Index: /trunk/archive/scripts/src/phase2/Makefile.am =================================================================== --- /trunk/archive/scripts/src/phase2/Makefile.am (revision 5785) +++ /trunk/archive/scripts/src/phase2/Makefile.am (revision 5786) @@ -12,5 +12,4 @@ pmFPAConceptsSet.c \ pmFPAConstruct.c \ - pmFPAMorph.c \ pmFPARead.c \ pmFPAWrite.c \ @@ -29,5 +28,4 @@ pmFPAConceptsSet.h \ pmFPAConstruct.h \ - pmFPAMorph.h \ pmFPARead.h \ pmFPAWrite.h \ Index: /trunk/archive/scripts/src/phase2/ipprc.config =================================================================== --- /trunk/archive/scripts/src/phase2/ipprc.config (revision 5785) +++ /trunk/archive/scripts/src/phase2/ipprc.config (revision 5786) @@ -23,7 +23,13 @@ pap S32 10 pmFPAPrint S32 10 + pmFPAWrite S32 10 # pmConfigRead S32 10 # pmFPAWriteMask S32 10 # pmFPAWriteWeight S32 10 pmChipMosaic S32 10 + pmFPAMorph S32 10 + spliceCells S32 10 + spliceImage S32 10 + setConceptItem S32 10 +# pap_psMetadataCopy S32 9 END Index: /trunk/archive/scripts/src/phase2/megacam_raw.config =================================================================== --- /trunk/archive/scripts/src/phase2/megacam_raw.config (revision 5785) +++ /trunk/archive/scripts/src/phase2/megacam_raw.config (revision 5786) @@ -133,5 +133,5 @@ # Default PS concepts that may be specified by value DEFAULTS METADATA - CELL.READDIR S32 2 # Cell is read in x direction + CELL.READDIR S32 1 # Cell is read in x direction CELL.BAD S32 0 CELL.TIMESYS STR UTC Index: /trunk/archive/scripts/src/phase2/nonlin.dat =================================================================== --- /trunk/archive/scripts/src/phase2/nonlin.dat (revision 5785) +++ /trunk/archive/scripts/src/phase2/nonlin.dat (revision 5786) @@ -103,4 +103,4 @@ 9900 99 10000 100 -10001 -1 -1e6 -1 +10001 100 +1e6 100 Index: /trunk/archive/scripts/src/phase2/papPhase2.c =================================================================== --- /trunk/archive/scripts/src/phase2/papPhase2.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/papPhase2.c (revision 5786) @@ -17,5 +17,5 @@ #include "pmSubtractBias.h" #include "pmChipMosaic.h" -#include "pmFPAMorph.h" +//#include "pmFPAMorph.h" // Phase 2 needs to: @@ -39,17 +39,17 @@ // phase2 INPUT.fits OUTPUT.fits -bias BIAS.fits -dark DARK.fits -flat FLAT.fits -chip CHIP -// +// // Most are self-explanatory. "-chip" says "only work on this particular chip". -#define RECIPE "PHASE2" // Name of the recipe to use +#define RECIPE "PHASE2" // Name of the recipe to use static psMemId memPrintAlloc(const psMemBlock *mb) { printf("Allocated memory block %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); return 0; } @@ -57,8 +57,8 @@ { printf("Freed memory block %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); return 0; } @@ -67,8 +67,8 @@ psMemBlock *mb = ((psMemBlock*)ptr) - 1; printf("Memory block %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); } @@ -102,12 +102,12 @@ psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", 0, "Chip number to process (if positive)", -1); if (! psArgumentParse(arguments, &argc, argv) || argc != 3) { - printf("\nPan-STARRS Phase 2 processing\n\n"); - printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]); - psArgumentHelp(arguments); - psFree(arguments); - exit(EXIT_FAILURE); - } - const char *inputName = argv[1]; // Name of input image - const char *outputName = argv[2]; // Name of output image + printf("\nPan-STARRS Phase 2 processing\n\n"); + printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]); + psArgumentHelp(arguments); + psFree(arguments); + exit(EXIT_FAILURE); + } + const char *inputName = argv[1]; // Name of input image + const char *outputName = argv[2]; // Name of output image const char *biasName = psMetadataLookupString(NULL, arguments, "-bias"); // Name of bias image const char *darkName = psMetadataLookupString(NULL, arguments, "-dark"); // Name of dark image @@ -123,42 +123,34 @@ // Open the input -#ifdef PRODUCTION psFits *inputFile = psFitsOpen(inputName, "r"); // File handle for FITS file -#else - psFits *inputFile = psFitsAlloc(inputName); // File handle for FITS file -#endif if (! inputFile) { - psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName); + psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName); exit(EXIT_FAILURE); } psMetadata *header = psFitsReadHeader(NULL, inputFile); // FITS header #if PRODUCTION - psDB *database = pmConfigDB(site); // Database handle + psDB *database = pmConfigDB(site); // Database handle #else - psDB *database = NULL; // Database handle + psDB *database = NULL; // Database handle #endif // Open the output and output mask // We do it here so that we don't process the whole lot and then find out we can't open the output file -#ifdef PRODUCTION psFits *outputFile = psFitsOpen(outputName, "w"); -#else - psFits *outputFile = psFitsAlloc(outputName); -#endif if (! outputFile) { - psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName); - exit(EXIT_FAILURE); + psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName); + exit(EXIT_FAILURE); } // Get camera configuration from header if not already defined if (! camera) { - camera = pmConfigCameraFromHeader(site, header); - if (! camera) { - psErrorStackPrint(stderr, "Can't find camera configuration!\n"); - exit(EXIT_FAILURE); - } + camera = pmConfigCameraFromHeader(site, header); + if (! camera) { + psErrorStackPrint(stderr, "Can't find camera configuration!\n"); + exit(EXIT_FAILURE); + } } else if (! pmConfigValidateCamera(camera, header)) { - psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName); - exit(EXIT_FAILURE); + psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName); + exit(EXIT_FAILURE); } if (! recipe && !(recipe = pmConfigRecipeFromCamera(camera, RECIPE))) { @@ -175,106 +167,106 @@ // Set various tasks - bool doMask = false; // Mask bad pixles - bool doNonLin = false; // Non-linearity correction - bool doBias = false; // Bias subtraction - bool doDark = false; // Dark subtraction - bool doOverscan = false; // Overscan subtraction - bool doFlat = false; // Flat-field normalisation - bool doFringe = false; // Fringe subtraction - bool doSource = false; // Source identification and photometry - bool doAstrom = false; // Astrometry + bool doMask = false; // Mask bad pixles + bool doNonLin = false; // Non-linearity correction + bool doBias = false; // Bias subtraction + bool doDark = false; // Dark subtraction + bool doOverscan = false; // Overscan subtraction + bool doFlat = false; // Flat-field normalisation + bool doFringe = false; // Fringe subtraction + bool doSource = false; // Source identification and photometry + bool doAstrom = false; // Astrometry pmOverscanAxis overscanMode = PM_OVERSCAN_NONE; // Axis for overscan pmFit overscanFitType = PM_FIT_NONE; // Fit type for overscan - int overscanBins = 1; // Number of pixels per bin for overscan - psStats *overscanStats = NULL; // Statistics for overscan - void *overscanFit = NULL; // Overscan fit (polynomial or spline) + int overscanBins = 1; // Number of pixels per bin for overscan + psStats *overscanStats = NULL; // Statistics for overscan + void *overscanFit = NULL; // Overscan fit (polynomial or spline) if (psMetadataLookupBool(NULL, recipe, "MASK")) { - if (strlen(maskName) > 0) { - doMask = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking " - "will be performed.\n"); - } + if (strlen(maskName) > 0) { + doMask = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking " + "will be performed.\n"); + } } if (psMetadataLookupBool(NULL, recipe, "NONLIN")) { - doNonLin = true; + doNonLin = true; } if (psMetadataLookupBool(NULL, recipe, "BIAS")) { - if (strlen(biasName) > 0) { - doBias = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- " - "no bias subtraction will be performed.\n"); - } + if (strlen(biasName) > 0) { + doBias = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- " + "no bias subtraction will be performed.\n"); + } } if (psMetadataLookupBool(NULL, recipe, "DARK")) { - if (strlen(darkName) > 0) { - doDark = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- " - "no dark subtraction will be performed.\n"); - } + if (strlen(darkName) > 0) { + doDark = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- " + "no dark subtraction will be performed.\n"); + } } if (psMetadataLookupBool(NULL, recipe, "OVERSCAN")) { - doOverscan = true; - psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE"); - if (strcasecmp(mode, "INDIVIDUAL") == 0) { - overscanMode = PM_OVERSCAN_ROWS; - // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip - } else if (strcasecmp(mode, "ALL") == 0) { - overscanMode = PM_OVERSCAN_ALL; - } else if (strcasecmp(mode, "NONE") != 0) { - psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:" - " assuming NONE.\n", mode); - } - psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT"); - if (strcasecmp(fit, "POLYNOMIAL") == 0) { - overscanFitType = PM_FIT_POLYNOMIAL; - int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit - overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); - } else if (strcasecmp(fit, "SPLINE") == 0) { - overscanFitType = PM_FIT_SPLINE; - int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit - overscanFit = NULL; -// overscanFit = psSpline1DAlloc(); - } else if (strcasecmp(fit, "NONE") != 0) { - psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:" - " assuming NONE.\n", fit); - } - overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN"); - if (overscanBins <= 0) { - psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins); - overscanBins = 1; - } - psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT"); - if (strcasecmp(stat, "MEAN")) { - overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - } else if (strcasecmp(stat, "MEDIAN")) { - overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); - } else { - psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat); - overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - } + doOverscan = true; + psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE"); + if (strcasecmp(mode, "INDIVIDUAL") == 0) { + overscanMode = PM_OVERSCAN_ROWS; + // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip + } else if (strcasecmp(mode, "ALL") == 0) { + overscanMode = PM_OVERSCAN_ALL; + } else if (strcasecmp(mode, "NONE") != 0) { + psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:" + " assuming NONE.\n", mode); + } + psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT"); + if (strcasecmp(fit, "POLYNOMIAL") == 0) { + overscanFitType = PM_FIT_POLYNOMIAL; + int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit + overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } else if (strcasecmp(fit, "SPLINE") == 0) { + overscanFitType = PM_FIT_SPLINE; + int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit + overscanFit = NULL; +// overscanFit = psSpline1DAlloc(); + } else if (strcasecmp(fit, "NONE") != 0) { + psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:" + " assuming NONE.\n", fit); + } + overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN"); + if (overscanBins <= 0) { + psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins); + overscanBins = 1; + } + psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT"); + if (strcasecmp(stat, "MEAN")) { + overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + } else if (strcasecmp(stat, "MEDIAN")) { + overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + } else { + psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat); + overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + } } if (psMetadataLookupBool(NULL, recipe, "FLAT")) { - if (strlen(flatName) > 0) { - doFlat = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- " - "no flat-fielding will be performed.\n"); - } + if (strlen(flatName) > 0) { + doFlat = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- " + "no flat-fielding will be performed.\n"); + } } // Chip selection if (chipNum >= 0) { - if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) || - ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) || - ! pmFPASelectChip(mask, chipNum)) { - psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum); + if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) || + ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) || + ! pmFPASelectChip(mask, chipNum)) { + psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum); } @@ -284,16 +276,8 @@ psLogMsg("phase2", PS_LOG_INFO, "Opening input image: %s\n", inputName); if (! pmFPARead(input, inputFile, header, database)) { - psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName); - exit(EXIT_FAILURE); - } - - pmFPAPrint(input); - exit(1); - -#ifdef PRODUCTION + psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName); + exit(EXIT_FAILURE); + } psFitsClose(inputFile); -#else - psFree(inputFile); -#endif #if 0 @@ -302,47 +286,47 @@ //#else { - // Generate mask and weight frame - p_pmHDU *hdu = NULL; - if (input->hdu) { - hdu = input->hdu; - } - psArray *chips = input->chips; - for (int chipNum = 0; chipNum < chips->n; chipNum++) { - pmChip *chip = chips->data[chipNum]; - if (chip->valid) { - if (chip->hdu) { - hdu = chip->hdu; - } - psArray *cells = chip->cells; - for (int cellNum = 0; cellNum < cells->n; cellNum++) { - pmCell *cell = cells->data[cellNum]; - if (cell->valid) { - if (cell->hdu) { - hdu = cell->hdu; - } - - hdu->masks = psArrayAlloc(hdu->images->n); - hdu->weights = psArrayAlloc(hdu->images->n); - psArray *readouts = cell->readouts; - for (int readNum = 0; readNum < hdu->images->n; readNum++) { - psImage *image = hdu->images->data[readNum]; - psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); - psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); - pmReadout *readout = readouts->data[readNum]; - for (int j = 0; j < image->numRows; j++) { - for (int i = 0; i < image->numCols; i++) { - mask->data.U8[j][i] = j + i; - weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]); - } - } - hdu->masks->data[readNum] = mask; - hdu->weights->data[readNum] = weight; - readout->mask = psMemIncrRefCounter(mask); - readout->weight = psMemIncrRefCounter(weight); - } - } - } - } - } + // Generate mask and weight frame + p_pmHDU *hdu = NULL; + if (input->hdu) { + hdu = input->hdu; + } + psArray *chips = input->chips; + for (int chipNum = 0; chipNum < chips->n; chipNum++) { + pmChip *chip = chips->data[chipNum]; + if (chip->valid) { + if (chip->hdu) { + hdu = chip->hdu; + } + psArray *cells = chip->cells; + for (int cellNum = 0; cellNum < cells->n; cellNum++) { + pmCell *cell = cells->data[cellNum]; + if (cell->valid) { + if (cell->hdu) { + hdu = cell->hdu; + } + + hdu->masks = psArrayAlloc(hdu->images->n); + hdu->weights = psArrayAlloc(hdu->images->n); + psArray *readouts = cell->readouts; + for (int readNum = 0; readNum < hdu->images->n; readNum++) { + psImage *image = hdu->images->data[readNum]; + psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); + psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); + pmReadout *readout = readouts->data[readNum]; + for (int j = 0; j < image->numRows; j++) { + for (int i = 0; i < image->numCols; i++) { + mask->data.U8[j][i] = j + i; + weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]); + } + } + hdu->masks->data[readNum] = mask; + hdu->weights->data[readNum] = weight; + readout->mask = psMemIncrRefCounter(mask); + readout->weight = psMemIncrRefCounter(weight); + } + } + } + } + } } #endif @@ -354,477 +338,455 @@ // Load the calibration frames, if required if (doBias) { + psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias + psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias + if (! pmConfigValidateCamera(camera, biasHeader)) { + psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName); + if (! pmFPARead(bias, biasFile, biasHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName); + exit(EXIT_FAILURE); + } + psFree(biasHeader); + psFitsClose(biasFile); + } + + if (doDark) { + psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark + psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark + if (! pmConfigValidateCamera(camera, darkHeader)) { + psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName); + if (! pmFPARead(dark, darkFile, darkHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName); + exit(EXIT_FAILURE); + } + psFree(darkHeader); + psFitsClose(darkFile); + } + + if (doMask) { + psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask + psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask + if (! pmConfigValidateCamera(camera, maskHeader)) { + psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName); + if (! pmFPARead(mask, maskFile, maskHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName); + exit(EXIT_FAILURE); + } + psFree(maskHeader); + psFitsClose(maskFile); + } + + if (doFlat) { + psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat + if (! flatFile) { + psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName); + exit(EXIT_FAILURE); + } + psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat + if (! pmConfigValidateCamera(camera, flatHeader)) { + psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName); + if (! pmFPARead(flat, flatFile, flatHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName); + exit(EXIT_FAILURE); + } + psFree(flatHeader); + psFitsClose(flatFile); + } + + psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2")); + + psArray *inputChips = input->chips; // Array of chips in input image + psArray *biasChips = bias->chips; // Array of chips in bias image + psArray *darkChips = dark->chips; // Array of chips in dark image + psArray *maskChips = mask->chips; // Array of chips in mask image + psArray *flatChips = flat->chips; // Array of chips in flat image + for (int i = 0; i < inputChips->n; i++) { + pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image + pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image + pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image + pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image + pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image + + if (! inputChip->valid) { + continue; + } + + psArray *inputCells = inputChip->cells; // Array of cells in input image + psArray *biasCells = biasChip->cells; // Array of cells in bias image + psArray *darkCells = darkChip->cells; // Array of cells in dark image + psArray *maskCells = maskChip->cells; // Array of cells in mask image + psArray *flatCells = flatChip->cells; // Array of cells in flat image + + for (int j = 0; j < inputCells->n; j++) { + pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image + pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image + pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag + pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image + pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image + + if (! inputCell->valid) { + continue; + } + + psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image + if (doBias && biasCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" + " used."); + } + if (doDark && darkCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" + " used."); + } + if (doMask && maskCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" + " used."); + } + if (doFlat && flatCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be" + " used."); + } + + pmReadout *biasReadout = NULL; // Bias readout + pmReadout *maskReadout = NULL; // Mask readout + pmReadout *flatReadout = NULL; // Flat readout + psImage *biasImage = NULL; // Bias pixels + psImage *darkImage = NULL; // Dark pixels + float darkTime = 1.0; // Dark time for dark image + psImage *maskImage = NULL; // Mask pixels + + if (doBias) { + biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image + biasImage = biasReadout->image; + } + if (doDark) { + pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image + darkImage = darkReadout->image; + darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME"); + if (darkTime <= 0.0) { + psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime); + exit(EXIT_FAILURE); + } + } + if (doMask) { + maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image + maskImage = maskReadout->image; + } + if (doFlat) { + flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image + } + + for (int k = 0; k < inputReadouts->n; k++) { + pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image + psImage *inputImage = inputReadout->image; // The actual image + + // Mask bad pixels + if (doMask) { + float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION"); + + // Need to change this later to grow the mask by the size of the convolution kernel + (void)pmMaskBadPixels(inputReadout, maskImage, + PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation, + PM_MASK_TRAP, 0); + } + + // Non-linearity correction + if (doNonLin) { + psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA"); + if (! dataItem) { + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to " + "find NONLIN.DATA in recipe --- ignored.\n"); + } else { + switch (dataItem->type) { + case PS_DATA_VECTOR: + { + // These are the polynomial coefficients + psVector *coeff = dataItem->data.V; // The coefficient vector + if (coeff->type.type != PS_TYPE_F64) { + psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version + coeff = temp; + } + psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, + PS_POLYNOMIAL_ORD); + psFree(correction->coeff); + correction->coeff = psMemIncrRefCounter(coeff->data.F64); + (void)pmNonLinearityPolynomial(inputReadout, correction); + psFree(coeff); + psFree(correction); + } + break; + case PS_DATA_STRING: + { + // This is a filename + psString name = dataItem->data.V; // Filename + psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0); + if (psLookupTableRead(table) <= 0) { + psErrorStackPrint(stderr, "Unable to read non-linearity correction file " + "%s --- ignored\n", name); + } else { #ifdef PRODUCTION - psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias + (void)pmNonLinearityLookup(inputReadout, table); #else - psFits *biasFile = psFitsAlloc(biasName); -#endif - psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias - if (! pmConfigValidateCamera(camera, biasHeader)) { - psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName); - if (! pmFPARead(bias, biasFile, biasHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName); - exit(EXIT_FAILURE); - } - psFree(biasHeader); + psVector *influx = table->values->data[0]; + psVector *outflux = table->values->data[1]; + (void)pmNonLinearityLookup(inputReadout, table->values->data[0], + table->values->data[1]); +#endif + } + psFree(table); + } + break; + case PS_DATA_METADATA: + { + // This is a menu + psMetadata *options = dataItem->data.V; // Options with concept values as keys + bool mdok = false; // Success of MD lookup + psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE"); + if (! mdok || ! concept) { + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " + "unable to find NONLIN.SOURCE in recipe --- ignored.\n"); + } else { + psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept); + if (! conceptValueItem) { + psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s " + "for non-linearity correction --- ignored.\n", concept); + } else if (conceptValueItem->type != PS_DATA_STRING) { + psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as" + " expected for non-linearity correction --- ignored.\n", + concept); + } else { + psString conceptValue = conceptValueItem->data.V; + // Get the value of the concept + psMetadataItem *optionItem = psMetadataLookup(options, conceptValue); + if (!optionItem) { + psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA" + " --- ignored.\n", conceptValue); + } else if (optionItem->type == PS_DATA_VECTOR) { + // These are the polynomial coefficients + psVector *coeff = optionItem->data.V; // The coefficient vector + if (coeff->type.type != PS_TYPE_F64) { + psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); + coeff = temp; + } + // Polynomial correction + psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, + PS_POLYNOMIAL_ORD); + psFree(correction->coeff); + correction->coeff = psMemIncrRefCounter(coeff->data.F64); + (void)pmNonLinearityPolynomial(inputReadout, correction); + psFree(coeff); + psFree(correction); + } else if (optionItem->type == PS_DATA_STRING) { + // This is a filename + psString tableName = optionItem->data.V; // The filename + psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0); + int numLines = 0; // Number of lines read from table + if ((numLines = psLookupTableRead(table)) <= 0) { + psErrorStackPrint(stderr, "Unable to read non-linearity " + "correction file %s --- ignored\n", + tableName); + } else { #ifdef PRODUCTION - psFitsClose(biasFile); + (void)pmNonLinearityLookup(inputReadout, table); #else - psFree(biasFile); -#endif - } - - if (doDark) { + printf("XXX: Non-linearity correction from lookup table not " + "yet implemented.\n"); +#endif + } + psFree(table); + } else { + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction " + "desired but unable to interpret NONLIN.DATA for %s" + " --- ignored\n", conceptValue); + } + } + } + } + break; + default: + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " + "NONLIN.DATA is of invalid type --- ignored.\n"); + } + } + } // Non-linearity correction + + // Bias, dark and overscan subtraction are all merged. + + // Changed the definition of pmSubtractBias to do the dark subtraction as well, but + // it's not in there yet. Once it is, we can get rid of "pedestal". + #ifdef PRODUCTION - psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark + psImage *pedestal = NULL; // Pedestal image (bias + scaled dark) + if (doDark) { + // Dark time for input image + float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME"); + if (inputTime <= 0) { + psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n", + inputTime); + exit(EXIT_FAILURE); + } + + pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*", + psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); + } + if (doBias) { + if (pedestal) { + if (biasImage->numRows != darkImage->numRows || + biasImage->numCols != darkImage->numCols) { + psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: " + "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows, + darkImage->numCols, darkImage->numRows); + exit(EXIT_FAILURE); + } + (void)psBinaryOp(pedestal, pedestal, "+", biasImage); + } else { + pedestal = psMemIncrRefCounter(biasImage); + } + } +#endif + if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { + // Need to get the read direction + int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR"); + if (readdir == 1) { +// psmodule-0.8.0 has confused PM_OVERSCAN_ROWS and PM_OVERSCAN_COLUMNS; bug 608. +#ifdef PRODUCTION + overscanMode = PM_OVERSCAN_ROWS; #else - psFits *darkFile = psFitsAlloc(darkName); -#endif - psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark - if (! pmConfigValidateCamera(camera, darkHeader)) { - psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName); - if (! pmFPARead(dark, darkFile, darkHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName); - exit(EXIT_FAILURE); - } - psFree(darkHeader); + overscanMode = PM_OVERSCAN_COLUMNS; +#endif + } else if (readdir == 2) { #ifdef PRODUCTION - psFitsClose(darkFile); + overscanMode = PM_OVERSCAN_COLUMNS; #else - psFree(darkFile); -#endif - } - - if (doMask) { + overscanMode = PM_OVERSCAN_ROWS; +#endif + } else { + psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n", + readdir); + overscanMode = PM_OVERSCAN_ROWS; + } + } + + if (doBias || doOverscan || doDark) { + psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions #ifdef PRODUCTION - psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask + (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, + overscanStats, overscanBins, overscanFitType, pedestal); #else - psFits *maskFile = psFitsAlloc(maskName); -#endif - psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask - if (! pmConfigValidateCamera(camera, maskHeader)) { - psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName); - if (! pmFPARead(mask, maskFile, maskHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName); - exit(EXIT_FAILURE); - } - psFree(maskHeader); - -#ifdef PRODUCTION - psFitsClose(maskFile); -#else - psFree(maskFile); -#endif - } - - if (doFlat) { -#ifdef PRODUCTION - psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat -#else - psFits *flatFile = psFitsAlloc(flatName); -#endif - if (! flatFile) { - psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName); - exit(EXIT_FAILURE); - } - psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat - if (! pmConfigValidateCamera(camera, flatHeader)) { - psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName); - if (! pmFPARead(flat, flatFile, flatHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName); - exit(EXIT_FAILURE); - } - psFree(flatHeader); -#ifdef PRODUCTION - psFitsClose(flatFile); -#else - psFree(flatFile); -#endif - } - - psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2")); - - psArray *inputChips = input->chips; // Array of chips in input image - psArray *biasChips = bias->chips; // Array of chips in bias image - psArray *darkChips = dark->chips; // Array of chips in dark image - psArray *maskChips = mask->chips; // Array of chips in mask image - psArray *flatChips = flat->chips; // Array of chips in flat image - for (int i = 0; i < inputChips->n; i++) { - pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image - pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image - pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image - pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image - pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image - - if (! inputChip->valid) { - continue; - } - - psArray *inputCells = inputChip->cells; // Array of cells in input image - psArray *biasCells = biasChip->cells; // Array of cells in bias image - psArray *darkCells = darkChip->cells; // Array of cells in dark image - psArray *maskCells = maskChip->cells; // Array of cells in mask image - psArray *flatCells = flatChip->cells; // Array of cells in flat image - - for (int j = 0; j < inputCells->n; j++) { - pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image - pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image - pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag - pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image - pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image - - if (! inputCell->valid) { - continue; - } - - psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image - if (doBias && biasCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" - " used."); - } - if (doDark && darkCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" - " used."); - } - if (doMask && maskCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" - " used."); - } - if (doFlat && flatCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be" - " used."); - } - - pmReadout *biasReadout = NULL; // Bias readout - pmReadout *maskReadout = NULL; // Mask readout - pmReadout *flatReadout = NULL; // Flat readout - psImage *biasImage = NULL; // Bias pixels - psImage *darkImage = NULL; // Dark pixels - float darkTime = 1.0; // Dark time for dark image - psImage *maskImage = NULL; // Mask pixels - - if (doBias) { - biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image - biasImage = biasReadout->image; - } - if (doDark) { - pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image - darkImage = darkReadout->image; - darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME"); - if (darkTime <= 0.0) { - psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime); - exit(EXIT_FAILURE); - } - } - if (doMask) { - maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image - maskImage = maskReadout->image; - } - if (doFlat) { - flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image - } - - for (int k = 0; k < inputReadouts->n; k++) { - pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image - psImage *inputImage = inputReadout->image; // The actual image - - // Mask bad pixels - if (doMask) { - float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION"); - - // Need to change this later to grow the mask by the size of the convolution kernel - (void)pmMaskBadPixels(inputReadout, maskImage, - PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation, - PM_MASK_TRAP, 0); - } - - // Non-linearity correction - if (doNonLin) { - psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA"); - if (! dataItem) { - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to " - "find NONLIN.DATA in recipe --- ignored.\n"); - } else { - switch (dataItem->type) { - case PS_DATA_VECTOR: - { - // These are the polynomial coefficients - psVector *coeff = dataItem->data.V; // The coefficient vector - if (coeff->type.type != PS_TYPE_F64) { - psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version - coeff = temp; - } - psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, - PS_POLYNOMIAL_ORD); - psFree(correction->coeff); - correction->coeff = psMemIncrRefCounter(coeff->data.F64); - (void)pmNonLinearityPolynomial(inputReadout, correction); - psFree(coeff); - psFree(correction); - } - break; - case PS_DATA_STRING: - { - // This is a filename - psString name = dataItem->data.V; // Filename - psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0); - if (psLookupTableRead(table) <= 0) { - psErrorStackPrint(stderr, "Unable to read non-linearity correction file " - "%s --- ignored\n", name); - } else { -#ifdef PRODUCTION - (void)pmNonLinearityLookup(inputReadout, table); -#else - psVector *influx = table->values->data[0]; - psVector *outflux = table->values->data[1]; - (void)pmNonLinearityLookup(inputReadout, table->values->data[0], - table->values->data[1]); -#endif - } - psFree(table); - } - break; - case PS_DATA_METADATA: - { - // This is a menu - psMetadata *options = dataItem->data.V; // Options with concept values as keys - bool mdok = false; // Success of MD lookup - psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE"); - if (! mdok || ! concept) { - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " - "unable to find NONLIN.SOURCE in recipe --- ignored.\n"); - } else { - psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept); - if (! conceptValueItem) { - psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s " - "for non-linearity correction --- ignored.\n", concept); - } else if (conceptValueItem->type != PS_DATA_STRING) { - psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as" - " expected for non-linearity correction --- ignored.\n", - concept); - } else { - psString conceptValue = conceptValueItem->data.V; - // Get the value of the concept - psMetadataItem *optionItem = psMetadataLookup(options, conceptValue); - if (!optionItem) { - psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA" - " --- ignored.\n", conceptValue); - } else if (optionItem->type == PS_DATA_VECTOR) { - // These are the polynomial coefficients - psVector *coeff = optionItem->data.V; // The coefficient vector - if (coeff->type.type != PS_TYPE_F64) { - psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); - coeff = temp; - } - // Polynomial correction - psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, - PS_POLYNOMIAL_ORD); - psFree(correction->coeff); - correction->coeff = psMemIncrRefCounter(coeff->data.F64); - (void)pmNonLinearityPolynomial(inputReadout, correction); - psFree(coeff); - psFree(correction); - } else if (optionItem->type == PS_DATA_STRING) { - // This is a filename - psString tableName = optionItem->data.V; // The filename - psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0); - int numLines = 0; // Number of lines read from table - if ((numLines = psLookupTableRead(table)) <= 0) { - psErrorStackPrint(stderr, "Unable to read non-linearity " - "correction file %s --- ignored\n", - tableName); - } else { -#ifdef PRODUCTION - (void)pmNonLinearityLookup(inputReadout, table); -#else - printf("XXX: Non-linearity correction from lookup table not " - "yet implemented.\n"); -#endif - } - psFree(table); - } else { - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction " - "desired but unable to interpret NONLIN.DATA for %s" - " --- ignored\n", conceptValue); - } - } - } - } - break; - default: - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " - "NONLIN.DATA is of invalid type --- ignored.\n"); - } - } - } // Non-linearity correction - - // Bias, dark and overscan subtraction are all merged. - - // Changed the definition of pmSubtractBias to do the dark subtraction as well, but - // it's not in there yet. Once it is, we can get rid of "pedestal". - -#ifdef PRODUCTION - psImage *pedestal = NULL; // Pedestal image (bias + scaled dark) - if (doDark) { - // Dark time for input image - float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME"); - if (inputTime <= 0) { - psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n", - inputTime); - exit(EXIT_FAILURE); - } - - pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*", - psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); - } - if (doBias) { - if (pedestal) { - if (biasImage->numRows != darkImage->numRows || - biasImage->numCols != darkImage->numCols) { - psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: " - "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows, - darkImage->numCols, darkImage->numRows); - exit(EXIT_FAILURE); - } - (void)psBinaryOp(pedestal, pedestal, "+", biasImage); - } else { - pedestal = psMemIncrRefCounter(biasImage); - } - } -#endif - if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { - // Need to get the read direction - int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR"); - if (readdir == 1) { - overscanMode = PM_OVERSCAN_ROWS; - } else if (readdir == 2) { - overscanMode = PM_OVERSCAN_COLUMNS; - } else { - psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n", - readdir); - overscanMode = PM_OVERSCAN_ROWS; - } - } - - if (doBias || doOverscan || doDark) { - psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions -#ifdef PRODUCTION - (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, - overscanStats, overscanBins, overscanFitType, pedestal); -#else - (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, - overscanStats, overscanBins, overscanFitType, biasReadout); -#endif - psFree(inputOverscans); - - // Output overscan fit results, if required - if (doOverscan) { - if (! overscanFit) { - psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n"); - } else { - switch (overscanFitType) { - case PM_FIT_POLYNOMIAL: - { - psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial - psString coeffs = NULL; // String containing the coefficients - for (int i = 0; i < poly->nX; i++) { - psStringAppend(&coeffs, "%e ", poly->coeff[i]); - } - psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n", - coeffs); - psFree(coeffs); - } - break; - case PM_FIT_SPLINE: - { - psSpline1D *spline = (psSpline1D*)overscanFit; // The spline - psString coeffs = NULL; // String containing the coefficients - for (int i = 0; i < spline->n; i++) { - psPolynomial1D *poly = spline->spline[i]; // i-th polynomial - psStringAppend(&coeffs, "%d: ", i); - for (int j = 0; j < poly->nX; j++) { - psStringAppend(&coeffs, "%e ", poly->coeff[i]); - } - psStringAppend(&coeffs, "\n"); - } - psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n", - coeffs); - psFree(coeffs); - } - break; - case PM_FIT_NONE: - break; - default: - psAbort(__func__, "Should never get here!!!\n"); - } - } - } - } - - // XXX: Temporary: until the other functions are altered to do this themselves. - // Trim, so that flat, fringe etc computations are faster. + (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, + overscanStats, overscanBins, overscanFitType, biasReadout); +#endif + psFree(inputOverscans); + + // Output overscan fit results, if required + if (doOverscan) { + if (! overscanFit) { + psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n"); + } else { + switch (overscanFitType) { + case PM_FIT_POLYNOMIAL: + { + psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial + psString coeffs = NULL; // String containing the coefficients + for (int i = 0; i < poly->nX; i++) { + psStringAppend(&coeffs, "%e ", poly->coeff[i]); + } + psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n", + coeffs); + psFree(coeffs); + } + break; + case PM_FIT_SPLINE: + { + psSpline1D *spline = (psSpline1D*)overscanFit; // The spline + psString coeffs = NULL; // String containing the coefficients + for (int i = 0; i < spline->n; i++) { + psPolynomial1D *poly = spline->spline[i]; // i-th polynomial + psStringAppend(&coeffs, "%d: ", i); + for (int j = 0; j < poly->nX; j++) { + psStringAppend(&coeffs, "%e ", poly->coeff[i]); + } + psStringAppend(&coeffs, "\n"); + } + psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n", + coeffs); + psFree(coeffs); + } + break; + case PM_FIT_NONE: + break; + default: + psAbort(__func__, "Should never get here!!!\n"); + } + } + } + } + + // XXX: Temporary: until the other functions are altered to do this themselves. + // Trim, so that flat, fringe etc computations are faster. #if 0 #ifndef PRODUCTION - psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC"); - psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion); - psFree(inputReadout->image); - inputReadout->image = trimmed; -#endif -#endif - // Flat-field correction - if (doFlat) { - psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n"); + psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC"); + psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion); + psFree(inputReadout->image); + inputReadout->image = trimmed; +#endif +#endif + // Flat-field correction + if (doFlat) { + psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n"); #ifndef PRODUCTION - psImage *dummyImage = psImageAlloc(inputReadout->image->numCols, - inputReadout->image->numRows, - PS_TYPE_U8); - pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1); - (void)pmFlatField(inputReadout, dummyMask, flatReadout); - psFree(dummyMask); - psFree(dummyImage); -#endif - } - - // Fringe subtraction - if (doFringe) { + psImage *dummyImage = psImageAlloc(inputReadout->image->numCols, + inputReadout->image->numRows, + PS_TYPE_U8); + pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1); + (void)pmFlatField(inputReadout, dummyMask, flatReadout); + psFree(dummyMask); + psFree(dummyImage); +#endif + } + + // Fringe subtraction + if (doFringe) { #if 0 - // Placeholder - pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, - psU8 maskVal, int binFactor, psStats *stats, float clipSD); -#endif - } - - // Source identification and photometry - if (doSource) { - } - - // Astrometry - if (doAstrom) { - } - - - } // Iterating over readouts - } // Iterating over cells - - // Testing pmChipMosaic - { - psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip - psFits *mosaicFile = psFitsAlloc("mosaic.fits"); // FITS file to which to write - psFitsWriteImage(mosaicFile, NULL, mosaic, 0); - psFree(mosaicFile); - psFree(mosaic); - } - + // Placeholder + pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, + psU8 maskVal, int binFactor, psStats *stats, float clipSD); +#endif + } + + // Source identification and photometry + if (doSource) { + } + + // Astrometry + if (doAstrom) { + } + + + } // Iterating over readouts + } // Iterating over cells + +#if 0 + { + // Testing pmChipMosaic + psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip + psFits *mosaicFile = psFitsOpen("mosaic.fits", "w"); // FITS file to which to write + psFitsWriteImage(mosaicFile, NULL, mosaic, 0); + psFitsClose(mosaicFile); + psFree(mosaic); + } +#endif + } // Iterating over chips @@ -832,16 +794,20 @@ +#if 0 + { + // Testing pmFPAMorph + int nBad = 0; + psMetadata *newCamera = psMetadataConfigParse(NULL, &nBad, "megacam_splice.config", true); + pmFPA *newFPA = pmFPAConstruct(newCamera); + pmFPAMorph(newFPA, input, true, 0, 0); + //pmFPAPrint(newFPA); + psFits *newFile = psFitsOpen("morph.fits", "w"); + pmFPAWrite(newFile, newFPA, database); + psFitsClose(newFile); + psFree(newFPA); + } +#endif + #if 1 - // Testing pmFPAMorph - psMetadata *newCamera = psMetadataConfigParse(NULL, NULL, "megacam_splice.config", true); - pmFPA *newFPA = pmFPAConstruct(newCamera); - pmFPAMorph(newFPA, input, true, 0, 0); - psFits *newFile = psFitsAlloc("morph.fits"); - pmFPAWrite(newFile, newFPA, database); - psFree(newFile); - psFree(newFPA); -#endif - -#if 0 // Write the output pmFPAWrite(outputFile, input, database); @@ -852,9 +818,5 @@ psLogMsg("phase2", PS_LOG_INFO, "Output completed after %f sec.\n", psTimerMark("phase2")); -#ifdef PRODUCTION psFitsClose(outputFile); -#else - psFree(outputFile); -#endif psFree(arguments); @@ -874,5 +836,5 @@ psTimerStop(); -#if 1 +#if 0 psMemCheckCorruption(true); psMemBlock **leaks = NULL; // List of leaks @@ -881,10 +843,10 @@ #if 1 for (int i = 0; i < nLeaks; i++) { - psMemBlock *mb = leaks[i]; - printf("Memory leak %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + psMemBlock *mb = leaks[i]; + printf("Memory leak %ld (%ld):\n" + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); } #endif Index: /trunk/archive/scripts/src/phase2/phase2.config =================================================================== --- /trunk/archive/scripts/src/phase2/phase2.config (revision 5785) +++ /trunk/archive/scripts/src/phase2/phase2.config (revision 5786) @@ -3,9 +3,9 @@ # List of tasks to perform MASK BOOL FALSE # Mask bad pixels -NONLIN BOOL FALSE # Non-linearity correction -BIAS BOOL FALSE # Bias subtraction +NONLIN BOOL TRUE # Non-linearity correction +BIAS BOOL TRUE # Bias subtraction DARK BOOL FALSE # Dark subtraction -OVERSCAN BOOL FALSE # Overscan subtraction -FLAT BOOL FALSE # Flat-field normalisation +OVERSCAN BOOL TRUE # Overscan subtraction +FLAT BOOL TRUE # Flat-field normalisation FRINGE BOOL FALSE # Fringe subtraction SOURCE BOOL FALSE # Source identification and photometry Index: /trunk/archive/scripts/src/phase2/pmChipMosaic.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmChipMosaic.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmChipMosaic.c (revision 5786) @@ -1,3 +1,4 @@ #include +#include #include "pslib.h" @@ -14,16 +15,20 @@ } -psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic - int xBinChip, int yBinChip // Binning of mosaic image in x and y +// Mosaic multiple images, with flips, binning and offsets +psImage *p_pmImageMosaic(const psArray *source, // Images to splice in + const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? + const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of + // source images + int xBinTarget, int yBinTarget, // Binning in x and y of target images + const psVector *x0, const psVector *y0 // Offsets for source images on target ) { - psArray *cells = chip->cells; // The array of cells - psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates - psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates - psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x - psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y - psVector *xParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in x - psVector *yParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in y - psArray *trimsec = psArrayAlloc(cells->n); // Trim section + assert(source); + assert(xFlip && xFlip->type.type == PS_TYPE_U8); + assert(yFlip && yFlip->type.type == PS_TYPE_U8); + assert(xBinSource && xBinSource->type.type == PS_TYPE_S32); + assert(yBinSource && yBinSource->type.type == PS_TYPE_S32); + assert(x0 && x0->type.type == PS_TYPE_S32); + assert(y0 && y0->type.type == PS_TYPE_S32); // Get the maximum extent of the mosaic image @@ -32,93 +37,138 @@ int yMin = INT_MAX; int yMax = 0; - for (int i = 0; i < cells->n; i++) { - pmCell *cell = cells->data[i]; // The cell of interest - x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); - yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); - xParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - yParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - trimsec->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC")); - psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, ((psRegion*)trimsec->data[i])->x0, - ((psRegion*)trimsec->data[i])->x1, ((psRegion*)trimsec->data[i])->y0, - ((psRegion*)trimsec->data[i])->y1); - // Size of cell in x and y - int nx = (int)(((psRegion*)trimsec->data[i])->x1 - (int)((psRegion*)trimsec->data[i])->x0); - int ny = (int)(((psRegion*)trimsec->data[i])->y1 - (int)((psRegion*)trimsec->data[i])->y0); - psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i], - x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * nx, y0->data.S32[i], - y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny); + for (int i = 0; i < source->n; i++) { + psImage *image = source->data[i]; // The image of interest - COMPARE(x0->data.S32[i], xMin, xMax); - COMPARE(y0->data.S32[i], yMin, yMax); - // Subtract the parity to get the inclusive limit (not exclusive) - COMPARE(x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * nx - xParity->data.S32[i], xMin, xMax); - COMPARE(y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny - yParity->data.S32[i], yMin, yMax); + assert(image->type.type == PS_TYPE_F32); // Only implemented for F32 images so far. + + // Size of cell in x and y + int xParity = xFlip->data.U8[i] ? -1 : 1; + int yParity = yFlip->data.U8[i] ? -1 : 1; + psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i], + x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols, y0->data.S32[i], + y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows); + + COMPARE(x0->data.S32[i], xMin, xMax); + COMPARE(y0->data.S32[i], yMin, yMax); + // Subtract the parity to get the inclusive limit (not exclusive) + COMPARE(x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols - xParity, xMin, xMax); + COMPARE(y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows - yParity, yMin, yMax); } // Set up the image // Since both upper and lower values are inclusive, we need to add one to the size - float xSize = (float)(xMax - xMin + 1) / (float)xBinChip; + float xSize = (float)(xMax - xMin + 1) / (float)xBinTarget; if (xSize - (int)xSize > 0) { - xSize += 1; + xSize += 1; } - float ySize = (float)(yMax - yMin + 1) / (float)yBinChip; + float ySize = (float)(yMax - yMin + 1) / (float)yBinTarget; if (ySize - (int)ySize > 0) { - ySize += 1; + ySize += 1; } - psTrace(__func__, 3, "Mosaicked chip will be %dx%d\n", (int)xSize, (int)ySize); + psTrace(__func__, 3, "Spliced image will be %dx%d\n", (int)xSize, (int)ySize); psImage *mosaic = psImageAlloc((int)xSize, (int)ySize, PS_TYPE_F32); // The mosaic image psImageInit(mosaic, 0.0); - // Set the image + // Next pass through the images to do the mosaicking + for (int i = 0; i < source->n; i++) { + psImage *image = source->data[i]; // The image of interest + if (xBinSource->data.S32[i] == xBinTarget && yBinSource->data.S32[i] == yBinTarget && + xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) { + // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason + psImageOverlaySection(mosaic, image, x0->data.S32[i], y0->data.S32[i], "+"); + } else { + // We have to do the hard work ourself + for (int y = 0; y < image->numRows; y++) { + int yParity = yFlip->data.U8[i] ? -1 : 1; + float yTargetBase = (y0->data.S32[i] + yParity * yBinSource->data.S32[i] * y) / yBinTarget; + for (int x = 0; x < image->numCols; x++) { + int xParity = xFlip->data.U8[i] ? -1 : 1; + float xTargetBase = (x0->data.S32[i] + xParity * xBinSource->data.S32[i] * x) / + xBinTarget; + + // In case the original image is binned but the mosaic is not, we need to fill in the + // values in the mosaic. + for (int j = 0; j < yBinSource->data.S32[i]; j++) { + int yTarget = (int)(yTargetBase + yParity * (float)j / (float)yBinTarget); + for (int i = 0; i < xBinSource->data.S32[i]; i++) { + int xTarget = (int)(xTargetBase + xParity * (float)i / (float)xBinTarget); + + mosaic->data.F32[yTarget][xTarget] += image->data.F32[y][x]; + } + } // Iterating over mosaic image for binned input image + } + } // Iterating over input image + } + } + + return mosaic; +} + +psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic + int xBinChip, int yBinChip // Binning of mosaic image in x and y + ) +{ + psArray *cells = chip->cells; // The array of cells + psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked + psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates + psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates + psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x + psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y + psVector *xFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in x? + psVector *yFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in y? + + // Set up the required inputs for (int i = 0; i < cells->n; i++) { - pmCell *cell = cells->data[i]; // The cell of interest - psArray *readouts = cell->readouts; // The array of readouts - if (readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " - "be mosaicked.\n", i); - } - psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic - psImage *trimmed = psImageSubset(image, *(psRegion*)(trimsec->data[i])); // Trimmed image (no overscan) + pmCell *cell = cells->data[i]; // The cell of interest + x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); + yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); + int xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + if (xParity == 1) { + xFlip->data.U8[i] = 0; + } else if (xParity == -1) { + xFlip->data.U8[i] = 1; + } else { + psLogMsg(__func__, PS_LOG_WARN, "The x parity of cell %d is not +/- 1 (it's %d) --- " + "assuming +1.\n", i, xParity); + xFlip->data.U8[i] = 0; + } + if (yParity == 1) { + yFlip->data.U8[i] = 0; + } else if (yParity == -1) { + yFlip->data.U8[i] = 1; + } else { + psLogMsg(__func__, PS_LOG_WARN, "The y parity of cell %d is not +/- 1 (it's %d) --- " + "assuming +1.\n", i, yParity); + yFlip->data.U8[i] = 0; + } - if (xBin->data.S32[i] == xBinChip && yBin->data.S32[i] == yBinChip && xParity->data.S32[i] == 1 && - yParity->data.S32[i] == 1) { - // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason - psImageOverlaySection(mosaic, trimmed, x0->data.S32[i], y0->data.S32[i], "+"); - } else { - // We have to do the hard work ourself - for (int y = 0; y < trimmed->numRows; y++) { - float yTargetBase = (y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * y) / - yBinChip; - for (int x = 0; x < trimmed->numCols; x++) { - float xTargetBase = (x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * x) / - xBinChip; + // Trim the image to get rid of the overscan + psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); + psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, trimsec->x0, trimsec->x1, + trimsec->y0, trimsec->y1); + psArray *readouts = cell->readouts; // The array of readouts + if (readouts->n > 1) { + psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " + "be mosaicked.\n", i); + } + psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic + images->data[i] = psImageSubset(image, *trimsec); // Trimmed image + } - // In case the original image is binned but the mosaic is not, we need to fill in the - // values in the mosaic. - for (int j = 0; j < yBin->data.S32[i]; j++) { - int yTarget = (int)(yTargetBase + yParity->data.S32[i] * (float)j / (float)yBinChip); - for (int i = 0; i < xBin->data.S32[i]; i++) { - int xTarget = (int)(xTargetBase + - xParity->data.S32[i] * (float)i / (float)xBinChip); + // Mosaic the images together and we're done + psImage *mosaic = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0); - mosaic->data.F32[yTarget][xTarget] += trimmed->data.F32[y][x]; - } - } // Iterating over mosaic image for binned input image - } - } // Iterating over input image - } - psFree(trimmed); - } // Iterating over cells - + // Clean up psFree(x0); psFree(y0); psFree(xBin); psFree(yBin); - psFree(xParity); - psFree(yParity); - psFree(trimsec); + psFree(xFlip); + psFree(yFlip); + psFree(images); return mosaic; Index: /trunk/archive/scripts/src/phase2/pmChipMosaic.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmChipMosaic.h (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmChipMosaic.h (revision 5786) @@ -5,6 +5,16 @@ #include "pmFPA.h" -psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic - int xBinChip, int yBinChip // Binning of mosaic image in x and y +// Mosaic multiple images, with flips, binning and offsets +psImage *p_pmImageMosaic(const psArray *source, // Images to splice in + const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? + const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of + // source images + int xBinTarget, int yBinTarget, // Binning in x and y of target images + const psVector *x0, const psVector *y0 // Offsets for source images on target + ); + +// Mosaic all the cells in a chip together (neglecting the overscans) +psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic + int xBinChip, int yBinChip // Binning of mosaic image in x and y ); Index: /trunk/archive/scripts/src/phase2/pmFPA.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPA.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPA.c (revision 5786) @@ -77,4 +77,5 @@ fpa->chips = psArrayAlloc(0); + fpa->phu = NULL; fpa->hdu = NULL; @@ -89,7 +90,8 @@ psFree(fpa->concepts); - psFree((psPtr)fpa->camera); + psFree(fpa->camera); psFree(fpa->chips); + psFree(fpa->phu); psFree(fpa->hdu); } Index: /trunk/archive/scripts/src/phase2/pmFPA.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPA.h (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPA.h (revision 5786) @@ -19,29 +19,29 @@ typedef struct { // Astrometric transformations - psPlaneDistort* fromTangentPlane; // Transformation from tangent plane to focal plane - psPlaneDistort* toTangentPlane; // Transformation from focal plane to tangent plane - psProjection *projection; // Projection from tangent plane to sky + psPlaneDistort* fromTangentPlane; // Transformation from tangent plane to focal plane + psPlaneDistort* toTangentPlane; // Transformation from focal plane to tangent plane + psProjection *projection; // Projection from tangent plane to sky // Information - psMetadata *concepts; // Cache for PS concepts - psMetadata *analysis; // FPA-level analysis metadata - const psMetadata *camera; // Camera configuration - psArray *chips; // The chips - p_pmHDU *hdu; // FITS data - psMetadata *phu; // Primary Header + psMetadata *concepts; // Cache for PS concepts + psMetadata *analysis; // FPA-level analysis metadata + const psMetadata *camera; // Camera configuration + psArray *chips; // The chips + p_pmHDU *hdu; // FITS data + psMetadata *phu; // Primary Header } pmFPA; typedef struct { // Offset specifying position on focal plane - int col0; // Offset from the left of FPA. - int row0; // Offset from the bottom of FPA. + int col0; // Offset from the left of FPA. + int row0; // Offset from the bottom of FPA. // Astrometric transformations - psPlaneTransform* toFPA; // Transformation from chip to FPA coordinates - psPlaneTransform* fromFPA; // Transformation from FPA to chip coordinates + psPlaneTransform* toFPA; // Transformation from chip to FPA coordinates + psPlaneTransform* fromFPA; // Transformation from FPA to chip coordinates // Information - psMetadata *concepts; // Cache for PS concepts - psMetadata *analysis; // Chip-level analysis metadata - psArray *cells; // The cells (referred to by name) - pmFPA *parent; // Parent FPA - bool valid; // Do we bother about reading and working with this chip? + psMetadata *concepts; // Cache for PS concepts + psMetadata *analysis; // Chip-level analysis metadata + psArray *cells; // The cells (referred to by name) + pmFPA *parent; // Parent FPA + bool valid; // Do we bother about reading and working with this chip? p_pmHDU *hdu; // FITS data } pmChip; @@ -49,35 +49,35 @@ typedef struct { // Offset specifying position on chip - int col0; // Offset from the left of chip. - int row0; // Offset from the bottom of chip. + int col0; // Offset from the left of chip. + int row0; // Offset from the bottom of chip. // Astrometric transformations - psPlaneTransform* toChip; // Transformations from cell to chip coordinates - psPlaneTransform* toFPA; // Transformations from cell to FPA coordinates - psPlaneTransform* toSky; // Transformations from cell to sky coordinates + psPlaneTransform* toChip; // Transformations from cell to chip coordinates + psPlaneTransform* toFPA; // Transformations from cell to FPA coordinates + psPlaneTransform* toSky; // Transformations from cell to sky coordinates // Information - psMetadata *concepts; // Cache for PS concepts - psMetadata *camera; // Camera information - psMetadata *analysis; // Cell-level analysis metadata - psArray *readouts; // The readouts (referred to by number) - pmChip *parent; // Parent chip - bool valid; // Do we bother about reading and working with this cell? - p_pmHDU *hdu; // FITS data + psMetadata *concepts; // Cache for PS concepts + psMetadata *camera; // Camera information + psMetadata *analysis; // Cell-level analysis metadata + psArray *readouts; // The readouts (referred to by number) + pmChip *parent; // Parent chip + bool valid; // Do we bother about reading and working with this cell? + p_pmHDU *hdu; // FITS data } pmCell; typedef struct { // Position on the cell - int col0; // Offset from the left of cell. - int row0; // Offset from the bottom of cell. - int colBins; // Amount of binning in x-dimension and parity (from sign) - int rowBins; // Amount of binning in y-dimension and parity (from sign) + int col0; // Offset from the left of cell. + int row0; // Offset from the bottom of cell. + int colBins; // Amount of binning in x-dimension and parity (from sign) + int rowBins; // Amount of binning in y-dimension and parity (from sign) // Information - psImage *image; // Imaging area of readout - psImage *mask; // Mask for image - psImage *weight; // Weight for image + psImage *image; // Imaging area of readout + psImage *mask; // Mask for image + psImage *weight; // Weight for image #if 0 - psList *bias; // List of bias section (sub-)images + psList *bias; // List of bias section (sub-)images #endif - psMetadata *analysis; // Readout-level analysis metadata - pmCell *parent; // Parent cell + psMetadata *analysis; // Readout-level analysis metadata + pmCell *parent; // Parent cell } pmReadout; Index: /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.c (revision 5786) @@ -14,30 +14,30 @@ psMetadataItem *p_pmFPAConceptGetFromCamera(pmCell *cell, // The cell - const char *concept // Name of concept + const char *concept // Name of concept ) { if (cell) { - psMetadata *camera = cell->camera; // Camera data - // Need the CELL.NAME first, which should be in the cell->concepts since creation of the cell - bool mdStatus = true; // Status of MD lookup - // Finally, the info that we want - psMetadataItem *item = psMetadataLookup(camera, concept); - return item; + psMetadata *camera = cell->camera; // Camera data + // Need the CELL.NAME first, which should be in the cell->concepts since creation of the cell + bool mdStatus = true; // Status of MD lookup + // Finally, the info that we want + psMetadataItem *item = psMetadataLookup(camera, concept); + return item; } return NULL; -} +} psMetadataItem *p_pmFPAConceptGetFromHeader(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept - ) -{ - bool mdStatus = true; // Status of MD lookup + pmChip *chip, // The chip that contains the cell + pmCell *cell, // The cell + const char *concept // Name of concept + ) +{ + bool mdStatus = true; // Status of MD lookup psMetadata *translation = psMetadataLookupMD(&mdStatus, fpa->camera, "TRANSLATION"); // FITS translation if (! mdStatus) { - psError(PS_ERR_IO, false, "Unable to find TRANSLATION in camera configuration.\n"); - return NULL; + psError(PS_ERR_IO, false, "Unable to find TRANSLATION in camera configuration.\n"); + return NULL; } @@ -45,36 +45,36 @@ const char *keyword = psMetadataLookupString(&mdStatus, translation, concept); if (mdStatus && strlen(keyword) > 0) { - // We have a FITS header to look up --- search each level - if (cell && cell->hdu) { - psMetadataItem *cellItem = psMetadataLookup(cell->hdu->header, keyword); - if (cellItem) { - // XXX: Need to clean up before returning - return cellItem; - } - } - - if (chip && chip->hdu) { - psMetadataItem *chipItem = psMetadataLookup(chip->hdu->header, keyword); - if (chipItem) { - // XXX: Need to clean up before returning - return chipItem; - } - } - - if (fpa->hdu) { - psMetadataItem *fpaItem = psMetadataLookup(fpa->hdu->header, keyword); - if (fpaItem) { - // XXX: Need to clean up before returning - return fpaItem; - } - } - - if (fpa->phu) { - psMetadataItem *fpaItem = psMetadataLookup(fpa->phu, keyword); - if (fpaItem) { - // XXX: Need to clean up before returning - return fpaItem; - } - } + // We have a FITS header to look up --- search each level + if (cell && cell->hdu) { + psMetadataItem *cellItem = psMetadataLookup(cell->hdu->header, keyword); + if (cellItem) { + // XXX: Need to clean up before returning + return cellItem; + } + } + + if (chip && chip->hdu) { + psMetadataItem *chipItem = psMetadataLookup(chip->hdu->header, keyword); + if (chipItem) { + // XXX: Need to clean up before returning + return chipItem; + } + } + + if (fpa->hdu) { + psMetadataItem *fpaItem = psMetadataLookup(fpa->hdu->header, keyword); + if (fpaItem) { + // XXX: Need to clean up before returning + return fpaItem; + } + } + + if (fpa->phu) { + psMetadataItem *fpaItem = psMetadataLookup(fpa->phu, keyword); + if (fpaItem) { + // XXX: Need to clean up before returning + return fpaItem; + } + } } @@ -86,51 +86,51 @@ // Look for a default psMetadataItem *p_pmFPAConceptGetFromDefault(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept - ) -{ - bool mdOK = true; // Status of MD lookup + pmChip *chip, // The chip that contains the cell + pmCell *cell, // The cell + const char *concept // Name of concept + ) +{ + bool mdOK = true; // Status of MD lookup psMetadata *defaults = psMetadataLookupMD(&mdOK, fpa->camera, "DEFAULTS"); if (! mdOK) { - psError(PS_ERR_IO, false, "Unable to find DEFAULTS in camera configuration.\n"); - return NULL; + psError(PS_ERR_IO, false, "Unable to find DEFAULTS in camera configuration.\n"); + return NULL; } psMetadataItem *defItem = psMetadataLookup(defaults, concept); if (defItem) { - if (defItem->type == PS_DATA_METADATA) { - // A dependent default - psTrace(__func__, 7, "Evaluating dependent default....\n"); - psMetadata *dependents = defItem->data.V; // The list of dependents - // Find out what it depends on - psString dependName = psStringCopy(concept); - psStringAppend(&dependName, ".DEPEND"); - psString dependsOn = psMetadataLookupString(&mdOK, defaults, dependName); - if (! mdOK) { - psError(PS_ERR_IO, false, "Unable to find %s in camera configuration for dependent default" - " --- ignored\n", dependName); - // XXX: Need to clean up before returning - return NULL; - } - psFree(dependName); - // Find the value of the dependent concept - psMetadataItem *depItem = p_pmFPAConceptGetFromHeader(fpa, chip, cell, dependsOn); - if (! depItem) { - psError(PS_ERR_IO, true, "Unable to find value for %s (required for %s)\n", dependsOn, - concept); - return NULL; - } - if (depItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "Value of %s is not of type string, as required for dependency" - " --- ignored.\n", dependsOn); - } - - defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after - } + if (defItem->type == PS_DATA_METADATA) { + // A dependent default + psTrace(__func__, 7, "Evaluating dependent default....\n"); + psMetadata *dependents = defItem->data.V; // The list of dependents + // Find out what it depends on + psString dependName = psStringCopy(concept); + psStringAppend(&dependName, ".DEPEND"); + psString dependsOn = psMetadataLookupString(&mdOK, defaults, dependName); + if (! mdOK) { + psError(PS_ERR_IO, false, "Unable to find %s in camera configuration for dependent default" + " --- ignored\n", dependName); + // XXX: Need to clean up before returning + return NULL; + } + psFree(dependName); + // Find the value of the dependent concept + psMetadataItem *depItem = p_pmFPAConceptGetFromHeader(fpa, chip, cell, dependsOn); + if (! depItem) { + psError(PS_ERR_IO, true, "Unable to find value for %s (required for %s)\n", dependsOn, + concept); + return NULL; + } + if (depItem->type != PS_DATA_STRING) { + psError(PS_ERR_IO, true, "Value of %s is not of type string, as required for dependency" + " --- ignored.\n", dependsOn); + } + + defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after + } } // XXX: Need to clean up before returning - return defItem; // defItem is either NULL or points to what was desired + return defItem; // defItem is either NULL or points to what was desired } @@ -139,93 +139,93 @@ // XXX: Not tested psMetadataItem *p_pmFPAConceptGetFromDB(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Name of concept + pmChip *chip, // The chip that contains the cell + pmCell *cell, // The cell + psDB *db, // DB handle + const char *concept // Name of concept ) { if (! db) { - // No database initialised - return NULL; - } - - bool mdStatus = true; // Status of MD lookup + // No database initialised + return NULL; + } + + bool mdStatus = true; // Status of MD lookup psMetadata *database = psMetadataLookupMD(&mdStatus, fpa->camera, "DATABASE"); if (! mdStatus) { - // No error, because not everyone needs to use the DB - return NULL; + // No error, because not everyone needs to use the DB + return NULL; } psMetadata *dbLookup = psMetadataLookupMD(&mdStatus, database, concept); if (dbLookup) { - const char *tableName = psMetadataLookupString(&mdStatus, dbLookup, "TABLE"); // Name of the table - const char *colName = psMetadataLookupString(&mdStatus, dbLookup, "COLUMN"); // Name of the column - const char *givenCols = psMetadataLookupString(&mdStatus, dbLookup, "GIVENDBCOL"); // Name of "where" - // columns - const char *givenPS = psMetadataLookupString(&mdStatus, dbLookup, "GIVENPS"); // Values for "where" - // columns - - // Now, need to get the "given"s - if (strlen(givenCols) || strlen(givenPS)) { - psList *cols = papSplit(givenCols, ",;"); // List of column names - psList *values = papSplit(givenPS, ",;"); // List of value names for the columns - psMetadata *selection = psMetadataAlloc(); // The stuff to select in the DB - if (cols->n != values->n) { - psLogMsg(__func__, PS_LOG_WARN, "The GIVENDBCOL and GIVENPS entries for %s do not have " - "the same number of entries --- ignored.\n", concept); - } else { - // Iterators for the lists - psListIterator *colsIter = psListIteratorAlloc(cols, PS_LIST_HEAD, false); - psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false); - char *column = NULL; // Name of the column - while (column = psListGetAndIncrement(colsIter)) { - char *name = psListGetAndIncrement(valuesIter); // Name for the value - if (!strlen(column) || !strlen(name)) { - psLogMsg(__func__, PS_LOG_WARN, "One of the columns or value names for %s is " - " empty --- ignored.\n", concept); - } else { - // Search for the value name - psMetadataItem *item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, name); - if (! item) { - item = p_pmFPAConceptGetFromDefault(fpa, chip, cell, name); - } - if (! item) { - psLogMsg(__func__, PS_LOG_ERROR, "Unable to find the value name %s for DB " - " lookup on %s --- ignored.\n", name, concept); - } else { - // We need to create a new psMetadataItem. I don't think we can't simply hack - // the existing one, since that could conceivably cause memory leaks - psMetadataItem *newItem = psMetadataItemAlloc(concept, item->type, - item->comment, item->data.V); - psMetadataAddItem(selection, newItem, PS_LIST_TAIL, PS_META_REPLACE); - psFree(newItem); - } - } - psFree(name); - psFree(column); - } // Iterating through the columns - psFree(colsIter); - psFree(valuesIter); - - psArray *dbResult = psDBSelectRows(db, tableName, selection, 2); // Lookup result - // Note that we use limit=2 in order to test if there are multiple rows returned - - psMetadataItem *result = NULL; // The final result of the DB lookup - if (dbResult->n == 0) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find any rows in DB for %s --- ignored\n", - concept); - } else { - if (dbResult-> n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Multiple rows returned in DB lookup for %s --- " - " using the first one only.\n", concept); - } - result = (psMetadataItem*)dbResult->data[0]; - } - // XXX: Need to clean up before returning - return result; - } - psFree(cols); - psFree(values); - } + const char *tableName = psMetadataLookupString(&mdStatus, dbLookup, "TABLE"); // Name of the table + const char *colName = psMetadataLookupString(&mdStatus, dbLookup, "COLUMN"); // Name of the column + const char *givenCols = psMetadataLookupString(&mdStatus, dbLookup, "GIVENDBCOL"); // Name of "where" + // columns + const char *givenPS = psMetadataLookupString(&mdStatus, dbLookup, "GIVENPS"); // Values for "where" + // columns + + // Now, need to get the "given"s + if (strlen(givenCols) || strlen(givenPS)) { + psList *cols = papSplit(givenCols, ",;"); // List of column names + psList *values = papSplit(givenPS, ",;"); // List of value names for the columns + psMetadata *selection = psMetadataAlloc(); // The stuff to select in the DB + if (cols->n != values->n) { + psLogMsg(__func__, PS_LOG_WARN, "The GIVENDBCOL and GIVENPS entries for %s do not have " + "the same number of entries --- ignored.\n", concept); + } else { + // Iterators for the lists + psListIterator *colsIter = psListIteratorAlloc(cols, PS_LIST_HEAD, false); + psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false); + char *column = NULL; // Name of the column + while (column = psListGetAndIncrement(colsIter)) { + char *name = psListGetAndIncrement(valuesIter); // Name for the value + if (!strlen(column) || !strlen(name)) { + psLogMsg(__func__, PS_LOG_WARN, "One of the columns or value names for %s is " + " empty --- ignored.\n", concept); + } else { + // Search for the value name + psMetadataItem *item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, name); + if (! item) { + item = p_pmFPAConceptGetFromDefault(fpa, chip, cell, name); + } + if (! item) { + psLogMsg(__func__, PS_LOG_ERROR, "Unable to find the value name %s for DB " + " lookup on %s --- ignored.\n", name, concept); + } else { + // We need to create a new psMetadataItem. I don't think we can't simply hack + // the existing one, since that could conceivably cause memory leaks + psMetadataItem *newItem = psMetadataItemAlloc(concept, item->type, + item->comment, item->data.V); + psMetadataAddItem(selection, newItem, PS_LIST_TAIL, PS_META_REPLACE); + psFree(newItem); + } + } + psFree(name); + psFree(column); + } // Iterating through the columns + psFree(colsIter); + psFree(valuesIter); + + psArray *dbResult = psDBSelectRows(db, tableName, selection, 2); // Lookup result + // Note that we use limit=2 in order to test if there are multiple rows returned + + psMetadataItem *result = NULL; // The final result of the DB lookup + if (dbResult->n == 0) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to find any rows in DB for %s --- ignored\n", + concept); + } else { + if (dbResult-> n > 1) { + psLogMsg(__func__, PS_LOG_WARN, "Multiple rows returned in DB lookup for %s --- " + " using the first one only.\n", concept); + } + result = (psMetadataItem*)dbResult->data[0]; + } + // XXX: Need to clean up before returning + return result; + } + psFree(cols); + psFree(values); + } } // Doing the "given"s. @@ -236,8 +236,8 @@ // Concept lookup psMetadataItem *p_pmFPAConceptGet(pmFPA *fpa, // The FPA - pmChip *chip,// The chip - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Concept name + pmChip *chip,// The chip + pmCell *cell, // The cell + psDB *db, // DB handle + const char *concept // Concept name ) { @@ -245,5 +245,5 @@ psMetadataItem *item = p_pmFPAConceptGetFromCamera(cell, concept); if (! item) { - item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, concept); + item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, concept); } if (! item) { @@ -257,11 +257,11 @@ -void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept +void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ) { @@ -269,35 +269,35 @@ float value = NAN; if (item) { - switch (item->type) { - case PS_DATA_F32: - value = item->data.F32; - break; - case PS_DATA_F64: - // Assume it's OK to truncate to floating point from double - value = (float)item->data.F64; - break; - case PS_DATA_S32: - // Promote to float - value = (float)item->data.S32; - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of floating point type (%x) --- treating as " - "undefined.\n", name, comment, item->type); - } + switch (item->type) { + case PS_DATA_F32: + value = item->data.F32; + break; + case PS_DATA_F64: + // Assume it's OK to truncate to floating point from double + value = (float)item->data.F64; + break; + case PS_DATA_S32: + // Promote to float + value = (float)item->data.S32; + break; + default: + psError(PS_ERR_IO, true, "Concept %s (%s) is not of floating point type (%x) --- treating as " + "undefined.\n", name, comment, item->type); + } } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); + psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); } psTrace(__func__, 7, "Adding %s (%s): %f\n", name, comment, value); - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F32, comment, value); -} - -void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept + psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F32 | PS_META_REPLACE, comment, value); +} + +void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ) { @@ -305,35 +305,35 @@ double value = NAN; if (item) { - switch (item->type) { - case PS_TYPE_F64: - value = item->data.F64; - break; - case PS_TYPE_F32: - // Promote to double - value = (double)item->data.F32; - break; - case PS_TYPE_S32: - // Promote to double - value = (double)item->data.S32; - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of double-precision floating point type (%x) " - "--- treating as undefined.\n", name, comment, item->type); - } + switch (item->type) { + case PS_TYPE_F64: + value = item->data.F64; + break; + case PS_TYPE_F32: + // Promote to double + value = (double)item->data.F32; + break; + case PS_TYPE_S32: + // Promote to double + value = (double)item->data.S32; + break; + default: + psError(PS_ERR_IO, true, "Concept %s (%s) is not of double-precision floating point type (%x) " + "--- treating as undefined.\n", name, comment, item->type); + } } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); + psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); } psTrace(__func__, 7, "Adding %s (%s): %f\n", name, comment, value); - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F64, comment, value); -} - -void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept + psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F64 | PS_META_REPLACE, comment, value); +} + +void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ) { @@ -341,37 +341,37 @@ int value = 0; if (item) { - switch (item->type) { - case PS_TYPE_S32: - value = item->data.S32; - break; - case PS_TYPE_F32: - psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F32 --- converting.\n", - name, comment); - value = (int)item->data.F32; - break; - case PS_TYPE_F64: - psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F64 --- converting.\n", - name, comment); - value = (int)item->data.F64; - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of integer type (%x) --- treating as " - "undefined.\n", name, comment, item->type); - } + switch (item->type) { + case PS_TYPE_S32: + value = item->data.S32; + break; + case PS_TYPE_F32: + psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F32 --- converting.\n", + name, comment); + value = (int)item->data.F32; + break; + case PS_TYPE_F64: + psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F64 --- converting.\n", + name, comment); + value = (int)item->data.F64; + break; + default: + psError(PS_ERR_IO, true, "Concept %s (%s) is not of integer type (%x) --- treating as " + "undefined.\n", name, comment, item->type); + } } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); + psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); } psTrace(__func__, 7, "Adding %s (%s): %d\n", name, comment, value); - - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_S32, comment, value); + + psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_S32 | PS_META_REPLACE, comment, value); } void p_pmFPAConceptGetString(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ) { @@ -379,25 +379,25 @@ psString value = NULL; if (item) { - switch (item->type) { - case PS_DATA_STRING: - value = psMemIncrRefCounter(item->data.V); - break; - case PS_DATA_F32: - psStringAppend(&value, "%f", item->data.F32); - break; - case PS_DATA_S32: - psStringAppend(&value, "%d", item->data.S32); - break; - default: - psError(PS_ERR_IO, true, "Concept %s (%s) is not of string type (%x) --- treating as " - "undefined.\n", name, comment, item->type); - } + switch (item->type) { + case PS_DATA_STRING: + value = psMemIncrRefCounter(item->data.V); + break; + case PS_DATA_F32: + psStringAppend(&value, "%f", item->data.F32); + break; + case PS_DATA_S32: + psStringAppend(&value, "%d", item->data.S32); + break; + default: + psError(PS_ERR_IO, true, "Concept %s (%s) is not of string type (%x) --- treating as " + "undefined.\n", name, comment, item->type); + } } else { - psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); - value = psStringCopy(""); + psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); + value = psStringCopy(""); } psTrace(__func__, 7, "Adding %s (%s): %s\n", name, comment, value); - psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_DATA_STRING, comment, value); + psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_DATA_STRING | PS_META_REPLACE, comment, value); psFree(value); } @@ -409,10 +409,10 @@ // Ingest concepts for the FPA -void pmFPAIngestConcepts(pmFPA *fpa, // The FPA - psDB *db // DB handle +void pmFPAIngestConcepts(pmFPA *fpa, // The FPA + psDB *db // DB handle ) { if (! fpa->concepts) { - fpa->concepts = psMetadataAlloc(); + fpa->concepts = psMetadataAlloc(); } @@ -436,68 +436,68 @@ // FPA.RA { - double ra = NAN; // The RA - psMetadataItem *raItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.RA"); // The FPA.RA item - if (raItem) { - switch (raItem->type) { - case PS_TYPE_F32: - ra = raItem->data.F32; - break; - case PS_TYPE_F64: - ra = raItem->data.F64; - break; - case PS_DATA_STRING: - // Sexagesimal format - { - int big, medium; - float small; - // XXX: Upgrade path is to allow dd:mm.mmm - if (sscanf(raItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && - sscanf(raItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { - psError(PS_ERR_IO, true, "Cannot interpret FPA.RA: %s\n", raItem->data.V); - break; - } - ra = abs(big) + (float)medium/60.0 + small/3600.0; - if (big < 0) { - ra *= -1.0; - } - } - break; - default: - psError(PS_ERR_IO, true, "FPA.RA is of an unexpected type: %x\n", raItem->type); - } - - // How to interpret the RA - const psMetadata *camera = fpa->camera; // Camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString raFormat = psMetadataLookupString(&mdok, formats, "FPA.RA"); - if (mdok && strlen(raFormat) > 0) { - if (strcasecmp(raFormat, "HOURS") == 0) { - ra *= M_PI / 12.0; - } else if (strcasecmp(raFormat, "DEGREES") == 0) { - ra *= M_PI / 180.0; - } else if (strcasecmp(raFormat, "RADIANS") == 0) { - // No action required - } else { - psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.RA in FORMATS (%s) --- assuming" - " HOURS.\n"); - ra *= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FPA.RA in FORMATS --- assuming HOURS.\n"); - ra *= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FORMAT metadata in camera configuration --- " - "assuming format for FPA.RA is HOURS.\n"); - ra *= M_PI / 12.0; - } - } else { - psError(PS_ERR_IO, false, "Couldn't find FPA.RA.\n"); - } - - psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_DATA_F64, - "Right Ascension of the boresight (radians)", ra); + double ra = NAN; // The RA + psMetadataItem *raItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.RA"); // The FPA.RA item + if (raItem) { + switch (raItem->type) { + case PS_TYPE_F32: + ra = raItem->data.F32; + break; + case PS_TYPE_F64: + ra = raItem->data.F64; + break; + case PS_DATA_STRING: + // Sexagesimal format + { + int big, medium; + float small; + // XXX: Upgrade path is to allow dd:mm.mmm + if (sscanf(raItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && + sscanf(raItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { + psError(PS_ERR_IO, true, "Cannot interpret FPA.RA: %s\n", raItem->data.V); + break; + } + ra = abs(big) + (float)medium/60.0 + small/3600.0; + if (big < 0) { + ra *= -1.0; + } + } + break; + default: + psError(PS_ERR_IO, true, "FPA.RA is of an unexpected type: %x\n", raItem->type); + } + + // How to interpret the RA + const psMetadata *camera = fpa->camera; // Camera configuration data + bool mdok = true; // Status of MD lookup + psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); + if (mdok && formats) { + psString raFormat = psMetadataLookupString(&mdok, formats, "FPA.RA"); + if (mdok && strlen(raFormat) > 0) { + if (strcasecmp(raFormat, "HOURS") == 0) { + ra *= M_PI / 12.0; + } else if (strcasecmp(raFormat, "DEGREES") == 0) { + ra *= M_PI / 180.0; + } else if (strcasecmp(raFormat, "RADIANS") == 0) { + // No action required + } else { + psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.RA in FORMATS (%s) --- assuming" + " HOURS.\n"); + ra *= M_PI / 12.0; + } + } else { + psError(PS_ERR_IO, false, "Unable to find FPA.RA in FORMATS --- assuming HOURS.\n"); + ra *= M_PI / 12.0; + } + } else { + psError(PS_ERR_IO, false, "Unable to find FORMAT metadata in camera configuration --- " + "assuming format for FPA.RA is HOURS.\n"); + ra *= M_PI / 12.0; + } + } else { + psError(PS_ERR_IO, false, "Couldn't find FPA.RA.\n"); + } + + psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_DATA_F64 | PS_META_REPLACE, + "Right Ascension of the boresight (radians)", ra); } @@ -505,68 +505,68 @@ // FPA.DEC { - double dec = NAN; // The DEC - psMetadataItem *decItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.DEC"); // The FPA.DEC item - if (decItem) { - switch (decItem->type) { - case PS_TYPE_F32: - dec = decItem->data.F32; - break; - case PS_TYPE_F64: - dec = decItem->data.F64; - break; - case PS_DATA_STRING: - // Sexagesimal format - { - int big, medium; - float small; - // XXX: Upgrade path is to allow dd:mm.mmm - if (sscanf(decItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && - sscanf(decItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { - psError(PS_ERR_IO, true, "Cannot interpret FPA.DEC: %s\n", decItem->data.V); - break; - } - dec = abs(big) + (float)medium/60.0 + small/3600.0; - if (big < 0) { - dec *= -1.0; - } - } - break; - default: - psError(PS_ERR_IO, true, "FPA.DEC is of an unexpected type: %x\n", decItem->type); - } - - // How to interpret the DEC - const psMetadata *camera = fpa->camera; // Camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString decFormat = psMetadataLookupString(&mdok, formats, "FPA.DEC"); - if (mdok && strlen(decFormat) > 0) { - if (strcasecmp(decFormat, "HOURS") == 0) { - dec *= M_PI / 12.0; - } else if (strcasecmp(decFormat, "DEGREES") == 0) { - dec *= M_PI / 180.0; - } else if (strcasecmp(decFormat, "RADIANS") == 0) { - // No action required - } else { - psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.DEC in FORMATS (%s) --- " - "assuming DEGREES.\n"); - dec *= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FPA.DEC in FORMATS --- assuming DEGREES.\n"); - dec *= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Unable to find FORMATS metadata in camera configuration --- " - "assuming format for FPA.DEC is DEGREES.\n"); - dec *= M_PI / 180.0; - } - } else { - psError(PS_ERR_IO, false, "Couldn't find FPA.DEC.\n"); - } - - psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_DATA_F64, - "Declination of the boresight (radians)", dec); + double dec = NAN; // The DEC + psMetadataItem *decItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.DEC"); // The FPA.DEC item + if (decItem) { + switch (decItem->type) { + case PS_TYPE_F32: + dec = decItem->data.F32; + break; + case PS_TYPE_F64: + dec = decItem->data.F64; + break; + case PS_DATA_STRING: + // Sexagesimal format + { + int big, medium; + float small; + // XXX: Upgrade path is to allow dd:mm.mmm + if (sscanf(decItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && + sscanf(decItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { + psError(PS_ERR_IO, true, "Cannot interpret FPA.DEC: %s\n", decItem->data.V); + break; + } + dec = abs(big) + (float)medium/60.0 + small/3600.0; + if (big < 0) { + dec *= -1.0; + } + } + break; + default: + psError(PS_ERR_IO, true, "FPA.DEC is of an unexpected type: %x\n", decItem->type); + } + + // How to interpret the DEC + const psMetadata *camera = fpa->camera; // Camera configuration data + bool mdok = true; // Status of MD lookup + psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); + if (mdok && formats) { + psString decFormat = psMetadataLookupString(&mdok, formats, "FPA.DEC"); + if (mdok && strlen(decFormat) > 0) { + if (strcasecmp(decFormat, "HOURS") == 0) { + dec *= M_PI / 12.0; + } else if (strcasecmp(decFormat, "DEGREES") == 0) { + dec *= M_PI / 180.0; + } else if (strcasecmp(decFormat, "RADIANS") == 0) { + // No action required + } else { + psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.DEC in FORMATS (%s) --- " + "assuming DEGREES.\n"); + dec *= M_PI / 180.0; + } + } else { + psError(PS_ERR_IO, false, "Unable to find FPA.DEC in FORMATS --- assuming DEGREES.\n"); + dec *= M_PI / 180.0; + } + } else { + psError(PS_ERR_IO, false, "Unable to find FORMATS metadata in camera configuration --- " + "assuming format for FPA.DEC is DEGREES.\n"); + dec *= M_PI / 180.0; + } + } else { + psError(PS_ERR_IO, false, "Couldn't find FPA.DEC.\n"); + } + + psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_DATA_F64 | PS_META_REPLACE, + "Declination of the boresight (radians)", dec); } @@ -577,12 +577,12 @@ // Ingest concepts for the chip -bool pmChipIngestConcepts(pmChip *chip, // The chip - psDB *db // DB handle - ) -{ - pmFPA *fpa = chip->parent; // The parent FPA +bool pmChipIngestConcepts(pmChip *chip, // The chip + psDB *db // DB handle + ) +{ + pmFPA *fpa = chip->parent; // The parent FPA if (! chip->concepts) { - chip->concepts = psMetadataAlloc(); + chip->concepts = psMetadataAlloc(); } @@ -591,42 +591,430 @@ // Pau. } + + +// Add corrective to a position --- in case the user wants FORTRAN indexing (like the FITS standard...) +static bool correctPosition(pmFPA *fpa, // FPA, contains the camera configuration + pmCell *cell, // Cell containing the concept to correct + const char *conceptName // Name of concept to correct + ) +{ + bool mdok = false; // Result of MD lookup + psMetadata *formats = psMetadataLookupMD(&mdok, fpa->camera, "FORMATS"); + if (mdok && formats) { + psString format = psMetadataLookupString(&mdok, formats, conceptName); + if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { + psMetadataItem *valueItem = psMetadataLookup(cell->concepts, conceptName); + valueItem->data.S32 -= 1; + } + } + return true; +} + +bool p_pmCellIngestConcept(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + const char *concept // Name of the concept + ) +{ + if (strcmp(concept, "CELL.GAIN") == 0) { + p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.GAIN", "CCD gain (e/count)"); + } else if (strcmp(concept, "CELL.READNOISE") == 0) { + p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.READNOISE", "CCD read noise (e)"); + } else if (strcmp(concept, "CELL.SATURATION") == 0) { + p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.SATURATION", + "Saturation level (ADU)"); + } else if (strcmp(concept, "CELL.BAD") == 0) { + p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.BAD", "Bad level (ADU)"); + } else if (strcmp(concept, "CELL.XPARITY") == 0) { + p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.XPARITY", + "Orientation in x compared to the rest of the FPA"); + } else if (strcmp(concept, "CELL.YPARITY") == 0) { + p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.YPARITY", + "Orientation in y compared to the rest of the FPA"); + } else if (strcmp(concept, "CELL.READDIR") == 0) { + p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.READDIR", + "Read direction: 1=row, 2=col"); + } else if (strcmp(concept, "CELL.EXPOSURE") == 0) { // used to be READOUT.EXPOSURE + p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.EXPOSURE", "Exposure time (sec)"); + } else if (strcmp(concept, "CELL.DARKTIME") == 0) { // used to be READOUT.DARKTIME + p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.DARKTIME", + "Time since CCD flush (sec)"); + // These take some extra work + } else if (strcmp(concept, "CELL.TRIMSEC") == 0) { + psRegion *trimsec = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed by value) + + psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC"); + if (! secItem) { + psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.\n"); + *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); + } else if (secItem->type != PS_DATA_STRING) { + psError(PS_ERR_IO, true, "CELL.TRIMSEC is not of type STR (%x)\n", secItem->type); + *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); + } else { + psString section = secItem->data.V; // The section string + + psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC.SOURCE"); + if (! sourceItem) { + psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.SOURCE.\n"); + *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); + } else if (sourceItem->type != PS_DATA_STRING) { + psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); + *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); + } else { + psString source = sourceItem->data.V; // The source string + + if (strcasecmp(source, "VALUE") == 0) { + *trimsec = psRegionFromString(section); + } else if (strcasecmp(source, "HEADER") == 0) { + psMetadata *header = NULL; // The FITS header + if (cell->hdu) { + header = cell->hdu->header; + } else if (chip->hdu) { + header = chip->hdu->header; + } else if (fpa->hdu) { + header = fpa->hdu->header; + } + if (! header) { + psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); + *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); + } else { + bool mdok = true; // Status of MD lookup + psString secValue = psMetadataLookupString(&mdok, header, section); + if (! mdok || ! secValue) { + psError(PS_ERR_IO, false, "Unable to locate header %s\n", section); + *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); + } else { + *trimsec = psRegionFromString(secValue); + } + } + } else { + psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE (%s) is not HEADER or VALUE --- trying " + "VALUE.\n", source); + *trimsec = psRegionFromString(section); + } // Value of CELL.TRIMSEC.SOURCE + } // Looking up CELL.TRIMSEC.SOURCE + } // Looking up CELL.TRIMSEC + + psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE, + "Trim section", trimsec); + psFree(trimsec); + + } else if (strcmp(concept, "CELL.BIASSEC") == 0) { + psList *biassecs = psListAlloc(NULL); // List of bias sections + + psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC"); + if (! secItem) { + psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.\n"); + } else if (secItem->type != PS_DATA_STRING) { + psError(PS_ERR_IO, true, "CELL.BIASSEC is not of type STR (%x)\n", secItem->type); + } else { + psString sections = secItem->data.V; // The section string + + psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC.SOURCE"); + if (! sourceItem) { + psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.SOURCE.\n"); + } else if (sourceItem->type != PS_DATA_STRING) { + psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); + } else { + psString source = sourceItem->data.V; // The source string + + psList *secList = papSplit(sections, " ;"); // List of sections + psListIterator *secIter = psListIteratorAlloc(secList, PS_LIST_HEAD, false); // Iterator over + // sections + psString aSection = NULL; // A section from the list + while (aSection = psListGetAndIncrement(secIter)) { + psRegion *region = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed + // by value) + + if (strcasecmp(source, "VALUE") == 0) { + *region = psRegionFromString(aSection); + } else if (strcasecmp(source, "HEADER") == 0) { + psMetadata *header = NULL; // The FITS header + if (cell->hdu) { + header = cell->hdu->header; + } else if (chip->hdu) { + header = chip->hdu->header; + } else if (fpa->hdu) { + header = fpa->hdu->header; + } + if (! header) { + psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); + *region = psRegionSet(0.0,0.0,0.0,0.0); + } else { + bool mdok = true; // Status of MD lookup + psString secValue = psMetadataLookupString(&mdok, header, aSection); + if (! mdok || ! secValue) { + psError(PS_ERR_IO, false, "Unable to locate header %s\n", aSection); + *region = psRegionSet(0.0,0.0,0.0,0.0); + } else { + *region = psRegionFromString(secValue); + } + } + } else { + psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE (%s) is not HEADER or VALUE --- trying " + "VALUE.\n", source); + *region = psRegionFromString(aSection); + } // Value of CELL.BIASSEC.SOURCE + + psListAdd(biassecs, PS_LIST_TAIL, region); + psFree(region); + } // Iterating over multiple sections + psFree(secIter); + psFree(secList); + } // Looking up CELL.BIASSEC.SOURCE + } // Looking up CELL.BIASSEC + + psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE, + "Bias sections", biassecs); + psFree(biassecs); + + } else if (strcmp(concept, "CELL.XBIN") == 0) { + int xBin = 1; // Binning factor in x + psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.XBIN"); + if (! binItem) { + psError(PS_ERR_IO, false, "Couldn't find CELL.XBIN.\n"); + } else if (binItem->type == PS_DATA_STRING) { + psString binString = binItem->data.V; // The string containing the binning + if (sscanf(binString, "%d %*d", &xBin) != 1 && + sscanf(binString, "%d,%*d", &xBin) != 1) { + psError(PS_ERR_IO, true, "Unable to read string to get x binning: %s\n", binString); + } + } else if (binItem->type == PS_TYPE_S32) { + xBin = binItem->data.S32; + } else { + psError(PS_ERR_IO, true, "Note sure how to interpret CELL.XBIN of type %x --- assuming 1.\n", + binItem->type); + } + + psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_TYPE_S32 | PS_META_REPLACE, + "Binning in x", xBin); + + } else if (strcmp(concept, "CELL.YBIN") == 0) { + int yBin = 1; // Binning factor in y + psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.YBIN"); + if (! binItem) { + psError(PS_ERR_IO, false, "Couldn't find CELL.YBIN.\n"); + } else if (binItem->type == PS_DATA_STRING) { + psString binString = binItem->data.V; // The string containing the binning + if (sscanf(binString, "%*d %d", &yBin) != 1 && + sscanf(binString, "%*d,%d", &yBin) != 1) { + psError(PS_ERR_IO, true, "Unable to read string to get y binning: %s\n", binString); + } + } else if (binItem->type == PS_TYPE_S32) { + yBin = binItem->data.S32; + } else { + psError(PS_ERR_IO, true, "Note sure how to interpret CELL.YBIN of type %x --- assuming 1.\n", + binItem->type); + } + + psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_TYPE_S32 | PS_META_REPLACE, + "Binning in y", yBin); + + } else if (strcmp(concept, "CELL.TIME") == 0 || strcmp(concept, "CELL.TIMESYS") == 0) { + // Do CELL.TIME and CELL.TIMESYS together + psTime *time = NULL; // The time + psTimeType timeSys = PS_TIME_UTC; // The time system + + // CELL.TIMESYS + psMetadataItem *sysItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIMESYS"); + if (! sysItem) { + psError(PS_ERR_IO, true, "Couldn't find CELL.TIMESYS --- assuming UTC.\n"); + } else if (sysItem->type != PS_DATA_STRING) { + psError(PS_ERR_IO, true, "CELL.TIMESYS isn't of type STRING --- assuming UTC.\n"); + } else { + psString sys = sysItem->data.V; // The time system string + if (strcasecmp(sys, "TAI") == 0) { + timeSys = PS_TIME_TAI; + } else if (strcasecmp(sys, "UTC") == 0) { + timeSys = PS_TIME_UTC; + } else if (strcasecmp(sys, "UT1") == 0) { + timeSys = PS_TIME_UT1; + } else if (strcasecmp(sys, "TT") == 0) { + timeSys = PS_TIME_TT; + } else { + psError(PS_ERR_IO, true, "Can't interpret CELL.TIMESYS --- assuming UTC.\n"); + } + } + psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIMESYS", PS_TYPE_S32 | PS_META_REPLACE, + "Time system", timeSys); + + psMetadataItem *timeItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIME"); + if (! timeItem) { + psError(PS_ERR_IO, false, "Couldn't find CELL.TIME.\n"); + } else { + // Get format + const psMetadata *camera = fpa->camera; // The camera configuration data + bool mdok = true; // Status of MD lookup + psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); + if (mdok && formats) { + psString timeFormat = psMetadataLookupString(&mdok, formats, "CELL.TIME"); + if (mdok && strlen(timeFormat) > 0) { + switch (timeItem->type) { + case PS_DATA_STRING: + { + psString timeString = timeItem->data.V; // String with the time + if (strcasecmp(timeFormat, "ISO") == 0) { + // timeString contains an ISO time + time = psTimeFromISO(timeString, timeSys); + } else if (strstr(timeFormat, "SEPARATE")) { + // timeString contains headers for the date and time + psMetadata *header = NULL; // The FITS header + if (cell->hdu) { + header = cell->hdu->header; + } else if (chip->hdu) { + header = chip->hdu->header; + } else if (fpa->hdu) { + header = fpa->hdu->header; + } + if (! header) { + psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); + } else { + // Get the headers + char *stuff1 = strpbrk(timeString, " ,;"); + psString dateName = psStringNCopy(timeString, + strlen(timeString) - strlen(stuff1)); + char *stuff2 = strpbrk(stuff1, "abcdefghijklmnopqrstuvwxyz" + "ABCDEFGHIJKLMNOPQRSTUVWXYZ"); + psString timeName = psStringCopy(stuff2); + + bool mdok = true; // Status of MD lookup + psString dateString = psMetadataLookupString(&mdok, header, dateName); + psFree(dateName); + int day = 0, month = 0, year = 0; + if (sscanf(dateString, "%d-%d-%d", &day, &month, &year) != 3 && + sscanf(dateString, "%d/%d/%d", &day, &month, &year) != 3) { + psError(PS_ERR_IO, true, "Unable to read date: %s\n", dateString); + } else { + if (strstr(timeFormat, "BACKWARDS")) { + int temp = day; + day = year; + year = temp; + } + if (strstr(timeFormat, "PRE2000") || year < 2000) { + year += 2000; + } + + psMetadataItem *timeItem = psMetadataLookup(header, timeName); + if (! timeItem) { + psError(PS_ERR_IO, false, "Unable to find time header: %s\n", + timeName); + } else if (timeItem->type == PS_DATA_STRING) { + // Time is a string, in the usual way: + psStringAppend(&dateString, "T%s", timeItem->data.V); + } else { + // Assume that time is specified in Second of Day + double seconds = NAN; + switch (timeItem->type) { + case PS_TYPE_S32: + seconds = timeItem->data.S32; + break; + case PS_TYPE_F32: + seconds = timeItem->data.F32; + break; + case PS_TYPE_F64: + seconds = timeItem->data.F64; + break; + default: + psError(PS_ERR_IO, true, "Time header (%s) is not of an " + "expected type: %x\n", timeName, timeItem->type); + } + // Now print to timeString as "hh:mm:ss.ss" + int hours = seconds / 3600; + seconds -= (double)hours * 3600.0; + int minutes = seconds / 60; + seconds -= (double)minutes * 60.0; + psStringAppend(&dateString, "T%02d:%02d:%02f", hours, minutes, + seconds); + } + time = psTimeFromISO(dateString, timeSys); + } // Reading date and time + psFree(timeName); + } // Reading headers + } else { + psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%s) --- trying " + "ISO\n", timeString); + time = psTimeFromISO(timeString, timeSys); + } // Interpreting the time string + } + break; + case PS_TYPE_F32: + { + double timeValue = (double)timeItem->data.F32; + if (strcasecmp(timeFormat, "JD") == 0) { + time = psTimeFromJD(timeValue); + } else if (strcasecmp(timeFormat, "MJD") == 0) { + time = psTimeFromMJD(timeValue); + } else { + psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " + "JD\n", timeValue); + time = psTimeFromJD(timeValue); + } + } + break; + case PS_TYPE_F64: + { + double timeValue = (double)timeItem->data.F64; + if (strcasecmp(timeFormat, "JD") == 0) { + time = psTimeFromJD(timeValue); + } else if (strcasecmp(timeFormat, "MJD") == 0) { + time = psTimeFromMJD(timeValue); + } else { + psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " + "JD\n", timeValue); + time = psTimeFromJD(timeValue); + } + } + break; + default: + psError(PS_ERR_IO, true, "Unable to parse CELL.TIME.\n"); + } + } else { + psError(PS_ERR_IO, false, "Unable to find CELL.TIME in FORMATS.\n"); + } // Getting the format + } else { + psError(PS_ERR_IO, false, "Unable to find FORMATS in camera configuration.\n"); + } // Getting the formats + } // Getting CELL.TIME + + psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIME", PS_DATA_TIME | PS_META_REPLACE, + "Time of exposure", time); + psFree(time); + + // These are new and experimental concepts: CELL.X0 and CELL.Y0 + } else if (strcmp(concept, "CELL.X0") == 0) { + p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.X0", "Position of (0,0) on the chip"); + correctPosition(fpa, cell, "CELL.X0"); + } else if (strcmp(concept, "CELL.Y0") == 0) { + p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.Y0", "Position of (0,0) on the chip"); + correctPosition(fpa, cell, "CELL.Y0"); + } + + return true; +} + // Ingest concepts for the cell bool pmCellIngestConcepts(pmCell *cell, // The cell - psDB *db // DB handle - ) -{ - pmChip *chip = cell->parent; // The parent chip - pmFPA *fpa = chip->parent; // The parent FPA + psDB *db // DB handle + ) +{ + pmChip *chip = cell->parent; // The parent chip + pmFPA *fpa = chip->parent; // The parent FPA if (! cell->concepts) { - cell->concepts = psMetadataAlloc(); + cell->concepts = psMetadataAlloc(); } // CELL.NAME --- added by pmFPAConstruct - // CELL.GAIN - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.GAIN", "CCD gain (e/count)"); - - // CELL.READNOISE - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.READNOISE", "CCD read noise (e)"); - - // CELL.SATURATION - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.SATURATION", "Saturation level (ADU)"); - - // CELL.BAD - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.BAD", "Bad level (ADU)"); - - // CELL.XPARITY - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.XPARITY", "Orientation in x compared to the " - "rest of the FPA"); - - // CELL.YPARITY - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.YPARITY", "Orientation in y compared to the " - "rest of the FPA"); - - // CELL.READDIR - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.READDIR", "Read direction: 1=row, 2=col"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.GAIN"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.READNOISE"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.SATURATION"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.BAD"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.XPARITY"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.YPARITY"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.READDIR"); // These used to be pmReadoutGetExposure and pmReadoutGetDarkTime, but that doesn't really make sense at @@ -634,383 +1022,16 @@ // REALLY derived? They're not in the FITS headers, because a readout is a plane in a 3D image. We'll // have to dream up some additional suffix to specify these, but for now.... - - // CELL.EXPOSURE (used to be READOUT.EXPOSURE) - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.EXPOSURE", "Exposure time (sec)"); - - // CELL.DARKTIME (used to be READOUT.DARKTIME) - p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.DARKTIME", "Time since CCD flush (sec)"); - - // These take some extra work - - // CELL.TRIMSEC - { - psRegion *trimsec = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed by value) - - psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC"); - if (! secItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.\n"); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else if (secItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.TRIMSEC is not of type STR (%x)\n", secItem->type); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - psString section = secItem->data.V; // The section string - - psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC.SOURCE"); - if (! sourceItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.SOURCE.\n"); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else if (sourceItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - psString source = sourceItem->data.V; // The source string - - if (strcasecmp(source, "VALUE") == 0) { - *trimsec = psRegionFromString(section); - } else if (strcasecmp(source, "HEADER") == 0) { - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - bool mdok = true; // Status of MD lookup - psString secValue = psMetadataLookupString(&mdok, header, section); - if (! mdok || ! secValue) { - psError(PS_ERR_IO, false, "Unable to locate header %s\n", section); - *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); - } else { - *trimsec = psRegionFromString(secValue); - } - } - } else { - psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE (%s) is not HEADER or VALUE --- trying " - "VALUE.\n", source); - *trimsec = psRegionFromString(section); - } // Value of CELL.TRIMSEC.SOURCE - } // Looking up CELL.TRIMSEC.SOURCE - } // Looking up CELL.TRIMSEC - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN, - "Trim section", trimsec); - psFree(trimsec); - } - - // CELL.BIASSEC - { - psList *biassecs = psListAlloc(NULL); // List of bias sections - - psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC"); - if (! secItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.\n"); - } else if (secItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.BIASSEC is not of type STR (%x)\n", secItem->type); - } else { - psString sections = secItem->data.V; // The section string - - psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC.SOURCE"); - if (! sourceItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.SOURCE.\n"); - } else if (sourceItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); - } else { - psString source = sourceItem->data.V; // The source string - - psList *secList = papSplit(sections, " ;"); // List of sections - psListIterator *secIter = psListIteratorAlloc(secList, PS_LIST_HEAD, false); // Iterator over - // sections - psString aSection = NULL; // A section from the list - while (aSection = psListGetAndIncrement(secIter)) { - psRegion *region = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed - // by value) - - if (strcasecmp(source, "VALUE") == 0) { - *region = psRegionFromString(aSection); - } else if (strcasecmp(source, "HEADER") == 0) { - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - *region = psRegionSet(0.0,0.0,0.0,0.0); - } else { - bool mdok = true; // Status of MD lookup - psString secValue = psMetadataLookupString(&mdok, header, aSection); - if (! mdok || ! secValue) { - psError(PS_ERR_IO, false, "Unable to locate header %s\n", aSection); - *region = psRegionSet(0.0,0.0,0.0,0.0); - } else { - *region = psRegionFromString(secValue); - } - } - } else { - psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE (%s) is not HEADER or VALUE --- trying " - "VALUE.\n", source); - *region = psRegionFromString(aSection); - } // Value of CELL.BIASSEC.SOURCE - - psListAdd(biassecs, PS_LIST_TAIL, region); - psFree(region); - } // Iterating over multiple sections - psFree(secIter); - psFree(secList); - } // Looking up CELL.BIASSEC.SOURCE - } // Looking up CELL.BIASSEC - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST, "Bias sections", biassecs); - psFree(biassecs); - } - - // CELL.XBIN - { - int xBin = 1; // Binning factor in x - psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.XBIN"); - if (! binItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.XBIN.\n"); - } else if (binItem->type == PS_DATA_STRING) { - psString binString = binItem->data.V; // The string containing the binning - if (sscanf(binString, "%d %*d", &xBin) != 1 && - sscanf(binString, "%d,%*d", &xBin) != 1) { - psError(PS_ERR_IO, true, "Unable to read string to get x binning: %s\n", binString); - } - } else if (binItem->type == PS_TYPE_S32) { - xBin = binItem->data.S32; - } else { - psError(PS_ERR_IO, true, "Note sure how to interpret CELL.XBIN of type %x --- assuming 1.\n", - binItem->type); - } - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_TYPE_S32, "Binning in x", xBin); - } - - // CELL.XBIN - { - int yBin = 1; // Binning factor in y - psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.YBIN"); - if (! binItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.YBIN.\n"); - } else if (binItem->type == PS_DATA_STRING) { - psString binString = binItem->data.V; // The string containing the binning - if (sscanf(binString, "%*d %d", &yBin) != 1 && - sscanf(binString, "%*d,%d", &yBin) != 1) { - psError(PS_ERR_IO, true, "Unable to read string to get y binning: %s\n", binString); - } - } else if (binItem->type == PS_TYPE_S32) { - yBin = binItem->data.S32; - } else { - psError(PS_ERR_IO, true, "Note sure how to interpret CELL.YBIN of type %x --- assuming 1.\n", - binItem->type); - } - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_TYPE_S32, "Binning in y", yBin); - } - - // CELL.TIME and CELL.TIMESYS - { - psTime *time = NULL; // The time - psTimeType timeSys = PS_TIME_UTC; // The time system - - // CELL.TIMESYS - psMetadataItem *sysItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIMESYS"); - if (! sysItem) { - psError(PS_ERR_IO, true, "Couldn't find CELL.TIMESYS --- assuming UTC.\n"); - } else if (sysItem->type != PS_DATA_STRING) { - psError(PS_ERR_IO, true, "CELL.TIMESYS isn't of type STRING --- assuming UTC.\n"); - } else { - psString sys = sysItem->data.V; // The time system string - if (strcasecmp(sys, "TAI") == 0) { - timeSys = PS_TIME_TAI; - } else if (strcasecmp(sys, "UTC") == 0) { - timeSys = PS_TIME_UTC; - } else if (strcasecmp(sys, "UT1") == 0) { - timeSys = PS_TIME_UT1; - } else if (strcasecmp(sys, "TT") == 0) { - timeSys = PS_TIME_TT; - } else { - psError(PS_ERR_IO, true, "Can't interpret CELL.TIMESYS --- assuming UTC.\n"); - } - } - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIMESYS", PS_TYPE_S32, "Time system", timeSys); - - psMetadataItem *timeItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIME"); - if (! timeItem) { - psError(PS_ERR_IO, false, "Couldn't find CELL.TIME.\n"); - } else { - // Get format - const psMetadata *camera = fpa->camera; // The camera configuration data - bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString timeFormat = psMetadataLookupString(&mdok, formats, "CELL.TIME"); - if (mdok && strlen(timeFormat) > 0) { - switch (timeItem->type) { - case PS_DATA_STRING: - { - psString timeString = timeItem->data.V; // String with the time - if (strcasecmp(timeFormat, "ISO") == 0) { - // timeString contains an ISO time - time = psTimeFromISO(timeString, timeSys); - } else if (strstr(timeFormat, "SEPARATE")) { - // timeString contains headers for the date and time - psMetadata *header = NULL; // The FITS header - if (cell->hdu) { - header = cell->hdu->header; - } else if (chip->hdu) { - header = chip->hdu->header; - } else if (fpa->hdu) { - header = fpa->hdu->header; - } - if (! header) { - psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); - } else { - // Get the headers - char *stuff1 = strpbrk(timeString, " ,;"); - psString dateName = psStringNCopy(timeString, - strlen(timeString) - strlen(stuff1)); - char *stuff2 = strpbrk(stuff1, - "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"); - psString timeName = psStringCopy(stuff2); - - bool mdok = true; // Status of MD lookup - psString dateString = psMetadataLookupString(&mdok, header, dateName); - psFree(dateName); - int day = 0, month = 0, year = 0; - if (sscanf(dateString, "%d-%d-%d", &day, &month, &year) != 3 && - sscanf(dateString, "%d/%d/%d", &day, &month, &year) != 3) { - psError(PS_ERR_IO, true, "Unable to read date: %s\n", dateString); - } else { - if (strstr(timeFormat, "BACKWARDS")) { - int temp = day; - day = year; - year = temp; - } - if (strstr(timeFormat, "PRE2000") || year < 2000) { - year += 2000; - } - - psMetadataItem *timeItem = psMetadataLookup(header, timeName); - if (! timeItem) { - psError(PS_ERR_IO, false, "Unable to find time header: %s\n", - timeName); - } else if (timeItem->type == PS_DATA_STRING) { - // Time is a string, in the usual way: - psStringAppend(&dateString, "T%s", timeItem->data.V); - } else { - // Assume that time is specified in Second of Day - double seconds = NAN; - switch (timeItem->type) { - case PS_TYPE_S32: - seconds = timeItem->data.S32; - break; - case PS_TYPE_F32: - seconds = timeItem->data.F32; - break; - case PS_TYPE_F64: - seconds = timeItem->data.F64; - break; - default: - psError(PS_ERR_IO, true, "Time header (%s) is not of an " - "expected type: %x\n", timeName, timeItem->type); - } - // Now print to timeString as "hh:mm:ss.ss" - int hours = seconds / 3600; - seconds -= (double)hours * 3600.0; - int minutes = seconds / 60; - seconds -= (double)minutes * 60.0; - psStringAppend(&dateString, "T%02d:%02d:%02f", hours, minutes, - seconds); - } - time = psTimeFromISO(dateString, timeSys); - } // Reading date and time - psFree(timeName); - } // Reading headers - } else { - psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%s) --- trying " - "ISO\n", timeString); - time = psTimeFromISO(timeString, timeSys); - } // Interpreting the time string - } - break; - case PS_TYPE_F32: - { - double timeValue = (double)timeItem->data.F32; - if (strcasecmp(timeFormat, "JD") == 0) { - time = psTimeFromJD(timeValue); - } else if (strcasecmp(timeFormat, "MJD") == 0) { - time = psTimeFromMJD(timeValue); - } else { - psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " - "JD\n", timeValue); - time = psTimeFromJD(timeValue); - } - } - break; - case PS_TYPE_F64: - { - double timeValue = (double)timeItem->data.F64; - if (strcasecmp(timeFormat, "JD") == 0) { - time = psTimeFromJD(timeValue); - } else if (strcasecmp(timeFormat, "MJD") == 0) { - time = psTimeFromMJD(timeValue); - } else { - psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " - "JD\n", timeValue); - time = psTimeFromJD(timeValue); - } - } - break; - default: - psError(PS_ERR_IO, true, "Unable to parse CELL.TIME.\n"); - } - } else { - psError(PS_ERR_IO, false, "Unable to find CELL.TIME in FORMATS.\n"); - } // Getting the format - } else { - psError(PS_ERR_IO, false, "Unable to find FORMATS in camera configuration.\n"); - } // Getting the formats - } // Getting CELL.TIME - - psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIME", PS_DATA_UNKNOWN, "Time of exposure", time); - psFree(time); - } - - // These are new and experimental concepts: CELL.X0 and CELL.Y0 - - // CELL.X0 - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.X0", "Position of (0,0) on the chip "); - // CELL.Y0 - p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.Y0", "Position of (0,0) on the chip"); - // Add corrective - { - const psMetadata *camera = fpa->camera; - bool mdok = false; // Result of MD lookup - psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); - if (mdok && formats) { - psString format = psMetadataLookupString(&mdok, formats, "CELL.X0"); - if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { - psMetadataItem *cellx0 = psMetadataLookup(cell->concepts, "CELL.X0"); - cellx0->data.S32 -= 1; - } - format = psMetadataLookupString(&mdok, formats, "CELL.Y0"); - if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { - psMetadataItem *celly0 = psMetadataLookup(cell->concepts, "CELL.Y0"); - celly0->data.S32 -= 1; - } - } - } - - // Pau. + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.EXPOSURE"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.DARKTIME"); + + + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.TRIMSEC"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.BIASSEC"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.TIME"); // This also does CELL.TIMESYS + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.XBIN"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.YBIN"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.X0"); + p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.Y0"); + } @@ -1018,16 +1039,16 @@ // Retrieve a concept psMetadataItem *pmCellGetConcept(pmCell *cell, // The cell - const char *concept // The concept + const char *concept // The concept ) { psMetadataItem *item = psMetadataLookup(cell->concepts, concept); if (! item) { - pmChip *chip = cell->parent; - item = psMetadataLookup(chip->concepts, concept); - if (! item) { - pmFPA *fpa = chip->parent; - item = psMetadataLookup(fpa->concepts, concept); - } - } - return item; // item is either NULL or is what we want. -} + pmChip *chip = cell->parent; + item = psMetadataLookup(chip->concepts, concept); + if (! item) { + pmFPA *fpa = chip->parent; + item = psMetadataLookup(fpa->concepts, concept); + } + } + return item; // item is either NULL or is what we want. +} Index: /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.h (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPAConceptsGet.h (revision 5786) @@ -7,78 +7,86 @@ psMetadataItem *p_pmFPAConceptGetFromCamera(pmCell *cell, // The cell - const char *concept // Name of concept + const char *concept // Name of concept ); psMetadataItem *p_pmFPAConceptGetFromHeader(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept + pmChip *chip, // The chip that contains the cell + pmCell *cell, // The cell + const char *concept // Name of concept ); psMetadataItem *p_pmFPAConceptGetFromDefault(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - const char *concept // Name of concept + pmChip *chip, // The chip that contains the cell + pmCell *cell, // The cell + const char *concept // Name of concept ); psMetadataItem *p_pmFPAConceptGetFromDB(pmFPA *fpa, // The FPA that contains the chip - pmChip *chip, // The chip that contains the cell - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Name of concept + pmChip *chip, // The chip that contains the cell + pmCell *cell, // The cell + psDB *db, // DB handle + const char *concept // Name of concept ); psMetadataItem *p_pmFPAConceptGet(pmFPA *fpa, // The FPA - pmChip *chip,// The chip - pmCell *cell, // The cell - psDB *db, // DB handle - const char *concept // Concept name + pmChip *chip,// The chip + pmCell *cell, // The cell + psDB *db, // DB handle + const char *concept // Concept name ); -void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept +void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ); -void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept +void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ); -void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept +void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ); void p_pmFPAConceptGetString(pmFPA *fpa, // The FPA - pmChip *chip, // The chip - pmCell *cell, // The cell - psDB *db, // DB handle - psMetadata *concepts, // The concepts MD - const char *name, // Name of the concept - const char *comment // Comment for concept + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + psMetadata *concepts, // The concepts MD + const char *name, // Name of the concept + const char *comment // Comment for concept ); +// Ingest a single cell concept; this is the workhorse for ingesting cell concepts +bool p_pmCellIngestConcept(pmFPA *fpa, // The FPA + pmChip *chip, // The chip + pmCell *cell, // The cell + psDB *db, // DB handle + const char *concept // Name of the concept + ); + // Ingest concepts for the FPA -void pmFPAIngestConcepts(pmFPA *fpa, // The FPA - psDB *db // DB handle +void pmFPAIngestConcepts(pmFPA *fpa, // The FPA + psDB *db // DB handle ); // Ingest concepts for the chip -bool pmChipIngestConcepts(pmChip *chip, // The chip - psDB *db // DB handle +bool pmChipIngestConcepts(pmChip *chip, // The chip + psDB *db // DB handle ); // Ingest concepts for the cell bool pmCellIngestConcepts(pmCell *cell, // The cell - psDB *db // DB handle + psDB *db // DB handle ); // Retrieve a concept psMetadataItem *pmCellGetConcept(pmCell *cell, // The cell - const char *concept // The concept + const char *concept // The concept ); Index: /trunk/archive/scripts/src/phase2/pmFPAConceptsSet.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConceptsSet.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPAConceptsSet.c (revision 5786) @@ -109,4 +109,5 @@ const char *keyword = psMetadataLookupString(&mdok, translation, concept->name); if (mdok && strlen(keyword) > 0) { + psTrace(__func__, 6, "It's in keyword %s\n", keyword); psMetadataItem *headerItem = NULL; // Item to add to header // XXX: Need to expand range of types @@ -131,31 +132,23 @@ // We have a FITS header to look up --- search each level if (cell && cell->hdu) { - psMetadataItem *cellItem = psMetadataLookup(cell->hdu->header, keyword); - if (cellItem) { - // XXX: Need to clean up before returning - psMetadataAddItem(cell->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } + psTrace(__func__, 7, "Adding to the cell level header...\n"); + psMetadataAddItem(cell->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); + status = true; } else if (chip && chip->hdu) { - psMetadataItem *chipItem = psMetadataLookup(chip->hdu->header, keyword); - if (chipItem) { - // XXX: Need to clean up before returning - psMetadataAddItem(chip->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } + psTrace(__func__, 7, "Adding to the chip level header...\n"); + psMetadataAddItem(chip->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); + status = true; } else if (fpa->hdu) { - psMetadataItem *fpaItem = psMetadataLookup(fpa->hdu->header, keyword); - if (fpaItem) { - // XXX: Need to clean up before returning - psMetadataAddItem(fpa->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } - } else if (fpa->phu) { - psMetadataItem *fpaItem = psMetadataLookup(fpa->phu, keyword); - if (fpaItem) { - // XXX: Need to clean up before returning - psMetadataAddItem(fpa->phu, headerItem, PS_LIST_TAIL, PS_META_REPLACE); - status = true; - } + psTrace(__func__, 7, "Adding to the FPA level header...\n"); + psMetadataAddItem(fpa->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); + status = true; + } else { + // In desperation, add to the PHU --- it HAS to be in the header somewhere! + if (! fpa->phu) { + fpa->phu = psMetadataAlloc(); + } + psTrace(__func__, 7, "Adding to the PHU...\n"); + psMetadataAddItem(fpa->phu, headerItem, PS_LIST_TAIL, PS_META_REPLACE); + status = true; } psFree(headerItem); @@ -182,4 +175,5 @@ psMetadataItem *defItem = psMetadataLookup(defaults, concept->name); + bool status = false; // Result of checking the database if (defItem) { if (defItem->type == PS_DATA_METADATA) { @@ -209,11 +203,17 @@ " --- ignored.\n", dependsOn); } - - defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after + + defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after + } + + status = compareConcepts(defItem, concept); + if (! status) { + psError(PS_ERR_IO, true, "Concept %s is specified by default in the camera configuration, " + "but doesn't match the actual value.\n", concept->name); } } // XXX: Need to clean up before returning - return compareConcepts(defItem, concept); + return status; } @@ -330,12 +330,16 @@ { // Try headers, database, defaults in order + psTrace(__func__, 3, "Trying to set concept %s...\n", concept->name); bool status = setConceptInCamera(cell, concept); // Status for return if (! status) { + psTrace(__func__, 5, "Trying header....\n"); status = setConceptInHeader(fpa, chip, cell, concept); } if (! status) { + psTrace(__func__, 5, "Trying database....\n"); status = setConceptInDB(fpa, chip, cell, db, concept); } if (! status) { + psTrace(__func__, 5, "Checking defaults....\n"); status = setConceptInDefault(fpa, chip, cell, concept); } Index: /trunk/archive/scripts/src/phase2/pmFPAConstruct.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAConstruct.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPAConstruct.c (revision 5786) @@ -67,5 +67,4 @@ const char *phuType = psMetadataLookupString(&mdStatus, camera, "PHU"); // What is the PHU? const char *extType = psMetadataLookupString(&mdStatus, camera, "EXTENSIONS"); // What's in the extns? - if (strcasecmp(phuType, "FPA") == 0) { Index: /trunk/archive/scripts/src/phase2/pmFPAMorph.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAMorph.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPAMorph.c (revision 5786) @@ -14,5 +14,5 @@ // Flip the chip in x static psImage *xFlipImage(psImage *image, // Image to be flipped - int size // Size of image in flip dimension + int size // Size of image in flip dimension ) { @@ -20,19 +20,19 @@ switch (image->type.type) { case PS_TYPE_U16: - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - flipped->data.U16[y][size - x - 1] = image->data.U16[y][x]; - } - } - break; + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + flipped->data.U16[y][size - x - 1] = image->data.U16[y][x]; + } + } + break; case PS_TYPE_F32: - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - flipped->data.F32[y][size - x - 1] = image->data.F32[y][x]; - } - } - break; + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + flipped->data.F32[y][size - x - 1] = image->data.F32[y][x]; + } + } + break; default: - psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type); + psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type); } @@ -42,5 +42,5 @@ // Flip the chip in y static psImage *yFlipImage(psImage *image, // Image to be flipped - int size // Size of image in flip dimension + int size // Size of image in flip dimension ) { @@ -48,12 +48,12 @@ switch (image->type.type) { case PS_TYPE_U16: - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - flipped->data.U16[size - y - 1][x] = image->data.U16[y][x]; - } - } - break; + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + flipped->data.U16[size - y - 1][x] = image->data.U16[y][x]; + } + } + break; default: - psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type); + psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type); } @@ -65,14 +65,14 @@ // Return a list of the region strings and the method for interpreting them static psList *getRegions(psString *method, // Method for evaluating the regions - psString methodValues // The METHOD:VALUES string + psString methodValues // The METHOD:VALUES string ) { char *values = strchr(methodValues, ':'); // The VALUES *method = psStringNCopy(methodValues, strlen(methodValues) - strlen(values)); // The METHOD - psList *regionStrings = NULL; // List of the region strings + psList *regionStrings = NULL; // List of the region strings if (strncmp(*method, "HEADER", 6) == 0 || strncmp(*method, "HD", 2) == 0) { - regionStrings = papSplit(values, ", "); + regionStrings = papSplit(values, ", "); } else if (strncmp(*method, "VALUE", 5) == 0 || strncmp(*method, "VAL", 3) == 0) { - regionStrings = papSplit(values, "; "); + regionStrings = papSplit(values, "; "); } @@ -84,41 +84,43 @@ // Splice an image into another image static psRegion *spliceImage(psImage *target, // Target image onto which to splice - int *x0, int *y0, // Starting coordinates to splice to - psImage *source, // Source image from which to splice - const psRegion *from, // Region to splice from - int readdir, // Read direction - bool xFlip, bool yFlip // Flip x and y axes? - ) + int *x0, int *y0, // Starting coordinates to splice to + psImage *source, // Source image from which to splice + const psRegion *from, // Region to splice from + int readdir, // Read direction + bool xFlip, bool yFlip // Flip x and y axes? + ) { psImage *toSplice = psImageSubset(source, *from); // Subimage, to splice into target if (xFlip) { - int size = 0; // Size of flipped image - if (readdir == 1) { - size = toSplice->numCols; - } else if (readdir == 2) { - size = target->numCols; - } - psImage *temp = xFlipImage(toSplice, size); - psFree(toSplice); - toSplice = temp; + int size = 0; // Size of flipped image + if (readdir == 1) { + size = toSplice->numCols; + } else if (readdir == 2) { + size = target->numCols; + } + psImage *temp = xFlipImage(toSplice, size); + psFree(toSplice); + toSplice = temp; } if (yFlip) { - int size = 0; // Size of flipped image - if (readdir == 1) { - size = target->numRows; - } else if (readdir == 2) { - size = toSplice->numRows; - } - psImage *temp = yFlipImage(toSplice, size); - psFree(toSplice); - toSplice = temp; - } - + int size = 0; // Size of flipped image + if (readdir == 1) { + size = target->numRows; + } else if (readdir == 2) { + size = toSplice->numRows; + } + psImage *temp = yFlipImage(toSplice, size); + psFree(toSplice); + toSplice = temp; + } + + psTrace(__func__, 5, "Splicing %dx%d image in at %d,%d\n", toSplice->numCols, toSplice->numRows, + *x0, *y0); (void)psImageOverlaySection(target, toSplice, *x0, *y0, "="); if (readdir == 1) { - *x0 += toSplice->numCols; + *x0 += toSplice->numCols; } else if (readdir == 2) { - *y0 += toSplice->numRows; + *y0 += toSplice->numRows; } @@ -136,12 +138,12 @@ // Splice bias regions (overscans) -static bool spliceBias(psImage *splice, // Image to which to splice - int *x0, int *y0, // Starting position for splice - const pmCell *in, // Input cell - pmCell *out, // Output cell - int xFlip, int yFlip, // Flip the image? - int readNum, // Number of readout - int readdir // Read direction - ) +static bool spliceBias(psImage *splice, // Image to which to splice + int *x0, int *y0, // Starting position for splice + const pmCell *in, // Input cell + pmCell *out, // Output cell + int xFlip, int yFlip, // Flip the image? + int readNum, // Number of readout + int readdir // Read direction + ) { psArray *readouts = in->readouts; // The readouts; @@ -152,20 +154,30 @@ psList *outBiassecs = psListAlloc(NULL); // List of biassecs for output while (inBiassec = psListGetAndIncrement(inBiassecsIter)) { - psRegion *outBiassec = spliceImage(splice, x0, y0, readout->image, inBiassec, readdir, xFlip, - yFlip); - psListAdd(outBiassecs, PS_LIST_TAIL, outBiassec); + psRegion *outBiassec = spliceImage(splice, x0, y0, readout->image, inBiassec, readdir, xFlip, + yFlip); + psListAdd(outBiassecs, PS_LIST_TAIL, outBiassec); } psFree(inBiassecsIter); psMetadataAdd(out->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE, - "BIASSEC from pmFPAMorph", outBiassecs); - + "BIASSEC from pmFPAMorph", outBiassecs); + return true; } +psImage *p_pmImageMosaic(const psArray *source, // Images to splice in + const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? + const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of + // source images + int xBinTarget, int yBinTarget, // Binning in x and y of target images + const psVector *x0, const psVector *y0 // Offsets for source images on target + ); // Splice a list of cells together +// We need to do the following: +// 1. Copy the contents of the cells (i.e., the readouts) over from the in to the out +// 2. Splice together everything in the bucket for the image that we will write out static psArray *spliceCells(psList *outCells, // List of target cells (required for parity info) - psList *inCells, // List of cells to splice together - bool posDep // Position dependent placement of overscans? + psList *inCells, // List of cells to splice together + bool posDep // Position dependent placement of overscans? ) { @@ -176,130 +188,292 @@ if (inCells->n != outCells->n) { - psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n", - inCells->n, outCells->n); - return NULL; - } - - int numCells = inCells->n; // Number of cells - int readdir = 0; // Read direction for CCD; currently unknown - int numReadouts = 0; // Number of readouts in a cell - psArray *spliced = NULL; // Array of spliced readouts - psElemType type = 0; // Image type (U16, S32, F32, etc) - - psVector *xSize = NULL; // Size of spliced image in x - psVector *ySize = NULL; // Size of spliced image in y + psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n", + inCells->n, outCells->n); + return NULL; + } + + int numCells = inCells->n; // Number of cells + int readdir = 0; // Read direction for CCD; currently unknown + int numReadouts = 0; // Number of readouts in a cell + psArray *spliced = NULL; // Array of spliced readouts + psElemType type = 0; // Image type (U16, S32, F32, etc) psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in x? psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in y? - + psVector *xBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in x + psVector *yBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in y + psVector *x0 = psVectorAlloc(numCells, PS_TYPE_S32); // Offset in x + psVector *y0 = psVectorAlloc(numCells, PS_TYPE_S32); // Offset in y + int xBinOut = 0; // Binning in x for output + int yBinOut = 0; // Binning in y for output // We go through the cells to make sure everything's kosher, and to get the sizes. psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false); // Iterator for cells - pmCell *inCell = NULL; // Cell from list of input cells psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator for cells - int cellNum = 0; // Cell number; + pmCell *inCell = NULL; // Cell from list of input cells + pmCell *outCell = NULL; // Cell from list of output cells + int cellNum = -1; // Cell number + while (inCell = psListGetAndIncrement(inCellsIter) && outCell = psListGetAndIncrement(outCellsIter)) { + cellNum++; + psArray *inReadouts = inCell->readouts; // The readouts comprising the cell + + // Get and check the read direction + int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR"); + psTrace(__func__, 10, "Checking read direction for cell %d: %d\n", cellNum, cellReaddir); + if (cellNum == 0) { + // First run through; set it. + readdir = cellReaddir; + } else if (readdir != cellReaddir) { + psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I " + "don't know how to do that!\n", readdir, cellReaddir); + // XXX Clean up before returning + return NULL; + } + + // Get and check the number of readouts + if (! spliced) { + numReadouts = inReadouts->n; + spliced = psArrayAlloc(numReadouts); + xSize = psVectorAlloc(numReadouts, PS_TYPE_S32); + ySize = psVectorAlloc(numReadouts, PS_TYPE_S32); + for (int i = 0; i < numReadouts; i++) { + xSize->data.S32[i] = 0; + ySize->data.S32[i] = 0; + } + } else if (inReadouts->n != numReadouts) { + psError(PS_ERR_IO, true, "Trying to splice cells with different number of reads (%d vs %d)\n", + numReadouts, readouts->n); + // XXX Clean up before returning + return NULL; + } + + // Get and check the binning --- we can't splice to cells that have different binnings, because + // the spliced image can't be defined. Hence all the "to" cells must have the same binning. + xBinIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XBIN"); + yBinIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YBIN"); + int xBinCheck = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XBIN"); + int yBinCheck = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YBIN"); + if (xBinOut == 0) { + xBinOut = xBinCheck; + } else if (xBinCheck != xBinOut) { + psError(PS_ERR_IO, true, "Trying to splice to cells with different x binnings (%d vs %d)\n", + xBinOut, xBinCheck); + // XXX Clean up before returning + return NULL; + } + if (yBinOut == 0) { + yBinOut = yBinCheck; + } else if (yBinCheck != yBinOut) { + psError(PS_ERR_IO, true, "Trying to splice to cells with different y binnings (%d vs %d)\n", + yBinOut, yBinCheck); + // XXX Clean up before returning + return NULL; + } + + // Copy the readouts over + psArray *outReadouts = psArrayAlloc(inReadouts->n); // Array of readouts for the output cell + for (int i = 0; i < inReadouts->n; i++) { + outReadouts->data[i] = psMemIncrRefCounter(inReadouts->data[i]); + } + + // Get the offsets + x0->data.S32[i] = psMetadataLookupS32(NULL, inCell->concepts, "CELL.X0"); + y0->data.S32[i] = psMetadataLookupS32(NULL, inCell->concepts, "CELL.Y0"); + + // Get and check the parity + int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY"); + int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY"); + int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY"); + int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YPARITY"); + if (xParityIn == 0 || xParityOut == 0 || yParityIn == 0 || yParityOut == 0) { + psError(PS_ERR_IO, false, "Unable to determine parity of cell!\n"); + // XXX Clean up before returning + return NULL; + } + if (xParityIn != xParityOut) { + for (int i = 0; i < outReadouts->n; i++) { + psImage *image = outReadouts->data[i]; + outReadouts->data[i] = xFlipImage(image); + psFree(image); + } + // Correct CELL.X0 for the parity shift + psMetadataItem *x0item = psMetadataLookup(inCell->concepts, "CELL.X0"); + x0item->data.S32 -= xSize->data.S32[cellNum]; + } + if (yParityIn != yParityOut) { + for (int i = 0; i < outReadouts->n; i++) { + psImage *image = outReadouts->data[i]; + outReadouts->data[i] = xFlipImage(image); + psFree(image); + } + // Correct CELL.Y 0 for the parity shift + psMetadataItem *y0item = psMetadataLookup(inCell->concepts, "CELL.Y0"); + y0item->data.S32 -= ySize->data.S32[cellNum]; + } + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // DO I WANT TO DO THIS???? + // Wouldn't it be better to splice all the images together first, and then divvy up the components? + ////////////////////////////////////////////////////////////////////////////////////////////////////// + + + // Manipulate the images, if required + for (int i = 0; i < outReadouts->n; i++) { + psImage *image = manipulate(outReadouts + + psImage *image = psMemIncrRefCounter(outReadouts->data[i]); + if (xParityIn != xParityOut) { + psImage *image = xFlipImage(image); + } + + + // THIS ISN'T SO SIMPLE: we need to SET the output x0,y0 based on where the image actually gets pasted + // Get the offsets + x0->data.S32[i] = psMetadataLookupS32(NULL, outCell->concepts, "CELL.X0"); + y0->data.S32[i] = psMetadataLookupS32(NULL, outCell->concepts, "CELL.Y0"); + + + + +// Splice a list of cells together +static psArray *spliceCells(psList *outCells, // List of target cells (required for parity info) + psList *inCells, // List of cells to splice together + bool posDep // Position dependent placement of overscans? + ) +{ + assert(outCells); + assert(inCells); + + psTrace(__func__, 5, "Splicing together all cells in the bucket.\n"); + + if (inCells->n != outCells->n) { + psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n", + inCells->n, outCells->n); + return NULL; + } + + int numCells = inCells->n; // Number of cells + int readdir = 0; // Read direction for CCD; currently unknown + int numReadouts = 0; // Number of readouts in a cell + psArray *spliced = NULL; // Array of spliced readouts + psElemType type = 0; // Image type (U16, S32, F32, etc) + + psVector *xSize = NULL; // Size of spliced image in x + psVector *ySize = NULL; // Size of spliced image in y + + psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in x? + psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in y? + + + // We go through the cells to make sure everything's kosher, and to get the sizes. + psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false); // Iterator for cells + psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator for cells + pmCell *inCell = NULL; // Cell from list of input cells + int cellNum = 0; // Cell number; while (inCell = psListGetAndIncrement(inCellsIter)) { - psArray *readouts = inCell->readouts; // The readouts comprising the cell - - // Get and check the read direction - int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR"); - if (readdir == 0) { - // Zero means unknown; set it. - readdir = cellReaddir; - } else if (readdir != cellReaddir) { - psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I " - "don't know how to do that!\n", readdir, cellReaddir); - // XXX Clean up before returning - return NULL; - } - - // Get and check the number of readouts - if (! spliced) { - numReadouts = readouts->n; - spliced = psArrayAlloc(numReadouts); - xSize = psVectorAlloc(numReadouts, PS_TYPE_S32); - ySize = psVectorAlloc(numReadouts, PS_TYPE_S32); - for (int i = 0; i < numReadouts; i++) { - xSize->data.S32[i] = 0; - ySize->data.S32[i] = 0; - } - } else if (readouts->n != numReadouts) { - psError(PS_ERR_IO, true, "Trying to splice cells with different number of reads (%d vs %d)\n", - numReadouts, readouts->n); - // XXX Clean up before returning - return NULL; - } - - // Get and check the cell parities - pmCell *outCell = psListGetAndIncrement(outCellsIter); // Corresponding output cell - int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY"); - int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY"); - int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY"); - int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YPARITY"); - if (xParityIn == 0 || xParityOut == 0 || yParityIn == 0 || yParityOut == 0) { - psError(PS_ERR_IO, false, "Unable to determine parity of cell!\n"); - // XXX Clean up before returning - return NULL; - } - - if (yParityIn == yParityOut) { - yFlip->data.U8[cellNum] = 0; - } else { - psTrace(__func__, 7, "Need to flip y.\n"); - yFlip->data.U8[cellNum] = 1; - } - if (xParityIn == xParityOut) { - xFlip->data.U8[cellNum] = 0; - } else { - psTrace(__func__, 7, "Need to flip x.\n"); - xFlip->data.U8[cellNum] = 1; - } - cellNum++; - - // Calculate the sizes of the spliced images - for (int i = 0; i < numReadouts; i++) { - pmReadout *readout = readouts->data[i]; // The readout of interest - psImage *image = readout->image; // The image pixels - - // Check image type - if (type == 0) { - type = image->type.type; - } else if (image->type.type != type) { - psError(PS_ERR_IO, true, "Trying to splice readouts with different image types (%d vs %d)\n", - type, image->type.type); - // XXX Clean up before returning - return NULL; - } - - psRegion *trimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // Trim section - psList *biassecs = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.BIASSEC"); // Bias section - psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator - psRegion *biassec = NULL; // Bias section from list iteration - - // Get the size of the spliced image - if (readdir == 1) { - xSize->data.S32[i] += trimsec->x1 - trimsec->x0; - ySize->data.S32[i] = MAX(ySize->data.S32[i], trimsec->y1 - trimsec->y0); - while (biassec = psListGetAndIncrement(biassecsIter)) { - xSize->data.S32[i] += biassec->x1 - biassec->x0; - ySize->data.S32[i] = MAX(ySize->data.S32[i], biassec->y1 - biassec->y0); - } - } else if (readdir == 2) { - ySize->data.S32[i] += trimsec->y1 - trimsec->y0; - xSize->data.S32[i] = MAX(xSize->data.S32[i], trimsec->x1 - trimsec->x0); - while (biassec = psListGetAndIncrement(biassecsIter)) { - ySize->data.S32[i] += biassec->y1 - biassec->y0; - xSize->data.S32[i] = MAX(ySize->data.S32[i], biassec->x1 - biassec->x0); - } - } else { - psError(PS_ERR_IO, true, "Invalid read direction: %d\n", readdir); - // XXX Clean up before returning - return NULL; - } - psFree(biassecsIter); - } - - // Copy the concepts over - psMetadataCopy(outCell->concepts, inCell->concepts); + psArray *readouts = inCell->readouts; // The readouts comprising the cell + + // Get and check the read direction + int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR"); + psTrace(__func__, 10, "Checking read direction for cell %d: %d\n", cellNum, cellReaddir); + if (cellNum == 0) { + // First run through; set it. + readdir = cellReaddir; + } else if (readdir != cellReaddir) { + psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I " + "don't know how to do that!\n", readdir, cellReaddir); + // XXX Clean up before returning + return NULL; + } + + // Get and check the number of readouts + if (! spliced) { + numReadouts = readouts->n; + spliced = psArrayAlloc(numReadouts); + xSize = psVectorAlloc(numReadouts, PS_TYPE_S32); + ySize = psVectorAlloc(numReadouts, PS_TYPE_S32); + for (int i = 0; i < numReadouts; i++) { + xSize->data.S32[i] = 0; + ySize->data.S32[i] = 0; + } + } else if (readouts->n != numReadouts) { + psError(PS_ERR_IO, true, "Trying to splice cells with different number of reads (%d vs %d)\n", + numReadouts, readouts->n); + // XXX Clean up before returning + return NULL; + } + + // Get and check the cell parities + pmCell *outCell = psListGetAndIncrement(outCellsIter); // Corresponding output cell + int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY"); + int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY"); + int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY"); + int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YPARITY"); + if (xParityIn == 0 || xParityOut == 0 || yParityIn == 0 || yParityOut == 0) { + psError(PS_ERR_IO, false, "Unable to determine parity of cell!\n"); + // XXX Clean up before returning + return NULL; + } + + if (yParityIn == yParityOut) { + yFlip->data.U8[cellNum] = 0; + } else { + psTrace(__func__, 7, "Need to flip y.\n"); + yFlip->data.U8[cellNum] = 1; + } + if (xParityIn == xParityOut) { + xFlip->data.U8[cellNum] = 0; + } else { + psTrace(__func__, 7, "Need to flip x.\n"); + xFlip->data.U8[cellNum] = 1; + } + cellNum++; + + // Calculate the sizes of the spliced images + for (int i = 0; i < numReadouts; i++) { + pmReadout *readout = readouts->data[i]; // The readout of interest + psImage *image = readout->image; // The image pixels + + // Check image type + if (type == 0) { + type = image->type.type; + } else if (image->type.type != type) { + psError(PS_ERR_IO, true, "Trying to splice readouts with different image types (%d vs %d)\n", + type, image->type.type); + // XXX Clean up before returning + return NULL; + } + + psRegion *trimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // Trim section + psList *biassecs = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.BIASSEC"); // Bias section + psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator + psRegion *biassec = NULL; // Bias section from list iteration + + // Get the size of the spliced image + if (readdir == 1) { + psTrace(__func__, 7, "TRIMSEC is %.0fx%.0f\n", trimsec->x1 - trimsec->x0, + trimsec->y1 - trimsec->y0); + xSize->data.S32[i] += trimsec->x1 - trimsec->x0; + ySize->data.S32[i] = MAX(ySize->data.S32[i], trimsec->y1 - trimsec->y0); + while (biassec = psListGetAndIncrement(biassecsIter)) { + psTrace(__func__, 9, "BIASSEC is %.0fx%.0f\n", biassec->x1 - biassec->x0, + biassec->y1 - biassec->y0); + xSize->data.S32[i] += biassec->x1 - biassec->x0; + ySize->data.S32[i] = MAX(ySize->data.S32[i], biassec->y1 - biassec->y0); + } + } else if (readdir == 2) { + ySize->data.S32[i] += trimsec->y1 - trimsec->y0; + xSize->data.S32[i] = MAX(xSize->data.S32[i], trimsec->x1 - trimsec->x0); + while (biassec = psListGetAndIncrement(biassecsIter)) { + ySize->data.S32[i] += biassec->y1 - biassec->y0; + xSize->data.S32[i] = MAX(ySize->data.S32[i], biassec->x1 - biassec->x0); + } + } else { + psError(PS_ERR_IO, true, "Invalid read direction: %d\n", readdir); + // XXX Clean up before returning + return NULL; + } + psFree(biassecsIter); + } } psFree(inCellsIter); @@ -308,11 +482,11 @@ // Make sure all the readouts have the same size - int numRows = ySize->data.S32[0]; // Number of rows for spliced image - int numCols = xSize->data.S32[0]; // Number of columns for spliced image + int numRows = ySize->data.S32[0]; // Number of rows for spliced image + int numCols = xSize->data.S32[0]; // Number of columns for spliced image for (int i = 1; i < numReadouts; i++) { - if (xSize->data.S32[i] != numCols || ySize->data.S32[i] != numRows) { - psError(PS_ERR_IO, true, "Spliced readouts would have different sizes: %dx%d vs %dx%d\n", - numCols, numRows, xSize->data.S32[i], ySize->data.S32[i]); - } + if (xSize->data.S32[i] != numCols || ySize->data.S32[i] != numRows) { + psError(PS_ERR_IO, true, "Spliced readouts would have different sizes: %dx%d vs %dx%d\n", + numCols, numRows, xSize->data.S32[i], ySize->data.S32[i]); + } } @@ -321,57 +495,61 @@ // Now we have done the requisite checks, and know the sizes; we just go through and splice together for (int i = 0; i < numReadouts; i++) { - psImage *splice = psImageAlloc(numCols, numRows, type); // The spliced image - int x0 = 0, y0 = 0; // Position at which the splice begins - - // Position-dependent overscans first - if (posDep) { - for (int cellNum = 0; cellNum < inCells->n / 2; cellNum++) { - pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest - pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest - spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, - readdir); - } - } - - // Then the images - psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false);// Iterator for cells - pmCell *inCell = NULL; // Cell from the list of cells - psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator - int cellNum = 0; // Cell number - while (inCell = psListGetAndIncrement(inCellsIter)) { - psArray *readouts = inCell->readouts; // The readouts comprising the cell - pmReadout *readout = readouts->data[i]; // The specific readout of interest - pmCell *outCell = psListGetAndIncrement(outCellsIter); // Output cell - psRegion *inTrimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // TRIMSEC - psRegion *outTrimsec = spliceImage(splice, &x0, &y0, readout->image, inTrimsec, readdir, - xFlip->data.U8[cellNum], yFlip->data.U8[cellNum]); - - // Update the output TRIMSEC - psMetadataAdd(outCell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE, - "TRIMSEC from pmFPAMorph", outTrimsec); - } - psFree(inCellsIter); - psFree(outCellsIter); - - // Then the biases - if (! posDep) { - // Need to only do half the biases - for (int cellNum = inCells->n / 2; cellNum < inCells->n; cellNum++) { - pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest - pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest - spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, - readdir); - } - } else { - // Need to do all the biases - for (int cellNum = 0; cellNum < inCells->n; cellNum++) { - pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest - pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest - spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, - readdir); - } - } - - spliced->data[i] = splice; + psImage *splice = psImageAlloc(numCols, numRows, type); // The spliced image + int x0 = 0, y0 = 0; // Position at which the splice begins + + // Position-dependent overscans first + if (posDep) { + psTrace(__func__, 8, "Doing first position-dependent biases...\n"); + for (int cellNum = 0; cellNum < inCells->n / 2; cellNum++) { + psTrace(__func__, 9, "Cell %d...\n", cellNum); + pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest + pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest + spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, + readdir); + } + } + + // Then the images + psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false);// Iterator for cells + pmCell *inCell = NULL; // Cell from the list of cells + psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator + int cellNum = 0; // Cell number + while (inCell = psListGetAndIncrement(inCellsIter)) { + psArray *readouts = inCell->readouts; // The readouts comprising the cell + pmReadout *readout = readouts->data[i]; // The specific readout of interest + pmCell *outCell = psListGetAndIncrement(outCellsIter); // Output cell + psRegion *inTrimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // TRIMSEC + psRegion *outTrimsec = spliceImage(splice, &x0, &y0, readout->image, inTrimsec, readdir, + xFlip->data.U8[cellNum], yFlip->data.U8[cellNum]); + + // Update the output TRIMSEC + psMetadataAdd(outCell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE, + "TRIMSEC from pmFPAMorph", outTrimsec); + } + psFree(inCellsIter); + psFree(outCellsIter); + + // Then the biases + if (posDep) { + // Need to only do half the biases + psTrace(__func__, 8, "Doing second position-dependent biases...\n"); + for (int cellNum = inCells->n / 2; cellNum < inCells->n; cellNum++) { + psTrace(__func__, 9, "Cell %d...\n", cellNum); + pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest + pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest + spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, + readdir); + } + } else { + // Need to do all the biases + for (int cellNum = 0; cellNum < inCells->n; cellNum++) { + pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest + pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest + spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, + readdir); + } + } + + spliced->data[i] = splice; } // Iterating over readouts @@ -384,9 +562,9 @@ -bool pmFPAMorph(pmFPA *toFPA, // FPA structure to which to morph - pmFPA *fromFPA, // FPA structure from which to morph - bool positionDependent, // Is the position of the overscan dependent on the position? - int chipNum, // Chip number, in case the scopes are different - int cellNum // Cell number, in case the scopes are different +bool pmFPAMorph(pmFPA *toFPA, // FPA structure to which to morph + pmFPA *fromFPA, // FPA structure from which to morph + bool positionDependent, // Is the position of the overscan dependent on the position? + int chipNum, // Chip number, in case the scopes are different + int cellNum // Cell number, in case the scopes are different ) { @@ -394,94 +572,154 @@ psList *sourceCells = psListAlloc(NULL); // List of source cells - psArray *toChips = toFPA->chips; // Array of chips - psArray *fromChips = fromFPA->chips;// Array of chips + psArray *toChips = toFPA->chips; // Array of chips + psArray *fromChips = fromFPA->chips;// Array of chips + + psTrace(__func__, 1, "Copying FPA concepts over...\n"); + pap_psMetadataCopy(toFPA->concepts, fromFPA->concepts); + pap_psMetadataCopy(toFPA->phu, fromFPA->phu); for (int i = 0; i < toChips->n; i++) { - pmChip *toChip = toChips->data[i]; - pmChip *fromChip = NULL; - // Select the correct chip - if (toChips->n == 1 && fromChips->n > chipNum) { - fromChip = fromChips->data[chipNum]; - } else if (fromChips->n == 1 && toChips->n > chipNum) { - fromChip = fromChips->data[0]; - } else if (toChips->n == fromChips->n) { - fromChip = fromChips->data[i]; - } else { - psError(PS_ERR_IO, true, "Unable to discern intended chip.\n"); - return false; - } - - psArray *toCells = toChip->cells; // Array of cells - psArray *fromCells = fromChip->cells; // Array of cells - - for (int j = 0; j < toCells->n; j++) { - pmCell *toCell = toCells->data[j]; - pmCell *fromCell = NULL; - // Select the correct cell - if (toCells->n == 1 && fromCells->n > cellNum) { - fromCell = fromCells->data[cellNum]; - } else if (fromCells->n == 1 && toCells->n > cellNum) { - fromCell = fromCells->data[0]; - } else if (toCells->n == fromCells->n) { - fromCell = fromCells->data[j]; - } else { - psError(PS_ERR_IO, true, "Unable to discern intended cell.\n"); - return false; - } - - psTrace(__func__, 5, "Putting chip %d, cell %d in the bucket.\n", i, j); - psListAdd(targetCells, PS_LIST_TAIL, toCell); - psListAdd(sourceCells, PS_LIST_TAIL, fromCell); - -#if 1 - int xParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.XPARITY"); - int yParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.YPARITY"); - int xParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.XPARITY"); - int yParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.YPARITY"); - psTrace(__func__, 3, "Chip %d, cell %d: in (%d,%d) out (%d,%d)\n", i, j, xParityIn, yParityIn, - xParityOut, yParityOut); + pmChip *toChip = toChips->data[i]; + pmChip *fromChip = NULL; + // Select the correct chip + if (toChips->n == 1 && fromChips->n > chipNum) { + fromChip = fromChips->data[chipNum]; + } else if (fromChips->n == 1 && toChips->n > chipNum) { + fromChip = fromChips->data[0]; + } else if (toChips->n == fromChips->n) { + fromChip = fromChips->data[i]; + } else { + psError(PS_ERR_IO, true, "Unable to discern intended chip.\n"); + return false; + } + + if (! fromChip->valid) { + toChip->valid = false; + continue; + } + + psTrace(__func__, 2, "Copying chip %d concepts over...\n", i); + pap_psMetadataCopy(toChip->concepts, fromChip->concepts); + + psArray *toCells = toChip->cells; // Array of cells + psArray *fromCells = fromChip->cells; // Array of cells + + for (int j = 0; j < toCells->n; j++) { + pmCell *toCell = toCells->data[j]; + pmCell *fromCell = NULL; + // Select the correct cell + if (toCells->n == 1 && fromCells->n > cellNum) { + fromCell = fromCells->data[cellNum]; + } else if (fromCells->n == 1 && toCells->n > cellNum) { + fromCell = fromCells->data[0]; + } else if (toCells->n == fromCells->n) { + fromCell = fromCells->data[j]; + } else { + psError(PS_ERR_IO, true, "Unable to discern intended cell.\n"); + return false; + } + + if (! fromCell->valid) { + toCell->valid = false; + continue; + } + +#ifdef TESTING + int xParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.XPARITY"); + int yParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.YPARITY"); + int xParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.XPARITY"); + int yParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.YPARITY"); + psTrace(__func__, 3, "Chip %d, cell %d: in (%d,%d) out (%d,%d)\n", i, j, xParityIn, yParityIn, + xParityOut, yParityOut); #endif - // Copy the cell contents over - toCell->readouts = psMemIncrRefCounter(fromCell->readouts); - - if (toCell->hdu && strlen(toCell->hdu->extname) > 0) { - // Splice the component cells - p_pmHDU *hdu = toCell->hdu; - if (! hdu->header) { - hdu->header = psMetadataAlloc(); - } - hdu->images = spliceCells(targetCells, sourceCells, positionDependent); - // Purge the lists - while (psListRemove(targetCells, PS_LIST_TAIL)); - while (psListRemove(sourceCells, PS_LIST_TAIL)); - } - - } - - if (toChip->hdu && strlen(toChip->hdu->extname) > 0) { - // Splice the component cells - p_pmHDU *hdu = toChip->hdu; - if (! hdu->header) { - hdu->header = psMetadataAlloc(); - } - hdu->images = spliceCells(targetCells, sourceCells, positionDependent); - // Purge the lists - while (psListRemove(targetCells, PS_LIST_TAIL)); - while (psListRemove(sourceCells, PS_LIST_TAIL)); - } +#ifdef TESTING + psMetadataPrint(toCell->concepts, 7); +#endif + + // Copy the concepts over + + // Need to be a little tricky here --- some concepts are already set from the camera configuration + // file (e.g., CELL.NAME), and we want to preserve these, not replace it with the old value. + psTrace(__func__, 3, "Copying cell %d concepts over...\n", j); + psMetadata *newConcepts = pap_psMetadataCopy(NULL, fromCell->concepts); + pap_psMetadataCopy(newConcepts, toCell->concepts); // This preserves any already existing concepts + psFree(toCell->concepts); + toCell->concepts = newConcepts; + + // Now, we need to check if the following concepts for the target cell are defined statically: + // CELL.XPARITY, CELL.YPARITY, CELL.XBIN, CELL.YBIN, CELL.X0, CELL.Y0 + // 1. If they are specified by the header, we can do what we want (because we can set the header). + // 2. If they are not specified by the header, we use those values. + const char *checkConcepts[4] = { "CELL.XPARITY", + "CELL.YPARITY", + "CELL.XBIN", + "CELL.YBIN", + "CELL.X0", + "CELL.Y0" }; // Concepts that we need to check + bool mdok = false; // Result of MD lookup + psMetadata *translation = psMetadataLookupMD(&mdok, toFPA->camera, "TRANSLATION"); // Header + // translation table + if (!mdok || ! translation) { + psError(PS_ERR_IO, true, "Unable to find TRANSLATION in camera configuration!\n"); + return false; + } + for (int c = 0; c < 4; c++) { + psString headerName = psMetadataLookupString(&mdok, translation, checkConcepts[c]); // Name of + // header holding the concept, or NULL + if (!mdok || strlen(headerName) <= 0) { + p_pmCellIngestConcept(toFPA, toChip, toCell, db, checkConcepts[c]); + } + } + +#ifdef TESTING + psMetadataPrint(toCell->concepts, 7); +#endif + + psTrace(__func__, 5, "Putting chip %d, cell %d in the bucket.\n", i, j); + psListAdd(targetCells, PS_LIST_TAIL, toCell); + psListAdd(sourceCells, PS_LIST_TAIL, fromCell); + + // Copy the cell contents over + toCell->readouts = psMemIncrRefCounter(fromCell->readouts); + + if (toCell->hdu && strlen(toCell->hdu->extname) > 0) { + // Splice the component cells + p_pmHDU *hdu = toCell->hdu; + if (! hdu->header) { + hdu->header = psMetadataAlloc(); + } + hdu->images = spliceCells(targetCells, sourceCells, positionDependent); + // Purge the lists + while (psListRemove(targetCells, PS_LIST_TAIL)); + while (psListRemove(sourceCells, PS_LIST_TAIL)); + } + + } + + if (toChip->hdu && strlen(toChip->hdu->extname) > 0) { + // Splice the component cells + p_pmHDU *hdu = toChip->hdu; + if (! hdu->header) { + hdu->header = psMetadataAlloc(); + } + hdu->images = spliceCells(targetCells, sourceCells, positionDependent); + // Purge the lists + while (psListRemove(targetCells, PS_LIST_TAIL)); + while (psListRemove(sourceCells, PS_LIST_TAIL)); + } } if (toFPA->hdu && strlen(toFPA->hdu->extname) > 0) { - // Splice the component cells - p_pmHDU *hdu = toFPA->hdu; - if (! hdu->header) { - hdu->header = psMetadataAlloc(); - } - hdu->images = spliceCells(targetCells, sourceCells, positionDependent); - // Purge the lists - while (psListRemove(targetCells, PS_LIST_TAIL)); - while (psListRemove(sourceCells, PS_LIST_TAIL)); + // Splice the component cells + p_pmHDU *hdu = toFPA->hdu; + if (! hdu->header) { + hdu->header = psMetadataAlloc(); + } + hdu->images = spliceCells(targetCells, sourceCells, positionDependent); + // Purge the lists + while (psListRemove(targetCells, PS_LIST_TAIL)); + while (psListRemove(sourceCells, PS_LIST_TAIL)); } Index: /trunk/archive/scripts/src/phase2/pmFPARead.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPARead.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPARead.c (revision 5786) @@ -7,5 +7,5 @@ #include "pmFPARead.h" -#include "papStuff.h" // For "split" +#include "papStuff.h" // For "split" // NOTE: Need to deal with header inheritance @@ -18,5 +18,5 @@ // 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 + psFits *fits // The FITS file from which to read ) { @@ -25,17 +25,17 @@ 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; - } + 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; + 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; + psError(PS_ERR_IO, false, "Unable to read FITS header!\n"); + return false; } @@ -44,22 +44,22 @@ 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); + 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"); + 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 + 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 + 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 @@ -67,29 +67,29 @@ 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); + 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); #ifndef PRODUCTION - // 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"); + // 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"); #endif - psFree(image); - } else { - pixels->data[i] = image; - } - psTrace(__func__, 10, "Done\n"); + 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; - } + psError(PS_ERR_IO, false, "Unable to read image for extension %s, plane %d\n", extName, i); + return false; + } } @@ -105,21 +105,21 @@ // 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 + p_pmHDU *hdu // Pixel data, containing image, mask, weights ) { - psArray *images = hdu->images; // Array of images (each of which is a readout) - psArray *masks = hdu->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++) { - psImage *image = images->data[i]; // The i-th plane - psImage *mask = NULL; // The mask - if (masks) { - mask = masks->data[i]; - } - - 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 - // pmReadoutAlloc to pmReadoutAdd? + psImage *image = images->data[i]; // The i-th plane + psImage *mask = NULL; // The mask + if (masks) { + mask = masks->data[i]; + } + + 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 + // pmReadoutAlloc to pmReadoutAdd? } @@ -131,85 +131,86 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -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 +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 + 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"); - psMetadata *phu = psFitsReadHeader(NULL, fits); // Primary header - } - fpa->phu = 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 { + fpa->phu = psMemIncrRefCounter(phu); + } // Read in.... psTrace(__func__, 1, "Working on 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; - } + 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; + } } pmFPAIngestConcepts(fpa, db); - psArray *chips = fpa->chips; // Array of chips + 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 "valid" - if (! chip->valid) { - 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; - } - } - pmChipIngestConcepts(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 "valid" - if (! cell->valid) { - 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; - } - } - pmCellIngestConcepts(cell, db); - - psTrace(__func__, 5, "Allocating readouts for chip %d cell %d...\n", - i, j); - generateReadouts(cell, hdu); - } + pmChip *chip = chips->data[i]; // The chip + + // Only read chips marked "valid" + if (! chip->valid) { + 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; + } + } + pmChipIngestConcepts(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 "valid" + if (! cell->valid) { + 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; + } + } + pmCellIngestConcepts(cell, db); + + psTrace(__func__, 5, "Allocating readouts for chip %d cell %d...\n", + i, j); + generateReadouts(cell, hdu); + } } @@ -219,5 +220,5 @@ // 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 + pmCell *cell // The cell for which to translate ) { @@ -230,106 +231,106 @@ 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 + 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. + // 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. + // 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. + // 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. + // 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. + // 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. + // 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. } @@ -338,144 +339,146 @@ // Read a mask into the FPA -bool pmFPAReadMask(pmFPA *fpa, // FPA to read into - psFits *source // Source FITS file (for the original data) +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 + 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; + 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; + 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; + 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 + 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 + 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 = psMemIncrRefCounter(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); + } + + 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); + 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 + psFree(maskSource); + } // Valid cells + } // Iterating over cells + } // Valid chips } // Iterating over chips @@ -486,144 +489,146 @@ // 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) +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 + 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; + 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; + 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; + 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 + 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 + 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 = psMemIncrRefCounter(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); + } + + 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); + 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 + psFree(weightSource); + } // Valid cells + } // Iterating over cells + } // Valid chips } // Iterating over chips Index: /trunk/archive/scripts/src/phase2/pmFPAWrite.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFPAWrite.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmFPAWrite.c (revision 5786) @@ -7,12 +7,12 @@ #include "pmFPARead.h" -static bool writeHDU(psFits *fits, // FITS file to which to write - p_pmHDU *hdu // Pixel data to write - ) -{ - bool status = true; // Status of write, to return +static bool writeHDU(psFits *fits, // FITS file to which to write + p_pmHDU *hdu // Pixel data to write + ) +{ + bool status = true; // Status of write, to return for (int i = 0; i < hdu->images->n; i++) { - status &= psFitsWriteImage(fits, hdu->header, hdu->images->data[i], i); - // XXX: Insert here the writing on mask and weight images + status &= psFitsWriteImage(fits, hdu->header, hdu->images->data[i], i); + // XXX: Insert here the writing on mask and weight images } @@ -21,46 +21,51 @@ -bool pmFPAWrite(psFits *fits, // FITS file to which to write - pmFPA *fpa, // FPA to write - psDB *db // Database to update - ) -{ - bool status = true; // Status of writing, to return - +bool pmFPAWrite(psFits *fits, // FITS file to which to write + pmFPA *fpa, // FPA to write + psDB *db // Database to update + ) +{ + bool status = true; // Status of writing, to return + + psTrace(__func__, 7, "Outgesting FPA concepts...\n"); pmFPAOutgestConcepts(fpa, db); // Write the primary header if (fpa->phu) { - status &= psFitsMoveExtNum(fits, 0, false); - status &= psFitsWriteHeader(fpa->phu, fits); - } - - psArray *chips = fpa->chips; // Array of component chips + status &= psFitsMoveExtNum(fits, 0, false); + status &= psFitsWriteHeader(fpa->phu, fits); + } + + psArray *chips = fpa->chips; // Array of component chips for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The component chip - if (chip->valid) { - pmChipOutgestConcepts(chip, db); - - psArray *cells = chip->cells; // Array of component cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // The component cell - if (cell->valid) { - pmCellOutgestConcepts(cell, db); - - if (cell->hdu && strlen(cell->hdu->extname) > 0) { - status &= writeHDU(fits, cell->hdu); - } - } - } - - if (chip->hdu && strlen(chip->hdu->extname) > 0) { - status &= writeHDU(fits, chip->hdu); - } - } + pmChip *chip = chips->data[i]; // The component chip + if (chip->valid) { + psTrace(__func__, 1, "Writing out chip %d...\n", i); + + pmChipOutgestConcepts(chip, db); + + psArray *cells = chip->cells; // Array of component cells + for (int j = 0; j < cells->n; j++) { + pmCell *cell = cells->data[j]; // The component cell + if (cell->valid) { + psTrace(__func__, 2, "Writing out cell, %d...\n", j); + + pmCellOutgestConcepts(cell, db); + + if (cell->hdu && strlen(cell->hdu->extname) > 0) { + status &= writeHDU(fits, cell->hdu); + } + } + } + + if (chip->hdu && strlen(chip->hdu->extname) > 0) { + status &= writeHDU(fits, chip->hdu); + } + } } if (fpa->hdu && strlen(fpa->hdu->extname) > 0) { - status &= writeHDU(fits, fpa->hdu); + status &= writeHDU(fits, fpa->hdu); } @@ -71,101 +76,101 @@ -bool pmFPAWriteMask(pmFPA *fpa, // FPA containing mask to write - psFits *fits // FITS file for image +bool pmFPAWriteMask(pmFPA *fpa, // FPA containing mask to write + psFits *fits // FITS file for image ) { const psMetadata *camera = fpa->camera; // Camera configuration for FPA - bool mdok = false; // Status of MD lookup + 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; + 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; + 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; + 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 + 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 write the pixels - psFits *maskDest = psMemIncrRefCounter(fits); // Destination of mask image - psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name - 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); - } - - psFree(maskDest); - maskDest = psFitsAlloc(filename); - if (extname) { - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - } - 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); - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - psFree(extname); - } - - // We've arrived where the pixels are. Now we need to write them out. - psArray *readouts = cell->readouts; // The array of readouts - for (int readNum = 0; readNum < readouts->n; readNum++) { - pmReadout *readout = readouts->data[readNum]; // The readout of interest - if (! readout->mask) { - psLogMsg(__func__, PS_LOG_WARN, "No mask to write out in %d,%d,%d\n", - chipNum, cellNum, readNum); - } else { - // XXX: Need to add the extname to the existing header - if (! psFitsWriteImage(maskDest, header, readout->mask, readNum)) { - psError(PS_ERR_IO, false, "Unable to write mask plane %d in extension %s\n", - readNum, hdu->extname); - return false; - } - } - } // Iterating over readouts - psFree(header); - psFree(maskDest); - } // Valid cells - } // Iterating over cells - } // Valid chips + 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 write the pixels + psFits *maskDest = psMemIncrRefCounter(fits); // Destination of mask image + psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name + 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); + } + + psFree(maskDest); + maskDest = psFitsOpen(filename, "rw"); + if (extname) { + psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); + } + 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); + psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); + psFree(extname); + } + + // We've arrived where the pixels are. Now we need to write them out. + psArray *readouts = cell->readouts; // The array of readouts + for (int readNum = 0; readNum < readouts->n; readNum++) { + pmReadout *readout = readouts->data[readNum]; // The readout of interest + if (! readout->mask) { + psLogMsg(__func__, PS_LOG_WARN, "No mask to write out in %d,%d,%d\n", + chipNum, cellNum, readNum); + } else { + // XXX: Need to add the extname to the existing header + if (! psFitsWriteImage(maskDest, header, readout->mask, readNum)) { + psError(PS_ERR_IO, false, "Unable to write mask plane %d in extension %s\n", + readNum, hdu->extname); + return false; + } + } + } // Iterating over readouts + psFree(header); + psFree(maskDest); + } // Valid cells + } // Iterating over cells + } // Valid chips } // Iterating over chips @@ -174,100 +179,100 @@ -bool pmFPAWriteWeight(pmFPA *fpa, // FPA containing mask to write - psFits *fits // FITS file for image +bool pmFPAWriteWeight(pmFPA *fpa, // FPA containing mask to write + psFits *fits // FITS file for image ) { const psMetadata *camera = fpa->camera; // Camera configuration for FPA - bool mdok = false; // Status of MD lookup + 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; + 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; + 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; + 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 + 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 write the pixels - psFits *weightDest = psMemIncrRefCounter(fits); // Destination of weight image - psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name - 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); - } - - psFree(weightDest); - weightDest = psFitsAlloc(filename); - if (extname) { - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - } - 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); - psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); - psFree(extname); - } - - // We've arrived where the pixels are. Now we need to write them out. - psArray *readouts = cell->readouts; // The array of readouts - for (int readNum = 0; readNum < readouts->n; readNum++) { - pmReadout *readout = readouts->data[readNum]; // The readout of interest - if (! readout->weight) { - psLogMsg(__func__, PS_LOG_WARN, "No weight image to write out in %d,%d,%d\n", - chipNum, cellNum, readNum); - } else { - if (! psFitsWriteImage(weightDest, header, readout->weight, readNum)) { - psError(PS_ERR_IO, false, "Unable to write weight plane %d in extension %s\n", - readNum, hdu->extname); - return false; - } - } - } // Iterating over readouts - psFree(header); - psFree(weightDest); - } // Valid cells - } // Iterating over cells - } // Valid chips + 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 write the pixels + psFits *weightDest = psMemIncrRefCounter(fits); // Destination of weight image + psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name + 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); + } + + psFree(weightDest); + weightDest = psFitsOpen(filename, "rw"); + if (extname) { + psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); + } + 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); + psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); + psFree(extname); + } + + // We've arrived where the pixels are. Now we need to write them out. + psArray *readouts = cell->readouts; // The array of readouts + for (int readNum = 0; readNum < readouts->n; readNum++) { + pmReadout *readout = readouts->data[readNum]; // The readout of interest + if (! readout->weight) { + psLogMsg(__func__, PS_LOG_WARN, "No weight image to write out in %d,%d,%d\n", + chipNum, cellNum, readNum); + } else { + if (! psFitsWriteImage(weightDest, header, readout->weight, readNum)) { + psError(PS_ERR_IO, false, "Unable to write weight plane %d in extension %s\n", + readNum, hdu->extname); + return false; + } + } + } // Iterating over readouts + psFree(header); + psFree(weightDest); + } // Valid cells + } // Iterating over cells + } // Valid chips } // Iterating over chips Index: /trunk/archive/scripts/src/phase2/pmSubtractBias.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmSubtractBias.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/pmSubtractBias.c (revision 5786) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-03 01:30:32 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-14 03:47:35 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -214,5 +214,4 @@ psF32 x; psS32 i; - printf("Got here\n"); if (fit == PM_FIT_POLYNOMIAL) { // Fit a polynomial to the old overscan vector. Index: /trunk/archive/scripts/src/phase2/psAdditionals.c =================================================================== --- /trunk/archive/scripts/src/phase2/psAdditionals.c (revision 5785) +++ /trunk/archive/scripts/src/phase2/psAdditionals.c (revision 5786) @@ -5,6 +5,6 @@ -psMetadata *psMetadataCopy(psMetadata *out, - const psMetadata *in) +psMetadata *pap_psMetadataCopy(psMetadata *out, + const psMetadata *in) { PS_ASSERT_PTR_NON_NULL(in,NULL); @@ -20,13 +20,19 @@ // Need to look for MULTI, which won't be picked up using the iterator. psMetadataItem *multiCheckItem = psMetadataLookup(in, inItem->name); - int multiFlag = 0; // Flag to indicate MULTI or not + unsigned int flag = PS_META_REPLACE; // Flag to indicate MULTI; otherwise, replace if (multiCheckItem->type == PS_DATA_METADATA_MULTI) { - multiFlag = PS_DATA_METADATA_MULTI; - } + psTrace(__func__, 10, "MULTI: %s (%s)\n", inItem->name, inItem->comment); + flag = PS_DATA_METADATA_MULTI; + } + + psTrace(__func__, 5, "Copying %s (%s)...\n", inItem->name, inItem->comment); #define PS_METADATA_COPY_CASE(NAME,TYPE) \ case PS_TYPE_##NAME: \ - psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_TYPE_##NAME | multiFlag, inItem->comment, \ - inItem->data.TYPE); \ + if (! psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_TYPE_##NAME | flag, inItem->comment, \ + inItem->data.TYPE)) { \ + psErrorStackPrint(stderr, "Error copying %s (%s) in the metadata\n", inItem->name, \ + inItem->comment); \ + } \ break; @@ -34,16 +40,16 @@ // Numerical types PS_METADATA_COPY_CASE(BOOL,B); - PS_METADATA_COPY_CASE(S8, S8); - PS_METADATA_COPY_CASE(S16, S16); - PS_METADATA_COPY_CASE(S32, S32); - PS_METADATA_COPY_CASE(U8, U8); - PS_METADATA_COPY_CASE(U16, U16); - PS_METADATA_COPY_CASE(U32, U32); - PS_METADATA_COPY_CASE(F32, F32); - PS_METADATA_COPY_CASE(F64, F64); + PS_METADATA_COPY_CASE(S8,S8); + PS_METADATA_COPY_CASE(S16,S16); + PS_METADATA_COPY_CASE(S32,S32); + PS_METADATA_COPY_CASE(U8,U8); + PS_METADATA_COPY_CASE(U16,U16); + PS_METADATA_COPY_CASE(U32,U32); + PS_METADATA_COPY_CASE(F32,F32); + PS_METADATA_COPY_CASE(F64,F64); // String: relying on the fact that this will copy the string, not point at it. case PS_DATA_STRING: - psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_STRING | multiFlag, inItem->comment, + psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_STRING | flag, inItem->comment, inItem->data.V); break; @@ -52,6 +58,6 @@ case PS_DATA_METADATA: { - psMetadata *metadata = psMetadataCopy(NULL, inItem->data.md); - psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_METADATA | multiFlag, inItem->comment, + psMetadata *metadata = pap_psMetadataCopy(NULL, inItem->data.md); + psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_METADATA | flag, inItem->comment, metadata); break; @@ -60,5 +66,7 @@ default: numPointers++; - psMetadataItemAlloc(inItem->name, inItem->type, inItem->comment, inItem->data.V); + psTrace(__func__, 10, "Copying a pointer in the metadata: %x\n", inItem->type); + psMetadataAdd(out, PS_LIST_TAIL, inItem->name, inItem->type | flag, inItem->comment, + inItem->data.V); break; } @@ -68,5 +76,5 @@ if (numPointers > 0) { psLogMsg(__func__, PS_LOG_WARN, "Forced to copy %d pointers when copying metadata. Updating the " - "copied psMetadata will affect the original!\n"); + "copied psMetadata will affect the original!\n", numPointers); } Index: /trunk/archive/scripts/src/phase2/psAdditionals.h =================================================================== --- /trunk/archive/scripts/src/phase2/psAdditionals.h (revision 5785) +++ /trunk/archive/scripts/src/phase2/psAdditionals.h (revision 5786) @@ -5,5 +5,5 @@ // Deep copy of metadata -psMetadata *psMetadataCopy(psMetadata *out, const psMetadata *in); +psMetadata *pap_psMetadataCopy(psMetadata *out, const psMetadata *in); // Get a value from the metadata that we believe should be metadata.