Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26347) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26373) @@ -21,4 +21,5 @@ #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" +#include "pmSubtractionVisual.h" #include "pmSubtractionThreads.h" @@ -525,4 +526,7 @@ int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version + + // pmSubtractionVisualShowSubtraction(image->image, kernel->image, conv); + psFree(conv); return convolved; Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26347) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26373) @@ -831,4 +831,21 @@ } +// private functions used on pmSubtractionSolveEquation +bool psVectorWriteFile (char *filename, const psVector *vector); +bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header); + +psImage *p_pmSubSolve_wUt (psVector *w, psImage *U); +psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt); + +psVector *p_pmSubSolve_UtB (psImage *U, psVector *B); +psVector *p_pmSubSolve_wUtB (psVector *w, psVector *UtB); +psVector *p_pmSubSolve_VwUtB (psImage *V, psVector *wUtB); + +psVector *p_pmSubSolve_Ax (psImage *A, psVector *x); + +double p_pmSubSolve_y2 (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); +double p_pmSubSolve_VdV (psVector *x, psVector *y); +double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); + bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps, @@ -943,7 +960,5 @@ { psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); - psFits *fits = psFitsOpen("matrixInv.fits", "w"); - psFitsWriteImage(fits, NULL, inverse, 0, NULL); - psFitsClose(fits); + psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); psFree(inverse); } @@ -952,7 +967,5 @@ psImage *Xt = psMatrixTranspose(NULL, X); psImage *XtX = psMatrixMultiply(NULL, Xt, X); - psFits *fits = psFitsOpen("matrixErr.fits", "w"); - psFitsWriteImage(fits, NULL, XtX, 0, NULL); - psFitsClose(fits); + psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); psFree(X); psFree(Xt); @@ -961,25 +974,22 @@ #endif - // XXX test: save the matrix A and vector b: - { - psFits *fits = psFitsOpen("matrix.fits", "w"); - psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL); - psFitsClose(fits); - - FILE *f = fopen ("vector.dat", "w"); - int fd = fileno(f); - p_psVectorPrint (fd, sumVector, "B"); - fclose (f); - } - # define SVD_ANALYSIS 1 - // re-do this section with SVD: +# define COEFF_SIG 0.0 +# define SVD_TOL 1e-9 + + psVector *solution = NULL; + psVector *solutionErr = NULL; + + // Use SVD to determine the kernel coeffs (and validate) if (SVD_ANALYSIS) { - // We have sumVector and sumMatrix. we are trying to solve the equation: sumMatrix * x = sumVector. + + // 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 determiend or unstable. If not numerically ill-conditioned, the + // 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 + // 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. @@ -993,79 +1003,169 @@ // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 - // I get confused by the index values between the image vs matrix usage: In terms - // of the elements of an image A(x,y) = A->data.F32[y][x] = A_x,y, a matrix - // multiplication is: A_k,j * B_i,k = C_i,j + // XXX TEST with known matrix: + // psFree (sumMatrix); + // sumMatrix = psImageAlloc (2, 2, PS_TYPE_F64); + // sumMatrix->data.F64[0][0] = 1.0; + // sumMatrix->data.F64[0][1] = 0.5; + // + // sumMatrix->data.F64[1][0] = 0.5; + // sumMatrix->data.F64[1][1] = 2.0; + + // sumMatrix->data.F64[2][0] = 0.0; + // sumMatrix->data.F64[2][1] = 1.0; + + // TEST psFitsWriteImageSimple("A.fits", sumMatrix, NULL); + // TEST psVectorWriteFile("B.dat", sumVector);; + + // If I use the SVD trick to re-condition the matrix, I need to break out the + // kernel terms from the normalization and background terms. + + // XX int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + + // solve for the normalization first (independent of all other terms) + // XX double sumRR = sumMatrix->data.F64[normIndex][normIndex]; + // XX double sumIR = sumVector->data.F64[normIndex]; + // XX double norm = sumIR / sumRR; + + // 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); + + // note that we need to subtract norm*sumRC from the sumVector elements (since we + // are fitting to the residual image of I - norm*Ref) + for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { + // if (ix == normIndex) continue; + if (ix == bgIndex) continue; + for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { + // if (iy == normIndex) continue; + if (iy == bgIndex) continue; + kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; + ky++; + } + // kernelVector->data.F64[kx] = sumVector->data.F64[ix] - norm*sumMatrix->data.F64[normIndex][ix]; + kernelVector->data.F64[kx] = sumVector->data.F64[ix]; + kx++; + } + + psFitsWriteImageSimple("A.fits", kernelMatrix, NULL); psImage *U = NULL; psImage *V = NULL; psVector *w = NULL; - if (!psMatrixSVD (&U, &w, &V, sumMatrix)) { + if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); return NULL; } - // calculate A_inverse and x_svd - // A'_i,j = \sum_k (V_k,j * w_k * Ut_i,k) - // but U^T_i,j = U_j,i, so: - // A'_i,j = \sum_k (V_k,j * w_k * U_k,i) - - // v1_j = \sum_k (UT_k,j * B_k) - // v1_j = w_j * \sum_k (U_j,k * B_k) - // x_i = \sum_j V_j,i * w_j * \sum_k (U_j,k * B_k) - - // A'_i,j = \sum_k (V_k,j * w_k * U_k,i) - - psImage *Ainv = psImageAlloc(A->numCols, A->numRows, PS_TYPE_F32); - psVector *Xsvd = psVectorAlloc(A->numCols, PS_TYPE_F32); - for (int iy = 0; iy < A->numRows; iy++) { - double vsum = 0.0; - for (int ix = 0; ix < A->numCols; ix++) { - double msum = 0.0; - double vsum1 = 0.0; - for (int k = 0; k < A->numCols; k++) { - if (fabs(w->data.F32[k]) < FLT_EPSILON) continue; - msum += V->data.F32[iy][k] * U->data.F32[ix][k] / w->data.F32[k]; - vsum1 += U->data.F32[k][iy] * B->data.F32[k]; - } - Ainv->data.F32[iy][ix] = sum; - if (fabs(w->data.F32[ix]) < FLT_EPSILON) continue; - vsum += V->data.F32[iy][ix] * vsum1 / w->data.F32[ix]; - } - Xsvd->data.F32[iy] = vsum; - } - - // compare Xsvd[k] to dXsvd and zero out w entries that are not significant - # define COEFF_SIG 1.0 - for (int k = 0; k < A->numRows; k++) { - float dXsvd = sqrt(Ainv->data.F32[k][k]); - if (fabs(Xsvd->data.F32[k] / dXsvd) < COEFF_SIG) { - w->data.F32[k] = 0.0; - } - } - - // - - - } - - psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); - psFree(sumMatrix); - if (!luMatrix) { - psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - return NULL; - } - - psVector *solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - if (!solution) { - psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); - return NULL; - } + // TEST psFitsWriteImageSimple("U.fits", U, NULL); + // TEST psFitsWriteImageSimple("V.fits", V, NULL); + // TEST psVectorWriteFile("w.dat", w); + + // calculate A_inverse: + // Ainv = V * w * U^T + psImage *wUt = p_pmSubSolve_wUt (w, U); + psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); + + // TEST psFitsWriteImageSimple("AInv.fits", Ainv, NULL); + // TEST psImage *Io = p_pmSubSolve_VwUt (Ainv, sumMatrix); + // TEST psFitsWriteImageSimple("I.fits", Io, NULL); + + // kernel terms: + for (int i = 0; i < w->n; i++) { + fprintf (stderr, "w: %f\n", w->data.F64[i]); + } + + // XXX loop over w adding in more and more of the values until chisquare is no + // longer dropping significantly + + // zero-out ill-conditioned SVD components: + // largest w value is in w->data.F64[0]: + double Wmax = w->data.F64[0]; + for (int k = 1; k < w->n; k++) { + if (fabs(w->data.F64[k] / Wmax) < SVD_TOL) { + fprintf (stderr, "masking basis function %d : %f\n", k, w->data.F64[k] / Wmax); + w->data.F64[k] = NAN; + psTrace("psModules.imcombine", 5, "masking basis function %d : %f", k, w->data.F64[k] / Wmax); + } + } + + // solve for x: + // x = V * w * (U^T * B) + psVector *UtB = p_pmSubSolve_UtB (U, kernelVector); + psVector *wUtB = p_pmSubSolve_wUtB (w, UtB); + psVector *Xsvd = p_pmSubSolve_VwUtB (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 + psVector *Ax = p_pmSubSolve_Ax (kernelMatrix, Xsvd); + double Axx = p_pmSubSolve_VdV (Ax, Xsvd); + double Bx = p_pmSubSolve_VdV (kernelVector, Xsvd); + double y2 = p_pmSubSolve_y2 (kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + fprintf (stderr, "chi square = %f\n", chiSquare); + + // 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++) { + // XX if (ix == normIndex) { + // XX // solution->data.F64[ix] = norm; // + // XX solution->data.F64[ix] = 1.0; // XXX a CHEAT + // XX solutionErr->data.F64[ix] = 0.001; + // XX continue; + // XX } + 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++; + } + + // TEST psVectorWriteFile("X.dat", Xsvd); + // TEST psVector *Bo = p_pmSubSolve_Ax (sumMatrix, Xsvd); + // TEST psVectorWriteFile("Bo.dat", Bo); + + psFree (kernelMatrix); + psFree (kernelVector); + + psFree (U); + psFree (V); + psFree (w); + + psFree (wUt); + psFree (Ainv); + + psFree (UtB); + psFree (wUtB); + psFree (Xsvd); + } else { + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); + psFree(sumMatrix); + if (!luMatrix) { + psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + + solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + if (!solution) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + return NULL; + } + } if (!kernels->solution1) { @@ -1088,7 +1188,17 @@ int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms for (int i = 0; i < numKernels * numPoly; i++) { - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } + 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) { + fprintf (stderr, "drop\n"); + kernels->solution1->data.F64[i] = 0.0; + } else { + 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); @@ -1141,9 +1251,5 @@ #ifdef TESTING - { - psFits *fits = psFitsOpen("sumMatrix.fits", "w"); - psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL); - psFitsClose(fits); - } + psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); { psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); @@ -1345,155 +1451,155 @@ 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; - } - } - } + 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; - } - } - } - 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; + 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; + } + } + } + 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; #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); - } + { + 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 @@ -1516,2 +1622,282 @@ return deviations; } + +// we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w) +psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) { + + psAssert (w->n == U->numCols, "w and U dimensions do not match"); + + // wUt has dimensions transposed relative to Ut. + psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64); + psImageInit (wUt, 0.0); + + 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]; + } + } + return wUt; +} + +// XXX this is just standard matrix multiplication: use psMatrixMultiply? +psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) { + + psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match"); + + psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64); + + 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; + } + } + return Ainv; +} + +// we are supplied U, not Ut +psVector *p_pmSubSolve_UtB (psImage *U, psVector *B) { + + psAssert (U->numRows == B->n, "U and B dimensions do not match"); + + psVector *UtB = psVectorAlloc (U->numCols, PS_TYPE_F64); + + 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->data.F64[i] = sum; + } + return UtB; +} + +// w is diagonal +psVector *p_pmSubSolve_wUtB (psVector *w, psVector *UtB) { + + psAssert (w->n == UtB->n, "w and UtB dimensions do not match"); + + // wUt has dimensions transposed relative to Ut. + psVector *wUtB = psVectorAlloc (w->n, PS_TYPE_F64); + psVectorInit (wUtB, 0.0); + + for (int i = 0; i < w->n; i++) { + if (!isfinite(w->data.F64[i])) continue; + if (w->data.F64[i] == 0.0) continue; + wUtB->data.F64[i] = UtB->data.F64[i] / w->data.F64[i]; + } + return wUtB; +} + +// this is basically matrix * vector +psVector *p_pmSubSolve_VwUtB (psImage *V, psVector *wUtB) { + + psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match"); + + psVector *VwUtB = psVectorAlloc (V->numRows, PS_TYPE_F64); + + 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->data.F64[j] = sum; + } + return VwUtB; +} + +// this is basically matrix * vector +psVector *p_pmSubSolve_Ax (psImage *A, psVector *x) { + + psAssert (A->numCols == x->n, "A and x dimensions do not match"); + + psVector *B = psVectorAlloc (A->numRows, PS_TYPE_F64); + + 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->data.F64[j] = sum; + } + return B; +} + +// this is basically Vector * vector +double p_pmSubSolve_VdV (psVector *x, psVector *y) { + + psAssert (x->n == y->n, "x and y dimensions do not match"); + + double sum = 0.0; + for (int i = 0; i < x->n; i++) { + sum += x->data.F64[i] * y->data.F64[i]; + } + return sum; +} + +double p_pmSubSolve_y2 (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { + + int footprint = stamps->footprint; // Half-size of stamps + + double sum = 0.0; + 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; + +#ifdef USE_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; + } + } + } + return sum; +} + +double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { + + int footprint = stamps->footprint; // Half-size of stamps + int numKernels = kernels->num; // Number of kernels + + double sum = 0.0; + + psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + psImageInit(residual->image, 0.0); + + psImage *polyValues = NULL; // Polynomial values + + 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; + +#ifdef USE_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; + } + } + } + psFree (polyValues); + psFree (residual); + + return sum; +} + +// I get confused by the index values between the image vs matrix usage: In terms +// of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix +// multiplication is: A_k,j * B_i,k = C_i,j + + +bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) { + + psFits *fits = psFitsOpen(filename, "w"); + psFitsWriteImage(fits, header, image, 0, NULL); + psFitsClose(fits); + + return true; +} + +bool psVectorWriteFile (char *filename, const psVector *vector) { + + FILE *f = fopen (filename, "w"); + int fd = fileno(f); + p_psVectorPrint (fd, vector, "unnamed"); + fclose (f); + + return true; +} + Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26347) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26373) @@ -646,4 +646,6 @@ psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); + double maxValue = 0.0; + // generate the window pixels for (int y = -size; y <= size; y++) { @@ -666,4 +668,14 @@ } 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; } } Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26347) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26373) @@ -308,4 +308,6 @@ // if (!pmVisualIsVisual()) return true; + double sum; + int NXoff = index % NX; int NYoff = index / NX; @@ -315,42 +317,57 @@ // insert the (target) kernel into the (target) image: + 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]; + } + } + targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; // insert the (source) kernel into the (source) image: + 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]; + } + } + sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; + // insert the (convolution) kernel into the (convolution) image: + 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]; + } + } + 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]; + } + } + differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; // insert the (residual) kernel into the (residual) image: + 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]; + } + } + residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; // 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] = PS_SQR(residualImage->data.F32[y + NYpix][x + NXpix]) / (target->kernel[y][x] + 25.0); + fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); } }