/** @file  pmSourceSky.c
 *
 *  Functions to measure the local sky and sky variance for sources on images
 *
 *  @author GLG, MHPCC
 *  @author EAM, IfA: significant modifications.
 *
 *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2008-07-15 20:25:00 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

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

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAMaskWeight.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmTrend2D.h"
#include "pmPSF.h"
#include "pmModel.h"
#include "pmSource.h"
#include "pmSourceSky.h"

/******************************************************************************
pmSource *pmSourceLocalSky(source, statsOptions, Radius): this
routine creates a new pmSource.moments element if needed and sets pmSource.pmMoments.sky

The sky value is set from the pixels in the square annulus surrounding the
peak pixel.

The source.pixels and source.mask must already exist

This function modifies the source mask; it should only be called before the object aperture is defined
*****************************************************************************/

bool pmSourceLocalSky(
    pmSource *source,
    psStatsOptions statsOptions,
    psF32 Radius,
    psMaskType maskVal,
    psMaskType markVal)
{
    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(source, false);
    PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);
    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    PS_ASSERT_INT_POSITIVE(Radius, false);

    psStatsOptions statistic = psStatsSingleOption(statsOptions);
    if (statistic == 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics specified: %x\n", statsOptions);
        return NULL;
    }

    psImage *image = source->pixels;
    psImage *mask  = source->maskObj;
    pmPeak *peak  = source->peak;
    psRegion srcRegion;

    // maskVal is used to test for rejected pixels, and must include markVal
    maskVal |= markVal;

    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
    srcRegion = psRegionForImage(mask, srcRegion);

    psImageMaskRegion(mask, srcRegion, "OR", markVal);
    psStats *myStats = psStatsAlloc(statsOptions);
    if (!psImageStats(myStats, image, mask, maskVal)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
        psFree(myStats);
        return false;
    }
    psImageMaskRegion(mask, srcRegion, "AND", PS_NOT_U8(markVal));
    double value = psStatsGetValue(myStats, statistic);
    psFree(myStats);

    if (isnan(value)) {
        psTrace("psModules.objects", 3, "---- %s(false) end ----\n", __func__);
        return(false);
    }
    if (source->moments == NULL) {
        source->moments = pmMomentsAlloc();
    }
    source->moments->Sky = value;
    psTrace("psModules.objects", 3, "---- %s(true) end ----\n", __func__);
    return (true);
}

// A complementary function to pmSourceLocalSky: calculate the local median variance
bool pmSourceLocalSkyVariance(
    pmSource *source,
    psStatsOptions statsOptions,
    psF32 Radius,
    psMaskType maskVal,
    psMaskType markVal
)
{
    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(source, false);
    PS_ASSERT_IMAGE_NON_NULL(source->weight, false);
    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    PS_ASSERT_INT_POSITIVE(Radius, false);

    // maskVal is used to test for rejected pixels, and must include markVal
    maskVal |= markVal;

    psStatsOptions statistic = psStatsSingleOption(statsOptions);
    if (statistic == 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics specified: %x\n", statsOptions);
        return NULL;
    }

    psImage *image = source->weight;
    psImage *mask  = source->maskObj;
    pmPeak *peak  = source->peak;
    psRegion srcRegion;

    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
    srcRegion = psRegionForImage(mask, srcRegion);

    psImageMaskRegion(mask, srcRegion, "OR", markVal);
    psStats *myStats = psStatsAlloc(statsOptions);
    if (!psImageStats(myStats, image, mask, maskVal)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
        psFree(myStats);
        return false;
    }
    psImageMaskRegion(mask, srcRegion, "AND", PS_NOT_U8(markVal));
    double value = psStatsGetValue(myStats, statistic);
    psFree(myStats);

    if (isnan(value)) {
        psTrace("psModules.objects", 3, "---- %s(false) end ----\n", __func__);
        return(false);
    }
    if (source->moments == NULL) {
        source->moments = pmMomentsAlloc();
    }
    source->moments->dSky = value;
    psTrace("psModules.objects", 3, "---- %s(true) end ----\n", __func__);
    return (true);
}
