Index: trunk/psLib/src/imageops/psImageMapFit.c
===================================================================
--- trunk/psLib/src/imageops/psImageMapFit.c	(revision 24091)
+++ trunk/psLib/src/imageops/psImageMapFit.c	(revision 25753)
@@ -33,4 +33,5 @@
 #include "psStats.h"
 #include "psImageBinning.h"
+#include "psImageStructManip.h"
 #include "psImageMap.h"
 // #include "psImagePixelInterpolate.h"
@@ -752,2 +753,65 @@
 }
 
+// this function repairs an image with NAN pixels (only valid for a small-scale map -- no robust mean)
+bool psImageMapRepair (psImage *image) {
+
+    // we are going to repair the image by:
+    // 1) finding NAN pixels
+    // 2) if any of the neighbors are valid,
+    //    replace with the mean of the neighbors
+    // 3) otherwise, replace with the image mean
+
+    // copy the image so the repaired pixels do not affect the input
+    psImage *fixed = psImageCopy (NULL, image, PS_TYPE_F32);
+
+    // find the global mean
+    float mean = 0.0;
+    float npix = 0.0;
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    if (!isfinite(image->data.F32[iy][ix])) continue;
+	    mean += image->data.F32[iy][ix];
+	    npix += 1.0;
+	}
+    }
+    mean /= npix;
+
+    // find the NAN pixels:
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    if (isfinite(image->data.F32[iy][ix])) {
+		fixed->data.F32[iy][ix] = image->data.F32[iy][ix];
+		continue;
+	    }
+
+	    // find mean of all possible neighbors
+	    float meanLocal = 0.0;
+	    float npixLocal = 0.0;
+	    for (int jy = -1; jy <= +1; jy++) {
+		int ny = iy + jy;
+		if (ny < 0) continue;
+		if (ny >= image->numRows) continue;
+		for (int jx = -1; jx <= +1; jx++) {
+		    int nx = ix + jx;
+		    if (nx < 0) continue;
+		    if (nx >= image->numCols) continue;
+		    if (!isfinite(image->data.F32[ny][nx])) continue;
+		    meanLocal += image->data.F32[ny][nx];
+		    npixLocal += 1.0;
+		}
+	    }
+	    meanLocal = (npixLocal > 0.0) ? meanLocal / npixLocal : mean;
+	    fixed->data.F32[iy][ix] = meanLocal;
+	}
+    }
+    
+    // find the NAN pixels:
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    image->data.F32[iy][ix] = fixed->data.F32[iy][ix];
+	}
+    }
+    psFree (fixed);
+
+    return true;
+}