# include "psphotInternal.h"

// In this function, we smooth the image, then search for the peaks
psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe,
                          bool returnFootprints,
                          const int pass, psMaskType maskVal) {

    float SIGMA_SMTH, NSIGMA_SMTH, NSIGMA_PEAK;
    bool  status = false;

    // smooth the image and weight map
    psTimerStart ("psphot");

    if (pass == 1) {
        SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    } else {
        bool status_x, status_y;
        float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
        float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
        if (!status_x | !status_y) {
            psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined");
            return false;
        }
        SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    }

    // smooth the image, applying the mask as we go
    psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
    psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH);
    psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot"));

    // smooth the weight, applying the mask as we go
    psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
    psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
    psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot"));

    psImage *mask = readout->mask;

    // optionally save example images under trace
    if (psTraceGetLevel("psphot") > 5) {
        char name[64];
        sprintf (name, "imsmooth.v%d.fits", pass);
        psphotSaveImage (NULL, smooth_im, name);
        sprintf (name, "wtsmooth.v%d.fits", pass);
        psphotSaveImage (NULL, smooth_wt, name);
    }

    // build the significance image on top of smooth_im
    for (int j = 0; j < smooth_im->numRows; j++) {
        for (int i = 0; i < smooth_im->numCols; i++) {
            float value = smooth_im->data.F32[j][i];
            if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.U8[j][i] & maskVal)) {
                smooth_im->data.F32[j][i] = 0.0;
            } else {
                smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i];
            }
        }
    }
    psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot"));

    // optionally save example images under trace
    if (psTraceGetLevel("psphot") > 5) {
        char name[64];
        sprintf (name, "snsmooth.v%d.fits", pass);
        psphotSaveImage (NULL, smooth_im, name);
    }

    psTimerStart ("psphot");
    // set peak threshold

    // signal/noise limit for the detected peaks
    if (pass == 1) {
        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
    } else {
        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
    }

    // we need to define the threshold based on the value of NSIGMA_PEAK and the applied smoothing
    // gaussian SIGMA.  a peak in the significance image has an effective S/N for faint sources
    // of (S = Io*2pi*sigma_eff^2) / sqrt(Var), where sigma_eff^2 = sigma_obs^2 + sigma_sm^2.
    // the smoothing of the varience map does not affect the faint-source S/N since the noise
    // is constant under a faint source.

    // if sigma_sm = sigma_obs, then sigma_eff^2 = 2 sigma_sm^2. in this case, the threshold in
    // terms of smooth_im peak counts Io, for a desired S/N limit corresponds to
    // S/N = sqrt(Io)*4*pi*sigma_sm^2
    // thus, the threshold is:
    float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH);
    if (effArea < 1) {
        effArea = 1;                    // never less than a pixel
    }
    float threshold = PS_SQR(NSIGMA_PEAK) / effArea;

    // find the peaks in the smoothed image
    psArray *peaks = pmFindImagePeaks (smooth_im, threshold);
    if (peaks == NULL) {
        // XXX this may also be due to a programming or config error
        // XXX do we need to set something in the readout->analysis to indicate that
        // we tried and failed to find peaks (something in the header data)
        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
        return false;
    }

    // correct the peak values to S/N = sqrt(value*effArea)
    // get the peak flux from the unsmoothed image
    // the peak pixel coords are guaranteed to be on the image
    int row0 = readout->image->row0;
    int col0 = readout->image->col0;
    for (int i = 0; i < peaks->n; i++) {
        pmPeak *peak = peaks->data[i];
        peak->SN = sqrt(peak->value*effArea);
        peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
    }

    // limit the total number of returned peak as specified
    if (pass == 1) {
        psArraySort (peaks, pmPeakSortBySN);
        int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX");
        if (NMAX && (peaks->n > NMAX)) {
            psArray *tmpPeaks = psArrayAllocEmpty (NMAX);
            for (int i = 0; i < NMAX; i++) {
                psArrayAdd (tmpPeaks, 100, peaks->data[i]);
            }
            psFree (peaks);
            peaks = tmpPeaks;
        }
    }

    // optional dump of all peak data
    char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE");
    if (status && (output != NULL) && (output[0])) {
        pmPeaksWriteText (peaks, output);
    }
    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot"));

    // If they asked us to return a set of pmFootprints, not a raw pmPeak list,
    // go ahead and find the footprints, and assign them their peaks.
    //
    // N.b. We're not culling this list; call pmFootprintCullPeaks if you
    // want to do that
    //
    if (returnFootprints) {     // We want an array of pmFootprint, not pmPeak
        int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN");
        if (!status) {
            npixMin = 1;
        }
        float FOOTPRINT_NSIGMA_LIMIT =
            psMetadataLookupS32(&status, recipe,
                                (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2");
        if (!status) {
            FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK;
        }
        threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea;

        psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin);
        pmPeaksAssignToFootprints(footprints, peaks);

        psFree(peaks);
        peaks = footprints;             // well, you know what I mean
    }

    psFree (smooth_im);
    psFree (smooth_wt);

    return (peaks);
}


/************************************************************************************************************/
/*
 * Cull a set of peaks contained in a psArray of pmFootprints
 */
psErrorCode
psphotCullPeaks(const psImage *image,   // the image wherein lives the footprint
                const psImage *weight,  // corresponding variance image
                const psMetadata *recipe,
                psArray *footprints) {  // array of pmFootprints
    bool status = false;
    float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA");
    if (!status) {
        nsigma_delta = 0; // min.
    }
    float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN");
    if (!status) {
        nsigma_min = 0;
    }
    const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV");

    return pmFootprintArrayCullPeaks(image, weight, footprints,
                                     nsigma_delta, nsigma_min*skyStdev);
}
