Index: trunk/pswarp/src/pswarp.c =================================================================== --- trunk/pswarp/src/pswarp.c (revision 34089) +++ trunk/pswarp/src/pswarp.c (revision 34800) @@ -46,4 +46,9 @@ goto DIE; } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + if (!pswarpDefineBackground(config)) { + goto DIE; + } + } // Open the statistics file @@ -66,4 +71,10 @@ if (!pswarpLoop(config, stats)) { goto DIE; + } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + if (!pswarpLoopBackground(config, stats)) { + fprintf(stderr,"Dying!\n"); + goto DIE; + } } Index: trunk/pswarp/src/pswarp.h =================================================================== --- trunk/pswarp/src/pswarp.h (revision 34089) +++ trunk/pswarp/src/pswarp.h (revision 34800) @@ -62,4 +62,9 @@ psImage *region; + /** values which are needed to control the background model warping. */ + bool background_warping; + double offset_x; + double offset_y; + /** input values for this tile */ int gridX; @@ -81,5 +86,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: trunk/pswarp/src/pswarpArguments.c =================================================================== --- trunk/pswarp/src/pswarpArguments.c (revision 34089) +++ trunk/pswarp/src/pswarpArguments.c (revision 34800) @@ -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 = psMetadataLookupS32(&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 = psMetadataLookupS32(&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: trunk/pswarp/src/pswarpDefine.c =================================================================== --- trunk/pswarp/src/pswarpDefine.c (revision 34089) +++ trunk/pswarp/src/pswarpDefine.c (revision 34800) @@ -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: trunk/pswarp/src/pswarpFileNames.h =================================================================== --- trunk/pswarp/src/pswarpFileNames.h (revision 34089) +++ trunk/pswarp/src/pswarpFileNames.h (revision 34800) @@ -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: trunk/pswarp/src/pswarpLoop.c =================================================================== --- trunk/pswarp/src/pswarpLoop.c (revision 34089) +++ trunk/pswarp/src/pswarpLoop.c (revision 34800) @@ -179,4 +179,6 @@ // read WCS data from the corresponding header pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); + + if (bilevelAstrometry) { if (!pmAstromReadBilevelChip (chip, hdu->header)) { @@ -195,5 +197,5 @@ } } - + pmCell *cell; while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { @@ -227,5 +229,5 @@ pswarpTransformReadout(output, readout, config); - + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { psError(psErrorCodeLast(), false, "Unable to write files."); @@ -255,5 +257,5 @@ goto DONE; } - + pmCell *outCell = output->parent; ///< Output cell pmChip *outChip = outCell->parent; ///< Output chip @@ -450,2 +452,354 @@ 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; +} Index: trunk/pswarp/src/pswarpParseCamera.c =================================================================== --- trunk/pswarp/src/pswarpParseCamera.c (revision 34089) +++ trunk/pswarp/src/pswarpParseCamera.c (revision 34800) @@ -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: trunk/pswarp/src/pswarpTransformReadout.c =================================================================== --- trunk/pswarp/src/pswarpTransformReadout.c (revision 34089) +++ trunk/pswarp/src/pswarpTransformReadout.c (revision 34800) @@ -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 @@ -60,5 +61,6 @@ // output coordinates to input coordinates pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY); - + // if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { + // XXX optionally modify the grid based on this result and force the maxError < XXX double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid @@ -130,4 +132,11 @@ args->goodPixels = 0; + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { + args->background_warping = true; + args->offset_x = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"); + args->offset_y = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"); + } + + // allocate a job psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE"); @@ -171,5 +180,4 @@ } else { pswarpTransformTileArgs *args = job->args->data[0]; - // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n); goodPixels += args->goodPixels; xMin = PS_MIN(args->xMin, xMin); @@ -183,4 +191,6 @@ jacobian += args->jacobian * args->goodPixels; } + + } psFree(job); @@ -201,8 +211,10 @@ if (goodPixels > 0 && psMetadataLookupBool(&mdok, recipe, "SOURCES")) { + if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { if (!pswarpTransformSources(output, input, config)) { - psError(psErrorCodeLast(), false, "Unable to interpolate image."); - return false; + psError(psErrorCodeLast(), false, "Unable to interpolate image."); + return false; } + } } Index: trunk/pswarp/src/pswarpTransformSources.c =================================================================== --- trunk/pswarp/src/pswarpTransformSources.c (revision 34089) +++ trunk/pswarp/src/pswarpTransformSources.c (revision 34800) @@ -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; Index: trunk/pswarp/src/pswarpTransformTile.c =================================================================== --- trunk/pswarp/src/pswarpTransformTile.c (revision 34089) +++ trunk/pswarp/src/pswarpTransformTile.c (revision 34800) @@ -46,4 +46,8 @@ args->jacobian = NAN; + args->background_warping = false; + args->offset_x = 0.0; + args->offset_y = 0.0; + return args; } @@ -83,16 +87,32 @@ // Iterate over the output image pixels (parent frame) long goodPixels = 0; ///< Number of input pixels landing on the output image + for (int y = yMin; y < yMax; y++) { for (int x = xMin; x < xMax; x++) { - // Only transform those pixels requested if (region && region->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) { continue; } - // pswarpMapApply converts the output coordinate (x,y) to the input coordinate. // both are in the parent frames of the input and output images. double xIn, yIn; // Input pixel coordinates pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5); + + // This needs to use a more reliable method to do this offset and limiting + // if (args->interp->mode == 8) { + if (args->background_warping) { + // double xOffset = 177.0 / 400.0; // (modelsize * modelbinning - xsize) / 2.0 + // double yOffset = 166.0 / 400.0; // (modelsize * modelbinning - ysize) / 2.0 + xIn += args->offset_x; + yIn += args->offset_y; + + if ((xIn > inNumCols - args->offset_x)|| + (yIn > inNumRows - args->offset_y)|| + (xIn < args->offset_x)|| + (yIn < args->offset_y)) { + continue; + } + } + if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) { continue;