/** @file  pmFringeStats.c
 *
 *  @author Eugene Magnier, IfA
 *
 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-06-21 03:21:16 $
 *
 *  Copyright 2004 IfA
 */

#include <stdio.h>
#include <assert.h>
#include <pslib.h>
#include "pmFPA.h"
#include "pmFringeStats.h"

#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))


// Future optimisations for speed:
//
// 1. Clipping --- don't re-do the matrix setup again, but carry matrix and vector around, subtract
// contributions from clipped data points.
// 2. Faster psImageStats (use memcpy?)



//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// pmFringeRegions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, const psImage *image)
{
    PS_ASSERT_PTR_NON_NULL(fringe, false);
    PS_ASSERT_PTR_NON_NULL(image, false);
    PS_ASSERT_IMAGE_NON_EMPTY(image, false);

    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);
    fringe->x->n = fringe->y->n = fringe->mask->n = fringe->nRequested;
    psVectorInit(fringe->mask, 0);

    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
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions)
{
    PS_ASSERT_PTR_NON_NULL(regions, false);

    pmFringeStats *stats = psAlloc(sizeof(pmFringeStats));
    (void)psMemSetDeallocator(stats, (psFreeFunc)fringeStatsFree);

    int numRegions = regions->x->n;     // Number of regions
    stats->regions = psMemIncrRefCounter(regions);
    stats->f = psVectorAlloc(numRegions, PS_TYPE_F32);
    stats->df = psVectorAlloc(numRegions, PS_TYPE_F32);
    stats->f->n = stats->df->n = numRegions;

    return stats;
}

pmFringeStats *pmFringeStatsMeasure(pmFringeRegions *fringe, pmReadout *readout, psMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(fringe, NULL);
    PS_ASSERT_PTR_NON_NULL(readout, NULL);
    PS_ASSERT_PTR_NON_NULL(readout->image, NULL);
    PS_ASSERT_IMAGE_NON_EMPTY(readout->image, NULL);

    if (!fringe->x || !fringe->y) {
        // create the fringe vectors for this image
        pmFringeRegionsCreatePoints(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 *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Median statistics only
    psStats *medianSd = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); // Median and SD

    // Measure the sky over the image
    psImage *sky = psImageAlloc(fringe->nX, fringe->nY, PS_TYPE_F32);
    for (int i = 0; i < fringe->nY; i++) {
        int y0 = (float)i * (float)image->numRows / (float)fringe->nY;
        int y1 = (float)(i + 1) * (float)image->numRows / (float)fringe->nY;
        for (int j = 0; j < fringe->nX; j++) {
            int x0 = (float)j * (float)image->numCols / (float)fringe->nX;
            int x1 = (float)(j + 1) * (float)image->numCols / (float)fringe->nX;
            psRegion region = psRegionSet(x0, x1, y0, y1);
            psImage *subImage = psImageSubset(image, region); // Subimage of the sky region
            psImage *subMask = psImageSubset(mask, region); // Subimage of the sky region
            psImageStats(median, subImage, subMask, maskVal);
            sky->data.F32[i][j] = median->sampleMedian;
        }
    }

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

        int xSky = xPt[i] / (float)image->numCols * (float)sky->numCols;
        int ySky = yPt[i] / (float)image->numRows * (float)sky->numRows;

        fPt[i] = medianSd->sampleMedian - sky->data.F32[ySky][xSky];
        dfPt[i] = 1.0 / medianSd->sampleStdev;

        psTrace(__func__, 7, "[%d:%d,%d:%d]: %f %f\n", (int)region.x0, (int)region.x1,
                (int)region.y0, (int)region.y1, fPt[i], dfPt[i]);
    }

    return measurements;
}


bool pmFringeStatsWriteFits(psFits *fits, const pmFringeStats *fringe, const char *extname)
{
    // Make sure the input is well-behaved
    PS_ASSERT_PTR_NON_NULL(fits, false);
    PS_ASSERT_PTR_NON_NULL(fringe, false);
    pmFringeRegions *regions = fringe->regions; // The fringe regions
    PS_ASSERT_PTR_NON_NULL(regions, false);
    psVector *x = regions->x;           // The x positions
    psVector *y = regions->y;           // The y positions
    psVector *mask = regions->mask;     // The region mask
    int numRows = x->n;                 // Number of rows in the table
    PS_ASSERT_INT_POSITIVE(regions->nRequested, false);
    PS_ASSERT_VECTOR_NON_NULL(x, false);
    PS_ASSERT_VECTOR_NON_NULL(y, false);
    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
    PS_ASSERT_VECTOR_SIZE(y, numRows, false);
    if (mask) {
        PS_ASSERT_VECTOR_NON_NULL(mask, false);
        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
        PS_ASSERT_VECTOR_SIZE(mask, numRows, false);
    }
    psVector *f = fringe->f;            // The fringe measurements
    psVector *df = fringe->df;      // The fringe standard deviatiations
    PS_ASSERT_VECTOR_NON_NULL(f, false);
    PS_ASSERT_VECTOR_NON_NULL(df, false);
    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, false);
    PS_ASSERT_VECTOR_TYPE(df, PS_TYPE_F32, false);
    PS_ASSERT_VECTOR_SIZE(f, numRows, false);
    PS_ASSERT_VECTOR_SIZE(df, numRows, false);

    if (extname && strlen(extname) > 0) {
        if (!psFitsMoveExtName(fits, extname)) {
            psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname);
            return false;
        }
    } else if (!psFitsMoveExtNum(fits, 0, false)) {
        psError(PS_ERR_IO, false, "Unable to move to PHU\n");
        return false;
    }

    // We need to write:
    // Scalars: dX, dY, nX, nY
    // Vectors: x, y, mask, f, df

    psMetadata *scalars = psMetadataAlloc(); // Metadata to hold the scalars; will be the header
    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width",
                     regions->dX);
    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height",
                     regions->dY);
    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x",
                     regions->nX);
    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y",
                     regions->nY);

    psArray *table = psArrayAlloc(numRows); // The table
    table->n = numRows;
    // Translate the vectors into the required format for psFitsWriteTable()
    for (long i = 0; i < numRows; i++) {
        psMetadata *row = psMetadataAlloc();
        psArraySet(table, i, row);
        psMetadataAddF32(row, PS_LIST_TAIL, "x", PS_META_REPLACE, "Fringe position in x", x->data.F32[i]);
        psMetadataAddF32(row, PS_LIST_TAIL, "y", PS_META_REPLACE, "Fringe position in y", y->data.F32[i]);
        if (mask) {
            psMetadataAddU8(row, PS_LIST_TAIL, "mask", PS_META_REPLACE, "Mask", mask->data.F32[i]);
        } else {
            psMetadataAddU8(row, PS_LIST_TAIL, "mask", PS_META_REPLACE, "Mask", 0);
        }
        psMetadataAddF32(row, PS_LIST_TAIL, "f", PS_META_REPLACE, "Fringe measurement", f->data.F32[i]);
        psMetadataAddF32(row, PS_LIST_TAIL, "df", PS_META_REPLACE, "Fringe stdev", df->data.F32[i]);
        psFree(row);                    // Drop reference
    }

    bool success;                       // Success of operation
    if (!(success = psFitsWriteTable(fits, scalars, table, extname))) {
        psError(PS_ERR_IO, false, "Unable to write fringe data to extension %s\n", extname);
    }
    psFree(scalars);
    psFree(table);

    return success;
}

pmFringeStats *pmFringeStatsReadFits(const psFits *fits, const char *extname)
{
    PS_ASSERT_PTR_NON_NULL(fits, NULL);

    if (extname && strlen(extname) > 0) {
        if (!psFitsMoveExtName(fits, extname)) {
            psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname);
            return NULL;
        }
    } else if (!psFitsMoveExtNum(fits, 0, false)) {
        psError(PS_ERR_IO, false, "Unable to move to PHU\n");
        return NULL;
    }

    psMetadata *header = psFitsReadHeader(NULL, fits); // The FITS header
    if (!header) {
        psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname);
        return NULL;
    }

    // Read the scalars from the header
    #define READ_SCALAR(SCALAR, NAME) \
    int SCALAR = psMetadataLookupS32(&mdok, header, NAME); \
    if (!mdok || SCALAR <= 0) { \
        psError(PS_ERR_IO, true, "Unable to find " NAME " in header of extension %s.\n", extname); \
        return NULL; \
    }

    // Need to retrieve the scalars: dX, dY, nX, nY
    bool mdok = true;                   // Status of MD lookup
    READ_SCALAR(dX, "PSFRNGDX");
    READ_SCALAR(dY, "PSFRNGDX");
    READ_SCALAR(nX, "PSFRNGDX");
    READ_SCALAR(nY, "PSFRNGDX");
    psFree(header);

    // Now the vectors: x, y, mask, f, df
    psArray *table = psFitsReadTable(fits); // The table
    long numRows = table->n;            // Number of rows

    pmFringeRegions *regions = pmFringeRegionsAlloc(numRows, dX, dY, nX, nY); // The fringe regions
    psVector *x = psVectorAlloc(numRows, PS_TYPE_F32); // x position
    psVector *y = psVectorAlloc(numRows, PS_TYPE_F32); // y position
    psVector *mask = psVectorAlloc(numRows, PS_TYPE_U8); // mask
    x->n = y->n = mask->n = numRows;
    regions->x = x;
    regions->y = y;
    regions->mask = mask;
    pmFringeStats *fringes = pmFringeStatsAlloc(regions); // The fringe measurements
    psFree(regions);                    // Drop reference
    psVector *f = fringes->f;           // fringe measurement
    psVector *df = fringes->df;     // fringe stdev

    #define READ_ROW(VECTOR, TYPE, NAME, DESCRIPTION) \
    VECTOR->data.TYPE[i] = psMetadataLookup##TYPE(&mdok, row, NAME); \
    if (!mdok) { \
        psError(PS_ERR_IO, true, "Unable to find " #DESCRIPTION " .\n"); \
        psFree(table); \
        psFree(fringes); \
        return NULL; \
    }

    // Translate the table into vectors
    for (long i = 0; i < numRows; i++) {
        psMetadata *row = table->data[i]; // Table row
        READ_ROW(x, F32, "x", "x position");
        READ_ROW(y, F32, "y", "y position");
        READ_ROW(mask, U8, "mask", "mask value");
        READ_ROW(f, F32, "f", "fringe measurement");
        READ_ROW(df, F32, "df", "fringe standard deviation");
    }
    psFree(table);

    return fringes;
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// pmFringeScale
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

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);
    scale->coeff->n = nFringeFrames + 1;
    scale->coeff->n = nFringeFrames + 1;

    return scale;
}

// Determine the fringe scales through solving the least-squares problem
static bool scaleMeasure(pmFringeScale *scale, // Scale to return
                         pmFringeStats *science, // The fringe measurements for the science image
                         psArray *fringes // Array of fringe measurements for the templates
                        )
{
    assert(scale);
    assert(science);
    assert(fringes);
    assert(scale->nFringeFrames == fringes->n);

    psVector *mask = science->regions->mask; // The region mask

    int numCoeffs = fringes->n + 1;     // Number of coefficients: scales for the templates plus a background
    int numPoints = science->regions->nRequested; // Number of points (i.e., fringe measurements)

    psImage *A = psImageAlloc(numCoeffs, numCoeffs, PS_TYPE_F64); // The least-squares matrix
    psVector *B = psVectorAlloc(numCoeffs, PS_TYPE_F64); // The least-squares vector
    B->n = numCoeffs;

    // Generate the least-squares matrix and vector
    for (int i = 0; i < numCoeffs; i++) {
        psVector *fringe1 = NULL;       // A fringe measurement
        if (i != 0) {
            pmFringeStats *fringe = fringes->data[i - 1];
            fringe1 = fringe->f;
        }

        // Fill in the upper part of the matrix
        for (int j = i; j < numCoeffs; j++) {
            psVector *fringe2 = NULL;   // Another fringe measurement
            if (j != 0) {
                pmFringeStats *fringe = fringes->data[j - 1];
                fringe2 = fringe->f;
            }

            double matrix = 0.0;        // The matrix sum
            for (int k = 0; k < numPoints; k++) {
                if (!mask->data.U8[k]) {
                    psF32 f1 = (fringe1) ? fringe1->data.F32[k] : 1.0; // Contribution from i fringe
                    psF32 f2 = (fringe2) ? fringe2->data.F32[k] : 1.0; // Contribution from j fringe
                    psF32 dsInv = science->df->data.F32[k]; // 1 / sigma
                    matrix += f1 * f2 * dsInv * dsInv;
                }
            }
            A->data.F64[i][j] = matrix;
        }

        // Use symmetry to fill in the lower part of the matrix
        for (int j = 0; j < i; j++) {
            A->data.F64[i][j] = A->data.F64[j][i];
        }

        double vector = 0.0;            // The vector sum
        for (int k = 0; k < numPoints; k++) {
            if (!mask->data.U8[k]) {
                psF32 f1 = (fringe1) ? fringe1->data.F32[k] : 1.0; // Contribution from fringe 1
                psF32 s = science->f->data.F32[k]; // Contribution from science measurement
                psF32 dsInv = science->df->data.F32[k]; // 1 / sigma
                vector += f1 * s * dsInv * dsInv;
            }
        }
        B->data.F64[i] = vector;
    }

    if (psTraceGetLevel(__func__) >= 5) {
        printf("From %d points:\n", numPoints);
        for (int i = 0; i < numCoeffs; i++) {
            for (int j = 0; j < numCoeffs; j++) {
                printf("%.2e ", A->data.F64[i][j]);
            }
            printf("\n");
        }
    }

    // Solve the least-squares equation
    if (! psMatrixGJSolve(A, B)) {
        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
        return false;
    }

    // Copy the results over
    for (int i = 0; i < numCoeffs; i++) {
        scale->coeff->data.F32[i] = B->data.F64[i];
        scale->coeffErr->data.F32[i] = sqrt(A->data.F64[i][i]);
    }

    return true;
}

// Measure the fringe differences for each region
static bool fringeScaleDiffs(psVector *diff, // Vector of differences
                             pmFringeStats *science, // Science fringe measurements
                             psArray *fringes, // Template fringe measurements
                             pmFringeScale *scale // Fringe scales
                            )
{
    assert(diff);
    assert(diff->type.type == PS_TYPE_F32);
    assert(science);
    assert(fringes);
    assert(scale);
    assert(diff->n == science->regions->nRequested);
    assert(fringes->n == scale->nFringeFrames);

    psVector *mask = science->regions->mask; // The region mask

    for (int i = 0; i < diff->n; i++) {
        if (!mask->data.U8[i]) {
            float difference = science->f->data.F32[i] - scale->coeff->data.F32[0];
            for (int j = 0; j < fringes->n; j++) {
                pmFringeStats *fringe = fringes->data[j]; // The fringe of interest
                difference -= scale->coeff->data.F32[j + 1] * fringe->f->data.F32[i];
            }
            diff->data.F32[i] = difference * difference * science->df->data.F32[i] * science->df->data.F32[i];
        }
    }

    return true;
}

// Clip regions based on the differences; return the number masked
static int clipRegions(psVector *diffs, // Differences
                       psVector *mask,  // Region mask
                       float rej        // Rejection limit in standard deviations
                      )
{
    assert(diffs);
    assert(diffs->type.type == PS_TYPE_F32);
    assert(mask);
    assert(mask->type.type == PS_TYPE_U8);
    assert(diffs->n == mask->n);

    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics
    psVectorStats(stats, diffs, NULL, mask, 1);
    float middle = stats->sampleMedian; // The middle of the distribution
    float thresh = rej * 0.74 * (stats->sampleUQ - stats->sampleLQ); // The rejection threshold
    psFree(stats);

    int numClipped = 0;                 // Number clipped
    for (int i = 0; i < diffs->n; i++) {
        psTrace(__func__, 10, "Region %d (%d): %f\n", i, mask->data.U8[i], diffs->data.F32[i]);
        if (!mask->data.U8[i] && fabs(diffs->data.F32[i]) > middle + thresh) {
            psTrace(__func__, 5, "Masking %d: %f\n", i, diffs->data.F32[i]);
            mask->data.U8[i] = 1;
            numClipped++;
        }
    }

    return numClipped;
}


// XXX include the fringe error (fringe->df) in the fit?
pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes, float rej,
                                    unsigned int nIter, float keepFrac)
{
    PS_ASSERT_PTR_NON_NULL(science, NULL);
    PS_ASSERT_PTR_NON_NULL(fringes, NULL);
    PS_ASSERT_INT_POSITIVE(fringes->n, NULL);
    PS_ASSERT_INT_POSITIVE(nIter, NULL);

    pmFringeRegions *regions = science->regions; // The fringe regions
    int numRegions = regions->nRequested; // Number of regions

    // Ensure we are dealing with the SAME fringe points for all the inputs.
    // Otherwise, we're going to get crazy results.
    for (long i = 0; i < numRegions; i++) {
        float xScience = regions->x->data.F32[i]; // The x position for the science image
        float yScience = regions->y->data.F32[i]; // The y position for the science image
        for (long j = 0; j < fringes->n; j++) {
            pmFringeStats *fringe = fringes->data[j]; // The fringe statistics from a fringe image
            pmFringeRegions *fringeRegions = fringe->regions; // The fringe regions for that fringe image
            if (fringeRegions->x->data.F32[i] != xScience ||
                    fringeRegions->y->data.F32[i] != yScience) {
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Science and fringe measurement regions "
                        "don't match.\n");
                return NULL;
            }
        }
    }

    // Set up the mask
    if (!regions->mask) {
        regions->mask = psVectorAlloc(numRegions, PS_TYPE_U8);
        regions->mask->n = numRegions;
        psVectorInit(regions->mask, 0);
    }
    psVector *mask = regions->mask;     // The region mask
    psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Median statistics
    unsigned int numClipped = 0;        // Total number clipped
    psVector *diff = psVectorAlloc(numRegions, PS_TYPE_F32); // The differences between obs. and pred.
    diff->n = numRegions;

    pmFringeScale *scale = pmFringeScaleAlloc(fringes->n); // The fringe scales

    // Get rid of bad data points
    for (int i = 0; i < fringes->n; i++) {
        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
        for (int j = 0; j < numRegions; j++) {
            if (!finite(fringe->f->data.F32[j])) {
                mask->data.U8[j] = 1;
                psTrace(__func__, 9, "Masking region %d because not finite in fringe %d.\n", j, i);
            }
        }
    }

    // Get rid of the extreme outliers by assuming most of the points are somewhat clustered
    psVectorStats(median, science->f, NULL, NULL, 0);
    scale->coeff->data.F32[0] = median->sampleMedian;
    for (int i = 0; i < fringes->n; i++) {
        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
        psVectorStats(median, fringe->f, NULL, NULL, 0);
        scale->coeff->data.F32[0] -= median->sampleMedian;
        scale->coeff->data.F32[i] = 0.0;
    }
    psFree(median);
    fringeScaleDiffs(diff, science, fringes, scale);
    numClipped = clipRegions(diff, mask, 3.0*rej);
    psTrace(__func__, 4, "%d regions clipped in initial pass.\n", numClipped);

    unsigned int iter = 0;              // Iteration number
    unsigned int iterClip = 0;          // Number clipped in this iteration
    do {
        iter++;
        scaleMeasure(scale, science, fringes); // The scales
        psTrace(__func__, 1, "Fringe scales after iteration %d:\n", iter);
        psTrace(__func__, 1, "Background: %f %f\n", scale->coeff->data.F32[0], scale->coeffErr->data.F32[0]);
        for (int i = 0; i < scale->nFringeFrames; i++) {
            psTrace(__func__, 1, "%d: %f %f\n", i, scale->coeff->data.F32[i + 1],
                    scale->coeffErr->data.F32[i + 1]);
        }

        fringeScaleDiffs(diff, science, fringes, scale);
        iterClip = clipRegions(diff, mask, rej); // Number clipped
        numClipped += iterClip;
        psTrace(__func__, 9, "Clipped: %d\tFrac: %f\n", iterClip, (float)numClipped/(float)numRegions);
    } while (iterClip > 0 && iter < nIter && (float)numClipped/(float)numRegions <= 1.0 - keepFrac);
    psFree(diff);

    // A final iteration with the last clipping
    scaleMeasure(scale, science, fringes);

    return scale;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Fringe correction
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

// XXX note that this modifies the input fringe images
psImage *pmFringeCorrect(pmReadout *readout, pmFringeRegions *fringes, psArray *fringeImages,
                         psArray *fringeStats, psMaskType maskVal, float rej,
                         unsigned int nIter, float keepFrac)
{
    PS_ASSERT_PTR_NON_NULL(readout, NULL);
    PS_ASSERT_PTR_NON_NULL(readout->image, NULL);
    PS_ASSERT_IMAGE_NON_EMPTY(readout->image, NULL);
    PS_ASSERT_PTR_NON_NULL(fringes, NULL);
    PS_ASSERT_PTR_NON_NULL(fringeImages, NULL);
    PS_ASSERT_PTR_NON_NULL(fringeStats, NULL);
    PS_ASSERT_INT_EQUAL(fringeImages->n, fringeStats->n, NULL);
    PS_ASSERT_INT_POSITIVE(nIter, NULL);

    // measure the fringe stats for the science frame and solve for the scales
    pmFringeStats *scienceStats = pmFringeStatsMeasure(fringes, readout, maskVal);

    if (psTraceGetLevel(__func__) > 9) {
        for (int i = 0; i < fringes->nRequested; i++) {
            printf("%f", scienceStats->f->data.F32[i]);
            for (int j = 0; j < fringeStats->n; j++) {
                pmFringeStats *fringe = fringeStats->data[j];
                printf("\t%f", fringe->f->data.F32[i]);
            }
            printf("\n");
        }
    }

    pmFringeScale *scale = pmFringeScaleMeasure(scienceStats, fringeStats, rej, nIter, keepFrac);
    psFree(scienceStats);

    psTrace(__func__, 7, "Fringe solution:\n");
    for (int i = 0; i < fringeImages->n + 1; i++) {
        psTrace(__func__, 7, "%d: %f %f\n", i, scale->coeff->data.F32[i],
                scale->coeffErr->data.F32[i]);
    }

    // build the fringe correction image
    // XXX we could save data space by making the first image the output image
    psImage *sumFringe = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
    //psBinaryOp(sumFringe, sumFringe, "+", psScalarAlloc(scale->coeff->data.F32[0], PS_TYPE_F32));
    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;
}

