Index: trunk/psastro/src/pmAstromGrid.c
===================================================================
--- trunk/psastro/src/pmAstromGrid.c	(revision 5510)
+++ trunk/psastro/src/pmAstromGrid.c	(revision 5560)
@@ -1,5 +1,5 @@
 # include "pmAstrom.h"
 
-static double maxOffset;  // maximum allowed offset between lists, in raw pixels
+static double maxOffpix;  // maximum allowed offset between lists, in raw pixels
 static double Scale;      // grid pixel scale
 static double Offset;     // deltas to pixels
@@ -7,113 +7,75 @@
 // local function to convert x,y coords to grid bins
 // it requires the globals defined above
-static bool AstromGridBin (int *pOut, int *qOut, double dP, int dQ) {
-
-    if (fabs(dP) > maxOffset) return false;
-    if (fabs(dQ) > maxOffset) return false;
-
-    *pOut = dP / Scale + Offset;
-    *qOut = dQ / Scale + Offset;
+static bool AstromGridBin (int *dx, int *dy, double dX, double dY) {
+
+    if (fabs(dX) > maxOffpix) return false;
+    if (fabs(dY) > maxOffpix) return false;
+
+    *dx = dX / Scale + Offset;
+    *dy = dY / Scale + Offset;
     return true;
 }
 
-// match two star lists
-pmAstromGridMatchStats pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config) {
-
-    double xMin, xMax, yMin, yMax;
-    pmAstromObj *ob1, *ob2;
-
-    pmAstromGridMatchStat minStat, newStat;
-    psPlane center;
-
-    // find center of the st2 field (focal-plane coords)
-    xMin = yMin = +1e10; 
-    xMax = yMax = -1e10;
-    for (int i = 0; i < st2->n; i++) {
-	ob2 = (pmAstromObj *)st2->data[i];
-	xMin = PS_MIN (ob2->FP.x, xMin);
-	xMax = PS_MAX (ob2->FP.x, xMax);
-	yMin = PS_MIN (ob2->FP.y, yMin);
-	yMax = PS_MAX (ob2->FP.y, yMax);
-    }
-    center.x = 0.5*(xMin + xMax);
-    center.y = 0.5*(yMin + yMax);
-
-    double minAngle = psMetadataLookupF32 (&status, config, "GRID.MIN.ANGLE");
-    double maxAngle = psMetadataLookupF32 (&status, config, "GRID.MAX.ANGLE");
-    double delAngle = psMetadataLookupF32 (&status, config, "GRID.DEL.ANGLE");
-
-    minStat.minMetric = 1e10;
-    for (angle = minAngle; angle <= maxAngle; angle += delAngle) {
-	st2r = pmAstromRotateObj (st2, center, angle);
-	newStat = pmAstromGridMatchAngle (st1, st2r, config);
-	newStat.angle  = angle;
-	newStat.center = center;
-	if (newStat.minMetric < minStat.minMetric) {
-	    minStat = newStat;
-	}
-	psFree (st2r);
-    }
-    return (minStat);
-}
-
 // match the two lists using the binned delta-delta max
-pmAstromGridMatchStats pmAstromGridMatchAngle (psArray *st1, psArray *st2, psMetadata *config) {
-
+pmAstromStats pmAstromGridAngle (psArray *raw, psArray *ref, psMetadata *config) {
+
+    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
-    pmAstromGridMatchStats matchStats; // output match statistics
+    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, "GRID.OFFSET"); 
+    double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET"); 
 
     // sampling scale of the grid
-    double gridScale  = psMetadataLookupF32 (&status, config, "GRID.SCALE");
+    double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
 
     // set the static scaling factors
-    nPix = (int)(gridOffset / gridScale + 0.5);  // half-grid
-    nPix = 2*nPix + 1;
+    nPixHalf = (int)(gridOffset / gridScale + 0.5);  // half-grid
+    nPix = 2*nPixHalf + 1;                           // full grid width
 
     // these are globals used by p_pmAstromGridBin
-    maxOffpix = gridScale * (nPix + 0.5);
+    maxOffpix = gridScale * (nPixHalf + 0.5);            // max offset from true center
     Offset    = maxOffpix / gridScale; 
     Scale     = gridScale;
 
-    // XXX EAM : can we assume the allocated image is init-ed?
-    psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_S32);
+    // 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);
-    psImageInit (gridDX);
-    psImageInit (gridDY);
-    psImageInit (gridD2);
+    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 = gridDP->data.F32;
-    psF32 *DY = gridDQ->data.F32;
-    psF32 *D2 = gridD2->data.F32;
+    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 < st1->n; i++) {
-	ob1 = (pmAstromObj *)st1->data[i];
-	for (int j = 0; j < st2->n; j++) {
-	    ob2 = (pmAstromObj *)st2->data[i];
+    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 (!p_pmAstromGridBin (&iX, &iY, dX, dY)) {
+	    if (!AstromGridBin (&iX, &iY, dX, dY)) {
 		continue; // matched pair is too far offset
-
 	    }
 
 	    // accumulate bin stats
 	    NP[iY][iX] ++;
-	    DP[iY][iX] += dX;
-	    DQ[iY][iX] += dY;
+	    DX[iY][iX] += dX;
+	    DY[iY][iX] += dY;
 	    D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY);
 	}
@@ -124,18 +86,19 @@
 	double minMetric = 1e10;
 	double minVar = 1e10;
-	double minX = -1;
-	double minY = -1;
+	int minX = -1;
+	int minY = -1;
 	double metric, var;
 
 	// find the max pixel
-	psStats *stats = psStatsAlloc (PS_STAT_MAX);
-	stats = psImageStats (stats, gridNP);
+	psStats *imStats = psStatsAlloc (PS_STAT_MAX);
+	imStats = psImageStats (imStats, gridNP, NULL, 0);
 
 	// only check bins with at least 1/2 of max bin
-	int minNpt = 0.5*stats->max;
+	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->nY; j++) {
-	    for (int i = 0; i < gridNP->nX; i++) {
+	for (int j = 0; j < gridNP->numRows; j++) {
+	    for (int i = 0; i < gridNP->numCols; i++) {
 
 		if (NP[j][i] < minNpts) continue;
@@ -145,4 +108,6 @@
 		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;
@@ -155,20 +120,84 @@
 
 	// convert the bin to delta-delta
-	matchStats.offset.x  = DP[minY][minX] / NP[minY][minX];
-	matchStats.offset.y  = DQ[minY][minX] / NP[minY][minX];
-	matchStats.minMetric = minMetric;
-	matchStats.minVar    = minVar;
-	matchStats.nMatch    = NP[minY][minX];
-    }
-    return (matchStats);
+	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);
+}
+
+// match two star lists
+pmAstromStats pmAstromGridMatch (psArray *raw, psArray *ref, psMetadata *config) {
+
+    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);
 }
 
 // apply the measured FPA offset and rotation (stat) to the fpa astrom structures
-psFPA *pmAstromGridApply (psPlaneTransform *map, pmAstromGridMatchStat stat) {
-
-    // stat.angle, stat.center, stat.offse
-    // I think i need to know the center reference....
-
-    return (fpa);
+psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat) {
+
+    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);
 }