#include <stdio.h>
#include <string.h>
#include "pslib.h"
#include "papmodule.h"
#include "papStuff.h"
#include "papFocalPlane.h"
#include "pmFPARead.h"

// NOTE: Need to deal with header inheritance

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// File-static functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

// Read a FITS extension into a chip
static bool readExtension(psArray **images, // Array of images in which to read pixels
			  psMetadata **header, // Metadata in which to read header
			  psFits *fits,	// The FITS file from which to read
			  const char *extName // The name of the extension
    )
{
    psTrace(__func__, 7, "Moving to extension %s...\n", extName);
    if (strncmp(extName, "PHU", 3) == 0) {
	if (!psFitsMoveExtNum(fits, 0, false)) {
	    psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n");
	    return false;
	}
    } else if (! psFitsMoveExtName(fits, extName)) {
	psError(PS_ERR_IO, false, "Unable to find extension %s in FITS file!\n", extName);
	return false;
    }
    psTrace(__func__, 7, "Reading header....\n");
    *header = psFitsReadHeader(*header, fits);
    if (! *header) {
	psError(PS_ERR_IO, false, "Unable to read FITS header!\n");
	return false;
    }

    psTrace(__func__, 7, "Checking NAXIS....\n");
    bool mdStatus = false;
    int nAxis = psMetadataLookupS32(&mdStatus, *header, "NAXIS");
    if (!mdStatus) {
	psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header of extension %s!\n",
		 extName);
    }
    if (nAxis != 2 && nAxis != 3) {
	psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into a single image "
		 "anyway.\n");
    }
    psTrace(__func__, 9, "NAXIS = %d\n", nAxis);

    int numPlanes = 1;			// Number of planes
    if (nAxis == 3) {
	numPlanes = psMetadataLookupS32(&mdStatus, *header, "NAXIS3");
	if (!mdStatus) {
	    psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n");
	    return false;
	}
    }
    psRegion region = {0, 0, 0, 0};	// Region to read is everything

    // Read each plane into the array
    psTrace(__func__, 7, "Reading %d planes into array....\n", numPlanes);
    *images = psArrayAlloc(numPlanes); // Array of images
    for (int i = 0; i < numPlanes; i++) {
	psTrace(__func__, 9, "Reading plane %d\n", i);
	psMemCheckCorruption(true);
	(*images)->data[i] = psFitsReadImage(NULL, fits, region, i);
	psTrace(__func__, 10, "Done\n");
        if (! (*images)->data[i]) {
	    psError(PS_ERR_IO, false, "Unable to read image for extension %s, plane %d\n", extName, i);
	    return false;
	}
    }

    return true;
}

// Read a list of regions (e.g., [x0:x1,y0:y1];[x2:x3,y2:y3]...), extract subimages and put them on a list
static int readMultipleRegions(psList *list, // List to which to add images
			       const char *string, // String from which to read regions.
			       psImage *image // Image from which to get the region
    )
{
    int numFound = 0;			// Number of regions found
    char *stringCopy = psStringCopy(string); // Copy of the string (to preserve it intact)
    char *startOfRegion = stringCopy; // The string containing all the regions yet to be parsed
    while (startOfRegion) {
	psRegion region = psRegionFromString(startOfRegion); // The region
	if (isnan(region.x0) || isnan(region.x1) || isnan(region.y0) || isnan(region.y1)) {
	    psLogMsg(__func__, PS_LOG_WARN, "Unable to read region (%s) --- ignored.\nOriginal string was "
		     "%s.", startOfRegion, string);
	} else {
	    psTrace(__func__, 9, "region found: %f,%f --> %f,%f\n", region.x0, region.y0, region.x1,
		    region.y1);
	    numFound++;
	    // Convert from FITS standard to PS standard
	    region.x0 -= 1;
	    region.y0 -= 1;
	    // We don't touch the x1 and y1 values because they aren't included by psImageSubset.
	    psImage *subImage = psImageSubset(image, region);
	    psListAdd(list, PS_LIST_TAIL, subImage);
	}

	startOfRegion = strchr(startOfRegion, ';');
	if (startOfRegion) {
	    // Eat the semi-colon
	    startOfRegion += 1;
	}
    }
    psFree(stringCopy);

    return numFound;
}


// Return a list of portions
static psList *portions(psImage *image,	// The image to portion up
			const char *methodValue, // How to determine the portions
			psMetadata *header // FITS header from the file
    )
{
    psList *subImages = psListAlloc(NULL); // The list of subimages, to be returned

    psTrace(__func__, 7, "Parsing %s to get subimages....\n", methodValue);
    // Parse the methodValue string into the METHOD:VALUE (e.g., HEADER:TRIMSEC, or VALUE:[1:100,1:100])
    char *value = strchr(methodValue, ':') + 1; // Find the first colon, which is the split character
    if (! value) {
	psError(PS_ERR_IO, true, "Unable to separate string (%s) into METHOD:VALUE\n", methodValue);
	return false;
    }
    psTrace(__func__, 7, "Value is %s\n", value);

    if (strncmp(methodValue, "HEADER:", 7) == 0 || strncmp(methodValue, "HD:", 3) == 0) {
	// Read the portions from the header
	psList *keywords = papSplit(value, ",;"); // List of headers to look at
	psListIterator *keywordsIter = psListIteratorAlloc(keywords, PS_LIST_HEAD, false); // Iterator
	psString keyword = NULL;	// The header keyword
	bool mdStatus = false;		// Status for MD lookup
	while (keyword = psListGetAndIncrement(keywordsIter)) {
	    psString regionString = psMetadataLookupString(&mdStatus, header, keyword);
	    if (! mdStatus) {
		psLogMsg(__func__, PS_LOG_WARN, "Unable to find header %s to determine image region --- "
			 "ignored.\n", keyword);
	    } else {
		// We have the region: [x0:x1,y0:y1]
		psRegion region = psRegionFromString(regionString);
		// Convert from FITS standard to PS standard
		region.x0 -= 1;
		region.y0 -= 1;
		// We don't touch the x1 and y1 values because they aren't included by psImageSubset.
		psImage *subImage = psImageSubset(image, region);
		psListAdd(subImages, PS_LIST_TAIL, subImage);
	    }
	}
    } else if (strncmp(methodValue, "VALUE:", 6) == 0 || strncmp(methodValue, "VAL:", 4) == 0) {
	// Read the portions right here
	if (readMultipleRegions(subImages, value, image) == 0) {
	    psLogMsg(__func__, PS_LOG_WARN, "No value specified: %s --- ignored.\n", methodValue);
	}
    } else {
	// Try to interpret the methodValue as a "value" only (ie, no "HEADER:" or "VALUE:")
	// Otherwise, we give up
	if (readMultipleRegions(subImages, methodValue, image) == 0) {
	    psLogMsg(__func__, PS_LOG_WARN, "Can't parse %s into a region --- ignored.\n", methodValue);
	}
    }

    psTrace(__func__, 7, "%d regions found.\n", subImages->size);

    return subImages;
}

// Portion out an image into the cell
static bool portionCell(papCell *cell,	// The cell that gets its bits
			psArray *images, // Array of images from which the readouts are assigned
			psMetadata *header // FITS header from the file
    )
{
    bool mdStatus = false;		// Status of MD lookup

    // How do we get the trim section?
    psString trimsecString = psMetadataLookupString(&mdStatus, cell->values, "CELL.TRIMSEC");
    if (! mdStatus) {
	psError(PS_ERR_IO, false, "Unable to find CELL.TRIMSEC specified --- unable to read cell.\n");
	return false;
    }

    // How do we get the bias section(s)?
    psString biassecString = psMetadataLookupString(&mdStatus, cell->values, "CELL.BIASSEC");
    if (! mdStatus) {
	psError(PS_ERR_IO, false, "Unable to find CELL.BIASSEC specified --- unable to read cell.\n");
	return false;
    }

    // Iterate over each of the image planes
    psArray *readouts = psArrayAlloc(images->n);
    for (int i = 0; i < images->n; i++) {
	// Get the TRIMSEC region
	psList *trimsecs = portions(images->data[i], trimsecString, header);
	if (trimsecs->size == 0) {
	    psError(PS_ERR_IO, true, "No TRIMSEC region found for image.  Specified: %s\n", trimsecString);
	    readouts->data[i] = NULL;
	} else {
	    // There can be only ONE!
	    if (trimsecs->size > 1) {
		psLogMsg(__func__, PS_LOG_WARN, "Multiple TRIMSECs were specified --- using only the first."
			 "\n");
	    }
	    psImage *trimsec = psListGet(trimsecs, PS_LIST_HEAD);
	    psFree(trimsecs);

	    // Get the BIASSEC regions
	    psList *biassecs = portions(images->data[i], biassecString, header);
	    if (biassecs->size == 0) {
		psLogMsg(__func__, PS_LOG_WARN, "No BIASSEC region found for image.  Specified: %s\n",
			 biassecString);
	    }

	    papReadout *readout = papReadoutAlloc(cell, i, trimsec, biassecs, 0,0,0,0,1,1);
	    readouts->data[i] = readout;
	}

    }

    // Stuff the array of readouts into the cell
    cell->readouts = readouts;

    return true;

}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Public functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool pmFPARead(papFPA *fpa,		// FPA to read into
	       psFits *fits		// FITS file from which to read
    )
{
    psArray *pixels = NULL;		// Current pixels from FITS file (array of images)
    psMetadata *header = NULL;		// Current header from FITS file

    psTrace(__func__, 1, "Working on FPA...\n");
    if (fpa->extname) {
	psTrace(__func__, 3, "Reading pixels from extension %s into FPA.\n", fpa->extname);
	if (! readExtension(&fpa->pixels, &fpa->header, fits, fpa->extname)) {
	    psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", fpa->extname);
	    return false;
	}
	pixels = fpa->pixels;
	header = fpa->header;
    }

    psArray *chips = fpa->chips;	// Array of chips
    // Iterate over the FPA
    for (int i = 0; i < chips->n; i++) {
	papChip *chip = chips->data[i]; // The chip
	psTrace(__func__, 2, "Working on chip %d...\n", i);
	if (chip->extname) {
	    psTrace(__func__, 3, "Reading pixels from extension %s into chip %d.\n", chip->extname, i);
	    if (! readExtension(&chip->pixels, &chip->header, fits, chip->extname)) {
		psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", chip->extname);
		return false;
	    }
	    pixels = chip->pixels;
	    header = chip->header;
	}
	// Iterate over the chip
	psArray *cells = chip->cells;	// Array of cells
	for (int j = 0; j < cells->n; j++) {
	    papCell *cell = cells->data[j]; // The cell
	    psTrace(__func__, 3, "Working on cell %d...\n", j);
	    if (cell->extname) {
		psTrace(__func__, 5, "Reading pixels from extension %s into cell %d.\n", cell->extname, j);
		if (! readExtension(&cell->pixels, &cell->header, fits, cell->extname)) {
		    psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", cell->extname);
		    return false;
		}
		pixels = cell->pixels;
		header = cell->header;
	    }

	    psTrace(__func__, 5, "Allocating readouts and extracting overscans etc for chip %d cell %d...\n",
		    i, j);
	    portionCell(cell, pixels, header);
	}
    }

    return true;
}
