/** @file ppSubReadoutSubtract.c
 *
 *  @brief
 *
 *  @ingroup ppSub
 *
 *  @author IfA
 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-18 00:31:20 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>

#include "ppSub.h"

bool ppSubReadoutSubtract(pmConfig *config)
{
    psAssert(config, "Require configuration");

    // Look up recipe values
    bool mdok = false;                  // Status of MD lookup
    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
    psAssert(recipe, "We checked this earlier, so it should be here.");

    bool reverse = psMetadataLookupBool(&mdok, config->arguments, "REVERSE"); // Reverse sense of subtraction?

    pmFPAview *view = ppSubViewReadout(); // View to readout

    // Subtraction is: minuend - subtrahend
    pmReadout *minuend = NULL;          // Positive image
    pmReadout *subtrahend = NULL;       // Negative image
    if (reverse) {
        minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
        subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
    } else {
        minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
        subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
    }

// XXX this should be replaced by SAVE options to the inConv and refConf pmFPAfiles
#ifdef TESTING
    {
        pmReadoutMaskApply(minuend, maskVal);
        psFits *fits = psFitsOpen("minuend.fits", "w");
        psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
        psFitsClose(fits);
    }
    {
        pmReadoutMaskApply(subtrahend, maskVal);
        psFits *fits = psFitsOpen("subtrahend.fits", "w");
        psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
            psFitsClose(fits);
    }
#endif

    // Do the actual subtraction
    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
    outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
    if (minuend->variance && subtrahend->variance) {
        outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
    }
    outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);

    psArray *covars = psArrayAlloc(2);  // Covariance pseudo-matrices
    covars->data[0] = psMemIncrRefCounter(minuend->covariance);
    covars->data[1] = psMemIncrRefCounter(subtrahend->covariance);
    outRO->covariance = psImageCovarianceAverage(covars);
    psFree(covars);

    outRO->data_exists = true;
    outRO->parent->data_exists = true;
    outRO->parent->parent->data_exists = true;

    pmSubtractionVisualShowSubtraction(minuend->image, subtrahend->image, outRO->image);

    // Copy concepts from the input to the output
    pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT"); // Input file
    pmFPA *inFPA = inFile->fpa;         // Input FPA
    pmFPAfile *outFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); // Output file
    pmFPA *outFPA = outFile->fpa;       // Output FPA
    if (!pmConceptsCopyFPA(outFPA, inFPA, true, true)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
        psFree(outRO);
        psFree(view);
        return false;
    }

    // Copy astrometry over
    // It should find its way into the output images and photometry
    pmHDU *inHDU = inFPA->hdu;          // Input HDU
    pmHDU *outHDU = outFPA->hdu;        // Output HDU
    pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSUB.OUTPUT"); // Output chip
    psFree(view);
    if (!outHDU || !inHDU) {
        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU at FPA level to copy astrometry.");
        return false;
    }
    if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA.");
        return false;
    }
    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
        return false;
    }

    return true;
}
