Changeset 7256 for trunk/psModules/src/imcombine/pmReadoutCombine.c

Timestamp:

May 31, 2006, 4:26:32 PM (20 years ago)

Author:

Paul Price

Message:

pmReadoutCombine is working, but it's kind of slow.

File:

• trunk/psModules/src/imcombine/pmReadoutCombine.c (modified) (11 diffs)

trunk/psModules/src/imcombine/pmReadoutCombine.c

@@ -5 +5 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-05-31 22:51:15 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-01 02:26:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -25 +25 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-#if 1
 // Return the statistic of interest
-static double getStat(const psStats *stats, // Statistics structure
-                      psStatsOptions option // Statistics option
-                     )
+static inline double getStat(const psStats *stats, // Statistics structure
+                             psStatsOptions option // Statistics option
+                            )
 {
     switch (option) {
@@ -59 +58 @@
     return NAN;
 }
-#endif
-
-// Mask for combination --- used only locally
-typedef enum {
-    PM_READOUT_COMBINE_CLEAR      = 0x00, // No problem
-    PM_READOUT_COMBINE_NO_IMAGE   = 0x01, // No image available
-    PM_READOUT_COMBINE_MASKED     = 0x02, // Pixel is masked
-    PM_READOUT_COMBINE_BAD        = 0x03, // Pixel is bad (NO_IMAGE | MASKED)
-    PM_READOUT_COMBINE_CLIPPED    = 0x04  // Pixel has been clipped
-} pmReadoutCombineMask;
-
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -96 +84 @@
  *****************************************************************************/
 bool pmReadoutCombine(pmReadout *output,
-                      const psArray *inputs,
+                      psArray *inputs,
                       const psVector *zero,
                       const psVector *scale,
@@ -144 +132 @@
     psVector *mask   = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for stack
     mask->n = inputs->n;
-    psVectorInit(mask, 0);
+
     bool haveWeights = false;           // Do we have weight images?
     bool valid = false;                 // Do we have a single valid input?
@@ -150 +138 @@
         pmReadout *readout = inputs->data[i]; // Readout of interest
         if (!readout || !readout->image) {
-            mask->data.U8[i] = PM_READOUT_COMBINE_NO_IMAGE;
-            continue;
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Input readout %d is NULL or has NULL image.\n", i);
+            return false;
         }
 
@@ -179 +167 @@
         psError(PS_ERR_IO, true, "No valid inputs.\n");
         return NULL;
+    }
+
+    // Apply scale and zero
+    if (scale || zero) {
+        for (int i = 0; i < inputs->n; i++) {
+            pmReadout *readout = inputs->data[i]; // The particular readout
+            if (zero) {
+                psBinaryOp(readout->image, readout->image, "-",
+                           psScalarAlloc(zero->data.F32[i], PS_TYPE_F32));
+            }
+            if (scale) {
+                float scaleFactor = 1.0 / scale->data.F32[i]; // Scaling
+                psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(scaleFactor, PS_TYPE_F32));
+                if (haveWeights) {
+                    psBinaryOp(readout->weight, readout->weight, "*",
+                               psScalarAlloc(scaleFactor * scaleFactor, PS_TYPE_F32));
+                }
+            }
+        }
     }
 
@@ -277 +284 @@
         }
         for (int j = minInputCols; j < maxInputCols; j++) {
-
             int numValid = 0;           // Number of valid pixels in the stack
+            memset(mask->data.U8, 0, mask->n * sizeof(psU8)); // Reset the mask
             for (int r = 0; r < inputs->n; r++) {
                 pmReadout *readout = inputs->data[r]; // Input readout
@@ -284 +291 @@
                 int xIn = j - readout->col0; // x position on input readout
 
-                // Check existence and bounds
-                if (mask->data.U8[r] & PM_READOUT_COMBINE_NO_IMAGE ||
-                        xIn < 0 || xIn >= readout->image->numCols ||
-                        yIn < 0 || yIn >= readout->image->numRows) {
+                psImage *image = readout->image; // The readout image
+
+                // Check bounds
+                if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) {
                     continue;
                 }
+
+                pixels->data.F32[r] = image->data.F32[yIn][xIn];
 
                 // Check mask
                 if (readout->mask && readout->mask->data.U8[yIn][xIn] & maskVal) {
-                    mask->data.U8[r] &= PM_READOUT_COMBINE_MASKED;
+                    mask->data.U8[r] = 1;
                     continue;
                 }
 
-                pixels->data.F32[r] = readout->image->data.F32[yIn][xIn];
-
-                // Apply zero and scale
-                if (zero) {
-                    pixels->data.F32[r] -= zero->data.F32[r];
-                }
-                if (scale) {
-                    pixels->data.F32[r] /= scale->data.F32[r];
-                }
-
-                if (haveWeights) {
-                    weights->data.F32[r] = readout->weight->data.F32[yIn][xIn];
-                    if (scale) {
-                        weights->data.F32[r] /= scale->data.F32[r] * scale->data.F32[r];
-                    }
-                }
                 numValid++;
             }
@@ -331 +324 @@
                 int numHigh = numValid * params->fracHigh; // Number of high pixels to clip
                 // Low pixels
-                for (int k = 0, numMasked = 0; numMasked < numLow; k++) {
+                for (int k = 0, numMasked = 0; numMasked < numLow && k < index->n; k++) {
                     // Don't count the ones that are already masked
-                    if (! mask->data.U8[index->data.S32[k]] & PM_READOUT_COMBINE_BAD) {
-                        mask->data.U8[index->data.S32[k]] |= PM_READOUT_COMBINE_CLIPPED;
+                    if (!mask->data.U8[index->data.S32[k]]) {
+                        mask->data.U8[index->data.S32[k]] = 1;
                         numMasked++;
                     }
                 }
                 // High pixels
-                for (int k = pixels->n, numMasked = 0; numMasked < numHigh; k++) {
+                for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k < index->n; k++) {
                     // Don't count the ones that are already masked
-                    if (! mask->data.U8[index->data.S32[k]] & PM_READOUT_COMBINE_BAD) {
-                        mask->data.U8[index->data.S32[k]] |= PM_READOUT_COMBINE_CLIPPED;
+                    if (! mask->data.U8[index->data.S32[k]]) {
+                        mask->data.U8[index->data.S32[k]] = 1;
                         numMasked++;
                     }
@@ -349 +342 @@
 
             // Combination
-            psVectorStats(stats, pixels, weights, mask, PM_READOUT_COMBINE_BAD | PM_READOUT_COMBINE_CLIPPED);
+            psVectorStats(stats, pixels, weights, mask, 1);
             output->image->data.F32[yOut][xOut] = getStat(stats, params->combine);
             if (haveWeights) {
-                output->weight->data.F32[yOut][xOut] = PS_SQR(getStat(stats, combineStdev)); // Variance
-            }
-
-            // Clear the clipping
-            psBinaryOp(mask, mask, "&", psScalarAlloc(~PM_READOUT_COMBINE_CLIPPED, PS_TYPE_U8));
-
+                float stdev = getStat(stats, combineStdev);
+                output->weight->data.F32[yOut][xOut] = PS_SQR(stdev); // Variance
+            }
         }
     }