IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 16, 2020, 1:54:47 PM (6 years ago)
Author:
tdeboer
Message:

revert to working Ohana build

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

Legend:

Unmodified
Added
Removed

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

    r41322 r41340  
    1919$(SRC)/init.$(ARCH).o          \
    2020$(SRC)/accum.$(ARCH).o          \
    21 $(SRC)/antialias.$(ARCH).o              \
    2221$(SRC)/applyfit1d.$(ARCH).o     \
    2322$(SRC)/applyfit2d.$(ARCH).o     \
    … …  
    7978$(SRC)/join.$(ARCH).o           \
    8079$(SRC)/jpeg.$(ARCH).o           \
    81 $(SRC)/kapamemory.$(ARCH).o             \
    8280$(SRC)/kern.$(ARCH).o           \
    8381$(SRC)/keyword.$(ARCH).o        \
    … …  
    9189$(SRC)/lookup.$(ARCH).o \
    9290$(SRC)/matrix.$(ARCH).o \
    93 $(SRC)/match1d.$(ARCH).o        \
    9491$(SRC)/match2d.$(ARCH).o        \
    9592$(SRC)/mkrgb.$(ARCH).o  \
    9693$(SRC)/mcreate.$(ARCH).o        \
    9794$(SRC)/medacc.$(ARCH).o \
    98 $(SRC)/mgaussdev.$(ARCH).o      \
    9995$(SRC)/mget.$(ARCH).o           \
    10096$(SRC)/mget3d.$(ARCH).o         \
    … …  
    10399$(SRC)/medimage.$(ARCH).o       \
    104100$(SRC)/medimage_commands.$(ARCH).o \
    105 $(SRC)/medimage_calc.$(ARCH).o \
    106101$(SRC)/mset.$(ARCH).o           \
    107102$(SRC)/needles.$(ARCH).o        \
    … …  
    117112$(SRC)/point.$(ARCH).o          \
    118113$(SRC)/ps.$(ARCH).o             \
    119 $(SRC)/pdf.$(ARCH).o            \
    120114$(SRC)/print_vectors.$(ARCH).o  \
    121115$(SRC)/mprint.$(ARCH).o         \
    … …  
    130124$(SRC)/queuesubstr.$(ARCH).o    \
    131125$(SRC)/queueinit.$(ARCH).o      \
    132 $(SRC)/queue2book.$(ARCH).o     \
    133126$(SRC)/radial.$(ARCH).o \
    134127$(SRC)/rd.$(ARCH).o             \
    … …  
    176169$(SRC)/vgauss.$(ARCH).o            \
    177170$(SRC)/vlorentz.$(ARCH).o          \
    178 $(SRC)/vsigmoid.$(ARCH).o          \
    179171$(SRC)/vellipse.$(ARCH).o          \
    180172$(SRC)/vmaxwell.$(ARCH).o          \
    … …  
    184176$(SRC)/vzload.$(ARCH).o            \
    185177$(SRC)/vpeaks.$(ARCH).o            \
    186 $(SRC)/vtransitions.$(ARCH).o     \
    187178$(SRC)/vpop.$(ARCH).o              \
    188179$(SRC)/vroll.$(ARCH).o             \
    … …  
    191182$(SRC)/vsmooth.$(ARCH).o           \
    192183$(SRC)/vstats.$(ARCH).o            \
    193 $(SRC)/virls.$(ARCH).o             \
    194 $(SRC)/vwtmean.$(ARCH).o \
    195184$(SRC)/xsection.$(ARCH).o          \
    196185$(SRC)/vsh.$(ARCH).o               \

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

    r41269 r41340  
    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 
    369  INT = FALSE;
    3710  if ((N = get_argument (argc, argv, "-int"))) {
    … …  
    4316    gprint (GP_ERR, "USAGE: create vector start end [delta] [-int]\n");
    4417    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");
    4618    return (FALSE);
    4719  }

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

    r41164 r41340  
    11# include "data.h"
    22
    3 enum {SUM, MEAN, MEDIAN, INNER, NGOOD};
    4 
    5 double 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 }
     3enum {SUM, MEAN, MEDIAN};
    654
    665int cut (int argc, char **argv) {
    676 
    687  int i, j, N, Nx, Ny, Mode;
    69   float *Vin, *Vbuf;
     8  float *Vin, *Vbuf, value;
    709  int sx, sy, nx, ny;
    7110  Vector *xvec, *yvec;
    … …  
    7312
    7413  Mode = SUM;
    75   if ((N = get_argument (argc, argv, "-sum"))) {
    76     remove_argument (N, &argc, argv);
    77     Mode = SUM;
    78   }
    7914  if ((N = get_argument (argc, argv, "-median"))) {
    8015    remove_argument (N, &argc, argv);
    … …  
    8419    remove_argument (N, &argc, argv);
    8520    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;
    9421  }
    9522
    … …  
    10936  Ny = buf[0].matrix.Naxis[1];
    11037
    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);
     38  if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) {
     39    gprint (GP_ERR, "region out of range\n");
    11440    return (FALSE);
    11541  }
    … …  
    13157
    13258    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;
     59      /* accumulate values */
     60      Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i;
     61      for (j = 0; j < ny; j++, Vin += Nx) {
     62        Vbuf[j] = *Vin;
    13763      }
    138       if (i + sx >= Nx) {
    139         yvec[0].elements.Flt[i] = NAN;
    140         continue;
     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; }
    14174      }
    142 
    143       /* 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 
    152       Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i;
    153       int Npix = 0;
    154       for (j = 0; j < ny; j++, Vin += Nx) {
    155         if (!isfinite(*Vin)) continue;
    156         Vbuf[Npix] = *Vin;
    157         Npix ++;
    158       }
    159       yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
     75      yvec[0].elements.Flt[i] = value;
    16076    }
    16177    free (Vbuf);
    … …  
    17389
    17490    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;
     91      /* accumulate values */
     92      Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx;
     93      for (j = 0; j < nx; j++, Vin ++) {
     94        Vbuf[j] = *Vin;
    17995      }
    180       if (i + sy >= Ny) {
    181         yvec[0].elements.Flt[i] = NAN;
    182         continue;
     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; }
    183106      }
    184 
    185       /* 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 
    194       Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx;
    195       int Npix = 0;
    196       for (j = 0; j < nx; j++, Vin ++) {
    197         if (!isfinite(*Vin)) continue;
    198         Vbuf[Npix] = *Vin;
    199         Npix ++;
    200       }
    201       yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
     107      yvec[0].elements.Flt[i] = value;
    202108    }
    203109    free (Vbuf);

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

    r41283 r41340  
    8080      case FIELD_TYPE_DOUBLE:
    8181      case FIELD_TYPE_TIME:
    82       case FIELD_TYPE_TIMESTAMP:
    8382      case FIELD_TYPE_DATE:
    8483      case FIELD_TYPE_DATETIME:
    … …  
    116115            break;
    117116          case FIELD_TYPE_TIME:
    118           case FIELD_TYPE_TIMESTAMP:
    119117          case FIELD_TYPE_DATE:
    120118          case FIELD_TYPE_DATETIME:

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

    r41164 r41340  
    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  */
    102
    113int extract (int argc, char **argv) {
    124 
    13   int N;
     5  int i, j;
     6  float *Vin, *Vout;
     7  int sx, sy, nx, ny, NX, NY;
     8  int Sx, Sy, Nx, Ny;
    149  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   }
    2410
    2511  if (argc != 11) {
    … …  
    2915
    3016  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    31   int NX = in[0].matrix.Naxis[0];
    32   int NY = in[0].matrix.Naxis[1];
     17  NX = in[0].matrix.Naxis[0];
     18  NY = in[0].matrix.Naxis[1];
    3319
    34   int sx = atof (argv[3]);
    35   int sy = atof (argv[4]);
    36   int nx = atof (argv[5]);
    37   int ny = atof (argv[6]);
     20  sx = atof (argv[3]);
     21  sy = atof (argv[4]);
     22  nx = atof (argv[5]);
     23  ny = atof (argv[6]);
    3824
    39   int Sx = atof (argv[7]);
    40   int Sy = atof (argv[8]);
    41   int Nx = atof (argv[9]);
    42   int Ny = atof (argv[10]);
     25  Sx = atof (argv[7]);
     26  Sy = atof (argv[8]);
     27  Nx = atof (argv[9]);
     28  Ny = atof (argv[10]);
    4329
    4430  if ((Sy + ny > Ny) || (Sx + nx > Nx)) {
    45     gprint (GP_ERR, "source pixels extend beyond target pixels\n");
     31    gprint (GP_ERR, "mismatch between source and dest regions\n");
    4632    gprint (GP_ERR, "%d + %d > %d or %d + %d > %d\n", Sy, ny, Ny, Sx, nx, Nx);
    47     // return (FALSE);
     33    return (FALSE);
    4834  }
    49 
    50   // XXX : allow source region to fall outside source image
    5135
    5236  /* region is not on first image */
    … …  
    5539      (sy > in[0].matrix.Naxis[1])) {
    5640    gprint (GP_ERR, "region outside of source image\n");
    57     // return (FALSE);
     41    return (FALSE);
    5842  }
    5943
     44  if ((Sx + nx > Nx) || (Sy + ny > Ny)) {
     45    gprint (GP_ERR, "source region larger than dest region\n");
     46    return (FALSE);
     47  }
    6048  if ((Sx < 0) || (Sy < 0)) {
    6149    gprint (GP_ERR, "dest region out of range\n");
    62     // return (FALSE);
     50    return (FALSE);
    6351  }
    6452
    … …  
    8876  }
    8977
    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++) {
     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++) {
    13284      if (i + sx < 0) continue;
    13385      if (i + sx >= NX) continue;
    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];
     86      *Vout = *Vin;
    14787    }
    14888  }
    … …  
    15191
    15292}
     93

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

    r41322 r41340  
    33
    44int accum            PROTO((int, char **));
    5 int antialias        PROTO((int, char **));
    65int applyfit         PROTO((int, char **));
    76int applyfit1d       PROTO((int, char **));
    … …  
    6867int join             PROTO((int, char **));
    6968int jpeg             PROTO((int, char **));
    70 int kapamemory       PROTO((int, char **));
    7169int kern             PROTO((int, char **));
    7270int keyword          PROTO((int, char **));
    … …  
    8179int lookup           PROTO((int, char **));
    8280int matrix           PROTO((int, char **));
    83 int match1d          PROTO((int, char **));
    8481int match2d          PROTO((int, char **));
    8582int mkrgb            PROTO((int, char **));
    8683int mcreate          PROTO((int, char **));
    8784int medacc           PROTO((int, char **));
    88 int mgaussdev        PROTO((int, char **));
    8985int mget             PROTO((int, char **));
    9086int mget3d           PROTO((int, char **));
    … …  
    10399int parity           PROTO((int, char **));
    104100int point            PROTO((int, char **));
    105 int pdf              PROTO((int, char **));
    106101int ps               PROTO((int, char **));
    107102int vprint           PROTO((int, char **));
    … …  
    116111int queuesubstr      PROTO((int, char **));
    117112int queuesize        PROTO((int, char **));
    118 int queue2book       PROTO((int, char **));
    119113int rd               PROTO((int, char **));
    120114int rdseg            PROTO((int, char **));
    … …  
    162156int vgauss           PROTO((int, char **));
    163157int vlorentz         PROTO((int, char **));
    164 int vsigmoid         PROTO((int, char **));
    165158int vellipse         PROTO((int, char **));
    166159int vmaxwell         PROTO((int, char **));
    … …  
    170163int vzload           PROTO((int, char **));
    171164int vstats           PROTO((int, char **));
    172 int vstats           PROTO((int, char **));
    173 int virls            PROTO((int, char **));
    174 int vwtmean          PROTO((int, char **));
    175165int vroll            PROTO((int, char **));
    176166int vshift           PROTO((int, char **));
    177167int vpeaks           PROTO((int, char **));
    178 int vtransitions     PROTO((int, char **));
    179168int vpop             PROTO((int, char **));
    180169int vsmooth          PROTO((int, char **));
    … …  
    198187static Command cmds[] = { 
    199188  {1, "accum",        accum,            "accumulate vector values in another vector"},
    200   {1, "antialias",    antialias,        "set anti-alias sigma value for display"},
    201189  {1, "applyfit",     applyfit1d,       "apply 1-d fit to new vector"},
    202190  {1, "applyfit1d",   applyfit1d,       "apply 1-d fit to new vector"},
    … …  
    267255  {1, "join",         join,             "find the join of two ID vectors"},
    268256  {1, "jpeg",         jpeg,             "convert display image to JPEG"},
    269   {1, "kapamemory",   kapamemory,       "manage kapa memory dump options"},
    270257  {1, "kern",         kern,             "convolve with 3x3 kernel"},
    271258  {1, "keyword",      keyword,          "extract a FITS keyword from image header"},
    … …  
    278265  {1, "imcreate",     mcreate,          "create an image"},
    279266  {1, "medacc",       medacc,           "accumulate vector values in another vector"},
    280   {1, "mgaussdev",    mgaussdev,        "generate a gaussian deviate image"},
    281267  {1, "mget",         mget,             "extract a vector from an image"},
    282268  {1, "mget3d",       mget3d,           "extract a vector from a 3D image"},
    … …  
    287273  {1, "medimage",     medimage_command, "median image manipulation"},
    288274  {1, "matrix",       matrix,           "matrix math operations"},
    289   {1, "match1d",      match1d,          "match 2 vectors and return matched indexes"},
    290275  {1, "match2d",      match2d,          "match 2 pairs of X,Y vectors and return matched indexes"},
    291276  {1, "mkrgb",        mkrgb,            "convert 3 images to rgb jpeg (use Kapa for better control)"},
    … …  
    304289  {1, "ppm",          jpeg,             "convert display graphic to PPM"},
    305290  {1, "ps",           ps,               "convert display to PostScript"},
    306   {1, "pdf",          pdf,              "convert display to PDF"},
    307291  {1, "print_vectors", vprint,          "print a set of vectors"},
    308292  {1, "vprint",       vprint,           "print a set of vectors"},
    … …  
    317301  {1, "queuesize",    queuesize,        "show queue size"},
    318302  {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"},
    321303  {1, "rd",           rd,               "load fits image"},
    322304  {1, "rdseg",        rdseg,            "read a segment of an image from a file"},
    … …  
    365347  {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
    366348  {1, "vlorentz",     vlorentz,         "fit a Lorentzian to a vector"},
    367   {1, "vsigmoid",     vsigmoid,         "fit a Sigmoid to a vector"},
    368349  {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
    369350  {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
    … …  
    373354  {1, "vload",        vload,            "load vectors as overlay on image display"},
    374355  {1, "vmaxwell",     vmaxwell,         "fit a Maxwellian to a 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"},
     356  {1, "vpeaks",       vpeaks,           "fine coord and flux of peaks in vector"},
    377357  {1, "vpop",         vpop,             "remove first element of a vector"},
    378358  {1, "vroll",        vroll,            "roll vector elements by 1 entry"},
    … …  
    380360  {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
    381361  {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"},
    384362  {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
    385363  {1, "vsh",          vsh,              "Vector Spherical Harmonics"},

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

    r41164 r41340  
    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   }
    169
    1710  /** check basic syntax **/
    … …  
    4336  y0 = yin[0].elements.Flt[0];
    4437 
    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   }
     38  /* in and out vectors are sorted */
     39  j = 0;
    4940
    5041  // every output pixel should get a value unless
    5142  // we are below the lowest in or above the highest in
    52   for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
     43  for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
    5344
    5445    yout[0].elements.Flt[i] = NAN;
    … …  
    7566
    7667  // fill in the rest with NANs
    77   int NinLast = xin[0].Nelements - 1;
    7868  while (i < yout[0].Nelements) {
    79     yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN;
     69    yout[0].elements.Flt[i] = NAN;
    8070    i++;
    8171  }
    8272
    8373  return (TRUE);
     74   
    8475}

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

    r41269 r41340  
    3232  }
    3333
    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;
     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    }
    5646  }
    5747  if (Variable) free (Variable);

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

    r41164 r41340  
    1515int medimage_add (int argc, char **argv) {
    1616
    17   int N;
    1817  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   }
    26 
    27   if (argc != 3) {
    28     gprint (GP_ERR, "USAGE: medimage add (name) (image) [-variance variance]\n");
     18
     19  if (argc != 3) {
     20    gprint (GP_ERR, "USAGE: medimage add (name) (image)\n");
    2921    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");
    3122    return FALSE;
    3223  }
    3324
    3425  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   }
    4326
    4427  MedImageType *median = FindMedImage (argv[1]);
    … …  
    6346  // new image should match existing medimage dimensions
    6447
    65   // AddMedImage (median, image, var);
    6648  int Ninput = median->Ninput;
    6749  median->Ninput ++;
    68   REALLOCATE (median->flx, float *, median->Ninput);
    69   REALLOCATE (median->var, float *, median->Ninput);
    70 
    71   ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
    72   memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    73 
    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);
     50  REALLOCATE (median->buffers, float *, median->Ninput);
     51
     52  ALLOCATE (median->buffers[Ninput], float, median->Nx*median->Ny);
     53  memcpy (median->buffers[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
     54
     55  return TRUE;
     56}
     57
     58int 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    }
    78120  }
    79121
    … …  
    82124
    83125/*
    84    int medimage_save (int argc, char **argv) {
    85 
    86    int N;
    87 
    88    int APPEND = FALSE;
    89    if ((N = get_argument (argc, argv, "-append"))) {
    90    APPEND = TRUE;
    91    remove_argument (N, &argc, argv);
    92    }
    93 
    94    if (argc != 3) {
    95    gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
    96    return FALSE;
    97    }
    98 
    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    }
    105 
    106    int medimage_load (int argc, char **argv) {
    107 
    108    if (argc != 3) {
    109    gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
    110    return FALSE;
    111    }
    112 
    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    }
     126int medimage_save (int argc, char **argv) {
     127
     128  int N;
     129
     130  int APPEND = FALSE;
     131  if ((N = get_argument (argc, argv, "-append"))) {
     132    APPEND = TRUE;
     133    remove_argument (N, &argc, argv);
     134  }
     135
     136  if (argc != 3) {
     137    gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
     138    return FALSE;
     139  }
     140
     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}
     147
     148int medimage_load (int argc, char **argv) {
     149
     150  if (argc != 3) {
     151    gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
     152    return FALSE;
     153  }
     154
     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}
    119161*/
    120162

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

    r41164 r41340  
    3939      if (*X < start) continue;
    4040      if (*X > end) continue;
    41       if (!isfinite(*Y)) continue;
    4241      if (*Y < ymax) continue;
    4342      xmax = *X;
    … …  
    6463      if (*X < start) continue;
    6564      if (*X > end) continue;
    66       if (!isfinite(*Y)) continue;
    6765      if (*Y < ymax) continue;
    6866      xmax = *X;

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

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

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

    r41305 r41340  
    3535
    3636// vector types
    37 enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_DATE, COLTYPE_CHAR, COLTYPE_STR, COLTYPE_HMS};
     37enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_CHAR};
    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, str, time, date, hms\n");
     93    gprint (GP_ERR, "       type is int, float, char, time\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");
    9895    gprint (GP_ERR, "       for time, values are human-readable date/time strings YYYY/MM/DD,hh:mm:ss\n");
    99     gprint (GP_ERR, "         date & time values are stored as a float using the TIMEFORMAT, TIMEREF values for the conversion\n");
     96    gprint (GP_ERR, "       the resulting vector is a float based on the TIMEFORMAT, TIMEREF values\n");
    10097
    10198    gprint (GP_ERR, "     -fits options:\n");
    … …  
    155152      if (!strcasecmp(ptr, "int"))   { coltype[i] = COLTYPE_INT; }
    156153      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; }
    159154      if (!strcasecmp(ptr, "time"))  { coltype[i] = COLTYPE_TIME; }
    160       if (!strcasecmp(ptr, "hms"))   { coltype[i] = COLTYPE_HMS; }
    161155      if (!coltype[i]) goto bad_colname;
    162156    }
    … …  
    216210        break;
    217211      case COLTYPE_FLT:
    218       case COLTYPE_HMS:
    219212      case COLTYPE_TIME:
    220       case COLTYPE_DATE:
     213        // note that COLTYPE_TIME is a type of float
    221214        ResetVector (vec[i], OPIHI_FLT, NELEM);
    222         break;
    223       case COLTYPE_STR:
    224         ResetVector (vec[i], OPIHI_STR, NELEM);
    225215        break;
    226216      case COLTYPE_CHAR:
    … …  
    297287        time_t tvalue;
    298288        int readStatus = FALSE;
    299         int dataStatus = FALSE;
    300289        // need to make the if cases for coltype[i]
    301290        switch (coltype[i]) {
    … …  
    322311              break;
    323312            }
    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             }
    340313          case COLTYPE_FLT:
    341314            readStatus = IsCSV ? dparse_csv (&dvalue, col[i], c0) : dparse (&dvalue, col[i], c0);
    … …  
    347320            vec[i][0].elements.Flt[Nelem] = readStatus ? dvalue : NAN;
    348321            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           }
    374322        }
    375323        if (!readStatus && VERBOSE) {
    … …  
    397345              break;
    398346            case COLTYPE_FLT:
    399             case COLTYPE_HMS:
    400347            case COLTYPE_TIME:
    401             case COLTYPE_DATE:
    402348              ResetVector (vec[i], OPIHI_FLT, NELEM);
    403               break;
    404             case COLTYPE_STR:
    405               ResetVector (vec[i], OPIHI_STR, NELEM);
    406349              break;
    407350            case COLTYPE_CHAR:
    … …  
    423366        break;
    424367      case COLTYPE_FLT:
    425       case COLTYPE_HMS:
    426368      case COLTYPE_TIME:
    427       case COLTYPE_DATE:
    428369        ResetVector (vec[i], OPIHI_FLT, Nelem);
    429         break;
    430       case COLTYPE_STR:
    431         ResetVector (vec[i], OPIHI_STR, Nelem);
    432370        break;
    433371      case COLTYPE_CHAR:
    … …  
    531469  }
    532470
    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 
     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
    539473  // if (Ny > 10000), force CharAsVectors
    540474  int CharAsVectors = FALSE;
    … …  
    542476    remove_argument (N, &argc, argv);
    543477    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;
    549478  }
    550479
    … …  
    722651    if (!FITS_TRANSPOSE) {
    723652      // read string column into a list rather than a vector
    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);
     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, '$');
    741665          }
    742           sprintf (varname, "%s:n", fieldName);
    743           set_int_variable (varname, Ny);
    744           continue;
     666          set_str_variable (varname, value);
     667          free (value);
    745668        }
    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         }
     669        sprintf (varname, "%s:n", fieldName);
     670        set_int_variable (varname, Ny);
     671        continue;
    762672      }
    763673
    764674      // define the multifield vector names (Nval vectors x Ny elements)
    765       // CharAsVectors is handled below automatically
    766675      ALLOCATE (vec, Vector *, Nval);
    767676      for (j = 0; j < Nval; j++) {

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

    r41164 r41340  
    33int roll (int argc, char **argv) {
    44 
     5  int Nbytes, Nextra;
     6  int dX, dx, dy, nx, ny;
    57  Buffer *buf;
    6 
    7   // this function is probably not needed
    8   gprint (GP_ERR, "ERROR: use 'shift' instead of 'roll'\n");
    9   return (FALSE);
    108
    119  if (argc != 4) {
    … …  
    1614  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    1715
    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;
     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");
    4320  }
    4421
    45   if (dx == 0) {
    46     // moves in just dy can be done row-by-row
     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];
    4727
    48     ALLOCATE_PTR (output, float, Nx*Ny);
     28  dX = abs(dx);
     29  Nbytes = nx * ny * sizeof (float);
     30  Nextra = (nx * ny + dX) * sizeof (float);
    4931
    50     float *input = (float *) buf[0].matrix.buffer;
     32  REALLOCATE (buf[0].matrix.buffer, char, Nextra);
    5133
    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;
     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);
    6738  }
    6839
    69   ALLOCATE_PTR (output, float, Nx*Ny);
    70   float *input = (float *) buf[0].matrix.buffer;
     40  REALLOCATE (buf[0].matrix.buffer, char, Nbytes);
    7141
    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;
     42  return (TRUE);
     43
    8444}
    8545
    86 

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

    r41269 r41340  
    33int rotate (int argc, char **argv) {
    44 
    5   int i, j, NX, NY;
     5  int i, j, NX, NY, X, Y, Lx, Ly, N;
    66  float *in_buff, *out_buff, *c;
    7   double Xo, Yo, x, y, X1, Y1;
     7  double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, 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   double angle = atof (argv[2]);
    224   double CosAngle = cos (angle*RAD_DEG);
    225   double SinAngle = sin (angle*RAD_DEG);
     223  angle = atof (argv[2]);
     224  CosAngle = cos (angle*RAD_DEG);
     225  SinAngle = sin (angle*RAD_DEG);
    226226 
    227227  gprint (GP_ERR, "rotating: %f %f %f\n", angle, CosAngle, SinAngle);
    228228
    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);
     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  */
    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   X1 = (Xo - NX/2)*CosAngle + (Yo - NY/2)*SinAngle + Lx/2;
    244   Y1 = (NX/2 - Xo)*SinAngle + (Yo - NY/2)*CosAngle + Ly/2;
     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;
    245249  gfits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    246250  gfits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    … …  
    251255  gfits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    252256  gfits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    253   PC11 = pc11*CosAngle + pc21*SinAngle;
    254   PC12 = pc12*CosAngle + pc22*SinAngle;
    255   PC21 = pc21*CosAngle - pc11*SinAngle;
    256   PC22 = pc22*CosAngle - pc12*SinAngle;
     257  PC11 = pc11*CosAngle - pc21*SinAngle;
     258  PC12 = pc12*CosAngle - pc22*SinAngle;
     259  PC21 = pc21*CosAngle + pc11*SinAngle;
     260  PC22 = pc22*CosAngle + pc12*SinAngle;
    257261  gfits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    258262  gfits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    … …  
    266270  gfits_create_matrix (&buf[0].header, &buf[0].matrix);
    267271  gfits_print_alt (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
    268 
    269272  out_buff = (float *)buf[0].matrix.buffer;
    270 
    271273  for (j = 0; j < Ly; j++) {
    272274    for (i = 0; i < Lx; i++, out_buff++) {
    273275
    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 
     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;
    295289      *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];
    297290    }
    298291  }

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

    r41164 r41340  
    77int spline_load (int argc, char **argv);
    88int spline_save (int argc, char **argv);
    9 int spline_print (int argc, char **argv);
    109int spline_delete (int argc, char **argv);
    1110int spline_rename (int argc, char **argv);
    … …  
    1918  {1, "load",       spline_load,       "write a spline to a FITS file"},
    2019  {1, "save",       spline_save,       "read a spline from a FITS file"},
    21   {1, "print",      spline_print,      "print a spline to stdout"},
    2220  {1, "delete",     spline_delete,     "delete a spline"},
    2321  {1, "rename",     spline_rename,     "rename a spline"},
    … …  
    3735  gprint (GP_ERR, "    spline load   (spline) (filename)            : load a spline from a FITS file\n");
    3836  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");
    4037
    4138  return FALSE;

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

    r41164 r41340  
    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");
    3530    return FALSE;
    3631  }
    … …  
    9893}
    9994
    100 int 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 
    12295int spline_load (int argc, char **argv) {
    12396
    … …  
    136109int spline_delete (int argc, char **argv) {
    137110
    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   }
     111  int status;
     112  Spline *spline;
    145113
    146114  if (argc != 2) {
    … …  
    149117  }
    150118
    151   Spline *spline = FindSpline (argv[1]);
     119  spline = FindSpline (argv[1]);
    152120  if (spline == NULL) {
    153     if (QUIET) return TRUE;
    154121    gprint (GP_ERR, "spline %s not found\n", argv[1]);
    155122    return FALSE;
    … …  
    163130int spline_rename (int argc, char **argv) {
    164131
     132  Spline *spline;
     133
    165134  if (argc != 3) {
    166135    gprint (GP_ERR, "USAGE: spline rename (spline) (newname)\n");
    … …  
    168137  }
    169138
    170   Spline *spline = FindSpline (argv[1]);
     139  spline = FindSpline (argv[1]);
    171140  if (spline == NULL) {
    172141    gprint (GP_ERR, "spline %s not found\n", argv[1]);
    173142    return FALSE;
    174   }
    175 
    176   Spline *newname = FindSpline (argv[2]);
    177   if (newname != NULL) {
    178     DeleteSpline (newname);
    179143  }
    180144
    … …  
    183147  return TRUE;
    184148}
     149
     150/*
     151int 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

    r41164 r41340  
    44 memtest1
    55end
    6 
    7 macro 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  
    16 end
    17 
    18 macro 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
    46 end
    47 
    486
    497# Test if cut works
    … …  
    7937 end
    8038end
    81 
    82 macro 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    
    120 end

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

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

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

    r41164 r41340  
    11
    2 macro 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
     2macro go
     3 mcreate a 30 30
     4 for i 0 40
     5  set a$i = zero(a) + $i
    176 end
    187
    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
     8 for i 0 40
     9  medimage add t1 a$i
     10 end
    3011end
    31 
    32 macro 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
    61 end
    62  
    63 macro 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
    112 end
    113 
    114 macro 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
    165 end
    166 
    167 
    168 ###################33
    169 
    170 
    171 macro 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)}"
    205 end
    206 
    207 macro 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)}"
    243 end
    244  
    245 macro 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)}
    296 end
    297 
    298 macro 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)})"
    357 end
    358 
    359 macro 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)})"
    425 end
    426 
    427 ##############################
    428 macro 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)})"
    465 end
    466 
    467 macro 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
    540 end

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

    r41164 r41340  
    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]);
    3727    remove_argument (N, &argc, argv);
    3828  }
    … …  
    7464  }
    7565
    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 
    8866  // use bisection to find the value
    8967
    … …  
    9371
    9472  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     }
    10473    while (Nhi - Nlo > 10) {
    10574      N = 0.5*(Nlo + Nhi);
    … …  
    12089    }
    12190  } 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     }
    13191    while (Nhi - Nlo > 10) {
    13292      N = 0.5*(Nlo + Nhi);

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

    r41269 r41340  
    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: if colormap name is given as file:path/to/file, the colormap is read from the file\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");
    4045    gprint (GP_ERR, "   the colormap file must contain 4 columns: f R G B; each line defines a color transition.\n");
    4146    gprint (GP_ERR, "   f: 0 - 1 defines the scale value for the transition\n");
    4247    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");
    4356    return (FALSE);
    4457  }

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

    r41164 r41340  
    33/* local private functions */
    44opihi_flt fgaussOD (opihi_flt, opihi_flt *, int, opihi_flt *);
    5 int vgauss_apply (int argc, char **argv);
    65
    76# define GET_VAR(V,A) \
    … …  
    2120  Vector *xvec, *yvec, *svec, *ovec;
    2221  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   }
    2922
    3023  Quiet = FALSE;
    … …  
    111104  } 
    112105  if (!Quiet) gprint (GP_ERR, "%d iterations\n", i);
    113   set_variable ("Chisq", chisq);
    114   set_variable ("ChisqNu", chisq / (float) (Npts - Npar));
    115106
    116107  for (i = 0; i < Npts; i++) {
    … …  
    134125}
    135126
    136 int 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 
    170127/* pars: x_o, sigma, I, back */
    171128opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
    … …  
    178135  f = par[2]*r + par[3];
    179136
    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   }
     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;
    186141 
    187142  return (f);

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

    r41269 r41340  
    3232  }
    3333
    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
     34  /* option generate a FITS output table */
    3735  char *FITS = NULL;
    38   Header *fitsheader = NULL;
    39   Buffer *headbuffer = NULL;
    4036  if ((N = get_argument (argc, argv, "-fits"))) {
    4137    remove_argument (N, &argc, argv);
    4238    FITS = strcreate (argv[N]);
    4339    remove_argument (N, &argc, argv);
    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  
     40  }
     41
    5342  /* option generate a FITS output table */
    5443  int CSV = FALSE;
    … …  
    135124
    136125  if (FITS) {
    137     int status = WriteVectorTableFITS (argv[1], FITS, fitsheader, vec, Nvec, append, compress, format, Ntile);
     126    int status = WriteVectorTableFITS (argv[1], FITS, vec, Nvec, append, compress, format, Ntile);
    138127    free (vec);
    139128    return status;
    … …  
    154143  if (ADD_HEADER) {
    155144    for (j = 0; j < Nvec; j++) {
    156       if (j == 0) fprintf (f, "# ");
    157145      if (CSV) {
    158146        fprintf (f, "%s,", vec[j][0].name);
    159147      } else {
     148        if (j == 0) fprintf (f, "# ");
    160149        fprintf (f, "%s ", vec[j][0].name);
    161150      }
    … …  
    166155  /* default output format */
    167156  if (format == (char *) NULL) {
    168     char padChar = CSV ? ',' : ' ';
    169157    for (i = 0; i < vec[0][0].Nelements; i++) {
    170158      for (j = 0; j < Nvec; j++) {
    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;
     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          }
    181171        }
    182172      }
    … …  
    225215        fmttype[j] = 'd';
    226216        break;
    227       case 's':
    228         fmttype[j] = 's';
    229         break;
    230217      default:
    231         gprint (GP_ERR, "syntax error in format (only e,f,d,c,x,s allowed)\n");
     218        gprint (GP_ERR, "syntax error in format (only e,f,d,c,x allowed)\n");
    232219        return (FALSE);
    233220    }
    … …  
    236223  strcat (fmtlist[Nvec-1], p0);
    237224 
    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 
    257225  for (i = 0; i < vec[0][0].Nelements; i++) {
    258226    for (j = 0; j < Nvec; j++) {
    … …  
    271239        }
    272240      }
    273       if (fmttype[j] == 's') {
    274         fprintf (f, fmtlist[j], vec[j][0].elements.Str[i]);
    275       }
    276241    }
    277242    fprintf (f, "\n");
    … …  
    295260    gprint (GP_ERR, "  -append : write to the end of the existing file\n");
    296261    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");
    299262    gprint (GP_ERR, "  -csv : write a comma-separated values file (eg, to read in excel)\n");
    300263    gprint (GP_ERR, "  -f \"format\" : provide formatting codes for output:\n");
Note: See TracChangeset for help on using the changeset viewer.