Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 21351) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 21363) @@ -36,6 +36,6 @@ #define PEAK_FLUX 1.0e4 // Peak flux for each source #define SKY_VALUE 0.0e0 // Sky value for fake image -#define WEIGHT_VAL 3.0 // Weighting for image -#define WEIGHT_FACTOR 10.0 // Factor to multiply image by to get weighting +#define VARIANCE_VAL 3.0 // Variance for image +#define VARIANCE_FACTOR 10.0 // Factor to multiply image by to get variance #define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF #define SOURCE_FIT_ITERATIONS 100 // Number of iterations for source fitting @@ -149,5 +149,5 @@ pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source psAssert (model, "failed to generate model: should this be an error or not?"); - float srcRadius = model->modelRadius(model->params, PS_SQR(WEIGHT_VAL)); // Radius for source + float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source if (srcRadius > maxRadius) { maxRadius = srcRadius; @@ -216,11 +216,12 @@ psFree(envelope); - // XXX Setting the weight seems to be an art + // XXX Setting the variance seems to be an art // Can't set it too high so that pixels are rejected as insignificant // Can't set it too low so that it's hard to get to the minimum // Have also tried: - // *** readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "*", readout->image); - // *** readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "*", psScalarAlloc(WEIGHT_FACTOR, PS_TYPE_F32)); - readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "+", psScalarAlloc(WEIGHT_VAL, PS_TYPE_F32)); + // *** readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "*", readout->image); + // *** readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "*", psScalarAlloc(VARIANCE_FACTOR, PS_TYPE_F32)); + readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "+", + psScalarAlloc(VARIANCE_VAL, PS_TYPE_F32)); readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); psImageInit(readout->mask, 0); @@ -242,9 +243,9 @@ psFree(source->pixels); - psFree(source->weight); + psFree(source->variance); psFree(source->maskView); psFree(source->maskObj); source->pixels = NULL; - source->weight = NULL; + source->variance = NULL; source->maskView = NULL; source->maskObj = NULL; @@ -280,5 +281,5 @@ options->psfFieldYo = 0; - pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, WEIGHT_VAL, true); + pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true); pmPSFtry *try = pmPSFtryModel(fakes, modelName, options, 0, 0xff); Index: trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 21351) +++ trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 21363) @@ -37,5 +37,5 @@ params->iter = 1; params->rej = INFINITY; - params->weights = false; + params->variances = false; return params; @@ -75,5 +75,5 @@ psStringAppend(&comment, "Combining using statistic: %x", params->combine); if (!hdu->header) { - hdu->header = psMetadataAlloc(); + hdu->header = psMetadataAlloc(); } psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); @@ -82,14 +82,14 @@ // note the clipping parameters, if used if (params->combine == PS_STAT_CLIPPED_MEAN) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - - // note the use of weights - if (params->weights) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + // note the use of variances + if (params->variances) { psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - "Using input weights to combine images", ""); + "Using input variances to combine images", ""); } @@ -121,5 +121,5 @@ // generate the required output images based on the specified sizes - pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->weights, params->blank); + pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->variances, params->blank); psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); @@ -177,7 +177,7 @@ return false; } - if (params->weights && !readout->weight) { + if (params->variances && !readout->variance) { psError(PS_ERR_UNEXPECTED_NULL, true, - "Rejection based on weights requested, but no weights supplied for image %d.\n", i); + "Rejection based on variances requested, but no variances supplied for image %d.\n", i); return false; } @@ -190,10 +190,12 @@ } - // pthread_t id = pthread_self(); - // char name[64]; - // sprintf (name, "%x", (unsigned int) id); - // psTimerStart (name); - - psStatsOptions combineStdev = 0; // Statistics option for weights +#if 0 + pthread_t id = pthread_self(); + char name[64]; + sprintf(name, "%x", (unsigned int)id); + psTimerStart(name); +#endif + + psStatsOptions combineStdev = 0; // Statistics option for variances switch (params->combine) { case PS_STAT_SAMPLE_MEAN: @@ -250,10 +252,10 @@ psVectorMaskType *maskData = mask->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference mask - psVector *weights = NULL; // Stack of weights - psVector *errors = NULL; // Stack of errors (sqrt of variance/weights), for psVectorStats - psF32 *weightsData = NULL; // Dereference weights - if (params->weights) { - weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights - weightsData = weights->data.F32; + psVector *variances = NULL; // Stack of variances + psVector *errors = NULL; // Stack of errors (sqrt of variance), for psVectorStats + psF32 *variancesData = NULL; // Dereference variances + if (params->variances) { + variances = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of variances + variancesData = variances->data.F32; } psVector *index = NULL; // The indices to sort the pixels @@ -279,7 +281,7 @@ psF32 **outputImage = output->image->data.F32; // Output image psImageMaskType **outputMask = output->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Output mask - psF32 **outputWeight = NULL; // Output weight map - if (output->weight) { - outputWeight = output->weight->data.F32; + psF32 **outputVariance = NULL; // Output variance map + if (output->variance) { + outputVariance = output->variance->data.F32; } @@ -323,6 +325,6 @@ } - if (params->weights) { - weightsData[r] = readout->weight->data.F32[yIn][xIn]; + if (params->variances) { + variancesData[r] = readout->variance->data.F32[yIn][xIn]; } @@ -332,6 +334,6 @@ if (scale) { pixelsData[r] *= invScale->data.F32[r]; - if (params->weights) { - weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; + if (params->variances) { + variancesData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; } } @@ -376,6 +378,6 @@ // XXXXX this step probably is very expensive : convert errors to variance everywhere? - if (params->weights) { - errors = (psVector*)psUnaryOp(errors, weights, "sqrt"); + if (params->variances) { + errors = (psVector*)psUnaryOp(errors, variances, "sqrt"); } @@ -386,6 +388,6 @@ outputImage[yOut][xOut] = NAN; outputMask[yOut][xOut] = params->blank; - if (params->weights) { - outputWeight[yOut][xOut] = NAN; + if (params->variances) { + outputVariance[yOut][xOut] = NAN; } sigma->data.F32[yOut][xOut] = NAN; @@ -393,7 +395,7 @@ outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); outputMask[yOut][xOut] = isfinite(outputImage[yOut][xOut]) ? 0 : params->blank; - if (params->weights) { + if (params->variances) { float stdev = psStatsGetValue(stats, combineStdev); - outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance + outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance // XXXX this is not the correct formal error. // also, the weighted mean is not obviously the correct thing here @@ -412,5 +414,5 @@ psFree(pixels); psFree(mask); - psFree(weights); + psFree(variances); psFree(errors); psFree(stats); Index: trunk/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 21351) +++ trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 21363) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004-2006 Institute for Astronomy, University of Hawaii */ @@ -22,6 +22,6 @@ typedef struct { psStatsOptions combine; ///< Statistic to use when performing the combination - psImageMaskType maskVal; ///< Mask value - psImageMaskType blank; ///< Mask value to give blank (i.e., no data) pixels + psImageMaskType maskVal; ///< Mask value + psImageMaskType blank; ///< Mask value to give blank (i.e., no data) pixels int nKeep; ///< Mimimum number of pixels to keep float fracHigh; ///< Fraction of high pixels to immediately throw @@ -29,5 +29,5 @@ int iter; ///< Number of iterations for clipping (for CLIPPED_MEAN only) float rej; ///< Rejection threshould for clipping (for CLIPPED_MEAN only) - bool weights; ///< Use the supplied weights (instead of calculated stdev)? + bool variances; ///< Use the supplied variances (instead of calculated stdev)? } pmCombineParams; @@ -41,5 +41,5 @@ /// Combine multiple readouts, applying zero and scale, with optional minmax clipping bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned - const psArray *inputs, ///< Array of input readouts (F32 image and weight, U8 mask) + const psArray *inputs, ///< Array of input readouts const psVector *zero, ///< Zero corrections to subtract from input, or NULL const psVector *scale, ///< Scale corrections to divide into input, or NULL Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 21351) +++ trunk/psModules/src/imcombine/pmStack.c (revision 21363) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -30,5 +30,6 @@ #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time #define PIXEL_MAP_BUFFER 2 // Number of entries to add to pixel map at a time -//#define VARIANCE_FACTORS // Use variance factors when calculating the variances? +#define VARIANCE_FACTORS // Use variance factors when calculating the variances? +//#define ADD_VARIANCE // Allow additional variance (besides variance factor)? #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly @@ -263,5 +264,5 @@ psImage *image = data->readout->image; // Image of interest - psImage *variance = data->readout->weight; // Variance ("weight") map of interest + psImage *variance = data->readout->variance; // Variance map of interest pixelData->data.F32[num] = image->data.F32[yIn][xIn]; if (variance) { @@ -442,9 +443,9 @@ *numCols = data->readout->image->numCols; *numRows = data->readout->image->numRows; - if (data->readout->weight) { + if (data->readout->variance) { *haveVariances = true; - PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); - PS_ASSERT_IMAGE_TYPE(data->readout->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false); + PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } *haveRejects = (data->reject != NULL); @@ -473,7 +474,7 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false); if (*haveVariances) { - PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); - PS_ASSERT_IMAGE_TYPE(data->readout->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false); + PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } } @@ -606,14 +607,11 @@ weights->data.F32[i] = data->weight; stack->data[i] = psMemIncrRefCounter(data->readout); - // Variance factor - float vf = psMetadataLookupF32(NULL, data->readout->parent->concepts, "CELL.VARFACTOR"); // Var factor - if (!isfinite(vf)) { - psWarning("Non-finite CELL.VARFACTOR for image %d --- setting to unity.", i); - vf = 1.0; - } - varFactors->data.F32[i] = vf; + varFactors->data.F32[i] = data->readout->covariance ? + psImageCovarianceFactor(data->readout->covariance) : 1.0; +#if 0 if (isfinite(data->addVariance)) { varFactors->data.F32[i] *= data->addVariance; } +#endif if (!haveRejects && !data->inspect) { data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); @@ -649,5 +647,5 @@ psImage *combinedImage = combined->image; // Combined image psImage *combinedMask = combined->mask; // Combined mask - psImage *combinedVariance = combined->weight; // Combined variance map + psImage *combinedVariance = combined->variance; // Combined variance map psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, @@ -702,8 +700,8 @@ } - psImage *combinedVariance = combined->weight; // Combined variance map + psImage *combinedVariance = combined->variance; // Combined variance map if (haveVariances && !combinedVariance) { - combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedVariance = combined->weight; + combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedVariance = combined->variance; } 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; Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 21363) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -70,4 +70,10 @@ float sigmaRej, ///< Number of RMS deviations above zero at which to reject int footprint ///< Half-size of stamp + ); + +/// Generate the convolution kernel +psKernel *pmSubtractionKernel(const pmSubtractionKernels *kernels, ///< Kernel parameters + float x, float y, ///< Normalised position [-1,1] for which to generate image + bool wantDual ///< Calculate for the dual kernel? ); Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 21363) @@ -24,5 +24,5 @@ static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication const psKernel *image2, // Second image in multiplication - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image int footprint // (Half-)Size of stamp ) @@ -31,5 +31,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -42,5 +42,5 @@ const psKernel *image1, // First image in multiplication const psKernel *image2, // Second image in multiplication - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image const psImage *polyValues, // Spatial polynomial values int numKernels, // Number of kernel basis functions @@ -50,5 +50,5 @@ ) { - double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products if (!isfinite(sum)) { return false; @@ -79,5 +79,5 @@ const psKernel *image1, // First image in multiplication const psKernel *image2, // Second image in multiplication - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image const psImage *polyValues, // Spatial polynomial values int numKernels, // Number of kernel basis functions @@ -87,5 +87,5 @@ ) { - double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products if (!isfinite(sum)) { return false; @@ -120,5 +120,5 @@ const psArray *convolutions1, // Convolutions for element 1 const psArray *convolutions2, // Convolutions for element 2 - const psKernel *weight, // Weight image + const psKernel *variance, // Variance image const psImage *polyValues, // Polynomial values int numKernels, // Number of kernel basis functions @@ -135,5 +135,5 @@ psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element - if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, weight, polyValues, numKernels, + if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, variance, polyValues, numKernels, footprint, spatialOrder, symmetric)) { psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j); @@ -151,5 +151,5 @@ const psArray *convolutions, // Convolutions of source with kernels const psKernel *input, // Input stamp, or NULL - const psKernel *weight, // Weight stamp + const psKernel *variance, // Variance stamp const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp @@ -171,5 +171,5 @@ // Square part of the matrix (convolution-convolution products) - if (!calculateMatrixSquare(matrix, convolutions, convolutions, weight, polyValues, numKernels, + if (!calculateMatrixSquare(matrix, convolutions, convolutions, variance, polyValues, numKernels, spatialOrder, footprint)) { return false; @@ -185,5 +185,5 @@ // Normalisation-convolution terms - if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels, + if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, variance, polyValues, numKernels, footprint, spatialOrder, true)) { psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); @@ -195,5 +195,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -218,5 +218,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x]; + double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x]; double value = input->kernel[y][x] * invNoise2; sumI += value; @@ -253,5 +253,5 @@ const psKernel *input, // Input stamp, or NULL if !normAndBG const psKernel *target, // Target stamp - const psKernel *weight, // Weight stamp + const psKernel *variance, // Variance stamp const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp @@ -277,5 +277,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -297,5 +297,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x]; + float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x]; sumIT += value * input->kernel[y][x]; sumT += value; @@ -329,5 +329,5 @@ const psArray *convolutions2, // Convolutions of image 2 const psKernel *image1, // Image 1 stamp - const psKernel *weight, // Weight stamp + const psKernel *variance, // Variance stamp const psImage *polyValues, // Spatial polynomial values int footprint // (Half-)Size of stamp @@ -348,5 +348,5 @@ PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); - if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, weight, polyValues, numKernels, + if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, variance, polyValues, numKernels, spatialOrder, footprint)) { return false; @@ -356,5 +356,5 @@ // Normalisation psKernel *conv = convolutions2->data[i]; // Convolution - if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, weight, polyValues, numKernels, + if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, variance, polyValues, numKernels, footprint, spatialOrder, false)) { psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); @@ -366,5 +366,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; + sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; } } @@ -559,13 +559,13 @@ case PM_SUBTRACTION_MODE_1: status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, - stamp->weight, polyValues, footprint, true); + stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint, true); + stamp->image2, stamp->variance, polyValues, footprint, true); break; case PM_SUBTRACTION_MODE_2: status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, - stamp->weight, polyValues, footprint, true); + stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, - stamp->image1, stamp->weight, polyValues, footprint, true); + stamp->image1, stamp->variance, polyValues, footprint, true); break; case PM_SUBTRACTION_MODE_DUAL: @@ -581,14 +581,14 @@ #endif status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, - stamp->weight, polyValues, footprint, true); + stamp->variance, polyValues, footprint, true); status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, - stamp->weight, polyValues, footprint, false); + stamp->variance, polyValues, footprint, false); status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, - stamp->convolutions2, stamp->image1, stamp->weight, polyValues, + stamp->convolutions2, stamp->image1, stamp->variance, polyValues, footprint); status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, - stamp->image2, stamp->weight, polyValues, footprint, true); + stamp->image2, stamp->variance, polyValues, footprint, true); status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, - stamp->image2, stamp->weight, polyValues, footprint, false); + stamp->image2, stamp->variance, polyValues, footprint, false); break; default: @@ -1033,5 +1033,5 @@ // Calculate residuals - psKernel *weight = stamp->weight; // Weight postage stamp + psKernel *variance = stamp->variance; // Variance postage stamp psImageInit(residual->image, 0.0); if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { @@ -1098,5 +1098,5 @@ 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]; + double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x]; deviation += dev; #ifdef TESTING @@ -1141,10 +1141,10 @@ psFitsClose(fits); } - if (stamp->weight) { + if (stamp->variance) { psString filename = NULL; - psStringAppend(&filename, "stamp_weight_%03d.fits", i); + psStringAppend(&filename, "stamp_variance_%03d.fits", i); psFits *fits = psFitsOpen(filename, "w"); psFree(filename); - psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); + psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL); psFitsClose(fits); } Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 21363) @@ -60,5 +60,5 @@ const pmReadout *ro2, // Readout 2 const psImage *subMask, // Mask for subtraction, or NULL - psImage *weight, // Weight map + psImage *variance, // Variance map const psRegion *region, // Region of interest, or NULL float threshold, // Threshold for stamp finding @@ -80,5 +80,5 @@ psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, size)) { + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -111,7 +111,7 @@ conv1->mask = NULL; } - if (conv1->weight) { - psFree(conv1->weight); - conv1->weight = NULL; + if (conv1->variance) { + psFree(conv1->variance); + conv1->variance = NULL; } } @@ -126,7 +126,7 @@ conv2->mask = NULL; } - if (conv2->weight) { - psFree(conv2->weight); - conv2->weight = NULL; + if (conv2->variance) { + psFree(conv2->variance); + conv2->variance = NULL; } } @@ -191,15 +191,15 @@ } - // Where does our weight map come from? - // Getting the weight exactly right is not necessary --- it's just used for weighting. - psImage *weight = NULL; // Weight image to use - if (ro1->weight && ro2->weight) { - weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight); - } else if (ro1->weight) { - weight = psMemIncrRefCounter(ro1->weight); - } else if (ro2->weight) { - weight = psMemIncrRefCounter(ro2->weight); + // Where does our variance map come from? + // Getting the variance exactly right is not necessary --- it's just used for weighting. + psImage *variance = NULL; // Variance image to use + if (ro1->variance && ro2->variance) { + variance = (psImage*)psBinaryOp(NULL, ro1->variance, "+", ro2->variance); + } else if (ro1->variance) { + variance = psMemIncrRefCounter(ro1->variance); + } else if (ro2->variance) { + variance = psMemIncrRefCounter(ro2->variance); } else { - weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); + variance = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); } @@ -274,5 +274,5 @@ // We get the stamps here; we will also attempt to get stamps at the first iteration, but it // doesn't matter. - if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing, + if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, threshold, stampSpacing, size, footprint, subMode)) { goto MATCH_ERROR; @@ -336,5 +336,5 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing, + if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, threshold, stampSpacing, size, footprint, subMode)) { goto MATCH_ERROR; @@ -439,8 +439,8 @@ psFree(subMask); subMask = NULL; - psFree(weight); - weight = NULL; - - if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBad)) { + psFree(variance); + variance = NULL; + + if (!pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) { psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); goto MATCH_ERROR; @@ -482,5 +482,5 @@ psFree(kernels); psFree(stamps); - psFree(weight); + psFree(variance); psFree(rng); return false; Index: trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 21363) @@ -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 weight + const pmSubtractionStamp *stamp, // Stamp with variance const psKernel *target, // Target stamp int kernelIndex, // Index for kernel component @@ -78,10 +78,10 @@ ) { - psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *variance = stamp->variance; // Variance, sigma(x)^2 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 = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { @@ -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 weight + const pmSubtractionStamp *stamp, // Stamp with input and variance int kernelIndex, // Index for kernel component int footprint, // Size of region of interest @@ -104,9 +104,9 @@ ) { - psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *variance = stamp->variance; // Variance, sigma(x)^2 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, wt++, conv++) { @@ -120,5 +120,5 @@ static double accumulateChi2(const psKernel *target, // Target stamp - pmSubtractionStamp *stamp, // Stamp with weight + pmSubtractionStamp *stamp, // Stamp with variance int kernelIndex, // Index for kernel component double coeff, // Coefficient of convolution @@ -129,10 +129,10 @@ { double chi2 = 0.0; - psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *variance = stamp->variance; // Variance, sigma(x)^2 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 = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { @@ -146,10 +146,10 @@ // Return the initial value of chi^2 static double initialChi2(const psKernel *target, // Target stamp - const pmSubtractionStamp *stamp, // Stamp with weight + const pmSubtractionStamp *stamp, // Stamp with variance int footprint, // Size of convolution pmSubtractionMode mode // Mode of subtraction ) { - psKernel *weight = stamp->weight; // Weight map + psKernel *variance = stamp->variance; // Variance map psKernel *source; // Source stamp switch (mode) { @@ -167,5 +167,5 @@ for (int y = -footprint; y <= footprint; y++) { psF32 *in = &target->kernel[y][-footprint]; // Dereference input - psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance psF32 *ref = &source->kernel[y][-footprint]; // Derference reference for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { @@ -180,5 +180,5 @@ // Subtract a convolution from the input static void subtractConvolution(psKernel *target, // Target stamp - const pmSubtractionStamp *stamp, // Stamp with weight + const pmSubtractionStamp *stamp, // Stamp with variance int kernelIndex, // Index for kernel component float coeff, // Coefficient of subtraction @@ -288,7 +288,7 @@ // This sum is invariant to the kernel - psKernel *weight = stamp->weight; // Weight map for stamp + psKernel *variance = stamp->variance; // Variance map for stamp for (int v = -footprint; v <= footprint; v++) { - psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight map + psF32 *wt = &variance->kernel[v][-footprint]; // Dereference variance map for (int u = -footprint; u <= footprint; u++, wt++) { sum1 += 1.0 / *wt; Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 21363) @@ -54,5 +54,5 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->weight); + psFree(stamp->variance); psFree(stamp->convolutions1); psFree(stamp->convolutions2); @@ -211,5 +211,5 @@ stamp->image1 = NULL; stamp->image2 = NULL; - stamp->weight = NULL; + stamp->variance = NULL; stamp->convolutions1 = NULL; stamp->convolutions2 = NULL; @@ -333,8 +333,8 @@ psFree(stamp->image1); psFree(stamp->image2); - psFree(stamp->weight); + psFree(stamp->variance); psFree(stamp->convolutions1); psFree(stamp->convolutions2); - stamp->image1 = stamp->image2 = stamp->weight = NULL; + stamp->image1 = stamp->image2 = stamp->variance = NULL; stamp->convolutions1 = stamp->convolutions2 = NULL; @@ -480,5 +480,5 @@ bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, - psImage *weight, int kernelSize) + psImage *variance, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -490,7 +490,7 @@ PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); } - PS_ASSERT_IMAGE_NON_NULL(weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false); - PS_ASSERT_IMAGE_TYPE(weight, 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); @@ -520,5 +520,5 @@ assert(stamp->image1 == NULL); assert(stamp->image2 == NULL); - assert(stamp->weight == NULL); + assert(stamp->variance == NULL); psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest @@ -534,6 +534,6 @@ } - psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp - stamp->weight = psKernelAllocFromImage(wtSub, size, size); + psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp + stamp->variance = psKernelAllocFromImage(wtSub, size, size); psFree(wtSub); // Drop reference Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 21363) @@ -60,5 +60,5 @@ psKernel *image1; ///< Reference image postage stamp psKernel *image2; ///< Input image postage stamp - psKernel *weight; ///< Weight image postage stamp, or NULL + psKernel *variance; ///< Variance image 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 @@ -120,5 +120,5 @@ psImage *image1, ///< Reference image psImage *image2, ///< Input image (or NULL) - psImage *weight, ///< Weight (variance) map + psImage *variance, ///< Variance map int kernelSize ///< Kernel half-size ); Index: trunk/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 21351) +++ trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 21363) @@ -27,6 +27,6 @@ // XXX psMutexInit(ro->image); // XXX } - // XXX if (ro->weight) { - // XXX psMutexInit(ro->weight); + // XXX if (ro->variance) { + // XXX psMutexInit(ro->variance); // XXX } @@ -43,6 +43,6 @@ // XXX psMutexDestroy(ro->image); // XXX } - // XXX if (ro->weight) { - // XXX psMutexDestroy(ro->weight); + // XXX if (ro->variance) { + // XXX psMutexDestroy(ro->variance); // XXX }