Changeset 30031 for trunk/psModules/src/objects

Timestamp:

Dec 14, 2010, 9:29:23 AM (15 years ago)

Author:

eugene

Message:

modifications to psImageMapFit to allow for unfittable data (at a given resolution)

Location:

trunk/psModules/src/objects

Files:

• . (modified) (1 prop)
• pmPSF_IO.c (modified) (3 diffs)
• pmPSFtry.h (modified) (2 diffs)
• pmPSFtryMakePSF.c (modified) (9 diffs)
• pmPSFtryModel.c (modified) (2 diffs)
• pmTrend2D.c (modified) (4 diffs)
• pmTrend2D.h (modified) (1 diff)

trunk/psModules/src/objects/pmPSF_IO.c

@@ -62 +62 @@
 #include "pmSourceIO.h"
 
+bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header);
+
 bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file)
 {
@@ -851 +853 @@
 
     // read the psf clump data for each region
+    status = false;
     if (roAnalysis) {
-        int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN");
-        if (!status) {
-            // read old-style psf clump data
-
-            char regionName[64];
-            snprintf (regionName, 64, "PSF.CLUMP.REGION.000");
-            psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
-
-            psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
-            if (!regionMD) {
-                regionMD = psMetadataAlloc();
-                psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
-                psFree (regionMD);
-            }
-
-            // psf clump data
-            pmPSFClump psfClump;
-
-            psfClump.X  = psMetadataLookupF32 (&status, header, "PSF_CLX" );  assert(status);
-            psfClump.Y  = psMetadataLookupF32 (&status, header, "PSF_CLY" );  assert(status);
-            psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX");  assert(status);
-            psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY");  assert(status);
-
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
-        } else {
-            psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
-
-            for (int i = 0; i < nRegions; i++) {
-                char key[10];
-                char regionName[64];
-                snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
-
-                psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
-                if (!regionMD) {
-                    regionMD = psMetadataAlloc();
-                    psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
-                    psFree (regionMD);
-                }
-
-                // psf clump data
-                pmPSFClump psfClump;
-
-                snprintf (key, 9, "CLX_%03d", i);
-                psfClump.X  = psMetadataLookupF32 (&status, header, key);  assert(status);
-                snprintf (key, 9, "CLY_%03d", i);
-                psfClump.Y  = psMetadataLookupF32 (&status, header, key);  assert(status);
-                snprintf (key, 9, "CLDX_%03d", i);
-                psfClump.dX = psMetadataLookupF32 (&status, header, key);  assert(status);
-                snprintf (key, 9, "CLDY_%03d", i);
-                psfClump.dY = psMetadataLookupF32 (&status, header, key);  assert(status);
-
-                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
-                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
-                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
-                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
-            }
-        }
+        status = pmPSFmodelReadPSFClump (roAnalysis, header);
+        if (!status) {
+            psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 0);
+        }
+    } 
+    if (!roAnalysis || !status) {
+        psWarning ("no PSF.CLUMP data available for PSF model");
     }
 
@@ -1123 +1073 @@
 }
 
-// XXX pmPSF to/from Metadata functions were defined for 1.22 and earlier, but were dropped
+bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header) {
+
+    bool status = false;;
+
+    int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN");
+    if (!status) {
+        // read old-style psf clump data
+
+        char regionName[64];
+        snprintf (regionName, 64, "PSF.CLUMP.REGION.000");
+        psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+
+        if (!regionMD) {
+            regionMD = psMetadataAlloc();
+            psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
+            psFree (regionMD);
+        }
+
+        // psf clump data
+        pmPSFClump psfClump;
+        psfClump.X  = psMetadataLookupF32 (&status, header, "PSF_CLX" );  if (!status) return false;
+        psfClump.Y  = psMetadataLookupF32 (&status, header, "PSF_CLY" );  if (!status) return false;
+        psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX");  if (!status) return false;
+        psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY");  if (!status) return false;
+
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
+        psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
+    } else {
+        for (int i = 0; i < nRegions; i++) {
+            char key[10];
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
+
+            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+            if (!regionMD) {
+                regionMD = psMetadataAlloc();
+                psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
+                psFree (regionMD);
+            }
+
+            // psf clump data
+            pmPSFClump psfClump;
+
+            snprintf (key, 9, "CLX_%03d", i);
+            psfClump.X  = psMetadataLookupF32 (&status, header, key);  if (!status) return false;
+            snprintf (key, 9, "CLY_%03d", i);
+            psfClump.Y  = psMetadataLookupF32 (&status, header, key);  if (!status) return false;
+            snprintf (key, 9, "CLDX_%03d", i);
+            psfClump.dX = psMetadataLookupF32 (&status, header, key);  if (!status) return false;
+            snprintf (key, 9, "CLDY_%03d", i);
+            psfClump.dY = psMetadataLookupF32 (&status, header, key);  if (!status) return false;
+
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
+        }
+        psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
+    }
+    return true;
+}

trunk/psModules/src/objects/pmPSFtry.h

@@ -100 +100 @@
 bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
 
-bool pmPSFtryMakePSF (pmPSFtry *psfTry);
+bool pmPSFtryMakePSF (bool *goodFit, pmPSFtry *psfTry);
 
 bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
@@ -123 +123 @@
 );
 
-bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask);
+bool pmPSFFitShapeParams (bool *goodFit, pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask);
 
 float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);

trunk/psModules/src/objects/pmPSFtryMakePSF.c

@@ -50 +50 @@
 Note: some of the array entries may be NULL (failed fits); ignore them.
  *****************************************************************************/
-bool pmPSFtryMakePSF (pmPSFtry *psfTry)
+bool pmPSFtryMakePSF (bool *goodFit, pmPSFtry *psfTry)
 {
     PS_ASSERT_PTR_NON_NULL(psfTry, false);
@@ -74 +74 @@
 
     // fit the shape parameters (SXX, SYY, SXY) as a function of position
-    if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, srcMask)) {
+    if (!pmPSFFitShapeParams (goodFit, psf, psfTry->sources, x, y, srcMask)) {
         psFree(x);
         psFree(y);
         return false;
+    }
+    if (!goodFit) {
+        psWarning ("poor fit to PSF shape parameters for trend order %d, %d, skipping\n", psf->trendNx, psf->trendNy);
+        psFree(x);
+        psFree(y);
+        return true;
     }
 
@@ -115 +121 @@
         // the mask is carried from previous steps and updated with this operation
         // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
-        if (!pmTrend2DFit (trend, srcMask, 0xff, x, y, z, NULL)) {
+        if (!pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, z, NULL)) {
             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
             psFree(x);
@@ -122 +128 @@
             return false;
         }
+        if (!goodFit) {
+            // if we do not get a good fit (but do not actually hit an error), 
+            // tell the calling program to try something else
+            psWarning ("poor fit to PSF parameter %d for trend order %d, %d, skipping\n", i, psf->trendNx, psf->trendNy);
+            psFree(x);
+            psFree(y);
+            psFree(z);
+            return true;
+        }
         if (trend->mode == PM_TREND_MAP) {
             // p_psImagePrint (2, trend->map->map, "param N Before"); // XXX TEST:
@@ -163 +178 @@
 
 // fit the shape parameters using the supplied order (pmPSF->trendNx,trendNy)
-bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
+bool pmPSFFitShapeParams (bool *goodFit, pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
 
     // we are doing a robust fit.  after each pass, we drop points which are more deviant than
@@ -219 +234 @@
         trend = psf->params->data[PM_PAR_E0];
         trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here
-        status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
+        status &= pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, e0, NULL);
+        if (!goodFit) {
+            psFree (e0);
+            psFree (e1);
+            psFree (e2);
+            return true;
+        }
         mean = psStatsGetValue (trend->stats, meanOption);
         stdev = psStatsGetValue (trend->stats, stdevOption);
@@ -228 +249 @@
         trend = psf->params->data[PM_PAR_E1];
         trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here
-        status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
+        status &= pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, e1, NULL);
+        if (!goodFit) {
+            psFree (e0);
+            psFree (e1);
+            psFree (e2);
+            return true;
+        }
         mean = psStatsGetValue (trend->stats, meanOption);
         stdev = psStatsGetValue (trend->stats, stdevOption);
@@ -237 +264 @@
         trend = psf->params->data[PM_PAR_E2];
         trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here
-        status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
+        status &= pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, e2, NULL);
+        if (!goodFit) {
+            psFree (e0);
+            psFree (e1);
+            psFree (e2);
+            return true;
+        }
         mean = psStatsGetValue (trend->stats, meanOption);
         stdev = psStatsGetValue (trend->stats, stdevOption);
@@ -246 +279 @@
         if (!status) {
             psError (PS_ERR_UNKNOWN, true, "failed to fit PSF shape params");
+            psFree (e0);
+            psFree (e1);
+            psFree (e2);
             return false;
         }

trunk/psModules/src/objects/pmPSFtryModel.c

@@ -136 +136 @@
 
         // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
-        if (!pmPSFtryMakePSF (psfTry)) {
+        bool goodFit = false;
+        if (!pmPSFtryMakePSF (&goodFit, psfTry)) {
             psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
             psFree(psfTry);
             return NULL;
         }
+        if (!goodFit) {
+            psWarning ("poor psf fit for order %d, skipping\n", i);
+            continue;
+        }
 
         // stage 3: refit with fixed shape parameters, measure pmPSFtry->metric
@@ -169 +174 @@
     }
     psFree (srcMask);
+
+    if (!minPSF) {
+        psError(PS_ERR_UNKNOWN, false, "failed to construct a valid psf model from the sources");
+        psFree(psfTry);
+        return NULL;
+    }
 
     // keep the ones matching the min systematic error:

trunk/psModules/src/objects/pmTrend2D.c

@@ -179 +179 @@
 }
 
-bool pmTrend2DFit(pmTrend2D *trend, psVector *mask, psVectorMaskType maskVal, const psVector *x,
+bool pmTrend2DFit(bool *goodFit, pmTrend2D *trend, psVector *mask, psVectorMaskType maskVal, const psVector *x,
                   const psVector *y, const psVector *f, const psVector *df)
 {
@@ -188 +188 @@
     PS_ASSERT_VECTOR_NON_NULL(f, false);
 
-    bool status;
+    bool status = false;
+    *goodFit = false;
+    // for the psImageMap fit, it is possible to have valid data but no valid solution for
+    // example, an isolated cell may not be reached from other cells, making the solution
+    // degenerate.  psImageMapFit should probably handle this case, but until it does, we allow
+    // it to fail on the result, but not yield an error (goodFit = false).
+    // psVectorClipFitPolynomial2D can not fail in this way (really?), so goodFit is always
+    // true
+
     switch (trend->mode) {
       case PM_TREND_POLY_ORD:
@@ -196 +204 @@
         // of points in the image, and potentially based on the fractional range of the
         // data?
+        *goodFit = true;
         break;
 
@@ -201 +210 @@
         // XXX supply fraction from trend elements
         // XXX need to add the API which adjusts the scale
-        status = psImageMapClipFit(trend->map, trend->stats, mask, maskVal, x, y, f, df);
+        status = psImageMapClipFit(goodFit, trend->map, trend->stats, mask, maskVal, x, y, f, df);
         break;
 

trunk/psModules/src/objects/pmTrend2D.h

@@ -76 +76 @@
     );
 
-bool pmTrend2DFit(pmTrend2D *trend,
+bool pmTrend2DFit(bool *goodFit,
+                  pmTrend2D *trend,
                   psVector *mask,       // Warning: mask is modified!
                   psVectorMaskType maskVal,