Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 25998) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 25999) @@ -29,4 +29,5 @@ #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time #define MIN_SAMPLE_STATS 7 // Minimum number to use sample statistics; otherwise use quartiles +#define USE_SYS_ERR // Use systematic error image? ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -36,5 +37,5 @@ // Generate the kernel to apply to the variance from the normal kernel static psKernel *varianceKernel(psKernel *out, // Output kernel - psKernel *normalKernel // Normal kernel + const psKernel *normalKernel // Normal kernel ) { @@ -48,12 +49,8 @@ // Take the square of the normal kernel - double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels + double sumVariance = 0.0; // Sum of the variance kernels for (int v = yMin; v <= yMax; v++) { for (int u = xMin; u <= xMax; u++) { - float value = normalKernel->kernel[v][u]; // Value of interest - float value2 = PS_SQR(value); // Value squared - sumNormal += value; - sumVariance += value2; - out->kernel[v][u] = value2; + sumVariance += out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]); } } @@ -65,4 +62,40 @@ return out; } + +// Contribute to an image of the solved kernel component for ISIS +static void solvedKernelISIS(psKernel *kernel, // Kernel, updated + const pmSubtractionKernels *kernels, // Kernel basis functions + float value, // Normalisation value for basis function + int index // Index of basis function of interest + ) +{ + int size = kernels->size; // Kernel half-size + psArray *preCalc = kernels->preCalc->data[index]; // 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] += yValue * xKernel->data.F32[x]; + } + } + // Photometric scaling for even kernels only + if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { + 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 + + return; +} + // Generate an image of the solved kernel @@ -70,4 +103,5 @@ const pmSubtractionKernels *kernels, // Kernel basis functions const psImage *polyValues, // Spatial polynomial values + bool normalise, // Add normalisation? bool wantDual // Want the dual (second) kernel? ) @@ -115,13 +149,5 @@ case PM_SUBTRACTION_KERNEL_GUNK: { if (i < kernels->inner) { - // Using pre-calculated function - psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values - // Iterating over the kernel - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - kernel->kernel[v][u] += preCalc->kernel[v][u] * value; - // Photometric scaling is built into the preCalc kernel --- no subtraction! - } - } + solvedKernelISIS(kernel, kernels, value, i); } else { // Using delta function @@ -134,17 +160,5 @@ } case PM_SUBTRACTION_KERNEL_ISIS: { - psArray *preCalc = kernels->preCalc->data[i]; // Precalculated values - 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++) { - for (int x = 0, u = -size; u <= size; x++, u++) { - kernel->kernel[v][u] += value * xKernel->data.F32[x] * yKernel->data.F32[y]; - } - } - // Photometric scaling for even kernels only - if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { - kernel->kernel[0][0] -= value; - } + solvedKernelISIS(kernel, kernels, value, i); break; } @@ -168,6 +182,8 @@ } - // Put in the normalisation component - kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels)); + if (normalise) { + // Put in the normalisation component + kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels)); + } return kernel; @@ -339,5 +355,5 @@ ) { - *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual); + *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, true, wantDual); if (variance || subMask) { *kernelVariance = varianceKernel(*kernelVariance, *kernelImage); @@ -413,4 +429,5 @@ } +#ifdef USE_SYS_ERR // Generate an image that can be used to track systematic errors static psImage *subtractionSysErrImage(const psImage *image, // Image from which to make sys err image @@ -433,4 +450,64 @@ return sys; +} +#endif + +// Convolve a stamp using an ISIS kernel basis function +static psKernel *convolveStampISIS(const psKernel *image, // Image to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int index, // Index of basis function of interest + int footprint // Half-size of stamp + ) +{ + psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data +#if 1 + // Convolving using separable convolution + psVector *xKernel = preCalc->data[0]; // Kernel in x + psVector *yKernel = preCalc->data[1]; // Kernel in y + int size = kernels->size; // Size of kernel + + // Convolve in x + // Need to convolve a bit more than the footprint, for the y convolution + int yMin = -size - footprint, yMax = size + footprint; // Range for y + psKernel *temp = psKernelAlloc(yMin, yMax, + -footprint, footprint); // Temporary convolution; NOTE: wrong way! + for (int y = yMin; y <= yMax; y++) { + for (int x = -footprint; x <= footprint; x++) { + float value = 0.0; // Value of convolved pixel + int uMin = x - size, uMax = x + size; // Range for u + psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values + psF32 *imageData = &image->kernel[y][uMin]; // Image values + for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) { + value += *xKernelData * *imageData; + } + temp->kernel[x][y] = value; // NOTE: putting in wrong way! + } + } + + // Convolve in y + psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint);// Convolved image + for (int x = -footprint; x <= footprint; x++) { + for (int y = -footprint; y <= footprint; y++) { + float value = 0.0; // Value of convolved pixel + int vMin = y - size, vMax = y + size; // Range for v + psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values + psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way! + for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) { + value += *yKernelData * *imageData; + } + convolved->kernel[y][x] = value; + } + } + psFree(temp); + + // Photometric scaling for even kernels only + if (kernels->u->data.S32[index] % 2 == 0 && kernels->v->data.S32[index] % 2 == 0) { + convolveSub(convolved, image, footprint); + } + return convolved; +#else + // Convolving using precalculated kernel + return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]); +#endif } @@ -456,5 +533,5 @@ psKernel *kernel; // Kernel to use if (!preKernel) { - kernel = solvedKernel(NULL, kernels, polyValues, wantDual); + kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual); } else { kernel = psMemIncrRefCounter(preKernel); @@ -586,5 +663,5 @@ if (index < kernels->inner) { // Photometric scaling is already built in to the precalculated kernel - return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]); + return convolveStampISIS(image, kernels, index, footprint); } // Using delta function @@ -596,48 +673,5 @@ } case PM_SUBTRACTION_KERNEL_ISIS: { - psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values - psVector *xKernel = preCalc->data[0]; // Kernel in x - psVector *yKernel = preCalc->data[1]; // Kernel in y - int size = kernels->size; // Size of kernel - - // Convolve in x - // Need to convolve a bit more than the footprint, for the y convolution - int yMin = -size - footprint, yMax = size + footprint; // Range for y - psKernel *temp = psKernelAlloc(yMin, yMax, - -footprint, footprint); // Temporary convolution; NOTE: wrong way! - for (int y = yMin; y <= yMax; y++) { - for (int x = -footprint; x <= footprint; x++) { - float value = 0.0; // Value of convolved pixel - int uMin = x - size, uMax = x + size; // Range for u - psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values - psF32 *imageData = &image->kernel[y][uMin]; // Image values - for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) { - value += *xKernelData * *imageData; - } - temp->kernel[x][y] = value; // NOTE: putting in wrong way! - } - } - - // Convolve in y - psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint);// Convolved image - for (int x = -footprint; x <= footprint; x++) { - for (int y = -footprint; y <= footprint; y++) { - float value = 0.0; // Value of convolved pixel - int vMin = y - size, vMax = y + size; // Range for v - psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values - psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way! - for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) { - value += *yKernelData * *imageData; - } - convolved->kernel[y][x] = value; - } - } - psFree(temp); - - // Photometric scaling for even kernels only - if (kernels->u->data.S32[index] % 2 == 0 && kernels->v->data.S32[index] % 2 == 0) { - convolveSub(convolved, image, footprint); - } - return convolved; + return convolveStampISIS(image, kernels, index, footprint); } case PM_SUBTRACTION_KERNEL_RINGS: { @@ -901,5 +935,5 @@ psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial - psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual); // The appropriate kernel psFree(polyValues); @@ -932,5 +966,5 @@ psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); - psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual); // The appropriate kernel psFree(polyValues); @@ -949,5 +983,5 @@ } -#if 0 +#if 1 psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) { @@ -957,15 +991,41 @@ PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - psArray *images = psArrayAlloc(kernels->solution1->n - 1); // Images of each kernel to return - psVector *fakeSolution = psVectorAlloc(kernels->solution1->n, PS_TYPE_F64); // Fake solution vector - psVectorInit(fakeSolution, 0.0); - - for (int i = 0; i < kernels->solution1->n - 1; i++) { - fakeSolution->data.F64[i] = kernels->solution1->data.F64[i]; - images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual); - fakeSolution->data.F64[i] = 0.0; - } - - psFree(fakeSolution); + psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution of interest + psVector *backup = psVectorCopy(NULL, solution, PS_TYPE_F64); // Backup version + + int num = wantDual ? solution->n - 1 : solution->n; // Number of kernel basis functions + + psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial + psArray *images = psArrayAlloc(num + 1); // Images of each kernel to return + + // The whole kernel + { + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual); // The appropriate kernel + images->data[0] = psMemIncrRefCounter(kernel->image); + psFree(kernel); + } + + // The parts + psVectorInit(solution, 0.0); + for (int i = 0; i < num; i++) { + solution->data.F64[i] = backup->data.F64[i]; + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, false, wantDual); // The appropriate kernel +#if 0 + int size = kernels->size; + double sum = 0.0; + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += kernel->kernel[v][u]; + } + } + fprintf(stderr, "Kernel %d: %lf\n", i, sum); +#endif + images->data[i + 1] = psMemIncrRefCounter(kernel->image); + psFree(kernel); + solution->data.F64[i] = 0.0; + } + psFree(polyValues); + psVectorCopy(solution, backup, PS_TYPE_F64); + psFree(backup); return images; @@ -1129,14 +1189,8 @@ if (!out1->image) { out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - // XXX if (threaded) { - // XXX psMutexInit(out1->image); - // XXX } } if (ro1->variance) { if (!out1->variance) { out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - // XXX if (threaded) { - // XXX psMutexInit(out1->variance); - // XXX } } psImageInit(out1->variance, 0.0); @@ -1146,14 +1200,8 @@ if (!out2->image) { out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - // XXX if (threaded) { - // XXX psMutexInit(out2->image); - // XXX } } if (ro2->variance) { if (!out2->variance) { out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - // XXX if (threaded) { - // XXX psMutexInit(out2->variance); - // XXX } } psImageInit(out2->variance, 0.0); @@ -1162,7 +1210,4 @@ psImage *convMask = NULL; // Convolved mask image (common to inputs 1 and 2) if (subMask) { - // XXX if (threaded) { - // XXX psMutexInit(subMask); - // XXX } if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { if (!out1->mask) { @@ -1188,16 +1233,12 @@ psImage *sys1 = NULL, *sys2 = NULL; // Systematic error images +#ifdef USE_SYS_ERR if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { sys1 = subtractionSysErrImage(ro1->image, sysError); - // XXX if (threaded && sys1) { - // XXX psMutexInit(sys1); - // XXX } } if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { sys2 = subtractionSysErrImage(ro2->image, sysError); - // XXX if (threaded && sys2) { - // XXX psMutexInit(sys2); - // XXX } - } + } +#endif int size = kernels->size; // Half-size of kernel @@ -1248,26 +1289,7 @@ psArrayAdd(args, 1, (pmReadout*)ro1); // Casting away const psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const - // Since adding to the array can impact the reference count, we need to lock - // XXX if (sys1) { - // XXX psMutexLock(sys1); - // XXX } psArrayAdd(args, 1, sys1); - // XXX if (sys1) { - // XXX psMutexUnlock(sys1); - // XXX } - // XXX if (sys2) { - // XXX psMutexLock(sys2); - // XXX } psArrayAdd(args, 1, sys2); - // XXX if (sys2) { - // XXX psMutexUnlock(sys2); - // XXX } - // XXX if (subMask) { - // XXX psMutexLock(subMask); - // XXX } psArrayAdd(args, 1, subMask); - // XXX if (subMask) { - // XXX psMutexUnlock(subMask); - // XXX } PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(args, maskPoor, PS_TYPE_IMAGE_MASK); @@ -1308,14 +1330,4 @@ psFree(job); } - - // XXX if (subMask) { - // XXX psMutexDestroy(subMask); - // XXX } - // XXX if (sys1) { - // XXX psMutexDestroy(sys1); - // XXX } - // XXX if (sys2) { - // XXX psMutexDestroy(sys2); - // XXX } } psImageConvolveSetThreads(oldThreads); Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25998) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25999) @@ -16,4 +16,5 @@ #define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image +//#define TESTING bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header, @@ -117,5 +118,5 @@ -#if 0 +#ifdef TESTING // Generate images of the kernel components { @@ -128,5 +129,20 @@ } psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, false); - psFits *kernelFile = psFitsOpen("kernels.fits", "w"); + psFits *kernelFile = psFitsOpen("kernels1.fits", "w"); + (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL); + psFitsClose(kernelFile); + psFree(kernelImages); + psFree(header); + } + if (kernels->solution2) { + psMetadata *header = psMetadataAlloc(); // Header + for (int i = 0; i < kernels->solution2->n; i++) { + psString name = NULL; // Header keyword + psStringAppend(&name, "SOLN%04d", i); + psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, kernels->solution2->data.F64[i]); + psFree(name); + } + psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, true); + psFits *kernelFile = psFitsOpen("kernels2.fits", "w"); (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL); psFitsClose(kernelFile); Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 25998) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 25999) @@ -15,7 +15,8 @@ #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? + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -23,127 +24,414 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Calculate the sum over a stamp product -static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication - const psKernel *image2, // Second image in multiplication - const psKernel *variance, // Variance image - int footprint // (Half-)Size of stamp - ) +// 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 + ) { - double sum = 0.0; // Sum of the image products + // (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 value = image1->kernel[y][x] * image2->kernel[y][x]; + double in = input->kernel[y][x]; + double ref = reference->kernel[y][x]; + double ir = in * ref; + double rr = PS_SQR(ref); + double one = 1.0; #ifdef USE_VARIANCE - value /= variance->kernel[y][x]; -#endif - sum += value; - } - } - return sum; -} - -// Calculate a single element of the least-squares matrix, with the polynomial expansions in one direction -static inline bool calculateMatrixElement1(psImage *matrix, // Matrix to calculate - int i, int j, // Coordinates of element - const psKernel *image1, // First image in multiplication - const psKernel *image2, // Second image in multiplication - const psKernel *variance, // Variance image - const psImage *polyValues, // Spatial polynomial values - int numKernels, // Number of kernel basis functions - int footprint, // (Half-)Size of stamp - int spatialOrder, // Maximum order of spatial variation - bool symmetric // Is the matrix symmetric? - ) -{ - double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products - if (!isfinite(sum)) { - return false; - } - - // Generate the pseudo-convolutions from the spatial polynomial terms - for (int iyOrder = 0, iIndex = i; iyOrder <= spatialOrder; iyOrder++) { - for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { - double convPoly = sum * polyValues->data.F64[iyOrder][ixOrder]; - - assert(iIndex < matrix->numRows && j < matrix->numCols); - - matrix->data.F64[iIndex][j] = convPoly; - if (symmetric) { - - assert(iIndex < matrix->numCols && j < matrix->numRows); - - matrix->data.F64[j][iIndex] = convPoly; - } - } - } + float varVal = 1.0 / variance->kernel[y][x]; + rr *= varVal; + ir *= varVal; + in *= varVal; + ref *= varVal; + one *= varVal; +#endif + sumRR += rr; + sumIR += ir; + sumR += ref; + sumI += in; + sum1 += one; + } + } + 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; } -// Calculate a single element of the least-squares matrix, with the polynomial expansions in both directions -static inline bool calculateMatrixElement2(psImage *matrix, // Matrix to calculate - int i, int j, // Coordinates of element - const psKernel *image1, // First image in multiplication - const psKernel *image2, // Second image in multiplication - const psKernel *variance, // Variance image - const psImage *polyValues, // Spatial polynomial values - int numKernels, // Number of kernel basis functions - int footprint, // (Half-)Size of stamp - int spatialOrder, // Maximum order of spatial variation - bool symmetric // Is the matrix symmetric? - ) +// 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 *image1, // Image 1 + const psKernel *image2, // Image 2 + const psKernel *variance, // Variance image + const psArray *convolutions1, // Convolutions of image 1 for each kernel + const psArray *convolutions2, // Convolutions of image 2 for each kernel + const pmSubtractionKernels *kernels, // Kernels + const psImage *polyValues, // Spatial polynomial values + int footprint // (Half-)Size of stamp + ) { - double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products - if (!isfinite(sum)) { - return false; - } - - // Generate the pseudo-convolutions from the spatial polynomial terms - for (int iyOrder = 0, iIndex = i; iyOrder <= spatialOrder; iyOrder++) { - for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { + // A_ij = A_i A_j + // B_ij = B_i B_j + // C_ij = A_i B_j + // d_i = A_i I_2 + // e_i = B_i I_2 + + // A_i = I_1 * k_i + // B_i = I_2 * k_i + + // Background: A_i = 1.0 + // Normalisation: A_i = I_1 + + 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 - for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) { - for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex += numKernels) { - double convPoly = sum * iPoly * polyValues->data.F64[jyOrder][jxOrder]; - - assert(iIndex < matrix->numRows && jIndex < matrix->numCols); - - matrix->data.F64[iIndex][jIndex] = convPoly; - if (symmetric) { - - assert(iIndex < matrix->numCols && jIndex < matrix->numRows); - - matrix->data.F64[jIndex][iIndex] = convPoly; - } - } - } - } - } + 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 *iConv1 = convolutions1->data[i]; // Convolution 1 for index i + psKernel *iConv2 = convolutions2->data[i]; // Convolution 2 for index i + for (int j = i; j < numKernels; j++) { + psKernel *jConv1 = convolutions1->data[j]; // Convolution 1 for index j + psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j + + double sumAA = 0.0; // Sum of convolution products for matrix A + double sumBB = 0.0; // Sum of convolution products for matrix B + double sumAB = 0.0; // Sum of convolution products for matrix C + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x]; + double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; +#ifdef USE_VARIANCE + aa /= variance->kernel[y][x]; + bb /= variance->kernel[y][x]; + ab /= variance->kernel[y][x]; +#endif + sumAA += aa; + sumBB += bb; + sumAB += ab; + } + } + + // 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 aa = sumAA * poly2[iTerm][jTerm]; + double bb = sumBB * poly2[iTerm][jTerm]; + double ab = sumAB * poly2[iTerm][jTerm]; + matrix1->data.F64[iIndex][jIndex] = aa; + matrix1->data.F64[jIndex][iIndex] = aa; + matrix2->data.F64[iIndex][jIndex] = bb; + matrix2->data.F64[jIndex][iIndex] = bb; + matrixX->data.F64[iIndex][jIndex] = ab; + } + } + } + for (int j = 0; j < i; j++) { + psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j + double sumAB = 0.0; // Sum of convolution products for matrix C + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; +#ifdef USE_VARIANCE + ab /= variance->kernel[y][x]; +#endif + sumAB += ab; + } + } + + // 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 ab = sumAB * poly2[iTerm][jTerm]; + matrixX->data.F64[iIndex][jIndex] = ab; + } + } + } + + double sumAI2 = 0.0; // Sum of A.I_2 products (for vector 1) + double sumBI2 = 0.0; // Sum of B.I_2 products (for vector 2) + double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix 1, normalisation) + double sumA = 0.0; // Sum of A (for matrix 1, background) + double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix X, normalisation) + double sumB = 0.0; // Sum of B products (for matrix X, background) + double sumI2 = 0.0; // Sum of I_2 (for vector 1, background) + double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector 1, normalisation) + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float a = iConv1->kernel[y][x]; + float b = iConv2->kernel[y][x]; + float i1 = image1->kernel[y][x]; + float i2 = image2->kernel[y][x]; + + double ai2 = a * i2; + double bi2 = b * i2; + double ai1 = a * i1; + double bi1 = b * i1; + double i1i2 = i1 * i2; + +#ifdef USE_VARIANCE + float varVal = 1.0 / variance->kernel[y][x]; + ai2 *= varVal; + bi2 *= varVal; + ai1 *= varVal; + bi1 *= varVal; + i1i2 *= varVal; + a *= varVal; + b *= varVal; + i2 *= varVal; +#endif + + sumAI2 += ai2; + sumBI2 += bi2; + sumAI1 += ai1; + sumA += a; + sumBI1 += bi1; + sumB += b; + sumI2 += i2; + sumI1I2 += i1i2; + } + } + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + double ai2 = sumAI2 * poly[iTerm]; + double bi2 = sumBI2 * poly[iTerm]; + double ai1 = sumAI1 * poly[iTerm]; + double a = sumA * poly[iTerm]; + double bi1 = sumBI1 * poly[iTerm]; + double b = sumB * poly[iTerm]; + + matrix1->data.F64[iIndex][normIndex] = ai1; + matrix1->data.F64[normIndex][iIndex] = ai1; + matrix1->data.F64[iIndex][bgIndex] = a; + matrix1->data.F64[bgIndex][iIndex] = a; + vector1->data.F64[iIndex] = ai2; + vector2->data.F64[iIndex] = bi2; + matrixX->data.F64[iIndex][normIndex] = bi1; + matrixX->data.F64[iIndex][bgIndex] = b; + } + } + + double sumI1 = 0.0; // Sum of I_1 (for matrix 1, background-normalisation) + double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix 1, normalisation-normalisation) + double sum1 = 0.0; // Sum of 1 (for matrix 1, background-background) + double sumI2 = 0.0; // Sum of I_2 (for vector 1, background) + double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector 1, normalisation) + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float i1 = image1->kernel[y][x]; + float i2 = image2->kernel[y][x]; + + double i1i1 = i1 * i1; + double one = 1.0; + double i1i2 = i1 * i2; + +#ifdef USE_VARIANCE + float varVal = 1.0 / variance->kernel[y][x]; + i1 *= varVal; + i1i1 *= varVal; + one *= varVal; + i2 *= varVal; + i1i2 *= varVal; +#endif + + sumI1 += i1; + sumI1I1 += i1i1; + sum1 += one; + sumI2 += i2; + sumI1I2 += i1i2; + } + } + matrix1->data.F64[bgIndex][normIndex] = sumI1; + matrix1->data.F64[normIndex][bgIndex] = sumI1; + matrix1->data.F64[normIndex][normIndex] = sumI1I1; + matrix1->data.F64[bgIndex][bgIndex] = sum1; + vector1->data.F64[bgIndex] = sumI2; + vector1->data.F64[normIndex] = sumI1I2; + return true; } -// Calculate the square part of the matrix derived from multiplying convolutions -static bool calculateMatrixSquare(psImage *matrix, // Matrix to calculate - const psArray *convolutions1, // Convolutions for element 1 - const psArray *convolutions2, // Convolutions for element 2 - const psKernel *variance, // Variance image - const psImage *polyValues, // Polynomial values - int numKernels, // Number of kernel basis functions - int spatialOrder, // Order of spatial variation - int footprint // Half-size of stamp +// Merge dual matrices and vectors into single matrix equation +// Have: Aa = Ct.b + d +// Have: Ca = Bb + e +// Set: F = ( A -Ct ; C -B ) +// Set: g = ( a ; b ) +// Set: h = ( d ; e ) +// So that we combine the above two equations: Fg = h +static bool calculateEquationDual(psImage **outMatrix, + psVector **outVector, + const psImage *sumMatrix1, + const psImage *sumMatrix2, + const psImage *sumMatrixX, + const psVector *sumVector1, + const psVector *sumVector2 ) { - bool symmetric = (convolutions1 == convolutions2 ? true : false); // Is matrix symmetric? - - for (int i = 0; i < numKernels; i++) { - psKernel *iConv = convolutions1->data[i]; // Convolution for i-th element - - for (int j = (symmetric ? i : 0); j < numKernels; j++) { - psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element - - if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, variance, polyValues, numKernels, - footprint, spatialOrder, symmetric)) { - psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j); - return false; - } + psAssert(sumMatrix1 && sumMatrix2 && sumMatrixX, "Require input matrices"); + psAssert(sumVector1 && sumVector2, "Require input vectors"); + int num1 = sumVector1->n; // Number of parameters in first set + int num2 = sumVector2->n; // Number of parameters in second set + int num = num1 + num2; // Number of parameters in new set + + psAssert(sumMatrix1->type.type == PS_TYPE_F64 && + sumMatrix2->type.type == PS_TYPE_F64 && + sumMatrixX->type.type == PS_TYPE_F64 && + sumVector1->type.type == PS_TYPE_F64 && + sumVector2->type.type == PS_TYPE_F64, + "Require input type is F64"); + + psAssert(outMatrix, "Require output matrix"); + psAssert(outVector, "Require output vector"); + if (!*outMatrix) { + *outMatrix = psImageAlloc(num, num, PS_TYPE_F64); + } + if (!*outVector) { + *outVector = psVectorAlloc(num, PS_TYPE_F64); + } + psImage *matrix = *outMatrix; + psVector *vector = *outVector; + + psAssert(sumMatrix1->numCols == num1 && sumMatrix1->numRows == num1, "Require size NxN"); + psAssert(sumMatrix2->numCols == num2 && sumMatrix2->numRows == num2, "Require size MxM"); + psAssert(sumMatrixX->numCols == num1 && sumMatrixX->numRows == num2, "Require size MxN"); + + memcpy(vector->data.F64, sumVector1->data.F64, num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + memcpy(&vector->data.F64[num1], sumVector2->data.F64, num2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + + for (int i = 0; i < num1; i++) { + memcpy(matrix->data.F64[i], sumMatrix1->data.F64[i], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + for (int j = 0, k = num1; j < num2; j++, k++) { + matrix->data.F64[i][k] = - sumMatrixX->data.F64[j][i]; + } + } + for (int i1 = 0, i2 = num1; i1 < num2; i1++, i2++) { + memcpy(matrix->data.F64[i2], sumMatrixX->data.F64[i1], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + for (int j = 0, k = num1; j < num2; j++, k++) { + matrix->data.F64[i2][k] = - sumMatrix2->data.F64[i1][j]; } } @@ -152,263 +440,10 @@ } -// Calculate least-squares matrix and vector -static bool calculateMatrix(psImage *matrix, // Matrix to calculate - const pmSubtractionKernels *kernels, // Kernel components - const psArray *convolutions, // Convolutions of source with kernels - const psKernel *input, // Input stamp, or NULL - const psKernel *variance, // Variance stamp - const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - bool normAndBG // Calculate normalisation and background terms? - ) -{ - int numKernels = kernels->num; // Number of kernel components - int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation - int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms - int bgOrder = kernels->bgOrder; // Maximum order of background fit - int numBackground = normAndBG ? PM_SUBTRACTION_POLYTERMS(bgOrder) : 0; // Number of background terms - int numTerms = numKernels * numSpatial + (normAndBG ? 1 + numBackground : 0); // Total number of terms - assert(matrix); - assert(matrix->numCols == matrix->numRows); - assert(matrix->numCols == numTerms); - assert(convolutions && convolutions->n == numKernels); - assert(polyValues); - assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image - - // Square part of the matrix (convolution-convolution products) - if (!calculateMatrixSquare(matrix, convolutions, convolutions, variance, polyValues, numKernels, - spatialOrder, footprint)) { - return false; - } - - // XXX To support higher-order background model than simply constant, the below code needs to be updated. - if (normAndBG) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - - for (int i = 0; i < numKernels; i++) { - psKernel *conv = convolutions->data[i]; // Convolution for i-th element - - // Normalisation-convolution terms - if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, variance, polyValues, numKernels, - footprint, spatialOrder, true)) { - psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); - return false; - } - - // Background-convolution terms - double sumC = 0.0; // Sum of the convolution - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double value = conv->kernel[y][x]; -#ifdef USE_VARIANCE - value /= variance->kernel[y][x]; -#endif - sumC += value; - } - } - if (!isfinite(sumC)) { - psTrace("psModules.imcombine", 2, "Bad sumC at %d", i); - return false; - } - - for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { - for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { - double value = sumC * polyValues->data.F64[yOrder][xOrder]; - matrix->data.F64[index][bgIndex] = value; - matrix->data.F64[bgIndex][index] = value; - } - } - } - - // Background only, normalisation only, and background-normalisation terms - double sum1 = 0.0; // Sum of the weighting - double sumI = 0.0; // Sum of the input - double sumII = 0.0; // Sum of the input squared - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double invNoise2 = 1.0; -#ifdef USE_VARIANCE - invNoise2 /= variance->kernel[y][x]; -#endif - double value = input->kernel[y][x] * invNoise2; - sumI += value; - sumII += value * input->kernel[y][x]; - sum1 += invNoise2; - } - } - if (!isfinite(sumI)) { - psTrace("psModules.imcombine", 2, "Bad sumI detected"); - return false; - } - if (!isfinite(sumII)) { - psTrace("psModules.imcombine", 2, "Bad sumII detected"); - return false; - } - if (!isfinite(sum1)) { - psTrace("psModules.imcombine", 2, "Bad sum1 detected"); - return false; - } - matrix->data.F64[normIndex][normIndex] = sumII; - matrix->data.F64[bgIndex][bgIndex] = sum1; - matrix->data.F64[normIndex][bgIndex] = sumI; - matrix->data.F64[bgIndex][normIndex] = sumI; - } - - return true; -} - - -// Calculate least-squares matrix and vector -static bool calculateVector(psVector *vector, // Vector to calculate, or NULL - const pmSubtractionKernels *kernels, // Kernel components - const psArray *convolutions, // Convolutions of source with kernels - const psKernel *input, // Input stamp, or NULL if !normAndBG - const psKernel *target, // Target stamp - const psKernel *variance, // Variance stamp - const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - bool normAndBG // Calculate normalisation and background terms? - ) -{ - int numKernels = kernels->num; // Number of kernel components - int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation - int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms - int bgOrder = kernels->bgOrder; // Maximum order of background fit - int numBackground = normAndBG ? PM_SUBTRACTION_POLYTERMS(bgOrder) : 0; // Number of background terms - int numTerms = numKernels * numSpatial + (normAndBG ? 1 + numBackground : 0); // Total number of terms - assert(vector && vector->n == numTerms); - assert(convolutions && convolutions->n == numKernels); - assert(target); - assert(polyValues); - assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image - - // Convolution terms - for (int i = 0; i < numKernels; i++) { - psKernel *conv = convolutions->data[i]; // Convolution for i-th element - double sumTC = 0.0; // Sum of the target and convolution - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double value = target->kernel[y][x] * conv->kernel[y][x]; -#ifdef USE_VARIANCE - value /= variance->kernel[y][x]; -#endif - sumTC += value; - } - } - if (!isfinite(sumTC)) { - psTrace("psModules.imcombine", 2, "Bad sumTC at %d", i); - return false; - } - for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { - for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { - vector->data.F64[index] = sumTC * polyValues->data.F64[yOrder][xOrder]; - } - } - } - - if (normAndBG) { - // Background terms - double sumT = 0.0; // Sum of the target - double sumIT = 0.0; // Sum of the input-target product - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double value = target->kernel[y][x]; -#ifdef USE_VARIANCE - value /= variance->kernel[y][x]; -#endif - sumIT += value * input->kernel[y][x]; - sumT += value; - } - } - if (!isfinite(sumT)) { - psTrace("psModules.imcombine", 2, "Bad sumI detected"); - return false; - } - if (!isfinite(sumIT)) { - psTrace("psModules.imcombine", 2, "Bad sumIT detected"); - return false; - } - - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation term - vector->data.F64[normIndex] = sumIT; - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background term - vector->data.F64[bgIndex] = sumT; - } - - return true; -} - - - -// Calculate the cross-matrix, composed of convolutions of each image -// Note that the cross-matrix is NOT square -static bool calculateMatrixCross(psImage *matrix, // Matrix to calculate - const pmSubtractionKernels *kernels, // Kernel components - const psArray *convolutions1, // Convolutions of image 1 - const psArray *convolutions2, // Convolutions of image 2 - const psKernel *image1, // Image 1 stamp - const psKernel *variance, // Variance stamp - const psImage *polyValues, // Spatial polynomial values - int footprint // (Half-)Size of stamp - ) -{ - assert(matrix); - int numKernels = kernels->num; // Number of kernel components - int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation - int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial polynomial terms - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - int numCols = numKernels * numSpatial + 1 + numBackground; // Number of columns - int numRows = numKernels * numSpatial; // Number of rows - assert(matrix->numCols == numCols && matrix->numRows == numRows); - assert(convolutions1 && convolutions1->n == numKernels); - assert(convolutions2 && convolutions2->n == numKernels); - - int normIndex, bgIndex; // Indices in matrix for normalisation and background terms - PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); - - if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, variance, polyValues, numKernels, - spatialOrder, footprint)) { - return false; - } - - for (int i = 0; i < numKernels; i++) { - // Normalisation - psKernel *conv = convolutions2->data[i]; // Convolution - if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, variance, polyValues, numKernels, - footprint, spatialOrder, false)) { - psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); - return false; - } - - // Background - double sumC = 0.0; // Sum of the weighting - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double value = conv->kernel[y][x]; -#ifdef USE_VARIANCE - value /= variance->kernel[y][x]; -#endif - sumC += value; - } - } - if (!isfinite(sumC)) { - psTrace("psModules.imcombine", 2, "Bad sumC detected at %d", i); - return false; - } - for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { - for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { - matrix->data.F64[index][bgIndex] = sumC * polyValues->data.F64[yOrder][xOrder]; - } - } - } - - return true; -} - // Add in penalty term to least-squares vector -static bool calculatePenalty(psVector *vector, // Vector to which to add in penalty term - const pmSubtractionKernels *kernels // Kernel parameters +static bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term + psVector *vector, // Vector to which to add in penalty term + const pmSubtractionKernels *kernels, // Kernel parameters + float norm // Normalisation ) { @@ -423,5 +458,6 @@ for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { - vector->data.F64[index] -= penalties->data.F32[i]; + // Contribution to chi^2: a_i^2 P_i + matrix->data.F64[index][index] -= norm * penalties->data.F32[i]; } } @@ -582,14 +618,12 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: - 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); + status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, + stamp->variance, stamp->convolutions1, kernels, polyValues, + footprint); break; case PM_SUBTRACTION_MODE_2: - 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); + status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2, + stamp->variance, stamp->convolutions2, kernels, polyValues, + footprint); break; case PM_SUBTRACTION_MODE_DUAL: @@ -604,15 +638,8 @@ psVectorInit(stamp->vector2, NAN); #endif - status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, - stamp->variance, polyValues, footprint, true); - status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, - stamp->variance, polyValues, footprint, false); - status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, - stamp->convolutions2, stamp->image1, stamp->variance, polyValues, - footprint); - status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->variance, polyValues, footprint, true); - status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, - stamp->image2, stamp->variance, polyValues, footprint, false); + status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, + stamp->matrixX, stamp->image1, stamp->image2, stamp->variance, + stamp->convolutions1, stamp->convolutions2, kernels, polyValues, + footprint); break; default: @@ -629,5 +656,5 @@ #ifdef TESTING - if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { + { psString matrixName = NULL; psStringAppend(&matrixName, "matrix1_%d.fits", index); @@ -825,5 +852,8 @@ #endif - calculatePenalty(sumVector, kernels); +#if 0 + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); +#endif #ifdef TESTING @@ -909,160 +939,177 @@ } } - calculatePenalty(sumVector1, kernels); - calculatePenalty(sumVector2, kernels); - - // Pure matrix operations - - // A * a = Ct * b + d - // C * a = B * b + e - // - // a = (Ct * Bi * C - A)i (Ct * Bi * e - d) - // b = Bi * (C * a - e) - psVector *a = psVectorRecycle(kernels->solution1, numParams, PS_TYPE_F64); - psVector *b = psVectorRecycle(kernels->solution2, numParams2, PS_TYPE_F64); -#ifdef TESTING - psVectorInit(a, NAN); - psVectorInit(b, NAN); -#endif - psImage *A = sumMatrix1; - psImage *B = sumMatrix2; - psImage *C = sumMatrixX; - psVector *d = sumVector1; - psVector *e = sumVector2; - - assert(a->n == numParams); - assert(b->n == numParams2); - assert(A->numRows == numParams && A->numCols == numParams); - assert(B->numRows == numParams2 && B->numCols == numParams2); - assert(C->numRows == numParams2 && C->numCols == numParams); - assert(d->n == numParams); - assert(e->n == numParams2); - - psImage *Bi = psMatrixInvert(NULL, B, NULL); - assert(Bi->numRows == numParams2 && Bi->numCols == numParams2); - psImage *Ct = psMatrixTranspose(NULL, C); - assert(Ct->numRows == numParams && Ct->numCols == numParams2); - - psImage *BiC = psMatrixMultiply(NULL, Bi, C); - assert(BiC->numRows == numParams2 && BiC->numCols == numParams); - psImage *CtBi = psMatrixMultiply(NULL, Ct, Bi); - assert(CtBi->numRows == numParams && CtBi->numCols == numParams2); - - psImage *CtBiC = psMatrixMultiply(NULL, Ct, BiC); - assert(CtBiC->numRows == numParams && CtBiC->numCols == numParams); - - psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A); - assert(F->numRows == numParams && F->numCols == numParams); - float det = 0.0; - psImage *Fi = psMatrixInvert(NULL, F, &det); - assert(Fi->numRows == numParams && Fi->numCols == numParams); - psTrace("psModules.imcombine", 4, "Determinant of F: %f\n", det); - - psVector *g = psVectorAlloc(numParams, PS_TYPE_F64); -#ifdef TESTING - psVectorInit(g, NAN); -#endif - assert(CtBi->numRows == numParams && CtBi->numCols == numParams2); - assert(e->n == numParams2); - assert(d->n == numParams); - for (int i = 0; i < numParams; i++) { - double value = 0.0; - for (int j = 0; j < numParams2; j++) { - value += CtBi->data.F64[i][j] * e->data.F64[j]; - } - g->data.F64[i] = value - d->data.F64[i]; - } - - assert(Fi->numRows == numParams && Fi->numCols == numParams); - assert(g->n == numParams); - for (int i = 0; i < numParams; i++) { - double value = 0.0; - for (int j = 0; j < numParams; j++) { - value += Fi->data.F64[i][j] * g->data.F64[j]; - } - a->data.F64[i] = value; - } - - psVector *h = psVectorAlloc(numParams2, PS_TYPE_F64); -#ifdef TESTING - psVectorInit(h, NAN); -#endif - assert(C->numRows == numParams2 && C->numCols == numParams); - assert(a->n == numParams); - assert(e->n == numParams2); - for (int i = 0; i < numParams2; i++) { - double value = 0.0; - for (int j = 0; j < numParams; j++) { - value += C->data.F64[i][j] * a->data.F64[j]; - } - h->data.F64[i] = value - e->data.F64[i]; - } - - assert(Bi->numRows == numParams2 && Bi->numCols == numParams2); - assert(h->n == numParams2); - for (int i = 0; i < numParams2; i++) { - double value = 0.0; - for (int j = 0; j < numParams2; j++) { - value += Bi->data.F64[i][j] * h->data.F64[j]; - } - b->data.F64[i] = value; - } - - -#if 0 - for (int i = 0; i < numParams; i++) { - double aVal1 = 0.0, bVal1 = 0.0; - for (int j = 0; j < numParams2; j++) { - aVal1 += A->data.F64[i][j] * a->data.F64[j]; - bVal1 += Ct->data.F64[i][j] * b->data.F64[j]; - } - bVal1 += d->data.F64[i]; - for (int j = numParams2; j < numParams; j++) { - aVal1 += A->data.F64[i][j] * a->data.F64[j]; - } - printf("%d: %lf\n", i, aVal1 - bVal1); - } - - for (int i = 0; i < numParams2; i++) { - double aVal2 = 0.0, bVal2 = 0.0; - for (int j = 0; j < numParams2; j++) { - aVal2 += C->data.F64[i][j] * a->data.F64[j]; - bVal2 += B->data.F64[i][j] * b->data.F64[j]; - } - bVal2 += e->data.F64[i]; - for (int j = numParams2; j < numParams; j++) { - aVal2 += C->data.F64[i][j] * a->data.F64[j]; - } - printf("%d: %lf\n", i, aVal2 - bVal2); - } -#endif + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + calculatePenalty(sumMatrix1, sumVector1, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]); + calculatePenalty(sumMatrix2, sumVector2, kernels, -sumMatrix1->data.F64[bgIndex][bgIndex]); + + psImage *sumMatrix = NULL; // Combined matrix + psVector *sumVector = NULL; // Combined vector + calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2, + sumMatrixX, sumVector1, sumVector2); #ifdef TESTING { - psFits *fits = psFitsOpen("sumMatrix1.fits", "w"); - psFitsWriteImage(fits, NULL, sumMatrix1, 0, NULL); + psFits *fits = psFitsOpen("sumMatrix.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL); psFitsClose(fits); } { - psFits *fits = psFitsOpen("sumMatrix2.fits", "w"); - psFitsWriteImage(fits, NULL, sumMatrix2, 0, NULL); + psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); + psFits *fits = psFitsOpen("sumVector.fits", "w"); + for (int i = 0; i < numParams + numParams2; i++) { + image->data.F64[0][i] = sumVector->data.F64[i]; + } + psFitsWriteImage(fits, NULL, image, 0, NULL); + psFree(image); psFitsClose(fits); } +#endif + + psVector *solution = NULL; // Solution to equation! { - psFits *fits = psFitsOpen("sumMatrixX.fits", "w"); - psFitsWriteImage(fits, NULL, sumMatrixX, 0, NULL); + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector); + if (!solution) { + psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n"); + psFree(sumMatrix); + psFree(sumVector); + return NULL; + } + + int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations + int numKernels = kernels->num; // Number of kernel basis functions + + // Remove a kernel basis for image 1 from the equation +#define MASK_BASIS_1(INDEX) \ + { \ + for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \ + for (int k = 0; k < numParams2; k++) { \ + sumMatrix1->data.F64[k][index] = 0.0; \ + sumMatrix1->data.F64[index][k] = 0.0; \ + sumMatrixX->data.F64[k][index] = 0.0; \ + } \ + sumMatrix1->data.F64[bgIndex][index] = 0.0; \ + sumMatrix1->data.F64[index][bgIndex] = 0.0; \ + sumMatrix1->data.F64[normIndex][index] = 0.0; \ + sumMatrix1->data.F64[index][normIndex] = 0.0; \ + sumMatrix1->data.F64[index][index] = 1.0; \ + sumVector1->data.F64[index] = 0.0; \ + } \ + } + + // Remove a kernel basis for image 2 from the equation +#define MASK_BASIS_2(INDEX) \ + { \ + for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \ + for (int k = 0; k < numParams2; k++) { \ + sumMatrix2->data.F64[k][index] = 0.0; \ + sumMatrix2->data.F64[index][k] = 0.0; \ + sumMatrixX->data.F64[index][k] = 0.0; \ + } \ + sumMatrix2->data.F64[index][index] = 1.0; \ + sumMatrixX->data.F64[index][normIndex] = 0.0; \ + sumMatrixX->data.F64[index][bgIndex] = 0.0; \ + sumVector2->data.F64[index] = 0.0; \ + } \ + } + + #define TOL 1.0e-5 + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + double norm = solution->data.F64[normIndex]; // Normalisation + double thresh = norm * TOL; // Threshold for low parameters + for (int i = 0; i < numKernels; i++) { + // Getting 0th order parameter value. In the presence of spatial variation, the actual value + // of the parameter will vary over the image. We are in effect getting the value in the + // centre of the image. If we use different polynomial functions (e.g., Chebyshev), we may + // have to change this to properly determine the value of the parameter at the centre. + double param1 = solution->data.F64[i], + param2 = solution->data.F64[numParams + i]; // Parameters of interest + bool mask1 = false, mask2 = false; // Masked the parameter? + if (fabs(param1) < thresh) { + psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i); + MASK_BASIS_1(i); + mask1 = true; + } + if (fabs(param2) < thresh) { + psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i); + MASK_BASIS_2(i); + mask2 = true; + } + + if (!mask1 && !mask2) { + if (fabs(param1) < fabs(param2)) { + psTrace("psModules.imcombine", 7, "Parameter %d: 1 < 2\n", i); + MASK_BASIS_1(i); + } else { + psTrace("psModules.imcombine", 7, "Parameter %d: 2 < 1\n", i); + MASK_BASIS_2(i); + } + } + } + } + + calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2, + sumMatrixX, sumVector1, sumVector2); + +#ifdef TESTING + { + psFits *fits = psFitsOpen("sumMatrixFix.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL); psFitsClose(fits); } { - psFits *fits = psFitsOpen("sumFinverse.fits", "w"); - psFitsWriteImage(fits, NULL, Fi, 0, NULL); + psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); + psFits *fits = psFitsOpen("sumVectorFix.fits", "w"); + for (int i = 0; i < numParams + numParams2; i++) { + image->data.F64[0][i] = sumVector->data.F64[i]; + } + psFitsWriteImage(fits, NULL, image, 0, NULL); + psFree(image); psFitsClose(fits); } #endif - kernels->solution1 = a; - kernels->solution2 = b; - - // XXXXX Free temporary matrices and vectors + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector); + if (!solution) { + psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n"); + psFree(sumMatrix); + psFree(sumVector); + return NULL; + } + + psFree(sumMatrix1); + psFree(sumMatrix2); + psFree(sumMatrixX); + psFree(sumVector1); + psFree(sumVector2); + + psFree(sumMatrix); + psFree(sumVector); + +#ifdef TESTING + { + psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); + psFits *fits = psFitsOpen("solnVector.fits", "w"); + for (int i = 0; i < numParams + numParams2; i++) { + image->data.F64[0][i] = solution->data.F64[i]; + } + psFitsWriteImage(fits, NULL, image, 0, NULL); + psFree(image); + psFitsClose(fits); + } +#endif + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc(numParams, PS_TYPE_F64); + } + if (!kernels->solution2) { + kernels->solution2 = psVectorAlloc(numParams2, PS_TYPE_F64); + } + + memcpy(kernels->solution1->data.F64, solution->data.F64, numParams * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + memcpy(kernels->solution2->data.F64, &solution->data.F64[numParams], + numParams2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + + psFree(solution); } Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 25998) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 25999) @@ -81,4 +81,5 @@ kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew); kernels->inner = start; + kernels->num += numNew; // Generate a set of kernels for each (u,v) @@ -94,9 +95,15 @@ kernels->v->data.S32[index] = v; kernels->preCalc->data[index] = NULL; - kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v)); + kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); } } + + kernels->widths->n = start + numNew; + kernels->u->n = start + numNew; + kernels->v->n = start + numNew; + kernels->preCalc->n = start + numNew; + kernels->penalties->n = start + numNew; return true; @@ -140,7 +147,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,5 +157,6 @@ for (int u = -size, x = 0; u <= size; u++, x++) { double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel - moment += value * (PS_SQR(u) + PS_SQR(v)); + kernel->kernel[v][u] = value; + moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); } } @@ -164,6 +173,19 @@ 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); } + + +#if 0 + double sum = 0.0; // Sum of kernel component + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += kernel->kernel[v][u]; + } + } + fprintf(stderr, "%d sum: %lf\n", index, sum); +#endif kernels->widths->data.F32[index] = fwhms->data.F32[i]; @@ -456,8 +478,7 @@ PS_ASSERT_INT_LESS_THAN(inner, size, NULL); - // 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 + kernels->type = PM_SUBTRACTION_KERNEL_GUNK; psStringPrepend(&kernels->description, "GUNK="); psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); @@ -577,5 +598,5 @@ poly->data.F32[j] = polyVal; norm += polyVal; - moment += polyVal * (PS_SQR(u) + PS_SQR(v)); + moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v)); psVectorExtend(uCoords, RINGS_BUFFER, 1); @@ -603,5 +624,5 @@ psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); } -// moment /= norm; + moment /= norm; }