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

Timestamp:

Mar 31, 2019, 8:47:40 AM (7 years ago)

Author:

eugene

Message:

add vwtmean, virls, medimage wtmean & irls options, spline print

File:

• branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c (modified) (6 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) {
@@ -15 +20 @@
 int medimage_add (int argc, char **argv) {
 
+  int N;
   Buffer *image;
-
-  Buffer *var = NULL;
-  if ((N = get_argument (argc, argv, "-var"))) {
-    remove_argument (N, &argc, argv);
-    if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  Buffer *sig = NULL;
+
+  if ((N = get_argument (argc, argv, "-sigma"))) {
+    remove_argument (N, &argc, argv);
+    if ((sig = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
   }
 
   if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-var var]\n");
+    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-sigma sigma]\n");
     gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
     gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
@@ -33 +39 @@
   if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
 
-XXX: match dimensions wit var image
+  if (sig) {
+    if ((sig->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
+        (sig->matrix.Naxis[1] != image->matrix.Naxis[1])) {
+      gprint (GP_ERR, "sigma buffer does not match image buffer dimensions\n");
+      return FALSE;
+    }
+  }
 
   MedImageType *median = FindMedImage (argv[1]);
@@ -56 +68 @@
   // new image should match existing medimage dimensions
 
+  // AddMedImage (median, image, sig);
   int Ninput = median->Ninput;
   median->Ninput ++;
+  REALLOCATE (median->flx, float *, median->Ninput);
   REALLOCATE (median->sig, float *, median->Ninput);
-  REALLOCATE (median->var, float *, median->Ninput);
-
-  ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
-  memcpy (median->sig[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
-
-  median->var[Ninput] = NULL;
-  if (var) {
-    ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
-    memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
-  }
-
-  return TRUE;
-}
-
+
+  ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
+  memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+
+  median->sig[Ninput] = NULL;
+  if (sig) {
+    ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
+    memcpy (median->sig[Ninput], sig->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+  }
+
+  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 CALC_MEAN = FALSE;
+  
+  int mode = CALC_MEDIAN;
   if ((N = get_argument (argc, argv, "-mean"))) {
-    CALC_MEAN = TRUE;
+    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);
   }
 
   if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean]\n");
+    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;
@@ -102 +126 @@
   int Ny = median->Ny;
 
-  float *value = NULL;
-  ALLOCATE (value, float, Ninput);
+  ALLOCATE_PTR (val, float, Ninput);
+  ALLOCATE_PTR (wgt, float, Ninput);
 
   gfits_free_matrix (&output->matrix);
@@ -117 +141 @@
       int Npix = ix + Nx*iy;
       for (n = 0; n < Ninput; n++) {
-        float v = median->sig[n][Npix];
+        float v = median->flx[n][Npix];
         if (!isfinite(v)) continue;
-        value[N] = v;
+        val[N] = v;
+        wgt[N] = 1.0;
+        if (median->sig[n]) {
+          float s = median->sig[n][Npix];
+          if (!isfinite(s)) continue;
+          if (fabs(s) < 2*FLT_MIN) s = 2*FLT_MIN;
+          wgt[N] = 1.0 / SQ(s);
+        }
         N++;
       }
       if (N == 0) continue;
 
-    XXXX: add in calc options:
-      weighted mean
-      irls
-
-      if (CALC_MEAN) {
-        float sum = 0.0;
-        for (n = 0; n < N; n++) {
-          sum += value[n];
+      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;
         }
-        outvalue[Npix] = sum / (float) N;
-      } else {
-        fsort (value, N);
-        outvalue[Npix] = value[(int)(0.5*N)];
+        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;
-}
+  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));
+}
+
 
 /* 
-int medimage_save (int argc, char **argv) {
-
-  int N;
-
-  int APPEND = FALSE;
-  if ((N = get_argument (argc, argv, "-append"))) {
-    APPEND = TRUE;
-    remove_argument (N, &argc, argv);
-  }
-
-  if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
-    return FALSE;
-  }
-
-  if (!SaveMedImage(argv[2], argv[1], APPEND)) {
-    gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
-    return (FALSE);
-  }
-  return TRUE;
-}
-
-int medimage_load (int argc, char **argv) {
-
-  if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
-    return FALSE;
-  }
-
-  if (!LoadMedImage(argv[2], argv[1])) {
-    gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
-    return (FALSE);
-  }
-  return TRUE;
-}
+   int medimage_save (int argc, char **argv) {
+
+   int N;
+
+   int APPEND = FALSE;
+   if ((N = get_argument (argc, argv, "-append"))) {
+   APPEND = TRUE;
+   remove_argument (N, &argc, argv);
+   }
+
+   if (argc != 3) {
+   gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
+   return FALSE;
+   }
+
+   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
+   gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
+   return (FALSE);
+   }
+   return TRUE;
+   }
+
+   int medimage_load (int argc, char **argv) {
+
+   if (argc != 3) {
+   gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
+   return FALSE;
+   }
+
+   if (!LoadMedImage(argv[2], argv[1])) {
+   gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
+   return (FALSE);
+   }
+   return TRUE;
+   }
 */