Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 18287) +++ 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;