Index: /trunk/pswarp/doc/notes.20130406.txt =================================================================== --- /trunk/pswarp/doc/notes.20130406.txt (revision 35563) +++ /trunk/pswarp/doc/notes.20130406.txt (revision 35563) @@ -0,0 +1,108 @@ + +20130406 + +examining pswarp to handle inverse transformations or more general transformations: + +pswarp: + * setup / init functions (psLibInit, psphotInit, etc) + * pswarpArguments : parse things like -file (name), -list (name), etc + * pswarpParseCamera : examine minimal files, define cameras & pmFPAfile + pmFPA + this function reads headers for files + + uses pmFPAfileDefineFromArgs / pmFPAfileBindFromArgs for inputs + uses pswarpDefineSkycell for outputs + + pmFPAfileDefineSkycell is used to generate the output + FPAs. This forces the file to be FPA-level only, which is wrong. + Not certain of the impact of replacing with pmFPAfileDefineOutput + + pmFPAfileDefineSkycell: + * appends -SKYCELL to the camera name + + * pswarpOptions : check the recipe options, reconcile with cmd-line options + + * pswarpDefine : create output image pixels, update related concepts, + associate output astrometry (skycell) with output pmFPA structures + + + +in pswarpParseCamera: + * pswarpDefineSkycell : this function is based on pmFPAfileDefineFromArgs + * reads fileset from config->arguments (a fileset is an array of filenames) + * select first filename & translate to UNIX namespace + * open fits file + * read fits header (PHU) + * inherit some config elements from input config (user, system, files, arguments) + * determine format from header + (NOTE: currently does not examine or set fileLevel) + * build the pmFPA structure based on camera & name + * bind filerule to fpa, generate the pmFPAfile + * sets file->fileLevel based on format, but slightly out of sync wrt pmFPAfileDefineFromArgs + (NOTE: does not expect / examine any other file) + + +in pswarpLoop (major blocks): + + * find the needed recipe / filerules / etc + * pswarpIOChecksBefore() : generate the outputs? (isn't this done in pswarpDefine?) + * read input astrometry headers (also loads sources or not?) + * apply astrometry headers to input->fpa + * pswarpTransformReadout : transform pixels, etc + * pswarpPixelsLit : measure filling & masks stats + * modify the covariance + * modify the ZP & trimsec (seems like these should be encapsulated) + * update headers + * measure the PSF + * measure pixel stats on output + * add MD5 to header + + + +---- + +pswarp issues: + +1) operations need to loop over the output astrometry / targets for fpa/chip/cell levels + a) I replaced pswarpDefineSkycell with an update to pmFPAfileDefineXX which does the looping + b) is it identical in behavior to pswarpDefineSkycell? + * YES : I can drop pswarpDefineSkycell + +2) how do we load the SOURCES? + + * I need to load the headers for the astrometry, and I do this by + setting the type of the input astrometry to be WCS. is this the + same as CMF, but the latter reads the sources?? + + * by changing the filerules entry CMF <-> WCS pswarpDefine.c either does or does not read the sources + + (for testing, I want to use a skycell image or template to define the astrometry; this would not have sources to be read) + +3) does the I/O loop in pswarpLoop (pswarpIOChecksBefore) do anything useful? + ** I don't think so... this operation is done in pswarpDefine. + +4) how do I determine the overlaps correctly? (I'd like to do this before I load sources or pixels) + + ** pswarpMapSetLocalModel does not seem to depend on astrometry model style (GOOD!) + +----- + +Changes to support arbitrary fpa to arbitrary fpa: + +in pswarpDefine: + * load astrometry headers (source and target) + * determine overlap (source -> target) + * loop over target chips + * determine o + * activate output targets which will get data. + * generate output images + +in pswarpLoop + * need to loop over output targets as well as input sources... + +*** looks like I can select the correct output chips and write them. + Currently, the output format is inherited from the output + astrometry description. that means we get a MEF output file with + one chip per extension. that seems poor, but I think we can fix + this with the filerules. + +** trouble running warp -> chip using astrometry from skycell (ok using astrom from warp) Index: /trunk/pswarp/src/Makefile.am =================================================================== --- /trunk/pswarp/src/Makefile.am (revision 35562) +++ /trunk/pswarp/src/Makefile.am (revision 35563) @@ -28,12 +28,21 @@ pswarp.c \ pswarpArguments.c \ - pswarpExit.c \ + pswarpExitCode.c \ pswarpLoop.c \ + pswarpLoopBackground.c \ + pswarpCleanup.c \ + pswarpStatsFile.c \ pswarpDefine.c \ - pswarpDefineSkycell.c \ + pswarpDefineBackground.c \ + pswarpDefineSkycell.c \ + pswarpDefineLayout.c \ pswarpErrorCodes.c \ + pswarpLoadAstrometry.c \ + pswarpMakePSF.c \ pswarpMapGrid.c \ pswarpMaskStats.c \ pswarpMatchRange.c \ + pswarpOverlaps.c \ + pswarpOptions.c \ pswarpParseCamera.c \ pswarpPixelsLit.c \ @@ -43,4 +52,6 @@ pswarpTransformSources.c \ pswarpTransformTile.c \ + pswarpUpdateStatistics.c \ + pswarpUpdateMetadata.c \ pswarpVersion.c \ pswarpFiles.c Index: /trunk/pswarp/src/pswarp.c =================================================================== --- /trunk/pswarp/src/pswarp.c (revision 35562) +++ /trunk/pswarp/src/pswarp.c (revision 35563) @@ -12,5 +12,4 @@ #include "pswarp.h" -#include "pswarpFileNames.h" int main (int argc, char **argv) @@ -22,12 +21,10 @@ psphotInit(); - const char *statsName = NULL; // Filename for statistics - psMetadata *stats = NULL; // Container for statistics - FILE *statsFile = NULL; // File stream for statistics - pmConfig *config = pswarpArguments(argc, argv); if (!config) { - goto DIE; + pswarpCleanup(config, NULL); } + + pswarpStatsFile *statsFile = pswarpStatsFileOpen (config); pswarpVersionPrint(); @@ -35,111 +32,29 @@ // load identify the data sources if (!pswarpParseCamera(config)) { - goto DIE; + pswarpCleanup(config, statsFile); } if (!pswarpOptions(config)) { - goto DIE; + pswarpCleanup(config, statsFile); } - // load the skycell layout information - if (!pswarpDefine(config)) { - goto DIE; - } - if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { - if (!pswarpDefineBackground(config)) { - goto DIE; - } + // load the input & output astrometry, find the output overlaps, generate the output pixels + if (!pswarpDefineLayout(config)) { + pswarpCleanup(config, statsFile); } - // Open the statistics file - bool mdok; ///< Status of MD lookup - statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); ///< Filename for statistics - if (mdok && statsName && strlen(statsName) > 0) { - psString resolved = pmConfigConvertFilename(statsName, config, true, true); - statsFile = fopen(resolved, "w"); - if (!statsFile) { - psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); - psFree(resolved); - goto DIE; - } - psFree(resolved); - stats = psMetadataAlloc(); - psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0); + // load input pixels & and warp + if (!pswarpLoop(config, statsFile->md)) { + pswarpCleanup(config, statsFile); } - // load and warp - if (!pswarpLoop(config, stats)) { - goto DIE; - } - if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { - if (!pswarpLoopBackground(config, stats)) { - fprintf(stderr,"Dying!\n"); - goto DIE; - } + // load input pixels & and warp + if (!pswarpLoopBackground(config, statsFile->md)) { + pswarpCleanup(config, statsFile); } + // output and free + // NOTE: pswarpCleanup calls exit psLogMsg("pswarp", PS_LOG_INFO, "complete pswarp run: %f sec\n", psTimerMark("pswarp")); - -DIE: - { - psExit exitValue = pswarpExitCode(PS_EXIT_SUCCESS); // Exit code - - // Ensure everything is written out, at every level - pswarpFileActivation(config, detectorFiles, true); - pswarpFileActivation(config, skycellFiles, true); - pswarpFileActivation(config, photFiles, true); - pswarpFileActivation(config, independentFiles, true); - if (!pswarpIOChecksAfter(config)) { - psError(psErrorCodeLast(), false, "Unable to write files."); - exitValue = pswarpExitCode(exitValue); - pmFPAfileFreeSetStrict(false); - } - - // Write out summary statistics - if (stats && statsFile) { - psMetadataAddF32(stats, PS_LIST_TAIL, "DT_WARP", 0, "Time for warp completion", - psTimerMark("pswarp")); - - const char *statsMDC = psMetadataConfigFormat(stats); - if (!statsMDC) { - psError(psErrorCodeLast(), false, "Unable to get statistics file."); - exitValue = pswarpExitCode(exitValue); - } - psFree(stats); - if (fprintf(statsFile, "%s", statsMDC) != strlen(statsMDC)) { - psError(PSWARP_ERR_IO, true, "Unable to write statistics file."); - exitValue = pswarpExitCode(exitValue); - } - psFree(statsMDC); - if (fclose(statsFile) == EOF) { - psError(PSWARP_ERR_IO, true, "Unable to close statistics file."); - exitValue = pswarpExitCode(exitValue); - } - pmConfigRunFilenameAddWrite(config, "STATS", statsName); - } - - // Dump configuration - psString dump_file = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); - if (dump_file) { - if (!pmConfigDump(config, dump_file)) { - psError(psErrorCodeLast(), false, "Unable to dump configuration"); - exitValue = pswarpExitCode(exitValue); - } - } - - psThreadPoolFinalize(); - psMemCheckCorruption(stderr, true); - - psFree(config); - - psTimerStop(); - pmVisualClose(); - pmModelClassCleanup(); - pmConceptsDone(); - pmConfigDone(); - psLibFinalize(); - - exitValue = pswarpExitCode(exitValue); - return exitValue; - } + pswarpCleanup(config, statsFile); } Index: /trunk/pswarp/src/pswarp.h =================================================================== --- /trunk/pswarp/src/pswarp.h (revision 35562) +++ /trunk/pswarp/src/pswarp.h (revision 35563) @@ -22,4 +22,5 @@ #include #include +#include #define THREADED 1 @@ -78,4 +79,10 @@ } pswarpTransformTileArgs; +typedef struct { + FILE *f; // File stream for statistics + char *name; // Filename for statistics + psMetadata *md; // Container for statistics +} pswarpStatsFile; + pswarpTransformTileArgs *pswarpTransformTileArgsAlloc(void); bool pswarpTransformTile (pswarpTransformTileArgs *args); @@ -88,7 +95,8 @@ bool pswarpDefineBackground (pmConfig *config); bool pswarpLoop (pmConfig *config, psMetadata *stats); +bool pswarpLoopSkycell (pmConfig *config, psMetadata *stats); bool pswarpLoopBackground (pmConfig *config, psMetadata *stats); psExit pswarpExitCode(psExit exitValue); -bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config); +bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config, bool backgroundWarp); bool pswarpTransformSources(pmReadout *output, pmReadout *input, pmConfig *config); @@ -160,2 +168,53 @@ ); +pswarpStatsFile *pswarpStatsFileAlloc (); +pswarpStatsFile *pswarpStatsFileOpen (pmConfig *config); +bool pswarpStatsFileSave (pmConfig *config, pswarpStatsFile *statsFile); + +// cleanup memory and exit +void pswarpCleanup (pmConfig *config, pswarpStatsFile *statsFile) PS_ATTR_NORETURN; + +// load the astrometry header info and generate the tranformations +bool pswarpDefineLayout (pmConfig *config); + +// XXX function for testing +bool pswarpDumpOutput (pmConfig *config); + +// structure to describe approximate bounds for each fpa element (eg, chip) +// P,Q are a locally linear projection +typedef struct { + psVector *Pmin; + psVector *Pmax; + psVector *Qmin; + psVector *Qmax; +} pswarpBounds; + +bool pswarpFindOverlap (pmFPA *input, pmFPA *output, pswarpBounds *src, pswarpBounds *tgt); +psProjection *pswarpLocalFrame (pmFPA *fpa); +pswarpBounds *pswarpMakeBounds (pmFPA *fpa, psProjection *frame); +bool pswarpBoundsAppend(pswarpBounds *bounds, float Pmin, float Pmax, float Qmin, float Qmax); +pswarpBounds *pswarpBoundsAlloc(); + +bool pswarpLoadAstrometry (pmFPAfile *target, pmFPAfile *astrom, pmConfig *config); + +bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config, bool backgroundWarp); +bool pswarpMakePSF (pmConfig *config, pmFPAfile *output, psMetadata *stats); +bool pswarpUpdateStatistics (pmFPA *output, psMetadata *stats, pmFPA *input, pmFPA *astrom, pmConfig *config); +bool pswarpUpdateMetadata (pmFPA *output, pmFPA *skycell, pmFPA *input, pmFPA *astrom, pmConfig *config, bool fullImage); + +// XXX functions in pswarpParseCamera + +bool AddStringAsArray (psMetadata *md, char *string, char *name); + +pmFPAfile *pswarpDefineInputFile(pmConfig *config,// Configuration + pmFPAfile *bind, // File to which to bind, or NULL + char *filerule, // Name of file rule + char *argname, // Argument name + pmFPAfileType fileType // Type of file + ); + +bool pswarpParseSingleInput (pmConfig *config); +bool pswarpParseMultiInput (pmConfig *config, psMetadata *fileListMD); + +bool pswarpModifyChipAstrom (pmConfig *config, pmFPAview *view, pmChip *chip, pmFPAfile *astrom, bool bilevelAstrometry, double xBin, double yBin); +bool pswarpGetInputScales (double *xBin, double *yBin, pmConfig *config, pmFPAview *view, pmChip *chip); Index: /trunk/pswarp/src/pswarpArguments.c =================================================================== --- /trunk/pswarp/src/pswarpArguments.c (revision 35562) +++ /trunk/pswarp/src/pswarpArguments.c (revision 35563) @@ -12,10 +12,10 @@ # include "pswarp.h" -# include - static void usage (void) { - fprintf(stderr, "USAGE: pswarp [-file image(s)] [-list imagelist] [options] (output) (skycell)\n"); + fprintf(stderr, "USAGE: pswarp [-input input.mdc] [-file image(s)] [-list imagelist] [options] (output) (skycell)\n"); fprintf(stderr, " options:\n"); + fprintf(stderr, " [-input input.mdc] : input image information in a metadata file\n"); + fprintf(stderr, " [-file input.fits[,input.fits]] : input image to be warped\n"); fprintf(stderr, " [-astrom astrom.cmp] : provide an alternative astrometry calibration\n"); fprintf(stderr, " [-mask mask.fits] : provide a corresponding mask image\n"); @@ -24,5 +24,4 @@ exit(PS_EXIT_CONFIG_ERROR); } - pmConfig *pswarpArguments (int argc, char **argv) { @@ -55,12 +54,11 @@ } + // XXX move to the single group below? pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", "-astromlist"); - { - int arg; ///< Argument Number - if ((arg = psArgumentGet(argc, argv, "-psphot-visual"))) { - psArgumentRemove(arg, &argc, argv); - pmVisualSetVisual(true); - } + // turn on psphot visualization + if ((N = psArgumentGet(argc, argv, "-psphot-visual"))) { + psArgumentRemove(N, &argc, argv); + pmVisualSetVisual(true); } @@ -83,14 +81,50 @@ pswarpSetThreads(); - pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); - pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist"); - pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist"); - pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist"); + // there are three mutually exclusive ways of providing the input + // 1) supply -file (filename) [-mask .. -variance ..] on the command line + // 2) supply -input (input.mdc) on the command line + // 3) load inputs from RUN config info + + // below, we check first for -file then for -input. failure to find either implies use of + // the configuration metadata file. + + bool singleInput = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); + if (singleInput) { + if (psArgumentGet(argc, argv, "-input")) { + psErrorStackPrint(stderr, "error in arguments : -input and -file / -list are mutually exclusive"); + exit(PS_EXIT_CONFIG_ERROR); + } + pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist"); + pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist"); + pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist"); + } else { + // find the input data file (an mdc file) + if ((N = psArgumentGet(argc, argv, "-input"))) { + if (argc <= N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove(N, &argc, argv); + + unsigned int numBad = 0; // Number of bad lines + psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[N], false); // Input file info + if (!inputs || numBad > 0) { + psErrorStackPrint(stderr, "Unable to cleanly read MDC file with inputs."); + exit(PS_EXIT_CONFIG_ERROR); + } + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "INPUTS", 0, "Metadata with input details", inputs); + psFree(inputs); + + psArgumentRemove(N, &argc, argv); + } + } + if (argc != 3) { + usage(); + } + if (psErrorCodeLast() != PS_ERR_NONE) { + psErrorStackPrint(stderr, "error in arguments"); + exit(PS_EXIT_CONFIG_ERROR); + } - if (argc != 3) { - usage(); - } - - psArray *array; @@ -107,114 +141,2 @@ return config; } - -/** - * Parse the recipe and format into the arguments - */ -bool pswarpOptions(pmConfig *config) -{ - // Select the appropriate recipe - psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PSWARP_RECIPE); - if (!recipe) { - psError(PSWARP_ERR_CONFIG, true, "Can't find %s recipe!\n", PSWARP_RECIPE); - return false; - } - - // Get grid size - bool status; ///< Status of MD lookup - int nGridX = psMetadataLookupS32(&status, recipe, "GRID.NX"); - if (!status || nGridX <= 0) { - nGridX = 128; - psWarning("GRID.NX is not set in the recipe --- defaulting to %d", nGridX); - } - int nGridY = psMetadataLookupS32(&status, recipe, "GRID.NY"); - if (!status) { - nGridY = 128; - psWarning("GRID.NY is not set in the recipe --- defaulting to %d", nGridY); - } - - // Get interpolation mode - const char *name = psMetadataLookupStr (&status, recipe, "INTERPOLATION.MODE"); ///< Name of interp mode - if (!name) { - name = "BILINEAR"; - psLogMsg("pswarp", 3, "defaulting to %s interpolation", name); - } - psImageInterpolateMode interpolationMode = psImageInterpolateModeFromString(name); ///< Mode for interp. - if (interpolationMode == PS_INTERPOLATE_NONE) { - interpolationMode = PS_INTERPOLATE_BILINEAR; - psLogMsg ("pswarp", 3, - "Unknown interpolation mode %s, defaulting to bilinear interpolation\n", name); - name = "BILINEAR"; - } - - int numKernels = psMetadataLookupS32(&status, recipe, "INTERPOLATION.NUM"); - if (!status) { - numKernels = 0; - psWarning("INTERPOLATION.NUM is not set in the recipe --- defaulting to %d", numKernels); - } - - float poorFrac = psMetadataLookupF32(&status, recipe, "POOR.FRAC"); ///< Frac of bad flux for a "poor" - if (!status) { - poorFrac = 0.0; - psWarning("POOR.FRAC is not set in the %s recipe --- defaulting to %f.", PSWARP_RECIPE, poorFrac); - } - - bool PSF = psMetadataLookupBool(&status, recipe, "PSF"); ///< Generate a PSF model? - if (!status) { - PSF = true; - psWarning("PSF is not set in the %s recipe --- defaulting to TRUE.", PSWARP_RECIPE); - } - - 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, - "Iso-astrom grid spacing in x", nGridX); - psMetadataAddS32(recipe, PS_LIST_TAIL, "GRID.NY", PS_META_REPLACE, - "Iso-astrom grid spacing in y", nGridY); - psMetadataAddStr(recipe, PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, - "Interpolation mode", name); - psMetadataAddS32(recipe, PS_LIST_TAIL, "INTERPOLATION.NUM", PS_META_REPLACE, - "Interpolation pre-calculated kernels", numKernels); - psMetadataAddF32(recipe, PS_LIST_TAIL, "POOR.FRAC", PS_META_REPLACE, - "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, - "Iso-astrom grid spacing in x", nGridX); - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "GRID.NY", 0, - "Iso-astrom grid spacing in y", nGridY); - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.MODE", 0, - "Interpolation mode", interpolationMode); - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.NUM", 0, - "Interpolation pre-calculated kernels", numKernels); - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "POOR.FRAC", 0, - "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"); - - return (config); -} Index: /trunk/pswarp/src/pswarpCleanup.c =================================================================== --- /trunk/pswarp/src/pswarpCleanup.c (revision 35563) +++ /trunk/pswarp/src/pswarpCleanup.c (revision 35563) @@ -0,0 +1,120 @@ +/** @file pswarp.c + * + * @brief cleanup the open files & memory then exit + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-11 23:27:58 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pswarp.h" + +// Lists of file rules which we may need to save here +static char *outputFiles[] = { + "PSWARP.OUTPUT", + "PSWARP.OUTPUT.MASK", + "PSWARP.OUTPUT.VARIANCE", + "PSWARP.OUTPUT.BKGMODEL", + "PSWARP.OUTPUT.SOURCES", + "PSPHOT.OUTPUT", + "PSPHOT.RESID", + "PSPHOT.BACKMDL", + "PSPHOT.BACKMDL.STDEV", + "PSPHOT.BACKGND", + "PSPHOT.BACKSUB", + "PSPHOT.PSF.SAVE", + "SOURCE.PLOT.MOMENTS", + "SOURCE.PLOT.PSFMODEL", + "SOURCE.PLOT.APRESID", + NULL +}; + +void pswarpCleanup (pmConfig *config, pswarpStatsFile *statsFile) +{ + psExit exitValue = pswarpExitCode(PS_EXIT_SUCCESS); // Exit code + + // activate all of the relevant output files + pmFPAfileActivate(config->files, false, NULL); + for (int i = 0; outputFiles[i] != NULL; i++) { + pmFPAfileActivate(config->files, true, outputFiles[i]); + } + + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT"); + if (!output) { + psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n"); + pmFPAfileFreeSetStrict(false); + goto DONE; + } + + pmFPAview *view = pmFPAviewAlloc(0); + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + pmFPAfileFreeSetStrict(false); + goto DONE; + } + pmChip *chip; + while ((chip = pmFPAviewNextChip (view, output->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."); + pmFPAfileFreeSetStrict(false); + goto DONE; + } + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, output->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."); + pmFPAfileFreeSetStrict(false); + goto DONE; + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + pmFPAfileFreeSetStrict(false); + goto DONE; + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + pmFPAfileFreeSetStrict(false); + goto DONE; + } + psFree(view); + + if (!pswarpStatsFileSave (config, statsFile)) { + psError(psErrorCodeLast(), false, "Unable to save stats file."); + exitValue = pswarpExitCode(exitValue); + } + + // Dump configuration + bool mdok = false; + psString dump_file = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); + if (dump_file) { + if (!pmConfigDump(config, dump_file)) { + psError(psErrorCodeLast(), false, "Unable to dump configuration"); + exitValue = pswarpExitCode(exitValue); + } + } + + psThreadPoolFinalize(); + psMemCheckCorruption(stderr, true); + + psFree(config); + + psTimerStop(); + pmVisualClose(); + pmModelClassCleanup(); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + +DONE: + exitValue = pswarpExitCode(exitValue); + exit (exitValue); +} Index: /trunk/pswarp/src/pswarpDefine.c =================================================================== --- /trunk/pswarp/src/pswarpDefine.c (revision 35562) +++ /trunk/pswarp/src/pswarpDefine.c (revision 35563) @@ -1,6 +1,5 @@ /** @file pswarpDefine.c * - * @brief - * + * @brief generates the output images * @ingroup pswarp * @@ -13,7 +12,4 @@ # include "pswarp.h" -/** - * loads the skycell metadata - */ bool pswarpDefine (pmConfig *config) { @@ -42,53 +38,90 @@ } - // 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 the output astrometry + // XXX this section loops over the input skycell to load the header info + { + pmFPAview *view = pmFPAviewAlloc(0); - // 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; - } + pmChip *chip; + if (!pmFPAfileRead (skycell, view, config)) { + psError(PS_ERR_IO, false, "failed READ at FPA %s", skycell->name); + psFree(view); + return false; + } + while ((chip = pmFPAviewNextChip (view, skycell->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 (!pmFPAfileRead (skycell, view, config)) { + psError(psErrorCodeLast(), false, "failed READ at CHIP %s", skycell->name); + psFree(view); + return false; + } + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, skycell->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 (!pmFPAfileRead (skycell, view, config)) { + psError(psErrorCodeLast(), false, "failed READ at CELL %s", skycell->name); + psFree(view); + return false; + } + if (!pmFPAfileClose(skycell, view)) { + psError(psErrorCodeLast(), false, "failed CLOSE at CELL %s", skycell->name); + psFree (view); + return false; + } - pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell - pmReadout *readout = pmReadoutAlloc(target); ///< Target readout - readout->image = psImageAlloc(numCols, numRows, 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"); - psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE"); - psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE"); - 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"); - - + // we've got the output astrom header + pmHDU *hdu = pmHDUFromCell(cell); ///< 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; + } + + // XXX generate the output pixels + pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell + pmReadout *readout = pmReadoutAlloc(target); ///< Target readout + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->image, NAN); + psFree(readout); // Drop reference + + bool status = false; + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XBIN"); + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YBIN"); + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XSIZE"); + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YSIZE"); + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XPARITY"); + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YPARITY"); + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.X0"); + psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.Y0"); + } + if (!pmFPAfileClose(skycell, view)) { + psError(psErrorCodeLast(), false, "failed CLOSE at CHIP %s", skycell->name); + psFree (view); + return false; + } + } + if (!pmFPAfileClose(skycell, view)) { + psError(psErrorCodeLast(), false, "failed CLOSE at FPA %s", skycell->name); + psFree (view); + return false; + } + psFree(view); } - // XXX this is not a sufficient test + // the section below converts the header astrometry info (in skycell) to the pmFPA + // astrometry structures, saving them on the output fpa structure. + + // XXX this block should be modified to loop over all chips + + pmFPAview *view = pmFPAviewAlloc (0); view->chip = 0; view->cell = 0; @@ -134,144 +167,2 @@ } -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/pswarpDefineBackground.c =================================================================== --- /trunk/pswarp/src/pswarpDefineBackground.c (revision 35563) +++ /trunk/pswarp/src/pswarpDefineBackground.c (revision 35563) @@ -0,0 +1,155 @@ +/** @file pswarpDefineBackground.c + * + * @brief generates the output background model + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pswarp.h" + +bool pswarpDefineBackground (pmConfig *config) { + + if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) return true; + + // 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/pswarpDefineLayout.c =================================================================== --- /trunk/pswarp/src/pswarpDefineLayout.c (revision 35563) +++ /trunk/pswarp/src/pswarpDefineLayout.c (revision 35563) @@ -0,0 +1,196 @@ +/** @file pswarpDefineLayout.c + * + * @brief load input & output astrometry, determine overlaps + * @ingroup pswarp + * + * @author IfA + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pswarp.h" + +bool pswarpGenerateOutputCell (pmChip *chip, pmCell *cell, pmCell *refcell, float xBin, float yBin); + +// load the astrometry header info and generate the tranformations +bool pswarpDefineLayout (pmConfig *config) { + + bool status = false; + + pmFPAfile *output = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT"); + if (!output) { + psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n"); + return false; + } + // select the input data sources + pmFPAfile *skycell = psMetadataLookupPtr (&status, config->files, "PSWARP.SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_CONFIG, false, "Can't find skycell data!\n"); + return false; + } + if (!pswarpLoadAstrometry (output, skycell, config)) { + psError(PSWARP_ERR_CONFIG, false, "problem loading output astrometry\n"); + return false; + } + + // if a background model is supplied and the output is requested, then + // this will exist. if not, skip related features + pmFPAfile *bkgModel = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT.BKGMODEL"); + + // chips are not processed unless we have determined they overlap the inputs + pmFPAExcludeChips (output->fpa); + + 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; + } + + for (int i = 0; i < nInputs; i++) { + // place input astrometry transformations in 'input' + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSWARP.INPUT", i); + if (!input) { + psError(PSWARP_ERR_CONFIG, false, "Can't find input data!\n"); + return false; + } + + // input astrometry may be embedded in 'input' or supplied separately + pmFPAfile *astrom = pmFPAfileSelectSingle(config->files, "PSWARP.ASTROM", i); + if (!astrom) { + astrom = input; + } + if (astrom->camera != input->camera) { + // XXX this may cause an error for output / skycell + psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match."); + return false; + } + if (!pswarpLoadAstrometry (input, astrom, config)) { + psError(PSWARP_ERR_CONFIG, false, "problem loading input astrometry\n"); + return false; + } + + // given the input data, determine which output elements should be generated + // and generate the appropriate images + + // find the R,D center for the skycell, use for common projection + psProjection *frame = pswarpLocalFrame (skycell->fpa); + + // generate Lmin,max, Mmin,max for both datasets + pswarpBounds *srcBounds = pswarpMakeBounds (input->fpa, frame); + pswarpBounds *tgtBounds = pswarpMakeBounds (skycell->fpa, frame); + + // find the output (tgt) chips which overlap the input (src) chips + pswarpFindOverlap (input->fpa, output->fpa, srcBounds, tgtBounds); + psFree (srcBounds); + psFree (tgtBounds); + } + + // The loop below generates the output pixels. XXX Should this be deferred until we + // actually need them? + + // Generate the output chips (pixels on output->fpa, concepts from skycell->fpa) + pmFPAview *view = pmFPAviewAlloc(0); + + // check input astrometry mode + bool bilevelAstrometry = psMetadataLookupBool (NULL, skycell->fpa->analysis, "ASTROMETRY.BILEVEL"); + + // XXX this is a bit muddled : we are setting up the astrometry elements for the output + // background model here instead of in pswarpLoadAstrometry, where it would be more natural + // the only explanation I have is that 1) I need to have the binning factors for the input image, but + // 2) at this point I don't yet have the input background models loaded and 3) the source + // for the binning factor of the input and output background models is different + // (header->CCDSUMi for input, recipe / pmFPAfile.xbin.ybin for the other) + + // if we have an output background model, generate hdu, etc and supply astrometry + if (bkgModel && bilevelAstrometry) { + // top-level elements of the bkgModel astrometry + bkgModel->fpa->toTPA = psMemIncrRefCounter (skycell->fpa->toTPA); + bkgModel->fpa->fromTPA = psMemIncrRefCounter (skycell->fpa->fromTPA); + bkgModel->fpa->toSky = psMemIncrRefCounter (skycell->fpa->toSky); + } + + pmChip *chip; + while ((chip = pmFPAviewNextChip (view, output->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process) { continue; } + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, output->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process) { continue; } + + pmCell *refcell = pmFPAviewThisCell(view, skycell->fpa); ///< Target cell + + if (!pswarpGenerateOutputCell (chip, cell, refcell, output->xBin, output->yBin)){ + psError(PSWARP_ERR_DATA, false, "failed to generate output cell"); + psFree(view); + return false; + } + + if (bkgModel) { + pmChip *bkgChip = pmFPAviewThisChip(view, bkgModel->fpa); ///< Target cell + pmCell *bkgCell = pmFPAviewThisCell(view, bkgModel->fpa); ///< Target cell + if (!pswarpGenerateOutputCell (bkgChip, bkgCell, refcell, bkgModel->xBin, bkgModel->yBin)) { + psError(PSWARP_ERR_DATA, false, "failed to generate output cell"); + psFree(view); + return false; + } + if (!pswarpModifyChipAstrom (config, view, bkgChip, skycell, bilevelAstrometry, bkgModel->xBin, bkgModel->yBin)) { + psError(PSWARP_ERR_DATA, false, "problem with output astrometry"); + psFree(view); + return false; + } + } + } + } + psFree (view); + return true; +} + +bool pswarpGenerateOutputCell (pmChip *chip, pmCell *cell, pmCell *refcell, float xBin, float yBin) { + + bool status = false; + + // we've got the output astrom header + pmHDU *hdu = pmHDUFromCell(refcell); ///< HDU for source + if (!hdu || !hdu->header) { + psError(PM_ERR_PROG, false, "Unable to find header for sky cell."); + return false; + } + int numCols = psMetadataLookupS32(&status, hdu->header, "NAXIS1") / xBin; ///< Number of output columns + int numRows = psMetadataLookupS32(&status, hdu->header, "NAXIS2") / yBin; ///< Number of output rows + if ((numCols == 0) || (numRows == 0)) { + psError(PSWARP_ERR_DATA, false, "output cell has invalid dimensions %d x %d", numCols, numRows); + return false; + } + + // generate the pixels for output->fpa + pmReadout *readout = pmReadoutAlloc(cell); ///< output readout + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + + psImageInit(readout->image, NAN); + psFree(readout); // Drop reference + + // copy the image concepts from the skycell + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XBIN"); + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YBIN"); + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XSIZE"); + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YSIZE"); + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XPARITY"); + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YPARITY"); + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.X0"); + psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.Y0"); + + chip->data_exists = true; + chip->file_exists = true; + chip->process = true; + + cell->data_exists = true; + cell->file_exists = true; + cell->process = true; + + // copy the basic headers across from astrom ref to output + pmHDU *outHDU = pmHDUFromCell (cell); ///< HDU for the output warped image + outHDU->header = psMetadataCopy(outHDU->header, hdu->header); + pswarpVersionHeader(outHDU->header); + return true; +} Index: /trunk/pswarp/src/pswarpDefineSkycell.c =================================================================== --- /trunk/pswarp/src/pswarpDefineSkycell.c (revision 35562) +++ /trunk/pswarp/src/pswarpDefineSkycell.c (revision 35563) @@ -73,15 +73,18 @@ psFitsClose(fits); - // We need to force the format for the skycell to be equivalent to SIMPLE. Determine - // the current format from the header; Determine camera if not specified already - // XXX EAM : this operation should be defined as a pmConfig function (pmConfigCopy?) - skyConfig = pmConfigAlloc(); - skyConfig->user = psMemIncrRefCounter(config->user); - skyConfig->system = psMemIncrRefCounter(config->system); + { // this bit is unique to pswarpDefineSkycell: - psFree (skyConfig->files); - skyConfig->files = psMemIncrRefCounter(config->files); - psFree (skyConfig->arguments); - skyConfig->arguments = psMemIncrRefCounter(config->arguments); + // We need to force the format for the skycell to be equivalent to SIMPLE. Determine + // the current format from the header; Determine camera if not specified already + // XXX EAM : this operation should be defined as a pmConfig function (pmConfigCopy?) + skyConfig = pmConfigAlloc(); + skyConfig->user = psMemIncrRefCounter(config->user); + skyConfig->system = psMemIncrRefCounter(config->system); + + psFree (skyConfig->files); + skyConfig->files = psMemIncrRefCounter(config->files); + psFree (skyConfig->arguments); + skyConfig->arguments = psMemIncrRefCounter(config->arguments); + } format = pmConfigCameraFormatFromHeader (NULL, NULL, NULL, skyConfig, phu, false); @@ -94,6 +97,6 @@ // Record the camera name of the skycell, so we can save its configuration - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SKYCELL.CAMERA", 0, "Name of camera for skycell", - skyConfig->cameraName); + // XXX this can be put outside of the f + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SKYCELL.CAMERA", 0, "Name of camera for skycell", skyConfig->cameraName); // build the template fpa, set up the basic view Index: unk/pswarp/src/pswarpExit.c =================================================================== --- /trunk/pswarp/src/pswarpExit.c (revision 35562) +++ (revision ) @@ -1,82 +1,0 @@ -/** @file pswarpExit.c - * - * @brief - * - * @ingroup pswarp - * - * @author IfA - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-05 20:44:04 $ - * Copyright 2009 Institute for Astronomy, University of Hawaii - */ - -# include "pswarp.h" - - -psExit pswarpExitCode(psExit exitValue) -{ - if (exitValue != PS_EXIT_SUCCESS) { - return exitValue; - } - - // gcc -Wswitch complains here if err is declared as type psErrorCode - // the collection of ps*ErrorCode values are enums defined separately for - // each module (psphot, pswarp, etc). the lowest type, psErrorCode is only the base set and does - // not include the possible psphot values - - // for now, to get around this, we just use an int for the switch - - // psErrorCode errorCode = psErrorCodeLast(); // Error code - int errorCode = psErrorCodeLast(); // Error code - - if (errorCode != PS_ERR_NONE) { - pmFPAfileFreeSetStrict(false); - psErrorStackPrint(stderr, "Unable to perform warp."); - switch (errorCode) { - case PSWARP_ERR_UNKNOWN: - case PS_ERR_UNKNOWN: - psLogMsg("pswarp", PS_LOG_WARN, "Unknown error code: %x", errorCode); - return PS_EXIT_UNKNOWN_ERROR; - case PM_ERR_CONFIG: - // We get this when the config file is not readable. Set to data like error - psLogMsg("pswarp", PS_LOG_WARN, "Config Error: %x", errorCode); - exitValue = PS_EXIT_SYS_ERROR; - break; - case PS_ERR_IO: - case PS_ERR_DB_CLIENT: - case PS_ERR_DB_SERVER: - case PS_ERR_BAD_FITS: - case PS_ERR_OS_CALL_FAILED: - case PM_ERR_SYS: - case PSWARP_ERR_IO: - psLogMsg("pswarp", PS_LOG_WARN, "I/O error code: %x", errorCode); - return PS_EXIT_SYS_ERROR; - case PS_ERR_BAD_PARAMETER_VALUE: - case PS_ERR_BAD_PARAMETER_TYPE: - case PS_ERR_BAD_PARAMETER_NULL: - case PS_ERR_BAD_PARAMETER_SIZE: - case PSWARP_ERR_ARGUMENTS: - case PSWARP_ERR_CONFIG: - psLogMsg("pswarp", PS_LOG_WARN, "Configuration error code: %x", errorCode); - return PS_EXIT_CONFIG_ERROR; - case PSPHOT_ERR_PSF: - case PSWARP_ERR_DATA: - case PSWARP_ERR_NO_OVERLAP: - psLogMsg("pswarp", PS_LOG_WARN, "Data error code: %x", errorCode); - return PS_EXIT_DATA_ERROR; - case PS_ERR_UNEXPECTED_NULL: - case PS_ERR_PROGRAMMING: - case PSWARP_ERR_NOT_IMPLEMENTED: - psLogMsg("pswarp", PS_LOG_WARN, "Programming error code: %x", errorCode); - exitValue = PS_EXIT_PROG_ERROR; - break; - default: - // It's a programming error if we're not dealing with the error correctly - psLogMsg("pswarp", PS_LOG_WARN, "Unrecognised error code: %x", errorCode); - return PS_EXIT_PROG_ERROR; - } - } - - return exitValue; -} - Index: /trunk/pswarp/src/pswarpExitCode.c =================================================================== --- /trunk/pswarp/src/pswarpExitCode.c (revision 35563) +++ /trunk/pswarp/src/pswarpExitCode.c (revision 35563) @@ -0,0 +1,76 @@ +/** @file pswarpExit.c + * + * @brief + * + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pswarp.h" + +psExit pswarpExitCode(psExit exitValue) +{ + if (exitValue != PS_EXIT_SUCCESS) { + return exitValue; + } + + // gcc -Wswitch complains here if err is declared as type psErrorCode + // the collection of ps*ErrorCode values are enums defined separately for + // each module (psphot, pswarp, etc). the lowest type, psErrorCode is only the base set and does + // not include the possible psphot values + + // for now, to get around this, we just use an int for the switch + + // psErrorCode errorCode = psErrorCodeLast(); // Error code + int errorCode = psErrorCodeLast(); // Error code + + if (errorCode != PS_ERR_NONE) { + pmFPAfileFreeSetStrict(false); + psErrorStackPrint(stderr, "Unable to perform warp."); + switch (errorCode) { + case PSWARP_ERR_UNKNOWN: + case PS_ERR_UNKNOWN: + psLogMsg("pswarp", PS_LOG_WARN, "Unknown error code: %x", errorCode); + return PS_EXIT_UNKNOWN_ERROR; + case PS_ERR_IO: + case PS_ERR_DB_CLIENT: + case PS_ERR_DB_SERVER: + case PS_ERR_BAD_FITS: + case PS_ERR_OS_CALL_FAILED: + case PM_ERR_SYS: + case PSWARP_ERR_IO: + psLogMsg("pswarp", PS_LOG_WARN, "I/O error code: %x", errorCode); + return PS_EXIT_SYS_ERROR; + case PS_ERR_BAD_PARAMETER_VALUE: + case PS_ERR_BAD_PARAMETER_TYPE: + case PS_ERR_BAD_PARAMETER_NULL: + case PS_ERR_BAD_PARAMETER_SIZE: + case PSWARP_ERR_ARGUMENTS: + case PSWARP_ERR_CONFIG: + psLogMsg("pswarp", PS_LOG_WARN, "Configuration error code: %x", errorCode); + return PS_EXIT_CONFIG_ERROR; + case PSPHOT_ERR_PSF: + case PSWARP_ERR_DATA: + case PSWARP_ERR_NO_OVERLAP: + psLogMsg("pswarp", PS_LOG_WARN, "Data error code: %x", errorCode); + return PS_EXIT_DATA_ERROR; + case PS_ERR_UNEXPECTED_NULL: + case PS_ERR_PROGRAMMING: + case PSWARP_ERR_NOT_IMPLEMENTED: + psLogMsg("pswarp", PS_LOG_WARN, "Programming error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; + default: + // It's a programming error if we're not dealing with the error correctly + psLogMsg("pswarp", PS_LOG_WARN, "Unrecognised error code: %x", errorCode); + return PS_EXIT_PROG_ERROR; + } + } + + return exitValue; +} + Index: /trunk/pswarp/src/pswarpLoadAstrometry.c =================================================================== --- /trunk/pswarp/src/pswarpLoadAstrometry.c (revision 35563) +++ /trunk/pswarp/src/pswarpLoadAstrometry.c (revision 35563) @@ -0,0 +1,152 @@ +/** @file pswarpLoadAstrometry.c + * + * @brief load the headers and construct astrometric transformations + * @ingroup pswarp + * + * @author IfA + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pswarp.h" + +// we have a pmFPAfile which defines the astrometry via headers. we force it to have type +// WCS and read the headers. We place the resulting info on the target astrometry +// containers +bool pswarpLoadAstrometry (pmFPAfile *target, pmFPAfile *astrom, pmConfig *config) { + + pmChip *chip = NULL; + pmCell *cell = NULL; + + // XXX assert that astrom be WCS or CMF? + + // XXX set the type to be WCS + pmFPAfileType saveType = astrom->type; + astrom->type = PM_FPA_FILE_WCS; + + // **** read in all of the headers from the astrometry source + + pmFPAview *view = pmFPAviewAlloc(0); + + if (!pmFPAfileRead (astrom, view, config)) { + psError(PS_ERR_IO, false, "failed READ at FPA %s", astrom->name); + psFree(view); + return false; + } + + while ((chip = pmFPAviewNextChip (view, astrom->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 (!pmFPAfileRead (astrom, view, config)) { + psError(psErrorCodeLast(), false, "failed READ at CHIP %s", astrom->name); + psFree(view); + return false; + } + while ((cell = pmFPAviewNextCell (view, astrom->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process) { continue; } + if (!pmFPAfileRead (astrom, view, config)) { + psError(psErrorCodeLast(), false, "failed READ at CELL %s", astrom->name); + psFree(view); + return false; + } + if (!pmFPAfileClose(astrom, view)) { + psError(psErrorCodeLast(), false, "failed CLOSE at CELL %s", astrom->name); + psFree (view); + return false; + } + + // we've got the output astrom header + pmHDU *hdu = pmHDUFromCell(cell); ///< HDU for source + if (!hdu || !hdu->header) { + psError(PM_ERR_PROG, false, "Unable to find header for sky cell."); + psFree(view); + return false; + } + } + if (!pmFPAfileClose(astrom, view)) { + psError(psErrorCodeLast(), false, "failed CLOSE at CHIP %s", astrom->name); + psFree (view); + return false; + } + } + if (!pmFPAfileClose(astrom, view)) { + psError(psErrorCodeLast(), false, "failed CLOSE at FPA %s", astrom->name); + psFree (view); + return false; + } + + // **** convert the header astrometry info (in astrom) to the pmFPA astrometry + // structures, saving them on the target fpa structure. + + pmFPAviewReset (view); + + bool bilevelAstrometry = false; + pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa); ///< Astrometry PHU + if (phu) { + char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp(&ctype[4], "-DIS"); + } + } + // apply the bilevel astrometry elements to the target + if (bilevelAstrometry) { + if (!pmAstromReadBilevelMosaic(astrom->fpa, phu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell."); + psFree(view); + return false; + } + } + psMetadataAddBool (astrom->fpa->analysis, PS_LIST_TAIL, "ASTROMETRY.BILEVEL", PS_META_REPLACE, "bilevel astrometry?", bilevelAstrometry); + + // for pswarpLocalFrame, I need transformations and HDUs on a single fpa (so set on astrom as well as target) + // target->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA); + // target->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA); + // target->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky); + + while ((chip = pmFPAviewNextChip (view, astrom->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; } + + // we've got the astrom header + pmHDU *hdu = pmHDUFromChip(chip); ///< HDU for source + if (!hdu || !hdu->header) { + psError(PM_ERR_PROG, false, "Unable to find header for sky cell."); + psFree(view); + return false; + } + + // We read from the astrometry source into the target. + pmChip *targetChip = pmFPAviewThisChip(view, target->fpa); ///< Chip in the output + if (!targetChip) continue; // only load astrometry into output chips which exist! + + if (bilevelAstrometry) { + if (!pmAstromReadBilevelChip(chip, 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(astrom->fpa, chip, hdu->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell."); + psFree(view); + return false; + } + } + targetChip->toFPA = psMemIncrRefCounter (chip->toFPA); + targetChip->fromFPA = psMemIncrRefCounter (chip->fromFPA); + } + + // for pswarpLocalFrame, I need transformations and HDUs on a single fpa (so set on astrom as well as target) + target->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA); + target->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA); + target->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky); + + // reset the type to the original value + astrom->type = saveType; + + psFree(view); + return true; +} + Index: /trunk/pswarp/src/pswarpLoop.c =================================================================== --- /trunk/pswarp/src/pswarpLoop.c (revision 35562) +++ /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; -} Index: /trunk/pswarp/src/pswarpLoopBackground.c =================================================================== --- /trunk/pswarp/src/pswarpLoopBackground.c (revision 35563) +++ /trunk/pswarp/src/pswarpLoopBackground.c (revision 35563) @@ -0,0 +1,276 @@ +/** @file pswarpLoopBackground.c + * + * @brief background model processing loop for pswarp + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pswarp.h" + +#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. +bool pswarpLoopBackground(pmConfig *config, psMetadata *stats) +{ + bool status; + bool mdok; // Status of MD lookup + + if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) return true; + + int xGrid = psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID"); + int yGrid = psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID"); + + // WHAT IS THIS?? should it be disabled at the end? + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8); + + // 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; + } + + // select the input data sources + pmFPAfile *output = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT.BKGMODEL"); + if (!output) { + psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n"); + return false; + } + + // use the external astrometry source if supplied + pmFPAfile *skycell = psMetadataLookupPtr(&status, config->files, "PSWARP.SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry."); + return false; + } + + pmFPAview *view = pmFPAviewAlloc(0); + + 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; + } + + // we are only reading the background models, not writing anything + pmFPAfileActivate(config->files, false, NULL); + pmFPAfileActivate(config->files, true, "PSWARP.BKGMODEL"); + + // loop over this section once per input group + for (int i = 0; i < nInputs; i++) { + + // select the input data source : we are reading & transforming pixels for the background model + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSWARP.BKGMODEL", 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; + } + + // check input astrometry mode + bool bilevelAstrometry = psMetadataLookupBool (NULL, astrom->fpa->analysis, "ASTROMETRY.BILEVEL"); + if (bilevelAstrometry) { + // top-level elements of the input model astrometry + input->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA); + input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA); + input->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky); + } + + pmFPAviewReset (view); + + pmFPAview *view = pmFPAviewAlloc(0); + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto FAIL; + } + + pmChip *chip; + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + if (!chip->process || !chip->file_exists) { continue; } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto FAIL; + } + + // adjust output astrometry + double xBin = NAN; + double yBin = NAN; + pswarpGetInputScales (&xBin, &yBin, config, view, chip); + double fluxRatio = xBin * yBin / (float) (xGrid * yGrid); + + if (!pswarpModifyChipAstrom (config, view, chip, astrom, bilevelAstrometry, xBin, yBin)) { + psError(psErrorCodeLast(), false, "failed to set BKGMODEL astrometry."); + goto FAIL; + } + + 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)) { + 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; + } + + // re-normalize the BKGMODEL pixels by modified astrometry + for (int x = 0; x < readout->image->numCols; x++) { + for (int y = 0; y < readout->image->numRows; y++) { + readout->image->data.F32[y][x] *= fluxRatio; + } + } + + // transform the actual BKGMODEL pixels from readout to output + pswarpTransformToTarget (output->fpa, readout, config, true); + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Readout + psError(psErrorCodeLast(), false, "Unable to close files."); + goto FAIL; + } + } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Cell + psError(psErrorCodeLast(), false, "Unable to close files."); + goto FAIL; + } + } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Chip + psError(psErrorCodeLast(), false, "Unable to close files."); + goto FAIL; + } + } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // FPA + psError(psErrorCodeLast(), false, "Unable to close files."); + goto FAIL; + } + + if (!pswarpUpdateMetadata (output->fpa, skycell->fpa, input->fpa, astrom->fpa, config, false)) { + psError(psErrorCodeLast(), false, "problem generating statistics."); + goto FAIL; + } + } + + psFree (view); + return true; + +FAIL: + psFree (view); + return false; +} + +// **** input == (small) BKGMODEL, astrom == (fullsize) ASTROM +// use astrometry from input ref astrometry and BINNING info from input model to modify the astrometry +// elements of input chip. scale is the area of an input pixel is raw image pixels +bool pswarpGetInputScales (double *xBin, double *yBin, pmConfig *config, pmFPAview *view, pmChip *chip) { + + // Pull information from the header of the background files so we can use it to set values. + pmHDU *hdu = pmFPAviewThisHDU(view, chip->parent); + 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")); + + *xBin = CCDSUM1; + *yBin = CCDSUM2; + + return true; +} + +// **** input == (small) BKGMODEL, astrom == (fullsize) ASTROM +// use astrometry from input ref astrometry and BINNING info from input model to modify the astrometry +// elements of input chip. scale is the area of an input pixel is raw image pixels +bool pswarpModifyChipAstrom (pmConfig *config, pmFPAview *view, pmChip *chip, pmFPAfile *astrom, bool bilevelAstrometry, double xBin, double yBin) { + + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); + + // generate a WCS structure from the ASTROM header keywords + pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); + + // re-scale the terms of the WCS to match the BKGMODEL scale + double cd1f = xBin; + double cd2f = yBin; + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; + + // WCS->trans->x->nX/nY + // adjust the scale of the WCS coeffs + 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); + } + } + + // this is a bit crude : use a locally temporary header to store the WCS information + // which is then applied to the input chip (why not modify the input chip header?) + psMetadata *tmpHeader = psMetadataAlloc(); + + // write the modified WCS to the header + // XXX this should probably not modify the astrom header + pmAstromWCStoHeader (tmpHeader, WCS); + + // use the modified header to update the WCS elements of the BKGMODEL fpa/chip structures + if (bilevelAstrometry) { + if (!pmAstromReadBilevelChip (chip, tmpHeader)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA."); + return false; + } + } else { + // we use a default FPA pixel scale of 1.0 + if (!pmAstromReadWCS (chip->parent, chip, tmpHeader, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); + return false; + } + } + + return true; +} Index: /trunk/pswarp/src/pswarpLoopSkycell.c =================================================================== --- /trunk/pswarp/src/pswarpLoopSkycell.c (revision 35563) +++ /trunk/pswarp/src/pswarpLoopSkycell.c (revision 35563) @@ -0,0 +1,447 @@ +/** @file pswarpLoop.c + * + * @brief mail processing loop for pswarp + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#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. +bool pswarpLoopSkycell(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; + } + + // 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"); + 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"); + if (!output) { + psError(PSWARP_ERR_CONFIG, true, "Can't find output 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. + // XXX why is this being done here, doesn't it duplicate the work in pswarpDefine.c?? + pswarpFileActivation(config, detectorFiles, false); + pswarpFileActivation(config, photFiles, false); + pswarpFileActivation(config, independentFiles, false); + pswarpFileActivation(config, skycellFiles, false); + + // pswarpFileActivation(config, skycellFiles, true); + // if (0) { + // if (!pswarpIOChecksBefore(config)) { + // psError(psErrorCodeLast(), false, "Unable to read files."); + // goto DONE; + // } + // } + + // 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"); + } + + // 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); + 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; + } + + // *** main transformation block + 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 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); + 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); + 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)) { + 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; + } + } + + 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(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(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."); + 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: /trunk/pswarp/src/pswarpMakePSF.c =================================================================== --- /trunk/pswarp/src/pswarpMakePSF.c (revision 35563) +++ /trunk/pswarp/src/pswarpMakePSF.c (revision 35563) @@ -0,0 +1,112 @@ +/** @file pswarpLoop.c + * + * @brief mail processing loop for pswarp + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pswarp.h" + +// 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. + +bool pswarpMakePSF (pmConfig *config, pmFPAfile *output, psMetadata *stats) { + + bool status = false; + + // 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 (!psMetadataLookupBool(&status, recipe, "PSF")) { + psLogMsg("pswarp", PS_LOG_INFO, "Skipping PSF measurement"); + return true; + } + + // XXX move above the loop + pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT); + + // supply the readout and fpa of interest to psphot + pmFPAfile *photFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.INPUT"); + pmFPACopy(photFile->fpa, output->fpa); + + pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell + + pmChip *chip; + while ((chip = pmFPAviewNextChip (view, photFile->fpa, 1)) != NULL) { + psTrace ("pswarpMakePSF", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, photFile->fpa, 1)) != NULL) { + psTrace ("pswarpMakePSF", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, photFile->fpa, 1)) != NULL) { + if (!readout->data_exists) { + continue; + } + + // 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"); + psLogMsg ("psphot", 3, "no psf sources for this readout (%d %d %d)", view->chip, view->cell, view->readout); + continue; + } + + // measure the PSF using these sources + if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) { + // A failure is likely a data quality issue. raise a quality flag, but do not skip other readouts + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to determine PSF for readout (%d %d %d)", view->chip, view->cell, view->readout); + psWarning("Unable to determine PSF --- suspect bad data quality : readout (%d %d %d)", view->chip, view->cell, view->readout); + if (stats) { + int currentQuality = psMetadataLookupS32(&status, stats, "QUALITY"); + if (currentQuality == 0) { + psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, "Unable to determine PSF", psErrorCodeLast()); + } + } + psErrorClear(); + } + } + } + + // Ensure seeing is carried over from the PSPHOT.INPUT chip to the output chip + pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Chip with seeing + psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing + item->data.F32 = psMetadataLookupF32(&status, chip->concepts, "CHIP.SEEING"); + } + + psFree(view); + return true; +} + +# if (0) +bool pswarpMakePSF_test () { + +#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); + return true; +} +# endif Index: /trunk/pswarp/src/pswarpMaskStats.c =================================================================== --- /trunk/pswarp/src/pswarpMaskStats.c (revision 35562) +++ /trunk/pswarp/src/pswarpMaskStats.c (revision 35563) @@ -1,4 +1,3 @@ #include "pswarp.h" -#include bool pswarpMaskStats(const pmReadout *readout, psMetadata *stats, const pmConfig *config) @@ -31,16 +30,13 @@ dynamicMaskVal,advisoryMaskVal)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate masks for readout."); - return(false); + return false; } - psMetadataAddS32(stats, PS_LIST_TAIL,"MASKFRAC_NPIX", 0, - "Number of valid pixels", Npix_valid); - psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_STATIC", 0, - "Fraction of pixels statically masked", (float) Npix_static / Npix_valid); - psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_DYNAMIC", 0, - "Fraction of pixels dynamically masked", (float) Npix_dynamic / Npix_valid); - psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAGIC", 0, - "Fraction of pixels magically masked", (float) Npix_magic / Npix_valid); - psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_ADVISORY", 0, - "Fraction of pixels masked as an advisory", (float) Npix_advisory / Npix_valid); - return(true); + + // XXX with multiple inputs (eg, output stacks -> exposure), these only represent the last input + psMetadataAddS32(stats, PS_LIST_TAIL,"MASKFRAC_NPIX", PS_META_REPLACE, "Number of valid pixels", Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_STATIC", PS_META_REPLACE, "Fraction of pixels statically masked", (float) Npix_static / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_DYNAMIC", PS_META_REPLACE, "Fraction of pixels dynamically masked", (float) Npix_dynamic / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAGIC", PS_META_REPLACE, "Fraction of pixels magically masked", (float) Npix_magic / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_ADVISORY", PS_META_REPLACE, "Fraction of pixels masked as an advisory", (float) Npix_advisory / Npix_valid); + return true; } Index: /trunk/pswarp/src/pswarpMatchRange.c =================================================================== --- /trunk/pswarp/src/pswarpMatchRange.c (revision 35562) +++ /trunk/pswarp/src/pswarpMatchRange.c (revision 35563) @@ -115,5 +115,5 @@ psPlaneTransformApply (srcFP, fpaSrc->fromTPA, srcTP); psPlaneTransformApply (srcPix, chipSrc->fromFPA, srcFP); - fprintf (stderr, "-> %f,%f ", srcPix->x, srcPix->y); + fprintf (stderr, "-> %f,%f\n", srcPix->x, srcPix->y); # endif Index: /trunk/pswarp/src/pswarpOptions.c =================================================================== --- /trunk/pswarp/src/pswarpOptions.c (revision 35563) +++ /trunk/pswarp/src/pswarpOptions.c (revision 35563) @@ -0,0 +1,120 @@ +/** @file pswarpOptionss.c + * + * @brief + * + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 03:10:36 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pswarp.h" + +/** + * Parse the recipe and format into the arguments + */ +bool pswarpOptions(pmConfig *config) +{ + // Select the appropriate recipe + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, true, "Can't find %s recipe!\n", PSWARP_RECIPE); + return false; + } + + // Get grid size + bool status; ///< Status of MD lookup + int nGridX = psMetadataLookupS32(&status, recipe, "GRID.NX"); + if (!status || nGridX <= 0) { + nGridX = 128; + psWarning("GRID.NX is not set in the recipe --- defaulting to %d", nGridX); + } + int nGridY = psMetadataLookupS32(&status, recipe, "GRID.NY"); + if (!status) { + nGridY = 128; + psWarning("GRID.NY is not set in the recipe --- defaulting to %d", nGridY); + } + + // Get interpolation mode + const char *name = psMetadataLookupStr (&status, recipe, "INTERPOLATION.MODE"); ///< Name of interp mode + if (!name) { + name = "BILINEAR"; + psLogMsg("pswarp", 3, "defaulting to %s interpolation", name); + } + psImageInterpolateMode interpolationMode = psImageInterpolateModeFromString(name); ///< Mode for interp. + if (interpolationMode == PS_INTERPOLATE_NONE) { + interpolationMode = PS_INTERPOLATE_BILINEAR; + psLogMsg ("pswarp", 3, + "Unknown interpolation mode %s, defaulting to bilinear interpolation\n", name); + name = "BILINEAR"; + } + + int numKernels = psMetadataLookupS32(&status, recipe, "INTERPOLATION.NUM"); + if (!status) { + numKernels = 0; + psWarning("INTERPOLATION.NUM is not set in the recipe --- defaulting to %d", numKernels); + } + + float poorFrac = psMetadataLookupF32(&status, recipe, "POOR.FRAC"); ///< Frac of bad flux for a "poor" + if (!status) { + poorFrac = 0.0; + psWarning("POOR.FRAC is not set in the %s recipe --- defaulting to %f.", PSWARP_RECIPE, poorFrac); + } + + bool PSF = psMetadataLookupBool(&status, recipe, "PSF"); ///< Generate a PSF model? + if (!status) { + PSF = true; + psWarning("PSF is not set in the %s recipe --- defaulting to TRUE.", PSWARP_RECIPE); + } + + // BACKGROUND.MODEL gets set in config->arguments (to false) if no input model is found + bool doBKG = psMetadataLookupBool(&status,config->arguments, "BACKGROUND.MODEL"); ///< Generate the warped background model? + if (!status) { + // if it is not in config->arguments, look in recipe + 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, "Iso-astrom grid spacing in x", nGridX); + psMetadataAddS32(recipe, PS_LIST_TAIL, "GRID.NY", PS_META_REPLACE, "Iso-astrom grid spacing in y", nGridY); + psMetadataAddStr(recipe, PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "Interpolation mode", name); + psMetadataAddS32(recipe, PS_LIST_TAIL, "INTERPOLATION.NUM", PS_META_REPLACE, "Interpolation pre-calculated kernels", numKernels); + psMetadataAddF32(recipe, PS_LIST_TAIL, "POOR.FRAC", PS_META_REPLACE, "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, "Iso-astrom grid spacing in x", nGridX); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "GRID.NY", 0, "Iso-astrom grid spacing in y", nGridY); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.MODE", 0, "Interpolation mode", interpolationMode); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.NUM", 0, "Interpolation pre-calculated kernels", numKernels); + psMetadataAddF32(config->arguments, PS_LIST_TAIL, "POOR.FRAC", 0, "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 pswarpOptions...\n"); + + return (config); +} Index: /trunk/pswarp/src/pswarpOverlaps.c =================================================================== --- /trunk/pswarp/src/pswarpOverlaps.c (revision 35563) +++ /trunk/pswarp/src/pswarpOverlaps.c (revision 35563) @@ -0,0 +1,238 @@ +# include "pswarp.h" + +static void pswarpBoundsFree (pswarpBounds *bounds) { + psFree (bounds->Pmin); + psFree (bounds->Pmax); + psFree (bounds->Qmin); + psFree (bounds->Qmax); +} + +pswarpBounds *pswarpBoundsAlloc() { + + pswarpBounds *bounds = psAlloc(sizeof(pswarpBounds)); + psMemSetDeallocator(bounds, (psFreeFunc)pswarpBoundsFree); + + bounds->Pmin = psVectorAllocEmpty(32, PS_DATA_F32); + bounds->Pmax = psVectorAllocEmpty(32, PS_DATA_F32); + bounds->Qmin = psVectorAllocEmpty(32, PS_DATA_F32); + bounds->Qmax = psVectorAllocEmpty(32, PS_DATA_F32); + return bounds; +} + +bool pswarpBoundsAppend(pswarpBounds *bounds, float Pmin, float Pmax, float Qmin, float Qmax) { + psVectorAppend (bounds->Pmin, Pmin); + psVectorAppend (bounds->Pmax, Pmax); + psVectorAppend (bounds->Qmin, Qmin); + psVectorAppend (bounds->Qmax, Qmax); + return true; +} + +pswarpBounds *pswarpMakeBounds (pmFPA *fpa, psProjection *frame) { + + pswarpBounds *bounds = pswarpBoundsAlloc (); + + psPlane *CH = psPlaneAlloc(); + psPlane *FP = psPlaneAlloc(); + psPlane *TP = psPlaneAlloc(); + psPlane *FR = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + + pmFPAview *view = pmFPAviewAlloc (0); + + pmChip *chip = NULL; + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { + // this ensures one entry per chip (regardless of existence) + pswarpBoundsAppend (bounds, NAN, NAN, NAN, NAN); + continue; + } + + // we've got the output astrom header + pmHDU *hdu = pmHDUFromChip(chip); ///< HDU for source + assert (hdu); + assert (hdu->header); + 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, "astrom has invalid dimensions %d x %d", numCols, numRows); + psFree(view); + return false; + } + + float Pmin = +1e9; + float Pmax = -1e9; + float Qmin = +1e9; + float Qmax = -1e9; + + // find P,Q at corners and center of edges (and center just for ease of code) + for (float ix = 0.0; ix <= 1.0; ix += 0.5) { + for (float iy = 0.0; iy <= 1.0; iy += 0.5) { + + CH->x = ix * numCols; + CH->y = iy * numRows; + psPlaneTransformApply(FP, chip->toFPA, CH); + psPlaneTransformApply(TP, fpa->toTPA, FP); + psDeproject (sky, TP, fpa->toSky); + psProject (FR, sky, frame); + + Pmin = PS_MIN(Pmin, FR->x); + Pmax = PS_MAX(Pmax, FR->x); + Qmin = PS_MIN(Qmin, FR->y); + Qmax = PS_MAX(Qmax, FR->y); + } + } + pswarpBoundsAppend (bounds, Pmin, Pmax, Qmin, Qmax); + } + psFree(view); + psFree(CH); + psFree(FP); + psFree(TP); + psFree(FR); + psFree(sky); + + return bounds; +} + +// find R,D for the center of the output pmFPA +psProjection *pswarpLocalFrame (pmFPA *fpa) { + + // find R,D for each corner of each chip and take the average + psVector *Rvec = psVectorAllocEmpty (256, PS_DATA_F32); + psVector *Dvec = psVectorAllocEmpty (256, PS_DATA_F32); + + psPlane *CH = psPlaneAlloc(); + psPlane *FP = psPlaneAlloc(); + psPlane *TP = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + + pmFPAview *view = pmFPAviewAlloc (0); + + // we are looping over the chips of the OUTPUT fpa. chip->file_exists is false for all chips. + // at this point, I cannot exclude any of these chips (I do not know which output chips will not be created). + // NOTE : I either have to use SKYCELL fpa and ensure SKYCELL->chip->toFPA exists or grab hdu from SKYCELL + + pmChip *chip = NULL; + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + // if (!chip->file_exists) { continue; } // this fpa refers to the output file, so no chip exist + assert (chip->toFPA); + chip->process = true; // start with process true, we will de-activate later (pswarpFindOverlaps) + + // we've got the output astrom header + pmHDU *hdu = pmHDUFromChip(chip); ///< HDU for source + assert (hdu); + assert (hdu->header); + 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, "astrom has invalid dimensions %d x %d", numCols, numRows); + psFree(view); + return false; + } + + // get R,D for all corners + for (float ix = 0.0; ix <= 1.0; ix += 1.0) { + for (float iy = 0.0; iy <= 1.0; iy += 1.0) { + CH->x = ix * numCols; + CH->y = iy * numRows; + psPlaneTransformApply(FP, chip->toFPA, CH); + psPlaneTransformApply(TP, fpa->toTPA, FP); + psDeproject (sky, TP, fpa->toSky); + psVectorAppend (Rvec, sky->r); + psVectorAppend (Dvec, sky->d); + } + } + } + psFree(view); + psFree(CH); + psFree(FP); + psFree(TP); + psFree(sky); + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + psVectorStats (stats, Rvec, NULL, NULL, 0); + double Rmid = stats->sampleMean; + + psStatsInit (stats); + psVectorStats (stats, Dvec, NULL, NULL, 0); + double Dmid = stats->sampleMean; + + psProjection *frame = psProjectionAlloc (Rmid, Dmid, PS_RAD_DEG/3600.0, PS_RAD_DEG/3600.0, PS_PROJ_TAN); + return frame; +} + +bool pswarpFindOverlap (pmFPA *input, pmFPA *output, pswarpBounds *src, pswarpBounds *tgt) { + + bool status = false; + + // we have the source and target bounds. loop over all output elements and check if any of the source + // elements overlap it. + + // also add the list if input (src) cells which land on each output (tgt) chip we do this + // by only identifying the input and output chips and adding all input chip cells to each + // output cell + + assert (output->chips->n == tgt->Pmin->n); + assert (input->chips->n == src->Pmin->n); + + for (int j = 0; j < output->chips->n; j++) { + pmChip *chip = output->chips->data[j]; + + // we have src bounds + float Pmin = tgt->Pmin->data.F32[j]; + float Pmax = tgt->Pmax->data.F32[j]; + float Qmin = tgt->Qmin->data.F32[j]; + float Qmax = tgt->Qmax->data.F32[j]; + + psArray *inputChips = psArrayAllocEmpty(8); + + bool hasOverlap = false; + for (int i = 0; i < src->Pmin->n; i++) { + // overlaps in P? + if (!isfinite(src->Pmin->data.F32[i])) continue; + if (!isfinite(src->Pmax->data.F32[i])) continue; + if (!isfinite(src->Qmin->data.F32[i])) continue; + if (!isfinite(src->Qmax->data.F32[i])) continue; + + if (Pmin > src->Pmax->data.F32[i]) continue; + if (Pmax < src->Pmin->data.F32[i]) continue; + + if (Qmin > src->Qmax->data.F32[i]) continue; + if (Qmax < src->Qmin->data.F32[i]) continue; + + hasOverlap = true; + psArrayAdd (inputChips, 1, input->chips->data[i]); + } + chip->process |= hasOverlap; + if (hasOverlap) { + pmChipSelectCells (chip); + // get the complete list of input cells for this set of input chips + // pmConceptsAverageCells needs a psList (not a psArray) + psArray *inputCells = psArrayAllocEmpty(8); + for (int nChip = 0; nChip < inputChips->n; nChip++) { + pmChip *inputChip = inputChips->data[nChip]; + for (int nCell = 0; nCell < inputChip->cells->n; nCell++) { + psArrayAdd (inputCells, 8, inputChip->cells->data[nCell]); + } + } + + // place the inputCells on each of the output cells + for (int nCell = 0; nCell < chip->cells->n; nCell++) { + pmCell *outputCell = chip->cells->data[nCell]; + psList *outList = psMetadataLookupPtr(&status, outputCell->analysis, "INPUT.CELLS"); + if (!outList) { + outList = psListAlloc(NULL); + psMetadataAddPtr(outputCell->analysis, PS_LIST_TAIL, "INPUT.CELLS", PS_DATA_LIST , "input cells touching this output cell", outList); + psFree (outList); + } + for (int nCell = 0; nCell < inputCells->n; nCell ++) { + psListAdd (outList, PS_LIST_TAIL, inputCells->data[nCell]); + } + } + psFree (inputCells); + } + psFree (inputChips); + } + return true; +} Index: /trunk/pswarp/src/pswarpParseCamera.c =================================================================== --- /trunk/pswarp/src/pswarpParseCamera.c (revision 35562) +++ /trunk/pswarp/src/pswarpParseCamera.c (revision 35563) @@ -13,90 +13,47 @@ #include "pswarp.h" -// Define an input file -static pmFPAfile *defineInputFile(pmConfig *config,// Configuration - pmFPAfile *bind, // File to which to bind, or NULL - char *filerule, // Name of file rule - char *argname, // Argument name - pmFPAfileType fileType // Type of file - ) -{ - bool status; - - pmFPAfile *file = NULL; - // look for the file on the argument list - if (bind) { - file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); - } else { - file = pmFPAfileDefineFromArgs(&status, config, filerule, argname); - } - if (!status) { - psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); - return false; - } - if (!file) { - // look for the file on the RUN metadata - file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return - if (!status) { - psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); - return NULL; - } - } - - if (!file) { - return NULL; - } - - if (file->type != fileType) { - psError(PSWARP_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); - return NULL; - } - - return file; -} - - - +static bool foundAstrom = false; +static bool foundVariance = false; +static bool foundBackground = false; bool pswarpParseCamera(pmConfig *config) { psAssert(config, "Require configuration"); - bool mdok; // Status of MD lookup - - // The input image(s) is required: it defines the camera - pmFPAfile *input = defineInputFile(config, NULL, "PSWARP.INPUT", "INPUT", PM_FPA_FILE_IMAGE); - if (!input) { - psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.INPUT"); - return false; - } - - pmFPAfile *astrom = defineInputFile(config, NULL, "PSWARP.ASTROM", "ASTROM", PM_FPA_FILE_CMF); - psLogMsg("pswarp", PS_LOG_INFO, "Astrometry source: %s", astrom ? "supplied" : "header"); - - pmFPAfile *inMask = defineInputFile(config, input, "PSWARP.MASK", "MASK", PM_FPA_FILE_MASK); - if (!inMask) { - psLogMsg("pswarp", PS_LOG_INFO, "No mask supplied"); - } - - pmFPAfile *inVariance = defineInputFile(config, input, "PSWARP.VARIANCE", "VARIANCE", - PM_FPA_FILE_VARIANCE); - if (!inVariance) { - psLogMsg("pswarp", PS_LOG_INFO, "No variance supplied"); - } - - pmFPAfile *inBackground = defineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND", - PM_FPA_FILE_IMAGE); - if (!inBackground) { - psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied"); - } - + bool status; // Status of MD lookup + + // *** parse the input information (either from -file or from -input) + + // if INPUTS exists, we have a metadata file to provide input, variance, mask, etc + psMetadata *inputFile = psMetadataLookupPtr(&status, config->arguments, "INPUTS"); + if (inputFile) { + if (!pswarpParseMultiInput (config, inputFile)) { + psError(PSWARP_ERR_CONFIG, false, "failed to parse multi-input file"); + return false; + } + } else { + if (!pswarpParseSingleInput (config)) { + psError(PSWARP_ERR_CONFIG, false, "failed to parse multi-input file"); + return false; + } + } + + // once we have read the input mask headers (if any), we can use that to set the mask bits + if (!pswarpSetMaskBits(config)) { + psError(psErrorCodeLast(), false, "failed to set mask bits"); + return false; + } + + // *** parse the output information (output target, output astrometry [skycell]) + // The input skycell is a required argument: it defines the output image - // XXX we may need a different skycell structure here + pmConfig *skyConfig = NULL; pmFPAfile *skycell = NULL; - pmConfig *skyConfig = NULL; - bool status = pswarpDefineSkycell(&skycell, &skyConfig, config, "PSWARP.SKYCELL", "SKYCELL"); - if (!status) { - psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.SKYCELL"); - return false; - } + + skycell = pmFPAfileDefineNewConfig(&status, &skyConfig, config, "PSWARP.SKYCELL", "SKYCELL"); + if (!status || !skycell) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.SKYCELL"); + return false; + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SKYCELL.CAMERA", 0, "Name of camera for skycell", skyConfig->cameraName); psFree(skyConfig); @@ -107,6 +64,11 @@ } + // use the skycell camera or the input camera? + bool useInputCamera = !strcmp (skycell->cameraName, "SIMPLE"); + // The output skycell - pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT"); + pmFPAfile *output = useInputCamera ? + pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT") : + pmFPAfileDefineOutputForFormat(config, NULL, "PSWARP.OUTPUT", skycell->cameraName, skycell->formatName); if (!output) { psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT"); @@ -114,6 +76,9 @@ } output->save = true; - - pmFPAfile *outMask = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.MASK"); + pmFPAAddSourceFromFormat(output->fpa, output->format); // ** builds the HDUs, is this OK? + + pmFPAfile *outMask = useInputCamera ? + pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.MASK"): + pmFPAfileDefineOutputForFormat(config, output->fpa, "PSWARP.OUTPUT.MASK", skycell->cameraName, skycell->formatName); if (!outMask) { psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.MASK"); @@ -122,6 +87,9 @@ outMask->save = true; - if (inVariance) { - pmFPAfile *outVariance = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.VARIANCE"); + // only create an output variance in we supply an input variance + if (foundVariance) { + pmFPAfile *outVariance = useInputCamera ? + pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.VARIANCE"): + pmFPAfileDefineOutputForFormat(config, output->fpa, "PSWARP.OUTPUT.VARIANCE", skycell->cameraName, skycell->formatName); if (!outVariance) { psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.VARIANCE"); @@ -131,6 +99,9 @@ } - if (astrom && psMetadataLookupBool(&mdok, recipe, "SOURCES")) { - pmFPAfile *outSources = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.SOURCES"); + // XXX we assume input sources come from the input astrom description, but this need not be true (we could define an input sources file) + if (foundAstrom && psMetadataLookupBool(&status, recipe, "SOURCES")) { + pmFPAfile *outSources = useInputCamera ? + pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.SOURCES"): + pmFPAfileDefineOutputForFormat(config, output->fpa, "PSWARP.OUTPUT.SOURCES", skycell->cameraName, skycell->formatName); if (!outSources) { psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.SOURCES"); @@ -140,24 +111,25 @@ } - 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 - // psphotDefineFiles - pmFPAfile *psphotInput = pmFPAfileDefineSkycell(config, NULL, "PSPHOT.INPUT"); + // only create an output background if an input background is supplied + if (foundBackground && psMetadataLookupBool(&status, recipe, "BACKGROUND.MODEL")) { + // pmFPAfile *outBackground = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL"); + pmFPAfile *outBackground = useInputCamera ? + pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL"): + pmFPAfileDefineOutputForFormat(config, NULL, "PSWARP.OUTPUT.BKGMODEL", skycell->cameraName, skycell->formatName); + if (!outBackground) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL"); + return false; + } + outBackground->xBin = psMetadataLookupS32(&status, recipe, "BKG.XGRID"); + outBackground->yBin = psMetadataLookupS32(&status, recipe, "BKG.YGRID"); + outBackground->save = true; + pmFPAAddSourceFromFormat(outBackground->fpa, outBackground->format); // ** builds the HDUs, is this OK? + } + + if (psMetadataLookupBool(&status, recipe, "PSF")) { + // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by psphotDefineFiles + pmFPAfile *psphotInput = useInputCamera ? + pmFPAfileDefineSkycell(config, NULL, "PSPHOT.INPUT"): + pmFPAfileDefineOutputForFormat(config, NULL, "PSPHOT.INPUT", skycell->cameraName, skycell->formatName); if (!psphotInput) { psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.INPUT")); @@ -165,8 +137,18 @@ } psphotInput->src = psMemIncrRefCounter(output->fpa); - // specify the number of psphot input images + pmFPAAddSourceFromFormat(psphotInput->fpa, psphotInput->format); // ** builds the HDUs, is this OK? + + // specify the number of psphot input images (psphotReadout loops over all input images) psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1); - pmFPAfile *psphotInSources = pmFPAfileDefineSkycell(config, output->fpa, "PSPHOT.INPUT.CMF"); + // the input sources (read from the input astrometry file) are transformed (in pswarpLoop) to the readout->analysis + // entries of the output file PSWARP.OUTPUT.SOURCES, associated with the main output pmFPAfile PSWARP.OUTPUT + + // the PSPHOT.INPUT.CMF is used to supply those sources to psphot (specifically psphotLoadPSFSources). + + // pmFPAfile *psphotInSources = pmFPAfileDefineSkycell(config, output->fpa, "PSPHOT.INPUT.CMF"); + pmFPAfile *psphotInSources = useInputCamera ? + pmFPAfileDefineSkycell(config, output->fpa, "PSPHOT.INPUT.CMF"): + pmFPAfileDefineOutputForFormat(config, output->fpa, "PSPHOT.INPUT.CMF", skycell->cameraName, skycell->formatName); if (!psphotInSources) { psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.INPUT.CMF")); @@ -195,8 +177,69 @@ psphotSources->save = false; } + + // only keep the necessary (output) camera configs and relevant recipes + const char *skyCamera = psMetadataLookupStr(NULL, config->arguments, "SKYCELL.CAMERA"); + pmConfigCamerasCull(config, skyCamera); + pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); + + psTrace("pswarp", 1, "Done with pswarpParseCamera...\n"); + return true; +} + +bool pswarpParseSingleInput (pmConfig *config) { + + // The input image(s) is required: it defines the camera + pmFPAfile *input = pswarpDefineInputFile(config, NULL, "PSWARP.INPUT", "INPUT", PM_FPA_FILE_IMAGE); + if (!input) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.INPUT"); + return false; + } + + // Define the input astrometry file(s) (may be either WCS or CMF) + bool status = false; + pmFPAfile *astrom = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM", "ASTROM"); + if (!astrom) { + // look for the file on the RUN metadata + astrom = pmFPAfileDefineFromRun(&status, NULL, config, "PSWARP.ASTROM"); // File to return + if (!status) { + psLogMsg("pswarp", PS_LOG_INFO, "Failed to load file definition for %s", "PSWARP.ASTROM"); + } + } + if (astrom) { + if ((astrom->type != PM_FPA_FILE_CMF) && (astrom->type != PM_FPA_FILE_WCS)) { + psLogMsg("pswarp", PS_LOG_INFO, "%s is neither CMF nor WCS", "PSWARP.ASTROM"); + // XXX raise an error here?? + } else { + foundAstrom = true; + } + } + psLogMsg("pswarp", PS_LOG_INFO, "Astrometry source: %s", astrom ? "supplied" : "header"); + + // Define the input mask file(s) + pmFPAfile *inMask = pswarpDefineInputFile(config, input, "PSWARP.MASK", "MASK", PM_FPA_FILE_MASK); + if (!inMask) { + psLogMsg("pswarp", PS_LOG_INFO, "No mask supplied"); + } + + // Define the input variance file(s) + pmFPAfile *inVariance = pswarpDefineInputFile(config, input, "PSWARP.VARIANCE", "VARIANCE", PM_FPA_FILE_VARIANCE); + if (inVariance) { + foundVariance = true; + } else { + psLogMsg("pswarp", PS_LOG_INFO, "No variance supplied"); + } + + pmFPAfile *inBackground = pswarpDefineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND", PM_FPA_FILE_IMAGE); + if (inBackground) { + foundBackground = true; + } else { + // cannot do the background model if an input is not supplied. + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "no input background supplied", false); + psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied"); + } // Chip selection: turn on only the chips specified - char *chipLine = psMetadataLookupStr(&mdok, config->arguments, "CHIP_SELECTIONS"); - if (mdok) { + char *chipLine = psMetadataLookupStr(&status, config->arguments, "CHIP_SELECTIONS"); + if (status) { psArray *chips = psStringSplitArray (chipLine, ",", false); if (chips->n > 0) { @@ -213,5 +256,140 @@ } - psTrace("pswarp", 1, "Done with pswarpParseCamera...\n"); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NUM_INPUTS", PS_META_REPLACE, "input file sets", 1); return true; } + +// read input file information from a metadata file +// XXX for now, in the multi-input format, we disable PSF, BACKGROUND, and CHIP_SELECTIONS +bool pswarpParseMultiInput (pmConfig *config, psMetadata *fileListMD) { + + // the multi-input file consists of a set of metadata blocks, each with entries + // for INPUT and (possibly) ASTROM, MASK, VARIANCE, ??? + + bool status = false; + + for (int i = 0; i < fileListMD->list->n; i++) { + psMetadataItem *item = psMetadataGet(fileListMD, i); + if (item->type != PS_DATA_METADATA) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type METADATA", item->name); + return false; + } + + // next input metadata block of interest + psMetadata *fileBlockMD = item->data.md; + + psString inputName = psMetadataLookupStr(&status, fileBlockMD, "INPUT"); // Name of image + if (!inputName || !strlen(inputName)) { + psError(PS_ERR_UNKNOWN, false, "input file not found in metadata block %d (%s)", i, item->name); + return false; + } + AddStringAsArray (config->arguments, inputName, "FILENAMES"); + + pmFPAfile *input = pswarpDefineInputFile (config, NULL, "PSWARP.INPUT", "FILENAMES", PM_FPA_FILE_IMAGE); + if (!input) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.INPUT"); + return false; + } + + psString astromName = psMetadataLookupStr(&status, fileBlockMD, "ASTROM"); // Name of astrom file + if (astromName && strlen(astromName)) { + AddStringAsArray (config->arguments, astromName, "FILENAMES"); + + pmFPAfile *astrom = pmFPAfileDefineFromArgs (&status, config, "PSWARP.ASTROM", "FILENAMES"); + if (!status) { + psLogMsg("pswarp", PS_LOG_INFO, "Failed to load file definition for %s", "PSWARP.ASTROM"); + } + if (astrom) { + if ((astrom->type != PM_FPA_FILE_CMF) && (astrom->type != PM_FPA_FILE_WCS)) { + psLogMsg("pswarp", PS_LOG_INFO, "%s is neither CMF nor WCS", "PSWARP.ASTROM"); + } else { + foundAstrom = true; + } + } + psLogMsg("pswarp", PS_LOG_INFO, "Astrometry source: %s", astrom ? "supplied" : "header"); + } + + pmFPAfile *mask = NULL; + psString maskName = psMetadataLookupStr(&status, fileBlockMD, "MASK"); // Name of mask + if (maskName && strlen(maskName)) { + AddStringAsArray (config->arguments, maskName, "FILENAMES"); + + mask = pswarpDefineInputFile (config, input, "PSWARP.MASK", "FILENAMES", PM_FPA_FILE_MASK); + } + if (!mask) { + psLogMsg("pswarp", PS_LOG_INFO, "No mask supplied"); + } + + pmFPAfile *variance = NULL; + psString varianceName = psMetadataLookupStr(&status, fileBlockMD, "VARIANCE"); // Name of variance + if (varianceName && strlen(varianceName)) { + AddStringAsArray (config->arguments, varianceName, "FILENAMES"); + + variance = pswarpDefineInputFile (config, input, "PSWARP.VARIANCE", "FILENAMES", PM_FPA_FILE_VARIANCE); + } + if (variance) { + foundVariance = true; + } else { + psLogMsg("pswarp", PS_LOG_INFO, "No variance supplied"); + } + + // XXX for now we do not have an option to supply a background model set in multi input mode + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "no input background supplied", false); + } + + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NUM_INPUTS", PS_META_REPLACE, "input file sets", fileListMD->list->n); + return true; +} + +// Define an input file based on name in config->arguments or config RUN block +pmFPAfile *pswarpDefineInputFile(pmConfig *config,// Configuration + pmFPAfile *bind, // File to which to bind, or NULL + char *filerule, // Name of file rule + char *argname, // Argument name + pmFPAfileType fileType // Type of file + ) +{ + bool status; + + pmFPAfile *file = NULL; + // look for the file on the argument list + if (bind) { + file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); + } else { + file = pmFPAfileDefineFromArgs(&status, config, filerule, argname); + } + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return false; + } + if (!file) { + // look for the file on the RUN metadata + file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + } + + if (!file) { + return NULL; + } + + if (file->type != fileType) { + psError(PSWARP_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + return file; +} + +bool AddStringAsArray (psMetadata *md, char *string, char *name) { + + psArray *dummy = psArrayAlloc(1); // Dummy array of filenames for this FPA + dummy->data[0] = psStringCopy(string); + + psMetadataAddArray(md, PS_LIST_TAIL, name, PS_META_REPLACE, "string added as array", dummy); + psFree(dummy); + return true; +} + Index: /trunk/pswarp/src/pswarpPixelsLit.c =================================================================== --- /trunk/pswarp/src/pswarpPixelsLit.c (revision 35562) +++ /trunk/pswarp/src/pswarpPixelsLit.c (revision 35563) @@ -64,26 +64,18 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue)) { - if (y > yMax) { - yMax = y; - } - if (y < yMin) { - yMin = y; - } - if (x > xMax) { - xMax = x; - } - if (x < xMin) { - xMin = x; - } - } + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue) { continue; } + xMin = PS_MIN (xMin, x); + xMax = PS_MAX (xMax, x); + yMin = PS_MIN (yMin, y); + yMax = PS_MAX (yMax, y); } } if (stats) { - psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMIN", 0, "Minimum valid x value", xMin); - psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMAX", 0, "Maximum valid x value", xMax); - psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMIN", 0, "Minimum valid y value", yMin); - psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMAX", 0, "Maximum valid y value", yMax); + // XXX with multiple inputs (eg, output stacks -> exposure), these only represent the last input + psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMIN", PS_META_REPLACE, "Minimum valid x value", xMin); + psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMAX", PS_META_REPLACE, "Maximum valid x value", xMax); + psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMIN", PS_META_REPLACE, "Minimum valid y value", yMin); + psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMAX", PS_META_REPLACE, "Maximum valid y value", yMax); } Index: /trunk/pswarp/src/pswarpStatsFile.c =================================================================== --- /trunk/pswarp/src/pswarpStatsFile.c (revision 35563) +++ /trunk/pswarp/src/pswarpStatsFile.c (revision 35563) @@ -0,0 +1,94 @@ +/** @file pswarpStatsFile.c + * + * @brief functions for managing the stats file (or return NULL, or exit with an error) + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pswarp.h" + +static void pswarpStatsFileFree (pswarpStatsFile *statsFile) { + psFree (statsFile->name); + psFree (statsFile->md); +} + +pswarpStatsFile *pswarpStatsFileAlloc () { + + pswarpStatsFile *statsFile = psAlloc(sizeof(pswarpStatsFile)); + psMemSetDeallocator(statsFile, (psFreeFunc)pswarpStatsFileFree); + + statsFile->f = NULL; + statsFile->name = NULL; + statsFile->md = NULL; + + return statsFile; +} + +// Open the statistics file +pswarpStatsFile *pswarpStatsFileOpen (pmConfig *config) { + + bool mdok; + + char *name = psMetadataLookupStr(&mdok, config->arguments, "STATS"); ///< Filename for statistics + + pswarpStatsFile *statsFile = pswarpStatsFileAlloc (); + statsFile->name = psMemIncrRefCounter(name); + + if (!mdok) return statsFile; // XXX this is probably a config error, but treat as ok for now + if (!name) return statsFile; + if (strlen(name) == 0) return statsFile; + + psString resolved = pmConfigConvertFilename(name, config, true, true); + + statsFile->f = fopen(resolved, "w"); + if (!statsFile->f) { + psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); + psFree(resolved); + pswarpCleanup(config, statsFile); + } + psFree(resolved); + + statsFile->md = psMetadataAlloc(); + psMetadataAddS32(statsFile->md, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0); + + return statsFile; +} + +bool pswarpStatsFileSave (pmConfig *config, pswarpStatsFile *statsFile) { + + if (!statsFile) return true; + if (!statsFile->md) return true; + + // Write out summary statistics + psMetadataAddF32(statsFile->md, PS_LIST_TAIL, "DT_WARP", 0, "Time for warp completion", psTimerMark("pswarp")); + + // convert the psMetadata to a string block + const char *statsMDC = psMetadataConfigFormat(statsFile->md); + if (!statsMDC) { + psError(psErrorCodeLast(), false, "Unable to get statistics file."); + return false; + } + + if (fprintf(statsFile->f, "%s", statsMDC) != strlen(statsMDC)) { + psError(PSWARP_ERR_IO, true, "Unable to write statistics file."); + return false; + } + psFree(statsMDC); + + if (fclose(statsFile->f) == EOF) { + psError(PSWARP_ERR_IO, true, "Unable to close statistics file."); + statsFile->f = NULL; + return false; + } + statsFile->f = NULL; + + pmConfigRunFilenameAddWrite(config, "STATS", statsFile->name); + psFree(statsFile); + + return true; +} + Index: /trunk/pswarp/src/pswarpTransformReadout.c =================================================================== --- /trunk/pswarp/src/pswarpTransformReadout.c (revision 35562) +++ /trunk/pswarp/src/pswarpTransformReadout.c (revision 35563) @@ -16,5 +16,5 @@ * NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT */ -bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config) +bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config, bool backgroundWarp) { // XXX this implementation currently ignores the use of the region @@ -52,8 +52,4 @@ } - // XXX unused int nThreads = psMetadataLookupS32(&mdok, config->arguments, "NTHREADS"); // Number of threads - // XXX unused if (!mdok) { - // XXX unused nThreads = 0; - // XXX unused } float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); ///< Flux fraction for "poor" @@ -61,5 +57,4 @@ // 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 @@ -117,4 +112,6 @@ psAssert (xGridMax < grid->nXpts, "xGridMax too big\n"); psAssert (yGridMax < grid->nYpts, "yGridMax too big\n"); + + // fprintf (stderr, "warp %d,%d - %d,%d\n", xGridMin, yGridMin, xGridMax, yGridMax); // create jobs and supply them to the threads @@ -132,10 +129,9 @@ args->goodPixels = 0; - if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { + if (backgroundWarp) { 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 @@ -210,11 +206,9 @@ psFree(interp); - 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; + if (goodPixels > 0 && !backgroundWarp && psMetadataLookupBool(&mdok, recipe, "SOURCES")) { + if (!pswarpTransformSources(output, input, config)) { + psError(psErrorCodeLast(), false, "Unable to transform sources."); + return false; } - } } Index: /trunk/pswarp/src/pswarpTransformTile.c =================================================================== --- /trunk/pswarp/src/pswarpTransformTile.c (revision 35562) +++ /trunk/pswarp/src/pswarpTransformTile.c (revision 35563) @@ -89,4 +89,7 @@ for (int y = yMin; y < yMax; y++) { + + int yOut = y - outRow0; ///< Position on output image + for (int x = xMin; x < xMax; x++) { // Only transform those pixels requested @@ -94,4 +97,11 @@ continue; } + + int xOut = x - outCol0; + + // XXX the existing image may already have valid data -- probably should keep + // the best, but for the moment, just keep the pixel which is not NAN + if (isfinite(outImageData[yOut][xOut])) continue; + // pswarpMapApply converts the output coordinate (x,y) to the input coordinate. // both are in the parent frames of the input and output images. @@ -127,8 +137,7 @@ } - int xOut = x - outCol0, yOut = y - outRow0; ///< Position on output image - if (outImageData) { - outImageData[yOut][xOut] = imageValue * jacobian; + // XXX TEST outImageData[yOut][xOut] = value; + outImageData[yOut][xOut] = imageValue * jacobian; } if (outVarData) { Index: /trunk/pswarp/src/pswarpUpdateMetadata.c =================================================================== --- /trunk/pswarp/src/pswarpUpdateMetadata.c (revision 35563) +++ /trunk/pswarp/src/pswarpUpdateMetadata.c (revision 35563) @@ -0,0 +1,145 @@ +/** @file pswarpUpdateMetadata.c + * + * @brief update generic metadata info for warp + * @ingroup pswarp + * + * @author IfA + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pswarp.h" +#define WCS_NONLIN_TOL 0.001 // Non-linear tolerance for header WCS + +// once the output fpa elements have been built, loop over the fpa and generate metadata +bool pswarpUpdateMetadata (pmFPA *output, pmFPA *skycell, pmFPA *input, pmFPA *astrom, pmConfig *config, bool fullImage) { + + pmFPAview *view = pmFPAviewAlloc(0); + + // check output astrometry mode + bool bilevelAstrometry = psMetadataLookupBool (NULL, skycell->analysis, "ASTROMETRY.BILEVEL"); + + pmChip *chip; + 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; } + + // use this in output metadata info (MD5 sums) + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + + pmCell *cell; + 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; } + + // use this in output metadata info (MD5 sums) + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) { + // skip empty output readouts + if (!readout->data_exists) continue; + + if (fullImage) { + // Set covariance matrix for output + bool status = false; + psList *covariances = psMetadataLookupPtr(&status, readout->analysis, PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices + psAssert(covariances, "Should be there"); + psArray *covars = psListToArray(covariances); // Array of covariance matrices + if (covars->n) { + psKernel *covar = psImageCovarianceAverage(covars); + psMetadataRemoveKey(readout->analysis, PSWARP_ANALYSIS_COVARIANCES); + + // Correct covariance matrix scale for the mean (square root of the) Jacobian + double jacobian = psMetadataLookupF64(NULL, readout->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian + int goodPixels = psMetadataLookupS32(NULL, readout->analysis, PSWARP_ANALYSIS_GOODPIX); // Good pixels + jacobian /= goodPixels; + readout->covariance = psImageCovarianceScale(covar, jacobian); + psFree(covar); + + if (readout->variance) { + psImageCovarianceTransfer(readout->variance, readout->covariance); + } + + psFree(covars); + } + } + + // Add MD5 information for readout + psString headerName = NULL; ///< Header name for MD5 + psVector *md5 = psImageMD5(readout->image); ///< md5 hash + psString md5string = psMD5toString(md5); ///< String + + pmHDU *hdu = pmHDUFromReadout(readout); + psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, "Image MD5", md5string); + psFree(md5); + psFree(md5string); + psFree(headerName); + } + + if (fullImage) { + psList *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "INPUT.CELLS"); + if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) { + psError(psErrorCodeLast(), false, "Unable to average cell concepts."); + psFree(view); + return false; + } + } + + // XXX Is this too ad-hoc? + psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); ///< Trim section + trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels + } + + // update astrometry headers + pmHDU *hdu = pmHDUFromChip(chip); ///< HDU for source + if (!hdu || !hdu->header) { + psError(PM_ERR_PROG, false, "Unable to find header for output."); + psFree(view); + return false; + } + + if (bilevelAstrometry) { + if (!pmAstromWriteBilevelChip(hdu->header, chip, WCS_NONLIN_TOL)) { + 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 (!pmAstromWriteWCS(hdu->header, output, chip, WCS_NONLIN_TOL)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell."); + psFree(view); + return false; + } + } + } + + if (!psMetadataCopy(output->concepts, input->concepts)) { + psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output."); + return false; + } + + // Update ZP from the astrometry + if (fullImage) { + psMetadataItem *item = psMetadataLookup(output->concepts, "FPA.ZP"); + item->data.F32 = psMetadataLookupF32(NULL, astrom->concepts, "FPA.ZP"); + } + + // apply the bilevel astrometry elements to the target + if (bilevelAstrometry) { + pmHDU *phu = pmFPAviewThisPHU(view, output); ///< Astrometry PHU + if (!phu->header) { + phu->header = psMetadataAlloc (); + } + if (!pmAstromWriteBilevelMosaic(phu->header, output, WCS_NONLIN_TOL)) { + psError(psErrorCodeLast(), false, "Unable to write bilevel mosaic astrometry for skycell."); + psFree(view); + return false; + } + } + + return true; +} Index: /trunk/pswarp/src/pswarpUpdateStatistics.c =================================================================== --- /trunk/pswarp/src/pswarpUpdateStatistics.c (revision 35563) +++ /trunk/pswarp/src/pswarpUpdateStatistics.c (revision 35563) @@ -0,0 +1,121 @@ +/** @file pswarpUpdateStatistics.c + * + * @brief generate output statistics + * @ingroup pswarp + * + * @author IfA + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pswarp.h" + +// once the output fpa elements have been built, loop over the fpa and generate stats +// for each readout +bool pswarpUpdateStatistics (pmFPA *output, psMetadata *fpaStats, pmFPA *input, pmFPA *astrom, pmConfig *config) { + + if (!fpaStats) { + psLogMsg("pswarp", PS_LOG_INFO, "stats not requested, skipping"); + return true; + } + + // load the recipe + bool status = false; + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE); + return false; + } + + // output mask bits + psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT"); + psAssert (status, "MASK.OUTPUT was not defined"); + + int Nreadout = 0; // count the number of readouts with data + + bool doMaskStats = psMetadataLookupBool(&status, recipe, "MASK.STATS"); + + pmFPAview *view = pmFPAviewAlloc(0); + + pmChip *chip; + 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; } + + // use this in output metadata info (MD5 sums) + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + + psMetadata *chipStats = psMetadataLookupPtr (&status, fpaStats, chipName); + if (!chipStats) { + chipStats = psMetadataAlloc (); + psMetadataAdd (fpaStats, PS_LIST_TAIL, chipName, PS_DATA_METADATA, "chip stats folder", chipStats); + psFree (chipStats); + } + + pmCell *cell; + 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; } + + // use this in output metadata info (MD5 sums) + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); + + psMetadata *cellStats = psMetadataLookupPtr (&status, chipStats, cellName); + if (!cellStats) { + cellStats = psMetadataAlloc (); + psMetadataAdd (chipStats, PS_LIST_TAIL, cellName, PS_DATA_METADATA, "cell stats folder", cellStats); + psFree (cellStats); + } + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) { + // skip empty output readouts + if (!readout->data_exists) { + Nreadout ++; + continue; + } + + psString readoutName = NULL; + psStringAppend (&readoutName, "READOUT.%02d", view->readout); + + psMetadata *readoutStats = psMetadataLookupPtr (&status, cellStats, readoutName); + if (!readoutStats) { + readoutStats = psMetadataAlloc (); + psMetadataAdd (cellStats, PS_LIST_TAIL, readoutName, PS_DATA_METADATA, "readout stats folder", readoutStats); + psFree (readoutStats); + } + + if (!pswarpPixelsLit(readout, readoutStats, config)) { + psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); + psFree(view); + return false; + } + + if (doMaskStats && !pswarpMaskStats(readout, readoutStats, config)) { + psError(psErrorCodeLast(), false, "Unable to calculate mask stats."); + psFree(view); + return false; + } + + psFree (readoutName); + } + } + } + + // Perform statistics on the output image (ppStatsFPA loops down the hierarchy) + pmFPAviewReset (view); + if (!ppStatsFPA(fpaStats, output, view, maskValue, config)) { + psWarning("Unable to perform statistics on warped image."); + } + psFree(view); + + if (Nreadout == 0) { + psWarning("No overlap between input and skycell."); + if (fpaStats) { + psMetadataAddS32(fpaStats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP); + } + } + + return true; +} Index: /trunk/pswarp/src/pswarpVersion.c =================================================================== --- /trunk/pswarp/src/pswarpVersion.c (revision 35562) +++ /trunk/pswarp/src/pswarpVersion.c (revision 35563) @@ -13,10 +13,4 @@ #include #endif - -#include -#include -#include -#include -#include #include "pswarp.h"