/** @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.1.2.3 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-04-19 06:31:55 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "pslib.h"
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmMaskBadPixels.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(__func__, 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);

    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_EXT2));


    psF64 tmpF64;
    p_psGetStatValue(myStats, &tmpF64);
    psFree(myStats);

    if (isnan(tmpF64)) {
        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
        return(false);
    }
    if (source->moments == NULL) {
        source->moments = pmMomentsAlloc();
    }
    source->moments->Sky = (psF32) tmpF64;
    psTrace(__func__, 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(__func__, 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);

    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));

    psF64 tmpF64;
    p_psGetStatValue(myStats, &tmpF64);
    psFree(myStats);

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