IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30322


Ignore:
Timestamp:
Jan 20, 2011, 11:53:03 AM (15 years ago)
Author:
eugene
Message:

moved subtractionGetStamps to pmSubtractionStamps and renamed as pmSubtractionStampsSelect; added return value to pmSubtractionStampsGetWindow to specify a size retry; if desired normalization window is too large, retry with a larger stamp size; move around structures to prevent annoying build collisions

Location:
branches/eam_branches/ipp-20101205/psModules/src
Files:
1 added
22 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101205/psModules/src/extras/pmVisual.c

    r30266 r30322  
    1616#include <pslib.h>
    1717
    18 #include "pmHDU.h"
    19 #include "pmFPA.h"
    20 #include "pmFPAfile.h"
    21 #include "pmAstrometryObjects.h"
    22 #include "pmSubtractionStamps.h"
    23 
    24 #include "pmTrend2D.h"
    25 #include "pmResiduals.h"
    26 #include "pmGrowthCurve.h"
    27 #include "pmSpan.h"
    28 #include "pmFootprintSpans.h"
    29 #include "pmFootprint.h"
    30 #include "pmPeaks.h"
    31 #include "pmMoments.h"
    32 #include "pmModelFuncs.h"
    33 #include "pmModel.h"
    34 #include "pmSourceMasks.h"
    35 #include "pmSourceExtendedPars.h"
    36 #include "pmSourceDiffStats.h"
    37 #include "pmSource.h"
    38 #include "pmSourceFitModel.h"
    39 #include "pmPSF.h"
    40 #include "pmPSFtry.h"
    41 
    42 #include "pmFPAExtent.h"
    43 
    44 #include "pmAstrometryVisual.h"
    45 #include "pmSubtractionVisual.h"
    46 #include "pmStackVisual.h"
    47 #include "pmSourceVisual.h"
     18bool pmSubtractionVisualClose(void);
     19bool pmAstromVisualClose(void);
     20bool pmSubtractionVisualClose(void);
     21bool pmStackVisualClose(void);
     22bool pmSourceVisualClose(void);
     23
     24// #include "pmHDU.h"
     25// #include "pmFPA.h"
     26// #include "pmFPAfile.h"
     27// #include "pmAstrometryObjects.h"
     28// #include "pmSubtractionStamps.h"
     29// #include "pmTrend2D.h"
     30// #include "pmResiduals.h"
     31// #include "pmGrowthCurve.h"
     32// #include "pmSpan.h"
     33// #include "pmFootprintSpans.h"
     34// #include "pmFootprint.h"
     35// #include "pmPeaks.h"
     36// #include "pmMoments.h"
     37// #include "pmModelFuncs.h"
     38// #include "pmModel.h"
     39// #include "pmSourceMasks.h"
     40// #include "pmSourceExtendedPars.h"
     41// #include "pmSourceDiffStats.h"
     42// #include "pmSource.h"
     43// #include "pmSourceFitModel.h"
     44// #include "pmPSF.h"
     45// #include "pmPSFtry.h"
     46// #include "pmFPAExtent.h"
     47// #include "pmAstrometryVisual.h"
     48// #include "pmSubtractionVisual.h"
     49// #include "pmStackVisual.h"
     50// #include "pmSourceVisual.h"
    4851
    4952# if (HAVE_KAPA)

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/Makefile.am

    r26893 r30322  
    3030        pmStackReject.h         \
    3131        pmSubtraction.h         \
     32        pmSubtractionTypes.h            \
    3233        pmSubtractionAnalysis.h \
    3334        pmSubtractionEquation.h \

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmStackReject.c

    r28405 r30322  
    77#include <pslib.h>
    88
     9#include "pmFPA.h"
     10#include "pmSubtractionTypes.h"
    911#include "pmSubtraction.h"
    1012#include "pmSubtractionThreads.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c

    r30303 r30322  
    2020#include "pmHDU.h"                      // Required for pmFPA.h
    2121#include "pmFPA.h"
     22#include "pmSubtractionTypes.h"
    2223#include "pmSubtractionStamps.h"
    2324#include "pmSubtractionEquation.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h

    r30288 r30322  
    1414#define PM_SUBTRACTION_H
    1515
    16 #include <pslib.h>
    17 
    18 #include <pmHDU.h>
    19 #include <pmFPA.h>
    20 #include <pmSubtractionKernels.h>
    21 #include <pmSubtractionStamps.h>
     16// #include <pslib.h>
     17// #include <pmHDU.h>
     18// #include <pmFPA.h>
     19// #include <pmSubtractionKernels.h>
     20// #include <pmSubtractionStamps.h>
    2221
    2322// if we use the original ppSub implementation, we subtract a central delta-function for all
    … …  
    3029/// @addtogroup imcombine Image Combinations
    3130/// @{
    32 
    33 /// Mask values for the subtraction mask
    34 typedef enum {
    35     PM_SUBTRACTION_MASK_CLEAR          = 0x00, // No masking
    36     PM_SUBTRACTION_MASK_BAD_1          = 0x01, // Image 1 is bad
    37     PM_SUBTRACTION_MASK_BAD_2          = 0x02, // Image 2 is bad
    38     PM_SUBTRACTION_MASK_CONVOLVE_1     = 0x04, // If image 1 is convolved, would be poor or bad
    39     PM_SUBTRACTION_MASK_CONVOLVE_2     = 0x08, // If image 2 is convolved, would be poor or bad
    40     PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad
    41     PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad
    42     PM_SUBTRACTION_MASK_BORDER         = 0x40, // Image border
    43     PM_SUBTRACTION_MASK_REJ            = 0x80, // Previously tried as a stamp, and rejected
    44 } pmSubtractionMasks;
    45 
    46 typedef struct {
    47     double score;
    48     pmSubtractionMode mode;
    49     int spatialOrder;
    50     int nGood;
    51     psVector *fluxes;
    52     psVector *chisq;
    53     psVector *moments;
    54     psVector *stampMask;
    55 } pmSubtractionQuality;
    5631
    5732/// Number of terms in a polynomial

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionAnalysis.c

    r29595 r30322  
    77
    88#include "pmFPA.h"
     9#include "pmSubtractionTypes.h"
    910#include "pmSubtraction.h"
    1011#include "pmSubtractionKernels.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionDeconvolve.c

    r26893 r30322  
    1010
    1111#include "pmFPA.h"
     12#include "pmSubtractionTypes.h"
    1213#include "pmSubtractionKernels.h"
    1314#include "pmSubtractionDeconvolve.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c

    r30308 r30322  
    99#include "pmErrorCodes.h"
    1010#include "pmVisual.h"
     11#include "pmFPA.h"
     12#include "pmSubtractionTypes.h"
    1113#include "pmSubtraction.h"
    1214#include "pmSubtractionKernels.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h

    r30289 r30322  
    55#include "pmSubtractionKernels.h"
    66#include "pmSubtraction.h"
    7 
    8 // typedef enum {
    9 //     PM_SUBTRACTION_EQUATION_NONE    = 0x00,
    10 //     PM_SUBTRACTION_EQUATION_NORM    = 0x01,
    11 //     PM_SUBTRACTION_EQUATION_BG      = 0x02,
    12 //     PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
    13 //     PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
    14 // } pmSubtractionEquationCalculationMode;
    157
    168/// Execute a thread job to calculate the least-squares equation for a stamp

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.v0.c

    r29543 r30322  
    99#include "pmErrorCodes.h"
    1010#include "pmSubtraction.h"
     11#include "pmSubtractionTypes.h"
    1112#include "pmSubtractionKernels.h"
    1213#include "pmSubtractionStamps.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionHermitian.c

    r26893 r30322  
    88#include <pslib.h>
    99
     10#include "pmSubtractionTypes.h"
    1011#include "pmSubtractionHermitian.h"
    1112

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionIO.c

    r27321 r30322  
    1212#include "pmConceptsRead.h"
    1313
     14#include "pmSubtractionTypes.h"
    1415#include "pmSubtraction.h"
    1516#include "pmSubtractionKernels.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.c

    r30266 r30322  
    88#include <pslib.h>
    99
     10#include "pmFPA.h"
     11#include "pmSubtractionTypes.h"
    1012#include "pmSubtraction.h"
    1113#include "pmSubtractionKernels.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.h

    r30266 r30322  
    22#define PM_SUBTRACTION_KERNELS_H
    33
    4 #include <string.h>
    5 #include <pslib.h>
    6 
    7 /// Type of subtraction kernel
    8 typedef enum {
    9     PM_SUBTRACTION_KERNEL_NONE,         ///< Nothing --- an error
    10     PM_SUBTRACTION_KERNEL_POIS,         ///< Pan-STARRS Optimal Image Subtraction --- delta functions
    11     PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
    12     PM_SUBTRACTION_KERNEL_ISIS_RADIAL,  ///< ISIS + higher-order radial Hermitians
    13     PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
    14     PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
    15     PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
    16     PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
    17     PM_SUBTRACTION_KERNEL_GUNK,         ///< Grid United with Normal Kernel --- POIS and ISIS hybrid
    18     PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
    19 } pmSubtractionKernelsType;
    20 
    21 /// Modes --- specifies which image to convolve
    22 typedef enum {
    23     PM_SUBTRACTION_MODE_ERR,            // Error in the mode
    24     PM_SUBTRACTION_MODE_1,              // Convolve image 1
    25     PM_SUBTRACTION_MODE_2,              // Convolve image 2
    26     PM_SUBTRACTION_MODE_UNSURE,         // Not sure yet which image to convolve so try to satisfy both
    27     PM_SUBTRACTION_MODE_DUAL,           // Dual convolution
    28 } pmSubtractionMode;
    29 
    30 /// Kernels specification
    31 typedef struct {
    32     pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
    33     psString description;               ///< Description of the kernel parameters
    34     int xMin, xMax, yMin, yMax;         ///< Bounds of image (for normalisation)
    35     long num;                           ///< Number of kernel components (not including the spatial ones)
    36     psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM)
    37     psVector *widths;                   ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM)
    38     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    39     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM)
    40     float penalty;                      ///< Penalty for wideness
    41     psVector *penalties1;               ///< Penalty for each kernel component
    42     psVector *penalties2;               ///< Penalty for each kernel component
    43     bool havePenalties;                 ///< flag to test if we have already calculated the penalties or not.
    44     int size;                           ///< The half-size of the kernel
    45     int inner;                          ///< The size of an inner region
    46     int spatialOrder;                   ///< The spatial order of the kernels
    47     int bgOrder;                        ///< The order for the background fitting
    48     pmSubtractionMode mode;             ///< Mode for subtraction
    49     psVector *solution1, *solution2;    ///< Solution for the PSF matching
    50     psVector *solution1err, *solution2err; ///< error in solution for the PSF matching
    51     // Quality information
    52     float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
    53     int numStamps;                      ///< Number of good stamps
    54     float fResSigmaMean;                ///< mean fractional stdev of residuals
    55     float fResSigmaStdev;               ///< stdev of fractional stdev of residuals
    56     float fResOuterMean;                ///< mean fractional positive swing in residuals
    57     float fResOuterStdev;               ///< stdev of fractional positive swing in residuals
    58     float fResTotalMean;                ///< mean fractional negative swing in residuals
    59     float fResTotalStdev;               ///< stdev of fractional negative swing in residuals
    60     psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
    61 } pmSubtractionKernels;
    62 
    63 // pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures
    64 typedef struct {
    65     psVector *uCoords;                  // used by RINGS
    66     psVector *vCoords;                  // used by RINGS
    67     psVector *poly;                     // used by RINGS
    68 
    69     psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
    70     psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
    71     psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
    72 } pmSubtractionKernelPreCalc;
     4// #include <string.h>
     5// #include <pslib.h>
    736
    747// Assertion to check pmSubtractionKernels

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMask.c

    r27402 r30322  
    77
    88#include "pmErrorCodes.h"
     9#include "pmFPA.h"
     10#include "pmSubtractionTypes.h"
    911#include "pmSubtraction.h"
    1012#include "pmSubtractionKernels.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c

    r30304 r30322  
    1111#include "pmFPA.h"
    1212#include "pmHDUUtils.h"
     13#include "pmSubtractionTypes.h"
     14#include "pmSubtraction.h"
    1315#include "pmSubtractionParams.h"
    1416#include "pmSubtractionKernels.h"
    1517#include "pmSubtractionStamps.h"
    1618#include "pmSubtractionEquation.h"
    17 #include "pmSubtraction.h"
    1819#include "pmSubtractionAnalysis.h"
    1920#include "pmSubtractionMask.h"
    … …  
    5758}
    5859
    59 
    60 static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read
    61                                  const pmReadout *ro1, // Readout 1
    62                                  const pmReadout *ro2, // Readout 2
    63                                  const psImage *subMask, // Mask for subtraction, or NULL
    64                                  psImage *variance,  // Variance map
    65                                  const psRegion *region, // Region of interest
    66                                  float thresh1,  // Threshold for stamp finding on readout 1
    67                                  float thresh2,  // Threshold for stamp finding on readout 2
    68                                  float stampSpacing, // Spacing between stamps
    69                                  float normFrac,     // Fraction of flux in window for normalisation window
    70                                  float sysError,     // Relative systematic error in images
    71                                  float skyError,     // Relative systematic error in images
    72                                  int size,         // Kernel half-size
    73                                  int footprint,     // Convolution footprint for stamps
    74                                  pmSubtractionMode mode // Mode for subtraction
    75     )
    76 {
    77     PS_ASSERT_PTR_NON_NULL(stamps, false);
    78     PM_ASSERT_READOUT_NON_NULL(ro1, false);
    79     PM_ASSERT_READOUT_NON_NULL(ro2, false);
    80     PS_ASSERT_IMAGE_NON_NULL(subMask, false);
    81     PS_ASSERT_IMAGE_NON_NULL(variance, false);
    82     PS_ASSERT_PTR_NON_NULL(region, false);
    83 
    84     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    85 
    86     psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
    87 
    88     *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
    89                                       size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
    90     if (!*stamps) {
    91         psError(psErrorCodeLast(), false, "Unable to find stamps.");
    92         return false;
    93     }
    94 
    95     memCheck("  find stamps");
    96 
    97     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    98     if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
    99         psError(psErrorCodeLast(), false, "Unable to extract stamps.");
    100         return false;
    101     }
    102 
    103     memCheck("   extract stamps");
    104     pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1);
    105     return true;
    106 }
    107 
    10860// Check input arguments
    10961static bool subtractionMatchCheck(pmReadout *conv1, pmReadout *conv2, // Convolved images
    … …  
    12173                                  float badFrac,   // Maximum fraction of bad input pixels to accept
    12274                                  pmSubtractionMode subMode // Mode of subtraction
    123                                   )
     75    )
    12476{
    12577    if (subMode != PM_SUBTRACTION_MODE_2) {
    … …  
    700652            regionString = psRegionToString(*region);
    701653            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
    702                     regionString, numCols, numRows);
     654                     regionString, numCols, numRows);
    703655
    704656            if (stampsName && strlen(stampsName) > 0) {
    … …  
    712664            }
    713665
    714             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    715             // doesn't matter.
    716             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
    717                                       stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
    718                 goto MATCH_ERROR;
    719             }
    720 
    721 
    722             // generate the window function from the set of stamps
    723             if (!pmSubtractionStampsGetWindow(stamps, size)) {
    724                 psError(psErrorCodeLast(), false, "Unable to get stamp window.");
    725                 goto MATCH_ERROR;
    726             }
    727 
    728             // Define kernel basis functions
    729             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    730                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
    731                                                           optFWHMs, optOrder, stamps, footprint,
    732                                                           optThreshold, penalty, bounds, subMode);
    733                 if (!kernels) {
    734                     psErrorClear();
    735                     psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
    736                 }
    737             }
    738             if (kernels == NULL) {
    739                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    740                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    741                                                        inner, binning, ringsOrder, penalty, bounds, subMode);
    742             }
    743 
    744             memCheck("kernels");
    745 
    746             if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
    747                 pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
    748                 switch (newMode) {
    749                   case PM_SUBTRACTION_MODE_1:
    750                     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
    751                     break;
    752                   case PM_SUBTRACTION_MODE_2:
    753                     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
    754                     break;
    755                   default:
    756                     psError(psErrorCodeLast(), false, "Unable to determine subtraction order.");
    757                     goto MATCH_ERROR;
    758                 }
    759                 subMode = newMode;
    760             }
    761 
    762             int numRejected = -1;       // Number of rejected stamps in each iteration
    763             for (int k = 0; (k < iter) && (numRejected != 0); k++) {
    764                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    765 
    766                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    767                                           stampThresh1, stampThresh2, stampSpacing, normFrac,
    768                                           sysError, skyError, size, footprint, subMode)) {
    769                     goto MATCH_ERROR;
    770                 }
    771 
    772                 // generate the window function from the set of stamps
    773                 if (!pmSubtractionStampsGetWindow(stamps, size)) {
    774                     psError(psErrorCodeLast(), false, "Unable to get stamps window.");
    775                     goto MATCH_ERROR;
    776                 }
     666            bool tryAgain = true;
     667            while (tryAgain) {
     668                // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
     669                // doesn't matter.
     670                if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
     671                                               stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
     672                    goto MATCH_ERROR;
     673                }
     674
     675                // generate the window function from the set of stamps
     676                if (!pmSubtractionStampsGetWindow(&tryAgain, stamps, size)) {
     677                    // if we failed, it might be due to the desired normWindow being larger than the current size.
     678                    // in this case, just adjust the size and try again.
     679                    if (tryAgain) {
     680                        // keep the same footprint-size buffer:
     681                        int boundary = footprint - size;
     682                        size = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2;
     683                        footprint = size + boundary;
     684
     685                        // we need to reconstruct everything, so just free the stamps here and retry
     686                        psFree(stamps);
     687                    } else {
     688                        // unrecoverable error
     689                        psError(psErrorCodeLast(), false, "Unable to get stamp window.");
     690                        goto MATCH_ERROR;
     691                    }
     692                }
     693            }
     694
     695            // Define kernel basis functions
     696            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
     697                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
     698                                                          optFWHMs, optOrder, stamps, footprint,
     699                                                          optThreshold, penalty, bounds, subMode);
     700                if (!kernels) {
     701                    psErrorClear();
     702                    psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
     703                }
     704            }
     705            if (kernels == NULL) {
     706                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
     707                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
     708                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
     709            }
     710
     711            memCheck("kernels");
     712
     713            if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
     714                pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
     715                switch (newMode) {
     716                  case PM_SUBTRACTION_MODE_1:
     717                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
     718                    break;
     719                  case PM_SUBTRACTION_MODE_2:
     720                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
     721                    break;
     722                  default:
     723                    psError(psErrorCodeLast(), false, "Unable to determine subtraction order.");
     724                    goto MATCH_ERROR;
     725                }
     726                subMode = newMode;
     727            }
     728
     729            int numRejected = -1;       // Number of rejected stamps in each iteration
     730            for (int k = 0; (k < iter) && (numRejected != 0); k++) {
     731                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
     732
     733                if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region,
     734                                               stampThresh1, stampThresh2, stampSpacing, normFrac,
     735                                               sysError, skyError, size, footprint, subMode)) {
     736                    goto MATCH_ERROR;
     737                }
     738
     739                // Generate the window function from the set of stamps.  Since we succeeded
     740                // above to define a normWindow, if we fail here, it is probably unrecoverable.
     741                if (!pmSubtractionStampsGetWindow(NULL, stamps, size)) {
     742                    psError(psErrorCodeLast(), false, "Unable to get stamps window.");
     743                    goto MATCH_ERROR;
     744                }
    777745
    778746                // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue
    779                 psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
    780                 if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) {
    781                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    782                     goto MATCH_ERROR;
    783                 }
     747                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
     748                if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) {
     749                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     750                    goto MATCH_ERROR;
     751                }
    784752
    785753                // on each iteration, we start from scratch
    … …  
    797765               
    798766                // reject the deviant stamps based on the stats of the best match
    799                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
    800                 numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej);
    801                 if (numRejected < 0) {
    802                     psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    803                     goto MATCH_ERROR;
    804                 }
    805                 memCheck("  reject stamps");
    806             }
     767                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
     768                numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej);
     769                if (numRejected < 0) {
     770                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
     771                    goto MATCH_ERROR;
     772                }
     773                memCheck("  reject stamps");
     774            }
    807775
    808776            // apply the best fit so we are ready to roll
    … …  
    811779                goto MATCH_ERROR;
    812780            }
    813             psFree(stamps);
     781            psFree(stamps);
    814782            psFree(bestMatch);
    815             memCheck("solution");
    816 
    817             if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
    818                 psError(psErrorCodeLast(), false, "Unable to generate QA data");
    819                 goto MATCH_ERROR;
    820             }
    821             memCheck("diag outputs");
    822 
    823             psTrace("psModules.imcombine", 2, "Convolving...\n");
    824             if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
    825                                        kernelError, covarFrac, region, kernels, true, useFFT)) {
    826                 psError(psErrorCodeLast(), false, "Unable to convolve image.");
    827                 goto MATCH_ERROR;
    828             }
    829 
    830             psFree(kernels);
    831             kernels = NULL;
    832         }
     783            memCheck("solution");
     784
     785            if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
     786                psError(psErrorCodeLast(), false, "Unable to generate QA data");
     787                goto MATCH_ERROR;
     788            }
     789            memCheck("diag outputs");
     790
     791            psTrace("psModules.imcombine", 2, "Convolving...\n");
     792            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
     793                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
     794                psError(psErrorCodeLast(), false, "Unable to convolve image.");
     795                goto MATCH_ERROR;
     796            }
     797
     798            psFree(kernels);
     799            kernels = NULL;
     800        }
    833801    }
    834802    psFree(rng);
    … …  
    844812
    845813    if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
    846         psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
    847         goto MATCH_ERROR;
     814        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
     815        goto MATCH_ERROR;
    848816    }
    849817    if (conv2 && !pmSubtractionBorder(conv2->image, conv2->variance, conv2->mask, size, maskBad)) {
    850         psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
    851         goto MATCH_ERROR;
     818        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
     819        goto MATCH_ERROR;
    852820    }
    853821
    … …  
    860828#ifdef TESTING
    861829    {
    862         if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    863             psFits *fits = psFitsOpen("convolved1.fits", "w");
    864             psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
    865             psFitsClose(fits);
    866         }
    867 
    868         if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    869             psFits *fits = psFitsOpen("convolved2.fits", "w");
    870             psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);
    871             psFitsClose(fits);
    872         }
     830        if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
     831            psFits *fits = psFitsOpen("convolved1.fits", "w");
     832            psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
     833            psFitsClose(fits);
     834        }
     835
     836        if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
     837            psFits *fits = psFitsOpen("convolved2.fits", "w");
     838            psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);
     839            psFitsClose(fits);
     840        }
    873841    }
    874842#endif
    … …  
    895863// increment).
    896864static int subtractionOrderWidth(const psKernel *kernel, // Image
    897                                 float bg, // Background in image
    898                                 int size, // Maximum size
    899                                 const psArray *models, // Buffer of models
    900                                 const psVector *modelSums // Buffer of model sums
     865                                float bg, // Background in image
     866                                int size, // Maximum size
     867                                const psArray *models, // Buffer of models
     868                                const psVector *modelSums // Buffer of model sums
    901869    )
    902870{
    … …  
    911879    psVector *chi2 = psVectorAlloc(size, PS_TYPE_F32); // chi^2 as a function of radius
    912880    for (int sigma = 0; sigma < size; sigma++) {
    913         double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian
    914         psKernel *model = models->data[sigma]; // Model of interest
    915         for (int y = yMin; y <= yMax; y++) {
    916             for (int x = xMin; x <= xMax; x++) {
    917                 sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg);
    918             }
    919         }
    920         float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian
    921         double sumDev2 = 0.0;           // Sum of square deviations
    922         for (int y = yMin; y <= yMax; y++) {
    923             for (int x = xMin; x <= xMax; x++) {
    924                 float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation
    925                 sumDev2 += PS_SQR(dev);
    926             }
    927         }
    928         chi2->data.F32[sigma] = sumDev2;
     881        double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian
     882        psKernel *model = models->data[sigma]; // Model of interest
     883        for (int y = yMin; y <= yMax; y++) {
     884            for (int x = xMin; x <= xMax; x++) {
     885                sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg);
     886            }
     887        }
     888        float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian
     889        double sumDev2 = 0.0;           // Sum of square deviations
     890        for (int y = yMin; y <= yMax; y++) {
     891            for (int x = xMin; x <= xMax; x++) {
     892                float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation
     893                sumDev2 += PS_SQR(dev);
     894            }
     895        }
     896        chi2->data.F32[sigma] = sumDev2;
    929897    }
    930898
    … …  
    933901    float bestChi2 = INFINITY;          // Best chi^2
    934902    for (int i = 0; i < size; i++) {
    935         if (chi2->data.F32[i] < bestChi2) {
    936             bestChi2 = chi2->data.F32[i];
    937             bestIndex = i;
    938         }
     903        if (chi2->data.F32[i] < bestChi2) {
     904            bestChi2 = chi2->data.F32[i];
     905            bestIndex = i;
     906        }
    939907    }
    940908    psFree(chi2);
    … …  
    945913
    946914bool pmSubtractionOrderStamp(psVector *ratios, psVector *mask, const pmSubtractionStampList *stamps,
    947                              const psArray *models, const psVector *modelSums,
    948                              int index, float bg1, float bg2)
     915                             const psArray *models, const psVector *modelSums,
     916                             int index, float bg1, float bg2)
    949917{
    950918    PS_ASSERT_VECTOR_NON_NULL(ratios, false);
    … …  
    960928    pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
    961929    psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE || stamp->status == PM_SUBTRACTION_STAMP_USED,
    962              "We checked this earlier.");
     930             "We checked this earlier.");
    963931
    964932    // Widths of stars
    … …  
    967935
    968936    if (width1 == 0 || width2 == 0) {
    969         ratios->data.F32[index] = NAN;
    970         mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff;
     937        ratios->data.F32[index] = NAN;
     938        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff;
    971939    } else {
    972         ratios->data.F32[index] = (float)width1 / (float)width2;
    973         mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0;
    974         psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n",
    975                 index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]);
     940        ratios->data.F32[index] = (float)width1 / (float)width2;
     941        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0;
     942        psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n",
     943                index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]);
    976944    }
    977945
    … …  
    1007975    psVector *modelSums = psVectorAlloc(size, PS_TYPE_F64); // Gaussian model sums
    1008976    for (int sigma = 0; sigma < size; sigma++) {
    1009         psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model
    1010         float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared
    1011         double sumGG = 0.0;         // Sum of square of Gaussian
    1012         for (int y = -size; y <= size; y++) {
    1013             int y2 = PS_SQR(y);     // y squared
    1014             for (int x = -size; x <= size; x++) {
    1015                 float rad2 = PS_SQR(x) + y2; // Radius squared
    1016                 float value = expf(-rad2 * invSigma2); // Model value
    1017                 model->kernel[y][x] = value;
    1018                 sumGG += PS_SQR(value);
    1019             }
    1020         }
    1021         models->data[sigma] = model;
    1022         modelSums->data.F64[sigma] = 1.0 / sumGG;
     977        psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model
     978        float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared
     979        double sumGG = 0.0;         // Sum of square of Gaussian
     980        for (int y = -size; y <= size; y++) {
     981            int y2 = PS_SQR(y);     // y squared
     982            for (int x = -size; x <= size; x++) {
     983                float rad2 = PS_SQR(x) + y2; // Radius squared
     984                float value = expf(-rad2 * invSigma2); // Model value
     985                model->kernel[y][x] = value;
     986                sumGG += PS_SQR(value);
     987            }
     988        }
     989        models->data[sigma] = model;
     990        modelSums->data.F64[sigma] = 1.0 / sumGG;
    1023991    }
    1024992
    1025993    // Fit models to stamps
    1026994    for (int i = 0; i < stamps->num; i++) {
    1027         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    1028         if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
    1029             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
    1030             continue;
    1031         }
    1032 
    1033         if (pmSubtractionThreaded()) {
    1034             psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER");
    1035             psArrayAdd(job->args, 1, ratios);
    1036             psArrayAdd(job->args, 1, mask);
    1037             psArrayAdd(job->args, 1, stamps);
    1038             psArrayAdd(job->args, 1, models);
    1039             psArrayAdd(job->args, 1, modelSums);
    1040             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
    1041             PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32);
    1042             PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
    1043             if (!psThreadJobAddPending(job)) {
    1044                 return false;
    1045             }
    1046         } else {
    1047             if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
    1048                 psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
    1049                 psFree(models);
    1050                 psFree(modelSums);
    1051                 psFree(ratios);
    1052                 psFree(mask);
    1053                 return false;
    1054             }
    1055         }
     995        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     996        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
     997            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
     998            continue;
     999        }
     1000
     1001        if (pmSubtractionThreaded()) {
     1002            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER");
     1003            psArrayAdd(job->args, 1, ratios);
     1004            psArrayAdd(job->args, 1, mask);
     1005            psArrayAdd(job->args, 1, stamps);
     1006            psArrayAdd(job->args, 1, models);
     1007            psArrayAdd(job->args, 1, modelSums);
     1008            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
     1009            PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32);
     1010            PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
     1011            if (!psThreadJobAddPending(job)) {
     1012                return false;
     1013            }
     1014        } else {
     1015            if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
     1016                psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
     1017                psFree(models);
     1018                psFree(modelSums);
     1019                psFree(ratios);
     1020                psFree(mask);
     1021                return false;
     1022            }
     1023        }
    10561024    }
    10571025
    10581026    if (!psThreadPoolWait(true)) {
    1059         psError(psErrorCodeLast(), false, "Error waiting for threads.");
    1060         psFree(models);
    1061         psFree(modelSums);
    1062         psFree(ratios);
    1063         psFree(mask);
    1064             return false;
     1027        psError(psErrorCodeLast(), false, "Error waiting for threads.");
     1028        psFree(models);
     1029        psFree(modelSums);
     1030        psFree(ratios);
     1031        psFree(mask);
     1032        return false;
    10651033    }
    10661034
    … …  
    10701038    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    10711039    if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) {
    1072         psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
    1073         psFree(mask);
    1074         psFree(ratios);
    1075         psFree(stats);
    1076         return PM_SUBTRACTION_MODE_ERR;
     1040        psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
     1041        psFree(mask);
     1042        psFree(ratios);
     1043        psFree(stats);
     1044        return PM_SUBTRACTION_MODE_ERR;
    10771045    }
    10781046    psFree(ratios);
    … …  
    10811049    // XXX raise an error here or not?
    10821050    if (isnan(stats->robustMedian)) {
    1083         psFree(stats);
    1084         return PM_SUBTRACTION_MODE_ERR;
     1051        psFree(stats);
     1052        return PM_SUBTRACTION_MODE_ERR;
    10851053    }
    10861054
    … …  
    10951063// Test a subtraction mode by performing a single iteration
    10961064static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode
    1097                                 pmSubtractionKernels *kernels, // Kernel description
    1098                                 const char *description, // Description for trace
    1099                                 psImage *subMask,  // Subtraction mask
    1100                                 float rej               // Rejection threshold
    1101                                 )
     1065                                pmSubtractionKernels *kernels, // Kernel description
     1066                                const char *description, // Description for trace
     1067                                psImage *subMask,  // Subtraction mask
     1068                                float rej               // Rejection threshold
     1069    )
    11021070{
    11031071    assert(stamps);
    … …  
    11141082    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    11151083    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    1116         psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1117         return false;
     1084        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1085        return false;
    11181086    }
    11191087
    11201088    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
    11211089    if (!pmSubtractionSolveEquation(kernels, stamps)) {
    1122         psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1123         return false;
     1090        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1091        return false;
    11241092    }
    11251093
    … …  
    11271095    psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    11281096    if (!deviations) {
    1129         psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    1130         return false;
     1097        psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
     1098        return false;
    11311099    }
    11321100
    … …  
    11351103    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    11361104    if (numRejected < 0) {
    1137         psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    1138         psFree(deviations);
    1139         return false;
     1105        psError(psErrorCodeLast(), false, "Unable to reject stamps.");
     1106        psFree(deviations);
     1107        return false;
    11401108    }
    11411109    psFree(deviations);
    11421110
    11431111    if (numRejected > 0) {
    1144         // Allow re-fit with reduced stamps set
    1145         psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1146         if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    1147             psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1148             return false;
    1149         }
    1150 
    1151         psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    1152         if (!pmSubtractionSolveEquation(kernels, stamps)) {
    1153             psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1154             return false;
    1155         }
    1156         psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
    1157 
    1158         psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    1159         if (!deviations) {
    1160             psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    1161             return false;
    1162         }
    1163         psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description);
    1164         long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
    1165         if (numRejected < 0) {
    1166             psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    1167             psFree(deviations);
    1168             return false;
    1169         }
    1170         psFree(deviations);
     1112        // Allow re-fit with reduced stamps set
     1113        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1114        if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     1115            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1116            return false;
     1117        }
     1118
     1119        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
     1120        if (!pmSubtractionSolveEquation(kernels, stamps)) {
     1121            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1122            return false;
     1123        }
     1124        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
     1125
     1126        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
     1127        if (!deviations) {
     1128            psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
     1129            return false;
     1130        }
     1131        psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description);
     1132        long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
     1133        if (numRejected < 0) {
     1134            psError(psErrorCodeLast(), false, "Unable to reject stamps.");
     1135            psFree(deviations);
     1136            return false;
     1137        }
     1138        psFree(deviations);
    11711139    }
    11721140# endif
    … …  
    11761144
    11771145pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels,
    1178                                         const psImage *subMask, float rej)
     1146                                        const psImage *subMask, float rej)
    11791147{
    11801148    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR);
    … …  
    11881156
    11891157    if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) {
    1190         psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1");
    1191         psFree(stamps1);
    1192         psFree(kernels1);
    1193         psFree(subMask1);
    1194         return PM_SUBTRACTION_MODE_ERR;
     1158        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1");
     1159        psFree(stamps1);
     1160        psFree(kernels1);
     1161        psFree(subMask1);
     1162        return PM_SUBTRACTION_MODE_ERR;
    11951163    }
    11961164    psFree(subMask1);
    … …  
    12031171
    12041172    if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) {
    1205         psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2");
    1206         psFree(stamps2);
    1207         psFree(kernels2);
    1208         psFree(subMask2);
    1209         psFree(stamps1);
    1210         psFree(kernels1);
    1211         return PM_SUBTRACTION_MODE_ERR;
     1173        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2");
     1174        psFree(stamps2);
     1175        psFree(kernels2);
     1176        psFree(subMask2);
     1177        psFree(stamps1);
     1178        psFree(kernels1);
     1179        return PM_SUBTRACTION_MODE_ERR;
    12121180    }
    12131181    psFree(subMask2);
    … …  
    12171185    pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels
    12181186    psLogMsg("psModules.imcombine", PS_LOG_INFO,
    1219              "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n",
    1220              kernels1->mean, kernels1->rms, kernels1->numStamps,
    1221              kernels2->mean, kernels2->rms, kernels2->numStamps);
     1187             "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n",
     1188             kernels1->mean, kernels1->rms, kernels1->numStamps,
     1189             kernels2->mean, kernels2->rms, kernels2->numStamps);
    12221190
    12231191    if (kernels1->mean < kernels2->mean) {
    1224         bestStamps = stamps1;
    1225         bestKernels = kernels1;
     1192        bestStamps = stamps1;
     1193        bestKernels = kernels1;
    12261194    } else {
    1227         bestStamps = stamps2;
    1228         bestKernels = kernels2;
     1195        bestStamps = stamps2;
     1196        bestKernels = kernels2;
    12291197    }
    12301198
    … …  
    12441212
    12451213bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
    1246                               float scaleRef, float scaleMin, float scaleMax)
     1214                              float scaleRef, float scaleMin, float scaleMax)
    12471215{
    12481216    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
    … …  
    12661234
    12671235    if (isfinite(scaleMin) && scale < scaleMin) {
    1268         scale = scaleMin;
     1236        scale = scaleMin;
    12691237    }
    12701238    if (isfinite(scaleMax) && scale > scaleMax) {
    1271         scale = scaleMax;
     1239        scale = scaleMax;
    12721240    }
    12731241
    12741242    for (int i = 0; i < widths->n; i++) {
    1275         widths->data.F32[i] *= scale;
     1243        widths->data.F32[i] *= scale;
    12761244    }
    12771245    *kernelSize = *kernelSize * scale + 0.5;
    … …  
    12791247
    12801248    psLogMsg("psModules.imcombine", PS_LOG_INFO,
    1281              "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
     1249             "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
    12821250
    12831251    return true;

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionParams.c

    r30286 r30322  
    99
    1010#include "pmErrorCodes.h"
     11#include "pmFPA.h"
     12#include "pmSubtractionTypes.h"
     13#include "pmSubtraction.h"
    1114#include "pmSubtractionStamps.h"
    12 #include "pmSubtraction.h"
    1315#include "pmSubtractionParams.h"
    1416

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c

    r30289 r30322  
    2727#include "pmSource.h"
    2828
     29#include "pmSubtractionTypes.h"
    2930#include "pmSubtraction.h"
    3031#include "pmSubtractionStamps.h"
     32#include "pmSubtractionVisual.h"
    3133
    3234#define STAMP_LIST_BUFFER 20            // Number of stamps to add to list at a time
    … …  
    367369}
    368370
    369 
    370 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1,
    371                                                 const psImage *image2, const psImage *subMask,
    372                                                 const psRegion *region, float thresh1, float thresh2,
    373                                                 int size, int footprint, float spacing, float normFrac,
    374                                                 float sysErr, float skyErr, pmSubtractionMode mode)
     371bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read
     372                               const pmReadout *ro1, // Readout 1
     373                               const pmReadout *ro2, // Readout 2
     374                               const psImage *subMask, // Mask for subtraction, or NULL
     375                               psImage *variance,  // Variance map
     376                               const psRegion *region, // Region of interest
     377                               float thresh1,  // Threshold for stamp finding on readout 1
     378                               float thresh2,  // Threshold for stamp finding on readout 2
     379                               float stampSpacing, // Spacing between stamps
     380                               float normFrac,     // Fraction of flux in window for normalisation window
     381                               float sysError,     // Relative systematic error in images
     382                               float skyError,     // Relative systematic error in images
     383                               int size,         // Kernel half-size
     384                               int footprint,     // Convolution footprint for stamps
     385                               pmSubtractionMode mode // Mode for subtraction
     386    )
     387{
     388    PS_ASSERT_PTR_NON_NULL(stamps, false);
     389    PM_ASSERT_READOUT_NON_NULL(ro1, false);
     390    PM_ASSERT_READOUT_NON_NULL(ro2, false);
     391    PS_ASSERT_IMAGE_NON_NULL(subMask, false);
     392    PS_ASSERT_IMAGE_NON_NULL(variance, false);
     393    PS_ASSERT_PTR_NON_NULL(region, false);
     394
     395    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
     396
     397    psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
     398
     399    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
     400                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
     401    if (!*stamps) {
     402        psError(psErrorCodeLast(), false, "Unable to find stamps.");
     403        return false;
     404    }
     405
     406    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
     407    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
     408        psError(psErrorCodeLast(), false, "Unable to extract stamps.");
     409        return false;
     410    }
     411
     412    pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1);
     413    return true;
     414}
     415
     416pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps,
     417                                                const psImage *image1,
     418                                                const psImage *image2,
     419                                                const psImage *subMask,
     420                                                const psRegion *region,
     421                                                float thresh1,
     422                                                float thresh2,
     423                                                int size,
     424                                                int footprint,
     425                                                float spacing,
     426                                                float normFrac,
     427                                                float sysErr,
     428                                                float skyErr,
     429                                                pmSubtractionMode mode)
    375430{
    376431    if (!image1 && !image2) {
    … …  
    687742}
    688743
    689 
    690 bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize)
     744// we are essentially using aperture photometry to determine the photometric match between the
     745// images.  we need to choose an appropriate-sized aperture for this analysis.  If it is too
     746// large, the measurement will be noisy (and possibly biased) due to the sky noise.  If it is
     747// too small, or inconsistent, the measurement will be biased.  We use Kron-mag like aperture
     748// scaled by the first radial moment.
     749bool pmSubtractionStampsGetWindow(bool *tryAgain, pmSubtractionStampList *stamps, int kernelSize)
    691750{
    692751    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    693752    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    694753
     754    // if we succeed, or fail with an unrecoverable error, do not try again
     755    if (tryAgain) {
     756        *tryAgain = false;
     757    }
     758
    695759    int size = stamps->footprint; // Size of postage stamps
    696760
     761    // window for moments calculations downstream
     762    psFree (stamps->window);
     763    stamps->window = psKernelAlloc(-size, size, -size, size);
     764    psImageInit(stamps->window->image, 0.0);
     765
     766    // window1 and window2 are mean stamp images used here to measure the
     767    // first radial moment, and thus the normalization window
    697768    psFree (stamps->window1);
    698769    stamps->window1 = psKernelAlloc(-size, size, -size, size);
    … …  
    703774    psImageInit(stamps->window2->image, 0.0);
    704775
    705     // Generate a weighting window based on the fwhms (20% larger than the largest)
    706     {
    707         float fwhm1, fwhm2;
    708 
    709         // XXX this is annoyingly hack-ish
    710         pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
    711 
    712         float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
    713 
    714         psFree (stamps->window);
    715         stamps->window = psKernelAlloc(-size, size, -size, size);
    716         psImageInit(stamps->window->image, 0.0);
    717 
    718         for (int y = -size; y <= size; y++) {
    719             for (int x = -size; x <= size; x++) {
    720                 stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
    721             }
    722         }
     776    // Generate an initial weighting window based on the fwhms (50% larger than the largest)
     777    float fwhm1, fwhm2;
     778
     779    // XXX this is annoyingly hack-ish
     780    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     781   
     782    float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
     783   
     784    for (int y = -size; y <= size; y++) {
     785        for (int x = -size; x <= size; x++) {
     786            stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
     787        }
    723788    }
    724789
    … …  
    810875    psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2);
    811876
    812 # if (1)
    813     // this block attempts to calculate the radius based on the first radial moment
     877    // attempt to calculate the normalization window based on the first radial moment
    814878    double Sr1 = 0.0;
    815879    double Sr2 = 0.0;
    … …  
    833897    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
    834898
     899    // if the calculated normWindows are too large, we will fall off the stamps.  In this case, we need to try again.
     900    if ((stamps->normWindow1 > size) || (stamps->normWindow2 > size)) {
     901        if (tryAgain) {
     902            *tryAgain = true;
     903        }
     904        psFree (stats);
     905        psFree (flux1);
     906        psFree (flux2);
     907        psFree (norm1);
     908        psFree (norm2);
     909        return false;
     910    }
     911
     912    // this is an unrecoverable error : something really bogus in the data
     913    if (stamps->normWindow1 == 0) {
     914        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1).");
     915        psFree (stats);
     916        psFree (flux1);
     917        psFree (flux2);
     918        psFree (norm1);
     919        psFree (norm2);
     920        return false;
     921    }
     922    if (stamps->normWindow2 == 0) {
     923        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2).");
     924        psFree (stats);
     925        psFree (flux1);
     926        psFree (flux2);
     927        psFree (norm1);
     928        psFree (norm2);
     929        return false;
     930    }
     931
    835932    // Generate a weighting window based on the kron radii
    836     {
    837         float radius = 2.0 * PS_MAX(R1, R2);
    838 
    839         psImageInit(stamps->window->image, 0.0);
    840 
    841         for (int y = -size; y <= size; y++) {
    842             for (int x = -size; x <= size; x++) {
    843                 if (hypot(x,y) > radius) continue;
    844                 stamps->window->kernel[y][x] = 1.0;
    845             }
    846         }
    847     }
    848 
    849     // complain if normWindow1 or normWindow2 are larger than size
    850     psAssert(stamps->normWindow1 < size, "norm Window 1 too large");
    851     psAssert(stamps->normWindow2 < size, "norm Window 2 too large");
    852 
    853 # else
    854     // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent).
    855     // It did not do very well (though a true curve-of-growth analysis might be better...)
    856     bool done1 = false;
    857     bool done2 = false;
    858     double prior1 = 0.0;
    859     double prior2 = 0.0;
    860     double delta1o = 1.0;
    861     double delta2o = 1.0;
    862     for (int radius = 1; radius <= size && !(done1 && done2); radius++) {
    863         double within1 = 0.0;
    864         double within2 = 0.0;
    865         for (int y = -radius; y <= radius; y++) {
    866             for (int x = -radius; x <= radius; x++) {
    867                 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
    868                     within1 += stamps->window1->kernel[y][x];
    869                 }
    870                 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
    871                     within2 += stamps->window2->kernel[y][x];
    872                 }
    873             }
    874         }
    875         double delta1 = (within1 - prior1) / within1;
    876         if (!done1 && (fabs(delta1) < stamps->normFrac)) {
    877             // interpolate to the radius at which delta2 is normFrac:
    878             stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
    879             psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
    880             done1 = true;
    881         }
    882         double delta2 = (within2 - prior2) / within2;
    883         if (!done2 && (fabs(delta2) < stamps->normFrac)) {
    884             // interpolate to the radius at which delta2 is normFrac:
    885             stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
    886             psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
    887             done2 = true;
    888         }
    889         psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2);
    890 
    891         prior1 = within1;
    892         prior2 = within2;
    893         delta1o = delta1;
    894         delta2o = delta2;
    895 
    896         // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) {
    897         //     stamps->normWindow1 = radius;
    898         //     done1 = true;
    899         // }
    900         // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) {
    901         //     stamps->normWindow2 = radius;
    902         //     done2 = true;
    903         // }
    904 
    905     }
    906 # endif
    907 
    908     psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2);
    909     if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) {
    910         psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1).");
    911         return false;
    912     }
    913     if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) {
    914         psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2).");
    915         return false;
     933    float radius = 2.0 * PS_MAX(R1, R2);
     934    psImageInit(stamps->window->image, 0.0);
     935
     936    // we use a top-hat window for the moments analysis
     937    for (int y = -size; y <= size; y++) {
     938        for (int x = -size; x <= size; x++) {
     939            if (hypot(x,y) > radius) continue;
     940            stamps->window->kernel[y][x] = 1.0;
     941        }
    916942    }
    917943
    … …  
    928954        }
    929955    }
    930 
    931 #if 0
    932     {
    933         psFits *fits = NULL;
    934         fits = psFitsOpen ("window1.norm.fits", "w");
    935         psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
    936         psFitsClose (fits);
    937         fits = psFitsOpen ("window2.norm.fits", "w");
    938         psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
    939         psFitsClose (fits);
    940     }
    941 #endif
    942956
    943957    psFree (stats);

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.h

    r29543 r30322  
    11#ifndef PM_SUBTRACTION_STAMPS_H
    22#define PM_SUBTRACTION_STAMPS_H
    3 
    4 #include <pslib.h>
    5 
    6 #include "pmSubtractionKernels.h"
    7 
    8 /// Status of stamp
    9 typedef enum {
    10     PM_SUBTRACTION_STAMP_INIT,          ///< Initial state
    11     PM_SUBTRACTION_STAMP_FOUND,         ///< Found a suitable source for this stamp
    12     PM_SUBTRACTION_STAMP_CALCULATE,     ///< Calculate matrix and vector values for this stamp
    13     PM_SUBTRACTION_STAMP_USED,          ///< Use this stamp
    14     PM_SUBTRACTION_STAMP_REJECTED,      ///< This stamp has been rejected
    15     PM_SUBTRACTION_STAMP_NONE           ///< No stamp in this region
    16 } pmSubtractionStampStatus;
    17 
    18 /// A list of stamps
    19 typedef struct {
    20     long num;                           ///< Number of stamps
    21     psArray *stamps;                    ///< The stamps
    22     psArray *regions;                   ///< Regions for each stamp
    23     psArray *x, *y;                     ///< Coordinates for possible stamps (or NULL)
    24     psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
    25     int footprint;                      ///< Half-size of stamps
    26     float normFrac;                     ///< Fraction of flux in window for normalisation window
    27     float normValue;                    ///< calculated normalization
    28     float normValue2;                   ///< calculated normalization
    29     psKernel *window1;                  ///< window function generated from ensemble of stamps (input 1)
    30     psKernel *window2;                  ///< window function generated from ensemble of stamps (input 2)
    31     psKernel *window;                   ///< weighting window function (sigma = 1.1 * MAX(fwhm))
    32     float normWindow1;                  ///< Size of window for measuring normalisation
    33     float normWindow2;                  ///< Size of window for measuring normalisation
    34     float sysErr;                       ///< Systematic error
    35     float skyErr;                       ///< increase effective readnoise
    36 } pmSubtractionStampList;
    373
    384/// Allocate a list of stamps
    … …  
    7440    );
    7541
    76 
    77 /// A stamp for image subtraction
    78 typedef struct {
    79     float x, y;                         ///< Position
    80     float flux;                         ///< Flux
    81     float xNorm, yNorm;                 ///< Normalised position
    82     psKernel *image1;                   ///< Reference image postage stamp
    83     psKernel *image2;                   ///< Input image postage stamp
    84     psKernel *weight;                   ///< Weight image (1/variance) postage stamp, or NULL
    85     psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
    86     psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
    87     psImage *matrix;                    ///< Least-squares matrix, or NULL
    88     psVector *vector;                   ///< Least-squares vector, or NULL
    89     double norm;                        ///< Normalisation difference
    90     double normI1;                       ///< Sum(flux) for image 1
    91     double normI2;                       ///< Sum(flux) for image 2
    92     double normSquare1;                 ///< Sum(flux^2) for image 1 (used for penalty)
    93     double normSquare2;                 ///< Sum(flux^2) for image 2 (used for penalty)
    94     pmSubtractionStampStatus status;    ///< Status of stamp
    95     psVector *MxxI1;                    ///< second moments of convolution images
    96     psVector *MyyI1;                    ///< second moments of convolution images
    97     psVector *MxxI2;                    ///< second moments of convolution images
    98     psVector *MyyI2;                    ///< second moments of convolution images
    99     double MxxI1raw;
    100     double MyyI1raw;
    101     double MxxI2raw;
    102     double MyyI2raw;
    103 } pmSubtractionStamp;
    104 
    10542/// Allocate a stamp
    10643pmSubtractionStamp *pmSubtractionStampAlloc(void);
     44
     45// find and extract the stamps
     46bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read
     47                               const pmReadout *ro1, // Readout 1
     48                               const pmReadout *ro2, // Readout 2
     49                               const psImage *subMask, // Mask for subtraction, or NULL
     50                               psImage *variance,  // Variance map
     51                               const psRegion *region, // Region of interest
     52                               float thresh1,  // Threshold for stamp finding on readout 1
     53                               float thresh2,  // Threshold for stamp finding on readout 2
     54                               float stampSpacing, // Spacing between stamps
     55                               float normFrac,     // Fraction of flux in window for normalisation window
     56                               float sysError,     // Relative systematic error in images
     57                               float skyError,     // Relative systematic error in images
     58                               int size,         // Kernel half-size
     59                               int footprint,     // Convolution footprint for stamps
     60                               pmSubtractionMode mode // Mode for subtraction
     61    );
     62
    10763
    10864/// Find stamps on an image
    … …  
    172128/// Calculate the window and normalisation window from the stamps
    173129bool pmSubtractionStampsGetWindow(
     130    bool *tryAgain,                     ///< re-try with new stamp size?
    174131    pmSubtractionStampList *stamps,     ///< List of stamps
    175132    int kernelSize                      ///< Half-size of kernel

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionThreads.c

    r30287 r30322  
    66#include <pslib.h>
    77
     8#include "pmSubtractionTypes.h"
    89#include "pmSubtractionMatch.h"
    910#include "pmSubtractionEquation.h"

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c

    r30289 r30322  
    1616
    1717#include "pmKapaPlots.h"
     18#include "pmFPA.h"
     19#include "pmSubtractionTypes.h"
    1820#include "pmSubtraction.h"
    1921#include "pmSubtractionStamps.h"

  • branches/eam_branches/ipp-20101205/psModules/src/psmodules.h

    r29388 r30322  
    9898// the following headers are from psModule:imcombine
    9999#include <pmStack.h>
     100#include <pmSubtractionTypes.h>
    100101#include <pmStackReject.h>
    101102#include <pmSubtraction.h>
Note: See TracChangeset for help on using the changeset viewer.