Index: /branches/czw_branch/20120906/pswarp/src/pswarp.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarp.c (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarp.c (revision 34411) @@ -46,4 +46,10 @@ goto DIE; } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + fprintf(stderr,"Here!\n"); + if (!pswarpDefineBackground(config)) { + goto DIE; + } + } // Open the statistics file @@ -66,4 +72,11 @@ if (!pswarpLoop(config, stats)) { goto DIE; + } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + fprintf(stderr,"There!\n"); + if (!pswarpLoopBackground(config, stats)) { + fprintf(stderr,"Dying!\n"); + goto DIE; + } } Index: /branches/czw_branch/20120906/pswarp/src/pswarp.h =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarp.h (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarp.h (revision 34411) @@ -81,5 +81,7 @@ bool pswarpParseCamera (pmConfig *config); bool pswarpDefine (pmConfig *config); +bool pswarpDefineBackground (pmConfig *config); bool pswarpLoop (pmConfig *config, psMetadata *stats); +bool pswarpLoopBackground (pmConfig *config, psMetadata *stats); psExit pswarpExitCode(psExit exitValue); bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config); Index: /branches/czw_branch/20120906/pswarp/src/pswarpArguments.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpArguments.c (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarpArguments.c (revision 34411) @@ -86,5 +86,6 @@ pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist"); - + pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist"); + if (argc != 3) { usage(); @@ -164,4 +165,21 @@ } + bool doBKG = psMetadataLookupBool(&status,recipe, "BACKGROUND.MODEL"); ///< Generate the warped background model? + if (!status) { + doBKG = false; + psWarning("BACKGROUND.MODEL is not set in the %s recipe -- defaulting to FALSE.", PSWARP_RECIPE); + } + int bkgXgrid = psMetadataLookupBool(&status,recipe, "BKG.XGRID"); ///< Xsize of background model + if (!status) { + bkgXgrid = 10; + psWarning("BKG.XGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgXgrid); + } + int bkgYgrid = psMetadataLookupBool(&status,recipe, "BKG.YGRID"); ///< Xsize of background model + if (!status) { + bkgYgrid = 10; + psWarning("BKG.YGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgYgrid); + } + + // Set recipe values in the recipe (since we've possibly altered some) psMetadataAddS32(recipe, PS_LIST_TAIL, "GRID.NX", PS_META_REPLACE, @@ -176,5 +194,9 @@ "Fraction of bad flux for a pixel to be marked as poor", poorFrac); psMetadataAddBool(recipe, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF); - + psMetadataAddBool(recipe, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG); + psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid); + psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid); + + // Set recipe values in the arguments psMetadataAddS32(config->arguments, PS_LIST_TAIL, "GRID.NX", 0, @@ -189,4 +211,7 @@ "Fraction of bad flux for a pixel to be marked as poor", poorFrac); psMetadataAddBool(config->arguments, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid); psTrace("pswarp", 1, "Done with pswarpArguments...\n"); Index: /branches/czw_branch/20120906/pswarp/src/pswarpDefine.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpDefine.c (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarpDefine.c (revision 34411) @@ -41,5 +41,5 @@ return false; } - + // open the full skycell file; no need to defer different depths. only load the header data pmFPAview *view = pmFPAviewAlloc (0); @@ -86,4 +86,6 @@ psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0"); psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0"); + + } @@ -123,9 +125,153 @@ } } - + view->chip = view->cell = view->readout = -1; pmFPAAddSourceFromView(output->fpa, view, output->format); + psFree (view); return true; } + +bool pswarpDefineBackground (pmConfig *config) { + + // load the PSWARP recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n"); + return false; + } + + // select the input data sources + pmFPAfile *skycell = psMetadataLookupPtr (NULL, config->files, "PSWARP.SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_CONFIG, false, "Can't find skycell data!\n"); + return false; + } + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT.BKGMODEL"); + if (!output) { + psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n"); + return false; + } + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT"); + if (!input) { + psError(PSWARP_ERR_CONFIG, false, "Can't find input data!\n"); + return false; + } + + // open the full skycell file; no need to defer different depths. only load the header data + pmFPAview *view = pmFPAviewAlloc (0); + pmFPAfileOpen (skycell, view, config); + + // Read header and create target + { + if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) { + psError(psErrorCodeLast(), false, "Unable to read headers for skycell."); + psFree(view); + return false; + } + view->chip = 0; + view->cell = 0; + view->readout = 0; + pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell + pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source + if (!hdu || !hdu->header) { + psError(PM_ERR_PROG, false, "Unable to find header for sky cell."); + psFree(view); + return false; + } + int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns + int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows + if ((numCols == 0) || (numRows == 0)) { + psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows); + psFree(view); + return false; + } + + pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell + pmReadout *readout = pmReadoutAlloc(target); ///< Target readout + readout->image = psImageAlloc(numCols / output->xBin, numRows / output->yBin, PS_TYPE_F32); + psImageInit(readout->image, NAN); + psFree(readout); // Drop reference + + bool status = false; + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN"); +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XBIN",PS_META_REPLACE,"",output->xBin); */ +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YBIN",PS_META_REPLACE,"",output->yBin); */ + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE"); +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XSIZE",PS_META_REPLACE,"",numCols / output->xBin); */ +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YSIZE",PS_META_REPLACE,"",numRows / output->yBin); */ + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0"); + + + } + + // XXX this is not a sufficient test + view->chip = 0; + view->cell = 0; + view->readout = -1; + pmHDU *phu = pmFPAviewThisPHU(view, skycell->fpa); ///< Skycell PHU + pmHDU *hdu = pmFPAviewThisHDU(view, skycell->fpa); ///< Skycell header + + pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); + + double cd1f = 1.0 * output->xBin; + double cd2f = 1.0 * output->yBin; + + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + + WCS->trans->x->coeff[1][0] *= cd1f; + WCS->trans->x->coeff[0][1] *= cd2f; + WCS->trans->y->coeff[1][0] *= cd1f; + WCS->trans->y->coeff[0][1] *= cd2f; + + + pmAstromWCStoHeader (hdu->header,WCS); + + bool bilevelAstrometry = false; + if (phu) { + char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp(&ctype[4], "-DIS"); + } + } + + // We read from the skycell into the output. i.e., the output receives the desired astrometry. + pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output + if (bilevelAstrometry) { + if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell."); + psFree(view); + return false; + } + if (!pmAstromReadBilevelChip(outputChip, hdu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell."); + psFree(view); + return false; + } + } else { + // we use a default FPA pixel scale of 1.0 + if (!pmAstromReadWCS(output->fpa, outputChip, hdu->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell."); + psFree(view); + return false; + } + } + + + view->chip = view->cell = view->readout = -1; + pmFPAAddSourceFromView(output->fpa, view, output->format); + + + psFree (view); + return true; +} Index: /branches/czw_branch/20120906/pswarp/src/pswarpFileNames.h =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpFileNames.h (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarpFileNames.h (revision 34411) @@ -12,4 +12,5 @@ "PSWARP.MASK", "PSWARP.VARIANCE", + "PSWARP.BKGMODEL", NULL }; @@ -20,4 +21,5 @@ "PSWARP.OUTPUT.MASK", "PSWARP.OUTPUT.VARIANCE", + "PSWARP.OUTPUT.BKGMODEL", NULL }; Index: /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c (revision 34411) @@ -179,4 +179,6 @@ // read WCS data from the corresponding header pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); + + if (bilevelAstrometry) { if (!pmAstromReadBilevelChip (chip, hdu->header)) { @@ -196,4 +198,5 @@ } + pmCell *cell; while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { @@ -227,5 +230,5 @@ pswarpTransformReadout(output, readout, config); - + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { psError(psErrorCodeLast(), false, "Unable to write files."); @@ -255,5 +258,5 @@ goto DONE; } - + pmCell *outCell = output->parent; ///< Output cell pmChip *outChip = outCell->parent; ///< Output chip @@ -450,2 +453,469 @@ 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) { + psTrace ("pswarp", 4, "AChip %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, "ACell %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"); + } + +/* // 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"); */ + + // 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); + 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; + } + + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); + psMetadataAddS32(config->arguments,PS_LIST_TAIL, + "INTERPOLATION.MODE",PS_META_REPLACE,"", PS_INTERPOLATE_BILINEAR); ///< Mode + psMetadataAddS32(config->arguments,PS_LIST_TAIL, + "INTERPOLATION.NUM",PS_META_REPLACE,"", 64); ///< Mode + + 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; + } + } +#define CORRECT_INPUT_WCS 1 +#if CORRECT_INPUT_WCS + // Correct the input WCS + pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); + + double cd1f = 1.0 * 400; + double cd2f = 1.0 * 400; + + WCS->crpix1 = WCS->crpix1 / cd1f + 0.5; + WCS->crpix2 = WCS->crpix2 / cd2f - 0.5; + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + + WCS->trans->x->coeff[1][0] *= cd1f; + WCS->trans->x->coeff[0][1] *= cd2f; + WCS->trans->y->coeff[1][0] *= cd1f; + WCS->trans->y->coeff[0][1] *= cd2f; + + + pmAstromWCStoHeader (hdu->header,WCS); + // End WCS work. +#endif + + 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; + } + + // 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); + } + + 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, "", 2); +/* psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NX",PS_META_REPLACE,"",100); ///< Number of grid points in x */ +/* psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NY",PS_META_REPLACE,"",100); ///< Number of grid points in y */ + + 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; + } + } + + 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; +} Index: /branches/czw_branch/20120906/pswarp/src/pswarpParseCamera.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpParseCamera.c (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarpParseCamera.c (revision 34411) @@ -84,4 +84,10 @@ } + pmFPAfile *inBackground = defineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND", + PM_FPA_FILE_IMAGE); + if (!inBackground) { + psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied"); + } + // The input skycell is a required argument: it defines the output image // XXX we may need a different skycell structure here @@ -134,4 +140,20 @@ } + if (inBackground && psMetadataLookupBool(&mdok, recipe, "BACKGROUND.MODEL")) { + pmFPAfile *outBackground = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL"); +/* pmFPAfile *outBackground = pmFPAfileDefineFromFPA(config,output->fpa, */ +/* psMetadataLookupS32(&mdok,recipe,"BKG.XGRID"), */ +/* psMetadataLookupS32(&mdok,recipe,"BKG.YGRID"), */ +/* "PSWARP.OUTPUT.BKGMODEL"); */ + outBackground->xBin = psMetadataLookupS32(&mdok, recipe, "BKG.XGRID"); + outBackground->yBin = psMetadataLookupS32(&mdok, recipe, "BKG.YGRID"); + + if (!outBackground) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL"); + return false; + } + outBackground->save = true; + } + if (psMetadataLookupBool(&mdok, recipe, "PSF")) { // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by Index: /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c (revision 34411) @@ -29,4 +29,5 @@ psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE"); ///< Mode for interp + int numKernels = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.NUM"); ///< Number of kernels Index: /branches/czw_branch/20120906/pswarp/src/pswarpTransformSources.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpTransformSources.c (revision 34410) +++ /branches/czw_branch/20120906/pswarp/src/pswarpTransformSources.c (revision 34411) @@ -41,5 +41,5 @@ if (!outDetections) { outDetections = pmDetectionsAlloc(); - psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Warped sources", outDetections); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY , "Warped sources", outDetections); } psArray *outSources = outDetections->allSources;