IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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;
Note: See TracChangeset for help on using the changeset viewer.