ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.29
Committed: Wed Jul 9 11:23:07 1997 UTC (27 years, 3 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 2.28: +2 -2 lines
Log Message:
changed rtrace behavior to send bogus record when ray dir==0

File Contents

# Content
1 /* Copyright (c) 1995 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * Interpolate and extrapolate pictures with different view parameters.
9 *
10 * Greg Ward 09Dec89
11 */
12
13 #include "standard.h"
14
15 #include <ctype.h>
16
17 #include "view.h"
18
19 #include "color.h"
20
21 #include "resolu.h"
22
23 #define LOG2 0.69314718055994530942
24
25 #define pscan(y) (ourpict+(y)*hresolu)
26 #define sscan(y) (ourspict+(y)*hresolu)
27 #define wscan(y) (ourweigh+(y)*hresolu)
28 #define zscan(y) (ourzbuf+(y)*hresolu)
29 #define bscan(y) (ourbpict+(y)*hresolu)
30 #define averaging (ourweigh != NULL)
31 #define blurring (ourbpict != NULL)
32 #define usematrix (hasmatrix & !averaging)
33 #define zisnorm (!usematrix | ourview.type != VT_PER)
34
35 #define MAXWT 1000. /* maximum pixel weight (averaging) */
36
37 #define F_FORE 1 /* fill foreground */
38 #define F_BACK 2 /* fill background */
39
40 #define PACKSIZ 256 /* max. calculation packet size */
41
42 #define RTCOM "rtrace -h- -ovl -fff "
43
44 #define ABS(x) ((x)>0?(x):-(x))
45
46 struct position {int x,y; float z;};
47
48 #define NSTEPS 64 /* number steps in overlap prescan */
49 #define MINSTEP 4 /* minimum worthwhile preview step */
50
51 struct bound {int min,max;};
52
53 VIEW ourview = STDVIEW; /* desired view */
54 int hresolu = 512; /* horizontal resolution */
55 int vresolu = 512; /* vertical resolution */
56 double pixaspect = 1.0; /* pixel aspect ratio */
57
58 double zeps = .02; /* allowed z epsilon */
59
60 COLR *ourpict; /* output picture (COLR's) */
61 COLOR *ourspict; /* output pixel sums (averaging) */
62 float *ourweigh = NULL; /* output pixel weights (averaging) */
63 float *ourzbuf; /* corresponding z-buffer */
64 COLOR *ourbpict = NULL; /* blurred picture (view averaging) */
65
66 VIEW avgview; /* average view for -B option */
67 int nvavg; /* number of views averaged */
68
69 char *progname;
70
71 int fillo = F_FORE|F_BACK; /* selected fill options */
72 int fillsamp = 0; /* sample separation (0 == inf) */
73 extern int backfill(), rcalfill(); /* fill functions */
74 int (*fillfunc)() = backfill; /* selected fill function */
75 COLR backcolr = BLKCOLR; /* background color */
76 COLOR backcolor = BLKCOLOR; /* background color (float) */
77 double backz = 0.0; /* background z value */
78 int normdist = 1; /* i/o normalized distance? */
79 char ourfmt[LPICFMT+1] = PICFMT; /* original picture format */
80 double ourexp = -1; /* original picture exposure */
81 int expadj = 0; /* exposure adjustment (f-stops) */
82 double rexpadj = 1; /* real exposure adjustment */
83
84 VIEW theirview; /* input view */
85 int gotview; /* got input view? */
86 int wrongformat = 0; /* input in another format? */
87 RESOLU tresolu; /* input resolution */
88 double theirexp; /* input picture exposure */
89 MAT4 theirs2ours; /* transformation matrix */
90 int hasmatrix = 0; /* has transformation matrix */
91
92 int PDesc[3] = {-1,-1,-1}; /* rtrace process descriptor */
93 #define childpid (PDesc[2])
94 unsigned short queue[PACKSIZ][2]; /* pending pixels */
95 int packsiz; /* actual packet size */
96 int queuesiz = 0; /* number of pixels pending */
97
98 extern double movepixel();
99
100
101 main(argc, argv) /* interpolate pictures */
102 int argc;
103 char *argv[];
104 {
105 #define check(ol,al) if (argv[an][ol] || \
106 badarg(argc-an-1,argv+an+1,al)) \
107 goto badopt
108 int gotvfile = 0;
109 int doavg = -1;
110 int doblur = 0;
111 char *zfile = NULL;
112 char *expcomp = NULL;
113 int i, an, rval;
114
115 progname = argv[0];
116
117 for (an = 1; an < argc && argv[an][0] == '-'; an++) {
118 rval = getviewopt(&ourview, argc-an, argv+an);
119 if (rval >= 0) {
120 an += rval;
121 continue;
122 }
123 switch (argv[an][1]) {
124 case 'e': /* exposure */
125 check(2,"f");
126 expcomp = argv[++an];
127 break;
128 case 't': /* threshold */
129 check(2,"f");
130 zeps = atof(argv[++an]);
131 break;
132 case 'a': /* average */
133 check(2,NULL);
134 doavg = 1;
135 break;
136 case 'B': /* blur views */
137 check(2,NULL);
138 doblur = 1;
139 break;
140 case 'q': /* quick (no avg.) */
141 check(2,NULL);
142 doavg = 0;
143 break;
144 case 'n': /* dist. normalized? */
145 check(2,NULL);
146 normdist = !normdist;
147 break;
148 case 'f': /* fill type */
149 switch (argv[an][2]) {
150 case '0': /* none */
151 check(3,NULL);
152 fillo = 0;
153 break;
154 case 'f': /* foreground */
155 check(3,NULL);
156 fillo = F_FORE;
157 break;
158 case 'b': /* background */
159 check(3,NULL);
160 fillo = F_BACK;
161 break;
162 case 'a': /* all */
163 check(3,NULL);
164 fillo = F_FORE|F_BACK;
165 break;
166 case 's': /* sample */
167 check(3,"i");
168 fillsamp = atoi(argv[++an]);
169 break;
170 case 'c': /* color */
171 check(3,"fff");
172 fillfunc = backfill;
173 setcolor(backcolor, atof(argv[an+1]),
174 atof(argv[an+2]), atof(argv[an+3]));
175 setcolr(backcolr, colval(backcolor,RED),
176 colval(backcolor,GRN),
177 colval(backcolor,BLU));
178 an += 3;
179 break;
180 case 'z': /* z value */
181 check(3,"f");
182 fillfunc = backfill;
183 backz = atof(argv[++an]);
184 break;
185 case 'r': /* rtrace */
186 check(3,"s");
187 fillfunc = rcalfill;
188 calstart(RTCOM, argv[++an]);
189 break;
190 default:
191 goto badopt;
192 }
193 break;
194 case 'z': /* z file */
195 check(2,"s");
196 zfile = argv[++an];
197 break;
198 case 'x': /* x resolution */
199 check(2,"i");
200 hresolu = atoi(argv[++an]);
201 break;
202 case 'y': /* y resolution */
203 check(2,"i");
204 vresolu = atoi(argv[++an]);
205 break;
206 case 'p': /* pixel aspect */
207 if (argv[an][2] != 'a')
208 goto badopt;
209 check(3,"f");
210 pixaspect = atof(argv[++an]);
211 break;
212 case 'v': /* view file */
213 if (argv[an][2] != 'f')
214 goto badopt;
215 check(3,"s");
216 gotvfile = viewfile(argv[++an], &ourview, 0, 0);
217 if (gotvfile < 0)
218 syserror(argv[an]);
219 else if (gotvfile == 0) {
220 fprintf(stderr, "%s: bad view file\n",
221 argv[an]);
222 exit(1);
223 }
224 break;
225 default:
226 badopt:
227 fprintf(stderr, "%s: command line error at '%s'\n",
228 progname, argv[an]);
229 goto userr;
230 }
231 }
232 /* check arguments */
233 if ((argc-an)%2)
234 goto userr;
235 if (fillsamp == 1)
236 fillo &= ~F_BACK;
237 if (doavg < 0)
238 doavg = (argc-an) > 2;
239 if (expcomp != NULL)
240 if (expcomp[0] == '+' | expcomp[0] == '-') {
241 expadj = atof(expcomp) + (expcomp[0]=='+' ? .5 : -.5);
242 if (doavg | doblur)
243 rexpadj = pow(2.0, atof(expcomp));
244 else
245 rexpadj = pow(2.0, (double)expadj);
246 } else {
247 if (!isflt(expcomp))
248 goto userr;
249 rexpadj = atof(expcomp);
250 expadj = log(rexpadj)/LOG2 + (rexpadj>1 ? .5 : -.5);
251 if (!(doavg | doblur))
252 rexpadj = pow(2.0, (double)expadj);
253 }
254 /* set view */
255 if (nextview(doblur ? stdin : (FILE *)NULL) == EOF) {
256 fprintf(stderr, "%s: no view on standard input!\n",
257 progname);
258 exit(1);
259 }
260 normaspect(viewaspect(&ourview), &pixaspect, &hresolu, &vresolu);
261 /* allocate frame */
262 if (doavg) {
263 ourspict = (COLOR *)bmalloc(hresolu*vresolu*sizeof(COLOR));
264 ourweigh = (float *)bmalloc(hresolu*vresolu*sizeof(float));
265 if (ourspict == NULL | ourweigh == NULL)
266 syserror(progname);
267 } else {
268 ourpict = (COLR *)bmalloc(hresolu*vresolu*sizeof(COLR));
269 if (ourpict == NULL)
270 syserror(progname);
271 }
272 if (doblur) {
273 ourbpict = (COLOR *)bmalloc(hresolu*vresolu*sizeof(COLOR));
274 if (ourbpict == NULL)
275 syserror(progname);
276 }
277 ourzbuf = (float *)bmalloc(hresolu*vresolu*sizeof(float));
278 if (ourzbuf == NULL)
279 syserror(progname);
280 /* new header */
281 newheader("RADIANCE", stdout);
282 /* run pictures */
283 do {
284 bzero((char *)ourzbuf, hresolu*vresolu*sizeof(float));
285 for (i = an; i < argc; i += 2)
286 addpicture(argv[i], argv[i+1]);
287 if (fillo&F_BACK) /* fill in spaces */
288 backpicture(fillfunc, fillsamp);
289 else
290 fillpicture(fillfunc);
291 /* aft clipping */
292 clipaft();
293 } while (addblur() && nextview(stdin) != EOF);
294 /* close calculation */
295 caldone();
296 /* add to header */
297 printargs(argc, argv, stdout);
298 compavgview();
299 if (doblur | gotvfile) {
300 fputs(VIEWSTR, stdout);
301 fprintview(&avgview, stdout);
302 putc('\n', stdout);
303 }
304 if (pixaspect < .99 | pixaspect > 1.01)
305 fputaspect(pixaspect, stdout);
306 if (ourexp > 0)
307 ourexp *= rexpadj;
308 else
309 ourexp = rexpadj;
310 if (ourexp < .995 | ourexp > 1.005)
311 fputexpos(ourexp, stdout);
312 if (strcmp(ourfmt, PICFMT)) /* print format if known */
313 fputformat(ourfmt, stdout);
314 putc('\n', stdout);
315 /* write picture */
316 writepicture();
317 /* write z file */
318 if (zfile != NULL)
319 writedistance(zfile);
320
321 exit(0);
322 userr:
323 fprintf(stderr,
324 "Usage: %s [view opts][-t eps][-z zout][-e spec][-B][-a|-q][-fT][-n] pfile zspec ..\n",
325 progname);
326 exit(1);
327 #undef check
328 }
329
330
331 headline(s) /* process header string */
332 char *s;
333 {
334 char fmt[32];
335
336 if (isheadid(s))
337 return;
338 if (formatval(fmt, s)) {
339 if (globmatch(ourfmt, fmt)) {
340 wrongformat = 0;
341 strcpy(ourfmt, fmt);
342 } else
343 wrongformat = 1;
344 return;
345 }
346 if (nvavg < 2) {
347 putc('\t', stdout);
348 fputs(s, stdout);
349 }
350 if (isexpos(s)) {
351 theirexp *= exposval(s);
352 return;
353 }
354 if (isview(s) && sscanview(&theirview, s) > 0)
355 gotview++;
356 }
357
358
359 nextview(fp) /* get and set next view */
360 FILE *fp;
361 {
362 char linebuf[256];
363 char *err;
364 register int i;
365
366 if (fp != NULL) {
367 do /* get new view */
368 if (fgets(linebuf, sizeof(linebuf), fp) == NULL)
369 return(EOF);
370 while (!isview(linebuf) || !sscanview(&ourview, linebuf));
371 }
372 /* set new view */
373 if ((err = setview(&ourview)) != NULL) {
374 fprintf(stderr, "%s: %s\n", progname, err);
375 exit(1);
376 }
377 if (!nvavg) { /* first view */
378 copystruct(&avgview, &ourview);
379 return(nvavg++);
380 }
381 /* add to average view */
382 for (i = 0; i < 3; i++) {
383 avgview.vp[i] += ourview.vp[i];
384 avgview.vdir[i] += ourview.vdir[i];
385 avgview.vup[i] += ourview.vup[i];
386 }
387 avgview.horiz += ourview.horiz;
388 avgview.vert += ourview.vert;
389 avgview.hoff += ourview.hoff;
390 avgview.voff += ourview.voff;
391 avgview.vfore += ourview.vfore;
392 avgview.vaft += ourview.vaft;
393 return(nvavg++);
394 }
395
396
397 compavgview() /* compute average view */
398 {
399 register int i;
400 double f;
401
402 if (nvavg < 2)
403 return;
404 f = 1.0/nvavg;
405 for (i = 0; i < 3; i++) {
406 avgview.vp[i] *= f;
407 avgview.vdir[i] *= f;
408 avgview.vup[i] *= f;
409 }
410 avgview.horiz *= f;
411 avgview.vert *= f;
412 avgview.hoff *= f;
413 avgview.voff *= f;
414 avgview.vfore *= f;
415 avgview.vaft *= f;
416 if (setview(&avgview) != NULL) /* in case of emergency... */
417 copystruct(&avgview, &ourview);
418 pixaspect = viewaspect(&avgview) * hresolu / vresolu;
419 }
420
421
422 addpicture(pfile, zspec) /* add picture to output */
423 char *pfile, *zspec;
424 {
425 FILE *pfp;
426 int zfd;
427 char *err;
428 COLR *scanin;
429 float *zin;
430 struct position *plast;
431 struct bound *xlim, ylim;
432 int y;
433 /* open picture file */
434 if ((pfp = fopen(pfile, "r")) == NULL)
435 syserror(pfile);
436 /* get header with exposure and view */
437 theirexp = 1.0;
438 copystruct(&theirview, &stdview);
439 gotview = 0;
440 if (nvavg < 2)
441 printf("%s:\n", pfile);
442 getheader(pfp, headline, NULL);
443 if (wrongformat || !gotview || !fgetsresolu(&tresolu, pfp)) {
444 fprintf(stderr, "%s: picture format error\n", pfile);
445 exit(1);
446 }
447 if (ourexp <= 0)
448 ourexp = theirexp;
449 else if (ABS(theirexp-ourexp) > .01*ourexp)
450 fprintf(stderr, "%s: different exposure (warning)\n", pfile);
451 if (err = setview(&theirview)) {
452 fprintf(stderr, "%s: %s\n", pfile, err);
453 exit(1);
454 }
455 /* compute transformation */
456 hasmatrix = pixform(theirs2ours, &theirview, &ourview);
457 /* get z specification or file */
458 zin = (float *)malloc(scanlen(&tresolu)*sizeof(float));
459 if (zin == NULL)
460 syserror(progname);
461 if ((zfd = open(zspec, O_RDONLY)) == -1) {
462 double zvalue;
463 register int x;
464 if (!isflt(zspec) || (zvalue = atof(zspec)) <= 0.0)
465 syserror(zspec);
466 for (x = scanlen(&tresolu); x-- > 0; )
467 zin[x] = zvalue;
468 }
469 /* compute transferrable perimeter */
470 xlim = (struct bound *)malloc(numscans(&tresolu)*sizeof(struct bound));
471 if (xlim == NULL)
472 syserror(progname);
473 if (!getperim(xlim, &ylim, zin, zfd)) { /* overlapping area? */
474 free((char *)zin);
475 free((char *)xlim);
476 if (zfd != -1)
477 close(zfd);
478 fclose(pfp);
479 return;
480 }
481 /* allocate scanlines */
482 scanin = (COLR *)malloc(scanlen(&tresolu)*sizeof(COLR));
483 plast = (struct position *)calloc(scanlen(&tresolu),
484 sizeof(struct position));
485 if (scanin == NULL | plast == NULL)
486 syserror(progname);
487 /* skip to starting point */
488 for (y = 0; y < ylim.min; y++)
489 if (freadcolrs(scanin, scanlen(&tresolu), pfp) < 0) {
490 fprintf(stderr, "%s: read error\n", pfile);
491 exit(1);
492 }
493 if (zfd != -1 && lseek(zfd,
494 (long)ylim.min*scanlen(&tresolu)*sizeof(float), 0) < 0)
495 syserror(zspec);
496 /* load image */
497 for (y = ylim.min; y <= ylim.max; y++) {
498 if (freadcolrs(scanin, scanlen(&tresolu), pfp) < 0) {
499 fprintf(stderr, "%s: read error\n", pfile);
500 exit(1);
501 }
502 if (zfd != -1 && read(zfd, (char *)zin,
503 scanlen(&tresolu)*sizeof(float))
504 < scanlen(&tresolu)*sizeof(float))
505 syserror(zspec);
506 addscanline(xlim+y, y, scanin, zin, plast);
507 }
508 /* clean up */
509 free((char *)xlim);
510 free((char *)scanin);
511 free((char *)zin);
512 free((char *)plast);
513 fclose(pfp);
514 if (zfd != -1)
515 close(zfd);
516 }
517
518
519 pixform(xfmat, vw1, vw2) /* compute view1 to view2 matrix */
520 register MAT4 xfmat;
521 register VIEW *vw1, *vw2;
522 {
523 double m4t[4][4];
524
525 if (vw1->type != VT_PER & vw1->type != VT_PAR)
526 return(0);
527 if (vw2->type != VT_PER & vw2->type != VT_PAR)
528 return(0);
529 setident4(xfmat);
530 xfmat[0][0] = vw1->hvec[0];
531 xfmat[0][1] = vw1->hvec[1];
532 xfmat[0][2] = vw1->hvec[2];
533 xfmat[1][0] = vw1->vvec[0];
534 xfmat[1][1] = vw1->vvec[1];
535 xfmat[1][2] = vw1->vvec[2];
536 xfmat[2][0] = vw1->vdir[0];
537 xfmat[2][1] = vw1->vdir[1];
538 xfmat[2][2] = vw1->vdir[2];
539 xfmat[3][0] = vw1->vp[0];
540 xfmat[3][1] = vw1->vp[1];
541 xfmat[3][2] = vw1->vp[2];
542 setident4(m4t);
543 m4t[0][0] = vw2->hvec[0]/vw2->hn2;
544 m4t[1][0] = vw2->hvec[1]/vw2->hn2;
545 m4t[2][0] = vw2->hvec[2]/vw2->hn2;
546 m4t[3][0] = -DOT(vw2->vp,vw2->hvec)/vw2->hn2;
547 m4t[0][1] = vw2->vvec[0]/vw2->vn2;
548 m4t[1][1] = vw2->vvec[1]/vw2->vn2;
549 m4t[2][1] = vw2->vvec[2]/vw2->vn2;
550 m4t[3][1] = -DOT(vw2->vp,vw2->vvec)/vw2->vn2;
551 m4t[0][2] = vw2->vdir[0];
552 m4t[1][2] = vw2->vdir[1];
553 m4t[2][2] = vw2->vdir[2];
554 m4t[3][2] = -DOT(vw2->vp,vw2->vdir);
555 multmat4(xfmat, xfmat, m4t);
556 return(1);
557 }
558
559
560 addscanline(xl, y, pline, zline, lasty) /* add scanline to output */
561 struct bound *xl;
562 int y;
563 COLR *pline;
564 float *zline;
565 struct position *lasty; /* input/output */
566 {
567 FVECT pos;
568 struct position lastx, newpos;
569 double wt;
570 register int x;
571
572 lastx.z = 0;
573 for (x = xl->max; x >= xl->min; x--) {
574 pix2loc(pos, &tresolu, x, y);
575 pos[2] = zline[x];
576 if ((wt = movepixel(pos)) <= FTINY) {
577 lasty[x].z = lastx.z = 0; /* mark invalid */
578 continue;
579 }
580 /* add pixel to our image */
581 newpos.x = pos[0] * hresolu;
582 newpos.y = pos[1] * vresolu;
583 newpos.z = zline[x];
584 addpixel(&newpos, &lastx, &lasty[x], pline[x], wt, pos[2]);
585 lasty[x].x = lastx.x = newpos.x;
586 lasty[x].y = lastx.y = newpos.y;
587 lasty[x].z = lastx.z = newpos.z;
588 }
589 }
590
591
592 addpixel(p0, p1, p2, pix, w, z) /* fill in pixel parallelogram */
593 struct position *p0, *p1, *p2;
594 COLR pix;
595 double w;
596 double z;
597 {
598 double zt = 2.*zeps*p0->z; /* threshold */
599 COLOR pval; /* converted+weighted pixel */
600 int s1x, s1y, s2x, s2y; /* step sizes */
601 int l1, l2, c1, c2; /* side lengths and counters */
602 int p1isy; /* p0p1 along y? */
603 int x1, y1; /* p1 position */
604 register int x, y; /* final position */
605
606 /* compute vector p0p1 */
607 if (fillo&F_FORE && ABS(p1->z-p0->z) <= zt) {
608 s1x = p1->x - p0->x;
609 s1y = p1->y - p0->y;
610 l1 = ABS(s1x);
611 if (p1isy = (ABS(s1y) > l1))
612 l1 = ABS(s1y);
613 else if (l1 < 1)
614 l1 = 1;
615 } else {
616 l1 = s1x = s1y = 1;
617 p1isy = -1;
618 }
619 /* compute vector p0p2 */
620 if (fillo&F_FORE && ABS(p2->z-p0->z) <= zt) {
621 s2x = p2->x - p0->x;
622 s2y = p2->y - p0->y;
623 if (p1isy == 1)
624 l2 = ABS(s2x);
625 else {
626 l2 = ABS(s2y);
627 if (p1isy != 0 && ABS(s2x) > l2)
628 l2 = ABS(s2x);
629 }
630 if (l2 < 1)
631 l2 = 1;
632 } else
633 l2 = s2x = s2y = 1;
634 /* fill the parallelogram */
635 if (averaging) {
636 colr_color(pval, pix);
637 scalecolor(pval, w);
638 }
639 for (c1 = l1; c1-- > 0; ) {
640 x1 = p0->x + c1*s1x/l1;
641 y1 = p0->y + c1*s1y/l1;
642 for (c2 = l2; c2-- > 0; ) {
643 x = x1 + c2*s2x/l2;
644 if (x < 0 | x >= hresolu)
645 continue;
646 y = y1 + c2*s2y/l2;
647 if (y < 0 | y >= vresolu)
648 continue;
649 if (averaging) {
650 if (zscan(y)[x] <= 0 || zscan(y)[x]-z
651 > zeps*zscan(y)[x]) {
652 copycolor(sscan(y)[x], pval);
653 wscan(y)[x] = w;
654 zscan(y)[x] = z;
655 } else if (z-zscan(y)[x] <= zeps*zscan(y)[x]) {
656 addcolor(sscan(y)[x], pval);
657 wscan(y)[x] += w;
658 }
659 } else if (zscan(y)[x] <= 0 || zscan(y)[x]-z
660 > zeps*zscan(y)[x]) {
661 copycolr(pscan(y)[x], pix);
662 zscan(y)[x] = z;
663 }
664 }
665 }
666 }
667
668
669 double
670 movepixel(pos) /* reposition image point */
671 register FVECT pos;
672 {
673 FVECT pt, tdir, odir;
674 double d;
675 register int i;
676
677 if (pos[2] <= 0) /* empty pixel */
678 return(0);
679 if (usematrix) {
680 pos[0] += theirview.hoff - .5;
681 pos[1] += theirview.voff - .5;
682 if (normdist & theirview.type == VT_PER)
683 d = sqrt(1. + pos[0]*pos[0]*theirview.hn2
684 + pos[1]*pos[1]*theirview.vn2);
685 else
686 d = 1.;
687 pos[2] += d*theirview.vfore;
688 if (theirview.type == VT_PER) {
689 pos[2] /= d;
690 pos[0] *= pos[2];
691 pos[1] *= pos[2];
692 }
693 multp3(pos, pos, theirs2ours);
694 if (pos[2] <= ourview.vfore)
695 return(0);
696 if (ourview.type == VT_PER) {
697 pos[0] /= pos[2];
698 pos[1] /= pos[2];
699 }
700 pos[0] += .5 - ourview.hoff;
701 pos[1] += .5 - ourview.voff;
702 pos[2] -= ourview.vfore;
703 } else {
704 if (viewray(pt, tdir, &theirview, pos[0], pos[1]) < -FTINY)
705 return(0);
706 if (!normdist & theirview.type == VT_PER) /* adjust */
707 pos[2] *= sqrt(1. + pos[0]*pos[0]*theirview.hn2
708 + pos[1]*pos[1]*theirview.vn2);
709 pt[0] += tdir[0]*pos[2];
710 pt[1] += tdir[1]*pos[2];
711 pt[2] += tdir[2]*pos[2];
712 viewloc(pos, &ourview, pt);
713 if (pos[2] <= 0)
714 return(0);
715 }
716 if (pos[0] < 0 | pos[0] >= 1-FTINY | pos[1] < 0 | pos[1] >= 1-FTINY)
717 return(0);
718 if (!averaging)
719 return(1);
720 if (ourview.type == VT_PAR) /* compute our direction */
721 VCOPY(odir, ourview.vdir);
722 else
723 for (i = 0; i < 3; i++)
724 odir[i] = (pt[i] - ourview.vp[i])/pos[2];
725 d = DOT(odir,tdir); /* compute pixel weight */
726 if (d >= 1.-1./MAXWT/MAXWT)
727 return(MAXWT); /* clip to maximum weight */
728 return(1./sqrt(1.-d));
729 }
730
731
732 getperim(xl, yl, zline, zfd) /* compute overlapping image area */
733 register struct bound *xl;
734 struct bound *yl;
735 float *zline;
736 int zfd;
737 {
738 int step;
739 FVECT pos;
740 register int x, y;
741 /* set up step size */
742 if (scanlen(&tresolu) < numscans(&tresolu))
743 step = scanlen(&tresolu)/NSTEPS;
744 else
745 step = numscans(&tresolu)/NSTEPS;
746 if (step < MINSTEP) { /* not worth cropping? */
747 yl->min = 0;
748 yl->max = numscans(&tresolu) - 1;
749 x = scanlen(&tresolu) - 1;
750 for (y = numscans(&tresolu); y--; ) {
751 xl[y].min = 0;
752 xl[y].max = x;
753 }
754 return(1);
755 }
756 yl->min = 32000; yl->max = 0; /* search for points on image */
757 for (y = step - 1; y < numscans(&tresolu); y += step) {
758 if (zfd != -1) {
759 if (lseek(zfd, (long)y*scanlen(&tresolu)*sizeof(float),
760 0) < 0)
761 syserror("lseek");
762 if (read(zfd, (char *)zline,
763 scanlen(&tresolu)*sizeof(float))
764 < scanlen(&tresolu)*sizeof(float))
765 syserror("read");
766 }
767 xl[y].min = 32000; xl[y].max = 0; /* x max */
768 for (x = scanlen(&tresolu); (x -= step) > 0; ) {
769 pix2loc(pos, &tresolu, x, y);
770 pos[2] = zline[x];
771 if (movepixel(pos) > FTINY) {
772 xl[y].max = x + step - 1;
773 xl[y].min = x - step + 1; /* x min */
774 if (xl[y].min < 0)
775 xl[y].min = 0;
776 for (x = step - 1; x < xl[y].max; x += step) {
777 pix2loc(pos, &tresolu, x, y);
778 pos[2] = zline[x];
779 if (movepixel(pos) > FTINY) {
780 xl[y].min = x - step + 1;
781 break;
782 }
783 }
784 if (y < yl->min) /* y limits */
785 yl->min = y - step + 1;
786 yl->max = y + step - 1;
787 break;
788 }
789 }
790 /* fill in between */
791 if (y < step) {
792 xl[y-1].min = xl[y].min;
793 xl[y-1].max = xl[y].max;
794 } else {
795 if (xl[y].min < xl[y-step].min)
796 xl[y-1].min = xl[y].min;
797 else
798 xl[y-1].min = xl[y-step].min;
799 if (xl[y].max > xl[y-step].max)
800 xl[y-1].max = xl[y].max;
801 else
802 xl[y-1].max = xl[y-step].max;
803 }
804 for (x = 2; x < step; x++)
805 copystruct(xl+y-x, xl+y-1);
806 }
807 if (yl->max >= numscans(&tresolu))
808 yl->max = numscans(&tresolu) - 1;
809 y -= step;
810 for (x = numscans(&tresolu) - 1; x > y; x--) /* fill bottom rows */
811 copystruct(xl+x, xl+y);
812 return(yl->max >= yl->min);
813 }
814
815
816 backpicture(fill, samp) /* background fill algorithm */
817 int (*fill)();
818 int samp;
819 {
820 int *yback, xback;
821 int y;
822 register int x, i;
823 /* get back buffer */
824 yback = (int *)malloc(hresolu*sizeof(int));
825 if (yback == NULL)
826 syserror(progname);
827 for (x = 0; x < hresolu; x++)
828 yback[x] = -2;
829 /*
830 * Xback and yback are the pixel locations of suitable
831 * background values in each direction.
832 * A value of -2 means unassigned, and -1 means
833 * that there is no suitable background in this direction.
834 */
835 /* fill image */
836 for (y = 0; y < vresolu; y++) {
837 xback = -2;
838 for (x = 0; x < hresolu; x++)
839 if (zscan(y)[x] <= 0) { /* empty pixel */
840 /*
841 * First, find background from above or below.
842 * (farthest assigned pixel)
843 */
844 if (yback[x] == -2) {
845 for (i = y+1; i < vresolu; i++)
846 if (zscan(i)[x] > 0)
847 break;
848 if (i < vresolu
849 && (y <= 0 || zscan(y-1)[x] < zscan(i)[x]))
850 yback[x] = i;
851 else
852 yback[x] = y-1;
853 }
854 /*
855 * Next, find background from left or right.
856 */
857 if (xback == -2) {
858 for (i = x+1; i < hresolu; i++)
859 if (zscan(y)[i] > 0)
860 break;
861 if (i < hresolu
862 && (x <= 0 || zscan(y)[x-1] < zscan(y)[i]))
863 xback = i;
864 else
865 xback = x-1;
866 }
867 /*
868 * If we have no background for this pixel,
869 * use the given fill function.
870 */
871 if (xback < 0 && yback[x] < 0)
872 goto fillit;
873 /*
874 * Compare, and use the background that is
875 * farther, unless one of them is next to us.
876 * If the background is too distant, call
877 * the fill function.
878 */
879 if ( yback[x] < 0
880 || (xback >= 0 && ABS(x-xback) <= 1)
881 || ( ABS(y-yback[x]) > 1
882 && zscan(yback[x])[x]
883 < zscan(y)[xback] ) ) {
884 if (samp > 0 && ABS(x-xback) >= samp)
885 goto fillit;
886 if (averaging) {
887 copycolor(sscan(y)[x],
888 sscan(y)[xback]);
889 wscan(y)[x] = wscan(y)[xback];
890 } else
891 copycolr(pscan(y)[x],
892 pscan(y)[xback]);
893 zscan(y)[x] = zscan(y)[xback];
894 } else {
895 if (samp > 0 && ABS(y-yback[x]) > samp)
896 goto fillit;
897 if (averaging) {
898 copycolor(sscan(y)[x],
899 sscan(yback[x])[x]);
900 wscan(y)[x] =
901 wscan(yback[x])[x];
902 } else
903 copycolr(pscan(y)[x],
904 pscan(yback[x])[x]);
905 zscan(y)[x] = zscan(yback[x])[x];
906 }
907 continue;
908 fillit:
909 (*fill)(x,y);
910 if (fill == rcalfill) { /* use it */
911 clearqueue();
912 xback = x;
913 yback[x] = y;
914 }
915 } else { /* full pixel */
916 yback[x] = -2;
917 xback = -2;
918 }
919 }
920 free((char *)yback);
921 }
922
923
924 fillpicture(fill) /* paint in empty pixels using fill */
925 int (*fill)();
926 {
927 register int x, y;
928
929 for (y = 0; y < vresolu; y++)
930 for (x = 0; x < hresolu; x++)
931 if (zscan(y)[x] <= 0)
932 (*fill)(x,y);
933 if (fill == rcalfill)
934 clearqueue();
935 }
936
937
938 clipaft() /* perform aft clipping as indicated */
939 {
940 register int x, y;
941 int adjtest = ourview.type == VT_PER & zisnorm;
942 double tstdist;
943 double yzn2, vx;
944
945 if (ourview.vaft <= FTINY)
946 return(0);
947 tstdist = ourview.vaft - ourview.vfore;
948 for (y = 0; y < vresolu; y++) {
949 if (adjtest) { /* adjust test */
950 yzn2 = (y+.5)/vresolu + ourview.voff - .5;
951 yzn2 = 1. + yzn2*yzn2*ourview.vn2;
952 tstdist = (ourview.vaft - ourview.vfore)*sqrt(yzn2);
953 }
954 for (x = 0; x < hresolu; x++)
955 if (zscan(y)[x] > tstdist) {
956 if (adjtest) {
957 vx = (x+.5)/hresolu + ourview.hoff - .5;
958 if (zscan(y)[x] <= (ourview.vaft -
959 ourview.vfore) *
960 sqrt(vx*vx*ourview.hn2 + yzn2))
961 continue;
962 }
963 if (averaging)
964 bzero(sscan(y)[x], sizeof(COLOR));
965 else
966 bzero(pscan(y)[x], sizeof(COLR));
967 zscan(y)[x] = 0.0;
968 }
969 }
970 return(1);
971 }
972
973
974 addblur() /* add to blurred picture */
975 {
976 COLOR cval;
977 double d;
978 register int i;
979
980 if (!blurring)
981 return(0);
982 i = hresolu*vresolu;
983 if (nvavg < 2)
984 if (averaging)
985 while (i--) {
986 copycolor(ourbpict[i], ourspict[i]);
987 d = 1.0/ourweigh[i];
988 scalecolor(ourbpict[i], d);
989 }
990 else
991 while (i--)
992 colr_color(ourbpict[i], ourpict[i]);
993 else
994 if (averaging)
995 while (i--) {
996 copycolor(cval, ourspict[i]);
997 d = 1.0/ourweigh[i];
998 scalecolor(cval, d);
999 addcolor(ourbpict[i], cval);
1000 }
1001 else
1002 while (i--) {
1003 colr_color(cval, ourpict[i]);
1004 addcolor(ourbpict[i], cval);
1005 }
1006 /* print view */
1007 printf("VIEW%d:", nvavg);
1008 fprintview(&ourview, stdout);
1009 putchar('\n');
1010 return(1);
1011 }
1012
1013
1014 writepicture() /* write out picture (alters buffer) */
1015 {
1016 int y;
1017 register int x;
1018 double d;
1019
1020 fprtresolu(hresolu, vresolu, stdout);
1021 for (y = vresolu-1; y >= 0; y--)
1022 if (blurring) {
1023 for (x = 0; x < hresolu; x++) { /* compute avg. */
1024 d = rexpadj/nvavg;
1025 scalecolor(bscan(y)[x], d);
1026 }
1027 if (fwritescan(bscan(y), hresolu, stdout) < 0)
1028 syserror(progname);
1029 } else if (averaging) {
1030 for (x = 0; x < hresolu; x++) { /* average pixels */
1031 d = rexpadj/wscan(y)[x];
1032 scalecolor(sscan(y)[x], d);
1033 }
1034 if (fwritescan(sscan(y), hresolu, stdout) < 0)
1035 syserror(progname);
1036 } else {
1037 if (expadj)
1038 shiftcolrs(pscan(y), hresolu, expadj);
1039 if (fwritecolrs(pscan(y), hresolu, stdout) < 0)
1040 syserror(progname);
1041 }
1042 }
1043
1044
1045 writedistance(fname) /* write out z file (alters buffer) */
1046 char *fname;
1047 {
1048 int donorm = normdist & !zisnorm ? 1 :
1049 ourview.type == VT_PER & !normdist & zisnorm ? -1 : 0;
1050 int fd;
1051 int y;
1052
1053 if ((fd = open(fname, O_WRONLY|O_CREAT|O_TRUNC, 0666)) == -1)
1054 syserror(fname);
1055 for (y = vresolu-1; y >= 0; y--) {
1056 if (donorm) {
1057 double vx, yzn2, d;
1058 register int x;
1059 yzn2 = (y+.5)/vresolu + ourview.voff - .5;
1060 yzn2 = 1. + yzn2*yzn2*ourview.vn2;
1061 for (x = 0; x < hresolu; x++) {
1062 vx = (x+.5)/hresolu + ourview.hoff - .5;
1063 d = sqrt(vx*vx*ourview.hn2 + yzn2);
1064 if (donorm > 0)
1065 zscan(y)[x] *= d;
1066 else
1067 zscan(y)[x] /= d;
1068 }
1069 }
1070 if (write(fd, (char *)zscan(y), hresolu*sizeof(float))
1071 < hresolu*sizeof(float))
1072 syserror(fname);
1073 }
1074 close(fd);
1075 }
1076
1077
1078 backfill(x, y) /* fill pixel with background */
1079 int x, y;
1080 {
1081 if (averaging) {
1082 copycolor(sscan(y)[x], backcolor);
1083 wscan(y)[x] = 1;
1084 } else
1085 copycolr(pscan(y)[x], backcolr);
1086 zscan(y)[x] = backz;
1087 }
1088
1089
1090 calstart(prog, args) /* start fill calculation */
1091 char *prog, *args;
1092 {
1093 char combuf[512];
1094 char *argv[64];
1095 int rval;
1096 register char **wp, *cp;
1097
1098 if (childpid != -1) {
1099 fprintf(stderr, "%s: too many calculations\n", progname);
1100 exit(1);
1101 }
1102 strcpy(combuf, prog);
1103 strcat(combuf, args);
1104 cp = combuf;
1105 wp = argv;
1106 for ( ; ; ) {
1107 while (isspace(*cp)) /* nullify spaces */
1108 *cp++ = '\0';
1109 if (!*cp) /* all done? */
1110 break;
1111 *wp++ = cp; /* add argument to list */
1112 while (*++cp && !isspace(*cp))
1113 ;
1114 }
1115 *wp = NULL;
1116 /* start process */
1117 if ((rval = open_process(PDesc, argv)) < 0)
1118 syserror(progname);
1119 if (rval == 0) {
1120 fprintf(stderr, "%s: command not found\n", argv[0]);
1121 exit(1);
1122 }
1123 packsiz = rval/(6*sizeof(float)) - 1;
1124 if (packsiz > PACKSIZ)
1125 packsiz = PACKSIZ;
1126 queuesiz = 0;
1127 }
1128
1129
1130 caldone() /* done with calculation */
1131 {
1132 if (childpid == -1)
1133 return;
1134 clearqueue();
1135 close_process(PDesc);
1136 childpid = -1;
1137 }
1138
1139
1140 rcalfill(x, y) /* fill with ray-calculated pixel */
1141 int x, y;
1142 {
1143 if (queuesiz >= packsiz) /* flush queue if needed */
1144 clearqueue();
1145 /* add position to queue */
1146 queue[queuesiz][0] = x;
1147 queue[queuesiz][1] = y;
1148 queuesiz++;
1149 }
1150
1151
1152 clearqueue() /* process queue */
1153 {
1154 FVECT orig, dir;
1155 float fbuf[6*(PACKSIZ+1)];
1156 register float *fbp;
1157 register int i;
1158 double vx, vy;
1159
1160 if (queuesiz == 0)
1161 return;
1162 fbp = fbuf;
1163 for (i = 0; i < queuesiz; i++) {
1164 viewray(orig, dir, &ourview,
1165 (queue[i][0]+.5)/hresolu,
1166 (queue[i][1]+.5)/vresolu);
1167 *fbp++ = orig[0]; *fbp++ = orig[1]; *fbp++ = orig[2];
1168 *fbp++ = dir[0]; *fbp++ = dir[1]; *fbp++ = dir[2];
1169 }
1170 /* mark end and get results */
1171 bzero((char *)fbp, 6*sizeof(float));
1172 if (process(PDesc, fbuf, fbuf, 4*sizeof(float)*(queuesiz+1),
1173 6*sizeof(float)*(queuesiz+1)) !=
1174 4*sizeof(float)*(queuesiz+1)) {
1175 fprintf(stderr, "%s: error reading from rtrace process\n",
1176 progname);
1177 exit(1);
1178 }
1179 fbp = fbuf;
1180 for (i = 0; i < queuesiz; i++) {
1181 if (ourexp > 0 && ourexp != 1.0) {
1182 fbp[0] *= ourexp;
1183 fbp[1] *= ourexp;
1184 fbp[2] *= ourexp;
1185 }
1186 if (averaging) {
1187 setcolor(sscan(queue[i][1])[queue[i][0]],
1188 fbp[0], fbp[1], fbp[2]);
1189 wscan(queue[i][1])[queue[i][0]] = 1;
1190 } else
1191 setcolr(pscan(queue[i][1])[queue[i][0]],
1192 fbp[0], fbp[1], fbp[2]);
1193 if (zisnorm)
1194 zscan(queue[i][1])[queue[i][0]] = fbp[3];
1195 else {
1196 vx = (queue[i][0]+.5)/hresolu + ourview.hoff - .5;
1197 vy = (queue[i][1]+.5)/vresolu + ourview.voff - .5;
1198 zscan(queue[i][1])[queue[i][0]] = fbp[3] / sqrt(1. +
1199 vx*vx*ourview.hn2 + vy*vy*ourview.vn2);
1200 }
1201 fbp += 4;
1202 }
1203 queuesiz = 0;
1204 }
1205
1206
1207 syserror(s) /* report error and exit */
1208 char *s;
1209 {
1210 perror(s);
1211 exit(1);
1212 }