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

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

File:

• branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c (modified) (2 diffs)

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));
-}
-
 
 /*