/** @file  pmSource.c
 *
 *  Functions to define and manipulate sources on images
 *
 *  @author GLG, MHPCC
 *  @author EAM, IfA: significant modifications.
 *
 *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-04-20 19:19:05 $
 *
 *  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"

static void sourceFree(pmSource *tmp)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psFree(tmp->peak);
    psFree(tmp->pixels);
    psFree(tmp->weight);
    psFree(tmp->mask);
    psFree(tmp->moments);
    psFree(tmp->modelPSF);
    psFree(tmp->modelEXT);
    psFree(tmp->blends);
    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
}

/******************************************************************************
pmSourceAlloc(): Allocate the pmSource structure and initialize its members
to NULL.
*****************************************************************************/
pmSource *pmSourceAlloc()
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource));
    tmp->peak = NULL;
    tmp->pixels = NULL;
    tmp->weight = NULL;
    tmp->mask = NULL;
    tmp->moments = NULL;
    tmp->blends = NULL;
    tmp->modelPSF = NULL;
    tmp->modelEXT = NULL;
    tmp->type = PM_SOURCE_TYPE_UNKNOWN;
    tmp->mode = PM_SOURCE_MODE_DEFAULT;
    psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);

    tmp->psfMag = 32.0;
    tmp->extMag = 32.0;
    tmp->errMag = 32.0;
    tmp->apMag  = 32.0;
    tmp->pixWeight = 0.0;

    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmp);
}

// x,y are defined in the parent image coords of readout->image
bool pmSourceDefinePixels(pmSource *mySource,
                          const pmReadout *readout,
                          psF32 x,
                          psF32 y,
                          psF32 Radius)
{
    psRegion srcRegion;

    // Grab a subimage of the original image of size (2 * outerRadius).
    srcRegion = psRegionForSquare (x, y, Radius);
    srcRegion = psRegionForImage (readout->image, srcRegion);

    mySource->pixels = psImageSubset(readout->image, srcRegion);
    mySource->weight = psImageSubset(readout->weight, srcRegion);
    mySource->mask   = psImageSubset(readout->mask,  srcRegion);
    mySource->region = srcRegion;

    return true;
}

bool pmSourceRedefinePixels(pmSource *mySource,
                            const pmReadout *readout,
                            psF32 x,
                            psF32 y,
                            psF32 Radius)
{
    bool extend;
    psRegion newRegion;

    if (Radius == 0)
        return false;

    // check to see if new region is completely contained within old region
    newRegion = psRegionForSquare (x, y, Radius);
    newRegion = psRegionForImage (readout->image, newRegion);

    extend = false;
    extend |= (int)(newRegion.x0) < (int)(mySource->region.x0);
    extend |= (int)(newRegion.x1) > (int)(mySource->region.x1);
    extend |= (int)(newRegion.y0) < (int)(mySource->region.y0);
    extend |= (int)(newRegion.y1) > (int)(mySource->region.y1);

    extend |= (mySource->pixels == NULL);
    extend |= (mySource->weight == NULL);
    extend |= (mySource->mask == NULL);

    // extend = true;
    if (extend) {
        // Grab a new subimage
        psFree (mySource->pixels);
        psFree (mySource->weight);
        psFree (mySource->mask);

        mySource->pixels = psImageSubset(readout->image,  newRegion);
        mySource->weight = psImageSubset(readout->weight, newRegion);
        mySource->mask   = psImageSubset(readout->mask,   newRegion);
        mySource->region = newRegion;
    }

    return extend;
}

/******************************************************************************
    pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 
    sigma-x, sigma-y plane. return 0,0 clump in case of error. 
*****************************************************************************/

// XXX EAM include a S/N cutoff in selecting the sources?
pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);

    psArray *peaks  = NULL;
    pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};
    pmPSFClump psfClump = emptyClump;

    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
    PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);

    psF32 SCALE, SN_MIN;
    psF32 SX_MAX, SY_MAX, AR_MAX, AR_MIN;

    // find the sigmaX, sigmaY clump
    {
        psStats *stats  = NULL;
        psImage *splane = NULL;
        int binX, binY;
        bool status;

        SCALE = psMetadataLookupF32 (&status, metadata, "MOMENTS_SCALE");
        if (!status)
            SCALE = 0.1;
        SN_MIN = psMetadataLookupF32 (&status, metadata, "MOMENTS_SN_MIN");
        if (!status)
            SN_MIN = 30.0;
        SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
        if (!status)
            SX_MAX = 3.0;
        SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
        if (!status)
            SY_MAX = 3.0;
        AR_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_AR_MAX");
        if (!status)
            AR_MAX = 3.0;
        AR_MIN = 1.0 / AR_MAX;

        // construct a sigma-plane image
        // psImageAlloc does zero the data
        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
        for (int i = 0; i < splane->numRows; i++)
        {
            memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
        }

        // place the sources in the sigma-plane image (ignore 0,0 values?)
        for (psS32 i = 0 ; i < sources->n ; i++)
        {
            pmSource *tmpSrc = (pmSource *) sources->data[i];
            if (tmpSrc == NULL) {
                continue;
            }
            if (tmpSrc->moments == NULL) {
                continue;
            }

            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
            if (tmpSrc->moments->Sx < 0.05) {
                continue;
            }
            if (tmpSrc->moments->Sy < 0.05) {
                continue;
            }
            if (tmpSrc->moments->Sx > SX_MAX) {
                continue;
            }
            if (tmpSrc->moments->Sy > SY_MAX) {
                continue;
            }
            if ((tmpSrc->moments->Sx / tmpSrc->moments->Sy) > AR_MAX) {
                continue;
            }
            if ((tmpSrc->moments->Sx / tmpSrc->moments->Sy) < AR_MIN) {
                continue;
            }
            if (tmpSrc->moments->SN < 30) {
                continue;
            }

            // for the moment, force splane dimensions to be 10x10 image pix
            binX = tmpSrc->moments->Sx/SCALE;
            if (binX < 0)
                continue;
            if (binX >= splane->numCols)
                continue;

            binY = tmpSrc->moments->Sy/SCALE;
            if (binY < 0)
                continue;
            if (binY >= splane->numRows)
                continue;

            splane->data.F32[binY][binX] += 1.0;
        }

        // find the peak in this image
        stats = psStatsAlloc (PS_STAT_MAX);
        stats = psImageStats (stats, splane, NULL, 0);
        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);

        psFree (splane);
        psFree (stats);

    }
    // XXX EAM : possible errors:
    //           1) no peak in splane
    //           2) no significant peak in splane

    // measure statistics on Sx, Sy if Sx, Sy within range of clump
    {
        pmPeak *clump;
        psF32 minSx, maxSx;
        psF32 minSy, maxSy;
        psVector *tmpSx = NULL;
        psVector *tmpSy = NULL;
        psStats *stats  = NULL;

        // XXX EAM : this lets us takes the single highest peak
        psArraySort (peaks, pmPeaksCompareDescend);
        clump = peaks->data[0];
        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);

        // define section window for clump
        minSx = clump->x * SCALE - 0.2;
        maxSx = clump->x * SCALE + 0.2;
        minSy = clump->y * SCALE - 0.2;
        maxSy = clump->y * SCALE + 0.2;

        tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);
        tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);
        tmpSx->n = 0;
        tmpSy->n = 0;

        // XXX clip sources based on flux?
        // create vectors with Sx, Sy values in window
        for (psS32 i = 0 ; i < sources->n ; i++)
        {
            pmSource *tmpSrc = (pmSource *) sources->data[i];

            if (tmpSrc->moments->Sx < minSx)
                continue;
            if (tmpSrc->moments->Sx > maxSx)
                continue;
            if (tmpSrc->moments->Sy < minSy)
                continue;
            if (tmpSrc->moments->Sy > maxSy)
                continue;
            tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;
            tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;
            tmpSx->n++;
            tmpSy->n++;
            if (tmpSx->n == tmpSx->nalloc) {
                psVectorRealloc (tmpSx, tmpSx->nalloc + 100);
                psVectorRealloc (tmpSy, tmpSy->nalloc + 100);
            }
        }

        // measures stats of Sx, Sy
        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);

        stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
        psfClump.X  = stats->clippedMean;
        psfClump.dX = stats->clippedStdev;

        stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
        psfClump.Y  = stats->clippedMean;
        psfClump.dY = stats->clippedStdev;

        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
        // these values should be pushed on the metadata somewhere

        psFree (stats);
        psFree (peaks);
        psFree (tmpSx);
        psFree (tmpSy);
    }

    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return (psfClump);
}

/******************************************************************************
    pmSourceRoughClass(source, metadata): make a guess at the source
    classification.
     
    XXX: push the clump info into the metadata?
     
    XXX: How can this function ever return FALSE?
     
    EAM: I moved S/N calculation to pmSourceMoments, using weight image
*****************************************************************************/

bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);

    psBool rc = true;

    int Nsat     = 0;
    int Next     = 0;
    int Nstar    = 0;
    int Npsf     = 0;
    int Ncr      = 0;
    int Nsatstar = 0;
    // psRegion allArray = psRegionSet (0, 0, 0, 0);
    psRegion inner;

    // report stats on S/N values for star-like objects
    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    starsn->n = 0;

    // check return status value (do these exist?)
    bool status;
    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");

    // XXX allow clump size to be scaled relative to sigmas?
    // make rough IDs based on clumpX,Y,DX,DY
    for (psS32 i = 0 ; i < sources->n ; i++) {

        pmSource *tmpSrc = (pmSource *) sources->data[i];

        tmpSrc->peak->type = 0;

        psF32 sigX = tmpSrc->moments->Sx;
        psF32 sigY = tmpSrc->moments->Sy;

        // XXX EAM : can we use the value of SATURATE if mask is NULL?
        // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
        inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2);
        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PM_MASK_SAT);

        // saturated star (size consistent with PSF or larger)
        // Nsigma should be user-configured parameter
        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
        big = true;
        if ((Nsatpix > 1) && big) {
            tmpSrc->type = PM_SOURCE_TYPE_STAR;
            tmpSrc->mode = PM_SOURCE_MODE_SATSTAR;
            Nsatstar ++;
            continue;
        }

        // saturated object (not a star, eg bleed trails, hot pixels)
        if (Nsatpix > 1) {
            tmpSrc->type = PM_SOURCE_TYPE_SATURATED;
            tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
            Nsat ++;
            continue;
        }

        // likely defect (too small to be stellar) (push out to 3 sigma)
        // low S/N objects which are small are probably stellar
        // only set candidate defects if
        if ((sigX < 0.05) || (sigY < 0.05)) {
            tmpSrc->type = PM_SOURCE_TYPE_DEFECT;
            tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
            Ncr ++;
            continue;
        }

        // likely unsaturated extended source (too large to be stellar)
        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
            tmpSrc->type = PM_SOURCE_TYPE_EXTENDED;
            tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
            Next ++;
            continue;
        }

        // the rest are probable stellar objects
        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
        starsn->n ++;
        Nstar ++;

        // PSF star (within 1.5 sigma of clump center, S/N > limit)
        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
            tmpSrc->type = PM_SOURCE_TYPE_STAR;
            tmpSrc->mode = PM_SOURCE_MODE_PSFSTAR;
            Npsf ++;
            continue;
        }

        // random type of star
        tmpSrc->type = PM_SOURCE_TYPE_STAR;
        tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
    }

    {
        psStats *stats  = NULL;
        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
        stats = psVectorStats (stats, starsn, NULL, NULL, 0);
        psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
        psFree (stats);
        psFree (starsn);
    }

    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
    psTrace (".pmObjects.pmSourceRoughClass", 2, "Next:     %3d\n", Next);
    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);

    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    return(rc);
}

/******************************************************************************
pmSourceMoments(source, radius): this function takes a subImage defined in the
pmSource data structure, along with the peak location, and determines the
various moments associated with that peak.
 
Requires the following to have been created:
    pmSource
    pmSource->peak
    pmSource->pixels
    pmSource->weight
    pmSource->mask
 
XXX: The peak calculations are done in image coords, not subImage coords.
 
XXX EAM : this version clips input pixels on S/N
XXX EAM : this version returns false for several reasons
*****************************************************************************/
# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)

bool pmSourceMoments(pmSource *source,
                     psF32 radius)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(source, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    PS_ASSERT_PTR_NON_NULL(source->mask, false);
    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);

    //
    // XXX: Verify the setting for sky if source->moments == NULL.
    //
    psF32 sky = 0.0;
    if (source->moments == NULL) {
        source->moments = pmMomentsAlloc();
    } else {
        sky = source->moments->Sky;
    }

    //
    // Sum = SUM (z - sky)
    // X1  = SUM (x - xc)*(z - sky)
    // X2  = SUM (x - xc)^2 * (z - sky)
    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    //
    psF32 peakPixel = -PS_MAX_F32;
    psS32 numPixels = 0;
    psF32 Sum = 0.0;
    psF32 Var = 0.0;
    psF32 X1 = 0.0;
    psF32 Y1 = 0.0;
    psF32 X2 = 0.0;
    psF32 Y2 = 0.0;
    psF32 XY = 0.0;
    psF32 x  = 0;
    psF32 y  = 0;
    psF32 R2 = PS_SQR(radius);

    psF32 xPeak = source->peak->x;
    psF32 yPeak = source->peak->y;
    psF32 xOff = source->pixels->col0 - source->peak->x;
    psF32 yOff = source->pixels->row0 - source->peak->y;

    // XXX why do I get different results for these two methods of finding Sx?
    // XXX Sx, Sy would be better measured if we clip pixels close to sky
    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
    // We loop through all pixels in this subimage (source->pixels), and for each
    // pixel that is not masked, AND within the radius of the peak pixel, we
    // proceed with the moments calculation.  need to do two loops for a
    // numerically stable result.  first loop: get the sums.
    // XXX EAM : mask == 0 is valid

    for (psS32 row = 0; row < source->pixels->numRows ; row++) {

        psF32 *vPix = source->pixels->data.F32[row];
        psF32 *vWgt = source->weight->data.F32[row];
        psU8  *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row];

        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
            if ((vMsk != NULL) && *vMsk) {
                vMsk++;
                continue;
            }

            psF32 xDiff = col + xOff;
            psF32 yDiff = row + yOff;

            // radius is just a function of (xDiff, yDiff)
            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
                if (vMsk != NULL)
                    vMsk++;
                continue;
            }

            psF32 pDiff = *vPix - sky;
            psF32 wDiff = *vWgt;

            // XXX EAM : check for valid S/N in pixel
            // XXX EAM : should this limit be user-defined?
            if (PS_SQR(pDiff) < wDiff) {
                if (vMsk != NULL)
                    vMsk++;
                continue;
            }

            Var += wDiff;
            Sum += pDiff;

            psF32 xWght = xDiff * pDiff;
            psF32 yWght = yDiff * pDiff;

            X1  += xWght;
            Y1  += yWght;

            XY  += xDiff * yWght;
            X2  += xDiff * xWght;
            Y2  += yDiff * yWght;

            peakPixel = PS_MAX (*vPix, peakPixel);
            numPixels++;
            if (vMsk != NULL)
                vMsk++;
        }
    }

    // if we have less than (1/4) of the possible pixels, force a retry
    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
        return (false);
    }

    psTrace (".psModules.pmSourceMoments", 5,
             "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);

    //
    // first moment X  = X1/Sum + xc
    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    // Sxy             = XY / Sum
    //
    x = X1/Sum;
    y = Y1/Sum;
    if ((fabs(x) > radius) || (fabs(y) > radius)) {
        psTrace (".psModules.pmSourceMoments", 3,
                 "large centroid swing; invalid peak %d, %d\n",
                 source->peak->x, source->peak->y);
        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
        return (false);
    }

    source->moments->x = x + xPeak;
    source->moments->y = y + yPeak;

    // XXX EAM : Sxy needs to have x*y subtracted
    source->moments->Sxy = XY/Sum - x*y;
    source->moments->Sum = Sum;
    source->moments->SN  = Sum / sqrt(Var);
    source->moments->Peak = peakPixel;
    source->moments->nPixels = numPixels;

    // XXX EAM : these values can be negative, so we need to limit the range
    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));

    psTrace (".psModules.pmSourceMoments", 4,
             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
             sky, Sum, source->moments->x, source->moments->y,
             source->moments->Sx, source->moments->Sy, source->moments->Sxy);

    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
    return(true);
}

pmModel *pmSourceSelectModel (pmSource *source)
{
    switch (source->type) {
    case PM_SOURCE_TYPE_STAR:
        return source->modelPSF;

    case PM_SOURCE_TYPE_EXTENDED:
        return source->modelEXT;

    default:
        return NULL;
    }
    return NULL;
}
