IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18248


Ignore:
Timestamp:
Jun 20, 2008, 3:27:33 PM (18 years ago)
Author:
Paul Price
Message:

Dual convolution seems to be working. Removed division by 'weight' when generating equation --- seems to bias the fits.

Location:
branches/pap_branch_080617/psModules/src/imcombine
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtraction.c

    r18163 r18248  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.94 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2008-06-17 22:16:38 $
     6 *  @version $Revision: 1.94.2.1 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2008-06-21 01:27:33 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
     
    512512    PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1);
    513513
    514     double totalSquareDev = 0.0;        // Total square deviation from zero
     514    // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but
     515    // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like
     516    // with large N.  Therefore, let's just treat this as a Gaussian distribution.
     517
     518    double mean = 0.0;                  // Mean deviation
    515519    int numStamps = 0;                  // Number of used stamps
    516520    for (int i = 0; i < stamps->num; i++) {
     
    519523            continue;
    520524        }
    521         totalSquareDev += PS_SQR(deviations->data.F32[i]);
     525        mean += deviations->data.F32[i];
    522526        numStamps++;
    523527    }
    524 
    525     float rms = sqrt(totalSquareDev / (double)numStamps); // Convert to RMS
     528    mean /= numStamps;
     529
     530    double rms = 0.0;                   // Standard deviation
     531    for (int i = 0; i < stamps->num; i++) {
     532        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     533        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
     534            continue;
     535        }
     536        rms += PS_SQR(deviations->data.F32[i] - mean);
     537    }
     538    rms = sqrt(rms / (numStamps - 1));
    526539
    527540    if (rmsPtr) {
     
    549562    int numRejected = 0;                // Number of stamps rejected
    550563    int numGood = 0;                    // Number of good stamps
    551     double newSquareDev = 0.0;          // New square deviation
     564    double newMean = 0.0;               // New mean
    552565    for (int i = 0; i < stamps->num; i++) {
    553566        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    554567        if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
    555             if (deviations->data.F32[i] > limit) {
     568            // Should we reject stars with low deviation?  Well, if this is really a Gaussian-like
     569            // distribution and they're low, then we have the right to ask why.  Isn't it suspicious that
     570            // they're anomalously low, compared to the rest of the population which (we hope) is indicative
     571            // of normality?  Besides, the standard deviation is going to be blown up by stars that didn't
     572            // subtract well, in which case very few (if any) stars will be legitimately rejected for being
     573            // low.
     574            if (deviations->data.F32[i] - mean > limit) {
    556575                // Mask out the stamp in the image so you it's not found again
    557576                psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
     
    587606            } else {
    588607                numGood++;
    589                 newSquareDev += PS_SQR(deviations->data.F32[i]);
     608                newMean += deviations->data.F32[i];
    590609            }
    591610        }
    592611    }
     612    newMean /= numGood;
    593613
    594614    if (numRejected > 0) {
    595615        psLogMsg("psModules.imcombine", PS_LOG_INFO,
    596                  "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n",
    597                  numGood, numRejected, rms,
    598                  sqrt(newSquareDev / (double)numGood));
     616                 "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n",
     617                 numGood, numRejected, mean, newMean);
    599618    } else {
    600619        psLogMsg("psModules.imcombine", PS_LOG_INFO,
    601                  "%d good stamps; 0 rejected.\nRMS deviation: %f\n",
    602                  numGood, rms);
     620                 "%d good stamps; 0 rejected.\nMean deviation: %lf\n",
     621                 numGood, mean);
    603622    }
    604623
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c

    r18239 r18248  
    3030    for (int y = - footprint; y <= footprint; y++) {
    3131        for (int x = - footprint; x <= footprint; x++) {
    32             sum += image1->kernel[y][x] * image2->kernel[y][x] / weight->kernel[y][x];
     32            sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x];
    3333        }
    3434    }
     
    194194            for (int y = - footprint; y <= footprint; y++) {
    195195                for (int x = - footprint; x <= footprint; x++) {
    196                     sumC += conv->kernel[y][x] / weight->kernel[y][x];
     196                    sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
    197197                }
    198198            }
     
    217217        for (int y = - footprint; y <= footprint; y++) {
    218218            for (int x = - footprint; x <= footprint; x++) {
    219                 double invNoise2 = 1.0 / weight->kernel[y][x];
     219                double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x];
    220220                double value = input->kernel[y][x] * invNoise2;
    221221                sumI += value;
     
    276276        for (int y = - footprint; y <= footprint; y++) {
    277277            for (int x = - footprint; x <= footprint; x++) {
    278                     sumTC += target->kernel[y][x] * conv->kernel[y][x] / weight->kernel[y][x];
     278                sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
    279279            }
    280280        }
     
    296296        for (int y = - footprint; y <= footprint; y++) {
    297297            for (int x = - footprint; x <= footprint; x++) {
    298                 float value = target->kernel[y][x] / weight->kernel[y][x];
     298                float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x];
    299299                sumIT += value * input->kernel[y][x];
    300300                sumT += value;
     
    365365        for (int y = - footprint; y <= footprint; y++) {
    366366            for (int x = - footprint; x <= footprint; x++) {
    367                 sumC += conv->kernel[y][x] / weight->kernel[y][x];
     367                sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
    368368            }
    369369        }
     
    384384
    385385// Add in penalty term to least-squares vector
    386 static bool calculatePenalty(int numPixels, // Number of pixels; for normalisation
    387                              psVector *vector, // Vector to which to add in penalty term
     386static bool calculatePenalty(psVector *vector, // Vector to which to add in penalty term
    388387                             const pmSubtractionKernels *kernels // Kernel parameters
    389388    )
     
    393392    }
    394393
    395     float penalty = numPixels * kernels->penalty; // Penalty value
    396394    psVector *penalties = kernels->penalties; // Penalties for each kernel component
    397395    int spatialOrder = kernels->spatialOrder; // Order of spatial variations
     
    400398        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
    401399            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    402                 vector->data.F64[index] -= penalty * penalties->data.F32[i];
     400                vector->data.F64[index] -= penalties->data.F32[i];
    403401            }
    404402        }
     
    711709            }
    712710        }
    713         calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector, kernels);
     711        calculatePenalty(sumVector, kernels);
    714712
    715713        psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
     
    758756            }
    759757        }
    760         calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector1, kernels);
    761         calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector2, kernels);
     758        calculatePenalty(sumVector1, kernels);
     759        calculatePenalty(sumVector2, kernels);
    762760
    763761#if 0
     
    812810        psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A);
    813811        assert(F->numRows == numParams && F->numCols == numParams);
    814         float det = NAN;
     812        float det = 0.0;
    815813        psImage *Fi = psMatrixInvert(NULL, F, &det);
    816814        assert(Fi->numRows == numParams && Fi->numCols == numParams);
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.c

    r18240 r18248  
    7676            kernels->v->data.S32[index] = v;
    7777            kernels->preCalc->data[index] = NULL;
    78             kernels->penalties->data.F32[index] = PS_SQR(u) + PS_SQR(v);
     78            kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v));
    7979
    8080            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
     
    8787pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
    8888                                                    const psVector *fwhms, const psVector *orders,
    89                                                     pmSubtractionMode mode)
     89                                                    float penalty, pmSubtractionMode mode)
    9090{
    9191    PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
     
    107107    }
    108108
    109     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS,
    110                                                               size, spatialOrder, mode); // The kernels
    111     psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder);
     109    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
     110                                                              spatialOrder, penalty, mode); // The kernels
     111    psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    112112
    113113    psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements",
     
    130130                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
    131131                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
    132                         moment += preCalc->kernel[v][u] * sqrtf(PS_SQR(u) + PS_SQR(v));
     132                        moment += preCalc->kernel[v][u] * (PS_SQR(u) + PS_SQR(v));
    133133                    }
    134134                }
     
    151151                }
    152152                kernels->preCalc->data[index] = preCalc;
    153                 kernels->penalties->data.F32[index] = fabsf(moment);
    154 
    155                 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
    156                         fwhms->data.F32[i], uOrder, vOrder);
     153                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
     154
     155                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
     156                        fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
    157157            }
    158158        }
     
    167167
    168168pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
    169                                                 int size, int spatialOrder, pmSubtractionMode mode)
     169                                                int size, int spatialOrder, float penalty,
     170                                                pmSubtractionMode mode)
    170171{
    171172    pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
     
    196197}
    197198
    198 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode)
     199pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
     200                                               pmSubtractionMode mode)
    199201{
    200202    PS_ASSERT_INT_POSITIVE(size, NULL);
     
    203205    int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions
    204206
    205     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS,
    206                                                               size, spatialOrder, mode); // The kernels
    207     psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder);
     207    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, size,
     208                                                              spatialOrder, penalty, mode); // The kernels
     209    psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
    208210    psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
    209211             size, spatialOrder, num);
     
    219221pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
    220222                                               const psVector *fwhms, const psVector *orders,
    221                                                pmSubtractionMode mode)
    222 {
    223     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
    224                                                                   fwhms, orders, mode); // Kernels
     223                                               float penalty, pmSubtractionMode mode)
     224{
     225    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
     226                                                                  penalty, mode); // Kernels
    225227    if (!kernels) {
    226228        return NULL;
     
    250252
    251253pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
    252                                                pmSubtractionMode mode)
     254                                               float penalty, pmSubtractionMode mode)
    253255{
    254256    PS_ASSERT_INT_POSITIVE(size, NULL);
     
    270272    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
    271273
    272     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM,
    273                                                               size, spatialOrder, mode); // The kernels
    274     psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder);
     274    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
     275                                                              spatialOrder, penalty, mode); // The kernels
     276    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
     277                   penalty);
    275278
    276279    psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements",
     
    339342
    340343
    341 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode)
     344pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
     345                                                pmSubtractionMode mode)
    342346{
    343347    PS_ASSERT_INT_POSITIVE(size, NULL);
     
    365369    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
    366370
    367     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES,
    368                                                               size, spatialOrder, mode); // The kernels
    369     psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder);
     371    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
     372                                                              spatialOrder, penalty, mode); // The kernels
     373    psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
    370374
    371375    psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements",
     
    433437// Grid United with Normal Kernel
    434438pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
    435                                                const psVector *orders, int inner, pmSubtractionMode mode)
     439                                               const psVector *orders, int inner, float penalty,
     440                                               pmSubtractionMode mode)
    436441{
    437442    PS_ASSERT_INT_POSITIVE(size, NULL);
     
    445450    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
    446451
    447     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
    448                                                                   fwhms, orders, mode); // Kernels
     452    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
     453                                                                  penalty, mode); // Kernels
    449454    psStringPrepend(&kernels->description, "GUNK=");
    450455    psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
     
    461466// RINGS --- just what it says
    462467pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
    463                                                 pmSubtractionMode mode)
     468                                                float penalty, pmSubtractionMode mode)
    464469{
    465470    PS_ASSERT_INT_POSITIVE(size, NULL);
     
    491496    int num = numRings * numPoly; // Total number of basis functions
    492497
    493     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS,
    494                                                               size, spatialOrder, mode); // The kernels
    495     psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder);
     498    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
     499                                                              spatialOrder, penalty, mode); // The kernels
     500    psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
     501                   penalty);
    496502
    497503    psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements",
     
    562568                                    poly->data.F32[j] = polyVal;
    563569                                    norm += polyVal;
    564                                     moment += polyVal * sqrtf(PS_SQR(u) + PS_SQR(v));
     570                                    moment += polyVal * (PS_SQR(u) + PS_SQR(v));
    565571
    566572                                    psVectorExtend(uCoords, RINGS_BUFFER, 1);
     
    596602                kernels->u->data.S32[index] = uOrder;
    597603                kernels->v->data.S32[index] = vOrder;
    598                 kernels->penalties->data.F32[index] = fabsf(moment);
     604                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    599605
    600606                psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,
     
    609615pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
    610616                                                   const psVector *fwhms, const psVector *orders, int inner,
    611                                                    int binning, int ringsOrder, pmSubtractionMode mode)
     617                                                   int binning, int ringsOrder, float penalty,
     618                                                   pmSubtractionMode mode)
    612619{
    613620    switch (type) {
    614621      case PM_SUBTRACTION_KERNEL_POIS:
    615         return pmSubtractionKernelsPOIS(size, spatialOrder, mode);
     622        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode);
    616623      case PM_SUBTRACTION_KERNEL_ISIS:
    617         return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode);
     624        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
    618625      case PM_SUBTRACTION_KERNEL_SPAM:
    619         return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode);
     626        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
    620627      case PM_SUBTRACTION_KERNEL_FRIES:
    621         return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode);
     628        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode);
    622629      case PM_SUBTRACTION_KERNEL_GUNK:
    623         return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode);
     630        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode);
    624631      case PM_SUBTRACTION_KERNEL_RINGS:
    625         return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode);
     632        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode);
    626633      default:
    627634        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
     
    676683    int binning = 0;                    // Binning to use
    677684    int ringsOrder = 0;                 // Polynomial order for rings
     685    float penalty = 0.0;                // Penalty for wideness
    678686
    679687    if (strncmp(description, "ISIS", 4) == 0) {
     
    703711
    704712            ptr++;                      // Eat ','
    705             spatialOrder = parseStringInt(ptr);
     713            PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
     714            penalty = parseStringFloat(ptr);
    706715        }
    707716    } else if (strncmp(description, "RINGS", 5) == 0) {
     
    711720        PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt);
    712721        PARSE_STRING_NUMBER(ringsOrder, ptr, ',', parseStringInt);
    713         PARSE_STRING_NUMBER(spatialOrder, ptr, ')', parseStringInt);
     722        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
     723        PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
    714724    } else {
    715725        psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported.");
     
    718728
    719729    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders,
    720                                         inner, binning, ringsOrder, mode);
     730                                        inner, binning, ringsOrder, penalty, mode);
    721731}
    722732
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h

    r18239 r18248  
    3333    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    3434    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
    35     float penalty;                      ///< Penalty value
     35    float penalty;                      ///< Penalty for wideness
    3636    psVector *penalties;                ///< Penalty for each kernel component
    3737    int size;                           ///< The half-size of the kernel
     
    109109                                                int size, ///< Half-size of kernel
    110110                                                int spatialOrder, ///< Order of spatial variations
     111                                                float penalty, ///< Penalty for wideness
    111112                                                pmSubtractionMode mode ///< Mode for subtraction
    112113    );
     
    115116pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims)
    116117                                               int spatialOrder, ///< Order of spatial variations
     118                                               float penalty, ///< Penalty for wideness
    117119                                               pmSubtractionMode mode ///< Mode for subtraction
    118120    );
     
    123125                                                    const psVector *fwhms, ///< Gaussian FWHMs
    124126                                                    const psVector *orders, ///< Polynomial order of gaussians
     127                                                    float penalty, ///< Penalty for wideness
    125128                                                    pmSubtractionMode mode ///< Mode for subtraction
    126129    );
     
    131134                                               const psVector *fwhms, ///< Gaussian FWHMs
    132135                                               const psVector *orders, ///< Polynomial order of gaussians
     136                                               float penalty, ///< Penalty for wideness
    133137                                               pmSubtractionMode mode ///< Mode for subtraction
    134138                                               );
     
    139143                                               int inner, ///< Inner radius to preserve unbinned
    140144                                               int binning, ///< Kernel binning factor
     145                                               float penalty, ///< Penalty for wideness
    141146                                               pmSubtractionMode mode ///< Mode for subtraction
    142147    );
     
    146151                                                int spatialOrder, ///< Order of spatial variations
    147152                                                int inner, ///< Inner radius to preserve unbinned
     153                                                float penalty, ///< Penalty for wideness
    148154                                                pmSubtractionMode mode ///< Mode for subtraction
    149155    );
     
    155161                                               const psVector *orders, ///< Polynomial order of gaussians
    156162                                               int inner, ///< Inner radius containing grid of delta functions
     163                                               float penalty, ///< Penalty for wideness
    157164                                               pmSubtractionMode mode ///< Mode for subtraction
    158165    );
     
    163170                                                int inner, ///< Inner radius to preserve unbinned
    164171                                                int ringsOrder, ///< Polynomial order
     172                                                float penalty, ///< Penalty for wideness
    165173                                                pmSubtractionMode mode ///< Mode for subtraction
    166174    );
     
    176184                                                   int binning, ///< Kernel binning factor
    177185                                                   int ringsOrder, ///< Polynomial order for RINGS
     186                                                   float penalty, ///< Penalty for wideness
    178187                                                   pmSubtractionMode mode ///< Mode for subtraction
    179188    );
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c

    r18239 r18248  
    9797                        pmSubtractionKernelsType type, int size, int spatialOrder,
    9898                        const psVector *isisWidths, const psVector *isisOrders,
    99                         int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs,
    100                         int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
    101                         psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
     99                        int inner, int ringsOrder, int binning, float penalty,
     100                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
     101                        int iter, float rej, psMaskType maskBad, psMaskType maskBlank,
     102                        float badFrac, pmSubtractionMode mode)
    102103{
    103104    if (mode != PM_SUBTRACTION_MODE_2) {
     
    300301            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    301302                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    302                                                           stamps, footprint, optThreshold, mode);
     303                                                          stamps, footprint, optThreshold, penalty, mode);
    303304                if (!kernels) {
    304305                    psErrorClear();
     
    309310                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    310311                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    311                                                        inner, binning, ringsOrder, mode);
     312                                                       inner, binning, ringsOrder, penalty, mode);
    312313            }
    313314
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.h

    r17783 r18248  
    3636                        int ringsOrder, ///< RINGS polynomial order
    3737                        int binning,    ///< SPAM kernel binning
     38                        float penalty,  ///< Penalty for wideness
    3839                        bool optimum,   ///< Search for optimum ISIS kernel?
    3940                        const psVector *optFWHMs, ///< FWHMs for optimum search
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.c

    r15756 r18248  
    204204                                                      int spatialOrder, const psVector *fwhms, int maxOrder,
    205205                                                      const pmSubtractionStampList *stamps, int footprint,
    206                                                       float tolerance, pmSubtractionMode mode)
     206                                                      float tolerance, float penalty, pmSubtractionMode mode)
    207207{
    208208    if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
     
    231231    psVector *orders = psVectorAlloc(numGaussians, PS_TYPE_S32); // Polynomial orders
    232232    psVectorInit(orders, maxOrder);
    233     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
    234                                                                   fwhms, orders, mode); // Kernels
     233    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
     234                                                                  penalty, mode); // Kernels
    235235    psFree(orders);
    236236    psFree(kernels->description);
  • branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.h

    r15443 r18248  
    1616                                                      int footprint, ///< Convolution footprint for stamps
    1717                                                      float tolerance, ///< Maximum difference in chi^2
     18                                                      float penalty, ///< Penalty for wideness
    1819                                                      pmSubtractionMode mode // Mode for subtraction
    1920    );
Note: See TracChangeset for help on using the changeset viewer.