/** @file  pmFringeStats.c
 *
 *  @author Eugene Magnier, IfA
 *
 *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-04-12 23:50:45 $
 *
 *  Copyright 2004 IfA
 */

#include "pslib.h"
#include "pmFPA.h"
#include "pmFringeStats.h"

static void fringeRegionsFree(pmFringeRegions *fringe)
{
    psFree(fringe->x);
    psFree(fringe->y);
    psFree(fringe->mask);
    return;
}

pmFringeRegions *pmFringeRegionsAlloc(int nPts, int dX, int dY, int nX, int nY)
{
    pmFringeRegions *fringe = psAlloc(sizeof(pmFringeRegions));
    (void)psMemSetDeallocator(fringe, (psFreeFunc)fringeRegionsFree);

    fringe->x = NULL;
    fringe->y = NULL;
    fringe->mask = NULL;

    fringe->nRequested = nPts;
    fringe->nAccepted = 0;

    fringe->dX = dX;
    fringe->dY = dY;
    fringe->nX = nX;
    fringe->nY = nY;

    return fringe;
}


static void fringeStatsFree(pmFringeStats *stats)
{
    psFree(stats->regions);
    psFree(stats->f);
    psFree(stats->df);
}

pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions)
{
    pmFringeStats *stats = psAlloc(sizeof(pmFringeStats));
    (void)psMemSetDeallocator(stats, (psFreeFunc)fringeStatsFree);

    stats->regions = psMemIncrRefCounter(regions);
    stats->f = psVectorAlloc(regions->x->n, PS_TYPE_F32);
    stats->df = psVectorAlloc(regions->x->n, PS_TYPE_F32);

    return stats;
}


static void fringeScaleFree(pmFringeScale *scale)
{
    psFree(scale->coeff);
    psFree(scale->coeffErr);
    return;
}

pmFringeScale *pmFringeScaleAlloc(int nFringeFrames)
{
    pmFringeScale *scale = psAlloc(sizeof(pmFringeScale));
    (void)psMemSetDeallocator(scale, (psFreeFunc)fringeScaleFree);

    scale->nFringeFrames = nFringeFrames;
    scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
    scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);

    return scale;
}

bool pmFringeStatsCreatePoints(pmFringeRegions *fringe, psImage *image)
{

    double frnd;
    // create fringe->nRequested

    psRandom *rnd = psRandomAlloc(PS_RANDOM_TAUS, 0);

    fringe->x = psVectorRecycle(fringe->x, fringe->nRequested, PS_TYPE_F32);
    fringe->y = psVectorRecycle(fringe->y, fringe->nRequested, PS_TYPE_F32);
    fringe->mask = psVectorRecycle(fringe->mask, fringe->nRequested, PS_TYPE_U8);

    int nX = image->numCols;
    int nY = image->numRows;

    psF32 *xPt = fringe->x->data.F32;
    psF32 *yPt = fringe->y->data.F32;

    int dX = fringe->dX;
    int dY = fringe->dY;

    // generate random points located within image bounds
    for (int i = 0; i < fringe->nRequested; i++) {
        frnd = psRandomUniform(rnd);
        xPt[i] = (nX - 2*dX)* frnd + dX;
        frnd = psRandomUniform(rnd);
        yPt[i] = (nY - 2*dY)* frnd + dY;
    }
    return true;
}

pmFringeStats *pmFringeStatsMeasure(pmFringeRegions *fringe, pmReadout *readout, psMaskType maskVal)
{
    if (!fringe->x || !fringe->y) {
        // create the fringe vectors for this image
        pmFringeStatsCreatePoints(fringe, readout->image);
    }

    PS_ASSERT_PTR_NON_NULL(fringe->x, false);
    PS_ASSERT_PTR_NON_NULL(fringe->y, false);

    pmFringeStats *measurements = pmFringeStatsAlloc(fringe);

    psF32 *xPt = fringe->x->data.F32;
    psF32 *yPt = fringe->y->data.F32;
    psF32 *fPt = measurements->f->data.F32;
    psF32 *dfPt = measurements->df->data.F32;

    int dX = fringe->dX;
    int dY = fringe->dY;

    psImage *image = readout->image;
    psImage *mask  = readout->mask;

    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);

    for (int i = 0; i < fringe->x->n; i++) {
        psRegion region = psRegionSet(xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, yPt[i] + dY + 1);
        region = psRegionForImage(image, region);
        psImage *subImage = psImageSubset(image, region);
        psImage *subMask = psImageSubset(mask, region);
        psImageStats(stats, subImage, subMask, maskVal);

        fPt[i] = stats->sampleMedian;
        dfPt[i] = stats->sampleStdev;
    }

    return measurements;
}

// XXX include the fringe error (fringe->df) in the fit?
pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes)
{
    double sum;
    psVector *v1;
    psVector *v2;

    pmFringeScale *scale = pmFringeScaleAlloc(fringes->n);

    int nCof = fringes->n + 1;
    int nPts = science->f->n;

    psImage *A = psImageAlloc(nCof, nCof, PS_TYPE_F64);
    psVector *B = psVectorAlloc(nCof, PS_TYPE_F64);

    for (int i = 0; i < nCof; i++) {
        psF32 *p1 = NULL;
        if (i != 0) {
            v1 = fringes->data[i - 1];
            p1 = v1->data.F32;
        }
        for (int j = 0; j < nCof; j++) {
            psF32 *p2 = NULL;
            if (j != 0) {
                v2 = fringes->data[j - 1];
                p2 = v2->data.F32;
            }

            sum = 0;
            for (int k = 0; k < nPts; k++) {
                psF32 f1 = (p1 == NULL) ? 1.0 : p1[k];
                psF32 f2 = (p2 == NULL) ? 1.0 : p2[k];
                sum += f1*f2;
            }
            A->data.F32[i][j] = sum;
        }

        sum = 0;
        psF32 *p2 = science->f->data.F32;
        for (int j = 0; j < nPts; j++) {
            psF32 f1 = (p1 == NULL) ? 1.0 : p1[j];
            psF32 f2 = p2[j];
            sum += f1*f2;
        }
        B->data.F32[i] = sum;
    }

    if (! psGaussJordan(A, B)) {
        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    }


    for (int i = 0; i < nCof; i++) {
        scale->coeff->data.F32[i] = B->data.F32[i];
        scale->coeffErr->data.F32[i] = sqrt(A->data.F32[i][i]);
    }

    return scale;
}

// XXX note that this modifies the input fringe images
psImage *pmFringeCorrect(pmReadout *readout, pmFringeRegions *fringes, psArray *fringeImages, psArray *fringeStats, psMaskType maskVal)
{
    // measure the fringe stats for the science frame and solve for the scales
    pmFringeStats *scienceStats = pmFringeStatsMeasure(fringes, readout, maskVal);
    pmFringeScale *scale = pmFringeScaleMeasure(scienceStats, fringeStats);
    psFree(scienceStats);

    // build the fringe correction image
    // XXX we could save data space by making the first image the output image
    psImage *sumFringe = NULL;
    for (int i = 0; i < fringeImages->n; i++) {

        // rescale the fringe image
        psBinaryOp(fringeImages->data[i], fringeImages->data[i], "*",
                   psScalarAlloc(scale->coeff->data.F32[i+1], PS_TYPE_F32));

        // sum together
        sumFringe = (psImage*)psBinaryOp(sumFringe, sumFringe, "+", fringeImages->data[i]);
    }
    psFree(scale);

    // subtract the resulting fringe frame
    readout->image = (psImage*)psBinaryOp(readout->image, readout->image, "-", sumFringe);

    return sumFringe;
}

