Index: trunk/psLib/src/math/Makefile.am
===================================================================
--- trunk/psLib/src/math/Makefile.am	(revision 7375)
+++ trunk/psLib/src/math/Makefile.am	(revision 7380)
@@ -8,15 +8,19 @@
 	psBinaryOp.c \
 	psCompare.c \
+	psEllipse.c \
 	psMatrix.c \
 	psMinimizeLMM.c \
 	psMinimizePowell.c \
 	psMinimizePolyFit.c \
+	psPolynomial.c \
+	psPolynomialUtils.c \
 	psRandom.c \
 	psRegion.c \
 	psRegionForImage.c \
-	psPolynomial.c \
+	psSparse.c \
 	psSpline.c \
 	psStats.c \
-	psMathUtils.c
+	psMathUtils.c \
+	psVectorSmooth.c
 
 EXTRA_DIST = math.i
@@ -28,15 +32,19 @@
 	psCompare.h \
 	psConstants.h \
+	psEllipse.h \
 	psMatrix.h \
 	psMinimizeLMM.h \
 	psMinimizePowell.h \
 	psMinimizePolyFit.h \
+	psPolynomial.h \
+	psPolynomialUtils.h \
 	psRandom.h \
 	psRegion.h \
 	psRegionForImage.h \
-	psPolynomial.h \
+	psSparse.h \
 	psSpline.h \
 	psStats.h \
-	psMathUtils.h
+	psMathUtils.h \
+	psVectorSmooth.h
 
 CLEANFILES = *~
Index: trunk/psLib/src/math/psEllipse.c
===================================================================
--- trunk/psLib/src/math/psEllipse.c	(revision 7380)
+++ trunk/psLib/src/math/psEllipse.c	(revision 7380)
@@ -0,0 +1,61 @@
+#include <stdio.h>
+
+#include "psConstants.h"
+#include "psEllipse.h"
+
+psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments)
+{
+    psEllipseAxes axes;
+
+    double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.xy));
+    if (f > (moments.x2 + moments.y2) / 2.0) {
+        f = 0.98*(moments.x2 + moments.y2) / 2.0;
+    }
+
+    axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f);
+    axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f);
+    axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2;
+    // theta in radians
+
+    return axes;
+}
+
+psEllipseShape psEllipseAxesToShape(psEllipseAxes axes)
+{
+    psEllipseShape shape;
+
+    double r1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor);
+    double r2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor);
+
+    double sxr = r1 + r2*cos(2*axes.theta);
+    double syr = r1 - r2*cos(2*axes.theta);
+
+    shape.sx = 1.0 / sqrt(sxr);
+    shape.sy = 1.0 / sqrt(syr);
+    shape.sxy = r2*sin(2*axes.theta);
+
+    return shape;
+}
+
+psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape)
+{
+    psEllipseAxes axes;
+
+    double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy);
+    double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy);
+
+    // force the axis ratio to be less than 10
+    double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2));
+
+    shape.sxy = PS_MIN(PS_MAX(shape.sxy, -r1), r1);
+
+    axes.theta = atan2 (-2.0*shape.sxy, f2) / 2.0;
+
+    double Ar = 0.25*f1 + 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
+    double Br = 0.25*f1 - 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
+
+    axes.minor = 1.0 / sqrt (Ar);
+    axes.major = 1.0 / sqrt (Br);
+
+    return axes;
+}
Index: trunk/psLib/src/math/psEllipse.h
===================================================================
--- trunk/psLib/src/math/psEllipse.h	(revision 7380)
+++ trunk/psLib/src/math/psEllipse.h	(revision 7380)
@@ -0,0 +1,40 @@
+/** @file  psEllipse.h
+ *
+ * functions to manipulate sparse matrices equations
+ *
+ */
+
+#ifndef PS_ELLIPSE_H
+#define PS_ELLIPSE_H
+
+// strucures to define elliptical shape parameters
+typedef struct
+{
+    double major;                       // Major axis
+    double minor;                       // Minor axis
+    double theta;                       // Position angle
+}
+psEllipseAxes;
+
+typedef struct
+{
+    double x2;                          // Moment of xx
+    double y2;                          // Moment of yy
+    double xy;                          // Moment of xy
+}
+psEllipseMoments;
+
+typedef struct
+{
+    double sx;                          // Shape parameter in x
+    double sy;                          // Shape parameter in y
+    double sxy;                         // Shape parameter in xy
+}
+psEllipseShape;
+
+// conversions between elliptical shape representations
+psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments);
+psEllipseShape psEllipseAxesToShape(psEllipseAxes axes);
+psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape);
+
+#endif
Index: trunk/psLib/src/math/psPolynomialUtils.c
===================================================================
--- trunk/psLib/src/math/psPolynomialUtils.c	(revision 7380)
+++ trunk/psLib/src/math/psPolynomialUtils.c	(revision 7380)
@@ -0,0 +1,173 @@
+#include <stdio.h>
+#include "psMemory.h"
+#include "psError.h"
+#include "psLogMsg.h"
+#include "psVector.h"
+#include "psImage.h"
+#include "psBinaryOp.h"
+#include "psStats.h"
+#include "psConstants.h"
+#include "psPolynomial.h"
+#include "psMinimizePolyFit.h"
+#include "psCoord.h"
+#include "psPolynomialUtils.h"
+
+psPolynomial4D *psVectorChiClipFitPolynomial4D(
+    psPolynomial4D *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)
+{
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, 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, x, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(fErr, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+
+    // clipping range defined by min and max and/or clipSigma
+    float minClipSigma;
+    float maxClipSigma;
+    if (isfinite(stats->max)) {
+        maxClipSigma = +fabs(stats->max);
+    } else {
+        maxClipSigma = +fabs(stats->clipSigma);
+    }
+    if (isfinite(stats->min)) {
+        minClipSigma = -fabs(stats->min);
+    } else {
+        minClipSigma = -fabs(stats->clipSigma);
+    }
+    psVector *fit   = NULL;
+    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    // eventual expansion: user supplies one of various stats option pairs,
+    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
+    // evaluate the clipping sigma
+    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < stats->clipIter; N++) {
+        int Nkeep = 0;
+
+        poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
+        fit = psPolynomial4DEvalVector (poly, x, y, z, t);
+        resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit);
+
+        stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
+        psTrace (__func__, 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+        // set mask if pts are not valid
+        // we are masking out any point which is out of range
+        // recovery is not allowed with this scheme
+        for (int i = 0; i < resid->n; i++) {
+            if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
+                continue;
+            }
+            float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev);
+            if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            Nkeep ++;
+        }
+
+        psTrace (__func__, 4, "keeping %d of %d pts for fit\n",
+                 Nkeep, x->n);
+
+        stats->clippedNvalues = Nkeep;
+        psFree (fit);
+    }
+    // Free local temporary variables
+    psFree (resid);
+
+    if (poly == NULL) {
+        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
+        return(NULL);
+    }
+    return(poly);
+}
+
+psPolynomial2D *psImageBicubeFit(const psImage *image, long x, long y)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_INT_WITHIN_RANGE(x, 0, image->numCols - 1, NULL);
+    PS_ASSERT_INT_WITHIN_RANGE(y, 0, image->numRows - 1, NULL);
+
+    long ix = x - image->col0;
+    long iy = y - image->row0;
+
+    psF32 *Fm = &image->data.F32[iy - 1][ix];
+    psF32 *Fo = &image->data.F32[iy + 0][ix];
+    psF32 *Fp = &image->data.F32[iy + 1][ix];
+
+    double Fxm = Fm[-1] + Fo[-1] + Fp[-1];
+    double Fxp = Fm[+1] + Fo[+1] + Fp[+1];
+    double Fym = Fm[-1] + Fm[+0] + Fm[+1];
+    double Fyp = Fp[-1] + Fp[+0] + Fp[+1];
+    double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1];
+
+    psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    poly->mask[2][2] = 1;
+    poly->mask[1][2] = 1;
+    poly->mask[2][1] = 1;
+
+    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
+
+    poly->coeff[1][0] = (Fxp - Fxm)/6.0;
+    poly->coeff[0][1] = (Fyp - Fym)/6.0;
+
+    poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0;
+    poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0;
+
+    poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0;
+
+    return poly;
+}
+
+psPlane psImageBicubeMin(const psPolynomial2D *poly)
+{
+    psPlane min = { NAN, NAN, 0, 0 };   // Minimum value to return
+
+    PS_ASSERT_PTR_NON_NULL(poly, min);
+    PS_ASSERT_INT_EQUAL(poly->nX, 2, min);
+    PS_ASSERT_INT_EQUAL(poly->nY, 2, min);
+
+    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]); // Determinant
+    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
+    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
+
+    return min;
+}
Index: trunk/psLib/src/math/psPolynomialUtils.h
===================================================================
--- trunk/psLib/src/math/psPolynomialUtils.h	(revision 7380)
+++ trunk/psLib/src/math/psPolynomialUtils.h	(revision 7380)
@@ -0,0 +1,38 @@
+/** @file  psPolynomialUtils.h
+ *
+ * extra psPolynomial-related functions
+ */
+
+#ifndef PS_POLYNOMIAL_UTILS_H
+#define PS_POLYNOMIAL_UTILS_H
+
+#include "psVector.h"
+#include "psImage.h"
+#include "psCoord.h"
+#include "psStats.h"
+#include "psPolynomial.h"
+
+// perform vector clip-fit based on significance of deviations
+psPolynomial4D *psVectorChiClipFitPolynomial4D(
+    psPolynomial4D *poly,               // Polynomial to fit
+    psStats *stats,                     // Statistics to use in clipping
+    const psVector *mask,               // Mask for input values
+    psMaskType maskValue,               // Mask value
+    const psVector *f,                  // Value of the function, f(x,y,z,t)
+    const psVector *fErr,               // Error in the value
+    const psVector *x,                  // x ordinate
+    const psVector *y,                  // y ordinate
+    const psVector *z,                  // z ordinate
+    const psVector *t                   // t ordinate
+);
+
+// fit a 2D 2nd order polynomial to the 9 pixels centered on (x,y)
+psPolynomial2D *psImageBicubeFit(const psImage *image, // Image to fit
+                                 long x, long y // Pixels of centre of fit
+                                );
+
+// detemine the min(max) of the special 2D 2nd order polynomial
+psPlane psImageBicubeMin(const psPolynomial2D *poly // The polynomial
+                        );
+
+#endif /* PS_POLY_UTILS_H */
Index: trunk/psLib/src/math/psRegion.c
===================================================================
--- trunk/psLib/src/math/psRegion.c	(revision 7375)
+++ trunk/psLib/src/math/psRegion.c	(revision 7380)
@@ -84,2 +84,7 @@
 }
 
+bool inline psRegionIsNaN(const psRegion region)
+{
+    return isnan(region.x0) || isnan(region.x1) || isnan(region.y0) || isnan(region.y1);
+}
+
Index: trunk/psLib/src/math/psRegion.h
===================================================================
--- trunk/psLib/src/math/psRegion.h	(revision 7375)
+++ trunk/psLib/src/math/psRegion.h	(revision 7380)
@@ -70,4 +70,7 @@
 );
 
+// Test if any element of the region is NaN
+bool inline psRegionIsNaN(const psRegion region// Region to check
+                         );
 
 #endif
Index: trunk/psLib/src/math/psSparse.c
===================================================================
--- trunk/psLib/src/math/psSparse.c	(revision 7380)
+++ trunk/psLib/src/math/psSparse.c	(revision 7380)
@@ -0,0 +1,216 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "psMemory.h"
+#include "psError.h"
+#include "psLogMsg.h"
+#include "psVector.h"
+#include "psConstants.h"
+#include "psSparse.h"
+
+
+#define BUFFER 100                      // Size to increment at each go
+
+static void sparseFree(psSparse *sparse)
+{
+    if (!sparse) {
+        return;
+    }
+    psFree(sparse->Aij);
+    psFree(sparse->Bfj);
+    psFree(sparse->Qii);
+    psFree(sparse->Si);
+    psFree(sparse->Sj);
+    return;
+}
+
+// allocate a sparse matrix container for Nrows, with Nelem slots allocated
+psSparse *psSparseAlloc(long Nrows, long Nelem)
+{
+    psSparse *sparse = (psSparse *)psAlloc(sizeof(psSparse));
+    psMemSetDeallocator(sparse, (psFreeFunc)sparseFree);
+
+    sparse->Aij = psVectorAlloc(Nelem, PS_DATA_F32);
+    sparse->Si  = psVectorAlloc(Nelem, PS_DATA_S32);
+    sparse->Sj  = psVectorAlloc(Nelem, PS_DATA_S32);
+
+    sparse->Aij->n = 0;
+    sparse->Si->n  = 0;
+    sparse->Sj->n  = 0;
+    sparse->Nelem = 0;
+
+    sparse->Bfj = psVectorAlloc(Nrows, PS_DATA_F32);
+    sparse->Qii = psVectorAlloc(Nrows, PS_DATA_F32);
+    sparse->Bfj->n = Nrows;
+    sparse->Qii->n = Nrows;
+
+    sparse->Nrows = Nrows;
+
+    return sparse;
+}
+
+// user should only add elements above the diagonal, but we don't check this
+bool psSparseMatrixElement(psSparse *sparse, long i, long j, float value)
+{
+    PS_ASSERT_PTR_NON_NULL(sparse, false);
+    PS_ASSERT_INT_NONNEGATIVE(i, false);
+    PS_ASSERT_INT_NONNEGATIVE(j, false);
+
+    if (i < j) {
+        psLogMsg(__func__, PS_LOG_WARN, "i=%ld, j=%ld refers to a sub-diagonal element; values switched.\n");
+        long temp = i;
+        i = j;
+        j = temp;
+    }
+
+    if (i == j) {
+        // add to the diagonal
+        sparse->Qii->data.F32[i] = value;
+
+        // check vectors lengths and extend if needed
+        if (sparse->Nelem >= sparse->Aij->nalloc) {
+            psVectorRealloc(sparse->Aij, sparse->Aij->nalloc + BUFFER);
+            psVectorRealloc(sparse->Si,  sparse->Si->nalloc + BUFFER);
+            psVectorRealloc(sparse->Sj,  sparse->Sj->nalloc + BUFFER);
+        }
+
+        long k = sparse->Nelem;         // Index at which to add
+        sparse->Aij->data.F32[k] = value;
+        sparse->Si->data.S32[k]  = i;
+        sparse->Sj->data.S32[k]  = j;
+
+        sparse->Nelem ++;
+        sparse->Aij->n ++;
+        sparse->Si->n ++;
+        sparse->Sj->n ++;
+    } else {
+        // check vectors lengths and extend if needed
+        if (sparse->Nelem >= sparse->Aij->nalloc - 1) {
+            psVectorRealloc(sparse->Aij, sparse->Aij->nalloc + BUFFER);
+            psVectorRealloc(sparse->Si,  sparse->Si->nalloc + BUFFER);
+            psVectorRealloc(sparse->Sj,  sparse->Sj->nalloc + BUFFER);
+        }
+
+        long k = sparse->Nelem;         // Index at which to add
+        sparse->Aij->data.F32[k] = value;
+        sparse->Si->data.S32[k]  = i;
+        sparse->Sj->data.S32[k]  = j;
+        k++;
+
+        sparse->Aij->data.F32[k] = value;
+        sparse->Si->data.S32[k]  = j;
+        sparse->Sj->data.S32[k]  = i;
+
+        sparse->Nelem  += 2;
+        sparse->Aij->n += 2;
+        sparse->Si->n  += 2;
+        sparse->Sj->n  += 2;
+    }
+
+    return true;
+}
+
+void inline psSparseVectorElement(psSparse *sparse, long i, float value)
+{
+
+    sparse->Bfj->data.F32[i] = value;
+    return;
+}
+
+// multiply A * x
+psVector *psSparseMatrixTimesVector(psVector *output, const psSparse *matrix, const psVector *vector)
+{
+    PS_ASSERT_PTR_NON_NULL(matrix, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(vector, NULL);
+
+    output = psVectorRecycle(output, vector->n, PS_TYPE_F32);
+
+    long Nelem = 0;                     // Number of elements
+    for (long j = 0; j < vector->n; j++) {
+        double F = 0;                    // Running total
+        while (matrix->Sj->data.S32[Nelem] == j) {
+            long i = matrix->Si->data.S32[Nelem];
+            F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem];
+            Nelem++;
+        }
+        output->data.F32[j] = F;
+    }
+    return output;
+}
+
+psVector *psSparseSolve(psVector *output, psSparseConstraint constraint, const psSparse *sparse, int Niter)
+{
+    PS_ASSERT_PTR_NON_NULL(sparse, NULL);
+    PS_ASSERT_INT_POSITIVE(Niter, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(constraint.paramDelta, 0.0, NULL);
+
+    // Dereference some vectors
+    psVector *Qii = sparse->Qii;
+    psVector *Bfj = sparse->Bfj;
+
+    output = psVectorCopy(output, Bfj, PS_DATA_F32);
+
+    // temporary storage for intermediate results
+    psVector *dQ = psVectorAlloc(output->n, PS_DATA_F32);
+    dQ->n = output->n;
+
+    for (long j = 0; j < Niter; j++) {
+        dQ = psSparseMatrixTimesVector(dQ, sparse, output);
+        for (long i = 0; i < dQ->n; i++) {
+            psF32 dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
+            if (fabs (dG) > constraint.paramDelta) {
+                if (dG > 0) {
+                    dG = +constraint.paramDelta;
+                } else {
+                    dG = -constraint.paramDelta;
+                }
+            }
+            output->data.F32[i] -= dG;
+            output->data.F32[i] = PS_MAX(output->data.F32[i], constraint.paramMin);
+            output->data.F32[i] = PS_MIN(output->data.F32[i], constraint.paramMax);
+        }
+    }
+    psFree(dQ);
+
+    return output;
+}
+
+bool psSparseResort(psSparse *sparse)
+{
+    long Nelem = sparse->Nelem;
+
+    psVector *index = psVectorSortIndex(NULL, sparse->Sj); // Index key for sorting
+    if (!index) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to sort sparse matrix.\n");
+        return false;
+    }
+    psVector *Aij = sparse->Aij;
+    psVector *Si = sparse->Si;
+    psVector *Sj = sparse->Sj;
+
+    // allocate new temporary vectors
+    psVector *tAij = psVectorAlloc(Nelem, PS_DATA_F32);
+    psVector *tSi  = psVectorAlloc(Nelem, PS_DATA_S32);
+    psVector *tSj  = psVectorAlloc(Nelem, PS_DATA_S32);
+    tAij->n = tSi->n = tSj->n = Nelem;
+
+    for (long i = 0; i < Nelem; i++) {
+        long j = index->data.U32[i];
+        tAij->data.F32[i] = Aij->data.F32[j];
+        tSi->data.S32[i]  = Si->data.S32[j];
+        tSj->data.S32[i]  = Sj->data.S32[j];
+    }
+    psFree(index);
+    psFree(Aij);
+    psFree(Si);
+    psFree(Sj);
+
+    sparse->Aij = tAij;
+    sparse->Si = tSi;
+    sparse->Sj = tSj;
+
+    return true;
+}
Index: trunk/psLib/src/math/psSparse.h
===================================================================
--- trunk/psLib/src/math/psSparse.h	(revision 7380)
+++ trunk/psLib/src/math/psSparse.h	(revision 7380)
@@ -0,0 +1,67 @@
+/** @file  psSparse.h
+ *
+ * functions to manipulate sparse matrices equations
+ *
+ */
+
+#ifndef PS_SPARSE_H
+#define PS_SPARSE_H
+
+// constraints to limit the range of the matrix equation solution
+typedef struct
+{
+    double paramDelta;
+    double paramMin;
+    double paramMax;
+}
+psSparseConstraint;
+
+// A sparse matrix equation: A x = Bf
+typedef struct
+{
+    psVector *Aij;                      // Aij contains the populated elements of the matrix
+    psVector *Bfj;                      // Bfj contains the elements of the vector Bf
+    psVector *Qii;                      // Qii contains the diagonal elements of Aij
+    psVector *Si;                       // Si contains the i-index values of Aij
+    psVector *Sj;                       // Sj contains the j-index values of Aij
+    long Nelem;                         // Number of elements
+    long Nrows;                         // Number of rows
+}
+psSparse;
+
+// allocate a sparse matrix structure
+psSparse *psSparseAlloc(long Nrows, long Nelem);
+
+// add a new matrix element
+// user should only add elements above the diagonal
+bool psSparseMatrixElement(psSparse *sparse, // Matrix to which to add
+                           long i, long j, // Matrix indices at which to add
+                           float value  // Value to add
+                          );
+
+// define a new sparse matrix equation vector element
+void inline psSparseVectorElement(psSparse *sparse, // Matrix to which to add
+                                  long i,      // Index to add
+                                  float value  // Value to add
+                                 );
+
+// perform the operation matrix * vector on a sparse matrix and a vector
+psVector *psSparseMatrixTimesVector(psVector *output, // Output vector, or NULL
+                                    const psSparse *matrix, // Sparse matrix
+                                    const psVector *vector // Corresponding vector
+                                   );
+
+// re-sort a sparse matrix to have all elements in index order rather than insertion order
+// call this before solving, but after populating matrix and vector
+bool psSparseResort(psSparse *sparse    // Matrix to re-sort
+                   );
+
+// solve the equation A x = Bf for the value of x
+// a good starting guess is the vector Bf
+psVector *psSparseSolve(psVector *output,// The output vector, or NULL
+                        psSparseConstraint constraint, // Constraint to limit the range of the solution
+                        const psSparse *sparse, // Sparse matrix
+                        int Niter       // Number of iterations
+                       );
+
+#endif /* PS_SPARSE_H */
Index: trunk/psLib/src/math/psVectorSmooth.c
===================================================================
--- trunk/psLib/src/math/psVectorSmooth.c	(revision 7380)
+++ trunk/psLib/src/math/psVectorSmooth.c	(revision 7380)
@@ -0,0 +1,103 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "psMemory.h"
+#include "psError.h"
+#include "psLogMsg.h"
+#include "psVector.h"
+#include "psVectorSmooth.h"
+#include "psConstants.h"
+
+
+psVector *psVectorSmooth(psVector *output,
+                         const psVector *input,
+                         double sigma,
+                         double Nsigma
+                        )
+{
+    PS_ASSERT_VECTOR_NON_NULL(input, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(sigma, 0.0, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(Nsigma, 0.0, NULL);
+
+    // relevant terms
+    long Nrange = sigma*Nsigma + 0.5;   // Extent of smoothing
+    long Npixel = 2*Nrange + 1;         // Total size of smoothing kernel
+    long Nbin = input->n;               // Number of elements
+    double factor = -0.5/(sigma*sigma); // Factor for Gaussian
+
+    #define VECTOR_SMOOTH_CASE(TYPE) \
+case PS_TYPE_##TYPE: { \
+        output = psVectorRecycle(output, Nbin, PS_TYPE_##TYPE); \
+        /* generate normalized gaussian */ \
+        psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_##TYPE); \
+        double sum = 0.0; \
+        for (long i = -Nrange; i < Nrange + 1; i++) { \
+            gaussnorm->data.TYPE[i+Nrange] = exp(factor*i*i); \
+            sum += gaussnorm->data.TYPE[i+Nrange]; \
+        } \
+        for (long i = -Nrange; i < Nrange + 1; i++) { \
+            gaussnorm->data.TYPE[i+Nrange] /= sum; \
+        } \
+        ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \
+        \
+        /* smooth vector */ \
+        psVector *temp = psVectorAlloc(Nbin, PS_TYPE_##TYPE); \
+        ps##TYPE *vi = input->data.TYPE; \
+        ps##TYPE *vo = temp->data.TYPE; \
+        /* smooth first Nrange pixels, with renorm */ \
+        /* XXX need to check that this does not run over end for narrow vectors */ \
+        for (long i = 0; i < Nrange; i++, vi++, vo++) { \
+            ps##TYPE *vr = vi - i; \
+            ps##TYPE *vg = gauss - i; \
+            double g = 0; \
+            double s = 0; \
+            for (int n = -i; n < Nrange + 1; n++, vr++, vg++) { \
+                s += *vg * *vr; \
+                g += *vg; \
+            } \
+            *vo = s / g; \
+        } \
+        /* smooth middle pixels */ \
+        for (long i = Nrange; i < Nbin - Nrange; i++, vi++, vo++) { \
+            ps##TYPE *vr = vi - Nrange; \
+            ps##TYPE *vg = gauss - Nrange; \
+            double s = 0; \
+            for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { \
+                s += *vg * *vr; \
+            } \
+            *vo = s; \
+        } \
+        /* smooth last Nrange pixels, with renorm */ \
+        /* XXX does this miss the last column? */ \
+        for (long i = Nbin - Nrange; i < Nbin; i++, vi++, vo++) { \
+            ps##TYPE *vr = vi - Nrange; \
+            ps##TYPE *vg = gauss - Nrange; \
+            double g = 0; \
+            double s = 0; \
+            for (int n = -Nrange; n < Nbin - i - 1; n++, vr++, vg++) { \
+                s += *vg * *vr; \
+                g += *vg; \
+            } \
+            *vo = s / g; \
+        } \
+        memcpy(output->data.TYPE, temp->data.TYPE, Nbin*sizeof(ps##TYPE)); \
+        psFree(temp); \
+        psFree(gaussnorm); \
+        break; \
+    }
+
+    switch (input->type.type) {
+        VECTOR_SMOOTH_CASE(F32);
+        VECTOR_SMOOTH_CASE(F64);
+    default: {
+            char* typeStr;
+            PS_TYPE_NAME(typeStr, input->type.type);
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type %s is not valid.", typeStr);
+            return NULL;
+        }
+    }
+    return output;
+}
Index: trunk/psLib/src/math/psVectorSmooth.h
===================================================================
--- trunk/psLib/src/math/psVectorSmooth.h	(revision 7380)
+++ trunk/psLib/src/math/psVectorSmooth.h	(revision 7380)
@@ -0,0 +1,17 @@
+/** @file  psVectorSmooth.h
+ *
+ * smooth the input vector
+ *
+ */
+
+#ifndef PS_VECTOR_SMOOTH_H
+#define PS_VECTOR_SMOOTH_H
+
+// Smooth a vector with a Gaussian
+psVector *psVectorSmooth(psVector *output, // Output vector, or NULL
+                         const psVector *input, // Input vector (F32 or F64 only)
+                         double sigma,  // Gausian width (standard deviations)
+                         double Nsigma  // Number of standard deviations for Gaussian to extend (either side)
+                        );
+
+#endif