Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c	(revision 14773)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c	(revision 14774)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-03 20:29:07 $
+ *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:19:35 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -20,6 +20,10 @@
 #include "psError.h"
 #include "psAbort.h"
+
+#include "psFits.h"
 #include "psAssert.h"
 #include "psRegion.h"
+#include "psFitsImage.h"
+
 #include "psVector.h"
 #include "psImage.h"
@@ -28,4 +32,5 @@
 #include "psImageMap.h"
 #include "psImagePixelInterpolate.h"
+#include "psImageUnbin.h"
 
 static void psImageMapFree(psImageMap *map) {
@@ -34,5 +39,5 @@
 
     psFree (map->map);
-    psFree (map->field);
+    // psFree (map->field); XXX ??? this should be freed here, but that causes an error...
     psFree (map->stats);
     psFree (map->binning);
@@ -55,10 +60,27 @@
 
     psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
+    psImageBinningSetSkip (map->binning, map->field);
 
     return map;
 }
 
+bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff) {
+
+    assert (map);
+
+    map->binning->nXruff = nXruff;
+    map->binning->nYruff = nYruff;
+
+    psImageRecycle (map->map, nXruff, nYruff, PS_TYPE_F32);
+
+    psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
+    psImageBinningSetSkip (map->binning, map->field);
+
+    return true;
+}
+
 // generate a psImageMap (or NULL) with the given number of superpixels in X and Y
-bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, float badFrac) {
+// this function returns an error if the output map has impossible holes
+bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
 
     psImage *mask = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_MASK);
@@ -70,19 +92,23 @@
     // we can do this by accumulating a vector of pixel indexes for each cell
     psArray *pixelSets = psArrayAlloc (map->map->numCols*map->map->numRows);
-    for (int i = 0; i < pixelsSets->n; i++) {
+    for (int i = 0; i < pixelSets->n; i++) {
 	pixelSets->data[i] = psVectorAllocEmpty (4, PS_TYPE_S32);
     }
     // associate each value with a cell
     for (int i = 0; i < x->n; i++) {
-	// XXX use the psImageUnbin command to convert from x,y, to xRuff,yRuff?
-	int xRuff = x->data.F32[i] / map->binning->nXbin;
-	int yRuff = y->data.F32[i] / map->binning->nYbin;
+	int xRuff = psImageBinningGetRuffX (map->binning, x->data.F32[i]);
+	int yRuff = psImageBinningGetRuffY (map->binning, y->data.F32[i]);
 
 	int bin = xRuff + yRuff*map->map->numCols;
+	assert (bin >= 0);
+	assert (bin < pixelSets->n);
 	
 	psVector *pixels = pixelSets->data[bin];
 	pixels->data.S32[pixels->n] = i;
-	psVectorExtend (vector, 4, 1);
-    }
+	psVectorExtend (pixels, 4, 1);
+    }
+
+    // stats structure for getting the position centers
+    psStats *meanStat = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
 
     // accumulate the x,y coords for each point to calculate the mean position.
@@ -92,17 +118,28 @@
 	for (int ix = 0; ix < Nx; ix++) {
     
-	    // select the vector.  use psImageUnbin command to convert xRuff,yRuff to
-	    // xFine,yFine
+	    // pixel index for this cell
 	    psVector *pixels = pixelSets->data[ix + iy*Nx];
 
-	    psVector xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
-	    psVector yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
-	    psVector fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
-
+	    // storage vectors
+	    psVector *xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+	    psVector *yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+	    psVector *fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+
+	    // error vector, if needed
+	    psVector *dfCell = NULL;
+	    if (df) {
+		dfCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+	    }		
+
+	    // collect data for this cell
 	    for (int i = 0; i < pixels->n; i++) {
-		bin = pixels->data.S32[i];
-		xCell->data.F32[i] = x->data.F32[bin];
-		yCell->data.F32[i] = y->data.F32[bin];
-		fCell->data.F32[i] = f->data.F32[bin];
+		int bin = pixels->data.S32[i];
+		// convert x,y in the fine image to the ruff image
+		xCell->data.F32[i]  = psImageBinningGetRuffX (map->binning, x->data.F32[bin]);
+		yCell->data.F32[i]  = psImageBinningGetRuffY (map->binning, y->data.F32[bin]);
+		fCell->data.F32[i]  = f->data.F32[bin];
+		if (df) {
+		    dfCell->data.F32[i] = df->data.F32[bin];
+		}
 	    }
 
@@ -113,5 +150,5 @@
 	    // XXX need to supply a mask and skip the masked pixels when calculating the centroid
 	    // this will not in general be properly weighted...
-	    if (psVectorStats (map->stats, fCell, NULL, NULL, 0)) {
+	    if (psVectorStats (map->stats, fCell, dfCell, NULL, 0)) {
 		mask->data.U8[iy][ix] = 0;
 		// XXX ensure only one option is selected, or save both position and width
@@ -119,6 +156,8 @@
 
 		// calculate the mean position and save:
+		psStatsInit (meanStat);
 		psVectorStats (meanStat, xCell, NULL, NULL, 0);
 		xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); 
+		psStatsInit (meanStat);
 		psVectorStats (meanStat, yCell, NULL, NULL, 0);
 		yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); 
@@ -130,12 +169,24 @@
 	    psFree (yCell);
 	    psFree (fCell);
-	}
-    }
-
-    // at this point, for each map pixel, we have (f,x,y), or the pixel is masked.
-
+	    psFree (dfCell);
+	}
+    }
+    psFree (pixelSets);
+    psFree (meanStat);
+   // at this point, for each map pixel, we have (f,x,y), or the pixel is masked.
+
+    psFits *fits = NULL;
+
+    fits = psFitsOpen ("imageMap.raw.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
+    // did this analysis succeed?  (enough good or OK pixels?)
     psImage *state = psImagePixelInterpolateState (&map->nBad, &map->nPoor, mask, 0xff);
     map->nGood = mask->numCols * mask->numRows - map->nBad - map->nPoor;
     if (map->nBad > badFrac * mask->numCols * mask->numRows) {
+	psFree (xCoord);
+	psFree (yCoord);
+	psFree (mask);
 	return false;
     }
@@ -143,10 +194,89 @@
     // fit the valid pixels to (0,1,2) order polynomials, interpolate values to the pixel center
     // XXX I need to be careful about the pixel coordinates: center is 0,0 or 0.5, 0.5?
-    psImagePixelInterpolateCenters (map->map, xCoord, yCoord, state, mask, 0xff);
+    psImagePixelInterpolateCenter (map->map, xCoord, yCoord, state, mask, 0xff);
     psFree (xCoord);
     psFree (yCoord);
 
+    fits = psFitsOpen ("imageMap.ref.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
     psImagePixelInterpolatePoor (map->map, state, mask, 0xff);
 
+    fits = psFitsOpen ("imageMap.fix.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
+    psFree (state);
+    psFree (mask);
     return true;
 }
+
+// using the points given, generate a map with maximum resolution that yields only good and ok pixels
+bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
+    
+    int nXruff, nYruff;
+	
+    while (!psImageMapGenerate (map, x, y, f, df, badFrac)) {
+	// if we failed to build an acceptable map, decrease nXruff, nYruff as appropriate, and
+	// try again...  try to keep the aspect ratio.
+
+	// if both axes are at 1, give up
+	if ((map->binning->nXruff == 1) && (map->binning->nYruff == 1)) {
+	    return false;
+	}
+
+	// if one axis is at 1, decrement the other
+	if (map->binning->nXruff == 1) {
+	    nXruff = map->binning->nXruff;
+	    nYruff = map->binning->nYruff - 1;
+	    psImageMapModifyScale (map, nXruff, nYruff);
+	    continue;
+	}
+	if (map->binning->nYruff == 1) {
+	    nYruff = map->binning->nYruff;
+	    nXruff = map->binning->nXruff - 1;
+	    psImageMapModifyScale (map, nXruff, nYruff);
+	    continue;
+	}
+
+	// otherwise, decrement the larger axis, and set the smaller based
+	// on the aspect ratio
+	float aRatio = map->binning->nXruff / map->binning->nYruff;
+	if (map->binning->nXruff > map->binning->nYruff) {
+	    nXruff = map->binning->nXruff - 1;
+	    nYruff = (int)(0.5 + (nXruff / aRatio));
+	} else {
+	    nYruff = map->binning->nYruff - 1;
+	    nXruff = (int)(0.5 + (nYruff * aRatio));
+	}
+
+	psImageMapModifyScale (map, nXruff, nYruff);
+    }
+    return true;
+}
+
+double psImageMapEval (psImageMap *map, float x, float y) {
+
+    double result;
+
+    result = psImageUnbinPixel(x, y, map->map, map->binning);
+
+    return result; 
+}
+
+psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y) {
+
+    assert (x);
+    assert (y);
+    assert (x->n == y->n);
+    assert (map);
+
+    psVector *result = psVectorAlloc (x->n, PS_TYPE_F32);
+
+    for (int i = 0; i < x->n; i++) {
+	result->data.F32[i] = psImageUnbinPixel(x->data.F32[i], y->data.F32[i], map->map, map->binning);
+    }
+
+    return result; 
+}
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h	(revision 14773)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h	(revision 14774)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-02 02:03:58 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:19:35 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -35,6 +35,16 @@
 psImageMap *psImageMapAlloc(psImage *field, psImageBinning *binning, psStats *stats) PS_ATTR_MALLOC;
 
+bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff);
+
 // generate a psImageMap (or NULL) with the given number of superpixels in X and Y
-bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, float badFrac);
+bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
+
+bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
+
+// apply the psImageMap to the given coordinate (fine image pixels)
+double psImageMapEval (psImageMap *map, float x, float y);
+
+// apply the psImageMap to the given coordinate vectors (fine image pixels)
+psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y);
 
 /// @}