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

#include "pmFPA.h"
#include "pmFPARead.h"
#include "pmConcepts.h"

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


// Read a FITS extension into a chip
static bool readExtension(p_pmHDU *hdu, // Pixel data into which to read
                          psFits *fits  // The FITS file from which to read
                         )
{
    const char *extName = hdu->extname; // Extension name

    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");
    psMetadata *header = psFitsReadHeader(NULL, fits); // Header
    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);
    psArray *pixels = psArrayAlloc(numPlanes); // Array of images
    for (int i = 0; i < numPlanes; i++) {
        psTrace(__func__, 9, "Reading plane %d\n", i);
        psMemCheckCorruption(true);
        psImage *image = psFitsReadImage(NULL, fits, region, i);

        // XXX: Type conversion here to support the modules, which don't have multiple type support yet
        if (image->type.type != PS_TYPE_F32) {
            pixels->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);

            // XXX: Temporary fix for writing images, until psFits gets cleaned up
            psMetadataItem *bitpixItem = psMetadataLookup(header, "BITPIX");
            bitpixItem->data.S32 = -32;
            psMetadataRemove(header, 0, "BZERO");
            psMetadataRemove(header, 0, "BSCALE");

            psFree(image);
        } else {
            pixels->data[i] = image;
        }
        psTrace(__func__, 10, "Done\n");
        if (! pixels->data[i]) {
            psError(PS_ERR_IO, false, "Unable to read image for extension %s, plane %d\n", extName, i);
            return false;
        }
    }

    // Update the HDU with the new data
    hdu->images = pixels;
    hdu->header = header;

    // Mask and weight not set yet
    hdu->weights = NULL;
    hdu->masks = NULL;

    return true;
}

// Portion out an image into the cell
static bool generateReadouts(pmCell *cell, // The cell that gets its bits
                             p_pmHDU *hdu // Pixel data, containing image, mask, weights
                            )
{
    psArray *images = hdu->images;      // Array of images (each of which is a readout)

    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC");
    psList *biassecs = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC");

    // Iterate over each of the image planes
    for (int i = 0; i < images->n; i++) {
        psImage *image = images->data[i]; // The i-th plane
        pmReadout *readout = pmReadoutAlloc(cell);

        readout->image = psMemIncrRefCounter(psImageSubset(image, *trimsec)); // The image corresponding to the trim region

        // Get the list of overscans
        psListIterator *iter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator
        psRegion *biassec = NULL;       // A BIASSEC region from the list
        while ((biassec = psListGetAndIncrement(iter))) {
            psImage *overscan = psMemIncrRefCounter(psImageSubset(image, *biassec));
            psListAdd(readout->bias, PS_LIST_TAIL, overscan);
            psFree(overscan);
        }
        psFree(iter);

        readout->mask = NULL;
        readout->weight = NULL;

        psFree(readout);                // Drop reference
    }

    return true;
}

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

bool pmFPARead(pmFPA *fpa,              // FPA to read into
               psFits *fits,            // FITS file from which to read
               psMetadata *phu,         // Primary header
               psDB *db                 // Database handle, for concept ingest
              )
{
    p_pmHDU *hdu = NULL;        // Pixel data from FITS file

    // Read the PHU, if required
    if (! phu) {
        if (! psFitsMoveExtNum(fits, 0, false)) {
            psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n");
            return false;
        }
        psTrace(__func__, 7, "Reading PHU....\n");
        fpa->phu = psFitsReadHeader(NULL, fits); // Primary header
    } else if (! fpa->phu) {
        fpa->phu = psMemIncrRefCounter(phu);
    }

    // Read in....
    psTrace(__func__, 1, "Reading FPA...\n");
    if (fpa->hdu) {
        hdu = fpa->hdu;
        psTrace(__func__, 3, "Reading pixels from extension %s into FPA.\n", hdu->extname);
        if (! readExtension(hdu, fits)) {
            psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname);
            return false;
        }
    }
    pmConceptsReadFPA(fpa, db);

    psArray *chips = fpa->chips;        // Array of chips
    // Iterate over the FPA
    for (int i = 0; i < chips->n; i++) {
        pmChip *chip = chips->data[i]; // The chip

        // Only read chips marked to "read"
        if (! chip || ! chip->process || chip->exists) {
            psTrace(__func__, 2, "Ignoring chip %d...\n", i);
            continue;
        }
        psTrace(__func__, 2, "Reading in chip %d...\n", i);

        if (chip->hdu) {
            hdu = chip->hdu;
            psTrace(__func__, 3, "Reading pixels from extension %s into chip %d.\n", hdu->extname, i);
            if (! readExtension(hdu, fits)) {
                psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname);
                return false;
            }
        }

        pmConceptsReadChip(chip, db);

        // Iterate over the chip
        psArray *cells = chip->cells;   // Array of cells
        for (int j = 0; j < cells->n; j++) {
            pmCell *cell = cells->data[j]; // The cell

            // Only read cells marked to "read"
            if (!cell || ! cell->process || cell->exists) {
                psTrace(__func__, 3, "Ignoring chip %d...\n", i);
                continue;
            }
            psTrace(__func__, 3, "Reading in cell %d...\n", j);

            if (cell->hdu) {
                hdu = cell->hdu;
                psTrace(__func__, 5, "Reading pixels from extension %s into cell %d.\n", hdu->extname,
                        j);
                if (! readExtension(hdu, fits)) {
                    psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n",
                            hdu->extname);
                    return false;
                }
            }
            psTrace(__func__, 5, "Reading concepts for cell %d...\n", j);
            pmConceptsReadCell(cell, db);

            psTrace(__func__, 5, "Allocating readouts for chip %d cell %d...\n", i, j);
            generateReadouts(cell, hdu);

            cell->exists = true;
        }
        chip->exists = true;
    }

    psTrace(__func__, 1, "Done reading FPA...\n");

    return true;
}

// Translate a name from the configuration file, containing something like "%a_%d" to a real value
psString p_pmFPATranslateName(const psString name, // The name to translate
                              const pmCell *cell // The cell for which to translate
                             )
{
    // %a is the FPA.NAME
    // %b is the CHIP.NAME
    // %c is the CELL.NAME
    // %d is the chip number
    // %e if the cell number
    // %f is the extension name

    psString translation = psMemIncrRefCounter(name); // The translated string
    char *temp = NULL;                  // Temporary string

    pmChip *chip = cell->parent;        // Chip of interest
    pmFPA *fpa = chip->parent;          // FPA of interest

    // FPA.NAME
    if ((temp = strstr(translation, "%a"))) {
        // This is not particularly friendly to the memory allocation, but the cache should make it OK,
        // and there's not much of it anyway...
        psFree(translation);
        translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string
        psString fpaName = psMetadataLookupStr(NULL, fpa->concepts, "FPA.NAME"); // The value of FPA.NAME
        psStringAppend(&translation, "%s%s", fpaName, temp + 2);
        // So "translation" now contains the first part, the replaced string, and the last part.
    }

    // CHIP.NAME
    if ((temp = strstr(translation, "%b"))) {
        // This is not particularly friendly to the memory allocation, but the cache should make it OK,
        // and there's not much of it anyway...
        psFree(translation);
        translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string
        psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // CHIP.NAME's value
        psStringAppend(&translation, "%s%s", chipName, temp + 2);
        // So "translation" now contains the first part, the replaced string, and the last part.
    }

    // CELL.NAME
    if ((temp = strstr(translation, "%c"))) {
        // This is not particularly friendly to the memory allocation, but the cache should make it OK,
        // and there's not much of it anyway...
        psFree(translation);
        translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string
        psString cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // CELL.NAME's value
        psStringAppend(&translation, "%s%s", cellName, temp + 2);
        // So "translation" now contains the first part, the replaced string, and the last part.
    }

    // Chip number
    if ((temp = strstr(translation, "%d"))) {
        // This is not particularly friendly to the memory allocation, but the cache should make it OK,
        // and there's not much of it anyway...
        psFree(translation);
        translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string
        // Search for the pointer to get the chip number
        int chipNum = -1;
        psArray *chips = fpa->chips;    // The array of chips
        for (int i = 0; i < chips->n && chipNum < 0; i++) {
            if (chips->data[i] == chip) {
                chipNum = i;
            }
        }
        if (chipNum < 0) {
            psError(PS_ERR_IO, true, "Unable to find chip to get name: %s\n", name);
            // Try to muddle on by leaving the number out
            psStringAppend(&translation, "%s", temp + 2);
        } else {
            psStringAppend(&translation, "%d%s", chipNum, temp + 2);
        }
        // So "translation" now contains the first part, the replaced string, and the last part.
    }

    // Cell number
    if ((temp = strstr(translation, "%e"))) {
        // This is not particularly friendly to the memory allocation, but the cache should make it OK,
        // and there's not much of it anyway...
        psFree(translation);
        translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string
        // Search for the pointer to get the cell number
        int cellNum = -1;
        psArray *cells = chip->cells;   // The array of cells
        for (int i = 0; i < cells->n && cellNum < 0; i++) {
            if (cells->data[i] == cell) {
                cellNum = i;
            }
        }
        if (cellNum < 0) {
            psError(PS_ERR_IO, true, "Unable to find cell to get name: %s\n", name);
            // Try to muddle on by leaving the number out
            psStringAppend(&translation, "%s", temp + 2);
        } else {
            psStringAppend(&translation, "%d%s", cellNum, temp + 2);
        }
        // So "translation" now contains the first part, the replaced string, and the last part.
    }

    // Extension name
    if ((temp = strstr(translation, "%f"))) {
        // This is not particularly friendly to the memory allocation, but the cache should make it OK,
        // and there's not much of it anyway...
        psFree(translation);
        translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string
        const char *extname = NULL;
        if (cell->hdu) {
            extname = cell->hdu->extname;
        } else if (chip->hdu) {
            extname = chip->hdu->extname;
        } else if (fpa->hdu) {
            extname = fpa->hdu->extname;
        }
        psStringAppend(&translation, "%s%s", extname, temp + 2);
        // So "translation" now contains the first part, the replaced string, and the last part.
    }

    return translation;
}

// Get filename and extension out of "somefile.fits:ext"
psString p_pmFPATranslateFileExt(psString *extName, // Extension name, to be returned
                                 const psString name, // The string to be translated and then parsed
                                 const pmCell *cell // The cell
                                )
{
    psString filenameExt = p_pmFPATranslateName(name, cell);
    const char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn
    psString filename = NULL;           // The filename
    if (extName && *extName) {
        psFree(*extName);
    }
    if (extName) {
        *extName = NULL;
    }

    if (colon) {
        filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon));
        if (strlen(colon) > 1) {
            if (extName) {
                *extName = psStringCopy(colon + 1);
            }
        }
    } else {
        filename = psMemIncrRefCounter(filenameExt);
    }

    return filename;
}

// Read a mask into the FPA
bool pmFPAReadMask(pmFPA *fpa,          // FPA to read into
                   psFits *source       // Source FITS file (for the original data)
                  )
{
    const psMetadata *camera = fpa->camera; // Camera configuration for FPA
    bool mdok = false;                  // Status of MD lookup

    // Get the required information from the camera configuration
    psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data
    if (! mdok || ! supps) {
        psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n");
        return false;
    }
    psString sourceType = psMetadataLookupStr(&mdok, supps, "MASK.SOURCE"); // Type of source: EXT | FILE
    if (! mdok || strlen(sourceType) <= 0) {
        psError(PS_ERR_IO, false, "Unable to find MASK.SOURCE in SUPPLEMENTARY section of camera "
                "configuration!\n");
        return false;
    }
    psString name = psMetadataLookupStr(&mdok, supps, "MASK.NAME"); // Name of mask
    if (! mdok || strlen(sourceType) <= 0) {
        psError(PS_ERR_IO, false, "Unable to find MASK.NAME in SUPPLEMENTARY section of camera "
                "configuration!\n");
        return false;
    }

    // Go through the FPA to each cell/readout to get the mask
    p_pmHDU *hdu = fpa->hdu;            // The HDU into which we will read the mask
    psArray *chips = fpa->chips;        // Array of chips
    for (int chipNum = 0; chipNum < chips->n; chipNum++) {
        pmChip *chip = chips->data[chipNum]; // The current chip of interest
        if (chip->process && !chip->exists) {
            if (chip->hdu) {
                hdu = chip->hdu;
            }
            psArray *cells = chip->cells;       // Array of cells
            for (int cellNum = 0; cellNum < cells->n; cellNum++) {
                pmCell *cell = cells->data[cellNum]; // The current cell of interest
                if (cell->process && !cell->exists) {
                    if (cell->hdu) {
                        hdu = cell->hdu;
                    }

                    // Now, need to find out where to get the pixels
                    psFits *maskSource = psMemIncrRefCounter(source); // Source of mask image
                    if (strcasecmp(sourceType, "FILE") == 0) {
                        // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt"
                        psString extname = NULL; // Extension name
                        psString filename = p_pmFPATranslateFileExt(&extname, name, cell); // The filename

                        psFree(maskSource);
                        maskSource = psFitsOpen(filename, "r");
                        if (extname) {
                            if (! psFitsMoveExtName(maskSource, extname)) {
                                psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read mask.\n",
                                         extname);
                                return false;
                            }
                        }
                        psFree(filename);
                        psFree(extname);
                    } else if (strncasecmp(sourceType, "EXT", 3) == 0) {
                        // Source is an extension in the original file
                        psString extname = p_pmFPATranslateName(name, cell);
                        psFitsMoveExtName(maskSource, extname);
                    }

                    // We've arrived where the pixels are.  Now we need to read them in.
                    psMetadata *header = psFitsReadHeader(NULL, maskSource); // The header
                    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 for "
                                 "mask (%s)!\n", name);
                    }
                    if (nAxis != 2 && nAxis != 3) {
                        psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into "
                                 "a single image anyway.\n");
                    }

                    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");
                            // Try to proceed by taking only the first plane
                            numPlanes = 1;
                        }
                        if (numPlanes != 1 && numPlanes != cell->readouts->n) {
                            psError(PS_ERR_IO, false, "Number of masks (%d) does not match number of "
                                    "readouts (%d)\n", numPlanes, cell->readouts->n);
                            // Try to proceed by taking only the first plane
                            numPlanes = 1;
                        }
                    }

                    hdu->masks = psArrayAlloc(hdu->images->n);

                    // Read each plane into the array
                    psArray *readouts = cell->readouts; // The array of readouts
                    for (int i = 0; i < hdu->masks->n; i++) {
                        psImage *mask = NULL; // The mask to be added
                        if (i < numPlanes) {
                            // Read the mask from the file
                            psTrace(__func__, 9, "Reading plane %d\n", i);
                            psRegion region = {0, 0, 0, 0};
                            mask = psFitsReadImage(NULL, maskSource, region, i);
                        } else {
                            // One mask in the file is provided for all planes in the original image
                            psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i);
                            psImage *original = hdu->masks->data[0];
                            mask = psImageCopy(NULL, original, original->type.type);
                        }
                        hdu->masks->data[0] = mask;
                        pmReadout *readout = readouts->data[i];
                        readout->mask = mask;
                        // Check the dimensions
                        // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size?
                        if (mask->numCols < readout->image->numCols ||
                                mask->numRows < readout->image->numRows) {
                            psError(PS_ERR_IO, false, "Mask size (%dx%d) not compatible with image (%dx%d)\n",
                                    mask->numCols, mask->numRows, readout->image->numCols,
                                    readout->image->numRows);
                            return false;
                        }
                    } // Iterating over readouts
                    psFree(maskSource);
                } // Valid cells
            } // Iterating over cells
        } // Valid chips
    } // Iterating over chips

    return true;
}


// Read a mask into the FPA
// This is just a copy of the above pmFPAReadMask, replacing "mask" with "weight" throughout.
bool pmFPAReadWeight(pmFPA *fpa,        // FPA to read into
                     psFits *source     // Source FITS file (for the original data)
                    )
{
    const psMetadata *camera = fpa->camera; // Camera configuration for FPA
    bool mdok = false;                  // Status of MD lookup

    // Get the required information from the camera configuration
    psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data
    if (! mdok || ! supps) {
        psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n");
        return false;
    }
    psString sourceType = psMetadataLookupStr(&mdok, supps, "WEIGHT.SOURCE"); // Type of source: EXT | FILE
    if (! mdok || strlen(sourceType) <= 0) {
        psError(PS_ERR_IO, false, "Unable to find WEIGHT.SOURCE in SUPPLEMENTARY section of camera "
                "configuration!\n");
        return false;
    }
    psString name = psMetadataLookupStr(&mdok, supps, "WEIGHT.NAME"); // Name of weight
    if (! mdok || strlen(sourceType) <= 0) {
        psError(PS_ERR_IO, false, "Unable to find WEIGHT.NAME in SUPPLEMENTARY section of camera "
                "configuration!\n");
        return false;
    }

    // Go through the FPA to each cell/readout to get the weight
    p_pmHDU *hdu = fpa->hdu;            // The HDU into which we will read the weight
    psArray *chips = fpa->chips;        // Array of chips
    for (int chipNum = 0; chipNum < chips->n; chipNum++) {
        pmChip *chip = chips->data[chipNum]; // The current chip of interest
        if (chip->process && !chip->exists) {
            if (chip->hdu) {
                hdu = chip->hdu;
            }
            psArray *cells = chip->cells;       // Array of cells
            for (int cellNum = 0; cellNum < cells->n; cellNum++) {
                pmCell *cell = cells->data[cellNum]; // The current cell of interest
                if (cell->process && !cell->exists) {
                    if (cell->hdu) {
                        hdu = cell->hdu;
                    }

                    // Now, need to find out where to get the pixels
                    psFits *weightSource = psMemIncrRefCounter(source); // Source of weight image
                    if (strcasecmp(sourceType, "FILE") == 0) {
                        // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt"
                        psString extname = NULL; // Extension name
                        psString filename = p_pmFPATranslateFileExt(&extname, name, cell); // The filename

                        psFree(weightSource);
                        weightSource = psFitsOpen(filename, "r");
                        if (extname) {
                            if (! psFitsMoveExtName(weightSource, extname)) {
                                psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read "
                                         "weight.\n", extname);
                                return false;
                            }
                        }
                        psFree(filename);
                        psFree(extname);
                    } else if (strncasecmp(sourceType, "EXT", 3) == 0) {
                        // Source is an extension in the original file
                        psString extname = p_pmFPATranslateName(name, cell);
                        psFitsMoveExtName(weightSource, extname);
                    }

                    // We've arrived where the pixels are.  Now we need to read them in.
                    psMetadata *header = psFitsReadHeader(NULL, weightSource); // The header
                    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 for "
                                 "weight (%s)!\n", name);
                    }
                    if (nAxis != 2 && nAxis != 3) {
                        psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into "
                                 "a single image anyway.\n");
                    }

                    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");
                            // Try to proceed by taking only the first plane
                            numPlanes = 1;
                        }
                        if (numPlanes != 1 && numPlanes != cell->readouts->n) {
                            psError(PS_ERR_IO, false, "Number of weights (%d) does not match number of "
                                    "readouts (%d)\n", numPlanes, cell->readouts->n);
                            // Try to proceed by taking only the first plane
                            numPlanes = 1;
                        }
                    }

                    hdu->weights = psArrayAlloc(hdu->images->n);

                    // Read each plane into the array
                    psArray *readouts = cell->readouts; // The array of readouts
                    for (int i = 0; i < hdu->weights->n; i++) {
                        psImage *weight = NULL; // The weight to be added
                        if (i < numPlanes) {
                            // Read the weight from the file
                            psTrace(__func__, 9, "Reading plane %d\n", i);
                            psRegion region = {0, 0, 0, 0};
                            weight = psFitsReadImage(NULL, weightSource, region, i);
                        } else {
                            // One weight in the file is provided for all planes in the original image
                            psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i);
                            psImage *original = hdu->weights->data[0];
                            weight = psImageCopy(NULL, original, original->type.type);
                        }
                        hdu->weights->data[0] = weight;
                        pmReadout *readout = readouts->data[i];
                        readout->weight = weight;
                        // Check the dimensions
                        // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size?
                        if (weight->numCols < readout->image->numCols ||
                                weight->numRows < readout->image->numRows) {
                            psError(PS_ERR_IO, false, "Weight size (%dx%d) not compatible with image "
                                    "(%dx%d)\n", weight->numCols, weight->numRows, readout->image->numCols,
                                    readout->image->numRows);
                            return false;
                        }
                    } // Iterating over readouts
                    psFree(weightSource);
                } // Valid cells
            } // Iterating over cells
        } // Valid chips
    } // Iterating over chips

    return true;
}
