Index: trunk/psModules/src/imcombine/README =================================================================== --- trunk/psModules/src/imcombine/README (revision 29543) +++ trunk/psModules/src/imcombine/README (revision 29543) @@ -0,0 +1,9 @@ + +pmSubtractionEquation.v0.c + + This version of the code includes the old-style form of the matrix + equations, in which the normalization and background were calculated + as part of the matrix inversion. We've decided that this method is + too sensitive to errors in the normalization and no longer use this + approach. + Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 29391) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 29543) @@ -380,6 +380,5 @@ // measure the source moments: tophat windowing, no pixel S/N cutoff - // XXX probably should be passing the maskVal to this function so we can pass it along here... - if (!pmSourceMoments(source, maxRadius, 0.0, 1.0, maskVal)) { + if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, maskVal)) { // Can't do anything about it; limp along as best we can psErrorClear(); Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 29391) +++ trunk/psModules/src/imcombine/pmStack.c (revision 29543) @@ -19,4 +19,5 @@ #include +#include // for memset #include @@ -34,9 +35,18 @@ -//#define TESTING // Enable test output -//#define TEST_X 843-1 // x coordinate to examine -//#define TEST_Y 813-1 // y coordinate to examine -//#define TEST_RADIUS 0 // Radius to examine - +# if (0) +#define TESTING // Enable test output +#define TEST_X 4963 // x coordinate to examine +#define TEST_Y 5353 // y coordinate to examine +#define TEST_X 40 // x coordinate to examine +#define TEST_Y 40 // y coordinate to examine +#define TEST_RADIUS 0.5 // Radius to examine +# endif + +# ifdef TESTING +# define CHECKPIX(XPIX,YPIX,MSG,...) { if (PS_SQR(XPIX - TEST_X) + PS_SQR(YPIX - TEST_Y) <= PS_SQR(TEST_RADIUS)) { fprintf(stderr,MSG,__VA_ARGS__); } } +# else +# define CHECKPIX(XPIX,YPIX,MSG,...) { } +# endif // Data structure for use as a buffer when combining pixels @@ -250,9 +260,5 @@ ) { -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Marking image %d, pixel %d,%d for inspection\n", source, x, y); - } -#endif + CHECKPIX(x, y, "Marking image %d, pixel %d,%d for inspection\n", source, x, y); pmStackData *data = inputs->data[source]; // Stack data of interest if (!data) { @@ -272,9 +278,5 @@ ) { -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Marking pixel image %d, pixel %d,%d for rejection\n", source, x, y); - } -#endif + CHECKPIX(x, y, "Marking pixel image %d, pixel %d,%d for rejection\n", source, x, y); pmStackData *data = inputs->data[source]; // Stack data of interest if (!data) { @@ -290,4 +292,6 @@ static void combineExtract(int *num, // Number of good pixels bool *suspect, // Any suspect pixels? + psImageMaskType *badMask, // OR of all bad (masked) pixels + psImageMaskType *goodMask, // OR of all good (unmasked) pixels combineBuffer *buffer, // Buffer with vectors psImage *image, // Combined image, for output @@ -300,6 +304,6 @@ const psVector *reject, // Indices of pixels to reject, or NULL int x, int y, // Coordinates of interest; frame of output image - psImageMaskType maskVal, // Value to mask - psImageMaskType maskSuspect // Value to suspect + psImageMaskType badMaskBits, // Value to mask as 'bad' + psImageMaskType suspectMaskBits // Value to mask as 'suspect' ) { @@ -322,4 +326,12 @@ } + // mask values to store possible mask combinations + *badMask = 0; + *goodMask = 0xffff; + + int nGoodBits[16]; // accumulate the good pixel bits here for fuzzy logic + psAssert (sizeof(psImageMaskType) == 2, "psImageMaskType is not the expected size"); + memset (nGoodBits, 0, 16*sizeof(int)); + // Extract the pixel and mask data int numGood = 0; // Number of good pixels @@ -328,5 +340,22 @@ // should be because of how pixelMapGenerate works if (reject && reject->data.U16[j] == i) { + // pixels can be rejected because: + // 1) only 1 input pixel (and 'safe' is true) + // 2) only 2 input pixels were available and they were mutually inconsistent (or variance info was missing) + // 3) NOTE : ifdef'ed out code for 3 inputs case + // 4) outlier from sample of N pixels + // 5) some of these may have been suspect. + // XXX raise a specific mask bit for these (currently results in BLANK) j++; + + pmStackData *data = inputs->data[i]; // Stack data of interest + if (data) { + int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout + psImage *mask = data->readout->mask; // Mask of interest + *badMask |= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the bad bits used + CHECKPIX(x, y, "reject: adding bit to mask: %d : %x (badMask = %x)\n", i, mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn], *badMask); + } else { + CHECKPIX(x, y, "reject: no item in data (badMask = %x)\n", *badMask); + } continue; } @@ -334,4 +363,5 @@ pmStackData *data = inputs->data[i]; // Stack data of interest if (!data) { + CHECKPIX(x, y, "skip: no item in data (badMask = %x)\n", *badMask); continue; } @@ -339,10 +369,23 @@ int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout psImage *mask = data->readout->mask; // Mask of interest - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits) { + *badMask |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits); // save the bad bits used + CHECKPIX(x, y, "skip: adding bit to mask: %d : %x (badMask = %x)\n", i, mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn], *badMask); continue; } - pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ? - true : false; + pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & suspectMaskBits ? true : false; + + // *goodMask &= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the mask bits still used + // check for set bits and increment counter as appropriate + { + psImageMaskType value = 0x0001; + for (int nbit = 0; nbit < 16; nbit ++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) { + nGoodBits[nbit] ++; + } + value <<= 1; + } + } psImage *image = data->readout->image; // Image of interest @@ -356,4 +399,6 @@ pixelSources->data.U16[numGood] = i; numGood++; + + CHECKPIX(x, y, "keep: %d : %x (badMask = %x)\n", i, mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn], *badMask); } pixelData->n = numGood; @@ -367,11 +412,24 @@ *num = numGood; + // set the mask bits if nGoodBits[i] > f*numGood + { + # define SUSPECT_FRACTION 0.65 + *goodMask = 0x0000; + psImageMaskType value = 0x0001; + for (int nbit = 0; nbit < 16; nbit ++) { + if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) { + *goodMask |= value; + } + value <<= 1; + } + } + #ifdef TESTING if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { for (int i = 0; i < numGood; i++) { - fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d\n", + fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n", i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i], - pixelSuspects->data.U8[i]); + pixelSuspects->data.U8[i], badMaskBits, suspectMaskBits, *badMask, *goodMask); } } @@ -380,5 +438,4 @@ return; } - // Combine pixels @@ -392,5 +449,7 @@ combineBuffer *buffer, // Buffer with vectors int x, int y, // Coordinates of interest; frame of output image - psImageMaskType bad, // Value for bad pixels + psImageMaskType blankMask, // Value for empty pixels + psImageMaskType badMask, // Value for bad pixels + psImageMaskType goodMask, // Value for good pixels bool safe, // Safe combination? float invTotalWeight // Inverse of total weight for all inputs @@ -404,15 +463,17 @@ // Default option is that the pixel is bad float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map - psImageMaskType maskValue = bad; // Value for combined mask float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted) + + // default output mask value. badMask is the OR of all unused input pixels. + // if there are no input pixels, this value will be 0, in which case we want to set the output pixel to BLANK. + // if there are only good input pixels, and they do not result in a valid pixel, we still want to set this to BLANK. + psImageMaskType maskValue = badMask ? badMask : blankMask; // Value for combined mask + + CHECKPIX(x, y, "bad vs good : %x %x %x\n", maskValue, badMask, blankMask); switch (num) { case 0: { // Nothing to combine: it's bad -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "No inputs to combine, pixel %d,%d is bad.\n", x, y); - } -#endif + CHECKPIX(x, y, "No inputs to combine, pixel %d,%d is bad.\n", x, y); break; } @@ -428,17 +489,9 @@ expWeightValue = pixelExps->data.F32[0]; } - maskValue = 0; -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Single input to combine, safety off, pixel %d,%d --> %f\n", - x, y, imageValue); - } -#endif - } -#ifdef TESTING - else if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Single input to combine, safety on, pixel %d,%d is bad.\n", x, y); - } -#endif + maskValue = goodMask; + CHECKPIX(x, y, "Single input to combine, safety off, pixel %d,%d --> %f\n", x, y, imageValue); + } else { + CHECKPIX(x, y, "Single input to combine, safety on, pixel %d,%d is bad.\n", x, y); + } break; } @@ -446,22 +499,11 @@ // Automatically accept the mean of the pixels only if we're not playing safe if (!safe) { - if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue, - pixelData, pixelVariances, pixelExps, pixelWeights)) { - maskValue = 0; -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n", - x, y, imageValue, varianceValue); - } -#endif + if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue, pixelData, pixelVariances, pixelExps, pixelWeights)) { + maskValue = goodMask; + CHECKPIX(x, y, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue); } + } else { + CHECKPIX(x, y, "Two inputs to combine, safety on, pixel %d,%d is bad\n", x, y); } -#ifdef TESTING - else { - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Two inputs to combine, safety on, pixel %d,%d is bad\n", x, y); - } - } -#endif break; } @@ -472,10 +514,6 @@ break; } - maskValue = 0; -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue); - } -#endif + maskValue = goodMask; + CHECKPIX(x, y, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue); break; } @@ -522,10 +560,5 @@ int numIter = PS_MAX(iter * num, 1); // Number of iterations -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Testing pixel %d,%d: %d %f %f %f %d %d\n", - x, y, numIter, rej, sys, olympic, useVariance, safe); - } -#endif + CHECKPIX(x, y, "Testing pixel %d,%d: %d %f %f %f %d %d\n", x, y, numIter, rej, sys, olympic, useVariance, safe); psVector *pixelData = buffer->pixels; // Values for the pixel of interest @@ -548,12 +581,7 @@ // Systematic error contributes to the rejection level float sysVar = PS_SQR(sys * pixelData->data.F32[i]); -#ifdef TESTING - // Correct variance for comparison against weighted mean including itself - float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights; - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Variance %d (%d), pixel %d,%d: %f %f %f\n", i, pixelSources->data.U16[i], - x, y, pixelVariances->data.F32[i], sysVar, compare); - } -#endif + CHECKPIX(x, y, "Variance %d (%d), pixel %d,%d: %f %f %f\n", + i, pixelSources->data.U16[i], x, y, + pixelVariances->data.F32[i], sysVar, 1.0 - pixelWeights->data.F32[i] / sumWeights); pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar); } @@ -593,10 +621,5 @@ combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); } -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Flagged both inputs for pixel %d,%d (%f > %f x %f\n)", - x, y, diff, rej, sqrtf(sigma2)); - } -#endif + CHECKPIX(x, y, "Flagged both inputs for pixel %d,%d (%f > %f x %f\n)", x, y, diff, rej, sqrtf(sigma2)); } } else if (i == 0 && safe) { @@ -708,17 +731,9 @@ float median = combinationWeightedOlympic(pixelData, pixelWeights, olympic, buffer->sort); // Median for stack -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Flag with variance pixel %d,%d: median = %f\n", x, y, median); - } -#endif + CHECKPIX(x, y, "Flag with variance pixel %d,%d: median = %f\n", x, y, median); float worst = -INFINITY; // Largest deviation for (int j = 0; j < num; j++) { float diff = pixelData->data.F32[j] - median; // Difference from expected -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff); - } -#endif + CHECKPIX(x, y, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff); // Comparing squares --- cheaper than lots of sqrts @@ -737,11 +752,5 @@ combinationMedianStdev(&median, &stdev, pixelData, buffer->sort); float limit = rej * stdev; // Rejection limit -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, - "Flag without variance pixel %d,%d; median = %f, stdev = %f, limit = %f\n", - x, y, median, stdev, limit); - } -#endif + CHECKPIX(x, y, "Flag without variance pixel %d,%d; median = %f, stdev = %f, limit = %f\n", x, y, median, stdev, limit); float worst = -INFINITY; // Largest deviation for (int j = 0; j < num; j++) { @@ -763,9 +772,5 @@ if (maskIndex >= 0) { if (suspect) { -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Throwing out all suspect pixels for %d,%d\n", x, y); - } -#endif + CHECKPIX(x, y, "Throwing out all suspect pixels for %d,%d\n", x, y); // Throw out all suspect pixels int numGood = 0; // Number of good pixels @@ -796,9 +801,5 @@ } else { // Throw out masked pixel -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Throwing out input %d for pixel %d,%d\n", maskIndex, x, y); - } -#endif + CHECKPIX(x, y, "Throwing out input %d for pixel %d,%d\n", maskIndex, x, y); combineMarkInspect(inputs, x, y, pixelSources->data.U16[maskIndex]); int numGood = 0; // Number of good pixels @@ -999,9 +1000,19 @@ /// Stack input images -bool pmStackCombine(pmReadout *combined, pmReadout *expmaps, psArray *input, - psImageMaskType maskVal, psImageMaskType maskSuspect, - psImageMaskType bad, int kernelSize, - float iter, float rej, float sys, float olympic, - bool useVariance, bool safe, bool rejection) +bool pmStackCombine( + pmReadout *combined, // output stacked readout + pmReadout *expmaps, // output exposure map information + psArray *input, // input exposures + psImageMaskType badMaskBits, // treat these bits as 'bad' + psImageMaskType suspectMaskBits, // treat these bits as 'suspect' + psImageMaskType blankMaskBits, // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?) + int kernelSize, + float iter, + float rej, + float sys, + float olympic, + bool useVariance, + bool safe, + bool rejection) { bool haveVariances; // Do we have the variance maps? @@ -1012,5 +1023,5 @@ } PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - PS_ASSERT_INT_POSITIVE(bad, false); + PS_ASSERT_INT_POSITIVE(blankMaskBits, false); if (isnan(rej)) { PS_ASSERT_FLOAT_EQUAL(iter, 0, false); @@ -1070,5 +1081,5 @@ } psFree(stack); - pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, bad); + pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, blankMaskBits); psTrace("psModules.imcombine", 1, "Have for combination [%d:%d,%d:%d] (%dx%d)\n", minInputCols, maxInputCols, minInputRows, maxInputRows, xSize, ySize); @@ -1131,10 +1142,12 @@ for (int y = minInputRows; y < maxInputRows; y++) { for (int x = minInputCols; x < maxInputCols; x++) { + CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection); + #ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf(stderr, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", - x, y, maskVal, bad, iter, rej, sys, olympic, useVariance, safe, rejection); - } -#endif + if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { + fprintf (stderr, "combine\n"); + } +# endif + psVector *reject = NULL; // Images to reject for this pixel if (rejection) { @@ -1154,11 +1167,11 @@ #endif } - - int num; // Number of good pixels - bool suspect; // Suspect pixels in stack? - combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance, - input, weights, exps, addVariance, reject, x, y, maskVal, maskSuspect); - combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight, - num, buffer, x, y, bad, safe, totalExpWeight); + + int num; // Number of good pixels + bool suspect; // Suspect pixels in stack? + psImageMaskType badMask = 0; // OR of mask bits in all bad input pixels + psImageMaskType goodMask = 0; // OR of mask bits in all good input pixels + combineExtract(&num, &suspect, &badMask, &goodMask, buffer, combinedImage, combinedMask, combinedVariance, input, weights, exps, addVariance, reject, x, y, badMaskBits, suspectMaskBits); + combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight, num, buffer, x, y, blankMaskBits, badMask, goodMask, safe, totalExpWeight); if (iter > 0) { Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 29543) @@ -540,4 +540,8 @@ psImage *conv = psImageConvolveFFT(NULL, image->image, NULL, 0, kernel); // Convolved image + + // note: do not attempt to renormalize kernels here: cannot have different stars with + // different kernel ratios + int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 29543) @@ -299,17 +299,17 @@ kernels->rms); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, 0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, 0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, 0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); + + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, 0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, 0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, 0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); } Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 29543) @@ -24,10 +24,10 @@ #define PM_SUBTRACTION_ANALYSIS_DECONV_MAX "SUBTRACTION.DECONV.MAX" // Maximum deconvolution fraction -#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN "SUBTRACTION.RES.FSIGMA.MEAN" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV "SUBTRACTION.RES.FSIGMA.STDEV" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN "SUBTRACTION.RES.FMIN.MEAN" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV "SUBTRACTION.RES.FMIN.STDEV" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN "SUBTRACTION.RES.FMAX.MEAN" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV "SUBTRACTION.RES.FMAX.STDEV" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN "SUBTRACTION.FRES.MEAN" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV "SUBTRACTION.FRES.STDEV" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN "SUBTRACTION.FRES.OUTER.MEAN" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV "SUBTRACTION.FRES.OUTER.STDEV" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN "SUBTRACTION.FRES.TOTAL.MEAN" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV "SUBTRACTION.FRES.TOTAL.STDEV" // RMS stamp deviation // Derive QA information about the subtraction Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 29543) @@ -19,6 +19,10 @@ //#define USE_WEIGHT // Include weight (1/variance) in equation? -//#define USE_WINDOW // Include weight (1/variance) in equation? - + +// XXX TEST: +//# define USE_WINDOW // window to avoid neighbor contamination + +# define PENALTY false +# define MOMENTS (!PENALTY) ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -27,7 +31,7 @@ // Calculate the least-squares matrix and vector -static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated - psVector *vector, // Least-squares vector, updated - double *norm, // Normalisation, updated +static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated + double normValue, // Normalisation, supplied const psKernel *input, // Input image (target) const psKernel *reference, // Reference image (convolution source) @@ -37,8 +41,5 @@ const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - int normWindow1, // Window (half-)size for normalisation measurement - int normWindow2, // Window (half-)size for normalisation measurement - const pmSubtractionEquationCalculationMode mode + int footprint // (Half-)Size of stamp ) { @@ -50,13 +51,21 @@ int numKernels = kernels->num; // Number of kernels - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background int spatialOrder = kernels->spatialOrder; // Order of spatial variation int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms double poly[numPoly]; // Polynomial terms double poly2[numPoly][numPoly]; // Polynomial-polynomial values + int numParams = numKernels * numPoly; + + psAssert(matrix && + matrix->type.type == PS_TYPE_F64 && + matrix->numCols == numParams && + matrix->numRows == numParams, + "Least-squares matrix is bad."); + psAssert(vector && + vector->type.type == PS_TYPE_F64 && + vector->n == numParams, + "Least-squares vector is bad."); // Evaluate polynomial-polynomial terms - // XXX we can skip this if we are not calculating kernel coeffs for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { @@ -84,6 +93,4 @@ // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] - // normalization - // bg 0, bg 1, bg 2 (only 0 is currently used?) for (int i = 0; i < numKernels; i++) { @@ -107,18 +114,15 @@ // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double value = sumCC * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex] = value; - matrix->data.F64[jIndex][iIndex] = value; - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double value = sumCC * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex] = value; + matrix->data.F64[jIndex][iIndex] = value; + } + } } double sumRC = 0.0; // Sum of the reference-convolution products double sumIC = 0.0; // Sum of the input-convolution products - double sumC = 0.0; // Sum of the convolution for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -128,10 +132,8 @@ double ic = in * conv; double rc = ref * conv; - double c = conv; if (weight) { float wtVal = weight->kernel[y][x]; ic *= wtVal; rc *= wtVal; - c *= wtVal; } if (window) { @@ -139,97 +141,14 @@ ic *= winVal; rc *= winVal; - c *= winVal; } sumIC += ic; sumRC += rc; - sumC += c; - } - } + } + } + // Spatial variation for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - double normTerm = sumRC * poly[iTerm]; - double bgTerm = sumC * poly[iTerm]; - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = normTerm; - matrix->data.F64[normIndex][iIndex] = normTerm; - } - if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = bgTerm; - matrix->data.F64[bgIndex][iIndex] = bgTerm; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = sumIC * poly[iTerm]; - if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { - // subtract norm * sumRC * poly[iTerm] - psAssert (kernels->solution1, "programming error: define solution first!"); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation - vector->data.F64[iIndex] -= norm * normTerm; - } - } - } - } - - double sumRR = 0.0; // Sum of the reference product - double sumIR = 0.0; // Sum of the input-reference product - double sum1 = 0.0; // Sum of the background - double sumR = 0.0; // Sum of the reference - double sumI = 0.0; // Sum of the input - double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double in = input->kernel[y][x]; - double ref = reference->kernel[y][x]; - double ir = in * ref; - double rr = PS_SQR(ref); - double one = 1.0; - - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { - normI1 += ref; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { - normI2 += in; - } - - if (weight) { - float wtVal = weight->kernel[y][x]; - rr *= wtVal; - ir *= wtVal; - in *= wtVal; - ref *= wtVal; - one *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - rr *= winVal; - ir *= winVal; - in *= winVal; - ref *= winVal; - one *= winVal; - } - sumRR += rr; - sumIR += ir; - sumR += ref; - sumI += in; - sum1 += one; - } - } - - *norm = normI2 / normI1; - - fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); - - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumRR; - vector->data.F64[normIndex] = sumIR; - // subtract sum over kernels * kernel solution - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI; - } - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[normIndex][bgIndex] = sumR; - matrix->data.F64[bgIndex][normIndex] = sumR; + vector->data.F64[iIndex] = (sumIC - normValue*sumRC) * poly[iTerm]; + } } @@ -255,33 +174,31 @@ // Calculate the least-squares matrix and vector for dual convolution -static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated - psVector *vector, // Least-squares vector, updated - double *norm, // Normalisation, updated - const psKernel *image1, // Image 1 - const psKernel *image2, // Image 2 +// XXX we could avoid calculating these values on successive passes *if* the stamp has not changed. +static bool calculateDualMatrixVector(pmSubtractionStamp *stamp, // stamp of interest + double normValue, // Normalisation, updated + double normValue2, // Normalisation, updated 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 const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - int normWindow1, // Window (half-)size for normalisation measurement - int normWindow2, // Window (half-)size for normalisation measurement - const pmSubtractionEquationCalculationMode mode + int footprint // (Half-)Size of stamp ) { - int numKernels = kernels->num; // Number of kernels - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numKernels = kernels->num; // Number of kernels + int spatialOrder = kernels->spatialOrder; // Order of spatial variation int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms double poly[numPoly]; // Polynomial terms double poly2[numPoly][numPoly]; // Polynomial-polynomial values - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters - int numParams2 = numKernels * numPoly; // Number of additional parameters for dual - int numDual = numParams + numParams2; // Total number of parameters for dual + int numParams = numKernels * numPoly; // Number of regular parameters + int numParams2 = numKernels * numPoly; // Number of additional parameters for dual + int numDual = numParams + numParams2; // Total number of parameters for dual + + psImage *matrix = stamp->matrix; // Least-squares matrix, updated + psVector *vector = stamp->vector; // Least-squares vector, updated + psKernel *image1 = stamp->image1; // Image 1 + psKernel *image2 = stamp->image2; // Image 2 + psArray *convolutions1 = stamp->convolutions1; // Convolutions of image 1 for each kernel + psArray *convolutions2 = stamp->convolutions2; // Convolutions of image 2 for each kernel psAssert(matrix && @@ -290,4 +207,5 @@ matrix->numRows == numDual, "Least-squares matrix is bad."); + psAssert(vector && vector->type.type == PS_TYPE_F64 && @@ -318,4 +236,18 @@ } + // the order of the elements in the matrix and vector is: + // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0] + // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] + // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] + + // for DUAL convolution analysis, we apply the normalization to I1 as follows: + // norm = I2 / I1 + // + // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i + // I2c = I2 + \Sum_i b_i I2 \cross k_i + + // we cannot absorb the normalization into a_i until the analysis is complete, or the + // second moment terms are incorrectly calculated. + for (int i = 0; i < numKernels; i++) { psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i @@ -328,9 +260,14 @@ double sumBB = 0.0; // Sum of convolution products between image 2 double sumAB = 0.0; // Sum of convolution products across images 1 and 2 + + double MxxAA = 0.0; + double MyyAA = 0.0; + double MxxBB = 0.0; + double MyyBB = 0.0; for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x]; + 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]; + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; if (weight) { float wtVal = weight->kernel[y][x]; @@ -348,27 +285,76 @@ sumBB += bb; sumAB += ab; - } - } + + if (MOMENTS) { + MxxAA += x*x*aa; + MyyAA += y*y*aa; + MxxBB += x*x*bb; + MyyBB += y*y*bb; + } + } + } + + if (MOMENTS) { + MxxAA /= stamp->normSquare1 * PS_SQR(normValue); + MyyAA /= stamp->normSquare1 * PS_SQR(normValue); + MxxBB /= stamp->normSquare2; + MyyBB /= stamp->normSquare2; + } + + // XXX this makes the Chisq portion independent of the normalization and star flux + // but may be mis-scaling between stars of different fluxes + sumAA /= PS_SQR(stamp->normI2); + sumAB /= PS_SQR(stamp->normI2); + sumBB /= PS_SQR(stamp->normI2); + + // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB); // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double aa = sumAA * poly2[iTerm][jTerm]; - double bb = sumBB * poly2[iTerm][jTerm]; - double ab = sumAB * poly2[iTerm][jTerm]; - - matrix->data.F64[iIndex][jIndex] = aa; - matrix->data.F64[jIndex][iIndex] = aa; - - matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; - matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; - - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double aa = sumAA * poly2[iTerm][jTerm]; + double bb = sumBB * poly2[iTerm][jTerm]; + double ab = sumAB * poly2[iTerm][jTerm]; + + matrix->data.F64[iIndex][jIndex] = aa; + matrix->data.F64[jIndex][iIndex] = aa; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; + matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; + + // add in second moments + if (MOMENTS) { + matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA; + matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA; + + matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA; + matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB; + matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB; + + matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB; + matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB; + + // XXX this uses the single moments, first try + // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j]; + // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j]; + // + // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j]; + // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j]; + // + // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j]; + // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j]; + // + // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j]; + // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j]; + } + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } + + // we need to calculate the lower-diagonal AB elements since they are not symmetric for A <-> B for (int j = 0; j < i; j++) { psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j @@ -376,5 +362,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; if (weight) { ab *= weight->kernel[y][x]; @@ -387,13 +373,15 @@ } + // XXX this makes the Chisq portion independent of the normalization and star flux + // but may be mis-scaling between stars of different fluxes + sumAB /= PS_SQR(stamp->normI2); + // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double ab = sumAB * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double ab = sumAB * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } } } @@ -402,8 +390,10 @@ double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) - double sumA = 0.0; // Sum of A (for matrix, background) double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) - double sumB = 0.0; // Sum of B products (for matrix, background) - double sumI2 = 0.0; // Sum of I_2 (for vector, background) + + double MxxAI1 = 0.0; + double MyyAI1 = 0.0; + double MxxBI2 = 0.0; + double MyyBI2 = 0.0; for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -413,8 +403,8 @@ float i2 = image2->kernel[y][x]; - double ai2 = a * i2; + double ai2 = a * i2 * normValue; double bi2 = b * i2; - double ai1 = a * i1; - double bi1 = b * i1; + double ai1 = a * i1 * PS_SQR(normValue); + double bi1 = b * i1 * normValue; if (weight) { @@ -424,7 +414,4 @@ ai1 *= wtVal; bi1 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; } if (window) { @@ -434,17 +421,35 @@ ai1 *= wtVal; bi1 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; } sumAI2 += ai2; sumBI2 += bi2; sumAI1 += ai1; - sumA += a; sumBI1 += bi1; - sumB += b; - sumI2 += i2; - } - } + + if (MOMENTS) { + MxxAI1 += x*x*ai1; + MyyAI1 += y*y*ai1; + MxxBI2 += x*x*bi2; + MyyBI2 += y*y*bi2; + } + } + } + + if (MOMENTS) { + MxxAI1 /= stamp->normSquare1 * PS_SQR(normValue); + MyyAI1 /= stamp->normSquare1 * PS_SQR(normValue); + MxxBI2 /= stamp->normSquare2; + MyyBI2 /= stamp->normSquare2; + } + + // fprintf (stderr, "i : %d : M(xx,yy)(AI1,BI2) : %f %f %f %f\n", i, MxxAI1, MyyAI1, MxxBI2, MyyBI2); + + // XXX this makes the Chisq portion independent of the normalization and star flux + // but may be mis-scaling between stars of different fluxes + sumAI1 /= PS_SQR(stamp->normI2); + sumBI1 /= PS_SQR(stamp->normI2); + sumAI2 /= PS_SQR(stamp->normI2); + sumBI2 /= PS_SQR(stamp->normI2); + // Spatial variation for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { @@ -452,95 +457,19 @@ double bi2 = sumBI2 * poly[iTerm]; double ai1 = sumAI1 * poly[iTerm]; - double a = sumA * poly[iTerm]; double bi1 = sumBI1 * poly[iTerm]; - double b = sumB * poly[iTerm]; - - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = ai1; - matrix->data.F64[normIndex][iIndex] = ai1; - matrix->data.F64[iIndex + numParams][normIndex] = bi1; - matrix->data.F64[normIndex][iIndex + numParams] = bi1; - } - if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = a; - matrix->data.F64[bgIndex][iIndex] = a; - matrix->data.F64[iIndex + numParams][bgIndex] = b; - matrix->data.F64[bgIndex][iIndex + numParams] = b; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = ai2; - vector->data.F64[iIndex + numParams] = bi2; - if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { - // subtract norm * sumRC * poly[iTerm] - psAssert (kernels->solution1, "programming error: define solution first!"); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation - vector->data.F64[iIndex] -= norm * ai1; - vector->data.F64[iIndex + numParams] -= norm * bi1; - } - } - } - } - - double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) - double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) - double sum1 = 0.0; // Sum of 1 (for matrix, background-background) - double sumI2 = 0.0; // Sum of I_2 (for vector, background) - double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) - double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double i1 = image1->kernel[y][x]; - double i2 = image2->kernel[y][x]; - - double i1i1 = i1 * i1; - double one = 1.0; - double i1i2 = i1 * i2; - - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { - normI1 += i1; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { - normI2 += i2; - } - - if (weight) { - float wtVal = weight->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } - sumI1 += i1; - sumI1I1 += i1i1; - sum1 += one; - sumI2 += i2; - sumI1I2 += i1i2; - } - } - - *norm = normI2 / normI1; - fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); - - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumI1I1; - vector->data.F64[normIndex] = sumI1I2; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI2; - } - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[bgIndex][normIndex] = sumI1; - matrix->data.F64[normIndex][bgIndex] = sumI1; + vector->data.F64[iIndex] = ai2 - ai1; + vector->data.F64[iIndex + numParams] = bi2 - bi1; + + // fprintf (stderr, "i : %d : V(I1,I2) : %f %f\n", i, vector->data.F64[iIndex], vector->data.F64[iIndex + numParams]); + + // add in second moments + if (MOMENTS) { + vector->data.F64[iIndex] -= kernels->penalty * MxxAI1; + vector->data.F64[iIndex] -= kernels->penalty * MyyAI1; + + vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2; + vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2; + } + } } @@ -565,12 +494,14 @@ } -#if 1 // Add in penalty term to least-squares vector bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term psVector *vector, // Vector to which to add in penalty term const pmSubtractionKernels *kernels, // Kernel parameters - float norm // Normalisation + float normSquare1, // Normalisation for image 1 + float normSquare2 // Normalisation for image 2 ) { + psAssert (kernels->mode == PM_SUBTRACTION_MODE_DUAL, "only use penalties for dual convolution"); + if (kernels->penalty == 0.0) { return true; @@ -589,6 +520,5 @@ // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] - // [norm] - // [bg] + // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] @@ -600,21 +530,19 @@ // Contribution to chi^2: a_i^2 P_i psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty"); - fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]); - matrix->data.F64[index][index] += norm * penalties1->data.F32[i]; - if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], - matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]); - matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i]; - // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) - // penalties scale with second moments - // - } - } - } - } - - return true; -} -# endif + // fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]); + matrix->data.F64[index][index] += normSquare1 * penalties1->data.F32[i]; + + // fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], + // matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]); + matrix->data.F64[index + numParams][index + numParams] += normSquare2 * penalties2->data.F32[i]; + // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) + // penalties scale with second moments + // + } + } + } + + return true; +} ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -647,12 +575,5 @@ int index, bool wantDual) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); - PS_ASSERT_INT_POSITIVE(index, NAN); -#endif - psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, @@ -662,11 +583,5 @@ double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); -#endif - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation return kernels->solution1->data.F64[normIndex]; @@ -676,11 +591,5 @@ const psImage *polyValues) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); -#endif - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); @@ -690,4 +599,208 @@ // Public functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmSubtractionCalculateMoments( + pmSubtractionKernels *kernels, // Kernels + pmSubtractionStampList *stamps) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + // XXX skip this, right? + return true; + + // these are only used by DUAL mode + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true; + + psTimerStart("pmSubtractionCalculateMoments"); + + int footprint = stamps->footprint; // Half-size of stamps + + // Loop over each stamp and calculate its normalization factor + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; + if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; + + pmSubtractionCalculateMomentsStamp(kernels, stamp, footprint, stamps->normWindow2, stamps->normWindow1); + } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate moments: %f sec", psTimerClear("pmSubtractionCalculateMoments")); + + return true; +} + +bool pmSubtractionCalculateMomentsStamp( + pmSubtractionKernels *kernels, // Kernels + pmSubtractionStamp *stamp, // stamp on which to save normalization) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ) +{ + double Mxx, Myy; + + int numKernels = kernels->num; + + // Generate convolutions: these are generated once and saved + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + psError(psErrorCodeLast(), false, "Unable to convolve stamp"); + return false; + } + + if (!stamp->MxxI1) { + stamp->MxxI1 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + if (!stamp->MyyI1) { + stamp->MyyI1 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + if (!stamp->MxxI2) { + stamp->MxxI2 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + if (!stamp->MyyI2) { + stamp->MyyI2 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + + for (int i = 0; i < numKernels; i++) { + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions1->data[i], footprint, normWindow1); + stamp->MxxI1->data.F32[i] = Mxx / stamp->normI1; + stamp->MyyI1->data.F32[i] = Myy / stamp->normI1; + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions2->data[i], footprint, normWindow2); + stamp->MxxI2->data.F32[i] = Mxx / stamp->normI2; + stamp->MyyI2->data.F32[i] = Myy / stamp->normI2; + } + + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image1, footprint, normWindow1); + stamp->MxxI1raw = Mxx / stamp->normI1; + stamp->MyyI1raw = Myy / stamp->normI1; + + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image2, footprint, normWindow2); + stamp->MxxI2raw = Mxx / stamp->normI2; + stamp->MyyI2raw = Myy / stamp->normI2; + + // fprintf (stderr, "Mxx I1: %f, Myy I1: %f, Mxx I2: %f, Myy I2: %f\n", stamp->MxxI1raw, stamp->MyyI1raw, stamp->MxxI2raw, stamp->MyyI2raw); + + return true; +} + +bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window) { + + double Sxx = 0.0; + double Syy = 0.0; + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + if (PS_SQR(x) + PS_SQR(y) > PS_SQR(window)) continue; + + double flux = image->kernel[y][x]; + + Sxx += PS_SQR(x) * flux; + Syy += PS_SQR(y) * flux; + } + } + *Mxx = Sxx; + *Myy = Syy; + return true; +} + +///--------- + +bool pmSubtractionCalculateNormalization( + pmSubtractionStampList *stamps, + const pmSubtractionMode mode) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + psTimerStart("pmSubtractionCalculateNormalization"); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + psVector *norm2 = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + + int footprint = stamps->footprint; // Half-size of stamps + + // Loop over each stamp and calculate its normalization factor + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; + if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; + + // XXX skip this if we have already calculated it? (stamp->norm does not change, just the median statistic) + // XXX maybe not: the star may have changed for a given stamp -- only if norm is reset to NAN can we do this + if (mode == PM_SUBTRACTION_MODE_2) { + pmSubtractionCalculateNormalizationStamp(stamp, stamp->image2, stamp->image1, footprint, stamps->normWindow2, stamps->normWindow1); + } else { + pmSubtractionCalculateNormalizationStamp(stamp, stamp->image1, stamp->image2, footprint, stamps->normWindow1, stamps->normWindow2); + } + psVectorAppend(norms, stamp->norm); + psVectorAppend(norm2, stamp->normSquare2 / stamp->normSquare1); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(norms); + psFree(norm2); + return false; + } + stamps->normValue = stats->robustMedian; + + psStatsInit(stats); + if (!psVectorStats(stats, norm2, NULL, NULL, 0)) { + psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(norms); + psFree(norm2); + return false; + } + stamps->normValue2 = stats->robustMedian; + + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2); + + psFree(stats); + psFree(norms); + psFree(norm2); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate normalization: %f sec", psTimerClear("pmSubtractionCalculateNormalization")); + + return true; +} + +bool pmSubtractionCalculateNormalizationStamp( + pmSubtractionStamp *stamp, // stamp on which to save normalization) + const psKernel *image1, // Input image (target) + const psKernel *image2, // Reference image (convolution source) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ) +{ + double normI1 = 0.0; // Sum of I_1 within the normalisation window (aperture) + double normI2 = 0.0; // Sum of I_2 within the normalisation window (aperture) + double normSquare1 = 0.0; // Sum of (I_1)^2 within the normalisation window (aperture) + double normSquare2 = 0.0; // Sum of (I_2)^2 within the normalisation window (aperture) + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double im1 = image1->kernel[y][x]; + double im2 = image2->kernel[y][x]; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { + normI1 += im1; + normSquare1 += PS_SQR(im1); + } + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { + normI2 += im2; + normSquare2 += PS_SQR(im2); + } + } + } + stamp->norm = normI2 / normI1; + stamp->normI1 = normI1; + stamp->normI2 = normI2; + stamp->normSquare1 = normSquare1; + stamp->normSquare2 = normSquare2; + + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); + + return true; +} bool pmSubtractionCalculateEquationThread(psThreadJob *job) @@ -715,5 +828,4 @@ int numKernels = kernels->num; // Number of kernel basis functions int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). @@ -722,5 +834,6 @@ // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the // number of coefficients for the spatial polynomial, normalisation and a constant background offset. - int numParams = numKernels * numSpatial + 1 + numBackground; + // XXX we no longer solve for the normalization and background in the matrix inversion + int numParams = numKernels * numSpatial; if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { // An additional image is convolved @@ -790,18 +903,13 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: - status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, - weight, window, stamp->convolutions1, kernels, - polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image2, stamp->image1, + weight, window, stamp->convolutions1, kernels, polyValues, footprint); break; case PM_SUBTRACTION_MODE_2: - status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, - weight, window, stamp->convolutions2, kernels, - polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode); + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2, + weight, window, stamp->convolutions2, kernels, polyValues, footprint); break; case PM_SUBTRACTION_MODE_DUAL: - status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, - stamp->image1, stamp->image2, - weight, window, stamp->convolutions1, stamp->convolutions2, - kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + status = calculateDualMatrixVector(stamp, stamps->normValue, stamps->normValue2, weight, window, kernels, polyValues, footprint); break; default: @@ -928,7 +1036,8 @@ int numKernels = kernels->num; // Number of kernel basis functions int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for - int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution + // XXX int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + // XXX int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for + int numParams = numKernels * numSpatial; // Number of parameters being solved for + int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { // An additional image is convolved @@ -960,5 +1069,6 @@ if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { - // Accumulate the least-squares matricies and vectors + // Accumulate the least-squares matrices and vectors. These are generated for the + // kernel elements, excluding the background and normalization. psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector @@ -966,7 +1076,83 @@ psImageInit(sumMatrix, 0.0); - psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations - int numStamps = 0; // Number of good stamps + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); + numStamps++; + } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); + } + } + +#if 0 + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); + psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); +#endif + + psVector *solution = NULL; // Solution to equation! + solution = psVectorAlloc(numParams, PS_TYPE_F64); + psVectorInit(solution, 0); + + // XXX TEST: try some constraint on the svd solution + // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + // SINGLE solution + if (1) { + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + } else { + psVector *PERM = NULL; + psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix); + solution = psMatrixLUSolution(solution, LU, sumVector, PERM); + psFree (LU); + psFree (PERM); + } + +# if (0) + for (int i = 0; i < solution->n; i++) { + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[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); + } + + 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]; + } + + // Apply the normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + kernels->solution1->data.F64[normIndex] = stamps->normValue; + kernels->solution1->data.F64[bgIndex] = 0.0; + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + + } else { + // Dual convolution solution + // Accumulate the least-squares matrices and vectors. These are generated for the + // kernel elements, excluding the background and normalization. + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix, 0.0); + psVectorInit(sumVector, 0.0); + + int numStamps = 0; // Number of good stamps + double normSquare1 = 0.0; // Sum of (I_1)^2 over stamps + double normSquare2 = 0.0; // Sum of (I_2)^2 over stamps for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -974,5 +1160,8 @@ (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - psVectorAppend(norms, stamp->norm); + + normSquare1 += stamp->normSquare1; + normSquare2 += stamp->normSquare2; + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; @@ -983,7 +1172,13 @@ #if 0 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); + psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); #endif + + if (PENALTY) { + calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2); + } psVector *solution = NULL; // Solution to equation! @@ -991,212 +1186,8 @@ psVectorInit(solution, 0); -#if 0 - // Regular, straight-forward solution - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); -#else - { - // Solve normalisation and background separately - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm - if (!psVectorStats(stats, norms, NULL, NULL, 0)) { - psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); - psFree(stats); - psFree(sumMatrix); - psFree(sumVector); - psFree(norms); - return false; - } - - // double normValue = 1.0; - double normValue = stats->robustMedian; - // double bgValue = 0.0; - - psFree(stats); - -#ifdef TESTING - fprintf(stderr, "Norm: %lf\n", normValue); -#endif - // Solve kernel components - for (int i = 0; i < numSolution1; i++) { - sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; - - sumMatrix->data.F64[i][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i] = 0.0; - } - sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; - sumMatrix->data.F64[bgIndex][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][bgIndex] = 0.0; - - sumMatrix->data.F64[normIndex][normIndex] = 1.0; - sumVector->data.F64[normIndex] = 0.0; - - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - - solution->data.F64[normIndex] = normValue; - } -# endif - -#if (1) - for (int i = 0; i < solution->n; i++) { - fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]); - } -#endif - - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64); - psVectorInit(kernels->solution1, 0.0); - } - - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - int spatialOrder = kernels->spatialOrder; // Order of spatial variation - int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms - for (int i = 0; i < numKernels * numPoly; i++) { - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } - - psFree(norms); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - -#ifdef TESTING - // XXX double-check for NAN in data: - for (int ix = 0; ix < kernels->solution1->n; ix++) { - if (!isfinite(kernels->solution1->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - - } else { - // Dual convolution solution - - // Accumulation of stamp matrices/vectors - psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); - psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); - psImageInit(sumMatrix, 0.0); - psVectorInit(sumVector, 0.0); - - psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations - - int numStamps = 0; // Number of good stamps - for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_USED) { - (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); - (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - - psVectorAppend(norms, stamp->norm); - - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); - numStamps++; - } - } - -#ifdef TESTING - psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); - psVectorWriteFile("sumVector.dat", sumVector); -#endif - -#if 1 - // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); - - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0); -#endif - - psVector *solution = NULL; // Solution to equation! - solution = psVectorAlloc(numParams, PS_TYPE_F64); - psVectorInit(solution, 0); - -#if 0 - // Regular, straight-forward solution - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); -#else - { - // Solve normalisation and background separately - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - -#if 0 - psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); - psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); - - normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; - normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; - normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; - - normVector->data.F64[0] = sumVector->data.F64[normIndex]; - normVector->data.F64[1] = sumVector->data.F64[bgIndex]; - - psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); - - double normValue = normSolution->data.F64[0]; - double bgValue = normSolution->data.F64[1]; - - psFree(normMatrix); - psFree(normVector); - psFree(normSolution); -#endif - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm - if (!psVectorStats(stats, norms, NULL, NULL, 0)) { - psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); - psFree(stats); - psFree(sumMatrix); - psFree(sumVector); - psFree(norms); - return false; - } - - double normValue = stats->robustMedian; - - psFree(stats); - -#ifdef TESTING - fprintf(stderr, "Norm: %lf\n", normValue); -#endif - - // Solve kernel components - for (int i = 0; i < numSolution2; i++) { - sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; - sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; - - sumMatrix->data.F64[i][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i] = 0.0; - - sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; - } - sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; - sumMatrix->data.F64[bgIndex][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][bgIndex] = 0.0; - - sumMatrix->data.F64[normIndex][normIndex] = 1.0; - - sumVector->data.F64[normIndex] = 0.0; - - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - - solution->data.F64[normIndex] = normValue; - } -#endif - - -#if (1) + // DUAL solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + +#if (0) for (int i = 0; i < solution->n; i++) { fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); @@ -1204,11 +1195,6 @@ #endif - psFree(sumMatrix); - psFree(sumVector); - - psFree(norms); - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64); + kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); psVectorInit (kernels->solution1, 0.0); } @@ -1218,30 +1204,27 @@ } - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - for (int i = 0; i < numKernels * numSpatial; i++) { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; - } - } - - - memcpy(kernels->solution1->data.F64, solution->data.F64, - numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], - numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - + // for DUAL convolution analysis, we apply the normalization to I1 as follows: + // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i + // I2c = I2 + \Sum_i b_i I2 \cross k_i + + // We absorb the normalization into a_i after the analysis is complete to be consistent + // with the SINGLE definitions of the convolutions + + int numKernels = kernels->num; + for (int i = 0; i < numKernels * numSpatial; i++) { + // we solve for coefficients + kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue; + kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; + } + + // Apply the normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + kernels->solution1->data.F64[normIndex] = stamps->normValue; + kernels->solution1->data.F64[bgIndex] = 0.0; + + psFree(sumMatrix); + psFree(sumVector); psFree(solution); - } @@ -1265,7 +1248,10 @@ } -bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { - - // XXX measure some useful stats on the residuals +// measure some useful stats on the stamp residuals: +// fResSigma : the residual stdev / total flux +// fResOuter : the residual fabs / total flux for R > 2 pix +// fResTotal : the residual fabs / total flux for R > 0 pix +bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { + float sum = 0.0; float peak = 0.0; @@ -1277,21 +1263,19 @@ } - // only count pixels with more than X% of the source flux - // calculate stdev(dflux) + // init counters + int npix = 0; float dflux1 = 0.0; float dflux2 = 0.0; - int npix = 0; - - float dmax = 0.0; - float dmin = 0.0; + float dOuter = 0.0; + float dTotal = 0.0; for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm); - if (dflux < 0.02*sum) continue; dflux1 += residual->kernel[y][x]; dflux2 += PS_SQR(residual->kernel[y][x]); - dmax = PS_MAX(residual->kernel[y][x], dmax); - dmin = PS_MIN(residual->kernel[y][x], dmin); + dTotal += fabs(residual->kernel[y][x]); + if (hypot(x,y) > 2.0) { + dOuter += fabs(residual->kernel[y][x]); + } npix ++; } @@ -1299,12 +1283,12 @@ float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix)); if (!isfinite(sum)) return false; - if (!isfinite(dmax)) return false; - if (!isfinite(dmin)) return false; if (!isfinite(peak)) return false; - - // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak); - psVectorAppend(fSigRes, sigma/sum); - psVectorAppend(fMaxRes, dmax/peak); - psVectorAppend(fMinRes, dmin/peak); + if (!isfinite(dOuter)) return false; + if (!isfinite(dTotal)) return false; + + // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum); + psVectorAppend(fResSigma, sigma/sum); + psVectorAppend(fResOuter, dOuter/sum); + psVectorAppend(fResTotal, dTotal/sum); return true; } @@ -1329,7 +1313,7 @@ pmSubtractionVisualShowFitInit (stamps); - psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); - psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); - psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); // we want to save the residual images for the 9 brightest stamps. @@ -1443,5 +1427,5 @@ } - pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint); + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint); } else { @@ -1479,5 +1463,5 @@ } - pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint); + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint); } @@ -1558,27 +1542,27 @@ psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psVectorStats (stats, fSigRes, NULL, NULL, 0); - kernels->fSigResMean = stats->robustMedian; - kernels->fSigResStdev = stats->robustStdev; + psVectorStats (stats, fResSigma, NULL, NULL, 0); + kernels->fResSigmaMean = stats->robustMedian; + kernels->fResSigmaStdev = stats->robustStdev; psStatsInit (stats); - psVectorStats (stats, fMaxRes, NULL, NULL, 0); - kernels->fMaxResMean = stats->robustMedian; - kernels->fMaxResStdev = stats->robustStdev; + psVectorStats (stats, fResOuter, NULL, NULL, 0); + kernels->fResOuterMean = stats->robustMedian; + kernels->fResOuterStdev = stats->robustStdev; psStatsInit (stats); - psVectorStats (stats, fMinRes, NULL, NULL, 0); - kernels->fMinResMean = stats->robustMedian; - kernels->fMinResStdev = stats->robustStdev; + psVectorStats (stats, fResTotal, NULL, NULL, 0); + kernels->fResTotalMean = stats->robustMedian; + kernels->fResTotalStdev = stats->robustStdev; // XXX save these values somewhere - psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f", - kernels->fSigResMean, kernels->fSigResStdev, - kernels->fMaxResMean, kernels->fMaxResStdev, - kernels->fMinResMean, kernels->fMinResStdev); - - psFree (fSigRes); - psFree (fMaxRes); - psFree (fMinRes); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f", + kernels->fResSigmaMean, kernels->fResSigmaStdev, + kernels->fResOuterMean, kernels->fResOuterStdev, + kernels->fResTotalMean, kernels->fResTotalStdev); + + psFree (fResSigma); + psFree (fResOuter); + psFree (fResTotal); psFree (stats); } @@ -1586,5 +1570,4 @@ psFree(residual); psFree(polyValues); - return deviations; @@ -1912,312 +1895,2 @@ return true; } - - -# if 0 - -#ifdef TESTING - psFitsWriteImageSimple("A.fits", sumMatrix, NULL); - psVectorWriteFile ("B.dat", sumVector); -#endif - -# define SVD_ANALYSIS 0 -# define COEFF_SIG 0.0 -# define SVD_TOL 0.0 - - // Use SVD to determine the kernel coeffs (and validate) - if (SVD_ANALYSIS) { - - // We have sumVector and sumMatrix. we are trying to solve the following equation: - // sumMatrix * x = sumVector. - - // we can use any standard matrix inversion to solve this. However, the basis - // functions in general have substantial correlation, so that the solution may be - // somewhat poorly determined or unstable. If not numerically ill-conditioned, the - // system of equations may be statistically ill-conditioned. Noise in the image - // will drive insignificant, but correlated, terms in the solution. To avoid these - // problems, we can use SVD to identify numerically unconstrained values and to - // avoid statistically badly determined value. - - // A = sumMatrix, B = sumVector - // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T - // x = V (1/w) (U^T B) - // \sigma_x = sqrt(diag(A^{-1})) - // solve for x and A^{-1} to get x & dx - // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 - // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 - - // If I use the SVD trick to re-condition the matrix, I need to break out the - // kernel and normalization terms from the background term. - // XXX is this true? or was this due to an error in the analysis? - - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - - // now pull out the kernel elements into their own square matrix - psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); - psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { - if (ix == bgIndex) continue; - for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { - if (iy == bgIndex) continue; - kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; - ky++; - } - kernelVector->data.F64[kx] = sumVector->data.F64[ix]; - kx++; - } - - psImage *U = NULL; - psImage *V = NULL; - psVector *w = NULL; - if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { - psError(psErrorCodeLast(), false, "failed to perform SVD on sumMatrix\n"); - return NULL; - } - - // calculate A_inverse: - // Ainv = V * w * U^T - psImage *wUt = p_pmSubSolve_wUt (w, U); - psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); - psImage *Xvar = NULL; - psFree (wUt); - -# ifdef TESTING - // kernel terms: - for (int i = 0; i < w->n; i++) { - fprintf (stderr, "w: %f\n", w->data.F64[i]); - } -# endif - // loop over w adding in more and more of the values until chisquare is no longer - // dropping significantly. - // XXX this does not seem to work very well: we seem to need all terms even for - // simple cases... - - psVector *Xsvd = NULL; - { - psVector *Ax = NULL; - psVector *UtB = NULL; - psVector *wUtB = NULL; - - psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); - psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); - psVectorInit (wMask, 1); // start by masking everything - - double chiSquareLast = NAN; - int maxWeight = 0; - - double Axx, Bx, y2; - - // XXX this is an attempt to exclude insignificant modes. - // it was not successful with the ISIS kernel set: removing even - // the least significant mode leaves additional ringing / noise - // because the terms are so coupled. - for (int k = 0; false && (k < w->n); k++) { - - // unmask the k-th weight - wMask->data.U8[k] = 0; - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; - chiSquareLast = chiSquare; - - // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); - if (k && !maxWeight && (deltaChi < 1.0)) { - maxWeight = k; - } - } - - // keep all terms or we get extra ringing - maxWeight = w->n; - psVectorInit (wMask, 1); - for (int k = 0; k < maxWeight; k++) { - wMask->data.U8[k] = 0; - } - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); - - // re-calculate A^{-1} to get new variances: - // Ainv = V * w * U^T - // XXX since we keep all terms, this is identical to Ainv - psImage *wUt = p_pmSubSolve_wUt (wApply, U); - Xvar = p_pmSubSolve_VwUt (V, wUt); - psFree (wUt); - - psFree (Ax); - psFree (UtB); - psFree (wUtB); - psFree (wApply); - psFree (wMask); - } - - // copy the kernel solutions to the full solution vector: - solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { - if (ix == bgIndex) { - solution->data.F64[ix] = 0; - solutionErr->data.F64[ix] = 0.001; - continue; - } - solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); - solution->data.F64[ix] = Xsvd->data.F64[kx]; - kx++; - } - - psFree (kernelMatrix); - psFree (kernelVector); - - psFree (U); - psFree (V); - psFree (w); - - psFree (Ainv); - psFree (Xsvd); - } else { - psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); - if (!luMatrix) { - psError(PM_ERR_DATA, true, "LU Decomposition of least-squares matrix failed.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - psFree(luMatrix); - psFree(permutation); - return NULL; - } - - solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); - psFree(luMatrix); - psFree(permutation); - if (!solution) { - psError(PM_ERR_DATA, true, "Failed to solve the least-squares system.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - return NULL; - } - - // XXX LUD does not provide A^{-1}? fake the error for now - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - for (int ix = 0; ix < sumVector->n; ix++) { - solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; - } - } - - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); - psVectorInit (kernels->solution1, 0.0); - } - - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - int spatialOrder = kernels->spatialOrder; // Order of spatial variation - int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms - for (int i = 0; i < numKernels * numPoly; i++) { - // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); - if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { - // XXX fprintf (stderr, "drop\n"); - kernels->solution1->data.F64[i] = 0.0; - } else { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } - } - // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); - // fprintf (stderr, "chi square Brute = %f\n", chiSquare); - - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); -# endif - -#ifdef TESTING - // XXX double-check for NAN in data: - for (int iy = 0; iy < stamp->matrix->numRows; iy++) { - for (int ix = 0; ix < stamp->matrix->numCols; ix++) { - if (!isfinite(stamp->matrix->data.F64[iy][ix])) { - fprintf (stderr, "WARNING: NAN in matrix\n"); - } - } - } - for (int ix = 0; ix < stamp->vector->n; ix++) { - if (!isfinite(stamp->vector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } - { - psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); - psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); - psFree(inverse); - } - { - psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); - psImage *Xt = psMatrixTranspose(NULL, X); - psImage *XtX = psMatrixMultiply(NULL, Xt, X); - psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); - psFree(X); - psFree(Xt); - psFree(XtX); - } -#endif - Index: trunk/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 29543) @@ -65,3 +65,30 @@ ); +bool pmSubtractionCalculateNormalization( + pmSubtractionStampList *stamps, + const pmSubtractionMode mode); + +bool pmSubtractionCalculateNormalizationStamp( + pmSubtractionStamp *stamp, // stamp on which to save normalization) + const psKernel *input, // Input image (target) + const psKernel *reference, // Reference image (convolution source) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ); + +bool pmSubtractionCalculateMoments( + pmSubtractionKernels *kernels, // Kernels + pmSubtractionStampList *stamps); + +bool pmSubtractionCalculateMomentsStamp( + pmSubtractionKernels *kernels, // Kernels + pmSubtractionStamp *stamp, // stamp on which to save normalization) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ); + +bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window); + #endif Index: trunk/psModules/src/imcombine/pmSubtractionEquation.v0.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 29543) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 29543) @@ -0,0 +1,2238 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmErrorCodes.h" +#include "pmSubtraction.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionStamps.h" +#include "pmSubtractionThreads.h" + +#include "pmSubtractionEquation.h" +#include "pmSubtractionVisual.h" + +//#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? + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Calculate the least-squares matrix and vector +static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated + double *norm, // Normalisation, updated + const psKernel *input, // Input image (target) + const psKernel *reference, // Reference image (convolution source) + const psKernel *weight, // Weight image + const psKernel *window, // Window image + const psArray *convolutions, // Convolutions for each kernel + const pmSubtractionKernels *kernels, // Kernels + const psImage *polyValues, // Spatial polynomial values + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2, // Window (half-)size for normalisation measurement + const pmSubtractionEquationCalculationMode mode + ) +{ + // (I - R * sum_i a_i k_i - g) (R * k_j) = 0 + // I C_j = sum_i C_i C_j + + // Background: C_i = 1.0 + // Normalisation: C_i = R + + int numKernels = kernels->num; // Number of kernels + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + double poly[numPoly]; // Polynomial terms + double poly2[numPoly][numPoly]; // Polynomial-polynomial values + + // Evaluate polynomial-polynomial terms + // XXX we can skip this if we are not calculating kernel coeffs + for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { + double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial + poly[iIndex] = iPoly; + for (int jyOrder = 0, jIndex = 0; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex++) { + double jPoly = polyValues->data.F64[jyOrder][jxOrder]; + poly2[iIndex][jIndex] = iPoly * jPoly; + } + } + } + } + + // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we + // choose to calculate + psImageInit(matrix, 0.0); + psVectorInit(vector, 1.0); + for (int i = 0; i < matrix->numCols; i++) { + matrix->data.F64[i][i] = 1.0; + } + + // the order of the elements in the matrix and vector is: + // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0] + // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] + // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] + // normalization + // bg 0, bg 1, bg 2 (only 0 is currently used?) + + for (int i = 0; i < numKernels; i++) { + psKernel *iConv = convolutions->data[i]; // Convolution for index i + for (int j = i; j < numKernels; j++) { + psKernel *jConv = convolutions->data[j]; // Convolution for index j + + double sumCC = 0.0; // Sum of convolution products + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; + if (weight) { + cc *= weight->kernel[y][x]; + } + if (window) { + cc *= window->kernel[y][x]; + } + sumCC += cc; + } + } + + // Spatial variation of kernel coeffs + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double value = sumCC * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex] = value; + matrix->data.F64[jIndex][iIndex] = value; + } + } + } + } + + double sumRC = 0.0; // Sum of the reference-convolution products + double sumIC = 0.0; // Sum of the input-convolution products + double sumC = 0.0; // Sum of the convolution + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float conv = iConv->kernel[y][x]; + float in = input->kernel[y][x]; + float ref = reference->kernel[y][x]; + double ic = in * conv; + double rc = ref * conv; + double c = conv; + if (weight) { + float wtVal = weight->kernel[y][x]; + ic *= wtVal; + rc *= wtVal; + c *= wtVal; + } + if (window) { + float winVal = window->kernel[y][x]; + ic *= winVal; + rc *= winVal; + c *= winVal; + } + sumIC += ic; + sumRC += rc; + sumC += c; + } + } + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + double normTerm = sumRC * poly[iTerm]; + double bgTerm = sumC * poly[iTerm]; + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][normIndex] = normTerm; + matrix->data.F64[normIndex][iIndex] = normTerm; + } + if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][bgIndex] = bgTerm; + matrix->data.F64[bgIndex][iIndex] = bgTerm; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + vector->data.F64[iIndex] = sumIC * poly[iTerm]; + if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { + // subtract norm * sumRC * poly[iTerm] + psAssert (kernels->solution1, "programming error: define solution first!"); + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation + vector->data.F64[iIndex] -= norm * normTerm; + } + } + } + } + + double sumRR = 0.0; // Sum of the reference product + double sumIR = 0.0; // Sum of the input-reference product + double sum1 = 0.0; // Sum of the background + double sumR = 0.0; // Sum of the reference + double sumI = 0.0; // Sum of the input + double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double in = input->kernel[y][x]; + double ref = reference->kernel[y][x]; + double ir = in * ref; + double rr = PS_SQR(ref); + double one = 1.0; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { + normI1 += ref; + } + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { + normI2 += in; + } + + if (weight) { + float wtVal = weight->kernel[y][x]; + rr *= wtVal; + ir *= wtVal; + in *= wtVal; + ref *= wtVal; + one *= wtVal; + } + if (window) { + float winVal = window->kernel[y][x]; + rr *= winVal; + ir *= winVal; + in *= winVal; + ref *= winVal; + one *= winVal; + } + sumRR += rr; + sumIR += ir; + sumR += ref; + sumI += in; + sum1 += one; + } + } + + *norm = normI2 / normI1; + + fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); + + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + matrix->data.F64[normIndex][normIndex] = sumRR; + vector->data.F64[normIndex] = sumIR; + // subtract sum over kernels * kernel solution + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + matrix->data.F64[bgIndex][bgIndex] = sum1; + vector->data.F64[bgIndex] = sumI; + } + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { + matrix->data.F64[normIndex][bgIndex] = sumR; + matrix->data.F64[bgIndex][normIndex] = sumR; + } + + // check for any NAN values in the result, skip if found: + for (int iy = 0; iy < matrix->numRows; iy++) { + for (int ix = 0; ix < matrix->numCols; ix++) { + if (!isfinite(matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + return false; + } + } + } + for (int ix = 0; ix < vector->n; ix++) { + if (!isfinite(vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + return false; + } + } + + return true; +} + + +// Calculate the least-squares matrix and vector for dual convolution +static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated + double *norm, // Normalisation, updated + const psKernel *image1, // Image 1 + 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 + const pmSubtractionKernels *kernels, // Kernels + const psImage *polyValues, // Spatial polynomial values + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2, // Window (half-)size for normalisation measurement + const pmSubtractionEquationCalculationMode mode + ) +{ + int numKernels = kernels->num; // Number of kernels + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + double poly[numPoly]; // Polynomial terms + double poly2[numPoly][numPoly]; // Polynomial-polynomial values + + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters + int numParams2 = numKernels * numPoly; // Number of additional parameters for dual + int numDual = numParams + numParams2; // Total number of parameters for dual + + psAssert(matrix && + matrix->type.type == PS_TYPE_F64 && + matrix->numCols == numDual && + matrix->numRows == numDual, + "Least-squares matrix is bad."); + psAssert(vector && + vector->type.type == PS_TYPE_F64 && + vector->n == numDual, + "Least-squares vector is bad."); + + // Evaluate polynomial-polynomial terms + for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { + double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial + poly[iIndex] = iPoly; + for (int jyOrder = 0, jIndex = 0; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex++) { + double jPoly = polyValues->data.F64[jyOrder][jxOrder]; + poly2[iIndex][jIndex] = iPoly * jPoly; + } + } + } + } + + + // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we + // choose to calculate + psImageInit(matrix, 0.0); + psVectorInit(vector, 1.0); + for (int i = 0; i < matrix->numCols; i++) { + matrix->data.F64[i][i] = 1.0; + } + + for (int i = 0; i < numKernels; i++) { + psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i + psKernel *iConv2 = convolutions2->data[i]; // Convolution 2 for index i + for (int j = i; j < numKernels; j++) { + psKernel *jConv1 = convolutions1->data[j]; // Convolution 1 for index j + psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j + + double sumAA = 0.0; // Sum of convolution products between image 1 + double sumBB = 0.0; // Sum of convolution products between image 2 + double sumAB = 0.0; // Sum of convolution products across images 1 and 2 + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x]; + double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; + if (weight) { + float wtVal = weight->kernel[y][x]; + aa *= wtVal; + bb *= wtVal; + ab *= wtVal; + } + if (window) { + float wtVal = window->kernel[y][x]; + aa *= wtVal; + bb *= wtVal; + ab *= wtVal; + } + sumAA += aa; + sumBB += bb; + sumAB += ab; + } + } + + // Spatial variation of kernel coeffs + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double aa = sumAA * poly2[iTerm][jTerm]; + double bb = sumBB * poly2[iTerm][jTerm]; + double ab = sumAB * poly2[iTerm][jTerm]; + + matrix->data.F64[iIndex][jIndex] = aa; + matrix->data.F64[jIndex][iIndex] = aa; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; + matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; + + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } + } + for (int j = 0; j < i; j++) { + psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j + double sumAB = 0.0; // Sum of convolution products for matrix C + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; + if (weight) { + ab *= weight->kernel[y][x]; + } + if (window) { + ab *= window->kernel[y][x]; + } + sumAB += ab; + } + } + + // Spatial variation of kernel coeffs + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double ab = sumAB * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } + } + + double sumAI2 = 0.0; // Sum of A.I_2 products (for vector) + double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) + double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) + double sumA = 0.0; // Sum of A (for matrix, background) + double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) + double sumB = 0.0; // Sum of B products (for matrix, background) + double sumI2 = 0.0; // Sum of I_2 (for vector, background) + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double a = iConv1->kernel[y][x]; + double b = iConv2->kernel[y][x]; + float i1 = image1->kernel[y][x]; + float i2 = image2->kernel[y][x]; + + double ai2 = a * i2; + double bi2 = b * i2; + double ai1 = a * i1; + double bi1 = b * i1; + + if (weight) { + float wtVal = weight->kernel[y][x]; + ai2 *= wtVal; + bi2 *= wtVal; + ai1 *= wtVal; + bi1 *= wtVal; + a *= wtVal; + b *= wtVal; + i2 *= wtVal; + } + if (window) { + float wtVal = window->kernel[y][x]; + ai2 *= wtVal; + bi2 *= wtVal; + ai1 *= wtVal; + bi1 *= wtVal; + a *= wtVal; + b *= wtVal; + i2 *= wtVal; + } + sumAI2 += ai2; + sumBI2 += bi2; + sumAI1 += ai1; + sumA += a; + sumBI1 += bi1; + sumB += b; + sumI2 += i2; + } + } + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + double ai2 = sumAI2 * poly[iTerm]; + double bi2 = sumBI2 * poly[iTerm]; + double ai1 = sumAI1 * poly[iTerm]; + double a = sumA * poly[iTerm]; + double bi1 = sumBI1 * poly[iTerm]; + double b = sumB * poly[iTerm]; + + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][normIndex] = ai1; + matrix->data.F64[normIndex][iIndex] = ai1; + matrix->data.F64[iIndex + numParams][normIndex] = bi1; + matrix->data.F64[normIndex][iIndex + numParams] = bi1; + } + if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][bgIndex] = a; + matrix->data.F64[bgIndex][iIndex] = a; + matrix->data.F64[iIndex + numParams][bgIndex] = b; + matrix->data.F64[bgIndex][iIndex + numParams] = b; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + vector->data.F64[iIndex] = ai2; + vector->data.F64[iIndex + numParams] = bi2; + if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { + // subtract norm * sumRC * poly[iTerm] + psAssert (kernels->solution1, "programming error: define solution first!"); + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation + vector->data.F64[iIndex] -= norm * ai1; + vector->data.F64[iIndex + numParams] -= norm * bi1; + } + } + } + } + + double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) + double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) + double sum1 = 0.0; // Sum of 1 (for matrix, background-background) + double sumI2 = 0.0; // Sum of I_2 (for vector, background) + double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) + double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double i1 = image1->kernel[y][x]; + double i2 = image2->kernel[y][x]; + + double i1i1 = i1 * i1; + double one = 1.0; + double i1i2 = i1 * i2; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { + normI1 += i1; + } + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { + normI2 += i2; + } + + if (weight) { + float wtVal = weight->kernel[y][x]; + i1 *= wtVal; + i1i1 *= wtVal; + one *= wtVal; + i2 *= wtVal; + i1i2 *= wtVal; + } + if (window) { + float wtVal = window->kernel[y][x]; + i1 *= wtVal; + i1i1 *= wtVal; + one *= wtVal; + i2 *= wtVal; + i1i2 *= wtVal; + } + sumI1 += i1; + sumI1I1 += i1i1; + sum1 += one; + sumI2 += i2; + sumI1I2 += i1i2; + } + } + + *norm = normI2 / normI1; + fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); + + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + matrix->data.F64[normIndex][normIndex] = sumI1I1; + vector->data.F64[normIndex] = sumI1I2; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + matrix->data.F64[bgIndex][bgIndex] = sum1; + vector->data.F64[bgIndex] = sumI2; + } + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { + matrix->data.F64[bgIndex][normIndex] = sumI1; + matrix->data.F64[normIndex][bgIndex] = sumI1; + } + + // check for any NAN values in the result, skip if found: + for (int iy = 0; iy < matrix->numRows; iy++) { + for (int ix = 0; ix < matrix->numCols; ix++) { + if (!isfinite(matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + return false; + } + } + } + for (int ix = 0; ix < vector->n; ix++) { + if (!isfinite(vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + return false; + } + } + + + return true; +} + +#if 1 +// Add in penalty term to least-squares vector +bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term + psVector *vector, // Vector to which to add in penalty term + const pmSubtractionKernels *kernels, // Kernel parameters + float norm // Normalisation + ) +{ + if (kernels->penalty == 0.0) { + return true; + } + + psVector *penalties1 = kernels->penalties1; // Penalties for each kernel component (input) + psVector *penalties2 = kernels->penalties2; // Penalties for each kernel component (ref) + + int spatialOrder = kernels->spatialOrder; // Order of spatial variations + int numKernels = kernels->num; // Number of kernel components + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations + int numParams = numKernels * numSpatial; // Number of kernel parameters + + // order is : + // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0] + // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] + // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] + // [norm] + // [bg] + // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] + // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] + // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1] + + for (int i = 0; i < numKernels; i++) { + for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { + // Contribution to chi^2: a_i^2 P_i + psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty"); + fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]); + matrix->data.F64[index][index] += norm * penalties1->data.F32[i]; + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], + matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]); + matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i]; + // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) + // penalties scale with second moments + // + } + } + } + } + + return true; +} +# endif + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Semi-public functions +// XXX We might like to define these functions as "extern inline" but gcc currently doesn't handle this in c99 +// mode. See http://gcc.gnu.org/ml/gcc/2006-11/msg00006.html +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Calculate the value of a polynomial, specified by coefficients and polynomial values +double p_pmSubtractionCalculatePolynomial(const psVector *coeff, // Coefficients + const psImage *polyValues, // Polynomial values + int order, // Order of polynomials + int index, // Index at which to begin + int step // Step between subsequent indices + ) +{ + double sum = 0.0; // Value of the polynomial sum + for (int yOrder = 0; yOrder <= order; yOrder++) { + for (int xOrder = 0; xOrder <= order - yOrder; xOrder++, index += step) { + + assert(index < coeff->n); + + sum += coeff->data.F64[index] * polyValues->data.F64[yOrder][xOrder]; + } + } + return sum; +} + +double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, const psImage *polyValues, + int index, bool wantDual) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); + PS_ASSERT_INT_POSITIVE(index, NAN); +#endif + + psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector + return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, + kernels->num); +} + +double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); +#endif + + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + return kernels->solution1->data.F64[normIndex]; +} + +double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, + const psImage *polyValues) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); +#endif + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmSubtractionCalculateEquationThread(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + + pmSubtractionStampList *stamps = job->args->data[0]; // List of stamps + 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) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PS_ASSERT_INT_NONNEGATIVE(index, false); + PS_ASSERT_INT_LESS_THAN(index, stamps->num, false); + + int footprint = stamps->footprint; // Half-size of stamps + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numKernels = kernels->num; // Number of kernel basis functions + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + + // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). + // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3 + + // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the + // number of coefficients for the spatial polynomial, normalisation and a constant background offset. + int numParams = numKernels * numSpatial + 1 + numBackground; + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + // An additional image is convolved + numParams += numKernels * numSpatial; + } + + pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest + 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? + if (new) { + stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64); + } +#ifdef TESTING + psImageInit(stamp->matrix, NAN); + psVectorInit(stamp->vector, NAN); +#endif + + bool status; // Status of least-squares matrix/vector calculation + + psKernel *weight = NULL; + psKernel *window = NULL; + +#ifdef USE_WEIGHT + weight = stamp->weight; +#endif +#ifdef USE_WINDOW + window = stamps->window; +#endif + + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, + weight, window, stamp->convolutions1, kernels, + polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + break; + case PM_SUBTRACTION_MODE_2: + status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, + weight, window, stamp->convolutions2, kernels, + polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode); + break; + case PM_SUBTRACTION_MODE_DUAL: + status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, + stamp->image1, stamp->image2, + weight, window, stamp->convolutions1, stamp->convolutions2, + kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + if (!status) { + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + psWarning("Rejecting stamp %d (%d,%d) because of bad equation", + index, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); + } else { + stamp->status = PM_SUBTRACTION_STAMP_USED; + } + +#ifdef TESTING + { + psString matrixName = NULL; + psStringAppend(&matrixName, "matrix_%d.fits", index); + psFits *matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stamp->matrix, 0, NULL); + psFitsClose(matrixFile); + + matrixName = NULL; + psStringAppend(&matrixName, "vector_%d.fits", index); + psImage *dummy = psImageAlloc(stamp->vector->n, 1, PS_TYPE_F64); + memcpy(dummy->data.F64[0], stamp->vector->data.F64, + PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector->n); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL); + psFree(dummy); + psFitsClose(matrixFile); + } +#endif + + psFree(polyValues); + + return true; +} + +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, + const pmSubtractionEquationCalculationMode mode) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + + psTimerStart("pmSubtractionCalculateEquation"); + + // We iterate over each stamp, allocate the matrix and vectors if + // necessary, and then calculate those matrix/vectors. + 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; + } + + if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) { + psAbort ("bad stamp"); + } + if (!isfinite(stamp->x) && !isfinite(stamp->y)) { + psAbort ("bad stamp"); + } + + if (pmSubtractionThreaded()) { + psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION"); + psArrayAdd(job->args, 1, stamps); + 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); + } + } + + if (!psThreadPoolWait(true)) { + psError(psErrorCodeLast(), false, "Error waiting for threads."); + return false; + } + + pmSubtractionVisualPlotLeastSquares(stamps); + pmSubtractionVisualShowKernels((pmSubtractionKernels *)kernels); + pmSubtractionVisualShowBasis(stamps); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", + psTimerClear("pmSubtractionCalculateEquation")); + + + return true; +} + +// 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); + +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) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + // Check inputs + int numKernels = kernels->num; // Number of kernel basis functions + int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for + int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + // An additional image is convolved + numSolution2 = numKernels * numSpatial; + numParams += numSolution2; + } + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + PS_ASSERT_PTR_NON_NULL(stamp, false); + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + + PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); + PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); + PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix, false); + PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, false); + } + + psString ds9name = NULL; // Filename for ds9 region file + static int ds9num = 0; // File number for ds9 region file + psStringAppend(&ds9name, "stamps_solution_%d.ds9", ds9num); + FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "solution stamps"); + psFree(ds9name); + ds9num++; + + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { + // Accumulate the least-squares matricies and vectors + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector + psVectorInit(sumVector, 0.0); + psImageInit(sumMatrix, 0.0); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + + int numStamps = 0; // Number of good stamps + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + + psVectorAppend(norms, stamp->norm); + + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); + numStamps++; + } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); + } + } + +#if 0 + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); +#endif + + psVector *solution = NULL; // Solution to equation! + solution = psVectorAlloc(numParams, PS_TYPE_F64); + psVectorInit(solution, 0); + +#if 0 + // Regular, straight-forward solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); +#else + { + // Solve normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(sumMatrix); + psFree(sumVector); + psFree(norms); + return false; + } + + // double normValue = 1.0; + double normValue = stats->robustMedian; + // double bgValue = 0.0; + + psFree(stats); + +#ifdef TESTING + fprintf(stderr, "Norm: %lf\n", normValue); +#endif + // Solve kernel components + for (int i = 0; i < numSolution1; i++) { + sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; + + sumMatrix->data.F64[i][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i] = 0.0; + } + sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; + sumMatrix->data.F64[bgIndex][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][bgIndex] = 0.0; + + sumMatrix->data.F64[normIndex][normIndex] = 1.0; + sumVector->data.F64[normIndex] = 0.0; + + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + + solution->data.F64[normIndex] = normValue; + } +# endif + +#if (1) + for (int i = 0; i < solution->n; i++) { + fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]); + } +#endif + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64); + psVectorInit(kernels->solution1, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + + psFree(norms); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int ix = 0; ix < kernels->solution1->n; ix++) { + if (!isfinite(kernels->solution1->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + + } else { + // Dual convolution solution + + // Accumulation of stamp matrices/vectors + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix, 0.0); + psVectorInit(sumVector, 0.0); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + + int numStamps = 0; // Number of good stamps + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + + psVectorAppend(norms, stamp->norm); + + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); + numStamps++; + } + } + +#if 0 + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); + psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); +#endif + +#if 1 + // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); + + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0); +#endif + + psVector *solution = NULL; // Solution to equation! + solution = psVectorAlloc(numParams, PS_TYPE_F64); + psVectorInit(solution, 0); + +#if 0 + // Regular, straight-forward solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); +#else + { + // Solve normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + +#if 0 + psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); + psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); + + normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; + normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; + normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; + + normVector->data.F64[0] = sumVector->data.F64[normIndex]; + normVector->data.F64[1] = sumVector->data.F64[bgIndex]; + + psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); + + double normValue = normSolution->data.F64[0]; + double bgValue = normSolution->data.F64[1]; + + psFree(normMatrix); + psFree(normVector); + psFree(normSolution); +#endif + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(sumMatrix); + psFree(sumVector); + psFree(norms); + return false; + } + + double normValue = stats->robustMedian; + + psFree(stats); + +#ifdef TESTING + fprintf(stderr, "Norm: %lf\n", normValue); +#endif + + // Solve kernel components + for (int i = 0; i < numSolution2; i++) { + sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; + sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; + + sumMatrix->data.F64[i][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i] = 0.0; + + sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; + } + sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; + sumMatrix->data.F64[bgIndex][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][bgIndex] = 0.0; + + sumMatrix->data.F64[normIndex][normIndex] = 1.0; + + sumVector->data.F64[normIndex] = 0.0; + +// save the matrix and vector after the NULLs have been set +#if 0 + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); + psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); +#endif + + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6); + // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 3e-4); + // psVectorCopy (solution, sumVector, PS_TYPE_F64); + // psMatrixGJSolve(sumMatrix, solution); + solution->data.F64[normIndex] = normValue; + } +#endif + + +#if (1) + for (int i = 0; i < solution->n; i++) { + fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); + } +#endif + + psFree(sumMatrix); + psFree(sumVector); + + psFree(norms); + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64); + psVectorInit (kernels->solution1, 0.0); + } + if (!kernels->solution2) { + kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64); + psVectorInit (kernels->solution2, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + for (int i = 0; i < numKernels * numSpatial; i++) { + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; + } + } + + + memcpy(kernels->solution1->data.F64, solution->data.F64, + numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], + numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + + psFree(solution); + + } + + if (ds9) { + fclose(ds9); + } + + if (psTraceGetLevel("psModules.imcombine") >= 7) { + for (int i = 0; i < kernels->solution1->n; i++) { + psTrace("psModules.imcombine", 7, "Solution 1 %d: %f\n", i, kernels->solution1->data.F64[i]); + } + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + for (int i = 0; i < kernels->solution2->n; i++) { + psTrace("psModules.imcombine", 7, "Solution 2 %d: %f\n", i, kernels->solution2->data.F64[i]); + } + } + } + + // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const + return true; +} + +// measure some useful stats on the stamp residuals: +// fResSigma : the residual stdev / total flux +// fResOuter : the residual fabs / total flux for R > 2 pix +// fResTotal : the residual fabs / total flux for R > 0 pix +bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { + + float sum = 0.0; + float peak = 0.0; + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm); + peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm)); + } + } + + // init counters + int npix = 0; + float dflux1 = 0.0; + float dflux2 = 0.0; + float dOuter = 0.0; + float dTotal = 0.0; + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + dflux1 += residual->kernel[y][x]; + dflux2 += PS_SQR(residual->kernel[y][x]); + dTotal += fabs(residual->kernel[y][x]); + if (hypot(x,y) > 2.0) { + dOuter += fabs(residual->kernel[y][x]); + } + npix ++; + } + } + float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix)); + if (!isfinite(sum)) return false; + if (!isfinite(peak)) return false; + if (!isfinite(dOuter)) return false; + if (!isfinite(dTotal)) return false; + + fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum); + psVectorAppend(fResSigma, sigma/sum); + psVectorAppend(fResOuter, dOuter/sum); + psVectorAppend(fResTotal, dTotal/sum); + return true; +} + +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); + + psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps + int footprint = stamps->footprint; // Half-size of stamps + long numPixels = PS_SQR(2 * footprint + 1); // Number of pixels in footprint + double devNorm = 1.0 / (double)numPixels; // Normalisation for deviations + int numKernels = kernels->num; // Number of kernels + + 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); + psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + + // we want to save the residual images for the 9 brightest stamps. + // identify the 9 brightest stamps + psVector *keepStamps = psVectorAlloc(stamps->num, PS_TYPE_S32); + psVectorInit (keepStamps, 0); + { + psVector *flux = psVectorAlloc(stamps->num, PS_TYPE_F32); + psVectorInit (flux, 0.0); + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (!isfinite(stamp->flux)) continue; + flux->data.F32[i] = stamp->flux; + } + + psVector *index = psVectorSortIndex(NULL, flux); + for (int i = 0; (i < stamps->num) && (i < 9); i++) { + int n = stamps->num - i - 1; + keepStamps->data.S32[index->data.S32[n]] = 1; + } + psFree (flux); + psFree (index); + + // this function is called multiple times in the iteration, but + // we only know after the interation is done if we will try again. + // therefore we must save the sample each time, and blow away the old one + // if it exists. + psFree (kernels->sampleStamps); + kernels->sampleStamps = psArrayAllocEmpty(9); + } + + psString log = psStringCopy("Deviations:\n"); // Log message with deviations + 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(psErrorCodeLast(), false, "Unable to generate image of kernel."); + return false; + } + double sumKernel = 0; // Sum of kernel, for normalising residual + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1; // Full size of kernel + for (int y = 0; y < fullSize; y++) { + for (int x = 0; x < fullSize; x++) { + sumKernel += image->data.F32[y][x]; + } + } + psFree(image); + devNorm = 1.0 / sumKernel / numPixels; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + source = stamp->image2; + convolutions = stamp->convolutions2; + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution + double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, + false); // Coefficient + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += convolution->kernel[y][x] * coefficient; + } + } + } + + // XXX visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i); + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x]; + } + } + + if (keepStamps->data.S32[i]) { + psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32); + psArrayAdd (kernels->sampleStamps, 9, sample); + psFree (sample); + } + + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint); + + } else { + // Dual convolution + psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image + psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image + psKernel *image1 = stamp->image1; // The first image + psKernel *image2 = stamp->image2; // The second image + + for (int j = 0; j < numKernels; j++) { + psKernel *conv1 = convolutions1->data[j]; // Convolution of first image + psKernel *conv2 = convolutions2->data[j]; // Convolution of second image + double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1 + double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2 + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1; + } + } + } + + // XXX visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i); + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x]; + } + } + if (keepStamps->data.S32[i]) { + psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32); + psArrayAdd (kernels->sampleStamps, 9, sample); + psFree (sample); + } + + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint); + } + + 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]); + 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)); + 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 + + } + + psFree(keepStamps); + + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log); + psFree(log); + + // calculate and report the normalization and background for the image center + { + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0); + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background + psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background); + + pmSubtractionVisualShowFit(norm); + pmSubtractionVisualPlotFit(kernels); + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psVectorStats (stats, fResSigma, NULL, NULL, 0); + kernels->fResSigmaMean = stats->robustMedian; + kernels->fResSigmaStdev = stats->robustStdev; + + psStatsInit (stats); + psVectorStats (stats, fResOuter, NULL, NULL, 0); + kernels->fResOuterMean = stats->robustMedian; + kernels->fResOuterStdev = stats->robustStdev; + + psStatsInit (stats); + psVectorStats (stats, fResTotal, NULL, NULL, 0); + kernels->fResTotalMean = stats->robustMedian; + kernels->fResTotalStdev = stats->robustStdev; + + // XXX save these values somewhere + psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f", + kernels->fResSigmaMean, kernels->fResSigmaStdev, + kernels->fResOuterMean, kernels->fResOuterStdev, + kernels->fResTotalMean, kernels->fResTotalStdev); + + psFree (fResSigma); + psFree (fResOuter); + psFree (fResTotal); + psFree (stats); + } + + psFree(residual); + psFree(polyValues); + + 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; +} + +// 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; +} + + +# if 0 + +#ifdef TESTING + psFitsWriteImageSimple("A.fits", sumMatrix, NULL); + psVectorWriteFile ("B.dat", sumVector); +#endif + +# define SVD_ANALYSIS 0 +# define COEFF_SIG 0.0 +# define SVD_TOL 0.0 + + // Use SVD to determine the kernel coeffs (and validate) + if (SVD_ANALYSIS) { + + // We have sumVector and sumMatrix. we are trying to solve the following equation: + // sumMatrix * x = sumVector. + + // we can use any standard matrix inversion to solve this. However, the basis + // functions in general have substantial correlation, so that the solution may be + // somewhat poorly determined or unstable. If not numerically ill-conditioned, the + // system of equations may be statistically ill-conditioned. Noise in the image + // will drive insignificant, but correlated, terms in the solution. To avoid these + // problems, we can use SVD to identify numerically unconstrained values and to + // avoid statistically badly determined value. + + // A = sumMatrix, B = sumVector + // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T + // x = V (1/w) (U^T B) + // \sigma_x = sqrt(diag(A^{-1})) + // solve for x and A^{-1} to get x & dx + // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 + // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 + + // If I use the SVD trick to re-condition the matrix, I need to break out the + // kernel and normalization terms from the background term. + // XXX is this true? or was this due to an error in the analysis? + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + + // now pull out the kernel elements into their own square matrix + psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); + psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { + if (ix == bgIndex) continue; + for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { + if (iy == bgIndex) continue; + kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; + ky++; + } + kernelVector->data.F64[kx] = sumVector->data.F64[ix]; + kx++; + } + + psImage *U = NULL; + psImage *V = NULL; + psVector *w = NULL; + if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { + psError(psErrorCodeLast(), false, "failed to perform SVD on sumMatrix\n"); + return NULL; + } + + // calculate A_inverse: + // Ainv = V * w * U^T + psImage *wUt = p_pmSubSolve_wUt (w, U); + psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); + psImage *Xvar = NULL; + psFree (wUt); + +# ifdef TESTING + // kernel terms: + for (int i = 0; i < w->n; i++) { + fprintf (stderr, "w: %f\n", w->data.F64[i]); + } +# endif + // loop over w adding in more and more of the values until chisquare is no longer + // dropping significantly. + // XXX this does not seem to work very well: we seem to need all terms even for + // simple cases... + + psVector *Xsvd = NULL; + { + psVector *Ax = NULL; + psVector *UtB = NULL; + psVector *wUtB = NULL; + + psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); + psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); + psVectorInit (wMask, 1); // start by masking everything + + double chiSquareLast = NAN; + int maxWeight = 0; + + double Axx, Bx, y2; + + // XXX this is an attempt to exclude insignificant modes. + // it was not successful with the ISIS kernel set: removing even + // the least significant mode leaves additional ringing / noise + // because the terms are so coupled. + for (int k = 0; false && (k < w->n); k++) { + + // unmask the k-th weight + wMask->data.U8[k] = 0; + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; + chiSquareLast = chiSquare; + + // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); + if (k && !maxWeight && (deltaChi < 1.0)) { + maxWeight = k; + } + } + + // keep all terms or we get extra ringing + maxWeight = w->n; + psVectorInit (wMask, 1); + for (int k = 0; k < maxWeight; k++) { + wMask->data.U8[k] = 0; + } + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); + + // re-calculate A^{-1} to get new variances: + // Ainv = V * w * U^T + // XXX since we keep all terms, this is identical to Ainv + psImage *wUt = p_pmSubSolve_wUt (wApply, U); + Xvar = p_pmSubSolve_VwUt (V, wUt); + psFree (wUt); + + psFree (Ax); + psFree (UtB); + psFree (wUtB); + psFree (wApply); + psFree (wMask); + } + + // copy the kernel solutions to the full solution vector: + solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { + if (ix == bgIndex) { + solution->data.F64[ix] = 0; + solutionErr->data.F64[ix] = 0.001; + continue; + } + solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); + solution->data.F64[ix] = Xsvd->data.F64[kx]; + kx++; + } + + psFree (kernelMatrix); + psFree (kernelVector); + + psFree (U); + psFree (V); + psFree (w); + + psFree (Ainv); + psFree (Xsvd); + } else { + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); + if (!luMatrix) { + psError(PM_ERR_DATA, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + + solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); + psFree(luMatrix); + psFree(permutation); + if (!solution) { + psError(PM_ERR_DATA, true, "Failed to solve the least-squares system.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + return NULL; + } + + // XXX LUD does not provide A^{-1}? fake the error for now + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + for (int ix = 0; ix < sumVector->n; ix++) { + solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; + } + } + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); + psVectorInit (kernels->solution1, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); + if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { + // XXX fprintf (stderr, "drop\n"); + kernels->solution1->data.F64[i] = 0.0; + } else { + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + } + // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); + // fprintf (stderr, "chi square Brute = %f\n", chiSquare); + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); +# endif + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int iy = 0; iy < stamp->matrix->numRows; iy++) { + for (int ix = 0; ix < stamp->matrix->numCols; ix++) { + if (!isfinite(stamp->matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + } + } + } + for (int ix = 0; ix < stamp->vector->n; ix++) { + if (!isfinite(stamp->vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + +#ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + +#ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } + { + psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); + psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); + psFree(inverse); + } + { + psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); + psImage *Xt = psMatrixTranspose(NULL, X); + psImage *XtX = psMatrixMultiply(NULL, Xt, X); + psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); + psFree(X); + psFree(Xt); + psFree(XtX); + } +#endif + Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 29543) @@ -183,4 +183,10 @@ } +# define CENTRAL_DELTA 0 + +# if (CENTRAL_DELTA) + +// XXX *** this code used the central pixel to force zero net flux, +// Alard actually uses kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0]) static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int uOrder, int vOrder, float fwhm, @@ -220,10 +226,10 @@ // Re-normalize // scale2D = 1.0 / fabs(sum); - scale2D = 1.0 / sqrt(sum2); + scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm); zeroNull = true; } else { // Odd functions: choose normalisation so that parameters have about the same strength as for even // functions, no subtraction of null pixel because the sum is already (near) zero - scale2D = 1.0 / sqrt(sum2); + scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm); zeroNull = false; } @@ -235,10 +241,10 @@ if (forceZeroNull) { // Force rescaling and subtraction of null pixel even though the order doesn't indicate it's even - scale2D = 1.0 / fabs(sum); + scale2D = 1.0 / fabs(sum) / PS_SQR(fwhm); zeroNull = true; } if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) { // Odd function - scale2D = 1.0 / sqrt(sum2); + scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm); } @@ -255,8 +261,8 @@ if (zeroNull) { // preCalc->kernel->kernel[0][0] -= 1.0; - preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2); - } - -#if 1 + preCalc->kernel->kernel[0][0] -= sum * scale2D; + } + +#if 0 { double Sum = 0.0; // Sum of kernel component @@ -287,4 +293,102 @@ return true; } + +# else /* CENTRAL_DELTA */ + +static bool zeroIsNormal = false; + +// this code uses kernel(0) to force zero flux, and is invalid for other kinds of normalizations +static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, + int index, int uOrder, int vOrder, float fwhm, + bool AlardLuptonStyle, bool forceZeroNull) +{ + // 1) for odd functions: no renormalization + // 2) for even functions, normalize the kernel to unity + // 3) for even functions & index > 0, subtract kernel(0) + + // Calculate moments + double sum = 0.0, sum2 = 0.0; // Sum of kernel component + float min = INFINITY, max = -INFINITY; // Minimum and maximum kernel value + + for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) { + for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) { + double value = preCalc->kernel->kernel[v][u]; + double value2 = PS_SQR(value); + sum += value; + sum2 += value2; + min = PS_MIN(value, min); + max = PS_MAX(value, max); + } + } + +#if 0 + fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max); +#endif + + float scale2D = 1.0; // Scaling for 2-D kernels + float scale1D = 1.0; // Scaling for 1-D kernels + + if (uOrder % 2 == 0 && vOrder % 2 == 0) { + + scale2D = 1.0 / sum; // Scaling for 2-D kernels + scale1D = sqrtf(scale2D); // Scaling for 1-D kernels + + for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) { + for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) { + preCalc->kernel->kernel[v][u] *= scale2D; + } + } + if (index == 0) { + zeroIsNormal = true; + } else { + psAssert(zeroIsNormal, "failed to normalize zero kernel first"); + pmSubtractionKernelPreCalc *zeroKernel = kernels->preCalc->data[0]; + psAssert(zeroKernel, "failed to supply zero kernel"); + for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) { + for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) { + preCalc->kernel->kernel[v][u] -= zeroKernel->kernel->kernel[v][u]; + } + } + } + + // XXX why do we bother renormlizing the 1D kernels? I don't think we use them again... + if (preCalc->xKernel) { + psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32)); + } + if (preCalc->yKernel) { + psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32)); + } + } + +#if 0 + { + double Sum = 0.0; // Sum of kernel component + double Sum2 = 0.0; // Sum of kernel component + float min = INFINITY, max = -INFINITY; // Minimum and maximum kernel value + for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) { + for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) { + double value = preCalc->kernel->kernel[v][u]; + Sum += value; + Sum2 += PS_SQR(value); + min = PS_MIN(preCalc->kernel->kernel[v][u], min); + max = PS_MAX(preCalc->kernel->kernel[v][u], max); + } + } + fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D); + } +#endif + + kernels->widths->data.F32[index] = fwhm; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + if (kernels->preCalc->data[index]) { + psFree(kernels->preCalc->data[index]); + } + kernels->preCalc->data[index] = preCalc; + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder); + + return true; +} +# endif /* Central Delta */ pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, @@ -603,10 +707,10 @@ kernels->sampleStamps = NULL; - kernels->fSigResMean = NAN; - kernels->fSigResStdev = NAN; - kernels->fMaxResMean = NAN; - kernels->fMaxResStdev = NAN; - kernels->fMinResMean = NAN; - kernels->fMinResStdev = NAN; + kernels->fResSigmaMean = NAN; + kernels->fResSigmaStdev = NAN; + kernels->fResOuterMean = NAN; + kernels->fResOuterStdev = NAN; + kernels->fResTotalMean = NAN; + kernels->fResTotalStdev = NAN; return kernels; Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 29543) @@ -51,10 +51,10 @@ float mean, rms; ///< Mean and RMS of chi^2 from stamps int numStamps; ///< Number of good stamps - float fSigResMean; ///< mean fractional stdev of residuals - float fSigResStdev; ///< stdev of fractional stdev of residuals - float fMaxResMean; ///< mean fractional positive swing in residuals - float fMaxResStdev; ///< stdev of fractional positive swing in residuals - float fMinResMean; ///< mean fractional negative swing in residuals - float fMinResStdev; ///< stdev of fractional negative swing in residuals + float fResSigmaMean; ///< mean fractional stdev of residuals + float fResSigmaStdev; ///< stdev of fractional stdev of residuals + float fResOuterMean; ///< mean fractional positive swing in residuals + float fResOuterStdev; ///< stdev of fractional positive swing in residuals + float fResTotalMean; ///< mean fractional negative swing in residuals + float fResTotalStdev; ///< stdev of fractional negative swing in residuals psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations } pmSubtractionKernels; Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 29543) @@ -28,10 +28,5 @@ static bool useFFT = true; // Do convolutions using FFT -# define SEPARATE 0 -# if (SEPARATE) -# define SUBMODE PM_SUBTRACTION_EQUATION_NORM -# else # define SUBMODE PM_SUBTRACTION_EQUATION_ALL -# endif //#define TESTING @@ -280,26 +275,26 @@ } -bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) { - - if (!readout) return true; - - psImage *image = readout->image; - psImage *mask = readout->mask; - psImage *variance = readout->variance; - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue; - bool valid = false; - valid = isfinite(image->data.F32[y][x]); - if (variance) { - valid &= isfinite(variance->data.F32[y][x]); - } - if (valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal; - } - } - - return true; -} +// bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) { +// +// if (!readout) return true; +// +// psImage *image = readout->image; +// psImage *mask = readout->mask; +// psImage *variance = readout->variance; +// for (int y = 0; y < image->numRows; y++) { +// for (int x = 0; x < image->numCols; x++) { +// if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue; +// bool valid = false; +// valid = isfinite(image->data.F32[y][x]); +// if (variance) { +// valid &= isfinite(variance->data.F32[y][x]); +// } +// if (valid) continue; +// mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal; +// } +// } +// +// return true; +// } bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, @@ -392,6 +387,7 @@ } - pmSubtractionMaskInvalid(ro1, maskVal); - pmSubtractionMaskInvalid(ro2, maskVal); + // XXX this is done before calling this function + // pmSubtractionMaskInvalid(ro1, maskVal); + // pmSubtractionMaskInvalid(ro2, maskVal); // General background subtraction, since this is done in pmSubtractionMatch @@ -565,6 +561,6 @@ memCheck("start"); - pmSubtractionMaskInvalid(ro1, maskVal); - pmSubtractionMaskInvalid(ro2, maskVal); + // pmSubtractionMaskInvalid(ro1, maskVal); + // pmSubtractionMaskInvalid(ro2, maskVal); psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels @@ -672,5 +668,4 @@ kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, inner, binning, ringsOrder, penalty, bounds, subMode); - // pmSubtractionVisualShowKernels(kernels); } @@ -678,27 +673,4 @@ if (subMode == PM_SUBTRACTION_MODE_UNSURE) { -#if 0 - // Get backgrounds - psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background - psVector *buffer = NULL;// Buffer for stats - if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) { - psError(PM_ERR_DATA, false, "Unable to measure background of image 1."); - psFree(bgStats); - psFree(buffer); - goto MATCH_ERROR; - } - float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1 - if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) { - psError(PM_ERR_DATA, false, "Unable to measure background of image 2."); - psFree(bgStats); - psFree(buffer); - goto MATCH_ERROR; - } - float bg2 = psStatsGetValue(bgStats, BG_STAT); // Background for image 2 - psFree(bgStats); - psFree(buffer); - - pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use -#endif pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); switch (newMode) { @@ -732,6 +704,21 @@ } - // XXX step 1: calculate normalization - psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); + // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue + psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); + if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + // step 0b : calculate the moments, pass along to the next steps via stamps->normValue + // XXX this step is now skipped -- delete + psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); + if (!pmSubtractionCalculateMoments(kernels, stamps)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + // step 1: generate the elements of the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); @@ -739,5 +726,6 @@ } - psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n"); + // step 2: solve the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Solving kernel equations...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); @@ -746,23 +734,4 @@ memCheck(" solve equation"); -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters - psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" calculate equation"); - - psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" solve equation"); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -770,5 +739,4 @@ goto MATCH_ERROR; } - memCheck(" calculate deviations"); @@ -787,5 +755,6 @@ // if we hit the max number of iterations and we have rejected stamps, re-solve if (numRejected > 0) { - // XXX step 1: calculate normalization + + // step 1: generate the elements of the matrix equation Ax = B psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { @@ -794,5 +763,5 @@ } - // solve normalization + // step 2: solve the matrix equation Ax = B psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { @@ -801,23 +770,4 @@ } -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters - psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - // solve kernel parameters - psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" solve equation"); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -837,5 +787,4 @@ goto MATCH_ERROR; } - memCheck("diag outputs"); @@ -1134,21 +1083,4 @@ } -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } -# endif - psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations @@ -1181,22 +1113,4 @@ } psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); - -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations @@ -1288,5 +1202,5 @@ bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, - float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax) + float scaleRef, float scaleMin, float scaleMax) { PS_ASSERT_PTR_NON_NULL(kernelSize, false); @@ -1294,6 +1208,4 @@ PS_ASSERT_VECTOR_NON_NULL(widths, false); PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false); - PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false); PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false); PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false); @@ -1301,9 +1213,13 @@ PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); -// float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference + float fwhm1; + float fwhm2; + + pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + psAssert(isfinite(fwhm1), "fwhm 1 not set"); + psAssert(isfinite(fwhm2), "fwhm 2 not set"); + + // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference float scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor - - // XXX save these values in a static for later use - pmSubtractionSetFWHMs(fwhm1, fwhm2); if (isfinite(scaleMin) && scale < scaleMin) { Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 29543) @@ -104,5 +104,4 @@ int *stampSize, ///< Half-size of the stamp (footprint) psVector *widths, ///< ISIS widths - float fwhm1, float fwhm2, ///< FWHMs for inputs float scaleRef, ///< Reference width for scaling float scaleMin, ///< Minimum scaling ratio, or NAN Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 29543) @@ -51,4 +51,5 @@ psFree(list->y); psFree(list->flux); + psFree(list->window); psFree(list->window1); psFree(list->window2); @@ -96,4 +97,6 @@ // fprintf (stderr, "xMin, xMax: %d, %d -> ", xMin, xMax); + // float xRaw = xMin; + // float yRaw = yMin; // Ensure we're not going to go outside the bounds of the image @@ -103,6 +106,12 @@ yMax = PS_MIN(numRows - border - 1, yMax); - if (xMax < xMin) return false; - if (yMax < yMin) return false; + if (xMax < xMin) { + // fprintf (stderr, "%f,%f : x-border\n", xRaw, yRaw); + return false; + } + if (yMax < yMin) { + // fprintf (stderr, "%f,%f : y-border\n", xRaw, yRaw); + return false; + } psAssert (xMin <= xMax, "x mismatch?"); @@ -118,4 +127,5 @@ if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) { + // fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw); return false; } @@ -133,4 +143,8 @@ } + // if (!found) { + // fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw); + // } + return found; } @@ -232,4 +246,6 @@ list->flux = NULL; list->normFrac = normFrac; + list->normValue = NAN; + list->window = NULL; list->window1 = NULL; list->window2 = NULL; @@ -256,4 +272,5 @@ out->y = NULL; out->flux = NULL; + out->window = psMemIncrRefCounter(in->window); out->window1 = psMemIncrRefCounter(in->window1); out->window2 = psMemIncrRefCounter(in->window2); @@ -331,4 +348,18 @@ stamp->vector = NULL; stamp->norm = NAN; + stamp->normI1 = NAN; + stamp->normI2 = NAN; + stamp->normSquare1 = NAN; + stamp->normSquare2 = NAN; + + stamp->MxxI1 = NULL; + stamp->MyyI1 = NULL; + stamp->MxxI2 = NULL; + stamp->MyyI2 = NULL; + + stamp->MxxI1raw = NAN; + stamp->MyyI1raw = NAN; + stamp->MxxI2raw = NAN; + stamp->MyyI2raw = NAN; return stamp; @@ -431,4 +462,5 @@ // Take stamps off the top of the (sorted) list for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) { + // fprintf (stderr, "%d : xList: %ld elements\n", i, xList->n); // Chop off the top of the list xList->n = j; @@ -459,9 +491,4 @@ yList->data.F32[j], yList->data.F32[j], numCols, numRows, border); #endif - if (0 && goodStamp) { - fprintf (stderr, "find: %6.1f < %6.1f < %6.1f, %6.1f < %6.1f %6.1f\n", - region->x0 + size, xStamp, region->x1 - size, - region->y0 + size, yStamp, region->y1 - size); - } } } else { @@ -656,4 +683,24 @@ psImageInit(stamps->window2->image, 0.0); + // Generate a weighting window based on the fwhms (20% larger than the largest) + { + float fwhm1, fwhm2; + + // XXX this is annoyingly hack-ish + pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + + float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35; + + psFree (stamps->window); + stamps->window = psKernelAlloc(-size, size, -size, size); + psImageInit(stamps->window->image, 0.0); + + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma)); + } + } + } + // generate normalizations for each stamp psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32); @@ -678,5 +725,5 @@ // storage vector for flux data - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); @@ -706,5 +753,5 @@ psAbort ("failed to generate stats"); } - float f1 = stats->robustMedian; + float f1 = stats->sampleMedian; psStatsInit (stats); @@ -712,5 +759,5 @@ psAbort ("failed to generate stats"); } - float f2 = stats->robustMedian; + float f2 = stats->sampleMedian; stamps->window1->kernel[y][x] = f1; @@ -728,13 +775,25 @@ } +#if 0 + { + psFits *fits = NULL; + fits = psFitsOpen ("window1.raw.fits", "w"); + psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL); + psFitsClose (fits); + fits = psFitsOpen ("window2.raw.fits", "w"); + psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL); + psFitsClose (fits); + } +#endif + psTrace("psModules.imcombine", 3, "Window total (1): %f, threshold: %f\n", sum1, (1.0 - stamps->normFrac) * sum1); psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2); -# if (0) +# if (1) // this block attempts to calculate the radius based on the first radial moment - bool done1 = false; - bool done2 = false; - double prior1 = 0.0; - double prior2 = 0.0; + double Sr1 = 0.0; + double Sr2 = 0.0; + double Sf1 = 0.0; + double Sf2 = 0.0; for (int y = -size; y <= size; y++) { for (int x = -size; x <= size; x++) { @@ -750,6 +809,8 @@ float R2 = Sr2 / Sf2; - stamps->normWindow1 = 2.5*R1; - stamps->normWindow1 = 2.5*R2; + stamps->normWindow1 = 2.0*R1; + stamps->normWindow2 = 2.0*R2; + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2); + # else // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent). @@ -778,5 +839,5 @@ // interpolate to the radius at which delta2 is normFrac: stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1); - fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1); + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1); done1 = true; } @@ -785,5 +846,5 @@ // interpolate to the radius at which delta2 is normFrac: stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2); - fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2); + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2); done2 = true; } @@ -833,8 +894,8 @@ { psFits *fits = NULL; - fits = psFitsOpen ("window1.fits", "w"); + fits = psFitsOpen ("window1.norm.fits", "w"); psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL); psFitsClose (fits); - fits = psFitsOpen ("window2.fits", "w"); + fits = psFitsOpen ("window2.norm.fits", "w"); psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL); psFitsClose (fits); @@ -906,6 +967,13 @@ float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull); - penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment - penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment + if (1) { + penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment + penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment + // penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment + // penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment + } else { + penalty1 = PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment + penalty2 = PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment + } } kernels->penalties1->data.F32[index] = kernels->penalty * penalty1; @@ -915,5 +983,5 @@ psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty"); - fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2); + // fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2); return true; @@ -942,6 +1010,6 @@ penalty *= 1.0 / sum2; - if (1) { - fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2); + if (0) { + // fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2); // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty); } @@ -983,6 +1051,7 @@ if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) { - psError(PM_ERR_PROG, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y); - return false; + psLogMsg("psModules.imcombine", 3, "Stamp %d (%d,%d) is within the region border.\n", i, x, y); + stamp->status = PM_SUBTRACTION_STAMP_NONE; + continue; } Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 29543) @@ -25,8 +25,11 @@ int footprint; ///< Half-size of stamps float normFrac; ///< Fraction of flux in window for normalisation window + float normValue; ///< calculated normalization + float normValue2; ///< calculated normalization psKernel *window1; ///< window function generated from ensemble of stamps (input 1) psKernel *window2; ///< window function generated from ensemble of stamps (input 2) - float normWindow1; ///< Size of window for measuring normalisation - float normWindow2; ///< Size of window for measuring normalisation + psKernel *window; ///< weighting window function (sigma = 1.1 * MAX(fwhm)) + float normWindow1; ///< Size of window for measuring normalisation + float normWindow2; ///< Size of window for measuring normalisation float sysErr; ///< Systematic error float skyErr; ///< increase effective readnoise @@ -85,5 +88,17 @@ psVector *vector; ///< Least-squares vector, or NULL double norm; ///< Normalisation difference + double normI1; ///< Sum(flux) for image 1 + double normI2; ///< Sum(flux) for image 2 + double normSquare1; ///< Sum(flux^2) for image 1 (used for penalty) + double normSquare2; ///< Sum(flux^2) for image 2 (used for penalty) pmSubtractionStampStatus status; ///< Status of stamp + psVector *MxxI1; ///< second moments of convolution images + psVector *MyyI1; ///< second moments of convolution images + psVector *MxxI2; ///< second moments of convolution images + psVector *MyyI2; ///< second moments of convolution images + double MxxI1raw; + double MyyI1raw; + double MxxI2raw; + double MyyI2raw; } pmSubtractionStamp; Index: trunk/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 29391) +++ trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 29543) @@ -146,4 +146,5 @@ } pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true); + psFree(canvas); pmVisualAskUser(&plotStamps); @@ -152,5 +153,5 @@ /** Plot the least-squares matrix of each stamp */ -bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) { +bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps) { if (!pmVisualTestLevel("ppsub.chisq", 1)) return true; @@ -209,4 +210,14 @@ pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true); + if (0) { + static int count = 0; + char filename[64]; + sprintf (filename, "chisq.%02d.fits", count); + count ++; + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, canvas32, 0, NULL); + psFitsClose (fits); + } + pmVisualAskUser(&plotLeastSquares); psFree(canvas); @@ -271,8 +282,9 @@ } } - fprintf (stderr, "kernel %d, sum %f\n", i, sum); + // fprintf (stderr, "kernel %d, sum %f\n", i, sum); } pmVisualScaleImage(kapa1, output, "Image", 0, true); + psFree(output); pmVisualAskUser(&plotImage); return true; @@ -286,4 +298,6 @@ return false; } + + // XXX clear the overlay(s) (red at least!) // get the kernel sizes @@ -297,8 +311,12 @@ if (!isfinite(stamp->flux)) continue; if (!stamp->convolutions1 && !stamp->convolutions2) continue; + // fprintf (stderr, "flux: %f, maxFlux: %f ", stamp->flux, maxFlux); if (!maxStamp) { maxFlux = stamp->flux; maxStamp = stamp; + // fprintf (stderr, "maxStamp %d\n", i); continue; + } else { + // fprintf (stderr, "\n"); } if (stamp->flux > maxFlux) { @@ -337,14 +355,22 @@ 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\n", i, sum); + // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); } pmVisualScaleImage(kapa2, output, "Image", 0, true); - } + + if (0) { + psFits *fits = psFitsOpen("basis.1.fits", "w"); + psFitsWriteImage(fits, NULL, output, 0, NULL); + psFitsClose(fits); + } + } if (maxStamp->convolutions2) { @@ -371,13 +397,22 @@ 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\n", i, sum); + // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); } pmVisualScaleImage(kapa2, output, "Image", 1, true); + + if (0) { + psFits *fits = psFitsOpen("basis.2.fits", "w"); + psFitsWriteImage(fits, NULL, output, 0, NULL); + psFitsClose(fits); + } + psFree(output); } @@ -676,4 +711,5 @@ psFree (x); psFree (y); + psFree (polyValues); pmVisualAskUser(NULL);