Changeset 4322 for trunk/psModules/src/pmImageSubtract.c

Timestamp:

Jun 20, 2005, 1:21:05 PM (21 years ago)

Author:

gusciora

Message:

The code is in better shape than the release version. CalculateEquation
generates a resonably looking image kernel (to me) for POIS. However, for
ISIS kernels, we are generating NANs.

File:

• trunk/psModules/src/pmImageSubtract.c (modified) (20 diffs)

trunk/psModules/src/pmImageSubtract.c

@@ -1 +1 @@
 /** @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<stdio.h>
@@ -60 +66 @@
     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 +319 @@
                                  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 +359 @@
     // 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 +385 @@
                 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 +503 @@
     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 +540 @@
         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 +622 @@
 
             psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s);
+
             //
             // Iterate over all pixels surrounding this stamp.
@@ -572 +628 @@
             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 +674 @@
                                 // 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 +707 @@
 
                     } 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 +752 @@
             }
             stamp->status = PM_STAMP_USED;
+        } else {
+            // Stamp is ignored since it's not PM_STAMP_RECALC
         }
     }
@@ -721 +770 @@
 {
     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 +798 @@
     }
 
+    //
+    // 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 +808 @@
     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 +830 @@
     }
 
+    // XXX: Check output from these routines.
     psVector *permutation = NULL;
     psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
@@ -930 +1008 @@
                                             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 +1100 @@
         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 +1261 @@
     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 +1301 @@
             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];
         }
     }