Commit af837613 authored by Leo Pound Singer's avatar Leo Pound Singer

Finishing touches to new multi-order HEALPix functions

Original: f0551ce4104f392a36075c9c6e3f59f4fa894ed2
parent 5cece41e
......@@ -155,9 +155,9 @@ static void nest2uniq_loop(
for (npy_intp i = 0; i < n; i ++)
{
*(int64_t *) &args[2][i * steps[2]] = nest2uniq64(
*(int8_t *) &args[0][i * steps[0]],
*(int64_t *) &args[1][i * steps[1]]);
*(uint64_t *) &args[2][i * steps[2]] = nest2uniq64(
*(int8_t *) &args[0][i * steps[0]],
*(uint64_t *) &args[1][i * steps[1]]);
}
}
......@@ -169,9 +169,9 @@ static void uniq2nest_loop(
for (npy_intp i = 0; i < n; i ++)
{
*(int64_t *) &args[2][i * steps[2]] = *(int64_t *) &args[0][i * steps[0]];
*(int8_t *) &args[1][i * steps[1]] = uniq2nest64(
(int64_t *) &args[2][i * steps[2]]);
*(int8_t *) &args[1][i * steps[1]] = uniq2nest64(
*(uint64_t *) &args[0][i * steps[0]],
(uint64_t *) &args[2][i * steps[2]]);
}
}
......@@ -183,8 +183,8 @@ static void uniq2order_loop(
for (npy_intp i = 0; i < n; i ++)
{
*(int8_t *) &args[1][i * steps[1]] = uniq2order64(
*(int64_t *) &args[0][i * steps[0]]);
*(int8_t *) &args[1][i * steps[1]] = uniq2order64(
*(uint64_t *) &args[0][i * steps[0]]);
}
}
......
......@@ -33,9 +33,9 @@ uint64_t nest2uniq64(uint8_t order, uint64_t nest)
}
int8_t uniq2order64(uint64_t n)
int8_t uniq2order64(uint64_t uniq)
{
if (n < 4)
if (uniq < 4)
return -1;
int8_t order;
......@@ -43,35 +43,36 @@ int8_t uniq2order64(uint64_t n)
int64_t o;
asm("bsrq %1, %0\n\t"
: "=r" (o)
: "rm" (n));
: "rm" (uniq));
order = o;
#else
for (order = -1; n; n >>= 1, order ++)
for (order = -1; uniq; uniq >>= 1, order ++)
/* noop */;
#endif
return (order >> 1) - 1;
}
double uniq2pixarea64(uint64_t n)
double uniq2pixarea64(uint64_t uniq)
{
return ldexp(M_PI / 3, -2 * uniq2order64(n));
return ldexp(M_PI / 3, -2 * uniq2order64(uniq));
}
int8_t uniq2nest64(uint64_t *n)
int8_t uniq2nest64(uint64_t uniq, uint64_t *nest)
{
int8_t order = uniq2order64(*n);
*n -= (uint64_t) 1 << 2 * (order + 1);
int8_t order = uniq2order64(uniq);
*nest = uniq - ((uint64_t) 1 << 2 * (order + 1));
return order;
}
void pix2ang_uniq64(uint64_t ipix, double *theta, double *phi)
void pix2ang_uniq64(uint64_t uniq, double *theta, double *phi)
{
int8_t order = uniq2nest64(&ipix);
uint64_t nest;
int8_t order = uniq2nest64(uniq, &nest);
int64_t nside = (int64_t) 1 << order;
pix2ang_nest64(nside, ipix, theta, phi);
pix2ang_nest64(nside, nest, theta, phi);
}
......@@ -79,8 +80,7 @@ void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t
{
const size_t pixelsize = offset + itemsize;
const void *pixel = (const char *) pixels + (len - 1) * pixelsize;
uint64_t ipix = *(const uint64_t *) pixel;
const int8_t max_order = uniq2order64(ipix);
const int8_t max_order = uniq2order64(*(const uint64_t *) pixel);
*npix = 12 * ((size_t) 1 << 2 * max_order);
void *ret = calloc(*npix, itemsize);
if (!ret)
......@@ -89,11 +89,11 @@ void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t
for (size_t i = 0; i < len; i ++)
{
pixel = (const char *) pixels + i * pixelsize;
ipix = *(const uint64_t *) pixel;
const int8_t order = uniq2nest64(&ipix);
uint64_t nest;
const int8_t order = uniq2nest64(*(const uint64_t *) pixel, &nest);
const size_t reps = (size_t) 1 << 2 * (max_order - order);
for (size_t j = 0; j < reps; j ++)
memcpy((char *) ret + (ipix * reps + j) * itemsize, (const char *) pixel + offset, itemsize);
memcpy((char *) ret + (nest * reps + j) * itemsize, (const char *) pixel + offset, itemsize);
}
return ret;
......
......@@ -33,14 +33,14 @@
/* Convert a NESTED pixel index to NUNIQ ordering. */
uint64_t nest2uniq64(uint8_t order, uint64_t nest);
int8_t uniq2order64(uint64_t n);
int8_t uniq2order64(uint64_t uniq);
double uniq2pixarea64(uint64_t n);
double uniq2pixarea64(uint64_t uniq);
/* Convert a NUNIQ pixel index to NESTED ordering. */
int8_t uniq2nest64(uint64_t *n);
int8_t uniq2nest64(uint64_t uniq, uint64_t *nest);
void pix2ang_uniq64(uint64_t ipix, double *theta, double *phi);
void pix2ang_uniq64(uint64_t uniq, double *theta, double *phi);
void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t len, size_t *npix);
......
......@@ -1455,8 +1455,8 @@ static void test_nest2uniq64(uint8_t order, uint64_t nest, uint64_t uniq)
"expected nest2uniq64(%u, %llu) = %llu, got %llu",
(unsigned) order, nest, uniq, uniq_result);
uint64_t nest_result = uniq;
const uint8_t order_result = uniq2nest64(&nest_result);
uint64_t nest_result;
const uint8_t order_result = uniq2nest64(uniq, &nest_result);
gsl_test(!(nest_result == nest && order_result == order),
"expected uniq2nest64(%llu) = (%u, %llu), got (%u, %llu)",
uniq, (unsigned) order, nest, order_result, nest_result);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment