ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.22
Committed: Tue Dec 27 18:04:29 1994 UTC (29 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.21: +28 -21 lines
Log Message:
bug fixes related to -n, -z and -fr options

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 (theirview.type == VT_PER) {
593 if (normdist) /* adjust distance */
594 pos[2] /= sqrt(1. + pos[0]*pos[0]*theirview.hn2
595 + pos[1]*pos[1]*theirview.vn2);
596 pos[0] *= pos[2];
597 pos[1] *= pos[2];
598 }
599 multp3(pos, pos, theirs2ours);
600 if (pos[2] <= 0)
601 return(0);
602 if (ourview.type == VT_PER) {
603 pos[0] /= pos[2];
604 pos[1] /= pos[2];
605 }
606 pos[0] += .5 - ourview.hoff;
607 pos[1] += .5 - ourview.voff;
608 } else {
609 if (viewray(pt, tdir, &theirview, pos[0], pos[1]) < -FTINY)
610 return(0);
611 if (!normdist & theirview.type == VT_PER) /* adjust */
612 pos[2] *= sqrt(1. + pos[0]*pos[0]*theirview.hn2
613 + pos[1]*pos[1]*theirview.vn2);
614 pt[0] += tdir[0]*pos[2];
615 pt[1] += tdir[1]*pos[2];
616 pt[2] += tdir[2]*pos[2];
617 viewloc(pos, &ourview, pt);
618 if (pos[2] <= 0)
619 return(0);
620 }
621 if (pos[0] < 0 || pos[0] >= 1-FTINY || pos[1] < 0 || pos[1] >= 1-FTINY)
622 return(0);
623 if (!averaging)
624 return(1);
625 if (ourview.type == VT_PAR) /* compute our direction */
626 VCOPY(odir, ourview.vdir);
627 else
628 for (i = 0; i < 3; i++)
629 odir[i] = (pt[i] - ourview.vp[i])/pos[2];
630 d = DOT(odir,tdir); /* compute pixel weight */
631 if (d >= 1.-1./MAXWT/MAXWT)
632 return(MAXWT); /* clip to maximum weight */
633 return(1./sqrt(1.-d));
634 }
635
636
637 getperim(xl, yl, zline, zfd) /* compute overlapping image area */
638 register struct bound *xl;
639 struct bound *yl;
640 float *zline;
641 int zfd;
642 {
643 int step;
644 FVECT pos;
645 register int x, y;
646 /* set up step size */
647 if (scanlen(&tresolu) < numscans(&tresolu))
648 step = scanlen(&tresolu)/NSTEPS;
649 else
650 step = numscans(&tresolu)/NSTEPS;
651 if (step < MINSTEP) { /* not worth cropping? */
652 yl->min = 0;
653 yl->max = numscans(&tresolu) - 1;
654 x = scanlen(&tresolu) - 1;
655 for (y = numscans(&tresolu); y--; ) {
656 xl[y].min = 0;
657 xl[y].max = x;
658 }
659 return(1);
660 }
661 yl->min = 32000; yl->max = 0; /* search for points on image */
662 for (y = step - 1; y < numscans(&tresolu); y += step) {
663 if (zfd != -1) {
664 if (lseek(zfd, (long)y*scanlen(&tresolu)*sizeof(float),
665 0) < 0)
666 syserror("lseek");
667 if (read(zfd, (char *)zline,
668 scanlen(&tresolu)*sizeof(float))
669 < scanlen(&tresolu)*sizeof(float))
670 syserror("read");
671 }
672 xl[y].min = 32000; xl[y].max = 0; /* x max */
673 for (x = scanlen(&tresolu); (x -= step) > 0; ) {
674 pix2loc(pos, &tresolu, x, y);
675 pos[2] = zline[x];
676 if (movepixel(pos) > FTINY) {
677 xl[y].max = x + step - 1;
678 xl[y].min = x - step + 1; /* x min */
679 if (xl[y].min < 0)
680 xl[y].min = 0;
681 for (x = step - 1; x < xl[y].max; x += step) {
682 pix2loc(pos, &tresolu, x, y);
683 pos[2] = zline[x];
684 if (movepixel(pos) > FTINY) {
685 xl[y].min = x - step + 1;
686 break;
687 }
688 }
689 if (y < yl->min) /* y limits */
690 yl->min = y - step + 1;
691 yl->max = y + step - 1;
692 break;
693 }
694 }
695 /* fill in between */
696 if (y < step) {
697 xl[y-1].min = xl[y].min;
698 xl[y-1].max = xl[y].max;
699 } else {
700 if (xl[y].min < xl[y-step].min)
701 xl[y-1].min = xl[y].min;
702 else
703 xl[y-1].min = xl[y-step].min;
704 if (xl[y].max > xl[y-step].max)
705 xl[y-1].max = xl[y].max;
706 else
707 xl[y-1].max = xl[y-step].max;
708 }
709 for (x = 2; x < step; x++)
710 copystruct(xl+y-x, xl+y-1);
711 }
712 if (yl->max >= numscans(&tresolu))
713 yl->max = numscans(&tresolu) - 1;
714 y -= step;
715 for (x = numscans(&tresolu) - 1; x > y; x--) /* fill bottom rows */
716 copystruct(xl+x, xl+y);
717 return(yl->max >= yl->min);
718 }
719
720
721 backpicture(fill, samp) /* background fill algorithm */
722 int (*fill)();
723 int samp;
724 {
725 int *yback, xback;
726 int y;
727 register int x, i;
728 /* get back buffer */
729 yback = (int *)malloc(hresolu*sizeof(int));
730 if (yback == NULL)
731 syserror(progname);
732 for (x = 0; x < hresolu; x++)
733 yback[x] = -2;
734 /*
735 * Xback and yback are the pixel locations of suitable
736 * background values in each direction.
737 * A value of -2 means unassigned, and -1 means
738 * that there is no suitable background in this direction.
739 */
740 /* fill image */
741 for (y = 0; y < vresolu; y++) {
742 xback = -2;
743 for (x = 0; x < hresolu; x++)
744 if (zscan(y)[x] <= 0) { /* empty pixel */
745 /*
746 * First, find background from above or below.
747 * (farthest assigned pixel)
748 */
749 if (yback[x] == -2) {
750 for (i = y+1; i < vresolu; i++)
751 if (zscan(i)[x] > 0)
752 break;
753 if (i < vresolu
754 && (y <= 0 || zscan(y-1)[x] < zscan(i)[x]))
755 yback[x] = i;
756 else
757 yback[x] = y-1;
758 }
759 /*
760 * Next, find background from left or right.
761 */
762 if (xback == -2) {
763 for (i = x+1; i < hresolu; i++)
764 if (zscan(y)[i] > 0)
765 break;
766 if (i < hresolu
767 && (x <= 0 || zscan(y)[x-1] < zscan(y)[i]))
768 xback = i;
769 else
770 xback = x-1;
771 }
772 /*
773 * If we have no background for this pixel,
774 * use the given fill function.
775 */
776 if (xback < 0 && yback[x] < 0)
777 goto fillit;
778 /*
779 * Compare, and use the background that is
780 * farther, unless one of them is next to us.
781 * If the background is too distant, call
782 * the fill function.
783 */
784 if ( yback[x] < 0
785 || (xback >= 0 && ABS(x-xback) <= 1)
786 || ( ABS(y-yback[x]) > 1
787 && zscan(yback[x])[x]
788 < zscan(y)[xback] ) ) {
789 if (samp > 0 && ABS(x-xback) >= samp)
790 goto fillit;
791 if (averaging) {
792 copycolor(sscan(y)[x],
793 sscan(y)[xback]);
794 wscan(y)[x] = wscan(y)[xback];
795 } else
796 copycolr(pscan(y)[x],
797 pscan(y)[xback]);
798 zscan(y)[x] = zscan(y)[xback];
799 } else {
800 if (samp > 0 && ABS(y-yback[x]) > samp)
801 goto fillit;
802 if (averaging) {
803 copycolor(sscan(y)[x],
804 sscan(yback[x])[x]);
805 wscan(y)[x] =
806 wscan(yback[x])[x];
807 } else
808 copycolr(pscan(y)[x],
809 pscan(yback[x])[x]);
810 zscan(y)[x] = zscan(yback[x])[x];
811 }
812 continue;
813 fillit:
814 (*fill)(x,y);
815 if (fill == rcalfill) { /* use it */
816 clearqueue();
817 xback = x;
818 yback[x] = y;
819 }
820 } else { /* full pixel */
821 yback[x] = -2;
822 xback = -2;
823 }
824 }
825 free((char *)yback);
826 }
827
828
829 fillpicture(fill) /* paint in empty pixels using fill */
830 int (*fill)();
831 {
832 register int x, y;
833
834 for (y = 0; y < vresolu; y++)
835 for (x = 0; x < hresolu; x++)
836 if (zscan(y)[x] <= 0)
837 (*fill)(x,y);
838 }
839
840
841 clipaft() /* perform aft clipping as indicated */
842 {
843 register int x, y;
844 int adjtest = ourview.type == VT_PER & zisnorm;
845 double tstdist;
846 double yzn2, vx;
847
848 if (ourview.vaft <= FTINY)
849 return;
850 tstdist = ourview.vaft - ourview.vfore;
851 for (y = 0; y < vresolu; y++) {
852 if (adjtest) { /* adjust test */
853 yzn2 = (y+.5)/vresolu + ourview.voff - .5;
854 yzn2 = 1. + yzn2*yzn2*ourview.vn2;
855 tstdist = (ourview.vaft - ourview.vfore)*sqrt(yzn2);
856 }
857 for (x = 0; x < hresolu; x++)
858 if (zscan(y)[x] > tstdist) {
859 if (adjtest) {
860 vx = (x+.5)/hresolu + ourview.hoff - .5;
861 if (zscan(y)[x] <= (ourview.vaft -
862 ourview.vfore) *
863 sqrt(vx*vx*ourview.hn2 + yzn2))
864 continue;
865 }
866 if (averaging)
867 bzero(sscan(y)[x], sizeof(COLOR));
868 else
869 bzero(pscan(y)[x], sizeof(COLR));
870 zscan(y)[x] = 0.0;
871 }
872 }
873 }
874
875
876 writepicture() /* write out picture (alters buffer) */
877 {
878 int y;
879 register int x;
880 double d;
881
882 fprtresolu(hresolu, vresolu, stdout);
883 for (y = vresolu-1; y >= 0; y--)
884 if (averaging) {
885 for (x = 0; x < hresolu; x++) { /* average pixels */
886 d = rexpadj/wscan(y)[x];
887 scalecolor(sscan(y)[x], d);
888 }
889 if (fwritescan(sscan(y), hresolu, stdout) < 0)
890 syserror(progname);
891 } else {
892 if (expadj)
893 shiftcolrs(pscan(y), hresolu, expadj);
894 if (fwritecolrs(pscan(y), hresolu, stdout) < 0)
895 syserror(progname);
896 }
897 }
898
899
900 writedistance(fname) /* write out z file (alters buffer) */
901 char *fname;
902 {
903 int donorm = normdist & !zisnorm ? 1 :
904 ourview.type == VT_PER & !normdist & zisnorm ? -1 : 0;
905 int fd;
906 int y;
907
908 if ((fd = open(fname, O_WRONLY|O_CREAT|O_TRUNC, 0666)) == -1)
909 syserror(fname);
910 for (y = vresolu-1; y >= 0; y--) {
911 if (donorm) {
912 double vx, yzn2, d;
913 register int x;
914 yzn2 = (y+.5)/vresolu + ourview.voff - .5;
915 yzn2 = 1. + yzn2*yzn2*ourview.vn2;
916 for (x = 0; x < hresolu; x++) {
917 vx = (x+.5)/hresolu + ourview.hoff - .5;
918 d = sqrt(vx*vx*ourview.hn2 + yzn2);
919 if (donorm > 0)
920 zscan(y)[x] *= d;
921 else
922 zscan(y)[x] /= d;
923 }
924 }
925 if (write(fd, (char *)zscan(y), hresolu*sizeof(float))
926 < hresolu*sizeof(float))
927 syserror(fname);
928 }
929 close(fd);
930 }
931
932
933 backfill(x, y) /* fill pixel with background */
934 int x, y;
935 {
936 if (averaging) {
937 copycolor(sscan(y)[x], backcolor);
938 wscan(y)[x] = 1;
939 } else
940 copycolr(pscan(y)[x], backcolr);
941 zscan(y)[x] = backz;
942 }
943
944
945 calstart(prog, args) /* start fill calculation */
946 char *prog, *args;
947 {
948 char combuf[512];
949 char *argv[64];
950 int rval;
951 register char **wp, *cp;
952
953 if (childpid != -1) {
954 fprintf(stderr, "%s: too many calculations\n", progname);
955 exit(1);
956 }
957 strcpy(combuf, prog);
958 strcat(combuf, args);
959 cp = combuf;
960 wp = argv;
961 for ( ; ; ) {
962 while (isspace(*cp)) /* nullify spaces */
963 *cp++ = '\0';
964 if (!*cp) /* all done? */
965 break;
966 *wp++ = cp; /* add argument to list */
967 while (*++cp && !isspace(*cp))
968 ;
969 }
970 *wp = NULL;
971 /* start process */
972 if ((rval = open_process(PDesc, argv)) < 0)
973 syserror(progname);
974 if (rval == 0) {
975 fprintf(stderr, "%s: command not found\n", argv[0]);
976 exit(1);
977 }
978 packsiz = rval/(6*sizeof(float)) - 1;
979 if (packsiz > PACKSIZ)
980 packsiz = PACKSIZ;
981 queuesiz = 0;
982 }
983
984
985 caldone() /* done with calculation */
986 {
987 if (childpid == -1)
988 return;
989 clearqueue();
990 close_process(PDesc);
991 childpid = -1;
992 }
993
994
995 rcalfill(x, y) /* fill with ray-calculated pixel */
996 int x, y;
997 {
998 if (queuesiz >= packsiz) /* flush queue if needed */
999 clearqueue();
1000 /* add position to queue */
1001 queue[queuesiz][0] = x;
1002 queue[queuesiz][1] = y;
1003 queuesiz++;
1004 }
1005
1006
1007 clearqueue() /* process queue */
1008 {
1009 FVECT orig, dir;
1010 float fbuf[6*(PACKSIZ+1)];
1011 register float *fbp;
1012 register int i;
1013 double vx, vy;
1014
1015 if (queuesiz == 0)
1016 return;
1017 fbp = fbuf;
1018 for (i = 0; i < queuesiz; i++) {
1019 viewray(orig, dir, &ourview,
1020 (queue[i][0]+.5)/hresolu,
1021 (queue[i][1]+.5)/vresolu);
1022 *fbp++ = orig[0]; *fbp++ = orig[1]; *fbp++ = orig[2];
1023 *fbp++ = dir[0]; *fbp++ = dir[1]; *fbp++ = dir[2];
1024 }
1025 /* mark end and get results */
1026 bzero((char *)fbp, 6*sizeof(float));
1027 if (process(PDesc, fbuf, fbuf, 4*sizeof(float)*queuesiz,
1028 6*sizeof(float)*(queuesiz+1)) !=
1029 4*sizeof(float)*queuesiz) {
1030 fprintf(stderr, "%s: error reading from rtrace process\n",
1031 progname);
1032 exit(1);
1033 }
1034 fbp = fbuf;
1035 for (i = 0; i < queuesiz; i++) {
1036 if (ourexp > 0 && ourexp != 1.0) {
1037 fbp[0] *= ourexp;
1038 fbp[1] *= ourexp;
1039 fbp[2] *= ourexp;
1040 }
1041 if (averaging) {
1042 setcolor(sscan(queue[i][1])[queue[i][0]],
1043 fbp[0], fbp[1], fbp[2]);
1044 wscan(queue[i][1])[queue[i][0]] = 1;
1045 } else
1046 setcolr(pscan(queue[i][1])[queue[i][0]],
1047 fbp[0], fbp[1], fbp[2]);
1048 if (zisnorm)
1049 zscan(queue[i][1])[queue[i][0]] = fbp[3];
1050 else {
1051 vx = (queue[i][0]+.5)/hresolu + ourview.hoff - .5;
1052 vy = (queue[i][1]+.5)/vresolu + ourview.voff - .5;
1053 zscan(queue[i][1])[queue[i][0]] = fbp[3] / sqrt(1. +
1054 vx*vx*ourview.hn2 + vy*vy*ourview.vn2);
1055 }
1056 fbp += 4;
1057 }
1058 queuesiz = 0;
1059 }
1060
1061
1062 syserror(s) /* report error and exit */
1063 char *s;
1064 {
1065 perror(s);
1066 exit(1);
1067 }