Index: /branches/eam_branches/ipp-20130419/pswarp/doc/notes.20130406.txt =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/doc/notes.20130406.txt (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/doc/notes.20130406.txt (revision 35434) @@ -10,4 +10,11 @@ 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 @@ -91,2 +98,8 @@ * 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. + Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarp.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarp.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarp.c (revision 35434) @@ -39,5 +39,5 @@ } -# if (0) +# if (1) // XXX I'm going to add the new code before removing the old code. // load the input & output astrometry, find the output overlaps, generate the output pixels @@ -46,8 +46,17 @@ } +# if (0) pswarpDumpOutput (config); exit (0); # endif + // load and warp + if (!pswarpLoop(config, statsFile->md)) { + pswarpCleanup(config, statsFile); + } + pswarpDumpOutput (config); + exit (0); + +# else // load the skycell layout information: // read the output astrometry description @@ -60,7 +69,6 @@ pswarpCleanup(config, statsFile); } - // load and warp - if (!pswarpLoop(config, statsFile->md)) { + if (!pswarpLoopSkycell(config, statsFile->md)) { pswarpCleanup(config, statsFile); } @@ -69,5 +77,5 @@ pswarpCleanup(config, statsFile); } - +# endif psLogMsg("pswarp", PS_LOG_INFO, "complete pswarp run: %f sec\n", psTimerMark("pswarp")); pswarpCleanup(config, statsFile); Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarp.h =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarp.h (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarp.h (revision 35434) @@ -194,3 +194,7 @@ pswarpBounds *pswarpBoundsAlloc(); -bool pswarpLoadAstrometry (pmConfig *config, pmFPAfile *astrom, pmFPAfile *target); +bool pswarpLoadAstrometry (pmFPAfile *target, pmFPAfile *astrom, pmConfig *config); + +bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config); +bool pswarpMakePSF (pmFPA *output); +bool pswarpUpdateStatistics (pmFPA *output, pmFPA *input, pmConfig *config); Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpDefineLayout.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpDefineLayout.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpDefineLayout.c (revision 35434) @@ -25,5 +25,5 @@ astrom = input; } - if (!pswarpLoadAstrometry (config, input, astrom)) { + if (!pswarpLoadAstrometry (input, astrom, config)) { psError(PSWARP_ERR_CONFIG, false, "problem loading input astrometry\n"); return false; @@ -41,5 +41,5 @@ return false; } - if (!pswarpLoadAstrometry (config, output, skycell)) { + if (!pswarpLoadAstrometry (output, skycell, config)) { psError(PSWARP_ERR_CONFIG, false, "problem loading output astrometry\n"); return false; @@ -50,11 +50,18 @@ // find the R,D center for the skycell, use for common projection - psProjection *frame = pswarpLocalFrame (output->fpa); + psProjection *frame = pswarpLocalFrame (skycell->fpa); // generate Lmin,max, Mmin,max for both datasets + fprintf (stderr, "srcBounds:\n"); + // pswarpBounds *srcBounds = pswarpMakeBounds (astrom->fpa, frame); pswarpBounds *srcBounds = pswarpMakeBounds (input->fpa, frame); - pswarpBounds *tgtBounds = pswarpMakeBounds (output->fpa, frame); + + fprintf (stderr, "srcBounds:\n"); + pswarpBounds *tgtBounds = pswarpMakeBounds (skycell->fpa, frame); + pswarpFindOverlap (output->fpa, srcBounds, tgtBounds); + // XXX the loop below generates the output pixels. 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); @@ -63,12 +70,14 @@ 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 (!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 || !cell->file_exists) { continue; } + if (!cell->process) { continue; } + + pmCell *refcell = pmFPAviewThisCell(view, skycell->fpa); ///< Target cell // we've got the output astrom header - pmHDU *hdu = pmHDUFromCell(cell); ///< HDU for source + pmHDU *hdu = pmHDUFromCell(refcell); ///< HDU for source if (!hdu || !hdu->header) { psError(PM_ERR_PROG, false, "Unable to find header for sky cell."); @@ -92,5 +101,4 @@ // copy the image concepts from the skycell bool status = false; - pmCell *refcell = pmFPAviewThisCell(view, skycell->fpa); ///< Target cell psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XBIN"); psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YBIN"); @@ -101,4 +109,16 @@ psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.X0"); psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.Y0"); + + chip->data_exists = true; + chip->file_exists = true; + cell->data_exists = true; + cell->file_exists = true; + +# if (0) + // ******* this chunk comes from pswarpLoop, but should be done when the output chip is generated + pmHDU *outHDU = pmHDUFromCell (cell); ///< HDU for the output warped image + outHDU->header = psMetadataCopy(outHDU->header, hdu->header); + pswarpVersionHeader(outHDU->header); +# endif } } Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoadAstrometry.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoadAstrometry.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoadAstrometry.c (revision 35434) @@ -14,8 +14,10 @@ // WCS and read the headers. We place the resulting info on the target astrometry // containers -bool pswarpLoadAstrometry (pmConfig *config, pmFPAfile *astrom, pmFPAfile *target) { +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 @@ -42,5 +44,5 @@ 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 || !cell->file_exists) { continue; } + if (!cell->process) { continue; } if (!pmFPAfileRead (astrom, view, config)) { psError(psErrorCodeLast(), false, "failed READ at CELL %s", astrom->name); @@ -89,5 +91,5 @@ // apply the bilevel astrometry elements to the target if (bilevelAstrometry) { - if (!pmAstromReadBilevelMosaic(target->fpa, phu->header)) { + if (!pmAstromReadBilevelMosaic(astrom->fpa, phu->header)) { psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell."); psFree(view); @@ -95,4 +97,8 @@ } } + // 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) { @@ -110,6 +116,8 @@ // 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(targetChip, hdu->header)) { + if (!pmAstromReadBilevelChip(chip, hdu->header)) { psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell."); psFree(view); @@ -118,5 +126,5 @@ } else { // we use a default FPA pixel scale of 1.0 - if (!pmAstromReadWCS(target->fpa, targetChip, hdu->header, 1.0)) { + if (!pmAstromReadWCS(target->fpa, chip, hdu->header, 1.0)) { psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell."); psFree(view); @@ -124,4 +132,6 @@ } } + targetChip->toFPA = psMemIncrRefCounter (chip->toFPA); + targetChip->fromFPA = psMemIncrRefCounter (chip->fromFPA); } Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoop.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoop.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoop.c (revision 35434) @@ -20,7 +20,4 @@ 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 @@ -29,4 +26,5 @@ // load the recipe + bool status = false; psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE); if (!recipe) { @@ -35,12 +33,9 @@ } + // XXX this should probably happen elsewhere (logically it belongs to file I/O, right?) 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 @@ -62,20 +57,21 @@ } - // 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"); + // select the input data sources + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "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?? + + // use the external astrometry source if supplied + pmFPAfile *skycell = psMetadataLookupPtr(NULL, config->files, "PSWARP.SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry."); + return false; + } + + // XXX NOTE : 'output' was a readout, now we need to loop over all output targets + + // ensure everyone is off except what we need below pswarpFileActivation(config, detectorFiles, false); pswarpFileActivation(config, photFiles, false); @@ -83,80 +79,16 @@ 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 + // XXX NOTE : do not tackle source issues yet... + // XXX pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES"); + + // Done with the detector side of things + pswarpFileActivation(config, detectorFiles, true); + pswarpFileActivation(config, independentFiles, true); + + pmFPAview *view = pmFPAviewAlloc(0); + + // XXX : need to generate a list of input cells for each output cell + // XXX : psList *cells = psListAlloc(NULL); // List of cells, for concepts averaging // files associated with the science image @@ -166,5 +98,7 @@ } - // *** main transformation block + // *** 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) { @@ -176,23 +110,4 @@ } - // 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) { @@ -203,6 +118,4 @@ goto DONE; } - - psListAdd(cells, PS_LIST_TAIL, cell); // process each of the readouts @@ -217,4 +130,6 @@ } +/// XXX ignore sources for now +# if (0) // Copy the detections from the astrometry carrier to the input, so they can be accessed by // pswarpTransformReadout @@ -224,6 +139,7 @@ psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections); } - - pswarpTransformReadout(output, readout, config); +# endif + + pswarpTransformToTarget(output->fpa, readout, config); if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { @@ -242,105 +158,4 @@ } } - - 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."); @@ -352,7 +167,164 @@ 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. + // XXXX turn these on later .. pswarpUpdateStatistics (output); + // XXXX turn these on later .. pswarpMakePSF (output); + + // XXX have not yet written the output to disk + + psFree(view); + + DONE: + return true; +} + +# if (0) +// once the output fpa elements have been built, loop over the fpa and generate stats +// for each readout +bool pswarpUpdateStatistics (pmFPA *output, pmFPA *input, pmConfig *config) { + + // output mask bits + psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT"); + psAssert (status, "MASK.OUTPUT was not defined"); + + 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; } + + 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; } + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, output->fpa, 1)) != NULL) { + if (!readout->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; + } + + if (!pswarpPixelsLit(readout, stats, config)) { + psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); + psFree(view); + goto DONE; + } + bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS"); + if (doStats) { + if (!pswarpMaskStats(readout, 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); + } + } + } + + // XXX need to get the input Cells for this output cell (save on the analysis somewhere + // XXX skip this for now +# if (0) + psArray *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "something"); + if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) { + psError(psErrorCodeLast(), false, "Unable to average cell concepts."); + psFree(view); + goto DONE; + } +# endif + // 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 + + } + + // Perform statistics on the output image + if (stats) { + // XXX this should be done at each level... + if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) { + psWarning("Unable to perform statistics on warped image."); + } + } + } + if (!psMetadataCopy(output->concepts, input->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"); + } + + // need to update the chip and fpa level astrometry appropriate, see + // code in psastro + +# if (0) + // 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); +# endif + +# if (0) + if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) { + psError(psErrorCodeLast(), false, "Unable to generate WCS header."); + goto DONE; + } +# endif + return true; +} +# endif + +# if (0) +// 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) { + + // 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(&mdok, recipe, "PSF")) { @@ -419,29 +391,32 @@ 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; } +# endif + +// 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) { + + 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; } + + 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; } + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) { + pswarpTransformReadout (readout, input, config); + } + } + } + return true; +} + Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoopSkycell.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoopSkycell.c (revision 35434) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoopSkycell.c (revision 35434) @@ -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 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 + 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: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpMatchRange.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpMatchRange.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpMatchRange.c (revision 35434) @@ -98,5 +98,5 @@ // demo forward and backward transformation -# if (0) +# if (1) srcPix->x = *minX; srcPix->y = *minY; @@ -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: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpOverlaps.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpOverlaps.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpOverlaps.c (revision 35434) @@ -83,4 +83,5 @@ } } + fprintf (stderr, "bounds %d : %f - %f, %f - %f\n", view->chip, Pmin, Pmax, Qmin, Qmax); pswarpBoundsAppend (bounds, Pmin, Pmax, Qmin, Qmax); } @@ -109,8 +110,14 @@ 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->process || !chip->file_exists) { continue; } + // 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 @@ -153,5 +160,5 @@ double Dmid = stats->sampleMean; - psProjection *frame = psProjectionAlloc (Rmid, Dmid, 1.0, 1.0, PS_PROJ_TAN); + psProjection *frame = psProjectionAlloc (Rmid, Dmid, PS_RAD_DEG/3600.0, PS_RAD_DEG/3600.0, PS_PROJ_TAN); return frame; } @@ -162,17 +169,14 @@ // elements overlap it. - pmFPAview *view = pmFPAviewAlloc (0); + assert (fpa->chips->n == tgt->Pmin->n); - pmChip *chip = NULL; - int Nout = 0; - 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) { continue; } + for (int j = 0; j < fpa->chips->n; j++) { + pmChip *chip = fpa->chips->data[j]; // we have src bounds - float Pmin = tgt->Pmin->data.F32[Nout]; - float Pmax = tgt->Pmax->data.F32[Nout]; - float Qmin = tgt->Qmin->data.F32[Nout]; - float Qmax = tgt->Qmax->data.F32[Nout]; + 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]; bool valid = false; @@ -183,4 +187,9 @@ valid = Pvalid && Qvalid; } + // if (valid) { + // fprintf (stderr, "chip %d is valid\n", j); + // } else { + // fprintf (stderr, "chip %d is NOT valid\n", j); + // } chip->process = valid; } Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpParseCamera.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpParseCamera.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpParseCamera.c (revision 35434) @@ -128,5 +128,10 @@ // The output skycell - pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT"); + // XXX The output needs to be in a more generic format, not just skycell. + // I need to understand how badly things depend on SKYCELL being a special format + // (is it just pswarp, or is it everything downstream as well??) + // XXX for a test, generate a generic output (based on output target astrometry) + // XXX pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT"); + pmFPAfile *output = pmFPAfileDefineOutputForFormat(config, NULL, "PSWARP.OUTPUT", skycell->cameraName, skycell->formatName); if (!output) { psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT"); @@ -135,4 +140,6 @@ output->save = true; + // XXX for a test, skip these +# if (0) pmFPAfile *outMask = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.MASK"); if (!outMask) { @@ -159,4 +166,5 @@ outSources->save = true; } +# endif if (inBackground && psMetadataLookupBool(&mdok, recipe, "BACKGROUND.MODEL")) { @@ -172,5 +180,6 @@ } - if (psMetadataLookupBool(&mdok, recipe, "PSF")) { + /// XXX TEST + if (false && psMetadataLookupBool(&mdok, recipe, "PSF")) { // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by // psphotDefineFiles Index: /branches/eam_branches/ipp-20130419/pswarp/src/pswarpTransformReadout.c =================================================================== --- /branches/eam_branches/ipp-20130419/pswarp/src/pswarpTransformReadout.c (revision 35433) +++ /branches/eam_branches/ipp-20130419/pswarp/src/pswarpTransformReadout.c (revision 35434) @@ -117,4 +117,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