Changeset 17155 for branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c

Timestamp:

Mar 27, 2008, 1:28:38 PM (18 years ago)

Author:

Paul Price

Message:

Adding supplementary output: count and stdev images.

File:

• branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c (modified) (10 diffs)

branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c

@@ -18 +18 @@
 
 //#define SHOW_BUSY 1                   // Show that the function is busy
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -94 +95 @@
     }
 
-    psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination
+    psStatsOptions combineStdev = 0; // Statistics option for weights
+    switch (params->combine) {
+      case PS_STAT_SAMPLE_MEAN:
+      case PS_STAT_SAMPLE_MEDIAN:
+        combineStdev = PS_STAT_SAMPLE_STDEV;
+        break;
+      case PS_STAT_ROBUST_MEDIAN:
+        combineStdev = PS_STAT_ROBUST_STDEV;
+        break;
+      case PS_STAT_FITTED_MEAN:
+        combineStdev = PS_STAT_FITTED_STDEV;
+        break;
+      case PS_STAT_CLIPPED_MEAN:
+        combineStdev = PS_STAT_CLIPPED_STDEV;
+        break;
+      default:
+        psAbort("Should never get here --- checked params->combine before.\n");
+    }
+
+    psStats *stats = psStatsAlloc(params->combine | combineStdev); // The statistics to use in the combination
     if (params->combine == PS_STAT_CLIPPED_MEAN) {
         stats->clipSigma = params->rej;
@@ -107 +127 @@
         }
     }
+    if (params->weights && first) {
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
+                         "Using input weights to combine images", "");
+    }
 
     int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
@@ -120 +144 @@
     psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
 
-    psStatsOptions combineStdev = 0; // Statistics option for weights
-    if (params->weights) {
-        if (first) {
-            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
-                             "Using input weights to combine images", "");
-        }
-
-        // Get the correct statistics option for weights
-        switch (params->combine) {
-        case PS_STAT_SAMPLE_MEAN:
-        case PS_STAT_SAMPLE_MEDIAN:
-            combineStdev = PS_STAT_SAMPLE_STDEV;
-            break;
-        case PS_STAT_ROBUST_MEDIAN:
-            combineStdev = PS_STAT_ROBUST_STDEV;
-            break;
-        case PS_STAT_FITTED_MEAN:
-            combineStdev = PS_STAT_FITTED_STDEV;
-            break;
-        case PS_STAT_CLIPPED_MEAN:
-            combineStdev = PS_STAT_CLIPPED_STDEV;
-            break;
-        default:
-            psAbort("Should never get here --- checked params->combine before.\n");
-        }
-        stats->options |= combineStdev;
-    }
+    psImage *counts = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize,
+                                             PS_TYPE_U16, 0);
+    if (!counts) {
+        return false;
+    }
+    psImage *sigma = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize,
+                                            PS_TYPE_F32, NAN);
+    if (!sigma) {
+        psFree(counts);
+        return false;
+    }
+
+    stats->options |= combineStdev;
 
     // We loop through each pixel in the output image.  We loop through each input readout.  We determine if
@@ -273 +283 @@
                 outputMask[yOut][xOut] = params->blank;
                 outputImage[yOut][xOut] = NAN;
+                counts->data.U16[yOut][xOut] = 0;
+                sigma->data.F32[yOut][xOut] = NAN;
                 continue;
             }
@@ -288 +300 @@
                         maskData[indexData[k]] = 1;
                         numMasked++;
+                        numValid--;
                     }
                 }
@@ -293 +306 @@
                 for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) {
                     // Don't count the ones that are already masked
-                    if (! maskData[indexData[k]]) {
+                    if (!maskData[indexData[k]]) {
                         maskData[indexData[k]] = 1;
                         numMasked++;
+                        numValid--;
                     }
                 }
             }
+            counts->data.U16[yOut][xOut] = numValid;
 
             // XXXXX this step probably is very expensive : convert errors to variance everywhere?
@@ -314 +329 @@
                     outputWeight[yOut][xOut] = NAN;
                 }
+                sigma->data.F32[yOut][xOut] = NAN;
             } else {
                 outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
@@ -323 +339 @@
                     // also, the weighted mean is not obviously the correct thing here
                 }
+                sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev);
             }
         }
@@ -338 +355 @@
     psFree(stats);
     psFree(invScale);
+    psFree(counts);
+    psFree(sigma);
 
     // Update the "concepts"