Changeset 40575

Timestamp:

Dec 5, 2018, 2:13:03 PM (8 years ago)

Author:

eugene

Message:

fix avselect function to cross dvo catalog boundaries

Location:

trunk/Ohana/src/opihi

Files:

• dvo/avselect.c (modified) (1 diff)
• dvo/skyregion.c (modified) (1 diff)
• include/dvoshell.h (modified) (1 diff)

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

@@ -102 +102 @@
   }
 
+  RADIUS = atof (argv[1]);
+
   /* load regions which contain all supplied RA,DEC coordinates */
-  if ((skylist = SelectRegionsByCoordVectors (RAvec, DECvec)) == NULL) goto escape;
+  if ((skylist = SelectRegionsByCoordVectorsAndRadius (RAvec, DECvec, RADIUS/3600.0)) == NULL) goto escape;
 
   // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results

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

@@ -147 +147 @@
 }
 
+void sort_skylist_by_index (SkyList *list) {
+
+# define SWAPFUNC(A,B){ \
+    SkyRegion *tmp; tmp = list[0].regions[A]; list[0].regions[A] = list[0].regions[B]; list[0].regions[B] = tmp; \
+    char *tmpname; tmpname = list[0].filename[A]; list[0].filename[A] = list[0].filename[B]; list[0].filename[B] = tmpname; \
+}
+# define COMPARE(A,B)(list[0].regions[A][0].index < list[0].regions[B][0].index)
+
+  OHANA_SORT (list[0].Nregions, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+SkyList *SkyListUniqueSubset (SkyList *input) {
+
+  sort_skylist_by_index (input);
+  for (int i = 0; i < input->Nregions; i++) {
+    fprintf (stderr, "idx: %s : %s : %d\n", input->regions[i][0].name, input->filename[i], input->regions[i][0].index);
+  }
+
+  SkyList *list;
+  
+  int Nout = 0;
+  int NOUT = input->Nregions;
+
+  ALLOCATE (list, SkyList, 1);
+  ALLOCATE (list[0].regions, SkyRegion *, NOUT);
+  ALLOCATE (list[0].filename, char *, NOUT);
+  list[0].ownElements = FALSE; // this list is only holding a view to the elements
+  
+  if (input->Nregions == 0) {
+    list->Nregions = 0;
+    return list;
+  }
+
+  list->regions[Nout] = input->regions[0];
+  list->filename[Nout] = input->filename[0];
+
+  for (int i = 0; i < input->Nregions; i++) {
+    if (list->regions[Nout][0].index == input->regions[i][0].index) continue;
+    Nout ++;
+    list->regions[Nout] = input->regions[i];
+    list->filename[Nout] = input->filename[i];
+  }
+  list->Nregions = Nout;
+  REALLOCATE (list[0].regions, SkyRegion *, Nout);
+  REALLOCATE (list[0].filename, char *, Nout);
+
+  // free the input list?
+  return list;
+}
+
+// return a list of all regions covered by this list of coordinates and the given radius
+SkyList *SelectRegionsByCoordVectorsAndRadius (Vector *RA, Vector *DEC, float Radius) {
+  
+  int Npts, Nout, NOUT, *found;
+  double ra, dec;
+  SkyList *list;
+  SkyTable *sky;
+
+  sky = GetSkyTable ();
+
+  Npts = RA->Nelements;
+  ALLOCATE (found, int, Npts);
+  memset (found, 0, Npts*sizeof(int));
+
+  // output list
+  Nout = 0;
+  NOUT = 100;
+  ALLOCATE (list, SkyList, 1);
+  ALLOCATE (list[0].regions,  SkyRegion *, NOUT);
+  ALLOCATE (list[0].filename,  char *, NOUT);
+  list[0].Nregions = 0;
+  list[0].ownElements = FALSE; // this list is only holding a view to the elements
+  strcpy (list[0].hosts, sky[0].hosts);
+
+  for (int i = 0; i < Npts; i++) {
+    if (found[i]) continue;
+
+    SkyList *new = SkyListByRadius (sky, -1, RA->elements.Flt[i], DEC->elements.Flt[i], Radius);
+    if (new->Nregions == 0) continue;
+
+    // append new regions to output list
+    if (Nout + new->Nregions >= NOUT) {
+        NOUT += new->Nregions + 100;
+        REALLOCATE (list[0].regions, SkyRegion *, NOUT);
+        REALLOCATE (list[0].filename, char *, NOUT);
+    }
+    
+    for (int n = 0; n < new->Nregions; n++) {
+      list[0].regions[Nout] = new[0].regions[n];
+      list[0].filename[Nout] = new[0].filename[n];
+      Nout ++;
+    }
+    found[i] = TRUE;
+
+    // scan over the remaining points to find any that lie within these regions
+    // if we sorted the ra,dec inputs we could break out of this more quickly
+    for (int j = i + 1; j < Npts; j++) {
+      ra = RA->elements.Flt[j];
+      dec = DEC->elements.Flt[j];
+      for (int n = 0; !found[j] && (n < new->Nregions); n++) {
+        if (ra  < new[0].regions[n][0].Rmin) continue;
+        if (ra  > new[0].regions[n][0].Rmax) continue;
+        if (dec < new[0].regions[n][0].Dmin) continue;
+        if (dec > new[0].regions[n][0].Dmax) continue;
+        found[j] = TRUE;
+      }
+    }
+    SkyListFree (new);
+  }
+  list[0].Nregions = Nout;
+
+  // before I go any further, I need to strip down to a unique subset
+  SkyList *uniq = SkyListUniqueSubset(list);
+
+  free (found);
+
+  return (uniq);
+}
+

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

@@ -103 +103 @@
 int           SkyRegionByPoint_r    PROTO((SkyTable *table, SkyList *list, int depth, double ra, double dec));
 SkyList      *SelectRegionsByCoordVectors PROTO((Vector *RA, Vector *DEC));
+SkyList      *SelectRegionsByCoordVectorsAndRadius (Vector *RA, Vector *DEC, float Radius);
 
 int             find_matches_by_vectors_closest PROTO((SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index));