IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 27, 2019, 11:22:02 AM (7 years ago)
Author:
eugene
Message:

extend cut function to apply alternate stats to the output vector; extend extract to allow source region to fall beyond bounds of input image; add commands antialias, kapamemory, match1d, mgaussdev, vsigmoid (fit a sigmoid), virls (irls on vector), vwtmean (weighted mean), vtransitions; extend medimage to report variance; fix error in rotate; update some tests; add vgauss -apply option

Location:
trunk/Ohana/src/opihi/cmd.data
Files:
17 edited
14 copied

Legend:

Unmodified
Added
Removed

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

    r40637 r41164  
    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        \
    … …  
    169174$(SRC)/vgauss.$(ARCH).o            \
    170175$(SRC)/vlorentz.$(ARCH).o          \
     176$(SRC)/vsigmoid.$(ARCH).o          \
    171177$(SRC)/vellipse.$(ARCH).o          \
    172178$(SRC)/vmaxwell.$(ARCH).o          \
    … …  
    176182$(SRC)/vzload.$(ARCH).o            \
    177183$(SRC)/vpeaks.$(ARCH).o            \
     184$(SRC)/vtransitions.$(ARCH).o     \
    178185$(SRC)/vpop.$(ARCH).o              \
    179186$(SRC)/vroll.$(ARCH).o             \
    … …  
    182189$(SRC)/vsmooth.$(ARCH).o           \
    183190$(SRC)/vstats.$(ARCH).o            \
     191$(SRC)/virls.$(ARCH).o             \
     192$(SRC)/vwtmean.$(ARCH).o \
    184193$(SRC)/xsection.$(ARCH).o          \
    185194$(SRC)/vsh.$(ARCH).o               \

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

    r33662 r41164  
    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/extract.c

    r28241 r41164  
    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

    r40637 r41164  
    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 **));
    … …  
    156160int vgauss           PROTO((int, char **));
    157161int vlorentz         PROTO((int, char **));
     162int vsigmoid         PROTO((int, char **));
    158163int vellipse         PROTO((int, char **));
    159164int vmaxwell         PROTO((int, char **));
    … …  
    163168int vzload           PROTO((int, char **));
    164169int vstats           PROTO((int, char **));
     170int vstats           PROTO((int, char **));
     171int virls            PROTO((int, char **));
     172int vwtmean          PROTO((int, char **));
    165173int vroll            PROTO((int, char **));
    166174int vshift           PROTO((int, char **));
    167175int vpeaks           PROTO((int, char **));
     176int vtransitions     PROTO((int, char **));
    168177int vpop             PROTO((int, char **));
    169178int vsmooth          PROTO((int, char **));
    … …  
    187196static Command cmds[] = { 
    188197  {1, "accum",        accum,            "accumulate vector values in another vector"},
     198  {1, "antialias",    antialias,        "set anti-alias sigma value for display"},
    189199  {1, "applyfit",     applyfit1d,       "apply 1-d fit to new vector"},
    190200  {1, "applyfit1d",   applyfit1d,       "apply 1-d fit to new vector"},
    … …  
    255265  {1, "join",         join,             "find the join of two ID vectors"},
    256266  {1, "jpeg",         jpeg,             "convert display image to JPEG"},
     267  {1, "kapamemory",   kapamemory,       "manage kapa memory dump options"},
    257268  {1, "kern",         kern,             "convolve with 3x3 kernel"},
    258269  {1, "keyword",      keyword,          "extract a FITS keyword from image header"},
    … …  
    265276  {1, "imcreate",     mcreate,          "create an image"},
    266277  {1, "medacc",       medacc,           "accumulate vector values in another vector"},
     278  {1, "mgaussdev",    mgaussdev,        "generate a gaussian deviate image"},
    267279  {1, "mget",         mget,             "extract a vector from an image"},
    268280  {1, "mget3d",       mget3d,           "extract a vector from a 3D image"},
    … …  
    273285  {1, "medimage",     medimage_command, "median image manipulation"},
    274286  {1, "matrix",       matrix,           "matrix math operations"},
     287  {1, "match1d",      match1d,          "match 2 vectors and return matched indexes"},
    275288  {1, "match2d",      match2d,          "match 2 pairs of X,Y vectors and return matched indexes"},
    276289  {1, "mkrgb",        mkrgb,            "convert 3 images to rgb jpeg (use Kapa for better control)"},
    … …  
    347360  {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
    348361  {1, "vlorentz",     vlorentz,         "fit a Lorentzian to a vector"},
     362  {1, "vsigmoid",     vsigmoid,         "fit a Sigmoid to a vector"},
    349363  {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
    350364  {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
    … …  
    354368  {1, "vload",        vload,            "load vectors as overlay on image display"},
    355369  {1, "vmaxwell",     vmaxwell,         "fit a Maxwellian to a vector"},
    356   {1, "vpeaks",       vpeaks,           "fine coord and flux of peaks in vector"},
     370  {1, "vpeaks",       vpeaks,           "find coord and flux of peaks in vector"},
     371  {1, "vtransitions", vtransitions,     "find points in vector that cross the transition value"},
    357372  {1, "vpop",         vpop,             "remove first element of a vector"},
    358373  {1, "vroll",        vroll,            "roll vector elements by 1 entry"},
    … …  
    360375  {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
    361376  {1, "vstats",       vstats,           "statistics on a vector"},
     377  {1, "vwtmean",      vwtmean,          "weighted mean of a vector"},
     378  {1, "virls",        virls,            "IRLS mean of a vector"},
    362379  {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
    363380  {1, "vsh",          vsh,              "Vector Spherical Harmonics"},

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

    r40007 r41164  
    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/medimage_commands.c

    r40654 r41164  
    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

    r20936 r41164  
    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/roll.c

    r7917 r41164  
    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

    r34088 r41164  
    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;
    273270  for (j = 0; j < Ly; j++) {
    274271    for (i = 0; i < Lx; i++, out_buff++) {
    275272
    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;
     273      float xo = (i - Lx/2);
     274      float yo = (j - Ly/2);
     275
     276      x =  xo*CosAngle + yo*SinAngle;
     277      y = -xo*SinAngle + yo*CosAngle;
     278
     279      int I = x + NX/2;
     280      int J = y + NY/2;
     281
     282      if (I < 0) continue;
     283      if (I >= NX - 1) continue;
     284      if (J < 0) continue;
     285      if (J >= NY - 1) continue;
     286
     287      c = &in_buff[I + NX*J];
     288
     289      int X = (int) x;
     290      int Y = (int) y;
     291      double fx = x - X;
     292      double fy = y - Y;
     293
    289294      *out_buff = (c[0]*(1-fx) + c[1]*fx)*(1-fy) + (c[NX+1]*fx + c[NX]*(1-fx))*fy;
     295//    *out_buff = c[0];
    290296    }
    291297  }

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

    r39457 r41164  
    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

    r39457 r41164  
    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

    r16056 r41164  
    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

    r14176 r41164  
    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

    r36679 r41164  
    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

    r40623 r41164  
    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/vgauss.c

    r39457 r41164  
    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

    r40574 r41164  
    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;
    … …  
    260271    gprint (GP_ERR, "  -append : write to the end of the existing file\n");
    261272    gprint (GP_ERR, "  -fits NAME : write a fits table (extention name is NAME, column names match vector names)\n");
     273    gprint (GP_ERR, "     in FITS output context, -header takes an additional argument which is interpretted\n");
     274    gprint (GP_ERR, "     as a buffer containing header keywords to supplement the FITS table header\n");
    262275    gprint (GP_ERR, "  -csv : write a comma-separated values file (eg, to read in excel)\n");
    263276    gprint (GP_ERR, "  -f \"format\" : provide formatting codes for output:\n");
Note: See TracChangeset for help on using the changeset viewer.