IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 16, 2020, 2:04:27 PM (6 years ago)
Author:
tdeboer
Message:

quick undo

Location:
trunk/Ohana/src/opihi/cmd.data
Files:
23 edited
16 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.data/Makefile

    r41340 r41341  
    1919$(SRC)/init.$(ARCH).o          \
    2020$(SRC)/accum.$(ARCH).o          \
     21$(SRC)/antialias.$(ARCH).o              \
    2122$(SRC)/applyfit1d.$(ARCH).o     \
    2223$(SRC)/applyfit2d.$(ARCH).o     \
    … …  
    7879$(SRC)/join.$(ARCH).o           \
    7980$(SRC)/jpeg.$(ARCH).o           \
     81$(SRC)/kapamemory.$(ARCH).o             \
    8082$(SRC)/kern.$(ARCH).o           \
    8183$(SRC)/keyword.$(ARCH).o        \
    … …  
    8991$(SRC)/lookup.$(ARCH).o \
    9092$(SRC)/matrix.$(ARCH).o \
     93$(SRC)/match1d.$(ARCH).o        \
    9194$(SRC)/match2d.$(ARCH).o        \
    9295$(SRC)/mkrgb.$(ARCH).o  \
    9396$(SRC)/mcreate.$(ARCH).o        \
    9497$(SRC)/medacc.$(ARCH).o \
     98$(SRC)/mgaussdev.$(ARCH).o      \
    9599$(SRC)/mget.$(ARCH).o           \
    96100$(SRC)/mget3d.$(ARCH).o         \
    … …  
    99103$(SRC)/medimage.$(ARCH).o       \
    100104$(SRC)/medimage_commands.$(ARCH).o \
     105$(SRC)/medimage_calc.$(ARCH).o \
    101106$(SRC)/mset.$(ARCH).o           \
    102107$(SRC)/needles.$(ARCH).o        \
    … …  
    112117$(SRC)/point.$(ARCH).o          \
    113118$(SRC)/ps.$(ARCH).o             \
     119$(SRC)/pdf.$(ARCH).o            \
    114120$(SRC)/print_vectors.$(ARCH).o  \
    115121$(SRC)/mprint.$(ARCH).o         \
    … …  
    124130$(SRC)/queuesubstr.$(ARCH).o    \
    125131$(SRC)/queueinit.$(ARCH).o      \
     132$(SRC)/queue2book.$(ARCH).o     \
    126133$(SRC)/radial.$(ARCH).o \
    127134$(SRC)/rd.$(ARCH).o             \
    … …  
    169176$(SRC)/vgauss.$(ARCH).o            \
    170177$(SRC)/vlorentz.$(ARCH).o          \
     178$(SRC)/vsigmoid.$(ARCH).o          \
    171179$(SRC)/vellipse.$(ARCH).o          \
    172180$(SRC)/vmaxwell.$(ARCH).o          \
    … …  
    176184$(SRC)/vzload.$(ARCH).o            \
    177185$(SRC)/vpeaks.$(ARCH).o            \
     186$(SRC)/vtransitions.$(ARCH).o     \
    178187$(SRC)/vpop.$(ARCH).o              \
    179188$(SRC)/vroll.$(ARCH).o             \
    … …  
    182191$(SRC)/vsmooth.$(ARCH).o           \
    183192$(SRC)/vstats.$(ARCH).o            \
     193$(SRC)/virls.$(ARCH).o             \
     194$(SRC)/vwtmean.$(ARCH).o \
    184195$(SRC)/xsection.$(ARCH).o          \
    185196$(SRC)/vsh.$(ARCH).o               \

  • trunk/Ohana/src/opihi/cmd.data/create.c

    r41340 r41341  
    77  Vector *vec;
    88 
     9  // create a vector of empty strings
     10  if ((N = get_argument (argc, argv, "-str"))) {
     11    remove_argument (N, &argc, argv);
     12
     13    char *stringValue = NULL;
     14    if ((N = get_argument (argc, argv, "-value"))) {
     15      remove_argument (N, &argc, argv);
     16      stringValue = strcreate (argv[N]);
     17      remove_argument (N, &argc, argv);
     18    }
     19
     20    if (argc != 3) {
     21      gprint (GP_ERR, "USAGE: create vector Nelements -value word\n");
     22      return (FALSE);
     23    }
     24
     25    if ((vec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     26    int Nelements = atoi (argv[2]);
     27
     28    // a newly reset vector has NULL-valued pointers
     29    ResetVector (vec, OPIHI_STR, Nelements);
     30    for (int i = 0; i < Nelements; i++) {
     31      vec[0].elements.Str[i] = stringValue ? strcreate (stringValue) : strcreate ("");
     32    }
     33    return TRUE;
     34  }
     35
    936  INT = FALSE;
    1037  if ((N = get_argument (argc, argv, "-int"))) {
    … …  
    1643    gprint (GP_ERR, "USAGE: create vector start end [delta] [-int]\n");
    1744    gprint (GP_ERR, " -int : resulting vector is integer type (delta must be integer)\n");
     45    gprint (GP_ERR, " -str : resulting vector is string type (only give number of elements)\n");
    1846    return (FALSE);
    1947  }

  • trunk/Ohana/src/opihi/cmd.data/cut.c

    r41340 r41341  
    11# include "data.h"
    22
    3 enum {SUM, MEAN, MEDIAN};
     3enum {SUM, MEAN, MEDIAN, INNER, NGOOD};
     4
     5double cutstat (float *value, int Nvalue, int mode) {
     6
     7  double output = 0.0;
     8
     9  if (Nvalue == 0) return NAN;
     10
     11  if (mode == MEDIAN) {
     12    fsort (value, Nvalue);
     13    if (Nvalue % 2) {
     14      int Ncenter = Nvalue / 2;
     15      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue), "cutstat");
     16      output = value[Ncenter];
     17    } else {
     18      int Ncenter = Nvalue / 2 - 1;
     19      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue - 1), "cutstat");
     20      output = 0.5*(value[Ncenter] + value[Ncenter + 1]);
     21    }
     22    return output;
     23  }
     24  if (mode == INNER) {
     25    fsort (value, Nvalue);
     26
     27    int Ns = 0, Ne = 0;
     28    if (Nvalue % 2) {
     29      // for an odd number of points take the same number below
     30      // and above the center value
     31      int Ncenter = Nvalue / 2;
     32      int Nquarter = 0.25*Nvalue;
     33      Ns = Ncenter - Nquarter;
     34      Ne = Ncenter + Nquarter;
     35    } else {
     36      // for an even number, the middle lies between two points
     37      // take the same number below as above
     38      int Ncenter = (int)(Nvalue / 2) - 1;
     39      int Nquarter = 0.25*Nvalue;
     40      Ns = Ncenter     - Nquarter;
     41      Ne = Ncenter + 1 + Nquarter;
     42    }
     43    int Nv = 0;
     44    for (int i = Ns; i <= Ne; i++) {
     45      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
     46      output += value[i];
     47      Nv ++;
     48    }
     49    output /= Nv;
     50    return output;
     51  }
     52  if (mode == NGOOD) {
     53    return Nvalue;
     54  }
     55  if ((mode == MEAN) || (mode == SUM)) {
     56    for (int i = 0; i < Nvalue; i++) {
     57      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
     58      output += value[i];
     59    }
     60    if (mode == MEAN) { output /= Nvalue; }
     61    return output;
     62  }
     63  return NAN;
     64}
    465
    566int cut (int argc, char **argv) {
    667 
    768  int i, j, N, Nx, Ny, Mode;
    8   float *Vin, *Vbuf, value;
     69  float *Vin, *Vbuf;
    970  int sx, sy, nx, ny;
    1071  Vector *xvec, *yvec;
    … …  
    1273
    1374  Mode = SUM;
     75  if ((N = get_argument (argc, argv, "-sum"))) {
     76    remove_argument (N, &argc, argv);
     77    Mode = SUM;
     78  }
    1479  if ((N = get_argument (argc, argv, "-median"))) {
    1580    remove_argument (N, &argc, argv);
    … …  
    1984    remove_argument (N, &argc, argv);
    2085    Mode = MEAN;
     86  }
     87  if ((N = get_argument (argc, argv, "-inner"))) {
     88    remove_argument (N, &argc, argv);
     89    Mode = INNER;
     90  }
     91  if ((N = get_argument (argc, argv, "-ngood"))) {
     92    remove_argument (N, &argc, argv);
     93    Mode = NGOOD;
    2194  }
    2295
    … …  
    36109  Ny = buf[0].matrix.Naxis[1];
    37110
    38   if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) {
    39     gprint (GP_ERR, "region out of range\n");
     111  // ny & nx do not need to be constrained by the buffer, but they do need to be sensible
     112  if ((nx < 0) || (ny < 0) || (nx > 1e8) || (ny > 1e8)) {
     113    gprint (GP_ERR, "warning : extraction size is crazy: %d,%d\n", nx, ny);
    40114    return (FALSE);
    41115  }
    … …  
    57131
    58132    for (i = 0; i < nx; i++) {
     133      // for out-of-range areas, set the output pixels to NAN
     134      if (i + sx < 0) {
     135        yvec[0].elements.Flt[i] = NAN;
     136        continue;
     137      }
     138      if (i + sx >= Nx) {
     139        yvec[0].elements.Flt[i] = NAN;
     140        continue;
     141      }
     142
    59143      /* accumulate values */
     144      // skip out-of-range areas in y
     145      if (sy < 0) {
     146        ny += sy;
     147        sy = 0;
     148      }
     149      if (sy > Ny) sy = Ny;
     150      if (sy + ny >= Ny) ny = Ny - sy;
     151
    60152      Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i;
     153      int Npix = 0;
    61154      for (j = 0; j < ny; j++, Vin += Nx) {
    62         Vbuf[j] = *Vin;
    63       }
    64       /* apply stat of choice */
    65       if (Mode == MEDIAN) {
    66         fsort (Vbuf, ny);
    67         value = Vbuf[(int)(0.5*ny)];
    68       } else {
    69         value = 0;
    70         for (j = 0; j < ny; j++) {
    71           value += Vbuf[j];
    72         }
    73         if (Mode == MEAN) { value /= ny; }
    74       }
    75       yvec[0].elements.Flt[i] = value;
     155        if (!isfinite(*Vin)) continue;
     156        Vbuf[Npix] = *Vin;
     157        Npix ++;
     158      }
     159      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
    76160    }
    77161    free (Vbuf);
    … …  
    89173
    90174    for (i = 0; i < ny; i++) {
     175      // for out-of-range areas, set the output pixels to NAN
     176      if (i + sy < 0) {
     177        yvec[0].elements.Flt[i] = NAN;
     178        continue;
     179      }
     180      if (i + sy >= Ny) {
     181        yvec[0].elements.Flt[i] = NAN;
     182        continue;
     183      }
     184
    91185      /* accumulate values */
     186      // skip out-of-range areas in x
     187      if (sx < 0) {
     188        nx += sx;
     189        sx = 0;
     190      }
     191      if (sx > Nx) sx = Nx;
     192      if (sx + nx >= Nx) nx = Nx - sx;
     193
    92194      Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx;
     195      int Npix = 0;
    93196      for (j = 0; j < nx; j++, Vin ++) {
    94         Vbuf[j] = *Vin;
    95       }
    96       /* apply stat of choice */
    97       if (Mode == MEDIAN) {
    98         fsort (Vbuf, nx);
    99         value = Vbuf[(int)(0.5*nx)];
    100       } else {
    101         value = 0;
    102         for (j = 0; j < nx; j++) {
    103           value += Vbuf[j];
    104         }
    105         if (Mode == MEAN) { value /= nx; }
    106       }
    107       yvec[0].elements.Flt[i] = value;
     197        if (!isfinite(*Vin)) continue;
     198        Vbuf[Npix] = *Vin;
     199        Npix ++;
     200      }
     201      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
    108202    }
    109203    free (Vbuf);

  • trunk/Ohana/src/opihi/cmd.data/dbselect.c

    r41340 r41341  
    8080      case FIELD_TYPE_DOUBLE:
    8181      case FIELD_TYPE_TIME:
     82      case FIELD_TYPE_TIMESTAMP:
    8283      case FIELD_TYPE_DATE:
    8384      case FIELD_TYPE_DATETIME:
    … …  
    115116            break;
    116117          case FIELD_TYPE_TIME:
     118          case FIELD_TYPE_TIMESTAMP:
    117119          case FIELD_TYPE_DATE:
    118120          case FIELD_TYPE_DATETIME:

  • trunk/Ohana/src/opihi/cmd.data/extract.c

    r41340 r41341  
    11# include "data.h"
     2
     3/* <from> : source image, must exist
     4   <to>   : target image -- if it does not exist, it is created of size (Nx,Ny)
     5   sx, sy : source starting coordinate -- need not be on a valid image pixel
     6   nx, ny : number of source pixels -- currently must not go out of bounds -> allow out-of-bounds
     7   Sx, Sy : target starting coordinate -- need not be on a valid image pixel
     8   Nx, Ny : redundant information UNLESS <to> does exist (in which case it is required information) -> make Nx,Ny optional if <to> exists?
     9 */
    210
    311int extract (int argc, char **argv) {
    412 
    5   int i, j;
    6   float *Vin, *Vout;
    7   int sx, sy, nx, ny, NX, NY;
    8   int Sx, Sy, Nx, Ny;
     13  int N;
    914  Buffer *in, *out;
     15
     16  float initValue = 0.0;
     17  int initOutput = FALSE;
     18  if ((N = get_argument (argc, argv, "-init"))) {
     19    remove_argument (N, &argc, argv);
     20    initValue = atof (argv[N]);
     21    remove_argument (N, &argc, argv);
     22    initOutput = TRUE;
     23  }
    1024
    1125  if (argc != 11) {
    … …  
    1529
    1630  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    17   NX = in[0].matrix.Naxis[0];
    18   NY = in[0].matrix.Naxis[1];
     31  int NX = in[0].matrix.Naxis[0];
     32  int NY = in[0].matrix.Naxis[1];
    1933
    20   sx = atof (argv[3]);
    21   sy = atof (argv[4]);
    22   nx = atof (argv[5]);
    23   ny = atof (argv[6]);
     34  int sx = atof (argv[3]);
     35  int sy = atof (argv[4]);
     36  int nx = atof (argv[5]);
     37  int ny = atof (argv[6]);
    2438
    25   Sx = atof (argv[7]);
    26   Sy = atof (argv[8]);
    27   Nx = atof (argv[9]);
    28   Ny = atof (argv[10]);
     39  int Sx = atof (argv[7]);
     40  int Sy = atof (argv[8]);
     41  int Nx = atof (argv[9]);
     42  int Ny = atof (argv[10]);
    2943
    3044  if ((Sy + ny > Ny) || (Sx + nx > Nx)) {
    31     gprint (GP_ERR, "mismatch between source and dest regions\n");
     45    gprint (GP_ERR, "source pixels extend beyond target pixels\n");
    3246    gprint (GP_ERR, "%d + %d > %d or %d + %d > %d\n", Sy, ny, Ny, Sx, nx, Nx);
    33     return (FALSE);
     47    // return (FALSE);
    3448  }
     49
     50  // XXX : allow source region to fall outside source image
    3551
    3652  /* region is not on first image */
    … …  
    3955      (sy > in[0].matrix.Naxis[1])) {
    4056    gprint (GP_ERR, "region outside of source image\n");
    41     return (FALSE);
     57    // return (FALSE);
    4258  }
    4359
    44   if ((Sx + nx > Nx) || (Sy + ny > Ny)) {
    45     gprint (GP_ERR, "source region larger than dest region\n");
    46     return (FALSE);
    47   }
    4860  if ((Sx < 0) || (Sy < 0)) {
    4961    gprint (GP_ERR, "dest region out of range\n");
    50     return (FALSE);
     62    // return (FALSE);
    5163  }
    5264
    … …  
    7688  }
    7789
    78   for (j = 0; j < ny; j++) {
    79     if (j + sy < 0) continue;
    80     if (j + sy >= NY) continue;
    81     Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx; 
    82     Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;   
    83     for (i = 0; i < nx; i++, Vin++, Vout++) {
     90  // (NX,NY) : source image dimensions
     91  // (Nx,Ny) : target image dimensions
     92
     93  // allow (sx, sy), (Sx, Sy) to be off image
     94  // allow (sx+nx, sy+ny) to be off image
     95
     96  // if (sx < 0) {
     97  //   nx += sx;
     98  //   sx = 0;
     99  // }
     100  // if (sx > NX) sx = NX;
     101  //
     102  // if (Sx < 0) {
     103  //   Nx += sx;
     104  //   sx = 0;
     105  // }
     106  // if (sx > NX) sx = NX;
     107 
     108  float *Vin = (float *)(in[0].matrix.buffer);
     109  float *Vout = (float *)(out[0].matrix.buffer);
     110
     111  if (initOutput) {
     112    for (int j = 0; j < Ny; j++) {
     113      for (int i = 0; i < Nx; i++) {
     114        Vout[i + Nx*j] = initValue;
     115      }
     116    }
     117  }
     118
     119  int Nps = NX*NY;
     120  int Npt = Nx*Ny;
     121
     122  for (int j = 0; j < ny; j++) {
     123    if (j + sy < 0) continue; // not yet on source image
     124    if (j + Sy < 0) continue; // not yet on target image
     125    if (j + sy >= NY) continue; // past edge of source image
     126    if (j + Sy >= Ny) continue; // past edge of target image
     127
     128    // Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx; 
     129    // Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;   
     130
     131    for (int i = 0; i < nx; i++) {
    84132      if (i + sx < 0) continue;
    85133      if (i + sx >= NX) continue;
    86       *Vout = *Vin;
     134      if (i + Sx < 0) continue;
     135      if (i + Sx >= Nx) continue;
     136
     137      int ps = i + sx + (j + sy)*NX;
     138      int pt = i + Sx + (j + Sy)*Nx;
     139
     140      myAssert (pt >= 0, "oops 1");
     141      myAssert (pt < Npt, "oops 2");
     142
     143      myAssert (ps >= 0, "oops 3");
     144      myAssert (ps < Nps, "oops 4");
     145
     146      Vout[pt] = Vin[ps];
    87147    }
    88148  }
    … …  
    91151
    92152}
    93 

  • trunk/Ohana/src/opihi/cmd.data/init.c

    r41340 r41341  
    33
    44int accum            PROTO((int, char **));
     5int antialias        PROTO((int, char **));
    56int applyfit         PROTO((int, char **));
    67int applyfit1d       PROTO((int, char **));
    … …  
    6768int join             PROTO((int, char **));
    6869int jpeg             PROTO((int, char **));
     70int kapamemory       PROTO((int, char **));
    6971int kern             PROTO((int, char **));
    7072int keyword          PROTO((int, char **));
    … …  
    7981int lookup           PROTO((int, char **));
    8082int matrix           PROTO((int, char **));
     83int match1d          PROTO((int, char **));
    8184int match2d          PROTO((int, char **));
    8285int mkrgb            PROTO((int, char **));
    8386int mcreate          PROTO((int, char **));
    8487int medacc           PROTO((int, char **));
     88int mgaussdev        PROTO((int, char **));
    8589int mget             PROTO((int, char **));
    8690int mget3d           PROTO((int, char **));
    … …  
    99103int parity           PROTO((int, char **));
    100104int point            PROTO((int, char **));
     105int pdf              PROTO((int, char **));
    101106int ps               PROTO((int, char **));
    102107int vprint           PROTO((int, char **));
    … …  
    111116int queuesubstr      PROTO((int, char **));
    112117int queuesize        PROTO((int, char **));
     118int queue2book       PROTO((int, char **));
    113119int rd               PROTO((int, char **));
    114120int rdseg            PROTO((int, char **));
    … …  
    156162int vgauss           PROTO((int, char **));
    157163int vlorentz         PROTO((int, char **));
     164int vsigmoid         PROTO((int, char **));
    158165int vellipse         PROTO((int, char **));
    159166int vmaxwell         PROTO((int, char **));
    … …  
    163170int vzload           PROTO((int, char **));
    164171int vstats           PROTO((int, char **));
     172int vstats           PROTO((int, char **));
     173int virls            PROTO((int, char **));
     174int vwtmean          PROTO((int, char **));
    165175int vroll            PROTO((int, char **));
    166176int vshift           PROTO((int, char **));
    167177int vpeaks           PROTO((int, char **));
     178int vtransitions     PROTO((int, char **));
    168179int vpop             PROTO((int, char **));
    169180int vsmooth          PROTO((int, char **));
    … …  
    187198static Command cmds[] = { 
    188199  {1, "accum",        accum,            "accumulate vector values in another vector"},
     200  {1, "antialias",    antialias,        "set anti-alias sigma value for display"},
    189201  {1, "applyfit",     applyfit1d,       "apply 1-d fit to new vector"},
    190202  {1, "applyfit1d",   applyfit1d,       "apply 1-d fit to new vector"},
    … …  
    255267  {1, "join",         join,             "find the join of two ID vectors"},
    256268  {1, "jpeg",         jpeg,             "convert display image to JPEG"},
     269  {1, "kapamemory",   kapamemory,       "manage kapa memory dump options"},
    257270  {1, "kern",         kern,             "convolve with 3x3 kernel"},
    258271  {1, "keyword",      keyword,          "extract a FITS keyword from image header"},
    … …  
    265278  {1, "imcreate",     mcreate,          "create an image"},
    266279  {1, "medacc",       medacc,           "accumulate vector values in another vector"},
     280  {1, "mgaussdev",    mgaussdev,        "generate a gaussian deviate image"},
    267281  {1, "mget",         mget,             "extract a vector from an image"},
    268282  {1, "mget3d",       mget3d,           "extract a vector from a 3D image"},
    … …  
    273287  {1, "medimage",     medimage_command, "median image manipulation"},
    274288  {1, "matrix",       matrix,           "matrix math operations"},
     289  {1, "match1d",      match1d,          "match 2 vectors and return matched indexes"},
    275290  {1, "match2d",      match2d,          "match 2 pairs of X,Y vectors and return matched indexes"},
    276291  {1, "mkrgb",        mkrgb,            "convert 3 images to rgb jpeg (use Kapa for better control)"},
    … …  
    289304  {1, "ppm",          jpeg,             "convert display graphic to PPM"},
    290305  {1, "ps",           ps,               "convert display to PostScript"},
     306  {1, "pdf",          pdf,              "convert display to PDF"},
    291307  {1, "print_vectors", vprint,          "print a set of vectors"},
    292308  {1, "vprint",       vprint,           "print a set of vectors"},
    … …  
    301317  {1, "queuesize",    queuesize,        "show queue size"},
    302318  {1, "queuesubstr",  queuesubstr,      "bulk replace strings in queue"},
     319  {1, "queue2book",   queue2book,       "convert queue with ipptool output to book"},
     320  {1, "ipptool2book", queue2book,       "convert queue with ipptool output to book"},
    303321  {1, "rd",           rd,               "load fits image"},
    304322  {1, "rdseg",        rdseg,            "read a segment of an image from a file"},
    … …  
    347365  {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
    348366  {1, "vlorentz",     vlorentz,         "fit a Lorentzian to a vector"},
     367  {1, "vsigmoid",     vsigmoid,         "fit a Sigmoid to a vector"},
    349368  {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
    350369  {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
    … …  
    354373  {1, "vload",        vload,            "load vectors as overlay on image display"},
    355374  {1, "vmaxwell",     vmaxwell,         "fit a Maxwellian to a vector"},
    356   {1, "vpeaks",       vpeaks,           "fine coord and flux of peaks in vector"},
     375  {1, "vpeaks",       vpeaks,           "find coord and flux of peaks in vector"},
     376  {1, "vtransitions", vtransitions,     "find points in vector that cross the transition value"},
    357377  {1, "vpop",         vpop,             "remove first element of a vector"},
    358378  {1, "vroll",        vroll,            "roll vector elements by 1 entry"},
    … …  
    360380  {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
    361381  {1, "vstats",       vstats,           "statistics on a vector"},
     382  {1, "vwtmean",      vwtmean,          "weighted mean of a vector"},
     383  {1, "virls",        virls,            "IRLS mean of a vector"},
    362384  {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
    363385  {1, "vsh",          vsh,              "Vector Spherical Harmonics"},

  • trunk/Ohana/src/opihi/cmd.data/interpolate_presort.c

    r41340 r41341  
    77  double x0, dx, dy, y0;
    88  Vector *xout, *yout, *xin, *yin;
     9
     10  int N;
     11  int FillEnds = FALSE;
     12  if ((N = get_argument (argc, argv, "-fill-ends"))) {
     13    FillEnds = TRUE;
     14    remove_argument (N, &argc, argv);
     15  }
    916
    1017  /** check basic syntax **/
    … …  
    3643  y0 = yin[0].elements.Flt[0];
    3744 
    38   /* in and out vectors are sorted */
    39   j = 0;
     45  // fill in the start with NANs
     46  for (i = 0; (i < xout[0].Nelements) && (xout[0].elements.Flt[i] < xin[0].elements.Flt[0]); i++) {
     47    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[0] : NAN;
     48  }
    4049
    4150  // every output pixel should get a value unless
    4251  // we are below the lowest in or above the highest in
    43   for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
     52  for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
    4453
    4554    yout[0].elements.Flt[i] = NAN;
    … …  
    6675
    6776  // fill in the rest with NANs
     77  int NinLast = xin[0].Nelements - 1;
    6878  while (i < yout[0].Nelements) {
    69     yout[0].elements.Flt[i] = NAN;
     79    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN;
    7080    i++;
    7181  }
    7282
    7383  return (TRUE);
    74    
    7584}

  • trunk/Ohana/src/opihi/cmd.data/list_vectors.c

    r41340 r41341  
    3232  }
    3333
    34   if (vec->type == OPIHI_FLT) {
    35     if (Variable) {
    36       set_str_variable (Variable, "FLT");
    37     } else {
    38       gprint (GP_LOG, "%s : FLT\n", argv[1]);
    39     }
    40   } else {
    41     if (Variable) {
    42       set_str_variable (Variable, "INT");
    43     } else {
    44       gprint (GP_LOG, "%s : INT\n", argv[1]);
    45     }
     34  switch (vec->type) {
     35    case OPIHI_FLT:
     36      if (Variable) {
     37        set_str_variable (Variable, "FLT");
     38      } else {
     39        gprint (GP_LOG, "%s : FLT\n", argv[1]);
     40      }
     41      break;
     42    case OPIHI_INT:
     43      if (Variable) {
     44        set_str_variable (Variable, "INT");
     45      } else {
     46        gprint (GP_LOG, "%s : INT\n", argv[1]);
     47      }
     48      break;
     49    case OPIHI_STR:
     50      if (Variable) {
     51        set_str_variable (Variable, "STR");
     52      } else {
     53        gprint (GP_LOG, "%s : STR\n", argv[1]);
     54      }
     55      break;
    4656  }
    4757  if (Variable) free (Variable);

  • trunk/Ohana/src/opihi/cmd.data/medimage_commands.c

    r41340 r41341  
    1515int medimage_add (int argc, char **argv) {
    1616
     17  int N;
    1718  Buffer *image;
     19  Buffer *var = NULL;
     20
     21  if ((N = get_argument (argc, argv, "-variance"))) {
     22    remove_argument (N, &argc, argv);
     23    if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     24    remove_argument (N, &argc, argv);
     25  }
    1826
    1927  if (argc != 3) {
    20     gprint (GP_ERR, "USAGE: medimage add (name) (image)\n");
     28    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-variance variance]\n");
    2129    gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
     30    gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
    2231    return FALSE;
    2332  }
    2433
    2534  if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     35
     36  if (var) {
     37    if ((var->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
     38        (var->matrix.Naxis[1] != image->matrix.Naxis[1])) {
     39      gprint (GP_ERR, "variance buffer does not match image buffer dimensions\n");
     40      return FALSE;
     41    }
     42  }
    2643
    2744  MedImageType *median = FindMedImage (argv[1]);
    … …  
    4663  // new image should match existing medimage dimensions
    4764
     65  // AddMedImage (median, image, var);
    4866  int Ninput = median->Ninput;
    4967  median->Ninput ++;
    50   REALLOCATE (median->buffers, float *, median->Ninput);
     68  REALLOCATE (median->flx, float *, median->Ninput);
     69  REALLOCATE (median->var, float *, median->Ninput);
    5170
    52   ALLOCATE (median->buffers[Ninput], float, median->Nx*median->Ny);
    53   memcpy (median->buffers[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
     71  ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
     72  memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    5473
    55   return TRUE;
    56 }
    57 
    58 int medimage_calc (int argc, char **argv) {
    59 
    60   int ix, iy, n, N;
    61   Buffer *output;
    62 
    63   int CALC_MEAN = FALSE;
    64   if ((N = get_argument (argc, argv, "-mean"))) {
    65     CALC_MEAN = TRUE;
    66     remove_argument (N, &argc, argv);
    67   }
    68 
    69   if (argc != 3) {
    70     gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean]\n");
    71     gprint (GP_ERR, "       calculate the median image for the median image set\n");
    72     return FALSE;
    73   }
    74 
    75   MedImageType *median = FindMedImage (argv[1]);
    76   if (!median) {
    77     gprint (GP_ERR, "median image %s not found\n", argv[1]);
    78     return FALSE;
    79   }
    80 
    81   if ((output = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    82 
    83   int Ninput = median->Ninput;
    84   int Nx = median->Nx;
    85   int Ny = median->Ny;
    86 
    87   float *value = NULL;
    88   ALLOCATE (value, float, Ninput);
    89 
    90   gfits_free_matrix (&output->matrix);
    91   gfits_free_header (&output->header);
    92   if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    93 
    94   float *outvalue = (float *) output->matrix.buffer;
    95 
    96   for (iy = 0; iy < Ny; iy++) {
    97     for (ix = 0; ix < Nx; ix++) {
    98 
    99       int N = 0;
    100       int Npix = ix + Nx*iy;
    101       for (n = 0; n < Ninput; n++) {
    102         float v = median->buffers[n][Npix];
    103         if (!isfinite(v)) continue;
    104         value[N] = v;
    105         N++;
    106       }
    107       if (N == 0) continue;
    108 
    109       if (CALC_MEAN) {
    110         float sum = 0.0;
    111         for (n = 0; n < N; n++) {
    112           sum += value[n];
    113         }
    114         outvalue[Npix] = sum / (float) N;
    115       } else {
    116         fsort (value, N);
    117         outvalue[Npix] = value[(int)(0.5*N)];
    118       }
    119     }
     74  median->var[Ninput] = NULL;
     75  if (var) {
     76    ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
     77    memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    12078  }
    12179
    … …  
    12482
    12583/*
    126 int medimage_save (int argc, char **argv) {
     84   int medimage_save (int argc, char **argv) {
    12785
    128   int N;
     86   int N;
    12987
    130   int APPEND = FALSE;
    131   if ((N = get_argument (argc, argv, "-append"))) {
    132     APPEND = TRUE;
    133     remove_argument (N, &argc, argv);
    134   }
     88   int APPEND = FALSE;
     89   if ((N = get_argument (argc, argv, "-append"))) {
     90   APPEND = TRUE;
     91   remove_argument (N, &argc, argv);
     92   }
    13593
    136   if (argc != 3) {
    137     gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
    138     return FALSE;
    139   }
     94   if (argc != 3) {
     95   gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
     96   return FALSE;
     97   }
    14098
    141   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
    142     gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
    143     return (FALSE);
    144   }
    145   return TRUE;
    146 }
     99   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
     100   gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
     101   return (FALSE);
     102   }
     103   return TRUE;
     104   }
    147105
    148 int medimage_load (int argc, char **argv) {
     106   int medimage_load (int argc, char **argv) {
    149107
    150   if (argc != 3) {
    151     gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
    152     return FALSE;
    153   }
     108   if (argc != 3) {
     109   gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
     110   return FALSE;
     111   }
    154112
    155   if (!LoadMedImage(argv[2], argv[1])) {
    156     gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
    157     return (FALSE);
    158   }
    159   return TRUE;
    160 }
     113   if (!LoadMedImage(argv[2], argv[1])) {
     114   gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
     115   return (FALSE);
     116   }
     117   return TRUE;
     118   }
    161119*/
    162120

  • trunk/Ohana/src/opihi/cmd.data/peak.c

    r41340 r41341  
    3939      if (*X < start) continue;
    4040      if (*X > end) continue;
     41      if (!isfinite(*Y)) continue;
    4142      if (*Y < ymax) continue;
    4243      xmax = *X;
    … …  
    6364      if (*X < start) continue;
    6465      if (*X > end) continue;
     66      if (!isfinite(*Y)) continue;
    6567      if (*Y < ymax) continue;
    6668      xmax = *X;

  • trunk/Ohana/src/opihi/cmd.data/print_vectors.c

    r41340 r41341  
    5656        if (vec[i][0].type == OPIHI_FLT) {
    5757          gprint (GP_LOG, "%f ", vec[i][0].elements.Flt[j]);
    58         } else {
     58        }
     59        if (vec[i][0].type == OPIHI_INT) {
    5960          gprint (GP_LOG, OPIHI_INT_FMT" ", vec[i][0].elements.Int[j]);
    60         }
     61        }
     62        if (vec[i][0].type == OPIHI_STR) {
     63          gprint (GP_LOG, "%s ", vec[i][0].elements.Str[j]);
     64        }
    6165      }
    6266    }

  • trunk/Ohana/src/opihi/cmd.data/read_vectors.c

    r41340 r41341  
    3535
    3636// vector types
    37 enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_CHAR};
     37enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_DATE, COLTYPE_CHAR, COLTYPE_STR, COLTYPE_HMS};
    3838
    3939static int      Nvec     = 0;
    … …  
    9191    gprint (GP_ERR, "     -skip N : skip N lines before reading\n");
    9292    gprint (GP_ERR, "     column names may include a type: name:type\n");
    93     gprint (GP_ERR, "       type is int, float, char, time\n");
     93    gprint (GP_ERR, "       type is int, float, char, str, time, date, hms\n");
    9494    gprint (GP_ERR, "       for char, values are placed into a list $name:0 - $name:n\n");
     95    gprint (GP_ERR, "       for str, values are placed into a string-typed vector\n");
     96    gprint (GP_ERR, "       for hms, values are sexigesimal values HH:MM:SS (good for degrees as well)\n");
     97    gprint (GP_ERR, "       for date, values are human-readable date strings YYYY/MM/DD\n");
    9598    gprint (GP_ERR, "       for time, values are human-readable date/time strings YYYY/MM/DD,hh:mm:ss\n");
    96     gprint (GP_ERR, "       the resulting vector is a float based on the TIMEFORMAT, TIMEREF values\n");
     99    gprint (GP_ERR, "         date & time values are stored as a float using the TIMEFORMAT, TIMEREF values for the conversion\n");
    97100
    98101    gprint (GP_ERR, "     -fits options:\n");
    … …  
    152155      if (!strcasecmp(ptr, "int"))   { coltype[i] = COLTYPE_INT; }
    153156      if (!strcasecmp(ptr, "char"))  { coltype[i] = COLTYPE_CHAR; }
     157      if (!strcasecmp(ptr, "str"))   { coltype[i] = COLTYPE_STR; }
     158      if (!strcasecmp(ptr, "date"))  { coltype[i] = COLTYPE_DATE; }
    154159      if (!strcasecmp(ptr, "time"))  { coltype[i] = COLTYPE_TIME; }
     160      if (!strcasecmp(ptr, "hms"))   { coltype[i] = COLTYPE_HMS; }
    155161      if (!coltype[i]) goto bad_colname;
    156162    }
    … …  
    210216        break;
    211217      case COLTYPE_FLT:
     218      case COLTYPE_HMS:
    212219      case COLTYPE_TIME:
    213         // note that COLTYPE_TIME is a type of float
     220      case COLTYPE_DATE:
    214221        ResetVector (vec[i], OPIHI_FLT, NELEM);
     222        break;
     223      case COLTYPE_STR:
     224        ResetVector (vec[i], OPIHI_STR, NELEM);
    215225        break;
    216226      case COLTYPE_CHAR:
    … …  
    287297        time_t tvalue;
    288298        int readStatus = FALSE;
     299        int dataStatus = FALSE;
    289300        // need to make the if cases for coltype[i]
    290301        switch (coltype[i]) {
    … …  
    311322              break;
    312323            }
     324          case COLTYPE_STR:
     325            {
     326              // I need to get an isolated word in 'value' with the string value of this field
     327              char *ptr = IsCSV ? ptrparse_csv (col[i], c0) : ptrparse (col[i], c0);
     328              char *value = NULL;
     329              if (IsCSV) {
     330                char *end = parse_nextword_csv (ptr);
     331                if (end) {
     332                  value = end ? strncreate (ptr, end - ptr) : strcreate (ptr);
     333                }
     334              } else {
     335                value = thisword(ptr);
     336              }
     337              vec[i][0].elements.Str[Nelem] = value; // needs to be freed later.
     338              break;
     339            }
    313340          case COLTYPE_FLT:
    314341            readStatus = IsCSV ? dparse_csv (&dvalue, col[i], c0) : dparse (&dvalue, col[i], c0);
    … …  
    320347            vec[i][0].elements.Flt[Nelem] = readStatus ? dvalue : NAN;
    321348            break;
     349          case COLTYPE_DATE: {
     350            char *string = NULL;
     351            if (IsCSV) {
     352              char *ptr = ptrparse_csv (col[i], c0);
     353              string = strcreate (ptr); // create separate string (NULL-safe)
     354              ptr = (string == NULL) ? NULL : strchr (string, ',');
     355              if (ptr != NULL) *ptr = 0; // place an EOL here so parsing does not go past the comma
     356            } else {
     357              char *ptr = ptrparse (col[i], c0); // NULL-safe
     358              string = getword (ptr); // NULL-safe
     359            }         
     360            tvalue = ohana_date_to_sec (string); // NULL-safe (returns 0)
     361            dvalue = TimeValue (tvalue, TimeReference, TimeFormat);
     362            vec[i][0].elements.Flt[Nelem] = string ? dvalue : NAN;
     363            FREE (string);
     364            break;
     365          }
     366          case COLTYPE_HMS: {
     367            char *ptr = IsCSV ? ptrparse_csv (col[i], c0) : ptrparse (col[i], c0);
     368            char *string = strcreate (ptr); // make a copy so we can munge it in ohana_dms_to_ddd
     369            dataStatus = ohana_dms_to_ddd (&dvalue, string);
     370            vec[i][0].elements.Flt[Nelem] = ((string != NULL) && dataStatus) ? dvalue : NAN;
     371            FREE (string);
     372            break;
     373          }
    322374        }
    323375        if (!readStatus && VERBOSE) {
    … …  
    345397              break;
    346398            case COLTYPE_FLT:
     399            case COLTYPE_HMS:
    347400            case COLTYPE_TIME:
     401            case COLTYPE_DATE:
    348402              ResetVector (vec[i], OPIHI_FLT, NELEM);
     403              break;
     404            case COLTYPE_STR:
     405              ResetVector (vec[i], OPIHI_STR, NELEM);
    349406              break;
    350407            case COLTYPE_CHAR:
    … …  
    366423        break;
    367424      case COLTYPE_FLT:
     425      case COLTYPE_HMS:
    368426      case COLTYPE_TIME:
     427      case COLTYPE_DATE:
    369428        ResetVector (vec[i], OPIHI_FLT, Nelem);
     429        break;
     430      case COLTYPE_STR:
     431        ResetVector (vec[i], OPIHI_STR, Nelem);
    370432        break;
    371433      case COLTYPE_CHAR:
    … …  
    469531  }
    470532
    471   // if CharAsVectors, char fields will be saved as vectors NAME:0 -- NAME:n for n characters
    472   // else char fields will be saved as $NAME:0 - $NAME:m for m rows
     533  // by default, we now store a string-type field as a string-type vector.
     534
     535  // if CharAsList is selected (and My < 10000), char fields will be saved as $NAME:0 - $NAME:m for m rows
     536
     537  // if CharAsVectors is selected, char fields will be saved as vectors NAME:0 -- NAME:n for n characters
     538
    473539  // if (Ny > 10000), force CharAsVectors
    474540  int CharAsVectors = FALSE;
    … …  
    476542    remove_argument (N, &argc, argv);
    477543    CharAsVectors = TRUE;
     544  }
     545  int CharAsList = FALSE;
     546  if ((N = get_argument (argc, argv, "-char-list"))) {
     547    remove_argument (N, &argc, argv);
     548    CharAsList = TRUE;
    478549  }
    479550
    … …  
    651722    if (!FITS_TRANSPOSE) {
    652723      // read string column into a list rather than a vector
    653       if (!strcmp (type, "char") && !CharAsVectors && (Ny < 3000)) {
    654         char *fieldName = argv[i];
    655         char *Ptr = data;
    656         char varname[1024];  // used as a buffer for the names of string fields
    657         for (j = 0; j < Ny; j++) {
    658           set_list_varname (varname, fieldName, j, FALSE);
    659           char *value = strncreate (&Ptr[j*Nval], Nval);
    660           // replace instances of $ with _
    661           char *p = strchr (value, '$');
    662           while (p) {
    663             *p = '_';
    664             p = strchr (p, '$');
     724      if (!strcmp (type, "char")) {
     725        // save char-type field as a List:
     726        if (CharAsList && (Ny < 3000)) {
     727          char *fieldName = argv[i];
     728          char *Ptr = data;
     729          char varname[1024];  // used as a buffer for the names of string fields
     730          for (j = 0; j < Ny; j++) {
     731            set_list_varname (varname, fieldName, j, FALSE);
     732            char *value = strncreate (&Ptr[j*Nval], Nval);
     733            // replace instances of $ with _
     734            char *p = strchr (value, '$');
     735            while (p) {
     736              *p = '_';
     737              p = strchr (p, '$');
     738            }
     739            set_str_variable (varname, value);
     740            free (value);
    665741          }
    666           set_str_variable (varname, value);
    667           free (value);
     742          sprintf (varname, "%s:n", fieldName);
     743          set_int_variable (varname, Ny);
     744          continue;
    668745        }
    669         sprintf (varname, "%s:n", fieldName);
    670         set_int_variable (varname, Ny);
    671         continue;
     746        // save char-type field as a string-type vector:
     747        if (!CharAsVectors) {
     748          Vector *myVector = NULL;
     749          sprintf (name, "%s", argv[i]);
     750          if ((myVector = SelectVector (name, ANYVECTOR, TRUE)) == NULL) ESCAPE ("bad vector name\n");
     751          ResetVector (myVector, OPIHI_STR, Ny);
     752
     753          char *Ptr = data;
     754          for (j = 0; j < Ny; j++) {
     755            myVector[0].elements.Str[j] = strncreate (&Ptr[j*Nval], Nval);
     756            // replace instances of $ with _
     757            char *p = strchr (myVector[0].elements.Str[j], '$');
     758            while (p) { *p = '_'; p = strchr (p, '$'); }
     759          }
     760          continue;
     761        }
    672762      }
    673763
    674764      // define the multifield vector names (Nval vectors x Ny elements)
     765      // CharAsVectors is handled below automatically
    675766      ALLOCATE (vec, Vector *, Nval);
    676767      for (j = 0; j < Nval; j++) {

  • trunk/Ohana/src/opihi/cmd.data/roll.c

    r41340 r41341  
    33int roll (int argc, char **argv) {
    44 
    5   int Nbytes, Nextra;
    6   int dX, dx, dy, nx, ny;
    75  Buffer *buf;
     6
     7  // this function is probably not needed
     8  gprint (GP_ERR, "ERROR: use 'shift' instead of 'roll'\n");
     9  return (FALSE);
    810
    911  if (argc != 4) {
    … …  
    1416  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    1517
    16   dx = atof (argv[2]);
    17   dy = atof (argv[3]);
    18   if (dy != 0) {
    19     gprint (GP_ERR, "only x rolls implemented for now\n");
     18  int dx = atoi (argv[2]);
     19  int dy = atoi (argv[3]);
     20
     21  /* if (dx,dy < 0), we are moving the start position backwards,
     22     if (dx,dy > 0), we are moving the start position forward */
     23 
     24  int Nx = buf[0].matrix.Naxis[0];
     25  int Ny = buf[0].matrix.Naxis[1];
     26
     27  if (dy == 0) {
     28    // moves in just dx can be done row-by-row
     29
     30    int dX = abs(dx);
     31    int Nbytes = (Nx * Ny - dX) * sizeof (float);
     32    int Nextra = dX * sizeof (float);
     33
     34    if (dx < 0) {
     35      memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes);
     36      memset (&buf[0].matrix.buffer[Nbytes], 0, Nextra);
     37    } else {
     38      memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes);
     39      memset (buf[0].matrix.buffer, 0, Nextra);
     40    }
     41
     42    return TRUE;
    2043  }
    2144
    22   /* if (dx < 0), we are moving the start position back by dx pixels,
    23      if (dx > 0), we are moving the start position forward by dx pixels */
    24  
    25   nx = buf[0].matrix.Naxis[0];
    26   ny = buf[0].matrix.Naxis[1];
     45  if (dx == 0) {
     46    // moves in just dy can be done row-by-row
    2747
    28   dX = abs(dx);
    29   Nbytes = nx * ny * sizeof (float);
    30   Nextra = (nx * ny + dX) * sizeof (float);
     48    ALLOCATE_PTR (output, float, Nx*Ny);
    3149
    32   REALLOCATE (buf[0].matrix.buffer, char, Nextra);
     50    float *input = (float *) buf[0].matrix.buffer;
    3351
    34   if (dx < 0) {
    35     memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes);
    36   } else {
    37     memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes);
     52    for (int iy = 0; iy < Ny; iy++) {
     53      if (iy + dy < 0) continue;
     54      if (iy + dy >= Ny) continue;
     55      memcpy (&output[iy*Nx], &input[(iy + dy)*Nx], Nx);
     56    }
     57    if (dy > 0) {
     58      memset (output, 0, dy*Nx*sizeof(float));
     59    } else {
     60      int Nbytes = Nx * (Ny - abs(dy)) * sizeof (float);
     61      int Nextra = abs(dy) * Nx * sizeof (float);
     62      memset (&output[Nbytes], 0, Nextra);
     63    }
     64    free (buf[0].matrix.buffer);
     65    buf[0].matrix.buffer = (char *) output;
     66    return TRUE;
    3867  }
    3968
    40   REALLOCATE (buf[0].matrix.buffer, char, Nbytes);
     69  ALLOCATE_PTR (output, float, Nx*Ny);
     70  float *input = (float *) buf[0].matrix.buffer;
    4171
    42   return (TRUE);
    43 
     72  for (int iy = 0; iy < Ny; iy++) {
     73    if (iy + dy < 0) continue;
     74    if (iy + dy >= Ny) continue;
     75    for (int ix = 0; ix < Nx; ix++) {
     76      if (ix + dx < 0) continue;
     77      if (ix + dx >= Nx) continue;
     78      output[ix + iy*Nx] = input[(ix + dx) + (iy + dy)*Nx];
     79    }
     80  }
     81  free (buf[0].matrix.buffer);
     82  buf[0].matrix.buffer = (char *) output;
     83  return TRUE;
    4484}
    4585
     86

  • trunk/Ohana/src/opihi/cmd.data/rotate.c

    r41340 r41341  
    33int rotate (int argc, char **argv) {
    44 
    5   int i, j, NX, NY, X, Y, Lx, Ly, N;
     5  int i, j, NX, NY;
    66  float *in_buff, *out_buff, *c;
    7   double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, x, y, X1, Y1;
     7  double Xo, Yo, x, y, X1, Y1;
    88  double pc11, pc12, pc21, pc22, PC11, PC12, PC21, PC22;
    99  Buffer *buf;
    1010
    11   Xo = 0;
    12   Yo = 0;
    13   if ((N = get_argument (argc, argv, "-center"))) {
    14     remove_argument (N, &argc, argv);
    15     Xo  = atof(argv[N]);
    16     remove_argument (N, &argc, argv);
    17     Yo  = atof(argv[N]);
    18     remove_argument (N, &argc, argv);
    19   }
     11//  Xo = 0;
     12//  Yo = 0;
     13//  if ((N = get_argument (argc, argv, "-center"))) {
     14//    remove_argument (N, &argc, argv);
     15//    Xo  = atof(argv[N]);
     16//    remove_argument (N, &argc, argv);
     17//    Yo  = atof(argv[N]);
     18//    remove_argument (N, &argc, argv);
     19//  }
    2020
    2121  if (argc != 3) {
    … …  
    221221  }
    222222
    223   angle = atof (argv[2]);
    224   CosAngle = cos (angle*RAD_DEG);
    225   SinAngle = sin (angle*RAD_DEG);
     223  double angle = atof (argv[2]);
     224  double CosAngle = cos (angle*RAD_DEG);
     225  double SinAngle = sin (angle*RAD_DEG);
    226226 
    227227  gprint (GP_ERR, "rotating: %f %f %f\n", angle, CosAngle, SinAngle);
    228228
    229   Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
    230   Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
    231   dX = MAX(0,NY*SinAngle);
    232   dY = MAX(0,-NX*SinAngle);
    233   /*
    234   gprint (GP_ERR, "%f %f  -->  ", Xo, Yo);
    235   X1 = Xo*CosAngle - Yo*SinAngle + dX;
    236   Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
    237   gprint (GP_ERR, "%f %f\n", X1, Y1);
    238   */
     229  // we are rotating about the center pixel, (NX/2, NY/2),
     230  // but we are then putting the result in a new image
     231  // of size (Lx,Ly).
     232
     233  // (x,y) = (i - Nx/2),(j - Ny/2)
     234  // (x',y') = R(theta) (x,y)
     235  // (I,J) = (x' + Lx/2),(y' + Ly/2)
     236
     237  int Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
     238  int Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
    239239
    240240  /* fix reference pixel */
    241241  gfits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
    242242  gfits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
    243   /*
    244   X1 = (Xo - dX)*CosAngle - (Yo - dY)*SinAngle;
    245   Y1 = (dX - Xo)*SinAngle + (Yo - dY)*CosAngle;
    246   */
    247   X1 = Xo*CosAngle - Yo*SinAngle + dX;
    248   Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
     243  X1 = (Xo - NX/2)*CosAngle + (Yo - NY/2)*SinAngle + Lx/2;
     244  Y1 = (NX/2 - Xo)*SinAngle + (Yo - NY/2)*CosAngle + Ly/2;
    249245  gfits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    250246  gfits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    … …  
    255251  gfits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    256252  gfits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    257   PC11 = pc11*CosAngle - pc21*SinAngle;
    258   PC12 = pc12*CosAngle - pc22*SinAngle;
    259   PC21 = pc21*CosAngle + pc11*SinAngle;
    260   PC22 = pc22*CosAngle + pc12*SinAngle;
     253  PC11 = pc11*CosAngle + pc21*SinAngle;
     254  PC12 = pc12*CosAngle + pc22*SinAngle;
     255  PC21 = pc21*CosAngle - pc11*SinAngle;
     256  PC22 = pc22*CosAngle - pc12*SinAngle;
    261257  gfits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    262258  gfits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    … …  
    270266  gfits_create_matrix (&buf[0].header, &buf[0].matrix);
    271267  gfits_print_alt (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
     268
    272269  out_buff = (float *)buf[0].matrix.buffer;
     270
    273271  for (j = 0; j < Ly; j++) {
    274272    for (i = 0; i < Lx; i++, out_buff++) {
    275273
    276       x = (i - dX)*CosAngle + (j - dY)*SinAngle;
    277       y = (dX - i)*SinAngle + (j - dY)*CosAngle;
    278       X = (int) x;
    279       Y = (int) y;
    280 
    281       if (X < 0) continue;
    282       if (X >= NX - 1) continue;
    283       if (Y < 0) continue;
    284       if (Y >= NY - 1) continue;
    285 
    286       c = &in_buff[X + NX*Y];
    287       fx = x - X;
    288       fy = y - Y;
     274      float xo = (i - Lx/2);
     275      float yo = (j - Ly/2);
     276
     277      x =  xo*CosAngle + yo*SinAngle;
     278      y = -xo*SinAngle + yo*CosAngle;
     279
     280      int I = x + NX/2;
     281      int J = y + NY/2;
     282
     283      if (I < 0) continue;
     284      if (I >= NX - 1) continue;
     285      if (J < 0) continue;
     286      if (J >= NY - 1) continue;
     287
     288      c = &in_buff[I + NX*J];
     289
     290      int X = (int) x;
     291      int Y = (int) y;
     292      double fx = x - X;
     293      double fy = y - Y;
     294
    289295      *out_buff = (c[0]*(1-fx) + c[1]*fx)*(1-fy) + (c[NX+1]*fx + c[NX]*(1-fx))*fy;
     296//    *out_buff = c[0];
    290297    }
    291298  }

  • trunk/Ohana/src/opihi/cmd.data/spline.c

    r41340 r41341  
    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;

  • trunk/Ohana/src/opihi/cmd.data/spline_commands.c

    r41340 r41341  
    2828  if (xvec->Nelements != yvec->Nelements) {
    2929    gprint (GP_ERR, "x and y vectors must be the same length\n");
     30    return FALSE;
     31  }
     32
     33  if (xvec->Nelements < 3) {
     34    gprint (GP_ERR, "cannot make a spline with fewer than 3 knots\n");
    3035    return FALSE;
    3136  }
    … …  
    9398}
    9499
     100int spline_print (int argc, char **argv) {
     101
     102  int i;
     103
     104  if (argc != 2) {
     105    gprint (GP_ERR, "USAGE: spline print (name)\n");
     106    return FALSE;
     107  }
     108
     109  Spline *myspline = FindSpline (argv[1]);
     110  if (!myspline) {
     111    gprint (GP_ERR, "spline %s not found\n", argv[1]);
     112    return (FALSE);
     113  }
     114
     115  for (i = 0; i < myspline->Nknots; i++) {
     116    gprint (GP_LOG, "%lf : %lf : %lf\n", myspline->xk[i], myspline->yk[i], myspline->y2[i]);
     117  }   
     118
     119  return TRUE;
     120}
     121
    95122int spline_load (int argc, char **argv) {
    96123
    … …  
    109136int spline_delete (int argc, char **argv) {
    110137
    111   int status;
    112   Spline *spline;
     138  int N, status;
     139
     140  int QUIET = FALSE;
     141  if ((N = get_argument (argc, argv, "-q"))) {
     142    QUIET = TRUE;
     143    remove_argument (N, &argc, argv);
     144  }
    113145
    114146  if (argc != 2) {
    … …  
    117149  }
    118150
    119   spline = FindSpline (argv[1]);
     151  Spline *spline = FindSpline (argv[1]);
    120152  if (spline == NULL) {
     153    if (QUIET) return TRUE;
    121154    gprint (GP_ERR, "spline %s not found\n", argv[1]);
    122155    return FALSE;
    … …  
    130163int spline_rename (int argc, char **argv) {
    131164
    132   Spline *spline;
    133 
    134165  if (argc != 3) {
    135166    gprint (GP_ERR, "USAGE: spline rename (spline) (newname)\n");
    … …  
    137168  }
    138169
    139   spline = FindSpline (argv[1]);
     170  Spline *spline = FindSpline (argv[1]);
    140171  if (spline == NULL) {
    141172    gprint (GP_ERR, "spline %s not found\n", argv[1]);
    142173    return FALSE;
     174  }
     175
     176  Spline *newname = FindSpline (argv[2]);
     177  if (newname != NULL) {
     178    DeleteSpline (newname);
    143179  }
    144180
    … …  
    147183  return TRUE;
    148184}
    149 
    150 /*
    151 int spline_listspline (int argc, char **argv) {
    152 
    153   Spline *spline;
    154 
    155   if (argc != 2) {
    156     gprint (GP_ERR, "USAGE: spline listspline (spline)\n");
    157     return FALSE;
    158   }
    159 
    160   spline = FindSpline (argv[1]);
    161   if (spline == NULL) {
    162     gprint (GP_ERR, "spline %s not found\n", argv[1]);
    163     return FALSE;
    164   }
    165 
    166   ListPages (spline);
    167   return TRUE;
    168 }
    169 */
    170 

  • trunk/Ohana/src/opihi/cmd.data/test/cut.sh

    r41340 r41341  
    44 memtest1
    55end
     6
     7macro test2
     8
     9  mgaussdev z 50 100 0.0 1.0
     10
     11  cut z x y_sum    x 0 0 z[][0] z[0][]
     12  cut z x y_median x 0 0 z[][0] z[0][] -median
     13  cut z x y_mean   x 0 0 z[][0] z[0][] -mean
     14  cut z x y_inner  x 0 0 z[][0] z[0][] -inner
     15 
     16end
     17
     18macro checkrange
     19  if ($0 != 2)
     20    echo "USAGE: checkrange (Npts)"
     21    break
     22  end
     23
     24  $Npts = $1
     25  if ($Npts % 2)
     26    $Ncenter = int($Npts / 2)
     27    $Nquarter = int(0.25 * $Npts)
     28    $Ns = $Ncenter - $Nquarter
     29    $Ne = $Ncenter + $Nquarter
     30  else
     31    $Ncenter = int($Npts / 2) - 1
     32    $Nquarter = int(0.25 * $Npts)
     33    $Ns = $Ncenter - $Nquarter
     34    $Ne = $Ncenter + $Nquarter + 1
     35  end
     36  if ($Ns < 0)
     37    echo "error: Ns < 0: $Ns $Ncenter $Nquarter"
     38  end
     39  if ($Ne >= $Npts)
     40    echo "error: Ne >= Npts: $Ne $Ncenter $Nquarter"
     41  end
     42  echo "$Npts : $Ncenter : $Nquarter : $Ns $Ne"
     43  for i 0 $Npts
     44    echo $i {($i >= $Ns) && ($i <= $Ne)}
     45  end
     46end
     47
    648
    749# Test if cut works
    … …  
    3779 end
    3880end
     81
     82macro memtest1
     83
     84 $i = 0
     85 mcreate tim 100 100
     86 cut tim xdir imx X 40 4 60 6
     87 # do one to set up memory that should stay used
     88
     89 memory check
     90 for i 0 100
     91   cut tim xdir imx X 40 4 60 6
     92 end
     93 memory check
     94   
     95 for i 0 100
     96   cut tim xdir imx y 40 4 60 6
     97 end
     98 memory check
     99   
     100 for i 0 100
     101   cut -median tim xdir imx X 40 4 60 6
     102 end
     103 memory check
     104   
     105 for i 0 100
     106   cut -median tim xdir imx y 40 4 60 6
     107 end
     108 memory check
     109   
     110 for i 0 100
     111   cut -mean tim xdir imx X 40 4 60 6
     112 end
     113 memory check
     114   
     115 for i 0 100
     116   cut -mean tim xdir imx y 40 4 60 6
     117 end
     118 memory check
     119   
     120end

  • trunk/Ohana/src/opihi/cmd.data/test/delete.sh

    r41340 r41341  
    4949
    5050end
     51
     52# Test delete
     53macro test2
     54
     55 $PASS = 1
     56
     57 for i 0 1000
     58   create v 0 200
     59   delete v
     60 end
     61end

  • trunk/Ohana/src/opihi/cmd.data/test/medimage.sh

    r41340 r41341  
    11
    2 macro go
    3  mcreate a 30 30
    4  for i 0 40
    5   set a$i = zero(a) + $i
    6  end
    7 
    8  for i 0 40
    9   medimage add t1 a$i
    10  end
    11 end
     2macro test.mean
     3
     4 $Nsample = 16
     5 medimage delete -q t1
     6 for i 0 $Nsample
     7  mgaussdev t 50 50 0.0 1.0
     8  medimage add t1 t
     9  imhist -q t x y -range -10 10 -delta 0.1
     10
     11  $C0 = 0
     12  $C1 = 1.5
     13  $C2 = 400
     14  $C3 = 0
     15  vgauss -q x y con yf
     16  # echo $C1
     17 end
     18
     19 medimage calc t1 T -mean
     20
     21 imhist -q T x y -range -10 10 -delta 0.1
     22 lim -n 1 x y; clear; box; plot -x hist x y
     23 $C0 = 0
     24 $C1 = 1.5
     25 $C2 = 400
     26 $C3 = 0
     27 vgauss -q x y con yf
     28 echo "expect {1/sqrt($Nsample)} : $C1"
     29 plot -c red -x line x yf
     30end
     31
     32macro test.median
     33
     34 $Nsample = 16
     35 medimage delete -q t1
     36 for i 0 $Nsample
     37  mgaussdev t 50 50 0.0 1.0
     38  medimage add t1 t
     39  imhist -q t x y -range -10 10 -delta 0.1
     40
     41  $C0 = 0
     42  $C1 = 1.5
     43  $C2 = 400
     44  $C3 = 0
     45  vgauss -q x y con yf
     46 end
     47
     48 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     49 # (somewhat higher scatter)
     50 medimage calc t1 T
     51
     52 imhist -q T x y -range -10 10 -delta 0.1
     53 lim -n 1 x y; clear; box; plot -x hist x y
     54 $C0 = 0
     55 $C1 = 1.5
     56 $C2 = 400
     57 $C3 = 0
     58 vgauss -q x y con yf
     59 echo "expect {1/sqrt($Nsample)} : $C1 (actually should be a bit higher)"
     60 plot -c red -x line x yf
     61end
     62 
     63macro test.wtmean
     64
     65 $Nsample = 8
     66 $sig1 = 1.0
     67 $sig2 = 3.0
     68
     69 medimage delete -q t1
     70 for i 0 $Nsample
     71  mgaussdev t 50 50 0.0 $sig1
     72  set v = $sig1^2 + zero(t)     
     73  medimage add t1 t -variance v
     74  imhist -q t x y -range -10 10 -delta 0.1
     75
     76  $C0 = 0
     77  $C1 = 1.5
     78  $C2 = 400
     79  $C3 = 0
     80  vgauss -q x y con yf
     81  # echo $C1
     82 end
     83
     84 for i 0 $Nsample
     85  mgaussdev t 50 50 0.0 $sig2
     86  set v = $sig2^2 + zero(t)     
     87  medimage add t1 t -variance v
     88  imhist -q t x y -range -10 10 -delta 0.1
     89
     90  $C0 = 0
     91  $C1 = 1.5
     92  $C2 = 400
     93  $C3 = 0
     94  vgauss -q x y con yf
     95  # echo $C1
     96 end
     97
     98 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     99 # (somewhat higher scatter)
     100 medimage calc t1 T -wtmean
     101
     102 imhist -q T x y -range -10 10 -delta 0.1
     103 lim -n 1 x y; clear; box; plot -x hist x y
     104 $C0 = 0
     105 $C1 = 1.5
     106 $C2 = 400
     107 $C3 = 0
     108 vgauss -q x y con yf
     109 $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
     110 echo "expect {1/sqrt($S1)} : $C1"
     111 plot -c red -x line x yf
     112end
     113
     114macro test.irls
     115 medimage delete -q t1
     116 $Nsample = 16
     117 $sig = 1.0
     118 for i 0 $Nsample
     119  mgaussdev t 50 50 0.0 $sig
     120  set v = $sig^2 + zero(t)     
     121
     122  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     123  set ts = t + bad
     124
     125  medimage add t1 ts -variance v
     126  imhist -q t x y -range -10 10 -delta 0.1
     127
     128  $C0 = 0
     129  $C1 = 1.5
     130  $C2 = 400
     131  $C3 = 0
     132  vgauss -q x y con yf
     133  # echo $C1
     134 end
     135
     136 # get stats for straight mean:
     137 medimage calc t1 Tm -mean
     138
     139 imhist -q Tm x y -range -10 10 -delta 0.1
     140 lim -n 1 x y; clear; box; plot -x hist x y
     141 $C0 = 0
     142 $C1 = 1.5
     143 $C2 = 400
     144 $C3 = 0
     145 vgauss -q x y con yf
     146 echo "sigma from straight stdev: $C1"
     147 # stats Tm
     148 
     149 plot -c red -x line x yf
     150
     151 # get stats for irls
     152 medimage calc t1 Ti -irls
     153
     154 imhist -q Ti x y -range -10 10 -delta 0.1
     155 lim -n 2 x y; clear; box; plot -x hist x y
     156 $C0 = 0
     157 $C1 = 1.5
     158 $C2 = 400
     159 $C3 = 0
     160 vgauss -q x y con yf
     161 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
     162 # stats Ti
     163 
     164 plot -c red -x line x yf
     165end
     166
     167
     168###################33
     169
     170
     171macro test.mean.var
     172
     173 $Nsample = 64
     174 $sig = 2.0
     175 medimage delete -q t1
     176 for i 0 $Nsample
     177  mgaussdev t 100 100 0.0 $sig
     178  medimage add t1 t
     179  imhist -q t x y -range -10 10 -delta 0.1
     180
     181  $C0 = 0
     182  $C1 = 1.5
     183  $C2 = 400
     184  $C3 = 0
     185  vgauss -q x y con yf
     186  # echo $C1
     187 end
     188
     189 medimage calc t1 T -mean -variance Tv
     190
     191 imhist -q T x y -range -10 10 -delta 0.1
     192 lim -n 1 x y; clear; box; plot -x hist x y
     193 $C0 = 0
     194 $C1 = 1.5
     195 $C2 = 400
     196 $C3 = 0
     197 vgauss -q x y con yf
     198 plot -c red -x line x yf
     199
     200 imhist Tv xv yv -range -1 4 -delta 0.1
     201 lim -n 2 xv yv; clear; box; plot xv yv -x hist
     202
     203 stat Tv
     204 echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
     205end
     206
     207macro test.median.var
     208
     209 $Nsample = 64
     210 $sig = 2.0
     211 medimage delete -q t1
     212 for i 0 $Nsample
     213  mgaussdev t 50 50 0.0 $sig
     214  medimage add t1 t
     215  imhist -q t x y -range -10 10 -delta 0.1
     216
     217  $C0 = 0
     218  $C1 = 1.5
     219  $C2 = 400
     220  $C3 = 0
     221  vgauss -q x y con yf
     222  # echo $C1
     223 end
     224
     225 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     226 # (somewhat higher scatter)
     227 medimage calc t1 T -variance Tv
     228
     229 imhist -q T x y -range -10 10 -delta 0.1
     230 lim -n 1 x y; clear; box; plot -x hist x y
     231 $C0 = 0
     232 $C1 = 1.5
     233 $C2 = 400
     234 $C3 = 0
     235 vgauss -q x y con yf
     236 plot -c red -x line x yf
     237
     238 imhist Tv xv yv -range -1 4 -delta 0.1
     239 lim -n 2 xv yv; clear; box; plot xv yv -x hist
     240
     241 stat Tv
     242 echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
     243end
     244 
     245macro test.wtmean.var
     246
     247 $Nsample = 32
     248 $sig1 = 1.0
     249 $sig2 = 1.0
     250
     251 medimage delete -q t1
     252 for i 0 $Nsample
     253  mgaussdev t 50 50 0.0 $sig1
     254  set v = $sig1^2 + zero(t)     
     255  medimage add t1 t -variance v
     256  imhist -q t x y -range -10 10 -delta 0.1
     257
     258  $C0 = 0
     259  $C1 = 1.5
     260  $C2 = 400
     261  $C3 = 0
     262  vgauss -q x y con yf
     263  # echo $C1
     264 end
     265
     266 for i 0 $Nsample
     267  mgaussdev t 50 50 0.0 $sig2
     268  set v = $sig2^2 + zero(t)     
     269  medimage add t1 t -variance v
     270  imhist -q t x y -range -10 10 -delta 0.1
     271
     272  $C0 = 0
     273  $C1 = 1.5
     274  $C2 = 400
     275  $C3 = 0
     276  vgauss -q x y con yf
     277  # echo $C1
     278 end
     279
     280 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     281 # (somewhat higher scatter)
     282 medimage calc t1 T -wtmean -variance Tv
     283
     284 imhist -q T x y -range -10 10 -delta 0.1
     285 lim -n 1 x y; clear; box; plot -x hist x y
     286 $C0 = 0
     287 $C1 = 1.5
     288 $C2 = 400
     289 $C3 = 0
     290 vgauss -q x y con yf
     291 plot -c red -x line x yf
     292
     293 stat -q Tv
     294 $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
     295 echo $C1 vs {sqrt($MEDIAN)} : expect {1/sqrt($S1)}
     296end
     297
     298macro test.irls.var
     299
     300 $Nsample = 16
     301 $sig = 1.0
     302
     303 medimage delete -q t1
     304 for i 0 $Nsample
     305  mgaussdev t 50 50 0.0 $sig
     306  set v = $sig^2 + zero(t)     
     307
     308  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     309  set ts = t + bad
     310
     311  medimage add t1 ts -variance v
     312  imhist -q t x y -range -10 10 -delta 0.1
     313
     314  $C0 = 0
     315  $C1 = 1.5
     316  $C2 = 400
     317  $C3 = 0
     318  vgauss -q x y con yf
     319  # echo $C1
     320 end
     321
     322 # get stats for straight mean:
     323 medimage calc t1 Tm -mean -variance Tv
     324
     325 imhist -q Tm x y -range -10 10 -delta 0.1
     326 lim -n 1 x y; clear; box; plot -x hist x y
     327 $C0 = 0
     328 $C1 = 1.5
     329 $C2 = 400
     330 $C3 = 0
     331 vgauss -q x y con yf
     332 echo "sigma from straight stdev: $C1"
     333 # stats Tm
     334 
     335 plot -c red -x line x yf
     336
     337 # get stats for irls
     338 medimage calc t1 Ti -irls -variance Tv
     339
     340 imhist -q Ti x y -range -10 10 -delta 0.1
     341 lim -n 2 x y; clear; box; plot -x hist x y
     342 $C0 = 0
     343 $C1 = 1.5
     344 $C2 = 400
     345 $C3 = 0
     346 vgauss -q x y con yf
     347 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
     348 # stats Ti
     349 
     350 plot -c red -x line x yf
     351
     352 set dTv = sqrt(Tv)
     353 imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
     354
     355 stat -q Tv
     356 echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
     357end
     358
     359macro test.irls.boot.var
     360
     361 $Nsample = 64
     362 $sig = 1.0
     363
     364 medimage delete -q t1
     365 for i 0 $Nsample
     366  mgaussdev t 200 200 0.0 $sig
     367  set v = $sig^2 + zero(t)     
     368
     369  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     370  set ts = t + bad
     371
     372  mgaussdev noise 200 200 0.0 0.5
     373  set ts = ts + noise
     374
     375  medimage add t1 ts -variance v
     376  imhist -q t x y -range -10 10 -delta 0.1
     377
     378  $C0 = 0
     379  $C1 = 1.5
     380  $C2 = 400
     381  $C3 = 0
     382  vgauss -q x y con yf
     383  # echo $C1
     384 end
     385
     386 # get stats for straight mean:
     387 medimage calc t1 Tm -mean -variance Tv
     388
     389 imhist -q Tm x y -range -10 10 -delta 0.02
     390 lim -n 1 x y; clear; box; plot -x hist x y
     391 peak -q x y
     392 $C0 = $peakpos
     393 $C1 = 1.5*$sig / sqrt($Nsample)
     394 $C2 = $peakval
     395 $C3 = 0
     396 vgauss -q x y con yf
     397 echo "sigma from straight stdev: $C1"
     398 # stats Tm
     399 
     400 plot -c red -x line x yf
     401
     402 # get stats for irls
     403 date
     404 medimage calc t1 Ti -irls -variance Tv -bootstrap-iter 100
     405 date
     406
     407 imhist -q Ti x y -range -10 10 -delta 0.02
     408 lim -n 2 x y; clear; box; plot -x hist x y
     409 peak -q x y
     410 $C0 = $peakpos
     411 $C1 = 1.5*$sig / sqrt($Nsample)
     412 $C2 = $peakval
     413 $C3 = 0
     414 vgauss -q x y con yf
     415 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
     416 # stats Ti
     417 
     418 plot -c red -x line x yf
     419
     420 set dTv = sqrt(Tv)
     421 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
     422
     423 stat -q Tv
     424 echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
     425end
     426
     427##############################
     428macro test.irls.boot.test
     429
     430 $Nsample = 100
     431 $sig1 = 1.0
     432
     433 medimage delete -q t1
     434 for i 0 $Nsample
     435  mgaussdev t 100 100 0.0 $sig1
     436  set v = $sig1^2 + zero(t)     
     437
     438  medimage add t1 t -variance v
     439 end
     440
     441 # get stats for irls
     442 medimage calc t1 Ti -irls -variance Tv -bootstrap
     443
     444 imhist -q Ti x y -range {-10*$sig1/sqrt($Nsample)} {10*$sig1/sqrt($Nsample)} -delta 0.01
     445 lim -n 2 x y; clear; box; plot -x hist x y
     446 peak -q x y
     447 $C0 = $peakpos
     448 $C1 = 1.5*$sig1/sqrt($Nsample)
     449 $C2 = $peakval
     450 $C3 = 0
     451 vgauss x y con yf
     452 echo "sigma from irls: $C1 (ideal is {$sig1/sqrt($Nsample)})"
     453 # stats Ti
     454 
     455 plot -c red -x line x yf
     456
     457 set dTv = sqrt(Tv)
     458 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
     459
     460 stat -q irls_npt
     461 $Npix = $MEAN
     462
     463 stat -q Tv
     464 echo "sigma of irls average: $C1, sqrt(mean) of irls variance: {sqrt($MEAN)}, (ideal is {$sig1/sqrt($Npix)})"
     465end
     466
     467macro test.irls.range.var
     468 medimage delete -q t1
     469 for i 0 8
     470  mgaussdev t 50 50 0.0 1.0
     471  set v = 1.0 + zero(t)     
     472
     473  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     474  set ts = t + bad
     475
     476  medimage add t1 ts -variance v
     477  imhist -q t x y -range -10 10 -delta 0.1
     478
     479  $C0 = 0
     480  $C1 = 1.5
     481  $C2 = 400
     482  $C3 = 0
     483  vgauss -q x y con yf
     484  echo $C1
     485 end
     486
     487 for i 0 8
     488  mgaussdev t 50 50 0.0 3.0
     489  set v = 3.0 + zero(t)     
     490
     491  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     492  set ts = t + bad
     493
     494  medimage add t1 ts -variance v
     495  imhist -q t x y -range -10 10 -delta 0.1
     496
     497  $C0 = 0
     498  $C1 = 1.5
     499  $C2 = 400
     500  $C3 = 0
     501  vgauss -q x y con yf
     502  echo $C1
     503 end
     504
     505 # get stats for straight mean:
     506 medimage calc t1 Tm -mean -variance Tv
     507
     508 imhist -q Tm x y -range -10 10 -delta 0.1
     509 lim -n 1 x y; clear; box; plot -x hist x y
     510 $C0 = 0
     511 $C1 = 1.5
     512 $C2 = 400
     513 $C3 = 0
     514 vgauss -q x y con yf
     515 echo $C1
     516 stats Tm
     517 
     518 plot -c red -x line x yf
     519
     520 # get stats for irls
     521 medimage calc t1 Ti -irls -variance Tv
     522
     523 imhist -q Ti x y -range -10 10 -delta 0.1
     524 lim -n 2 x y; clear; box; plot -x hist x y
     525 $C0 = 0
     526 $C1 = 1.5
     527 $C2 = 400
     528 $C3 = 0
     529 vgauss -q x y con yf
     530 echo $C1
     531 stats Ti
     532 
     533 plot -c red -x line x yf
     534
     535 stat -q Tv
     536 echo $C1 vs {sqrt($MEDIAN)}
     537
     538 set dTv = sqrt(Tv)
     539 imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
     540end

  • trunk/Ohana/src/opihi/cmd.data/threshold.c

    r41340 r41341  
    2525    remove_argument (N, &argc, argv);
    2626    BinMax = atoi (argv[N]);
     27    remove_argument (N, &argc, argv);
     28  }
     29
     30  float ValMin = NAN;
     31  float ValMax = NAN;
     32  if ((N = get_argument (argc, argv, "-vrange"))) {
     33    remove_argument (N, &argc, argv);
     34    ValMin = atof (argv[N]);
     35    remove_argument (N, &argc, argv);
     36    ValMax = atof (argv[N]);
    2737    remove_argument (N, &argc, argv);
    2838  }
    … …  
    6474  }
    6575
     76  // ValMin, ValMax not compatible with BinMin, BinMax
     77  if (isfinite(ValMin)) {
     78    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMin);
     79    BinMin = (vecx[0].elements.Flt[binTest+1] == ValMin) ? binTest : binTest + 1;
     80    // XXX dangerous: check that this does not run off the end, etc
     81  }
     82  if (isfinite(ValMax)) {
     83    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMax);
     84    BinMax = (vecx[0].elements.Flt[binTest+1] == ValMax) ? binTest : binTest + 1;
     85    // XXX dangerous: check that this does not run off the end, etc
     86  }
     87
    6688  // use bisection to find the value
    6789
    … …  
    7193
    7294  if (!REVERSE) {
     95    // this algorithm assumes vecy[BinMax] > threshold, vecy[BinMin] < threshold
     96    if (vecy[0].elements.Flt[BinMin] > value) {
     97      gprint (GP_ERR, "ERROR: all values above threshold\n");
     98      return FALSE;
     99    }
     100    if (vecy[0].elements.Flt[BinMax] < value) {
     101      gprint (GP_ERR, "ERROR: all values below threshold\n");
     102      return FALSE;
     103    }
    73104    while (Nhi - Nlo > 10) {
    74105      N = 0.5*(Nlo + Nhi);
    … …  
    89120    }
    90121  } else {
     122    // this algorithm assumes vecy[BinMin] > threshold, vecy[BinMax] < threshold
     123    if (vecy[0].elements.Flt[BinMin] < value) {
     124      gprint (GP_ERR, "ERROR: all values below threshold\n");
     125      return FALSE;
     126    }
     127    if (vecy[0].elements.Flt[BinMax] > value) {
     128      gprint (GP_ERR, "ERROR: all values above threshold\n");
     129      return FALSE;
     130    }
    91131    while (Nhi - Nlo > 10) {
    92132      N = 0.5*(Nlo + Nhi);

  • trunk/Ohana/src/opihi/cmd.data/tvcolors.c

    r41340 r41341  
    3737    gprint (GP_ERR, "USAGE: tvcolors (colormap) [-nan red green blue]\n");
    3838    gprint (GP_ERR, " colormap options : greyscale, -greyscale, rainbow, heat, fullcolor, ruffcolor (also grayscale, -grayscale)\n");
    39     gprint (GP_ERR, " colormap file name of the for WORD:PATH where WORD = (file,csvf,lgcy,cetf)\n");
    40     gprint (GP_ERR, " lgcy:\n");
    41     gprint (GP_ERR, "   the colormap file must contain 4 columns: f R B G; each line defines a color transition.\n");
    42     gprint (GP_ERR, "   f: 0 - 1 defines the scale value for the transition\n");
    43     gprint (GP_ERR, "   R,B,G: 0 - 1 define the value of the color at the transition point\n");
    44     gprint (GP_ERR, " file:\n");
     39    gprint (GP_ERR, " colormap file: if colormap name is given as file:path/to/file, the colormap is read from the file\n");
    4540    gprint (GP_ERR, "   the colormap file must contain 4 columns: f R G B; each line defines a color transition.\n");
    4641    gprint (GP_ERR, "   f: 0 - 1 defines the scale value for the transition\n");
    4742    gprint (GP_ERR, "   R,G,B: 0 - 1 define the value of the color at the transition point\n");
    48     gprint (GP_ERR, " csvf:\n");
    49     gprint (GP_ERR, "   the colormap file must contain 4 columns: f,R,G,B separate by commas; each line defines a color transition.\n");
    50     gprint (GP_ERR, "   f: 0 - 1 defines the scale value for the transition\n");
    51     gprint (GP_ERR, "   R,G,B: 0 - 1 define the value of the color at the transition point\n");
    52     gprint (GP_ERR, " cetf:\n");
    53     gprint (GP_ERR, "   the colormap file must 256 rows contain 3 columns: R,G,B separate by commas; each line defines a color transition.\n");
    54     gprint (GP_ERR, "   R,G,B: 0 - 1 define the value of the color at the transition point\n");
    55     gprint (GP_ERR, "   CETF-format files can be found at : https://peterkovesi.com/projects/colourmaps\n");
    5643    return (FALSE);
    5744  }

  • trunk/Ohana/src/opihi/cmd.data/vgauss.c

    r41340 r41341  
    33/* local private functions */
    44opihi_flt fgaussOD (opihi_flt, opihi_flt *, int, opihi_flt *);
     5int vgauss_apply (int argc, char **argv);
    56
    67# define GET_VAR(V,A) \
    … …  
    2021  Vector *xvec, *yvec, *svec, *ovec;
    2122  char *c, name[16];
     23
     24  if ((N = get_argument (argc, argv, "-apply"))) {
     25    remove_argument (N, &argc, argv);
     26    int status = vgauss_apply (argc, argv);
     27    return status;
     28  }
    2229
    2330  Quiet = FALSE;
    … …  
    104111  } 
    105112  if (!Quiet) gprint (GP_ERR, "%d iterations\n", i);
     113  set_variable ("Chisq", chisq);
     114  set_variable ("ChisqNu", chisq / (float) (Npts - Npar));
    106115
    107116  for (i = 0; i < Npts; i++) {
    … …  
    125134}
    126135
     136int vgauss_apply (int argc, char **argv) {
     137
     138  opihi_flt par[4];
     139  Vector *xvec, *ovec;
     140  char *c;
     141
     142  // NOTE: -apply has already been stripped by vgauss call
     143  if (argc != 3) {
     144    gprint (GP_ERR, "USAGE: vgauss -apply <x> <out>\n");
     145    return (FALSE);
     146  }
     147 
     148  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     149  if ((ovec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     150
     151  CastVector (xvec, OPIHI_FLT);
     152
     153  GET_VAR (par[0], "C0");
     154  GET_VAR (par[1], "C1");
     155  GET_VAR (par[2], "C2");
     156  GET_VAR (par[3], "C3");
     157  int Npar = 4;
     158
     159  int Npts = xvec[0].Nelements;
     160  ResetVector (ovec, OPIHI_FLT, Npts);
     161
     162  for (int i = 0; i < Npts; i++) {
     163    ovec[0].elements.Flt[i] = fgaussOD (xvec[0].elements.Flt[i], par, Npar, NULL);
     164  }
     165  ovec[0].Nelements = Npts;
     166
     167  return TRUE;
     168}
     169
    127170/* pars: x_o, sigma, I, back */
    128171opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
    … …  
    135178  f = par[2]*r + par[3];
    136179
    137   dpar[0] = par[2]*r*z/par[1];
    138   dpar[1] = par[2]*r*z*z/par[1];
    139   dpar[2] = r;
    140   dpar[3] = 1;
     180  if (dpar) {
     181    dpar[0] = par[2]*r*z/par[1];
     182    dpar[1] = par[2]*r*z*z/par[1];
     183    dpar[2] = r;
     184    dpar[3] = 1;
     185  }
    141186 
    142187  return (f);

  • trunk/Ohana/src/opihi/cmd.data/write_vectors.c

    r41340 r41341  
    3232  }
    3333
    34   /* option generate a FITS output table */
     34  // option generate a FITS output table (FITS holds the filename)
     35  // in FITS output context, -header is interpretted as a buffer containing
     36  // header keywords to supplement the FITS table header
    3537  char *FITS = NULL;
     38  Header *fitsheader = NULL;
     39  Buffer *headbuffer = NULL;
    3640  if ((N = get_argument (argc, argv, "-fits"))) {
    3741    remove_argument (N, &argc, argv);
    3842    FITS = strcreate (argv[N]);
    3943    remove_argument (N, &argc, argv);
    40   }
    41 
     44
     45    if ((N = get_argument (argc, argv, "-header"))) {
     46      remove_argument (N, &argc, argv);
     47      if ((headbuffer = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     48      fitsheader = &headbuffer->header;
     49      remove_argument (N, &argc, argv);
     50    }
     51  }
     52 
    4253  /* option generate a FITS output table */
    4354  int CSV = FALSE;
    … …  
    124135
    125136  if (FITS) {
    126     int status = WriteVectorTableFITS (argv[1], FITS, vec, Nvec, append, compress, format, Ntile);
     137    int status = WriteVectorTableFITS (argv[1], FITS, fitsheader, vec, Nvec, append, compress, format, Ntile);
    127138    free (vec);
    128139    return status;
    … …  
    143154  if (ADD_HEADER) {
    144155    for (j = 0; j < Nvec; j++) {
     156      if (j == 0) fprintf (f, "# ");
    145157      if (CSV) {
    146158        fprintf (f, "%s,", vec[j][0].name);
    147159      } else {
    148         if (j == 0) fprintf (f, "# ");
    149160        fprintf (f, "%s ", vec[j][0].name);
    150161      }
    … …  
    155166  /* default output format */
    156167  if (format == (char *) NULL) {
     168    char padChar = CSV ? ',' : ' ';
    157169    for (i = 0; i < vec[0][0].Nelements; i++) {
    158170      for (j = 0; j < Nvec; j++) {
    159         if (vec[j][0].type == OPIHI_FLT) {
    160           if (CSV) {
    161             fprintf (f, "%.12g,", vec[j][0].elements.Flt[i]);
    162           } else {
    163             fprintf (f, "%.12g ", vec[j][0].elements.Flt[i]);
    164           }
    165         } else {
    166           if (CSV) {
    167             fprintf (f, OPIHI_INT_FMT",", vec[j][0].elements.Int[i]);
    168           } else {
    169             fprintf (f, OPIHI_INT_FMT" ", vec[j][0].elements.Int[i]);
    170           }
     171        switch (vec[j][0].type) {
     172          case OPIHI_FLT:
     173            fprintf (f, "%.12g%c", vec[j][0].elements.Flt[i], padChar);
     174            break;
     175          case OPIHI_INT:
     176            fprintf (f, OPIHI_INT_FMT"%c", vec[j][0].elements.Int[i], padChar);
     177            break;
     178          case OPIHI_STR:
     179            fprintf (f, "%s%c", vec[j][0].elements.Str[i], padChar);
     180            break;
    171181        }
    172182      }
    … …  
    215225        fmttype[j] = 'd';
    216226        break;
     227      case 's':
     228        fmttype[j] = 's';
     229        break;
    217230      default:
    218         gprint (GP_ERR, "syntax error in format (only e,f,d,c,x allowed)\n");
     231        gprint (GP_ERR, "syntax error in format (only e,f,d,c,x,s allowed)\n");
    219232        return (FALSE);
    220233    }
    … …  
    223236  strcat (fmtlist[Nvec-1], p0);
    224237 
     238  // check format types against vector types:
     239  for (j = 0; j < Nvec; j++) {
     240    switch (vec[j][0].type) {
     241      case OPIHI_FLT:
     242      case OPIHI_INT:
     243        if (fmttype[j] == 's') {
     244          gprint (GP_ERR, "mismatch between string format and numerical vector for %s\n", vec[j][0].name);
     245          return FALSE;
     246        }
     247        break;
     248      case OPIHI_STR:
     249        if (fmttype[j] != 's') {
     250          gprint (GP_ERR, "mismatch between numerical format and string vector for %s\n", vec[j][0].name);
     251          return FALSE;
     252        }
     253        break;
     254    }
     255  }
     256
    225257  for (i = 0; i < vec[0][0].Nelements; i++) {
    226258    for (j = 0; j < Nvec; j++) {
    … …  
    239271        }
    240272      }
     273      if (fmttype[j] == 's') {
     274        fprintf (f, fmtlist[j], vec[j][0].elements.Str[i]);
     275      }
    241276    }
    242277    fprintf (f, "\n");
    … …  
    260295    gprint (GP_ERR, "  -append : write to the end of the existing file\n");
    261296    gprint (GP_ERR, "  -fits NAME : write a fits table (extention name is NAME, column names match vector names)\n");
     297    gprint (GP_ERR, "     in FITS output context, -header takes an additional argument which is interpretted\n");
     298    gprint (GP_ERR, "     as a buffer containing header keywords to supplement the FITS table header\n");
    262299    gprint (GP_ERR, "  -csv : write a comma-separated values file (eg, to read in excel)\n");
    263300    gprint (GP_ERR, "  -f \"format\" : provide formatting codes for output:\n");
Note: See TracChangeset for help on using the changeset viewer.