Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 14420) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 14455) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-07 19:02:55 $ + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-09 20:25:52 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -30,4 +30,11 @@ // Private (file-static) functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Return the number of polynomial terms there are for a given order +static inline long polyTerms(int order // Order of polynomial + ) +{ + return (order + 1) * (order + 2) / 2; +} // Given a stamp coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j @@ -160,12 +167,8 @@ case PM_SUBTRACTION_KERNEL_ISIS: { \ psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \ - psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \ for (int v = -size; v <= size; v++) { \ for (int u = -size; u <= size; u++) { \ (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \ - /* Subtract (0,0) kernel from other kernels to preserve photometric scaling */ \ - if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \ - (TARGET)->kernel[v][u] += subValue * FUNC(subKernel->kernel[v][u]); \ - } \ + /* Photometric scaling is already built in to the precalculated kernel */ \ } \ } \ @@ -207,17 +210,17 @@ ) { - int numKernels = solution->n - 1; // Number of kernel basis functions + int numKernels = kernels->num; // Number of kernel basis functions assert(kernels->u->n == numKernels); assert(kernels->v->n == numKernels); - assert(kernels->xOrder->n == numKernels); - assert(kernels->yOrder->n == numKernels); + int spatialOrder = kernels->spatialOrder; // Maximum spatial polynomial order + assert(solution->n == numKernels * polyTerms(spatialOrder) + polyTerms(kernels->bgOrder)); // Ensure the subIndex for POIS kernels is what is expected assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS && kernels->type != PM_SUBTRACTION_KERNEL_SPAM && - kernels->type != PM_SUBTRACTION_KERNEL_FRIES) || - (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->xOrder->data.S32[kernels->subIndex] == 0 && - kernels->yOrder->data.S32[kernels->subIndex] == 0)); + kernels->type != PM_SUBTRACTION_KERNEL_FRIES && + kernels->type != PM_SUBTRACTION_KERNEL_GUNK && + kernels->type != PM_SUBTRACTION_KERNEL_RINGS) || + (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0)); int size = kernels->size; // Kernel half-size @@ -228,17 +231,19 @@ for (int i = 0; i < numKernels; i++) { - int xOrder = kernels->xOrder->data.S32[i]; // Polynomial order in x - int yOrder = kernels->yOrder->data.S32[i]; // Polynomial order in y - float value = weightFunc(polyValues->data.F64[yOrder][xOrder] * - solution->data.F64[i]); // Value to sum - float subValue = weightFunc(-solution->data.F64[i]); // Value to subtract (actually add, because of -) + double value = 0.0; // Value of kernel coefficient + for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { + value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index]; + } + } + switch (kernels->type) { case PM_SUBTRACTION_KERNEL_POIS: { int u = kernels->u->data.S32[i]; // Offset in x int v = kernels->v->data.S32[i]; // Offset in y - kernel->kernel[v][u] += value; + kernel->kernel[v][u] += weightFunc(value); if (kernels->spatialOrder > 0 && i != kernels->subIndex) { // The (0,0) element is subtracted from most kernels to preserve photometric scaling - kernel->kernel[0][0] += subValue; + kernel->kernel[0][0] += weightFunc(- value); } break; @@ -253,14 +258,14 @@ // Normalising sum of kernel component to unity - value /= weightFunc((uStop - uStart) * (vStop - vStart)); + float norm = 1.0 / weightFunc((uStop - uStart) * (vStop - vStart)); for (int v = vStart; v <= vStop; v++) { for (int u = uStart; u <= uStop; u++) { - kernel->kernel[v][u] += value; + kernel->kernel[v][u] += norm * weightFunc(value); + if (kernels->spatialOrder > 0 && i != kernels->subIndex) { + // The (0,0) element is subtracted from most kernels to preserve photometric scaling + kernel->kernel[0][0] += weightFunc(- value); + } } - } - if (kernels->spatialOrder > 0 && i != kernels->subIndex) { - /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ - kernel->kernel[0][0] += subValue; } break; @@ -270,38 +275,38 @@ // Using pre-calculated function psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values + // Iterating over the kernel for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - kernel->kernel[v][u] += value * weightFunc(preCalc->kernel[v][u]); + kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value); + // Photometric scaling is built into the preCalc kernel --- no subtraction! } } } else { // Using delta function + bool subtract = (kernels->spatialOrder > 0 && i != kernels->subIndex); // Subtract (0,0)? int u = kernels->u->data.S32[i]; // Offset in x int v = kernels->v->data.S32[i]; // Offset in y - kernel->kernel[v][u] += value; - } - // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling - if (kernels->spatialOrder > 0 && i != kernels->subIndex) { - kernel->kernel[0][0] += subValue; - } - break; - } - case PM_SUBTRACTION_KERNEL_ISIS: { - psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values - psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - kernel->kernel[v][u] += value * weightFunc(preCalc->kernel[v][u]); - // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling - if (kernels->spatialOrder > 0 && i != kernels->subIndex) { - kernel->kernel[v][u] += subValue * weightFunc(subKernel->kernel[v][u]); - } + kernel->kernel[v][u] += weightFunc(value); + if (subtract) { + // The (0,0) element is subtracted from most kernels to preserve photometric scaling + kernel->kernel[0][0] += weightFunc(- value); } } break; } + case PM_SUBTRACTION_KERNEL_ISIS: { + psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values + // Iterating over the kernel + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value); + // Photometric scaling is built into the preCalc kernel --- no subtraction! + } + } + break; + } case PM_SUBTRACTION_KERNEL_RINGS: { if (i == kernels->subIndex) { - kernel->kernel[0][0] += value; + kernel->kernel[0][0] += weightFunc(value); break; } @@ -311,14 +316,10 @@ psVector *poly = preCalc->data[2]; // Polynomial values int num = uCoords->n; // Number of pixels - value /= weightFunc(num); // Normalising sum of kernel component to unity + for (int j = 0; j < num; j++) { int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates - kernel->kernel[v][u] += value * weightFunc(poly->data.F32[j]); + kernel->kernel[v][u] += weightFunc(value * poly->data.F32[j]); } - - // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling - if (kernels->spatialOrder > 0) { - kernel->kernel[0][0] += subValue; - } + kernel->kernel[0][0] += weightFunc(- value * num); break; } @@ -336,12 +337,7 @@ int index, // Kernel basis function index int x, int y, // Pixel around which to convolve - const psImage *image, // Image to convolve - const psImage *polyValues // Spatial polynomial values + const psImage *image // Image to convolve ) { - int xOrder = kernels->xOrder->data.S32[index]; // Polynomial order in x - int yOrder = kernels->yOrder->data.S32[index]; // Polynomial order in y - double polyValue = polyValues->data.F64[yOrder][xOrder]; // Value of spatial polynomial - switch (kernels->type) { case PM_SUBTRACTION_KERNEL_POIS: { @@ -349,5 +345,5 @@ int u = kernels->u->data.S32[index]; // Offset in x int v = kernels->v->data.S32[index]; // Offset in y - double value = polyValue * image->data.F32[y + v][x + u]; // Value of convolution + float value = image->data.F32[y + v][x + u]; // Value of convolution if (kernels->spatialOrder > 0 && index != kernels->subIndex) { // The (0,0) element is subtracted from most kernels to preserve photometric scaling @@ -366,5 +362,5 @@ for (int v = vStart; v <= vStop; v++) { for (int u = uStart; u <= uStop; u++) { - sum += polyValue * image->data.F32[y + v][x + u]; + sum += image->data.F32[y + v][x + u]; } } @@ -377,5 +373,4 @@ } case PM_SUBTRACTION_KERNEL_GUNK: { - double value; // The value to return if (index < kernels->inner) { // Using pre-calculated function @@ -388,11 +383,10 @@ } } - value = polyValue * sum; - } else { - // Using delta function - int u = kernels->u->data.S32[index]; // Offset in x - int v = kernels->v->data.S32[index]; // Offset in y - value = polyValue * image->data.F32[y + v][x + u]; // Value of convolution - } + return sum; + } + // Using delta function + int u = kernels->u->data.S32[index]; // Offset in x + int v = kernels->v->data.S32[index]; // Offset in y + float value = image->data.F32[y + v][x + u]; // Value of convolution // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling if (kernels->spatialOrder > 0 && index != kernels->subIndex) { @@ -403,22 +397,17 @@ case PM_SUBTRACTION_KERNEL_ISIS: { psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel - psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract int size = kernels->size; // Kernel half-size double sum = 0.0; // Accumulated sum from convolution - double sub = 0.0; // Accumulated sum to subtract for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u]; - // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling - if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - sub += subKernel->kernel[v][u] * image->data.F32[y + v][x + u]; - } + // Photometric scaling is already built in to the precalculated kernel } } - return polyValue * sum - sub; + return sum; } case PM_SUBTRACTION_KERNEL_RINGS: { if (index == kernels->subIndex) { - return image->data.F32[0][0]; + return image->data.F32[y][x]; } psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data @@ -432,7 +421,6 @@ sum += image->data.F32[y + v][x + u] * poly->data.F32[j]; } - sum /= (double)num; // Normalising sum of kernel component to unity // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling - return polyValue * sum - image->data.F32[0][0]; + return sum - num * image->data.F32[y][x]; } default: @@ -592,6 +580,4 @@ PS_ASSERT_PTR_NON_NULL(kernels, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); if (weight) { PS_ASSERT_IMAGE_NON_NULL(weight, false); @@ -599,7 +585,4 @@ PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); } - if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); - } PS_ASSERT_INT_NONNEGATIVE(footprint, false); @@ -609,9 +592,13 @@ int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation - int numKernels = kernels->u->n; // Number of kernel basis functions - int numParams = numKernels + 1; // Total number of parameters to solve for: coefficient of each kernel - // basis function, and a constant background offset. - int bgIndex = numKernels; // Index in matrix for the background - psVector *convolutions = psVectorAlloc(numKernels, PS_TYPE_F64); // Convolutions for each kernel + int numKernels = kernels->num; // Number of kernel basis functions + int numSpatial = polyTerms(spatialOrder); // Number of spatial variations + int numBackground = polyTerms(kernels->bgOrder); // Number of background terms + + // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the + // number of coefficients for the spatial polynomial, and a constant background offset. + int numParams = numKernels * numSpatial + numBackground; + int bgIndex = numParams - numBackground; // Index in matrix for the background + psVector *convolutions = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); // Convolutions // We iterate over each stamp, allocate the matrix and vectors if @@ -637,8 +624,7 @@ psVectorInit(stampVector, 0.0); - // Pre-evaluated spatial polynomial values - psImage *polyValues = spatialPolyValues(spatialOrder, - 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols, - 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows); + float xNorm = 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols; // Normalised x coord + float yNorm = 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows; // Normalised y coord + psImage *polyValues = spatialPolyValues(spatialOrder, xNorm, yNorm); // Spatial polynomial terms for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { @@ -648,12 +634,19 @@ // Generate the convolutions - for (int i = 0; i < numKernels; i++) { - convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues); + for (int j = 0; j < numKernels; j++) { + double value = convolvePixel(kernels, j, x, y, reference); // Value from convolution + // Generate the pseudo-convolutions from the spatial polynomial terms + for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; + xOrder++, index += numKernels) { + convolutions->data.F64[index] = value * polyValues->data.F64[yOrder][xOrder]; + } + } } // Generate the least-squares matrix and vector // Upper diagonal only - for (int i = 0; i < numKernels; i++) { - for (int j = i; j < numKernels; j++) { + for (int i = 0; i < bgIndex; i++) { + for (int j = i; j < bgIndex; j++) { stampMatrix->data.F64[i][j] += convolutions->data.F64[i] * convolutions->data.F64[j] * invNoise2; @@ -674,5 +667,5 @@ // Fill in lower diagonal of symmetric matrix, while checking for bad values bool bad = false; // Are there bad values? - for (int i = 0; i < numKernels; i++) { + for (int i = 0; i < bgIndex; i++) { for (int j = 0; j < i; j++) { stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i]; @@ -694,5 +687,6 @@ if (bad) { stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - psTrace("psModules.imcombine", 3, "Rejecting stamp %d because of bad equation\n", i); + psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n", + i, stamp->x, stamp->y); } else { stamp->status = PM_SUBTRACTION_STAMP_USED; @@ -807,12 +801,8 @@ PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1); PS_ASSERT_PTR_NON_NULL(kernels, -1); - PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, -1); + PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + + polyTerms(kernels->bgOrder), -1); + PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, -1); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, -1); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, -1); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, -1); - if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - PS_ASSERT_ARRAY_NON_EMPTY(kernels->preCalc, -1); - PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, -1); - } // Measure deviations @@ -836,8 +826,7 @@ int xStamp = stamp->x, yStamp = stamp->y; // Coordinates of stamp - // Precalulate polynomial values - psImage *polyValues = spatialPolyValues(kernels->spatialOrder, - 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols, - 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows); + float xNorm = 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols; // Normalised x coord + float yNorm = 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows; // Normalised y coord + psImage *polyValues = spatialPolyValues(kernels->spatialOrder, xNorm, yNorm); // Polynomial terms #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS @@ -878,5 +867,6 @@ } deviations->data.F32[i] = sqrt(stats->sampleMean / 2.0); - psTrace("psModules.imcombine", 1, "Deviation for stamp %d: %f\n", i, deviations->data.F32[i]); + psTrace("psModules.imcombine", 1, "Deviation for stamp %d (%d,%d): %f\n", + i, stamp->x, stamp->y, deviations->data.F32[i]); totalSquareDev += PS_SQR(deviations->data.F32[i]); numStamps++; @@ -896,5 +886,5 @@ if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) { // Mask out the stamp in the image so you it's not found again - psTrace("psModules.imcombine", 3, "Rejecting stamp %d\n", i); + psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, stamp->x, stamp->y); numRejected++; for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { @@ -922,14 +912,10 @@ PS_ASSERT_VECTOR_NON_NULL(solution, NULL); PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + + polyTerms(kernels->bgOrder), NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false); + PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, NULL); - PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, NULL); - } // Precalulate polynomial values @@ -958,14 +944,10 @@ PS_ASSERT_VECTOR_NON_NULL(solution, NULL); PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + + polyTerms(kernels->bgOrder), NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false); + PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, NULL); - PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, NULL); - } psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return @@ -987,5 +969,6 @@ bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, const psImage *inImage, const psImage *inWeight, const psImage *subMask, - psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels) + psMaskType blank, const psRegion *region, const psVector *solution, + const pmSubtractionKernels *kernels) { PS_ASSERT_IMAGE_NON_NULL(inImage, false); @@ -1022,14 +1005,25 @@ PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false); PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false); + PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + + polyTerms(kernels->bgOrder), -1); + PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); - if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false); - PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, false); - } - - int size = kernels->size; // Half-size of kernel + if (region) { + if (psRegionIsNaN(*region)) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string); + psFree(string); + return false; + } + if (region->x0 < 0 || region->x1 > inImage->numCols || + region->y0 < 0 || region->y1 > inImage->numRows) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)", + string, inImage->numCols, inImage->numRows); + psFree(string); + return false; + } + } + float background = solution->data.F64[solution->n-1]; // The difference in background int numCols = inImage->numCols, numRows = inImage->numRows; // Image dimensions @@ -1057,12 +1051,23 @@ } + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1; // Full size of kernel psImage *polyValues = NULL; // Pre-calculated polynomial values - int fullSize = 2 * size + 1; // Full size of kernel psKernel *kernelImage = NULL; // Kernel for the image psKernel *kernelWeight = NULL; // Kernel for the weight map - for (int j = size; j < inImage->numRows - size; j += fullSize) { - for (int i = size; i < inImage->numCols - size; i += fullSize) { + // Get region for convolution: [xMin:xMax,yMin:yMax] + int xMin = size, xMax = numCols - size; + int yMin = size, yMax = numRows - size; + if (region) { + xMin = PS_MAX(region->x0, xMin); + xMax = PS_MIN(region->x1, xMax); + yMin = PS_MAX(region->y0, yMin); + yMax = PS_MIN(region->y1, yMax); + } + + for (int j = yMin; j < yMax; j += fullSize) { + for (int i = xMin; i < xMax; i += fullSize) { // Only generate polynomial values every kernel footprint, since we have already assumed @@ -1085,6 +1090,6 @@ #endif - for (int y = j; y < PS_MIN(j + fullSize, numRows - size); y++) { - for (int x = i; x < PS_MIN(i + fullSize, numCols - size); x++) { + for (int y = j; y < PS_MIN(j + fullSize, yMax); y++) { + for (int x = i; x < PS_MIN(i + fullSize, xMax); x++) { // Check and propagate the kernel footprint, if required if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] & @@ -1123,21 +1128,42 @@ psFree(polyValues); - // Mark the non-convolved part as blank - for (int y = size; y < inImage->numRows - size; y++) { + return true; +} + +bool pmSubtractionBorder(psImage *image, psImage *weight, psImage *mask, + const pmSubtractionKernels *kernels, psMaskType blank) +{ + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + if (mask) { + PS_ASSERT_IMAGE_NON_NULL(mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false); + } + if (weight) { + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); + } + + int numCols = image->numCols, numRows = image->numRows; // Image dimensions + + int size = kernels->size; // Half-size of kernel + for (int y = size; y < numRows - size; y++) { for (int x = 0; x < size; x++) { - markBlank(convImage, convMask, convWeight, x, y, blank); + markBlank(image, mask, weight, x, y, blank); } for (int x = numCols - size; x < numCols; x++) { - markBlank(convImage, convMask, convWeight, x, y, blank); + markBlank(image, mask, weight, x, y, blank); } } for (int y = 0; y < size; y++) { for (int x = 0; x < numCols; x++) { - markBlank(convImage, convMask, convWeight, x, y, blank); + markBlank(image, mask, weight, x, y, blank); } } for (int y = numRows - size; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - markBlank(convImage, convMask, convWeight, x, y, blank); + markBlank(image, mask, weight, x, y, blank); } } Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 14420) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 14455) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-20 01:53:52 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-09 20:25:52 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -83,8 +83,17 @@ const psImage *subMask, ///< Subtraction mask (or NULL) psMaskType blank, ///< Mask value for blank regions + const psRegion *region, ///< Region to convolve (or NULL) const psVector *solution, ///< The solution vector const pmSubtractionKernels *kernels ///< Kernel parameters ); +/// Mark the non-convolved part of the image as blank +bool pmSubtractionBorder(psImage *image,///< Image + psImage *weight, ///< Weight map (or NULL) + psImage *mask, ///< Mask (or NULL) + const pmSubtractionKernels *kernels, ///< Kernel parameters (for kernel size) + psMaskType blank ///< Mask value for blank regions + ); + /// @} #endif Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14420) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14455) @@ -19,6 +19,4 @@ psFree(kernels->uStop); psFree(kernels->vStop); - psFree(kernels->xOrder); - psFree(kernels->yOrder); psFree(kernels->preCalc); } @@ -50,4 +48,5 @@ kernels->type = type; + kernels->num = numBasisFunctions; kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); @@ -55,6 +54,4 @@ kernels->uStop = NULL; kernels->vStop = NULL; - kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); - kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->subIndex = 0; kernels->preCalc = NULL; @@ -62,4 +59,5 @@ kernels->inner = 0; kernels->spatialOrder = spatialOrder; + kernels->bgOrder = 0; return kernels; @@ -71,5 +69,5 @@ PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); - int num = PS_SQR(2 * size + 1) * (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions + int num = PS_SQR(2 * size + 1); // Number of basis functions pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, @@ -81,26 +79,15 @@ // Generate a set of kernels for each (u,v) for (int v = - size, index = 0; v <= size; v++) { - for (int u = - size; u <= size; u++) { - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { - for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->u->data.S32[index] = u; - kernels->v->data.S32[index] = v; - kernels->xOrder->data.S32[index] = xOrder; - kernels->yOrder->data.S32[index] = yOrder; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, - u, v, xOrder, yOrder); - } - } - } - } - - kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; + for (int u = - size; u <= size; u++, index++) { + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); + } + } + + kernels->subIndex = num / 2; assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->xOrder->data.S32[kernels->subIndex] == 0 && - kernels->yOrder->data.S32[kernels->subIndex] == 0); + kernels->v->data.S32[kernels->subIndex] == 0); return kernels; @@ -127,5 +114,4 @@ num += (gaussOrder + 1) * (gaussOrder + 2) / 2; } - num *= (spatialOrder + 1) * (spatialOrder + 2) / 2; pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, @@ -138,10 +124,12 @@ kernels->widths = psVectorAlloc(num, PS_TYPE_F32); kernels->preCalc = psArrayAlloc(num); + psKernel *subtract = NULL; // Kernel to subtract to maintain flux scaling // Set the kernel parameters for (int i = 0, index = 0; i < numGaussians; i++) { + float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian // Iterate over (u,v) order for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { - for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) { + for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { // Set the pre-calculated kernel @@ -150,37 +138,36 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) * + sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i])); } } - // Normalise sum of kernel component to unity - psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32)); - - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { - for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->widths->data.F32[index] = sigmas->data.F32[i]; - kernels->u->data.S32[index] = uOrder; - kernels->v->data.S32[index] = vOrder; - kernels->xOrder->data.S32[index] = xOrder; - kernels->yOrder->data.S32[index] = yOrder; - kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc); - - psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index, - sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder); + if (index == 0) { + subtract = preCalc; + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + preCalc->kernel[v][u] /= sum; + } + } + } else if (uOrder % 2 == 0 && vOrder % 2 == 0) { + // Normalise sum of kernel component to unity for even functions + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u]; + } } } - psFree(preCalc); // Drop reference + kernels->widths->data.F32[index] = sigmas->data.F32[i]; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + kernels->preCalc->data[index] = preCalc; + + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, + sigmas->data.F32[i], uOrder, vOrder); } } } - kernels->subIndex = 0; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->xOrder->data.S32[kernels->subIndex] == 0 && - kernels->yOrder->data.S32[kernels->subIndex] == 0); + kernels->subIndex = -1; if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { @@ -220,6 +207,5 @@ psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - int num = PS_SQR(2 * numTotal + 1) * - (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions + int num = PS_SQR(2 * numTotal + 1); // Number of basis functions psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); @@ -264,33 +250,23 @@ int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel - for (int j = - numTotal; j <= numTotal; j++) { + for (int j = - numTotal; j <= numTotal; j++, index++) { int v = locations->data.S32[numTotal + j]; // Location of pixel int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { - for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->u->data.S32[index] = u; - kernels->v->data.S32[index] = v; - kernels->uStop->data.S32[index] = uStop; - kernels->vStop->data.S32[index] = vStop; - kernels->xOrder->data.S32[index] = xOrder; - kernels->yOrder->data.S32[index] = yOrder; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, - u, uStop, v, vStop, xOrder, yOrder); - } - } - } - } - - kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, + u, uStop, v, vStop); + } + } + + kernels->subIndex = num / 2; assert(kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0 && kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0 && - kernels->xOrder->data.S32[kernels->subIndex] == 0 && - kernels->yOrder->data.S32[kernels->subIndex] == 0); + kernels->vStop->data.S32[kernels->subIndex] == 0); psFree(locations); @@ -327,6 +303,5 @@ psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - int num = PS_SQR(2 * numTotal + 1) * - (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions + int num = PS_SQR(2 * numTotal + 1); // Number of basis functions psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); @@ -369,33 +344,23 @@ int u = start->data.S32[numTotal + i]; // Location of pixel int uStop = stop->data.S32[numTotal + i]; // Width of pixel - for (int j = - numTotal; j <= numTotal; j++) { + for (int j = - numTotal; j <= numTotal; j++, index++) { int v = start->data.S32[numTotal + j]; // Location of pixel int vStop = stop->data.S32[numTotal + j]; // Width of pixel - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { - for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->u->data.S32[index] = u; - kernels->v->data.S32[index] = v; - kernels->uStop->data.S32[index] = uStop; - kernels->vStop->data.S32[index] = vStop; - kernels->xOrder->data.S32[index] = xOrder; - kernels->yOrder->data.S32[index] = yOrder; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, - u, uStop, v, vStop, xOrder, yOrder); - } - } - } - } - - kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, + u, uStop, v, vStop); + } + } + + kernels->subIndex = num / 2; assert(kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0 && kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0 && - kernels->xOrder->data.S32[kernels->subIndex] == 0 && - kernels->yOrder->data.S32[kernels->subIndex] == 0); + kernels->vStop->data.S32[kernels->subIndex] == 0); psFree(start); @@ -429,7 +394,6 @@ int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements - int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel - - int num = (numGaussianVars + numInner) * numSpatial; // Total number of basis functions + + int num = numGaussianVars + numInner; // Total number of basis functions pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK, @@ -440,15 +404,15 @@ psFree(params); - kernels->widths = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32); - kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial); - kernels->inner = numGaussianVars * numSpatial; + kernels->widths = psVectorAlloc(numGaussianVars, PS_TYPE_F32); + kernels->preCalc = psArrayAlloc(numGaussianVars); + kernels->inner = numGaussianVars; + psKernel *subtract = NULL; // Kernel to subtract to maintain flux scaling // Set the Gaussian kernel parameters for (int i = 0, index = 0; i < numGaussians; i++) { + float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian // Iterate over (u,v) order for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { - for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) { - - + for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { // Set the pre-calculated kernel psKernel *preCalc = psKernelAlloc(-size, size, -size, size); @@ -456,57 +420,48 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) * + sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i])); } } - // Normalise sum of kernel component to unity - psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32)); - - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { - for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->widths->data.F32[index] = sigmas->data.F32[i]; - kernels->u->data.S32[index] = uOrder; - kernels->v->data.S32[index] = vOrder; - kernels->xOrder->data.S32[index] = xOrder; - kernels->yOrder->data.S32[index] = yOrder; - kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc); - - psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index, - sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder); + if (index == 0) { + subtract = preCalc; + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + preCalc->kernel[v][u] /= sum; + } + } + } else if (uOrder % 2 == 0 && vOrder % 2 == 0) { + // Normalise sum of kernel component to unity for even functions + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u]; + } } } - psFree(preCalc); // Drop reference + kernels->widths->data.F32[index] = sigmas->data.F32[i]; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + kernels->preCalc->data[index] = preCalc; + + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, + sigmas->data.F32[i], uOrder, vOrder); } } } - // Generate a grid set of kernels for each (u,v) for (int v = - inner, index = kernels->inner; v <= inner; v++) { - for (int u = - inner; u <= inner; u++) { - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { - for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->u->data.S32[index] = u; - kernels->v->data.S32[index] = v; - kernels->xOrder->data.S32[index] = xOrder; - kernels->yOrder->data.S32[index] = yOrder; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, - u, v, xOrder, yOrder); - } - } - } - } - - kernels->subIndex = kernels->inner + (numInner - 1) * numSpatial / 2; + for (int u = - inner; u <= inner; u++, index++) { + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); + } + } + + kernels->subIndex = inner + num / 2; assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->xOrder->data.S32[kernels->subIndex] == 0 && - kernels->yOrder->data.S32[kernels->subIndex] == 0); + kernels->v->data.S32[kernels->subIndex] == 0); return kernels; @@ -543,7 +498,6 @@ int numRings = numOuter + numInner; // Number of rings (not including the central pixel) int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring - int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel - - int num = (numRings * numPoly + 1) * numSpatial; // Total number of basis functions + + int num = numRings * numPoly + 1; // Total number of basis functions pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, @@ -589,5 +543,5 @@ // Iterate over (u,v) order for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) { - for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++) { + for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++, index++) { psArray *data = psArrayAlloc(3); // Container for data @@ -604,4 +558,5 @@ } else { int j = 0; // Index for data + double norm = 0.0; // Normalisation for (int v = -size; v <= size; v++) { int v2 = PS_SQR(v); // Square of v @@ -616,5 +571,5 @@ uCoords->data.S32[j] = u; vCoords->data.S32[j] = v; - poly->data.F32[j] = uPoly * vPoly; + norm += poly->data.F32[j] = uPoly * vPoly; psVectorExtend(uCoords, RINGS_BUFFER, 1); @@ -629,23 +584,17 @@ } } + // Normalise kernel component to unit sum + psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); + } psTrace("psModules.imcombine", 8, "%ld pixels in ring\n", uCoords->n); - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { - for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->preCalc->data[index] = psMemIncrRefCounter(data); - kernels->u->data.S32[index] = uOrder; - kernels->v->data.S32[index] = vOrder; - kernels->xOrder->data.S32[index] = xOrder; - kernels->yOrder->data.S32[index] = yOrder; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d\n", index, - i, uOrder, vOrder, xOrder, yOrder); - } - } - psFree(data); + kernels->preCalc->data[index] = data; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, + i, uOrder, vOrder); } } @@ -653,6 +602,4 @@ kernels->subIndex = 0; - assert(kernels->xOrder->data.S32[kernels->subIndex] == 0 && - kernels->yOrder->data.S32[kernels->subIndex] == 0); return kernels; Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14420) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14455) @@ -17,8 +17,8 @@ typedef struct { pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels + long num; ///< Number of kernel components (not including the spatial ones) psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) psVector *widths; ///< Gaussian widths (ISIS) or ring radius (RINGS) psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - psVector *xOrder, *yOrder; ///< Spatial Polynomial order (for all) int subIndex; ///< Index of kernel to be subtracted (to maintain flux conservation) psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) @@ -26,4 +26,5 @@ int inner; ///< The size of an inner region int spatialOrder; ///< The spatial order of the kernels + int bgOrder; ///< The order for the background fitting } pmSubtractionKernels; Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14420) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14455) @@ -42,5 +42,5 @@ psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, - float threshold, float spacing) + const psRegion *region, float threshold, float spacing) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -52,12 +52,36 @@ } PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); + if (region) { + if (psRegionIsNaN(*region)) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string); + psFree(string); + return false; + } + if (region->x0 < 0 || region->x1 > image->numCols || + region->y0 < 0 || region->y1 > image->numRows) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)", + string, image->numCols, image->numRows); + psFree(string); + return false; + } + } - // Size of image - int numRows = image->numRows; - int numCols = image->numCols; + int numRows = image->numRows, numCols = image->numCols; // Size of image + + // Get region in which to find stamps: [xMin:xMax,yMin:yMax] + int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows; + if (region) { + xMin = PS_MAX(region->x0, xMin); + xMax = PS_MIN(region->x1, xMax); + yMin = PS_MAX(region->y0, yMin); + yMax = PS_MIN(region->y1, yMax); + } + int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest // Number of stamps - int xNumStamps = numCols / spacing + 1; - int yNumStamps = numRows / spacing + 1; + int xNumStamps = xSize / spacing + 1; + int yNumStamps = ySize / spacing + 1; int numFound = 0; // Number of stamps found @@ -92,12 +116,12 @@ // Bounds of region to search for stamp - int yMin = j * numRows / (yNumStamps + 1); - int yMax = (j + 1) * numRows / (yNumStamps + 1) - 1; - int xMin = i * numCols / (xNumStamps + 1); - int xMax = (i + 1) * numCols / (xNumStamps + 1) - 1; - assert(yMax < image->numRows && xMax < image->numCols && yMin >= 0 && xMin >= 0); + int yStart = yMin + j * ySize / (yNumStamps + 1); + int yStop = yMin + (j + 1) * ySize / (yNumStamps + 1) - 1; + int xStart = xMin + i * xSize / (xNumStamps + 1); + int xStop = xMin + (i + 1) * xSize / (xNumStamps + 1) - 1; + assert(yStop < numRows && xStop < numCols && yStart >= 0 && xStart >= 0); - for (int y = yMin; y <= yMax ; y++) { - for (int x = xMin; x <= xMax ; x++) { + for (int y = yStart; y <= yStop ; y++) { + for (int x = xStart; x <= xStop ; x++) { if ((!subMask || !(subMask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT | Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14420) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14455) @@ -26,4 +26,5 @@ const psImage *image, ///< Image for which to find stamps const psImage *mask, ///< Mask + const psRegion *region, ///< Region in which to search for stamps float threshold, ///< Threshold for stamps in the image float spacing ///< Rough spacing for stamps