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

#include "psvideophot.h"

// XXX: derive these from a recipe perhaps basing them on seeing
#define DX 10
#define DY 5
#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);
static psArray *readVideoTable(pmFPAfile *input, psString video_extension_name);

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

    psString class_id = psMetadataLookupStr(&status, config->arguments, "CLASS_ID");
    psString cell_id = psMetadataLookupStr(&status, config->arguments, "CELL_ID");


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

    psArray *output = NULL;

    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);
        psString chipName = psMetadataLookupStr(&status, chip->concepts, "CHIP.NAME");
        if (strcmp(chipName, class_id)) {
            continue;
        }
        if (!chip->process || !chip->file_exists) {
            ESCAPE("target chip has no data");
        }
        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
            ESCAPE("load failure for Chip");
        }
        (void) class_id;

        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);
            psString cellName = psMetadataLookupStr(&status, cell->concepts, "CELL.NAME");
            if (strcmp(cellName, cell_id)) {
                continue;
            }
            if (!cell->process || !cell->file_exists) {
                ESCAPE("target cell has no data");
            }
            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                ESCAPE("load failure for Cell");
            }

	    psLogMsg("psvideophotLoop", 1, "Found target cell %s %d. %ld readouts\n", cellName, view->cell, cell->readouts->n);

            // 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) {
              psLogMsg("psvideophotLoop", 1, "cell %s has only %ld readouts\n", cellName, cell->readouts->n);
              ESCAPE("invalid video Cell");
	    }

            // get the position of this cell on the chip 
            int cellX0 = psMetadataLookupS32(&status, cell->concepts, "CELL.X0");
            int cellY0 = psMetadataLookupS32(&status, cell->concepts, "CELL.Y0");
            int xParity = psMetadataLookupS32(&status, cell->concepts, "CELL.XPARITY");
            int yParity = psMetadataLookupS32(&status, cell->concepts, "CELL.YPARITY");
            int xBin = psMetadataLookupS32(&status, cell->concepts, "CELL.XBIN");
            int yBin = psMetadataLookupS32(&status, cell->concepts, "CELL.YBIN");

            // avoid errors about unused variables
            (void) xParity; (void) yParity; (void) cellX0; (void)cellY0;

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

            // go find the video extension data
            psString video_extension_name = psMetadataLookupStr(&status, config->arguments, "VIDEOEXTNAME");
            if (!video_extension_name) {
                psStringAppend(&video_extension_name, "video_table_%s", cellName);
            }
            psArray *video_table = readVideoTable(input, video_extension_name);
            if (!video_table) {
                psError(PS_ERR_UNKNOWN, false, "failed to find video extension:%s", video_extension_name);
                ESCAPE("invalid image file");
            }
            psFree(video_extension_name);
            // psMetadataAddArray(cell->analysis, PS_LIST_TAIL, "VIDEO_DATA", 0, "", video_table);

            output = (psArray *) psMemIncrRefCounter(video_table);

            if (cell->readouts->n != video_table->n) {
                psError(PS_ERR_UNKNOWN, true, "number of readouts %ld does not match length of video extension: %ld", 
                    cell->readouts->n, video_table->n);
                ESCAPE("invalid data");
            }

            // process each of the readouts
            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");
	        }
                int frame = view->readout;
                if (frame < SKIP - 1) {
                    continue;
                }
                if (!readout->data_exists) {
                    continue;
                }

                psMetadata *row = video_table->data[frame];
                if (!row) {
                    psError(PS_ERR_UNKNOWN, true, "video_table row for frame %d is null", frame);
                    // XXX: rather than fail, should we just skip this frame
                    ESCAPE(" ");
                }

                // extract the camera's measurements from the video table
                // We should probably get the column names from a configuration file so that we can
                // handle any changes more gracefully
                // Time at which frame was completely received from controller
                double frame_start = psMetadataLookupF64(&status, row, "frame_start_time");
                // Time at which frame was completely received from controller
                double frame_complete = psMetadataLookupF64(&status, row, "frame_complete_time");
                // star centriod, x axis
                double centroid_x = psMetadataLookupF32(&status, row, "centroid_x");
                // star centriod, y axis
                double centroid_y = psMetadataLookupF64(&status, row, "centroid_y");
                // star full width at half maximum
                double fwhm = psMetadataLookupF64(&status, row, "fwhm");
                // star full width at half maximum, x axis
                double fwhm_x = psMetadataLookupF64(&status, row, "fwhm_x");
                // star full width at half maximum, y axis
                double fwhm_y = psMetadataLookupF64(&status, row, "fwhm_y");
                // calculated sky level
                double sky = psMetadataLookupF64(&status, row, "sky");
                // flux from star
                double flux = psMetadataLookupF64(&status, row, "flux");
                // signal to noise ration of video frame
                double snr = psMetadataLookupF64(&status, row, "snr");
                // number of parallel row shits before integration (cleans)
                double rowpre = psMetadataLookupS32(&status, row, "rowpre");
                // number of cols skipped before video readout
                double cnpix1 = psMetadataLookupS32(&status, row, "cnpix1");
                // number of rows skipped before video readout
                double cnpix2 = psMetadataLookupS32(&status, row, "cnpix2");

                // avoid errors from unused variables
                (void) frame_start; (void) frame_complete;
                (void) rowpre; (void) cnpix1; (void) cnpix2; (void) fwhm_x; (void) fwhm_y; (void) fwhm;

                double star_chip_x = cellX0 + (centroid_x * xParity);
                double star_chip_y = cellY0 + (centroid_y * yParity);
                (void) star_chip_x; (void) star_chip_y;

                // measured position of star in readout coordinates
                double star_readout_x = (centroid_x - cnpix1) / xBin;
                double star_readout_y = (centroid_y - cnpix2) / yBin;

#ifdef notdef
                printf("%4d %9.2f %9.2f %9.2f %.1f %.1f %9.6f %6.3f %6.3f %6.3f\n",
                        frame, flux, snr, sky, star_readout_x, star_readout_y, frame_complete - frame_start,
                            fwhm, fwhm_x, fwhm_y);
#endif
                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);
                    // XXX: rather than fail, should we just skip this frame
                    ESCAPE("failed to generate stats");
                }

                double stddev = stats->sampleStdev;
                double median = stats->sampleMedian;
                double vpFlux = 0;
                double vpSum = 0;
                double xy[2] = {0, 0};
                if (stddev > threshold) {

                    Centering(readout->image, EDGE, &vpSum, xy, median);

                    if (isfinite(xy[0]) && isfinite(xy[1])) {
                        vpFlux = SubSum(readout->image, xy, median);
                    }

//                    printf("%4d %9.2f %9.2f %9.2f %.1f %.1f %9.2f\n",
//                             frame, vpFlux, vpFlux/median, median, xy[0], xy[1], stddev);
                } else {
//                    printf("%4d %9.2f %9.2f %8.2f\n",
//                            frame, 0.0,   stddev, median);
                }
                if (frame % 60 == 0) {
                    printf("frame  flux        sky       snr     X    Y     dt        fwhm   fwhm_x fwhm_y  vp_flux   vp_median vp_sum    vp_stddev vp_x vp_y\n");
                }
                printf("%4d   %9.2f %9.2f %9.2f %.1f %.1f %9.6f %6.3f %6.3f %6.3f   %9.2f %9.2f %9.2f %9.2f %.1f %.1f\n",
                        frame, flux, sky, snr, star_readout_x, star_readout_y, frame_complete - frame_start,
                            fwhm, fwhm_x, fwhm_y,
                            vpFlux, median, vpSum, stddev, xy[0], xy[1]);
                psMetadataAddF32(row, PS_LIST_TAIL, "vp_flux", 0, "vp flux", vpFlux);
                psMetadataAddF32(row, PS_LIST_TAIL, "vp_median", 0, "vp meidan", median);
                psMetadataAddF32(row, PS_LIST_TAIL, "vp_stddev", 0, "vp meidan", stddev);
                psMetadataAddF32(row, PS_LIST_TAIL, "vp_sum", 0, "vp sum of source", vpSum);
                psMetadataAddF32(row, PS_LIST_TAIL, "vp_x", 0, "centroid x from vp", xy[0]);
                psMetadataAddF32(row, PS_LIST_TAIL, "vp_y", 0, "centroid y from vp", xy[1]);
            }
            psLogMsg ("psvideophotLoop", 5, "Done processing video cell");
        }

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

    psString outputName = psMetadataLookupStr(&status, config->arguments, "OUTPUT");
    psString realName = pmConfigConvertFilename(outputName, config, true, true);
    if (!realName) {
        psError(PS_ERR_UNKNOWN, false, "failed to convert filename for %s", outputName);
        ESCAPE("cannot save output");
    }
    psFits *outFits = psFitsOpen(realName, "w");
    if (!outFits) {
        psError(PS_ERR_UNKNOWN, false, "failed to open output filename for %s", outputName);
        ESCAPE(" ");
    }

    psMetadata *outputHeader = lastHDU ? lastHDU->header : NULL;
    if (!psFitsWriteTable(outFits, outputHeader, output, "videophotemetry")) {
        psError(PS_ERR_UNKNOWN, false, "failed to write videophotometry table");
        ESCAPE("cannot write output");
    }
    psFitsClose(outFits);

    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-2+sxf/sum;
    xy[1] = edge-2+syf/sum;
    xy[0] = sxf/sum;
    xy[1] = 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);
}
static psArray *readVideoTable(pmFPAfile *input, psString video_extension_name)
{
    if (!psFitsMoveExtName(input->fits, video_extension_name)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to move to video extension: %s", video_extension_name);
        return NULL;
    }

    return psFitsReadTable(input->fits);
}
