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

#include <stdio.h>
#include <strings.h>
#include <assert.h>
#include <string.h>
#include <pslib.h>

#include "pmFPA.h"
#include "pmFPAFlags.h"
#include "pmHDU.h"
#include "pmHDUUtils.h"
#include "pmConcepts.h"
#include "pmFPAHeader.h"

#include "pmFPARead.h"

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Definitions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

// Specify what to read
typedef enum {
    FPA_READ_TYPE_IMAGE,                // Read image
    FPA_READ_TYPE_MASK,                 // Read mask
    FPA_READ_TYPE_WEIGHT,               // Read weight map
    FPA_READ_TYPE_HEADER                // Read header
} fpaReadType;

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

// Carve a readout from the image pixels
static bool readoutCarve(pmReadout *readout, // Readout to be carved up
                         psImage *image, // Image that will be carved
                         const psRegion *trimsec, // Trim section
                         const psList *biassecs, // Bias sections
			 fpaReadType type
                        )
{
    assert(readout);
    assert(image);
    assert(trimsec);
    assert(biassecs);

    // The image corresponding to the trim region
    if (psRegionIsNaN(*trimsec)) {
        psString regionString = psRegionToString(*trimsec);
        psError(PS_ERR_UNKNOWN, true, "Invalid trim section: %s\n", regionString);
        psFree(regionString);
        psFree(readout);
        return false;
    }
    psRegion region = psRegionSet(PS_MAX(trimsec->x0 - readout->col0, 0), // x0
                                  PS_MIN(trimsec->x1 - readout->col0, image->numCols), // x1
                                  PS_MAX(trimsec->y0 - readout->row0, 0), // y0
                                  PS_MIN(trimsec->y1 - readout->row0, image->numRows) // y1
                                 );

    // place the image subset in the appropriate target location, freeing if needed
    // XXX why psMemIncrRefCounter on psImageSubset??
    switch (type) {
      case FPA_READ_TYPE_IMAGE:
	if (readout->image) {
	    psFree (readout->image);
	}
	readout->image = psMemIncrRefCounter(psImageSubset(image, region));	
	break;
      case FPA_READ_TYPE_MASK:
	if (readout->mask) {
	    psFree (readout->mask);
	}
	readout->mask = psMemIncrRefCounter(psImageSubset(image, region));	
	break;
      case FPA_READ_TYPE_WEIGHT:
	if (readout->weight) {
	    psFree (readout->weight);
	}
	readout->weight = psMemIncrRefCounter(psImageSubset(image, region));	
	break;
      default:
	psAbort("Unknown read type: %x\n", type);
    }
    
    // Get the list of overscans
    // XXX should this step only be performed for IMAGE, not MASK and WEIGHT types?
    // XXX that would allow us to overlay a MASK and WEIGHT which have been trimmed...
    if (readout->bias->n != 0) {
        // Make way!
        psFree(readout->bias);
        readout->bias = psListAlloc(NULL);
    }
    psListIterator *iter = psListIteratorAlloc((psList*)biassecs, PS_LIST_HEAD, false); // Iterator
    psRegion *biassec = NULL;       // A BIASSEC region from the list
    while ((biassec = psListGetAndIncrement(iter))) {
        if (psRegionIsNaN(*biassec)) {
            psString regionString = psRegionToString(*biassec);
            psError(PS_ERR_IO, true, "Invalid bias section: %s\n", regionString);
            psFree(regionString);
            psFree(readout);
            return false;
        }
        psRegion region = psRegionSet(PS_MAX(biassec->x0 - readout->col0, 0), // x0
                                      PS_MIN(biassec->x1 - readout->col0, image->numCols), // x1
                                      PS_MAX(biassec->y0 - readout->row0, 0), // y0
                                      PS_MIN(biassec->y1 - readout->row0, image->numRows) // y1
                                     );
        psImage *overscan = psMemIncrRefCounter(psImageSubset(image, region));
        psListAdd(readout->bias, PS_LIST_TAIL, overscan);
        psFree(overscan);
    }
    psFree(iter);

    return true;
}

// Read a component of a readout
psImage *readoutReadComponent(psImage *image, // Image into which to read
                              psFits *fits, // FITS file from which to read
                              const psRegion *region, // Region to read
                              int readdir, // Read direction (1=rows, 2=cols)
                              int min,  // Minimum row/col number to read
                              int max,   // Maximum row/col number to read
                              int z,     // Image plane to read
                              float bad // Bad value
                             )
{
    assert(fits);
    assert(region);

    bool resize = false;                // Do we need to resize the image once read?
    psRegion toRead = psRegionSet(region->x0, region->x1, region->y0, region->y1); // Region to read
    psRegion fullSize = toRead;         // Full sized region
    if (readdir == 1) {
        if (toRead.y0 <= min) {
            toRead.y0 = min;
        } else {
            fullSize.y0 = min;
            resize = true;
        }
        if (toRead.y1 >= max) {
            toRead.y1 = max;
        } else {
            fullSize.y1 = max;
            resize = false;
        }
    } else if (readdir == 2) {
        if (toRead.x0 <= min) {
            toRead.x0 = min;
        } else {
            fullSize.x0 = min;
            resize = true;
        }
        if (toRead.x1 >= max) {
            toRead.x1 = max;
        } else {
            fullSize.x1 = max;
            resize = false;
        }
    } else {
        psAbort("Read direction can only be 1 (rows) or 2 (cols).\n");
    }

    psTrace("psModules.camera", 5, "Reading section [%.0f:%.0f,%.0f:%.0f]\n",
            toRead.x0, toRead.x1, toRead.y0, toRead.y1);
    image = psFitsReadImageBuffer(image, fits, toRead, z); // Desired pixels
    psTrace("psModules.camera", 7, "Image is %dx%d\n", image->numCols, image->numRows);

    // XXX: We only support F32 for now
    if (image->type.type != PS_TYPE_F32) {
        psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
        psFree(image);
        image = temp;
    }

    if (resize) {
        // For some reason, the region of interest is smaller than the number of pixels we want.
        psTrace("psModules.camera", 5, "Resizing image to %.0fx%.0f\n",
                fullSize.x1 - fullSize.x0, fullSize.y1 - fullSize.y0);
        psImage *temp = psImageAlloc(fullSize.x1 - fullSize.x0, fullSize.y1 - fullSize.y0, image->type.type);
        psImageInit(temp, bad);
        psImageOverlaySection(temp, image, toRead.x0, toRead.y0, "=");
        psFree(image);
        image = temp;
    }

    return image;
}


// Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadWeight
// Does most of the work for the reading --- reads the HDU, and portions the HDU into readouts.
static bool cellRead(pmCell *cell,      // Cell into which to read
                     psFits *fits,      // FITS file from which to read
                     psDB *db,          // Database handle, for concepts ingest
                     fpaReadType type   // Type to read
                    )
{
    assert(cell);
    assert(fits);

    pmHDU *hdu = pmHDUFromCell(cell);   // The HDU
    if (!hdu) {
        return true;                    // We read everything we could
    }

    // check if we have read the desired data, read it if needed
    bool (*hduReadFunc)(pmHDU*, psFits*) = NULL; // Function to use to read the HDU
    void *dataPointer = NULL; 		// pointer to location of desired data
    switch (type) {
      case FPA_READ_TYPE_IMAGE:
        hduReadFunc = pmHDURead;
        dataPointer = hdu->images;
        break;
      case FPA_READ_TYPE_HEADER:
        hduReadFunc = pmHDUReadHeader;
        dataPointer = hdu->header;
        break;
      case FPA_READ_TYPE_MASK:
        hduReadFunc = pmHDUReadMask;
        dataPointer = hdu->masks;
        break;
      case FPA_READ_TYPE_WEIGHT:
        hduReadFunc = pmHDUReadWeight;
        dataPointer = hdu->weights;
        break;
      default:
        psAbort("Unknown read type: %x\n", type);
    }

    // do we have the data we want (image, header, or etc).
    if (!dataPointer) {
	// attempt to read in the desired data
	if (!hduReadFunc(hdu, fits)) {
	    psError(PS_ERR_UNKNOWN, false, "Unable to read HDU for cell.\n");
	    return false;
	}
    }

    // load in the concept information for this cell
    if (!pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, true, NULL)) {
        //psError(PS_ERR_UNKNOWN, false, "Failed to read concepts for cell");
        //return false;
        psWarning("Difficulty reading concepts for cell; attempting to proceed.");
    }

    // skip the image arrays completely for the header-only files
    if (type == FPA_READ_TYPE_HEADER) {
	pmCellSetDataStatus(cell, true);
	return true;
    }

    // set up pointers for the different possible image arrays
    psArray *imageArray = NULL;	// Array of images in the HDU
    psElemType imageType = PS_TYPE_F32; // Expected type for image
    switch (type) {
      case FPA_READ_TYPE_IMAGE:
        imageArray = hdu->images;
        imageType = PS_TYPE_F32;
        break;
      case FPA_READ_TYPE_MASK:
	imageArray = hdu->masks;
        imageType = PS_TYPE_MASK;
        break;
      case FPA_READ_TYPE_WEIGHT:
        imageArray = hdu->weights;
        imageType = PS_TYPE_F32;
        break;
      default:
        psAbort("Unknown read type: %x\n", type);
    }

    // Having read the cell, we now have to cut it up
    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC");
    psList *biassecs = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC");
    if (psRegionIsNaN(*trimsec)) {
        psError(PS_ERR_IO, false, "CELL.TRIMSEC is not set --- can't read cell.\n");
        return false;
    }

    // Iterate over each of the image planes, converting type if necessary, and extracting the bits that
    // matter (CELL.TRIMSEC, CELL.BIASSEC) into readouts with readoutCarve.
    for (int i = 0; i < imageArray->n; i++) {
        psImage *source = imageArray->data[i]; // Source image, from the i-th plane

        // Type conversion here to support the modules, which don't have multiple type support yet
        if (source->type.type != imageType) {
            psImage *temp = psImageCopy(NULL, source, imageType); // Temporary image
            psFree(imageArray->data[i]);
            imageArray->data[i] = temp;
            source = temp;
        }

        pmReadout *readout;             // Readout into which to read
        if (cell->readouts->n > i && cell->readouts->data[i]) {
            readout = psMemIncrRefCounter(cell->readouts->data[i]);
        } else {
            readout = pmReadoutAlloc(cell);
        }

        if (!readoutCarve(readout, source, trimsec, biassecs, type)) {
            psError(PS_ERR_UNEXPECTED_NULL, false,
                    "Unable to carve readout into image and bias sections for %d the plane.", i);
            return NULL;
        }
        psFree(readout);                // Drop reference
    }

    pmCellSetDataStatus(cell, true);
    return true;
}


// Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadWeight
// Iterates over component cells, reading each
static bool chipRead(pmChip *chip,      // Chip into which to read
                     psFits *fits,      // FITS file from which to read
                     psDB *db,          // Database handle, for concepts ingest
                     fpaReadType type   // Type to read
                    )
{
    assert(chip);
    assert(fits);

    bool success = false;               // Were we able to read at least one HDU?
    psArray *cells = chip->cells;       // Array of cells
    for (int i = 0; i < cells->n; i++) {
        pmCell *cell = cells->data[i];  // The cell of interest
        success |= cellRead(cell, fits, db, type);
    }
    if (success) {
        if (!pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_HEADER, true, true, NULL)) {
            psError(PS_ERR_IO, false, "Failed to read concepts for FPA.\n");
            return false;
        }
        // XXX probably could just use chip->data_exists
        pmChipSetDataStatus(chip, true);
    }

    return success;
}


// Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadWeight
// Iterates over component chips, reading each
static bool fpaRead(pmFPA *fpa,         // FPA into which to read
                    psFits *fits,       // FITS file from which to read
                    psDB *db,           // Database handle, for concepts ingest
                    fpaReadType type    // Type to read
                   )
{
    assert(fpa);
    assert(fits);

    bool success = false;               // Were we able to read at least one HDU?
    psArray *chips = fpa->chips;        // Array of chips
    for (int i = 0; i < chips->n; i++) {
        pmChip *chip = chips->data[i];  // The cell of interest
        success |= chipRead(chip, fits, db, type);
    }
    if (success) {
        if (!pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_HEADER, true, NULL)) {
            psError(PS_ERR_IO, false, "Failed to read concepts for FPA.\n");
            return false;
        }
    } else {
        psError(PS_ERR_UNKNOWN, false, "Unable to read any chips in FPA");
    }

    return success;
}

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


//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Reading images
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool pmReadoutReadNext(pmReadout *readout, psFits *fits, int z, int numScans)
{
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);
    PS_ASSERT_INT_NONNEGATIVE(z, false);
    PS_ASSERT_INT_NONNEGATIVE(numScans, false);

    // Get the HDU and read the header
    pmCell *cell = readout->parent;     // The parent cell

    pmHDU *hdu = pmHDUFromCell(cell);   // The HDU
    if (!hdu || hdu->blankPHU) {
        return false;
    }

    if (!pmCellReadHeader(cell, fits)) {
        psError(PS_ERR_IO, false, "Unable to read header for cell!\n");
        return false;
    }

    // Make sure we have the information we need
    if (!pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS |
                            PM_CONCEPT_SOURCE_DEFAULTS, true, NULL)) {
        psError(PS_ERR_IO, false, "Failed to read concepts for cell.\n");
        return false;
    }

    // Get the trim and bias sections
    bool mdok = true;                   // Status of MD lookup
    psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim sections
    if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) {
        psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set.\n");
        return false;
    }
    psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections
    if (!mdok || !biassecs) {
        psError(PS_ERR_IO, true, "CELL.BIASSEC is not set.\n");
        return false;
    }
    int readdir = psMetadataLookupS32(&mdok, cell->concepts, "CELL.READDIR"); // Read direction
    if (!mdok || readdir == 0 || (readdir != 1 && readdir != 2)) {
        psError(PS_ERR_IO, true, "CELL.READDIR is not set to -1 or +1.\n");
        return false;
    }
    float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
    if (!mdok) {
        psLogMsg(__func__, PS_LOG_WARN, "CELL.BAD is not set --- assuming zero.\n");
        bad = 0.0;
    }

    // Check the third dimension
    int naxis = psMetadataLookupS32(&mdok, hdu->header, "NAXIS"); // The number of axes
    if (!mdok) {
        psError(PS_ERR_IO, true, "Unable to find NAXIS in header for extension %s\n", hdu->extname);
        return false;
    }
    if (naxis == 0) {
        // No pixels to read, as for a PHU.
        return false;
    }
    if (naxis < 2 || naxis > 3) {
        psError(PS_ERR_IO, true, "NAXIS in header of extension %s (= %d) is not valid.\n",
                hdu->extname, naxis);
        return false;
    }
    int naxis3 = 1;                     // The number of image planes
    if (naxis == 3) {
        naxis3 = psMetadataLookupS32(&mdok, hdu->header, "NAXIS3");
        if (!mdok) {
            psError(PS_ERR_IO, true, "Unable to find NAXIS3 in header for extension %s\n", hdu->extname);
            return false;
        }
    }
    if (z >= naxis3) {
        // Nothing to see here.  Move along.
        return false;
    }

    // Get the size of the image plane
    int naxis1 = psMetadataLookupS32(&mdok, hdu->header, "NAXIS1"); // The number of columns
    if (!mdok) {
        psError(PS_ERR_IO, true, "Unable to find NAXIS1 in header for extension %s\n", hdu->extname);
        return false;
    }
    int naxis2 = psMetadataLookupS32(&mdok, hdu->header, "NAXIS2"); // The number of columns
    if (!mdok) {
        psError(PS_ERR_IO, true, "Unable to find NAXIS2 in header for extension %s\n", hdu->extname);
        return false;
    }

    // rationalize trimsec against naxis1, naxis2
    // valid range for trimsec is 1-Nx,1-Ny
    // if (trimsec->x0 == 0) trimsec->x0 = 1;
    if (trimsec->x1 <  1)
        trimsec->x1 = naxis1 + trimsec->x1;
    // if (trimsec->y0 == 0) trimsec->y0 = 1;
    if (trimsec->y1 <  1)
        trimsec->y1 = naxis2 + trimsec->y1;

    int maxSize;                        // Number of cols,rows in image
    if (readdir == 1) {
        maxSize = PS_MIN(naxis2, trimsec->y1);
    } else {
        maxSize = PS_MIN(naxis1, trimsec->x1);
    }

    // Read the FITS image
    int offset = 0;                     // The offset from the start of the image to where we are now
    if (readout->image) {
        if (readdir == 1) {
            // Reading rows
            if (numScans == 0) {
                numScans = naxis2 - readout->row0;
            }
            readout->row0 += readout->image->numRows;
            offset = readout->row0;
        } else {
            // Reading columns
            if (numScans == 0) {
                numScans = naxis1 - readout->col0;
            }
            readout->col0 += readout->image->numCols;
            offset = readout->col0;
        }
    } else {
        if (numScans == 0) {
            if (readdir == 1) {
                numScans = naxis2;
            } else
                numScans = naxis1;
        }
    }

    // Blow away existing data.
    // Do this before returning, so that we're not returning data from a previous read
    psFree(readout->image);
    readout->image = NULL;

    while (readout->bias->n > 0) {
        psListRemove(readout->bias, PS_LIST_HEAD);
    }

    if (offset >= maxSize) {
        // We've read everything there is
        psTrace("psModules.camera", 7, "Read everything.\n");
        psFree(readout->image);
        return false;
    }
    int upper = PS_MIN(offset + numScans, maxSize); // The upper limit for the pixel read
    psTrace("psModules.camera", 7, "offset=%d, upper = %d, image is %dx%d, trimsec [%.0f:%.0f,%.0f:%.0f]\n",
            offset, upper, naxis1, naxis2, trimsec->x0, trimsec->x1, trimsec->y0, trimsec->y1);

    // Get the new the trim section
    readout->image = readoutReadComponent(readout->image, fits, trimsec, readdir,
                                          offset, upper, z, bad); // The image

    // Get the new bias sections
    psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator for BIASSEC
    psRegion *biassec = NULL;           // Bias section from iteration
    while ((biassec = psListGetAndIncrement(biassecsIter))) {
        psImage *bias = readoutReadComponent(NULL, fits, biassec, readdir, offset, upper, z, bad); // The bias
        psListAdd(readout->bias, PS_LIST_TAIL, bias);
        psFree(bias);                   // Drop reference
    }
    psFree(biassecsIter);

    return true;
}


bool pmCellRead(pmCell *cell, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(cell, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return cellRead(cell, fits, db, FPA_READ_TYPE_IMAGE);
}

bool pmChipRead(pmChip *chip, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return chipRead(chip, fits, db, FPA_READ_TYPE_IMAGE);
}

bool pmFPARead(pmFPA *fpa, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return fpaRead(fpa, fits, db, FPA_READ_TYPE_IMAGE);
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Reading the mask
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool pmCellReadMask(pmCell *cell, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(cell, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return cellRead(cell, fits, db, FPA_READ_TYPE_MASK);
}

bool pmChipReadMask(pmChip *chip, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return chipRead(chip, fits, db, FPA_READ_TYPE_MASK);
}

bool pmFPAReadMask(pmFPA *fpa, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return fpaRead(fpa, fits, db, FPA_READ_TYPE_MASK);
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Reading the weight map
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool pmCellReadWeight(pmCell *cell, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(cell, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return cellRead(cell, fits, db, FPA_READ_TYPE_WEIGHT);
}

bool pmChipReadWeight(pmChip *chip, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return chipRead(chip, fits, db, FPA_READ_TYPE_WEIGHT);
}

bool pmFPAReadWeight(pmFPA *fpa, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return fpaRead(fpa, fits, db, FPA_READ_TYPE_WEIGHT);
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Reading the image header
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool pmCellReadHeaderSet(pmCell *cell, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(cell, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return cellRead(cell, fits, db, FPA_READ_TYPE_HEADER);
}

bool pmChipReadHeaderSet(pmChip *chip, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return chipRead(chip, fits, db, FPA_READ_TYPE_HEADER);
}

bool pmFPAReadHeaderSet(pmFPA *fpa, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);
    PS_ASSERT_PTR_NON_NULL(fits, false);

    return fpaRead(fpa, fits, db, FPA_READ_TYPE_HEADER);
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Reading FITS tables
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

int pmCellReadTable(pmCell *cell, psFits *fits, const char *name)
{
    PS_ASSERT_PTR_NON_NULL(cell, 0);
    PS_ASSERT_PTR_NON_NULL(fits, 0);
    PS_ASSERT_STRING_NON_EMPTY(name, 0);

    const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip
    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
    psString extname = NULL;            // Extension name
    psStringAppend(&extname, "%s_%s_%s", name, chipName, cellName);

    // XXX Could do a table lookup from the camera format, in case the input file isn't laid out with
    // NAME_CHIP_CELL extension names --- use these as keys, and the value as the proper extension name.
    // Allow interpolation of concepts, e.g., "{CHIP.NAME}" --> "ccd13".

    if (!psFitsMoveExtName(fits, extname)) {
        psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname);
        psFree(extname);
        return 0;
    }

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

    psString headerName = NULL;         // Name for header
    psStringAppend(&headerName, "%s.HEADER", name);
    if (!psMetadataAdd(cell->analysis, PS_LIST_TAIL, headerName, PS_DATA_METADATA,
                       "FITS table header", header)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to add header from extension %s to analysis metadata "
                "for chip %s, cell %s\n", extname, chipName, cellName);
        psFree(headerName);
        psFree(header);
        psFree(extname);
        return 0;
    }
    psFree(headerName);
    psFree(header);

    psArray *table = psFitsReadTable(fits); // The table
    if (!psMetadataAdd(cell->analysis, PS_LIST_TAIL, name, PS_DATA_ARRAY, "FITS table", table)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to add table from extension %s to analysis metadata "
                "for chip %s, cell %s\n", extname, chipName, cellName);
        psFree(table);
        psFree(extname);
        return 0;
    }

    psFree(extname);
    psFree(table);                      // Dropping reference

    return 1;
}


int pmChipReadTable(pmChip *chip, psFits *fits, const char *name)
{
    PS_ASSERT_PTR_NON_NULL(chip, 0);
    PS_ASSERT_PTR_NON_NULL(fits, 0);
    PS_ASSERT_STRING_NON_EMPTY(name, 0);

    int numRead = 0;                    // Number of reads
    psArray *cells = chip->cells;       // Array of cells
    for (int i = 0; i < cells->n; i++) {
        pmCell *cell = cells->data[i];  // Cell of interest
        numRead += pmCellReadTable(cell, fits, name);
    }

    return numRead;
}


int pmFPAReadTable(pmFPA *fpa, psFits *fits, const char *name)
{
    PS_ASSERT_PTR_NON_NULL(fpa, 0);
    PS_ASSERT_PTR_NON_NULL(fits, 0);
    PS_ASSERT_STRING_NON_EMPTY(name, 0);

    int numRead = 0;                    // Number of reads
    psArray *chips = fpa->chips;        // Array of chips
    for (int i = 0; i < chips->n; i++) {
        pmChip *chip = chips->data[i];  // Chip of interest
        numRead += pmChipReadTable(chip, fits, name);
    }

    return numRead;
}
