Changeset 21173 for trunk/psModules/src/objects/pmPSFtry.c

Timestamp:

Jan 26, 2009, 2:01:33 PM (17 years ago)

Author:

eugene

Message:

keep input mask values, pass along to pmTrend functions; detect a few failure modes

File:

• trunk/psModules/src/objects/pmPSFtry.c (modified) (10 diffs)

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.67 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-12-08 02:51:14 $
+ *  @version $Revision: 1.68 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-01-27 00:01:33 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -175 +175 @@
 
         pmSource *source = psfTry->sources->data[i];
+        if (!source->moments) {
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
+            continue;
+        }
+        if (!source->moments->nPixels) {
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
+            continue;
+        }
+
         source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
         if (source->modelEXT == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "failed to build model");
-            psFree(psfTry);
-            return NULL;
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
+            continue;
         }
 
@@ -591 +600 @@
     psVector *e2   = psVectorAlloc (sources->n, PS_TYPE_F32);
     psVector *mag  = psVectorAlloc (sources->n, PS_TYPE_F32);
+
     for (int i = 0; i < sources->n; i++) {
+        // skip any masked sources (failed to fit one of the model steps or get a magnitude)
+        if (srcMask->data.U8[i]) continue;
+        
         pmSource *source = sources->data[i];
-        if (source->modelEXT == NULL) continue;
-        // XXX I am relying on the fact that none of the masked sources
-        // have modelEXT set here.  perhaps use the value of psfTry->mask instead?
+        assert (source->modelEXT); // all unmasked sources should have modelEXT
 
         psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
@@ -619 +630 @@
         for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
 
-            psVectorInit (mask, 0);
+            // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
+            for (int i = 0; i < mask->n; i++) {
+                mask->data.U8[i] = srcMask->data.U8[i];
+            }
             if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
                 break;
@@ -635 +649 @@
         }
 
-        // XXX supply the resulting mask values back to srcMask
-        psVectorInit (mask, 0);
+        // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
+        for (int i = 0; i < mask->n; i++) {
+            mask->data.U8[i] = srcMask->data.U8[i];
+        }
         if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
             psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
@@ -836 +852 @@
     // the mask marks the values not used to calculate the ApTrend
     psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_U8);
-    psVectorInit (fitMask, 0);
+    // copy mask values to fitMask as a starting point
+    for (int i = 0; i < fitMask->n; i++) {
+        fitMask->data.U8[i] = mask->data.U8[i];
+    }
 
     // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
@@ -884 +903 @@
 
     // construct the fitted values and the residuals
-    psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], x_tst, y_tst);
-    psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], x_tst, y_tst);
-    psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], x_tst, y_tst);
+    psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst);
+    psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst);
+    psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst);
 
     psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit);
@@ -895 +914 @@
     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10);
     // psTraceSetLevel ("psLib.math.vectorClippedStats", 10);
-    pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, psStatsStdevOption(psf->psfTrendStats->options));
+    pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options));
     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0);
     // psTraceSetLevel ("psLib.math.vectorClippedStats", 0);
@@ -932 +951 @@
 
 // calculate the scatter of the parameters
-bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psStatsOptions stdevOpt)
+bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psMaskType maskValue, psStatsOptions stdevOpt)
 {
 
@@ -941 +960 @@
     float dEsquare = 0.0;
     psStatsInit (stats);
-    psVectorStats (stats, e0res, NULL, NULL, 0xff);
+    psVectorStats (stats, e0res, NULL, mask, maskValue);
     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
 
     psStatsInit (stats);
-    psVectorStats (stats, e1res, NULL, NULL, 0xff);
+    psVectorStats (stats, e1res, NULL, mask, maskValue);
     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
 
     psStatsInit (stats);
-    psVectorStats (stats, e2res, NULL, NULL, 0xff);
+    psVectorStats (stats, e2res, NULL, mask, maskValue);
     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));