Changeset 10574

Timestamp:

Dec 8, 2006, 10:46:16 AM (20 years ago)

Author:

Paul Price

Message:

Adding additional metadata to output.

File:

• trunk/psModules/src/detrend/pmBias.c (modified) (11 diffs)

trunk/psModules/src/detrend/pmBias.c

@@ -110 +110 @@
 
 // Produce an overscan vector from an array of pixels
-static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options
+static psVector *overscanVector(float *chi2, // chi^2 from fit
+                                pmOverscanOptions *overscanOpts, // Overscan options
                                 const psArray *pixels, // Array of vectors containing the pixel values
                                 psStats *myStats // Statistic to use in reducing the overscan
@@ -145 +146 @@
 
     // Fit the overscan, if required
+    psVector *fitted;                   // Fitted overscan values
     switch (overscanOpts->fitType) {
     case PM_FIT_NONE:
         // No fitting --- that's easy.
+        fitted = psMemIncrRefCounter(reduced);
         break;
     case PM_FIT_POLY_ORD:
@@ -159 +162 @@
         }
         psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
-        psFree(reduced);
-        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
         break;
     case PM_FIT_POLY_CHEBY:
@@ -172 +174 @@
         }
         psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
-        psFree(reduced);
-        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
         break;
     case PM_FIT_SPLINE:
@@ -179 +180 @@
         // input spline
         overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
-        psFree(reduced);
-        reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate);
+        fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate);
         break;
     default:
@@ -188 +188 @@
     }
 
+    if (chi2) {
+        *chi2 = 0.0;                    // chi^2 (sort of)
+        for (int i = 0; i < reduced->n; i++) {
+            *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]);
+        }
+    }
+
+    psFree(reduced);
     psFree(ordinate);
     psFree(mask);
 
-    return reduced;
+    return fitted;
 }
 
@@ -273 +281 @@
             // We do the regular overscan subtraction
             bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction
+            float chi2 = NAN;           // chi^2 from fit
 
             if (readRows) {
@@ -301 +310 @@
 
                 // Reduce the overscans
-                psVector *reduced = overscanVector(overscanOpts, pixels, stats);
+                psVector *reduced = overscanVector(&chi2, overscanOpts, pixels, stats);
                 psFree(pixels);
                 if (! reduced) {
@@ -344 +353 @@
 
                 // Reduce the overscans
-                psVector *reduced = overscanVector(overscanOpts, pixels, stats);
+                psVector *reduced = overscanVector(&chi2, overscanOpts, pixels, stats);
                 psFree(pixels);
                 if (! reduced) {
@@ -364 +373 @@
             case PM_FIT_POLY_CHEBY: {
                     psString comment = NULL;    // Comment to add
-                    psStringAppend(&comment, "Overscan fit: ");
+                    psStringAppend(&comment, "Overscan fit (chi^2: %.2f):", chi2);
                     psPolynomial1D *poly = overscanOpts->poly; // The polynomial
                     for (int i = 0; i < poly->nX; i++) {
@@ -376 +385 @@
                     psSpline1D *spline = overscanOpts->spline; // The spline
                     for (int i = 0; i < spline->n; i++) {
-                        psStringAppend(&comment, "Overscan fit %d: ", i);
+                        psStringAppend(&comment, "Overscan fit (chi^2: %.2f) %d :", chi2, i);
                         psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
                         for (int j = 0; j < poly->nX; j++) {