Index: trunk/psModules/src/astrom/pmAstrometryDistortion.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryDistortion.c (revision 10827) +++ trunk/psModules/src/astrom/pmAstrometryDistortion.c (revision 10829) @@ -7,6 +7,6 @@ * @author EAM, IfA * -* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-22 21:23:06 $ +* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-25 01:49:49 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -50,9 +50,9 @@ } -psArray *pmAstromMeasureGradients(psArray *grads, psArray *rawstars, psArray *refstars, psArray *matches, psMetadata *config) +psArray *pmAstromMeasureGradients(psArray *gradients, psArray *rawstars, psArray *refstars, psArray *matches, psMetadata *config) { - if (grads == NULL) { - grads = psArrayAllocEmpty (100); + if (gradients == NULL) { + gradients = psArrayAllocEmpty (100); } @@ -141,11 +141,11 @@ grad->dTPdM.y = local->coeff[0][1]; - psArrayAdd (grads, 100, grad); - } - } - return grads; + psArrayAdd (gradients, 100, grad); + } + } + return gradients; } -bool pmAstromFitDistortion(pmFPA *fpa, psArray *grads, psMetadata *config) +bool pmAstromFitDistortion(pmFPA *fpa, psArray *gradients, psMetadata *config) { @@ -154,14 +154,14 @@ // assign the gradient elements to psVectors for fitting - psVector *dPdL = psVectorAlloc (grads->n, PS_TYPE_F32); - psVector *dQdL = psVectorAlloc (grads->n, PS_TYPE_F32); - psVector *dPdM = psVectorAlloc (grads->n, PS_TYPE_F32); - psVector *dQdM = psVectorAlloc (grads->n, PS_TYPE_F32); - psVector *L = psVectorAlloc (grads->n, PS_TYPE_F32); - psVector *M = psVectorAlloc (grads->n, PS_TYPE_F32); - - for (int i = 0; i < grads->n; i++) { - - pmAstromGradient *grad = grads->data[i]; + psVector *dPdL = psVectorAlloc (gradients->n, PS_TYPE_F32); + psVector *dQdL = psVectorAlloc (gradients->n, PS_TYPE_F32); + psVector *dPdM = psVectorAlloc (gradients->n, PS_TYPE_F32); + psVector *dQdM = psVectorAlloc (gradients->n, PS_TYPE_F32); + psVector *L = psVectorAlloc (gradients->n, PS_TYPE_F32); + psVector *M = psVectorAlloc (gradients->n, PS_TYPE_F32); + + for (int i = 0; i < gradients->n; i++) { + + pmAstromGradient *grad = gradients->data[i]; dPdL->data.F32[i] = grad->dTPdL.x; @@ -176,5 +176,5 @@ psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - psVector *mask = psVectorAlloc (grads->n, PS_TYPE_MASK); + psVector *mask = psVectorAlloc (gradients->n, PS_TYPE_MASK); psVectorInit (mask, 0); @@ -185,9 +185,13 @@ // set masks based on fpa->toTP - for (int i = 0; i < fpa->toTPA->x->nX; i++) { - for (int j = 0; j < fpa->toTPA->x->nY; j++) { + for (int i = 0; i <= fpa->toTPA->x->nX; i++) { + for (int j = 0; j <= fpa->toTPA->x->nY; j++) { if (fpa->toTPA->x->mask[i][j]) { - localX->mask[i-1][j] = 1; - localY->mask[i][j-1] = 1; + if ((i > 0) && (i < fpa->toTPA->x->nX)) { + localX->mask[i-1][j] = 1; + } + if ((j > 0) && (j < fpa->toTPA->y->nY)) { + localY->mask[i][j-1] = 1; + } } } @@ -247,5 +251,5 @@ psVector *dQdMf = psPolynomial2DEvalVector (localY, L, M); - FILE *f = fopen ("grads.dat", "w"); + FILE *f = fopen ("gradients.dat", "w"); for (int i = 0; i < dPdLf->n; i++) { fprintf (f, "%f %f | %f %f %f %f | %f %f %f %f\n", @@ -290,4 +294,7 @@ psFree (mask); + // XXX need to reset the fromFPA terms here + + return true; } Index: trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 10827) +++ 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); } Index: trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 10827) +++ trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 10829) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-17 09:45:20 $ +* @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-25 01:49:49 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -79,14 +79,23 @@ 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. +typedef struct +{ + psStats *xStats; + psStats *yStats; + int nMatch; ///< + double nSigma; ///< +} +pmAstromFitResults; + +/* + * + * If the two sets of coordinates are expected to agree very well (ie, the current best-guess + * astrometric solution is quite close to reality), perform a match based on a simple radius + * test. The following functions accept two sets of pmAstromObj sources and determines the + * matched objects between the two lists using coordinates of the desired depth (depending on + * the function). The input and reference sources must have been projected to the desired depth + * (eg, for Focal Plane coordinates, to pmAstromObj.FP). The specified radius must be in the + * units for the matching depth (chip: pixels, focal plane: microns, tangent plane: + * degrees. The output consists an array of pmAstromMatch values, defined above. * */ @@ -94,15 +103,20 @@ const psArray *st1, const psArray *st2, - const psMetadata *config + double RADIUS +); +psArray *pmAstromRadiusMatchFP( + const psArray *st1, + const psArray *st2, + double RADIUS ); psArray *pmAstromRadiusMatchTP( const psArray *st1, const psArray *st2, - const psMetadata *config + double RADIUS ); psArray *pmAstromRadiusMatchChip( const psArray *st1, const psArray *st2, - const psMetadata *config + double RADIUS ); @@ -270,4 +284,5 @@ */ +pmAstromFitResults *pmAstromFitResultsAlloc(void); /* @@ -312,11 +327,10 @@ * */ -bool pmAstromMatchFit( +pmAstromFitResults *pmAstromMatchFit( psPlaneTransform *map, - psArray *st1, - psArray *st2, + psArray *raw, + psArray *ref, psArray *match, - psMetadata *config, - psMetadata *updates + psStats *stats ); Index: trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10827) +++ trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10829) @@ -7,6 +7,6 @@ * @author EAM, IfA * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-22 21:23:06 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-25 01:49:49 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -32,6 +32,6 @@ // interpret header WCS (only handles traditional WCS for the moment) -// plateScale is nominal physical scale on tangent plane (radians / TPA physical units) -bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale) +// pixelScale is microns per pixel +bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double pixelScale) { pmAstromWCS *wcs = pmAstromWCSfromHeader (header); @@ -40,5 +40,5 @@ } - bool status = pmAstromWCStoFPA (fpa, chip, wcs, plateScale); + bool status = pmAstromWCStoFPA (fpa, chip, wcs, pixelScale); psFree (wcs); @@ -67,5 +67,4 @@ // interpret chip header WCS as bilevel chip components -// plateScale is nominal physical scale on tangent plane (radians / TPA physical units) bool pmAstromReadBilevelChip (pmChip *chip, const psMetadata *header) { @@ -82,5 +81,4 @@ // convert toFPA / toSky components to traditional WCS -// plateScale is nominal physical scale on tangent plane (microns / arcsecond) bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header) { @@ -359,6 +357,6 @@ // interpret header WCS (only handles traditional WCS for the moment) -// plateScale is nominal physical scale on tangent plane (radians / TPA physical units) -bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale) +// pixelScale is the pixel size in microns per pixel +bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double pixelScale) { psPlaneTransform *toFPA; @@ -367,37 +365,37 @@ * wcs->trans, which will convert X,Y in pixels to L,M in degrees. we also have the following * elements defined: - * type (CTYPE) + * type (projection type) * crval1,2 (in RA,DEC degrees) * crpix1,2 * cdelt1,2 (in degrees / pixel) - * plateScale (radians / physical TPA units) + * pixelScale (microns / pixel) * * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways: - * 1) the output is in pixels (not degrees): divide by cdelt1,2 raised to an appropriate power + * 1) the output is in microns (not degrees): divide by cdelt1,2 * 2) X,Y are applied directly, without an applied Xo,Yo offset * 3) there is an allowed Lo,Mo term ([0][0] coefficients) */ - // convert wcs->trans to a matrix which yields L,M in pixels - double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]); - double cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]); - for (int i = 0; i <= wcs->trans->x->nX; i++) { - for (int j = 0; j <= wcs->trans->x->nX; j++) { - wcs->trans->x->coeff[i][j] /= cdelt1; - wcs->trans->y->coeff[i][j] /= cdelt2; - } - } - // create transformation with 0,0 reference pixel toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2); - // scale from FPA to TPA (microns / pixel) - double pdelt1 = cdelt1*PM_RAD_DEG / plateScale; - double pdelt2 = cdelt2*PM_RAD_DEG / plateScale; + // modify scale of toFPA to yield L,M in microns + double cdelt1 = hypot (toFPA->x->coeff[1][0], toFPA->x->coeff[0][1]); + double cdelt2 = hypot (toFPA->y->coeff[1][0], toFPA->y->coeff[0][1]); + for (int i = 0; i <= toFPA->x->nX; i++) { + for (int j = 0; j <= toFPA->x->nX; j++) { + toFPA->x->coeff[i][j] *= pixelScale/cdelt1; + toFPA->y->coeff[i][j] *= pixelScale/cdelt2; + } + } + + // scale from FPA to TPA (degrees / micron) + double pdelt1 = cdelt1 / pixelScale; + double pdelt2 = cdelt2 / pixelScale; float rX = 1.0; float rY = 1.0; - // projection from TPA to SKY - psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, plateScale, plateScale, wcs->toSky->type); + // projection from TPA ("linear" degrees) to SKY (radians) + psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, PM_RAD_DEG, PM_RAD_DEG, wcs->toSky->type); if (fpa->toSky == NULL) { @@ -558,6 +556,6 @@ // cdelt1,2 convert from pixels->degrees - double cdelt1 = fpa->toTPA->x->coeff[1][0]*fpa->toSky->Xs*PM_DEG_RAD; - double cdelt2 = fpa->toTPA->y->coeff[0][1]*fpa->toSky->Ys*PM_DEG_RAD; + double cdelt1 = fpa->toTPA->x->coeff[1][0]; + double cdelt2 = fpa->toTPA->y->coeff[0][1]; wcs->cdelt1 = cdelt1; wcs->cdelt2 = cdelt2; Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 10827) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 10829) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-03 18:48:10 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-25 01:49:50 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -256,4 +256,11 @@ if (phu) { outhead = psMetadataCopy (NULL, phu->header); + + // if we have mosaic-level astrometry information, add it here: + updates = psMetadataLookupPtr (&status, file->fpa->analysis, "PSASTRO.HEADER"); + if (updates) { + psMetadataCopy (outhead, updates); + } + psMetadataAddBool (outhead, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true); // XXX what is the EXTNAME??