Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 21363) @@ -60,5 +60,5 @@ // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel - // This is required to keep the relative scaling between the image and the weight map + // This is required to keep the relative scaling between the image and the variance map psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32)); @@ -287,6 +287,6 @@ // Convolve an image using FFT -static void convolveWeightFFT(psImage *target,// Place the result in here - psImage *weight, // Weight map to convolve +static void convolveVarianceFFT(psImage *target,// Place the result in here + psImage *variance, // Variance map to convolve psImage *sys, // Systematic error image psImage *mask, // Mask image @@ -302,13 +302,13 @@ bool threaded = pmSubtractionThreaded(); // Are we running threaded? - psImage *subWeight = convolveSubsetAlloc(weight, border, threaded); // Weight map + psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once - psImage *convWeight = psImageConvolveFFT(NULL, subWeight, subMask, maskVal, kernel); // Convolved weight + psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys - convolveSubsetFree(weight, subWeight, threaded); + convolveSubsetFree(variance, subVariance, threaded); convolveSubsetFree(sys, subSys, threaded); convolveSubsetFree(mask, subMask, threaded); @@ -319,5 +319,5 @@ for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) { - target->data.F32[yTarget][xTarget] = convWeight->data.F32[ySource][xSource] + + target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] + convSys->data.F32[ySource][xSource]; } @@ -326,9 +326,9 @@ int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { - memcpy(&target->data.F32[yTarget][xMin], &convWeight->data.F32[ySource][size], numBytes); - } - } - - psFree(convWeight); + memcpy(&target->data.F32[yTarget][xMin], &convVariance->data.F32[ySource][size], numBytes); + } + } + + psFree(convVariance); psFree(convSys); @@ -361,10 +361,10 @@ // Convolve a region of an image static inline void convolveRegion(psImage *convImage, // Convolved image (output) - psImage *convWeight, // Convolved weight map (output), or NULL + psImage *convVariance, // Convolved variance map (output), or NULL psImage *convMask, // Convolve mask (output), or NULL psKernel **kernelImage, // Convolution kernel for the image - psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL + psKernel **kernelVariance, // Convolution kernel for the variance map, or NULL psImage *image, // Image to convolve - psImage *weight, // Weight map to convolve, or NULL + psImage *variance, // Variance map to convolve, or NULL psImage *sys, // Systematic error image, or NULL psImage *subMask, // Subtraction mask @@ -381,6 +381,6 @@ { *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual); - if (weight || subMask) { - *kernelWeight = varianceKernel(*kernelWeight, *kernelImage); + if (variance || subMask) { + *kernelVariance = varianceKernel(*kernelVariance, *kernelImage); } @@ -398,17 +398,17 @@ } - // Convolve the image and weight + // Convolve the image and variance if (useFFT) { // Use Fast Fourier Transform to do the convolution // This provides a big speed-up for large kernels convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size); - if (weight) { - convolveWeightFFT(convWeight, weight, sys, subMask, subBad, *kernelWeight, region, kernels->size); + if (variance) { + convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size); } } else { // XXX Direct convolution doesn't account for bad pixels yet convolveDirect(convImage, image, *kernelImage, region, background, kernels->size); - if (weight) { - convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size); + if (variance) { + convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size); } } @@ -885,6 +885,6 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->weight); - stamp->image1 = stamp->image2 = stamp->weight = NULL; + psFree(stamp->variance); + stamp->image1 = stamp->image2 = stamp->variance = NULL; psFree(stamp->matrix1); psFree(stamp->matrix2); @@ -915,5 +915,5 @@ } -psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +psKernel *pmSubtractionKernel(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); @@ -922,12 +922,19 @@ PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - // Precalulate polynomial values - psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); - - // The appropriate kernel - psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); - + psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel psFree(polyValues); + return kernel; +} + +psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); + + psKernel *kernel = pmSubtractionKernel(kernels, x, y, wantDual); // Convolution kernel psImage *image = psMemIncrRefCounter(kernel->image); // Image of the kernel psFree(kernel); @@ -989,5 +996,5 @@ -// XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data +// XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image int x0, int y0, // Offsets for image @@ -1010,5 +1017,5 @@ psKernel *kernelImage = NULL; // Kernel for the images - psKernel *kernelWeight = NULL; // Kernel for the weight maps + psKernel *kernelVariance = NULL; // Kernel for the variance maps // Only generate polynomial values every kernel footprint, since we have already assumed @@ -1020,16 +1027,16 @@ if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight, - ro1->image, ro1->weight, sys1, subMask, kernels, polyValues, background, *region, + convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance, + ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT, false); } if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - convolveRegion(out2->image, out2->weight, convMask, &kernelImage, &kernelWeight, - ro2->image, ro2->weight, sys2, subMask, kernels, polyValues, background, *region, + convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance, + ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL); } psFree(kernelImage); - psFree(kernelWeight); + psFree(kernelVariance); psFree(polyValues); @@ -1148,12 +1155,12 @@ // XXX } } - if (ro1->weight) { - if (!out1->weight) { - out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + if (ro1->variance) { + if (!out1->variance) { + out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // XXX if (threaded) { - // XXX psMutexInit(out1->weight); + // XXX psMutexInit(out1->variance); // XXX } } - psImageInit(out1->weight, 0.0); + psImageInit(out1->variance, 0.0); } } @@ -1165,12 +1172,12 @@ // XXX } } - if (ro2->weight) { - if (!out2->weight) { - out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + if (ro2->variance) { + if (!out2->variance) { + out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // XXX if (threaded) { - // XXX psMutexInit(out2->weight); + // XXX psMutexInit(out2->variance); // XXX } } - psImageInit(out2->weight, 0.0); + psImageInit(out2->variance, 0.0); } } @@ -1232,5 +1239,5 @@ psImage *polyValues = NULL; // Pre-calculated polynomial values psKernel *kernelImage = NULL; // Kernel for the images - psKernel *kernelWeight = NULL; // Kernel for the weight maps + psKernel *kernelVariance = NULL; // Kernel for the variance maps #endif @@ -1299,7 +1306,7 @@ psFree(job); } else { - subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, - sys1, sys2, subMask, maskBad, maskPoor, poorFrac, subRegion, - kernels, doBG, useFFT); + subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2, + subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG, + useFFT); } psFree(subRegion); @@ -1334,9 +1341,23 @@ // XXX } } - psImageConvolveSetThreads(oldThreads); psFree(sys1); psFree(sys2); + + // Calculate covariances + // This can take a while, so we only do it for a single instance + // XXX psImageCovarianceCalculate could be multithreaded + if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance); + psFree(kernel); + } + if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, + kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel + out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance); + psFree(kernel); + } // Copy anything that wasn't convolved @@ -1345,6 +1366,7 @@ if (out2) { out2->image = psMemIncrRefCounter(ro2->image); - out2->weight = psMemIncrRefCounter(ro2->weight); + out2->variance = psMemIncrRefCounter(ro2->variance); out2->mask = psMemIncrRefCounter(ro2->mask); + out2->covariance = psMemIncrRefCounter(ro2->covariance); } break; @@ -1352,6 +1374,7 @@ if (out1) { out1->image = psMemIncrRefCounter(ro1->image); - out1->weight = psMemIncrRefCounter(ro1->weight); + out1->variance = psMemIncrRefCounter(ro1->variance); out1->mask = psMemIncrRefCounter(ro1->mask); + out1->covariance = psMemIncrRefCounter(ro1->covariance); } break;