Index: trunk/archive/scripts/src/phase2/pmFPARead.c =================================================================== --- trunk/archive/scripts/src/phase2/pmFPARead.c (revision 5371) +++ trunk/archive/scripts/src/phase2/pmFPARead.c (revision 5564) @@ -17,9 +17,9 @@ // Read a FITS extension into a chip -static bool readExtension(pmPixelData *pixelData, // Pixel data into which to read +static bool readExtension(p_pmHDU *hdu, // Pixel data into which to read psFits *fits // The FITS file from which to read ) { - const char *extName = pixelData->extname; // Extension name + const char *extName = hdu->extname; // Extension name psTrace(__func__, 7, "Moving to extension %s...\n", extName); @@ -95,6 +95,6 @@ // Update the HDU with the new data - pixelData->images = pixels; - pixelData->header = header; + hdu->images = pixels; + hdu->header = header; // XXX: Insert mask and weight reading here @@ -105,9 +105,9 @@ // Portion out an image into the cell static bool generateReadouts(pmCell *cell, // The cell that gets its bits - pmPixelData *data // Pixel data, containing image, mask, weights + p_pmHDU *hdu // Pixel data, containing image, mask, weights ) { - psArray *images = data->images; // Array of images (each of which is a readout) - psArray *masks = data->masks; // Array of masks (one for each readout) + psArray *images = hdu->images; // Array of images (each of which is a readout) + psArray *masks = hdu->masks; // Array of masks (one for each readout) // Iterate over each of the image planes for (int i = 0; i < images->n; i++) { @@ -118,5 +118,5 @@ } - pmReadout *readout = pmReadoutAlloc(cell, image, mask, 0,0,0,0,1,1); + pmReadout *readout = pmReadoutAlloc(cell, image, mask, 0,0,1,1); psFree(readout); // This seems silly, but the alloc registers it with the cell, so we // don't have to do anything with it. Perhaps we should change @@ -137,5 +137,5 @@ ) { - pmPixelData *pixelData = NULL; // Pixel data from FITS file + p_pmHDU *hdu = NULL; // Pixel data from FITS file // Read the PHU, if required @@ -152,9 +152,9 @@ // Read in.... psTrace(__func__, 1, "Working on FPA...\n"); - if (fpa->data) { - pixelData= fpa->data; - psTrace(__func__, 3, "Reading pixels from extension %s into FPA.\n", pixelData->extname); - if (! readExtension(pixelData, fits)) { - psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", pixelData->extname); + 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; } @@ -174,9 +174,9 @@ psTrace(__func__, 2, "Reading in chip %d...\n", i); - if (chip->data) { - pixelData = chip->data; - psTrace(__func__, 3, "Reading pixels from extension %s into chip %d.\n", pixelData->extname, i); - if (! readExtension(pixelData, fits)) { - psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", pixelData->extname); + 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; } @@ -196,11 +196,11 @@ psTrace(__func__, 3, "Reading in cell %d...\n", j); - if (cell->data) { - pixelData = cell->data; - psTrace(__func__, 5, "Reading pixels from extension %s into cell %d.\n", pixelData->extname, + if (cell->hdu) { + hdu = cell->hdu; + psTrace(__func__, 5, "Reading pixels from extension %s into cell %d.\n", hdu->extname, j); - if (! readExtension(pixelData, fits)) { + if (! readExtension(hdu, fits)) { psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", - pixelData->extname); + hdu->extname); return false; } @@ -210,5 +210,5 @@ psTrace(__func__, 5, "Allocating readouts for chip %d cell %d...\n", i, j); - generateReadouts(cell, pixelData); + generateReadouts(cell, hdu); } } @@ -216,2 +216,416 @@ return true; } + +// Translate a name from the configuration file, containing something like "%a_%d" to a real value +psString p_pmFPATranslateName(psString name, // The name to translate + 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 = psMetadataLookupString(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 = psMetadataLookupString(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 = psMetadataLookupString(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; +} + +// 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 = psMetadataLookupString(&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 = psMetadataLookupString(&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->valid) { + 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->valid) { + if (cell->hdu) { + hdu = cell->hdu; + } + + // Now, need to find out where to get the pixels + psFits *maskSource = source; // Source of mask image + if (strcasecmp(sourceType, "FILE") == 0) { + // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt" + psString filenameExt = p_pmFPATranslateName(name, cell); + char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn + psString filename = NULL; // The filename + psString extname = NULL;// The extenstion name + if (colon) { + filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); + if (strlen(colon) > 1) { + extname = psStringCopy(colon + 1); + } + } else { + filename = psMemIncrRefCounter(filenameExt); + } + + maskSource = psFitsAlloc(filename); + 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); + psFree(filenameExt); + } 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 + } // 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 = psMetadataLookupString(&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 = psMetadataLookupString(&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->valid) { + 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->valid) { + if (cell->hdu) { + hdu = cell->hdu; + } + + // Now, need to find out where to get the pixels + psFits *weightSource = source; // Source of weight image + if (strcasecmp(sourceType, "FILE") == 0) { + // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt" + psString filenameExt = p_pmFPATranslateName(name, cell); + char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn + psString filename = NULL; // The filename + psString extname = NULL;// The extenstion name + if (colon) { + filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); + if (strlen(colon) > 1) { + extname = psStringCopy(colon + 1); + } + } else { + filename = psMemIncrRefCounter(filenameExt); + } + + weightSource = psFitsAlloc(filename); + 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); + psFree(filenameExt); + } 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 + } // Valid cells + } // Iterating over cells + } // Valid chips + } // Iterating over chips + + return true; +}