Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 17294)
+++ /trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 17295)
@@ -19,4 +19,5 @@
 
 #define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
+#define BG_STAT PS_STAT_ROBUST_MEDIAN   // Statistic to use for background
 
 static bool useFFT = true;              // Do convolutions using FFT
@@ -272,5 +273,29 @@
 
             if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) {
-                pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode
+                // Get backgrounds
+                psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background
+                psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+                psVector *buffer = NULL;// Buffer for stats
+                if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskBad, rng)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 1.");
+                    psFree(bgStats);
+                    psFree(rng);
+                    psFree(buffer);
+                    goto MATCH_ERROR;
+                }
+                float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1
+                if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskBad, rng)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 2.");
+                    psFree(bgStats);
+                    psFree(rng);
+                    psFree(buffer);
+                    goto MATCH_ERROR;
+                }
+                float bg2 = psStatsGetValue(bgStats, BG_STAT); // Background for image 2
+                psFree(bgStats);
+                psFree(rng);
+                psFree(buffer);
+
+                pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use
                 switch (newMode) {
                   case PM_SUBTRACTION_MODE_1:
@@ -529,98 +554,88 @@
 }
 
-// Calculate the second order moments for an image
-static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments
-                                    int radius   // Maximum radius
-                                    )
+
+// Determine a rough width (integer value) of the star in the image
+// XXX Could improve this by using a user-provided list of floating-point widths (or an end point and
+// increment).
+static int subtractionOrderWidth(const psKernel *kernel, // Image
+                                 float bg, // Background in image
+                                 int size, // Maximum size
+                                 psArray **models, // Buffer of models
+                                 psVector **modelSums // Buffer of model sums
+    )
 {
-    assert(kernel && kernel->kernel);
-
-    int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x
-    int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y
-
-    float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment)
-    float sum = 0.0;       // Sum (zero-th moment)
-    for (int y = yMin; y <= yMax; y++) {
-        for (int x = xMin; x <= xMax; x++) {
-            xCentroid += kernel->kernel[y][x] * x;
-            yCentroid += kernel->kernel[y][x] * y;
-            sum += kernel->kernel[y][x];
-        }
-    }
-    xCentroid /= sum;
-    yCentroid /= sum;
-
-    float eta20 = 0.0, eta02 = 0.0;     // Second moments
-    for (int y = yMin; y <= yMax; y++) {
-        float yDiff = y - yCentroid;
-        for (int x = xMin; x <= xMax; x++) {
-            float xDiff = x - xCentroid;
-            eta20 += PS_SQR(xDiff) * kernel->kernel[y][x];
-            eta02 += PS_SQR(yDiff) * kernel->kernel[y][x];
-        }
-    }
-
-    // Normalise to calculate the scale-invariant
-    float sum2 = PS_SQR(sum);
-    eta20 /= sum2;
-    eta02 /= sum2;
-    // eta11 /= sum2;
-
-    return eta20 + eta02;
+    assert(kernel);
+    assert(models);
+    assert(modelSums);
+
+    int xMin = kernel->xMin, xMax = kernel->xMax; // Bounds in x
+    int yMin = kernel->yMin, yMax = kernel->yMax; // Bounds in y
+
+    // Generate models
+    if (!*models) {
+        assert(!*modelSums);
+        *models = psArrayAlloc(size);
+        *modelSums = psVectorAlloc(size, PS_TYPE_F64);
+        for (int sigma = 0; sigma < size; sigma++) {
+            psKernel *model = psKernelAlloc(xMin, xMax, yMin, yMax); // Gaussian model
+            float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared
+            double sumGG = 0.0;         // Sum of square of Gaussian
+            for (int y = yMin; y <= yMax; y++) {
+                int y2 = PS_SQR(y);     // y squared
+                for (int x = xMin; x <= xMax; x++) {
+                    float rad2 = PS_SQR(x) + y2; // Radius squared
+                    float value = expf(-rad2 * invSigma2); // Model value
+                    model->kernel[y][x] = value;
+                    sumGG += PS_SQR(value);
+                }
+            }
+            (*models)->data[sigma] = model;
+            (*modelSums)->data.F64[sigma] = sumGG;
+        }
+    }
+
+    // Fit gaussians of varying widths to the image, record the chi^2
+    psVector *chi2 = psVectorAlloc(size, PS_TYPE_F32); // chi^2 as a function of radius
+    for (int sigma = 0; sigma < size; sigma++) {
+        double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian
+        psKernel *model = (*models)->data[sigma]; // Model of interest
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg);
+            }
+        }
+        float norm = sumFG / (*modelSums)->data.F64[sigma]; // Normalisation for Gaussian
+        double sumDev2 = 0.0;           // Sum of square deviations
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation
+                sumDev2 += PS_SQR(dev);
+            }
+        }
+        chi2->data.F32[sigma] = sumDev2;
+    }
+
+    // Find the minimum chi^2
+    int bestIndex = -1;                 // Index of best chi^2
+    float bestChi2 = INFINITY;          // Best chi^2
+    for (int i = 0; i < size; i++) {
+        if (chi2->data.F32[i] < bestChi2) {
+            bestChi2 = chi2->data.F32[i];
+            bestIndex = i;
+        }
+    }
+    psFree(chi2);
+
+    return bestIndex + 1;
 }
 
-#if 0
-// Calculate the deviations for a particular subtraction order
-static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel
-                                           pmSubtractionStampList *stamps, // Stamps to convolve
-                                           const pmSubtractionKernels *kernels, // Kernel basis functions
-                                           int footprint, // Stamp footprint
-                                           pmSubtractionMode mode // Mode of subtraction
-                                           )
+pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, float bg1, float bg2)
 {
-    assert(stamps);
-    assert(footprint >= 0);
-    assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2);
-
-    if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-        return NULL;
-    }
-
-    psVector *solution = pmSubtractionSolveEquation(NULL, stamps);
-    if (!solution) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-        return NULL;
-    }
-
-    if (sumKernel) {
-        float sum = 0.0;                // Sum of the kernel
-        psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel
-        for (int y = 0; y < image->numRows; y++) {
-            for (int x = 0; x < image->numCols; x++) {
-                sum += image->data.F32[y][x];
-            }
-        }
-        psFree(image);
-        *sumKernel = sum;
-    }
-
-    psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode);
-    psFree(solution);
-    if (!deviations) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
-        return NULL;
-    }
-
-    return deviations;
-}
-#endif
-
-pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius)
-{
-    PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR);
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, PM_SUBTRACTION_MODE_ERR);
 
     psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps
-    psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments
+    psVector *ratios = psVectorAlloc(stamps->num, PS_TYPE_F32); // Ratios of widths
+    psArray *models = NULL;             // Gaussian models
+    psVector *modelSums = NULL;         // Gaussian model sums
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -630,20 +645,31 @@
         }
         mask->data.PS_TYPE_MASK_DATA[i] = 0;
-        moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) /
-            subtractionOrderMoment(stamp->image2, radius);
-    }
+
+        // Widths of stars
+        int width1 = subtractionOrderWidth(stamp->image1, bg1, stamps->footprint, &models, &modelSums);
+        int width2 = subtractionOrderWidth(stamp->image2, bg2, stamps->footprint, &models, &modelSums);
+
+        if (width1 == -1 || width2 == -2) {
+            ratios->data.F32[i] = NAN;
+            mask->data.PS_TYPE_MASK_DATA[i] = 0xff;
+        } else {
+            ratios->data.F32[i] = (float)width1 / (float)width2;
+        }
+    }
+    psFree(models);
+    psFree(modelSums);
 
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-    if (!psVectorStats(stats, moments, NULL, mask, 0xff)) {
+    if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio.");
         psFree(mask);
-        psFree(moments);
+        psFree(ratios);
         psFree(stats);
         return PM_SUBTRACTION_MODE_ERR;
     }
-    psFree(moments);
+    psFree(ratios);
     psFree(mask);
 
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median moment: %lf", stats->robustMedian);
     pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2);
     psFree(stats);