Changeset 33652 for trunk/Ohana/src/relastro/src/UpdateObjects.c

Timestamp:

Apr 1, 2012, 3:00:19 PM (14 years ago)

Author:

eugene

Message:

updates to relastro to support distributed dvo; some optmimization (eg, image matching to area); cleanup names of static helper arrays in ImageOps, MosaicOps; cleanup rules for -reset options (relastro does not reset photom bits); skip objects with NAN error bars (can be used to exclude some photcodes from astrometry); enum for mask bits in FitChip and related; choose max order using same rules as psastro; plug leaks; reduce verbosity; big re-org of main (for parallel code); only use accepted measurements for Tmin,Tmax,Trange,Tmean; require a min parallax factor to fit parallaxes; use MeasureTiny where possible; use new ParsePhotcodeList function; testparallax function; threaded version of UpdateChips; bcatalog limit density selecting objects by number of measurements; high-speed now copies support tables to output database; high-speed may write to a parallel database; test for unsorted database (& fail if true)

File:

• trunk/Ohana/src/relastro/src/UpdateObjects.c (modified) (9 diffs)

trunk/Ohana/src/relastro/src/UpdateObjects.c

@@ -38 +38 @@
 }  
 
+// This function operates on both Measure and MeasureTiny.  In the big stages, this should
+// be called with just MeasureTiny set and Measure == NULL
 int UpdateObjects (Catalog *catalog, int Ncatalog) {
 
@@ -72 +74 @@
   T2000 = ohana_date_to_sec ("2000/01/01");
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
+
   Nsecfilt = GetPhotcodeNsecfilt ();
-  assert (catalog[0].Nsecfilt == Nsecfilt);
+  if (Ncatalog) {
+    assert (catalog[0].Nsecfilt == Nsecfilt);
+  }
 
   NaveSum = NparSum = NpmSum = NskipSum = 0;
   for (i = 0; i < Ncatalog; i++) {
 
-    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
+    if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
 
     Nave = Npar = Npm = Nskip = 0;
@@ -86 +91 @@
       XVERB = FALSE;
 
-      // skip objects which are known to be problematic
-      // XXX include this code or not?
-# if (0)
-      if (catalog[i].average[j].code & STAR_BAD) {
-        Nskip ++;
-        continue;  
-      }
-# endif
-
       if (catalog[i].average[j].Nmeasure == 0) {
           continue;
@@ -101 +97 @@
       N = 0;
       m = catalog[i].average[j].measureOffset;
-      Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25);
+      MeasureTiny *measure = &catalog[i].measureT[m];
+      Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL;
+      // when we update the output measure values, we need to do it here
+
       mode = FIT_MODE;
 
       // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
-      Tmean = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-
-        //does the measurement pass the supplied filtering constraints?
-        if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) {
-          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+
+        // does the measurement pass the supplied filtering constraints?
+        if (!MeasFilterTestTiny(&measure[k], FALSE)) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
           continue;
         }
 
         //outlier rejection
-        if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
-          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+        if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
           continue;
         }
 
-        // exclude measurements by previous outlier detection
-        // XXX include this code or not?
-# if (0)
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) { 
-          catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
+        measure[k].dbFlags |= ID_MEAS_USED_OBJ;
+        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
+
+        R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        D[N] = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+
+        // XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects
+        // to get Time in years relative to UNIX Tzero (1970/01/01)
+        T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
+
+        // dX, dY : error in arcsec -- 
+        dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
+        dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
+
+        // allow a given photcode or measurement to be
+        // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
+        if (isnan(dX[N])) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
           continue;
         }
-# endif
-
-        catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
-
-        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-        D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-        T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
-
-        Tmin = MIN(Tmin, T[N]);
-        Tmax = MAX(Tmax, T[N]);
-        Tmean += T[N];
-
-        // dX, dY : error in arcsec -- 
-        dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
-        dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
+        if (isnan(dY[N])) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          continue;
+        }
 
         // add systematic error in quadrature, if desired
         // only do this after the fit has converged (or you will never improve the poor images)
         // if (INCLUDE_SYS_ERR) {
-        // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys);
+        // float dRsys = FromShortPixels(measure[k].dRsys);
         // dX[N] = hypot(dX[N], dRsys);
         // dY[N] = hypot(dY[N], dRsys);
@@ -154 +156 @@
         // dY[N] = 0.1;
 
-        dT[N] = catalog[i].measure[m].dt;
+        dT[N] = measure[k].dt;
 
         // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
@@ -172 +174 @@
       catalog[i].average[j].flags &= ~ID_STAR_FEW;
 
+      // find Tmin & Tmax from the list of accepted measurements
+      Tmean = 0;
+      Tmin = Tmax = T[0];
+      for (k = 0; k < N; k++) {
+        Tmin = MIN(Tmin, T[k]);
+        Tmax = MAX(Tmax, T[k]);
+        Tmean += T[k];
+      }
       // XXX add the parallax factor range as a criterion as well
       Trange = Tmax - Tmin;
       if (Trange < PM_DT_MIN) mode = FIT_AVERAGE;
       if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;
+
+      if (FIT_TARGET == TARGET_HIGH_SPEED) {
+          Tmean = 0.5*(Tmax - Tmin);
+      } else {
+          Tmean /= (float) N;
+      }
 
       // too few measurements for average position (require 2 values)
@@ -192 +208 @@
       coords.crval2 = D[0];
 
-      if (FIT_TARGET == TARGET_HIGH_SPEED) {
-          Tmean = 0.5*(Tmax - Tmin);
-      } else {
-          Tmean /= (float) N;
-      }
-      
       // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e);
       XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
       XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
-      // XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
 
       // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
@@ -232 +241 @@
         // fprintf (stderr, "parallax fitting is still untested (%s, %d)\n", __FILE__, __LINE__);
 
+        float pXmin = +2.0;
+        float pXmax = -2.0;
+        float pYmin = +2.0;
+        float pYmax = -2.0;
         for (k = 0; k < N; k++) {
           ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
-        }
-        FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
-        XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
-        catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
-        Npar ++;
+          pXmin = MIN (pXmin, pX[k]);
+          pXmax = MAX (pXmax, pX[k]);
+          pYmin = MIN (pYmin, pY[k]);
+          pYmax = MAX (pYmax, pY[k]);
+        }
+        float dXRange = pXmax - pXmin;
+        float dYRange = pYmax - pYmin;
+        float parRange = hypot (dXRange, dYRange);
+        
+        if (parRange >= PAR_FACTOR_MIN) {
+          FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
+          XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
+          catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
+          Npar ++;
+        }
       }   
 
@@ -307 +330 @@
 
       // the measure fields must be updated before the average fields
-      m = catalog[i].average[j].measureOffset;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-        setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+        setMeanR (fit.Ro, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        setMeanD (fit.Do, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        if (measureBig) {
+          setMeanR_Big (fit.Ro, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+          setMeanD_Big (fit.Do, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        }
       }