Index: /trunk/psLib/configure.ac
===================================================================
--- /trunk/psLib/configure.ac	(revision 7379)
+++ /trunk/psLib/configure.ac	(revision 7380)
@@ -26,5 +26,5 @@
 
 SRCPATH='${top_srcdir}/src'
-SRCDIRS="sys astro db fft fits imageops math mathtypes types xml"
+SRCDIRS="sys astro db fft fits imageops jpeg math mathtypes types xml"
 # escape two escapes at this level so \\ gets passed to the shell and \ to perl
 SRCINC=`echo "${SRCDIRS=}" | ${PERL} -pe "s|(\w+)|-I\\\\${SRCPATH=}/\1|g"`
@@ -241,4 +241,8 @@
 AC_SUBST([GSL_CFLAGS])
 
+dnl libjpeg -------------------------------------------------------------------
+
+AC_CHECK_LIB(jpeg,jpeg_CreateCompress,[],[AC_MSG_ERROR([jpeg library not found.])])
+
 dnl ------------------- XML2 options ---------------------
 AC_ARG_WITH(xml2-config,
@@ -316,5 +320,4 @@
 AM_CONDITIONAL(DOXYGEN, test x$doxygen = xtrue)
 
-
 dnl ------- restore CFLAGS / LDFLAGS (tests done) --------
 CFLAGS="${CFLAGS} -Wall -Werror"
@@ -338,4 +341,5 @@
   src/fits/Makefile
   src/imageops/Makefile
+  src/jpeg/Makefile
   src/math/Makefile
   src/mathtypes/Makefile
@@ -349,4 +353,5 @@
   test/fits/Makefile
   test/imageops/Makefile
+  test/jpeg/Makefile
   test/math/Makefile
   test/mathtypes/Makefile
Index: /trunk/psLib/src/imageops/Makefile.am
===================================================================
--- /trunk/psLib/src/imageops/Makefile.am	(revision 7379)
+++ /trunk/psLib/src/imageops/Makefile.am	(revision 7380)
@@ -5,4 +5,5 @@
 libpslibimageops_la_CPPFLAGS = $(SRCINC)
 libpslibimageops_la_SOURCES = \
+	psImageBackground.c \
 	psImageConvolve.c \
 	psImageGeomManip.c \
@@ -11,5 +12,6 @@
 	psImageStats.c \
 	psImageStructManip.c \
-    psImageMaskOps.c
+	psImageMaskOps.c \
+	psImageUnbin.c
 
 EXTRA_DIST = imageops.i
@@ -17,4 +19,5 @@
 pslibincludedir = $(includedir)
 pslibinclude_HEADERS = \
+	psImageBackground.h \
 	psImageConvolve.h \
 	psImageGeomManip.h \
@@ -23,5 +26,6 @@
 	psImageStats.h \
 	psImageStructManip.h \
-    psImageMaskOps.h
+	psImageMaskOps.h \
+	psImageUnbin.h
 
 CLEANFILES = *~
Index: /trunk/psLib/src/imageops/psImageBackground.c
===================================================================
--- /trunk/psLib/src/imageops/psImageBackground.c	(revision 7380)
+++ /trunk/psLib/src/imageops/psImageBackground.c	(revision 7380)
@@ -0,0 +1,87 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "psMemory.h"
+#include "psImage.h"
+#include "psVector.h"
+#include "psStats.h"
+#include "psType.h"
+#include "psConstants.h"
+#include "psRandom.h"
+#include "psError.h"
+
+
+psStats *psImageBackground(const psImage *image, const psImage *mask, psMaskType maskValue,
+                           double fmin, double fmax, long nMax, psRandom *rng)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    if (mask) {
+        PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, NULL);
+    }
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(fmin, 0.0, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(fmax, 0.0, NULL);
+    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(fmin, 1.0, NULL);
+    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(fmax, 1.0, NULL);
+    PS_ASSERT_INT_POSITIVE(nMax, NULL);
+
+    // Size of image
+    long nx = image->numCols;
+    long ny = image->numRows;
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE); // Statistics, for return
+
+    int Nsubset = PS_MIN(nMax, nx*ny);  // Number of pixels in nubset
+    int Npixels = nx*ny;                // Total number of pixels
+
+    psVector *values = psVectorAlloc(Nsubset, PS_TYPE_F32); // Vector containing subsample
+
+    // Minimum and maximum values
+    float min = values->data.F32[0];
+    float max = values->data.F32[0];
+
+    long n = 0;                         // Number of actual pixels in subset
+    for (long i = 0; i < Nsubset; i++) {
+        double frnd = psRandomUniform(rng);
+        int pixel = Npixels * frnd;
+        int ix = pixel % nx;
+        int iy = pixel / nx;
+
+        if (mask && mask->data.U8[iy][ix] & maskValue) {
+            continue;
+        }
+
+        float value = image->data.F32[iy][ix];
+        min = PS_MIN(value, min);
+        max = PS_MIN(value, max);
+        values->data.F32[n] = value;
+        n++;
+    }
+    values->n = n;
+
+    int imin = fmin * n;
+    int imax = fmax * n;
+    int npts = imax - imin + 1;
+
+    if (!psVectorSort(values, values)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
+        psFree(stats);
+        psFree(values);
+        return NULL;
+    }
+
+    float value = 0;
+    for (long i = imin; (i <= imax) && (i < n); i++) {
+        value += values->data.F32[i];
+    }
+    value = value / npts;
+
+    stats->robustMedian = value;
+    stats->robustUQ = values->data.F32[imax];
+    stats->robustLQ = values->data.F32[imin];
+
+    psFree(values);
+    return stats;
+    // XXX correct for selection bias??
+}
Index: /trunk/psLib/src/imageops/psImageBackground.h
===================================================================
--- /trunk/psLib/src/imageops/psImageBackground.h	(revision 7380)
+++ /trunk/psLib/src/imageops/psImageBackground.h	(revision 7380)
@@ -0,0 +1,18 @@
+#ifndef PS_IMAGE_BACKGROUND_H
+#define PS_IMAGE_BACKGROUND_H
+
+#include "psStats.h"
+#include "psImage.h"
+#include "psType.h"
+
+// Get the background for an image
+psStats *psImageBackground(const psImage *image, // Image for which to get the background
+                           const psImage *mask, // Mask image
+                           psMaskType maskValue, // Mask pixels which this mask value
+                           double fmin, // Fraction to return in the lower quartile field (0.25 for LQ)
+                           double fmax, // Fraction to return in the upper quartile field (0.75 for LQ)
+                           long nMax,   // Maximum number of pixels to subsample
+                           psRandom *rng // Random number generator (for pixel selection)
+                          );
+
+#endif // #ifndef PS_IMAGE_BACKGROUND_H
Index: /trunk/psLib/src/imageops/psImageGeomManip.c
===================================================================
--- /trunk/psLib/src/imageops/psImageGeomManip.c	(revision 7379)
+++ /trunk/psLib/src/imageops/psImageGeomManip.c	(revision 7380)
@@ -10,6 +10,6 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-06 22:55:18 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-07 03:22:06 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -23,4 +23,5 @@
 #include "psImageGeomManip.h"
 
+#include "psAbort.h"
 #include "psError.h"
 #include "psImage.h"
@@ -878,2 +879,89 @@
 }
 
+
+#define FLIP_X_CASE(TYPENAME,TYPE) \
+case TYPENAME: { \
+    long numRows = input->numRows; \
+    long numCols = input->numCols; \
+    for (long i = 0; i < numRows; i++) { \
+        for (long j = 0; j < numCols; j++) { \
+            output->data.TYPE[i][j] = input->data.TYPE[i][numCols - j - 1]; \
+        } \
+    } \
+    break; \
+}
+
+#define FLIP_Y_CASE(TYPENAME,TYPE) \
+case TYPENAME: { \
+    long numRows = input->numRows; \
+    long numCols = input->numCols; \
+    for (long i = 0; i < numRows; i++) { \
+        for (long j = 0; j < numCols; j++) { \
+            output->data.TYPE[i][j] = input->data.TYPE[numRows - i - 1][j]; \
+        } \
+    } \
+    break; \
+}
+
+psImage *psImageFlip(psImage *output, const psImage *input, bool xFlip, bool yFlip)
+{
+    PS_ASSERT_IMAGE_NON_NULL(input, NULL);
+
+    if (xFlip && yFlip) {
+        // This is equivalent to a 180 degree rotation;
+        return psImageRotate(output, input, M_PI, NAN, PS_INTERPOLATE_BILINEAR);
+    }
+
+    if (!xFlip && !yFlip) {
+        // They want something, so let's give it to them
+        return psImageCopy(output, input, input->type.type);
+    }
+
+    output = psImageRecycle(output, input->numCols, input->numRows, input->type.type);
+
+    if (xFlip) {
+        switch (input->type.type) {
+            FLIP_X_CASE(PS_TYPE_U8,  U8);
+            FLIP_X_CASE(PS_TYPE_U16, U16);
+            FLIP_X_CASE(PS_TYPE_U32, U32);
+            FLIP_X_CASE(PS_TYPE_U64, U64);
+            FLIP_X_CASE(PS_TYPE_S8,  S8);
+            FLIP_X_CASE(PS_TYPE_S16, S16);
+            FLIP_X_CASE(PS_TYPE_S32, S32);
+            FLIP_X_CASE(PS_TYPE_S64, S64);
+            FLIP_X_CASE(PS_TYPE_F32, F32);
+            FLIP_X_CASE(PS_TYPE_F64, F64);
+            FLIP_X_CASE(PS_TYPE_C32, C32);
+            FLIP_X_CASE(PS_TYPE_C64, C64);
+        default:
+            psFree(output);
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown type for input image: %x\n", input->type.type);
+            return NULL;
+        }
+        return output;
+    }
+    if (yFlip) {
+        switch (input->type.type) {
+            FLIP_Y_CASE(PS_TYPE_U8,  U8);
+            FLIP_Y_CASE(PS_TYPE_U16, U16);
+            FLIP_Y_CASE(PS_TYPE_U32, U32);
+            FLIP_Y_CASE(PS_TYPE_U64, U64);
+            FLIP_Y_CASE(PS_TYPE_S8,  S8);
+            FLIP_Y_CASE(PS_TYPE_S16, S16);
+            FLIP_Y_CASE(PS_TYPE_S32, S32);
+            FLIP_Y_CASE(PS_TYPE_S64, S64);
+            FLIP_Y_CASE(PS_TYPE_F32, F32);
+            FLIP_Y_CASE(PS_TYPE_F64, F64);
+            FLIP_Y_CASE(PS_TYPE_C32, C32);
+            FLIP_Y_CASE(PS_TYPE_C64, C64);
+        default:
+            psFree(output);
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown type for input image: %x\n", input->type.type);
+            return NULL;
+        }
+        return output;
+    }
+
+    psAbort(__func__, "Should never get here.\n");
+    return NULL;
+}
Index: /trunk/psLib/src/imageops/psImageGeomManip.h
===================================================================
--- /trunk/psLib/src/imageops/psImageGeomManip.h	(revision 7379)
+++ /trunk/psLib/src/imageops/psImageGeomManip.h	(revision 7380)
@@ -8,6 +8,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-01 02:43:57 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-07 03:22:06 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -151,9 +151,16 @@
     psRegion region,                   ///< the size of the transformed image
     const psPixels* pixels,            /**< if not NULL, consists of psPixelCoords and specifies
-                                            * which pixels in output image shall be transformed;
-                                            * otherwise, entire image generated*/
+                                                * which pixels in output image shall be transformed;
+                                                * otherwise, entire image generated*/
     psImageInterpolateMode mode,       ///< the interpolation scheme to be used
     double exposedValue                ///< Exposed value to which non-corresponding pixels are set
 );
 
+// Flip the input image
+psImage *psImageFlip(psImage *output,   // Output image, or NULL
+                     const psImage *input, // Input image
+                     bool xFlip,        // Flip x axis?
+                     bool yFlip         // Flip y axis?
+                    );
+
 #endif // #ifndef PS_IMAGE_GEOM_MANIP_H
Index: /trunk/psLib/src/imageops/psImagePixelManip.c
===================================================================
--- /trunk/psLib/src/imageops/psImagePixelManip.c	(revision 7379)
+++ /trunk/psLib/src/imageops/psImagePixelManip.c	(revision 7380)
@@ -10,6 +10,6 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-05-16 03:27:13 $
+ *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-07 03:22:06 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -413,3 +413,2 @@
     return numClipped;
 }
-
Index: /trunk/psLib/src/jpeg/.cvsignore
===================================================================
--- /trunk/psLib/src/jpeg/.cvsignore	(revision 7380)
+++ /trunk/psLib/src/jpeg/.cvsignore	(revision 7380)
@@ -0,0 +1,7 @@
+.deps
+.libs
+Makefile
+Makefile.in
+*.la
+*.lo
+*.loT
Index: /trunk/psLib/src/jpeg/Makefile.am
===================================================================
--- /trunk/psLib/src/jpeg/Makefile.am	(revision 7380)
+++ /trunk/psLib/src/jpeg/Makefile.am	(revision 7380)
@@ -0,0 +1,13 @@
+#Makefile for image jpeg functions of psLib
+#
+noinst_LTLIBRARIES = libpslibjpeg.la
+
+libpslibjpeg_la_CPPFLAGS = $(SRCINC)
+libpslibjpeg_la_SOURCES = \
+	psImageJpeg.c
+
+pslibincludedir = $(includedir)
+pslibinclude_HEADERS = \
+	psImageJpeg.h
+
+CLEANFILES = *~
Index: /trunk/psLib/src/jpeg/psImageJpeg.c
===================================================================
--- /trunk/psLib/src/jpeg/psImageJpeg.c	(revision 7380)
+++ /trunk/psLib/src/jpeg/psImageJpeg.c	(revision 7380)
@@ -0,0 +1,198 @@
+#include <stdio.h>
+#include <strings.h>  // for strcasecmp
+#include <string.h>
+
+#include <jpeglib.h>
+
+#include "psMemory.h"
+#include "psImage.h"
+#include "psVector.h"
+#include "psError.h"
+#include "psConstants.h"
+#include "psImageJpeg.h"
+
+static void imageJpegColormapFree(psImageJpegColormap *map)
+{
+
+    if (!map) {
+        return;
+    }
+
+    psFree(map->red);
+    psFree(map->green);
+    psFree(map->blue);
+    return;
+}
+
+psImageJpegColormap *psImageJpegColormapAlloc ()
+{
+    psImageJpegColormap *map = psAlloc(sizeof(psImageJpegColormap));
+    psMemSetDeallocator(map, (psFreeFunc)imageJpegColormapFree);
+
+    map->red   = psVectorAlloc(256, PS_TYPE_U8);
+    map->blue  = psVectorAlloc(256, PS_TYPE_U8);
+    map->green = psVectorAlloc(256, PS_TYPE_U8);
+
+    return map;
+}
+
+psImageJpegColormap *psImageJpegColormapSet(psImageJpegColormap *map, const char *name)
+{
+
+    if (!map) {
+        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
+    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Invalid colormap name: %s --- using greyscale\n", name);
+    return psImageJpegColormapSet (map, "greyscale");
+}
+
+bool psImageJpeg(const psImageJpegColormap *map, const psImage *image, const char *filename,
+                 float min, float max)
+{
+    PS_ASSERT_PTR_NON_NULL(map, false);
+    PS_ASSERT_IMAGE_NON_NULL(image, false);
+    PS_ASSERT_VECTOR_NON_NULL(map->red, false);
+    PS_ASSERT_VECTOR_NON_NULL(map->green, false);
+    PS_ASSERT_VECTOR_NON_NULL(map->blue, false);
+    PS_ASSERT_PTR_NON_NULL(filename, false);
+    PS_ASSERT_INT_POSITIVE(strlen(filename), false);
+
+    struct jpeg_compress_struct cinfo;
+    struct jpeg_error_mgr jerr;
+
+    long pixel;
+    JSAMPLE *jpegLine;   // Points to data for current line
+    JSAMPROW jpegLineList[1];  // pointer to JSAMPLE row[s]
+    JSAMPLE *outPix;
+
+    /* 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) {
+        psError(PS_ERR_IO, true, "failed to open %s for output\n", filename);
+        return false;
+    }
+    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: /trunk/psLib/src/jpeg/psImageJpeg.h
===================================================================
--- /trunk/psLib/src/jpeg/psImageJpeg.h	(revision 7380)
+++ /trunk/psLib/src/jpeg/psImageJpeg.h	(revision 7380)
@@ -0,0 +1,39 @@
+/** @file  psImageJpeg.h
+ *
+ * the functions to generate JPEG images from psImage
+ */
+
+#ifndef PS_IMAGE_JPEG_H
+#define PS_IMAGE_JPEG_H
+
+#include "psImage.h"
+
+typedef struct
+{
+    psVector *red;                      // Red colormap
+    psVector *green;                    // Green colormap
+    psVector *blue;                     // Blue colormap
+}
+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
+psImageJpegColormap *psImageJpegColormapSet(psImageJpegColormap *map, // Colormap to set
+        const char *name // Name of colormap
+                                           );
+
+// write out a JPEG file using the supplied image and colormap
+// output goes to the specified filename
+bool psImageJpeg(const psImageJpegColormap *map, // Color map
+                 const psImage *image,  // Image to write
+                 const char *filename,  // Filename of JPEG
+                 float min, float max   // Minimum and maximum values
+                );
+
+#endif /* PS_IMAGE_JPEG_H */
Index: /trunk/psLib/src/math/Makefile.am
===================================================================
--- /trunk/psLib/src/math/Makefile.am	(revision 7379)
+++ /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 7379)
+++ /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 7379)
+++ /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
Index: /trunk/psLib/src/pslib_strict.h
===================================================================
--- /trunk/psLib/src/pslib_strict.h	(revision 7379)
+++ /trunk/psLib/src/pslib_strict.h	(revision 7380)
@@ -9,6 +9,6 @@
 *  @author Eric Van Alst, MHPCC
 *
-*  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-05-25 21:19:45 $
+*  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-06-07 03:22:06 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -42,4 +42,5 @@
 
 #include "psXML.h"
+
 #include "psImageConvolve.h"
 #include "psImageGeomManip.h"
@@ -49,4 +50,7 @@
 #include "psImageStructManip.h"
 #include "psImageMaskOps.h"
+#include "psImageUnbin.h"
+
+#include "psImageJpeg.h"
 
 #include "psBinaryOp.h"
@@ -61,4 +65,5 @@
 #include "psRegionForImage.h"
 #include "psPolynomial.h"
+#include "psPolynomialUtils.h"
 #include "psSpline.h"
 #include "psStats.h"
@@ -84,6 +89,12 @@
 #include "psMetadata.h"
 #include "psMetadataConfig.h"
+#include "psMetadataItemParse.h"
+#include "psMetadataItemCompare.h"
 #include "psPixels.h"
 #include "psArguments.h"
+#include "psVectorSmooth.h"
+#include "psImageBackground.h"
+#include "psEllipse.h"
+#include "psSparse.h"
 
 #endif // #ifndef PS_LIB_STRICT_H
Index: /trunk/psLib/src/sys/Makefile.am
===================================================================
--- /trunk/psLib/src/sys/Makefile.am	(revision 7379)
+++ /trunk/psLib/src/sys/Makefile.am	(revision 7380)
@@ -9,4 +9,5 @@
 	psError.c      \
 	psErrorCodes.c \
+	psLine.c       \
 	psLogMsg.c     \
 	psMemory.c     \
@@ -32,4 +33,5 @@
 	psError.h      \
 	psErrorCodes.h \
+	psLine.h       \
 	psLogMsg.h     \
 	psMemory.h     \
Index: /trunk/psLib/src/sys/psLine.c
===================================================================
--- /trunk/psLib/src/sys/psLine.c	(revision 7380)
+++ /trunk/psLib/src/sys/psLine.c	(revision 7380)
@@ -0,0 +1,59 @@
+#include <stdio.h>
+#include <stdarg.h>
+
+#include "psMemory.h"
+#include "psConstants.h"
+#include "psLine.h"
+
+static void lineFree(psLine *line)
+{
+
+    if (!line) {
+        return;
+    }
+
+    psFree (line->line);
+    return;
+}
+
+// allocate a psLine structrue
+psLine *psLineAlloc(long Nline)
+{
+    psLine *line = psAlloc(sizeof(psLine));
+    psMemSetDeallocator(line, (psFreeFunc)lineFree);
+    line->Nline = 0;
+    line->NLINE = Nline;
+    line->line = psAlloc(Nline);
+
+    return line;
+}
+
+bool psLineInit(psLine *line)
+{
+    if (!line) {
+        return false;
+    }
+
+    line->Nline = 0;
+    return true;
+}
+
+bool psLineAdd(psLine *line, const char *format, ...)
+{
+    if (!line) {
+        return false;
+    }
+
+    long nMax = line->NLINE - line->Nline;
+
+    va_list ap;
+    va_start(ap, format);
+    long 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: /trunk/psLib/src/sys/psLine.h
===================================================================
--- /trunk/psLib/src/sys/psLine.h	(revision 7380)
+++ /trunk/psLib/src/sys/psLine.h	(revision 7380)
@@ -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
+{
+    long NLINE;                 // allocated length
+    long Nline;                 // current length
+    char *line;                 // character string data
+}
+psLine;
+
+// allocate a line object of length Nline
+psLine *psLineAlloc(long Nline);
+
+// (re-)init the line
+bool psLineInit(psLine *line);
+
+// add the new string segment to the line
+bool psLineAdd(psLine *line, const char *format, ...);
+
+#endif /* PS_LINE_H */
Index: /trunk/psLib/src/sys/psString.c
===================================================================
--- /trunk/psLib/src/sys/psString.c	(revision 7379)
+++ /trunk/psLib/src/sys/psString.c	(revision 7380)
@@ -13,6 +13,6 @@
  *  @author David Robbins, MHPCC
  *
- *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-05-31 20:56:37 $
+ *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-07 03:22:06 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -31,9 +31,9 @@
 psString psStringCopy(const char *string)
 {
-    PS_ASSERT_PTR_NON_NULL(string, NULL);
+    // Pass through NULL values!
+
     // Allocate memory using psAlloc function
     // Copy input string to memory just allocated
     // Return the copy
-    // Pass through NULL values
     return string ? strcpy(psAlloc(strlen(string) + 1), string) : NULL;
 }
@@ -218,18 +218,15 @@
 // given the input string, search for all copies of the key, and replace with the replacement value
 // the input string may be freed if not needed
-char *psStringSubstitute (char *input, char *replace, char *key)
-{
-
-    char *p;
-
-    if (key == NULL)
+char *psStringSubstitute(char *input, const char *replace, const char *key)
+{
+    if (key == NULL || strlen(key) == 0) {
         return input;
-    if (strlen(key) == 0)
-        return input;
+    }
 
     while (true) {
-        p = strstr (input, key);
-        if (p == NULL)
+        char *p = strstr (input, key);
+        if (!p) {
             return input;
+        }
 
         // we have input = xxxkeyxxx
@@ -241,15 +238,14 @@
 
         // copy the first segement into 'output'
-        strncpy (output, input, Nc);
+        strncpy(output, input, Nc);
 
         // copy the key replacement to the start of the key
-
-        strcpy (&output[Nc], replace);
+        strcpy(&output[Nc], replace);
         Nc += strlen (replace);
 
         // copy the remainder to the end of the replacement
-        strcpy (&output[Nc], p + strlen(key));
-
-        psFree (input);
+        strcpy(&output[Nc], p + strlen(key));
+
+        psFree(input);
         input = output;
     }
@@ -257,2 +253,35 @@
 }
 
+psArray *psStringSplitArray(const char *string, const char *splitters, bool multi)
+{
+
+    psList *list = psStringSplit(string, splitters, multi);
+    psArray *array = psListToArray(list);
+    psFree (list);
+    return array;
+}
+
+#ifndef whitespace
+#define whitespace(c) (((c) == ' ') || ((c) == '\t'))
+#endif
+
+/* Strip whitespace from the start and end of STRING. */
+size_t psStringStrip(char *string)
+{
+    long i;
+
+    if (!string || strlen(string) == 0) {
+        return 0;
+    }
+
+    for (i = 0; i < strlen(string) && whitespace(string[i]); i++)
+        ; // No action
+    if (i) {
+        memmove (string, string + i, strlen(string+i)+1);
+    }
+    for (i = strlen (string) - 1; (i > 0) && whitespace(string[i]); i--)
+        ; // No action
+    string[++i] = 0;
+
+    return i;
+}
Index: /trunk/psLib/src/sys/psString.h
===================================================================
--- /trunk/psLib/src/sys/psString.h	(revision 7379)
+++ /trunk/psLib/src/sys/psString.h	(revision 7380)
@@ -14,6 +14,6 @@
  *  @author David Robbins, MHPCC
  *
- *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-06-02 21:33:34 $
+ *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-07 03:22:06 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -107,11 +107,20 @@
 );
 
-// given the input string, search for all copies of the key, and replace with the replacement value
+// Same as psStringSplit, but return result as an array
+psArray *psStringSplitArray(const char *string, // String to split
+                            const char *splitters, // Characters on which to split
+                            bool multipleAreSignificant // Are multiple occurences significant?
+                           );
+
+// Given the input string, search for all copies of the key, and replace with the replacement value
 // the input string may be freed if not needed
-char *psStringSubstitute (
-    char *input,    ///< input string to be modified
-    char *replace,    ///< replacement value
-    char *key    ///< string to be replaced in input
-);
+char *psStringSubstitute (char *input,  ///< input string to be modified
+                          const char *replace, ///< replacement value
+                          const char *key ///< string to be replaced in input
+                         );
+
+// strip whitespace from head and tail of string
+size_t psStringStrip(char *string);
+
 
 /** @} */// Doxygen - End of SystemGroup Functions
Index: /trunk/psLib/src/types/Makefile.am
===================================================================
--- /trunk/psLib/src/types/Makefile.am	(revision 7379)
+++ /trunk/psLib/src/types/Makefile.am	(revision 7380)
@@ -12,4 +12,6 @@
 	psMetadata.c \
 	psMetadataConfig.c \
+	psMetadataItemParse.c \
+	psMetadataItemCompare.c \
 	psPixels.c \
 	psArguments.c
@@ -26,4 +28,6 @@
 	psMetadata.h \
 	psMetadataConfig.h \
+	psMetadataItemParse.h \
+	psMetadataItemCompare.h \
 	psPixels.h \
 	psArguments.h
Index: /trunk/psLib/src/types/psMetadata.c
===================================================================
--- /trunk/psLib/src/types/psMetadata.c	(revision 7379)
+++ /trunk/psLib/src/types/psMetadata.c	(revision 7380)
@@ -12,6 +12,6 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.107 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-06-06 22:54:22 $
+ *  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-07 03:22:06 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -445,4 +445,18 @@
     return out;
 }
+
+
+// may need to extend this to change the keyname in the copy
+bool psMetadataItemTransfer(psMetadata *out, const psMetadata *in, const char *key)
+{
+
+    psMetadataItem *item = psMetadataLookup(in, key);
+    if (!item) {
+        return false;
+    }
+
+    return psMetadataAddItem(out, item, PS_LIST_TAIL, PS_META_REPLACE);
+}
+
 
 bool psMetadataAddItem(psMetadata *md,
Index: /trunk/psLib/src/types/psMetadata.h
===================================================================
--- /trunk/psLib/src/types/psMetadata.h	(revision 7379)
+++ /trunk/psLib/src/types/psMetadata.h	(revision 7380)
@@ -11,6 +11,6 @@
 *  @author Ross Harman, MHPCC
 *
-*  @version $Revision: 1.76 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-06-06 22:54:22 $
+*  @version $Revision: 1.77 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-06-07 03:22:06 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -361,4 +361,10 @@
 );
 
+// Copy a metadata item from one psMetadata to another
+bool psMetadataItemTransfer(psMetadata *out, // Destination: copy is placed here
+                            const psMetadata *in, // Source: item comes from here
+                            const char *key // key to identify the metadata item
+                           );
+
 /** Add existing metadata item to metadata collection.
  *
Index: /trunk/psLib/src/types/psMetadataItemCompare.c
===================================================================
--- /trunk/psLib/src/types/psMetadataItemCompare.c	(revision 7380)
+++ /trunk/psLib/src/types/psMetadataItemCompare.c	(revision 7380)
@@ -0,0 +1,83 @@
+#include <stdio.h>
+#include <string.h>
+#include <strings.h>
+
+#include "psType.h"
+#include "psError.h"
+#include "psAbort.h"
+#include "psTrace.h"
+#include "psMetadata.h"
+#include "psMetadataItemCompare.h"
+
+bool psMetadataItemCompare(const psMetadataItem *compare, // Item to compare to the template
+                           const psMetadataItem *template // The template
+                          )
+{
+
+    // First order checks
+    if (! compare || ! template) {
+        return false;
+    }
+    if (strcmp(compare->name, template->name)
+            != 0) {
+        return false;
+    }
+
+
+    #define COMPARE_CASE(TEMPLATENAME, TEMPLATETYPE, COMPARENAME, COMPARETYPENAME) \
+case PS_TYPE_##COMPARETYPENAME: { \
+        TEMPLATETYPE value = (TEMPLATETYPE)compare->data.COMPARENAME; \
+        if (value != compare->data.COMPARENAME) { \
+            return false; \
+        } \
+        return (template->data.TEMPLATENAME == value) ? true : false; \
+    }
+
+    #define TEMPLATE_CASE(TYPENAME, NAME, TYPE) \
+case PS_TYPE_##TYPENAME: \
+    switch(compare->type) { \
+        COMPARE_CASE(NAME, TYPE, B  , BOOL); \
+        COMPARE_CASE(NAME, TYPE, U8 , U8 ); \
+        COMPARE_CASE(NAME, TYPE, U16, U16); \
+        COMPARE_CASE(NAME, TYPE, U32, U32); \
+        COMPARE_CASE(NAME, TYPE, S8 , S8 ); \
+        COMPARE_CASE(NAME, TYPE, S16, S16); \
+        COMPARE_CASE(NAME, TYPE, S32, S32); \
+        COMPARE_CASE(NAME, TYPE, F32, F32); \
+        COMPARE_CASE(NAME, TYPE, F64, F64); \
+    default: \
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Don't know how to compare types %x and %x\n", \
+                compare->type, template->type); \
+        return false; \
+    }
+
+
+    switch (template->type) {
+        TEMPLATE_CASE(BOOL, B  , psBool)
+        ;
+        TEMPLATE_CASE(U8,   U8 , psU8 );
+        TEMPLATE_CASE(U16,  U16, psU16);
+        TEMPLATE_CASE(U32,  U32, psU32);
+        TEMPLATE_CASE(S8 ,  S8 , psS8 );
+        TEMPLATE_CASE(S16,  S16, psS16);
+        TEMPLATE_CASE(S32,  S32, psS32);
+        TEMPLATE_CASE(F32,  F32, psF32);
+        TEMPLATE_CASE(F64,  F64, psF64);
+    case PS_DATA_STRING:
+        if (compare->type != PS_DATA_STRING) {
+            return false;
+        }
+        psTrace(__func__, 10, "Comparing '%s' with '%s'\n", compare->data.V, template->
+                data.V);
+        return (strcasecmp(compare->data.V, template->
+                           data.V) == 0) ? true : false;
+    default:
+        // Simply don't know how to compare more complex types.
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Don't know how to compare types %x and %x\n",
+                compare->type, template->type);
+        return false;
+    }
+
+    psAbort(__func__, "Should never get here.\n");
+    return false;
+}
Index: /trunk/psLib/src/types/psMetadataItemCompare.h
===================================================================
--- /trunk/psLib/src/types/psMetadataItemCompare.h	(revision 7380)
+++ /trunk/psLib/src/types/psMetadataItemCompare.h	(revision 7380)
@@ -0,0 +1,12 @@
+#ifndef PS_METADATA_ITEM_COMPARE_H
+#define PS_METADATA_ITEM_COMPARE_H
+
+#include "psMetadata.h"
+
+bool psMetadataItemCompare(const psMetadataItem *compare, // Item to compare to the template
+                           const psMetadataItem *template // The template
+                          )
+;
+
+
+#endif
Index: /trunk/psLib/src/types/psMetadataItemParse.c
===================================================================
--- /trunk/psLib/src/types/psMetadataItemParse.c	(revision 7380)
+++ /trunk/psLib/src/types/psMetadataItemParse.c	(revision 7380)
@@ -0,0 +1,277 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <strings.h>
+
+#include "psMemory.h"
+#include "psError.h"
+#include "psMetadata.h"
+#include "psMetadataItemParse.h"
+
+# define PS_METADATA_ITEM_PARSE_NUMBER(INTYPE,OUTTYPE) \
+case PS_DATA_##INTYPE: \
+return (ps##OUTTYPE)(item->data.INTYPE); \
+
+
+// NOTE: This function flows through so that errors may be handled by the "default" case.
+#define PS_METADATA_ITEM_PARSE_STRING_FLOAT(OUTTYPE,FUNCTION) \
+case PS_DATA_STRING: { \
+    char *end = NULL; \
+    ps##OUTTYPE number = FUNCTION(item->data.V, &end); \
+    if (end != item->data.V) { \
+        return number; \
+    } \
+}
+
+// NOTE: This function flows through so that errors may be handled by the "default" case.
+#define PS_METADATA_ITEM_PARSE_STRING_INT(OUTTYPE,FUNCTION) \
+case PS_DATA_STRING: { \
+    char *end = NULL; \
+    ps##OUTTYPE number = FUNCTION(item->data.V, &end, 10); \
+    if (end != item->data.V) { \
+        return number; \
+    } \
+}
+
+psBool psMetadataItemParseBool(const psMetadataItem *item
+                              )
+{
+    PS_ASSERT_PTR_NON_NULL(item, false);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, Bool);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, Bool);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  Bool);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, Bool);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, Bool);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  Bool);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, Bool);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, Bool);
+    case PS_DATA_STRING:
+        if (strcasecmp(item->data.V, "true") == 0) {
+            return true;
+        } else if (strcasecmp(item->data.V, "false") == 0) {
+            return false;
+        }
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of boolean type (%x) --- "
+                "treating as false.\n", item->name, item->comment, item->type);
+        return false;
+    }
+}
+
+psF32 psMetadataItemParseF32(const psMetadataItem *item
+                            )
+{
+    PS_ASSERT_PTR_NON_NULL(item, NAN);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, F32);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, F32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  F32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, F32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, F32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  F32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, F32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, F32);
+        PS_METADATA_ITEM_PARSE_STRING_FLOAT(F32, strtof);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of floating point type (%x) --- "
+                "treating as NaN.\n", item->name, item->comment, item->type);
+        return NAN;
+    }
+}
+
+psF64 psMetadataItemParseF64(const psMetadataItem *item
+                            )
+{
+    PS_ASSERT_PTR_NON_NULL(item, NAN);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, F64);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, F64);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  F64);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, F64);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, F64);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  F64);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, F64);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, F64);
+        PS_METADATA_ITEM_PARSE_STRING_FLOAT(F32, strtod);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of double-precision floating "
+                "point type (%x) --- treating as NaN.\n", item->name, item->comment, item->type);
+        return NAN;
+    }
+}
+
+psU8 psMetadataItemParseU8(const psMetadataItem *item
+                          )
+{
+    PS_ASSERT_PTR_NON_NULL(item, 0);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, U8);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, U8);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  U8);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, U8);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, U8);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  U8);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, U8);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, U8);
+        PS_METADATA_ITEM_PARSE_STRING_INT(U8,strtoul);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of integer type (%x) --- "
+                "treating as zero.\n", item->name, item->comment, item->type);
+        return 0;
+    }
+}
+
+psU16 psMetadataItemParseU16(const psMetadataItem *item
+                            )
+{
+    PS_ASSERT_PTR_NON_NULL(item, 0);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, U16);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, U16);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  U16);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, U16);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, U16);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  U16);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, U16);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, U16);
+        PS_METADATA_ITEM_PARSE_STRING_INT(U16,strtoul);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of integer type (%x) --- "
+                "treating as zero.\n", item->name, item->comment, item->type);
+        return 0;
+    }
+}
+
+psU32 psMetadataItemParseU32(const psMetadataItem *item
+                            )
+{
+    PS_ASSERT_PTR_NON_NULL(item, 0);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, U32);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, U32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  U32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, U32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, U32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  U32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, U32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, U32);
+        PS_METADATA_ITEM_PARSE_STRING_INT(U32,strtoul);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of integer type (%x) --- "
+                "treating as zero.\n", item->name, item->comment, item->type);
+        return 0;
+    }
+}
+
+psS8 psMetadataItemParseS8(const psMetadataItem *item
+                          )
+{
+    PS_ASSERT_PTR_NON_NULL(item, 0);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, S8);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, S8);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  S8);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, S8);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, S8);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  S8);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, S8);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, S8);
+        PS_METADATA_ITEM_PARSE_STRING_INT(S8,strtol);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of integer type (%x) --- "
+                "treating as zero.\n", item->name, item->comment, item->type);
+        return 0;
+    }
+}
+
+psS16 psMetadataItemParseS16(const psMetadataItem *item
+                            )
+{
+    PS_ASSERT_PTR_NON_NULL(item, 0);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, S16);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, S16);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  S16);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, S16);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, S16);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  S16);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, S16);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, S16);
+        PS_METADATA_ITEM_PARSE_STRING_INT(S16,strtol);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of integer type (%x) --- "
+                "treating as zero.\n", item->name, item->comment, item->type);
+        return 0;
+    }
+}
+
+psS32 psMetadataItemParseS32(const psMetadataItem *item
+                            )
+{
+    PS_ASSERT_PTR_NON_NULL(item, 0);
+
+    switch (item->type) {
+        PS_METADATA_ITEM_PARSE_NUMBER(F32, S32);
+        PS_METADATA_ITEM_PARSE_NUMBER(F64, S32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S8,  S32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S16, S32);
+        PS_METADATA_ITEM_PARSE_NUMBER(S32, S32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U8,  S32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U16, S32);
+        PS_METADATA_ITEM_PARSE_NUMBER(U32, S32);
+        PS_METADATA_ITEM_PARSE_STRING_INT(S32,strtol);
+        // Flow through
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of integer type (%x) --- "
+                "treating as zero.\n", item->name, item->comment, item->type);
+        return 0;
+    }
+}
+
+# define PS_METADATA_ITEM_PARSE_STRING(TYPE,MODE) \
+case PS_DATA_##TYPE: { \
+    psString value = NULL; \
+    psStringAppend(&value, MODE, item->data.TYPE); \
+    return value; } \
+
+psString psMetadataItemParseString(const psMetadataItem *item
+                                  )
+{
+    PS_ASSERT_PTR_NON_NULL(item, psStringCopy(""));
+
+    switch (item->type) {
+    case PS_DATA_STRING:
+        return psMemIncrRefCounter(item->data.V);
+
+        PS_METADATA_ITEM_PARSE_STRING(F32, "%f");
+        PS_METADATA_ITEM_PARSE_STRING(F64, "%f");
+        PS_METADATA_ITEM_PARSE_STRING(S8,  "%d");
+        PS_METADATA_ITEM_PARSE_STRING(S16, "%d");
+        PS_METADATA_ITEM_PARSE_STRING(S32, "%d");
+        PS_METADATA_ITEM_PARSE_STRING(U8,  "%d");
+        PS_METADATA_ITEM_PARSE_STRING(U16, "%d");
+        PS_METADATA_ITEM_PARSE_STRING(U32, "%d");
+
+    default:
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Item %s (%s) is not of string type (%x) --- treating as "
+                "undefined.\n", item->name, item->comment, item->type);
+        return psStringCopy("");
+    }
+}
+
Index: /trunk/psLib/src/types/psMetadataItemParse.h
===================================================================
--- /trunk/psLib/src/types/psMetadataItemParse.h	(revision 7380)
+++ /trunk/psLib/src/types/psMetadataItemParse.h	(revision 7380)
@@ -0,0 +1,19 @@
+#ifndef PS_METADATA_ITEM_PARSE_H
+#define PS_METADATA_ITEM_PARSE_H
+
+#include "psType.h"
+#include "psMetadata.h"
+
+// Parse a psMetadataItem as a particular type
+psBool psMetadataItemParseBool(const psMetadataItem *item);
+psF32 psMetadataItemParseF32(const psMetadataItem *item);
+psF64 psMetadataItemParseF64(const psMetadataItem *item);
+psS8  psMetadataItemParseS8(const psMetadataItem *item);
+psS16 psMetadataItemParseS16(const psMetadataItem *item);
+psS32 psMetadataItemParseS32(const psMetadataItem *item);
+psU8  psMetadataItemParseU8(const psMetadataItem *item);
+psU16 psMetadataItemParseU16(const psMetadataItem *item);
+psU32 psMetadataItemParseU32(const psMetadataItem *item);
+psString psMetadataItemParseString(const psMetadataItem *item);
+
+#endif
Index: /trunk/psLib/test/jpeg/Makefile.am
===================================================================
--- /trunk/psLib/test/jpeg/Makefile.am	(revision 7380)
+++ /trunk/psLib/test/jpeg/Makefile.am	(revision 7380)
@@ -0,0 +1,20 @@
+#Makefile for image jpeg tests of psLib
+#
+AM_CPPFLAGS = $(SRCINC) $(PSLIB_CFLAGS) -DVERIFIED_DIR=\"$(srcdir)/verified\"
+AM_LDFLAGS = -L$(top_builddir)/src -lpslib $(PSLIB_LIBS)
+
+TESTS = 
+
+check_DATA =
+
+check_PROGRAMS = $(TESTS)
+
+TESTS_ENVIRONMENT = perl $(top_srcdir)/test/runTest -verified=$(srcdir)/verified
+
+EXTRA_DIST = verified
+
+CLEANFILES = $(TESTS) $(check_DATA) core core.* *~
+
+tests: $(TESTS) $(check_DATA)
+
+test: check
Index: /trunk/psLib/test/math/tst_psSparse.c
===================================================================
--- /trunk/psLib/test/math/tst_psSparse.c	(revision 7380)
+++ /trunk/psLib/test/math/tst_psSparse.c	(revision 7380)
@@ -0,0 +1,59 @@
+
+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;
+    x->n = 3;
+    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;
+}
+