Changeset 4115 for trunk/psphot/src/pspsf.c

Timestamp:

Jun 4, 2005, 5:22:13 PM (21 years ago)

Author:

eugene

Message:

running sources in brightness order, using a mask

File:

• trunk/psphot/src/pspsf.c (modified) (8 diffs)

trunk/psphot/src/pspsf.c

@@ -36 +36 @@
     test->metricStats = NULL;
     test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
 
     for (int i = 0; i < test->modelFLT->n; i++) {
         test->modelFLT->data[i] = NULL;
         test->modelPSF->data[i] = NULL;
+        test->mask->data.U8[i]  = 0;
     }   
     return (test);
@@ -57 +59 @@
         psModel  *model  = pmSourceModelGuess (source, test->modelType); 
 
-        // fit model as FLT, not PSF (skip poor fits)
-        if (!pmSourceFitModel (source, model, false)) continue;
+        // fit model as FLT, not PSF (mask & skip poor fits)
+        if (!pmSourceFitModel (source, model, false)) {
+          test->mask->data.U8[i] = 1;
+          continue;
+        }
         test->modelFLT->data[i] = model;
         Nflt ++;
@@ -66 +71 @@
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    pmPSFFromModels (test->psf, test->modelFLT);
+    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
 
     // stage 3: refit with fixed shape parameters
@@ -73 +78 @@
         psSource *source = test->sources->data[i];
         psModel  *modelFLT = test->modelFLT->data[i];
-        if (modelFLT == NULL) continue;
+
+        // masked for: bad model fit, outlier in parameters
+        if (test->mask->data.U8[i]) continue;
+        if (modelFLT == NULL) {
+          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
+          continue;
+        }
         psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
 
         // fit model as PSF, not FLT (skip poor fits)
-        if (!pmSourceFitModel (source, modelPSF, true)) continue;
+        if (!pmSourceFitModel (source, modelPSF, true)) {
+          test->mask->data.U8[i] = 1;
+          continue;
+        }
         test->modelPSF->data[i] = modelPSF;
         Npsf ++;
@@ -92 +106 @@
         float fitSum = 0;
       
+        // is this metricMask redundant with test->mask ?
         metricMask->data.U8[i] = 1;
         test->metric->data.F64[i] = 0;
 
+        if (test->mask->data.U8[i]) continue;
+
         psModel *model = test->modelPSF->data[i];
-        if (model == NULL) continue;
+        if (model == NULL) {
+          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
+          continue;
+        }
 
         psImage *image = ((psSource *)test->sources->data[i])->pixels;
         psImage *mask  = ((psSource *)test->sources->data[i])->mask;
 
+        // this metric is Ap-Fit
         float sky = model->params->data.F32[0];
-
         for (int ix = 0; ix < image->numCols; ix++) {
             for (int iy = 0; iy < image->numRows; iy++) {
@@ -110 +130 @@
             }
         }
-        // this metric is Ap-Fit
         // fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
+        // keep the good metrics
         if ((fitSum > 0) && (obsSum > 0)) {
             metricMask->data.U8[i] = 0;
             test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
-        }       
+        } else {
+          test->mask->data.U8[i] = 1;
+        }
     }
 
@@ -128 +150 @@
 // input: an array of psModels, pre-allocated psf
 // some of the array entries may be NULL, ignore them
-bool pmPSFFromModels (pmPSF *psf, psArray *models) {
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
 
     int n;
 
-    // construct the x and y vectors from the collection of objects
-    // count the number of valid model measurements
+    // construct the fit vectors from the collection of objects
+    // use the mask to ignore missing fits 
     psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
     psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
-    n = 0;
+    psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
+
     for (int i = 0; i < models->n; i++) {
         psModel *model = models->data[i];
@@ -142 +166 @@
 
         // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
-        x->data.F64[n] = model->params->data.F32[2];
-        y->data.F64[n] = model->params->data.F32[3];
-        n++;
-    }
-    x->n = y->n = n;
+        x->data.F64[i] = model->params->data.F32[2];
+        y->data.F64[i] = model->params->data.F32[3];
+    }
 
     // we are doing a robust fit.  after each pass, we drop points which are 
-    // more deviant than three sigma.
-    psVector *z = psVectorAlloc (x->n, PS_TYPE_F64);
-    psVector *dz = psVectorAlloc (x->n, PS_TYPE_F64);
+    // more deviant than three sigma. 
+    // mask is currently updated for each pass. this is inconsistent
+
     for (int i = 0; i < psf->params->n; i++) {
-        n = 0;
         for (int j = 0; j < models->n; j++) {
             psModel *model = models->data[j];
             if (model == NULL) continue;
-            z->data.F64[n] = model->params->data.F32[i + 4];
-            dz->data.F64[n] = 1;
-            n++;
+            z->data.F64[j] = model->params->data.F32[i + 4];
+            dz->data.F64[j] = 1;
             // XXX EAM : need to use actual errors?
             // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
         }
-        if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
-        z->n = n;
-        dz->n = n;
-
-        psf->params->data[i] = RobustFit2D (psf->params->data[i], x, y, z, dz);
+        // if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
+
+        psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
         // psPolynomial2DDump (psf->params->data[i]);
     }