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

    // Check that the gain is set (this is used by both pmReadoutGenerateMask and pmReadoutGenerateVariance)
    { 
      float gain = psMetadataLookupF32(NULL, input->parent->concepts, "CELL.GAIN"); // Gain for cell
      if (!isfinite(gain)) {
	psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe
	psAssert(recipe, "Should be there!");
	bool override = psMetadataLookupBool(NULL, recipe, "GAIN.OVERRIDE"); // Override the bad gain?
	if (override) {
	  psWarning("CELL.GAIN is not set for readout (%d,%d,%d) --- setting to unity.", view->chip, view->cell, view->readout);
	  psMetadataItem *item = psMetadataLookup(input->parent->concepts, "CELL.GAIN"); // Gain item
	  psAssert(item, "Should be there!");
	  item->data.F32 = 1.0;

	  // for unity gain, there is no modification for the readnoise, note that it has (effectively) been updated
	  psMetadataRemoveKey(input->parent->concepts, "CELL.READNOISE.UPDATE");
	} else {
	  psWarning("CELL.GAIN is NAN for readout (%d,%d,%d), image will be masked.", view->chip, view->cell, view->readout);
	}
      }
    }

    // Masking on the basis of pixel value needs to be done before anything else, so the values are pristine.
    if (options->doMaskBuild) {
        psImageMaskType satMask = options->doMaskSat ? options->satMask : 0;
        psImageMaskType lowMask = options->doMaskLow ? options->lowMask : 0;
        pmReadoutGenerateMask(input, satMask, lowMask);
    }
    // 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 (options->doMaskBurntool) {
      ppImageBurntoolMask(config,options,view,input);
    }


# 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.");
            psFree(detview);
            return false;
        }
    }

    // Weight on the basis of pixel value needs to be done after the overscan has been subtracted
    if (options->doVarianceBuild) {
        // create the target mask and variance images
        psImage *noiseImage = NULL;
        if (options->doNoiseMap) {
            // XXX convert the noiseMap image to a binned image
            pmReadout *noiseMap = NULL;
            noiseMap = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.NOISEMAP");
            noiseImage = psImageCopy (NULL, input->image, PS_TYPE_F32);
            psImageInit (noiseImage, 0.0);

            // XXX this works, but is not really quite right: the model shoud include the
            // offset information, we are not really getting exactly the right mapping from the
            // original file.
            psImageBinning *binning = psImageBinningAlloc();
            binning->nXruff = noiseMap->image->numCols;
            binning->nYruff = noiseMap->image->numRows;
            binning->nXfine = input->image->numCols;
            binning->nYfine = input->image->numRows;
            psImageBinningSetScale(binning, PS_IMAGE_BINNING_LEFT);

            psImageUnbin (noiseImage, noiseMap->image, binning);
            psFree (binning);
        }

        pmReadoutGenerateVariance(input, noiseImage, true);
        psFree (noiseImage);
    }

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

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

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

/*     // Pattern noise correction */
/*     if (options->doPattern) { */
/*         if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */
/*                           options->patternThresh, options->patternMean, options->patternStdev, */
/*                           options->maskValue, options->darkMask)) { */
/*             psFree(detview); */
/*             return false; */
/*         } */
/*     } */

    // Normalization by a single (known) constant
    bool mdok;                          // Status of MD lookup
    float norm = psMetadataLookupF32(&mdok, config->arguments, "NORMALIZATION");
    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 (1)
    // Normalization by per-class values
    psMetadata *normlist = psMetadataLookupMetadata(&mdok, config->arguments, "NORMALIZATION.TABLE");
    if (normlist) {
        pmFPAfile *inputFile = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.INPUT");

        // get the menu of class IDs
        psMetadata *menu = psMetadataLookupMetadata(&mdok, inputFile->camera, "CLASSID");
        if (!menu) {
            psError(PS_ERR_IO, false, "Unable to find CLASSID metadata in camera configuration");
            psFree(detview);
            return false;
        }
        // get the rule for class_id for the desired class
        const char *rule = psMetadataLookupStr(&mdok, menu, options->normClass);
        if (!rule) {
            psError(PS_ERR_IO, false, "Unable to find NORM.CLASS value %s in CLASSID in camera configuration", options->normClass);
            psFree(detview);
            return false;
        }
        // get the class_id from the rule
        char *classID = pmFPAfileNameFromRule(rule, inputFile, view);
        if (!classID) {
            psError(PS_ERR_IO, false, "error converting CLASSID rule %s to name\n", rule);
            psFree(detview);
            return false;
        }

        // get normalization from the class_id
        float norm = psMetadataLookupF32 (&mdok, normlist, classID);
        if (!mdok) {
            psError(PS_ERR_IO, false, "failed to find class ID %s in normalization table\n", classID);
            psFree(detview);
            return false;
        }

        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);

        // apply the normalization
        psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32));

        psFree (classID);
    }
# endif

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

    ppImageMemoryDump("detrend");

    psFree(detview);
    return true;
}
