Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c (revision 18238) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c (revision 18239) @@ -382,7 +382,34 @@ } + +// 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 + const pmSubtractionKernels *kernels // Kernel parameters + ) +{ + if (kernels->penalty == 0.0) { + return true; + } + + float penalty = numPixels * kernels->penalty; // Penalty value + psVector *penalties = kernels->penalties; // Penalties for each kernel component + int spatialOrder = kernels->spatialOrder; // Order of spatial variations + int numKernels = kernels->num; // Number of kernel components + for (int i = 0; i < numKernels; i++) { + 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]; + } + } + } + + return true; +} + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Semi-public functions -// XXX EAM these cannot be inline :: talk to Josh +// XXX We might like to define these functions as "extern inline" but gcc currently doesn't handle this in c99 +// mode. See http://gcc.gnu.org/ml/gcc/2006-11/msg00006.html ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -675,4 +702,5 @@ psVectorInit(sumVector, 0.0); psImageInit(sumMatrix, 0.0); + int numStamps = 0; // Number of good stamps for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -680,6 +708,8 @@ (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); - } - } + numStamps++; + } + } + calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector, kernels); psVector *permutation = NULL; // Permutation vector, required for LU decomposition @@ -716,4 +746,5 @@ psVectorInit(sumVector2, 0.0); + int numStamps = 0; // Number of good stamps for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -724,6 +755,9 @@ (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1); (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); - } - } + numStamps++; + } + } + calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector1, kernels); + calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector2, kernels); #if 0 Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.c (revision 18238) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.c (revision 18239) @@ -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] = PS_SQR(u) + PS_SQR(v); psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); @@ -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] * sqrtf(PS_SQR(u) + PS_SQR(v)); } } @@ -146,4 +151,5 @@ } kernels->preCalc->data[index] = preCalc; + kernels->penalties->data.F32[index] = moment; psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, @@ -173,4 +179,6 @@ kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); kernels->preCalc = psArrayAlloc(numBasisFunctions); + kernels->penalty = 0.0; + kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); kernels->uStop = NULL; kernels->vStop = NULL; @@ -324,4 +332,7 @@ psFree(widths); + psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels."); + psVectorInit(kernels->penalties, 0.0); + return kernels; } @@ -413,4 +424,7 @@ psFree(start); psFree(stop); + + psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels."); + psVectorInit(kernels->penalties, 0.0); return kernels; @@ -520,4 +534,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 +542,5 @@ uCoords->n = vCoords->n = poly->n = 1; radiusLast = 0; + moment = 0.0; } else { int j = 0; // Index for data @@ -546,4 +562,5 @@ poly->data.F32[j] = polyVal; norm += polyVal; + moment += polyVal * sqrtf(PS_SQR(u) + PS_SQR(v)); psVectorExtend(uCoords, RINGS_BUFFER, 1); @@ -571,4 +588,5 @@ psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); } + moment /= norm; } @@ -578,4 +596,5 @@ kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; + kernels->penalties->data.F32[index] = moment; psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h (revision 18238) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h (revision 18239) @@ -33,4 +33,6 @@ 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 + psVector *penalties; ///< Penalty for each kernel component int size; ///< The half-size of the kernel int inner; ///< The size of an inner region Index: /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c (revision 18238) +++ /branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c (revision 18239) @@ -430,5 +430,9 @@ // Generate image with convolution kernels int fullSize = 2 * size + 1 + 1; // Full size of kernel - psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32); + int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize; + psImage *convKernels = psImageAlloc((mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) * + imageSize - 1 + + (mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0), + imageSize - 1, PS_TYPE_F32); psImageInit(convKernels, NAN); for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) { @@ -451,4 +455,27 @@ } psFree(kernel); + + if (mode == PM_SUBTRACTION_MODE_DUAL) { + kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, + (float)j / (float)KERNEL_MOSAIC, + true); // Image of the kernel + if (!kernel) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); + psFree(convKernels); + goto MATCH_ERROR; + } + + if (psImageOverlaySection(convKernels, kernel, + (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize + + 4, + (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image."); + psFree(kernel); + psFree(convKernels); + goto MATCH_ERROR; + } + psFree(kernel); + } + } }