Index: trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 7254) +++ trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 7256) @@ -5,6 +5,6 @@ * @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,9 +25,8 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -#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,15 +58,4 @@ 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,5 +84,5 @@ *****************************************************************************/ bool pmReadoutCombine(pmReadout *output, - const psArray *inputs, + psArray *inputs, const psVector *zero, const psVector *scale, @@ -144,5 +132,5 @@ 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,6 +138,6 @@ 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,4 +167,23 @@ 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,6 +284,6 @@ } 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,33 +291,19 @@ 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,16 +324,16 @@ 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,13 +342,10 @@ // 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 + } } }