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

#include "ppImage.h"

bool ppImageDetrendReadout(pmConfig *config, ppImageOptions *options, pmFPAview *view)
{
    // construct a view for the detrend images (which have only one readout)
    pmFPAview *detview = pmFPAviewAlloc(0);
    *detview = *view;
    detview->readout = 0;

    // find the currently selected readout
    pmReadout *input = pmFPAfileThisReadout(config->files, view, "PPIMAGE.INPUT");

    // Masking on the basis of pixel value needs to be done before anything else, so the values are pristine.
    if (options->doMaskBuild) {
        pmReadoutGenerateMask(input, options->satMask, options->badMask);
    }
    // apply the externally supplied mask to the input->mask pixels
    if (options->doMask) {
        pmReadout *mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK");
        pmMaskBadPixels(input, mask, options->maskValue);
    }

# if 0
    // Non-linearity correction
    if (options->doNonLin) {
        ppImageDetrendNonLinear(detrend->input, input, options);
    }
# endif

    // set up the dark and bias
    pmCell *dark = NULL;                // Multi-dark
    pmReadout *oldDark = NULL;          // Old-fashioned dark
    pmReadout  *bias = NULL;
    if (options->doBias) {
        bias = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.BIAS");
    }
    if (options->doDark) {
        bool mdok;                      // Status of MD lookup
        psMetadata *recipe = psMetadataLookupPtr (&mdok, config->recipes, RECIPE_NAME);
        assert(mdok && recipe);
        if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) {
            oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK");
        } else {
            dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK");
        }
    }

    // Bias, dark and overscan subtraction are all merged.
    if (options->doBias || options->doOverscan) {
        if (!pmBiasSubtract(input, options->overscan, bias, oldDark, view)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to subtract bias.");
            return false;
        }
    }

    // Weight on the basis of pixel value needs to be done after the overscan has been subtracted
    if (options->doWeightBuild) {
        // create the target mask and weight images
        pmReadoutGenerateWeight(input, true);
    }

    if (options->doDark && dark) {
        if (!pmDarkApply(input, dark, options->maskValue)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to subtract dark.");
            return false;
        }
    }

    if (options->doRemnance) {
        if (!pmRemnance(input, options->maskValue, options->badMask,
                        options->remnanceSize, options->remnanceThresh)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to mask remnance.");
            return false;
        }
    }

    // Shutter correction
    if (options->doShutter) {
        pmReadout *shutter = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.SHUTTER");
        if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config))) {
            return false;
        }
    }

    // Flat-field correction (no options used?)
    if (options->doFlat) {
        pmReadout *flat = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.FLAT");
        if (!pmFlatField(input, flat, options->flatMask)) {
            return false;
        }
    }

    // Normalisation by a (known) constant
    bool mdok;                          // Status of MD lookup
    float norm = psMetadataLookupF32(&mdok, config->arguments, "NORMALISATION");
    if (mdok && isfinite(norm) && norm != 1.0) {
        pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
        psString comment = NULL;        // Comment to add
        psStringAppend(&comment, "Normalization: %f", norm);
        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
        psFree(comment);

        psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    }

    if (options->doFringe) {
        pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE");
        if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) {
            return false;
        }
    }

    ppImageMemoryDump("detrend");

    psFree(detview);
    return true;
}
