/** @file  pmGrowthCurveGenerate.c
 *
 *  Generate the curve-of-growth
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-01-27 06:39:38 $
 *
 *  Copyright 2004 Institute for Astronomy, University of Hawaii
 *
 */

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

/*****************************************************************************/
/* INCLUDE FILES                                                             */
/*****************************************************************************/

#include <strings.h>  // for strcasecmp
#include <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmSpan.h"
#include "pmFootprint.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 "pmModelClass.h"
#include "pmModelUtils.h"
#include "pmSourcePhotometry.h"
#include "pmFPAMaskWeight.h"
#include "psVectorBracket.h"
#include "pmErrorCodes.h"

pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc);

/*****************************************************************************/
/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
/*****************************************************************************/

// we generate the growth curve for the center of the image with the specified psf model
bool pmGrowthCurveGenerate (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal)
{
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_PTR_NON_NULL(readout->image, false);

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

    // XXX something of a hack: measure the growth curve at a number of points in the field and
    // average them together

    psArray *growths = psArrayAllocEmpty (100);

    for (float ix = -0.4; ix <= +0.4; ix += 0.2) {
	for (float iy = -0.4; iy <= +0.4; iy += 0.2) {

	    // use the center of the center pixel of the image
	    // 0.5 PIX: is this offset needed? probably -- the psf model uses 0.5,0.5 as the center, double check
	    float xc = (int)(ix*readout->image->numCols + 0.5*readout->image->numCols) + readout->image->col0 + 0.5;
	    float yc = (int)(iy*readout->image->numRows + 0.5*readout->image->numRows) + readout->image->row0 + 0.5;

	    pmGrowthCurve *growth = pmGrowthCurveForPosition (readout->image, psf, ignore, maskVal, markVal, xc, yc);
	    if (!growth) continue;

	    psArrayAdd (growths, 100, growth);
	    psFree (growth);
	}
    }
    psAssert (growths->n, "cannot build growth curve (psf model is invalid everywhere)");

    // just use a simple sample median to get the 'best' value from each growth curve...
    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);

    psVector *values = psVectorAlloc (growths->n, PS_DATA_F32);

    // median the values for the fitMags
    values->n = 0;
    for (int j = 0; j < growths->n; j++) {
	pmGrowthCurve *growth = growths->data[j];
	if (!isfinite(growth->fitMag)) continue;
	psVectorAppend (values, growth->fitMag);
    }
    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
	psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
	return false;
    }
    psf->growth->fitMag = stats->sampleMedian;

    // loop over a range of source fluxes
    // no need to interpolate since we have forced the object center
    // to 0.5, 0.5 above
    for (int i = 0; i < psf->growth->radius->n; i++) {

	// median the values for each radial bin
	values->n = 0;
	for (int j = 0; j < growths->n; j++) {
	    pmGrowthCurve *growth = growths->data[j];
	    if (!isfinite(growth->apMag->data.F32[i])) continue;
	    psVectorAppend (values, growth->apMag->data.F32[i]);
	}
	if (!psVectorStats (stats, values, NULL, NULL, 0)) {
	    psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
	    return false;
	}
	psf->growth->apMag->data.F32[i] = stats->sampleMedian;
    }
    psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius);
    psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef;

    psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f (fitMag - apMag @ ref : %f - %f)\n", psf->growth->apLoss, psf->growth->fitMag, psf->growth->apRef);

    psFree (growths);
    psFree (stats);
    psFree (values);

    return true;
}

pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc) {

    float fitMag, apMag;
    float radius;

    assert (psf->growth);

    float minRadius = psf->growth->radius->data.F32[0];
    pmGrowthCurve *growth = pmGrowthCurveAlloc (minRadius, psf->growth->maxRadius, psf->growth->refRadius);

    float dx = growth->maxRadius + 1;
    float dy = growth->maxRadius + 1;

    // create template model
    pmModel *modelRef = pmModelAlloc(psf->type);

    // assign the x and y coords to the image center
    // create an object with center intensity of 1000
    modelRef->params->data.F32[PM_PAR_SKY] = 0;
    modelRef->params->data.F32[PM_PAR_I0] = 1000;
    modelRef->params->data.F32[PM_PAR_XPOS] = xc;
    modelRef->params->data.F32[PM_PAR_YPOS] = yc;

    // create modelPSF from this model
    pmModel *model = pmModelFromPSF (modelRef, psf);
    if (!model) {
	psFree (growth);
	return NULL;
    }

    // measure the fitMag for this model
    pmSourcePhotometryModel (&fitMag, NULL, model);
    growth->fitMag = fitMag;

    // generate working image for this source
    psRegion region = {xc - dx, xc + dx, yc - dy, yc + dy};

    // force region to stop at dimensions of image
    region = psRegionForImage (image, region);

    // the view, image, and mask retain col0,row0
    psImage *view = psImageSubset (image, region);
    psImage *pixels = psImageCopy (NULL, view, PS_TYPE_F32);
    psImage *mask = psImageCopy (NULL, view, PS_TYPE_IMAGE_MASK);

    psImageInit (pixels, 0.0);
    psImageInit (mask, 0);

    // place the reference object in the image center
    // no need to mask the source here
    // XXX should we measure this for the analytical model only or the full model?
    pmModelAdd (pixels, NULL, model, PM_MODEL_OP_FULL, maskVal);

    // Loop over a range of radii.  No need to interpolate since we have forced the object
    // center to 0.5, 0.5 above
    for (int i = 0; i < growth->radius->n; i++) {

        radius = growth->radius->data.F32[i];

        // mask the given aperture and measure the apMag
        psImageKeepCircle (mask, xc, yc, radius, "OR", markVal);
        if (!pmSourcePhotometryAper (&apMag, model, pixels, mask, maskVal)) {
	    psFree (growth);
	    psFree (view);
	    psFree (pixels);
	    psFree (mask);
	    psFree (model);
	    psFree (modelRef);
	    return NULL;
        }
	psImageMaskPixels (mask, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask

        // the 'ignore' mode is for testing
        if (ignore) {
            growth->apMag->data.F32[i] = fitMag;
        } else {
            growth->apMag->data.F32[i] = apMag;
        }
    }
    
    psFree (view);
    psFree (pixels);
    psFree (mask);
    psFree (model);
    psFree (modelRef);

    // psLogMsg ("psModules", 4, "GrowthCurve for %f,%f\n", xc, yc);

    return growth;
}
