Index: /trunk/psModules/src/pmReadoutCombine.c =================================================================== --- /trunk/psModules/src/pmReadoutCombine.c (revision 2583) +++ /trunk/psModules/src/pmReadoutCombine.c (revision 2584) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-05 04:49:47 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-01 21:20:24 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -55,5 +55,4 @@ PS_PTR_CHECK_NULL(inputs, NULL); PS_PTR_CHECK_NULL(params, NULL); - unsigned int maskVal = params->maskVal; psStats *stats = params->stats; int i; @@ -64,8 +63,7 @@ int minInputRows = MAX_INT; psListElem *tmpInput = NULL; + psReadout *tmpReadout = NULL; int numInputs = 0; - psReadout *tmpReadout = NULL; - float tmpF; - + int tmpI; // @@ -91,15 +89,15 @@ minInputRows = PS_MIN(minInputRows, (tmpReadout->row0 + tmpReadout->image->row0)); - tmpF = tmpReadout->row0 + + tmpI = tmpReadout->row0 + tmpReadout->image->row0 + tmpReadout->image->numRows; - maxInputRows = PS_MAX(maxInputRows, tmpF); + maxInputRows = PS_MAX(maxInputRows, tmpI); minInputCols = PS_MIN(minInputCols, (tmpReadout->col0 + tmpReadout->image->col0)); - tmpF = tmpReadout->col0 + + tmpI = tmpReadout->col0 + tmpReadout->image->col0 + tmpReadout->image->numCols; - maxInputCols = PS_MAX(maxInputCols, tmpF); + maxInputCols = PS_MAX(maxInputCols, tmpI); tmpInput = tmpInput->next; numInputs++; @@ -170,5 +168,5 @@ // For each input readout, we create a pointer to that readout in // "tmpReadouts[]", and we store the min/max pixel indices for that - // readout, in putput image coordinates, in the psVectors + // readout, in output image coordinates, in the psVectors // (outRowLower, outColLower, outRowUpper, outColUpper). i = 0; @@ -184,8 +182,4 @@ tmpReadouts[i]->image->col0 + tmpReadouts[i]->image->numCols; - - // printf("Bounds: [%d][%d] to [%d][%d]\n", outRowLower->data.U32[i], outColLower->data.U32[i], outRowUpper->data.U32[i], outColUpper->data.U32[i]); - - tmpInput = tmpInput->next; i++; @@ -203,5 +197,5 @@ for (j = output->col0; j < (output->col0 + output->numCols) ; j++) { for (int r = 0; r < numInputs ; r++) { - // printf("[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outRowLower->data.U32[r], outColLower->data.U32[r], outRowUpper->data.U32[r], outColUpper->data.U32[r]); + // printf("[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outRowLower->data.U32[r], outColLower->data.U32[r], outRowUpper->data.U32[r], outColUpper->data.U32[r]); if ((outRowLower->data.U32[r] <= i) && (outColLower->data.U32[r] <= j) && @@ -215,5 +209,5 @@ if ((NULL == tmpReadouts[r]->mask) || - !(maskVal && tmpReadouts[r]->mask->data.U8[imageRow][imageCol])) { + !(params->maskVal && tmpReadouts[r]->mask->data.U8[imageRow][imageCol])) { tmpPixels->data.F32[r] = tmpReadouts[r]->image->data.F32[imageRow][imageCol]; tmpPixelMask->data.U8[r] = 0; @@ -226,5 +220,5 @@ tmpPixelMask->data.U8[r] = 1; } - // printf("readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[r]); + // printf("readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[r]); } // At this point, we have scanned all input readouts for this Index: /trunk/psModules/src/pmSubtractSky.c =================================================================== --- /trunk/psModules/src/pmSubtractSky.c (revision 2584) +++ /trunk/psModules/src/pmSubtractSky.c (revision 2584) @@ -0,0 +1,244 @@ +/** @file pmSubtractSky.c + * + * This file will contain a module which will create a model of the + * background sky and subtract that from the input image. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-01 21:20:24 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ +#include +#include +#include "pslib.h" +#include "psConstants.h" + +int p_psDetermineNumBits(unsigned int data) +{ + int i; + unsigned int tmpData = data; + int numBits = 0; + + for (i=0;i> 1; + } + return(numBits); +} + +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); + } Index: /trunk/psModules/src/pmSubtractSky.h =================================================================== --- /trunk/psModules/src/pmSubtractSky.h (revision 2584) +++ /trunk/psModules/src/pmSubtractSky.h (revision 2584) @@ -0,0 +1,34 @@ +/** @file pmSubtractSky.h + * + * This file will contain a module which will create a model of the + * background sky and subtract that from the input image. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-01 21:20:24 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if !defined(PM_SUBTRACT_SKY_H) +#define PM_SUBTRACT_SKY_H + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +// XXX: have Desonia do this correctly. +#include "../../psLib/src/dataManip/psConstants.h" + +psImage *pmSubtractSky(psReadout *in, + void *fitSpec, + psFit fit, + int binFactor, + psStats *stats, + float clipSD); +#endif