ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.21
Committed: Mon Dec 26 21:02:31 1994 UTC (29 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.20: +42 -10 lines
Log Message:
add -e exposure option and fixed small bug

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