Index: /branches/rel10_ifa/psModules/src/pslib/Makefile.am
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/Makefile.am	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/Makefile.am	(revision 6530)
@@ -0,0 +1,18 @@
+noinst_LTLIBRARIES = libpsmodulepslib.la
+
+libpsmodulepslib_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS)
+libpsmodulepslib_la_LDFLAGS  = -release $(PACKAGE_VERSION)
+libpsmodulepslib_la_SOURCES  = \
+    psAdditionals.c \
+    psImageJpeg.c \
+    psLine.c \
+    psPolynomialUtils.c \
+    psSparse.c
+
+psmoduleincludedir = $(includedir)
+psmoduleinclude_HEADERS = \
+    psAdditionals.h \
+    psImageJpeg.h \
+    psLine.h \
+    psPolynomialUtils.h \
+    psSparse.h
Index: /branches/rel10_ifa/psModules/src/pslib/psAdditionals.c
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psAdditionals.c	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psAdditionals.c	(revision 6530)
@@ -0,0 +1,202 @@
+#include <stdio.h>
+#include <strings.h>
+#include "pslib.h"
+
+#include "psAdditionals.h"
+
+
+psMetadata *pap_psMetadataCopy(psMetadata *out,
+                               const psMetadata *in)
+{
+    PS_ASSERT_PTR_NON_NULL(in,NULL);
+    if (out ==  NULL) {
+        out = psMetadataAlloc();
+    }
+    psMetadataItem *inItem = NULL;
+    psMetadataIterator *iter = psMetadataIteratorAlloc(*(psMetadata**)&in, PS_LIST_HEAD, NULL);
+    unsigned long numPointers = 0;      // Number of pointers we were forced to copy
+    while ((inItem = psMetadataGetAndIncrement(iter))) {
+        // Need to look for MULTI, which won't be picked up using the iterator.
+        psMetadataItem *multiCheckItem = psMetadataLookup(in, inItem->name);
+        unsigned int flag = PS_META_REPLACE; // Flag to indicate MULTI; otherwise, replace
+        if (multiCheckItem->type == PS_DATA_METADATA_MULTI) {
+            psTrace(__func__, 10, "MULTI: %s (%s)\n", inItem->name, inItem->comment);
+            flag = PS_DATA_METADATA_MULTI;
+        }
+
+        psTrace(__func__, 5, "Copying %s (%s)...\n", inItem->name, inItem->comment);
+
+        #define PS_METADATA_COPY_CASE(NAME,TYPE) \
+    case PS_TYPE_##NAME: \
+        if (! psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_TYPE_##NAME | flag, inItem->comment, \
+                            inItem->data.TYPE)) { \
+            psErrorStackPrint(stderr, "Error copying %s (%s) in the metadata\n", inItem->name, \
+                              inItem->comment); \
+        } \
+        break;
+
+        switch (inItem->type) {
+            // Numerical types
+            PS_METADATA_COPY_CASE(BOOL,B);
+            PS_METADATA_COPY_CASE(S8,S8);
+            PS_METADATA_COPY_CASE(S16,S16);
+            PS_METADATA_COPY_CASE(S32,S32);
+            PS_METADATA_COPY_CASE(U8,U8);
+            PS_METADATA_COPY_CASE(U16,U16);
+            PS_METADATA_COPY_CASE(U32,U32);
+            PS_METADATA_COPY_CASE(F32,F32);
+            PS_METADATA_COPY_CASE(F64,F64);
+
+            // String: relying on the fact that this will copy the string, not point at it.
+        case PS_DATA_STRING:
+            psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_STRING | flag, inItem->comment,
+                          inItem->data.V);
+            break;
+
+            // Metadata: copy the next level and stuff that in too
+        case PS_DATA_METADATA: {
+                psMetadata *metadata = pap_psMetadataCopy(NULL, inItem->data.md);
+                psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_METADATA | flag, inItem->comment,
+                              metadata);
+                break;
+            }
+            // Other kinds of pointers
+        default:
+            numPointers++;
+            psTrace(__func__, 10, "Copying a pointer in the metadata: %x\n", inItem->type);
+            psMetadataAdd(out, PS_LIST_TAIL, inItem->name, inItem->type | flag, inItem->comment,
+                          inItem->data.V);
+            break;
+        }
+    }
+    psFree(iter);
+
+    if (numPointers > 0) {
+        psLogMsg(__func__, PS_LOG_WARN, "Forced to copy %d pointers when copying metadata.  Updating the "
+                 "copied psMetadata will affect the original!\n", numPointers);
+    }
+    return out;
+}
+
+// may need to extend this to change the keyname in the copy
+bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key)
+{
+
+    psMetadataItem *item = psMetadataLookup (in, key);
+    if (item == NULL)
+        return false;
+
+    psMetadataAddItem (out, item, PS_LIST_TAIL, PS_META_REPLACE);
+    return true;
+}
+
+void psMetadataPrint(psMetadata *md, int level)
+{
+    psMetadataIterator *iter = psMetadataIteratorAlloc(md, PS_LIST_HEAD, NULL); // Iterator
+    psMetadataItem *item = NULL;        // Item from metadata
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        // Indent...
+        for (int i = 0; i < level; i++) {
+            printf(" ");
+        }
+        printf("%s", item->name);
+        if (item->comment && strlen(item->comment) > 0) {
+            printf(" (%s)", item->comment);
+        }
+        printf(": ");
+        switch (item->type) {
+        case PS_DATA_STRING:
+            printf("%s", (char*)item->data.V);
+            break;
+        case PS_DATA_BOOL:
+            if (item->data.B) {
+                printf("True");
+            } else {
+                printf("False");
+            }
+            break;
+        case PS_DATA_S32:
+            printf("%d", item->data.S32);
+            break;
+        case PS_DATA_F32:
+            printf("%f", item->data.F32);
+            break;
+        case PS_DATA_F64:
+            printf("%f", item->data.F64);
+            break;
+        case PS_DATA_METADATA:
+            printf("\n");
+            psMetadataPrint(item->data.V, level + 1);
+            break;
+        default:
+            printf("\n");
+            psError(PS_ERR_IO, false, "Non-printable metadata type: %x\n", item->type);
+        }
+        printf("\n");
+    }
+    psFree(iter);
+
+    return;
+}
+
+// XXX: This should probably be implemented using strpbrk
+psList *psStringSplit(const char *string,
+                      const char *splitters)
+{
+    psList *values = psListAlloc(NULL); // The list of values to return
+    unsigned int length = strlen(string); // The length of the string
+    unsigned int numSplitters = strlen(splitters); // Number of characters that might split
+    unsigned int start = 0;             // The position of the start of a word
+    for (int i = 1; i < length; i++) {
+        bool split = false;             // Is this character a splitter?
+        for (int j = 0; j < numSplitters && ! split; j++) {
+            if (string[i] == splitters[j]) {
+                split = true;
+            }
+        }
+        if (split) {
+            if (i == start) {
+                // Some idiot put in two spaces, or two commas or something
+                start++;
+            } else {
+                // We're at the end of the word
+                psString word = psStringNCopy(&string[start], i - start);
+                (void)psListAdd(values, PS_LIST_TAIL, word);
+                start = i + 1;
+                psFree(word);
+            }
+        }
+    }
+    if (start < length) {
+        // Copy the last word
+        psString word = psStringNCopy(&string[start], length - start);
+        (void)psListAdd(values, PS_LIST_TAIL, word);
+        psFree(word);
+    }
+
+    return values;
+}
+
+#ifndef whitespace
+#define whitespace(c) (((c) == ' ') || ((c) == '\t'))
+#endif
+
+/* Strip whitespace from the start and end of STRING. */
+int psStringStrip (char *string)
+{
+
+    int i;
+
+    if (string == (char *) NULL)
+        return (FALSE);
+
+    for (i = 0; whitespace (string[i]); i++)
+        ;
+    if (i)
+        memmove (string, string + i, strlen(string+i)+1);
+    for (i = strlen (string) - 1; (i > 0) && whitespace (string[i]); i--)
+        ;
+    string[++i] = 0;
+    return (i);
+
+}
Index: /branches/rel10_ifa/psModules/src/pslib/psAdditionals.h
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psAdditionals.h	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psAdditionals.h	(revision 6530)
@@ -0,0 +1,35 @@
+// Functions that should go into psLib.
+
+
+#ifndef PS_ADDITIONALS_H
+#define PS_ADDITIONALS_H
+
+#include "pslib.h"
+
+// Deep copy of metadata
+// Corrected version of MHPCC code in psLib at the moment
+psMetadata *pap_psMetadataCopy(psMetadata *out, // Target, to which the copy is made
+                               const psMetadata *in // Source, from which the copy is made
+                              );
+
+// copy a metadata item from one psMetadata to another
+bool psMetadataItemTransfer (psMetadata *out, // Destination: copy is placed here
+                             psMetadata *in,  // Source: item comes from here
+                             char *key        // key to identify the metadata item
+                            );
+
+
+// Print out the metadata
+void psMetadataPrint(psMetadata *md,    // Metadata to print
+                     int level          // Indent level
+                    );
+
+// Split string on given characters
+psList *psStringSplit(const char *string, // String to split
+                      const char *splitters // Characters on which to split
+                     );
+
+// strip whitespace from head and tail of string
+int          psStringStrip (char *string);
+
+#endif
Index: /branches/rel10_ifa/psModules/src/pslib/psImageJpeg.c
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psImageJpeg.c	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psImageJpeg.c	(revision 6530)
@@ -0,0 +1,191 @@
+# include <stdio.h>
+# include <strings.h>  // for strcasecmp
+# include <unistd.h>   // for unlink
+# include <pslib.h>
+# include "jpeglib.h"
+# include "psImageJpeg.h"
+
+static void psImageJpegColormapFree (psImageJpegColormap *map)
+{
+
+    if (map == NULL)
+        return;
+
+    psFree (map->red);
+    psFree (map->green);
+    psFree (map->blue);
+    return;
+}
+
+psImageJpegColormap *psImageJpegColormapAlloc ()
+{
+
+    psImageJpegColormap *map;
+    map = psAlloc (sizeof(psImageJpegColormap));
+    psMemSetDeallocator(map, (psFreeFunc) psImageJpegColormapFree);
+
+    map->red   = psVectorAlloc (256, PS_TYPE_U8);
+    map->blue  = psVectorAlloc (256, PS_TYPE_U8);
+    map->green = psVectorAlloc (256, PS_TYPE_U8);
+
+    return (map);
+}
+
+psImageColormap *psImageJpegColormapSet (psImageJpegColormap *map, char *name)
+{
+
+    if (map == NULL) {
+        map = psImageJpegColormapAlloc ();
+    }
+
+    /* grayscale */
+    if ((!strcasecmp (name, "grayscale")) || (!strcasecmp (name, "greyscale"))) {
+        for (int i = 0; i < map->red->n; i++) {
+            map->red->data.U8[i]   = PS_JPEG_RANGELIM(i);
+            map->green->data.U8[i] = PS_JPEG_RANGELIM(i);
+            map->blue->data.U8[i]  = PS_JPEG_RANGELIM(i);
+        }
+        return map;
+    }
+
+    /* -grayscale */
+    if ((!strcasecmp (name, "-grayscale")) || (!strcasecmp (name, "-greyscale"))) {
+        for (int i = 0; i < map->red->n; i++) {
+            map->red->data.U8[i]   = PS_JPEG_RANGELIM(256 - i);
+            map->green->data.U8[i] = PS_JPEG_RANGELIM(256 - i);
+            map->blue->data.U8[i]  = PS_JPEG_RANGELIM(256 - i);
+        }
+        return map;
+    }
+
+    /* rainbow */
+    if (!strcasecmp (name, "rainbow")) {
+        int I1 = 0.25*map->red->n;
+        int I2 = 0.50*map->red->n;
+        int I3 = 0.75*map->red->n;
+        for (int i = 0; i < I1; i++) {
+            map->red->data.U8[i]   = 0;
+            map->green->data.U8[i] = 0;
+            map->blue->data.U8[i]  = PS_JPEG_RANGELIM(4*i);
+        }
+        for (int i = I1; i < I2; i++) {
+            map->red->data.U8[i]   = PS_JPEG_RANGELIM(4*(i - I1));
+            map->green->data.U8[i] = 0;
+            map->blue->data.U8[i]  = PS_JPEG_RANGELIM(4*(I2 - i));
+        }
+        for (int i = I2; i < I3; i++) {
+            map->red->data.U8[i]   = 255;
+            map->green->data.U8[i] = 4*(i - I2);
+            map->blue->data.U8[i]  = 0;
+        }
+        for (int i = I3; i < map->red->n; i++) {
+            map->red->data.U8[i]   = 255;
+            map->green->data.U8[i] = 255;
+            map->blue->data.U8[i]  = PS_JPEG_RANGELIM(4*(i - I3));
+        }
+        return map;
+    }
+
+    /* heat */
+    if (!strcasecmp (name, "heat")) {
+        int I1 = 0.25*map->red->n;
+        int I2 = 0.50*map->red->n;
+        int I3 = 0.75*map->red->n;
+        for (int i = 0; i < I1; i++) {
+            map->red->data.U8[i]   = PS_JPEG_RANGELIM(2*i);
+            map->green->data.U8[i] = 0;
+            map->blue->data.U8[i]  = 0;
+        }
+        for (int i = I1; i < I2; i++) {
+            map->red->data.U8[i]   = PS_JPEG_RANGELIM(2*i);
+            map->green->data.U8[i] = PS_JPEG_RANGELIM(2*(i - I1));
+            map->blue->data.U8[i]  = 0;
+        }
+        for (int i = I2; i < I3; i++) {
+            map->red->data.U8[i]   = 255;
+            map->green->data.U8[i] = PS_JPEG_RANGELIM(2*(i - I1));
+            map->blue->data.U8[i]  = PS_JPEG_RANGELIM(2*(i - I2));
+        }
+        for (int i = I3; i < map->red->n; i++) {
+            map->red->data.U8[i]   = 255;
+            map->green->data.U8[i] = 255;
+            map->blue->data.U8[i]  = PS_JPEG_RANGELIM(2*(i - I2));
+        }
+        return map;
+    }
+
+    // invalid colormap : warn user
+    map = psImageJpegColormapSet (map, "greyscale");
+    return map;
+}
+
+bool psImageJpeg (psImageJpegColormap *map, psImage *image, char *filename, float min, float max)
+{
+
+    struct jpeg_compress_struct cinfo;
+    struct jpeg_error_mgr jerr;
+
+    int pixel;
+    JSAMPLE *jpegLine;   // Points to data for current line
+    JSAMPROW jpegLineList[1];  // pointer to JSAMPLE row[s]
+    JSAMPLE *outPix;
+
+    if (map->red == NULL)
+        return false;
+    if (map->green == NULL)
+        return false;
+    if (map->blue == NULL)
+        return false;
+
+    /* JPEG init calls */
+    cinfo.err = jpeg_std_error (&jerr);
+    jpeg_create_compress (&cinfo);
+
+    /* open file, prep for jpeg */
+    FILE *f = fopen (filename, "w");
+    if (f == NULL) {
+        fprintf (stderr, "failed to open %s for output\n", filename);
+        return (TRUE);
+    }
+    jpeg_stdio_dest(&cinfo, f);
+
+    /* set up color jpeg buffers */
+    int quality = 75;
+    cinfo.image_width = image->numCols; // image width and height, in pixels
+    cinfo.image_height = image->numRows;
+    cinfo.input_components = 3;
+    cinfo.in_color_space = JCS_RGB;
+    jpeg_set_defaults (&cinfo);
+    jpeg_set_quality (&cinfo, quality, true); // limit to baseline-JPEG values
+    jpeg_start_compress (&cinfo, true);
+
+    jpegLine = psAlloc (3*image->numCols*sizeof(JSAMPLE));
+    jpegLineList[0] = jpegLine;
+
+    psU8 *Rpix = map->red->data.U8;
+    psU8 *Gpix = map->green->data.U8;
+    psU8 *Bpix = map->blue->data.U8;
+
+    float zero = min;
+    float scale = 256.0/(max - min);
+
+    for (int j = 0; j < image->numRows; j++) {
+        psF32 *row = image->data.F32[j];
+
+        outPix = jpegLine;
+        for (int i = 0; i < image->numCols; i++, outPix += 3) {
+            pixel = PS_JPEG_SCALEVALUE(row[i],zero,scale);
+            outPix[0] = Rpix[pixel];
+            outPix[1] = Gpix[pixel];
+            outPix[2] = Bpix[pixel];
+        }
+        jpeg_write_scanlines (&cinfo, jpegLineList, 1);
+    }
+
+    jpeg_finish_compress (&cinfo);
+    fclose (f);
+    jpeg_destroy_compress (&cinfo);
+
+    psFree (jpegLine);
+    return true;
+}
Index: /branches/rel10_ifa/psModules/src/pslib/psImageJpeg.h
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psImageJpeg.h	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psImageJpeg.h	(revision 6530)
@@ -0,0 +1,33 @@
+/** @file  psImageJpeg.h
+ *
+ * the functions to generate JPEG images from psImage 
+ */
+
+# ifndef PS_IMAGE_JPEG_H
+# define PS_IMAGE_JPEG_H
+
+typedef struct
+{
+    psString name;
+    psVector *red;
+    psVector *green;
+    psVector *blue;
+}
+psImageJpegColormap;
+
+# define PS_JPEG_RANGELIM(A)(PS_MAX(0,PS_MIN(255,(A))))
+
+# define PS_JPEG_SCALEVALUE(VALUE,ZERO,SCALE)(PS_MAX(0,PS_MIN(255,(SCALE*(VALUE-ZERO)))))
+
+// allocate a colormap (does not define the map values)
+psImageJpegColormap *psImageJpegColormapAlloc ()
+{
+
+    // set the colormap values using the supplied name
+    psImageColormap *psImageJpegColormapSet (psImageJpegColormap *map, char *name) {
+
+        // write out a JPEG file using the supplied image and colormap
+        // output goes to the specified filename
+        bool psImageJpeg (psImageJpegColormap *map, psImage *image, char *filename, float min, float max) {
+
+            # endif /* PS_IMAGE_JPEG_H */
Index: /branches/rel10_ifa/psModules/src/pslib/psLine.c
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psLine.c	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psLine.c	(revision 6530)
@@ -0,0 +1,55 @@
+# include <pslib.h>
+# include "psLine.h"
+
+static void psLineFree (psLine *line)
+{
+
+    if (line == NULL)
+        return;
+
+    psFree (line->line);
+    return;
+}
+
+// allocate a psLine structrue
+psLine *psLineAlloc (int Nline)
+{
+
+    psLine *line;
+    line = psAlloc (sizeof(psLine));
+    psMemSetDeallocator(line, (psFreeFunc) psLineFree);
+
+    line->Nline = 0;
+    line->NLINE = Nline;
+    line->line = psAlloc (Nline);
+    return (line);
+}
+
+bool psLineInit (psLine *line)
+{
+    if (line == NULL)
+        return (false);
+    line->Nline = 0;
+    return (true);
+}
+
+bool psLineAdd (psLine *line, char *format, ...)
+{
+
+    int Nchar;
+    va_list ap;
+
+    if (line == NULL)
+        return (false);
+
+    int nMax = line->NLINE - line->Nline;
+
+    va_start (ap, format);
+    Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap);
+    line->Nline += PS_MIN (nMax - 1, Nchar);
+    va_end (ap);
+
+    if (Nchar >= nMax)
+        return (false);
+    return (true);
+}
Index: /branches/rel10_ifa/psModules/src/pslib/psLine.h
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psLine.h	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psLine.h	(revision 6530)
@@ -0,0 +1,27 @@
+/** @file  psLine.h
+ *
+ * the psLine functions allow manipulation of fixed-length lines
+ */
+
+# ifndef PS_LINE_H
+# define PS_LINE_H
+
+// structure to carry a dynamic string
+typedef struct
+{
+    int NLINE;
+    int Nline;
+    char *line;
+}
+psLine;
+
+// allocate a line object of length Nline
+psLine      *psLineAlloc (int Nline);
+
+// (re-)init the line
+bool      psLineInit (psLine *line);
+
+// add the new string segment to the line
+bool      psLineAdd (psLine *line, char *format, ...);
+
+# endif /* PS_LINE_H */
Index: /branches/rel10_ifa/psModules/src/pslib/psPolynomialUtils.c
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psPolynomialUtils.c	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psPolynomialUtils.c	(revision 6530)
@@ -0,0 +1,161 @@
+# include <stdio.h>
+# include <strings.h>  // for strcasecmp
+# include <unistd.h>   // for unlink
+# include <pslib.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 (psImage *image, int x, int y)
+{
+
+    int ix = x - image->col0;
+    int 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 (psPolynomial2D *poly)
+{
+
+    psPlane min;
+
+    min.xErr = min.yErr = 0;
+
+    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]);
+
+    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: /branches/rel10_ifa/psModules/src/pslib/psPolynomialUtils.h
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psPolynomialUtils.h	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psPolynomialUtils.h	(revision 6530)
@@ -0,0 +1,28 @@
+/** @file  psPolynomialUtils.h
+ *
+ * extra psPolynomial-related functions
+ */
+
+# ifndef PS_POLY_UTILS_H
+# define PS_POLY_UTILS_H
+
+// perform vector clip-fit based on significance of deviations
+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);
+
+// fit a 2D 2nd order polynomial to the 9 pixels centered on (x,y)
+psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
+
+// detemine the min(max) of the special 2D 2nd order polynomial
+psPlane psImageBicubeMin (psPolynomial2D *poly);
+
+# endif /* PS_POLY_UTILS_H */
Index: /branches/rel10_ifa/psModules/src/pslib/psSparse.c
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psSparse.c	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psSparse.c	(revision 6530)
@@ -0,0 +1,245 @@
+# include <pslib.h>
+# include "psSparse.h"
+
+static void psSparseFree (psSparse *sparse)
+{
+    if (sparse == NULL)
+        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 (int Nrows, int Nelem)
+{
+
+    psSparse *sparse = (psSparse *) psAlloc (sizeof(psSparse));
+    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->Nrows = Nrows;
+    psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree);
+    return (sparse);
+}
+
+// user should only add elements above the diagonal, but we don't check this
+void psSparseMatrixElement (psSparse *sparse, int i, int j, float value)
+{
+
+    int k;
+
+    if (i < j) {
+        fprintf (stderr, "*** error: subdiagonal element ***\n");
+        return;
+    }
+
+    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 + 100);
+            psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
+            psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
+        }
+
+        k = sparse->Nelem;
+        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 + 100);
+            psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
+            psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
+        }
+
+        k = sparse->Nelem;
+        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;
+}
+
+void psSparseVectorElement (psSparse *sparse, int i, float value)
+{
+
+    sparse->Bfj->data.F32[i] = value;
+    return;
+}
+
+// multiple A * x
+psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector)
+{
+
+    int i, Nelem;
+    float F;
+
+    if (output == NULL) {
+        output = psVectorAlloc (vector->n, PS_DATA_F32);
+    }
+
+    Nelem = 0;
+    for (int j = 0; j < vector->n; j++) {
+        F = 0;
+        while (matrix->Sj->data.S32[Nelem] == j) {
+            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 *guess, psSparseConstraint constraint, psSparse *sparse, int Niter)
+{
+
+    psF32 dG;
+
+    psVector *Qii = sparse->Qii;
+    psVector *Bfj = sparse->Bfj;
+
+    guess = psVectorCopy (guess, Bfj, PS_DATA_F32);
+
+    // temporary storage for intermediate results
+    psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32);
+
+    for (int j = 0; j < Niter; j++) {
+        dQ = psSparseMatrixTimesVector (dQ, sparse, guess);
+        for (int i = 0; i < dQ->n; i++) {
+            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;
+                }
+            }
+            guess->data.F32[i] -= dG;
+            guess->data.F32[i] = PS_MAX (guess->data.F32[i], constraint.paramMin);
+            guess->data.F32[i] = PS_MIN (guess->data.F32[i], constraint.paramMax);
+        }
+    }
+    psFree (dQ);
+    return (guess);
+}
+
+void psSparseResort (psSparse *sparse)
+{
+
+    int Nelem = sparse->Nelem;
+
+    psVector *index = psVectorSortIndex (NULL, sparse->Sj);
+    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);
+    for (int i = 0; i < Nelem; i++) {
+        int 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;
+}
+
+void psSparseMatrixTest ()
+{
+
+    // build a sparse matrix
+    psSparse *sparse = psSparseAlloc (3, 9);
+
+    psSparseMatrixElement (sparse, 0, 0, 3.0);
+    psSparseMatrixElement (sparse, 1, 1, 2.0);
+    psSparseMatrixElement (sparse, 2, 2, 1.0);
+
+    psSparseMatrixElement (sparse, 1, 0, 0.1);
+    psSparseMatrixElement (sparse, 2, 0, -0.1);
+
+    psSparseResort (sparse);
+    for (int i = 0; i < sparse->Nelem; i++) {
+        fprintf (stderr, "%d %d %f\n",
+                 sparse->Si->data.S32[i],
+                 sparse->Sj->data.S32[i],
+                 sparse->Aij->data.F32[i]);
+    }
+
+    psVector *x = psVectorAlloc (3, PS_DATA_F32);
+    x->data.F32[0] = 3;
+    x->data.F32[1] = 5;
+    x->data.F32[2] = 7;
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    psVector *B = psSparseMatrixTimesVector (NULL, sparse, x);
+    fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]);
+
+    sparse->Bfj->data.F32[0] = B->data.F32[0];
+    sparse->Bfj->data.F32[1] = B->data.F32[1];
+    sparse->Bfj->data.F32[2] = B->data.F32[2];
+
+    psSparseConstraint constraint;
+    constraint.paramMin   = -1e8;
+    constraint.paramMax   = +1e8;
+    constraint.paramDelta = +1e8;
+
+    x = psSparseSolve (x, constraint, sparse, 0);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 1);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 2);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 3);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 4);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+    return;
+}
+
Index: /branches/rel10_ifa/psModules/src/pslib/psSparse.h
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psSparse.h	(revision 6530)
+++ /branches/rel10_ifa/psModules/src/pslib/psSparse.h	(revision 6530)
@@ -0,0 +1,62 @@
+
+/** @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
+// Aij contains the populated elements of the matrix
+// Bfj contains the elements of the vector Bf
+// Qii contains the diagonal elements of Aij
+// Si contains the i-index values of Aij
+// Sj contains the j-index values of Aij
+typedef struct
+{
+    psVector *Aij;
+    psVector *Bfj;
+    psVector *Qii;
+    psVector *Si;
+    psVector *Sj;
+    int Nelem;
+    int Nrows;
+}
+psSparse;
+
+// allocate a sparse matrix structure
+psSparse *psSparseAlloc (int Nrows, int Nelem);
+
+// add a new matrix element
+// user should only add elements above the diagonal
+void psSparseMatrixElement (psSparse *sparse, int i, int j, float value);
+
+// define a new sparse matrix equation vector element
+void psSparseVectorElement (psSparse *sparse, int i, float value);
+
+// perform the operation matrix * vector on a sparse matrix and a vector
+psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *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
+void psSparseResort (psSparse *sparse);
+
+// solve the equation A x = Bf for the value of x
+// a good starting guess is the vector Bf
+psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter);
+
+// test of the sparse matrix solutions (move to test suite)
+void psSparseMatrixTest ();
+
+# endif /* PS_SPARSE_H */