# include "psphotInternal.h"

bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view, bool firstPass)
{
    bool status = true;

    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");

    // loop over the available readouts
    for (int i = 0; i < num; i++) {
	if (!psphotStackMatchPSFsReadout (config, view, i)) {
            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
	    return false;
	}
    }
    return true;
}

// convolve the image to match desired PSF
bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) {

    bool status;
    int pass;
    float NSIGMA_PEAK = 25.0;
    int NMAX = 0;

    // find the currently selected readout
    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
    psAssert (fileRaw, "missing file?");

    pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CONV", index); // File of interest
    psAssert (fileCnv, "missing file?");

    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
    psAssert (readoutRaw, "missing readout?");

    pmReadout *readoutCnv = pmFPAviewThisReadout(view, fileCnv->fpa);
    psAssert (readoutCnv, "missing readout?");

    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    psAssert (maskVal, "missing mask value?");

    /***** set up recipe options *****/

    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    psAssert(stackRecipe, "We've thrown an error on this before.");

    // Look up appropriate values from the ppSub recipe
    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    psAssert(subRecipe, "recipe missing");

    int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size

    float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction

    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor
    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels

    bool mdok;                          // Status of MD lookup
    float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads

    if (!pmReadoutMaskNonfinite(readout, maskVal)) {
        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout.");
        return false;
    }

    // Image Matching (PSFs or just flux)
    if (options->convolve) {
      // Full match of PSFs
        pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily

        // Read previously produced kernel
        if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
	  loadKernel();
        } else {
	  matchKernel();
        } // !DEBUG.STACK

	saveMatchData();

	saveChiSquare();

	renormKernel();

        // Reject image completely if the maximum deconvolution fraction exceeds the limit
        float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
        if (deconv > deconvLimit) {
            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
            psFree(conv);
            return NULL;
        }

        readout->analysis = psMetadataCopy(readout->analysis, conv->analysis);

        psFree(conv);
    } else {
        // only match the flux
        float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
        psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
        psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
    }

    // Ensure the background value is zero
    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
        psWarning("Can't measure background for image.");
        psErrorClear();
    } else {
        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
            (void)psBinaryOp(readout->image, readout->image, "-",
                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
        }
    }

    if (!stackRenormaliseReadout(config, readout)) {
        psFree(rng);
        psFree(bg);
        return false;
    }

    // Measure the variance level for the weighting
    if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
            psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
            psFree(rng);
            psFree(bg);
            return false;
        }
        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
    } else {
        options->weightings->data.F32[index] = 1.0;
    }
    psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
             index, options->weightings->data.F32[index]);

    psFree(rng);
    psFree(bg);

    dumpImage3();

    return true;
}
