IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 13, 2011, 12:19:53 PM (15 years ago)
Author:
eugene
Message:

some code reorganzation (move all types into pmSubtractionTypes.h); add stage to check on different mode and order options, choosing best based on chisq of subtraction; careful with the window and the kernel sizes: if the measured kron radius is too large, allow the window to grow as needed; use the 1st radial moment measured on the stamps to set the kernel scaling; remove some cruft from the code; calling program now passes in kernel scaling options to be used when 1st radial moment is measured; new function to update the kernel description string used for kernel I/O; update the kernel description after the spatial order and direction is chosen; need to carry the kernel fwhms and spatial order to allow for update; calculate the psf-match solution errors; psf solution now uses the weights to get sensible chisq values, window is used to calculate the moments; penalty scale is now adjusted to be consistent with sensible reduced chisq; improved visualizations; use range-limiter in SVD of 1e-10; calculate chisq and moments for the solution and use in the evaluation; split out the stamp selection / convolution steps; reject sources after fitting a model to the chisq assuming non-zero systematic error

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtractionMatch.c

    r29596 r30622  
    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"
    … …  
    2728
    2829static bool useFFT = true;              // Do convolutions using FFT
    29 
    30 # define SUBMODE PM_SUBTRACTION_EQUATION_ALL
    3130
    3231//#define TESTING
    … …  
    5958}
    6059
    61 
    62 static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read
    63                                  const pmReadout *ro1, // Readout 1
    64                                  const pmReadout *ro2, // Readout 2
    65                                  const psImage *subMask, // Mask for subtraction, or NULL
    66                                  psImage *variance,  // Variance map
    67                                  const psRegion *region, // Region of interest
    68                                  float thresh1,  // Threshold for stamp finding on readout 1
    69                                  float thresh2,  // Threshold for stamp finding on readout 2
    70                                  float stampSpacing, // Spacing between stamps
    71                                  float normFrac,     // Fraction of flux in window for normalisation window
    72                                  float sysError,     // Relative systematic error in images
    73                                  float skyError,     // Relative systematic error in images
    74                                  int size,         // Kernel half-size
    75                                  int footprint,     // Convolution footprint for stamps
    76                                  pmSubtractionMode mode // Mode for subtraction
    77     )
    78 {
    79     PS_ASSERT_PTR_NON_NULL(stamps, false);
    80     PM_ASSERT_READOUT_NON_NULL(ro1, false);
    81     PM_ASSERT_READOUT_NON_NULL(ro2, false);
    82     PS_ASSERT_IMAGE_NON_NULL(subMask, false);
    83     PS_ASSERT_IMAGE_NON_NULL(variance, false);
    84     PS_ASSERT_PTR_NON_NULL(region, false);
    85 
    86     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    87 
    88     psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
    89 
    90     *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
    91                                       size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
    92     if (!*stamps) {
    93         psError(psErrorCodeLast(), false, "Unable to find stamps.");
    94         return false;
    95     }
    96 
    97     memCheck("  find stamps");
    98 
    99     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    100     if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
    101         psError(psErrorCodeLast(), false, "Unable to extract stamps.");
    102         return false;
    103     }
    104 
    105     memCheck("   extract stamps");
    106     pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1);
    107     return true;
    108 }
    109 
    11060// Check input arguments
    11161static bool subtractionMatchCheck(pmReadout *conv1, pmReadout *conv2, // Convolved images
    … …  
    12373                                  float badFrac,   // Maximum fraction of bad input pixels to accept
    12474                                  pmSubtractionMode subMode // Mode of subtraction
    125                                   )
     75    )
    12676{
    12777    if (subMode != PM_SUBTRACTION_MODE_2) {
    … …  
    473423}
    474424
     425bool pmSubtractionMatchAttempt(pmSubtractionQuality **bestMatch, pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, pmSubtractionMode mode, int spatialOrder, bool final) {
     426
     427    pmSubtractionMode nativeMode = kernels->mode;
     428    pmSubtractionMode nativeOrder = kernels->spatialOrder;
     429
     430    kernels->mode = mode;
     431    kernels->spatialOrder = spatialOrder;
     432
     433    // we always need to recalculate the matrix equation elements...
     434    pmSubtractionStampsResetStatus(stamps);
     435
     436    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
     437    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
     438        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
     439        return false;
     440    }
     441
     442    // step 1: generate the elements of the matrix equation Ax = B
     443    psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
     444    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     445        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     446        return false;
     447    }
     448                   
     449    // step 2: solve the matrix equation Ax = B
     450    psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
     451    if (!pmSubtractionSolveEquation(kernels, stamps)) {
     452        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     453        return false;
     454    }
     455    memCheck("  solve equation");
     456
     457    // calculate the score for this model fit attempt
     458    // XXX store the chisq, flux and moments for stamp rejection
     459    pmSubtractionCalculateChisqAndMoments(bestMatch, stamps, kernels); // Stamp deviations
     460
     461    // display the input and model stamps
     462    pmSubtractionVisualShowFit(stamps, kernels);
     463    pmSubtractionVisualPlotFit(kernels);
     464    pmSubtractionVisualPlotConvKernels(kernels);
     465
     466    // reset the kernel if desired (on final pass, do not reset)
     467    if (!final) {
     468        kernels->mode = nativeMode;
     469        kernels->spatialOrder = nativeOrder;
     470    } else {
     471      pmSubtractionKernelsMakeDescription(kernels);
     472      psLogMsg("psModules.imcombine", PS_LOG_INFO, "final kernel: %s", kernels->description);
     473    }
     474    return true;
     475}
    475476
    476477bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
    … …  
    478479                        const psArray *sources, const char *stampsName,
    479480                        pmSubtractionKernelsType type, int size, int spatialOrder,
    480                         const psVector *isisWidths, const psVector *isisOrders,
     481                        psVector *isisWidths, const psVector *isisOrders,
    481482                        int inner, int ringsOrder, int binning, float penalty,
    482483                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    … …  
    559560    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    560561
     562    pmSubtractionQuality *bestMatch = NULL;
     563
     564    int N_TEST_MODES;
     565    int N_TEST_ORDER = spatialOrder;
     566
     567    pmSubtractionMode TestModes[3];
     568    switch (subMode) {
     569      case PM_SUBTRACTION_MODE_1:
     570        N_TEST_MODES = 1;
     571        TestModes[0] = PM_SUBTRACTION_MODE_1;
     572        break;
     573      case PM_SUBTRACTION_MODE_2:
     574        N_TEST_MODES = 1;
     575        TestModes[0] = PM_SUBTRACTION_MODE_2;
     576        break;
     577      case PM_SUBTRACTION_MODE_DUAL:
     578        N_TEST_MODES = 3;
     579        TestModes[0] = PM_SUBTRACTION_MODE_1;
     580        TestModes[1] = PM_SUBTRACTION_MODE_2;
     581        TestModes[2] = PM_SUBTRACTION_MODE_DUAL;
     582        break;
     583      default:
     584        psError(psErrorCodeLast(), false, "For now, only modes 1, 2, and DUAL are supported.");
     585        goto MATCH_ERROR;
     586    }
     587   
    561588    memCheck("start");
    562589
    … …  
    628655            regionString = psRegionToString(*region);
    629656            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
    630                     regionString, numCols, numRows);
     657                     regionString, numCols, numRows);
    631658
    632659            if (stampsName && strlen(stampsName) > 0) {
    … …  
    640667            }
    641668
    642             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    643             // doesn't matter.
    644             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
    645                                       stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
    646                 goto MATCH_ERROR;
    647             }
    648 
    649 
    650             // generate the window function from the set of stamps
    651             if (!pmSubtractionStampsGetWindow(stamps, size)) {
    652                 psError(psErrorCodeLast(), false, "Unable to get stamp window.");
    653                 goto MATCH_ERROR;
    654             }
    655 
    656             // Define kernel basis functions
    657             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    658                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
    659                                                           optFWHMs, optOrder, stamps, footprint,
    660                                                           optThreshold, penalty, bounds, subMode);
    661                 if (!kernels) {
    662                     psErrorClear();
    663                     psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
    664                 }
    665             }
    666             if (kernels == NULL) {
    667                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    668                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    669                                                        inner, binning, ringsOrder, penalty, bounds, subMode);
    670             }
    671 
    672             memCheck("kernels");
    673 
    674             if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
    675                 pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
    676                 switch (newMode) {
    677                   case PM_SUBTRACTION_MODE_1:
    678                     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
    679                     break;
    680                   case PM_SUBTRACTION_MODE_2:
    681                     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
    682                     break;
    683                   default:
    684                     psError(psErrorCodeLast(), false, "Unable to determine subtraction order.");
    685                     goto MATCH_ERROR;
    686                 }
    687                 subMode = newMode;
    688             }
    689 
    690             int numRejected = -1;       // Number of rejected stamps in each iteration
    691             for (int k = 0; k < iter && numRejected != 0; k++) {
    692                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    693 
    694                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    695                                           stampThresh1, stampThresh2, stampSpacing, normFrac,
    696                                           sysError, skyError, size, footprint, subMode)) {
    697                     goto MATCH_ERROR;
    698                 }
    699 
    700                 // generate the window function from the set of stamps
    701                 if (!pmSubtractionStampsGetWindow(stamps, size)) {
    702                     psError(psErrorCodeLast(), false, "Unable to get stamps window.");
    703                     goto MATCH_ERROR;
    704                 }
     669            bool tryAgain = true;
     670            while (tryAgain) {
     671                // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
     672                // doesn't matter.
     673                if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
     674                                               stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
     675                    goto MATCH_ERROR;
     676                }
     677
     678                // generate the window function from the set of stamps
     679                if (!pmSubtractionStampsGetWindow(&tryAgain, stamps, size)) {
     680                    // if we failed, it might be due to the desired normWindow being larger than the current footprint.
     681                    // in this case, just adjust the footprint and try again.
     682                    if (tryAgain) {
     683                        // keep the border constant
     684                        int border = footprint - size;
     685                        size = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2;
     686                        footprint = size + border;
     687
     688                        // we need to reconstruct everything, so just free the stamps here and retry
     689                        psFree(stamps);
     690                    } else {
     691                        // unrecoverable error
     692                        psError(psErrorCodeLast(), false, "Unable to get stamp window.");
     693                        goto MATCH_ERROR;
     694                    }
     695                }
     696            }
     697
     698            // check on the kernel scaling -- if the kron-based radial moments are very different, adjust to match them
     699            {
     700                // float fwhm1;
     701                // float fwhm2;
     702                // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     703                // psAssert(isfinite(fwhm1), "fwhm 1 not set");
     704                // psAssert(isfinite(fwhm2), "fwhm 2 not set");
     705
     706                // XXX this is BAD: depends on the relationship below:
     707                // stamps->normWindow1 = 2.75*R1;
     708                // stamps->normWindow2 = 2.75*R2;
     709                float radMoment1 = stamps->normWindow1 / 2.75;
     710                float radMoment2 = stamps->normWindow2 / 2.75;
     711                pmSubtractionParamsScale(NULL, NULL, isisWidths, radMoment1, radMoment2);
     712
     713                // float maxFWHM = PS_MAX(fwhm1, fwhm2);
     714                // float maxRadial = PS_MAX(radMoment1, radMoment2);
     715               
     716                // if (fabs(2.0*(maxFWHM - maxRadial)/(maxFWHM + maxRadial)) > 0.25) {
     717                // if (1) {
     718                //
     719                //     float scale = maxRadial / maxFWHM;
     720                //     psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale);
     721                //     
     722                //     for (int i = 0; i < isisWidths->n; i++) {
     723                //      isisWidths->data.F32[i] *= scale;
     724                //     }
     725                // }
     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                bool tryAgain = true;
     767                while (tryAgain) {
     768                    // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
     769                    // doesn't matter.
     770                    if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
     771                                                   stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
     772                        goto MATCH_ERROR;
     773                    }
     774
     775                    // generate the window function from the set of stamps
     776                    if (!pmSubtractionStampsGetWindow(&tryAgain, stamps, size)) {
     777                        // if we failed, it might be due to the desired normWindow being larger than the current footprint.
     778                        // in this case, just adjust the footprint and try again.
     779                        if (tryAgain) {
     780                            footprint = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2;
     781
     782                            // we need to reconstruct everything, so just free the stamps here and retry
     783                            psFree(stamps);
     784                        } else {
     785                            // unrecoverable error
     786                            psError(psErrorCodeLast(), false, "Unable to get stamp window.");
     787                            goto MATCH_ERROR;
     788                        }
     789                    }
     790                }
    705791
    706792                // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue
    707                 psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
    708                 if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) {
    709                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    710                     goto MATCH_ERROR;
    711                 }
    712 
    713                 // step 0b : calculate the moments, pass along to the next steps via stamps->normValue
    714                 // XXX this step is now skipped -- delete
    715                 psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
    716                 if (!pmSubtractionCalculateMoments(kernels, stamps)) {
    717                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    718                     goto MATCH_ERROR;
    719                 }
    720 
    721                 // step 1: generate the elements of the matrix equation Ax = B
    722                 psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
    723                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    724                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    725                     goto MATCH_ERROR;
    726                 }
    727 
    728                 // step 2: solve the matrix equation Ax = B
    729                 psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
    730                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    731                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    732                     goto MATCH_ERROR;
    733                 }
    734                 memCheck("  solve equation");
    735 
    736                 pmSubtractionVisualPlotConvKernels(kernels);
    737 
    738                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    739                 if (!deviations) {
    740                     psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    741                     goto MATCH_ERROR;
    742                 }
    743                 memCheck("   calculate deviations");
    744 
    745                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
    746                 numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    747                 if (numRejected < 0) {
    748                     psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    749                     psFree(deviations);
    750                     goto MATCH_ERROR;
    751                 }
    752                 psFree(deviations);
    753 
    754                 memCheck("  reject stamps");
    755             }
    756 
    757             // if we hit the max number of iterations and we have rejected stamps, re-solve
    758             if (numRejected > 0) {
    759 
    760                 // step 1: generate the elements of the matrix equation Ax = B
    761                 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
    762                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    763                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    764                     goto MATCH_ERROR;
    765                 }
    766 
    767                 // step 2: solve the matrix equation Ax = B
    768                 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
    769                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    770                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    771                     goto MATCH_ERROR;
    772                 }
    773 
    774                 pmSubtractionVisualPlotConvKernels(kernels);
    775 
    776                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    777                 if (!deviations) {
    778                     psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    779                     goto MATCH_ERROR;
    780                 }
    781                 pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
    782                 psFree(deviations);
    783             }
    784             psFree(stamps);
    785             stamps = NULL;
    786 
    787             memCheck("solution");
    788 
    789             if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
    790                 psError(psErrorCodeLast(), false, "Unable to generate QA data");
    791                 goto MATCH_ERROR;
    792             }
    793             memCheck("diag outputs");
    794 
    795             psTrace("psModules.imcombine", 2, "Convolving...\n");
    796             if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
    797                                        kernelError, covarFrac, region, kernels, true, useFFT)) {
    798                 psError(psErrorCodeLast(), false, "Unable to convolve image.");
    799                 goto MATCH_ERROR;
    800             }
    801 
    802             psFree(kernels);
    803             kernels = NULL;
    804         }
     793                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
     794                if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) {
     795                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     796                    goto MATCH_ERROR;
     797                }
     798
     799                // on each iteration, we start from scratch
     800                psFree(bestMatch);
     801
     802                // choose the spatial order and subtraction direction (1, 2, dual)
     803                // XXX need to make these respect recipe somewhat
     804                for (int order = 0; order <= N_TEST_ORDER; order++) {
     805                    for (int j = 0; j < N_TEST_MODES; j++) {
     806                        if (!pmSubtractionMatchAttempt(&bestMatch, kernels, stamps, TestModes[j], order, false)) {
     807                            goto MATCH_ERROR;
     808                        }
     809                    }
     810                }
     811               
     812                // reject the deviant stamps based on the stats of the best match
     813                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
     814                numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej);
     815                if (numRejected < 0) {
     816                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
     817                    goto MATCH_ERROR;
     818                }
     819                memCheck("  reject stamps");
     820            }
     821
     822            // apply the best fit so we are ready to roll
     823            psLogMsg("psModules.imcombine", PS_LOG_INFO, "applying order: %d, mode: %d\n", bestMatch->spatialOrder, bestMatch->mode);
     824            if (!pmSubtractionMatchAttempt(NULL, kernels, stamps, bestMatch->mode, bestMatch->spatialOrder, true)) {
     825                goto MATCH_ERROR;
     826            }
     827            psFree(stamps);
     828            psFree(bestMatch);
     829            memCheck("solution");
     830
     831            if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
     832                psError(psErrorCodeLast(), false, "Unable to generate QA data");
     833                goto MATCH_ERROR;
     834            }
     835            memCheck("diag outputs");
     836
     837            psTrace("psModules.imcombine", 2, "Convolving...\n");
     838            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
     839                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
     840                psError(psErrorCodeLast(), false, "Unable to convolve image.");
     841                goto MATCH_ERROR;
     842            }
     843
     844            psFree(kernels);
     845            kernels = NULL;
     846        }
    805847    }
    806848    psFree(rng);
    … …  
    816858
    817859    if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
    818         psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
    819         goto MATCH_ERROR;
     860        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
     861        goto MATCH_ERROR;
    820862    }
    821863    if (conv2 && !pmSubtractionBorder(conv2->image, conv2->variance, conv2->mask, size, maskBad)) {
    822         psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
    823         goto MATCH_ERROR;
     864        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
     865        goto MATCH_ERROR;
    824866    }
    825867
    … …  
    832874#ifdef TESTING
    833875    {
    834         if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    835             psFits *fits = psFitsOpen("convolved1.fits", "w");
    836             psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
    837             psFitsClose(fits);
    838         }
    839 
    840         if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    841             psFits *fits = psFitsOpen("convolved2.fits", "w");
    842             psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);
    843             psFitsClose(fits);
    844         }
     876        if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
     877            psFits *fits = psFitsOpen("convolved1.fits", "w");
     878            psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
     879            psFitsClose(fits);
     880        }
     881
     882        if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
     883            psFits *fits = psFitsOpen("convolved2.fits", "w");
     884            psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);
     885            psFitsClose(fits);
     886        }
    845887    }
    846888#endif
    … …  
    858900    psFree(variance);
    859901    psFree(rng);
     902    psFree(bestMatch);
    860903    return false;
    861904}
    … …  
    866909// increment).
    867910static int subtractionOrderWidth(const psKernel *kernel, // Image
    868                                 float bg, // Background in image
    869                                 int size, // Maximum size
    870                                 const psArray *models, // Buffer of models
    871                                 const psVector *modelSums // Buffer of model sums
     911                                float bg, // Background in image
     912                                int size, // Maximum size
     913                                const psArray *models, // Buffer of models
     914                                const psVector *modelSums // Buffer of model sums
    872915    )
    873916{
    … …  
    882925    psVector *chi2 = psVectorAlloc(size, PS_TYPE_F32); // chi^2 as a function of radius
    883926    for (int sigma = 0; sigma < size; sigma++) {
    884         double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian
    885         psKernel *model = models->data[sigma]; // Model of interest
    886         for (int y = yMin; y <= yMax; y++) {
    887             for (int x = xMin; x <= xMax; x++) {
    888                 sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg);
    889             }
    890         }
    891         float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian
    892         double sumDev2 = 0.0;           // Sum of square deviations
    893         for (int y = yMin; y <= yMax; y++) {
    894             for (int x = xMin; x <= xMax; x++) {
    895                 float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation
    896                 sumDev2 += PS_SQR(dev);
    897             }
    898         }
    899         chi2->data.F32[sigma] = sumDev2;
     927        double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian
     928        psKernel *model = models->data[sigma]; // Model of interest
     929        for (int y = yMin; y <= yMax; y++) {
     930            for (int x = xMin; x <= xMax; x++) {
     931                sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg);
     932            }
     933        }
     934        float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian
     935        double sumDev2 = 0.0;           // Sum of square deviations
     936        for (int y = yMin; y <= yMax; y++) {
     937            for (int x = xMin; x <= xMax; x++) {
     938                float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation
     939                sumDev2 += PS_SQR(dev);
     940            }
     941        }
     942        chi2->data.F32[sigma] = sumDev2;
    900943    }
    901944
    … …  
    904947    float bestChi2 = INFINITY;          // Best chi^2
    905948    for (int i = 0; i < size; i++) {
    906         if (chi2->data.F32[i] < bestChi2) {
    907             bestChi2 = chi2->data.F32[i];
    908             bestIndex = i;
    909         }
     949        if (chi2->data.F32[i] < bestChi2) {
     950            bestChi2 = chi2->data.F32[i];
     951            bestIndex = i;
     952        }
    910953    }
    911954    psFree(chi2);
    … …  
    916959
    917960bool pmSubtractionOrderStamp(psVector *ratios, psVector *mask, const pmSubtractionStampList *stamps,
    918                              const psArray *models, const psVector *modelSums,
    919                              int index, float bg1, float bg2)
     961                             const psArray *models, const psVector *modelSums,
     962                             int index, float bg1, float bg2)
    920963{
    921964    PS_ASSERT_VECTOR_NON_NULL(ratios, false);
    … …  
    931974    pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
    932975    psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE || stamp->status == PM_SUBTRACTION_STAMP_USED,
    933              "We checked this earlier.");
     976             "We checked this earlier.");
    934977
    935978    // Widths of stars
    … …  
    938981
    939982    if (width1 == 0 || width2 == 0) {
    940         ratios->data.F32[index] = NAN;
    941         mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff;
     983        ratios->data.F32[index] = NAN;
     984        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff;
    942985    } else {
    943         ratios->data.F32[index] = (float)width1 / (float)width2;
    944         mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0;
    945         psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n",
    946                 index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]);
     986        ratios->data.F32[index] = (float)width1 / (float)width2;
     987        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0;
     988        psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n",
     989                index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]);
    947990    }
    948991
    … …  
    9781021    psVector *modelSums = psVectorAlloc(size, PS_TYPE_F64); // Gaussian model sums
    9791022    for (int sigma = 0; sigma < size; sigma++) {
    980         psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model
    981         float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared
    982         double sumGG = 0.0;         // Sum of square of Gaussian
    983         for (int y = -size; y <= size; y++) {
    984             int y2 = PS_SQR(y);     // y squared
    985             for (int x = -size; x <= size; x++) {
    986                 float rad2 = PS_SQR(x) + y2; // Radius squared
    987                 float value = expf(-rad2 * invSigma2); // Model value
    988                 model->kernel[y][x] = value;
    989                 sumGG += PS_SQR(value);
    990             }
    991         }
    992         models->data[sigma] = model;
    993         modelSums->data.F64[sigma] = 1.0 / sumGG;
     1023        psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model
     1024        float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared
     1025        double sumGG = 0.0;         // Sum of square of Gaussian
     1026        for (int y = -size; y <= size; y++) {
     1027            int y2 = PS_SQR(y);     // y squared
     1028            for (int x = -size; x <= size; x++) {
     1029                float rad2 = PS_SQR(x) + y2; // Radius squared
     1030                float value = expf(-rad2 * invSigma2); // Model value
     1031                model->kernel[y][x] = value;
     1032                sumGG += PS_SQR(value);
     1033            }
     1034        }
     1035        models->data[sigma] = model;
     1036        modelSums->data.F64[sigma] = 1.0 / sumGG;
    9941037    }
    9951038
    9961039    // Fit models to stamps
    9971040    for (int i = 0; i < stamps->num; i++) {
    998         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    999         if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
    1000             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
    1001             continue;
    1002         }
    1003 
    1004         if (pmSubtractionThreaded()) {
    1005             psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER");
    1006             psArrayAdd(job->args, 1, ratios);
    1007             psArrayAdd(job->args, 1, mask);
    1008             psArrayAdd(job->args, 1, stamps);
    1009             psArrayAdd(job->args, 1, models);
    1010             psArrayAdd(job->args, 1, modelSums);
    1011             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
    1012             PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32);
    1013             PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
    1014             if (!psThreadJobAddPending(job)) {
    1015                 return false;
    1016             }
    1017         } else {
    1018             if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
    1019                 psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
    1020                 psFree(models);
    1021                 psFree(modelSums);
    1022                 psFree(ratios);
    1023                 psFree(mask);
    1024                 return false;
    1025             }
    1026         }
     1041        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     1042        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
     1043            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
     1044            continue;
     1045        }
     1046
     1047        if (pmSubtractionThreaded()) {
     1048            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER");
     1049            psArrayAdd(job->args, 1, ratios);
     1050            psArrayAdd(job->args, 1, mask);
     1051            psArrayAdd(job->args, 1, stamps);
     1052            psArrayAdd(job->args, 1, models);
     1053            psArrayAdd(job->args, 1, modelSums);
     1054            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
     1055            PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32);
     1056            PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
     1057            if (!psThreadJobAddPending(job)) {
     1058                return false;
     1059            }
     1060        } else {
     1061            if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
     1062                psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
     1063                psFree(models);
     1064                psFree(modelSums);
     1065                psFree(ratios);
     1066                psFree(mask);
     1067                return false;
     1068            }
     1069        }
    10271070    }
    10281071
    10291072    if (!psThreadPoolWait(true)) {
    1030         psError(psErrorCodeLast(), false, "Error waiting for threads.");
    1031         psFree(models);
    1032         psFree(modelSums);
    1033         psFree(ratios);
    1034         psFree(mask);
    1035             return false;
     1073        psError(psErrorCodeLast(), false, "Error waiting for threads.");
     1074        psFree(models);
     1075        psFree(modelSums);
     1076        psFree(ratios);
     1077        psFree(mask);
     1078        return false;
    10361079    }
    10371080
    … …  
    10411084    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    10421085    if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) {
    1043         psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
    1044         psFree(mask);
    1045         psFree(ratios);
    1046         psFree(stats);
    1047         return PM_SUBTRACTION_MODE_ERR;
     1086        psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
     1087        psFree(mask);
     1088        psFree(ratios);
     1089        psFree(stats);
     1090        return PM_SUBTRACTION_MODE_ERR;
    10481091    }
    10491092    psFree(ratios);
    … …  
    10521095    // XXX raise an error here or not?
    10531096    if (isnan(stats->robustMedian)) {
    1054         psFree(stats);
    1055         return PM_SUBTRACTION_MODE_ERR;
     1097        psFree(stats);
     1098        return PM_SUBTRACTION_MODE_ERR;
    10561099    }
    10571100
    … …  
    10661109// Test a subtraction mode by performing a single iteration
    10671110static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode
    1068                                 pmSubtractionKernels *kernels, // Kernel description
    1069                                 const char *description, // Description for trace
    1070                                 psImage *subMask,  // Subtraction mask
    1071                                 float rej               // Rejection threshold
    1072                                 )
     1111                                pmSubtractionKernels *kernels, // Kernel description
     1112                                const char *description, // Description for trace
     1113                                psImage *subMask,  // Subtraction mask
     1114                                float rej               // Rejection threshold
     1115    )
    10731116{
    10741117    assert(stamps);
    10751118    assert(kernels);
    10761119
     1120    psAbort("this function is not working");
     1121# if (0)
     1122    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
     1123    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
     1124        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
     1125        return false;
     1126    }
     1127
    10771128    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1078     if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    1079         psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1080         return false;
     1129    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     1130        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1131        return false;
    10811132    }
    10821133
    10831134    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
    1084     if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    1085         psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1086         return false;
     1135    if (!pmSubtractionSolveEquation(kernels, stamps)) {
     1136        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1137        return false;
    10871138    }
    10881139
    … …  
    10901141    psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    10911142    if (!deviations) {
    1092         psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    1093         return false;
    1094     }
    1095 
     1143        psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
     1144        return false;
     1145    }
     1146
     1147    // XXX this needs to be made consistent with the modified 'reject stamps' function
    10961148    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
    10971149    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    10981150    if (numRejected < 0) {
    1099         psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    1100         psFree(deviations);
    1101         return false;
     1151        psError(psErrorCodeLast(), false, "Unable to reject stamps.");
     1152        psFree(deviations);
     1153        return false;
    11021154    }
    11031155    psFree(deviations);
    11041156
    11051157    if (numRejected > 0) {
    1106         // Allow re-fit with reduced stamps set
    1107         psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1108         if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
    1109             psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1110             return false;
    1111         }
    1112 
    1113         psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    1114         if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
    1115             psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    1116             return false;
    1117         }
    1118         psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
    1119 
    1120         psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    1121         if (!deviations) {
    1122             psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    1123             return false;
    1124         }
    1125         psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description);
    1126         long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
    1127         if (numRejected < 0) {
    1128             psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    1129             psFree(deviations);
    1130             return false;
    1131         }
    1132         psFree(deviations);
    1133     }
    1134 
     1158        // Allow re-fit with reduced stamps set
     1159        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1160        if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     1161            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1162            return false;
     1163        }
     1164
     1165        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
     1166        if (!pmSubtractionSolveEquation(kernels, stamps)) {
     1167            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1168            return false;
     1169        }
     1170        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
     1171
     1172        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
     1173        if (!deviations) {
     1174            psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
     1175            return false;
     1176        }
     1177        psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description);
     1178        long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
     1179        if (numRejected < 0) {
     1180            psError(psErrorCodeLast(), false, "Unable to reject stamps.");
     1181            psFree(deviations);
     1182            return false;
     1183        }
     1184        psFree(deviations);
     1185    }
     1186# endif
    11351187    return true;
    11361188}
    … …  
    11381190
    11391191pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels,
    1140                                         const psImage *subMask, float rej)
     1192                                        const psImage *subMask, float rej)
    11411193{
    11421194    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR);
    … …  
    11501202
    11511203    if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) {
    1152         psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1");
    1153         psFree(stamps1);
    1154         psFree(kernels1);
    1155         psFree(subMask1);
    1156         return PM_SUBTRACTION_MODE_ERR;
     1204        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1");
     1205        psFree(stamps1);
     1206        psFree(kernels1);
     1207        psFree(subMask1);
     1208        return PM_SUBTRACTION_MODE_ERR;
    11571209    }
    11581210    psFree(subMask1);
    … …  
    11651217
    11661218    if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) {
    1167         psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2");
    1168         psFree(stamps2);
    1169         psFree(kernels2);
    1170         psFree(subMask2);
    1171         psFree(stamps1);
    1172         psFree(kernels1);
    1173         return PM_SUBTRACTION_MODE_ERR;
     1219        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2");
     1220        psFree(stamps2);
     1221        psFree(kernels2);
     1222        psFree(subMask2);
     1223        psFree(stamps1);
     1224        psFree(kernels1);
     1225        return PM_SUBTRACTION_MODE_ERR;
    11741226    }
    11751227    psFree(subMask2);
    … …  
    11791231    pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels
    11801232    psLogMsg("psModules.imcombine", PS_LOG_INFO,
    1181              "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n",
    1182              kernels1->mean, kernels1->rms, kernels1->numStamps,
    1183              kernels2->mean, kernels2->rms, kernels2->numStamps);
     1233             "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n",
     1234             kernels1->mean, kernels1->rms, kernels1->numStamps,
     1235             kernels2->mean, kernels2->rms, kernels2->numStamps);
    11841236
    11851237    if (kernels1->mean < kernels2->mean) {
    1186         bestStamps = stamps1;
    1187         bestKernels = kernels1;
     1238        bestStamps = stamps1;
     1239        bestKernels = kernels1;
    11881240    } else {
    1189         bestStamps = stamps2;
    1190         bestKernels = kernels2;
     1241        bestStamps = stamps2;
     1242        bestKernels = kernels2;
    11911243    }
    11921244
    … …  
    12041256}
    12051257
    1206 
    1207 bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
    1208                               float scaleRef, float scaleMin, float scaleMax)
     1258static float scaleRefOption = NAN;
     1259static float scaleMinOption = NAN;
     1260static float scaleMaxOption = NAN;
     1261static bool  scaleOption = false;
     1262
     1263bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax) {
     1264
     1265    if (scale) {
     1266        PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
     1267        PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
     1268        PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
     1269        PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
     1270    }
     1271
     1272    scaleRefOption = scaleRef;
     1273    scaleMinOption = scaleMin;
     1274    scaleMaxOption = scaleMax;
     1275    scaleOption = scale;
     1276   
     1277    return true;
     1278}
     1279
     1280bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2)
    12091281{
    1210     PS_ASSERT_PTR_NON_NULL(kernelSize, false);
    1211     PS_ASSERT_PTR_NON_NULL(stampSize, false);
     1282    // PS_ASSERT_PTR_NON_NULL(kernelSize, false);
     1283    // PS_ASSERT_PTR_NON_NULL(stampSize, false);
    12121284    PS_ASSERT_VECTOR_NON_NULL(widths, false);
    12131285    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
    1214     PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
    1215     PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
    1216     PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
    1217     PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
    1218 
    1219     float fwhm1;
    1220     float fwhm2;
    1221 
    1222     pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
    1223     psAssert(isfinite(fwhm1), "fwhm 1 not set");
    1224     psAssert(isfinite(fwhm2), "fwhm 2 not set");
     1286
     1287    if (!scaleOption) return true;
     1288
     1289    // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     1290    // psAssert(isfinite(fwhm1), "fwhm 1 not set");
     1291    // psAssert(isfinite(fwhm2), "fwhm 2 not set");
    12251292   
    12261293    // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
    1227     float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
    1228 
    1229     if (isfinite(scaleMin) && scale < scaleMin) {
    1230         scale = scaleMin;
    1231     }
    1232     if (isfinite(scaleMax) && scale > scaleMax) {
    1233         scale = scaleMax;
     1294    float scale = PS_MAX(fwhm1, fwhm2) / scaleRefOption;      // Scaling factor
     1295
     1296    if (isfinite(scaleMinOption) && scale < scaleMinOption) {
     1297        scale = scaleMinOption;
     1298    }
     1299    if (isfinite(scaleMaxOption) && scale > scaleMaxOption) {
     1300        scale = scaleMaxOption;
    12341301    }
    12351302
    12361303    for (int i = 0; i < widths->n; i++) {
    1237         widths->data.F32[i] *= scale;
    1238     }
    1239     *kernelSize = *kernelSize * scale + 0.5;
    1240     *stampSize = *stampSize * scale + 0.5;
    1241 
    1242     psLogMsg("psModules.imcombine", PS_LOG_INFO,
    1243              "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
     1304        widths->data.F32[i] *= scale;
     1305    }
     1306    if (kernelSize) {
     1307        *kernelSize = *kernelSize * scale + 0.5;
     1308    }
     1309    if (stampSize) {
     1310        *stampSize = *stampSize * scale + 0.5;
     1311    }
     1312
     1313    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Scaling kernel parameters by %f", scale);
     1314    if (kernelSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified kernel size %d", *kernelSize);
     1315    if (stampSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified stamp size %d", *stampSize);
    12441316
    12451317    return true;
Note: See TracChangeset for help on using the changeset viewer.