Index: /branches/pap/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 25833) +++ /branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 25834) @@ -964,5 +964,5 @@ } -#if 0 +#if 1 psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) { @@ -972,15 +972,19 @@ 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]; + 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 + + psArray *images = psArrayAlloc(num); // Images of each kernel to return + psVectorInit(solution, 0.0); + + for (int i = 0; i < num; i++) { + solution->data.F64[i] = backup->data.F64[i]; images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual); - fakeSolution->data.F64[i] = 0.0; - } - - psFree(fakeSolution); + solution->data.F64[i] = 0.0; + } + psVectorCopy(solution, backup, PS_TYPE_F64); + psFree(backup); return images; Index: /branches/pap/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25833) +++ /branches/pap/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25834) @@ -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: /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 25833) +++ /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 25834) @@ -19,9 +19,11 @@ #define USE_VARIANCE // Include variance in equation? +//#define OLD_FUNCTIONS // Use old functions + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Private (file-static) functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -#if 1 +#ifndef OLD_FUNCTIONS // Calculate the least-squares matrix and vector static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated @@ -136,4 +138,5 @@ double ir = in * ref; double rr = PS_SQR(ref); + double one = 1.0; #ifdef USE_VARIANCE float varVal = 1.0 / variance->kernel[y][x]; @@ -142,4 +145,5 @@ in *= varVal; ref *= varVal; + one *= varVal; #endif sumRR += rr; @@ -147,5 +151,5 @@ sumR += ref; sumI += in; - sum1 += 1.0; + sum1 += one; } } @@ -159,5 +163,4 @@ } -#if 0 // Calculate the least-squares matrix and vector for dual convolution static bool calculateDualMatrixVector(psImage *matrix1, // Least-squares matrix, updated @@ -166,8 +169,9 @@ 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 *image1, // Image 1 + const psKernel *image2, // Image 2 const psKernel *variance, // Variance image - const psArray *convolutions, // Convolutions for each kernel + 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 @@ -175,9 +179,15 @@ ) { - // (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 + // 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 @@ -205,16 +215,26 @@ for (int i = 0; i < numKernels; i++) { - psKernel *iConv = convolutions->data[i]; // Convolution for index 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 *jConv = convolutions->data[j]; // Convolution for index j - - double sumCC = 0.0; // Sum of convolution products + 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 cc = iConv->kernel[y][x] * jConv->kernel[y][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 - cc /= variance->kernel[y][x]; -#endif - sumCC += cc; + aa /= variance->kernel[y][x]; + bb /= variance->kernel[y][x]; + ab /= variance->kernel[y][x]; +#endif + sumAA += aa; + sumBB += bb; + sumAB += ab; } } @@ -223,79 +243,139 @@ 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 + 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 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; + 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]; - ic *= varVal; - rc *= varVal; - c *= varVal; -#endif - sumIC += ic; - sumRC += rc; - sumC += c; + 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 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 + 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++) { - double in = input->kernel[y][x]; - double ref = reference->kernel[y][x]; - double ir = in * ref; - double rr = PS_SQR(ref); + 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]; - 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; + 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; } -#endif #endif @@ -321,4 +401,5 @@ } +#ifdef OLD_FUNCTIONS // 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 @@ -684,5 +765,5 @@ return true; } - +#endif // Add in penalty term to least-squares vector @@ -860,5 +941,5 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: -#if 0 +#ifdef OLD_FUNCTIONS status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, stamp->variance, polyValues, footprint, true); @@ -872,5 +953,5 @@ break; case PM_SUBTRACTION_MODE_2: -#if 0 +#ifdef OLD_FUNCTIONS status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, stamp->variance, polyValues, footprint, true); @@ -894,4 +975,5 @@ psVectorInit(stamp->vector2, NAN); #endif +#ifdef OLD_FUNCTIONS status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, stamp->variance, polyValues, footprint, true); @@ -905,4 +987,10 @@ status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, stamp->image2, stamp->variance, polyValues, footprint, false); +#else + 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); +#endif break; default: