Changeset 4958

Timestamp:

Sep 7, 2005, 11:35:50 AM (21 years ago)

Author:

gusciora

Message:

Split of psFunctions.c into psSpline.c and psPolynomial.c.

Location:

trunk/psLib/src/math

Files:

• Makefile.am (modified) (2 diffs)
• psFunctions.c (modified) (2 diffs)
• psMinimize.c (modified) (15 diffs)
• psMinimize.h (modified) (3 diffs)
• psPolynomial.c (added)
• psPolynomial.h (added)
• psSpline.c (added)
• psSpline.h (added)
• psStats.c (modified) (15 diffs)

trunk/psLib/src/math/Makefile.am

@@ -9 +9 @@
         psMinimize.c \
         psRandom.c \
-        psFunctions.c \
+        psPolynomial.c \
+        psSpline.c \
         psStats.c \
         psUnaryOp.c
@@ -23 +24 @@
         psMinimize.h \
         psRandom.h \
-        psFunctions.h \
+        psPolynomial.h \
+        psSpline.h \
         psStats.h \
         psUnaryOp.h

trunk/psLib/src/math/psFunctions.c

@@ -7 +7 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-08-31 22:28:35 $
+*  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-07 21:35:50 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -1307 +1307 @@
 /*****************************************************************************
 p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
+ 
+XXX: Assert that the psVector and psScalar have the same type.
  *****************************************************************************/
 psS32 p_psVectorBinDisect(psVector *bins,

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.133 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-02 21:32:06 $
+ *  @version $Revision: 1.134 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-07 21:35:50 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -29 +29 @@
 #include "psImageStructManip.h"
 #include "psBinaryOp.h"
+#include "psLogMsg.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
@@ -112 +113 @@
 
 /******************************************************************************
-p_psBuildSums1D(x, polyOrder, sums): this routine calculates the powers of
-input parameter "x" between 0 and input parameter polyOrder.  The result is
-returned as a psVector sums.
+buildSums1D(x, polyOrder, sums): this routine calculates the powers of input
+parameter "x" between 0 and input parameter polyOrder.  The result is returned
+as a psVector sums.
  
 XXX: Use a static vector.
+ 
+XXX: should the argument be polyOrder or numTerms?
  *****************************************************************************/
 static void buildSums1D(psF64 x,
@@ -124 +127 @@
     psS32 i = 0;
     psF64 xSum = 0.0;
+    psS32 numTerms = polyOrder + 1;
 
     if (sums == NULL) {
-        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
-    }
-    if (polyOrder > sums->n) {
-        sums = psVectorRealloc(sums, polyOrder);
+        sums = psVectorAlloc(numTerms, PS_TYPE_F64);
+    } else if (numTerms > sums->n) {
+        sums = psVectorRealloc(sums, numTerms);
     }
 
     xSum = 1.0;
-    for (i = 0; i <= polyOrder; i++) {
+    for (i = 0; i < numTerms; i++) {
         sums->data.F64[i] = xSum;
         xSum *= x;
@@ -708 +711 @@
 those coords.
  
-// XXX EAM this is my re-implementation of MinLM
+XXX EAM this is my re-implementation of MinLM
  
 XXX: This must work for both F32 and F64.  F32 is currently implemented.
@@ -1306 +1309 @@
  
 XXX: type F32 is done via vector conversion only.
+ 
+XXX: Get rid of this.  Use new argument list.
  *****************************************************************************/
 psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* poly,
@@ -2143 +2148 @@
 
 /******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
 EAM Code:
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
  *****************************************************************************/
 
-// XXX EAM : my alternate EAMBuildSums1D
 static psVector *EAMBuildSums1D(
     psVector* sums,
@@ -2155 +2165 @@
     psF64 xSum = 0.0;
 
+    //
+    // XXX: Why do we multiply by 2 here?  Better to do it outside and have the
+    // definition of this function remain sensible.
+    //
     nSum = 2*nTerm;
     if (sums == NULL) {
         sums = psVectorAlloc(nSum, PS_TYPE_F64);
-    }
-    if (nSum > sums->n) {
+    } else if (nSum > sums->n) {
         sums = psVectorRealloc(sums, nSum);
     }
@@ -2177 +2190 @@
     const psVector *x,
     const psVector *y,
-    const psVector *yErr
-)
+    const psVector *yErr)
 {
     // I think this is 1 dimension down
@@ -2277 +2289 @@
     const psPolynomial2D *myPoly,
     const psVector *x,
-    const psVector *y
-)
+    const psVector *y)
 {
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
@@ -2312 +2323 @@
     psF64 y,
     psS32 nXterm,
-    psS32 nYterm
-)
+    psS32 nYterm)
 {
     psS32 nXsum = 0;
@@ -2349 +2359 @@
     const psVector* y,
     const psVector* z,
-    const psVector* zErr
-)
+    const psVector* zErr)
 {
     // I think this is 1 dimension down
@@ -2439 +2448 @@
     const psVector* y,
     const psVector* z,
-    const psVector* dz
-)
+    const psVector* dz)
 {
     psVector *X;
@@ -2487 +2495 @@
 
 // XXX EAM : be careful here with F32 vs F64 vectors
+/*
+ Basically, you repetitively fit a polynomial to a set of data points,
+ reject the points that did not fit well, then refit the polynomial.
+ 
+ Basically, simply fit the polynomial to the data.  They compare the
+ fit (by evaluating the x data with that polynomial and subtracting
+ from the original f data).  That's the residual.  Loop through all
+ data and if the ((residual - mean) > 3*stDev), mask that data point,
+ and fit the polynomial again.  Do this 3 times.
+*/
+
 psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
                             psVector* mask,
@@ -2531 +2550 @@
 }
 
-
+/******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+NEW Code:
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+ *****************************************************************************/
+
+#define PS_VECTOR_GEN_INDEX_F32(VEC, SIZE) \
+VEC = psVectorAlloc(SIZE, PS_TYPE_F32); \
+for (psS32 i = 0 ; i < SIZE ; i++) { \
+    VEC->data.F32[i] = (psF32) i; \
+} \
+
+psPolynomial1D *psVectorFitPolynomial1D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *x32 = NULL;
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    } else {
+        PS_VECTOR_GEN_INDEX_F32(x32, f->n);
+    }
+
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        // XXX: Call EAM code
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        // XXX: Call my code
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+
+    // Free psVectors that were created for NULL arguments.
+    if (x == NULL) {
+        psFree(x32);
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+
+    return(NULL);
+}
+
+
+psPolynomial2D *psVectorFitPolynomial2D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial3D *psVectorFitPolynomial3D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial4D *psVectorFitPolynomial4D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+
+psPolynomial1D *psVectorClipFitPolynomial1D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *x32 = NULL;
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    // Free psVectors that were created for NULL arguments.
+    if (x == NULL) {
+        psFree(x32);
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial2D *psVectorClipFitPolynomial2D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial3D *psVectorClipFitPolynomial3D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial4D *psVectorClipFitPolynomial4D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+

trunk/psLib/src/math/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-08-30 01:14:13 $
+ *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-07 21:35:50 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -91 +91 @@
 );
 
+
+
+
+
+psPolynomial1D *psVectorFitPolynomial1D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x
+);
+
+psPolynomial2D *psVectorFitPolynomial2D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y
+);
+
+psPolynomial3D *psVectorFitPolynomial3D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z
+);
+
+psPolynomial4D *psVectorFitPolynomial4D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t
+);
+
+
+psPolynomial1D *psVectorClipFitPolynomial1D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x
+);
+
+psPolynomial2D *psVectorClipFitPolynomial2D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y
+);
+
+psPolynomial3D *psVectorClipFitPolynomial3D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z
+);
+
+psPolynomial4D *psVectorClipFitPolynomial4D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t
+);
+
+
+
+
 /** Derive a one-dimensional spline fit.
  *
@@ -230 +326 @@
 /* \} */// End of MathGroup Functions
 
+
+
+
 #endif // #ifndef PS_MINIMIZE_H
 

trunk/psLib/src/math/psStats.c

@@ -14 +14 @@
  *      stats->binsize
  *
- *  @version $Revision: 1.139 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-08-30 01:14:13 $
+ *  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-07 21:35:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -46 +46 @@
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-#define PS_GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
+#define PS_GAUSS_WIDTH 5                // The width of the Gaussian smoothing.
+// This corresponds to N in the ADD.
 #define PS_CLIPPED_NUM_ITER_LB 1
 #define PS_CLIPPED_NUM_ITER_UB 10
@@ -292 +293 @@
 max of the input vector.  If there was a problem with the max calculation,
 this routine sets stats->max to NAN.
+ 
+XXX: Do we need to factor errors into it?
  *****************************************************************************/
 psS32 p_psVectorMax(const psVector* myVector,
@@ -614 +617 @@
 /******************************************************************************
 p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
-robustHistogram with a Gaussian of width sigma.
+robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
+smoothed data.
  
 XXX: Only PS_TYPE_F32 is supported.
@@ -621 +625 @@
 call that.  Is that possible?
  *****************************************************************************/
-psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
+psVector *p_psVectorSmoothHistGaussian(psHistogram *histogram,
                                        psF32 sigma)
 {
-    PS_ASSERT_PTR_NON_NULL(robustHistogram, NULL);
-    PS_ASSERT_PTR_NON_NULL(robustHistogram->bounds, NULL);
-
-    psS32 i = 0;                  // Loop index variable
-    psS32 j = 0;                  // Loop index variable
-    psF32 iMid;
-    psF32 jMid;
-    psS32 numBins = robustHistogram->nums->n;
-    psS32 numBounds = robustHistogram->bounds->n;
-    psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32);
+    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
+
+    psS32 numBins = histogram->nums->n;
+    psS32 numBounds = histogram->bounds->n;
+    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32);
+    psF32 firstBound = histogram->bounds->data.F32[0];
+    psF32 lastBound = histogram->bounds->data.F32[numBounds-1];
+    psScalar x;
+    x.type.type = PS_TYPE_F32;
     psS32 jMin = 0;
     psS32 jMax = 0;
-    psF32 firstBound = robustHistogram->bounds->data.F32[0];
-    psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1];
+
+    if (histogram->uniform == false) {
+        //
+        // We get here if the histogram is non-uniform.
+        //
+
+        for (psS32 i = 0; i < numBins; i++) {
+            // Determine the midpoint of bin i.
+            psS32 iMid = PS_BIN_MIDPOINT(histogram, i);
+
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            //
+            x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
+            if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+                jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x);
+                if (jMin < 0) {
+                    psError(PS_ERR_UNEXPECTED_NULL,
+                            false,
+                            PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
+                    return(NULL);
+                }
+            } else if (x.data.F32 <= firstBound) {
+                jMin = 0;
+            } else if (x.data.F32 >= lastBound) {
+                jMin = histogram->bounds->n - 1;
+            }
+
+            x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
+            if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+                jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x);
+                if (jMax < 0) {
+                    psError(PS_ERR_UNEXPECTED_NULL,
+                            false,
+                            PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
+                    return(NULL);
+                }
+            } else if (x.data.F32 <= firstBound) {
+                jMax = 0;
+            } else if (x.data.F32 >= lastBound) {
+                jMax = histogram->bounds->n - 1;
+            }
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            for (psS32 j = jMin ; j <= jMax ; j++) {
+                psS32 jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
+        }
+    } else {
+        //
+        // We get here if the histogram is uniform.
+        //
+
+        for (psS32 i = 0; i < numBins; i++) {
+            psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0];
+            psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize);
+
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            //
+            psS32 jMin = i - gaussWidth;
+            if (jMin < 0 ) {
+                jMin = 0;
+            }
+            psS32 jMax = i + gaussWidth;
+            if (jMax > (histogram->bounds->n - 1)) {
+                jMax = (histogram->bounds->n - 1);
+            }
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            psS32 iMid = PS_BIN_MIDPOINT(histogram, i);
+            for (psS32 j = jMin ; j <= jMax ; j++) {
+                psS32 jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
+        }
+    }
+
+    return(smooth);
+}
+/******************************************************************************
+p_psVectorSmoothHistGaussianNEW(): This routine smoothes the data in the input
+robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
+smoothed data.
+ 
+XXX: Only PS_TYPE_F32 is supported.
+ 
+XXX: Write a general routine which smoothes a psVector.  This routine should
+call that.  Is that possible?
+ *****************************************************************************/
+psVector *p_psVectorSmoothHistGaussianNEW(psHistogram *histogram,
+        psF32 sigma)
+{
+    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
+
+    psS32 numBins = histogram->nums->n;
+    psS32 numBounds = histogram->bounds->n;
+    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32);
+    psF32 firstBound = histogram->bounds->data.F32[0];
+    psF32 lastBound = histogram->bounds->data.F32[numBounds-1];
     psScalar x;
-
     x.type.type = PS_TYPE_F32;
-    for (i = 0; i < numBins; i++) {
-        // Determine the midpoint of bin i.
-        iMid = (robustHistogram->bounds->data.F32[i] +
-                robustHistogram->bounds->data.F32[i+1]) / 2.0;
-
-
-        // We determine the bin numbers corresponding to a range of data
-        // values surrounding iMid.  The ranges is of size
-        // s*PS_GAUSS_WIDTH*sigma
-
-        // YYY: The p_psVectorBinDisect() routine does much of the work of
-        // the following conditionals, however, it also reports a warning
-        // message.  I don't want the warning message so I reproduce the
-        // conditionals here.  Maybe p_psVectorBinDisect() should not produce
-        // warnings?
-
-        x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
-        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
-            jMin = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x);
-            if (jMin < 0) {
-                psError(PS_ERR_UNEXPECTED_NULL,
-                        false,
-                        PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
-                return(NULL);
-            }
-        } else if (x.data.F32 <= firstBound) {
-            jMin = 0;
-        } else if (x.data.F32 >= lastBound) {
-            jMin = robustHistogram->bounds->n - 1;
-        }
-
-        x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
-        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
-            jMax = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x);
-            if (jMax < 0) {
-                psError(PS_ERR_UNEXPECTED_NULL,
-                        false,
-                        PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
-                return(NULL);
-            }
-        } else if (x.data.F32 <= firstBound) {
-            jMax = 0;
-        } else if (x.data.F32 >= lastBound) {
-            jMax = robustHistogram->bounds->n - 1;
-        }
-
-        smooth->data.F32[i] = 0.0;
-        for (j = jMin ; j <= jMax ; j++) {
-            jMid = (robustHistogram->bounds->data.F32[j] +
-                    robustHistogram->bounds->data.F32[j+1]) / 2.0;
-            smooth->data.F32[i] +=
-                robustHistogram->nums->data.F32[j] *
-                psGaussian(jMid, iMid, sigma, true);
+    psS32 jMin = 0;
+    psS32 jMax = 0;
+
+    if (histogram->uniform == false) {
+        //
+        // We get here if the histogram is non-uniform.
+        //
+
+        for (psS32 i = 0; i < numBins; i++) {
+            // Determine the midpoint of bin i.
+            psS32 iMid = PS_BIN_MIDPOINT(histogram, i);
+
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            //
+            x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
+            if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+                jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x);
+                if (jMin < 0) {
+                    psError(PS_ERR_UNEXPECTED_NULL,
+                            false,
+                            PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
+                    return(NULL);
+                }
+            } else if (x.data.F32 <= firstBound) {
+                jMin = 0;
+            } else if (x.data.F32 >= lastBound) {
+                jMin = histogram->bounds->n - 1;
+            }
+
+            x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
+            if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+                jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x);
+                if (jMax < 0) {
+                    psError(PS_ERR_UNEXPECTED_NULL,
+                            false,
+                            PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
+                    return(NULL);
+                }
+            } else if (x.data.F32 <= firstBound) {
+                jMax = 0;
+            } else if (x.data.F32 >= lastBound) {
+                jMax = histogram->bounds->n - 1;
+            }
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            for (psS32 j = jMin ; j <= jMax ; j++) {
+                psS32 jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
+        }
+    } else {
+        //
+        // We get here if the histogram is uniform.
+        //
+
+        for (psS32 i = 0; i < numBins; i++) {
+            psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0];
+            psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize);
+
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            //
+            psS32 jMin = i - gaussWidth;
+            if (jMin < 0 ) {
+                jMin = 0;
+            }
+            psS32 jMax = i + gaussWidth;
+            if (jMax > (histogram->bounds->n - 1)) {
+                jMax = (histogram->bounds->n - 1);
+            }
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            psS32 iMid = PS_BIN_MIDPOINT(histogram, i);
+            for (psS32 j = jMin ; j <= jMax ; j++) {
+                psS32 jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
         }
     }
@@ -1142 +1296 @@
 These macros and functions define the following functions:
  
-<    p_psNormalizeVectorRange(myData, low, high)
+    p_psNormalizeVectorRange(myData, low, high)
  
 That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
@@ -1251 +1405 @@
 
 /******************************************************************************
-p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
-as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
-polynomials here) and a range of x-values for which it is defined:
-[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
-that f(x) == midpoint.  This functions uses a binary-search algorithm on the
-range and assumes that the polynomial is monotonically increasing or
-decreasing within that range.
+p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, getThisValue): This routine
+takes as input a 1-D polynomial of arbitrary order and a range of x-values for
+which it is defined:  [rangeLow, rangeHigh].  It determines the x-value of
+that polynomial such that f(x) == getThisValue.  This function uses a
+binary-search algorithm on the range and assumes that the polynomial is
+monotonically increasing or decreasing within that range.
  
 XXX: Terminate when f(x)-getThisValue is within some error tolerance.
@@ -1269 +1422 @@
 {
     PS_ASSERT_POLY_NON_NULL(myPoly, NAN);
-    PS_FLOAT_COMPARE(rangeLow, rangeHigh, NAN);
+    PS_ASSERT_FLOAT_LARGER_THAN(rangeHigh, rangeLow, NAN);
     // We ensure that the requested f(y) value, which is getThisValue, is
     // falls within the range of y-values of the polynomial "myPoly" in the
     // specified x-range (rangeLow:rangeHigh).
-    psF32 fLo = psPolynomial1DEval(
-                    myPoly,
-                    rangeLow
-                );
-    psF32 fHi = psPolynomial1DEval(
-                    myPoly,
-                    rangeHigh
-                );
+    psF32 fLo = psPolynomial1DEval(myPoly, rangeLow);
+    psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh);
     if (!((fLo <= getThisValue) && (fHi >= getThisValue))) {
         psError(PS_ERR_UNKNOWN,
@@ -1300 +1447 @@
         oldMidpoint = midpoint;
 
-        f = psPolynomial1DEval(
-                myPoly,
-                midpoint
-            );
+        f = psPolynomial1DEval(myPoly, midpoint);
         if (fabs(f - getThisValue) <= FLT_EPSILON) {
             return (midpoint);
@@ -1331 +1475 @@
 XXX: the vectors do not have to be the same length.  Must insert the proper
 tests to ensure that binNum is within acceptable ranges for both vectors.
+ 
+XXX: This currently assumes that the three points are monotonically increasing
+or decreasing: so, it works for the cumulative histogram vectors, but not for
+arbitrary vectors.  We should probably test that condition.
 *****************************************************************************/
 psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
@@ -1345 +1493 @@
     PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (yVec->n - 1), NAN);
 
-    //    PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);
-    //    PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);
-    //    PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);
-    //    PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD);
     psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     psVector *y = psVectorAlloc(3, PS_TYPE_F64);
@@ -1364 +1508 @@
         y->data.F64[1] = yVec->data.F32[binNum];
         y->data.F64[2] = yVec->data.F32[binNum + 1];
+
+        //
+        // Ensure that the y values are monotonic.
+        //
+        // XXX: This routine should probably be rewritten in a more general fashion
+        // so that the folloiwng checks are not necessary.
+        //
+        if (y->data.F64[0] < y->data.F64[1]) {
+            if (!(y->data.F64[1] <= y->data.F64[2])) {
+                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
+                psFree(myPoly);
+                psFree(x);
+                psFree(y);
+                psFree(yErr);
+                return(NAN);
+            }
+        } else {
+            if (!(y->data.F64[1] >= y->data.F64[2])) {
+                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
+                psFree(myPoly);
+                psFree(x);
+                psFree(y);
+                psFree(yErr);
+                return(NAN);
+            }
+        }
 
         // Ensure that yVal is within the range of the bins we are using.
@@ -1637 +1807 @@
         // code no longer produces sensible results.
         // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably
-        // remove some of theabove code that calculated the initial estimate
+        // remove some of the above code that calculated the initial estimate
         // for the mean and sigma.
 
@@ -1645 +1815 @@
             int index = i - modeBinNum + dL;
             // XXX: Should this be the natural log?
-            y->data.F32[index] = robustHistogramVector->data.F32[i];
-            //            y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);
+            //            y->data.F32[index] = robustHistogramVector->data.F32[i];
+            y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);
             x->data.F32[index] = (psF32) index;
         }
@@ -1742 +1912 @@
     psFree(robustHistogramVector);
     psFree(cumulativeRobustSums);
+
     return(0);
 }
+
+
+/*****************************************************************************
+XXX: Is there a psLib function for this?
+ *****************************************************************************/
+psVector *PsVectorDup(psVector *in)
+{
+    psVector *out = psVectorAlloc(in->n, in->type.type);
+
+    if (in->type.type == PS_TYPE_F32) {
+        for (psS32 i = 0 ; i < in->n ; i++) {
+            out->data.F32[i] = in->data.F32[i];
+        }
+    } else if (in->type.type == PS_TYPE_F64) {
+        for (psS32 i = 0 ; i < in->n ; i++) {
+            out->data.F64[i] = in->data.F64[i];
+        }
+    } else {
+        printf("XXX: Generate an error here.\n");
+        return(NULL);
+    }
+    return(out);
+}
+
+/******************************************************************************
+XXX: This function need to be written.  Actually, it simply needs to be
+retrieved from the CVS repository, since it was written earlier, then
+discarded.  At present, it was deleted from the CVS repository, so we might
+have to retrieve it from tape.
+*****************************************************************************/
+psVector *Fit1DGaussian(psVector *x, psVector*y)
+{
+    printf("XXX: Generate an error here.\n");
+    printf("XXX: Error: This function was previously part of psStats.c, was removed, was purged from CVS, and now needs to be retrieved from tape.\n");
+    return(NULL);
+}
+
+/******************************************************************************
+ 
+p_psVectorRobustStatsNew(myVector, maskVector, maskVal, stats): This is the new
+version of the robust stats routine.
+ 
+XXX: MUST DO: If the errors in the input values are known, then the same
+approach is used, except that the histograms become probability density
+functions (PDFs). In this case, the input values are spread out, so that they
+do not simply contribute a single unit to the histogram, but rather contribute
+a fraction of a value, equivalent to the weight. In the interests of speed, a
+boxcar PDF may be used to represent each input value (as opposed to a
+Gaussian), where the boxcar width is equal to 2p2 ln 2 times the error and
+each input value contributes constant area. Then the robust median and
+standard deviation are estimated in the same manner as above.
+ 
+XXX: Check for errors in psLib routines that we call.
+*****************************************************************************/
+psS32 p_psVectorRobustStatsNew(const psVector* myVector,
+                               const psVector* errors,
+                               const psVector* maskVector,
+                               psU32 maskVal,
+                               psStats* stats)
+{
+    psHistogram *robustHistogram = NULL;
+    psHistogram *cumulativeRobustHistogram = NULL;
+    psS32 numBins = 0;
+    psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
+    tmpScalar->type.type = PS_TYPE_F32;
+    psS32 totalDataPoints = 0;
+    psS32 rc = 0;
+    psVector *tmpMaskVec = PsVectorDup((psVector *) maskVector);
+
+    while (1) {
+        //
+        // Determine the bin size of the robust histogram.  This is done
+        // by computing the total range of data values and dividing by 1000.0.
+        //
+        psStats* tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX);
+        rc = p_psVectorMin(myVector, tmpMaskVec, maskVal, tmpStatsMinMax);
+        rc|= p_psVectorMax(myVector, tmpMaskVec, maskVal, tmpStatsMinMax);
+        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(tmpStatsMinMax);
+            psFree(tmpMaskVec);
+            psFree(tmpScalar);
+            return(1);
+        }
+        psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f;
+
+        //
+        // If all data points have the same value, then we set the appropiate
+        // members of stats and return.
+        //
+        if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) {
+            if (stats->options & PS_STAT_ROBUST_MEDIAN) {
+                stats->robustMedian = tmpStatsMinMax->min;
+            }
+            if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+                stats->robustUQ = tmpStatsMinMax->min;
+                stats->robustLQ = tmpStatsMinMax->min;
+            }
+            // XXX: Set these to the number of unmasked data points?
+            stats->robustNfit = 0.0;
+            stats->robustN50 = 0.0;
+            psFree(tmpStatsMinMax);
+            psFree(tmpMaskVec);
+            psFree(tmpScalar);
+
+            return(0);
+        }
+
+        //
+        // ADD: Step 0.
+        // Construct the histogram with the specified bin size.
+        //
+        // NOTE: we can not specify the bin size precisely since the argument
+        // to psHistogramAlloc() is the number of bins, not the binSize.
+        // If we get here, we know that binSize != 0.0.
+        //
+        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize);
+        robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
+        cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
+
+        // Populate the histogram array.
+        psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal);
+
+        //
+        // ADD: Step 1.
+        // Construct the cumulative histogram from the specific histogram
+        //
+        cumulativeRobustHistogram->nums->data.F32[0] = robustHistogram->nums->data.F32[0];
+        for (psS32 i = 1 ; i < robustHistogram->nums->n ; i++) {
+            cumulativeRobustHistogram->nums->data.F32[i] = cumulativeRobustHistogram->nums->data.F32[i-1] +
+                    robustHistogram->nums->data.F32[i];
+        }
+
+        //
+        // ADD: Step 2.
+        // Find the bin which contains the 50% data point.
+        //
+        totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1];
+        tmpScalar->data.F32 = totalDataPoints/2.0;
+        psS32 binMedian = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+        if (binMedian != 0) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50% data point.\n");
+            psFree(tmpStatsMinMax);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpScalar);
+            return(1);
+        }
+
+        //
+        // ADD: Step 3.
+        // Interpolate to the exact 50% position: this is the robust histogram median.
+        // XXX: Check for errors here!
+        //
+        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
+                                  *(psVector* *)&robustHistogram->bounds,
+                                  *(psVector* *)&robustHistogram->nums,
+                                  binMedian,
+                                  totalDataPoints/2.0);
+
+        //
+        // ADD: Step 4.
+        // Find the bins which contains the 15.8655% and 84.1345% data points.
+        //
+        tmpScalar->data.F32 = totalDataPoints * 0.158655f;
+        psS32 binLo = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+        tmpScalar->data.F32 = totalDataPoints * 0.841345f;
+        psS32 binHi = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+        if ((binLo != 0) || (binHi != 0)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the15.8655% and 84.1345% data point\n");
+            psFree(tmpStatsMinMax);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpScalar);
+            return(1);
+        }
+
+        //
+        // ADD: Step 4b.
+        // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
+        //
+        psF32 binLoF32 = fitQuadraticSearchForYThenReturnX(
+                             *(psVector* *)&robustHistogram->bounds,
+                             *(psVector* *)&robustHistogram->nums,
+                             binLo,
+                             totalDataPoints * 0.158655f);
+        psF32 binHiF32 = fitQuadraticSearchForYThenReturnX(
+                             *(psVector* *)&robustHistogram->bounds,
+                             *(psVector* *)&robustHistogram->nums,
+                             binHi,
+                             totalDataPoints * 0.841345f);
+
+        //
+        // ADD: Step 5.
+        // Determine SIGMA as 1/2 of the distance between these positions.
+        //
+        psF32 sigma = (binHiF32 - binLoF32) / 2.0;
+
+        //
+        // ADD: Step 6.
+        // If the measured SIGMA is less than 2 times the bin size, exclude
+        // points which are more than 25 bins from the median,
+        // recalculate the bin size, and perform the algorithm again.
+        //
+        if (sigma < (2 * binSize)) {
+            psF32 medianLo = robustHistogram->bounds->data.F32[binMedian - 25];
+            psF32 medianHi = robustHistogram->bounds->data.F32[binMedian + 25];
+            for (psS32 i = 0 ; i < myVector->n ; i++) {
+                if ((myVector->data.F32[i] < medianLo) ||
+                        (myVector->data.F32[i] > medianHi)) {
+                    tmpMaskVec->data.U8[i] = 1;
+                }
+            }
+        } else {
+            //
+            // ADD: Step 7.
+            // Find the bins which contains the 25% and 75% data points.
+            //
+            tmpScalar->data.F32 = totalDataPoints * 0.25f;
+            psS32 binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            tmpScalar->data.F32 = totalDataPoints * 0.75f;
+            psS32 binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            if ((binLo25 != 0) || (binHi25 != 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25% and 75% data points\n");
+                psFree(tmpStatsMinMax);
+                psFree(robustHistogram);
+                psFree(cumulativeRobustHistogram);
+                psFree(tmpScalar);
+                return(1);
+            }
+
+            //
+            // ADD: Step 8.
+            // Interpolate to find these two positions exactly: these are the upper
+            // and lower quartile positions.
+            // XXX: Check for errors.
+            //
+            psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
+                                   *(psVector* *)&robustHistogram->bounds,
+                                   *(psVector* *)&robustHistogram->nums,
+                                   binLo25,
+                                   totalDataPoints * 0.25f);
+            psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
+                                   *(psVector* *)&robustHistogram->bounds,
+                                   *(psVector* *)&robustHistogram->nums,
+                                   binHi25,
+                                   totalDataPoints * 0.75f);
+
+            stats->robustLQ = binLo25F32;
+            stats->robustUQ = binHi25F32;
+            // XXX: No idea how to calculate stats->stdev
+
+            // PS_BIN_MIDPOINT(robustHistogram, modeBinNum);
+
+            // XXX: I think sumNfit == sumN50 here.
+            stats->robustNfit = -1;
+            stats->robustN50 = -1;
+
+            //
+            // Perform the Robust Histogram Statistics algorithm above
+            //
+
+            //
+            // Smooth the resulting histogram with a Gaussian with SIGMA_s = 1
+            // bin.
+            //
+            // XXX: SIGMA_s is defined nowhere in the ADD.
+            //
+            psF32 SIGMA_S = 1.0;
+            p_psVectorSmoothHistGaussian(robustHistogram, SIGMA_S);
+
+            //
+            // Find the bin with the peak value in the range 2 SIGMA of the
+            // robust histogram median.
+            //
+            // XXX: SIGMA is defined nowhere in the ADD.
+            //
+            psF32 SIGMA = 2.0;
+            psS32 binMin = binMedian - (SIGMA * PS_GAUSS_WIDTH);
+            if (binMin < 0) {
+                binMin = 0;
+            }
+            psS32 binMax = binMedian + (2 * PS_GAUSS_WIDTH);
+            if (binMin > (robustHistogram->nums->n - 1)) {
+                binMin = (robustHistogram->nums->n - 1);
+            }
+            psS32 binNum = binNum;
+            psF32 binMaxNums = robustHistogram->nums->data.F32[binNum];
+            for (psS32 i = binNum+1 ; i <= binMax ; i++) {
+                if (robustHistogram->nums->data.F32[i] > binMaxNums) {
+                    binNum = i;
+                    binMaxNums = robustHistogram->nums->data.F32[i];
+                }
+            }
+
+            //
+            // Fit a Gaussian to the bins in the range 2 SIGMA of the robust
+            // histogram median.
+            //
+            // XXX: SIGMA is defined nowhere in the ADD.
+            //
+            psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
+            psVector *x = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
+            psS32 j = 0;
+            for (psS32 i = binNum ; i <= binMax ; i++) {
+                y->data.F32[j] = robustHistogram->nums->data.F32[i];
+                x->data.F32[j] = PS_BIN_MIDPOINT(robustHistogram, i);
+            }
+            //
+            // XXX: This function need to be written.  Actually, it simply
+            // needs to be retrieved from the CVS repository, since it was
+            // written earlier, then discarded.  At present, it was deleted
+            // from the CVS repository, so we might have to retrieve it from
+            // tape.
+            //
+            psVector *params = Fit1DGaussian(x, y);
+
+            //
+            // The robust mean mean_r is derived directly from the fitted
+            // Gaussian mean.
+            //
+            stats->robustMean = params->data.F32[0];
+
+            //
+            // The robust standard deviation, SIGMA_r is determined by
+            // subtracting the smoothing scale in quadrature: SIGMA_r^2 = SIGMA^2
+            // - SIGMA_s^2
+            //
+            // XXX: SIGMA and SIGMA_s are defined nowhere in the ADD.  We must figure
+            // out what they are.
+            //
+            stats->robustStdev = sqrt(PS_SQR(SIGMA) - PS_SQR(SIGMA_S));
+
+            psFree(tmpStatsMinMax);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpScalar);
+            psFree(params);
+
+            return(0);
+        }
+
+        psFree(tmpStatsMinMax);
+        psFree(robustHistogram);
+        psFree(cumulativeRobustHistogram);
+    }
+    return(1);
+}
+
+
+
+
+
+
+
+
 
 /*****************************************************************************/