Index: /trunk/ppSub/.cvsignore
===================================================================
--- /trunk/ppSub/.cvsignore	(revision 13341)
+++ /trunk/ppSub/.cvsignore	(revision 13341)
@@ -0,0 +1,17 @@
+Makefile
+Makefile.in
+aclocal.m4
+autom4te.cache
+compile
+config.log
+config.status
+configure
+depcomp
+install-sh
+missing
+config.guess
+libtool
+ltmain.sh
+stamp-h1
+config.sub
+test
Index: /trunk/ppSub/Makefile.am
===================================================================
--- /trunk/ppSub/Makefile.am	(revision 13341)
+++ /trunk/ppSub/Makefile.am	(revision 13341)
@@ -0,0 +1,3 @@
+SUBDIRS = src
+
+CLEANFILES = *.pyc *~ core core.*
Index: /trunk/ppSub/autogen.sh
===================================================================
--- /trunk/ppSub/autogen.sh	(revision 13341)
+++ /trunk/ppSub/autogen.sh	(revision 13341)
@@ -0,0 +1,103 @@
+#!/bin/sh
+# Run this to generate all the initial makefiles, etc.
+
+srcdir=`dirname $0`
+test -z "$srcdir" && srcdir=.
+
+ORIGDIR=`pwd`
+cd $srcdir
+
+PROJECT=ppSub
+TEST_TYPE=-f
+# change this to be a unique filename in the top level dir
+FILE=autogen.sh
+
+DIE=0
+
+LIBTOOLIZE=libtoolize
+ACLOCAL="aclocal $ACLOCAL_FLAGS"
+AUTOHEADER=autoheader
+AUTOMAKE=automake
+AUTOCONF=autoconf
+
+($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $LIBTOOlIZE installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/"
+        DIE=1
+}
+
+($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $ACLOCAL installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOHEADER installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+        DIE=1
+}
+
+($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOMAKE installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOCONF installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+        DIE=1
+}
+
+if test "$DIE" -eq 1; then
+        exit 1
+fi
+
+test $TEST_TYPE $FILE || {
+        echo "You must run this script in the top-level $PROJECT directory"
+        exit 1
+}
+
+if test -z "$*"; then
+        echo "I am going to run ./configure with no arguments - if you wish "
+        echo "to pass any to it, please specify them on the $0 command line."
+fi
+
+$LIBTOOLIZE --copy --force || echo "$LIBTOOlIZE failed"
+$ACLOCAL || echo "$ACLOCAL failed"
+$AUTOHEADER || echo "$AUTOHEADER failed"
+$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE  failed"
+$AUTOCONF || echo "$AUTOCONF failed"
+
+cd $ORIGDIR
+
+run_configure=true
+for arg in $*; do
+    case $arg in
+        --no-configure)
+            run_configure=false
+            ;;
+        *)
+            ;;
+    esac
+done
+
+if $run_configure; then
+    $srcdir/configure --enable-maintainer-mode "$@"
+    echo
+    echo "Now type 'make' to compile $PROJECT."
+else
+    echo
+    echo "Now run 'configure' and 'make' to compile $PROJECT."
+fi
Index: /trunk/ppSub/configure.ac
===================================================================
--- /trunk/ppSub/configure.ac	(revision 13341)
+++ /trunk/ppSub/configure.ac	(revision 13341)
@@ -0,0 +1,38 @@
+dnl Process this file with autoconf to produce a configure script.
+AC_PREREQ(2.59)
+
+AC_INIT([ppSub], [0.1.1], [ipp-support@ifa.hawaii.edu])
+AC_CONFIG_SRCDIR([src])
+
+AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2])
+AM_CONFIG_HEADER([src/config.h])
+AM_MAINTAINER_MODE
+
+IPP_STDCFLAGS
+
+AC_LANG(C)
+AC_GNU_SOURCE
+AC_PROG_CC
+AC_PROG_INSTALL
+AC_PROG_LIBTOOL
+AC_SYS_LARGEFILE
+
+PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0])
+PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0])
+PKG_CHECK_MODULES([PPSTATS], [ppStats >= 1.0.0]) 
+
+dnl Set CFLAGS for build
+IPP_STDOPTS
+CFLAGS="${CFLAGS} -Wall -Werror -std=c99"
+echo "PPSUB_CFLAGS: $PPSUB_CFLAGS"
+echo "PPSUB_LIBS: $PPSUB_LIBS"
+
+AC_SUBST([PPSUB_CFLAGS])
+AC_SUBST([PPSUB_LIBS])
+
+AC_CONFIG_FILES([
+  Makefile
+  src/Makefile
+])
+
+AC_OUTPUT
Index: /trunk/ppSub/src/.cvsignore
===================================================================
--- /trunk/ppSub/src/.cvsignore	(revision 13341)
+++ /trunk/ppSub/src/.cvsignore	(revision 13341)
@@ -0,0 +1,10 @@
+*.o
+*.lo
+.libs
+.deps
+Makefile
+Makefile.in
+ppSub
+config.h
+config.h.in
+stamp-h1
Index: /trunk/ppSub/src/Makefile.am
===================================================================
--- /trunk/ppSub/src/Makefile.am	(revision 13341)
+++ /trunk/ppSub/src/Makefile.am	(revision 13341)
@@ -0,0 +1,22 @@
+bin_PROGRAMS = ppSub
+ppSub_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PPSUB_CFLAGS)
+ppSub_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PPSUB_LIBS)
+
+ppSub_SOURCES =			\
+	ppSub.c			\
+	ppSubArguments.c	\
+	ppSubCamera.c		\
+	ppSubLoop.c		\
+	ppSubReadout.c		\
+	ppSubVersion.c            
+
+noinst_HEADERS = \
+	ppSub.h
+
+clean-local:
+	-rm -f TAGS
+
+# Tags for emacs
+tags:
+	etags `find . -name \*.[ch] -print`
+
Index: /trunk/ppSub/src/ppSub.c
===================================================================
--- /trunk/ppSub/src/ppSub.c	(revision 13341)
+++ /trunk/ppSub/src/ppSub.c	(revision 13341)
@@ -0,0 +1,51 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppSub.h"
+
+int main(int argc, char *argv[])
+{
+    psExit exitValue = PS_EXIT_SUCCESS; // Exit value
+    psTimerStart("ppSub");
+    psLibInit(NULL);
+
+    pmConfig *config = pmConfigRead(&argc, argv, PPSUB_RECIPE); // Configuration
+    if (!config) {
+        psErrorStackPrint(stderr, "Error reading configuration.");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    if (!ppSubArguments(argc, argv, config)) {
+        psErrorStackPrint(stderr, "Error reading arguments.\n");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    if (!ppSubCamera(config)) {
+        psErrorStackPrint(stderr, "Error reading configuration.\n");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    if (!ppSubLoop(config)) {
+        psErrorStackPrint(stderr, "Error performing subtraction.\n");
+        exitValue = PS_EXIT_PROG_ERROR;
+        goto die;
+    }
+
+ die:
+    psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub"));
+    psTimerStop();
+
+    psFree(config);
+    pmConfigDone();
+    psLibFinalize();
+
+    exit(exitValue);
+}
Index: /trunk/ppSub/src/ppSub.h
===================================================================
--- /trunk/ppSub/src/ppSub.h	(revision 13341)
+++ /trunk/ppSub/src/ppSub.h	(revision 13341)
@@ -0,0 +1,29 @@
+#ifndef PP_SUB_H
+#define PP_SUB_H
+
+#define PPSUB_RECIPE "PPSUB"            /// Name of the recipe to use
+
+/// Parse the arguments
+bool ppSubArguments(int argc, char *argv[], ///< Command-line arguments
+                    pmConfig *config    ///< Configuration
+    );
+
+/// Parse the camera input
+bool ppSubCamera(pmConfig *config       ///< Configuration
+    );
+
+/// Loop over the FPA hierarchy
+bool ppSubLoop(pmConfig *config         ///< Configuration
+    );
+
+/// Perform PSF-matched image subtraction on the readout
+bool ppSubReadout(pmConfig *config,     ///< Configuration
+                  const pmFPAview *view ///< View of readout to subtract
+    );
+
+/// Put the program version information into a metadata
+void ppSubVersionMetadata(psMetadata *metadata ///< Metadata to populate
+    );
+
+
+#endif
Index: /trunk/ppSub/src/ppSubArguments.c
===================================================================
--- /trunk/ppSub/src/ppSubArguments.c	(revision 13341)
+++ /trunk/ppSub/src/ppSubArguments.c	(revision 13341)
@@ -0,0 +1,231 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppSub.h"
+
+// Print usage information and die
+static void usage(const char *program,  // Name of the program
+                  psMetadata *arguments, // Command-line arguments
+                  pmConfig *config      // Configuration
+    )
+{
+    fprintf(stderr, "\nPan-STARRS PSF-matched image subtraction\n\n");
+    fprintf(stderr, "Usage: %s INPUT.fits REFERENCE.fits OUTPUT.fits\n", program);
+    fprintf(stderr, "\n");
+    psArgumentHelp(arguments);
+    psFree(arguments);
+    psFree(config);
+    pmConfigDone();
+    psLibFinalize();
+    exit(PS_EXIT_CONFIG_ERROR);
+}
+
+// Get a float-point value from the command-line or recipe, and add it to the arguments
+#define VALUE_ARG_RECIPE_FLOAT(ARGNAME, RECIPENAME, TYPE) { \
+    ps##TYPE value = psMetadataLookup##TYPE(NULL, arguments, ARGNAME); \
+    if (isnan(value)) { \
+        bool mdok; \
+        value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \
+        if (!mdok) { \
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
+                RECIPENAME, PPSUB_RECIPE); \
+            goto ERROR; \
+        } \
+    } \
+    psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \
+}
+
+// Get an integer value from the command-line or recipe, and add it to the arguments
+#define VALUE_ARG_RECIPE_INT(ARGNAME, RECIPENAME, TYPE, UNSET) { \
+    ps##TYPE value = psMetadataLookup##TYPE(NULL, arguments, ARGNAME); \
+    if (value == UNSET) { \
+        bool mdok; \
+        value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \
+        if (!mdok) { \
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
+                RECIPENAME, PPSUB_RECIPE); \
+            goto ERROR; \
+        } \
+    } \
+    psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \
+}
+
+// Get a string value from the command-line and add it to the target
+static bool valueArgStr(pmConfig *config,      // Configuration
+                        psMetadata *arguments, // Command-line arguments
+                        const char *argName, // Argument name in the command-line arguments
+                        const char *mdName, // Name for value in the metadata
+                        psMetadata *target // Target metadata to which to add value
+                        )
+{
+    psString value = psMetadataLookupStr(NULL, arguments, argName); // Value of interest
+    if (value && strlen(value) > 0) {
+        return psMetadataAddStr(target, PS_LIST_TAIL, mdName, 0, NULL, value);
+    }
+    return false;
+}
+
+// Get a vector from the command-line or recipe, and add it to the target
+static bool vectorArgRecipe(pmConfig *config,   // Configuration
+                            psMetadata *arguments, // Command-line arguments
+                            const char *argName, // Argument name in the command-line arguments
+                            const psMetadata *recipe, // Recipe
+                            const char *recipeName, // Name for value in the recipe
+                            psMetadata *target, // Target to which to add value
+                            psElemType type // Type for vector
+    )
+{
+    psVector *vector;                   // Vector
+    psString string = psMetadataLookupStr(NULL, arguments, argName); // String from arguments
+    if (string) {
+        psArray *array = psStringSplitArray(string, ", ", false); // Array of strings
+        vector = psVectorAlloc(array->n, PS_TYPE_F32);
+        for (int i = 0; i < array->n; i++) {
+            const char *subString = array->data[i]; // String with a value
+            char *end;                  // Ptr to end of string parsed
+
+            switch (type) {
+              case PS_TYPE_F32:
+                vector->data.F32[i] = strtof(subString, &end);
+                break;
+              case PS_TYPE_S32: {
+                  long value = strtol(subString, &end, 10);
+                  if (value > PS_MAX_S32 || value < PS_MIN_S32) {
+                      psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                              "%s list includes value beyond S32 representation: %s",
+                              argName, string);
+                      psFree(vector);
+                      return false;
+                  }
+                  vector->data.S32[i] = value;
+                  break;
+              }
+              default:
+                psAbort("Unsupported type: %x\n", type);
+            }
+            if (end == subString) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to decipher %s list: %s",
+                        argName, string);
+                psFree(vector);
+                return false;
+            }
+        }
+        psFree(array);
+    } else {
+        vector = psMetadataLookupPtr(NULL, recipe, recipeName);
+        if (!psMemCheckVector(vector) || vector->type.type != type) {
+            psError(PS_ERR_BAD_PARAMETER_TYPE, false, "%s in recipe %s is not a vector of type F32.",
+                    recipeName, PPSUB_RECIPE);
+            return false;
+        }
+    }
+
+    psMetadataAddVector(target, PS_LIST_TAIL, recipeName, 0, NULL, vector);
+    psFree(vector);                     // Drop reference
+
+    return true;
+}
+
+
+bool ppSubArguments(int argc, char *argv[], pmConfig *config)
+{
+    assert(config);
+
+    psMetadata *arguments = psMetadataAlloc(); // Command-line arguments
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-inmask", 0, "Input mask image", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-inweight", 0, "Input weight image", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-refmask", 0, "Referencemask image", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-refweight", 0, "Referenceweight image", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stat", 0, "Statistics file", NULL);
+
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-size", 0, "Kernel half-size (pixels)", 0);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-order", 0, "Spatial polynomial order", 0);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Kernel type (POIS or ISIS)", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-widths", 0, "ISIS Gaussian widths (comma-separated)", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-orders", 0, "ISIS polynomial orders (comma-separated)", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-spacing", 0, "Typical stamp spacing (pixels)", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-footprint", 0, "Stamp footprint half-size (pixels)", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold", 0, "Minimum threshold for stamps (ADU)", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-rej", 0, "Rejection thresold (sigma)", NAN);
+    psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-bad", 0, "Mask value for bad pixels", 0);
+    psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-blank", 0, "Mask value for blank region", 0);
+
+    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 4) {
+        usage(argv[0], arguments, config);
+    }
+
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT", 0, "Name of the input image", argv[1]);
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "REF", 0, "Name of the reference image", argv[2]);
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[3]);
+
+    valueArgStr(config, arguments, "-inmask",    "INPUT.MASK",    config->arguments);
+    valueArgStr(config, arguments, "-inweight",  "INPUT.WEIGHT",  config->arguments);
+    valueArgStr(config, arguments, "-refmask",   "REF.MASK",      config->arguments);
+    valueArgStr(config, arguments, "-refweight", "REF.WEIGHT",    config->arguments);
+    valueArgStr(config, arguments, "-stat",      "STATS",         config->arguments);
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
+    if (!recipe) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE);
+        goto ERROR;
+    }
+
+    VALUE_ARG_RECIPE_INT("-size",        "KERNEL.SIZE",     S32, 0);
+    VALUE_ARG_RECIPE_INT("-order",       "SPATIAL.ORDER",   S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-spacing",   "STAMP.SPACING",   F32);
+    VALUE_ARG_RECIPE_INT("-footprint",   "STAMP.FOOTPRINT", S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-threshold", "STAMP.THRESHOLD", F32);
+    VALUE_ARG_RECIPE_INT("-iter",        "ITER",            S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-rej",       "REJ",             F32);
+    VALUE_ARG_RECIPE_INT("-mask-bad",    "MASK.BAD",        U8, 0);
+    VALUE_ARG_RECIPE_INT("-mask-blank",  "MASK.BLANK",      U8, 0);
+
+    vectorArgRecipe(config, arguments, "-isis-widths", recipe, "ISIS.WIDTHS", config->arguments, PS_TYPE_F32);
+    vectorArgRecipe(config, arguments, "-isis-orders", recipe, "ISIS.ORDERS", config->arguments, PS_TYPE_S32);
+
+    psVector *widths = psMetadataLookupPtr(NULL, config->arguments, "ISIS.WIDTHS"); // ISIS Gaussian widths
+    psVector *orders = psMetadataLookupPtr(NULL, config->arguments, "ISIS.ORDERS"); // ISIS Polynomial orders
+    if (widths->n != orders->n) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Size of vectors for ISIS widths and orders do not match.");
+        goto ERROR;
+    }
+
+    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "REVERSE", 0, "Reverse sense of subtraction",
+                      psMetadataLookupBool(NULL, arguments, "-reverse"));
+
+    // Translate the kernel type
+    psString type = psMetadataLookupStr(NULL, arguments, "-type"); // Name of kernel type
+    if (!type || strlen(type) == 0) {
+        type = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE");
+        if (!type || strlen(type) == 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find KERNEL.TYPE specified.");
+            goto ERROR;
+        }
+    }
+    pmSubtractionKernelsType kernelType; // Type of kernel
+    if (strcasecmp(type, "POIS") == 0) {
+        kernelType = PM_SUBTRACTION_KERNEL_POIS;
+    } else if (strcasecmp(type, "ISIS") == 0) {
+        kernelType = PM_SUBTRACTION_KERNEL_ISIS;
+    } else {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", type);
+        goto ERROR;
+    }
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "KERNEL.TYPE", 0, "Type of kernel", kernelType);
+
+    psTrace("ppSub", 1, "Done reading command-line arguments\n");
+    return true;
+
+ERROR:
+    psFree(arguments);
+    return false;
+}
+
+
Index: /trunk/ppSub/src/ppSubCamera.c
===================================================================
--- /trunk/ppSub/src/ppSubCamera.c	(revision 13341)
+++ /trunk/ppSub/src/ppSubCamera.c	(revision 13341)
@@ -0,0 +1,100 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppSub.h"
+
+bool ppSubCamera(pmConfig *config)
+{
+    bool status = false;                // Status of definition
+
+    // Input image
+    pmFPAfile *input = pmFPAfileDefineFromArgs(&status, config, "PPSUB.INPUT", "INPUT");
+    if (!status || !input) {
+        psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.INPUT");
+        return false;
+    }
+    if (input->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_IO, true, "PPSUB.INPUT is not of type IMAGE");
+        return false;
+    }
+
+    // Input mask
+    pmFPAfile *inputMask = pmFPAfileBindFromArgs(NULL, input, config, "PPSUB.INPUT.MASK", "INPUT.MASK");
+    if (inputMask && inputMask->type != PM_FPA_FILE_MASK) {
+        psError(PS_ERR_IO, true, "PPSUB.INPUT.MASK is not of type MASK");
+        return false;
+    }
+
+    // Input weight map
+    pmFPAfile *inputWeight = pmFPAfileBindFromArgs(NULL, input, config, "PPSUB.INPUT.WEIGHT", "INPUT.WEIGHT");
+    if (inputWeight && inputWeight->type != PM_FPA_FILE_WEIGHT) {
+        psError(PS_ERR_IO, true, "PPSUB.INPUT.WEIGHT is not of type WEIGHT");
+        return false;
+    }
+
+    // Reference image
+    status = false;
+    pmFPAfile *ref = pmFPAfileDefineFromArgs(&status, config, "PPSUB.REF", "REF");
+    if (!ref) {
+        psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.REF");
+        return false;
+    }
+    if (ref->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_IO, true, "PPSUB.REF is not of type IMAGE");
+        return false;
+    }
+
+    // Reference mask
+    pmFPAfile *refMask = pmFPAfileBindFromArgs(NULL, input, config, "PPSUB.REF.MASK", "REF.MASK");
+    if (refMask && refMask->type != PM_FPA_FILE_MASK) {
+        psError(PS_ERR_IO, true, "PPSUB.REF.MASK is not of type MASK");
+        return false;
+    }
+
+    // Reference weight map
+    pmFPAfile *refWeight = pmFPAfileBindFromArgs(NULL, input, config, "PPSUB.REF.WEIGHT", "REF.WEIGHT");
+    if (refWeight && refWeight->type != PM_FPA_FILE_WEIGHT) {
+        psError(PS_ERR_IO, true, "PPSUB.REF.WEIGHT is not of type WEIGHT");
+        return false;
+    }
+
+    // Output image
+    pmFPAfile *output = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT");
+    if (!output) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT"));
+        return false;
+    }
+    if (output->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_IO, true, "PPSUB.OUTPUT is not of type IMAGE");
+        return false;
+    }
+
+    // Output mask
+    pmFPAfile *outMask = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.OUTPUT.MASK");
+    if (!outMask) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.MASK"));
+        return false;
+    }
+    if (outMask->type != PM_FPA_FILE_MASK) {
+        psError(PS_ERR_IO, true, "PPSUB.OUTPUT.MASK is not of type MASK");
+        return false;
+    }
+
+    // Output weight
+    pmFPAfile *outWeight = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.OUTPUT.WEIGHT");
+    if (!outWeight) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.WEIGHT"));
+        return false;
+    }
+    if (outWeight->type != PM_FPA_FILE_WEIGHT) {
+        psError(PS_ERR_IO, true, "PPSUB.OUTPUT.WEIGHT is not of type WEIGHT");
+        return false;
+    }
+
+    return true;
+}
Index: /trunk/ppSub/src/ppSubLoop.c
===================================================================
--- /trunk/ppSub/src/ppSubLoop.c	(revision 13341)
+++ /trunk/ppSub/src/ppSubLoop.c	(revision 13341)
@@ -0,0 +1,171 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <ppStats.h>
+
+#include "ppSub.h"
+
+bool ppSubLoop(pmConfig *config)
+{
+    const char *statsName = psMetadataLookupStr(NULL, config->arguments, "STATS"); // Filename for statistics
+    psMetadata *stats = NULL;           // Container for statistics
+    FILE *statsFile = NULL;             // File stream for statistics
+    if (statsName && strlen(statsName) > 0) {
+        psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename
+        statsFile = fopen(resolved, "w");
+        if (!statsFile) {
+            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
+            psFree(resolved);
+            return false;
+        } else {
+            stats = psMetadataAlloc();
+        }
+        psFree(resolved);
+    }
+
+    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT");
+    if (!input) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n");
+        return false;
+    }
+
+    pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF");
+    if (!reference) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find reference data!\n");
+        return false;
+    }
+
+    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT");
+    if (!output) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n");
+        return false;
+    }
+
+    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
+    pmHDU *lastHDU = NULL;              // Last HDU that was updated
+
+    // Iterate over the FPA hierarchy
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return false;
+    }
+
+    pmChip *inChip;                    // Input chip of interest
+    while ((inChip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
+        pmChip *refChip = pmFPAviewThisChip(view, reference->fpa); // Reference chip of interest
+        if ((!inChip->file_exists && refChip->file_exists) ||
+            (inChip->file_exists && !refChip->file_exists)) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "FPA format discrepency between input and reference");
+            psFree(view);
+            return false;
+        }
+
+        if (!inChip->file_exists) {
+            continue;
+        }
+
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            return false;
+        }
+
+        pmCell *inCell;                // Cell of interest
+        while ((inCell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
+            pmCell *refCell = pmFPAviewThisCell(view, reference->fpa); // Reference cell of interest
+            if ((!inCell->file_exists && refCell->file_exists) ||
+                (inCell->file_exists && !refCell->file_exists)) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "FPA format discrepency between input and reference");
+                psFree(view);
+                return false;
+            }
+            if (!inCell->file_exists) {
+                continue;
+            }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                return false;
+            }
+
+            // Put version information into the header
+            pmHDU *hdu = pmHDUFromCell(inCell);
+            if (hdu && hdu != lastHDU) {
+                if (!hdu->header) {
+                    hdu->header = psMetadataAlloc();
+                }
+                ppSubVersionMetadata(hdu->header);
+                lastHDU = hdu;
+            }
+
+            pmReadout *inRO;           // Readin of interest
+            while ((inRO = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    return false;
+                }
+                pmReadout *refRO = pmFPAviewThisReadout(view, reference->fpa);// Reference readout of interest
+                if ((!inRO->data_exists && refRO->data_exists) ||
+                    (inRO->data_exists && !refRO->data_exists)) {
+                    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                            "FPA format discrepency between input and reference");
+                    psFree(view);
+                    return false;
+                }
+                if (!inRO->data_exists) {
+                    continue;
+                }
+
+                // Perform the analysis
+                if (!ppSubReadout(config, view)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to detrend inRO.\n");
+                    return false;
+                }
+
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    return false;
+                }
+            }
+
+            // Perform statistics on the cell
+            if (stats) {
+                pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); // Output file
+                if (!output) {
+                    psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPSUB.OUTPUT.\n");
+                    return false;
+                }
+                ppStats(stats, output->fpa, view, config);
+            }
+
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                return false;
+            }
+        }
+
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            return false;
+        }
+    }
+
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+
+    psFree(view);
+
+    // Write out summary statistics
+    if (stats) {
+        const char *statsMDC = psMetadataConfigFormat(stats);
+        if (!statsMDC || strlen(statsMDC) == 0) {
+            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
+        } else {
+            fprintf(statsFile, "%s", statsMDC);
+        }
+        psFree((void *)statsMDC);
+        fclose(statsFile);
+
+        psFree(stats);
+    }
+
+    return true;
+}
Index: /trunk/ppSub/src/ppSubReadout.c
===================================================================
--- /trunk/ppSub/src/ppSubReadout.c	(revision 13341)
+++ /trunk/ppSub/src/ppSubReadout.c	(revision 13341)
@@ -0,0 +1,174 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppSub.h"
+
+#define MASK_BAD      0x01              // Mask value for bad pixel
+#define MASK_STAMP    0x02              // Mask value for bad stamp (and bad stamp region)
+
+bool ppSubReadout(pmConfig *config, const pmFPAview *view)
+{
+    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
+    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REFERENCE"); // Reference readout
+    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
+    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+
+    psImage *input = inRO->image;       // Input image
+    psImage *reference = refRO->image;  // Reference image
+    PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
+
+    // Look up appropriate values
+    int size = psMetadataLookupS32(NULL, config->arguments, "KERNEL.SIZE"); // Kernel half-size
+    int order = psMetadataLookupS32(NULL, config->arguments, "SPATIAL.ORDER"); // Spatial polynomial order
+    float spacing = psMetadataLookupF32(NULL, config->arguments, "STAMP.SPACING"); // Typical stamp spacing
+    int footprint = psMetadataLookupS32(NULL, config->arguments, "STAMP.FOOTPRINT"); // Stamp half-size
+    float threshold = psMetadataLookupF32(NULL, config->arguments, "STAMP.THRESHOLD"); // Threshold for stmps
+    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
+    float rej = psMetadataLookupF32(NULL, config->arguments, "REJ"); // Rejection threshold
+    pmSubtractionKernelsType type = psMetadataLookupBool(NULL, config->arguments,
+                                                         "KERNEL.TYPE"); // Kernel type
+    bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction?
+    psVector *widths = psMetadataLookupPtr(NULL, config->arguments, "ISIS.WIDTHS"); // ISIS Gaussian widths
+    psVector *orders = psMetadataLookupPtr(NULL, config->arguments, "ISIS.ORDERS"); // ISIS Polynomial orders
+    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
+    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
+
+    if (!inRO->mask && !pmReadoutGenerateMask(inRO)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for input image");
+        return false;
+    }
+    if (!inRO->weight && !pmReadoutGenerateWeight(inRO, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for input image");
+        return false;
+    }
+    if (!refRO->mask && !pmReadoutGenerateMask(refRO)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for reference image");
+        return false;
+    }
+    if (!refRO->weight && !pmReadoutGenerateWeight(refRO, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for reference image");
+        return false;
+    }
+
+    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
+    int numCols = input->numCols, numRows = input->numRows; // Image dimensions
+    psImage *stampMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // Mask to use for stamps
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            stampMask->data.PS_TYPE_MASK_DATA[y][x] =
+                (refRO->mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) ? MASK_BAD : 0;
+        }
+    }
+
+#if 0
+    if (!inRO->weight) {
+        // Need to ensure input image is completely above zero, so we're not weighting by negative numbers
+        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+        if (!psImageBackground(stats, inRO->image, inRO->mask, maskBad, NULL)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine statistics of input image.");
+            return false;
+        }
+        float diff = stats->robustMedian - 5.0 * stats->robustStdev;
+        if (diff < 0.0) {
+            psLogMsg("ppSub", PS_LOG_INFO, "Adding %f to image background to get positive pixels.");
+            (void)psBinaryOp(in->image, inRO->image, "+", psScalarAlloc(diff, PS_TYPE_F32));
+        }
+    }
+#endif
+
+    pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order,
+                                                                 widths, orders); // Kernel basis functions
+    psArray *stamps = NULL;             // Stamps for matching PSF
+    psVector *solution = NULL;          // Solution to match PSF
+
+    for (int i = 0, numRejected = 0; i < iter && numRejected != 0; i++) {
+        stamps = pmSubtractionFindStamps(stamps, refRO->image, stampMask, MASK_BAD, MASK_STAMP,
+                                         threshold, spacing, footprint);
+        if (!stamps) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image.");
+            goto ERROR;
+        }
+
+        if (!pmSubtractionCalculateEquation(stamps, refRO->image, inRO->image, inRO->weight,
+                                            kernels, footprint)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+            goto ERROR;
+         }
+
+        solution = pmSubtractionSolveEquation(solution, stamps);
+        if (!solution) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+            goto ERROR;
+        }
+
+        numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, stampMask, MASK_STAMP,
+                                                solution, footprint, rej, kernels);
+        if (numRejected < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+            goto ERROR;
+        }
+        psLogMsg("ppSub", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, i);
+    }
+
+    if (!pmSubtractionConvolve(&outRO->image, &outRO->weight, &outRO->mask,
+                               refRO->image, refRO->weight, refRO->mask,
+                               MASK_BAD, maskBlank, solution, kernels)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
+        goto ERROR;
+    }
+
+    // Do the subtraction
+    if (reverse) {
+        (void)psBinaryOp(outRO->image, outRO->image, "-", inRO->image);
+    } else {
+        (void)psBinaryOp(outRO->image, inRO->image, "-", outRO->image);
+    }
+    (void)psBinaryOp(outRO->mask, outRO->mask, "|", inRO->mask);
+    (void)psBinaryOp(outRO->weight, outRO->weight, "+", inRO->weight);
+
+#if 0
+    // Generate image with convolution kernels
+    int fullSize = 2 * size + 1;        // Full size of kernel
+    psImage *convKernels = psImageAlloc(5 * fullSize, 5 * fullSize, PS_TYPE_F32);
+    for (int j = -2; j <= 2; j++) {
+        for (int i = -2; i <= 2; i++) {
+            psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
+                                                       (float)j / 2.0); // Image of the kernel
+            if (!kernel) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+                psFree(convKernels);
+                goto ERROR;
+            }
+            if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
+                                      (j + 2) * fullSize, "=") == 0) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
+                psFree(kernel);
+                psFree(convKernels);
+                goto ERROR;
+            }
+        }
+    }
+
+    // XXX What do we do with this image?
+    psFree(convKernels);
+#endif
+
+
+        psFree(kernels);
+        psFree(stamps);
+        psFree(solution);
+        psFree(outRO);
+        return true;
+
+    ERROR:
+        psFree(kernels);
+        psFree(stamps);
+        psFree(solution);
+        psFree(outRO);
+        return false;
+}
Index: /trunk/ppSub/src/ppSubVersion.c
===================================================================
--- /trunk/ppSub/src/ppSubVersion.c	(revision 13341)
+++ /trunk/ppSub/src/ppSubVersion.c	(revision 13341)
@@ -0,0 +1,60 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <ppStats.h>
+
+#include "ppSub.h"
+
+static const char *cvsTag = "$Name: not supported by cvs2svn $";// CVS tag name
+
+psString ppSubVersion(void)
+{
+    psString version = NULL;            // Version, to return
+    psStringAppend(&version, "%s-%s",PACKAGE_NAME,PACKAGE_VERSION);
+    return version;
+}
+
+psString ppSubVersionLong(void)
+{
+    psString version = ppSubVersion(); // Version, to return
+    psString tag = psStringStripCVS(cvsTag, "Name"); // CVS tag
+    psStringAppend(&version, " (cvs tag %s) %s, %s", tag, __DATE__, __TIME__);
+    psFree(tag);
+    return version;
+}
+
+
+void ppSubVersionMetadata(psMetadata *metadata)
+{
+    PS_ASSERT_METADATA_NON_NULL(metadata,);
+
+    psString pslib = psLibVersionLong();// psLib version
+    psString psmodules = psModulesVersionLong(); // psModules version
+    psString ppStats = ppStatsVersionLong(); // ppStats version
+    psString ppSub = ppSubVersionLong(); // ppSub version
+
+    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
+    psString timeString = psTimeToISO(time); // The time in an ISO string
+    psFree(time);
+    psString head = NULL;               // Head string
+    psStringAppend(&head, "ppSub processing at %s. Component information:", timeString);
+    psFree(timeString);
+
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, head, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, pslib, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, psmodules, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, ppStats, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, ppSub, "");
+
+    psFree(head);
+    psFree(pslib);
+    psFree(psmodules);
+    psFree(ppStats);
+    psFree(ppSub);
+
+    return;
+}