Index: /trunk/psModules/src/astrom/Makefile.am =================================================================== --- /trunk/psModules/src/astrom/Makefile.am (revision 5673) +++ /trunk/psModules/src/astrom/Makefile.am (revision 5674) @@ -3,7 +3,10 @@ libpsmoduleastrom_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) libpsmoduleastrom_la_LDFLAGS = -release $(PACKAGE_VERSION) -libpsmoduleastrom_la_SOURCES = pmAstrometry.c +libpsmoduleastrom_la_SOURCES = \ + pmAstrometry.c \ + pmAstrometryObjects.c psmoduleincludedir = $(includedir) psmoduleinclude_HEADERS = \ - pmAstrometry.h + pmAstrometry.h \ + pmAstrometryObjects.h Index: /trunk/psModules/src/astrom/pmAstrometry.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometry.c (revision 5673) +++ /trunk/psModules/src/astrom/pmAstrometry.c (revision 5674) @@ -10,7 +10,9 @@ * XXX: We should review the extent of the warning messages on these functions * when the transformations are not successful. -* -* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-27 19:28:56 $ +* +* XXX: Should we implement non-linear cell->chip transforms? +* +* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-05 20:49:30 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -23,5 +25,4 @@ #include #include "pslib.h" -#include "psDB.h" #include "pmAstrometry.h" @@ -34,7 +35,8 @@ XXX: Macro this for speed. *****************************************************************************/ -static psS32 checkValidImageCoords(double x, - double y, - psImage* tmpImage) +static psS32 checkValidImageCoords( + double x, + double y, + psImage* tmpImage) { PS_ASSERT_IMAGE_NON_NULL(tmpImage, 0); @@ -51,4 +53,5 @@ /****************************************************************************** *****************************************************************************/ +#define PS_FREE_HIERARCHY 1 static void readoutFree(pmReadout *readout) { @@ -77,5 +80,4 @@ } -#define PS_FREE_HIERARCHY 1 static void cellFree(pmCell *cell) { @@ -255,5 +257,5 @@ } -bool cellCheckParents(pmCell *cell) +static psBool cellCheckParents(pmCell *cell) { if (cell == NULL) { @@ -273,5 +275,5 @@ } -bool chipCheckParents(pmChip *chip) +static psBool chipCheckParents(pmChip *chip) { if (chip == NULL) { @@ -293,5 +295,5 @@ } -bool pmFPACheckParents(pmFPA *fpa) +psBool pmFPACheckParents(pmFPA *fpa) { if (fpa == NULL) { @@ -617,8 +619,5 @@ /***************************************************************************** -XXX: What about units for the (x,y) coords? - -XXX: This has not been tested. - *****************************************************************************/ + *****************************************************************************/ psSphere* pmCoordTPToSky( psSphere* outSphere, @@ -635,8 +634,6 @@ return(rc); } + /***************************************************************************** -XXX: What about units for the (x,y) coords? - -XXX: This has not been tested. *****************************************************************************/ psPlane* pmCoordSkyToTP( @@ -655,7 +652,6 @@ } - -// XXX: Should we check return codes of intermediate transforms? -// XXX: This has not been tested. +/***************************************************************************** + *****************************************************************************/ psSphere* pmCoordCellToSky( psSphere* skyCoord, @@ -700,4 +696,6 @@ } +/***************************************************************************** + *****************************************************************************/ psPlane* pmCoordSkyToCell( psPlane* cellCoord, @@ -749,6 +747,6 @@ } -// XXX: This has not been tested. -// XXX: Verify the plane coords to sphere coord code. +/***************************************************************************** + *****************************************************************************/ psSphere* pmCoordCellToSkyQuick( psSphere* outSphere, @@ -778,9 +776,4 @@ /***************************************************************************** -XXX: What about units for the (x,y) coords? - -XXX: This has not been tested. - -XXX: Verify the plane coords to sphere coord code. *****************************************************************************/ psPlane* pmCoordSkyToCellQuick( Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 5674) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 5674) @@ -0,0 +1,625 @@ +/** @file pmAstrometryObjects.c +* +* @brief This file defines the basic types for matching objects +* based on their astrometry. +* +* @ingroup AstroImage +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-05 20:49:30 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include "pslib.h" +#include "pmAstrometry.h" +#include "pmAstrometryObjects.h" + + +#define PM_ASTROMETRYOBJECTS_DEBUG 1 +/****************************************************************************** +pmAstromObjSortByFPX(**a, **b): sort by mag (descending) + +Is this a private routine? +Should we do the early asserts? + ******************************************************************************/ + +int pmAstromObjSortByFPX( + const void **a, + const void **b) +{ + if (PM_ASTROMETRYOBJECTS_DEBUG) { + PS_ASSERT_PTR_NON_NULL(a, 0); + PS_ASSERT_PTR_NON_NULL(*a, 0); + PS_ASSERT_PTR_NON_NULL(b, 0); + PS_ASSERT_PTR_NON_NULL(*b, 0); + } + + pmAstromObj *A = *(pmAstromObj **)a; + pmAstromObj *B = *(pmAstromObj **)b; + + psF32 diff = A->FP.x - B->FP.x; + if (diff > FLT_EPSILON) { + return (+1); + } + + if (diff < FLT_EPSILON) { + return (-1); + } + return (0); +} + + + +/****************************************************************************** +pmAstromObjSortByMag(**a, **b): sort by mag (descending) + +Is this a private routine? +Should we do the early asserts? + ******************************************************************************/ +int pmAstromObjSortByMag( + const void **a, + const void **b) +{ + if (PM_ASTROMETRYOBJECTS_DEBUG) { + PS_ASSERT_PTR_NON_NULL(a, 0); + PS_ASSERT_PTR_NON_NULL(*a, 0); + PS_ASSERT_PTR_NON_NULL(b, 0); + PS_ASSERT_PTR_NON_NULL(*b, 0); + } + + pmAstromObj *A = *(pmAstromObj **)a; + pmAstromObj *B = *(pmAstromObj **)b; + + psF32 diff = A->Mag - B->Mag; + if (diff > FLT_EPSILON) { + return (-1); + } + + if (diff < FLT_EPSILON) { + return (+1); + } + + return (0); +} + + +/****************************************************************************** +pmAstromRadiusMatch(st1, st2, config) + ******************************************************************************/ +psArray *pmAstromRadiusMatch( + psArray *st1, + psArray *st2, + psMetadata *config) +{ + PS_ASSERT_PTR_NON_NULL(st1, NULL); + PS_ASSERT_PTR_NON_NULL(st2, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + bool status; + int jStart; + double dX, dY, dR; + + double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS"); + double RADIUS_SQR = PS_SQR (RADIUS); + + // sort both lists by Focal Plane X coord + st1 = psArraySort (st1, pmAstromObjSortByFPX); + st2 = psArraySort (st2, pmAstromObjSortByFPX); + + psArray *matches = psArrayAlloc (100); + matches->n = 0; + + int i = 0; + int j = 0; + while ((i < st1->n) && (j < st2->n)) { + + dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x; + if (dX < -RADIUS) { + i++; + continue; + } + if (dX > +RADIUS) { + j++; + continue; + } + + jStart = j; + while ((dX > -RADIUS) && (j < st2->n)) { + + dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x; + dY = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y; + dR = dX*dX + dY*dY; + + if (dR > RADIUS_SQR) { + j++; + continue; + } + + // got a match; add to output list + pmAstromMatch *match = pmAstromMatchAlloc (i, j); + psArrayAdd (matches, 100, match); + + j++; + } + j = jStart; + i++; + } + return (matches); +} + + + +/****************************************************************************** +pmAstromMatchFit(map, raw, ref, match, config): take two matched star lists +and fit a psPlaneTransform between them + ******************************************************************************/ +psPlaneTransform *pmAstromMatchFit( + psPlaneTransform *map, + psArray *raw, + psArray *ref, + psArray *match, + psMetadata *config) +{ + PS_ASSERT_PTR_NON_NULL(raw, NULL); + 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; + + if (map == NULL) { + int nX = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NX"); + if (!status) { + nX = 1; + } + int nY = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NY"); + if (!status) { + nY = 1; + } + map = psPlaneTransformAlloc (nX, nY); + // XXX EAM : not clear that we are allowed to use orders > 1 + } + + psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + stats->clipSigma = psMetadataLookupF32 (&status, config, "PSASTRO.CHIP.NSIGMA"); + stats->clipIter = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NITER"); + + // XXX EAM : clip fit seems to only work for F32! + // XXX EAM : clip fit fails to handle NULL error + psVector *X = psVectorAlloc (match->n, PS_TYPE_F32); + psVector *Y = psVectorAlloc (match->n, PS_TYPE_F32); + psVector *x = psVectorAlloc (match->n, PS_TYPE_F32); + psVector *y = psVectorAlloc (match->n, PS_TYPE_F32); + psVector *wt = psVectorAlloc (match->n, PS_TYPE_F32); + + // take the matched stars, first fit + for (int i = 0; i < match->n; i++) { + + pair = match->data[i]; + rawStar = raw->data[pair->i1]; + refStar = ref->data[pair->i2]; + + X->data.F32[i] = rawStar->chip.x; + Y->data.F32[i] = rawStar->chip.y; + + x->data.F32[i] = refStar->FP.x; + y->data.F32[i] = refStar->FP.y; + + wt->data.F32[i] = 1.0; + } + + // no masking, no errors + psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, wt, x, y); + psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, wt, x, y); + psFree (x); + psFree (y); + psFree (X); + psFree (Y); + psFree (wt); + psFree (stats); + + return (map); +} + + + +/****************************************************************************** +pmAstromRotateObj(old, center, angle): rotate the focal-plane coordinates +about the center coordinate angle specified in radians + ******************************************************************************/ +psArray *pmAstromRotateObj( + psArray *old, + psPlane center, + double angle) +{ + PS_ASSERT_PTR_NON_NULL(old, NULL); + + double X, Y; + pmAstromObj *newObj; + pmAstromObj *oldObj; + + psArray *new = psArrayAlloc (old->n); + + double cs = cos(angle); + double sn = sin(angle); + double xCenter = center.x; + double yCenter = center.y; + + for (int i = 0; i < old->n; i++) { + + oldObj = (pmAstromObj *)old->data[i]; + newObj = pmAstromObjCopy (oldObj); + + X = oldObj->FP.x - xCenter; + Y = oldObj->FP.y - yCenter; + + newObj->FP.x = X*cs + Y*sn + xCenter; + newObj->FP.y = Y*cs - X*sn + yCenter; + + new->data[i] = newObj; + } + return (new); +} + + + +/****************************************************************************** +pmAstromObjFree(obj) + ******************************************************************************/ +static void pmAstromObjFree(pmAstromObj *obj) +{ + if (obj == NULL) { + return; + } + + return; +} + + +/****************************************************************************** +pmAstromObjAlloc() + ******************************************************************************/ +pmAstromObj *pmAstromObjAlloc(void) +{ + pmAstromObj *obj = psAlloc (sizeof(pmAstromObj)); + psMemSetDeallocator (obj, (psFreeFunc) pmAstromObjFree); + + obj->pix.x = 0; + obj->pix.y = 0; + + obj->FP.x = 0; + obj->FP.y = 0; + + obj->TP.x = 0; + obj->TP.y = 0; + + obj->sky.r = 0; + obj->sky.d = 0; + + obj->Mag = 0; + obj->dMag = 0; + + return (obj); +} + + +/****************************************************************************** +pmAstromObjCopy(old) + ******************************************************************************/ +pmAstromObj *pmAstromObjCopy(pmAstromObj *old) +{ + PS_ASSERT_PTR_NON_NULL(old, NULL); + pmAstromObj *obj = pmAstromObjAlloc(); + + obj[0] = old[0]; + + return(obj); +} + + +/****************************************************************************** + ******************************************************************************/ +static void pmAstromMatchFree (pmAstromMatch *match) +{ + if (match == NULL) + return; + return; +} + + +/****************************************************************************** + ******************************************************************************/ +pmAstromMatch *pmAstromMatchAlloc( + int i1, + int i2) +{ + pmAstromMatch *match = psAlloc (sizeof(pmAstromMatch)); + psMemSetDeallocator(match, (psFreeFunc) pmAstromMatchFree); + + match->i1 = i1; + match->i2 = i2; + + return (match); +} + + +static double maxOffpix; // maximum allowed offset between lists, in raw pixels +static double Scale; // grid pixel scale static +double Offset; // deltas to pixels +/****************************************************************************** +AstromGridBin(*dx, *dy, dX, dY): local function to convert x,y coords to grid +bins it requires the globals defined above. + + ******************************************************************************/ +static bool AstromGridBin( + int *dx, + int *dy, + double dX, + double dY) +{ + if (PM_ASTROMETRYOBJECTS_DEBUG) { + PS_ASSERT_PTR_NON_NULL(dx, false); + PS_ASSERT_PTR_NON_NULL(dy, false); + } + + if (fabs(dX) > maxOffpix) { + return false; + } + + if (fabs(dY) > maxOffpix) { + return false; + } + + *dx = dX / Scale + Offset; + *dy = dY / Scale + Offset; + return true; +} + + +/****************************************************************************** +pmAstromGridAngle(raw, ref, config): match the two lists using the binned +delta-delta max. + ******************************************************************************/ +pmAstromStats pmAstromGridAngle( + psArray *raw, + psArray *ref, + 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); + + bool status; + int nPix; // size of matching grid + int nPixHalf; // half-size of matching grid + double dX, dY; // offset between a possible matched pair + int iX, iY; // corresponding grid bin + + pmAstromObj *ob1, *ob2; // short-cut pointers to the objects + pmAstromStats stats; // output match statistics + + // max allowed offset in either X or Y directions + double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET"); + + // sampling scale of the grid + double gridScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE"); + + // set the static scaling factors + nPixHalf = (int)(gridOffset / gridScale + 0.5); // half-grid + nPix = 2*nPixHalf + 1; // full grid width + + // these are globals used by p_pmAstromGridBin + maxOffpix = gridScale * (nPixHalf + 0.5); // max offset from true center + Offset = maxOffpix / gridScale; + Scale = gridScale; + + // images used as accumulators for the loop below + psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_U32); + psImage *gridDX = psImageAlloc (nPix, nPix, PS_TYPE_F32); + psImage *gridDY = psImageAlloc (nPix, nPix, PS_TYPE_F32); + psImage *gridD2 = psImageAlloc (nPix, nPix, PS_TYPE_F32); + psImageInit (gridNP, 0); + psImageInit (gridDX, 0); + psImageInit (gridDY, 0); + psImageInit (gridD2, 0); + + // short-cut names for grid images + psS32 **NP = gridNP->data.S32; + psF32 **DX = gridDX->data.F32; + psF32 **DY = gridDY->data.F32; + psF32 **D2 = gridD2->data.F32; + + // accumulate grids for focal plane (L,M) matches + for (int i = 0; i < raw->n; i++) { + ob1 = (pmAstromObj *)raw->data[i]; + for (int j = 0; j < ref->n; j++) { + ob2 = (pmAstromObj *)ref->data[j]; + dX = ob1->FP.x - ob2->FP.x; + dY = ob1->FP.y - ob2->FP.y; + + // fprintf (stderr, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP.x, ob2->FP.x, ob1->FP.y, ob2->FP.y); + // find bin coordinates for this delta-delta + if (!AstromGridBin (&iX, &iY, dX, dY)) { + continue; // matched pair is too far offset + } + + // accumulate bin stats + NP[iY][iX] ++; + DX[iY][iX] += dX; + DY[iY][iX] += dY; + D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY); + } + } + + // now assess the grid images + { + double minMetric = 1e10; + double minVar = 1e10; + int minX = -1; + int minY = -1; + double metric, var; + + // find the max pixel + psStats *imStats = psStatsAlloc (PS_STAT_MAX); + imStats = psImageStats (imStats, gridNP, NULL, 0); + + // only check bins with at least 1/2 of max bin + int minNpts = 0.5*imStats->max; + fprintf (stderr, "minNpts: %d, max: %d\n", minNpts, (int)(imStats->max)); + + // find the 'best' bin + for (int j = 0; j < gridNP->numRows; j++) + { + for (int i = 0; i < gridNP->numCols; i++) { + + if (NP[j][i] < minNpts) + continue; + + // this metric emphasize a narrow peak with lots of sources over one with few. + var = fabs((D2[j][i]/NP[j][i]) - PS_SQR(DX[j][i]/NP[j][i]) - PS_SQR(DY[j][i]/NP[j][i])); + metric = var / PS_SQR(PS_SQR(NP[j][i])); + + fprintf (stderr, "try : %f %f (%d pts, %f var, %f met)\n", DX[j][i]/NP[j][i], DY[j][i]/NP[j][i], NP[j][i], var, metric); + + if (metric < minMetric) { + minMetric = metric; + minVar = var; + minX = i; + minY = j; + } + } + } + + // 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]; + + // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better + } + return (stats); +} + + + +/****************************************************************************** +pmAstromGridMatch(*raw, *ref, *config): match two star lists. + ******************************************************************************/ + +pmAstromStats pmAstromGridMatch( + psArray *raw, + psArray *ref, + 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); + + bool status; + double xMin, xMax, yMin, yMax; + pmAstromObj *obj; + psArray *rot; + + pmAstromStats minStat, newStat; + psPlane center; + + // find center of the raw field (focal-plane coords) + xMin = yMin = +1e10; + xMax = yMax = -1e10; + for (int i = 0; i < raw->n; i++) { + obj = (pmAstromObj *)raw->data[i]; + xMin = PS_MIN (obj->FP.x, xMin); + xMax = PS_MAX (obj->FP.x, xMax); + yMin = PS_MIN (obj->FP.y, yMin); + yMax = PS_MAX (obj->FP.y, yMax); + } + center.x = 0.5*(xMin + xMax); + 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; + for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { + rot = pmAstromRotateObj (raw, center, angle); + newStat = pmAstromGridAngle (rot, ref, config); + newStat.angle = angle; + newStat.center = center; + if (newStat.minMetric < minStat.minMetric) { + minStat = newStat; + } + psFree (rot); + } + fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric); + return (minStat); +} + + +/****************************************************************************** +pmAstromGridApply(*map, stat): apply the measured FPA offset and rotation +(stat) to the fpa astrom structures. + ******************************************************************************/ +psPlaneTransform *pmAstromGridApply( + psPlaneTransform *map, + pmAstromStats stat) +{ + PS_ASSERT_PTR_NON_NULL(map, NULL); + PS_ASSERT_POLY_NON_NULL(map->x, NULL); + 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); + + // 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; + + // original rotation matrix + double pc1_1 = map->x->coeff[1][0]; + double pc1_2 = map->x->coeff[0][1]; + double pc2_1 = map->y->coeff[1][0]; + double pc2_2 = map->y->coeff[0][1]; + + // new rotation matrix + map->x->coeff[1][0] = +cs*pc1_1 + sn*pc2_1; + map->x->coeff[0][1] = +cs*pc1_2 + sn*pc2_2; + map->y->coeff[1][0] = -sn*pc1_1 + cs*pc2_1; + map->y->coeff[0][1] = -sn*pc1_2 + cs*pc2_2; + + return (map); +} + +/* Illustration of the grid bins + dX bin + -35:-25 -> 0 bin = dX / Scale + Offset + -25:-15 -> 1 Scale = 10 + -15:-05 -> 2 Offset = 3.5 + -05:+05 -> 3 nPix = 3 (maxOffset = 35 = (nPix + 0.5)*dXix + +05:+15 -> 4 dPix = 10 + +15:+25 -> 5 + +25:+35 -> 6 + + maxOffsetRequest = 30 + nPix = (int) (maxOffset / dPix + 0.5); + maxOffset = (nPix + 0.5)*Scale; +*/ + Index: /trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 5674) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 5674) @@ -0,0 +1,390 @@ +/** @file pmAstrometryObjects.h +* +* @brief This file defines the basic types for matching objects +* based on their astrometry. +* +* @ingroup AstroImage +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-05 20:49:30 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +#if !defined(PM_ASTROMETRYOBJECTS_H) +#define PM_PM_ASTROMETRYOBJECTS_H_H + +# include +# include // for strcasecmp +# include // for unlink +# include +# include + +/* + * + * This structure specifies the coordinate of the detection in each of the + * four necessary coordinate frames: pix defines the position in the psReadout + * frame, FP defines the position in the Focal Plane frame, TP defines the + * position in the Tangent Plane frame, sky defines the position on the Celestial + * Sphere. In addition, a measurement of the brightness is given by the element + * Mag. Such a data structure should be used for both the raw and the reference + * stars. In astrometric processing, the raw detections will be projected using + * the best available information to each of these coordinate frames from the pix + * coordinates, while the reference detections will be projected to the other + * frames from the sky coordinates. + * + * XXX: There are more members here than in the SDRS. + * + */ +typedef struct +{ + psPlane pix; ///< the position in the pmReadout frame + psPlane cell; ///< the position in the pmCell frame + psPlane chip; ///< the position in the pmChip frame + psPlane FP; ///< the position in the pmFPA frame + psPlane TP; ///< the position in the tangent plane + psSphere sky; ///< the position on the Celestial Sphere. + double Mag; ///< A measurement of the brightness. + double dMag; ///< What is this? +} +pmAstromObj; + + + +/* + * + * The pmAstromMatch structure defines the cross-correlation between two + * arrays. A single such data item specifies that item number pmAstromMatch.idx1 + * in the first list corresponds to pmAstromMatch.idx2 in the second list. + * + */ +typedef struct +{ + int i1; ///< What is this? + int i2; ///< What is this? +} +pmAstromMatch; + + +/* + * + * XXX: Not in SDRS. + * + */ +typedef struct +{ + psPlane center; ///< + psPlane offset; ///< + double angle; ///< + double minMetric; ///< + double minVar; ///< + int nMatch; ///< +} +pmAstromStats; + + + +/* + * + * If the two sets of coordinates are expected to agree very well (ie, the + * current best-guess astrometric solution is quite close to the radius. The + * following function accepts two sets of pmAstromObj sources and determines the + * matched objects between the two lists. The input sources must have been + * projected to the Focal Plane coordinates (pmAstromObj.FP), and the supplied + * options entry must contain the desired match radius (keyword: + * ASTROM.MATCH.RADIUS). The output consists an array of pmAstromMatch values, + * defined below. + * + */ +psArray *pmAstromRadiusMatch( + psArray *st1, + psArray *st2, + psMetadata *config +); + + + +/* + * + * This function accepts an array of pmAstromObj objects and rotates them by + * the given angle about the given center coordinate pCenter,qCenter in the Focal + * Plane Array coordinates. + * + * XXX: This differs from the SDRS + * + */ +/* SDRS +psArray *pmAstromRotateObj( + psArray *old, + double angle, + double pCenter, + double qCenter +); +*/ +psArray *pmAstromRotateObj( + psArray *old, + psPlane center, + double angle +); + + +/* + * + * If the two sets of coordinates are not known to agree well, but the + * relative scale and approximate relative rotation is known, then a much faster + * match can be found using pair-pair displacements. In such a case, the two + * lists can be considered as having the same coordinate system, with an unknown + * relative displacement. In this algorithm, all possible pair-wise differences + * between the source positions in the two lists are constructed and accumulated + * in a grid of possible offset values. The resulting grid is searched for a + * cluster representing the offset between the two input lists. This algorithm + * can only tolerate a small error in the relative scale or the relative rotation + * of the two coordinate lists. However, this process is naturally O(N2), and is + * thus advantageous over triangle matching in some circumstances. This process + * can be extended to allow a larger uncertainty in the relative rotation by + * allowing the procedure to scan over a range of rotations. We define the + * following function to apply this matching algorithm: + * + * XXX: In the SDRS, this function is a pointer. + * + */ +pmAstromStats pmAstromGridMatch( + psArray *st1, + psArray *st2, + psMetadata *config +); + + +/* + * + * The result of a pmAstromGridMatch may be used to modify the astrometry + * transformation information for a pmFPA image hierarchy structure. The result + * of pmAstromGridMatch defines the adjustments which should be made to the + * reference coordinate of the projection (pmFPA.projection.R,D) and the + * effective rotation of the Focal Plane. The rotation implies modification of + * the linear terms of the pmFPA.toTangentPlane transformation. These two + * adjustments are made using the function: + * + * XXX: This function name is different in the SDRS. + * + */ +psPlaneTransform *pmAstromGridApply( + psPlaneTransform *map, + pmAstromStats stat +); + + +/* + * + * This function is identical to pmAstromGridMatch, but is valid for only a + * single relative rotation. The input config information need not contain any of + * the GRID.*.ANGLE entries (they will be ignored). + * + * XXX: This function name is different in the SDRS. + * + */ +/* in pmAstromGrid.c */ +pmAstromStats pmAstromGridAngle( + psArray *st1, + psArray *st2, + psMetadata *config); + + + +/* + * + * This function accepts the raw and reference source lists and the list of + * matched entries. It uses the matched list to determine a polynomial + * transformation between the two coordinate systems. The fitting uses clipping + * to exclude outliers, likely representing poor matches. The config element must + * contain the information ASTROM.NSIGMA (specifying the number of sigma used in + * the clipping) and ASTROM.NCLIP (specifying the number of clipping iterations + * must be performed). The config element must also specify the order of the + * polynomial fit (keyword: ASTROM.ORDER). The result of this fit is a set of + * modifications of the components of the pmFPA.toTangentPlane transformation, + * and the modifications of the reference coordinate of the projection + * (pmFPA.projection.R,D) and the projection scale (pmFPA.projection.Xs,Ys). The + * modifications to pmFPA.toTangentPlane incorporate the rotation component of + * the linear terms and the higher-order terms of the polynomial fits. + * + * XXX: No prototype code. + * + */ +bool pmAstromFitFPA( + pmFPA *fpa, + psArray *st1, + psArray *st2, + psArray *match, + psMetadata *config +); + + + +/* + * + * This function accepts the raw and reference source lists for a single chip + * and the list of matched entries. It uses the matched list to determine a + * polynomial transformation between the two coordinate systems. The fitting + * uses clipping to exclude outliers, likely representing poor matches. The + * config element must contain the information ASTROM.NSIGMA + *(specifying the number of sigma used in the clipping) and ASTROM.NCLIP + *(specifying the number of clipping iterations must be performed). The config + *element must also specify the order of the polynomial fit (keyword: + *ASTROM.ORDER). The result of this fit is a set of modifications of the + *components of the pmChip.toFPA transformation. + * + * XXX: No prototype code. + * + */ +bool pmAstromFitChip( + pmFPA *fpa, + psArray *st1, + psArray *st2, + psArray *match, + psMetadata *config +); + + +/* + * + * The following function determines the position residual, in the tangent + * plane, as a function of position in the focal plane, for a collection of raw + * measurements and matched reference stars. The configuration data must include + * the bin size over which the gradient is measured (keyword: ASTROM.GRAD.BOX). + * The function returns an array of pmAstromGradient structures, defined below. + * + * XXX: No prototype code. + * + */ +psArray pmAstromMeasureGradients( + psArray *starlist1, + psArray *starlist2, + psArray *match, + psMetadata *config +); + + + +/* + * + * The following data structure carries the information about the residual + * gradient of source positions in the tangent plane (pmAstromObj.TP) as a + * function of position in the focal plane (pmAstromObj.FP). + * + */ +typedef struct +{ + psPlane FP; + psPlane dTPdL; + psPlane dTPdM; +} +pmAstromGradient; + + +/* + * + * The gradient set measured above can be fitted with a pair of 2D + * polynomials. The resulting fits can then be related back to the implied + * polynomials which represent the distortion. The following function performs + * the fit and applies the result to the distortion transformation of the + * supplied pmFPA structure. The configuration variable supplies the polynomial + * order (keyword: ASTROM.DISTORT.ORDER). + * + * XXX: No prototype code. + * + */ +psArray pmAstromFitDistortion( + pmFPA *fpa, + psArray *gradients, + psMetadata *config); + + + + + +/******************************************************************************* + The following functions and structs were in the prototype code, but not the + SDRS. + ******************************************************************************/ +/* + * + * + * + * + */ + + +/* + * + * Allocates a pmAstromObj struct. + * + */ +pmAstromObj *pmAstromObjAlloc (void); + + + +/* + * + * Copies a pmAstromObj struct. + * + */ +pmAstromObj *pmAstromObjCopy( + pmAstromObj *old +); + + + +/* + * + * + * + */ +pmAstromMatch *pmAstromMatchAlloc( + int i1, + int i2 +); + + + + +/* + * + * + * + */ +psPlaneTransform *pmAstromMatchFit( + psPlaneTransform *map, + psArray *st1, + psArray *st2, + psArray *match, + psMetadata *config +); + + + +/* + * + * + * + */ +int pmAstromObjSortByFPX( + const void **a, + const void **b +); + + + +/* + * + * + * + */ +int pmAstromObjSortByMag( + const void **a, + const void **b +); + +#endif +