Index: /trunk/psModules/src/detrend/Makefile.am =================================================================== --- /trunk/psModules/src/detrend/Makefile.am (revision 10298) +++ /trunk/psModules/src/detrend/Makefile.am (revision 10299) @@ -9,5 +9,5 @@ pmMaskBadPixels.c \ pmNonLinear.c \ - pmSubtractBias.c \ + pmBias.c \ pmDetrendDB.c \ pmShutterCorrection.c \ @@ -15,5 +15,5 @@ psIOBuffer.c -# pmSubtractSky.c +# pmSkySubtract.c pkginclude_HEADERS = \ @@ -23,5 +23,5 @@ pmMaskBadPixels.h \ pmNonLinear.h \ - pmSubtractBias.h \ + pmBias.h \ pmDetrendDB.h \ pmShutterCorrection.h \ @@ -29,5 +29,5 @@ psIOBuffer.h -# pmSubtractSky.h +# pmSkySubtract.h CLEANFILES = *~ Index: /trunk/psModules/src/detrend/pmBias.c =================================================================== --- /trunk/psModules/src/detrend/pmBias.c (revision 10299) +++ /trunk/psModules/src/detrend/pmBias.c (revision 10299) @@ -0,0 +1,374 @@ +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmFPA.h" +#include "pmSubtractBias.h" + + +static 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 + + // Check parities + int xIpar = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.XPARITY"); + int xSpar = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.XPARITY"); + if (xIpar != xSpar) { + psError(PS_ERR_UNKNOWN, true, "images for subtraction do not have the same " + "CELL.XPARITY (%d vs %d).\n pmSubtractBias must be upgraded to handle this situation\n", + xIpar, xSpar); + return false; + } + + int yIpar = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.YPARITY"); + int ySpar = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.YPARITY"); + if (yIpar != ySpar) { + psError(PS_ERR_UNKNOWN, true, "images for subtraction do not have the same " + "CELL.YPARITY (%d vs %d).\n pmSubtractBias must be upgraded to handle this situation\n", + xIpar, xSpar); + return false; + } + + // 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; +} + +// 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 + ) +{ + assert(overscanOpts); + assert(pixels); + assert(myStats); + + psStatsOptions statistic = psStatsSingleOption(myStats->options); // Statistic to use + assert(statistic != 0); + + // 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); + reduced->data.F32[i] = psStatsGetValue(myStats, statistic); + } 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; +} + + +bool pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, + const pmReadout *bias, const pmReadout *dark) +{ + psTrace("psModules.detrend", 4, + "---- pmSubtractBias() begin ----\n"); + PS_ASSERT_PTR_NON_NULL(in, NULL); + PS_ASSERT_IMAGE_NON_NULL(in->image, NULL); + PS_ASSERT_IMAGE_TYPE(in->image, PS_TYPE_F32, NULL); + if (bias) { + PS_ASSERT_IMAGE_NON_NULL(bias->image, NULL); + PS_ASSERT_IMAGE_TYPE(bias->image, PS_TYPE_F32, NULL); + } + if (dark) { + PS_ASSERT_IMAGE_NON_NULL(dark->image, NULL); + PS_ASSERT_IMAGE_TYPE(dark->image, PS_TYPE_F32, NULL); + } + + // psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now, used for reporting + + 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 false; + } + + psList *overscans = in->bias; // List of the overscan images + + psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use + if (statistic == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", + overscanOpts->stat); + return false; + } + psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original + + // Reduce all overscan pixels to a single value + if (overscanOpts->single) { + psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); + 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); + 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(stats, pixels, NULL, NULL, 0); + psFree(pixels); + double reduced = psStatsGetValue(stats, statistic); // Result of statistics + (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); + 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 = PS_MAX(0,diff); i < PS_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, stats); + psFree(pixels); + if (! reduced) { + return false; + } + + // 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); + 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 = PS_MAX(0,diff); i < PS_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, stats); + psFree(pixels); + if (! reduced) { + return false; + } + + // 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(stats); + } // End of overscan subtraction + + // Bias frame subtraction + if (bias) { + if (!SubtractFrame(in, bias, 1.0)) + return false; + } + + if (dark) { + // Get the scaling + float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); + float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); + if (!SubtractFrame(in, dark, inTime/darkTime)) + return false; + } + + return true; +} + + Index: /trunk/psModules/src/detrend/pmBias.h =================================================================== --- /trunk/psModules/src/detrend/pmBias.h (revision 10299) +++ /trunk/psModules/src/detrend/pmBias.h (revision 10299) @@ -0,0 +1,68 @@ +/// @file pmSubtractBias.h +/// +/// @brief Subtract the overscan, bias and dark +/// +/// @ingroup Detrend +/// +/// @author George Gusciora, MHPCC +/// @author Paul Price, IfA +/// +/// @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2006-11-30 00:52:21 $ +/// +/// Copyright 2004--2006 Institute for Astronomy, University of Hawaii +/// + +#ifndef PM_SUBTRACT_BIAS_H +#define PM_SUBTRACT_BIAS_H + +#include + +#include "pmFPA.h" + +/// Type of fit to perform +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; + +/// Options for overscan subtraction +/// +/// The overscan subtraction may be performed by reducing all overscan regions to a single value (e.g., if +/// there is no structure); or the overscan may be fit perpendicular to the read direction (usually the +/// columns) with a particular functional form; or a single value may be subtracted for each read/scan without +/// fitting (if the structure defies characterisation). In any case, statistics are required to reduce +/// multiple values to a single value (either for the scan, or for the entire overscan regions). +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; + +/// Allocator for overscan options +pmOverscanOptions *pmOverscanOptionsAlloc(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 splines + psStats *stat ///< Statistic to use + ); + +/// Subtract the overscan, bias and/or dark +/// +/// Subtracts the overscan, as measured from the bias member of the input readout (if options are non-NULL), +/// bias (if non-NULL) and dark (if non-NULL) scaled by the CELL.DARKTIME concept. +bool pmSubtractBias(pmReadout *in, ///< Input readout, to be overscan/bias/dark corrected + pmOverscanOptions *overscanOpts, ///< Options for overscan subtraction, or NULL + const pmReadout *bias, ///< Bias image to subtract, or NULL + const pmReadout *dark ///< Dark image to scale and subtract, or NULL + ); + +#endif Index: /trunk/psModules/src/detrend/pmSkySubtract.c =================================================================== --- /trunk/psModules/src/detrend/pmSkySubtract.c (revision 10299) +++ /trunk/psModules/src/detrend/pmSkySubtract.c (revision 10299) @@ -0,0 +1,725 @@ +/** @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-11-30 00:52:21 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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. + *****************************************************************************/ +#if 0 +static psImage *binImage(psImage *origImage, + int binFactor, + psStatsOptions statOptions) +{ + psTrace("psModules.detrend", 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("psModules.detrend", 4, "Exiting binImage(%d)\n", binFactor); + return(origImage); +} +#endif + +/****************************************************************************** +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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend") >= 10) { + for (i=0; i < localPolyTerms ; i++) { + printf("x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]); + } + } + + psTrace("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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); + psVector *outPerm = NULL; + + // + // Initialize A matrix and B vector. + // + psImageInit(A, 0.0); + psVectorInit(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); + 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("psModules.detrend", 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("psModules.detrend", 4, + "Exiting ImageFitPolynomial()\n"); + // psTrace("psModules.detrend", 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("psModules.detrend", 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); + psImageInit(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); + psImageInit(binnedMaskImage, 0); + } + psTrace("psModules.detrend", 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. + psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + if (!psImageStats(myStats, binnedImage, NULL, 0)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Couldn't get statistics for image.\n"); + return NULL; + } + psF64 binnedMean = myStats->sampleMean; + psF64 binnedStdev = myStats->sampleStdev; + psFree(myStats); + psTrace("psModules.detrend", 6, + "binned StDev is %f\n", binnedStdev); + + // Clip all pixels which are more than clipSD sigmas from the mean. + psTrace("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 4, + "---- pmSubtractSky() exit successfully ----\n"); + return(in); +} Index: /trunk/psModules/src/detrend/pmSkySubtract.h =================================================================== --- /trunk/psModules/src/detrend/pmSkySubtract.h (revision 10299) +++ /trunk/psModules/src/detrend/pmSkySubtract.h (revision 10299) @@ -0,0 +1,36 @@ +/** @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-11-30 00:52:21 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifndef PM_SUBTRACT_SKY_H +#define PM_SUBTRACT_SKY_H + +#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 Index: unk/psModules/src/detrend/pmSubtractBias.c =================================================================== --- /trunk/psModules/src/detrend/pmSubtractBias.c (revision 10298) +++ (revision ) @@ -1,374 +1,0 @@ -#if HAVE_CONFIG_H -#include -#endif - -#include -#include -#include - -#include "pmFPA.h" -#include "pmSubtractBias.h" - - -static 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 - - // Check parities - int xIpar = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.XPARITY"); - int xSpar = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.XPARITY"); - if (xIpar != xSpar) { - psError(PS_ERR_UNKNOWN, true, "images for subtraction do not have the same " - "CELL.XPARITY (%d vs %d).\n pmSubtractBias must be upgraded to handle this situation\n", - xIpar, xSpar); - return false; - } - - int yIpar = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.YPARITY"); - int ySpar = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.YPARITY"); - if (yIpar != ySpar) { - psError(PS_ERR_UNKNOWN, true, "images for subtraction do not have the same " - "CELL.YPARITY (%d vs %d).\n pmSubtractBias must be upgraded to handle this situation\n", - xIpar, xSpar); - return false; - } - - // 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; -} - -// 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 - ) -{ - assert(overscanOpts); - assert(pixels); - assert(myStats); - - psStatsOptions statistic = psStatsSingleOption(myStats->options); // Statistic to use - assert(statistic != 0); - - // 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); - reduced->data.F32[i] = psStatsGetValue(myStats, statistic); - } 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; -} - - -bool pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, - const pmReadout *bias, const pmReadout *dark) -{ - psTrace("psModules.detrend", 4, - "---- pmSubtractBias() begin ----\n"); - PS_ASSERT_PTR_NON_NULL(in, NULL); - PS_ASSERT_IMAGE_NON_NULL(in->image, NULL); - PS_ASSERT_IMAGE_TYPE(in->image, PS_TYPE_F32, NULL); - if (bias) { - PS_ASSERT_IMAGE_NON_NULL(bias->image, NULL); - PS_ASSERT_IMAGE_TYPE(bias->image, PS_TYPE_F32, NULL); - } - if (dark) { - PS_ASSERT_IMAGE_NON_NULL(dark->image, NULL); - PS_ASSERT_IMAGE_TYPE(dark->image, PS_TYPE_F32, NULL); - } - - // psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now, used for reporting - - 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 false; - } - - psList *overscans = in->bias; // List of the overscan images - - psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use - if (statistic == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", - overscanOpts->stat); - return false; - } - psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original - - // Reduce all overscan pixels to a single value - if (overscanOpts->single) { - psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); - 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); - 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(stats, pixels, NULL, NULL, 0); - psFree(pixels); - double reduced = psStatsGetValue(stats, statistic); // Result of statistics - (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); - 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 = PS_MAX(0,diff); i < PS_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, stats); - psFree(pixels); - if (! reduced) { - return false; - } - - // 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); - 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 = PS_MAX(0,diff); i < PS_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, stats); - psFree(pixels); - if (! reduced) { - return false; - } - - // 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(stats); - } // End of overscan subtraction - - // Bias frame subtraction - if (bias) { - if (!SubtractFrame(in, bias, 1.0)) - return false; - } - - if (dark) { - // Get the scaling - float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); - float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); - if (!SubtractFrame(in, dark, inTime/darkTime)) - return false; - } - - return true; -} - - Index: unk/psModules/src/detrend/pmSubtractBias.h =================================================================== --- /trunk/psModules/src/detrend/pmSubtractBias.h (revision 10298) +++ (revision ) @@ -1,68 +1,0 @@ -/// @file pmSubtractBias.h -/// -/// @brief Subtract the overscan, bias and dark -/// -/// @ingroup Detrend -/// -/// @author George Gusciora, MHPCC -/// @author Paul Price, IfA -/// -/// @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ -/// @date $Date: 2006-10-24 03:24:33 $ -/// -/// Copyright 2004--2006 Institute for Astronomy, University of Hawaii -/// - -#ifndef PM_SUBTRACT_BIAS_H -#define PM_SUBTRACT_BIAS_H - -#include - -#include "pmFPA.h" - -/// Type of fit to perform -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; - -/// Options for overscan subtraction -/// -/// The overscan subtraction may be performed by reducing all overscan regions to a single value (e.g., if -/// there is no structure); or the overscan may be fit perpendicular to the read direction (usually the -/// columns) with a particular functional form; or a single value may be subtracted for each read/scan without -/// fitting (if the structure defies characterisation). In any case, statistics are required to reduce -/// multiple values to a single value (either for the scan, or for the entire overscan regions). -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; - -/// Allocator for overscan options -pmOverscanOptions *pmOverscanOptionsAlloc(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 splines - psStats *stat ///< Statistic to use - ); - -/// Subtract the overscan, bias and/or dark -/// -/// Subtracts the overscan, as measured from the bias member of the input readout (if options are non-NULL), -/// bias (if non-NULL) and dark (if non-NULL) scaled by the CELL.DARKTIME concept. -bool pmSubtractBias(pmReadout *in, ///< Input readout, to be overscan/bias/dark corrected - pmOverscanOptions *overscanOpts, ///< Options for overscan subtraction, or NULL - const pmReadout *bias, ///< Bias image to subtract, or NULL - const pmReadout *dark ///< Dark image to scale and subtract, or NULL - ); - -#endif Index: unk/psModules/src/detrend/pmSubtractSky.c =================================================================== --- /trunk/psModules/src/detrend/pmSubtractSky.c (revision 10298) +++ (revision ) @@ -1,725 +1,0 @@ -/** @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.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-24 22:55:05 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - * - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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. - *****************************************************************************/ -#if 0 -static psImage *binImage(psImage *origImage, - int binFactor, - psStatsOptions statOptions) -{ - psTrace("psModules.detrend", 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("psModules.detrend", 4, "Exiting binImage(%d)\n", binFactor); - return(origImage); -} -#endif - -/****************************************************************************** -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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend") >= 10) { - for (i=0; i < localPolyTerms ; i++) { - printf("x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]); - } - } - - psTrace("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 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); - psVector *outPerm = NULL; - - // - // Initialize A matrix and B vector. - // - psImageInit(A, 0.0); - psVectorInit(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); - 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("psModules.detrend", 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("psModules.detrend", 4, - "Exiting ImageFitPolynomial()\n"); - // psTrace("psModules.detrend", 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("psModules.detrend", 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); - psImageInit(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); - psImageInit(binnedMaskImage, 0); - } - psTrace("psModules.detrend", 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. - psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - if (!psImageStats(myStats, binnedImage, NULL, 0)) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Couldn't get statistics for image.\n"); - return NULL; - } - psF64 binnedMean = myStats->sampleMean; - psF64 binnedStdev = myStats->sampleStdev; - psFree(myStats); - psTrace("psModules.detrend", 6, - "binned StDev is %f\n", binnedStdev); - - // Clip all pixels which are more than clipSD sigmas from the mean. - psTrace("psModules.detrend", 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("psModules.detrend", 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("psModules.detrend", 4, - "---- pmSubtractSky() exit successfully ----\n"); - return(in); -} Index: unk/psModules/src/detrend/pmSubtractSky.h =================================================================== --- /trunk/psModules/src/detrend/pmSubtractSky.h (revision 10298) +++ (revision ) @@ -1,36 +1,0 @@ -/** @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.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-18 02:29:15 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifndef PM_SUBTRACT_SKY_H -#define PM_SUBTRACT_SKY_H - -#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 Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 10298) +++ /trunk/psModules/src/psmodules.h (revision 10299) @@ -43,10 +43,10 @@ #include #include -#include +#include #include #include #include #include -// #include +// #include // the following headers are from psModule:astrom