/** @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.8 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-09-15 09:49:01 $
 *
 *  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 "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
*****************************************************************************/

bool pmSourceLocalSky(
    pmSource *source,
    psStatsOptions statsOptions,
    psF32 Radius)
{
    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->mask, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    PS_ASSERT_INT_POSITIVE(Radius, false);
    PS_ASSERT_INT_NONNEGATIVE(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->mask;
    pmPeak *peak  = source->peak;
    psRegion srcRegion;

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

    psImageMaskRegion(mask, srcRegion, "OR", PM_MASK_MARK);
    psStats *myStats = psStatsAlloc(statsOptions);
    myStats = psImageStats(myStats, image, mask, 0xff);
    psImageMaskRegion(mask, srcRegion, "AND", NOT_U8(PM_MASK_MARK));
    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)
{
    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->mask, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    PS_ASSERT_INT_POSITIVE(Radius, false);
    PS_ASSERT_INT_NONNEGATIVE(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->weight;
    psImage *mask  = source->mask;
    pmPeak *peak  = source->peak;
    psRegion srcRegion;

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

    psImageMaskRegion(mask, srcRegion, "OR", PM_MASK_MARK);
    psStats *myStats = psStatsAlloc(statsOptions);
    myStats = psImageStats(myStats, image, mask, 0xff);
    psImageMaskRegion(mask, srcRegion, "AND", NOT_U8(PM_MASK_MARK));
    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);
}
