# include "psphotInternal.h"

// In this function, we smooth the image and variance, then generate the significance image :
// (S/N)^2.  If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the
// smoothing kernel.
psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal) {

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

    // smooth the image and variance map
    psTimerStart ("psphot.smooth");
    bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve

    // XXX we can a) choose fft to convolve if needed and b) multithread fftw

    float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
    if (!isfinite(minGauss)) {
        psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
        minGauss = 0.5;
    }

    // if we have already determined the PSF model, then we have a better idea how to smooth this image
    bool statusMajor, statusMinor;
    float fwhmMajor = psMetadataLookupF32(&statusMajor, readout->analysis, "FWHM_MAJ");
    float fwhmMinor = psMetadataLookupF32(&statusMinor, readout->analysis, "FWHM_MIN");
    if (statusMajor && statusMinor) {
        // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?)
        if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) {
            psWarning("fwhmMajor (%f) or fwhmMinor (%f) is bad!", fwhmMajor, fwhmMinor);
        }
        SIGMA_SMTH  = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrt(2.0*log(2.0)));
        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
        guess = false;
    } else {
        // if we do not know the FWHM, use the guess smoothing kernel supplied.
        // it is a configuration error if these are not supplied
        SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
        PS_ASSERT (status, NULL);
        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
        PS_ASSERT (status, NULL);
        guess = true;
    }
    // record the actual smoothing sigma
    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections", SIGMA_SMTH);

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

    // Smooth the variance, applying the mask as we go.  The variance is smoothed by the PSF^2,
    // renomalized to maintain the input level of the variance.  We achieve this by smoothing
    // with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with unity normalization.  Note that
    // this process yields a smoothed image with correlated errors.  The pixel-to-pixel
    // variations in smooth_im will be decreased by a factor of 4*pi*SIGMA_SMTH^2, but for
    // measurements based on apertures comparable to or larger than the smoothing kernel, the
    // effective per-pixel variance is maintained.
    psImage *smooth_wt = psImageCopy(NULL, readout->variance, PS_TYPE_F32);
    psImageSmoothMask_Threaded(smooth_wt, smooth_wt, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss);
    psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth"));

    psImage *mask = readout->mask;

    // optionally save example images under trace
    // XXX change these to recipe value checks
    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);
    }

    // Calculate correction factor for the covariance produced by the (potentially multiple) smoothing
    psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing
    float factor = 1.0 / psImageCovarianceCalculateFactor(kernel, readout->covariance);
    psFree(kernel);

    // record the effective area and significance scaling factor
    float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH);
    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);

    // XXX multithread this if needed
    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.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal)) {
                smooth_im->data.F32[j][i] = 0.0;
            } else {
                smooth_im->data.F32[j][i] = factor * PS_SQR(value) / smooth_wt->data.F32[j][i];
            }
        }
    }
    psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark("psphot.smooth"));

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

    psFree(smooth_wt);

    psImageConvolveSetThreads(oldThreads);

    return smooth_im;
}

# if (0)
{
    // threadingdemo
    // smooth the image, applying the mask as we go
    psImage *smooth_im = psImageCopy(NULL, readout->image, PS_TYPE_F32);
    psImageSmoothMask(smooth_im, smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH, minGauss);

    psLogMsg("psphot", PS_LOG_INFO, "smooth image: %f sec\n", psTimerMark("psphot.smooth"));

    psImageConvolveSetThreads(true);

    psTimerStart ("psphot.smooth");

    // XXX a quick test of the threaded version of the function:
    psImage *smooth_test = psImageCopy(NULL, readout->image, PS_TYPE_F32);
    psImageSmoothMask_Threaded(smooth_test, smooth_test, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH, minGauss);

    psLogMsg("psphot", PS_LOG_INFO, "smooth_threaded image: %f sec\n", psTimerMark("psphot.smooth"));

    psphotSaveImage (NULL, smooth_im, "smooth_im.fits");
    psphotSaveImage (NULL, smooth_test, "smooth_test.fits");
    exit (0);
}
# endif
