Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 15799) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 15809) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.74 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-12-07 01:57:15 $ + * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-13 21:32:44 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -165,9 +165,6 @@ } - if (!wantDual) { - // Put in the normalisation component - kernel->kernel[0][0] += p_pmSubtractionSolutionNorm(kernels); - } - + // Put in the normalisation component + kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels)); return kernel; @@ -279,8 +276,9 @@ float background, // Background value to apply psRegion region, // Region to convolve - bool useFFT // Use FFT to convolve? + bool useFFT, // Use FFT to convolve? + bool wantDual // Want the dual convolution? ) { - *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, false); + *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual); if (weight) { *kernelWeight = varianceKernel(*kernelWeight, *kernelImage); @@ -559,4 +557,6 @@ stamp->x = 0; stamp->y = 0; + stamp->xNorm = NAN; + stamp->yNorm = NAN; stamp->status = PM_SUBTRACTION_STAMP_REJECTED; // Recalculate convolutions @@ -568,4 +568,11 @@ psFree(stamp->weight); stamp->image1 = stamp->image2 = stamp->weight = NULL; + psFree(stamp->matrix1); + psFree(stamp->matrix2); + psFree(stamp->matrixX); + stamp->matrix1 = stamp->matrix2 = stamp->matrixX = NULL; + psFree(stamp->vector1); + psFree(stamp->vector2); + stamp->vector1 = stamp->vector2 = NULL; } else { numGood++; @@ -651,5 +658,5 @@ PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false); - if (kernels->mode != PM_SUBTRACTION_MODE_1) { + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { PS_ASSERT_PTR_NON_NULL(out2, false); PS_ASSERT_PTR_NON_NULL(ro2, false); @@ -812,9 +819,9 @@ convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image, weight, - subMask, maskSource, kernels, polyValues, background, subRegion, useFFT); + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, false); if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, ro2->image, ro2->weight, - subMask, maskSource, kernels, polyValues, background, subRegion, useFFT); + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, true); } Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 15799) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 15809) @@ -58,7 +58,12 @@ for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { double convPoly = sum * polyValues->data.F64[iyOrder][ixOrder]; - // XXX Not sure I've got this the right way around + + 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; } @@ -93,6 +98,12 @@ 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; } @@ -139,5 +150,4 @@ const psArray *convolutions, // Convolutions of source with kernels const psKernel *input, // Input stamp, or NULL - const psKernel *target, // Target stamp, or NULL const psKernel *weight, // Weight stamp const psImage *polyValues, // Spatial polynomial values @@ -156,5 +166,4 @@ assert(matrix->numCols == 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 @@ -168,11 +177,11 @@ // XXX To support higher-order background model than simply constant, the below code needs to be updated. if (normAndBG) { - // Normalisation and background terms - int normIndex, bgIndex; // Indices in matrix for normalisation and background - PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); + 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 + + // Normalisation-convolution terms if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels, footprint, spatialOrder, true)) { @@ -181,5 +190,5 @@ } - // Background + // Background-convolution terms double sumC = 0.0; // Sum of the convolution for (int y = - footprint; y <= footprint; y++) { @@ -192,13 +201,15 @@ 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] = matrix->data.F64[bgIndex][index] = - sumC * polyValues->data.F64[yOrder][xOrder]; + double value = sumC * polyValues->data.F64[yOrder][xOrder]; + matrix->data.F64[index][bgIndex] = value; + matrix->data.F64[bgIndex][index] = value; } } } - // Normalisation only and background only terms in the matrix + // Background only, normalisation only, and background-normalisation terms double sum1 = 0.0; // Sum of the weighting double sumI = 0.0; // Sum of the input @@ -206,17 +217,17 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float invNoise2 = 1.0 / weight->kernel[y][x]; - float value = input->kernel[y][x] * invNoise2; + double invNoise2 = 1.0 / weight->kernel[y][x]; + double value = input->kernel[y][x] * invNoise2; sumI += value; - sumII += input->kernel[y][x] * 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(sumI)) { - psTrace("psModules.imcombine", 2, "Bad sumI detected"); return false; } @@ -239,9 +250,10 @@ const pmSubtractionKernels *kernels, // Kernel components const psArray *convolutions, // Convolutions of source with kernels - const psKernel *input, // Input stamp, or NULL - const psKernel *target, // Target stamp, or NULL + const psKernel *input, // Input stamp, or NULL if !normAndBG + const psKernel *target, // Target stamp const psKernel *weight, // Weight stamp const psImage *polyValues, // Spatial polynomial values - int footprint // (Half-)Size of stamp + int footprint, // (Half-)Size of stamp + bool normAndBG // Calculate normalisation and background terms? ) { @@ -250,11 +262,11 @@ int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms int bgOrder = kernels->bgOrder; // Maximum order of background fit - int numBackground = PM_SUBTRACTION_POLYTERMS(bgOrder); // Number of background terms - int numTerms = numKernels * numSpatial + 1 + numBackground; // Total number of terms + 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(input); assert(target); assert(polyValues); + assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image // Convolution terms @@ -278,27 +290,29 @@ } - // Normalisation and background terms - double sumIT = 0.0; // Sum of the target and input - double sumT = 0.0; // Sum of the target - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - float value = target->kernel[y][x] / weight->kernel[y][x]; - sumIT += value * input->kernel[y][x]; - sumT += value; - } - } - if (!isfinite(sumIT)) { - psTrace("psModules.imcombine", 2, "Bad sumIT detected"); - return false; - } - if (!isfinite(sumT)) { - psTrace("psModules.imcombine", 2, "Bad sumI detected"); - return false; - } - - int normIndex, bgIndex; // Indices in vector for normalisation and background terms - PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); - vector->data.F64[normIndex] = sumIT; - vector->data.F64[bgIndex] = sumT; + 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++) { + float value = target->kernel[y][x] / weight->kernel[y][x]; + 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; @@ -383,4 +397,7 @@ for (int yOrder = 0; yOrder <= order; yOrder++) { for (int xOrder = 0; xOrder <= order - yOrder; xOrder++, index += step) { + + assert(index < coeff->n); + sum += coeff->data.F64[index] * polyValues->data.F64[yOrder][xOrder]; } @@ -468,9 +485,36 @@ } +#ifdef TESTING + for (int j = 0; j < numKernels; j++) { + if (stamp->convolutions1) { + psString convName = NULL; + psStringAppend(&convName, "conv1_%03d_%03d.fits", i, j); + psFits *fits = psFitsOpen(convName, "w"); + psFree(convName); + psKernel *conv = stamp->convolutions1->data[j]; + psFitsWriteImage(fits, NULL, conv->image, 0, NULL); + psFitsClose(fits); + } + + if (stamp->convolutions2) { + psString convName = NULL; + psStringAppend(&convName, "conv2_%03d_%03d.fits", i, j); + psFits *fits = psFitsOpen(convName, "w"); + psFree(convName); + psKernel *conv = stamp->convolutions2->data[j]; + psFitsWriteImage(fits, NULL, conv->image, 0, NULL); + psFitsClose(fits); + } + } +#endif + polyValues = p_pmSubtractionPolynomial(polyValues, spatialOrder, stamp->xNorm, stamp->yNorm); - stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64); + + stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); + stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64); #ifdef TESTING - psVectorInit(stamp->vector, NAN); + psImageInit(stamp->matrix1, NAN); + psVectorInit(stamp->vector1, NAN); #endif @@ -479,42 +523,35 @@ case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: - stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); -#ifdef TESTING - psImageInit(stamp->matrix1, NAN); -#endif status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint, true); - status &= calculateVector(stamp->vector, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint); + stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint, true); break; case PM_SUBTRACTION_MODE_2: - stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); -#ifdef TESTING - psImageInit(stamp->matrix1, NAN); -#endif status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, - stamp->image1, stamp->weight, polyValues, footprint, true); - status &= calculateVector(stamp->vector, kernels, stamp->convolutions2, stamp->image2, - stamp->image1, stamp->weight, polyValues, footprint); + stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, + stamp->image1, stamp->weight, polyValues, footprint, true); break; case PM_SUBTRACTION_MODE_DUAL: - stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64); + stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); #ifdef TESTING - psImageInit(stamp->matrix1, NAN); psImageInit(stamp->matrix2, NAN); psImageInit(stamp->matrixX, NAN); -#endif - status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint, true); - status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, stamp->image2, - stamp->image1, stamp->weight, polyValues, footprint, false); - // XXX Not sure I'm passing the correct stamp->image to calculateMatrixCross: + psVectorInit(stamp->vector2, NAN); +#endif + status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, + stamp->weight, polyValues, footprint, true); + status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, + stamp->weight, polyValues, footprint, false); status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, - stamp->convolutions2, stamp->image2, stamp->weight, polyValues, + stamp->convolutions2, stamp->image1, stamp->weight, polyValues, footprint); - status &= calculateVector(stamp->vector, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint); + status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, + stamp->image2, stamp->weight, polyValues, footprint, false); break; default: @@ -540,8 +577,8 @@ matrixName = NULL; - psStringAppend(&matrixName, "vector_%d.fits", i); - psImage *dummy = psImageAlloc(stamp->vector->n, 1, PS_TYPE_F64); - memcpy(dummy->data.F64[0], stamp->vector->data.F64, - PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector->n); + psStringAppend(&matrixName, "vector1_%d.fits", i); + psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64); + memcpy(dummy->data.F64[0], stamp->vector1->data.F64, + PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n); matrixFile = psFitsOpen(matrixName, "w"); psFree(matrixName); @@ -550,4 +587,16 @@ psFitsClose(matrixFile); + if (stamp->vector2) { + matrixName = NULL; + psStringAppend(&matrixName, "vector2_%d.fits", i); + dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64); + memcpy(dummy->data.F64[0], stamp->vector2->data.F64, + PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL); + psFree(dummy); + psFitsClose(matrixFile); + } if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { @@ -575,12 +624,8 @@ } -bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps, - float tolerance) +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); - if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, false); - } // Check inputs @@ -594,31 +639,27 @@ } + PS_ASSERT_VECTOR_NON_NULL(stamp->vector1, false); + if (numParams == -1) { + numParams = stamp->vector1->n; + } + PS_ASSERT_VECTOR_SIZE(stamp->vector1, (long)numParams, false); + PS_ASSERT_VECTOR_TYPE(stamp->vector1, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); + PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false); PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false); - if (numParams == -1) { - numParams = stamp->matrix1->numCols; + if (numParams2 == 0) { numParams2 = stamp->matrix2->numCols; } - PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false); PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, false); - PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false); PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false); - PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); - PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); - PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); - } else { - PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); - if (numParams == -1) { - numParams = stamp->vector->n; - } - PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); - PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); - PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); - PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); - PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); + PS_ASSERT_VECTOR_NON_NULL(stamp->vector2, false); + PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false); + PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false); } } @@ -638,5 +679,5 @@ if (stamp->status == PM_SUBTRACTION_STAMP_USED) { (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); - (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); } } @@ -667,9 +708,11 @@ psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64); psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64); - psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64); psImageInit(sumMatrix1, 0.0); psImageInit(sumMatrix2, 0.0); psImageInit(sumMatrixX, 0.0); - psVectorInit(sumVector, 0.0); + psVectorInit(sumVector1, 0.0); + psVectorInit(sumVector2, 0.0); for (int i = 0; i < stamps->num; i++) { @@ -679,113 +722,171 @@ (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2); (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX); - (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - } - } - - psVector *permutation1 = NULL, *permutation2 = NULL; // Permutation vectors for LUD - psImage *lu1 = psMatrixLUD(NULL, &permutation1, sumMatrix1); // LUD of matrix 1 - psImage *lu2 = psMatrixLUD(NULL, &permutation2, sumMatrix2); // LUD of matrix 2 - psFree(sumMatrix1); - psFree(sumMatrix2); - if (!lu1 || !lu2) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate LU decomposed matrices."); - psFree(permutation1); - psFree(permutation2); - psFree(lu1); - psFree(lu2); - psFree(sumMatrixX); - psFree(sumVector); - return NULL; - } - - psVector *lastSol2 = psVectorAlloc(numParams2, PS_TYPE_F64); // Value of sol2 in last iteration - psVectorInit(lastSol2, 0.0); - - double diff = NAN; // Difference between iterations - psVector *vector = NULL; // RHS for equations - int iter = 0; // Iteration number - do { - if (!vector) { - // Start with traditional Alard-Lupton solution - vector = psVectorCopy(NULL, sumVector, PS_TYPE_F64); - } else { - // Generate vector for RHS of equation 1 - vector->n = numParams; - psVector *sol2 = kernels->solution2; - assert(sol2); - for (int i = 0; i < numParams; i++) { - double value = 0.0; // Value of vector element - for (int j = 0; j < numParams2; j++) { - // XXX Not sure I've got this the right way around - value += sol2->data.F64[j] * sumMatrixX->data.F64[j][i]; - } - vector->data.F64[i] = value; - } - } - - kernels->solution1 = psMatrixLUSolve(kernels->solution1, lu1, vector, permutation1); - if (!kernels->solution1) { - psError(PS_ERR_UNKNOWN, false, "Unable to solve initial matrix equation."); - psFree(permutation1); - psFree(permutation2); - psFree(lu1); - psFree(lu2); - psFree(sumMatrixX); - psFree(sumVector); - psFree(vector); - psFree(lastSol2); - return NULL; - } - - // Generate vector for RHS of equation 2 - vector->n = numParams2; - psVector *sol1 = kernels->solution1; // Solution of first equation - for (int i = 0; i < numParams2; i++) { - double value = 0.0; // Value of vector element - for (int j = 0; j < numParams; j++) { - // XXX Not sure I've got this the right way around - value += sol1->data.F64[j] * sumMatrixX->data.F64[i][j]; - } - vector->data.F64[i] = value; - } - - kernels->solution2 = psMatrixLUSolve(kernels->solution2, lu2, vector, permutation2); - if (!kernels->solution2) { - psError(PS_ERR_UNKNOWN, false, "Unable to solve matrix equation 2."); - psFree(permutation1); - psFree(permutation2); - psFree(lu1); - psFree(lu2); - psFree(sumMatrixX); - psFree(sumVector); - psFree(vector); - psFree(lastSol2); - return NULL; - } - - // Get difference produced by iteration - diff = 0.0; + (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1); + (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); + } + } + #if 0 - psVector *sol2 = kernels->solution2; // Solution of second equation - for (int i = 0; i < numParams2; i++) { - diff += fabs(sol2->data.F64[i] - lastSol2->data.F64[i]); - } - lastSol2 = psVectorCopy(lastSol2, sol2, PS_TYPE_F64); - psTrace("psModules.imcombine", 4, "Solution iteration %d: %lf\n", iter, diff); -#endif - - psVectorInit(kernels->solution2, 0.0); - - iter++; - } while (diff > tolerance); - - psFree(permutation1); - psFree(permutation2); - psFree(lu1); - psFree(lu2); - psFree(sumMatrixX); - psFree(sumVector); - psFree(vector); - psFree(lastSol2); + // Apply weighting to maximise the difference between solution coefficients for the same functions + double fudge = PS_SQR(2 * stamps->footprint + 1); // Fudge factor + for (int i = 0; i < kernels->num; i++) { + sumMatrix1->data.F64[i][i] -= fudge; + sumMatrix2->data.F64[i][i] += fudge; + } +#endif + + // 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 = NAN; + 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 + + { + psFits *fits = psFitsOpen("sumMatrix1.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrix1, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("sumMatrix2.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrix2, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("sumMatrixX.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrixX, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("sumFinverse.fits", "w"); + psFitsWriteImage(fits, NULL, Fi, 0, NULL); + psFitsClose(fits); + } + + + kernels->solution1 = a; + kernels->solution2 = b; + + // XXXXX Free temporary matrices and vectors + } @@ -817,7 +918,4 @@ int numKernels = kernels->num; // Number of kernels - psVector *coeff1 = psVectorAlloc(numKernels, PS_TYPE_F64); // Coefficients for first image - psVector *coeff2 = kernels->mode == PM_SUBTRACTION_MODE_DUAL ? - psVectorAlloc(numKernels, PS_TYPE_F64) : NULL; // Coefficients for second image psImage *polyValues = NULL; // Polynomial values psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image @@ -832,7 +930,4 @@ // Calculate coefficients of the kernel basis functions polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); - for (int j = 0; j < numKernels; j++) { - coeff1->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); - } double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background @@ -863,5 +958,6 @@ for (int j = 0; j < numKernels; j++) { psKernel *convolution = convolutions->data[j]; // Convolution - double coefficient = coeff1->data.F64[j]; // Coefficient + double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, + false); // Coefficient for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -880,20 +976,15 @@ psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image psKernel *image1 = stamp->image1; // The first image - - for (int j = 0; j < numKernels; j++) { - coeff2->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); - } + psKernel *image2 = stamp->image2; // The second image for (int j = 0; j < numKernels; j++) { psKernel *conv1 = convolutions1->data[j]; // Convolution of first image psKernel *conv2 = convolutions2->data[j]; // Convolution of second image - - double coefficient1 = coeff1->data.F64[j]; // Coefficient for convolution 1 - double coefficient2 = coeff2->data.F64[j]; // Coefficient for convolution 2 + double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1 + double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2 for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] += conv1->kernel[y][x] * coefficient1 + - conv2->kernel[y][x] * coefficient2; + residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1; } } @@ -901,5 +992,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] += image1->kernel[y][x] * norm - background; + residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm; } } @@ -909,5 +1000,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - deviation += PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x]; + double dev = PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x]; + deviation += dev; +#ifdef TESTING + residual->kernel[y][x] = dev; +#endif } } @@ -932,4 +1027,20 @@ psFitsClose(fits); } + { + psString filename = NULL; + psStringAppend(&filename, "stamp_image_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL); + psFitsClose(fits); + } + { + psString filename = NULL; + psStringAppend(&filename, "stamp_weight_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); + psFitsClose(fits); + } #endif @@ -937,5 +1048,4 @@ psFree(residual); psFree(polyValues); - psFree(coeff1); return deviations; Index: trunk/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 15799) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 15809) @@ -12,6 +12,5 @@ /// Solve the least-squares equation to match the image quality bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters - const pmSubtractionStampList *stamps, ///< Stamps - float tolerance ///< Tolerance for iteration (dual mode only) + const pmSubtractionStampList *stamps ///< Stamps ); Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 15799) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 15809) @@ -136,4 +136,5 @@ } } + preCalc->kernel[0][0] -= 1.0; } @@ -216,19 +217,4 @@ } - // Subtract normalisation from all the others to maintain flux scaling - if (spatialOrder > 0) { - int numGaussians = fwhms->n; // Number of Gaussians - for (int i = 0, index = 0; i < numGaussians; i++) { - for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { - for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { - if (uOrder % 2 == 0 && vOrder % 2 == 0) { - psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract - kernel->kernel[0][0] -= 1.0; - } - } - } - } - } - if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { for (int i = 0; i < kernels->num; i++) { @@ -246,5 +232,5 @@ } } - psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %lf\n", i, sum); + psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %le\n", i, sum); } } @@ -447,19 +433,4 @@ psStringPrepend(&kernels->description, "GUNK="); psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); - - // Subtract unity from the kernels to maintain photometric flux scaling - int numGaussians = fwhms->n; // Number of Gaussians - if (spatialOrder > 0) { - for (int i = 0, index = 0; i < numGaussians; i++) { - for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { - for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { - if (uOrder % 2 == 0 && vOrder % 2 == 0) { - psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract - kernel->kernel[0][0] -= 1.0; - } - } - } - } - } int numISIS = kernels->num; // Number of ISIS kernels Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15799) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15809) @@ -18,5 +18,4 @@ #include "pmSubtractionMatch.h" -#define TOL 1.0e-3 // Tolerance for subtraction solution #define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image @@ -24,10 +23,11 @@ -//#define TESTING +#define TESTING +//#define TESTING_MEMORY // Output memory usage information static void memCheck(const char *where) { -#ifdef TESTING +#ifdef TESTING_MEMORY psMemBlock **leaks = NULL; int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true); @@ -216,5 +216,4 @@ psString regionString = NULL; // String for region pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF - psVector *solution = NULL; // Solution to match PSF pmSubtractionKernels *kernels = NULL; // Kernel basis functions int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions @@ -261,5 +260,5 @@ stampSpacing, mode); } else if (stampsName && strlen(stampsName) > 0) { - stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, footprint, + stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, footprint, stampSpacing, mode); } @@ -331,5 +330,5 @@ psTrace("psModules.imcombine", 3, "Solving equation...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) { + if (!pmSubtractionSolveEquation(kernels, stamps)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); goto MATCH_ERROR; @@ -360,5 +359,5 @@ if (numRejected > 0) { psTrace("psModules.imcombine", 3, "Solving equation...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) { + if (!pmSubtractionSolveEquation(kernels, stamps)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); goto MATCH_ERROR; @@ -413,5 +412,5 @@ } -#ifdef TESTING +#if 0 { // Generate images of the kernel components @@ -442,4 +441,5 @@ kernels = NULL; +#if 0 // Put the solution on the metadata { @@ -460,7 +460,5 @@ } } - - psFree(solution); - solution = NULL; +#endif // There is data in the readout now @@ -496,4 +494,18 @@ +#ifdef TESTING + { + psFits *fits = psFitsOpen("convolved1.fits", "w"); + psFitsWriteImage(fits, NULL, conv1->image, 0, NULL); + psFitsClose(fits); + + if (mode == PM_SUBTRACTION_MODE_DUAL) { + psFits *fits = psFitsOpen("convolved2.fits", "w"); + psFitsWriteImage(fits, NULL, conv2->image, 0, NULL); + psFitsClose(fits); + } + } +#endif + return true; @@ -504,5 +516,4 @@ psFree(kernels); psFree(stamps); - psFree(solution); psFree(weight); return false; Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 15799) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 15809) @@ -63,7 +63,7 @@ psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL - psImage *matrix1, *matrix2; ///< Matrices for each image, or NULL + psImage *matrix1, *matrix2; ///< Least-squares matrices for each image, or NULL psImage *matrixX; ///< Cross-matrix (for mode DUAL), or NULL - psVector *vector; ///< Associated vector (when mode not DUAL), or NULL + psVector *vector1, *vector2; ///< Least-squares vectors for each image, or NULL pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp; @@ -108,4 +108,5 @@ pmSubtractionStampList *pmSubtractionStampsSetFromFile( const char *filename, ///< Filename of file containing x,y (or x,y,flux) on each line + const psImage *image, ///< Image for flux of stamp const psImage *subMask, ///< Mask, or NULL const psRegion *region, ///< Region to search, or NULL