|
|
1.1 root 1: #include <stdio.h>
2: #include <math.h>
3: #include "view2d.h"
4: #define MAXPTS 4096
5:
6: int NX = 16, NY = 16;
7: char *progname;
8:
9: typedef struct Vector
10: {
11: float x, y, z;
12: } Vector;
13:
14: typedef struct Frame
15: {
16: Vector data[MAXPTS];
17: double time;
18: int npts;
19: } Frame;
20: Frame frame[40];
21: int interp = 1;
22:
23: float grid[200][200];
24: float xmin, xmax, ymin, ymax;
25: float xstep, ystep;
26: short *data;
27: int verbose;
28:
29: main(argc, argv)
30: char **argv;
31: {
32: register i, j;
33: int nframes;
34: float x;
35:
36: verbose = 0;
37: progname = *argv;
38: for(argc--, argv++; **argv == '-'; argv++)
39: switch(argv[0][1])
40: {
41: case 'n':
42: NY = -1;
43: sscanf(&argv[0][2], "%d,%d", &NX, &NY);
44: if(NY < 0) NY = NX;
45: break;
46: case 'c': /* piecewise constant */
47: interp = 2;
48: break;
49: case 's': /* scatterplot on black background */
50: interp = 3;
51: break;
52: case 'i':
53: interp = 0;
54: break;
55: case 'v':
56: verbose = 1;
57: break;
58: default:
59: error("usage: vdata -c -s -i -v -nNX,NY file");
60: }
61: if(*argv)
62: if(freopen(*argv, "r", stdin)==NULL)
63: error("freopen(%s) failed",*argv);
64: if((NX>200)||(NY>200)) error("NX, NY must be <= 200");
65: data = (short *)Malloc(NX*NY*sizeof(short));
66:
67: if(interp)
68: {
69: xmax = ymax = -HUGE;
70: xmin = ymin = HUGE;
71: }
72: for(nframes = 0; readf(&frame[nframes]); nframes++){
73: if(nframes>38) error("too many frames");
74: }
75: if(nframes<=0) error("empty input");
76: if(verbose)fprintf(stderr, "nx=%d ny=%d nframes=%d\n", NX, NY, nframes);
77:
78: if(interp)
79: {
80: xstep = (xmax-xmin)/(NX-1);
81: ystep = (ymax-ymin)/(NY-1);
82: if(verbose)fprintf(stderr, "%.5g<=x<=%.5g %.5g<=y<=%.5g\n",
83: xmin, xmax, ymin, ymax);
84: }
85: else
86: {
87: xmin = 0; xstep = 1; xmax = xmin + (NX-1);
88: ymin = 0; ystep = 1; ymax = ymin + (NY-1);
89: }
90: for(i = 0; i < nframes; i++){
91: fit(&frame[i]);
92: out(frame[i].time,frame[i].npts);
93: }
94: exit(0);
95: }
96:
97: readf(s)
98: Frame *s;
99: {
100: register Vector *v;
101: register i;
102: float f;
103:
104: if(scanf("%d %E", &s->npts, &s->time) != 2) return(0);
105: if(verbose) fprintf(stderr,"%d %g\n",s->npts,s->time);
106: if(s->npts>MAXPTS) error("too many data points");
107: if(interp)
108: for(v = s->data, i = 0; i < s->npts; i++, v++)
109: {
110: scanf("%f %f %f", &v->x, &v->y, &v->z);
111: if(v->x < xmin) xmin = v->x;
112: if(v->x > xmax) xmax = v->x;
113: if(v->y < ymin) ymin = v->y;
114: if(v->y > ymax) ymax = v->y;
115: }
116: else
117: for(v = s->data, i = 0; i < s->npts; i++, v++)
118: {
119: scanf("%f", &v->z);
120: }
121: return(1);
122: }
123:
124: fit(s)
125: Frame *s;
126: {
127: register Vector *v;
128: Vector h;
129: float g, x, y;
130: register i, j, k;
131: Vector *c1, *c2, *c3;
132: float d, d1, d2, d3;
133:
134: if( (interp==1) || (interp==2) ) {
135: for(j = 0; j < NY; j++)
136: for(i = 0; i < NX; i++) {
137: x = xmin + i*xstep;
138: y = ymin + j*ystep;
139: g = grid[i][j];
140: d1 = d2 = d3 = HUGE;
141: for(v = s->data, k = 0; k < s->npts; k++, v++) {
142: d = (v->x - x)*(v->x - x) +
143: (v->y - y)*(v->y - y);
144: if(d < d1) {
145: c3 = c2; d3 = d2;
146: c2 = c1; d2 = d1;
147: c1 = v; d1 = d;
148: } else if(d < d2) {
149: c3 = c2; d3 = d2;
150: c2 = v; d2 = d;
151: } else if(d < d3) {
152: c3 = v; d3 = d;
153: }
154: }
155: /* c1 c2 c3 three closest points */
156: if(interp==1) {
157: h.x = x; h.y = y; h.z = g;
158: interp_(&h, c1, c2, c3);
159: grid[i][j] = h.z;
160: } else if(interp==2)
161: grid[i][j] = c1->z;
162: }
163: } else if( interp==3 ) {
164: for(j = 0; j < NY; j++)
165: for(i = 0; i < NX; i++)
166: grid[i][j] = -1e26;
167: for(v = s->data, k = 0; k < s->npts; k++, v++) {
168: i = (NX-1)*(v->x-xmin)/(xmax-xmin);
169: j = (NY-1)*(v->y-ymin)/(ymax-ymin);
170: if((i<0)||(i>=NX)||(j<0)||(j>=NY))
171: error("can't: i,j=%d,%d\n",i,j);
172: grid[i][j] = v->z;
173: }
174: } else { /* option -i */
175: for(v = s->data, j = 0; j < NY; j++)
176: for(i = 0; i < NX; i++, v++)
177: grid[i][j] = v->z;
178: }
179: }
180:
181: out(tim,npt)
182: double tim;
183: int npt;
184: {
185: register short *d, i, j;
186: register float zmax, zmin, zrange, f, v2;
187: short u, v;
188: extern double pow2();
189:
190: zmin = HUGE;
191: zmax = -HUGE;
192: for(j = 0; j < NY; j++)
193: for(i = 0; i < NX; i++){
194: if(grid[i][j] < -1e25 ) continue;
195: if(grid[i][j] < zmin) zmin = grid[i][j];
196: if(grid[i][j] > zmax) zmax = grid[i][j];
197: }
198:
199: if(verbose) fprintf(stderr, "time=%g npts=%d %g<=z<=%g\n",
200: tim, npt, zmin, zmax);
201: zmax = (zmax<0)?0:zmax;
202: zmin = (zmin>0)?0:zmin;
203: zrange = zmax - zmin;
204: if(zrange<=0)
205: v=0, u=0, v2=1, f=0;
206: else{
207: f = 2*BIG-2;
208: v = log(zrange)/log(2.)-14;
209: v2 = pow2(v);
210: while( 2*zrange < f*v2 ){ v--; v2 /= 2; }
211: f = -(zmin+zmax)/(2*v2);
212: u = f;
213: if( (f<-BIG)||(f>BIG) ){
214: if(verbose)fprintf(stderr,"zmin=%g zmax=%g u=%d v=%d\n",
215: zmin,zmax,u,v);
216: error("u out of bounds f=%g",f);
217: }
218: }
219: d = data;
220: for(j = 0; j < NY; j++)
221: for(i = 0; i < NX; i++)
222: if(grid[i][j] < -1e25){ *d++ = -BIG-1; }
223: else{ *d++ = f + grid[i][j]/v2; }
224: view2d(1, NX, NY, tim, u, v, 0, 0, 0, data);
225: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.