Index: /trunk/archive/modules/src/phase4.c =================================================================== --- /trunk/archive/modules/src/phase4.c (revision 2185) +++ /trunk/archive/modules/src/phase4.c (revision 2185) @@ -0,0 +1,212 @@ +#include +#include + +// Configuration for Phase 4 +typedef struct { + // Input parameters + int nInputs; // Number of input files + const char *camera; // Name of the camera configuration file + const psList *inNames; // Input file names + unsigned long skyId; // Sky tangent plane identification number + // Output parameters + const char *combined; // Combined image filename + const char *subtracted; // Subtracted image filename + const char *variables; // Output variable candidates filename + const char *recipeFile; // Recipe filename +} phase4Config; + + +// A tangent plane (for combined image and static sky) +// Assume we're dealing with a rectangle only for now. +typedef struct { + long int id; // Identification number --- key in MDDB + double ra0, dec0; // Tangent point --- from MDDB + double sizeRA, sizeDec; // Size in RA and Dec --- from MDDB + double scale; // Pixel scale --- from MDDB + char *name; // File name --- from MDDB or image server + psImage *pixels; // The pixels that comprise the tangent plane + psImage *mask; // Mask image + psSomeCollectionOfObjects *objects; // Objects in the image +} psTangentPlane; + +void main (int argc, char **argv) +{ + // Bury the command line parsing into some function that just works + phase4Config *config = parseCommandLine(argc, argv); // Parse the command line + + // Read the recipe metadata + psMetadata *recipe = psMetadataParseConfig(NULL, config->recipeFile, false); + + // Don't particularly care about the camera and site details --- Phase 4 doesn't care where + // the data comes from, but it does need the following for each image: + // * Gain and read noise for noise model --- from the camera configuration + // * Background value for noise model (assume constant over the region of interest) --- from image MD + // * Photometric calibration for relative and absolute scalings --- from image MD + // * Astrometric calibration to map to the sky --- from file specified in image MD + // We need to be told which camera; assume all data is from the same camera (easy to generalise + // if this is not true). + + // Read the camera configuration + psMetadata *camera = psMetadataParseConfig(NULL, config->camera, false); + // What are the keywords for gain, readnoise, photometric calibration, background? + // This basically treats the "camera" configuration as a translation table. + const char *gainHeader = psMetadataLookup(camera, "GAIN"); + const char *readnoiseHeader = psMetadataLookup(camera, "READNOISE"); + const char *photcalHeader = psMetadataLookup(camera, "PHOTCAL"); + const char *backgroundHeader = psMetadataLookup(camera, "BACKGRND"); + const char *astromHeader = psMetadataLookup(camera, "ASTROM"); + + // Initialise the various arrays/vectors we need + psArray *headers = psArrayAlloc(config->nInputs); // The image headers + psArray *inputs = psArrayAlloc(config->nInputs); // The images + psArray *astrometry = psArrayAlloc(config->nInputs); // The astrometry data for each image + psVector *gain = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Gains + psVector *readnoise = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Read noises + psVector *photcal = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Photometric calibrations + psVector *background = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Background flux values + + // This function queries the MDDB to get sizeRA,sizeDec,ra0,dec0,scale,name for the tangent plane of + // interest + psTangentPlane *staticSky = psTangentPlaneFromDB(config->sky); + + // Read in each of the inputs + for (psListSetIterator(config->inNames, PS_LIST_HEAD), // Set the iterator and get the first one + char *inName = psListGetNext(input), + int i = 0; + inName != NULL; // Keep going until nothing left in the list + inName = psListGetNext(input), // Get the next one on the list + i++) { + + // Read the header + psMetadata *header = psFITSReadHeader(inName); + // Put it in the array of headers for future reference + headers->data[i] = header; + + // Check the validity of the given header based on expectactions for this camera + status = pmCameraValidateHeaders(headers, camera); + // Need to act if status is bad --- generate a warning and ignore this image + + // Get the gain, readnoise, photometric calibration, background + // I think this is not strictly correct --- doesn't psMetadataLookup return a void* ? + gain->data.F32[i] = (psF32) psMetadataLookup(header, gainHeader); + readnoise->data.F32[i] = (psF32) psMetadataLookup(header, readnoiseHeader); + photcal->data.F32[i] = (psF32) psMetadataLookup(header, photcalHeader); + background->data.F32[i] = (psF32) psMetadataLookup(header, backgroundHeader); + + // Get the name of the astrometry file (XML) and generate the FPA structure from it (no pixels) + const char *astromFile = psMetadataLookup(header, astromHeader); + psFPA *fpa = psFPAFromAstrometry(astromFile); + astrometry->data[i] = fpa; + + // Iterate over the FPA structure + for (int chipnum = 0, int nChips = fpa->chips->n; chipnum < nChips; chipnum++) { + psChip *chip = fpa->chips->data[chipnum]; // The chip of interest + for (int cellnum = 0, int nCells = chip->cells->n; cellnum < nCells; cellnum++) { + psCell *cell = chip->cells->data[cellnum]; // The cell of interest + + // QUESTION: Do we want to bother about the cells that were used for fast guiding? + // Assume that we care --- may as well. + + // QUESTION: Are multiple readouts stacked in phase 4 or earlier? + // Assume that they should be stacked here. + + // QUESTION: How do we know the names of the chip, cell, readout in the FITS file? + // Assume that each of psFPA, psChip, psCell, psReadout contains a "const char *source" + // which contains some sort of URI, e.g.: + // FPA: fpa12345_12.fits --- This particular FPA contains only a single OTA + // Chip: fpa12345_12.fits --- The OTA is entirely contained in the FITS file + // Cell: fpa12345_12.fits|cell34 --- The cell has extension name "cell34" + // Readout: fpa12345_12.fits|cell34|56 --- The readout is the 56th slice in the 3rd axis + // Assume further that, given the "source", there exists a module that loads the pixels. + + + // If there are multiple readouts, then stack them + if (cell->readouts->n > 1) { + psList *readouts = psListAlloc(NULL); // List of readouts for stacking + + // Read the images in + for (int readoutNum = 0; readoutNum < cell->readouts->n; readoutNum++) { + (void)psReadoutRead(cell->readouts->data[readoutNum]); // Read the readout into memory + psListAdd(readouts, PS_LIST_TAIL, cell->readouts->data[readoutNum]); // Put on list + } + + // QUESTION: Do we want to bother about combining these "properly", or just sum them? + // Assume that we just sum them, which saves measuring the zeros and scales. CR rejection + // can be done when we combine multiple images (consecutive or from multiple telescopes). + + // Perhaps we want a pmReadoutStack module, which would convert an array of readouts from + // fast guiding to a single pseudo-readout. + + // Just do a simple mean + psStats *combineStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics for combination + pmCombineParams *combine = pmCombineParamsAlloc(combineStats, 0, 0.0, 0.0, 0); // How to + // combine + psImage *stack = pmReadoutCombine(NULL, readouts, combine, 0.0, 1.0, false, 1.0, 0.0); + // And now multiply by the number to get the sum + (void)psImageBinaryOp(stack, stack, "*", psScalarAlloc(cell->readouts->n, PS_TYPE_U16)); + + // Probably also have to multiply by some exposure time ratio to scale with the rest of + // the focal plane --- neglected here. + + // Clean up the temporary steps + psFree(readouts); + psFree(combine); + psFree(combineStats); + + // Stuff the stack back into the cell --- probably should play with the metadata, but I'm + // going to ignore that for now. + + // Clean out the old structure + for (int readoutNum = 0; readoutNum < cell->readouts->n; readoutNum++) { + psFree(cell->readouts->data[readoutNum]); + } + // Resize and stuff the stack in. + psArrayRealloc(cell->readouts, 1); + cell->readouts->data[0] = stack; + } else { + // Otherwise, just read the cell in + cell->readouts->data[0] = psReadoutRead(cell->readouts->data[0]); + } // Reading in readouts + } // Cells + } // Chips + } // FPAs + + // Now all the data is in memory. Combine the multiple images + psTangentPlane *combined = pmFPACombine(NULL, staticSky, inputs, gain, readnoise, photcal, background); + + // Write the combined image out + psTangentPlaneWriteFITS(combined, config->combined); + + // Run the object detection/measurement --- this is not well thought out yet + (void)psDetectAndMeasureObjects(combined, otherParameters); + + // Load the static sky into the structure already defined + (void)psTangentPlaneRead(staticSky); + + // Get the seeing ratio --- involves matching stars and measuring the mean ratio of FWHM + float seeingRatio = pmSeeingRatio(combined, staticSky, otherParameters); + + // Do the subtraction --- involves measuring the convolution kernel, applying the kernel and subtracting + psTangentPlane *subtracted = NULL; + if (seeingRatio < 1.0) { + subtracted = pmImageSubtraction(NULL, combined, staticSky, otherParameters); + } else { + subtracted = pmImageSubtraction(NULL, staticSky, combined, otherParameters); + } + + // Write out the subtracted image + psTangentPlaneWriteFITS(subtraction, config->subtracted); + + // Find and measure objects on the subtracted image + (void)psDetectAndMeasureObjects(subtracted, otherParameters); + + // Filter the list of candidate variable sources + psList *variables = pmFilterVariables(subtracted, combined, staticSky, otherParameters); + + // Output the list of variable sources + (void)pmOutputVariables(variables, config->variables); + + // All done. +} + +