Changeset 1022

Timestamp:

Jun 14, 2004, 9:40:15 AM (22 years ago)

Author:

gusciora

Message:

...

Location:

trunk/psLib

Files:

• src/dataManip/psFunctions.c (modified) (5 diffs)
• src/dataManip/psMinimize.c (modified) (2 diffs)
• src/dataManip/psStats.c (modified) (11 diffs)
• src/image/psImageStats.c (modified) (20 diffs)
• src/image/psImageStats.h (modified) (1 diff)
• src/imageops/psImageStats.c (modified) (20 diffs)
• src/imageops/psImageStats.h (modified) (1 diff)
• src/math/psMinimize.c (modified) (2 diffs)
• src/math/psPolynomial.c (modified) (5 diffs)
• src/math/psSpline.c (modified) (5 diffs)
• src/math/psStats.c (modified) (11 diffs)
• src/sys/psLogMsg.c (modified) (2 diffs)
• src/sysUtils/psLogMsg.c (modified) (2 diffs)
• test/image/tst_psImageStats02.c (modified) (4 diffs)
• test/sysUtils/tst_psHash02.c (modified) (1 diff)
• test/sysUtils/tst_psHash03.c (modified) (2 diffs)
• test/sysUtils/tst_psLogMsg02.c (modified) (1 diff)

trunk/psLib/src/dataManip/psFunctions.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-14 19:32:42 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-14 19:40:14 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -73 +73 @@
  Gaussian distribution as well.
  
- GUS: There is no way to seed the random generator.
+ NOTE: There is no way to seed the random generator.
  *****************************************************************************/
 psVector *psGaussianDev(float mean,
@@ -94 +94 @@
     }
 
-    // GUS: Should I free r as well?
+    // NOTE: Should I free r as well?
     return(gauss);
 }
@@ -320 +320 @@
     float xSum = 1.0;
 
-    // GUS: Do we want to flag this case?
+    // NOTE: Do we want to flag this case?
     if (myPoly->n == 0) {
         return(1.0);
@@ -646 +646 @@
     double xSum = 1.0;
 
-    // GUS: Do we want to flag this case?
+    // NOTE: Do we want to flag this case?
     if (myPoly->n == 0) {
         return(1.0);

trunk/psLib/src/dataManip/psMinimize.c

@@ -15 +15 @@
 #include "float.h"
 #include <math.h>
-// GUS: rewrite so there is no maximum order for the polynomials.
+// NOTE: rewrite so there is no maximum order for the polynomials.
 #define MAX_POLY_ORDER 10
 #define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
@@ -243 +243 @@
         }
         col[j] = 1.0;
-        // GUS: substitue the LUD rotine
+        // NOTE: substitue the LUD rotine
         //        lubksb(A, N, indx, col);
         for(i=1;i<=N;i++) {

trunk/psLib/src/dataManip/psStats.c

@@ -95 +95 @@
     float binSize = 0.0;                // The histogram bin size
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (n == 0) {
         return(NULL);
     }
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (lower > upper) {
         return(NULL);
@@ -146 +146 @@
     int i;                              // Loop index variable
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (bounds == NULL) {
         //        psAbort(__func__, "psHistogram requested with NULL bounds");
@@ -152 +152 @@
     }
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (bounds->n <= 1) {
         //        psAbort(__func__, "psHistogram requested with NULL bounds");
@@ -215 +215 @@
     int numBins = 0;                            // The total number of bins
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (out == NULL) {
         return(NULL);
     }
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (in == NULL) {
         return(out);
@@ -240 +240 @@
         }
     }
-    // GUS: determine the correct action for a variety of other cases:
+    // NOTE: determine the correct action for a variety of other cases:
     // in vector has 0 elements, and histogram structure has zero bins.
 
@@ -268 +268 @@
                     // bin number requires a bit more work.
                 } else {
-                    // GUS: This is slow.  Put a smarter algorithm here to
+                    // NOTE: This is slow.  Put a smarter algorithm here to
                     // find the correct bin number (bin search, probably)
                     for (j=0;j<(out->bounds->n)-1;j++) {
@@ -841 +841 @@
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
-    // GUS: This is the
+    // NOTE: This is the
 }
 
@@ -950 +950 @@
 
     // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
-    // GUS: This step is dependent on the functions in psMinimize.c being
+    // NOTE: This step is dependent on the functions in psMinimize.c being
     // implemented.  Currently, they are not.
 
@@ -1122 +1122 @@
 
     // 1. Compute the sample median.
-    // GUS: This seems odd.  Verify with IfA that we want to calculate the
+    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
     // median here, not the mean.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
@@ -1198 +1198 @@
                        unsigned int maskVal)
 {
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (in == NULL) {
         return(stats);
@@ -1237 +1237 @@
 
     // ************************************************************************
-    // GUS: The Stdev calculation requires the mean.  Should we assume the
+    // NOTE: The Stdev calculation requires the mean.  Should we assume the
     // mean has already been calculated?  Or should we always calculate it?
     if (stats->options & PS_STAT_SAMPLE_STDEV) {

trunk/psLib/src/image/psImageStats.c

@@ -19 +19 @@
 /// This routine must determine the various statistics for the image.
 /*****************************************************************************
-    GUS: verify that image/mask have the correct types, and sizes.
+    NOTE: verify that image/mask have the correct types, and sizes.
  *****************************************************************************/
 psStats *psImageStats(psStats *stats,
@@ -32 +32 @@
     int j = 0;
 
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if ((stats == NULL) ||
             (in == NULL)) {
         return(NULL);
     }
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if (stats->options == 0) {
         return(stats);
@@ -52 +52 @@
 
     if (mask != NULL) {
-        // GUS: verify that mask data is PS_TYPE_U8.
-        // GUS: figure out mask types
+        // NOTE: verify that mask data is PS_TYPE_U8.
+        // NOTE: figure out mask types
 
         junkMask = psVectorAlloc(mask->numRows * mask->numCols,
                                  mask->type.type);
         junkMask->n = junkMask->nalloc;
-        // GUS: Is there a more efficient way to do this?
+        // NOTE: Is there a more efficient way to do this?
         ptr=0;
         for (i=0;i<mask->numRows;i++) {
@@ -78 +78 @@
     NOTE: We assume that the psHistogram structure out has already been
     allocated and initialized.
-    GUS: verify that image/mask have the, correct types and  sizes.
+    NOTE: verify that image/mask have the, correct types and  sizes.
  *****************************************************************************/
 psHistogram *psImageHistogram(psHistogram *out,
@@ -91 +91 @@
     int j = 0;
 
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if ((out == NULL) ||
             (in == NULL)) {
@@ -100 +100 @@
     junkData->n = junkData->nalloc;
 
-    // GUS: Is there a more efficient way to do this?  memcopy() won't work.
+    // NOTE: Is there a more efficient way to do this?  memcopy() won't work.
     // Can we trick the junkData structure to use the image buffer, then
     // untrick it before we deallocate it (so we won't deallocate that buffer
     // twice?
 
-    // GUS: Make sure you have the numRows/NumCols in the right place.
+    // NOTE: Make sure you have the numRows/NumCols in the right place.
     ptr=0;
     for (i=0;i<in->numRows;i++) {
@@ -114 +114 @@
 
     if (mask != NULL) {
-        // GUS: verify that mask data is PS_TYPE_U8.
-        // GUS: figure out mask types
+        // NOTE: verify that mask data is PS_TYPE_U8.
+        // NOTE: figure out mask types
 
         junkMask = psVectorAlloc(mask->numRows * mask->numCols,
@@ -121 +121 @@
         junkMask->n = junkMask->nalloc;
 
-        // GUS: Is there a more efficient way to do this?
+        // NOTE: Is there a more efficient way to do this?
         ptr=0;
         for (i=0;i<mask->numRows;i++) {
@@ -137 +137 @@
     return(out);
 }
+
+float *p_psCalcScaleFactorsFit(int n)
+{
+    int i = 0;
+    float tmp = 0.0;
+    float *scalingFactors = (float *) psAlloc(n * sizeof(float));
+
+    for (i=0;i<n;i++) {
+        //     ((2.0 * (float) i) / ((float) (n-1))) - 1.0;
+        //        tmp = (float) (i + 1);
+        tmp = (float) (n - i);
+        tmp = (M_PI * (tmp - 0.5)) / ((float) n);
+        scalingFactors[i] = cos(tmp);
+    }
+
+    return(scalingFactors);
+}
+
+float *p_psCalcScaleFactorsEval(int n)
+{
+    return p_psCalcScaleFactorsFit(n);
+
+    int i = 0;
+    //    float tmp = 0.0;
+    float *scalingFactors = (float *) psAlloc(n * sizeof(float));
+    for (i=0;i<n;i++) {
+        //          scalingFactors[i] = ((2.0 * (float) i) / ((float) (n-1))) - 1.0;
+        tmp = (float) (i + 1);
+        tmp = (M_PI * (tmp - 0.5)) / ((float) n);
+        scalingFactors[i] = cos(tmp);
+    }
+    return(scalingFactors);
+}
+
+psPolynomial1D **p_psCreateChebyshevPolys(int maxChebyPoly)
+{
+    psPolynomial1D **chebPolys = NULL;
+    int i = 0;
+    int j = 0;
+
+    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
+                                            sizeof(psPolynomial1D *));
+    for (i=0;i<maxChebyPoly;i++) {
+        chebPolys[i] = psPolynomial1DAlloc(i+1);
+    }
+
+    // Create the Chebyshev polynomials.  Polynomial i has i-th order.
+    chebPolys[0]->coeff[0] = 1;
+    chebPolys[1]->coeff[1] = 1;
+    for (i=2;i<maxChebyPoly;i++) {
+        for (j=0;j<chebPolys[i-1]->n;j++) {
+            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
+        }
+        for (j=0;j<chebPolys[i-2]->n;j++) {
+            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
+        }
+    }
+
+    return(chebPolys);
+}
+
 
 /*****************************************************************************
@@ -146 +207 @@
     Output:
     Internal Data Structures:
+ 
  chebPolys[i][j] 
-        sums[i][j]
+ 
+ sums[i][j]: This will contain the sum of 
+  input->data.F32[x][y] *
+                psEvalPolynomial1D((float) x, chebPolys[i]) *
+                psEvalPolynomial1D((float) y, chebPolys[j]);
+      over all pixels (x,y) in the image.
+ 
+ 
  *****************************************************************************/
 psPolynomial2D *
@@ -153 +222 @@
                      psPolynomial2D *coeffs)
 {
+    int k = 0;
+    int l = 0;
     int x = 0;
     int y = 0;
@@ -160 +231 @@
     psPolynomial1D **chebPolys = NULL;
     int maxChebyPoly = 0;
-
+    float *cScalingFactors = NULL;
+    float *rScalingFactors = NULL;
+    float tmp = 0.0;
 
     // Create the sums[][] data structure.  This will hold the LHS of equation
     // 29 in the ADD: sums[k][l] = SUM { image(x,y) * Tk(x) * Tl(y) }
-
     sums = (float **) psAlloc(coeffs->nX * sizeof(float *));
     for (i=0;i<coeffs->nX;i++) {
         sums[i] = (float *) psAlloc(coeffs->nY * sizeof(float));
     }
-    for (i=0;i<coeffs->nX;i++) {
-        for (j=0;j<coeffs->nY;j++) {
-            sums[i][j] = 0.0;
-        }
-    }
+
+    // We scale the pixel positions to values between -1.0 and 1.0
+    cScalingFactors = p_psCalcScaleFactorsFit(input->numRows);
+    rScalingFactors = p_psCalcScaleFactorsFit(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -180 +251 @@
         maxChebyPoly = coeffs->nY;
     }
-
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
-                                            sizeof(psPolynomial1D *));
-    for (i=0;i<maxChebyPoly;i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i);
-    }
-
-    // Create the Chebyshev polynomials
-    chebPolys[1]->coeff[0] = 1;
-    for (i=2;i<maxChebyPoly;i++) {
-        for (j=0;j<chebPolys[i-1]->n;j++) {
-            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
-        }
-        for (j=0;j<chebPolys[i-2]->n;j++) {
-            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
-        }
-    }
-
-    for (x=0;x<input->numRows;x++) {
-        for (y=0;y<input->numCols;y++) {
-            for (i=0;i<coeffs->nX;i++) {
-                for (j=0;j<coeffs->nY;j++) {
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
+
+    // Sanity check for the Chebyshevs
+    for (i=0;i<coeffs->nX;i++) {
+        for (j=0;j<coeffs->nY;j++) {
+            tmp = 0.0;
+            for (x=0;x<input->numRows;x++) {
+                tmp+= psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                      psEvalPolynomial1D(rScalingFactors[x], chebPolys[j]);
+
+            }
+        }
+    }
+
+    for (i=0;i<coeffs->nX;i++) {
+        for (j=0;j<coeffs->nY;j++) {
+            sums[i][j] = 0.0;
+            for (x=0;x<input->numRows;x++) {
+                for (y=0;y<input->numCols;y++) {
                     sums[i][j]+= input->data.F32[x][y] *
-                                 psEvalPolynomial1D((float) x, chebPolys[i]) *
-                                 psEvalPolynomial1D((float) y, chebPolys[j]);
-
+                                 psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                                 psEvalPolynomial1D(cScalingFactors[y], chebPolys[j]);
                 }
             }
@@ -211 +278 @@
     }
 
-    // GUS: Check the math on this at a later date.
-    coeffs->coeff[0][0] = sums[0][0] / (coeffs->nX * coeffs->nY);
-    for (i=0;i<coeffs->nX;i++) {
-        coeffs->coeff[i][0] = 0.5 *
-                              ((sums[i][0] * 4.0) / (coeffs->nX * coeffs->nY)) -
-                              (2.0 * coeffs->coeff[0][0]);
+    // NOTE: Check the math on this at a later date.
+    coeffs->coeff[0][0] = sums[0][0] / ((float) (coeffs->nX * coeffs->nY));
+    for (i=0;i<coeffs->nX;i++) {
+        coeffs->coeff[i][0] = ((sums[i][0] * 2.0) /
+                               ((float) (coeffs->nX * coeffs->nY))) -
+                              coeffs->coeff[0][0];
     }
     for (j=0;j<coeffs->nY;j++) {
-        coeffs->coeff[0][j] = 0.5 *
-                              ((sums[0][j] * 4.0) / (coeffs->nX * coeffs->nY)) -
-                              (2.0 * coeffs->coeff[0][0]);
+        coeffs->coeff[0][j] =
+            ((sums[0][j] * 2.0) /
+             ((float) (coeffs->nX * coeffs->nY))) -
+            coeffs->coeff[0][0];
     }
     for (i=1;i<coeffs->nX;i++) {
         for (j=1;j<coeffs->nY;j++) {
             coeffs->coeff[i][j] =
-                (sums[i][0] * 4.0) / (coeffs->nX * coeffs->nY) -
-                coeffs->coeff[0][0] -
-                coeffs->coeff[i][0] -
-                coeffs->coeff[0][j];
-        }
-    }
+                ((sums[i][0] * 4.0) / ((float) (coeffs->nX * coeffs->nY))) -
+                (coeffs->coeff[0][0] +
+                 coeffs->coeff[i][0] +
+                 coeffs->coeff[0][j]);
+        }
+    }
+
+    for (k=0;k<coeffs->nX;k++) {
+        for (l=0;l<coeffs->nY;l++) {
+            tmp = 0.0;
+            for (i=0;i<coeffs->nX;i++) {
+                for (j=0;j<coeffs->nY;j++) {
+                    for (x=0;x<input->numRows;x++) {
+                        for (y=0;y<input->numCols;y++) {
+                            tmp+= (coeffs->coeff[i][j] *
+                                   psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                                   psEvalPolynomial1D(rScalingFactors[y], chebPolys[j]) *
+                                   psEvalPolynomial1D(rScalingFactors[x], chebPolys[k]) *
+                                   psEvalPolynomial1D(rScalingFactors[y], chebPolys[l]));
+
+                        }
+                    }
+                }
+            }
+        }
+    }
+
 
     // Free the Chebyshev polynomials that were created in this routine.
@@ -244 +333 @@
     }
     psFree(sums);
+    psFree(cScalingFactors);
+    psFree(rScalingFactors);
 
     return(coeffs);
@@ -255 +346 @@
                       const psPolynomial2D *coeffs)
 {
-    int i = 0;
-    int j = 0;
     int x = 0;
     int y = 0;
+    int i = 0;
+    int j = 0;
     float **sums = NULL;
     psPolynomial1D **chebPolys = NULL;
     int maxChebyPoly = 0;
+    float *cScalingFactors = NULL;
+    float *rScalingFactors = NULL;
     float polySum = 0.0;
 
@@ -275 +368 @@
         }
     }
+
+    // We scale the pixel positions to values between -1.0 and 1.0
+    cScalingFactors = p_psCalcScaleFactorsEval(input->numRows);
+    rScalingFactors = p_psCalcScaleFactorsEval(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -282 +379 @@
     }
 
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
-                                            sizeof(psPolynomial1D *));
-    for (i=0;i<maxChebyPoly;i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i);
-    }
-
-    // Create the Chebyshev polynomials
-    chebPolys[1]->coeff[0] = 1;
-    for (i=2;i<maxChebyPoly;i++) {
-        for (j=0;j<chebPolys[i-1]->n;j++) {
-            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
-        }
-        for (j=0;j<chebPolys[i-2]->n;j++) {
-            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
-        }
-    }
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
 
     for (x=0;x<input->numRows;x++) {
@@ -304 +386 @@
             for (i=0;i<coeffs->nX;i++) {
                 for (j=0;j<coeffs->nY;j++) {
-                    polySum+= psEvalPolynomial1D((float)x, chebPolys[i]) *
-                              psEvalPolynomial1D((float)y, chebPolys[j]) *
+                    polySum+= psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                              psEvalPolynomial1D(cScalingFactors[y], chebPolys[j]) *
                               coeffs->coeff[i][j];
+
                 }
             }
@@ -324 +407 @@
     }
     psFree(sums);
+    psFree(cScalingFactors);
+    psFree(rScalingFactors);
 
     return(0);

trunk/psLib/src/image/psImageStats.h

@@ +1 @@
+/** @file psImageStats.h
+ *  \brief Routines for calculating statistics on images.
+ *  \ingroup Stats
+ *  @ingroup Stats
+ *
+ *  This file will hold the prototypes for procedures which calculate
+ *  statistic on images, histograms on images, and fit/evaluate Chebyshev
+ *  polynomials to images.
+ *
+ *  @author George Gusciora, MHPCC
+ *
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-14 19:40:14 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
 #if !defined(PS_IMAGE_STATS_H)
 #define PS_IMAGE_STATS_H
+
+
 #include "psVector.h"
 #include "psImage.h"

trunk/psLib/src/imageops/psImageStats.c

@@ -19 +19 @@
 /// This routine must determine the various statistics for the image.
 /*****************************************************************************
-    GUS: verify that image/mask have the correct types, and sizes.
+    NOTE: verify that image/mask have the correct types, and sizes.
  *****************************************************************************/
 psStats *psImageStats(psStats *stats,
@@ -32 +32 @@
     int j = 0;
 
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if ((stats == NULL) ||
             (in == NULL)) {
         return(NULL);
     }
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if (stats->options == 0) {
         return(stats);
@@ -52 +52 @@
 
     if (mask != NULL) {
-        // GUS: verify that mask data is PS_TYPE_U8.
-        // GUS: figure out mask types
+        // NOTE: verify that mask data is PS_TYPE_U8.
+        // NOTE: figure out mask types
 
         junkMask = psVectorAlloc(mask->numRows * mask->numCols,
                                  mask->type.type);
         junkMask->n = junkMask->nalloc;
-        // GUS: Is there a more efficient way to do this?
+        // NOTE: Is there a more efficient way to do this?
         ptr=0;
         for (i=0;i<mask->numRows;i++) {
@@ -78 +78 @@
     NOTE: We assume that the psHistogram structure out has already been
     allocated and initialized.
-    GUS: verify that image/mask have the, correct types and  sizes.
+    NOTE: verify that image/mask have the, correct types and  sizes.
  *****************************************************************************/
 psHistogram *psImageHistogram(psHistogram *out,
@@ -91 +91 @@
     int j = 0;
 
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if ((out == NULL) ||
             (in == NULL)) {
@@ -100 +100 @@
     junkData->n = junkData->nalloc;
 
-    // GUS: Is there a more efficient way to do this?  memcopy() won't work.
+    // NOTE: Is there a more efficient way to do this?  memcopy() won't work.
     // Can we trick the junkData structure to use the image buffer, then
     // untrick it before we deallocate it (so we won't deallocate that buffer
     // twice?
 
-    // GUS: Make sure you have the numRows/NumCols in the right place.
+    // NOTE: Make sure you have the numRows/NumCols in the right place.
     ptr=0;
     for (i=0;i<in->numRows;i++) {
@@ -114 +114 @@
 
     if (mask != NULL) {
-        // GUS: verify that mask data is PS_TYPE_U8.
-        // GUS: figure out mask types
+        // NOTE: verify that mask data is PS_TYPE_U8.
+        // NOTE: figure out mask types
 
         junkMask = psVectorAlloc(mask->numRows * mask->numCols,
@@ -121 +121 @@
         junkMask->n = junkMask->nalloc;
 
-        // GUS: Is there a more efficient way to do this?
+        // NOTE: Is there a more efficient way to do this?
         ptr=0;
         for (i=0;i<mask->numRows;i++) {
@@ -137 +137 @@
     return(out);
 }
+
+float *p_psCalcScaleFactorsFit(int n)
+{
+    int i = 0;
+    float tmp = 0.0;
+    float *scalingFactors = (float *) psAlloc(n * sizeof(float));
+
+    for (i=0;i<n;i++) {
+        //     ((2.0 * (float) i) / ((float) (n-1))) - 1.0;
+        //        tmp = (float) (i + 1);
+        tmp = (float) (n - i);
+        tmp = (M_PI * (tmp - 0.5)) / ((float) n);
+        scalingFactors[i] = cos(tmp);
+    }
+
+    return(scalingFactors);
+}
+
+float *p_psCalcScaleFactorsEval(int n)
+{
+    return p_psCalcScaleFactorsFit(n);
+
+    int i = 0;
+    //    float tmp = 0.0;
+    float *scalingFactors = (float *) psAlloc(n * sizeof(float));
+    for (i=0;i<n;i++) {
+        //          scalingFactors[i] = ((2.0 * (float) i) / ((float) (n-1))) - 1.0;
+        tmp = (float) (i + 1);
+        tmp = (M_PI * (tmp - 0.5)) / ((float) n);
+        scalingFactors[i] = cos(tmp);
+    }
+    return(scalingFactors);
+}
+
+psPolynomial1D **p_psCreateChebyshevPolys(int maxChebyPoly)
+{
+    psPolynomial1D **chebPolys = NULL;
+    int i = 0;
+    int j = 0;
+
+    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
+                                            sizeof(psPolynomial1D *));
+    for (i=0;i<maxChebyPoly;i++) {
+        chebPolys[i] = psPolynomial1DAlloc(i+1);
+    }
+
+    // Create the Chebyshev polynomials.  Polynomial i has i-th order.
+    chebPolys[0]->coeff[0] = 1;
+    chebPolys[1]->coeff[1] = 1;
+    for (i=2;i<maxChebyPoly;i++) {
+        for (j=0;j<chebPolys[i-1]->n;j++) {
+            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
+        }
+        for (j=0;j<chebPolys[i-2]->n;j++) {
+            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
+        }
+    }
+
+    return(chebPolys);
+}
+
 
 /*****************************************************************************
@@ -146 +207 @@
     Output:
     Internal Data Structures:
+ 
  chebPolys[i][j] 
-        sums[i][j]
+ 
+ sums[i][j]: This will contain the sum of 
+  input->data.F32[x][y] *
+                psEvalPolynomial1D((float) x, chebPolys[i]) *
+                psEvalPolynomial1D((float) y, chebPolys[j]);
+      over all pixels (x,y) in the image.
+ 
+ 
  *****************************************************************************/
 psPolynomial2D *
@@ -153 +222 @@
                      psPolynomial2D *coeffs)
 {
+    int k = 0;
+    int l = 0;
     int x = 0;
     int y = 0;
@@ -160 +231 @@
     psPolynomial1D **chebPolys = NULL;
     int maxChebyPoly = 0;
-
+    float *cScalingFactors = NULL;
+    float *rScalingFactors = NULL;
+    float tmp = 0.0;
 
     // Create the sums[][] data structure.  This will hold the LHS of equation
     // 29 in the ADD: sums[k][l] = SUM { image(x,y) * Tk(x) * Tl(y) }
-
     sums = (float **) psAlloc(coeffs->nX * sizeof(float *));
     for (i=0;i<coeffs->nX;i++) {
         sums[i] = (float *) psAlloc(coeffs->nY * sizeof(float));
     }
-    for (i=0;i<coeffs->nX;i++) {
-        for (j=0;j<coeffs->nY;j++) {
-            sums[i][j] = 0.0;
-        }
-    }
+
+    // We scale the pixel positions to values between -1.0 and 1.0
+    cScalingFactors = p_psCalcScaleFactorsFit(input->numRows);
+    rScalingFactors = p_psCalcScaleFactorsFit(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -180 +251 @@
         maxChebyPoly = coeffs->nY;
     }
-
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
-                                            sizeof(psPolynomial1D *));
-    for (i=0;i<maxChebyPoly;i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i);
-    }
-
-    // Create the Chebyshev polynomials
-    chebPolys[1]->coeff[0] = 1;
-    for (i=2;i<maxChebyPoly;i++) {
-        for (j=0;j<chebPolys[i-1]->n;j++) {
-            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
-        }
-        for (j=0;j<chebPolys[i-2]->n;j++) {
-            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
-        }
-    }
-
-    for (x=0;x<input->numRows;x++) {
-        for (y=0;y<input->numCols;y++) {
-            for (i=0;i<coeffs->nX;i++) {
-                for (j=0;j<coeffs->nY;j++) {
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
+
+    // Sanity check for the Chebyshevs
+    for (i=0;i<coeffs->nX;i++) {
+        for (j=0;j<coeffs->nY;j++) {
+            tmp = 0.0;
+            for (x=0;x<input->numRows;x++) {
+                tmp+= psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                      psEvalPolynomial1D(rScalingFactors[x], chebPolys[j]);
+
+            }
+        }
+    }
+
+    for (i=0;i<coeffs->nX;i++) {
+        for (j=0;j<coeffs->nY;j++) {
+            sums[i][j] = 0.0;
+            for (x=0;x<input->numRows;x++) {
+                for (y=0;y<input->numCols;y++) {
                     sums[i][j]+= input->data.F32[x][y] *
-                                 psEvalPolynomial1D((float) x, chebPolys[i]) *
-                                 psEvalPolynomial1D((float) y, chebPolys[j]);
-
+                                 psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                                 psEvalPolynomial1D(cScalingFactors[y], chebPolys[j]);
                 }
             }
@@ -211 +278 @@
     }
 
-    // GUS: Check the math on this at a later date.
-    coeffs->coeff[0][0] = sums[0][0] / (coeffs->nX * coeffs->nY);
-    for (i=0;i<coeffs->nX;i++) {
-        coeffs->coeff[i][0] = 0.5 *
-                              ((sums[i][0] * 4.0) / (coeffs->nX * coeffs->nY)) -
-                              (2.0 * coeffs->coeff[0][0]);
+    // NOTE: Check the math on this at a later date.
+    coeffs->coeff[0][0] = sums[0][0] / ((float) (coeffs->nX * coeffs->nY));
+    for (i=0;i<coeffs->nX;i++) {
+        coeffs->coeff[i][0] = ((sums[i][0] * 2.0) /
+                               ((float) (coeffs->nX * coeffs->nY))) -
+                              coeffs->coeff[0][0];
     }
     for (j=0;j<coeffs->nY;j++) {
-        coeffs->coeff[0][j] = 0.5 *
-                              ((sums[0][j] * 4.0) / (coeffs->nX * coeffs->nY)) -
-                              (2.0 * coeffs->coeff[0][0]);
+        coeffs->coeff[0][j] =
+            ((sums[0][j] * 2.0) /
+             ((float) (coeffs->nX * coeffs->nY))) -
+            coeffs->coeff[0][0];
     }
     for (i=1;i<coeffs->nX;i++) {
         for (j=1;j<coeffs->nY;j++) {
             coeffs->coeff[i][j] =
-                (sums[i][0] * 4.0) / (coeffs->nX * coeffs->nY) -
-                coeffs->coeff[0][0] -
-                coeffs->coeff[i][0] -
-                coeffs->coeff[0][j];
-        }
-    }
+                ((sums[i][0] * 4.0) / ((float) (coeffs->nX * coeffs->nY))) -
+                (coeffs->coeff[0][0] +
+                 coeffs->coeff[i][0] +
+                 coeffs->coeff[0][j]);
+        }
+    }
+
+    for (k=0;k<coeffs->nX;k++) {
+        for (l=0;l<coeffs->nY;l++) {
+            tmp = 0.0;
+            for (i=0;i<coeffs->nX;i++) {
+                for (j=0;j<coeffs->nY;j++) {
+                    for (x=0;x<input->numRows;x++) {
+                        for (y=0;y<input->numCols;y++) {
+                            tmp+= (coeffs->coeff[i][j] *
+                                   psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                                   psEvalPolynomial1D(rScalingFactors[y], chebPolys[j]) *
+                                   psEvalPolynomial1D(rScalingFactors[x], chebPolys[k]) *
+                                   psEvalPolynomial1D(rScalingFactors[y], chebPolys[l]));
+
+                        }
+                    }
+                }
+            }
+        }
+    }
+
 
     // Free the Chebyshev polynomials that were created in this routine.
@@ -244 +333 @@
     }
     psFree(sums);
+    psFree(cScalingFactors);
+    psFree(rScalingFactors);
 
     return(coeffs);
@@ -255 +346 @@
                       const psPolynomial2D *coeffs)
 {
-    int i = 0;
-    int j = 0;
     int x = 0;
     int y = 0;
+    int i = 0;
+    int j = 0;
     float **sums = NULL;
     psPolynomial1D **chebPolys = NULL;
     int maxChebyPoly = 0;
+    float *cScalingFactors = NULL;
+    float *rScalingFactors = NULL;
     float polySum = 0.0;
 
@@ -275 +368 @@
         }
     }
+
+    // We scale the pixel positions to values between -1.0 and 1.0
+    cScalingFactors = p_psCalcScaleFactorsEval(input->numRows);
+    rScalingFactors = p_psCalcScaleFactorsEval(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -282 +379 @@
     }
 
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
-                                            sizeof(psPolynomial1D *));
-    for (i=0;i<maxChebyPoly;i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i);
-    }
-
-    // Create the Chebyshev polynomials
-    chebPolys[1]->coeff[0] = 1;
-    for (i=2;i<maxChebyPoly;i++) {
-        for (j=0;j<chebPolys[i-1]->n;j++) {
-            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
-        }
-        for (j=0;j<chebPolys[i-2]->n;j++) {
-            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
-        }
-    }
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
 
     for (x=0;x<input->numRows;x++) {
@@ -304 +386 @@
             for (i=0;i<coeffs->nX;i++) {
                 for (j=0;j<coeffs->nY;j++) {
-                    polySum+= psEvalPolynomial1D((float)x, chebPolys[i]) *
-                              psEvalPolynomial1D((float)y, chebPolys[j]) *
+                    polySum+= psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                              psEvalPolynomial1D(cScalingFactors[y], chebPolys[j]) *
                               coeffs->coeff[i][j];
+
                 }
             }
@@ -324 +407 @@
     }
     psFree(sums);
+    psFree(cScalingFactors);
+    psFree(rScalingFactors);
 
     return(0);

trunk/psLib/src/imageops/psImageStats.h

@@ +1 @@
+/** @file psImageStats.h
+ *  \brief Routines for calculating statistics on images.
+ *  \ingroup Stats
+ *  @ingroup Stats
+ *
+ *  This file will hold the prototypes for procedures which calculate
+ *  statistic on images, histograms on images, and fit/evaluate Chebyshev
+ *  polynomials to images.
+ *
+ *  @author George Gusciora, MHPCC
+ *
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-14 19:40:14 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
 #if !defined(PS_IMAGE_STATS_H)
 #define PS_IMAGE_STATS_H
+
+
 #include "psVector.h"
 #include "psImage.h"

trunk/psLib/src/math/psMinimize.c

@@ -15 +15 @@
 #include "float.h"
 #include <math.h>
-// GUS: rewrite so there is no maximum order for the polynomials.
+// NOTE: rewrite so there is no maximum order for the polynomials.
 #define MAX_POLY_ORDER 10
 #define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
@@ -243 +243 @@
         }
         col[j] = 1.0;
-        // GUS: substitue the LUD rotine
+        // NOTE: substitue the LUD rotine
         //        lubksb(A, N, indx, col);
         for(i=1;i<=N;i++) {

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-14 19:32:42 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-14 19:40:14 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -73 +73 @@
  Gaussian distribution as well.
  
- GUS: There is no way to seed the random generator.
+ NOTE: There is no way to seed the random generator.
  *****************************************************************************/
 psVector *psGaussianDev(float mean,
@@ -94 +94 @@
     }
 
-    // GUS: Should I free r as well?
+    // NOTE: Should I free r as well?
     return(gauss);
 }
@@ -320 +320 @@
     float xSum = 1.0;
 
-    // GUS: Do we want to flag this case?
+    // NOTE: Do we want to flag this case?
     if (myPoly->n == 0) {
         return(1.0);
@@ -646 +646 @@
     double xSum = 1.0;
 
-    // GUS: Do we want to flag this case?
+    // NOTE: Do we want to flag this case?
     if (myPoly->n == 0) {
         return(1.0);

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-14 19:32:42 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-14 19:40:14 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -73 +73 @@
  Gaussian distribution as well.
  
- GUS: There is no way to seed the random generator.
+ NOTE: There is no way to seed the random generator.
  *****************************************************************************/
 psVector *psGaussianDev(float mean,
@@ -94 +94 @@
     }
 
-    // GUS: Should I free r as well?
+    // NOTE: Should I free r as well?
     return(gauss);
 }
@@ -320 +320 @@
     float xSum = 1.0;
 
-    // GUS: Do we want to flag this case?
+    // NOTE: Do we want to flag this case?
     if (myPoly->n == 0) {
         return(1.0);
@@ -646 +646 @@
     double xSum = 1.0;
 
-    // GUS: Do we want to flag this case?
+    // NOTE: Do we want to flag this case?
     if (myPoly->n == 0) {
         return(1.0);

trunk/psLib/src/math/psStats.c

@@ -95 +95 @@
     float binSize = 0.0;                // The histogram bin size
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (n == 0) {
         return(NULL);
     }
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (lower > upper) {
         return(NULL);
@@ -146 +146 @@
     int i;                              // Loop index variable
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (bounds == NULL) {
         //        psAbort(__func__, "psHistogram requested with NULL bounds");
@@ -152 +152 @@
     }
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (bounds->n <= 1) {
         //        psAbort(__func__, "psHistogram requested with NULL bounds");
@@ -215 +215 @@
     int numBins = 0;                            // The total number of bins
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (out == NULL) {
         return(NULL);
     }
 
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (in == NULL) {
         return(out);
@@ -240 +240 @@
         }
     }
-    // GUS: determine the correct action for a variety of other cases:
+    // NOTE: determine the correct action for a variety of other cases:
     // in vector has 0 elements, and histogram structure has zero bins.
 
@@ -268 +268 @@
                     // bin number requires a bit more work.
                 } else {
-                    // GUS: This is slow.  Put a smarter algorithm here to
+                    // NOTE: This is slow.  Put a smarter algorithm here to
                     // find the correct bin number (bin search, probably)
                     for (j=0;j<(out->bounds->n)-1;j++) {
@@ -841 +841 @@
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
-    // GUS: This is the
+    // NOTE: This is the
 }
 
@@ -950 +950 @@
 
     // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
-    // GUS: This step is dependent on the functions in psMinimize.c being
+    // NOTE: This step is dependent on the functions in psMinimize.c being
     // implemented.  Currently, they are not.
 
@@ -1122 +1122 @@
 
     // 1. Compute the sample median.
-    // GUS: This seems odd.  Verify with IfA that we want to calculate the
+    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
     // median here, not the mean.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
@@ -1198 +1198 @@
                        unsigned int maskVal)
 {
-    // GUS: Verify that this is the correct action.
+    // NOTE: Verify that this is the correct action.
     if (in == NULL) {
         return(stats);
@@ -1237 +1237 @@
 
     // ************************************************************************
-    // GUS: The Stdev calculation requires the mean.  Should we assume the
+    // NOTE: The Stdev calculation requires the mean.  Should we assume the
     // mean has already been calculated?  Or should we always calculate it?
     if (stats->options & PS_STAT_SAMPLE_STDEV) {

trunk/psLib/src/sys/psLogMsg.c

@@ -11 +11 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-12 05:50:01 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-14 19:40:15 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -105 +105 @@
 
     default:
-        // GUS: Should you log this error properly?
+        // NOTE: Should you log this error properly?
         fprintf(stderr,"Unknown p_psGlobalLogDest: %d (ignored)\n", dest);
         break;

trunk/psLib/src/sysUtils/psLogMsg.c

@@ -11 +11 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-12 05:50:01 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-14 19:40:15 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -105 +105 @@
 
     default:
-        // GUS: Should you log this error properly?
+        // NOTE: Should you log this error properly?
         fprintf(stderr,"Unknown p_psGlobalLogDest: %d (ignored)\n", dest);
         break;

trunk/psLib/test/image/tst_psImageStats02.c

@@ -4 +4 @@
     thouroughly tested elsewhere, we will only test psImageStats() with
     the PS_STAT_SAMPLE_MEAN here.
+ 
+    NOTE: After you debug, set CHEBY_X_DIM != CHEBY_Y_DIM.
+    Also, uses non-square images.
  *****************************************************************************/
 #include <stdio.h>
@@ -11 +14 @@
 #include "psImage.h"
 #include "psImageStats.h"
-#define NUM_BINS 20
-#define IMAGE_SIZE 10
-#define CHEBY_X_DIM 5
+#define IMAGE_SIZE 3
+#define CHEBY_X_DIM 3
 #define CHEBY_Y_DIM 3
 
@@ -38 +40 @@
         for (j=0;j<IMAGE_SIZE;j++) {
             tmpImage->data.F32[i][j] = (float) (i + j);
+            tmpImage->data.F32[i][j] = 4.0;
             outImage->data.F32[i][j] = 0.0;
-        }
-    }
-    for (i=0;i<IMAGE_SIZE;i++) {
-        for (j=0;j<IMAGE_SIZE;j++) {
-            printf("(%f, %f)\n", tmpImage->data.F32[i][j],
-                   outImage->data.F32[i][j]);
         }
     }
@@ -72 +69 @@
     for (i=0;i<IMAGE_SIZE;i++) {
         for (j=0;j<IMAGE_SIZE;j++) {
-            printf("(%f, %f)\n", tmpImage->data.F32[i][j],
+            printf("pixel[%d][%d] is (%f, %f)\n", i, j,
+                   tmpImage->data.F32[i][j],
                    outImage->data.F32[i][j]);
         }

trunk/psLib/test/sysUtils/tst_psHash02.c

@@ -3 +3 @@
     and retrieved correctly.
  
-    GUS: Add code to test whether duplicates are handled correctly (use a
+    NOTE: Add code to test whether duplicates are handled correctly (use a
     small hash table and lots of keys).
  *****************************************************************************/

trunk/psLib/test/sysUtils/tst_psHash03.c

@@ -2 +2 @@
     This code will test whether hash tables entries can be removed correctly.
  
-    GUS: Add code to test whether duplicates are handled correctly.
+    NOTE: Add code to test whether duplicates are handled correctly.
  *****************************************************************************/
 #include <stdio.h>
@@ -102 +102 @@
     while (myKeys[i] != NULL) {
         id = psHashRemove(myHashTable, myKeys[i], (void (*)(void *))IdFree);
-        // GUS: The psHashRemove() procedure removes the entry from the hash
+        // The psHashRemove() procedure removes the entry from the hash
         // table, but does not delete the data.  The following is necessary
         // to delete the data as well (which is returned from the call).

trunk/psLib/test/sysUtils/tst_psLogMsg02.c

@@ -67 +67 @@
                 testStatus);
 
+
+
+
+    printPositiveTestHeader(stdout,
+                            "psLogMsg functions",
+                            "Output Format");
+
+    psLogMsg("Under 15 chars", 0, "Hello World!\n");
+    psLogMsg("This string is more than 15 chars", 0, "Hello World!\n");
+    psLogMsg(__func__, 0, "Line #1\n");
+    psLogMsg(__func__, 0, "Line #2\n");
+    psLogMsg(__func__, 0, "Line #3");
+    psLogMsg(__func__, 0, "Line #4");
+
+    printFooter(stdout,
+                "psLogMsg functions",
+                "Output Format",
+                testStatus);
+
+
     if (psMemCheckLeaks(currentId,NULL,NULL) != 0) {
         psAbort(__func__,"Memory Leaks!");