Index: /trunk/psModules/src/pmImageSubtract.c =================================================================== --- /trunk/psModules/src/pmImageSubtract.c (revision 4035) +++ /trunk/psModules/src/pmImageSubtract.c (revision 4036) @@ -5,8 +5,26 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-05-25 23:08:33 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-25 23:59:08 $ * * 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). * */ @@ -18,15 +36,5 @@ #include "pmImageSubtract.h" -// XXX: Put this is psConstants.h -#define PS_ASSERT_VECTOR_SIZE(VEC, SIZE, RVAL) \ -if (VEC->n != SIZE) { \ - psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ - "psVector %s has size %d, should be %d." \ - #VEC, VEC->n, SIZE); \ - return(RVAL); \ -} - /******************************************************************************* -XXX: Store the size and spatialOrder in a private member of the psSubtractionKernels struct. ******************************************************************************/ @@ -67,4 +75,6 @@ // for convenience. // + // XXX: replace (i, j) with (xOrder, yOrder) for clarity. + // psS32 ptr = 0; for (psS32 order = 0; order <= spatialOrder; order++) { @@ -89,4 +99,6 @@ // Iterate over spatial order. This loop creates the terms for // x^i * y^j such that (i+j) <= spatialOrder. + // + // XXX: replace (i, j) with (xOrder, yOrder) for clarity. // for (psS32 order = 0; order <= spatialOrder; order++) { @@ -132,5 +144,5 @@ psS32 numSigmas = sigmas->n; psS32 nBF = 0; - for (psS32 i = 0 ; i < numSigmas ; i++) { + for (psS32 s = 0 ; s < numSigmas ; s++) { nBF+= ((orders->n + 1) * (orders->n + 2) / 2); } @@ -146,6 +158,6 @@ // XXX: There should be a macro or a psLib vectorDuplicate(). tmp->sigma = psVectorAlloc(numSigmas, PS_TYPE_F32); - for (psS32 i = 0 ; i < numSigmas ; i++) { - tmp->sigma->data.F32[i] = sigmas->data.F32[i]; + for (psS32 s = 0 ; s < numSigmas ; s++) { + tmp->sigma->data.F32[s] = sigmas->data.F32[s]; } tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); @@ -165,7 +177,8 @@ // This loop creates // - // - // + // XXX: Put the v loop before the u loop. + // XXX: Put the orderYTerm loop before the orderXTerm loop. // XXX: must loop on the sigmas. + // psS32 sigPtr = 0; // XXX: Review ptr. @@ -175,5 +188,5 @@ for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { psS32 orderYTerm = order - orderXTerm; - // psS32 imgPtr = orderYTerm + (orderXTerm * spatialOrder); + // psS32 imgPtr = orderYTerm + (orderXTerm * spatialOrder); psImage *currKernel = (psImage *) psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); @@ -194,7 +207,9 @@ } - // // sigPtr loops through each sigma value in the *sigmas argument. + // XXX: Put the v loop before the u loop. + // XXX: Put the orderYTerm loop before the orderXTerm loop. + // psS32 ptr = 0; for (psS32 sigPtr = 0 ; sigPtr < numSigmas ; sigPtr++) { @@ -218,4 +233,6 @@ for (psS32 order = 0; order <= spatialOrder; order++) { for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { + psS32 orderYTerm = order - orderXTerm; + // XXX: Should simply add a pointer to an existing image and increase // the memory reference counter for equivalent spatial orders. @@ -227,5 +244,4 @@ } - psS32 orderYTerm = order - orderXTerm; tmp->u->data.F32[ptr] = xOrderP; tmp->v->data.F32[ptr] = yOrderP; @@ -331,6 +347,6 @@ // psS32 num = 0; - for (psS32 i = 0; i < xNum; i++) { - for (psS32 j = 0; j < yNum; j++) { + for (psS32 j = 0; j < yNum; j++) { + for (psS32 i = 0; i < xNum; i++) { pmStamp *stamp = (pmStamp *) stamps->data[num]; @@ -351,16 +367,17 @@ // border of pixels from the image edge. // - // XXX: must set these - psS32 numX = 0; - psS32 numY = 0; - for (psS32 x = border + i * (numRows - 2.0 * border) / numX; - x < border + (i + 1) * (numRows - 2.0 * border) / numX; x++) { - for (psS32 y = border + j * (numCols - 2.0 * border) / numY; - y < border + (j + 1) * (numCols - 2.0 * border) / numY; y++) { + // XXX: verify this (numX, numY), then get rid of it. + // + psS32 numX = xNum; + psS32 numY = yNum; + for (psS32 y = border + j * (numCols - 2.0 * border) / numY; + y < border + (j + 1) * (numCols - 2.0 * border) / numY; y++) { + for (psS32 x = border + i * (numRows - 2.0 * border) / numX; + x < border + (i + 1) * (numRows - 2.0 * border) / numX; x++) { // Determine if this pixel is larger than the max, and unmasked. - if ((image->data.F32[x][y] > max) && - !(mask->data.U8[x][y] & mask->data.U8[x][y])) { - max = image->data.F32[x][y]; + 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; @@ -373,5 +390,5 @@ // Otherwise, mark the stamp as PM_STAMP_NONE // - if (image->data.F32[bestx][besty] > threshold) { + if (image->data.F32[besty][bestx] > threshold) { stamp->x = bestx; stamp->y = besty; @@ -394,4 +411,6 @@ XXX: Modify loop so that terms higher than spatialOrder are not computed. + +XXX: Modify this so that [i][j] location is calculated as (x^j * y^i)? ******************************************************************************/ static psImage *GenSpatialOrder(psS32 spatialOrder, @@ -490,14 +509,16 @@ // psImage *polyValues = GenSpatialOrder(spatialOrder, - ((psF64) (stamp->y - numHalfCols)) / ((psF64) numHalfCols), - ((psF64) (stamp->x - numHalfRows)) / ((psF64) numHalfRows)); + ((psF64) (stamp->x - numHalfCols)) / ((psF64) numHalfCols), + ((psF64) (stamp->y - numHalfRows)) / ((psF64) numHalfRows)); // // Iterate over all pixels surrounding this stamp. // + // XXX: I modified the y,x loop order, and an [x][y] reference (1.0/reference->data.F32[y][x]). + // Verify it. // - for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { - for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { - psF32 invNoise2 = 1.0/reference->data.F32[x][y]; // The inverse of the noise, squared. + 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 @@ -711,6 +732,4 @@ psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position - - // XXX: Ensure you have the (x, y) order correct. psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); @@ -722,5 +741,5 @@ // XXX: What's the story with this? #if 0 - // XXX: x/y order correct? + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; @@ -777,6 +796,4 @@ psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position - - // XXX: Ensure you have the (x, y) order correct. psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); @@ -785,5 +802,4 @@ // XXX: What's the story with this? #if 0 - // XXX: x/y order correct? psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; @@ -798,6 +814,4 @@ for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { // XXX: Should I do something special for the k==0 kernel? - // The POIS code does. - conv += solution->data.F64[k] * input->data.F32[yy+row][xx+col] * @@ -943,6 +957,5 @@ psS32 y = stamp->y; // Stamp y coord if (stamp->status == PM_STAMP_USED) { - // XXX: Verify proper x,y order - psRegion myReg = psRegionSet(y - ySize, y + ySize, x - xSize, x + xSize); + psRegion myReg = psRegionSet(x - xSize, x + xSize, y - ySize, y + ySize); psImage *refStamp = psImageSubset((psImage *) refImage, myReg); psImage *inStamp = psImageSubset((psImage *) inImage, myReg); @@ -954,5 +967,4 @@ (void)psBinaryOp(subStamp, subStamp, "/", inStamp); (void)psBinaryOp(subStamp, subStamp, "*", subStamp); - // XXX: Verify proper x,y order myReg = psRegionSet(kernelSize, kernelSize + 2 * footprint, kernelSize, kernelSize + 2 * footprint); @@ -1069,4 +1081,5 @@ /******************************************************************************* +XXX: There are definite (x, y) (row, col) problems here. ******************************************************************************/ psImage *pmSubtractionKernelImage(psImage *out, @@ -1107,6 +1120,6 @@ // centered at the specified (x, y) position. // - numRows = (psS32) ceilf( ((psF32) kernelSize) / ((x + 1.0)/2.0)); - numRows = PS_MAX(numRows, (((psF32) kernelSize) / (1.0 - ((x + 1.0)/2.0)))); + 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)))); @@ -1119,4 +1132,5 @@ PS_IMAGE_SET_F32(out, 0.0); + // XXX: switch (i, j) so they correspond to (x, y). if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { for (psS32 k = 0 ; k < nBF ; k++) { @@ -1143,6 +1157,4 @@ psS32 xOrder = kernels->xOrder->data.F32[nBF]; psS32 yOrder = kernels->yOrder->data.F32[nBF]; - // XXX: Verify that (u, v) correspond to (x, y) and that we have the row,col stuff correct. - // XXX: Verify the same for xOrder,yOrder. // XXX: Verify that this is correct. out->data.F32[yCenter - (psS32) v][xCenter - (psS32) u]+= solution->data.F64[nBF] * @@ -1154,3 +1166,2 @@ return(out); } -// This code is