Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtraction.c (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtraction.c (revision 18248) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.94 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-06-17 22:16:38 $ + * @version $Revision: 1.94.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-06-21 01:27:33 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -512,5 +512,9 @@ PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); - double totalSquareDev = 0.0; // Total square deviation from zero + // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but + // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like + // with large N. Therefore, let's just treat this as a Gaussian distribution. + + double mean = 0.0; // Mean deviation int numStamps = 0; // Number of used stamps for (int i = 0; i < stamps->num; i++) { @@ -519,9 +523,18 @@ continue; } - totalSquareDev += PS_SQR(deviations->data.F32[i]); + mean += deviations->data.F32[i]; numStamps++; } - - float rms = sqrt(totalSquareDev / (double)numStamps); // Convert to RMS + mean /= numStamps; + + double rms = 0.0; // Standard deviation + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + rms += PS_SQR(deviations->data.F32[i] - mean); + } + rms = sqrt(rms / (numStamps - 1)); if (rmsPtr) { @@ -549,9 +562,15 @@ int numRejected = 0; // Number of stamps rejected int numGood = 0; // Number of good stamps - double newSquareDev = 0.0; // New square deviation + double newMean = 0.0; // New mean for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (stamp->status == PM_SUBTRACTION_STAMP_USED) { - if (deviations->data.F32[i] > limit) { + // Should we reject stars with low deviation? Well, if this is really a Gaussian-like + // distribution and they're low, then we have the right to ask why. Isn't it suspicious that + // they're anomalously low, compared to the rest of the population which (we hope) is indicative + // of normality? Besides, the standard deviation is going to be blown up by stars that didn't + // subtract well, in which case very few (if any) stars will be legitimately rejected for being + // low. + if (deviations->data.F32[i] - mean > limit) { // Mask out the stamp in the image so you it's not found again psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, @@ -587,18 +606,18 @@ } else { numGood++; - newSquareDev += PS_SQR(deviations->data.F32[i]); + newMean += deviations->data.F32[i]; } } } + newMean /= numGood; if (numRejected > 0) { psLogMsg("psModules.imcombine", PS_LOG_INFO, - "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n", - numGood, numRejected, rms, - sqrt(newSquareDev / (double)numGood)); + "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n", + numGood, numRejected, mean, newMean); } else { psLogMsg("psModules.imcombine", PS_LOG_INFO, - "%d good stamps; 0 rejected.\nRMS deviation: %f\n", - numGood, rms); + "%d good stamps; 0 rejected.\nMean deviation: %lf\n", + numGood, mean); } Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c (revision 18248) @@ -30,5 +30,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sum += image1->kernel[y][x] * image2->kernel[y][x] / weight->kernel[y][x]; + sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x]; } } @@ -194,5 +194,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / weight->kernel[y][x]; + sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; } } @@ -217,5 +217,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double invNoise2 = 1.0 / weight->kernel[y][x]; + double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x]; double value = input->kernel[y][x] * invNoise2; sumI += value; @@ -276,5 +276,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumTC += target->kernel[y][x] * conv->kernel[y][x] / weight->kernel[y][x]; + sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; } } @@ -296,5 +296,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float value = target->kernel[y][x] / weight->kernel[y][x]; + float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x]; sumIT += value * input->kernel[y][x]; sumT += value; @@ -365,5 +365,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / weight->kernel[y][x]; + sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; } } @@ -384,6 +384,5 @@ // Add in penalty term to least-squares vector -static bool calculatePenalty(int numPixels, // Number of pixels; for normalisation - psVector *vector, // Vector to which to add in penalty term +static bool calculatePenalty(psVector *vector, // Vector to which to add in penalty term const pmSubtractionKernels *kernels // Kernel parameters ) @@ -393,5 +392,4 @@ } - float penalty = numPixels * kernels->penalty; // Penalty value psVector *penalties = kernels->penalties; // Penalties for each kernel component int spatialOrder = kernels->spatialOrder; // Order of spatial variations @@ -400,5 +398,5 @@ for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { - vector->data.F64[index] -= penalty * penalties->data.F32[i]; + vector->data.F64[index] -= penalties->data.F32[i]; } } @@ -711,5 +709,5 @@ } } - calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector, kernels); + calculatePenalty(sumVector, kernels); psVector *permutation = NULL; // Permutation vector, required for LU decomposition @@ -758,6 +756,6 @@ } } - calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector1, kernels); - calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector2, kernels); + calculatePenalty(sumVector1, kernels); + calculatePenalty(sumVector2, kernels); #if 0 @@ -812,5 +810,5 @@ psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A); assert(F->numRows == numParams && F->numCols == numParams); - float det = NAN; + float det = 0.0; psImage *Fi = psMatrixInvert(NULL, F, &det); assert(Fi->numRows == numParams && Fi->numCols == numParams); Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.c (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.c (revision 18248) @@ -76,5 +76,5 @@ kernels->v->data.S32[index] = v; kernels->preCalc->data[index] = NULL; - kernels->penalties->data.F32[index] = PS_SQR(u) + PS_SQR(v); + kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v)); psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); @@ -87,5 +87,5 @@ pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, const psVector *fwhms, const psVector *orders, - pmSubtractionMode mode) + float penalty, pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); @@ -107,7 +107,7 @@ } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, - size, spatialOrder, mode); // The kernels - psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, + spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements", @@ -130,5 +130,5 @@ sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma)); - moment += preCalc->kernel[v][u] * sqrtf(PS_SQR(u) + PS_SQR(v)); + moment += preCalc->kernel[v][u] * (PS_SQR(u) + PS_SQR(v)); } } @@ -151,8 +151,8 @@ } kernels->preCalc->data[index] = preCalc; - kernels->penalties->data.F32[index] = fabsf(moment); - - psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, - fwhms->data.F32[i], uOrder, vOrder); + kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); + + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, + fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); } } @@ -167,5 +167,6 @@ pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, - int size, int spatialOrder, pmSubtractionMode mode) + int size, int spatialOrder, float penalty, + pmSubtractionMode mode) { pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return @@ -196,5 +197,6 @@ } -pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode) +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, + pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -203,7 +205,7 @@ int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, - size, spatialOrder, mode); // The kernels - psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, size, + spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty); psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", size, spatialOrder, num); @@ -219,8 +221,8 @@ pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, const psVector *fwhms, const psVector *orders, - pmSubtractionMode mode) -{ - pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders, mode); // Kernels + float penalty, pmSubtractionMode mode) +{ + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, + penalty, mode); // Kernels if (!kernels) { return NULL; @@ -250,5 +252,5 @@ pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, - pmSubtractionMode mode) + float penalty, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -270,7 +272,8 @@ psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, - size, spatialOrder, mode); // The kernels - psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, + spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder, + penalty); psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements", @@ -339,5 +342,6 @@ -pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode) +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty, + pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -365,7 +369,7 @@ psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, - size, spatialOrder, mode); // The kernels - psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, + spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty); psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements", @@ -433,5 +437,6 @@ // Grid United with Normal Kernel pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, - const psVector *orders, int inner, pmSubtractionMode mode) + const psVector *orders, int inner, float penalty, + pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -445,6 +450,6 @@ PS_ASSERT_INT_LESS_THAN(inner, size, NULL); - pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders, mode); // Kernels + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, + penalty, mode); // Kernels psStringPrepend(&kernels->description, "GUNK="); psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); @@ -461,5 +466,5 @@ // RINGS --- just what it says pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, - pmSubtractionMode mode) + float penalty, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -491,7 +496,8 @@ int num = numRings * numPoly; // Total number of basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, - size, spatialOrder, mode); // The kernels - psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, + spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder, + penalty); psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements", @@ -562,5 +568,5 @@ poly->data.F32[j] = polyVal; norm += polyVal; - moment += polyVal * sqrtf(PS_SQR(u) + PS_SQR(v)); + moment += polyVal * (PS_SQR(u) + PS_SQR(v)); psVectorExtend(uCoords, RINGS_BUFFER, 1); @@ -596,5 +602,5 @@ kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; - kernels->penalties->data.F32[index] = fabsf(moment); + kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, @@ -609,19 +615,20 @@ pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *fwhms, const psVector *orders, int inner, - int binning, int ringsOrder, pmSubtractionMode mode) + int binning, int ringsOrder, float penalty, + pmSubtractionMode mode) { switch (type) { case PM_SUBTRACTION_KERNEL_POIS: - return pmSubtractionKernelsPOIS(size, spatialOrder, mode); + return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode); case PM_SUBTRACTION_KERNEL_ISIS: - return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode); + return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode); case PM_SUBTRACTION_KERNEL_SPAM: - return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode); + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); case PM_SUBTRACTION_KERNEL_FRIES: - return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode); + return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode); case PM_SUBTRACTION_KERNEL_GUNK: - return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode); + return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode); case PM_SUBTRACTION_KERNEL_RINGS: - return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode); + return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); @@ -676,4 +683,5 @@ int binning = 0; // Binning to use int ringsOrder = 0; // Polynomial order for rings + float penalty = 0.0; // Penalty for wideness if (strncmp(description, "ISIS", 4) == 0) { @@ -703,5 +711,6 @@ ptr++; // Eat ',' - spatialOrder = parseStringInt(ptr); + PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); + penalty = parseStringFloat(ptr); } } else if (strncmp(description, "RINGS", 5) == 0) { @@ -711,5 +720,6 @@ PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt); PARSE_STRING_NUMBER(ringsOrder, ptr, ',', parseStringInt); - PARSE_STRING_NUMBER(spatialOrder, ptr, ')', parseStringInt); + PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); + PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); } else { psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported."); @@ -718,5 +728,5 @@ return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, - inner, binning, ringsOrder, mode); + inner, binning, ringsOrder, penalty, mode); } Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h (revision 18248) @@ -33,5 +33,5 @@ psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) - float penalty; ///< Penalty value + float penalty; ///< Penalty for wideness psVector *penalties; ///< Penalty for each kernel component int size; ///< The half-size of the kernel @@ -109,4 +109,5 @@ int size, ///< Half-size of kernel int spatialOrder, ///< Order of spatial variations + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -115,4 +116,5 @@ pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims) int spatialOrder, ///< Order of spatial variations + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -123,4 +125,5 @@ const psVector *fwhms, ///< Gaussian FWHMs const psVector *orders, ///< Polynomial order of gaussians + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -131,4 +134,5 @@ const psVector *fwhms, ///< Gaussian FWHMs const psVector *orders, ///< Polynomial order of gaussians + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -139,4 +143,5 @@ int inner, ///< Inner radius to preserve unbinned int binning, ///< Kernel binning factor + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -146,4 +151,5 @@ int spatialOrder, ///< Order of spatial variations int inner, ///< Inner radius to preserve unbinned + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -155,4 +161,5 @@ const psVector *orders, ///< Polynomial order of gaussians int inner, ///< Inner radius containing grid of delta functions + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -163,4 +170,5 @@ int inner, ///< Inner radius to preserve unbinned int ringsOrder, ///< Polynomial order + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); @@ -176,4 +184,5 @@ int binning, ///< Kernel binning factor int ringsOrder, ///< Polynomial order for RINGS + float penalty, ///< Penalty for wideness pmSubtractionMode mode ///< Mode for subtraction ); Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c (revision 18248) @@ -97,7 +97,8 @@ pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *isisWidths, const psVector *isisOrders, - int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs, - int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad, - psMaskType maskBlank, float badFrac, pmSubtractionMode mode) + int inner, int ringsOrder, int binning, float penalty, + bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, + int iter, float rej, psMaskType maskBad, psMaskType maskBlank, + float badFrac, pmSubtractionMode mode) { if (mode != PM_SUBTRACTION_MODE_2) { @@ -300,5 +301,5 @@ if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold, mode); + stamps, footprint, optThreshold, penalty, mode); if (!kernels) { psErrorClear(); @@ -309,5 +310,5 @@ // Not an ISIS/GUNK kernel, or the optimum kernel search failed kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder, mode); + inner, binning, ringsOrder, penalty, mode); } Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.h (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.h (revision 18248) @@ -36,4 +36,5 @@ int ringsOrder, ///< RINGS polynomial order int binning, ///< SPAM kernel binning + float penalty, ///< Penalty for wideness bool optimum, ///< Search for optimum ISIS kernel? const psVector *optFWHMs, ///< FWHMs for optimum search Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.c (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.c (revision 18248) @@ -204,5 +204,5 @@ int spatialOrder, const psVector *fwhms, int maxOrder, const pmSubtractionStampList *stamps, int footprint, - float tolerance, pmSubtractionMode mode) + float tolerance, float penalty, pmSubtractionMode mode) { if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { @@ -231,6 +231,6 @@ psVector *orders = psVectorAlloc(numGaussians, PS_TYPE_S32); // Polynomial orders psVectorInit(orders, maxOrder); - pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders, mode); // Kernels + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, + penalty, mode); // Kernels psFree(orders); psFree(kernels->description); Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.h =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.h (revision 18247) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionParams.h (revision 18248) @@ -16,4 +16,5 @@ int footprint, ///< Convolution footprint for stamps float tolerance, ///< Maximum difference in chi^2 + float penalty, ///< Penalty for wideness pmSubtractionMode mode // Mode for subtraction );