ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.34
Committed: Thu Jun 26 00:58:10 2003 UTC (20 years, 10 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.33: +10 -11 lines
Log Message:
Abstracted process and path handling for Windows.
Renamed FLOAT to RREAL because of conflict on Windows.
Added conditional compiles for some signal handlers.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: pinterp.c,v 2.33 2003/02/22 02:07:27 greg Exp $";
3 #endif
4 /*
5 * Interpolate and extrapolate pictures with different view parameters.
6 *
7 * Greg Ward 09Dec89
8 */
9
10 #include "copyright.h"
11
12 #include <ctype.h>
13
14 #include "standard.h"
15 #include "rtprocess.h" /* Windows: must come before color.h */
16 #include "view.h"
17 #include "color.h"
18
19 #define LOG2 0.69314718055994530942
20
21 #define pscan(y) (ourpict+(y)*hresolu)
22 #define sscan(y) (ourspict+(y)*hresolu)
23 #define wscan(y) (ourweigh+(y)*hresolu)
24 #define zscan(y) (ourzbuf+(y)*hresolu)
25 #define bscan(y) (ourbpict+(y)*hresolu)
26 #define averaging (ourweigh != NULL)
27 #define blurring (ourbpict != NULL)
28 #define usematrix (hasmatrix & !averaging)
29 #define zisnorm (!usematrix | ourview.type != VT_PER)
30
31 #define MAXWT 1000. /* maximum pixel weight (averaging) */
32
33 #define F_FORE 1 /* fill foreground */
34 #define F_BACK 2 /* fill background */
35
36 #define PACKSIZ 256 /* max. calculation packet size */
37
38 #define RTCOM "rtrace -h- -ovl -fff -ld- -i- -I- "
39
40 #define ABS(x) ((x)>0?(x):-(x))
41
42 struct position {int x,y; float z;};
43
44 #define NSTEPS 64 /* number steps in overlap prescan */
45 #define MINSTEP 4 /* minimum worthwhile preview step */
46
47 struct bound {int min,max;};
48
49 VIEW ourview = STDVIEW; /* desired view */
50 int hresolu = 512; /* horizontal resolution */
51 int vresolu = 512; /* vertical resolution */
52 double pixaspect = 1.0; /* pixel aspect ratio */
53
54 double zeps = .02; /* allowed z epsilon */
55
56 COLR *ourpict; /* output picture (COLR's) */
57 COLOR *ourspict; /* output pixel sums (averaging) */
58 float *ourweigh = NULL; /* output pixel weights (averaging) */
59 float *ourzbuf; /* corresponding z-buffer */
60 COLOR *ourbpict = NULL; /* blurred picture (view averaging) */
61
62 VIEW avgview; /* average view for -B option */
63 int nvavg; /* number of views averaged */
64
65 char *progname;
66
67 int fillo = F_FORE|F_BACK; /* selected fill options */
68 int fillsamp = 0; /* sample separation (0 == inf) */
69 extern int backfill(), rcalfill(); /* fill functions */
70 int (*fillfunc)() = backfill; /* selected fill function */
71 COLR backcolr = BLKCOLR; /* background color */
72 COLOR backcolor = BLKCOLOR; /* background color (float) */
73 double backz = 0.0; /* background z value */
74 int normdist = 1; /* i/o normalized distance? */
75 char ourfmt[LPICFMT+1] = PICFMT; /* original picture format */
76 double ourexp = -1; /* original picture exposure */
77 int expadj = 0; /* exposure adjustment (f-stops) */
78 double rexpadj = 1; /* real exposure adjustment */
79
80 VIEW theirview; /* input view */
81 int gotview; /* got input view? */
82 int wrongformat = 0; /* input in another format? */
83 RESOLU tresolu; /* input resolution */
84 double theirexp; /* input picture exposure */
85 MAT4 theirs2ours; /* transformation matrix */
86 int hasmatrix = 0; /* has transformation matrix */
87
88 static SUBPROC PDesc = SP_INACTIVE; /* rtrace process descriptor */
89 unsigned short queue[PACKSIZ][2]; /* pending pixels */
90 int packsiz; /* actual packet size */
91 int queuesiz = 0; /* number of pixels pending */
92
93 extern double movepixel();
94
95
96 main(argc, argv) /* interpolate pictures */
97 int argc;
98 char *argv[];
99 {
100 #define check(ol,al) if (argv[an][ol] || \
101 badarg(argc-an-1,argv+an+1,al)) \
102 goto badopt
103 int gotvfile = 0;
104 int doavg = -1;
105 int doblur = 0;
106 char *zfile = NULL;
107 char *expcomp = NULL;
108 int i, an, rval;
109
110 progname = argv[0];
111
112 for (an = 1; an < argc && argv[an][0] == '-'; an++) {
113 rval = getviewopt(&ourview, argc-an, argv+an);
114 if (rval >= 0) {
115 an += rval;
116 continue;
117 }
118 switch (argv[an][1]) {
119 case 'e': /* exposure */
120 check(2,"f");
121 expcomp = argv[++an];
122 break;
123 case 't': /* threshold */
124 check(2,"f");
125 zeps = atof(argv[++an]);
126 break;
127 case 'a': /* average */
128 check(2,NULL);
129 doavg = 1;
130 break;
131 case 'B': /* blur views */
132 check(2,NULL);
133 doblur = 1;
134 break;
135 case 'q': /* quick (no avg.) */
136 check(2,NULL);
137 doavg = 0;
138 break;
139 case 'n': /* dist. normalized? */
140 check(2,NULL);
141 normdist = !normdist;
142 break;
143 case 'f': /* fill type */
144 switch (argv[an][2]) {
145 case '0': /* none */
146 check(3,NULL);
147 fillo = 0;
148 break;
149 case 'f': /* foreground */
150 check(3,NULL);
151 fillo = F_FORE;
152 break;
153 case 'b': /* background */
154 check(3,NULL);
155 fillo = F_BACK;
156 break;
157 case 'a': /* all */
158 check(3,NULL);
159 fillo = F_FORE|F_BACK;
160 break;
161 case 's': /* sample */
162 check(3,"i");
163 fillsamp = atoi(argv[++an]);
164 break;
165 case 'c': /* color */
166 check(3,"fff");
167 fillfunc = backfill;
168 setcolor(backcolor, atof(argv[an+1]),
169 atof(argv[an+2]), atof(argv[an+3]));
170 setcolr(backcolr, colval(backcolor,RED),
171 colval(backcolor,GRN),
172 colval(backcolor,BLU));
173 an += 3;
174 break;
175 case 'z': /* z value */
176 check(3,"f");
177 fillfunc = backfill;
178 backz = atof(argv[++an]);
179 break;
180 case 'r': /* rtrace */
181 check(3,"s");
182 fillfunc = rcalfill;
183 calstart(RTCOM, argv[++an]);
184 break;
185 default:
186 goto badopt;
187 }
188 break;
189 case 'z': /* z file */
190 check(2,"s");
191 zfile = argv[++an];
192 break;
193 case 'x': /* x resolution */
194 check(2,"i");
195 hresolu = atoi(argv[++an]);
196 break;
197 case 'y': /* y resolution */
198 check(2,"i");
199 vresolu = atoi(argv[++an]);
200 break;
201 case 'p': /* pixel aspect */
202 if (argv[an][2] != 'a')
203 goto badopt;
204 check(3,"f");
205 pixaspect = atof(argv[++an]);
206 break;
207 case 'v': /* view file */
208 if (argv[an][2] != 'f')
209 goto badopt;
210 check(3,"s");
211 gotvfile = viewfile(argv[++an], &ourview, NULL);
212 if (gotvfile < 0)
213 syserror(argv[an]);
214 else if (gotvfile == 0) {
215 fprintf(stderr, "%s: bad view file\n",
216 argv[an]);
217 exit(1);
218 }
219 break;
220 default:
221 badopt:
222 fprintf(stderr, "%s: command line error at '%s'\n",
223 progname, argv[an]);
224 goto userr;
225 }
226 }
227 /* check arguments */
228 if ((argc-an)%2)
229 goto userr;
230 if (fillsamp == 1)
231 fillo &= ~F_BACK;
232 if (doavg < 0)
233 doavg = (argc-an) > 2;
234 if (expcomp != NULL)
235 if (expcomp[0] == '+' | expcomp[0] == '-') {
236 expadj = atof(expcomp) + (expcomp[0]=='+' ? .5 : -.5);
237 if (doavg | doblur)
238 rexpadj = pow(2.0, atof(expcomp));
239 else
240 rexpadj = pow(2.0, (double)expadj);
241 } else {
242 if (!isflt(expcomp))
243 goto userr;
244 rexpadj = atof(expcomp);
245 expadj = log(rexpadj)/LOG2 + (rexpadj>1 ? .5 : -.5);
246 if (!(doavg | doblur))
247 rexpadj = pow(2.0, (double)expadj);
248 }
249 /* set view */
250 if (nextview(doblur ? stdin : (FILE *)NULL) == EOF) {
251 fprintf(stderr, "%s: no view on standard input!\n",
252 progname);
253 exit(1);
254 }
255 normaspect(viewaspect(&ourview), &pixaspect, &hresolu, &vresolu);
256 /* allocate frame */
257 if (doavg) {
258 ourspict = (COLOR *)bmalloc(hresolu*vresolu*sizeof(COLOR));
259 ourweigh = (float *)bmalloc(hresolu*vresolu*sizeof(float));
260 if (ourspict == NULL | ourweigh == NULL)
261 syserror(progname);
262 } else {
263 ourpict = (COLR *)bmalloc(hresolu*vresolu*sizeof(COLR));
264 if (ourpict == NULL)
265 syserror(progname);
266 }
267 if (doblur) {
268 ourbpict = (COLOR *)bmalloc(hresolu*vresolu*sizeof(COLOR));
269 if (ourbpict == NULL)
270 syserror(progname);
271 }
272 ourzbuf = (float *)bmalloc(hresolu*vresolu*sizeof(float));
273 if (ourzbuf == NULL)
274 syserror(progname);
275 /* new header */
276 newheader("RADIANCE", stdout);
277 fputnow(stdout);
278 /* run pictures */
279 do {
280 bzero((char *)ourzbuf, hresolu*vresolu*sizeof(float));
281 for (i = an; i < argc; i += 2)
282 addpicture(argv[i], argv[i+1]);
283 if (fillo&F_BACK) /* fill in spaces */
284 backpicture(fillfunc, fillsamp);
285 else
286 fillpicture(fillfunc);
287 /* aft clipping */
288 clipaft();
289 } while (addblur() && nextview(stdin) != EOF);
290 /* close calculation */
291 caldone();
292 /* add to header */
293 printargs(argc, argv, stdout);
294 compavgview();
295 if (doblur | gotvfile) {
296 fputs(VIEWSTR, stdout);
297 fprintview(&avgview, stdout);
298 putc('\n', stdout);
299 }
300 if (pixaspect < .99 | pixaspect > 1.01)
301 fputaspect(pixaspect, stdout);
302 if (ourexp > 0)
303 ourexp *= rexpadj;
304 else
305 ourexp = rexpadj;
306 if (ourexp < .995 | ourexp > 1.005)
307 fputexpos(ourexp, stdout);
308 if (strcmp(ourfmt, PICFMT)) /* print format if known */
309 fputformat(ourfmt, stdout);
310 putc('\n', stdout);
311 /* write picture */
312 writepicture();
313 /* write z file */
314 if (zfile != NULL)
315 writedistance(zfile);
316
317 exit(0);
318 userr:
319 fprintf(stderr,
320 "Usage: %s [view opts][-t eps][-z zout][-e spec][-B][-a|-q][-fT][-n] pfile zspec ..\n",
321 progname);
322 exit(1);
323 #undef check
324 }
325
326
327 int
328 headline(s) /* process header string */
329 char *s;
330 {
331 char fmt[32];
332
333 if (isheadid(s))
334 return(0);
335 if (formatval(fmt, s)) {
336 if (globmatch(ourfmt, fmt)) {
337 wrongformat = 0;
338 strcpy(ourfmt, fmt);
339 } else
340 wrongformat = 1;
341 return(0);
342 }
343 if (nvavg < 2) {
344 putc('\t', stdout);
345 fputs(s, stdout);
346 }
347 if (isexpos(s)) {
348 theirexp *= exposval(s);
349 return(0);
350 }
351 if (isview(s) && sscanview(&theirview, s) > 0)
352 gotview++;
353 return(0);
354 }
355
356
357 nextview(fp) /* get and set next view */
358 FILE *fp;
359 {
360 char linebuf[256];
361 char *err;
362 register int i;
363
364 if (fp != NULL) {
365 do /* get new view */
366 if (fgets(linebuf, sizeof(linebuf), fp) == NULL)
367 return(EOF);
368 while (!isview(linebuf) || !sscanview(&ourview, linebuf));
369 }
370 /* set new view */
371 if ((err = setview(&ourview)) != NULL) {
372 fprintf(stderr, "%s: %s\n", progname, err);
373 exit(1);
374 }
375 if (!nvavg) { /* first view */
376 copystruct(&avgview, &ourview);
377 return(nvavg++);
378 }
379 /* add to average view */
380 for (i = 0; i < 3; i++) {
381 avgview.vp[i] += ourview.vp[i];
382 avgview.vdir[i] += ourview.vdir[i];
383 avgview.vup[i] += ourview.vup[i];
384 }
385 avgview.horiz += ourview.horiz;
386 avgview.vert += ourview.vert;
387 avgview.hoff += ourview.hoff;
388 avgview.voff += ourview.voff;
389 avgview.vfore += ourview.vfore;
390 avgview.vaft += ourview.vaft;
391 return(nvavg++);
392 }
393
394
395 compavgview() /* compute average view */
396 {
397 register int i;
398 double f;
399
400 if (nvavg < 2)
401 return;
402 f = 1.0/nvavg;
403 for (i = 0; i < 3; i++) {
404 avgview.vp[i] *= f;
405 avgview.vdir[i] *= f;
406 avgview.vup[i] *= f;
407 }
408 avgview.horiz *= f;
409 avgview.vert *= f;
410 avgview.hoff *= f;
411 avgview.voff *= f;
412 avgview.vfore *= f;
413 avgview.vaft *= f;
414 if (setview(&avgview) != NULL) /* in case of emergency... */
415 copystruct(&avgview, &ourview);
416 pixaspect = viewaspect(&avgview) * hresolu / vresolu;
417 }
418
419
420 addpicture(pfile, zspec) /* add picture to output */
421 char *pfile, *zspec;
422 {
423 FILE *pfp;
424 int zfd;
425 char *err;
426 COLR *scanin;
427 float *zin;
428 struct position *plast;
429 struct bound *xlim, ylim;
430 int y;
431 /* open picture file */
432 if ((pfp = fopen(pfile, "r")) == NULL)
433 syserror(pfile);
434 /* get header with exposure and view */
435 theirexp = 1.0;
436 copystruct(&theirview, &stdview);
437 gotview = 0;
438 if (nvavg < 2)
439 printf("%s:\n", pfile);
440 getheader(pfp, headline, NULL);
441 if (wrongformat || !gotview || !fgetsresolu(&tresolu, pfp)) {
442 fprintf(stderr, "%s: picture format error\n", pfile);
443 exit(1);
444 }
445 if (ourexp <= 0)
446 ourexp = theirexp;
447 else if (ABS(theirexp-ourexp) > .01*ourexp)
448 fprintf(stderr, "%s: different exposure (warning)\n", pfile);
449 if (err = setview(&theirview)) {
450 fprintf(stderr, "%s: %s\n", pfile, err);
451 exit(1);
452 }
453 /* compute transformation */
454 hasmatrix = pixform(theirs2ours, &theirview, &ourview);
455 /* get z specification or file */
456 zin = (float *)malloc(scanlen(&tresolu)*sizeof(float));
457 if (zin == NULL)
458 syserror(progname);
459 if ((zfd = open(zspec, O_RDONLY)) == -1) {
460 double zvalue;
461 register int x;
462 if (!isflt(zspec) || (zvalue = atof(zspec)) <= 0.0)
463 syserror(zspec);
464 for (x = scanlen(&tresolu); x-- > 0; )
465 zin[x] = zvalue;
466 }
467 /* compute transferrable perimeter */
468 xlim = (struct bound *)malloc(numscans(&tresolu)*sizeof(struct bound));
469 if (xlim == NULL)
470 syserror(progname);
471 if (!getperim(xlim, &ylim, zin, zfd)) { /* overlapping area? */
472 free((void *)zin);
473 free((void *)xlim);
474 if (zfd != -1)
475 close(zfd);
476 fclose(pfp);
477 return;
478 }
479 /* allocate scanlines */
480 scanin = (COLR *)malloc(scanlen(&tresolu)*sizeof(COLR));
481 plast = (struct position *)calloc(scanlen(&tresolu),
482 sizeof(struct position));
483 if (scanin == NULL | plast == NULL)
484 syserror(progname);
485 /* skip to starting point */
486 for (y = 0; y < ylim.min; y++)
487 if (freadcolrs(scanin, scanlen(&tresolu), pfp) < 0) {
488 fprintf(stderr, "%s: read error\n", pfile);
489 exit(1);
490 }
491 if (zfd != -1 && lseek(zfd,
492 (off_t)ylim.min*scanlen(&tresolu)*sizeof(float), 0) < 0)
493 syserror(zspec);
494 /* load image */
495 for (y = ylim.min; y <= ylim.max; y++) {
496 if (freadcolrs(scanin, scanlen(&tresolu), pfp) < 0) {
497 fprintf(stderr, "%s: read error\n", pfile);
498 exit(1);
499 }
500 if (zfd != -1 && read(zfd, (char *)zin,
501 scanlen(&tresolu)*sizeof(float))
502 < scanlen(&tresolu)*sizeof(float))
503 syserror(zspec);
504 addscanline(xlim+y, y, scanin, zin, plast);
505 }
506 /* clean up */
507 free((void *)xlim);
508 free((void *)scanin);
509 free((void *)zin);
510 free((void *)plast);
511 fclose(pfp);
512 if (zfd != -1)
513 close(zfd);
514 }
515
516
517 pixform(xfmat, vw1, vw2) /* compute view1 to view2 matrix */
518 register MAT4 xfmat;
519 register VIEW *vw1, *vw2;
520 {
521 double m4t[4][4];
522
523 if (vw1->type != VT_PER & vw1->type != VT_PAR)
524 return(0);
525 if (vw2->type != VT_PER & vw2->type != VT_PAR)
526 return(0);
527 setident4(xfmat);
528 xfmat[0][0] = vw1->hvec[0];
529 xfmat[0][1] = vw1->hvec[1];
530 xfmat[0][2] = vw1->hvec[2];
531 xfmat[1][0] = vw1->vvec[0];
532 xfmat[1][1] = vw1->vvec[1];
533 xfmat[1][2] = vw1->vvec[2];
534 xfmat[2][0] = vw1->vdir[0];
535 xfmat[2][1] = vw1->vdir[1];
536 xfmat[2][2] = vw1->vdir[2];
537 xfmat[3][0] = vw1->vp[0];
538 xfmat[3][1] = vw1->vp[1];
539 xfmat[3][2] = vw1->vp[2];
540 setident4(m4t);
541 m4t[0][0] = vw2->hvec[0]/vw2->hn2;
542 m4t[1][0] = vw2->hvec[1]/vw2->hn2;
543 m4t[2][0] = vw2->hvec[2]/vw2->hn2;
544 m4t[3][0] = -DOT(vw2->vp,vw2->hvec)/vw2->hn2;
545 m4t[0][1] = vw2->vvec[0]/vw2->vn2;
546 m4t[1][1] = vw2->vvec[1]/vw2->vn2;
547 m4t[2][1] = vw2->vvec[2]/vw2->vn2;
548 m4t[3][1] = -DOT(vw2->vp,vw2->vvec)/vw2->vn2;
549 m4t[0][2] = vw2->vdir[0];
550 m4t[1][2] = vw2->vdir[1];
551 m4t[2][2] = vw2->vdir[2];
552 m4t[3][2] = -DOT(vw2->vp,vw2->vdir);
553 multmat4(xfmat, xfmat, m4t);
554 return(1);
555 }
556
557
558 addscanline(xl, y, pline, zline, lasty) /* add scanline to output */
559 struct bound *xl;
560 int y;
561 COLR *pline;
562 float *zline;
563 struct position *lasty; /* input/output */
564 {
565 FVECT pos;
566 struct position lastx, newpos;
567 double wt;
568 register int x;
569
570 lastx.z = 0;
571 for (x = xl->max; x >= xl->min; x--) {
572 pix2loc(pos, &tresolu, x, y);
573 pos[2] = zline[x];
574 if ((wt = movepixel(pos)) <= FTINY) {
575 lasty[x].z = lastx.z = 0; /* mark invalid */
576 continue;
577 }
578 /* add pixel to our image */
579 newpos.x = pos[0] * hresolu;
580 newpos.y = pos[1] * vresolu;
581 newpos.z = zline[x];
582 addpixel(&newpos, &lastx, &lasty[x], pline[x], wt, pos[2]);
583 lasty[x].x = lastx.x = newpos.x;
584 lasty[x].y = lastx.y = newpos.y;
585 lasty[x].z = lastx.z = newpos.z;
586 }
587 }
588
589
590 addpixel(p0, p1, p2, pix, w, z) /* fill in pixel parallelogram */
591 struct position *p0, *p1, *p2;
592 COLR pix;
593 double w;
594 double z;
595 {
596 double zt = 2.*zeps*p0->z; /* threshold */
597 COLOR pval; /* converted+weighted pixel */
598 int s1x, s1y, s2x, s2y; /* step sizes */
599 int l1, l2, c1, c2; /* side lengths and counters */
600 int p1isy; /* p0p1 along y? */
601 int x1, y1; /* p1 position */
602 register int x, y; /* final position */
603
604 /* compute vector p0p1 */
605 if (fillo&F_FORE && ABS(p1->z-p0->z) <= zt) {
606 s1x = p1->x - p0->x;
607 s1y = p1->y - p0->y;
608 l1 = ABS(s1x);
609 if (p1isy = (ABS(s1y) > l1))
610 l1 = ABS(s1y);
611 else if (l1 < 1)
612 l1 = 1;
613 } else {
614 l1 = s1x = s1y = 1;
615 p1isy = -1;
616 }
617 /* compute vector p0p2 */
618 if (fillo&F_FORE && ABS(p2->z-p0->z) <= zt) {
619 s2x = p2->x - p0->x;
620 s2y = p2->y - p0->y;
621 if (p1isy == 1)
622 l2 = ABS(s2x);
623 else {
624 l2 = ABS(s2y);
625 if (p1isy != 0 && ABS(s2x) > l2)
626 l2 = ABS(s2x);
627 }
628 if (l2 < 1)
629 l2 = 1;
630 } else
631 l2 = s2x = s2y = 1;
632 /* fill the parallelogram */
633 if (averaging) {
634 colr_color(pval, pix);
635 scalecolor(pval, w);
636 }
637 for (c1 = l1; c1-- > 0; ) {
638 x1 = p0->x + c1*s1x/l1;
639 y1 = p0->y + c1*s1y/l1;
640 for (c2 = l2; c2-- > 0; ) {
641 x = x1 + c2*s2x/l2;
642 if (x < 0 | x >= hresolu)
643 continue;
644 y = y1 + c2*s2y/l2;
645 if (y < 0 | y >= vresolu)
646 continue;
647 if (averaging) {
648 if (zscan(y)[x] <= 0 || zscan(y)[x]-z
649 > zeps*zscan(y)[x]) {
650 copycolor(sscan(y)[x], pval);
651 wscan(y)[x] = w;
652 zscan(y)[x] = z;
653 } else if (z-zscan(y)[x] <= zeps*zscan(y)[x]) {
654 addcolor(sscan(y)[x], pval);
655 wscan(y)[x] += w;
656 }
657 } else if (zscan(y)[x] <= 0 || zscan(y)[x]-z
658 > zeps*zscan(y)[x]) {
659 copycolr(pscan(y)[x], pix);
660 zscan(y)[x] = z;
661 }
662 }
663 }
664 }
665
666
667 double
668 movepixel(pos) /* reposition image point */
669 register FVECT pos;
670 {
671 FVECT pt, tdir, odir;
672 double d;
673
674 if (pos[2] <= 0) /* empty pixel */
675 return(0);
676 if (usematrix) {
677 pos[0] += theirview.hoff - .5;
678 pos[1] += theirview.voff - .5;
679 if (normdist & theirview.type == VT_PER)
680 d = sqrt(1. + pos[0]*pos[0]*theirview.hn2
681 + pos[1]*pos[1]*theirview.vn2);
682 else
683 d = 1.;
684 pos[2] += d*theirview.vfore;
685 if (theirview.type == VT_PER) {
686 pos[2] /= d;
687 pos[0] *= pos[2];
688 pos[1] *= pos[2];
689 }
690 multp3(pos, pos, theirs2ours);
691 if (pos[2] <= ourview.vfore)
692 return(0);
693 if (ourview.type == VT_PER) {
694 pos[0] /= pos[2];
695 pos[1] /= pos[2];
696 }
697 pos[0] += .5 - ourview.hoff;
698 pos[1] += .5 - ourview.voff;
699 pos[2] -= ourview.vfore;
700 } else {
701 if (viewray(pt, tdir, &theirview, pos[0], pos[1]) < -FTINY)
702 return(0);
703 if (!normdist & theirview.type == VT_PER) /* adjust */
704 pos[2] *= sqrt(1. + pos[0]*pos[0]*theirview.hn2
705 + pos[1]*pos[1]*theirview.vn2);
706 pt[0] += tdir[0]*pos[2];
707 pt[1] += tdir[1]*pos[2];
708 pt[2] += tdir[2]*pos[2];
709 viewloc(pos, &ourview, pt);
710 if (pos[2] <= 0)
711 return(0);
712 }
713 if (pos[0] < 0 | pos[0] >= 1-FTINY | pos[1] < 0 | pos[1] >= 1-FTINY)
714 return(0);
715 if (!averaging)
716 return(1);
717 /* compute pixel weight */
718 if (ourview.type == VT_PAR) {
719 d = DOT(ourview.vdir,tdir);
720 d = 1. - d*d;
721 } else {
722 VSUB(odir, pt, ourview.vp);
723 d = DOT(odir,tdir);
724 d = 1. - d*d/DOT(odir,odir);
725 }
726 if (d <= 1./MAXWT/MAXWT)
727 return(MAXWT); /* clip to maximum weight */
728 return(1./sqrt(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, (off_t)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((void *)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 (PDesc.running) {
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 (!PDesc.running)
1133 return;
1134 clearqueue();
1135 close_process(&PDesc);
1136 }
1137
1138
1139 rcalfill(x, y) /* fill with ray-calculated pixel */
1140 int x, y;
1141 {
1142 if (queuesiz >= packsiz) /* flush queue if needed */
1143 clearqueue();
1144 /* add position to queue */
1145 queue[queuesiz][0] = x;
1146 queue[queuesiz][1] = y;
1147 queuesiz++;
1148 }
1149
1150
1151 clearqueue() /* process queue */
1152 {
1153 FVECT orig, dir;
1154 float fbuf[6*(PACKSIZ+1)];
1155 register float *fbp;
1156 register int i;
1157 double vx, vy;
1158
1159 if (queuesiz == 0)
1160 return;
1161 fbp = fbuf;
1162 for (i = 0; i < queuesiz; i++) {
1163 viewray(orig, dir, &ourview,
1164 (queue[i][0]+.5)/hresolu,
1165 (queue[i][1]+.5)/vresolu);
1166 *fbp++ = orig[0]; *fbp++ = orig[1]; *fbp++ = orig[2];
1167 *fbp++ = dir[0]; *fbp++ = dir[1]; *fbp++ = dir[2];
1168 }
1169 /* mark end and get results */
1170 bzero((char *)fbp, 6*sizeof(float));
1171 if (process(&PDesc, (char *)fbuf, (char *)fbuf,
1172 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 }