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); } }