Index: /trunk/psModules/src/pmImageSubtract.c =================================================================== --- /trunk/psModules/src/pmImageSubtract.c (revision 4416) +++ /trunk/psModules/src/pmImageSubtract.c (revision 4417) @@ -1,39 +1,40 @@ /** @file ImageSubtract.c -* -* This file will ... -* -* @author Paul Price, IfA (original prototype) -* @author GLG, MHPCC -* -* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-06-20 23:21:05 $ -* -* Copyright 2004 Maui High Performance Computing Center, University of Hawaii -* -* XXX: sync with iFa on this: -* The (x, y) (row, col) issue is becoming a problem. In this file, and I -* think, the rest of psLib and psModules, the following conventions are used: -* -* 1) x will correspond to the column, and y will correspond to the row. -* 2) When used in function prototypes, the column (and hence x) appears -* first. -* 3) When used to index 2-D arrays, obviously, the row (and hence, y) -* appears first. (2 and 3 are the source of confusion). -* 4) When (u, v) are used in certain structures. -* u corresponds to x -* v corresponds to y -* 5) When element (a, b) is casually referred to in comments, or -* documentation it is unclear where a is the row, or the column. -* 6) A convention on loop index variables (i, j) would be convenient. -* Currently, sometimes i corresponds to the column (x), -* usually it corresponds to the row (y). -* -* XXX: The following variables are used an interpreted this way: -* kernelSize: Means that the actual kernel is a square (1 + 2 * kernelSize) per side. -* border: When accessing an image, a swath of pixels this wide is ignored. -* footprint: When accessing a stample, a square of pixels, footprint pixels per side, -* are looked at. We must ensure that (footprint+kernelSize) pixels exist -* around the center. -*/ + * + * This file will contain code which creates a set of kernel basis + * functions, solves for their solution, and applies them to an image. + * + * @author Paul Price, IfA (original prototype) + * @author GLG, MHPCC + * + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-06-28 23:23:55 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * XXX: sync with IfA on this: + * The (x, y) (row, col) issue is becoming a problem. In this file, and I + * think, the rest of psLib and psModules, the following conventions are used: + * + * 1) x will correspond to the column, and y will correspond to the row. + * 2) When used in function prototypes, the column (and hence x) appears + * first. + * 3) When used to index 2-D arrays, obviously, the row (and hence, y) + * appears first. (2 and 3 are the source of confusion). + * 4) When (u, v) are used in certain structures. + * u corresponds to x + * v corresponds to y + * 5) When element (a, b) is casually referred to in comments, or + * documentation it is unclear where a is the row, or the column. + * 6) A convention on loop index variables (i, j) would be convenient. + * Currently, sometimes i corresponds to the column (x), + * usually it corresponds to the row (y). + * + * XXX: The following variables are used an interpreted this way: + * kernelSize: Means that the actual kernel is a square (1 + 2 * kernelSize) per side. + * border: When accessing an image, a swath of pixels this wide is ignored. + * footprint: When accessing a stample, a square of pixels, footprint pixels per side, + * are looked at. We must ensure that (footprint+kernelSize) pixels exist + * around the center. + */ #include @@ -43,16 +44,7 @@ #include "pmImageSubtract.h" -// XXX: Put this is psConstants.h -#define PS_ASSERT_INT_LESS_THAN(VAR1, VAR2, RVAL) \ -if (!(VAR1 < VAR2)) { \ - psError(PS_ERR_UNKNOWN, true, \ - "Error: %s is not less than %s (%d, %d)", #VAR1, #VAR2, VAR1, VAR2); \ - return(RVAL); \ -} - /******************************************************************************* Private alloc/free functions. - -XXX: It's not clear if the SubtractionKernels alloc/free functions are doable. +XXX: It's not clear if the SubtractionKernels alloc/free functions are feasable. ******************************************************************************/ void p_pmStampFree(pmStamp *stamp) @@ -90,5 +82,5 @@ } -psSubtractionKernels *p_pmSubtractionKernelsAlloc(psS32 numBasisFunctions, +psSubtractionKernels *p_pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type) { @@ -106,4 +98,6 @@ int spatialOrder) { + psTrace("ImageSubtract.pmSubtractionKernelsAllocPOIS", 1, + "Calling pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder); PS_ASSERT_INT_POSITIVE(size, NULL); PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); @@ -152,5 +146,5 @@ // - // Iterate over (u,v). Generate a set of kernels foreach (u, v). + // Iterate over (u,v). Generate a set of kernels for each (u, v). // for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { @@ -176,4 +170,6 @@ } + psTrace("ImageSubtract.pmSubtractionKernelsAllocPOIS", 1, + "Exiting pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder); return(tmp); } @@ -181,8 +177,4 @@ /******************************************************************************* XXX: Get the types correct (u, v, xOrder, yOrder). - -XXX: For the various spatial orders, verify that the preCalc data does not -vary with the (x^i * y^j) since that comes from the image coordinates when the -kernel is convolved with it. XXX: Code review this. @@ -193,12 +185,15 @@ int spatialOrder) { + PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + psTrace("ImageSubtract.pmSubtractionKernelsAllocISIS", 1, + "Calling pmSubtractionKernelsAllocISIS(%d, %d, %d, %d)\n", + sigmas->n, orders->n, size, spatialOrder); PS_ASSERT_INT_POSITIVE(size, NULL); PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); - PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(orders, NULL); PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); // - // Calculate the number of basis functions (nBF), and others. + // Calculate the number of basis functions (nBF). // psS32 numSigmas = sigmas->n; @@ -283,5 +278,7 @@ (2.0 * PS_SQR(sigmas->data.F32[sigPtr])); currPreCalc->data.F32[v+size][u+size] = - exp(exponent) * pow(uScaled, orderXTerm) * pow(vScaled, orderYTerm); + exp(-exponent) * + pow(uScaled, orderXTerm) * + pow(vScaled, orderYTerm); } } @@ -298,4 +295,7 @@ } + psTrace("ImageSubtract.pmSubtractionKernelsAllocISIS", 1, + "Exiting pmSubtractionKernelsAllocISIS(%d, %d, %d, %d)\n", + sigmas->n, orders->n, size, spatialOrder); return(tmp); } @@ -310,5 +310,5 @@ XXX: Do we need to care about case where yNum/xNum does not evenly divide the -nnumber of rows/columns in the mage? +nnumber of rows/columns in the image? ******************************************************************************/ psArray *pmSubtractionFindStamps(psArray *stamps, ///< Output stamps, or NULL @@ -322,4 +322,7 @@ ) { + psTrace("ImageSubtract.pmSubtractionFindStamps", 1, + "Calling pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n", + maskVal, threshold, xNum, yNum, border); PS_ASSERT_IMAGE_NON_NULL(image, NULL); PS_ASSERT_IMAGE_NON_EMPTY(image, NULL); @@ -360,4 +363,5 @@ // // XXX: Must handle cases where image size is not an even multiple of xNum or yNum + // they are currently ignored. // psS32 num = 0; @@ -440,4 +444,7 @@ } } + psTrace("ImageSubtract.pmSubtractionFindStamps", 1, + "Exiting pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n", + maskVal, threshold, xNum, yNum, border); return(stamps); } @@ -455,4 +462,7 @@ psF32 y) { + psTrace("ImageSubtract.GenSpatialOrder", 5, + "Calling GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y); + psImage *polyValues = psImageAlloc(spatialOrder+1, spatialOrder+1, PS_TYPE_F64); @@ -468,4 +478,7 @@ } + psTrace("ImageSubtract.GenSpatialOrder", 5, + "Exiting GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y); + return(polyValues); } @@ -473,19 +486,17 @@ /******************************************************************************* -GeneralKernelConvolve(input, kernels, kernelID, col, row): This routine +IsisKernelConvolve(input, kernels, kernelID, col, row): This routine convolves a single kernel basis function with a pixel in an image. - -XXX: merge this code with the other convolution code in this file. - -XXX: Must we ensure that the pixels being accessed are inside the image? Is -there any garantee that the kernelSize, footprint stuff will be set correctly -to ensure no seg faults. - ******************************************************************************/ -static psF32 GeneralKernelConvolve(const psImage *input, - const psSubtractionKernels *kernels, - psS32 kernelID, - psS32 col, - psS32 row) -{ + ******************************************************************************/ +static psF32 IsisKernelConvolve(const psImage *input, + const psSubtractionKernels *kernels, + psS32 kernelID, + psS32 col, + psS32 row) +{ + + + psTrace("ImageSubtract.IsisKernelConvolve", 5, + "Calling IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row); psS32 spatialOrder = kernels->p_spatialOrder; psS32 kernelSize = kernels->p_size; @@ -496,8 +507,16 @@ psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position + psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); - psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; + + psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; psImage *preCalc = (psImage *) kernels->preCalc->data[kernelID]; + + // XXX: Are the following asserts really necessary? + PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(row-kernelSize, 0, NAN); + PS_ASSERT_INT_LESS_THAN(row+kernelSize, input->numRows, NAN); + PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(col-kernelSize, 0, NAN); + PS_ASSERT_INT_LESS_THAN(col+kernelSize, input->numCols, NAN); psF32 conv = 0.0; for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { @@ -510,13 +529,198 @@ psFree(polyValues); + psTrace("ImageSubtract.IsisKernelConvolve", 5, + "Exiting IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row); return(conv); } - - - /******************************************************************************* -XXX: Should we assert that the footprint is equal to the kernel size, or does -they have nothing in common. +ConvolvePixelPois(input, mask, badStampMaskVal, solution, kernels, col, row): + +This routine takes a single pixel in the psImage input and convolves it with +the set of kernel basis functions and their appropriate weights in solution. +It returns the value of the convolved pixel. + +XXX: Static structure for polyValues? + ******************************************************************************/ +static psF32 ConvolvePixelPois(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels, + psS32 col, + psS32 row) +{ + psTrace("ImageSubtract.ConvolvePixelPois", 5, + "Calling ConvolvePixelPois(%d, %d)\n", col, row); + psS32 nBF = kernels->u->n; + psF32 numColsHalf = 0.5 * (psF32) input->numCols; + psF32 numRowsHalf = 0.5 * (psF32) input->numRows; + psF32 background = solution->data.F64[solution->n-1]; + psS32 spatialOrder = kernels->p_spatialOrder; + psF32 conv = background; // Initial convolved value + + if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { + psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position + psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position + psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); + + // Iterate over the kernel basis functions + for (psS32 k = 0; k < nBF; k++) { + psS32 u = (psS32) kernels->u->data.F32[k]; + psS32 v = (psS32) kernels->v->data.F32[k]; + + // XXX: What's the story with this? + #if 0 + + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; + #else + + psF32 polyVal = 1.0; + #endif + + // XXX: Why this? + if (k == 0) { + conv += solution->data.F64[k] * input->data.F32[row - v][col - u] * polyVal; + } else { + conv += solution->data.F64[k] * + (input->data.F32[row - v][col - u] * polyVal - input->data.F32[row][col]); + } + } + psFree(polyValues); + } + + psTrace("ImageSubtract.ConvolvePixelPois", 5, + "Exiting ConvolvePixelPois(%d, %d)\n", col, row); + return(conv); +} + + + +/******************************************************************************* +ConvolvePixelIsis(input, mask, badStampMaskVal, solution, kernels, col, row): + +This routine takes a single pixel in the psImage input and convolves it with +the set of kernel basis functions and their appropriate weights in solution. +It returns the value of the convolved pixel. + +XXX: Static structure for polyValues? + ******************************************************************************/ +static psF32 ConvolvePixelIsis(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels, + psS32 col, + psS32 row) +{ + psTrace("ImageSubtract.ConvolvePixelIsis", 5, + "Calling ConvolvePixelIsis(%d, %d)\n", col, row); + psF32 background = solution->data.F64[solution->n-1]; + psF32 conv = background; // Initial convolved value + + if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { + // Iterate over the kernel basis functions + for (psS32 k = 0; k < kernels->u->n; k++) { + conv += IsisKernelConvolve(input, kernels, k, col, row); + } + } + + psTrace("ImageSubtract.ConvolvePixelIsis", 5, + "Exiting ConvolvePixelIsis(%d, %d)\n", col, row); + return(conv); +} + +/******************************************************************************* +ConvolveImage(input, mask, badStampMaskVal, solution, kernels): convolves an +arbitrary image with either an ISIS or POIS set of kernel basis functions. + ******************************************************************************/ +static psImage *ConvolveImage(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels) +{ + psTrace("ImageSubtract.ConvolveImage", 5, "Calling ConvolveImage()\n"); + PS_ASSERT_IMAGE_NON_NULL(input, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(input, NULL); + PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, NULL); + if (mask != NULL) { + PS_ASSERT_IMAGES_SIZE_EQUAL(input, mask, NULL); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); + } + PS_ASSERT_VECTOR_NON_NULL(solution, NULL); + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); + PS_ASSERT_PTR_NON_NULL(kernels, NULL); + 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->preCalc != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); + } else { + if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); + return(NULL); + } + } + psS32 nBF = kernels->u->n; + PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); + + psS32 numCols = input->numCols; + psS32 numRows = input->numRows; + psS32 kernelSize = kernels->p_size; + + psImage *convolved = psImageAlloc(numCols, numRows, PS_TYPE_F32); + + for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { + for (psS32 x = kernelSize; x < numCols - kernelSize; x++) { + if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + convolved->data.F32[y][x] = ConvolvePixelPois(input, mask, badStampMaskVal, + solution, kernels, x, y); + } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + convolved->data.F32[y][x] = ConvolvePixelIsis(input, mask, badStampMaskVal, + solution, kernels, x, y); + } else { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel type. Returning NULL\n"); + return(NULL); + } + } + } + + // + // Pad the rest of the convolved image with 0.0 + // + for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { + for (psS32 x = 0; x < kernelSize; x++) { + convolved->data.F32[y][x] = 0.0; + } + for (psS32 x = numCols - kernelSize; x < numCols; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + for (psS32 y = 0; y < kernelSize; y++) { + for (psS32 x = 0; x < numCols; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + for (psS32 y = numRows - kernelSize; y < numRows; y++) { + for (psS32 x = 0; x < numCols; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + + psTrace("ImageSubtract.ConvolveImage", 5, "Exiting ConvolveImage()\n"); + return convolved; +} + + + + + +/******************************************************************************* +XXX: We should assert that the (footprint, kernelSize, imageSize) stuff +ensures that all data is accessed in bounds? ******************************************************************************/ bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, @@ -527,4 +731,6 @@ ) { + psTrace("ImageSubtract.pmSubtractionCalculateEquation", 1, + "Calling pmSubtractionCalculateEquation()\n"); PS_ASSERT_PTR_NON_NULL(stamps, false); PS_ASSERT_IMAGE_NON_NULL(reference, false); @@ -623,17 +829,15 @@ psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s); - // - // Iterate over all pixels surrounding this stamp. - // - for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { - for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { - - psF32 invNoise2 = 1.0/reference->data.F32[y][x]; // The inverse of the noise, squared. - - // XXX: Think about a better decomposition of the computation - // based on whether the kernels are ISIS or POIS. - - psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): pixel (%d, %d).\n", y, x); - if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + // + // Iterate over all pixels surrounding this stamp. + // + for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { + for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { + psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): pixel (%d, %d).\n", y, x); + + // The inverse of the noise, squared. + psF32 invNoise2 = 1.0/reference->data.F32[y][x]; + // // Iterate over the first convolution */ @@ -647,5 +851,4 @@ // // First convolution. This will set the value for the stampVector. - // // // XXX: verify the [y-v2][x-u2] subscript. This generated errors in @@ -705,24 +908,50 @@ stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; - - } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + } + } + } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + // + // Iterate over all pixels surrounding this stamp. + // + // XXX: Why isn't there a polyValues term here? + // + + for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { + for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { + psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): pixel (%d, %d).\n", y, x); + psF32 invNoise2 = 1.0/reference->data.F32[y][x]; // The inverse of the noise, squared. + for (psS32 k1 = 0; k1 < numKernels; k1++) { - psF32 conv1 = GeneralKernelConvolve(reference, kernels, k1, x, y); + psF32 conv1 = IsisKernelConvolve(reference, kernels, k1, x, y); for (psS32 k2 = k1; k2 < numKernels; k2++) { //printf("(k1, k2) is (%d, %d)\n", k1, k2); - psF32 conv2 = GeneralKernelConvolve(reference, kernels, k2, x, y); - stampMatrix->data.F64[k1][k2] += conv1 * conv2 * invNoise2; + psF32 conv2 = IsisKernelConvolve(reference, kernels, k2, x, y); + stampMatrix->data.F64[k1][k2]+= conv1 * conv2 * invNoise2; } - stampVector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2; + stampVector->data.F64[k1]+= input->data.F32[y][x] * conv1 * invNoise2; stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2; } - } else { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel->type.\n"); - return(false); - } // if-else on kernel type. - } // y-loop for all pixels. - } // x-loop for all pixels. + stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; + stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; + } + } + } else { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel->type.\n"); + return(false); + } psFree(polyValues); + + // XXX: Generate psTrace() + if (0) { + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + if (stamp->status == PM_STAMP_RECALC) { + psVector *stampVector = stamp->vector; + printf("COOL: stamp %d vector:\n", s); + PS_VECTOR_PRINT_F64(stampVector); + } + } + } // @@ -756,4 +985,6 @@ } } + psTrace("ImageSubtract.pmSubtractionCalculateEquation", 1, + "Exiting pmSubtractionCalculateEquation()\n"); return(true); } @@ -763,5 +994,4 @@ /******************************************************************************* -XXX: Assert correct vector matrix sizes. ******************************************************************************/ psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL @@ -769,4 +999,6 @@ ) { + psTrace("ImageSubtract.pmSubtractionSolveEquation", 1, + "Calling pmSubtractionSolveEquation()\n"); PS_ASSERT_PTR_NON_NULL(stamps, NULL); psS32 size = -1; @@ -802,4 +1034,6 @@ // // XXX: Test these functions with size=-1. This caused seg faults during test. + // This should be done in the psImage.c and psVector.c test files. It + // should never occur here. // psImage *sumMatrix = psImageAlloc(size, size, PS_TYPE_F64); @@ -829,9 +1063,39 @@ } } - + psVector *permutation = NULL; // XXX: Check output from these routines. - psVector *permutation = NULL; + + // XXX: psTrace() + if (1) { + PS_IMAGE_PRINT_F64(sumMatrix); + } + psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); + if (luMatrix == NULL) { + psError(PS_ERR_UNKNOWN, true, "Failed to LU-Decompose the matrix.\n"); + psFree(sumMatrix); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return(NULL); + } + // XXX: psTrace() + if (1) { + PS_IMAGE_PRINT_F64(luMatrix); + } + solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation); + // XXX: psTrace() + if (1) { + PS_VECTOR_PRINT_F64(solution); + } + if (solution == NULL) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the matrix.\n"); + psFree(sumMatrix); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return(NULL); + } psFree(sumMatrix); @@ -840,201 +1104,7 @@ psFree(permutation); + psTrace("ImageSubtract.pmSubtractionSolveEquation", 1, + "Exiting pmSubtractionSolveEquation()\n"); return(solution); -} - - -/******************************************************************************* -ConvolvePixelPois(input, mask, badStampMaskVal, solution, kernels, col, row): - -This routine takes a single pixel in the psImage input and convolves it with -the set of kernel basis functions and their appropriate weights in solution. -It returns the value of the convolved pixel. - -XXX: Static structure for polyValues? - ******************************************************************************/ -static psF32 ConvolvePixelPois(const psImage *input, - const psImage *mask, - psU32 badStampMaskVal, - const psVector *solution, - const psSubtractionKernels *kernels, - psS32 col, - psS32 row) -{ - psS32 nBF = kernels->u->n; - psF32 numColsHalf = 0.5 * (psF32) input->numCols; - psF32 numRowsHalf = 0.5 * (psF32) input->numRows; - psF32 background = solution->data.F64[solution->n-1]; - psS32 spatialOrder = kernels->p_spatialOrder; - psF32 conv = background; // Initial convolved value - - if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { - psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position - psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position - psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); - - // Iterate over the kernel basis functions - for (psS32 k = 0; k < nBF; k++) { - psS32 u = (psS32) kernels->u->data.F32[k]; - psS32 v = (psS32) kernels->v->data.F32[k]; - - // XXX: What's the story with this? - #if 0 - - psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; - psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; - psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; - #else - - psF32 polyVal = 1.0; - #endif - - // XXX: Why this? - if (k == 0) { - conv += solution->data.F64[k] * input->data.F32[row - v][col - u] * polyVal; - } else { - conv += solution->data.F64[k] * - (input->data.F32[row - v][col - u] * polyVal - input->data.F32[row][col]); - } - } - psFree(polyValues); - } - - return(conv); -} - - - -/******************************************************************************* -ConvolvePixelIsis(input, mask, badStampMaskVal, solution, kernels, col, row): - -This routine takes a single pixel in the psImage input and convolves it with -the set of kernel basis functions and their appropriate weights in solution. -It returns the value of the convolved pixel. - -XXX: Static structure for polyValues? - ******************************************************************************/ -static psF32 ConvolvePixelIsis(const psImage *input, - const psImage *mask, - psU32 badStampMaskVal, - const psVector *solution, - const psSubtractionKernels *kernels, - psS32 col, - psS32 row) -{ - psS32 nBF = kernels->u->n; - psF32 numColsHalf = 0.5 * (psF32) input->numCols; - psF32 numRowsHalf = 0.5 * (psF32) input->numRows; - psF32 background = solution->data.F64[solution->n-1]; - psS32 spatialOrder = kernels->p_spatialOrder; - psS32 kernelSize = kernels->p_size; - psF32 conv = background; // Initial convolved value - - if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { - psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position - psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position - psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); - - // Iterate over the kernel basis functions - for (psS32 k = 0; k < nBF; k++) { - psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; - psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; - psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; - - psImage *preCalc = (psImage *) kernels->preCalc->data[k]; - for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { - for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { - // XXX: Should I do something special for the k==0 kernel? - conv += solution->data.F64[k] * - input->data.F32[yy+row][xx+col] * - preCalc->data.F32[yy+kernelSize][xx+kernelSize] * - polyVal; - } - } - } - psFree(polyValues); - } - - return(conv); -} - -/******************************************************************************* -ConvolveImage(input, mask, badStampMaskVal, solution, kernels): convolves an -arbitrary image with either an ISIS or POIS set of kernel basis functions. - ******************************************************************************/ -static psImage *ConvolveImage(const psImage *input, - const psImage *mask, - psU32 badStampMaskVal, - const psVector *solution, - const psSubtractionKernels *kernels) -{ - PS_ASSERT_IMAGE_NON_NULL(input, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(input, NULL); - PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, NULL); - if (mask != NULL) { - PS_ASSERT_IMAGES_SIZE_EQUAL(input, mask, NULL); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(solution, NULL); - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); - 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->preCalc != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); - } else { - if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); - return(NULL); - } - } - psS32 nBF = kernels->u->n; - PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); - - psS32 numCols = input->numCols; - psS32 numRows = input->numRows; - psS32 kernelSize = kernels->p_size; - - psImage *convolved = psImageAlloc(numCols, numRows, PS_TYPE_F32); - - for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { - for (psS32 x = kernelSize; x < numCols - kernelSize; x++) { - if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { - convolved->data.F32[y][x] = ConvolvePixelPois(input, mask, badStampMaskVal, - solution, kernels, x, y); - } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - convolved->data.F32[y][x] = ConvolvePixelIsis(input, mask, badStampMaskVal, - solution, kernels, x, y); - } else { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel type. Returning NULL\n"); - return(NULL); - } - } - } - - // - // Pad the rest of the convolved image with 0.0 - // - for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { - for (psS32 x = 0; x < kernelSize; x++) { - convolved->data.F32[y][x] = 0.0; - } - for (psS32 x = numCols - kernelSize; x < numCols; x++) { - convolved->data.F32[y][x] = 0.0; - } - } - for (psS32 y = 0; y < kernelSize; y++) { - for (psS32 x = 0; x < numCols; x++) { - convolved->data.F32[y][x] = 0.0; - } - } - for (psS32 y = numRows - kernelSize; y < numRows; y++) { - for (psS32 x = 0; x < numCols; x++) { - convolved->data.F32[y][x] = 0.0; - } - } - - return convolved; } @@ -1052,4 +1122,6 @@ const psVector *solution) { + psTrace("ImageSubtract.CalculateDeviations", 5, + "Calling CalculateDeviations()\n"); PS_ASSERT_PTR_NON_NULL(stamps, NULL); if (deviations != NULL) { @@ -1115,10 +1187,5 @@ psImage *subStampTrim = psImageSubset((psImage *) subStamp, myReg); psImage *maskStampTrim = psImageSubset((psImage *) maskStamp, myReg); - // Copy image to workaround bug 305 - // XXX: Modify. This bug is now fixed. - // psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32); - // psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8); - // (void)psImageStats(stats, tempImage, tempMask, badStampMaskVal); - (void)psImageStats(stats, subStampTrim, maskStampTrim, badStampMaskVal); + psImageStats(stats, subStampTrim, maskStampTrim, badStampMaskVal); deviations->data.F32[s] = stats->sampleMean * (psF32)footprint * (psF32)footprint * 4.0; @@ -1136,4 +1203,6 @@ psFree(subStamp); + psTrace("ImageSubtract.CalculateDeviations", 5, + "Exiting CalculateDeviations()\n"); return deviations; } @@ -1152,4 +1221,6 @@ ) { + psTrace("ImageSubtract.pmSubtractionRejectStamps", 5, + "Calling pmSubtractionRejectStamps()\n"); PS_ASSERT_PTR_NON_NULL(stamps, false); PS_ASSERT_IMAGE_NON_NULL(refImage, false); @@ -1225,4 +1296,6 @@ psFree(deviations); + psTrace("ImageSubtract.pmSubtractionRejectStamps", 5, + "Exiting pmSubtractionRejectStamps()\n"); return(true); } @@ -1237,4 +1310,6 @@ ) { + psTrace("ImageSubtract.pmSubtractionKernelImage", 5, + "Calling pmSubtractionKernelImage()\n"); PS_ASSERT_VECTOR_NON_NULL(solution, NULL); PS_ASSERT_PTR_NON_NULL(kernels, NULL); @@ -1280,5 +1355,5 @@ psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; - psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; + psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; // XXX: Verify that this is correct. @@ -1304,9 +1379,11 @@ out->data.F32[kernelSize - v][kernelSize - u]+= - solution->data.F64[k] * polyValues->data.F32[yOrder][xOrder]; + solution->data.F64[k] * polyValues->data.F64[yOrder][xOrder]; } } psFree(polyValues); + psTrace("ImageSubtract.pmSubtractionKernelImage", 5, + "Exiting pmSubtractionKernelImage()\n"); return(out); } Index: /trunk/psModules/src/pmImageSubtract.h =================================================================== --- /trunk/psModules/src/pmImageSubtract.h (revision 4416) +++ /trunk/psModules/src/pmImageSubtract.h (revision 4417) @@ -1,10 +1,12 @@ /** @file ImageSubtract.h * - * This file will ... + * This file will contain code which creates a set of kernel basis + * functions, solves for their solution, and applies them to an image. * + * @author Paul Price, IfA (original prototype) * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-20 23:21:05 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-06-28 23:23:55 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,5 +28,5 @@ typedef enum { PM_SUBTRACTION_KERNEL_POIS, ///< POIS kernel --- delta functions - PM_SUBTRACTION_KERNEL_ISIS ///< ISIS kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_ISIS ///< ISIS kernel --- gaussians modified by polynomials } pmSubtractionKernelsType; @@ -32,12 +34,12 @@ { pmSubtractionKernelsType type; ///< Type ofKernels --- allowing the use of multiple kernels - psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) - psVector *sigma; ///< Width of Gaussian (for ISIS) - 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 (to + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) + psVector *sigma; ///< Width of Gaussian (for ISIS) + 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 (to ///< accelerate ISIS; don't use for POIS) - psS32 p_size; ///< The halfsize of the kernel - psS32 p_spatialOrder; ///< The spatial order of the kernels + psS32 p_size; ///< The halfsize of the kernel + psS32 p_spatialOrder; ///< The spatial order of the kernels } psSubtractionKernels; @@ -52,9 +54,9 @@ typedef enum { - PM_STAMP_INIT, ///< Initial state - PM_STAMP_USED, ///< Use this stamp - PM_STAMP_REJECTED, ///< This stamp has been rejected - PM_STAMP_RECALC, ///< Having been reset, this stamp is to be recalculated - PM_STAMP_NONE ///< No stamp in this region + PM_STAMP_INIT, ///< Initial state + PM_STAMP_USED, ///< Use this stamp + PM_STAMP_REJECTED, ///< This stamp has been rejected + PM_STAMP_RECALC, ///< Having been reset, this stamp is to be recalculated + PM_STAMP_NONE ///< No stamp in this region } pmStampStatus;