Changeset 2719

Timestamp:

Dec 15, 2004, 2:46:22 PM (22 years ago)

Author:

gusciora

Message:

* empty log message *

File:

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

trunk/psModules/src/pmSubtractSky.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-01 21:46:06 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-16 00:46:22 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -17 +17 @@
 #include "psConstants.h"
 
-// XXX: this is pmFit in pmSubtractBias.c, and psFit here.
+// XXX: this is pmFit in pmSubtractBias.c, named psFit here.
 typedef enum {
     PM_FIT_NONE,                              ///< No fit
@@ -38 +38 @@
     return(numBits);
 }
-
-
 
 psU64 getHighestPriorityStatOption(psU64 statOptions)
@@ -123 +121 @@
 }
 
+/******************************************************************************
+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:
+    polyTerms[i][0] = the power to which X is raised in the i-th term of in an
+    poly-order sky background polynomial.
+ 
+    polyTerms[i][1] = the power to which Y is raised in the i-th term of in an
+    poly-order sky background polynomial.
+ *****************************************************************************/
+psS32 **buildPolyTerms(psS32 polyOrder)
+{
+    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 **polyTerms = (psS32 **) psAlloc(localPolyTerms * sizeof(psS32 *));
+    for (i=0; i < localPolyTerms ; i++) {
+        polyTerms[i] = (psS32 *) psAlloc(2 * sizeof(psS32));
+    }
+
+    i=0;
+    // This code segment loops through each term i in the polynomial and
+    // 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 (num=0;num<=order;num++) {
+            polyTerms[i][0] = order-num;
+            polyTerms[i][1] = num;
+            i++;
+        }
+    }
+
+    return(polyTerms);
+}
+
+#define PS_TWENTY 20
+
+// 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
+ */
+void buildSums(psF64 x,
+               psF64 y,
+               psF64 sums[PS_TWENTY][PS_TWENTY],
+               psS32 polyOrder)
+{
+    psS32 i = 0;
+    psS32 j = 0;
+    psF64 xSum = 0.0;
+    psF64 ySum = 0.0;
+
+    xSum = 1.0;
+    ySum = 1.0;
+    for(i=0;i<=polyOrder;i++) {
+        ySum = xSum;
+        for(j=0;j<=polyOrder;j++) {
+            sums[i][j] = ySum;
+            ySum*= y;
+        }
+        xSum*= x;
+    }
+}
+
+/******************************************************************************
+ 
+ *****************************************************************************/
+psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly,
+                                   psImage *binnedImage,
+                                   psImage *maskImage)
+{
+    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);
+    PS_IMAGE_CHECK_TYPE(binnedImage, 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);
+    psS32 i;
+    psS32 j;
+    psS32 x;
+    psS32 y;
+    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;
+
+    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;
+    }
+
+    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);
+
+                /************************************************************
+                Equation (7) from the 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
+                that polynomial.  The i-th equation is stored in row i of
+                matrix A[][] (matrix A[][] has origin (1,1), not (0,0)).  To
+                compute A[i][j] we simply multiply the j-th term of the Sky
+                Background Polynomial (SBP) by the i-th term of SBP.
+                ************************************************************/
+                for (aRow=0;aRow<localPolyTerms;aRow++) {
+                    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] ]);
+                    }
+                }
+
+                // 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] ];
+                }
+            }
+        }
+    }
+    psImage *ALUD = psMatrixLUD(NULL, outPerm, A);
+    psVector *C = psMatrixLUSolve(C, ALUD, B, outPerm);
+
+    for (i=0;i<localPolyTerms;i++) {
+        myPoly->coeff[ polyTerms[i][0] ][ polyTerms[i][1] ] = C->data.F64[i];
+    }
+
+    for (i=0;i<localPolyTerms;i++) {
+        psFree(polyTerms[i]);
+    }
+    psFree(polyTerms);
+    psFree(A);
+    psFree(ALUD);
+    psFree(B);
+    psFree(C);
+    psFree(outPerm);
+
+    return(myPoly);
+}
+
 
 /******************************************************************************
 psReadout pmSubtractSky():
  
-XXX: use static vectors for myStats.
+XXX: use static vectors for myStats, and the binned image
  *****************************************************************************/
 psReadout *pmSubtractSky(psReadout *in,
@@ -136 +292 @@
                          float clipSD)
 {
-    if (fitSpec == NULL) {
+    // Return the original input readout if the fit specs are poorly defined.
+    if ((fitSpec == NULL) ||
+            ((fit == PM_FIT_NONE) || (fit == PM_FIT_SPLINE))) {
         return(in);
     }
-    if ((fit != PM_FIT_POLYNOMIAL) && (fit != PM_FIT_SPLINE)) {
-        // No fit is specified.
-        return(in);
-    }
-
     psImage *origImage = in->image;
     psImage *maskImage = in->mask;
     psImage *binnedImage = NULL;
     psPolynomial2D *myPoly;
-    psSpline1D *mySpline;
 
     psStatsOptions statOptions = stats->options;
@@ -154 +306 @@
         //XXX  psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n");
         statOptions = getHighestPriorityStatOption(statOptions);
+        if (statOptions <= 0) {
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: pmSubtractSky(): unallowable stats->options was requested\n");
+            return(in);
+        }
     }
 
     // Bin the input image according to input parameters.
-    if ((binFactor <= 0) ||
-            (stats == NULL)) {
+    if ((binFactor <= 0) || (stats == NULL)) {
+        // Simply use the original image: no binning.
         if (binFactor <= 0) {
             psLogMsg(__func__, PS_LOG_WARN,
@@ -169 +326 @@
         binnedImage = origImage;
     } else {
+        // Call a private function to do the binning.
         binnedImage = binImage(origImage, binFactor, statOptions);
     }
@@ -175 +333 @@
     if (clipSD <= 0.0) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD)
-        ;
+                 "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD);
     } else {
+        // Determine the mean and standard deviation of the binned image.
         psF64 binnedMean;
         psF64 binnedStdev;
@@ -189 +347 @@
         psFree(myStats);
 
+        // Clip all pixels which are more than clipSD sigmas from the mean.
+        // XXX: Is this correct?  We simply set the mask.
+        // XXX: we must unset this later since we modify the image mask.
+        // XXX: Determine which pixels, mask or image, should be clipped.
         for (int row = 0; row < binnedImage->numRows ; row++) {
             for (int col = 0; col < binnedImage->numCols ; col++) {
@@ -201 +363 @@
     // XXX: fit the polynomial to the binned image
     if (fit == PM_FIT_POLYNOMIAL) {
-        // Fit a polynomial to the old overscan vector.
         myPoly = (psPolynomial2D *) fitSpec;
-        myPoly = psImageFitPolynomial(myPoly, origImage);
+
+        // 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);
-
-    } else if (fit == PM_FIT_SPLINE) {
-        // Fit a spline to the old overscan vector.
-        mySpline = (psSpline1D *) fitSpec;
-        // XXX: What do we do?  We don't have 2-D splines.
-        return(in);
     }
 
     //Subtract the polynomially fitted image from the original image
     if (binFactor <= 0) {
+        // The binned image is the same size as the original image.
         for (int row = 0; row < origImage->numRows ; row++) {
             for (int col = 0; col < origImage->numCols ; col++) {
@@ -231 +390 @@
 
     }
+
+    if (!((binFactor <= 0) || (stats == NULL))) {
+        psFree(binnedImage);
+    }
+
     return(in);
 }