Index: trunk/psModules/src/pmImageSubtract.c =================================================================== --- trunk/psModules/src/pmImageSubtract.c (revision 4223) +++ trunk/psModules/src/pmImageSubtract.c (revision 4322) @@ -1,33 +1,39 @@ /** @file ImageSubtract.c - * - * This file will ... - * - * @author Paul Price, IfA (original prototype) - * @author GLG, MHPCC - * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-13 20:09:53 $ - * - * 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). - * - */ +* +* 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. +*/ #include @@ -60,5 +66,7 @@ pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp)); stamp->x = 0; + stamp->p_xSize = 0; stamp->y = 0; + stamp->p_ySize = 0; stamp->matrix = NULL; stamp->vector = NULL; @@ -311,5 +319,5 @@ psS32 xNum, ///< Number of stamps in x psS32 yNum, ///< Number of stamps in y - psS32 border ///< Border around image to ignore (should be size of kernel) + psS32 border ///< Border around image to ignore (should be size of kernel or larger) ) { @@ -351,9 +359,10 @@ // Iterate over the image sections // + // XXX: Must handle cases where image size is not an even multiple of xNum or yNum + // psS32 num = 0; for (psS32 j = 0; j < yNum; j++) { for (psS32 i = 0; i < xNum; i++) { pmStamp *stamp = (pmStamp *) stamps->data[num]; - // // Only find a new stamp if we need to @@ -376,32 +385,54 @@ 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[y][x] > max) { - if ((mask == NULL) || !((mask->data.U8[y][x]) & maskVal)) { - max = image->data.F32[y][x]; - bestx = x; - besty = y; + psS32 yMin = border + j * (numCols - 2.0 * border) / numY; + psS32 yMax = (border + (j + 1) * (numCols - 2.0 * border) / numY) - 1; + psS32 xMin = border + i * (numRows - 2.0 * border) / numX; + psS32 xMax = (border + (i + 1) * (numRows - 2.0 * border) / numX) - 1; + + if ((yMax >= image->numRows) || + (xMax >= image->numCols) || + (yMin < 0) || + (xMin < 0)) { + // XXX: We skip this stamp since its borders extends beyond the image. + // XXX: This is here mainly as a safeguard. We need to redefine the above + // min/max pixels calculation to ensure that all stamps are legitimate. + + stamp->x = -1; + stamp->y = -1; + stamp->status = PM_STAMP_NONE; + } else { + stamp->p_xSize = 1 + (xMax - xMin); + stamp->p_ySize = 1 + (yMax - yMin); + stamp->p_xMin = xMin; + stamp->p_xMax = xMax; + stamp->p_yMin = yMin; + stamp->p_yMax = yMax; + + for (psS32 y = yMin; y <= yMax ; y++) { + for (psS32 x = xMin; x <= xMax ; x++) { + // Determine if this pixel is larger than the max, and unmasked. + 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; + } } } } - } - - // - // If the max pixel is larger than the threshold, we keep this stamp. - // Otherwise, mark the stamp as PM_STAMP_NONE - // - if (image->data.F32[besty][bestx] >= threshold) { - stamp->x = bestx; - stamp->y = besty; - stamp->status = PM_STAMP_RECALC; - } else { - stamp->x = bestx; - stamp->y = besty; - stamp->status = PM_STAMP_NONE; + + // + // If the max pixel is larger than the threshold, we keep this stamp. + // Otherwise, mark the stamp as PM_STAMP_NONE + // + if (image->data.F32[besty][bestx] >= threshold) { + stamp->x = bestx; + stamp->y = besty; + stamp->status = PM_STAMP_RECALC; + } else { + stamp->x = bestx; + stamp->y = besty; + stamp->status = PM_STAMP_NONE; + } } } @@ -472,8 +503,4 @@ 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] * @@ -513,5 +540,33 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); } + psS32 kernelSize = kernels->p_size; PS_ASSERT_INT_NONNEGATIVE(footprint, false); + // + // For each legitimate stamp, ensure that the footprint is small enough to perform + // the full calculation. + // + // XXX: Verify with IfA that this is a reasonable action. + // + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + if (stamp->status == PM_STAMP_RECALC) { + // XXX: trace message + // printf("stamp %d (x, y) is (%d, %d). Footprint is %d. kernelSize is %d.\n", s, stamp->x, stamp->y, footprint, kernelSize); + if (((stamp->y - (footprint + kernelSize)) < 0) || + ((stamp->x - (footprint + kernelSize)) < 0) || + ((stamp->y + footprint + kernelSize) >= input->numRows) || + ((stamp->x + footprint + kernelSize) >= input->numCols)) { + stamp->status = PM_STAMP_NONE; + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: stamp %d will be ignored. It exceeds image size: access columns (%d to %d) and rows (%d to %d)\n", + s, + stamp->x - (footprint + kernelSize), + (stamp->x + footprint + kernelSize) - 1, + stamp->y - (footprint + kernelSize), + (stamp->y + footprint + kernelSize) - 1); + } + } + } + psS32 numHalfRows = reference->numRows; psS32 numHalfCols = reference->numCols; @@ -567,4 +622,5 @@ psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s); + // // Iterate over all pixels surrounding this stamp. @@ -572,4 +628,5 @@ 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. @@ -617,14 +674,8 @@ // 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]; - // // Assuming that the first kernel component is 0 order in x and y, and 0 offset @@ -656,15 +707,11 @@ } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - // 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, 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; - } stampVector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2; @@ -705,4 +752,6 @@ } stamp->status = PM_STAMP_USED; + } else { + // Stamp is ignored since it's not PM_STAMP_RECALC } } @@ -721,7 +770,24 @@ { 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; + psS32 size = -1; + psS32 s = 0; + + // + // Determine the size of the stamp vectors and matrix. + // We iterate until we find the first acceptable stamp. + // + while ((size == -1) && (s < stamps->n)) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + PS_ASSERT_PTR_NON_NULL(stamp, NULL); + if (stamp->status == PM_STAMP_USED) { + size = ((pmStamp *) stamps->data[s])->vector->n; + PS_ASSERT_INT_POSITIVE(size, NULL); + } + s++; + } + if (size == -1) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: no acceptable stamps. Returning NULL\n"); + return(NULL); + } if (solution != NULL) { @@ -732,4 +798,9 @@ } + // + // Create the solution matrix and vector. + // + // XXX: Test these functions with size=-1. This caused seg faults during test. + // psImage *sumMatrix = psImageAlloc(size, size, PS_TYPE_F64); psVector *sumVector = psVectorAlloc(size, PS_TYPE_F64); @@ -737,15 +808,21 @@ PS_IMAGE_SET_F64(sumMatrix, 0.0); + // + // Verify that all stamps have similar sizes. + // Compute the sum matrix and vector. + // for (psS32 s = 0; s < stamps->n; s++) { pmStamp *stamp = (pmStamp *) stamps->data[s]; - psImage *stampMatrix = stamp->matrix; - psVector *stampVector = stamp->vector; - - PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_SIZE(stampVector, size, NULL); - PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, NULL); - PS_ASSERT_IMAGE_SIZE(stampMatrix, size, size, NULL); if (stamp->status == PM_STAMP_USED) { + PS_ASSERT_INT_EQUAL(((pmStamp *) stamps->data[s])->vector->n, size, NULL); + + psImage *stampMatrix = stamp->matrix; + psVector *stampVector = stamp->vector; + PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_SIZE(stampVector, size, NULL); + PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, NULL); + PS_ASSERT_IMAGE_SIZE(stampMatrix, size, size, NULL); + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stampMatrix); (void)psBinaryOp(sumVector, sumVector, "+", stampVector); @@ -753,4 +830,5 @@ } + // XXX: Check output from these routines. psVector *permutation = NULL; psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); @@ -930,5 +1008,5 @@ solution, kernels, x, y); } else { - printf("XXX: Generate WARNING: unknown kernel type\n"); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel type. Returning NULL\n"); return(NULL); } @@ -1022,4 +1100,5 @@ psS32 y = stamp->y; // Stamp y coord if (stamp->status == PM_STAMP_USED) { + psRegion myReg = psRegionSet(x - xSize, x + xSize, y - ySize, y + ySize); psImage *refStamp = psImageSubset((psImage *) refImage, myReg); @@ -1182,5 +1261,5 @@ if (out != NULL) { if ((out->numCols < (1+2*kernelSize)) || (out->numRows < (1+2*kernelSize))) { - printf("XXX: generate WARNING: out image is not large enough.\n"); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: out image is not large enough.\n"); return(out); } @@ -1222,12 +1301,8 @@ 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[kernelSize - v][kernelSize - u]+= - solution->data.F64[k] * - polyValues->data.F64[yOrder][xOrder] * - polyVal; - polyVal = polyVal; + solution->data.F64[k] * polyValues->data.F32[yOrder][xOrder]; } }