Changeset 18962

Timestamp:

Aug 8, 2008, 8:10:29 AM (18 years ago)

Author:

Paul Price

Message:

Adding threads to pmSubtractionCalculateEquation and pmSubtractionOrder. More work needs to be done in threading the final convolution.

Location:

trunk/psModules/src

Files:

• imcombine/Makefile.am (modified) (2 diffs)
• imcombine/pmSubtraction.c (modified) (1 diff)
• imcombine/pmSubtractionEquation.c (modified) (4 diffs)
• imcombine/pmSubtractionEquation.h (modified) (1 diff)
• imcombine/pmSubtractionMatch.c (modified) (8 diffs)
• imcombine/pmSubtractionMatch.h (modified) (2 diffs)
• imcombine/pmSubtractionThreads.c (added)
• imcombine/pmSubtractionThreads.h (added)
• psmodules.h (modified) (1 diff)

trunk/psModules/src/imcombine/Makefile.am

@@ -16 +16 @@
         pmSubtractionParams.c   \
         pmSubtractionStamps.c   \
+        pmSubtractionThreads.c  \
         pmPSFEnvelope.c
 
@@ -31 +32 @@
         pmSubtractionParams.h   \
         pmSubtractionStamps.h   \
+        pmSubtractionThreads.h  \
         pmPSFEnvelope.h
 

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -550 +550 @@
         rms = 0.74 * (stats->sampleUQ - stats->sampleLQ);
     }
+    psFree(stats);
 
     psTrace("psModules.imcombine", 1, "Mean: %f\n", mean);

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -10 +10 @@
 #include "pmSubtractionKernels.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionThreads.h"
 
 #include "pmSubtractionEquation.h"
@@ -479 +480 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels)
+bool pmSubtractionCalculateEquationThread(const psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    pmSubtractionStampList *stamps = job->args->data[0]; // List of stamps
+    const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
+    int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
+
+    return pmSubtractionCalculateEquationStamp(stamps, kernels, index);
+}
+
+bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
+                                         int index)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
+    PS_ASSERT_INT_NONNEGATIVE(index, false);
+    PS_ASSERT_INT_LESS_THAN(index, stamps->num, false);
 
     int footprint = stamps->footprint;  // Half-size of stamps
@@ -494 +509 @@
     int numParams = numKernels * numSpatial + 1 + numBackground;
 
-    psImage *polyValues = NULL;         // Polynomial terms
+    pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
+    psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state.");
+
+    // Generate convolutions
+    if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", index);
+        return NULL;
+    }
+
+#ifdef TESTING
+    for (int j = 0; j < numKernels; j++) {
+        if (stamp->convolutions1) {
+            psString convName = NULL;
+            psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
+            psFits *fits = psFitsOpen(convName, "w");
+            psFree(convName);
+            psKernel *conv = stamp->convolutions1->data[j];
+            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+            psFitsClose(fits);
+        }
+
+        if (stamp->convolutions2) {
+            psString convName = NULL;
+            psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
+            psFits *fits = psFitsOpen(convName, "w");
+            psFree(convName);
+            psKernel *conv = stamp->convolutions2->data[j];
+            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+            psFitsClose(fits);
+        }
+    }
+#endif
+
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
+                                                    stamp->xNorm, stamp->yNorm); // Polynomial terms
+
+    bool new = stamp->vector1 ? false : true; // Is this a new run?
+    if (new) {
+        stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+        stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64);
+    }
+#ifdef TESTING
+    psImageInit(stamp->matrix1, NAN);
+    psVectorInit(stamp->vector1, NAN);
+#endif
+
+    bool status;                    // Status of least-squares matrix/vector calculation
+    switch (kernels->mode) {
+      case PM_SUBTRACTION_MODE_1:
+        status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
+                                 stamp->weight, polyValues, footprint, true);
+        status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
+                                  stamp->image2, stamp->weight, polyValues, footprint, true);
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2,
+                                 stamp->weight, polyValues, footprint, true);
+        status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2,
+                                  stamp->image1, stamp->weight, polyValues, footprint, true);
+        break;
+      case PM_SUBTRACTION_MODE_DUAL:
+        if (new) {
+            stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
+            stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64);
+            stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64);
+        }
+#ifdef TESTING
+        psImageInit(stamp->matrix2, NAN);
+        psImageInit(stamp->matrixX, NAN);
+        psVectorInit(stamp->vector2, NAN);
+#endif
+        status  = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
+                                  stamp->weight, polyValues, footprint, true);
+        status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL,
+                                  stamp->weight, polyValues, footprint, false);
+        status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1,
+                                       stamp->convolutions2, stamp->image1, stamp->weight, polyValues,
+                                       footprint);
+        status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
+                                  stamp->image2, stamp->weight, polyValues, footprint, true);
+        status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL,
+                                  stamp->image2, stamp->weight, polyValues, footprint, false);
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", kernels->mode);
+    }
+
+    if (!status) {
+        stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+        psWarning("Rejecting stamp %d (%d,%d) because of bad equation",
+                  index, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+    } else {
+        stamp->status = PM_SUBTRACTION_STAMP_USED;
+    }
+
+#ifdef TESTING
+    if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
+        psString matrixName = NULL;
+        psStringAppend(&matrixName, "matrix1_%d.fits", index);
+        psFits *matrixFile = psFitsOpen(matrixName, "w");
+        psFree(matrixName);
+        psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL);
+        psFitsClose(matrixFile);
+
+        matrixName = NULL;
+        psStringAppend(&matrixName, "vector1_%d.fits", index);
+        psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64);
+        memcpy(dummy->data.F64[0], stamp->vector1->data.F64,
+               PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n);
+        matrixFile = psFitsOpen(matrixName, "w");
+        psFree(matrixName);
+        psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
+        psFree(dummy);
+        psFitsClose(matrixFile);
+
+        if (stamp->vector2) {
+            matrixName = NULL;
+            psStringAppend(&matrixName, "vector2_%d.fits", index);
+            dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);
+            memcpy(dummy->data.F64[0], stamp->vector2->data.F64,
+                   PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);
+            matrixFile = psFitsOpen(matrixName, "w");
+            psFree(matrixName);
+            psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
+            psFree(dummy);
+            psFitsClose(matrixFile);
+        }
+
+        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+            matrixName = NULL;
+            psStringAppend(&matrixName, "matrix2_%d.fits", index);
+            matrixFile = psFitsOpen(matrixName, "w");
+            psFree(matrixName);
+            psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL);
+            psFitsClose(matrixFile);
+
+            matrixName = NULL;
+            psStringAppend(&matrixName, "matrixX_%d.fits", index);
+            matrixFile = psFitsOpen(matrixName, "w");
+            psFree(matrixName);
+            psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL);
+            psFitsClose(matrixFile);
+        }
+    }
+#endif
+
+    psFree(polyValues);
+
+    return true;
+}
+
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
 
     // We iterate over each stamp, allocate the matrix and vectors if
@@ -504 +673 @@
         }
 
-        // Generate convolutions
-        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
-            psFree(polyValues);
-            return NULL;
-        }
-
-#ifdef TESTING
-        for (int j = 0; j < numKernels; j++) {
-            if (stamp->convolutions1) {
-                psString convName = NULL;
-                psStringAppend(&convName, "conv1_%03d_%03d.fits", i, j);
-                psFits *fits = psFitsOpen(convName, "w");
-                psFree(convName);
-                psKernel *conv = stamp->convolutions1->data[j];
-                psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-                psFitsClose(fits);
-            }
-
-            if (stamp->convolutions2) {
-                psString convName = NULL;
-                psStringAppend(&convName, "conv2_%03d_%03d.fits", i, j);
-                psFits *fits = psFitsOpen(convName, "w");
-                psFree(convName);
-                psKernel *conv = stamp->convolutions2->data[j];
-                psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-                psFitsClose(fits);
-            }
-        }
-#endif
-
-        polyValues = p_pmSubtractionPolynomial(polyValues, spatialOrder, stamp->xNorm, stamp->yNorm);
-
-
-            stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
-            stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64);
-#ifdef TESTING
-            psImageInit(stamp->matrix1, NAN);
-            psVectorInit(stamp->vector1, NAN);
-#endif
-
-        bool status;                    // Status of least-squares matrix/vector calculation
-        switch (kernels->mode) {
-          case PM_SUBTRACTION_MODE_1:
-            status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
-                                     stamp->weight, polyValues, footprint, true);
-            status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
-                                      stamp->image2, stamp->weight, polyValues, footprint, true);
-            break;
-          case PM_SUBTRACTION_MODE_2:
-            status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2,
-                                     stamp->weight, polyValues, footprint, true);
-            status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2,
-                                      stamp->image1, stamp->weight, polyValues, footprint, true);
-            break;
-          case PM_SUBTRACTION_MODE_DUAL:
-            stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
-            stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64);
-            stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64);
-#ifdef TESTING
-            psImageInit(stamp->matrix2, NAN);
-            psImageInit(stamp->matrixX, NAN);
-            psVectorInit(stamp->vector2, NAN);
-#endif
-            status  = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
-                                      stamp->weight, polyValues, footprint, true);
-            status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL,
-                                      stamp->weight, polyValues, footprint, false);
-            status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1,
-                                           stamp->convolutions2, stamp->image1, stamp->weight, polyValues,
-                                           footprint);
-            status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
-                                      stamp->image2, stamp->weight, polyValues, footprint, true);
-            status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL,
-                                      stamp->image2, stamp->weight, polyValues, footprint, false);
-            break;
-          default:
-            psAbort("Unsupported subtraction mode: %x", kernels->mode);
-        }
-
-        if (!status) {
-            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            psWarning("Rejecting stamp %d (%d,%d) because of bad equation",
-                      i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+        if (pmSubtractionThreaded()) {
+            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION");
+            psArrayAdd(job->args, 1, stamps);
+            psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
+            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                return false;
+            }
+            psFree(job);
         } else {
-            stamp->status = PM_SUBTRACTION_STAMP_USED;
-        }
-
-#ifdef TESTING
-        if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
-            psString matrixName = NULL;
-            psStringAppend(&matrixName, "matrix1_%d.fits", i);
-            psFits *matrixFile = psFitsOpen(matrixName, "w");
-            psFree(matrixName);
-            psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL);
-            psFitsClose(matrixFile);
-
-            matrixName = NULL;
-            psStringAppend(&matrixName, "vector1_%d.fits", i);
-            psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64);
-            memcpy(dummy->data.F64[0], stamp->vector1->data.F64,
-                   PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n);
-            matrixFile = psFitsOpen(matrixName, "w");
-            psFree(matrixName);
-            psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
-            psFree(dummy);
-            psFitsClose(matrixFile);
-
-            if (stamp->vector2) {
-                matrixName = NULL;
-                psStringAppend(&matrixName, "vector2_%d.fits", i);
-                dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);
-                memcpy(dummy->data.F64[0], stamp->vector2->data.F64,
-                       PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);
-                matrixFile = psFitsOpen(matrixName, "w");
-                psFree(matrixName);
-                psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
-                psFree(dummy);
-                psFitsClose(matrixFile);
-            }
-
-            if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-                matrixName = NULL;
-                psStringAppend(&matrixName, "matrix2_%d.fits", i);
-                matrixFile = psFitsOpen(matrixName, "w");
-                psFree(matrixName);
-                psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL);
-                psFitsClose(matrixFile);
-
-                matrixName = NULL;
-                psStringAppend(&matrixName, "matrixX_%d.fits", i);
-                matrixFile = psFitsOpen(matrixName, "w");
-                psFree(matrixName);
-                psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL);
-                psFitsClose(matrixFile);
-            }
-        }
-#endif
-
-    }
-    psFree(polyValues);
+            pmSubtractionCalculateEquationStamp(stamps, kernels, i);
+        }
+    }
+
+    if (!psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        return false;
+    }
 
     return true;

trunk/psModules/src/imcombine/pmSubtractionEquation.h

@@ -4 +4 @@
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionKernels.h"
+
+/// Execute a thread job to calculate the least-squares equation for a stamp
+bool pmSubtractionCalculateEquationThread(const psThreadJob *job ///< Job to execute
+    );
+
+/// Calculate the least-squares equation to match the image quality for a single stamp
+bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps
+                                         const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                         int index ///< Index of stamp
+                                    );
 
 /// Calculate the least-squares equation to match the image quality

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -16 +16 @@
 #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 +633 @@
                                  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 +641 @@
     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 +648 @@
     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 +654 @@
             }
         }
-        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 +679 @@
 }
 
+
+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 +735 @@
     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 +765 @@
         }
 
-        // 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);

trunk/psModules/src/imcombine/pmSubtractionMatch.h

@@ -15 +15 @@
 #define PM_SUBTRACTION_ANALYSIS_STAMPS_NUM "SUBTRACTION.NUM"// Name of the number of stamps in the analysis MD
 #define PM_SUBTRACTION_ANALYSIS_VARFACTOR "SUBTRACTION.VARFACTOR"// Name of variance factor in the analysis MD
+
 
 /// Match two images
@@ -51 +52 @@
     );
 
+/// Execute a thread job to measure the PSF width ratios
+bool pmSubtractionOrderThread(const psThreadJob *job ///< Job to execute
+    );
+
+/// Measure the PSF width ratio for a single stamp
+bool pmSubtractionOrderStamp(psVector *ratios, ///< PSF width ratios
+                             psVector *mask, ///< Mask for PSF width ratios
+                             const pmSubtractionStampList *stamps, ///< List of stamps
+                             const psArray *models, ///< Pre-calculated gaussian models
+                             const psVector *modelSums, ///< Pre-calculated gaussian model sums
+                             int index, ///< Index of stamp
+                             float bg1, ///< Background for image 1
+                             float bg2  ///< Background for image 2
+    );
+
 /// Determine which image to convolve
 pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, ///< Stamps that have been extracted
-                                     float bg1, float bg2 // Background for each image
+                                     float bg1, float bg2 ///< Background for each image
     );
 

trunk/psModules/src/psmodules.h

@@ -85 +85 @@
 #include <pmStackReject.h>
 #include <pmSubtraction.h>
+#include <pmSubtractionThreads.h>
 #include <pmSubtractionStamps.h>
 #include <pmSubtractionKernels.h>