Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 18348)
+++ trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 18962)
@@ -16,5 +16,7 @@
 #include "pmSubtraction.h"
 #include "pmSubtractionMask.h"
+#include "pmSubtractionThreads.h"
 #include "pmSubtractionMatch.h"
+
 
 #define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
@@ -631,6 +633,6 @@
                                  float bg, // Background in image
                                  int size, // Maximum size
-                                 psArray **models, // Buffer of models
-                                 psVector **modelSums // Buffer of model sums
+                                 const psArray *models, // Buffer of models
+                                 const psVector *modelSums // Buffer of model sums
     )
 {
@@ -639,29 +641,6 @@
     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;
-        }
-    }
+    int xMin = -size, xMax = size; // Bounds in x
+    int yMin = -size, yMax = size; // Bounds in y
 
     // Fit gaussians of varying widths to the image, record the chi^2
@@ -669,5 +648,5 @@
     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
+        psKernel *model = models->data[sigma]; // Model of interest
         for (int y = yMin; y <= yMax; y++) {
             for (int x = xMin; x <= xMax; x++) {
@@ -675,5 +654,5 @@
             }
         }
-        float norm = sumFG / (*modelSums)->data.F64[sigma]; // Normalisation for Gaussian
+        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++) {
@@ -700,4 +679,54 @@
 }
 
+
+bool pmSubtractionOrderStamp(psVector *ratios, psVector *mask, const pmSubtractionStampList *stamps,
+                             const psArray *models, const psVector *modelSums,
+                             int index, float bg1, float bg2)
+{
+    PS_ASSERT_VECTOR_NON_NULL(ratios, false);
+    PS_ASSERT_VECTOR_NON_NULL(mask, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(ratios, mask, false);
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PS_ASSERT_INT_NONNEGATIVE(index, false);
+    PS_ASSERT_INT_LESS_THAN(index, stamps->num, false);
+    PS_ASSERT_ARRAY_NON_NULL(models, false);
+    PS_ASSERT_VECTOR_NON_NULL(modelSums, false);
+    PS_ASSERT_VECTOR_SIZE(modelSums, models->n, false);
+
+    pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
+    psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE || stamp->status == PM_SUBTRACTION_STAMP_USED,
+             "We checked this earlier.");
+
+    // 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 == 0 || width2 == 0) {
+        ratios->data.F32[index] = NAN;
+        mask->data.PS_TYPE_MASK_DATA[index] = 0xff;
+    } else {
+        ratios->data.F32[index] = (float)width1 / (float)width2;
+        mask->data.PS_TYPE_MASK_DATA[index] = 0;
+    }
+
+    return true;
+}
+
+bool pmSubtractionOrderThread(const psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    psVector *ratios = job->args->data[0]; // Ratios of widths
+    psVector *mask = job->args->data[1]; // Mask for ratios
+    const pmSubtractionStampList *stamps = job->args->data[2]; // List of stamps
+    const psArray *models = job->args->data[3]; // Gaussian models
+    const psVector *modelSums = job->args->data[4]; // Gaussian model sums
+    int index = PS_SCALAR_VALUE(job->args->data[5], S32); // Stamp index
+    float bg1 = PS_SCALAR_VALUE(job->args->data[6], F32); // Background of image 1
+    float bg2 = PS_SCALAR_VALUE(job->args->data[7], F32); // Background of image 2
+
+    return pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, index, bg1, bg2);
+}
+
 pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, float bg1, float bg2)
 {
@@ -706,6 +735,27 @@
     psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps
     psVector *ratios = psVectorAlloc(stamps->num, PS_TYPE_F32); // Ratios of widths
-    psArray *models = NULL;             // Gaussian models
-    psVector *modelSums = NULL;         // Gaussian model sums
+
+    // Generate models
+    int size = stamps->footprint;       // Maximum size
+    psArray *models = psArrayAlloc(size); // Gaussian models
+    psVector *modelSums = psVectorAlloc(size, PS_TYPE_F64); // Gaussian model sums
+    for (int sigma = 0; sigma < size; sigma++) {
+        psKernel *model = psKernelAlloc(-size, size, -size, size); // 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 = -size; y <= size; y++) {
+            int y2 = PS_SQR(y);     // y squared
+            for (int x = -size; x <= size; 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] = 1.0 / sumGG;
+    }
+
+    // Fit models to stamps
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -715,16 +765,40 @@
         }
 
-        // 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 == 0 || width2 == 0) {
-            ratios->data.F32[i] = NAN;
-            mask->data.PS_TYPE_MASK_DATA[i] = 0xff;
+        if (pmSubtractionThreaded()) {
+            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER");
+            psArrayAdd(job->args, 1, ratios);
+            psArrayAdd(job->args, 1, mask);
+            psArrayAdd(job->args, 1, stamps);
+            psArrayAdd(job->args, 1, models);
+            psArrayAdd(job->args, 1, modelSums);
+            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                return false;
+            }
+            psFree(job);
         } else {
-            ratios->data.F32[i] = (float)width1 / (float)width2;
-            mask->data.PS_TYPE_MASK_DATA[i] = 0;
-        }
-    }
+            if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to measure PSF width for stamp %d", i);
+                psFree(models);
+                psFree(modelSums);
+                psFree(ratios);
+                psFree(mask);
+                return false;
+            }
+        }
+    }
+
+    if (!psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        psFree(models);
+        psFree(modelSums);
+        psFree(ratios);
+        psFree(mask);
+            return false;
+    }
+
     psFree(models);
     psFree(modelSums);