Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile	(revision 40757)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile	(revision 40758)
@@ -95,4 +95,5 @@
 $(SRC)/mcreate.$(ARCH).o	\
 $(SRC)/medacc.$(ARCH).o	\
+$(SRC)/mgaussdev.$(ARCH).o	\
 $(SRC)/mget.$(ARCH).o		\
 $(SRC)/mget3d.$(ARCH).o		\
@@ -101,4 +102,5 @@
 $(SRC)/medimage.$(ARCH).o	\
 $(SRC)/medimage_commands.$(ARCH).o \
+$(SRC)/medimage_calc.$(ARCH).o \
 $(SRC)/mset.$(ARCH).o		\
 $(SRC)/needles.$(ARCH).o	\
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c	(revision 40757)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c	(revision 40758)
@@ -85,4 +85,5 @@
 int mcreate          PROTO((int, char **));
 int medacc           PROTO((int, char **));
+int mgaussdev        PROTO((int, char **));
 int mget             PROTO((int, char **));
 int mget3d           PROTO((int, char **));
@@ -273,4 +274,5 @@
   {1, "imcreate",     mcreate,          "create an image"},
   {1, "medacc",       medacc,           "accumulate vector values in another vector"},
+  {1, "mgaussdev",    mgaussdev,        "generate a gaussian deviate image"},
   {1, "mget",         mget,             "extract a vector from an image"},
   {1, "mget3d",       mget3d,           "extract a vector from a 3D image"},
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_calc.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_calc.c	(revision 40758)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_calc.c	(revision 40758)
@@ -0,0 +1,365 @@
+# include "data.h"
+
+float weight_cauchy_square_flt (float x2);
+float irls_mean (float *val, float *wgt, int N, float *outvar, int *Npt);
+float irls_fraction_interpolate (float *values, float fraction, int Npts);
+float irls_robust_stdev (float *values, int Npts);
+float irls_wtmean (float *val, float *wgt, int *idx, int Npt, float *variance);
+void irls_bootstrap (int *idx, int Npts);
+void irls_init (int N);
+void irls_free (void);
+
+# define IRLS_TOLERANCE 1e-4
+int BOOTSTRAP = FALSE;
+int BOOTSTRAP_NITER = 100;
+
+static float *irls_valsub  = NULL;
+static float *irls_wgtsub  = NULL;
+static int   *irls_idx     = NULL;
+static float *irls_testval = NULL;
+
+enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN};
+
+int medimage_calc (int argc, char **argv) {
+
+  int ix, iy, n, N;
+  Buffer *output;
+  
+  BOOTSTRAP = FALSE;
+  if ((N = get_argument (argc, argv, "-bootstrap"))) {
+    BOOTSTRAP = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  BOOTSTRAP_NITER = 100;
+  if ((N = get_argument (argc, argv, "-bootstrap-iter"))) {
+    remove_argument (N, &argc, argv);
+    BOOTSTRAP_NITER = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+    BOOTSTRAP = TRUE;
+  }
+
+  int mode = CALC_MEDIAN;
+  if ((N = get_argument (argc, argv, "-mean"))) {
+    mode = CALC_MEAN;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-irls"))) {
+    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; }
+    mode = CALC_IRLS;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-wtmean"))) {
+    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; }
+    mode = CALC_WTMEAN;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-median"))) {
+    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; }
+    mode = CALC_MEDIAN;
+    remove_argument (N, &argc, argv);
+  }
+
+  Buffer *variance = NULL;
+  if ((N = get_argument (argc, argv, "-variance"))) {
+    remove_argument (N, &argc, argv);
+    if ((variance = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean,-irls,-wtmean,-median] [-variance output] [-bootstrap] [-bootstrap-iter (N)]\n");
+    gprint (GP_ERR, "       calculate the median image for the median image set\n");
+    gprint (GP_ERR, "  This function calculates the 'average' image for the supplied set of images.\n");
+    gprint (GP_ERR, "  Four options are available to calculate the average:\n");
+    gprint (GP_ERR, "   -mean (straight arithmetic mean)\n");
+    gprint (GP_ERR, "   -median\n");
+    gprint (GP_ERR, "   -wtmean (arithmetic mean, weighted by supplied variances)\n");
+    gprint (GP_ERR, "   -irls (iteratively-reweighted least-squares)\n");
+    gprint (GP_ERR, "  If -variance is supplied the returned image represents the estimate of the variance on the average image\n");
+    gprint (GP_ERR, "  For the four average methods above, the sqrt of the variance has the following meanings:\n");
+    gprint (GP_ERR, "   -mean   : standard deviation / sqrt(Npts)\n");
+    gprint (GP_ERR, "   -median : standard deviation / sqrt(Npts)\n");
+    gprint (GP_ERR, "   -wtmean : formal error on the weighted mean [sum (1 / variances)]\n");
+    gprint (GP_ERR, "   -irls   : formal error on the pixels used (excluding ~3 sigma outliers)\n");
+    gprint (GP_ERR, "   note that the irls formal error slightly underestimates the observed standard deviation\n");
+    return FALSE;
+  }
+
+  MedImageType *median = FindMedImage (argv[1]);
+  if (!median) {
+    gprint (GP_ERR, "median image %s not found\n", argv[1]);
+    return FALSE;
+  }
+
+  if ((output = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
+  int Ninput = median->Ninput;
+  int Nx = median->Nx;
+  int Ny = median->Ny;
+
+  ALLOCATE_PTR (val, float, Ninput);
+  ALLOCATE_PTR (wgt, float, Ninput);
+
+  if (mode == CALC_IRLS) irls_init (Ninput);
+
+  gfits_free_matrix (&output->matrix);
+  gfits_free_header (&output->header);
+  if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
+
+  float *outvalue = (float *) output->matrix.buffer;
+  float *varvalue = NULL;
+  if (variance) {
+    gfits_free_matrix (&variance->matrix);
+    gfits_free_header (&variance->header);
+    if (!CreateBuffer (variance, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
+
+    varvalue = (float *) variance->matrix.buffer;
+  }
+
+  // save the number of points per pixel
+  Buffer *nptbuf = SelectBuffer ("irls_npt", ANYBUFFER, TRUE);
+  gfits_free_matrix (&nptbuf->matrix);
+  gfits_free_header (&nptbuf->header);
+  if (!CreateBuffer (nptbuf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
+  float *NptVal = (float *) nptbuf->matrix.buffer;
+
+  for (iy = 0; iy < Ny; iy++) {
+    for (ix = 0; ix < Nx; ix++) {
+
+      int N = 0;
+      int Npix = ix + Nx*iy;
+      for (n = 0; n < Ninput; n++) {
+	float v = median->flx[n][Npix];
+	if (!isfinite(v)) continue;
+	val[N] = v;
+	wgt[N] = 1.0;
+	if (median->var[n]) {
+	  float s = median->var[n][Npix];
+	  if (!isfinite(s)) continue;
+	  if (fabs(s) < 2*FLT_MIN) s = 2*FLT_MIN;
+	  wgt[N] = 1.0 / s;
+	}
+	N++;
+      }
+      if (N == 0) continue;
+
+      switch (mode) {
+	case CALC_MEDIAN:
+	  fsort (val, N);
+	  outvalue[Npix] = val[(int)(0.5*N)];
+	  if (varvalue) {
+	    float sum = 0.0;
+	    for (n = 0; n < N; n++) {
+	      sum += SQ(val[n] - outvalue[Npix]);
+	    }
+	    // variance on the mean (stdev / sqrt(N))^2
+	    varvalue[Npix] = sum / (N - 1) / N;
+	  }
+	  break;
+	case CALC_MEAN: {
+	  float sum = 0.0;
+	  for (n = 0; n < N; n++) {
+	    sum += val[n];
+	  }
+	  outvalue[Npix] = sum / (float) N;
+
+	  if (varvalue) {
+	    float Sum = 0.0;
+	    for (n = 0; n < N; n++) {
+	      Sum += SQ(val[n] - outvalue[Npix]);
+	    }
+	    // variance on the mean (stdev / sqrt(N))^2
+	    varvalue[Npix] = Sum / (N - 1) / N;
+	  }
+	  break;
+	}
+	case CALC_WTMEAN: {
+	  float variance;
+	  outvalue[Npix] = irls_wtmean (val, wgt, NULL, N, &variance);
+	  if (varvalue) { varvalue[Npix] = variance; }
+	  break;
+	}
+	case CALC_IRLS: 
+	  if (varvalue) {
+	    int Npts = 0;
+	    outvalue[Npix] = irls_mean (val, wgt, N, &varvalue[Npix], &Npts);
+	    NptVal[Npix] = Npts;
+	  } else {
+	    int Npts = 0;
+	    outvalue[Npix] = irls_mean (val, wgt, N, NULL, &Npts);
+	    NptVal[Npix] = Npts;
+	  }
+      }
+    }
+  }
+  irls_free();
+  return TRUE;
+}
+
+// allocate temporary vectors
+void irls_init (int N) {
+
+  ALLOCATE (irls_valsub, float, N);
+  ALLOCATE (irls_wgtsub, float, N);
+  ALLOCATE (irls_idx, int, N);
+  ALLOCATE (irls_testval, float, BOOTSTRAP_NITER);
+}
+
+// free temporary vectors
+void irls_free (void) {
+  FREE (irls_valsub); 
+  FREE (irls_wgtsub);
+  FREE (irls_idx);
+  FREE (irls_testval);
+  irls_valsub = NULL;
+  irls_wgtsub = NULL;
+  irls_idx = NULL;
+  irls_testval = NULL;
+}
+
+float irls_mean (float *val, float *wgt, int N, float *outvar, int *Npt) {
+
+  // calculate weighted mean
+  float Value = irls_wtmean (val, wgt, NULL, N, NULL);
+
+  int converged = FALSE;
+  for (int i = 0; (i < 10) && !converged; i++) {
+    float ValueLast = Value;
+    float S1 = 0.0, S2 = 0.0;
+    
+    // calculate weight modification based on distances (squared).
+    // use modifier to calculate new weighted mean
+    for (int n = 0; n < N; n++) {
+      float dV = (val[n] - Value);
+      float d2 = SQ(dV) * wgt[n];
+      
+      float Mod = weight_cauchy_square_flt (d2);
+      S1 += Mod * wgt[n] * val[n];
+      S2 += Mod * wgt[n];
+    }
+    Value = S1 / S2;
+
+    float delta = fabs(Value - ValueLast);
+    if (delta < Value * IRLS_TOLERANCE) converged = TRUE;
+  }
+
+  if (outvar) {
+    if (BOOTSTRAP) {
+      // generate a subset vector of just the accepted points
+    
+      // calculate stdev of high-weight points
+      int npt = 0;
+      for (int n = 0; n < N; n++) {
+
+	float dV = (val[n] - Value);
+	float d2 = SQ(dV) * wgt[n];
+      
+	float Mod = weight_cauchy_square_flt (d2);
+	if (Mod < 0.1) continue; // totally ad-hoc number
+
+	irls_valsub[npt] = val[n];
+	irls_wgtsub[npt] = wgt[n];
+	npt ++;
+      }
+
+      for (int iter = 0; iter < BOOTSTRAP_NITER; iter++) {
+	irls_bootstrap (irls_idx, npt); // fill idx with the index of the resampled points
+	irls_testval[iter] = irls_wtmean (irls_valsub, irls_wgtsub, irls_idx, npt, NULL);
+      }
+
+      float sigma = irls_robust_stdev (irls_testval, BOOTSTRAP_NITER);
+
+      *outvar = SQ(sigma);
+      *Npt = npt;
+    } else {
+      // calculate stdev of high-weight points
+      float S1 = 0.0, S2 = 0.0;
+      int npt = 0;
+      for (int n = 0; n < N; n++) {
+
+	float dV = (val[n] - Value);
+	float d2 = SQ(dV) * wgt[n];
+      
+	float Mod = weight_cauchy_square_flt (d2);
+	if (Mod < 0.1) continue; // totally ad-hoc number
+
+	S1 += SQ(dV);
+	S2 += wgt[n];
+	npt ++;
+      }
+      *Npt = npt;
+      *outvar = 1 / S2;
+    }
+  }
+
+  return Value;
+}
+
+void irls_bootstrap (int *idx, int Npts) {
+  // generate an index (idx) containing the index values
+  // for the Npts randomly resampled points
+  for (int i = 0; i < Npts; i++) {
+    idx[i] = Npts * drand48();
+  }
+}
+
+float irls_wtmean (float *val, float *wgt, int *idx, int Npt, float *variance) {
+  float S1 = 0.0, S2 = 0.0;
+  for (int n = 0; n < Npt; n++) {
+    int N = idx ? idx[n] : n;
+    S1 += wgt[N] * val[N];
+    S2 += wgt[N];
+  }
+  float Value = S1 / S2;
+  if (variance) *variance = 1.0 / S2;
+  return Value;
+}
+
+float irls_robust_stdev (float *values, int Npts) {
+
+  fsort (values, Npts);
+
+  float Slo = irls_fraction_interpolate (values, 0.158655, Npts);
+  float Shi = irls_fraction_interpolate (values, 0.841345, Npts);
+  float sigma = (Shi - Slo) / 2.0;
+  
+  return sigma;
+}
+
+float irls_fraction_interpolate (float *values, float fraction, int Npts) {
+
+  float F = fraction * Npts;
+  int   N = fraction * Npts;
+
+  if (N < 0        ) return NAN;
+  if (N >= Npts - 2) return NAN;
+
+  // interpolate between N,N+1
+    
+  float S = (F - N) * (values[N+1] - values[N]) + values[N];
+  return S;
+}
+
+
+// exp(-(x^2/s^2)/2) = (1/2)
+//     -(x^2/s^2)/2  = ln(1/2)
+//      (x^2/s^2)/2  = ln(2)
+//      (x^2/s^2)    = 2ln(2)
+//      (x  /s)      = sqrt(2ln(2)) : half-width at half-max
+//       FWHM        = 2sqrt(2ln(2))
+
+// R2 = (X / 2.385)^2 = (X^2 / 2.385^2)
+
+# define CAUCY_FACTOR 2.0
+
+/*
+float weight_cauchy_square_flt (float x2) {
+  return (1.0);
+}
+*/
+
+float weight_cauchy_square_flt (float x2) {
+  float r2 = x2 / CAUCY_FACTOR;
+  return (1.0 / (1.0 + r2));
+}
+
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c	(revision 40757)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c	(revision 40758)
@@ -1,8 +1,3 @@
 # include "data.h"
-
-float weight_cauchy_square_flt (float x2);
-float irls_mean (float *val, float *wgt, int N);
-
-# define IRLS_TOLERANCE 1e-4
 
 int medimage_list (int argc, char **argv) {
@@ -85,161 +80,4 @@
   return TRUE;
 }
-
-enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN};
-int medimage_calc (int argc, char **argv) {
-
-  int ix, iy, n, N;
-  Buffer *output;
-  
-  int mode = CALC_MEDIAN;
-  if ((N = get_argument (argc, argv, "-mean"))) {
-    mode = CALC_MEAN;
-    remove_argument (N, &argc, argv);
-  }
-  if ((N = get_argument (argc, argv, "-irls"))) {
-    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -irls, -wtmean\n"); return FALSE; }
-    mode = CALC_IRLS;
-    remove_argument (N, &argc, argv);
-  }
-  if ((N = get_argument (argc, argv, "-wtmean"))) {
-    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -irls, -wtmean\n"); return FALSE; }
-    mode = CALC_WTMEAN;
-    remove_argument (N, &argc, argv);
-  }
-
-  Buffer *variance = NULL;
-  if ((N = get_argument (argc, argv, "-variance"))) {
-    remove_argument (N, &argc, argv);
-    if ((variance = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
-    remove_argument (N, &argc, argv);
-  }
-
-  if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean,-irls,-wtmean]\n");
-    gprint (GP_ERR, "       calculate the median image for the median image set\n");
-    return FALSE;
-  }
-
-  MedImageType *median = FindMedImage (argv[1]);
-  if (!median) {
-    gprint (GP_ERR, "median image %s not found\n", argv[1]);
-    return FALSE;
-  }
-
-  if ((output = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
-
-  int Ninput = median->Ninput;
-  int Nx = median->Nx;
-  int Ny = median->Ny;
-
-  ALLOCATE_PTR (val, float, Ninput);
-  ALLOCATE_PTR (wgt, float, Ninput);
-
-  gfits_free_matrix (&output->matrix);
-  gfits_free_header (&output->header);
-  if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
-
-  float *outvalue = (float *) output->matrix.buffer;
-
-  for (iy = 0; iy < Ny; iy++) {
-    for (ix = 0; ix < Nx; ix++) {
-
-      int N = 0;
-      int Npix = ix + Nx*iy;
-      for (n = 0; n < Ninput; n++) {
-	float v = median->flx[n][Npix];
-	if (!isfinite(v)) continue;
-	val[N] = v;
-	wgt[N] = 1.0;
-	if (median->var[n]) {
-	  float s = median->var[n][Npix];
-	  if (!isfinite(s)) continue;
-	  if (fabs(s) < 2*FLT_MIN) s = 2*FLT_MIN;
-	  wgt[N] = 1.0 / s;
-	}
-	N++;
-      }
-      if (N == 0) continue;
-
-      switch (mode) {
-	case CALC_MEDIAN:
-	  fsort (val, N);
-	  outvalue[Npix] = val[(int)(0.5*N)];
-	  break;
-	case CALC_MEAN: {
-	  float sum = 0.0;
-	  for (n = 0; n < N; n++) {
-	    sum += val[n];
-	  }
-	  outvalue[Npix] = sum / (float) N;
-	  break;
-	}
-	case CALC_WTMEAN: {
-	  float S1 = 0.0, S2 = 0.0;
-	  for (n = 0; n < N; n++) {
-	    S1 += wgt[n] * val[n];
-	    S2 += wgt[n];
-	  }
-	  outvalue[Npix] = S1 / S2;
-	  break;
-	}
-	case CALC_IRLS: 
-	  outvalue[Npix] = irls_mean (val, wgt, N);
-      }
-    }
-  }
-  return TRUE;
-}
-
-float irls_mean (float *val, float *wgt, int N) {
-
-  // calculate weighted mean
-  float S1 = 0.0, S2 = 0.0;
-  for (int n = 0; n < N; n++) {
-    S1 += wgt[n] * val[n];
-    S2 += wgt[n];
-  }
-  float Value = S1 / S2;
-  
-  int converged = FALSE;
-  for (int i = 0; (i < 10) && !converged; i++) {
-    float ValueLast = Value;
-
-    float S1 = 0.0, S2 = 0.0;
-    
-    // calculate weight modification based on distances (squared).
-    // use modifier to calculate new weighted mean
-    for (int n = 0; n < N; n++) {
-      float dV = (val[n] - Value);
-      float d2 = SQ(dV) * wgt[n];
-      
-      float Mod = weight_cauchy_square_flt (d2);
-      S1 += Mod * wgt[n] * val[n];
-      S2 += Mod * wgt[n];
-    }
-    Value = S1 / S2;
-
-    float delta = fabs(Value - ValueLast);
-    if (delta < Value * IRLS_TOLERANCE) converged = TRUE;
-  }
-  return Value;
-}
-
-// exp(-(x^2/s^2)/2) = (1/2)
-//     -(x^2/s^2)/2  = ln(1/2)
-//      (x^2/s^2)/2  = ln(2)
-//      (x^2/s^2)    = 2ln(2)
-//      (x  /s)      = sqrt(2ln(2)) : half-width at half-max
-//       FWHM        = 2sqrt(2ln(2))
-
-// R2 = (X / 2.385)^2 = (X^2 / 2.385^2)
-
-# define CAUCY_FACTOR 1.0
-
-float weight_cauchy_square_flt (float x2) {
-  float r2 = x2 / CAUCY_FACTOR;
-  return (1.0 / (1.0 + r2));
-}
-
 
 /* 
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/mgaussdev.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/mgaussdev.c	(revision 40758)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/mgaussdev.c	(revision 40758)
@@ -0,0 +1,34 @@
+# include "data.h"
+
+int mgaussdev (int argc, char **argv) {
+  
+  Buffer *buf;
+
+  if (argc != 6) goto usage;
+
+  if ((buf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  int Nx = atof (argv[2]);
+  int Ny = atof (argv[3]);
+
+  double mean = atof (argv[4]);
+  double sigma = atof (argv[5]);
+
+  /* I should encapsulate this in a create_default_buffer */
+  gfits_free_matrix (&buf[0].matrix);
+  gfits_free_header (&buf[0].header);
+
+  // 3D CUBE OPTION: if (!CreateBuffer3D (buf, Nx, Ny, Nz, -32, 1.0, 0.0)) return FALSE;
+  if (!CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
+
+  ohana_gaussdev_init ();
+  
+  float *v = (float *) buf[0].matrix.buffer;
+  for (int i = 0; i < Nx*Ny; i++, v++) {
+    *v = ohana_gaussdev_rnd (mean, sigma);
+  }
+  return (TRUE);
+
+ usage:
+  gprint (GP_ERR, "USAGE: mgaussdev (buff) Nx Ny mean sigma\n");
+  return (FALSE);
+}
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/test/medimage.sh
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/test/medimage.sh	(revision 40757)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/test/medimage.sh	(revision 40758)
@@ -1,11 +1,540 @@
 
-macro go
- mcreate a 30 30
- for i 0 40
-  set a$i = zero(a) + $i
- end
-
- for i 0 40
-  medimage add t1 a$i
- end
-end
+macro test.mean
+
+ $Nsample = 16
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 1.0
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ medimage calc t1 T -mean
+
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "expect {1/sqrt($Nsample)} : $C1"
+ plot -c red -x line x yf
+end
+
+macro test.median
+
+ $Nsample = 16
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 1.0
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+ end
+
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T
+
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "expect {1/sqrt($Nsample)} : $C1 (actually should be a bit higher)"
+ plot -c red -x line x yf
+end
+ 
+macro test.wtmean
+
+ $Nsample = 8
+ $sig1 = 1.0
+ $sig2 = 3.0
+
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig1
+  set v = $sig1^2 + zero(t)      
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig2
+  set v = $sig2^2 + zero(t)      
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T -wtmean
+
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
+ echo "expect {1/sqrt($S1)} : $C1"
+ plot -c red -x line x yf
+end
+
+macro test.irls
+ medimage delete -q t1
+ $Nsample = 16
+ $sig = 1.0
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig
+  set v = $sig^2 + zero(t)      
+
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean
+
+ imhist -q Tm x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from straight stdev: $C1"
+ # stats Tm
+ 
+ plot -c red -x line x yf
+
+ # get stats for irls
+ medimage calc t1 Ti -irls
+
+ imhist -q Ti x y -range -10 10 -delta 0.1
+ lim -n 2 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
+ # stats Ti
+ 
+ plot -c red -x line x yf
+end
+
+
+###################33
+
+
+macro test.mean.var
+
+ $Nsample = 64
+ $sig = 2.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 100 100 0.0 $sig
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ medimage calc t1 T -mean -variance Tv
+
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ plot -c red -x line x yf
+
+ imhist Tv xv yv -range -1 4 -delta 0.1
+ lim -n 2 xv yv; clear; box; plot xv yv -x hist
+
+ stat Tv
+ echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
+end
+
+macro test.median.var
+
+ $Nsample = 64
+ $sig = 2.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T -variance Tv
+
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ plot -c red -x line x yf
+
+ imhist Tv xv yv -range -1 4 -delta 0.1
+ lim -n 2 xv yv; clear; box; plot xv yv -x hist
+
+ stat Tv
+ echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
+end
+ 
+macro test.wtmean.var
+
+ $Nsample = 32
+ $sig1 = 1.0
+ $sig2 = 1.0
+
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig1
+  set v = $sig1^2 + zero(t)      
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig2
+  set v = $sig2^2 + zero(t)      
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T -wtmean -variance Tv
+
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ plot -c red -x line x yf
+
+ stat -q Tv
+ $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
+ echo $C1 vs {sqrt($MEDIAN)} : expect {1/sqrt($S1)}
+end
+
+macro test.irls.var
+
+ $Nsample = 16
+ $sig = 1.0
+
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig
+  set v = $sig^2 + zero(t)      
+
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean -variance Tv
+
+ imhist -q Tm x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from straight stdev: $C1"
+ # stats Tm
+ 
+ plot -c red -x line x yf
+
+ # get stats for irls
+ medimage calc t1 Ti -irls -variance Tv
+
+ imhist -q Ti x y -range -10 10 -delta 0.1
+ lim -n 2 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
+ # stats Ti
+ 
+ plot -c red -x line x yf
+
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+
+ stat -q Tv
+ echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
+end
+
+macro test.irls.boot.var
+
+ $Nsample = 64
+ $sig = 1.0
+
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 200 200 0.0 $sig
+  set v = $sig^2 + zero(t)      
+
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+
+  mgaussdev noise 200 200 0.0 0.5
+  set ts = ts + noise
+
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean -variance Tv
+
+ imhist -q Tm x y -range -10 10 -delta 0.02
+ lim -n 1 x y; clear; box; plot -x hist x y
+ peak -q x y
+ $C0 = $peakpos
+ $C1 = 1.5*$sig / sqrt($Nsample)
+ $C2 = $peakval
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from straight stdev: $C1"
+ # stats Tm
+ 
+ plot -c red -x line x yf
+
+ # get stats for irls
+ date
+ medimage calc t1 Ti -irls -variance Tv -bootstrap-iter 100
+ date
+
+ imhist -q Ti x y -range -10 10 -delta 0.02
+ lim -n 2 x y; clear; box; plot -x hist x y
+ peak -q x y
+ $C0 = $peakpos
+ $C1 = 1.5*$sig / sqrt($Nsample)
+ $C2 = $peakval
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
+ # stats Ti
+ 
+ plot -c red -x line x yf
+
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range 0 {5*$sig/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+
+ stat -q Tv
+ echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
+end
+
+##############################
+macro test.irls.boot.test
+
+ $Nsample = 100
+ $sig1 = 1.0
+
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 100 100 0.0 $sig1
+  set v = $sig1^2 + zero(t)      
+
+  medimage add t1 t -variance v
+ end
+
+ # get stats for irls
+ medimage calc t1 Ti -irls -variance Tv -bootstrap
+
+ imhist -q Ti x y -range {-10*$sig1/sqrt($Nsample)} {10*$sig1/sqrt($Nsample)} -delta 0.01
+ lim -n 2 x y; clear; box; plot -x hist x y
+ peak -q x y
+ $C0 = $peakpos
+ $C1 = 1.5*$sig1/sqrt($Nsample)
+ $C2 = $peakval
+ $C3 = 0
+ vgauss x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig1/sqrt($Nsample)})"
+ # stats Ti
+ 
+ plot -c red -x line x yf
+
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range 0 {5*$sig1/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+
+ stat -q irls_npt
+ $Npix = $MEAN
+
+ stat -q Tv
+ echo "sigma of irls average: $C1, sqrt(mean) of irls variance: {sqrt($MEAN)}, (ideal is {$sig1/sqrt($Npix)})"
+end
+
+macro test.irls.range.var
+ medimage delete -q t1
+ for i 0 8
+  mgaussdev t 50 50 0.0 1.0
+  set v = 1.0 + zero(t)      
+
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  echo $C1
+ end
+
+ for i 0 8
+  mgaussdev t 50 50 0.0 3.0
+  set v = 3.0 + zero(t)      
+
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  echo $C1
+ end
+
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean -variance Tv
+
+ imhist -q Tm x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo $C1
+ stats Tm
+ 
+ plot -c red -x line x yf
+
+ # get stats for irls
+ medimage calc t1 Ti -irls -variance Tv
+
+ imhist -q Ti x y -range -10 10 -delta 0.1
+ lim -n 2 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo $C1
+ stats Ti
+ 
+ plot -c red -x line x yf
+
+ stat -q Tv
+ echo $C1 vs {sqrt($MEDIAN)}
+
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+end