Index: /trunk/psModules/src/detrend/Makefile.am =================================================================== --- /trunk/psModules/src/detrend/Makefile.am (revision 7017) +++ /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 Index: /trunk/psModules/src/imcombine/pmImageSubtract.c =================================================================== --- /trunk/psModules/src/imcombine/pmImageSubtract.c (revision 7018) +++ /trunk/psModules/src/imcombine/pmImageSubtract.c (revision 7018) @@ -0,0 +1,1407 @@ +/** @file ImageSubtract.c + * + * This file will contain code which creates a set of kernel basis + * functions, solves for their solution, and applies them to an image. + * + * @author Paul Price, IfA (original prototype) + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-01 01:56:48 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * XXX: sync with IfA on this: + * The (x, y) (row, col) issue is becoming a problem. In this file, and I + * think, the rest of psLib and psModules, the following conventions are used: + * + * 1) x will correspond to the column, and y will correspond to the row. + * 2) When used in function prototypes, the column (and hence x) appears + * first. + * 3) When used to index 2-D arrays, obviously, the row (and hence, y) + * appears first. (2 and 3 are the source of confusion). + * 4) When (u, v) are used in certain structures. + * u corresponds to x + * v corresponds to y + * 5) When element (a, b) is casually referred to in comments, or + * documentation it is unclear where a is the row, or the column. + * 6) A convention on loop index variables (i, j) would be convenient. + * Currently, sometimes i corresponds to the column (x), + * usually it corresponds to the row (y). + * + * XXX: The following variables are used an interpreted this way: + * kernelSize: Means that the actual kernel is a square (1 + 2 * kernelSize) per side. + * border: When accessing an image, a swath of pixels this wide is ignored. + * footprint: When accessing a stample, a square of pixels, footprint pixels per side, + * are looked at. We must ensure that (footprint+kernelSize) pixels exist + * around the center. + */ + +#include +#include +#include "pslib.h" +#include "psConstants.h" +#include "pmImageSubtract.h" + +/******************************************************************************* + Private alloc/free functions. +XXX: It's not clear if the SubtractionKernels alloc/free functions are feasable. + ******************************************************************************/ +void p_pmStampFree(pmStamp *stamp) +{ + psFree(stamp->matrix); + psFree(stamp->vector); +} + +pmStamp *p_pmStampAlloc(pmStampStatus status) +{ + pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp)); + stamp->x = 0; + stamp->p_xSize = 0; + stamp->y = 0; + stamp->p_ySize = 0; + stamp->matrix = NULL; + stamp->vector = NULL; + stamp->status = status; + + psMemSetDeallocator(stamp, (psFreeFunc)p_pmStampFree); + + return(stamp); +} + +void p_pmSubtractionKernelsFree(psSubtractionKernels *kernels) +{ + psFree(kernels->u); + psFree(kernels->v); + psFree(kernels->sigma); + psFree(kernels->xOrder); + psFree(kernels->yOrder); + psFree(kernels->preCalc); + + psFree(kernels); +} + +psSubtractionKernels *p_pmSubtractionKernelsAlloc(int numBasisFunctions, + pmSubtractionKernelsType type) +{ + psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); + + tmp->type = type; + psMemSetDeallocator(tmp, (psFreeFunc) p_pmSubtractionKernelsFree); + return(tmp); +} + +/******************************************************************************* +psSubtractionKernels struct. + ******************************************************************************/ +psSubtractionKernels *pmSubtractionKernelsAllocPOIS(int size, + int spatialOrder) +{ + psTrace("ImageSubtract.pmSubtractionKernelsAllocPOIS", 3, + "Calling pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder); + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); + // + // Calculate the number of basis functions (nBF) + // + psS32 xKernelHalfSize = size; + psS32 yKernelHalfSize = size; + psS32 nBF = (2 * xKernelHalfSize + 1) * + (2 * yKernelHalfSize + 1) * + (spatialOrder + 1) * + (spatialOrder + 2) / 2; + + // + // Generate the new psSubtractionKernels data structure: + // + psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); + tmp->type = PM_SUBTRACTION_KERNEL_POIS; + tmp->u = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->v = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->u->n = tmp->u->nalloc; + tmp->v->n = tmp->v->nalloc; + tmp->sigma = NULL; + tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->xOrder->n = tmp->xOrder->nalloc; + tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->yOrder->n = tmp->yOrder->nalloc; + tmp->subIndex = 0; + tmp->preCalc = NULL; + tmp->size = size; + tmp->spatialOrder = spatialOrder; + + // + // This corresponds to the Kernel Basis Function with (u, v) = (0, 0) + // + // Put the (u,v) = (0,0) component right at the start of the list + // for convenience. + // + psS32 ptr = 0; + for (psS32 order = 0; order <= spatialOrder; order++) { + for (psS32 xOrder = 0; xOrder <= order; xOrder++) { + psS32 yOrder = order - xOrder; + tmp->u->data.F32[ptr] = 0; + tmp->v->data.F32[ptr] = 0; + tmp->xOrder->data.F32[ptr] = xOrder; + tmp->yOrder->data.F32[ptr] = yOrder; + ptr++; + } + } + + // + // Iterate over (u,v). Generate a set of kernels for each (u, v). + // + for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { + for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { + // Already did (u,v) = (0,0): it's at the start, so skip it now. + if ((u != 0) || (v != 0)) { + // + // Iterate over spatial order. This loop creates the terms for + // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. + // + for (psS32 order = 0; order <= spatialOrder; order++) { + for (psS32 xOrder = 0; xOrder <= order; xOrder++) { + psS32 yOrder = order - xOrder; + tmp->u->data.F32[ptr] = u; + tmp->v->data.F32[ptr] = v; + tmp->xOrder->data.F32[ptr] = xOrder; + tmp->yOrder->data.F32[ptr] = yOrder; + ptr++; + } + } + } + } + } + + psTrace("ImageSubtract.pmSubtractionKernelsAllocPOIS", 3, + "Exiting pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder); + return(tmp); +} + +/******************************************************************************* +XXX: Get the types correct (u, v, xOrder, yOrder). + +XXX: Code review this. + ******************************************************************************/ +psSubtractionKernels *pmSubtractionKernelsAllocISIS(const psVector *sigmas, + const psVector *orders, + int size, + int spatialOrder) +{ + PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + psTrace("ImageSubtract.pmSubtractionKernelsAllocISIS", 3, + "Calling pmSubtractionKernelsAllocISIS(%d, %d, %d, %d)\n", + sigmas->n, orders->n, size, spatialOrder); + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); + PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); + // + // Calculate the number of basis functions (nBF). + // + psS32 numSigmas = sigmas->n; + + // XXX: Get rid of the sigma loop? We merely multiple nBF by numSigmas? + // XXX: Verify that all this is correct. + psS32 nBF = 0; + for (psS32 s = 0 ; s < numSigmas ; s++) { + for (psS32 o = 0 ; o < orders->n ; o++) { + nBF+=((orders->data.S32[o] + 1) * (orders->data.S32[o] + 2) / 2); + } + } + nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2); + + // + // Generate the new psSubtractionKernels data structure: + // + psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); + tmp->type = PM_SUBTRACTION_KERNEL_ISIS; + tmp->u = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->v = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->u->n = tmp->u->nalloc; + tmp->v->n = tmp->v->nalloc; + tmp->sigma = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->subIndex = 0; + tmp->size = size; + tmp->spatialOrder = spatialOrder; + tmp->preCalc = psArrayAlloc(nBF); + tmp->u->n = tmp->u->nalloc; + tmp->v->n = tmp->v->nalloc; + tmp->sigma->n = tmp->sigma->nalloc; + tmp->xOrder->n = tmp->xOrder->nalloc; + tmp->yOrder->n = tmp->yOrder->nalloc; + tmp->preCalc->n = tmp->preCalc->nalloc; + + // + // We loop through all combinations of sigmas and polynomial orders + // creating the kernel basis functions. + // + psS32 ptr = 0; + for (psS32 sigPtr = 0 ; sigPtr < numSigmas ; sigPtr++) { + tmp->sigma->data.F32[sigPtr] = sigmas->data.F32[sigPtr]; + // + // (xOrderP, yOrderP) are the order of the polynomial that modify the + // gaussian in the kernel. They correspond to (j, k) in equation (5) + // from the psModules SDRS. + // + for (psS32 o = 0 ; o < orders->n ; o++) { + for (psS32 orderP = 0 ; orderP <= orders->data.S32[o] ; orderP++) { + for (psS32 xOrderP = 0 ; xOrderP <= orderP ; xOrderP++) { + psS32 yOrderP = orderP - xOrderP; + + psImage *currPreCalc = psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); + PS_IMAGE_SET_F32(currPreCalc, 0.0); + psBool setPreCalc = true; + // + // We loop through all spatial orders. Since they have no effect on + // the preCalc images, we only calculate them once, and store pointers + // in tmp->preCalc->data[ptr] for other spatial orders. + // + for (psS32 order = 0; order <= spatialOrder; order++) { + for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { + PS_ASSERT_INT_LESS_THAN(ptr, nBF, NULL); + + psS32 orderYTerm = order - orderXTerm; + + tmp->u->data.F32[ptr] = xOrderP; + tmp->v->data.F32[ptr] = yOrderP; + tmp->xOrder->data.F32[ptr] = orderXTerm; + tmp->yOrder->data.F32[ptr] = orderYTerm; + tmp->sigma->data.F32[ptr] = sigmas->data.F32[sigPtr]; + tmp->preCalc->data[ptr] = (psPtr *) currPreCalc; + + // + // We calculate the preCalc image only the first time through + // this loop. Otherwise, we increment the memory reference + // counter. + // + if (setPreCalc == true) { + for (psS32 v = -size; v <= size; v++) { + for (psS32 u = -size; u <= size; u++) { + // Scale the (u,v) coordinates in kernel space to [-1.0:1.0]. + psF32 uScaled = ((psF32) u) / ((psF32) size); + psF32 vScaled = ((psF32) v) / ((psF32) size); + + // Compute the value of the kernel at location (u, v): + psF32 exponent = (PS_SQR(uScaled) + PS_SQR(vScaled)) / + (2.0 * PS_SQR(sigmas->data.F32[sigPtr])); + currPreCalc->data.F32[v+size][u+size] = + exp(-exponent) * + pow(uScaled, orderXTerm) * + pow(vScaled, orderYTerm); + } + } + setPreCalc = false; + } else { + psMemIncrRefCounter(currPreCalc); + } + ptr++; + } + } + } + } + } + } + + psTrace("ImageSubtract.pmSubtractionKernelsAllocISIS", 3, + "Exiting pmSubtractionKernelsAllocISIS(%d, %d, %d, %d)\n", + sigmas->n, orders->n, size, spatialOrder); + return(tmp); +} + +/******************************************************************************* +pmSubtractionFindStamps(stamps, image, mask, maskVal, threshold, xNum, yNum, border) + +XXX: The SDRS and the proptotype code differ significantly. + Prototype: When a maximum pixel is found within a stamp, an area of size + 2*footprint is searched around that pixel looking for masked pixels. + SDRS: none of that is required. + +XXX: Do we need to care about case where yNum/xNum does not evenly divide the +nnumber of rows/columns in the image? + ******************************************************************************/ +psArray *pmSubtractionFindStamps(psArray *stamps, ///< Output stamps, or NULL + const psImage *image, ///< Image for which to find stamps + const psImage *mask, ///< Mask + psU32 maskVal, ///< Value for mask + psF32 threshold, ///< Threshold for stamps in the image + psS32 xNum, ///< Number of stamps in x + psS32 yNum, ///< Number of stamps in y + psS32 border ///< Border around image to ignore (should be size of kernel or larger) + ) +{ + psTrace("ImageSubtract.pmSubtractionFindStamps", 3, + "Calling pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n", + maskVal, threshold, xNum, yNum, border); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + if (mask != NULL) { + PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); + } + PS_ASSERT_INT_POSITIVE(xNum, NULL); + PS_ASSERT_INT_POSITIVE(yNum, NULL); + PS_ASSERT_INT_POSITIVE(border, NULL); + PS_ASSERT_INT_LARGER_THAN(image->numCols, xNum, NULL); + PS_ASSERT_INT_LARGER_THAN(image->numRows, yNum, NULL); + PS_ASSERT_INT_LARGER_THAN(image->numCols, (2 * border), NULL); + PS_ASSERT_INT_LARGER_THAN(image->numRows, (2 * border), NULL); + + if (stamps != NULL) { + PS_ASSERT_INT_EQUAL(stamps->n, (xNum * yNum), NULL); + // + // Ensure that a pmStamp struct exists at each psArray location. + // + for (psS32 s = 0 ; s < (xNum * yNum) ; s++) { + if (NULL == stamps->data[s]) { + stamps->data[s] = (psPtr *) p_pmStampAlloc(PM_STAMP_REJECTED); + } + } + } else { + stamps = (psArray *) psArrayAlloc(xNum * yNum); + stamps->n = stamps->nalloc; + for (psS32 s = 0 ; s < (xNum * yNum) ; s++) { + stamps->data[s] = (psPtr *) p_pmStampAlloc(PM_STAMP_REJECTED); + } + } + psS32 numRows = image->numRows; + psS32 numCols = image->numCols; + + // + // Iterate over the image sections + // + // XXX: Must handle cases where image size is not an even multiple of xNum or yNum + // they are currently ignored. + // + psS32 num = 0; + for (psS32 j = 0; j < yNum; j++) { + for (psS32 i = 0; i < xNum; i++) { + pmStamp *stamp = (pmStamp *) stamps->data[num]; + // + // Only find a new stamp if we need to + // + if (stamp->status == PM_STAMP_REJECTED) { + // + // Find maximum non-masked value in the image section, + // but don't include a footprint around the edge + // + psF32 max = -INFINITY; + psS32 bestx = 0; + psS32 besty = 0; + // + // The following nested loop iterates over every pixel in the mask + // associated with this (i, j). It ignores pixels within a + // border of pixels from the image edge. + // + // XXX: verify (numX, numY), then get rid of it. + // + psS32 numX = xNum; + psS32 numY = yNum; + psS32 yMin = border + j * (numCols - 2.0 * border) / numY; + psS32 yMax = (border + (j + 1) * (numCols - 2.0 * border) / numY) - 1; + psS32 xMin = border + i * (numRows - 2.0 * border) / numX; + psS32 xMax = (border + (i + 1) * (numRows - 2.0 * border) / numX) - 1; + + if ((yMax >= image->numRows) || + (xMax >= image->numCols) || + (yMin < 0) || + (xMin < 0)) { + // XXX: We skip this stamp since its borders extends beyond the image. + // XXX: This is here mainly as a safeguard. We need to redefine the above + // min/max pixels calculation to ensure that all stamps are legitimate. + + stamp->x = -1; + stamp->y = -1; + stamp->status = PM_STAMP_NONE; + } else { + stamp->p_xSize = 1 + (xMax - xMin); + stamp->p_ySize = 1 + (yMax - yMin); + stamp->p_xMin = xMin; + stamp->p_xMax = xMax; + stamp->p_yMin = yMin; + stamp->p_yMax = yMax; + + for (psS32 y = yMin; y <= yMax ; y++) { + for (psS32 x = xMin; x <= xMax ; x++) { + // Determine if this pixel is larger than the max, and unmasked. + if (image->data.F32[y][x] > max) { + if ((mask == NULL) || !((mask->data.U8[y][x]) & maskVal)) { + max = image->data.F32[y][x]; + bestx = x; + besty = y; + } + } + } + } + + // + // If the max pixel is larger than the threshold, we keep this stamp. + // Otherwise, mark the stamp as PM_STAMP_NONE + // + if (image->data.F32[besty][bestx] >= threshold) { + stamp->x = bestx; + stamp->y = besty; + stamp->status = PM_STAMP_RECALC; + } else { + stamp->x = bestx; + stamp->y = besty; + stamp->status = PM_STAMP_NONE; + } + } + } + num++; + } + } + psTrace("ImageSubtract.pmSubtractionFindStamps", 3, + "Exiting pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n", + maskVal, threshold, xNum, yNum, border); + return(stamps); +} + +/******************************************************************************* +GenSpatialOrder(spatialOrder, x, y): generates and returns a psImage in which +the [i][j] location is calculated as (x^i * y^j). + +XXX: Modify loop so that terms higher than spatialOrder are not computed. + +XXX: Modify this so that [i][j] location is calculated as (x^j * y^i)? + ******************************************************************************/ +static psImage *GenSpatialOrder(psS32 spatialOrder, + psF32 x, + psF32 y) +{ + psTrace("ImageSubtract.GenSpatialOrder", 4, + "Calling GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y); + + psImage *polyValues = psImageAlloc(spatialOrder+1, spatialOrder+1, PS_TYPE_F64); + + psF64 xSum = 1.0; + psF64 ySum = 1.0; + for (psS32 i = 0; i < spatialOrder + 1; i++) { + ySum = xSum; + for (psS32 j = 0; j < spatialOrder + 1; j++) { + polyValues->data.F64[i][j] = ySum; + ySum*= y; + } + xSum*= x; + } + + psTrace("ImageSubtract.GenSpatialOrder", 4, + "Exiting GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y); + + return(polyValues); +} + + +/******************************************************************************* +IsisKernelConvolve(input, kernels, kernelID, col, row): This routine +convolves a single kernel basis function with a pixel in an image. + ******************************************************************************/ +static psF32 IsisKernelConvolve(const psImage *input, + const psSubtractionKernels *kernels, + psS32 kernelID, + psS32 col, + psS32 row) +{ + + + psTrace("ImageSubtract.IsisKernelConvolve", 4, + "Calling IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row); + psS32 spatialOrder = kernels->spatialOrder; + psS32 kernelSize = kernels->size; + psS32 xOrder = (psS32) kernels->xOrder->data.F32[kernelID]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[kernelID]; + psF32 numColsHalf = 0.5 * (psF32) input->numCols; + psF32 numRowsHalf = 0.5 * (psF32) input->numRows; + psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position + psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position + + psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); + + psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; + + psImage *preCalc = (psImage *) kernels->preCalc->data[kernelID]; + + // XXX: Are the following asserts really necessary? + PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(row-kernelSize, 0, NAN); + PS_ASSERT_INT_LESS_THAN(row+kernelSize, input->numRows, NAN); + PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(col-kernelSize, 0, NAN); + PS_ASSERT_INT_LESS_THAN(col+kernelSize, input->numCols, NAN); + psF32 conv = 0.0; + for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { + for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { + conv += input->data.F32[yy+row][xx+col] * + preCalc->data.F32[yy+kernelSize][xx+kernelSize] * + polyVal; + } + } + psFree(polyValues); + + psTrace("ImageSubtract.IsisKernelConvolve", 4, + "Exiting IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row); + return(conv); +} + +/******************************************************************************* +ConvolvePixelPois(input, mask, badStampMaskVal, solution, kernels, col, row): + +This routine takes a single pixel in the psImage input and convolves it with +the set of kernel basis functions and their appropriate weights in solution. +It returns the value of the convolved pixel. + +XXX: Static structure for polyValues? + ******************************************************************************/ +static psF32 ConvolvePixelPois(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels, + psS32 col, + psS32 row) +{ + psTrace("ImageSubtract.ConvolvePixelPois", 4, + "Calling ConvolvePixelPois(%d, %d)\n", col, row); + psS32 nBF = kernels->u->n; + psF32 numColsHalf = 0.5 * (psF32) input->numCols; + psF32 numRowsHalf = 0.5 * (psF32) input->numRows; + psF32 background = solution->data.F64[solution->n-1]; + psS32 spatialOrder = kernels->spatialOrder; + psF32 conv = background; // Initial convolved value + + if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { + psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position + psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position + psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); + + // Iterate over the kernel basis functions + for (psS32 k = 0; k < nBF; k++) { + psS32 u = (psS32) kernels->u->data.F32[k]; + psS32 v = (psS32) kernels->v->data.F32[k]; + + // XXX: What's the story with this? + #if 0 + + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; + #else + + psF32 polyVal = 1.0; + #endif + + // XXX: Why this? + if (k == 0) { + conv += solution->data.F64[k] * input->data.F32[row - v][col - u] * polyVal; + } else { + conv += solution->data.F64[k] * + (input->data.F32[row - v][col - u] * polyVal - input->data.F32[row][col]); + } + } + psFree(polyValues); + } + + psTrace("ImageSubtract.ConvolvePixelPois", 4, + "Exiting ConvolvePixelPois(%d, %d)\n", col, row); + return(conv); +} + + + +/******************************************************************************* +ConvolvePixelIsis(input, mask, badStampMaskVal, solution, kernels, col, row): + +This routine takes a single pixel in the psImage input and convolves it with +the set of kernel basis functions and their appropriate weights in solution. +It returns the value of the convolved pixel. + +XXX: Static structure for polyValues? + ******************************************************************************/ +static psF32 ConvolvePixelIsis(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels, + psS32 col, + psS32 row) +{ + psTrace("ImageSubtract.ConvolvePixelIsis", 4, + "Calling ConvolvePixelIsis(%d, %d)\n", col, row); + psF32 background = solution->data.F64[solution->n-1]; + psF32 conv = background; // Initial convolved value + + if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { + // Iterate over the kernel basis functions + for (psS32 k = 0; k < kernels->u->n; k++) { + conv += IsisKernelConvolve(input, kernels, k, col, row); + } + } + + psTrace("ImageSubtract.ConvolvePixelIsis", 4, + "Exiting ConvolvePixelIsis(%d, %d)\n", col, row); + return(conv); +} + +/******************************************************************************* +ConvolveImage(input, mask, badStampMaskVal, solution, kernels): convolves an +arbitrary image with either an ISIS or POIS set of kernel basis functions. + ******************************************************************************/ +static psImage *ConvolveImage(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels) +{ + psTrace("ImageSubtract.ConvolveImage", 4, "Calling ConvolveImage()\n"); + PS_ASSERT_IMAGE_NON_NULL(input, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(input, NULL); + PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, NULL); + if (mask != NULL) { + PS_ASSERT_IMAGES_SIZE_EQUAL(input, mask, NULL); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); + } + PS_ASSERT_VECTOR_NON_NULL(solution, NULL); + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); + PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); + if (kernels->preCalc != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); + } else { + if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); + return(NULL); + } + } + psS32 nBF = kernels->u->n; + PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); + + psS32 numCols = input->numCols; + psS32 numRows = input->numRows; + psS32 kernelSize = kernels->size; + + psImage *convolved = psImageAlloc(numCols, numRows, PS_TYPE_F32); + + for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { + for (psS32 x = kernelSize; x < numCols - kernelSize; x++) { + if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + convolved->data.F32[y][x] = ConvolvePixelPois(input, mask, badStampMaskVal, + solution, kernels, x, y); + } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + convolved->data.F32[y][x] = ConvolvePixelIsis(input, mask, badStampMaskVal, + solution, kernels, x, y); + } else { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel type. Returning NULL\n"); + return(NULL); + } + } + } + + // + // Pad the rest of the convolved image with 0.0 + // + for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { + for (psS32 x = 0; x < kernelSize; x++) { + convolved->data.F32[y][x] = 0.0; + } + for (psS32 x = numCols - kernelSize; x < numCols; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + for (psS32 y = 0; y < kernelSize; y++) { + for (psS32 x = 0; x < numCols; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + for (psS32 y = numRows - kernelSize; y < numRows; y++) { + for (psS32 x = 0; x < numCols; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + + psTrace("ImageSubtract.ConvolveImage", 4, "Exiting ConvolveImage()\n"); + return convolved; +} + + + + + +/******************************************************************************* +XXX: We should assert that the (footprint, kernelSize, imageSize) stuff +ensures that all data is accessed in bounds? + ******************************************************************************/ +bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, + const psImage *reference, ///< Reference image + const psImage *input, ///< Input image + const psSubtractionKernels *kernels, ///< The kernel basis functions + psS32 footprint ///< Half-size of region over which to calculate equation + ) +{ + psTrace("ImageSubtract.pmSubtractionCalculateEquation", 3, + "Calling pmSubtractionCalculateEquation()\n"); + PS_ASSERT_PTR_NON_NULL(stamps, false); + PS_ASSERT_IMAGE_NON_NULL(reference, false); + PS_ASSERT_IMAGE_NON_EMPTY(reference, false); + PS_ASSERT_IMAGE_NON_NULL(input, false); + PS_ASSERT_IMAGE_NON_EMPTY(input, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(reference, input, false); + PS_ASSERT_PTR_NON_NULL(kernels, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); + if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); + } + psS32 kernelSize = kernels->size; + PS_ASSERT_INT_NONNEGATIVE(footprint, false); + // + // For each legitimate stamp, ensure that the footprint is small enough to perform + // the full calculation. + // + // XXX: Verify with IfA that this is a reasonable action. + // + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + if (stamp->status == PM_STAMP_RECALC) { + // XXX: trace message + // printf("stamp %d (x, y) is (%d, %d). Footprint is %d. kernelSize is %d.\n", s, stamp->x, stamp->y, footprint, kernelSize); + if (((stamp->y - (footprint + kernelSize)) < 0) || + ((stamp->x - (footprint + kernelSize)) < 0) || + ((stamp->y + footprint + kernelSize) >= input->numRows) || + ((stamp->x + footprint + kernelSize) >= input->numCols)) { + stamp->status = PM_STAMP_NONE; + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: stamp %d will be ignored. It exceeds image size: access columns (%d to %d) and rows (%d to %d)\n", + s, + stamp->x - (footprint + kernelSize), + (stamp->x + footprint + kernelSize) - 1, + stamp->y - (footprint + kernelSize), + (stamp->y + footprint + kernelSize) - 1); + } + } + } + + psS32 numHalfRows = reference->numRows; + psS32 numHalfCols = reference->numCols; + psS32 spatialOrder = kernels->spatialOrder; + + // + // The numSolveParams incorporates the additional parameter for the + // background value, which we must solve for. + // + psS32 numKernels = kernels->u->n; + int numSolveParams = numKernels + 1; + int bgIndex = numKernels; // Index in matrix for the background + + // + // We iterate over each stamp, allocate the matrix and vectors if + // necessary, and then calculate those matrix/vectors. + // + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + psTrace("pmSubtractionCalculateEquation", 6, "subCalcEqn(): stamp %d\n", s); + if (stamp->status == PM_STAMP_RECALC) { + psTrace("pmSubtractionCalculateEquation", 6, "subCalcEqn(): stamp %d: status is PM_STAMP_RECALC.\n", s); + psImage *stampMatrix = stamp->matrix; + psVector *stampVector = stamp->vector; + + if (stampMatrix == NULL) { + stampMatrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); + stamp->matrix = stampMatrix; + } else { + PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_SIZE(stampMatrix, numSolveParams, numSolveParams, false); + } + PS_IMAGE_SET_F64(stampMatrix, 0.0); + + if (stampVector == NULL) { + stampVector = psVectorAlloc(numSolveParams, PS_TYPE_F64); + stampVector->n = stampVector->nalloc; + stamp->vector = stampVector; + } else { + PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, false); + PS_ASSERT_VECTOR_SIZE(stampVector, numSolveParams, false); + } + PS_VECTOR_SET_F64(stampVector, 0.0); + psTrace("pmSubtractionCalculateEquation", 6, "subCalcEqn(): stamp %d: allocate matrix and vector.\n", s); + + // + // Evaluate the spatial-order polynomial. The [i][j]-th element of + // the psImage polyValues will hold (x^i * y^j) for the stamp. The + // (x, y) value are scaled to [-1:1] + // + psImage *polyValues = GenSpatialOrder(spatialOrder, + ((psF64) (stamp->x - numHalfCols)) / ((psF64) numHalfCols), + ((psF64) (stamp->y - numHalfRows)) / ((psF64) numHalfRows)); + + psTrace("pmSubtractionCalculateEquation", 6, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s); + + if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + // + // Iterate over all pixels surrounding this stamp. + // + for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { + for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { + psTrace("pmSubtractionCalculateEquation", 6, "subCalcEqn(): pixel (%d, %d).\n", y, x); + + // The inverse of the noise, squared. + psF32 invNoise2 = 1.0/reference->data.F32[y][x]; + + // + // Iterate over the first convolution */ + // + for (psS32 k1 = 0; k1 < numKernels; k1++) { + psS32 u1 = kernels->u->data.F32[k1]; // Offset in x + psS32 v1 = kernels->v->data.F32[k1]; // Offset in y + psS32 i1 = kernels->xOrder->data.F32[k1]; // Polynomial order in x + psS32 j1 = kernels->yOrder->data.F32[k1]; // Polynomial order in y + + // + // First convolution. This will set the value for the stampVector. + // + // XXX: verify the [y-v2][x-u2] subscript. This generated errors in + // testing, depending on kernel size and footprint. + // + psF32 conv1 = polyValues->data.F64[j1][i1] * reference->data.F32[y - v1][x - u1]; + + // + // Assuming that the first kernel component is 0 order in x and y, and 0 offset + // XXX: I don't understand this: + // + if (k1 != 0) { + conv1 -= reference->data.F32[y][x]; + } + + // + // Iterate over the second convolution + // + for (psS32 k2 = k1; k2 < numKernels; k2++) { + psS32 u2 = (psS32) kernels->u->data.F32[k2]; // Offset in x + psS32 v2 = (psS32) kernels->v->data.F32[k2]; // Offset in y + psS32 i2 = (psS32) kernels->xOrder->data.F32[k2]; // Polynomial order in x + psS32 j2 = (psS32) kernels->yOrder->data.F32[k2]; // Polynomial order in y + // + // XXX: verify the [y-v2][x-u2] subscript. This generated errors in + // testing, depending on kernel size and footprint. + // + // Second convolution + // + psF32 conv2 = polyValues->data.F64[j2][i2] * + reference->data.F32[y-v2][x-u2]; + // + // Assuming that the first kernel component is 0 order in x and y, and 0 offset + // XXX: I don't understand this: + if (k2 != 0) { + // + conv2 -= reference->data.F32[y][x]; + } + + // Add into the matrix element + stampMatrix->data.F64[k1][k2] += conv1 * conv2 * invNoise2; + + } // Iteration on second convolution + + // Add into the vector element + stampVector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2; + + /* Background term */ + stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2; + + } // Iteration on first convolution + + // + // Background only terms. + // XXX: understand this. + // + stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; + stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; + } + } + } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + // + // Iterate over all pixels surrounding this stamp. + // + // XXX: Why isn't there a polyValues term here? + // + + for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { + for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { + psTrace("pmSubtractionCalculateEquation", 6, "subCalcEqn(): pixel (%d, %d).\n", y, x); + psF32 invNoise2 = 1.0/reference->data.F32[y][x]; // The inverse of the noise, squared. + + for (psS32 k1 = 0; k1 < numKernels; k1++) { + psF32 conv1 = IsisKernelConvolve(reference, kernels, k1, x, y); + + for (psS32 k2 = k1; k2 < numKernels; k2++) { + //printf("(k1, k2) is (%d, %d)\n", k1, k2); + psF32 conv2 = IsisKernelConvolve(reference, kernels, k2, x, y); + stampMatrix->data.F64[k1][k2]+= conv1 * conv2 * invNoise2; + } + stampVector->data.F64[k1]+= input->data.F32[y][x] * conv1 * invNoise2; + stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2; + } + stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; + stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; + } + } + } else { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel->type.\n"); + return(false); + } + psFree(polyValues); + + // XXX: Generate psTrace() + if (0) { + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + if (stamp->status == PM_STAMP_RECALC) { + psVector *stampVector = stamp->vector; + printf("STAMP: stamp %d vector:\n", s); + PS_VECTOR_PRINT_F64(stampVector); + } + } + } + + // + // Fill in other side of symmetric matrix + // + // XXX: understand this. + // XXX: Why aren't they using numSolveParams instead of numKernels? + // XXX: is this POIS specific? + // + for (psS32 k1 = 0; k1 < numKernels; k1++) { + for (psS32 k2 = 0; k2 < k1; k2++) { + stampMatrix->data.F64[k1][k2] = stampMatrix->data.F64[k2][k1]; + } + stampMatrix->data.F64[bgIndex][k1] = stampMatrix->data.F64[k1][bgIndex]; + } + + // + // XXX: Why aren't they using numSolveParams instead of numKernels? + // XXX: is this POIS specific? + // + #define XXX_CONFIG_PENALTY 1.0 + for (psS32 k = 0; k < numKernels; k++) + { + psS32 u = kernels->u->data.F32[k]; // Offset in x + psS32 v = kernels->v->data.F32[k]; // Offset in y + stampMatrix->data.F64[k][k] += XXX_CONFIG_PENALTY * (psF32)(u*u + v*v); + } + stamp->status = PM_STAMP_USED; + } else { + // Stamp is ignored since it's not PM_STAMP_RECALC + } + } + psTrace("ImageSubtract.pmSubtractionCalculateEquation", 3, + "Exiting pmSubtractionCalculateEquation()\n"); + return(true); +} + + + + +/******************************************************************************* + ******************************************************************************/ +psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL + const psArray *stamps ///< Array of stamps + ) +{ + psTrace("ImageSubtract.pmSubtractionSolveEquation", 3, + "Calling pmSubtractionSolveEquation()\n"); + PS_ASSERT_PTR_NON_NULL(stamps, NULL); + psS32 size = -1; + psS32 s = 0; + + // + // Determine the size of the stamp vectors and matrix. + // We iterate until we find the first acceptable stamp. + // + while ((size == -1) && (s < stamps->n)) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + PS_ASSERT_PTR_NON_NULL(stamp, NULL); + if (stamp->status == PM_STAMP_USED) { + size = ((pmStamp *) stamps->data[s])->vector->n; + PS_ASSERT_INT_POSITIVE(size, NULL); + } + s++; + } + if (size == -1) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: no acceptable stamps. Returning NULL\n"); + return(NULL); + } + + if (solution != NULL) { + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_SIZE(solution, size, NULL); + } else { + solution = psVectorAlloc(size, PS_TYPE_F64); + solution->n = solution->nalloc; + } + + // + // Create the solution matrix and vector. + // + // XXX: Test these functions with size=-1. This caused seg faults during test. + // This should be done in the psImage.c and psVector.c test files. It + // should never occur here. + // + psImage *sumMatrix = psImageAlloc(size, size, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(size, PS_TYPE_F64); + sumVector->n = sumVector->nalloc; + PS_VECTOR_SET_F64(sumVector, 0.0); + PS_IMAGE_SET_F64(sumMatrix, 0.0); + + // + // Verify that all stamps have similar sizes. + // Compute the sum matrix and vector. + // + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + + if (stamp->status == PM_STAMP_USED) { + PS_ASSERT_INT_EQUAL(((pmStamp *) stamps->data[s])->vector->n, size, NULL); + + psImage *stampMatrix = stamp->matrix; + psVector *stampVector = stamp->vector; + PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_SIZE(stampVector, size, NULL); + PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, NULL); + PS_ASSERT_IMAGE_SIZE(stampMatrix, size, size, NULL); + + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stampMatrix); + (void)psBinaryOp(sumVector, sumVector, "+", stampVector); + } + } + psVector *permutation = NULL; + // XXX: Check output from these routines. + + // XXX: psTrace() + if (0) { + PS_IMAGE_PRINT_F64(sumMatrix); + } + + psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); + if (luMatrix == NULL) { + psError(PS_ERR_UNKNOWN, true, "Failed to LU-Decompose the matrix.\n"); + psFree(sumMatrix); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return(NULL); + } + // XXX: psTrace() + if (0) { + PS_IMAGE_PRINT_F64(luMatrix); + } + + solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation); + // XXX: psTrace() + // XXX: should we be checking for NAN's in the solution vector? + if (0) { + PS_VECTOR_PRINT_F64(solution); + } + if (solution == NULL) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the matrix.\n"); + psFree(sumMatrix); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return(NULL); + } + + psFree(sumMatrix); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + + psTrace("ImageSubtract.pmSubtractionSolveEquation", 3, + "Exiting pmSubtractionSolveEquation()\n"); + return(solution); +} + + +/******************************************************************************* + ******************************************************************************/ +static psVector *CalculateDeviations(psVector *deviations, + psArray *stamps, + psS32 footprint, + const psImage *refImage, + const psImage *inImage, + const psImage *mask, + psU32 badStampMaskVal, + const psSubtractionKernels *kernels, + const psVector *solution) +{ + psTrace("ImageSubtract.CalculateDeviations", 4, + "Calling CalculateDeviations()\n"); + PS_ASSERT_PTR_NON_NULL(stamps, NULL); + if (deviations != NULL) { + PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(deviations, stamps, NULL); + } else { + deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); + deviations->n = deviations->nalloc; + // XXX: Probably not necessary. + PS_VECTOR_SET_F32(deviations, 0.0); + } + PS_ASSERT_IMAGE_NON_NULL(refImage, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(refImage, NULL); + PS_ASSERT_IMAGE_TYPE(refImage, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGE_NON_NULL(inImage, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(inImage, NULL); + PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, NULL); + PS_ASSERT_IMAGE_NON_NULL(mask, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(mask, NULL); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, NULL); + PS_ASSERT_VECTOR_NON_NULL(solution, NULL); + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); + PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); + if (kernels->preCalc != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); + } else { + if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); + return(NULL); + } + } + psS32 nBF = kernels->u->n; + PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); + + psS32 kernelSize = kernels->size; + int xSize = footprint + kernelSize; + int ySize = footprint + kernelSize; + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics + psImage *subStamp = psImageAlloc(2 * xSize, 2 * ySize, PS_TYPE_F32); // Subtraction of stamp + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest + psS32 x = stamp->x; // Stamp x coord + psS32 y = stamp->y; // Stamp y coord + if (stamp->status == PM_STAMP_USED) { + + psRegion myReg = psRegionSet(x - xSize, x + xSize, y - ySize, y + ySize); + psImage *refStamp = psImageSubset((psImage *) refImage, myReg); + psImage *inStamp = psImageSubset((psImage *) inImage, myReg); + psImage *maskStamp = psImageSubset((psImage *) mask, myReg); + psImage *convRefStamp = ConvolveImage(refStamp, maskStamp, badStampMaskVal, solution, kernels); + + // Calculate chi^2 + (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp); + (void)psBinaryOp(subStamp, subStamp, "/", inStamp); + (void)psBinaryOp(subStamp, subStamp, "*", subStamp); + myReg = psRegionSet(kernelSize, kernelSize + 2 * footprint, + kernelSize, kernelSize + 2 * footprint); + psImage *subStampTrim = psImageSubset((psImage *) subStamp, myReg); + psImage *maskStampTrim = psImageSubset((psImage *) maskStamp, myReg); + psImageStats(stats, subStampTrim, maskStampTrim, badStampMaskVal); + + deviations->data.F32[s] = stats->sampleMean * (psF32)footprint * (psF32)footprint * 4.0; + // XXX: Allocate and free these elsewhere. + psFree(refStamp); + psFree(inStamp); + psFree(maskStamp); + psFree(convRefStamp); + psFree(subStampTrim); + psFree(maskStampTrim); + } + } + + psFree(stats); + psFree(subStamp); + + psTrace("ImageSubtract.CalculateDeviations", 4, + "Exiting CalculateDeviations()\n"); + return deviations; +} + +/******************************************************************************* + ******************************************************************************/ +bool pmSubtractionRejectStamps(psArray *stamps, ///< Array of stamps to check for rejection + psImage *mask, ///< Mask image + psU32 badStampMaskVal, ///< Value to use in mask for bad stamp + psS32 footprint, ///< Region to mask if stamp is bad + psF32 sigmaRej, ///< Number of RMS deviations above zero at which to reject + const psImage *refImage, ///< Reference image + const psImage *inImage, ///< Input image + const psVector *solution, ///< Solution vector + const psSubtractionKernels *kernels ///< Array of kernel parameters + ) +{ + psTrace("ImageSubtract.pmSubtractionRejectStamps", 3, + "Calling pmSubtractionRejectStamps()\n"); + PS_ASSERT_PTR_NON_NULL(stamps, false); + PS_ASSERT_IMAGE_NON_NULL(refImage, false); + PS_ASSERT_IMAGE_NON_EMPTY(refImage, false); + PS_ASSERT_IMAGE_TYPE(refImage, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(inImage, false); + PS_ASSERT_IMAGE_NON_EMPTY(inImage, false); + PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, false); + PS_ASSERT_IMAGE_NON_NULL(mask, false); + PS_ASSERT_IMAGE_NON_EMPTY(mask, false); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, false); + PS_ASSERT_VECTOR_NON_NULL(solution, false); + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false); + PS_ASSERT_PTR_NON_NULL(kernels, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); + if (kernels->preCalc != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); + } else { + if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); + return(false); + } + } + + psS32 nBF = kernels->u->n; + PS_ASSERT_VECTOR_SIZE(solution, nBF+1, false); + + psVector *deviations = CalculateDeviations(NULL, + stamps, + footprint, + refImage, + inImage, + mask, + badStampMaskVal, + kernels, + solution); + // + // Calculate the deviation from zero. + // + psF64 meanDev = 0.0; + psS32 numDev = 0; + for (psS32 i = 0; i < deviations->n; i++) { + pmStamp *stamp = stamps->data[i]; + if (stamp->status == PM_STAMP_USED) { + meanDev += PS_SQR(deviations->data.F32[i]); + numDev++; + } + } + psF32 rmsDev = sqrtf(meanDev / (psF64)(numDev - 1)); + psF32 limit = rmsDev * sigmaRej; + + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + if (stamp->status == PM_STAMP_USED && fabsf(deviations->data.F32[s]) > limit) { + // Mask out the stamp in the image so you don't find it again + for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { + for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { + mask->data.U8[y][x] |= badStampMaskVal; + } + } + + // Set stamp for replacement + stamp->x = 0; + stamp->y = 0; + stamp->status = PM_STAMP_REJECTED; + } + } + + psFree(deviations); + psTrace("ImageSubtract.pmSubtractionRejectStamps", 3, + "Exiting pmSubtractionRejectStamps()\n"); + return(true); +} + +/******************************************************************************* + ******************************************************************************/ +psImage *pmSubtractionKernelImage(psImage *out, + const psVector *solution, + const psSubtractionKernels *kernels, + psF32 x, + psF32 y + ) +{ + psTrace("ImageSubtract.pmSubtractionKernelImage", 3, + "Calling pmSubtractionKernelImage()\n"); + PS_ASSERT_VECTOR_NON_NULL(solution, NULL); + PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); + if (kernels->preCalc != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); + } else { + if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); + return(NULL); + } + } + PS_ASSERT_INT_EQUAL(1+kernels->u->n, solution->n, NULL); + + psS32 nBF = kernels->u->n; + psS32 spatialOrder = kernels->spatialOrder; + psS32 kernelSize = kernels->size; + + if (out != NULL) { + if ((out->numCols < (1+2*kernelSize)) || (out->numRows < (1+2*kernelSize))) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: out image is not large enough.\n"); + return(out); + } + } else { + out = psImageAlloc(1+2*kernelSize, 1+2*kernelSize, PS_TYPE_F32); + } + PS_IMAGE_SET_F32(out, 0.0); + + // + // Generate the spatial-order polynomial. The [i][j]-th element of + // the psImage polyValues will hold (x^i * y^j) for the stamp. + // + psImage *polyValues = GenSpatialOrder(spatialOrder, x, y); + + // XXX: switch (i, j) so they correspond to (x, y). + if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + for (psS32 k = 0 ; k < nBF ; k++) { + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; + + // XXX: Verify that this is correct. + for (psS32 i = -kernelSize ; i <= kernelSize ; i++) { + for (psS32 j = -kernelSize ; j <= kernelSize ; j++) { + psImage *preCalc = (psImage *) kernels->preCalc->data[k]; + out->data.F32[i+kernelSize][j+kernelSize]+= + solution->data.F64[k] * + preCalc->data.F32[i+kernelSize][j+kernelSize] * + polyVal; + } + } + } + } else if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + for (psS32 k = 0 ; k < nBF ; k++) { + // XXX: Why don't we have compilation warnings on type here (if + // we remove the (psS32) cast)? + psS32 u = (psS32) kernels->u->data.F32[k]; + psS32 v = (psS32) kernels->v->data.F32[k]; + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + // XXX: Verify that this is correct. + + out->data.F32[kernelSize - v][kernelSize - u]+= + solution->data.F64[k] * polyValues->data.F64[yOrder][xOrder]; + } + } + psFree(polyValues); + + psTrace("ImageSubtract.pmSubtractionKernelImage", 3, + "Exiting pmSubtractionKernelImage()\n"); + return(out); +} Index: /trunk/psModules/src/imcombine/pmImageSubtract.h =================================================================== --- /trunk/psModules/src/imcombine/pmImageSubtract.h (revision 7018) +++ /trunk/psModules/src/imcombine/pmImageSubtract.h (revision 7018) @@ -0,0 +1,127 @@ +/** @file ImageSubtract.h + * + * This file will contain code which creates a set of kernel basis + * functions, solves for their solution, and applies them to an image. + * + * @author Paul Price, IfA (original prototype) + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-01 01:56:48 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if !defined(PM_IMAGE_COMBINE_H) +#define PM_IMAGE_COMBINE_H + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "psConstants.h" + +typedef enum { + PM_SUBTRACTION_KERNEL_POIS, ///< POIS kernel --- delta functions + PM_SUBTRACTION_KERNEL_ISIS ///< ISIS kernel --- gaussians modified by polynomials +} pmSubtractionKernelsType; + +typedef struct +{ + pmSubtractionKernelsType type; ///< Type ofKernels --- allowing the use of multiple kernels + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) + psVector *sigma; ///< Width of Gaussian (for ISIS) + psVector *xOrder, *yOrder; ///< Spatial Polynomial order (for all) + int subIndex; ///< Index of kernel to be subtracted to maintain flux conservation + psArray *preCalc; ///< Array of images containing pre-calculated kernel (to + ///< accelerate ISIS; don't use for POIS) + psS32 size; ///< The halfsize of the kernel + psS32 spatialOrder; ///< The spatial order of the kernels +} +psSubtractionKernels; + +psSubtractionKernels *pmSubtractionKernelsAllocPOIS( + int size, + int SpatialOrder +); + +psSubtractionKernels *pmSubtractionKernelsAllocISIS( + const psVector *sigmas, + const psVector *orders, + int size, + int SpatialOrder +); + +typedef enum { + PM_STAMP_INIT, ///< Initial state + PM_STAMP_USED, ///< Use this stamp + PM_STAMP_REJECTED, ///< This stamp has been rejected + PM_STAMP_RECALC, ///< Having been reset, this stamp is to be recalculated + PM_STAMP_NONE ///< No stamp in this region +} pmStampStatus; + +typedef struct +{ + int x, y; ///< Position + int p_xSize; + int p_ySize; + int p_xMin; + int p_xMax; + int p_yMin; + int p_yMax; + psImage *matrix; ///< Associated matrix + psVector *vector; ///< Assoicated vector + pmStampStatus status; ///< Status ofstamp +} +pmStamp; + +psArray *pmSubtractionFindStamps( + psArray *stamps, ///< Output stamps, or NULL + const psImage *image, ///< Image for which to find stamps + const psImage *mask, ///< Mask + psU32 maskVal, ///< Value for mask + psF32 threshold, ///< Threshold for stamps in the image + psS32 xNum, ///< Number of stamps in x + psS32 yNum, ///< Number of stamps in y + psS32 border ///< Border around image to ignore (should be size of kernel) +); + +bool pmSubtractionCalculateEquation( + psArray *stamps, ///< The stamps for which to calculate the equation, + const psImage *reference, ///< Reference image + const psImage *input, ///< Input image + const psSubtractionKernels *kernels,///< The kernel basis functions + psS32 footprint ///< Half-size of region over which to calculate equation +); + + +psVector *pmSubtractionSolveEquation( + psVector *solution, ///< Solution vector, or NULL + const psArray *stamps ///< Array of stamps +); + +bool pmSubtractionRejectStamps( + psArray *stamps, ///< Array of stamps to check for rejection + psImage *mask, ///< Mask image + psU32 badStampMaskVal, ///< Value to use in mask for bad stamp + psS32 footprint, ///< Region to mask if stamp is bad + psF32 sigmaRej, ///< Number of RMS deviations above zero at which to reject + const psImage *refImage, ///< Reference image + const psImage *inImage, ///< Input image + const psVector *solution, ///< Solution vector + const psSubtractionKernels *kernels ///< Array of kernel parameters +); + +psImage *pmSubtractionKernelImage( + psImage *out, + const psVector *solution, + const psSubtractionKernels *kernels, + psF32 x, + psF32 y +); + +#endif