![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to regrid.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv8dc / cmd / view2d|
Return to regrid.c CVS log | Up to [Research Unix] / researchv8dc / cmd / view2d |
1.1 root 1: #include <stdio.h>
2: #include "view2d.h"
3: #define NMAX 1030
4: char *progname;
5:
6: short timewarp;
7: double ts, te;
8: extern Rd2d rd;
9: int verbose;
10:
11: typedef struct Frame {
12: double time;
13: short *p;
14: } Frame;
15: Frame frame[3];
16: Frame *here, *tween, *there, *tframe;
17:
18: short nx, ny; /* input grid */
19: short mx, my; /* output grid */
20: float x[NMAX], y[NMAX];
21: int detrend;
22:
23: main(argc, argv)
24: int argc;
25: char **argv;
26: {
27: int i, j;
28: short *in; /* input values */
29: short *mid; /* intermediate values */
30: int nfr; /* number of frames to be drawn */
31: int ifr;
32: double twtime, timestep, t;
33: short *a, *b, *c;
34: float *d; /* trend plane */
35: int fd;
36: int tinterp = 0;
37: int ulim = 0;
38: int timebar = 0;
39: float ufmin, ufmax;
40: char *Malloc();
41:
42: timewarp = 0;
43: detrend = 0;
44: verbose = 0;
45: progname = argv[0];
46: here = frame+0;
47: tween = frame+1;
48: there = frame+2;
49: nfr = -1; mx = -1; my = -1; /* flag uninitialized values */
50:
51: for(argc--, argv++; *argv; argv++){
52: if(**argv == '-' ){
53: switch(argv[0][1]) {
54: case 'b':
55: i = sscanf(&argv[0][2], "%d", &timebar);
56: if(i!=1) timebar = -1;
57: break;
58: case 'n':
59: i = sscanf(&argv[0][2], "%hd, %hd", &mx, &my);
60: if(i == 1) { my = mx; i = 2; }
61: if((i!=2)||(mx<0)||(mx>NMAX)||(my<0)||(my>NMAX)) error("bad NX, NY");
62: break;
63: case 'f':
64: i = sscanf(&argv[0][2], "%d", &nfr);
65: tinterp = 1;
66: break;
67: case 't':
68: i = sscanf(&argv[0][2], "%E, %E", &ts, &te);
69: if(i==0) error("bad TS,TE");
70: timewarp = i;
71: break;
72: case 'r':
73: detrend++;
74: break;
75: case 'm':
76: i = sscanf(&argv[0][2], "%e, %e", &ufmin, &ufmax);
77: if(i!=2) error("bad fmin,fmax");
78: ulim++;
79: break;
80: case 'v':
81: verbose++;
82: break;
83: }
84: }else{
85: if(fd) error("can only read one file");
86: fd = Open(*argv,0);
87: }
88: }
89:
90: rd2dh(fd,&nx,&ny);
91: if((timewarp>0)&&(verbose))
92: fprintf(stderr,"timewarp=%d ts=%g te=%g\n",timewarp,ts,te);
93: if(ulim){
94: rd.fmin = ufmin;
95: rd.fmax = ufmax;
96: g_rang2();
97: }
98: if(timebar&&(ts==te)){
99: fprintf(stderr,"couldn't determine boundary times for timebar\n");
100: timebar = 0;
101: }
102: if(nfr==-1) nfr = rd.nfr;
103: if((nx > NMAX) || (ny > NMAX)) error("n too large");
104: if(mx==-1){ mx = nx; my = ny; }
105: if(verbose){
106: fprintf(stderr,"nx=%d ny=%d nframes=%d\n",nx,ny,rd.nfr);
107: fprintf(stderr,"mx=%d my=%d\n",mx,my);
108: fprintf(stderr,"global starting_time=%g ending_time=%g\n",rd.ts,rd.te);
109: fprintf(stderr,"fmin=%g fmax=%g\n",rd.fmin,rd.fmax);
110: fprintf(stderr,"pmin=%d pmax=%d\n",rd.pmin,rd.pmax);
111: fprintf(stderr,"u=%d v=%d\n",rd.u,rd.v);
112: }
113: for(i=0; i<nx; i++) x[i] = i/(nx-1.);
114: for(j=0; j<ny; j++) y[j] = j/(ny-1.);
115:
116: in = (short *)Malloc(nx*ny*sizeof(short));
117: mid = (short *)Malloc(mx*ny*sizeof(short));
118: if(timebar){
119: if(timebar==-1) timebar = my/30;
120: if(timebar<1) timebar = 1;
121: i = mx*(my+3*timebar)*sizeof(short);
122: }else{
123: i = mx*my*sizeof(short);
124: }
125: j = (nx*ny*sizeof(short)+mx*ny*sizeof(short)+3*i)/1000;
126: if(detrend) j += (8*mx*my)*sizeof(float)/1000;
127: if(j>100) fprintf(stderr,"using %dkB workspace\n",j);
128: here->p = 3*timebar*mx+(short *)Malloc(i);
129: tween->p = 3*timebar*mx+(short *)Malloc(i);
130: there->p = 3*timebar*mx+(short *)Malloc(i);
131: if(detrend) d = (float *)Malloc(mx*my*sizeof(float));
132:
133: if(rdframe(here,in,mid)==0) error("unexpected empty first frame");
134: if(detrend){
135: fit(here->p,d);
136: rmfit(here->p,d);
137: }
138:
139: if(tinterp&&(rd.nfr>1)){
140: ifr = 1;
141: twtime = ts;
142: timestep = (te-ts)/((nfr-1.)+1e-20);
143: rdframe(there,in,mid);
144: if(detrend) rmfit(there->p,d);
145: while( ifr<=nfr ){
146: while( twtime > there->time ){
147: tframe = there; there = here; here = tframe; /* swap */
148: if( rdframe(there,in,mid)==0 ) error("unexpected EOF!");
149: if(detrend) rmfit(there->p,d);
150: }
151: t = (twtime-here->time)/((there->time-here->time)+1e-20);
152: if(verbose) fprintf(stderr,"twtime=%g t=%g\n",twtime,t);
153: if( (t<0) || (t>1) ){
154: fprintf(stderr,"here=%g there=%g\n",here->time,there->time);
155: error("can't happen: extrap t=%g twtime=%g",t,twtime);
156: }
157: for(i=mx*my, a=here->p, b=there->p, c=tween->p; i>0; i--){
158: if( (*a < -BIG)||(*b < -BIG) ){
159: *c++ = -BIG-1; a++; b++;
160: }else{
161: *c++ = (1-t) * *a++ + t * *b++;
162: }
163: }
164: wrframe(twtime,tween->p,timebar);
165: ifr++;
166: twtime = ts+(ifr-1)*timestep;
167: if(twtime>te) twtime=te;
168: if(twtime>rd.te) break;
169: }
170: }else{ /* no time interpolation */
171: wrframe(here->time,here->p,timebar);
172: /*----old version-----
173: while( rdframe(there,in,mid) ){
174: if(detrend) rmfit(there->p,d);
175: wrframe(there->time,there->p,timebar);
176: tframe = there; there = here; here = tframe;
177: }
178: */
179: while( rdframe(here,in,mid) ){
180: if(detrend) rmfit(here->p,d);
181: wrframe(here->time,here->p,timebar);
182: }
183: }
184:
185: exit(0);
186: }
187:
188:
189: wrframe(t,p,timebar)
190: double t;
191: short p[];
192: int timebar;
193: {
194: if(detrend){
195: view2d(1,mx,my,t,rd.u,rd.v,0,0,0,p);
196: }else{
197: if(timebar){
198: int i,j;
199: int k = mx*(t-ts)/(te-ts);
200: short *q;
201: short space = -BIG-1; /* black */
202: short bar = rd.pmin + .8*((int)rd.pmax-rd.pmin); /* bar color */
203: short base = rd.pmin; /* base-bar color */
204: q = &p[-3*timebar*mx];
205: for(i=0;i<timebar;i++){
206: for(j=0;j<k;j++){ *q++ = bar; }
207: for(j=k;j<mx;j++){ *q++ = base; }
208: }
209: for(i=0;i<2*timebar*mx;i++){ *q++ = space; }
210: }
211: view2d(1,mx,my+3*timebar,t,rd.u,rd.v,1,rd.pmin,rd.pmax,&p[-3*timebar*mx]);
212: }
213: }
214:
215:
216: int /* returns 0 on EOF */
217: rdframe ( h, u, v )
218: Frame *h;
219: short *u;
220: short *v;
221: {
222: int i, ii, j, k, m;
223: float s, t;
224: short *w;
225:
226: if( rd2di( &h->time, u ) == 0 ) return(0);
227: if(verbose) fprintf(stderr,"time=%g\n",h->time);
228:
229: if((mx==nx)&&(my==ny)){
230: w = h->p;
231: m=nx*ny;
232: for(i = 0; i < m; i++){ *w++ = *u++; }
233: }else{
234: /* interpolate in x */
235: ii = 0;
236: for(i = 0; i < mx; i++){
237: t = i/(mx-1.);
238: while( (ii<nx-1) && (t>=x[ii+1]) ) ii++;
239: if(t == x[ii]) {
240: for(j = 0; j < ny; j++) {
241: v[i+j*mx] = u[ii+j*nx];
242: }
243: }else{
244: s = (t-x[ii])/(x[ii+1]-x[ii]);
245: for(j = 0; j < ny; j++){
246: k = ii+j*nx;
247: if( (u[k] < -BIG)||(u[k+1] < -BIG) ){
248: v[i+j*mx] = -BIG-1;
249: }else{
250: v[i+j*mx] = (u[k] + s*(u[k+1]-u[k]));
251: }
252: }
253: }
254: }
255:
256: /* interpolate in y */
257: w = h->p;
258: ii = 0;
259: for(i = 0; i < my; i++){
260: t = i/(my-1.);
261: while( (ii<ny-1) && (t>=y[ii+1]) ) ii++;
262: if(t == y[ii]){
263: for(j = 0; j < mx; j++){
264: w[j+i*mx] = v[j+ii*mx];
265: }
266: }else{
267: s = (t-y[ii])/(y[ii+1]-y[ii]);
268: for(j = 0; j < mx; j++){
269: k = j+ii*mx;
270: if( (v[k] < -BIG)||(v[k+mx] < -BIG) ){
271: w[j+i*mx] = -BIG-1;
272: }else{
273: w[j+i*mx] = (v[k] + s*(v[k+mx]-v[k]));
274: }
275: }
276: }
277: }
278: }
279: return(1);
280: }
281:
282:
283: fit(p,d)
284: short *p;
285: float *d; /* l2 fit to p */
286: {
287: float *z; /* floating copy of first frame data */
288: float *zz; /* pointer into z */
289: short *l; /* shares storage with z */
290: float *work, plane[3]; /* workspace for tl2fit */
291: float x, y;
292: float u2, v2, pow2();
293: int i, j;
294: int imx, imy, imxmy;
295: char *Malloc();
296:
297: i = mx*my*sizeof(float);
298: z = (float *)Malloc(i);
299: work = (float *)Malloc(6*i);
300: imx = mx; imy = my; imxmy = mx*my;
301: u2 = rd.u; v2 = pow2(rd.v);
302: for( i=mx*my, zz=z; i>0; i--){
303: *zz++ = ((*p++)-u2)*v2;
304: }
305: tl2fit_(&imx,&imy,&imxmy,z,work,plane);
306: free(z); free(work);
307: if(verbose) fprintf(stderr,"l2 plane: %g + %g x + %g y for 0<=x,y<=1\n",
308: plane[0], plane[1], plane[2]);
309:
310: for(j=1; j<=my; j++){
311: y = (j-1)/(my-1.);
312: for(i=1; i<=mx; i++){
313: x = (i-1)/(mx-1.);
314: *d++ = (plane[0] + x*plane[1] + y*plane[2])/v2 + u2;
315: }
316: }
317: }
318:
319: rmfit(p,d) /* p -= d */
320: short *p;
321: float *d;
322: {
323: int i;
324: float f;
325: short *p0, *p1;
326: p0 = p;
327: p1 = p+mx*my-1;
328: for(i=mx*my; i>0; i--, p++, d++){
329: f = *p - *d;
330: *p = (f>BIG)? BIG: ( (f<-BIG)? -BIG: f );
331: }
332: }
333:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.