Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 26584) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 26893) @@ -21,4 +21,5 @@ #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" +#include "pmSubtractionVisual.h" #include "pmSubtractionThreads.h" @@ -63,21 +64,19 @@ } -// Contribute to an image of the solved kernel component for ISIS -static void solvedKernelISIS(psKernel *kernel, // Kernel, updated - const pmSubtractionKernels *kernels, // Kernel basis functions - float value, // Normalisation value for basis function - int index // Index of basis function of interest +// Contribute to an image of the solved kernel component using the preCalculated image +static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated + const pmSubtractionKernels *kernels, // Kernel basis functions + float value, // Normalisation value for basis function + int index // Index of basis function of interest ) { int size = kernels->size; // Kernel half-size - psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values #if 0 - psVector *xKernel = preCalc->data[0]; // Kernel in x - psVector *yKernel = preCalc->data[1]; // Kernel in y // Iterating over the kernel for (int y = 0, v = -size; v <= size; y++, v++) { - float yValue = value * yKernel->data.F32[y]; + float yValue = value * preCalc->yKernel->data.F32[y]; for (int x = 0, u = -size; u <= size; x++, u++) { - kernel->kernel[v][u] += yValue * xKernel->data.F32[x]; + kernel->kernel[v][u] += yValue * preCalc->xKernel->data.F32[x]; } } @@ -87,8 +86,7 @@ } #else - psKernel *k = preCalc->data[2]; // Kernel image for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - kernel->kernel[v][u] += value * k->kernel[v][u]; + kernel->kernel[v][u] += value * preCalc->kernel->kernel[v][u]; } } @@ -119,4 +117,9 @@ for (int i = 0; i < numKernels; i++) { double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value + if (wantDual) { + // The model is built with the dual convolution terms added, so to produce zero residual the + // equation results in negative coefficients which we must undo. + value *= -1.0; + } switch (kernels->type) { @@ -149,5 +152,5 @@ case PM_SUBTRACTION_KERNEL_GUNK: { if (i < kernels->inner) { - solvedKernelISIS(kernel, kernels, value, i); + solvedKernelPreCalc(kernel, kernels, value, i); } else { // Using delta function @@ -159,18 +162,19 @@ break; } - case PM_SUBTRACTION_KERNEL_ISIS: { - solvedKernelISIS(kernel, kernels, value, i); + case PM_SUBTRACTION_KERNEL_ISIS: + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + case PM_SUBTRACTION_KERNEL_HERM: + case PM_SUBTRACTION_KERNEL_DECONV_HERM: { + solvedKernelPreCalc(kernel, kernels, value, i); break; } case PM_SUBTRACTION_KERNEL_RINGS: { - psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data - psVector *uCoords = preCalc->data[0]; // u coordinates - psVector *vCoords = preCalc->data[1]; // v coordinates - psVector *poly = preCalc->data[2]; // Polynomial values - int num = uCoords->n; // Number of pixels + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Precalculated kernels + int num = preCalc->uCoords->n; // Number of pixels for (int j = 0; j < num; j++) { - int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates - kernel->kernel[v][u] += poly->data.F32[j] * value; + int u = preCalc->uCoords->data.S32[j]; + int v = preCalc->vCoords->data.S32[j]; // Kernel coordinates + kernel->kernel[v][u] += preCalc->poly->data.F32[j] * value; } // Photometric scaling is built into the kernel --- no subtraction! @@ -419,5 +423,8 @@ *target |= maskBad; } else if (*source & subConvPoor) { + *target &= ~maskBad; *target |= maskPoor; + } else { + *target &= ~maskBad & ~maskPoor; } } @@ -454,16 +461,14 @@ #endif -// Convolve a stamp using an ISIS kernel basis function -static psKernel *convolveStampISIS(const psKernel *image, // Image to convolve - const pmSubtractionKernels *kernels, // Kernel basis functions - int index, // Index of basis function of interest - int footprint // Half-size of stamp +// Convolve a stamp using a pre-calculated kernel basis function +static psKernel *convolveStampPreCalc(const psKernel *image, // Image to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int index, // Index of basis function of interest + int footprint // Half-size of stamp ) { - psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data #if 0 // Convolving using separable convolution - psVector *xKernel = preCalc->data[0]; // Kernel in x - psVector *yKernel = preCalc->data[1]; // Kernel in y int size = kernels->size; // Size of kernel @@ -477,5 +482,5 @@ float value = 0.0; // Value of convolved pixel int uMin = x - size, uMax = x + size; // Range for u - psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values + psF32 *xKernelData = &preCalc->xKernel->data.F32[xKernel->n - 1]; // Kernel values psF32 *imageData = &image->kernel[y][uMin]; // Image values for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) { @@ -492,5 +497,5 @@ float value = 0.0; // Value of convolved pixel int vMin = y - size, vMax = y + size; // Range for v - psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values + psF32 *yKernelData = &preCalc->yKernel->data.F32[yKernel->n - 1]; // Kernel values psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way! for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) { @@ -509,5 +514,5 @@ #else // Convolving using precalculated kernel - return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]); + return p_pmSubtractionConvolveStampPrecalc(image, preCalc->kernel); #endif } @@ -525,4 +530,7 @@ int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version + + // pmSubtractionVisualShowSubtraction(image->image, kernel->image, conv); + psFree(conv); return convolved; @@ -569,19 +577,21 @@ } +void p_pmSubtractionPolynomialNormCoords(float *xOut, float *yOut, float xIn, float yIn, + int xMin, int xMax, int yMin, int yMax) +{ + float xNormSize = xMax - xMin, yNormSize = yMax - yMin; // Size to use for normalisation + *xOut = 2.0 * (float)(xIn - xMin - xNormSize/2.0) / xNormSize; + *yOut = 2.0 * (float)(yIn - yMin - yNormSize/2.0) / yNormSize; + return; +} + psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels, - int numCols, int numRows, int x, int y) + int x, int y) { assert(kernels); - assert(numCols > 0 && numRows > 0); - - // Size to use when calculating normalised coordinates (different from actual size when convolving - // subimage) - int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols); - int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows); - - // Normalised coordinates - float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize; - float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize; - + + float xNorm, yNorm; // Normalised coordinates + p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, x, y, + kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax); return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm); } @@ -664,5 +674,5 @@ if (index < kernels->inner) { // Photometric scaling is already built in to the precalculated kernel - return convolveStampISIS(image, kernels, index, footprint); + return convolveStampPreCalc(image, kernels, index, footprint); } // Using delta function @@ -673,17 +683,18 @@ return convolved; } - case PM_SUBTRACTION_KERNEL_ISIS: { - return convolveStampISIS(image, kernels, index, footprint); - } + case PM_SUBTRACTION_KERNEL_ISIS: + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + case PM_SUBTRACTION_KERNEL_HERM: + case PM_SUBTRACTION_KERNEL_DECONV_HERM: { + return convolveStampPreCalc(image, kernels, index, footprint); + } case PM_SUBTRACTION_KERNEL_RINGS: { - psKernel *convolved = psKernelAlloc(-footprint, footprint, - -footprint, footprint); // Convolved image - psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data - psVector *uCoords = preCalc->data[0]; // u coordinates - psVector *vCoords = preCalc->data[1]; // v coordinates - psVector *poly = preCalc->data[2]; // Polynomial values - int num = uCoords->n; // Number of pixels - psS32 *uData = uCoords->data.S32, *vData = vCoords->data.S32; // Dereference u,v coordinates - psF32 *polyData = poly->data.F32; // Dereference polynomial values + psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data + + int num = preCalc->uCoords->n; // Number of pixels + psS32 *uData = preCalc->uCoords->data.S32; // Dereference v coordinate + psS32 *vData = preCalc->vCoords->data.S32; // Dereference u coordinate + psF32 *polyData = preCalc->poly->data.F32; // Dereference polynomial values psF32 **imageData = image->kernel; // Dereference image psF32 **convData = convolved->kernel; // Dereference convolved image @@ -765,6 +776,4 @@ - - int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, const psVector *deviations, psImage *subMask, float sigmaRej) @@ -860,4 +869,6 @@ int numGood = 0; // Number of good stamps double newMean = 0.0; // New mean + psString log = NULL; // Log message + psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit); for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -872,5 +883,8 @@ // Mask out the stamp in the image so you it's not found again psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, - (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); + psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i, + (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), + fabsf(deviations->data.F32[i] - mean)); numRejected++; for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { @@ -895,11 +909,8 @@ psFree(stamp->weight); stamp->image1 = stamp->image2 = stamp->weight = NULL; - psFree(stamp->matrix1); - psFree(stamp->matrix2); - psFree(stamp->matrixX); - stamp->matrix1 = stamp->matrix2 = stamp->matrixX = NULL; - psFree(stamp->vector1); - psFree(stamp->vector2); - stamp->vector1 = stamp->vector2 = NULL; + psFree(stamp->matrix); + stamp->matrix = NULL; + psFree(stamp->vector); + stamp->vector = NULL; } else { numGood++; @@ -910,4 +921,10 @@ } newMean /= numGood; + + if (numRejected == 0) { + psStringAppend(&log, "\n"); + } + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log); + psFree(log); if (ds9) { @@ -995,5 +1012,5 @@ psVector *backup = psVectorCopy(NULL, solution, PS_TYPE_F64); // Backup version - int num = wantDual ? solution->n - 1 : solution->n; // Number of kernel basis functions + int num = kernels->num; // Number of kernel basis functions psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial @@ -1060,7 +1077,6 @@ // Only generate polynomial values every kernel footprint, since we have already assumed // (with the stamps) that it does not vary rapidly on this scale. - psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows, - xMin + x0 + size + 1, - yMin + y0 + size + 1); + psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + x0 + size + 1, + yMin + y0 + size + 1); // Polynomial float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term @@ -1137,5 +1153,5 @@ bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, psImage *subMask, int stride, psImageMaskType maskBad, psImageMaskType maskPoor, - float poorFrac, float kernelError, const psRegion *region, + float poorFrac, float kernelError, float covarFrac, const psRegion *region, const pmSubtractionKernels *kernels, bool doBG, bool useFFT) { @@ -1146,4 +1162,5 @@ PM_ASSERT_READOUT_NON_NULL(ro1, false); PM_ASSERT_READOUT_IMAGE(ro1, false); + PM_ASSERT_READOUT_IMAGE(out1, false); numCols = ro1->image->numCols; numRows = ro1->image->numRows; @@ -1155,4 +1172,5 @@ PM_ASSERT_READOUT_NON_NULL(ro2, false); PM_ASSERT_READOUT_IMAGE(ro2, false); + PM_ASSERT_READOUT_IMAGE(out2, false); if (numCols == 0 && numRows == 0) { numCols = ro2->image->numCols; @@ -1177,4 +1195,6 @@ PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false); PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(kernelError, 1.0, false); + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(covarFrac, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN(covarFrac, 1.0, false); if (region && psRegionIsNaN(*region)) { psString string = psRegionToString(*region); @@ -1188,41 +1208,10 @@ bool threaded = pmSubtractionThreaded(); // Running threaded? - // Outputs - if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out1->image) { - out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - if (ro1->variance) { - if (!out1->variance) { - out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - psImageInit(out1->variance, 0.0); - } - } - if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out2->image) { - out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - if (ro2->variance) { - if (!out2->variance) { - out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - psImageInit(out2->variance, 0.0); - } - } psImage *convMask = NULL; // Convolved mask image (common to inputs 1 and 2) if (subMask) { if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out1->mask) { - out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - } - psImageInit(out1->mask, 0); convMask = out1->mask; } if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out2->mask) { - out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - } - psImageInit(out2->mask, 0); if (!convMask) { convMask = out2->mask; @@ -1244,6 +1233,6 @@ // Get region for convolution: [xMin:xMax,yMin:yMax] - int xMin = size, xMax = numCols - size; - int yMin = size, yMax = numRows - size; + int xMin = kernels->xMin + size, xMax = kernels->xMax - size; + int yMin = kernels->yMin + size, yMax = kernels->yMax - size; if (region) { xMin = PS_MAX(region->x0, xMin); @@ -1336,8 +1325,10 @@ // Calculate covariances - // This can take a while, so we only do it for a single instance - // XXX psImageCovarianceCalculate could be multithreaded + // This can be fairly involved, so we only do it for a single instance + // Enable threads for covariance calculation, since we're not threading on top of it. + oldThreads = psImageCovarianceSetThreads(true); if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + psKernelTruncate(kernel, covarFrac); out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance); psFree(kernel); @@ -1346,7 +1337,9 @@ psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel + psKernelTruncate(kernel, covarFrac); out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance); psFree(kernel); } + psImageCovarianceSetThreads(oldThreads); // Copy anything that wasn't convolved @@ -1388,5 +1381,4 @@ } - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve image: %f sec", psTimerClear("pmSubtractionConvolve"));