Changeset 36796

Timestamp:

Jun 2, 2014, 6:07:00 AM (12 years ago)

Author:

eugene

Message:

add basic bogus-detection code

Location:

branches/eam_branches/ipp-20140423/Ohana/src/relastro

Files:

• include/relastro.h (modified) (2 diffs)
• src/args.c (modified) (2 diffs)
• src/bcatalog.c (modified) (5 diffs)
• src/dvo_astrom_ops.c (modified) (2 diffs)
• src/load_catalogs.c (modified) (2 diffs)

branches/eam_branches/ipp-20140423/Ohana/src/relastro/include/relastro.h

@@ -233 +233 @@
 int    CLIP_THRESH;
 int USE_BASIC_CHECK;
+
+int ExcludeBogus;
+double ExcludeBogusRadius;
 
 FitMode FIT_MODE;
@@ -417 +420 @@
 int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
 int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
+double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt);
+double getMeanD_Big (Measure *measure, Average *average, SecFilt *secfilt);
 int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
 int setMeanD_Big (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);

branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/args.c

@@ -183 +183 @@
     remove_argument (N, &argc, argv);
     FlagOutlier = TRUE;
+  }
+
+  ExcludeBogus = FALSE;
+  ExcludeBogusRadius = 0.0;
+  if ((N = get_argument (argc, argv, "-exclude-bogus"))) {
+    remove_argument (N, &argc, argv);
+    ExcludeBogusRadius = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ExcludeBogus = TRUE;
   }
 
@@ -593 +602 @@
     remove_argument (N, &argc, argv);
     FlagOutlier = TRUE;
+  }
+
+  ExcludeBogus = FALSE;
+  ExcludeBogusRadius = 0.0;
+  if ((N = get_argument (argc, argv, "-exclude-bogus"))) {
+    remove_argument (N, &argc, argv);
+    ExcludeBogusRadius = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ExcludeBogus = TRUE;
   }
 

branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/bcatalog.c

@@ -6 +6 @@
 static int Nskip2 = 0;
 
+static int NskipBogus = 0;
+
 // test image: 2013/06/15,13:25:51, GPC1.r.XY50
 static int CHECK_TEST_IMAGE = TRUE;
@@ -21 +23 @@
   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
   int Nsecfilt;
+  Coords coords;
+
+  /* for outlier rejection, project coordinates to a plane centered on the object with units of arcsec */
+  coords.crval1 = 0;
+  coords.crval2 = 0;
+  coords.crpix1 = 0;
+  coords.crpix2 = 0;
+  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
+  coords.pc1_1  = coords.pc2_2 = 1.0;
+  coords.pc1_2  = coords.pc2_1 = 0.0;
+  coords.Npolyterms = 1;
+  strcpy (coords.ctype, "DEC--SIN");
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -79 +93 @@
       }
 
-      // filter out outliers
+      // filter out outliers - these are detections inconsistent with the offset distribution
       if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
         catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
@@ -88 +102 @@
       }
       catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP;
+
+      // exclude bogus 
+      if (ExcludeBogus) {
+          double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
+          double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
+          coords.crval1 = catalog[0].average[i].R;
+          coords.crval2 = catalog[0].average[i].D;
+          double Xi, Yi;
+          RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
+          double radius = hypot(Xi, Yi);
+          // XXX this radius needs to be configurable
+          if (radius > ExcludeBogusRadius) {
+              NskipBogus ++;
+              continue;
+          }
+      }
 
       // re-assess on each run of relastro if a measurement should be used
@@ -159 +189 @@
 
 void bcatalog_show_skips () {
-  fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
-  fprintf (stderr, "Nkeep: %d, %d\n", Nkeep1, Nkeep2);
+    fprintf (stderr, "NskipBogus: %d\n", NskipBogus);
+    // fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
+    // fprintf (stderr, "Nkeep: %d, %d\n", Nkeep1, Nkeep2);
 }
 

branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/dvo_astrom_ops.c

@@ -91 +91 @@
 }
 
+double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt) {
+
+  double ra;
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  ra = average[0].R - measure[0].dR / 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (ra);
+}
+
+double getMeanD_Big (Measure *measure, Average *average, SecFilt *secfilt) {
+
+  double dec;
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  dec = average[0].D - measure[0].dD / 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (dec);
+}
+
 int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
 
@@ -144 +184 @@
   return (TRUE);
 }
+

branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/load_catalogs.c

@@ -84 +84 @@
   }
 
-  // XXX TEST : bcatalog_show_skips();
+  bcatalog_show_skips();
 
   Nstar = 0;
@@ -162 +162 @@
     
     if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);       strcpy (command, tmpline); }
+    if (ExcludeBogus)  { snprintf (tmpline, 1024, "%s -exclude-bogus %f", command, ExcludeBogusRadius); strcpy (command, tmpline); }
     
     if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command);     strcpy (command, tmpline); }