Changeset 18860 for trunk/psModules/src/detrend/pmShutterCorrection.c

Timestamp:

Aug 1, 2008, 1:58:37 PM (18 years ago)

Author:

eugene

Message:

splitting off pmShutterCorrectionGeneratePending for threaded coding

File:

• trunk/psModules/src/detrend/pmShutterCorrection.c (modified) (5 diffs)

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -957 +957 @@
 }
 
+bool pmShutterCorrectionGeneratePrepare(pmReadout *shutter, pmReadout *pattern, const psArray *inputs, psMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(shutter, false);
+    PS_ASSERT_PTR_NON_NULL(pattern, false);
+    PS_ASSERT_ARRAY_NON_NULL(inputs, false);
+
+    // determine the output image size based on the input images
+    int row0, col0, numCols, numRows;
+    if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, inputs)) {
+        psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
+        return false;
+    }
+
+    // generate the required output image based on the specified sizes
+    pmReadoutStackDefineOutput(shutter, col0, row0, numCols, numRows, false, false, maskVal);
+    if (pattern) {
+        pmReadoutStackDefineOutput(pattern, col0, row0, numCols, numRows, false, false, maskVal);
+    }
+
+    psImage *nums = pmReadoutSetAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_COUNT, numCols, numRows, PS_TYPE_U16, 0); // Image with number fitted per pixel
+    if (!nums) {
+        return false;
+    }
+    psImage *sigma = pmReadoutSetAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_SIGMA, numCols, numRows, PS_TYPE_F32, NAN); // Image with stdev per pixel
+    if (!sigma) {
+        return false;
+    }
+
+    // Update the "concepts"
+    psList *inputCells = psListAlloc(NULL); // List of cells
+    for (long i = 0; i < inputs->n; i++) {
+        pmReadout *readout = inputs->data[i]; // Readout of interest
+        psListAdd(inputCells, PS_LIST_TAIL, readout->parent);
+    }
+    bool success = pmConceptsAverageCells(shutter->parent, inputCells, NULL, NULL, true);
+    psFree(inputCells);
+
+    // Correct the exposure times --- they don't make sense any more.
+    psMetadataItem *item = psMetadataLookup(shutter->parent->concepts, "CELL.EXPOSURE");
+    item->data.F32 = NAN;
+    item = psMetadataLookup(shutter->parent->concepts, "CELL.DARKTIME");
+    item->data.F32 = NAN;
+
+    shutter->data_exists = true;
+    shutter->parent->data_exists = true;
+    shutter->parent->parent->data_exists = true;
+
+    pattern->data_exists = true;
+    if (pattern->parent) { 
+        pattern->parent->data_exists = true;
+        if (pattern->parent->parent) {
+            pattern->parent->parent->data_exists = true;
+        }
+    }
+
+    return success;
+}
+
 bool pmShutterCorrectionGenerate(pmReadout *shutter, pmReadout *pattern, const psArray *inputs,
                                  float reference, const pmShutterCorrectionData *data,
@@ -962 +1020 @@
 {
     PS_ASSERT_PTR_NON_NULL(shutter, false);
+    PS_ASSERT_PTR_NON_NULL(pattern, false);
     PS_ASSERT_ARRAY_NON_NULL(inputs, false);
     PS_ASSERT_INT_EQUAL(data->num, inputs->n, false);
@@ -975 +1034 @@
     }
 
-    pmReadoutUpdateSize(shutter, minInputCols, minInputRows, xSize, ySize, false, false, maskVal);
-    if (pattern) {
-        pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false, false, maskVal);
-    }
-
-    psImage *nums = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize,
-                                           PS_TYPE_U16, 0); // Image with number fitted per pixel
+    psImage *nums = pmReadoutGetAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_COUNT);
     if (!nums) {
         return false;
     }
-    psImage *sigma = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize,
-                                            PS_TYPE_F32, NAN); // Image with stdev per pixel
+    psImage *sigma = pmReadoutGetAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_SIGMA);
     if (!sigma) {
-        psFree(nums);
         return false;
     }
 
     psImage *shutterImage = shutter->image; // Shutter correction image
-    psImage *patternImage; // Illumination pattern
-    if (pattern) {
-        patternImage = psMemIncrRefCounter(pattern->image);
-    } else {
-        patternImage = psImageAlloc(xSize, ySize, PS_TYPE_F32);
-    }
+    psImage *patternImage = pattern->image; // Illumination pattern
 
     int num = data->num;                // Number of images
@@ -1024 +1070 @@
                     errors->data.F32[r] = sqrtf(readout->weight->data.F32[yIn][xIn]) * ref;
                 } else {
-                    errors->data.F32[r] = sqrtf(image->data.F32[yIn][xIn]) * ref;
+                    // XXX guess that the input data is Poisson distributed; if we go negative, force high
+                    errors->data.F32[r] = sqrtf(fabs(image->data.F32[yIn][xIn])) * ref;
                 }
             }
@@ -1049 +1096 @@
     psFree(errors);
     psFree(counts);
-    psFree(nums);
-    psFree(sigma);
-
-    // Update the "concepts"
-    psList *inputCells = psListAlloc(NULL); // List of cells
-    for (long i = 0; i < inputs->n; i++) {
-        pmReadout *readout = inputs->data[i]; // Readout of interest
-        psListAdd(inputCells, PS_LIST_TAIL, readout->parent);
-    }
-    bool success = pmConceptsAverageCells(shutter->parent, inputCells, NULL, NULL, true);
-    psFree(inputCells);
-
-    // Correct the exposure times --- they don't make sense any more.
-    psMetadataItem *item = psMetadataLookup(shutter->parent->concepts, "CELL.EXPOSURE");
-    item->data.F32 = NAN;
-    item = psMetadataLookup(shutter->parent->concepts, "CELL.DARKTIME");
-    item->data.F32 = NAN;
-
-    shutter->data_exists = true;
-    shutter->parent->data_exists = true;
-    shutter->parent->parent->data_exists = true;
-
-    if (pattern) {
-        pattern->data_exists = true;
-        pattern->parent->data_exists = true;
-        pattern->parent->parent->data_exists = true;
-    }
-
-    return success;
-}
+
+    return true;
+}