Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c	(revision 37690)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c	(revision 37691)
@@ -5,4 +5,6 @@
 
 int dump_map_data (float *x, float *y, float *f, int Npts, char *filename);
+int AstromOffsetMapFit_Chisq (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir);
+int AstromOffsetMapFit_Median (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir);
 
 float AstromOffsetMapValue (AstromOffsetMap *map, float x, float y, int xdir) {
@@ -63,7 +65,13 @@
 }
 
+int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) {
+
+  int status = AstromOffsetMapFit_Median (map, x, y, f, Npts, xdir);
+  return status;
+}
+
 // given (x,y),value vector sets, choose the map that minimizes the difference between 
 // the values and bilinear interpolation of the map
-int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) {
+int AstromOffsetMapFit_Chisq (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) {
 
   int i, j, ix, iy, jx, jy;
@@ -361,4 +369,74 @@
 }
 
+// given (x,y),value vector sets, choose the map that minimizes the difference between 
+// the values and bilinear interpolation of the map
+int AstromOffsetMapFit_Median (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) {
+
+  int i, ix, iy;
+
+  // choose to map direction:
+  float **V = xdir ? map->dXv : map->dYv;
+
+  // measure clipped median in each bin
+  VStatsType stats;
+  stats.statmode = VSTATS_INNER_MEAN;
+
+  int Nx = map->Nx;
+  int Ny = map->Ny;
+  int Npix = Nx*Ny;
+
+  double **values;
+  int    *Nvalue;
+  ALLOCATE (Nvalue, int, Npix);
+  ALLOCATE (values, double *, Npix);
+  for (i = 0; i < Npix; i++) {
+    ALLOCATE (values[i], double, Npts);
+    Nvalue[i] = 0;
+  }
+
+  // assign the points to the map cells
+  for (i = 0; i < Npts; i++) {
+
+    // data value & weight for this point
+    if (!isfinite(f[i])) continue;
+
+    // if (mask && (mask[i] & maskValue)) continue;
+    
+    // base coordinate offset for this point (x,y) relative to this map element (n,m)
+    // double dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5);
+    // double dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
+
+    // bin for this point
+    ix = MAX(0, MIN(Nx, (int)(x[i] * map->dX)));
+    iy = MAX(0, MIN(Ny, (int)(y[i] * map->dY)));
+    int I = ix + Nx * iy;
+
+    int N = Nvalue[I];
+    values[I][N] = f[i];
+    Nvalue[I] ++;
+  }    
+
+  for (ix = 0; ix < Nx; ix++) {
+    for (iy = 0; iy < Ny; iy++) {
+      int I = ix + Nx * iy;
+      if (Nvalue[I] > 5) {
+	vstats_getstats (values[I], NULL, NULL, Nvalue[I], &stats);
+	V[ix][iy] = stats.mean;
+      } else {
+	V[ix][iy] = NAN;
+      }
+    }
+  }
+
+  // XXX free someone?
+  for (i = 0; i < Npix; i++) {
+    free (values[i]);
+  }    
+  free (values);
+  free (Nvalue);
+
+  return TRUE;
+}
+
 // this function repairs an image with NAN pixels (only valid for a small-scale map -- no robust mean)
 int AstromOffsetMapRepair (AstromOffsetMap *map, int xdir) {