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

#include "ppImage.h"

bool ppImageDetrendReadout(pmConfig *config, ppImageOptions *options, pmFPAview *view)
{
  // psTimerStart("detrend.readout");

    // 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);
	}
      }
    }
    // psLogMsg ("ppImage", 6, "check gain: %f sec\n", psTimerMark ("detrend.readout"));

    // Check to see if we're in a chip that contains video
    // NOTE: this is an esoteric test for a GPC1 feature : this is irrelevant & will not work for other cameras
    bool hasVideo = false;
    {
      // XXX test for GPC1? or CAN_HAVE_VIDEO in camera config?
      if (!input) goto done_video_check;
      if (!input->parent) goto done_video_check;
      if (!input->parent->parent) goto done_video_check;
      if (!input->parent->parent->hdu) goto done_video_check;
      if (!input->parent->parent->hdu->header) goto done_video_check;
      if (psMetadataLookupStr(NULL,input->parent->parent->hdu->header,"CELLMODE")) {
	psLogMsg ("ppImage.detrend", PS_LOG_DETAIL, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask);
	char *Vptr = strchr(psMetadataLookupStr(NULL,input->parent->parent->hdu->header,"CELLMODE"),'V');
	if (Vptr) {
	  hasVideo = options->hasVideo = true;
	  psLogMsg ("ppImage.detrend", PS_LOG_INFO, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask);
	}
      }
    }
    done_video_check:
    // psLogMsg ("ppImage", 6, "check video: %f sec\n", psTimerMark ("detrend.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;
        if (!pmReadoutGenerateMask(input, satMask, lowMask)) {
	  psError(PS_ERR_UNKNOWN, false, "Unable to generate a mask.");
	  psFree(detview);
	  return false;
	}
	// psLogMsg ("ppImage", 6, "generate mask: %f sec\n", psTimerMark ("detrend.readout"));
    }
    // apply the externally supplied mask to the input->mask pixels
    if (options->doMask) {
      pmReadout *mask;
      if ((options->useVideoMask)&&(hasVideo)) {
        mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEOMASK");
      }
      else {
	mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK");
      }
      pmMaskBadPixels(input, mask, options->maskValue);
      // psLogMsg ("ppImage", 6, "apply mask: %f sec\n", psTimerMark ("detrend.readout"));
    }

    if (options->doApplyBurntool)
    {
        // extern bool ppImageBurntoolApply(pmConfig *, ppImageOptions *, pmFPAview *, pmReadout *);
        ppImageBurntoolApply(config, options, view, input);
	// psLogMsg ("ppImage", 6, "apply burntool: %f sec\n", psTimerMark ("detrend.readout"));
    }
    
    if (options->doMaskBurntool) {
        if (options->doApplyBurntool) {
            // build burntool mask from data input burntool table
            ppImageBurntoolMaskFromTable(config,options,view,input);
        } else {
            // build burntool mask from data in the input image's fits extension
            ppImageBurntoolMask(config,options,view,input);
        }
	// psLogMsg ("ppImage", 6, "apply burntool mask: %f sec\n", psTimerMark ("detrend.readout"));
    }

    // Subtract the overscan
    if (options->doOverscan) {
      if (!pmOverscanSubtract (input, options->overscan)) {
	psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan.");
	psFree(detview);
	return false;
      }
      // psLogMsg ("ppImage", 6, "subtract overscan: %f sec\n", psTimerMark ("detrend.readout"));
    }

    // Non-linearity correction
    if (options->doNonLin) {
      if (!ppImageDetrendNonLinear(input,detview,config)) {
	psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity");
	psFree(detview);
	return(false);
      }
      // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
    }

    // 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 ((options->useVideoDark)&&(hasVideo)) {
	  if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) {
            oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEODARK");
	  } else {
            dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.VIDEODARK");
	  }
	}
	else {
	  if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) {
            oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK");
	  } else {
            dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK");
	  }
	}
    }

    // Bias and temperature-independent-dark subtraction are merged.
    if (options->doBias) {
        if (!pmBiasSubtract(input, bias, oldDark, view)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to subtract bias.");
            psFree(detview);
            return false;
        }
	// psLogMsg ("ppImage", 6, "apply bias: %f sec\n", psTimerMark ("detrend.readout"));
    }
    
    // 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.
	    // CZW 2012-03-21: I do not believe this is true anymore.  In any case, this seems to be what
	    //                 ppImageReplaceBackground does to do sky subtraction.
            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);
	    // Stolen from pmSkySubtract.c
	    // CZW: Unneeded, as psImageUnbin does the bilinear interpolation?  It still looks blocky. 
/* 	    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe */
/* 	    psAssert(recipe, "Should be there!"); */

/* 	    psS32 xBinFactor = psMetadataLookupS32(NULL,recipe,"NOISE.XBIN"); */
/* 	    psS32 yBinFactor = psMetadataLookupS32(NULL,recipe,"NOISE.YBIN"); */
/* 	    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, */
/* 									       noiseMap->image, */
/* 									       NULL, NULL, */
/* 									       0, 0.0, 0.0, 0, 0, 0.0); */
/* 	    for (psS32 row = 0; row < input->image->numRows; row++) { */
/* 	      for (psS32 col = 0; col < input->image->numCols; col++) { */
/* 		// We calculate the F32 value of the pixel coordinates in the */
/* 		// binned image and then use a pixel interpolation routine to */
/* 		// determine the value of the pixel at that location. */
/* 		psF32 binRowF64 = ((psF32) row) / ((psF32) yBinFactor); */
/* 		psF32 binColF64 = ((psF32) col) / ((psF32) xBinFactor); */
/* 		// We add 0.5 to the pixel locations since the pixel */
/* 		// interpolation routine defines the location of pixel */
/* 		// (i, j) as (i+0.5, j+0.5). */
/* 		binRowF64+= 0.5; */
/* 		binColF64+= 0.5; */

/* 		double binPixel; */
/* 		if (!psImagePixelInterpolate(&binPixel, NULL, NULL, binColF64, binRowF64, interp)) { */
/* 		  psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); */
/* 		  psFree(interp); */
/* 		  psFree(noiseImage); */
/* 		  return NULL; */
/* 		} */
/* 		noiseImage->data.F32[row][col] = binPixel; */
/* 	      } */
/* 	    } */
/* 	    psFree(interp); */
        }

        if (!pmReadoutGenerateVariance(input, noiseImage, true)) {
	    psError(PS_ERR_UNKNOWN, false, "Unable to generate a variance image.");
	    psFree(detview);
	    return false;
	}
        psFree (noiseImage);
	// psLogMsg ("ppImage", 6, "generate variance: %f sec\n", psTimerMark ("detrend.readout"));
    }

    if (options->doDark && dark) {
        if (!pmDarkApply(input, dark, options->darkMask)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to subtract dark.");
            psFree(detview);
            return false;
        }
	// psLogMsg ("ppImage", 6, "apply dark: %f sec\n", psTimerMark ("detrend.readout"));
    }

    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;
        }
	// psLogMsg ("ppImage", 6, "mask remnance: %f sec\n", psTimerMark ("detrend.readout"));
    }

    // Shutter correction
    if (options->doShutter) {
        pmReadout *shutter = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.SHUTTER");
        if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config))) {
            psFree(detview);
            return false;
        }
	// psLogMsg ("ppImage", 6, "shutter correction: %f sec\n", psTimerMark ("detrend.readout"));
    }

    // 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;
        }
	// psLogMsg ("ppImage", 6, "apply flat: %f sec\n", psTimerMark ("detrend.readout"));
    }

/*     // 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));
	// psLogMsg ("ppImage", 6, "renormalize: %f sec\n", psTimerMark ("detrend.readout"));
    }

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

	// psLogMsg ("ppImage", 6, "renormalize by class: %f sec\n", psTimerMark ("detrend.readout"));
    }
# endif

    if (options->doFringe) {
        pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE");
        if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) {
            psFree(detview);
            return false;
        }
	// psLogMsg ("ppImage", 6, "measure fringe: %f sec\n", psTimerMark ("detrend.readout"));
    }

    ppImageMemoryDump("detrend");

    // psLogMsg ("ppImage", 5, "detrend readout: %f sec\n", psTimerMark ("detrend.readout"));
    psFree(detview);
    return true;
}
