# include "pslib.h"

typedef struct {
    double xMin;
    double yMin;
    double xMax;
    double yMax;
    double delta;
    double midValue;
} pmFringePoint;

stats *pmFringeStats (psImage *image, psArray *fringePoints, psMetadata *config) {

    psRegion region;

    double RADIUS = psMetadataLookupR32 (&status, myHeader, "FRINGE_SQUARE_RADIUS");

    psStats *stats = psStatsAlloc (PS_STATS_ROBUST_MEDIAN | PS_STATS_ROBUST_STDEV);

    // clip out the stars?

    psVector *delta = psVectorAlloc (fringePoints->n, PS_DATA_F32);

    for (i = 0; i < fringePoints->n; i++) {

	pmFringePoint *point = fringePoints->data[i];
	
	region = psRegionForSquare (point->xMin, point->yMin, RADIUS);
	region = psRegionForImage (image, region);
	subimage = psImageSubset (image, region);
	stats = psImageStats (stats, subimage, NULL, 0);
	psImageFree (subimage);

	min = stats->robustMedian;

	region = psRegionForSquare (point->xMax, point->yMax, RADIUS);
	region = psRegionForImage (image, region);
	subimage = psImageSubset (image, region);
	stats = psImageStats (stats, subimage, NULL, 0);
	psImageFree (subimage);

	max = stats->robustMedian;

	// check for NaN, ignore 
	point->delta = max - min;
	point->midValue = 0.5*(min + max);

	delta->data.F32[i] = point->delta;
    }

    stats = psVectorStats (stats, delta, NULL, NULL, 0);
    return (value);
}
