Changeset 34800 for trunk/psModules/src/detrend/pmDark.c

Timestamp:

Dec 11, 2012, 2:04:31 PM (13 years ago)

Author:

watersc1

Message:

Merge from my branch including background restoration and stack stage projection cell binned images.

Location:

trunk

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/detrend/pmDark.c (modified) (3 diffs)

trunk/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;