Changeset 29543

Timestamp:

Oct 25, 2010, 2:45:41 PM (16 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20100823

Location:

trunk/psModules/src/imcombine

Files:

• README (copied) (copied from branches/eam_branches/ipp-20100823/psModules/src/imcombine/README )
• pmPSFEnvelope.c (modified) (1 diff)
• pmStack.c (modified) (30 diffs)
• pmSubtraction.c (modified) (1 diff)
• pmSubtractionAnalysis.c (modified) (1 diff)
• pmSubtractionAnalysis.h (modified) (1 diff)
• pmSubtractionEquation.c (modified) (47 diffs)
• pmSubtractionEquation.h (modified) (1 diff)
• pmSubtractionEquation.v0.c (copied) (copied from branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.v0.c )
• pmSubtractionKernels.c (modified) (6 diffs)
• pmSubtractionKernels.h (modified) (1 diff)
• pmSubtractionMatch.c (modified) (19 diffs)
• pmSubtractionMatch.h (modified) (1 diff)
• pmSubtractionStamps.c (modified) (23 diffs)
• pmSubtractionStamps.h (modified) (2 diffs)
• pmSubtractionVisual.c (modified) (9 diffs)

trunk/psModules/src/imcombine/pmPSFEnvelope.c

@@ -380 +380 @@
 
         // measure the source moments: tophat windowing, no pixel S/N cutoff
-        // XXX probably should be passing the maskVal to this function so we can pass it along here...
-        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0, maskVal)) {
+        if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, maskVal)) {
             // Can't do anything about it; limp along as best we can
             psErrorClear();

trunk/psModules/src/imcombine/pmStack.c

@@ -19 +19 @@
 
 #include <stdio.h>
+#include <string.h> // for memset
 #include <pslib.h>
 
@@ -34 +35 @@
 
 
-//#define TESTING                         // Enable test output
-//#define TEST_X 843-1                     // x coordinate to examine
-//#define TEST_Y 813-1                     // y coordinate to examine
-//#define TEST_RADIUS 0                    // Radius to examine
-
+# if (0)
+#define TESTING                         // Enable test output
+#define TEST_X 4963                       // x coordinate to examine
+#define TEST_Y 5353                       // y coordinate to examine
+#define TEST_X 40                       // x coordinate to examine
+#define TEST_Y 40                       // y coordinate to examine
+#define TEST_RADIUS 0.5                 // Radius to examine
+# endif
+
+# ifdef TESTING
+# define CHECKPIX(XPIX,YPIX,MSG,...) { if (PS_SQR(XPIX - TEST_X) + PS_SQR(YPIX - TEST_Y) <= PS_SQR(TEST_RADIUS)) { fprintf(stderr,MSG,__VA_ARGS__); } }
+# else
+# define CHECKPIX(XPIX,YPIX,MSG,...) { }
+# endif    
 
 // Data structure for use as a buffer when combining pixels
@@ -250 +260 @@
                                       )
 {
-#ifdef TESTING
-    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-        fprintf(stderr, "Marking image %d, pixel %d,%d for inspection\n", source, x, y);
-    }
-#endif
+    CHECKPIX(x, y, "Marking image %d, pixel %d,%d for inspection\n", source, x, y);
     pmStackData *data = inputs->data[source]; // Stack data of interest
     if (!data) {
@@ -272 +278 @@
                                      )
 {
-#ifdef TESTING
-    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-        fprintf(stderr, "Marking pixel image %d, pixel %d,%d for rejection\n", source, x, y);
-    }
-#endif
+    CHECKPIX(x, y, "Marking pixel image %d, pixel %d,%d for rejection\n", source, x, y);
     pmStackData *data = inputs->data[source]; // Stack data of interest
     if (!data) {
@@ -290 +292 @@
 static void combineExtract(int *num,                        // Number of good pixels
                            bool *suspect,                   // Any suspect pixels?
+                           psImageMaskType *badMask,        // OR of all bad (masked) pixels
+                           psImageMaskType *goodMask,       // OR of all good (unmasked) pixels
                            combineBuffer *buffer, // Buffer with vectors
                            psImage *image, // Combined image, for output
@@ -300 +304 @@
                            const psVector *reject, // Indices of pixels to reject, or NULL
                            int x, int y, // Coordinates of interest; frame of output image
-                           psImageMaskType maskVal, // Value to mask
-                           psImageMaskType maskSuspect // Value to suspect
+                           psImageMaskType badMaskBits, // Value to mask as 'bad'
+                           psImageMaskType suspectMaskBits // Value to mask as 'suspect'
                            )
 {
@@ -322 +326 @@
     }
 
+    // mask values to store possible mask combinations
+    *badMask = 0;
+    *goodMask = 0xffff;
+
+    int nGoodBits[16]; // accumulate the good pixel bits here for fuzzy logic
+    psAssert (sizeof(psImageMaskType) == 2, "psImageMaskType is not the expected size");
+    memset (nGoodBits, 0, 16*sizeof(int));
+
     // Extract the pixel and mask data
     int numGood = 0;                    // Number of good pixels
@@ -328 +340 @@
         // should be because of how pixelMapGenerate works
         if (reject && reject->data.U16[j] == i) {
+            // pixels can be rejected because:
+            // 1) only 1 input pixel (and 'safe' is true)
+            // 2) only 2 input pixels were available and they were mutually inconsistent (or variance info was missing)
+            // 3) NOTE : ifdef'ed out code for 3 inputs case 
+            // 4) outlier from sample of N pixels
+            // 5) some of these may have been suspect.
+            // XXX raise a specific mask bit for these (currently results in BLANK)
             j++;
+
+            pmStackData *data = inputs->data[i]; // Stack data of interest
+            if (data) {
+                int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
+                psImage *mask = data->readout->mask; // Mask of interest
+                *badMask |= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the bad bits used
+                CHECKPIX(x, y, "reject: adding bit to mask: %d : %x (badMask = %x)\n", i, mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn], *badMask);
+            } else {
+                CHECKPIX(x, y, "reject: no item in data (badMask = %x)\n", *badMask);
+            }
             continue;
         }
@@ -334 +363 @@
         pmStackData *data = inputs->data[i]; // Stack data of interest
         if (!data) {
+            CHECKPIX(x, y, "skip: no item in data (badMask = %x)\n", *badMask);
             continue;
         }
@@ -339 +369 @@
         int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
         psImage *mask = data->readout->mask; // Mask of interest
-        if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal) {
+        if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits) {
+            *badMask |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits); // save the bad bits used
+            CHECKPIX(x, y, "skip: adding bit to mask: %d : %x (badMask = %x)\n", i, mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn], *badMask);
             continue;
         }
 
-        pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ?
-            true : false;
+        pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & suspectMaskBits ? true : false;
+
+        // *goodMask &= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the mask bits still used
+        // check for set bits and increment counter as appropriate
+        { 
+            psImageMaskType value = 0x0001;
+            for (int nbit = 0; nbit < 16; nbit ++) {
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) {
+                    nGoodBits[nbit] ++;
+                }
+                value <<= 1;
+            }
+        }
 
         psImage *image = data->readout->image; // Image of interest
@@ -356 +399 @@
         pixelSources->data.U16[numGood] = i;
         numGood++;
+
+        CHECKPIX(x, y, "keep: %d : %x (badMask = %x)\n", i, mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn], *badMask);
     }
     pixelData->n = numGood;
@@ -367 +412 @@
     *num = numGood;
 
+    // set the mask bits if nGoodBits[i] > f*numGood
+    {
+        # define SUSPECT_FRACTION 0.65
+        *goodMask = 0x0000;
+        psImageMaskType value = 0x0001;
+        for (int nbit = 0; nbit < 16; nbit ++) {
+            if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) {
+                *goodMask |= value;
+            }
+            value <<= 1;
+        }
+    }
+
 #ifdef TESTING
     if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
         for (int i = 0; i < numGood; i++) {
-            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d\n",
+            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
                     i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
                     addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
-                    pixelSuspects->data.U8[i]);
+                    pixelSuspects->data.U8[i], badMaskBits, suspectMaskBits, *badMask, *goodMask);
         }
     }
@@ -380 +438 @@
     return;
 }
-
 
 // Combine pixels
@@ -392 +449 @@
                           combineBuffer *buffer, // Buffer with vectors
                           int x, int y, // Coordinates of interest; frame of output image
-                          psImageMaskType bad, // Value for bad pixels
+                          psImageMaskType blankMask, // Value for empty pixels
+                          psImageMaskType badMask, // Value for bad pixels
+                          psImageMaskType goodMask, // Value for good pixels
                           bool safe,           // Safe combination?
                           float invTotalWeight    // Inverse of total weight for all inputs
@@ -404 +463 @@
     // Default option is that the pixel is bad
     float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
-    psImageMaskType maskValue = bad;    // Value for combined mask
     float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted)
+
+    // default output mask value.  badMask is the OR of all unused input pixels.  
+    // if there are no input pixels, this value will be 0, in which case we want to set the output pixel to BLANK.
+    // if there are only good input pixels, and they do not result in a valid pixel, we still want to set this to BLANK.
+    psImageMaskType maskValue = badMask ? badMask : blankMask;    // Value for combined mask 
+
+    CHECKPIX(x, y, "bad vs good : %x %x %x\n", maskValue, badMask, blankMask);
 
     switch (num) {
       case 0: {
           // Nothing to combine: it's bad
-#ifdef TESTING
-          if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-              fprintf(stderr, "No inputs to combine, pixel %d,%d is bad.\n", x, y);
-          }
-#endif
+          CHECKPIX(x, y, "No inputs to combine, pixel %d,%d is bad.\n", x, y);
           break;
       }
@@ -428 +489 @@
                   expWeightValue = pixelExps->data.F32[0];
               }
-              maskValue = 0;
-#ifdef TESTING
-              if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                  fprintf(stderr, "Single input to combine, safety off, pixel %d,%d --> %f\n",
-                          x, y, imageValue);
-              }
-#endif
-          }
-#ifdef TESTING
-          else if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-              fprintf(stderr, "Single input to combine, safety on, pixel %d,%d is bad.\n", x, y);
-          }
-#endif
+              maskValue = goodMask;
+              CHECKPIX(x, y, "Single input to combine, safety off, pixel %d,%d --> %f\n", x, y, imageValue);
+          } else {
+              CHECKPIX(x, y, "Single input to combine, safety on, pixel %d,%d is bad.\n", x, y);
+          }
           break;
       }
@@ -446 +499 @@
           // Automatically accept the mean of the pixels only if we're not playing safe
           if (!safe) {
-              if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
-                                          pixelData, pixelVariances, pixelExps, pixelWeights)) {
-                  maskValue = 0;
-#ifdef TESTING
-                  if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                      fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n",
-                              x, y, imageValue, varianceValue);
-                  }
-#endif
+              if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue, pixelData, pixelVariances, pixelExps, pixelWeights)) {
+                  maskValue = goodMask;
+                  CHECKPIX(x, y, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
               }
+          } else {
+              CHECKPIX(x, y, "Two inputs to combine, safety on, pixel %d,%d is bad\n", x, y);
           }
-#ifdef TESTING
-          else {
-              if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                  fprintf(stderr, "Two inputs to combine, safety on, pixel %d,%d is bad\n", x, y);
-              }
-          }
-#endif
           break;
       }
@@ -472 +514 @@
               break;
           }
-          maskValue = 0;
-#ifdef TESTING
-          if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-              fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
-          }
-#endif
+          maskValue = goodMask;
+          CHECKPIX(x, y, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
           break;
       }
@@ -522 +560 @@
     int numIter = PS_MAX(iter * num, 1); // Number of iterations
 
-#ifdef TESTING
-    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-        fprintf(stderr, "Testing pixel %d,%d: %d %f %f %f %d %d\n",
-                x, y, numIter, rej, sys, olympic, useVariance, safe);
-    }
-#endif
+    CHECKPIX(x, y, "Testing pixel %d,%d: %d %f %f %f %d %d\n", x, y, numIter, rej, sys, olympic, useVariance, safe);
 
     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
@@ -548 +581 @@
             // Systematic error contributes to the rejection level
             float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
-#ifdef TESTING
-            // Correct variance for comparison against weighted mean including itself
-            float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
-            if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                fprintf(stderr, "Variance %d (%d), pixel %d,%d: %f %f %f\n", i, pixelSources->data.U16[i],
-                        x, y, pixelVariances->data.F32[i], sysVar, compare);
-            }
-#endif
+            CHECKPIX(x, y, "Variance %d (%d), pixel %d,%d: %f %f %f\n", 
+                     i, pixelSources->data.U16[i], x, y, 
+                     pixelVariances->data.F32[i], sysVar, 1.0 - pixelWeights->data.F32[i] / sumWeights);
             pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
         }
@@ -593 +621 @@
                           combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
                       }
-#ifdef TESTING
-                      if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                          fprintf(stderr, "Flagged both inputs for pixel %d,%d (%f > %f x %f\n)",
-                                  x, y, diff, rej, sqrtf(sigma2));
-                      }
-#endif
+                      CHECKPIX(x, y, "Flagged both inputs for pixel %d,%d (%f > %f x %f\n)", x, y, diff, rej, sqrtf(sigma2));
                   }
               } else if (i == 0 && safe) {
@@ -708 +731 @@
                   float median = combinationWeightedOlympic(pixelData, pixelWeights,
                                                             olympic, buffer->sort); // Median for stack
-#ifdef TESTING
-                  if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                      fprintf(stderr, "Flag with variance pixel %d,%d: median = %f\n", x, y, median);
-                  }
-#endif
+                  CHECKPIX(x, y, "Flag with variance pixel %d,%d: median = %f\n", x, y, median);
                   float worst = -INFINITY; // Largest deviation
                   for (int j = 0; j < num; j++) {
                       float diff = pixelData->data.F32[j] - median; // Difference from expected
-#ifdef TESTING
-                      if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                          fprintf(stderr, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff);
-                      }
-#endif
+                      CHECKPIX(x, y, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff);
 
                       // Comparing squares --- cheaper than lots of sqrts
@@ -737 +752 @@
                   combinationMedianStdev(&median, &stdev, pixelData, buffer->sort);
                   float limit = rej * stdev; // Rejection limit
-#ifdef TESTING
-                  if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                      fprintf(stderr,
-                              "Flag without variance pixel %d,%d; median = %f, stdev = %f, limit = %f\n",
-                              x, y, median, stdev, limit);
-                  }
-#endif
+                  CHECKPIX(x, y, "Flag without variance pixel %d,%d; median = %f, stdev = %f, limit = %f\n", x, y, median, stdev, limit);
                   float worst = -INFINITY; // Largest deviation
                   for (int j = 0; j < num; j++) {
@@ -763 +772 @@
         if (maskIndex >= 0) {
             if (suspect) {
-#ifdef TESTING
-                if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                    fprintf(stderr, "Throwing out all suspect pixels for %d,%d\n", x, y);
-                }
-#endif
+                CHECKPIX(x, y, "Throwing out all suspect pixels for %d,%d\n", x, y);
                 // Throw out all suspect pixels
                 int numGood = 0;        // Number of good pixels
@@ -796 +801 @@
             } else {
                 // Throw out masked pixel
-#ifdef TESTING
-                if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                    fprintf(stderr, "Throwing out input %d for pixel %d,%d\n", maskIndex, x, y);
-                }
-#endif
+                CHECKPIX(x, y, "Throwing out input %d for pixel %d,%d\n", maskIndex, x, y);
                 combineMarkInspect(inputs, x, y, pixelSources->data.U16[maskIndex]);
                 int numGood = 0;        // Number of good pixels
@@ -999 +1000 @@
 
 /// Stack input images
-bool pmStackCombine(pmReadout *combined, pmReadout *expmaps, psArray *input,
-                    psImageMaskType maskVal, psImageMaskType maskSuspect,
-                    psImageMaskType bad, int kernelSize,
-                    float iter, float rej, float sys, float olympic,
-                    bool useVariance, bool safe, bool rejection)
+bool pmStackCombine(
+    pmReadout *combined,                // output stacked readout
+    pmReadout *expmaps,                 // output exposure map information
+    psArray *input,                     // input exposures
+    psImageMaskType badMaskBits,        // treat these bits as 'bad'
+    psImageMaskType suspectMaskBits,    // treat these bits as 'suspect'
+    psImageMaskType blankMaskBits,      // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
+    int kernelSize,
+    float iter, 
+    float rej, 
+    float sys, 
+    float olympic,
+    bool useVariance, 
+    bool safe, 
+    bool rejection)
 {
     bool haveVariances;                 // Do we have the variance maps?
@@ -1012 +1023 @@
     }
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
-    PS_ASSERT_INT_POSITIVE(bad, false);
+    PS_ASSERT_INT_POSITIVE(blankMaskBits, false);
     if (isnan(rej)) {
         PS_ASSERT_FLOAT_EQUAL(iter, 0, false);
@@ -1070 +1081 @@
     }
     psFree(stack);
-    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, bad);
+    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, blankMaskBits);
     psTrace("psModules.imcombine", 1, "Have for combination [%d:%d,%d:%d] (%dx%d)\n",
             minInputCols, maxInputCols, minInputRows, maxInputRows, xSize, ySize);
@@ -1131 +1142 @@
     for (int y = minInputRows; y < maxInputRows; y++) {
         for (int x = minInputCols; x < maxInputCols; x++) {
+            CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection);
+
 #ifdef TESTING
-            if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                fprintf(stderr, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n",
-                        x, y, maskVal, bad, iter, rej, sys, olympic, useVariance, safe, rejection);
-            }
-#endif
+            if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
+                fprintf (stderr, "combine\n");
+            }
+# endif
+
             psVector *reject = NULL; // Images to reject for this pixel
             if (rejection) {
@@ -1154 +1167 @@
 #endif
             }
-
-            int num;                    // Number of good pixels
-            bool suspect;               // Suspect pixels in stack?
-            combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
-                           input, weights, exps, addVariance, reject, x, y, maskVal, maskSuspect);
-            combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight,
-                          num, buffer, x, y, bad, safe, totalExpWeight);
+ 
+            int num;                      // Number of good pixels
+            bool suspect;                 // Suspect pixels in stack?
+            psImageMaskType badMask = 0;  // OR of mask bits in all bad input pixels
+            psImageMaskType goodMask = 0; // OR of mask bits in all good input pixels
+            combineExtract(&num, &suspect, &badMask, &goodMask, buffer, combinedImage, combinedMask, combinedVariance, input, weights, exps, addVariance, reject, x, y, badMaskBits, suspectMaskBits);
+            combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight, num, buffer, x, y, blankMaskBits, badMask, goodMask, safe, totalExpWeight);
 
             if (iter > 0) {

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -540 +540 @@
 
     psImage *conv = psImageConvolveFFT(NULL, image->image, NULL, 0, kernel); // Convolved image
+
+    // note: do not attempt to renormalize kernels here: cannot have different stars with
+    // different kernel ratios
+
     int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
     psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version

trunk/psModules/src/imcombine/pmSubtractionAnalysis.c

@@ -299 +299 @@
                          kernels->rms);
 
-        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
-        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
-        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
-        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
-        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
-        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
-
-        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
-        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
-        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
-        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
-        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
-        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
+
+        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
+        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
+        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
+        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
+        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
+        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
     }
 

trunk/psModules/src/imcombine/pmSubtractionAnalysis.h

@@ -24 +24 @@
 #define PM_SUBTRACTION_ANALYSIS_DECONV_MAX   "SUBTRACTION.DECONV.MAX"   // Maximum deconvolution fraction
 
-#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN  "SUBTRACTION.RES.FSIGMA.MEAN"      // RMS stamp deviation
-#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV "SUBTRACTION.RES.FSIGMA.STDEV"      // RMS stamp deviation
-#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN    "SUBTRACTION.RES.FMIN.MEAN"      // RMS stamp deviation
-#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV   "SUBTRACTION.RES.FMIN.STDEV"      // RMS stamp deviation
-#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN    "SUBTRACTION.RES.FMAX.MEAN"      // RMS stamp deviation
-#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV   "SUBTRACTION.RES.FMAX.STDEV"      // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN  "SUBTRACTION.FRES.MEAN" // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV "SUBTRACTION.FRES.STDEV" // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN  "SUBTRACTION.FRES.OUTER.MEAN" // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV "SUBTRACTION.FRES.OUTER.STDEV" // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN  "SUBTRACTION.FRES.TOTAL.MEAN" // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV "SUBTRACTION.FRES.TOTAL.STDEV" // RMS stamp deviation
 
 // Derive QA information about the subtraction

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -19 +19 @@
 
 //#define USE_WEIGHT                      // Include weight (1/variance) in equation?
-//#define USE_WINDOW                      // Include weight (1/variance) in equation?
-
+
+// XXX TEST:
+//# define USE_WINDOW                      // window to avoid neighbor contamination
+
+# define PENALTY false
+# define MOMENTS (!PENALTY)
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -27 +31 @@
 
 // Calculate the least-squares matrix and vector
-static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
-                                  psVector *vector, // Least-squares vector, updated
-                                  double *norm,     // Normalisation, updated
+static bool calculateMatrixVector(psImage *matrix,       // Least-squares matrix, updated
+                                  psVector *vector,      // Least-squares vector, updated
+                                  double normValue,      // Normalisation, supplied
                                   const psKernel *input, // Input image (target)
                                   const psKernel *reference, // Reference image (convolution source)
@@ -37 +41 @@
                                   const pmSubtractionKernels *kernels, // Kernels
                                   const psImage *polyValues, // Spatial polynomial values
-                                  int footprint, // (Half-)Size of stamp
-                                  int normWindow1, // Window (half-)size for normalisation measurement
-                                  int normWindow2, // Window (half-)size for normalisation measurement
-                                  const pmSubtractionEquationCalculationMode mode
+                                  int footprint // (Half-)Size of stamp
                                   )
 {
@@ -50 +51 @@
 
     int numKernels = kernels->num;                      // Number of kernels
-    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-    int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     double poly[numPoly];                                 // Polynomial terms
     double poly2[numPoly][numPoly];                       // Polynomial-polynomial values
+    int numParams = numKernels * numPoly;
+
+    psAssert(matrix &&
+             matrix->type.type == PS_TYPE_F64 &&
+             matrix->numCols == numParams &&
+             matrix->numRows == numParams,
+             "Least-squares matrix is bad.");
+    psAssert(vector &&
+             vector->type.type == PS_TYPE_F64 &&
+             vector->n == numParams,
+             "Least-squares vector is bad.");
 
     // Evaluate polynomial-polynomial terms
-    // XXX we can skip this if we are not calculating kernel coeffs
     for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) {
         for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) {
@@ -84 +93 @@
     // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0]
     // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m]
-    // normalization
-    // bg 0, bg 1, bg 2 (only 0 is currently used?)
 
     for (int i = 0; i < numKernels; i++) {
@@ -107 +114 @@
 
             // Spatial variation of kernel coeffs
-            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
-                        double value = sumCC * poly2[iTerm][jTerm];
-                        matrix->data.F64[iIndex][jIndex] = value;
-                        matrix->data.F64[jIndex][iIndex] = value;
-                    }
-                }
-            }
+            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                    double value = sumCC * poly2[iTerm][jTerm];
+                    matrix->data.F64[iIndex][jIndex] = value;
+                    matrix->data.F64[jIndex][iIndex] = value;
+                }
+            }
         }
 
         double sumRC = 0.0;             // Sum of the reference-convolution products
         double sumIC = 0.0;             // Sum of the input-convolution products
-        double sumC = 0.0;              // Sum of the convolution
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
@@ -128 +132 @@
                 double ic = in * conv;
                 double rc = ref * conv;
-                double c = conv;
                 if (weight) {
                     float wtVal = weight->kernel[y][x];
                     ic *= wtVal;
                     rc *= wtVal;
-                    c *= wtVal;
                 }
                 if (window) {
@@ -139 +141 @@
                     ic *= winVal;
                     rc *= winVal;
-                    c  *= winVal;
                 }
                 sumIC += ic;
                 sumRC += rc;
-                sumC += c;
-            }
-        }
+            }
+        }
+
         // Spatial variation
         for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-            double normTerm = sumRC * poly[iTerm];
-            double bgTerm = sumC * poly[iTerm];
-            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
-                matrix->data.F64[iIndex][normIndex] = normTerm;
-                matrix->data.F64[normIndex][iIndex] = normTerm;
-            }
-            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
-                matrix->data.F64[iIndex][bgIndex] = bgTerm;
-                matrix->data.F64[bgIndex][iIndex] = bgTerm;
-            }
-            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-                vector->data.F64[iIndex] = sumIC * poly[iTerm];
-                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
-                    // subtract norm * sumRC * poly[iTerm]
-                    psAssert (kernels->solution1, "programming error: define solution first!");
-                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
-                    vector->data.F64[iIndex] -= norm * normTerm;
-                }
-            }
-        }
-    }
-
-    double sumRR = 0.0;                 // Sum of the reference product
-    double sumIR = 0.0;                 // Sum of the input-reference product
-    double sum1 = 0.0;                  // Sum of the background
-    double sumR = 0.0;                  // Sum of the reference
-    double sumI = 0.0;                  // Sum of the input
-    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
-    for (int y = - footprint; y <= footprint; y++) {
-        for (int x = - footprint; x <= footprint; x++) {
-            double in = input->kernel[y][x];
-            double ref = reference->kernel[y][x];
-            double ir = in * ref;
-            double rr = PS_SQR(ref);
-            double one = 1.0;
-
-            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) {
-                normI1 += ref;
-            }
-            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) {
-                normI2 += in;
-            }
-
-            if (weight) {
-                float wtVal = weight->kernel[y][x];
-                rr *= wtVal;
-                ir *= wtVal;
-                in *= wtVal;
-                ref *= wtVal;
-                one *= wtVal;
-            }
-            if (window) {
-                float  winVal = window->kernel[y][x];
-                rr      *= winVal;
-                ir      *= winVal;
-                in      *= winVal;
-                ref *= winVal;
-                one *= winVal;
-            }
-            sumRR += rr;
-            sumIR += ir;
-            sumR += ref;
-            sumI += in;
-            sum1 += one;
-        }
-    }
-
-    *norm = normI2 / normI1;
-
-    fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm);
-
-    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
-        matrix->data.F64[normIndex][normIndex] = sumRR;
-        vector->data.F64[normIndex] = sumIR;
-        // subtract sum over kernels * kernel solution
-    }
-    if (mode & PM_SUBTRACTION_EQUATION_BG) {
-        matrix->data.F64[bgIndex][bgIndex] = sum1;
-        vector->data.F64[bgIndex] = sumI;
-    }
-    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
-        matrix->data.F64[normIndex][bgIndex] = sumR;
-        matrix->data.F64[bgIndex][normIndex] = sumR;
+            vector->data.F64[iIndex] = (sumIC - normValue*sumRC) * poly[iTerm];
+        }
     }
 
@@ -255 +174 @@
 
 // Calculate the least-squares matrix and vector for dual convolution
-static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated
-                                      psVector *vector, // Least-squares vector, updated
-                                      double *norm,     // Normalisation, updated
-                                      const psKernel *image1, // Image 1
-                                      const psKernel *image2, // Image 2
+// XXX we could avoid calculating these values on successive passes *if* the stamp has not changed.
+static bool calculateDualMatrixVector(pmSubtractionStamp *stamp,              // stamp of interest
+                                      double normValue,       // Normalisation, updated
+                                      double normValue2,      // Normalisation, updated
                                       const psKernel *weight,  // Weight image
                                       const psKernel *window,  // Window image
-                                      const psArray *convolutions1, // Convolutions of image 1 for each kernel
-                                      const psArray *convolutions2, // Convolutions of image 2 for each kernel
                                       const pmSubtractionKernels *kernels, // Kernels
                                       const psImage *polyValues, // Spatial polynomial values
-                                      int footprint, // (Half-)Size of stamp
-                                      int normWindow1, // Window (half-)size for normalisation measurement
-                                      int normWindow2, // Window (half-)size for normalisation measurement
-                                      const pmSubtractionEquationCalculationMode mode
+                                      int footprint // (Half-)Size of stamp
                                       )
 {
-    int numKernels = kernels->num;                      // Number of kernels
-    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-    int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
-    int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
+    int numKernels = kernels->num;                        // Number of kernels
+    int spatialOrder = kernels->spatialOrder;             // Order of spatial variation
     int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     double poly[numPoly];                                 // Polynomial terms
     double poly2[numPoly][numPoly];                       // Polynomial-polynomial values
 
-    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
-    int numParams = numKernels * numPoly + 1 + numBackground;       // Number of regular parameters
-    int numParams2 = numKernels * numPoly;                          // Number of additional parameters for dual
-    int numDual = numParams + numParams2;                           // Total number of parameters for dual
+    int numParams = numKernels * numPoly;                 // Number of regular parameters
+    int numParams2 = numKernels * numPoly;         // Number of additional parameters for dual
+    int numDual = numParams + numParams2;          // Total number of parameters for dual
+
+    psImage *matrix = stamp->matrix;               // Least-squares matrix, updated
+    psVector *vector = stamp->vector;              // Least-squares vector, updated
+    psKernel *image1 = stamp->image1;              // Image 1
+    psKernel *image2 = stamp->image2;              // Image 2
+    psArray *convolutions1 = stamp->convolutions1; // Convolutions of image 1 for each kernel
+    psArray *convolutions2 = stamp->convolutions2; // Convolutions of image 2 for each kernel
 
     psAssert(matrix &&
@@ -290 +207 @@
              matrix->numRows == numDual,
              "Least-squares matrix is bad.");
+
     psAssert(vector &&
              vector->type.type == PS_TYPE_F64 &&
@@ -318 +236 @@
     }
 
+    // the order of the elements in the matrix and vector is:
+    // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0]
+    // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0]
+    // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m]
+
+    // for DUAL convolution analysis, we apply the normalization to I1 as follows:
+    // norm = I2 / I1
+    // 
+    // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i
+    // I2c =      I2 + \Sum_i b_i      I2 \cross k_i
+
+    // we cannot absorb the normalization into a_i until the analysis is complete, or the
+    // second moment terms are incorrectly calculated.
+
     for (int i = 0; i < numKernels; i++) {
         psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i
@@ -328 +260 @@
             double sumBB = 0.0;         // Sum of convolution products between image 2
             double sumAB = 0.0;         // Sum of convolution products across images 1 and 2
+
+            double MxxAA = 0.0;
+            double MyyAA = 0.0;
+            double MxxBB = 0.0;
+            double MyyBB = 0.0;
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
-                    double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x];
+                    double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue);
                     double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
-                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
+                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue;
                     if (weight) {
                         float wtVal = weight->kernel[y][x];
@@ -348 +285 @@
                     sumBB += bb;
                     sumAB += ab;
-                }
-            }
+
+                    if (MOMENTS) {
+                        MxxAA += x*x*aa;
+                        MyyAA += y*y*aa;
+                        MxxBB += x*x*bb;
+                        MyyBB += y*y*bb;
+                    }
+                }
+            }
+
+            if (MOMENTS) {
+                MxxAA /= stamp->normSquare1 * PS_SQR(normValue);
+                MyyAA /= stamp->normSquare1 * PS_SQR(normValue);
+                MxxBB /= stamp->normSquare2;
+                MyyBB /= stamp->normSquare2;
+            }
+
+            // XXX this makes the Chisq portion independent of the normalization and star flux
+            // but may be mis-scaling between stars of different fluxes
+            sumAA /= PS_SQR(stamp->normI2);
+            sumAB /= PS_SQR(stamp->normI2);
+            sumBB /= PS_SQR(stamp->normI2);
+
+            // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB);
 
             // Spatial variation of kernel coeffs
-            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
-                        double aa = sumAA * poly2[iTerm][jTerm];
-                        double bb = sumBB * poly2[iTerm][jTerm];
-                        double ab = sumAB * poly2[iTerm][jTerm];
-
-                        matrix->data.F64[iIndex][jIndex] = aa;
-                        matrix->data.F64[jIndex][iIndex] = aa;
-
-                        matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
-                        matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
-
-                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
-                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
-                    }
-                }
-            }
-        }
+            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                    double aa = sumAA * poly2[iTerm][jTerm];
+                    double bb = sumBB * poly2[iTerm][jTerm];
+                    double ab = sumAB * poly2[iTerm][jTerm];
+
+                    matrix->data.F64[iIndex][jIndex] = aa;
+                    matrix->data.F64[jIndex][iIndex] = aa;
+
+                    matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
+                    matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
+
+                    // add in second moments
+                    if (MOMENTS) {
+                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA;
+                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA;
+
+                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA;
+                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA;
+
+                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB;
+                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB;
+
+                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB;
+                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB;
+
+                        // XXX this uses the single moments, first try
+                        // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
+                        // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
+                        // 
+                        // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
+                        // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
+                        // 
+                        // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
+                        // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
+                        // 
+                        // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
+                        // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
+                    }
+                    matrix->data.F64[iIndex][jIndex + numParams] = ab;
+                    matrix->data.F64[jIndex + numParams][iIndex] = ab;
+                }
+            }
+        }
+
+        // we need to calculate the lower-diagonal AB elements since they are not symmetric for A <-> B
         for (int j = 0; j < i; j++) {
             psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j
@@ -376 +362 @@
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
-                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
+                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue;
                     if (weight) {
                         ab *= weight->kernel[y][x];
@@ -387 +373 @@
             }
 
+            // XXX this makes the Chisq portion independent of the normalization and star flux
+            // but may be mis-scaling between stars of different fluxes
+            sumAB /= PS_SQR(stamp->normI2);
+
             // Spatial variation of kernel coeffs
-            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
-                        double ab = sumAB * poly2[iTerm][jTerm];
-                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
-                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
-                    }
-                }
+            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                    double ab = sumAB * poly2[iTerm][jTerm];
+                    matrix->data.F64[iIndex][jIndex + numParams] = ab;
+                    matrix->data.F64[jIndex + numParams][iIndex] = ab;
+                }
             }
         }
@@ -402 +390 @@
         double sumBI2 = 0.0;            // Sum of B.I_2 products (for vector)
         double sumAI1 = 0.0;            // Sum of A.I_1 products (for matrix, normalisation)
-        double sumA = 0.0;              // Sum of A (for matrix, background)
         double sumBI1 = 0.0;            // Sum of B.I_1 products (for matrix, normalisation)
-        double sumB = 0.0;              // Sum of B products (for matrix, background)
-        double sumI2 = 0.0;             // Sum of I_2 (for vector, background)
+
+        double MxxAI1 = 0.0;
+        double MyyAI1 = 0.0;
+        double MxxBI2 = 0.0;
+        double MyyBI2 = 0.0;
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
@@ -413 +403 @@
                 float i2 = image2->kernel[y][x];
 
-                double ai2 = a * i2;
+                double ai2 = a * i2 * normValue;
                 double bi2 = b * i2;
-                double ai1 = a * i1;
-                double bi1 = b * i1;
+                double ai1 = a * i1 * PS_SQR(normValue);
+                double bi1 = b * i1 * normValue;
 
                 if (weight) {
@@ -424 +414 @@
                     ai1 *= wtVal;
                     bi1 *= wtVal;
-                    a *= wtVal;
-                    b *= wtVal;
-                    i2 *= wtVal;
                 }
                 if (window) {
@@ -434 +421 @@
                     ai1 *= wtVal;
                     bi1 *= wtVal;
-                    a *= wtVal;
-                    b *= wtVal;
-                    i2 *= wtVal;
                 }
                 sumAI2 += ai2;
                 sumBI2 += bi2;
                 sumAI1 += ai1;
-                sumA += a;
                 sumBI1 += bi1;
-                sumB += b;
-                sumI2 += i2;
-            }
-        }
+
+                if (MOMENTS) {
+                    MxxAI1 += x*x*ai1;
+                    MyyAI1 += y*y*ai1;
+                    MxxBI2 += x*x*bi2;
+                    MyyBI2 += y*y*bi2;
+                }
+            }
+        }
+
+        if (MOMENTS) {
+            MxxAI1 /= stamp->normSquare1 * PS_SQR(normValue);
+            MyyAI1 /= stamp->normSquare1 * PS_SQR(normValue);
+            MxxBI2 /= stamp->normSquare2;
+            MyyBI2 /= stamp->normSquare2;
+        }
+
+        // fprintf (stderr, "i : %d : M(xx,yy)(AI1,BI2) : %f %f %f %f\n", i, MxxAI1, MyyAI1, MxxBI2, MyyBI2);
+
+        // XXX this makes the Chisq portion independent of the normalization and star flux
+        // but may be mis-scaling between stars of different fluxes
+        sumAI1 /= PS_SQR(stamp->normI2);
+        sumBI1 /= PS_SQR(stamp->normI2);
+        sumAI2 /= PS_SQR(stamp->normI2);
+        sumBI2 /= PS_SQR(stamp->normI2);
+
         // Spatial variation
         for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
@@ -452 +457 @@
             double bi2 = sumBI2 * poly[iTerm];
             double ai1 = sumAI1 * poly[iTerm];
-            double a   = sumA * poly[iTerm];
             double bi1 = sumBI1 * poly[iTerm];
-            double b   = sumB * poly[iTerm];
-
-            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
-                matrix->data.F64[iIndex][normIndex] = ai1;
-                matrix->data.F64[normIndex][iIndex] = ai1;
-                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
-                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
-            }
-            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
-                matrix->data.F64[iIndex][bgIndex] = a;
-                matrix->data.F64[bgIndex][iIndex] = a;
-                matrix->data.F64[iIndex + numParams][bgIndex] = b;
-                matrix->data.F64[bgIndex][iIndex + numParams] = b;
-            }
-            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-                vector->data.F64[iIndex] = ai2;
-                vector->data.F64[iIndex + numParams] = bi2;
-                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
-                    // subtract norm * sumRC * poly[iTerm]
-                    psAssert (kernels->solution1, "programming error: define solution first!");
-                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
-                    vector->data.F64[iIndex] -= norm * ai1;
-                    vector->data.F64[iIndex + numParams] -= norm * bi1;
-                }
-            }
-        }
-    }
-
-    double sumI1 = 0.0;                 // Sum of I_1 (for matrix, background-normalisation)
-    double sumI1I1 = 0.0;               // Sum of I_1^2 (for matrix, normalisation-normalisation)
-    double sum1 = 0.0;                  // Sum of 1 (for matrix, background-background)
-    double sumI2 = 0.0;                 // Sum of I_2 (for vector, background)
-    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
-    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
-    for (int y = - footprint; y <= footprint; y++) {
-        for (int x = - footprint; x <= footprint; x++) {
-            double i1 = image1->kernel[y][x];
-            double i2 = image2->kernel[y][x];
-
-            double i1i1 = i1 * i1;
-            double one = 1.0;
-            double i1i2 = i1 * i2;
-
-            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) {
-                normI1 += i1;
-            }
-            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) {
-                normI2 += i2;
-            }
-
-            if (weight) {
-                float wtVal = weight->kernel[y][x];
-                i1 *= wtVal;
-                i1i1 *= wtVal;
-                one *= wtVal;
-                i2 *= wtVal;
-                i1i2 *= wtVal;
-            }
-            if (window) {
-                float wtVal = window->kernel[y][x];
-                i1 *= wtVal;
-                i1i1 *= wtVal;
-                one *= wtVal;
-                i2 *= wtVal;
-                i1i2 *= wtVal;
-            }
-            sumI1 += i1;
-            sumI1I1 += i1i1;
-            sum1 += one;
-            sumI2 += i2;
-            sumI1I2 += i1i2;
-        }
-    }
-
-    *norm = normI2 / normI1;
-    fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm);
-
-    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
-        matrix->data.F64[normIndex][normIndex] = sumI1I1;
-        vector->data.F64[normIndex] = sumI1I2;
-    }
-    if (mode & PM_SUBTRACTION_EQUATION_BG) {
-        matrix->data.F64[bgIndex][bgIndex] = sum1;
-        vector->data.F64[bgIndex] = sumI2;
-    }
-    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
-        matrix->data.F64[bgIndex][normIndex] = sumI1;
-        matrix->data.F64[normIndex][bgIndex] = sumI1;
+            vector->data.F64[iIndex]             = ai2 - ai1;
+            vector->data.F64[iIndex + numParams] = bi2 - bi1;
+
+            // fprintf (stderr, "i : %d : V(I1,I2) : %f %f\n", i, vector->data.F64[iIndex], vector->data.F64[iIndex + numParams]);
+
+            // add in second moments
+            if (MOMENTS) {
+                vector->data.F64[iIndex]             -= kernels->penalty * MxxAI1;
+                vector->data.F64[iIndex]             -= kernels->penalty * MyyAI1;
+
+                vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2;
+                vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2;
+            }
+        }
     }
 
@@ -565 +494 @@
 }
 
-#if 1
 // Add in penalty term to least-squares vector
 bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
                       psVector *vector,                    // Vector to which to add in penalty term
                       const pmSubtractionKernels *kernels, // Kernel parameters
-                      float norm                           // Normalisation
+                      float normSquare1,                   // Normalisation for image 1
+                      float normSquare2            // Normalisation for image 2
   )
 {
+    psAssert (kernels->mode == PM_SUBTRACTION_MODE_DUAL, "only use penalties for dual convolution");
+
     if (kernels->penalty == 0.0) {
         return true;
@@ -589 +520 @@
     // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0]
     // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1]
-    // [norm]
-    // [bg]
+
     // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0]
     // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0]
@@ -600 +530 @@
                 // Contribution to chi^2: a_i^2 P_i
                 psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty");
-                fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]);
-                matrix->data.F64[index][index] += norm * penalties1->data.F32[i];
-                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-                    fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 
-                             matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]);
-                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i];                        
-                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
-                    // penalties scale with second moments
-                    //
-                }
-            }
-        }
-    }
-
-    return true;
-}
-# endif
+                // fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]);
+                matrix->data.F64[index][index] += normSquare1 * penalties1->data.F32[i];
+
+                // fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 
+                // matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]);
+                matrix->data.F64[index + numParams][index + numParams] += normSquare2 * penalties2->data.F32[i];                             
+                // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
+                // penalties scale with second moments
+                //
+            }
+        }
+    }
+
+    return true;
+}
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -647 +575 @@
                                     int index, bool wantDual)
 {
-#if 0
     // This is probably in a tight loop, so don't check inputs
-    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);
-    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);
-    PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);
-    PS_ASSERT_INT_POSITIVE(index, NAN);
-#endif
-
     psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector
     return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index,
@@ -662 +583 @@
 double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels)
 {
-#if 0
     // This is probably in a tight loop, so don't check inputs
-    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);
-    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);
-    PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);
-#endif
-
     int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     return kernels->solution1->data.F64[normIndex];
@@ -676 +591 @@
                                                 const psImage *polyValues)
 {
-#if 0
     // This is probably in a tight loop, so don't check inputs
-    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);
-    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);
-    PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);
-#endif
-
     int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1);
@@ -690 +599 @@
 // Public functions
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+bool pmSubtractionCalculateMoments(
+    pmSubtractionKernels *kernels, // Kernels
+    pmSubtractionStampList *stamps)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+
+    // XXX skip this, right?
+    return true;
+
+    // these are only used by DUAL mode
+    if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true;
+
+    psTimerStart("pmSubtractionCalculateMoments");
+
+    int footprint = stamps->footprint;  // Half-size of stamps
+
+    // Loop over each stamp and calculate its normalization factor 
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue;
+        if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue;
+
+        pmSubtractionCalculateMomentsStamp(kernels, stamp, footprint, stamps->normWindow2, stamps->normWindow1);
+    }
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate moments: %f sec", psTimerClear("pmSubtractionCalculateMoments"));
+
+    return true;
+}
+
+bool pmSubtractionCalculateMomentsStamp(
+    pmSubtractionKernels *kernels, // Kernels
+    pmSubtractionStamp *stamp,          // stamp on which to save normalization)
+    int footprint,                      // (Half-)Size of stamp
+    int normWindow1,                    // Window (half-)size for normalisation measurement
+    int normWindow2                     // Window (half-)size for normalisation measurement
+    )
+{
+    double Mxx, Myy;
+
+    int numKernels = kernels->num;
+
+    // Generate convolutions: these are generated once and saved
+    if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+        psError(psErrorCodeLast(), false, "Unable to convolve stamp");
+        return false;
+    }
+
+    if (!stamp->MxxI1) {
+        stamp->MxxI1 = psVectorAlloc (numKernels, PS_TYPE_F32);
+    }
+    if (!stamp->MyyI1) {
+        stamp->MyyI1 = psVectorAlloc (numKernels, PS_TYPE_F32);
+    }
+    if (!stamp->MxxI2) {
+        stamp->MxxI2 = psVectorAlloc (numKernels, PS_TYPE_F32);
+    }
+    if (!stamp->MyyI2) {
+        stamp->MyyI2 = psVectorAlloc (numKernels, PS_TYPE_F32);
+    }
+
+    for (int i = 0; i < numKernels; i++) {
+        pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions1->data[i], footprint, normWindow1);
+        stamp->MxxI1->data.F32[i] = Mxx / stamp->normI1;
+        stamp->MyyI1->data.F32[i] = Myy / stamp->normI1;
+        pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions2->data[i], footprint, normWindow2);
+        stamp->MxxI2->data.F32[i] = Mxx / stamp->normI2;
+        stamp->MyyI2->data.F32[i] = Myy / stamp->normI2;
+    }
+
+    pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image1, footprint, normWindow1);
+    stamp->MxxI1raw = Mxx / stamp->normI1;
+    stamp->MyyI1raw = Myy / stamp->normI1;
+
+    pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image2, footprint, normWindow2);
+    stamp->MxxI2raw = Mxx / stamp->normI2;
+    stamp->MyyI2raw = Myy / stamp->normI2;
+
+    // fprintf (stderr, "Mxx I1: %f, Myy I1: %f, Mxx I2: %f, Myy I2: %f\n", stamp->MxxI1raw, stamp->MyyI1raw, stamp->MxxI2raw, stamp->MyyI2raw);
+
+    return true;
+}
+
+bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window) {
+
+    double Sxx = 0.0;
+    double Syy = 0.0;
+    for (int y = - footprint; y <= footprint; y++) {
+        for (int x = - footprint; x <= footprint; x++) {
+            if (PS_SQR(x) + PS_SQR(y) > PS_SQR(window)) continue;
+
+            double flux = image->kernel[y][x];
+
+            Sxx += PS_SQR(x) * flux;
+            Syy += PS_SQR(y) * flux;
+        }
+    }
+    *Mxx = Sxx;
+    *Myy = Syy;
+    return true;
+}
+
+///---------
+
+bool pmSubtractionCalculateNormalization(
+    pmSubtractionStampList *stamps,
+    const pmSubtractionMode mode)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+
+    psTimerStart("pmSubtractionCalculateNormalization");
+
+    psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
+    psVector *norm2 = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
+
+    int footprint = stamps->footprint;  // Half-size of stamps
+
+    // Loop over each stamp and calculate its normalization factor 
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue;
+        if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue;
+
+        // XXX skip this if we have already calculated it? (stamp->norm does not change, just the median statistic)
+        // XXX maybe not: the star may have changed for a given stamp -- only if norm is reset to NAN can we do this
+        if (mode == PM_SUBTRACTION_MODE_2) {
+            pmSubtractionCalculateNormalizationStamp(stamp, stamp->image2, stamp->image1, footprint, stamps->normWindow2, stamps->normWindow1);
+        } else {
+            pmSubtractionCalculateNormalizationStamp(stamp, stamp->image1, stamp->image2, footprint, stamps->normWindow1, stamps->normWindow2);
+        }
+        psVectorAppend(norms, stamp->norm);
+        psVectorAppend(norm2, stamp->normSquare2 / stamp->normSquare1);
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
+    if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
+        psError(PM_ERR_DATA, false, "Unable to determine median normalisation");
+        psFree(stats);
+        psFree(norms);
+        psFree(norm2);
+        return false;
+    }
+    stamps->normValue = stats->robustMedian;
+
+    psStatsInit(stats);
+    if (!psVectorStats(stats, norm2, NULL, NULL, 0)) {
+        psError(PM_ERR_DATA, false, "Unable to determine median normalisation");
+        psFree(stats);
+        psFree(norms);
+        psFree(norm2);
+        return false;
+    }
+    stamps->normValue2 = stats->robustMedian;
+
+    psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2);
+
+    psFree(stats);
+    psFree(norms);
+    psFree(norm2);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate normalization: %f sec", psTimerClear("pmSubtractionCalculateNormalization"));
+
+    return true;
+}
+
+bool pmSubtractionCalculateNormalizationStamp(
+    pmSubtractionStamp *stamp,          // stamp on which to save normalization)
+    const psKernel *image1,             // Input image (target)
+    const psKernel *image2,             // Reference image (convolution source)
+    int footprint,                      // (Half-)Size of stamp
+    int normWindow1,                    // Window (half-)size for normalisation measurement
+    int normWindow2                     // Window (half-)size for normalisation measurement
+    )
+{
+    double normI1 = 0.0;  // Sum of I_1 within the normalisation window (aperture)
+    double normI2 = 0.0;  // Sum of I_2 within the normalisation window (aperture)
+    double normSquare1 = 0.0;  // Sum of (I_1)^2 within the normalisation window (aperture)
+    double normSquare2 = 0.0;  // Sum of (I_2)^2 within the normalisation window (aperture)
+    for (int y = - footprint; y <= footprint; y++) {
+        for (int x = - footprint; x <= footprint; x++) {
+            double im1 = image1->kernel[y][x];
+            double im2 = image2->kernel[y][x];
+
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) {
+                normI1 += im1;
+                normSquare1 += PS_SQR(im1);
+            }
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) {
+                normI2 += im2;
+                normSquare2 += PS_SQR(im2);
+            }
+        }
+    }
+    stamp->norm = normI2 / normI1;
+    stamp->normI1 = normI1;
+    stamp->normI2 = normI2;
+    stamp->normSquare1 = normSquare1;
+    stamp->normSquare2 = normSquare2;
+
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f  (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2);
+
+    return true;
+}
 
 bool pmSubtractionCalculateEquationThread(psThreadJob *job)
@@ -715 +828 @@
     int numKernels = kernels->num;      // Number of kernel basis functions
     int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations
-    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
 
     // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial).
@@ -722 +834 @@
     // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the
     // number of coefficients for the spatial polynomial, normalisation and a constant background offset.
-    int numParams = numKernels * numSpatial + 1 + numBackground;
+    // XXX we no longer solve for the normalization and background in the matrix inversion
+    int numParams = numKernels * numSpatial;
     if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         // An additional image is convolved
@@ -790 +903 @@
     switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_1:
-        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
-                                       weight, window, stamp->convolutions1, kernels,
-                                       polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
+        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image2, stamp->image1,
+                                       weight, window, stamp->convolutions1, kernels, polyValues, footprint);
         break;
       case PM_SUBTRACTION_MODE_2:
-        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
-                                       weight, window, stamp->convolutions2, kernels,
-                                       polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode);
+        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2,
+                                       weight, window, stamp->convolutions2, kernels, polyValues, footprint);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
-        status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm,
-                                           stamp->image1, stamp->image2,
-                                           weight, window, stamp->convolutions1, stamp->convolutions2,
-                                           kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
+        status = calculateDualMatrixVector(stamp, stamps->normValue, stamps->normValue2, weight, window, kernels, polyValues, footprint);
         break;
       default:
@@ -928 +1036 @@
     int numKernels = kernels->num;      // Number of kernel basis functions
     int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
-    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
-    int numParams = numKernels * numSpatial + 1 + numBackground;    // Number of parameters being solved for
-    int numSolution1 = numParams, numSolution2 = 0;                 // Number of parameters for each solution
+    // XXX int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
+    // XXX int numParams = numKernels * numSpatial + 1 + numBackground;    // Number of parameters being solved for
+    int numParams = numKernels * numSpatial;        // Number of parameters being solved for
+    int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution
     if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         // An additional image is convolved
@@ -960 +1069 @@
 
     if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
-        // Accumulate the least-squares matricies and vectors
+        // Accumulate the least-squares matrices and vectors.  These are generated for the
+        // kernel elements, excluding the background and normalization.
         psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix
         psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector
@@ -966 +1076 @@
         psImageInit(sumMatrix, 0.0);
 
-        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
-
         int numStamps = 0;              // Number of good stamps
+        for (int i = 0; i < stamps->num; i++) {
+            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
+                
+                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
+                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
+
+                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
+                numStamps++;
+            } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
+                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
+            }
+        }
+
+#if 0
+        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
+        psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
+        psVectorWriteFile("sumVector.dat", sumVector);
+        psFree (save);
+#endif
+
+        psVector *solution = NULL;                       // Solution to equation!
+        solution = psVectorAlloc(numParams, PS_TYPE_F64);
+        psVectorInit(solution, 0);
+
+        // XXX TEST: try some constraint on the svd solution
+        // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+        // SINGLE solution
+        if (1) {
+            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+        } else {
+            psVector *PERM = NULL;
+            psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix);
+            solution = psMatrixLUSolution(solution, LU, sumVector, PERM);
+            psFree (LU);
+            psFree (PERM);
+        }
+
+# if (0)
+        for (int i = 0; i < solution->n; i++) {
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]);
+        }
+# endif
+
+        if (!kernels->solution1) {
+            kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+            psVectorInit(kernels->solution1, 0.0);
+        }
+
+        int numKernels = kernels->num;
+        int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
+        int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+        for (int i = 0; i < numKernels * numPoly; i++) {
+            kernels->solution1->data.F64[i] = solution->data.F64[i];
+        }
+
+        // Apply the normalisation and background separately
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
+        kernels->solution1->data.F64[normIndex] = stamps->normValue;
+        kernels->solution1->data.F64[bgIndex] = 0.0;
+
+        psFree(solution);
+        psFree(sumVector);
+        psFree(sumMatrix);
+
+    } else {
+        // Dual convolution solution
+        // Accumulate the least-squares matrices and vectors.  These are generated for the
+        // kernel elements, excluding the background and normalization.
+        psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
+        psImageInit(sumMatrix, 0.0);
+        psVectorInit(sumVector, 0.0);
+
+        int numStamps = 0;         // Number of good stamps
+        double normSquare1 = 0.0; // Sum of (I_1)^2 over stamps
+        double normSquare2 = 0.0; // Sum of (I_2)^2 over stamps
         for (int i = 0; i < stamps->num; i++) {
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -974 +1160 @@
                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
-                psVectorAppend(norms, stamp->norm);
+
+                normSquare1 += stamp->normSquare1;
+                normSquare2 += stamp->normSquare2;
+
                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
                 numStamps++;
@@ -983 +1172 @@
 
 #if 0
-        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
-        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
+        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
+        psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
+        psVectorWriteFile("sumVector.dat", sumVector);
+        psFree (save);
 #endif
+
+        if (PENALTY) {
+            calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2);
+        }
 
         psVector *solution = NULL;                       // Solution to equation!
@@ -991 +1186 @@
         psVectorInit(solution, 0);
 
-#if 0
-        // Regular, straight-forward solution
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-#else
-        {
-            // Solve normalisation and background separately
-            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
-
-            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
-            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
-                psError(PM_ERR_DATA, false, "Unable to determine median normalisation");
-                psFree(stats);
-                psFree(sumMatrix);
-                psFree(sumVector);
-                psFree(norms);
-                return false;
-            }
-
-            // double normValue = 1.0;
-            double normValue = stats->robustMedian;
-            // double bgValue = 0.0;
-
-            psFree(stats);
-
-#ifdef TESTING
-            fprintf(stderr, "Norm: %lf\n", normValue);
-#endif
-            // Solve kernel components
-            for (int i = 0; i < numSolution1; i++) {
-                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
-
-                sumMatrix->data.F64[i][normIndex] = 0.0;
-                sumMatrix->data.F64[normIndex][i] = 0.0;
-            }
-            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
-            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
-            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
-
-            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
-            sumVector->data.F64[normIndex] = 0.0;
-
-            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-
-            solution->data.F64[normIndex] = normValue;
-        }
-# endif
-
-#if (1)
-        for (int i = 0; i < solution->n; i++) {
-            fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]);
-        }
-#endif
-
-        if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64);
-            psVectorInit(kernels->solution1, 0.0);
-        }
-
-        // only update the solutions that we chose to calculate:
-        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
-            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
-        }
-        if (mode & PM_SUBTRACTION_EQUATION_BG) {
-            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
-            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
-        }
-        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-            int numKernels = kernels->num;
-            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
-            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
-            for (int i = 0; i < numKernels * numPoly; i++) {
-                kernels->solution1->data.F64[i] = solution->data.F64[i];
-            }
-        }
-
-        psFree(norms);
-        psFree(solution);
-        psFree(sumVector);
-        psFree(sumMatrix);
-
-#ifdef TESTING
-        // XXX double-check for NAN in data:
-        for (int ix = 0; ix < kernels->solution1->n; ix++) {
-            if (!isfinite(kernels->solution1->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector\n");
-            }
-        }
-#endif
-
-    } else {
-        // Dual convolution solution
-
-        // Accumulation of stamp matrices/vectors
-        psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
-        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
-        psImageInit(sumMatrix, 0.0);
-        psVectorInit(sumVector, 0.0);
-
-        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
-
-        int numStamps = 0;              // Number of good stamps
-        for (int i = 0; i < stamps->num; i++) {
-            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
-                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
-                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
-
-                psVectorAppend(norms, stamp->norm);
-
-                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
-                numStamps++;
-            }
-        }
-
-#ifdef TESTING
-        psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
-        psVectorWriteFile("sumVector.dat", sumVector);
-#endif
-
-#if 1
-        // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
-        // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
-
-        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0);
-#endif
-
-        psVector *solution = NULL;                       // Solution to equation!
-        solution = psVectorAlloc(numParams, PS_TYPE_F64);
-        psVectorInit(solution, 0);
-
-#if 0
-        // Regular, straight-forward solution
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-#else
-        {
-            // Solve normalisation and background separately
-            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
-
-#if 0
-            psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64);
-            psVector *normVector = psVectorAlloc(2, PS_TYPE_F64);
-
-            normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex];
-            normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex];
-            normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex];
-
-            normVector->data.F64[0] = sumVector->data.F64[normIndex];
-            normVector->data.F64[1] = sumVector->data.F64[bgIndex];
-
-            psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN);
-
-            double normValue = normSolution->data.F64[0];
-            double bgValue = normSolution->data.F64[1];
-
-            psFree(normMatrix);
-            psFree(normVector);
-            psFree(normSolution);
-#endif
-
-            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
-            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
-                psError(PM_ERR_DATA, false, "Unable to determine median normalisation");
-                psFree(stats);
-                psFree(sumMatrix);
-                psFree(sumVector);
-                psFree(norms);
-                return false;
-            }
-
-            double normValue = stats->robustMedian;
-
-            psFree(stats);
-
-#ifdef TESTING
-            fprintf(stderr, "Norm: %lf\n", normValue);
-#endif
-
-            // Solve kernel components
-            for (int i = 0; i < numSolution2; i++) {
-                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
-                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1];
-
-                sumMatrix->data.F64[i][normIndex] = 0.0;
-                sumMatrix->data.F64[normIndex][i] = 0.0;
-
-                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
-                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
-            }
-            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
-            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
-            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
-
-            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
-
-            sumVector->data.F64[normIndex] = 0.0;
-
-            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-
-            solution->data.F64[normIndex] = normValue;
-        }
-#endif
-
-
-#if (1)
+        // DUAL solution
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+
+#if (0)
         for (int i = 0; i < solution->n; i++) {
             fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
@@ -1204 +1195 @@
 #endif
 
-        psFree(sumMatrix);
-        psFree(sumVector);
-
-        psFree(norms);
-
         if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64);
+            kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64);
             psVectorInit (kernels->solution1, 0.0);
         }
@@ -1218 +1204 @@
         }
 
-        // only update the solutions that we chose to calculate:
-        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
-            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
-        }
-        if (mode & PM_SUBTRACTION_EQUATION_BG) {
-            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
-            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
-        }
-        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-            int numKernels = kernels->num;
-            for (int i = 0; i < numKernels * numSpatial; i++) {
-                // XXX fprintf (stderr, "keep\n");
-                kernels->solution1->data.F64[i] = solution->data.F64[i];
-                kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
-            }
-        }
-
-
-        memcpy(kernels->solution1->data.F64, solution->data.F64,
-               numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
-        memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1],
-               numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
-
+        // for DUAL convolution analysis, we apply the normalization to I1 as follows:
+        // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i
+        // I2c =      I2 + \Sum_i b_i      I2 \cross k_i
+
+        // We absorb the normalization into a_i after the analysis is complete to be consistent
+        // with the SINGLE definitions of the convolutions
+
+        int numKernels = kernels->num;
+        for (int i = 0; i < numKernels * numSpatial; i++) {
+            // we solve for coefficients 
+            kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue;
+            kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
+        }
+
+        // Apply the normalisation and background separately
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
+        kernels->solution1->data.F64[normIndex] = stamps->normValue;
+        kernels->solution1->data.F64[bgIndex] = 0.0;
+
+        psFree(sumMatrix);
+        psFree(sumVector);
         psFree(solution);
-
     }
 
@@ -1265 +1248 @@
 }
 
-bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
-
-    // XXX measure some useful stats on the residuals
+// measure some useful stats on the stamp residuals:
+// fResSigma : the residual stdev / total flux
+// fResOuter : the residual fabs / total flux for R > 2 pix
+// fResTotal : the residual fabs / total flux for R > 0 pix
+bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
+
     float sum = 0.0;
     float peak = 0.0;
@@ -1277 +1263 @@
     }
 
-    // only count pixels with more than X% of the source flux
-    // calculate stdev(dflux)
+    // init counters
+    int npix = 0;
     float dflux1 = 0.0;
     float dflux2 = 0.0;
-    int npix = 0;
-
-    float dmax = 0.0;
-    float dmin = 0.0;
+    float dOuter = 0.0;
+    float dTotal = 0.0;
 
     for (int y = - footprint; y <= footprint; y++) {
         for (int x = - footprint; x <= footprint; x++) {
-            float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
-            if (dflux < 0.02*sum) continue;
             dflux1 += residual->kernel[y][x];
             dflux2 += PS_SQR(residual->kernel[y][x]);
-            dmax = PS_MAX(residual->kernel[y][x], dmax);
-            dmin = PS_MIN(residual->kernel[y][x], dmin);
+            dTotal += fabs(residual->kernel[y][x]);
+            if (hypot(x,y) > 2.0) {
+              dOuter += fabs(residual->kernel[y][x]);
+            }
             npix ++;
         }
@@ -1299 +1283 @@
     float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
     if (!isfinite(sum))  return false;
-    if (!isfinite(dmax)) return false;
-    if (!isfinite(dmin)) return false;
     if (!isfinite(peak)) return false;
-
-    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
-    psVectorAppend(fSigRes, sigma/sum);
-    psVectorAppend(fMaxRes, dmax/peak);
-    psVectorAppend(fMinRes, dmin/peak);
+    if (!isfinite(dOuter)) return false;
+    if (!isfinite(dTotal)) return false;
+
+    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum);
+    psVectorAppend(fResSigma, sigma/sum);
+    psVectorAppend(fResOuter, dOuter/sum);
+    psVectorAppend(fResTotal, dTotal/sum);
     return true;
 }
@@ -1329 +1313 @@
     pmSubtractionVisualShowFitInit (stamps);
 
-    psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
-    psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
-    psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
 
     // we want to save the residual images for the 9 brightest stamps.
@@ -1443 +1427 @@
             }
 
-            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
+            pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint);
 
         } else {
@@ -1479 +1463 @@
             }
 
-            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
+            pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint);
         }
 
@@ -1558 +1542 @@
 
         psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-        psVectorStats (stats, fSigRes, NULL, NULL, 0);
-        kernels->fSigResMean = stats->robustMedian;
-        kernels->fSigResStdev = stats->robustStdev;
+        psVectorStats (stats, fResSigma, NULL, NULL, 0);
+        kernels->fResSigmaMean  = stats->robustMedian;
+        kernels->fResSigmaStdev = stats->robustStdev;
 
         psStatsInit (stats);
-        psVectorStats (stats, fMaxRes, NULL, NULL, 0);
-        kernels->fMaxResMean = stats->robustMedian;
-        kernels->fMaxResStdev = stats->robustStdev;
+        psVectorStats (stats, fResOuter, NULL, NULL, 0);
+        kernels->fResOuterMean  = stats->robustMedian;
+        kernels->fResOuterStdev = stats->robustStdev;
 
         psStatsInit (stats);
-        psVectorStats (stats, fMinRes, NULL, NULL, 0);
-        kernels->fMinResMean = stats->robustMedian;
-        kernels->fMinResStdev = stats->robustStdev;
+        psVectorStats (stats, fResTotal, NULL, NULL, 0);
+        kernels->fResTotalMean  = stats->robustMedian;
+        kernels->fResTotalStdev = stats->robustStdev;
 
         // XXX save these values somewhere
-        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",
-                 kernels->fSigResMean, kernels->fSigResStdev,
-                 kernels->fMaxResMean, kernels->fMaxResStdev,
-                 kernels->fMinResMean, kernels->fMinResStdev);
-
-        psFree (fSigRes);
-        psFree (fMaxRes);
-        psFree (fMinRes);
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f",
+                 kernels->fResSigmaMean, kernels->fResSigmaStdev,
+                 kernels->fResOuterMean, kernels->fResOuterStdev,
+                 kernels->fResTotalMean, kernels->fResTotalStdev);
+
+        psFree (fResSigma);
+        psFree (fResOuter);
+        psFree (fResTotal);
         psFree (stats);
     }
@@ -1586 +1570 @@
     psFree(residual);
     psFree(polyValues);
-
 
     return deviations;
@@ -1912 +1895 @@
     return true;
 }
-
-
-# if 0
-
-#ifdef TESTING
-        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
-        psVectorWriteFile ("B.dat", sumVector);
-#endif
-
-# define SVD_ANALYSIS 0
-# define COEFF_SIG 0.0
-# define SVD_TOL 0.0
-
-        // Use SVD to determine the kernel coeffs (and validate)
-        if (SVD_ANALYSIS) {
-
-            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
-            // sumMatrix * x = sumVector.
-
-            // we can use any standard matrix inversion to solve this.  However, the basis
-            // functions in general have substantial correlation, so that the solution may be
-            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
-            // system of equations may be statistically ill-conditioned.  Noise in the image
-            // will drive insignificant, but correlated, terms in the solution.  To avoid these
-            // problems, we can use SVD to identify numerically unconstrained values and to
-            // avoid statistically badly determined value.
-
-            // A = sumMatrix, B = sumVector
-            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
-            // x = V (1/w) (U^T B)
-            // \sigma_x = sqrt(diag(A^{-1}))
-            // solve for x and A^{-1} to get x & dx
-            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
-            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
-
-            // If I use the SVD trick to re-condition the matrix, I need to break out the
-            // kernel and normalization terms from the background term.
-            // XXX is this true?  or was this due to an error in the analysis?
-
-            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
-
-            // now pull out the kernel elements into their own square matrix
-            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
-            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
-
-            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
-                if (ix == bgIndex) continue;
-                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
-                    if (iy == bgIndex) continue;
-                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
-                    ky++;
-                }
-                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
-                kx++;
-            }
-
-            psImage *U = NULL;
-            psImage *V = NULL;
-            psVector *w = NULL;
-            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
-                psError(psErrorCodeLast(), false, "failed to perform SVD on sumMatrix\n");
-                return NULL;
-            }
-
-            // calculate A_inverse:
-            // Ainv = V * w * U^T
-            psImage *wUt  = p_pmSubSolve_wUt (w, U);
-            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
-            psImage *Xvar = NULL;
-            psFree (wUt);
-
-# ifdef TESTING
-            // kernel terms:
-            for (int i = 0; i < w->n; i++) {
-                fprintf (stderr, "w: %f\n", w->data.F64[i]);
-            }
-# endif
-            // loop over w adding in more and more of the values until chisquare is no longer
-            // dropping significantly.
-            // XXX this does not seem to work very well: we seem to need all terms even for
-            // simple cases...
-
-            psVector *Xsvd = NULL;
-            {
-                psVector *Ax = NULL;
-                psVector *UtB = NULL;
-                psVector *wUtB = NULL;
-
-                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
-                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
-                psVectorInit (wMask, 1); // start by masking everything
-
-                double chiSquareLast = NAN;
-                int maxWeight = 0;
-
-                double Axx, Bx, y2;
-
-                // XXX this is an attempt to exclude insignificant modes.
-                // it was not successful with the ISIS kernel set: removing even
-                // the least significant mode leaves additional ringing / noise
-                // because the terms are so coupled.
-                for (int k = 0; false && (k < w->n); k++) {
-
-                    // unmask the k-th weight
-                    wMask->data.U8[k] = 0;
-                    p_pmSubSolve_SetWeights(wApply, w, wMask);
-
-                    // solve for x:
-                    // x = V * w * (U^T * B)
-                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
-                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
-                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
-
-                    // chi-square for this system of equations:
-                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
-                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
-                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
-                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
-                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
-                    p_pmSubSolve_y2 (&y2, kernels, stamps);
-
-                    // apparently, this works (compare with the brute force value below
-                    double chiSquare = Axx - 2.0*Bx + y2;
-                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
-                    chiSquareLast = chiSquare;
-
-                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
-                    if (k && !maxWeight && (deltaChi < 1.0)) {
-                        maxWeight = k;
-                    }
-                }
-
-                // keep all terms or we get extra ringing
-                maxWeight = w->n;
-                psVectorInit (wMask, 1);
-                for (int k = 0; k < maxWeight; k++) {
-                    wMask->data.U8[k] = 0;
-                }
-                p_pmSubSolve_SetWeights(wApply, w, wMask);
-
-                // solve for x:
-                // x = V * w * (U^T * B)
-                p_pmSubSolve_UtB (&UtB, U, kernelVector);
-                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
-                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
-
-                // chi-square for this system of equations:
-                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
-                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
-                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
-                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
-                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
-                p_pmSubSolve_y2 (&y2, kernels, stamps);
-
-                // apparently, this works (compare with the brute force value below
-                double chiSquare = Axx - 2.0*Bx + y2;
-                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
-
-                // re-calculate A^{-1} to get new variances:
-                // Ainv = V * w * U^T
-                // XXX since we keep all terms, this is identical to Ainv
-                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
-                Xvar = p_pmSubSolve_VwUt (V, wUt);
-                psFree (wUt);
-
-                psFree (Ax);
-                psFree (UtB);
-                psFree (wUtB);
-                psFree (wApply);
-                psFree (wMask);
-            }
-
-            // copy the kernel solutions to the full solution vector:
-            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
-            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
-
-            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
-                if (ix == bgIndex) {
-                    solution->data.F64[ix] = 0;
-                    solutionErr->data.F64[ix] = 0.001;
-                    continue;
-                }
-                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
-                solution->data.F64[ix] = Xsvd->data.F64[kx];
-                kx++;
-            }
-
-            psFree (kernelMatrix);
-            psFree (kernelVector);
-
-            psFree (U);
-            psFree (V);
-            psFree (w);
-
-            psFree (Ainv);
-            psFree (Xsvd);
-        } else {
-            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
-            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
-            if (!luMatrix) {
-                psError(PM_ERR_DATA, true, "LU Decomposition of least-squares matrix failed.\n");
-                psFree(solution);
-                psFree(sumVector);
-                psFree(sumMatrix);
-                psFree(luMatrix);
-                psFree(permutation);
-                return NULL;
-            }
-
-            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
-            psFree(luMatrix);
-            psFree(permutation);
-            if (!solution) {
-                psError(PM_ERR_DATA, true, "Failed to solve the least-squares system.\n");
-                psFree(solution);
-                psFree(sumVector);
-                psFree(sumMatrix);
-                return NULL;
-            }
-
-            // XXX LUD does not provide A^{-1}?  fake the error for now
-            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
-            for (int ix = 0; ix < sumVector->n; ix++) {
-                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
-            }
-        }
-
-        if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
-            psVectorInit (kernels->solution1, 0.0);
-        }
-
-        // only update the solutions that we chose to calculate:
-        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
-            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
-        }
-        if (mode & PM_SUBTRACTION_EQUATION_BG) {
-            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
-            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
-        }
-        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
-            int numKernels = kernels->num;
-            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
-            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
-            for (int i = 0; i < numKernels * numPoly; i++) {
-                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
-                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
-                    // XXX fprintf (stderr, "drop\n");
-                    kernels->solution1->data.F64[i] = 0.0;
-                } else {
-                    // XXX fprintf (stderr, "keep\n");
-                    kernels->solution1->data.F64[i] = solution->data.F64[i];
-                }
-            }
-        }
-        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
-        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
-
-        psFree(solution);
-        psFree(sumVector);
-        psFree(sumMatrix);
-# endif
-
-#ifdef TESTING
-              // XXX double-check for NAN in data:
-                for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
-                    for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
-                        if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
-                            fprintf (stderr, "WARNING: NAN in matrix\n");
-                        }
-                    }
-                }
-                for (int ix = 0; ix < stamp->vector->n; ix++) {
-                    if (!isfinite(stamp->vector->data.F64[ix])) {
-                        fprintf (stderr, "WARNING: NAN in vector\n");
-                    }
-                }
-#endif
-
-#ifdef TESTING
-        for (int ix = 0; ix < sumVector->n; ix++) {
-            if (!isfinite(sumVector->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector\n");
-            }
-        }
-#endif
-
-#ifdef TESTING
-        for (int ix = 0; ix < sumVector->n; ix++) {
-            if (!isfinite(sumVector->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector\n");
-            }
-        }
-        {
-            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
-            psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
-            psFree(inverse);
-        }
-        {
-            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
-            psImage *Xt = psMatrixTranspose(NULL, X);
-            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
-            psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
-            psFree(X);
-            psFree(Xt);
-            psFree(XtX);
-        }
-#endif
-

trunk/psModules/src/imcombine/pmSubtractionEquation.h

@@ -65 +65 @@
     );
 
+bool pmSubtractionCalculateNormalization(
+  pmSubtractionStampList *stamps,
+  const pmSubtractionMode mode);
+
+bool pmSubtractionCalculateNormalizationStamp(
+    pmSubtractionStamp *stamp,          // stamp on which to save normalization)
+    const psKernel *input,              // Input image (target)
+    const psKernel *reference,          // Reference image (convolution source)
+    int footprint,                      // (Half-)Size of stamp
+    int normWindow1,                    // Window (half-)size for normalisation measurement
+    int normWindow2                     // Window (half-)size for normalisation measurement
+  );
+
+bool pmSubtractionCalculateMoments(
+    pmSubtractionKernels *kernels, // Kernels
+    pmSubtractionStampList *stamps);
+
+bool pmSubtractionCalculateMomentsStamp(
+    pmSubtractionKernels *kernels, // Kernels
+    pmSubtractionStamp *stamp,          // stamp on which to save normalization)
+    int footprint,                      // (Half-)Size of stamp
+    int normWindow1,                    // Window (half-)size for normalisation measurement
+    int normWindow2                     // Window (half-)size for normalisation measurement
+    );
+
+bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window);
+
 #endif

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -183 +183 @@
 }
 
+# define CENTRAL_DELTA 0
+
+# if (CENTRAL_DELTA)
+
+// XXX *** this code used the central pixel to force zero net flux,
+// Alard actually uses kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0])
 static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
                                                 int index, int uOrder, int vOrder, float fwhm,
@@ -220 +226 @@
             // Re-normalize 
             // scale2D  = 1.0 / fabs(sum);
-            scale2D  = 1.0 / sqrt(sum2);
+            scale2D  = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
             zeroNull = true;
         } else {
             // Odd functions: choose normalisation so that parameters have about the same strength as for even
             // functions, no subtraction of null pixel because the sum is already (near) zero
-            scale2D = 1.0 / sqrt(sum2);
+            scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
             zeroNull = false;
         }
@@ -235 +241 @@
     if (forceZeroNull) {
         // Force rescaling and subtraction of null pixel even though the order doesn't indicate it's even
-        scale2D = 1.0 / fabs(sum);
+        scale2D = 1.0 / fabs(sum) / PS_SQR(fwhm);
         zeroNull = true;
     }
     if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
         // Odd function
-        scale2D = 1.0 / sqrt(sum2);
+        scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
     }
 
@@ -255 +261 @@
     if (zeroNull) {
         // preCalc->kernel->kernel[0][0] -= 1.0;
-        preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
-    }
-
-#if 1
+        preCalc->kernel->kernel[0][0] -= sum * scale2D;
+    }
+
+#if 0
     {
         double Sum = 0.0;   // Sum of kernel component
@@ -287 +293 @@
     return true;
 }
+
+# else /* CENTRAL_DELTA */
+
+static bool zeroIsNormal = false;
+
+// this code uses kernel(0) to force zero flux, and is invalid for other kinds of normalizations
+static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
+                                                int index, int uOrder, int vOrder, float fwhm,
+                                                bool AlardLuptonStyle, bool forceZeroNull)
+{
+    // 1) for odd functions: no renormalization 
+    // 2) for even functions, normalize the kernel to unity
+    // 3) for even functions & index > 0, subtract kernel(0)
+
+    // Calculate moments
+    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
+    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+
+    for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+        for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+            double value = preCalc->kernel->kernel[v][u];
+            double value2 = PS_SQR(value);
+            sum += value;
+            sum2 += value2;
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+        }
+    }
+
+#if 0
+    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max);
+#endif
+
+    float scale2D = 1.0;                // Scaling for 2-D kernels
+    float scale1D = 1.0;                // Scaling for 1-D kernels
+
+    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+
+        scale2D = 1.0 / sum;            // Scaling for 2-D kernels
+        scale1D = sqrtf(scale2D);               // Scaling for 1-D kernels
+
+        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+                preCalc->kernel->kernel[v][u] *= scale2D;
+            }
+        }
+        if (index == 0) {
+            zeroIsNormal = true;
+        } else {
+            psAssert(zeroIsNormal, "failed to normalize zero kernel first");
+            pmSubtractionKernelPreCalc *zeroKernel = kernels->preCalc->data[0];
+            psAssert(zeroKernel, "failed to supply zero kernel");
+            for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+                for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+                    preCalc->kernel->kernel[v][u] -= zeroKernel->kernel->kernel[v][u];
+                }
+            }
+        }
+
+        // XXX why do we bother renormlizing the 1D kernels?  I don't think we use them again...
+        if (preCalc->xKernel) {
+            psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+        }
+        if (preCalc->yKernel) {
+            psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+        }
+    }
+
+#if 0
+    {
+        double Sum = 0.0;   // Sum of kernel component
+        double Sum2 = 0.0;   // Sum of kernel component
+        float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+                double value = preCalc->kernel->kernel[v][u];
+                Sum += value;
+                Sum2 += PS_SQR(value);
+                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
+                max = PS_MAX(preCalc->kernel->kernel[v][u], max);
+            }
+        }
+        fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D);
+    }
+#endif
+
+    kernels->widths->data.F32[index] = fwhm;
+    kernels->u->data.S32[index] = uOrder;
+    kernels->v->data.S32[index] = vOrder;
+    if (kernels->preCalc->data[index]) {
+        psFree(kernels->preCalc->data[index]);
+    }
+    kernels->preCalc->data[index] = preCalc;
+    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
+
+    return true;
+}
+# endif /* Central Delta */
 
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
@@ -603 +707 @@
     kernels->sampleStamps = NULL;
 
-    kernels->fSigResMean  = NAN;
-    kernels->fSigResStdev = NAN;
-    kernels->fMaxResMean  = NAN;
-    kernels->fMaxResStdev = NAN;
-    kernels->fMinResMean  = NAN;
-    kernels->fMinResStdev = NAN;
+    kernels->fResSigmaMean  = NAN;
+    kernels->fResSigmaStdev = NAN;
+    kernels->fResOuterMean  = NAN;
+    kernels->fResOuterStdev = NAN;
+    kernels->fResTotalMean  = NAN;
+    kernels->fResTotalStdev = NAN;
 
     return kernels;

trunk/psModules/src/imcombine/pmSubtractionKernels.h

@@ -51 +51 @@
     float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
     int numStamps;                      ///< Number of good stamps
-    float fSigResMean;                  ///< mean fractional stdev of residuals
-    float fSigResStdev;                 ///< stdev of fractional stdev of residuals
-    float fMaxResMean;                  ///< mean fractional positive swing in residuals
-    float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
-    float fMinResMean;                  ///< mean fractional negative swing in residuals
-    float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
+    float fResSigmaMean;                ///< mean fractional stdev of residuals
+    float fResSigmaStdev;               ///< stdev of fractional stdev of residuals
+    float fResOuterMean;                ///< mean fractional positive swing in residuals
+    float fResOuterStdev;               ///< stdev of fractional positive swing in residuals
+    float fResTotalMean;                ///< mean fractional negative swing in residuals
+    float fResTotalStdev;               ///< stdev of fractional negative swing in residuals
     psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
 } pmSubtractionKernels;

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -28 +28 @@
 static bool useFFT = true;              // Do convolutions using FFT
 
-# define SEPARATE 0
-# if (SEPARATE)
-# define SUBMODE PM_SUBTRACTION_EQUATION_NORM
-# else
 # define SUBMODE PM_SUBTRACTION_EQUATION_ALL
-# endif
 
 //#define TESTING
@@ -280 +275 @@
 }
 
-bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {
-
-    if (!readout) return true;
-
-    psImage *image = readout->image;
-    psImage *mask  = readout->mask;
-    psImage *variance = readout->variance;
-    for (int y = 0; y < image->numRows; y++) {
-        for (int x = 0; x < image->numCols; x++) {
-            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;
-            bool valid = false;
-            valid = isfinite(image->data.F32[y][x]);
-            if (variance) {
-                valid &= isfinite(variance->data.F32[y][x]);
-            }
-            if (valid) continue;
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;
-        }
-    }
-
-    return true;
-}
+// bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {
+// 
+//     if (!readout) return true;
+// 
+//     psImage *image = readout->image;
+//     psImage *mask  = readout->mask;
+//     psImage *variance = readout->variance;
+//     for (int y = 0; y < image->numRows; y++) {
+//         for (int x = 0; x < image->numCols; x++) {
+//             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;
+//             bool valid = false;
+//             valid = isfinite(image->data.F32[y][x]);
+//             if (variance) {
+//                 valid &= isfinite(variance->data.F32[y][x]);
+//             }
+//             if (valid) continue;
+//             mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;
+//         }
+//     }
+// 
+//     return true;
+// }
 
 bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
@@ -392 +387 @@
     }
 
-    pmSubtractionMaskInvalid(ro1, maskVal);
-    pmSubtractionMaskInvalid(ro2, maskVal);
+    // XXX this is done before calling this function
+    // pmSubtractionMaskInvalid(ro1, maskVal);
+    // pmSubtractionMaskInvalid(ro2, maskVal);
 
     // General background subtraction, since this is done in pmSubtractionMatch
@@ -565 +561 @@
     memCheck("start");
 
-    pmSubtractionMaskInvalid(ro1, maskVal);
-    pmSubtractionMaskInvalid(ro2, maskVal);
+    // pmSubtractionMaskInvalid(ro1, maskVal);
+    // pmSubtractionMaskInvalid(ro2, maskVal);
 
     psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
@@ -672 +668 @@
                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
                                                        inner, binning, ringsOrder, penalty, bounds, subMode);
-                // pmSubtractionVisualShowKernels(kernels);
             }
 
@@ -678 +673 @@
 
             if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
-#if 0
-                // Get backgrounds
-                psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background
-                psVector *buffer = NULL;// Buffer for stats
-                if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) {
-                    psError(PM_ERR_DATA, false, "Unable to measure background of image 1.");
-                    psFree(bgStats);
-                    psFree(buffer);
-                    goto MATCH_ERROR;
-                }
-                float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1
-                if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) {
-                    psError(PM_ERR_DATA, false, "Unable to measure background of image 2.");
-                    psFree(bgStats);
-                    psFree(buffer);
-                    goto MATCH_ERROR;
-                }
-                float bg2 = psStatsGetValue(bgStats, BG_STAT); // Background for image 2
-                psFree(bgStats);
-                psFree(buffer);
-
-                pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use
-#endif
                 pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
                 switch (newMode) {
@@ -732 +704 @@
                 }
 
-                // XXX step 1: calculate normalization
-                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue
+                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
+                if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+                // step 0b : calculate the moments, pass along to the next steps via stamps->normValue
+                // XXX this step is now skipped -- delete
+                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
+                if (!pmSubtractionCalculateMoments(kernels, stamps)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+                // step 1: generate the elements of the matrix equation Ax = B
+                psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
@@ -739 +726 @@
                 }
 
-                psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n");
+                // step 2: solve the matrix equation Ax = B
+                psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
@@ -746 +734 @@
                 memCheck("  solve equation");
 
-# if (SEPARATE)
-                // set USED -> CALCULATE
-                pmSubtractionStampsResetStatus (stamps);
-
-                // XXX step 2: calculate kernel parameters
-                psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-                memCheck("  calculate equation");
-
-                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
-                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-                memCheck("  solve equation");
-# endif
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
@@ -770 +739 @@
                     goto MATCH_ERROR;
                 }
-
                 memCheck("   calculate deviations");
 
@@ -787 +755 @@
             // if we hit the max number of iterations and we have rejected stamps, re-solve
             if (numRejected > 0) {
-                // XXX step 1: calculate normalization
+
+                // step 1: generate the elements of the matrix equation Ax = B
                 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
@@ -794 +763 @@
                 }
 
-                // solve normalization
+                // step 2: solve the matrix equation Ax = B
                 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
@@ -801 +770 @@
                 }
 
-# if (SEPARATE)
-                // set USED -> CALCULATE
-                pmSubtractionStampsResetStatus (stamps);
-
-                // XXX step 2: calculate kernel parameters
-                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-
-                // solve kernel parameters
-                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
-                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-                memCheck("  solve equation");
-# endif
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
@@ -837 +787 @@
                 goto MATCH_ERROR;
             }
-
             memCheck("diag outputs");
 
@@ -1134 +1083 @@
     }
 
-# if (SEPARATE)
-    // set USED -> CALCULATE
-    pmSubtractionStampsResetStatus (stamps);
-
-    psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description);
-    if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
-        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-        return false;
-    }
-
-    psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description);
-    if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
-        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-        return false;
-    }
-# endif
-
     psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
     psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
@@ -1181 +1113 @@
         }
         psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
-
-# if (SEPARATE)
-        // set USED -> CALCULATE
-        pmSubtractionStampsResetStatus (stamps);
-
-        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
-        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
-            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-            return false;
-        }
-
-        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
-        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
-            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-            return false;
-        }
-        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
-# endif
 
         psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
@@ -1288 +1202 @@
 
 bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
-                              float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax)
+                              float scaleRef, float scaleMin, float scaleMax)
 {
     PS_ASSERT_PTR_NON_NULL(kernelSize, false);
@@ -1294 +1208 @@
     PS_ASSERT_VECTOR_NON_NULL(widths, false);
     PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);
     PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
     PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
@@ -1301 +1213 @@
     PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
 
-//    float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
+    float fwhm1;
+    float fwhm2;
+
+    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+    psAssert(isfinite(fwhm1), "fwhm 1 not set");
+    psAssert(isfinite(fwhm2), "fwhm 2 not set");
+    
+    // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
     float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
-
-    // XXX save these values in a static for later use
-    pmSubtractionSetFWHMs(fwhm1, fwhm2);
 
     if (isfinite(scaleMin) && scale < scaleMin) {

trunk/psModules/src/imcombine/pmSubtractionMatch.h

@@ -104 +104 @@
     int *stampSize,                     ///< Half-size of the stamp (footprint)
     psVector *widths,                   ///< ISIS widths
-    float fwhm1, float fwhm2,           ///< FWHMs for inputs
     float scaleRef,                     ///< Reference width for scaling
     float scaleMin,                     ///< Minimum scaling ratio, or NAN

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -51 +51 @@
     psFree(list->y);
     psFree(list->flux);
+    psFree(list->window);
     psFree(list->window1);
     psFree(list->window2);
@@ -96 +97 @@
 
     // fprintf (stderr, "xMin, xMax: %d, %d -> ", xMin, xMax);
+    // float xRaw = xMin;
+    // float yRaw = yMin;
 
     // Ensure we're not going to go outside the bounds of the image
@@ -103 +106 @@
     yMax = PS_MIN(numRows - border - 1, yMax);
 
-    if (xMax < xMin) return false;
-    if (yMax < yMin) return false;
+    if (xMax < xMin) {
+        // fprintf (stderr, "%f,%f : x-border\n", xRaw, yRaw);
+        return false;
+    }
+    if (yMax < yMin) {
+        // fprintf (stderr, "%f,%f : y-border\n", xRaw, yRaw);
+        return false;
+    }
 
     psAssert (xMin <= xMax, "x mismatch?");
@@ -118 +127 @@
             if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] &
                 (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) {
+                // fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw);
                 return false;
             }
@@ -133 +143 @@
     }
 
+    // if (!found) {
+    //  fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw);
+    // }
+
     return found;
 }
@@ -232 +246 @@
     list->flux = NULL;
     list->normFrac = normFrac;
+    list->normValue = NAN;
+    list->window = NULL;
     list->window1 = NULL;
     list->window2 = NULL;
@@ -256 +272 @@
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->window1 = psMemIncrRefCounter(in->window1);
     out->window2 = psMemIncrRefCounter(in->window2);
@@ -331 +348 @@
     stamp->vector = NULL;
     stamp->norm = NAN;
+    stamp->normI1 = NAN;
+    stamp->normI2 = NAN;
+    stamp->normSquare1 = NAN;
+    stamp->normSquare2 = NAN;
+
+    stamp->MxxI1 = NULL;
+    stamp->MyyI1 = NULL;
+    stamp->MxxI2 = NULL;
+    stamp->MyyI2 = NULL;
+
+    stamp->MxxI1raw = NAN;
+    stamp->MyyI1raw = NAN;
+    stamp->MxxI2raw = NAN;
+    stamp->MyyI2raw = NAN;
 
     return stamp;
@@ -431 +462 @@
                 // Take stamps off the top of the (sorted) list
                 for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
+                    // fprintf (stderr, "%d : xList: %ld elements\n", i, xList->n);
                     // Chop off the top of the list
                     xList->n = j;
@@ -459 +491 @@
                                             yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
 #endif
-                    if (0 && goodStamp) {
-                        fprintf (stderr, "find: %6.1f < %6.1f < %6.1f, %6.1f < %6.1f %6.1f\n",
-                                 region->x0 + size, xStamp, region->x1 - size,
-                                 region->y0 + size, yStamp, region->y1 - size);
-                    }
                 }
             } else {
@@ -656 +683 @@
     psImageInit(stamps->window2->image, 0.0);
 
+    // Generate a weighting window based on the fwhms (20% larger than the largest)
+    { 
+        float fwhm1, fwhm2;
+
+        // XXX this is annoyingly hack-ish
+        pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+
+        float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
+
+        psFree (stamps->window);
+        stamps->window = psKernelAlloc(-size, size, -size, size);
+        psImageInit(stamps->window->image, 0.0);
+
+        for (int y = -size; y <= size; y++) {
+            for (int x = -size; x <= size; x++) {
+                stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
+            }
+        }
+    }
+
     // generate normalizations for each stamp
     psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32);
@@ -678 +725 @@
 
     // storage vector for flux data
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
     psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
@@ -706 +753 @@
                 psAbort ("failed to generate stats");
             }
-            float f1 = stats->robustMedian;
+            float f1 = stats->sampleMedian;
 
             psStatsInit (stats);
@@ -712 +759 @@
                 psAbort ("failed to generate stats");
             }
-            float f2 = stats->robustMedian;
+            float f2 = stats->sampleMedian;
 
             stamps->window1->kernel[y][x] = f1;
@@ -728 +775 @@
     }
 
+#if 0
+    {
+        psFits *fits = NULL;
+        fits = psFitsOpen ("window1.raw.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
+        psFitsClose (fits);
+        fits = psFitsOpen ("window2.raw.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
+        psFitsClose (fits);
+    }
+#endif
+
     psTrace("psModules.imcombine", 3, "Window total (1): %f, threshold: %f\n", sum1, (1.0 - stamps->normFrac) * sum1);
     psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2);
 
-# if (0)
+# if (1)
     // this block attempts to calculate the radius based on the first radial moment
-    bool done1 = false;
-    bool done2 = false;
-    double prior1 = 0.0;
-    double prior2 = 0.0;
+    double Sr1 = 0.0;
+    double Sr2 = 0.0;
+    double Sf1 = 0.0;
+    double Sf2 = 0.0;
     for (int y = -size; y <= size; y++) {
         for (int x = -size; x <= size; x++) {
@@ -750 +809 @@
     float R2 = Sr2 / Sf2;
 
-    stamps->normWindow1 = 2.5*R1;
-    stamps->normWindow1 = 2.5*R2;
+    stamps->normWindow1 = 2.0*R1;
+    stamps->normWindow2 = 2.0*R2;
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
+
 # else
     // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent).
@@ -778 +839 @@
             // interpolate to the radius at which delta2 is normFrac:
             stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
-            fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
+            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
             done1 = true;
         }
@@ -785 +846 @@
             // interpolate to the radius at which delta2 is normFrac:
             stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
-            fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
+            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
             done2 = true;
         }
@@ -833 +894 @@
     {
         psFits *fits = NULL;
-        fits = psFitsOpen ("window1.fits", "w");
+        fits = psFitsOpen ("window1.norm.fits", "w");
         psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
         psFitsClose (fits);
-        fits = psFitsOpen ("window2.fits", "w");
+        fits = psFitsOpen ("window2.norm.fits", "w");
         psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
         psFitsClose (fits);
@@ -906 +967 @@
         float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull);
 
-        penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
-        penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        if (1) {
+            penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            // penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            // penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        } else {
+            penalty1 = PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+            penalty2 = PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        }
     }
     kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;
@@ -915 +983 @@
     psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
 
-    fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
+    // fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
 
     return true;
@@ -942 +1010 @@
     penalty *= 1.0 / sum2;
 
-    if (1) {
-        fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
+    if (0) {
+        // fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
         // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     }
@@ -983 +1051 @@
 
         if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) {
-            psError(PM_ERR_PROG, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
-            return false;
+          psLogMsg("psModules.imcombine", 3, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
+          stamp->status = PM_SUBTRACTION_STAMP_NONE;
+          continue;
         }
 

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -25 +25 @@
     int footprint;                      ///< Half-size of stamps
     float normFrac;                     ///< Fraction of flux in window for normalisation window
+    float normValue;                    ///< calculated normalization
+    float normValue2;                   ///< calculated normalization
     psKernel *window1;                  ///< window function generated from ensemble of stamps (input 1)
     psKernel *window2;                  ///< window function generated from ensemble of stamps (input 2)
-    float normWindow1;                    ///< Size of window for measuring normalisation
-    float normWindow2;                    ///< Size of window for measuring normalisation
+    psKernel *window;                   ///< weighting window function (sigma = 1.1 * MAX(fwhm))
+    float normWindow1;                  ///< Size of window for measuring normalisation
+    float normWindow2;                  ///< Size of window for measuring normalisation
     float sysErr;                       ///< Systematic error
     float skyErr;                       ///< increase effective readnoise
@@ -85 +88 @@
     psVector *vector;                   ///< Least-squares vector, or NULL
     double norm;                        ///< Normalisation difference
+    double normI1;                       ///< Sum(flux) for image 1
+    double normI2;                       ///< Sum(flux) for image 2
+    double normSquare1;                 ///< Sum(flux^2) for image 1 (used for penalty)
+    double normSquare2;                 ///< Sum(flux^2) for image 2 (used for penalty)
     pmSubtractionStampStatus status;    ///< Status of stamp
+    psVector *MxxI1;                    ///< second moments of convolution images
+    psVector *MyyI1;                    ///< second moments of convolution images
+    psVector *MxxI2;                    ///< second moments of convolution images
+    psVector *MyyI2;                    ///< second moments of convolution images
+    double MxxI1raw;
+    double MyyI1raw;
+    double MxxI2raw;
+    double MyyI2raw;
 } pmSubtractionStamp;
 

trunk/psModules/src/imcombine/pmSubtractionVisual.c

@@ -146 +146 @@
     }
     pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);
+    psFree(canvas);
 
     pmVisualAskUser(&plotStamps);
@@ -152 +153 @@
 
 /** Plot the least-squares matrix of each stamp */
-bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) {
+bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps) {
 
     if (!pmVisualTestLevel("ppsub.chisq", 1)) return true;
@@ -209 +210 @@
     pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
 
+    if (0) {
+        static int count = 0;
+        char filename[64];
+        sprintf (filename, "chisq.%02d.fits", count);
+        count ++;
+        psFits *fits = psFitsOpen (filename, "w");
+        psFitsWriteImage (fits, NULL, canvas32, 0, NULL);
+        psFitsClose (fits);
+    }
+
     pmVisualAskUser(&plotLeastSquares);
     psFree(canvas);
@@ -271 +282 @@
             }
         }
-        fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+        // fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     }                                                    
         
     pmVisualScaleImage(kapa1, output, "Image", 0, true);
+    psFree(output);
     pmVisualAskUser(&plotImage);
     return true;
@@ -286 +298 @@
         return false;
     }
+
+    // XXX clear the overlay(s) (red at least!)
 
     // get the kernel sizes
@@ -297 +311 @@
         if (!isfinite(stamp->flux)) continue;
         if (!stamp->convolutions1 && !stamp->convolutions2) continue;
+        // fprintf (stderr, "flux: %f, maxFlux: %f  ", stamp->flux, maxFlux);
         if (!maxStamp) {
             maxFlux = stamp->flux;
             maxStamp = stamp;
+            // fprintf (stderr, "maxStamp %d\n", i);
             continue;
+        } else {
+            // fprintf (stderr, "\n");
         }
         if (stamp->flux > maxFlux) {
@@ -337 +355 @@
             
             double sum = 0.0;
+            double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
                     sum += kernel->kernel[y][x];
+                    sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
         }               
         pmVisualScaleImage(kapa2, output, "Image", 0, true);
-    }                                    
+
+        if (0) {
+            psFits *fits = psFitsOpen("basis.1.fits", "w");
+            psFitsWriteImage(fits, NULL, output, 0, NULL);
+            psFitsClose(fits);
+        }
+    }
         
     if (maxStamp->convolutions2) {
@@ -371 +397 @@
             
             double sum = 0.0;
+            double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
                     sum += kernel->kernel[y][x];
+                    sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
         }               
         pmVisualScaleImage(kapa2, output, "Image", 1, true);
+
+        if (0) {
+            psFits *fits = psFitsOpen("basis.2.fits", "w");
+            psFitsWriteImage(fits, NULL, output, 0, NULL);
+            psFitsClose(fits);
+        }
+        psFree(output);
     }                                    
         
@@ -676 +711 @@
     psFree (x);
     psFree (y);
+    psFree (polyValues);
 
     pmVisualAskUser(NULL);