IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 40662


Ignore:
Timestamp:
Mar 31, 2019, 8:47:40 AM (7 years ago)
Author:
eugene
Message:

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

Location:
branches/eam_branches/ohana.20190329/src/opihi
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

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

    r40637 r40662  
    182182$(SRC)/vsmooth.$(ARCH).o           \
    183183$(SRC)/vstats.$(ARCH).o            \
     184$(SRC)/virls.$(ARCH).o             \
     185$(SRC)/vwtmean.$(ARCH).o \
    184186$(SRC)/xsection.$(ARCH).o          \
    185187$(SRC)/vsh.$(ARCH).o               \

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

    r40637 r40662  
    163163int vzload           PROTO((int, char **));
    164164int vstats           PROTO((int, char **));
     165int vstats           PROTO((int, char **));
     166int virls            PROTO((int, char **));
     167int vwtmean          PROTO((int, char **));
    165168int vroll            PROTO((int, char **));
    166169int vshift           PROTO((int, char **));
    … …  
    360363  {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
    361364  {1, "vstats",       vstats,           "statistics on a vector"},
     365  {1, "vwtmean",      vwtmean,          "weighted mean of a vector"},
     366  {1, "virls",        virls,            "IRLS mean of a vector"},
    362367  {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
    363368  {1, "vsh",          vsh,              "Vector Spherical Harmonics"},

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

    r40661 r40662  
    11# include "data.h"
     2
     3float weight_cauchy_square_flt (float x2);
     4float irls_mean (float *val, float *wgt, int N);
     5
     6# define IRLS_TOLERANCE 1e-4
    27
    38int medimage_list (int argc, char **argv) {
    … …  
    1520int medimage_add (int argc, char **argv) {
    1621
     22  int N;
    1723  Buffer *image;
    18 
    19   Buffer *var = NULL;
    20   if ((N = get_argument (argc, argv, "-var"))) {
    21     remove_argument (N, &argc, argv);
    22     if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     24  Buffer *sig = NULL;
     25
     26  if ((N = get_argument (argc, argv, "-sigma"))) {
     27    remove_argument (N, &argc, argv);
     28    if ((sig = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    2329    remove_argument (N, &argc, argv);
    2430  }
    2531
    2632  if (argc != 3) {
    27     gprint (GP_ERR, "USAGE: medimage add (name) (image) [-var var]\n");
     33    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-sigma sigma]\n");
    2834    gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
    2935    gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
    … …  
    3339  if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    3440
    35 XXX: match dimensions wit var image
     41  if (sig) {
     42    if ((sig->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
     43        (sig->matrix.Naxis[1] != image->matrix.Naxis[1])) {
     44      gprint (GP_ERR, "sigma buffer does not match image buffer dimensions\n");
     45      return FALSE;
     46    }
     47  }
    3648
    3749  MedImageType *median = FindMedImage (argv[1]);
    … …  
    5668  // new image should match existing medimage dimensions
    5769
     70  // AddMedImage (median, image, sig);
    5871  int Ninput = median->Ninput;
    5972  median->Ninput ++;
     73  REALLOCATE (median->flx, float *, median->Ninput);
    6074  REALLOCATE (median->sig, float *, median->Ninput);
    61   REALLOCATE (median->var, float *, median->Ninput);
    62 
    63   ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
    64   memcpy (median->sig[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    65 
    66   median->var[Ninput] = NULL;
    67   if (var) {
    68     ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
    69     memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    70   }
    71 
    72   return TRUE;
    73 }
    74 
     75
     76  ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
     77  memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
     78
     79  median->sig[Ninput] = NULL;
     80  if (sig) {
     81    ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
     82    memcpy (median->sig[Ninput], sig->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
     83  }
     84
     85  return TRUE;
     86}
     87
     88enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN};
    7589int medimage_calc (int argc, char **argv) {
    7690
    7791  int ix, iy, n, N;
    7892  Buffer *output;
    79 
    80   int CALC_MEAN = FALSE;
     93 
     94  int mode = CALC_MEDIAN;
    8195  if ((N = get_argument (argc, argv, "-mean"))) {
    82     CALC_MEAN = TRUE;
     96    mode = CALC_MEAN;
     97    remove_argument (N, &argc, argv);
     98  }
     99  if ((N = get_argument (argc, argv, "-irls"))) {
     100    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -irls, -wtmean\n"); return FALSE; }
     101    mode = CALC_IRLS;
     102    remove_argument (N, &argc, argv);
     103  }
     104  if ((N = get_argument (argc, argv, "-wtmean"))) {
     105    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -irls, -wtmean\n"); return FALSE; }
     106    mode = CALC_WTMEAN;
    83107    remove_argument (N, &argc, argv);
    84108  }
    85109
    86110  if (argc != 3) {
    87     gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean]\n");
     111    gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean,-irls,-wtmean]\n");
    88112    gprint (GP_ERR, "       calculate the median image for the median image set\n");
    89113    return FALSE;
    … …  
    102126  int Ny = median->Ny;
    103127
    104   float *value = NULL;
    105   ALLOCATE (value, float, Ninput);
     128  ALLOCATE_PTR (val, float, Ninput);
     129  ALLOCATE_PTR (wgt, float, Ninput);
    106130
    107131  gfits_free_matrix (&output->matrix);
    … …  
    117141      int Npix = ix + Nx*iy;
    118142      for (n = 0; n < Ninput; n++) {
    119         float v = median->sig[n][Npix];
     143        float v = median->flx[n][Npix];
    120144        if (!isfinite(v)) continue;
    121         value[N] = v;
     145        val[N] = v;
     146        wgt[N] = 1.0;
     147        if (median->sig[n]) {
     148          float s = median->sig[n][Npix];
     149          if (!isfinite(s)) continue;
     150          if (fabs(s) < 2*FLT_MIN) s = 2*FLT_MIN;
     151          wgt[N] = 1.0 / SQ(s);
     152        }
    122153        N++;
    123154      }
    124155      if (N == 0) continue;
    125156
    126     XXXX: add in calc options:
    127       weighted mean
    128       irls
    129 
    130       if (CALC_MEAN) {
    131         float sum = 0.0;
    132         for (n = 0; n < N; n++) {
    133           sum += value[n];
     157      switch (mode) {
     158        case CALC_MEDIAN:
     159          fsort (val, N);
     160          outvalue[Npix] = val[(int)(0.5*N)];
     161          break;
     162        case CALC_MEAN: {
     163          float sum = 0.0;
     164          for (n = 0; n < N; n++) {
     165            sum += val[n];
     166          }
     167          outvalue[Npix] = sum / (float) N;
     168          break;
    134169        }
    135         outvalue[Npix] = sum / (float) N;
    136       } else {
    137         fsort (value, N);
    138         outvalue[Npix] = value[(int)(0.5*N)];
     170        case CALC_WTMEAN: {
     171          float S1 = 0.0, S2 = 0.0;
     172          for (n = 0; n < N; n++) {
     173            S1 += wgt[n] * val[n];
     174            S2 += wgt[n];
     175          }
     176          outvalue[Npix] = S1 / S2;
     177          break;
     178        }
     179        case CALC_IRLS:
     180          outvalue[Npix] = irls_mean (val, wgt, N);
    139181      }
    140182    }
    141183  }
    142 
    143   return TRUE;
    144 }
     184  return TRUE;
     185}
     186
     187float irls_mean (float *val, float *wgt, int N) {
     188
     189  // calculate weighted mean
     190  float S1 = 0.0, S2 = 0.0;
     191  for (int n = 0; n < N; n++) {
     192    S1 += wgt[n] * val[n];
     193    S2 += wgt[n];
     194  }
     195  float Value = S1 / S2;
     196 
     197  int converged = FALSE;
     198  for (int i = 0; (i < 10) && !converged; i++) {
     199    float ValueLast = Value;
     200
     201    float S1 = 0.0, S2 = 0.0;
     202   
     203    // calculate weight modification based on distances (squared).
     204    // use modifier to calculate new weighted mean
     205    for (int n = 0; n < N; n++) {
     206      float dV = (val[n] - Value);
     207      float d2 = SQ(dV) * wgt[n];
     208     
     209      float Mod = weight_cauchy_square_flt (d2);
     210      S1 += Mod * wgt[n] * val[n];
     211      S2 += Mod * wgt[n];
     212    }
     213    Value = S1 / S2;
     214
     215    float delta = fabs(Value - ValueLast);
     216    if (delta < Value * IRLS_TOLERANCE) converged = TRUE;
     217  }
     218  return Value;
     219}
     220
     221// exp(-(x^2/s^2)/2) = (1/2)
     222//     -(x^2/s^2)/2  = ln(1/2)
     223//      (x^2/s^2)/2  = ln(2)
     224//      (x^2/s^2)    = 2ln(2)
     225//      (x  /s)      = sqrt(2ln(2)) : half-width at half-max
     226//       FWHM        = 2sqrt(2ln(2))
     227
     228// R2 = (X / 2.385)^2 = (X^2 / 2.385^2)
     229
     230# define CAUCY_FACTOR 1.0
     231
     232float weight_cauchy_square_flt (float x2) {
     233  float r2 = x2 / CAUCY_FACTOR;
     234  return (1.0 / (1.0 + r2));
     235}
     236
    145237
    146238/*
    147 int medimage_save (int argc, char **argv) {
    148 
    149   int N;
    150 
    151   int APPEND = FALSE;
    152   if ((N = get_argument (argc, argv, "-append"))) {
    153     APPEND = TRUE;
    154     remove_argument (N, &argc, argv);
    155   }
    156 
    157   if (argc != 3) {
    158     gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
    159     return FALSE;
    160   }
    161 
    162   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
    163     gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
    164     return (FALSE);
    165   }
    166   return TRUE;
    167 }
    168 
    169 int medimage_load (int argc, char **argv) {
    170 
    171   if (argc != 3) {
    172     gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
    173     return FALSE;
    174   }
    175 
    176   if (!LoadMedImage(argv[2], argv[1])) {
    177     gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
    178     return (FALSE);
    179   }
    180   return TRUE;
    181 }
     239   int medimage_save (int argc, char **argv) {
     240
     241   int N;
     242
     243   int APPEND = FALSE;
     244   if ((N = get_argument (argc, argv, "-append"))) {
     245   APPEND = TRUE;
     246   remove_argument (N, &argc, argv);
     247   }
     248
     249   if (argc != 3) {
     250   gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
     251   return FALSE;
     252   }
     253
     254   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
     255   gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
     256   return (FALSE);
     257   }
     258   return TRUE;
     259   }
     260
     261   int medimage_load (int argc, char **argv) {
     262
     263   if (argc != 3) {
     264   gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
     265   return FALSE;
     266   }
     267
     268   if (!LoadMedImage(argv[2], argv[1])) {
     269   gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
     270   return (FALSE);
     271   }
     272   return TRUE;
     273   }
    182274*/
    183275

  • branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline.c

    r39457 r40662  
    77int spline_load (int argc, char **argv);
    88int spline_save (int argc, char **argv);
     9int spline_print (int argc, char **argv);
    910int spline_delete (int argc, char **argv);
    1011int spline_rename (int argc, char **argv);
    … …  
    1819  {1, "load",       spline_load,       "write a spline to a FITS file"},
    1920  {1, "save",       spline_save,       "read a spline from a FITS file"},
     21  {1, "print",      spline_print,      "print a spline to stdout"},
    2022  {1, "delete",     spline_delete,     "delete a spline"},
    2123  {1, "rename",     spline_rename,     "rename a spline"},
    … …  
    3537  gprint (GP_ERR, "    spline load   (spline) (filename)            : load a spline from a FITS file\n");
    3638  gprint (GP_ERR, "    spline save   (spline) (filename) [-append]  : save a spline in FITS format\n");
     39  gprint (GP_ERR, "    spline print  (spline)                       : print a spline to stdout\n");
    3740
    3841  return FALSE;

  • branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline_commands.c

    r39457 r40662  
    9393}
    9494
     95int spline_print (int argc, char **argv) {
     96
     97  int i;
     98
     99  if (argc != 2) {
     100    gprint (GP_ERR, "USAGE: spline print (name)\n");
     101    return FALSE;
     102  }
     103
     104  Spline *myspline = FindSpline (argv[1]);
     105  if (!myspline) {
     106    gprint (GP_ERR, "spline %s not found\n", argv[1]);
     107    return (FALSE);
     108  }
     109
     110  for (i = 0; i < myspline->Nknots; i++) {
     111    gprint (GP_LOG, "%lf : %lf : %lf\n", myspline->xk[i], myspline->yk[i], myspline->y2[i]);
     112  }   
     113
     114  return TRUE;
     115}
     116
    95117int spline_load (int argc, char **argv) {
    96118

  • branches/eam_branches/ohana.20190329/src/opihi/include/data.h

    r40661 r40662  
    4848  int Nx;
    4949  int Ny;
     50  float **flx;
    5051  float **sig;
    51   float **var;
    5252} MedImageType;
    5353

  • branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c

    r40661 r40662  
    3232  free (medimage[0].name);
    3333  for (i = 0; i < medimage[0].Ninput; i++) {
    34     free (medimage[0].sig[i]);
     34    free (medimage[0].flx[i]);
     35    FREE (medimage[0].sig[i]);
    3536  }
     37  free (medimage[0].flx);
    3638  free (medimage[0].sig);
    3739  free (medimage);
    … …  
    7577  medimage->Nx = Nx;
    7678  medimage->Ny = Ny;
     79  ALLOCATE (medimage->flx, float *, 1);
    7780  ALLOCATE (medimage->sig, float *, 1);
    7881

  • branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c

    r40661 r40662  
    190190    if (psf) Npof -= Npsf;
    191191
     192    // integral and fractional pixel offsets of PSF (none, by default)
     193    int dxi = 0;
     194    float dxf = 0.0, dxr = 1.0;
     195   
    192196    if (psf_trace) {
    193       float dx = -spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
    194       Npof += dx;
    195     }
     197      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
     198      dxi = floor(dx);
     199      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
     200      dxr = 1 - dxf;
     201      Npof += dxi;
     202      if (iy % 200 == 0) { fprintf (stderr, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
     203    }
     204
     205    // if fractional offset is small, do not interpolate
     206    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
    196207
    197208    // flux = obj * PSF + sky
    … …  
    205216        // only add in the flux if we are in range of the PSF
    206217        if ((n >= 0) && (n < NpsfFull)) {
    207           value += objVy * psfV[n];
     218          if (doInterp) {
     219            if ((n > 0) && (n < NpsfFull)) {
     220              value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     221            }
     222            if (n == 0) {
     223              value += objVy*psfV[n]*dxr;
     224            }
     225          } else {
     226            value += objVy*psfV[n];
     227          }
    208228        }
    209229      } else {
Note: See TracChangeset for help on using the changeset viewer.