Index: /trunk/ppImage/src/ppConfigLoad.c
===================================================================
--- /trunk/ppImage/src/ppConfigLoad.c	(revision 5857)
+++ /trunk/ppImage/src/ppConfigLoad.c	(revision 5857)
@@ -0,0 +1,34 @@
+# include "ppImage.h"
+
+bool ppConfigLoad (ppConfig *config, int argc, char **argv) {
+
+    ppConfig config;
+
+    // Parse the configurations
+    config->site = NULL;            // Site configuration
+    config->camera = NULL;          // Camera configuration
+    config->recipe = NULL;          // Recipe configuration
+
+    if (! pmConfigRead(&config->site, &config->camera, &config->recipe, &argc, argv, RECIPE)) {
+        psErrorStackPrint(stderr, "Can't find site configuration!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    // Parse other command-line arguments
+    config->arguments = psMetadataAlloc(); // The arguments, with default values
+
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "-bias", 0, "Name of the bias image", "");
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "-dark", 0, "Name of the dark image", "");
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "-flat", 0, "Name of the flat-field image", "");
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "-mask", 0, "Name of the mask image", "");
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "-chip", 0, "Chip number to process (if positive)", -1);
+    if (! psArgumentParse(config->arguments, &argc, argv) || argc != 3) {
+        printf("\nPan-STARRS Phase 2 processing\n\n");
+        printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]);
+        psArgumentHelp(config->arguments);
+        psFree(config->arguments);
+        exit(EXIT_FAILURE);
+    }
+    config->database = pmConfigDB(config->site);  // Database handle
+    return true;
+} 
Index: /trunk/ppImage/src/ppDataLoadMeta.c
===================================================================
--- /trunk/ppImage/src/ppDataLoadMeta.c	(revision 5857)
+++ /trunk/ppImage/src/ppDataLoadMeta.c	(revision 5857)
@@ -0,0 +1,54 @@
+# include "ppImage.h"
+
+bool ppDataLoadMeta (ppData *data, ppConfig *config, char **argv) {
+
+    const char *inputName = argv[1];    // Name of input image
+    const char *outputName = argv[2];   // Name of output image
+
+    // Open the input image
+    psLogMsg("phase2", PS_LOG_INFO, "Opening input image: %s\n", inputName);
+    data->inputFile = psFitsOpen(inputName, "r"); // File handle for FITS file
+    if (! data->inputFile) {
+        psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName);
+        exit(EXIT_FAILURE);
+    }
+    data->header = psFitsReadHeader(NULL, data->inputFile); // FITS header
+
+    // Open the output image
+    // We do it here so that we don't process the whole lot and then find out we can't open the output file
+    data->outputFile = psFitsOpen(outputName, "w");
+    if (! data->outputFile) {
+        psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName);
+        exit(EXIT_FAILURE);
+    }
+
+    // Open the output mask
+    // XXX EAM : not done
+
+    // Get camera configuration from header if not already defined
+    if (! config->camera) {
+        config->camera = pmConfigCameraFromHeader(config->site, config->header);
+        if (! config->camera) {
+            psErrorStackPrint(stderr, "Can't find camera configuration!\n");
+            exit(EXIT_FAILURE);
+        }
+   } else if (! pmConfigValidateCamera(config->camera, config->header)) {
+        psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName);
+        exit(EXIT_FAILURE);
+    }
+
+    // determine the correct recipe to use
+    if (! config->recipe && !(config->recipe = pmConfigRecipeFromCamera(config->camera, RECIPE))) {
+        psErrorStackPrint(stderr, "Can't find recipe configuration!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    // Construct camera in preparation for reading
+    data->input = pmFPAConstruct(camera);
+    data->mask  = pmFPAConstruct(camera);
+    data->bias  = pmFPAConstruct(camera);
+    data->dark  = pmFPAConstruct(camera);
+    data->flat  = pmFPAConstruct(camera);
+
+    return true;
+}
Index: /trunk/ppImage/src/ppDataLoadPixels.c
===================================================================
--- /trunk/ppImage/src/ppDataLoadPixels.c	(revision 5857)
+++ /trunk/ppImage/src/ppDataLoadPixels.c	(revision 5857)
@@ -0,0 +1,137 @@
+# include "ppImage.h"
+
+// XXX EAM : this function reads in all data from all chips in one pass
+bool ppDataLoadPixels (ppData *data, ppRecipe *options, ppConfig *config) {
+
+    // if we are using a single chip, select it
+    if (options->chipNum >= 0) {
+        if (! pmFPASelectChip(data->input, options->chipNum)) {
+            psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", options->chipNum);
+            exit(EXIT_FAILURE);
+        }
+	if (! pmFPASelectChip(data->bias, options->chipNum)) psAbort (__func__, "invalid missing bias chip");
+	if (! pmFPASelectChip(data->dark, options->chipNum)) psAbort (__func__, "invalid missing dark chip");
+	if (! pmFPASelectChip(data->mask, options->chipNum)) psAbort (__func__, "invalid missing mask chip");
+
+        psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", options->chipNum);
+    }
+    
+    // Read in the input pixels
+    if (! pmFPARead(data->input, data->inputFile, data->header, config->database)) {
+        psErrorStackPrint(stderr, "Unable to populate camera from input FITS file\n");
+        exit(EXIT_FAILURE);
+    }
+    psFitsClose(data->inputFile);
+
+    // Generate input mask and weight frame
+    // ** note the distinction between the input->mask and the external mask **
+    // XXX is this a sensible distinction? **
+    p_pmHDU *hdu = NULL;
+    if (data->input->hdu) { hdu = data->input->hdu; }
+    psArray *chips = data->input->chips;
+    for (int chipNum = 0; chipNum < chips->n; chipNum++) {
+	pmChip *chip = chips->data[chipNum];
+	if (! chip->valid) { continue; } 
+	if (chip->hdu) { hdu = chip->hdu; }
+
+	psArray *cells = chip->cells;
+	for (int cellNum = 0; cellNum < cells->n; cellNum++) {
+	    pmCell *cell = cells->data[cellNum];
+	    if (! cell->valid) { continue; }
+	    if (cell->hdu) { hdu = cell->hdu; }
+
+	    hdu->masks = psArrayAlloc(hdu->images->n);
+	    hdu->weights = psArrayAlloc(hdu->images->n);
+	    psArray *readouts = cell->readouts;
+	    for (int readNum = 0; readNum < hdu->images->n; readNum++) {
+		psImage *image = hdu->images->data[readNum];
+		psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
+		psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+		pmReadout *readout = readouts->data[readNum];
+
+		// XXX EAM : this needs a gain and read-noise for the weight
+		// XXX EAM : should weight be stdev or var ?
+		for (int j = 0; j < image->numRows; j++) {
+		    for (int i = 0; i < image->numCols; i++) {
+			mask->data.U8[j][i] = 0;
+			weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]);
+		    }
+		}
+		hdu->masks->data[readNum] = mask;
+		hdu->weights->data[readNum] = weight;
+		readout->mask = psMemIncrRefCounter(mask);
+		readout->weight = psMemIncrRefCounter(weight);
+	    }
+	}
+    }
+
+    // Load the calibration frames, if required
+    // XXX : by using the entries like data->biasFile we can read in the inner loop
+    if (options->doBias) {
+        psFits *biasFile = psFitsOpen(options->biasName, "r"); // File handle for bias
+        psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias
+        if (! pmConfigValidateCamera(config->camera, biasHeader)) {
+            psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", options->biasName);
+            exit(EXIT_FAILURE);
+        }
+        psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", options->biasName);
+        if (! pmFPARead(data->bias, biasFile, biasHeader, config->database)) {
+            psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", options->biasName);
+            exit(EXIT_FAILURE);
+        }
+        psFree(biasHeader);
+        psFitsClose(biasFile);
+    }
+
+    if (options->doDark) {
+        psFits *darkFile = psFitsOpen(options->darkName, "r"); // File handle for dark
+        psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark
+        if (! pmConfigValidateCamera(config->camera, darkHeader)) {
+            psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", options->darkName);
+            exit(EXIT_FAILURE);
+        }
+        psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", options->darkName);
+        if (! pmFPARead(data->dark, darkFile, darkHeader, config->database)) {
+            psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", options->darkName);
+            exit(EXIT_FAILURE);
+        }
+        psFree(darkHeader);
+        psFitsClose(darkFile);
+    }
+
+    if (options->doMask) {
+        psFits *maskFile = psFitsOpen(options->maskName, "r"); // File handle for mask
+        psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask
+        if (! pmConfigValidateCamera(config->camera, maskHeader)) {
+            psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", options->maskName);
+            exit(EXIT_FAILURE);
+        }
+        psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", options->maskName);
+        if (! pmFPARead(data->mask, maskFile, maskHeader, config->database)) {
+            psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", options->maskName);
+            exit(EXIT_FAILURE);
+        }
+        psFree(maskHeader);
+        psFitsClose(maskFile);
+    }
+
+    if (options->doFlat) {
+        psFits *flatFile = psFitsOpen(options->flatName, "r"); // File handle for flat
+        if (! flatFile) {
+            psErrorStackPrint(stderr, "Can't open flat image: %s\n", options->flatName);
+            exit(EXIT_FAILURE);
+        }
+        psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat
+        if (! pmConfigValidateCamera(config->camera, flatHeader)) {
+            psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", options->flatName);
+            exit(EXIT_FAILURE);
+        }
+        psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", options->flatName);
+        if (! pmFPARead(data->flat, flatFile, flatHeader, config->database)) {
+            psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", options->flatName);
+            exit(EXIT_FAILURE);
+        }
+        psFree(flatHeader);
+        psFitsClose(flatFile);
+    }
+
Index: /trunk/ppImage/src/ppImage.c
===================================================================
--- /trunk/ppImage/src/ppImage.c	(revision 5857)
+++ /trunk/ppImage/src/ppImage.c	(revision 5857)
@@ -0,0 +1,34 @@
+# include "ppImage.h"
+
+int main (int argc, char **argv) {
+
+    ppData data;
+    ppConfig config;
+    ppRecipe options;
+
+    psTimerStart("phase2");
+
+    // Parse the configuration and arguments
+    ppConfigLoad (&config, argc, argv);
+
+    // Open the input image, output image, output mask
+    // Get camera configuration from header if not already defined
+    // Construct camera in preparation for reading
+    ppDataLoadMeta (&data, &config, argv);
+
+    // Set various tasks (define optional operations)
+    ppRecipeInit (&options, &config)
+
+    // Read in the input pixels
+    // Generate mask and weight frame
+    // Read in detrend pixels
+    ppDataLoadPixels (&data, &options, &config);
+
+    // Image Arithmetic Loop
+    ppImageLoop (&data, &options, &config);
+
+    // output options
+    ppImageOutput (&data, &options, &config);
+
+    exit (0);
+}
Index: /trunk/ppImage/src/ppImage.h
===================================================================
--- /trunk/ppImage/src/ppImage.h	(revision 5857)
+++ /trunk/ppImage/src/ppImage.h	(revision 5857)
@@ -0,0 +1,73 @@
+#include <stdio.h>
+#include <strings.h>
+
+#include "pslib.h"
+
+#include "psAdditionals.h"
+
+#include "pmAstrometry.h"
+#include "pmReadout.h"
+#include "pmConfig.h"
+#include "pmFPAConstruct.h"
+#include "pmFPARead.h"
+#include "pmFPAConceptsGet.h"
+#include "pmFPAWrite.h"
+
+#include "pmFlatField.h"
+#include "pmMaskBadPixels.h"
+#include "pmNonLinear.h"
+#include "pmSubtractBias.h"
+#include "pmChipMosaic.h"
+//#include "pmFPAMorph.h"
+
+#define RECIPE "PHASE2"                 // Name of the recipe to use
+
+typedef struct {
+    psMetadata *site;
+    psMetadata *camera;
+    psMetadata *recipe;
+    psMetadata *arguments;
+    psDB       *database;
+} ppConfig;
+
+typedef struct {
+    bool doMask;			// Mask bad pixles
+    bool doNonLin;			// Non-linearity correction
+    bool doBias;			// Bias subtraction
+    bool doDark;			// Dark subtraction
+    bool doOverscan;			// Overscan subtraction
+    bool doFlat;			// Flat-field normalisation
+    bool doFringe;			// Fringe subtraction
+    bool doSource;			// Source identification and photometry
+    bool doAstrom;			// Astrometry
+    int overscanBins;			// Number of pixels per bin for overscan
+    psStats *overscanStats;		// Statistics for overscan
+    void *overscanFit;			// Overscan fit (polynomial or spline)
+    pmFit overscanFitType;		// Fit type for overscan
+    pmOverscanAxis overscanMode;	// Axis for overscan
+
+    char *maskName;
+    char *biasName;
+    char *darkName;
+    char *flatName;
+} ppRecipe;
+
+typedef struct {
+    psMetadata *header;
+
+    pmFPA *input;
+    pmFPA *weight;
+    pmFPA *mask;
+    pmFPA *bias;
+    pmFPA *dark;
+    pmFPA *flat;
+    pmFPA *fringe;
+
+    psFits *inputFile;
+    psFits *weightFile;
+    psFits *maskFile;
+    psFits *biasFile;
+    psFits *darkFile;
+    psFits *flatFile;
+    psFits *fringeFile;
+} ppData;
Index: /trunk/ppImage/src/ppImageDetrendBias.c
===================================================================
--- /trunk/ppImage/src/ppImageDetrendBias.c	(revision 5857)
+++ /trunk/ppImage/src/ppImageDetrendBias.c	(revision 5857)
@@ -0,0 +1,87 @@
+# include "ppImage.h"
+
+bool ppDetrendBias (pmCell *input, psImage *bias, psImage *dark, ppRecipe *options) {
+
+    psPolynomial1D *poly;
+    psSpline1D *spline;
+
+    // XXX EAM : this is an odd trigger for this
+    if (options->overscanMode == PM_OVERSCAN_ROWS || options->overscanMode == PM_OVERSCAN_COLUMNS) {
+	// Need to get the read direction
+	int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR");
+	if (readdir == 1) {
+// psmodule-0.8.0 has confused PM_OVERSCAN_ROWS and PM_OVERSCAN_COLUMNS; bug 608.
+#ifdef PRODUCTION
+	    options->overscanMode = PM_OVERSCAN_ROWS;
+#else
+	    options->overscanMode = PM_OVERSCAN_COLUMNS;
+#endif
+	} else if (readdir == 2) {
+#ifdef PRODUCTION
+	    options->overscanMode = PM_OVERSCAN_COLUMNS;
+#else
+	    options->overscanMode = PM_OVERSCAN_ROWS;
+#endif
+	} else {
+	    psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n",
+			      readdir);
+	    options->overscanMode = PM_OVERSCAN_ROWS;
+	}
+    }
+
+    psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions
+#ifdef PRODUCTION
+    (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode,
+			 overscanStats, overscanBins, overscanFitType, pedestal);
+#else
+    (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode,
+			 overscanStats, overscanBins, overscanFitType, biasReadout);
+#endif
+    psFree(inputOverscans);
+
+    // Output overscan fit results, if required
+    if (!doOverscan) { 
+	return true;
+    }
+
+    if (! overscanFit) {
+	psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n");
+	return true;
+    } 
+
+    switch (overscanFitType) {
+      case PM_FIT_POLYNOMIAL:
+	poly = (psPolynomial1D*)overscanFit; // The polynomial
+	psString coeffs = NULL;     // String containing the coefficients
+	for (int i = 0; i < poly->nX; i++) {
+	    psStringAppend(&coeffs, "%e ", poly->coeff[i]);
+	}
+	psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n",
+		 coeffs);
+	psFree(coeffs);
+	break;
+
+      case PM_FIT_SPLINE:
+	spline = (psSpline1D*)overscanFit; // The spline
+	psString coeffs = NULL;     // String containing the coefficients
+	for (int i = 0; i < spline->n; i++) {
+	    psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
+	    psStringAppend(&coeffs, "%d: ", i);
+	    for (int j = 0; j < poly->nX; j++) {
+		psStringAppend(&coeffs, "%e ", poly->coeff[i]);
+	    }
+	    psStringAppend(&coeffs, "\n");
+	}
+	psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n",
+		 coeffs);
+	psFree(coeffs);
+	break;
+
+      case PM_FIT_NONE:
+	break;
+
+      default:
+	psAbort(__func__, "Should never get here!!!\n");
+    }
+    return true;
+}
Index: /trunk/ppImage/src/ppImageDetrendFlat.c
===================================================================
--- /trunk/ppImage/src/ppImageDetrendFlat.c	(revision 5857)
+++ /trunk/ppImage/src/ppImageDetrendFlat.c	(revision 5857)
@@ -0,0 +1,15 @@
+# include "ppImage.h"
+
+bool ppDetrendFlat () {
+
+    psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n");
+#ifndef PRODUCTION
+    psImage *dummyImage = psImageAlloc(inputReadout->image->numCols,
+				       inputReadout->image->numRows,
+				       PS_TYPE_U8);
+    pmReadout *dummyMask = pmReadoutAlloc(NULL);
+    dummyMask->image = dummyImage;
+    (void)pmFlatField(inputReadout, dummyMask, flatReadout);
+    psFree(dummyMask);
+    psFree(dummyImage);
+#endif
Index: /trunk/ppImage/src/ppImageDetrendMask.c
===================================================================
--- /trunk/ppImage/src/ppImageDetrendMask.c	(revision 5857)
+++ /trunk/ppImage/src/ppImageDetrendMask.c	(revision 5857)
@@ -0,0 +1,14 @@
+# include "ppImage.h"
+
+bool ppDetrendMask (pmCell *cell, pmReadout *readout, pmImage *mask) {
+
+    float saturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
+
+    // Need to change this later to grow the mask by the size of the convolution kernel
+    (void)pmMaskBadPixels(readout, mask, 
+			  PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation,
+			  PM_MASK_TRAP, 0);
+
+    return true;
+}
+
Index: /trunk/ppImage/src/ppImageDetrendNonLinear.c
===================================================================
--- /trunk/ppImage/src/ppImageDetrendNonLinear.c	(revision 5857)
+++ /trunk/ppImage/src/ppImageDetrendNonLinear.c	(revision 5857)
@@ -0,0 +1,133 @@
+# include "ppImage.h"
+
+bool ppDetrendNonLinear (pmCell *cell, pmReadout *readout, psMetadata *recipe) {
+
+    psString name;
+    psVector *coeff;
+    psMetadata *options;
+
+    psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA");
+    if (! dataItem) {
+	psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to "
+		 "find NONLIN.DATA in recipe --- ignored.\n");
+	return false;
+    } 
+
+    switch (dataItem->type) {
+      case PS_DATA_VECTOR:
+	// These are the polynomial coefficients
+	coeff = dataItem->data.V; // The coefficient vector
+	if (coeff->type.type != PS_TYPE_F64) {
+	    psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version
+	    coeff = temp;
+	}
+	psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1,
+							 PS_POLYNOMIAL_ORD);
+	psFree(correction->coeff);
+	correction->coeff = psMemIncrRefCounter(coeff->data.F64);
+	(void)pmNonLinearityPolynomial(inputReadout, correction);
+	psFree(coeff);
+	psFree(correction);
+	return true;
+
+      case PS_DATA_STRING:
+	// This is a filename
+	name = dataItem->data.V;       // Filename
+	psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0);
+	if (psLookupTableRead(table) <= 0) {
+	    psErrorStackPrint(stderr, "Unable to read non-linearity correction file "
+			      "%s --- ignored\n", name);
+	    return false;
+	} 
+#ifdef PRODUCTION
+	(void)pmNonLinearityLookup(inputReadout, table);
+#else
+	psVector *influx = table->values->data[0];
+	psVector *outflux = table->values->data[1];
+	(void)pmNonLinearityLookup(inputReadout, table->values->data[0],
+				   table->values->data[1]);
+#endif
+	psFree(table);
+	return true;
+
+      case PS_DATA_METADATA:
+	// This is a menu
+	options = dataItem->data.V; // Options with concept values as keys
+	bool mdok = false; // Success of MD lookup
+	psString concept = psMetadataLookupStr(&mdok, recipe, "NONLIN.SOURCE");
+	if (! mdok || ! concept) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
+		     "unable to find NONLIN.SOURCE in recipe --- ignored.\n");
+	    return false;
+	} 
+	psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept);
+	if (! conceptValueItem) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s "
+		     "for non-linearity correction --- ignored.\n", concept);
+	    return false;
+	} 
+	if (conceptValueItem->type != PS_DATA_STRING) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as"
+		     " expected for non-linearity correction --- ignored.\n",
+		     concept);
+	    return false;
+	} 
+
+	psString conceptValue = conceptValueItem->data.V;
+	// Get the value of the concept
+	psMetadataItem *optionItem = psMetadataLookup(options, conceptValue);
+	if (!optionItem) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA"
+		     " --- ignored.\n", conceptValue);
+	    return false;
+	} 
+	if (optionItem->type == PS_DATA_VECTOR) {
+	    // These are the polynomial coefficients
+	    coeff = optionItem->data.V; // The coefficient vector
+	    if (coeff->type.type != PS_TYPE_F64) {
+		psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64);
+		coeff = temp;
+	    }
+	    // Polynomial correction
+	    psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1,
+							     PS_POLYNOMIAL_ORD);
+	    psFree(correction->coeff);
+	    correction->coeff = psMemIncrRefCounter(coeff->data.F64);
+	    (void)pmNonLinearityPolynomial(inputReadout, correction);
+	    psFree(coeff);
+	    psFree(correction);
+	    return true;
+	} 
+	if (optionItem->type == PS_DATA_STRING) {
+	    // This is a filename
+	    psString tableName = optionItem->data.V; // The filename
+	    psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0);
+	    int numLines = 0; // Number of lines read from table
+	    if ((numLines = psLookupTableRead(table)) <= 0) {
+		psErrorStackPrint(stderr, "Unable to read non-linearity "
+				  "correction file %s --- ignored\n",
+				  tableName);
+		return false;
+	    } 
+#ifdef PRODUCTION
+	    (void)pmNonLinearityLookup(inputReadout, table);
+#else
+	    printf("XXX: Non-linearity correction from lookup table not "
+		   "yet implemented.\n");
+#endif
+	    psFree(table);
+	    return true;
+	}
+
+	psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction "
+		 "desired but unable to interpret NONLIN.DATA for %s"
+		 " --- ignored\n", conceptValue);
+	return false;
+  
+      default:
+	psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
+		 "NONLIN.DATA is of invalid type --- ignored.\n");
+    }
+    return false;
+}
+
Index: /trunk/ppImage/src/ppImageDetrendPedestal.c
===================================================================
--- /trunk/ppImage/src/ppImageDetrendPedestal.c	(revision 5857)
+++ /trunk/ppImage/src/ppImageDetrendPedestal.c	(revision 5857)
@@ -0,0 +1,36 @@
+# include "ppImage.h"
+
+psImage *ppDetrendPedestal (pmCell *input, psImage *bias, psImage *dark, ppRecipe *options) {
+
+    psImage *pedestal = NULL;       // Pedestal image (bias + scaled dark)
+
+#ifdef PRODUCTION
+    if (options->doDark) {
+	// Dark time for input image
+	float inputTime = psMetadataLookupF32(NULL, input->concepts, "CELL.DARKTIME");
+	if (inputTime <= 0) {
+	    psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n",
+			      inputTime);
+	    exit(EXIT_FAILURE);
+	}
+	pedestal = (psImage*)psBinaryOp(NULL, dark, "*", psScalarAlloc(inputTime * darkTime, PS_TYPE_F32));
+    }
+    if (options->doBias) {
+	if (pedestal) {
+	    if (bias->numRows != dark->numRows ||
+		bias->numCols != dark->numCols) {
+		psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: "
+			"%dx%d vs %dx%d\n", bias->numCols, bias->numRows,
+			dark->numCols, dark->numRows);
+		exit(EXIT_FAILURE);
+	    }
+	    (void)psBinaryOp(pedestal, pedestal, "+", bias);
+	} else {
+	    pedestal = psMemIncrRefCounter(bias);
+	}
+    }
+#endif
+
+    return pedestal;
+}
+
Index: /trunk/ppImage/src/ppImageLoop.c
===================================================================
--- /trunk/ppImage/src/ppImageLoop.c	(revision 5857)
+++ /trunk/ppImage/src/ppImageLoop.c	(revision 5857)
@@ -0,0 +1,39 @@
+# include "ppImage.h"
+
+bool ppImageLoop (ppData *data, ppRecipe *options, ppConfig *config) {
+
+    psArray *inputChips = data->input->chips; // Array of chips in input image
+    psArray *biasChips 	= data->bias->chips;  // Array of chips in bias image
+    psArray *darkChips 	= data->dark->chips;  // Array of chips in dark image
+    psArray *maskChips 	= data->mask->chips;  // Array of chips in mask image
+    psArray *flatChips 	= data->flat->chips;  // Array of chips in flat image
+
+    for (int i = 0; i < inputChips->n; i++) {
+        pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image
+        pmChip *biasChip  = biasChips->data[i];  // Chip of interest in bias image
+        pmChip *darkChip  = darkChips->data[i];	 // Chip of interest in dark image
+        pmChip *maskChip  = maskChips->data[i];	 // Chip of interest in mask image
+        pmChip *flatChip  = flatChips->data[i];  // Chip of interest in flat image
+
+        if (! inputChip->valid) { continue; }
+
+        psArray *inputCells = inputChip->cells; // Array of cells in input image
+        psArray *biasCells  = biasChip->cells; 	// Array of cells in bias image
+        psArray *darkCells  = darkChip->cells; 	// Array of cells in dark image
+        psArray *maskCells  = maskChip->cells; 	// Array of cells in mask image
+        psArray *flatCells  = flatChip->cells; 	// Array of cells in flat image
+
+        for (int j = 0; j < inputCells->n; j++) {
+            pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image
+            pmCell *biasCell  = biasCells->data[j];  // Cell of interest in bias image
+            pmCell *darkCell  = darkCells->data[j];  // Cell of interest in dark imag
+            pmCell *maskCell  = maskCells->data[j];  // Cell of interest in mask image
+            pmCell *flatCell  = flatCells->data[j];  // Cell of interest in flat image
+
+            if (! inputCell->valid) { continue; }
+
+	    ppReadoutDetrend (inputCell, maskCell, biasCell, darkCell, flatCell);
+	}
+    }
+    return true;
+}
Index: /trunk/ppImage/src/ppImageOutput.c
===================================================================
--- /trunk/ppImage/src/ppImageOutput.c	(revision 5857)
+++ /trunk/ppImage/src/ppImageOutput.c	(revision 5857)
@@ -0,0 +1,29 @@
+# include "ppImage.h"
+
+bool ppImageLoop (ppData *data, ppRecipe *options, ppConfig *config) {
+
+    // Write the output
+    pmFPAWrite(outputFile, input, database);
+    pmFPAWriteMask(input, outputFile);
+    pmFPAWriteWeight(input, outputFile);
+
+    psLogMsg("phase2", PS_LOG_INFO, "Output completed after %f sec.\n", psTimerMark("phase2"));
+
+    psFitsClose(outputFile);
+
+    psFree(arguments);
+    psFree(site);
+    psFree(header);
+    psFree(camera);
+    psFree(recipe);
+    psFree(input);
+    psFree(mask);
+    psFree(bias);
+    psFree(dark);
+    psFree(flat);
+    psFree(overscanFit);
+    psFree(overscanStats);
+
+    psLogMsg("phase2", PS_LOG_INFO, "Output completed after %f sec.\n", psTimerMark("phase2"));
+    return true;
+}
Index: /trunk/ppImage/src/ppReadoutDetrend.c
===================================================================
--- /trunk/ppImage/src/ppReadoutDetrend.c	(revision 5857)
+++ /trunk/ppImage/src/ppReadoutDetrend.c	(revision 5857)
@@ -0,0 +1,81 @@
+# include "ppImage.h"
+
+bool ppReadoutDetrend (pmCell *input, pmCell *mask, pmCell *bias, pmCell *dark, pmCell *flat, ppRecipe *options, ppConfig *config) {
+
+    psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image
+
+    pmReadout *maskReadout = NULL; // Mask readout
+    psImage *maskImage = NULL;  // Mask pixels
+    if (doMask) {
+	maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image
+	maskImage = maskReadout->image;
+	if (maskCell->readouts->n > 1) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be used.");
+	}
+    }
+
+    pmReadout *biasReadout = NULL; // Bias readout
+    psImage *biasImage = NULL;  // Bias pixels
+    if (doBias) {
+	biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image
+	biasImage = biasReadout->image;
+	if (biasCell->readouts->n > 1) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be used.");
+	}
+    }
+
+    pmReadout *darkReadout = NULL; // Readout of interest in dark image
+    psImage *darkImage = NULL;  // Dark pixels
+    float darkTime = 1.0;       // Dark time for dark image
+    if (doDark) {
+	pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image
+	darkImage = darkReadout->image;
+	darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME");
+	if (darkTime <= 0.0) {
+	    psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime);
+	    exit(EXIT_FAILURE);
+	}
+	if (darkCell->readouts->n > 1) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be used.");
+	}
+    }
+
+    pmReadout *flatReadout = NULL; // Flat readout
+    psImage *flatImage = NULL;  // Flat pixels
+    if (doFlat) {
+	flatReadout = flatCell->readouts->data[0]; // Readout of interest in flat image
+	flatImage = flatReadout->image;
+	if (flatCell->readouts->n > 1) {
+	    psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be"
+		     " used.");
+	}
+    }
+
+    for (int k = 0; k < inputReadouts->n; k++) {
+	pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image
+	psImage *inputImage = inputReadout->image; // The actual image
+
+	// Mask bad pixels
+	if (options->doMask) {
+	    ppDetrendMask (inputCell, inputReadout, maskImage);
+	}
+
+	// Non-linearity correction
+	if (options->doNonLin) {
+	    ppDetrendNonLinear (inputCell, inputReadout, config->recipe);
+	}
+
+	// Changed the definition of pmSubtractBias to do the dark subtraction as well, but
+	// it's not in there yet.  Once it is, we can get rid of "pedestal".
+	// generate the pedestal image
+	psImage *pedestal = ppDetrendPedestal (inputCell, biasImage, darkImage, options);
+
+	// Bias, dark and overscan subtraction are all merged.
+	if (options->doBias || options->doOverscan || options->doDark) {
+	    ppDetrendBias (inputCell, pedestal, options);
+	}
+
+	// Flat-field correction
+	if (doFlat) {
+	    ppDetrendFlat (inputCell, flatImage, options);
+	}
Index: /trunk/ppImage/src/ppReciptInit.c
===================================================================
--- /trunk/ppImage/src/ppReciptInit.c	(revision 5857)
+++ /trunk/ppImage/src/ppReciptInit.c	(revision 5857)
@@ -0,0 +1,125 @@
+# include "ppImage.h"
+
+bool ppRecipeInit (ppRecipe *options, ppConfig *config) {
+
+    // default initial values
+    options->doMask = false;		// Mask bad pixles
+    options->doNonLin = false;		// Non-linearity correction
+    options->doBias = false;		// Bias subtraction
+    options->doDark = false;		// Dark subtraction
+    options->doOverscan = false;	// Overscan subtraction
+    options->doFlat = false;		// Flat-field normalisation
+    options->doFringe = false;		// Fringe subtraction
+    options->doSource = false;		// Source identification and photometry
+    options->doAstrom = false;		// Astrometry
+    options->overscanBins = 1;		// Number of pixels per bin for overscan
+    options->overscanStats = NULL;      // Statistics for overscan
+    options->overscanFit = NULL;	// Overscan fit (polynomial or spline)
+    options->overscanFitType = PM_FIT_NONE; // Fit type for overscan
+    options->overscanMode = PM_OVERSCAN_NONE; // Axis for overscan
+
+    // Chip selection
+    options->chipNum  = psMetadataLookupS32(NULL, config->arguments, "-chip"); // Chip number to work on
+    
+    // mask - recipe options
+    options->maskName = psMetadataLookupStr(NULL, config->arguments, "-mask");
+    if (psMetadataLookupBool(NULL, config->recipe, "MASK")) {
+        if (strlen(options->maskName) > 0) {
+            options->doMask = true;
+        } else {
+            psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied"
+		     " --- no masking will be performed.\n");
+        }
+    }
+
+    // non-linear - recipe options
+    if (psMetadataLookupBool(NULL, config->recipe, "NONLIN")) {
+        options->doNonLin = true;
+    }
+
+    // bias - recipe options
+    options->biasName = = psMetadataLookupStr(NULL, config->arguments, "-bias");
+    if (psMetadataLookupBool(NULL, config->recipe, "BIAS")) {
+        if (strlen(biasName) > 0) {
+            options->doBias = true;
+        } else {
+            psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- "
+                     "no bias subtraction will be performed.\n");
+        }
+    }
+
+    // dark - recipe options
+    options->darkName = psMetadataLookupStr(NULL, config->arguments, "-dark");
+    if (psMetadataLookupBool(NULL, config->recipe, "DARK")) {
+        if (strlen(options->darkName) > 0) {
+            options->doDark = true;
+        } else {
+            psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- "
+                     "no dark subtraction will be performed.\n");
+        }
+    }
+
+    // overscan - recipe options
+    if (psMetadataLookupBool(NULL, config->recipe, "OVERSCAN")) {
+        options->doOverscan = true;
+        psString mode = psMetadataLookupStr(NULL, config->recipe, "OVERSCAN.MODE");
+        if (strcasecmp(mode, "INDIVIDUAL") == 0) {
+            options->overscanMode = PM_OVERSCAN_ROWS;
+            // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip
+        } else if (strcasecmp(mode, "ALL") == 0) {
+            options->overscanMode = PM_OVERSCAN_ALL;
+        } else if (strcasecmp(mode, "NONE") != 0) {
+            psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:"
+                     " assuming NONE.\n", mode);
+        }
+        psString fit = psMetadataLookupStr(NULL, config->recipe, "OVERSCAN.FIT");
+        if (strcasecmp(fit, "POLYNOMIAL") == 0) {
+            options->overscanFitType = PM_FIT_POLYNOMIAL;
+            int order = psMetadataLookupS32(NULL, config->recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
+            options->overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
+        } else if (strcasecmp(fit, "SPLINE") == 0) {
+            options->overscanFitType = PM_FIT_SPLINE;
+            int order = psMetadataLookupS32(NULL, config->recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
+            options->overscanFit = NULL;
+	    // XXX : not in psLib yet : options->overscanFit = psSpline1DAlloc();
+        } else if (strcasecmp(fit, "NONE") != 0) {
+            psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:"
+                     " assuming NONE.\n", fit);
+        }
+        options->overscanBins = psMetadataLookupS32(NULL, config->recipe, "OVERSCAN.BIN");
+        if (options->overscanBins <= 0) {
+            psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", options->overscanBins);
+            options->overscanBins = 1;
+        }
+        psString stat = psMetadataLookupStr(NULL, config->recipe, "OVERSCAN.STAT");
+        if (strcasecmp(stat, "MEAN")) {
+            options->overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+        } else if (strcasecmp(stat, "MEDIAN")) {
+            options->overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+        } else {
+            psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat);
+            options->overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+        }
+    }
+
+    // flat-field - recipe options
+    options->flatName = psMetadataLookupStr(NULL, config->arguments, "-flat");
+    if (psMetadataLookupBool(NULL, config->recipe, "FLAT")) {
+        if (strlen(options->flatName) > 0) {
+            options->doFlat = true;
+        } else {
+            psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- "
+                     "no flat-fielding will be performed.\n");
+        }
+    }
+
+    // XXX need to add the following:
+
+    // fringe - recipe options
+
+    // photom - recipe options
+
+    // astrom - recipe options
+
+    return true;
+}
