Changeset 17284

Timestamp:

Apr 2, 2008, 12:37:15 PM (18 years ago)

Author:

eugene

Message:

add function to join multiple detection sets into single objects

Location:

trunk/Ohana/src/photdbc

Files:

• Makefile (modified) (2 diffs)
• include/photdbc.h (modified) (1 diff)
• src/ConfigInit.c (modified) (3 diffs)
• src/copy_images.c (modified) (4 diffs)
• src/join_stars.c (modified) (14 diffs)
• src/make_subcatalog.c (modified) (3 diffs)
• src/photdbc.c (modified) (1 diff)

trunk/Ohana/src/photdbc/Makefile

@@ -27 +27 @@
 $(SRC)/copy_images.$(ARCH).o       \
 $(SRC)/Shutdown.$(ARCH).o          \
+$(SRC)/join_stars.$(ARCH).o        \
 $(SRC)/make_subcatalog.$(ARCH).o        
 
@@ -37 +38 @@
 $(SRC)/load_catalogs.$(ARCH).o     \
 $(SRC)/gcatalog.$(ARCH).o          \
-$(SRC)/join_stars.$(ARCH).o        \
 $(SRC)/unique_measures.$(ARCH).o   \
 $(SRC)/flag_measures.$(ARCH).o     \

trunk/Ohana/src/photdbc/include/photdbc.h

@@ -112 +112 @@
 void SetProtect (int mode);
 int SetSignals ();
-
-
 int copy_images (char *outdir);

trunk/Ohana/src/photdbc/src/ConfigInit.c

@@ -26 +26 @@
   if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
 
-  // WarnConfig (config, "JOIN_RADIUS",            "%lf", 0, &JOIN_RADIUS);
+  WarnConfig (config, "PHOTDBC_JOIN_RADIUS",       "%lf", 0, &JOIN_RADIUS);
   // WarnConfig (config, "UNIQ_RADIUS",            "%lf", 0, &UNIQ_RADIUS);
 
@@ -59 +59 @@
   }
 
-  SetZeroPoint (ZERO_POINT);
-
   if (*CATMODE == 0) strcpy (CATMODE, "RAW");
   if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR");
@@ -70 +68 @@
     exit (1);
   }
+  SetZeroPoint (ZERO_POINT);
 
   free (config);

trunk/Ohana/src/photdbc/src/copy_images.c

@@ -10 +10 @@
   struct stat filestat;
   char *path;
+  char filename[1024];
 
   path = pathname (ImageCat);
@@ -28 +29 @@
   // load the input image data
   if (!dvo_image_load (&in, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", in.filename);
-  dvo_image_unlock (&in);
 
   // define output filename (replace CATDIR)
@@ -35 +35 @@
 
   // lock the output catalog
-  status = dvo_image_lock (&out, ImageOut, 60.0, LCK_SOFT);
+  status = dvo_image_lock (&out, ImageOut, 60.0, LCK_XCLD);
   if (!status) Shutdown ("ERROR: failure to lock image catalog %s", ImageOut);
   if (out.dbstate != LCK_EMPTY) Shutdown ("ERROR: image table exists %s", ImageOut);
@@ -46 +46 @@
   dvo_image_addrows (&out, image, Nimage);
 
-  dvo_image_save (&out, VERBOSE);
+  dvo_image_update (&out, VERBOSE);
   dvo_image_unlock (&out);
+
+  dvo_image_unlock (&in);
+
+  // create the photcode file
+  sprintf (filename, "%s/Photcodes.dat", outdir);
+  SavePhotcodesFITS (filename);
+
   return (TRUE);
 }

trunk/Ohana/src/photdbc/src/join_stars.c

@@ -8 +8 @@
   double *X, *Y, dX, dY, dR, RADIUS2;
   double Sr, Sd, Rmid, Dmid;
+  int basecode, baseNsec, Nsecfilt, Nfirst;
 
   Average *naverage, *average;
@@ -16 +17 @@
 
   if (VERBOSE) fprintf (stderr, "joining overlapping stars\n");
+  if (VERBOSE) fprintf (stderr, "require base star to have i-band\n");
+
+  basecode = GetPhotcodeCodebyName ("i");
+  baseNsec = GetPhotcodeNsec (basecode);
+  Nsecfilt = catalog[0].Nsecfilt;
 
   average = catalog[0].average;
@@ -26 +32 @@
 
   /* reference for coords is this region */
+  Naves = 0;
   Rmid = Dmid = 0;
   for (i = 0; i < Naverage; i++) {
+    // XXX we should be vigilant against R,D becoming nan : this must be due to pm...
+    if (isnan(average[i].R) || isnan(average[i].D)) {
+      average[i].R = 0.0;
+      average[i].D = 0.0;
+      continue;
+    }
     Rmid += average[i].R;
     Dmid += average[i].D;
+    Naves ++;
   }
   Rmid /= Naverage;
@@ -59 +73 @@
 
   ALLOCATE (naverage, Average, Naverage);
+  ALLOCATE (nsecfilt, SecFilt, Naverage*Nsecfilt);
   ALLOCATE (mpointer, Mpointer, Nmeasure);
   for (i = 0; i < Nmeasure; i++) mpointer[i].averef = -1;
@@ -69 +84 @@
     Nj = index[j];
 
+    // if ((average[Ni].R > 131.259) && (average[Ni].R < 131.267) && (average[Ni].D > 20.440) && (average[Ni].D < 20.450)) {
+    // fprintf (stderr, "outer: %f, %f - %f, %f (%f, %f) == (%f, %f)\n", average[Ni].R, average[Ni].D, average[Nj].R, average[Nj].D, 
+    // 3600.0*(average[Ni].R - average[Nj].R), 3600.0*(average[Ni].D - average[Nj].D), X[i] - X[j], Y[i] - Y[j]);
+    // }
+
+    // require base star to meet certain conditions:
+    if (isnan(secfilt[Ni*Nsecfilt + baseNsec].M)) {
+      i++;
+      continue;
+    }
+
     /* a new star, add it to naverage[] */
     if (!found[i]) {    
-      naverage[Naves]        = average[Ni];
-      naverage[Naves].offset = Nmeas;
-      for (k = 0; k < average[Ni].Nm; k++) {
-        m = average[Ni].offset + k;
+      naverage[Naves]               = average[Ni];
+      naverage[Naves].measureOffset = Nmeas;
+
+      for (k = 0; k < Nsecfilt; k++) {
+        nsecfilt[Naves*Nsecfilt + k] = secfilt[Ni*Nsecfilt + k];
+      }
+
+      for (k = 0; k < average[Ni].Nmeasure; k++) {
+        m = average[Ni].measureOffset + k;
         mpointer[Nmeas].measure = m;
         mpointer[Nmeas].averef  = Naves;
@@ -87 +118 @@
 
     if (found[j]) { j++; continue; }  // don't duplicate
-    if (j == i) {   j++; continue; }  // don't auto-correlate
+    if (j == i)   { j++; continue; }  // don't auto-correlate
 
     dX = X[i] - X[j];
@@ -102 +133 @@
     for (; (dX > -2*JOIN_RADIUS) && (j < Naverage); j++) {
       Nj = index[j];
+
+      // if ((average[Ni].R > 131.259) && (average[Ni].R < 131.267) && (average[Ni].D > 20.440) && (average[Ni].D < 20.450)) {
+      // fprintf (stderr, "inner: %f, %f - %f, %f (%f, %f) == (%f, %f)\n", average[Ni].R, average[Ni].D, average[Nj].R, average[Nj].D, 
+      // 3600.0*(average[Ni].R - average[Nj].R), 3600.0*(average[Ni].D - average[Nj].D), X[i] - X[j], Y[i] - Y[j]);
+      // }
+
       if (found[j]) continue;
       dX = X[i] - X[j];
@@ -108 +145 @@
       if (dR < RADIUS2) {  /* matched star, join to first */
 
+        // if ((average[Ni].R > 131.259) && (average[Ni].R < 131.267) && (average[Ni].D > 20.440) && (average[Ni].D < 20.450)) {
+        // fprintf (stderr, "match: %f, %f - %f, %f\n", average[Ni].R, average[Ni].D, average[Nj].R, average[Nj].D);
+        // }
+
         /* define pointers for new measures */
-        for (k = 0; k < average[Nj].Nm; k++) {
-          m = average[Nj].offset + k;
+        for (k = 0; k < average[Nj].Nmeasure; k++) {
+          m = average[Nj].measureOffset + k;
           mpointer[Nmeas].measure = m;
           mpointer[Nmeas].averef  = Ncurr;
@@ -117 +158 @@
           Nmeas ++;
         }
-        naverage[Ncurr].Nm += average[Nj].Nm;
+        Nfirst = average[Nj].Nmeasure;
+        naverage[Ncurr].Nmeasure += average[Nj].Nmeasure;
         found[j] = TRUE;
-        
+
+# if 0
         /* recalculate naverage[Ncurr].RA,DEC */
         /* this must be done here to keep the average position consistent
            for the next star found */
-        for (Sr = Sd = k = 0; k < naverage[Ncurr].Nm; k++) {
-          m = naverage[Ncurr].offset + k;
+        for (Sr = Sd = k = 0; k < naverage[Ncurr].Nmeasure; k++) {
+          m = naverage[Ncurr].measureOffset + k;
           Sr += mpointer[m].R;
           Sd += mpointer[m].D;
         }
-        Sr = Sr / naverage[Ncurr].Nm;
-        Sd = Sd / naverage[Ncurr].Nm;
+        Sr = Sr / naverage[Ncurr].Nmeasure;
+        Sd = Sd / naverage[Ncurr].Nmeasure;
         naverage[Ncurr].R = Sr;
         naverage[Ncurr].D = Sd;
 
-        /* update original measurement offsets */
-        for (k = 0; k < naverage[Ncurr].Nm; k++) {
-          m = naverage[Ncurr].offset + k;
+        /* update original measurement offsets for new detections */
+        for (k = Nfirst; k < naverage[Ncurr].Nmeasure; k++) {
+          m = naverage[Ncurr].measureOffset + k;
           M = mpointer[m].measure;
           measure[M].dR = 3600.0*(Sr - mpointer[m].R);
@@ -143 +186 @@
         /* update current reference star position */
         RD_to_XY (&X[i], &Y[i], Sr, Sd, &tcoords);
+# else
+        /* update original measurement offsets for new detections */
+        for (k = Nfirst; k < naverage[Ncurr].Nmeasure; k++) {
+          m = naverage[Ncurr].measureOffset + k;
+          M = mpointer[m].measure;
+          measure[M].dR = 3600.0*(naverage[Ncurr].R - mpointer[m].R);
+          measure[M].dD = 3600.0*(naverage[Ncurr].D - mpointer[m].D);
+        }
+# endif
       }
     }
@@ -148 +200 @@
     j = first_j;
   }
+
   if (Nmeas != Nmeasure) {
-    fprintf (stderr, "failure to match measures\n");
+    fprintf (stderr, "failure to match %d measures (%d of %d matched)\n", Nmeasure - Nmeas, Nmeas, Nmeasure);
   }
   
   /* create a new Measure table in the appropriate sequence */
-  ALLOCATE (nmeasure, Measure, Nmeasure);
-  for (i = 0; i < Nmeasure; i++) {
+  ALLOCATE (nmeasure, Measure, Nmeas);
+  for (i = 0; i < Nmeas; i++) {
     nmeasure[i]        = measure[mpointer[i].measure];
     nmeasure[i].averef = mpointer[i].averef;
@@ -160 +213 @@
 
   /* create an empty SecFilt table */
-  ALLOCATE (nsecfilt, SecFilt, Naves*catalog[0].Nsecfilt);
-  bzero (nsecfilt, Naves*catalog[0].Nsecfilt*sizeof(SecFilt));
+  // ALLOCATE (nsecfilt, SecFilt, Naves*catalog[0].Nsecfilt);
+  // bzero (nsecfilt, Naves*catalog[0].Nsecfilt*sizeof(SecFilt));
 
   free (catalog[0].average);
@@ -171 +224 @@
   catalog[0].secfilt = nsecfilt;
   
+  // allow output catalog to retain fewer measures
   catalog[0].Naverage = Naves;
+  catalog[0].Nmeasure = Nmeas;
   catalog[0].Nsecf_mem = Naves*catalog[0].Nsecfilt;
+  
+  return;
 }
 
@@ -185 +242 @@
    (measures are just moved around).  thus, catalog[0].Nmeasure does not change.
 */   
-   

trunk/Ohana/src/photdbc/src/make_subcatalog.c

@@ -24 +24 @@
 
     // exclude stars with too few measurements
-    if (NMEAS_MIN && (catalog[0].average[i].Nm < NMEAS_MIN)) continue; 
+    if (NMEAS_MIN && (catalog[0].average[i].Nmeasure < NMEAS_MIN)) continue; 
 
     /* assign average and secfilt values */
     subcatalog[0].average[Naverage] = catalog[0].average[i];
-    subcatalog[0].average[Naverage].offset = Nmeasure;
+    subcatalog[0].average[Naverage].measureOffset = Nmeasure;
     for (j = 0; j < Nsecfilt; j++) {
       subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j];
@@ -35 +35 @@
     minMag = 32;
     Nm = 0;
-    for (j = 0; j < catalog[0].average[i].Nm; j++) {
+    for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
 
-      offset = catalog[0].average[i].offset + j;
+      offset = catalog[0].average[i].measureOffset + j;
 
       # if 0
@@ -83 +83 @@
       continue; 
     }
-    subcatalog[0].average[Naverage].Nn = 0;
-    subcatalog[0].average[Naverage].Nm = Nm;
+    subcatalog[0].average[Naverage].Nmissing = 0;
+    subcatalog[0].average[Naverage].Nmeasure = Nm;
     Naverage ++;
     if (Naverage == NAVERAGE) {

trunk/Ohana/src/photdbc/src/photdbc.c

@@ -60 +60 @@
         
     // XXX add other filters here:
-    // join_stars (catalog);
+    join_stars (&outcatalog);
     // unique_measures (catalog);
     // flag_measures (&db, catalog);