/** @file  pmSourcePhotometry.c
 *
 *  @author EAM, IfA; GLG, MHPCC
 *
 *  @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 "pmGrowthCurve.h"
#include "pmModel.h"
#include "pmPSF.h"
#include "pmSource.h"
#include "pmModelGroup.h"
#include "pmSourcePhotometry.h"

# define DO_SKY 0

static float AP_MIN_SN = 0.0;

bool pmSourceMagnitudesInit (psMetadata *config)
{

    bool status;

    float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN");
    if (status) {
        AP_MIN_SN = limit;
    }
    return true;
}

/**
    this function is used to calculate the three defined source magnitudes:
    - apMag  : only if S/N > AP_MIN_SN
             : is optionally corrected for curve-of-growth if:
        - the source is a STAR (PSF)
        - the option is selected (how??)
    - psfMag : all sources with non-NULL modelPSF
             : is optionally corrected for aperture residual if:
        - the source is a STAR (PSF)
        - the option is selected (how??)
 
    - extMag : all sources with non-NULL modelEXT
**/

// XXX EAM : masked region should be (optionally) elliptical
bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode)
{

    int status = false;
    bool isPSF;
    float x, y;
    float rflux;
    float SN;
    pmModel *model;

    source->psfMag = 99.0;
    source->extMag = 99.0;
    source->errMag = 99.0;
    source->apMag  = 99.0;

    switch (source->type) {
    case PM_SOURCE_TYPE_STAR:
        model = source->modelPSF;
        if (model == NULL)
            return false;
        isPSF = true;
        break;

    case PM_SOURCE_TYPE_EXTENDED:
        model = source->modelEXT;
        if (model == NULL)
            return false;
        isPSF = false;
        break;

    default:
        return false;
    }

    if (model->dparams->data.F32[1] > 0) {
        SN = model->params->data.F32[1] / model->dparams->data.F32[1];
        source->errMag = 1.0 / SN;
    } else {
        SN = 0.0;
        source->errMag = 0.0;
    }
    x = model->params->data.F32[2];
    y = model->params->data.F32[3];

    // measure object model photometry
    status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf) {
        // convert to the equivalent 2D model?
        source->psfMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0);
    }

    if (SN < AP_MIN_SN)
        return false;

    // replace source flux
    // XXX test to see if source has been subtracted?
    pmModelAdd (source->pixels, source->mask, model, false, false);

    // set aperture mask circle to model radius
    psImageKeepCircle (source->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);

    // measure the weight of included pixels
    if (mode & PM_SOURCE_PHOT_WEIGHT) {
        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->mask);
    }

    // measure object aperture photometry
    status = pmSourcePhotometryAper  (&source->apMag, model, source->pixels, source->mask);

    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
    if (status && isPSF && psf) {
        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
            source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusTMP);
        }
        if (mode & PM_SOURCE_PHOT_APCORR) {
            rflux   = pow (10.0, 0.4*source->psfMag);
            source->apMag  -= PS_SQR(model->radiusTMP)*rflux * psf->skyBias + psf->skySat / rflux;
        }
    }

    // unmask aperture
    psImageKeepCircle (source->mask, x, y, model->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));

    // subtract object, leave local sky
    pmModelSub (source->pixels, source->mask, model, false, false);

    return status;
}

/*
aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y)
(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y)
(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y)
 
*/

// return source model magnitude
bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
{

    float fitSum = 0;
    *fitMag = 99.0;

    if (model == NULL) {
        return false;
    }

    // measure fitMag
    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
    fitSum = modelFluxFunc (model->params);
    if (fitSum <= 0)
        return false;
    if (!isfinite(fitSum))
        return false;
    *fitMag = -2.5*log10(fitSum);

    return (true);
}

// return source aperture magnitude
bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask)
{
    float apSum = 0;
    float sky = 0;
    *apMag = 99.0;

    if (model == NULL) {
        return false;
    }

    if (DO_SKY) {
        sky = model->params->data.F32[0];
    } else {
        sky = 0;
    }

    // measure apMag
    for (int ix = 0; ix < image->numCols; ix++) {
        for (int iy = 0; iy < image->numRows; iy++) {
            if (mask->data.U8[iy][ix])
                continue;
            apSum += image->data.F32[iy][ix] - sky;
        }
    }
    if (apSum <= 0)
        return false;

    *apMag = -2.5*log10(apSum);
    return (true);
}

// return source aperture magnitude
bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask)
{
    float modelSum = 0;
    float validSum = 0;
    float sky = 0;
    float value;

    int Xo, Yo, dP;
    int dX, DX, NX;
    int dY, DY, NY;

    *pixWeight = 0.0;

    if (model == NULL)
        return false;

    // we only care about the value of the object model, not the local sky
    if (DO_SKY) {
        sky = model->params->data.F32[0];
    } else {
        sky = 0;
    }

    // the model function returns the source flux at a position
    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    coord->n = 2;

    psVector *params = model->params;

    Xo = params->data.F32[2];
    Yo = params->data.F32[3];

    dX = Xo - image->col0;
    dP = image->numCols - dX;
    DX = PS_MAX(dX, dP);
    NX = image->numCols;

    dY = Yo - image->row0;
    dP = image->numRows - dY;
    DY = PS_MAX(dY, dP);
    NY = image->numRows;

    // measure modelSum and validSum
    // XXX this does not work for sources near the edge: we need to measure for
    // a full square region
    for (int ix = -DX; ix < DX + 1; ix++) {
        int mx = ix + dX;
        for (int iy = -DY; iy < DY + 1; iy++) {
            int my = iy + dY;

            coord->data.F32[0] = (psF32) (ix + Xo);
            coord->data.F32[1] = (psF32) (iy + Yo);

            // for the full model, add all points
            value = modelFunc (NULL, params, coord) - sky;
            modelSum += value;

            // include count only the unmasked pixels within the image area
            if (mx < 0)
                continue;
            if (my < 0)
                continue;
            if (mx >= NX)
                continue;
            if (my >= NY)
                continue;
            if (mask->data.U8[my][mx])
                continue;

            validSum += value;
        }
    }
    psFree (coord);

    if (validSum <= 0)
        return false;

    *pixWeight = validSum / modelSum;
    return (true);
}

float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj)
{

    int Xs, Xe, Ys, Ye;
    int xi, xj, yi, yj;
    int xIs, xJs, yIs, yJs;
    int xIe, yIe;
    float flux, wt;

    psImage *Pi = Mi->pixels;
    psImage *Pj = Mj->pixels;

    psImage *Wi = Mi->weight;

    psImage *Ti = Mi->mask;
    psImage *Tj = Mj->mask;

    Xs = PS_MAX (Pi->col0, Pj->col0);
    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);

    Ys = PS_MAX (Pi->row0, Pj->row0);
    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);

    xIs = Xs - Pi->col0;
    xJs = Xs - Pj->col0;
    yIs = Ys - Pi->row0;
    yJs = Ys - Pj->row0;

    xIe = Xe - Pi->col0;
    yIe = Ye - Pi->row0;

    // note that this is addressing the same image pixels,
    // though only if both are source not model images
    flux = 0;
    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
            if (Ti->data.U8[yi][xi])
                continue;
            if (Tj->data.U8[yj][xj])
                continue;
            wt = Wi->data.F32[yi][xi];
            if (wt > 0) {
                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
            }
        }
    }
    return (flux);
}

float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj)
{

    int Xs, Xe, Ys, Ye;
    int xi, xj, yi, yj;
    int xIs, xJs, yIs, yJs;
    int xIe, yIe;
    float flux, wt;

    psImage *Pi = Mi->pixels;
    psImage *Pj = Mj->pixels;

    psImage *Wi = Mi->weight;

    psImage *Ti = Mi->mask;
    psImage *Tj = Mj->mask;

    Xs = PS_MAX (Pi->col0, Pj->col0);
    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);

    Ys = PS_MAX (Pi->row0, Pj->row0);
    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);

    xIs = Xs - Pi->col0;
    xJs = Xs - Pj->col0;
    yIs = Ys - Pi->row0;
    yJs = Ys - Pj->row0;

    xIe = Xe - Pi->col0;
    yIe = Ye - Pi->row0;

    // note that this is addressing the same image pixels,
    // though only if both are source not model images
    flux = 0;
    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
            if (Ti->data.U8[yi][xi])
                continue;
            if (Tj->data.U8[yj][xj])
                continue;
            wt = Wi->data.F32[yi][xi];
            if (wt > 0) {
                flux += 1.0 / wt;
            }
        }
    }
    return (flux);
}

