ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.24
Committed: Thu Sep 14 13:40:17 1995 UTC (29 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.23: +187 -55 lines
Log Message:
added -B option for blurring images

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