Index: trunk/psModules/src/pmSubtractSky.c =================================================================== --- trunk/psModules/src/pmSubtractSky.c (revision 2584) +++ trunk/psModules/src/pmSubtractSky.c (revision 2585) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-01 21:20:24 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-01 21:26:17 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -32,213 +32,198 @@ } -PS_STAT_SAMPLE_MEAN = 0x000001, - PS_STAT_SAMPLE_MEDIAN = 0x000002, - PS_STAT_SAMPLE_STDEV = 0x000004, - PS_STAT_SAMPLE_QUARTILE = 0x000008, - PS_STAT_ROBUST_MEAN = 0x000010, - PS_STAT_ROBUST_MEDIAN = 0x000020, - PS_STAT_ROBUST_MODE = 0x000040, - PS_STAT_ROBUST_STDEV = 0x000080, - PS_STAT_ROBUST_QUARTILE = 0x000100, - PS_STAT_CLIPPED_MEAN = 0x000200, - PS_STAT_CLIPPED_STDEV = 0x000400, - PS_STAT_MAX = 0x000800, - PS_STAT_MIN = 0x001000, - PS_STAT_USE_RANGE = 0x002000, - PS_STAT_USE_BINSIZE = 0x004000, - PS_STAT_ROBUST_FOR_SAMPLE = 0x008000 - - psU64 getHighestPriorityStatOption(psU64 statOptions) - { - switch (statOptions) { - case PS_STAT_SAMPLE_MEAN: - return(PS_STAT_SAMPLE_MEAN); - case PS_STAT_SAMPLE_MEDIAN: - return(PS_STAT_SAMPLE_MEDIAN); - case PS_STAT_CLIPPED_MEAN: - return(PS_STAT_CLIPPED_MEAN); - case PS_STAT_ROBUST_MEAN: - return(PS_STAT_ROBUST_MEAN); - case PS_STAT_ROBUST_MEDIAN: - return(PS_STAT_ROBUST_MEDIAN); - case PS_STAT_ROBUST_MODE: - return(PS_STAT_ROBUST_MODE); - } - psError(XXX); - return(-1); - } - - - /****************************************************************************** - psImage *binImage(origImage, binFactor, statOptions): This routine takes an - input psImage and scales it smaller by a factor of binFactor. The statistic - used in combining input pixels is specified in statOptions. - - XXX: use static vectors for myStats, binVector and binMask. - *****************************************************************************/ - psImage *binImage(psImage *origImage, - int binFactor, - psStatsOptions statOptions) - { - psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32); - psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8); - psStats *myStats = psStatsAlloc(statOptions) - - if (binFactor <= 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: binImage(): binFactor is %d\n", binFactor); - return(origImage); - } - if (binFactor == 1) { - return(origImage); - } - - for (int row = 0; row < origImage->numRows ; row+=binFactor) { - for (int col = 0; col < origImage->numCols ; col+=binFactor) { - int count = 0; - for (int binRow = 0; binRow <= binFactor ; binRow++) { - for (int binCol = 0; binCol <= binFactor ; binCol++) { - if (((row + binRow) < origImage->numRows) && - ((col + binCol) < origImage->numCol)) { - binVector->data.F32[count] = - origImage->data.F32[row + binRow][col + binCol]; - binMask->data.U8[count] = 0; - } else { - binVector->data.F32[count] = 0.0; - binMask->data.U8[count] = 1; - } - count++; - } - } - myStats = psVectorStats(myStats, binVector, binMask, 1); - psF64 statValue; - psBool rc = p_psGetStatValue(myStats, &statValue); - - if (rc == true) { - origImage->data.F32[row][col] = (psF32) statValue; - } else { - origImage->data.F32[row][col] = 0.0; - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(), binImage(): p_psGetStatValue() was FALSE\n"); - } - } - } - psFree(myStats); - psFree(binVector); - psFree(binMask); - - return(origImage); - } - - - /****************************************************************************** - psReadout pmSubtractSky(): - - XXX: use static vectors for myStats. - *****************************************************************************/ - psReadout *pmSubtractSky(psReadout *in, - void *fitSpec, - psFit fit, - int binFactor, - psStats *stats, - float clipSD) - { - if (fitSpec == NULL) { - return(in); - } - if ((fit != PM_FIT_POLYNOMIAL) && (fit != PM_FIT_SPLINE)) { - // No fit is specified. - return(in); - } - - psImage *origImage = in->image; - psImage *binnedImage = NULL; - psPolynomial1D *myPoly; - psSpline1D *mySpline; - - psStatsOptions statOptions = stats->options; - if (1 < p_psDetermineNumBits(statOptions) { - //XXX psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n"); - statOptions = getHighestPriorityStatOption(statOptions) - ; - } - - // Bin the input image according to input parameters. - if ((binFactor <= 0) || - (stats == NULL)) { - if (binFactor <= 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor); - } - if (stats == NULL) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): input parameter stats is NULL\n"); - } - binnedImage = origImage; - } else { - binnedImage = binImage(origImage, binFactor, statOptions); - } - - // Clip pixels that are outside the acceptable range. - if (clipSD <= 0.0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD) - ; - } else { - psF64 binnedMean; - psF64 binnedStdev; - psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - myStats = psImageStats(myStats, binedImage, NULL, 0); - p_psGetStatValue(myStats, &binnedMean); - - myStats->options = PS_STAT_SAMPLE_STDEV; - myStats = psImageStats(myStats, binedImage, NULL, 0); - p_psGetStatValue(myStats, &binnedStdev); - psFree(myStats); - - for (int row = 0; row < binnedImage->numRows ; row++) { - for (int col = 0; col < binnedImage->numCols ; col++) { - if (fabs(binnedImage->data.F32[row][col] - binnedMean) > - (clipSD * binnedStdev)) { - binnedImage->mask.U8[row][col] = 1; - } - } - } - } - - // XXX: fit the polynomial to the binned image - if (fit == PM_FIT_POLYNOMIAL) { - // Fit a polynomial to the old overscan vector. - myPoly = (psPolynomial2D *) - fitSpec; - myPoly = psImageFitPolynomial(myPoly, origImage); - // XXX Do we need to do something with ordinate scaling if Chebyshev? - binnedImage = psImageEvalPolynomial(binnedImage, myPoly); - - } else if (fit == PM_FIT_SPLINE) { - // Fit a spline to the old overscan vector. - mySpline = (psSpline1D *) - fitSpec; - // XXX: What do we do? We don't have 2-D splines. - return(in); - } - - //Subtract the polynomially fitted image from the original image - if (binFactor <= 0) { - for (int row = 0; row < origImage->numRows ; row++) { - for (int col = 0; col < origImage->numCols ; col++) { - origImage->data.F32[row][col]-= binnedImage->data.F32[row][col]; - } - } - } else { - for (int row = 0; row < origImage->numRows ; row++) { - for (int col = 0; col < origImage->numCols ; col++) { - int binRow = row / binFactor; - int binCol = col / binFactor; - origImage->data.F32[row][col]-= binnedImage->data.F32[binRow][binCol]; - } - } - - } - return(in); - } + + +psU64 getHighestPriorityStatOption(psU64 statOptions) +{ + switch (statOptions) { + case PS_STAT_SAMPLE_MEAN: + return(PS_STAT_SAMPLE_MEAN); + case PS_STAT_SAMPLE_MEDIAN: + return(PS_STAT_SAMPLE_MEDIAN); + case PS_STAT_CLIPPED_MEAN: + return(PS_STAT_CLIPPED_MEAN); + case PS_STAT_ROBUST_MEAN: + return(PS_STAT_ROBUST_MEAN); + case PS_STAT_ROBUST_MEDIAN: + return(PS_STAT_ROBUST_MEDIAN); + case PS_STAT_ROBUST_MODE: + return(PS_STAT_ROBUST_MODE); + } + psError(XXX); + return(-1); +} + + +/****************************************************************************** +psImage *binImage(origImage, binFactor, statOptions): This routine takes an +input psImage and scales it smaller by a factor of binFactor. The statistic +used in combining input pixels is specified in statOptions. + +XXX: use static vectors for myStats, binVector and binMask. + *****************************************************************************/ +psImage *binImage(psImage *origImage, + int binFactor, + psStatsOptions statOptions) +{ + psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32); + psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8); + psStats *myStats = psStatsAlloc(statOptions) + + if (binFactor <= 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: binImage(): binFactor is %d\n", binFactor); + return(origImage); + } + if (binFactor == 1) { + return(origImage); + } + + for (int row = 0; row < origImage->numRows ; row+=binFactor) { + for (int col = 0; col < origImage->numCols ; col+=binFactor) { + int count = 0; + for (int binRow = 0; binRow <= binFactor ; binRow++) { + for (int binCol = 0; binCol <= binFactor ; binCol++) { + if (((row + binRow) < origImage->numRows) && + ((col + binCol) < origImage->numCol)) { + binVector->data.F32[count] = + origImage->data.F32[row + binRow][col + binCol]; + binMask->data.U8[count] = 0; + } else { + binVector->data.F32[count] = 0.0; + binMask->data.U8[count] = 1; + } + count++; + } + } + myStats = psVectorStats(myStats, binVector, binMask, 1); + psF64 statValue; + psBool rc = p_psGetStatValue(myStats, &statValue); + + if (rc == true) { + origImage->data.F32[row][col] = (psF32) statValue; + } else { + origImage->data.F32[row][col] = 0.0; + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(), binImage(): p_psGetStatValue() was FALSE\n"); + } + } + } + psFree(myStats); + psFree(binVector); + psFree(binMask); + + return(origImage); +} + + +/****************************************************************************** +psReadout pmSubtractSky(): + +XXX: use static vectors for myStats. + *****************************************************************************/ +psReadout *pmSubtractSky(psReadout *in, + void *fitSpec, + psFit fit, + int binFactor, + psStats *stats, + float clipSD) +{ + if (fitSpec == NULL) { + return(in); + } + if ((fit != PM_FIT_POLYNOMIAL) && (fit != PM_FIT_SPLINE)) { + // No fit is specified. + return(in); + } + + psImage *origImage = in->image; + psImage *binnedImage = NULL; + psPolynomial1D *myPoly; + psSpline1D *mySpline; + + psStatsOptions statOptions = stats->options; + if (1 < p_psDetermineNumBits(statOptions) { + //XXX psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n"); + statOptions = getHighestPriorityStatOption(statOptions) + ; + } + + // Bin the input image according to input parameters. + if ((binFactor <= 0) || + (stats == NULL)) { + if (binFactor <= 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor); + } + if (stats == NULL) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): input parameter stats is NULL\n"); + } + binnedImage = origImage; + } else { + binnedImage = binImage(origImage, binFactor, statOptions); + } + + // Clip pixels that are outside the acceptable range. + if (clipSD <= 0.0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD) + ; + } else { + psF64 binnedMean; + psF64 binnedStdev; + psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + myStats = psImageStats(myStats, binedImage, NULL, 0); + p_psGetStatValue(myStats, &binnedMean); + + myStats->options = PS_STAT_SAMPLE_STDEV; + myStats = psImageStats(myStats, binedImage, NULL, 0); + p_psGetStatValue(myStats, &binnedStdev); + psFree(myStats); + + for (int row = 0; row < binnedImage->numRows ; row++) { + for (int col = 0; col < binnedImage->numCols ; col++) { + if (fabs(binnedImage->data.F32[row][col] - binnedMean) > + (clipSD * binnedStdev)) { + binnedImage->mask.U8[row][col] = 1; + } + } + } + } + + // XXX: fit the polynomial to the binned image + if (fit == PM_FIT_POLYNOMIAL) { + // Fit a polynomial to the old overscan vector. + myPoly = (psPolynomial2D *) + fitSpec; + myPoly = psImageFitPolynomial(myPoly, origImage); + // XXX Do we need to do something with ordinate scaling if Chebyshev? + binnedImage = psImageEvalPolynomial(binnedImage, myPoly); + + } else if (fit == PM_FIT_SPLINE) { + // Fit a spline to the old overscan vector. + mySpline = (psSpline1D *) + fitSpec; + // XXX: What do we do? We don't have 2-D splines. + return(in); + } + + //Subtract the polynomially fitted image from the original image + if (binFactor <= 0) { + for (int row = 0; row < origImage->numRows ; row++) { + for (int col = 0; col < origImage->numCols ; col++) { + origImage->data.F32[row][col]-= binnedImage->data.F32[row][col]; + } + } + } else { + for (int row = 0; row < origImage->numRows ; row++) { + for (int col = 0; col < origImage->numCols ; col++) { + int binRow = row / binFactor; + int binCol = col / binFactor; + origImage->data.F32[row][col]-= binnedImage->data.F32[binRow][binCol]; + } + } + + } + return(in); + }