Index: trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 6872) +++ trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 6873) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-04 01:01:33 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -36,15 +36,287 @@ } +static void pmCombineParamsFree (pmCombineParams *params) +{ + + if (params == NULL) + return; + + psFree (params->stats); + return; +} + +pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions) +{ + + pmCombineParams *params = psAlloc (sizeof(pmCombineParams)); + psMemSetDeallocator(params, (psFreeFunc) pmCombineParamsFree); + + params->stats = psStatsAlloc (statsOptions); + params->maskVal = 0; + params->fracHigh = 0.25; + params->fracHigh = 0.25; + params->nKeep = 3; + + return (params); +} + /****************************************************************************** XXX: Must add support for S16 and S32 types. F32 currently supported. *****************************************************************************/ psImage *pmReadoutCombine(psImage *output, - const psList *inputs, - psCombineParams *params, + const psArray *inputs, const psVector *zero, const psVector *scale, + pmCombineParams *params, bool applyZeroScale, psF32 gain, psF32 readnoise) +{ + PS_ASSERT_PTR_NON_NULL(inputs, NULL); + PS_ASSERT_PTR_NON_NULL(params, NULL); + PS_ASSERT_PTR_NON_NULL(params->stats, NULL); + if (zero != NULL) { + PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL); + // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(zero, NULL); + } + if (scale != NULL) { + PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL); + // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL); + } + if ((zero != NULL) && (scale != NULL)) { + PS_ASSERT_VECTOR_TYPE_EQUAL(zero, scale, NULL); + // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL); + } + + psStats *stats = params->stats; + psS32 maxInputCols = 0; + psS32 maxInputRows = 0; + psS32 minInputCols = PS_MAX_S32; + psS32 minInputRows = PS_MAX_S32; + pmReadout *tmpReadout = NULL; + psS32 tmpI; + psElemType outputType = PS_TYPE_F32; + + if (DetermineNumBits(stats->options) != 1) { + psError(PS_ERR_UNKNOWN, true, + "Multiple statistical options have been requested. Returning NULL.\n"); + return(NULL); + } + + // We step through each readout in the input image list. If any readout is + // NULL, empty, or has the wrong type, we generate an error and return + // NULL. We determine how big of an output image is needed to combine + // these input images. We do this by taking the + // max(readout->col0 + readout->numCols + image->col0 + image->numCols) + // max(readout->row0 + readout->numRows + image->row0 + image->numRows) + // + for (int i = 0; i < inputs->n; i++) { + tmpReadout = inputs->data[i]; + PS_ASSERT_READOUT_NON_NULL(tmpReadout, output); + PS_ASSERT_READOUT_NON_EMPTY(tmpReadout, output); + PS_ASSERT_READOUT_TYPE(tmpReadout, PS_TYPE_F32, output); + + minInputRows = PS_MIN(minInputRows, (tmpReadout->row0 + tmpReadout->image->row0)); + tmpI = tmpReadout->row0 + + tmpReadout->image->row0 + + tmpReadout->image->numRows; + maxInputRows = PS_MAX(maxInputRows, tmpI); + + minInputCols = PS_MIN(minInputCols, (tmpReadout->col0 + tmpReadout->image->col0)); + tmpI = tmpReadout->col0 + + tmpReadout->image->col0 + + tmpReadout->image->numCols; + maxInputCols = PS_MAX(maxInputCols, tmpI); + } + + // We ensure that the zero vector is of the proper size. + if (zero != NULL) { + PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL); + if (zero->n < inputs->n) { + psError(PS_ERR_UNKNOWN, true, "zero vector has incorrect size (%d). Returning NULL.\n", zero->n); + return(NULL); + } else if (zero->n > inputs->n) { + // XXX EAM : abort on this condition? is probably an error + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: the zero vector too many elements (%d)\n", zero->n); + } + } + + // We ensure that the scale vector is of the proper size. + if (scale != NULL) { + PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL); + if (scale->n < inputs->n) { + psError(PS_ERR_UNKNOWN, true, "scale vector has incorrect size (%d). Returning NULL.\n", scale->n); + return(NULL); + } else if (scale->n > inputs->n) { + // XXX EAM : abort on this condition? is probably an error + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: the scale vector has too many elements (%d)\n", scale->n); + } + } + + // At this point, the following variables have been computed: + // maxInputRows: the largest input row value, in output image space. + // maxInputCols: the largest input column value, in output image space. + // minInputRows: the smallest input row value, in output image space. + // minInputCols: the smallest input column value, in output image space. + // + if (output == NULL) { + output = psImageAlloc(maxInputCols-minInputCols, maxInputRows-minInputRows, outputType); + *(psS32 *) &(output->col0) = minInputCols; + *(psS32 *) &(output->row0) = minInputRows; + } else { + + // XXX EAM : recycle the existing output image? why would we care about the existing pixels? + PS_ASSERT_IMAGE_TYPE(output, PS_TYPE_F32, NULL); + if (((output->col0 + output->numCols) < maxInputCols) || + ((output->row0 + output->numRows) < maxInputRows)) { + psError(PS_ERR_UNKNOWN, true, + "Output image (%d, %d) is too small to hold combined images. Returning NULL.\n", + output->row0 + output->numRows, + output->col0 + output->numCols); + return(NULL); + } + + // reset output origin using logic of above + *(psS32 *) &(output->col0) = minInputCols; + *(psS32 *) &(output->row0) = minInputRows; + } + + psVector *tmpPixels = psVectorAlloc(inputs->n, PS_TYPE_F32); + psVector *tmpPixelsKeep = psVectorAlloc(inputs->n, PS_TYPE_F32); + psVector *outRowLower = psVectorAlloc(inputs->n, PS_TYPE_U32); + psVector *outRowUpper = psVectorAlloc(inputs->n, PS_TYPE_U32); + psVector *outColLower = psVectorAlloc(inputs->n, PS_TYPE_U32); + psVector *outColUpper = psVectorAlloc(inputs->n, PS_TYPE_U32); + + // For each input readout, we store the min/max pixel indices for that readout, in detector coordinates, + // in the psVectors (outRowLower, outColLower, outRowUpper, outColUpper). + for (int i = 0; i < inputs->n; i++) { + tmpReadout = (pmReadout *) inputs->data[i]; + outRowLower->data.U32[i] = tmpReadout->row0 + tmpReadout->image->row0; + outColLower->data.U32[i] = tmpReadout->col0 + tmpReadout->image->col0; + outRowUpper->data.U32[i] = tmpReadout->row0 + + tmpReadout->image->row0 + + tmpReadout->image->numRows; + outColUpper->data.U32[i] = tmpReadout->col0 + + tmpReadout->image->col0 + + tmpReadout->image->numCols; + } + + // We loop through each pixel in the output image. We loop through each + // input readout. We determine if that output pixel is contained in the + // image from that readout. If so, we save it in psVector tmpPixels. + // If not, we set a mask for that element in tmpPixels. Then, we mask off + // pixels not between fracLow and fracHigh. Then we call the vector + // stats routine on those pixels/mask. Then we set the output pixel value + // to the result of the stats call. + + int nx, ny; + int nKeep, nMin; + float keepFrac = 1.0 - params->fracLow - params->fracHigh; + float value = 0; + psF32 *saveVector = tmpPixelsKeep->data.F32; + + for (int j = output->row0; j < (output->row0 + output->numRows) ; j++) { + if (j % 10 == 0) + fprintf (stderr, "."); + for (int i = output->col0; i < (output->col0 + output->numCols) ; i++) { + int nPix = 0; + for (int r = 0; r < inputs->n; r++) { + tmpReadout = (pmReadout *) inputs->data[r]; + + // psTrace (__func__, 6, "[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outColLower->data.U32[r], outRowLower->data.U32[r], outColUpper->data.U32[r], outRowUpper->data.U32[r]); + if (i < outColLower->data.U32[r]) + continue; + if (i >= outColUpper->data.U32[r]) + continue; + if (j < outRowLower->data.U32[r]) + continue; + if (j >= outRowUpper->data.U32[r]) + continue; + + nx = i - (tmpReadout->col0 + tmpReadout->image->col0); + ny = j - (tmpReadout->row0 + tmpReadout->image->row0); + + if (tmpReadout->mask != NULL) { + if (tmpReadout->mask->data.U8[ny][nx] && params->maskVal) + continue; + } + + tmpPixels->data.F32[nPix] = tmpReadout->image->data.F32[ny][nx]; + // psTrace (__func__, 6, "readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[nPix]); + nPix ++; + } + tmpPixels->n = nPix; + + // are there enough valid pixels to apply fracLow,fracHigh? + nKeep = nPix * keepFrac; + if (nKeep >= params->nKeep) { + psVectorSort (tmpPixels, tmpPixels); + nMin = nPix * params->fracLow; + tmpPixelsKeep->data.F32 = &tmpPixels->data.F32[nMin]; + tmpPixelsKeep->n = nKeep; + } else { + tmpPixelsKeep->data.F32 = tmpPixels->data.F32; + tmpPixelsKeep->n = nPix; + } + + // tmpPixelsKeep is already sorted. sample mean and median are very easy + if (stats->options & PS_STAT_SAMPLE_MEAN) { + value = 0; + for (int r = 0; r < tmpPixelsKeep->n; r++) { + value += tmpPixelsKeep->data.F32[r]; + } + if (tmpPixelsKeep->n == 0) { + value = 0; + } else { + value = value / tmpPixelsKeep->n; + } + } + if (stats->options & PS_STAT_SAMPLE_MEDIAN) { + int r = tmpPixelsKeep->n / 2; + if (tmpPixelsKeep->n == 0) { + value = 0; + goto got_value; + } + if (tmpPixelsKeep->n % 2 == 1) { + int r = 0.5*tmpPixelsKeep->n; + value = tmpPixelsKeep->data.F32[r]; + goto got_value; + } + if (tmpPixelsKeep->n % 2 == 0) { + value = 0.5*(tmpPixelsKeep->data.F32[r] + + tmpPixelsKeep->data.F32[r-1]); + goto got_value; + } + } +got_value: + output->data.F32[j-output->row0][i-output->col0] = value; + } + } + tmpPixelsKeep->data.F32 = saveVector; + + psFree(tmpPixels); + psFree(tmpPixelsKeep); + psFree(outRowLower); + psFree(outRowUpper); + psFree(outColLower); + psFree(outColUpper); + + return(output); +} + +/****************************************************************************** +XXX: Must add support for S16 and S32 types. F32 currently supported. + *****************************************************************************/ +psImage *pmReadoutCombine_OLD(psImage *output, + const psList *inputs, + pmCombineParams *params, + const psVector *zero, + const psVector *scale, + bool applyZeroScale, + psF32 gain, + psF32 readnoise) { PS_ASSERT_PTR_NON_NULL(inputs, NULL); @@ -418,6 +690,6 @@ psRegion minRegion; psRegion maxRegion; - psStats *minStats = psStatsAlloc(PS_STAT_FITTED_MEAN); - psStats *maxStats = psStatsAlloc(PS_STAT_FITTED_MEAN); + psStats *minStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + psStats *maxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); psStats *diffStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); psVector *diffs = psVectorAlloc(fringePoints->n, PS_TYPE_F32); @@ -454,6 +726,6 @@ } - fp->midValue = 0.5 * (maxStats->fittedMean + minStats->fittedMean); - fp->delta = maxStats->fittedMean - minStats->fittedMean; + fp->midValue = 0.5 * (maxStats->robustMedian + minStats->robustMedian); + fp->delta = maxStats->robustMedian - minStats->robustMedian; diffs->data.F32[i] = fp->delta; } @@ -464,11 +736,9 @@ psFree(diffs); if (diffStats == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not determine fitted median of the differences.\n"); + psError(PS_ERR_UNKNOWN, true, "Could not determine robust median of the differences.\n"); return(NULL); } return(diffStats); } - - /** Index: trunk/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 6872) +++ trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 6873) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-17 18:01:05 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii