Index: /trunk/psastro/src/Makefile.am =================================================================== --- /trunk/psastro/src/Makefile.am (revision 19509) +++ /trunk/psastro/src/Makefile.am (revision 19510) @@ -50,5 +50,4 @@ psastroConvert.c \ psastroChipAstrom.c \ - psastroOneChip.c \ psastroOneChipGrid.c \ psastroOneChipFit.c \ Index: unk/psastro/src/psastroEnforceChips.c =================================================================== --- /trunk/psastro/src/psastroEnforceChips.c (revision 19509) +++ (revision ) @@ -1,152 +1,0 @@ -# include "psastroInternal.h" -# define NONLIN_TOL 0.001 /* tolerance in pixels */ - -XXX add a function to supply for all entries, as opposed to fixing the poor astrometry ones. -bool psastroEnforceChips (pmConfig *config, psMetadata *recipe) { - - bool status; - - bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS"); - if (!status) { - fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); - } - if (!fixChips) return true; - - // identify reference astrometry table - // if not defined, correction was not requested; skip step - pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL"); - if (!astrom) return true; - - // select the input data sources - pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); - if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; - } - - // make sure the astrometry model is loaded - // de-activate all files except PSASTRO.MODEL - // XXX do I need to supply the input->fpa->concepts to the astrom->fpa->concepts? - psFree (astrom->fpa->concepts); - astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts); - pmFPAfileActivate (config->files, false, NULL); - pmFPAfileActivate (config->files, true, "PSASTRO.MODEL"); - - pmFPAview *view = pmFPAviewAlloc (0); - - // files associated with the science image - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError (PS_ERR_IO, false, "Can't load the astrometry model file"); - return false; - } - - // loop over all chips, replace input astrometry elements with those from astrom - pmChip *obsChip = NULL; - while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); - if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - - // set the chip astrometry using the astrom file - pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); - - psFree (obsChip->toFPA); - psFree (obsChip->fromFPA); - - // XXX can we do a straight copy? I think so, but we may need to apply a scaling factor based - // on the toSky or toTPA values - obsChip->toFPA = psMemIncrRefCounter (refChip->toFPA); - obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA); - - // XXX test to write out adjusted header - pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL); - } - - psFree (input->fpa->toSky); - psFree (input->fpa->toTPA); - psFree (input->fpa->fromTPA); - input->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky); - input->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA); - input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA); - - psMetadata *updates = psMetadataAlloc(); - pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL); - psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); - psFree (updates); - - psFree (view); - return true; -} - -# if (0) - - // physical pixel scale in microns per pixel - double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE"); - if (!status) { - psError(PS_ERR_IO, true, "Failed to lookup pixel tolerance"); - psFree (view); - return false; - } - double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE"); - if (!status) { - psError(PS_ERR_IO, true, "Failed to lookup angle tolerance"); - psFree (view); - return false; - } - - pmChip *obsChip = NULL; - - // loop over all chips - while ((obsChip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { - psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); - if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - - // set the chip astrometry using the astrom file - pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); - - // bad Astrometry test: ref pixel or angle outside nominal - - psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; - psPlane obsCoord, refCoord; - - // find location of 0,0 pixel in focal plane coords for this chip - psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); - - // find location of 0,0 pixel in focal plane coords for ref chip - psPlaneTransformApply (&refCoord, refChip->toFPA, &refPixel); - - psPlane offPixel = {0.0, 0.0, 0.0, 0.0}; - psPlane obsOffPt, refOffPt; - - // find location of 0,0 pixel in focal plane coords for this chip - psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); - - // find location of 0,0 pixel in focal plane coords for ref chip - psPlaneTransformApply (&refOffPt, refChip->toFPA, &offPixel); - - double obsAngle = atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); - double refAngle = atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x); - - bool badAstrom = false; - badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; - badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; - badAstrom |= fabs(obsAngle - refAngle) > angleTol; - - if (!badAstrom) continue; - - psFree (obsChip->toFPA); - psFree (obsChip->fromFPA); - - // XXX can we do a straight copy? I think so, but we may need to apply a scaling factor based - // on the toSky or toTPA values - obsChip->toFPA = psMemIncrRefCounter (refChip->toFPA); - obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA); - } - psFree (view); - return true; -} - -// XXX for this function to work, I need to: -// 1 *) define badAstrom (ref pixel too far from nominal?) (ref angle too far off?) -// 2) load the astrom fpa in the pmFPAfileIO stage -// 3) allow for this operation to be optional -# endif Index: unk/psastro/src/psastroOneChip.c =================================================================== --- /trunk/psastro/src/psastroOneChip.c (revision 19509) +++ (revision ) @@ -1,223 +1,0 @@ -# include "psastroInternal.h" - -# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE)\ - VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \ - if (!status) { \ - psAbort ("Failed to find %s in recipe", NAME); } - -bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) { - - bool status; - pmAstromStats *stats = NULL; - - // default value for match/fit : radius is in pixels - REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); - - // run the match/fit sequence NITER times - REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); - - // correct radius to FP units (physical pixel scale in microns per pixel) - REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); - RADIUS *= pixelScale; - - // select the desired chip order - REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32); - - // allowed limits for valid solutions - REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32); - REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32); - - // do we need to get a rough initial match? - REQUIRED_RECIPE_VALUE (bool gridSearch, "PSASTRO.GRID.SEARCH", Bool); - - // do we need to get a rough initial match? - REQUIRED_RECIPE_VALUE (int maxNstar, "PSASTRO.GRID.NSTAR.MAX", S32); - - if (gridSearch) { - // generate the bright subset of maxNstar entries (note: rawstars is sorted by S/N) - psArray *subset = psArrayAlloc (PS_MIN (maxNstar, rawstars->n)); - for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) { - subset->data[i] = psMemIncrRefCounter (rawstars->data[i]); - } - - // XXX set clump scale from recipe - psArray *gridStars = psastroRemoveClumps (subset, 150); - psFree (subset); - - psArray *refSubset = psastroRemoveClumps (refstars, 150); - - psLogMsg ("psastro", 3, "grid search using %ld raw vs %ld ref stars\n", gridStars->n, refSubset->n); - - // find initial offset / rotation / scale - pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refSubset, recipe); - if (gridStats == NULL) { - psLogMsg ("psastro", 3, "failed to find a grid match solution\n"); - psFree (gridStars); - psFree (refSubset); - return false; - } - psLogMsg ("psastro", 3, "basic grid search result - offset: %f,%f pixels, rotation: %f deg\n", gridStats->offset.x, gridStats->offset.y, DEG_RAD*gridStats->angle); - - // tweak the position by finding peak of matches stars - stats = pmAstromGridTweak (gridStars, refSubset, recipe, gridStats); - if (stats == NULL) { - psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n"); - psFree (gridStats); - psFree (gridStars); - psFree (refSubset); - return false; - } - psLogMsg ("psastro", 3, "tweak grid search result - offset: %f,%f pixels, rotation: %f deg\n", stats->offset.x, stats->offset.y, DEG_RAD*stats->angle); - - // adjust the chip.toFPA terms only - pmAstromGridApply (chip->toFPA, stats); - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - psFree (gridStats); - psFree (gridStars); - psFree (refSubset); - } - - psArray *match = NULL; - psStats *fitStats = NULL; - pmAstromFitResults *results = NULL; - - for (int iter = 0; iter < nIter; iter++) { - - char name[128]; - - sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter); - float radius = psMetadataLookupF32 (&status, recipe, name); - radius *= pixelScale; - if (!status || (radius == 0.0)) { - radius = RADIUS; - } - - - // use small radius to match stars - match = pmAstromRadiusMatchFP (rawstars, refstars, radius); - if (match == NULL) { - psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); - psFree (stats); - return false; - } - - // modify the order to correspond to the actual number of matched stars: - if ((match->n < 11) && (order >= 3)) order = 2; - if ((match->n < 7) && (order >= 2)) order = 1; - if ((match->n < 4) && (order >= 1)) order = 0; - if (order < 1) { - psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); - psFree (match); - psFree (stats); - return false; - } - - // create output toFPA; set masks appropriate to the Elixir DVO astrometry format - psFree (chip->toFPA); - chip->toFPA = psPlaneTransformAlloc (order, order); - for (int i = 0; i <= chip->toFPA->x->nX; i++) { - for (int j = 0; j <= chip->toFPA->x->nY; j++) { - if (i + j > order) { - chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET; - chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; - } - } - } - - // XXX allow statitic to be set by the user - // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA"); - fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER"); - - // improved fit for astrometric terms - results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); - if (!results) { - psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); - psFree (match); - psFree (stats); - psFree (fitStats); - return false; - } - - // determine fromFPA transformation and apply new transformation to raw & ref stars - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - - if (iter < nIter - 1) { - psFree (fitStats); - psFree (results); - psFree (match); - } - - // toSky converts from FPA & TPA units (microns) to sky units (radians) - float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; - // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; - float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; - int astNstar = results->yStats->clippedNvalues; - psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); - } - - - - // toSky converts from FPA & TPA units (microns) to sky units (radians) - float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; - - // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns) - // float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale; - float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale; - - // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns) - // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; - float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; - int astNstar = results->yStats->clippedNvalues; - - bool validSolution = true; - - // XXX should these result in errors or be handled another way? - psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar); - if (astError > maxError) { - psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); - validSolution = false; - } - if (astNstar < minNstar) { - psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); - validSolution = false; - } - - // DVO expects NASTRO = 0 if we fail to find a solution - psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR", PS_META_REPLACE, "astrometry error (pixels)", pixError); - psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "astrometry error (arcsec)", astError); - if (validSolution) { - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar)); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); - } else { - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); - } - psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars - - // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars - // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - - // XXX check if we correctly applied the new transformation: - if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRawstars (rawstars, fpa, chip); - psastroDumpMatchedStars ("match.dat", rawstars, refstars, match); - psastroDumpStars (refstars, "refstars.cal.dat"); - } - - if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotOneChipFit (rawstars, refstars, match, recipe); - } - - psFree (match); - psFree (stats); - psFree (results); - psFree (fitStats); - return validSolution; -} - -// psastroWriteStars ("raw.1.dat", rawstars); -// psastroWriteStars ("ref.1.dat", refstars); -// psastroWriteTransform (chip->toFPA); -