Changeset 4322

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.

Location:

trunk/psModules

Files:

• src/pmImageSubtract.c (modified) (20 diffs)
• src/pmImageSubtract.h (modified) (2 diffs)
• test/tst_pmImageSubtract.c (modified) (11 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];
         }
     }

trunk/psModules/src/pmImageSubtract.h

@@ -5 +5 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-05-25 23:00:50 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-06-20 23:21:05 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -62 +62 @@
 {
     int x, y;      ///< Position
+    int p_xSize;
+    int p_ySize;
+    int p_xMin;
+    int p_xMax;
+    int p_yMin;
+    int p_yMax;
     psImage *matrix;     ///< Associated matrix
     psVector *vector;     ///< Assoicated vector

trunk/psModules/test/tst_pmImageSubtract.c

@@ -10 +10 @@
  *  data.  More work need to be done to verify the results.
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-06-13 20:12:18 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-06-20 23:21:05 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -33 +33 @@
 int main(int argc, char* argv[])
 {
+    psLogSetFormat("HLNM");
     return !runTestSuite(stderr, "Test Point Driver", tests, argc, argv);
 }
@@ -495 +496 @@
 }
 
-#define TST03_THRESHOLD 3.0
-#define TST03_MASK_VAL 1
+
+#define TST03_THRESHOLD  3.0
+#define TST03_MASK_VAL  1
 #define TST03_KERNEL_SIZE 2
 #define TST03_SPATIAL_ORDER 3
 #define TST03_ORDER_LENGTH 2
 #define TST03_SIGMA_LENGTH 3
-#define TST03_FOOTPRINT 5
+#define TST03_PSF_MAX  10.0
+#define TST03_IMAGE_SIZE 40
+#define TST03_NUM_COLS  TST03_IMAGE_SIZE
+#define TST03_NUM_ROWS  TST03_IMAGE_SIZE
+#define TST03_NUM_STAMPS 2
+#define TST03_NUM_STAMPS_COLS TST03_NUM_STAMPS
+#define TST03_NUM_STAMPS_ROWS TST03_NUM_STAMPS
+#define TST03_BORDER  TST03_KERNEL_SIZE
+//#define TST03_FOOTPRINT (((TST03_IMAGE_SIZE - (2 * TST03_BORDER)) / TST03_NUM_STAMPS) - TST03_KERNEL_SIZE)
+#define TST03_FOOTPRINT  4
+#define TST03_PSF_WIDTH  (TST03_FOOTPRINT/2 - 1)
+
+psBool genObject(psImage *tstImg,
+                 psImage *refImg,
+                 psS32 xMin,
+                 psS32 xMax,
+                 psS32 yMin,
+                 psS32 yMax)
+{
+    if (0) {
+        for (psS32 y = yMin ; y <= yMax ; y++) {
+            for (psS32 x = xMin ; x <= xMax ; x++) {
+                refImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = TST03_THRESHOLD + (psF32) (y/10 + x/10);
+                tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = TST03_THRESHOLD + (psF32) (y/10 + x/10);
+                refImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = 100.0;
+                tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = 100.0;
+
+                // Add some noise for the test image.
+                tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2]+= genRanFloat(0.0, 1.0);
+
+            }
+        }
+
+    } else {
+        if (((1 + 2 * TST03_PSF_WIDTH) > (xMax - xMin)) ||
+                ((1 + 2 * TST03_PSF_WIDTH) > (yMax - yMin))) {
+            printf("INCORRECT TEST CONFIGURATION: TST03_PSF_WIDTH is too big.\n");
+            printf("TST03_PSF_WIDTH is %d: (xMin - xMax) is %d\n", TST03_PSF_WIDTH, (xMax - xMin));
+            return(FALSE);
+        }
+        //
+        // HEY
+        // This code basically creates a peak at the center of the stamp with
+        // a height of TST03_PSF_MAX
+        //
+        psS32 xCenter = (xMax + xMin) / 2;
+        psS32 yCenter = (yMax + yMin) / 2;
+        psF32 subImageWidth = 1.0 + sqrtf(PS_SQR(((psF32) ((yMax-yMin)/2))) + PS_SQR(((psF32) ((xMax-xMin)/2))));
+        for (psS32 y = yMin ; y <= yMax ; y++) {
+            for (psS32 x = xMin ; x <= xMax ; x++) {
+                psF32 dist = sqrtf(PS_SQR((psF32) (y - yCenter)) + PS_SQR((psF32) (x - xCenter)));
+                psF32 pixel = TST03_PSF_MAX * PS_SQR(((psF32) (subImageWidth - dist)) / ((psF32) subImageWidth));
+                if (pixel < 0.0) {
+                    pixel = 0.0;
+                }
+                refImg->data.F32[y][x] = pixel;
+                tstImg->data.F32[y][x] = pixel;
+                // Add some noise for the test image.
+                tstImg->data.F32[y][x]+= genRanFloat(0.0, 1.0);
+            }
+        }
+    }
+    return(TRUE);
+}
+
 /*******************************************************************************
 NOTE: This function returns FALSE if there were no errors.
@@ -548 +614 @@
             psS32 xMax = PS_MIN(numCols-1, (border + (i + 1) * (numCols - 2.0 * border) / xNum) - 1);
 
-            // Set the center pixel in the stamp.
-            refImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = TST03_THRESHOLD + (psF32) (i + j);
-            tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = TST03_THRESHOLD + (psF32) (i + j);
-            refImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = 100.0;
-            tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = 100.0;
-
-            // Add some noise for the test image.
-            tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2]+= genRanFloat(0.0, 1.0);
+            if (0) {
+                // Set the center pixel in the stamp.
+                refImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = TST03_THRESHOLD + (psF32) (i + j);
+                tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = TST03_THRESHOLD + (psF32) (i + j);
+                refImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = 100.0;
+                tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2] = 100.0;
+
+                // Add some noise for the test image.
+                tstImg->data.F32[(yMax+yMin)/2][(xMax+xMin)/2]+= genRanFloat(0.0, 1.0);
+            } else {
+                genObject(tstImg, refImg, xMin, xMax, yMin, yMax);
+            }
         }
     }
@@ -593 +663 @@
         //-------------------------------------------------------------------------
         printf("Calling with a NULL psArray stamps.  Should generate error, return FALSE.\n");
-        psBool rc = pmSubtractionCalculateEquation(NULL, refImg, tstImg, myKernels, TST03_KERNEL_SIZE);
+        psBool rc = pmSubtractionCalculateEquation(NULL, refImg, tstImg, myKernels, TST03_FOOTPRINT);
         if (rc == TRUE) {
             printf("TEST ERROR: pmSubtractionCalculateEquation() returned TRUE.\n");
@@ -601 +671 @@
         //-------------------------------------------------------------------------
         printf("Calling with a NULL reference images.  Should generate error, return FALSE.\n");
-        rc = pmSubtractionCalculateEquation(stamps, NULL, tstImg, myKernels, TST03_KERNEL_SIZE);
+        rc = pmSubtractionCalculateEquation(stamps, NULL, tstImg, myKernels, TST03_FOOTPRINT);
         if (rc == TRUE) {
             printf("TEST ERROR: pmSubtractionCalculateEquation() returned TRUE.\n");
@@ -609 +679 @@
         //-------------------------------------------------------------------------
         printf("Calling with a NULL input images.  Should generate error, return FALSE.\n");
-        rc = pmSubtractionCalculateEquation(stamps, refImg, NULL, myKernels, TST03_KERNEL_SIZE);
+        rc = pmSubtractionCalculateEquation(stamps, refImg, NULL, myKernels, TST03_FOOTPRINT);
         if (rc == TRUE) {
             printf("TEST ERROR: pmSubtractionCalculateEquation() returned TRUE.\n");
@@ -617 +687 @@
         //-------------------------------------------------------------------------
         printf("Calling with a NULL kernel basis functions.  Should generate error, return FALSE.\n");
-        rc = pmSubtractionCalculateEquation(stamps, refImg, tstImg, NULL, TST03_KERNEL_SIZE);
+        rc = pmSubtractionCalculateEquation(stamps, refImg, tstImg, NULL, TST03_FOOTPRINT);
         if (rc == TRUE) {
             printf("TEST ERROR: pmSubtractionCalculateEquation() returned TRUE.\n");
@@ -626 +696 @@
         printf("Calling with acceptable input parameters.  Should return TRUE.\n");
 
-        rc = pmSubtractionCalculateEquation(stamps, refImg, tstImg, myKernels, TST03_KERNEL_SIZE);
+        rc = pmSubtractionCalculateEquation(stamps, refImg, tstImg, myKernels, TST03_FOOTPRINT);
         if (rc != TRUE) {
             printf("TEST ERROR: pmSubtractionCalculateEquation() returned FALSE.\n");
@@ -641 +711 @@
 
             //-------------------------------------------------------------------------
-            printf("Calling with pmSubtractionSolveEquation() acceptable input parameters.  Should return non-NULL.\n");
+            printf("Calling pmSubtractionSolveEquation() with acceptable input parameters.  Should return non-NULL.\n");
             solution = pmSubtractionSolveEquation(NULL, stamps);
             if (solution == NULL) {
                 printf("TEST ERROR: pmSubtractionSolveEquation() returned NULL.\n");
                 testStatus = true;
+            } else {
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionRejectStamps() with acceptable input parameters.  Should return TRUE.\n");
+                rc = pmSubtractionRejectStamps(stamps, maskImg, 0xff, TST03_FOOTPRINT, 1.0, refImg,
+                                               tstImg, solution, myKernels);
+                if (rc != TRUE) {
+                    printf("TEST ERROR: pmSubtractionRejectStamps() returned FALSE.\n");
+                    testStatus = true;
+                } else {
+                    printf("The solution vector is:\n");
+                    for (psS32 i = 0 ; i < solution->n ; i++) {
+                        printf("(%.2f) ", solution->data.F32[i]);
+                    }
+                    printf("\n");
+                }
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionKernelImage() with NULL solution.  Should generate error, return NULL.\n");
+                psImage *kernelImg = pmSubtractionKernelImage(NULL, NULL, myKernels, 0.0, 0.0);
+                if (kernelImg != NULL) {
+                    printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
+                    testStatus = true;
+                }
+                free(kernelImg);
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionKernelImage() with NULL kernels.  Should generate error, return NULL.\n");
+                kernelImg = pmSubtractionKernelImage(NULL, solution, NULL, 0.0, 0.0);
+                if (kernelImg != NULL) {
+                    printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
+                    testStatus = true;
+                }
+                free(kernelImg);
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionKernelImage() unallowable x value.  Should generate error, return NULL.\n");
+                kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, -2.0, 0.0);
+                if (kernelImg != NULL) {
+                    printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
+                    testStatus = true;
+                }
+                free(kernelImg);
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionKernelImage() unallowable x value.  Should generate error, return NULL.\n");
+                kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 2.0, 0.0);
+                if (kernelImg != NULL) {
+                    printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
+                    testStatus = true;
+                }
+                free(kernelImg);
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionKernelImage() unallowable y value.  Should generate error, return NULL.\n");
+                kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 0.0, -2.0);
+                if (kernelImg != NULL) {
+                    printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
+                    testStatus = true;
+                }
+                free(kernelImg);
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionKernelImage() unallowable y value.  Should generate error, return NULL.\n");
+                kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 0.0, 2.0);
+                if (kernelImg != NULL) {
+                    printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
+                    testStatus = true;
+                }
+                free(kernelImg);
+
+                //-------------------------------------------------------------------------
+                printf("Calling pmSubtractionKernelImage() with acceptable input parameters.  Should return a psImage.\n");
+                kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 0.5, 0.5);
+                if (kernelImg == NULL) {
+                    printf("TEST ERROR: pmSubtractionKernelImage() returned NULL.\n");
+                    testStatus = true;
+                } else {
+                    for (psS32 row = 0 ; row < kernelImg->numRows; row++) {
+                        for (psS32 col = 0 ; col < kernelImg->numCols; col++) {
+                            printf("%f ", kernelImg->data.F32[row][col]);
+                        }
+                        printf("\n");
+                    }
+                }
+                free(kernelImg);
+
+                psFree(solution);
             }
-
-            //-------------------------------------------------------------------------
-            printf("Calling with pmSubtractionRejectStamps() acceptable input parameters.  Should return TRUE.\n");
-            rc = pmSubtractionRejectStamps(stamps, maskImg, 0xff, TST03_FOOTPRINT, 1.0, refImg, tstImg, solution, myKernels);
-            if (rc != TRUE) {
-                printf("TEST ERROR: pmSubtractionRejectStamps() returned FALSE.\n");
-                testStatus = true;
-            }
-
-            //-------------------------------------------------------------------------
-            printf("Calling pmSubtractionKernelImage() with NULL solution.  Should generate error, return NULL.\n");
-            psImage *kernelImg = pmSubtractionKernelImage(NULL, NULL, myKernels, 0.0, 0.0);
-            if (kernelImg != NULL) {
-                printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
-                testStatus = true;
-            }
-            free(kernelImg);
-
-            //-------------------------------------------------------------------------
-            printf("Calling pmSubtractionKernelImage() with NULL kernels.  Should generate error, return NULL.\n");
-            kernelImg = pmSubtractionKernelImage(NULL, solution, NULL, 0.0, 0.0);
-            if (kernelImg != NULL) {
-                printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
-                testStatus = true;
-            }
-            free(kernelImg);
-
-            //-------------------------------------------------------------------------
-            printf("Calling pmSubtractionKernelImage() unallowable x value.  Should generate error, return NULL.\n");
-            kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, -2.0, 0.0);
-            if (kernelImg != NULL) {
-                printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
-                testStatus = true;
-            }
-            free(kernelImg);
-
-            //-------------------------------------------------------------------------
-            printf("Calling pmSubtractionKernelImage() unallowable x value.  Should generate error, return NULL.\n");
-            kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 2.0, 0.0);
-            if (kernelImg != NULL) {
-                printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
-                testStatus = true;
-            }
-            free(kernelImg);
-
-            //-------------------------------------------------------------------------
-            printf("Calling pmSubtractionKernelImage() unallowable y value.  Should generate error, return NULL.\n");
-            kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 0.0, -2.0);
-            if (kernelImg != NULL) {
-                printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
-                testStatus = true;
-            }
-            free(kernelImg);
-
-            //-------------------------------------------------------------------------
-            printf("Calling pmSubtractionKernelImage() unallowable y value.  Should generate error, return NULL.\n");
-            kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 0.0, 2.0);
-            if (kernelImg != NULL) {
-                printf("TEST ERROR: pmSubtractionKernelImage() returned non-NULL.\n");
-                testStatus = true;
-            }
-            free(kernelImg);
-
-            //-------------------------------------------------------------------------
-            printf("Calling with pmSubtractionKernelImage() acceptable input parameters.  Should return a psImage.\n");
-            kernelImg = pmSubtractionKernelImage(NULL, solution, myKernels, 0.0, 0.0);
-            if (kernelImg == NULL) {
-                printf("TEST ERROR: pmSubtractionKernelImage() returned NULL.\n");
-                testStatus = true;
-            }
-            free(kernelImg);
-
-            psFree(solution);
         }
     }
@@ -739 +824 @@
     bool testStatus = false;
 
-    testStatus|= testSubCalcEqu(100, 100, 2, 2, 3, PM_SUBTRACTION_KERNEL_POIS);
-    testStatus|= testSubCalcEqu(100, 100, 2, 2, 3, PM_SUBTRACTION_KERNEL_ISIS);
+    testStatus|= testSubCalcEqu(TST03_NUM_COLS,
+                                TST03_NUM_ROWS,
+                                TST03_NUM_STAMPS_COLS,
+                                TST03_NUM_STAMPS_ROWS,
+                                TST03_BORDER,
+                                PM_SUBTRACTION_KERNEL_POIS);
+
+    testStatus|= testSubCalcEqu(TST03_NUM_COLS,
+                                TST03_NUM_ROWS,
+                                TST03_NUM_STAMPS_COLS,
+                                TST03_NUM_STAMPS_ROWS,
+                                TST03_BORDER,
+                                PM_SUBTRACTION_KERNEL_ISIS);
 
     return(!testStatus);
 }
 
+//This code will