Index: trunk/archive/scripts/src/phase2/papPhase2.c =================================================================== --- trunk/archive/scripts/src/phase2/papPhase2.c (revision 4820) +++ trunk/archive/scripts/src/phase2/papPhase2.c (revision 5104) @@ -4,10 +4,16 @@ #include "pslib.h" #include "psAdditionals.h" - -#include "pmSubtractBias.h" #include "pmFPA.h" #include "pmConfig.h" #include "pmFPAConstruct.h" +#include "pmFPARead.h" +#include "pmFPAConceptsGet.h" +#include "pmFPAWrite.h" + +#include "pmFlatField.h" +#include "pmMaskBadPixels.h" +#include "pmNonLinear.h" +#include "pmSubtractBias.h" // Phase 2 needs to: @@ -37,5 +43,5 @@ #define RECIPE "PHASE2" // Name of the recipe to use -static psMemoryId memPrintAlloc(const psMemBlock *mb) +static psMemId memPrintAlloc(const psMemBlock *mb) { printf("Allocated memory block %lld (%lld):\n" @@ -46,5 +52,5 @@ return 0; } -static psMemoryId memPrintFree(const psMemBlock *mb) +static psMemId memPrintFree(const psMemBlock *mb) { printf("Freed memory block %lld (%lld):\n" @@ -69,11 +75,11 @@ #if 0 // Hunting memory leaks - psMemAllocateCallbackSetID(3539); - psMemFreeCallbackSetID(3539); + psMemAllocateCallbackSetID(23289); + psMemFreeCallbackSetID(23289); psMemAllocateCallbackSet(memPrintAlloc); psMemFreeCallbackSet(memPrintFree); #endif -// psTraceSetLevel(".", 10); + psTraceSetLevel("pmConfigCameraFromHeader", 10); // Parse the configurations @@ -88,9 +94,9 @@ // Parse other command-line arguments psMetadata *arguments = psMetadataAlloc(); // The arguments, with default values - psMetadataAddStr(arguments, PS_LIST_TAIL, "-bias", "Name of the bias image", ""); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-dark", "Name of the dark image", ""); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-flat", "Name of the flat-field image", ""); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", "Name of the mask image", ""); - psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", "Chip number to process (if positive)", -1); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-bias", 0, "Name of the bias image", ""); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-dark", 0, "Name of the dark image", ""); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-flat", 0, "Name of the flat-field image", ""); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", 0, "Name of the mask image", ""); + 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"); @@ -113,5 +119,4 @@ printf("Mask: %s\n", maskName); printf("Chip: %d\n", chipNum); - psFree(arguments); // Open the input @@ -132,5 +137,4 @@ #endif -#if 0 // 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 @@ -144,4 +148,6 @@ exit(EXIT_FAILURE); } + +#if 0 psString outputMaskName = psStringCopy(outputName); (void)psStringAppend(&outputMaskName, ".mask"); @@ -158,5 +164,4 @@ #endif - // Get camera configuration from header if not already defined if (! camera) { @@ -176,9 +181,9 @@ // Construct camera in preparation for reading - pmFPA *input = pmFPAConstruct(camera, database); - pmFPA *mask = pmFPAConstruct(camera, database); - pmFPA *bias = pmFPAConstruct(camera, database); - pmFPA *dark = pmFPAConstruct(camera, database); - pmFPA *flat = pmFPAConstruct(camera, database); + pmFPA *input = pmFPAConstruct(camera); + pmFPA *mask = pmFPAConstruct(camera); + pmFPA *bias = pmFPAConstruct(camera); + pmFPA *dark = pmFPAConstruct(camera); + pmFPA *flat = pmFPAConstruct(camera); // Set various tasks @@ -198,5 +203,5 @@ if (psMetadataLookupBool(NULL, recipe, "MASK")) { - if (maskName) { + if (strlen(maskName) > 0) { doMask = true; } else { @@ -209,5 +214,5 @@ } if (psMetadataLookupBool(NULL, recipe, "BIAS")) { - if (biasName) { + if (strlen(biasName) > 0) { doBias = true; } else { @@ -217,5 +222,5 @@ } if (psMetadataLookupBool(NULL, recipe, "DARK")) { - if (darkName) { + if (strlen(darkName) > 0) { doDark = true; } else { @@ -260,6 +265,7 @@ } } + if (psMetadataLookupBool(NULL, recipe, "FLAT")) { - if (flatName) { + if (strlen(flatName) > 0) { doFlat = true; } else { @@ -269,5 +275,4 @@ } -#if 0 // Chip selection if (chipNum >= 0) { @@ -282,5 +287,5 @@ // Read in the input pixels - if (! pmFPARead(input, inputFile)) { + if (! pmFPARead(input, inputFile, header, database)) { psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName); exit(EXIT_FAILURE); @@ -292,67 +297,100 @@ #endif + pmFPAPrint(input); // Load the calibration frames, if required - if (biasName) { + if (doBias) { +#ifdef PRODUCTION psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias +#else + psFits *biasFile = psFitsAlloc(biasName); +#endif psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias if (! pmConfigValidateCamera(camera, biasHeader)) { - psError("phase2", true, "Bias (%s) does not seem to be from the same camera.\n", biasName); + psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName); exit(EXIT_FAILURE); } psFree(biasHeader); - if (! pmFPARead(bias, biasFile)) { + if (! pmFPARead(bias, biasFile, NULL, database)) { psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName); exit(EXIT_FAILURE); } +#ifdef PRODUCTION psFitsClose(biasFile); - } - - if (darkName) { +#else + psFree(biasFile); +#endif + } + + if (doDark) { +#ifdef PRODUCTION psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark +#else + psFits *darkFile = psFitsAlloc(darkName); +#endif psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark if (! pmConfigValidateCamera(camera, darkHeader)) { - psError("phase2", true, "Dark (%s) does not seem to be from the same camera.\n", darkName); + psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName); exit(EXIT_FAILURE); } psFree(darkHeader); - if (! pmFPARead(dark, darkFile)) { + if (! pmFPARead(dark, darkFile, NULL, database)) { psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName); exit(EXIT_FAILURE); } +#ifdef PRODUCTION psFitsClose(darkFile); - } - - if (maskName) { +#else + psFree(darkFile); +#endif + } + + if (doMask) { +#ifdef PRODUCTION psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask +#else + psFits *maskFile = psFitsAlloc(maskName); +#endif psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask if (! pmConfigValidateCamera(camera, maskHeader)) { - psError("phase2", true, "Mask (%s) does not seem to be from the same camera.\n", maskName); + psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName); exit(EXIT_FAILURE); } psFree(maskHeader); - if (! pmFPARead(mask, maskFile)) { + if (! pmFPARead(mask, maskFile, NULL, database)) { psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName); exit(EXIT_FAILURE); } +#ifdef PRODUCTION psFitsClose(maskFile); +#else + psFree(maskFile); +#endif } - if (flatName) { + if (doFlat) { +#ifdef PRODUCTION psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat +#else + psFits *flatFile = psFitsAlloc(flatName); +#endif psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat if (! pmConfigValidateCamera(camera, flatHeader)) { - psError("phase2", true, "Flat (%s) does not seem to be from the same camera.\n", flatName); + psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName); exit(EXIT_FAILURE); } psFree(flatHeader); - if (! pmFPARead(flat, flatFile)) { + if (! pmFPARead(flat, flatFile, NULL, database)) { psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName); exit(EXIT_FAILURE); } +#ifdef PRODUCTION psFitsClose(flatFile); +#else + psFree(flatFile); +#endif } - psArray *inputChips = inputs->chips; // Array of chips in input image + 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 @@ -377,9 +415,9 @@ for (int j = 0; j < inputCells->n; j++) { - pmCell *inputCell = inputCells->data[i]; // Cell of interest in input image - pmCell *biasCell = biasCells->data[i]; // Cell of interest in bias image - pmCell *darkCell = darkCells->data[i]; // Cell of interest in dark imag - pmCell *maskCell = maskCells->data[i]; // Cell of interest in mask image - pmCell *flatCell = flatCells->data[i]; // Cell of interest in flat image + 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) { @@ -388,33 +426,37 @@ psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image - if (biasCell->readouts->n != 1) { + if (doBias && biasCell->readouts->n > 1) { psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" " used."); } - if (darkCell->readouts->n != 1) { + if (doDark && darkCell->readouts->n > 1) { psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" " used."); } - if (maskCell->readouts->n != 1) { + if (doMask && maskCell->readouts->n > 1) { psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" " used."); } - if (flatCell->readouts->n != 1) { + 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 - - if (biasName) { - pmReadout *biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image + psImage *maskImage = NULL; // Mask pixels + + if (doBias) { + biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image biasImage = biasReadout->image; } - if (darkName) { + if (doDark) { pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image darkImage = darkReadout->image; - darkTime = pmReadoutGetDarkTime(darkReadout); + darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME"); if (darkTime <= 0.0) { psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime); @@ -422,9 +464,10 @@ } } - if (maskName) { - pmReadout *maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image - } - if (flatName) { - pmReadout *flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image + 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 } @@ -432,12 +475,14 @@ pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image psImage *inputImage = inputReadout->image; // The actual image - psList *inputBias = inputReadout->bias; // List of overscan bias regions + psList *inputOverscans = inputReadout->overscans; // List of overscan bias regions // 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, maskReadout, - PS_MASK_TRAP | PS_MASK_BADCOL | PS_MASK_SAT, saturation, - PS_MASK_TRAP, 0); + (void)pmMaskBadPixels(inputReadout, maskImage, + PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation, + PM_MASK_TRAP, 0); } @@ -450,91 +495,106 @@ } else { switch (dataItem->type) { - case PS_META_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 + case PS_META_VEC: + { + // 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 + psFree(coeff); + coeff = temp; + } + psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, + coeff->n + 1); + for (int i = 0; i < coeff->n; i++) { + correction->coeff[i] = coeff->data.F64[i]; + } + (void)pmNonLinearityPolynomial(inputReadout, correction); psFree(coeff); - coeff = temp; + psFree(correction); } - psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, - coeff->n + 1); - for (int i = 0; i < coeff->n; i++) { - correction->coeff[i] = coeff->data.F64[i]; - } - (void)pmNonLinearityPolynomial(inputReadout, correction); - psFree(coeffCopy); - psFree(coeff); - psFree(correction); break; case PS_META_STR: - // 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", tableName); - } else { - (void)pmNonLinearityLookup(inputReadout, table); + { + // 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 + printf("XXX: Non-linearity with lookup table not yet implemented.\n"); +#endif + } + psFree(table); } - psFree(table); break; case PS_META_META: - // This is a menu - psMetadata *options = dataItem->data.V; // Options with concept values as keys - bool mdok = false; // Success of MD lookup - psString concept = psMetadataLookupStr(&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_META_STR) { - psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as" - " expected for non-linearity correction --- ignored.\n", - concept); + { + // 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 { - 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_META_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); + 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_META_STR) { + 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_META_VEC) { + // 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); + psFree(coeff); + coeff = temp; + } + // Polynomial correction + psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, + coeff->n + 1); + for (int i = 0; i < coeff->n; i++) { + correction->coeff[i] = coeff->data.F64[i]; + } + (void)pmNonLinearityPolynomial(inputReadout, correction); psFree(coeff); - coeff = temp; + psFree(correction); + } else if (optionItem->type == PS_META_STR) { + // 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); } - // Polynomial correction - psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, - coeff->n + 1); - for (int i = 0; i < coeff->n; i++) { - correction->coeff[i] = coeff->data.F64[i]; - } - (void)pmNonLinearityPolynomial(inputReadout, correction); - psFree(coeffCopy); - psFree(coeff); - psFree(correction); - } else if (optionItem->type == PS_META_STR) { - // This is a filename - psString tableName = optionItem->data.V; // The filename - psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0); - if ((int numLines = psLookupTableRead(table)) <= 0) { - psErrorStackPrint(stderr, "Unable to read non-linearity " - "correction file %s --- ignored\n", tableName); - } else { - (void)pmNonLinearityLookup(inputReadout, table); - } - psFree(table); - } else { - psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction " - "desired but unable to interpret NONLIN.DATA for %s" - " --- ignored\n", conceptValue); } } @@ -553,7 +613,9 @@ // 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) { - float inputTime = pmReadoutGetDarkTime(inputReadout); // Dark time for input image + // 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", @@ -562,6 +624,6 @@ } - pedestal = psBinaryOp(NULL, darkImage, "*", - psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); + pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*", + psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); } if (doBias) { @@ -569,5 +631,5 @@ if (biasImage->numRows != darkImage->numRows || biasImage->numCols != darkImage->numCols) { - psError("phase2", true, "Bias and dark images have different dimensions: " + 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); @@ -579,8 +641,8 @@ } } - +#endif if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { // Need to get the read direction - int readdir = pmCellGetReaddir(inputCell); // Read direction + int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR"); if (readdir == 1) { overscanMode = PM_OVERSCAN_ROWS; @@ -594,39 +656,53 @@ } - void *overscanResult = NULL; // Result of overscan fit - (void)pmSubtractBias(inputImage, overscanResult, inputBias, overscanMode, overscanStat, - overscanBin, overscanFit, pedestal); - - // Output overscan fit results, if required - if (doOverscan) { - switch (overscanFit) { - case PM_FIT_POLYNOMIAL: - psPolynomial1D *poly = (psPolynomial1D*)overscanResult; // The polynomial - psString coeffs = NULL; // String containing the coefficients - for (int i = 0; i < poly->n; i++) { - psStringAppend(&coeffs, "%.2f ", poly->coeffs[i]); + printf("mode: %d\n", overscanMode); + + if (doBias || doOverscan || doDark) { + void *overscanResult = NULL; // Result of overscan fit +#ifdef PRODUCTION + (void)pmSubtractBias(inputReadout, overscanResult, inputOverscans, overscanMode, + overscanStats, overscanBins, overscanFit, pedestal); +#else + (void)pmSubtractBias(inputReadout, overscanResult, inputOverscans, overscanMode, + overscanStats, overscanBins, overscanFit, biasReadout); +#endif + + // Output overscan fit results, if required + if (doOverscan) { + switch (overscanFit) { + case PM_FIT_POLYNOMIAL: + { + psPolynomial1D *poly = (psPolynomial1D*)overscanResult; // The polynomial + psString coeffs = NULL; // String containing the coefficients + for (int i = 0; i < poly->n; i++) { + psStringAppend(&coeffs, "%.2f ", 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*)overscanResult; // 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->n; j++) { + psStringAppend(&coeffs, "%.2f ", 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"); } - psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n", - coeffs); - psFree(coeffs); - break; - case PM_FIT_SPLINE: - psSpline1D *spline = (psSpline1D*)overscanResult; // 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->n; j++) { - psStringAppend(&coeffs, "%.2f ", poly->coeffs[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("Should never get here!!!\n"); } } @@ -634,5 +710,5 @@ // Flat-field correction if (doFlat) { - (void)pmFlatField(inputReadout, maskReadout, flatReadout); + (void)pmFlatField(inputReadout, flatReadout); } @@ -659,16 +735,15 @@ // Write the output - pmFPAWrite(outputFile, input); - pmFPAWriteMask(outputMaskFile, input); + pmFPAWrite(outputFile, input, database); +// pmFPAWriteMask(outputMaskFile, input); #ifdef PRODUCTION psFitsClose(outputFile); - psFitsClose(outputMaskFile); +// psFitsClose(outputMaskFile); #else psFree(outputFile); - psFree(outputMaskFile); -#endif - -#endif - +// psFree(outputMaskFile); +#endif + + psFree(arguments); psFree(site); psFree(header); @@ -681,5 +756,4 @@ psFree(flat); psFree(overscanStats); - psFree(inputFile); #if 1