ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/image.c
Revision: 2.16
Committed: Sat Feb 22 02:07:22 2003 UTC (21 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.15: +75 -14 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

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