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

Timestamp:

Nov 19, 2006, 7:26:08 PM (20 years ago)

Author:

eugene

Message:

working on global linear fit with sky level

File:

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

trunk/psphot/src/psphotEnsemblePSF.c

@@ -10 +10 @@
 // the analysis is performed wrt the simulated pixel values
 
-static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope);
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope, bool constant_weights);
 
 bool psphotEnsemblePSF (pmReadout *readout, psArray *refSources, psMetadata *recipe, pmPSF *psf, bool final) {
@@ -162 +162 @@
     // set the sky, sky_x, sky_y components of border matrix
     // XXX ignore sky slope terms for now
-    SetBorderT (border, readout, refSources, fitSources, index, false);
+    SetBorderT (border, readout, refSources, fitSources, index, false, CONSTANT_PHOTOMETRIC_WEIGHTS);
 
     psSparseConstraint constraint;
@@ -170 +170 @@
 
     // solve for normalization terms (need include local sky?)
-    # if (0)
+    # if (1)
     psVector *norm = NULL;
     psVector *skyfit = NULL;
     psSparseBorderSolve (&norm, &skyfit, constraint, border, 3);
+    fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
     # else
     psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
@@ -226 +227 @@
 
 // calculate the weight terms for the sky fit component of the matrix
-static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope) {
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope, bool constant_weights) {
 
     // generate the image-wide weight terms
@@ -251 +252 @@
     psU8 **mask = readout->mask->data.U8;
 
-    double w, x, y, x2, xy, y2, xc, yc;
-    w = x = y = x2 = xy = y2 = 0;
+    double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
+    w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
     
     int col0 = readout->image->col0;
@@ -260 +261 @@
         for (int i = 0; i < readout->image->numCols; i++) {
             if (mask[j][i]) continue;
-            w  +=       image[j][i]/weight[j][i];
+            if (constant_weights) {
+                wt = 1.0;
+            } else {
+                wt = weight[j][i];
+            }
+            f = image[j][i];
+            w   += 1/wt;
+            fo  += f/wt;
+
             if (!fitSlope) continue;
             xc = i + col0;
             yc = j + row0;
-            x  +=    xc*image[j][i]/weight[j][i];
-            y  +=    yc*image[j][i]/weight[j][i];
-            x2 += xc*xc*image[j][i]/weight[j][i];
-            xy += xc*yc*image[j][i]/weight[j][i];
-            y2 += yc*yc*image[j][i]/weight[j][i];
+            x  +=    xc/wt;
+            y  +=    yc/wt;
+            x2 += xc*xc/wt;
+            xy += xc*yc/wt;
+            y2 += yc*yc/wt;
+            fx +=  f*xc/wt;
+            fy +=  f*yc/wt;
         }
     }
@@ -276 +287 @@
 
     // set the Border T elements
+    psSparseBorderElementG (border, 0, fo);
     psSparseBorderElementT (border, 0, 0, w);
     if (fitSlope) {
+        psSparseBorderElementG (border, 0, fx);
+        psSparseBorderElementG (border, 0, fy);
         psSparseBorderElementT (border, 1, 0, x);
         psSparseBorderElementT (border, 2, 0, y);