# include "psphotInternal.h"

// In this function, we smooth the image and weight, 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 weight 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;
    }

    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) {
      // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?)
      SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (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(recipe, 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 weight, 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->weight, 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 weight: %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);
    }

    // We have an input image with PSF size sigma_r.  We have smoothed it with a kernel of size
    // sigma_s.  The result is an image with PSF size sigma_o: sigma_o^2 = sigma_r^2 +
    // sigma_s^2.  Ideally, we are choosing sigma_s = sigma_r, in which case sigma_o^2 = 2
    // sigma_s^2.  If we do not know the input image PSF size (initial detection stage), then
    // we are assuming that sigma_r = sigma_s.

    // Build the significance image on top of smooth_im.  We need to correct the ratio im/wt by
    // 2 factors: 1) the relationship the peak value and the integrated flux, and 2) the
    // relationship between the per-pixel variance (var_i) and the total variance included in
    // the flux measurement (effective area).  These latter correction comes from: flux = Io *
    // 2\pi\sigma_o^2 and total variance = var_i * 4\pi\sigma_o^2, thus (S/N)^2 = flux^2 / var
    // = (Io/var_i) \pi sigma_o^2

    // in fact, because of the smoothing and the resulting correlated errors, the per-pixel
    // variance is a bit smaller than the value predicted here.  The total variance should be
    // \alpha 4\pi(\sigma_r_^2 + \sigma_s^2) where \alpha approaches 1 for \sigma_s / \sigma_r
    // << 1, and \alpha approaches 0.5 for \sigma_s / \sigma_r >> 1, with a value of 0.7
    // (actually 1/1.5) for \sigma_s = \sigma_r

    // if the smoothing kernel size is known to be the psf size (\sigma_s = \sigma_r), then we
    // can use the scaling factor of 1/1.5; otherwise, we want to be conservative and use \alpha
    // = 0.5.  Thus, for the guess smoothing kernel, (S/N)^2 = 3.0 (Io/var_i) \pi sigma_s^2 for
    // the guess case and (S/N)^2 = 2.0 (Io/var_i) \pi sigma_s^2 for the known case

    // XXX hmm, somehow the simulations are coming up with 4pi, not 3pi (as if \alpha is 1/2.0,
    // not 1/1.5 for sigma_s = sigma_r...

    float factor = guess ? 4.0 * M_PI * PS_SQR(SIGMA_SMTH) : 4.0 * M_PI * PS_SQR(SIGMA_SMTH);
    // record the effective area and significance scaling factor
    float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH);
    psMetadataAddF32(recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
    psMetadataAddF32(recipe, 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
