Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26553) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26562) @@ -907,11 +907,8 @@ 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; + psFree(stamp->matrix); + stamp->matrix = NULL; + psFree(stamp->vector); + stamp->vector = NULL; } else { numGood++; Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26553) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26562) @@ -15,5 +15,5 @@ #include "pmSubtractionVisual.h" -#define TESTING // TESTING output for debugging; may not work with threads! +// #define TESTING // TESTING output for debugging; may not work with threads! // #define USE_WEIGHT // Include weight (1/variance) in equation? @@ -26,5 +26,4 @@ // Calculate the least-squares matrix and vector -// XXX why are we calculating these values on each iteration? aren't they identical (no change in pixel values...) static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated psVector *vector, // Least-squares vector, updated @@ -221,10 +220,8 @@ } + // 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 +static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated const psKernel *image1, // Image 1 const psKernel *image2, // Image 2 @@ -238,16 +235,4 @@ ) { - // 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 @@ -257,4 +242,19 @@ double poly[numPoly]; // Polynomial terms double poly2[numPoly][numPoly]; // Polynomial-polynomial values + + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters + int numParams2 = numKernels * numPoly; // Number of additional parameters for dual + int numDual = numParams + numParams2; // Total number of parameters for dual + + psAssert(matrix && + matrix->type.type == PS_TYPE_F64 && + matrix->numCols == numDual && + matrix->numRows == numDual, + "Least-squares matrix is bad."); + psAssert(vector && + vector->type.type == PS_TYPE_F64 && + vector->n == numDual, + "Least-squares vector is bad."); // Evaluate polynomial-polynomial terms @@ -280,7 +280,7 @@ 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 + double sumAA = 0.0; // Sum of convolution products between image 1 + double sumBB = 0.0; // Sum of convolution products between image 2 + double sumAB = 0.0; // Sum of convolution products across images 1 and 2 for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -312,9 +312,13 @@ 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; + + matrix->data.F64[iIndex][jIndex] = aa; + matrix->data.F64[jIndex][iIndex] = aa; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; + matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; + + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; } } @@ -340,16 +344,17 @@ 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) + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } + + double sumAI2 = 0.0; // Sum of A.I_2 products (for vector) + double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) + double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) + double sumA = 0.0; // Sum of A (for matrix, background) + double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) + double sumB = 0.0; // Sum of B products (for matrix, background) + double sumI2 = 0.0; // Sum of I_2 (for vector, background) for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -402,20 +407,22 @@ 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) + matrix->data.F64[iIndex][normIndex] = ai1; + matrix->data.F64[normIndex][iIndex] = ai1; + matrix->data.F64[iIndex][bgIndex] = a; + matrix->data.F64[bgIndex][iIndex] = a; + matrix->data.F64[iIndex + numParams][normIndex] = bi1; + matrix->data.F64[normIndex][iIndex + numParams] = bi1; + matrix->data.F64[iIndex + numParams][bgIndex] = b; + matrix->data.F64[bgIndex][iIndex + numParams] = b; + vector->data.F64[iIndex] = ai2; + vector->data.F64[iIndex + numParams] = bi2; + } + } + + double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) + double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) + double sum1 = 0.0; // Sum of 1 (for matrix, background-background) + double sumI2 = 0.0; // Sum of I_2 (for vector, background) + double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -450,77 +457,13 @@ } } - 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; + matrix->data.F64[bgIndex][normIndex] = sumI1; + matrix->data.F64[normIndex][bgIndex] = sumI1; + matrix->data.F64[normIndex][normIndex] = sumI1I1; + matrix->data.F64[bgIndex][bgIndex] = sum1; + vector->data.F64[bgIndex] = sumI2; + vector->data.F64[normIndex] = sumI1I2; return true; } - -// 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 - ) -{ - 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]; - } - } - - return true; -} - #if 1 @@ -562,9 +505,9 @@ // Calculate the value of a polynomial, specified by coefficients and polynomial values double p_pmSubtractionCalculatePolynomial(const psVector *coeff, // Coefficients - const psImage *polyValues, // Polynomial values - int order, // Order of polynomials - int index, // Index at which to begin - int step // Step between subsequent indices - ) + const psImage *polyValues, // Polynomial values + int order, // Order of polynomials + int index, // Index at which to begin + int step // Step between subsequent indices + ) { double sum = 0.0; // Value of the polynomial sum @@ -581,5 +524,5 @@ double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, const psImage *polyValues, - int index, bool wantDual) + int index, bool wantDual) { #if 0 @@ -659,4 +602,8 @@ // number of coefficients for the spatial polynomial, normalisation and a constant background offset. int numParams = numKernels * numSpatial + 1 + numBackground; + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + // An additional image is convolved + numParams += numKernels * numSpatial; + } pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest @@ -698,12 +645,12 @@ stamp->xNorm, stamp->yNorm); // Polynomial terms - bool new = stamp->vector1 ? false : true; // Is this a new run? + bool new = stamp->vector ? false : true; // Is this a new run? if (new) { - stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); - stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64); + stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64); } #ifdef TESTING - psImageInit(stamp->matrix1, NAN); - psVectorInit(stamp->vector1, NAN); + psImageInit(stamp->matrix, NAN); + psVectorInit(stamp->vector, NAN); #endif @@ -722,28 +669,17 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: - status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1, weight, window, stamp->convolutions1, kernels, polyValues, footprint, mode); break; case PM_SUBTRACTION_MODE_2: - status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2, + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2, weight, window, stamp->convolutions2, kernels, polyValues, footprint, mode); break; case PM_SUBTRACTION_MODE_DUAL: - if (new) { - 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->matrix2, NAN); - psImageInit(stamp->matrixX, NAN); - psVectorInit(stamp->vector2, NAN); -#endif - status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, - stamp->matrixX, stamp->image1, stamp->image2, weight, window, - stamp->convolutions1, stamp->convolutions2, kernels, polyValues, - footprint); + status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2, + weight, window, stamp->convolutions1, stamp->convolutions2, + kernels, polyValues, footprint); break; default: @@ -762,15 +698,15 @@ { psString matrixName = NULL; - psStringAppend(&matrixName, "matrix1_%d.fits", index); + psStringAppend(&matrixName, "matrix_%d.fits", index); psFits *matrixFile = psFitsOpen(matrixName, "w"); psFree(matrixName); - psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL); + psFitsWriteImage(matrixFile, NULL, stamp->matrix, 0, NULL); psFitsClose(matrixFile); matrixName = NULL; - psStringAppend(&matrixName, "vector1_%d.fits", index); - 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); + psStringAppend(&matrixName, "vector_%d.fits", index); + 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); matrixFile = psFitsOpen(matrixName, "w"); psFree(matrixName); @@ -778,33 +714,4 @@ psFree(dummy); psFitsClose(matrixFile); - - if (stamp->vector2) { - matrixName = NULL; - psStringAppend(&matrixName, "vector2_%d.fits", index); - 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) { - matrixName = NULL; - psStringAppend(&matrixName, "matrix2_%d.fits", index); - matrixFile = psFitsOpen(matrixName, "w"); - psFree(matrixName); - psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL); - psFitsClose(matrixFile); - - matrixName = NULL; - psStringAppend(&matrixName, "matrixX_%d.fits", index); - matrixFile = psFitsOpen(matrixName, "w"); - psFree(matrixName); - psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL); - psFitsClose(matrixFile); - } } #endif @@ -898,6 +805,15 @@ // Check inputs - int numParams = -1; // Number of parameters - int numParams2 = 0; // Number of parameters for part solution (DUAL mode) + int numKernels = kernels->num; // Number of kernel basis functions + int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for + int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + // An additional image is convolved + numSolution2 = numKernels * numSpatial; + numParams += numSolution2; + } + for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -907,32 +823,10 @@ } - 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->matrix2, false); - PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false); - if (numParams2 == 0) { - numParams2 = stamp->matrix2->numCols; - } - PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false); - PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, 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->vector2, false); - PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false); - PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false); - } - } - if (numParams == -1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found."); - return NULL; + 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); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix, false); + PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, false); } @@ -958,20 +852,20 @@ #ifdef TESTING // XXX double-check for NAN in data: - for (int iy = 0; iy < stamp->matrix1->numRows; iy++) { - for (int ix = 0; ix < stamp->matrix1->numCols; ix++) { - if (!isfinite(stamp->matrix1->data.F64[iy][ix])) { - fprintf (stderr, "WARNING: NAN in matrix1\n"); + for (int iy = 0; iy < stamp->matrix->numRows; iy++) { + for (int ix = 0; ix < stamp->matrix->numCols; ix++) { + if (!isfinite(stamp->matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); } } } - for (int ix = 0; ix < stamp->vector1->n; ix++) { - if (!isfinite(stamp->vector1->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector1\n"); + for (int ix = 0; ix < stamp->vector->n; ix++) { + if (!isfinite(stamp->vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); } } #endif - (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); - (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; @@ -984,5 +878,5 @@ for (int ix = 0; ix < sumVector->n; ix++) { if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector1\n"); + fprintf (stderr, "WARNING: NAN in vector\n"); } } @@ -997,5 +891,5 @@ for (int ix = 0; ix < sumVector->n; ix++) { if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector1\n"); + fprintf (stderr, "WARNING: NAN in vector\n"); } } @@ -1016,5 +910,5 @@ #endif -# define SVD_ANALYSIS 1 +# define SVD_ANALYSIS 0 # define COEFF_SIG 0.0 # define SVD_TOL 0.0 @@ -1282,5 +1176,5 @@ for (int ix = 0; ix < kernels->solution1->n; ix++) { if (!isfinite(kernels->solution1->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector1\n"); + fprintf (stderr, "WARNING: NAN in vector\n"); } } @@ -1291,14 +1185,8 @@ // Accumulation of stamp matrices/vectors - psImage *sumMatrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); - psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64); - psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64); - psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64); - psVector *sumVector2 = psVectorAlloc(numParams2, PS_TYPE_F64); - psImageInit(sumMatrix1, 0.0); - psImageInit(sumMatrix2, 0.0); - psImageInit(sumMatrixX, 0.0); - psVectorInit(sumVector1, 0.0); - psVectorInit(sumVector2, 0.0); + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix, 0.0); + psVectorInit(sumVector, 0.0); int numStamps = 0; // Number of good stamps @@ -1306,9 +1194,7 @@ pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (stamp->status == PM_SUBTRACTION_STAMP_USED) { - (void)psBinaryOp(sumMatrix1, sumMatrix1, "+", stamp->matrix1); - (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2); - (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX); - (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1); - (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; @@ -1318,45 +1204,20 @@ #ifdef TESTING - psFitsWriteImageSimple ("sumMatrix1.fits", sumMatrix1, NULL); - psFitsWriteImageSimple ("sumMatrix2.fits", sumMatrix2, NULL); - psFitsWriteImageSimple ("sumMatrixX.fits", sumMatrixX, NULL); - psVectorWriteFile("sumVector1.dat", sumVector1); - psVectorWriteFile("sumVector2.dat", sumVector2); -#endif - + psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); + psVectorWriteFile("sumVector.dat", sumVector); +#endif #if 1 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]); -#endif - - psImage *sumMatrix = NULL; // Combined matrix - psVector *sumVector = NULL; // Combined vector - calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2, - sumMatrixX, sumVector1, sumVector2); - -#ifdef TESTING - psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, 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); - } - psVectorWriteFile("sumVector.dat", sumVector); + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); #endif psVector *solution = NULL; // Solution to equation! - solution = psVectorAlloc(numParams + numParams2, PS_TYPE_F64); + solution = psVectorAlloc(numParams, PS_TYPE_F64); psVectorInit(solution, 0); #if 0 { - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-5); + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); if (!solution) { psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n"); @@ -1367,14 +1228,6 @@ #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++) { - fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]); - image->data.F64[0][i] = solution->data.F64[i]; - } - psFitsWriteImage(fits, NULL, image, 0, NULL); - psFree(image); - psFitsClose(fits); + for (int i = 0; i < numParams; i++) { + fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]); } #endif @@ -1382,4 +1235,15 @@ int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations int numKernels = kernels->num; // Number of kernel basis functions + +#define MASK_BASIS(INDEX) \ + { \ + for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \ + for (int k = 0; k < numParams; k++) { \ + sumMatrix->data.F64[index][k] = sumMatrix->data.F64[k][index] = 0.0; \ + } \ + sumMatrix->data.F64[index][index] = 1.0; \ + sumVector->data.F64[index] = 0.0; \ + } \ + } // Remove a kernel basis for image 1 from the equation @@ -1427,14 +1291,14 @@ // 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 + param2 = solution->data.F64[numSolution1 + 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); + MASK_BASIS(i); mask1 = true; } if (fabs(param2) < thresh) { psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i); - MASK_BASIS_2(i); + MASK_BASIS(numSolution1 + i); mask2 = true; } @@ -1443,77 +1307,41 @@ if (fabs(param1) > fabs(param2)) { psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i); - MASK_BASIS_2(i); + MASK_BASIS(numSolution1 + i); } else { psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i); - MASK_BASIS_1(i); + MASK_BASIS(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); - } - { - 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 - } -#endif - - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-7); - 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); + psFitsWriteImageSimple ("sumMatrixFix.fits", sumMatrix, NULL); + psVectorWriteFile("sumVectorFix.dat", sumVector); +#endif + + } +#endif + + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + + + for (int i = 0; i < solution->n; i++) { + fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); + } psFree(sumMatrix); psFree(sumVector); - -#ifdef TESTING - { - psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); - psFits *fits = psFitsOpen("solnVectorFix.fits", "w"); - for (int i = 0; i < numParams + numParams2; i++) { - fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[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); + kernels->solution1 = psVectorAlloc(numSolution1, 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)); + kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64); + } + + memcpy(kernels->solution1->data.F64, solution->data.F64, + numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], + numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); psFree(solution); Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26553) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26562) @@ -58,11 +58,6 @@ psFree(stamp->convolutions1); psFree(stamp->convolutions2); - - psFree(stamp->matrix1); - psFree(stamp->matrix2); - psFree(stamp->matrixX); - psFree(stamp->vector1); - psFree(stamp->vector2); - + psFree(stamp->matrix); + psFree(stamp->vector); } @@ -284,9 +279,6 @@ } - outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL; - outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL; - outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL; - outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL; - outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL; + outStamp->matrix = inStamp->matrix ? psImageCopy(NULL, inStamp->matrix, PS_TYPE_F64) : NULL; + outStamp->vector = inStamp->vector ? psVectorCopy(NULL, inStamp->vector, PS_TYPE_F64) : NULL; out->stamps->data[i] = outStamp; @@ -314,9 +306,6 @@ stamp->convolutions2 = NULL; - stamp->matrix1 = NULL; - stamp->matrix2 = NULL; - stamp->matrixX = NULL; - stamp->vector1 = NULL; - stamp->vector2 = NULL; + stamp->matrix = NULL; + stamp->vector = NULL; return stamp; @@ -384,5 +373,5 @@ if (!stamps) { - stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr); + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr); } @@ -439,8 +428,8 @@ subMask, xMin, xMax, yMin, yMax, numCols, numRows, border); - // XXX reset for a test: - xStamp = xList->data.F32[j]; - yStamp = yList->data.F32[j]; - // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp); + // XXX reset for a test: + xStamp = xList->data.F32[j]; + yStamp = yList->data.F32[j]; + // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp); } } else { @@ -549,5 +538,5 @@ } - // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix); + // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix); bool found = false; @@ -635,17 +624,17 @@ pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (!stamp) continue; - if (!stamp->image1) continue; - if (!stamp->image2) continue; - - float sum1 = 0.0; - float sum2 = 0.0; - for (int y = -size; y <= size; y++) { - for (int x = -size; x <= size; x++) { - sum1 += stamp->image1->kernel[y][x]; - sum2 += stamp->image2->kernel[y][x]; - } - } - norm1->data.F32[i] = sum1; - norm2->data.F32[i] = sum2; + if (!stamp->image1) continue; + if (!stamp->image2) continue; + + float sum1 = 0.0; + float sum2 = 0.0; + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + sum1 += stamp->image1->kernel[y][x]; + sum2 += stamp->image2->kernel[y][x]; + } + } + norm1->data.F32[i] = sum1; + norm2->data.F32[i] = sum2; } @@ -658,40 +647,40 @@ // generate the window pixels for (int y = -size; y <= size; y++) { - for (int x = -size; x <= size; x++) { - - flux->n = 0; - for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (!stamp) continue; - if (!stamp->image1) continue; - if (!stamp->image2) continue; - - psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]); - psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]); - } - - psStatsInit (stats); - if (!psVectorStats (stats, flux, NULL, NULL, 0)) { - psAbort ("failed to generate stats"); - } - stamps->window->kernel[y][x] = stats->robustMedian; - if (maxValue < stats->robustMedian) { - maxValue = stats->robustMedian; - } - } + for (int x = -size; x <= size; x++) { + + flux->n = 0; + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (!stamp) continue; + if (!stamp->image1) continue; + if (!stamp->image2) continue; + + psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]); + psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]); + } + + psStatsInit (stats); + if (!psVectorStats (stats, flux, NULL, NULL, 0)) { + psAbort ("failed to generate stats"); + } + stamps->window->kernel[y][x] = stats->robustMedian; + if (maxValue < stats->robustMedian) { + maxValue = stats->robustMedian; + } + } } // re-normalize so chisquare values are sensible for (int y = -size; y <= size; y++) { - for (int x = -size; x <= size; x++) { - stamps->window->kernel[y][x] /= maxValue; - } + for (int x = -size; x <= size; x++) { + stamps->window->kernel[y][x] /= maxValue; + } } # ifdef TESTING { - psFits *fits = psFitsOpen ("window.fits", "w"); - psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL); - psFitsClose (fits); + psFits *fits = psFitsOpen ("window.fits", "w"); + psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL); + psFitsClose (fits); } # endif @@ -743,5 +732,5 @@ return false; } - // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y); + // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y); // Catch memory leaks --- these should have been freed and NULLed before @@ -768,7 +757,7 @@ if ((isfinite(stamps->skyErr) && (stamps->skyErr > 0)) || - (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) { + (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) { float sysErr = 0.25 * PS_SQR(stamps->sysErr); // Systematic error - float skyErr = stamps->skyErr; + float skyErr = stamps->skyErr; psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images for (int y = -size; y <= size; y++) { @@ -823,5 +812,5 @@ y->data.F32[numOut] = source->peak->yf; } - // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]); + // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]); numOut++; } @@ -873,8 +862,8 @@ for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (!stamp) continue; - if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; - stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (!stamp) continue; + if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; + stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; } return true; Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h (revision 26553) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h (revision 26562) @@ -24,5 +24,5 @@ psArray *flux; ///< Fluxes for possible stamps (or NULL) int footprint; ///< Half-size of stamps - psKernel *window; ///< window function generated from ensemble of stamps + psKernel *window; ///< window function generated from ensemble of stamps float sysErr; ///< Systematic error float skyErr; ///< Systematic error @@ -76,7 +76,6 @@ 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; ///< Least-squares matrices for each image, or NULL - psImage *matrixX; ///< Cross-matrix (for mode DUAL), or NULL - psVector *vector1, *vector2; ///< Least-squares vectors for each image, or NULL + psImage *matrix; ///< Least-squares matrix, or NULL + psVector *vector; ///< Least-squares vector, or NULL pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp; Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26553) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26562) @@ -159,7 +159,5 @@ if (stamp == NULL) continue; - psImage *im = stamp->matrix1; - if (im == NULL) im = stamp->matrix2; - if (im == NULL) im = stamp->matrixX; + psImage *im = stamp->matrix; if (im == NULL) continue; @@ -189,7 +187,5 @@ if (stamp == NULL) continue; - psImage *im = stamp->matrix1; - if (im == NULL) im = stamp->matrix2; - if (im == NULL) im = stamp->matrixX; + psImage *im = stamp->matrix; if (im == NULL) continue;