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

# include "psphot.h"

// 2006.02.02 : no leaks
psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe) {

    bool  status = false;
    float NSIGMA;
    float SIGMA;
    float threshold;
    float value;

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

    SIGMA  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
    NSIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");

    // XXX we are somewhat doing the wrong thing here with the mask.
    // it would be better to smooth the image using the mask information
    psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
    psImageSmooth (smooth_im, SIGMA, NSIGMA);
    psLogMsg ("psphot", 4, "smooth image: %f sec\n", psTimerMark ("psphot"));

    psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
    psImageSmooth (smooth_wt, SIGMA, NSIGMA);
    psLogMsg ("psphot", 4, "smooth weight: %f sec\n", psTimerMark ("psphot"));

    psImage *mask = readout->mask;

    // XXX make this a user-option?
    // psphotSaveImage (NULL, smooth_im, "imsmooth.fits");
    // psphotSaveImage (NULL, smooth_wt, "wtsmooth.fits");

    psTimerStart ("psphot");
    for (int j = 0; j < smooth_im->numRows; j++) {
	for (int i = 0; i < smooth_im->numCols; i++) {
	    value = smooth_im->data.F32[j][i];
	    smooth_im->data.F32[j][i] = 0.0;
	    if (mask->data.U8[j][i]) continue;
	    if (smooth_wt->data.F32[j][i] <= 0) continue;
	    smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i];
	}
    }
    psLogMsg ("psphot", 4, "built S/N map: %f sec\n", psTimerMark ("psphot"));

    psTimerStart ("psphot");

    // set peak threshold
    NSIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
    // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; 
    threshold = PS_SQR(NSIGMA);
    psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold);

    // find the peaks in the smoothed image 
    psArray *peaks = pmFindImagePeaks (smooth_im, threshold);
    if (peaks == NULL) psAbort ("find peaks", "no peaks found");
    psFree (smooth_im);
    psFree (smooth_wt);

    // 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", 3, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot"));

    return (peaks);
}

// In this function, we smooth the image, then search for the peaks 
// Should we also subtract a super-binned image? (as an option?)
// XXX : We need to gracefully handle no source detections
