Changeset 37116 for trunk/Ohana/src/relphot

Timestamp:

Jul 24, 2014, 5:10:56 PM (12 years ago)

Author:

eugene

Message:

merge relphot changes from eam_branches/ipp-20140717

Location:

trunk/Ohana

Files:

• . (modified) (1 prop)
• src/relphot/include/relphot.h (modified) (2 diffs)
• src/relphot/src/ImageOps.c (modified) (1 diff)
• src/relphot/src/StarOps.c (modified) (2 diffs)
• src/relphot/src/bcatalog.c (modified) (1 diff)
• src/relphot/src/setMrelCatalog.c (modified) (16 diffs)

trunk/Ohana/src/relphot/include/relphot.h

@@ -440 +440 @@
 int MatchImageName (off_t meas, int cat, char *name);
 int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID);
+int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID);
 
 int load_tess (char *treefile);
@@ -451 +452 @@
 int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
 int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found);
-int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell);
+int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure);
 int setGlobalObjStats (Average *average, Measure *measure);
 

trunk/Ohana/src/relphot/src/ImageOps.c

@@ -491 +491 @@
   if (projectID[i] != myProjectionID) return FALSE;
   if (skycellID[i] != mySkycellID) return FALSE;
+  
+  return TRUE;
+}
+
+// returns image.Mcal - ff(x,y)
+int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID) {
+
+  off_t i;
+
+  if (!MeasureToImage) return FALSE;
+
+  if (meas < 0) return FALSE;
+
+  i = MeasureToImage[cat][meas];
+  if (i == -1) return FALSE;
+
+  if (tessID[i]    == -1) return FALSE;
+  if (projectID[i] == -1) return FALSE;
+  if (skycellID[i] == -1) return FALSE;
+
+  *myTessID = tessID[i];
+  *myProjectionID = projectID[i];
+  *mySkycellID =  skycellID[i];
   
   return TRUE;

trunk/Ohana/src/relphot/src/StarOps.c

@@ -1086 +1086 @@
           if (ecode != thisCode) { continue; }
 
+          // SKIP gpc1 stack data
+          if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
+          
+          // SKIP gpc1 forced-warp data
+          if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
+          
           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
           Mcal  = getMcal  (m, i, flatcorr, catalog);
@@ -1131 +1137 @@
           if (ecode != thisCode) { continue; }
 
+          // SKIP gpc1 stack data
+          if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
+          
+          // SKIP gpc1 forced-warp data
+          if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
+          
           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
           Mcal  = getMcal  (m, i, flatcorr, catalog);

trunk/Ohana/src/relphot/src/bcatalog.c

@@ -78 +78 @@
         int Nsec = GetPhotcodeNsec(thisCode);
 
-        fprintf (stderr, "really use dvo_secfilt_init?");
+        // fprintf (stderr, "really use dvo_secfilt_init?");
         dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]);
 

trunk/Ohana/src/relphot/src/setMrelCatalog.c

@@ -50 +50 @@
   int isSetMrelFinal = (pass >= 0);
 
-  char *primaryCell = NULL;
-  if (isSetMrelFinal) {
-    ALLOCATE (primaryCell, char, DVO_MAX_PATH);
-  }
+  // XX char *primaryCell = NULL;
+  // XX if (isSetMrelFinal) {
+  // XX   ALLOCATE (primaryCell, char, DVO_MAX_PATH);
+  // XX }
 
   for (j = 0; j < catalog[Nc].Naverage; j++) {
@@ -64 +64 @@
     // only apply Stack operation on setMrelFinal in first pass 
     if (isSetMrelFinal && (pass == 0)) {
-      setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], primaryCell);
+      setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m]);
       setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
     }
@@ -74 +74 @@
     }
   }
-  if (primaryCell) free (primaryCell);
+  // if (primaryCell) free (primaryCell);
   return (TRUE);
 }
@@ -441 +441 @@
 
 // only apply Stack operation on setMrelFinal in first pass 
-int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell) {
-
-  off_t k, ID;
-
-  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
+// this function has 3 goals, not to be confused:
+// 1) find and mark the PRIMARY detections (regardless of the quality of the detection)
+// 2) select the BEST detections per filter (regardless of PRIMARY)
+// 3) apply the zero point and AB->Jy transformations
+int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure) {
+
+  off_t k;
+
+  float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0;
 
   // set the primary projection cell and skycell for this coordinate
@@ -467 +471 @@
     int haveStack = FALSE;
 
-    // need to find the measurement closest to the center of its skycell, as well as the
-    // closest for the subset of primary projection cells
-
-    float stackCenterOffsetMin = 1e9;
-    off_t stackCenterMeasureMin = -1;
-
-    float stackPrimaryOffsetMin = 1e9;
-    off_t stackPrimaryMeasureMin = -1;
+    float psfQFbest = 0.0;
+
+    off_t stackBestMeasure = -1;
+    off_t stackPrimaryMeasure = -1;
 
     int isBad = FALSE;
@@ -481 +481 @@
     off_t meas = measureOffset;
     for (k = 0; k < average[0].Nmeasure; k++, meas++) {
+
+      // only examine gpc1 stack data
+      if (!isGPC1stack(measure[k].photcode)) continue;
 
       // skip measurements that do not match the current photcode
@@ -489 +492 @@
       // clear this bit for all measurements
       measure[k].dbFlags &= ~ID_MEAS_STACK_PRIMARY;
-
-      // only examine gpc1 stack data
-      if (!isGPC1stack(measure[k].photcode)) continue;
+      measure[k].dbFlags &= ~ID_MEAS_STACK_PHOT_SRC;
 
       haveStack = TRUE;
       haveStackObject = TRUE;
       
+      // match measurement to its image (this is just a check, right?)
+      if (getImageEntry (meas, cat) < 0) {
+        fprintf (stderr, "*** warning : this should be impossible\n");
+        abort();
+        continue;
+      }
+
+      // find the PRIMARY measurement
+
+      // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the 
+      // primary skycell for this position
+      if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
+        stackPrimaryMeasure = k;
+        measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
+        secfilt[Nsec].stackPrmryOff = meas;
+      }
+
+      // now choose the BEST measurements (may also be PRIMARY)
+
       if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); 
 
-      if (getImageEntry (meas, cat) < 0) continue;
       // measurements without an image are either external reference photometry or
       // data for which the associated image has not been loaded (probably because of
@@ -510 +529 @@
       if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid);
 
-      Msys = PhotSys (&measure[k], &average[0], &secfilt[0], MAG_CLASS_PSF);
-      if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys);
-
-      unsigned int stackImageID;
-
-      // which stack image should we use for the mean value?
-      // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell
-      if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
-        float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
-        if (stackPrimaryOffset < stackPrimaryOffsetMin) {
-          stackPrimaryOffsetMin = stackPrimaryOffset;
-          // stackPrimaryIDmin = stackImageID;
-          stackPrimaryMeasureMin = k;
-        }
-      }
-
-      // get the center distance for the generic case:
-      float stackCenterOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
-      if (stackCenterOffset < stackCenterOffsetMin) {
-        stackCenterOffsetMin = stackCenterOffset;
-        // stackCenterIDmin = stackImageID;
-        stackCenterMeasureMin = k;
+      // NOTE: negative and insignificant vlues values are allowed, but not NAN flux values
+      Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF);
+      if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst);
+
+      if (measure[k].psfQFperf > psfQFbest) {
+        psfQFbest = measure[k].psfQFperf;
+        stackBestMeasure = k;
       }
 
@@ -554 +558 @@
 
     // measurements which are bad will not have a valid stack entry and are skipped
-    k = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
+    k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure;
     if (k < 0) continue;
-
 
     // we are now populating stackDetectID not stack Image ID in secfilt
     // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
-    ID = measure[k].extID; // for the stack, this is the stackDetectID
+    // ID = measure[k].extID; // for the stack, this is the stackDetectID
 
     // get the zero point for the selected image
@@ -605 +608 @@
     secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
 
+    // Jy to AB mags
     secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN;
     secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
     secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN;
 
-    secfilt[Nsec].stackDetectID = ID;
-
+    secfilt[Nsec].stackBestOff = k + measureOffset;
+
+    // this is the measurement used by secfilt[]
     measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
-    if (stackPrimaryMeasureMin >= 0) {
-      measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
+    if (k == stackPrimaryMeasure) {
       secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
     }
@@ -642 +646 @@
 
   off_t k;
-  float Fsys = 0, Mcal= 0;
+  float Mcal= 0;
 
   // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
@@ -664 +668 @@
 
   // option for a test print
-  if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
+  if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) {
     fprintf (stderr, "test obj\n");
     print_measure_set_alt (average, secfilt, measure);
@@ -680 +684 @@
     if (found[Nsec]) continue;  
       
+    // I need to restrict the measurements coming only from the same skycell as the secfilt values
+    // if we cannot make that association, skip this set of warps
+    int stkTessID, stkProjID, stkSkycellID;
+    if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue;
+
     off_t meas = measureOffset;
 
@@ -694 +703 @@
       if (!code) continue;
       if (code->equiv != thisCode) { continue; }
+
+      measure[k].dbFlags &= ~ID_MEAS_WARP_USED;
+
+      if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue;
 
       if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
@@ -720 +733 @@
       // are small (and not measured)
       float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
-      if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys);
+      if (isnan(Fpsf)) SKIP_THIS_MEAS(Nsys);
       // if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
       // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
 
       float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
+      if (isnan(dFpsf)) SKIP_THIS_MEAS (Nsys);
+
+      int isBad = (measure[k].psfQF < 0.85);
+      isBad |= isnan(measure[k].psfQF);
+      if (isBad) SKIP_THIS_MEAS (Nsys);
+
+      int isPoor = (measure[k].psfQFperf < 0.85);
+      if ((pass == 0) && isPoor) SKIP_THIS_MEAS (Nsys);
+
       dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error
       Fpsflist[Npsf] = Fpsf * Fcal;
@@ -730 +752 @@
       wpsflist[Npsf] = 1.0;
       Npsf ++;
+
+      measure[k].dbFlags |= ID_MEAS_WARP_USED;
 
       float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);