Index: trunk/psModules/src/pmImageSubtract.c =================================================================== --- trunk/psModules/src/pmImageSubtract.c (revision 4168) +++ trunk/psModules/src/pmImageSubtract.c (revision 4185) @@ -4,7 +4,8 @@ * * @author GLG, MHPCC + * @author Paul Price, IfA * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-09 00:59:53 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-06-09 06:18:19 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -52,5 +53,4 @@ psFree(stamp->matrix); psFree(stamp->vector); - psFree(stamp); } @@ -392,9 +392,10 @@ // Determine if this pixel is larger than the max, and unmasked. - if ((image->data.F32[y][x] > max) && - !(mask->data.U8[y][x] & mask->data.U8[y][x])) { - max = image->data.F32[y][x]; - bestx = x; - besty = y; + if (image->data.F32[y][x] > max) { + if ((mask == NULL) || !((mask->data.U8[y][x]) & maskVal)) { + max = image->data.F32[y][x]; + bestx = x; + besty = y; + } } } @@ -405,5 +406,5 @@ // Otherwise, mark the stamp as PM_STAMP_NONE // - if (image->data.F32[besty][bestx] > threshold) { + if (image->data.F32[besty][bestx] >= threshold) { stamp->x = bestx; stamp->y = besty; @@ -451,7 +452,58 @@ /******************************************************************************* +GeneralKernelConvolve(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. + ******************************************************************************/ +psF32 GeneralKernelConvolve(const psImage *input, + const psSubtractionKernels *kernels, + psS32 kernelID, + psS32 col, + psS32 row) +{ + psS32 spatialOrder = kernels->p_spatialOrder; + psS32 kernelSize = kernels->p_size; + psS32 xOrder = (psS32) kernels->xOrder->data.F32[kernelID]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[kernelID]; + psF32 numColsHalf = 0.5 * (psF32) input->numCols; + psF32 numRowsHalf = 0.5 * (psF32) input->numRows; + 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]; + + psImage *preCalc = (psImage *) kernels->preCalc->data[kernelID]; + psF32 conv = 0.0; + for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { + for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { + //printf("HERE: (%d, %d)\n", yy, xx); + //printf("HERE: (%d, %d)\n", yy+row, xx+col); + //printf("HERE: (%d, %d)\n", yy-kernelSize, xx-kernelSize); + //printf("KERNEL SIZE is %d\n", kernelSize); + conv += input->data.F32[yy+row][xx+col] * + preCalc->data.F32[yy+kernelSize][xx+kernelSize] * + polyVal; + } + } + psFree(polyValues); + + return(conv); +} + + + + +/******************************************************************************* XXX: How is the spatial order factor calculated? Is it simply a x^iy^j power evaluated at the (x, y) center for the stamp? Why bother with the 2-D polynomial in evaluating it? + +XXX: Should we assert that the footprint is equal to the kernel size, or does +they have nothing in common. ******************************************************************************/ bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, @@ -494,6 +546,7 @@ for (psS32 s = 0; s < stamps->n; s++) { pmStamp *stamp = (pmStamp *) stamps->data[s]; - + psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d\n", s); if (stamp->status == PM_STAMP_RECALC) { + psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: status is PM_STAMP_RECALC.\n", s); psImage *stampMatrix = stamp->matrix; psVector *stampVector = stamp->vector; @@ -501,4 +554,5 @@ if (stampMatrix == NULL) { stampMatrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); + stamp->matrix = stampMatrix; } else { PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, false); @@ -509,4 +563,5 @@ if (stampVector == NULL) { stampVector = psVectorAlloc(numSolveParams, PS_TYPE_F64); + stamp->vector = stampVector; } else { PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, false); @@ -514,4 +569,5 @@ } PS_VECTOR_SET_F64(stampVector, 0.0); + psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: allocate matrix and vector.\n", s); // @@ -524,4 +580,5 @@ ((psF64) (stamp->y - numHalfRows)) / ((psF64) numHalfRows)); + psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s); // // Iterate over all pixels surrounding this stamp. @@ -534,4 +591,5 @@ // 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) { // @@ -547,4 +605,8 @@ // First convolution. This will set the value for the stampVector. // + // + // XXX: verify the [y-v2][x-u2] subscript. This generated errors in + // testing, depending on kernel size and footprint. + // psF32 conv1 = polyValues->data.F64[j1][i1] * reference->data.F32[y - v1][x - u1]; @@ -561,13 +623,21 @@ // for (psS32 k2 = k1; k2 < numKernels; k2++) { - psS32 u2 = kernels->u->data.F32[k2]; // Offset in x - psS32 v2 = kernels->v->data.F32[k2]; // Offset in y - psS32 i2 = kernels->xOrder->data.F32[k2]; // Polynomial order in x - psS32 j2 = kernels->yOrder->data.F32[k2]; // Polynomial order in y - + psS32 u2 = (psS32) kernels->u->data.F32[k2]; // Offset in x + psS32 v2 = (psS32) kernels->v->data.F32[k2]; // Offset in y + psS32 i2 = (psS32) kernels->xOrder->data.F32[k2]; // Polynomial order in x + psS32 j2 = (psS32) kernels->yOrder->data.F32[k2]; // Polynomial order in y + // + // XXX: verify the [y-v2][x-u2] subscript. This generated errors in + // testing, depending on kernel size and footprint. + // + //printf("footprint is %d\n", footprint); + //printf("Stamp (%d, %d).\n", stamp->y, stamp->x); + //printf("HERE (y, x) is (%d, %d). (v2, u2) is (%d, %d).\n", y, x, v2, u2); + //printf("HERE 00 (%d %d) (%d %d)\n", j2, i2, y-v2, x-u2); // // Second convolution // - psF32 conv2 = polyValues->data.F64[j2][i2] * reference->data.F32[y - v2][x - u2]; + psF32 conv2 = polyValues->data.F64[j2][i2] * + reference->data.F32[y-v2][x-u2]; // @@ -600,21 +670,19 @@ } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - printf("XXX: put some warning message here (ISIS kernels not implemented).\n"); - return(false); + // printf("XXX: put some warning message here (ISIS kernels not implemented).\n"); + // return(false); // XXX: HEY: code this - /* for (psS32 k1 = 0; k1 < numKernels; k1++) { - psF32 conv1 = GeneralKernelConvolve(reference, kernels, k1); + psF32 conv1 = GeneralKernelConvolve(reference, kernels, k1, x, y); for (psS32 k2 = k1; k2 < numKernels; k2++) { - psF32 conv2 = GeneralKernelConvolve(reference, kernels, k2); + psF32 conv2 = GeneralKernelConvolve(reference, kernels, k2, x, y); stampMatrix->data.F64[k1][k2] += conv1 * conv2 * invNoise2; } - vector->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 { printf("XXX: put some warning message here (bad kernel->type).\n"); @@ -662,4 +730,5 @@ /******************************************************************************* +XXX: Assert correct vector matrix sizes. ******************************************************************************/ psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL @@ -668,4 +737,6 @@ { PS_ASSERT_PTR_NON_NULL(stamps, NULL); + PS_ASSERT_IMAGE_NON_NULL(((pmStamp *) stamps->data[0])->matrix, NULL); + PS_ASSERT_VECTOR_NON_NULL(((pmStamp *) stamps->data[0])->vector, NULL); psS32 size = ((pmStamp *) stamps->data[0])->vector->n; @@ -804,15 +875,9 @@ // Iterate over the kernel basis functions for (psS32 k = 0; k < nBF; 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 - - psImage *preCalc = (psImage *) kernels->preCalc->data[nBF]; + + psImage *preCalc = (psImage *) kernels->preCalc->data[k]; for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { @@ -820,5 +885,5 @@ conv += solution->data.F64[k] * input->data.F32[yy+row][xx+col] * - preCalc->data.F32[yy-kernelSize][xx-kernelSize] * + preCalc->data.F32[yy+kernelSize][xx+kernelSize] * polyVal; } @@ -854,5 +919,13 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, 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); @@ -943,5 +1016,13 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, 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); @@ -1026,5 +1107,14 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); + if (kernels->preCalc != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); + } 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(false); + } + } + psS32 nBF = kernels->u->n; PS_ASSERT_VECTOR_SIZE(solution, nBF+1, false); @@ -1091,5 +1181,13 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, 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); + } + } PS_ASSERT_INT_EQUAL(1+kernels->u->n, solution->n, NULL); @@ -1097,69 +1195,58 @@ psS32 spatialOrder = kernels->p_spatialOrder; psS32 kernelSize = kernels->p_size; - psS32 xCenter; // The pixel location for the center of the kernel in out img. - psS32 yCenter; // The pixel location for the center of the kernel in out img. - psS32 numRows; - psS32 numCols; - psS32 imgSize; if (out != NULL) { - xCenter = (psS32) ( ((x+1.0)/2.0) * (psF32) (out->numCols)); - yCenter = (psS32) ( ((y+1.0)/2.0) * (psF32) (out->numRows)); - if ( ((xCenter - kernelSize) < 0) || - ((xCenter + kernelSize) > out->numCols) || - ((yCenter - kernelSize) < 0) || - ((yCenter + kernelSize) > out->numRows)) { - + if ((out->numCols < (1+2*kernelSize)) || (out->numRows < (1+2*kernelSize))) { printf("XXX: generate WARNING: out image is not large enough.\n"); return(out); } } else { - // - // We calculate the minimize size image that can contain the kernel - // centered at the specified (x, y) position. - // - numRows = (psS32) ceilf( ((psF32) kernelSize) / ((y + 1.0)/2.0)); - numRows = PS_MAX(numRows, (((psF32) kernelSize) / (1.0 - ((y + 1.0)/2.0)))); - numCols = (psS32) ceilf( ((psF32) kernelSize) / ((x + 1.0)/2.0)); - numCols = PS_MAX(numRows, (((psF32) kernelSize) / (1.0 - ((x + 1.0)/2.0)))); - imgSize = PS_MAX(numRows, numCols); - - out = psImageAlloc(imgSize, imgSize, PS_TYPE_F32); - xCenter = (psS32) ( ((x+1.0)/2.0) * (psF32) (out->numCols)); - yCenter = (psS32) ( ((y+1.0)/2.0) * (psF32) (out->numRows)); + out = psImageAlloc(1+2*kernelSize, 1+2*kernelSize, PS_TYPE_F32); } PS_IMAGE_SET_F32(out, 0.0); + + // + // Generate the spatial-order polynomial. The [i][j]-th element of + // the psImage polyValues will hold (x^i * y^j) for the stamp. + // + psImage *polyValues = GenSpatialOrder(spatialOrder, x, y); // XXX: switch (i, j) so they correspond to (x, y). if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { for (psS32 k = 0 ; k < nBF ; k++) { - // XXX: verify these loop bounds with the size of the preCalc images. + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; + + // XXX: Verify that this is correct. for (psS32 i = -kernelSize ; i <= kernelSize ; i++) { for (psS32 j = -kernelSize ; j <= kernelSize ; j++) { - psImage *preCalc = (psImage *) kernels->preCalc->data[nBF]; - out->data.F32[yCenter + i][xCenter + j]+= solution->data.F64[nBF] * - preCalc->data.F32[i+kernelSize][j+kernelSize]; + psImage *preCalc = (psImage *) kernels->preCalc->data[k]; + out->data.F32[i+kernelSize][j+kernelSize]+= + solution->data.F64[k] * + preCalc->data.F32[i+kernelSize][j+kernelSize] * + polyVal; } } } } else if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { - // - // Generate the spatial-order polynomial. The [i][j]-th element of - // the psImage polyValues will hold (x^i * y^j) for the stamp. - // - psImage *polyValues = GenSpatialOrder(spatialOrder, x, y); - for (psS32 k = 0 ; k < nBF ; k++) { - // XXX: Why don't we have compilation warnings on type here? - psS32 u = kernels->u->data.F32[nBF]; - psS32 v = kernels->v->data.F32[nBF]; - psS32 xOrder = kernels->xOrder->data.F32[nBF]; - psS32 yOrder = kernels->yOrder->data.F32[nBF]; + // XXX: Why don't we have compilation warnings on type here (if + // we remove the (psS32) cast)? + psS32 u = (psS32) kernels->u->data.F32[k]; + psS32 v = (psS32) kernels->v->data.F32[k]; + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; // XXX: Verify that this is correct. - out->data.F32[yCenter - (psS32) v][xCenter - (psS32) u]+= solution->data.F64[nBF] * - polyValues->data.F64[yOrder][xOrder]; - } - psFree(polyValues); - } + + out->data.F32[kernelSize - v][kernelSize - u]+= + solution->data.F64[k] * + polyValues->data.F64[yOrder][xOrder] * + polyVal; + polyVal = polyVal; + } + } + psFree(polyValues); return(out);