ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.20
Committed: Mon Dec 26 19:57:22 1994 UTC (29 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.19: +19 -18 lines
Log Message:
bug fixes and improvements to last change (averaging)

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