Index: /trunk/psModules/src/imcombine/Makefile.am =================================================================== --- /trunk/psModules/src/imcombine/Makefile.am (revision 18961) +++ /trunk/psModules/src/imcombine/Makefile.am (revision 18962) @@ -16,4 +16,5 @@ pmSubtractionParams.c \ pmSubtractionStamps.c \ + pmSubtractionThreads.c \ pmPSFEnvelope.c @@ -31,4 +32,5 @@ pmSubtractionParams.h \ pmSubtractionStamps.h \ + pmSubtractionThreads.h \ pmPSFEnvelope.h Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 18961) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 18962) @@ -550,4 +550,5 @@ rms = 0.74 * (stats->sampleUQ - stats->sampleLQ); } + psFree(stats); psTrace("psModules.imcombine", 1, "Mean: %f\n", mean); Index: /trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 18961) +++ /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 18962) @@ -10,4 +10,5 @@ #include "pmSubtractionKernels.h" #include "pmSubtractionStamps.h" +#include "pmSubtractionThreads.h" #include "pmSubtractionEquation.h" @@ -479,8 +480,22 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -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,5 +509,159 @@ 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,145 +673,23 @@ } - // 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; Index: /trunk/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 18961) +++ /trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 18962) @@ -4,4 +4,14 @@ #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 Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 18961) +++ /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); Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 18961) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 18962) @@ -15,4 +15,5 @@ #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,7 +52,22 @@ ); +/// 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 ); Index: /trunk/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 18962) +++ /trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 18962) @@ -0,0 +1,42 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmSubtractionMatch.h" +#include "pmSubtractionEquation.h" + +bool threaded = false; // Run with threads? + +bool pmSubtractionThreaded(void) +{ + return threaded; +} + +void pmSubtractionThreadSetup(void) +{ + if (threaded) { + psAbort("%s called multiple times.", __func__); + } + + threaded = true; + + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_ORDER", 8); + task->function = &pmSubtractionOrderThread; + psThreadTaskAdd(task); + psFree(task); + } + + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 3); + task->function = &pmSubtractionCalculateEquationThread; + psThreadTaskAdd(task); + psFree(task); + } + + return; +} + Index: /trunk/psModules/src/imcombine/pmSubtractionThreads.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionThreads.h (revision 18962) +++ /trunk/psModules/src/imcombine/pmSubtractionThreads.h (revision 18962) @@ -0,0 +1,10 @@ +#ifndef PM_SUBTRACTION_THREADS_H +#define PM_SUBTRACTION_THREADS_H + +/// Return whether theads have been activated +bool pmSubtractionThreaded(void); + +/// Set up threading for image matching +void pmSubtractionThreadSetup(void); + +#endif Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 18961) +++ /trunk/psModules/src/psmodules.h (revision 18962) @@ -85,4 +85,5 @@ #include #include +#include #include #include