Changeset 10171

Timestamp:

Nov 22, 2006, 7:28:54 PM (19 years ago)

Author:

eugene

Message:

added code to fit global sky (and slopes), with config options

File:

• trunk/psphot/src/psphotEnsemblePSF.c (modified) (14 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, bool constant_weights);
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool constant_weights, int SKY_FIT_ORDER);
 
 bool psphotEnsemblePSF (pmReadout *readout, psArray *refSources, psMetadata *recipe, pmPSF *psf, bool final) {
@@ -19 +19 @@
     float y;
     float f;
-    float r;
+    // float r;
 
     psTimerStart ("psphot");
@@ -36 +36 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
 
-    const bool CONSTANT_PHOTOMETRIC_WEIGHTS =
+    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     if (!status) {
         psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
+    }
+    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
+    if (!status) {
+        SKY_FIT_ORDER = 0;
+    }
+    bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
+    if (!status) {
+        SKY_FIT_LINEAR = false;
     }
 
@@ -100 +108 @@
 
     // vectors to store stats for each object
-    psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32);
+    // psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32);
     psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32);
 
@@ -106 +114 @@
     // for just sky: 1 row; for x,y terms: 3 rows
     psSparse *sparse = psSparseAlloc (fitSources->n, 100);
-    psSparseBorder *border = psSparseBorderAlloc (sparse, 1);
+    int nBorder = (SKY_FIT_ORDER == 0) ? 1 : 3;
+    psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
 
     // fill out the sparse matrix elements and border elements (B)
@@ -116 +125 @@
         pmSource *Fi = fitSources->data[i];
 
-        // scale by diagonal element (auto-cross-product)
-        r = pmSourceCrossProduct (Fi, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        weight->data.F32[i] = r;
+        // diagonal elements of the sparse matrix (auto-cross-product)
+        f = pmSourceCrossProduct (Fi, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        psSparseMatrixElement (sparse, i, i, f);
+
+        // XXX being used? weight->data.F32[i] = r;
 
         // the formal error depends on the weighting scheme
@@ -125 +136 @@
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
-            errors->data.F32[i] = 1.0 / sqrt(r);
+            errors->data.F32[i] = 1.0 / sqrt(f);
         }
 
-        psSparseMatrixElement (sparse, i, i, 1.0);
 
         // find the image x model value
         f = pmSourceCrossProduct (Ri, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        psSparseVectorElement (sparse, i, f / r);
+        psSparseVectorElement (sparse, i, f);
 
         // add the per-source weights (border region)
-        float p = pmSourceWeight (Fi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        // float px = pmSourceWeight (Fi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        // float py = pmSourceWeight (Fi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        psSparseBorderElementB (border, i, 0, p);
-        // psSparseBorderElementB (border, i, 1, px);
-        // psSparseBorderElementB (border, i, 2, py);
+        switch (SKY_FIT_ORDER) {
+          case 1:
+            f = pmSourceWeight (Fi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            psSparseBorderElementB (border, i, 1, f);
+            f = pmSourceWeight (Fi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            psSparseBorderElementB (border, i, 2, f);
+
+          case 0:
+            f = pmSourceWeight (Fi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            psSparseBorderElementB (border, i, 0, f);
+            break;
+
+          default:
+            psAbort(PS_FILE_LINE, "Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER);
+            break;
+        }
 
         // loop over all other stars following this one
@@ -154 +174 @@
             // got an overlap; calculate cross-product and add to output array
             f = pmSourceCrossProduct (Fi, Fj, CONSTANT_PHOTOMETRIC_WEIGHTS);
-            psSparseMatrixElement (sparse, j, i, f / r);
+            psSparseMatrixElement (sparse, j, i, f);
         }
     }
@@ -161 +181 @@
 
     // 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, CONSTANT_PHOTOMETRIC_WEIGHTS);
+    SetBorderT (border, readout, refSources, fitSources, index, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER);
 
     psSparseConstraint constraint;
@@ -170 +189 @@
 
     // solve for normalization terms (need include local sky?)
-    # 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);
-    psVector *skyfit = NULL;
-    # endif
+    if (SKY_FIT_LINEAR) {
+        psSparseBorderSolve (&norm, &skyfit, constraint, border, 5);
+        fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
+    } else {
+        norm = psSparseSolve (NULL, constraint, sparse, 5);
+        skyfit = NULL;
+    }
 
     // adjust fitSources, set refSources and subtract
@@ -220 +239 @@
     psFree (norm);
     psFree (skyfit);
-    psFree (weight);
     psFree (errors);
     psFree (border);
@@ -229 +247 @@
 
 // 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, bool constant_weights) {
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool constant_weights, int SKY_FIT_ORDER) {
 
     // generate the image-wide weight terms
@@ -271 +289 @@
             w   += 1/wt;
             fo  += f/wt;
-
-            if (!fitSlope) continue;
+            if (SKY_FIT_ORDER == 0) continue;
+
             xc = i + col0;
             yc = j + row0;
@@ -291 +309 @@
     psSparseBorderElementG (border, 0, fo);
     psSparseBorderElementT (border, 0, 0, w);
-    if (fitSlope) {
+    if (SKY_FIT_ORDER > 0) {
         psSparseBorderElementG (border, 0, fx);
         psSparseBorderElementG (border, 0, fy);