Changeset 6481 for trunk/psphot/src/psphotEnsemblePSF.c

Timestamp:

Feb 23, 2006, 6:16:11 PM (20 years ago)

Author:

eugene

Message:

fixed some errors with double sources, added second-pass linear PSF fit, multiple-depths

File:

• trunk/psphot/src/psphotEnsemblePSF.c (modified) (12 diffs)

trunk/psphot/src/psphotEnsemblePSF.c

@@ -1 +1 @@
 # include "psphot.h"
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);
 
 // 2006.02.07 : no leaks!
 // fit all reasonable sources with the linear PSF model
-bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) { 
+bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final) { 
 
     bool  status;
@@ -9 +10 @@
     float y;
     float f;
+    float r;
+
+    // psRegion allArray = psRegionSet (0, 0, 0, 0);
 
     psTimerStart ("psphot");
@@ -29 +33 @@
 
     // option to limit analysis to a specific region
-    bool UseAnalysisRegion = false;
-    psRegion AnalysisRegion = {0, 0, 0, 0};
     char *region = psMetadataLookupStr (&status, config, "ANALYSIS_REGION");
-    if (status) {
-        UseAnalysisRegion = true;
-        AnalysisRegion = psRegionFromString (region);
-        psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n", 
-                  AnalysisRegion.x0, AnalysisRegion.y0, 
-                  AnalysisRegion.x1, AnalysisRegion.y1);
-    }
+    psRegion AnalysisRegion = psRegionFromString (region);
+    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
+    // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
 
     for (int i = 0; i < sources->n; i++) {
@@ -45 +44 @@
         // skip non-astronomical objects (very likely defects)
         // XXX EAM : should we try these anyway?
-        if (inSource->mode &  PM_SOURCE_BLEND) continue;
         if (inSource->type == PM_SOURCE_DEFECT) continue; 
         if (inSource->type == PM_SOURCE_SATURATED) continue;
-
-        if (UseAnalysisRegion) {
-            if (inSource->moments->x < AnalysisRegion.x0) continue;
-            if (inSource->moments->y < AnalysisRegion.y0) continue;
-            if (inSource->moments->x > AnalysisRegion.x1) continue;
-            if (inSource->moments->y > AnalysisRegion.y1) continue;
-        }
+        if (final) {
+            if (inSource->mode &  PM_SOURCE_SUBTRACTED) continue;
+        } else {
+            if (inSource->mode &  PM_SOURCE_BLEND) continue;
+        }
+
+        if (inSource->moments->x < AnalysisRegion.x0) continue;
+        if (inSource->moments->y < AnalysisRegion.y0) continue;
+        if (inSource->moments->x > AnalysisRegion.x1) continue;
+        if (inSource->moments->y > AnalysisRegion.y1) continue;
 
         pmSource *otSource = pmSourceAlloc ();
@@ -73 +74 @@
             // peak-up on peak (for non-sat objects)
 
-            int ix = inSource->peak->x;
-            int iy = inSource->peak->y;
+            // ix,iy must land on inSource->pixels 
+            int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
+            int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));
 
             psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
@@ -104 +106 @@
         otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
         otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
+        otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);
 
         // build the model image 
@@ -117 +120 @@
         psImageKeepCircle (mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
         pmSourceAddModel (flux, mask, model, false, false);
+
+        // calculate nDOF (nPix - 1)
+        // int Nmaskpix = psImageCountPixelMask (mask, allArray, PSPHOT_MASK_SATURATED);
+        // model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;
 
         // save source in list
@@ -128 +135 @@
     // fill out the sparse matrix
     psSparse *sparse = psSparseAlloc (models->n, 100);
+    psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
     for (int i = 0; i < models->n; i++) {
         int N = index->data.U32[i];
@@ -133 +141 @@
         pmSource *Mi = models->data[i];
 
-        // diagonal element (auto-cross-product)
-        f = pmSourceCrossProduct (Mi, Mi);
-        psSparseMatrixElement (sparse, i, i, f);
+        // scale by diagonal element (auto-cross-product)
+        r = pmSourceCrossProduct (Mi, Mi);
+        weight->data.F32[i] = r;
+
+        psSparseMatrixElement (sparse, i, i, 1.0);
 
         // find the image x model value
         f = pmSourceCrossProduct (Fi, Mi);
-        psSparseVectorElement (sparse, i, f);
+        psSparseVectorElement (sparse, i, f / r);
 
         // loop over all other stars following this one
@@ -153 +163 @@
             // got an overlap; calculate cross-product and add to output array
             f = pmSourceCrossProduct (Mi, Mj);
-            psSparseMatrixElement (sparse, j, i, f); 
+            psSparseMatrixElement (sparse, j, i, f / r); 
         }
     }
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
+
+    psSparseConstraint constraint;
+    constraint.paramMin   = 0;
+    constraint.paramMax   = 1e8;
+    constraint.paramDelta = 1e8;
 
     // solve for normalization terms (need include local sky?)
     psSparseResort (sparse);
-    psVector *norm = psSparseSolve (NULL, sparse, 3);
+    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
 
     // adjust models, set sources and subtract
@@ -168 +183 @@
         pmSource *Mi = models->data[i];
 
+        // if we already have a PSF model, free it.
+        psFree (Fi->modelPSF);
+
         // need to increment counter so we can free models here and sources above
         Fi->modelPSF = psMemCopy (Mi->modelPSF);
@@ -173 +191 @@
         // assign linearly-fitted normalization
         Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
+        Fi->modelPSF->dparams->data.F32[1] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
+        // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
 
         // subtract object
         pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
         Fi->mode |= PM_SOURCE_SUBTRACTED;
-        Fi->mode |= PM_SOURCE_TEMPSUB;
-    }
+        if (!final) Fi->mode |= PM_SOURCE_TEMPSUB;
+    }
+
+    // measure chisq for each source
+    for (int i = 0; final && (i < models->n); i++) {
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
+        pmModel *model = Fi->modelPSF;
+
+        x = model->params->data.F32[2];
+        y = model->params->data.F32[3];
+
+        psImageKeepCircle (Fi->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
+        pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
+        psImageKeepCircle (Fi->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
+    }
+
     psFree (index);
     psFree (sparse);
     psFree (models);
     psFree (norm);
+    psFree (weight);
 
     psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
     return true;
 }
+
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight) {
+
+    double dC = 0.0;
+    int Npix = 0;
+    for (int j = 0; j < image->numRows; j++) {
+        for (int i = 0; i < image->numCols; i++) {
+            if (mask->data.U8[j][i]) continue;
+            if (weight->data.F32[j][i] <= 0) continue;
+            dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
+            Npix ++;
+        }
+    }
+    model->nDOF = Npix - 1;
+    model->chisq = dC;
+
+    return (true);
+}