Index: /branches/pap_branch_080207/ppStack/src/ppStackCamera.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackCamera.c (revision 16366) +++ /branches/pap_branch_080207/ppStack/src/ppStackCamera.c (revision 16367) @@ -10,4 +10,61 @@ #include "ppStack.h" + + +// Define an output convolved image file +static pmFPAfile *defineOutputConvolved(const char *name, // FPA file name + pmFPA *fpa, // FPA to bind + const pmConfig *config, // Configuration + pmFPAfileType type // Expected type + ) +{ + pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, name); + if (!file) { + psError(PS_ERR_UNKNOWN, false, "Unable to define output convolved file %s", name); + return NULL; + } + if (file->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSTACK.OUTCONV is not of type %s", pmFPAfileStringFromType(type)); + return NULL; + } + + return file; +} + +// Define an input convolved image file, using the output as a basis +static pmFPAfile *defineInputConvolved(const char *inputName, // Input FPA file name + pmFPAfile *outFile, // Corresponding output FPA file + pmConfig *config, // Configuration + pmFPAfileType type // Expected type + ) +{ + pmFPAview *view = pmFPAviewAlloc(0); // View into sky cells + view->chip = view->cell = view->readout = 0; + + psString imageName = pmFPANameFromRule(outFile->filerule, outFile->fpa, view); + psArray *imageNames = psArrayAlloc(1); + imageNames->data[0] = imageName; + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INCONV.FILENAMES", PS_META_REPLACE, + "Filenames of input convolved image files", imageNames); + psFree(imageNames); + bool found = false; // Found the file? + pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPSTACK.INCONV", + "INCONV.FILENAMES"); + psMetadataRemoveKey(config->arguments, "INCONV.FILENAMES"); + if (!imageFile || !found) { + psError(PS_ERR_UNKNOWN, false, "Unable to define %s file", inputName); + return NULL; + } + if (imageFile->type != type) { + psError(PS_ERR_IO, true, "PPSTACK.INCONV is not of type %s", + pmFPAfileStringFromType(type)); + return NULL; + } + + return imageFile; +} + + + bool ppStackCamera(pmConfig *config) @@ -138,4 +195,27 @@ } } + + // Output convolved files + pmFPAfile *outconvImage = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config, + PM_FPA_FILE_IMAGE); + pmFPAfile *outconvMask = defineOutputConvolved("PPSTACK.OUTCONV.MASK", imageFile->fpa, config, + PM_FPA_FILE_MASK); + pmFPAfile *outconvWeight = defineOutputConvolved("PPSTACK.OUTCONV.WEIGHT", imageFile->fpa, config, + PM_FPA_FILE_WEIGHT); + if (!outconvImage || !outconvMask || !outconvWeight) { + return false; + } + + // Input convolved files + pmFPAfile *inconvImage = defineInputConvolved("PPSTACK.INCONV", outconvImage, config, + PM_FPA_FILE_IMAGE); + pmFPAfile *inconvMask = defineInputConvolved("PPSTACK.INCONV.MASK", outconvMask, config, + PM_FPA_FILE_MASK); + pmFPAfile *inconvWeight = defineInputConvolved("PPSTACK.INCONV.WEIGHT", outconvWeight, config, + PM_FPA_FILE_WEIGHT); + if (!inconvImage || !inconvMask || !inconvWeight) { + return false; + } + psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, Index: /branches/pap_branch_080207/ppStack/src/ppStackLoop.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackLoop.c (revision 16366) +++ /branches/pap_branch_080207/ppStack/src/ppStackLoop.c (revision 16367) @@ -10,4 +10,149 @@ #include "ppStack.h" + +// Here follows lists of files for activation/deactivation at various stages. Each must be NULL-terminated. + +// All files in the system +static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", + "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT", + "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", + "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", + "PSPHOT.PSF.SAVE", "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", + "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", "PSPHOT.BACKMDL.STDEV", + "PSPHOT.BACKGND", "PSPHOT.BACKSUB", "SOURCE.PLOT.MOMENTS", + "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", "PSPHOT.INPUT.CMF", + 0 }; + +// Files required in preparation for convolution +static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 }; + +// Files for convolution; these should be turned on one by one +static char *convolutionFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", + "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT", + 0 }; + +// Input files for the combination +static char *inCombineFiles = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 }; + +// Output files for the combination and photometry +static char *outCombineFiles = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", + "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", + "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB", + "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", + "PSPHOT.INPUT.CMF", 0 }; + + +// Activate/deactivate a list of files +static void fileActivation(pmConfig *config, // Configuration + char **files, // Files to turn on/off + bool state // Activation state + ) +{ + assert(config); + assert(files); + + for (int i = 0; files[i] != NULL; i++) { + pmFPAfileActivate(config->files, state, files[i]); + } + return; +} + +// Activate/deactivate a single element for a list; return array of files +static psArray *fileActivationSingle(pmConfig *config, // Configuration + char **files, // Files to turn on/off + bool state, // Activation state + int num // Number of file in sequence + ) +{ + assert(config); + assert(files); + + psList *list = psListAlloc(NULL); // List of files + + for (int i = 0; files[i] != NULL; i++) { + pmFPAfile *file = pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file + psListAdd(list, PS_LIST_TAIL, file); + } + + psArray *array = psListToArray(list); + psFree(list); + + return array; +} + +// Set the data level for a list of files +static void fileSetDataLevel(pmConfig *config, // Configuration + char **files, // Files for which to set level + pmFPALevel level // Level to set + ) +{ + assert(config); + assert(files); + + for (int i = 0; files[i] != NULL; i++) { + psArray *selected = pmFPAfileSelect(config->files, files[i]); // Selected files of interest + for (int j = 0; j < selected->n; j++) { + pmFPAfile *file = selected->data[j]; + assert(file); + file->dataLevel = level; + } + psFree(selected); + } + return; +} + + + +// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells +static pmFPAview *filesIterateDown(pmConfig *config // Configuration + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + return view; +} + +// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells +static bool filesIterateUp(pmConfig *config // Configuration + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + view->chip = view->cell = view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + return true; +} + bool ppStackLoop(pmConfig *config) @@ -32,79 +177,217 @@ } - pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Token input file - if (!input) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n"); - return false; - } pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file if (!output) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n"); + psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!"); return false; } - pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy - pmHDU *lastHDU = NULL; // Last HDU that was updated - - // Iterate over the FPA hierarchy - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return false; - } - - pmChip *chip; // Chip of interest - while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return false; - } - - pmCell *cell; // Cell of interest - while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + // Preparation iteration: Load the sources, and get a target PSF model + psArray *sources = NULL; // Sources to use for PSF matching + pmPSF *targetPSF = NULL; // Target PSF + { + pmFPAfileActivate(config->files, false, NULL); + fileActivation(config, prepareFiles, true); + pmFPAview *view = filesIterateDown(config); + if (!view) { + return false; + } + + pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Readout + // We want to hang on to the 'sources' even when its host FPA is blown away + sources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES")); + if (!sources) { + psError(PS_ERR_UNKNOWN, true, "Unable to find sources."); + psFree(view); + return false; + } + + // Generate list of PSFs + psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, + "^PPSTACK.INPUT$"); // Iterator + psMetadataItem *fileItem; // Item from iteration + int fileNum = 0; // Number of file + psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope + while ((fileItem = psMetadataGetAndIncrement(fileIter))) { + assert(fileItem->type == PS_DATA_UNKNOWN); + pmFPAfile *inputFile = fileItem->data.V; // An input file + pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF + pmPSF *psf = psMetadataLookupPtr(NULL, ro->parent->parent->analysis, "PSPHOT.PSF"); // PSF + if (!psf) { + psError(PS_ERR_UNKNOWN, false, "Unable to find PSF."); + psFree(view); + psFree(sources); return false; } - - // Put version information into the header - pmHDU *hdu = pmHDUFromCell(cell); - if (hdu && hdu != lastHDU) { - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - ppStackVersionMetadata(hdu->header); - lastHDU = hdu; - } - - pmReadout *readout; // Readout of interest - while ((readout = pmFPAviewNextReadout(view, input->fpa, 1))) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return false; - } - - // Perform the analysis - if (!ppStackReadout(config, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); - return false; - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return false; + psListAdd(psfList, PS_LIST_TAIL, psf); + } + psFree(fileIter); + psFree(view); + + psArray *psfArray = psListToArray(psfList); // Array of PSFs + psFree(psfList); + targetPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel, + psfOrder, psfOrder); + psFree(psfArray); + if (!targetPSF) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); + psFree(sources); + return false; + } + + filesIterateUp(config); + } + + + // Generate convolutions and write them to disk + psArray *cells = psArrayAlloc(numFiles); // Cells for convolved images --- a handle for reading again + for (int i = 0; i < numFiles; i++) { + pmFPAfileActivate(config->files, false, NULL); + psArray *files = fileActivationSingle(config, prepareFiles, true, i); + pmFPAview *view = filesIterateDown(config); + if (!view) { + return false; + } + + pmReadout *readout = pmFPAviewThisReadout(view, files->data[0]); // Input readout + psFree(view); + + // Background subtraction, scaling and normalisation is performed automatically by the image matching + if (!ppStackMatch(ro, sources, outPSF, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i); + psFree(convolved); + return false; + } + +#ifdef CONVOLUTION_FILES + if (readout->image) { + psString name = NULL; // Name of image + psStringAppend(&name, "convolved%03d_image.fits", i); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, readout->image, 0, NULL); + psFitsClose(fits); + } + if (readout->mask) { + psString name = NULL; // Name of image + psStringAppend(&name, "convolved%03d_mask.fits", i); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, readout->mask, 0, NULL); + psFitsClose(fits); + } + if (readout->weight) { + psString name = NULL; // Name of image + psStringAppend(&name, "convolved%03d_weight.fits", i); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, readout->weight, 0, NULL); + psFitsClose(fits); + } +#endif // CONVOLUTION_FILES + + cells->data[i] = psMemIncrRefCounter(readout->parent); + filesIterateUp(config); + } + + // Set files to read at the readout level --- then when we iterate down, the pmFPAfileIO stuff should take + // care of the basic things, and we can just worry about grabbing the readouts by chunks + fileSetDataLevel(config, inCombineFiles, PM_FPA_LEVEL_CHUNK); + pmFPAfileActivate(config->files, false, NULL); + fileActivation(inCombineFiles, true); + fileActivation(outCombineFiles, true); + + { + pmFPAview *view = filesIterateDown(config); + if (!view) { + psFree(cells); + return false; + } + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.OUTPUT"); // Output readout + psFree(view); + + psArray *readouts = psArrayAlloc(numFiles); // Readouts for convolved images + for (int i = 0; i < numFiles; i++) { + readouts->data[i] = pmReadoutAlloc(cells->data[i]); // Readout into which to read + } + + // Read convolutions by chunks + int numRead = 0; // Number of inputs read + int numChunk = 0; // Chunk number + bool more = true; // More to read? + for (int numChunk = 0; more; numChunk++) { + for (int i = 0; i < numFiles; i++) { + pmReadout *readout = readouts->data[i]; + assert(readout); + + // Files, containing the FITS handle + pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i); + pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i); + pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i); + + // FITS handles + psFits *fitsImage = imageFile->fits; + psFits *fitsMask = maskFile->fits; + psFits *fitsWeight = weightFile->fits; + + if (!pmReadoutReadChunk(readout, fitsImage, 0, numScans, overlap) || + !pmReadoutReadChunkMask(readout, fitsMask, 0, numScans, overlap) || + !pmReadoutReadChunkWeight(readout, fitsWeight, 0, numScans, overlap)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i); + psFree(readouts); + psFree(cells); } } - // Perform statistics on the cell - if (stats) { - ppStatsFPA(stats, output->fpa, view, maskBlank, config); - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + if (!ppStackReadout(config, outRO, readouts)) { + psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); return false; } - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; - } - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; + + for (int i = 0; i < numFiles && more; i++) { + pmReadout *readout = readouts->data[i]; + assert(readout); + + // Files, containing the FITS handle + pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i); + pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i); + pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i); + + // FITS handles + psFits *fitsImage = imageFile->fits; + psFits *fitsMask = maskFile->fits; + psFits *fitsWeight = weightFile->fits; + + more &= pmReadoutMore(readout, fitsImage, 0, numScans, overlap); + more &= pmReadoutMoreMask(readout, fitsMask, 0, numScans, overlap); + more &= pmReadoutMoreWeight(readout, fitsWeight, 0, numScans, overlap); + } + } + + // Get rid of the input files + fileActivation(outCombineFiles, false); + filesIterateUp(config); + + ppStackPhotometry(config, outRO); + + // Put version information into the header + pmHDU *hdu = pmHDUFromCell(outRO->parent); + if (!hdu) { + psError(PS_ERR_UNKNOWN, true, "Unable to find HDU for output."); + return false; + } + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + ppStackVersionMetadata(hdu->header); + + // Statistics on output + if (stats) { + ppStatsFPA(stats, output->fpa, view, maskBlank, config); + } + + // Write out the output files + filesActivation(outCombineFiles, true); + filesIterateUp(config); }