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

Timestamp:

Nov 16, 2006, 4:45:46 PM (20 years ago)

Author:

eugene

Message:

moved pmSourceChisq to psModules, adding sky correction

File:

• trunk/psphot/src/psphotEnsemblePSF.c (modified) (9 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, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
@@ -29 +27 @@
     psArray  *models = psArrayAllocEmpty (sources->n);
     psVector *index  = psVectorAllocEmpty (sources->n, PS_TYPE_U32);
-    // DROP models->n = index->n = 0;
 
     // option to limit analysis to a specific region
@@ -66 +63 @@
         // really saturated stars should be re-measured for a better centroid
         // XXX EAM : move this to a 'clear satstar function'
+        // XXX : place this in pmSourceMoments or wherever it is called?
         if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
             status = pmSourceMoments (inSource, INNER_RADIUS);
@@ -109 +107 @@
 
         // make temporary copies of the image pixels and mask
+        // we need to have a copy which will not be modified by changes to its neighbor
         otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
         otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
@@ -142 +141 @@
     psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
     psVector *errors = psVectorAlloc (models->n, PS_TYPE_F32);
+
+    // create the border matrix (includes the sparse matrix)
+    // for just sky: 1 row; for x,y terms: 3 rows
+    psSparseBorder *border = psSparseBorderAlloc (sparse, 1);
 
     for (int i = 0; i < models->n; i++) {
@@ -166 +169 @@
         psSparseVectorElement (sparse, i, f / r);
 
+        // add the per-source weights (border region)
+        p = pmSourceWeight (Mi, CONST, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        // px = pmSourceWeight (Mi, XTERM, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        // py = pmSourceWeight (Mi, YTERM, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        psSparseBorderElementB (border, i, 0, p);
+        // psSparseBorderElementB (border, i, 1, px);
+        // psSparseBorderElementB (border, i, 2, py);
+
         // loop over all other stars following this one
         for (int j = i + 1; j < models->n; j++) {
@@ -183 +194 @@
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
+    // generate the image-wide weight terms
+    // first, mask all object pixels
+    for (int i = 0; i < models->n; i++) {
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
+        pmSource *Mi = models->data[i];
+
+        // XXX need to select the correct mode (check this against the selections below
+        // - which model is appropriate for the given source?
+        // - which radius is appropriate for the given source?
+        // - where do x,y come from?
+        psImageKeepCircle (Mi->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
+    }   
+
+    // accumulate the image statistics from the masked regions
+    for (int j = 0; j < numRows; j++) {
+        for (int i = 0; i < numCols; i++) {
+            w  +=     image->data.F32[j][i]/weight->data.F32[j][i];
+            x  +=   i*image->data.F32[j][i]/weight->data.F32[j][i];
+            y  +=   j*image->data.F32[j][i]/weight->data.F32[j][i];
+            x2 += i*i*image->data.F32[j][i]/weight->data.F32[j][i];
+            xy += i*j*image->data.F32[j][i]/weight->data.F32[j][i];
+            y2 += j*j*image->data.F32[j][i]/weight->data.F32[j][i];
+        }
+    }
+    psSparseBorderElementT (border, 0, 0, w);
+    if (0) {
+        psSparseBorderElementT (border, 1, 0, x);
+        psSparseBorderElementT (border, 2, 0, y);
+        psSparseBorderElementT (border, 0, 1, x);
+        psSparseBorderElementT (border, 1, 1, x2);
+        psSparseBorderElementT (border, 2, 1, xy);
+        psSparseBorderElementT (border, 0, 2, y);
+        psSparseBorderElementT (border, 1, 2, xy);
+        psSparseBorderElementT (border, 2, 2, y2);
+    }    
+
     psSparseConstraint constraint;
     constraint.paramMin   = 0.0;
@@ -190 +238 @@
     // solve for normalization terms (need include local sky?)
     psSparseResort (sparse);
-    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
+    psSparseBorderSolve (&xFit, &yFit, constraint, border, 3);
+
+    // psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
 
     // adjust models, set sources and subtract
@@ -244 +294 @@
     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);
-}