ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/image.c
Revision: 2.34
Committed: Tue Mar 11 02:21:46 2008 UTC (16 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R9
Changes since 2.33: +47 -11 lines
Log Message:
Added planisphere view type (-vts option) as requested by Axel Jacobs

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: image.c,v 2.33 2006/11/19 01:14:31 greg Exp $";
3 #endif
4 /*
5 * image.c - routines for image generation.
6 *
7 * External symbols declared in view.h
8 */
9
10 #include "copyright.h"
11
12 #include <ctype.h>
13 #include "rtio.h"
14 #include "rtmath.h"
15 #include "paths.h"
16 #include "view.h"
17
18
19 #define FEQ(x,y) (fabs((x)-(y)) <= FTINY)
20 #define VEQ(v,w) (FEQ((v)[0],(w)[0]) && FEQ((v)[1],(w)[1]) \
21 && FEQ((v)[2],(w)[2]))
22
23 VIEW stdview = STDVIEW; /* default view parameters */
24
25 static gethfunc gethview;
26
27
28 char *
29 setview(v) /* set hvec and vvec, return message on error */
30 register VIEW *v;
31 {
32 static char ill_horiz[] = "illegal horizontal view size";
33 static char ill_vert[] = "illegal vertical view size";
34
35 if (v->vaft < -FTINY || (v->vaft > FTINY && v->vaft <= v->vfore))
36 return("illegal fore/aft clipping plane");
37
38 if (v->vdist <= FTINY)
39 return("illegal view distance");
40 v->vdist *= normalize(v->vdir); /* normalize direction */
41 if (v->vdist == 0.0)
42 return("zero view direction");
43
44 if (normalize(v->vup) == 0.0) /* normalize view up */
45 return("zero view up vector");
46
47 fcross(v->hvec, v->vdir, v->vup); /* compute horiz dir */
48
49 if (normalize(v->hvec) == 0.0)
50 return("view up parallel to view direction");
51
52 fcross(v->vvec, v->hvec, v->vdir); /* compute vert dir */
53
54 if (v->horiz <= FTINY)
55 return(ill_horiz);
56 if (v->vert <= FTINY)
57 return(ill_vert);
58
59 switch (v->type) {
60 case VT_PAR: /* parallel view */
61 v->hn2 = v->horiz;
62 v->vn2 = v->vert;
63 break;
64 case VT_PER: /* perspective view */
65 if (v->horiz >= 180.0-FTINY)
66 return(ill_horiz);
67 if (v->vert >= 180.0-FTINY)
68 return(ill_vert);
69 v->hn2 = 2.0 * tan(v->horiz*(PI/180.0/2.0));
70 v->vn2 = 2.0 * tan(v->vert*(PI/180.0/2.0));
71 break;
72 case VT_CYL: /* cylindrical panorama */
73 if (v->horiz > 360.0+FTINY)
74 return(ill_horiz);
75 if (v->vert >= 180.0-FTINY)
76 return(ill_vert);
77 v->hn2 = v->horiz * (PI/180.0);
78 v->vn2 = 2.0 * tan(v->vert*(PI/180.0/2.0));
79 break;
80 case VT_ANG: /* angular fisheye */
81 if (v->horiz > 360.0+FTINY)
82 return(ill_horiz);
83 if (v->vert > 360.0+FTINY)
84 return(ill_vert);
85 v->hn2 = v->horiz * (PI/180.0);
86 v->vn2 = v->vert * (PI/180.0);
87 break;
88 case VT_HEM: /* hemispherical fisheye */
89 if (v->horiz > 180.0+FTINY)
90 return(ill_horiz);
91 if (v->vert > 180.0+FTINY)
92 return(ill_vert);
93 v->hn2 = 2.0 * sin(v->horiz*(PI/180.0/2.0));
94 v->vn2 = 2.0 * sin(v->vert*(PI/180.0/2.0));
95 break;
96 case VT_PLS: /* planispheric fisheye */
97 if (v->horiz >= 360.0-FTINY)
98 return(ill_horiz);
99 if (v->vert >= 360.0-FTINY)
100 return(ill_vert);
101 v->hn2 = 2.*sin(v->horiz*(PI/180.0/2.0)) /
102 (1.0 + cos(v->horiz*(PI/180.0/2.0)));
103 v->vn2 = 2.*sin(v->vert*(PI/180.0/2.0)) /
104 (1.0 + cos(v->vert*(PI/180.0/2.0)));
105 break;
106 default:
107 return("unknown view type");
108 }
109 if (v->type != VT_ANG && v->type != VT_PLS) {
110 if (v->type != VT_CYL) {
111 v->hvec[0] *= v->hn2;
112 v->hvec[1] *= v->hn2;
113 v->hvec[2] *= v->hn2;
114 }
115 v->vvec[0] *= v->vn2;
116 v->vvec[1] *= v->vn2;
117 v->vvec[2] *= v->vn2;
118 }
119 v->hn2 *= v->hn2;
120 v->vn2 *= v->vn2;
121
122 return(NULL);
123 }
124
125
126 void
127 normaspect(va, ap, xp, yp) /* fix pixel aspect or resolution */
128 double va; /* view aspect ratio */
129 double *ap; /* pixel aspect in (or out if 0) */
130 int *xp, *yp; /* x and y resolution in (or out if *ap!=0) */
131 {
132 if (*ap <= FTINY)
133 *ap = va * *xp / *yp; /* compute pixel aspect */
134 else if (va * *xp > *ap * *yp)
135 *xp = *yp / va * *ap + .5; /* reduce x resolution */
136 else
137 *yp = *xp * va / *ap + .5; /* reduce y resolution */
138 }
139
140
141 double
142 viewray(orig, direc, v, x, y) /* compute ray origin and direction */
143 FVECT orig, direc;
144 register VIEW *v;
145 double x, y;
146 {
147 double d, z;
148
149 x += v->hoff - 0.5;
150 y += v->voff - 0.5;
151
152 switch(v->type) {
153 case VT_PAR: /* parallel view */
154 orig[0] = v->vp[0] + v->vfore*v->vdir[0]
155 + x*v->hvec[0] + y*v->vvec[0];
156 orig[1] = v->vp[1] + v->vfore*v->vdir[1]
157 + x*v->hvec[1] + y*v->vvec[1];
158 orig[2] = v->vp[2] + v->vfore*v->vdir[2]
159 + x*v->hvec[2] + y*v->vvec[2];
160 VCOPY(direc, v->vdir);
161 return(v->vaft > FTINY ? v->vaft - v->vfore : 0.0);
162 case VT_PER: /* perspective view */
163 direc[0] = v->vdir[0] + x*v->hvec[0] + y*v->vvec[0];
164 direc[1] = v->vdir[1] + x*v->hvec[1] + y*v->vvec[1];
165 direc[2] = v->vdir[2] + x*v->hvec[2] + y*v->vvec[2];
166 orig[0] = v->vp[0] + v->vfore*direc[0];
167 orig[1] = v->vp[1] + v->vfore*direc[1];
168 orig[2] = v->vp[2] + v->vfore*direc[2];
169 d = normalize(direc);
170 return(v->vaft > FTINY ? (v->vaft - v->vfore)*d : 0.0);
171 case VT_HEM: /* hemispherical fisheye */
172 z = 1.0 - x*x*v->hn2 - y*y*v->vn2;
173 if (z < 0.0)
174 return(-1.0);
175 z = sqrt(z);
176 direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0];
177 direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1];
178 direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2];
179 orig[0] = v->vp[0] + v->vfore*direc[0];
180 orig[1] = v->vp[1] + v->vfore*direc[1];
181 orig[2] = v->vp[2] + v->vfore*direc[2];
182 return(v->vaft > FTINY ? v->vaft - v->vfore : 0.0);
183 case VT_CYL: /* cylindrical panorama */
184 d = x * v->horiz * (PI/180.0);
185 z = cos(d);
186 x = sin(d);
187 direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0];
188 direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1];
189 direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2];
190 orig[0] = v->vp[0] + v->vfore*direc[0];
191 orig[1] = v->vp[1] + v->vfore*direc[1];
192 orig[2] = v->vp[2] + v->vfore*direc[2];
193 d = normalize(direc);
194 return(v->vaft > FTINY ? (v->vaft - v->vfore)*d : 0.0);
195 case VT_ANG: /* angular fisheye */
196 x *= (1.0/180.0)*v->horiz;
197 y *= (1.0/180.0)*v->vert;
198 d = x*x + y*y;
199 if (d > 1.0)
200 return(-1.0);
201 d = sqrt(d);
202 z = cos(PI*d);
203 d = d <= FTINY ? PI : sqrt(1.0 - z*z)/d;
204 x *= d;
205 y *= d;
206 direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0];
207 direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1];
208 direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2];
209 orig[0] = v->vp[0] + v->vfore*direc[0];
210 orig[1] = v->vp[1] + v->vfore*direc[1];
211 orig[2] = v->vp[2] + v->vfore*direc[2];
212 return(v->vaft > FTINY ? v->vaft - v->vfore : 0.0);
213 case VT_PLS: /* planispheric fisheye */
214 x *= sqrt(v->hn2);
215 y *= sqrt(v->vn2);
216 d = x*x + y*y;
217 z = (1. - d)/(1. + d);
218 d = d <= FTINY*FTINY ? PI : sqrt((1.0 - z*z)/d);
219 x *= d;
220 y *= d;
221 direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0];
222 direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1];
223 direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2];
224 orig[0] = v->vp[0] + v->vfore*direc[0];
225 orig[1] = v->vp[1] + v->vfore*direc[1];
226 orig[2] = v->vp[2] + v->vfore*direc[2];
227 return(v->vaft > FTINY ? v->vaft - v->vfore : 0.0);
228 }
229 return(-1.0);
230 }
231
232
233 void
234 viewloc(ip, v, p) /* find image location for point */
235 FVECT ip;
236 register VIEW *v;
237 FVECT p;
238 {
239 double d, d2;
240 FVECT disp;
241
242 VSUB(disp, p, v->vp);
243
244 switch (v->type) {
245 case VT_PAR: /* parallel view */
246 ip[2] = DOT(disp,v->vdir) - v->vfore;
247 break;
248 case VT_PER: /* perspective view */
249 d = DOT(disp,v->vdir);
250 ip[2] = VLEN(disp);
251 if (d < 0.0) { /* fold pyramid */
252 ip[2] = -ip[2];
253 d = -d;
254 }
255 if (d > FTINY) {
256 d = 1.0/d;
257 disp[0] *= d;
258 disp[1] *= d;
259 disp[2] *= d;
260 }
261 ip[2] *= (1.0 - v->vfore*d);
262 break;
263 case VT_HEM: /* hemispherical fisheye */
264 d = normalize(disp);
265 if (DOT(disp,v->vdir) < 0.0)
266 ip[2] = -d;
267 else
268 ip[2] = d;
269 ip[2] -= v->vfore;
270 break;
271 case VT_CYL: /* cylindrical panorama */
272 d = DOT(disp,v->hvec);
273 d2 = DOT(disp,v->vdir);
274 ip[0] = 180.0/PI * atan2(d,d2) / v->horiz + 0.5 - v->hoff;
275 d = 1.0/sqrt(d*d + d2*d2);
276 ip[1] = DOT(disp,v->vvec)*d/v->vn2 + 0.5 - v->voff;
277 ip[2] = VLEN(disp);
278 ip[2] *= (1.0 - v->vfore*d);
279 return;
280 case VT_ANG: /* angular fisheye */
281 ip[0] = 0.5 - v->hoff;
282 ip[1] = 0.5 - v->voff;
283 ip[2] = normalize(disp) - v->vfore;
284 d = DOT(disp,v->vdir);
285 if (d >= 1.0-FTINY)
286 return;
287 if (d <= -(1.0-FTINY)) {
288 ip[0] += 180.0/v->horiz;
289 return;
290 }
291 d = (180.0/PI)*acos(d) / sqrt(1.0 - d*d);
292 ip[0] += DOT(disp,v->hvec)*d/v->horiz;
293 ip[1] += DOT(disp,v->vvec)*d/v->vert;
294 return;
295 case VT_PLS: /* planispheric fisheye */
296 ip[0] = 0.5 - v->hoff;
297 ip[1] = 0.5 - v->voff;
298 ip[2] = normalize(disp) - v->vfore;
299 d = DOT(disp,v->vdir);
300 if (d >= 1.0-FTINY)
301 return;
302 if (d <= -(1.0-FTINY))
303 return; /* really an error */
304 d = sqrt(1.0 - d*d) / (1.0 + d);
305 ip[0] += DOT(disp,v->hvec)*d/sqrt(v->hn2);
306 ip[1] += DOT(disp,v->vvec)*d/sqrt(v->vn2);
307 return;
308 }
309 ip[0] = DOT(disp,v->hvec)/v->hn2 + 0.5 - v->hoff;
310 ip[1] = DOT(disp,v->vvec)/v->vn2 + 0.5 - v->voff;
311 }
312
313
314 void
315 pix2loc(loc, rp, px, py) /* compute image location from pixel pos. */
316 RREAL loc[2];
317 register RESOLU *rp;
318 int px, py;
319 {
320 register int x, y;
321
322 if (rp->rt & YMAJOR) {
323 x = px;
324 y = py;
325 } else {
326 x = py;
327 y = px;
328 }
329 if (rp->rt & XDECR)
330 x = rp->xr-1 - x;
331 if (rp->rt & YDECR)
332 y = rp->yr-1 - y;
333 loc[0] = (x+.5)/rp->xr;
334 loc[1] = (y+.5)/rp->yr;
335 }
336
337
338 void
339 loc2pix(pp, rp, lx, ly) /* compute pixel pos. from image location */
340 int pp[2];
341 register RESOLU *rp;
342 double lx, ly;
343 {
344 register int x, y;
345
346 x = lx * rp->xr;
347 y = ly * rp->yr;
348 if (rp->rt & XDECR)
349 x = rp->xr-1 - x;
350 if (rp->rt & YDECR)
351 y = rp->yr-1 - y;
352 if (rp->rt & YMAJOR) {
353 pp[0] = x;
354 pp[1] = y;
355 } else {
356 pp[0] = y;
357 pp[1] = x;
358 }
359 }
360
361
362 int
363 getviewopt(v, ac, av) /* process view argument */
364 register VIEW *v;
365 int ac;
366 register char *av[];
367 {
368 #define check(c,l) if ((av[0][c]&&av[0][c]!=' ') || \
369 badarg(ac-1,av+1,l)) return(-1)
370
371 if (ac <= 0 || av[0][0] != '-' || av[0][1] != 'v')
372 return(-1);
373 switch (av[0][2]) {
374 case 't': /* type */
375 if (!av[0][3] || av[0][3]==' ')
376 return(-1);
377 check(4,"");
378 v->type = av[0][3];
379 return(0);
380 case 'p': /* point */
381 check(3,"fff");
382 v->vp[0] = atof(av[1]);
383 v->vp[1] = atof(av[2]);
384 v->vp[2] = atof(av[3]);
385 return(3);
386 case 'd': /* direction */
387 check(3,"fff");
388 v->vdir[0] = atof(av[1]);
389 v->vdir[1] = atof(av[2]);
390 v->vdir[2] = atof(av[3]);
391 v->vdist = 1.;
392 return(3);
393 case 'u': /* up */
394 check(3,"fff");
395 v->vup[0] = atof(av[1]);
396 v->vup[1] = atof(av[2]);
397 v->vup[2] = atof(av[3]);
398 return(3);
399 case 'h': /* horizontal size */
400 check(3,"f");
401 v->horiz = atof(av[1]);
402 return(1);
403 case 'v': /* vertical size */
404 check(3,"f");
405 v->vert = atof(av[1]);
406 return(1);
407 case 'o': /* fore clipping plane */
408 check(3,"f");
409 v->vfore = atof(av[1]);
410 return(1);
411 case 'a': /* aft clipping plane */
412 check(3,"f");
413 v->vaft = atof(av[1]);
414 return(1);
415 case 's': /* shift */
416 check(3,"f");
417 v->hoff = atof(av[1]);
418 return(1);
419 case 'l': /* lift */
420 check(3,"f");
421 v->voff = atof(av[1]);
422 return(1);
423 default:
424 return(-1);
425 }
426 #undef check
427 }
428
429
430 int
431 sscanview(vp, s) /* get view parameters from string */
432 VIEW *vp;
433 register char *s;
434 {
435 int ac;
436 char *av[4];
437 int na;
438 int nvopts = 0;
439
440 while (isspace(*s))
441 if (!*s++)
442 return(0);
443 while (*s) {
444 ac = 0;
445 do {
446 if (ac || *s == '-')
447 av[ac++] = s;
448 while (*s && !isspace(*s))
449 s++;
450 while (isspace(*s))
451 s++;
452 } while (*s && ac < 4);
453 if ((na = getviewopt(vp, ac, av)) >= 0) {
454 if (na+1 < ac)
455 s = av[na+1];
456 nvopts++;
457 } else if (ac > 1)
458 s = av[1];
459 }
460 return(nvopts);
461 }
462
463
464 void
465 fprintview(vp, fp) /* write out view parameters */
466 register VIEW *vp;
467 FILE *fp;
468 {
469 fprintf(fp, " -vt%c", vp->type);
470 fprintf(fp, " -vp %.6g %.6g %.6g", vp->vp[0], vp->vp[1], vp->vp[2]);
471 fprintf(fp, " -vd %.6g %.6g %.6g", vp->vdir[0]*vp->vdist,
472 vp->vdir[1]*vp->vdist,
473 vp->vdir[2]*vp->vdist);
474 fprintf(fp, " -vu %.6g %.6g %.6g", vp->vup[0], vp->vup[1], vp->vup[2]);
475 fprintf(fp, " -vh %.6g -vv %.6g", vp->horiz, vp->vert);
476 fprintf(fp, " -vo %.6g -va %.6g", vp->vfore, vp->vaft);
477 fprintf(fp, " -vs %.6g -vl %.6g", vp->hoff, vp->voff);
478 }
479
480
481 char *
482 viewopt(vp) /* translate to minimal view string */
483 register VIEW *vp;
484 {
485 static char vwstr[128];
486 register char *cp = vwstr;
487
488 *cp = '\0';
489 if (vp->type != stdview.type) {
490 sprintf(cp, " -vt%c", vp->type);
491 cp += strlen(cp);
492 }
493 if (!VEQ(vp->vp,stdview.vp)) {
494 sprintf(cp, " -vp %.6g %.6g %.6g",
495 vp->vp[0], vp->vp[1], vp->vp[2]);
496 cp += strlen(cp);
497 }
498 if (!FEQ(vp->vdist,stdview.vdist) || !VEQ(vp->vdir,stdview.vdir)) {
499 sprintf(cp, " -vd %.6g %.6g %.6g",
500 vp->vdir[0]*vp->vdist,
501 vp->vdir[1]*vp->vdist,
502 vp->vdir[2]*vp->vdist);
503 cp += strlen(cp);
504 }
505 if (!VEQ(vp->vup,stdview.vup)) {
506 sprintf(cp, " -vu %.6g %.6g %.6g",
507 vp->vup[0], vp->vup[1], vp->vup[2]);
508 cp += strlen(cp);
509 }
510 if (!FEQ(vp->horiz,stdview.horiz)) {
511 sprintf(cp, " -vh %.6g", vp->horiz);
512 cp += strlen(cp);
513 }
514 if (!FEQ(vp->vert,stdview.vert)) {
515 sprintf(cp, " -vv %.6g", vp->vert);
516 cp += strlen(cp);
517 }
518 if (!FEQ(vp->vfore,stdview.vfore)) {
519 sprintf(cp, " -vo %.6g", vp->vfore);
520 cp += strlen(cp);
521 }
522 if (!FEQ(vp->vaft,stdview.vaft)) {
523 sprintf(cp, " -va %.6g", vp->vaft);
524 cp += strlen(cp);
525 }
526 if (!FEQ(vp->hoff,stdview.hoff)) {
527 sprintf(cp, " -vs %.6g", vp->hoff);
528 cp += strlen(cp);
529 }
530 if (!FEQ(vp->voff,stdview.voff)) {
531 sprintf(cp, " -vl %.6g", vp->voff);
532 cp += strlen(cp);
533 }
534 return(vwstr);
535 }
536
537
538 int
539 isview(s) /* is this a view string? */
540 char *s;
541 {
542 static char *altname[]={NULL,VIEWSTR,"rpict","rview","rvu","rpiece","pinterp",NULL};
543 extern char *progname;
544 register char *cp;
545 register char **an;
546 /* add program name to list */
547 if (altname[0] == NULL) {
548 for (cp = progname; *cp; cp++)
549 ;
550 while (cp > progname && !ISDIRSEP(cp[-1]))
551 cp--;
552 altname[0] = cp;
553 }
554 /* skip leading path */
555 cp = s;
556 while (*cp && *cp != ' ')
557 cp++;
558 while (cp > s && !ISDIRSEP(cp[-1]))
559 cp--;
560 for (an = altname; *an != NULL; an++)
561 if (!strncmp(*an, cp, strlen(*an)))
562 return(1);
563 return(0);
564 }
565
566
567 struct myview {
568 VIEW *hv;
569 int ok;
570 };
571
572
573 static int
574 gethview( /* get view from header */
575 char *s,
576 void *v
577 )
578 {
579 if (isview(s) && sscanview(((struct myview*)v)->hv, s) > 0)
580 ((struct myview*)v)->ok++;
581 return(0);
582 }
583
584
585 int
586 viewfile(fname, vp, rp) /* get view from file */
587 char *fname;
588 VIEW *vp;
589 RESOLU *rp;
590 {
591 struct myview mvs;
592 FILE *fp;
593
594 if (fname == NULL || !strcmp(fname, "-"))
595 fp = stdin;
596 else if ((fp = fopen(fname, "r")) == NULL)
597 return(-1);
598
599 mvs.hv = vp;
600 mvs.ok = 0;
601
602 getheader(fp, gethview, &mvs);
603
604 if (rp != NULL && !fgetsresolu(rp, fp))
605 mvs.ok = 0;
606
607 fclose(fp);
608
609 return(mvs.ok);
610 }