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

#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>
#include <psphot.h>

#include "ppStack.h"

#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
                          PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
                          PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources

bool ppStackInputPhotometry(const pmReadout *ro, const psArray *sources, const pmConfig *config)
{
    PM_ASSERT_READOUT_NON_NULL(ro, false);
    PS_ASSERT_ARRAY_NON_NULL(sources, false);

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

    bool mdok;                          // Status of MD lookup

    float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
    if (!mdok) {
        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
        return false;
    }
    float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
    if (!mdok) {
        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
        return false;
    }
    float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
    if (!mdok) {
        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
        return false;
    }

    psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
    if (!mdok || !maskValStr) {
        psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
        return false;
    }
    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask

    psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout

    // Measure background
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
    if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
        psFree(stats);
        psFree(rng);
        return false;
    }
    float bg = stats->robustMedian; // Background level
    psFree(stats);
    psFree(rng);

    // Photometer sources
    int numCols = image->numCols, numRows = image->numRows; // Size of image
    int numSources = sources->n; // Number of sources
    int numGood = 0;            // Number of good sources
    float maxMag = -INFINITY;   // Maximum magnitude
    for (int j = 0; j < numSources; j++) {
        pmSource *source = sources->data[j]; // Source of interest
        if (source->mode & PHOT_SOURCE_MASK || !isfinite(source->psfMag)) {
            source->psfMag = NAN;
            continue;
        }
        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
        int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
        // Range of integration
        int xMin = PS_MAX(0, xCentral - zpRadius), xMax = PS_MIN(numCols - 1, xCentral + zpRadius);
        int yMin = PS_MAX(0, yCentral - zpRadius), yMax = PS_MIN(numRows - 1, yCentral + zpRadius);

        // Integrate
        double sumImage = 0.0, sumKernel = 0.0; // Sums from integration
        int numBadPix = 0;         // Number of bad pixels
        for (int v = yMin; v <= yMax; v++) {
            float dy2 = PS_SQR(y - v); // Distance from centroid
            for (int u = xMin; u <= xMax; u++) {
                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[v][u] & maskVal) {
                    numBadPix++;
                    continue;
                }
                float dx2 = PS_SQR(x - u); // Distance from centroid
                double kernel = exp(-0.5 * (dx2 + dy2) / PS_SQR(zpSigma)); // Kernel value
                sumImage += (image->data.F32[v][u] - bg) * kernel;
                sumKernel += kernel;
            }
        }

        if (numBadPix > zpFrac * M_PI * PS_SQR(zpRadius) || !isfinite(sumImage)) {
            source->psfMag = NAN;
        } else {
            source->psfMag = - 2.5 * log10(sumImage / sumKernel * M_PI * PS_SQR(zpRadius));
            if (isfinite(source->psfMag)) {
                numGood++;
                maxMag = PS_MAX(maxMag, source->psfMag);
            }
        }
    }
    psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);

    return true;
}



bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view)
{
    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    pmFPACopy(photFile->fpa, readout->parent->parent->parent);

    if (psMetadataLookupBool(NULL, config->arguments, "-visual")) {
        pmVisualSetVisual(true);
    }

    // Need to ensure aperture residual is not calculated
    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
    psMetadataItem *item = psMetadataLookup(recipe, "MEASURE.APTREND"); // Item determining aptrend
    if (!item) {
        psWarning("Unable to find MEASURE.APTREND in psphot recipe");
        psErrorClear();
    } else {
        item->data.B = false;
    }

    // set maskValue and markValue in the psphot recipe
    psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for mark", markValue);

    // XXX replace with psphotReadoutKnownSources (config, view)
    // XXX this function requires that we construct the input source list and place it on PSPHOT.INPUT.CMF
    pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
    psArray *inSources = psMetadataLookupPtr (NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
    if (!inSources) {
        psError(PS_ERR_UNKNOWN, false, "Unable to find input sources");
        return false;
    }

    if (!psphotReadoutKnownSources(config, view, inSources)) {
        // Clear the error, so that the output files are written.
        psWarning("Unable to perform photometry on stacked image.");
        psErrorStackPrint(stderr, "Error stack from photometry:");
        psErrorClear();
    }

    if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
        !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
        !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
        psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
        return false;
    }

    pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");

    return true;
}
