/** @file ppSubMatchPSFs.c
 *
 *  @brief
 *
 *  @ingroup ppSub
 *
 *  @author IfA
 *  @version $Revision: 1.3 $ $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 <psphot.h>

#include "ppSub.h"

bool ppSubMatchPSFs(pmConfig *config, const pmFPAview *view)
{
    psAssert(config, "Require configuration");
    psAssert(view, "Require view");

    // Look up recipe values
    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    psAssert(recipe, "We checked this earlier, so it should be here.");

    // Input images
    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout

    // Output image holders
    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved
    if (!inConv) {
        pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Cell for convolved input
        inConv = pmReadoutAlloc(cell);
    }
    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved
    if (!refConv) {
        pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Cell for convolved ref.
        refConv = pmReadoutAlloc(cell);
    }

    bool mdok;                          // Status of MD lookup

    // Sources in image, used for stamps: these must be loaded from previous analysis stages
    pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
    psArray *sources = NULL;            // Sources in image; used for stamps
    if (sourcesRO) {
        sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
    }

    int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
    int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches
    float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
    float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
    float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
    const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps

    const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type
    psAssert(typeStr, "We put it here in ppSubArguments.c");
    pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel
    if (type == PM_SUBTRACTION_KERNEL_NONE) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr);
        return false;
    }

    int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size
    int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
    psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
    psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders

    int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
    int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
    int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
    float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness
    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
    float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
    float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel

    float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"

    // Options which control the generation of optimal parameters
    bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
    float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
    float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
    float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
    psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
    if (optimum) {
        optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
    }
    int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
    float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search

    // XXX Do we need/want to define different values for BAD and POOR subtraction vs BAD and POOR warp?
    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs
    psImageMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); // Bits to mask for poor pixels
    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels

    bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution?
    pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode

    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
    if (threads > 0) {
        pmSubtractionThreadsInit(inRO, refRO);
    }

    // Match the PSFs
    if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, spacing, threshold,
                            sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
                            binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, sys,
                            maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
        return false;
    }

    psFree(optWidths);

    pmSubtractionThreadsFinalize(inRO, refRO);

    psImageCovarianceTransfer(inConv->variance, inConv->covariance);
    psImageCovarianceTransfer(refConv->variance, refConv->covariance);

    // XXX drop the pixels associated with inRO and refRO (now that we have inConv and refConf)

    psphotSaveImage (NULL, inRO->image, "inRO.fits");
    psphotSaveImage (NULL, refRO->image, "refRO.fits");
    psphotSaveImage (NULL, inConv->image, "inConv.fits");
    psphotSaveImage (NULL, refConv->image, "refConv.fits");

    return true;
}
