Changeset 34442 for branches/czw_branch/20120906/psModules/src/detrend

Timestamp:

Sep 14, 2012, 6:10:32 PM (14 years ago)

Author:

watersc1

Message:

Still doesn't work correctly. I've simplified and corrected the problems I was having before, but I still have residuals of ~50 pixels in the warped background model relative to the science warp. I'm beginning to think that there's some issue with the interpolated background model, and that it's not quite the correct size? That's pretty much all I can think of at this point.

File:

• branches/czw_branch/20120906/psModules/src/detrend/pmDark.c (modified) (3 diffs)

branches/czw_branch/20120906/psModules/src/detrend/pmDark.c

@@ -353 +353 @@
 
     // retrieve the required parameter vectors
-    psArray *values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
-    psAssert(values, "values not supplied");
+    psArray *in_values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
+    psAssert(in_values, "values not supplied");
     psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK");
     psAssert(roMask, "roMask not supplied");
-    psVector *orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
-    psAssert(orders, "orders not supplied");
-
-    psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting
+    psVector *max_orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
+    psAssert(max_orders, "orders not supplied");
+
+    psArray *values_set = psArrayAlloc(max_orders->n);
+    psArray *poly_set = psArrayAlloc(max_orders->n);
+    psVector *logL = psVectorAlloc(max_orders->n,PS_TYPE_F64);
+
+    for (int i = 0; i < max_orders->n; i++) {
+      psVector *orders = psVectorAlloc(i+1,PS_TYPE_U8);
+      for (int j = 0; j < orders->n; j++) {
+        orders->data.U8[j] = max_orders->data.U8[j];
+      }
+      poly_set->data[i] =  psPolynomialMDAlloc(orders); // Polynomial for fitting
+      
+      psArray *values = psArrayAlloc(in_values->n);
+      
+      for (int j = 0; j < values->n; j++) {
+        psVector *these_values = psVectorAlloc(i+1,PS_TYPE_F32);
+        psVector *input_values = in_values->data[j];
+
+        for (int k = 0; k < orders->n; k++) {
+          these_values->data.F32[k] = input_values->data.F32[k];
+        }
+        values->data[j] = these_values;
+      }
+      values_set->data[i] = values;
+    }
 
     // retrieve the norm vector, if supplied
@@ -383 +406 @@
     }
 
-    pmDarkVisualInit(values);
+    pmDarkVisualInit(values_set->data[max_orders->n - 1]);
 
     pmReadout *outReadout = output->readouts->data[0];
@@ -423 +446 @@
             }
 
-            if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
+            int k_best = 0;
+            for (int k = 0; k < max_orders->n; k++) {
+              psPolynomialMD *poly = poly_set->data[k];
+              psArray *values = values_set->data[k];
+              
+              if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
                 psErrorClear();         // Nothing we can do about it
                 psVectorInit(poly->coeff, NAN);
-            }
-
-            pmDarkVisualPixelFit(pixels, mask);
-            pmDarkVisualPixelModel(poly, values);
-
-            for (int k = 0; k < poly->coeff->n; k++) {
+              }
+
+              pmDarkVisualPixelFit(pixels, mask);
+              pmDarkVisualPixelModel(poly, values);
+
+              // Insert math here to choose optimum model.
+              logL->data.F64[k] = 0.0;
+              psPolynomialMD *polySig = poly_set->data[0];
+              for (int m = 0; m < poly->deviations->n; m++) {
+                logL->data.F64[k] += pow(poly->deviations->data.F32[m] / polySig->stdevFit,2);
+/*              if ((xOut == 20) && (yOut == 256)) { */
+/*                psTrace("psModules.detrend",3,"pmDarkCombine DEV: %d %d: input %d models: Norders: %d logL: %g value: %g\n", */
+/*                        xOut,yOut,m,k,logL->data.F64[k],poly->deviations->data.F32[m]); */
+/*              } */
+              }
+              if (k > 0) {
+                if ( ( logL->data.F64[k - 1] - logL->data.F64[k] ) > 1) { // Hard coded criterion for a ~5% limit with one degree of freedom
+                  k_best = k;
+                }
+              }
+              if ((xOut <= 600) && (yOut <= 600)) {
+                psTrace("psModules.detrend",3,"pmDarkCombine: %d %d: models: Norders: %d logL: %g BestOrders: %d\n",
+                        xOut,yOut,k,logL->data.F64[k],k_best);
+              }
+            }
+            if (k_best > 1) {
+              k_best = 1;
+            }
+/*          k_best = 1; */
+            // Select the polynomial that seems best.
+            psPolynomialMD *poly = poly_set->data[k_best];
+              
+            //            for (int k = 0; k < poly->coeff->n; k++) {
+            for (int k = 0; k < max_orders->n + 1; k++) { // There is one more coefficient than is stored here.
                 pmReadout *ro = output->readouts->data[k]; // Readout of interest
-                ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+                if (k < poly->coeff->n) {
+                  ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+                }
+                else {
+                  ro->image->data.F32[yOut][xOut] = 0.0;
+                }
             }
             counts->data.U16[yOut][xOut] = poly->numFit;