V10/cmd/view2d/vdata.c

#include	<stdio.h>
#include	<math.h>
#include	"view2d.h"
#define MAXPTS 4096

int NX = 16, NY = 16;
char *progname;

typedef struct Vector
{
	float x, y, z;
} Vector;

typedef struct Frame
{
	Vector data[MAXPTS];
	double time;
        int npts;
} Frame;
Frame frame[40];
int interp = 1;

float grid[200][200];
float xmin, xmax, ymin, ymax;
float xstep, ystep;
short *data;
int verbose;

main(argc, argv)
	char **argv;
{
	register i, j;
	int nframes;
	float x;

	verbose = 0;
	progname = *argv;
	for(argc--, argv++; **argv == '-'; argv++)
		switch(argv[0][1])
		{
		case 'n':
			NY = -1;
			sscanf(&argv[0][2], "%d,%d", &NX, &NY);
			if(NY < 0) NY = NX;
			break;
		case 'c':   /* piecewise constant */
			interp = 2;
			break;
		case 's':   /* scatterplot on black background */
			interp = 3;
			break;
		case 'i':
			interp = 0;
			break;
		case 'v':
			verbose = 1;
			break;
		default:
			error("usage: vdata -c -s -i -v -nNX,NY file");
		}
	if(*argv)
		if(freopen(*argv, "r", stdin)==NULL)
			error("freopen(%s) failed",*argv);
	if((NX>200)||(NY>200)) error("NX, NY must be <= 200");
	data = (short *)Malloc(NX*NY*sizeof(short));

	if(interp)
	{
		xmax = ymax = -HUGE;
		xmin = ymin = HUGE;
	}
	for(nframes = 0; readf(&frame[nframes]); nframes++){
	  if(nframes>38) error("too many frames");
        }
	if(nframes<=0) error("empty input");
	if(verbose)fprintf(stderr, "nx=%d ny=%d nframes=%d\n", NX, NY, nframes);

	if(interp)
	{
		xstep = (xmax-xmin)/(NX-1);
		ystep = (ymax-ymin)/(NY-1);
		if(verbose)fprintf(stderr, "%.5g<=x<=%.5g  %.5g<=y<=%.5g\n",
			xmin, xmax, ymin, ymax);
	}
	else
	{
		xmin = 0; xstep = 1; xmax = xmin + (NX-1);
		ymin = 0; ystep = 1; ymax = ymin + (NY-1);
	}
	for(i = 0; i < nframes; i++){
		fit(&frame[i]);
		out(frame[i].time,frame[i].npts);
	}
	exit(0);
}

readf(s)
	Frame *s;
{
	register Vector *v;
	register i;
	float f;

	if(scanf("%d %E", &s->npts, &s->time) != 2) return(0);
	if(verbose) fprintf(stderr,"%d %g\n",s->npts,s->time);
	if(s->npts>MAXPTS) error("too many data points");
	if(interp)
		for(v = s->data, i = 0; i < s->npts; i++, v++)
		{
			scanf("%f %f %f", &v->x, &v->y, &v->z);
			if(v->x < xmin) xmin = v->x;
			if(v->x > xmax) xmax = v->x;
			if(v->y < ymin) ymin = v->y;
			if(v->y > ymax) ymax = v->y;
		}
	else
		for(v = s->data, i = 0; i < s->npts; i++, v++)
		{
			scanf("%f", &v->z);
		}
	return(1);
}

fit(s)
	Frame *s;
{
	register Vector *v;
	Vector h;
	float g, x, y;
	register i, j, k;
	Vector *c1, *c2, *c3;
	float d, d1, d2, d3;

	if( (interp==1) || (interp==2) ) {
		for(j = 0; j < NY; j++)
			for(i = 0; i < NX; i++) {
				x = xmin + i*xstep;
				y = ymin + j*ystep;
				g = grid[i][j];
				d1 = d2 = d3 = HUGE;
				for(v = s->data, k = 0; k < s->npts; k++, v++) {
					d = (v->x - x)*(v->x - x) +
						(v->y - y)*(v->y - y);
					if(d < d1) {
						c3 = c2; d3 = d2;
						c2 = c1; d2 = d1;
						c1 = v; d1 = d;
					} else if(d < d2) {
						c3 = c2; d3 = d2;
						c2 = v; d2 = d;
					} else if(d < d3) {
						c3 = v; d3 = d;
					}
				}
				/* c1 c2 c3 three closest points */
				if(interp==1) {
					h.x = x; h.y = y; h.z = g;
					interp_(&h, c1, c2, c3);
					grid[i][j] = h.z;
				} else if(interp==2)
					grid[i][j] = c1->z;
			}
	} else if( interp==3 ) {
		for(j = 0; j < NY; j++)
			for(i = 0; i < NX; i++)
				grid[i][j] = -1e26;
		for(v = s->data, k = 0; k < s->npts; k++, v++) {
			i = (NX-1)*(v->x-xmin)/(xmax-xmin);
			j = (NY-1)*(v->y-ymin)/(ymax-ymin);
			if((i<0)||(i>=NX)||(j<0)||(j>=NY))
				error("can't: i,j=%d,%d\n",i,j);
			grid[i][j] = v->z;
		}
	} else {  /* option -i */
		for(v = s->data, j = 0; j < NY; j++)
			for(i = 0; i < NX; i++, v++)
				grid[i][j] = v->z;
	}
}

out(tim,npt)
  double tim;
  int npt;
{
	register short *d, i, j;
	register float zmax, zmin, zrange, f, v2;
	short u, v;
	extern double pow2();

	zmin = HUGE;
	zmax = -HUGE;
	for(j = 0; j < NY; j++)
		for(i = 0; i < NX; i++){
			if(grid[i][j] < -1e25 ) continue;
			if(grid[i][j] < zmin) zmin = grid[i][j];
			if(grid[i][j] > zmax) zmax = grid[i][j];
		}

	if(verbose) fprintf(stderr, "time=%g npts=%d %g<=z<=%g\n",
		 tim, npt, zmin, zmax);
	zmax = (zmax<0)?0:zmax;
	zmin = (zmin>0)?0:zmin;
	zrange = zmax - zmin;
	if(zrange<=0)
		v=0, u=0, v2=1, f=0;
	else{
		f = 2*BIG-2;
		v = log(zrange)/log(2.)-14;
		v2 = pow2(v);
		while( 2*zrange < f*v2 ){ v--; v2 /= 2; }
		f = -(zmin+zmax)/(2*v2);
		u = f;
		if( (f<-BIG)||(f>BIG) ){
			if(verbose)fprintf(stderr,"zmin=%g zmax=%g u=%d v=%d\n",
				zmin,zmax,u,v);
			error("u out of bounds f=%g",f);
		}
	}
	d = data;
	for(j = 0; j < NY; j++)
		for(i = 0; i < NX; i++)
			if(grid[i][j] < -1e25){ *d++ = -BIG-1; }
			 else{ *d++ = f + grid[i][j]/v2; }
	view2d(1, NX, NY, tim, u, v, 0, 0, 0, data);
}