Index: /branches/eam_branches/ipp-20101205/psModules/src/extras/pmVisual.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/extras/pmVisual.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/extras/pmVisual.c (revision 30322) @@ -16,34 +16,37 @@ #include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPAfile.h" -#include "pmAstrometryObjects.h" -#include "pmSubtractionStamps.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmFPAExtent.h" - -#include "pmAstrometryVisual.h" -#include "pmSubtractionVisual.h" -#include "pmStackVisual.h" -#include "pmSourceVisual.h" +bool pmSubtractionVisualClose(void); +bool pmAstromVisualClose(void); +bool pmSubtractionVisualClose(void); +bool pmStackVisualClose(void); +bool pmSourceVisualClose(void); + +// #include "pmHDU.h" +// #include "pmFPA.h" +// #include "pmFPAfile.h" +// #include "pmAstrometryObjects.h" +// #include "pmSubtractionStamps.h" +// #include "pmTrend2D.h" +// #include "pmResiduals.h" +// #include "pmGrowthCurve.h" +// #include "pmSpan.h" +// #include "pmFootprintSpans.h" +// #include "pmFootprint.h" +// #include "pmPeaks.h" +// #include "pmMoments.h" +// #include "pmModelFuncs.h" +// #include "pmModel.h" +// #include "pmSourceMasks.h" +// #include "pmSourceExtendedPars.h" +// #include "pmSourceDiffStats.h" +// #include "pmSource.h" +// #include "pmSourceFitModel.h" +// #include "pmPSF.h" +// #include "pmPSFtry.h" +// #include "pmFPAExtent.h" +// #include "pmAstrometryVisual.h" +// #include "pmSubtractionVisual.h" +// #include "pmStackVisual.h" +// #include "pmSourceVisual.h" # if (HAVE_KAPA) Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/Makefile.am (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/Makefile.am (revision 30322) @@ -30,4 +30,5 @@ pmStackReject.h \ pmSubtraction.h \ + pmSubtractionTypes.h \ pmSubtractionAnalysis.h \ pmSubtractionEquation.h \ Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmStackReject.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmStackReject.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmStackReject.c (revision 30322) @@ -7,4 +7,6 @@ #include +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionThreads.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c (revision 30322) @@ -20,4 +20,5 @@ #include "pmHDU.h" // Required for pmFPA.h #include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h (revision 30322) @@ -14,10 +14,9 @@ #define PM_SUBTRACTION_H -#include - -#include -#include -#include -#include +// #include +// #include +// #include +// #include +// #include // if we use the original ppSub implementation, we subtract a central delta-function for all @@ -30,28 +29,4 @@ /// @addtogroup imcombine Image Combinations /// @{ - -/// Mask values for the subtraction mask -typedef enum { - PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking - PM_SUBTRACTION_MASK_BAD_1 = 0x01, // Image 1 is bad - PM_SUBTRACTION_MASK_BAD_2 = 0x02, // Image 2 is bad - PM_SUBTRACTION_MASK_CONVOLVE_1 = 0x04, // If image 1 is convolved, would be poor or bad - PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be poor or bad - PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad - PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad - PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border - PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected -} pmSubtractionMasks; - -typedef struct { - double score; - pmSubtractionMode mode; - int spatialOrder; - int nGood; - psVector *fluxes; - psVector *chisq; - psVector *moments; - psVector *stampMask; -} pmSubtractionQuality; /// Number of terms in a polynomial Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 30322) @@ -7,4 +7,5 @@ #include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionDeconvolve.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 30322) @@ -10,4 +10,5 @@ #include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtractionKernels.h" #include "pmSubtractionDeconvolve.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (revision 30322) @@ -9,4 +9,6 @@ #include "pmErrorCodes.h" #include "pmVisual.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h (revision 30322) @@ -5,12 +5,4 @@ #include "pmSubtractionKernels.h" #include "pmSubtraction.h" - -// typedef enum { -// PM_SUBTRACTION_EQUATION_NONE = 0x00, -// PM_SUBTRACTION_EQUATION_NORM = 0x01, -// PM_SUBTRACTION_EQUATION_BG = 0x02, -// PM_SUBTRACTION_EQUATION_KERNELS = 0x04, -// PM_SUBTRACTION_EQUATION_ALL = 0x07, // value should be NORM | BG | KERNELS -// } pmSubtractionEquationCalculationMode; /// Execute a thread job to calculate the least-squares equation for a stamp Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.v0.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 30322) @@ -9,4 +9,5 @@ #include "pmErrorCodes.h" #include "pmSubtraction.h" +#include "pmSubtractionTypes.h" #include "pmSubtractionKernels.h" #include "pmSubtractionStamps.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionHermitian.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionHermitian.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionHermitian.c (revision 30322) @@ -8,4 +8,5 @@ #include +#include "pmSubtractionTypes.h" #include "pmSubtractionHermitian.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionIO.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionIO.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionIO.c (revision 30322) @@ -12,4 +12,5 @@ #include "pmConceptsRead.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.c (revision 30322) @@ -8,4 +8,6 @@ #include +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.h (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.h (revision 30322) @@ -2,73 +2,6 @@ #define PM_SUBTRACTION_KERNELS_H -#include -#include - -/// Type of subtraction kernel -typedef enum { - PM_SUBTRACTION_KERNEL_NONE, ///< Nothing --- an error - PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions - PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials - PM_SUBTRACTION_KERNEL_ISIS_RADIAL, ///< ISIS + higher-order radial Hermitians - PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels - PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved Hermitian polynomial kernels - PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions - PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction - PM_SUBTRACTION_KERNEL_GUNK, ///< Grid United with Normal Kernel --- POIS and ISIS hybrid - PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction -} pmSubtractionKernelsType; - -/// Modes --- specifies which image to convolve -typedef enum { - PM_SUBTRACTION_MODE_ERR, // Error in the mode - PM_SUBTRACTION_MODE_1, // Convolve image 1 - PM_SUBTRACTION_MODE_2, // Convolve image 2 - PM_SUBTRACTION_MODE_UNSURE, // Not sure yet which image to convolve so try to satisfy both - PM_SUBTRACTION_MODE_DUAL, // Dual convolution -} pmSubtractionMode; - -/// Kernels specification -typedef struct { - pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels - psString description; ///< Description of the kernel parameters - int xMin, xMax, yMin, yMax; ///< Bounds of image (for normalisation) - long num; ///< Number of kernel components (not including the spatial ones) - psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM) - psVector *widths; ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM) - psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM) - float penalty; ///< Penalty for wideness - psVector *penalties1; ///< Penalty for each kernel component - psVector *penalties2; ///< Penalty for each kernel component - bool havePenalties; ///< flag to test if we have already calculated the penalties or not. - int size; ///< The half-size of the kernel - int inner; ///< The size of an inner region - int spatialOrder; ///< The spatial order of the kernels - int bgOrder; ///< The order for the background fitting - pmSubtractionMode mode; ///< Mode for subtraction - psVector *solution1, *solution2; ///< Solution for the PSF matching - psVector *solution1err, *solution2err; ///< error in solution for the PSF matching - // Quality information - float mean, rms; ///< Mean and RMS of chi^2 from stamps - int numStamps; ///< Number of good stamps - float fResSigmaMean; ///< mean fractional stdev of residuals - float fResSigmaStdev; ///< stdev of fractional stdev of residuals - float fResOuterMean; ///< mean fractional positive swing in residuals - float fResOuterStdev; ///< stdev of fractional positive swing in residuals - float fResTotalMean; ///< mean fractional negative swing in residuals - float fResTotalStdev; ///< stdev of fractional negative swing in residuals - psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations -} pmSubtractionKernels; - -// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures -typedef struct { - psVector *uCoords; // used by RINGS - psVector *vCoords; // used by RINGS - psVector *poly; // used by RINGS - - psVector *xKernel; // used by ISIS, HERM, DECONV_HERM - psVector *yKernel; // used by ISIS, HERM, DECONV_HERM - psKernel *kernel; // used by ISIS, HERM, DECONV_HERM -} pmSubtractionKernelPreCalc; +// #include +// #include // Assertion to check pmSubtractionKernels Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMask.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMask.c (revision 30322) @@ -7,4 +7,6 @@ #include "pmErrorCodes.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c (revision 30322) @@ -11,9 +11,10 @@ #include "pmFPA.h" #include "pmHDUUtils.h" +#include "pmSubtractionTypes.h" +#include "pmSubtraction.h" #include "pmSubtractionParams.h" #include "pmSubtractionKernels.h" #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" -#include "pmSubtraction.h" #include "pmSubtractionAnalysis.h" #include "pmSubtractionMask.h" @@ -57,53 +58,4 @@ } - -static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read - const pmReadout *ro1, // Readout 1 - const pmReadout *ro2, // Readout 2 - const psImage *subMask, // Mask for subtraction, or NULL - psImage *variance, // Variance map - const psRegion *region, // Region of interest - float thresh1, // Threshold for stamp finding on readout 1 - float thresh2, // Threshold for stamp finding on readout 2 - float stampSpacing, // Spacing between stamps - float normFrac, // Fraction of flux in window for normalisation window - float sysError, // Relative systematic error in images - float skyError, // Relative systematic error in images - int size, // Kernel half-size - int footprint, // Convolution footprint for stamps - pmSubtractionMode mode // Mode for subtraction - ) -{ - PS_ASSERT_PTR_NON_NULL(stamps, false); - PM_ASSERT_READOUT_NON_NULL(ro1, false); - PM_ASSERT_READOUT_NON_NULL(ro2, false); - PS_ASSERT_IMAGE_NON_NULL(subMask, false); - PS_ASSERT_IMAGE_NON_NULL(variance, false); - PS_ASSERT_PTR_NON_NULL(region, false); - - psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - - psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest - - *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, - size, footprint, stampSpacing, normFrac, sysError, skyError, mode); - if (!*stamps) { - psError(psErrorCodeLast(), false, "Unable to find stamps."); - return false; - } - - memCheck(" find stamps"); - - psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) { - psError(psErrorCodeLast(), false, "Unable to extract stamps."); - return false; - } - - memCheck(" extract stamps"); - pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1); - return true; -} - // Check input arguments static bool subtractionMatchCheck(pmReadout *conv1, pmReadout *conv2, // Convolved images @@ -121,5 +73,5 @@ float badFrac, // Maximum fraction of bad input pixels to accept pmSubtractionMode subMode // Mode of subtraction - ) + ) { if (subMode != PM_SUBTRACTION_MODE_2) { @@ -700,5 +652,5 @@ regionString = psRegionToString(*region); psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n", - regionString, numCols, numRows); + regionString, numCols, numRows); if (stampsName && strlen(stampsName) > 0) { @@ -712,74 +664,90 @@ } - // We get the stamps here; we will also attempt to get stamps at the first iteration, but it - // doesn't matter. - if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2, - stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) { - goto MATCH_ERROR; - } - - - // generate the window function from the set of stamps - if (!pmSubtractionStampsGetWindow(stamps, size)) { - psError(psErrorCodeLast(), false, "Unable to get stamp window."); - goto MATCH_ERROR; - } - - // Define kernel basis functions - if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, - optFWHMs, optOrder, stamps, footprint, - optThreshold, penalty, bounds, subMode); - if (!kernels) { - psErrorClear(); - psWarning("Unable to derive optimum ISIS kernel --- switching to default."); - } - } - if (kernels == NULL) { - // Not an ISIS/GUNK kernel, or the optimum kernel search failed - kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder, penalty, bounds, subMode); - } - - memCheck("kernels"); - - if (subMode == PM_SUBTRACTION_MODE_UNSURE) { - pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); - switch (newMode) { - case PM_SUBTRACTION_MODE_1: - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2."); - break; - case PM_SUBTRACTION_MODE_2: - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); - break; - default: - psError(psErrorCodeLast(), false, "Unable to determine subtraction order."); - goto MATCH_ERROR; - } - subMode = newMode; - } - - int numRejected = -1; // Number of rejected stamps in each iteration - for (int k = 0; (k < iter) && (numRejected != 0); k++) { - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - - if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, - stampThresh1, stampThresh2, stampSpacing, normFrac, - sysError, skyError, size, footprint, subMode)) { - goto MATCH_ERROR; - } - - // generate the window function from the set of stamps - if (!pmSubtractionStampsGetWindow(stamps, size)) { - psError(psErrorCodeLast(), false, "Unable to get stamps window."); - goto MATCH_ERROR; - } + bool tryAgain = true; + while (tryAgain) { + // We get the stamps here; we will also attempt to get stamps at the first iteration, but it + // doesn't matter. + if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2, + stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) { + goto MATCH_ERROR; + } + + // generate the window function from the set of stamps + if (!pmSubtractionStampsGetWindow(&tryAgain, stamps, size)) { + // if we failed, it might be due to the desired normWindow being larger than the current size. + // in this case, just adjust the size and try again. + if (tryAgain) { + // keep the same footprint-size buffer: + int boundary = footprint - size; + size = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2; + footprint = size + boundary; + + // we need to reconstruct everything, so just free the stamps here and retry + psFree(stamps); + } else { + // unrecoverable error + psError(psErrorCodeLast(), false, "Unable to get stamp window."); + goto MATCH_ERROR; + } + } + } + + // Define kernel basis functions + if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { + kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, + optFWHMs, optOrder, stamps, footprint, + optThreshold, penalty, bounds, subMode); + if (!kernels) { + psErrorClear(); + psWarning("Unable to derive optimum ISIS kernel --- switching to default."); + } + } + if (kernels == NULL) { + // Not an ISIS/GUNK kernel, or the optimum kernel search failed + kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, + inner, binning, ringsOrder, penalty, bounds, subMode); + } + + memCheck("kernels"); + + if (subMode == PM_SUBTRACTION_MODE_UNSURE) { + pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); + switch (newMode) { + case PM_SUBTRACTION_MODE_1: + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2."); + break; + case PM_SUBTRACTION_MODE_2: + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); + break; + default: + psError(psErrorCodeLast(), false, "Unable to determine subtraction order."); + goto MATCH_ERROR; + } + subMode = newMode; + } + + int numRejected = -1; // Number of rejected stamps in each iteration + for (int k = 0; (k < iter) && (numRejected != 0); k++) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); + + if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region, + stampThresh1, stampThresh2, stampSpacing, normFrac, + sysError, skyError, size, footprint, subMode)) { + goto MATCH_ERROR; + } + + // Generate the window function from the set of stamps. Since we succeeded + // above to define a normWindow, if we fail here, it is probably unrecoverable. + if (!pmSubtractionStampsGetWindow(NULL, stamps, size)) { + psError(psErrorCodeLast(), false, "Unable to get stamps window."); + goto MATCH_ERROR; + } // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue - psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); - if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } + psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); + if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } // on each iteration, we start from scratch @@ -797,12 +765,12 @@ // reject the deviant stamps based on the stats of the best match - psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej); - if (numRejected < 0) { - psError(psErrorCodeLast(), false, "Unable to reject stamps."); - goto MATCH_ERROR; - } - memCheck(" reject stamps"); - } + psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); + numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej); + if (numRejected < 0) { + psError(psErrorCodeLast(), false, "Unable to reject stamps."); + goto MATCH_ERROR; + } + memCheck(" reject stamps"); + } // apply the best fit so we are ready to roll @@ -811,24 +779,24 @@ goto MATCH_ERROR; } - psFree(stamps); + psFree(stamps); psFree(bestMatch); - memCheck("solution"); - - if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) { - psError(psErrorCodeLast(), false, "Unable to generate QA data"); - goto MATCH_ERROR; - } - memCheck("diag outputs"); - - psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, - kernelError, covarFrac, region, kernels, true, useFFT)) { - psError(psErrorCodeLast(), false, "Unable to convolve image."); - goto MATCH_ERROR; - } - - psFree(kernels); - kernels = NULL; - } + memCheck("solution"); + + if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) { + psError(psErrorCodeLast(), false, "Unable to generate QA data"); + goto MATCH_ERROR; + } + memCheck("diag outputs"); + + psTrace("psModules.imcombine", 2, "Convolving...\n"); + if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, + kernelError, covarFrac, region, kernels, true, useFFT)) { + psError(psErrorCodeLast(), false, "Unable to convolve image."); + goto MATCH_ERROR; + } + + psFree(kernels); + kernels = NULL; + } } psFree(rng); @@ -844,10 +812,10 @@ if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) { - psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); - goto MATCH_ERROR; + psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); + goto MATCH_ERROR; } if (conv2 && !pmSubtractionBorder(conv2->image, conv2->variance, conv2->mask, size, maskBad)) { - psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); - goto MATCH_ERROR; + psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); + goto MATCH_ERROR; } @@ -860,15 +828,15 @@ #ifdef TESTING { - if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { - psFits *fits = psFitsOpen("convolved1.fits", "w"); - psFitsWriteImage(fits, NULL, conv1->image, 0, NULL); - psFitsClose(fits); - } - - if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { - psFits *fits = psFitsOpen("convolved2.fits", "w"); - psFitsWriteImage(fits, NULL, conv2->image, 0, NULL); - psFitsClose(fits); - } + if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { + psFits *fits = psFitsOpen("convolved1.fits", "w"); + psFitsWriteImage(fits, NULL, conv1->image, 0, NULL); + psFitsClose(fits); + } + + if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { + psFits *fits = psFitsOpen("convolved2.fits", "w"); + psFitsWriteImage(fits, NULL, conv2->image, 0, NULL); + psFitsClose(fits); + } } #endif @@ -895,8 +863,8 @@ // increment). static int subtractionOrderWidth(const psKernel *kernel, // Image - float bg, // Background in image - int size, // Maximum size - const psArray *models, // Buffer of models - const psVector *modelSums // Buffer of model sums + float bg, // Background in image + int size, // Maximum size + const psArray *models, // Buffer of models + const psVector *modelSums // Buffer of model sums ) { @@ -911,20 +879,20 @@ psVector *chi2 = psVectorAlloc(size, PS_TYPE_F32); // chi^2 as a function of radius for (int sigma = 0; sigma < size; sigma++) { - double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian - psKernel *model = models->data[sigma]; // Model of interest - for (int y = yMin; y <= yMax; y++) { - for (int x = xMin; x <= xMax; x++) { - sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg); - } - } - float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian - double sumDev2 = 0.0; // Sum of square deviations - for (int y = yMin; y <= yMax; y++) { - for (int x = xMin; x <= xMax; x++) { - float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation - sumDev2 += PS_SQR(dev); - } - } - chi2->data.F32[sigma] = sumDev2; + double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian + psKernel *model = models->data[sigma]; // Model of interest + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg); + } + } + float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian + double sumDev2 = 0.0; // Sum of square deviations + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation + sumDev2 += PS_SQR(dev); + } + } + chi2->data.F32[sigma] = sumDev2; } @@ -933,8 +901,8 @@ float bestChi2 = INFINITY; // Best chi^2 for (int i = 0; i < size; i++) { - if (chi2->data.F32[i] < bestChi2) { - bestChi2 = chi2->data.F32[i]; - bestIndex = i; - } + if (chi2->data.F32[i] < bestChi2) { + bestChi2 = chi2->data.F32[i]; + bestIndex = i; + } } psFree(chi2); @@ -945,6 +913,6 @@ bool pmSubtractionOrderStamp(psVector *ratios, psVector *mask, const pmSubtractionStampList *stamps, - const psArray *models, const psVector *modelSums, - int index, float bg1, float bg2) + const psArray *models, const psVector *modelSums, + int index, float bg1, float bg2) { PS_ASSERT_VECTOR_NON_NULL(ratios, false); @@ -960,5 +928,5 @@ pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE || stamp->status == PM_SUBTRACTION_STAMP_USED, - "We checked this earlier."); + "We checked this earlier."); // Widths of stars @@ -967,11 +935,11 @@ if (width1 == 0 || width2 == 0) { - ratios->data.F32[index] = NAN; - mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff; + ratios->data.F32[index] = NAN; + mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff; } else { - ratios->data.F32[index] = (float)width1 / (float)width2; - mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0; - psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n", - index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]); + ratios->data.F32[index] = (float)width1 / (float)width2; + mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0; + psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n", + index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]); } @@ -1007,60 +975,60 @@ psVector *modelSums = psVectorAlloc(size, PS_TYPE_F64); // Gaussian model sums for (int sigma = 0; sigma < size; sigma++) { - psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model - float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared - double sumGG = 0.0; // Sum of square of Gaussian - for (int y = -size; y <= size; y++) { - int y2 = PS_SQR(y); // y squared - for (int x = -size; x <= size; x++) { - float rad2 = PS_SQR(x) + y2; // Radius squared - float value = expf(-rad2 * invSigma2); // Model value - model->kernel[y][x] = value; - sumGG += PS_SQR(value); - } - } - models->data[sigma] = model; - modelSums->data.F64[sigma] = 1.0 / sumGG; + psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model + float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared + double sumGG = 0.0; // Sum of square of Gaussian + for (int y = -size; y <= size; y++) { + int y2 = PS_SQR(y); // y squared + for (int x = -size; x <= size; x++) { + float rad2 = PS_SQR(x) + y2; // Radius squared + float value = expf(-rad2 * invSigma2); // Model value + model->kernel[y][x] = value; + sumGG += PS_SQR(value); + } + } + models->data[sigma] = model; + modelSums->data.F64[sigma] = 1.0 / sumGG; } // Fit models to stamps for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; - continue; - } - - if (pmSubtractionThreaded()) { - psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER"); - psArrayAdd(job->args, 1, ratios); - psArrayAdd(job->args, 1, mask); - psArrayAdd(job->args, 1, stamps); - psArrayAdd(job->args, 1, models); - psArrayAdd(job->args, 1, modelSums); - PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32); - if (!psThreadJobAddPending(job)) { - return false; - } - } else { - if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) { - psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i); - psFree(models); - psFree(modelSums); - psFree(ratios); - psFree(mask); - return false; - } - } + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; + continue; + } + + if (pmSubtractionThreaded()) { + psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER"); + psArrayAdd(job->args, 1, ratios); + psArrayAdd(job->args, 1, mask); + psArrayAdd(job->args, 1, stamps); + psArrayAdd(job->args, 1, models); + psArrayAdd(job->args, 1, modelSums); + PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32); + if (!psThreadJobAddPending(job)) { + return false; + } + } else { + if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) { + psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i); + psFree(models); + psFree(modelSums); + psFree(ratios); + psFree(mask); + return false; + } + } } if (!psThreadPoolWait(true)) { - psError(psErrorCodeLast(), false, "Error waiting for threads."); - psFree(models); - psFree(modelSums); - psFree(ratios); - psFree(mask); - return false; + psError(psErrorCodeLast(), false, "Error waiting for threads."); + psFree(models); + psFree(modelSums); + psFree(ratios); + psFree(mask); + return false; } @@ -1070,9 +1038,9 @@ psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) { - psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio."); - psFree(mask); - psFree(ratios); - psFree(stats); - return PM_SUBTRACTION_MODE_ERR; + psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio."); + psFree(mask); + psFree(ratios); + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; } psFree(ratios); @@ -1081,6 +1049,6 @@ // XXX raise an error here or not? if (isnan(stats->robustMedian)) { - psFree(stats); - return PM_SUBTRACTION_MODE_ERR; + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; } @@ -1095,9 +1063,9 @@ // Test a subtraction mode by performing a single iteration static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode - pmSubtractionKernels *kernels, // Kernel description - const char *description, // Description for trace - psImage *subMask, // Subtraction mask - float rej // Rejection threshold - ) + pmSubtractionKernels *kernels, // Kernel description + const char *description, // Description for trace + psImage *subMask, // Subtraction mask + float rej // Rejection threshold + ) { assert(stamps); @@ -1114,12 +1082,12 @@ psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); if (!pmSubtractionCalculateEquation(stamps, kernels)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; } psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description); if (!pmSubtractionSolveEquation(kernels, stamps)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; } @@ -1127,6 +1095,6 @@ psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { - psError(psErrorCodeLast(), false, "Unable to calculate deviations."); - return false; + psError(psErrorCodeLast(), false, "Unable to calculate deviations."); + return false; } @@ -1135,38 +1103,38 @@ long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); if (numRejected < 0) { - psError(psErrorCodeLast(), false, "Unable to reject stamps."); - psFree(deviations); - return false; + psError(psErrorCodeLast(), false, "Unable to reject stamps."); + psFree(deviations); + return false; } psFree(deviations); if (numRejected > 0) { - // Allow re-fit with reduced stamps set - psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); - - psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations - if (!deviations) { - psError(psErrorCodeLast(), false, "Unable to calculate deviations."); - return false; - } - psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description); - long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); - if (numRejected < 0) { - psError(psErrorCodeLast(), false, "Unable to reject stamps."); - psFree(deviations); - return false; - } - psFree(deviations); + // Allow re-fit with reduced stamps set + psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); + if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; + } + + psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; + } + psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); + + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(psErrorCodeLast(), false, "Unable to calculate deviations."); + return false; + } + psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description); + long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); + if (numRejected < 0) { + psError(psErrorCodeLast(), false, "Unable to reject stamps."); + psFree(deviations); + return false; + } + psFree(deviations); } # endif @@ -1176,5 +1144,5 @@ pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels, - const psImage *subMask, float rej) + const psImage *subMask, float rej) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR); @@ -1188,9 +1156,9 @@ if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) { - psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1"); - psFree(stamps1); - psFree(kernels1); - psFree(subMask1); - return PM_SUBTRACTION_MODE_ERR; + psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1"); + psFree(stamps1); + psFree(kernels1); + psFree(subMask1); + return PM_SUBTRACTION_MODE_ERR; } psFree(subMask1); @@ -1203,11 +1171,11 @@ if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) { - psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2"); - psFree(stamps2); - psFree(kernels2); - psFree(subMask2); - psFree(stamps1); - psFree(kernels1); - return PM_SUBTRACTION_MODE_ERR; + psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2"); + psFree(stamps2); + psFree(kernels2); + psFree(subMask2); + psFree(stamps1); + psFree(kernels1); + return PM_SUBTRACTION_MODE_ERR; } psFree(subMask2); @@ -1217,14 +1185,14 @@ pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels psLogMsg("psModules.imcombine", PS_LOG_INFO, - "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n", - kernels1->mean, kernels1->rms, kernels1->numStamps, - kernels2->mean, kernels2->rms, kernels2->numStamps); + "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n", + kernels1->mean, kernels1->rms, kernels1->numStamps, + kernels2->mean, kernels2->rms, kernels2->numStamps); if (kernels1->mean < kernels2->mean) { - bestStamps = stamps1; - bestKernels = kernels1; + bestStamps = stamps1; + bestKernels = kernels1; } else { - bestStamps = stamps2; - bestKernels = kernels2; + bestStamps = stamps2; + bestKernels = kernels2; } @@ -1244,5 +1212,5 @@ bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, - float scaleRef, float scaleMin, float scaleMax) + float scaleRef, float scaleMin, float scaleMax) { PS_ASSERT_PTR_NON_NULL(kernelSize, false); @@ -1266,12 +1234,12 @@ if (isfinite(scaleMin) && scale < scaleMin) { - scale = scaleMin; + scale = scaleMin; } if (isfinite(scaleMax) && scale > scaleMax) { - scale = scaleMax; + scale = scaleMax; } for (int i = 0; i < widths->n; i++) { - widths->data.F32[i] *= scale; + widths->data.F32[i] *= scale; } *kernelSize = *kernelSize * scale + 0.5; @@ -1279,5 +1247,5 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, - "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize); + "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize); return true; Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionParams.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionParams.c (revision 30322) @@ -9,6 +9,8 @@ #include "pmErrorCodes.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" +#include "pmSubtraction.h" #include "pmSubtractionStamps.h" -#include "pmSubtraction.h" #include "pmSubtractionParams.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c (revision 30322) @@ -27,6 +27,8 @@ #include "pmSource.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionStamps.h" +#include "pmSubtractionVisual.h" #define STAMP_LIST_BUFFER 20 // Number of stamps to add to list at a time @@ -367,10 +369,63 @@ } - -pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1, - const psImage *image2, const psImage *subMask, - const psRegion *region, float thresh1, float thresh2, - int size, int footprint, float spacing, float normFrac, - float sysErr, float skyErr, pmSubtractionMode mode) +bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 + const psImage *subMask, // Mask for subtraction, or NULL + psImage *variance, // Variance map + const psRegion *region, // Region of interest + float thresh1, // Threshold for stamp finding on readout 1 + float thresh2, // Threshold for stamp finding on readout 2 + float stampSpacing, // Spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysError, // Relative systematic error in images + float skyError, // Relative systematic error in images + int size, // Kernel half-size + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction + ) +{ + PS_ASSERT_PTR_NON_NULL(stamps, false); + PM_ASSERT_READOUT_NON_NULL(ro1, false); + PM_ASSERT_READOUT_NON_NULL(ro2, false); + PS_ASSERT_IMAGE_NON_NULL(subMask, false); + PS_ASSERT_IMAGE_NON_NULL(variance, false); + PS_ASSERT_PTR_NON_NULL(region, false); + + psTrace("psModules.imcombine", 3, "Finding stamps...\n"); + + psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest + + *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, + size, footprint, stampSpacing, normFrac, sysError, skyError, mode); + if (!*stamps) { + psError(psErrorCodeLast(), false, "Unable to find stamps."); + return false; + } + + psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) { + psError(psErrorCodeLast(), false, "Unable to extract stamps."); + return false; + } + + pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1); + return true; +} + +pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, + const psImage *image1, + const psImage *image2, + const psImage *subMask, + const psRegion *region, + float thresh1, + float thresh2, + int size, + int footprint, + float spacing, + float normFrac, + float sysErr, + float skyErr, + pmSubtractionMode mode) { if (!image1 && !image2) { @@ -687,12 +742,28 @@ } - -bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize) +// we are essentially using aperture photometry to determine the photometric match between the +// images. we need to choose an appropriate-sized aperture for this analysis. If it is too +// large, the measurement will be noisy (and possibly biased) due to the sky noise. If it is +// too small, or inconsistent, the measurement will be biased. We use Kron-mag like aperture +// scaled by the first radial moment. +bool pmSubtractionStampsGetWindow(bool *tryAgain, pmSubtractionStampList *stamps, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); + // if we succeed, or fail with an unrecoverable error, do not try again + if (tryAgain) { + *tryAgain = false; + } + int size = stamps->footprint; // Size of postage stamps + // window for moments calculations downstream + psFree (stamps->window); + stamps->window = psKernelAlloc(-size, size, -size, size); + psImageInit(stamps->window->image, 0.0); + + // window1 and window2 are mean stamp images used here to measure the + // first radial moment, and thus the normalization window psFree (stamps->window1); stamps->window1 = psKernelAlloc(-size, size, -size, size); @@ -703,22 +774,16 @@ psImageInit(stamps->window2->image, 0.0); - // Generate a weighting window based on the fwhms (20% larger than the largest) - { - float fwhm1, fwhm2; - - // XXX this is annoyingly hack-ish - pmSubtractionGetFWHMs(&fwhm1, &fwhm2); - - float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35; - - psFree (stamps->window); - stamps->window = psKernelAlloc(-size, size, -size, size); - psImageInit(stamps->window->image, 0.0); - - for (int y = -size; y <= size; y++) { - for (int x = -size; x <= size; x++) { - stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma)); - } - } + // Generate an initial weighting window based on the fwhms (50% larger than the largest) + float fwhm1, fwhm2; + + // XXX this is annoyingly hack-ish + pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + + float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35; + + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma)); + } } @@ -810,6 +875,5 @@ psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2); -# if (1) - // this block attempts to calculate the radius based on the first radial moment + // attempt to calculate the normalization window based on the first radial moment double Sr1 = 0.0; double Sr2 = 0.0; @@ -833,85 +897,47 @@ psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2); + // if the calculated normWindows are too large, we will fall off the stamps. In this case, we need to try again. + if ((stamps->normWindow1 > size) || (stamps->normWindow2 > size)) { + if (tryAgain) { + *tryAgain = true; + } + psFree (stats); + psFree (flux1); + psFree (flux2); + psFree (norm1); + psFree (norm2); + return false; + } + + // this is an unrecoverable error : something really bogus in the data + if (stamps->normWindow1 == 0) { + psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1)."); + psFree (stats); + psFree (flux1); + psFree (flux2); + psFree (norm1); + psFree (norm2); + return false; + } + if (stamps->normWindow2 == 0) { + psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2)."); + psFree (stats); + psFree (flux1); + psFree (flux2); + psFree (norm1); + psFree (norm2); + return false; + } + // Generate a weighting window based on the kron radii - { - float radius = 2.0 * PS_MAX(R1, R2); - - psImageInit(stamps->window->image, 0.0); - - for (int y = -size; y <= size; y++) { - for (int x = -size; x <= size; x++) { - if (hypot(x,y) > radius) continue; - stamps->window->kernel[y][x] = 1.0; - } - } - } - - // complain if normWindow1 or normWindow2 are larger than size - psAssert(stamps->normWindow1 < size, "norm Window 1 too large"); - psAssert(stamps->normWindow2 < size, "norm Window 2 too large"); - -# else - // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent). - // It did not do very well (though a true curve-of-growth analysis might be better...) - bool done1 = false; - bool done2 = false; - double prior1 = 0.0; - double prior2 = 0.0; - double delta1o = 1.0; - double delta2o = 1.0; - for (int radius = 1; radius <= size && !(done1 && done2); radius++) { - double within1 = 0.0; - double within2 = 0.0; - for (int y = -radius; y <= radius; y++) { - for (int x = -radius; x <= radius; x++) { - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { - within1 += stamps->window1->kernel[y][x]; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { - within2 += stamps->window2->kernel[y][x]; - } - } - } - double delta1 = (within1 - prior1) / within1; - if (!done1 && (fabs(delta1) < stamps->normFrac)) { - // interpolate to the radius at which delta2 is normFrac: - stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1); - psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1); - done1 = true; - } - double delta2 = (within2 - prior2) / within2; - if (!done2 && (fabs(delta2) < stamps->normFrac)) { - // interpolate to the radius at which delta2 is normFrac: - stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2); - psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2); - done2 = true; - } - psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2); - - prior1 = within1; - prior2 = within2; - delta1o = delta1; - delta2o = delta2; - - // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) { - // stamps->normWindow1 = radius; - // done1 = true; - // } - // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) { - // stamps->normWindow2 = radius; - // done2 = true; - // } - - } -# endif - - psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2); - if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) { - psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1)."); - return false; - } - if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) { - psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2)."); - return false; + float radius = 2.0 * PS_MAX(R1, R2); + psImageInit(stamps->window->image, 0.0); + + // we use a top-hat window for the moments analysis + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + if (hypot(x,y) > radius) continue; + stamps->window->kernel[y][x] = 1.0; + } } @@ -928,16 +954,4 @@ } } - -#if 0 - { - psFits *fits = NULL; - fits = psFitsOpen ("window1.norm.fits", "w"); - psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL); - psFitsClose (fits); - fits = psFitsOpen ("window2.norm.fits", "w"); - psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL); - psFitsClose (fits); - } -#endif psFree (stats); Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.h (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.h (revision 30322) @@ -1,38 +1,4 @@ #ifndef PM_SUBTRACTION_STAMPS_H #define PM_SUBTRACTION_STAMPS_H - -#include - -#include "pmSubtractionKernels.h" - -/// Status of stamp -typedef enum { - PM_SUBTRACTION_STAMP_INIT, ///< Initial state - PM_SUBTRACTION_STAMP_FOUND, ///< Found a suitable source for this stamp - PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp - PM_SUBTRACTION_STAMP_USED, ///< Use this stamp - PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected - PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region -} pmSubtractionStampStatus; - -/// A list of stamps -typedef struct { - long num; ///< Number of stamps - psArray *stamps; ///< The stamps - psArray *regions; ///< Regions for each stamp - psArray *x, *y; ///< Coordinates for possible stamps (or NULL) - psArray *flux; ///< Fluxes for possible stamps (or NULL) - int footprint; ///< Half-size of stamps - float normFrac; ///< Fraction of flux in window for normalisation window - float normValue; ///< calculated normalization - float normValue2; ///< calculated normalization - psKernel *window1; ///< window function generated from ensemble of stamps (input 1) - psKernel *window2; ///< window function generated from ensemble of stamps (input 2) - psKernel *window; ///< weighting window function (sigma = 1.1 * MAX(fwhm)) - float normWindow1; ///< Size of window for measuring normalisation - float normWindow2; ///< Size of window for measuring normalisation - float sysErr; ///< Systematic error - float skyErr; ///< increase effective readnoise -} pmSubtractionStampList; /// Allocate a list of stamps @@ -74,35 +40,25 @@ ); - -/// A stamp for image subtraction -typedef struct { - float x, y; ///< Position - float flux; ///< Flux - float xNorm, yNorm; ///< Normalised position - psKernel *image1; ///< Reference image postage stamp - psKernel *image2; ///< Input image postage stamp - psKernel *weight; ///< Weight image (1/variance) postage stamp, or NULL - psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL - psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL - psImage *matrix; ///< Least-squares matrix, or NULL - psVector *vector; ///< Least-squares vector, or NULL - double norm; ///< Normalisation difference - double normI1; ///< Sum(flux) for image 1 - double normI2; ///< Sum(flux) for image 2 - double normSquare1; ///< Sum(flux^2) for image 1 (used for penalty) - double normSquare2; ///< Sum(flux^2) for image 2 (used for penalty) - pmSubtractionStampStatus status; ///< Status of stamp - psVector *MxxI1; ///< second moments of convolution images - psVector *MyyI1; ///< second moments of convolution images - psVector *MxxI2; ///< second moments of convolution images - psVector *MyyI2; ///< second moments of convolution images - double MxxI1raw; - double MyyI1raw; - double MxxI2raw; - double MyyI2raw; -} pmSubtractionStamp; - /// Allocate a stamp pmSubtractionStamp *pmSubtractionStampAlloc(void); + +// find and extract the stamps +bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 + const psImage *subMask, // Mask for subtraction, or NULL + psImage *variance, // Variance map + const psRegion *region, // Region of interest + float thresh1, // Threshold for stamp finding on readout 1 + float thresh2, // Threshold for stamp finding on readout 2 + float stampSpacing, // Spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysError, // Relative systematic error in images + float skyError, // Relative systematic error in images + int size, // Kernel half-size + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction + ); + /// Find stamps on an image @@ -172,4 +128,5 @@ /// Calculate the window and normalisation window from the stamps bool pmSubtractionStampsGetWindow( + bool *tryAgain, ///< re-try with new stamp size? pmSubtractionStampList *stamps, ///< List of stamps int kernelSize ///< Half-size of kernel Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionThreads.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionThreads.c (revision 30322) @@ -6,4 +6,5 @@ #include +#include "pmSubtractionTypes.h" #include "pmSubtractionMatch.h" #include "pmSubtractionEquation.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionTypes.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionTypes.h (revision 30322) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionTypes.h (revision 30322) @@ -0,0 +1,168 @@ +/* @file pmSubtraction.h + * + * PSF-matched image subtraction, based on the Alard & Lupton (1998) and Alard (2000) methods. + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ + * Copyright 2004-207 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_SUBTRACTION_TYPES_H +#define PM_SUBTRACTION_TYPES_H + +/// @addtogroup imcombine Image Combinations +/// @{ + +/// Mask values for the subtraction mask +typedef enum { + PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking + PM_SUBTRACTION_MASK_BAD_1 = 0x01, // Image 1 is bad + PM_SUBTRACTION_MASK_BAD_2 = 0x02, // Image 2 is bad + PM_SUBTRACTION_MASK_CONVOLVE_1 = 0x04, // If image 1 is convolved, would be poor or bad + PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be poor or bad + PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad + PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad + PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border + PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected +} pmSubtractionMasks; + +/// Type of subtraction kernel +typedef enum { + PM_SUBTRACTION_KERNEL_NONE, ///< Nothing --- an error + PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions + PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_ISIS_RADIAL, ///< ISIS + higher-order radial Hermitians + PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels + PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved Hermitian polynomial kernels + PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions + PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction + PM_SUBTRACTION_KERNEL_GUNK, ///< Grid United with Normal Kernel --- POIS and ISIS hybrid + PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction +} pmSubtractionKernelsType; + +/// Modes --- specifies which image to convolve +typedef enum { + PM_SUBTRACTION_MODE_ERR, // Error in the mode + PM_SUBTRACTION_MODE_1, // Convolve image 1 + PM_SUBTRACTION_MODE_2, // Convolve image 2 + PM_SUBTRACTION_MODE_UNSURE, // Not sure yet which image to convolve so try to satisfy both + PM_SUBTRACTION_MODE_DUAL, // Dual convolution +} pmSubtractionMode; + +/// Status of stamp +typedef enum { + PM_SUBTRACTION_STAMP_INIT, ///< Initial state + PM_SUBTRACTION_STAMP_FOUND, ///< Found a suitable source for this stamp + PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp + PM_SUBTRACTION_STAMP_USED, ///< Use this stamp + PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected + PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region +} pmSubtractionStampStatus; + +typedef struct { + double score; + pmSubtractionMode mode; + int spatialOrder; + int nGood; + psVector *fluxes; + psVector *chisq; + psVector *moments; + psVector *stampMask; +} pmSubtractionQuality; + +/// Kernels specification +typedef struct { + pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels + psString description; ///< Description of the kernel parameters + int xMin, xMax, yMin, yMax; ///< Bounds of image (for normalisation) + long num; ///< Number of kernel components (not including the spatial ones) + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM) + psVector *widths; ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM) + psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM) + float penalty; ///< Penalty for wideness + psVector *penalties1; ///< Penalty for each kernel component + psVector *penalties2; ///< Penalty for each kernel component + bool havePenalties; ///< flag to test if we have already calculated the penalties or not. + int size; ///< The half-size of the kernel + int inner; ///< The size of an inner region + int spatialOrder; ///< The spatial order of the kernels + int bgOrder; ///< The order for the background fitting + pmSubtractionMode mode; ///< Mode for subtraction + psVector *solution1, *solution2; ///< Solution for the PSF matching + psVector *solution1err, *solution2err; ///< error in solution for the PSF matching + // Quality information + float mean, rms; ///< Mean and RMS of chi^2 from stamps + int numStamps; ///< Number of good stamps + float fResSigmaMean; ///< mean fractional stdev of residuals + float fResSigmaStdev; ///< stdev of fractional stdev of residuals + float fResOuterMean; ///< mean fractional positive swing in residuals + float fResOuterStdev; ///< stdev of fractional positive swing in residuals + float fResTotalMean; ///< mean fractional negative swing in residuals + float fResTotalStdev; ///< stdev of fractional negative swing in residuals + psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations +} pmSubtractionKernels; + +// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures +typedef struct { + psVector *uCoords; // used by RINGS + psVector *vCoords; // used by RINGS + psVector *poly; // used by RINGS + + psVector *xKernel; // used by ISIS, HERM, DECONV_HERM + psVector *yKernel; // used by ISIS, HERM, DECONV_HERM + psKernel *kernel; // used by ISIS, HERM, DECONV_HERM +} pmSubtractionKernelPreCalc; + +/// A list of stamps +typedef struct { + long num; ///< Number of stamps + psArray *stamps; ///< The stamps + psArray *regions; ///< Regions for each stamp + psArray *x, *y; ///< Coordinates for possible stamps (or NULL) + psArray *flux; ///< Fluxes for possible stamps (or NULL) + int footprint; ///< Half-size of stamps + float normFrac; ///< Fraction of flux in window for normalisation window + float normValue; ///< calculated normalization + float normValue2; ///< calculated normalization + psKernel *window1; ///< window function generated from ensemble of stamps (input 1) + psKernel *window2; ///< window function generated from ensemble of stamps (input 2) + psKernel *window; ///< weighting window function (sigma = 1.1 * MAX(fwhm)) + float normWindow1; ///< Size of window for measuring normalisation + float normWindow2; ///< Size of window for measuring normalisation + float sysErr; ///< Systematic error + float skyErr; ///< increase effective readnoise +} pmSubtractionStampList; + +/// A stamp for image subtraction +typedef struct { + float x, y; ///< Position + float flux; ///< Flux + float xNorm, yNorm; ///< Normalised position + psKernel *image1; ///< Reference image postage stamp + psKernel *image2; ///< Input image postage stamp + psKernel *weight; ///< Weight image (1/variance) postage stamp, or NULL + psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL + psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL + psImage *matrix; ///< Least-squares matrix, or NULL + psVector *vector; ///< Least-squares vector, or NULL + double norm; ///< Normalisation difference + double normI1; ///< Sum(flux) for image 1 + double normI2; ///< Sum(flux) for image 2 + double normSquare1; ///< Sum(flux^2) for image 1 (used for penalty) + double normSquare2; ///< Sum(flux^2) for image 2 (used for penalty) + pmSubtractionStampStatus status; ///< Status of stamp + psVector *MxxI1; ///< second moments of convolution images + psVector *MyyI1; ///< second moments of convolution images + psVector *MxxI2; ///< second moments of convolution images + psVector *MyyI2; ///< second moments of convolution images + double MxxI1raw; + double MyyI1raw; + double MxxI2raw; + double MyyI2raw; +} pmSubtractionStamp; + +#endif Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c (revision 30322) @@ -16,4 +16,6 @@ #include "pmKapaPlots.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionStamps.h" Index: /branches/eam_branches/ipp-20101205/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/psmodules.h (revision 30321) +++ /branches/eam_branches/ipp-20101205/psModules/src/psmodules.h (revision 30322) @@ -98,4 +98,5 @@ // the following headers are from psModule:imcombine #include +#include #include #include