IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9639


Ignore:
Timestamp:
Oct 18, 2006, 4:56:40 PM (20 years ago)
Author:
magnier
Message:

added various error checks and cleaned up output stats

Location:
trunk/psModules/src/astrom
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/pmAstrometryObjects.c

    r9626 r9639  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-10-18 18:40:19 $
     10*  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-10-19 02:56:40 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3030#include "pmAstrometryObjects.h"
    3131
     32# define DEG_RAD 57.295779513082322
     33# define RAD_DEG  0.017453292519943
    3234
    3335#define PM_ASTROMETRYOBJECTS_DEBUG 1
    … …  
    293295    }
    294296    fclose (f);
     297    psFree (xfit);
     298    psFree (yfit);
    295299    # endif
    296300
    … …  
    299303    // XXX this is a somewhat silly place to write these updates...
    300304    if (updates) {
    301         psMetadataAddF32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", PS_MIN(x->n, y->n));
     305        psMetadataAddF32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", stats->clippedNvalues);
    302306        psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "", hypot(dX, dY));
    303307        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", hypot(dX, dY)/sqrt(PS_MIN(x->n, y->n)));
    304308        psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: must be defined somewhere
    305309    }
     310
     311    float maxError = psMetadataLookupF32 (&status, config, "PSASTRO.MAX.ERROR");
     312    float astError = hypot(dX, dY);
     313    int minNstar = psMetadataLookupS32 (&status, config, "PSASTRO.MIN.NSTAR");
     314    int astNstar = stats->clippedNvalues;
     315
     316    if (astError > maxError) {
     317        psError(PS_ERR_UNKNOWN, true, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
     318        return NULL;
     319    }
     320    if (stats->clippedNvalues < minNstar) {
     321        psError(PS_ERR_UNKNOWN, true, "solution uses too few stars: %d < %d", astNstar, minNstar);
     322        return NULL;
     323    }
     324    psLogMsg ("psModules.astrom", 3, "astrometry solution: error: %f, Nstars: %d", astError, astNstar);
    306325
    307326    psFree (x);
    … …  
    355374}
    356375
    357 
     376/******************************************************************************
     377pmAstromStatsFree(stats)
     378 ******************************************************************************/
     379static void pmAstromStatsFree(pmAstromStats *stats)
     380{
     381    if (stats == NULL)
     382        return;
     383    return;
     384}
     385
     386/******************************************************************************
     387pmAstromStatsAlloc()
     388 ******************************************************************************/
     389pmAstromStats *pmAstromStatsAlloc(void)
     390{
     391    pmAstromStats *stats = psAlloc (sizeof(pmAstromStats));
     392    psMemSetDeallocator (stats, (psFreeFunc)pmAstromStatsFree);
     393
     394    //    stats->center = {0, 0, 0, 0};
     395    //    stats->offset = {0, 0, 0, 0};
     396    stats->angle     = 0.0;
     397    stats->minMetric = 0.0;
     398    stats->minVar    = 0.0;
     399    stats->nMatch    = 0;
     400    stats->nTest     = 0;
     401    stats->nSigma    = 0;
     402
     403    return (stats);
     404}
    358405
    359406/******************************************************************************
    … …  
    488535delta-delta max.
    489536 ******************************************************************************/
    490 pmAstromStats pmAstromGridAngle(
     537pmAstromStats *pmAstromGridAngle(
    491538    const psArray *raw,
    492539    const psArray *ref,
    493540    const psMetadata *config)
    494541{
    495     // XXX: What to do if input parameters are bad?
    496     pmAstromStats badStat;
    497     PS_ASSERT_PTR_NON_NULL(raw, badStat);
    498     PS_ASSERT_PTR_NON_NULL(ref, badStat);
    499     PS_ASSERT_PTR_NON_NULL(config, badStat);
     542    PS_ASSERT_PTR_NON_NULL(raw, NULL);
     543    PS_ASSERT_PTR_NON_NULL(ref, NULL);
     544    PS_ASSERT_PTR_NON_NULL(config, NULL);
    500545
    501546    bool status;
    … …  
    506551
    507552    const pmAstromObj *ob1, *ob2; // short-cut pointers to the objects
    508     pmAstromStats stats;    // output match statistics
     553
     554    pmAstromStats *stats = pmAstromStatsAlloc();    // output match statistics
    509555
    510556    // max allowed offset in either X or Y directions
    … …  
    534580
    535581    // short-cut names for grid images
    536     psS32 **NP = gridNP->data.S32;
     582    psU32 **NP = gridNP->data.U32;
    537583    psF32 **DX = gridDX->data.F32;
    538584    psF32 **DY = gridDY->data.F32;
    … …  
    570616
    571617        // find the max pixel
    572         psStats *imStats = psStatsAlloc (PS_STAT_MAX);
     618        psStats *imStats = psStatsAlloc (PS_STAT_MAX | PS_STAT_MAX | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    573619        imStats = psImageStats (imStats, gridNP, NULL, 0);
    574620
    … …  
    601647
    602648        // convert the bin to delta-delta
    603         stats.offset.x  = DX[minY][minX] / NP[minY][minX];
    604         stats.offset.y  = DY[minY][minX] / NP[minY][minX];
    605         stats.minMetric = minMetric;
    606         stats.minVar    = minVar;
    607         stats.nMatch    = NP[minY][minX];
     649        stats->offset.x  = DX[minY][minX] / NP[minY][minX];
     650        stats->offset.y  = DY[minY][minX] / NP[minY][minX];
     651        stats->minMetric = minMetric;
     652        stats->minVar    = minVar;
     653        stats->nMatch    = NP[minY][minX];
    608654
    609655        psFree (imStats);
    610656        // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better
    611657    }
     658
     659    // sort the NP values and choose
     660    psVector *listNP = psVectorAlloc (nPix*nPix, PS_TYPE_U32);
     661    listNP->n = nPix*nPix;
     662    int n = 0;
     663    for (int i = 0; i < nPix; i++) {
     664        for (int j = 0; j < nPix; j++) {
     665            listNP->data.U32[n] = gridNP->data.U32[j][i];
     666            n++;
     667        }
     668    }
     669    psVector *sort = psVectorSort (NULL, listNP);
     670    stats->nTest = sort->data.U32[sort->n - 5];
     671    stats->nSigma = (stats->nMatch - stats->nTest) / sqrt(stats->nTest);
     672
     673    psFree (sort);
     674    psFree (listNP);
    612675    psFree (gridNP);
    613676    psFree (gridDX);
    … …  
    623686 ******************************************************************************/
    624687
    625 pmAstromStats pmAstromGridMatch(const psArray *raw,
    626                                 const psArray *ref,
    627                                 const psMetadata *config)
    628 {
    629     // XXX: What to do if input parameters are bad?
    630     pmAstromStats badStat;
    631     PS_ASSERT_PTR_NON_NULL(raw, badStat);
    632     PS_ASSERT_PTR_NON_NULL(ref, badStat);
    633     PS_ASSERT_PTR_NON_NULL(config, badStat);
     688pmAstromStats *pmAstromGridMatch(const psArray *raw,
     689                                 const psArray *ref,
     690                                 const psMetadata *config)
     691{
     692    PS_ASSERT_PTR_NON_NULL(raw, NULL);
     693    PS_ASSERT_PTR_NON_NULL(ref, NULL);
     694    PS_ASSERT_PTR_NON_NULL(config, NULL);
    634695
    635696    bool status;
    … …  
    638699    psArray *rot;
    639700
    640     pmAstromStats minStat;
    641     pmAstromStats newStat = {{0, 0, 0, 0}, {0, 0, 0, 0}, 0, 0, 0, 0};
     701    pmAstromStats *minStat = pmAstromStatsAlloc ();
     702    pmAstromStats *newStat = NULL;
     703
    642704    psPlane center;
    643705
    … …  
    655717    center.y = 0.5*(yMin + yMax);
    656718
    657     double minAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE");
    658     double maxAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE");
    659     double delAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE");
    660 
    661     minStat.minMetric = 1e10;
     719    double minAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE");
     720    double maxAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE");
     721    double delAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE");
     722    double minSigma = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.SIGMA");
     723
     724    minStat->minMetric = 1e10;
    662725    for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
    663726        rot = pmAstromRotateObj (raw, center, angle);
    664727        newStat = pmAstromGridAngle (rot, ref, config);
    665         psTrace("psModule.astrom.grid.match", 5, "try %f : %f %f (%d pts, %f var, %f met (%f log))",
    666                 angle, newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric,
    667                 log10(newStat.minMetric));
    668 
    669         newStat.angle  = angle;
    670         newStat.center = center;
    671         if (newStat.minMetric < minStat.minMetric) {
    672             minStat = newStat;
    673             psTrace("psModule.astrom.grid.match", 5, "use %f : %f %f (%d pts, %f var, %f met (%f log))",
    674                     angle, minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, minStat.
    675                     minMetric, log10(minStat.minMetric));
    676         }
     728        newStat->angle  = angle;
     729        newStat->center = center;
     730
     731        if (newStat->minMetric < minStat->minMetric) {
     732            *minStat = *newStat;
     733            psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *",
     734                      minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
     735        } else {
     736            psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)",
     737                      newStat->offset.x, newStat->offset.y, DEG_RAD*newStat->angle, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric));
     738
     739        }
     740        psFree (newStat);
    677741        psFree (rot);
    678742    }
    679     psTrace("psModule.astrom.grid.match", 4, "best: %f %f (%d pts, %f var, %f met)",
    680             minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, log10(minStat.minMetric));
     743    psLogMsg ("psModule.astrom.grid.match", 4, "grid best - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)",
     744              minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
     745
     746    // I need to decide if a solution is likely to be a good solution or just a mis-match
     747    // one posibility: how significant is the peak relative to the 4th or 5th most significant pixel?
     748
     749    if (minStat->nSigma < minSigma) {
     750        psError(PS_ERR_UNKNOWN, true, "Failed to find a valid match (%f sigma for best)", minStat->nSigma);
     751        return NULL;
     752    }
    681753    return (minStat);
    682754}
    … …  
    685757pmAstromGridTweak(*raw, *ref, *recipe, stats): improve match for two star lists.
    686758 ******************************************************************************/
    687 pmAstromStats pmAstromGridTweak(
     759pmAstromStats *pmAstromGridTweak(
    688760    psArray *raw,
    689761    psArray *ref,
    690762    psMetadata *recipe,
    691     pmAstromStats stats)
     763    pmAstromStats *stats)
    692764{
    693765    bool status;
    … …  
    697769    int nBin, xBin, yBin;
    698770
    699     rot = pmAstromRotateObj (raw, stats.center, stats.angle);
     771    rot = pmAstromRotateObj (raw, stats->center, stats->angle);
    700772
    701773    // sampling scale of the grid
    … …  
    717789        for (int j = 0; j < ref->n; j++) {
    718790            ob2 = (pmAstromObj *)ref->data[j];
    719             dX = ob1->FP->x - ob2->FP->x - stats.offset.x;
    720             dY = ob1->FP->y - ob2->FP->y - stats.offset.y;
     791            dX = ob1->FP->x - ob2->FP->x - stats->offset.x;
     792            dY = ob1->FP->y - ob2->FP->y - stats->offset.y;
    721793
    722794            xBin = (dX + tweakRange) / tweakScale;
    … …  
    760832
    761833    // adjust offset by peak center
    762     pmAstromStats tweak = stats;
    763     tweak.offset.x += xPeak;
    764     tweak.offset.y += yPeak;
     834    pmAstromStats *tweak = pmAstromStatsAlloc();
     835    *tweak = *stats;
     836    tweak->offset.x += xPeak;
     837    tweak->offset.y += yPeak;
    765838
    766839    psFree (rot);
    … …  
    777850psPlaneTransform *pmAstromGridApply(
    778851    psPlaneTransform *map,
    779     pmAstromStats stat)
     852    pmAstromStats *stat)
    780853{
    781854    PS_ASSERT_PTR_NON_NULL(map, NULL);
    … …  
    783856    PS_ASSERT_POLY_NON_NULL(map->y, NULL);
    784857
    785     double cs = cos (stat.angle);
    786     double sn = sin (stat.angle);
    787 
    788     double dx = (map->x->coeff[0][0] - stat.center.x);
    789     double dy = (map->y->coeff[0][0] - stat.center.y);
     858    double cs = cos (stat->angle);
     859    double sn = sin (stat->angle);
     860
     861    double dx = (map->x->coeff[0][0] - stat->center.x);
     862    double dy = (map->y->coeff[0][0] - stat->center.y);
    790863
    791864    // new offset
    792     map->x->coeff[0][0] =  cs*dx + sn*dy - stat.offset.x + stat.center.x;
    793     map->y->coeff[0][0] = -sn*dx + cs*dy - stat.offset.y + stat.center.y;
     865    map->x->coeff[0][0] =  cs*dx + sn*dy - stat->offset.x + stat->center.x;
     866    map->y->coeff[0][0] = -sn*dx + cs*dy - stat->offset.y + stat->center.y;
    794867
    795868    // original rotation matrix

  • trunk/psModules/src/astrom/pmAstrometryObjects.h

    r9386 r9639  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-10-07 03:54:08 $
     10*  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-10-19 02:56:40 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    7474    double  minVar;                     ///<
    7575    int     nMatch;                     ///<
     76    int     nTest;                      ///<
     77    double  nSigma;                     ///<
    7678}
    7779pmAstromStats;
    78 
    79 
    8080
    8181/*
    … …  
    108108
    109109
     110pmAstromStats *pmAstromStatsAlloc(void);
    110111
    111112/*
    … …  
    153154 *
    154155 */
    155 pmAstromStats pmAstromGridMatch(
     156pmAstromStats *pmAstromGridMatch(
    156157    const psArray *st1,
    157158    const psArray *st2,
    … …  
    162163pmAstromGridTweak(*raw, *ref, *recipe, stats): improve match for two star lists.
    163164 ******************************************************************************/
    164 pmAstromStats pmAstromGridTweak(
     165pmAstromStats *pmAstromGridTweak(
    165166    psArray *raw,
    166167    psArray *ref,
    167168    psMetadata *recipe,
    168     pmAstromStats stats);
     169    pmAstromStats *stats);
    169170
    170171/*
    … …  
    183184psPlaneTransform *pmAstromGridApply(
    184185    psPlaneTransform *map,
    185     pmAstromStats stat
     186    pmAstromStats *stat
    186187);
    187188
    … …  
    197198 */
    198199/* in pmAstromGrid.c */
    199 pmAstromStats pmAstromGridAngle(
     200pmAstromStats *pmAstromGridAngle(
    200201    const psArray *st1,
    201202    const psArray *st2,
Note: See TracChangeset for help on using the changeset viewer.