ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.40
Committed: Fri Jan 2 12:47:01 2004 UTC (20 years, 4 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.39: +8 -4 lines
Log Message:
Fixed typing/prototype of getheader() and its callback.

File Contents

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