Changeset 2760

Timestamp:

Dec 20, 2004, 11:38:46 AM (22 years ago)

Author:

gusciora

Message:

Working version.

File:

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

trunk/psModules/src/pmSubtractSky.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-18 02:27:42 $
+ *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-20 21:38:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 
 /******************************************************************************
- *****************************************************************************/
-int p_psDetermineNumBits(unsigned int data)
-{
-    int i;
-    unsigned int tmpData = data;
-    int numBits = 0;
-
-    for (i=0;i<sizeof(unsigned int);i++) {
+p_psDetermineNumBits(data): This routine takes an insigned int as an argument
+and returns the number of non-zero bits.
+ *****************************************************************************/
+psS32 p_psDetermineNumBits(psU32 data)
+{
+    psS32 i;
+    psU32 tmpData = data;
+    psS32 numBits = 0;
+
+    for (i=0;i<sizeof(psU32);i++) {
         if (0x0001 && tmpData) {
             numBits++;
@@ -37 +39 @@
 
 /******************************************************************************
+getHighestPriorityStatOption(statOptions): this routine takes as in put a
+psStats->options with multiple options set and returns one with a single
+option set according to the precedence set in the SDRS.
  *****************************************************************************/
 psU64 getHighestPriorityStatOption(psU64 statOptions)
@@ -85 +90 @@
     psStats *myStats = psStatsAlloc(statOptions);
 
-    for (int row = 0; row < origImage->numRows ; row+=binFactor) {
-        for (int col = 0; col < origImage->numCols ; col+=binFactor) {
-            int count = 0;
-            for (int binRow = 0; binRow <= binFactor ; binRow++) {
-                for (int binCol = 0; binCol <= binFactor ; binCol++) {
+    for (psS32 row = 0; row < origImage->numRows ; row+=binFactor) {
+        for (psS32 col = 0; col < origImage->numCols ; col+=binFactor) {
+            psS32 count = 0;
+            for (psS32 binRow = 0; binRow <= binFactor ; binRow++) {
+                for (psS32 binCol = 0; binCol <= binFactor ; binCol++) {
                     if (((row + binRow) < origImage->numRows) &&
                             ((col + binCol) < origImage->numCols)) {
@@ -123 +128 @@
 
 /******************************************************************************
+CalculatePolyTerms(xOrder, yOrder): this routine will calculate the number of
+coefficients (or terms) in a polynomial of order (xOrder, yOrder).
  *****************************************************************************/
 psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder)
@@ -194 +201 @@
  
     sums[i][j] == x^i * y^j
- *****************************************************************************/
-// XXX: Use variable size arrays for polynomial sums.
+ 
+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.
+ *****************************************************************************/
 #define PS_MAX_POLYNOMIAL_ORDER 20
 void buildSums(psF64 x,
@@ -221 +232 @@
 
 /******************************************************************************
+ImageFitPolynomial(myPoly, binnedImage, 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
+symmetrical in orders of X and Y.  The polynomial must be type
+PS_POLYNOMIAL_ORD.  If there are not enough rows or columns in the input image
+for the order of the polynomial, then that order is reduced.  The algorithm
+used in this routine is based on that of the pilot project ADD, but is not
+documented anywhere.
  *****************************************************************************/
 psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly,
@@ -263 +282 @@
     psS32 aRow;
     psS32 aCol;
-    // XXX: Document this.  The myPoly->nX and ->nY terms are actual 1 larger
+    // XXX: Document this.  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);
@@ -273 +292 @@
     psVector *outPerm = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
 
+    //
+    // Initialize A matrix and B vector.
+    //
     for(i=0;i<A->numRows;i++) {
         for(j=0;j<A->numCols;j++) {
@@ -282 +304 @@
     }
 
+
+    //
+    // We build the A matrix and B vector.
+    //
     for (x=0;x<binnedImage->numRows;x++) {
         for (y=0;y<binnedImage->numCols;y++) {
@@ -315 +341 @@
     for (aRow=0;aRow<localPolyTerms;aRow++) {
         for (aCol=0;aCol<localPolyTerms;aCol++) {
-            psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6,
+            psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 8,
                     "A[%d][%d] is %f\n", aRow, aCol, A->data.F64[aRow][aCol]);
         }
@@ -325 +351 @@
     }
 
-
+    //
     // 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];
@@ -338 +367 @@
     }
 
+    //
     // Free data structures that were allocated in this module.
+    //
     for (i=0;i<localPolyTerms;i++) {
         psFree(polyTerms[i]);
@@ -349 +380 @@
     psFree(outPerm);
 
+    //
     // We restore the original size of the polynomial and set remaining
     // coefficients to 0.0.
+    //
     if (oldPolyX != -1) {
         myPoly->nX = oldPolyX;
@@ -382 +415 @@
                          void *fitSpec,
                          psFit fit,
-                         int binFactor,
+                         psS32 binFactor,
                          psStats *stats,
-                         float clipSD)
-{
+                         psF32 clipSD)
+{
+    PS_READOUT_CHECK_NULL(in, NULL);
+    PS_READOUT_CHECK_EMPTY(in, NULL);
+    PS_READOUT_CHECK_TYPE(in, PS_TYPE_F32, NULL);
     psTrace(".psModule.pmSubtractSky", 4,
             "---- pmSubtractSky() begin ----\n");
@@ -396 +432 @@
     psImage *origImage = in->image;
     psImage *binnedImage = NULL;
-    psPolynomial2D *myPoly;
+    psPolynomial2D *myPoly = NULL;
     psImage *binnedMaskImage = NULL;
-
+    psU32 oldStatOptions = 0;
+
+    //
+    // Determine which statistic to use when binning pixels, if any.
+    //
     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.
+        oldStatOptions = stats->options;
         stats->options = statOptions;
-        if (statOptions <= 0) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractSky(): unallowable stats->options was requested\n");
-            return(in);
-        }
+    }
+    if (0 >= p_psDetermineNumBits(statOptions)) {
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: pmSubtractSky(): no stats->options was requested\n");
+        return(in);
     }
 
@@ -430 +471 @@
     //
     // Bin the input image according to input parameters.
-    //
-    if ((binFactor <= 1) ||
-            (stats == NULL) ||
-            (0 == p_psDetermineNumBits(statOptions))) {
+    // Create a new binned image mask.
+    //
+    if ((binFactor <= 1) || (stats == NULL) || (0 == p_psDetermineNumBits(statOptions))) {
+        // No binning is required here.  Simply create a copy of the image
+        // and a mask.
         binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
+        if (in->mask != NULL) {
+            binnedMaskImage = psImageCopy(binnedMaskImage, in->mask, PS_TYPE_U8);
+        } else {
+            binnedMaskImage = psImageAlloc(binnedImage->numCols,
+                                           binnedImage->numRows,
+                                           PS_TYPE_U8);
+            PS_IMAGE_SET_U8(binnedMaskImage, 0);
+        }
     } else {
-        // Add the original image mask in here.
-        binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
-        binnedImage = psImageRebin(NULL, binnedImage, NULL, 0, binFactor, stats);
+        binnedImage = psImageRebin(NULL, origImage, in->mask, 0, binFactor, stats);
+        binnedMaskImage = psImageAlloc(binnedImage->numCols,
+                                       binnedImage->numRows,
+                                       PS_TYPE_U8);
+        PS_IMAGE_SET_U8(binnedMaskImage, 0);
     }
     psTrace(".psModule.pmSubtractSky", 4,
@@ -456 +508 @@
         myStats =  psImageStats(myStats, binnedImage, NULL, 0);
         p_psGetStatValue(myStats, &binnedMean);
-        psTrace(".psModule.pmSubtractSky", 8,
+        psTrace(".psModule.pmSubtractSky", 6,
                 "binned Mean is %f\n", binnedMean);
 
@@ -463 +515 @@
         p_psGetStatValue(myStats, &binnedStdev);
         psFree(myStats);
-        psTrace(".psModule.pmSubtractSky", 8,
+        psTrace(".psModule.pmSubtractSky", 6,
                 "binned StDev is %f\n", binnedStdev);
 
         // 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.
-
-        binnedMaskImage = psImageAlloc(binnedImage->numCols,
-                                       binnedImage->numRows,
-                                       PS_TYPE_U8);
-
-        psTrace(".psModule.pmSubtractSky", 8,
+        psTrace(".psModule.pmSubtractSky", 6,
                 "clipSD is %f\n", clipSD);
 
-        for (int row = 0; row < binnedImage->numRows ; row++) {
-            for (int col = 0; col < binnedImage->numCols ; col++) {
+        for (psS32 row = 0; row < binnedImage->numRows ; row++) {
+            for (psS32 col = 0; col < binnedImage->numCols ; col++) {
                 if (fabs(binnedImage->data.F32[row][col] - binnedMean) >
                         (clipSD * binnedStdev)) {
@@ -499 +543 @@
 
         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]);
-                }
-            }
+            // Set the pixels in the binned image to that of the polynomial.
+            binnedImage = psImageEvalPolynomial(binnedImage, myPoly);
         } else {
             psLogMsg(__func__, PS_LOG_WARN,
@@ -516 +552 @@
                 psFree(binnedImage);
             }
+            if (oldStatOptions != 0) {
+                stats->options = statOptions;
+            }
             return(in);
         }
@@ -527 +566 @@
     if (binFactor <= 1) {
         // 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++) {
+        for (psS32 row = 0; row < origImage->numRows ; row++) {
+            for (psS32 col = 0; col < origImage->numCols ; col++) {
                 origImage->data.F32[row][col]-= binnedImage->data.F32[row][col];
             }
         }
     } else {
-        for (int row = 0; row < origImage->numRows ; row++) {
-            for (int col = 0; col < origImage->numCols ; col++) {
+        for (psS32 row = 0; row < origImage->numRows ; row++) {
+            for (psS32 col = 0; col < origImage->numCols ; col++) {
                 // We calculate the F32 value of the pixel coordinates in the
                 // binned image and then use a pixel interpolation routine to
@@ -548 +587 @@
 
                 psF32 binPixel = (psF32) psImagePixelInterpolate(
-                                     binnedImage, binRowF64, binColF64,
+                                     binnedImage, binColF64, binRowF64,
                                      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]);
-            }
-        }
-
-    }
-
+                        "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f\n",
+                        row, col, binRowF64-0.5, binColF64-0.5, binPixel);
+            }
+        }
+
+    }
     psFree(binnedMaskImage);
-    if (!((binFactor <= 1) || (stats == NULL))) {
-        psFree(binnedImage);
+    psFree(binnedImage);
+    if (oldStatOptions != 0) {
+        stats->options = statOptions;
     }