Changeset 40523

Timestamp:

Aug 30, 2018, 5:07:23 PM (8 years ago)

Author:

eugene

Message:

add avselect function to dvo (returns all matched objects within a radius of a list of points)

Location:

trunk/Ohana/src/opihi

Files:

• dvo/Makefile (modified) (1 diff)
• dvo/avmatch.c (modified) (1 diff)
• dvo/avperiodomatch.c (modified) (1 diff)
• dvo/avselect.c (added)
• dvo/find_matches.c (modified) (2 diffs)
• dvo/init.c (modified) (2 diffs)
• dvo/mmatch.c (modified) (5 diffs)
• include/dvoshell.h (modified) (2 diffs)

trunk/Ohana/src/opihi/dvo/Makefile

@@ -46 +46 @@
 $(SRC)/avextract.$(ARCH).o              \
 $(SRC)/avmatch.$(ARCH).o                \
+$(SRC)/avselect.$(ARCH).o               \
 $(SRC)/avperiodogram.$(ARCH).o          \
 $(SRC)/avperiodomatch.$(ARCH).o         \

trunk/Ohana/src/opihi/dvo/avmatch.c

@@ -209 +209 @@
     }
 
-    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
+    find_matches_by_vectors_closest (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
 
     for (j = 0; (j < NPTS) && !interrupt; j++) {

trunk/Ohana/src/opihi/dvo/avperiodomatch.c

@@ -264 +264 @@
     dvo_catalog_unlock (&catalog);
 
-    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, idxValue);
+    find_matches_by_vectors_closest (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, idxValue);
 
     // Scan the catalog for objects which match the WHERE clause

trunk/Ohana/src/opihi/dvo/find_matches.c

@@ -5 +5 @@
 // attempt to match every RA,DEC entry with an entry in the catalog.  the result is stored in 'index'
 // index >= 0 : valid match; index == -2 : no valid match; index == -1 : not contained by this catalog
-int find_matches_by_vectors (SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index) {
+int find_matches_by_vectors_closest (SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index) {
 
   off_t i, j, J, Jmin, Npoints, Nmatch;
@@ -148 +148 @@
   return (TRUE);
 }
+
+// attempt to match every RA,DEC entry with an entry in the catalog.  the result is stored in 'result'
+AvselectResult *find_matches_by_vectors_allmatch (SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *nresult) {
+
+  off_t i, j, J;
+  double dX, dY, dR, Doff, Dmin, Dmax, Rmin, Rmax;
+  int status;
+  double RADIUS2;
+  Coords tcoords;
+
+  // XXX dvoconvert does not correctly maintain 'sorted' 
+  // assert(catalog[0].sorted);
+
+  off_t Npoints = RAvec->Nelements;
+
+  /** allocate local arrays (points) **/
+  ALLOCATE_PTR (X1, double, Npoints);
+  ALLOCATE_PTR (Y1, double, Npoints);
+  ALLOCATE_PTR (N1, off_t,  Npoints);
+
+  ALLOCATE_PTR (inCatalog, int,  Npoints);
+
+  /** allocate local arrays (catalog) **/
+  off_t Nave = catalog[0].Naverage;
+  ALLOCATE_PTR (X2, double, Nave);
+  ALLOCATE_PTR (Y2, double, Nave);
+  ALLOCATE_PTR (N2, off_t,  Nave);
+
+  off_t Nresult = 0;
+  off_t NRESULT = 1000;
+  ALLOCATE_PTR (result, AvselectResult, NRESULT);
+
+  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has
+   * the advantage that every point in R,D has a mapping to a unique X,Y.  However, note
+   * that not all possible X,Y points map back to R,D and the local plate scale changes
+   * far from the projection pole. We use the center of the region (catalog) for crval1,2.
+   */
+  InitCoords (&tcoords, "DEC--ARC");
+  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
+  if (region[0].Dmax < 90) {
+    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
+  } else {
+    tcoords.crval2 = 90.0;
+  }
+  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+
+  // this region includes a boundary layer of size RADIUS
+  if (fabs(region[0].Dmin) < fabs(region[0].Dmax)) {
+    Doff = RAD_DEG*region[0].Dmax;
+  } else {
+    Doff = RAD_DEG*region[0].Dmin;
+  }    
+  if (Doff < 80) {
+    Rmin = region[0].Rmin - RADIUS / 3600.0 / cos(Doff);
+    Rmax = region[0].Rmax + RADIUS / 3600.0 / cos(Doff);
+  } else {
+    Rmin = 0.0;
+    Rmax = 360.0;
+  }
+  Dmin = region[0].Dmin - RADIUS / 3600.0;
+  Dmax = region[0].Dmax + RADIUS / 3600.0;
+
+  // identify the entries contained by this catalog
+  for (i = 0; i < Npoints; i++) {
+    inCatalog[i] = FALSE;
+    if (RAvec->elements.Flt[i] < Rmin) continue;
+    if (RAvec->elements.Flt[i] > Rmax) continue;
+    if (DECvec->elements.Flt[i] < Dmin) continue;
+    if (DECvec->elements.Flt[i] > Dmax) continue;
+    inCatalog[i] = TRUE;
+  }
+
+  /* build spatial index (RA sort) referencing input array sequence */
+  for (i = 0; i < Npoints; i++) {
+    // we need to have a finite number for the sort below; inCatalog == 0 entries are skipped
+    // in the matching process
+    X1[i] = Y1[i] = 0.0; 
+    N1[i] = i;
+    if (!inCatalog[i]) continue;
+    status = RD_to_XY (&X1[i], &Y1[i], RAvec->elements.Flt[i], DECvec->elements.Flt[i], &tcoords);
+    assert (status);
+  }
+  if (Npoints > 1) sort_coords_index (X1, Y1, N1, Npoints);
+
+  /* build spatial index (RA sort) */
+  for (i = 0; i < Nave; i++) {
+    RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
+    N2[i] = i;
+  }
+  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
+
+  /* choose a radius for matches */
+  RADIUS2 = RADIUS*RADIUS;
+
+# define NEXTi { i++; continue; }
+# define NEXTj { j++; continue; }
+                                                                                                                                                                                 
+  /** find matched stars **/
+  for (i = j = 0; (i < Npoints) && (j < Nave); ) {
+    if (!inCatalog[N1[i]]) NEXTi;
+    if (!finite(X1[i]) || !finite(Y1[i])) NEXTi;
+    if (!finite(X2[j]) || !finite(Y2[j])) NEXTj;
+    
+    /* negative dX: j is too large */
+    dX = X1[i] - X2[j];
+    if (dX <= -1.02*RADIUS) NEXTi;
+
+    /* positive dX, i is too large */
+    if (dX >= 1.02*RADIUS) NEXTj;
+
+    /* within match range; look for matches */
+    for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) {
+      /* find closest match for this detection */
+      dX = X1[i] - X2[J];
+      dY = Y1[i] - Y2[J];
+      dR = dX*dX + dY*dY;
+      if (dR > RADIUS2) continue;
+
+      /*** a match is found, set the result for this entry ***/
+      result[Nresult].Ncat = N2[J];
+      result[Nresult].Nseq = N1[i];
+      result[Nresult].Roff = sqrt(dR);
+
+      Nresult ++;
+      CHECK_REALLOCATE (result, AvselectResult, NRESULT, Nresult, 1000);
+    }
+    NEXTi;
+  }
+
+  free (X1);
+  free (Y1);
+  free (N1);
+  free (N2);
+  free (X2);
+  free (Y2);
+
+  *nresult = Nresult;
+  return (result);
+}

trunk/Ohana/src/opihi/dvo/init.c

@@ -3 +3 @@
 int avextract       PROTO((int, char **));
 int avmatch         PROTO((int, char **));
+int avselect        PROTO((int, char **));
 int avperiodogram   PROTO((int, char **));
 int avperiodomatch  PROTO((int, char **));
@@ -69 +70 @@
   {1, "avextract",   avextract,    "extract average data values"},
   {1, "avmatch",     avmatch,      "extract average data values matched to RA,DEC points"},
+  {1, "avselect",    avselect,     "extract average data values within range of a list of RA,DEC points"},
   {1, "avperiodogram",  avperiodogram, "perform periodogram on objects based on restrictions"},
   {1, "avperiodomatch", avperiodomatch, "perform periodogram on objects based on ra,dec list"},

trunk/Ohana/src/opihi/dvo/mmatch.c

@@ -84 +84 @@
 
   // parse the fields to be extracted and returned
-  int first = 3;
-  if (CoordsFile) {
-    first = 1;
-  }
+  // this is a syntax check that we can perform before reading the input vectors
+  // NOTE: This block is missing in avmatch / avselect
+  int first = CoordsFile ? 1 : 3;
   fields = dbCmdlineFields (argc-first, &argv[first], DVO_TABLE_MEASURE, &last, &Nfields);
   if (fields == NULL) goto help;
@@ -166 +165 @@
     }      
 
+    // NOTE: image metadata is only needed to mmatch (not avmatch) since measures are tied to images, not averages
     if (loadImages) {
       Image *image;
@@ -245 +245 @@
     char hostfile[1024];
     snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name);
+
     dvo_catalog_init (&catalog, TRUE);
     catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
@@ -265 +266 @@
 
     // returns the matches to AVERAGE objects; what do we do for each MEASURE?
-    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
+    find_matches_by_vectors_closest (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
 
     for (j = 0; (j < Nelem) && !interrupt; j++) {
@@ -331 +332 @@
   ClearInterrupt (old_sigaction);
 
+  // XXXX ???? this seems odd
   for (n = 0; n < Nfields; n++) {
     ResetVector (vec[n], fields[n].type, Npts);

trunk/Ohana/src/opihi/include/dvoshell.h

@@ -52 +52 @@
 } PeriodogramResult;
 
+typedef struct {
+  off_t Ncat;
+  off_t Nseq;
+  float Roff;
+} AvselectResult;
+
 /*** Periodogram functions (avperiodogram and avperiodomatch) ***/
 void PeriodogramResultFree (PeriodogramResult *result);
@@ -98 +104 @@
 SkyList      *SelectRegionsByCoordVectors PROTO((Vector *RA, Vector *DEC));
 
-int           find_matches_by_vectors PROTO((SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index));
+int             find_matches_by_vectors_closest PROTO((SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index));
+AvselectResult *find_matches_by_vectors_allmatch PROTO((SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *Nresult));
 
 int           SetImageSelection     PROTO((int mode, SkyRegionSelection *selection));