#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"

// Phase 2 needs to:
// * Read configurations
// * Read in the images
// * Flag bad pixels
// * Non-linearity correction
// * Bias/dark subtraction (MHPCC)
// * Flat-field (MHPCC)
// * Fringe subtraction (Not for now)
// * Source identification and photometry (IPP prototype)
// * Astrometry (IPP prototype soon)
// * Write calibrated image
// * Write source catalogue, astrometry, etc.

// Currently neglecting different input types (F32, S32, U16, etc)

// Would like to fold in Nebulous at some stage as well.

// Currently neglect to convolve the reference frames with the orthogonal transfer kernel

// phase2 INPUT.fits OUTPUT.fits -bias BIAS.fits -dark DARK.fits -flat FLAT.fits -chip CHIP
//
// Most are self-explanatory.  "-chip" says "only work on this particular chip".


#define RECIPE "PHASE2"                 // Name of the recipe to use

static psMemId memPrintAlloc(const psMemBlock *mb)
{
    printf("Allocated memory block %ld (%ld):\n"
           "\tFile %s, line %d, size %d\n"
           "\tPosts: %x %x %x\n",
           mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
           *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
    return 0;
}
static psMemId memPrintFree(const psMemBlock *mb)
{
    printf("Freed memory block %ld (%ld):\n"
           "\tFile %s, line %d, size %d\n"
           "\tPosts: %x %x %x\n",
           mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
           *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
    return 0;
}
static void memPrint(const psPtr ptr)
{
    psMemBlock *mb = ((psMemBlock*)ptr) - 1;
    printf("Memory block %ld (%ld):\n"
           "\tFile %s, line %d, size %d\n"
           "\tPosts: %x %x %x\n",
           mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
           *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
}

int main(int argc, char *argv[])
{
#if 0
    // Hunting memory leaks
    psMemAllocCallbackSetID(22367);
    psMemFreeCallbackSetID(22367);
    psMemAllocCallbackSet(memPrintAlloc);
    psMemFreeCallbackSet(memPrintFree);
#endif

    psTimerStart("phase2");

    // Parse the configurations
    psMetadata *site = NULL;            // Site configuration
    psMetadata *camera = NULL;          // Camera configuration
    psMetadata *recipe = NULL;          // Recipe configuration
    if (! pmConfigRead(&site, &camera, &recipe, &argc, argv, RECIPE)) {
        psErrorStackPrint(stderr, "Can't find site configuration!\n");
        exit(EXIT_FAILURE);
    }

    // Parse other command-line arguments
    psMetadata *arguments = psMetadataAlloc(); // The arguments, with default values
    psMetadataAddStr(arguments, PS_LIST_TAIL, "-bias", 0, "Name of the bias image", "");
    psMetadataAddStr(arguments, PS_LIST_TAIL, "-dark", 0, "Name of the dark image", "");
    psMetadataAddStr(arguments, PS_LIST_TAIL, "-flat", 0, "Name of the flat-field image", "");
    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", 0, "Name of the mask image", "");
    psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", 0, "Chip number to process (if positive)", -1);
    if (! psArgumentParse(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(arguments);
        psFree(arguments);
        exit(EXIT_FAILURE);
    }
    const char *inputName = argv[1];    // Name of input image
    const char *outputName = argv[2];   // Name of output image
    const char *biasName = psMetadataLookupStr(NULL, arguments, "-bias"); // Name of bias image
    const char *darkName = psMetadataLookupStr(NULL, arguments, "-dark"); // Name of dark image
    const char *flatName = psMetadataLookupStr(NULL, arguments, "-flat"); // Name of flat-field image
    const char *maskName = psMetadataLookupStr(NULL, arguments, "-mask"); // Name of mask image
    const int chipNum = psMetadataLookupS32(NULL, arguments, "-chip"); // Chip number to work on
    printf("Input: %s\nOutput: %s\n", inputName, outputName);
    printf("Bias: %s\n", biasName);
    printf("Dark: %s\n", darkName);
    printf("Flat: %s\n", flatName);
    printf("Mask: %s\n", maskName);
    printf("Chip: %d\n", chipNum);

    // Open the input
    psFits *inputFile = psFitsOpen(inputName, "r"); // File handle for FITS file
    if (! inputFile) {
        psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName);
        exit(EXIT_FAILURE);
    }
    psMetadata *header = psFitsReadHeader(NULL, inputFile); // FITS header
#if PRODUCTION
    psDB *database = pmConfigDB(site);  // Database handle
#else
    psDB *database = NULL;              // Database handle
#endif

    // Open the output and output mask
    // We do it here so that we don't process the whole lot and then find out we can't open the output file
    psFits *outputFile = psFitsOpen(outputName, "w");
    if (! outputFile) {
        psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName);
        exit(EXIT_FAILURE);
    }

    // Get camera configuration from header if not already defined
    if (! camera) {
        camera = pmConfigCameraFromHeader(site, header);
        if (! camera) {
            psErrorStackPrint(stderr, "Can't find camera configuration!\n");
            exit(EXIT_FAILURE);
        }
   } else if (! pmConfigValidateCamera(camera, header)) {
        psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName);
        exit(EXIT_FAILURE);
    }
    if (! recipe && !(recipe = pmConfigRecipeFromCamera(camera, RECIPE))) {
        psErrorStackPrint(stderr, "Can't find recipe configuration!\n");
        exit(EXIT_FAILURE);
    }

    // Construct camera in preparation for reading
    pmFPA *input = pmFPAConstruct(camera);
    pmFPA *mask = pmFPAConstruct(camera);
    pmFPA *bias = pmFPAConstruct(camera);
    pmFPA *dark = pmFPAConstruct(camera);
    pmFPA *flat = pmFPAConstruct(camera);

    // Set various tasks
    bool doMask = false;                // Mask bad pixles
    bool doNonLin = false;              // Non-linearity correction
    bool doBias = false;                // Bias subtraction
    bool doDark = false;                // Dark subtraction
    bool doOverscan = false;            // Overscan subtraction
    bool doFlat = false;                // Flat-field normalisation
    bool doFringe = false;              // Fringe subtraction
    bool doSource = false;              // Source identification and photometry
    bool doAstrom = false;              // Astrometry
    pmOverscanAxis overscanMode = PM_OVERSCAN_NONE; // Axis for overscan
    pmFit overscanFitType = PM_FIT_NONE; // Fit type for overscan
    int overscanBins = 1;               // Number of pixels per bin for overscan
    psStats *overscanStats = NULL;      // Statistics for overscan
    void *overscanFit = NULL;           // Overscan fit (polynomial or spline)

    if (psMetadataLookupBool(NULL, recipe, "MASK")) {
        if (strlen(maskName) > 0) {
            doMask = true;
        } else {
            psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking "
                     "will be performed.\n");
        }
    }
    if (psMetadataLookupBool(NULL, recipe, "NONLIN")) {
        doNonLin = true;
    }
    if (psMetadataLookupBool(NULL, recipe, "BIAS")) {
        if (strlen(biasName) > 0) {
            doBias = true;
        } else {
            psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- "
                     "no bias subtraction will be performed.\n");
        }
    }
    if (psMetadataLookupBool(NULL, recipe, "DARK")) {
        if (strlen(darkName) > 0) {
            doDark = true;
        } else {
            psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- "
                     "no dark subtraction will be performed.\n");
        }
    }
    if (psMetadataLookupBool(NULL, recipe, "OVERSCAN")) {
        doOverscan = true;
        psString mode = psMetadataLookupStr(NULL, recipe, "OVERSCAN.MODE");
        if (strcasecmp(mode, "INDIVIDUAL") == 0) {
            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) {
            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, recipe, "OVERSCAN.FIT");
        if (strcasecmp(fit, "POLYNOMIAL") == 0) {
            overscanFitType = PM_FIT_POLYNOMIAL;
            int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
            overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
        } else if (strcasecmp(fit, "SPLINE") == 0) {
            overscanFitType = PM_FIT_SPLINE;
            int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
            overscanFit = NULL;
//          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);
        }
        overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN");
        if (overscanBins <= 0) {
            psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins);
            overscanBins = 1;
        }
        psString stat = psMetadataLookupStr(NULL, recipe, "OVERSCAN.STAT");
        if (strcasecmp(stat, "MEAN")) {
            overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
        } else if (strcasecmp(stat, "MEDIAN")) {
            overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
        } else {
            psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat);
            overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
        }
    }

    if (psMetadataLookupBool(NULL, recipe, "FLAT")) {
        if (strlen(flatName) > 0) {
            doFlat = true;
        } else {
            psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- "
                     "no flat-fielding will be performed.\n");
        }
    }

    // Chip selection
    if (chipNum >= 0) {
        if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) ||
            ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) ||
            ! pmFPASelectChip(mask, chipNum)) {
            psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum);
            exit(EXIT_FAILURE);
        }
        psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum);
    }

    psLogMsg("phase2", PS_LOG_INFO, "Setup completed after %f sec.\n", psTimerMark("phase2"));

    // Read in the input pixels
    psLogMsg("phase2", PS_LOG_INFO, "Opening input image: %s\n", inputName);
    if (! pmFPARead(input, inputFile, header, database)) {
        psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName);
        exit(EXIT_FAILURE);
    }
    psFitsClose(inputFile);

#if 0
    pmFPAReadMask(input, inputFile);
    pmFPAReadWeight(input, inputFile);
//#else
    {
        // Generate mask and weight frame
        p_pmHDU *hdu = NULL;
        if (input->hdu) {
            hdu = input->hdu;
        }
        psArray *chips = input->chips;
        for (int chipNum = 0; chipNum < chips->n; chipNum++) {
            pmChip *chip = chips->data[chipNum];
            if (chip->valid) {
                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) {
                        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];
                            for (int j = 0; j < image->numRows; j++) {
                                for (int i = 0; i < image->numCols; i++) {
                                    mask->data.U8[j][i] = j + i;
                                    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);
                        }
                    }
                }
            }
        }
    }
#endif

#if 0
    pmFPAPrint(input);
#endif

    // Load the calibration frames, if required
    if (doBias) {
        psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias
        psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias
        if (! pmConfigValidateCamera(camera, biasHeader)) {
            psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName);
            exit(EXIT_FAILURE);
        }
        psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName);
        if (! pmFPARead(bias, biasFile, biasHeader, database)) {
            psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName);
            exit(EXIT_FAILURE);
        }
        psFree(biasHeader);
        psFitsClose(biasFile);
    }

    if (doDark) {
        psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark
        psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark
        if (! pmConfigValidateCamera(camera, darkHeader)) {
            psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName);
            exit(EXIT_FAILURE);
        }
        psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName);
        if (! pmFPARead(dark, darkFile, darkHeader, database)) {
            psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName);
            exit(EXIT_FAILURE);
        }
        psFree(darkHeader);
        psFitsClose(darkFile);
    }

    if (doMask) {
        psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask
        psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask
        if (! pmConfigValidateCamera(camera, maskHeader)) {
            psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName);
            exit(EXIT_FAILURE);
        }
        psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName);
        if (! pmFPARead(mask, maskFile, maskHeader, database)) {
            psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName);
            exit(EXIT_FAILURE);
        }
        psFree(maskHeader);
        psFitsClose(maskFile);
    }

    if (doFlat) {
        psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat
        if (! flatFile) {
            psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName);
            exit(EXIT_FAILURE);
        }
        psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat
        if (! pmConfigValidateCamera(camera, flatHeader)) {
            psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName);
            exit(EXIT_FAILURE);
        }
        psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName);
        if (! pmFPARead(flat, flatFile, flatHeader, database)) {
            psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName);
            exit(EXIT_FAILURE);
        }
        psFree(flatHeader);
        psFitsClose(flatFile);
    }

    psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2"));

    psArray *inputChips = input->chips; // Array of chips in input image
    psArray *biasChips = bias->chips;   // Array of chips in bias image
    psArray *darkChips = dark->chips;   // Array of chips in dark image
    psArray *maskChips = mask->chips;   // Array of chips in mask image
    psArray *flatChips = 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;
            }

            psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image
            if (doBias && biasCell->readouts->n > 1) {
                psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be"
                         " used.");
            }
            if (doDark && darkCell->readouts->n > 1) {
                psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be"
                         " used.");
            }
            if (doMask && maskCell->readouts->n > 1) {
                psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be"
                         " used.");
            }
            if (doFlat && flatCell->readouts->n > 1) {
                psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be"
                         " used.");
            }

            pmReadout *biasReadout = NULL; // Bias readout
            pmReadout *maskReadout = NULL; // Mask readout
            pmReadout *flatReadout = NULL; // Flat readout
            psImage *biasImage = NULL;  // Bias pixels
            psImage *darkImage = NULL;  // Dark pixels
            float darkTime = 1.0;       // Dark time for dark image
            psImage *maskImage = NULL;  // Mask pixels

            if (doBias) {
                biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image
                biasImage = biasReadout->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 (doMask) {
                maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image
                maskImage = maskReadout->image;
            }
            if (doFlat) {
                flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image
            }

            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 (doMask) {
                    float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION");

                    // Need to change this later to grow the mask by the size of the convolution kernel
                    (void)pmMaskBadPixels(inputReadout, maskImage,
                                          PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation,
                                          PM_MASK_TRAP, 0);
                }

                // Non-linearity correction
                if (doNonLin) {
                    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");
                    } else {
                        switch (dataItem->type) {
                          case PS_DATA_VECTOR:
                            {
                                // These are the polynomial coefficients
                                psVector *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);
                            }
                            break;
                          case PS_DATA_STRING:
                            {
                                // This is a filename
                                psString 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);
                                } else {
#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);
                            }
                            break;
                          case PS_DATA_METADATA:
                            {
                                // This is a menu
                                psMetadata *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");
                                } else {
                                    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);
                                    } else 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);
                                    } else {
                                        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);
                                        } else if (optionItem->type == PS_DATA_VECTOR) {
                                            // These are the polynomial coefficients
                                            psVector *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);
                                        } else 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);
                                            } else {
#ifdef PRODUCTION
                                                (void)pmNonLinearityLookup(inputReadout, table);
#else
                                                printf("XXX: Non-linearity correction from lookup table not "
                                                       "yet implemented.\n");
#endif
                                            }
                                            psFree(table);
                                        } else {
                                            psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction "
                                                     "desired but unable to interpret NONLIN.DATA for %s"
                                                     " --- ignored\n", conceptValue);
                                        }
                                    }
                                }
                            }
                            break;
                          default:
                            psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
                                     "NONLIN.DATA is of invalid type --- ignored.\n");
                        }
                    }
                } // Non-linearity correction

                // Bias, dark and overscan subtraction are all merged.

                // 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".

#ifdef PRODUCTION
                psImage *pedestal = NULL;       // Pedestal image (bias + scaled dark)
                if (doDark) {
                    // Dark time for input image
                    float inputTime = psMetadataLookupF32(NULL, inputCell->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, darkImage, "*",
                                                    psScalarAlloc(inputTime * darkTime, PS_TYPE_F32));
                }
                if (doBias) {
                    if (pedestal) {
                        if (biasImage->numRows != darkImage->numRows ||
                            biasImage->numCols != darkImage->numCols) {
                            psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: "
                                    "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows,
                                    darkImage->numCols, darkImage->numRows);
                            exit(EXIT_FAILURE);
                        }
                        (void)psBinaryOp(pedestal, pedestal, "+", biasImage);
                    } else {
                        pedestal = psMemIncrRefCounter(biasImage);
                    }
                }
#endif
                if (overscanMode == PM_OVERSCAN_ROWS || 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
                        overscanMode = PM_OVERSCAN_ROWS;
#else
                        overscanMode = PM_OVERSCAN_COLUMNS;
#endif
                    } else if (readdir == 2) {
#ifdef PRODUCTION
                        overscanMode = PM_OVERSCAN_COLUMNS;
#else
                        overscanMode = PM_OVERSCAN_ROWS;
#endif
                    } else {
                        psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n",
                                          readdir);
                        overscanMode = PM_OVERSCAN_ROWS;
                    }
                }

                if (doBias || doOverscan || doDark) {
                    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) {
                        if (! overscanFit) {
                            psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n");
                        } else {
                            switch (overscanFitType) {
                              case PM_FIT_POLYNOMIAL:
                                {
                                    psPolynomial1D *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:
                                {
                                    psSpline1D *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");
                            }
                        }
                    }
                }

                // XXX: Temporary: until the other functions are altered to do this themselves.
                // Trim, so that flat, fringe etc computations are faster.
#if 0
#ifndef PRODUCTION
                psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC");
                psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion);
                psFree(inputReadout->image);
                inputReadout->image = trimmed;
#endif
#endif
                // Flat-field correction
                if (doFlat) {
                    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
                }

                // Fringe subtraction
                if (doFringe) {
#if 0
                    // Placeholder
                    pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask,
                                             psU8 maskVal, int binFactor, psStats *stats, float clipSD);
#endif
                }

                // Source identification and photometry
                if (doSource) {
                }

                // Astrometry
                if (doAstrom) {
                }


            } // Iterating over readouts
        } // Iterating over cells

#if 0
        {
            // Testing pmChipMosaic
            psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip
            psFits *mosaicFile = psFitsOpen("mosaic.fits", "w"); // FITS file to which to write
            psFitsWriteImage(mosaicFile, NULL, mosaic, 0);
            psFitsClose(mosaicFile);
            psFree(mosaic);
        }
#endif

    } // Iterating over chips

    psLogMsg("phase2", PS_LOG_INFO, "Processing completed after %f sec.\n", psTimerMark("phase2"));


#if 0
    {
        // Testing pmFPAMorph
        int nBad = 0;
        psMetadata *newCamera = psMetadataConfigParse(NULL, &nBad, "megacam_splice.config", true);
        pmFPA *newFPA = pmFPAConstruct(newCamera);
        pmFPAMorph(newFPA, input, true, 0, 0);
        //pmFPAPrint(newFPA);
        psFits *newFile = psFitsOpen("morph.fits", "w");
        pmFPAWrite(newFile, newFPA, database);
        psFitsClose(newFile);
        psFree(newFPA);
    }
#endif

#if 1
    // Write the output
    pmFPAWrite(outputFile, input, database);
    pmFPAWriteMask(input, outputFile);
    pmFPAWriteWeight(input, outputFile);
#endif

    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, "Cleanup completed after %f sec.\n", psTimerMark("phase2"));
    psTimerStop();

#if 0
    psMemCheckCorruption(true);
    psMemBlock **leaks = NULL;          // List of leaks
    int nLeaks = psMemCheckLeaks(0, &leaks, NULL, false); // Number of leaks
    printf("%d leaks found.\n", nLeaks);
#if 1
    for (int i = 0; i < nLeaks; i++) {
        psMemBlock *mb = leaks[i];
        printf("Memory leak %ld (%ld):\n"
               "\tFile %s, line %d, size %d\n"
               "\tPosts: %x %x %x\n",
               mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
               *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
    }
#endif
#endif

}
