Index: trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.data/Makefile (revision 21064) +++ trunk/Ohana/src/opihi/cmd.data/Makefile (revision 21508) @@ -20,5 +20,5 @@ $(SRC)/accum.$(ARCH).o \ $(SRC)/applyfit2d.$(ARCH).o \ -$(SRC)/applyfit.$(ARCH).o \ +$(SRC)/applyfit1d.$(ARCH).o \ $(SRC)/box.$(ARCH).o \ $(SRC)/book.$(ARCH).o \ @@ -26,5 +26,5 @@ $(SRC)/center.$(ARCH).o \ $(SRC)/clear.$(ARCH).o \ -$(SRC)/clip.$(ARCH).o \ +$(SRC)/imclip.$(ARCH).o \ $(SRC)/close.$(ARCH).o \ $(SRC)/concat.$(ARCH).o \ @@ -44,6 +44,6 @@ $(SRC)/fft1d.$(ARCH).o \ $(SRC)/fft2d.$(ARCH).o \ +$(SRC)/fit1d.$(ARCH).o \ $(SRC)/fit2d.$(ARCH).o \ -$(SRC)/fit.$(ARCH).o \ $(SRC)/gaussj.$(ARCH).o \ $(SRC)/gaussdeviate.$(ARCH).o \ @@ -107,5 +107,5 @@ $(SRC)/spline_apply.$(ARCH).o \ $(SRC)/spline_construct.$(ARCH).o \ -$(SRC)/stats.$(ARCH).o \ +$(SRC)/imstats.$(ARCH).o \ $(SRC)/style.$(ARCH).o \ $(SRC)/subraster.$(ARCH).o \ @@ -131,5 +131,5 @@ $(SRC)/vroll.$(ARCH).o \ $(SRC)/vsmooth.$(ARCH).o \ -$(SRC)/vstat.$(ARCH).o \ +$(SRC)/vstats.$(ARCH).o \ $(SRC)/wd.$(ARCH).o \ $(SRC)/write_vectors.$(ARCH).o \ Index: trunk/Ohana/src/opihi/cmd.data/applyfit.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/applyfit.c (revision 21064) +++ (revision ) @@ -1,66 +1,0 @@ -# include "data.h" - -int applyfit (int argc, char **argv) { - - int i, j, order; - char *c, name[64]; - double *C, X; - opihi_flt *y; - Vector *xvec, *yvec; - - if (argc != 3) { - gprint (GP_ERR, "USAGE: applyfit x y\n"); - return (FALSE); - } - - c = get_variable ("Cn"); - if (c == NULL) { - gprint (GP_ERR, "no fit available\n"); - return (FALSE); - } - order = atof (c); - free (c); - - if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((yvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); - - ALLOCATE (C, double, order+1); - for (i = 0; i < order + 1; i++) { - sprintf (name, "C%d", i); - c = get_variable (name); - if (c == NULL) { - gprint (GP_ERR, "missing fit term %d\n", i); - return (FALSE); - } - C[i] = atof (c); - free (c); - } - - ResetVector (yvec, OPIHI_FLT, xvec[0].Nelements); - bzero (yvec[0].elements.Flt, sizeof(opihi_flt)*yvec[0].Nelements); - y = yvec[0].elements.Flt; - - if (xvec[0].type == OPIHI_FLT) { - opihi_flt *x = xvec[0].elements.Flt; - for (j = 0; j < xvec[0].Nelements; j++, x++, y++) { - X = 1; - for (i = 0; i < order + 1; i++) { - *y += C[i]*X; - X = X * (*x); - } - } - } else { - opihi_int *x = xvec[0].elements.Int; - for (j = 0; j < xvec[0].Nelements; j++, x++, y++) { - X = 1; - for (i = 0; i < order + 1; i++) { - *y += C[i]*X; - X = X * (*x); - } - } - } - - free (C); - return (TRUE); - -} Index: trunk/Ohana/src/opihi/cmd.data/applyfit1d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/applyfit1d.c (revision 21508) +++ trunk/Ohana/src/opihi/cmd.data/applyfit1d.c (revision 21508) @@ -0,0 +1,66 @@ +# include "data.h" + +int applyfit1d (int argc, char **argv) { + + int i, j, order; + char *c, name[64]; + double *C, X; + opihi_flt *y; + Vector *xvec, *yvec; + + if (argc != 3) { + gprint (GP_ERR, "USAGE: applyfit x y\n"); + return (FALSE); + } + + c = get_variable ("Cn"); + if (c == NULL) { + gprint (GP_ERR, "no fit available\n"); + return (FALSE); + } + order = atof (c); + free (c); + + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + ALLOCATE (C, double, order+1); + for (i = 0; i < order + 1; i++) { + sprintf (name, "C%d", i); + c = get_variable (name); + if (c == NULL) { + gprint (GP_ERR, "missing fit term %d\n", i); + return (FALSE); + } + C[i] = atof (c); + free (c); + } + + ResetVector (yvec, OPIHI_FLT, xvec[0].Nelements); + bzero (yvec[0].elements.Flt, sizeof(opihi_flt)*yvec[0].Nelements); + y = yvec[0].elements.Flt; + + if (xvec[0].type == OPIHI_FLT) { + opihi_flt *x = xvec[0].elements.Flt; + for (j = 0; j < xvec[0].Nelements; j++, x++, y++) { + X = 1; + for (i = 0; i < order + 1; i++) { + *y += C[i]*X; + X = X * (*x); + } + } + } else { + opihi_int *x = xvec[0].elements.Int; + for (j = 0; j < xvec[0].Nelements; j++, x++, y++) { + X = 1; + for (i = 0; i < order + 1; i++) { + *y += C[i]*X; + X = X * (*x); + } + } + } + + free (C); + return (TRUE); + +} Index: trunk/Ohana/src/opihi/cmd.data/clip.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/clip.c (revision 21064) +++ (revision ) @@ -1,73 +1,0 @@ -# include "data.h" - -int clip (int argc, char **argv) { - - int i, Npix, DO_NAN, DO_INF, N; - double min, Vmin, max, Vmax, nan_val, inf_val; - float *in; - Buffer *buf; - - inf_val = nan_val = min = Vmin = max = Vmax = 0; - - DO_NAN = FALSE; - if ((N = get_argument (argc, argv, "-nan"))) { - remove_argument (N, &argc, argv); - nan_val = atof(argv[N]); - remove_argument (N, &argc, argv); - DO_NAN = TRUE; - } - - DO_INF = FALSE; - if ((N = get_argument (argc, argv, "-inf"))) { - remove_argument (N, &argc, argv); - inf_val = atof(argv[N]); - remove_argument (N, &argc, argv); - DO_INF = TRUE; - } - - if ((argc != 6) && (!(DO_INF || DO_NAN))) { - gprint (GP_ERR, "USAGE: clip (buffer) [min Vmin max Vmax] [-inf val] [-nan val]\n"); - return (FALSE); - } - - if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - - if (argc == 6) { - min = atof (argv[2]); - Vmin = atof (argv[3]); - max = atof (argv[4]); - Vmax = atof (argv[5]); - } - - Npix = buf[0].matrix.Naxis[0]*buf[0].matrix.Naxis[1]; - in = (float *) buf[0].matrix.buffer; - - if (argc == 6) { - for (i = 0; i < Npix; i++, in++) { - if (*in < min) - *in = Vmin; - if (*in > max) - *in = Vmax; - } - } - in = (float *) buf[0].matrix.buffer; - if (DO_NAN) { - for (i = 0; i < Npix; i++, in++) { - if (isnan (*in)) { - *in = nan_val; - } - } - } - in = (float *) buf[0].matrix.buffer; - if (DO_INF) { - for (i = 0; i < Npix; i++, in++) { - if (!finite (*in)) { - *in = inf_val; - } - } - } - return (TRUE); - -} - - Index: trunk/Ohana/src/opihi/cmd.data/fit.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/fit.c (revision 21064) +++ (revision ) @@ -1,234 +1,0 @@ -# include "data.h" - -int fit (int argc, char **argv) { - - double **c, **b, *s, X, Y, dY, dY2; - double ClipNSigma, mean, sigma, maxsigma; - int i, j, nterm, mterm, order, Npt, Nmask; - int N, Weight, Quiet, ClipNiter; - Vector *xvec, *yvec, *dyvec; - opihi_flt *x, *y, *dy, *yf, *yfit; - char name[64], *mask; - - 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); - } - - ClipNSigma = 0; - ClipNiter = 1; - if ((N = get_argument (argc, argv, "-clip"))) { - remove_argument (N, &argc, argv); - ClipNSigma = atof(argv[N]); - remove_argument (N, &argc, argv); - ClipNiter = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - dy = NULL; - dyvec = NULL; - Weight = FALSE; - if ((N = get_argument (argc, argv, "-dy"))) { - remove_argument (N, &argc, argv); - if ((dyvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - Weight = TRUE; - } - - if (argc != 4) { - gprint (GP_ERR, "USAGE: fit x y order [-dy wt] [-quiet/-q] [-clip Nsigma Niter]\n"); - return (FALSE); - } - - if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); - - if (xvec[0].Nelements != yvec[0].Nelements) { - gprint (GP_ERR, "vectors must have same length\n"); - return (FALSE); - } - REQUIRE_VECTOR_FLT (xvec, FALSE); - REQUIRE_VECTOR_FLT (yvec, FALSE); - - if (Weight) { - REQUIRE_VECTOR_FLT (dyvec, FALSE); - if (xvec[0].Nelements != dyvec[0].Nelements) { - gprint (GP_ERR, "vectors must have same length\n"); - return (FALSE); - } - } - - - /* nterm is number of polynomial terms, starting at x^0 */ - order = atof (argv[3]); - nterm = order + 1; - mterm = 2*nterm; - - ALLOCATE (yfit, opihi_flt, xvec[0].Nelements); - ALLOCATE (mask, char, xvec[0].Nelements); - memset (mask, 0, xvec[0].Nelements); - - ALLOCATE (s, double, mterm); - ALLOCATE (b, double *, nterm); - ALLOCATE (c, double *, nterm); - for (i = 0; i < nterm; i++) { - ALLOCATE (c[i], double, nterm); - ALLOCATE (b[i], double, 1); - } - - Nmask = 0; - sigma = 0.0; - - for (N = 0; N < ClipNiter; N++) { - - /* init registers for current pass */ - memset (s, 0, mterm*sizeof(double)); - for (i = 0; i < nterm; i++) { - memset (c[i], 0, nterm*sizeof(double)); - memset (b[i], 0, sizeof(double)); - } - - /* perform linear fit */ - x = xvec[0].elements.Flt; - y = yvec[0].elements.Flt; - if (Weight) dy = dyvec[0].elements.Flt; - - for (i = 0; i < xvec[0].Nelements; i++, x++, y++) { - if (mask[i]) continue; - if (!(finite(*x) && finite(*y))) continue; - dY = 1.0; - if (Weight) { - dY = 1.0 / SQ(*dy); - dy ++; - } - X = 1*dY; - Y = *y*dY; - for (j = 0; j < nterm; j++) { - s[j] += X; - b[j][0] += Y; - X = X * (*x); - Y = Y * (*x); - } - for (j = nterm; j < mterm; j++) { - s[j] += X; - X = X * (*x); - } - } - for (i = 0; i < nterm; i++) { - for (j = 0; j < nterm; j++) { - c[i][j] = s[i + j]; - } - } - if (!dgaussjordan (c, b, nterm, 1)) goto escape; - - /* generate fitted values */ - x = xvec[0].elements.Flt; - yf = yfit; - for (i = 0; i < xvec[0].Nelements; i++, x++, yf++) { - if (!finite(*x)) continue; - *yf = 0; - X = 1; - for (j = 0; j < order + 1; j++) { - *yf += b[j][0]*X; - X = X * (*x); - } - } - - /* measure fit residual scatter */ - x = xvec[0].elements.Flt; - y = yvec[0].elements.Flt; - yf = yfit; - dY = dY2 = 0; - for (i = Npt = 0; i < xvec[0].Nelements; i++, x++, y++, yf++) { - if (mask[i]) continue; - if (!finite(*x)) continue; - dY += (*y - *yf); - dY2 += SQ(*y - *yf); - Npt ++; - } - mean = dY / Npt; - sigma = sqrt (fabs(dY2/Npt - SQ(mean))); - maxsigma = ClipNSigma * sigma; - - /* mask outlier points */ - x = xvec[0].elements.Flt; - y = yvec[0].elements.Flt; - yf = yfit; - Nmask = 0; - for (i = 0; ClipNSigma && (i < xvec[0].Nelements); i++, x++, y++, yf++) { - dY = (*y - *yf); - if (fabs(dY) > maxsigma) { - mask[i] = TRUE; - Nmask ++; - } else { - mask[i] = FALSE; - } - } - } - - /* print & save basic fit parameters */ - if (!Quiet) gprint (GP_ERR, "y = "); - for (i = 0; i < nterm; i++) { - sprintf (name, "C%d", i); - set_variable (name, b[i][0]); - if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i); - } - sprintf (name, "Cn"); - set_variable (name, (double) order); - - /* print & save basic fit parameters */ - if (!Quiet) gprint (GP_ERR, "\n"); - if (!Quiet) gprint (GP_ERR, " "); - for (i = 0; i < nterm; i++) { - sprintf (name, "dC%d", i); - set_variable (name, sqrt(c[i][i])); - if (!Quiet) gprint (GP_ERR, "%f ", sqrt(c[i][i])); - } - if (!Quiet) gprint (GP_ERR, "\n"); - - set_variable ("dC", sigma); - set_variable ("Cnv", (xvec[0].Nelements - Nmask)); - - /* save mask and yfit for testing? */ - if (1) { - Vector *fvec, *mvec; - if ((fvec = SelectVector ("yfit", ANYVECTOR, TRUE)) == NULL) return (FALSE); - if ((mvec = SelectVector ("mask", ANYVECTOR, TRUE)) == NULL) return (FALSE); - free (fvec[0].elements.Flt); - fvec[0].elements.Flt = yfit; - fvec[0].Nelements = xvec[0].Nelements; - mvec[0].Nelements = xvec[0].Nelements; - - REALLOCATE (mvec[0].elements.Flt, opihi_flt, xvec[0].Nelements); - for (i = 0; i < xvec[0].Nelements; i++) { - mvec[0].elements.Flt[i] = mask[i]; - } - } else { - free (yfit); - } - free (mask); - - for (i = 0; i < nterm; i++) { - free (b[i]); - free (c[i]); - } - free (b); - free (c); - free (s); - return (TRUE); - -escape: - for (i = 0; i < nterm; i++) { - free (b[i]); - free (c[i]); - } - free (b); - free (c); - free (s); - return (FALSE); -} Index: trunk/Ohana/src/opihi/cmd.data/fit1d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/fit1d.c (revision 21508) +++ trunk/Ohana/src/opihi/cmd.data/fit1d.c (revision 21508) @@ -0,0 +1,234 @@ +# include "data.h" + +int fit1d (int argc, char **argv) { + + double **c, **b, *s, X, Y, dY, dY2; + double ClipNSigma, mean, sigma, maxsigma; + int i, j, nterm, mterm, order, Npt, Nmask; + int N, Weight, Quiet, ClipNiter; + Vector *xvec, *yvec, *dyvec; + opihi_flt *x, *y, *dy, *yf, *yfit; + char name[64], *mask; + + 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); + } + + ClipNSigma = 0; + ClipNiter = 1; + if ((N = get_argument (argc, argv, "-clip"))) { + remove_argument (N, &argc, argv); + ClipNSigma = atof(argv[N]); + remove_argument (N, &argc, argv); + ClipNiter = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + dy = NULL; + dyvec = NULL; + Weight = FALSE; + if ((N = get_argument (argc, argv, "-dy"))) { + remove_argument (N, &argc, argv); + if ((dyvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + Weight = TRUE; + } + + if (argc != 4) { + gprint (GP_ERR, "USAGE: fit x y order [-dy wt] [-quiet/-q] [-clip Nsigma Niter]\n"); + return (FALSE); + } + + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + if (xvec[0].Nelements != yvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + REQUIRE_VECTOR_FLT (xvec, FALSE); + REQUIRE_VECTOR_FLT (yvec, FALSE); + + if (Weight) { + REQUIRE_VECTOR_FLT (dyvec, FALSE); + if (xvec[0].Nelements != dyvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + } + + + /* nterm is number of polynomial terms, starting at x^0 */ + order = atof (argv[3]); + nterm = order + 1; + mterm = 2*nterm; + + ALLOCATE (yfit, opihi_flt, xvec[0].Nelements); + ALLOCATE (mask, char, xvec[0].Nelements); + memset (mask, 0, xvec[0].Nelements); + + ALLOCATE (s, double, mterm); + ALLOCATE (b, double *, nterm); + ALLOCATE (c, double *, nterm); + for (i = 0; i < nterm; i++) { + ALLOCATE (c[i], double, nterm); + ALLOCATE (b[i], double, 1); + } + + Nmask = 0; + sigma = 0.0; + + for (N = 0; N < ClipNiter; N++) { + + /* init registers for current pass */ + memset (s, 0, mterm*sizeof(double)); + for (i = 0; i < nterm; i++) { + memset (c[i], 0, nterm*sizeof(double)); + memset (b[i], 0, sizeof(double)); + } + + /* perform linear fit */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + if (Weight) dy = dyvec[0].elements.Flt; + + for (i = 0; i < xvec[0].Nelements; i++, x++, y++) { + if (mask[i]) continue; + if (!(finite(*x) && finite(*y))) continue; + dY = 1.0; + if (Weight) { + dY = 1.0 / SQ(*dy); + dy ++; + } + X = 1*dY; + Y = *y*dY; + for (j = 0; j < nterm; j++) { + s[j] += X; + b[j][0] += Y; + X = X * (*x); + Y = Y * (*x); + } + for (j = nterm; j < mterm; j++) { + s[j] += X; + X = X * (*x); + } + } + for (i = 0; i < nterm; i++) { + for (j = 0; j < nterm; j++) { + c[i][j] = s[i + j]; + } + } + if (!dgaussjordan (c, b, nterm, 1)) goto escape; + + /* generate fitted values */ + x = xvec[0].elements.Flt; + yf = yfit; + for (i = 0; i < xvec[0].Nelements; i++, x++, yf++) { + if (!finite(*x)) continue; + *yf = 0; + X = 1; + for (j = 0; j < order + 1; j++) { + *yf += b[j][0]*X; + X = X * (*x); + } + } + + /* measure fit residual scatter */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + yf = yfit; + dY = dY2 = 0; + for (i = Npt = 0; i < xvec[0].Nelements; i++, x++, y++, yf++) { + if (mask[i]) continue; + if (!finite(*x)) continue; + dY += (*y - *yf); + dY2 += SQ(*y - *yf); + Npt ++; + } + mean = dY / Npt; + sigma = sqrt (fabs(dY2/Npt - SQ(mean))); + maxsigma = ClipNSigma * sigma; + + /* mask outlier points */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + yf = yfit; + Nmask = 0; + for (i = 0; ClipNSigma && (i < xvec[0].Nelements); i++, x++, y++, yf++) { + dY = (*y - *yf); + if (fabs(dY) > maxsigma) { + mask[i] = TRUE; + Nmask ++; + } else { + mask[i] = FALSE; + } + } + } + + /* print & save basic fit parameters */ + if (!Quiet) gprint (GP_ERR, "y = "); + for (i = 0; i < nterm; i++) { + sprintf (name, "C%d", i); + set_variable (name, b[i][0]); + if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i); + } + sprintf (name, "Cn"); + set_variable (name, (double) order); + + /* print & save basic fit parameters */ + if (!Quiet) gprint (GP_ERR, "\n"); + if (!Quiet) gprint (GP_ERR, " "); + for (i = 0; i < nterm; i++) { + sprintf (name, "dC%d", i); + set_variable (name, sqrt(c[i][i])); + if (!Quiet) gprint (GP_ERR, "%f ", sqrt(c[i][i])); + } + if (!Quiet) gprint (GP_ERR, "\n"); + + set_variable ("dC", sigma); + set_variable ("Cnv", (xvec[0].Nelements - Nmask)); + + /* save mask and yfit for testing? */ + if (1) { + Vector *fvec, *mvec; + if ((fvec = SelectVector ("yfit", ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((mvec = SelectVector ("mask", ANYVECTOR, TRUE)) == NULL) return (FALSE); + free (fvec[0].elements.Flt); + fvec[0].elements.Flt = yfit; + fvec[0].Nelements = xvec[0].Nelements; + mvec[0].Nelements = xvec[0].Nelements; + + REALLOCATE (mvec[0].elements.Flt, opihi_flt, xvec[0].Nelements); + for (i = 0; i < xvec[0].Nelements; i++) { + mvec[0].elements.Flt[i] = mask[i]; + } + } else { + free (yfit); + } + free (mask); + + for (i = 0; i < nterm; i++) { + free (b[i]); + free (c[i]); + } + free (b); + free (c); + free (s); + return (TRUE); + +escape: + for (i = 0; i < nterm; i++) { + free (b[i]); + free (c[i]); + } + free (b); + free (c); + free (s); + return (FALSE); +} Index: trunk/Ohana/src/opihi/cmd.data/imclip.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/imclip.c (revision 21508) +++ trunk/Ohana/src/opihi/cmd.data/imclip.c (revision 21508) @@ -0,0 +1,73 @@ +# include "data.h" + +int imclip (int argc, char **argv) { + + int i, Npix, DO_NAN, DO_INF, N; + double min, Vmin, max, Vmax, nan_val, inf_val; + float *in; + Buffer *buf; + + inf_val = nan_val = min = Vmin = max = Vmax = 0; + + DO_NAN = FALSE; + if ((N = get_argument (argc, argv, "-nan"))) { + remove_argument (N, &argc, argv); + nan_val = atof(argv[N]); + remove_argument (N, &argc, argv); + DO_NAN = TRUE; + } + + DO_INF = FALSE; + if ((N = get_argument (argc, argv, "-inf"))) { + remove_argument (N, &argc, argv); + inf_val = atof(argv[N]); + remove_argument (N, &argc, argv); + DO_INF = TRUE; + } + + if ((argc != 6) && (!(DO_INF || DO_NAN))) { + gprint (GP_ERR, "USAGE: clip (buffer) [min Vmin max Vmax] [-inf val] [-nan val]\n"); + return (FALSE); + } + + if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + if (argc == 6) { + min = atof (argv[2]); + Vmin = atof (argv[3]); + max = atof (argv[4]); + Vmax = atof (argv[5]); + } + + Npix = buf[0].matrix.Naxis[0]*buf[0].matrix.Naxis[1]; + in = (float *) buf[0].matrix.buffer; + + if (argc == 6) { + for (i = 0; i < Npix; i++, in++) { + if (*in < min) + *in = Vmin; + if (*in > max) + *in = Vmax; + } + } + in = (float *) buf[0].matrix.buffer; + if (DO_NAN) { + for (i = 0; i < Npix; i++, in++) { + if (isnan (*in)) { + *in = nan_val; + } + } + } + in = (float *) buf[0].matrix.buffer; + if (DO_INF) { + for (i = 0; i < Npix; i++, in++) { + if (!finite (*in)) { + *in = inf_val; + } + } + } + return (TRUE); + +} + + Index: trunk/Ohana/src/opihi/cmd.data/imstats.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/imstats.c (revision 21508) +++ trunk/Ohana/src/opihi/cmd.data/imstats.c (revision 21508) @@ -0,0 +1,125 @@ +# include "data.h" + +int imstats (int argc, char **argv) { + + int i, j, Nmode, Imode; + double Npix, N1, N2, max, min, range, median, mode, IgnoreValue; + float *V; + int sx, sy, nx, ny, *hist, Nhist, bin; + int Ignore, Quiet, N; + Buffer *buf; + + IgnoreValue = 0; + Ignore = FALSE; + if ((N = get_argument (argc, argv, "-ignore"))) { + Ignore = TRUE; + remove_argument (N, &argc, argv); + IgnoreValue = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + 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 != 2) && (argc != 6)) { + gprint (GP_ERR, "USAGE: stats sx sy nx ny\n"); + gprint (GP_ERR, "OR: stats \n"); + return (FALSE); + } + + if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + if (argc == 6) { + sx = strcmp (argv[2], "-") ? atof (argv[2]) : 0; + sy = strcmp (argv[3], "-") ? atof (argv[3]) : 0; + nx = strcmp (argv[4], "-") ? atof (argv[4]) : buf[0].matrix.Naxis[0]; + ny = strcmp (argv[5], "-") ? atof (argv[5]) : buf[0].matrix.Naxis[1]; + } else { + sx = 0; + sy = 0; + nx = buf[0].matrix.Naxis[0]; + ny = buf[0].matrix.Naxis[1]; + } + + Npix = N1 = N2 = 0; + if ((sx < 0) || (sy < 0) || + (sx+nx > buf[0].matrix.Naxis[0]) || + (sy+ny > buf[0].matrix.Naxis[1])) { + gprint (GP_ERR, "region out of range\n"); + return (FALSE); + } + + max = -1e32; + min = +1e32; + for (j = sy; j < sy + ny; j++) { + V = (float *)(buf[0].matrix.buffer) + j*buf[0].matrix.Naxis[0] + sx; + for (i = 0; i < nx; i++, V++) { + if (Ignore && (fabs (*V - IgnoreValue) < 1e-8)) continue; + if (isnan(*V)) continue; + if (isinf(*V)) continue; + N1 += *V; + N2 += (*V)*(*V); + Npix += 1.0; + max = MAX (max, *V); + min = MIN (min, *V); + } + } + N1 = N1 / Npix; + +/* calculate mode, median */ + median = mode = 0.5*(max + min); + if ((max - min) != 0) { + range = 0xffff / (max - min); + ALLOCATE (hist, int, 0x10000); + bzero (hist, 0x10000*sizeof(int)); + for (j = sy; j < sy + ny; j++) { + V = (float *)(buf[0].matrix.buffer) + j*buf[0].matrix.Naxis[0] + sx; + for (i = 0; i < nx; i++, V++) { + if (Ignore && (fabs (*V - IgnoreValue) < 1e-8)) continue; + if (isnan(*V)) continue; + if (isinf(*V)) continue; + bin = MIN (MAX (0, (*V - min) * range), 0xffff); + hist[bin] ++; + } + } + Nhist = 0; + for (i = 0; (i < 0xffff) && (Nhist < 0.5*Npix); i++) + Nhist += hist[i]; + median = i / range + min; + Nmode = hist[0]; + Imode = 0; + for (i = 1; i < 0x10000; i++) { + if (hist[i] > Nmode) { + Nmode = hist[i]; + Imode = i; + } + } + mode = Imode / range + min; + free (hist); + } + + if (!Quiet) { + gprint (GP_LOG, " mean stdev min max median Npix Total\n"); + gprint (GP_LOG, "%7.4g %7.4g %7.4g %7.4g %7.4g %7.0f %7.4g\n", N1, sqrt (N2/Npix - N1*N1), + min, max, median, Npix, Npix*N1); + } + + set_variable ("MIN", min); + set_variable ("MAX", max); + set_variable ("MEDIAN", median); + set_variable ("MEAN", N1); + set_variable ("MODE", mode); + set_variable ("TOTAL", N1*Npix); + set_variable ("NPIX", Npix); + set_variable ("SIGMA", sqrt (N2/Npix - N1*N1)); + + return (TRUE); +} + Index: trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/init.c (revision 21064) +++ trunk/Ohana/src/opihi/cmd.data/init.c (revision 21508) @@ -3,4 +3,5 @@ int accum PROTO((int, char **)); int applyfit PROTO((int, char **)); +int applyfit1d PROTO((int, char **)); int applyfit2d PROTO((int, char **)); int box PROTO((int, char **)); @@ -11,4 +12,5 @@ int clear PROTO((int, char **)); int clip PROTO((int, char **)); +int imclip PROTO((int, char **)); int close_device PROTO((int, char **)); int concat PROTO((int, char **)); @@ -29,6 +31,7 @@ int fft1d PROTO((int, char **)); int fft2d PROTO((int, char **)); +int fit PROTO((int, char **)); +int fit1d PROTO((int, char **)); int fit2d PROTO((int, char **)); -int fit PROTO((int, char **)); int gaussjordan PROTO((int, char **)); int gaussdeviate PROTO((int, char **)); @@ -92,4 +95,5 @@ int spline_construct_cmd PROTO((int, char **)); int stats PROTO((int, char **)); +int imstats PROTO((int, char **)); int style PROTO((int, char **)); int subraster PROTO((int, char **)); @@ -113,5 +117,5 @@ int vload PROTO((int, char **)); int vzload PROTO((int, char **)); -int vstat PROTO((int, char **)); +int vstats PROTO((int, char **)); int vroll PROTO((int, char **)); int vpop PROTO((int, char **)); @@ -124,17 +128,20 @@ static Command cmds[] = { {1, "accum", accum, "accumulate vector values in another vector"}, - {1, "applyfit", applyfit, "apply fit to new vector"}, + {1, "applyfit", applyfit1d, "apply 1-d fit to new vector"}, + {1, "applyfit1d", applyfit1d, "apply 1-d fit to new vector"}, {1, "applyfit2d", applyfit2d, "apply 2-d fit to new vector"}, + {1, "book", book_command, "commands to manipulate book/page/word data"}, {1, "box", box, "draw a box on the plot"}, - {1, "book", book_command, "commands to manipulate book/page/word data"}, + {1, "buffers", list_buffers, "list the currently allocated buffers (images)"}, + {1, "images", list_buffers, "list the currently allocated images"}, {1, "center", center, "center image on coords"}, - {1, "parity", parity, "set image parity"}, {1, "circstats", circstats, "circular statistics"}, {1, "clear", clear, "erase plot"}, - {1, "clip", clip, "clip values in a buffer to be within a range"}, + {1, "clip", imclip, "clip values in an image to be within a range"}, {1, "close", close_device, "close the current display device"}, - {1, "concat", concat, "reduce vector dimension"}, + {1, "concat", concat, "append values to the end of a vector"}, {1, "contour", contour, "create contour from image"}, {1, "create", create, "create a new vector"}, + {1, "vcreate", create, "create a new vector"}, {1, "cumulative", cumulative, "build a cumulative histogram from a specific histogram"}, {1, "cursor", cursor, "get coords from cursor"}, @@ -143,81 +150,87 @@ {1, "dbconnect", dbconnect, "setup mysql db connection"}, {1, "dbselect", dbselect, "extract vectors from mysql database table"}, - {1, "delete", delete, "delete vectors or matrices"}, + {1, "delete", delete, "delete vectors or images"}, {1, "device", device, "set / get current graphics device"}, - {1, "dimendown", dimendown, "convert matrix to vector"}, - {1, "dimenup", dimenup, "convert vector to matrix"}, - {1, "erase", erase, "erase objects on an overlay"}, - {1, "extract", extract, "extract a portion of a buffer into another buffer"}, - {1, "fft1d", fft1d, "fft on the pixel-stream in an image"}, + {1, "dimendown", dimendown, "convert image to vector"}, + {1, "dimenup", dimenup, "convert vector to image"}, + {1, "dot", dot, "plot a single point"}, + {1, "erase", erase, "erase objects on an image overlay"}, + {1, "extract", extract, "extract a portion of a image into another image"}, + {1, "fft1d", fft1d, "fft on a vector"}, {1, "fft2d", fft2d, "fft on an image"}, - {1, "fit", fit, "fit polynomial to vector pair"}, + {1, "fit", fit1d, "fit polynomial to vector pair"}, + {1, "fit1d", fit1d, "fit polynomial to vector pair"}, {1, "fit2d", fit2d, "fit 2-d polynomial to vector triplet"}, - {1, "gaussj", gaussjordan, "solve Ax = B (N-D)"}, {1, "gaussdev", gaussdeviate, "generate a gaussian deviate vector"}, {1, "gaussint", gaussintegral, "return the integrated gaussian vector"}, - {1, "grid", grid, "plot cartesian grid"}, - {1, "gridify", gridify, "convert vector triplet to buffer"}, + {1, "gaussj", gaussjordan, "solve Ax = B (N-Dimensional)"}, + {1, "grid", grid, "plot cartesian grid on graph"}, + {1, "gridify", gridify, "convert vector triplet to image (same as vgrid?)"}, {1, "grow", grow, "grow a mask"}, - {1, "ungridify", ungridify, "convert buffer region to vector triplet"}, - {1, "header", header, "print buffer header"}, + {1, "header", header, "print image header"}, {1, "histogram", histogram, "generate histogram from vector"}, + {1, "imbin", rebin, "rebin image data by factor of N"}, + {1, "imclip", imclip, "clip values in an image to be within a range"}, {1, "imcut", imcut, "linear image cut between arbitrary coords"}, - {1, "imhist", imhist, "histogram of an image region"}, + {1, "imhistogram", imhist, "histogram of an image region"}, {1, "imsmooth", imsmooth, "circular gaussian smoothing"}, + {1, "imstats", imstats, "statistics on a portion of an image"}, {1, "integrate", integrate, "integrate a vector"}, {1, "interpolate", interpolate, "interpolate between vector pairs"}, - {1, "jpeg", jpeg, "write text line on graph"}, - {1, "png", jpeg, "write text line on graph"}, - {1, "ppm", jpeg, "write text line on graph"}, + {1, "jpeg", jpeg, "convert display image to JPEG"}, {1, "kern", kern, "convolve with 3x3 kernel"}, - {1, "keyword", keyword, "extract a FITS keyword from buffer header"}, + {1, "keyword", keyword, "extract a FITS keyword from image header"}, {1, "labels", labels, "define labels for plot"}, {1, "limits", limits, "define plot limits"}, - {1, "line", line, "plot line"}, - {1, "buffers", list_buffers, "list the currently allocated buffers"}, - {1, "vectors", list_vectors, "list vectors"}, - {1, "load", load, "load an SAOimage style overlay"}, + {1, "line", line, "plot a line"}, + {1, "load", load, "load an SAOimage/DS9 style overlay"}, {1, "lookup", lookup, "convert vector via lookup table (vector pair)"}, - {1, "mkrgb", mkrgb, "convert 3 images to rgb jpeg"}, - {1, "mcreate", mcreate, "create a matrix"}, + {1, "mcreate", mcreate, "create an image"}, + {1, "imcreate", mcreate, "create an image"}, {1, "medacc", medacc, "accumulate vector values in another vector"}, - {1, "mget", mget, "extract a vector from a matrix"}, + {1, "mget", mget, "extract a vector from an image"}, + {1, "imget", mget, "extract a vector from an image"}, {1, "minterp", minterp, "interpolate image pixels"}, - {1, "mset", mset, "insert a vector in a matrix"}, + {1, "iminterp", minterp, "interpolate image pixels"}, + {1, "mkrgb", mkrgb, "convert 3 images to rgb jpeg (use Kapa for better control)"}, + {1, "imset", mset, "insert a vector in an image"}, + {1, "parity", parity, "set image parity"}, {1, "peak", peak, "find vector peak in range"}, {1, "periodogram", periodogram, "measure periods in unevenly sampled data"}, {1, "plot", plot, "plot a pair of vectors"}, - {1, "dot", dot, "plot a single point"}, - {1, "point", point, "load overlay with single point"}, - {1, "ps", ps, "define labels for plot"}, + {1, "png", jpeg, "convert display graphic to PNG"}, + {1, "point", point, "load image overlay with a single point"}, + {1, "ppm", jpeg, "convert display graphic to PPM"}, + {1, "ps", ps, "convert display to PostScript"}, + {1, "queuedelete", queuedelete, "delete a queue"}, + {1, "queuedrop", queuedrop, "drop values from queue matching a key"}, + {1, "queueinit", queueinit, "create an empty queue"}, + {1, "queuelist", queuelist, "list defined queues"}, + {1, "queueload", queueload, "load queue from command"}, {1, "queuepop", queuepop, "pop value from queue to variable"}, - {1, "queuedrop", queuedrop, "drop values from queue matching a key"}, {1, "queueprint", queueprint, "print the contents of a queue"}, {1, "queuepush", queuepush, "push value onto queue"}, - {1, "queueinit", queueinit, "create an empty queue"}, - {1, "queuedelete", queuedelete, "delete a queue"}, - {1, "queuelist", queuelist, "list defined queues"}, - {1, "queueload", queueload, "load queue from command"}, + {1, "queuesize", queuesize, "show queue size"}, {1, "queuesubstr", queuesubstr, "bulk replace strings in queue"}, - {1, "queuesize", queuesize, "show queue size"}, {1, "rd", rd, "load fits image"}, {1, "rdseg", rdseg, "read a segment of an image from a file"}, {1, "read", read_vectors, "read vectors from datafile"}, - {1, "rebin", rebin, "rebin data by factor of N"}, + {1, "rebin", rebin, "rebin image data by factor of N"}, {1, "resize", resize, "set graphics/image window size"}, {1, "roll", roll, "roll image to new start point"}, {1, "rotate", rotate, "rotate image"}, - {1, "save", save, "save an SAOimage style overlay"}, + {1, "save", save, "save an SAOimage style image overlay"}, {1, "section", section, "define section of graph"}, - {1, "set", set, "vector math"}, + {1, "select", vect_select, "selective vector assignment"}, + {1, "set", set, "image and vector math"}, {1, "shift", shift, "shift data in an image"}, {1, "sort", sort_vectors, "sort list of vectors"}, {1, "spline.apply", spline_apply_cmd, "apply spline fit to generate an image"}, {1, "spline.const", spline_construct_cmd, "create spline 2nd deriv. terms"}, - {1, "stats", stats, "give statistics on a portion of a buffer"}, + {1, "stats", imstats, "statistics on a portion of an image"}, {1, "style", style, "set the style for graph plots"}, {1, "subraster", subraster, "subraster of fits image"}, - {1, "subset", subset, "expand vector dimension"}, - {1, "svd", svd, "singular value decomposition of a matrix"}, + {1, "subset", subset, "generate a vector from a portion of another vector"}, + {1, "svd", svd, "singular value decomposition of a matrix (image)"}, {1, "swapbytes", swapbytes, "byte swap thing"}, {1, "textline", textline, "write text line on graph"}, @@ -225,23 +238,25 @@ {1, "tvchannel", tvchannel, "set the current tv channel"}, {1, "tvcolors", tvcolors, "set the tv colormap"}, - {1, "tvcontour", tvcontour, "send contour to overlay"}, - {1, "tvgrid", tvgrid, "wait until return is typed"}, + {1, "tvcontour", tvcontour, "send contour to image overlay"}, + {1, "tvgrid", tvgrid, "RA/DEC grid on an image"}, + {1, "ungridify", ungridify, "convert image region to vector triplet"}, {1, "uniq", uniq, "create a uniq vector subset from a vector"}, {1, "unsign", unsign, "toggle the UNSIGN status"}, - {1, "vbin", vbin, "bin values in a vector to be within a range"}, + {1, "vbin", vbin, "rebin vector data by a fector of N"}, {1, "vclip", vclip, "clip values in a vector to be within a range"}, - {1, "select", vect_select, "selective vector assignment"}, - {1, "vgauss", vgauss, ""}, - {1, "vmaxwell", vmaxwell, ""}, - {1, "vgrid", vgrid, ""}, - {1, "vload", vload, "load vectors on Kii"}, - {1, "vzload", vzload, "load vectors on Kii"}, - {1, "vstat", vstat, "get info from imreg database"}, - {1, "vsmooth", vsmooth, "gaussian smooth of a vector"}, + {1, "vectors", list_vectors, "list vectors"}, + {1, "vgauss", vgauss, "fit a Gaussian to a vector"}, + {1, "vgrid", vgrid, "generate an image from a triplet of vectors"}, + {1, "vhistogram", histogram, "generate histogram from vector"}, + {1, "vload", vload, "load vectors as overlay on image display"}, + {1, "vmaxwell", vmaxwell, "fit a Maxwellian to a vector"}, + {1, "vpop", vpop, "remove first element of a vector"}, {1, "vroll", vroll, "roll vector elements"}, - {1, "vpop", vpop, "remove first element"}, + {1, "vsmooth", vsmooth, "Gaussian smooth of a vector"}, + {1, "vstats", vstats, "statistics on a vector"}, + {1, "vzload", vzload, "load vectors as overlay on image display (scaled points)"}, {1, "wd", wd, "write an image to a file"}, {1, "write", write_vectors, "write vectors to datafile"}, - {1, "zap", zap, "delete pixels"}, + {1, "zap", zap, "assign values to pixel regions"}, {1, "zplot", zplot, "plot x y with size scaled by z"}, }; Index: trunk/Ohana/src/opihi/cmd.data/stats.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/stats.c (revision 21064) +++ (revision ) @@ -1,125 +1,0 @@ -# include "data.h" - -int stats (int argc, char **argv) { - - int i, j, Nmode, Imode; - double Npix, N1, N2, max, min, range, median, mode, IgnoreValue; - float *V; - int sx, sy, nx, ny, *hist, Nhist, bin; - int Ignore, Quiet, N; - Buffer *buf; - - IgnoreValue = 0; - Ignore = FALSE; - if ((N = get_argument (argc, argv, "-ignore"))) { - Ignore = TRUE; - remove_argument (N, &argc, argv); - IgnoreValue = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - 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 != 2) && (argc != 6)) { - gprint (GP_ERR, "USAGE: stats sx sy nx ny\n"); - gprint (GP_ERR, "OR: stats \n"); - return (FALSE); - } - - if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - - if (argc == 6) { - sx = strcmp (argv[2], "-") ? atof (argv[2]) : 0; - sy = strcmp (argv[3], "-") ? atof (argv[3]) : 0; - nx = strcmp (argv[4], "-") ? atof (argv[4]) : buf[0].matrix.Naxis[0]; - ny = strcmp (argv[5], "-") ? atof (argv[5]) : buf[0].matrix.Naxis[1]; - } else { - sx = 0; - sy = 0; - nx = buf[0].matrix.Naxis[0]; - ny = buf[0].matrix.Naxis[1]; - } - - Npix = N1 = N2 = 0; - if ((sx < 0) || (sy < 0) || - (sx+nx > buf[0].matrix.Naxis[0]) || - (sy+ny > buf[0].matrix.Naxis[1])) { - gprint (GP_ERR, "region out of range\n"); - return (FALSE); - } - - max = -1e32; - min = +1e32; - for (j = sy; j < sy + ny; j++) { - V = (float *)(buf[0].matrix.buffer) + j*buf[0].matrix.Naxis[0] + sx; - for (i = 0; i < nx; i++, V++) { - if (Ignore && (fabs (*V - IgnoreValue) < 1e-8)) continue; - if (isnan(*V)) continue; - if (isinf(*V)) continue; - N1 += *V; - N2 += (*V)*(*V); - Npix += 1.0; - max = MAX (max, *V); - min = MIN (min, *V); - } - } - N1 = N1 / Npix; - -/* calculate mode, median */ - median = mode = 0.5*(max + min); - if ((max - min) != 0) { - range = 0xffff / (max - min); - ALLOCATE (hist, int, 0x10000); - bzero (hist, 0x10000*sizeof(int)); - for (j = sy; j < sy + ny; j++) { - V = (float *)(buf[0].matrix.buffer) + j*buf[0].matrix.Naxis[0] + sx; - for (i = 0; i < nx; i++, V++) { - if (Ignore && (fabs (*V - IgnoreValue) < 1e-8)) continue; - if (isnan(*V)) continue; - if (isinf(*V)) continue; - bin = MIN (MAX (0, (*V - min) * range), 0xffff); - hist[bin] ++; - } - } - Nhist = 0; - for (i = 0; (i < 0xffff) && (Nhist < 0.5*Npix); i++) - Nhist += hist[i]; - median = i / range + min; - Nmode = hist[0]; - Imode = 0; - for (i = 1; i < 0x10000; i++) { - if (hist[i] > Nmode) { - Nmode = hist[i]; - Imode = i; - } - } - mode = Imode / range + min; - free (hist); - } - - if (!Quiet) { - gprint (GP_LOG, " mean stdev min max median Npix Total\n"); - gprint (GP_LOG, "%7.4g %7.4g %7.4g %7.4g %7.4g %7.0f %7.4g\n", N1, sqrt (N2/Npix - N1*N1), - min, max, median, Npix, Npix*N1); - } - - set_variable ("MIN", min); - set_variable ("MAX", max); - set_variable ("MEDIAN", median); - set_variable ("MEAN", N1); - set_variable ("MODE", mode); - set_variable ("TOTAL", N1*Npix); - set_variable ("NPIX", Npix); - set_variable ("SIGMA", sqrt (N2/Npix - N1*N1)); - - return (TRUE); -} - Index: trunk/Ohana/src/opihi/cmd.data/vstat.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vstat.c (revision 21064) +++ (revision ) @@ -1,139 +1,0 @@ -# include "data.h" - -int vstat (int argc, char **argv) { - - int i, N; - double max, min, sum, var, dvar, mean, stdev; - float IgnoreValue; - int Ignore, Quiet; - - int *Nval, bin, Nmode, Nmed; - double dx, mode, median; - Vector *vec; - - IgnoreValue = 0; - Ignore = FALSE; - if ((N = get_argument (argc, argv, "-ignore"))) { - Ignore = TRUE; - remove_argument (N, &argc, argv); - IgnoreValue = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - 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 != 2) { - gprint (GP_ERR, "USAGE: vstat (vector)\n"); - return (FALSE); - } - - if ((vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - - /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */ - - /* calculate max, min, mean, sum, npix */ - max = -HUGE_VAL; - min = HUGE_VAL; - sum = N = 0; - if (vec[0].type == OPIHI_FLT) { - opihi_flt *X = vec[0].elements.Flt; - for (i = 0; i < vec[0].Nelements; i++, X++) { - if (!finite (*X)) continue; - if (Ignore && (*X == IgnoreValue)) continue; - max = MAX (*X, max); - min = MIN (*X, min); - sum += *X; - N++; - } - } else { - opihi_int *X = vec[0].elements.Int; - for (i = 0; i < vec[0].Nelements; i++, X++) { - if (!finite (*X)) continue; - if (Ignore && (*X == IgnoreValue)) continue; - max = MAX (*X, max); - min = MIN (*X, min); - sum += *X; - N++; - } - } - mean = sum / N; - - /* calculate median and mode with resolution of (max - min) / 1000 */ - dx = (max - min) / 1000; - if (dx == 0) { - median = mode = min; - stdev = 0.0; - goto skip; - } - - ALLOCATE (Nval, int, 1002); - bzero (Nval, 1000*sizeof(int)); - var = 0; - if (vec[0].type == OPIHI_FLT) { - opihi_flt *X = vec[0].elements.Flt; - for (i = 0; i < vec[0].Nelements; i++, X++) { - if (!finite (*X)) continue; - if (Ignore && (*X == IgnoreValue)) continue; - bin = MAX (0, MIN (1000, (*X - min) / dx)); - Nval[bin] ++; - dvar = (*X - mean); - var += dvar*dvar; - } - } else { - opihi_int *X = vec[0].elements.Int; - for (i = 0; i < vec[0].Nelements; i++, X++) { - if (!finite (*X)) continue; - if (Ignore && (*X == IgnoreValue)) continue; - bin = MAX (0, MIN (1000, (*X - min) / dx)); - Nval[bin] ++; - dvar = (*X - mean); - var += dvar*dvar; - } - } - stdev = sqrt (var / N); - - Nmode = 0; - mode = Nval[Nmode]; - median = 0; - Nmed = -1; - for (i = 0; i < 1001; i++) { - if (Nmed == -1) { - median += Nval[i]; - if (median >= N / 2.0) { - Nmed = i; - median = i * dx + min; - } - } - if (mode < Nval[i]) { - Nmode = i; - mode = Nval[Nmode]; - } - } - mode = Nmode * dx + min; - free (Nval); - -skip: - if (!Quiet) { - gprint (GP_ERR, "mean: %g, stdev: %g, min: %g, max: %g, median: %g, mode: %g, Npts: %d\n", - mean, stdev, min, max, median, mode, N); - } - - set_variable ("MIN", min); - set_variable ("MAX", max); - set_variable ("MEDIAN", median); - set_variable ("MEAN", mean); - set_variable ("MODE", mode); - set_variable ("TOTAL", sum); - set_int_variable ("NPIX", N); - set_variable ("SIGMA", stdev); - return (TRUE); -} - Index: trunk/Ohana/src/opihi/cmd.data/vstats.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vstats.c (revision 21508) +++ trunk/Ohana/src/opihi/cmd.data/vstats.c (revision 21508) @@ -0,0 +1,139 @@ +# include "data.h" + +int vstats (int argc, char **argv) { + + int i, N; + double max, min, sum, var, dvar, mean, stdev; + float IgnoreValue; + int Ignore, Quiet; + + int *Nval, bin, Nmode, Nmed; + double dx, mode, median; + Vector *vec; + + IgnoreValue = 0; + Ignore = FALSE; + if ((N = get_argument (argc, argv, "-ignore"))) { + Ignore = TRUE; + remove_argument (N, &argc, argv); + IgnoreValue = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + 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 != 2) { + gprint (GP_ERR, "USAGE: vstat (vector)\n"); + return (FALSE); + } + + if ((vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */ + + /* calculate max, min, mean, sum, npix */ + max = -HUGE_VAL; + min = HUGE_VAL; + sum = N = 0; + if (vec[0].type == OPIHI_FLT) { + opihi_flt *X = vec[0].elements.Flt; + for (i = 0; i < vec[0].Nelements; i++, X++) { + if (!finite (*X)) continue; + if (Ignore && (*X == IgnoreValue)) continue; + max = MAX (*X, max); + min = MIN (*X, min); + sum += *X; + N++; + } + } else { + opihi_int *X = vec[0].elements.Int; + for (i = 0; i < vec[0].Nelements; i++, X++) { + if (!finite (*X)) continue; + if (Ignore && (*X == IgnoreValue)) continue; + max = MAX (*X, max); + min = MIN (*X, min); + sum += *X; + N++; + } + } + mean = sum / N; + + /* calculate median and mode with resolution of (max - min) / 1000 */ + dx = (max - min) / 1000; + if (dx == 0) { + median = mode = min; + stdev = 0.0; + goto skip; + } + + ALLOCATE (Nval, int, 1002); + bzero (Nval, 1000*sizeof(int)); + var = 0; + if (vec[0].type == OPIHI_FLT) { + opihi_flt *X = vec[0].elements.Flt; + for (i = 0; i < vec[0].Nelements; i++, X++) { + if (!finite (*X)) continue; + if (Ignore && (*X == IgnoreValue)) continue; + bin = MAX (0, MIN (1000, (*X - min) / dx)); + Nval[bin] ++; + dvar = (*X - mean); + var += dvar*dvar; + } + } else { + opihi_int *X = vec[0].elements.Int; + for (i = 0; i < vec[0].Nelements; i++, X++) { + if (!finite (*X)) continue; + if (Ignore && (*X == IgnoreValue)) continue; + bin = MAX (0, MIN (1000, (*X - min) / dx)); + Nval[bin] ++; + dvar = (*X - mean); + var += dvar*dvar; + } + } + stdev = sqrt (var / N); + + Nmode = 0; + mode = Nval[Nmode]; + median = 0; + Nmed = -1; + for (i = 0; i < 1001; i++) { + if (Nmed == -1) { + median += Nval[i]; + if (median >= N / 2.0) { + Nmed = i; + median = i * dx + min; + } + } + if (mode < Nval[i]) { + Nmode = i; + mode = Nval[Nmode]; + } + } + mode = Nmode * dx + min; + free (Nval); + +skip: + if (!Quiet) { + gprint (GP_ERR, "mean: %g, stdev: %g, min: %g, max: %g, median: %g, mode: %g, Npts: %d\n", + mean, stdev, min, max, median, mode, N); + } + + set_variable ("MIN", min); + set_variable ("MAX", max); + set_variable ("MEDIAN", median); + set_variable ("MEAN", mean); + set_variable ("MODE", mode); + set_variable ("TOTAL", sum); + set_int_variable ("NPIX", N); + set_variable ("SIGMA", stdev); + return (TRUE); +} +