ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.23
Committed: Fri Jan 6 13:19:37 1995 UTC (29 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.22: +9 -4 lines
Log Message:
bug fix in view fore clipping plane with matrix code

File Contents

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