Index: trunk/pswarp/src/pswarpLoop.c =================================================================== --- trunk/pswarp/src/pswarpLoop.c (revision 34800) +++ trunk/pswarp/src/pswarpLoop.c (revision 35563) @@ -1,6 +1,5 @@ /** @file pswarpLoop.c * - * @brief - * + * @brief main processing loop for pswarp * @ingroup pswarp * @@ -12,24 +11,11 @@ #include "pswarp.h" -#include -#include "pswarpFileNames.h" // Lists of file rules used at different stages -#define WCS_NONLIN_TOL 0.001 // Non-linear tolerance for header WCS -#define TESTING 0 // Testing output? - - - -// Loop over the inputs, warp them to the output skycell and then write out the output. +// Loop over the inputs, warp them to the output skycell and then update metadata bool pswarpLoop(pmConfig *config, psMetadata *stats) { - bool status; - bool mdok; // Status of MD lookup - - const char *skyCamera = psMetadataLookupStr(NULL, config->arguments, - "SKYCELL.CAMERA"); // Name of camera for skycell - pmConfigCamerasCull(config, skyCamera); - pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); // load the recipe + bool status = false; psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE); if (!recipe) { @@ -38,768 +24,170 @@ } - if (!pswarpSetMaskBits(config)) { - psError(psErrorCodeLast(), false, "failed to set mask bits"); - return NULL; - } - - // output mask bits - psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT"); - psAssert (status, "MASK.OUTPUT was not defined"); - // select the input data sources - pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT"); - if (!input) { - psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n"); + pmFPAfile *output = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT"); + if (!output) { + psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n"); return false; } // use the external astrometry source if supplied - pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM"); - if (!astrom) { - astrom = input; - } - - if (astrom->camera != input->camera) { - psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match."); + pmFPAfile *skycell = psMetadataLookupPtr(&status, config->files, "PSWARP.SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry."); return false; } - // select the output readout pmFPAview *view = pmFPAviewAlloc(0); - view->chip = 0; - view->cell = 0; - view->readout = 0; - pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT"); - if (!output) { - psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n"); + + int nInputs = psMetadataLookupS32(&status, config->arguments, "NUM_INPUTS"); + if (!status) { + psError(PSWARP_ERR_DATA, true, "number of inputs is not defined (programming error)"); return false; } - psFree (view); + // load in the input pixel data (ex. background model) + pmFPAfileActivate(config->files, false, NULL); + pmFPAfileActivate(config->files, true, "PSWARP.INPUT"); + pmFPAfileActivate(config->files, true, "PSWARP.MASK"); + pmFPAfileActivate(config->files, true, "PSWARP.VARIANCE"); - // Turn all skycell files on to generate them, and then turn them off for the loop over the input images - // the input, which is in a different format. - { - pswarpFileActivation(config, detectorFiles, false); - pswarpFileActivation(config, photFiles, false); - pswarpFileActivation(config, independentFiles, false); - pswarpFileActivation(config, skycellFiles, true); - if (!pswarpIOChecksBefore(config)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - pswarpFileActivation(config, skycellFiles, false); + // We re-activate the CMF load so we can transform the sources as well as the pixels. + // We only need to read in these if the astrometry source is CMF. + pmFPAfileActivate(config->files, true, "PSWARP.ASTROM"); + + // loop over this section once per input group + for (int i = 0; i < nInputs; i++) { + // select the input data sources + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSWARP.INPUT", i); + if (!input) { + psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n"); + return false; + } + + // select the input data sources + pmFPAfile *astrom = pmFPAfileSelectSingle(config->files, "PSWARP.ASTROM", i); + if (!astrom) { + astrom = input; + } + + pmFPAviewReset (view); + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto FAIL; + } + + // *** main transformation block + // *** this section loops over the input chips/cells and reads them one at a time + // *** the output chips/cells are filled where appropriate, but not yet written to disk + pmChip *chip; + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto FAIL; + } + + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto FAIL; + } + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto FAIL; + } + if (!readout->data_exists) { + continue; + } + + // Copy the detections from the astrometry carrier to the input, so they can be accessed by + // pswarpTransformReadout + if (astrom != input) { + pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry + pmDetections *detections = psMetadataLookupPtr(&status, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry + if (detections) { + psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections); + } + } + + pswarpTransformToTarget(output->fpa, readout, config, false); + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto FAIL; + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto FAIL; + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto FAIL; + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto FAIL; + } + + if (!pswarpUpdateStatistics (output->fpa, stats, input->fpa, astrom->fpa, config)) { + psError(psErrorCodeLast(), false, "problem generating statistics."); + goto FAIL; + } + if (!pswarpUpdateMetadata (output->fpa, skycell->fpa, input->fpa, astrom->fpa, config, true)) { + psError(psErrorCodeLast(), false, "problem generating statistics."); + goto FAIL; + } } - // Read the input astrometry - // XXX rather than use the activations here, this should just explicitly loop over the desired filerule - { - pmFPAfileActivate(config->files, true, "PSWARP.ASTROM"); - - pmChip *chip; - pmFPAview *view = pmFPAviewAlloc(0); - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); - if (!chip->process || !chip->file_exists) { continue; } - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - pmCell *cell; - while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { - psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); - if (!cell->process || !cell->file_exists) { continue; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) || - !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - psFree(view); - - pswarpFileActivation(config, detectorFiles, true); - pmFPAfileActivate(config->files, false, "PSWARP.ASTROM"); + if (!pswarpMakePSF (config, output, stats)) { + psError(psErrorCodeLast(), false, "problem generating PSF."); + goto FAIL; } - // Turn on the source output --- we need to get rid of these so that we can measure the PSF - pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES"); + psFree(view); + return true; - // Don't care about the skycell anymore --- we've read it, and that's all we need to do. - pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL"); - view = pmFPAviewAlloc(0); + FAIL: + psFree (view); + return false; +} - // find the FPA phu - bool bilevelAstrometry = false; - pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa); - if (phu) { - char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); - if (ctype) { - bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } - } - if (bilevelAstrometry) { - if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) { - psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA."); - psFree(view); - psFree(stats); - goto DONE; - } - } +// once the output fpa elements have been built, loop over the fpa and generate stats +// for each readout +bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config, bool backgroundWarp) { - psList *cells = psListAlloc(NULL); // List of cells, for concepts averaging - - // files associated with the science image - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - + pmFPAview *view = pmFPAviewAlloc(0); + pmChip *chip; - while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) { psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists) { continue; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - // read WCS data from the corresponding header - pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); - - - if (bilevelAstrometry) { - if (!pmAstromReadBilevelChip (chip, hdu->header)) { - psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA."); - psFree(view); - psFree(stats); - goto DONE; - } - } else { - // we use a default FPA pixel scale of 1.0 - if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) { - psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); - psFree(view); - psFree(stats); - goto DONE; - } - } - pmCell *cell; - while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + while ((cell = pmFPAviewNextCell (view, output, 1)) != NULL) { psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); if (!cell->process || !cell->file_exists) { continue; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - - psListAdd(cells, PS_LIST_TAIL, cell); // process each of the readouts pmReadout *readout; - while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - if (!readout->data_exists) { - continue; - } - - // Copy the detections from the astrometry carrier to the input, so they can be accessed by - // pswarpTransformReadout - pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry - pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry - if (detections) { - psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections); - } - - pswarpTransformReadout(output, readout, config); - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } + while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) { + pswarpTransformReadout (readout, input, config, backgroundWarp); + } + } } - - if (!output->data_exists) { - psWarning("No overlap between input and skycell."); - if (stats) { - psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, - "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP); - } - psphotFilesActivate(config, false); - psFree(cells); - psFree(view); - goto DONE; - } - - pmCell *outCell = output->parent; ///< Output cell - pmChip *outChip = outCell->parent; ///< Output chip - pmFPA *outFPA = outChip->parent; ///< Output FP - - if (!pswarpPixelsLit(output, stats, config)) { - psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); - psFree(cells); - psFree(view); - goto DONE; - } - bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS"); - if (doStats) { - if (!pswarpMaskStats(output, stats, config)) { - psError(psErrorCodeLast(), false, "Unable to calculate mask stats."); - psFree(cells); - psFree(view); - goto DONE; - } - } - // Set covariance matrix for output - { - psList *covariances = psMetadataLookupPtr(&mdok, output->analysis, - PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices - psAssert(covariances, "Should be there"); - psArray *covars = psListToArray(covariances); // Array of covariance matrices - psKernel *covar = psImageCovarianceAverage(covars); - psFree(covars); - psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES); - - // Correct covariance matrix scale for the mean (square root of the) Jacobian - double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian - int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX); // Good pixels - jacobian /= goodPixels; - output->covariance = psImageCovarianceScale(covar, jacobian); - psFree(covar); - - if (output->variance) { - psImageCovarianceTransfer(output->variance, output->covariance); - } - } - - if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) { - psError(psErrorCodeLast(), false, "Unable to average cell concepts."); - psFree(stats); - psFree(cells); - psFree(view); - goto DONE; - } - psFree(cells); - - psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section - trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels - - if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) { - psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output."); - psFree(stats); - psFree(view); - goto DONE; - } - - // Update ZP from the astrometry - { - psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP"); - item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP"); - } - - pmHDU *hdu = outFPA->hdu; ///< HDU for the output warped image - - // Copy header from target - { - pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell - skyView->chip = skyView->cell = 0; - pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell - psFree(skyView); - pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU - if (!skyHDU) { - psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU."); - psFree(view); - goto DONE; - } - hdu->header = psMetadataCopy(hdu->header, skyHDU->header); - } - - pswarpVersionHeader(hdu->header); - - if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) { - psError(psErrorCodeLast(), false, "Unable to generate WCS header."); - psFree(stats); - goto DONE; - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - - // Done with the detector side of things - pswarpFileActivation(config, detectorFiles, false); - pswarpFileActivation(config, independentFiles, false); - - - // We need a new PSF model for the warped frame. It would be good to generate this analytically, but - // that's going to be tricky. We have a list of sources, so we use those to redetermine the PSF model. - - if (psMetadataLookupBool(&mdok, recipe, "PSF")) { - pswarpFileActivation(config, photFiles, true); - if (!pswarpIOChecksBefore(config)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - - // supply the readout and fpa of interest to psphot - pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); - pmFPACopy(photFile->fpa, outFPA); - - pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell - view->chip = view->cell = view->readout = 0; - - // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF) - psArray *sources = psphotLoadPSFSources (config, view); - if (!sources) { - psError(psErrorCodeLast(), false, "No sources supplied to measure PSF"); - goto DONE; - } - - pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT); - - // measure the PSF using these sources - if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) { - // This is likely a data quality issue - // XXX Split into multiple cases using error codes? - psErrorStackPrint(stderr, "Unable to determine PSF"); - psWarning("Unable to determine PSF --- suspect bad data quality."); - if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) { - psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, - "Unable to determine PSF", psErrorCodeLast()); - } - psErrorClear(); - psphotFilesActivate(config, false); - } - - // Ensure seeing is carried over - pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa); // Chip with seeing - psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing - item->data.F32 = psMetadataLookupF32(NULL, photChip->concepts, "CHIP.SEEING"); - -// XXX EAM : put this in a visualization function -#if (TESTING) - { - #define PSF_SIZE 20 ///< Half-size of PSF - #define PSF_FLUX 10000 ///< Central flux for PSF - pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa); - pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF"); - psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32); - psImageInit(image, 0); - pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX); - pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0); - psFree(model); - psFits *fits = psFitsOpen("psf.fits", "w"); - psFitsWriteImage(fits, NULL, image, 0, NULL); - psFitsClose(fits); - psFree(image); - } -#endif - - psFree(view); - } - - // Perform statistics on the output image - if (stats) { - if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) { - psWarning("Unable to perform statistics on warped image."); - } - } - - - // Add MD5 information for readout - const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME"); - const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME"); - psString headerName = NULL; ///< Header name for MD5 - psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout); - psVector *md5 = psImageMD5(output->image); ///< md5 hash - psString md5string = psMD5toString(md5); ///< String - psFree(md5); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, - "Image MD5", md5string); - psFree(md5string); - psFree(headerName); - psFree(view); - - DONE: - return true; } -// Loop over the inputs, warp them to the output skycell and then write out the output. -bool pswarpLoopBackground(pmConfig *config, psMetadata *stats) -{ - bool status; - bool mdok; // Status of MD lookup - const char *skyCamera = psMetadataLookupStr(NULL, config->arguments, - "SKYCELL.CAMERA"); // Name of camera for skycell - pmConfigCamerasCull(config, skyCamera); - pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); - - // load the recipe - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE); - if (!recipe) { - psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE); - return false; - } - - if (!pswarpSetMaskBits(config)) { - psError(psErrorCodeLast(), false, "failed to set mask bits"); - return NULL; - } - - // select the input data sources - pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.BKGMODEL"); - if (!input) { - psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n"); - return false; - } - - // use the external astrometry source if supplied - pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM"); - if (!astrom) { - astrom = input; - } - - if (astrom->camera != input->camera) { - psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match."); - return false; - } - - // select the output readout - pmFPAview *view = pmFPAviewAlloc(0); - view->chip = 0; - view->cell = 0; - view->readout = 0; - pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT.BKGMODEL"); - if (!output) { - psError(PSWARP_ERR_CONFIG, true, "Can't find output background data!\n"); - return false; - } - psFree (view); - // Turn all skycell files on to generate them, and then turn them off for the loop over the input images - // the input, which is in a different format. - { - pswarpFileActivation(config, detectorFiles, false); - pswarpFileActivation(config, photFiles, false); - pswarpFileActivation(config, independentFiles, false); - pswarpFileActivation(config, skycellFiles, true); - if (!pswarpIOChecksBefore(config)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - pswarpFileActivation(config, skycellFiles, false); - } - // Read the input astrometry - // XXX rather than use the activations here, this should just explicitly loop over the desired filerule - { - - pmFPAfileActivate(config->files, true, "PSWARP.ASTROM"); - - pmChip *chip; - pmFPAview *view = pmFPAviewAlloc(0); - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - - while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { -#if 0 - // This needs to be removed because otherwise it throws an error of duplicate PSPHOT.DETECTIONS. - if (!chip->process || !chip->file_exists) { continue; } - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } -#endif - pmCell *cell; - while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { - psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast()); - if (!cell->process || !cell->file_exists) { continue; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) || - !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - psFree(view); - pswarpFileActivation(config, detectorFiles, true); - pmFPAfileActivate(config->files, false, "PSWARP.ASTROM"); - } - - // Don't care about the skycell anymore --- we've read it, and that's all we need to do. - pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL"); - view = pmFPAviewAlloc(0); - - // find the FPA phu - bool bilevelAstrometry = false; - pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa); - - // pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header); - - if (phu) { - char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); - if (ctype) { - bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } - } - if (bilevelAstrometry) { - if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) { - psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA."); - psFree(view); - psFree(stats); - goto DONE; - } - } - - psList *cells = psListAlloc(NULL); // List of cells, for concepts averaging - - // files associated with the science image - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - pmChip *chip; - while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psTrace ("pswarp", 4, "DChip %d: %x %x\n", view->chip, chip->file_exists, chip->process); - if (!chip->process || !chip->file_exists) { continue; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - - // Pull information from the header of the background files so we can use it to set values. - pmHDU *hdu = pmFPAviewThisHDU(view,input->fpa); - psMetadata *header = hdu->header; - - int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1"); - int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2"); - int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1"); - int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2"); - char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM"); - int CCDSUM1 = atoi(strtok(CCDSUM," ")); - int CCDSUM2 = atoi(strtok(NULL," ")); - - psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE, - "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1)); - psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE, - "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2)); - psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d", - psMetadataLookupS32(NULL,header,"IMNAXIS1"), - psMetadataLookupS32(NULL,header,"IMNAXIS2"), - psMetadataLookupS32(NULL,header,"NAXIS1"), - psMetadataLookupS32(NULL,header,"NAXIS2"), - CCDSUM1,CCDSUM2, - psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"), - psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"), - psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"), - psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID")); - - - // read WCS data from the corresponding header - hdu = pmFPAviewThisHDU (view, astrom->fpa); - - pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); - - double cd1f = (1.0 * CCDSUM1);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID")); - double cd2f = (1.0 * CCDSUM2);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID")); - - WCS->cdelt1 *= cd1f; - WCS->cdelt2 *= cd2f; - WCS->crpix1 = WCS->crpix1 / cd1f; - WCS->crpix2 = WCS->crpix2 / cd2f; - - // WCS->trans->x->nX/nY - for (int q = 0; q <= WCS->trans->x->nX; q++) { - for (int r = 0; r <= WCS->trans->x->nY; r++) { - WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); - } - } - for (int q = 0; q <= WCS->trans->y->nX; q++) { - for (int r = 0; r <= WCS->trans->y->nY; r++) { - WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); - } - } - pmAstromWCStoHeader (hdu->header,WCS); - - if (bilevelAstrometry) { - if (!pmAstromReadBilevelChip (chip, hdu->header)) { - psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA."); - psFree(view); - psFree(stats); - goto DONE; - } - } else { - // we use a default FPA pixel scale of 1.0 - if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) { - psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); - psFree(view); - psFree(stats); - goto DONE; - } - } - // Modify structure here. - - - pmCell *cell; - while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { - psTrace ("pswarp", 4, "DCell %d: %x %x\n", view->cell, cell->file_exists, cell->process); - if (!cell->process || !cell->file_exists) { continue; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - - psListAdd(cells, PS_LIST_TAIL, cell); - - // process each of the readouts - pmReadout *readout; - while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Unable to read files."); - goto DONE; - } - if (!readout->data_exists) { - continue; - } - - for (int x = 0; x < readout->image->numCols; x++) { - for (int y = 0; y < readout->image->numRows; y++) { - readout->image->data.F32[y][x] = readout->image->data.F32[y][x] * (cd1f * cd2f) / - (psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID") * - psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID")); - } - } - - psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8); - - psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true); - pswarpTransformReadout(output, readout, config); - psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false); - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - } - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - } - if (!output->data_exists) { - psWarning("No overlap between input and skycell."); - psphotFilesActivate(config, false); - psFree(cells); - psFree(view); - goto DONE; - } - pmCell *outCell = output->parent; ///< Output cell - pmChip *outChip = outCell->parent; ///< Output chip - pmFPA *outFPA = outChip->parent; ///< Output FP - -/* if (!pswarpPixelsLit(output, stats, config)) { */ -/* psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); */ -/* psFree(cells); */ -/* psFree(view); */ -/* goto DONE; */ -/* } */ - psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section - trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels - - if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) { - psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output."); - psFree(stats); - psFree(view); - goto DONE; - } - - pmHDU *hdu = outFPA->hdu; ///< HDU for the output warped image - - // Copy header from target - { - pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell - skyView->chip = skyView->cell = 0; - pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell - psFree(skyView); - pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU - if (!skyHDU) { - psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU."); - psFree(view); - goto DONE; - } - hdu->header = psMetadataCopy(hdu->header, skyHDU->header); - } - pswarpVersionHeader(hdu->header); - - if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) { - psError(psErrorCodeLast(), false, "Unable to generate WCS header."); - psFree(stats); - goto DONE; - } - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - goto DONE; - } - // Done with the detector side of things - pswarpFileActivation(config, detectorFiles, false); - pswarpFileActivation(config, independentFiles, false); - - - // Add MD5 information for readout - const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME"); - const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME"); - psString headerName = NULL; ///< Header name for MD5 - psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout); - psVector *md5 = psImageMD5(output->image); ///< md5 hash - psString md5string = psMD5toString(md5); ///< String - psFree(md5); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, - "Image MD5", md5string); - psFree(md5string); - psFree(headerName); - psFree(view); - DONE: - - return true; -}