Changeset 14356

Timestamp:

03/17/2010 04:52:35 PM (6 months ago)

Author:

dstn

Message:

dynamic indexing script

Location:

trunk/src/astrometry

Files:

• blind/dynamic-index.py (added)
• util/Makefile (modified) (1 diff)
• util/healpix-utils.c (modified) (4 diffs)
• util/healpix-utils.h (modified) (1 diff)
• util/healpix.c (modified) (4 diffs)
• util/healpix.h (modified) (2 diffs)
• util/healpix.py (modified) (1 diff)
• util/svn.c (modified) (1 diff)

trunk/src/astrometry/util/Makefile

@@ -176 +176 @@
         $(CC) $(SHAREDLIBFLAGS) -o $@ $^
 
-_healpix.so: healpix.o permutedsort.o starutil.o mathutil.o errors.o bl.o log.o gnu-specific.o
+_healpix.so: healpix-utils.o healpix.o permutedsort.o starutil.o mathutil.o errors.o bl.o log.o gnu-specific.o
         $(CC) $(SHAREDLIBFLAGS) -o $@ $^
 

trunk/src/astrometry/util/healpix-utils.c

@@ -79 +79 @@
 
 
-il* healpix_approx_rangesearch(double* xyz, double radius, int Nside, il* hps) {
+static il* hp_rangesearch(const double* xyz, double radius, int Nside, il* hps, bool approx) {
         int hp;
         double hprad = arcmin2dist(healpix_side_length_arcmin(Nside)) * sqrt(2);
@@ -96 +96 @@
                 nn = healpix_get_neighbours(hp, neighbours, Nside);
                 for (i=0; i<nn; i++) {
+                        bool tst;
                         double nxyz[3];
                         if (il_contains(frontier, neighbours[i]))
@@ -103 +104 @@
                         if (il_contains(hps, neighbours[i]))
                                 continue;
-                        healpix_to_xyzarr(neighbours[i], Nside, 0.5, 0.5, nxyz);
-                        if (sqrt(distsq(xyz, nxyz, 3)) - hprad < radius) {
+                        if (approx) {
+                                healpix_to_xyzarr(neighbours[i], Nside, 0.5, 0.5, nxyz);
+                                tst = (sqrt(distsq(xyz, nxyz, 3)) - hprad <= radius);
+                        } else {
+                                tst = (healpix_distance_to_xyz(neighbours[i], Nside, xyz, NULL) <= radius);
+                        }
+                        if (tst) {
                                 // in range!
                                 il_append(frontier, neighbours[i]);
@@ -118 +124 @@
         return hps;
 }
+
+il* healpix_rangesearch_xyz_approx(const double* xyz, double radius, int Nside, il* hps) {
+        return hp_rangesearch(xyz, radius, Nside, hps, TRUE);
+}
+
+il* healpix_rangesearch_xyz(const double* xyz, double radius, int Nside, il* hps) {
+        return hp_rangesearch(xyz, radius, Nside, hps, FALSE);
+}
+
+il* healpix_rangesearch_radec_approx(double ra, double dec, double radius, int Nside, il* hps) {
+        double xyz[3];
+        radecdeg2xyzarr(ra, dec, xyz);
+        return hp_rangesearch(xyz, radius, Nside, hps, TRUE);
+}
+
+il* healpix_rangesearch_radec(double ra, double dec, double radius, int Nside, il* hps) {
+        double xyz[3];
+        radecdeg2xyzarr(ra, dec, xyz);
+        return hp_rangesearch(xyz, radius, Nside, hps, FALSE);
+}

trunk/src/astrometry/util/healpix-utils.h

@@ -23 +23 @@
 
 /**
- Returns healpixes that may be within range of the given point.
+ Returns healpixes that are / may be within range of the given point, resp.
  */
-il* healpix_approx_rangesearch(double* xyz, double radius, int Nside, il* hps);
+il* healpix_rangesearch_xyz(const double* xyz, double radius, int Nside, il* hps);
+il* healpix_rangesearch_xyz_approx(const double* xyz, double radius, int Nside, il* hps);
+il* healpix_rangesearch_radec_approx(double ra, double dec, double radius, int Nside, il* hps);
+il* healpix_rangesearch_radec(double ra, double dec, double radius, int Nside, il* hps);
 
 /**

trunk/src/astrometry/util/healpix.c

@@ -994 +994 @@
 }
 
-int xyzarrtohealpix(double* xyz, int Nside) {
+int xyzarrtohealpix(const double* xyz, int Nside) {
         return xyztohealpix(xyz[0], xyz[1], xyz[2], Nside);
 }
 
-int64_t xyzarrtohealpixl(double* xyz, int Nside) {
+int64_t xyzarrtohealpixl(const double* xyz, int Nside) {
         return xyztohealpixl(xyz[0], xyz[1], xyz[2], Nside);
 }
 
-int xyzarrtohealpixf(double* xyz,int Nside, double* p_dx, double* p_dy) {
+int xyzarrtohealpixf(const double* xyz,int Nside, double* p_dx, double* p_dy) {
     return xyztohealpixf(xyz[0], xyz[1], xyz[2], Nside, p_dx, p_dy);
 }
@@ -1320 +1320 @@
 }
 
-double healpix_distance_to_radec(int hp, int Nside, double ra, double dec,
-                                                                 double* closestxyz) {
+double healpix_distance_to_xyz(int hp, int Nside, const double* xyz,
+                                                           double* closestxyz) {
         int thehp;
-        double xyz[3];
         // corners
         double cdx[4], cdy[4];
@@ -1338 +1337 @@
 
         // If the point is actually inside the healpix, dist = 0.
-        thehp = radecdegtohealpix(ra, dec, Nside);
+        thehp = xyzarrtohealpix(xyz, Nside);
         if (thehp == hp) {
                 if (closestxyz)
-                        radecdeg2xyzarr(ra, dec, closestxyz);
+                        memcpy(closestxyz, xyz, 3*sizeof(double));
                 return 0;
         }
 
         // Find two nearest corners.
-        radecdeg2xyzarr(ra, dec, xyz);
         for (i=0; i<4; i++) {
                 double cxyz[3];
@@ -1407 +1405 @@
 }
 
+double healpix_distance_to_radec(int hp, int Nside, double ra, double dec,
+                                                                 double* closestxyz) {
+        double xyz[3];
+        radecdeg2xyzarr(ra, dec, xyz);
+        return healpix_distance_to_xyz(hp, Nside, xyz, closestxyz);
+}
+
 void healpix_radec_bounds(int hp, int nside,
                                                   double* p_ralo, double* p_rahi,

trunk/src/astrometry/util/healpix.h

@@ -267 +267 @@
    a healpix index.
 */
-int xyzarrtohealpix(double* xyz, int Nside);
-
-int64_t xyzarrtohealpixl(double* xyz, int Nside);
-
-int xyzarrtohealpixf(double* xyz,int Nside, double* p_dx, double* p_dy);
+int xyzarrtohealpix(const double* xyz, int Nside);
+
+int64_t xyzarrtohealpixl(const double* xyz, int Nside);
+
+int xyzarrtohealpixf(const double* xyz,int Nside, double* p_dx, double* p_dy);
 
 /**
@@ -362 +362 @@
 
 /**
+ Returns the minimum distance (in degrees) between the given healpix
+ and the given xyz (point on unit sphere).
+ */
+double healpix_distance_to_xyz(int hp, int Nside, const double* xyz,
+                                                           double* closestxyz);
+
+/**
  Computes the RA,Dec bounding-box of the given healpix.  Results are
  in degrees.  RA may be wacky for healpixes spanning RA=0.

trunk/src/astrometry/util/healpix.py

@@ -20 +20 @@
         raise IOError('_healpix.so library not found')
 
-
+def il_to_list(il):
+        _lib.il_size.argtypes = [c_void_p]
+        _lib.il_size.restype = c_int
+        _lib.il_get.argtypes = [c_void_p, c_int]
+        _lib.il_get.restype = c_int
+        N = _lib.il_size(il)
+        lst = []
+        for i in range(N):
+                x = _lib.il_get(il, c_int(i))
+                lst.append(int(x))
+        return lst
+
+def healpix_rangesearch(ra, dec, radius, nside):
+        _lib.healpix_rangesearch_radec.restype = c_void_p
+        hplist = _lib.healpix_rangesearch_radec(c_double(ra), c_double(dec), c_double(radius), c_int(nside), c_void_p(None))
+        lst = il_to_list(hplist)
+        _lib.il_free.argtypes = [c_void_p]
+        _lib.il_free(hplist)
+        return lst
 
 def healpix_nside_for_side_length_arcmin(arcmin):

trunk/src/astrometry/util/svn.c

@@ -139 +139 @@
 //
 //
+//
+//
+//
+//