Index: trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 10780) +++ trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 10829) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-17 09:45:20 $ +* @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-25 01:49:49 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -121,24 +121,18 @@ /************************************************************************************************************/ - +// macro to generate code for radius match function based on desired member +// radius is in units of matching member (eg, pixels for chip, microns for FP, etc) #define MAKE_ASTROM_RADIUS(FUNC, MEMBER) \ psArray *FUNC( \ const psArray *st1, \ const psArray *st2, \ - const psMetadata *config) \ + double RADIUS) \ { \ PS_ASSERT_PTR_NON_NULL(st1, NULL); \ PS_ASSERT_PTR_NON_NULL(st2, NULL); \ - PS_ASSERT_PTR_NON_NULL(config, NULL); \ \ assert(st1->n == 0 || pmIsAstromObj(st1->data[0])); \ assert(st2->n == 0 || pmIsAstromObj(st2->data[0])); \ \ - bool status = false; \ - double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS"); \ - if (!status) { \ - psError(PS_ERR_IO, false, "Failed to lookup matching radius"); \ - return NULL; \ - } \ /* sort both lists by X coord; st1 first */ \ psVector *x1 = psVectorAlloc(st1->n, PS_TYPE_F64); \ @@ -186,18 +180,18 @@ */ MAKE_ASTROM_RADIUS(pmAstromRadiusMatch, FP) +MAKE_ASTROM_RADIUS(pmAstromRadiusMatchFP, FP) MAKE_ASTROM_RADIUS(pmAstromRadiusMatchTP, TP) MAKE_ASTROM_RADIUS(pmAstromRadiusMatchChip, chip) /****************************************************************************** -pmAstromMatchFit(map, raw, ref, match, config): take two matched star lists +pmAstromMatchFit(map, raw, ref, match, stats): take two matched star lists and fit a psPlaneTransform between them ******************************************************************************/ -bool pmAstromMatchFit( +pmAstromFitResults *pmAstromMatchFit( psPlaneTransform *map, psArray *raw, psArray *ref, psArray *match, - psMetadata *config, - psMetadata *updates) + psStats *stats) { PS_ASSERT_PTR_NON_NULL(map, NULL); @@ -205,12 +199,8 @@ PS_ASSERT_PTR_NON_NULL(ref, NULL); PS_ASSERT_PTR_NON_NULL(match, NULL); - PS_ASSERT_PTR_NON_NULL(config, NULL); - - bool status; - pmAstromObj *rawStar, *refStar; - pmAstromMatch *pair; - - // XXX EAM : clip fit seems to only work for F32! - // XXX EAM : clip fit fails to handle NULL error + PS_ASSERT_PTR_NON_NULL(stats, NULL); + + // reassign values for clip fit + // XXX set wt based on mag error? psVector *X = psVectorAlloc (match->n, PS_TYPE_F32); psVector *Y = psVectorAlloc (match->n, PS_TYPE_F32); @@ -220,8 +210,7 @@ // take the matched stars, first fit for (int i = 0; i < match->n; i++) { - - pair = match->data[i]; - rawStar = raw->data[pair->raw]; - refStar = ref->data[pair->ref]; + pmAstromMatch *pair = match->data[i]; + pmAstromObj *rawStar = raw->data[pair->raw]; + pmAstromObj *refStar = ref->data[pair->ref]; X->data.F32[i] = rawStar->chip->x; @@ -238,65 +227,23 @@ psVectorInit (mask, 0); - double dX, dY; - - psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - stats->clipSigma = psMetadataLookupF32 (&status, config, "PSASTRO.CHIP.NSIGMA"); - stats->clipIter = 1; - int clipIter = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NITER"); - // we run 'clipIter' cycles clipping in each of x and y + pmAstromFitResults *results = pmAstromFitResultsAlloc(); + results->xStats = psStatsAlloc (PS_STAT_NONE); + results->yStats = psStatsAlloc (PS_STAT_NONE); + *results->xStats = *stats; + *results->yStats = *stats; + results->xStats->clipIter = 1; + results->yStats->clipIter = 1; // fit chip-to-FPA transformation - // XXX should force x->n and y->n to be the same - for (int i = 0; i < clipIter; i++) { - psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, x, wt, X, Y); - psLogMsg ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, x->n); - dX = stats->clippedStdev; - - psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, y, wt, X, Y); - psLogMsg ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, y->n); - dY = stats->clippedStdev; - } - - // XXXX test dump output code - # if (0) - psVector *xfit = psPolynomial2DEvalVector(map->x, X, Y); - psVector *yfit = psPolynomial2DEvalVector(map->y, X, Y); - FILE *f = fopen ("match.dat", "w"); - for (int i = 0; i < x->n; i++) { - fprintf (f, "%3d %1d %1d %f %f %f %f %f %f\n", - i, mask->data.U8[i], mask->data.U8[i], - X->data.F32[i], Y->data.F32[i], - xfit->data.F32[i], yfit->data.F32[i], - x->data.F32[i], y->data.F32[i]); - } - fclose (f); - psFree (xfit); - psFree (yfit); - # endif - - // XXX this is somewhat ad-hoc: we need to be sure these values are in arcsec - // and we need to set NASTRO = 0 if we have a failure, etc - // XXX this is a somewhat silly place to write these updates... - if (updates) { - psMetadataAddF32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", stats->clippedNvalues); - psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "", hypot(dX, dY)); - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", hypot(dX, dY)/sqrt(PS_MIN(x->n, y->n))); - psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "", 2000.0); // XXX this is bogus: must be defined somewhere - } - - float maxError = psMetadataLookupF32 (&status, config, "PSASTRO.MAX.ERROR"); - float astError = hypot(dX, dY); - int minNstar = psMetadataLookupS32 (&status, config, "PSASTRO.MIN.NSTAR"); - int astNstar = stats->clippedNvalues; - - if (astError > maxError) { - psError(PS_ERR_UNKNOWN, true, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); - return NULL; - } - if (stats->clippedNvalues < minNstar) { - psError(PS_ERR_UNKNOWN, true, "solution uses too few stars: %d < %d", astNstar, minNstar); - return NULL; - } - psLogMsg ("psModules.astrom", 3, "astrometry solution: error: %f, Nstars: %d", astError, astNstar); + // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each + for (int i = 0; i < stats->clipIter; i++) { + psVectorClipFitPolynomial2D (map->x, results->xStats, mask, 0xff, x, wt, X, Y); + psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, x->n); + + psVectorClipFitPolynomial2D (map->y, results->yStats, mask, 0xff, y, wt, X, Y); + psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, y->n); + } + results->xStats->clipIter = stats->clipIter; + results->yStats->clipIter = stats->clipIter; psFree (x); @@ -305,8 +252,7 @@ psFree (Y); psFree (wt); - psFree (stats); psFree (mask); - return (map); + return (results); } @@ -377,4 +323,32 @@ return (stats); +} + +/****************************************************************************** +pmAstromFitResultsFree(stats) + ******************************************************************************/ +static void pmAstromFitResultsFree(pmAstromFitResults *results) +{ + if (results == NULL) + return; + psFree (results->xStats); + psFree (results->yStats); + return; +} + +/****************************************************************************** +pmAstromFitResultsAlloc() + ******************************************************************************/ +pmAstromFitResults *pmAstromFitResultsAlloc(void) +{ + pmAstromFitResults *results = psAlloc (sizeof(pmAstromFitResults)); + psMemSetDeallocator (results, (psFreeFunc)pmAstromFitResultsFree); + + results->xStats = NULL; + results->yStats = NULL; + results->nMatch = 0; + results->nSigma = 0; + + return (results); }