/** @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(ppSubData *data)
{
    psAssert(data, "Require processing data");
    pmConfig *config = data->config;    // Configuration
    psAssert(config, "Require configuration");

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

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

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

    // Load pre-calculated kernel, if available
    pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel

    psFree(view);

    // Sources in image, used for stamps: these must be loaded from previous analysis stages
    pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES");
    pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES");
    psArray *inSources = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES") :
        NULL; // Source list from input image
    psArray *refSources = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES") :
        NULL ; // Source list from reference image

    psArray *sources = NULL;            // Merged list of sources
    if (inSources && refSources) {
        float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Matching radius
        psArray *lists = psArrayAlloc(2); // Source lists
        lists->data[0] = psMemIncrRefCounter(inSources);
        lists->data[1] = psMemIncrRefCounter(refSources);
        sources = pmSourceMatchMerge(lists, radius);
        psFree(lists);
        if (!sources) {
            psError(PS_ERR_UNKNOWN, false, "Unable to merge source lists");
            return false;
        }
    } else if (inSources) {
        sources = psMemIncrRefCounter(inSources);
    } else if (refSources) {
        sources = psMemIncrRefCounter(refSources);
    }

    psMetadataAddArray(inConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE,
                       "Merged source list", sources);
    psMetadataAddArray(refConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE,
                       "Merged source list", sources);
    psFree(sources);                    // Drop reference

    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 *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

    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs
    psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", 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
    bool success = false;               // Operation was successful?
    if (kernelRO) {
        success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis,
                                            stride, sys, maskVal, maskBad, maskPoor, poorFrac, badFrac);
    } else {
        success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize,
                                     spacing, threshold, sources, data->stamps, type, size, order,
                                     widths, orders, inner, ringsOrder, binning, penalty, optimum,
                                     optWidths, optOrder, optThresh, iter, rej, sys, maskVal,
                                     maskBad, maskPoor, poorFrac, badFrac, subMode);
    }

    psFree(optWidths);
    pmSubtractionThreadsFinalize(inRO, refRO);

    if (!success) {
        psErrorCode error = psErrorCodeLast(); // Error code
        if (error == PM_ERR_STAMPS) {
            psErrorStackPrint(stderr, "Unable to find stamps");
            psWarning("Unable to find stamps --- suspect bad data quality.");
            ppSubDataQuality(data, error, PPSUB_FILES_ALL);
            return true;
        } else {
            psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
            return false;
        }
    }

    pmConceptsCopyFPA(inConv->parent->parent->parent, inRO->parent->parent->parent, true, true);
    pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true);

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

    return true;
}
