Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 26034) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 26035) @@ -29,5 +29,5 @@ #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time #define MIN_SAMPLE_STATS 7 // Minimum number to use sample statistics; otherwise use quartiles -#define USE_SYS_ERR // Use systematic error image? +#define USE_KERNEL_ERR // Use kernel error image? ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -266,5 +266,5 @@ static void convolveVarianceFFT(psImage *target,// Place the result in here psImage *variance, // Variance map to convolve - psImage *sys, // Systematic error image + psImage *kernelErr, // Kernel error image psImage *mask, // Mask image psImageMaskType maskVal, // Value to mask @@ -278,22 +278,22 @@ psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map - psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image + psImage *subKE = kernelErr ? psImageSubset(kernelErr, border) : NULL; // Kernel error image psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance - psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys + psImage *convKE = subKE ? psImageConvolveFFT(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE psFree(subVariance); - psFree(subSys); + psFree(subKE); psFree(subMask); // Now, we have to stick it in where it belongs int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region - if (convSys) { + if (convKE) { 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] = convVariance->data.F32[ySource][xSource] + - convSys->data.F32[ySource][xSource]; + convKE->data.F32[ySource][xSource]; } } @@ -306,5 +306,5 @@ psFree(convVariance); - psFree(convSys); + psFree(convKE); return; @@ -342,5 +342,5 @@ psImage *image, // Image to convolve psImage *variance, // Variance map to convolve, or NULL - psImage *sys, // Systematic error image, or NULL + psImage *kernelErr, // Kernel error image, or NULL psImage *subMask, // Subtraction mask const pmSubtractionKernels *kernels, // Kernels @@ -379,5 +379,6 @@ convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size); if (variance) { - convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size); + convolveVarianceFFT(convVariance, variance, kernelErr, subMask, subBad, *kernelVariance, + region, kernels->size); } } else { @@ -429,25 +430,25 @@ } -#ifdef USE_SYS_ERR -// Generate an image that can be used to track systematic errors -static psImage *subtractionSysErrImage(const psImage *image, // Image from which to make sys err image - float sysError // Relative systematic error - ) -{ - if (!isfinite(sysError) || sysError == 0.0) { +#ifdef USE_KERNEL_ERR +// Generate an image that can be used to track systematic errors in the kernel +static psImage *subtractionKernelErrImage(const psImage *image, // Image from which to make kernel error image + float kernelError // Relative systematic error in kernel + ) +{ + if (!isfinite(kernelError) || kernelError == 0.0) { return NULL; } int numCols = image->numCols, numRows = image->numRows; // Size of image - psImage *sys = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Systematic error image - - float sysError2 = PS_SQR(sysError); // Square of the systematic error + psImage *kernelErr = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Kernel error image + + float kernelError2 = PS_SQR(kernelError); // Square of the kernel error for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - sys->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * sysError2; - } - } - - return sys; + kernelErr->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * kernelError2; + } + } + + return kernelErr; } #endif @@ -892,6 +893,6 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->variance); - stamp->image1 = stamp->image2 = stamp->variance = NULL; + psFree(stamp->weight); + stamp->image1 = stamp->image2 = stamp->weight = NULL; psFree(stamp->matrix1); psFree(stamp->matrix2); @@ -1040,5 +1041,5 @@ psImage *convMask, // Output convolved mask const pmReadout *ro1, const pmReadout *ro2, // Input readouts - psImage *sys1, psImage *sys2, // Systematic error images + psImage *kernelErr1, psImage *kernelErr2, // Kernel error images psImage *subMask, // Input subtraction mask psImageMaskType maskBad, // Mask value to give bad pixels @@ -1066,11 +1067,12 @@ if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance, - ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region, - maskBad, maskPoor, poorFrac, useFFT, false); + ro1->image, ro1->variance, kernelErr1, 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->variance, convMask, &kernelImage, &kernelVariance, - ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region, - maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL); + ro2->image, ro2->variance, kernelErr2, subMask, kernels, polyValues, background, + *region, maskBad, maskPoor, poorFrac, useFFT, + kernels->mode == PM_SUBTRACTION_MODE_DUAL); } @@ -1117,6 +1119,6 @@ const pmReadout *ro1 = args->data[7]; // Input readout 1 const pmReadout *ro2 = args->data[8]; // Input readout 2 - psImage *sys1 = args->data[9]; // Systematic error image 1 - psImage *sys2 = args->data[10]; // Systematic error image 2 + psImage *kernelErr1 = args->data[9]; // Kernel error image 1 + psImage *kernelErr2 = args->data[10]; // Kernel error image 2 psImage *subMask = args->data[11]; // Subtraction mask psImageMaskType maskBad = PS_SCALAR_VALUE(args->data[12], PS_TYPE_IMAGE_MASK_DATA); // Output mask value for bad pixels @@ -1128,11 +1130,12 @@ bool useFFT = PS_SCALAR_VALUE(args->data[18], U8); // Use FFT for convolution? - return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2, - subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT); + return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, kernelErr1, + kernelErr2, subMask, maskBad, maskPoor, poorFrac, region, kernels, + doBG, useFFT); } bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, psImage *subMask, int stride, psImageMaskType maskBad, psImageMaskType maskPoor, - float poorFrac, float sysError, const psRegion *region, + float poorFrac, float kernelError, const psRegion *region, const pmSubtractionKernels *kernels, bool doBG, bool useFFT) { @@ -1172,6 +1175,6 @@ PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(poorFrac, 0.0, false); PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, false); - PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false); - PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(sysError, 1.0, false); + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(kernelError, 1.0, false); if (region && psRegionIsNaN(*region)) { psString string = psRegionToString(*region); @@ -1232,11 +1235,11 @@ } - psImage *sys1 = NULL, *sys2 = NULL; // Systematic error images -#ifdef USE_SYS_ERR + psImage *kernelErr1 = NULL, *kernelErr2 = NULL; // Kernel error images +#ifdef USE_KERNEL_ERR if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - sys1 = subtractionSysErrImage(ro1->image, sysError); + kernelErr1 = subtractionKernelErrImage(ro1->image, kernelError); } if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - sys2 = subtractionSysErrImage(ro2->image, sysError); + kernelErr2 = subtractionKernelErrImage(ro2->image, kernelError); } #endif @@ -1289,6 +1292,6 @@ psArrayAdd(args, 1, (pmReadout*)ro1); // Casting away const psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const - psArrayAdd(args, 1, sys1); - psArrayAdd(args, 1, sys2); + psArrayAdd(args, 1, kernelErr1); + psArrayAdd(args, 1, kernelErr2); psArrayAdd(args, 1, subMask); PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_IMAGE_MASK); @@ -1306,7 +1309,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, + kernelErr1, kernelErr2, subMask, maskBad, maskPoor, poorFrac, + subRegion, kernels, doBG, useFFT); } psFree(subRegion); @@ -1333,6 +1336,6 @@ psImageConvolveSetThreads(oldThreads); - psFree(sys1); - psFree(sys2); + psFree(kernelErr1); + psFree(kernelErr2); // Calculate covariances Index: /trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 26034) +++ /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 26035) @@ -17,5 +17,5 @@ //#define TESTING // TESTING output for debugging; may not work with threads! -#define USE_VARIANCE // Include variance in equation? +#define USE_WEIGHT // Include weight (1/variance) in equation? @@ -29,5 +29,5 @@ const psKernel *input, // Input image (target) const psKernel *reference, // Reference image (convolution source) - const psKernel *variance, // Variance image + const psKernel *weight, // Weight image const psArray *convolutions, // Convolutions for each kernel const pmSubtractionKernels *kernels, // Kernels @@ -74,6 +74,6 @@ for (int x = - footprint; x <= footprint; x++) { double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; -#ifdef USE_VARIANCE - cc /= variance->kernel[y][x]; +#ifdef USE_WEIGHT + cc *= weight->kernel[y][x]; #endif sumCC += cc; @@ -102,9 +102,9 @@ double rc = ref * conv; double c = conv; -#ifdef USE_VARIANCE - float varVal = 1.0 / variance->kernel[y][x]; - ic *= varVal; - rc *= varVal; - c *= varVal; +#ifdef USE_WEIGHT + float wtVal = weight->kernel[y][x]; + ic *= wtVal; + rc *= wtVal; + c *= wtVal; #endif sumIC += ic; @@ -137,11 +137,11 @@ double rr = PS_SQR(ref); double one = 1.0; -#ifdef USE_VARIANCE - float varVal = 1.0 / variance->kernel[y][x]; - rr *= varVal; - ir *= varVal; - in *= varVal; - ref *= varVal; - one *= varVal; +#ifdef USE_WEIGHT + float wtVal = weight->kernel[y][x]; + rr *= wtVal; + ir *= wtVal; + in *= wtVal; + ref *= wtVal; + one *= wtVal; #endif sumRR += rr; @@ -169,5 +169,5 @@ const psKernel *image1, // Image 1 const psKernel *image2, // Image 2 - const psKernel *variance, // Variance image + const psKernel *weight, // Weight image const psArray *convolutions1, // Convolutions of image 1 for each kernel const psArray *convolutions2, // Convolutions of image 2 for each kernel @@ -227,8 +227,9 @@ double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; -#ifdef USE_VARIANCE - aa /= variance->kernel[y][x]; - bb /= variance->kernel[y][x]; - ab /= variance->kernel[y][x]; +#ifdef USE_WEIGHT + float wtVal = weight->kernel[y][x]; + aa *= wtVal; + bb *= wtVal; + ab *= wtVal; #endif sumAA += aa; @@ -258,6 +259,6 @@ for (int x = - footprint; x <= footprint; x++) { double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; -#ifdef USE_VARIANCE - ab /= variance->kernel[y][x]; +#ifdef USE_WEIGHT + ab *= weight->kernel[y][x]; #endif sumAB += ab; @@ -295,14 +296,14 @@ double i1i2 = i1 * i2; -#ifdef USE_VARIANCE - float varVal = 1.0 / variance->kernel[y][x]; - ai2 *= varVal; - bi2 *= varVal; - ai1 *= varVal; - bi1 *= varVal; - i1i2 *= varVal; - a *= varVal; - b *= varVal; - i2 *= varVal; +#ifdef USE_WEIGHT + float wtVal = weight->kernel[y][x]; + ai2 *= wtVal; + bi2 *= wtVal; + ai1 *= wtVal; + bi1 *= wtVal; + i1i2 *= wtVal; + a *= wtVal; + b *= wtVal; + i2 *= wtVal; #endif @@ -351,11 +352,11 @@ double i1i2 = i1 * i2; -#ifdef USE_VARIANCE - float varVal = 1.0 / variance->kernel[y][x]; - i1 *= varVal; - i1i1 *= varVal; - one *= varVal; - i2 *= varVal; - i1i2 *= varVal; +#ifdef USE_WEIGHT + float wtVal = weight->kernel[y][x]; + i1 *= wtVal; + i1i1 *= wtVal; + one *= wtVal; + i2 *= wtVal; + i1i2 *= wtVal; #endif @@ -619,10 +620,10 @@ case PM_SUBTRACTION_MODE_1: status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, - stamp->variance, stamp->convolutions1, kernels, polyValues, + stamp->weight, stamp->convolutions1, kernels, polyValues, footprint); break; case PM_SUBTRACTION_MODE_2: status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2, - stamp->variance, stamp->convolutions2, kernels, polyValues, + stamp->weight, stamp->convolutions2, kernels, polyValues, footprint); break; @@ -639,5 +640,5 @@ #endif status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, - stamp->matrixX, stamp->image1, stamp->image2, stamp->variance, + stamp->matrixX, stamp->image1, stamp->image2, stamp->weight, stamp->convolutions1, stamp->convolutions2, kernels, polyValues, footprint); @@ -1163,5 +1164,5 @@ // Calculate residuals - psKernel *variance = stamp->variance; // Variance postage stamp + psKernel *weight = stamp->weight; // Weight postage stamp psImageInit(residual->image, 0.0); if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { @@ -1247,5 +1248,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x]; + double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x]; deviation += dev; #ifdef TESTING @@ -1290,10 +1291,10 @@ psFitsClose(fits); } - if (stamp->variance) { + if (stamp->weight) { psString filename = NULL; - psStringAppend(&filename, "stamp_variance_%03d.fits", i); + psStringAppend(&filename, "stamp_weight_%03d.fits", i); psFits *fits = psFitsOpen(filename, "w"); psFree(filename); - psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL); + psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); psFitsClose(fits); } Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 26034) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 26035) @@ -68,4 +68,5 @@ float thresh2, // Threshold for stamp finding on readout 2 float stampSpacing, // Spacing between stamps + float sysError, // Relative systematic error in images int size, // Kernel half-size int footprint, // Convolution footprint for stamps @@ -78,5 +79,5 @@ *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, - size, footprint, stampSpacing, mode); + size, footprint, stampSpacing, sysError, mode); if (!*stamps) { psError(psErrorCodeLast(), false, "Unable to find stamps."); @@ -101,5 +102,6 @@ const pmReadout *ro1, const pmReadout *ro2, // Input images int stride, // Size for convolution patches - float sysError, // Relative systematic error + float sysError, // Systematic error in images + float kernelError, // Systematic error in kernel psImageMaskType maskVal, // Value to mask for input psImageMaskType maskBad, // Mask for output bad pixels @@ -151,4 +153,8 @@ PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false); } + if (isfinite(kernelError)) { + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN(kernelError, 1.0, false); + } // Don't care about maskVal // Don't care about maskBad @@ -194,5 +200,5 @@ bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, - psMetadata *analysis, int stride, float sysError, + psMetadata *analysis, int stride, float kernelError, psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac, float badFrac) @@ -264,6 +270,6 @@ } - if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor, - poorFrac, badFrac, mode)) { + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, kernelError, + maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) { psFree(kernels); psFree(regions); @@ -300,5 +306,5 @@ if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, - sysError, region, kernel, true, useFFT)) { + kernelError, region, kernel, true, useFFT)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); psFree(outAnalysis); @@ -330,9 +336,10 @@ int inner, int ringsOrder, int binning, float penalty, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, - int iter, float rej, float sysError, psImageMaskType maskVal, psImageMaskType maskBad, - psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode) -{ - if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor, - poorFrac, badFrac, subMode)) { + int iter, float rej, float sysError, float kernelError, psImageMaskType maskVal, + psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac, + float badFrac, pmSubtractionMode subMode) +{ + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, kernelError, + maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) { return false; } @@ -463,8 +470,8 @@ if (stampsName && strlen(stampsName) > 0) { stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size, - footprint, stampSpacing, subMode); + footprint, stampSpacing, sysError, subMode); } else if (sources) { stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size, - footprint, stampSpacing, subMode); + footprint, stampSpacing, sysError, subMode); } @@ -472,5 +479,5 @@ // doesn't matter. if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2, - stampSpacing, size, footprint, subMode)) { + stampSpacing, sysError, size, footprint, subMode)) { goto MATCH_ERROR; } @@ -538,5 +545,5 @@ if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, - stampThresh1, stampThresh2, stampSpacing, + stampThresh1, stampThresh2, stampSpacing, sysError, size, footprint, subMode)) { goto MATCH_ERROR; @@ -608,5 +615,5 @@ psTrace("psModules.imcombine", 2, "Convolving...\n"); if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, - sysError, region, kernels, true, useFFT)) { + kernelError, region, kernels, true, useFFT)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); goto MATCH_ERROR; Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 26034) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 26035) @@ -39,5 +39,6 @@ int iter, ///< Rejection iterations float rej, ///< Rejection threshold - float sysError, ///< Relative systematic error + float sysError, ///< Relative systematic error in images + float kernelError, ///< Relative systematic error in kernel psImageMaskType maskVal, ///< Value to mask for input psImageMaskType maskBad, ///< Mask for output bad pixels @@ -55,5 +56,5 @@ psMetadata *analysis, ///< Analysis metadata with pre-calculated kernel, region int stride, ///< Size for convolution patches - float sysError, ///< Relative systematic error + float kernelError, ///< Relative systematic error in kernel psImageMaskType maskVal, ///< Value to mask for input psImageMaskType maskBad, ///< Mask for output bad pixels Index: /trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 26034) +++ /trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 26035) @@ -71,5 +71,5 @@ double *sumII, // Sum of I(x)^2/sigma(x)^2 double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 - const pmSubtractionStamp *stamp, // Stamp with variance + const pmSubtractionStamp *stamp, // Stamp const psKernel *target, // Target stamp int kernelIndex, // Index for kernel component @@ -78,13 +78,13 @@ ) { - psKernel *variance = stamp->variance; // Variance, sigma(x)^2 + psKernel *weight = stamp->weight; // Weight image psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { psF32 *in = &target->kernel[y][-footprint]; // Dereference input - psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { - double temp = *in / *wt; // Temporary product + double temp = *in * *wt; // Temporary product *sumI += temp; *sumII += *in * temp; @@ -98,5 +98,5 @@ static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2 double *sumCC, // Sum of conv(x)^2/sigma(x)^2 - const pmSubtractionStamp *stamp, // Stamp with input and variance + const pmSubtractionStamp *stamp, // Stamp with input and weight int kernelIndex, // Index for kernel component int footprint, // Size of region of interest @@ -104,12 +104,12 @@ ) { - psKernel *variance = stamp->variance; // Variance, sigma(x)^2 + psKernel *weight = stamp->weight; // Weight image psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, wt++, conv++) { - double convNoise = *conv / *wt; // Temporary product + double convNoise = *conv * *wt; // Temporary product *sumC += convNoise; *sumCC += *conv * convNoise; @@ -120,5 +120,5 @@ static double accumulateChi2(const psKernel *target, // Target stamp - pmSubtractionStamp *stamp, // Stamp with variance + pmSubtractionStamp *stamp, // Stamp with weight int kernelIndex, // Index for kernel component double coeff, // Coefficient of convolution @@ -129,13 +129,13 @@ { double chi2 = 0.0; - psKernel *variance = stamp->variance; // Variance, sigma(x)^2 + psKernel *weight = stamp->weight; // Weight image psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { psF32 *in = &target->kernel[y][-footprint]; // Dereference input - psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { - chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt; + chi2 += PS_SQR(*in - bg - coeff * *conv) * *wt; } } @@ -146,10 +146,10 @@ // Return the initial value of chi^2 static double initialChi2(const psKernel *target, // Target stamp - const pmSubtractionStamp *stamp, // Stamp with variance + const pmSubtractionStamp *stamp, // Stamp int footprint, // Size of convolution pmSubtractionMode mode // Mode of subtraction ) { - psKernel *variance = stamp->variance; // Variance map + psKernel *weight = stamp->weight; // Weight image psKernel *source; // Source stamp switch (mode) { @@ -167,9 +167,9 @@ for (int y = -footprint; y <= footprint; y++) { psF32 *in = &target->kernel[y][-footprint]; // Dereference input - psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *ref = &source->kernel[y][-footprint]; // Derference reference for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { float diff = *in - *ref; // Temporary value - chi2 += PS_SQR(diff) / *wt; + chi2 += PS_SQR(diff) * *wt; } } @@ -180,5 +180,5 @@ // Subtract a convolution from the input static void subtractConvolution(psKernel *target, // Target stamp - const pmSubtractionStamp *stamp, // Stamp with variance + const pmSubtractionStamp *stamp, // Stamp int kernelIndex, // Index for kernel component float coeff, // Coefficient of subtraction @@ -288,9 +288,9 @@ // This sum is invariant to the kernel - psKernel *variance = stamp->variance; // Variance map for stamp + psKernel *weight = stamp->weight; // Weight image for (int v = -footprint; v <= footprint; v++) { - psF32 *wt = &variance->kernel[v][-footprint]; // Dereference variance map + psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight for (int u = -footprint; u <= footprint; u++, wt++) { - sum1 += 1.0 / *wt; + sum1 += 1.0 * *wt; } } Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 26034) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 26035) @@ -54,5 +54,5 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->variance); + psFree(stamp->weight); psFree(stamp->convolutions1); psFree(stamp->convolutions2); @@ -180,5 +180,5 @@ pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, - int footprint, float spacing) + int footprint, float spacing, float sysErr) { pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return @@ -221,4 +221,5 @@ list->flux = NULL; list->footprint = footprint; + list->sysErr = sysErr; return list; @@ -256,5 +257,5 @@ outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL; outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL; - outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL; + outStamp->weight = inStamp->weight ? psKernelCopy(inStamp->weight) : NULL; if (inStamp->convolutions1) { @@ -305,5 +306,5 @@ stamp->image1 = NULL; stamp->image2 = NULL; - stamp->variance = NULL; + stamp->weight = NULL; stamp->convolutions1 = NULL; stamp->convolutions2 = NULL; @@ -322,5 +323,5 @@ const psImage *image2, const psImage *subMask, const psRegion *region, float thresh1, float thresh2, - int size, int footprint, float spacing, + int size, int footprint, float spacing, float sysErr, pmSubtractionMode mode) { @@ -379,5 +380,5 @@ if (!stamps) { - stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing); + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr); } @@ -452,8 +453,8 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->variance); + psFree(stamp->weight); psFree(stamp->convolutions1); psFree(stamp->convolutions2); - stamp->image1 = stamp->image2 = stamp->variance = NULL; + stamp->image1 = stamp->image2 = stamp->weight = NULL; stamp->convolutions1 = stamp->convolutions2 = NULL; @@ -487,5 +488,5 @@ const psImage *image, const psImage *subMask, const psRegion *region, int size, int footprint, - float spacing, pmSubtractionMode mode) + float spacing, float sysErr, pmSubtractionMode mode) { @@ -507,5 +508,6 @@ int numStars = x->n; // Number of stars pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, - region, footprint, spacing); // Stamp list + region, footprint, sysErr, + spacing); // Stamp list int numStamps = stamps->num; // Number of stamps @@ -616,8 +618,10 @@ PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); } - PS_ASSERT_IMAGE_NON_NULL(variance, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false); - PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false); - PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); + if (variance) { + PS_ASSERT_IMAGE_NON_NULL(variance, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false); + PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false); + PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); + } int numCols = image1->numCols, numRows = image1->numRows; // Size of images @@ -646,5 +650,5 @@ assert(stamp->image1 == NULL); assert(stamp->image2 == NULL); - assert(stamp->variance == NULL); + assert(stamp->weight == NULL); psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest @@ -660,7 +664,29 @@ } - psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp - stamp->variance = psKernelAllocFromImage(wtSub, size, size); - psFree(wtSub); // Drop reference + psKernel *weight = stamp->weight = psKernelAlloc(-size, size, -size, size); // Weight = 1/variance + if (variance) { + psImage *varSub = psImageSubset(variance, region); // Subimage with stamp + psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp + if (isfinite(stamps->sysErr) && stamps->sysErr > 0) { + float sysErr = 0.5 * stamps->sysErr; // Systematic error + psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + weight->kernel[y][x] = 1.0 / (var->kernel[y][x] + + sysErr * (image1->kernel[y][x] + image2->kernel[y][x])); + } + } + } else { + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + weight->kernel[y][x] = 1.0 / var->kernel[y][x]; + } + } + } + psFree(var); + psFree(varSub); + } else { + psImageInit(weight->image, 1.0); + } stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; @@ -672,5 +698,5 @@ pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, + int size, int footprint, float spacing, float sysErr, pmSubtractionMode mode) { @@ -702,5 +728,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, - footprint, spacing, mode); // Stamps to return + footprint, spacing, sysErr, mode); // Stamps psFree(x); psFree(y); @@ -716,5 +742,5 @@ pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, + int size, int footprint, float spacing, float sysErr, pmSubtractionMode mode) { @@ -734,5 +760,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint, - spacing, mode); + spacing, sysErr, mode); psFree(data); Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 26034) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 26035) @@ -24,4 +24,5 @@ psArray *flux; ///< Fluxes for possible stamps (or NULL) int footprint; ///< Half-size of stamps + float sysErr; ///< Systematic error } pmSubtractionStampList; @@ -31,5 +32,6 @@ const psRegion *region, // Region for stamps, or NULL int footprint, // Half-size of stamps - float spacing // Rough average spacing between stamps + float spacing, // Rough average spacing between stamps + float sysErr // Relative systematic error or NAN ); @@ -68,5 +70,5 @@ psKernel *image1; ///< Reference image postage stamp psKernel *image2; ///< Input image postage stamp - psKernel *variance; ///< Variance image postage stamp, or NULL + psKernel *weight; ///< Weight image (1/variance) postage stamp, or NULL psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL @@ -91,4 +93,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float sysErr, ///< Relative systematic error in images pmSubtractionMode mode ///< Mode for subtraction ); @@ -103,4 +106,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float sysErr, ///< Systematic error in images pmSubtractionMode mode ///< Mode for subtraction ); @@ -115,4 +119,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float sysErr, ///< Systematic error in images pmSubtractionMode mode ///< Mode for subtraction ); @@ -127,4 +132,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float sysErr, ///< Systematic error in images pmSubtractionMode mode ///< Mode for subtraction );