Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 18146) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 18287) @@ -24,4 +24,5 @@ psFree(kernels->vStop); psFree(kernels->preCalc); + psFree(kernels->penalties); psFree(kernels->solution1); psFree(kernels->solution2); @@ -60,4 +61,5 @@ kernels->v = psVectorRealloc(kernels->v, start + numNew); kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew); + kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew); kernels->inner = start; @@ -74,4 +76,5 @@ kernels->v->data.S32[index] = v; kernels->preCalc->data[index] = NULL; + 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); @@ -84,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); @@ -104,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", @@ -122,8 +125,10 @@ psKernel *preCalc = psKernelAlloc(-size, size, -size, size); double sum = 0.0; // Normalisation + double moment = 0.0; // Moment, for penalty for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { 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] * (PS_SQR(u) + PS_SQR(v)); } } @@ -146,7 +151,8 @@ } kernels->preCalc->data[index] = preCalc; - - 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)); } } @@ -161,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 @@ -173,4 +180,6 @@ kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); kernels->preCalc = psArrayAlloc(numBasisFunctions); + kernels->penalty = penalty; + kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); kernels->uStop = NULL; kernels->vStop = NULL; @@ -188,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); @@ -195,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); @@ -211,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; @@ -242,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); @@ -262,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", @@ -324,9 +335,13 @@ psFree(widths); + psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels."); + psVectorInit(kernels->penalties, 0.0); + return kernels; } -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); @@ -354,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", @@ -414,4 +429,7 @@ psFree(stop); + psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels."); + psVectorInit(kernels->penalties, 0.0); + return kernels; } @@ -419,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); @@ -431,6 +450,8 @@ PS_ASSERT_INT_LESS_THAN(inner, size, NULL); - pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders, mode); // Kernels + // XXX GUNK doesn't seem to work --- doesn't add the POIS components, or at least, they're not noticed + + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, + penalty, mode); // Kernels psStringPrepend(&kernels->description, "GUNK="); psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); @@ -447,5 +468,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); @@ -477,7 +498,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", @@ -520,4 +542,5 @@ psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial + double moment = 0.0; // Moment, for penalty if (i == 0) { @@ -527,4 +550,5 @@ uCoords->n = vCoords->n = poly->n = 1; radiusLast = 0; + moment = 0.0; } else { int j = 0; // Index for data @@ -546,4 +570,5 @@ poly->data.F32[j] = polyVal; norm += polyVal; + moment += polyVal * (PS_SQR(u) + PS_SQR(v)); psVectorExtend(uCoords, RINGS_BUFFER, 1); @@ -571,4 +596,5 @@ psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); } +// moment /= norm; } @@ -578,4 +604,5 @@ kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; + kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, @@ -590,19 +617,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); @@ -657,4 +685,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) { @@ -684,5 +713,6 @@ ptr++; // Eat ',' - spatialOrder = parseStringInt(ptr); + PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); + penalty = parseStringFloat(ptr); } } else if (strncmp(description, "RINGS", 5) == 0) { @@ -692,5 +722,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."); @@ -699,5 +730,5 @@ return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, - inner, binning, ringsOrder, mode); + inner, binning, ringsOrder, penalty, mode); }