|
|
1.1 root 1: /****** routines for reading view2d files ******
2:
3: globals that must be declared by calling program:
4: double ts; starting time
5: double te; ending time
6: short timewarp; == -1 ts, te unused; avoids seeks if possible
7: == 0 ts, te set from file
8: == 1 ts set by calling program, te set from file
9: == 2 ts, te set by calling program
10: another global, declared here in rd.c (with structure definition in view2d.h)
11: Rd2d rd;
12: short rd.pmin == global min pixel value
13: float rd.pmin == global max
14: short rd.u == global u value (origin)
15: short rd.v == global v value (exponent)
16: int rd.nfr == number of frames (not set if timewarp==-1)
17:
18: if p is the pixel value and f is the corresponding user's function value,
19: f = (p-u) * 2**v p = u + f / 2**v
20: if p<(-BIG) then the function value is undefined, and the pixel is black
21:
22: rd2dh(fd,pny,pny)
23: int fd; (on input) a file descriptor open for read and seek
24: short *pnx, *pny; (on output) size of the images
25: initializes the global variables ts, te, rd; chooses global scale factor;
26: skips forward through file to time>=ts;
27:
28: rd2di(ptime,p)
29: double *ptime; frame time
30: short p[]; nx by ny array of pixel values
31: reads a frame of data;
32:
33: rd2dj(j)
34: int j;
35: jumps to j frames forward (or back, if j<0)
36:
37: For a sample program illustrating the use of rd.c, see blue.c.
38: */
39:
40:
41: #include <stdio.h>
42: #include "view2d.h"
43: extern short timewarp;
44: extern double ts, te;
45: Rd2d rd;
46:
47: Header hd;
48: Header2 hd2;
49: short filver; /* normally 3; old files may still be 2 */
50: double hdtime; /* set by rdhd() */
51: double prevtim; /* old hdtime, to check monotonicity */
52: long frsiz; /* size of header and image */
53: int hdr_rd; /* are we part way through reading a header? */
54: #define MAXN 1025
55:
56: rd2dh(fd,pnx,pny) /* initialize */
57: int fd;
58: short *pnx, *pny;
59: {
60: long firstfr, imsiz;
61: int i, j;
62: float v2;
63: int hdsiz;
64: long Lseek();
65: double pow2();
66:
67: /* read first header */
68: rd.fd = fd;
69: filver = 0;
70: hdr_rd = 0;
71: if(rdhd()==0) error("empty input\n");
72: if(filver==3){ hdsiz = sizeof(hd); }
73: else{ hdsiz = sizeof(hd2); }
74: if( (hd.nx>MAXN) || (hd.ny>MAXN) ) error("nx,ny too big");
75: *pnx = hd.nx;
76: *pny = hd.ny;
77: rd.siz = hd.nx * hd.ny;
78: imsiz = rd.siz*sizeof(short);
79: frsiz = hdsiz+imsiz;
80: rd.ts = hdtime;
81: rd.u = hd.u;
82: rd.v = hd.v;
83: rd.fixuv = hd.fixuv;
84: rd.pmin = hd.pmin;
85: rd.pmax = hd.pmax;
86: if(timewarp>=0){
87: Lseek(fd,-frsiz,2);
88: rdhd();
89: rd.te = hdtime;
90: }
91:
92: /* locate first active header */
93: prevtim = -1.e28;
94: firstfr = 0;
95: if(timewarp>0){
96: Lseek(fd,firstfr,0);
97: rdhd();
98: if(ts<hdtime) ts = hdtime; /* don't extrapolate back in time */
99: while(hdtime<ts){
100: firstfr = Lseek(fd,imsiz,1);
101: if( rdhd() == 0 ) error("starting time is after EOF");
102: }
103: if(ts!=hdtime){
104: firstfr = Lseek(fd,-frsiz-hdsiz,1);
105: rdhd();
106: }
107: prevtim = hdtime;
108: }
109: if (timewarp==0){ te = rd.te; ts = rd.ts; }
110: else if(timewarp==1){ te = rd.te; }
111:
112: if( (timewarp>=0) || (rd.fixuv!=1) ){
113: rd.nfr = Lseek(fd,0L,2)/frsiz;
114: Lseek(fd,firstfr,0);
115: }else if( rd.fixuv==1 ){
116: hdr_rd = 1; /* so that initial header needn't be reread */
117: }
118:
119: /* global scaling */
120: if(rd.fixuv!=1){
121: g_range();
122: Lseek(fd,firstfr,0);
123: }else{
124: v2 = pow2(rd.v);
125: rd.fmin = (rd.pmin-rd.u)*v2;
126: rd.fmax = (rd.pmax-rd.u)*v2;
127: }
128: }
129:
130:
131:
132: int /* returns 0 upon EOF */
133: rd2di(ptime,p) /* read an image */
134: double *ptime;
135: short p[];
136: {
137: short *q;
138: int shift;
139: int i, j, k;
140: int pmin=rd.pmin;
141: int pmax=rd.pmax;
142:
143: if(hdr_rd==1){ hdr_rd=0; }
144: else if(rdhd()==0) return(0);
145: if(hdtime<prevtim) fprintf(stderr," time decreased! new=%g old=%g\n",hdtime,prevtim);
146: prevtim = hdtime;
147: *ptime = hdtime;
148:
149: Read(rd.fd,p,rd.siz*sizeof(short));
150:
151: /* convert to global units */
152: if( (hd.u!=rd.u) || (hd.v!=rd.v) ){
153: shift = rd.v - hd.v;
154: if(shift>15){
155: for( i=rd.siz, q=p; i>0; i--, q++ ){
156: *q = rd.u;
157: }
158: }else if(shift>=0){
159: for(k=1; shift>0;){ shift--; k*=2; }
160: for( i=rd.siz, q=p; i>0; i--, q++ ){
161: if( *q < -BIG ) continue;
162: j = ((*q-hd.u)/k) + rd.u;
163: *q = (j<pmin)? pmin:
164: ((j>pmax)? pmax: j);
165: }
166: }else{ /* shift<0 */
167: for(k=1; shift<0;){ shift++; k*=2; }
168: for( i=rd.siz, q=p; i>0; i--, q++ ){
169: if( *q < -BIG ) continue;
170: j = ((*q-hd.u)*k) + rd.u;
171: *q = (j<pmin)? pmin:
172: ((j>pmax)? pmax: j);
173: }
174: }
175: }else{ /* hd.u=rd.u and hd.v==rd.v */
176: for( i=rd.siz, q=p; i>0; i--, q++ ){
177: if( *q < -BIG ) continue;
178: if( *q<rd.pmin ){ *q = rd.pmin; }
179: else if( *q>rd.pmax ){ *q = rd.pmax; }
180: }
181: }
182: return(1);
183: }
184:
185:
186:
187: int
188: rdhd() /* returns 0 upon EOF */
189: {
190: long Lseek();
191:
192: if( filver==0 ){ /* get version number first time through */
193: if( EOFRead(rd.fd,&hd,sizeof(hd)) == 0 ) return(0);
194: if(hd.magic!=MAGIC) error("bad magic number %o at start",hd.magic);
195: filver = hd.ver;
196: if(filver==2){ Lseek(rd.fd,0L,0); }
197: else if(filver!=3){ error("bad version number %d",filver); }
198: }else if(filver==3){
199: if( EOFRead(rd.fd,&hd,sizeof(hd)) == 0 ) return(0);
200: }
201:
202: if(filver==2){
203: if( EOFRead(rd.fd,&hd2,sizeof(hd2)) == 0 ) return(0);
204: hd.magic = hd2.magic;
205: hd.ver = hd2.ver;
206: hd.nx = hd2.nx;
207: hd.ny = hd2.ny;
208: hd.u = hd2.u;
209: hd.v = hd2.v;
210: hd.fixuv = hd2.fixuv;
211: hd.pmin = hd2.pmin;
212: hd.pmax = hd2.pmax;
213: sprintf(hd.time, "%15.8g", (double)hd2.time);
214: }
215:
216: if(hd.magic!=MAGIC){
217: fprintf(stderr,"bad magic number %o\n",hd.magic);
218: fprintf(stderr,"file pointer now at %ld\n",Lseek(rd.fd,0L,1));
219: exit(2);
220: }
221: if(hd.ver!=filver){ error("inconsistent version number"); }
222: if(sscanf(hd.time,"%16E",&hdtime)!=1) error("can't parse time: %s",hd.time);
223: return(1);
224: }
225:
226:
227:
228: g_range()
229: { /* scan until last active header, updating range */
230: float gmin, gmax;
231: f_range(&rd.fmin,&rd.fmax);
232: while( f_range(&gmin,&gmax) ){
233: if( (timewarp==2) && (hdtime<te) ) break;
234: if( rd.v > hd.v ){
235: rd.u = hd.u;
236: rd.v = hd.v;
237: }
238: if( gmin < rd.fmin ) rd.fmin = gmin;
239: if( gmax > rd.fmax ) rd.fmax = gmax;
240: rd.te = hdtime;
241: }
242: g_rang2();
243: }
244:
245: g_rang2()
246: { /* choose global scaling */
247: float gmin, gmax, v2;
248: double pow2();
249: gmax = (rd.fmax<0)?0:rd.fmax;
250: gmin = (rd.fmin>0)?0:rd.fmin;
251: v2 = pow2(rd.v);
252: if( (gmin < (-BIG-rd.u)*v2) || (gmax > (BIG-rd.u)*v2) ){
253: while( (gmax-gmin)/v2 >= 2*BIG-2 ){ rd.v++; v2 *= 2; }
254: rd.u = (short)(-(gmin+gmax)/(2*v2));
255: }
256: rd.pmin = rd.u + (short)(rd.fmin/v2);
257: rd.pmax = rd.u + (short)(rd.fmax/v2);
258: }
259:
260:
261: int
262: f_range ( pfmin, pfmax )
263: float *pfmin, *pfmax;
264: {
265: int i, j;
266: short p[MAXN], *q;
267: short pmin = BIG;
268: short pmax = -BIG;
269: double pow2();
270: float u2, v2;
271: if( rdhd()==0 ) return(0);
272: u2 = hd.u;
273: v2 = pow2(hd.v);
274: for( i=hd.ny; i>0; i-- ){
275: Read(rd.fd,p,hd.nx*sizeof(short));
276: for( j=hd.nx, q=p; j>0; j--, q++ ){
277: if( *q < -BIG ) continue;
278: if( *q < pmin ){ pmin = *q; }
279: if( *q > pmax ){ pmax = *q; }
280: }
281: }
282: *pfmin = (pmin-u2)*v2;
283: *pfmax = (pmax-u2)*v2;
284: return(1);
285: }
286:
287: rd2dj(j)
288: int j;
289: {
290: long Lseek();
291: Lseek(rd.fd,j*frsiz,1);
292: prevtim = ts;
293: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.