Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26491) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26502) @@ -15,8 +15,8 @@ #include "pmSubtractionVisual.h" -//#define TESTING // TESTING output for debugging; may not work with threads! +// #define TESTING // TESTING output for debugging; may not work with threads! #define USE_WEIGHT // Include weight (1/variance) in equation? -#define USE_WINDOW // Include weight (1/variance) in equation? +// #define USE_WINDOW // Include weight (1/variance) in equation? @@ -157,5 +157,7 @@ if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { vector->data.F64[iIndex] = sumIC * poly[iTerm]; - if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { + // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B), + // instead, calculate A - \sum(k x B), with full hermitians + if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) { // subtract norm * sumRC * poly[iTerm] psAssert (kernels->solution1, "programming error: define solution first!"); @@ -228,4 +230,5 @@ const psKernel *image2, // Image 2 const psKernel *weight, // Weight image + const psKernel *window, // Window image const psArray *convolutions1, // Convolutions of image 1 for each kernel const psArray *convolutions2, // Convolutions of image 2 for each kernel @@ -291,4 +294,10 @@ ab *= wtVal; } + if (window) { + float wtVal = window->kernel[y][x]; + aa *= wtVal; + bb *= wtVal; + ab *= wtVal; + } sumAA += aa; sumBB += bb; @@ -319,4 +328,7 @@ if (weight) { ab *= weight->kernel[y][x]; + } + if (window) { + ab *= window->kernel[y][x]; } sumAB += ab; @@ -365,4 +377,15 @@ i2 *= wtVal; } + if (window) { + float wtVal = window->kernel[y][x]; + ai2 *= wtVal; + bi2 *= wtVal; + ai1 *= wtVal; + bi1 *= wtVal; + i1i2 *= wtVal; + a *= wtVal; + b *= wtVal; + i2 *= wtVal; + } sumAI2 += ai2; sumBI2 += bi2; @@ -411,4 +434,12 @@ if (weight) { float wtVal = weight->kernel[y][x]; + i1 *= wtVal; + i1i1 *= wtVal; + one *= wtVal; + i2 *= wtVal; + i1i2 *= wtVal; + } + if (window) { + float wtVal = window->kernel[y][x]; i1 *= wtVal; i1i1 *= wtVal; @@ -516,4 +547,7 @@ for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { // Contribution to chi^2: a_i^2 P_i + if (!isfinite(penalties->data.F32[i])) { + psAbort ("invalid penalty"); + } matrix->data.F64[index][index] -= norm * penalties->data.F32[i]; } @@ -702,5 +736,4 @@ break; case PM_SUBTRACTION_MODE_DUAL: - psAbort ("dual is disabled: need to add window and calculation mode"); if (new) { stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); @@ -714,5 +747,5 @@ #endif status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, - stamp->matrixX, stamp->image1, stamp->image2, weight, + stamp->matrixX, stamp->image1, stamp->image2, weight, window, stamp->convolutions1, stamp->convolutions2, kernels, polyValues, footprint); @@ -1587,4 +1620,8 @@ } } + + // 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++) { Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26491) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26502) @@ -132,4 +132,7 @@ kernels->preCalc->data[index] = NULL; kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); + if (!isfinite(kernels->penalties->data.F32[index])) { + psAbort ("invalid penalty"); + } psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); @@ -146,5 +149,5 @@ } -bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm) { +bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm, bool AlardLuptonStyle) { // Calculate moments @@ -157,33 +160,67 @@ } - // Normalize sum of kernel component to unity for even functions - if (uOrder % 2 == 0 && vOrder % 2 == 0) { - double sum = 0.0; // Sum of kernel component - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - sum += preCalc->kernel->kernel[v][u]; - } - } - sum = 1.0 / sqrt(sum); - psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); - psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); - psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); - -#if 1 - fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum, preCalc->kernel->kernel[0][0]); -#endif - - preCalc->kernel->kernel[0][0] -= 1.0; - moment *= PS_SQR(sum); - } - -#if 1 + // we have 4 cases here: + // 1) for odd functions, normalize the kernel by the maximum swing / Npix + // 2) for even functions, normalize the kernel to unity + // 3) for alard-lupton style normalization, subtract 1 from the 0,0 pixel for all even functions + // 4) for deconvolved hermitians, subtract 1 from the 0,0 pixel for the 0,0 function(s) + double sum = 0.0; // Sum of kernel component + double min = FLT_MAX; + double max = FLT_MIN; + for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { sum += preCalc->kernel->kernel[v][u]; + min = PS_MIN(preCalc->kernel->kernel[v][u], min); + max = PS_MAX(preCalc->kernel->kernel[v][u], max); } } - fprintf(stderr, "%d sum: %lf\n", index, sum); +#if 1 + fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, moment); +#endif + + // only even terms have non-zero sums + if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) { + moment /= sum; + } else { + moment = 0.0; + } + + bool zeroNull = false; + float scale1D = 1.0 / sqrt(sum); + float scale2D = 1.0 / sum; + + if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) { + zeroNull = true; + } + if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) { + zeroNull = true; + } + if ((uOrder % 2) || (vOrder % 2)) { + // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max); + scale2D = 1.0 / max; + scale1D = sqrt(scale2D); + } + + psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32)); + psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32)); + psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32)); + if (zeroNull) { + preCalc->kernel->kernel[0][0] -= 1.0; + } + +#if 1 + sum = 0.0; // Sum of kernel component + min = FLT_MAX; + max = FLT_MIN; + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += preCalc->kernel->kernel[v][u]; + min = PS_MIN(preCalc->kernel->kernel[v][u], min); + max = PS_MAX(preCalc->kernel->kernel[v][u], max); + } + } + fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D); #endif @@ -196,4 +233,7 @@ kernels->preCalc->data[index] = preCalc; kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); + if (!isfinite(kernels->penalties->data.F32[index])) { + psAbort ("invalid penalty"); + } psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, @@ -250,5 +290,5 @@ pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); } } @@ -302,5 +342,5 @@ for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); } } @@ -338,5 +378,5 @@ int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); - num += (gaussOrder + 1) * (gaussOrder + 2) / 2; + num += PS_SQR(gaussOrder + 1); } @@ -349,6 +389,10 @@ // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest) // generate the Gaussian deconvolution kernel - # define DECONV_SIGMA 1.7 + # define DECONV_SIGMA 1.6 psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECONV_SIGMA); + +# if 1 + psArray *deconKernels = psArrayAllocEmpty(100); +# endif // Set the kernel parameters @@ -357,5 +401,5 @@ // Iterate over (u,v) order for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { - for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { + for (int vOrder = 0; vOrder <= orders->data.S32[i]; vOrder++, index++) { pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values @@ -365,12 +409,41 @@ preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false); + // XXXX test demo that deconvolved kernel is valid - // psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); - // pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); - - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]); +# if 1 + psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); + psArrayAdd (deconKernels, 100, kernelConv); + psFree (kernelConv); + + if (!uOrder && !vOrder){ + pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); + } +# endif } } } + +# if 1 + psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32); + for (int i = 0; i < deconKernels->n; i++) { + for (int j = 0; j <= i; j++) { + psImage *t1 = deconKernels->data[i]; + psImage *t2 = deconKernels->data[j]; + + double sum = 0.0; + for (int iy = 0; iy < t1->numRows; iy++) { + for (int ix = 0; ix < t1->numCols; ix++) { + sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix]; + } + } + dot->data.F32[j][i] = sum; + dot->data.F32[i][j] = sum; + } + } + pmSubtractionVisualShowSubtraction (dot, NULL, NULL); + psFree (dot); + psFree (deconKernels); +# endif return kernels; @@ -857,4 +930,7 @@ kernels->v->data.S32[index] = vOrder; kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); + if (!isfinite(kernels->penalties->data.F32[index])) { + psAbort ("invalid penalty"); + } psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26491) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26502) @@ -81,5 +81,5 @@ // Search a region for a suitable stamp -bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return +bool stampSearch(float *xStamp, float *yStamp, // Coordinates of stamp, to return float *fluxStamp, // Flux of stamp, to return const psImage *image1, const psImage *image2, // Images to search @@ -116,6 +116,6 @@ ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel if (flux > *fluxStamp) { - *xStamp = x; - *yStamp = y; + *xStamp = x + 0.5; + *yStamp = y + 0.5; *fluxStamp = flux; found = true; @@ -405,5 +405,5 @@ numSearch++; - int xStamp = 0, yStamp = 0; // Coordinates of stamp + float xStamp = 0.0, yStamp = 0.0; // Coordinates of stamp float fluxStamp = -INFINITY; // Flux of stamp bool goodStamp = false; // Found a good stamp? @@ -437,4 +437,9 @@ goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, subMask, xMin, xMax, yMin, yMax, numCols, numRows, border); + + // XXX reset for a test: + xStamp = xList->data.F32[j]; + yStamp = yList->data.F32[j]; + fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp); } } else { @@ -543,4 +548,6 @@ } + fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix); + bool found = false; for (int j = 0; j < numStamps && !found; j++) { @@ -735,4 +742,5 @@ return false; } + fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y); // Catch memory leaks --- these should have been freed and NULLed before @@ -811,4 +819,5 @@ y->data.F32[numOut] = source->peak->yf; } + fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]); numOut++; } Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26491) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26502) @@ -251,9 +251,12 @@ int yPix = ySub * (2*footprint + 1 + 3) + footprint; + double sum = 0.0; for (int y = -footprint; y <= footprint; y++) { for (int x = -footprint; x <= footprint; x++) { output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; + sum += kernel->kernel[y][x]; } } + fprintf (stderr, "kernel %d, sum %f\n", i, sum); } @@ -274,12 +277,23 @@ // choose the brightest stamp - pmSubtractionStamp *maxStamp = stamps->stamps->data[0]; - float maxFlux = maxStamp->flux; - for (int i = 1; i < stamps->num; i++) { + pmSubtractionStamp *maxStamp = NULL; + float maxFlux = NAN; + for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (!isfinite(stamp->flux)) continue; + if (!stamp->convolutions1 && !stamp->convolutions2) continue; + if (!maxStamp) { + maxFlux = stamp->flux; + maxStamp = stamp; + continue; + } if (stamp->flux > maxFlux) { maxFlux = stamp->flux; maxStamp = stamp; } + } + + if (!isfinite(maxStamp->flux)) { + fprintf (stderr, "no valid stamps?\n"); } @@ -308,9 +322,12 @@ int yPix = ySub * (2*footprint + 1 + 3) + footprint; + double sum = 0.0; for (int y = -footprint; y <= footprint; y++) { for (int x = -footprint; x <= footprint; x++) { output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; + sum += kernel->kernel[y][x]; } } + fprintf (stderr, "kernel %d, sum %f\n", i, sum); } pmVisualScaleImage(kapa2, output, "Image", 0, true); @@ -339,9 +356,12 @@ int yPix = ySub * (2*footprint + 1 + 3) + footprint; + double sum = 0.0; for (int y = -footprint; y <= footprint; y++) { for (int x = -footprint; x <= footprint; x++) { output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; + sum += kernel->kernel[y][x]; } } + fprintf (stderr, "kernel %d, sum %f\n", i, sum); } pmVisualScaleImage(kapa2, output, "Image", 1, true); @@ -375,9 +395,9 @@ overlay[Noverlay].type = KII_OVERLAY_BOX; if ((stamp->x < 1.0) && (stamp->y < 1.0)) { - fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y); + // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y); continue; } if (!isfinite(stamp->x) && !isfinite(stamp->y)) { - fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y); + // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y); continue; } @@ -409,5 +429,5 @@ bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) { - if (!pmVisualIsVisual()) return true; + // if (!pmVisualIsVisual()) return true; // generate 4 storage images large enough to hold the N stamps: @@ -446,5 +466,5 @@ bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) { - if (!pmVisualIsVisual()) return true; + // if (!pmVisualIsVisual()) return true; double sum; @@ -517,5 +537,5 @@ bool pmSubtractionVisualShowFit() { - if (!pmVisualIsVisual()) return true; + // if (!pmVisualIsVisual()) return true; if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false; if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;