Changeset 5786
- Timestamp:
- Dec 13, 2005, 5:47:35 PM (20 years ago)
- Location:
- trunk/archive/scripts/src/phase2
- Files:
-
- 20 edited
-
Makefile.am (modified) (2 diffs)
-
ipprc.config (modified) (1 diff)
-
megacam_raw.config (modified) (1 diff)
-
nonlin.dat (modified) (1 diff)
-
papPhase2.c (modified) (14 diffs)
-
phase2.config (modified) (1 diff)
-
pmChipMosaic.c (modified) (3 diffs)
-
pmChipMosaic.h (modified) (1 diff)
-
pmFPA.c (modified) (2 diffs)
-
pmFPA.h (modified) (2 diffs)
-
pmFPAConceptsGet.c (modified) (18 diffs)
-
pmFPAConceptsGet.h (modified) (1 diff)
-
pmFPAConceptsSet.c (modified) (5 diffs)
-
pmFPAConstruct.c (modified) (1 diff)
-
pmFPAMorph.c (modified) (13 diffs)
-
pmFPARead.c (modified) (11 diffs)
-
pmFPAWrite.c (modified) (4 diffs)
-
pmSubtractBias.c (modified) (2 diffs)
-
psAdditionals.c (modified) (6 diffs)
-
psAdditionals.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/archive/scripts/src/phase2/Makefile.am
r5651 r5786 12 12 pmFPAConceptsSet.c \ 13 13 pmFPAConstruct.c \ 14 pmFPAMorph.c \15 14 pmFPARead.c \ 16 15 pmFPAWrite.c \ … … 29 28 pmFPAConceptsSet.h \ 30 29 pmFPAConstruct.h \ 31 pmFPAMorph.h \32 30 pmFPARead.h \ 33 31 pmFPAWrite.h \ -
trunk/archive/scripts/src/phase2/ipprc.config
r5651 r5786 23 23 pap S32 10 24 24 pmFPAPrint S32 10 25 pmFPAWrite S32 10 25 26 # pmConfigRead S32 10 26 27 # pmFPAWriteMask S32 10 27 28 # pmFPAWriteWeight S32 10 28 29 pmChipMosaic S32 10 30 pmFPAMorph S32 10 31 spliceCells S32 10 32 spliceImage S32 10 33 setConceptItem S32 10 34 # pap_psMetadataCopy S32 9 29 35 END -
trunk/archive/scripts/src/phase2/megacam_raw.config
r5621 r5786 133 133 # Default PS concepts that may be specified by value 134 134 DEFAULTS METADATA 135 CELL.READDIR S32 2# Cell is read in x direction135 CELL.READDIR S32 1 # Cell is read in x direction 136 136 CELL.BAD S32 0 137 137 CELL.TIMESYS STR UTC -
trunk/archive/scripts/src/phase2/nonlin.dat
r5621 r5786 103 103 9900 99 104 104 10000 100 105 10001 -1106 1e6 -1105 10001 100 106 1e6 100 -
trunk/archive/scripts/src/phase2/papPhase2.c
r5651 r5786 17 17 #include "pmSubtractBias.h" 18 18 #include "pmChipMosaic.h" 19 #include "pmFPAMorph.h"19 //#include "pmFPAMorph.h" 20 20 21 21 // Phase 2 needs to: … … 39 39 40 40 // phase2 INPUT.fits OUTPUT.fits -bias BIAS.fits -dark DARK.fits -flat FLAT.fits -chip CHIP 41 // 41 // 42 42 // Most are self-explanatory. "-chip" says "only work on this particular chip". 43 43 44 44 45 #define RECIPE "PHASE2" // Name of the recipe to use45 #define RECIPE "PHASE2" // Name of the recipe to use 46 46 47 47 static psMemId memPrintAlloc(const psMemBlock *mb) 48 48 { 49 49 printf("Allocated memory block %ld (%ld):\n" 50 "\tFile %s, line %d, size %d\n"51 "\tPosts: %x %x %x\n",52 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,53 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));50 "\tFile %s, line %d, size %d\n" 51 "\tPosts: %x %x %x\n", 52 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, 53 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); 54 54 return 0; 55 55 } … … 57 57 { 58 58 printf("Freed memory block %ld (%ld):\n" 59 "\tFile %s, line %d, size %d\n"60 "\tPosts: %x %x %x\n",61 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,62 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));59 "\tFile %s, line %d, size %d\n" 60 "\tPosts: %x %x %x\n", 61 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, 62 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); 63 63 return 0; 64 64 } … … 67 67 psMemBlock *mb = ((psMemBlock*)ptr) - 1; 68 68 printf("Memory block %ld (%ld):\n" 69 "\tFile %s, line %d, size %d\n"70 "\tPosts: %x %x %x\n",71 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,72 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));69 "\tFile %s, line %d, size %d\n" 70 "\tPosts: %x %x %x\n", 71 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, 72 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); 73 73 } 74 74 … … 102 102 psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", 0, "Chip number to process (if positive)", -1); 103 103 if (! psArgumentParse(arguments, &argc, argv) || argc != 3) { 104 printf("\nPan-STARRS Phase 2 processing\n\n");105 printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]);106 psArgumentHelp(arguments);107 psFree(arguments);108 exit(EXIT_FAILURE);109 } 110 const char *inputName = argv[1]; // Name of input image111 const char *outputName = argv[2]; // Name of output image104 printf("\nPan-STARRS Phase 2 processing\n\n"); 105 printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]); 106 psArgumentHelp(arguments); 107 psFree(arguments); 108 exit(EXIT_FAILURE); 109 } 110 const char *inputName = argv[1]; // Name of input image 111 const char *outputName = argv[2]; // Name of output image 112 112 const char *biasName = psMetadataLookupString(NULL, arguments, "-bias"); // Name of bias image 113 113 const char *darkName = psMetadataLookupString(NULL, arguments, "-dark"); // Name of dark image … … 123 123 124 124 // Open the input 125 #ifdef PRODUCTION126 125 psFits *inputFile = psFitsOpen(inputName, "r"); // File handle for FITS file 127 #else128 psFits *inputFile = psFitsAlloc(inputName); // File handle for FITS file129 #endif130 126 if (! inputFile) { 131 psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName);127 psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName); 132 128 exit(EXIT_FAILURE); 133 129 } 134 130 psMetadata *header = psFitsReadHeader(NULL, inputFile); // FITS header 135 131 #if PRODUCTION 136 psDB *database = pmConfigDB(site); // Database handle132 psDB *database = pmConfigDB(site); // Database handle 137 133 #else 138 psDB *database = NULL; // Database handle134 psDB *database = NULL; // Database handle 139 135 #endif 140 136 141 137 // Open the output and output mask 142 138 // We do it here so that we don't process the whole lot and then find out we can't open the output file 143 #ifdef PRODUCTION144 139 psFits *outputFile = psFitsOpen(outputName, "w"); 145 #else146 psFits *outputFile = psFitsAlloc(outputName);147 #endif148 140 if (! outputFile) { 149 psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName);150 exit(EXIT_FAILURE);141 psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName); 142 exit(EXIT_FAILURE); 151 143 } 152 144 153 145 // Get camera configuration from header if not already defined 154 146 if (! camera) { 155 camera = pmConfigCameraFromHeader(site, header);156 if (! camera) {157 psErrorStackPrint(stderr, "Can't find camera configuration!\n");158 exit(EXIT_FAILURE);159 }147 camera = pmConfigCameraFromHeader(site, header); 148 if (! camera) { 149 psErrorStackPrint(stderr, "Can't find camera configuration!\n"); 150 exit(EXIT_FAILURE); 151 } 160 152 } else if (! pmConfigValidateCamera(camera, header)) { 161 psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName);162 exit(EXIT_FAILURE);153 psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName); 154 exit(EXIT_FAILURE); 163 155 } 164 156 if (! recipe && !(recipe = pmConfigRecipeFromCamera(camera, RECIPE))) { … … 175 167 176 168 // Set various tasks 177 bool doMask = false; // Mask bad pixles178 bool doNonLin = false; // Non-linearity correction179 bool doBias = false; // Bias subtraction180 bool doDark = false; // Dark subtraction181 bool doOverscan = false; // Overscan subtraction182 bool doFlat = false; // Flat-field normalisation183 bool doFringe = false; // Fringe subtraction184 bool doSource = false; // Source identification and photometry185 bool doAstrom = false; // Astrometry169 bool doMask = false; // Mask bad pixles 170 bool doNonLin = false; // Non-linearity correction 171 bool doBias = false; // Bias subtraction 172 bool doDark = false; // Dark subtraction 173 bool doOverscan = false; // Overscan subtraction 174 bool doFlat = false; // Flat-field normalisation 175 bool doFringe = false; // Fringe subtraction 176 bool doSource = false; // Source identification and photometry 177 bool doAstrom = false; // Astrometry 186 178 pmOverscanAxis overscanMode = PM_OVERSCAN_NONE; // Axis for overscan 187 179 pmFit overscanFitType = PM_FIT_NONE; // Fit type for overscan 188 int overscanBins = 1; // Number of pixels per bin for overscan189 psStats *overscanStats = NULL; // Statistics for overscan190 void *overscanFit = NULL; // Overscan fit (polynomial or spline)180 int overscanBins = 1; // Number of pixels per bin for overscan 181 psStats *overscanStats = NULL; // Statistics for overscan 182 void *overscanFit = NULL; // Overscan fit (polynomial or spline) 191 183 192 184 if (psMetadataLookupBool(NULL, recipe, "MASK")) { 193 if (strlen(maskName) > 0) {194 doMask = true;195 } else {196 psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking "197 "will be performed.\n");198 }185 if (strlen(maskName) > 0) { 186 doMask = true; 187 } else { 188 psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking " 189 "will be performed.\n"); 190 } 199 191 } 200 192 if (psMetadataLookupBool(NULL, recipe, "NONLIN")) { 201 doNonLin = true;193 doNonLin = true; 202 194 } 203 195 if (psMetadataLookupBool(NULL, recipe, "BIAS")) { 204 if (strlen(biasName) > 0) {205 doBias = true;206 } else {207 psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- "208 "no bias subtraction will be performed.\n");209 }196 if (strlen(biasName) > 0) { 197 doBias = true; 198 } else { 199 psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- " 200 "no bias subtraction will be performed.\n"); 201 } 210 202 } 211 203 if (psMetadataLookupBool(NULL, recipe, "DARK")) { 212 if (strlen(darkName) > 0) {213 doDark = true;214 } else {215 psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- "216 "no dark subtraction will be performed.\n");217 }204 if (strlen(darkName) > 0) { 205 doDark = true; 206 } else { 207 psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- " 208 "no dark subtraction will be performed.\n"); 209 } 218 210 } 219 211 if (psMetadataLookupBool(NULL, recipe, "OVERSCAN")) { 220 doOverscan = true;221 psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE");222 if (strcasecmp(mode, "INDIVIDUAL") == 0) {223 overscanMode = PM_OVERSCAN_ROWS;224 // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip225 } else if (strcasecmp(mode, "ALL") == 0) {226 overscanMode = PM_OVERSCAN_ALL;227 } else if (strcasecmp(mode, "NONE") != 0) {228 psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:"229 " assuming NONE.\n", mode);230 }231 psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT");232 if (strcasecmp(fit, "POLYNOMIAL") == 0) {233 overscanFitType = PM_FIT_POLYNOMIAL;234 int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit235 overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);236 } else if (strcasecmp(fit, "SPLINE") == 0) {237 overscanFitType = PM_FIT_SPLINE;238 int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit239 overscanFit = NULL;240 // overscanFit = psSpline1DAlloc();241 } else if (strcasecmp(fit, "NONE") != 0) {242 psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:"243 " assuming NONE.\n", fit);244 }245 overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN");246 if (overscanBins <= 0) {247 psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins);248 overscanBins = 1;249 }250 psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT");251 if (strcasecmp(stat, "MEAN")) {252 overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);253 } else if (strcasecmp(stat, "MEDIAN")) {254 overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);255 } else {256 psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat);257 overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 258 }212 doOverscan = true; 213 psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE"); 214 if (strcasecmp(mode, "INDIVIDUAL") == 0) { 215 overscanMode = PM_OVERSCAN_ROWS; 216 // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip 217 } else if (strcasecmp(mode, "ALL") == 0) { 218 overscanMode = PM_OVERSCAN_ALL; 219 } else if (strcasecmp(mode, "NONE") != 0) { 220 psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:" 221 " assuming NONE.\n", mode); 222 } 223 psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT"); 224 if (strcasecmp(fit, "POLYNOMIAL") == 0) { 225 overscanFitType = PM_FIT_POLYNOMIAL; 226 int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit 227 overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); 228 } else if (strcasecmp(fit, "SPLINE") == 0) { 229 overscanFitType = PM_FIT_SPLINE; 230 int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit 231 overscanFit = NULL; 232 // overscanFit = psSpline1DAlloc(); 233 } else if (strcasecmp(fit, "NONE") != 0) { 234 psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:" 235 " assuming NONE.\n", fit); 236 } 237 overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN"); 238 if (overscanBins <= 0) { 239 psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins); 240 overscanBins = 1; 241 } 242 psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT"); 243 if (strcasecmp(stat, "MEAN")) { 244 overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 245 } else if (strcasecmp(stat, "MEDIAN")) { 246 overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); 247 } else { 248 psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat); 249 overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 250 } 259 251 } 260 252 261 253 if (psMetadataLookupBool(NULL, recipe, "FLAT")) { 262 if (strlen(flatName) > 0) {263 doFlat = true;264 } else {265 psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- "266 "no flat-fielding will be performed.\n");267 }254 if (strlen(flatName) > 0) { 255 doFlat = true; 256 } else { 257 psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- " 258 "no flat-fielding will be performed.\n"); 259 } 268 260 } 269 261 270 262 // Chip selection 271 263 if (chipNum >= 0) { 272 if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) ||273 ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) ||274 ! pmFPASelectChip(mask, chipNum)) {275 psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum);276 exit(EXIT_FAILURE);277 }278 psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum);264 if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) || 265 ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) || 266 ! pmFPASelectChip(mask, chipNum)) { 267 psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum); 268 exit(EXIT_FAILURE); 269 } 270 psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum); 279 271 } 280 272 … … 284 276 psLogMsg("phase2", PS_LOG_INFO, "Opening input image: %s\n", inputName); 285 277 if (! pmFPARead(input, inputFile, header, database)) { 286 psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName); 287 exit(EXIT_FAILURE); 288 } 289 290 pmFPAPrint(input); 291 exit(1); 292 293 #ifdef PRODUCTION 278 psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName); 279 exit(EXIT_FAILURE); 280 } 294 281 psFitsClose(inputFile); 295 #else296 psFree(inputFile);297 #endif298 282 299 283 #if 0 … … 302 286 //#else 303 287 { 304 // Generate mask and weight frame305 p_pmHDU *hdu = NULL;306 if (input->hdu) {307 hdu = input->hdu;308 }309 psArray *chips = input->chips;310 for (int chipNum = 0; chipNum < chips->n; chipNum++) {311 pmChip *chip = chips->data[chipNum];312 if (chip->valid) {313 if (chip->hdu) {314 hdu = chip->hdu;315 }316 psArray *cells = chip->cells;317 for (int cellNum = 0; cellNum < cells->n; cellNum++) {318 pmCell *cell = cells->data[cellNum];319 if (cell->valid) {320 if (cell->hdu) {321 hdu = cell->hdu;322 }323 324 hdu->masks = psArrayAlloc(hdu->images->n);325 hdu->weights = psArrayAlloc(hdu->images->n);326 psArray *readouts = cell->readouts;327 for (int readNum = 0; readNum < hdu->images->n; readNum++) {328 psImage *image = hdu->images->data[readNum];329 psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);330 psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);331 pmReadout *readout = readouts->data[readNum];332 for (int j = 0; j < image->numRows; j++) {333 for (int i = 0; i < image->numCols; i++) {334 mask->data.U8[j][i] = j + i;335 weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]);336 }337 }338 hdu->masks->data[readNum] = mask;339 hdu->weights->data[readNum] = weight;340 readout->mask = psMemIncrRefCounter(mask);341 readout->weight = psMemIncrRefCounter(weight);342 }343 }344 }345 }346 }288 // Generate mask and weight frame 289 p_pmHDU *hdu = NULL; 290 if (input->hdu) { 291 hdu = input->hdu; 292 } 293 psArray *chips = input->chips; 294 for (int chipNum = 0; chipNum < chips->n; chipNum++) { 295 pmChip *chip = chips->data[chipNum]; 296 if (chip->valid) { 297 if (chip->hdu) { 298 hdu = chip->hdu; 299 } 300 psArray *cells = chip->cells; 301 for (int cellNum = 0; cellNum < cells->n; cellNum++) { 302 pmCell *cell = cells->data[cellNum]; 303 if (cell->valid) { 304 if (cell->hdu) { 305 hdu = cell->hdu; 306 } 307 308 hdu->masks = psArrayAlloc(hdu->images->n); 309 hdu->weights = psArrayAlloc(hdu->images->n); 310 psArray *readouts = cell->readouts; 311 for (int readNum = 0; readNum < hdu->images->n; readNum++) { 312 psImage *image = hdu->images->data[readNum]; 313 psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); 314 psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); 315 pmReadout *readout = readouts->data[readNum]; 316 for (int j = 0; j < image->numRows; j++) { 317 for (int i = 0; i < image->numCols; i++) { 318 mask->data.U8[j][i] = j + i; 319 weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]); 320 } 321 } 322 hdu->masks->data[readNum] = mask; 323 hdu->weights->data[readNum] = weight; 324 readout->mask = psMemIncrRefCounter(mask); 325 readout->weight = psMemIncrRefCounter(weight); 326 } 327 } 328 } 329 } 330 } 347 331 } 348 332 #endif … … 354 338 // Load the calibration frames, if required 355 339 if (doBias) { 340 psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias 341 psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias 342 if (! pmConfigValidateCamera(camera, biasHeader)) { 343 psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName); 344 exit(EXIT_FAILURE); 345 } 346 psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName); 347 if (! pmFPARead(bias, biasFile, biasHeader, database)) { 348 psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName); 349 exit(EXIT_FAILURE); 350 } 351 psFree(biasHeader); 352 psFitsClose(biasFile); 353 } 354 355 if (doDark) { 356 psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark 357 psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark 358 if (! pmConfigValidateCamera(camera, darkHeader)) { 359 psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName); 360 exit(EXIT_FAILURE); 361 } 362 psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName); 363 if (! pmFPARead(dark, darkFile, darkHeader, database)) { 364 psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName); 365 exit(EXIT_FAILURE); 366 } 367 psFree(darkHeader); 368 psFitsClose(darkFile); 369 } 370 371 if (doMask) { 372 psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask 373 psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask 374 if (! pmConfigValidateCamera(camera, maskHeader)) { 375 psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName); 376 exit(EXIT_FAILURE); 377 } 378 psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName); 379 if (! pmFPARead(mask, maskFile, maskHeader, database)) { 380 psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName); 381 exit(EXIT_FAILURE); 382 } 383 psFree(maskHeader); 384 psFitsClose(maskFile); 385 } 386 387 if (doFlat) { 388 psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat 389 if (! flatFile) { 390 psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName); 391 exit(EXIT_FAILURE); 392 } 393 psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat 394 if (! pmConfigValidateCamera(camera, flatHeader)) { 395 psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName); 396 exit(EXIT_FAILURE); 397 } 398 psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName); 399 if (! pmFPARead(flat, flatFile, flatHeader, database)) { 400 psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName); 401 exit(EXIT_FAILURE); 402 } 403 psFree(flatHeader); 404 psFitsClose(flatFile); 405 } 406 407 psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2")); 408 409 psArray *inputChips = input->chips; // Array of chips in input image 410 psArray *biasChips = bias->chips; // Array of chips in bias image 411 psArray *darkChips = dark->chips; // Array of chips in dark image 412 psArray *maskChips = mask->chips; // Array of chips in mask image 413 psArray *flatChips = flat->chips; // Array of chips in flat image 414 for (int i = 0; i < inputChips->n; i++) { 415 pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image 416 pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image 417 pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image 418 pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image 419 pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image 420 421 if (! inputChip->valid) { 422 continue; 423 } 424 425 psArray *inputCells = inputChip->cells; // Array of cells in input image 426 psArray *biasCells = biasChip->cells; // Array of cells in bias image 427 psArray *darkCells = darkChip->cells; // Array of cells in dark image 428 psArray *maskCells = maskChip->cells; // Array of cells in mask image 429 psArray *flatCells = flatChip->cells; // Array of cells in flat image 430 431 for (int j = 0; j < inputCells->n; j++) { 432 pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image 433 pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image 434 pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag 435 pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image 436 pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image 437 438 if (! inputCell->valid) { 439 continue; 440 } 441 442 psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image 443 if (doBias && biasCell->readouts->n > 1) { 444 psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" 445 " used."); 446 } 447 if (doDark && darkCell->readouts->n > 1) { 448 psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" 449 " used."); 450 } 451 if (doMask && maskCell->readouts->n > 1) { 452 psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" 453 " used."); 454 } 455 if (doFlat && flatCell->readouts->n > 1) { 456 psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be" 457 " used."); 458 } 459 460 pmReadout *biasReadout = NULL; // Bias readout 461 pmReadout *maskReadout = NULL; // Mask readout 462 pmReadout *flatReadout = NULL; // Flat readout 463 psImage *biasImage = NULL; // Bias pixels 464 psImage *darkImage = NULL; // Dark pixels 465 float darkTime = 1.0; // Dark time for dark image 466 psImage *maskImage = NULL; // Mask pixels 467 468 if (doBias) { 469 biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image 470 biasImage = biasReadout->image; 471 } 472 if (doDark) { 473 pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image 474 darkImage = darkReadout->image; 475 darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME"); 476 if (darkTime <= 0.0) { 477 psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime); 478 exit(EXIT_FAILURE); 479 } 480 } 481 if (doMask) { 482 maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image 483 maskImage = maskReadout->image; 484 } 485 if (doFlat) { 486 flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image 487 } 488 489 for (int k = 0; k < inputReadouts->n; k++) { 490 pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image 491 psImage *inputImage = inputReadout->image; // The actual image 492 493 // Mask bad pixels 494 if (doMask) { 495 float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION"); 496 497 // Need to change this later to grow the mask by the size of the convolution kernel 498 (void)pmMaskBadPixels(inputReadout, maskImage, 499 PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation, 500 PM_MASK_TRAP, 0); 501 } 502 503 // Non-linearity correction 504 if (doNonLin) { 505 psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA"); 506 if (! dataItem) { 507 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to " 508 "find NONLIN.DATA in recipe --- ignored.\n"); 509 } else { 510 switch (dataItem->type) { 511 case PS_DATA_VECTOR: 512 { 513 // These are the polynomial coefficients 514 psVector *coeff = dataItem->data.V; // The coefficient vector 515 if (coeff->type.type != PS_TYPE_F64) { 516 psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version 517 coeff = temp; 518 } 519 psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, 520 PS_POLYNOMIAL_ORD); 521 psFree(correction->coeff); 522 correction->coeff = psMemIncrRefCounter(coeff->data.F64); 523 (void)pmNonLinearityPolynomial(inputReadout, correction); 524 psFree(coeff); 525 psFree(correction); 526 } 527 break; 528 case PS_DATA_STRING: 529 { 530 // This is a filename 531 psString name = dataItem->data.V; // Filename 532 psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0); 533 if (psLookupTableRead(table) <= 0) { 534 psErrorStackPrint(stderr, "Unable to read non-linearity correction file " 535 "%s --- ignored\n", name); 536 } else { 356 537 #ifdef PRODUCTION 357 psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias 538 (void)pmNonLinearityLookup(inputReadout, table); 358 539 #else 359 psFits *biasFile = psFitsAlloc(biasName); 360 #endif 361 psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias 362 if (! pmConfigValidateCamera(camera, biasHeader)) { 363 psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName); 364 exit(EXIT_FAILURE); 365 } 366 psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName); 367 if (! pmFPARead(bias, biasFile, biasHeader, database)) { 368 psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName); 369 exit(EXIT_FAILURE); 370 } 371 psFree(biasHeader); 540 psVector *influx = table->values->data[0]; 541 psVector *outflux = table->values->data[1]; 542 (void)pmNonLinearityLookup(inputReadout, table->values->data[0], 543 table->values->data[1]); 544 #endif 545 } 546 psFree(table); 547 } 548 break; 549 case PS_DATA_METADATA: 550 { 551 // This is a menu 552 psMetadata *options = dataItem->data.V; // Options with concept values as keys 553 bool mdok = false; // Success of MD lookup 554 psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE"); 555 if (! mdok || ! concept) { 556 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " 557 "unable to find NONLIN.SOURCE in recipe --- ignored.\n"); 558 } else { 559 psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept); 560 if (! conceptValueItem) { 561 psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s " 562 "for non-linearity correction --- ignored.\n", concept); 563 } else if (conceptValueItem->type != PS_DATA_STRING) { 564 psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as" 565 " expected for non-linearity correction --- ignored.\n", 566 concept); 567 } else { 568 psString conceptValue = conceptValueItem->data.V; 569 // Get the value of the concept 570 psMetadataItem *optionItem = psMetadataLookup(options, conceptValue); 571 if (!optionItem) { 572 psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA" 573 " --- ignored.\n", conceptValue); 574 } else if (optionItem->type == PS_DATA_VECTOR) { 575 // These are the polynomial coefficients 576 psVector *coeff = optionItem->data.V; // The coefficient vector 577 if (coeff->type.type != PS_TYPE_F64) { 578 psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); 579 coeff = temp; 580 } 581 // Polynomial correction 582 psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, 583 PS_POLYNOMIAL_ORD); 584 psFree(correction->coeff); 585 correction->coeff = psMemIncrRefCounter(coeff->data.F64); 586 (void)pmNonLinearityPolynomial(inputReadout, correction); 587 psFree(coeff); 588 psFree(correction); 589 } else if (optionItem->type == PS_DATA_STRING) { 590 // This is a filename 591 psString tableName = optionItem->data.V; // The filename 592 psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0); 593 int numLines = 0; // Number of lines read from table 594 if ((numLines = psLookupTableRead(table)) <= 0) { 595 psErrorStackPrint(stderr, "Unable to read non-linearity " 596 "correction file %s --- ignored\n", 597 tableName); 598 } else { 372 599 #ifdef PRODUCTION 373 psFitsClose(biasFile);600 (void)pmNonLinearityLookup(inputReadout, table); 374 601 #else 375 psFree(biasFile); 376 #endif 377 } 378 379 if (doDark) { 602 printf("XXX: Non-linearity correction from lookup table not " 603 "yet implemented.\n"); 604 #endif 605 } 606 psFree(table); 607 } else { 608 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction " 609 "desired but unable to interpret NONLIN.DATA for %s" 610 " --- ignored\n", conceptValue); 611 } 612 } 613 } 614 } 615 break; 616 default: 617 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " 618 "NONLIN.DATA is of invalid type --- ignored.\n"); 619 } 620 } 621 } // Non-linearity correction 622 623 // Bias, dark and overscan subtraction are all merged. 624 625 // Changed the definition of pmSubtractBias to do the dark subtraction as well, but 626 // it's not in there yet. Once it is, we can get rid of "pedestal". 627 380 628 #ifdef PRODUCTION 381 psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark 629 psImage *pedestal = NULL; // Pedestal image (bias + scaled dark) 630 if (doDark) { 631 // Dark time for input image 632 float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME"); 633 if (inputTime <= 0) { 634 psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n", 635 inputTime); 636 exit(EXIT_FAILURE); 637 } 638 639 pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*", 640 psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); 641 } 642 if (doBias) { 643 if (pedestal) { 644 if (biasImage->numRows != darkImage->numRows || 645 biasImage->numCols != darkImage->numCols) { 646 psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: " 647 "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows, 648 darkImage->numCols, darkImage->numRows); 649 exit(EXIT_FAILURE); 650 } 651 (void)psBinaryOp(pedestal, pedestal, "+", biasImage); 652 } else { 653 pedestal = psMemIncrRefCounter(biasImage); 654 } 655 } 656 #endif 657 if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { 658 // Need to get the read direction 659 int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR"); 660 if (readdir == 1) { 661 // psmodule-0.8.0 has confused PM_OVERSCAN_ROWS and PM_OVERSCAN_COLUMNS; bug 608. 662 #ifdef PRODUCTION 663 overscanMode = PM_OVERSCAN_ROWS; 382 664 #else 383 psFits *darkFile = psFitsAlloc(darkName); 384 #endif 385 psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark 386 if (! pmConfigValidateCamera(camera, darkHeader)) { 387 psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName); 388 exit(EXIT_FAILURE); 389 } 390 psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName); 391 if (! pmFPARead(dark, darkFile, darkHeader, database)) { 392 psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName); 393 exit(EXIT_FAILURE); 394 } 395 psFree(darkHeader); 665 overscanMode = PM_OVERSCAN_COLUMNS; 666 #endif 667 } else if (readdir == 2) { 396 668 #ifdef PRODUCTION 397 psFitsClose(darkFile);669 overscanMode = PM_OVERSCAN_COLUMNS; 398 670 #else 399 psFree(darkFile); 400 #endif 401 } 402 403 if (doMask) { 671 overscanMode = PM_OVERSCAN_ROWS; 672 #endif 673 } else { 674 psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n", 675 readdir); 676 overscanMode = PM_OVERSCAN_ROWS; 677 } 678 } 679 680 if (doBias || doOverscan || doDark) { 681 psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions 404 682 #ifdef PRODUCTION 405 psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask 683 (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, 684 overscanStats, overscanBins, overscanFitType, pedestal); 406 685 #else 407 psFits *maskFile = psFitsAlloc(maskName); 408 #endif 409 psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask 410 if (! pmConfigValidateCamera(camera, maskHeader)) { 411 psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName); 412 exit(EXIT_FAILURE); 413 } 414 psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName); 415 if (! pmFPARead(mask, maskFile, maskHeader, database)) { 416 psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName); 417 exit(EXIT_FAILURE); 418 } 419 psFree(maskHeader); 420 421 #ifdef PRODUCTION 422 psFitsClose(maskFile); 423 #else 424 psFree(maskFile); 425 #endif 426 } 427 428 if (doFlat) { 429 #ifdef PRODUCTION 430 psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat 431 #else 432 psFits *flatFile = psFitsAlloc(flatName); 433 #endif 434 if (! flatFile) { 435 psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName); 436 exit(EXIT_FAILURE); 437 } 438 psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat 439 if (! pmConfigValidateCamera(camera, flatHeader)) { 440 psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName); 441 exit(EXIT_FAILURE); 442 } 443 psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName); 444 if (! pmFPARead(flat, flatFile, flatHeader, database)) { 445 psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName); 446 exit(EXIT_FAILURE); 447 } 448 psFree(flatHeader); 449 #ifdef PRODUCTION 450 psFitsClose(flatFile); 451 #else 452 psFree(flatFile); 453 #endif 454 } 455 456 psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2")); 457 458 psArray *inputChips = input->chips; // Array of chips in input image 459 psArray *biasChips = bias->chips; // Array of chips in bias image 460 psArray *darkChips = dark->chips; // Array of chips in dark image 461 psArray *maskChips = mask->chips; // Array of chips in mask image 462 psArray *flatChips = flat->chips; // Array of chips in flat image 463 for (int i = 0; i < inputChips->n; i++) { 464 pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image 465 pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image 466 pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image 467 pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image 468 pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image 469 470 if (! inputChip->valid) { 471 continue; 472 } 473 474 psArray *inputCells = inputChip->cells; // Array of cells in input image 475 psArray *biasCells = biasChip->cells; // Array of cells in bias image 476 psArray *darkCells = darkChip->cells; // Array of cells in dark image 477 psArray *maskCells = maskChip->cells; // Array of cells in mask image 478 psArray *flatCells = flatChip->cells; // Array of cells in flat image 479 480 for (int j = 0; j < inputCells->n; j++) { 481 pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image 482 pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image 483 pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag 484 pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image 485 pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image 486 487 if (! inputCell->valid) { 488 continue; 489 } 490 491 psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image 492 if (doBias && biasCell->readouts->n > 1) { 493 psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" 494 " used."); 495 } 496 if (doDark && darkCell->readouts->n > 1) { 497 psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" 498 " used."); 499 } 500 if (doMask && maskCell->readouts->n > 1) { 501 psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" 502 " used."); 503 } 504 if (doFlat && flatCell->readouts->n > 1) { 505 psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be" 506 " used."); 507 } 508 509 pmReadout *biasReadout = NULL; // Bias readout 510 pmReadout *maskReadout = NULL; // Mask readout 511 pmReadout *flatReadout = NULL; // Flat readout 512 psImage *biasImage = NULL; // Bias pixels 513 psImage *darkImage = NULL; // Dark pixels 514 float darkTime = 1.0; // Dark time for dark image 515 psImage *maskImage = NULL; // Mask pixels 516 517 if (doBias) { 518 biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image 519 biasImage = biasReadout->image; 520 } 521 if (doDark) { 522 pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image 523 darkImage = darkReadout->image; 524 darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME"); 525 if (darkTime <= 0.0) { 526 psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime); 527 exit(EXIT_FAILURE); 528 } 529 } 530 if (doMask) { 531 maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image 532 maskImage = maskReadout->image; 533 } 534 if (doFlat) { 535 flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image 536 } 537 538 for (int k = 0; k < inputReadouts->n; k++) { 539 pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image 540 psImage *inputImage = inputReadout->image; // The actual image 541 542 // Mask bad pixels 543 if (doMask) { 544 float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION"); 545 546 // Need to change this later to grow the mask by the size of the convolution kernel 547 (void)pmMaskBadPixels(inputReadout, maskImage, 548 PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation, 549 PM_MASK_TRAP, 0); 550 } 551 552 // Non-linearity correction 553 if (doNonLin) { 554 psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA"); 555 if (! dataItem) { 556 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to " 557 "find NONLIN.DATA in recipe --- ignored.\n"); 558 } else { 559 switch (dataItem->type) { 560 case PS_DATA_VECTOR: 561 { 562 // These are the polynomial coefficients 563 psVector *coeff = dataItem->data.V; // The coefficient vector 564 if (coeff->type.type != PS_TYPE_F64) { 565 psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version 566 coeff = temp; 567 } 568 psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, 569 PS_POLYNOMIAL_ORD); 570 psFree(correction->coeff); 571 correction->coeff = psMemIncrRefCounter(coeff->data.F64); 572 (void)pmNonLinearityPolynomial(inputReadout, correction); 573 psFree(coeff); 574 psFree(correction); 575 } 576 break; 577 case PS_DATA_STRING: 578 { 579 // This is a filename 580 psString name = dataItem->data.V; // Filename 581 psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0); 582 if (psLookupTableRead(table) <= 0) { 583 psErrorStackPrint(stderr, "Unable to read non-linearity correction file " 584 "%s --- ignored\n", name); 585 } else { 586 #ifdef PRODUCTION 587 (void)pmNonLinearityLookup(inputReadout, table); 588 #else 589 psVector *influx = table->values->data[0]; 590 psVector *outflux = table->values->data[1]; 591 (void)pmNonLinearityLookup(inputReadout, table->values->data[0], 592 table->values->data[1]); 593 #endif 594 } 595 psFree(table); 596 } 597 break; 598 case PS_DATA_METADATA: 599 { 600 // This is a menu 601 psMetadata *options = dataItem->data.V; // Options with concept values as keys 602 bool mdok = false; // Success of MD lookup 603 psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE"); 604 if (! mdok || ! concept) { 605 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " 606 "unable to find NONLIN.SOURCE in recipe --- ignored.\n"); 607 } else { 608 psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept); 609 if (! conceptValueItem) { 610 psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s " 611 "for non-linearity correction --- ignored.\n", concept); 612 } else if (conceptValueItem->type != PS_DATA_STRING) { 613 psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as" 614 " expected for non-linearity correction --- ignored.\n", 615 concept); 616 } else { 617 psString conceptValue = conceptValueItem->data.V; 618 // Get the value of the concept 619 psMetadataItem *optionItem = psMetadataLookup(options, conceptValue); 620 if (!optionItem) { 621 psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA" 622 " --- ignored.\n", conceptValue); 623 } else if (optionItem->type == PS_DATA_VECTOR) { 624 // These are the polynomial coefficients 625 psVector *coeff = optionItem->data.V; // The coefficient vector 626 if (coeff->type.type != PS_TYPE_F64) { 627 psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); 628 coeff = temp; 629 } 630 // Polynomial correction 631 psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, 632 PS_POLYNOMIAL_ORD); 633 psFree(correction->coeff); 634 correction->coeff = psMemIncrRefCounter(coeff->data.F64); 635 (void)pmNonLinearityPolynomial(inputReadout, correction); 636 psFree(coeff); 637 psFree(correction); 638 } else if (optionItem->type == PS_DATA_STRING) { 639 // This is a filename 640 psString tableName = optionItem->data.V; // The filename 641 psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0); 642 int numLines = 0; // Number of lines read from table 643 if ((numLines = psLookupTableRead(table)) <= 0) { 644 psErrorStackPrint(stderr, "Unable to read non-linearity " 645 "correction file %s --- ignored\n", 646 tableName); 647 } else { 648 #ifdef PRODUCTION 649 (void)pmNonLinearityLookup(inputReadout, table); 650 #else 651 printf("XXX: Non-linearity correction from lookup table not " 652 "yet implemented.\n"); 653 #endif 654 } 655 psFree(table); 656 } else { 657 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction " 658 "desired but unable to interpret NONLIN.DATA for %s" 659 " --- ignored\n", conceptValue); 660 } 661 } 662 } 663 } 664 break; 665 default: 666 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " 667 "NONLIN.DATA is of invalid type --- ignored.\n"); 668 } 669 } 670 } // Non-linearity correction 671 672 // Bias, dark and overscan subtraction are all merged. 673 674 // Changed the definition of pmSubtractBias to do the dark subtraction as well, but 675 // it's not in there yet. Once it is, we can get rid of "pedestal". 676 677 #ifdef PRODUCTION 678 psImage *pedestal = NULL; // Pedestal image (bias + scaled dark) 679 if (doDark) { 680 // Dark time for input image 681 float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME"); 682 if (inputTime <= 0) { 683 psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n", 684 inputTime); 685 exit(EXIT_FAILURE); 686 } 687 688 pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*", 689 psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); 690 } 691 if (doBias) { 692 if (pedestal) { 693 if (biasImage->numRows != darkImage->numRows || 694 biasImage->numCols != darkImage->numCols) { 695 psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: " 696 "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows, 697 darkImage->numCols, darkImage->numRows); 698 exit(EXIT_FAILURE); 699 } 700 (void)psBinaryOp(pedestal, pedestal, "+", biasImage); 701 } else { 702 pedestal = psMemIncrRefCounter(biasImage); 703 } 704 } 705 #endif 706 if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { 707 // Need to get the read direction 708 int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR"); 709 if (readdir == 1) { 710 overscanMode = PM_OVERSCAN_ROWS; 711 } else if (readdir == 2) { 712 overscanMode = PM_OVERSCAN_COLUMNS; 713 } else { 714 psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n", 715 readdir); 716 overscanMode = PM_OVERSCAN_ROWS; 717 } 718 } 719 720 if (doBias || doOverscan || doDark) { 721 psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions 722 #ifdef PRODUCTION 723 (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, 724 overscanStats, overscanBins, overscanFitType, pedestal); 725 #else 726 (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, 727 overscanStats, overscanBins, overscanFitType, biasReadout); 728 #endif 729 psFree(inputOverscans); 730 731 // Output overscan fit results, if required 732 if (doOverscan) { 733 if (! overscanFit) { 734 psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n"); 735 } else { 736 switch (overscanFitType) { 737 case PM_FIT_POLYNOMIAL: 738 { 739 psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial 740 psString coeffs = NULL; // String containing the coefficients 741 for (int i = 0; i < poly->nX; i++) { 742 psStringAppend(&coeffs, "%e ", poly->coeff[i]); 743 } 744 psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n", 745 coeffs); 746 psFree(coeffs); 747 } 748 break; 749 case PM_FIT_SPLINE: 750 { 751 psSpline1D *spline = (psSpline1D*)overscanFit; // The spline 752 psString coeffs = NULL; // String containing the coefficients 753 for (int i = 0; i < spline->n; i++) { 754 psPolynomial1D *poly = spline->spline[i]; // i-th polynomial 755 psStringAppend(&coeffs, "%d: ", i); 756 for (int j = 0; j < poly->nX; j++) { 757 psStringAppend(&coeffs, "%e ", poly->coeff[i]); 758 } 759 psStringAppend(&coeffs, "\n"); 760 } 761 psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n", 762 coeffs); 763 psFree(coeffs); 764 } 765 break; 766 case PM_FIT_NONE: 767 break; 768 default: 769 psAbort(__func__, "Should never get here!!!\n"); 770 } 771 } 772 } 773 } 774 775 // XXX: Temporary: until the other functions are altered to do this themselves. 776 // Trim, so that flat, fringe etc computations are faster. 686 (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, 687 overscanStats, overscanBins, overscanFitType, biasReadout); 688 #endif 689 psFree(inputOverscans); 690 691 // Output overscan fit results, if required 692 if (doOverscan) { 693 if (! overscanFit) { 694 psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n"); 695 } else { 696 switch (overscanFitType) { 697 case PM_FIT_POLYNOMIAL: 698 { 699 psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial 700 psString coeffs = NULL; // String containing the coefficients 701 for (int i = 0; i < poly->nX; i++) { 702 psStringAppend(&coeffs, "%e ", poly->coeff[i]); 703 } 704 psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n", 705 coeffs); 706 psFree(coeffs); 707 } 708 break; 709 case PM_FIT_SPLINE: 710 { 711 psSpline1D *spline = (psSpline1D*)overscanFit; // The spline 712 psString coeffs = NULL; // String containing the coefficients 713 for (int i = 0; i < spline->n; i++) { 714 psPolynomial1D *poly = spline->spline[i]; // i-th polynomial 715 psStringAppend(&coeffs, "%d: ", i); 716 for (int j = 0; j < poly->nX; j++) { 717 psStringAppend(&coeffs, "%e ", poly->coeff[i]); 718 } 719 psStringAppend(&coeffs, "\n"); 720 } 721 psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n", 722 coeffs); 723 psFree(coeffs); 724 } 725 break; 726 case PM_FIT_NONE: 727 break; 728 default: 729 psAbort(__func__, "Should never get here!!!\n"); 730 } 731 } 732 } 733 } 734 735 // XXX: Temporary: until the other functions are altered to do this themselves. 736 // Trim, so that flat, fringe etc computations are faster. 777 737 #if 0 778 738 #ifndef PRODUCTION 779 psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC");780 psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion);781 psFree(inputReadout->image);782 inputReadout->image = trimmed;783 #endif 784 #endif 785 // Flat-field correction786 if (doFlat) {787 psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n");739 psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC"); 740 psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion); 741 psFree(inputReadout->image); 742 inputReadout->image = trimmed; 743 #endif 744 #endif 745 // Flat-field correction 746 if (doFlat) { 747 psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n"); 788 748 #ifndef PRODUCTION 789 psImage *dummyImage = psImageAlloc(inputReadout->image->numCols,790 inputReadout->image->numRows,791 PS_TYPE_U8);792 pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1);793 (void)pmFlatField(inputReadout, dummyMask, flatReadout);794 psFree(dummyMask);795 psFree(dummyImage);796 #endif 797 }798 799 // Fringe subtraction800 if (doFringe) {749 psImage *dummyImage = psImageAlloc(inputReadout->image->numCols, 750 inputReadout->image->numRows, 751 PS_TYPE_U8); 752 pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1); 753 (void)pmFlatField(inputReadout, dummyMask, flatReadout); 754 psFree(dummyMask); 755 psFree(dummyImage); 756 #endif 757 } 758 759 // Fringe subtraction 760 if (doFringe) { 801 761 #if 0 802 // Placeholder 803 pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, 804 psU8 maskVal, int binFactor, psStats *stats, float clipSD); 805 #endif 806 } 807 808 // Source identification and photometry 809 if (doSource) { 810 } 811 812 // Astrometry 813 if (doAstrom) { 814 } 815 816 817 } // Iterating over readouts 818 } // Iterating over cells 819 820 // Testing pmChipMosaic 821 { 822 psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip 823 psFits *mosaicFile = psFitsAlloc("mosaic.fits"); // FITS file to which to write 824 psFitsWriteImage(mosaicFile, NULL, mosaic, 0); 825 psFree(mosaicFile); 826 psFree(mosaic); 827 } 828 762 // Placeholder 763 pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, 764 psU8 maskVal, int binFactor, psStats *stats, float clipSD); 765 #endif 766 } 767 768 // Source identification and photometry 769 if (doSource) { 770 } 771 772 // Astrometry 773 if (doAstrom) { 774 } 775 776 777 } // Iterating over readouts 778 } // Iterating over cells 779 780 #if 0 781 { 782 // Testing pmChipMosaic 783 psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip 784 psFits *mosaicFile = psFitsOpen("mosaic.fits", "w"); // FITS file to which to write 785 psFitsWriteImage(mosaicFile, NULL, mosaic, 0); 786 psFitsClose(mosaicFile); 787 psFree(mosaic); 788 } 789 #endif 790 829 791 } // Iterating over chips 830 792 … … 832 794 833 795 796 #if 0 797 { 798 // Testing pmFPAMorph 799 int nBad = 0; 800 psMetadata *newCamera = psMetadataConfigParse(NULL, &nBad, "megacam_splice.config", true); 801 pmFPA *newFPA = pmFPAConstruct(newCamera); 802 pmFPAMorph(newFPA, input, true, 0, 0); 803 //pmFPAPrint(newFPA); 804 psFits *newFile = psFitsOpen("morph.fits", "w"); 805 pmFPAWrite(newFile, newFPA, database); 806 psFitsClose(newFile); 807 psFree(newFPA); 808 } 809 #endif 810 834 811 #if 1 835 // Testing pmFPAMorph836 psMetadata *newCamera = psMetadataConfigParse(NULL, NULL, "megacam_splice.config", true);837 pmFPA *newFPA = pmFPAConstruct(newCamera);838 pmFPAMorph(newFPA, input, true, 0, 0);839 psFits *newFile = psFitsAlloc("morph.fits");840 pmFPAWrite(newFile, newFPA, database);841 psFree(newFile);842 psFree(newFPA);843 #endif844 845 #if 0846 812 // Write the output 847 813 pmFPAWrite(outputFile, input, database); … … 852 818 psLogMsg("phase2", PS_LOG_INFO, "Output completed after %f sec.\n", psTimerMark("phase2")); 853 819 854 #ifdef PRODUCTION855 820 psFitsClose(outputFile); 856 #else857 psFree(outputFile);858 #endif859 821 860 822 psFree(arguments); … … 874 836 psTimerStop(); 875 837 876 #if 1838 #if 0 877 839 psMemCheckCorruption(true); 878 840 psMemBlock **leaks = NULL; // List of leaks … … 881 843 #if 1 882 844 for (int i = 0; i < nLeaks; i++) { 883 psMemBlock *mb = leaks[i];884 printf("Memory leak %ld (%ld):\n"885 "\tFile %s, line %d, size %d\n"886 "\tPosts: %x %x %x\n",887 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,888 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));845 psMemBlock *mb = leaks[i]; 846 printf("Memory leak %ld (%ld):\n" 847 "\tFile %s, line %d, size %d\n" 848 "\tPosts: %x %x %x\n", 849 mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, 850 *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); 889 851 } 890 852 #endif -
trunk/archive/scripts/src/phase2/phase2.config
r5651 r5786 3 3 # List of tasks to perform 4 4 MASK BOOL FALSE # Mask bad pixels 5 NONLIN BOOL FALSE # Non-linearity correction6 BIAS BOOL FALSE # Bias subtraction5 NONLIN BOOL TRUE # Non-linearity correction 6 BIAS BOOL TRUE # Bias subtraction 7 7 DARK BOOL FALSE # Dark subtraction 8 OVERSCAN BOOL FALSE # Overscan subtraction9 FLAT BOOL FALSE # Flat-field normalisation8 OVERSCAN BOOL TRUE # Overscan subtraction 9 FLAT BOOL TRUE # Flat-field normalisation 10 10 FRINGE BOOL FALSE # Fringe subtraction 11 11 SOURCE BOOL FALSE # Source identification and photometry -
trunk/archive/scripts/src/phase2/pmChipMosaic.c
r5621 r5786 1 1 #include <stdio.h> 2 #include <assert.h> 2 3 3 4 #include "pslib.h" … … 14 15 } 15 16 16 psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic 17 int xBinChip, int yBinChip // Binning of mosaic image in x and y 17 // Mosaic multiple images, with flips, binning and offsets 18 psImage *p_pmImageMosaic(const psArray *source, // Images to splice in 19 const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? 20 const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of 21 // source images 22 int xBinTarget, int yBinTarget, // Binning in x and y of target images 23 const psVector *x0, const psVector *y0 // Offsets for source images on target 18 24 ) 19 25 { 20 psArray *cells = chip->cells; // The array of cells 21 psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates 22 psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates 23 psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x 24 psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y 25 psVector *xParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in x 26 psVector *yParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in y 27 psArray *trimsec = psArrayAlloc(cells->n); // Trim section 26 assert(source); 27 assert(xFlip && xFlip->type.type == PS_TYPE_U8); 28 assert(yFlip && yFlip->type.type == PS_TYPE_U8); 29 assert(xBinSource && xBinSource->type.type == PS_TYPE_S32); 30 assert(yBinSource && yBinSource->type.type == PS_TYPE_S32); 31 assert(x0 && x0->type.type == PS_TYPE_S32); 32 assert(y0 && y0->type.type == PS_TYPE_S32); 28 33 29 34 // Get the maximum extent of the mosaic image … … 32 37 int yMin = INT_MAX; 33 38 int yMax = 0; 34 for (int i = 0; i < cells->n; i++) { 35 pmCell *cell = cells->data[i]; // The cell of interest 36 x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); 37 y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); 38 xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); 39 yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); 40 xParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); 41 yParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); 42 trimsec->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC")); 43 psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, ((psRegion*)trimsec->data[i])->x0, 44 ((psRegion*)trimsec->data[i])->x1, ((psRegion*)trimsec->data[i])->y0, 45 ((psRegion*)trimsec->data[i])->y1); 46 // Size of cell in x and y 47 int nx = (int)(((psRegion*)trimsec->data[i])->x1 - (int)((psRegion*)trimsec->data[i])->x0); 48 int ny = (int)(((psRegion*)trimsec->data[i])->y1 - (int)((psRegion*)trimsec->data[i])->y0); 49 psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i], 50 x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * nx, y0->data.S32[i], 51 y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny); 39 for (int i = 0; i < source->n; i++) { 40 psImage *image = source->data[i]; // The image of interest 52 41 53 COMPARE(x0->data.S32[i], xMin, xMax); 54 COMPARE(y0->data.S32[i], yMin, yMax); 55 // Subtract the parity to get the inclusive limit (not exclusive) 56 COMPARE(x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * nx - xParity->data.S32[i], xMin, xMax); 57 COMPARE(y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny - yParity->data.S32[i], yMin, yMax); 42 assert(image->type.type == PS_TYPE_F32); // Only implemented for F32 images so far. 43 44 // Size of cell in x and y 45 int xParity = xFlip->data.U8[i] ? -1 : 1; 46 int yParity = yFlip->data.U8[i] ? -1 : 1; 47 psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i], 48 x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols, y0->data.S32[i], 49 y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows); 50 51 COMPARE(x0->data.S32[i], xMin, xMax); 52 COMPARE(y0->data.S32[i], yMin, yMax); 53 // Subtract the parity to get the inclusive limit (not exclusive) 54 COMPARE(x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols - xParity, xMin, xMax); 55 COMPARE(y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows - yParity, yMin, yMax); 58 56 } 59 57 60 58 // Set up the image 61 59 // Since both upper and lower values are inclusive, we need to add one to the size 62 float xSize = (float)(xMax - xMin + 1) / (float)xBin Chip;60 float xSize = (float)(xMax - xMin + 1) / (float)xBinTarget; 63 61 if (xSize - (int)xSize > 0) { 64 xSize += 1;62 xSize += 1; 65 63 } 66 float ySize = (float)(yMax - yMin + 1) / (float)yBin Chip;64 float ySize = (float)(yMax - yMin + 1) / (float)yBinTarget; 67 65 if (ySize - (int)ySize > 0) { 68 ySize += 1;66 ySize += 1; 69 67 } 70 68 71 psTrace(__func__, 3, " Mosaicked chipwill be %dx%d\n", (int)xSize, (int)ySize);69 psTrace(__func__, 3, "Spliced image will be %dx%d\n", (int)xSize, (int)ySize); 72 70 psImage *mosaic = psImageAlloc((int)xSize, (int)ySize, PS_TYPE_F32); // The mosaic image 73 71 psImageInit(mosaic, 0.0); 74 72 75 // Set the image 73 // Next pass through the images to do the mosaicking 74 for (int i = 0; i < source->n; i++) { 75 psImage *image = source->data[i]; // The image of interest 76 if (xBinSource->data.S32[i] == xBinTarget && yBinSource->data.S32[i] == yBinTarget && 77 xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) { 78 // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason 79 psImageOverlaySection(mosaic, image, x0->data.S32[i], y0->data.S32[i], "+"); 80 } else { 81 // We have to do the hard work ourself 82 for (int y = 0; y < image->numRows; y++) { 83 int yParity = yFlip->data.U8[i] ? -1 : 1; 84 float yTargetBase = (y0->data.S32[i] + yParity * yBinSource->data.S32[i] * y) / yBinTarget; 85 for (int x = 0; x < image->numCols; x++) { 86 int xParity = xFlip->data.U8[i] ? -1 : 1; 87 float xTargetBase = (x0->data.S32[i] + xParity * xBinSource->data.S32[i] * x) / 88 xBinTarget; 89 90 // In case the original image is binned but the mosaic is not, we need to fill in the 91 // values in the mosaic. 92 for (int j = 0; j < yBinSource->data.S32[i]; j++) { 93 int yTarget = (int)(yTargetBase + yParity * (float)j / (float)yBinTarget); 94 for (int i = 0; i < xBinSource->data.S32[i]; i++) { 95 int xTarget = (int)(xTargetBase + xParity * (float)i / (float)xBinTarget); 96 97 mosaic->data.F32[yTarget][xTarget] += image->data.F32[y][x]; 98 } 99 } // Iterating over mosaic image for binned input image 100 } 101 } // Iterating over input image 102 } 103 } 104 105 return mosaic; 106 } 107 108 psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic 109 int xBinChip, int yBinChip // Binning of mosaic image in x and y 110 ) 111 { 112 psArray *cells = chip->cells; // The array of cells 113 psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked 114 psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates 115 psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates 116 psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x 117 psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y 118 psVector *xFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in x? 119 psVector *yFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in y? 120 121 // Set up the required inputs 76 122 for (int i = 0; i < cells->n; i++) { 77 pmCell *cell = cells->data[i]; // The cell of interest 78 psArray *readouts = cell->readouts; // The array of readouts 79 if (readouts->n > 1) { 80 psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " 81 "be mosaicked.\n", i); 82 } 83 psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic 84 psImage *trimmed = psImageSubset(image, *(psRegion*)(trimsec->data[i])); // Trimmed image (no overscan) 123 pmCell *cell = cells->data[i]; // The cell of interest 124 x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); 125 y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); 126 xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); 127 yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); 128 int xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); 129 int yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); 130 if (xParity == 1) { 131 xFlip->data.U8[i] = 0; 132 } else if (xParity == -1) { 133 xFlip->data.U8[i] = 1; 134 } else { 135 psLogMsg(__func__, PS_LOG_WARN, "The x parity of cell %d is not +/- 1 (it's %d) --- " 136 "assuming +1.\n", i, xParity); 137 xFlip->data.U8[i] = 0; 138 } 139 if (yParity == 1) { 140 yFlip->data.U8[i] = 0; 141 } else if (yParity == -1) { 142 yFlip->data.U8[i] = 1; 143 } else { 144 psLogMsg(__func__, PS_LOG_WARN, "The y parity of cell %d is not +/- 1 (it's %d) --- " 145 "assuming +1.\n", i, yParity); 146 yFlip->data.U8[i] = 0; 147 } 85 148 86 if (xBin->data.S32[i] == xBinChip && yBin->data.S32[i] == yBinChip && xParity->data.S32[i] == 1 && 87 yParity->data.S32[i] == 1) { 88 // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason 89 psImageOverlaySection(mosaic, trimmed, x0->data.S32[i], y0->data.S32[i], "+");90 } else { 91 // We have to do the hard work ourself 92 for (int y = 0; y < trimmed->numRows; y++) { 93 float yTargetBase = (y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * y) / 94 yBinChip; 95 for (int x = 0; x < trimmed->numCols; x++) { 96 float xTargetBase = (x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * x) / 97 xBinChip; 149 // Trim the image to get rid of the overscan 150 psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); 151 psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, trimsec->x0, trimsec->x1, 152 trimsec->y0, trimsec->y1); 153 psArray *readouts = cell->readouts; // The array of readouts 154 if (readouts->n > 1) { 155 psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will " 156 "be mosaicked.\n", i); 157 } 158 psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic 159 images->data[i] = psImageSubset(image, *trimsec); // Trimmed image 160 } 98 161 99 // In case the original image is binned but the mosaic is not, we need to fill in the 100 // values in the mosaic. 101 for (int j = 0; j < yBin->data.S32[i]; j++) { 102 int yTarget = (int)(yTargetBase + yParity->data.S32[i] * (float)j / (float)yBinChip); 103 for (int i = 0; i < xBin->data.S32[i]; i++) { 104 int xTarget = (int)(xTargetBase + 105 xParity->data.S32[i] * (float)i / (float)xBinChip); 162 // Mosaic the images together and we're done 163 psImage *mosaic = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0); 106 164 107 mosaic->data.F32[yTarget][xTarget] += trimmed->data.F32[y][x]; 108 } 109 } // Iterating over mosaic image for binned input image 110 } 111 } // Iterating over input image 112 } 113 psFree(trimmed); 114 } // Iterating over cells 115 165 // Clean up 116 166 psFree(x0); 117 167 psFree(y0); 118 168 psFree(xBin); 119 169 psFree(yBin); 120 psFree(x Parity);121 psFree(y Parity);122 psFree( trimsec);170 psFree(xFlip); 171 psFree(yFlip); 172 psFree(images); 123 173 124 174 return mosaic; -
trunk/archive/scripts/src/phase2/pmChipMosaic.h
r5621 r5786 5 5 #include "pmFPA.h" 6 6 7 psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic 8 int xBinChip, int yBinChip // Binning of mosaic image in x and y 7 // Mosaic multiple images, with flips, binning and offsets 8 psImage *p_pmImageMosaic(const psArray *source, // Images to splice in 9 const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? 10 const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of 11 // source images 12 int xBinTarget, int yBinTarget, // Binning in x and y of target images 13 const psVector *x0, const psVector *y0 // Offsets for source images on target 14 ); 15 16 // Mosaic all the cells in a chip together (neglecting the overscans) 17 psImage *pmChipMosaic(pmChip *chip, // Chip to mosaic 18 int xBinChip, int yBinChip // Binning of mosaic image in x and y 9 19 ); 10 20 -
trunk/archive/scripts/src/phase2/pmFPA.c
r5633 r5786 77 77 fpa->chips = psArrayAlloc(0); 78 78 79 fpa->phu = NULL; 79 80 fpa->hdu = NULL; 80 81 … … 89 90 90 91 psFree(fpa->concepts); 91 psFree( (psPtr)fpa->camera);92 psFree(fpa->camera); 92 93 psFree(fpa->chips); 93 94 95 psFree(fpa->phu); 94 96 psFree(fpa->hdu); 95 97 } -
trunk/archive/scripts/src/phase2/pmFPA.h
r5651 r5786 19 19 typedef struct { 20 20 // Astrometric transformations 21 psPlaneDistort* fromTangentPlane; // Transformation from tangent plane to focal plane22 psPlaneDistort* toTangentPlane; // Transformation from focal plane to tangent plane23 psProjection *projection; // Projection from tangent plane to sky21 psPlaneDistort* fromTangentPlane; // Transformation from tangent plane to focal plane 22 psPlaneDistort* toTangentPlane; // Transformation from focal plane to tangent plane 23 psProjection *projection; // Projection from tangent plane to sky 24 24 // Information 25 psMetadata *concepts; // Cache for PS concepts26 psMetadata *analysis; // FPA-level analysis metadata27 const psMetadata *camera; // Camera configuration28 psArray *chips; // The chips29 p_pmHDU *hdu; // FITS data30 psMetadata *phu; // Primary Header25 psMetadata *concepts; // Cache for PS concepts 26 psMetadata *analysis; // FPA-level analysis metadata 27 const psMetadata *camera; // Camera configuration 28 psArray *chips; // The chips 29 p_pmHDU *hdu; // FITS data 30 psMetadata *phu; // Primary Header 31 31 } pmFPA; 32 32 33 33 typedef struct { 34 34 // Offset specifying position on focal plane 35 int col0; // Offset from the left of FPA.36 int row0; // Offset from the bottom of FPA.35 int col0; // Offset from the left of FPA. 36 int row0; // Offset from the bottom of FPA. 37 37 // Astrometric transformations 38 psPlaneTransform* toFPA; // Transformation from chip to FPA coordinates39 psPlaneTransform* fromFPA; // Transformation from FPA to chip coordinates38 psPlaneTransform* toFPA; // Transformation from chip to FPA coordinates 39 psPlaneTransform* fromFPA; // Transformation from FPA to chip coordinates 40 40 // Information 41 psMetadata *concepts; // Cache for PS concepts42 psMetadata *analysis; // Chip-level analysis metadata43 psArray *cells; // The cells (referred to by name)44 pmFPA *parent; // Parent FPA45 bool valid; // Do we bother about reading and working with this chip?41 psMetadata *concepts; // Cache for PS concepts 42 psMetadata *analysis; // Chip-level analysis metadata 43 psArray *cells; // The cells (referred to by name) 44 pmFPA *parent; // Parent FPA 45 bool valid; // Do we bother about reading and working with this chip? 46 46 p_pmHDU *hdu; // FITS data 47 47 } pmChip; … … 49 49 typedef struct { 50 50 // Offset specifying position on chip 51 int col0; // Offset from the left of chip.52 int row0; // Offset from the bottom of chip.51 int col0; // Offset from the left of chip. 52 int row0; // Offset from the bottom of chip. 53 53 // Astrometric transformations 54 psPlaneTransform* toChip; // Transformations from cell to chip coordinates55 psPlaneTransform* toFPA; // Transformations from cell to FPA coordinates56 psPlaneTransform* toSky; // Transformations from cell to sky coordinates54 psPlaneTransform* toChip; // Transformations from cell to chip coordinates 55 psPlaneTransform* toFPA; // Transformations from cell to FPA coordinates 56 psPlaneTransform* toSky; // Transformations from cell to sky coordinates 57 57 // Information 58 psMetadata *concepts; // Cache for PS concepts59 psMetadata *camera; // Camera information60 psMetadata *analysis; // Cell-level analysis metadata61 psArray *readouts; // The readouts (referred to by number)62 pmChip *parent; // Parent chip63 bool valid; // Do we bother about reading and working with this cell?64 p_pmHDU *hdu; // FITS data58 psMetadata *concepts; // Cache for PS concepts 59 psMetadata *camera; // Camera information 60 psMetadata *analysis; // Cell-level analysis metadata 61 psArray *readouts; // The readouts (referred to by number) 62 pmChip *parent; // Parent chip 63 bool valid; // Do we bother about reading and working with this cell? 64 p_pmHDU *hdu; // FITS data 65 65 } pmCell; 66 66 67 67 typedef struct { 68 68 // Position on the cell 69 int col0; // Offset from the left of cell.70 int row0; // Offset from the bottom of cell.71 int colBins; // Amount of binning in x-dimension and parity (from sign)72 int rowBins; // Amount of binning in y-dimension and parity (from sign)69 int col0; // Offset from the left of cell. 70 int row0; // Offset from the bottom of cell. 71 int colBins; // Amount of binning in x-dimension and parity (from sign) 72 int rowBins; // Amount of binning in y-dimension and parity (from sign) 73 73 // Information 74 psImage *image; // Imaging area of readout75 psImage *mask; // Mask for image76 psImage *weight; // Weight for image74 psImage *image; // Imaging area of readout 75 psImage *mask; // Mask for image 76 psImage *weight; // Weight for image 77 77 #if 0 78 psList *bias; // List of bias section (sub-)images78 psList *bias; // List of bias section (sub-)images 79 79 #endif 80 psMetadata *analysis; // Readout-level analysis metadata81 pmCell *parent; // Parent cell80 psMetadata *analysis; // Readout-level analysis metadata 81 pmCell *parent; // Parent cell 82 82 } pmReadout; 83 83 -
trunk/archive/scripts/src/phase2/pmFPAConceptsGet.c
r5633 r5786 14 14 15 15 psMetadataItem *p_pmFPAConceptGetFromCamera(pmCell *cell, // The cell 16 const char *concept // Name of concept16 const char *concept // Name of concept 17 17 ) 18 18 { 19 19 if (cell) { 20 psMetadata *camera = cell->camera;// Camera data21 // Need the CELL.NAME first, which should be in the cell->concepts since creation of the cell22 bool mdStatus = true;// Status of MD lookup23 // Finally, the info that we want24 psMetadataItem *item = psMetadataLookup(camera, concept);25 return item;20 psMetadata *camera = cell->camera; // Camera data 21 // Need the CELL.NAME first, which should be in the cell->concepts since creation of the cell 22 bool mdStatus = true; // Status of MD lookup 23 // Finally, the info that we want 24 psMetadataItem *item = psMetadataLookup(camera, concept); 25 return item; 26 26 } 27 27 return NULL; 28 } 28 } 29 29 30 30 31 31 psMetadataItem *p_pmFPAConceptGetFromHeader(pmFPA *fpa, // The FPA that contains the chip 32 pmChip *chip, // The chip that contains the cell33 pmCell *cell, // The cell34 const char *concept // Name of concept35 ) 36 { 37 bool mdStatus = true; // Status of MD lookup32 pmChip *chip, // The chip that contains the cell 33 pmCell *cell, // The cell 34 const char *concept // Name of concept 35 ) 36 { 37 bool mdStatus = true; // Status of MD lookup 38 38 psMetadata *translation = psMetadataLookupMD(&mdStatus, fpa->camera, "TRANSLATION"); // FITS translation 39 39 if (! mdStatus) { 40 psError(PS_ERR_IO, false, "Unable to find TRANSLATION in camera configuration.\n");41 return NULL;40 psError(PS_ERR_IO, false, "Unable to find TRANSLATION in camera configuration.\n"); 41 return NULL; 42 42 } 43 43 … … 45 45 const char *keyword = psMetadataLookupString(&mdStatus, translation, concept); 46 46 if (mdStatus && strlen(keyword) > 0) { 47 // We have a FITS header to look up --- search each level48 if (cell && cell->hdu) {49 psMetadataItem *cellItem = psMetadataLookup(cell->hdu->header, keyword);50 if (cellItem) {51 // XXX: Need to clean up before returning52 return cellItem;53 }54 }55 56 if (chip && chip->hdu) {57 psMetadataItem *chipItem = psMetadataLookup(chip->hdu->header, keyword);58 if (chipItem) {59 // XXX: Need to clean up before returning60 return chipItem;61 }62 }63 64 if (fpa->hdu) {65 psMetadataItem *fpaItem = psMetadataLookup(fpa->hdu->header, keyword);66 if (fpaItem) {67 // XXX: Need to clean up before returning68 return fpaItem;69 }70 }71 72 if (fpa->phu) {73 psMetadataItem *fpaItem = psMetadataLookup(fpa->phu, keyword);74 if (fpaItem) {75 // XXX: Need to clean up before returning76 return fpaItem;77 }78 }47 // We have a FITS header to look up --- search each level 48 if (cell && cell->hdu) { 49 psMetadataItem *cellItem = psMetadataLookup(cell->hdu->header, keyword); 50 if (cellItem) { 51 // XXX: Need to clean up before returning 52 return cellItem; 53 } 54 } 55 56 if (chip && chip->hdu) { 57 psMetadataItem *chipItem = psMetadataLookup(chip->hdu->header, keyword); 58 if (chipItem) { 59 // XXX: Need to clean up before returning 60 return chipItem; 61 } 62 } 63 64 if (fpa->hdu) { 65 psMetadataItem *fpaItem = psMetadataLookup(fpa->hdu->header, keyword); 66 if (fpaItem) { 67 // XXX: Need to clean up before returning 68 return fpaItem; 69 } 70 } 71 72 if (fpa->phu) { 73 psMetadataItem *fpaItem = psMetadataLookup(fpa->phu, keyword); 74 if (fpaItem) { 75 // XXX: Need to clean up before returning 76 return fpaItem; 77 } 78 } 79 79 } 80 80 … … 86 86 // Look for a default 87 87 psMetadataItem *p_pmFPAConceptGetFromDefault(pmFPA *fpa, // The FPA that contains the chip 88 pmChip *chip, // The chip that contains the cell89 pmCell *cell, // The cell90 const char *concept // Name of concept91 ) 92 { 93 bool mdOK = true; // Status of MD lookup88 pmChip *chip, // The chip that contains the cell 89 pmCell *cell, // The cell 90 const char *concept // Name of concept 91 ) 92 { 93 bool mdOK = true; // Status of MD lookup 94 94 psMetadata *defaults = psMetadataLookupMD(&mdOK, fpa->camera, "DEFAULTS"); 95 95 if (! mdOK) { 96 psError(PS_ERR_IO, false, "Unable to find DEFAULTS in camera configuration.\n");97 return NULL;96 psError(PS_ERR_IO, false, "Unable to find DEFAULTS in camera configuration.\n"); 97 return NULL; 98 98 } 99 99 100 100 psMetadataItem *defItem = psMetadataLookup(defaults, concept); 101 101 if (defItem) { 102 if (defItem->type == PS_DATA_METADATA) {103 // A dependent default104 psTrace(__func__, 7, "Evaluating dependent default....\n");105 psMetadata *dependents = defItem->data.V; // The list of dependents106 // Find out what it depends on107 psString dependName = psStringCopy(concept);108 psStringAppend(&dependName, ".DEPEND");109 psString dependsOn = psMetadataLookupString(&mdOK, defaults, dependName);110 if (! mdOK) {111 psError(PS_ERR_IO, false, "Unable to find %s in camera configuration for dependent default"112 " --- ignored\n", dependName);113 // XXX: Need to clean up before returning114 return NULL;115 }116 psFree(dependName);117 // Find the value of the dependent concept118 psMetadataItem *depItem = p_pmFPAConceptGetFromHeader(fpa, chip, cell, dependsOn);119 if (! depItem) {120 psError(PS_ERR_IO, true, "Unable to find value for %s (required for %s)\n", dependsOn,121 concept);122 return NULL;123 }124 if (depItem->type != PS_DATA_STRING) {125 psError(PS_ERR_IO, true, "Value of %s is not of type string, as required for dependency"126 " --- ignored.\n", dependsOn);127 }128 129 defItem = psMetadataLookup(dependents, depItem->data.V);// This is now what we were after130 }102 if (defItem->type == PS_DATA_METADATA) { 103 // A dependent default 104 psTrace(__func__, 7, "Evaluating dependent default....\n"); 105 psMetadata *dependents = defItem->data.V; // The list of dependents 106 // Find out what it depends on 107 psString dependName = psStringCopy(concept); 108 psStringAppend(&dependName, ".DEPEND"); 109 psString dependsOn = psMetadataLookupString(&mdOK, defaults, dependName); 110 if (! mdOK) { 111 psError(PS_ERR_IO, false, "Unable to find %s in camera configuration for dependent default" 112 " --- ignored\n", dependName); 113 // XXX: Need to clean up before returning 114 return NULL; 115 } 116 psFree(dependName); 117 // Find the value of the dependent concept 118 psMetadataItem *depItem = p_pmFPAConceptGetFromHeader(fpa, chip, cell, dependsOn); 119 if (! depItem) { 120 psError(PS_ERR_IO, true, "Unable to find value for %s (required for %s)\n", dependsOn, 121 concept); 122 return NULL; 123 } 124 if (depItem->type != PS_DATA_STRING) { 125 psError(PS_ERR_IO, true, "Value of %s is not of type string, as required for dependency" 126 " --- ignored.\n", dependsOn); 127 } 128 129 defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after 130 } 131 131 } 132 132 133 133 // XXX: Need to clean up before returning 134 return defItem; // defItem is either NULL or points to what was desired134 return defItem; // defItem is either NULL or points to what was desired 135 135 } 136 136 … … 139 139 // XXX: Not tested 140 140 psMetadataItem *p_pmFPAConceptGetFromDB(pmFPA *fpa, // The FPA that contains the chip 141 pmChip *chip, // The chip that contains the cell142 pmCell *cell, // The cell143 psDB *db,// DB handle144 const char *concept // Name of concept141 pmChip *chip, // The chip that contains the cell 142 pmCell *cell, // The cell 143 psDB *db, // DB handle 144 const char *concept // Name of concept 145 145 ) 146 146 { 147 147 if (! db) { 148 // No database initialised149 return NULL;150 } 151 152 bool mdStatus = true; // Status of MD lookup148 // No database initialised 149 return NULL; 150 } 151 152 bool mdStatus = true; // Status of MD lookup 153 153 psMetadata *database = psMetadataLookupMD(&mdStatus, fpa->camera, "DATABASE"); 154 154 if (! mdStatus) { 155 // No error, because not everyone needs to use the DB156 return NULL;155 // No error, because not everyone needs to use the DB 156 return NULL; 157 157 } 158 158 159 159 psMetadata *dbLookup = psMetadataLookupMD(&mdStatus, database, concept); 160 160 if (dbLookup) { 161 const char *tableName = psMetadataLookupString(&mdStatus, dbLookup, "TABLE"); // Name of the table162 const char *colName = psMetadataLookupString(&mdStatus, dbLookup, "COLUMN"); // Name of the column163 const char *givenCols = psMetadataLookupString(&mdStatus, dbLookup, "GIVENDBCOL"); // Name of "where"164 // columns165 const char *givenPS = psMetadataLookupString(&mdStatus, dbLookup, "GIVENPS"); // Values for "where"166 // columns167 168 // Now, need to get the "given"s169 if (strlen(givenCols) || strlen(givenPS)) {170 psList *cols = papSplit(givenCols, ",;"); // List of column names171 psList *values = papSplit(givenPS, ",;"); // List of value names for the columns172 psMetadata *selection = psMetadataAlloc(); // The stuff to select in the DB173 if (cols->n != values->n) {174 psLogMsg(__func__, PS_LOG_WARN, "The GIVENDBCOL and GIVENPS entries for %s do not have "175 "the same number of entries --- ignored.\n", concept);176 } else {177 // Iterators for the lists178 psListIterator *colsIter = psListIteratorAlloc(cols, PS_LIST_HEAD, false);179 psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false);180 char *column = NULL;// Name of the column181 while (column = psListGetAndIncrement(colsIter)) {182 char *name = psListGetAndIncrement(valuesIter); // Name for the value183 if (!strlen(column) || !strlen(name)) {184 psLogMsg(__func__, PS_LOG_WARN, "One of the columns or value names for %s is "185 " empty --- ignored.\n", concept);186 } else {187 // Search for the value name188 psMetadataItem *item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, name);189 if (! item) {190 item = p_pmFPAConceptGetFromDefault(fpa, chip, cell, name);191 }192 if (! item) {193 psLogMsg(__func__, PS_LOG_ERROR, "Unable to find the value name %s for DB "194 " lookup on %s --- ignored.\n", name, concept);195 } else {196 // We need to create a new psMetadataItem. I don't think we can't simply hack197 // the existing one, since that could conceivably cause memory leaks198 psMetadataItem *newItem = psMetadataItemAlloc(concept, item->type,199 item->comment, item->data.V);200 psMetadataAddItem(selection, newItem, PS_LIST_TAIL, PS_META_REPLACE);201 psFree(newItem);202 }203 }204 psFree(name);205 psFree(column);206 } // Iterating through the columns207 psFree(colsIter);208 psFree(valuesIter);209 210 psArray *dbResult = psDBSelectRows(db, tableName, selection, 2); // Lookup result211 // Note that we use limit=2 in order to test if there are multiple rows returned212 213 psMetadataItem *result = NULL; // The final result of the DB lookup214 if (dbResult->n == 0) {215 psLogMsg(__func__, PS_LOG_WARN, "Unable to find any rows in DB for %s --- ignored\n", 216 concept);217 } else {218 if (dbResult-> n > 1) {219 psLogMsg(__func__, PS_LOG_WARN, "Multiple rows returned in DB lookup for %s --- "220 " using the first one only.\n", concept);221 }222 result = (psMetadataItem*)dbResult->data[0];223 }224 // XXX: Need to clean up before returning225 return result;226 }227 psFree(cols);228 psFree(values);229 }161 const char *tableName = psMetadataLookupString(&mdStatus, dbLookup, "TABLE"); // Name of the table 162 const char *colName = psMetadataLookupString(&mdStatus, dbLookup, "COLUMN"); // Name of the column 163 const char *givenCols = psMetadataLookupString(&mdStatus, dbLookup, "GIVENDBCOL"); // Name of "where" 164 // columns 165 const char *givenPS = psMetadataLookupString(&mdStatus, dbLookup, "GIVENPS"); // Values for "where" 166 // columns 167 168 // Now, need to get the "given"s 169 if (strlen(givenCols) || strlen(givenPS)) { 170 psList *cols = papSplit(givenCols, ",;"); // List of column names 171 psList *values = papSplit(givenPS, ",;"); // List of value names for the columns 172 psMetadata *selection = psMetadataAlloc(); // The stuff to select in the DB 173 if (cols->n != values->n) { 174 psLogMsg(__func__, PS_LOG_WARN, "The GIVENDBCOL and GIVENPS entries for %s do not have " 175 "the same number of entries --- ignored.\n", concept); 176 } else { 177 // Iterators for the lists 178 psListIterator *colsIter = psListIteratorAlloc(cols, PS_LIST_HEAD, false); 179 psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false); 180 char *column = NULL; // Name of the column 181 while (column = psListGetAndIncrement(colsIter)) { 182 char *name = psListGetAndIncrement(valuesIter); // Name for the value 183 if (!strlen(column) || !strlen(name)) { 184 psLogMsg(__func__, PS_LOG_WARN, "One of the columns or value names for %s is " 185 " empty --- ignored.\n", concept); 186 } else { 187 // Search for the value name 188 psMetadataItem *item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, name); 189 if (! item) { 190 item = p_pmFPAConceptGetFromDefault(fpa, chip, cell, name); 191 } 192 if (! item) { 193 psLogMsg(__func__, PS_LOG_ERROR, "Unable to find the value name %s for DB " 194 " lookup on %s --- ignored.\n", name, concept); 195 } else { 196 // We need to create a new psMetadataItem. I don't think we can't simply hack 197 // the existing one, since that could conceivably cause memory leaks 198 psMetadataItem *newItem = psMetadataItemAlloc(concept, item->type, 199 item->comment, item->data.V); 200 psMetadataAddItem(selection, newItem, PS_LIST_TAIL, PS_META_REPLACE); 201 psFree(newItem); 202 } 203 } 204 psFree(name); 205 psFree(column); 206 } // Iterating through the columns 207 psFree(colsIter); 208 psFree(valuesIter); 209 210 psArray *dbResult = psDBSelectRows(db, tableName, selection, 2); // Lookup result 211 // Note that we use limit=2 in order to test if there are multiple rows returned 212 213 psMetadataItem *result = NULL; // The final result of the DB lookup 214 if (dbResult->n == 0) { 215 psLogMsg(__func__, PS_LOG_WARN, "Unable to find any rows in DB for %s --- ignored\n", 216 concept); 217 } else { 218 if (dbResult-> n > 1) { 219 psLogMsg(__func__, PS_LOG_WARN, "Multiple rows returned in DB lookup for %s --- " 220 " using the first one only.\n", concept); 221 } 222 result = (psMetadataItem*)dbResult->data[0]; 223 } 224 // XXX: Need to clean up before returning 225 return result; 226 } 227 psFree(cols); 228 psFree(values); 229 } 230 230 } // Doing the "given"s. 231 231 … … 236 236 // Concept lookup 237 237 psMetadataItem *p_pmFPAConceptGet(pmFPA *fpa, // The FPA 238 pmChip *chip,// The chip239 pmCell *cell,// The cell240 psDB *db, // DB handle241 const char *concept // Concept name238 pmChip *chip,// The chip 239 pmCell *cell, // The cell 240 psDB *db, // DB handle 241 const char *concept // Concept name 242 242 ) 243 243 { … … 245 245 psMetadataItem *item = p_pmFPAConceptGetFromCamera(cell, concept); 246 246 if (! item) { 247 item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, concept);247 item = p_pmFPAConceptGetFromHeader(fpa, chip, cell, concept); 248 248 } 249 249 if (! item) { … … 257 257 258 258 259 void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA260 pmChip *chip, // The chip261 pmCell *cell, // The cell262 psDB *db,// DB handle263 psMetadata *concepts, // The concepts MD264 const char *name, // Name of the concept265 const char *comment // Comment for concept259 void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA 260 pmChip *chip, // The chip 261 pmCell *cell, // The cell 262 psDB *db, // DB handle 263 psMetadata *concepts, // The concepts MD 264 const char *name, // Name of the concept 265 const char *comment // Comment for concept 266 266 ) 267 267 { … … 269 269 float value = NAN; 270 270 if (item) { 271 switch (item->type) {272 case PS_DATA_F32:273 value = item->data.F32;274 break;275 case PS_DATA_F64:276 // Assume it's OK to truncate to floating point from double277 value = (float)item->data.F64;278 break;279 case PS_DATA_S32:280 // Promote to float281 value = (float)item->data.S32;282 break;283 default:284 psError(PS_ERR_IO, true, "Concept %s (%s) is not of floating point type (%x) --- treating as "285 "undefined.\n", name, comment, item->type);286 }271 switch (item->type) { 272 case PS_DATA_F32: 273 value = item->data.F32; 274 break; 275 case PS_DATA_F64: 276 // Assume it's OK to truncate to floating point from double 277 value = (float)item->data.F64; 278 break; 279 case PS_DATA_S32: 280 // Promote to float 281 value = (float)item->data.S32; 282 break; 283 default: 284 psError(PS_ERR_IO, true, "Concept %s (%s) is not of floating point type (%x) --- treating as " 285 "undefined.\n", name, comment, item->type); 286 } 287 287 } else { 288 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment);288 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); 289 289 } 290 290 psTrace(__func__, 7, "Adding %s (%s): %f\n", name, comment, value); 291 291 292 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F32 , comment, value);293 } 294 295 void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA296 pmChip *chip, // The chip297 pmCell *cell, // The cell298 psDB *db,// DB handle299 psMetadata *concepts, // The concepts MD300 const char *name, // Name of the concept301 const char *comment // Comment for concept292 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F32 | PS_META_REPLACE, comment, value); 293 } 294 295 void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA 296 pmChip *chip, // The chip 297 pmCell *cell, // The cell 298 psDB *db, // DB handle 299 psMetadata *concepts, // The concepts MD 300 const char *name, // Name of the concept 301 const char *comment // Comment for concept 302 302 ) 303 303 { … … 305 305 double value = NAN; 306 306 if (item) { 307 switch (item->type) {308 case PS_TYPE_F64:309 value = item->data.F64;310 break;311 case PS_TYPE_F32:312 // Promote to double313 value = (double)item->data.F32;314 break;315 case PS_TYPE_S32:316 // Promote to double317 value = (double)item->data.S32;318 break;319 default:320 psError(PS_ERR_IO, true, "Concept %s (%s) is not of double-precision floating point type (%x) "321 "--- treating as undefined.\n", name, comment, item->type);322 }307 switch (item->type) { 308 case PS_TYPE_F64: 309 value = item->data.F64; 310 break; 311 case PS_TYPE_F32: 312 // Promote to double 313 value = (double)item->data.F32; 314 break; 315 case PS_TYPE_S32: 316 // Promote to double 317 value = (double)item->data.S32; 318 break; 319 default: 320 psError(PS_ERR_IO, true, "Concept %s (%s) is not of double-precision floating point type (%x) " 321 "--- treating as undefined.\n", name, comment, item->type); 322 } 323 323 } else { 324 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment);324 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); 325 325 } 326 326 psTrace(__func__, 7, "Adding %s (%s): %f\n", name, comment, value); 327 327 328 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F64 , comment, value);329 } 330 331 void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA332 pmChip *chip,// The chip333 pmCell *cell,// The cell334 psDB *db,// DB handle335 psMetadata *concepts, // The concepts MD336 const char *name, // Name of the concept337 const char *comment // Comment for concept328 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_F64 | PS_META_REPLACE, comment, value); 329 } 330 331 void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA 332 pmChip *chip, // The chip 333 pmCell *cell, // The cell 334 psDB *db, // DB handle 335 psMetadata *concepts, // The concepts MD 336 const char *name, // Name of the concept 337 const char *comment // Comment for concept 338 338 ) 339 339 { … … 341 341 int value = 0; 342 342 if (item) { 343 switch (item->type) {344 case PS_TYPE_S32:345 value = item->data.S32;346 break;347 case PS_TYPE_F32:348 psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F32 --- converting.\n",349 name, comment);350 value = (int)item->data.F32;351 break;352 case PS_TYPE_F64:353 psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F64 --- converting.\n",354 name, comment);355 value = (int)item->data.F64;356 break;357 default:358 psError(PS_ERR_IO, true, "Concept %s (%s) is not of integer type (%x) --- treating as "359 "undefined.\n", name, comment, item->type);360 }343 switch (item->type) { 344 case PS_TYPE_S32: 345 value = item->data.S32; 346 break; 347 case PS_TYPE_F32: 348 psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F32 --- converting.\n", 349 name, comment); 350 value = (int)item->data.F32; 351 break; 352 case PS_TYPE_F64: 353 psLogMsg(__func__, PS_LOG_WARN, "Concept %s (%s) should be S32, but is F64 --- converting.\n", 354 name, comment); 355 value = (int)item->data.F64; 356 break; 357 default: 358 psError(PS_ERR_IO, true, "Concept %s (%s) is not of integer type (%x) --- treating as " 359 "undefined.\n", name, comment, item->type); 360 } 361 361 } else { 362 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment);362 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); 363 363 } 364 364 psTrace(__func__, 7, "Adding %s (%s): %d\n", name, comment, value); 365 366 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_S32 , comment, value);365 366 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_TYPE_S32 | PS_META_REPLACE, comment, value); 367 367 } 368 368 369 369 void p_pmFPAConceptGetString(pmFPA *fpa, // The FPA 370 pmChip *chip, // The chip371 pmCell *cell, // The cell372 psDB *db,// DB handle373 psMetadata *concepts, // The concepts MD374 const char *name, // Name of the concept375 const char *comment // Comment for concept370 pmChip *chip, // The chip 371 pmCell *cell, // The cell 372 psDB *db, // DB handle 373 psMetadata *concepts, // The concepts MD 374 const char *name, // Name of the concept 375 const char *comment // Comment for concept 376 376 ) 377 377 { … … 379 379 psString value = NULL; 380 380 if (item) { 381 switch (item->type) {382 case PS_DATA_STRING:383 value = psMemIncrRefCounter(item->data.V);384 break;385 case PS_DATA_F32:386 psStringAppend(&value, "%f", item->data.F32);387 break;388 case PS_DATA_S32:389 psStringAppend(&value, "%d", item->data.S32);390 break;391 default:392 psError(PS_ERR_IO, true, "Concept %s (%s) is not of string type (%x) --- treating as "393 "undefined.\n", name, comment, item->type);394 }381 switch (item->type) { 382 case PS_DATA_STRING: 383 value = psMemIncrRefCounter(item->data.V); 384 break; 385 case PS_DATA_F32: 386 psStringAppend(&value, "%f", item->data.F32); 387 break; 388 case PS_DATA_S32: 389 psStringAppend(&value, "%d", item->data.S32); 390 break; 391 default: 392 psError(PS_ERR_IO, true, "Concept %s (%s) is not of string type (%x) --- treating as " 393 "undefined.\n", name, comment, item->type); 394 } 395 395 } else { 396 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment);397 value = psStringCopy("");396 psError(PS_ERR_IO, true, "Concept %s (%s) is not defined.\n", name, comment); 397 value = psStringCopy(""); 398 398 } 399 399 psTrace(__func__, 7, "Adding %s (%s): %s\n", name, comment, value); 400 400 401 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_DATA_STRING , comment, value);401 psMetadataAdd(concepts, PS_LIST_TAIL, name, PS_DATA_STRING | PS_META_REPLACE, comment, value); 402 402 psFree(value); 403 403 } … … 409 409 410 410 // Ingest concepts for the FPA 411 void pmFPAIngestConcepts(pmFPA *fpa, // The FPA412 psDB *db// DB handle411 void pmFPAIngestConcepts(pmFPA *fpa, // The FPA 412 psDB *db // DB handle 413 413 ) 414 414 { 415 415 if (! fpa->concepts) { 416 fpa->concepts = psMetadataAlloc();416 fpa->concepts = psMetadataAlloc(); 417 417 } 418 418 … … 436 436 // FPA.RA 437 437 { 438 double ra = NAN;// The RA439 psMetadataItem *raItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.RA"); // The FPA.RA item440 if (raItem) {441 switch (raItem->type) {442 case PS_TYPE_F32:443 ra = raItem->data.F32;444 break;445 case PS_TYPE_F64:446 ra = raItem->data.F64;447 break;448 case PS_DATA_STRING:449 // Sexagesimal format450 {451 int big, medium;452 float small;453 // XXX: Upgrade path is to allow dd:mm.mmm454 if (sscanf(raItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 &&455 sscanf(raItem->data.V, "%d %d %f", &big, &medium, &small) != 3) {456 psError(PS_ERR_IO, true, "Cannot interpret FPA.RA: %s\n", raItem->data.V);457 break;458 }459 ra = abs(big) + (float)medium/60.0 + small/3600.0;460 if (big < 0) {461 ra *= -1.0;462 }463 }464 break;465 default:466 psError(PS_ERR_IO, true, "FPA.RA is of an unexpected type: %x\n", raItem->type);467 }468 469 // How to interpret the RA470 const psMetadata *camera = fpa->camera; // Camera configuration data471 bool mdok = true;// Status of MD lookup472 psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS");473 if (mdok && formats) {474 psString raFormat = psMetadataLookupString(&mdok, formats, "FPA.RA");475 if (mdok && strlen(raFormat) > 0) {476 if (strcasecmp(raFormat, "HOURS") == 0) {477 ra *= M_PI / 12.0;478 } else if (strcasecmp(raFormat, "DEGREES") == 0) {479 ra *= M_PI / 180.0;480 } else if (strcasecmp(raFormat, "RADIANS") == 0) {481 // No action required482 } else {483 psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.RA in FORMATS (%s) --- assuming"484 " HOURS.\n");485 ra *= M_PI / 12.0;486 }487 } else {488 psError(PS_ERR_IO, false, "Unable to find FPA.RA in FORMATS --- assuming HOURS.\n");489 ra *= M_PI / 12.0;490 }491 } else {492 psError(PS_ERR_IO, false, "Unable to find FORMAT metadata in camera configuration --- "493 "assuming format for FPA.RA is HOURS.\n");494 ra *= M_PI / 12.0;495 }496 } else {497 psError(PS_ERR_IO, false, "Couldn't find FPA.RA.\n");498 }499 500 psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_DATA_F64, 501 "Right Ascension of the boresight (radians)", ra);438 double ra = NAN; // The RA 439 psMetadataItem *raItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.RA"); // The FPA.RA item 440 if (raItem) { 441 switch (raItem->type) { 442 case PS_TYPE_F32: 443 ra = raItem->data.F32; 444 break; 445 case PS_TYPE_F64: 446 ra = raItem->data.F64; 447 break; 448 case PS_DATA_STRING: 449 // Sexagesimal format 450 { 451 int big, medium; 452 float small; 453 // XXX: Upgrade path is to allow dd:mm.mmm 454 if (sscanf(raItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && 455 sscanf(raItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { 456 psError(PS_ERR_IO, true, "Cannot interpret FPA.RA: %s\n", raItem->data.V); 457 break; 458 } 459 ra = abs(big) + (float)medium/60.0 + small/3600.0; 460 if (big < 0) { 461 ra *= -1.0; 462 } 463 } 464 break; 465 default: 466 psError(PS_ERR_IO, true, "FPA.RA is of an unexpected type: %x\n", raItem->type); 467 } 468 469 // How to interpret the RA 470 const psMetadata *camera = fpa->camera; // Camera configuration data 471 bool mdok = true; // Status of MD lookup 472 psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); 473 if (mdok && formats) { 474 psString raFormat = psMetadataLookupString(&mdok, formats, "FPA.RA"); 475 if (mdok && strlen(raFormat) > 0) { 476 if (strcasecmp(raFormat, "HOURS") == 0) { 477 ra *= M_PI / 12.0; 478 } else if (strcasecmp(raFormat, "DEGREES") == 0) { 479 ra *= M_PI / 180.0; 480 } else if (strcasecmp(raFormat, "RADIANS") == 0) { 481 // No action required 482 } else { 483 psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.RA in FORMATS (%s) --- assuming" 484 " HOURS.\n"); 485 ra *= M_PI / 12.0; 486 } 487 } else { 488 psError(PS_ERR_IO, false, "Unable to find FPA.RA in FORMATS --- assuming HOURS.\n"); 489 ra *= M_PI / 12.0; 490 } 491 } else { 492 psError(PS_ERR_IO, false, "Unable to find FORMAT metadata in camera configuration --- " 493 "assuming format for FPA.RA is HOURS.\n"); 494 ra *= M_PI / 12.0; 495 } 496 } else { 497 psError(PS_ERR_IO, false, "Couldn't find FPA.RA.\n"); 498 } 499 500 psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_DATA_F64 | PS_META_REPLACE, 501 "Right Ascension of the boresight (radians)", ra); 502 502 503 503 } … … 505 505 // FPA.DEC 506 506 { 507 double dec = NAN;// The DEC508 psMetadataItem *decItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.DEC"); // The FPA.DEC item509 if (decItem) {510 switch (decItem->type) {511 case PS_TYPE_F32:512 dec = decItem->data.F32;513 break;514 case PS_TYPE_F64:515 dec = decItem->data.F64;516 break;517 case PS_DATA_STRING:518 // Sexagesimal format519 {520 int big, medium;521 float small;522 // XXX: Upgrade path is to allow dd:mm.mmm523 if (sscanf(decItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 &&524 sscanf(decItem->data.V, "%d %d %f", &big, &medium, &small) != 3) {525 psError(PS_ERR_IO, true, "Cannot interpret FPA.DEC: %s\n", decItem->data.V);526 break;527 }528 dec = abs(big) + (float)medium/60.0 + small/3600.0;529 if (big < 0) {530 dec *= -1.0;531 }532 }533 break;534 default:535 psError(PS_ERR_IO, true, "FPA.DEC is of an unexpected type: %x\n", decItem->type);536 }537 538 // How to interpret the DEC539 const psMetadata *camera = fpa->camera; // Camera configuration data540 bool mdok = true;// Status of MD lookup541 psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS");542 if (mdok && formats) {543 psString decFormat = psMetadataLookupString(&mdok, formats, "FPA.DEC");544 if (mdok && strlen(decFormat) > 0) {545 if (strcasecmp(decFormat, "HOURS") == 0) {546 dec *= M_PI / 12.0;547 } else if (strcasecmp(decFormat, "DEGREES") == 0) {548 dec *= M_PI / 180.0;549 } else if (strcasecmp(decFormat, "RADIANS") == 0) {550 // No action required551 } else {552 psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.DEC in FORMATS (%s) --- "553 "assuming DEGREES.\n");554 dec *= M_PI / 180.0;555 }556 } else {557 psError(PS_ERR_IO, false, "Unable to find FPA.DEC in FORMATS --- assuming DEGREES.\n");558 dec *= M_PI / 180.0;559 }560 } else {561 psError(PS_ERR_IO, false, "Unable to find FORMATS metadata in camera configuration --- "562 "assuming format for FPA.DEC is DEGREES.\n");563 dec *= M_PI / 180.0;564 }565 } else {566 psError(PS_ERR_IO, false, "Couldn't find FPA.DEC.\n");567 }568 569 psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_DATA_F64, 570 "Declination of the boresight (radians)", dec);507 double dec = NAN; // The DEC 508 psMetadataItem *decItem = p_pmFPAConceptGet(fpa, NULL, NULL, db, "FPA.DEC"); // The FPA.DEC item 509 if (decItem) { 510 switch (decItem->type) { 511 case PS_TYPE_F32: 512 dec = decItem->data.F32; 513 break; 514 case PS_TYPE_F64: 515 dec = decItem->data.F64; 516 break; 517 case PS_DATA_STRING: 518 // Sexagesimal format 519 { 520 int big, medium; 521 float small; 522 // XXX: Upgrade path is to allow dd:mm.mmm 523 if (sscanf(decItem->data.V, "%d:%d:%f", &big, &medium, &small) != 3 && 524 sscanf(decItem->data.V, "%d %d %f", &big, &medium, &small) != 3) { 525 psError(PS_ERR_IO, true, "Cannot interpret FPA.DEC: %s\n", decItem->data.V); 526 break; 527 } 528 dec = abs(big) + (float)medium/60.0 + small/3600.0; 529 if (big < 0) { 530 dec *= -1.0; 531 } 532 } 533 break; 534 default: 535 psError(PS_ERR_IO, true, "FPA.DEC is of an unexpected type: %x\n", decItem->type); 536 } 537 538 // How to interpret the DEC 539 const psMetadata *camera = fpa->camera; // Camera configuration data 540 bool mdok = true; // Status of MD lookup 541 psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); 542 if (mdok && formats) { 543 psString decFormat = psMetadataLookupString(&mdok, formats, "FPA.DEC"); 544 if (mdok && strlen(decFormat) > 0) { 545 if (strcasecmp(decFormat, "HOURS") == 0) { 546 dec *= M_PI / 12.0; 547 } else if (strcasecmp(decFormat, "DEGREES") == 0) { 548 dec *= M_PI / 180.0; 549 } else if (strcasecmp(decFormat, "RADIANS") == 0) { 550 // No action required 551 } else { 552 psLogMsg(__func__, PS_LOG_WARN, "Don't understand FPA.DEC in FORMATS (%s) --- " 553 "assuming DEGREES.\n"); 554 dec *= M_PI / 180.0; 555 } 556 } else { 557 psError(PS_ERR_IO, false, "Unable to find FPA.DEC in FORMATS --- assuming DEGREES.\n"); 558 dec *= M_PI / 180.0; 559 } 560 } else { 561 psError(PS_ERR_IO, false, "Unable to find FORMATS metadata in camera configuration --- " 562 "assuming format for FPA.DEC is DEGREES.\n"); 563 dec *= M_PI / 180.0; 564 } 565 } else { 566 psError(PS_ERR_IO, false, "Couldn't find FPA.DEC.\n"); 567 } 568 569 psMetadataAdd(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_DATA_F64 | PS_META_REPLACE, 570 "Declination of the boresight (radians)", dec); 571 571 572 572 } … … 577 577 578 578 // Ingest concepts for the chip 579 bool pmChipIngestConcepts(pmChip *chip, // The chip580 psDB *db// DB handle581 ) 582 { 583 pmFPA *fpa = chip->parent; // The parent FPA579 bool pmChipIngestConcepts(pmChip *chip, // The chip 580 psDB *db // DB handle 581 ) 582 { 583 pmFPA *fpa = chip->parent; // The parent FPA 584 584 585 585 if (! chip->concepts) { 586 chip->concepts = psMetadataAlloc();586 chip->concepts = psMetadataAlloc(); 587 587 } 588 588 … … 591 591 // Pau. 592 592 } 593 594 595 // Add corrective to a position --- in case the user wants FORTRAN indexing (like the FITS standard...) 596 static bool correctPosition(pmFPA *fpa, // FPA, contains the camera configuration 597 pmCell *cell, // Cell containing the concept to correct 598 const char *conceptName // Name of concept to correct 599 ) 600 { 601 bool mdok = false; // Result of MD lookup 602 psMetadata *formats = psMetadataLookupMD(&mdok, fpa->camera, "FORMATS"); 603 if (mdok && formats) { 604 psString format = psMetadataLookupString(&mdok, formats, conceptName); 605 if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { 606 psMetadataItem *valueItem = psMetadataLookup(cell->concepts, conceptName); 607 valueItem->data.S32 -= 1; 608 } 609 } 610 return true; 611 } 612 613 bool p_pmCellIngestConcept(pmFPA *fpa, // The FPA 614 pmChip *chip, // The chip 615 pmCell *cell, // The cell 616 psDB *db, // DB handle 617 const char *concept // Name of the concept 618 ) 619 { 620 if (strcmp(concept, "CELL.GAIN") == 0) { 621 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.GAIN", "CCD gain (e/count)"); 622 } else if (strcmp(concept, "CELL.READNOISE") == 0) { 623 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.READNOISE", "CCD read noise (e)"); 624 } else if (strcmp(concept, "CELL.SATURATION") == 0) { 625 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.SATURATION", 626 "Saturation level (ADU)"); 627 } else if (strcmp(concept, "CELL.BAD") == 0) { 628 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.BAD", "Bad level (ADU)"); 629 } else if (strcmp(concept, "CELL.XPARITY") == 0) { 630 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.XPARITY", 631 "Orientation in x compared to the rest of the FPA"); 632 } else if (strcmp(concept, "CELL.YPARITY") == 0) { 633 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.YPARITY", 634 "Orientation in y compared to the rest of the FPA"); 635 } else if (strcmp(concept, "CELL.READDIR") == 0) { 636 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.READDIR", 637 "Read direction: 1=row, 2=col"); 638 } else if (strcmp(concept, "CELL.EXPOSURE") == 0) { // used to be READOUT.EXPOSURE 639 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.EXPOSURE", "Exposure time (sec)"); 640 } else if (strcmp(concept, "CELL.DARKTIME") == 0) { // used to be READOUT.DARKTIME 641 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.DARKTIME", 642 "Time since CCD flush (sec)"); 643 // These take some extra work 644 } else if (strcmp(concept, "CELL.TRIMSEC") == 0) { 645 psRegion *trimsec = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed by value) 646 647 psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC"); 648 if (! secItem) { 649 psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.\n"); 650 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 651 } else if (secItem->type != PS_DATA_STRING) { 652 psError(PS_ERR_IO, true, "CELL.TRIMSEC is not of type STR (%x)\n", secItem->type); 653 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 654 } else { 655 psString section = secItem->data.V; // The section string 656 657 psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC.SOURCE"); 658 if (! sourceItem) { 659 psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.SOURCE.\n"); 660 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 661 } else if (sourceItem->type != PS_DATA_STRING) { 662 psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); 663 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 664 } else { 665 psString source = sourceItem->data.V; // The source string 666 667 if (strcasecmp(source, "VALUE") == 0) { 668 *trimsec = psRegionFromString(section); 669 } else if (strcasecmp(source, "HEADER") == 0) { 670 psMetadata *header = NULL; // The FITS header 671 if (cell->hdu) { 672 header = cell->hdu->header; 673 } else if (chip->hdu) { 674 header = chip->hdu->header; 675 } else if (fpa->hdu) { 676 header = fpa->hdu->header; 677 } 678 if (! header) { 679 psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); 680 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 681 } else { 682 bool mdok = true; // Status of MD lookup 683 psString secValue = psMetadataLookupString(&mdok, header, section); 684 if (! mdok || ! secValue) { 685 psError(PS_ERR_IO, false, "Unable to locate header %s\n", section); 686 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 687 } else { 688 *trimsec = psRegionFromString(secValue); 689 } 690 } 691 } else { 692 psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE (%s) is not HEADER or VALUE --- trying " 693 "VALUE.\n", source); 694 *trimsec = psRegionFromString(section); 695 } // Value of CELL.TRIMSEC.SOURCE 696 } // Looking up CELL.TRIMSEC.SOURCE 697 } // Looking up CELL.TRIMSEC 698 699 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE, 700 "Trim section", trimsec); 701 psFree(trimsec); 702 703 } else if (strcmp(concept, "CELL.BIASSEC") == 0) { 704 psList *biassecs = psListAlloc(NULL); // List of bias sections 705 706 psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC"); 707 if (! secItem) { 708 psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.\n"); 709 } else if (secItem->type != PS_DATA_STRING) { 710 psError(PS_ERR_IO, true, "CELL.BIASSEC is not of type STR (%x)\n", secItem->type); 711 } else { 712 psString sections = secItem->data.V; // The section string 713 714 psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC.SOURCE"); 715 if (! sourceItem) { 716 psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.SOURCE.\n"); 717 } else if (sourceItem->type != PS_DATA_STRING) { 718 psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); 719 } else { 720 psString source = sourceItem->data.V; // The source string 721 722 psList *secList = papSplit(sections, " ;"); // List of sections 723 psListIterator *secIter = psListIteratorAlloc(secList, PS_LIST_HEAD, false); // Iterator over 724 // sections 725 psString aSection = NULL; // A section from the list 726 while (aSection = psListGetAndIncrement(secIter)) { 727 psRegion *region = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed 728 // by value) 729 730 if (strcasecmp(source, "VALUE") == 0) { 731 *region = psRegionFromString(aSection); 732 } else if (strcasecmp(source, "HEADER") == 0) { 733 psMetadata *header = NULL; // The FITS header 734 if (cell->hdu) { 735 header = cell->hdu->header; 736 } else if (chip->hdu) { 737 header = chip->hdu->header; 738 } else if (fpa->hdu) { 739 header = fpa->hdu->header; 740 } 741 if (! header) { 742 psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); 743 *region = psRegionSet(0.0,0.0,0.0,0.0); 744 } else { 745 bool mdok = true; // Status of MD lookup 746 psString secValue = psMetadataLookupString(&mdok, header, aSection); 747 if (! mdok || ! secValue) { 748 psError(PS_ERR_IO, false, "Unable to locate header %s\n", aSection); 749 *region = psRegionSet(0.0,0.0,0.0,0.0); 750 } else { 751 *region = psRegionFromString(secValue); 752 } 753 } 754 } else { 755 psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE (%s) is not HEADER or VALUE --- trying " 756 "VALUE.\n", source); 757 *region = psRegionFromString(aSection); 758 } // Value of CELL.BIASSEC.SOURCE 759 760 psListAdd(biassecs, PS_LIST_TAIL, region); 761 psFree(region); 762 } // Iterating over multiple sections 763 psFree(secIter); 764 psFree(secList); 765 } // Looking up CELL.BIASSEC.SOURCE 766 } // Looking up CELL.BIASSEC 767 768 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE, 769 "Bias sections", biassecs); 770 psFree(biassecs); 771 772 } else if (strcmp(concept, "CELL.XBIN") == 0) { 773 int xBin = 1; // Binning factor in x 774 psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.XBIN"); 775 if (! binItem) { 776 psError(PS_ERR_IO, false, "Couldn't find CELL.XBIN.\n"); 777 } else if (binItem->type == PS_DATA_STRING) { 778 psString binString = binItem->data.V; // The string containing the binning 779 if (sscanf(binString, "%d %*d", &xBin) != 1 && 780 sscanf(binString, "%d,%*d", &xBin) != 1) { 781 psError(PS_ERR_IO, true, "Unable to read string to get x binning: %s\n", binString); 782 } 783 } else if (binItem->type == PS_TYPE_S32) { 784 xBin = binItem->data.S32; 785 } else { 786 psError(PS_ERR_IO, true, "Note sure how to interpret CELL.XBIN of type %x --- assuming 1.\n", 787 binItem->type); 788 } 789 790 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_TYPE_S32 | PS_META_REPLACE, 791 "Binning in x", xBin); 792 793 } else if (strcmp(concept, "CELL.YBIN") == 0) { 794 int yBin = 1; // Binning factor in y 795 psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.YBIN"); 796 if (! binItem) { 797 psError(PS_ERR_IO, false, "Couldn't find CELL.YBIN.\n"); 798 } else if (binItem->type == PS_DATA_STRING) { 799 psString binString = binItem->data.V; // The string containing the binning 800 if (sscanf(binString, "%*d %d", &yBin) != 1 && 801 sscanf(binString, "%*d,%d", &yBin) != 1) { 802 psError(PS_ERR_IO, true, "Unable to read string to get y binning: %s\n", binString); 803 } 804 } else if (binItem->type == PS_TYPE_S32) { 805 yBin = binItem->data.S32; 806 } else { 807 psError(PS_ERR_IO, true, "Note sure how to interpret CELL.YBIN of type %x --- assuming 1.\n", 808 binItem->type); 809 } 810 811 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_TYPE_S32 | PS_META_REPLACE, 812 "Binning in y", yBin); 813 814 } else if (strcmp(concept, "CELL.TIME") == 0 || strcmp(concept, "CELL.TIMESYS") == 0) { 815 // Do CELL.TIME and CELL.TIMESYS together 816 psTime *time = NULL; // The time 817 psTimeType timeSys = PS_TIME_UTC; // The time system 818 819 // CELL.TIMESYS 820 psMetadataItem *sysItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIMESYS"); 821 if (! sysItem) { 822 psError(PS_ERR_IO, true, "Couldn't find CELL.TIMESYS --- assuming UTC.\n"); 823 } else if (sysItem->type != PS_DATA_STRING) { 824 psError(PS_ERR_IO, true, "CELL.TIMESYS isn't of type STRING --- assuming UTC.\n"); 825 } else { 826 psString sys = sysItem->data.V; // The time system string 827 if (strcasecmp(sys, "TAI") == 0) { 828 timeSys = PS_TIME_TAI; 829 } else if (strcasecmp(sys, "UTC") == 0) { 830 timeSys = PS_TIME_UTC; 831 } else if (strcasecmp(sys, "UT1") == 0) { 832 timeSys = PS_TIME_UT1; 833 } else if (strcasecmp(sys, "TT") == 0) { 834 timeSys = PS_TIME_TT; 835 } else { 836 psError(PS_ERR_IO, true, "Can't interpret CELL.TIMESYS --- assuming UTC.\n"); 837 } 838 } 839 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIMESYS", PS_TYPE_S32 | PS_META_REPLACE, 840 "Time system", timeSys); 841 842 psMetadataItem *timeItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIME"); 843 if (! timeItem) { 844 psError(PS_ERR_IO, false, "Couldn't find CELL.TIME.\n"); 845 } else { 846 // Get format 847 const psMetadata *camera = fpa->camera; // The camera configuration data 848 bool mdok = true; // Status of MD lookup 849 psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); 850 if (mdok && formats) { 851 psString timeFormat = psMetadataLookupString(&mdok, formats, "CELL.TIME"); 852 if (mdok && strlen(timeFormat) > 0) { 853 switch (timeItem->type) { 854 case PS_DATA_STRING: 855 { 856 psString timeString = timeItem->data.V; // String with the time 857 if (strcasecmp(timeFormat, "ISO") == 0) { 858 // timeString contains an ISO time 859 time = psTimeFromISO(timeString, timeSys); 860 } else if (strstr(timeFormat, "SEPARATE")) { 861 // timeString contains headers for the date and time 862 psMetadata *header = NULL; // The FITS header 863 if (cell->hdu) { 864 header = cell->hdu->header; 865 } else if (chip->hdu) { 866 header = chip->hdu->header; 867 } else if (fpa->hdu) { 868 header = fpa->hdu->header; 869 } 870 if (! header) { 871 psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); 872 } else { 873 // Get the headers 874 char *stuff1 = strpbrk(timeString, " ,;"); 875 psString dateName = psStringNCopy(timeString, 876 strlen(timeString) - strlen(stuff1)); 877 char *stuff2 = strpbrk(stuff1, "abcdefghijklmnopqrstuvwxyz" 878 "ABCDEFGHIJKLMNOPQRSTUVWXYZ"); 879 psString timeName = psStringCopy(stuff2); 880 881 bool mdok = true; // Status of MD lookup 882 psString dateString = psMetadataLookupString(&mdok, header, dateName); 883 psFree(dateName); 884 int day = 0, month = 0, year = 0; 885 if (sscanf(dateString, "%d-%d-%d", &day, &month, &year) != 3 && 886 sscanf(dateString, "%d/%d/%d", &day, &month, &year) != 3) { 887 psError(PS_ERR_IO, true, "Unable to read date: %s\n", dateString); 888 } else { 889 if (strstr(timeFormat, "BACKWARDS")) { 890 int temp = day; 891 day = year; 892 year = temp; 893 } 894 if (strstr(timeFormat, "PRE2000") || year < 2000) { 895 year += 2000; 896 } 897 898 psMetadataItem *timeItem = psMetadataLookup(header, timeName); 899 if (! timeItem) { 900 psError(PS_ERR_IO, false, "Unable to find time header: %s\n", 901 timeName); 902 } else if (timeItem->type == PS_DATA_STRING) { 903 // Time is a string, in the usual way: 904 psStringAppend(&dateString, "T%s", timeItem->data.V); 905 } else { 906 // Assume that time is specified in Second of Day 907 double seconds = NAN; 908 switch (timeItem->type) { 909 case PS_TYPE_S32: 910 seconds = timeItem->data.S32; 911 break; 912 case PS_TYPE_F32: 913 seconds = timeItem->data.F32; 914 break; 915 case PS_TYPE_F64: 916 seconds = timeItem->data.F64; 917 break; 918 default: 919 psError(PS_ERR_IO, true, "Time header (%s) is not of an " 920 "expected type: %x\n", timeName, timeItem->type); 921 } 922 // Now print to timeString as "hh:mm:ss.ss" 923 int hours = seconds / 3600; 924 seconds -= (double)hours * 3600.0; 925 int minutes = seconds / 60; 926 seconds -= (double)minutes * 60.0; 927 psStringAppend(&dateString, "T%02d:%02d:%02f", hours, minutes, 928 seconds); 929 } 930 time = psTimeFromISO(dateString, timeSys); 931 } // Reading date and time 932 psFree(timeName); 933 } // Reading headers 934 } else { 935 psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%s) --- trying " 936 "ISO\n", timeString); 937 time = psTimeFromISO(timeString, timeSys); 938 } // Interpreting the time string 939 } 940 break; 941 case PS_TYPE_F32: 942 { 943 double timeValue = (double)timeItem->data.F32; 944 if (strcasecmp(timeFormat, "JD") == 0) { 945 time = psTimeFromJD(timeValue); 946 } else if (strcasecmp(timeFormat, "MJD") == 0) { 947 time = psTimeFromMJD(timeValue); 948 } else { 949 psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " 950 "JD\n", timeValue); 951 time = psTimeFromJD(timeValue); 952 } 953 } 954 break; 955 case PS_TYPE_F64: 956 { 957 double timeValue = (double)timeItem->data.F64; 958 if (strcasecmp(timeFormat, "JD") == 0) { 959 time = psTimeFromJD(timeValue); 960 } else if (strcasecmp(timeFormat, "MJD") == 0) { 961 time = psTimeFromMJD(timeValue); 962 } else { 963 psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " 964 "JD\n", timeValue); 965 time = psTimeFromJD(timeValue); 966 } 967 } 968 break; 969 default: 970 psError(PS_ERR_IO, true, "Unable to parse CELL.TIME.\n"); 971 } 972 } else { 973 psError(PS_ERR_IO, false, "Unable to find CELL.TIME in FORMATS.\n"); 974 } // Getting the format 975 } else { 976 psError(PS_ERR_IO, false, "Unable to find FORMATS in camera configuration.\n"); 977 } // Getting the formats 978 } // Getting CELL.TIME 979 980 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIME", PS_DATA_TIME | PS_META_REPLACE, 981 "Time of exposure", time); 982 psFree(time); 983 984 // These are new and experimental concepts: CELL.X0 and CELL.Y0 985 } else if (strcmp(concept, "CELL.X0") == 0) { 986 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.X0", "Position of (0,0) on the chip"); 987 correctPosition(fpa, cell, "CELL.X0"); 988 } else if (strcmp(concept, "CELL.Y0") == 0) { 989 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.Y0", "Position of (0,0) on the chip"); 990 correctPosition(fpa, cell, "CELL.Y0"); 991 } 992 993 return true; 994 } 995 593 996 594 997 595 998 // Ingest concepts for the cell 596 999 bool pmCellIngestConcepts(pmCell *cell, // The cell 597 psDB *db// DB handle598 ) 599 { 600 pmChip *chip = cell->parent; // The parent chip601 pmFPA *fpa = chip->parent; // The parent FPA1000 psDB *db // DB handle 1001 ) 1002 { 1003 pmChip *chip = cell->parent; // The parent chip 1004 pmFPA *fpa = chip->parent; // The parent FPA 602 1005 603 1006 if (! cell->concepts) { 604 cell->concepts = psMetadataAlloc();1007 cell->concepts = psMetadataAlloc(); 605 1008 } 606 1009 607 1010 // CELL.NAME --- added by pmFPAConstruct 608 1011 609 // CELL.GAIN 610 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.GAIN", "CCD gain (e/count)"); 611 612 // CELL.READNOISE 613 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.READNOISE", "CCD read noise (e)"); 614 615 // CELL.SATURATION 616 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.SATURATION", "Saturation level (ADU)"); 617 618 // CELL.BAD 619 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.BAD", "Bad level (ADU)"); 620 621 // CELL.XPARITY 622 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.XPARITY", "Orientation in x compared to the " 623 "rest of the FPA"); 624 625 // CELL.YPARITY 626 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.YPARITY", "Orientation in y compared to the " 627 "rest of the FPA"); 628 629 // CELL.READDIR 630 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.READDIR", "Read direction: 1=row, 2=col"); 1012 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.GAIN"); 1013 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.READNOISE"); 1014 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.SATURATION"); 1015 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.BAD"); 1016 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.XPARITY"); 1017 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.YPARITY"); 1018 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.READDIR"); 631 1019 632 1020 // These used to be pmReadoutGetExposure and pmReadoutGetDarkTime, but that doesn't really make sense at … … 634 1022 // REALLY derived? They're not in the FITS headers, because a readout is a plane in a 3D image. We'll 635 1023 // have to dream up some additional suffix to specify these, but for now.... 636 637 // CELL.EXPOSURE (used to be READOUT.EXPOSURE) 638 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.EXPOSURE", "Exposure time (sec)"); 639 640 // CELL.DARKTIME (used to be READOUT.DARKTIME) 641 p_pmFPAConceptGetF32(fpa, chip, cell, db, cell->concepts, "CELL.DARKTIME", "Time since CCD flush (sec)"); 642 643 // These take some extra work 644 645 // CELL.TRIMSEC 646 { 647 psRegion *trimsec = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed by value) 648 649 psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC"); 650 if (! secItem) { 651 psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.\n"); 652 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 653 } else if (secItem->type != PS_DATA_STRING) { 654 psError(PS_ERR_IO, true, "CELL.TRIMSEC is not of type STR (%x)\n", secItem->type); 655 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 656 } else { 657 psString section = secItem->data.V; // The section string 658 659 psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TRIMSEC.SOURCE"); 660 if (! sourceItem) { 661 psError(PS_ERR_IO, false, "Couldn't find CELL.TRIMSEC.SOURCE.\n"); 662 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 663 } else if (sourceItem->type != PS_DATA_STRING) { 664 psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); 665 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 666 } else { 667 psString source = sourceItem->data.V; // The source string 668 669 if (strcasecmp(source, "VALUE") == 0) { 670 *trimsec = psRegionFromString(section); 671 } else if (strcasecmp(source, "HEADER") == 0) { 672 psMetadata *header = NULL; // The FITS header 673 if (cell->hdu) { 674 header = cell->hdu->header; 675 } else if (chip->hdu) { 676 header = chip->hdu->header; 677 } else if (fpa->hdu) { 678 header = fpa->hdu->header; 679 } 680 if (! header) { 681 psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); 682 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 683 } else { 684 bool mdok = true; // Status of MD lookup 685 psString secValue = psMetadataLookupString(&mdok, header, section); 686 if (! mdok || ! secValue) { 687 psError(PS_ERR_IO, false, "Unable to locate header %s\n", section); 688 *trimsec = psRegionSet(0.0, 0.0, 0.0, 0.0); 689 } else { 690 *trimsec = psRegionFromString(secValue); 691 } 692 } 693 } else { 694 psError(PS_ERR_IO, true, "CELL.TRIMSEC.SOURCE (%s) is not HEADER or VALUE --- trying " 695 "VALUE.\n", source); 696 *trimsec = psRegionFromString(section); 697 } // Value of CELL.TRIMSEC.SOURCE 698 } // Looking up CELL.TRIMSEC.SOURCE 699 } // Looking up CELL.TRIMSEC 700 701 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN, 702 "Trim section", trimsec); 703 psFree(trimsec); 704 } 705 706 // CELL.BIASSEC 707 { 708 psList *biassecs = psListAlloc(NULL); // List of bias sections 709 710 psMetadataItem *secItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC"); 711 if (! secItem) { 712 psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.\n"); 713 } else if (secItem->type != PS_DATA_STRING) { 714 psError(PS_ERR_IO, true, "CELL.BIASSEC is not of type STR (%x)\n", secItem->type); 715 } else { 716 psString sections = secItem->data.V; // The section string 717 718 psMetadataItem *sourceItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.BIASSEC.SOURCE"); 719 if (! sourceItem) { 720 psError(PS_ERR_IO, false, "Couldn't find CELL.BIASSEC.SOURCE.\n"); 721 } else if (sourceItem->type != PS_DATA_STRING) { 722 psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE is not of type STR (%x)\n", sourceItem->type); 723 } else { 724 psString source = sourceItem->data.V; // The source string 725 726 psList *secList = papSplit(sections, " ;"); // List of sections 727 psListIterator *secIter = psListIteratorAlloc(secList, PS_LIST_HEAD, false); // Iterator over 728 // sections 729 psString aSection = NULL; // A section from the list 730 while (aSection = psListGetAndIncrement(secIter)) { 731 psRegion *region = psAlloc(sizeof(psRegion)); // Make space for a psRegion (usually passed 732 // by value) 733 734 if (strcasecmp(source, "VALUE") == 0) { 735 *region = psRegionFromString(aSection); 736 } else if (strcasecmp(source, "HEADER") == 0) { 737 psMetadata *header = NULL; // The FITS header 738 if (cell->hdu) { 739 header = cell->hdu->header; 740 } else if (chip->hdu) { 741 header = chip->hdu->header; 742 } else if (fpa->hdu) { 743 header = fpa->hdu->header; 744 } 745 if (! header) { 746 psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); 747 *region = psRegionSet(0.0,0.0,0.0,0.0); 748 } else { 749 bool mdok = true; // Status of MD lookup 750 psString secValue = psMetadataLookupString(&mdok, header, aSection); 751 if (! mdok || ! secValue) { 752 psError(PS_ERR_IO, false, "Unable to locate header %s\n", aSection); 753 *region = psRegionSet(0.0,0.0,0.0,0.0); 754 } else { 755 *region = psRegionFromString(secValue); 756 } 757 } 758 } else { 759 psError(PS_ERR_IO, true, "CELL.BIASSEC.SOURCE (%s) is not HEADER or VALUE --- trying " 760 "VALUE.\n", source); 761 *region = psRegionFromString(aSection); 762 } // Value of CELL.BIASSEC.SOURCE 763 764 psListAdd(biassecs, PS_LIST_TAIL, region); 765 psFree(region); 766 } // Iterating over multiple sections 767 psFree(secIter); 768 psFree(secList); 769 } // Looking up CELL.BIASSEC.SOURCE 770 } // Looking up CELL.BIASSEC 771 772 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST, "Bias sections", biassecs); 773 psFree(biassecs); 774 } 775 776 // CELL.XBIN 777 { 778 int xBin = 1; // Binning factor in x 779 psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.XBIN"); 780 if (! binItem) { 781 psError(PS_ERR_IO, false, "Couldn't find CELL.XBIN.\n"); 782 } else if (binItem->type == PS_DATA_STRING) { 783 psString binString = binItem->data.V; // The string containing the binning 784 if (sscanf(binString, "%d %*d", &xBin) != 1 && 785 sscanf(binString, "%d,%*d", &xBin) != 1) { 786 psError(PS_ERR_IO, true, "Unable to read string to get x binning: %s\n", binString); 787 } 788 } else if (binItem->type == PS_TYPE_S32) { 789 xBin = binItem->data.S32; 790 } else { 791 psError(PS_ERR_IO, true, "Note sure how to interpret CELL.XBIN of type %x --- assuming 1.\n", 792 binItem->type); 793 } 794 795 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_TYPE_S32, "Binning in x", xBin); 796 } 797 798 // CELL.XBIN 799 { 800 int yBin = 1; // Binning factor in y 801 psMetadataItem *binItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.YBIN"); 802 if (! binItem) { 803 psError(PS_ERR_IO, false, "Couldn't find CELL.YBIN.\n"); 804 } else if (binItem->type == PS_DATA_STRING) { 805 psString binString = binItem->data.V; // The string containing the binning 806 if (sscanf(binString, "%*d %d", &yBin) != 1 && 807 sscanf(binString, "%*d,%d", &yBin) != 1) { 808 psError(PS_ERR_IO, true, "Unable to read string to get y binning: %s\n", binString); 809 } 810 } else if (binItem->type == PS_TYPE_S32) { 811 yBin = binItem->data.S32; 812 } else { 813 psError(PS_ERR_IO, true, "Note sure how to interpret CELL.YBIN of type %x --- assuming 1.\n", 814 binItem->type); 815 } 816 817 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_TYPE_S32, "Binning in y", yBin); 818 } 819 820 // CELL.TIME and CELL.TIMESYS 821 { 822 psTime *time = NULL; // The time 823 psTimeType timeSys = PS_TIME_UTC; // The time system 824 825 // CELL.TIMESYS 826 psMetadataItem *sysItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIMESYS"); 827 if (! sysItem) { 828 psError(PS_ERR_IO, true, "Couldn't find CELL.TIMESYS --- assuming UTC.\n"); 829 } else if (sysItem->type != PS_DATA_STRING) { 830 psError(PS_ERR_IO, true, "CELL.TIMESYS isn't of type STRING --- assuming UTC.\n"); 831 } else { 832 psString sys = sysItem->data.V; // The time system string 833 if (strcasecmp(sys, "TAI") == 0) { 834 timeSys = PS_TIME_TAI; 835 } else if (strcasecmp(sys, "UTC") == 0) { 836 timeSys = PS_TIME_UTC; 837 } else if (strcasecmp(sys, "UT1") == 0) { 838 timeSys = PS_TIME_UT1; 839 } else if (strcasecmp(sys, "TT") == 0) { 840 timeSys = PS_TIME_TT; 841 } else { 842 psError(PS_ERR_IO, true, "Can't interpret CELL.TIMESYS --- assuming UTC.\n"); 843 } 844 } 845 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIMESYS", PS_TYPE_S32, "Time system", timeSys); 846 847 psMetadataItem *timeItem = p_pmFPAConceptGet(fpa, chip, cell, db, "CELL.TIME"); 848 if (! timeItem) { 849 psError(PS_ERR_IO, false, "Couldn't find CELL.TIME.\n"); 850 } else { 851 // Get format 852 const psMetadata *camera = fpa->camera; // The camera configuration data 853 bool mdok = true; // Status of MD lookup 854 psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); 855 if (mdok && formats) { 856 psString timeFormat = psMetadataLookupString(&mdok, formats, "CELL.TIME"); 857 if (mdok && strlen(timeFormat) > 0) { 858 switch (timeItem->type) { 859 case PS_DATA_STRING: 860 { 861 psString timeString = timeItem->data.V; // String with the time 862 if (strcasecmp(timeFormat, "ISO") == 0) { 863 // timeString contains an ISO time 864 time = psTimeFromISO(timeString, timeSys); 865 } else if (strstr(timeFormat, "SEPARATE")) { 866 // timeString contains headers for the date and time 867 psMetadata *header = NULL; // The FITS header 868 if (cell->hdu) { 869 header = cell->hdu->header; 870 } else if (chip->hdu) { 871 header = chip->hdu->header; 872 } else if (fpa->hdu) { 873 header = fpa->hdu->header; 874 } 875 if (! header) { 876 psError(PS_ERR_IO, true, "Unable to find FITS header!\n"); 877 } else { 878 // Get the headers 879 char *stuff1 = strpbrk(timeString, " ,;"); 880 psString dateName = psStringNCopy(timeString, 881 strlen(timeString) - strlen(stuff1)); 882 char *stuff2 = strpbrk(stuff1, 883 "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"); 884 psString timeName = psStringCopy(stuff2); 885 886 bool mdok = true; // Status of MD lookup 887 psString dateString = psMetadataLookupString(&mdok, header, dateName); 888 psFree(dateName); 889 int day = 0, month = 0, year = 0; 890 if (sscanf(dateString, "%d-%d-%d", &day, &month, &year) != 3 && 891 sscanf(dateString, "%d/%d/%d", &day, &month, &year) != 3) { 892 psError(PS_ERR_IO, true, "Unable to read date: %s\n", dateString); 893 } else { 894 if (strstr(timeFormat, "BACKWARDS")) { 895 int temp = day; 896 day = year; 897 year = temp; 898 } 899 if (strstr(timeFormat, "PRE2000") || year < 2000) { 900 year += 2000; 901 } 902 903 psMetadataItem *timeItem = psMetadataLookup(header, timeName); 904 if (! timeItem) { 905 psError(PS_ERR_IO, false, "Unable to find time header: %s\n", 906 timeName); 907 } else if (timeItem->type == PS_DATA_STRING) { 908 // Time is a string, in the usual way: 909 psStringAppend(&dateString, "T%s", timeItem->data.V); 910 } else { 911 // Assume that time is specified in Second of Day 912 double seconds = NAN; 913 switch (timeItem->type) { 914 case PS_TYPE_S32: 915 seconds = timeItem->data.S32; 916 break; 917 case PS_TYPE_F32: 918 seconds = timeItem->data.F32; 919 break; 920 case PS_TYPE_F64: 921 seconds = timeItem->data.F64; 922 break; 923 default: 924 psError(PS_ERR_IO, true, "Time header (%s) is not of an " 925 "expected type: %x\n", timeName, timeItem->type); 926 } 927 // Now print to timeString as "hh:mm:ss.ss" 928 int hours = seconds / 3600; 929 seconds -= (double)hours * 3600.0; 930 int minutes = seconds / 60; 931 seconds -= (double)minutes * 60.0; 932 psStringAppend(&dateString, "T%02d:%02d:%02f", hours, minutes, 933 seconds); 934 } 935 time = psTimeFromISO(dateString, timeSys); 936 } // Reading date and time 937 psFree(timeName); 938 } // Reading headers 939 } else { 940 psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%s) --- trying " 941 "ISO\n", timeString); 942 time = psTimeFromISO(timeString, timeSys); 943 } // Interpreting the time string 944 } 945 break; 946 case PS_TYPE_F32: 947 { 948 double timeValue = (double)timeItem->data.F32; 949 if (strcasecmp(timeFormat, "JD") == 0) { 950 time = psTimeFromJD(timeValue); 951 } else if (strcasecmp(timeFormat, "MJD") == 0) { 952 time = psTimeFromMJD(timeValue); 953 } else { 954 psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " 955 "JD\n", timeValue); 956 time = psTimeFromJD(timeValue); 957 } 958 } 959 break; 960 case PS_TYPE_F64: 961 { 962 double timeValue = (double)timeItem->data.F64; 963 if (strcasecmp(timeFormat, "JD") == 0) { 964 time = psTimeFromJD(timeValue); 965 } else if (strcasecmp(timeFormat, "MJD") == 0) { 966 time = psTimeFromMJD(timeValue); 967 } else { 968 psError(PS_ERR_IO, true, "Not sure how to parse CELL.TIME (%f) --- trying " 969 "JD\n", timeValue); 970 time = psTimeFromJD(timeValue); 971 } 972 } 973 break; 974 default: 975 psError(PS_ERR_IO, true, "Unable to parse CELL.TIME.\n"); 976 } 977 } else { 978 psError(PS_ERR_IO, false, "Unable to find CELL.TIME in FORMATS.\n"); 979 } // Getting the format 980 } else { 981 psError(PS_ERR_IO, false, "Unable to find FORMATS in camera configuration.\n"); 982 } // Getting the formats 983 } // Getting CELL.TIME 984 985 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.TIME", PS_DATA_UNKNOWN, "Time of exposure", time); 986 psFree(time); 987 } 988 989 // These are new and experimental concepts: CELL.X0 and CELL.Y0 990 991 // CELL.X0 992 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.X0", "Position of (0,0) on the chip "); 993 // CELL.Y0 994 p_pmFPAConceptGetS32(fpa, chip, cell, db, cell->concepts, "CELL.Y0", "Position of (0,0) on the chip"); 995 // Add corrective 996 { 997 const psMetadata *camera = fpa->camera; 998 bool mdok = false; // Result of MD lookup 999 psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); 1000 if (mdok && formats) { 1001 psString format = psMetadataLookupString(&mdok, formats, "CELL.X0"); 1002 if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { 1003 psMetadataItem *cellx0 = psMetadataLookup(cell->concepts, "CELL.X0"); 1004 cellx0->data.S32 -= 1; 1005 } 1006 format = psMetadataLookupString(&mdok, formats, "CELL.Y0"); 1007 if (mdok && strlen(format) > 0 && strcasecmp(format, "FORTRAN") == 0) { 1008 psMetadataItem *celly0 = psMetadataLookup(cell->concepts, "CELL.Y0"); 1009 celly0->data.S32 -= 1; 1010 } 1011 } 1012 } 1013 1014 // Pau. 1024 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.EXPOSURE"); 1025 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.DARKTIME"); 1026 1027 1028 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.TRIMSEC"); 1029 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.BIASSEC"); 1030 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.TIME"); // This also does CELL.TIMESYS 1031 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.XBIN"); 1032 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.YBIN"); 1033 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.X0"); 1034 p_pmCellIngestConcept(fpa, chip, cell, db, "CELL.Y0"); 1035 1015 1036 } 1016 1037 … … 1018 1039 // Retrieve a concept 1019 1040 psMetadataItem *pmCellGetConcept(pmCell *cell, // The cell 1020 const char *concept // The concept1041 const char *concept // The concept 1021 1042 ) 1022 1043 { 1023 1044 psMetadataItem *item = psMetadataLookup(cell->concepts, concept); 1024 1045 if (! item) { 1025 pmChip *chip = cell->parent;1026 item = psMetadataLookup(chip->concepts, concept);1027 if (! item) {1028 pmFPA *fpa = chip->parent;1029 item = psMetadataLookup(fpa->concepts, concept);1030 }1031 } 1032 return item; // item is either NULL or is what we want.1033 } 1046 pmChip *chip = cell->parent; 1047 item = psMetadataLookup(chip->concepts, concept); 1048 if (! item) { 1049 pmFPA *fpa = chip->parent; 1050 item = psMetadataLookup(fpa->concepts, concept); 1051 } 1052 } 1053 return item; // item is either NULL or is what we want. 1054 } -
trunk/archive/scripts/src/phase2/pmFPAConceptsGet.h
r5105 r5786 7 7 8 8 psMetadataItem *p_pmFPAConceptGetFromCamera(pmCell *cell, // The cell 9 const char *concept // Name of concept9 const char *concept // Name of concept 10 10 ); 11 11 psMetadataItem *p_pmFPAConceptGetFromHeader(pmFPA *fpa, // The FPA that contains the chip 12 pmChip *chip, // The chip that contains the cell13 pmCell *cell, // The cell14 const char *concept // Name of concept12 pmChip *chip, // The chip that contains the cell 13 pmCell *cell, // The cell 14 const char *concept // Name of concept 15 15 ); 16 16 psMetadataItem *p_pmFPAConceptGetFromDefault(pmFPA *fpa, // The FPA that contains the chip 17 pmChip *chip, // The chip that contains the cell18 pmCell *cell, // The cell19 const char *concept // Name of concept17 pmChip *chip, // The chip that contains the cell 18 pmCell *cell, // The cell 19 const char *concept // Name of concept 20 20 ); 21 21 psMetadataItem *p_pmFPAConceptGetFromDB(pmFPA *fpa, // The FPA that contains the chip 22 pmChip *chip, // The chip that contains the cell23 pmCell *cell, // The cell24 psDB *db,// DB handle25 const char *concept // Name of concept22 pmChip *chip, // The chip that contains the cell 23 pmCell *cell, // The cell 24 psDB *db, // DB handle 25 const char *concept // Name of concept 26 26 ); 27 27 psMetadataItem *p_pmFPAConceptGet(pmFPA *fpa, // The FPA 28 pmChip *chip,// The chip29 pmCell *cell,// The cell30 psDB *db, // DB handle31 const char *concept // Concept name28 pmChip *chip,// The chip 29 pmCell *cell, // The cell 30 psDB *db, // DB handle 31 const char *concept // Concept name 32 32 ); 33 void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA34 pmChip *chip, // The chip35 pmCell *cell, // The cell36 psDB *db,// DB handle37 psMetadata *concepts, // The concepts MD38 const char *name, // Name of the concept39 const char *comment // Comment for concept33 void p_pmFPAConceptGetF32(pmFPA *fpa, // The FPA 34 pmChip *chip, // The chip 35 pmCell *cell, // The cell 36 psDB *db, // DB handle 37 psMetadata *concepts, // The concepts MD 38 const char *name, // Name of the concept 39 const char *comment // Comment for concept 40 40 ); 41 void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA42 pmChip *chip, // The chip43 pmCell *cell, // The cell44 psDB *db,// DB handle45 psMetadata *concepts, // The concepts MD46 const char *name, // Name of the concept47 const char *comment // Comment for concept41 void p_pmFPAConceptGetF64(pmFPA *fpa, // The FPA 42 pmChip *chip, // The chip 43 pmCell *cell, // The cell 44 psDB *db, // DB handle 45 psMetadata *concepts, // The concepts MD 46 const char *name, // Name of the concept 47 const char *comment // Comment for concept 48 48 ); 49 void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA50 pmChip *chip,// The chip51 pmCell *cell,// The cell52 psDB *db,// DB handle53 psMetadata *concepts, // The concepts MD54 const char *name, // Name of the concept55 const char *comment // Comment for concept49 void p_pmFPAConceptGetS32(pmFPA *fpa, // The FPA 50 pmChip *chip, // The chip 51 pmCell *cell, // The cell 52 psDB *db, // DB handle 53 psMetadata *concepts, // The concepts MD 54 const char *name, // Name of the concept 55 const char *comment // Comment for concept 56 56 ); 57 57 void p_pmFPAConceptGetString(pmFPA *fpa, // The FPA 58 pmChip *chip, // The chip59 pmCell *cell, // The cell60 psDB *db,// DB handle61 psMetadata *concepts, // The concepts MD62 const char *name, // Name of the concept63 const char *comment // Comment for concept58 pmChip *chip, // The chip 59 pmCell *cell, // The cell 60 psDB *db, // DB handle 61 psMetadata *concepts, // The concepts MD 62 const char *name, // Name of the concept 63 const char *comment // Comment for concept 64 64 ); 65 65 66 66 67 // Ingest a single cell concept; this is the workhorse for ingesting cell concepts 68 bool p_pmCellIngestConcept(pmFPA *fpa, // The FPA 69 pmChip *chip, // The chip 70 pmCell *cell, // The cell 71 psDB *db, // DB handle 72 const char *concept // Name of the concept 73 ); 74 67 75 // Ingest concepts for the FPA 68 void pmFPAIngestConcepts(pmFPA *fpa, // The FPA69 psDB *db// DB handle76 void pmFPAIngestConcepts(pmFPA *fpa, // The FPA 77 psDB *db // DB handle 70 78 ); 71 79 // Ingest concepts for the chip 72 bool pmChipIngestConcepts(pmChip *chip, // The chip73 psDB *db// DB handle80 bool pmChipIngestConcepts(pmChip *chip, // The chip 81 psDB *db // DB handle 74 82 ); 75 83 // Ingest concepts for the cell 76 84 bool pmCellIngestConcepts(pmCell *cell, // The cell 77 psDB *db// DB handle85 psDB *db // DB handle 78 86 ); 79 87 80 88 // Retrieve a concept 81 89 psMetadataItem *pmCellGetConcept(pmCell *cell, // The cell 82 const char *concept // The concept90 const char *concept // The concept 83 91 ); 84 92 -
trunk/archive/scripts/src/phase2/pmFPAConceptsSet.c
r5633 r5786 109 109 const char *keyword = psMetadataLookupString(&mdok, translation, concept->name); 110 110 if (mdok && strlen(keyword) > 0) { 111 psTrace(__func__, 6, "It's in keyword %s\n", keyword); 111 112 psMetadataItem *headerItem = NULL; // Item to add to header 112 113 // XXX: Need to expand range of types … … 131 132 // We have a FITS header to look up --- search each level 132 133 if (cell && cell->hdu) { 133 psMetadataItem *cellItem = psMetadataLookup(cell->hdu->header, keyword); 134 if (cellItem) { 135 // XXX: Need to clean up before returning 136 psMetadataAddItem(cell->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 137 status = true; 138 } 134 psTrace(__func__, 7, "Adding to the cell level header...\n"); 135 psMetadataAddItem(cell->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 136 status = true; 139 137 } else if (chip && chip->hdu) { 140 psMetadataItem *chipItem = psMetadataLookup(chip->hdu->header, keyword); 141 if (chipItem) { 142 // XXX: Need to clean up before returning 143 psMetadataAddItem(chip->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 144 status = true; 145 } 138 psTrace(__func__, 7, "Adding to the chip level header...\n"); 139 psMetadataAddItem(chip->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 140 status = true; 146 141 } else if (fpa->hdu) { 147 psMetadataItem *fpaItem = psMetadataLookup(fpa->hdu->header, keyword); 148 if (fpaItem) { 149 // XXX: Need to clean up before returning 150 psMetadataAddItem(fpa->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 151 status = true; 152 } 153 } else if (fpa->phu) { 154 psMetadataItem *fpaItem = psMetadataLookup(fpa->phu, keyword); 155 if (fpaItem) { 156 // XXX: Need to clean up before returning 157 psMetadataAddItem(fpa->phu, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 158 status = true; 159 } 142 psTrace(__func__, 7, "Adding to the FPA level header...\n"); 143 psMetadataAddItem(fpa->hdu->header, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 144 status = true; 145 } else { 146 // In desperation, add to the PHU --- it HAS to be in the header somewhere! 147 if (! fpa->phu) { 148 fpa->phu = psMetadataAlloc(); 149 } 150 psTrace(__func__, 7, "Adding to the PHU...\n"); 151 psMetadataAddItem(fpa->phu, headerItem, PS_LIST_TAIL, PS_META_REPLACE); 152 status = true; 160 153 } 161 154 psFree(headerItem); … … 182 175 183 176 psMetadataItem *defItem = psMetadataLookup(defaults, concept->name); 177 bool status = false; // Result of checking the database 184 178 if (defItem) { 185 179 if (defItem->type == PS_DATA_METADATA) { … … 209 203 " --- ignored.\n", dependsOn); 210 204 } 211 212 defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after 205 206 defItem = psMetadataLookup(dependents, depItem->data.V); // This is now what we were after 207 } 208 209 status = compareConcepts(defItem, concept); 210 if (! status) { 211 psError(PS_ERR_IO, true, "Concept %s is specified by default in the camera configuration, " 212 "but doesn't match the actual value.\n", concept->name); 213 213 } 214 214 } 215 215 216 216 // XXX: Need to clean up before returning 217 return compareConcepts(defItem, concept);217 return status; 218 218 } 219 219 … … 330 330 { 331 331 // Try headers, database, defaults in order 332 psTrace(__func__, 3, "Trying to set concept %s...\n", concept->name); 332 333 bool status = setConceptInCamera(cell, concept); // Status for return 333 334 if (! status) { 335 psTrace(__func__, 5, "Trying header....\n"); 334 336 status = setConceptInHeader(fpa, chip, cell, concept); 335 337 } 336 338 if (! status) { 339 psTrace(__func__, 5, "Trying database....\n"); 337 340 status = setConceptInDB(fpa, chip, cell, db, concept); 338 341 } 339 342 if (! status) { 343 psTrace(__func__, 5, "Checking defaults....\n"); 340 344 status = setConceptInDefault(fpa, chip, cell, concept); 341 345 } -
trunk/archive/scripts/src/phase2/pmFPAConstruct.c
r5651 r5786 67 67 const char *phuType = psMetadataLookupString(&mdStatus, camera, "PHU"); // What is the PHU? 68 68 const char *extType = psMetadataLookupString(&mdStatus, camera, "EXTENSIONS"); // What's in the extns? 69 70 69 71 70 if (strcasecmp(phuType, "FPA") == 0) { -
trunk/archive/scripts/src/phase2/pmFPAMorph.c
r5651 r5786 14 14 // Flip the chip in x 15 15 static psImage *xFlipImage(psImage *image, // Image to be flipped 16 int size// Size of image in flip dimension16 int size // Size of image in flip dimension 17 17 ) 18 18 { … … 20 20 switch (image->type.type) { 21 21 case PS_TYPE_U16: 22 for (int y = 0; y < image->numRows; y++) {23 for (int x = 0; x < image->numCols; x++) {24 flipped->data.U16[y][size - x - 1] = image->data.U16[y][x];25 }26 }27 break;22 for (int y = 0; y < image->numRows; y++) { 23 for (int x = 0; x < image->numCols; x++) { 24 flipped->data.U16[y][size - x - 1] = image->data.U16[y][x]; 25 } 26 } 27 break; 28 28 case PS_TYPE_F32: 29 for (int y = 0; y < image->numRows; y++) {30 for (int x = 0; x < image->numCols; x++) {31 flipped->data.F32[y][size - x - 1] = image->data.F32[y][x];32 }33 }34 break;29 for (int y = 0; y < image->numRows; y++) { 30 for (int x = 0; x < image->numCols; x++) { 31 flipped->data.F32[y][size - x - 1] = image->data.F32[y][x]; 32 } 33 } 34 break; 35 35 default: 36 psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type);36 psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type); 37 37 } 38 38 … … 42 42 // Flip the chip in y 43 43 static psImage *yFlipImage(psImage *image, // Image to be flipped 44 int size// Size of image in flip dimension44 int size // Size of image in flip dimension 45 45 ) 46 46 { … … 48 48 switch (image->type.type) { 49 49 case PS_TYPE_U16: 50 for (int y = 0; y < image->numRows; y++) {51 for (int x = 0; x < image->numCols; x++) {52 flipped->data.U16[size - y - 1][x] = image->data.U16[y][x];53 }54 }55 break;50 for (int y = 0; y < image->numRows; y++) { 51 for (int x = 0; x < image->numCols; x++) { 52 flipped->data.U16[size - y - 1][x] = image->data.U16[y][x]; 53 } 54 } 55 break; 56 56 default: 57 psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type);57 psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type); 58 58 } 59 59 … … 65 65 // Return a list of the region strings and the method for interpreting them 66 66 static psList *getRegions(psString *method, // Method for evaluating the regions 67 psString methodValues // The METHOD:VALUES string67 psString methodValues // The METHOD:VALUES string 68 68 ) 69 69 { 70 70 char *values = strchr(methodValues, ':'); // The VALUES 71 71 *method = psStringNCopy(methodValues, strlen(methodValues) - strlen(values)); // The METHOD 72 psList *regionStrings = NULL; // List of the region strings72 psList *regionStrings = NULL; // List of the region strings 73 73 if (strncmp(*method, "HEADER", 6) == 0 || strncmp(*method, "HD", 2) == 0) { 74 regionStrings = papSplit(values, ", ");74 regionStrings = papSplit(values, ", "); 75 75 } else if (strncmp(*method, "VALUE", 5) == 0 || strncmp(*method, "VAL", 3) == 0) { 76 regionStrings = papSplit(values, "; ");76 regionStrings = papSplit(values, "; "); 77 77 } 78 78 … … 84 84 // Splice an image into another image 85 85 static psRegion *spliceImage(psImage *target, // Target image onto which to splice 86 int *x0, int *y0, // Starting coordinates to splice to87 psImage *source, // Source image from which to splice88 const psRegion *from, // Region to splice from89 int readdir, // Read direction90 bool xFlip, bool yFlip // Flip x and y axes?91 )86 int *x0, int *y0, // Starting coordinates to splice to 87 psImage *source, // Source image from which to splice 88 const psRegion *from, // Region to splice from 89 int readdir, // Read direction 90 bool xFlip, bool yFlip // Flip x and y axes? 91 ) 92 92 { 93 93 psImage *toSplice = psImageSubset(source, *from); // Subimage, to splice into target 94 94 if (xFlip) { 95 int size = 0;// Size of flipped image96 if (readdir == 1) {97 size = toSplice->numCols;98 } else if (readdir == 2) {99 size = target->numCols;100 }101 psImage *temp = xFlipImage(toSplice, size);102 psFree(toSplice);103 toSplice = temp;95 int size = 0; // Size of flipped image 96 if (readdir == 1) { 97 size = toSplice->numCols; 98 } else if (readdir == 2) { 99 size = target->numCols; 100 } 101 psImage *temp = xFlipImage(toSplice, size); 102 psFree(toSplice); 103 toSplice = temp; 104 104 } 105 105 if (yFlip) { 106 int size = 0; // Size of flipped image 107 if (readdir == 1) { 108 size = target->numRows; 109 } else if (readdir == 2) { 110 size = toSplice->numRows; 111 } 112 psImage *temp = yFlipImage(toSplice, size); 113 psFree(toSplice); 114 toSplice = temp; 115 } 116 106 int size = 0; // Size of flipped image 107 if (readdir == 1) { 108 size = target->numRows; 109 } else if (readdir == 2) { 110 size = toSplice->numRows; 111 } 112 psImage *temp = yFlipImage(toSplice, size); 113 psFree(toSplice); 114 toSplice = temp; 115 } 116 117 psTrace(__func__, 5, "Splicing %dx%d image in at %d,%d\n", toSplice->numCols, toSplice->numRows, 118 *x0, *y0); 117 119 (void)psImageOverlaySection(target, toSplice, *x0, *y0, "="); 118 120 119 121 if (readdir == 1) { 120 *x0 += toSplice->numCols;122 *x0 += toSplice->numCols; 121 123 } else if (readdir == 2) { 122 *y0 += toSplice->numRows;124 *y0 += toSplice->numRows; 123 125 } 124 126 … … 136 138 137 139 // Splice bias regions (overscans) 138 static bool spliceBias(psImage *splice, // Image to which to splice139 int *x0, int *y0, // Starting position for splice140 const pmCell *in, // Input cell141 pmCell *out,// Output cell142 int xFlip, int yFlip, // Flip the image?143 int readNum,// Number of readout144 int readdir// Read direction145 )140 static bool spliceBias(psImage *splice, // Image to which to splice 141 int *x0, int *y0, // Starting position for splice 142 const pmCell *in, // Input cell 143 pmCell *out, // Output cell 144 int xFlip, int yFlip, // Flip the image? 145 int readNum, // Number of readout 146 int readdir // Read direction 147 ) 146 148 { 147 149 psArray *readouts = in->readouts; // The readouts; … … 152 154 psList *outBiassecs = psListAlloc(NULL); // List of biassecs for output 153 155 while (inBiassec = psListGetAndIncrement(inBiassecsIter)) { 154 psRegion *outBiassec = spliceImage(splice, x0, y0, readout->image, inBiassec, readdir, xFlip,155 yFlip);156 psListAdd(outBiassecs, PS_LIST_TAIL, outBiassec);156 psRegion *outBiassec = spliceImage(splice, x0, y0, readout->image, inBiassec, readdir, xFlip, 157 yFlip); 158 psListAdd(outBiassecs, PS_LIST_TAIL, outBiassec); 157 159 } 158 160 psFree(inBiassecsIter); 159 161 psMetadataAdd(out->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE, 160 "BIASSEC from pmFPAMorph", outBiassecs);161 162 "BIASSEC from pmFPAMorph", outBiassecs); 163 162 164 return true; 163 165 } 164 166 167 psImage *p_pmImageMosaic(const psArray *source, // Images to splice in 168 const psVector *xFlip, const psVector *yFlip, // Need to flip x and y? 169 const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of 170 // source images 171 int xBinTarget, int yBinTarget, // Binning in x and y of target images 172 const psVector *x0, const psVector *y0 // Offsets for source images on target 173 ); 165 174 166 175 // Splice a list of cells together 176 // We need to do the following: 177 // 1. Copy the contents of the cells (i.e., the readouts) over from the in to the out 178 // 2. Splice together everything in the bucket for the image that we will write out 167 179 static psArray *spliceCells(psList *outCells, // List of target cells (required for parity info) 168 psList *inCells, // List of cells to splice together169 bool posDep// Position dependent placement of overscans?180 psList *inCells, // List of cells to splice together 181 bool posDep // Position dependent placement of overscans? 170 182 ) 171 183 { … … 176 188 177 189 if (inCells->n != outCells->n) { 178 psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n", 179 inCells->n, outCells->n); 180 return NULL; 181 } 182 183 int numCells = inCells->n; // Number of cells 184 int readdir = 0; // Read direction for CCD; currently unknown 185 int numReadouts = 0; // Number of readouts in a cell 186 psArray *spliced = NULL; // Array of spliced readouts 187 psElemType type = 0; // Image type (U16, S32, F32, etc) 188 189 psVector *xSize = NULL; // Size of spliced image in x 190 psVector *ySize = NULL; // Size of spliced image in y 190 psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n", 191 inCells->n, outCells->n); 192 return NULL; 193 } 194 195 int numCells = inCells->n; // Number of cells 196 int readdir = 0; // Read direction for CCD; currently unknown 197 int numReadouts = 0; // Number of readouts in a cell 198 psArray *spliced = NULL; // Array of spliced readouts 199 psElemType type = 0; // Image type (U16, S32, F32, etc) 191 200 192 201 psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in x? 193 202 psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in y? 194 203 psVector *xBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in x 204 psVector *yBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in y 205 psVector *x0 = psVectorAlloc(numCells, PS_TYPE_S32); // Offset in x 206 psVector *y0 = psVectorAlloc(numCells, PS_TYPE_S32); // Offset in y 207 int xBinOut = 0; // Binning in x for output 208 int yBinOut = 0; // Binning in y for output 195 209 196 210 // We go through the cells to make sure everything's kosher, and to get the sizes. 197 211 psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false); // Iterator for cells 198 pmCell *inCell = NULL; // Cell from list of input cells199 212 psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator for cells 200 int cellNum = 0; // Cell number; 213 pmCell *inCell = NULL; // Cell from list of input cells 214 pmCell *outCell = NULL; // Cell from list of output cells 215 int cellNum = -1; // Cell number 216 while (inCell = psListGetAndIncrement(inCellsIter) && outCell = psListGetAndIncrement(outCellsIter)) { 217 cellNum++; 218 psArray *inReadouts = inCell->readouts; // The readouts comprising the cell 219 220 // Get and check the read direction 221 int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR"); 222 psTrace(__func__, 10, "Checking read direction for cell %d: %d\n", cellNum, cellReaddir); 223 if (cellNum == 0) { 224 // First run through; set it. 225 readdir = cellReaddir; 226 } else if (readdir != cellReaddir) { 227 psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I " 228 "don't know how to do that!\n", readdir, cellReaddir); 229 // XXX Clean up before returning 230 return NULL; 231 } 232 233 // Get and check the number of readouts 234 if (! spliced) { 235 numReadouts = inReadouts->n; 236 spliced = psArrayAlloc(numReadouts); 237 xSize = psVectorAlloc(numReadouts, PS_TYPE_S32); 238 ySize = psVectorAlloc(numReadouts, PS_TYPE_S32); 239 for (int i = 0; i < numReadouts; i++) { 240 xSize->data.S32[i] = 0; 241 ySize->data.S32[i] = 0; 242 } 243 } else if (inReadouts->n != numReadouts) { 244 psError(PS_ERR_IO, true, "Trying to splice cells with different number of reads (%d vs %d)\n", 245 numReadouts, readouts->n); 246 // XXX Clean up before returning 247 return NULL; 248 } 249 250 // Get and check the binning --- we can't splice to cells that have different binnings, because 251 // the spliced image can't be defined. Hence all the "to" cells must have the same binning. 252 xBinIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XBIN"); 253 yBinIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YBIN"); 254 int xBinCheck = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XBIN"); 255 int yBinCheck = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YBIN"); 256 if (xBinOut == 0) { 257 xBinOut = xBinCheck; 258 } else if (xBinCheck != xBinOut) { 259 psError(PS_ERR_IO, true, "Trying to splice to cells with different x binnings (%d vs %d)\n", 260 xBinOut, xBinCheck); 261 // XXX Clean up before returning 262 return NULL; 263 } 264 if (yBinOut == 0) { 265 yBinOut = yBinCheck; 266 } else if (yBinCheck != yBinOut) { 267 psError(PS_ERR_IO, true, "Trying to splice to cells with different y binnings (%d vs %d)\n", 268 yBinOut, yBinCheck); 269 // XXX Clean up before returning 270 return NULL; 271 } 272 273 // Copy the readouts over 274 psArray *outReadouts = psArrayAlloc(inReadouts->n); // Array of readouts for the output cell 275 for (int i = 0; i < inReadouts->n; i++) { 276 outReadouts->data[i] = psMemIncrRefCounter(inReadouts->data[i]); 277 } 278 279 // Get the offsets 280 x0->data.S32[i] = psMetadataLookupS32(NULL, inCell->concepts, "CELL.X0"); 281 y0->data.S32[i] = psMetadataLookupS32(NULL, inCell->concepts, "CELL.Y0"); 282 283 // Get and check the parity 284 int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY"); 285 int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY"); 286 int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY"); 287 int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YPARITY"); 288 if (xParityIn == 0 || xParityOut == 0 || yParityIn == 0 || yParityOut == 0) { 289 psError(PS_ERR_IO, false, "Unable to determine parity of cell!\n"); 290 // XXX Clean up before returning 291 return NULL; 292 } 293 if (xParityIn != xParityOut) { 294 for (int i = 0; i < outReadouts->n; i++) { 295 psImage *image = outReadouts->data[i]; 296 outReadouts->data[i] = xFlipImage(image); 297 psFree(image); 298 } 299 // Correct CELL.X0 for the parity shift 300 psMetadataItem *x0item = psMetadataLookup(inCell->concepts, "CELL.X0"); 301 x0item->data.S32 -= xSize->data.S32[cellNum]; 302 } 303 if (yParityIn != yParityOut) { 304 for (int i = 0; i < outReadouts->n; i++) { 305 psImage *image = outReadouts->data[i]; 306 outReadouts->data[i] = xFlipImage(image); 307 psFree(image); 308 } 309 // Correct CELL.Y 0 for the parity shift 310 psMetadataItem *y0item = psMetadataLookup(inCell->concepts, "CELL.Y0"); 311 y0item->data.S32 -= ySize->data.S32[cellNum]; 312 } 313 314 ////////////////////////////////////////////////////////////////////////////////////////////////////// 315 // DO I WANT TO DO THIS???? 316 // Wouldn't it be better to splice all the images together first, and then divvy up the components? 317 ////////////////////////////////////////////////////////////////////////////////////////////////////// 318 319 320 // Manipulate the images, if required 321 for (int i = 0; i < outReadouts->n; i++) { 322 psImage *image = manipulate(outReadouts 323 324 psImage *image = psMemIncrRefCounter(outReadouts->data[i]); 325 if (xParityIn != xParityOut) { 326 psImage *image = xFlipImage(image); 327 } 328 329 330 // THIS ISN'T SO SIMPLE: we need to SET the output x0,y0 based on where the image actually gets pasted 331 // Get the offsets 332 x0->data.S32[i] = psMetadataLookupS32(NULL, outCell->concepts, "CELL.X0"); 333 y0->data.S32[i] = psMetadataLookupS32(NULL, outCell->concepts, "CELL.Y0"); 334 335 336 337 338 // Splice a list of cells together 339 static psArray *spliceCells(psList *outCells, // List of target cells (required for parity info) 340 psList *inCells, // List of cells to splice together 341 bool posDep // Position dependent placement of overscans? 342 ) 343 { 344 assert(outCells); 345 assert(inCells); 346 347 psTrace(__func__, 5, "Splicing together all cells in the bucket.\n"); 348 349 if (inCells->n != outCells->n) { 350 psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n", 351 inCells->n, outCells->n); 352 return NULL; 353 } 354 355 int numCells = inCells->n; // Number of cells 356 int readdir = 0; // Read direction for CCD; currently unknown 357 int numReadouts = 0; // Number of readouts in a cell 358 psArray *spliced = NULL; // Array of spliced readouts 359 psElemType type = 0; // Image type (U16, S32, F32, etc) 360 361 psVector *xSize = NULL; // Size of spliced image in x 362 psVector *ySize = NULL; // Size of spliced image in y 363 364 psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in x? 365 psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in y? 366 367 368 // We go through the cells to make sure everything's kosher, and to get the sizes. 369 psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false); // Iterator for cells 370 psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator for cells 371 pmCell *inCell = NULL; // Cell from list of input cells 372 int cellNum = 0; // Cell number; 201 373 while (inCell = psListGetAndIncrement(inCellsIter)) { 202 psArray *readouts = inCell->readouts; // The readouts comprising the cell 203 204 // Get and check the read direction 205 int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR"); 206 if (readdir == 0) { 207 // Zero means unknown; set it. 208 readdir = cellReaddir; 209 } else if (readdir != cellReaddir) { 210 psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I " 211 "don't know how to do that!\n", readdir, cellReaddir); 212 // XXX Clean up before returning 213 return NULL; 214 } 215 216 // Get and check the number of readouts 217 if (! spliced) { 218 numReadouts = readouts->n; 219 spliced = psArrayAlloc(numReadouts); 220 xSize = psVectorAlloc(numReadouts, PS_TYPE_S32); 221 ySize = psVectorAlloc(numReadouts, PS_TYPE_S32); 222 for (int i = 0; i < numReadouts; i++) { 223 xSize->data.S32[i] = 0; 224 ySize->data.S32[i] = 0; 225 } 226 } else if (readouts->n != numReadouts) { 227 psError(PS_ERR_IO, true, "Trying to splice cells with different number of reads (%d vs %d)\n", 228 numReadouts, readouts->n); 229 // XXX Clean up before returning 230 return NULL; 231 } 232 233 // Get and check the cell parities 234 pmCell *outCell = psListGetAndIncrement(outCellsIter); // Corresponding output cell 235 int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY"); 236 int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY"); 237 int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY"); 238 int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YPARITY"); 239 if (xParityIn == 0 || xParityOut == 0 || yParityIn == 0 || yParityOut == 0) { 240 psError(PS_ERR_IO, false, "Unable to determine parity of cell!\n"); 241 // XXX Clean up before returning 242 return NULL; 243 } 244 245 if (yParityIn == yParityOut) { 246 yFlip->data.U8[cellNum] = 0; 247 } else { 248 psTrace(__func__, 7, "Need to flip y.\n"); 249 yFlip->data.U8[cellNum] = 1; 250 } 251 if (xParityIn == xParityOut) { 252 xFlip->data.U8[cellNum] = 0; 253 } else { 254 psTrace(__func__, 7, "Need to flip x.\n"); 255 xFlip->data.U8[cellNum] = 1; 256 } 257 cellNum++; 258 259 // Calculate the sizes of the spliced images 260 for (int i = 0; i < numReadouts; i++) { 261 pmReadout *readout = readouts->data[i]; // The readout of interest 262 psImage *image = readout->image; // The image pixels 263 264 // Check image type 265 if (type == 0) { 266 type = image->type.type; 267 } else if (image->type.type != type) { 268 psError(PS_ERR_IO, true, "Trying to splice readouts with different image types (%d vs %d)\n", 269 type, image->type.type); 270 // XXX Clean up before returning 271 return NULL; 272 } 273 274 psRegion *trimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // Trim section 275 psList *biassecs = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.BIASSEC"); // Bias section 276 psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator 277 psRegion *biassec = NULL; // Bias section from list iteration 278 279 // Get the size of the spliced image 280 if (readdir == 1) { 281 xSize->data.S32[i] += trimsec->x1 - trimsec->x0; 282 ySize->data.S32[i] = MAX(ySize->data.S32[i], trimsec->y1 - trimsec->y0); 283 while (biassec = psListGetAndIncrement(biassecsIter)) { 284 xSize->data.S32[i] += biassec->x1 - biassec->x0; 285 ySize->data.S32[i] = MAX(ySize->data.S32[i], biassec->y1 - biassec->y0); 286 } 287 } else if (readdir == 2) { 288 ySize->data.S32[i] += trimsec->y1 - trimsec->y0; 289 xSize->data.S32[i] = MAX(xSize->data.S32[i], trimsec->x1 - trimsec->x0); 290 while (biassec = psListGetAndIncrement(biassecsIter)) { 291 ySize->data.S32[i] += biassec->y1 - biassec->y0; 292 xSize->data.S32[i] = MAX(ySize->data.S32[i], biassec->x1 - biassec->x0); 293 } 294 } else { 295 psError(PS_ERR_IO, true, "Invalid read direction: %d\n", readdir); 296 // XXX Clean up before returning 297 return NULL; 298 } 299 psFree(biassecsIter); 300 } 301 302 // Copy the concepts over 303 psMetadataCopy(outCell->concepts, inCell->concepts); 374 psArray *readouts = inCell->readouts; // The readouts comprising the cell 375 376 // Get and check the read direction 377 int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR"); 378 psTrace(__func__, 10, "Checking read direction for cell %d: %d\n", cellNum, cellReaddir); 379 if (cellNum == 0) { 380 // First run through; set it. 381 readdir = cellReaddir; 382 } else if (readdir != cellReaddir) { 383 psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I " 384 "don't know how to do that!\n", readdir, cellReaddir); 385 // XXX Clean up before returning 386 return NULL; 387 } 388 389 // Get and check the number of readouts 390 if (! spliced) { 391 numReadouts = readouts->n; 392 spliced = psArrayAlloc(numReadouts); 393 xSize = psVectorAlloc(numReadouts, PS_TYPE_S32); 394 ySize = psVectorAlloc(numReadouts, PS_TYPE_S32); 395 for (int i = 0; i < numReadouts; i++) { 396 xSize->data.S32[i] = 0; 397 ySize->data.S32[i] = 0; 398 } 399 } else if (readouts->n != numReadouts) { 400 psError(PS_ERR_IO, true, "Trying to splice cells with different number of reads (%d vs %d)\n", 401 numReadouts, readouts->n); 402 // XXX Clean up before returning 403 return NULL; 404 } 405 406 // Get and check the cell parities 407 pmCell *outCell = psListGetAndIncrement(outCellsIter); // Corresponding output cell 408 int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY"); 409 int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY"); 410 int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY"); 411 int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YPARITY"); 412 if (xParityIn == 0 || xParityOut == 0 || yParityIn == 0 || yParityOut == 0) { 413 psError(PS_ERR_IO, false, "Unable to determine parity of cell!\n"); 414 // XXX Clean up before returning 415 return NULL; 416 } 417 418 if (yParityIn == yParityOut) { 419 yFlip->data.U8[cellNum] = 0; 420 } else { 421 psTrace(__func__, 7, "Need to flip y.\n"); 422 yFlip->data.U8[cellNum] = 1; 423 } 424 if (xParityIn == xParityOut) { 425 xFlip->data.U8[cellNum] = 0; 426 } else { 427 psTrace(__func__, 7, "Need to flip x.\n"); 428 xFlip->data.U8[cellNum] = 1; 429 } 430 cellNum++; 431 432 // Calculate the sizes of the spliced images 433 for (int i = 0; i < numReadouts; i++) { 434 pmReadout *readout = readouts->data[i]; // The readout of interest 435 psImage *image = readout->image; // The image pixels 436 437 // Check image type 438 if (type == 0) { 439 type = image->type.type; 440 } else if (image->type.type != type) { 441 psError(PS_ERR_IO, true, "Trying to splice readouts with different image types (%d vs %d)\n", 442 type, image->type.type); 443 // XXX Clean up before returning 444 return NULL; 445 } 446 447 psRegion *trimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // Trim section 448 psList *biassecs = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.BIASSEC"); // Bias section 449 psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator 450 psRegion *biassec = NULL; // Bias section from list iteration 451 452 // Get the size of the spliced image 453 if (readdir == 1) { 454 psTrace(__func__, 7, "TRIMSEC is %.0fx%.0f\n", trimsec->x1 - trimsec->x0, 455 trimsec->y1 - trimsec->y0); 456 xSize->data.S32[i] += trimsec->x1 - trimsec->x0; 457 ySize->data.S32[i] = MAX(ySize->data.S32[i], trimsec->y1 - trimsec->y0); 458 while (biassec = psListGetAndIncrement(biassecsIter)) { 459 psTrace(__func__, 9, "BIASSEC is %.0fx%.0f\n", biassec->x1 - biassec->x0, 460 biassec->y1 - biassec->y0); 461 xSize->data.S32[i] += biassec->x1 - biassec->x0; 462 ySize->data.S32[i] = MAX(ySize->data.S32[i], biassec->y1 - biassec->y0); 463 } 464 } else if (readdir == 2) { 465 ySize->data.S32[i] += trimsec->y1 - trimsec->y0; 466 xSize->data.S32[i] = MAX(xSize->data.S32[i], trimsec->x1 - trimsec->x0); 467 while (biassec = psListGetAndIncrement(biassecsIter)) { 468 ySize->data.S32[i] += biassec->y1 - biassec->y0; 469 xSize->data.S32[i] = MAX(ySize->data.S32[i], biassec->x1 - biassec->x0); 470 } 471 } else { 472 psError(PS_ERR_IO, true, "Invalid read direction: %d\n", readdir); 473 // XXX Clean up before returning 474 return NULL; 475 } 476 psFree(biassecsIter); 477 } 304 478 } 305 479 psFree(inCellsIter); … … 308 482 309 483 // Make sure all the readouts have the same size 310 int numRows = ySize->data.S32[0]; // Number of rows for spliced image311 int numCols = xSize->data.S32[0]; // Number of columns for spliced image484 int numRows = ySize->data.S32[0]; // Number of rows for spliced image 485 int numCols = xSize->data.S32[0]; // Number of columns for spliced image 312 486 for (int i = 1; i < numReadouts; i++) { 313 if (xSize->data.S32[i] != numCols || ySize->data.S32[i] != numRows) {314 psError(PS_ERR_IO, true, "Spliced readouts would have different sizes: %dx%d vs %dx%d\n",315 numCols, numRows, xSize->data.S32[i], ySize->data.S32[i]);316 }487 if (xSize->data.S32[i] != numCols || ySize->data.S32[i] != numRows) { 488 psError(PS_ERR_IO, true, "Spliced readouts would have different sizes: %dx%d vs %dx%d\n", 489 numCols, numRows, xSize->data.S32[i], ySize->data.S32[i]); 490 } 317 491 } 318 492 … … 321 495 // Now we have done the requisite checks, and know the sizes; we just go through and splice together 322 496 for (int i = 0; i < numReadouts; i++) { 323 psImage *splice = psImageAlloc(numCols, numRows, type); // The spliced image 324 int x0 = 0, y0 = 0; // Position at which the splice begins 325 326 // Position-dependent overscans first 327 if (posDep) { 328 for (int cellNum = 0; cellNum < inCells->n / 2; cellNum++) { 329 pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest 330 pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest 331 spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, 332 readdir); 333 } 334 } 335 336 // Then the images 337 psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false);// Iterator for cells 338 pmCell *inCell = NULL; // Cell from the list of cells 339 psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator 340 int cellNum = 0; // Cell number 341 while (inCell = psListGetAndIncrement(inCellsIter)) { 342 psArray *readouts = inCell->readouts; // The readouts comprising the cell 343 pmReadout *readout = readouts->data[i]; // The specific readout of interest 344 pmCell *outCell = psListGetAndIncrement(outCellsIter); // Output cell 345 psRegion *inTrimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // TRIMSEC 346 psRegion *outTrimsec = spliceImage(splice, &x0, &y0, readout->image, inTrimsec, readdir, 347 xFlip->data.U8[cellNum], yFlip->data.U8[cellNum]); 348 349 // Update the output TRIMSEC 350 psMetadataAdd(outCell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE, 351 "TRIMSEC from pmFPAMorph", outTrimsec); 352 } 353 psFree(inCellsIter); 354 psFree(outCellsIter); 355 356 // Then the biases 357 if (! posDep) { 358 // Need to only do half the biases 359 for (int cellNum = inCells->n / 2; cellNum < inCells->n; cellNum++) { 360 pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest 361 pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest 362 spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, 363 readdir); 364 } 365 } else { 366 // Need to do all the biases 367 for (int cellNum = 0; cellNum < inCells->n; cellNum++) { 368 pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest 369 pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest 370 spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, 371 readdir); 372 } 373 } 374 375 spliced->data[i] = splice; 497 psImage *splice = psImageAlloc(numCols, numRows, type); // The spliced image 498 int x0 = 0, y0 = 0; // Position at which the splice begins 499 500 // Position-dependent overscans first 501 if (posDep) { 502 psTrace(__func__, 8, "Doing first position-dependent biases...\n"); 503 for (int cellNum = 0; cellNum < inCells->n / 2; cellNum++) { 504 psTrace(__func__, 9, "Cell %d...\n", cellNum); 505 pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest 506 pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest 507 spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, 508 readdir); 509 } 510 } 511 512 // Then the images 513 psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false);// Iterator for cells 514 pmCell *inCell = NULL; // Cell from the list of cells 515 psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator 516 int cellNum = 0; // Cell number 517 while (inCell = psListGetAndIncrement(inCellsIter)) { 518 psArray *readouts = inCell->readouts; // The readouts comprising the cell 519 pmReadout *readout = readouts->data[i]; // The specific readout of interest 520 pmCell *outCell = psListGetAndIncrement(outCellsIter); // Output cell 521 psRegion *inTrimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // TRIMSEC 522 psRegion *outTrimsec = spliceImage(splice, &x0, &y0, readout->image, inTrimsec, readdir, 523 xFlip->data.U8[cellNum], yFlip->data.U8[cellNum]); 524 525 // Update the output TRIMSEC 526 psMetadataAdd(outCell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE, 527 "TRIMSEC from pmFPAMorph", outTrimsec); 528 } 529 psFree(inCellsIter); 530 psFree(outCellsIter); 531 532 // Then the biases 533 if (posDep) { 534 // Need to only do half the biases 535 psTrace(__func__, 8, "Doing second position-dependent biases...\n"); 536 for (int cellNum = inCells->n / 2; cellNum < inCells->n; cellNum++) { 537 psTrace(__func__, 9, "Cell %d...\n", cellNum); 538 pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest 539 pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest 540 spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, 541 readdir); 542 } 543 } else { 544 // Need to do all the biases 545 for (int cellNum = 0; cellNum < inCells->n; cellNum++) { 546 pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest 547 pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest 548 spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i, 549 readdir); 550 } 551 } 552 553 spliced->data[i] = splice; 376 554 377 555 } // Iterating over readouts … … 384 562 385 563 386 bool pmFPAMorph(pmFPA *toFPA, // FPA structure to which to morph387 pmFPA *fromFPA,// FPA structure from which to morph388 bool positionDependent, // Is the position of the overscan dependent on the position?389 int chipNum,// Chip number, in case the scopes are different390 int cellNum// Cell number, in case the scopes are different564 bool pmFPAMorph(pmFPA *toFPA, // FPA structure to which to morph 565 pmFPA *fromFPA, // FPA structure from which to morph 566 bool positionDependent, // Is the position of the overscan dependent on the position? 567 int chipNum, // Chip number, in case the scopes are different 568 int cellNum // Cell number, in case the scopes are different 391 569 ) 392 570 { … … 394 572 psList *sourceCells = psListAlloc(NULL); // List of source cells 395 573 396 psArray *toChips = toFPA->chips; // Array of chips 397 psArray *fromChips = fromFPA->chips;// Array of chips 574 psArray *toChips = toFPA->chips; // Array of chips 575 psArray *fromChips = fromFPA->chips;// Array of chips 576 577 psTrace(__func__, 1, "Copying FPA concepts over...\n"); 578 pap_psMetadataCopy(toFPA->concepts, fromFPA->concepts); 579 pap_psMetadataCopy(toFPA->phu, fromFPA->phu); 398 580 399 581 for (int i = 0; i < toChips->n; i++) { 400 pmChip *toChip = toChips->data[i]; 401 pmChip *fromChip = NULL; 402 // Select the correct chip 403 if (toChips->n == 1 && fromChips->n > chipNum) { 404 fromChip = fromChips->data[chipNum]; 405 } else if (fromChips->n == 1 && toChips->n > chipNum) { 406 fromChip = fromChips->data[0]; 407 } else if (toChips->n == fromChips->n) { 408 fromChip = fromChips->data[i]; 409 } else { 410 psError(PS_ERR_IO, true, "Unable to discern intended chip.\n"); 411 return false; 412 } 413 414 psArray *toCells = toChip->cells; // Array of cells 415 psArray *fromCells = fromChip->cells; // Array of cells 416 417 for (int j = 0; j < toCells->n; j++) { 418 pmCell *toCell = toCells->data[j]; 419 pmCell *fromCell = NULL; 420 // Select the correct cell 421 if (toCells->n == 1 && fromCells->n > cellNum) { 422 fromCell = fromCells->data[cellNum]; 423 } else if (fromCells->n == 1 && toCells->n > cellNum) { 424 fromCell = fromCells->data[0]; 425 } else if (toCells->n == fromCells->n) { 426 fromCell = fromCells->data[j]; 427 } else { 428 psError(PS_ERR_IO, true, "Unable to discern intended cell.\n"); 429 return false; 430 } 431 432 psTrace(__func__, 5, "Putting chip %d, cell %d in the bucket.\n", i, j); 433 psListAdd(targetCells, PS_LIST_TAIL, toCell); 434 psListAdd(sourceCells, PS_LIST_TAIL, fromCell); 435 436 #if 1 437 int xParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.XPARITY"); 438 int yParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.YPARITY"); 439 int xParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.XPARITY"); 440 int yParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.YPARITY"); 441 psTrace(__func__, 3, "Chip %d, cell %d: in (%d,%d) out (%d,%d)\n", i, j, xParityIn, yParityIn, 442 xParityOut, yParityOut); 582 pmChip *toChip = toChips->data[i]; 583 pmChip *fromChip = NULL; 584 // Select the correct chip 585 if (toChips->n == 1 && fromChips->n > chipNum) { 586 fromChip = fromChips->data[chipNum]; 587 } else if (fromChips->n == 1 && toChips->n > chipNum) { 588 fromChip = fromChips->data[0]; 589 } else if (toChips->n == fromChips->n) { 590 fromChip = fromChips->data[i]; 591 } else { 592 psError(PS_ERR_IO, true, "Unable to discern intended chip.\n"); 593 return false; 594 } 595 596 if (! fromChip->valid) { 597 toChip->valid = false; 598 continue; 599 } 600 601 psTrace(__func__, 2, "Copying chip %d concepts over...\n", i); 602 pap_psMetadataCopy(toChip->concepts, fromChip->concepts); 603 604 psArray *toCells = toChip->cells; // Array of cells 605 psArray *fromCells = fromChip->cells; // Array of cells 606 607 for (int j = 0; j < toCells->n; j++) { 608 pmCell *toCell = toCells->data[j]; 609 pmCell *fromCell = NULL; 610 // Select the correct cell 611 if (toCells->n == 1 && fromCells->n > cellNum) { 612 fromCell = fromCells->data[cellNum]; 613 } else if (fromCells->n == 1 && toCells->n > cellNum) { 614 fromCell = fromCells->data[0]; 615 } else if (toCells->n == fromCells->n) { 616 fromCell = fromCells->data[j]; 617 } else { 618 psError(PS_ERR_IO, true, "Unable to discern intended cell.\n"); 619 return false; 620 } 621 622 if (! fromCell->valid) { 623 toCell->valid = false; 624 continue; 625 } 626 627 #ifdef TESTING 628 int xParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.XPARITY"); 629 int yParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.YPARITY"); 630 int xParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.XPARITY"); 631 int yParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.YPARITY"); 632 psTrace(__func__, 3, "Chip %d, cell %d: in (%d,%d) out (%d,%d)\n", i, j, xParityIn, yParityIn, 633 xParityOut, yParityOut); 443 634 #endif 444 635 445 // Copy the cell contents over 446 toCell->readouts = psMemIncrRefCounter(fromCell->readouts); 447 448 if (toCell->hdu && strlen(toCell->hdu->extname) > 0) { 449 // Splice the component cells 450 p_pmHDU *hdu = toCell->hdu; 451 if (! hdu->header) { 452 hdu->header = psMetadataAlloc(); 453 } 454 hdu->images = spliceCells(targetCells, sourceCells, positionDependent); 455 // Purge the lists 456 while (psListRemove(targetCells, PS_LIST_TAIL)); 457 while (psListRemove(sourceCells, PS_LIST_TAIL)); 458 } 459 460 } 461 462 if (toChip->hdu && strlen(toChip->hdu->extname) > 0) { 463 // Splice the component cells 464 p_pmHDU *hdu = toChip->hdu; 465 if (! hdu->header) { 466 hdu->header = psMetadataAlloc(); 467 } 468 hdu->images = spliceCells(targetCells, sourceCells, positionDependent); 469 // Purge the lists 470 while (psListRemove(targetCells, PS_LIST_TAIL)); 471 while (psListRemove(sourceCells, PS_LIST_TAIL)); 472 } 636 #ifdef TESTING 637 psMetadataPrint(toCell->concepts, 7); 638 #endif 639 640 // Copy the concepts over 641 642 // Need to be a little tricky here --- some concepts are already set from the camera configuration 643 // file (e.g., CELL.NAME), and we want to preserve these, not replace it with the old value. 644 psTrace(__func__, 3, "Copying cell %d concepts over...\n", j); 645 psMetadata *newConcepts = pap_psMetadataCopy(NULL, fromCell->concepts); 646 pap_psMetadataCopy(newConcepts, toCell->concepts); // This preserves any already existing concepts 647 psFree(toCell->concepts); 648 toCell->concepts = newConcepts; 649 650 // Now, we need to check if the following concepts for the target cell are defined statically: 651 // CELL.XPARITY, CELL.YPARITY, CELL.XBIN, CELL.YBIN, CELL.X0, CELL.Y0 652 // 1. If they are specified by the header, we can do what we want (because we can set the header). 653 // 2. If they are not specified by the header, we use those values. 654 const char *checkConcepts[4] = { "CELL.XPARITY", 655 "CELL.YPARITY", 656 "CELL.XBIN", 657 "CELL.YBIN", 658 "CELL.X0", 659 "CELL.Y0" }; // Concepts that we need to check 660 bool mdok = false; // Result of MD lookup 661 psMetadata *translation = psMetadataLookupMD(&mdok, toFPA->camera, "TRANSLATION"); // Header 662 // translation table 663 if (!mdok || ! translation) { 664 psError(PS_ERR_IO, true, "Unable to find TRANSLATION in camera configuration!\n"); 665 return false; 666 } 667 for (int c = 0; c < 4; c++) { 668 psString headerName = psMetadataLookupString(&mdok, translation, checkConcepts[c]); // Name of 669 // header holding the concept, or NULL 670 if (!mdok || strlen(headerName) <= 0) { 671 p_pmCellIngestConcept(toFPA, toChip, toCell, db, checkConcepts[c]); 672 } 673 } 674 675 #ifdef TESTING 676 psMetadataPrint(toCell->concepts, 7); 677 #endif 678 679 psTrace(__func__, 5, "Putting chip %d, cell %d in the bucket.\n", i, j); 680 psListAdd(targetCells, PS_LIST_TAIL, toCell); 681 psListAdd(sourceCells, PS_LIST_TAIL, fromCell); 682 683 // Copy the cell contents over 684 toCell->readouts = psMemIncrRefCounter(fromCell->readouts); 685 686 if (toCell->hdu && strlen(toCell->hdu->extname) > 0) { 687 // Splice the component cells 688 p_pmHDU *hdu = toCell->hdu; 689 if (! hdu->header) { 690 hdu->header = psMetadataAlloc(); 691 } 692 hdu->images = spliceCells(targetCells, sourceCells, positionDependent); 693 // Purge the lists 694 while (psListRemove(targetCells, PS_LIST_TAIL)); 695 while (psListRemove(sourceCells, PS_LIST_TAIL)); 696 } 697 698 } 699 700 if (toChip->hdu && strlen(toChip->hdu->extname) > 0) { 701 // Splice the component cells 702 p_pmHDU *hdu = toChip->hdu; 703 if (! hdu->header) { 704 hdu->header = psMetadataAlloc(); 705 } 706 hdu->images = spliceCells(targetCells, sourceCells, positionDependent); 707 // Purge the lists 708 while (psListRemove(targetCells, PS_LIST_TAIL)); 709 while (psListRemove(sourceCells, PS_LIST_TAIL)); 710 } 473 711 474 712 } 475 713 476 714 if (toFPA->hdu && strlen(toFPA->hdu->extname) > 0) { 477 // Splice the component cells478 p_pmHDU *hdu = toFPA->hdu;479 if (! hdu->header) {480 hdu->header = psMetadataAlloc();481 }482 hdu->images = spliceCells(targetCells, sourceCells, positionDependent);483 // Purge the lists484 while (psListRemove(targetCells, PS_LIST_TAIL));485 while (psListRemove(sourceCells, PS_LIST_TAIL));715 // Splice the component cells 716 p_pmHDU *hdu = toFPA->hdu; 717 if (! hdu->header) { 718 hdu->header = psMetadataAlloc(); 719 } 720 hdu->images = spliceCells(targetCells, sourceCells, positionDependent); 721 // Purge the lists 722 while (psListRemove(targetCells, PS_LIST_TAIL)); 723 while (psListRemove(sourceCells, PS_LIST_TAIL)); 486 724 } 487 725 -
trunk/archive/scripts/src/phase2/pmFPARead.c
r5633 r5786 7 7 #include "pmFPARead.h" 8 8 9 #include "papStuff.h" // For "split"9 #include "papStuff.h" // For "split" 10 10 11 11 // NOTE: Need to deal with header inheritance … … 18 18 // Read a FITS extension into a chip 19 19 static bool readExtension(p_pmHDU *hdu, // Pixel data into which to read 20 psFits *fits// The FITS file from which to read20 psFits *fits // The FITS file from which to read 21 21 ) 22 22 { … … 25 25 psTrace(__func__, 7, "Moving to extension %s...\n", extName); 26 26 if (strncmp(extName, "PHU", 3) == 0) { 27 if (!psFitsMoveExtNum(fits, 0, false)) {28 psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n");29 return false;30 }27 if (!psFitsMoveExtNum(fits, 0, false)) { 28 psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n"); 29 return false; 30 } 31 31 } else if (! psFitsMoveExtName(fits, extName)) { 32 psError(PS_ERR_IO, false, "Unable to find extension %s in FITS file!\n", extName);33 return false;32 psError(PS_ERR_IO, false, "Unable to find extension %s in FITS file!\n", extName); 33 return false; 34 34 } 35 35 psTrace(__func__, 7, "Reading header....\n"); 36 36 psMetadata *header = psFitsReadHeader(NULL, fits); // Header 37 37 if (! header) { 38 psError(PS_ERR_IO, false, "Unable to read FITS header!\n");39 return false;38 psError(PS_ERR_IO, false, "Unable to read FITS header!\n"); 39 return false; 40 40 } 41 41 … … 44 44 int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); 45 45 if (!mdStatus) { 46 psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header of extension %s!\n",47 extName);46 psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header of extension %s!\n", 47 extName); 48 48 } 49 49 if (nAxis != 2 && nAxis != 3) { 50 psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into a single image "51 "anyway.\n");50 psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into a single image " 51 "anyway.\n"); 52 52 } 53 53 psTrace(__func__, 9, "NAXIS = %d\n", nAxis); 54 54 55 int numPlanes = 1; // Number of planes55 int numPlanes = 1; // Number of planes 56 56 if (nAxis == 3) { 57 numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3");58 if (!mdStatus) {59 psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n");60 return false;61 }62 } 63 psRegion region = {0, 0, 0, 0}; // Region to read is everything57 numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); 58 if (!mdStatus) { 59 psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); 60 return false; 61 } 62 } 63 psRegion region = {0, 0, 0, 0}; // Region to read is everything 64 64 65 65 // Read each plane into the array … … 67 67 psArray *pixels = psArrayAlloc(numPlanes); // Array of images 68 68 for (int i = 0; i < numPlanes; i++) { 69 psTrace(__func__, 9, "Reading plane %d\n", i);70 psMemCheckCorruption(true);71 psImage *image = psFitsReadImage(NULL, fits, region, i);72 73 // XXX: Type conversion here to support the modules, which don't have multiple type support yet74 if (image->type.type != PS_TYPE_F32) {75 pixels->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);69 psTrace(__func__, 9, "Reading plane %d\n", i); 70 psMemCheckCorruption(true); 71 psImage *image = psFitsReadImage(NULL, fits, region, i); 72 73 // XXX: Type conversion here to support the modules, which don't have multiple type support yet 74 if (image->type.type != PS_TYPE_F32) { 75 pixels->data[i] = psImageCopy(NULL, image, PS_TYPE_F32); 76 76 77 77 #ifndef PRODUCTION 78 // XXX: Temporary fix for writing images, until psFits gets cleaned up79 psMetadataItem *bitpixItem = psMetadataLookup(header, "BITPIX");80 bitpixItem->data.S32 = -32;81 psMetadataRemove(header, 0, "BZERO");82 psMetadataRemove(header, 0, "BSCALE");78 // XXX: Temporary fix for writing images, until psFits gets cleaned up 79 psMetadataItem *bitpixItem = psMetadataLookup(header, "BITPIX"); 80 bitpixItem->data.S32 = -32; 81 psMetadataRemove(header, 0, "BZERO"); 82 psMetadataRemove(header, 0, "BSCALE"); 83 83 #endif 84 84 85 psFree(image);86 } else {87 pixels->data[i] = image;88 }89 psTrace(__func__, 10, "Done\n");85 psFree(image); 86 } else { 87 pixels->data[i] = image; 88 } 89 psTrace(__func__, 10, "Done\n"); 90 90 if (! pixels->data[i]) { 91 psError(PS_ERR_IO, false, "Unable to read image for extension %s, plane %d\n", extName, i);92 return false;93 }91 psError(PS_ERR_IO, false, "Unable to read image for extension %s, plane %d\n", extName, i); 92 return false; 93 } 94 94 } 95 95 … … 105 105 // Portion out an image into the cell 106 106 static bool generateReadouts(pmCell *cell, // The cell that gets its bits 107 p_pmHDU *hdu // Pixel data, containing image, mask, weights107 p_pmHDU *hdu // Pixel data, containing image, mask, weights 108 108 ) 109 109 { 110 psArray *images = hdu->images; // Array of images (each of which is a readout)111 psArray *masks = hdu->masks; // Array of masks (one for each readout)110 psArray *images = hdu->images; // Array of images (each of which is a readout) 111 psArray *masks = hdu->masks; // Array of masks (one for each readout) 112 112 // Iterate over each of the image planes 113 113 for (int i = 0; i < images->n; i++) { 114 psImage *image = images->data[i]; // The i-th plane115 psImage *mask = NULL;// The mask116 if (masks) {117 mask = masks->data[i];118 }119 120 pmReadout *readout = pmReadoutAlloc(cell, image, mask, 0,0,1,1);121 psFree(readout);// This seems silly, but the alloc registers it with the cell, so we122 // don't have to do anything with it. Perhaps we should change123 // pmReadoutAlloc to pmReadoutAdd?114 psImage *image = images->data[i]; // The i-th plane 115 psImage *mask = NULL; // The mask 116 if (masks) { 117 mask = masks->data[i]; 118 } 119 120 pmReadout *readout = pmReadoutAlloc(cell, image, mask, 0,0,1,1); 121 psFree(readout); // This seems silly, but the alloc registers it with the cell, so we 122 // don't have to do anything with it. Perhaps we should change 123 // pmReadoutAlloc to pmReadoutAdd? 124 124 } 125 125 … … 131 131 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 132 132 133 bool pmFPARead(pmFPA *fpa, // FPA to read into134 psFits *fits,// FITS file from which to read135 psMetadata *phu,// Primary header136 psDB *db// Database handle, for concept ingest133 bool pmFPARead(pmFPA *fpa, // FPA to read into 134 psFits *fits, // FITS file from which to read 135 psMetadata *phu, // Primary header 136 psDB *db // Database handle, for concept ingest 137 137 ) 138 138 { 139 p_pmHDU *hdu = NULL; // Pixel data from FITS file139 p_pmHDU *hdu = NULL; // Pixel data from FITS file 140 140 141 141 // Read the PHU, if required 142 142 if (! phu) { 143 if (! psFitsMoveExtNum(fits, 0, false)) { 144 psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n"); 145 return false; 146 } 147 psTrace(__func__, 7, "Reading PHU....\n"); 148 psMetadata *phu = psFitsReadHeader(NULL, fits); // Primary header 149 } 150 fpa->phu = phu; 143 if (! psFitsMoveExtNum(fits, 0, false)) { 144 psError(PS_ERR_IO, false, "Unable to find PHU in FITS file!\n"); 145 return false; 146 } 147 psTrace(__func__, 7, "Reading PHU....\n"); 148 fpa->phu = psFitsReadHeader(NULL, fits); // Primary header 149 } else { 150 fpa->phu = psMemIncrRefCounter(phu); 151 } 151 152 152 153 // Read in.... 153 154 psTrace(__func__, 1, "Working on FPA...\n"); 154 155 if (fpa->hdu) { 155 hdu = fpa->hdu;156 psTrace(__func__, 3, "Reading pixels from extension %s into FPA.\n", hdu->extname);157 if (! readExtension(hdu, fits)) {158 psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname);159 return false;160 }156 hdu = fpa->hdu; 157 psTrace(__func__, 3, "Reading pixels from extension %s into FPA.\n", hdu->extname); 158 if (! readExtension(hdu, fits)) { 159 psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname); 160 return false; 161 } 161 162 } 162 163 pmFPAIngestConcepts(fpa, db); 163 164 164 psArray *chips = fpa->chips; // Array of chips165 psArray *chips = fpa->chips; // Array of chips 165 166 // Iterate over the FPA 166 167 for (int i = 0; i < chips->n; i++) { 167 pmChip *chip = chips->data[i]; // The chip168 169 // Only read chips marked "valid"170 if (! chip->valid) {171 psTrace(__func__, 2, "Ignoring chip %d...\n", i);172 continue;173 }174 psTrace(__func__, 2, "Reading in chip %d...\n", i);175 176 if (chip->hdu) {177 hdu = chip->hdu;178 psTrace(__func__, 3, "Reading pixels from extension %s into chip %d.\n", hdu->extname, i);179 if (! readExtension(hdu, fits)) {180 psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname);181 return false;182 }183 }184 pmChipIngestConcepts(chip, db);185 186 // Iterate over the chip187 psArray *cells = chip->cells;// Array of cells188 for (int j = 0; j < cells->n; j++) {189 pmCell *cell = cells->data[j]; // The cell190 191 // Only read cells marked "valid"192 if (! cell->valid) {193 psTrace(__func__, 3, "Ignoring chip %d...\n", i);194 continue;195 }196 psTrace(__func__, 3, "Reading in cell %d...\n", j);197 198 if (cell->hdu) {199 hdu = cell->hdu;200 psTrace(__func__, 5, "Reading pixels from extension %s into cell %d.\n", hdu->extname,201 j);202 if (! readExtension(hdu, fits)) {203 psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n",204 hdu->extname);205 return false;206 }207 }208 pmCellIngestConcepts(cell, db);209 210 psTrace(__func__, 5, "Allocating readouts for chip %d cell %d...\n",211 i, j);212 generateReadouts(cell, hdu);213 }168 pmChip *chip = chips->data[i]; // The chip 169 170 // Only read chips marked "valid" 171 if (! chip->valid) { 172 psTrace(__func__, 2, "Ignoring chip %d...\n", i); 173 continue; 174 } 175 psTrace(__func__, 2, "Reading in chip %d...\n", i); 176 177 if (chip->hdu) { 178 hdu = chip->hdu; 179 psTrace(__func__, 3, "Reading pixels from extension %s into chip %d.\n", hdu->extname, i); 180 if (! readExtension(hdu, fits)) { 181 psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", hdu->extname); 182 return false; 183 } 184 } 185 pmChipIngestConcepts(chip, db); 186 187 // Iterate over the chip 188 psArray *cells = chip->cells; // Array of cells 189 for (int j = 0; j < cells->n; j++) { 190 pmCell *cell = cells->data[j]; // The cell 191 192 // Only read cells marked "valid" 193 if (! cell->valid) { 194 psTrace(__func__, 3, "Ignoring chip %d...\n", i); 195 continue; 196 } 197 psTrace(__func__, 3, "Reading in cell %d...\n", j); 198 199 if (cell->hdu) { 200 hdu = cell->hdu; 201 psTrace(__func__, 5, "Reading pixels from extension %s into cell %d.\n", hdu->extname, 202 j); 203 if (! readExtension(hdu, fits)) { 204 psError(PS_ERR_IO, false, "Unable to read pixels from extension %s\n", 205 hdu->extname); 206 return false; 207 } 208 } 209 pmCellIngestConcepts(cell, db); 210 211 psTrace(__func__, 5, "Allocating readouts for chip %d cell %d...\n", 212 i, j); 213 generateReadouts(cell, hdu); 214 } 214 215 } 215 216 … … 219 220 // Translate a name from the configuration file, containing something like "%a_%d" to a real value 220 221 psString p_pmFPATranslateName(psString name, // The name to translate 221 pmCell *cell // The cell for which to translate222 pmCell *cell // The cell for which to translate 222 223 ) 223 224 { … … 230 231 231 232 psString translation = psMemIncrRefCounter(name); // The translated string 232 char *temp = NULL; // Temporary string233 234 pmChip *chip = cell->parent; // Chip of interest235 pmFPA *fpa = chip->parent; // FPA of interest233 char *temp = NULL; // Temporary string 234 235 pmChip *chip = cell->parent; // Chip of interest 236 pmFPA *fpa = chip->parent; // FPA of interest 236 237 237 238 // FPA.NAME 238 239 if (temp = strstr(translation, "%a")) { 239 // This is not particularly friendly to the memory allocation, but the cache should make it OK,240 // and there's not much of it anyway...241 psFree(translation);242 translation = psStringNCopy(name, strlen(name) - strlen(temp));// Copy first part of string243 psString fpaName = psMetadataLookupString(NULL, fpa->concepts, "FPA.NAME"); // The value of FPA.NAME244 psStringAppend(&translation, "%s%s", fpaName, temp + 2);245 // So "translation" now contains the first part, the replaced string, and the last part.240 // This is not particularly friendly to the memory allocation, but the cache should make it OK, 241 // and there's not much of it anyway... 242 psFree(translation); 243 translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string 244 psString fpaName = psMetadataLookupString(NULL, fpa->concepts, "FPA.NAME"); // The value of FPA.NAME 245 psStringAppend(&translation, "%s%s", fpaName, temp + 2); 246 // So "translation" now contains the first part, the replaced string, and the last part. 246 247 } 247 248 248 249 // CHIP.NAME 249 250 if (temp = strstr(translation, "%b")) { 250 // This is not particularly friendly to the memory allocation, but the cache should make it OK,251 // and there's not much of it anyway...252 psFree(translation);253 translation = psStringNCopy(name, strlen(name) - strlen(temp));// Copy first part of string254 psString chipName = psMetadataLookupString(NULL, chip->concepts, "CHIP.NAME"); // CHIP.NAME's value255 psStringAppend(&translation, "%s%s", chipName, temp + 2);256 // So "translation" now contains the first part, the replaced string, and the last part.251 // This is not particularly friendly to the memory allocation, but the cache should make it OK, 252 // and there's not much of it anyway... 253 psFree(translation); 254 translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string 255 psString chipName = psMetadataLookupString(NULL, chip->concepts, "CHIP.NAME"); // CHIP.NAME's value 256 psStringAppend(&translation, "%s%s", chipName, temp + 2); 257 // So "translation" now contains the first part, the replaced string, and the last part. 257 258 } 258 259 259 260 // CELL.NAME 260 261 if (temp = strstr(translation, "%c")) { 261 // This is not particularly friendly to the memory allocation, but the cache should make it OK,262 // and there's not much of it anyway...263 psFree(translation);264 translation = psStringNCopy(name, strlen(name) - strlen(temp));// Copy first part of string265 psString cellName = psMetadataLookupString(NULL, cell->concepts, "CELL.NAME"); // CELL.NAME's value266 psStringAppend(&translation, "%s%s", cellName, temp + 2);267 // So "translation" now contains the first part, the replaced string, and the last part.262 // This is not particularly friendly to the memory allocation, but the cache should make it OK, 263 // and there's not much of it anyway... 264 psFree(translation); 265 translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string 266 psString cellName = psMetadataLookupString(NULL, cell->concepts, "CELL.NAME"); // CELL.NAME's value 267 psStringAppend(&translation, "%s%s", cellName, temp + 2); 268 // So "translation" now contains the first part, the replaced string, and the last part. 268 269 } 269 270 270 271 // Chip number 271 272 if (temp = strstr(translation, "%d")) { 272 // This is not particularly friendly to the memory allocation, but the cache should make it OK,273 // and there's not much of it anyway...274 psFree(translation);275 translation = psStringNCopy(name, strlen(name) - strlen(temp));// Copy first part of string276 // Search for the pointer to get the chip number277 int chipNum = -1;278 psArray *chips = fpa->chips;// The array of chips279 for (int i = 0; i < chips->n && chipNum < 0; i++) {280 if (chips->data[i] == chip) {281 chipNum = i;282 }283 }284 if (chipNum < 0) {285 psError(PS_ERR_IO, true, "Unable to find chip to get name: %s\n", name);286 // Try to muddle on by leaving the number out287 psStringAppend(&translation, "%s", temp + 2);288 } else {289 psStringAppend(&translation, "%d%s", chipNum, temp + 2);290 }291 // So "translation" now contains the first part, the replaced string, and the last part.273 // This is not particularly friendly to the memory allocation, but the cache should make it OK, 274 // and there's not much of it anyway... 275 psFree(translation); 276 translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string 277 // Search for the pointer to get the chip number 278 int chipNum = -1; 279 psArray *chips = fpa->chips; // The array of chips 280 for (int i = 0; i < chips->n && chipNum < 0; i++) { 281 if (chips->data[i] == chip) { 282 chipNum = i; 283 } 284 } 285 if (chipNum < 0) { 286 psError(PS_ERR_IO, true, "Unable to find chip to get name: %s\n", name); 287 // Try to muddle on by leaving the number out 288 psStringAppend(&translation, "%s", temp + 2); 289 } else { 290 psStringAppend(&translation, "%d%s", chipNum, temp + 2); 291 } 292 // So "translation" now contains the first part, the replaced string, and the last part. 292 293 } 293 294 294 295 // Cell number 295 296 if (temp = strstr(translation, "%e")) { 296 // This is not particularly friendly to the memory allocation, but the cache should make it OK,297 // and there's not much of it anyway...298 psFree(translation);299 translation = psStringNCopy(name, strlen(name) - strlen(temp));// Copy first part of string300 // Search for the pointer to get the cell number301 int cellNum = -1;302 psArray *cells = chip->cells;// The array of cells303 for (int i = 0; i < cells->n && cellNum < 0; i++) {304 if (cells->data[i] == cell) {305 cellNum = i;306 }307 }308 if (cellNum < 0) {309 psError(PS_ERR_IO, true, "Unable to find cell to get name: %s\n", name);310 // Try to muddle on by leaving the number out311 psStringAppend(&translation, "%s", temp + 2);312 } else {313 psStringAppend(&translation, "%d%s", cellNum, temp + 2);314 }315 // So "translation" now contains the first part, the replaced string, and the last part.297 // This is not particularly friendly to the memory allocation, but the cache should make it OK, 298 // and there's not much of it anyway... 299 psFree(translation); 300 translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string 301 // Search for the pointer to get the cell number 302 int cellNum = -1; 303 psArray *cells = chip->cells; // The array of cells 304 for (int i = 0; i < cells->n && cellNum < 0; i++) { 305 if (cells->data[i] == cell) { 306 cellNum = i; 307 } 308 } 309 if (cellNum < 0) { 310 psError(PS_ERR_IO, true, "Unable to find cell to get name: %s\n", name); 311 // Try to muddle on by leaving the number out 312 psStringAppend(&translation, "%s", temp + 2); 313 } else { 314 psStringAppend(&translation, "%d%s", cellNum, temp + 2); 315 } 316 // So "translation" now contains the first part, the replaced string, and the last part. 316 317 } 317 318 318 319 // Extension name 319 320 if (temp = strstr(translation, "%f")) { 320 // This is not particularly friendly to the memory allocation, but the cache should make it OK,321 // and there's not much of it anyway...322 psFree(translation);323 translation = psStringNCopy(name, strlen(name) - strlen(temp));// Copy first part of string324 const char *extname = NULL;325 if (cell->hdu) {326 extname = cell->hdu->extname;327 } else if (chip->hdu) {328 extname = chip->hdu->extname;329 } else if (fpa->hdu) {330 extname = fpa->hdu->extname;331 }332 psStringAppend(&translation, "%s%s", extname, temp + 2);333 // So "translation" now contains the first part, the replaced string, and the last part.321 // This is not particularly friendly to the memory allocation, but the cache should make it OK, 322 // and there's not much of it anyway... 323 psFree(translation); 324 translation = psStringNCopy(name, strlen(name) - strlen(temp)); // Copy first part of string 325 const char *extname = NULL; 326 if (cell->hdu) { 327 extname = cell->hdu->extname; 328 } else if (chip->hdu) { 329 extname = chip->hdu->extname; 330 } else if (fpa->hdu) { 331 extname = fpa->hdu->extname; 332 } 333 psStringAppend(&translation, "%s%s", extname, temp + 2); 334 // So "translation" now contains the first part, the replaced string, and the last part. 334 335 } 335 336 … … 338 339 339 340 // Read a mask into the FPA 340 bool pmFPAReadMask(pmFPA *fpa, // FPA to read into341 psFits *source// Source FITS file (for the original data)341 bool pmFPAReadMask(pmFPA *fpa, // FPA to read into 342 psFits *source // Source FITS file (for the original data) 342 343 ) 343 344 { 344 345 const psMetadata *camera = fpa->camera; // Camera configuration for FPA 345 bool mdok = false; // Status of MD lookup346 bool mdok = false; // Status of MD lookup 346 347 347 348 // Get the required information from the camera configuration 348 349 psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data 349 350 if (! mdok || ! supps) { 350 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n");351 return false;351 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); 352 return false; 352 353 } 353 354 psString sourceType = psMetadataLookupString(&mdok, supps, "MASK.SOURCE"); // Type of source: EXT | FILE 354 355 if (! mdok || strlen(sourceType) <= 0) { 355 psError(PS_ERR_IO, false, "Unable to find MASK.SOURCE in SUPPLEMENTARY section of camera "356 "configuration!\n");357 return false;356 psError(PS_ERR_IO, false, "Unable to find MASK.SOURCE in SUPPLEMENTARY section of camera " 357 "configuration!\n"); 358 return false; 358 359 } 359 360 psString name = psMetadataLookupString(&mdok, supps, "MASK.NAME"); // Name of mask 360 361 if (! mdok || strlen(sourceType) <= 0) { 361 psError(PS_ERR_IO, false, "Unable to find MASK.NAME in SUPPLEMENTARY section of camera "362 "configuration!\n");363 return false;362 psError(PS_ERR_IO, false, "Unable to find MASK.NAME in SUPPLEMENTARY section of camera " 363 "configuration!\n"); 364 return false; 364 365 } 365 366 366 367 // Go through the FPA to each cell/readout to get the mask 367 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the mask368 psArray *chips = fpa->chips; // Array of chips368 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the mask 369 psArray *chips = fpa->chips; // Array of chips 369 370 for (int chipNum = 0; chipNum < chips->n; chipNum++) { 370 pmChip *chip = chips->data[chipNum]; // The current chip of interest 371 if (chip->valid) { 372 if (chip->hdu) { 373 hdu = chip->hdu; 374 } 375 psArray *cells = chip->cells; // Array of cells 376 for (int cellNum = 0; cellNum < cells->n; cellNum++) { 377 pmCell *cell = cells->data[cellNum]; // The current cell of interest 378 if (cell->valid) { 379 if (cell->hdu) { 380 hdu = cell->hdu; 381 } 382 383 // Now, need to find out where to get the pixels 384 psFits *maskSource = source; // Source of mask image 385 if (strcasecmp(sourceType, "FILE") == 0) { 386 // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt" 387 psString filenameExt = p_pmFPATranslateName(name, cell); 388 char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn 389 psString filename = NULL; // The filename 390 psString extname = NULL;// The extenstion name 391 if (colon) { 392 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); 393 if (strlen(colon) > 1) { 394 extname = psStringCopy(colon + 1); 395 } 396 } else { 397 filename = psMemIncrRefCounter(filenameExt); 398 } 399 400 maskSource = psFitsAlloc(filename); 401 if (extname) { 402 if (! psFitsMoveExtName(maskSource, extname)) { 403 psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read mask.\n", 404 extname); 405 return false; 406 } 407 } 408 psFree(filename); 409 psFree(extname); 410 psFree(filenameExt); 411 } else if (strncasecmp(sourceType, "EXT", 3) == 0) { 412 // Source is an extension in the original file 413 psString extname = p_pmFPATranslateName(name, cell); 414 psFitsMoveExtName(maskSource, extname); 415 } 416 417 // We've arrived where the pixels are. Now we need to read them in. 418 psMetadata *header = psFitsReadHeader(NULL, maskSource); // The header 419 bool mdStatus = false; 420 int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); 421 if (!mdStatus) { 422 psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header for " 423 "mask (%s)!\n", name); 424 } 425 if (nAxis != 2 && nAxis != 3) { 426 psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into " 427 "a single image anyway.\n"); 428 } 429 430 int numPlanes = 1; // Number of planes 431 if (nAxis == 3) { 432 numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); 433 if (!mdStatus) { 434 psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); 435 // Try to proceed by taking only the first plane 436 numPlanes = 1; 437 } 438 if (numPlanes != 1 && numPlanes != cell->readouts->n) { 439 psError(PS_ERR_IO, false, "Number of masks (%d) does not match number of " 440 "readouts (%d)\n", numPlanes, cell->readouts->n); 441 // Try to proceed by taking only the first plane 442 numPlanes = 1; 443 } 444 } 445 446 hdu->masks = psArrayAlloc(hdu->images->n); 447 448 // Read each plane into the array 449 psArray *readouts = cell->readouts; // The array of readouts 450 for (int i = 0; i < hdu->masks->n; i++) { 451 psImage *mask = NULL; // The mask to be added 452 if (i < numPlanes) { 453 // Read the mask from the file 454 psTrace(__func__, 9, "Reading plane %d\n", i); 455 psRegion region = {0, 0, 0, 0}; 456 mask = psFitsReadImage(NULL, maskSource, region, i); 457 } else { 458 // One mask in the file is provided for all planes in the original image 459 psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i); 460 psImage *original = hdu->masks->data[0]; 461 mask = psImageCopy(NULL, original, original->type.type); 462 } 463 hdu->masks->data[0] = mask; 464 pmReadout *readout = readouts->data[i]; 465 readout->mask = mask; 466 // Check the dimensions 467 // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size? 468 if (mask->numCols < readout->image->numCols || 469 mask->numRows < readout->image->numRows) 470 { 471 psError(PS_ERR_IO, false, "Mask size (%dx%d) not compatible with image (%dx%d)\n", 472 mask->numCols, mask->numRows, readout->image->numCols, 473 readout->image->numRows); 474 return false; 475 } 476 } // Iterating over readouts 477 } // Valid cells 478 } // Iterating over cells 479 } // Valid chips 371 pmChip *chip = chips->data[chipNum]; // The current chip of interest 372 if (chip->valid) { 373 if (chip->hdu) { 374 hdu = chip->hdu; 375 } 376 psArray *cells = chip->cells; // Array of cells 377 for (int cellNum = 0; cellNum < cells->n; cellNum++) { 378 pmCell *cell = cells->data[cellNum]; // The current cell of interest 379 if (cell->valid) { 380 if (cell->hdu) { 381 hdu = cell->hdu; 382 } 383 384 // Now, need to find out where to get the pixels 385 psFits *maskSource = psMemIncrRefCounter(source); // Source of mask image 386 if (strcasecmp(sourceType, "FILE") == 0) { 387 // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt" 388 psString filenameExt = p_pmFPATranslateName(name, cell); 389 char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn 390 psString filename = NULL; // The filename 391 psString extname = NULL;// The extenstion name 392 if (colon) { 393 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); 394 if (strlen(colon) > 1) { 395 extname = psStringCopy(colon + 1); 396 } 397 } else { 398 filename = psMemIncrRefCounter(filenameExt); 399 } 400 401 psFree(maskSource); 402 maskSource = psFitsOpen(filename, "r"); 403 if (extname) { 404 if (! psFitsMoveExtName(maskSource, extname)) { 405 psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read mask.\n", 406 extname); 407 return false; 408 } 409 } 410 psFree(filename); 411 psFree(extname); 412 psFree(filenameExt); 413 } else if (strncasecmp(sourceType, "EXT", 3) == 0) { 414 // Source is an extension in the original file 415 psString extname = p_pmFPATranslateName(name, cell); 416 psFitsMoveExtName(maskSource, extname); 417 } 418 419 // We've arrived where the pixels are. Now we need to read them in. 420 psMetadata *header = psFitsReadHeader(NULL, maskSource); // The header 421 bool mdStatus = false; 422 int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); 423 if (!mdStatus) { 424 psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header for " 425 "mask (%s)!\n", name); 426 } 427 if (nAxis != 2 && nAxis != 3) { 428 psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into " 429 "a single image anyway.\n"); 430 } 431 432 int numPlanes = 1; // Number of planes 433 if (nAxis == 3) { 434 numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); 435 if (!mdStatus) { 436 psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); 437 // Try to proceed by taking only the first plane 438 numPlanes = 1; 439 } 440 if (numPlanes != 1 && numPlanes != cell->readouts->n) { 441 psError(PS_ERR_IO, false, "Number of masks (%d) does not match number of " 442 "readouts (%d)\n", numPlanes, cell->readouts->n); 443 // Try to proceed by taking only the first plane 444 numPlanes = 1; 445 } 446 } 447 448 hdu->masks = psArrayAlloc(hdu->images->n); 449 450 // Read each plane into the array 451 psArray *readouts = cell->readouts; // The array of readouts 452 for (int i = 0; i < hdu->masks->n; i++) { 453 psImage *mask = NULL; // The mask to be added 454 if (i < numPlanes) { 455 // Read the mask from the file 456 psTrace(__func__, 9, "Reading plane %d\n", i); 457 psRegion region = {0, 0, 0, 0}; 458 mask = psFitsReadImage(NULL, maskSource, region, i); 459 } else { 460 // One mask in the file is provided for all planes in the original image 461 psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i); 462 psImage *original = hdu->masks->data[0]; 463 mask = psImageCopy(NULL, original, original->type.type); 464 } 465 hdu->masks->data[0] = mask; 466 pmReadout *readout = readouts->data[i]; 467 readout->mask = mask; 468 // Check the dimensions 469 // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size? 470 if (mask->numCols < readout->image->numCols || 471 mask->numRows < readout->image->numRows) 472 { 473 psError(PS_ERR_IO, false, "Mask size (%dx%d) not compatible with image (%dx%d)\n", 474 mask->numCols, mask->numRows, readout->image->numCols, 475 readout->image->numRows); 476 return false; 477 } 478 } // Iterating over readouts 479 psFree(maskSource); 480 } // Valid cells 481 } // Iterating over cells 482 } // Valid chips 480 483 } // Iterating over chips 481 484 … … 486 489 // Read a mask into the FPA 487 490 // This is just a copy of the above pmFPAReadMask, replacing "mask" with "weight" throughout. 488 bool pmFPAReadWeight(pmFPA *fpa, // FPA to read into489 psFits *source// Source FITS file (for the original data)491 bool pmFPAReadWeight(pmFPA *fpa, // FPA to read into 492 psFits *source // Source FITS file (for the original data) 490 493 ) 491 494 { 492 495 const psMetadata *camera = fpa->camera; // Camera configuration for FPA 493 bool mdok = false; // Status of MD lookup496 bool mdok = false; // Status of MD lookup 494 497 495 498 // Get the required information from the camera configuration 496 499 psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data 497 500 if (! mdok || ! supps) { 498 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n");499 return false;501 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); 502 return false; 500 503 } 501 504 psString sourceType = psMetadataLookupString(&mdok, supps, "WEIGHT.SOURCE"); // Type of source: EXT | FILE 502 505 if (! mdok || strlen(sourceType) <= 0) { 503 psError(PS_ERR_IO, false, "Unable to find WEIGHT.SOURCE in SUPPLEMENTARY section of camera "504 "configuration!\n");505 return false;506 psError(PS_ERR_IO, false, "Unable to find WEIGHT.SOURCE in SUPPLEMENTARY section of camera " 507 "configuration!\n"); 508 return false; 506 509 } 507 510 psString name = psMetadataLookupString(&mdok, supps, "WEIGHT.NAME"); // Name of weight 508 511 if (! mdok || strlen(sourceType) <= 0) { 509 psError(PS_ERR_IO, false, "Unable to find WEIGHT.NAME in SUPPLEMENTARY section of camera "510 "configuration!\n");511 return false;512 psError(PS_ERR_IO, false, "Unable to find WEIGHT.NAME in SUPPLEMENTARY section of camera " 513 "configuration!\n"); 514 return false; 512 515 } 513 516 514 517 // Go through the FPA to each cell/readout to get the weight 515 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the weight516 psArray *chips = fpa->chips; // Array of chips518 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the weight 519 psArray *chips = fpa->chips; // Array of chips 517 520 for (int chipNum = 0; chipNum < chips->n; chipNum++) { 518 pmChip *chip = chips->data[chipNum]; // The current chip of interest 519 if (chip->valid) { 520 if (chip->hdu) { 521 hdu = chip->hdu; 522 } 523 psArray *cells = chip->cells; // Array of cells 524 for (int cellNum = 0; cellNum < cells->n; cellNum++) { 525 pmCell *cell = cells->data[cellNum]; // The current cell of interest 526 if (cell->valid) { 527 if (cell->hdu) { 528 hdu = cell->hdu; 529 } 530 531 // Now, need to find out where to get the pixels 532 psFits *weightSource = source; // Source of weight image 533 if (strcasecmp(sourceType, "FILE") == 0) { 534 // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt" 535 psString filenameExt = p_pmFPATranslateName(name, cell); 536 char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn 537 psString filename = NULL; // The filename 538 psString extname = NULL;// The extenstion name 539 if (colon) { 540 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); 541 if (strlen(colon) > 1) { 542 extname = psStringCopy(colon + 1); 543 } 544 } else { 545 filename = psMemIncrRefCounter(filenameExt); 546 } 547 548 weightSource = psFitsAlloc(filename); 549 if (extname) { 550 if (! psFitsMoveExtName(weightSource, extname)) { 551 psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read " 552 "weight.\n", extname); 553 return false; 554 } 555 } 556 psFree(filename); 557 psFree(extname); 558 psFree(filenameExt); 559 } else if (strncasecmp(sourceType, "EXT", 3) == 0) { 560 // Source is an extension in the original file 561 psString extname = p_pmFPATranslateName(name, cell); 562 psFitsMoveExtName(weightSource, extname); 563 } 564 565 // We've arrived where the pixels are. Now we need to read them in. 566 psMetadata *header = psFitsReadHeader(NULL, weightSource); // The header 567 bool mdStatus = false; 568 int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); 569 if (!mdStatus) { 570 psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header for " 571 "weight (%s)!\n", name); 572 } 573 if (nAxis != 2 && nAxis != 3) { 574 psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into " 575 "a single image anyway.\n"); 576 } 577 578 int numPlanes = 1; // Number of planes 579 if (nAxis == 3) { 580 numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); 581 if (!mdStatus) { 582 psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); 583 // Try to proceed by taking only the first plane 584 numPlanes = 1; 585 } 586 if (numPlanes != 1 && numPlanes != cell->readouts->n) { 587 psError(PS_ERR_IO, false, "Number of weights (%d) does not match number of " 588 "readouts (%d)\n", numPlanes, cell->readouts->n); 589 // Try to proceed by taking only the first plane 590 numPlanes = 1; 591 } 592 } 593 594 hdu->weights = psArrayAlloc(hdu->images->n); 595 596 // Read each plane into the array 597 psArray *readouts = cell->readouts; // The array of readouts 598 for (int i = 0; i < hdu->weights->n; i++) { 599 psImage *weight = NULL; // The weight to be added 600 if (i < numPlanes) { 601 // Read the weight from the file 602 psTrace(__func__, 9, "Reading plane %d\n", i); 603 psRegion region = {0, 0, 0, 0}; 604 weight = psFitsReadImage(NULL, weightSource, region, i); 605 } else { 606 // One weight in the file is provided for all planes in the original image 607 psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i); 608 psImage *original = hdu->weights->data[0]; 609 weight = psImageCopy(NULL, original, original->type.type); 610 } 611 hdu->weights->data[0] = weight; 612 pmReadout *readout = readouts->data[i]; 613 readout->weight = weight; 614 // Check the dimensions 615 // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size? 616 if (weight->numCols < readout->image->numCols || 617 weight->numRows < readout->image->numRows) 618 { 619 psError(PS_ERR_IO, false, "Weight size (%dx%d) not compatible with image " 620 "(%dx%d)\n", weight->numCols, weight->numRows, readout->image->numCols, 621 readout->image->numRows); 622 return false; 623 } 624 } // Iterating over readouts 625 } // Valid cells 626 } // Iterating over cells 627 } // Valid chips 521 pmChip *chip = chips->data[chipNum]; // The current chip of interest 522 if (chip->valid) { 523 if (chip->hdu) { 524 hdu = chip->hdu; 525 } 526 psArray *cells = chip->cells; // Array of cells 527 for (int cellNum = 0; cellNum < cells->n; cellNum++) { 528 pmCell *cell = cells->data[cellNum]; // The current cell of interest 529 if (cell->valid) { 530 if (cell->hdu) { 531 hdu = cell->hdu; 532 } 533 534 // Now, need to find out where to get the pixels 535 psFits *weightSource = psMemIncrRefCounter(source); // Source of weight image 536 if (strcasecmp(sourceType, "FILE") == 0) { 537 // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt" 538 psString filenameExt = p_pmFPATranslateName(name, cell); 539 char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn 540 psString filename = NULL; // The filename 541 psString extname = NULL;// The extenstion name 542 if (colon) { 543 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); 544 if (strlen(colon) > 1) { 545 extname = psStringCopy(colon + 1); 546 } 547 } else { 548 filename = psMemIncrRefCounter(filenameExt); 549 } 550 551 psFree(weightSource); 552 weightSource = psFitsOpen(filename, "r"); 553 if (extname) { 554 if (! psFitsMoveExtName(weightSource, extname)) { 555 psLogMsg(__func__, PS_LOG_WARN, "Unable to find extension %s to read " 556 "weight.\n", extname); 557 return false; 558 } 559 } 560 psFree(filename); 561 psFree(extname); 562 psFree(filenameExt); 563 } else if (strncasecmp(sourceType, "EXT", 3) == 0) { 564 // Source is an extension in the original file 565 psString extname = p_pmFPATranslateName(name, cell); 566 psFitsMoveExtName(weightSource, extname); 567 } 568 569 // We've arrived where the pixels are. Now we need to read them in. 570 psMetadata *header = psFitsReadHeader(NULL, weightSource); // The header 571 bool mdStatus = false; 572 int nAxis = psMetadataLookupS32(&mdStatus, header, "NAXIS"); 573 if (!mdStatus) { 574 psLogMsg(__func__, PS_LOG_WARN, "There is no NAXIS keyword in the FITS header for " 575 "weight (%s)!\n", name); 576 } 577 if (nAxis != 2 && nAxis != 3) { 578 psLogMsg(__func__, PS_LOG_WARN, "Image is not 2- or 3-dimensional --- reading into " 579 "a single image anyway.\n"); 580 } 581 582 int numPlanes = 1; // Number of planes 583 if (nAxis == 3) { 584 numPlanes = psMetadataLookupS32(&mdStatus, header, "NAXIS3"); 585 if (!mdStatus) { 586 psError(PS_ERR_IO, false, "Unable to read NAXIS3 for 3-dimensional image!\n"); 587 // Try to proceed by taking only the first plane 588 numPlanes = 1; 589 } 590 if (numPlanes != 1 && numPlanes != cell->readouts->n) { 591 psError(PS_ERR_IO, false, "Number of weights (%d) does not match number of " 592 "readouts (%d)\n", numPlanes, cell->readouts->n); 593 // Try to proceed by taking only the first plane 594 numPlanes = 1; 595 } 596 } 597 598 hdu->weights = psArrayAlloc(hdu->images->n); 599 600 // Read each plane into the array 601 psArray *readouts = cell->readouts; // The array of readouts 602 for (int i = 0; i < hdu->weights->n; i++) { 603 psImage *weight = NULL; // The weight to be added 604 if (i < numPlanes) { 605 // Read the weight from the file 606 psTrace(__func__, 9, "Reading plane %d\n", i); 607 psRegion region = {0, 0, 0, 0}; 608 weight = psFitsReadImage(NULL, weightSource, region, i); 609 } else { 610 // One weight in the file is provided for all planes in the original image 611 psTrace(__func__, 9, "Copying plane 0 into plane %d\n", i); 612 psImage *original = hdu->weights->data[0]; 613 weight = psImageCopy(NULL, original, original->type.type); 614 } 615 hdu->weights->data[0] = weight; 616 pmReadout *readout = readouts->data[i]; 617 readout->weight = weight; 618 // Check the dimensions 619 // XXX: Perhaps this check should be against the CELL.TRIMSEC, not the image size? 620 if (weight->numCols < readout->image->numCols || 621 weight->numRows < readout->image->numRows) 622 { 623 psError(PS_ERR_IO, false, "Weight size (%dx%d) not compatible with image " 624 "(%dx%d)\n", weight->numCols, weight->numRows, readout->image->numCols, 625 readout->image->numRows); 626 return false; 627 } 628 } // Iterating over readouts 629 psFree(weightSource); 630 } // Valid cells 631 } // Iterating over cells 632 } // Valid chips 628 633 } // Iterating over chips 629 634 -
trunk/archive/scripts/src/phase2/pmFPAWrite.c
r5633 r5786 7 7 #include "pmFPARead.h" 8 8 9 static bool writeHDU(psFits *fits, // FITS file to which to write10 p_pmHDU *hdu// Pixel data to write11 ) 12 { 13 bool status = true; // Status of write, to return9 static bool writeHDU(psFits *fits, // FITS file to which to write 10 p_pmHDU *hdu // Pixel data to write 11 ) 12 { 13 bool status = true; // Status of write, to return 14 14 for (int i = 0; i < hdu->images->n; i++) { 15 status &= psFitsWriteImage(fits, hdu->header, hdu->images->data[i], i);16 // XXX: Insert here the writing on mask and weight images15 status &= psFitsWriteImage(fits, hdu->header, hdu->images->data[i], i); 16 // XXX: Insert here the writing on mask and weight images 17 17 } 18 18 … … 21 21 22 22 23 bool pmFPAWrite(psFits *fits, // FITS file to which to write 24 pmFPA *fpa, // FPA to write 25 psDB *db // Database to update 26 ) 27 { 28 bool status = true; // Status of writing, to return 29 23 bool pmFPAWrite(psFits *fits, // FITS file to which to write 24 pmFPA *fpa, // FPA to write 25 psDB *db // Database to update 26 ) 27 { 28 bool status = true; // Status of writing, to return 29 30 psTrace(__func__, 7, "Outgesting FPA concepts...\n"); 30 31 pmFPAOutgestConcepts(fpa, db); 31 32 32 33 // Write the primary header 33 34 if (fpa->phu) { 34 status &= psFitsMoveExtNum(fits, 0, false);35 status &= psFitsWriteHeader(fpa->phu, fits);36 } 37 38 psArray *chips = fpa->chips; // Array of component chips35 status &= psFitsMoveExtNum(fits, 0, false); 36 status &= psFitsWriteHeader(fpa->phu, fits); 37 } 38 39 psArray *chips = fpa->chips; // Array of component chips 39 40 for (int i = 0; i < chips->n; i++) { 40 pmChip *chip = chips->data[i]; // The component chip 41 if (chip->valid) { 42 pmChipOutgestConcepts(chip, db); 43 44 psArray *cells = chip->cells; // Array of component cells 45 for (int j = 0; j < cells->n; j++) { 46 pmCell *cell = cells->data[j]; // The component cell 47 if (cell->valid) { 48 pmCellOutgestConcepts(cell, db); 49 50 if (cell->hdu && strlen(cell->hdu->extname) > 0) { 51 status &= writeHDU(fits, cell->hdu); 52 } 53 } 54 } 55 56 if (chip->hdu && strlen(chip->hdu->extname) > 0) { 57 status &= writeHDU(fits, chip->hdu); 58 } 59 } 41 pmChip *chip = chips->data[i]; // The component chip 42 if (chip->valid) { 43 psTrace(__func__, 1, "Writing out chip %d...\n", i); 44 45 pmChipOutgestConcepts(chip, db); 46 47 psArray *cells = chip->cells; // Array of component cells 48 for (int j = 0; j < cells->n; j++) { 49 pmCell *cell = cells->data[j]; // The component cell 50 if (cell->valid) { 51 psTrace(__func__, 2, "Writing out cell, %d...\n", j); 52 53 pmCellOutgestConcepts(cell, db); 54 55 if (cell->hdu && strlen(cell->hdu->extname) > 0) { 56 status &= writeHDU(fits, cell->hdu); 57 } 58 } 59 } 60 61 if (chip->hdu && strlen(chip->hdu->extname) > 0) { 62 status &= writeHDU(fits, chip->hdu); 63 } 64 } 60 65 61 66 } 62 67 63 68 if (fpa->hdu && strlen(fpa->hdu->extname) > 0) { 64 status &= writeHDU(fits, fpa->hdu);69 status &= writeHDU(fits, fpa->hdu); 65 70 } 66 71 … … 71 76 72 77 73 bool pmFPAWriteMask(pmFPA *fpa, // FPA containing mask to write74 psFits *fits// FITS file for image78 bool pmFPAWriteMask(pmFPA *fpa, // FPA containing mask to write 79 psFits *fits // FITS file for image 75 80 ) 76 81 { 77 82 const psMetadata *camera = fpa->camera; // Camera configuration for FPA 78 bool mdok = false; // Status of MD lookup83 bool mdok = false; // Status of MD lookup 79 84 80 85 // Get the required information from the camera configuration 81 86 psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data 82 87 if (! mdok || ! supps) { 83 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n");84 return false;88 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); 89 return false; 85 90 } 86 91 psString sourceType = psMetadataLookupString(&mdok, supps, "MASK.SOURCE"); // Type of source: EXT | FILE 87 92 if (! mdok || strlen(sourceType) <= 0) { 88 psError(PS_ERR_IO, false, "Unable to find MASK.SOURCE in SUPPLEMENTARY section of camera "89 "configuration!\n");90 return false;93 psError(PS_ERR_IO, false, "Unable to find MASK.SOURCE in SUPPLEMENTARY section of camera " 94 "configuration!\n"); 95 return false; 91 96 } 92 97 psString name = psMetadataLookupString(&mdok, supps, "MASK.NAME"); // Name of mask 93 98 if (! mdok || strlen(sourceType) <= 0) { 94 psError(PS_ERR_IO, false, "Unable to find MASK.NAME in SUPPLEMENTARY section of camera "95 "configuration!\n");96 return false;99 psError(PS_ERR_IO, false, "Unable to find MASK.NAME in SUPPLEMENTARY section of camera " 100 "configuration!\n"); 101 return false; 97 102 } 98 103 99 104 // Go through the FPA to each cell/readout to get the mask 100 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the mask101 psArray *chips = fpa->chips; // Array of chips105 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the mask 106 psArray *chips = fpa->chips; // Array of chips 102 107 for (int chipNum = 0; chipNum < chips->n; chipNum++) { 103 pmChip *chip = chips->data[chipNum]; // The current chip of interest104 if (chip->valid) {105 if (chip->hdu) {106 hdu = chip->hdu;107 }108 psArray *cells = chip->cells;// Array of cells109 for (int cellNum = 0; cellNum < cells->n; cellNum++) {110 pmCell *cell = cells->data[cellNum]; // The current cell of interest111 if (cell->valid) {112 if (cell->hdu) {113 hdu = cell->hdu;114 }115 116 // Now, need to find out where to write the pixels117 psFits *maskDest = psMemIncrRefCounter(fits); // Destination of mask image118 psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name119 if (strcasecmp(sourceType, "FILE") == 0) {120 // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt"121 psString filenameExt = p_pmFPATranslateName(name, cell);122 char *colon = strchr(filenameExt, ':');// Pointer to a colon in the filename-extn123 psString filename = NULL; // The filename124 psString extname = NULL;// The extenstion name125 if (colon) {126 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon));127 if (strlen(colon) > 1) {128 extname = psStringCopy(colon + 1);129 }130 } else {131 filename = psMemIncrRefCounter(filenameExt);132 }133 134 psFree(maskDest);135 maskDest = psFitsAlloc(filename);136 if (extname) {137 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname);138 }139 psFree(filename);140 psFree(extname);141 psFree(filenameExt);142 } else if (strncasecmp(sourceType, "EXT", 3) == 0) {143 // Source is an extension in the original file144 psString extname = p_pmFPATranslateName(name, cell);145 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname);146 psFree(extname);147 }148 149 // We've arrived where the pixels are. Now we need to write them out.150 psArray *readouts = cell->readouts; // The array of readouts151 for (int readNum = 0; readNum < readouts->n; readNum++) {152 pmReadout *readout = readouts->data[readNum]; // The readout of interest153 if (! readout->mask) {154 psLogMsg(__func__, PS_LOG_WARN, "No mask to write out in %d,%d,%d\n",155 chipNum, cellNum, readNum);156 } else {157 // XXX: Need to add the extname to the existing header158 if (! psFitsWriteImage(maskDest, header, readout->mask, readNum)) {159 psError(PS_ERR_IO, false, "Unable to write mask plane %d in extension %s\n",160 readNum, hdu->extname);161 return false;162 }163 }164 } // Iterating over readouts165 psFree(header);166 psFree(maskDest);167 } // Valid cells168 } // Iterating over cells169 } // Valid chips108 pmChip *chip = chips->data[chipNum]; // The current chip of interest 109 if (chip->valid) { 110 if (chip->hdu) { 111 hdu = chip->hdu; 112 } 113 psArray *cells = chip->cells; // Array of cells 114 for (int cellNum = 0; cellNum < cells->n; cellNum++) { 115 pmCell *cell = cells->data[cellNum]; // The current cell of interest 116 if (cell->valid) { 117 if (cell->hdu) { 118 hdu = cell->hdu; 119 } 120 121 // Now, need to find out where to write the pixels 122 psFits *maskDest = psMemIncrRefCounter(fits); // Destination of mask image 123 psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name 124 if (strcasecmp(sourceType, "FILE") == 0) { 125 // Source is a file (with optional extension, e.g., "myMaskFile.fits:thisExt" 126 psString filenameExt = p_pmFPATranslateName(name, cell); 127 char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn 128 psString filename = NULL; // The filename 129 psString extname = NULL;// The extenstion name 130 if (colon) { 131 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); 132 if (strlen(colon) > 1) { 133 extname = psStringCopy(colon + 1); 134 } 135 } else { 136 filename = psMemIncrRefCounter(filenameExt); 137 } 138 139 psFree(maskDest); 140 maskDest = psFitsOpen(filename, "rw"); 141 if (extname) { 142 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); 143 } 144 psFree(filename); 145 psFree(extname); 146 psFree(filenameExt); 147 } else if (strncasecmp(sourceType, "EXT", 3) == 0) { 148 // Source is an extension in the original file 149 psString extname = p_pmFPATranslateName(name, cell); 150 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); 151 psFree(extname); 152 } 153 154 // We've arrived where the pixels are. Now we need to write them out. 155 psArray *readouts = cell->readouts; // The array of readouts 156 for (int readNum = 0; readNum < readouts->n; readNum++) { 157 pmReadout *readout = readouts->data[readNum]; // The readout of interest 158 if (! readout->mask) { 159 psLogMsg(__func__, PS_LOG_WARN, "No mask to write out in %d,%d,%d\n", 160 chipNum, cellNum, readNum); 161 } else { 162 // XXX: Need to add the extname to the existing header 163 if (! psFitsWriteImage(maskDest, header, readout->mask, readNum)) { 164 psError(PS_ERR_IO, false, "Unable to write mask plane %d in extension %s\n", 165 readNum, hdu->extname); 166 return false; 167 } 168 } 169 } // Iterating over readouts 170 psFree(header); 171 psFree(maskDest); 172 } // Valid cells 173 } // Iterating over cells 174 } // Valid chips 170 175 } // Iterating over chips 171 176 … … 174 179 175 180 176 bool pmFPAWriteWeight(pmFPA *fpa, // FPA containing mask to write177 psFits *fits// FITS file for image181 bool pmFPAWriteWeight(pmFPA *fpa, // FPA containing mask to write 182 psFits *fits // FITS file for image 178 183 ) 179 184 { 180 185 const psMetadata *camera = fpa->camera; // Camera configuration for FPA 181 bool mdok = false; // Status of MD lookup186 bool mdok = false; // Status of MD lookup 182 187 183 188 // Get the required information from the camera configuration 184 189 psMetadata *supps = psMetadataLookupMD(&mdok, camera, "SUPPLEMENTARY"); // Rules for supplementary data 185 190 if (! mdok || ! supps) { 186 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n");187 return false;191 psError(PS_ERR_IO, false, "Unable to find SUPPLEMENTARY in camera configuration!\n"); 192 return false; 188 193 } 189 194 psString sourceType = psMetadataLookupString(&mdok, supps, "WEIGHT.SOURCE"); // Type of source: EXT | FILE 190 195 if (! mdok || strlen(sourceType) <= 0) { 191 psError(PS_ERR_IO, false, "Unable to find WEIGHT.SOURCE in SUPPLEMENTARY section of camera "192 "configuration!\n");193 return false;196 psError(PS_ERR_IO, false, "Unable to find WEIGHT.SOURCE in SUPPLEMENTARY section of camera " 197 "configuration!\n"); 198 return false; 194 199 } 195 200 psString name = psMetadataLookupString(&mdok, supps, "WEIGHT.NAME"); // Name of weight 196 201 if (! mdok || strlen(sourceType) <= 0) { 197 psError(PS_ERR_IO, false, "Unable to find WEIGHT.NAME in SUPPLEMENTARY section of camera "198 "configuration!\n");199 return false;202 psError(PS_ERR_IO, false, "Unable to find WEIGHT.NAME in SUPPLEMENTARY section of camera " 203 "configuration!\n"); 204 return false; 200 205 } 201 206 202 207 // Go through the FPA to each cell/readout to get the weight 203 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the weight204 psArray *chips = fpa->chips; // Array of chips208 p_pmHDU *hdu = fpa->hdu; // The HDU into which we will read the weight 209 psArray *chips = fpa->chips; // Array of chips 205 210 for (int chipNum = 0; chipNum < chips->n; chipNum++) { 206 pmChip *chip = chips->data[chipNum]; // The current chip of interest207 if (chip->valid) {208 if (chip->hdu) {209 hdu = chip->hdu;210 }211 psArray *cells = chip->cells;// Array of cells212 for (int cellNum = 0; cellNum < cells->n; cellNum++) {213 pmCell *cell = cells->data[cellNum]; // The current cell of interest214 if (cell->valid) {215 if (cell->hdu) {216 hdu = cell->hdu;217 }218 219 // Now, need to find out where to write the pixels220 psFits *weightDest = psMemIncrRefCounter(fits); // Destination of weight image221 psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name222 if (strcasecmp(sourceType, "FILE") == 0) {223 // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt"224 psString filenameExt = p_pmFPATranslateName(name, cell);225 char *colon = strchr(filenameExt, ':');// Pointer to a colon in the filename-extn226 psString filename = NULL; // The filename227 psString extname = NULL;// The extenstion name228 if (colon) {229 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon));230 if (strlen(colon) > 1) {231 extname = psStringCopy(colon + 1);232 }233 } else {234 filename = psMemIncrRefCounter(filenameExt);235 }236 237 psFree(weightDest);238 weightDest = psFitsAlloc(filename);239 if (extname) {240 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname);241 }242 psFree(filename);243 psFree(extname);244 psFree(filenameExt);245 } else if (strncasecmp(sourceType, "EXT", 3) == 0) {246 // Source is an extension in the original file247 psString extname = p_pmFPATranslateName(name, cell);248 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname);249 psFree(extname);250 }251 252 // We've arrived where the pixels are. Now we need to write them out.253 psArray *readouts = cell->readouts; // The array of readouts254 for (int readNum = 0; readNum < readouts->n; readNum++) {255 pmReadout *readout = readouts->data[readNum]; // The readout of interest256 if (! readout->weight) {257 psLogMsg(__func__, PS_LOG_WARN, "No weight image to write out in %d,%d,%d\n",258 chipNum, cellNum, readNum);259 } else {260 if (! psFitsWriteImage(weightDest, header, readout->weight, readNum)) {261 psError(PS_ERR_IO, false, "Unable to write weight plane %d in extension %s\n",262 readNum, hdu->extname);263 return false;264 }265 }266 } // Iterating over readouts267 psFree(header);268 psFree(weightDest);269 } // Valid cells270 } // Iterating over cells271 } // Valid chips211 pmChip *chip = chips->data[chipNum]; // The current chip of interest 212 if (chip->valid) { 213 if (chip->hdu) { 214 hdu = chip->hdu; 215 } 216 psArray *cells = chip->cells; // Array of cells 217 for (int cellNum = 0; cellNum < cells->n; cellNum++) { 218 pmCell *cell = cells->data[cellNum]; // The current cell of interest 219 if (cell->valid) { 220 if (cell->hdu) { 221 hdu = cell->hdu; 222 } 223 224 // Now, need to find out where to write the pixels 225 psFits *weightDest = psMemIncrRefCounter(fits); // Destination of weight image 226 psMetadata *header = psMetadataAlloc(); // A dummy header, containing the extension name 227 if (strcasecmp(sourceType, "FILE") == 0) { 228 // Source is a file (with optional extension, e.g., "myWeightFile.fits:thisExt" 229 psString filenameExt = p_pmFPATranslateName(name, cell); 230 char *colon = strchr(filenameExt, ':'); // Pointer to a colon in the filename-extn 231 psString filename = NULL; // The filename 232 psString extname = NULL;// The extenstion name 233 if (colon) { 234 filename = psStringNCopy(filenameExt, strlen(filenameExt) - strlen(colon)); 235 if (strlen(colon) > 1) { 236 extname = psStringCopy(colon + 1); 237 } 238 } else { 239 filename = psMemIncrRefCounter(filenameExt); 240 } 241 242 psFree(weightDest); 243 weightDest = psFitsOpen(filename, "rw"); 244 if (extname) { 245 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); 246 } 247 psFree(filename); 248 psFree(extname); 249 psFree(filenameExt); 250 } else if (strncasecmp(sourceType, "EXT", 3) == 0) { 251 // Source is an extension in the original file 252 psString extname = p_pmFPATranslateName(name, cell); 253 psMetadataAddStr(header, PS_LIST_TAIL, "EXTNAME", 0, "Extension name", extname); 254 psFree(extname); 255 } 256 257 // We've arrived where the pixels are. Now we need to write them out. 258 psArray *readouts = cell->readouts; // The array of readouts 259 for (int readNum = 0; readNum < readouts->n; readNum++) { 260 pmReadout *readout = readouts->data[readNum]; // The readout of interest 261 if (! readout->weight) { 262 psLogMsg(__func__, PS_LOG_WARN, "No weight image to write out in %d,%d,%d\n", 263 chipNum, cellNum, readNum); 264 } else { 265 if (! psFitsWriteImage(weightDest, header, readout->weight, readNum)) { 266 psError(PS_ERR_IO, false, "Unable to write weight plane %d in extension %s\n", 267 readNum, hdu->extname); 268 return false; 269 } 270 } 271 } // Iterating over readouts 272 psFree(header); 273 psFree(weightDest); 274 } // Valid cells 275 } // Iterating over cells 276 } // Valid chips 272 277 } // Iterating over chips 273 278 -
trunk/archive/scripts/src/phase2/pmSubtractBias.c
r5462 r5786 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1. 2$ $Name: not supported by cvs2svn $9 * @date $Date: 2005-1 1-03 01:30:32$8 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2005-12-14 03:47:35 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 214 214 psF32 x; 215 215 psS32 i; 216 printf("Got here\n");217 216 if (fit == PM_FIT_POLYNOMIAL) { 218 217 // Fit a polynomial to the old overscan vector. -
trunk/archive/scripts/src/phase2/psAdditionals.c
r5651 r5786 5 5 6 6 7 psMetadata *p sMetadataCopy(psMetadata *out,8 const psMetadata *in)7 psMetadata *pap_psMetadataCopy(psMetadata *out, 8 const psMetadata *in) 9 9 { 10 10 PS_ASSERT_PTR_NON_NULL(in,NULL); … … 20 20 // Need to look for MULTI, which won't be picked up using the iterator. 21 21 psMetadataItem *multiCheckItem = psMetadataLookup(in, inItem->name); 22 int multiFlag = 0; // Flag to indicate MULTI or not22 unsigned int flag = PS_META_REPLACE; // Flag to indicate MULTI; otherwise, replace 23 23 if (multiCheckItem->type == PS_DATA_METADATA_MULTI) { 24 multiFlag = PS_DATA_METADATA_MULTI; 25 } 24 psTrace(__func__, 10, "MULTI: %s (%s)\n", inItem->name, inItem->comment); 25 flag = PS_DATA_METADATA_MULTI; 26 } 27 28 psTrace(__func__, 5, "Copying %s (%s)...\n", inItem->name, inItem->comment); 26 29 27 30 #define PS_METADATA_COPY_CASE(NAME,TYPE) \ 28 31 case PS_TYPE_##NAME: \ 29 psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_TYPE_##NAME | multiFlag, inItem->comment, \ 30 inItem->data.TYPE); \ 32 if (! psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_TYPE_##NAME | flag, inItem->comment, \ 33 inItem->data.TYPE)) { \ 34 psErrorStackPrint(stderr, "Error copying %s (%s) in the metadata\n", inItem->name, \ 35 inItem->comment); \ 36 } \ 31 37 break; 32 38 … … 34 40 // Numerical types 35 41 PS_METADATA_COPY_CASE(BOOL,B); 36 PS_METADATA_COPY_CASE(S8, S8);37 PS_METADATA_COPY_CASE(S16, S16);38 PS_METADATA_COPY_CASE(S32, S32);39 PS_METADATA_COPY_CASE(U8, U8);40 PS_METADATA_COPY_CASE(U16, U16);41 PS_METADATA_COPY_CASE(U32, U32);42 PS_METADATA_COPY_CASE(F32, F32);43 PS_METADATA_COPY_CASE(F64, F64);42 PS_METADATA_COPY_CASE(S8,S8); 43 PS_METADATA_COPY_CASE(S16,S16); 44 PS_METADATA_COPY_CASE(S32,S32); 45 PS_METADATA_COPY_CASE(U8,U8); 46 PS_METADATA_COPY_CASE(U16,U16); 47 PS_METADATA_COPY_CASE(U32,U32); 48 PS_METADATA_COPY_CASE(F32,F32); 49 PS_METADATA_COPY_CASE(F64,F64); 44 50 45 51 // String: relying on the fact that this will copy the string, not point at it. 46 52 case PS_DATA_STRING: 47 psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_STRING | multiFlag, inItem->comment,53 psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_STRING | flag, inItem->comment, 48 54 inItem->data.V); 49 55 break; … … 52 58 case PS_DATA_METADATA: 53 59 { 54 psMetadata *metadata = p sMetadataCopy(NULL, inItem->data.md);55 psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_METADATA | multiFlag, inItem->comment,60 psMetadata *metadata = pap_psMetadataCopy(NULL, inItem->data.md); 61 psMetadataAdd(out, PS_LIST_TAIL, inItem->name, PS_DATA_METADATA | flag, inItem->comment, 56 62 metadata); 57 63 break; … … 60 66 default: 61 67 numPointers++; 62 psMetadataItemAlloc(inItem->name, inItem->type, inItem->comment, inItem->data.V); 68 psTrace(__func__, 10, "Copying a pointer in the metadata: %x\n", inItem->type); 69 psMetadataAdd(out, PS_LIST_TAIL, inItem->name, inItem->type | flag, inItem->comment, 70 inItem->data.V); 63 71 break; 64 72 } … … 68 76 if (numPointers > 0) { 69 77 psLogMsg(__func__, PS_LOG_WARN, "Forced to copy %d pointers when copying metadata. Updating the " 70 "copied psMetadata will affect the original!\n" );78 "copied psMetadata will affect the original!\n", numPointers); 71 79 } 72 80 -
trunk/archive/scripts/src/phase2/psAdditionals.h
r5651 r5786 5 5 6 6 // Deep copy of metadata 7 psMetadata *p sMetadataCopy(psMetadata *out, const psMetadata *in);7 psMetadata *pap_psMetadataCopy(psMetadata *out, const psMetadata *in); 8 8 9 9 // Get a value from the metadata that we believe should be metadata.
Note:
See TracChangeset
for help on using the changeset viewer.
