Index: trunk/archive/scripts/src/phase2/papPhase2.c =================================================================== --- trunk/archive/scripts/src/phase2/papPhase2.c (revision 5651) +++ trunk/archive/scripts/src/phase2/papPhase2.c (revision 5786) @@ -17,5 +17,5 @@ #include "pmSubtractBias.h" #include "pmChipMosaic.h" -#include "pmFPAMorph.h" +//#include "pmFPAMorph.h" // Phase 2 needs to: @@ -39,17 +39,17 @@ // phase2 INPUT.fits OUTPUT.fits -bias BIAS.fits -dark DARK.fits -flat FLAT.fits -chip CHIP -// +// // Most are self-explanatory. "-chip" says "only work on this particular chip". -#define RECIPE "PHASE2" // Name of the recipe to use +#define RECIPE "PHASE2" // Name of the recipe to use static psMemId memPrintAlloc(const psMemBlock *mb) { printf("Allocated memory block %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); return 0; } @@ -57,8 +57,8 @@ { printf("Freed memory block %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); return 0; } @@ -67,8 +67,8 @@ psMemBlock *mb = ((psMemBlock*)ptr) - 1; printf("Memory block %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); } @@ -102,12 +102,12 @@ psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", 0, "Chip number to process (if positive)", -1); if (! psArgumentParse(arguments, &argc, argv) || argc != 3) { - printf("\nPan-STARRS Phase 2 processing\n\n"); - printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]); - psArgumentHelp(arguments); - psFree(arguments); - exit(EXIT_FAILURE); - } - const char *inputName = argv[1]; // Name of input image - const char *outputName = argv[2]; // Name of output image + printf("\nPan-STARRS Phase 2 processing\n\n"); + printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]); + psArgumentHelp(arguments); + psFree(arguments); + exit(EXIT_FAILURE); + } + const char *inputName = argv[1]; // Name of input image + const char *outputName = argv[2]; // Name of output image const char *biasName = psMetadataLookupString(NULL, arguments, "-bias"); // Name of bias image const char *darkName = psMetadataLookupString(NULL, arguments, "-dark"); // Name of dark image @@ -123,42 +123,34 @@ // Open the input -#ifdef PRODUCTION psFits *inputFile = psFitsOpen(inputName, "r"); // File handle for FITS file -#else - psFits *inputFile = psFitsAlloc(inputName); // File handle for FITS file -#endif if (! inputFile) { - psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName); + psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName); exit(EXIT_FAILURE); } psMetadata *header = psFitsReadHeader(NULL, inputFile); // FITS header #if PRODUCTION - psDB *database = pmConfigDB(site); // Database handle + psDB *database = pmConfigDB(site); // Database handle #else - psDB *database = NULL; // Database handle + psDB *database = NULL; // Database handle #endif // Open the output and output mask // We do it here so that we don't process the whole lot and then find out we can't open the output file -#ifdef PRODUCTION psFits *outputFile = psFitsOpen(outputName, "w"); -#else - psFits *outputFile = psFitsAlloc(outputName); -#endif if (! outputFile) { - psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName); - exit(EXIT_FAILURE); + psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName); + exit(EXIT_FAILURE); } // Get camera configuration from header if not already defined if (! camera) { - camera = pmConfigCameraFromHeader(site, header); - if (! camera) { - psErrorStackPrint(stderr, "Can't find camera configuration!\n"); - exit(EXIT_FAILURE); - } + camera = pmConfigCameraFromHeader(site, header); + if (! camera) { + psErrorStackPrint(stderr, "Can't find camera configuration!\n"); + exit(EXIT_FAILURE); + } } else if (! pmConfigValidateCamera(camera, header)) { - psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName); - exit(EXIT_FAILURE); + psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName); + exit(EXIT_FAILURE); } if (! recipe && !(recipe = pmConfigRecipeFromCamera(camera, RECIPE))) { @@ -175,106 +167,106 @@ // Set various tasks - bool doMask = false; // Mask bad pixles - bool doNonLin = false; // Non-linearity correction - bool doBias = false; // Bias subtraction - bool doDark = false; // Dark subtraction - bool doOverscan = false; // Overscan subtraction - bool doFlat = false; // Flat-field normalisation - bool doFringe = false; // Fringe subtraction - bool doSource = false; // Source identification and photometry - bool doAstrom = false; // Astrometry + bool doMask = false; // Mask bad pixles + bool doNonLin = false; // Non-linearity correction + bool doBias = false; // Bias subtraction + bool doDark = false; // Dark subtraction + bool doOverscan = false; // Overscan subtraction + bool doFlat = false; // Flat-field normalisation + bool doFringe = false; // Fringe subtraction + bool doSource = false; // Source identification and photometry + bool doAstrom = false; // Astrometry pmOverscanAxis overscanMode = PM_OVERSCAN_NONE; // Axis for overscan pmFit overscanFitType = PM_FIT_NONE; // Fit type for overscan - int overscanBins = 1; // Number of pixels per bin for overscan - psStats *overscanStats = NULL; // Statistics for overscan - void *overscanFit = NULL; // Overscan fit (polynomial or spline) + int overscanBins = 1; // Number of pixels per bin for overscan + psStats *overscanStats = NULL; // Statistics for overscan + void *overscanFit = NULL; // Overscan fit (polynomial or spline) if (psMetadataLookupBool(NULL, recipe, "MASK")) { - if (strlen(maskName) > 0) { - doMask = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking " - "will be performed.\n"); - } + if (strlen(maskName) > 0) { + doMask = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking " + "will be performed.\n"); + } } if (psMetadataLookupBool(NULL, recipe, "NONLIN")) { - doNonLin = true; + doNonLin = true; } if (psMetadataLookupBool(NULL, recipe, "BIAS")) { - if (strlen(biasName) > 0) { - doBias = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- " - "no bias subtraction will be performed.\n"); - } + if (strlen(biasName) > 0) { + doBias = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- " + "no bias subtraction will be performed.\n"); + } } if (psMetadataLookupBool(NULL, recipe, "DARK")) { - if (strlen(darkName) > 0) { - doDark = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- " - "no dark subtraction will be performed.\n"); - } + if (strlen(darkName) > 0) { + doDark = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- " + "no dark subtraction will be performed.\n"); + } } if (psMetadataLookupBool(NULL, recipe, "OVERSCAN")) { - doOverscan = true; - psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE"); - if (strcasecmp(mode, "INDIVIDUAL") == 0) { - overscanMode = PM_OVERSCAN_ROWS; - // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip - } else if (strcasecmp(mode, "ALL") == 0) { - overscanMode = PM_OVERSCAN_ALL; - } else if (strcasecmp(mode, "NONE") != 0) { - psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:" - " assuming NONE.\n", mode); - } - psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT"); - if (strcasecmp(fit, "POLYNOMIAL") == 0) { - overscanFitType = PM_FIT_POLYNOMIAL; - int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit - overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); - } else if (strcasecmp(fit, "SPLINE") == 0) { - overscanFitType = PM_FIT_SPLINE; - int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit - overscanFit = NULL; -// overscanFit = psSpline1DAlloc(); - } else if (strcasecmp(fit, "NONE") != 0) { - psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:" - " assuming NONE.\n", fit); - } - overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN"); - if (overscanBins <= 0) { - psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins); - overscanBins = 1; - } - psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT"); - if (strcasecmp(stat, "MEAN")) { - overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - } else if (strcasecmp(stat, "MEDIAN")) { - overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); - } else { - psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat); - overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - } + doOverscan = true; + psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE"); + if (strcasecmp(mode, "INDIVIDUAL") == 0) { + overscanMode = PM_OVERSCAN_ROWS; + // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip + } else if (strcasecmp(mode, "ALL") == 0) { + overscanMode = PM_OVERSCAN_ALL; + } else if (strcasecmp(mode, "NONE") != 0) { + psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:" + " assuming NONE.\n", mode); + } + psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT"); + if (strcasecmp(fit, "POLYNOMIAL") == 0) { + overscanFitType = PM_FIT_POLYNOMIAL; + int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit + overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } else if (strcasecmp(fit, "SPLINE") == 0) { + overscanFitType = PM_FIT_SPLINE; + int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit + overscanFit = NULL; +// overscanFit = psSpline1DAlloc(); + } else if (strcasecmp(fit, "NONE") != 0) { + psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:" + " assuming NONE.\n", fit); + } + overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN"); + if (overscanBins <= 0) { + psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins); + overscanBins = 1; + } + psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT"); + if (strcasecmp(stat, "MEAN")) { + overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + } else if (strcasecmp(stat, "MEDIAN")) { + overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + } else { + psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat); + overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + } } if (psMetadataLookupBool(NULL, recipe, "FLAT")) { - if (strlen(flatName) > 0) { - doFlat = true; - } else { - psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- " - "no flat-fielding will be performed.\n"); - } + if (strlen(flatName) > 0) { + doFlat = true; + } else { + psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- " + "no flat-fielding will be performed.\n"); + } } // Chip selection if (chipNum >= 0) { - if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) || - ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) || - ! pmFPASelectChip(mask, chipNum)) { - psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum); + if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) || + ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) || + ! pmFPASelectChip(mask, chipNum)) { + psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum); } @@ -284,16 +276,8 @@ psLogMsg("phase2", PS_LOG_INFO, "Opening input image: %s\n", inputName); if (! pmFPARead(input, inputFile, header, database)) { - psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName); - exit(EXIT_FAILURE); - } - - pmFPAPrint(input); - exit(1); - -#ifdef PRODUCTION + psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName); + exit(EXIT_FAILURE); + } psFitsClose(inputFile); -#else - psFree(inputFile); -#endif #if 0 @@ -302,47 +286,47 @@ //#else { - // Generate mask and weight frame - p_pmHDU *hdu = NULL; - if (input->hdu) { - hdu = input->hdu; - } - psArray *chips = input->chips; - for (int chipNum = 0; chipNum < chips->n; chipNum++) { - pmChip *chip = chips->data[chipNum]; - if (chip->valid) { - if (chip->hdu) { - hdu = chip->hdu; - } - psArray *cells = chip->cells; - for (int cellNum = 0; cellNum < cells->n; cellNum++) { - pmCell *cell = cells->data[cellNum]; - if (cell->valid) { - if (cell->hdu) { - hdu = cell->hdu; - } - - hdu->masks = psArrayAlloc(hdu->images->n); - hdu->weights = psArrayAlloc(hdu->images->n); - psArray *readouts = cell->readouts; - for (int readNum = 0; readNum < hdu->images->n; readNum++) { - psImage *image = hdu->images->data[readNum]; - psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); - psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); - pmReadout *readout = readouts->data[readNum]; - for (int j = 0; j < image->numRows; j++) { - for (int i = 0; i < image->numCols; i++) { - mask->data.U8[j][i] = j + i; - weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]); - } - } - hdu->masks->data[readNum] = mask; - hdu->weights->data[readNum] = weight; - readout->mask = psMemIncrRefCounter(mask); - readout->weight = psMemIncrRefCounter(weight); - } - } - } - } - } + // Generate mask and weight frame + p_pmHDU *hdu = NULL; + if (input->hdu) { + hdu = input->hdu; + } + psArray *chips = input->chips; + for (int chipNum = 0; chipNum < chips->n; chipNum++) { + pmChip *chip = chips->data[chipNum]; + if (chip->valid) { + if (chip->hdu) { + hdu = chip->hdu; + } + psArray *cells = chip->cells; + for (int cellNum = 0; cellNum < cells->n; cellNum++) { + pmCell *cell = cells->data[cellNum]; + if (cell->valid) { + if (cell->hdu) { + hdu = cell->hdu; + } + + hdu->masks = psArrayAlloc(hdu->images->n); + hdu->weights = psArrayAlloc(hdu->images->n); + psArray *readouts = cell->readouts; + for (int readNum = 0; readNum < hdu->images->n; readNum++) { + psImage *image = hdu->images->data[readNum]; + psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); + psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); + pmReadout *readout = readouts->data[readNum]; + for (int j = 0; j < image->numRows; j++) { + for (int i = 0; i < image->numCols; i++) { + mask->data.U8[j][i] = j + i; + weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]); + } + } + hdu->masks->data[readNum] = mask; + hdu->weights->data[readNum] = weight; + readout->mask = psMemIncrRefCounter(mask); + readout->weight = psMemIncrRefCounter(weight); + } + } + } + } + } } #endif @@ -354,477 +338,455 @@ // Load the calibration frames, if required if (doBias) { + psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias + psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias + if (! pmConfigValidateCamera(camera, biasHeader)) { + psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName); + if (! pmFPARead(bias, biasFile, biasHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName); + exit(EXIT_FAILURE); + } + psFree(biasHeader); + psFitsClose(biasFile); + } + + if (doDark) { + psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark + psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark + if (! pmConfigValidateCamera(camera, darkHeader)) { + psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName); + if (! pmFPARead(dark, darkFile, darkHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName); + exit(EXIT_FAILURE); + } + psFree(darkHeader); + psFitsClose(darkFile); + } + + if (doMask) { + psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask + psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask + if (! pmConfigValidateCamera(camera, maskHeader)) { + psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName); + if (! pmFPARead(mask, maskFile, maskHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName); + exit(EXIT_FAILURE); + } + psFree(maskHeader); + psFitsClose(maskFile); + } + + if (doFlat) { + psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat + if (! flatFile) { + psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName); + exit(EXIT_FAILURE); + } + psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat + if (! pmConfigValidateCamera(camera, flatHeader)) { + psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName); + exit(EXIT_FAILURE); + } + psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName); + if (! pmFPARead(flat, flatFile, flatHeader, database)) { + psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName); + exit(EXIT_FAILURE); + } + psFree(flatHeader); + psFitsClose(flatFile); + } + + psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2")); + + psArray *inputChips = input->chips; // Array of chips in input image + psArray *biasChips = bias->chips; // Array of chips in bias image + psArray *darkChips = dark->chips; // Array of chips in dark image + psArray *maskChips = mask->chips; // Array of chips in mask image + psArray *flatChips = flat->chips; // Array of chips in flat image + for (int i = 0; i < inputChips->n; i++) { + pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image + pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image + pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image + pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image + pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image + + if (! inputChip->valid) { + continue; + } + + psArray *inputCells = inputChip->cells; // Array of cells in input image + psArray *biasCells = biasChip->cells; // Array of cells in bias image + psArray *darkCells = darkChip->cells; // Array of cells in dark image + psArray *maskCells = maskChip->cells; // Array of cells in mask image + psArray *flatCells = flatChip->cells; // Array of cells in flat image + + for (int j = 0; j < inputCells->n; j++) { + pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image + pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image + pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag + pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image + pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image + + if (! inputCell->valid) { + continue; + } + + psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image + if (doBias && biasCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" + " used."); + } + if (doDark && darkCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" + " used."); + } + if (doMask && maskCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" + " used."); + } + if (doFlat && flatCell->readouts->n > 1) { + psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be" + " used."); + } + + pmReadout *biasReadout = NULL; // Bias readout + pmReadout *maskReadout = NULL; // Mask readout + pmReadout *flatReadout = NULL; // Flat readout + psImage *biasImage = NULL; // Bias pixels + psImage *darkImage = NULL; // Dark pixels + float darkTime = 1.0; // Dark time for dark image + psImage *maskImage = NULL; // Mask pixels + + if (doBias) { + biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image + biasImage = biasReadout->image; + } + if (doDark) { + pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image + darkImage = darkReadout->image; + darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME"); + if (darkTime <= 0.0) { + psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime); + exit(EXIT_FAILURE); + } + } + if (doMask) { + maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image + maskImage = maskReadout->image; + } + if (doFlat) { + flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image + } + + for (int k = 0; k < inputReadouts->n; k++) { + pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image + psImage *inputImage = inputReadout->image; // The actual image + + // Mask bad pixels + if (doMask) { + float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION"); + + // Need to change this later to grow the mask by the size of the convolution kernel + (void)pmMaskBadPixels(inputReadout, maskImage, + PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation, + PM_MASK_TRAP, 0); + } + + // Non-linearity correction + if (doNonLin) { + psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA"); + if (! dataItem) { + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to " + "find NONLIN.DATA in recipe --- ignored.\n"); + } else { + switch (dataItem->type) { + case PS_DATA_VECTOR: + { + // These are the polynomial coefficients + psVector *coeff = dataItem->data.V; // The coefficient vector + if (coeff->type.type != PS_TYPE_F64) { + psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version + coeff = temp; + } + psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, + PS_POLYNOMIAL_ORD); + psFree(correction->coeff); + correction->coeff = psMemIncrRefCounter(coeff->data.F64); + (void)pmNonLinearityPolynomial(inputReadout, correction); + psFree(coeff); + psFree(correction); + } + break; + case PS_DATA_STRING: + { + // This is a filename + psString name = dataItem->data.V; // Filename + psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0); + if (psLookupTableRead(table) <= 0) { + psErrorStackPrint(stderr, "Unable to read non-linearity correction file " + "%s --- ignored\n", name); + } else { #ifdef PRODUCTION - psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias + (void)pmNonLinearityLookup(inputReadout, table); #else - psFits *biasFile = psFitsAlloc(biasName); -#endif - psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias - if (! pmConfigValidateCamera(camera, biasHeader)) { - psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName); - if (! pmFPARead(bias, biasFile, biasHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName); - exit(EXIT_FAILURE); - } - psFree(biasHeader); + psVector *influx = table->values->data[0]; + psVector *outflux = table->values->data[1]; + (void)pmNonLinearityLookup(inputReadout, table->values->data[0], + table->values->data[1]); +#endif + } + psFree(table); + } + break; + case PS_DATA_METADATA: + { + // This is a menu + psMetadata *options = dataItem->data.V; // Options with concept values as keys + bool mdok = false; // Success of MD lookup + psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE"); + if (! mdok || ! concept) { + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " + "unable to find NONLIN.SOURCE in recipe --- ignored.\n"); + } else { + psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept); + if (! conceptValueItem) { + psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s " + "for non-linearity correction --- ignored.\n", concept); + } else if (conceptValueItem->type != PS_DATA_STRING) { + psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as" + " expected for non-linearity correction --- ignored.\n", + concept); + } else { + psString conceptValue = conceptValueItem->data.V; + // Get the value of the concept + psMetadataItem *optionItem = psMetadataLookup(options, conceptValue); + if (!optionItem) { + psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA" + " --- ignored.\n", conceptValue); + } else if (optionItem->type == PS_DATA_VECTOR) { + // These are the polynomial coefficients + psVector *coeff = optionItem->data.V; // The coefficient vector + if (coeff->type.type != PS_TYPE_F64) { + psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); + coeff = temp; + } + // Polynomial correction + psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, + PS_POLYNOMIAL_ORD); + psFree(correction->coeff); + correction->coeff = psMemIncrRefCounter(coeff->data.F64); + (void)pmNonLinearityPolynomial(inputReadout, correction); + psFree(coeff); + psFree(correction); + } else if (optionItem->type == PS_DATA_STRING) { + // This is a filename + psString tableName = optionItem->data.V; // The filename + psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0); + int numLines = 0; // Number of lines read from table + if ((numLines = psLookupTableRead(table)) <= 0) { + psErrorStackPrint(stderr, "Unable to read non-linearity " + "correction file %s --- ignored\n", + tableName); + } else { #ifdef PRODUCTION - psFitsClose(biasFile); + (void)pmNonLinearityLookup(inputReadout, table); #else - psFree(biasFile); -#endif - } - - if (doDark) { + printf("XXX: Non-linearity correction from lookup table not " + "yet implemented.\n"); +#endif + } + psFree(table); + } else { + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction " + "desired but unable to interpret NONLIN.DATA for %s" + " --- ignored\n", conceptValue); + } + } + } + } + break; + default: + psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " + "NONLIN.DATA is of invalid type --- ignored.\n"); + } + } + } // Non-linearity correction + + // Bias, dark and overscan subtraction are all merged. + + // Changed the definition of pmSubtractBias to do the dark subtraction as well, but + // it's not in there yet. Once it is, we can get rid of "pedestal". + #ifdef PRODUCTION - psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark + psImage *pedestal = NULL; // Pedestal image (bias + scaled dark) + if (doDark) { + // Dark time for input image + float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME"); + if (inputTime <= 0) { + psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n", + inputTime); + exit(EXIT_FAILURE); + } + + pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*", + psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); + } + if (doBias) { + if (pedestal) { + if (biasImage->numRows != darkImage->numRows || + biasImage->numCols != darkImage->numCols) { + psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: " + "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows, + darkImage->numCols, darkImage->numRows); + exit(EXIT_FAILURE); + } + (void)psBinaryOp(pedestal, pedestal, "+", biasImage); + } else { + pedestal = psMemIncrRefCounter(biasImage); + } + } +#endif + if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { + // Need to get the read direction + int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR"); + if (readdir == 1) { +// psmodule-0.8.0 has confused PM_OVERSCAN_ROWS and PM_OVERSCAN_COLUMNS; bug 608. +#ifdef PRODUCTION + overscanMode = PM_OVERSCAN_ROWS; #else - psFits *darkFile = psFitsAlloc(darkName); -#endif - psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark - if (! pmConfigValidateCamera(camera, darkHeader)) { - psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName); - if (! pmFPARead(dark, darkFile, darkHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName); - exit(EXIT_FAILURE); - } - psFree(darkHeader); + overscanMode = PM_OVERSCAN_COLUMNS; +#endif + } else if (readdir == 2) { #ifdef PRODUCTION - psFitsClose(darkFile); + overscanMode = PM_OVERSCAN_COLUMNS; #else - psFree(darkFile); -#endif - } - - if (doMask) { + overscanMode = PM_OVERSCAN_ROWS; +#endif + } else { + psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n", + readdir); + overscanMode = PM_OVERSCAN_ROWS; + } + } + + if (doBias || doOverscan || doDark) { + psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions #ifdef PRODUCTION - psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask + (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, + overscanStats, overscanBins, overscanFitType, pedestal); #else - psFits *maskFile = psFitsAlloc(maskName); -#endif - psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask - if (! pmConfigValidateCamera(camera, maskHeader)) { - psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName); - if (! pmFPARead(mask, maskFile, maskHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName); - exit(EXIT_FAILURE); - } - psFree(maskHeader); - -#ifdef PRODUCTION - psFitsClose(maskFile); -#else - psFree(maskFile); -#endif - } - - if (doFlat) { -#ifdef PRODUCTION - psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat -#else - psFits *flatFile = psFitsAlloc(flatName); -#endif - if (! flatFile) { - psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName); - exit(EXIT_FAILURE); - } - psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat - if (! pmConfigValidateCamera(camera, flatHeader)) { - psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName); - exit(EXIT_FAILURE); - } - psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName); - if (! pmFPARead(flat, flatFile, flatHeader, database)) { - psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName); - exit(EXIT_FAILURE); - } - psFree(flatHeader); -#ifdef PRODUCTION - psFitsClose(flatFile); -#else - psFree(flatFile); -#endif - } - - psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2")); - - psArray *inputChips = input->chips; // Array of chips in input image - psArray *biasChips = bias->chips; // Array of chips in bias image - psArray *darkChips = dark->chips; // Array of chips in dark image - psArray *maskChips = mask->chips; // Array of chips in mask image - psArray *flatChips = flat->chips; // Array of chips in flat image - for (int i = 0; i < inputChips->n; i++) { - pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image - pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image - pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image - pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image - pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image - - if (! inputChip->valid) { - continue; - } - - psArray *inputCells = inputChip->cells; // Array of cells in input image - psArray *biasCells = biasChip->cells; // Array of cells in bias image - psArray *darkCells = darkChip->cells; // Array of cells in dark image - psArray *maskCells = maskChip->cells; // Array of cells in mask image - psArray *flatCells = flatChip->cells; // Array of cells in flat image - - for (int j = 0; j < inputCells->n; j++) { - pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image - pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image - pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag - pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image - pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image - - if (! inputCell->valid) { - continue; - } - - psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image - if (doBias && biasCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" - " used."); - } - if (doDark && darkCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" - " used."); - } - if (doMask && maskCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" - " used."); - } - if (doFlat && flatCell->readouts->n > 1) { - psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be" - " used."); - } - - pmReadout *biasReadout = NULL; // Bias readout - pmReadout *maskReadout = NULL; // Mask readout - pmReadout *flatReadout = NULL; // Flat readout - psImage *biasImage = NULL; // Bias pixels - psImage *darkImage = NULL; // Dark pixels - float darkTime = 1.0; // Dark time for dark image - psImage *maskImage = NULL; // Mask pixels - - if (doBias) { - biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image - biasImage = biasReadout->image; - } - if (doDark) { - pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image - darkImage = darkReadout->image; - darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME"); - if (darkTime <= 0.0) { - psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime); - exit(EXIT_FAILURE); - } - } - if (doMask) { - maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image - maskImage = maskReadout->image; - } - if (doFlat) { - flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image - } - - for (int k = 0; k < inputReadouts->n; k++) { - pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image - psImage *inputImage = inputReadout->image; // The actual image - - // Mask bad pixels - if (doMask) { - float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION"); - - // Need to change this later to grow the mask by the size of the convolution kernel - (void)pmMaskBadPixels(inputReadout, maskImage, - PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation, - PM_MASK_TRAP, 0); - } - - // Non-linearity correction - if (doNonLin) { - psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA"); - if (! dataItem) { - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to " - "find NONLIN.DATA in recipe --- ignored.\n"); - } else { - switch (dataItem->type) { - case PS_DATA_VECTOR: - { - // These are the polynomial coefficients - psVector *coeff = dataItem->data.V; // The coefficient vector - if (coeff->type.type != PS_TYPE_F64) { - psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version - coeff = temp; - } - psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, - PS_POLYNOMIAL_ORD); - psFree(correction->coeff); - correction->coeff = psMemIncrRefCounter(coeff->data.F64); - (void)pmNonLinearityPolynomial(inputReadout, correction); - psFree(coeff); - psFree(correction); - } - break; - case PS_DATA_STRING: - { - // This is a filename - psString name = dataItem->data.V; // Filename - psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0); - if (psLookupTableRead(table) <= 0) { - psErrorStackPrint(stderr, "Unable to read non-linearity correction file " - "%s --- ignored\n", name); - } else { -#ifdef PRODUCTION - (void)pmNonLinearityLookup(inputReadout, table); -#else - psVector *influx = table->values->data[0]; - psVector *outflux = table->values->data[1]; - (void)pmNonLinearityLookup(inputReadout, table->values->data[0], - table->values->data[1]); -#endif - } - psFree(table); - } - break; - case PS_DATA_METADATA: - { - // This is a menu - psMetadata *options = dataItem->data.V; // Options with concept values as keys - bool mdok = false; // Success of MD lookup - psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE"); - if (! mdok || ! concept) { - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " - "unable to find NONLIN.SOURCE in recipe --- ignored.\n"); - } else { - psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept); - if (! conceptValueItem) { - psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s " - "for non-linearity correction --- ignored.\n", concept); - } else if (conceptValueItem->type != PS_DATA_STRING) { - psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as" - " expected for non-linearity correction --- ignored.\n", - concept); - } else { - psString conceptValue = conceptValueItem->data.V; - // Get the value of the concept - psMetadataItem *optionItem = psMetadataLookup(options, conceptValue); - if (!optionItem) { - psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA" - " --- ignored.\n", conceptValue); - } else if (optionItem->type == PS_DATA_VECTOR) { - // These are the polynomial coefficients - psVector *coeff = optionItem->data.V; // The coefficient vector - if (coeff->type.type != PS_TYPE_F64) { - psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); - coeff = temp; - } - // Polynomial correction - psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1, - PS_POLYNOMIAL_ORD); - psFree(correction->coeff); - correction->coeff = psMemIncrRefCounter(coeff->data.F64); - (void)pmNonLinearityPolynomial(inputReadout, correction); - psFree(coeff); - psFree(correction); - } else if (optionItem->type == PS_DATA_STRING) { - // This is a filename - psString tableName = optionItem->data.V; // The filename - psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0); - int numLines = 0; // Number of lines read from table - if ((numLines = psLookupTableRead(table)) <= 0) { - psErrorStackPrint(stderr, "Unable to read non-linearity " - "correction file %s --- ignored\n", - tableName); - } else { -#ifdef PRODUCTION - (void)pmNonLinearityLookup(inputReadout, table); -#else - printf("XXX: Non-linearity correction from lookup table not " - "yet implemented.\n"); -#endif - } - psFree(table); - } else { - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction " - "desired but unable to interpret NONLIN.DATA for %s" - " --- ignored\n", conceptValue); - } - } - } - } - break; - default: - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but " - "NONLIN.DATA is of invalid type --- ignored.\n"); - } - } - } // Non-linearity correction - - // Bias, dark and overscan subtraction are all merged. - - // Changed the definition of pmSubtractBias to do the dark subtraction as well, but - // it's not in there yet. Once it is, we can get rid of "pedestal". - -#ifdef PRODUCTION - psImage *pedestal = NULL; // Pedestal image (bias + scaled dark) - if (doDark) { - // Dark time for input image - float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME"); - if (inputTime <= 0) { - psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n", - inputTime); - exit(EXIT_FAILURE); - } - - pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*", - psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); - } - if (doBias) { - if (pedestal) { - if (biasImage->numRows != darkImage->numRows || - biasImage->numCols != darkImage->numCols) { - psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: " - "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows, - darkImage->numCols, darkImage->numRows); - exit(EXIT_FAILURE); - } - (void)psBinaryOp(pedestal, pedestal, "+", biasImage); - } else { - pedestal = psMemIncrRefCounter(biasImage); - } - } -#endif - if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { - // Need to get the read direction - int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR"); - if (readdir == 1) { - overscanMode = PM_OVERSCAN_ROWS; - } else if (readdir == 2) { - overscanMode = PM_OVERSCAN_COLUMNS; - } else { - psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n", - readdir); - overscanMode = PM_OVERSCAN_ROWS; - } - } - - if (doBias || doOverscan || doDark) { - psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions -#ifdef PRODUCTION - (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, - overscanStats, overscanBins, overscanFitType, pedestal); -#else - (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, - overscanStats, overscanBins, overscanFitType, biasReadout); -#endif - psFree(inputOverscans); - - // Output overscan fit results, if required - if (doOverscan) { - if (! overscanFit) { - psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n"); - } else { - switch (overscanFitType) { - case PM_FIT_POLYNOMIAL: - { - psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial - psString coeffs = NULL; // String containing the coefficients - for (int i = 0; i < poly->nX; i++) { - psStringAppend(&coeffs, "%e ", poly->coeff[i]); - } - psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n", - coeffs); - psFree(coeffs); - } - break; - case PM_FIT_SPLINE: - { - psSpline1D *spline = (psSpline1D*)overscanFit; // The spline - psString coeffs = NULL; // String containing the coefficients - for (int i = 0; i < spline->n; i++) { - psPolynomial1D *poly = spline->spline[i]; // i-th polynomial - psStringAppend(&coeffs, "%d: ", i); - for (int j = 0; j < poly->nX; j++) { - psStringAppend(&coeffs, "%e ", poly->coeff[i]); - } - psStringAppend(&coeffs, "\n"); - } - psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n", - coeffs); - psFree(coeffs); - } - break; - case PM_FIT_NONE: - break; - default: - psAbort(__func__, "Should never get here!!!\n"); - } - } - } - } - - // XXX: Temporary: until the other functions are altered to do this themselves. - // Trim, so that flat, fringe etc computations are faster. + (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode, + overscanStats, overscanBins, overscanFitType, biasReadout); +#endif + psFree(inputOverscans); + + // Output overscan fit results, if required + if (doOverscan) { + if (! overscanFit) { + psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n"); + } else { + switch (overscanFitType) { + case PM_FIT_POLYNOMIAL: + { + psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial + psString coeffs = NULL; // String containing the coefficients + for (int i = 0; i < poly->nX; i++) { + psStringAppend(&coeffs, "%e ", poly->coeff[i]); + } + psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n", + coeffs); + psFree(coeffs); + } + break; + case PM_FIT_SPLINE: + { + psSpline1D *spline = (psSpline1D*)overscanFit; // The spline + psString coeffs = NULL; // String containing the coefficients + for (int i = 0; i < spline->n; i++) { + psPolynomial1D *poly = spline->spline[i]; // i-th polynomial + psStringAppend(&coeffs, "%d: ", i); + for (int j = 0; j < poly->nX; j++) { + psStringAppend(&coeffs, "%e ", poly->coeff[i]); + } + psStringAppend(&coeffs, "\n"); + } + psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n", + coeffs); + psFree(coeffs); + } + break; + case PM_FIT_NONE: + break; + default: + psAbort(__func__, "Should never get here!!!\n"); + } + } + } + } + + // XXX: Temporary: until the other functions are altered to do this themselves. + // Trim, so that flat, fringe etc computations are faster. #if 0 #ifndef PRODUCTION - psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC"); - psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion); - psFree(inputReadout->image); - inputReadout->image = trimmed; -#endif -#endif - // Flat-field correction - if (doFlat) { - psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n"); + psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC"); + psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion); + psFree(inputReadout->image); + inputReadout->image = trimmed; +#endif +#endif + // Flat-field correction + if (doFlat) { + psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n"); #ifndef PRODUCTION - psImage *dummyImage = psImageAlloc(inputReadout->image->numCols, - inputReadout->image->numRows, - PS_TYPE_U8); - pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1); - (void)pmFlatField(inputReadout, dummyMask, flatReadout); - psFree(dummyMask); - psFree(dummyImage); -#endif - } - - // Fringe subtraction - if (doFringe) { + psImage *dummyImage = psImageAlloc(inputReadout->image->numCols, + inputReadout->image->numRows, + PS_TYPE_U8); + pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1); + (void)pmFlatField(inputReadout, dummyMask, flatReadout); + psFree(dummyMask); + psFree(dummyImage); +#endif + } + + // Fringe subtraction + if (doFringe) { #if 0 - // Placeholder - pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, - psU8 maskVal, int binFactor, psStats *stats, float clipSD); -#endif - } - - // Source identification and photometry - if (doSource) { - } - - // Astrometry - if (doAstrom) { - } - - - } // Iterating over readouts - } // Iterating over cells - - // Testing pmChipMosaic - { - psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip - psFits *mosaicFile = psFitsAlloc("mosaic.fits"); // FITS file to which to write - psFitsWriteImage(mosaicFile, NULL, mosaic, 0); - psFree(mosaicFile); - psFree(mosaic); - } - + // Placeholder + pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, + psU8 maskVal, int binFactor, psStats *stats, float clipSD); +#endif + } + + // Source identification and photometry + if (doSource) { + } + + // Astrometry + if (doAstrom) { + } + + + } // Iterating over readouts + } // Iterating over cells + +#if 0 + { + // Testing pmChipMosaic + psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip + psFits *mosaicFile = psFitsOpen("mosaic.fits", "w"); // FITS file to which to write + psFitsWriteImage(mosaicFile, NULL, mosaic, 0); + psFitsClose(mosaicFile); + psFree(mosaic); + } +#endif + } // Iterating over chips @@ -832,16 +794,20 @@ +#if 0 + { + // Testing pmFPAMorph + int nBad = 0; + psMetadata *newCamera = psMetadataConfigParse(NULL, &nBad, "megacam_splice.config", true); + pmFPA *newFPA = pmFPAConstruct(newCamera); + pmFPAMorph(newFPA, input, true, 0, 0); + //pmFPAPrint(newFPA); + psFits *newFile = psFitsOpen("morph.fits", "w"); + pmFPAWrite(newFile, newFPA, database); + psFitsClose(newFile); + psFree(newFPA); + } +#endif + #if 1 - // Testing pmFPAMorph - psMetadata *newCamera = psMetadataConfigParse(NULL, NULL, "megacam_splice.config", true); - pmFPA *newFPA = pmFPAConstruct(newCamera); - pmFPAMorph(newFPA, input, true, 0, 0); - psFits *newFile = psFitsAlloc("morph.fits"); - pmFPAWrite(newFile, newFPA, database); - psFree(newFile); - psFree(newFPA); -#endif - -#if 0 // Write the output pmFPAWrite(outputFile, input, database); @@ -852,9 +818,5 @@ psLogMsg("phase2", PS_LOG_INFO, "Output completed after %f sec.\n", psTimerMark("phase2")); -#ifdef PRODUCTION psFitsClose(outputFile); -#else - psFree(outputFile); -#endif psFree(arguments); @@ -874,5 +836,5 @@ psTimerStop(); -#if 1 +#if 0 psMemCheckCorruption(true); psMemBlock **leaks = NULL; // List of leaks @@ -881,10 +843,10 @@ #if 1 for (int i = 0; i < nLeaks; i++) { - psMemBlock *mb = leaks[i]; - printf("Memory leak %ld (%ld):\n" - "\tFile %s, line %d, size %d\n" - "\tPosts: %x %x %x\n", - mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, - *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); + psMemBlock *mb = leaks[i]; + printf("Memory leak %ld (%ld):\n" + "\tFile %s, line %d, size %d\n" + "\tPosts: %x %x %x\n", + mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock, + *(void**)((int8_t *)(mb + 1) + mb->userMemorySize)); } #endif