Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26546) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26547) @@ -117,4 +117,9 @@ for (int i = 0; i < numKernels; i++) { double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value + if (wantDual) { + // The model is built with the dual convolution terms added, so to produce zero residual the + // equation results in negative coefficients which we must undo. + value *= -1.0; + } switch (kernels->type) { Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 26546) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 26547) @@ -16,5 +16,5 @@ #define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image -//#define TESTING +#define TESTING bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header, Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26546) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26547) @@ -15,7 +15,7 @@ #include "pmSubtractionVisual.h" -// #define TESTING // TESTING output for debugging; may not work with threads! - -#define USE_WEIGHT // Include weight (1/variance) in equation? +#define TESTING // TESTING output for debugging; may not work with threads! + +// #define USE_WEIGHT // Include weight (1/variance) in equation? // #define USE_WINDOW // Include weight (1/variance) in equation? @@ -37,5 +37,5 @@ const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp - const pmSubtractionEquationCalculationMode mode + const pmSubtractionEquationCalculationMode mode ) { @@ -69,10 +69,10 @@ } - // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we + // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we // choose to calculate psImageInit(matrix, 0.0); psVectorInit(vector, 1.0); for (int i = 0; i < matrix->numCols; i++) { - matrix->data.F64[i][i] = 1.0; + matrix->data.F64[i][i] = 1.0; } @@ -85,87 +85,87 @@ 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]; - if (weight) { - cc *= weight->kernel[y][x]; - } - if (window) { - cc *= window->kernel[y][x]; - } - sumCC += cc; - } - } - - // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - 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; - if (weight) { - float wtVal = weight->kernel[y][x]; - ic *= wtVal; - rc *= wtVal; - c *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - ic *= winVal; - rc *= winVal; - c *= winVal; - } - 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]; - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = normTerm; - matrix->data.F64[normIndex][iIndex] = normTerm; - } - if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = bgTerm; - matrix->data.F64[bgIndex][iIndex] = bgTerm; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = sumIC * poly[iTerm]; - // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B), - // instead, calculate A - \sum(k x B), with full hermitians - if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) { - // subtract norm * sumRC * poly[iTerm] - psAssert (kernels->solution1, "programming error: define solution first!"); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation - vector->data.F64[iIndex] -= norm * normTerm; - } - } - } + 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]; + if (weight) { + cc *= weight->kernel[y][x]; + } + if (window) { + cc *= window->kernel[y][x]; + } + sumCC += cc; + } + } + + // Spatial variation of kernel coeffs + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + 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; + if (weight) { + float wtVal = weight->kernel[y][x]; + ic *= wtVal; + rc *= wtVal; + c *= wtVal; + } + if (window) { + float winVal = window->kernel[y][x]; + ic *= winVal; + rc *= winVal; + c *= winVal; + } + 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]; + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][normIndex] = normTerm; + matrix->data.F64[normIndex][iIndex] = normTerm; + } + if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][bgIndex] = bgTerm; + matrix->data.F64[bgIndex][iIndex] = bgTerm; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + vector->data.F64[iIndex] = sumIC * poly[iTerm]; + // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B), + // instead, calculate A - \sum(k x B), with full hermitians + if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) { + // subtract norm * sumRC * poly[iTerm] + psAssert (kernels->solution1, "programming error: define solution first!"); + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation + vector->data.F64[iIndex] -= norm * normTerm; + } + } + } } @@ -182,20 +182,20 @@ double rr = PS_SQR(ref); double one = 1.0; - if (weight) { - float wtVal = weight->kernel[y][x]; - rr *= wtVal; - ir *= wtVal; - in *= wtVal; - ref *= wtVal; - one *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - rr *= winVal; - ir *= winVal; - in *= winVal; - ref *= winVal; - one *= winVal; - } + if (weight) { + float wtVal = weight->kernel[y][x]; + rr *= wtVal; + ir *= wtVal; + in *= wtVal; + ref *= wtVal; + one *= wtVal; + } + if (window) { + float winVal = window->kernel[y][x]; + rr *= winVal; + ir *= winVal; + in *= winVal; + ref *= winVal; + one *= winVal; + } sumRR += rr; sumIR += ir; @@ -206,15 +206,15 @@ } if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumRR; - vector->data.F64[normIndex] = sumIR; - // subtract sum over kernels * kernel solution + matrix->data.F64[normIndex][normIndex] = sumRR; + vector->data.F64[normIndex] = sumIR; + // subtract sum over kernels * kernel solution } if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI; + matrix->data.F64[bgIndex][bgIndex] = sum1; + vector->data.F64[bgIndex] = sumI; } if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[normIndex][bgIndex] = sumR; - matrix->data.F64[bgIndex][normIndex] = sumR; + matrix->data.F64[normIndex][bgIndex] = sumR; + matrix->data.F64[bgIndex][normIndex] = sumR; } return true; @@ -288,16 +288,16 @@ double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; - if (weight) { - float wtVal = weight->kernel[y][x]; - aa *= wtVal; - bb *= wtVal; - ab *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - aa *= wtVal; - bb *= wtVal; - ab *= wtVal; - } + if (weight) { + float wtVal = weight->kernel[y][x]; + aa *= wtVal; + bb *= wtVal; + ab *= wtVal; + } + if (window) { + float wtVal = window->kernel[y][x]; + aa *= wtVal; + bb *= wtVal; + ab *= wtVal; + } sumAA += aa; sumBB += bb; @@ -326,10 +326,10 @@ for (int x = - footprint; x <= footprint; x++) { double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; - if (weight) { - ab *= weight->kernel[y][x]; - } - if (window) { - ab *= window->kernel[y][x]; - } + if (weight) { + ab *= weight->kernel[y][x]; + } + if (window) { + ab *= window->kernel[y][x]; + } sumAB += ab; } @@ -352,9 +352,8 @@ 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]; + double a = iConv1->kernel[y][x]; + double b = iConv2->kernel[y][x]; float i1 = image1->kernel[y][x]; float i2 = image2->kernel[y][x]; @@ -364,28 +363,25 @@ double ai1 = a * i1; double bi1 = b * i1; - double i1i2 = i1 * i2; - - if (weight) { - float wtVal = weight->kernel[y][x]; - ai2 *= wtVal; - bi2 *= wtVal; - ai1 *= wtVal; - bi1 *= wtVal; - i1i2 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - ai2 *= wtVal; - bi2 *= wtVal; - ai1 *= wtVal; - bi1 *= wtVal; - i1i2 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; - } + + if (weight) { + float wtVal = weight->kernel[y][x]; + ai2 *= wtVal; + bi2 *= wtVal; + ai1 *= wtVal; + bi1 *= wtVal; + a *= wtVal; + b *= wtVal; + i2 *= wtVal; + } + if (window) { + float wtVal = window->kernel[y][x]; + ai2 *= wtVal; + bi2 *= wtVal; + ai1 *= wtVal; + bi1 *= wtVal; + a *= wtVal; + b *= wtVal; + i2 *= wtVal; + } sumAI2 += ai2; sumBI2 += bi2; @@ -395,5 +391,4 @@ sumB += b; sumI2 += i2; - sumI1I2 += i1i2; } } @@ -425,6 +420,6 @@ 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 i1 = image1->kernel[y][x]; + double i2 = image2->kernel[y][x]; double i1i1 = i1 * i1; @@ -432,20 +427,20 @@ double i1i2 = i1 * i2; - if (weight) { - float wtVal = weight->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } + if (weight) { + float wtVal = weight->kernel[y][x]; + i1 *= wtVal; + i1i1 *= wtVal; + one *= wtVal; + i2 *= wtVal; + i1i2 *= wtVal; + } + if (window) { + float wtVal = window->kernel[y][x]; + i1 *= wtVal; + i1i1 *= wtVal; + one *= wtVal; + i2 *= wtVal; + i1i2 *= wtVal; + } sumI1 += i1; sumI1I1 += i1i1; @@ -515,5 +510,5 @@ 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]; + matrix->data.F64[i][k] = sumMatrixX->data.F64[j][i]; } } @@ -521,5 +516,5 @@ 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]; + matrix->data.F64[i2][k] = sumMatrix2->data.F64[i1][j]; } } @@ -529,4 +524,5 @@ +#if 1 // Add in penalty term to least-squares vector static bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term @@ -547,7 +543,7 @@ for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { // Contribution to chi^2: a_i^2 P_i - if (!isfinite(penalties->data.F32[i])) { - psAbort ("invalid penalty"); - } + if (!isfinite(penalties->data.F32[i])) { + psAbort ("invalid penalty"); + } matrix->data.F64[index][index] -= norm * penalties->data.F32[i]; } @@ -557,4 +553,6 @@ return true; } +#endif + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -657,5 +655,5 @@ int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). + // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3 @@ -728,10 +726,10 @@ status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, weight, window, stamp->convolutions1, kernels, - polyValues, footprint, mode); + polyValues, footprint, mode); break; case PM_SUBTRACTION_MODE_2: status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2, weight, window, stamp->convolutions2, kernels, - polyValues, footprint, mode); + polyValues, footprint, mode); break; case PM_SUBTRACTION_MODE_DUAL: @@ -834,10 +832,10 @@ } - if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) { - psAbort ("bad stamp"); - } - if (!isfinite(stamp->x) && !isfinite(stamp->y)) { - psAbort ("bad stamp"); - } + if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) { + psAbort ("bad stamp"); + } + if (!isfinite(stamp->x) && !isfinite(stamp->y)) { + psAbort ("bad stamp"); + } if (pmSubtractionThreaded()) { @@ -894,7 +892,7 @@ double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); -bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, - const pmSubtractionStampList *stamps, - const pmSubtractionEquationCalculationMode mode) +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, + const pmSubtractionStampList *stamps, + const pmSubtractionEquationCalculationMode mode) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); @@ -1024,261 +1022,261 @@ # define SVD_TOL 0.0 - psVector *solution = NULL; - psVector *solutionErr = NULL; + psVector *solution = NULL; + psVector *solutionErr = NULL; #ifdef TESTING - psFitsWriteImageSimple("A.fits", sumMatrix, NULL); - psVectorWriteFile ("B.dat", sumVector); -#endif - - // Use SVD to determine the kernel coeffs (and validate) - if (SVD_ANALYSIS) { - - // We have sumVector and sumMatrix. we are trying to solve the following equation: - // sumMatrix * x = sumVector. - - // we can use any standard matrix inversion to solve this. However, the basis - // functions in general have substantial correlation, so that the solution may be - // somewhat poorly determined or unstable. If not numerically ill-conditioned, the - // system of equations may be statistically ill-conditioned. Noise in the image - // will drive insignificant, but correlated, terms in the solution. To avoid these - // problems, we can use SVD to identify numerically unconstrained values and to - // avoid statistically badly determined value. - - // A = sumMatrix, B = sumVector - // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T - // x = V (1/w) (U^T B) - // \sigma_x = sqrt(diag(A^{-1})) - // solve for x and A^{-1} to get x & dx - // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 - // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 - - // If I use the SVD trick to re-condition the matrix, I need to break out the - // kernel and normalization terms from the background term. - // XXX is this true? or was this due to an error in the analysis? - - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - - // now pull out the kernel elements into their own square matrix - psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); - psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { - if (ix == bgIndex) continue; - for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { - if (iy == bgIndex) continue; - kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; - ky++; - } - kernelVector->data.F64[kx] = sumVector->data.F64[ix]; - kx++; - } - - psImage *U = NULL; - psImage *V = NULL; - psVector *w = NULL; - if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { - psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); - return NULL; - } - - // calculate A_inverse: - // Ainv = V * w * U^T - psImage *wUt = p_pmSubSolve_wUt (w, U); - psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); - psImage *Xvar = NULL; - psFree (wUt); + psFitsWriteImageSimple("A.fits", sumMatrix, NULL); + psVectorWriteFile("B.dat", sumVector); +#endif + + // Use SVD to determine the kernel coeffs (and validate) + if (SVD_ANALYSIS) { + + // We have sumVector and sumMatrix. we are trying to solve the following equation: + // sumMatrix * x = sumVector. + + // we can use any standard matrix inversion to solve this. However, the basis + // functions in general have substantial correlation, so that the solution may be + // somewhat poorly determined or unstable. If not numerically ill-conditioned, the + // system of equations may be statistically ill-conditioned. Noise in the image + // will drive insignificant, but correlated, terms in the solution. To avoid these + // problems, we can use SVD to identify numerically unconstrained values and to + // avoid statistically badly determined value. + + // A = sumMatrix, B = sumVector + // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T + // x = V (1/w) (U^T B) + // \sigma_x = sqrt(diag(A^{-1})) + // solve for x and A^{-1} to get x & dx + // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 + // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 + + // If I use the SVD trick to re-condition the matrix, I need to break out the + // kernel and normalization terms from the background term. + // XXX is this true? or was this due to an error in the analysis? + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + + // now pull out the kernel elements into their own square matrix + psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); + psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { + if (ix == bgIndex) continue; + for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { + if (iy == bgIndex) continue; + kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; + ky++; + } + kernelVector->data.F64[kx] = sumVector->data.F64[ix]; + kx++; + } + + psImage *U = NULL; + psImage *V = NULL; + psVector *w = NULL; + if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { + psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); + return NULL; + } + + // calculate A_inverse: + // Ainv = V * w * U^T + psImage *wUt = p_pmSubSolve_wUt (w, U); + psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); + psImage *Xvar = NULL; + psFree (wUt); # ifdef TESTING - // kernel terms: - for (int i = 0; i < w->n; i++) { - fprintf (stderr, "w: %f\n", w->data.F64[i]); - } + // kernel terms: + for (int i = 0; i < w->n; i++) { + fprintf (stderr, "w: %f\n", w->data.F64[i]); + } # endif - // loop over w adding in more and more of the values until chisquare is no longer - // dropping significantly. - // XXX this does not seem to work very well: we seem to need all terms even for - // simple cases... - - psVector *Xsvd = NULL; - { - psVector *Ax = NULL; - psVector *UtB = NULL; - psVector *wUtB = NULL; - - psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); - psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); - psVectorInit (wMask, 1); // start by masking everything - - double chiSquareLast = NAN; - int maxWeight = 0; - - double Axx, Bx, y2; - - // XXX this is an attempt to exclude insignificant modes. - // it was not successful with the ISIS kernel set: removing even - // the least significant mode leaves additional ringing / noise - // because the terms are so coupled. - for (int k = 0; false && (k < w->n); k++) { - - // unmask the k-th weight - wMask->data.U8[k] = 0; - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; - chiSquareLast = chiSquare; - - // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); - if (k && !maxWeight && (deltaChi < 1.0)) { - maxWeight = k; - } - } - - // keep all terms or we get extra ringing - maxWeight = w->n; - psVectorInit (wMask, 1); - for (int k = 0; k < maxWeight; k++) { - wMask->data.U8[k] = 0; - } - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); - - // re-calculate A^{-1} to get new variances: - // Ainv = V * w * U^T - // XXX since we keep all terms, this is identical to Ainv - psImage *wUt = p_pmSubSolve_wUt (wApply, U); - Xvar = p_pmSubSolve_VwUt (V, wUt); - psFree (wUt); - - psFree (Ax); - psFree (UtB); - psFree (wUtB); - psFree (wApply); - psFree (wMask); - } - - // copy the kernel solutions to the full solution vector: - solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { - if (ix == bgIndex) { - solution->data.F64[ix] = 0; - solutionErr->data.F64[ix] = 0.001; - continue; - } - solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); - solution->data.F64[ix] = Xsvd->data.F64[kx]; - kx++; - } - - psFree (kernelMatrix); - psFree (kernelVector); - - psFree (U); - psFree (V); - psFree (w); - - psFree (Ainv); - psFree (Xsvd); - } else { - psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); - if (!luMatrix) { - psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - psFree(luMatrix); - psFree(permutation); - return NULL; - } - - solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); - psFree(luMatrix); - psFree(permutation); - if (!solution) { - psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - return NULL; - } - - // XXX LUD does not provide A^{-1}? fake the error for now - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - for (int ix = 0; ix < sumVector->n; ix++) { - solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; - } - } - - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); - psVectorInit (kernels->solution1, 0.0); - } - - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - int spatialOrder = kernels->spatialOrder; // Order of spatial variation - int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms - for (int i = 0; i < numKernels * numPoly; i++) { - // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); - if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { - // XXX fprintf (stderr, "drop\n"); - kernels->solution1->data.F64[i] = 0.0; - } else { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } - } - // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); - // fprintf (stderr, "chi square Brute = %f\n", chiSquare); - - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); + // loop over w adding in more and more of the values until chisquare is no longer + // dropping significantly. + // XXX this does not seem to work very well: we seem to need all terms even for + // simple cases... + + psVector *Xsvd = NULL; + { + psVector *Ax = NULL; + psVector *UtB = NULL; + psVector *wUtB = NULL; + + psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); + psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); + psVectorInit (wMask, 1); // start by masking everything + + double chiSquareLast = NAN; + int maxWeight = 0; + + double Axx, Bx, y2; + + // XXX this is an attempt to exclude insignificant modes. + // it was not successful with the ISIS kernel set: removing even + // the least significant mode leaves additional ringing / noise + // because the terms are so coupled. + for (int k = 0; false && (k < w->n); k++) { + + // unmask the k-th weight + wMask->data.U8[k] = 0; + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; + chiSquareLast = chiSquare; + + // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); + if (k && !maxWeight && (deltaChi < 1.0)) { + maxWeight = k; + } + } + + // keep all terms or we get extra ringing + maxWeight = w->n; + psVectorInit (wMask, 1); + for (int k = 0; k < maxWeight; k++) { + wMask->data.U8[k] = 0; + } + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); + + // re-calculate A^{-1} to get new variances: + // Ainv = V * w * U^T + // XXX since we keep all terms, this is identical to Ainv + psImage *wUt = p_pmSubSolve_wUt (wApply, U); + Xvar = p_pmSubSolve_VwUt (V, wUt); + psFree (wUt); + + psFree (Ax); + psFree (UtB); + psFree (wUtB); + psFree (wApply); + psFree (wMask); + } + + // copy the kernel solutions to the full solution vector: + solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { + if (ix == bgIndex) { + solution->data.F64[ix] = 0; + solutionErr->data.F64[ix] = 0.001; + continue; + } + solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); + solution->data.F64[ix] = Xsvd->data.F64[kx]; + kx++; + } + + psFree (kernelMatrix); + psFree (kernelVector); + + psFree (U); + psFree (V); + psFree (w); + + psFree (Ainv); + psFree (Xsvd); + } else { + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); + if (!luMatrix) { + psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + + solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); + psFree(luMatrix); + psFree(permutation); + if (!solution) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + return NULL; + } + + // XXX LUD does not provide A^{-1}? fake the error for now + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + for (int ix = 0; ix < sumVector->n; ix++) { + solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; + } + } + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); + psVectorInit (kernels->solution1, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); + if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { + // XXX fprintf (stderr, "drop\n"); + kernels->solution1->data.F64[i] = 0.0; + } else { + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + } + // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); + // fprintf (stderr, "chi square Brute = %f\n", chiSquare); + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); #ifdef TESTING @@ -1299,5 +1297,5 @@ psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64); psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64); - psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *sumVector2 = psVectorAlloc(numParams2, PS_TYPE_F64); psImageInit(sumMatrix1, 0.0); psImageInit(sumMatrix2, 0.0); @@ -1320,7 +1318,19 @@ } + +#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 + + +#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]); + calculatePenalty(sumMatrix2, sumVector2, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]); +#endif psImage *sumMatrix = NULL; // Combined matrix @@ -1330,5 +1340,5 @@ #ifdef TESTING - psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); + psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); { psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); @@ -1341,7 +1351,12 @@ psFitsClose(fits); } + psVectorWriteFile("sumVector.dat", sumVector); #endif psVector *solution = NULL; // Solution to equation! + solution = psVectorAlloc(numParams + numParams2, PS_TYPE_F64); + psVectorInit(solution, 0); + +#if 0 { solution = psMatrixSolveSVD(solution, sumMatrix, sumVector); @@ -1423,26 +1438,28 @@ } } - } - - 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 + { + 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 + + calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2, + sumMatrixX, sumVector1, sumVector2); + } +#endif + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector); @@ -1463,4 +1480,5 @@ psFree(sumVector); + #ifdef TESTING { @@ -1468,4 +1486,5 @@ 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]; } @@ -1530,172 +1549,172 @@ for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest - if (stamp->status != PM_SUBTRACTION_STAMP_USED) { - deviations->data.F32[i] = NAN; - continue; - } - - // Calculate coefficients of the kernel basis functions - polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); - double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation - double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background - - // Calculate residuals - psKernel *weight = stamp->weight; // Weight postage stamp - psImageInit(residual->image, 0.0); - if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { - psKernel *target; // Target postage stamp - psKernel *source; // Source postage stamp - psArray *convolutions; // Convolution postage stamps for each kernel basis function - switch (kernels->mode) { - case PM_SUBTRACTION_MODE_1: - target = stamp->image2; - source = stamp->image1; - convolutions = stamp->convolutions1; - - // Having convolved image1 and changed its normalisation, we need to renormalise the residual - // so that it is on the scale of image1. - psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm, - false); // Kernel image - if (!image) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel."); - return false; - } - double sumKernel = 0; // Sum of kernel, for normalising residual - int size = kernels->size; // Half-size of kernel - int fullSize = 2 * size + 1; // Full size of kernel - for (int y = 0; y < fullSize; y++) { - for (int x = 0; x < fullSize; x++) { - sumKernel += image->data.F32[y][x]; - } - } - psFree(image); - devNorm = 1.0 / sumKernel / numPixels; - break; - case PM_SUBTRACTION_MODE_2: - target = stamp->image1; - source = stamp->image2; - convolutions = stamp->convolutions2; - break; - default: - psAbort("Unsupported subtraction mode: %x", kernels->mode); - } - - for (int j = 0; j < numKernels; j++) { - psKernel *convolution = convolutions->data[j]; // Convolution - double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, - false); // Coefficient - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient; - } - } - } - - // XXX visualize the target, source, convolution and residual - pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i); - - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm; - } - } - } else { - // Dual convolution - psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image - psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image - psKernel *image1 = stamp->image1; // The first image - 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 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] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1; - } - } - } - - // XXX visualize the target, source, convolution and residual - pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i); - - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm; - } - } - } - - double deviation = 0.0; // Sum of differences - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x]; - deviation += dev; + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + deviations->data.F32[i] = NAN; + continue; + } + + // Calculate coefficients of the kernel basis functions + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background + + // Calculate residuals + psKernel *weight = stamp->weight; // Weight postage stamp + psImageInit(residual->image, 0.0); + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { + psKernel *target; // Target postage stamp + psKernel *source; // Source postage stamp + psArray *convolutions; // Convolution postage stamps for each kernel basis function + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + source = stamp->image1; + convolutions = stamp->convolutions1; + + // Having convolved image1 and changed its normalisation, we need to renormalise the residual + // so that it is on the scale of image1. + psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm, + false); // Kernel image + if (!image) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel."); + return false; + } + double sumKernel = 0; // Sum of kernel, for normalising residual + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1; // Full size of kernel + for (int y = 0; y < fullSize; y++) { + for (int x = 0; x < fullSize; x++) { + sumKernel += image->data.F32[y][x]; + } + } + psFree(image); + devNorm = 1.0 / sumKernel / numPixels; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + source = stamp->image2; + convolutions = stamp->convolutions2; + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution + double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, + false); // Coefficient + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += convolution->kernel[y][x] * coefficient; + } + } + } + + // XXX visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i); + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x]; + } + } + } else { + // Dual convolution + psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image + psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image + psKernel *image1 = stamp->image1; // The first image + 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 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] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1; + } + } + } + + // XXX visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i); + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x]; + } + } + } + + double deviation = 0.0; // Sum of differences + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x]; + deviation += dev; #ifdef TESTING - residual->kernel[y][x] = dev; -#endif - } - } - deviations->data.F32[i] = devNorm * deviation; - psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", - i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); - if (!isfinite(deviations->data.F32[i])) { - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - psTrace("psModules.imcombine", 5, - "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", - i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); - continue; - } + residual->kernel[y][x] = dev; +#endif + } + } + deviations->data.F32[i] = devNorm * deviation; + psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); + if (!isfinite(deviations->data.F32[i])) { + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + psTrace("psModules.imcombine", 5, + "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + continue; + } #ifdef TESTING - { - psString filename = NULL; - psStringAppend(&filename, "resid_%03d.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, residual->image, 0, NULL); - psFitsClose(fits); - } - if (stamp->image1) { - psString filename = NULL; - psStringAppend(&filename, "stamp_image1_%03d.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL); - psFitsClose(fits); - } - if (stamp->image2) { - psString filename = NULL; - psStringAppend(&filename, "stamp_image2_%03d.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL); - psFitsClose(fits); - } - if (stamp->weight) { - 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 - + { + psString filename = NULL; + psStringAppend(&filename, "resid_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, residual->image, 0, NULL); + psFitsClose(fits); + } + if (stamp->image1) { + psString filename = NULL; + psStringAppend(&filename, "stamp_image1_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL); + psFitsClose(fits); + } + if (stamp->image2) { + psString filename = NULL; + psStringAppend(&filename, "stamp_image2_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL); + psFitsClose(fits); + } + if (stamp->weight) { + 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 + } // calculate and report the normalization and background for the image center - { - polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0); - double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation - double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background - psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background); - - pmSubtractionVisualShowFit(norm); - pmSubtractionVisualPlotFit(kernels); + { + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0); + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background + psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background); + + pmSubtractionVisualShowFit(norm); + pmSubtractionVisualPlotFit(kernels); } @@ -1716,9 +1735,9 @@ for (int i = 0; i < wUt->numCols; i++) { - for (int j = 0; j < wUt->numRows; j++) { - if (!isfinite(w->data.F64[j])) continue; - if (w->data.F64[j] == 0.0) continue; - wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j]; - } + for (int j = 0; j < wUt->numRows; j++) { + if (!isfinite(w->data.F64[j])) continue; + if (w->data.F64[j] == 0.0) continue; + wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j]; + } } return wUt; @@ -1733,11 +1752,11 @@ for (int i = 0; i < Ainv->numCols; i++) { - for (int j = 0; j < Ainv->numRows; j++) { - double sum = 0.0; - for (int k = 0; k < V->numCols; k++) { - sum += V->data.F64[j][k] * wUt->data.F64[k][i]; - } - Ainv->data.F64[j][i] = sum; - } + for (int j = 0; j < Ainv->numRows; j++) { + double sum = 0.0; + for (int k = 0; k < V->numCols; k++) { + sum += V->data.F64[j][k] * wUt->data.F64[k][i]; + } + Ainv->data.F64[j][i] = sum; + } } return Ainv; @@ -1752,9 +1771,9 @@ for (int i = 0; i < U->numCols; i++) { - double sum = 0.0; - for (int j = 0; j < U->numRows; j++) { - sum += B->data.F64[j] * U->data.F64[j][i]; - } - UtB[0]->data.F64[i] = sum; + double sum = 0.0; + for (int j = 0; j < U->numRows; j++) { + sum += B->data.F64[j] * U->data.F64[j][i]; + } + UtB[0]->data.F64[i] = sum; } return true; @@ -1771,7 +1790,7 @@ for (int i = 0; i < w->n; i++) { - if (!isfinite(w->data.F64[i])) continue; - if (w->data.F64[i] == 0.0) continue; - wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i]; + if (!isfinite(w->data.F64[i])) continue; + if (w->data.F64[i] == 0.0) continue; + wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i]; } return true; @@ -1786,9 +1805,9 @@ for (int j = 0; j < V->numRows; j++) { - double sum = 0.0; - for (int i = 0; i < V->numCols; i++) { - sum += V->data.F64[j][i] * wUtB->data.F64[i]; - } - VwUtB[0]->data.F64[j] = sum; + double sum = 0.0; + for (int i = 0; i < V->numCols; i++) { + sum += V->data.F64[j][i] * wUtB->data.F64[i]; + } + VwUtB[0]->data.F64[j] = sum; } return true; @@ -1803,9 +1822,9 @@ for (int j = 0; j < A->numRows; j++) { - double sum = 0.0; - for (int i = 0; i < A->numCols; i++) { - sum += A->data.F64[j][i] * x->data.F64[i]; - } - B[0]->data.F64[j] = sum; + double sum = 0.0; + for (int i = 0; i < A->numCols; i++) { + sum += A->data.F64[j][i] * x->data.F64[i]; + } + B[0]->data.F64[j] = sum; } return true; @@ -1819,5 +1838,5 @@ double sum = 0.0; for (int i = 0; i < x->n; i++) { - sum += x->data.F64[i] * y->data.F64[i]; + sum += x->data.F64[i] * y->data.F64[i]; } *value = sum; @@ -1832,45 +1851,45 @@ for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; - if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; - - psKernel *weight = NULL; - psKernel *window = NULL; - psKernel *input = NULL; + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; + + psKernel *weight = NULL; + psKernel *window = NULL; + psKernel *input = NULL; #ifdef USE_WEIGHT - weight = stamp->weight; + weight = stamp->weight; #endif #ifdef USE_WINDOW - window = stamps->window; -#endif - - switch (kernels->mode) { - // MODE_1 : convolve image 1 to match image 2 (and vice versa) - case PM_SUBTRACTION_MODE_1: - input = stamp->image2; - break; - case PM_SUBTRACTION_MODE_2: - input = stamp->image1; - break; - default: - psAbort ("programming error"); - } - - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double in = input->kernel[y][x]; - double value = in*in; - if (weight) { - float wtVal = weight->kernel[y][x]; - value *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - value *= winVal; - } - sum += value; - } - } + window = stamps->window; +#endif + + switch (kernels->mode) { + // MODE_1 : convolve image 1 to match image 2 (and vice versa) + case PM_SUBTRACTION_MODE_1: + input = stamp->image2; + break; + case PM_SUBTRACTION_MODE_2: + input = stamp->image1; + break; + default: + psAbort ("programming error"); + } + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double in = input->kernel[y][x]; + double value = in*in; + if (weight) { + float wtVal = weight->kernel[y][x]; + value *= wtVal; + } + if (window) { + float winVal = window->kernel[y][x]; + value *= winVal; + } + sum += value; + } + } } *y2 = sum; @@ -1892,68 +1911,68 @@ for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; - if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; - - psKernel *weight = NULL; - psKernel *window = NULL; - psKernel *target = NULL; - psKernel *source = NULL; - - psArray *convolutions = NULL; + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; + + psKernel *weight = NULL; + psKernel *window = NULL; + psKernel *target = NULL; + psKernel *source = NULL; + + psArray *convolutions = NULL; #ifdef USE_WEIGHT - weight = stamp->weight; + weight = stamp->weight; #endif #ifdef USE_WINDOW - window = stamps->window; -#endif - - switch (kernels->mode) { - // MODE_1 : convolve image 1 to match image 2 (and vice versa) - case PM_SUBTRACTION_MODE_1: - target = stamp->image2; - source = stamp->image1; - convolutions = stamp->convolutions1; - break; - case PM_SUBTRACTION_MODE_2: - target = stamp->image1; - source = stamp->image2; - convolutions = stamp->convolutions2; - break; - default: - psAbort ("programming error"); - } - - // Calculate coefficients of the kernel basis functions - polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); - double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation - double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background - - psImageInit(residual->image, 0.0); - for (int j = 0; j < numKernels; j++) { - psKernel *convolution = convolutions->data[j]; // Convolution - double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient; - } - } - } - - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x]; - double value = PS_SQR(resid); - if (weight) { - float wtVal = weight->kernel[y][x]; - value *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - value *= winVal; - } - sum += value; - } - } + window = stamps->window; +#endif + + switch (kernels->mode) { + // MODE_1 : convolve image 1 to match image 2 (and vice versa) + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + source = stamp->image1; + convolutions = stamp->convolutions1; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + source = stamp->image2; + convolutions = stamp->convolutions2; + break; + default: + psAbort ("programming error"); + } + + // Calculate coefficients of the kernel basis functions + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background + + psImageInit(residual->image, 0.0); + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution + double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient; + } + } + } + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x]; + double value = PS_SQR(resid); + if (weight) { + float wtVal = weight->kernel[y][x]; + value *= wtVal; + } + if (window) { + float winVal = window->kernel[y][x]; + value *= winVal; + } + sum += value; + } + } } psFree (polyValues); @@ -1966,5 +1985,5 @@ for (int i = 0; i < w->n; i++) { - wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i]; + wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i]; } return true; @@ -1981,9 +2000,9 @@ // generate Vn = V * w^{-1} for (int j = 0; j < Vn->numRows; j++) { - for (int i = 0; i < Vn->numCols; i++) { - if (!isfinite(w->data.F64[i])) continue; - if (w->data.F64[i] == 0.0) continue; - Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i]; - } + for (int i = 0; i < Vn->numCols; i++) { + if (!isfinite(w->data.F64[i])) continue; + if (w->data.F64[i] == 0.0) continue; + Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i]; + } } @@ -1993,11 +2012,11 @@ // generate Xvar = Vn * Vn^T for (int j = 0; j < Vn->numRows; j++) { - for (int i = 0; i < Vn->numCols; i++) { - double sum = 0.0; - for (int k = 0; k < Vn->numCols; k++) { - sum += Vn->data.F64[k][i]*Vn->data.F64[k][j]; - } - Xvar->data.F64[j][i] = sum; - } + for (int i = 0; i < Vn->numCols; i++) { + double sum = 0.0; + for (int k = 0; k < Vn->numCols; k++) { + sum += Vn->data.F64[k][i]*Vn->data.F64[k][j]; + } + Xvar->data.F64[j][i] = sum; + } } return Xvar; @@ -2014,5 +2033,5 @@ psFitsWriteImage(fits, header, image, 0, NULL); psFitsClose(fits); - + return true; } Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26546) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26547) @@ -90,6 +90,6 @@ for (int i = 0, x = -size; x <= size; i++, x++) { - float xf = x / sigma; - float z = -0.25*xf*xf; + float xf = x / sigma; + float z = -0.25*xf*xf; kernel->data.F32[i] = norm * p_pmSubtractionHermitianPolynomial(xf, order) * exp(z); } @@ -132,7 +132,5 @@ kernels->preCalc->data[index] = NULL; kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); - if (!isfinite(kernels->penalties->data.F32[index])) { - psAbort ("invalid penalty"); - } + psAssert(isfinite(kernels->penalties->data.F32[index]), "Invalid penalty"); psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); @@ -154,8 +152,8 @@ double moment = 0.0; // Moment, for penalty for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - double value = preCalc->kernel->kernel[v][u]; - moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); - } + for (int u = -size; u <= size; u++) { + double value = preCalc->kernel->kernel[v][u]; + moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); + } } @@ -171,9 +169,9 @@ for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - sum += preCalc->kernel->kernel[v][u]; - min = PS_MIN(preCalc->kernel->kernel[v][u], min); - max = PS_MAX(preCalc->kernel->kernel[v][u], max); - } + for (int u = -size; u <= size; u++) { + sum += preCalc->kernel->kernel[v][u]; + min = PS_MIN(preCalc->kernel->kernel[v][u], min); + max = PS_MAX(preCalc->kernel->kernel[v][u], max); + } } #if 0 @@ -183,7 +181,7 @@ // only even terms have non-zero sums if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) { - moment /= sum; + moment /= sum; } else { - moment = 0.0; + moment = 0.0; } @@ -193,13 +191,13 @@ if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) { - zeroNull = true; + zeroNull = true; } if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) { - zeroNull = true; + zeroNull = true; } if ((uOrder % 2) || (vOrder % 2)) { - // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max); - scale2D = 1.0 / max; - scale1D = sqrt(scale2D); + // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max); + scale2D = 1.0 / max; + scale1D = sqrt(scale2D); } @@ -208,5 +206,5 @@ psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32)); if (zeroNull) { - preCalc->kernel->kernel[0][0] -= 1.0; + preCalc->kernel->kernel[0][0] -= 1.0; } @@ -216,9 +214,9 @@ max = FLT_MIN; for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - sum += preCalc->kernel->kernel[v][u]; - min = PS_MIN(preCalc->kernel->kernel[v][u], min); - max = PS_MAX(preCalc->kernel->kernel[v][u], max); - } + for (int u = -size; u <= size; u++) { + sum += preCalc->kernel->kernel[v][u]; + min = PS_MIN(preCalc->kernel->kernel[v][u], min); + max = PS_MAX(preCalc->kernel->kernel[v][u], max); + } } fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D); @@ -229,14 +227,14 @@ kernels->v->data.S32[index] = vOrder; if (kernels->preCalc->data[index]) { - psFree(kernels->preCalc->data[index]); + psFree(kernels->preCalc->data[index]); } kernels->preCalc->data[index] = preCalc; kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); if (!isfinite(kernels->penalties->data.F32[index])) { - psAbort ("invalid penalty"); + psAbort ("invalid penalty"); } psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, - fwhm, uOrder, vOrder, fabsf(moment)); + fwhm, uOrder, vOrder, fabsf(moment)); return true; @@ -259,7 +257,7 @@ psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32); for (int i = 0; i < fwhmsIN->n; i++) { - if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; - psVectorAppend(fwhms, fwhmsIN->data.F32[i]); - psVectorAppend(orders, ordersIN->data.S32[i]); + if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; + psVectorAppend(fwhms, fwhmsIN->data.F32[i]); + psVectorAppend(orders, ordersIN->data.S32[i]); } @@ -288,7 +286,6 @@ for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { - pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); } } @@ -299,6 +296,6 @@ pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder, - const psVector *fwhmsIN, const psVector *ordersIN, - float penalty, pmSubtractionMode mode) + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); @@ -314,7 +311,7 @@ psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32); for (int i = 0; i < fwhmsIN->n; i++) { - if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; - psVectorAppend(fwhms, fwhmsIN->data.F32[i]); - psVectorAppend(orders, ordersIN->data.S32[i]); + if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; + psVectorAppend(fwhms, fwhmsIN->data.F32[i]); + psVectorAppend(orders, ordersIN->data.S32[i]); } @@ -342,5 +339,5 @@ for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); } } @@ -351,6 +348,6 @@ pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder, - const psVector *fwhmsIN, const psVector *ordersIN, - float penalty, pmSubtractionMode mode) + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); @@ -366,7 +363,7 @@ psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32); for (int i = 0; i < fwhmsIN->n; i++) { - if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; - psVectorAppend(fwhms, fwhmsIN->data.F32[i]); - psVectorAppend(orders, ordersIN->data.S32[i]); + if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; + psVectorAppend(fwhms, fwhmsIN->data.F32[i]); + psVectorAppend(orders, ordersIN->data.S32[i]); } @@ -405,20 +402,20 @@ pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values - // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel - psKernel *kernelTarget = preCalc->kernel; + // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel + psKernel *kernelTarget = preCalc->kernel; preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel - // XXX do we use Alard-Lupton normalization (last param true) or not? - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); - - // XXXX test demo that deconvolved kernel is valid + // XXX do we use Alard-Lupton normalization (last param true) or not? + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); + + // XXXX test demo that deconvolved kernel is valid # if 1 - psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); - psArrayAdd (deconKernels, 100, kernelConv); - psFree (kernelConv); - - if (!uOrder && !vOrder){ - pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); - } + psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); + psArrayAdd (deconKernels, 100, kernelConv); + psFree (kernelConv); + + if (!uOrder && !vOrder){ + pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); + } # endif } @@ -429,17 +426,17 @@ psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32); for (int i = 0; i < deconKernels->n; i++) { - for (int j = 0; j <= i; j++) { - psImage *t1 = deconKernels->data[i]; - psImage *t2 = deconKernels->data[j]; - - double sum = 0.0; - for (int iy = 0; iy < t1->numRows; iy++) { - for (int ix = 0; ix < t1->numCols; ix++) { - sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix]; - } - } - dot->data.F32[j][i] = sum; - dot->data.F32[i][j] = sum; - } + for (int j = 0; j <= i; j++) { + psImage *t1 = deconKernels->data[i]; + psImage *t2 = deconKernels->data[j]; + + double sum = 0.0; + for (int iy = 0; iy < t1->numRows; iy++) { + for (int ix = 0; ix < t1->numCols; ix++) { + sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix]; + } + } + dot->data.F32[j][i] = sum; + dot->data.F32[i][j] = sum; + } } pmSubtractionVisualShowSubtraction (dot, NULL, NULL); @@ -494,36 +491,36 @@ switch (type) { case PM_SUBTRACTION_KERNEL_ISIS: - preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size); - preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size); - preCalc->uCoords = NULL; - preCalc->vCoords = NULL; - preCalc->poly = NULL; - break; + preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size); + preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size); + preCalc->uCoords = NULL; + preCalc->vCoords = NULL; + preCalc->poly = NULL; + break; case PM_SUBTRACTION_KERNEL_HERM: - preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); - preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); - preCalc->uCoords = NULL; - preCalc->vCoords = NULL; - preCalc->poly = NULL; - break; + preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); + preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); + preCalc->uCoords = NULL; + preCalc->vCoords = NULL; + preCalc->poly = NULL; + break; case PM_SUBTRACTION_KERNEL_DECONV_HERM: - preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); - preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); - preCalc->uCoords = NULL; - preCalc->vCoords = NULL; - preCalc->poly = NULL; - break; + preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); + preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); + preCalc->uCoords = NULL; + preCalc->vCoords = NULL; + preCalc->poly = NULL; + break; case PM_SUBTRACTION_KERNEL_RINGS: - // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure - // we allocate these vectors here, but leave the kernel generation to the main function - preCalc->xKernel = NULL; - preCalc->yKernel = NULL; - preCalc->kernel = NULL; - preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords - preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords - preCalc->poly = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial - return preCalc; + // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure + // we allocate these vectors here, but leave the kernel generation to the main function + preCalc->xKernel = NULL; + preCalc->yKernel = NULL; + preCalc->kernel = NULL; + preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords + preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords + preCalc->poly = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial + return preCalc; default: - psAbort("programming error: invalid type for PreCalc kernel"); + psAbort("programming error: invalid type for PreCalc kernel"); } @@ -532,9 +529,9 @@ // generate 2D kernel from 1D realizations for (int v = -size, y = 0; v <= size; v++, y++) { - for (int u = -size, x = 0; u <= size; u++, x++) { - preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel - } - } - + for (int u = -size, x = 0; u <= size; u++, x++) { + preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel + } + } + return preCalc; } @@ -864,5 +861,5 @@ for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) { - pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0); + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0); double moment = 0.0; // Moment, for penalty @@ -870,9 +867,9 @@ // Central pixel is easy preCalc->uCoords->data.S32[0] = 0; - preCalc->vCoords->data.S32[0] = 0; + preCalc->vCoords->data.S32[0] = 0; preCalc->poly->data.F32[0] = 1.0; preCalc->uCoords->n = 1; - preCalc->vCoords->n = 1; - preCalc->poly->n = 1; + preCalc->vCoords->n = 1; + preCalc->poly->n = 1; radiusLast = 0; moment = 0.0; @@ -931,7 +928,7 @@ kernels->v->data.S32[index] = vOrder; kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); - if (!isfinite(kernels->penalties->data.F32[index])) { - psAbort ("invalid penalty"); - } + if (!isfinite(kernels->penalties->data.F32[index])) { + psAbort ("invalid penalty"); + } psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, @@ -1026,44 +1023,44 @@ type = PM_SUBTRACTION_KERNEL_GUNK; psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description); - } - - type = pmSubtractionKernelsTypeFromString (description); - psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must be ISIS, HERM or DECONV_HERM"); - - char *ptr = NULL; - switch (type) { - case PM_SUBTRACTION_KERNEL_ISIS: - case PM_SUBTRACTION_KERNEL_HERM: - ptr = (char*) description + 5; // Eat "ISIS(" or "HERM(" - break; - case PM_SUBTRACTION_KERNEL_DECONV_HERM: - ptr = (char*) description + 12; // Eat "DECONV_HERM(" - break; - default: - psAbort("programming error: invalid kernel type"); - } - PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); - - // Count the number of Gaussians - int numGauss = 0; - for (char *string = ptr; string; string = strchr(string + 1, '(')) { - numGauss++; - } - - fwhms = psVectorAlloc(numGauss, PS_TYPE_F32); - orders = psVectorAlloc(numGauss, PS_TYPE_S32); - - for (int i = 0; i < numGauss; i++) { - ptr++; // Eat the '(' - PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234," - PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)" - } - - ptr++; // Eat ',' - PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); - penalty = parseStringFloat(ptr); - - return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); - } + } + + type = pmSubtractionKernelsTypeFromString (description); + psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must be ISIS, HERM or DECONV_HERM"); + + char *ptr = NULL; + switch (type) { + case PM_SUBTRACTION_KERNEL_ISIS: + case PM_SUBTRACTION_KERNEL_HERM: + ptr = (char*) description + 5; // Eat "ISIS(" or "HERM(" + break; + case PM_SUBTRACTION_KERNEL_DECONV_HERM: + ptr = (char*) description + 12; // Eat "DECONV_HERM(" + break; + default: + psAbort("programming error: invalid kernel type"); + } + PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); + + // Count the number of Gaussians + int numGauss = 0; + for (char *string = ptr; string; string = strchr(string + 1, '(')) { + numGauss++; + } + + fwhms = psVectorAlloc(numGauss, PS_TYPE_F32); + orders = psVectorAlloc(numGauss, PS_TYPE_S32); + + for (int i = 0; i < numGauss; i++) { + ptr++; // Eat the '(' + PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234," + PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)" + } + + ptr++; // Eat ',' + PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); + penalty = parseStringFloat(ptr); + + return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); + } if (strncmp(description, "RINGS", 5) == 0) { @@ -1075,6 +1072,6 @@ PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); - return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); - } + return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); + } psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported."); Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26546) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26547) @@ -242,23 +242,23 @@ for (int i = 0; i < kernels->num; i++) { - pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; - psKernel *kernel = preCalc->kernel; - - int xSub = i % NXsub; - int ySub = i / NXsub; - - int xPix = xSub * (2*footprint + 1 + 3) + footprint; - int yPix = ySub * (2*footprint + 1 + 3) + footprint; - - double sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; - sum += kernel->kernel[y][x]; - } - } - fprintf (stderr, "kernel %d, sum %f\n", i, sum); - } - + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; + psKernel *kernel = preCalc->kernel; + + int xSub = i % NXsub; + int ySub = i / NXsub; + + int xPix = xSub * (2*footprint + 1 + 3) + footprint; + int yPix = ySub * (2*footprint + 1 + 3) + footprint; + + double sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; + sum += kernel->kernel[y][x]; + } + } + fprintf (stderr, "kernel %d, sum %f\n", i, sum); + } + pmVisualScaleImage(kapa1, output, "Image", 0, true); pmVisualAskUser(&plotImage); @@ -280,20 +280,20 @@ float maxFlux = NAN; for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; - if (!isfinite(stamp->flux)) continue; - if (!stamp->convolutions1 && !stamp->convolutions2) continue; - if (!maxStamp) { - maxFlux = stamp->flux; - maxStamp = stamp; - continue; - } - if (stamp->flux > maxFlux) { - maxFlux = stamp->flux; - maxStamp = stamp; - } + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (!isfinite(stamp->flux)) continue; + if (!stamp->convolutions1 && !stamp->convolutions2) continue; + if (!maxStamp) { + maxFlux = stamp->flux; + maxStamp = stamp; + continue; + } + if (stamp->flux > maxFlux) { + maxFlux = stamp->flux; + maxStamp = stamp; + } } if (!isfinite(maxStamp->flux)) { - fprintf (stderr, "no valid stamps?\n"); + fprintf (stderr, "no valid stamps?\n"); } @@ -301,71 +301,71 @@ if (maxStamp->convolutions1) { - // output image is a grid of NXsub by NYsub sub-images - nKernels = maxStamp->convolutions1->n; - int NXsub = sqrt(nKernels); - int NYsub = nKernels / NXsub; - if (nKernels % NXsub) NYsub++; - - int NXpix = NXsub * (2*footprint + 1 + 3); - int NYpix = NYsub * (2*footprint + 1 + 3); - - psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); - psImageInit (output, 0.0); - - for (int i = 0; i < nKernels; i++) { + // output image is a grid of NXsub by NYsub sub-images + nKernels = maxStamp->convolutions1->n; + int NXsub = sqrt(nKernels); + int NYsub = nKernels / NXsub; + if (nKernels % NXsub) NYsub++; + + int NXpix = NXsub * (2*footprint + 1 + 3); + int NYpix = NYsub * (2*footprint + 1 + 3); + + psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); + psImageInit (output, 0.0); + + for (int i = 0; i < nKernels; i++) { psKernel *kernel = maxStamp->convolutions1->data[i]; - - int xSub = i % NXsub; - int ySub = i / NXsub; - - int xPix = xSub * (2*footprint + 1 + 3) + footprint; - int yPix = ySub * (2*footprint + 1 + 3) + footprint; - - double sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; - sum += kernel->kernel[y][x]; - } - } - fprintf (stderr, "kernel %d, sum %f\n", i, sum); - } - pmVisualScaleImage(kapa2, output, "Image", 0, true); - } - + + int xSub = i % NXsub; + int ySub = i / NXsub; + + int xPix = xSub * (2*footprint + 1 + 3) + footprint; + int yPix = ySub * (2*footprint + 1 + 3) + footprint; + + double sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; + sum += kernel->kernel[y][x]; + } + } + fprintf (stderr, "kernel %d, sum %f\n", i, sum); + } + pmVisualScaleImage(kapa2, output, "Image", 0, true); + } + if (maxStamp->convolutions2) { - // output image is a grid of NXsub by NYsub sub-images - nKernels = maxStamp->convolutions2->n; - int NXsub = sqrt(nKernels); - int NYsub = nKernels / NXsub; - if (nKernels % NXsub) NYsub++; - - int NXpix = NXsub * (2*footprint + 1 + 3); - int NYpix = NYsub * (2*footprint + 1 + 3); - - psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); - psImageInit (output, 0.0); - - for (int i = 0; i < nKernels; i++) { + // output image is a grid of NXsub by NYsub sub-images + nKernels = maxStamp->convolutions2->n; + int NXsub = sqrt(nKernels); + int NYsub = nKernels / NXsub; + if (nKernels % NXsub) NYsub++; + + int NXpix = NXsub * (2*footprint + 1 + 3); + int NYpix = NYsub * (2*footprint + 1 + 3); + + psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); + psImageInit (output, 0.0); + + for (int i = 0; i < nKernels; i++) { psKernel *kernel = maxStamp->convolutions2->data[i]; - - int xSub = i % NXsub; - int ySub = i / NXsub; - - int xPix = xSub * (2*footprint + 1 + 3) + footprint; - int yPix = ySub * (2*footprint + 1 + 3) + footprint; - - double sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; - sum += kernel->kernel[y][x]; - } - } - fprintf (stderr, "kernel %d, sum %f\n", i, sum); - } - pmVisualScaleImage(kapa2, output, "Image", 1, true); - } - + + int xSub = i % NXsub; + int ySub = i / NXsub; + + int xPix = xSub * (2*footprint + 1 + 3) + footprint; + int yPix = ySub * (2*footprint + 1 + 3) + footprint; + + double sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; + sum += kernel->kernel[y][x]; + } + } + fprintf (stderr, "kernel %d, sum %f\n", i, sum); + } + pmVisualScaleImage(kapa2, output, "Image", 1, true); + } + pmVisualAskUser(&plotImage); return true; @@ -394,12 +394,12 @@ overlay[Noverlay].type = KII_OVERLAY_BOX; - if ((stamp->x < 1.0) && (stamp->y < 1.0)) { - // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y); - continue; - } - if (!isfinite(stamp->x) && !isfinite(stamp->y)) { - // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y); - continue; - } + if ((stamp->x < 1.0) && (stamp->y < 1.0)) { + // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y); + continue; + } + if (!isfinite(stamp->x) && !isfinite(stamp->y)) { + // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y); + continue; + } overlay[Noverlay].x = stamp->x; overlay[Noverlay].y = stamp->y; @@ -437,5 +437,5 @@ float NXf = sqrt(stamps->num); NX = (int) NXf == NXf ? NXf : NXf + 1.0; - + float NYf = stamps->num / NX; NY = (int) NYf == NY ? NYf : NYf + 1.0; @@ -453,5 +453,5 @@ differenceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - + psImageInit (sourceImage, 0.0); psImageInit (targetImage, 0.0); @@ -479,8 +479,8 @@ sum = 0.0; for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x]; - sum += targetImage->data.F32[y + NYpix][x + NXpix]; - } + for (int x = -footprint; x <= footprint; x++) { + targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x]; + sum += targetImage->data.F32[y + NYpix][x + NXpix]; + } } targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; @@ -489,8 +489,8 @@ sum = 0.0; for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x]; - sum += sourceImage->data.F32[y + NYpix][x + NXpix]; - } + for (int x = -footprint; x <= footprint; x++) { + sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x]; + sum += sourceImage->data.F32[y + NYpix][x + NXpix]; + } } sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; @@ -499,18 +499,18 @@ sum = 0.0; for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x]; - sum += convolutionImage->data.F32[y + NYpix][x + NXpix]; - } + for (int x = -footprint; x <= footprint; x++) { + convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x]; + sum += convolutionImage->data.F32[y + NYpix][x + NXpix]; + } } convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; - + // insert the (difference) kernel into the (difference) image: sum = 0.0; for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm; - sum += differenceImage->data.F32[y + NYpix][x + NXpix]; - } + for (int x = -footprint; x <= footprint; x++) { + differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm; + sum += differenceImage->data.F32[y + NYpix][x + NXpix]; + } } differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; @@ -519,8 +519,8 @@ sum = 0.0; for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm + convolution->kernel[y][x]; - sum += residualImage->data.F32[y + NYpix][x + NXpix]; - } + for (int x = -footprint; x <= footprint; x++) { + residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm + convolution->kernel[y][x]; + sum += residualImage->data.F32[y + NYpix][x + NXpix]; + } } residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; @@ -528,7 +528,7 @@ // insert the (fresidual) kernel into the (fresidual) image: for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); - } + for (int x = -footprint; x <= footprint; x++) { + fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); + } } return true; @@ -536,16 +536,17 @@ bool pmSubtractionVisualShowFit(double norm) { + + if (!pmVisualIsVisual()) return true; // XXX a dumb test : dump the residual image and exit { - psMetadata *header = psMetadataAlloc(); + psMetadata *header = psMetadataAlloc(); psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm); - psFits *fits = psFitsOpen("resid.fits", "w"); - psFitsWriteImage(fits, header, residualImage, 0, NULL); - psFitsClose(fits); - exit (0); - } - - // if (!pmVisualIsVisual()) return true; + psFits *fits = psFitsOpen("resid.fits", "w"); + psFitsWriteImage(fits, header, residualImage, 0, NULL); + psFitsClose(fits); + exit (0); + } + if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false; if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false; @@ -606,5 +607,5 @@ for (int i = 0; i < kernels->num; i++) { x->data.F32[i] = i; - y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false); + y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false); graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]); graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);