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