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

Timestamp:

Dec 17, 2004, 4:27:42 PM (22 years ago)

Author:

gusciora

Message:

This code appears to work correctly.

File:

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

trunk/psModules/src/pmSubtractSky.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-16 00:47:01 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-18 02:27:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -17 +17 @@
 #include "pslib.h"
 #include "psConstants.h"
-
-// XXX: this is pmFit in pmSubtractBias.c, named psFit here.
-typedef enum {
-    PM_FIT_NONE,                              ///< No fit
-    PM_FIT_POLYNOMIAL,                        ///< Fit polynomial
-    PM_FIT_SPLINE                             ///< Fit cubic splines
-} psFit;
-
+#include "pmSubtractSky.h"
+
+/******************************************************************************
+ *****************************************************************************/
 int p_psDetermineNumBits(unsigned int data)
 {
@@ -40 +36 @@
 }
 
+/******************************************************************************
+ *****************************************************************************/
 psU64 getHighestPriorityStatOption(psU64 statOptions)
 {
@@ -60 +58 @@
 }
 
-
 /******************************************************************************
 psImage *binImage(origImage, binFactor, statOptions): This routine takes an
@@ -67 +64 @@
  
 XXX: use static vectors for myStats, binVector and binMask.
+XXX: I coded this before I was aware of a psLib reBin function.  I don't
+use this function in this module.  I'm keeping it here in the event that
+requirements change and we might need a custom reBin function.
  *****************************************************************************/
 psImage *binImage(psImage *origImage,
@@ -72 +72 @@
                   psStatsOptions statOptions)
 {
+    if (binFactor <= 0) {
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: binImage(): binFactor is %d\n", binFactor);
+        return(origImage);
+    }
+    if (binFactor == 1) {
+        return(origImage);
+    }
+
     psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32);
     psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8);
     psStats *myStats = psStatsAlloc(statOptions);
-
-    if (binFactor <= 0) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: binImage(): binFactor is %d\n", binFactor);
-        return(origImage);
-    }
-    if (binFactor == 1) {
-        return(origImage);
-    }
 
     for (int row = 0; row < origImage->numRows ; row+=binFactor) {
@@ -115 +115 @@
         }
     }
-    psFree(myStats);
     psFree(binVector);
     psFree(binMask);
+    psFree(myStats);
 
     return(origImage);
@@ -123 +123 @@
 
 /******************************************************************************
-buildPolyTerms(): this routine computes a 2-D array polyTerms[] that holds
+ *****************************************************************************/
+psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder)
+{
+    psS32 maxOrder = PS_MAX(xOrder, yOrder);
+    psS32 localPolyTerms = 0;
+    psS32 order = 0;
+    psS32 num=0;
+
+    for (order=0;order<=maxOrder;order++) {
+        for (num=0;num<=order;num++) {
+            if (((order-num) <= xOrder) && (num <= yOrder)) {
+                localPolyTerms++;
+            }
+        }
+    }
+    return(localPolyTerms);
+}
+
+/******************************************************************************
+buildPolyTerms(): this routine computes a 2-D array polyTerms[][] that holds
 terms for the polynomial that is used to model the sky background.  We use
-this array primarily for convenience in many computations involving that sky
-model polynomials.  It is defined as:
+this array primarily for convenience in computations involving sky model
+polynomials.  It is defined as:
     polyTerms[i][0] = the power to which X is raised in the i-th term of in an
     poly-order sky background polynomial.
@@ -133 +152 @@
     poly-order sky background polynomial.
  *****************************************************************************/
-psS32 **buildPolyTerms(psS32 polyOrder)
+psS32 **buildPolyTerms(psS32 xOrder, psS32 yOrder)
 {
     psS32 i=0;
     psS32 order = 0;
     psS32 num=0;
-    psS32 localPolyTerms= (((polyOrder+1) * (polyOrder + 2)) / 2);
-
-    // We create the data structure which we hold the xy order of each coeff.
+    psS32 localPolyTerms = CalculatePolyTerms(xOrder, yOrder);
+    psS32 maxOrder = PS_MAX(xOrder, yOrder);
+
+    // Create the data structure which we hold the xy order of each coeff.
     psS32 **polyTerms = (psS32 **) psAlloc(localPolyTerms * sizeof(psS32 *));
     for (i=0; i < localPolyTerms ; i++) {
@@ -150 +170 @@
     // calculates the power to which x/y are raised in that i-th term.
     // We first do the 0-order terms, then the 1-order terms, etc.
-    for (order=0;order<=polyOrder;order++) {
+    for (order=0;order<=maxOrder;order++) {
         for (num=0;num<=order;num++) {
-            polyTerms[i][0] = order-num;
-            polyTerms[i][1] = num;
-            i++;
-        }
+            if (((order-num) <= xOrder) && (num <= yOrder)) {
+                polyTerms[i][0] = order-num;
+                polyTerms[i][1] = num;
+                i++;
+            }
+        }
+    }
+
+    // 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]);
     }
 
@@ -161 +189 @@
 }
 
-#define PS_TWENTY 20
-
+/******************************************************************************
+This procedure calculates various combinations of powers of x and y and stores
+them in the data structure sums[][].  After it completes:
+ 
+    sums[i][j] == x^i * y^j
+ *****************************************************************************/
 // XXX: Use variable size arrays for polynomial sums.
-/** @brief This procedure calculates various combinations of powers of x and y
- *   and stores them in the data structure sums[][].  After it completes:
- *          sums[i][j] == x^i * y^j
- */
+#define PS_MAX_POLYNOMIAL_ORDER 20
 void buildSums(psF64 x,
                psF64 y,
-               psF64 sums[PS_TWENTY][PS_TWENTY],
-               psS32 polyOrder)
+               psF64 sums[PS_MAX_POLYNOMIAL_ORDER][PS_MAX_POLYNOMIAL_ORDER],
+               psS32 xOrder,
+               psS32 yOrder)
 {
     psS32 i = 0;
@@ -180 +210 @@
     xSum = 1.0;
     ySum = 1.0;
-    for(i=0;i<=polyOrder;i++) {
+    for(i=0;i<=xOrder;i++) {
         ySum = xSum;
-        for(j=0;j<=polyOrder;j++) {
+        for(j=0;j<=yOrder;j++) {
             sums[i][j] = ySum;
             ySum*= y;
@@ -191 +221 @@
 
 /******************************************************************************
- 
  *****************************************************************************/
 psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly,
@@ -197 +226 @@
                                    psImage *maskImage)
 {
+    psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 4,
+            "---- ImageFitPolynomial() begin ----\n");
     PS_POLY_CHECK_NULL(myPoly, NULL);
     PS_POLY_CHECK_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL);
-    PS_INT_CHECK_NON_EQUALS(myPoly->nX, myPoly->nY, NULL);
     PS_IMAGE_CHECK_NULL(binnedImage, NULL);
     PS_IMAGE_CHECK_EMPTY(binnedImage, NULL);
@@ -207 +237 @@
     PS_IMAGE_CHECK_TYPE(maskImage, PS_TYPE_U8, NULL);
     PS_IMAGE_CHECK_SIZE_EQUAL(binnedImage, maskImage, NULL);
+    psS32 oldPolyX = -1;
+    psS32 oldPolyY = -1;
+
+    // The matrix equations become singular if there are more powers of X
+    // in myPoly then there are rows of the image.  I think.  Similarly for
+    // powers of Y and columns.  So.  Here we reduce the complexity of the
+    // polynomial if there are not enough rows/columns in the input image.
+
+    if (myPoly->nX > binnedImage->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) {
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in y-dimension.\n");
+        oldPolyY = myPoly->nY;
+        myPoly->nY = binnedImage->numCols;
+    }
     psS32 i;
     psS32 j;
@@ -213 +263 @@
     psS32 aRow;
     psS32 aCol;
-    psS32 polyOrder = myPoly->nX;
-    psS32 localPolyTerms= (((polyOrder+1) * (polyOrder + 2)) / 2);
-    psS32 **polyTerms = buildPolyTerms(polyOrder);
-    psF64 sums[PS_TWENTY][PS_TWENTY];
-    psImage *A = psImageAlloc(localPolyTerms+1, localPolyTerms+1, PS_TYPE_F64);
-    psVector *B = psVectorAlloc(localPolyTerms+1, PS_TYPE_F64);
-    psVector *outPerm = NULL;
+    // XXX: Document this.  The myPoly->nX and ->nY terms are actual 1 larger
+    // than the order of the polynomial.
+    psS32 **polyTerms = buildPolyTerms(myPoly->nX-1, myPoly->nY-1);
+    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);
+    psVector *B = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
+    psVector *outPerm = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
 
     for(i=0;i<A->numRows;i++) {
@@ -232 +284 @@
     for (x=0;x<binnedImage->numRows;x++) {
         for (y=0;y<binnedImage->numCols;y++) {
-            if (maskImage->data.U8[x][y] != 0) {
-                buildSums((psF64) x, (psF64) y, sums, polyOrder);
+            if (maskImage->data.U8[x][y] == 0) {
+                buildSums((psF64) x, (psF64) y, sums, myPoly->nX-1, myPoly->nY-1);
 
                 /************************************************************
-                Equation (7) from the ADD describes 16 linear equations.
+                Equation (7) from the pilot ADD describes 16 linear equations.
                 The i-th equation is simply the partial derivative of the
                 sky background polynomial (1) w.r.t. to the i-th term in
@@ -251 +303 @@
                     }
                 }
-
                 // Build the B[] vector, which is the right-hand side of (7).
-                for (i=0;i<=localPolyTerms;i++) {
+                for (i=0;i<localPolyTerms;i++) {
                     B->data.F64[i]+= binnedImage->data.F32[x][y] *
                                      sums[ polyTerms[i][0] ][ polyTerms[i][1] ];
@@ -260 +311 @@
         }
     }
-    psImage *ALUD = psMatrixLUD(NULL, outPerm, A);
-    psVector *C = psMatrixLUSolve(C, ALUD, B, outPerm);
-
+
+    // 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", 6,
+                    "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]);
+    }
+
+
+    // Solve the matrix equations for the polynomial coefficients C.
+    Aout = psMatrixLUD(Aout, outPerm, A);
+    psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
+    psMatrixLUSolve(C, Aout, B, outPerm);
+
+    // Set the appropriate coefficients in the myPoly structure.
     for (i=0;i<localPolyTerms;i++) {
         myPoly->coeff[ polyTerms[i][0] ][ polyTerms[i][1] ] = C->data.F64[i];
-    }
-
+        psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6,
+                "myPoly->coeff[%d][%d] is %f\n", polyTerms[i][0], polyTerms[i][1], myPoly->coeff[ polyTerms[i][0] ][ polyTerms[i][1] ]);
+    }
+
+    // Free data structures that were allocated in this module.
     for (i=0;i<localPolyTerms;i++) {
         psFree(polyTerms[i]);
@@ -272 +344 @@
     psFree(polyTerms);
     psFree(A);
-    psFree(ALUD);
+    psFree(Aout);
     psFree(B);
     psFree(C);
     psFree(outPerm);
 
+    // We restore the original size of the polynomial and set remaining
+    // coefficients to 0.0.
+    if (oldPolyX != -1) {
+        myPoly->nX = oldPolyX;
+        for (i=oldPolyX ; i < myPoly->nX ; i++) {
+            for (j=0;j<myPoly->nY ; j++) {
+                myPoly->coeff[i][j] = 0.0;
+            }
+        }
+    }
+    if (oldPolyY != -1) {
+        myPoly->nY = oldPolyY;
+        for (i=0 ; i < myPoly->nX ; i++) {
+            for (j=oldPolyY;j<myPoly->nY ; j++) {
+                myPoly->coeff[i][j] = 0.0;
+            }
+        }
+    }
+
+    psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 4,
+            "---- ImageFitPolynomial() end successfully ----\n");
     return(myPoly);
 }
@@ -293 +386 @@
                          float clipSD)
 {
+    psTrace(".psModule.pmSubtractSky", 4,
+            "---- pmSubtractSky() begin ----\n");
+
     // Return the original input readout if the fit specs are poorly defined.
     if ((fitSpec == NULL) ||
@@ -299 +395 @@
     }
     psImage *origImage = in->image;
-    psImage *maskImage = in->mask;
     psImage *binnedImage = NULL;
     psPolynomial2D *myPoly;
+    psImage *binnedMaskImage = NULL;
 
     psStatsOptions statOptions = stats->options;
-    if (1 < p_psDetermineNumBits(statOptions)) {
+    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.
+        stats->options = statOptions;
         if (statOptions <= 0) {
             psLogMsg(__func__, PS_LOG_WARN,
@@ -314 +412 @@
     }
 
+    //
+    // 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,
+                 "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor);
+    }
+    if (stats == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: pmSubtractSky(): input parameter stats is NULL\n");
+    }
+
+    //
     // Bin the input image according to input parameters.
-    if ((binFactor <= 0) || (stats == NULL)) {
-        // Simply use the original image: no binning.
-        if (binFactor <= 0) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor);
-        }
-        if (stats == NULL) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractSky(): input parameter stats is NULL\n");
-        }
-        binnedImage = origImage;
+    //
+    if ((binFactor <= 1) ||
+            (stats == NULL) ||
+            (0 == p_psDetermineNumBits(statOptions))) {
+        binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
     } else {
-        // Call a private function to do the binning.
-        binnedImage = binImage(origImage, binFactor, statOptions);
-    }
-
+        // Add the original image mask in here.
+        binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
+        binnedImage = psImageRebin(NULL, binnedImage, NULL, 0, binFactor, stats);
+    }
+    psTrace(".psModule.pmSubtractSky", 4,
+            "binnedImage size is (%d, %d)\n", binnedImage->numRows, binnedImage->numCols);
+
+    //
     // Clip pixels that are outside the acceptable range.
+    //
     if (clipSD <= 0.0) {
         psLogMsg(__func__, PS_LOG_WARN,
@@ -342 +456 @@
         myStats =  psImageStats(myStats, binnedImage, NULL, 0);
         p_psGetStatValue(myStats, &binnedMean);
+        psTrace(".psModule.pmSubtractSky", 8,
+                "binned Mean is %f\n", binnedMean);
 
         myStats->options = PS_STAT_SAMPLE_STDEV;
@@ -347 +463 @@
         p_psGetStatValue(myStats, &binnedStdev);
         psFree(myStats);
+        psTrace(".psModule.pmSubtractSky", 8,
+                "binned StDev is %f\n", binnedStdev);
 
         // Clip all pixels which are more than clipSD sigmas from the mean.
@@ -352 +470 @@
         // XXX: we must unset this later since we modify the image mask.
         // XXX: Determine which pixels, mask or image, should be clipped.
+
+        binnedMaskImage = psImageAlloc(binnedImage->numCols,
+                                       binnedImage->numRows,
+                                       PS_TYPE_U8);
+
+        psTrace(".psModule.pmSubtractSky", 8,
+                "clipSD is %f\n", clipSD);
+
         for (int row = 0; row < binnedImage->numRows ; row++) {
             for (int col = 0; col < binnedImage->numCols ; col++) {
                 if (fabs(binnedImage->data.F32[row][col] - binnedMean) >
                         (clipSD * binnedStdev)) {
-                    maskImage->data.U8[row][col] = 1;
+                    binnedMaskImage->data.U8[row][col] = 1;
                 }
             }
@@ -362 +488 @@
     }
 
-    // XXX: fit the polynomial to the binned image
+    //
+    // Fit the polynomial to the binned image
+    //
     if (fit == PM_FIT_POLYNOMIAL) {
         myPoly = (psPolynomial2D *) fitSpec;
-
-        // XXX Ensure that the polynomial is of type Chebyshev.
-        myPoly = ImageFitPolynomial(myPoly, binnedImage, maskImage);
-
-        // XXX Do we need to do something with ordinate scaling if Chebyshev?
-        binnedImage = psImageEvalPolynomial(binnedImage, myPoly);
-    }
-
+        PS_POLY_CHECK_NULL(myPoly, NULL);
+        PS_POLY_CHECK_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL);
+
+        myPoly = ImageFitPolynomial(myPoly, binnedImage, binnedMaskImage);
+
+        if (myPoly != NULL) {
+            // XXX:Add ordinary polynomials to psImageEvalPolynomial()
+            for (int row = 0; row < binnedImage->numRows ; row++) {
+                for (int col = 0; col < binnedImage->numCols ; col++) {
+                    binnedImage->data.F32[row][col] =
+                        psPolynomial2DEval(myPoly, (psF32) row, (psF32) col);
+                    psTrace(".psModule.pmSubtractSky", 8,
+                            "binned Image[%d][%d] is %f\n",
+                            row, col, binnedImage->data.F32[row][col]);
+                }
+            }
+        } else {
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: pmSubtractSky(): could not model sky with a polynomial.\n");
+            psFree(binnedMaskImage);
+            if (!((binFactor <= 1) || (stats == NULL))) {
+                psFree(binnedImage);
+            }
+            return(in);
+        }
+    } else {
+        //XXX: error
+    }
+
+    //
     //Subtract the polynomially fitted image from the original image
-    if (binFactor <= 0) {
+    //
+    if (binFactor <= 1) {
         // The binned image is the same size as the original image.
         for (int row = 0; row < origImage->numRows ; row++) {
@@ -384 +535 @@
         for (int row = 0; row < origImage->numRows ; row++) {
             for (int col = 0; col < origImage->numCols ; col++) {
-                int binRow = row / binFactor;
-                int binCol = col / binFactor;
-                origImage->data.F32[row][col]-= binnedImage->data.F32[binRow][binCol];
-            }
-        }
-
-    }
-
-    if (!((binFactor <= 0) || (stats == NULL))) {
+                // We calculate the F32 value of the pixel coordinates in the
+                // binned image and then use a pixel interpolation routine to
+                // determine the value of the pixel at that location.
+                psF32 binRowF64 = ((psF32) row) / ((psF32) binFactor);
+                psF32 binColF64 = ((psF32) col) / ((psF32) binFactor);
+
+                // We add 0.5 to the pixel locations since the pixel
+                // interpolation routine defines the location of pixel
+                // (i, j) as (i+0.5, j+0.5).
+                binRowF64+= 0.5;
+                binColF64+= 0.5;
+
+                psF32 binPixel = (psF32) psImagePixelInterpolate(
+                                     binnedImage, binRowF64, binColF64,
+                                     NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
+                origImage->data.F32[row][col]-= binPixel;
+
+                psTrace(".psModule.pmSubtractSky", 8,
+                        "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f (%f)\n",
+                        row, col, binRowF64-0.5, binColF64-0.5, binPixel,
+                        binnedImage->data.F32[(psS32)binRowF64][(psS32)binColF64]);
+            }
+        }
+
+    }
+
+    psFree(binnedMaskImage);
+    if (!((binFactor <= 1) || (stats == NULL))) {
         psFree(binnedImage);
     }
 
+    psTrace(".psModule.pmSubtractSky", 4,
+            "---- pmSubtractSky() exit successfully ----\n");
     return(in);
 }