Index: trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.c	(revision 42823)
+++ trunk/psLib/src/math/psMinimizePolyFit.c	(revision 42824)
@@ -69,4 +69,11 @@
 if ((ORIG != NULL) && (ORIG->type.type != PS_TYPE_F64)) { psFree(TEMP); }
 
+psVector *psVector_GetModifiedErrors_Caucy
+( const psVector *f,
+  const psVector *fEval,
+  const psVector *fErr,
+  const psVector *mask,
+  psVectorMaskType maskValue);
+
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -90,8 +97,8 @@
 returned as a psVector sums.
 *****************************************************************************/
-static psVector *BuildSums1D(
-    psVector* sums,
-    psF64 x,
-    psS32 nTerm)
+static psVector *BuildSums1D
+( psVector* sums,
+  psF64 x,
+  psS32 nTerm)
 {
     psS32 nSum = 0;
@@ -710,5 +717,7 @@
     switch (poly->type) {
     case PS_POLYNOMIAL_ORD:
-      if ((f64->n < 10000)&&(poly->nX > 1)) {
+      // if x is NULL, the domain is the index of the vector, in which case we do not
+      // need to rescale (range is guaranteed to be < 10000)
+      if ((x != NULL) && (f64->n < 10000) && (poly->nX > 1)) {
 	scale = true;
 
@@ -752,6 +761,5 @@
 	result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
 	psFree(z64); // Done with this.
-      }
-      else {
+      } else {
 	result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
       }
@@ -1030,4 +1038,102 @@
 }
 
+// These should probably be tunable:
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+# define WEIGHT_THRESHOLD 0.3
+
+// This function accepts F32 and F64 input vectors.
+bool psVectorIRLSFitPolynomial1D(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *xIn)
+{
+    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    if (mask != NULL) {
+	PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
+	PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
+    }
+    if (xIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
+    }
+
+    // Internal pointers for possibly NULL vectors.
+    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+
+    // initial fit with nominal errors
+    if (!psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x)) {
+	psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+	return false;
+    }
+
+    // use polyOld to save the last fit
+    psPolynomial1D *polyOld = NULL;
+
+    // use clipIter as max number of iterations
+    bool converged = false;
+    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
+        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n", N);
+
+	// evaluate the fit at the input positions
+	psVector *fEval = psPolynomial1DEvalVector (poly, x);
+
+	// calculate modified errors based on the deviation from the fit
+	psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
+	psFree (fEval);
+
+	// save the last fit (recycle the structure once allocated)
+	polyOld = psPolynomial1DCopy (polyOld, poly);
+
+	// calculate a new fit with modified errors:
+        if (!psVectorFitPolynomial1D(poly, mask, maskValue, f, modErr, x)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+            psFree(x);
+	    psFree(modErr);
+            return false;
+        }
+
+	// has the solution converged?
+	converged = true;
+	for (int ix = 0; ix <= poly->nX; ix++) {
+	  if ((fabs(poly->coeff[ix] - polyOld->coeff[ix]) > FIT_TOLERANCE * fabs(poly->coeff[ix])) && 
+	      (fabs(poly->coeff[ix] - polyOld->coeff[ix]) > FLT_TOLERANCE))
+	    converged = false;
+	}
+
+# if (0)	
+	// XXX test:
+	FILE *ftest = fopen ("irls.wt.dat", "w");
+	for (int i = 0; i < modErr->n; i++) {
+	    if (modErr->type.type == PS_TYPE_F64) {
+		fprintf (ftest, "%d %f\n", i, modErr->data.F64[i]);
+	    } else {
+		fprintf (ftest, "%d %f\n", i, modErr->data.F32[i]);
+	    }
+	}
+	fclose (ftest);
+# endif
+	psFree (modErr);
+    }
+
+    // Free local temporary variables
+    psFree(x);
+    psFree(polyOld);
+
+    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
+    return true;
+}
 
 /******************************************************************************
@@ -1613,4 +1719,97 @@
 }
 
+// This function accepts F32 and F64 input vectors.
+bool psVectorIRLSFitPolynomial2D(
+    psPolynomial2D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *xIn,
+    const psVector *yIn)
+{
+    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
+
+    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    if (mask != NULL) {
+	PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
+	PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
+    }
+    if (xIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
+    }
+    if (yIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(yIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(yIn, f->type.type, false);
+    }
+
+    // Internal pointers for possibly NULL vectors.  
+    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+    psVector *y = (yIn != NULL) ? psMemIncrRefCounter((psVector *) yIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+
+    // initial fit with nominal errors
+    if (!psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y)) {
+	psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+	psFree(x);
+	psFree(y);
+	return false;
+    }
+
+    // use polyOld to save the last fit
+    psPolynomial2D *polyOld = NULL;
+
+    // use clipIter as max number of iterations
+    bool converged = false;
+    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
+        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial2D()\n", N);
+
+	// evaluate the fit at the input positions
+	psVector *fEval = psPolynomial2DEvalVector (poly, x, y);
+
+	// calculate modified errors based on the deviation from the fit
+	psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
+	psFree (fEval);
+
+	// save the last fit (recycle the structure once allocated)
+	polyOld = psPolynomial2DCopy (polyOld, poly);
+
+	// calculate a new fit with modified errors:
+        if (!psVectorFitPolynomial2D(poly, mask, maskValue, f, modErr, x, y)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+            psFree(x);
+            psFree(y);
+	    psFree(modErr);
+            return false;
+        }
+
+	// has the solution converged?
+	converged = true;
+	for (int ix = 0; ix <= poly->nX; ix++) {
+	    for (int iy = 0; iy <= poly->nY; iy++) {
+		if ((fabs(poly->coeff[ix][iy] - polyOld->coeff[ix][iy]) > FIT_TOLERANCE * fabs(poly->coeff[ix][iy])) && 
+		    (fabs(poly->coeff[ix][iy] - polyOld->coeff[ix][iy]) > FLT_TOLERANCE))
+		    converged = false;
+	    }
+	}
+	psFree (modErr);
+    }
+
+    // Free local temporary variables
+    psFree(x);
+    psFree(y);
+    psFree(polyOld);
+
+    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
+    return true;
+}
 
 /******************************************************************************
@@ -2016,4 +2215,109 @@
     // Free local temporary variables
     psFree(resid);
+
+    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
+    return true;
+}
+
+// This function accepts F32 and F64 input vectors.
+bool psVectorIRLSFitPolynomial3D(
+    psPolynomial3D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *xIn,
+    const psVector *yIn,
+    const psVector *zIn)
+{
+    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
+
+    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    if (mask != NULL) {
+	PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
+	PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
+    }
+    if (xIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
+    }
+    if (yIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(yIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(yIn, f->type.type, false);
+    }
+    if (zIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(zIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(zIn, f->type.type, false);
+    }
+
+    // Internal pointers for possibly NULL vectors.  
+    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+    psVector *y = (yIn != NULL) ? psMemIncrRefCounter((psVector *) yIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+    psVector *z = (zIn != NULL) ? psMemIncrRefCounter((psVector *) zIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+
+    // initial fit with nominal errors
+    if (!psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z)) {
+	psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+	psFree(x);
+	psFree(y);
+	psFree(z);
+	return false;
+    }
+
+    // use polyOld to save the last fit
+    psPolynomial3D *polyOld = NULL;
+
+    // use clipIter as max number of iterations
+    bool converged = false;
+    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
+        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial3D()\n", N);
+
+	// evaluate the fit at the input positions
+	psVector *fEval = psPolynomial3DEvalVector (poly, x, y, z);
+
+	// calculate modified errors based on the deviation from the fit
+	psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
+	psFree (fEval);
+
+	// save the last fit (recycle the structure once allocated)
+	polyOld = psPolynomial3DCopy (polyOld, poly);
+
+	// calculate a new fit with modified errors:
+        if (!psVectorFitPolynomial3D(poly, mask, maskValue, f, modErr, x, y, z)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+            psFree(x);
+            psFree(y);
+            psFree(z);
+	    psFree(modErr);
+            return false;
+        }
+
+	// has the solution converged?
+	converged = true;
+	for (int ix = 0; ix <= poly->nX; ix++) {
+	    for (int iy = 0; iy <= poly->nY; iy++) {
+		for (int iz = 0; iz <= poly->nZ; iz++) {
+		    if ((fabs(poly->coeff[ix][iy][iz] - polyOld->coeff[ix][iy][iz]) > FIT_TOLERANCE * fabs(poly->coeff[ix][iy][iz])) && 
+			(fabs(poly->coeff[ix][iy][iz] - polyOld->coeff[ix][iy][iz]) > FLT_TOLERANCE))
+			converged = false;
+		}
+	    }
+	}
+	psFree (modErr);
+    }
+
+    // Free local temporary variables
+    psFree(x);
+    psFree(y);
+    psFree(z);
+    psFree(polyOld);
 
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
@@ -2455,2 +2759,169 @@
     return true;
 }
+
+// This function accepts F32 and F64 input vectors.
+bool psVectorIRLSFitPolynomial4D(
+    psPolynomial4D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *xIn,
+    const psVector *yIn,
+    const psVector *zIn,
+    const psVector *tIn)
+{
+    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
+
+    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
+    PS_ASSERT_POLY_NON_NULL(poly, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
+    if (mask != NULL) {
+	PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
+	PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
+        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
+    }
+    if (xIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
+    }
+    if (yIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(yIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(yIn, f->type.type, false);
+    }
+    if (zIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(zIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(zIn, f->type.type, false);
+    }
+    if (tIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(tIn, f, false);
+        PS_ASSERT_VECTOR_TYPE(tIn, f->type.type, false);
+    }
+
+    // Internal pointers for possibly NULL vectors.  
+    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+    psVector *y = (yIn != NULL) ? psMemIncrRefCounter((psVector *) yIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+    psVector *z = (zIn != NULL) ? psMemIncrRefCounter((psVector *) zIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+    psVector *t = (tIn != NULL) ? psMemIncrRefCounter((psVector *) tIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
+
+    // initial fit with nominal errors
+    if (!psVectorFitPolynomial4D(poly, mask, maskValue, f, fErr, x, y, z, t)) {
+	psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+	psFree(x);
+	psFree(y);
+	psFree(z);
+	psFree(t);
+	return false;
+    }
+
+    // use polyOld to save the last fit
+    psPolynomial4D *polyOld = NULL;
+
+    // use clipIter as max number of iterations
+    bool converged = false;
+    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
+        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial4D()\n", N);
+
+	// evaluate the fit at the input positions
+	psVector *fEval = psPolynomial4DEvalVector (poly, x, y, z, t);
+
+	// calculate modified errors based on the deviation from the fit
+	psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
+	psFree (fEval);
+
+	// save the last fit (recycle the structure once allocated)
+	polyOld = psPolynomial4DCopy (polyOld, poly);
+
+	// calculate a new fit with modified errors:
+        if (!psVectorFitPolynomial4D(poly, mask, maskValue, f, modErr, x, y, z, t)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
+            psFree(x);
+            psFree(y);
+            psFree(z);
+            psFree(t);
+	    psFree(modErr);
+            return false;
+        }
+
+	// has the solution converged?
+	converged = true;
+	for (int ix = 0; ix <= poly->nX; ix++) {
+	    for (int iy = 0; iy <= poly->nY; iy++) {
+		for (int iz = 0; iz <= poly->nZ; iz++) {
+		    for (int it = 0; it <= poly->nT; it++) {
+			if ((fabs(poly->coeff[ix][iy][iz][it] - polyOld->coeff[ix][iy][iz][it]) > FIT_TOLERANCE * fabs(poly->coeff[ix][iy][iz][it])) && 
+			    (fabs(poly->coeff[ix][iy][iz][it] - polyOld->coeff[ix][iy][iz][it]) > FLT_TOLERANCE))
+			    converged = false;
+		    }
+		}
+	    }
+	}
+	psFree (modErr);
+    }
+
+    // Free local temporary variables
+    psFree(x);
+    psFree(y);
+    psFree(z);
+    psFree(t);
+    psFree(polyOld);
+
+    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
+    return true;
+}
+
+// ######################## utilities ###################
+
+// Used by IRLS fitting
+// This function assumes the input vectors (f, fEval, fErr) all have the same type
+// This requirement is already enforced in the calling function (psVectorIRLSFitPolynomial1D)
+psVector *psVector_GetModifiedErrors_Caucy (
+    const psVector *f,
+    const psVector *fEval,
+    const psVector *fErr,
+    const psVector *mask,
+    psVectorMaskType maskValue) {
+
+  psVector *mErr = psVectorAlloc (f->n, f->type.type);
+
+  switch (f->type.type) {
+    case PS_TYPE_F64:
+      {
+	  psF64 *fPtr     = f->data.F64;
+	  psF64 *fErrPtr  = fErr->data.F64;
+	  psF64 *mErrPtr  = mErr->data.F64;
+	  psF64 *fEvalPtr = fEval->data.F64;
+	  for (int i = 0; i < f->n; i++, fPtr++, fEvalPtr++, fErrPtr++, mErrPtr++) {
+	      if (mask && (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValue)) continue;
+	      double dF = (*fPtr - *fEvalPtr) / 2.385;
+	      double dV = PS_SQR(*fErrPtr) + PS_SQR(dF);
+	      *mErrPtr = sqrt(dV);
+	  }
+      }
+      break;
+    case PS_TYPE_F32:
+      {
+	  psF32 *fPtr     = f->data.F32;
+	  psF32 *fErrPtr  = fErr->data.F32;
+	  psF32 *mErrPtr  = mErr->data.F32;
+	  psF32 *fEvalPtr = fEval->data.F32;
+	  for (int i = 0; i < f->n; i++, fPtr++, fEvalPtr++, fErrPtr++, mErrPtr++) {
+	      if (mask && (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValue)) continue;
+	      double dF = (*fPtr - *fEvalPtr) / 2.385;
+	      double dV = PS_SQR(*fErrPtr) + PS_SQR(dF);
+	      *mErrPtr = sqrt(dV);
+	  }
+      }
+      break;
+    default:
+        psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n");
+	return NULL;
+  }
+  return mErr;
+}
+
Index: trunk/psLib/src/math/psMinimizePolyFit.h
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.h	(revision 42823)
+++ trunk/psLib/src/math/psMinimizePolyFit.h	(revision 42824)
@@ -134,4 +134,50 @@
 );
 
+bool psVectorIRLSFitPolynomial1D(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x
+);
+
+bool psVectorIRLSFitPolynomial2D(
+    psPolynomial2D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y
+);
+
+bool psVectorIRLSFitPolynomial3D(
+    psPolynomial3D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z
+);
+
+bool psVectorIRLSFitPolynomial4D(
+    psPolynomial4D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psVectorMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t
+);
+
 /// @}
 #endif // #ifndef PS_MINIMIZE_POLYFIT_H