Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 25120)
+++ trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 26035)
@@ -54,5 +54,5 @@
     psFree(stamp->image1);
     psFree(stamp->image2);
-    psFree(stamp->variance);
+    psFree(stamp->weight);
     psFree(stamp->convolutions1);
     psFree(stamp->convolutions2);
@@ -180,5 +180,5 @@
 
 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
-                                                    int footprint, float spacing)
+                                                    int footprint, float spacing, float sysErr)
 {
     pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
@@ -221,4 +221,5 @@
     list->flux = NULL;
     list->footprint = footprint;
+    list->sysErr = sysErr;
 
     return list;
@@ -256,5 +257,5 @@
         outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL;
         outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL;
-        outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL;
+        outStamp->weight = inStamp->weight ? psKernelCopy(inStamp->weight) : NULL;
 
         if (inStamp->convolutions1) {
@@ -305,5 +306,5 @@
     stamp->image1 = NULL;
     stamp->image2 = NULL;
-    stamp->variance = NULL;
+    stamp->weight = NULL;
     stamp->convolutions1 = NULL;
     stamp->convolutions2 = NULL;
@@ -322,5 +323,5 @@
                                                 const psImage *image2, const psImage *subMask,
                                                 const psRegion *region, float thresh1, float thresh2,
-                                                int size, int footprint, float spacing,
+                                                int size, int footprint, float spacing, float sysErr,
                                                 pmSubtractionMode mode)
 {
@@ -379,5 +380,5 @@
 
     if (!stamps) {
-        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing);
+        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr);
     }
 
@@ -452,8 +453,8 @@
                 psFree(stamp->image1);
                 psFree(stamp->image2);
-                psFree(stamp->variance);
+                psFree(stamp->weight);
                 psFree(stamp->convolutions1);
                 psFree(stamp->convolutions2);
-                stamp->image1 = stamp->image2 = stamp->variance = NULL;
+                stamp->image1 = stamp->image2 = stamp->weight = NULL;
                 stamp->convolutions1 = stamp->convolutions2 = NULL;
 
@@ -487,5 +488,5 @@
                                                const psImage *image, const psImage *subMask,
                                                const psRegion *region, int size, int footprint,
-                                               float spacing, pmSubtractionMode mode)
+                                               float spacing, float sysErr, pmSubtractionMode mode)
 
 {
@@ -507,5 +508,6 @@
     int numStars = x->n;                // Number of stars
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
-                                                                 region, footprint, spacing); // Stamp list
+                                                                 region, footprint, sysErr,
+                                                                 spacing); // Stamp list
     int numStamps = stamps->num;        // Number of stamps
 
@@ -616,8 +618,10 @@
         PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
     }
-    PS_ASSERT_IMAGE_NON_NULL(variance, false);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false);
-    PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false);
-    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
+    if (variance) {
+        PS_ASSERT_IMAGE_NON_NULL(variance, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false);
+        PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false);
+        PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
+    }
 
     int numCols = image1->numCols, numRows = image1->numRows; // Size of images
@@ -646,5 +650,5 @@
         assert(stamp->image1 == NULL);
         assert(stamp->image2 == NULL);
-        assert(stamp->variance == NULL);
+        assert(stamp->weight == NULL);
 
         psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
@@ -660,7 +664,29 @@
         }
 
-        psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp
-        stamp->variance = psKernelAllocFromImage(wtSub, size, size);
-        psFree(wtSub);                  // Drop reference
+        psKernel *weight = stamp->weight = psKernelAlloc(-size, size, -size, size); // Weight = 1/variance
+        if (variance) {
+            psImage *varSub = psImageSubset(variance, region); // Subimage with stamp
+            psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp
+            if (isfinite(stamps->sysErr) && stamps->sysErr > 0) {
+                float sysErr = 0.5 * stamps->sysErr; // Systematic error
+                psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images
+                for (int y = -size; y <= size; y++) {
+                    for (int x = -size; x <= size; x++) {
+                        weight->kernel[y][x] = 1.0 / (var->kernel[y][x] +
+                                                      sysErr * (image1->kernel[y][x] + image2->kernel[y][x]));
+                    }
+                }
+            } else {
+                for (int y = -size; y <= size; y++) {
+                    for (int x = -size; x <= size; x++) {
+                        weight->kernel[y][x] = 1.0 / var->kernel[y][x];
+                    }
+                }
+            }
+            psFree(var);
+            psFree(varSub);
+        } else {
+            psImageInit(weight->image, 1.0);
+        }
 
         stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
@@ -672,5 +698,5 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
                                                           const psImage *subMask, const psRegion *region,
-                                                          int size, int footprint, float spacing,
+                                                          int size, int footprint, float spacing, float sysErr,
                                                           pmSubtractionMode mode)
 {
@@ -702,5 +728,5 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
-                                                            footprint, spacing, mode); // Stamps to return
+                                                            footprint, spacing, sysErr, mode); // Stamps
     psFree(x);
     psFree(y);
@@ -716,5 +742,5 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
                                                        const psImage *subMask, const psRegion *region,
-                                                       int size, int footprint, float spacing,
+                                                       int size, int footprint, float spacing, float sysErr,
                                                        pmSubtractionMode mode)
 {
@@ -734,5 +760,5 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
-                                                            spacing, mode);
+                                                            spacing, sysErr, mode);
     psFree(data);