Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 19058) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 19059) @@ -21,4 +21,5 @@ #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" +#include "pmSubtractionThreads.h" #include "pmSubtraction.h" @@ -28,4 +29,5 @@ #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,14 +215,31 @@ 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,5 +252,7 @@ psFree(subConv); psFree(convolved); + psMutexLock(target); psFree(subTarget); + psMutexUnlock(target); return; } @@ -285,4 +306,5 @@ // 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,4 +733,89 @@ +// 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,4 +874,6 @@ } + bool threaded = pmSubtractionThreaded(); // Running threaded? + // Outputs psImage *convImage1 = NULL, *convWeight1 = NULL; // Convolved image and weight from input 1 @@ -773,8 +882,14 @@ 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,8 +902,14 @@ 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,9 +954,4 @@ 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,50 +974,55 @@ } +#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 Index: /trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.h (revision 19058) +++ /trunk/psModules/src/imcombine/pmSubtraction.h (revision 19059) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-06-19 00:54:09 $ + * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-14 03:25:55 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -94,4 +94,8 @@ ); +/// Execute a thread job to convolve a patch of the image +bool pmSubtractionConvolveThread(const psThreadJob *job ///< Job to execute + ); + /// Convolve image in preparation for subtraction bool pmSubtractionConvolve(pmReadout *out1, ///< Output image 1 Index: /trunk/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 19058) +++ /trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 19059) @@ -8,4 +8,5 @@ #include "pmSubtractionMatch.h" #include "pmSubtractionEquation.h" +#include "pmSubtraction.h" bool threaded = false; // Run with threads? @@ -16,11 +17,49 @@ } -void pmSubtractionThreadSetup(void) +// Initialise a mutex in each of the input +static void subtractionMutexInit(pmReadout *ro) +{ + if (!ro) { + return; + } + + if (ro->image) { + psMutexInit(ro->image); + } + if (ro->weight) { + psMutexInit(ro->weight); + } + + return; +} + +static void subtractionMutexDestroy(pmReadout *ro) +{ + if (!ro) { + return; + } + + if (ro->image) { + psMutexDestroy(ro->image); + } + if (ro->weight) { + psMutexDestroy(ro->weight); + } + + return; +} + +void pmSubtractionThreadsInit(pmReadout *out1, pmReadout *out2, pmReadout *in1, pmReadout *in2) { if (threaded) { - psAbort("%s called multiple times.", __func__); + psAbort("Already running threaded."); } threaded = true; + + subtractionMutexInit(out1); + subtractionMutexInit(out2); + subtractionMutexInit(in1); + subtractionMutexInit(in2); { @@ -38,5 +77,30 @@ } + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CONVOLVE", 17); + task->function = &pmSubtractionConvolveThread; + psThreadTaskAdd(task); + psFree(task); + } + return; } + +void pmSubtractionThreadsFinalize(pmReadout *out1, pmReadout *out2, pmReadout *in1, pmReadout *in2) +{ + if (!threaded) { + return; + } + + threaded = false; + psThreadTaskRemove("PSMODULES_SUBTRACTION_ORDER"); + psThreadTaskRemove("PSMODULES_SUBTRACTION_CALCULATE_EQUATION"); + psThreadTaskRemove("PSMODULES_SUBTRACTION_CONVOLVE"); + + subtractionMutexDestroy(out1); + subtractionMutexDestroy(out2); + subtractionMutexDestroy(in1); + subtractionMutexDestroy(in1); + return; +} Index: /trunk/psModules/src/imcombine/pmSubtractionThreads.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionThreads.h (revision 19058) +++ /trunk/psModules/src/imcombine/pmSubtractionThreads.h (revision 19059) @@ -6,5 +6,15 @@ /// Set up threading for image matching -void pmSubtractionThreadSetup(void); +/// +/// Sets up thread tasks, and initialises mutexes in readouts +void pmSubtractionThreadsInit(pmReadout *out1, pmReadout *out2, // Output readouts + pmReadout *in1, pmReadout *in2 // Input readouts + ); + + +/// Take down threading for image matching +/// +/// Destroys thread tasks, and initialises mutexes in readouts +void pmSubtractionThreadsFinalize(pmReadout *out1, pmReadout *out2, pmReadout *in1, pmReadout *in2); #endif