Index: trunk/psModules/src/imcombine/pmStackReject.c
===================================================================
--- trunk/psModules/src/imcombine/pmStackReject.c	(revision 19282)
+++ trunk/psModules/src/imcombine/pmStackReject.c	(revision 19345)
@@ -8,4 +8,5 @@
 
 #include "pmSubtraction.h"
+#include "pmSubtractionThreads.h"
 #include "pmSubtractionKernels.h"
 
@@ -13,4 +14,103 @@
 
 //#define TESTING                         // Testing output
+
+// Mask values
+typedef enum {
+    PM_STACK_MASK_BAD      = 0x01,      // Bad pixel
+    PM_STACK_MASK_CONVOLVE = 0x02,      // Touching a bad pixel
+} pmStackMask;
+
+static bool threaded = false;           // Running threaded?
+
+
+// Grow the rejection mask
+static inline bool stackRejectGrow(psImage *target,   // Target mask image (product)
+                                   psImage *source, // Source mask image (to be grown)
+                                   const pmSubtractionKernels *kernels, // Subtraction kernels
+                                   int numCols, int numRows, // Size of image
+                                   int xMin, int xMax, int yMin, int yMax, // Bounds of convolution
+                                   float poorFrac       // Fraction for "poor"
+    )
+{
+    int size = kernels->size;           // Half-size of convolution kernel
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
+                                                              xMin + size + 1, yMin + size + 1); // Polynomial
+    int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues, false, poorFrac); // Radius of bad box
+    psFree(polyValues);
+
+    if (box > 0) {
+        // Convolve a subimage, then stick it in the target
+        if (threaded) {
+            psMutexLock(source);
+        }
+        psImage *mask = psImageSubset(source, psRegionSet(xMin - box, xMax + box,
+                                                          yMin - box, yMax + box)); // Mask to convolve
+        if (threaded) {
+            psMutexUnlock(source);
+        }
+        psImage *convolved = psImageConvolveMask(NULL, mask, PM_STACK_MASK_BAD, PM_STACK_MASK_CONVOLVE,
+                                                 -box, box, -box, box); // Convolved mask
+        if (threaded) {
+            psMutexLock(source);
+        }
+        psFree(mask);
+        if (threaded) {
+            psMutexUnlock(source);
+        }
+
+        int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy
+        psAssert(convolved->numCols - 2 * box == xMax - xMin, "Bad number of columns");
+        psAssert(convolved->numRows - 2 * box == yMax - yMin, "Bad number of rows");
+
+        for (int yTarget = yMin, ySource = box; yTarget < yMax; yTarget++, ySource++) {
+            memcpy(&target->data.PS_TYPE_MASK_DATA[yTarget][xMin],
+                   &convolved->data.PS_TYPE_MASK_DATA[ySource][box], numBytes);
+        }
+        psFree(convolved);
+    }
+    return true;
+}
+
+// Thread entry for stackRejectGrow
+static bool stackRejectGrowThread(psThreadJob *job // Job to execute
+    )
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    psArray *args = job->args;          // Job arguments
+    psImage *target = args->data[0];    // Target mask image
+    psImage *source = args->data[1];    // Source mask image
+    const pmSubtractionKernels *kernels = args->data[2]; // Subtraction kernels
+    int numCols = PS_SCALAR_VALUE(args->data[3], S32); // Number of columns
+    int numRows = PS_SCALAR_VALUE(args->data[4], S32); // Number of rows
+    int xMin = PS_SCALAR_VALUE(args->data[5], S32); // Minimum x value
+    int xMax = PS_SCALAR_VALUE(args->data[6], S32); // Maximum x value
+    int yMin = PS_SCALAR_VALUE(args->data[7], S32); // Minimum y value
+    int yMax = PS_SCALAR_VALUE(args->data[8], S32); // Maximum y value
+    float poorFrac = PS_SCALAR_VALUE(args->data[9], F32); // Fraction for "poor"
+
+    return stackRejectGrow(target, source, kernels, numCols, numRows, xMin, xMax, yMin, yMax, poorFrac);
+}
+
+bool pmStackRejectThreadsInit(void)
+{
+    if (threaded) {
+        psAbort("Already running threaded.");
+    }
+    threaded = true;
+
+    if (!pmSubtractionThreaded()) {
+        pmSubtractionThreadsInit(NULL, NULL);
+    }
+
+    {
+        psThreadTask *task = psThreadTaskAlloc("PSMODULES_STACK_REJECT_GROW", 10);
+        task->function = &stackRejectGrowThread;
+        psThreadTaskAdd(task);
+        psFree(task);
+    }
+
+    return true;
+}
 
 
@@ -52,4 +152,7 @@
     pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
     inRO->image = image;
+    if (threaded) {
+        psMutexInit(image);
+    }
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = subRegions->data[i]; // Region of interest
@@ -70,6 +173,6 @@
         float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) /
             (float)kernels->numRows;
-        psImage *image = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
-        if (!image) {
+        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
+        if (!kernel) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
             psFree(convRO);
@@ -78,10 +181,10 @@
         }
         float sum = 0.0;
-        for (int y = 0; y < image->numRows; y++) {
-            for (int x = 0; x < image->numCols; x++) {
-                sum += image->data.F32[y][x];
+        for (int y = 0; y < kernel->numRows; y++) {
+            for (int x = 0; x < kernel->numCols; x++) {
+                sum += kernel->data.F32[y][x];
             }
         }
-        psFree(image);
+        psFree(kernel);
 
         // Range for normalisation
@@ -97,4 +200,7 @@
         }
     }
+    if (threaded) {
+        psMutexDestroy(image);
+    }
     psFree(inRO);
     psImage *convolved = psMemIncrRefCounter(convRO->image);
@@ -123,11 +229,12 @@
         }
     }
-    psFree(convolved);
 
     // Now, grow the mask to include everything that touches a bad pixel in the convolution
-    // Mask values:
-    // 0x0f: we think this is bad
-    // 0xf0: this is within the convolution kernel of a bad pixel
-    mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0x0f);
+    psImage *source = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1),
+                                     PM_STACK_MASK_BAD); // Mask image to grow
+    psImage *target = psImageRecycle(convolved, numCols, numRows, PS_TYPE_U8); // Grown image
+    if (threaded) {
+        psMutexInit(source);
+    }
     for (int i = 0; i < subRegions->n; i++) {
         psRegion *region = subRegions->data[i]; // Subtraction region
@@ -141,5 +248,4 @@
         int yMin = PS_MAX(region->y0, size), yMax = PS_MIN(region->y1, numRows - size);
 
-        psImage *polyValues = NULL;     // Pre-calculated polynomial values
         for (int j = yMin; j < yMax; j += fullSize) {
             int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
@@ -147,31 +253,49 @@
                 int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
 
-                polyValues = p_pmSubtractionPolynomialFromCoords(polyValues, kernels, numCols, numRows,
-                                                                 i + size + 1, j + size + 1);
-                int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues,
-                                                   false, poorFrac); // Radius of bad box
-                if (box > 0) {
-                    // Convolve a subimage, then stick it in the original
-                    psImage *subMask = psImageSubset(mask, psRegionSet(i - box, xSubMax + box,
-                                                                       j - box, ySubMax + box)); // Subimage
-                    psImage *convolved = psImageConvolveMask(NULL, subMask, 0x0f, 0xf0,
-                                                             -box, box, -box, box); // Convolved mask
-                    psFree(subMask);
-
-                    int numBytes = (xSubMax - i) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy
-                    psAssert(convolved->numCols - 2 * box == xSubMax - i, "Bad number of columns");
-                    psAssert(convolved->numRows - 2 * box == ySubMax - j, "Bad number of rows");
-
-                    for (int yTarget = j, ySource = box; yTarget < ySubMax; yTarget++, ySource++) {
-                        memcpy(&mask->data.PS_TYPE_MASK_DATA[yTarget][i],
-                               &convolved->data.PS_TYPE_MASK_DATA[ySource][box], numBytes);
+                if (threaded) {
+                    psThreadJob *job = psThreadJobAlloc("PSMODULES_STACK_REJECT_GROW"); // Job to execute
+                    psArray *args = job->args; // Job arguments
+                    psArrayAdd(args, 1, target);
+                    psMutexLock(source);
+                    psArrayAdd(args, 1, source);
+                    psMutexUnlock(source);
+                    psArrayAdd(args, 1, kernels);
+                    PS_ARRAY_ADD_SCALAR(args, numCols, PS_TYPE_S32);
+                    PS_ARRAY_ADD_SCALAR(args, numRows, PS_TYPE_S32);
+                    PS_ARRAY_ADD_SCALAR(args, i, PS_TYPE_S32);
+                    PS_ARRAY_ADD_SCALAR(args, xSubMax, PS_TYPE_S32);
+                    PS_ARRAY_ADD_SCALAR(args, j, PS_TYPE_S32);
+                    PS_ARRAY_ADD_SCALAR(args, ySubMax, PS_TYPE_S32);
+                    PS_ARRAY_ADD_SCALAR(args, poorFrac, PS_TYPE_F32);
+                    if (!psThreadJobAddPending(job)) {
+                        psFree(job);
+                        psFree(source);
+                        psFree(target);
+                        return false;
                     }
-                    psFree(convolved);
+                    psFree(job);
+                } else if (!stackRejectGrow(target, source, kernels, numCols, numRows,
+                                              i, xSubMax, j, ySubMax, poorFrac)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to grow bad pixels.");
+                    psFree(source);
+                    psFree(target);
+                    return NULL;
                 }
             }
         }
-        psFree(polyValues);
-    }
-    bad = psPixelsFromMask(bad, mask, 0xff);
+    }
+
+    if (!psThreadPoolWait(false)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to grow bad pixels.");
+        psFree(source);
+        psFree(target);
+        return NULL;
+    }
+
+    if (threaded) {
+        psMutexDestroy(source);
+    }
+    psFree(source);
+    bad = psPixelsFromMask(bad, target, 0xff);
 
     return bad;