Index: trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 7004) +++ trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 7005) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-04-17 18:01:04 $ +* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-04-30 22:03:58 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,5 @@ #include #include "pslib.h" +#include "psVectorSmooth.h" #include "pmFPA.h" #include "pmAstrometryObjects.h" @@ -46,5 +47,5 @@ pmAstromObj *B = *(pmAstromObj **)b; - psF32 diff = A->FP.x - B->FP.x; + psF32 diff = A->FP->x - B->FP->x; if (diff > FLT_EPSILON) { return (+1); @@ -121,5 +122,5 @@ while ((i < st1->n) && (j < st2->n)) { - dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x; + dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x; if (dX < -RADIUS) { i++; @@ -134,6 +135,6 @@ 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; + 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; @@ -146,4 +147,5 @@ pmAstromMatch *match = pmAstromMatchAlloc (i, j); psArrayAdd (matches, 100, match); + psFree (match); j++; @@ -152,4 +154,5 @@ i++; } + psLogMsg (__func__, 3, "radius match: %d pairs\n", matches->n); return (matches); } @@ -213,16 +216,24 @@ 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; + 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); + // constant errors + psVector *xMask = psVectorAlloc (match->n, PS_TYPE_U8); + psVector *yMask = psVectorAlloc (match->n, PS_TYPE_U8); + xMask->n = yMask->n = match->n; + psVectorInit (xMask, 0); + psVectorInit (yMask, 0); + + // fit chip-to-FPA transformation + psVectorClipFitPolynomial2D (map->x, stats, xMask, 0, x, wt, X, Y); + psVectorClipFitPolynomial2D (map->y, stats, yMask, 0, y, wt, X, Y); + psFree (x); psFree (y); @@ -231,8 +242,9 @@ psFree (wt); psFree (stats); + psFree (xMask); + psFree (yMask); return (map); } - @@ -264,9 +276,9 @@ 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; + 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; @@ -286,4 +298,11 @@ } + psFree(obj->pix); + psFree(obj->cell); + psFree(obj->chip); + psFree(obj->FP); + psFree(obj->TP); + psFree(obj->sky); + return; } @@ -298,17 +317,11 @@ 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->pix = psPlaneAlloc(); + obj->cell = psPlaneAlloc(); + obj->chip = psPlaneAlloc(); + obj->FP = psPlaneAlloc(); + obj->TP = psPlaneAlloc(); + obj->sky = psSphereAlloc(); + obj->Mag = 0; obj->dMag = 0; @@ -325,5 +338,10 @@ pmAstromObj *obj = pmAstromObjAlloc(); - obj[0] = old[0]; + *obj->pix = *old->pix; + *obj->cell = *old->cell; + *obj->chip = *old->chip; + *obj->FP = *old->FP; + *obj->TP = *old->TP; + *obj->sky = *old->sky; return(obj); @@ -450,8 +468,8 @@ 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); + 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)) { @@ -491,9 +509,9 @@ continue; - // this metric emphasize a narrow peak with lots of sources over one with few. + // this metric emphasizes 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); + metric = var / PS_SQR(NP[j][i]) / 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) { @@ -513,6 +531,11 @@ stats.nMatch = NP[minY][minX]; + psFree (imStats); // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better } + psFree (gridNP); + psFree (gridDX); + psFree (gridDY); + psFree (gridD2); return (stats); } @@ -549,8 +572,8 @@ 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); + 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); @@ -565,15 +588,106 @@ rot = pmAstromRotateObj (raw, center, angle); newStat = pmAstromGridAngle (rot, ref, config); + fprintf (stderr, "try %f : %f %f (%d pts, %f var, %f met (%f log))\n", 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; + fprintf (stderr, "use %f : %f %f (%d pts, %f var, %f met (%f log))\n", angle, minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, minStat.minMetric, log10(minStat.minMetric)); } 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); + fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, log10(minStat.minMetric)); return (minStat); } +/****************************************************************************** +pmAstromGridTweak(*raw, *ref, *recipe, stats): improve match for two star lists. + ******************************************************************************/ +pmAstromStats pmAstromGridTweak( + psArray *raw, + psArray *ref, + psMetadata *recipe, + pmAstromStats stats) +{ + bool status; + pmAstromObj *ob1, *ob2; // short-cut pointers to the objects + double dX, dY; // offset between a possible matched pair + psArray *rot; + int nBin, xBin, yBin; + + rot = pmAstromRotateObj (raw, stats.center, stats.angle); + + // sampling scale of the grid + double tweakScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.SCALE"); + double tweakRange = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.RANGE"); + double tweakSmooth = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.SMOOTH"); + double tweakNsigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.NSIGMA"); + + nBin = 2*tweakRange / tweakScale; + psVector *xHist = psVectorAlloc (nBin, PS_TYPE_F32); + psVector *yHist = psVectorAlloc (nBin, PS_TYPE_F32); + xHist->n = yHist->n = nBin; + psVectorInit (xHist, 0); + psVectorInit (yHist, 0); + + // accumulate grids for focal plane (L,M) matches + for (int i = 0; i < rot->n; i++) { + ob1 = (pmAstromObj *)rot->data[i]; + 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; + + xBin = (dX + tweakRange) / tweakScale; + yBin = (dY + tweakRange) / tweakScale; + + if (xBin < 0) + continue; + if (yBin < 0) + continue; + if (xBin >= nBin) + continue; + if (yBin >= nBin) + continue; + + xHist->data.F32[xBin] += 1.0; + yHist->data.F32[yBin] += 1.0; + } + } + + // smooth histgram vector with gaussian of 1sigma = radius + psVectorSmooth (xHist, tweakSmooth, tweakNsigma); + psVectorSmooth (yHist, tweakSmooth, tweakNsigma); + + // select peak in x and in y + xBin = yBin = 0; + double xMax = 0; + double yMax = 0; + for (int i = 0; i < nBin; i++) { + if (xHist->data.F32[i] > xMax) { + xBin = i; + xMax = xHist->data.F32[i]; + } + if (yHist->data.F32[i] > yMax) { + yBin = i; + yMax = yHist->data.F32[i]; + } + } + double xPeak = xBin*tweakScale - tweakRange; + double yPeak = yBin*tweakScale - tweakRange; + psLogMsg (__func__, 3, "tweak peak by %f,%f\n", xPeak, yPeak); + + // adjust offset by peak center + pmAstromStats tweak = stats; + tweak.offset.x += xPeak; + tweak.offset.y += yPeak; + + psFree (rot); + psFree (xHist); + psFree (yHist); + + return tweak; +} /****************************************************************************** Index: trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 7004) +++ trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 7005) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-04-17 18:10:08 $ +* @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-04-30 22:03:58 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -41,16 +41,14 @@ 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? + 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; ///< object magnitude XXX what filter? + double dMag; ///< error on object magnitude } pmAstromObj; - - /* @@ -158,4 +156,13 @@ +/****************************************************************************** +pmAstromGridTweak(*raw, *ref, *recipe, stats): improve match for two star lists. + ******************************************************************************/ +pmAstromStats pmAstromGridTweak( + psArray *raw, + psArray *ref, + psMetadata *recipe, + pmAstromStats stats); + /* *