Changeset 19059 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Aug 13, 2008, 5:25:55 PM (18 years ago)

Author:

Paul Price

Message:

Adding threading for pmSubtractionConvolve. This was more difficult than for the others because the threads are working on a common image. Turns out that making child images using psImageSubset is not thread-safe: the image (which all threads are using) is modified. Wrapped those calls up in mutexes, and now threading seems to work fine.

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (11 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -21 +21 @@
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionEquation.h"
+#include "pmSubtractionThreads.h"
 
 #include "pmSubtraction.h"
@@ -28 +29 @@
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
 #define MIN_SAMPLE_STATS    7           // Minimum number to use sample statistics; otherwise use quartiles
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -213 +215 @@
                                   region.y0 - size, region.y1 + size); // Add a border
 
-    // Casting away const so psImageSubset can add the child
+    // Casting away const
+    psMutexLock((psImage*)image);
     psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve
+    psMutexUnlock((psImage*)image);
+
+    // Casting away const
+    psMutexLock((psImage*)mask);
     psImage *subMask = mask ? psImageSubset((psImage*)mask, border) : NULL; // Subimage mask
+    psMutexUnlock((psImage*)mask);
+
     psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution
+
+    psMutexLock((psImage*)image);
     psFree(subImage);
+    psMutexUnlock((psImage*)image);
+
+    psMutexLock((psImage*)mask);
     psFree(subMask);
+    psMutexUnlock((psImage*)mask);
 
     // Now, we have to chop off the borders, and stick it in where it belongs
+    // No locking because we own this one
     psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges
+
+    psMutexLock(target);
     psImage *subTarget = psImageSubset(target, region); // Target for this subregion
+    psMutexUnlock(target);
     if (background != 0.0) {
         psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
@@ -233 +252 @@
     psFree(subConv);
     psFree(convolved);
+    psMutexLock(target);
     psFree(subTarget);
+    psMutexUnlock(target);
     return;
 }
@@ -285 +306 @@
         // Use Fast Fourier Transform to do the convolution
         // This provides a big speed-up for large kernels
+
         convolveFFT(convImage, image, subMask, maskVal, *kernelImage, region, background, kernels->size);
         if (weight) {
@@ -711 +733 @@
 
 
+// XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data
+static bool subtractionConvolvePatch(int numCols, int numRows, int x0, int y0,
+                                     pmReadout *out1, pmReadout *out2, psImage *convMask,
+                                     const pmReadout *ro1, const pmReadout *ro2, const psImage *subMask,
+                                     psMaskType blank, psMaskType maskSource, psMaskType maskTarget,
+                                     const psRegion *region, const pmSubtractionKernels *kernels,
+                                     bool doBG, bool useFFT)
+{
+    int size = kernels->size;           // Half-size of kernel
+    int xMin = region->x0, xMax = region->x1, yMin = region->y0, yMax = region->y1; // Bounds of patch
+
+    // Size to use when calculating normalised coordinates (different from actual size when convolving
+    // subimage)
+    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
+    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
+
+    // Normalised coordinates
+    float yNorm = 2.0 * (float)(yMin + y0 + size + 1 - yNormSize/2.0) / (float)yNormSize; // Normalised coord
+    float xNorm = 2.0 * (float)(xMin + x0 + size + 1 - xNormSize/2.0) / (float)xNormSize; // Normalised coord
+
+    psKernel *kernelImage = NULL;       // Kernel for the images
+    psKernel *kernelWeight = NULL;      // Kernel for the weight maps
+
+    // Only generate polynomial values every kernel footprint, since we have already assumed
+    // (with the stamps) that it does not vary rapidly on this scale.
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder,
+                                                    xNorm, yNorm); // Pre-calculated polynomial values
+    float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term
+
+    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        convolveRegion(out1->image, out1->weight, &kernelImage, &kernelWeight, ro1->image, ro1->weight,
+                       subMask, maskSource, kernels, polyValues, background, *region, useFFT,
+                       false);
+    }
+    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        convolveRegion(out2->image, out2->weight, &kernelImage, &kernelWeight, ro2->image, ro2->weight,
+                       subMask, maskSource, kernels, polyValues, background, *region, useFFT,
+                       kernels->mode == PM_SUBTRACTION_MODE_DUAL);
+    }
+
+    psFree(kernelImage);
+    psFree(kernelWeight);
+    psFree(polyValues);
+
+    // Propagate the mask
+    if (subMask) {
+        for (int y = yMin; y < yMax; y++) {
+            for (int x = xMin; x < xMax; x++) {
+                if (subMask->data.PS_TYPE_MASK_DATA[y][x] & maskTarget) {
+                    convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
+                }
+            }
+        }
+    }
+
+    return true;
+}
+
+bool pmSubtractionConvolveThread(const psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    psArray *args = job->args;          // Arguments
+    int numCols = PS_SCALAR_VALUE(args->data[0], S32); // Number of columns
+    int numRows = PS_SCALAR_VALUE(args->data[1], S32); // Number of rows
+    int x0 = PS_SCALAR_VALUE(args->data[2], S32); // Offset in x
+    int y0 = PS_SCALAR_VALUE(args->data[3], S32); // Offset in x
+    pmReadout *out1 = args->data[4];    // Output readout 1
+    pmReadout *out2 = args->data[5];    // Output readout 2
+    psImage *convMask = args->data[6];  // Output convolved mask
+    const pmReadout *ro1 = args->data[7]; // Input readout 1
+    const pmReadout *ro2 = args->data[8]; // Input readout 2
+    const psImage *subMask = args->data[9]; // Subtraction mask
+    psMaskType blank = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value
+    psMaskType maskSource = PS_SCALAR_VALUE(args->data[11], U8); // Mask value corresponding to source image
+    psMaskType maskTarget = PS_SCALAR_VALUE(args->data[12], U8); // Mask value corresponding to target image
+    const psRegion *region = args->data[13]; // Region to convolve
+    const pmSubtractionKernels *kernels = args->data[14]; // Kernels
+    bool doBG = PS_SCALAR_VALUE(args->data[15], U8); // Do background subtraction?
+    bool useFFT = PS_SCALAR_VALUE(args->data[16], U8); // Use FFT for convolution?
+
+    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask,
+                                    ro1, ro2, subMask, blank, maskSource, maskTarget,
+                                    region, kernels, doBG, useFFT);
+}
 
 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2,
@@ -767 +874 @@
     }
 
+    bool threaded = pmSubtractionThreaded(); // Running threaded?
+
     // Outputs
     psImage *convImage1 = NULL, *convWeight1 = NULL; // Convolved image and weight from input 1
@@ -773 +882 @@
         if (!convImage1) {
             convImage1 = out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            if (threaded) {
+                psMutexInit(convImage1);
+            }
         }
         if (weight1) {
             if (!out1->weight) {
                 out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+                if (threaded) {
+                    psMutexInit(out1->weight);
+                }
             }
             convWeight1 = out1->weight;
@@ -787 +902 @@
         if (!convImage2) {
             convImage2 = out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            if (threaded) {
+                psMutexInit(convImage2);
+            }
         }
         if (weight2) {
             if (!out2->weight) {
                 out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+                if (threaded) {
+                    psMutexInit(out2->weight);
+                }
             }
             convWeight2 = out2->weight;
@@ -833 +954 @@
         yMax = PS_MIN(region->y1, yMax);
     }
-
-    // Size to use when calculating normalised coordinates (different from actual size when convolving
-    // subimage)
-    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
-    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
 
     psMaskType maskSource;              // Mask these pixels when convolving
@@ -858 +974 @@
     }
 
+#if 0
+    // XXX Use thread-specific data to store these
     psImage *polyValues = NULL;         // Pre-calculated polynomial values
     psKernel *kernelImage = NULL;       // Kernel for the images
     psKernel *kernelWeight = NULL;      // Kernel for the weight maps
+#endif
 
     for (int j = yMin; j < yMax; j += fullSize) {
         int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
-        float yNorm = 2.0 * (float)(j + y0 + size + 1 - yNormSize/2.0) /
-            (float)yNormSize; // Normalised coordinate
         for (int i = xMin; i < xMax; i += fullSize) {
             int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
-            float xNorm = 2.0 * (float)(i + x0 + size + 1 - xNormSize/2.0) /
-                (float)xNormSize; // Normalised coordinate
-
-            // Only generate polynomial values every kernel footprint, since we have already assumed
-            // (with the stamps) that it does not vary rapidly on this scale.
-            polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, xNorm, yNorm);
-            float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) :
-                0.0; // Background term
-            psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve
-
-            if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-                convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image1, weight1,
-                               subMask, maskSource, kernels, polyValues, background, subRegion, useFFT,
-                               false);
+
+            psRegion *subRegion = psRegionAlloc(i, xSubMax, j, ySubMax); // Bounds of subtraction
+            if (threaded) {
+                psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CONVOLVE");
+                psArray *args = job->args;
+                PS_ARRAY_ADD_SCALAR(args, numCols, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(args, numRows, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(args, x0, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(args, y0, PS_TYPE_S32);
+                psArrayAdd(args, 1, out1);
+                psArrayAdd(args, 1, out2);
+                psArrayAdd(args, 1, convMask);
+                psArrayAdd(args, 1, (pmReadout*)ro1); // Casting away const
+                psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const
+                psArrayAdd(args, 1, (psImage*)subMask); // Casting away const
+                PS_ARRAY_ADD_SCALAR(args, blank, PS_TYPE_U8);
+                PS_ARRAY_ADD_SCALAR(args, maskSource, PS_TYPE_U8);
+                PS_ARRAY_ADD_SCALAR(args, maskTarget, PS_TYPE_U8);
+                psArrayAdd(args, 1, subRegion);
+                psArrayAdd(args, 1, (pmSubtractionKernels*)kernels); // Casting away const
+                PS_ARRAY_ADD_SCALAR(args, doBG, PS_TYPE_U8);
+                PS_ARRAY_ADD_SCALAR(args, useFFT, PS_TYPE_U8);
+
+                if (!psThreadJobAddPending(job)) {
+                    psFree(job);
+                    return false;
+                }
+                psFree(job);
+            } else {
+                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask,
+                                         blank, maskSource, maskTarget, subRegion, kernels, doBG, useFFT);
             }
-            if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-                convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, image2, weight2,
-                               subMask, maskSource, kernels, polyValues, background, subRegion, useFFT,
-                               kernels->mode == PM_SUBTRACTION_MODE_DUAL);
-            }
-
-            // Propagate the mask
-            if (subMask) {
-                for (int y = j; y < ySubMax; y++) {
-                    for (int x = i; x < xSubMax; x++) {
-                        if (subMask->data.PS_TYPE_MASK_DATA[y][x] & maskTarget) {
-                            convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
-                        }
-                    }
-                }
-            }
-        }
-    }
-    psFree(kernelImage);
-    psFree(kernelWeight);
-    psFree(polyValues);
+            psFree(subRegion);
+        }
+    }
+
+    if (!psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        return false;
+    }
 
     // Copy anything that wasn't convolved