Index: /branches/pap/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 25832) +++ /branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 25833) @@ -135,10 +135,12 @@ case PM_SUBTRACTION_KERNEL_ISIS: { psArray *preCalc = kernels->preCalc->data[i]; // Precalculated values +#if 0 psVector *xKernel = preCalc->data[0]; // Kernel in x psVector *yKernel = preCalc->data[1]; // Kernel in y // Iterating over the kernel for (int y = 0, v = -size; v <= size; y++, v++) { + float yValue = value * yKernel->data.F32[y]; for (int x = 0, u = -size; u <= size; x++, u++) { - kernel->kernel[v][u] += value * xKernel->data.F32[x] * yKernel->data.F32[y]; + kernel->kernel[v][u] += yValue * xKernel->data.F32[x]; } } @@ -147,4 +149,12 @@ kernel->kernel[0][0] -= value; } +#else + psKernel *k = preCalc->data[2]; // Kernel image + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + kernel->kernel[v][u] += value * k->kernel[v][u]; + } + } +#endif break; } @@ -597,4 +607,5 @@ case PM_SUBTRACTION_KERNEL_ISIS: { psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values +#if 1 psVector *xKernel = preCalc->data[0]; // Kernel in x psVector *yKernel = preCalc->data[1]; // Kernel in y @@ -640,4 +651,8 @@ } return convolved; +#else + return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]); +#endif + } case PM_SUBTRACTION_KERNEL_RINGS: { Index: /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 25832) +++ /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 25833) @@ -15,5 +15,5 @@ #include "pmSubtractionVisual.h" -// #define TESTING // TESTING output for debugging; may not work with threads! +//#define TESTING // TESTING output for debugging; may not work with threads! #define USE_VARIANCE // Include variance in equation? @@ -22,4 +22,282 @@ // Private (file-static) functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#if 1 +// Calculate the least-squares matrix and vector +static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated + const psKernel *input, // Input image (target) + const psKernel *reference, // Reference image (convolution source) + const psKernel *variance, // Variance image + const psArray *convolutions, // Convolutions for each kernel + const pmSubtractionKernels *kernels, // Kernels + const psImage *polyValues, // Spatial polynomial values + int footprint // (Half-)Size of stamp + ) +{ + // (I - R * sum_i a_i k_i - g) (R * k_j) = 0 + // I C_j = sum_i C_i C_j + + // Background: C_i = 1.0 + // Normalisation: C_i = R + + int numKernels = kernels->num; // Number of kernels + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + double poly[numPoly]; // Polynomial terms + double poly2[numPoly][numPoly]; // Polynomial-polynomial values + + // Evaluate polynomial-polynomial terms + for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { + double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial + poly[iIndex] = iPoly; + for (int jyOrder = 0, jIndex = 0; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex++) { + double jPoly = polyValues->data.F64[jyOrder][jxOrder]; + poly2[iIndex][jIndex] = iPoly * jPoly; + } + } + } + } + + + for (int i = 0; i < numKernels; i++) { + psKernel *iConv = convolutions->data[i]; // Convolution for index i + for (int j = i; j < numKernels; j++) { + psKernel *jConv = convolutions->data[j]; // Convolution for index j + + double sumCC = 0.0; // Sum of convolution products + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; +#ifdef USE_VARIANCE + cc /= variance->kernel[y][x]; +#endif + sumCC += cc; + } + } + + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double value = sumCC * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex] = value; + matrix->data.F64[jIndex][iIndex] = value; + } + } + } + + double sumRC = 0.0; // Sum of the reference-convolution products + double sumIC = 0.0; // Sum of the input-convolution products + double sumC = 0.0; // Sum of the convolution + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float conv = iConv->kernel[y][x]; + float in = input->kernel[y][x]; + float ref = reference->kernel[y][x]; + double ic = in * conv; + double rc = ref * conv; + double c = conv; +#ifdef USE_VARIANCE + float varVal = 1.0 / variance->kernel[y][x]; + ic *= varVal; + rc *= varVal; + c *= varVal; +#endif + sumIC += ic; + sumRC += rc; + sumC += c; + } + } + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + double normTerm = sumRC * poly[iTerm]; + double bgTerm = sumC * poly[iTerm]; + matrix->data.F64[iIndex][normIndex] = normTerm; + matrix->data.F64[normIndex][iIndex] = normTerm; + matrix->data.F64[iIndex][bgIndex] = bgTerm; + matrix->data.F64[bgIndex][iIndex] = bgTerm; + vector->data.F64[iIndex] = sumIC * poly[iTerm]; + } + } + + double sumRR = 0.0; // Sum of the reference product + double sumIR = 0.0; // Sum of the input-reference product + double sum1 = 0.0; // Sum of the background + double sumR = 0.0; // Sum of the reference + double sumI = 0.0; // Sum of the input + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double in = input->kernel[y][x]; + double ref = reference->kernel[y][x]; + double ir = in * ref; + double rr = PS_SQR(ref); +#ifdef USE_VARIANCE + float varVal = 1.0 / variance->kernel[y][x]; + rr *= varVal; + ir *= varVal; + in *= varVal; + ref *= varVal; +#endif + sumRR += rr; + sumIR += ir; + sumR += ref; + sumI += in; + sum1 += 1.0; + } + } + matrix->data.F64[normIndex][normIndex] = sumRR; + matrix->data.F64[bgIndex][bgIndex] = sum1; + matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR; + vector->data.F64[normIndex] = sumIR; + vector->data.F64[bgIndex] = sumI; + + return true; +} + +#if 0 +// Calculate the least-squares matrix and vector for dual convolution +static bool calculateDualMatrixVector(psImage *matrix1, // Least-squares matrix, updated + psVector *vector1, // Least-squares vector, updated + psImage *matrix2, // Least-squares matrix, updated + psVector *vector2, // Least-squares vector, updated + psImage *matrixX, // Cross-matrix + const psKernel *input, // Input image (target) + const psKernel *reference, // Reference image (convolution source) + const psKernel *variance, // Variance image + const psArray *convolutions, // Convolutions for each kernel + const pmSubtractionKernels *kernels, // Kernels + const psImage *polyValues, // Spatial polynomial values + int footprint // (Half-)Size of stamp + ) +{ + // (I - R * sum_i a_i k_i - g) (R * k_j) = 0 + // I C_j = sum_i C_i C_j + + // Background: C_i = 1.0 + // Normalisation: C_i = R + + int numKernels = kernels->num; // Number of kernels + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + double poly[numPoly]; // Polynomial terms + double poly2[numPoly][numPoly]; // Polynomial-polynomial values + + // Evaluate polynomial-polynomial terms + for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { + double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial + poly[iIndex] = iPoly; + for (int jyOrder = 0, jIndex = 0; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex++) { + double jPoly = polyValues->data.F64[jyOrder][jxOrder]; + poly2[iIndex][jIndex] = iPoly * jPoly; + } + } + } + } + + + for (int i = 0; i < numKernels; i++) { + psKernel *iConv = convolutions->data[i]; // Convolution for index i + for (int j = i; j < numKernels; j++) { + psKernel *jConv = convolutions->data[j]; // Convolution for index j + + double sumCC = 0.0; // Sum of convolution products + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; +#ifdef USE_VARIANCE + cc /= variance->kernel[y][x]; +#endif + sumCC += cc; + } + } + + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double value = sumCC * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex] = value; + matrix->data.F64[jIndex][iIndex] = value; + } + } + } + + double sumRC = 0.0; // Sum of the reference-convolution products + double sumIC = 0.0; // Sum of the input-convolution products + double sumC = 0.0; // Sum of the convolution + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float conv = iConv->kernel[y][x]; + float in = input->kernel[y][x]; + float ref = reference->kernel[y][x]; + double ic = in * conv; + double rc = ref * conv; + double c = conv; +#ifdef USE_VARIANCE + float varVal = 1.0 / variance->kernel[y][x]; + ic *= varVal; + rc *= varVal; + c *= varVal; +#endif + sumIC += ic; + sumRC += rc; + sumC += c; + } + } + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + double normTerm = sumRC * poly[iTerm]; + double bgTerm = sumC * poly[iTerm]; + matrix->data.F64[iIndex][normIndex] = normTerm; + matrix->data.F64[normIndex][iIndex] = normTerm; + matrix->data.F64[iIndex][bgIndex] = bgTerm; + matrix->data.F64[bgIndex][iIndex] = bgTerm; + vector->data.F64[iIndex] = sumIC * poly[iTerm]; + } + } + + double sumRR = 0.0; // Sum of the reference product + double sumIR = 0.0; // Sum of the input-reference product + double sum1 = 0.0; // Sum of the background + double sumR = 0.0; // Sum of the reference + double sumI = 0.0; // Sum of the input + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double in = input->kernel[y][x]; + double ref = reference->kernel[y][x]; + double ir = in * ref; + double rr = PS_SQR(ref); +#ifdef USE_VARIANCE + float varVal = 1.0 / variance->kernel[y][x]; + rr *= varVal; + ir *= varVal; + in *= varVal; + ref *= varVal; +#endif + sumRR += rr; + sumIR += ir; + sumR += ref; + sumI += in; + sum1 += 1.0; + } + } + matrix->data.F64[normIndex][normIndex] = sumRR; + matrix->data.F64[bgIndex][bgIndex] = sum1; + matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR; + vector->data.F64[normIndex] = sumIR; + vector->data.F64[bgIndex] = sumI; + + return true; +} +#endif +#endif + // Calculate the sum over a stamp product @@ -582,14 +860,26 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: +#if 0 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, stamp->image2, stamp->variance, polyValues, footprint, true); +#else + status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, + stamp->variance, stamp->convolutions1, kernels, polyValues, + footprint); +#endif break; case PM_SUBTRACTION_MODE_2: +#if 0 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, stamp->image1, stamp->variance, polyValues, footprint, true); +#else + status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2, + stamp->variance, stamp->convolutions2, kernels, polyValues, + footprint); +#endif break; case PM_SUBTRACTION_MODE_DUAL: @@ -629,5 +919,5 @@ #ifdef TESTING - if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { + { psString matrixName = NULL; psStringAppend(&matrixName, "matrix1_%d.fits", index); Index: /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c (revision 25832) +++ /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c (revision 25833) @@ -140,7 +140,8 @@ for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { - psArray *preCalc = psArrayAlloc(2); // Array to hold precalculated values + psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel + psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size); // Kernel // Calculate moments @@ -149,4 +150,5 @@ for (int u = -size, x = 0; u <= size; u++, x++) { double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel + kernel->kernel[v][u] = value; moment += value * (PS_SQR(u) + PS_SQR(v)); } @@ -164,4 +166,6 @@ psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); + psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); + kernel->kernel[0][0] -= 1.0; moment *= PS_SQR(sum); }