Changeset 21363 for trunk/psModules/src/imcombine

Timestamp:

Feb 5, 2009, 4:31:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

Location:

trunk/psModules/src/imcombine

Files:

• pmPSFEnvelope.c (modified) (5 diffs)
• pmReadoutCombine.c (modified) (14 diffs)
• pmReadoutCombine.h (modified) (4 diffs)
• pmStack.c (modified) (8 diffs)
• pmSubtraction.c (modified) (21 diffs)
• pmSubtraction.h (modified) (2 diffs)
• pmSubtractionEquation.c (modified) (25 diffs)
• pmSubtractionMatch.c (modified) (9 diffs)
• pmSubtractionParams.c (modified) (10 diffs)
• pmSubtractionStamps.c (modified) (7 diffs)
• pmSubtractionStamps.h (modified) (2 diffs)
• pmSubtractionThreads.c (modified) (2 diffs)

trunk/psModules/src/imcombine/pmPSFEnvelope.c

@@ -36 +36 @@
 #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 +149 @@
             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 +216 @@
     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 +243 @@
 
         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 +281 @@
     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);

trunk/psModules/src/imcombine/pmReadoutCombine.c

@@ -37 +37 @@
     params->iter = 1;
     params->rej = INFINITY;
-    params->weights = false;
+    params->variances = false;
 
     return params;
@@ -75 +75 @@
     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 +82 @@
     // 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 +121 @@
 
     // 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 +177 @@
             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 +190 @@
     }
 
-    // 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 +252 @@
     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 +281 @@
     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 +325 @@
                 }
 
-                if (params->weights) {
-                    weightsData[r] = readout->weight->data.F32[yIn][xIn];
+                if (params->variances) {
+                    variancesData[r] = readout->variance->data.F32[yIn][xIn];
                 }
 
@@ -332 +334 @@
                 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 +378 @@
 
             // 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 +388 @@
                 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 +395 @@
                 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 +414 @@
     psFree(pixels);
     psFree(mask);
-    psFree(weights);
+    psFree(variances);
     psFree(errors);
     psFree(stats);

trunk/psModules/src/imcombine/pmReadoutCombine.h

@@ -5 +5 @@
  * @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 +22 @@
 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 +29 @@
     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 +41 @@
 /// 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

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @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 +30 @@
 #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 +264 @@
 
         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 +443 @@
     *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 +474 @@
         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 +607 @@
         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 +647 @@
         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 +700 @@
         }
 
-        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;
         }
 

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -60 +60 @@
 
     // 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 +287 @@
 
 // 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 +302 @@
     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 +319 @@
         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 +326 @@
         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 +361 @@
 // 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 +381 @@
 {
     *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
-    if (weight || subMask) {
-        *kernelWeight = varianceKernel(*kernelWeight, *kernelImage);
+    if (variance || subMask) {
+        *kernelVariance = varianceKernel(*kernelVariance, *kernelImage);
     }
 
@@ -398 +398 @@
     }
 
-    // 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 +885 @@
                 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 +915 @@
 }
 
-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 +922 @@
     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 +996 @@
 
 
-// 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 +1017 @@
 
     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 +1027 @@
 
     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 +1155 @@
             // 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 +1172 @@
             // 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 +1239 @@
     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 +1306 @@
                 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 +1341 @@
         // 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 +1366 @@
         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 +1374 @@
         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;

trunk/psModules/src/imcombine/pmSubtraction.h

@@ -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 +70 @@
                               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?
     );
 

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -24 +24 @@
 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 +31 @@
     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 +42 @@
                                            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 +50 @@
     )
 {
-    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 +79 @@
                                            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 +87 @@
     )
 {
-    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 +120 @@
                                   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 +135 @@
             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 +151 @@
                             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 +171 @@
 
     // 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 +185 @@
 
             // 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 +195 @@
             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 +218 @@
         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 +253 @@
                             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 +277 @@
         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 +297 @@
         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 +329 @@
                                  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 +348 @@
     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 +356 @@
         // 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 +366 @@
         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 +559 @@
       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 +581 @@
 #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 +1033 @@
 
         // 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 +1098 @@
         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 +1141 @@
             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);
         }

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -60 +60 @@
                       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 +80 @@
 
     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 +111 @@
             conv1->mask = NULL;
         }
-        if (conv1->weight) {
-            psFree(conv1->weight);
-            conv1->weight = NULL;
+        if (conv1->variance) {
+            psFree(conv1->variance);
+            conv1->variance = NULL;
         }
     }
@@ -126 +126 @@
             conv2->mask = NULL;
         }
-        if (conv2->weight) {
-            psFree(conv2->weight);
-            conv2->weight = NULL;
+        if (conv2->variance) {
+            psFree(conv2->variance);
+            conv2->variance = NULL;
         }
     }
@@ -191 +191 @@
     }
 
-    // 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 +274 @@
             // 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 +336 @@
                 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 +439 @@
     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 +482 @@
     psFree(kernels);
     psFree(stamps);
-    psFree(weight);
+    psFree(variance);
     psFree(rng);
     return false;

trunk/psModules/src/imcombine/pmSubtractionParams.c

@@ -71 +71 @@
                             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 +78 @@
     )
 {
-    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 +98 @@
 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 +104 @@
     )
 {
-    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 +120 @@
 
 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 +129 @@
 {
     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 +146 @@
 // 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 +167 @@
     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 +180 @@
 // 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 +288 @@
 
         // 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;

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -54 +54 @@
     psFree(stamp->image1);
     psFree(stamp->image2);
-    psFree(stamp->weight);
+    psFree(stamp->variance);
     psFree(stamp->convolutions1);
     psFree(stamp->convolutions2);
@@ -211 +211 @@
     stamp->image1 = NULL;
     stamp->image2 = NULL;
-    stamp->weight = NULL;
+    stamp->variance = NULL;
     stamp->convolutions1 = NULL;
     stamp->convolutions2 = NULL;
@@ -333 +333 @@
                 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 +480 @@
 
 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 +490 @@
         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 +520 @@
         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 +534 @@
         }
 
-        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
 

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -60 +60 @@
     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 +120 @@
                                 psImage *image1, ///< Reference image
                                 psImage *image2, ///< Input image (or NULL)
-                                psImage *weight, ///< Weight (variance) map
+                                psImage *variance, ///< Variance map
                                 int kernelSize ///< Kernel half-size
     );

trunk/psModules/src/imcombine/pmSubtractionThreads.c

@@ -27 +27 @@
     // XXX     psMutexInit(ro->image);
     // XXX }
-    // XXX if (ro->weight) {
-    // XXX     psMutexInit(ro->weight);
+    // XXX if (ro->variance) {
+    // XXX     psMutexInit(ro->variance);
     // XXX }
 
@@ -43 +43 @@
     // XXX     psMutexDestroy(ro->image);
     // XXX }
-    // XXX if (ro->weight) {
-    // XXX     psMutexDestroy(ro->weight);
+    // XXX if (ro->variance) {
+    // XXX     psMutexDestroy(ro->variance);
     // XXX }