ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/data.c
Revision: 2.8
Committed: Thu Apr 14 04:50:27 1994 UTC (30 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +3 -4 lines
Log Message:
changed extern char *libpath to extern char *getlibpath()

File Contents

# User Rev Content
1 greg 2.7 /* Copyright (c) 1993 Regents of the University of California */
2 greg 1.1
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * data.c - routines dealing with interpolated data.
9     *
10     * 6/4/86
11     */
12    
13     #include "standard.h"
14    
15     #include "color.h"
16    
17 greg 1.11 #include "resolu.h"
18    
19 greg 1.1 #include "data.h"
20    
21 greg 2.7 /* picture memory usage before warning */
22     #ifndef PSIZWARN
23     #ifdef BIGMEM
24     #define PSIZWARN 3000000
25     #else
26     #define PSIZWARN 1000000
27     #endif
28     #endif
29 greg 1.1
30 greg 2.7 #ifndef TABSIZ
31 greg 2.6 #define TABSIZ 97 /* table size (prime) */
32 greg 2.7 #endif
33 greg 2.6
34     #define hash(s) (shash(s)%TABSIZ)
35    
36    
37 greg 1.10 extern char *fgetword();
38    
39 greg 2.8 extern char *getlibpath(); /* library search path */
40 greg 1.1
41 greg 2.6 static DATARRAY *dtab[TABSIZ]; /* data array list */
42 greg 1.1
43 greg 2.6 static DATARRAY *ptab[TABSIZ]; /* picture list */
44 greg 1.1
45    
46     DATARRAY *
47     getdata(dname) /* get data array dname */
48     char *dname;
49     {
50 greg 1.10 char word[64];
51 greg 1.1 char *dfname;
52     FILE *fp;
53     int asize;
54 greg 1.4 register int i, j;
55 greg 1.1 register DATARRAY *dp;
56     /* look for array in list */
57 greg 2.6 for (dp = dtab[hash(dname)]; dp != NULL; dp = dp->next)
58 greg 1.1 if (!strcmp(dname, dp->name))
59     return(dp); /* found! */
60    
61     /*
62     * If we haven't loaded the data already, we will look
63 greg 2.5 * for it in the directories specified by the library path.
64 greg 1.1 *
65     * The file has the following format:
66     *
67 greg 1.4 * N
68 greg 1.1 * beg0 end0 n0
69     * beg1 end1 n1
70     * . . .
71 greg 1.4 * begN endN nN
72 greg 1.1 * data, later dimensions changing faster
73     * . . .
74     *
75 greg 1.4 * For irregularly spaced points, the following can be
76     * substituted for begi endi ni:
77     *
78 greg 1.6 * 0 0 ni p0i p1i .. pni
79 greg 1.1 */
80    
81 greg 2.8 if ((dfname = getpath(dname, getlibpath(), R_OK)) == NULL) {
82 greg 1.1 sprintf(errmsg, "cannot find data file \"%s\"", dname);
83     error(USER, errmsg);
84     }
85     if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
86     goto memerr;
87    
88     dp->name = savestr(dname);
89    
90     if ((fp = fopen(dfname, "r")) == NULL) {
91     sprintf(errmsg, "cannot open data file \"%s\"", dfname);
92     error(SYSTEM, errmsg);
93     }
94     /* get dimensions */
95 greg 1.10 if (fgetword(word, sizeof(word), fp) == NULL || !isint(word))
96 greg 1.1 goto scanerr;
97 greg 1.10 dp->nd = atoi(word);
98 greg 1.9 if (dp->nd <= 0 || dp->nd > MAXDDIM) {
99 greg 1.1 sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
100     error(USER, errmsg);
101     }
102     asize = 1;
103     for (i = 0; i < dp->nd; i++) {
104 greg 1.10 if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
105 greg 1.6 goto scanerr;
106 greg 1.10 dp->dim[i].org = atof(word);
107     if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
108     goto scanerr;
109     dp->dim[i].siz = atof(word);
110     if (fgetword(word, sizeof(word), fp) == NULL || !isint(word))
111     goto scanerr;
112     dp->dim[i].ne = atoi(word);
113 greg 1.6 if (dp->dim[i].ne < 2)
114     goto scanerr;
115     asize *= dp->dim[i].ne;
116     if ((dp->dim[i].siz -= dp->dim[i].org) == 0) {
117 greg 1.4 dp->dim[i].p = (double *)malloc(dp->dim[i].ne*sizeof(double));
118     if (dp->dim[i].p == NULL)
119     goto memerr;
120 greg 1.10 for (j = 0; j < dp->dim[i].ne; j++) {
121     if (fgetword(word, sizeof(word), fp) == NULL ||
122     !isflt(word))
123 greg 1.4 goto scanerr;
124 greg 1.10 dp->dim[i].p[j] = atof(word);
125     }
126 greg 1.4 for (j = 1; j < dp->dim[i].ne-1; j++)
127     if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
128     (dp->dim[i].p[j] < dp->dim[i].p[j+1]))
129     goto scanerr;
130     dp->dim[i].org = dp->dim[i].p[0];
131     dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
132     - dp->dim[i].p[0];
133     } else
134 greg 1.6 dp->dim[i].p = NULL;
135 greg 1.1 }
136     if ((dp->arr = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
137     goto memerr;
138    
139 greg 1.10 for (i = 0; i < asize; i++) {
140     if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
141 greg 1.1 goto scanerr;
142 greg 1.10 dp->arr[i] = atof(word);
143     }
144 greg 1.1 fclose(fp);
145 greg 2.6 i = hash(dname);
146     dp->next = dtab[i];
147     return(dtab[i] = dp);
148 greg 1.1
149     memerr:
150     error(SYSTEM, "out of memory in getdata");
151     scanerr:
152     sprintf(errmsg, "%s in data file \"%s\"",
153     feof(fp) ? "unexpected EOF" : "bad format", dfname);
154     error(USER, errmsg);
155     }
156    
157    
158 greg 1.5 static
159 greg 2.4 headaspect(s, iap) /* check string for aspect ratio */
160 greg 1.5 char *s;
161 greg 2.4 double *iap;
162 greg 1.5 {
163     if (isaspect(s))
164 greg 2.4 *iap *= aspectval(s);
165 greg 1.5 }
166    
167    
168 greg 1.1 DATARRAY *
169     getpict(pname) /* get picture pname */
170     char *pname;
171     {
172 greg 2.4 double inpaspect;
173 greg 1.1 char *pfname;
174     FILE *fp;
175     COLOR *scanin;
176 greg 1.11 int sl, ns;
177 greg 2.3 RESOLU inpres;
178 greg 1.11 FLOAT loc[2];
179     int y;
180     register int x, i;
181 greg 1.1 register DATARRAY *pp;
182     /* look for array in list */
183 greg 2.6 for (pp = ptab[hash(pname)]; pp != NULL; pp = pp->next)
184 greg 1.1 if (!strcmp(pname, pp->name))
185     return(pp); /* found! */
186    
187 greg 2.8 if ((pfname = getpath(pname, getlibpath(), R_OK)) == NULL) {
188 greg 1.1 sprintf(errmsg, "cannot find picture file \"%s\"", pname);
189     error(USER, errmsg);
190     }
191     if ((pp = (DATARRAY *)calloc(3, sizeof(DATARRAY))) == NULL)
192     goto memerr;
193    
194     pp[0].name =
195     pp[1].name =
196     pp[2].name = savestr(pname);
197    
198     if ((fp = fopen(pfname, "r")) == NULL) {
199     sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
200     error(SYSTEM, errmsg);
201     }
202 greg 2.3 #ifdef MSDOS
203     setmode(fileno(fp), O_BINARY);
204     #endif
205 greg 1.1 /* get dimensions */
206 greg 1.5 inpaspect = 1.0;
207 greg 2.4 getheader(fp, headaspect, &inpaspect);
208 greg 1.11 if (!fgetsresolu(&inpres, fp))
209 greg 1.1 goto readerr;
210 greg 2.7 #if PSIZWARN
211     /* check memory usage */
212     i = 3*sizeof(DATATYPE)*inpres.xr*inpres.yr;
213     if (i > PSIZWARN) {
214     sprintf(errmsg, "picture file \"%s\" using %d bytes of memory",
215     pname, i);
216     error(WARNING, errmsg);
217     }
218     #endif
219 greg 1.1 for (i = 0; i < 3; i++) {
220     pp[i].nd = 2;
221 greg 1.11 pp[i].dim[0].ne = inpres.yr;
222     pp[i].dim[1].ne = inpres.xr;
223 greg 1.1 pp[i].dim[0].org =
224     pp[i].dim[1].org = 0.0;
225 greg 1.11 if (inpres.xr <= inpres.yr*inpaspect) {
226     pp[i].dim[0].siz = inpaspect *
227     (double)inpres.yr/inpres.xr;
228 greg 1.7 pp[i].dim[1].siz = 1.0;
229     } else {
230 greg 1.1 pp[i].dim[0].siz = 1.0;
231 greg 1.11 pp[i].dim[1].siz = (double)inpres.xr/inpres.yr /
232     inpaspect;
233 greg 1.1 }
234 greg 1.7 pp[i].dim[0].p = pp[i].dim[1].p = NULL;
235 greg 1.11 pp[i].arr = (DATATYPE *)
236     malloc(inpres.xr*inpres.yr*sizeof(DATATYPE));
237 greg 1.1 if (pp[i].arr == NULL)
238     goto memerr;
239     }
240     /* load picture */
241 greg 1.11 sl = scanlen(&inpres);
242     ns = numscans(&inpres);
243     if ((scanin = (COLOR *)malloc(sl*sizeof(COLOR))) == NULL)
244 greg 1.1 goto memerr;
245 greg 1.11 for (y = 0; y < ns; y++) {
246     if (freadscan(scanin, sl, fp) < 0)
247 greg 1.1 goto readerr;
248 greg 1.11 for (x = 0; x < sl; x++) {
249     pix2loc(loc, &inpres, x, y);
250     i = (int)(loc[1]*inpres.yr)*inpres.xr +
251     (int)(loc[0]*inpres.xr);
252     pp[0].arr[i] = colval(scanin[x],RED);
253     pp[1].arr[i] = colval(scanin[x],GRN);
254     pp[2].arr[i] = colval(scanin[x],BLU);
255     }
256 greg 1.1 }
257     free((char *)scanin);
258     fclose(fp);
259 greg 2.6 i = hash(pname);
260 greg 1.1 pp[0].next =
261     pp[1].next =
262 greg 2.6 pp[2].next = ptab[i];
263     return(ptab[i] = pp);
264 greg 1.1
265     memerr:
266     error(SYSTEM, "out of memory in getpict");
267     readerr:
268     sprintf(errmsg, "bad picture file \"%s\"", pfname);
269     error(USER, errmsg);
270     }
271    
272    
273     freedata(dname) /* free memory associated with dname */
274     char *dname;
275     {
276 greg 2.5 DATARRAY head;
277 greg 2.6 int hval, nents;
278 greg 1.1 register DATARRAY *dp, *dpl;
279 greg 1.4 register int i;
280 greg 1.1
281 greg 2.6 if (dname == NULL) { /* free all if NULL */
282     hval = 0; nents = TABSIZ;
283     } else {
284     hval = hash(dname); nents = 1;
285     }
286     while (nents--) {
287     head.next = dtab[hval];
288     dpl = &head;
289     while ((dp = dpl->next) != NULL)
290     if (dname == NULL || !strcmp(dname, dp->name)) {
291     dpl->next = dp->next;
292     free((char *)dp->arr);
293     for (i = 0; i < dp->nd; i++)
294     if (dp->dim[i].p != NULL)
295     free((char *)dp->dim[i].p);
296     freestr(dp->name);
297     free((char *)dp);
298     } else
299     dpl = dp;
300     dtab[hval++] = head.next;
301     }
302 greg 1.1 }
303    
304    
305     freepict(pname) /* free memory associated with pname */
306     char *pname;
307     {
308 greg 2.5 DATARRAY head;
309 greg 2.6 int hval, nents;
310 greg 1.1 register DATARRAY *pp, *ppl;
311    
312 greg 2.6 if (pname == NULL) { /* free all if NULL */
313     hval = 0; nents = TABSIZ;
314     } else {
315     hval = hash(pname); nents = 1;
316     }
317     while (nents--) {
318     head.next = ptab[hval];
319     ppl = &head;
320     while ((pp = ppl->next) != NULL)
321     if (pname == NULL || !strcmp(pname, pp->name)) {
322     ppl->next = pp->next;
323     free((char *)pp[0].arr);
324     free((char *)pp[1].arr);
325     free((char *)pp[2].arr);
326     freestr(pp[0].name);
327     free((char *)pp);
328     } else
329     ppl = pp;
330     ptab[hval++] = head.next;
331     }
332 greg 1.1 }
333    
334    
335     double
336     datavalue(dp, pt) /* interpolate data value at a point */
337     register DATARRAY *dp;
338 greg 2.3 double *pt;
339 greg 1.1 {
340     DATARRAY sd;
341     int asize;
342 greg 1.6 int lower, upper;
343 greg 1.1 register int i;
344 greg 2.3 double x, y, y0, y1;
345 greg 1.4 /* set up dimensions for recursion */
346 greg 1.1 sd.nd = dp->nd - 1;
347     asize = 1;
348     for (i = 0; i < sd.nd; i++) {
349     sd.dim[i].org = dp->dim[i+1].org;
350     sd.dim[i].siz = dp->dim[i+1].siz;
351 greg 1.4 sd.dim[i].p = dp->dim[i+1].p;
352 greg 1.1 asize *= sd.dim[i].ne = dp->dim[i+1].ne;
353     }
354 greg 1.4 /* get independent variable */
355     if (dp->dim[0].p == NULL) { /* evenly spaced points */
356     x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
357     x = x * (dp->dim[0].ne - 1);
358     i = x;
359     if (i < 0)
360     i = 0;
361     else if (i > dp->dim[0].ne - 2)
362     i = dp->dim[0].ne - 2;
363     } else { /* unevenly spaced points */
364 greg 1.6 if (dp->dim[0].siz > 0.0) {
365     lower = 0;
366     upper = dp->dim[0].ne;
367     } else {
368     lower = dp->dim[0].ne;
369     upper = 0;
370     }
371     do {
372     i = (lower + upper) >> 1;
373     if (pt[0] >= dp->dim[0].p[i])
374     lower = i;
375 greg 1.8 else
376 greg 1.6 upper = i;
377     } while (i != (lower + upper) >> 1);
378 greg 1.8 if (i > dp->dim[0].ne - 2)
379 greg 1.4 i = dp->dim[0].ne - 2;
380     x = i + (pt[0] - dp->dim[0].p[i]) /
381     (dp->dim[0].p[i+1] - dp->dim[0].p[i]);
382     }
383     /* get dependent variable */
384 greg 1.1 if (dp->nd == 1) {
385     y0 = dp->arr[i];
386     y1 = dp->arr[i+1];
387     } else {
388     sd.arr = &dp->arr[i*asize];
389     y0 = datavalue(&sd, pt+1);
390     sd.arr = &dp->arr[(i+1)*asize];
391     y1 = datavalue(&sd, pt+1);
392     }
393     /*
394     * Extrapolate as far as one division, then
395     * taper off harmonically to zero.
396     */
397     if (x > i+2)
398     y = (2*y1-y0)/(x-i-1);
399     else if (x < i-1)
400     y = (2*y0-y1)/(i-x);
401     else
402     y = y0*((i+1)-x) + y1*(x-i);
403    
404     return(y);
405     }