Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 21060)
+++ /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 21061)
@@ -30,4 +30,5 @@
 $(SRC)/flux.$(ARCH).o		   \
 $(SRC)/fixwrap.$(ARCH).o	   \
+$(SRC)/fixcols.$(ARCH).o	   \
 $(SRC)/gauss.$(ARCH).o		   \
 $(SRC)/getvel.$(ARCH).o	   \
Index: /trunk/Ohana/src/opihi/cmd.astro/fixcols.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/fixcols.c	(revision 21061)
+++ /trunk/Ohana/src/opihi/cmd.astro/fixcols.c	(revision 21061)
@@ -0,0 +1,246 @@
+# include "astro.h"
+
+// this function is a hack to mask bad columns / rows in an image.  we assume that there
+// are columns or rows which are, as a whole, significant deviant from the others.
+// in fact, we measure this by segments
+
+// SIGMA is a robust measure of the sigma, STDEV is the formal stdev of the values in range
+int vector_stats (float *vect, int Nvect, float *MEDIAN, float *SIGMA, float *STDEV);
+
+int fixrows (int argc, char **argv) {
+  
+  int ix, iy, Nx, Ny, Nvect, start, stop;
+  float *Vin, *vect, *meds, median, stdev, sigma, Nsigma, value;
+  Buffer *in;
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: fixrows (buffer) (col_start) (col_stop) (Nsigma) (value)\n");
+    return (FALSE);
+  }
+
+  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  start = atoi (argv[2]);
+  stop = atoi (argv[3]);
+  Nsigma = atof (argv[4]);
+  value = atof (argv[5]);
+
+  Nx = in[0].matrix.Naxis[0];
+  Ny = in[0].matrix.Naxis[1];
+
+  ALLOCATE (meds, float, Ny);
+  ALLOCATE (vect, float, stop - start);
+
+  Vin  = (float *)in[0].matrix.buffer;
+
+  // measure the median for each column segment, accumulate the stats
+  for (iy = 0; iy < Ny; iy++) {
+
+    Nvect = 0;
+    for (ix = start; ix < stop; ix++) {
+      vect[Nvect] = Vin[ix + iy*Nx];
+      Nvect ++;
+    }
+
+    // get the median, robust sigma, and formal stdev of the column segment
+    vector_stats (vect, Nvect, &median, &sigma, &stdev);
+    meds[iy] = median;
+    // fprintf (stderr, "median : %f +/- %f (%f)\n", median, sigma, stdev);
+  }
+
+  // get the median, robust sigma, and formal stdev of the medians
+  vector_stats (meds, Nx, &median, &sigma, &stdev);
+  // fprintf (stderr, "median of medians : %f +/- %f (%f)\n", median, sigma, stdev);
+
+  // mask values which are more than Nsigma out of range
+  for (iy = 0; iy < Ny; iy++) {
+    if (fabs(meds[iy] - median) < Nsigma * sigma) continue;
+    for (ix = start; ix < stop; ix++) {
+      Vin[ix + iy*Nx] = value;
+    }
+  }
+  
+  free (meds);
+  free (vect);
+
+  return (TRUE);
+}
+
+int fixcols (int argc, char **argv) {
+  
+  int ix, iy, Nx, Ny, Nvect, start, stop;
+  float *Vin, *vect, *meds, median, stdev, sigma, Nsigma, value;
+  Buffer *in;
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: fixcols (buffer) (row_start) (row_stop) (Nsigma) (value)\n");
+    return (FALSE);
+  }
+
+  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  start = atoi (argv[2]);
+  stop = atoi (argv[3]);
+  Nsigma = atof (argv[4]);
+  value = atof (argv[5]);
+
+  Nx = in[0].matrix.Naxis[0];
+  Ny = in[0].matrix.Naxis[1];
+
+  ALLOCATE (meds, float, Nx);
+  ALLOCATE (vect, float, stop - start);
+
+  Vin  = (float *)in[0].matrix.buffer;
+
+  // measure the median for each column segment, accumulate the stats
+  for (ix = 0; ix < Nx; ix++) {
+
+    Nvect = 0;
+    for (iy = start; iy < stop; iy++) {
+      vect[Nvect] = Vin[ix + iy*Nx];
+      Nvect ++;
+    }
+
+    // get the median, robust sigma, and formal stdev of the column segment
+    vector_stats (vect, Nvect, &median, &sigma, &stdev);
+    meds[ix] = median;
+    // fprintf (stderr, "median : %f +/- %f (%f)\n", median, sigma, stdev);
+  }
+
+  // get the median, robust sigma, and formal stdev of the medians
+  vector_stats (meds, Nx, &median, &sigma, &stdev);
+  // fprintf (stderr, "median of medians : %f +/- %f (%f)\n", median, sigma, stdev);
+
+  // mask values which are more than Nsigma out of range
+  for (ix = 0; ix < Nx; ix++) {
+    if (fabs(meds[ix] - median) < Nsigma * sigma) continue;
+    for (iy = start; iy < stop; iy++) {
+      Vin[ix + iy*Nx] = value;
+    }
+  }
+  
+  free (meds);
+  free (vect);
+
+  return (TRUE);
+}
+
+// this macro is accumulating the histogram until the sum passes the desired threshold
+// point the fractional position of that threshold is found by interpolation
+# define SUMPOINT(POINT, NPOINT, FRACTION)				\
+  if ((NPOINT) == -1) {							\
+    (POINT) += Nval[i];							\
+    if ((POINT) >= FRACTION*N) {					\
+      (NPOINT) = i;							\
+      value = i * dx + min;						\
+      if (Nval[i] != 0) {						\
+	/* POINT for the bin center is not exactly FRACTION*N */	\
+	/* interpolate to the true value, but limit the scale of the */	\
+	/* interpolation to a single bin */				\
+	adjust = dx * (FRACTION*N - (POINT)) / Nval[i];			\
+	value += MAX (MIN (adjust, dx), -dx);				\
+      }									\
+      (POINT) = value;							\
+    }									\
+  }
+
+// for testing, this is a useful line:
+// fprintf (stderr, "for point %f (%f): bin: %d = %f, dx: %f, adjust: %f, value: %f, POINT: %f, Nval[i]: %d, Nval[i-1]: %d\n",
+// FRACTION, FRACTION*N, i, i*dx + min, dx, adjust, value, POINT, Nval[i], Nval[i-1]);
+
+int vector_stats (float *vect, int Nvect, float *MEDIAN, float *SIGMA, float *STDEV) {
+  
+  int i, iter, N, done;
+  double max, min, sum, var, dvar, mean, stdev, sigma, adjust, value;
+  int Nval[1001], bin, Nmode, N_MED, N_LQ, N_UQ, N_SM, N_SP;
+  double dx, mode, MED, UQ, LQ, SM, SP;
+  float *X;
+
+  /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */
+
+  /* calculate max, min, mean, sum, npix */
+  max = -HUGE_VAL;
+  min = HUGE_VAL;
+  sum = N = 0;
+  X = vect;
+  for (i = 0; i < Nvect; i++, X++) {
+    if (!finite (*X)) continue;
+    max = MAX (*X, max);
+    min = MIN (*X, min);
+    sum += *X;
+    N++;
+  }      
+  mean = sum / N;
+
+  // calculate median and mode with initial resolution dx = (max - min) / 1000 
+  // we adjust the resolution and retry if the resulting robust sigma < dx
+  dx = (max - min) / 1000;
+  if (dx == 0) {
+    MED = mode = min;
+    stdev = 0.0;
+    sigma = 0.0;
+    goto skip;
+  }
+
+  done = FALSE;
+  for (iter = 0; (iter < 3) && !done; iter ++) {
+
+    bzero (Nval, 1001*sizeof(int));
+    var = 0;
+    X = vect;
+    for (i = 0; i < Nvect; i++, X++) {
+      if (!finite (*X)) continue;
+      bin = MAX (0, MIN (1000, (*X - min) / dx));
+      Nval[bin] ++;
+      if (bin == 0) continue;
+      if (bin == 1000) continue;
+      // only count the values in range for the stdev
+      dvar = (*X - mean);
+      var += dvar*dvar;
+    }      
+    stdev = sqrt (var / N);
+
+    Nmode = 0;
+    mode = Nval[Nmode];
+    MED = SM = SP = LQ = UQ = 0;
+    N_MED = N_SM = N_SP = N_LQ = N_UQ = -1;
+    for (i = 0; i < 1001; i++) {
+      SUMPOINT (MED, N_MED, 0.5000);
+      SUMPOINT (LQ,  N_LQ,  0.2500);
+      SUMPOINT (UQ,  N_UQ,  0.7500);
+      SUMPOINT (SM,  N_SM,  0.1587);
+      SUMPOINT (SP,  N_SP,  0.8413);
+
+      if (mode < Nval[i]) {
+	Nmode = i;
+	mode = Nval[Nmode];
+      }
+    }
+
+    sigma = 0.5*(SP - SM);
+    mode = Nmode * dx + min;
+
+    // if (sigma < dx), we adjust the range to be centered on the median
+    if (sigma < dx) {
+      dx *= 0.1;
+      min = MED - 500*dx;
+    } else {
+      done = TRUE;
+    }
+  }
+
+skip:
+
+  *MEDIAN = MED;
+  *SIGMA = sigma;
+  *STDEV = stdev;
+
+  // set_variable ("MIN",      min);
+  // set_variable ("MAX",      max);
+  // set_variable ("MEAN",     mean);
+  // set_variable ("MODE",     mode);
+  // set_variable ("TOTAL",    sum);
+  // set_int_variable ("NPIX", N);
+  // set_variable ("SIGMA",    stdev);
+
+  return (TRUE);
+}
+
Index: /trunk/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 21060)
+++ /trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 21061)
@@ -13,4 +13,6 @@
 int flux                    PROTO((int, char **));
 int fixwrap                 PROTO((int, char **));
+int fixcols                 PROTO((int, char **));
+int fixrows                 PROTO((int, char **));
 int gauss                   PROTO((int, char **));
 int gaussfit                PROTO((int, char **));
@@ -50,4 +52,6 @@
   {1, "flux",        flux,         "flux in a convex contour"},
   {1, "fixwrap",     fixwrap,      "fix megacam over-wrapped pixels"},
+  {1, "fixcols",     fixcols,      "fix bad columns by comparing with others"},
+  {1, "fixrows",     fixrows,      "fix bad rows by comparing with others"},
   {1, "gauss",       gauss,        "get statistics on a star, assuming gaussian profile"},
   {1, "getvel",      getvel,       "rotcurve to velocities"},