# include "psphotInternal.h"

// Find peaks in the significance image above a threshold significance level.  The significance
// image must be constructed to represent (S/N)^2.  If nMax is non-zero, only return a maximum
// of nMax peaks
psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks) {

    bool status = false;

    psTimerStart ("psphot.peaks");

    // find the peaks in the smoothed image
    // NOTE : significance->variance actually carries the detection S/N image
    psArray *peaks = pmPeaksInImage (significance->variance, threshold);
    if (peaks == NULL) {
	// we only get a NULL peaks array due to a programming or config error. 
	// this will result in a failure.
        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
        return NULL;
    }
    // return the total number of peaks found before the nMax limit is applied
    if (totalPeaks) {
        *totalPeaks = peaks->n;
    }

    if (peaks->n == 0) {
        // 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)
	psLogMsg ("psphot", PS_LOG_INFO, "no peaks found in this image");
	return peaks;
    }

    // Convert the peak values to S/N = sqrt(significance).
    // Get the peak flux from the unsmoothed image.
    // Rescale the peak position errors using the peak variance
    // 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->rawFlux = readout->image->data.F32[peak->y-row0][peak->x-col0];
	peak->rawFluxStdev = sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
        peak->smoothFlux = significance->image->data.F32[peak->y-row0][peak->x-col0];
	peak->smoothFluxStdev = peak->smoothFlux / sqrt(significance->variance->data.F32[peak->y-row0][peak->x-col0]);
	// NOTE smoothFluxStdev is actually (sqrt(variance) / covar_factor)

	// do we need this or not?
        // peak->SN = sqrt(peak->detValue);

	if (readout->variance && isfinite (peak->dx)) {
	    peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
	}
	if (readout->variance && isfinite (peak->dy)) {
	    peak->dy *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
	}
    }

    // limit the total number of returned peaks as specified
    psArraySort (peaks, pmPeaksSortByRawFluxDescend);
    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 (output && strcasecmp (output, "NONE")) {
        pmPeaksWriteText (peaks, output);
    }
    psLogMsg ("psphot", PS_LOG_WARN, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot.peaks"));

    return peaks;

}
