Index: trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 9626) +++ trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 9639) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-10-18 18:40:19 $ +* @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-10-19 02:56:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -30,4 +30,6 @@ #include "pmAstrometryObjects.h" +# define DEG_RAD 57.295779513082322 +# define RAD_DEG 0.017453292519943 #define PM_ASTROMETRYOBJECTS_DEBUG 1 @@ -293,4 +295,6 @@ } fclose (f); + psFree (xfit); + psFree (yfit); # endif @@ -299,9 +303,24 @@ // XXX this is a somewhat silly place to write these updates... if (updates) { - psMetadataAddF32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", PS_MIN(x->n, y->n)); + 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); psFree (x); @@ -355,5 +374,33 @@ } - +/****************************************************************************** +pmAstromStatsFree(stats) + ******************************************************************************/ +static void pmAstromStatsFree(pmAstromStats *stats) +{ + if (stats == NULL) + return; + return; +} + +/****************************************************************************** +pmAstromStatsAlloc() + ******************************************************************************/ +pmAstromStats *pmAstromStatsAlloc(void) +{ + pmAstromStats *stats = psAlloc (sizeof(pmAstromStats)); + psMemSetDeallocator (stats, (psFreeFunc)pmAstromStatsFree); + + // stats->center = {0, 0, 0, 0}; + // stats->offset = {0, 0, 0, 0}; + stats->angle = 0.0; + stats->minMetric = 0.0; + stats->minVar = 0.0; + stats->nMatch = 0; + stats->nTest = 0; + stats->nSigma = 0; + + return (stats); +} /****************************************************************************** @@ -488,14 +535,12 @@ delta-delta max. ******************************************************************************/ -pmAstromStats pmAstromGridAngle( +pmAstromStats *pmAstromGridAngle( const psArray *raw, const psArray *ref, const psMetadata *config) { - // XXX: What to do if input parameters are bad? - pmAstromStats badStat; - PS_ASSERT_PTR_NON_NULL(raw, badStat); - PS_ASSERT_PTR_NON_NULL(ref, badStat); - PS_ASSERT_PTR_NON_NULL(config, badStat); + PS_ASSERT_PTR_NON_NULL(raw, NULL); + PS_ASSERT_PTR_NON_NULL(ref, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); bool status; @@ -506,5 +551,6 @@ const pmAstromObj *ob1, *ob2; // short-cut pointers to the objects - pmAstromStats stats; // output match statistics + + pmAstromStats *stats = pmAstromStatsAlloc(); // output match statistics // max allowed offset in either X or Y directions @@ -534,5 +580,5 @@ // short-cut names for grid images - psS32 **NP = gridNP->data.S32; + psU32 **NP = gridNP->data.U32; psF32 **DX = gridDX->data.F32; psF32 **DY = gridDY->data.F32; @@ -570,5 +616,5 @@ // find the max pixel - psStats *imStats = psStatsAlloc (PS_STAT_MAX); + psStats *imStats = psStatsAlloc (PS_STAT_MAX | PS_STAT_MAX | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); imStats = psImageStats (imStats, gridNP, NULL, 0); @@ -601,13 +647,30 @@ // convert the bin to delta-delta - stats.offset.x = DX[minY][minX] / NP[minY][minX]; - stats.offset.y = DY[minY][minX] / NP[minY][minX]; - stats.minMetric = minMetric; - stats.minVar = minVar; - stats.nMatch = NP[minY][minX]; + stats->offset.x = DX[minY][minX] / NP[minY][minX]; + stats->offset.y = DY[minY][minX] / NP[minY][minX]; + stats->minMetric = minMetric; + stats->minVar = minVar; + stats->nMatch = NP[minY][minX]; psFree (imStats); // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better } + + // sort the NP values and choose + psVector *listNP = psVectorAlloc (nPix*nPix, PS_TYPE_U32); + listNP->n = nPix*nPix; + int n = 0; + for (int i = 0; i < nPix; i++) { + for (int j = 0; j < nPix; j++) { + listNP->data.U32[n] = gridNP->data.U32[j][i]; + n++; + } + } + psVector *sort = psVectorSort (NULL, listNP); + stats->nTest = sort->data.U32[sort->n - 5]; + stats->nSigma = (stats->nMatch - stats->nTest) / sqrt(stats->nTest); + + psFree (sort); + psFree (listNP); psFree (gridNP); psFree (gridDX); @@ -623,13 +686,11 @@ ******************************************************************************/ -pmAstromStats pmAstromGridMatch(const psArray *raw, - const psArray *ref, - const psMetadata *config) -{ - // XXX: What to do if input parameters are bad? - pmAstromStats badStat; - PS_ASSERT_PTR_NON_NULL(raw, badStat); - PS_ASSERT_PTR_NON_NULL(ref, badStat); - PS_ASSERT_PTR_NON_NULL(config, badStat); +pmAstromStats *pmAstromGridMatch(const psArray *raw, + const psArray *ref, + const psMetadata *config) +{ + PS_ASSERT_PTR_NON_NULL(raw, NULL); + PS_ASSERT_PTR_NON_NULL(ref, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); bool status; @@ -638,6 +699,7 @@ psArray *rot; - pmAstromStats minStat; - pmAstromStats newStat = {{0, 0, 0, 0}, {0, 0, 0, 0}, 0, 0, 0, 0}; + pmAstromStats *minStat = pmAstromStatsAlloc (); + pmAstromStats *newStat = NULL; + psPlane center; @@ -655,28 +717,38 @@ center.y = 0.5*(yMin + yMax); - double minAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE"); - double maxAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE"); - double delAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE"); - - minStat.minMetric = 1e10; + double minAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE"); + double maxAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE"); + double delAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE"); + double minSigma = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.SIGMA"); + + minStat->minMetric = 1e10; for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { rot = pmAstromRotateObj (raw, center, angle); newStat = pmAstromGridAngle (rot, ref, config); - psTrace("psModule.astrom.grid.match", 5, "try %f : %f %f (%d pts, %f var, %f met (%f log))", - angle, newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric, - log10(newStat.minMetric)); - - newStat.angle = angle; - newStat.center = center; - if (newStat.minMetric < minStat.minMetric) { - minStat = newStat; - psTrace("psModule.astrom.grid.match", 5, "use %f : %f %f (%d pts, %f var, %f met (%f log))", - angle, minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, minStat. - minMetric, log10(minStat.minMetric)); - } + newStat->angle = angle; + newStat->center = center; + + if (newStat->minMetric < minStat->minMetric) { + *minStat = *newStat; + 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) *", + minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); + } else { + 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)", + newStat->offset.x, newStat->offset.y, DEG_RAD*newStat->angle, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric)); + + } + psFree (newStat); psFree (rot); } - psTrace("psModule.astrom.grid.match", 4, "best: %f %f (%d pts, %f var, %f met)", - minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, log10(minStat.minMetric)); + 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)", + minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); + + // I need to decide if a solution is likely to be a good solution or just a mis-match + // one posibility: how significant is the peak relative to the 4th or 5th most significant pixel? + + if (minStat->nSigma < minSigma) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a valid match (%f sigma for best)", minStat->nSigma); + return NULL; + } return (minStat); } @@ -685,9 +757,9 @@ pmAstromGridTweak(*raw, *ref, *recipe, stats): improve match for two star lists. ******************************************************************************/ -pmAstromStats pmAstromGridTweak( +pmAstromStats *pmAstromGridTweak( psArray *raw, psArray *ref, psMetadata *recipe, - pmAstromStats stats) + pmAstromStats *stats) { bool status; @@ -697,5 +769,5 @@ int nBin, xBin, yBin; - rot = pmAstromRotateObj (raw, stats.center, stats.angle); + rot = pmAstromRotateObj (raw, stats->center, stats->angle); // sampling scale of the grid @@ -717,6 +789,6 @@ for (int j = 0; j < ref->n; j++) { ob2 = (pmAstromObj *)ref->data[j]; - dX = ob1->FP->x - ob2->FP->x - stats.offset.x; - dY = ob1->FP->y - ob2->FP->y - stats.offset.y; + dX = ob1->FP->x - ob2->FP->x - stats->offset.x; + dY = ob1->FP->y - ob2->FP->y - stats->offset.y; xBin = (dX + tweakRange) / tweakScale; @@ -760,7 +832,8 @@ // adjust offset by peak center - pmAstromStats tweak = stats; - tweak.offset.x += xPeak; - tweak.offset.y += yPeak; + pmAstromStats *tweak = pmAstromStatsAlloc(); + *tweak = *stats; + tweak->offset.x += xPeak; + tweak->offset.y += yPeak; psFree (rot); @@ -777,5 +850,5 @@ psPlaneTransform *pmAstromGridApply( psPlaneTransform *map, - pmAstromStats stat) + pmAstromStats *stat) { PS_ASSERT_PTR_NON_NULL(map, NULL); @@ -783,13 +856,13 @@ PS_ASSERT_POLY_NON_NULL(map->y, NULL); - double cs = cos (stat.angle); - double sn = sin (stat.angle); - - double dx = (map->x->coeff[0][0] - stat.center.x); - double dy = (map->y->coeff[0][0] - stat.center.y); + double cs = cos (stat->angle); + double sn = sin (stat->angle); + + double dx = (map->x->coeff[0][0] - stat->center.x); + double dy = (map->y->coeff[0][0] - stat->center.y); // new offset - map->x->coeff[0][0] = cs*dx + sn*dy - stat.offset.x + stat.center.x; - map->y->coeff[0][0] = -sn*dx + cs*dy - stat.offset.y + stat.center.y; + map->x->coeff[0][0] = cs*dx + sn*dy - stat->offset.x + stat->center.x; + map->y->coeff[0][0] = -sn*dx + cs*dy - stat->offset.y + stat->center.y; // original rotation matrix