IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 10, 2010, 7:34:39 PM (16 years ago)
Author:
eugene
Message:

updates from eam_branches/20091201 (substantially changes to the kernel matching code; updates to pmDetection container, added pmPattern, etc)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26035 r26893  
    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
    7077                                 float sysError,     // Relative systematic error in images
     78                                 float skyError,     // Relative systematic error in images
    7179                                 int size,         // Kernel half-size
    7280                                 int footprint,     // Convolution footprint for stamps
    … …  
    7482    )
    7583{
     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
    7691    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    7792
    … …  
    7994
    8095    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
    81                                       size, footprint, stampSpacing, sysError, mode);
     96                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
    8297    if (!*stamps) {
    8398        psError(psErrorCodeLast(), false, "Unable to find stamps.");
    … …  
    88103
    89104    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    90     if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) {
     105    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
    91106        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
    92107        return false;
    … …  
    102117                                  const pmReadout *ro1, const pmReadout *ro2, // Input images
    103118                                  int stride, // Size for convolution patches
     119                                  float normFrac,           // Fraction of window for normalisation window
    104120                                  float sysError,           // Systematic error in images
     121                                  float skyError,           // Systematic error in images
    105122                                  float kernelError, // Systematic error in kernel
     123                                  float covarFrac,   // Fraction for kernel truncation before covariance
    106124                                  psImageMaskType maskVal, // Value to mask for input
    107125                                  psImageMaskType maskBad, // Mask for output bad pixels
    … …  
    149167    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
    150168    PS_ASSERT_INT_NONNEGATIVE(stride, false);
     169    if (isfinite(normFrac)) {
     170        PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false);
     171        PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false);
     172    }
    151173    if (isfinite(sysError)) {
    152174        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
    153175        PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
    154176    }
     177    if (isfinite(sysError)) {
     178        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.0, false);
     179    }
    155180    if (isfinite(kernelError)) {
    156181        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false);
    157182        PS_ASSERT_FLOAT_LESS_THAN(kernelError, 1.0, false);
    158183    }
     184    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(covarFrac, 0.0, false);
     185    PS_ASSERT_FLOAT_LESS_THAN(covarFrac, 1.0, false);
    159186    // Don't care about maskVal
    160187    // Don't care about maskBad
    … …  
    198225}
    199226
     227bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {
     228
     229    if (!readout) return true;
     230
     231    psImage *image = readout->image;
     232    psImage *mask  = readout->mask;
     233    psImage *variance = readout->variance;
     234    for (int y = 0; y < image->numRows; y++) {
     235        for (int x = 0; x < image->numCols; x++) {
     236            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;
     237            bool valid = false;
     238            valid = isfinite(image->data.F32[y][x]);
     239            if (variance) {
     240                valid &= isfinite(variance->data.F32[y][x]);
     241            }
     242            if (valid) continue;
     243            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;
     244        }
     245    }
     246
     247    return true;
     248}
    200249
    201250bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
    202                                psMetadata *analysis, int stride, float kernelError,
     251                               psMetadata *analysis, int stride, float kernelError, float covarFrac,
    203252                               psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor,
    204253                               float poorFrac, float badFrac)
    … …  
    270319    }
    271320
    272     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, kernelError,
     321    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError, covarFrac,
    273322                               maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
    274323        psFree(kernels);
    … …  
    277326    }
    278327
    279     psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
     328    pmSubtractionMaskInvalid(ro1, maskVal);
     329    pmSubtractionMaskInvalid(ro2, maskVal);
     330
     331    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     332
     333    psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0,
    280334                                         badFrac, mode); // Subtraction mask
    281335    if (!subMask) {
    … …  
    306360
    307361        if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
    308                                    kernelError, region, kernel, true, useFFT)) {
     362                                   kernelError, covarFrac, region, kernel, true, useFFT)) {
    309363            psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
    310364            psFree(outAnalysis);
    … …  
    336390                        int inner, int ringsOrder, int binning, float penalty,
    337391                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    338                         int iter, float rej, float sysError, float kernelError, psImageMaskType maskVal,
    339                         psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac,
    340                         float badFrac, pmSubtractionMode subMode)
     392                        int iter, float rej, float normFrac, float sysError, float skyError,
     393                        float kernelError, float covarFrac, psImageMaskType maskVal, psImageMaskType maskBad,
     394                        psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
    341395{
    342     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, kernelError,
    343                                maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
     396    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError,
     397                               covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
    344398        return false;
    345399    }
    … …  
    399453    // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
    400454    psImage *subMask = NULL;            // Mask for subtraction
    401     psRegion *region = NULL;            // Iso-kernel region
     455    psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region
    402456    psString regionString = NULL;       // String for region
    403457    pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
    … …  
    412466    memCheck("start");
    413467
    414     subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
    415                                 badFrac, subMode);
     468    pmSubtractionMaskInvalid(ro1, maskVal);
     469    pmSubtractionMaskInvalid(ro2, maskVal);
     470
     471    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     472
     473    subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode);
    416474    if (!subMask) {
    417475        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    … …  
    423481    // Get region of interest
    424482    int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
    425     float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
     483    float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions
    426484    if (isfinite(regionSize) && regionSize != 0.0) {
    427         xRegions = numCols / regionSize + 1;
    428         yRegions = numRows / regionSize + 1;
    429         xRegionSize = (float)numCols / (float)xRegions;
    430         yRegionSize = (float)numRows / (float)yRegions;
    431         region = psRegionAlloc(NAN, NAN, NAN, NAN);
    432     }
    433 
     485        xRegions = (bounds.x1 - bounds.x0) / regionSize + 1;
     486        yRegions = (bounds.y1 - bounds.y0) / regionSize + 1;
     487        xRegionSize = (float)(bounds.x1 - bounds.x0) / (float)xRegions;
     488        yRegionSize = (float)(bounds.y1 - bounds.y0) / (float)yRegions;
     489    } else {
     490        xRegionSize = bounds.x1 - bounds.x0;
     491        yRegionSize = bounds.y1 - bounds.y0;
     492    }
     493
     494    // General background subtraction and measurement of stamp threshold
    434495    float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
    435496    {
    436         psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
     497        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    437498        if (ro1) {
     499            psStatsInit(bg);
    438500            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
    439501                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
    … …  
    441503                goto MATCH_ERROR;
    442504            }
    443             stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
     505            stampThresh1 = threshold * bg->robustStdev;
     506            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
    444507        }
    445508        if (ro2) {
     509            psStatsInit(bg);
    446510            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
    447511                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
    … …  
    449513                goto MATCH_ERROR;
    450514            }
    451             stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
    452         }
     515            stampThresh2 = threshold * bg->robustStdev;
     516            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
     517       }
    453518        psFree(bg);
    454519    }
     520
     521    // Just in case the iso-kernel region doesn't cover the entire image...
     522    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
     523        subMode == PM_SUBTRACTION_MODE_DUAL) {
     524        if (!conv1->image) {
     525            conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     526        }
     527        psImageInit(conv1->image, NAN);
     528        if (ro1->variance) {
     529            if (!conv1->variance) {
     530                conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     531            }
     532            psImageInit(conv1->variance, NAN);
     533        }
     534        if (subMask) {
     535            if (!conv1->mask) {
     536                conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     537            }
     538            psImageInit(conv1->mask, maskBad);
     539        }
     540    }
     541    if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
     542        subMode == PM_SUBTRACTION_MODE_DUAL) {
     543        if (!conv2->image) {
     544            conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     545        }
     546        psImageInit(conv2->image, NAN);
     547        if (ro2->variance) {
     548            if (!conv2->variance) {
     549                conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     550            }
     551            psImageInit(conv2->variance, NAN);
     552        }
     553        if (subMask) {
     554            if (!conv2->mask) {
     555                conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     556            }
     557            psImageInit(conv2->mask, maskBad);
     558        }
     559    }
     560
    455561
    456562    // Iterate over iso-kernel regions
    … …  
    459565            psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
    460566                    j * xRegions + i + 1, xRegions * yRegions);
    461             if (region) {
    462                 *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
    463                                       (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
    464                 psFree(regionString);
    465                 regionString = psRegionToString(*region);
    466                 psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
    467                         regionString, numCols, numRows);
    468             }
     567            *region = psRegionSet(bounds.x0 + (int)(i * xRegionSize),
     568                                  bounds.x0 + (int)((i + 1) * xRegionSize),
     569                                  bounds.y0 + (int)(j * yRegionSize),
     570                                  bounds.y0 + (int)((j + 1) * yRegionSize));
     571            psFree(regionString);
     572            regionString = psRegionToString(*region);
     573            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
     574                    regionString, numCols, numRows);
    469575
    470576            if (stampsName && strlen(stampsName) > 0) {
    471577                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
    472                                                         footprint, stampSpacing, sysError, subMode);
     578                                                        footprint, stampSpacing, normFrac,
     579                                                        sysError, skyError, subMode);
    473580            } else if (sources) {
    474581                stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
    475                                                            footprint, stampSpacing, sysError, subMode);
     582                                                           footprint, stampSpacing, normFrac,
     583                                                           sysError, skyError, subMode);
    476584            }
    477585
    478586            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    479587            // doesn't matter.
    480             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
    481                                       stampSpacing, sysError, size, footprint, subMode)) {
     588            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
     589                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
    482590                goto MATCH_ERROR;
    483591            }
    484592
     593
     594            // generate the window function from the set of stamps
     595            if (!pmSubtractionStampsGetWindow(stamps, size)) {
     596                psError(PS_ERR_UNKNOWN, false, "Unable to get stamp window.");
     597                goto MATCH_ERROR;
     598            }
    485599
    486600            // Define kernel basis functions
    487601            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    488                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    489                                                           stamps, footprint, optThreshold, penalty, subMode);
     602                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
     603                                                          optFWHMs, optOrder, stamps, footprint,
     604                                                          optThreshold, penalty, bounds, subMode);
    490605                if (!kernels) {
    491606                    psErrorClear();
    … …  
    496611                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    497612                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    498                                                        inner, binning, ringsOrder, penalty, subMode);
     613                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
     614                // pmSubtractionVisualShowKernels(kernels);
    499615            }
    500616
    … …  
    545661
    546662                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    547                                           stampThresh1, stampThresh2, stampSpacing, sysError,
    548                                           size, footprint, subMode)) {
    549                     goto MATCH_ERROR;
    550                 }
    551 
    552                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
    553                 if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     663                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
     664                                          sysError, skyError, size, footprint, subMode)) {
     665                    goto MATCH_ERROR;
     666                }
     667
     668                // generate the window function from the set of stamps
     669                if (!pmSubtractionStampsGetWindow(stamps, size)) {
     670                    psError(PS_ERR_UNKNOWN, false, "Unable to get stamps window.");
     671                    goto MATCH_ERROR;
     672                }
     673
     674                // XXX step 1: calculate normalization
     675                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
     676                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    554677                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    555678                    goto MATCH_ERROR;
    556679                }
    557680
     681                psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n");
     682                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
     683                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     684                    goto MATCH_ERROR;
     685                }
     686                memCheck("  solve equation");
     687
     688# if (SEPARATE)
     689                // set USED -> CALCULATE
     690                pmSubtractionStampsResetStatus (stamps);
     691
     692                // XXX step 2: calculate kernel parameters
     693                psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");
     694                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     695                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     696                    goto MATCH_ERROR;
     697                }
    558698                memCheck("  calculate equation");
    559699
    560                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
    561 
    562                 if (!pmSubtractionSolveEquation(kernels, stamps)) {
     700                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
     701                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
    563702                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    564703                    goto MATCH_ERROR;
    565704                }
    566 
    567705                memCheck("  solve equation");
    568 
     706# endif
    569707                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    570708                if (!deviations) {
    … …  
    587725            }
    588726
     727            // if we hit the max number of iterations and we have rejected stamps, re-solve
    589728            if (numRejected > 0) {
    590                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
    591                 if (!pmSubtractionSolveEquation(kernels, stamps)) {
     729                // XXX step 1: calculate normalization
     730                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
     731                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    592732                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    593733                    goto MATCH_ERROR;
    594734                }
     735
     736                // solve normalization
     737                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
     738                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
     739                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     740                    goto MATCH_ERROR;
     741                }
     742
     743# if (SEPARATE)
     744                // set USED -> CALCULATE
     745                pmSubtractionStampsResetStatus (stamps);
     746
     747                // XXX step 2: calculate kernel parameters
     748                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
     749                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     750                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     751                    goto MATCH_ERROR;
     752                }
     753
     754                // solve kernel parameters
     755                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
     756                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
     757                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     758                    goto MATCH_ERROR;
     759                }
     760                memCheck("  solve equation");
     761# endif
    595762                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    596763                if (!deviations) {
    … …  
    615782            psTrace("psModules.imcombine", 2, "Convolving...\n");
    616783            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
    617                                        kernelError, region, kernels, true, useFFT)) {
     784                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
    618785                psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
    619786                goto MATCH_ERROR;
    … …  
    8971064    assert(kernels);
    8981065
    899     psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description);
    900     if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     1066    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1067    if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    9011068        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    9021069        return false;
    9031070    }
    9041071
    905     psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description);
    906     if (!pmSubtractionSolveEquation(kernels, stamps)) {
     1072    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
     1073    if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    9071074        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    9081075        return false;
    9091076    }
     1077
     1078# if (SEPARATE)
     1079    // set USED -> CALCULATE
     1080    pmSubtractionStampsResetStatus (stamps);
     1081
     1082    psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description);
     1083    if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1084        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1085        return false;
     1086    }
     1087
     1088    psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description);
     1089    if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1090        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1091        return false;
     1092    }
     1093# endif
    9101094
    9111095    psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
    … …  
    9271111    if (numRejected > 0) {
    9281112        // Allow re-fit with reduced stamps set
    929         psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    930         if (!pmSubtractionSolveEquation(kernels, stamps)) {
     1113        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1114        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
    9311115            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    9321116            return false;
    9331117        }
     1118
     1119        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
     1120        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
     1121            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1122            return false;
     1123        }
    9341124        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
     1125
     1126# if (SEPARATE)
     1127        // set USED -> CALCULATE
     1128        pmSubtractionStampsResetStatus (stamps);
     1129
     1130        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1131        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1132            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1133            return false;
     1134        }
     1135
     1136        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
     1137        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1138            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1139            return false;
     1140        }
     1141        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
     1142# endif
     1143
    9351144        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    9361145        if (!deviations) {
    … …  
    10191228}
    10201229
     1230
     1231bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
     1232                              float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax)
     1233{
     1234    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
     1235    PS_ASSERT_PTR_NON_NULL(stampSize, false);
     1236    PS_ASSERT_VECTOR_NON_NULL(widths, false);
     1237    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
     1238    PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);
     1239    PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);
     1240    PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
     1241    PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
     1242    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
     1243    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
     1244
     1245//    float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
     1246    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
     1247
     1248    if (isfinite(scaleMin) && scale < scaleMin) {
     1249        scale = scaleMin;
     1250    }
     1251    if (isfinite(scaleMax) && scale > scaleMax) {
     1252        scale = scaleMax;
     1253    }
     1254
     1255    for (int i = 0; i < widths->n; i++) {
     1256        widths->data.F32[i] *= scale;
     1257    }
     1258    *kernelSize = *kernelSize * scale + 0.5;
     1259    *stampSize = *stampSize * scale + 0.5;
     1260
     1261    psLogMsg("psModules.imcombine", PS_LOG_INFO,
     1262             "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
     1263
     1264    return true;
     1265}
Note: See TracChangeset for help on using the changeset viewer.