Index: /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractAnalysis.c
===================================================================
--- /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractAnalysis.c	(revision 23406)
+++ /branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractAnalysis.c	(revision 23407)
@@ -13,9 +13,9 @@
 # include "psastroInternal.h"
 
-# define ESCAPE { \
-  psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate"); \
-  psFree (view); \
-  return false; \
-}
+# define ESCAPE {							\
+	psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate"); \
+	psFree (view);							\
+	return false;							\
+    }
 
 # define TEST_OUTPUT 1
@@ -106,11 +106,11 @@
     char *filename = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
     if (!filename) {
-      psLogMsg ("psastro", PS_LOG_INFO, "output filename not supplieda");
-      return false;
+	psLogMsg ("psastro", PS_LOG_INFO, "output filename not supplieda");
+	return false;
     }
     psFits *outStars = psFitsOpen (filename, "w");
     if (!outStars) {
-      psError(PS_ERR_IO, false, "error opening file %s\n", filename);
-      return false;
+	psError(PS_ERR_IO, false, "error opening file %s\n", filename);
+	return false;
     }
 
@@ -131,13 +131,13 @@
         FILE *f = NULL;
         if (TEST_OUTPUT) {
-          char *filename = NULL;
-          char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
-          psStringAppend (&filename, "bright.%s.dat", chipname);
-          f = fopen (filename, "w");
-          if (!f) {
-            psWarning ("cannot create test output file %s\n", filename);
-            continue;
-          }
-          psFree (filename);
+	    char *filename = NULL;
+	    char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
+	    psStringAppend (&filename, "bright.%s.dat", chipname);
+	    f = fopen (filename, "w");
+	    if (!f) {
+		psWarning ("cannot create test output file %s\n", filename);
+		continue;
+	    }
+	    psFree (filename);
         }
 
@@ -160,5 +160,5 @@
 
                     if (TEST_OUTPUT) {
-                      fprintf (f, "stars %f %f  %f %f   %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->Mag);
+			fprintf (f, "stars %f %f  %f %f   %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->Mag);
                     }
 
@@ -190,22 +190,40 @@
 psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY) {
 
-  // skip if no pixels are on the image
-  // skip if center is not on the image
-  // if bounds fall off image, paste into a full-size image.
-
-  if (x < 0) return NULL;
-  if (y < 0) return NULL;
-  if (x >= input->numCols) return NULL;
-  if (y >= input->numRows) return NULL;
-
-  psRegion fullRegion = psRegionSet (x - dX, x + dX, y - dY, y + dY);
-  psRegion realRegion = psRegionForImage (input, fullRegion);
-  psImage *subset = psImageSubset (input, realRegion);
-
-  // if ((subset->numCols < dX) || (subset->numRows < dY)) {
-  //   psImage *fullsize = psImageAlloc (dX, dY, PS_TYPE_F32);
-  //   // fill in the fullsize image with the values from the subset
-  //   // I must already have done this elsewhere...
-  // }
-  return subset;
+    // skip if no pixels are on the image
+    // skip if center is not on the image
+    // if bounds fall off image, paste into a full-size image.
+
+    if (x < 0) return NULL;
+    if (y < 0) return NULL;
+    if (x >= input->numCols) return NULL;
+    if (y >= input->numRows) return NULL;
+
+    int xo = x; // index of center pixel (eg, 6.00 to 6.99 -> 6)
+    int yo = y;
+
+    // build an output image of size 2dX + 1, 2dY + 1
+    // psRegion fullRegion = psRegionSet (xo - dX, x + dX + 1, y - dY, y + dY + 1);
+    // psRegion realRegion = psRegionForImage (input, fullRegion);
+    // psImage *subset = psImageSubset (input, realRegion);
+
+    psImage *subset = psImageAlloc (2*dX+1, 2*dY+1, PS_TYPE_F32);
+    psImageInit (subset, 0.0);
+
+    // fill in the subset image with the values from the subset
+    // I must already have done this elsewhere...
+
+    for (int iy = yo - dY; iy < yo + dY + 1; iy++) {
+	for (int ix = xo - dX; ix < xo + dX + 1; ix++) {
+
+	    if (ix < input->col0) continue;
+	    if (iy < input->row0) continue;
+	    if (ix >= input->numCols) continue;
+	    if (iy >= input->numRows) continue;
+	
+	    int jx = ix + dX - xo; // runs from 0 to 2dX
+	    int jy = iy + dY - yo;
+	    subset->data.F32[jy][jx] = input->data.F32[iy][ix];
+	}
+    }
+    return subset;
 }