Index: trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.data/Makefile (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/Makefile (revision 41164) @@ -19,4 +19,5 @@ $(SRC)/init.$(ARCH).o \ $(SRC)/accum.$(ARCH).o \ +$(SRC)/antialias.$(ARCH).o \ $(SRC)/applyfit1d.$(ARCH).o \ $(SRC)/applyfit2d.$(ARCH).o \ @@ -78,4 +79,5 @@ $(SRC)/join.$(ARCH).o \ $(SRC)/jpeg.$(ARCH).o \ +$(SRC)/kapamemory.$(ARCH).o \ $(SRC)/kern.$(ARCH).o \ $(SRC)/keyword.$(ARCH).o \ @@ -89,8 +91,10 @@ $(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 \ @@ -99,4 +103,5 @@ $(SRC)/medimage.$(ARCH).o \ $(SRC)/medimage_commands.$(ARCH).o \ +$(SRC)/medimage_calc.$(ARCH).o \ $(SRC)/mset.$(ARCH).o \ $(SRC)/needles.$(ARCH).o \ @@ -169,4 +174,5 @@ $(SRC)/vgauss.$(ARCH).o \ $(SRC)/vlorentz.$(ARCH).o \ +$(SRC)/vsigmoid.$(ARCH).o \ $(SRC)/vellipse.$(ARCH).o \ $(SRC)/vmaxwell.$(ARCH).o \ @@ -176,4 +182,5 @@ $(SRC)/vzload.$(ARCH).o \ $(SRC)/vpeaks.$(ARCH).o \ +$(SRC)/vtransitions.$(ARCH).o \ $(SRC)/vpop.$(ARCH).o \ $(SRC)/vroll.$(ARCH).o \ @@ -182,4 +189,6 @@ $(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 41164) +++ trunk/Ohana/src/opihi/cmd.data/antialias.c (revision 41164) @@ -0,0 +1,26 @@ +# 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/cut.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/cut.c (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/cut.c (revision 41164) @@ -1,10 +1,71 @@ # include "data.h" -enum {SUM, MEAN, MEDIAN}; +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; +} int cut (int argc, char **argv) { int i, j, N, Nx, Ny, Mode; - float *Vin, *Vbuf, value; + float *Vin, *Vbuf; int sx, sy, nx, ny; Vector *xvec, *yvec; @@ -12,4 +73,8 @@ 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); @@ -19,4 +84,12 @@ 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; } @@ -36,6 +109,7 @@ Ny = buf[0].matrix.Naxis[1]; - if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) { - gprint (GP_ERR, "region out of range\n"); + // 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); return (FALSE); } @@ -57,21 +131,31 @@ 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; + } + if (i + sx >= Nx) { + yvec[0].elements.Flt[i] = NAN; + continue; + } + /* 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) { - Vbuf[j] = *Vin; - } - /* 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; } - } - yvec[0].elements.Flt[i] = value; + if (!isfinite(*Vin)) continue; + Vbuf[Npix] = *Vin; + Npix ++; + } + yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode); } free (Vbuf); @@ -89,21 +173,31 @@ 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; + } + if (i + sy >= Ny) { + yvec[0].elements.Flt[i] = NAN; + continue; + } + /* 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 ++) { - Vbuf[j] = *Vin; - } - /* 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; } - } - yvec[0].elements.Flt[i] = value; + if (!isfinite(*Vin)) continue; + Vbuf[Npix] = *Vin; + Npix ++; + } + yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode); } free (Vbuf); Index: trunk/Ohana/src/opihi/cmd.data/extract.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/extract.c (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/extract.c (revision 41164) @@ -1,11 +1,25 @@ # 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 i, j; - float *Vin, *Vout; - int sx, sy, nx, ny, NX, NY; - int Sx, Sy, Nx, Ny; + int N; 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) { @@ -15,22 +29,24 @@ if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - NX = in[0].matrix.Naxis[0]; - NY = in[0].matrix.Naxis[1]; + int NX = in[0].matrix.Naxis[0]; + int NY = in[0].matrix.Naxis[1]; - sx = atof (argv[3]); - sy = atof (argv[4]); - nx = atof (argv[5]); - ny = atof (argv[6]); + int sx = atof (argv[3]); + int sy = atof (argv[4]); + int nx = atof (argv[5]); + int ny = atof (argv[6]); - Sx = atof (argv[7]); - Sy = atof (argv[8]); - Nx = atof (argv[9]); - Ny = atof (argv[10]); + int Sx = atof (argv[7]); + int Sy = atof (argv[8]); + int Nx = atof (argv[9]); + int Ny = atof (argv[10]); if ((Sy + ny > Ny) || (Sx + nx > Nx)) { - gprint (GP_ERR, "mismatch between source and dest regions\n"); + gprint (GP_ERR, "source pixels extend beyond target pixels\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 */ @@ -39,14 +55,10 @@ (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); } @@ -76,13 +88,61 @@ } - 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++) { + // (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++) { if (i + sx < 0) continue; if (i + sx >= NX) continue; - *Vout = *Vin; + 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]; } } @@ -91,3 +151,2 @@ } - Index: trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/init.c (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/init.c (revision 41164) @@ -3,4 +3,5 @@ int accum PROTO((int, char **)); +int antialias PROTO((int, char **)); int applyfit PROTO((int, char **)); int applyfit1d PROTO((int, char **)); @@ -67,4 +68,5 @@ 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 **)); @@ -79,8 +81,10 @@ 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 **)); @@ -156,4 +160,5 @@ 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 **)); @@ -163,7 +168,11 @@ 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 **)); @@ -187,4 +196,5 @@ 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"}, @@ -255,4 +265,5 @@ {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"}, @@ -265,4 +276,5 @@ {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"}, @@ -273,4 +285,5 @@ {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)"}, @@ -347,4 +360,5 @@ {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"}, @@ -354,5 +368,6 @@ {1, "vload", vload, "load vectors as overlay on image display"}, {1, "vmaxwell", vmaxwell, "fit a Maxwellian to a vector"}, - {1, "vpeaks", vpeaks, "fine coord and flux of peaks in 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, "vpop", vpop, "remove first element of a vector"}, {1, "vroll", vroll, "roll vector elements by 1 entry"}, @@ -360,4 +375,6 @@ {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 40654) +++ trunk/Ohana/src/opihi/cmd.data/interpolate_presort.c (revision 41164) @@ -7,4 +7,11 @@ 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 **/ @@ -36,10 +43,12 @@ y0 = yin[0].elements.Flt[0]; - /* in and out vectors are sorted */ - j = 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; + } // every output pixel should get a value unless // we are below the lowest in or above the highest in - for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) { + for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) { yout[0].elements.Flt[i] = NAN; @@ -66,10 +75,10 @@ // fill in the rest with NANs + int NinLast = xin[0].Nelements - 1; while (i < yout[0].Nelements) { - yout[0].elements.Flt[i] = NAN; + yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN; i++; } return (TRUE); - } Index: trunk/Ohana/src/opihi/cmd.data/kapamemory.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/kapamemory.c (revision 41164) +++ trunk/Ohana/src/opihi/cmd.data/kapamemory.c (revision 41164) @@ -0,0 +1,56 @@ +# 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/match1d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/match1d.c (revision 41164) +++ trunk/Ohana/src/opihi/cmd.data/match1d.c (revision 41164) @@ -0,0 +1,248 @@ +# 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 41164) +++ trunk/Ohana/src/opihi/cmd.data/medimage_calc.c (revision 41164) @@ -0,0 +1,365 @@ +# 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 40654) +++ trunk/Ohana/src/opihi/cmd.data/medimage_commands.c (revision 41164) @@ -15,13 +15,30 @@ 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)\n"); + gprint (GP_ERR, "USAGE: medimage add (name) (image) [-variance variance]\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]); @@ -46,76 +63,17 @@ // new image should match existing medimage dimensions + // AddMedImage (median, image, var); int Ninput = median->Ninput; median->Ninput ++; - REALLOCATE (median->buffers, float *, median->Ninput); + REALLOCATE (median->flx, float *, median->Ninput); + REALLOCATE (median->var, 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); + ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny); + memcpy (median->flx[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)]; - } - } + 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); } @@ -124,39 +82,39 @@ /* -int medimage_save (int argc, char **argv) { + int medimage_save (int argc, char **argv) { - int N; + int N; - int APPEND = FALSE; - if ((N = get_argument (argc, argv, "-append"))) { - APPEND = TRUE; - remove_argument (N, &argc, argv); - } + 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 (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; -} + 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) { + int medimage_load (int argc, char **argv) { - if (argc != 3) { - gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n"); - return FALSE; - } + 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; -} + 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 41164) +++ trunk/Ohana/src/opihi/cmd.data/mgaussdev.c (revision 41164) @@ -0,0 +1,34 @@ +# 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/peak.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/peak.c (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/peak.c (revision 41164) @@ -39,4 +39,5 @@ if (*X < start) continue; if (*X > end) continue; + if (!isfinite(*Y)) continue; if (*Y < ymax) continue; xmax = *X; @@ -63,4 +64,5 @@ 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/roll.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/roll.c (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/roll.c (revision 41164) @@ -3,7 +3,9 @@ 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) { @@ -14,32 +16,71 @@ if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - dx = atof (argv[2]); - dy = atof (argv[3]); - if (dy != 0) { - gprint (GP_ERR, "only x rolls implemented for now\n"); + 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; } - /* 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]; + if (dx == 0) { + // moves in just dy can be done row-by-row - dX = abs(dx); - Nbytes = nx * ny * sizeof (float); - Nextra = (nx * ny + dX) * sizeof (float); + ALLOCATE_PTR (output, float, Nx*Ny); - REALLOCATE (buf[0].matrix.buffer, char, Nextra); + float *input = (float *) buf[0].matrix.buffer; - 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); + 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; } - REALLOCATE (buf[0].matrix.buffer, char, Nbytes); + ALLOCATE_PTR (output, float, Nx*Ny); + float *input = (float *) buf[0].matrix.buffer; - return (TRUE); - + 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; } + Index: trunk/Ohana/src/opihi/cmd.data/rotate.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/rotate.c (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/rotate.c (revision 41164) @@ -3,19 +3,19 @@ int rotate (int argc, char **argv) { - int i, j, NX, NY, X, Y, Lx, Ly, N; + int i, j, NX, NY; float *in_buff, *out_buff, *c; - double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, x, y, X1, Y1; + double Xo, Yo, 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,30 +221,26 @@ } - angle = atof (argv[2]); - CosAngle = cos (angle*RAD_DEG); - SinAngle = sin (angle*RAD_DEG); + double angle = atof (argv[2]); + double CosAngle = cos (angle*RAD_DEG); + double SinAngle = sin (angle*RAD_DEG); gprint (GP_ERR, "rotating: %f %f %f\n", angle, CosAngle, SinAngle); - 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); - */ + // 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); /* fix reference pixel */ gfits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo); gfits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo); - /* - 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; + X1 = (Xo - NX/2)*CosAngle + (Yo - NY/2)*SinAngle + Lx/2; + Y1 = (NX/2 - Xo)*SinAngle + (Yo - NY/2)*CosAngle + Ly/2; gfits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1); gfits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1); @@ -255,8 +251,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); @@ -270,22 +266,32 @@ 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++) { - 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; + 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; + *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 40654) +++ trunk/Ohana/src/opihi/cmd.data/spline.c (revision 41164) @@ -7,4 +7,5 @@ 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); @@ -18,4 +19,5 @@ {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"}, @@ -35,4 +37,5 @@ 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 40654) +++ trunk/Ohana/src/opihi/cmd.data/spline_commands.c (revision 41164) @@ -28,4 +28,9 @@ 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; } @@ -93,4 +98,26 @@ } +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) { @@ -109,6 +136,11 @@ int spline_delete (int argc, char **argv) { - int status; - Spline *spline; + int N, status; + + int QUIET = FALSE; + if ((N = get_argument (argc, argv, "-q"))) { + QUIET = TRUE; + remove_argument (N, &argc, argv); + } if (argc != 2) { @@ -117,6 +149,7 @@ } - spline = FindSpline (argv[1]); + Spline *spline = FindSpline (argv[1]); if (spline == NULL) { + if (QUIET) return TRUE; gprint (GP_ERR, "spline %s not found\n", argv[1]); return FALSE; @@ -130,6 +163,4 @@ int spline_rename (int argc, char **argv) { - Spline *spline; - if (argc != 3) { gprint (GP_ERR, "USAGE: spline rename (spline) (newname)\n"); @@ -137,8 +168,13 @@ } - spline = FindSpline (argv[1]); + Spline *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); } @@ -147,24 +183,2 @@ 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 40654) +++ trunk/Ohana/src/opihi/cmd.data/test/cut.sh (revision 41164) @@ -4,4 +4,46 @@ 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 @@ -37,2 +79,42 @@ 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 40654) +++ trunk/Ohana/src/opihi/cmd.data/test/delete.sh (revision 41164) @@ -49,2 +49,13 @@ 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 41164) +++ trunk/Ohana/src/opihi/cmd.data/test/headtest.sh (revision 41164) @@ -0,0 +1,15 @@ + +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 40654) +++ trunk/Ohana/src/opihi/cmd.data/test/medimage.sh (revision 41164) @@ -1,11 +1,540 @@ -macro go - mcreate a 30 30 - for i 0 40 - set a$i = zero(a) + $i - end - - for i 0 40 - medimage add t1 a$i - end -end +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 + 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 +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 41164) +++ trunk/Ohana/src/opihi/cmd.data/test/set.sh (revision 41164) @@ -0,0 +1,43 @@ + +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 41164) +++ trunk/Ohana/src/opihi/cmd.data/test/spline.sh (revision 41164) @@ -0,0 +1,11 @@ + +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 41164) +++ trunk/Ohana/src/opihi/cmd.data/test/vsigmoid.sh (revision 41164) @@ -0,0 +1,263 @@ + +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 41164) +++ trunk/Ohana/src/opihi/cmd.data/test/vstat.sh (revision 41164) @@ -0,0 +1,28 @@ + +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 40654) +++ trunk/Ohana/src/opihi/cmd.data/threshold.c (revision 41164) @@ -25,4 +25,14 @@ 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); } @@ -64,4 +74,16 @@ } + // 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 @@ -71,4 +93,13 @@ 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); @@ -89,4 +120,13 @@ } } 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/vgauss.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vgauss.c (revision 40654) +++ trunk/Ohana/src/opihi/cmd.data/vgauss.c (revision 41164) @@ -3,4 +3,5 @@ /* 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) \ @@ -20,4 +21,10 @@ 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; @@ -104,4 +111,6 @@ } 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++) { @@ -125,4 +134,38 @@ } +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) { @@ -135,8 +178,10 @@ f = par[2]*r + par[3]; - dpar[0] = par[2]*r*z/par[1]; - dpar[1] = par[2]*r*z*z/par[1]; - dpar[2] = r; - dpar[3] = 1; + 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; + } return (f); Index: trunk/Ohana/src/opihi/cmd.data/virls.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/virls.c (revision 41164) +++ trunk/Ohana/src/opihi/cmd.data/virls.c (revision 41164) @@ -0,0 +1,103 @@ +# 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 41164) +++ trunk/Ohana/src/opihi/cmd.data/vsigmoid.c (revision 41164) @@ -0,0 +1,192 @@ +# 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 41164) +++ trunk/Ohana/src/opihi/cmd.data/vtransitions.c (revision 41164) @@ -0,0 +1,175 @@ +# 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 41164) +++ trunk/Ohana/src/opihi/cmd.data/vwtmean.c (revision 41164) @@ -0,0 +1,70 @@ +# 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 40654) +++ trunk/Ohana/src/opihi/cmd.data/write_vectors.c (revision 41164) @@ -32,12 +32,23 @@ } - /* option generate a FITS output table */ + // 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 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; @@ -124,5 +135,5 @@ if (FITS) { - int status = WriteVectorTableFITS (argv[1], FITS, vec, Nvec, append, compress, format, Ntile); + int status = WriteVectorTableFITS (argv[1], FITS, fitsheader, vec, Nvec, append, compress, format, Ntile); free (vec); return status; @@ -260,4 +271,6 @@ 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");