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

#include "psvideophot.h"

// XXX: get these from a recipe
#define DX 20
#define DY 10
#define SKIP 0
#define EDGE 3
#define THRESHOLD 100

#define ESCAPE(MESSAGE) { \
  psError(PS_ERR_UNKNOWN, false, MESSAGE); \
  psFree(view); \
  pmFPAfileFreeSetStrict(false); \
  return false; \
}

// XXX: move these to another file 
static double SubSum (psImage *image, double xy[2], double median);
static void Centering(psImage *image, long edge, double *pSum, double xy[], double median);

bool psvideophotLoop(pmConfig *config, psvideophotOptions *options)
{
    psMetadata *stats = NULL;           // Statistics to output
    float timePhot = 0;                 // Amount of time spent in photometry

    if (options->doStats) {
        psMetadata *statsMD = psMetadataAlloc();
        psMetadataAddS32(statsMD, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
    }

    bool status;                        // Status of MD lookup

    pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PSVIDEOPHOT.INPUT");
    
    if (!status) {
        psErrorStackPrint(stderr, "Can't find input data!\n");
        psvideophotCleanup(config, options);
        exit(PS_EXIT_PROG_ERROR);
    }


    pmConfigCamerasCull(config, NULL);
    pmConfigRecipesCull(config, "PSVIDEOPHOT,PPSTATS,PSPHOT,MASKS,PSASTRO,JPEG");

    pmFPAview *view = pmFPAviewAlloc(0);// View for level of interest
    pmHDU *lastHDU = NULL;              // Last HDU that was updated

    // files associated with the science image
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        ESCAPE("load failure for FPA");
    }

    pmChip *chip;                       // Chip from FPA
    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
        psLogMsg ("psvideophotLoop", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
        if (!chip->process || !chip->file_exists) {
            continue;
        }
        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
            ESCAPE("load failure for Chip");
        }

        double threshold = THRESHOLD;

        psTimerStart(TIMER_PHOT);
        pmCell *cell;                   // Cell from chip
        while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
            psLogMsg ("psvideophotLoop", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
            if (!cell->process || !cell->file_exists) {
                continue;
            }
            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                ESCAPE("load failure for Cell");
            }

            // Put version information into the header
            pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
            if (hdu && hdu != lastHDU) {
                psvideophotVersionHeader(hdu->header);
                lastHDU = hdu;
            }

            if (cell->readouts->n <2) {
	      continue; 
	    }

            psString cellName = psMetadataLookupStr(&status, cell->concepts, "CELL.NAME");
	    psWarning("Found video cell %s %d\n", cellName, view->cell);

            psStats *stats = psStatsAlloc(
                PS_STAT_SAMPLE_MEDIAN | 
                PS_STAT_SAMPLE_STDEV |
                PS_STAT_ROBUST_MEDIAN | 
                PS_STAT_ROBUST_STDEV |
                PS_STAT_MIN | PS_STAT_MAX | PS_STAT_SAMPLE_QUARTILE 
                | PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS
                );

            // process each of the readouts
            int frame = 0;
            pmReadout *readout;         // Readout from cell
            while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
	        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                    ESCAPE("load failure for Readout");
	        }
                if (++frame < SKIP) {
                    continue;
                }
                if (!readout->data_exists) {
                    continue;
                }

                psStatsInit(stats);

                if (!psImageStats(stats, readout->image, NULL, 0)) {
                    psError(PS_ERR_UNKNOWN, false, "unable to generate image stats for cell %d readout %d\n",
                        view->cell, view->readout);
                    ESCAPE("failed to generate stats");
                }

                double stddev = stats->sampleStdev;
                double median = stats->sampleMedian;
                if (stddev > threshold) {

                    double sum = 0;
                    // XXX: use a point data type
                    double xy[2];
                    Centering(readout->image, EDGE, &sum, xy, median);
                    double sum2;
                    sum2 = SubSum(readout->image, xy, median);

                    printf("%4d %9.2f %9.2f %8.2f %.1f %.1f %9.2f\n",
                            frame, sum, stddev, median, xy[0], xy[1], sum2);
                } else {
                    printf("%4d %9.2f %9.2f %8.2f\n",
                            frame, 0.0,   stddev, median);
                }
            }
        }

        // we perform photometry on the readouts of this chip in the output

        timePhot += psTimerClear(TIMER_PHOT);

        view->cell = -1;

        // Close chip
        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
            ESCAPE("save failure for Chip");
        }
    }
    psFree (stats);
    
    if (psTraceGetLevel("psvideophot") >= 3) {
        psvideophotFileCheck(config);
    }
    
    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
        ESCAPE("save failure for FPA");
    }
    psFree(view);

    // Dump configuration
    psString dump_file = psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
    if (dump_file) {
        if (!pmConfigDump(config, dump_file)) {
            ESCAPE("Unable to dump configuration.");
        }
    }

    return true;
}

// calculate the sum of pixel values inside the window
static double Sum(psImage *image, long edge, double median) {
    long	x, y;
    double	sum=0;
    for (y=edge; y<=image->numRows-edge+1; y++) {
        // TYPE!
        psF32 *row = image->data.F32[y];
        for (x=edge; x<= image->numCols - edge+1; x++) {
            sum = sum + row[x] - median;
        }
    }	
    return(sum);
}

// find the star
static void Centering(psImage *image, long edge, double *pSum, double xy[], double median) {
    long	x, y;
    double	sum=0, sxf=0, syf=0;
	
    sum=Sum(image, edge, median);

    for (y=edge; y<= image->numRows -edge+1; y++) {
        // TYPE!!
        psF32 *row = image->data.F32[y];
        for (x=edge; x<=image->numCols -edge+1; x++) {
            sxf=sxf+x*(row[x] - median);
            syf=syf+y*(row[x] - median);
        }
    }
    xy[0] = edge-1+sxf/sum;
    xy[1] = edge-2+syf/sum;
    *pSum = sum;
    //printf("%f %f\n", xc, yc);
    //return(0);
}

// calculate the sum of pixel values in an aperture around a point
double SubSum (psImage *image, double xy[2], double median)
{
    double	sum=0;
    long        x, y;

    for (y = xy[1] - DY; y<= xy[1] + DY; y++) {
        if (y < image->row0) {
            continue;
        }
        if (y >= image->numRows) {
            break;
        }
        psF32 *row = image->data.F32[y];
        for (x = xy[0] - DX; x<= xy[0] + DX; x++) {
            if (x < image->col0) {
                continue;
            }
            if (x >= image->numCols) {
                break;
            }
            sum += row[x] - median;
        }
    }	
    return(sum);
}
