Index: /trunk/archive/scripts/src/papPhase2.c =================================================================== --- /trunk/archive/scripts/src/papPhase2.c (revision 4755) +++ /trunk/archive/scripts/src/papPhase2.c (revision 4755) @@ -0,0 +1,599 @@ +#include +#include "pslib.h" +#include "papmodule.h" + +// Phase 2 needs to: +// * Read configurations +// * Read in the images +// * Flag bad pixels +// * Non-linearity correction +// * Bias/dark subtraction (MHPCC) +// * Flat-field (MHPCC) +// * Fringe subtraction (Not for now) +// * Source identification and photometry (IPP prototype) +// * Astrometry (IPP prototype soon) +// * Write calibrated image +// * Write source catalogue, astrometry, etc. + +// Currently neglecting different input types (F32, S32, U16, etc) + +// Would like to fold in Nebulous at some stage as well. + +// Currently neglect to convolve the reference frames with the orthogonal transfer kernel + +// 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". + +void help(psMetadata *arguments, // Command-line arguments + const char *programName // Name of the program = argv[0] + ) +{ + printf("Pan-STARRS Phase 2 processing\n\n"); + printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", programName); + psArgumentHelp(arguments); +} + + + +int main(int argc, char *argv[]) +{ + // Parse the configurations + psMetadata *site = NULL; // Site configuration + psMetadata *camera = NULL; // Camera configuration + psMetadata *recipe = NULL; // Recipe configuration + if (! pmConfigRead(&site, &camera, &recipe, &argc, argv, "moduleName")) { + psErrorStackPrint("Can't find site configuration!\n"); + exit(EXIT_FAILURE); + } + + // Parse other command-line arguments + psMetadata *arguments = psMetadataAlloc(); // The arguments, with default values + psMetadataAddStr(arguments, PS_LIST_TAIL, "bias", "Name of the bias image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "dark", "Name of the dark image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "flat", "Name of the flat-field image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "mask", "Name of the mask image", NULL); + psMetadataAddS32(arguments, PS_LIST_TAIL, "chip", "Chip number to process (if positive)", -1); + if (! psArgumentParse(arguments, &argc, argv) || argc != 2) { + help(arguments, argv[0]); + 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 + const char *flatName = psMetadataLookupString(NULL, arguments, "flat"); // Name of flat-field image + const char *maskName = psMetadataLookupString(NULL, arguments, "mask"); // Name of mask image + const int chipNum = psMetadataLookupS32(NULL, arguments, "chip"); // Chip number to work on + + // Open the input + psFits *inputFile = psFitsOpen(inputName, "r"); // File handle for FITS file + if (! inputFile) { + psErrorStackPrint("Can't open input image: %s\n", inputName); + exit(EXIT_FAILURE); + } + psMetadata *header = psFitsReadHeader(NULL, inputFile); // FITS header + psDB *database = pmConfigDB(site); // Database handle + + // 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 + psFits *outputFile = psFitsOpen(outputName, "w"); + if (! outputFile) { + psErrorStackPrint("Can't open output image: %s\n", outputName); + exit(EXIT_FAILURE); + } + psString outputMaskName = psStringCopy(outputName); + (void)psStringAppend(outputMaskName, ".mask"); + psFits *outputMaskFile = psFits(outputMaskName, "w"); + if (! outputMaskFile) { + psErrorStackPrint("Can't open output mask image: %s\n", outputMaskName); + exit(EXIT_FAILURE); + } + psFree(outputMaskName); + + // Get camera configuration from header if not already defined + if (! camera) { + camera = pmConfigCameraFromHeader(site, header); + if (! camera) { + psErrorStackPrint("Can't find camera configuration!\n"); + exit(EXIT_FAILURE); + } + } else if (! pmConfigValidateCamera(camera, inputHeader)) { + psError("phase2", true, "%s does not seem to be from the camera.\n", inputName); + exit(EXIT_FAILURE); + } + if (! recipe && !(recipe = pmConfigRecipeFromCamera(camera, "moduleName"))) { + psErrorStackPrint("Can't find recipe configuration!\n"); + exit(EXIT_FAILURE); + } + + // Construct camera in preparation for reading + pmFPA *input = pmFPAConstruct(camera, db); + pmFPA *mask = pmFPAConstruct(camera, db); + pmFPA *bias = pmFPAConstruct(camera, db); + pmFPA *dark = pmFPAConstruct(camera, db); + pmFPA *flat = pmFPAConstruct(camera, db); + + // 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 + pmOverscanAxis overscanMode = PM_OVERSCAN_NONE; // Axis for overscan + pmFit overscanFit = PM_FIT_NONE; // Fit type for overscan + int overscanBins = 1; // Number of pixels per bin for overscan + psStats *overscanStats = NULL; // Statistics for overscan + + if (psMetadataLookupBool(NULL, recipe, "MASK")) { + if (maskName) { + 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; + } + if (psMetadataLookupBool(NULL, recipe, "BIAS")) { + if (biasName) { + 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 (darkName) { + 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); + } + psFree(mode); + psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT"); + if (strcasecmp(fit, "POLYNOMIAL") == 0) { + overscanFit = PM_FIT_POLYNOMIAL; + } else if (strcasecmp(fit, "SPLINE") == 0) { + overscanFit = PM_FIT_SPLINE; + } 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); + } + psFree(fit); + overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN"); + if (overscanBins <= 0) { + psErrorStackPrint("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("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 (flatName) { + 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("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); + } + + // Read in the input pixels + if (! pmFPARead(input, inputFile)) { + psErrorStackPrint("Unable to populate camera from FITS file: %s\n", inputName); + exit(EXIT_FAILURE); + } + psFitsClose(inputFile); + + // Load the calibration frames, if required + if (biasName) { + psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias + 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); + exit(EXIT_FAILURE); + } + psFree(biasHeader); + if (! pmFPARead(bias, biasFile)) { + psErrorStackPrint("Unable to populate bias camera from fits FITS: %s\n", biasName); + exit(EXIT_FAILURE); + } + psFitsClose(biasFile); + } + + if (darkName) { + psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark + 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); + exit(EXIT_FAILURE); + } + psFree(darkHeader); + if (! pmFPARead(dark, darkFile)) { + psErrorStackPrint("Unable to populate dark camera from fits FITS: %s\n", darkName); + exit(EXIT_FAILURE); + } + psFitsClose(darkFile); + } + + if (maskName) { + psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask + 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); + exit(EXIT_FAILURE); + } + psFree(maskHeader); + if (! pmFPARead(mask, maskFile)) { + psErrorStackPrint("Unable to populate mask camera from fits FITS: %s\n", maskName); + exit(EXIT_FAILURE); + } + psFitsClose(maskFile); + } + + if (flatName) { + psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat + 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); + exit(EXIT_FAILURE); + } + psFree(flatHeader); + if (! pmFPARead(flat, flatFile)) { + psErrorStackPrint("Unable to populate flat camera from fits FITS: %s\n", flatName); + exit(EXIT_FAILURE); + } + psFitsClose(flatFile); + } + + psArray *inputChips = inputs->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[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 + + if (! inputCell->valid) { + continue; + } + + psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image + if (biasCell->readouts->n != 1) { + psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be" + " used."); + } + if (darkCell->readouts->n != 1) { + psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be" + " used."); + } + if (maskCell->readouts->n != 1) { + psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be" + " used."); + } + if (flatCell->readouts->n != 1) { + psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be" + " used."); + } + + 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 + biasImage = biasReadout->image; + } + if (darkName) { + pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image + darkImage = darkReadout->image; + darkTime = pmReadoutGetDarkTime(darkReadout); + if (darkTime <= 0.0) { + psErrorStackPrint("DARKTIME for dark image (%f) is non-positive.\n", darkTime); + exit(EXIT_FAILURE); + } + } + 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 + } + + 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 + psList *inputBias = inputReadout->bias; // List of overscan bias regions + + // Mask bad pixels + if (doMask) { + // 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); + } + + // 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_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 + 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(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("Unable to read non-linearity correction file %s --- " + "ignored\n", tableName); + } else { + (void)pmNonLinearityLookup(inputReadout, 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); + } 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); + 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(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("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); + } + } + } + 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". + + psImage *pedestal = NULL; // Pedestal image (bias + scaled dark) + if (doDark) { + float inputTime = pmReadoutGetDarkTime(inputReadout); // Dark time for input image + if (inputTime <= 0) { + psErrorStackPrint("DARKTIME for input image (%f) is non-positive.\n", inputTime); + exit(EXIT_FAILURE); + } + + pedestal = psBinaryOp(NULL, darkImage, "*", + psScalarAlloc(inputTime * darkTime, PS_TYPE_F32)); + } + if (doBias) { + if (pedestal) { + if (biasImage->numRows != darkImage->numRows || + biasImage->numCols != darkImage->numCols) { + psError("phase2", 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); + } + } + + if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) { + // Need to get the read direction + int readdir = pmCellGetReaddir(inputCell); // Read direction + if (readdir == 1) { + overscanMode = PM_OVERSCAN_ROWS; + } else if (readdir == 2) { + overscanMode = PM_OVERSCAN_COLUMNS; + } else { + psErrorStackPrint("CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n", readdir); + overscanMode = PM_OVERSCAN_ROWS; + } + } + + 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]); + } + 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"); + } + } + + // Flat-field correction + if (doFlat) { + (void)pmFlatField(inputReadout, maskReadout, flatReadout); + } + + // 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 + } // Iterating over chips + + // Write the output + pmFPAWrite(outputFile, input); + pmFPAWriteMask(outputMaskFile, input); + psFitsClose(outputFile); + psFitsClose(outputMaskFile); + + psFree(input); + psFree(mask); + psFree(bias); + psFree(dark); + psFree(flat); + +}