ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pinterp.c
Revision: 2.19
Committed: Fri Dec 23 22:35:12 1994 UTC (29 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.18: +143 -62 lines
Log Message: 
added pixel averaging (-a and -q options)

File Contents

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