Index: trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.data/Makefile (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/Makefile (revision 41340) @@ -19,5 +19,4 @@ $(SRC)/init.$(ARCH).o \ $(SRC)/accum.$(ARCH).o \ -$(SRC)/antialias.$(ARCH).o \ $(SRC)/applyfit1d.$(ARCH).o \ $(SRC)/applyfit2d.$(ARCH).o \ @@ -79,5 +78,4 @@ $(SRC)/join.$(ARCH).o \ $(SRC)/jpeg.$(ARCH).o \ -$(SRC)/kapamemory.$(ARCH).o \ $(SRC)/kern.$(ARCH).o \ $(SRC)/keyword.$(ARCH).o \ @@ -91,10 +89,8 @@ $(SRC)/lookup.$(ARCH).o \ $(SRC)/matrix.$(ARCH).o \ -$(SRC)/match1d.$(ARCH).o \ $(SRC)/match2d.$(ARCH).o \ $(SRC)/mkrgb.$(ARCH).o \ $(SRC)/mcreate.$(ARCH).o \ $(SRC)/medacc.$(ARCH).o \ -$(SRC)/mgaussdev.$(ARCH).o \ $(SRC)/mget.$(ARCH).o \ $(SRC)/mget3d.$(ARCH).o \ @@ -103,5 +99,4 @@ $(SRC)/medimage.$(ARCH).o \ $(SRC)/medimage_commands.$(ARCH).o \ -$(SRC)/medimage_calc.$(ARCH).o \ $(SRC)/mset.$(ARCH).o \ $(SRC)/needles.$(ARCH).o \ @@ -117,5 +112,4 @@ $(SRC)/point.$(ARCH).o \ $(SRC)/ps.$(ARCH).o \ -$(SRC)/pdf.$(ARCH).o \ $(SRC)/print_vectors.$(ARCH).o \ $(SRC)/mprint.$(ARCH).o \ @@ -130,5 +124,4 @@ $(SRC)/queuesubstr.$(ARCH).o \ $(SRC)/queueinit.$(ARCH).o \ -$(SRC)/queue2book.$(ARCH).o \ $(SRC)/radial.$(ARCH).o \ $(SRC)/rd.$(ARCH).o \ @@ -176,5 +169,4 @@ $(SRC)/vgauss.$(ARCH).o \ $(SRC)/vlorentz.$(ARCH).o \ -$(SRC)/vsigmoid.$(ARCH).o \ $(SRC)/vellipse.$(ARCH).o \ $(SRC)/vmaxwell.$(ARCH).o \ @@ -184,5 +176,4 @@ $(SRC)/vzload.$(ARCH).o \ $(SRC)/vpeaks.$(ARCH).o \ -$(SRC)/vtransitions.$(ARCH).o \ $(SRC)/vpop.$(ARCH).o \ $(SRC)/vroll.$(ARCH).o \ @@ -191,6 +182,4 @@ $(SRC)/vsmooth.$(ARCH).o \ $(SRC)/vstats.$(ARCH).o \ -$(SRC)/virls.$(ARCH).o \ -$(SRC)/vwtmean.$(ARCH).o \ $(SRC)/xsection.$(ARCH).o \ $(SRC)/vsh.$(ARCH).o \ Index: trunk/Ohana/src/opihi/cmd.data/antialias.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/antialias.c (revision 41323) +++ (revision ) @@ -1,26 +1,0 @@ -# include "data.h" - -int antialias (int argc, char **argv) { - - int N, kapa; - Graphdata graphmode; - - char *name = NULL; - if ((N = get_argument (argc, argv, "-n"))) { - remove_argument (N, &argc, argv); - name = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if (!GetGraph (&graphmode, &kapa, name)) return (FALSE); - FREE (name); - - if (argc != 2) { - gprint (GP_ERR, "USAGE: antialias (sigma)\n"); - return (FALSE); - } - - float sigma = atof (argv[1]); - - KapaSetSmoothSigma (kapa, sigma); - return (TRUE); -} Index: trunk/Ohana/src/opihi/cmd.data/create.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/create.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/create.c (revision 41340) @@ -7,31 +7,4 @@ Vector *vec; - // create a vector of empty strings - if ((N = get_argument (argc, argv, "-str"))) { - remove_argument (N, &argc, argv); - - char *stringValue = NULL; - if ((N = get_argument (argc, argv, "-value"))) { - remove_argument (N, &argc, argv); - stringValue = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - - if (argc != 3) { - gprint (GP_ERR, "USAGE: create vector Nelements -value word\n"); - return (FALSE); - } - - if ((vec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) return (FALSE); - int Nelements = atoi (argv[2]); - - // a newly reset vector has NULL-valued pointers - ResetVector (vec, OPIHI_STR, Nelements); - for (int i = 0; i < Nelements; i++) { - vec[0].elements.Str[i] = stringValue ? strcreate (stringValue) : strcreate (""); - } - return TRUE; - } - INT = FALSE; if ((N = get_argument (argc, argv, "-int"))) { @@ -43,5 +16,4 @@ gprint (GP_ERR, "USAGE: create vector start end [delta] [-int]\n"); gprint (GP_ERR, " -int : resulting vector is integer type (delta must be integer)\n"); - gprint (GP_ERR, " -str : resulting vector is string type (only give number of elements)\n"); return (FALSE); } Index: trunk/Ohana/src/opihi/cmd.data/cut.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/cut.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/cut.c (revision 41340) @@ -1,71 +1,10 @@ # include "data.h" -enum {SUM, MEAN, MEDIAN, INNER, NGOOD}; - -double cutstat (float *value, int Nvalue, int mode) { - - double output = 0.0; - - if (Nvalue == 0) return NAN; - - if (mode == MEDIAN) { - fsort (value, Nvalue); - if (Nvalue % 2) { - int Ncenter = Nvalue / 2; - myAssert ((Ncenter >= 0) && (Ncenter < Nvalue), "cutstat"); - output = value[Ncenter]; - } else { - int Ncenter = Nvalue / 2 - 1; - myAssert ((Ncenter >= 0) && (Ncenter < Nvalue - 1), "cutstat"); - output = 0.5*(value[Ncenter] + value[Ncenter + 1]); - } - return output; - } - if (mode == INNER) { - fsort (value, Nvalue); - - int Ns = 0, Ne = 0; - if (Nvalue % 2) { - // for an odd number of points take the same number below - // and above the center value - int Ncenter = Nvalue / 2; - int Nquarter = 0.25*Nvalue; - Ns = Ncenter - Nquarter; - Ne = Ncenter + Nquarter; - } else { - // for an even number, the middle lies between two points - // take the same number below as above - int Ncenter = (int)(Nvalue / 2) - 1; - int Nquarter = 0.25*Nvalue; - Ns = Ncenter - Nquarter; - Ne = Ncenter + 1 + Nquarter; - } - int Nv = 0; - for (int i = Ns; i <= Ne; i++) { - myAssert ((i >= 0) && (i < Nvalue), "cutstat"); - output += value[i]; - Nv ++; - } - output /= Nv; - return output; - } - if (mode == NGOOD) { - return Nvalue; - } - if ((mode == MEAN) || (mode == SUM)) { - for (int i = 0; i < Nvalue; i++) { - myAssert ((i >= 0) && (i < Nvalue), "cutstat"); - output += value[i]; - } - if (mode == MEAN) { output /= Nvalue; } - return output; - } - return NAN; -} +enum {SUM, MEAN, MEDIAN}; int cut (int argc, char **argv) { int i, j, N, Nx, Ny, Mode; - float *Vin, *Vbuf; + float *Vin, *Vbuf, value; int sx, sy, nx, ny; Vector *xvec, *yvec; @@ -73,8 +12,4 @@ Mode = SUM; - if ((N = get_argument (argc, argv, "-sum"))) { - remove_argument (N, &argc, argv); - Mode = SUM; - } if ((N = get_argument (argc, argv, "-median"))) { remove_argument (N, &argc, argv); @@ -84,12 +19,4 @@ remove_argument (N, &argc, argv); Mode = MEAN; - } - if ((N = get_argument (argc, argv, "-inner"))) { - remove_argument (N, &argc, argv); - Mode = INNER; - } - if ((N = get_argument (argc, argv, "-ngood"))) { - remove_argument (N, &argc, argv); - Mode = NGOOD; } @@ -109,7 +36,6 @@ Ny = buf[0].matrix.Naxis[1]; - // ny & nx do not need to be constrained by the buffer, but they do need to be sensible - if ((nx < 0) || (ny < 0) || (nx > 1e8) || (ny > 1e8)) { - gprint (GP_ERR, "warning : extraction size is crazy: %d,%d\n", nx, ny); + if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) { + gprint (GP_ERR, "region out of range\n"); return (FALSE); } @@ -131,31 +57,21 @@ for (i = 0; i < nx; i++) { - // for out-of-range areas, set the output pixels to NAN - if (i + sx < 0) { - yvec[0].elements.Flt[i] = NAN; - continue; + /* accumulate values */ + Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i; + for (j = 0; j < ny; j++, Vin += Nx) { + Vbuf[j] = *Vin; } - if (i + sx >= Nx) { - yvec[0].elements.Flt[i] = NAN; - continue; + /* apply stat of choice */ + if (Mode == MEDIAN) { + fsort (Vbuf, ny); + value = Vbuf[(int)(0.5*ny)]; + } else { + value = 0; + for (j = 0; j < ny; j++) { + value += Vbuf[j]; + } + if (Mode == MEAN) { value /= ny; } } - - /* accumulate values */ - // skip out-of-range areas in y - if (sy < 0) { - ny += sy; - sy = 0; - } - if (sy > Ny) sy = Ny; - if (sy + ny >= Ny) ny = Ny - sy; - - Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i; - int Npix = 0; - for (j = 0; j < ny; j++, Vin += Nx) { - if (!isfinite(*Vin)) continue; - Vbuf[Npix] = *Vin; - Npix ++; - } - yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode); + yvec[0].elements.Flt[i] = value; } free (Vbuf); @@ -173,31 +89,21 @@ for (i = 0; i < ny; i++) { - // for out-of-range areas, set the output pixels to NAN - if (i + sy < 0) { - yvec[0].elements.Flt[i] = NAN; - continue; + /* accumulate values */ + Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx; + for (j = 0; j < nx; j++, Vin ++) { + Vbuf[j] = *Vin; } - if (i + sy >= Ny) { - yvec[0].elements.Flt[i] = NAN; - continue; + /* apply stat of choice */ + if (Mode == MEDIAN) { + fsort (Vbuf, nx); + value = Vbuf[(int)(0.5*nx)]; + } else { + value = 0; + for (j = 0; j < nx; j++) { + value += Vbuf[j]; + } + if (Mode == MEAN) { value /= nx; } } - - /* accumulate values */ - // skip out-of-range areas in x - if (sx < 0) { - nx += sx; - sx = 0; - } - if (sx > Nx) sx = Nx; - if (sx + nx >= Nx) nx = Nx - sx; - - Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx; - int Npix = 0; - for (j = 0; j < nx; j++, Vin ++) { - if (!isfinite(*Vin)) continue; - Vbuf[Npix] = *Vin; - Npix ++; - } - yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode); + yvec[0].elements.Flt[i] = value; } free (Vbuf); Index: trunk/Ohana/src/opihi/cmd.data/dbselect.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/dbselect.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/dbselect.c (revision 41340) @@ -80,5 +80,4 @@ case FIELD_TYPE_DOUBLE: case FIELD_TYPE_TIME: - case FIELD_TYPE_TIMESTAMP: case FIELD_TYPE_DATE: case FIELD_TYPE_DATETIME: @@ -116,5 +115,4 @@ break; case FIELD_TYPE_TIME: - case FIELD_TYPE_TIMESTAMP: case FIELD_TYPE_DATE: case FIELD_TYPE_DATETIME: Index: trunk/Ohana/src/opihi/cmd.data/extract.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/extract.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/extract.c (revision 41340) @@ -1,25 +1,11 @@ # include "data.h" - -/* : source image, must exist - : target image -- if it does not exist, it is created of size (Nx,Ny) - sx, sy : source starting coordinate -- need not be on a valid image pixel - nx, ny : number of source pixels -- currently must not go out of bounds -> allow out-of-bounds - Sx, Sy : target starting coordinate -- need not be on a valid image pixel - Nx, Ny : redundant information UNLESS does exist (in which case it is required information) -> make Nx,Ny optional if exists? - */ int extract (int argc, char **argv) { - int N; + int i, j; + float *Vin, *Vout; + int sx, sy, nx, ny, NX, NY; + int Sx, Sy, Nx, Ny; Buffer *in, *out; - - float initValue = 0.0; - int initOutput = FALSE; - if ((N = get_argument (argc, argv, "-init"))) { - remove_argument (N, &argc, argv); - initValue = atof (argv[N]); - remove_argument (N, &argc, argv); - initOutput = TRUE; - } if (argc != 11) { @@ -29,24 +15,22 @@ if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - int NX = in[0].matrix.Naxis[0]; - int NY = in[0].matrix.Naxis[1]; + NX = in[0].matrix.Naxis[0]; + NY = in[0].matrix.Naxis[1]; - int sx = atof (argv[3]); - int sy = atof (argv[4]); - int nx = atof (argv[5]); - int ny = atof (argv[6]); + sx = atof (argv[3]); + sy = atof (argv[4]); + nx = atof (argv[5]); + ny = atof (argv[6]); - int Sx = atof (argv[7]); - int Sy = atof (argv[8]); - int Nx = atof (argv[9]); - int Ny = atof (argv[10]); + Sx = atof (argv[7]); + Sy = atof (argv[8]); + Nx = atof (argv[9]); + Ny = atof (argv[10]); if ((Sy + ny > Ny) || (Sx + nx > Nx)) { - gprint (GP_ERR, "source pixels extend beyond target pixels\n"); + gprint (GP_ERR, "mismatch between source and dest regions\n"); gprint (GP_ERR, "%d + %d > %d or %d + %d > %d\n", Sy, ny, Ny, Sx, nx, Nx); - // return (FALSE); + return (FALSE); } - - // XXX : allow source region to fall outside source image /* region is not on first image */ @@ -55,10 +39,14 @@ (sy > in[0].matrix.Naxis[1])) { gprint (GP_ERR, "region outside of source image\n"); - // return (FALSE); + return (FALSE); } + if ((Sx + nx > Nx) || (Sy + ny > Ny)) { + gprint (GP_ERR, "source region larger than dest region\n"); + return (FALSE); + } if ((Sx < 0) || (Sy < 0)) { gprint (GP_ERR, "dest region out of range\n"); - // return (FALSE); + return (FALSE); } @@ -88,61 +76,13 @@ } - // (NX,NY) : source image dimensions - // (Nx,Ny) : target image dimensions - - // allow (sx, sy), (Sx, Sy) to be off image - // allow (sx+nx, sy+ny) to be off image - - // if (sx < 0) { - // nx += sx; - // sx = 0; - // } - // if (sx > NX) sx = NX; - // - // if (Sx < 0) { - // Nx += sx; - // sx = 0; - // } - // if (sx > NX) sx = NX; - - float *Vin = (float *)(in[0].matrix.buffer); - float *Vout = (float *)(out[0].matrix.buffer); - - if (initOutput) { - for (int j = 0; j < Ny; j++) { - for (int i = 0; i < Nx; i++) { - Vout[i + Nx*j] = initValue; - } - } - } - - int Nps = NX*NY; - int Npt = Nx*Ny; - - for (int j = 0; j < ny; j++) { - if (j + sy < 0) continue; // not yet on source image - if (j + Sy < 0) continue; // not yet on target image - if (j + sy >= NY) continue; // past edge of source image - if (j + Sy >= Ny) continue; // past edge of target image - - // Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx; - // Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx; - - for (int i = 0; i < nx; i++) { + for (j = 0; j < ny; j++) { + if (j + sy < 0) continue; + if (j + sy >= NY) continue; + Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx; + Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx; + for (i = 0; i < nx; i++, Vin++, Vout++) { if (i + sx < 0) continue; if (i + sx >= NX) continue; - if (i + Sx < 0) continue; - if (i + Sx >= Nx) continue; - - int ps = i + sx + (j + sy)*NX; - int pt = i + Sx + (j + Sy)*Nx; - - myAssert (pt >= 0, "oops 1"); - myAssert (pt < Npt, "oops 2"); - - myAssert (ps >= 0, "oops 3"); - myAssert (ps < Nps, "oops 4"); - - Vout[pt] = Vin[ps]; + *Vout = *Vin; } } @@ -151,2 +91,3 @@ } + Index: trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/init.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/init.c (revision 41340) @@ -3,5 +3,4 @@ int accum PROTO((int, char **)); -int antialias PROTO((int, char **)); int applyfit PROTO((int, char **)); int applyfit1d PROTO((int, char **)); @@ -68,5 +67,4 @@ int join PROTO((int, char **)); int jpeg PROTO((int, char **)); -int kapamemory PROTO((int, char **)); int kern PROTO((int, char **)); int keyword PROTO((int, char **)); @@ -81,10 +79,8 @@ int lookup PROTO((int, char **)); int matrix PROTO((int, char **)); -int match1d PROTO((int, char **)); int match2d PROTO((int, char **)); int mkrgb PROTO((int, char **)); int mcreate PROTO((int, char **)); int medacc PROTO((int, char **)); -int mgaussdev PROTO((int, char **)); int mget PROTO((int, char **)); int mget3d PROTO((int, char **)); @@ -103,5 +99,4 @@ int parity PROTO((int, char **)); int point PROTO((int, char **)); -int pdf PROTO((int, char **)); int ps PROTO((int, char **)); int vprint PROTO((int, char **)); @@ -116,5 +111,4 @@ int queuesubstr PROTO((int, char **)); int queuesize PROTO((int, char **)); -int queue2book PROTO((int, char **)); int rd PROTO((int, char **)); int rdseg PROTO((int, char **)); @@ -162,5 +156,4 @@ int vgauss PROTO((int, char **)); int vlorentz PROTO((int, char **)); -int vsigmoid PROTO((int, char **)); int vellipse PROTO((int, char **)); int vmaxwell PROTO((int, char **)); @@ -170,11 +163,7 @@ int vzload PROTO((int, char **)); int vstats PROTO((int, char **)); -int vstats PROTO((int, char **)); -int virls PROTO((int, char **)); -int vwtmean PROTO((int, char **)); int vroll PROTO((int, char **)); int vshift PROTO((int, char **)); int vpeaks PROTO((int, char **)); -int vtransitions PROTO((int, char **)); int vpop PROTO((int, char **)); int vsmooth PROTO((int, char **)); @@ -198,5 +187,4 @@ static Command cmds[] = { {1, "accum", accum, "accumulate vector values in another vector"}, - {1, "antialias", antialias, "set anti-alias sigma value for display"}, {1, "applyfit", applyfit1d, "apply 1-d fit to new vector"}, {1, "applyfit1d", applyfit1d, "apply 1-d fit to new vector"}, @@ -267,5 +255,4 @@ {1, "join", join, "find the join of two ID vectors"}, {1, "jpeg", jpeg, "convert display image to JPEG"}, - {1, "kapamemory", kapamemory, "manage kapa memory dump options"}, {1, "kern", kern, "convolve with 3x3 kernel"}, {1, "keyword", keyword, "extract a FITS keyword from image header"}, @@ -278,5 +265,4 @@ {1, "imcreate", mcreate, "create an image"}, {1, "medacc", medacc, "accumulate vector values in another vector"}, - {1, "mgaussdev", mgaussdev, "generate a gaussian deviate image"}, {1, "mget", mget, "extract a vector from an image"}, {1, "mget3d", mget3d, "extract a vector from a 3D image"}, @@ -287,5 +273,4 @@ {1, "medimage", medimage_command, "median image manipulation"}, {1, "matrix", matrix, "matrix math operations"}, - {1, "match1d", match1d, "match 2 vectors and return matched indexes"}, {1, "match2d", match2d, "match 2 pairs of X,Y vectors and return matched indexes"}, {1, "mkrgb", mkrgb, "convert 3 images to rgb jpeg (use Kapa for better control)"}, @@ -304,5 +289,4 @@ {1, "ppm", jpeg, "convert display graphic to PPM"}, {1, "ps", ps, "convert display to PostScript"}, - {1, "pdf", pdf, "convert display to PDF"}, {1, "print_vectors", vprint, "print a set of vectors"}, {1, "vprint", vprint, "print a set of vectors"}, @@ -317,6 +301,4 @@ {1, "queuesize", queuesize, "show queue size"}, {1, "queuesubstr", queuesubstr, "bulk replace strings in queue"}, - {1, "queue2book", queue2book, "convert queue with ipptool output to book"}, - {1, "ipptool2book", queue2book, "convert queue with ipptool output to book"}, {1, "rd", rd, "load fits image"}, {1, "rdseg", rdseg, "read a segment of an image from a file"}, @@ -365,5 +347,4 @@ {1, "vgauss", vgauss, "fit a Gaussian to a vector"}, {1, "vlorentz", vlorentz, "fit a Lorentzian to a vector"}, - {1, "vsigmoid", vsigmoid, "fit a Sigmoid to a vector"}, {1, "vellipse", vellipse, "fit a Ellipse to a vector pair"}, {1, "vgrid", vgrid, "generate an image from a triplet of vectors"}, @@ -373,6 +354,5 @@ {1, "vload", vload, "load vectors as overlay on image display"}, {1, "vmaxwell", vmaxwell, "fit a Maxwellian to a vector"}, - {1, "vpeaks", vpeaks, "find coord and flux of peaks in vector"}, - {1, "vtransitions", vtransitions, "find points in vector that cross the transition value"}, + {1, "vpeaks", vpeaks, "fine coord and flux of peaks in vector"}, {1, "vpop", vpop, "remove first element of a vector"}, {1, "vroll", vroll, "roll vector elements by 1 entry"}, @@ -380,6 +360,4 @@ {1, "vsmooth", vsmooth, "Gaussian smooth of a vector"}, {1, "vstats", vstats, "statistics on a vector"}, - {1, "vwtmean", vwtmean, "weighted mean of a vector"}, - {1, "virls", virls, "IRLS mean of a vector"}, {1, "vzload", vzload, "load vectors as overlay on image display (scaled points)"}, {1, "vsh", vsh, "Vector Spherical Harmonics"}, Index: trunk/Ohana/src/opihi/cmd.data/interpolate_presort.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/interpolate_presort.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/interpolate_presort.c (revision 41340) @@ -7,11 +7,4 @@ double x0, dx, dy, y0; Vector *xout, *yout, *xin, *yin; - - int N; - int FillEnds = FALSE; - if ((N = get_argument (argc, argv, "-fill-ends"))) { - FillEnds = TRUE; - remove_argument (N, &argc, argv); - } /** check basic syntax **/ @@ -43,12 +36,10 @@ y0 = yin[0].elements.Flt[0]; - // fill in the start with NANs - for (i = 0; (i < xout[0].Nelements) && (xout[0].elements.Flt[i] < xin[0].elements.Flt[0]); i++) { - yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[0] : NAN; - } + /* in and out vectors are sorted */ + j = 0; // every output pixel should get a value unless // we are below the lowest in or above the highest in - for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) { + for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) { yout[0].elements.Flt[i] = NAN; @@ -75,10 +66,10 @@ // fill in the rest with NANs - int NinLast = xin[0].Nelements - 1; while (i < yout[0].Nelements) { - yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN; + yout[0].elements.Flt[i] = NAN; i++; } return (TRUE); + } Index: trunk/Ohana/src/opihi/cmd.data/kapamemory.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/kapamemory.c (revision 41323) +++ (revision ) @@ -1,56 +1,0 @@ -# include "data.h" - -int kapamemory (int argc, char **argv) { - - int N, kapa; - - /* display source */ - char *name = NULL; - if ((N = get_argument (argc, argv, "-n"))) { - remove_argument (N, &argc, argv); - name = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if (!GetImage (NULL, &kapa, name)) return (FALSE); - FREE (name); - - if (argc < 2) { - gprint (GP_ERR, "USAGE: kapamemory (mode)\n"); - gprint (GP_ERR, " mode = dump, nlines, exit\n"); - return (FALSE); - } - - if (!strcasecmp(argv[1], "dump")) { - KapaMemoryDump (kapa); - return TRUE; - } - - if (!strcasecmp(argv[1], "exit")) { - if (argc != 3) goto exit_usage; - if (!strcasecmp(argv[2], "on")) { - KapaMemoryDumpOnExit (kapa, TRUE); - return TRUE; - } - if (!strcasecmp(argv[2], "off")) { - KapaMemoryDumpOnExit (kapa, FALSE); - return TRUE; - } - exit_usage: - gprint (GP_ERR, "USAGE: kapamemory exit (status)\n"); - gprint (GP_ERR, " statue = 'on' or 'off'\n"); - return FALSE; - } - - if (!strcasecmp(argv[1], "nlines")) { - if (argc != 3) { - gprint (GP_ERR, "USAGE: kapamemory nlines (N)\n"); - return FALSE; - } - int Nlines = atoi (argv[2]); - KapaMemoryDumpLines (kapa, Nlines); - return TRUE; - } - gprint (GP_ERR, "ERROR: unknown kapamemory mode %s\n", argv[1]); - return FALSE; -} - Index: trunk/Ohana/src/opihi/cmd.data/list_vectors.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/list_vectors.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/list_vectors.c (revision 41340) @@ -32,26 +32,16 @@ } - switch (vec->type) { - case OPIHI_FLT: - if (Variable) { - set_str_variable (Variable, "FLT"); - } else { - gprint (GP_LOG, "%s : FLT\n", argv[1]); - } - break; - case OPIHI_INT: - if (Variable) { - set_str_variable (Variable, "INT"); - } else { - gprint (GP_LOG, "%s : INT\n", argv[1]); - } - break; - case OPIHI_STR: - if (Variable) { - set_str_variable (Variable, "STR"); - } else { - gprint (GP_LOG, "%s : STR\n", argv[1]); - } - break; + if (vec->type == OPIHI_FLT) { + if (Variable) { + set_str_variable (Variable, "FLT"); + } else { + gprint (GP_LOG, "%s : FLT\n", argv[1]); + } + } else { + if (Variable) { + set_str_variable (Variable, "INT"); + } else { + gprint (GP_LOG, "%s : INT\n", argv[1]); + } } if (Variable) free (Variable); Index: trunk/Ohana/src/opihi/cmd.data/match1d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/match1d.c (revision 41323) +++ (revision ) @@ -1,248 +1,0 @@ -# include "data.h" - -int find_matches1d (Vector *X1, Vector *X2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch); -int find_matches1d_closest (Vector *X1, Vector *X2, double Radius, Vector *index); - -// match1d (X1) (X2) (Radius) [-index1 (index1)] [-index2 (index2)] [-nomatch1 nomatch1] [-nomatch2 nomatch2] -int match1d (int argc, char **argv) { - - int N, CLOSEST; - double Radius; - char *endptr; - Vector *X1vec, *X2vec; - Vector *index1, *index2; - - if ((N = get_argument (argc, argv, "-h"))) goto usage; - if ((N = get_argument (argc, argv, "--help"))) goto usage; - - CLOSEST = FALSE; - if ((N = get_argument (argc, argv, "-closest"))) { - remove_argument (N, &argc, argv); - CLOSEST = TRUE; - } - - if ((N = get_argument (argc, argv, "-index1"))) { - remove_argument (N, &argc, argv); - if ((index1 = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { - if ((index1 = SelectVector ("index1", ANYVECTOR, TRUE)) == NULL) return (FALSE); - } - - if ((N = get_argument (argc, argv, "-index2"))) { - remove_argument (N, &argc, argv); - if ((index2 = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { - if ((index2 = SelectVector ("index2", ANYVECTOR, TRUE)) == NULL) return (FALSE); - } - - Vector *radiusMatch = NULL; - if ((N = get_argument (argc, argv, "-radius"))) { - if (CLOSEST) { - gprint (GP_ERR, "error: -radius and -closest are currently incompatible\n"); - return (FALSE); - } - remove_argument (N, &argc, argv); - if ((radiusMatch = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } - - if (argc != 4) { - gprint (GP_ERR, "USAGE: match1d X1 X2 Radius [-index1 (index1)] [-index2 (index2)] [-closest]\n"); - gprint (GP_ERR, " use -h or --help for more detail\n"); - return (FALSE); - } - - /* - we have two modes of operation: - - without -closest, we are finding all matched pairs within the match radius. in this - case, the two index vectors have the same length, one entry per matched pair. - x1[index1],y1[index1] matches to x2[index2],y2[index2]. - - with -closest selected, we are finding the closest element of set 1 to each of set 2 - and vice versa. in this case, index1 is always the same length as x1,y1, while index2 - is the same lengths as x2,y2. x2[index1],y2[index1] matches x1,y1 while - x1[index2],y1[index2] matches x2,y2 - - */ - - if ((X1vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((X2vec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); - - REQUIRE_VECTOR_FLT (X1vec, FALSE); - REQUIRE_VECTOR_FLT (X2vec, FALSE); - - Radius = strtod (argv[3], &endptr); - if (*endptr) { - gprint (GP_ERR, "Radius must be numerical (%s)\n", argv[3]); - return (FALSE); - } - - if (CLOSEST) { - find_matches1d_closest (X1vec, X2vec, Radius, index1); - find_matches1d_closest (X2vec, X1vec, Radius, index2); - } else { - find_matches1d (X1vec, X2vec, Radius, index1, index2, radiusMatch); - } - return (TRUE); - -usage: - gprint (GP_ERR, "we have two modes of operation:\n\n"); - - gprint (GP_ERR, "without -closest, we are finding all matched pairs within the match radius. in this\n"); - gprint (GP_ERR, "case, the two index vectors have the same length, one entry per matched pair.\n"); - gprint (GP_ERR, "x1[index1] matches to x2[index2].\n\n"); - - gprint (GP_ERR, "with -closest selected, we are finding the closest element of set 1 to each of set 2\n"); - gprint (GP_ERR, "and vice versa. in this case, index1 is always the same length as x1, while index2\n"); - gprint (GP_ERR, "is the same lengths as x2. x2[index1] matches x1 while\n"); - gprint (GP_ERR, "x1[index2] matches x2\n\n"); - - gprint (GP_ERR, "if -index1 or -index2 is not supplied, the vectors are created with names index1 or index2\n"); - gprint (GP_ERR, "use 'reindex' to generate new vectors based on these index vectors\n"); - - gprint (GP_ERR, "if -radius (vector) is supplied, the vector will be filled with the distance between the matched pairs\n"); - gprint (GP_ERR, " not valid with -closest\n"); - return FALSE; -} - -// Radius is a 1-D separation -int find_matches1d (Vector *X1, Vector *X2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch) { - - off_t i, j, first_j, I, J, *N1, *N2, Nmatch, NMATCH, DMATCH; - double dX, dR; - - NMATCH = MAX(MAX(0.01*X1->Nelements, 0.01*X2->Nelements), 100); - DMATCH = NMATCH; - - ResetVector (index1, OPIHI_INT, NMATCH); - ResetVector (index2, OPIHI_INT, NMATCH); - if (radiusMatch) ResetVector (radiusMatch, OPIHI_FLT, NMATCH); - - ALLOCATE (N1, off_t, X1->Nelements); - ALLOCATE (N2, off_t, X2->Nelements); - - for (i = 0; i < X1->Nelements; i++) { N1[i] = i; } - for (i = 0; i < X2->Nelements; i++) { N2[i] = i; } - - dsort_indexonly (X1->elements.Flt, N1, X1->Nelements); - dsort_indexonly (X2->elements.Flt, N2, X2->Nelements); - - Nmatch = 0; - for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) { - I = N1[i]; - J = N2[j]; - - if (!isfinite(X1->elements.Flt[I])) { i++; continue; } - if (!isfinite(X2->elements.Flt[J])) { j++; continue; } - - dX = X1->elements.Flt[I] - X2->elements.Flt[J]; - - if (dX <= -1.02*Radius) { i++; continue; } - if (dX >= +1.02*Radius) { j++; continue; } - - // look for all matches of list2() to list1(i) - first_j = j; - for (j = first_j; (dX > -1.02*Radius) && (j < X2->Nelements); j++) { - J = N2[j]; - dX = X1->elements.Flt[I] - X2->elements.Flt[J]; - dR = fabs(dX); - if (dR < Radius) { - index1->elements.Int[Nmatch] = I; - index2->elements.Int[Nmatch] = J; - if (radiusMatch) radiusMatch->elements.Flt[Nmatch] = dR; - - // XXX track matches 1 and 2 with internal vector, save new nomatch index vectors - // after this loop - - Nmatch ++; - if (Nmatch >= NMATCH) { - NMATCH += DMATCH; - REALLOCATE (index1->elements.Int, opihi_int, NMATCH); - REALLOCATE (index2->elements.Int, opihi_int, NMATCH); - if (radiusMatch) { REALLOCATE (radiusMatch->elements.Flt, opihi_flt, NMATCH); } - } - } - } - j = first_j; - i++; - } - index1->Nelements = Nmatch; - index2->Nelements = Nmatch; - if (radiusMatch) radiusMatch->Nelements = Nmatch; - - free (N1); - free (N2); - - return (TRUE); -} - -// find the elements of X2,Y2 which are closest to each element of X1,Y1 (-1 if no match) -int find_matches1d_closest (Vector *X1, Vector *X2, double Radius, Vector *index) { - - off_t i, j, Jmin, Ji, I, J, *N1, *N2, NMATCH; - double dX, dR, Rmin; - - NMATCH = X1->Nelements; - ResetVector (index, OPIHI_INT, NMATCH); - - for (i = 0; i < index->Nelements; i++) { index->elements.Int[i] = -1; } - - ALLOCATE (N1, off_t, X1->Nelements); - ALLOCATE (N2, off_t, X2->Nelements); - - for (i = 0; i < X1->Nelements; i++) { N1[i] = i; } - for (i = 0; i < X2->Nelements; i++) { N2[i] = i; } - - dsort_indexonly (X1->elements.Flt, N1, X1->Nelements); - dsort_indexonly (X2->elements.Flt, N2, X2->Nelements); - - // find the closest entry in list 2 to the current entry in list 1: - for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) { - I = N1[i]; - J = N2[j]; - - dX = X1->elements.Flt[I] - X2->elements.Flt[J]; - - if (dX <= -1.02*Radius) { - // no match in list 2 to this entry - index->elements.Int[I] = -1; - i++; - continue; - } - if (dX >= +1.02*Radius) { j++; continue; } - - // look for closest matches of list2() to list1(i) - Jmin = -1; - for (Ji = j; (dX > -1.02*Radius) && (Ji < X2->Nelements); Ji++) { - J = N2[Ji]; - dX = X1->elements.Flt[I] - X2->elements.Flt[J]; - dR = fabs(dX); - if (dR > Radius) continue; - if (dR < Rmin) { - Rmin = dR; - Jmin = J; - } - } - - // no match in list 2 to this entry - if (Jmin == -1) { - index->elements.Int[I] = -1; - i++; - continue; - } - index->elements.Int[I] = Jmin; - i++; - } - index->Nelements = NMATCH; - - free (N1); - free (N2); - - return (TRUE); -} - - - Index: trunk/Ohana/src/opihi/cmd.data/medimage_calc.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/medimage_calc.c (revision 41323) +++ (revision ) @@ -1,365 +1,0 @@ -# include "data.h" - -float weight_cauchy_square_flt (float x2); -float irls_mean (float *val, float *wgt, int N, float *outvar, int *Npt); -float irls_fraction_interpolate (float *values, float fraction, int Npts); -float irls_robust_stdev (float *values, int Npts); -float irls_wtmean (float *val, float *wgt, int *idx, int Npt, float *variance); -void irls_bootstrap (int *idx, int Npts); -void irls_init (int N); -void irls_free (void); - -# define IRLS_TOLERANCE 1e-4 -int BOOTSTRAP = FALSE; -int BOOTSTRAP_NITER = 100; - -static float *irls_valsub = NULL; -static float *irls_wgtsub = NULL; -static int *irls_idx = NULL; -static float *irls_testval = NULL; - -enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN}; - -int medimage_calc (int argc, char **argv) { - - int ix, iy, n, N; - Buffer *output; - - BOOTSTRAP = FALSE; - if ((N = get_argument (argc, argv, "-bootstrap"))) { - BOOTSTRAP = TRUE; - remove_argument (N, &argc, argv); - } - BOOTSTRAP_NITER = 100; - if ((N = get_argument (argc, argv, "-bootstrap-iter"))) { - remove_argument (N, &argc, argv); - BOOTSTRAP_NITER = atoi (argv[N]); - remove_argument (N, &argc, argv); - BOOTSTRAP = TRUE; - } - - int mode = CALC_MEDIAN; - if ((N = get_argument (argc, argv, "-mean"))) { - mode = CALC_MEAN; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-irls"))) { - if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; } - mode = CALC_IRLS; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-wtmean"))) { - if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; } - mode = CALC_WTMEAN; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-median"))) { - if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; } - mode = CALC_MEDIAN; - remove_argument (N, &argc, argv); - } - - Buffer *variance = NULL; - if ((N = get_argument (argc, argv, "-variance"))) { - remove_argument (N, &argc, argv); - if ((variance = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } - - if (argc != 3) { - gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean,-irls,-wtmean,-median] [-variance output] [-bootstrap] [-bootstrap-iter (N)]\n"); - gprint (GP_ERR, " calculate the median image for the median image set\n"); - gprint (GP_ERR, " This function calculates the 'average' image for the supplied set of images.\n"); - gprint (GP_ERR, " Four options are available to calculate the average:\n"); - gprint (GP_ERR, " -mean (straight arithmetic mean)\n"); - gprint (GP_ERR, " -median\n"); - gprint (GP_ERR, " -wtmean (arithmetic mean, weighted by supplied variances)\n"); - gprint (GP_ERR, " -irls (iteratively-reweighted least-squares)\n"); - gprint (GP_ERR, " If -variance is supplied the returned image represents the estimate of the variance on the average image\n"); - gprint (GP_ERR, " For the four average methods above, the sqrt of the variance has the following meanings:\n"); - gprint (GP_ERR, " -mean : standard deviation / sqrt(Npts)\n"); - gprint (GP_ERR, " -median : standard deviation / sqrt(Npts)\n"); - gprint (GP_ERR, " -wtmean : formal error on the weighted mean [sum (1 / variances)]\n"); - gprint (GP_ERR, " -irls : formal error on the pixels used (excluding ~3 sigma outliers)\n"); - gprint (GP_ERR, " note that the irls formal error slightly underestimates the observed standard deviation\n"); - return FALSE; - } - - MedImageType *median = FindMedImage (argv[1]); - if (!median) { - gprint (GP_ERR, "median image %s not found\n", argv[1]); - return FALSE; - } - - if ((output = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE); - - int Ninput = median->Ninput; - int Nx = median->Nx; - int Ny = median->Ny; - - ALLOCATE_PTR (val, float, Ninput); - ALLOCATE_PTR (wgt, float, Ninput); - - if (mode == CALC_IRLS) irls_init (Ninput); - - gfits_free_matrix (&output->matrix); - gfits_free_header (&output->header); - if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE; - - float *outvalue = (float *) output->matrix.buffer; - float *varvalue = NULL; - if (variance) { - gfits_free_matrix (&variance->matrix); - gfits_free_header (&variance->header); - if (!CreateBuffer (variance, Nx, Ny, -32, 0.0, 1.0)) return FALSE; - - varvalue = (float *) variance->matrix.buffer; - } - - // save the number of points per pixel - Buffer *nptbuf = SelectBuffer ("irls_npt", ANYBUFFER, TRUE); - gfits_free_matrix (&nptbuf->matrix); - gfits_free_header (&nptbuf->header); - if (!CreateBuffer (nptbuf, Nx, Ny, -32, 0.0, 1.0)) return FALSE; - float *NptVal = (float *) nptbuf->matrix.buffer; - - for (iy = 0; iy < Ny; iy++) { - for (ix = 0; ix < Nx; ix++) { - - int N = 0; - int Npix = ix + Nx*iy; - for (n = 0; n < Ninput; n++) { - float v = median->flx[n][Npix]; - if (!isfinite(v)) continue; - val[N] = v; - wgt[N] = 1.0; - if (median->var[n]) { - float s = median->var[n][Npix]; - if (!isfinite(s)) continue; - if (fabs(s) < 2*FLT_MIN) s = 2*FLT_MIN; - wgt[N] = 1.0 / s; - } - N++; - } - if (N == 0) continue; - - switch (mode) { - case CALC_MEDIAN: - fsort (val, N); - outvalue[Npix] = val[(int)(0.5*N)]; - if (varvalue) { - float sum = 0.0; - for (n = 0; n < N; n++) { - sum += SQ(val[n] - outvalue[Npix]); - } - // variance on the mean (stdev / sqrt(N))^2 - varvalue[Npix] = sum / (N - 1) / N; - } - break; - case CALC_MEAN: { - float sum = 0.0; - for (n = 0; n < N; n++) { - sum += val[n]; - } - outvalue[Npix] = sum / (float) N; - - if (varvalue) { - float Sum = 0.0; - for (n = 0; n < N; n++) { - Sum += SQ(val[n] - outvalue[Npix]); - } - // variance on the mean (stdev / sqrt(N))^2 - varvalue[Npix] = Sum / (N - 1) / N; - } - break; - } - case CALC_WTMEAN: { - float variance; - outvalue[Npix] = irls_wtmean (val, wgt, NULL, N, &variance); - if (varvalue) { varvalue[Npix] = variance; } - break; - } - case CALC_IRLS: - if (varvalue) { - int Npts = 0; - outvalue[Npix] = irls_mean (val, wgt, N, &varvalue[Npix], &Npts); - NptVal[Npix] = Npts; - } else { - int Npts = 0; - outvalue[Npix] = irls_mean (val, wgt, N, NULL, &Npts); - NptVal[Npix] = Npts; - } - } - } - } - irls_free(); - return TRUE; -} - -// allocate temporary vectors -void irls_init (int N) { - - ALLOCATE (irls_valsub, float, N); - ALLOCATE (irls_wgtsub, float, N); - ALLOCATE (irls_idx, int, N); - ALLOCATE (irls_testval, float, BOOTSTRAP_NITER); -} - -// free temporary vectors -void irls_free (void) { - FREE (irls_valsub); - FREE (irls_wgtsub); - FREE (irls_idx); - FREE (irls_testval); - irls_valsub = NULL; - irls_wgtsub = NULL; - irls_idx = NULL; - irls_testval = NULL; -} - -float irls_mean (float *val, float *wgt, int N, float *outvar, int *Npt) { - - // calculate weighted mean - float Value = irls_wtmean (val, wgt, NULL, N, NULL); - - int converged = FALSE; - for (int i = 0; (i < 10) && !converged; i++) { - float ValueLast = Value; - float S1 = 0.0, S2 = 0.0; - - // calculate weight modification based on distances (squared). - // use modifier to calculate new weighted mean - for (int n = 0; n < N; n++) { - float dV = (val[n] - Value); - float d2 = SQ(dV) * wgt[n]; - - float Mod = weight_cauchy_square_flt (d2); - S1 += Mod * wgt[n] * val[n]; - S2 += Mod * wgt[n]; - } - Value = S1 / S2; - - float delta = fabs(Value - ValueLast); - if (delta < Value * IRLS_TOLERANCE) converged = TRUE; - } - - if (outvar) { - if (BOOTSTRAP) { - // generate a subset vector of just the accepted points - - // calculate stdev of high-weight points - int npt = 0; - for (int n = 0; n < N; n++) { - - float dV = (val[n] - Value); - float d2 = SQ(dV) * wgt[n]; - - float Mod = weight_cauchy_square_flt (d2); - if (Mod < 0.1) continue; // totally ad-hoc number - - irls_valsub[npt] = val[n]; - irls_wgtsub[npt] = wgt[n]; - npt ++; - } - - for (int iter = 0; iter < BOOTSTRAP_NITER; iter++) { - irls_bootstrap (irls_idx, npt); // fill idx with the index of the resampled points - irls_testval[iter] = irls_wtmean (irls_valsub, irls_wgtsub, irls_idx, npt, NULL); - } - - float sigma = irls_robust_stdev (irls_testval, BOOTSTRAP_NITER); - - *outvar = SQ(sigma); - *Npt = npt; - } else { - // calculate stdev of high-weight points - float S1 = 0.0, S2 = 0.0; - int npt = 0; - for (int n = 0; n < N; n++) { - - float dV = (val[n] - Value); - float d2 = SQ(dV) * wgt[n]; - - float Mod = weight_cauchy_square_flt (d2); - if (Mod < 0.1) continue; // totally ad-hoc number - - S1 += SQ(dV); - S2 += wgt[n]; - npt ++; - } - *Npt = npt; - *outvar = 1 / S2; - } - } - - return Value; -} - -void irls_bootstrap (int *idx, int Npts) { - // generate an index (idx) containing the index values - // for the Npts randomly resampled points - for (int i = 0; i < Npts; i++) { - idx[i] = Npts * drand48(); - } -} - -float irls_wtmean (float *val, float *wgt, int *idx, int Npt, float *variance) { - float S1 = 0.0, S2 = 0.0; - for (int n = 0; n < Npt; n++) { - int N = idx ? idx[n] : n; - S1 += wgt[N] * val[N]; - S2 += wgt[N]; - } - float Value = S1 / S2; - if (variance) *variance = 1.0 / S2; - return Value; -} - -float irls_robust_stdev (float *values, int Npts) { - - fsort (values, Npts); - - float Slo = irls_fraction_interpolate (values, 0.158655, Npts); - float Shi = irls_fraction_interpolate (values, 0.841345, Npts); - float sigma = (Shi - Slo) / 2.0; - - return sigma; -} - -float irls_fraction_interpolate (float *values, float fraction, int Npts) { - - float F = fraction * Npts; - int N = fraction * Npts; - - if (N < 0 ) return NAN; - if (N >= Npts - 2) return NAN; - - // interpolate between N,N+1 - - float S = (F - N) * (values[N+1] - values[N]) + values[N]; - return S; -} - - -// exp(-(x^2/s^2)/2) = (1/2) -// -(x^2/s^2)/2 = ln(1/2) -// (x^2/s^2)/2 = ln(2) -// (x^2/s^2) = 2ln(2) -// (x /s) = sqrt(2ln(2)) : half-width at half-max -// FWHM = 2sqrt(2ln(2)) - -// R2 = (X / 2.385)^2 = (X^2 / 2.385^2) - -# define CAUCY_FACTOR 2.0 - -/* -float weight_cauchy_square_flt (float x2) { - return (1.0); -} -*/ - -float weight_cauchy_square_flt (float x2) { - float r2 = x2 / CAUCY_FACTOR; - return (1.0 / (1.0 + r2)); -} - Index: trunk/Ohana/src/opihi/cmd.data/medimage_commands.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/medimage_commands.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/medimage_commands.c (revision 41340) @@ -15,30 +15,13 @@ int medimage_add (int argc, char **argv) { - int N; Buffer *image; - Buffer *var = NULL; - - if ((N = get_argument (argc, argv, "-variance"))) { - remove_argument (N, &argc, argv); - if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } - - if (argc != 3) { - gprint (GP_ERR, "USAGE: medimage add (name) (image) [-variance variance]\n"); + + if (argc != 3) { + gprint (GP_ERR, "USAGE: medimage add (name) (image)\n"); gprint (GP_ERR, " add the given image to the set of images to be medianed\n"); - gprint (GP_ERR, " optionally supply variance image (for weighted calculations)\n"); return FALSE; } if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE); - - if (var) { - if ((var->matrix.Naxis[0] != image->matrix.Naxis[0]) || - (var->matrix.Naxis[1] != image->matrix.Naxis[1])) { - gprint (GP_ERR, "variance buffer does not match image buffer dimensions\n"); - return FALSE; - } - } MedImageType *median = FindMedImage (argv[1]); @@ -63,17 +46,76 @@ // new image should match existing medimage dimensions - // AddMedImage (median, image, var); int Ninput = median->Ninput; median->Ninput ++; - REALLOCATE (median->flx, float *, median->Ninput); - REALLOCATE (median->var, float *, median->Ninput); - - ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny); - memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny); - - median->var[Ninput] = NULL; - if (var) { - ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny); - memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny); + REALLOCATE (median->buffers, float *, median->Ninput); + + ALLOCATE (median->buffers[Ninput], float, median->Nx*median->Ny); + memcpy (median->buffers[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny); + + return TRUE; +} + +int medimage_calc (int argc, char **argv) { + + int ix, iy, n, N; + Buffer *output; + + int CALC_MEAN = FALSE; + if ((N = get_argument (argc, argv, "-mean"))) { + CALC_MEAN = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 3) { + gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean]\n"); + gprint (GP_ERR, " calculate the median image for the median image set\n"); + return FALSE; + } + + MedImageType *median = FindMedImage (argv[1]); + if (!median) { + gprint (GP_ERR, "median image %s not found\n", argv[1]); + return FALSE; + } + + if ((output = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE); + + int Ninput = median->Ninput; + int Nx = median->Nx; + int Ny = median->Ny; + + float *value = NULL; + ALLOCATE (value, float, Ninput); + + gfits_free_matrix (&output->matrix); + gfits_free_header (&output->header); + if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE; + + float *outvalue = (float *) output->matrix.buffer; + + for (iy = 0; iy < Ny; iy++) { + for (ix = 0; ix < Nx; ix++) { + + int N = 0; + int Npix = ix + Nx*iy; + for (n = 0; n < Ninput; n++) { + float v = median->buffers[n][Npix]; + if (!isfinite(v)) continue; + value[N] = v; + N++; + } + if (N == 0) continue; + + if (CALC_MEAN) { + float sum = 0.0; + for (n = 0; n < N; n++) { + sum += value[n]; + } + outvalue[Npix] = sum / (float) N; + } else { + fsort (value, N); + outvalue[Npix] = value[(int)(0.5*N)]; + } + } } @@ -82,39 +124,39 @@ /* - int medimage_save (int argc, char **argv) { - - int N; - - int APPEND = FALSE; - if ((N = get_argument (argc, argv, "-append"))) { - APPEND = TRUE; - remove_argument (N, &argc, argv); - } - - if (argc != 3) { - gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n"); - return FALSE; - } - - if (!SaveMedImage(argv[2], argv[1], APPEND)) { - gprint (GP_ERR, "failed to save medimage %s\n", argv[1]); - return (FALSE); - } - return TRUE; - } - - int medimage_load (int argc, char **argv) { - - if (argc != 3) { - gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n"); - return FALSE; - } - - if (!LoadMedImage(argv[2], argv[1])) { - gprint (GP_ERR, "failed to load medimage %s\n", argv[1]); - return (FALSE); - } - return TRUE; - } +int medimage_save (int argc, char **argv) { + + int N; + + int APPEND = FALSE; + if ((N = get_argument (argc, argv, "-append"))) { + APPEND = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 3) { + gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n"); + return FALSE; + } + + if (!SaveMedImage(argv[2], argv[1], APPEND)) { + gprint (GP_ERR, "failed to save medimage %s\n", argv[1]); + return (FALSE); + } + return TRUE; +} + +int medimage_load (int argc, char **argv) { + + if (argc != 3) { + gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n"); + return FALSE; + } + + if (!LoadMedImage(argv[2], argv[1])) { + gprint (GP_ERR, "failed to load medimage %s\n", argv[1]); + return (FALSE); + } + return TRUE; +} */ Index: trunk/Ohana/src/opihi/cmd.data/mgaussdev.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/mgaussdev.c (revision 41323) +++ (revision ) @@ -1,34 +1,0 @@ -# include "data.h" - -int mgaussdev (int argc, char **argv) { - - Buffer *buf; - - if (argc != 6) goto usage; - - if ((buf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE); - int Nx = atof (argv[2]); - int Ny = atof (argv[3]); - - double mean = atof (argv[4]); - double sigma = atof (argv[5]); - - /* I should encapsulate this in a create_default_buffer */ - gfits_free_matrix (&buf[0].matrix); - gfits_free_header (&buf[0].header); - - // 3D CUBE OPTION: if (!CreateBuffer3D (buf, Nx, Ny, Nz, -32, 1.0, 0.0)) return FALSE; - if (!CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0)) return FALSE; - - ohana_gaussdev_init (); - - float *v = (float *) buf[0].matrix.buffer; - for (int i = 0; i < Nx*Ny; i++, v++) { - *v = ohana_gaussdev_rnd (mean, sigma); - } - return (TRUE); - - usage: - gprint (GP_ERR, "USAGE: mgaussdev (buff) Nx Ny mean sigma\n"); - return (FALSE); -} Index: trunk/Ohana/src/opihi/cmd.data/pdf.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/pdf.c (revision 41323) +++ (revision ) @@ -1,73 +1,0 @@ -# include "data.h" - -int pdf (int argc, char **argv) { - - char filename[1024], pagename[1024], *name; - int N, kapa, scaleMode, pageMode; - - if ((N = get_argument (argc, argv, "--help"))) goto help; - if ((N = get_argument (argc, argv, "-h"))) goto help; - - pageMode = KAPA_PS_NEWPLOT; - - /* new page? */ - /* - if ((N = get_argument (argc, argv, "-newpage"))) { - remove_argument (N, &argc, argv); - pageMode = KAPA_PS_NEWPAGE; - } - if ((N = get_argument (argc, argv, "-raw"))) { - remove_argument (N, &argc, argv); - pageMode = KAPA_PS_RAWPAGE; - } - */ - - /* scale image? */ - scaleMode = TRUE; - if ((N = get_argument (argc, argv, "-noscale"))) { - remove_argument (N, &argc, argv); - scaleMode = FALSE; - } - - /* what file? */ - filename[0] = 0; - if ((N = get_argument (argc, argv, "-name"))) { - remove_argument (N, &argc, argv); - strcpy (filename, argv[N]); - remove_argument (N, &argc, argv); - } - /* pagename ? */ - strcpy (pagename, "default"); - if ((N = get_argument (argc, argv, "-pagename"))) { - remove_argument (N, &argc, argv); - strcpy (pagename, argv[N]); - remove_argument (N, &argc, argv); - } - - /* which tool */ - name = NULL; - if ((N = get_argument (argc, argv, "-n"))) { - remove_argument (N, &argc, argv); - name = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - - if ((argc > 1) && filename[0]) goto help; - - if (argc > 1) strcpy (filename, argv[1]); - - // get the connection to kapa, false if none available - if (!GetGraphdata (NULL, &kapa, name)) return (FALSE); - FREE (name); - - if (!filename[0]) strcpy (filename, "kapa.ps"); - - /* tell Ximage/Xgraph to ps the image */ - KapaPDF (kapa, filename, scaleMode, pageMode, pagename); - return (TRUE); - -help: - gprint (GP_ERR, "USAGE: ps [-name file.ps] [-g | -i] [-n device] [-raw] [-noscale] [-newpage] [-pagename (name]\n"); - return (FALSE); -} - Index: trunk/Ohana/src/opihi/cmd.data/peak.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/peak.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/peak.c (revision 41340) @@ -39,5 +39,4 @@ if (*X < start) continue; if (*X > end) continue; - if (!isfinite(*Y)) continue; if (*Y < ymax) continue; xmax = *X; @@ -64,5 +63,4 @@ if (*X < start) continue; if (*X > end) continue; - if (!isfinite(*Y)) continue; if (*Y < ymax) continue; xmax = *X; Index: trunk/Ohana/src/opihi/cmd.data/print_vectors.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/print_vectors.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/print_vectors.c (revision 41340) @@ -56,11 +56,7 @@ if (vec[i][0].type == OPIHI_FLT) { gprint (GP_LOG, "%f ", vec[i][0].elements.Flt[j]); - } - if (vec[i][0].type == OPIHI_INT) { + } else { gprint (GP_LOG, OPIHI_INT_FMT" ", vec[i][0].elements.Int[j]); - } - if (vec[i][0].type == OPIHI_STR) { - gprint (GP_LOG, "%s ", vec[i][0].elements.Str[j]); - } + } } } Index: trunk/Ohana/src/opihi/cmd.data/queue2book.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/queue2book.c (revision 41323) +++ (revision ) @@ -1,248 +1,0 @@ -# include "pantasks.h" - -# define FREEKEYS \ - if (keys) { for (i = 0; i < Nkeys; i++) { free (keys[i]); } free (keys); } \ - for (i = 0; i < setWordN; i++) { free (setWordList[i]); free (setWordValue[i]); } \ - free (setWordValue); free (setWordList); - -// convert the named queue with ipptool metadata output into book format -int queue2book (int argc, char **argv) { - - int i, N, onPage, found, Unique, Nkeys, NKEYS; - char *line, *tmpword, *tmpvalue; - char pagename[512]; // XXX this should be made dynamic, though it is an unlikey problem - char *bookName, **keys, *p, *q; - char **setWordList; - char **setWordValue; - int setWordN, setWordNalloc; - - Page *page = NULL; - Book *book = NULL; - Queue *queue = NULL; - - /* supply additional constant words */ - setWordN = 0; - setWordNalloc = 10; - ALLOCATE (setWordList, char *, setWordNalloc); - ALLOCATE (setWordValue, char *, setWordNalloc); - while ((N = get_argument (argc, argv, "-setword"))) { - remove_argument (N, &argc, argv); - setWordList[setWordN] = strcreate (argv[N]); - remove_argument (N, &argc, argv); - setWordValue[setWordN] = strcreate (argv[N]); - remove_argument (N, &argc, argv); - setWordN ++; - if (setWordN >= setWordNalloc) { - setWordNalloc += 10; - REALLOCATE (setWordList, char *, setWordNalloc); - REALLOCATE (setWordValue, char *, setWordNalloc); - } - } - - Unique = FALSE; - if ((N = get_argument (argc, argv, "-uniq"))) { - remove_argument (N, &argc, argv); - Unique = TRUE; - } - - Nkeys = 0; - keys = NULL; - if ((N = get_argument (argc, argv, "-key"))) { - remove_argument (N, &argc, argv); - - /* parse the key (word:word:word) into constituents */ - NKEYS = 10; - ALLOCATE (keys, char *, NKEYS); - - p = q = argv[N]; - while (q != NULL) { - q = strchr (p, ':'); - keys[Nkeys] = (q == NULL) ? strcreate (p) : strncreate (p, q - p); - p = q + 1; - Nkeys ++; - CHECK_REALLOCATE (keys, char *, NKEYS, Nkeys, 10); - } - remove_argument (N, &argc, argv); - } - - if (argc != 3) { - gprint (GP_ERR, "USAGE: queue2book (queue) (book) [-uniq] [-key key[:key..]]\n"); - gprint (GP_ERR, " converts the named queue with ipptool metadata output into book format\n"); - FREEKEYS; - return (TRUE); - } - - queue = FindQueue (argv[1]); - if (queue == NULL) { - gprint (GP_ERR, "queue %s not found\n", argv[1]); - FREEKEYS; - return (TRUE); - } - - book = CreateBook (argv[2]); - - /* the first non-whitespace line of the queue should contain: - BookName MULTI - */ - - /* scan for first non-emtpy line */ - while (TRUE) { - line = PopQueue (queue); - if (!line) { - FREEKEYS; - return (TRUE); // no output in queue - not an error? - } - stripwhite (line); - if (line[0] == '#') { - free (line); - continue; - } - if (*line) break; - free (line); - } - bookName = thisword (line); - tmpword = nextword (line); - - if (!tmpword || strcmp(tmpword, "MULTI")) { - gprint (GP_ERR, "ERROR: missing metadata output name on first line\n"); - free (bookName); - free (line); - FREEKEYS; - return (TRUE); - } - free (line); - - onPage = FALSE; - while (TRUE) { - line = PopQueue (queue); - if (!line) { - free (bookName); - if (onPage) { - gprint (GP_ERR, "ERROR: unterminated metadata\n"); - FREEKEYS; - return (TRUE); - } - FREEKEYS; - return (TRUE); - } - stripwhite (line); - if (!*line) { - free (line); - continue; - } - if (line[0] == '#') { - free (line); - continue; - } - - tmpword = thisword (line); - if (tmpword == NULL) abort(); // should not happen if line exists - - if (!onPage) { - if (strcmp(tmpword, bookName)) { - gprint (GP_ERR, "ERROR: unexpected metadata name %s\n", tmpword); - free (line); - free (tmpword); - free (bookName); - FREEKEYS; - return (TRUE); - } - free (tmpword); - tmpword = thisword (nextword (line)); - if (strcmp(tmpword, "METADATA")) { - gprint (GP_ERR, "ERROR: missing METADATA tag\n"); - free (line); - free (tmpword); - free (bookName); - FREEKEYS; - return TRUE; - } - onPage = TRUE; - // we don't know the page name until we load its data? - sprintf (pagename, "page.%03d", book[0].Npages); - page = CreatePage (book, pagename); - free (line); - free (tmpword); - continue; - } - - /* close out the page */ - if (!strcmp(tmpword, "END")) { - // XXX set the page name based on the contents - free (tmpword); - free (line); - - if (keys) { - memset (pagename, 0, 512); - for (i = 0; i < Nkeys; i++) { - tmpword = BookGetWord (page, keys[i]); - if (tmpword == NULL) { - gprint (GP_ERR, "ERROR: missing key %s for ID\n", keys[i]); - FREEKEYS; - return (TRUE); - } - if (i > 0) strcat (pagename, ":"); - strcat (pagename, tmpword); - } - free (page[0].name); - page[0].name = strcreate (pagename); - } - - // add the fixed words to the page - for (i = 0; i < setWordN; i++) { - BookSetWord (page, setWordList[i], setWordValue[i]); - } - - if (Unique) { - /* search for an existing page with this name - delete this one if one is found */ - /* XXX don't attach this page to the book until this step is passed */ - - found = FALSE; - for (i = 0; !found && (i < book[0].Npages); i++) { - if (book[0].pages[i] == page) continue; - if (!strcmp (book[0].pages[i][0].name, page[0].name)) { - found = TRUE; - } - } - if (found) { - DeletePage (book, page); - } - } - - onPage = FALSE; - continue; - } - - tmpvalue = thisword (nextword(nextword(line))); - if (tmpvalue == NULL) { - gprint (GP_ERR, "ERROR: missing value for %s\n", tmpword); - free (tmpword); - free (line); - free (bookName); - FREEKEYS; - return (TRUE); - } - BookSetWord (page, tmpword, tmpvalue); - free (tmpword); - free (tmpvalue); - free (line); - } - FREEKEYS; - return (TRUE); -} - -/* scan through the queue lines. these should look like: - -bookName METADATA -key type value -key type value -... -END - -bookName METADATA -key type value -key type value -... -END -*/ Index: trunk/Ohana/src/opihi/cmd.data/read_vectors.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/read_vectors.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/read_vectors.c (revision 41340) @@ -35,5 +35,5 @@ // vector types -enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_DATE, COLTYPE_CHAR, COLTYPE_STR, COLTYPE_HMS}; +enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_CHAR}; static int Nvec = 0; @@ -91,11 +91,8 @@ gprint (GP_ERR, " -skip N : skip N lines before reading\n"); gprint (GP_ERR, " column names may include a type: name:type\n"); - gprint (GP_ERR, " type is int, float, char, str, time, date, hms\n"); + gprint (GP_ERR, " type is int, float, char, time\n"); gprint (GP_ERR, " for char, values are placed into a list $name:0 - $name:n\n"); - gprint (GP_ERR, " for str, values are placed into a string-typed vector\n"); - gprint (GP_ERR, " for hms, values are sexigesimal values HH:MM:SS (good for degrees as well)\n"); - gprint (GP_ERR, " for date, values are human-readable date strings YYYY/MM/DD\n"); gprint (GP_ERR, " for time, values are human-readable date/time strings YYYY/MM/DD,hh:mm:ss\n"); - gprint (GP_ERR, " date & time values are stored as a float using the TIMEFORMAT, TIMEREF values for the conversion\n"); + gprint (GP_ERR, " the resulting vector is a float based on the TIMEFORMAT, TIMEREF values\n"); gprint (GP_ERR, " -fits options:\n"); @@ -155,8 +152,5 @@ if (!strcasecmp(ptr, "int")) { coltype[i] = COLTYPE_INT; } if (!strcasecmp(ptr, "char")) { coltype[i] = COLTYPE_CHAR; } - if (!strcasecmp(ptr, "str")) { coltype[i] = COLTYPE_STR; } - if (!strcasecmp(ptr, "date")) { coltype[i] = COLTYPE_DATE; } if (!strcasecmp(ptr, "time")) { coltype[i] = COLTYPE_TIME; } - if (!strcasecmp(ptr, "hms")) { coltype[i] = COLTYPE_HMS; } if (!coltype[i]) goto bad_colname; } @@ -216,11 +210,7 @@ break; case COLTYPE_FLT: - case COLTYPE_HMS: case COLTYPE_TIME: - case COLTYPE_DATE: + // note that COLTYPE_TIME is a type of float ResetVector (vec[i], OPIHI_FLT, NELEM); - break; - case COLTYPE_STR: - ResetVector (vec[i], OPIHI_STR, NELEM); break; case COLTYPE_CHAR: @@ -297,5 +287,4 @@ time_t tvalue; int readStatus = FALSE; - int dataStatus = FALSE; // need to make the if cases for coltype[i] switch (coltype[i]) { @@ -322,20 +311,4 @@ break; } - case COLTYPE_STR: - { - // I need to get an isolated word in 'value' with the string value of this field - char *ptr = IsCSV ? ptrparse_csv (col[i], c0) : ptrparse (col[i], c0); - char *value = NULL; - if (IsCSV) { - char *end = parse_nextword_csv (ptr); - if (end) { - value = end ? strncreate (ptr, end - ptr) : strcreate (ptr); - } - } else { - value = thisword(ptr); - } - vec[i][0].elements.Str[Nelem] = value; // needs to be freed later. - break; - } case COLTYPE_FLT: readStatus = IsCSV ? dparse_csv (&dvalue, col[i], c0) : dparse (&dvalue, col[i], c0); @@ -347,29 +320,4 @@ vec[i][0].elements.Flt[Nelem] = readStatus ? dvalue : NAN; break; - case COLTYPE_DATE: { - char *string = NULL; - if (IsCSV) { - char *ptr = ptrparse_csv (col[i], c0); - string = strcreate (ptr); // create separate string (NULL-safe) - ptr = (string == NULL) ? NULL : strchr (string, ','); - if (ptr != NULL) *ptr = 0; // place an EOL here so parsing does not go past the comma - } else { - char *ptr = ptrparse (col[i], c0); // NULL-safe - string = getword (ptr); // NULL-safe - } - tvalue = ohana_date_to_sec (string); // NULL-safe (returns 0) - dvalue = TimeValue (tvalue, TimeReference, TimeFormat); - vec[i][0].elements.Flt[Nelem] = string ? dvalue : NAN; - FREE (string); - break; - } - case COLTYPE_HMS: { - char *ptr = IsCSV ? ptrparse_csv (col[i], c0) : ptrparse (col[i], c0); - char *string = strcreate (ptr); // make a copy so we can munge it in ohana_dms_to_ddd - dataStatus = ohana_dms_to_ddd (&dvalue, string); - vec[i][0].elements.Flt[Nelem] = ((string != NULL) && dataStatus) ? dvalue : NAN; - FREE (string); - break; - } } if (!readStatus && VERBOSE) { @@ -397,11 +345,6 @@ break; case COLTYPE_FLT: - case COLTYPE_HMS: case COLTYPE_TIME: - case COLTYPE_DATE: ResetVector (vec[i], OPIHI_FLT, NELEM); - break; - case COLTYPE_STR: - ResetVector (vec[i], OPIHI_STR, NELEM); break; case COLTYPE_CHAR: @@ -423,11 +366,6 @@ break; case COLTYPE_FLT: - case COLTYPE_HMS: case COLTYPE_TIME: - case COLTYPE_DATE: ResetVector (vec[i], OPIHI_FLT, Nelem); - break; - case COLTYPE_STR: - ResetVector (vec[i], OPIHI_STR, Nelem); break; case COLTYPE_CHAR: @@ -531,10 +469,6 @@ } - // by default, we now store a string-type field as a string-type vector. - - // if CharAsList is selected (and My < 10000), char fields will be saved as $NAME:0 - $NAME:m for m rows - - // if CharAsVectors is selected, char fields will be saved as vectors NAME:0 -- NAME:n for n characters - + // if CharAsVectors, char fields will be saved as vectors NAME:0 -- NAME:n for n characters + // else char fields will be saved as $NAME:0 - $NAME:m for m rows // if (Ny > 10000), force CharAsVectors int CharAsVectors = FALSE; @@ -542,9 +476,4 @@ remove_argument (N, &argc, argv); CharAsVectors = TRUE; - } - int CharAsList = FALSE; - if ((N = get_argument (argc, argv, "-char-list"))) { - remove_argument (N, &argc, argv); - CharAsList = TRUE; } @@ -722,46 +651,26 @@ if (!FITS_TRANSPOSE) { // read string column into a list rather than a vector - if (!strcmp (type, "char")) { - // save char-type field as a List: - if (CharAsList && (Ny < 3000)) { - char *fieldName = argv[i]; - char *Ptr = data; - char varname[1024]; // used as a buffer for the names of string fields - for (j = 0; j < Ny; j++) { - set_list_varname (varname, fieldName, j, FALSE); - char *value = strncreate (&Ptr[j*Nval], Nval); - // replace instances of $ with _ - char *p = strchr (value, '$'); - while (p) { - *p = '_'; - p = strchr (p, '$'); - } - set_str_variable (varname, value); - free (value); + if (!strcmp (type, "char") && !CharAsVectors && (Ny < 3000)) { + char *fieldName = argv[i]; + char *Ptr = data; + char varname[1024]; // used as a buffer for the names of string fields + for (j = 0; j < Ny; j++) { + set_list_varname (varname, fieldName, j, FALSE); + char *value = strncreate (&Ptr[j*Nval], Nval); + // replace instances of $ with _ + char *p = strchr (value, '$'); + while (p) { + *p = '_'; + p = strchr (p, '$'); } - sprintf (varname, "%s:n", fieldName); - set_int_variable (varname, Ny); - continue; + set_str_variable (varname, value); + free (value); } - // save char-type field as a string-type vector: - if (!CharAsVectors) { - Vector *myVector = NULL; - sprintf (name, "%s", argv[i]); - if ((myVector = SelectVector (name, ANYVECTOR, TRUE)) == NULL) ESCAPE ("bad vector name\n"); - ResetVector (myVector, OPIHI_STR, Ny); - - char *Ptr = data; - for (j = 0; j < Ny; j++) { - myVector[0].elements.Str[j] = strncreate (&Ptr[j*Nval], Nval); - // replace instances of $ with _ - char *p = strchr (myVector[0].elements.Str[j], '$'); - while (p) { *p = '_'; p = strchr (p, '$'); } - } - continue; - } + sprintf (varname, "%s:n", fieldName); + set_int_variable (varname, Ny); + continue; } // define the multifield vector names (Nval vectors x Ny elements) - // CharAsVectors is handled below automatically ALLOCATE (vec, Vector *, Nval); for (j = 0; j < Nval; j++) { Index: trunk/Ohana/src/opihi/cmd.data/roll.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/roll.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/roll.c (revision 41340) @@ -3,9 +3,7 @@ int roll (int argc, char **argv) { + int Nbytes, Nextra; + int dX, dx, dy, nx, ny; Buffer *buf; - - // this function is probably not needed - gprint (GP_ERR, "ERROR: use 'shift' instead of 'roll'\n"); - return (FALSE); if (argc != 4) { @@ -16,71 +14,32 @@ if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - int dx = atoi (argv[2]); - int dy = atoi (argv[3]); - - /* if (dx,dy < 0), we are moving the start position backwards, - if (dx,dy > 0), we are moving the start position forward */ - - int Nx = buf[0].matrix.Naxis[0]; - int Ny = buf[0].matrix.Naxis[1]; - - if (dy == 0) { - // moves in just dx can be done row-by-row - - int dX = abs(dx); - int Nbytes = (Nx * Ny - dX) * sizeof (float); - int Nextra = dX * sizeof (float); - - if (dx < 0) { - memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes); - memset (&buf[0].matrix.buffer[Nbytes], 0, Nextra); - } else { - memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes); - memset (buf[0].matrix.buffer, 0, Nextra); - } - - return TRUE; + dx = atof (argv[2]); + dy = atof (argv[3]); + if (dy != 0) { + gprint (GP_ERR, "only x rolls implemented for now\n"); } - if (dx == 0) { - // moves in just dy can be done row-by-row + /* if (dx < 0), we are moving the start position back by dx pixels, + if (dx > 0), we are moving the start position forward by dx pixels */ + + nx = buf[0].matrix.Naxis[0]; + ny = buf[0].matrix.Naxis[1]; - ALLOCATE_PTR (output, float, Nx*Ny); + dX = abs(dx); + Nbytes = nx * ny * sizeof (float); + Nextra = (nx * ny + dX) * sizeof (float); - float *input = (float *) buf[0].matrix.buffer; + REALLOCATE (buf[0].matrix.buffer, char, Nextra); - for (int iy = 0; iy < Ny; iy++) { - if (iy + dy < 0) continue; - if (iy + dy >= Ny) continue; - memcpy (&output[iy*Nx], &input[(iy + dy)*Nx], Nx); - } - if (dy > 0) { - memset (output, 0, dy*Nx*sizeof(float)); - } else { - int Nbytes = Nx * (Ny - abs(dy)) * sizeof (float); - int Nextra = abs(dy) * Nx * sizeof (float); - memset (&output[Nbytes], 0, Nextra); - } - free (buf[0].matrix.buffer); - buf[0].matrix.buffer = (char *) output; - return TRUE; + if (dx < 0) { + memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes); + } else { + memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes); } - ALLOCATE_PTR (output, float, Nx*Ny); - float *input = (float *) buf[0].matrix.buffer; + REALLOCATE (buf[0].matrix.buffer, char, Nbytes); - for (int iy = 0; iy < Ny; iy++) { - if (iy + dy < 0) continue; - if (iy + dy >= Ny) continue; - for (int ix = 0; ix < Nx; ix++) { - if (ix + dx < 0) continue; - if (ix + dx >= Nx) continue; - output[ix + iy*Nx] = input[(ix + dx) + (iy + dy)*Nx]; - } - } - free (buf[0].matrix.buffer); - buf[0].matrix.buffer = (char *) output; - return TRUE; + return (TRUE); + } - Index: trunk/Ohana/src/opihi/cmd.data/rotate.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/rotate.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/rotate.c (revision 41340) @@ -3,19 +3,19 @@ int rotate (int argc, char **argv) { - int i, j, NX, NY; + int i, j, NX, NY, X, Y, Lx, Ly, N; float *in_buff, *out_buff, *c; - double Xo, Yo, x, y, X1, Y1; + double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, x, y, X1, Y1; double pc11, pc12, pc21, pc22, PC11, PC12, PC21, PC22; Buffer *buf; -// Xo = 0; -// Yo = 0; -// if ((N = get_argument (argc, argv, "-center"))) { -// remove_argument (N, &argc, argv); -// Xo = atof(argv[N]); -// remove_argument (N, &argc, argv); -// Yo = atof(argv[N]); -// remove_argument (N, &argc, argv); -// } + Xo = 0; + Yo = 0; + if ((N = get_argument (argc, argv, "-center"))) { + remove_argument (N, &argc, argv); + Xo = atof(argv[N]); + remove_argument (N, &argc, argv); + Yo = atof(argv[N]); + remove_argument (N, &argc, argv); + } if (argc != 3) { @@ -221,26 +221,30 @@ } - double angle = atof (argv[2]); - double CosAngle = cos (angle*RAD_DEG); - double SinAngle = sin (angle*RAD_DEG); + angle = atof (argv[2]); + CosAngle = cos (angle*RAD_DEG); + SinAngle = sin (angle*RAD_DEG); gprint (GP_ERR, "rotating: %f %f %f\n", angle, CosAngle, SinAngle); - // we are rotating about the center pixel, (NX/2, NY/2), - // but we are then putting the result in a new image - // of size (Lx,Ly). - - // (x,y) = (i - Nx/2),(j - Ny/2) - // (x',y') = R(theta) (x,y) - // (I,J) = (x' + Lx/2),(y' + Ly/2) - - int Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle); - int Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle); + Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle); + Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle); + dX = MAX(0,NY*SinAngle); + dY = MAX(0,-NX*SinAngle); + /* + gprint (GP_ERR, "%f %f --> ", Xo, Yo); + X1 = Xo*CosAngle - Yo*SinAngle + dX; + Y1 = Xo*SinAngle + Yo*CosAngle + dY; + gprint (GP_ERR, "%f %f\n", X1, Y1); + */ /* fix reference pixel */ gfits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo); gfits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo); - X1 = (Xo - NX/2)*CosAngle + (Yo - NY/2)*SinAngle + Lx/2; - Y1 = (NX/2 - Xo)*SinAngle + (Yo - NY/2)*CosAngle + Ly/2; + /* + X1 = (Xo - dX)*CosAngle - (Yo - dY)*SinAngle; + Y1 = (dX - Xo)*SinAngle + (Yo - dY)*CosAngle; + */ + X1 = Xo*CosAngle - Yo*SinAngle + dX; + Y1 = Xo*SinAngle + Yo*CosAngle + dY; gfits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1); gfits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1); @@ -251,8 +255,8 @@ gfits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21); gfits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22); - PC11 = pc11*CosAngle + pc21*SinAngle; - PC12 = pc12*CosAngle + pc22*SinAngle; - PC21 = pc21*CosAngle - pc11*SinAngle; - PC22 = pc22*CosAngle - pc12*SinAngle; + PC11 = pc11*CosAngle - pc21*SinAngle; + PC12 = pc12*CosAngle - pc22*SinAngle; + PC21 = pc21*CosAngle + pc11*SinAngle; + PC22 = pc22*CosAngle + pc12*SinAngle; gfits_modify (&buf[0].header, "PC001001", "%le", 1, PC11); gfits_modify (&buf[0].header, "PC001002", "%le", 1, PC12); @@ -266,33 +270,22 @@ gfits_create_matrix (&buf[0].header, &buf[0].matrix); gfits_print_alt (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!"); - out_buff = (float *)buf[0].matrix.buffer; - for (j = 0; j < Ly; j++) { for (i = 0; i < Lx; i++, out_buff++) { - float xo = (i - Lx/2); - float yo = (j - Ly/2); - - x = xo*CosAngle + yo*SinAngle; - y = -xo*SinAngle + yo*CosAngle; - - int I = x + NX/2; - int J = y + NY/2; - - if (I < 0) continue; - if (I >= NX - 1) continue; - if (J < 0) continue; - if (J >= NY - 1) continue; - - c = &in_buff[I + NX*J]; - - int X = (int) x; - int Y = (int) y; - double fx = x - X; - double fy = y - Y; - + x = (i - dX)*CosAngle + (j - dY)*SinAngle; + y = (dX - i)*SinAngle + (j - dY)*CosAngle; + X = (int) x; + Y = (int) y; + + if (X < 0) continue; + if (X >= NX - 1) continue; + if (Y < 0) continue; + if (Y >= NY - 1) continue; + + c = &in_buff[X + NX*Y]; + fx = x - X; + fy = y - Y; *out_buff = (c[0]*(1-fx) + c[1]*fx)*(1-fy) + (c[NX+1]*fx + c[NX]*(1-fx))*fy; -// *out_buff = c[0]; } } Index: trunk/Ohana/src/opihi/cmd.data/spline.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/spline.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/spline.c (revision 41340) @@ -7,5 +7,4 @@ int spline_load (int argc, char **argv); int spline_save (int argc, char **argv); -int spline_print (int argc, char **argv); int spline_delete (int argc, char **argv); int spline_rename (int argc, char **argv); @@ -19,5 +18,4 @@ {1, "load", spline_load, "write a spline to a FITS file"}, {1, "save", spline_save, "read a spline from a FITS file"}, - {1, "print", spline_print, "print a spline to stdout"}, {1, "delete", spline_delete, "delete a spline"}, {1, "rename", spline_rename, "rename a spline"}, @@ -37,5 +35,4 @@ gprint (GP_ERR, " spline load (spline) (filename) : load a spline from a FITS file\n"); gprint (GP_ERR, " spline save (spline) (filename) [-append] : save a spline in FITS format\n"); - gprint (GP_ERR, " spline print (spline) : print a spline to stdout\n"); return FALSE; Index: trunk/Ohana/src/opihi/cmd.data/spline_commands.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/spline_commands.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/spline_commands.c (revision 41340) @@ -28,9 +28,4 @@ if (xvec->Nelements != yvec->Nelements) { gprint (GP_ERR, "x and y vectors must be the same length\n"); - return FALSE; - } - - if (xvec->Nelements < 3) { - gprint (GP_ERR, "cannot make a spline with fewer than 3 knots\n"); return FALSE; } @@ -98,26 +93,4 @@ } -int spline_print (int argc, char **argv) { - - int i; - - if (argc != 2) { - gprint (GP_ERR, "USAGE: spline print (name)\n"); - return FALSE; - } - - Spline *myspline = FindSpline (argv[1]); - if (!myspline) { - gprint (GP_ERR, "spline %s not found\n", argv[1]); - return (FALSE); - } - - for (i = 0; i < myspline->Nknots; i++) { - gprint (GP_LOG, "%lf : %lf : %lf\n", myspline->xk[i], myspline->yk[i], myspline->y2[i]); - } - - return TRUE; -} - int spline_load (int argc, char **argv) { @@ -136,11 +109,6 @@ int spline_delete (int argc, char **argv) { - int N, status; - - int QUIET = FALSE; - if ((N = get_argument (argc, argv, "-q"))) { - QUIET = TRUE; - remove_argument (N, &argc, argv); - } + int status; + Spline *spline; if (argc != 2) { @@ -149,7 +117,6 @@ } - Spline *spline = FindSpline (argv[1]); + spline = FindSpline (argv[1]); if (spline == NULL) { - if (QUIET) return TRUE; gprint (GP_ERR, "spline %s not found\n", argv[1]); return FALSE; @@ -163,4 +130,6 @@ int spline_rename (int argc, char **argv) { + Spline *spline; + if (argc != 3) { gprint (GP_ERR, "USAGE: spline rename (spline) (newname)\n"); @@ -168,13 +137,8 @@ } - Spline *spline = FindSpline (argv[1]); + spline = FindSpline (argv[1]); if (spline == NULL) { gprint (GP_ERR, "spline %s not found\n", argv[1]); return FALSE; - } - - Spline *newname = FindSpline (argv[2]); - if (newname != NULL) { - DeleteSpline (newname); } @@ -183,2 +147,24 @@ return TRUE; } + +/* +int spline_listspline (int argc, char **argv) { + + Spline *spline; + + if (argc != 2) { + gprint (GP_ERR, "USAGE: spline listspline (spline)\n"); + return FALSE; + } + + spline = FindSpline (argv[1]); + if (spline == NULL) { + gprint (GP_ERR, "spline %s not found\n", argv[1]); + return FALSE; + } + + ListPages (spline); + return TRUE; +} +*/ + Index: trunk/Ohana/src/opihi/cmd.data/test/cut.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/cut.sh (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/test/cut.sh (revision 41340) @@ -4,46 +4,4 @@ memtest1 end - -macro test2 - - mgaussdev z 50 100 0.0 1.0 - - cut z x y_sum x 0 0 z[][0] z[0][] - cut z x y_median x 0 0 z[][0] z[0][] -median - cut z x y_mean x 0 0 z[][0] z[0][] -mean - cut z x y_inner x 0 0 z[][0] z[0][] -inner - -end - -macro checkrange - if ($0 != 2) - echo "USAGE: checkrange (Npts)" - break - end - - $Npts = $1 - if ($Npts % 2) - $Ncenter = int($Npts / 2) - $Nquarter = int(0.25 * $Npts) - $Ns = $Ncenter - $Nquarter - $Ne = $Ncenter + $Nquarter - else - $Ncenter = int($Npts / 2) - 1 - $Nquarter = int(0.25 * $Npts) - $Ns = $Ncenter - $Nquarter - $Ne = $Ncenter + $Nquarter + 1 - end - if ($Ns < 0) - echo "error: Ns < 0: $Ns $Ncenter $Nquarter" - end - if ($Ne >= $Npts) - echo "error: Ne >= Npts: $Ne $Ncenter $Nquarter" - end - echo "$Npts : $Ncenter : $Nquarter : $Ns $Ne" - for i 0 $Npts - echo $i {($i >= $Ns) && ($i <= $Ne)} - end -end - # Test if cut works @@ -79,42 +37,2 @@ end end - -macro memtest1 - - $i = 0 - mcreate tim 100 100 - cut tim xdir imx X 40 4 60 6 - # do one to set up memory that should stay used - - memory check - for i 0 100 - cut tim xdir imx X 40 4 60 6 - end - memory check - - for i 0 100 - cut tim xdir imx y 40 4 60 6 - end - memory check - - for i 0 100 - cut -median tim xdir imx X 40 4 60 6 - end - memory check - - for i 0 100 - cut -median tim xdir imx y 40 4 60 6 - end - memory check - - for i 0 100 - cut -mean tim xdir imx X 40 4 60 6 - end - memory check - - for i 0 100 - cut -mean tim xdir imx y 40 4 60 6 - end - memory check - -end Index: trunk/Ohana/src/opihi/cmd.data/test/delete.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/delete.sh (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/test/delete.sh (revision 41340) @@ -49,13 +49,2 @@ end - -# Test delete -macro test2 - - $PASS = 1 - - for i 0 1000 - create v 0 200 - delete v - end -end Index: trunk/Ohana/src/opihi/cmd.data/test/headtest.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/headtest.sh (revision 41323) +++ (revision ) @@ -1,15 +1,0 @@ - -macro init.test - create x 0 10 - set y = x^2 - mcreate blank 0 0 - keyword blank TE0 -wd 56 - keyword blank TE1 -w "TEST LINE" - keyword blank TE2 -wf 5.12342341 - keyword blank TE3 -wb T - keyword blank TE4 -wc "not sure what do say" - keyword blank TE4 -wf 4.5 - keyword blank TE5 -ws "this is a long comment" - keyword blank TE2 -wc "this is a short comment" - write -fits DATA test2.fits x y -header blank -end Index: trunk/Ohana/src/opihi/cmd.data/test/medimage.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/medimage.sh (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/test/medimage.sh (revision 41340) @@ -1,540 +1,11 @@ -macro test.mean - - $Nsample = 16 - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 50 50 0.0 1.0 - medimage add t1 t - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 +macro go + mcreate a 30 30 + for i 0 40 + set a$i = zero(a) + $i end - medimage calc t1 T -mean - - imhist -q T x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo "expect {1/sqrt($Nsample)} : $C1" - plot -c red -x line x yf + for i 0 40 + medimage add t1 a$i + end end - -macro test.median - - $Nsample = 16 - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 50 50 0.0 1.0 - medimage add t1 t - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - end - - # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N) - # (somewhat higher scatter) - medimage calc t1 T - - imhist -q T x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo "expect {1/sqrt($Nsample)} : $C1 (actually should be a bit higher)" - plot -c red -x line x yf -end - -macro test.wtmean - - $Nsample = 8 - $sig1 = 1.0 - $sig2 = 3.0 - - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 50 50 0.0 $sig1 - set v = $sig1^2 + zero(t) - medimage add t1 t -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - for i 0 $Nsample - mgaussdev t 50 50 0.0 $sig2 - set v = $sig2^2 + zero(t) - medimage add t1 t -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N) - # (somewhat higher scatter) - medimage calc t1 T -wtmean - - imhist -q T x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2 - echo "expect {1/sqrt($S1)} : $C1" - plot -c red -x line x yf -end - -macro test.irls - medimage delete -q t1 - $Nsample = 16 - $sig = 1.0 - for i 0 $Nsample - mgaussdev t 50 50 0.0 $sig - set v = $sig^2 + zero(t) - - set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t) - set ts = t + bad - - medimage add t1 ts -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - # get stats for straight mean: - medimage calc t1 Tm -mean - - imhist -q Tm x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo "sigma from straight stdev: $C1" - # stats Tm - - plot -c red -x line x yf - - # get stats for irls - medimage calc t1 Ti -irls - - imhist -q Ti x y -range -10 10 -delta 0.1 - lim -n 2 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})" - # stats Ti - - plot -c red -x line x yf -end - - -###################33 - - -macro test.mean.var - - $Nsample = 64 - $sig = 2.0 - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 100 100 0.0 $sig - medimage add t1 t - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - medimage calc t1 T -mean -variance Tv - - imhist -q T x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - plot -c red -x line x yf - - imhist Tv xv yv -range -1 4 -delta 0.1 - lim -n 2 xv yv; clear; box; plot xv yv -x hist - - stat Tv - echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}" -end - -macro test.median.var - - $Nsample = 64 - $sig = 2.0 - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 50 50 0.0 $sig - medimage add t1 t - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N) - # (somewhat higher scatter) - medimage calc t1 T -variance Tv - - imhist -q T x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - plot -c red -x line x yf - - imhist Tv xv yv -range -1 4 -delta 0.1 - lim -n 2 xv yv; clear; box; plot xv yv -x hist - - stat Tv - echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}" -end - -macro test.wtmean.var - - $Nsample = 32 - $sig1 = 1.0 - $sig2 = 1.0 - - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 50 50 0.0 $sig1 - set v = $sig1^2 + zero(t) - medimage add t1 t -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - for i 0 $Nsample - mgaussdev t 50 50 0.0 $sig2 - set v = $sig2^2 + zero(t) - medimage add t1 t -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N) - # (somewhat higher scatter) - medimage calc t1 T -wtmean -variance Tv - - imhist -q T x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - plot -c red -x line x yf - - stat -q Tv - $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2 - echo $C1 vs {sqrt($MEDIAN)} : expect {1/sqrt($S1)} -end - -macro test.irls.var - - $Nsample = 16 - $sig = 1.0 - - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 50 50 0.0 $sig - set v = $sig^2 + zero(t) - - set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t) - set ts = t + bad - - medimage add t1 ts -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - # get stats for straight mean: - medimage calc t1 Tm -mean -variance Tv - - imhist -q Tm x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo "sigma from straight stdev: $C1" - # stats Tm - - plot -c red -x line x yf - - # get stats for irls - medimage calc t1 Ti -irls -variance Tv - - imhist -q Ti x y -range -10 10 -delta 0.1 - lim -n 2 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})" - # stats Ti - - plot -c red -x line x yf - - set dTv = sqrt(Tv) - imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist - - stat -q Tv - echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})" -end - -macro test.irls.boot.var - - $Nsample = 64 - $sig = 1.0 - - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 200 200 0.0 $sig - set v = $sig^2 + zero(t) - - set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t) - set ts = t + bad - - mgaussdev noise 200 200 0.0 0.5 - set ts = ts + noise - - medimage add t1 ts -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - # echo $C1 - end - - # get stats for straight mean: - medimage calc t1 Tm -mean -variance Tv - - imhist -q Tm x y -range -10 10 -delta 0.02 - lim -n 1 x y; clear; box; plot -x hist x y - peak -q x y - $C0 = $peakpos - $C1 = 1.5*$sig / sqrt($Nsample) - $C2 = $peakval - $C3 = 0 - vgauss -q x y con yf - echo "sigma from straight stdev: $C1" - # stats Tm - - plot -c red -x line x yf - - # get stats for irls - date - medimage calc t1 Ti -irls -variance Tv -bootstrap-iter 100 - date - - imhist -q Ti x y -range -10 10 -delta 0.02 - lim -n 2 x y; clear; box; plot -x hist x y - peak -q x y - $C0 = $peakpos - $C1 = 1.5*$sig / sqrt($Nsample) - $C2 = $peakval - $C3 = 0 - vgauss -q x y con yf - echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})" - # stats Ti - - plot -c red -x line x yf - - set dTv = sqrt(Tv) - 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 - - stat -q Tv - echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})" -end - -############################## -macro test.irls.boot.test - - $Nsample = 100 - $sig1 = 1.0 - - medimage delete -q t1 - for i 0 $Nsample - mgaussdev t 100 100 0.0 $sig1 - set v = $sig1^2 + zero(t) - - medimage add t1 t -variance v - end - - # get stats for irls - medimage calc t1 Ti -irls -variance Tv -bootstrap - - imhist -q Ti x y -range {-10*$sig1/sqrt($Nsample)} {10*$sig1/sqrt($Nsample)} -delta 0.01 - lim -n 2 x y; clear; box; plot -x hist x y - peak -q x y - $C0 = $peakpos - $C1 = 1.5*$sig1/sqrt($Nsample) - $C2 = $peakval - $C3 = 0 - vgauss x y con yf - echo "sigma from irls: $C1 (ideal is {$sig1/sqrt($Nsample)})" - # stats Ti - - plot -c red -x line x yf - - set dTv = sqrt(Tv) - 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 - - stat -q irls_npt - $Npix = $MEAN - - stat -q Tv - echo "sigma of irls average: $C1, sqrt(mean) of irls variance: {sqrt($MEAN)}, (ideal is {$sig1/sqrt($Npix)})" -end - -macro test.irls.range.var - medimage delete -q t1 - for i 0 8 - mgaussdev t 50 50 0.0 1.0 - set v = 1.0 + zero(t) - - set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t) - set ts = t + bad - - medimage add t1 ts -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo $C1 - end - - for i 0 8 - mgaussdev t 50 50 0.0 3.0 - set v = 3.0 + zero(t) - - set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t) - set ts = t + bad - - medimage add t1 ts -variance v - imhist -q t x y -range -10 10 -delta 0.1 - - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo $C1 - end - - # get stats for straight mean: - medimage calc t1 Tm -mean -variance Tv - - imhist -q Tm x y -range -10 10 -delta 0.1 - lim -n 1 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo $C1 - stats Tm - - plot -c red -x line x yf - - # get stats for irls - medimage calc t1 Ti -irls -variance Tv - - imhist -q Ti x y -range -10 10 -delta 0.1 - lim -n 2 x y; clear; box; plot -x hist x y - $C0 = 0 - $C1 = 1.5 - $C2 = 400 - $C3 = 0 - vgauss -q x y con yf - echo $C1 - stats Ti - - plot -c red -x line x yf - - stat -q Tv - echo $C1 vs {sqrt($MEDIAN)} - - set dTv = sqrt(Tv) - imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist -end Index: trunk/Ohana/src/opihi/cmd.data/test/set.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/set.sh (revision 41323) +++ (revision ) @@ -1,43 +1,0 @@ - -macro memtest1 - - create x -10 10 0.1 - set y = zero(x) - $i = 0 - - memory check - for i 0 1000 - set y = exp(x) - end - memory check - -end - -macro memtest2 - - create x -10 10 0.1 - set y = zero(x) - $i = 0 - - memory check - for i 0 1000 - set y = (x < 0) ? x : exp(x) - end - memory check - -end - -macro memtest3 - - create x -10 10 0.1 - set y = zero(x) - $i = 0 - - memory check - for i 0 1000 - set y = (x < 5 + 2*$i) ? x : exp(x) - end - memory check - -end - Index: trunk/Ohana/src/opihi/cmd.data/test/spline.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/spline.sh (revision 41323) +++ (revision ) @@ -1,11 +1,0 @@ - -macro test1 - - for i 0 100 - create x 0 50 - set y = x^2 - spline create t$i x y - spline rename t$i test - spline delete test - end -end Index: trunk/Ohana/src/opihi/cmd.data/test/vsigmoid.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/vsigmoid.sh (revision 41323) +++ (revision ) @@ -1,263 +1,0 @@ - -list tests - test1 - test2 - test3 -end - -macro test1 - $PASS = 1 - break -auto off - - create x -10 10 0.1 - set y = 0 + 5 / (1 + exp((x - 3)/4)) - - $C0 = 1 - $C1 = 2 - $C2 = 10 - $C3 = 1 - set dy = sqrt(y) - - vsigmoid -q x y dy yfit - lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit - # cursor - - if (abs($C0 - 0.0) > 0.01) - $PASS = 0 - end - if (abs($C1 - 3.0) > 0.01) - $PASS = 0 - end - if (abs($C2 - 5.0) > 0.01) - $PASS = 0 - end - if (abs($C3 - 0.0) > 0.01) - $PASS = 0 - end -end - -# noise of 0.1 -macro test2 - $PASS = 1 - break -auto off - - create x -10 10 0.1 - set y = 50 + 1000 / (1 + exp((x - 3)/2)) - set dy = (rnd(y) - 0.5)/0.5 - set y = y + dy - - $C0 = 2 - $C1 = 3 - $C2 = 900 - $C3 = 1 - - vsigmoid -q x y dy yfit - lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit - # cursor - - if (abs($C0 - 0.0) > 0.01) - $PASS = 0 - end - if (abs($C1 - 3.0) > 0.01) - $PASS = 0 - end - if (abs($C2 - 1000.0) > 0.1) - $PASS = 0 - end - if (abs($C3 - 0.0) > 0.1) - $PASS = 0 - end -end - -# poisson-distributed noise -macro test3 - $PASS = 1 - break -auto off - - $C0o = 3 - $C1o = 2 - $C2o = 20 - $C3o = -2 - - create x -10 10 0.1 - set y = $C3o + $C2o / (1 + exp((x - $C0o)/$C1o)) - - gaussdev dY y[] 0.0 1.0 - set dy = dY * sqrt(abs(y)) - set y = y + dy - - $C0 = $C0o + 2 - $C1 = $C1o + 2 - $C2 = $C2o + 50 - $C3 = $C3o + 10 - - vsigmoid x y dy yfit - lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit - # cursor - - $dS = 3.0 / sqrt(y[]) - - if (abs($C0 - $C0o) > $dS) - $PASS = 0 - end - if (abs($C1 - $C1o)/$C1o > $dS) - $PASS = 0 - end - if (abs($C2 - $C2o)/$C2o > $dS) - $PASS = 0 - end - if (abs($C3 - $C3o) > $dS) - $PASS = 0 - end -end - -macro test1rev - $PASS = 1 - break -auto off - - create x -10 10 0.1 - set y = 5 / (1 + exp(-(x - 3)/4)) - - $C0 = 1 - $C1 = 2 - $C2 = 10 - $C3 = 1 - set dy = sqrt(abs(y)) - - vsigmoid -r -q x y dy yfit - lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit - # cursor - - if (abs($C0 - 0.0) > 0.01) - $PASS = 0 - end - if (abs($C1 - 3.0) > 0.01) - $PASS = 0 - end - if (abs($C2 - 5.0) > 0.01) - $PASS = 0 - end - if (abs($C3 - 0.0) > 0.01) - $PASS = 0 - end -end - -# noise of 0.1 -macro test2rev - $PASS = 1 - break -auto off - - create x -10 10 0.1 - set y = 50 + 1000 / (1 + exp(-(x - 3)/2)) - set dy = (rnd(y) - 0.5)/0.5 - set y = y + dy - - $C0 = 2 - $C1 = 3 - $C2 = 900 - $C3 = 1 - - vsigmoid -r -q x y dy yfit - lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit - # cursor - - if (abs($C0 - 0.0) > 0.01) - $PASS = 0 - end - if (abs($C1 - 3.0) > 0.01) - $PASS = 0 - end - if (abs($C2 - 1000.0) > 0.1) - $PASS = 0 - end - if (abs($C3 - 0.0) > 0.1) - $PASS = 0 - end -end - -# poisson-distributed noise -macro test3rev - $PASS = 1 - break -auto off - - $C0o = 3 - $C1o = 2 - $C2o = 20 - $C3o = -2 - - create x -10 10 0.1 - set y = $C3o + $C2o / (1 + exp(-(x - $C0o)/$C1o)) - - gaussdev dY y[] 0.0 1.0 - set dy = dY * sqrt(abs(y)) - set y = y + dy - - $C0 = $C0o + 2 - $C1 = $C1o + 2 - $C2 = $C2o + 50 - $C3 = $C3o + 10 - - vsigmoid -r -q x y dy yfit - lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit - # cursor - - $dS = 3.0 / sqrt(y[]) - - if (abs($C0 - $C0o) > $dS) - $PASS = 0 - end - if (abs($C1 - $C1o)/$C1o > $dS) - $PASS = 0 - end - if (abs($C2 - $C2o)/$C2o > $dS) - $PASS = 0 - end - if (abs($C3 - $C3o) > $dS) - $PASS = 0 - end -end - -macro memtest1 - - create x -10 10 0.1 - set y = 0 + 5 / (1 + exp((x - 3)/4)) - - $C0 = 1 - $C1 = 2 - $C2 = 10 - $C3 = 1 - set dy = sqrt(y) - set yfit = zero(y) - - $i = 0 - - memory check - for i 0 1000 - vsigmoid -q x y dy yfit - end - memory check - -end - -macro memtest2 - - create x -10 10 0.1 - set y = 0 + 5 / (1 + exp((x - 3)/4)) - - $C0 = 1 - $C1 = 2 - $C2 = 10 - $C3 = 1 - set dy = sqrt(y) - set yfit = zero(y) - - $i = 0 - - memory check - for i 0 1000 - vsigmoid -q -r x y dy yfit - end - memory check - -end - Index: trunk/Ohana/src/opihi/cmd.data/test/vstat.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/vstat.sh (revision 41323) +++ (revision ) @@ -1,28 +1,0 @@ - -macro memtest1 - - # preset variables so memory usage does not change - $PMIN = -10 - $PMAX = 9.9 - $MIN = -10 - $MAX = 9.9 - $MEDIAN = -0.1097 - $MEAN = -0.05 - $MODE = -10 - $TOTAL = -10 - $NPIX = 200 - $NPTS = 200 - $NUSED = 200 - $SIGMA = 5.7879184514 - - create x -10 10 0.1 - $i = 0 - - memory check - for i 0 1000 - vstat -q x - end - memory check - -end - Index: trunk/Ohana/src/opihi/cmd.data/threshold.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/threshold.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/threshold.c (revision 41340) @@ -25,14 +25,4 @@ remove_argument (N, &argc, argv); BinMax = atoi (argv[N]); - remove_argument (N, &argc, argv); - } - - float ValMin = NAN; - float ValMax = NAN; - if ((N = get_argument (argc, argv, "-vrange"))) { - remove_argument (N, &argc, argv); - ValMin = atof (argv[N]); - remove_argument (N, &argc, argv); - ValMax = atof (argv[N]); remove_argument (N, &argc, argv); } @@ -74,16 +64,4 @@ } - // ValMin, ValMax not compatible with BinMin, BinMax - if (isfinite(ValMin)) { - int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMin); - BinMin = (vecx[0].elements.Flt[binTest+1] == ValMin) ? binTest : binTest + 1; - // XXX dangerous: check that this does not run off the end, etc - } - if (isfinite(ValMax)) { - int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMax); - BinMax = (vecx[0].elements.Flt[binTest+1] == ValMax) ? binTest : binTest + 1; - // XXX dangerous: check that this does not run off the end, etc - } - // use bisection to find the value @@ -93,13 +71,4 @@ if (!REVERSE) { - // this algorithm assumes vecy[BinMax] > threshold, vecy[BinMin] < threshold - if (vecy[0].elements.Flt[BinMin] > value) { - gprint (GP_ERR, "ERROR: all values above threshold\n"); - return FALSE; - } - if (vecy[0].elements.Flt[BinMax] < value) { - gprint (GP_ERR, "ERROR: all values below threshold\n"); - return FALSE; - } while (Nhi - Nlo > 10) { N = 0.5*(Nlo + Nhi); @@ -120,13 +89,4 @@ } } else { - // this algorithm assumes vecy[BinMin] > threshold, vecy[BinMax] < threshold - if (vecy[0].elements.Flt[BinMin] < value) { - gprint (GP_ERR, "ERROR: all values below threshold\n"); - return FALSE; - } - if (vecy[0].elements.Flt[BinMax] > value) { - gprint (GP_ERR, "ERROR: all values above threshold\n"); - return FALSE; - } while (Nhi - Nlo > 10) { N = 0.5*(Nlo + Nhi); Index: trunk/Ohana/src/opihi/cmd.data/tvcolors.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/tvcolors.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/tvcolors.c (revision 41340) @@ -37,8 +37,21 @@ gprint (GP_ERR, "USAGE: tvcolors (colormap) [-nan red green blue]\n"); gprint (GP_ERR, " colormap options : greyscale, -greyscale, rainbow, heat, fullcolor, ruffcolor (also grayscale, -grayscale)\n"); - gprint (GP_ERR, " colormap file: if colormap name is given as file:path/to/file, the colormap is read from the file\n"); + gprint (GP_ERR, " colormap file name of the for WORD:PATH where WORD = (file,csvf,lgcy,cetf)\n"); + gprint (GP_ERR, " lgcy:\n"); + gprint (GP_ERR, " the colormap file must contain 4 columns: f R B G; each line defines a color transition.\n"); + gprint (GP_ERR, " f: 0 - 1 defines the scale value for the transition\n"); + gprint (GP_ERR, " R,B,G: 0 - 1 define the value of the color at the transition point\n"); + gprint (GP_ERR, " file:\n"); gprint (GP_ERR, " the colormap file must contain 4 columns: f R G B; each line defines a color transition.\n"); gprint (GP_ERR, " f: 0 - 1 defines the scale value for the transition\n"); gprint (GP_ERR, " R,G,B: 0 - 1 define the value of the color at the transition point\n"); + gprint (GP_ERR, " csvf:\n"); + gprint (GP_ERR, " the colormap file must contain 4 columns: f,R,G,B separate by commas; each line defines a color transition.\n"); + gprint (GP_ERR, " f: 0 - 1 defines the scale value for the transition\n"); + gprint (GP_ERR, " R,G,B: 0 - 1 define the value of the color at the transition point\n"); + gprint (GP_ERR, " cetf:\n"); + 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"); + gprint (GP_ERR, " R,G,B: 0 - 1 define the value of the color at the transition point\n"); + gprint (GP_ERR, " CETF-format files can be found at : https://peterkovesi.com/projects/colourmaps\n"); return (FALSE); } Index: trunk/Ohana/src/opihi/cmd.data/vgauss.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vgauss.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/vgauss.c (revision 41340) @@ -3,5 +3,4 @@ /* local private functions */ opihi_flt fgaussOD (opihi_flt, opihi_flt *, int, opihi_flt *); -int vgauss_apply (int argc, char **argv); # define GET_VAR(V,A) \ @@ -21,10 +20,4 @@ Vector *xvec, *yvec, *svec, *ovec; char *c, name[16]; - - if ((N = get_argument (argc, argv, "-apply"))) { - remove_argument (N, &argc, argv); - int status = vgauss_apply (argc, argv); - return status; - } Quiet = FALSE; @@ -111,6 +104,4 @@ } if (!Quiet) gprint (GP_ERR, "%d iterations\n", i); - set_variable ("Chisq", chisq); - set_variable ("ChisqNu", chisq / (float) (Npts - Npar)); for (i = 0; i < Npts; i++) { @@ -134,38 +125,4 @@ } -int vgauss_apply (int argc, char **argv) { - - opihi_flt par[4]; - Vector *xvec, *ovec; - char *c; - - // NOTE: -apply has already been stripped by vgauss call - if (argc != 3) { - gprint (GP_ERR, "USAGE: vgauss -apply \n"); - return (FALSE); - } - - if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((ovec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); - - CastVector (xvec, OPIHI_FLT); - - GET_VAR (par[0], "C0"); - GET_VAR (par[1], "C1"); - GET_VAR (par[2], "C2"); - GET_VAR (par[3], "C3"); - int Npar = 4; - - int Npts = xvec[0].Nelements; - ResetVector (ovec, OPIHI_FLT, Npts); - - for (int i = 0; i < Npts; i++) { - ovec[0].elements.Flt[i] = fgaussOD (xvec[0].elements.Flt[i], par, Npar, NULL); - } - ovec[0].Nelements = Npts; - - return TRUE; -} - /* pars: x_o, sigma, I, back */ opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) { @@ -178,10 +135,8 @@ f = par[2]*r + par[3]; - if (dpar) { - dpar[0] = par[2]*r*z/par[1]; - dpar[1] = par[2]*r*z*z/par[1]; - dpar[2] = r; - dpar[3] = 1; - } + dpar[0] = par[2]*r*z/par[1]; + dpar[1] = par[2]*r*z*z/par[1]; + dpar[2] = r; + dpar[3] = 1; return (f); Index: trunk/Ohana/src/opihi/cmd.data/virls.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/virls.c (revision 41323) +++ (revision ) @@ -1,103 +1,0 @@ -# include "data.h" - -double opihi_irls_mean_dbl (double *val, double *wgt, int N, int *Converged); -double weight_cauchy_square_dbl (double x2); - -# define IRLS_TOLERANCE 1e-4 - -int virls (int argc, char **argv) { - - if (argc != 3) { - gprint (GP_ERR, "USAGE: virls (values) (sigmas)\n"); - gprint (GP_ERR, "NOTE: sigmas with abs values < %e will be truncated to that value\n", FLT_MIN); - return (FALSE); - } - - Vector *value, *sigma; - - if ((value = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((sigma = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); - - if (value->Nelements != sigma->Nelements) { - gprint (GP_ERR, "vectors are not the same length\n"); - return FALSE; - } - if (value->type != OPIHI_FLT) { - gprint (GP_ERR, "%s is not a floating point vector\n", argv[1]); - return FALSE; - } - if (sigma->type != OPIHI_FLT) { - gprint (GP_ERR, "%s is not a floating point vector\n", argv[2]); - return FALSE; - } - - ALLOCATE_PTR (val, opihi_flt, value->Nelements); - ALLOCATE_PTR (wgt, opihi_flt, value->Nelements); - - int N = 0; - for (int i = 0; i < value->Nelements; i++) { - if (!isfinite(value->elements.Flt[i])) continue; - - double s = sigma->elements.Flt[i]; - if (!isfinite(s)) continue; - - if (fabs(s) < FLT_MIN) s = FLT_MIN; - - val[N] = value->elements.Flt[i]; - wgt[N] = 1.0 / SQ(s); - N ++; - } - - int converged = TRUE; - double Value = opihi_irls_mean_dbl (val, wgt, N, &converged); - - set_variable ("WTMEAN", Value); - set_int_variable ("NUSED", N); - set_int_variable ("CONVERGED", converged); - return (TRUE); -} - -double opihi_irls_mean_dbl (double *val, double *wgt, int N, int *Converged) { - - // calculate weighted mean - double S1 = 0.0, S2 = 0.0; - for (int n = 0; n < N; n++) { - S1 += wgt[n] * val[n]; - S2 += wgt[n]; - } - double Value = S1 / S2; - - int converged = FALSE; - for (int i = 0; (i < 10) && !converged; i++) { - double ValueLast = Value; - - double S1 = 0.0, S2 = 0.0; - - // calculate weight modification based on distances (squared). - // use modifier to calculate new weighted mean - for (int n = 0; n < N; n++) { - double dV = (val[n] - Value); - double d2 = SQ(dV) * wgt[n]; - - double Mod = weight_cauchy_square_dbl (d2); - S1 += Mod * wgt[n] * val[n]; - S2 += Mod * wgt[n]; - fprintf (stderr, "%f %f : %f %f : %f\n", val[n], wgt[n], dV, sqrt(d2), Mod); - } - Value = S1 / S2; - - double delta = fabs(Value - ValueLast); - if (delta < Value * IRLS_TOLERANCE) converged = TRUE; - - // XXX if the answer is close to zero, we might not converge - } - *Converged = converged; - return Value; -} - -# define CAUCY_FACTOR 1.0 - -double weight_cauchy_square_dbl (double x2) { - double r2 = x2 / CAUCY_FACTOR; - return (1.0 / (1.0 + r2)); -} Index: trunk/Ohana/src/opihi/cmd.data/vsigmoid.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vsigmoid.c (revision 41323) +++ (revision ) @@ -1,192 +1,0 @@ -# include "data.h" - -typedef opihi_flt FitFunc (opihi_flt, opihi_flt *, int, opihi_flt *); - -/* local private functions */ -opihi_flt fsigmoidOD (opihi_flt, opihi_flt *, int, opihi_flt *); -opihi_flt rsigmoidOD (opihi_flt, opihi_flt *, int, opihi_flt *); - -# define GET_VAR(V,A) { \ - char *c = get_variable (A); \ - if (c == NULL) { \ - gprint (GP_ERR, "missing fit parameter A\n"); \ - return (FALSE); \ - } \ - V = atof (c); \ - free (c); } - -// XXX this fitting function works but fails if the input vectors have NANs -int vsigmoid (int argc, char **argv) { - - int N; - Vector *xvec, *yvec, *svec, *ovec; - - FitFunc *FUNC = fsigmoidOD; - if ((N = get_argument (argc, argv, "-r"))) { - FUNC = rsigmoidOD; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-rev"))) { - FUNC = rsigmoidOD; - remove_argument (N, &argc, argv); - } - - int Quiet = FALSE; - if ((N = get_argument (argc, argv, "-q"))) { - Quiet = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-quiet"))) { - Quiet = TRUE; - remove_argument (N, &argc, argv); - } - - if (argc != 5) { - gprint (GP_ERR, "USAGE: vsigmoid [-q] [-r] [-rev] (out)\n"); - gprint (GP_ERR, " may be the words 'con[stant]' or 'poi[sson]', in which case the error is constructed'\n"); - gprint (GP_ERR, " uses guesses: C0 (center), C1 (sigma), C2 (top), C3 (bottom)\n"); - gprint (GP_ERR, " forward sigmoid : bottom + top / (1 + exp(-z)), z = (x - center) / sigma\n"); - gprint (GP_ERR, " -r / -rev : reverse sigmoid : bottom + top / (1 + exp(z))\n"); - return (FALSE); - } - - if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((ovec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE); - - CastVector (xvec, OPIHI_FLT); - CastVector (yvec, OPIHI_FLT); - - int Nsvec = strlen(argv[3]); - - if ((Nsvec >= 3) && (!strncasecmp ("poisson", argv[3], Nsvec) || !strncasecmp ("constant", argv[3], Nsvec))) { - svec = SelectVector ("vsigmoid_err", ANYVECTOR, TRUE); - if (svec == NULL) return (FALSE); - MatchVector (svec, xvec, OPIHI_FLT); - opihi_flt *v1 = svec[0].elements.Flt; - opihi_flt *v2 = yvec[0].elements.Flt; - if (!strncasecmp ("poisson", argv[3], Nsvec)) { - for (int i = 0; i < svec[0].Nelements; i++) { - v1[i] = (isfinite(v2[i]) && (v2[i] > 0)) ? sqrt(v2[i]) : 1.0; - } - } else { - for (int i = 0; i < svec[0].Nelements; i++) { - v1[i] = 1.0; - } - } - } else { - if ((svec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); - } - - CastVector (svec, OPIHI_FLT); - // XXX Cast is failing. - - int Npts = xvec[0].Nelements; - ResetVector (ovec, OPIHI_FLT, Npts); - - ALLOCATE_PTR (dy, opihi_flt, Npts); - - opihi_flt par[4]; - GET_VAR (par[0], "C0"); - GET_VAR (par[1], "C1"); - GET_VAR (par[2], "C2"); - GET_VAR (par[3], "C3"); - int Npar = 4; - - // mrqmin takes the inverse variance (do not generate NANs) - opihi_flt *v1 = svec[0].elements.Flt; - opihi_flt *v2 = dy; - for (int i = 0; i < Npts; i++, v1++, v2++) { - *v2 = (*v1 == 0.0) ? 0.0 : 1.0 / (*v1 * *v1); - } - - opihi_flt ochisq = mrqinit (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, FUNC, !Quiet); - opihi_flt dchisq = ochisq + 2*Npts; - - int Niter; - for (Niter = 0; (Niter < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); Niter++) { - opihi_flt chisq = mrqmin (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, FUNC, !Quiet); - dchisq = ochisq - chisq; - ochisq = chisq; - if (!Quiet) gprint (GP_ERR, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar); - } - if (!Quiet) gprint (GP_ERR, "%d iterations\n", Niter); - - for (int i = 0; i < Npts; i++) { - ovec[0].elements.Flt[i] = FUNC (xvec[0].elements.Flt[i], par, Npar, dy); - } - ovec[0].Nelements = Npts; - /* set output *before* variable renomalization */ - - opihi_flt **covar = mrqcovar (Npar); - for (int i = 0; i < Npar; i++) { - char name[16]; - sprintf (name, "C%d", i); - set_variable (name, par[i]); - if (!Quiet) gprint (GP_ERR, "%d %f %f\n", i, par[i], sqrt(covar[i][i])); - sprintf (name, "dC%d", i); - set_variable (name, sqrt(covar[i][i])); - } - - free (dy); - mrqfree (Npar); - return (TRUE); -} - -/* pars: x_o, sigma, top, bottom */ -opihi_flt rsigmoidOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) { - OHANA_UNUSED_PARAM(Npar); - - double z = (x - par[0])/par[1]; - double r = exp (z); - double s = 1 + r; - double q = 1 / s; - double f = par[2]*q + par[3]; - - // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0 - // df/dq = par[2] - // dq/ds = -s^-2 = -q/s - // ds/dr = 1 - // dr/dz = r - // dz/dp0 = -1 / par[1] - // dz/dp1 = -z / par[1] - - // df/dp0 = par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1]) - // df/dp1 = par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1]) - - dpar[0] = par[2]*q*r/(s*par[1]); - dpar[1] = dpar[0]*z; - dpar[2] = q; - dpar[3] = 1; - - return (f); -} - -/* pars: x_o, sigma, top, bottom */ -opihi_flt fsigmoidOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) { - OHANA_UNUSED_PARAM(Npar); - - double z = (x - par[0])/par[1]; - double r = exp (-z); - double s = 1 + r; - double q = 1 / s; - double f = par[2]*q + par[3]; - - // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0 - // df/dq = par[2] - // dq/ds = -s^-2 = -q/s - // ds/dr = 1 - // dr/dz = -r - // dz/dp0 = -1 / par[1] - // dz/dp1 = -z / par[1] - - // df/dp0 = -par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1]) - // df/dp1 = -par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1]) - - dpar[0] = -par[2]*q*r/(s*par[1]); - dpar[1] = dpar[0]*z; - dpar[2] = q; - dpar[3] = 1; - - return (f); -} Index: trunk/Ohana/src/opihi/cmd.data/vtransitions.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vtransitions.c (revision 41323) +++ (revision ) @@ -1,175 +1,0 @@ -# include "data.h" - -static int isFltX = FALSE; -static int isFltY = FALSE; -static double interpolate_position (Vector *vecx, Vector *vecy, int bin, float value); - -int vtransitions (int argc, char **argv) { - - int N; - Vector *vecx, *vecy, *outbin, *outval, *outpos; - - int QUIET = FALSE; - if ((N = get_argument (argc, argv, "-q"))) { - QUIET = TRUE; - remove_argument (N, &argc, argv); - } - - int UP_TRANSITIONS = TRUE; - if ((N = get_argument (argc, argv, "-dn"))) { - UP_TRANSITIONS = FALSE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-down"))) { - UP_TRANSITIONS = FALSE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-up"))) { - UP_TRANSITIONS = TRUE; - remove_argument (N, &argc, argv); - } - - int BinMin = -1; - int BinMax = -1; - if ((N = get_argument (argc, argv, "-range"))) { - remove_argument (N, &argc, argv); - BinMin = atoi (argv[N]); - remove_argument (N, &argc, argv); - BinMax = atoi (argv[N]); - remove_argument (N, &argc, argv); - } - - if (argc != 4) { - gprint (GP_ERR, "USAGE: vtransitions (value) (-up|-down,-dn) [-q] [-range BinMin BinMax]\n"); - gprint (GP_ERR, " find the x coordinates at which we pass the specified value going up\n"); - gprint (GP_ERR, " -q : quiet mode\n"); - gprint (GP_ERR, " -r : reverse (find downward transition)\n"); - return (FALSE); - } - - if ((vecx = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((vecy = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if (vecx[0].Nelements != vecy[0].Nelements) return FALSE; - int Nelements = vecx[0].Nelements; - - if ((outbin = SelectVector ("vtr_bin", ANYVECTOR, TRUE)) == NULL) return (FALSE); // bins containing transitions - if ((outval = SelectVector ("vtr_val", ANYVECTOR, TRUE)) == NULL) return (FALSE); // value of bins containing transitions - if ((outpos = SelectVector ("vtr_pos", ANYVECTOR, TRUE)) == NULL) return (FALSE); // interpolated coordinate of transitions - - double value = atof (argv[3]); - - isFltX = (vecx[0].type == OPIHI_FLT); - isFltY = (vecy[0].type == OPIHI_FLT); - - // BinMin = -1 -> range not provide - if (BinMin == -1) { - BinMin = 0; - BinMax = Nelements - 1; - } else { - if ((BinMin >= vecx[0].Nelements) || (BinMin < -1*vecx[0].Nelements)) { - gprint (GP_ERR, "vector subscript out of range\n"); - return FALSE; - } - if (BinMin < 0) BinMin += vecx[0].Nelements; - if ((BinMax >= vecx[0].Nelements) || (BinMax < -1*vecx[0].Nelements)) { - gprint (GP_ERR, "vector subscript out of range\n"); - return FALSE; - } - if (BinMax < 0) BinMax += vecx[0].Nelements; - } - - int Nout = 0; - int NOUT = 100; - ResetVector (outbin, OPIHI_INT, NOUT); - ResetVector (outval, vecy[0].type, NOUT); - ResetVector (outpos, OPIHI_FLT, NOUT); - - // starting disposition - int isUp = isFltY ? (vecy[0].elements.Flt[BinMin] > value) : (vecy[0].elements.Int[BinMin] > value); - - for (int i = BinMin; i <= BinMax; i++) { - double testval = isFltY ? vecy[0].elements.Flt[i] : vecy[0].elements.Int[i]; - if (isUp && (testval < value)) { - isUp = FALSE; - if (!UP_TRANSITIONS) { - // save this value - outbin->elements.Int[Nout] = i; - if (isFltY) { - outval->elements.Flt[Nout] = testval; - } else { - outval->elements.Int[Nout] = testval; - } - outpos->elements.Flt[Nout] = interpolate_position (vecx, vecy, i, value); - Nout ++; - if (Nout == NOUT) { - NOUT += 100; - ResetVector (outbin, OPIHI_INT, NOUT); - ResetVector (outval, vecy[0].type, NOUT); - ResetVector (outpos, OPIHI_FLT, NOUT); - } - } - continue; - } - if (!isUp && (testval > value)) { - isUp = TRUE; - if (UP_TRANSITIONS) { - // save this value - outbin->elements.Int[Nout] = i; - if (isFltY) { - outval->elements.Flt[Nout] = testval; - } else { - outval->elements.Int[Nout] = testval; - } - outpos->elements.Flt[Nout] = interpolate_position (vecx, vecy, i, value); - Nout ++; - if (Nout == NOUT) { - NOUT += 100; - ResetVector (outbin, OPIHI_INT, NOUT); - ResetVector (outval, vecy[0].type, NOUT); - ResetVector (outpos, OPIHI_FLT, NOUT); - } - } - } - } - - ResetVector (outbin, OPIHI_INT, Nout); - ResetVector (outval, vecy[0].type, Nout); - ResetVector (outpos, OPIHI_FLT, Nout); - - if (!QUIET) gprint (GP_LOG, "found %d transitions\n", Nout); - - return (TRUE); -} - -static double interpolate_position (Vector *vecx, Vector *vecy, int bin, float value) { - - double x0, x1, y0, y1, Xvalue; - - if (vecy[0].Nelements == 0) return NAN; - if (vecy[0].Nelements == 1) return isFltY ? vecy[0].elements.Flt[0] : vecy[0].elements.Int[0]; - - if (bin == 0) { - // interpolate to value: - y0 = isFltY ? vecy[0].elements.Flt[bin] : vecy[0].elements.Int[bin]; - y1 = isFltY ? vecy[0].elements.Flt[bin+1] : vecy[0].elements.Int[bin+1]; - x0 = isFltX ? vecx[0].elements.Flt[bin] : vecy[0].elements.Int[bin]; - x1 = isFltX ? vecx[0].elements.Flt[bin+1] : vecy[0].elements.Int[bin+1]; - if (y0 == y1) { - Xvalue = 0.5*(x0 + x1); - } else { - Xvalue = (value - y0) * (x1 - x0) / (y1 - y0) + x0; - } - } else { - // interpolate to value: - y0 = isFltY ? vecy[0].elements.Flt[bin-1] : vecy[0].elements.Int[bin-1]; - y1 = isFltY ? vecy[0].elements.Flt[bin] : vecy[0].elements.Int[bin]; - x0 = isFltX ? vecx[0].elements.Flt[bin-1] : vecy[0].elements.Int[bin-1]; - x1 = isFltX ? vecx[0].elements.Flt[bin] : vecy[0].elements.Int[bin]; - if (y0 == y1) { - Xvalue = 0.5*(x0 + x1); - } else { - Xvalue = (value - y0) * (x1 - x0) / (y1 - y0) + x0; - } - } - return Xvalue; -} Index: trunk/Ohana/src/opihi/cmd.data/vwtmean.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vwtmean.c (revision 41323) +++ (revision ) @@ -1,70 +1,0 @@ -# include "data.h" - -double opihi_wt_mean_dbl (double *val, double *wgt, int N, double *Sigma); - -int vwtmean (int argc, char **argv) { - - if (argc != 3) { - gprint (GP_ERR, "USAGE: vwtmean (values) (sigmas)\n"); - gprint (GP_ERR, "NOTE: sigmas with abs values < %e will be truncated to that value\n", FLT_MIN); - return (FALSE); - } - - Vector *value, *sigma; - - if ((value = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((sigma = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); - - if (value->Nelements != sigma->Nelements) { - gprint (GP_ERR, "vectors are not the same length\n"); - return FALSE; - } - if (value->type != OPIHI_FLT) { - gprint (GP_ERR, "%s is not a floating point vector\n", argv[1]); - return FALSE; - } - if (sigma->type != OPIHI_FLT) { - gprint (GP_ERR, "%s is not a floating point vector\n", argv[2]); - return FALSE; - } - - ALLOCATE_PTR (val, opihi_flt, value->Nelements); - ALLOCATE_PTR (wgt, opihi_flt, value->Nelements); - - int N = 0; - for (int i = 0; i < value->Nelements; i++) { - if (!isfinite(value->elements.Flt[i])) continue; - - double s = sigma->elements.Flt[i]; - if (!isfinite(s)) continue; - - if (fabs(s) < FLT_MIN) s = FLT_MIN; - - val[N] = value->elements.Flt[i]; - wgt[N] = 1.0 / SQ(s); - N ++; - } - - double Sigma = NAN; - double Value = opihi_wt_mean_dbl (val, wgt, N, &Sigma); - - set_variable ("WTMEAN", Value); - set_variable ("WTSIGMA", Sigma); - set_int_variable ("NUSED", N); - return (TRUE); -} - -double opihi_wt_mean_dbl (double *val, double *wgt, int N, double *Sigma) { - - // calculate weighted mean - double S1 = 0.0, S2 = 0.0; - for (int n = 0; n < N; n++) { - S1 += wgt[n] * val[n]; - S2 += wgt[n]; - } - double Value = S1 / S2; - *Sigma = 1.0 / S2; - - return Value; -} - Index: trunk/Ohana/src/opihi/cmd.data/write_vectors.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/write_vectors.c (revision 41323) +++ trunk/Ohana/src/opihi/cmd.data/write_vectors.c (revision 41340) @@ -32,23 +32,12 @@ } - // option generate a FITS output table (FITS holds the filename) - // in FITS output context, -header is interpretted as a buffer containing - // header keywords to supplement the FITS table header + /* option generate a FITS output table */ char *FITS = NULL; - Header *fitsheader = NULL; - Buffer *headbuffer = NULL; if ((N = get_argument (argc, argv, "-fits"))) { remove_argument (N, &argc, argv); FITS = strcreate (argv[N]); remove_argument (N, &argc, argv); - - if ((N = get_argument (argc, argv, "-header"))) { - remove_argument (N, &argc, argv); - if ((headbuffer = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE); - fitsheader = &headbuffer->header; - remove_argument (N, &argc, argv); - } - } - + } + /* option generate a FITS output table */ int CSV = FALSE; @@ -135,5 +124,5 @@ if (FITS) { - int status = WriteVectorTableFITS (argv[1], FITS, fitsheader, vec, Nvec, append, compress, format, Ntile); + int status = WriteVectorTableFITS (argv[1], FITS, vec, Nvec, append, compress, format, Ntile); free (vec); return status; @@ -154,8 +143,8 @@ if (ADD_HEADER) { for (j = 0; j < Nvec; j++) { - if (j == 0) fprintf (f, "# "); if (CSV) { fprintf (f, "%s,", vec[j][0].name); } else { + if (j == 0) fprintf (f, "# "); fprintf (f, "%s ", vec[j][0].name); } @@ -166,17 +155,18 @@ /* default output format */ if (format == (char *) NULL) { - char padChar = CSV ? ',' : ' '; for (i = 0; i < vec[0][0].Nelements; i++) { for (j = 0; j < Nvec; j++) { - switch (vec[j][0].type) { - case OPIHI_FLT: - fprintf (f, "%.12g%c", vec[j][0].elements.Flt[i], padChar); - break; - case OPIHI_INT: - fprintf (f, OPIHI_INT_FMT"%c", vec[j][0].elements.Int[i], padChar); - break; - case OPIHI_STR: - fprintf (f, "%s%c", vec[j][0].elements.Str[i], padChar); - break; + if (vec[j][0].type == OPIHI_FLT) { + if (CSV) { + fprintf (f, "%.12g,", vec[j][0].elements.Flt[i]); + } else { + fprintf (f, "%.12g ", vec[j][0].elements.Flt[i]); + } + } else { + if (CSV) { + fprintf (f, OPIHI_INT_FMT",", vec[j][0].elements.Int[i]); + } else { + fprintf (f, OPIHI_INT_FMT" ", vec[j][0].elements.Int[i]); + } } } @@ -225,9 +215,6 @@ fmttype[j] = 'd'; break; - case 's': - fmttype[j] = 's'; - break; default: - gprint (GP_ERR, "syntax error in format (only e,f,d,c,x,s allowed)\n"); + gprint (GP_ERR, "syntax error in format (only e,f,d,c,x allowed)\n"); return (FALSE); } @@ -236,23 +223,4 @@ strcat (fmtlist[Nvec-1], p0); - // check format types against vector types: - for (j = 0; j < Nvec; j++) { - switch (vec[j][0].type) { - case OPIHI_FLT: - case OPIHI_INT: - if (fmttype[j] == 's') { - gprint (GP_ERR, "mismatch between string format and numerical vector for %s\n", vec[j][0].name); - return FALSE; - } - break; - case OPIHI_STR: - if (fmttype[j] != 's') { - gprint (GP_ERR, "mismatch between numerical format and string vector for %s\n", vec[j][0].name); - return FALSE; - } - break; - } - } - for (i = 0; i < vec[0][0].Nelements; i++) { for (j = 0; j < Nvec; j++) { @@ -271,7 +239,4 @@ } } - if (fmttype[j] == 's') { - fprintf (f, fmtlist[j], vec[j][0].elements.Str[i]); - } } fprintf (f, "\n"); @@ -295,6 +260,4 @@ gprint (GP_ERR, " -append : write to the end of the existing file\n"); gprint (GP_ERR, " -fits NAME : write a fits table (extention name is NAME, column names match vector names)\n"); - gprint (GP_ERR, " in FITS output context, -header takes an additional argument which is interpretted\n"); - gprint (GP_ERR, " as a buffer containing header keywords to supplement the FITS table header\n"); gprint (GP_ERR, " -csv : write a comma-separated values file (eg, to read in excel)\n"); gprint (GP_ERR, " -f \"format\" : provide formatting codes for output:\n");