ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/data.c
Revision: 2.31
Committed: Tue Jul 8 18:25:00 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad4R2, rad4R2P1
Changes since 2.30: +13 -13 lines
Log Message:
Eliminated unnecessary "extern" and "register" modifiers

File Contents

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