Changeset 40758

Timestamp:

May 30, 2019, 6:05:51 AM (7 years ago)

Author:

eugene

Message:

update medimage to return the variance, finish irls in medimage, add mgaussdev

Location:

branches/eam_branches/ohana.20190329/src/opihi/cmd.data

Files:

• Makefile (modified) (2 diffs)
• init.c (modified) (2 diffs)
• medimage_calc.c (added)
• medimage_commands.c (modified) (2 diffs)
• mgaussdev.c (added)
• test/medimage.sh (modified) (1 diff)

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile

@@ -95 +95 @@
 $(SRC)/mcreate.$(ARCH).o        \
 $(SRC)/medacc.$(ARCH).o \
+$(SRC)/mgaussdev.$(ARCH).o      \
 $(SRC)/mget.$(ARCH).o           \
 $(SRC)/mget3d.$(ARCH).o         \
@@ -101 +102 @@
 $(SRC)/medimage.$(ARCH).o       \
 $(SRC)/medimage_commands.$(ARCH).o \
+$(SRC)/medimage_calc.$(ARCH).o \
 $(SRC)/mset.$(ARCH).o           \
 $(SRC)/needles.$(ARCH).o        \

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c

@@ -85 +85 @@
 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 +274 @@
   {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"},

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c

@@ -1 +1 @@
 # 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 +80 @@
   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));
-}
-
 
 /* 

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/test/medimage.sh

@@ -1 +1 @@
 
-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