Changeset 10829 for trunk/psModules/src/astrom/pmAstrometryObjects.c

Timestamp:

Dec 24, 2006, 3:49:50 PM (20 years ago)

Author:

magnier

Message:

some API work to reduce dependencies on the recipes, support for mosaic astrometry

File:

• trunk/psModules/src/astrom/pmAstrometryObjects.c (modified) (8 diffs)

trunk/psModules/src/astrom/pmAstrometryObjects.c

@@ -8 +8 @@
 *  @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 +121 @@
 
 /************************************************************************************************************/
-
+// 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 +180 @@
  */
 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 +199 @@
     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 +210 @@
     // 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 +227 @@
     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 +252 @@
     psFree (Y);
     psFree (wt);
-    psFree (stats);
     psFree (mask);
 
-    return (map);
+    return (results);
 }
 
@@ -377 +323 @@
 
     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);
 }