Changeset 2810 for trunk/psModules/src/pmSubtractSky.c

Timestamp:

Dec 23, 2004, 10:27:51 AM (22 years ago)

Author:

gusciora

Message:

Fixed a few bugs in pmSubtractSky.c. Added tests for various input parameter
error conditions.

File:

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

trunk/psModules/src/pmSubtractSky.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-21 20:41:29 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-23 20:27:51 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -39 +39 @@
 
 /******************************************************************************
-getHighestPriorityStatOption(statOptions): this routine takes as in put a
+getHighestPriorityStatOption(statOptions): this routine takes as input a
 psStats->options with multiple options set and returns one with a single
 option set according to the precedence set in the SDRS.
@@ -45 +45 @@
 psU64 getHighestPriorityStatOption(psU64 statOptions)
 {
-    switch (statOptions) {
-    case PS_STAT_SAMPLE_MEAN:
+    if (statOptions & PS_STAT_SAMPLE_MEAN) {
         return(PS_STAT_SAMPLE_MEAN);
-    case PS_STAT_SAMPLE_MEDIAN:
+    } else if (statOptions & PS_STAT_SAMPLE_MEDIAN) {
         return(PS_STAT_SAMPLE_MEDIAN);
-    case PS_STAT_CLIPPED_MEAN:
+    } else if (statOptions & PS_STAT_CLIPPED_MEAN) {
         return(PS_STAT_CLIPPED_MEAN);
-    case PS_STAT_ROBUST_MEAN:
+    } else if (statOptions & PS_STAT_ROBUST_MEAN) {
         return(PS_STAT_ROBUST_MEAN);
-    case PS_STAT_ROBUST_MEDIAN:
+    } else if (statOptions & PS_STAT_ROBUST_MEDIAN) {
         return(PS_STAT_ROBUST_MEDIAN);
-    case PS_STAT_ROBUST_MODE:
+    } else if (statOptions & PS_STAT_ROBUST_MODE) {
         return(PS_STAT_ROBUST_MODE);
     }
-    //    psError(XXX);
+    psError(PS_ERR_UNKNOWN, true, "Unallowable option requested for statistically binning image pixels.\n");
     return(-1);
 }
@@ -129 +128 @@
 /******************************************************************************
 CalculatePolyTerms(xOrder, yOrder): this routine will calculate the number of
-coefficients (or terms) in a polynomial of order (xOrder, yOrder).
+coefficients (or terms) in a 2-D polynomial of order (xOrder, yOrder).
  *****************************************************************************/
 psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder)
@@ -187 +186 @@
     }
 
-    // XXX: Get rid of this.
-    for (i=0; i < localPolyTerms ; i++) {
-        psTrace(".psModule.pmSubtractSky.buildPolyTerms", 6,
-                "x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]);
+    if (psTraceGetLevel(".psModule.pmSubtractSky.buildPolyTerms") >= 10) {
+        for (i=0; i < localPolyTerms ; i++) {
+            printf("x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]);
+        }
     }
 
@@ -198 +197 @@
 /******************************************************************************
 This procedure calculates various combinations of powers of x and y and stores
-them in the data structure sums[][].  After it completes:
+them in the data structure p_psPolySums[][].  After it completes:
  
-    sums[i][j] == x^i * y^j
+    p_psPolySums[i][j] == x^i * y^j
  
-XXX: Use a psImage for the sums data structure?
-XXX: Check for non-NULL sums argument?
-XXX: Check for positive x- and yOrder.
-XXX: Use variable size arrays for polynomial sums.
+XXX: Use a psImage for the p_psPolySums data structure?
+XXX: Check for positive x- and y-Order.
+XXX: Use variable size arrays for p_psPolySums[][].
+XXX" Must initialize p_psPolySums[][]?
  *****************************************************************************/
 #define PS_MAX_POLYNOMIAL_ORDER 20
+psF64 p_psPolySums[PS_MAX_POLYNOMIAL_ORDER+1][PS_MAX_POLYNOMIAL_ORDER+1];
 void buildSums(psF64 x,
                psF64 y,
-               psF64 sums[PS_MAX_POLYNOMIAL_ORDER][PS_MAX_POLYNOMIAL_ORDER],
                psS32 xOrder,
                psS32 yOrder)
@@ -224 +223 @@
         ySum = xSum;
         for(j=0;j<=yOrder;j++) {
-            sums[i][j] = ySum;
+            p_psPolySums[i][j] = ySum;
             ySum*= y;
         }
@@ -232 +231 @@
 
 /******************************************************************************
-ImageFitPolynomial(myPoly, binnedImage, maskImage): this private routine takes
+ImageFitPolynomial(myPoly, dataImage, maskImage): this private routine takes
 an input image along with a mask and fits a polynomial to it.  The degree of
 the polynomial is specified by input parameter myPoly, and need not be
@@ -242 +241 @@
  *****************************************************************************/
 psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly,
-                                   psImage *binnedImage,
+                                   psImage *dataImage,
                                    psImage *maskImage)
 {
@@ -249 +248 @@
     PS_POLY_CHECK_NULL(myPoly, NULL);
     PS_POLY_CHECK_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL);
-    PS_IMAGE_CHECK_NULL(binnedImage, NULL);
-    PS_IMAGE_CHECK_EMPTY(binnedImage, NULL);
-    PS_IMAGE_CHECK_TYPE(binnedImage, PS_TYPE_F32, NULL);
+    PS_IMAGE_CHECK_NULL(dataImage, NULL);
+    PS_IMAGE_CHECK_EMPTY(dataImage, NULL);
+    PS_IMAGE_CHECK_TYPE(dataImage, PS_TYPE_F32, NULL);
     PS_IMAGE_CHECK_NULL(maskImage, NULL);
     PS_IMAGE_CHECK_EMPTY(maskImage, NULL);
     PS_IMAGE_CHECK_TYPE(maskImage, PS_TYPE_U8, NULL);
-    PS_IMAGE_CHECK_SIZE_EQUAL(binnedImage, maskImage, NULL);
+    PS_IMAGE_CHECK_SIZE_EQUAL(dataImage, maskImage, NULL);
     psS32 oldPolyX = -1;
     psS32 oldPolyY = -1;
@@ -264 +263 @@
     // polynomial if there are not enough rows/columns in the input image.
 
-    if (myPoly->nX > binnedImage->numRows) {
+    if (myPoly->nX > dataImage->numRows) {
         psLogMsg(__func__, PS_LOG_WARN,
                  "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in x-dimension.\n");
         oldPolyX = myPoly->nX;
-        myPoly->nX = binnedImage->numRows;
-    }
-    if (myPoly->nY > binnedImage->numCols) {
+        myPoly->nX = dataImage->numRows;
+    }
+    if (myPoly->nY > dataImage->numCols) {
         psLogMsg(__func__, PS_LOG_WARN,
                  "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in y-dimension.\n");
         oldPolyY = myPoly->nY;
-        myPoly->nY = binnedImage->numCols;
+        myPoly->nY = dataImage->numCols;
     }
     psS32 i;
@@ -282 +281 @@
     psS32 aRow;
     psS32 aCol;
-    // XXX: Document this.  The myPoly->nX and ->nY terms are actually 1 larger
+    // The myPoly->nX and ->nY terms are actually 1 larger
     // than the order of the polynomial.
     psS32 **polyTerms = buildPolyTerms(myPoly->nX-1, myPoly->nY-1);
+    // We determine how many coefficients will be in the polynomial that we
+    // are fitting to this image.
     psS32 localPolyTerms = CalculatePolyTerms(myPoly->nX-1, myPoly->nY-1);
-    psF64 sums[PS_MAX_POLYNOMIAL_ORDER][PS_MAX_POLYNOMIAL_ORDER];
     psImage *A = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64);
     psImage *Aout = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64);
@@ -295 +295 @@
     // Initialize A matrix and B vector.
     //
-    for(i=0;i<A->numRows;i++) {
-        for(j=0;j<A->numCols;j++) {
-            A->data.F64[i][j] = 0.0;
-        }
-    }
-    for(i=0;i<B->n;i++) {
-        B->data.F64[i] = 0.0;
-    }
-
+    PS_IMAGE_SET_F64(A, 0.0);
+    PS_VECTOR_SET_F64(B, 0.0);
 
     //
     // We build the A matrix and B vector.
     //
-    for (x=0;x<binnedImage->numRows;x++) {
-        for (y=0;y<binnedImage->numCols;y++) {
+    for (x=0;x<dataImage->numRows;x++) {
+        for (y=0;y<dataImage->numCols;y++) {
             if (maskImage->data.U8[x][y] == 0) {
-                buildSums((psF64) x, (psF64) y, sums, myPoly->nX-1, myPoly->nY-1);
+                buildSums((psF64) x, (psF64) y, myPoly->nX-1, myPoly->nY-1);
 
                 /************************************************************
+                This code dervies from equation (7) of the pilot ADD.  However,
+                it is not exactly the same in that the order of the polynomial
+                may be different in X And Y.
+
                 Equation (7) from the pilot ADD describes 16 linear equations.
                 The i-th equation is simply the partial derivative of the
@@ -325 +322 @@
                     for (aCol=0;aCol<localPolyTerms;aCol++) {
                         A->data.F64[aRow][aCol]+=
-                            (sums[ polyTerms[aCol][0] ][ polyTerms[aCol][1] ] *
-                             sums[ polyTerms[aRow][0] ][ polyTerms[aRow][1] ]);
+                            (p_psPolySums[ polyTerms[aCol][0] ][ polyTerms[aCol][1] ] *
+                             p_psPolySums[ polyTerms[aRow][0] ][ polyTerms[aRow][1] ]);
                     }
                 }
                 // Build the B[] vector, which is the right-hand side of (7).
                 for (i=0;i<localPolyTerms;i++) {
-                    B->data.F64[i]+= binnedImage->data.F32[x][y] *
-                                     sums[ polyTerms[i][0] ][ polyTerms[i][1] ];
+                    B->data.F64[i]+= dataImage->data.F32[x][y] *
+                                     p_psPolySums[ polyTerms[i][0] ][ polyTerms[i][1] ];
                 }
             }
@@ -338 +335 @@
     }
 
-    // XXX: Put this loop inside a psTrace conditional, somehow.
-    for (aRow=0;aRow<localPolyTerms;aRow++) {
-        for (aCol=0;aCol<localPolyTerms;aCol++) {
-            psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 8,
-                    "A[%d][%d] is %f\n", aRow, aCol, A->data.F64[aRow][aCol]);
-        }
-    }
-    // XXX: Put this loop inside a psTrace conditional, somehow.
-    for (i=0;i<=localPolyTerms;i++) {
-        psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6,
-                "B[%d] is %f\n", i, B->data.F64[i]);
+    if (psTraceGetLevel(".psModule.pmSubtractSky.ImageFitPolynomial") >= 8) {
+        for (aRow=0;aRow<localPolyTerms;aRow++) {
+            for (aCol=0;aCol<localPolyTerms;aCol++) {
+                printf("A[%d][%d] is %f\n", aRow, aCol,
+                       A->data.F64[aRow][aCol]);
+            }
+        }
+
+        for (i=0;i<=localPolyTerms;i++) {
+            printf("B[%d] is %f\n", i, B->data.F64[i]);
+        }
     }
 
     //
     // Solve the matrix equations for the polynomial coefficients C.
+    // XXX: How do we know if these matrix operations were successful?
     //
     Aout = psMatrixLUD(Aout, outPerm, A);
+    PS_IMAGE_CHECK_NULL(Aout, NULL);
+    PS_IMAGE_CHECK_EMPTY(Aout, NULL);
     psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
     psMatrixLUSolve(C, Aout, B, outPerm);
@@ -382 +382 @@
     //
     // We restore the original size of the polynomial and set remaining
-    // coefficients to 0.0.
+    // coefficients to 0.0, if necessary.
     //
     if (oldPolyX != -1) {
@@ -424 +424 @@
     psTrace(".psModule.pmSubtractSky", 4,
             "---- pmSubtractSky() begin ----\n");
+    psStatsOptions statOptions = 0;
 
     // Return the original input readout if the fit specs are poorly defined.
@@ -439 +440 @@
     // Determine which statistic to use when binning pixels, if any.
     //
-    psStatsOptions statOptions = stats->options;
-    if (1 < p_psDetermineNumBits(statOptions)) {
-        //XXX  psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n");
-        statOptions = getHighestPriorityStatOption(statOptions);
-        // XXX: Don't modify input parameter.
-        oldStatOptions = stats->options;
-        stats->options = statOptions;
-    }
-    if (0 >= p_psDetermineNumBits(statOptions)) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: pmSubtractSky(): no stats->options was requested\n");
-        return(in);
+    if (stats != NULL) {
+        statOptions = stats->options;
+        if (1 < p_psDetermineNumBits(statOptions)) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Multiple statistical options have been requested.\n");
+            statOptions = getHighestPriorityStatOption(statOptions);
+
+            // Save old input "stats" parameter.
+            oldStatOptions = stats->options;
+            stats->options = statOptions;
+        }
+        if (0 == p_psDetermineNumBits(statOptions)) {
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: pmSubtractSky(): no stats->options was requested\n");
+        }
     }
 
@@ -456 +459 @@
     // Generate required warning messages.
     //
-    if (0 == p_psDetermineNumBits(statOptions)) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: pmSubtractSky(): no stats->options was requested\n");
-    }
     if (binFactor <= 0) {
         psLogMsg(__func__, PS_LOG_WARN,
@@ -558 +557 @@
         }
     } else {
-        //XXX: error
+        // We shouldn't get here since we check this above.
+        psError(PS_ERR_UNKNOWN, true, "Unallowable fit type.");
     }