Skip to content
This repository has been archived by the owner on Dec 31, 2018. It is now read-only.

Commit

Permalink
New API for a number of pixel-related functions
Browse files Browse the repository at this point in the history 
All the functions that convert a pixel index into a direction and
vice versa do no longer accept a NSIDE parameter, but instead
a hpix_resolution_t structure. This allows to avoid repeating some
calculations a lot of time.
  • Loading branch information
@ziotom78
ziotom78 committed Dec 5, 2013
1 parent 785821a commit 858c990
Show file tree
Hide file tree
Showing 5 changed files with 145 additions and 110 deletions.
9 changes: 6 additions & 3 deletions src/bitmap.c
View file
Expand Up @@ -201,7 +201,7 @@ hpix_bmp_projection_trace(const hpix_bmp_projection_t * proj,
}

hpix_pixel_num_t pixel_idx =
angles_to_pixel_fn(nside, theta, phi);
angles_to_pixel_fn(hpix_map_resolution(map), theta, phi);
if(pixels[pixel_idx] > -1.6e+30)
*line_ptr = pixels[pixel_idx];
else
Expand All @@ -214,14 +214,17 @@ hpix_bmp_projection_trace(const hpix_bmp_projection_t * proj,
if(min_value || max_value)
{
if(min_value)
*min_value = FLT_MAX;
*min_value = DBL_MAX;

if(max_value)
*max_value = FLT_MIN;
*max_value = -DBL_MAX;

double * bitmap_ptr = bitmap;
for(size_t idx = 0; idx < num_of_pixels; ++idx, ++bitmap_ptr)
{
if(isnan(*bitmap_ptr) || isinf(*bitmap_ptr))
continue;

if(min_value && *min_value > *bitmap_ptr)
*min_value = *bitmap_ptr;
if(max_value && *max_value < *bitmap_ptr)
Expand Down
45 changes: 25 additions & 20 deletions src/hpixlib/hpix.h
View file
Expand Up @@ -147,6 +147,9 @@ double * hpix_map_pixels(const hpix_map_t * map);

size_t hpix_map_num_of_pixels(const hpix_map_t * map);

void hpix_init_resolution_from_nside(hpix_nside_t nside,
hpix_resolution_t * resolution);

const hpix_resolution_t * hpix_map_resolution(const hpix_map_t * map);

/* Functions implemented in integer_functions.c */
Expand Down Expand Up @@ -225,48 +228,50 @@ hpix_save_fits_pol_map(const char * file_name,
void hpix_angles_to_vector(double theta, double phi,
hpix_vector_t * vector);

typedef hpix_pixel_num_t hpix_angles_to_pixel_fn_t(hpix_nside_t,
double, double);
typedef hpix_pixel_num_t hpix_angles_to_pixel_fn_t(const hpix_resolution_t *,
double, double);

hpix_pixel_num_t hpix_angles_to_ring_pixel(hpix_nside_t nside,
double theta,
double phi);
hpix_pixel_num_t hpix_angles_to_ring_pixel(const hpix_resolution_t * resolution,
double theta,
double phi);

hpix_pixel_num_t hpix_angles_to_nest_pixel(hpix_nside_t nside,
double theta,
double phi);
hpix_pixel_num_t hpix_angles_to_nest_pixel(const hpix_resolution_t * resolution,
double theta,
double phi);

void hpix_vector_to_angles(const hpix_vector_t * vector,
double * theta, double * phi);

typedef hpix_pixel_num_t hpix_vector_to_pixel_fn_t(hpix_nside_t,
const hpix_vector_t *);

hpix_pixel_num_t hpix_vector_to_ring_pixel(hpix_nside_t nside,
hpix_pixel_num_t hpix_vector_to_ring_pixel(const hpix_resolution_t * resolution,
const hpix_vector_t * vector);

hpix_pixel_num_t hpix_vector_to_nest_pixel(hpix_nside_t nside,
hpix_pixel_num_t hpix_vector_to_nest_pixel(const hpix_resolution_t * resolution,
const hpix_vector_t * vector);

typedef void hpix_pixel_to_angles(hpix_nside_t, hpix_pixel_num_t,
double *, double *);

void hpix_ring_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,
double * theta, double * phi);
void hpix_ring_pixel_to_angles(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel,
double * theta, double * phi);

void hpix_nest_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,
double * theta, double * phi);
void hpix_nest_pixel_to_angles(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel,
double * theta, double * phi);

typedef void hpix_pixel_to_vector(hpix_nside_t, hpix_pixel_num_t,
hpix_vector_t * vector);

void hpix_ring_pixel_to_vector(hpix_nside_t nside,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector);
void hpix_ring_pixel_to_vector(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector);

void hpix_nest_pixel_to_vector(hpix_nside_t nside,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector);
void hpix_nest_pixel_to_vector(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector);

/* Functions implemented in bitmap.c */

Expand Down
78 changes: 45 additions & 33 deletions src/positions.c
View file
Expand Up @@ -37,7 +37,7 @@

void
hpix_angles_to_vector(double theta, double phi,
hpix_vector_t * vector)
hpix_vector_t * vector)
{
assert(vector);

Expand All @@ -51,16 +51,19 @@ hpix_angles_to_vector(double theta, double phi,


hpix_pixel_num_t
hpix_angles_to_ring_pixel(hpix_nside_t nside,
hpix_angles_to_ring_pixel(const hpix_resolution_t * resolution,
double theta,
double phi)
{
assert(resolution != NULL);

int jp, jm, ipix1;
int ir, ip;

int nl2 = 2*nside;
int nl4 = 4*nside;
int ncap = nl2*(nside-1);
hpix_nside_t nside = resolution->nside;
hpix_nside_t nl2 = resolution->nside_times_two;
hpix_nside_t nl4 = resolution->nside_times_four;
unsigned int ncap = resolution->ncap;

double z = cos(theta);
double z_abs = fabs(z);
Expand Down Expand Up @@ -97,7 +100,7 @@ hpix_angles_to_ring_pixel(hpix_nside_t nside,
ipix1 = 2 * ir * (ir - 1) + ip;
if (z <= 0.)
{
ipix1 = 12 * nside * nside - 2 * ir * (ir + 1) + ip;
ipix1 = resolution->num_of_pixels - 2 * ir * (ir + 1) + ip;
}
}
return ipix1 - 1;
Expand All @@ -107,10 +110,13 @@ hpix_angles_to_ring_pixel(hpix_nside_t nside,


hpix_pixel_num_t
hpix_angles_to_nest_pixel(hpix_nside_t nside,
hpix_angles_to_nest_pixel(const hpix_resolution_t * resolution,
double theta,
double phi)
{
assert(resolution != NULL);

hpix_nside_t nside = resolution->nside;
double z, z_abs, tt, tp, tmp;
int face_num,jp,jm;
long ifp, ifm;
Expand Down Expand Up @@ -193,44 +199,47 @@ hpix_vector_to_angles(const hpix_vector_t * vector,


hpix_pixel_num_t
hpix_vector_to_ring_pixel(hpix_nside_t nside,
const hpix_vector_t * vector)
hpix_vector_to_ring_pixel(const hpix_resolution_t * resolution,
const hpix_vector_t * vector)
{
double theta, phi;
hpix_vector_to_angles(vector, &theta, &phi);
return hpix_angles_to_ring_pixel(nside, theta, phi);
return hpix_angles_to_ring_pixel(resolution, theta, phi);
}

/**********************************************************************/


hpix_pixel_num_t
hpix_vector_to_nest_pixel(hpix_nside_t nside,
const hpix_vector_t * vector)
hpix_vector_to_nest_pixel(const hpix_resolution_t * resolution,
const hpix_vector_t * vector)
{
double theta, phi;
hpix_vector_to_angles(vector, &theta, &phi);
return hpix_angles_to_nest_pixel(nside, theta, phi);
return hpix_angles_to_nest_pixel(resolution, theta, phi);
}

/**********************************************************************/


void
hpix_ring_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,
hpix_ring_pixel_to_angles(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel,
double * theta, double * phi)
{
assert(resolution);
assert(theta && phi);

int nl2, nl4, ncap, iring, iphi, ip, ipix1;
double fact1, fact2, fodd, hip, fihip;
hpix_nside_t nside = resolution->nside;
hpix_nside_t nl2 = resolution->nside_times_two;
hpix_nside_t nl4 = resolution->nside_times_four;
unsigned int ncap = resolution->ncap;
int iring, iphi, ip, ipix1;
double fact1, fact2, fodd, hip, fihip;

ipix1 = pixel + 1; // in {1, npix}
nl2 = 2*nside;
nl4 = 4*nside;
ncap = 2*nside*(nside-1);
fact1 = 1.5*nside;
fact2 = 3.0*nside*nside;
fact1 = 1.5 * nside;
fact2 = 3.0 * resolution->pixels_per_face;

if(ipix1 <= ncap)
{
Expand All @@ -257,7 +266,7 @@ hpix_ring_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,
}
else {//! South Polar cap -----------------------------------

ip = 12 * nside * nside - ipix1 + 1;
ip = resolution->num_of_pixels - ipix1 + 1;
hip = ip / 2.;
/* bug corrige floor instead of 1.* */
fihip = floor(hip);
Expand All @@ -273,9 +282,11 @@ hpix_ring_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,


void
hpix_nest_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,
hpix_nest_pixel_to_angles(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel,
double * theta, double * phi)
{
assert(resolution);
assert(theta && phi);

int ix, iy, jrt, jr, nr, jpt, jp, kshift;
Expand All @@ -285,12 +296,13 @@ hpix_nest_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,
int jrll[12] = { 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 };
int jpll[12] = { 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 };

hpix_nside_t nside = resolution->nside;
double fn = 1.*nside;
double fact1 = 1./(3.*fn*fn);
double fact2 = 2./(3.*fn);
int nl4 = 4*nside;
hpix_nside_t nl4 = resolution->nside_times_four;

int npface = nside*nside;
unsigned int npface = resolution->pixels_per_face;

int face_num = pixel/npface; // face number in {0,11}
int ipf = (int)fmod(pixel,npface); // pixel number in the face {0,npface-1}
Expand Down Expand Up @@ -337,28 +349,28 @@ hpix_nest_pixel_to_angles(hpix_nside_t nside, hpix_pixel_num_t pixel,


void
hpix_ring_pixel_to_vector(hpix_nside_t nside,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector)
hpix_ring_pixel_to_vector(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector)
{
assert(vector);

double theta, phi;
hpix_ring_pixel_to_angles(nside, pixel_index, &theta, &phi);
hpix_ring_pixel_to_angles(resolution, pixel_index, &theta, &phi);
hpix_angles_to_vector(theta, phi, vector);
}

/**********************************************************************/


void
hpix_nest_pixel_to_vector(hpix_nside_t nside,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector)
hpix_nest_pixel_to_vector(const hpix_resolution_t * resolution,
hpix_pixel_num_t pixel_index,
hpix_vector_t * vector)
{
assert(vector);

double theta, phi;
hpix_nest_pixel_to_angles(nside, pixel_index, &theta, &phi);
hpix_nest_pixel_to_angles(resolution, pixel_index, &theta, &phi);
hpix_angles_to_vector(theta, phi, vector);
}
1 change: 1 addition & 0 deletions test/test_bmp_projection.c
View file
Expand Up @@ -79,6 +79,7 @@ START_TEST(mollweide_projection)
{
fail_unless(mollweide_proj_min == 0.0,
"hpix_bmp_to_mollweide_proj failed to set the minimum value");
fprintf(stderr, "***************** maximum = %f\n", mollweide_proj_max);
fail_unless(mollweide_proj_max == 767.0,
"hpix_bmp_to_mollweide_proj failed to set the maximum value");

Expand Down

0 comments on commit 858c990

Please sign in to comment.