Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (revision 30288) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (revision 30289) @@ -20,5 +20,9 @@ //# define USE_WEIGHT // Include weight (1/variance) in equation? # define USE_WINDOW // window to avoid neighbor contamination -// XXX if we want to have a weight and window, we'll need to pass through to the chisq function + +/* I believe we want to apply the WEIGHT to the chisq portions of the calculation (but not the WINDOW), + * and the WINDOW to the moments portiosn of the calculations (but not the WEIGHT) + * + */ # define PENALTY false @@ -106,5 +110,6 @@ cc *= weight->kernel[y][x]; } - if (window) { + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation + if (false && window) { cc *= window->kernel[y][x]; } @@ -137,5 +142,6 @@ rc *= wtVal; } - if (window) { + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation + if (false && window) { float winVal = window->kernel[y][x]; ic *= winVal; @@ -172,4 +178,5 @@ } +# define RENORM_BY_FLUX 0 // Calculate the least-squares matrix and vector for dual convolution @@ -267,32 +274,34 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { + + // XXX NOTE: clipping low S/N pixels does not seem to work very well + if (false && weight) { + float i1 = image1->kernel[y][x]; + float i2 = image2->kernel[y][x]; + float sn = (i1 + i2) / sqrt (weight->kernel[y][x]); + if (sn < 0.5) continue; + } + double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue); double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; - 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; - sumAB += ab; + + float wtVal = (weight) ? weight->kernel[y][x] : 1.0; + sumAA += wtVal*aa; + sumBB += wtVal*bb; + sumAB += wtVal*ab; if (MOMENTS) { - MxxAA += x*x*aa; - MyyAA += y*y*aa; - MxxBB += x*x*bb; - MyyBB += y*y*bb; + float winVal = (window) ? window->kernel[y][x] : 1.0; + MxxAA += winVal*x*x*aa; + MyyAA += winVal*y*y*aa; + MxxBB += winVal*x*x*bb; + MyyBB += winVal*y*y*bb; } } } + // XXX does normSquare1,2 mess up the relative scaling? + // XXX no: normSquare1,2 is the sum of the flux^2 for the source if (MOMENTS) { MxxAA /= stamp->normSquare1 * PS_SQR(normValue); @@ -304,7 +313,9 @@ // XXX this makes the Chisq portion independent of the normalization and star flux // but may be mis-scaling between stars of different fluxes +# if (RENORM_BY_FLUX) sumAA /= PS_SQR(stamp->normI2); sumAB /= PS_SQR(stamp->normI2); sumBB /= PS_SQR(stamp->normI2); +# endif // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB); @@ -327,5 +338,4 @@ matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE; matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE; - matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE; matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE; @@ -333,5 +343,4 @@ matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE; matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE; - matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE; matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE; @@ -353,5 +362,6 @@ ab *= weight->kernel[y][x]; } - if (window) { + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation + if (false && window) { ab *= window->kernel[y][x]; } @@ -362,5 +372,7 @@ // XXX this makes the Chisq portion independent of the normalization and star flux // but may be mis-scaling between stars of different fluxes +# if (RENORM_BY_FLUX) sumAB /= PS_SQR(stamp->normI2); +# endif // Spatial variation of kernel coeffs @@ -390,4 +402,10 @@ float i2 = image2->kernel[y][x]; + // XXX NOTE: clipping low S/N pixels does not seem to work very well + if (false && weight) { + float sn = (i1 + i2) / sqrt (weight->kernel[y][x]); + if (sn < 0.5) continue; + } + double ai2 = a * i2 * normValue; double bi2 = b * i2; @@ -395,28 +413,16 @@ double bi1 = b * i1 * normValue; - if (weight) { - float wtVal = weight->kernel[y][x]; - ai2 *= wtVal; - bi2 *= wtVal; - ai1 *= wtVal; - bi1 *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - ai2 *= wtVal; - bi2 *= wtVal; - ai1 *= wtVal; - bi1 *= wtVal; - } - sumAI2 += ai2; - sumBI2 += bi2; - sumAI1 += ai1; - sumBI1 += bi1; + float wtVal = (weight) ? weight->kernel[y][x] : 1.0; + sumAI2 += wtVal*ai2; + sumBI2 += wtVal*bi2; + sumAI1 += wtVal*ai1; + sumBI1 += wtVal*bi1; if (MOMENTS) { - MxxAI1 += x*x*ai1; - MyyAI1 += y*y*ai1; - MxxBI2 += x*x*bi2; - MyyBI2 += y*y*bi2; + float winVal = (window) ? window->kernel[y][x] : 1.0; + MxxAI1 += winVal*x*x*ai1; + MyyAI1 += winVal*y*y*ai1; + MxxBI2 += winVal*x*x*bi2; + MyyBI2 += winVal*y*y*bi2; } } @@ -434,8 +440,10 @@ // XXX this makes the Chisq portion independent of the normalization and star flux // but may be mis-scaling between stars of different fluxes +# if (RENORM_BY_FLUX) sumAI1 /= PS_SQR(stamp->normI2); sumBI1 /= PS_SQR(stamp->normI2); sumAI2 /= PS_SQR(stamp->normI2); sumBI2 /= PS_SQR(stamp->normI2); +# endif // Spatial variation @@ -799,11 +807,9 @@ pmSubtractionKernels *kernels = job->args->data[1]; // Kernels int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index - pmSubtractionEquationCalculationMode mode = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model - - return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode); -} - -bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, - int index, const pmSubtractionEquationCalculationMode mode) + + return pmSubtractionCalculateEquationStamp(stamps, kernels, index); +} + +bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, int index) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -832,40 +838,14 @@ psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state."); - // Generate convolutions: these are generated once and saved - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { - psError(psErrorCodeLast(), false, "Unable to convolve stamp %d.", index); - return NULL; - } - -#ifdef TESTING - for (int j = 0; j < numKernels; j++) { - if (stamp->convolutions1) { - psString convName = NULL; - psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j); - psFits *fits = psFitsOpen(convName, "w"); - psFree(convName); - psKernel *conv = stamp->convolutions1->data[j]; - psFitsWriteImage(fits, NULL, conv->image, 0, NULL); - psFitsClose(fits); - } - - if (stamp->convolutions2) { - psString convName = NULL; - psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j); - psFits *fits = psFitsOpen(convName, "w"); - psFree(convName); - psKernel *conv = stamp->convolutions2->data[j]; - psFitsWriteImage(fits, NULL, conv->image, 0, NULL); - psFitsClose(fits); - } - } -#endif - - // XXX visualize the set of convolved stamps - - psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, - stamp->xNorm, stamp->yNorm); // Polynomial terms - - bool new = stamp->vector ? false : true; // Is this a new run? + psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms + + // Is this a new run? Have we allocated the correct sized vector/matrix? + bool new = stamp->vector ? false : true; + if (!new && (stamp->vector->n != numParams)) { + psFree (stamp->vector); + psFree (stamp->matrix); + new = true; + } + if (new) { stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); @@ -940,6 +920,5 @@ } -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, - const pmSubtractionEquationCalculationMode mode) +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -952,4 +931,5 @@ for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { continue; @@ -968,10 +948,9 @@ psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32); if (!psThreadJobAddPending(job)) { return false; } } else { - pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode); + pmSubtractionCalculateEquationStamp(stamps, kernels, i); } } @@ -986,7 +965,5 @@ pmSubtractionVisualPlotLeastSquares(stamps); - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", - psTimerClear("pmSubtractionCalculateEquation")); - + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", psTimerClear("pmSubtractionCalculateEquation")); return true; @@ -997,27 +974,10 @@ bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header); -psImage *p_pmSubSolve_wUt (psVector *w, psImage *U); -psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt); - -bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask); - -bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B); -bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB); -bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB); - -bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x); -bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y); -bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); - -psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w); - -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) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + psTimerStart("pmSubtractionSolveEquation"); // Check inputs @@ -1041,4 +1001,7 @@ } + if (stamp->vector->n != numParams) { + fprintf (stderr, "mismatch length\n"); + } PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); @@ -1088,4 +1051,14 @@ #endif + // XXX TEST : print the matrix & vector + if (0) { + for (int iy = 0; iy < sumMatrix->numRows; iy++) { + for (int ix = 0; ix < sumMatrix->numCols; ix++) { + fprintf (stderr, "%e ", sumMatrix->data.F64[iy][ix]); + } + fprintf (stderr, " : %e\n", sumVector->data.F64[iy]); + } + } + psImage *invMatrix = NULL; psVector *solution = NULL; // Solution to equation! @@ -1096,31 +1069,35 @@ // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); // SINGLE solution - if (1) { - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - invMatrix = psMatrixInvert(NULL, sumMatrix, NULL); - } else { - psVector *PERM = NULL; - psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix); - solution = psMatrixLUSolution(solution, LU, sumVector, PERM); - psFree (LU); - psFree (PERM); - } - # if (1) + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + invMatrix = psMatrixInvert(NULL, sumMatrix, NULL); +# endif +# if (0) + psMatrixLUSolve(sumMatrixLU, sumVector); + solution = psMemIncrRefCounter(sumVector); + invMatrix = psMemIncrRefCounter(sumMatrix); +# endif +# if (0) + psMatrixGJSolve(sumMatrix, sumVector); + invMatrix = psMemIncrRefCounter(sumMatrix); + solution = psMemIncrRefCounter(sumVector); +# endif + +# if (0) for (int i = 0; i < solution->n; i++) { - psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i])); + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(fabs(invMatrix->data.F64[i][i]))); } # endif - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg - psVectorInit(kernels->solution1, 0.0); - kernels->solution1err = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg - psVectorInit(kernels->solution1err, 0.0); - } + // ensure we have a solution vector of the right size + kernels->solution1 = psVectorRecycle(kernels->solution1, sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution1err = psVectorRecycle(kernels->solution1err, sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + psVectorInit(kernels->solution1, 0.0); + psVectorInit(kernels->solution1err, 0.0); 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++) { kernels->solution1->data.F64[i] = solution->data.F64[i]; @@ -1180,4 +1157,14 @@ } + // XXX TEST : print the matrix & vector + if (0) { + for (int iy = 0; iy < sumMatrix->numRows; iy++) { + for (int ix = 0; ix < sumMatrix->numCols; ix++) { + fprintf (stderr, "%e ", sumMatrix->data.F64[iy][ix]); + } + fprintf (stderr, " : %e\n", sumVector->data.F64[iy]); + } + } + psImage *invMatrix = NULL; psVector *solution = NULL; // Solution to equation! @@ -1187,14 +1174,14 @@ // DUAL solution # if (1) - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-7); + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); invMatrix = psMatrixInvert(NULL, sumMatrix, NULL); # endif # if (0) psMatrixLUSolve(sumMatrix, sumVector); - solution = sumVector; - invMatrix = sumMatrix; + solution = psMemIncrRefCounter(sumVector); + invMatrix = psMemIncrRefCounter(sumMatrix); # endif -#if (1) +#if (0) for (int i = 0; i < solution->n; i++) { fprintf(stderr, "Dual solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i])); @@ -1202,16 +1189,20 @@ #endif - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); - psVectorInit (kernels->solution1, 0.0); - kernels->solution1err = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg - psVectorInit(kernels->solution1err, 0.0); - } - if (!kernels->solution2) { - kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64); - psVectorInit (kernels->solution2, 0.0); - kernels->solution2err = psVectorAlloc(numSolution2, PS_TYPE_F64); - psVectorInit(kernels->solution2err, 0.0); - } + // XXX TEST: manually set the coeffs to a desired solution + // solution->data.F64[0] = +1.826; + // solution->data.F64[1] = -0.115; + // solution->data.F64[2] = 0.0; + // solution->data.F64[3] = 0.0; + + // ensure we have solution vectors of the right size + kernels->solution1 = psVectorRecycle(kernels->solution1, numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution1err = psVectorRecycle(kernels->solution1err, numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution2 = psVectorRecycle(kernels->solution2, numSolution2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution2err = psVectorRecycle(kernels->solution2err, numSolution2, PS_TYPE_F64); // 1 for norm, 1 for bg + + psVectorInit(kernels->solution1, 0.0); + psVectorInit(kernels->solution1err, 0.0); + psVectorInit(kernels->solution2, 0.0); + psVectorInit(kernels->solution2err, 0.0); // for DUAL convolution analysis, we apply the normalization to I1 as follows: @@ -1259,4 +1250,6 @@ } } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Solve equation: %f sec", psTimerClear("pmSubtractionSolveEquation")); // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const @@ -1311,5 +1304,12 @@ // given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq -bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) { +bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) { + +# ifndef USE_WEIGHT + psAssert(weight == NULL, "impossible!"); +# endif +# ifndef USE_WINDOW + psAssert(window == NULL, "impossible!"); +# endif int npix = 0; @@ -1319,5 +1319,13 @@ for (int y = residual->yMin; y <= residual->yMax; y++) { for (int x = residual->xMin; x <= residual->xMax; x++) { - chisq += PS_SQR(residual->kernel[y][x]); + float value = PS_SQR(residual->kernel[y][x]); + if (weight) { + value *= weight->kernel[y][x]; + } + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation + if (false && window) { + value *= window->kernel[y][x]; + } + chisq += value; npix ++; } @@ -1345,5 +1353,39 @@ value2 = PS_SQR(value1); - if (weight) { + if (window) { + value1 *= window->kernel[y][x]; + value2 *= window->kernel[y][x]; + } + + fluxC1 += value1; + flux2 += value2; + fluxX += x * value2; + fluxY += y * value2; + // fluxX2 += PS_SQR(x) * value2; + // fluxY2 += PS_SQR(y) * value2; + fluxX2 += PS_SQR(x) * value1; + fluxY2 += PS_SQR(y) * value1; + } + } + // float Mx = fluxX / flux2; + // float My = fluxY / flux2; + // float Mxx = fluxX2 / flux2; + // float Myy = fluxY2 / flux2; + float Mxx = fluxX2 / fluxC1; + float Myy = fluxY2 / fluxC1; + + // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix); + moment += Mxx + Myy; + } + + // get the moments from convolved1 + if (convolved2) { + for (int y = residual->yMin; y <= residual->yMax; y++) { + for (int x = residual->xMin; x <= residual->xMax; x++) { + value1 = convolved2->kernel[y][x]; + value2 = PS_SQR(value1); + + // XXX NOTE: do NOT apply the weight to the moments calculation + if (false && weight) { value2 *= weight->kernel[y][x]; } @@ -1353,48 +1395,20 @@ } - fluxC1 += value1; + fluxC2 += value1; flux2 += value2; fluxX += x * value2; fluxY += y * value2; - fluxX2 += PS_SQR(x) * value2; - fluxY2 += PS_SQR(y) * value2; + // fluxX2 += PS_SQR(x) * value2; + // fluxY2 += PS_SQR(y) * value2; + fluxX2 += PS_SQR(x) * value1; + fluxY2 += PS_SQR(y) * value1; } } // float Mx = fluxX / flux2; // float My = fluxY / flux2; - float Mxx = fluxX2 / flux2; - float Myy = fluxY2 / flux2; - - // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix); - moment += Mxx + Myy; - } - - // get the moments from convolved1 - if (convolved2) { - for (int y = residual->yMin; y <= residual->yMax; y++) { - for (int x = residual->xMin; x <= residual->xMax; x++) { - value1 = convolved2->kernel[y][x]; - value2 = PS_SQR(value1); - - if (weight) { - value2 *= weight->kernel[y][x]; - } - if (window) { - value1 *= window->kernel[y][x]; - value2 *= window->kernel[y][x]; - } - - fluxC2 += value1; - flux2 += value2; - fluxX += x * value2; - fluxY += y * value2; - fluxX2 += PS_SQR(x) * value2; - fluxY2 += PS_SQR(y) * value2; - } - } - // float Mx = fluxX / flux2; - // float My = fluxY / flux2; - float Mxx = fluxX2 / flux2; - float Myy = fluxY2 / flux2; + // float Mxx = fluxX2 / flux2; + // float Myy = fluxY2 / flux2; + float Mxx = fluxX2 / fluxC2; + float Myy = fluxY2 / fluxC2; // fprintf (stderr, "conv2, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix); @@ -1410,13 +1424,10 @@ psVectorAppend(momentVector, moment); psVectorAppend(fluxesVector, flux); + psVectorAppend(stampMask, 0); return true; } -// typedef struct { -// float moment; -// float chisq; -// } diffStat; - -bool pmSubtractionCalculateChisqAndMoments(pmSubtractionStampList *stamps, +bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, + pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) { @@ -1425,7 +1436,11 @@ PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + psTimerStart("pmSubtractionCalculateChisqAndMoments"); + + // XXX need to save these somewhere psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); psVector *chisq = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK); int footprint = stamps->footprint; // Half-size of stamps @@ -1434,15 +1449,21 @@ psImage *polyValues = NULL; // Polynomial values + // storage for the image (convolved2 is not used in SINGLE mode) psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image - - // storage for the convolved source image (only convolved1 is used in SINGLE mode) psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + int nGood = 0; 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) { + // mark this stamp as unused (note that we have to append NANs to the other vectors to keep the lengths in sync) + psVectorAppend(moments, NAN); + psVectorAppend(fluxes, NAN); + psVectorAppend(chisq, NAN); + psVectorAppend(stampMask, 0x01); continue; } + nGood ++; // Calculate coefficients of the kernel basis functions @@ -1454,7 +1475,7 @@ psImageInit(residual->image, 0.0); - psKernel *weight = NULL; - psKernel *window = NULL; - + psKernel *weight = NULL; + psKernel *window = NULL; + #ifdef USE_WEIGHT weight = stamp->weight; @@ -1511,5 +1532,5 @@ } // XXX if we want to have a weight and window, we'll need to pass through to here - pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, NULL, residual, weight, window); + pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, NULL, residual, weight, window); } else { @@ -1555,20 +1576,55 @@ } - pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, convolved2, residual, weight, window); - } - } - - pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments); + pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, convolved2, residual, weight, window); + } + } // find the mean chisq and mean moment psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - psVectorStats (stats, chisq, NULL, NULL, 0); + psVectorStats (stats, chisq, NULL, stampMask, 0xff); float chisqValue = stats->sampleMean; psStatsInit(stats); - psVectorStats (stats, moments, NULL, NULL, 0); + psVectorStats (stats, moments, NULL, stampMask, 0xff); float momentValue = stats->sampleMean; - fprintf (stderr, "chisq: %f, moment: %f\n", chisqValue, momentValue); + // XXX this is probably wrong : I want to score against higher moments **compared with the raw momements** + // score : chisq * dMoments * modeFactor + float modeFactor = kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 2.0 : 1.0; + float orderFactor = 0.5*(kernels->spatialOrder + 1) * (kernels->spatialOrder + 2); + float score = chisqValue * momentValue * modeFactor * orderFactor; + + fprintf (stderr, "chisq: %f, moment: %f, score: %f\n", chisqValue, momentValue, score); + + // save this result if it is the first or the best (skip if bestMatch is NULL) + if (bestMatch) { + pmSubtractionQuality *match = *bestMatch; + bool keep = false; + if (match == NULL) { + *bestMatch = match = pmSubtractionQualityAlloc(); + keep = true; + } else { + if (score < match->score) { + psFree(match->fluxes); + psFree(match->chisq); + psFree(match->moments); + psFree(match->stampMask); + keep = true; + } + } + if (keep) { + fprintf (stderr, "keeping order: %d, mode: %d, score: %f\n", kernels->spatialOrder, kernels->mode, score); + match->score = score; + match->spatialOrder = kernels->spatialOrder; + match->mode = kernels->mode; + match->nGood = nGood; + match->fluxes = psMemIncrRefCounter(fluxes); + match->chisq = psMemIncrRefCounter(chisq); + match->moments = psMemIncrRefCounter(moments); + match->stampMask = psMemIncrRefCounter(stampMask); + } + } + + pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments); psFree(stats); @@ -1576,4 +1632,5 @@ psFree(fluxes); psFree(moments); + psFree(stampMask); psFree(residual); @@ -1582,13 +1639,17 @@ psFree(polyValues); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Chisq and Moments: %f sec", psTimerClear("pmSubtractionCalculateChisqAndMoments")); + return true; } -psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, - pmSubtractionKernels *kernels) +// XXX for now, let's not use this, and let's instead just use values from pmSubtractionCalculateChisqAndMoments +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + + psTimerStart("pmSubtractionCalculateDeviations"); psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps @@ -1600,7 +1661,4 @@ psImage *polyValues = NULL; // Polynomial values psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image - - // set up holding images for the visualization - pmSubtractionVisualShowFitInit (stamps); psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); @@ -1703,7 +1761,4 @@ } - // 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++) { @@ -1740,7 +1795,4 @@ } - // 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++) { @@ -1757,4 +1809,5 @@ } + double flux = 0.0; double deviation = 0.0; // Sum of differences for (int y = - footprint; y <= footprint; y++) { @@ -1762,57 +1815,17 @@ double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x]; deviation += dev; -#ifdef TESTING - residual->kernel[y][x] = dev; -#endif + flux += stamp->image1->kernel[y][x] + stamp->image2->kernel[y][x]; } } 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]); + psTrace("psModules.imcombine", 5, "Deviation and Flux for stamp %d (%d,%d): %f %f\n", + i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i], flux); psStringAppend(&log, "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)); + 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 - } @@ -1829,7 +1842,4 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background); - pmSubtractionVisualShowFit(norm); - pmSubtractionVisualPlotFit(kernels); - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psVectorStats (stats, fResSigma, NULL, NULL, 0); @@ -1862,303 +1872,7 @@ psFree(polyValues); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Deviations: %f sec", psTimerClear("pmSubtractionCalculateDeviations")); + 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 -bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) { - - psAssert (U->numRows == B->n, "U and B dimensions do not match"); - - UtB[0] = psVectorRecycle (UtB[0], 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[0]->data.F64[i] = sum; - } - return true; -} - -// w is diagonal -bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) { - - psAssert (w->n == UtB->n, "w and UtB dimensions do not match"); - - // wUt has dimensions transposed relative to Ut. - wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64); - psVectorInit (wUtB[0], 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[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i]; - } - return true; -} - -// this is basically matrix * vector -bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) { - - psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match"); - - VwUtB[0] = psVectorRecycle (*VwUtB, 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[0]->data.F64[j] = sum; - } - return true; -} - -// this is basically matrix * vector -bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) { - - psAssert (A->numCols == x->n, "A and x dimensions do not match"); - - B[0] = psVectorRecycle (*B, 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[0]->data.F64[j] = sum; - } - return true; -} - -// this is basically Vector * vector -bool p_pmSubSolve_VdV (double *value, 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]; - } - *value = sum; - return true; -} - -bool p_pmSubSolve_y2 (double *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; - } - } - } - *y2 = sum; - return true; -} - -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; -} - -bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) { - - for (int i = 0; i < w->n; i++) { - wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i]; - } - return true; -} - -// we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w) -psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) { - - psAssert (w->n == V->numCols, "w and U dimensions do not match"); - - psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64); - psImageInit (Vn, 0.0); - - // 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]; - } - } - - psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64); - psImageInit (Xvar, 0.0); - - // 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; - } - } - return Xvar; } @@ -2166,5 +1880,4 @@ // 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) { Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h (revision 30288) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h (revision 30289) @@ -4,12 +4,13 @@ #include "pmSubtractionStamps.h" #include "pmSubtractionKernels.h" +#include "pmSubtraction.h" -typedef enum { - PM_SUBTRACTION_EQUATION_NONE = 0x00, - PM_SUBTRACTION_EQUATION_NORM = 0x01, - PM_SUBTRACTION_EQUATION_BG = 0x02, - PM_SUBTRACTION_EQUATION_KERNELS = 0x04, - PM_SUBTRACTION_EQUATION_ALL = 0x07, // value should be NORM | BG | KERNELS -} pmSubtractionEquationCalculationMode; +// typedef enum { +// PM_SUBTRACTION_EQUATION_NONE = 0x00, +// PM_SUBTRACTION_EQUATION_NORM = 0x01, +// PM_SUBTRACTION_EQUATION_BG = 0x02, +// PM_SUBTRACTION_EQUATION_KERNELS = 0x04, +// PM_SUBTRACTION_EQUATION_ALL = 0x07, // value should be NORM | BG | KERNELS +// } pmSubtractionEquationCalculationMode; /// Execute a thread job to calculate the least-squares equation for a stamp @@ -20,18 +21,15 @@ bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps pmSubtractionKernels *kernels, ///< Kernel parameters - int index, ///< Index of stamp - const pmSubtractionEquationCalculationMode mode + int index ///< Index of stamp ); /// Calculate the least-squares equation to match the image quality bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps - pmSubtractionKernels *kernels, ///< Kernel parameters - const pmSubtractionEquationCalculationMode mode + pmSubtractionKernels *kernels ///< Kernel parameters ); /// Solve the least-squares equation to match the image quality bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters - const pmSubtractionStampList *stamps, ///< Stamps - const pmSubtractionEquationCalculationMode mode + const pmSubtractionStampList *stamps ///< Stamps ); @@ -92,8 +90,6 @@ bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window); -bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window); +bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window); -bool pmSubtractionCalculateChisqAndMoments(pmSubtractionStampList *stamps, - pmSubtractionKernels *kernels); - +bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, pmSubtractionStampList *stamps, pmSubtractionKernels *kernels); #endif Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c (revision 30288) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c (revision 30289) @@ -122,4 +122,5 @@ if ((image1 && image1->data.F32[y][x] < thresh1) || (image2 && image2->data.F32[y][x] < thresh2)) { + // fprintf (stderr, "%f,%f : thresh\n", xRaw, yRaw); continue; } @@ -429,4 +430,19 @@ stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, normFrac, sysErr, skyErr); + } + + // XXX TEST : dump all stars in the stamps here + if (0) { + FILE *f = fopen ("stamp.dat", "w"); + for (int i = 0; i < stamps->num; i++) { + psVector *xList = stamps->x->data[i]; + psVector *yList = stamps->y->data[i]; // Coordinate lists + psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes + + for (int j = 0; j < xList->n; j++) { + fprintf (f, "%d %d %f %f %f\n", i, j, xList->data.F32[j], yList->data.F32[j], fluxList->data.F32[j]); + } + } + fclose (f); } @@ -636,4 +652,6 @@ } + int nTotal = 0; + // Sort the list by flux, with the brightest last for (int i = 0; i < numStamps; i++) { @@ -662,6 +680,8 @@ stamps->y->data[i] = ySorted; stamps->flux->data[i] = fluxSorted; - } - + nTotal += num; + } + // fprintf (stderr, "nTotal %d\n", nTotal); + return stamps; } @@ -809,11 +829,11 @@ float R2 = Sr2 / Sf2; - stamps->normWindow1 = 2.0*R1; - stamps->normWindow2 = 2.0*R2; + stamps->normWindow1 = 2.75*R1; + stamps->normWindow2 = 2.75*R2; psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2); // Generate a weighting window based on the kron radii { - float radius = 1.5 * PS_MAX(R1, R2); + float radius = 2.0 * PS_MAX(R1, R2); psImageInit(stamps->window->image, 0.0); @@ -826,4 +846,8 @@ } } + + // complain if normWindow1 or normWindow2 are larger than size + psAssert(stamps->normWindow1 < size, "norm Window 1 too large"); + psAssert(stamps->normWindow2 < size, "norm Window 2 too large"); # else @@ -1142,4 +1166,8 @@ continue; } + + // XXX TEST + if (source->errMag > 0.1) continue; + if (source->modelPSF) { x->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_XPOS]; Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c (revision 30288) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c (revision 30289) @@ -20,4 +20,5 @@ #include "pmSubtractionEquation.h" #include "pmSubtractionKernels.h" +#include "pmSubtractionVisual.h" #include "pmVisual.h" @@ -423,5 +424,6 @@ } - // XXX clear the overlay(s) (red at least!) + // clear the overlay (red at least!) + KiiEraseOverlay (kapa2, "red"); // get the kernel sizes @@ -456,4 +458,5 @@ int nKernels = 0; + // paste in the kernel images, scaled by sum2 if (maxStamp->convolutions1) { // output image is a grid of NXsub by NYsub sub-images @@ -478,14 +481,16 @@ int yPix = ySub * (2*footprint + 1 + 3) + footprint; - double sum = 0.0; double sum2 = 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]; sum2 += PS_SQR(kernel->kernel[y][x]); } } - // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); + float scale = sqrt(sum2) / PS_SQR(2*footprint + 1); + 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] / scale; + } + } } pmVisualScaleImage(kapa2, output, "Image", 0, true); @@ -520,14 +525,16 @@ int yPix = ySub * (2*footprint + 1 + 3) + footprint; - double sum = 0.0; double sum2 = 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]; sum2 += PS_SQR(kernel->kernel[y][x]); } } - // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); + float scale = sqrt(sum2) / PS_SQR(2*footprint + 1); + 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] / scale; + } + } } pmVisualScaleImage(kapa2, output, "Image", 1, true); @@ -587,122 +594,4 @@ KiiLoadOverlay (kapa2, overlay, Noverlay, "red"); FREE (overlay); - return true; -} - -static int footprint = 0; -static int NX = 0; -static int NY = 0; -static psImage *sourceImage = NULL; -static psImage *targetImage = NULL; -static psImage *residualImage = NULL; -static psImage *fresidualImage = NULL; -static psImage *differenceImage = NULL; -static psImage *convolutionImage = NULL; - -bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) { - - if (!pmVisualTestLevel("ppsub.fit", 1)) return true; - - // generate 4 storage images large enough to hold the N stamps: - - footprint = stamps->footprint; - - 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; - - int NXpix = (2*footprint + 1) * NX; - NXpix += (NX > 1) ? 3 * NX : 0; - - int NYpix = (2*footprint + 1) * NY; - NYpix += (NY > 1) ? 3 * NY : 0; - - sourceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - targetImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - residualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - fresidualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - differenceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - - psImageInit (sourceImage, 0.0); - psImageInit (targetImage, 0.0); - psImageInit (residualImage, 0.0); - psImageInit (fresidualImage, 0.0); - psImageInit (differenceImage, 0.0); - psImageInit (convolutionImage, 0.0); - - return true; -} - -bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) { - - if (!pmVisualTestLevel("ppsub.stamp", 1)) return true; - - double sum; - - int NXoff = index % NX; - int NYoff = index / NX; - - int NXpix = NXoff * (2*footprint + 1 + 3) + footprint; - int NYpix = NYoff * (2*footprint + 1 + 3) + footprint; - - // 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] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); - } - } return true; } @@ -730,15 +619,15 @@ } -bool pmSubtractionVisualShowFit(double norm) { - - // for testing, dump the residual image and exit - if (0) { - 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); - } +static int footprint = 0; +static int NX = 0; +static int NY = 0; +static psImage *sourceImage = NULL; +static psImage *targetImage = NULL; +static psImage *residualImage = NULL; +static psImage *fresidualImage = NULL; +static psImage *differenceImage = NULL; +static psImage *convolutionImage = NULL; + +bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) { if (!pmVisualTestLevel("ppsub.fit", 1)) return true; @@ -746,4 +635,187 @@ if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false; if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false; + + // set up holding images for the visualization + pmSubtractionVisualShowFitInit (stamps); + + int numKernels = kernels->num; // Number of kernels + + psImage *polyValues = NULL; // Polynomial values + psKernel *residual = psKernelAlloc(-stamps->footprint, stamps->footprint, -stamps->footprint, stamps->footprint); // Residual image + + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + + 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) { continue; } + + // Calculate coefficients of the kernel basis functions + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); + double background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Difference in background + + 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; + 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; + } + } + } + // visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i); + } 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; + } + } + } + // visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i); + } + } + pmSubtractionVisualShowFitImage(norm); + + return true; +} + +// generate 4 storage images large enough to hold the stamps: +bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) { + + footprint = stamps->footprint; + + 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; + + int NXpix = (2*footprint + 1) * NX; + NXpix += (NX > 1) ? 3 * NX : 0; + + int NYpix = (2*footprint + 1) * NY; + NYpix += (NY > 1) ? 3 * NY : 0; + + sourceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + targetImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + residualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + fresidualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + differenceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + + psImageInit (sourceImage, 0.0); + psImageInit (targetImage, 0.0); + psImageInit (residualImage, 0.0); + psImageInit (fresidualImage, 0.0); + psImageInit (differenceImage, 0.0); + psImageInit (convolutionImage, 0.0); + + return true; +} + +bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) { + + double sum; + + int NXoff = index % NX; + int NYoff = index / NX; + + int NXpix = NXoff * (2*footprint + 1 + 3) + footprint; + int NYpix = NYoff * (2*footprint + 1 + 3) + footprint; + + // 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] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); + } + } + return true; +} + +bool pmSubtractionVisualShowFitImage(double norm) { KiiEraseOverlay (kapa1, "red"); @@ -792,4 +864,5 @@ psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32); psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32); + psVector *dy = psVectorAllocEmpty (kernels->num, PS_TYPE_F32); graphdata.xmin = -1.0; @@ -804,8 +877,9 @@ x->data.F32[i] = i; y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false); + dy->data.F32[i] = kernels->solution1err->data.F64[i]; graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]); graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]); } - x->n = y->n = kernels->num; + x->n = y->n = dy->n = kernels->num; float range; @@ -828,11 +902,15 @@ graphdata.size = 0.5; graphdata.style = 2; + graphdata.etype |= 0x01; KapaPrepPlot (kapa3, x->n, &graphdata); KapaPlotVector (kapa3, x->n, x->data.F32, "x"); KapaPlotVector (kapa3, x->n, y->data.F32, "y"); + KapaPlotVector (kapa3, x->n, dy->data.F32, "dym"); + KapaPlotVector (kapa3, x->n, dy->data.F32, "dyp"); psFree (x); psFree (y); + psFree (dy); psFree (polyValues); Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.h (revision 30288) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.h (revision 30289) @@ -8,7 +8,10 @@ bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed); bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub); + +bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels); bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps); bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index); -bool pmSubtractionVisualShowFit(double norm); +bool pmSubtractionVisualShowFitImage(double norm); + bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels); bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels);