IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:41:49 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/tap_branches
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/tap_branches

  • branches/tap_branches/psModules

  • branches/tap_branches/psModules/src/imcombine/pmSubtractionMatch.c

    r25464 r27838  
    2828static bool useFFT = true;              // Do convolutions using FFT
    2929
     30# define SEPARATE 0
     31# if (SEPARATE)
     32# define SUBMODE PM_SUBTRACTION_EQUATION_NORM
     33# else
     34# define SUBMODE PM_SUBTRACTION_EQUATION_ALL
     35# endif
    3036
    3137//#define TESTING
    … …  
    6470                                 const psImage *subMask, // Mask for subtraction, or NULL
    6571                                 psImage *variance,  // Variance map
    66                                  const psRegion *region, // Region of interest, or NULL
     72                                 const psRegion *region, // Region of interest
    6773                                 float thresh1,  // Threshold for stamp finding on readout 1
    6874                                 float thresh2,  // Threshold for stamp finding on readout 2
    6975                                 float stampSpacing, // Spacing between stamps
     76                                 float normFrac,     // Fraction of flux in window for normalisation window
     77                                 float sysError,     // Relative systematic error in images
     78                                 float skyError,     // Relative systematic error in images
    7079                                 int size,         // Kernel half-size
    7180                                 int footprint,     // Convolution footprint for stamps
    … …  
    7382    )
    7483{
     84    PS_ASSERT_PTR_NON_NULL(stamps, false);
     85    PM_ASSERT_READOUT_NON_NULL(ro1, false);
     86    PM_ASSERT_READOUT_NON_NULL(ro2, false);
     87    PS_ASSERT_IMAGE_NON_NULL(subMask, false);
     88    PS_ASSERT_IMAGE_NON_NULL(variance, false);
     89    PS_ASSERT_PTR_NON_NULL(region, false);
     90
    7591    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    7692
    … …  
    7894
    7995    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
    80                                       size, footprint, stampSpacing, mode);
     96                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
    8197    if (!*stamps) {
    8298        psError(psErrorCodeLast(), false, "Unable to find stamps.");
    … …  
    87103
    88104    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    89     if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) {
    90         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
     105    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
     106        psError(psErrorCodeLast(), false, "Unable to extract stamps.");
    91107        return false;
    92108    }
    … …  
    101117                                  const pmReadout *ro1, const pmReadout *ro2, // Input images
    102118                                  int stride, // Size for convolution patches
    103                                   float sysError, // Relative systematic error
     119                                  float normFrac,           // Fraction of window for normalisation window
     120                                  float sysError,           // Systematic error in images
     121                                  float skyError,           // Systematic error in images
     122                                  float kernelError, // Systematic error in kernel
     123                                  float covarFrac,   // Fraction for kernel truncation before covariance
    104124                                  psImageMaskType maskVal, // Value to mask for input
    105125                                  psImageMaskType maskBad, // Mask for output bad pixels
    … …  
    112132    if (subMode != PM_SUBTRACTION_MODE_2) {
    113133        PM_ASSERT_READOUT_NON_NULL(conv1, false);
     134        PM_ASSERT_READOUT_NON_NULL(ro1, false);
     135        PM_ASSERT_READOUT_IMAGE(ro1, false);
    114136        if (conv1->image) {
    115137            psFree(conv1->image);
    … …  
    127149    if (subMode != PM_SUBTRACTION_MODE_1) {
    128150        PM_ASSERT_READOUT_NON_NULL(conv2, false);
     151        PM_ASSERT_READOUT_NON_NULL(ro2, false);
     152        PM_ASSERT_READOUT_IMAGE(ro2, false);
    129153        if (conv2->image) {
    130154            psFree(conv2->image);
    … …  
    141165    }
    142166
    143     PM_ASSERT_READOUT_NON_NULL(ro1, false);
    144     PM_ASSERT_READOUT_NON_NULL(ro2, false);
    145     PM_ASSERT_READOUT_IMAGE(ro1, false);
    146     PM_ASSERT_READOUT_IMAGE(ro2, false);
    147     PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
     167    if (ro1 && ro2) {
     168        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
     169    }
    148170    PS_ASSERT_INT_NONNEGATIVE(stride, false);
     171    if (isfinite(normFrac)) {
     172        PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false);
     173        PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false);
     174    }
    149175    if (isfinite(sysError)) {
    150176        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
    151177        PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
    152178    }
     179    if (isfinite(sysError)) {
     180        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.0, false);
     181    }
     182    if (isfinite(kernelError)) {
     183        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false);
     184        PS_ASSERT_FLOAT_LESS_THAN(kernelError, 1.0, false);
     185    }
     186    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(covarFrac, 0.0, false);
     187    PS_ASSERT_FLOAT_LESS_THAN(covarFrac, 1.0, false);
    153188    // Don't care about maskVal
    154189    // Don't care about maskBad
    … …  
    164199}
    165200
     201
     202/// Allocate images, as required
     203static void subtractionMatchAlloc(pmReadout *conv1, pmReadout *conv2, // Output readouts
     204                                  const pmReadout *ro1, const pmReadout *ro2, // Input readouts
     205                                  const psImage *subMask,                     // Subtraction mask
     206                                  psImageMaskType maskBad,                    // Mask value for bad pixels
     207                                  pmSubtractionMode subMode,          // Subtraction mode
     208                                  int numCols, int numRows            // Size of image
     209    )
     210{
     211    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
     212        subMode == PM_SUBTRACTION_MODE_DUAL) {
     213        if (!conv1->image) {
     214            conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     215        }
     216        psImageInit(conv1->image, NAN);
     217        if (ro1->variance) {
     218            if (!conv1->variance) {
     219                conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     220            }
     221            psImageInit(conv1->variance, NAN);
     222        }
     223        if (subMask) {
     224            if (!conv1->mask) {
     225                conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     226            }
     227            psImageInit(conv1->mask, maskBad);
     228        }
     229    }
     230    if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
     231        subMode == PM_SUBTRACTION_MODE_DUAL) {
     232        if (!conv2->image) {
     233            conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     234        }
     235        psImageInit(conv2->image, NAN);
     236        if (ro2->variance) {
     237            if (!conv2->variance) {
     238                conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     239            }
     240            psImageInit(conv2->variance, NAN);
     241        }
     242        if (subMask) {
     243            if (!conv2->mask) {
     244                conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     245            }
     246            psImageInit(conv2->mask, maskBad);
     247        }
     248    }
     249
     250    return;
     251}
     252
     253
    166254static void subtractionAnalysisUpdate(pmReadout *conv1, pmReadout *conv2, // Convolved images
    167255                                      const psMetadata *analysis, // Analysis metadata
    … …  
    192280}
    193281
     282bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {
     283
     284    if (!readout) return true;
     285
     286    psImage *image = readout->image;
     287    psImage *mask  = readout->mask;
     288    psImage *variance = readout->variance;
     289    for (int y = 0; y < image->numRows; y++) {
     290        for (int x = 0; x < image->numCols; x++) {
     291            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;
     292            bool valid = false;
     293            valid = isfinite(image->data.F32[y][x]);
     294            if (variance) {
     295                valid &= isfinite(variance->data.F32[y][x]);
     296            }
     297            if (valid) continue;
     298            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;
     299        }
     300    }
     301
     302    return true;
     303}
    194304
    195305bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
    196                                psMetadata *analysis, int stride, float sysError,
     306                               psMetadata *analysis, int stride, float kernelError, float covarFrac,
    197307                               psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor,
    198308                               float poorFrac, float badFrac)
    … …  
    210320        while ((item = psMetadataGetAndIncrement(iter))) {
    211321            if (item->type != PS_DATA_UNKNOWN) {
    212                 psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unexpected type for kernel.");
     322                psError(PM_ERR_PROG, true, "Unexpected type for kernel.");
    213323                psFree(iter);
    214324                psFree(kernelList);
    … …  
    230340    }
    231341    if (psListLength(kernelList) == 0) {
    232         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find kernels");
     342        psError(PM_ERR_PROG, true, "Unable to find kernels");
    233343        psFree(kernelList);
    234344        return false;
    … …  
    245355        while ((item = psMetadataGetAndIncrement(iter))) {
    246356            if (item->type != PS_DATA_REGION) {
    247                 psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unexpected type for region.");
     357                psError(PM_ERR_PROG, true, "Unexpected type for region.");
    248358                psFree(iter);
    249359                psFree(kernels);
    … …  
    257367    }
    258368    if (regions->n != kernels->n) {
    259         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Differing number of kernels (%ld) and regions (%ld)",
     369        psError(PM_ERR_PROG, true, "Differing number of kernels (%ld) and regions (%ld)",
    260370                kernels->n, regions->n);
    261371        psFree(regions);
    … …  
    264374    }
    265375
    266     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,
    267                                poorFrac, badFrac, mode)) {
     376    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError, covarFrac,
     377                               maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
    268378        psFree(kernels);
    269379        psFree(regions);
    … …  
    271381    }
    272382
    273     psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
     383    int numCols, numRows;       // Size of image
     384    if (ro1) {
     385        numCols = ro1->image->numCols;
     386        numRows = ro1->image->numRows;
     387    } else if (ro2) {
     388        numCols = ro2->image->numCols;
     389        numRows = ro2->image->numRows;
     390    } else {
     391        psAbort("No input image provided.");
     392    }
     393
     394    pmSubtractionMaskInvalid(ro1, maskVal);
     395    pmSubtractionMaskInvalid(ro2, maskVal);
     396
     397    // General background subtraction, since this is done in pmSubtractionMatch
     398    {
     399        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     400        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for background
     401        if (ro1) {
     402            psStatsInit(bg);
     403            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
     404                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
     405                psFree(bg);
     406                psFree(rng);
     407                return false;
     408            }
     409            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
     410        }
     411        if (ro2) {
     412            psStatsInit(bg);
     413            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
     414                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
     415                psFree(bg);
     416                psFree(rng);
     417                return false;
     418            }
     419            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
     420        }
     421        psFree(bg);
     422        psFree(rng);
     423    }
     424
     425    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     426
     427    psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0,
    274428                                         badFrac, mode); // Subtraction mask
    275429    if (!subMask) {
    … …  
    283437    psMetadata *outHeader = psMetadataAlloc(); // Output header values
    284438
     439    subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, mode, numCols, numRows);
     440
    285441    psTrace("psModules.imcombine", 2, "Convolving...\n");
    286442    for (int i = 0; i < kernels->n; i++) {
    … …  
    288444        psRegion *region = regions->data[i]; // Region of interest
    289445
    290         if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region,
    291                                    ro1->image->numCols, ro1->image->numRows)) {
    292             psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
     446        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region, numCols, numRows)) {
     447            psError(psErrorCodeLast(), false, "Unable to generate QA data");
    293448            psFree(outAnalysis);
    294449            psFree(outHeader);
    … …  
    300455
    301456        if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
    302                                    sysError, region, kernel, true, useFFT)) {
    303             psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
     457                                   kernelError, covarFrac, region, kernel, true, useFFT)) {
     458            psError(psErrorCodeLast(), false, "Unable to convolve image.");
    304459            psFree(outAnalysis);
    305460            psFree(outHeader);
    … …  
    330485                        int inner, int ringsOrder, int binning, float penalty,
    331486                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    332                         int iter, float rej, float sysError, psImageMaskType maskVal, psImageMaskType maskBad,
     487                        int iter, float rej, float normFrac, float sysError, float skyError,
     488                        float kernelError, float covarFrac, psImageMaskType maskVal, psImageMaskType maskBad,
    333489                        psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
    334490{
    335     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,
    336                                poorFrac, badFrac, subMode)) {
     491    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError,
     492                               covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
    337493        return false;
    338494    }
     495
     496    // We need both inputs
     497    PM_ASSERT_READOUT_NON_NULL(ro1, false);
     498    PM_ASSERT_READOUT_NON_NULL(ro2, false);
    339499
    340500    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
    … …  
    392552    // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
    393553    psImage *subMask = NULL;            // Mask for subtraction
    394     psRegion *region = NULL;            // Iso-kernel region
     554    psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region
    395555    psString regionString = NULL;       // String for region
    396556    pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
    … …  
    405565    memCheck("start");
    406566
    407     subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
    408                                 badFrac, subMode);
     567    pmSubtractionMaskInvalid(ro1, maskVal);
     568    pmSubtractionMaskInvalid(ro2, maskVal);
     569
     570    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     571
     572    subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode);
    409573    if (!subMask) {
    410574        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    … …  
    416580    // Get region of interest
    417581    int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
    418     float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
     582    float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions
    419583    if (isfinite(regionSize) && regionSize != 0.0) {
    420         xRegions = numCols / regionSize + 1;
    421         yRegions = numRows / regionSize + 1;
    422         xRegionSize = (float)numCols / (float)xRegions;
    423         yRegionSize = (float)numRows / (float)yRegions;
    424         region = psRegionAlloc(NAN, NAN, NAN, NAN);
    425     }
    426 
     584        xRegions = (bounds.x1 - bounds.x0) / regionSize + 1;
     585        yRegions = (bounds.y1 - bounds.y0) / regionSize + 1;
     586        xRegionSize = (float)(bounds.x1 - bounds.x0) / (float)xRegions;
     587        yRegionSize = (float)(bounds.y1 - bounds.y0) / (float)yRegions;
     588    } else {
     589        xRegionSize = bounds.x1 - bounds.x0;
     590        yRegionSize = bounds.y1 - bounds.y0;
     591    }
     592
     593    // General background subtraction and measurement of stamp threshold
    427594    float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
    428595    {
    429         psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
     596        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    430597        if (ro1) {
     598            psStatsInit(bg);
    431599            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
    432                 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
     600                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
    433601                psFree(bg);
    434602                goto MATCH_ERROR;
    435603            }
    436             stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
     604            stampThresh1 = threshold * bg->robustStdev;
     605            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
    437606        }
    438607        if (ro2) {
     608            psStatsInit(bg);
    439609            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
    440                 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
     610                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
    441611                psFree(bg);
    442612                goto MATCH_ERROR;
    443613            }
    444             stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
     614            stampThresh2 = threshold * bg->robustStdev;
     615            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
    445616        }
    446617        psFree(bg);
    447618    }
     619
     620    subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, subMode, numCols, numRows);
    448621
    449622    // Iterate over iso-kernel regions
    … …  
    452625            psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
    453626                    j * xRegions + i + 1, xRegions * yRegions);
    454             if (region) {
    455                 *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
    456                                       (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
    457                 psFree(regionString);
    458                 regionString = psRegionToString(*region);
    459                 psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
    460                         regionString, numCols, numRows);
    461             }
     627            *region = psRegionSet(bounds.x0 + (int)(i * xRegionSize),
     628                                  bounds.x0 + (int)((i + 1) * xRegionSize),
     629                                  bounds.y0 + (int)(j * yRegionSize),
     630                                  bounds.y0 + (int)((j + 1) * yRegionSize));
     631            psFree(regionString);
     632            regionString = psRegionToString(*region);
     633            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
     634                    regionString, numCols, numRows);
    462635
    463636            if (stampsName && strlen(stampsName) > 0) {
    464637                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
    465                                                         footprint, stampSpacing, subMode);
     638                                                        footprint, stampSpacing, normFrac,
     639                                                        sysError, skyError, subMode);
    466640            } else if (sources) {
    467641                stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
    468                                                            footprint, stampSpacing, subMode);
     642                                                           footprint, stampSpacing, normFrac,
     643                                                           sysError, skyError, subMode);
    469644            }
    470645
    471646            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    472647            // doesn't matter.
    473             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
    474                            stampSpacing, size, footprint, subMode)) {
     648            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
     649                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
    475650                goto MATCH_ERROR;
    476651            }
    477652
     653
     654            // generate the window function from the set of stamps
     655            if (!pmSubtractionStampsGetWindow(stamps, size)) {
     656                psError(psErrorCodeLast(), false, "Unable to get stamp window.");
     657                goto MATCH_ERROR;
     658            }
    478659
    479660            // Define kernel basis functions
    480661            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    481                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    482                                                           stamps, footprint, optThreshold, penalty, subMode);
     662                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
     663                                                          optFWHMs, optOrder, stamps, footprint,
     664                                                          optThreshold, penalty, bounds, subMode);
    483665                if (!kernels) {
    484666                    psErrorClear();
    … …  
    489671                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    490672                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    491                                                        inner, binning, ringsOrder, penalty, subMode);
     673                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
     674                // pmSubtractionVisualShowKernels(kernels);
    492675            }
    493676
    … …  
    500683                psVector *buffer = NULL;// Buffer for stats
    501684                if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) {
    502                     psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 1.");
     685                    psError(PM_ERR_DATA, false, "Unable to measure background of image 1.");
    503686                    psFree(bgStats);
    504687                    psFree(buffer);
    … …  
    507690                float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1
    508691                if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) {
    509                     psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 2.");
     692                    psError(PM_ERR_DATA, false, "Unable to measure background of image 2.");
    510693                    psFree(bgStats);
    511694                    psFree(buffer);
    … …  
    527710                    break;
    528711                  default:
    529                     psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order.");
     712                    psError(psErrorCodeLast(), false, "Unable to determine subtraction order.");
    530713                    goto MATCH_ERROR;
    531714                }
    … …  
    538721
    539722                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    540                                           stampThresh1, stampThresh2, stampSpacing,
    541                                           size, footprint, subMode)) {
    542                     goto MATCH_ERROR;
    543                 }
    544 
    545                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
    546                 if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    547                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    548                     goto MATCH_ERROR;
    549                 }
    550 
     723                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
     724                                          sysError, skyError, size, footprint, subMode)) {
     725                    goto MATCH_ERROR;
     726                }
     727
     728                // generate the window function from the set of stamps
     729                if (!pmSubtractionStampsGetWindow(stamps, size)) {
     730                    psError(psErrorCodeLast(), false, "Unable to get stamps window.");
     731                    goto MATCH_ERROR;
     732                }
     733
     734                // XXX step 1: calculate normalization
     735                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
     736                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
     737                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     738                    goto MATCH_ERROR;
     739                }
     740
     741                psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n");
     742                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
     743                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     744                    goto MATCH_ERROR;
     745                }
     746                memCheck("  solve equation");
     747
     748# if (SEPARATE)
     749                // set USED -> CALCULATE
     750                pmSubtractionStampsResetStatus (stamps);
     751
     752                // XXX step 2: calculate kernel parameters
     753                psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");
     754                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     755                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     756                    goto MATCH_ERROR;
     757                }
    551758                memCheck("  calculate equation");
    552759
    553                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
    554 
    555                 if (!pmSubtractionSolveEquation(kernels, stamps)) {
    556                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    557                     goto MATCH_ERROR;
    558                 }
    559 
     760                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
     761                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
     762                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     763                    goto MATCH_ERROR;
     764                }
    560765                memCheck("  solve equation");
    561 
     766# endif
    562767                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    563768                if (!deviations) {
    564                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
     769                    psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    565770                    goto MATCH_ERROR;
    566771                }
    … …  
    571776                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    572777                if (numRejected < 0) {
    573                     psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
     778                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    574779                    psFree(deviations);
    575780                    goto MATCH_ERROR;
    … …  
    580785            }
    581786
     787            // if we hit the max number of iterations and we have rejected stamps, re-solve
    582788            if (numRejected > 0) {
    583                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
    584                 if (!pmSubtractionSolveEquation(kernels, stamps)) {
    585                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    586                     goto MATCH_ERROR;
    587                 }
     789                // XXX step 1: calculate normalization
     790                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
     791                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
     792                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     793                    goto MATCH_ERROR;
     794                }
     795
     796                // solve normalization
     797                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
     798                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
     799                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     800                    goto MATCH_ERROR;
     801                }
     802
     803# if (SEPARATE)
     804                // set USED -> CALCULATE
     805                pmSubtractionStampsResetStatus (stamps);
     806
     807                // XXX step 2: calculate kernel parameters
     808                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
     809                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     810                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     811                    goto MATCH_ERROR;
     812                }
     813
     814                // solve kernel parameters
     815                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
     816                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
     817                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     818                    goto MATCH_ERROR;
     819                }
     820                memCheck("  solve equation");
     821# endif
    588822                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    589823                if (!deviations) {
    590                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
     824                    psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    591825                    goto MATCH_ERROR;
    592826                }
    … …  
    600834
    601835            if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
    602                 psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
     836                psError(psErrorCodeLast(), false, "Unable to generate QA data");
    603837                goto MATCH_ERROR;
    604838            }
    … …  
    608842            psTrace("psModules.imcombine", 2, "Convolving...\n");
    609843            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
    610                                        sysError, region, kernels, true, useFFT)) {
    611                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
     844                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
     845                psError(psErrorCodeLast(), false, "Unable to convolve image.");
    612846                goto MATCH_ERROR;
    613847            }
    … …  
    629863
    630864    if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
    631         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
     865        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
    632866        goto MATCH_ERROR;
    633867    }
    634868    if (conv2 && !pmSubtractionBorder(conv2->image, conv2->variance, conv2->mask, size, maskBad)) {
    635         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
     869        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
    636870        goto MATCH_ERROR;
    637871    }
    … …  
    8321066        } else {
    8331067            if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
    834                 psError(PS_ERR_UNKNOWN, false, "Unable to measure PSF width for stamp %d", i);
     1068                psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
    8351069                psFree(models);
    8361070                psFree(modelSums);
    … …  
    8431077
    8441078    if (!psThreadPoolWait(true)) {
    845         psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
     1079        psError(psErrorCodeLast(), false, "Error waiting for threads.");
    8461080        psFree(models);
    8471081        psFree(modelSums);
    … …  
    8561090    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    8571091    if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) {
    858         psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio.");
     1092        psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
    8591093        psFree(mask);
    8601094        psFree(ratios);
    … …  
    8901124    assert(kernels);
    8911125
    892     psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description);
    893     if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    894         psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1126    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1127    if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
     1128        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    8951129        return false;
    8961130    }
    8971131
    898     psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description);
    899     if (!pmSubtractionSolveEquation(kernels, stamps)) {
    900         psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1132    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
     1133    if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
     1134        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    9011135        return false;
    9021136    }
     1137
     1138# if (SEPARATE)
     1139    // set USED -> CALCULATE
     1140    pmSubtractionStampsResetStatus (stamps);
     1141
     1142    psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description);
     1143    if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1144        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1145        return false;
     1146    }
     1147
     1148    psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description);
     1149    if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1150        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1151        return false;
     1152    }
     1153# endif
    9031154
    9041155    psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
    9051156    psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    9061157    if (!deviations) {
    907         psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
     1158        psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    9081159        return false;
    9091160    }
    … …  
    9121163    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    9131164    if (numRejected < 0) {
    914         psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
     1165        psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    9151166        psFree(deviations);
    9161167        return false;
    … …  
    9201171    if (numRejected > 0) {
    9211172        // Allow re-fit with reduced stamps set
     1173        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1174        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
     1175            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1176            return false;
     1177        }
     1178
    9221179        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    923         if (!pmSubtractionSolveEquation(kernels, stamps)) {
    924             psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1180        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
     1181            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    9251182            return false;
    9261183        }
    9271184        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
     1185
     1186# if (SEPARATE)
     1187        // set USED -> CALCULATE
     1188        pmSubtractionStampsResetStatus (stamps);
     1189
     1190        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1191        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1192            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1193            return false;
     1194        }
     1195
     1196        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
     1197        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1198            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     1199            return false;
     1200        }
     1201        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
     1202# endif
     1203
    9281204        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    9291205        if (!deviations) {
    930             psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
     1206            psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    9311207            return false;
    9321208        }
    … …  
    9341210        long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
    9351211        if (numRejected < 0) {
    936             psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
     1212            psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    9371213            psFree(deviations);
    9381214            return false;
    … …  
    9581234
    9591235    if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) {
    960         psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 1");
     1236        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1");
    9611237        psFree(stamps1);
    9621238        psFree(kernels1);
    … …  
    9731249
    9741250    if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) {
    975         psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 2");
     1251        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2");
    9761252        psFree(stamps2);
    9771253        psFree(kernels2);
    … …  
    10121288}
    10131289
     1290
     1291bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
     1292                              float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax)
     1293{
     1294    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
     1295    PS_ASSERT_PTR_NON_NULL(stampSize, false);
     1296    PS_ASSERT_VECTOR_NON_NULL(widths, false);
     1297    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
     1298    PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);
     1299    PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);
     1300    PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
     1301    PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
     1302    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
     1303    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
     1304
     1305//    float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
     1306    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
     1307
     1308    if (isfinite(scaleMin) && scale < scaleMin) {
     1309        scale = scaleMin;
     1310    }
     1311    if (isfinite(scaleMax) && scale > scaleMax) {
     1312        scale = scaleMax;
     1313    }
     1314
     1315    for (int i = 0; i < widths->n; i++) {
     1316        widths->data.F32[i] *= scale;
     1317    }
     1318    *kernelSize = *kernelSize * scale + 0.5;
     1319    *stampSize = *stampSize * scale + 0.5;
     1320
     1321    psLogMsg("psModules.imcombine", PS_LOG_INFO,
     1322             "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
     1323
     1324    return true;
     1325}
Note: See TracChangeset for help on using the changeset viewer.