Index: trunk/psModules/src/detrend/Makefile.am =================================================================== --- trunk/psModules/src/detrend/Makefile.am (revision 6999) +++ trunk/psModules/src/detrend/Makefile.am (revision 7018) @@ -8,5 +8,7 @@ pmFringeStats.c \ pmMaskBadPixels.c \ - pmNonLinear.c + pmNonLinear.c \ + pmSubtractBias.c \ + pmSubtractSky.c psmoduleincludedir = $(includedir) @@ -18,5 +20,7 @@ pmMaskBadPixelsErrors.h \ pmMaskBadPixels.h \ - pmNonLinear.h + pmNonLinear.h \ + pmSubtractBias.c \ + pmSubtractSky.c EXTRA_DIST = pmFlatFieldErrors.dat pmMaskBadPixelsErrors.dat Index: trunk/psModules/src/detrend/pmSubtractBias.c =================================================================== --- trunk/psModules/src/detrend/pmSubtractBias.c (revision 7018) +++ trunk/psModules/src/detrend/pmSubtractBias.c (revision 7018) @@ -0,0 +1,577 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** @file pmSubtractBias.c + * + * This file will contain a module which will subtract the detector bias + * in place from an input image. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-01 01:56:29 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include "pmSubtractBias.h" + +#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 +#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 + + +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +#define MIN(a,b) ((a) < (b) ? (a) : (b)) + + +// XXX: put these in psConstants.h +void PS_POLY1D_PRINT(psPolynomial1D *poly) +{ + printf("-------------- PS_POLY1D_PRINT() --------------\n"); + printf("poly->nX is %d\n", poly->nX); + for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) { + printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]); + } +} + +void PS_PRINT_SPLINE(psSpline1D *mySpline) +{ + printf("-------------- PS_PRINT_SPLINE() --------------\n"); + printf("mySpline->n is %d\n", mySpline->n); + for (psS32 i = 0 ; i < mySpline->n ; i++) { + PS_POLY1D_PRINT(mySpline->spline[i]); + } + PS_VECTOR_PRINT_F32(mySpline->knots); +} + +#define PS_IMAGE_PRINT_F32_HIDEF(NAME) \ +printf("======== printing %s ========\n", #NAME); \ +for (int i = 0 ; i < (NAME)->numRows ; i++) { \ + for (int j = 0 ; j < (NAME)->numCols ; j++) { \ + printf("%.5f ", (NAME)->data.F32[i][j]); \ + } \ + printf("\n"); \ +}\ + + +void overscanOptionsFree(pmOverscanOptions *options) +{ + psFree(options->stat); + psFree(options->poly); + psFree(options->spline); +} + +pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat) +{ + pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); + psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree); + + // Inputs + opts->single = single; + opts->fitType = fitType; + opts->order = order; + opts->stat = psMemIncrRefCounter(stat); + + // Outputs + opts->poly = NULL; + opts->spline = NULL; + + return opts; +} + + +/****************************************************************************** +psSubtractFrame(): this routine will take as input a readout for the input +image and a readout for the bias image. The bias image is subtracted in +place from the input image. +*****************************************************************************/ +static bool SubtractFrame(pmReadout *in,// Input readout + const pmReadout *sub, // Readout to be subtracted from input + float scale // Scale to apply before subtracting + ) +{ + assert(in); + assert(sub); + + psImage *inImage = in->image; // The input image + psImage *inMask = in->mask; // The input mask + psImage *subImage = sub->image; // The image to be subtracted + psImage *subMask = sub->mask; // The mask for the subtraction image + + // Offsets of the cells + int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0"); + int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0"); + int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0"); + int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0"); + + if ((inImage->numCols + x0in - x0sub) > subImage->numCols) { + psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n"); + return false; + } + if ((inImage->numRows + y0in - y0sub) > subImage->numRows) { + psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n"); + return false; + } + + if (scale == 1.0) { + for (int i = 0; i < inImage->numRows; i++) { + for (int j = 0; j < inImage->numCols; j++) { + inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub]; + if (inMask && subMask) { + inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; + } + } + } + } else { + for (int i = 0; i < inImage->numRows; i++) { + for (int j = 0; j < inImage->numCols; j++) { + inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale; + if (inMask && subMask) { + inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; + } + } + } + } + + return true; +} + + +#if 0 +/****************************************************************************** +ImageSubtractScalar(): subtract a scalar from the input image. + +XXX: Use a psLib function for this. + +XXX: This should + *****************************************************************************/ +static psImage *ImageSubtractScalar(psImage *image, + psF32 scalar) +{ + for (psS32 i=0;inumRows;i++) { + for (psS32 j=0;jnumCols;j++) { + image->data.F32[i][j]-= scalar; + } + } + return(image); +} +#endif + +/****************************************************************************** +GenNewStatOptions(): this routine will take as input the options member of the +stat data structure, determine if multiple options have been specified, issue +a warning message if so, and return the highest priority option (according to +the order of the if-statements in this code). The higher priority options are +listed lower in the code. + *****************************************************************************/ +static psStatsOptions GenNewStatOptions(const psStats *stat) +{ + psS32 numOptions = 0; + psStatsOptions opt = 0; + + if (stat->options & PS_STAT_ROBUST_MEDIAN) { + if (numOptions == 0) { + opt = PS_STAT_ROBUST_MEDIAN; + } + numOptions++; + } + + if (stat->options & PS_STAT_CLIPPED_MEAN) { + if (numOptions == 0) { + opt = PS_STAT_CLIPPED_MEAN; + } + numOptions++; + } + + if (stat->options & PS_STAT_SAMPLE_MEDIAN) { + if (numOptions == 0) { + opt = PS_STAT_SAMPLE_MEDIAN; + } + numOptions++; + } + + if (stat->options & PS_STAT_SAMPLE_MEAN) { + numOptions++; + opt = PS_STAT_SAMPLE_MEAN; + } + + + if (numOptions == 0) { + psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n"); + } + if (numOptions != 1) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.c: GenNewStatOptions(): Too many statistics options have been specified\n"); + } + return(opt); +} + + + +#if 0 +/****************************************************************************** +ScaleOverscanVector(): this routine takes as input an arbitrary vector, +creates a new vector of length n, and fills the new vector with the +interpolated values of the old vector. The type of interpolation is: + PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data. + PM_FIT_SPLINE: fit splines to the input vector data. +XXX: Doesn't it make more sense to do polynomial interpolation on a few +elements of the input vector, rather than fit a polynomial to the entire +vector? + *****************************************************************************/ +static psVector *ScaleOverscanVector(psVector *overscanVector, + psS32 n, + void *fitSpec, + pmFit fit) +{ + psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, + "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n); + // PS_VECTOR_PRINT_F32(overscanVector); + + if (NULL == overscanVector) { + return(overscanVector); + } + + // Allocate the new vector. + psVector *newVec = psVectorAlloc(n, PS_TYPE_F32); + + // + // If the new vector is the same size as the old, simply copy the data. + // + if (n == overscanVector->n) { + for (psS32 i = 0 ; i < n ; i++) { + newVec->data.F32[i] = overscanVector->data.F32[i]; + } + return(newVec); + } + psPolynomial1D *myPoly; + psSpline1D *mySpline; + psF32 x; + psS32 i; + if (fit == PM_FIT_POLYNOMIAL) { + // Fit a polynomial to the old overscan vector. + myPoly = (psPolynomial1D *) fitSpec; + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); + if (myPoly == NULL) { + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n"); + return(NULL); + } + + // For each element of the new vector, convert the x-ordinate to that + // of the old vector, use the fitted polynomial to determine the + // interpolated value at that point, and set the new vector. + for (i=0;in) / ((psF32) n); + newVec->data.F32[i] = psPolynomial1DEval(myPoly, x); + } + } else if (fit == PM_FIT_SPLINE) { + psS32 mustFreeSpline = 0; + // Fit a spline to the old overscan vector. + mySpline = (psSpline1D *) fitSpec; + // XXX: Does it make any sense to have a psSpline argument? + if (mySpline == NULL) { + mustFreeSpline = 1; + } + + // + // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL, + // splines endpoints will be from 0.0 to overscanVector->n-1. Must scale + // properly when doing the spline eval. + // + // mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); + mySpline = psVectorFitSpline1D(NULL, overscanVector); + if (mySpline == NULL) { + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n"); + return(NULL); + } + // PS_PRINT_SPLINE(mySpline); + + // For each element of the new vector, convert the x-ordinate to that + // of the old vector, use the fitted polynomial to determine the + // interpolated value at that point, and set the new vector. + for (i=0;in - 1] + x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n); + newVec->data.F32[i] = psSpline1DEval(mySpline, x); + } + if (mustFreeSpline ==1) { + psFree(mySpline); + } + // PS_VECTOR_PRINT_F32(newVec); + + + } else { + psError(PS_ERR_UNKNOWN, true, "unknown fit type. Returning NULL.\n"); + psFree(newVec); + return(NULL); + } + + psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, + "---- ScaleOverscanVector() exit ----\n"); + return(newVec); +} + +#endif + +// Produce an overscan vector from an array of pixels +static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options + const psArray *pixels, // Array of vectors containing the pixel values + psStats *myStats // Statistic to use in reducing the overscan + ) +{ + // Reduce the overscans + psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row + psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate + psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting + for (int i = 0; i < pixels->n; i++) { + psVector *values = pixels->data[i]; // Vector with overscan values + if (values->n > 0) { + mask->data.U8[i] = 0; + ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] + psVectorStats(myStats, values, NULL, NULL, 0); + double reducedVal = NAN; // Result of statistics + if (! p_psGetStatValue(myStats, &reducedVal)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result " + "of statistics on row %d.\n", i); + return NULL; + } + reduced->data.F32[i] = reducedVal; + } else if (overscanOpts->fitType == PM_FIT_NONE) { + psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the " + "image, and no fit is requested.\n"); + return NULL; + } else { + // We'll fit this one out + mask->data.U8[i] = 1; + } + } + + // Fit the overscan, if required + switch (overscanOpts->fitType) { + case PM_FIT_NONE: + // No fitting --- that's easy. + break; + case PM_FIT_POLY_ORD: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (! overscanOpts->poly) { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + psFree(reduced); + reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_POLY_CHEBY: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (! overscanOpts->poly) { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + psFree(reduced); + reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_SPLINE: + // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an + // input spline + overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); + psFree(reduced); + reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); + return NULL; + break; + } + + psFree(ordinate); + psFree(mask); + + return reduced; +} + + + +/****************************************************************************** +XXX: The SDRS does not specify type support. F32 is implemented here. + *****************************************************************************/ +pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, + const pmReadout *bias, const pmReadout *dark) +{ + psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, + "---- pmSubtractBias() begin ----\n"); + PS_ASSERT_READOUT_NON_NULL(in, NULL); + PS_ASSERT_READOUT_NON_EMPTY(in, NULL); + PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); + + psImage *image = in->image; // The input image + + // Overscan processing + if (overscanOpts) { + // Check for an unallowable pmFit. + if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD && + overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) { + psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", overscanOpts->fitType); + return(in); + } + + psList *overscans = in->bias; // List of the overscan images + + psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original + myStats->options = GenNewStatOptions(overscanOpts->stat); + + // Reduce all overscan pixels to a single value + if (overscanOpts->single) { + psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); + pixels->n = 0; + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) { + int index = pixels->n; // Index + pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); + // XXX Reimplement with memcpy + for (int i = 0; i < overscan->numRows; i++) { + for (int j = 0; j < overscan->numCols; j++) { + pixels->data.F32[index++] = overscan->data.F32[i][j]; + } + } + + } + psFree(iter); + + (void)psVectorStats(myStats, pixels, NULL, NULL, 0); + double reduced = NAN; // Result of statistics + if (! p_psGetStatValue(myStats, &reduced)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning input image.\n"); + return(in); + } + (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32)); + } else { + + // We do the regular overscan subtraction + + bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction + + if (readRows) { + // The read direction is rows + psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) { + psVector *values = psVectorAlloc(0, PS_TYPE_F32); + values->n = 0; + pixels->data[i] = values; + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) { + int diff = image->row0 - overscan->row0; // Offset between the two regions + for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) { + // i is row on overscan + // XXX Reimplement with memcpy + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numCols); + for (int j = 0; j < overscan->numCols; j++) { + values->data.F32[index++] = overscan->data.F32[i][j]; + } + values->n += overscan->numCols; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = overscanVector(overscanOpts, pixels, myStats); + psFree(pixels); + if (! reduced) { + return in; + } + + // Subtract row by row + for (int i = 0; i < image->numRows; i++) { + for (int j = 0; j < image->numCols; j++) { + image->data.F32[i][j] -= reduced->data.F32[i]; + } + } + psFree(reduced); + + } else { + // The read direction is columns + psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) { + psVector *values = psVectorAlloc(0, PS_TYPE_F32); + values->n = 0; + pixels->data[i] = values; + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) { + int diff = image->col0 - overscan->col0; // Offset between the two regions + for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) { + // i is column on overscan + // XXX Reimplement with memcpy + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numRows); + for (int j = 0; j < overscan->numRows; j++) { + values->data.F32[index++] = overscan->data.F32[i][j]; + } + values->n += overscan->numRows; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = overscanVector(overscanOpts, pixels, myStats); + psFree(pixels); + if (! reduced) { + return in; + } + + // Subtract column by column + for (int i = 0; i < image->numCols; i++) { + for (int j = 0; j < image->numRows; j++) { + image->data.F32[j][i] -= reduced->data.F32[i]; + } + } + psFree(reduced); + } + } + psFree(myStats); + } // End of overscan subtraction + + // Bias frame subtraction + if (bias) { + SubtractFrame(in, bias, 1.0); + } + + if (dark) { + // Get the scaling + float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); + float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); + SubtractFrame(in, dark, inTime/darkTime); + } + + return in; +} + + Index: trunk/psModules/src/detrend/pmSubtractBias.h =================================================================== --- trunk/psModules/src/detrend/pmSubtractBias.h (revision 7018) +++ trunk/psModules/src/detrend/pmSubtractBias.h (revision 7018) @@ -0,0 +1,77 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** @file pmSubtractBias.h + * + * This file will contain a module which will subtract the detector bias + * in place from an input image. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-01 01:56:29 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if !defined(PM_SUBTRACT_BIAS_H) +#define PM_SUBTRACT_BIAS_H + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" + +#include "pmFPA.h" + +typedef enum { + PM_OVERSCAN_NONE, ///< No overscan subtraction + PM_OVERSCAN_EDGE, ///< Subtract the statistic of pixels along the to-be-determined readout direction + PM_OVERSCAN_ROWS, ///< Subtract rows + PM_OVERSCAN_COLUMNS, ///< Subtract columns + PM_OVERSCAN_ALL ///< Subtract the statistic of all pixels in overscan region +} pmOverscanAxis; + +typedef enum { + PM_FIT_NONE, ///< No fit + PM_FIT_POLY_ORD, ///< Fit ordinary polynomial + PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial + PM_FIT_SPLINE ///< Fit cubic splines +} pmFit; + +typedef struct +{ + // Inputs + bool single; // Reduce all overscan regions to a single value? + pmFit fitType; // Type of fit to overscan + unsigned int order; // Order of polynomial, or number of spline pieces + psStats *stat; // Statistic to use when reducing the minor direction + // Outputs + psPolynomial1D *poly; // Result of polynomial fit + psSpline1D *spline; // Result of spline fit +} +pmOverscanOptions; + +pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat); + +pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, + const pmReadout *bias, const pmReadout *dark); + +#if 0 +pmReadout *pmSubtractBias(pmReadout *in, ///< The input pmReadout image + void *fitSpec, ///< A polynomial or spline, defining the fit type. + const psList *overscans, ///< A psList of overscan images + pmOverscanAxis overScanAxis, ///< Defines how overscans are applied + psStats *stat, ///< The statistic to be used in combining overscan data + int nBin, ///< The amount of binning to be done image pixels. + pmFit fit, ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE + const pmReadout *bias); ///< A possibly NULL bias pmReadout which is to be subtracted +#endif + +#endif Index: trunk/psModules/src/detrend/pmSubtractSky.c =================================================================== --- trunk/psModules/src/detrend/pmSubtractSky.c (revision 7018) +++ trunk/psModules/src/detrend/pmSubtractSky.c (revision 7018) @@ -0,0 +1,741 @@ +/** @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: 2006-05-01 01:56:29 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * + * + */ + +#include +#include +#include "pslib.h" +#include "pmSubtractSky.h" + +// XXX: Get rid of the. Create pmUtils.h +psImage *p_psDetermineTrimmedImage( + pmReadout *in +); + +/****************************************************************************** +DetermineNumBits(data): This routine takes an enum psStatsOptions as an +argument and returns the number of non-zero bits. + +XXX: This code is duplicated in the ReadoutCombine file. + *****************************************************************************/ +static psS32 DetermineNumBits(psStatsOptions data) +{ + psTrace("SubtractSky.DetermineNumBits", 4, "Calling DetermineNumBits(0x%x)\n", data); + + psS32 i; + psU64 tmpData = data; + psS32 numBits = 0; + + for (i=0;i<(8 * sizeof(psStatsOptions));i++) { + if (0x0001 & tmpData) { + numBits++; + } + tmpData = tmpData >> 1; + } + + psTrace("SubtractSky.DetermineNumBits", 4, + "Calling DetermineNumBits(0x%x) -> %d\n", data, numBits); + return(numBits); +} + +/****************************************************************************** +getHighestPriorityStatOption(statOptions): this routine takes as input a +psStats->options with multiple options set and returns one with a single +option set according to the precedence set in the SDRS. + *****************************************************************************/ +static psU64 getHighestPriorityStatOption(psU64 statOptions) +{ + psTrace("SubtractSky.getHighestPriorityStatOption", 4, + "Calling getHighestPriorityStatOption(0x%x)\n", statOptions); + + if (statOptions & PS_STAT_SAMPLE_MEAN) { + return(PS_STAT_SAMPLE_MEAN); + } else if (statOptions & PS_STAT_SAMPLE_MEDIAN) { + return(PS_STAT_SAMPLE_MEDIAN); + } else if (statOptions & PS_STAT_CLIPPED_MEAN) { + return(PS_STAT_CLIPPED_MEAN); + } else if (statOptions & PS_STAT_FITTED_MEAN) { + return(PS_STAT_FITTED_MEAN); + } else if (statOptions & PS_STAT_ROBUST_MEDIAN) { + return(PS_STAT_ROBUST_MEDIAN); + } + psError(PS_ERR_UNKNOWN, true, "Unallowable option requested for statistically binning image pixels.\n"); + return(-1); + // XXX + //else if (statOptions & PS_STAT_ROBUST_MODE) { + // return(PS_STAT_ROBUST_MODE); + //} +} + +/****************************************************************************** +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. +XXX: I coded this before I was aware of a psLib reBin function. I don't +use this function in this module. I'm keeping it here in the event that +requirements change and we might need a custom reBin function. + *****************************************************************************/ +/* +static psImage *binImage(psImage *origImage, + int binFactor, + psStatsOptions statOptions) +{ + psTrace("SubtractSky.binImage", 4, "Calling binImage(%d)\n", binFactor); + + if (binFactor <= 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: binImage(): binFactor is %d\n", binFactor); + return(origImage); + } + if (binFactor == 1) { + return(origImage); + } + + psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32); + psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8); + psStats *myStats = psStatsAlloc(statOptions); + + for (psS32 row = 0; row < origImage->numRows ; row+=binFactor) { + for (psS32 col = 0; col < origImage->numCols ; col+=binFactor) { + psS32 count = 0; + for (psS32 binRow = 0; binRow <= binFactor ; binRow++) { + for (psS32 binCol = 0; binCol <= binFactor ; binCol++) { + if (((row + binRow) < origImage->numRows) && + ((col + binCol) < origImage->numCols)) { + 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++; + } + } + psStats *rc1 = psVectorStats(myStats, binVector, NULL, binMask, 1); + if (rc1 == NULL) { + psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); + return(origImage); + } + psF64 statValue; + psBool rc = p_psGetStatValue(rc1, &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(binVector); + psFree(binMask); + psFree(myStats); + + psTrace("SubtractSky.binImage", 4, "Exiting binImage(%d)\n", binFactor); + return(origImage); +} +*/ + +/****************************************************************************** +CalculatePolyTerms(xOrder, yOrder): this routine will calculate the number of +coefficients (or terms) in a 2-D polynomial of order (xOrder, yOrder). + +XXX: Use your brain and figure out the analytical expression. + +XXX: Why isn't it simply (xOrder+1) * (yOrder+1)? + *****************************************************************************/ +static psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder) +{ + psTrace("SubtractSky.CalculatePolyTerms", 4, + "Calling CalculatePolyTerms(%d, %d)\n", xOrder, yOrder); + + psS32 maxOrder = PS_MAX(xOrder, yOrder); + psS32 localPolyTerms = 0; + psS32 order = 0; + psS32 num=0; + + for (order=0;order<=maxOrder;order++) { + for (num=0;num<=order;num++) { + if (((order-num) <= xOrder) && (num <= yOrder)) { + localPolyTerms++; + } + } + } + psTrace("SubtractSky.CalculatePolyTerms", 4, + "Exiting CalculatePolyTerms(%d, %d) -> %d\n", xOrder, yOrder, localPolyTerms); + return(localPolyTerms); + + // return((xOrder+1) * (yOrder+1)); +} + +/****************************************************************************** +buildPolyTerms(): this routine computes a 2-D array polyTerms[][] that holds +terms for the polynomial that is used to model the sky background. We use +this array primarily for convenience in computations involving sky model +polynomials. It is defined as: + polyTerms[i][0] = the power to which X is raised in the i-th term of in an + poly-order sky background polynomial. + + polyTerms[i][1] = the power to which Y is raised in the i-th term of in an + poly-order sky background polynomial. + *****************************************************************************/ +static psS32 **buildPolyTerms(psS32 xOrder, psS32 yOrder) +{ + psTrace("SubtractSky.buildPolyTerms", 4, + "Calling buildPolyTerms(%d, %d)\n", xOrder, yOrder); + + psS32 i=0; + psS32 order = 0; + psS32 num=0; + psS32 localPolyTerms = CalculatePolyTerms(xOrder, yOrder); + psS32 maxOrder = PS_MAX(xOrder, yOrder); + + // Create the data structure which we hold the xy order of each coeff. + psS32 **polyTerms = (psS32 **) psAlloc(localPolyTerms * sizeof(psS32 *)); + for (i=0; i < localPolyTerms ; i++) { + polyTerms[i] = (psS32 *) psAlloc(2 * sizeof(psS32)); + } + + i=0; + // This code segment loops through each term i in the polynomial and + // calculates the power to which x/y are raised in that i-th term. + // We first do the 0-order terms, then the 1-order terms, etc. + for (order=0;order<=maxOrder;order++) { + for (num=0;num<=order;num++) { + if (((order-num) <= xOrder) && (num <= yOrder)) { + polyTerms[i][0] = order-num; + polyTerms[i][1] = num; + i++; + } + } + } + + if (psTraceGetLevel(".psModule.pmSubtractSky.buildPolyTerms") >= 10) { + for (i=0; i < localPolyTerms ; i++) { + printf("x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]); + } + } + + psTrace("SubtractSky.buildPolyTerms", 4, + "Exiting buildPolyTerms(%d, %d)\n", xOrder, yOrder); + return(polyTerms); +} + +/****************************************************************************** +This procedure calculates various combinations of powers of x and y and stores +them in the data structure p_psPolySums[][]. After it completes: + + p_psPolySums[i][j] == x^i * y^j + +XXX: Use a psImage for the p_psPolySums data structure? +XXX: p_psPolySums: should this be a global? Did you get the storage classifier + and name correct? +XXX: Use variable size arrays for p_psPolySums[][]. +XXX: Must initialize p_psPolySums[][]? + *****************************************************************************/ +#define PS_MAX_POLYNOMIAL_ORDER 20 + +psF64 p_psPolySums[PS_MAX_POLYNOMIAL_ORDER+1][PS_MAX_POLYNOMIAL_ORDER+1]; +static void buildSums(psF64 x, + psF64 y, + psS32 xOrder, + psS32 yOrder) +{ + psTrace("SubtractSky.buildPolyTerms", 4, + "Calling buildPolyTerms(%d, %d)\n", xOrder, yOrder); + + psS32 i = 0; + psS32 j = 0; + psF64 xSum = 0.0; + psF64 ySum = 0.0; + + xSum = 1.0; + ySum = 1.0; + for(i=0;i<=xOrder;i++) { + ySum = xSum; + for(j=0;j<=yOrder;j++) { + p_psPolySums[i][j] = ySum; + ySum*= y; + } + xSum*= x; + } + psTrace("SubtractSky.buildPolyTerms", 4, + "Exiting buildPolyTerms(%d, %d)\n", xOrder, yOrder); +} + +/****************************************************************************** +ImageFitPolynomial(myPoly, dataImage, maskImage): this private routine takes +an input image along with a mask and fits a polynomial to it. The degree of +the polynomial is specified by input parameter myPoly, and need not be +symmetrical in orders of X and Y. The polynomial must be type +PS_POLYNOMIAL_ORD. If there are not enough rows or columns in the input image +for the order of the polynomial, then that order is reduced. The algorithm +used in this routine is based on that of the pilot project ADD, but is not +documented anywhere. + +XXX: Different trace message facilities in use here. + *****************************************************************************/ +static psPolynomial2D *ImageFitPolynomial( + psPolynomial2D *myPoly, + psImage *dataImage, + psImage *maskImage) +{ + psTrace("SubtractSky.ImageFitPolynomial", 4, + "Calling ImageFitPolynomial()\n"); + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_POLY_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL); + PS_ASSERT_IMAGE_NON_NULL(dataImage, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(dataImage, NULL); + PS_ASSERT_IMAGE_TYPE(dataImage, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGE_NON_NULL(maskImage, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(maskImage, NULL); + PS_ASSERT_IMAGE_TYPE(maskImage, PS_TYPE_U8, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(dataImage, maskImage, NULL); + psS32 oldPolyX = -1; + psS32 oldPolyY = -1; + + // The matrix equations become singular if there are more powers of X + // in myPoly then there are rows of the image. I think. Similarly for + // powers of Y and columns. So. Here we reduce the complexity of the + // polynomial if there are not enough rows/columns in the input image. + + if ((myPoly->nX + 1) > dataImage->numRows) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in x-dimension.\n"); + oldPolyX = myPoly->nX; + myPoly->nX = dataImage->numRows - 1; + } + if ((myPoly->nY + 1) > dataImage->numCols) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in y-dimension.\n"); + oldPolyY = myPoly->nY; + myPoly->nY = dataImage->numCols - 1; + } + psS32 i; + psS32 j; + psS32 x; + psS32 y; + psS32 aRow; + psS32 aCol; + psS32 **polyTerms = buildPolyTerms(myPoly->nX, myPoly->nY); + // We determine how many coefficients will be in the polynomial that we + // are fitting to this image. + psS32 localPolyTerms = CalculatePolyTerms(myPoly->nX, myPoly->nY); + psImage *A = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64); + psImage *Aout = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64); + psVector *B = psVectorAlloc(localPolyTerms, PS_TYPE_F64); + B->n = B->nalloc; + psVector *outPerm = NULL; + + // + // Initialize A matrix and B vector. + // + PS_IMAGE_SET_F64(A, 0.0); + PS_VECTOR_SET_F64(B, 0.0); + + // + // We build the A matrix and B vector. + // + for (x=0;xnumRows;x++) { + for (y=0;ynumCols;y++) { + if (maskImage->data.U8[x][y] == 0) { + buildSums((psF64) x, (psF64) y, myPoly->nX, myPoly->nY); + + /************************************************************ + This code dervies from equation (7) of the pilot ADD. However, + it is not exactly the same in that the order of the polynomial + may be different in X And Y. + + Equation (7) from the pilot ADD describes 16 linear equations. + The i-th equation is simply the partial derivative of the + sky background polynomial (1) w.r.t. to the i-th term in + that polynomial. The i-th equation is stored in row i of + matrix A[][] (matrix A[][] has origin (1,1), not (0,0)). To + compute A[i][j] we simply multiply the j-th term of the Sky + Background Polynomial (SBP) by the i-th term of SBP. + ************************************************************/ + for (aRow=0;aRowdata.F64[aRow][aCol]+= + (p_psPolySums[ polyTerms[aCol][0] ][ polyTerms[aCol][1] ] * + p_psPolySums[ polyTerms[aRow][0] ][ polyTerms[aRow][1] ]); + } + } + // Build the B[] vector, which is the right-hand side of (7). + for (i=0;idata.F64[i]+= dataImage->data.F32[x][y] * + p_psPolySums[ polyTerms[i][0] ][ polyTerms[i][1] ]; + } + } + } + } + + if (psTraceGetLevel(".psModule.pmSubtractSky.ImageFitPolynomial") >= 8) { + for (aRow=0;aRowdata.F64[aRow][aCol]); + } + } + + for (i=0;i<=localPolyTerms;i++) { + printf("B[%d] is %f\n", i, B->data.F64[i]); + } + } + + // + // Solve the matrix equations for the polynomial coefficients C. + // XXX: How do we know if these matrix operations were successful? + // + Aout = psMatrixLUD(Aout, &outPerm, A); + PS_ASSERT_IMAGE_NON_NULL(Aout, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(Aout, NULL); + psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64); + C->n = C->nalloc; + psMatrixLUSolve(C, Aout, B, outPerm); + + // + // Set the appropriate coefficients in the myPoly structure. + // + for (i=0;icoeff[ polyTerms[i][0] ][ polyTerms[i][1] ] = C->data.F64[i]; + psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6, + "myPoly->coeff[%d][%d] is %f\n", polyTerms[i][0], polyTerms[i][1], myPoly->coeff[ polyTerms[i][0] ][ polyTerms[i][1] ]); + } + + // + // Free data structures that were allocated in this module. + // + for (i=0;inX = oldPolyX; + for (i=oldPolyX ; i < (1 + myPoly->nX) ; i++) { + for (j=0;j<(1 + myPoly->nY) ; j++) { + myPoly->coeff[i][j] = 0.0; + } + } + } + if (oldPolyY != -1) { + myPoly->nY = oldPolyY; + for (i=0 ; i < (1 + myPoly->nX) ; i++) { + for (j=oldPolyY;j < (1 + myPoly->nY) ; j++) { + myPoly->coeff[i][j] = 0.0; + } + } + } + + psTrace("SubtractSky.ImageFitPolynomial", 4, + "Exiting ImageFitPolynomial()\n"); + // psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 4, + // "---- ImageFitPolynomial() end successfully ----\n"); + return(myPoly); +} + + +/****************************************************************************** +pmReadout pmSubtractSky(): + +XXX: use static vectors for myStats, and the binned image + +XXX: The SDR is silent about types. PS_TYPE_F32 is implemented here. + +XXX: Sync the psTrace message facilities. + *****************************************************************************/ +pmReadout *pmSubtractSky(pmReadout *in, + void *fitSpec, + psFit fit, + psS32 binFactor, + psStats *stats, + psF32 clipSD) +{ + PS_ASSERT_READOUT_NON_NULL(in, NULL); + PS_ASSERT_READOUT_NON_EMPTY(in, NULL); + PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); + PS_WARN_PTR_NON_NULL(in->parent); + if (in->parent != NULL) { + PS_WARN_PTR_NON_NULL(in->parent->concepts); + } + psTrace(".psModule.pmSubtractSky", 4, + "---- pmSubtractSky() begin ----\n"); + + if ((fit != PM_FIT_NONE) && + (fit != PM_FIT_POLYNOMIAL) && + (fit != PM_FIT_SPLINE)) { + psError(PS_ERR_UNKNOWN, true, "psFit is unallowable (%d). Returning in image.\n", fit); + return(in); + } + + psStatsOptions statOptions = 0; + + // + // Return the original input readout if the fit specs are poorly defined. + // No warning or error messages should be generated. + // + if ((fitSpec == NULL) || + ((fit == PM_FIT_NONE) || (fit == PM_FIT_SPLINE))) { + // psLogMsg(__func__, PS_LOG_WARN, "Fit specs are poorly defined. Returning in image.\n"); + return(in); + } + + // + // Determine trimmed image from metadata. + // + + psImage *trimmedImg = p_psDetermineTrimmedImage(in); + psImage *binnedImage = NULL; + psPolynomial2D *myPoly = NULL; + psImage *binnedMaskImage = NULL; + psU32 oldStatOptions = 0; + + // + // Determine which statistic to use when binning pixels, if any. + // + if (stats != NULL) { + statOptions = stats->options; + if (1 < DetermineNumBits(statOptions)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Multiple statistical options have been requested.\n"); + statOptions = getHighestPriorityStatOption(statOptions); + if (statOptions == -1) { + psError(PS_ERR_UNKNOWN, true, "Not allowable stats->option was specified. Returning in image.\n"); + return(in); + } + // Save old input "stats" parameter. + oldStatOptions = stats->options; + stats->options = statOptions; + } + if (0 == DetermineNumBits(statOptions)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): no stats->options was requested\n"); + } + } + + // + // Generate required warning messages. + // + 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"); + } + + // + // Bin the input image according to input parameters. + // Create a new binned image mask. + // + if ((binFactor <= 1) || (stats == NULL) || (0 == DetermineNumBits(statOptions))) { + // No binning is required here. Simply create a copy of the image + // and a mask. + binnedImage = psImageCopy(binnedImage, trimmedImg, PS_TYPE_F32); + if (binnedImage == NULL) { + psError(PS_ERR_UNKNOWN, false, "psImageCopy() returned NULL. Returning in image.\n"); + return(in); + } + + if (in->mask != NULL) { + binnedMaskImage = psImageCopy(binnedMaskImage, in->mask, PS_TYPE_U8); + if (binnedMaskImage == NULL) { + psError(PS_ERR_UNKNOWN, false, "psImageCopy() returned NULL. Returning in image.\n"); + psFree(binnedImage); + return(in); + } + } else { + binnedMaskImage = psImageAlloc(binnedImage->numCols, + binnedImage->numRows, + PS_TYPE_U8); + PS_IMAGE_SET_U8(binnedMaskImage, 0); + } + } else { + binnedImage = psImageRebin(NULL, trimmedImg, in->mask, 0, binFactor, stats); + if (binnedImage == NULL) { + psError(PS_ERR_UNKNOWN, false, "psImageRebin() returned NULL. Returning in image.\n"); + return(in); + } + binnedMaskImage = psImageAlloc(binnedImage->numCols, + binnedImage->numRows, + PS_TYPE_U8); + PS_IMAGE_SET_U8(binnedMaskImage, 0); + } + psTrace(".psModule.pmSubtractSky", 4, + "binnedImage size is (%d, %d)\n", binnedImage->numRows, binnedImage->numCols); + + // + // 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 { + // Determine the mean and standard deviation of the binned image. + psF64 binnedMean; + psF64 binnedStdev; + psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + psStats *rc = psImageStats(myStats, binnedImage, NULL, 0); + if (rc == NULL) { + psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); + return(in); + } + if (false == p_psGetStatValue(rc, &binnedMean)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine requested statistical operation. Returning in image.\n"); + return(in); + } + psTrace(".psModule.pmSubtractSky", 6, + "binned Mean is %f\n", binnedMean); + + myStats->options = PS_STAT_SAMPLE_STDEV; + rc = psImageStats(myStats, binnedImage, NULL, 0); + if (rc == NULL) { + psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); + return(in); + } + if (false == p_psGetStatValue(myStats, &binnedStdev)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine requested statistical operation. Returning in image.\n"); + return(in); + } + psFree(myStats); + psTrace(".psModule.pmSubtractSky", 6, + "binned StDev is %f\n", binnedStdev); + + // Clip all pixels which are more than clipSD sigmas from the mean. + psTrace(".psModule.pmSubtractSky", 6, + "clipSD is %f\n", clipSD); + + for (psS32 row = 0; row < binnedImage->numRows ; row++) { + for (psS32 col = 0; col < binnedImage->numCols ; col++) { + if (fabs(binnedImage->data.F32[row][col] - binnedMean) > + (clipSD * binnedStdev)) { + binnedMaskImage->data.U8[row][col] = 1; + } + } + } + } + + // + // Fit the polynomial to the binned image + // + if (fit == PM_FIT_POLYNOMIAL) { + myPoly = (psPolynomial2D *) fitSpec; + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_POLY_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL); + + myPoly = ImageFitPolynomial(myPoly, binnedImage, binnedMaskImage); + + if (myPoly != NULL) { + // Set the pixels in the binned image to that of the polynomial. + binnedImage = psImageEvalPolynomial(binnedImage, myPoly); + if (binnedImage == NULL) { + psError(PS_ERR_UNKNOWN, false, "psImageEvalPolynomial() returned NULL. Returning in image.\n"); + psFree(binnedMaskImage); + if (!((binFactor <= 1) || (stats == NULL))) { + psFree(binnedImage); + } + if (oldStatOptions != 0) { + stats->options = statOptions; + } + return(in); + } + } else { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): could not model sky with a polynomial.\n"); + psFree(binnedMaskImage); + if (!((binFactor <= 1) || (stats == NULL))) { + psFree(binnedImage); + } + if (oldStatOptions != 0) { + stats->options = statOptions; + } + return(in); + } + } else { + // We shouldn't get here since we check this above. + psError(PS_ERR_UNKNOWN, true, "Unallowable fit type. Returning in image.\n"); + psFree(binnedMaskImage); + if (!((binFactor <= 1) || (stats == NULL))) { + psFree(binnedImage); + } + if (oldStatOptions != 0) { + stats->options = statOptions; + } + return(in); + } + + // + //Subtract the polynomially fitted image from the original image + // + if (binFactor <= 1) { + // The binned image is the same size as the original image. + for (psS32 row = 0; row < trimmedImg->numRows ; row++) { + for (psS32 col = 0; col < trimmedImg->numCols ; col++) { + trimmedImg->data.F32[row][col]-= binnedImage->data.F32[row][col]; + } + } + } else { + for (psS32 row = 0; row < trimmedImg->numRows ; row++) { + for (psS32 col = 0; col < trimmedImg->numCols ; col++) { + // We calculate the F32 value of the pixel coordinates in the + // binned image and then use a pixel interpolation routine to + // determine the value of the pixel at that location. + psF32 binRowF64 = ((psF32) row) / ((psF32) binFactor); + psF32 binColF64 = ((psF32) col) / ((psF32) binFactor); + + // We add 0.5 to the pixel locations since the pixel + // interpolation routine defines the location of pixel + // (i, j) as (i+0.5, j+0.5). + binRowF64+= 0.5; + binColF64+= 0.5; + + psF32 binPixel = (psF32) psImagePixelInterpolate( + binnedImage, binColF64, binRowF64, + NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + trimmedImg->data.F32[row][col]-= binPixel; + + psTrace(".psModule.pmSubtractSky", 8, + "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f\n", + row, col, binRowF64-0.5, binColF64-0.5, binPixel); + } + } + + } + psFree(binnedMaskImage); + psFree(binnedImage); + if (oldStatOptions != 0) { + stats->options = statOptions; + } + + psTrace(".psModule.pmSubtractSky", 4, + "---- pmSubtractSky() exit successfully ----\n"); + return(in); +} Index: trunk/psModules/src/detrend/pmSubtractSky.h =================================================================== --- trunk/psModules/src/detrend/pmSubtractSky.h (revision 7018) +++ trunk/psModules/src/detrend/pmSubtractSky.h (revision 7018) @@ -0,0 +1,40 @@ +/** @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: 2006-05-01 01:56:29 $ + * + * 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" +#include "pmFPA.h" + +// XXX: this is pmFit in pmSubtractBias.c, named psFit here. +typedef enum { + PM_FIT_NONE, ///< No fit + PM_FIT_POLYNOMIAL, ///< Fit polynomial + PM_FIT_SPLINE ///< Fit cubic splines +} psFit; + +pmReadout *pmSubtractSky(pmReadout *in, + void *fitSpec, + psFit fit, + int binFactor, + psStats *stats, + float clipSD); +#endif