--- ray/src/px/pinterp.c 1990/01/18 23:58:29 1.23 +++ ray/src/px/pinterp.c 1998/10/27 09:08:26 2.31 @@ -1,3 +1,5 @@ +/* Copyright (c) 1995 Regents of the University of California */ + #ifndef lint static char SCCSid[] = "$SunId$ LBL"; #endif @@ -10,30 +12,44 @@ static char SCCSid[] = "$SunId$ LBL"; #include "standard.h" -#include +#include #include "view.h" #include "color.h" -#ifndef BSD -#define vfork fork -#endif +#include "resolu.h" +#define LOG2 0.69314718055994530942 + #define pscan(y) (ourpict+(y)*hresolu) +#define sscan(y) (ourspict+(y)*hresolu) +#define wscan(y) (ourweigh+(y)*hresolu) #define zscan(y) (ourzbuf+(y)*hresolu) +#define bscan(y) (ourbpict+(y)*hresolu) +#define averaging (ourweigh != NULL) +#define blurring (ourbpict != NULL) +#define usematrix (hasmatrix & !averaging) +#define zisnorm (!usematrix | ourview.type != VT_PER) +#define MAXWT 1000. /* maximum pixel weight (averaging) */ + #define F_FORE 1 /* fill foreground */ #define F_BACK 2 /* fill background */ -#define PACKSIZ 42 /* calculation packet size */ +#define PACKSIZ 256 /* max. calculation packet size */ -#define RTCOM "rtrace -h -ovl -fff -x %d %s" +#define RTCOM "rtrace -h- -ovl -fff -ld- -i- -I- " #define ABS(x) ((x)>0?(x):-(x)) struct position {int x,y; float z;}; +#define NSTEPS 64 /* number steps in overlap prescan */ +#define MINSTEP 4 /* minimum worthwhile preview step */ + +struct bound {int min,max;}; + VIEW ourview = STDVIEW; /* desired view */ int hresolu = 512; /* horizontal resolution */ int vresolu = 512; /* vertical resolution */ @@ -41,171 +57,261 @@ double pixaspect = 1.0; /* pixel aspect ratio */ double zeps = .02; /* allowed z epsilon */ -COLR *ourpict; /* output picture */ +COLR *ourpict; /* output picture (COLR's) */ +COLOR *ourspict; /* output pixel sums (averaging) */ +float *ourweigh = NULL; /* output pixel weights (averaging) */ float *ourzbuf; /* corresponding z-buffer */ +COLOR *ourbpict = NULL; /* blurred picture (view averaging) */ +VIEW avgview; /* average view for -B option */ +int nvavg; /* number of views averaged */ + char *progname; -int fill = F_FORE|F_BACK; /* selected fill algorithm */ +int fillo = F_FORE|F_BACK; /* selected fill options */ +int fillsamp = 0; /* sample separation (0 == inf) */ extern int backfill(), rcalfill(); /* fill functions */ -int (*deffill)() = backfill; /* selected fill function */ +int (*fillfunc)() = backfill; /* selected fill function */ COLR backcolr = BLKCOLR; /* background color */ +COLOR backcolor = BLKCOLOR; /* background color (float) */ double backz = 0.0; /* background z value */ -int normdist = 1; /* normalized distance? */ -double ourexp = -1; /* output picture exposure */ +int normdist = 1; /* i/o normalized distance? */ +char ourfmt[LPICFMT+1] = PICFMT; /* original picture format */ +double ourexp = -1; /* original picture exposure */ +int expadj = 0; /* exposure adjustment (f-stops) */ +double rexpadj = 1; /* real exposure adjustment */ -VIEW theirview = STDVIEW; /* input view */ +VIEW theirview; /* input view */ int gotview; /* got input view? */ -int thresolu, tvresolu; /* input resolution */ +int wrongformat = 0; /* input in another format? */ +RESOLU tresolu; /* input resolution */ double theirexp; /* input picture exposure */ -double theirs2ours[4][4]; /* transformation matrix */ +MAT4 theirs2ours; /* transformation matrix */ +int hasmatrix = 0; /* has transformation matrix */ -int childpid = -1; /* id of fill process */ -FILE *psend, *precv; /* pipes to/from fill calculation */ -int queue[PACKSIZ][2]; /* pending pixels */ -int queuesiz; /* number of pixels pending */ +int PDesc[3] = {-1,-1,-1}; /* rtrace process descriptor */ +#define childpid (PDesc[2]) +unsigned short queue[PACKSIZ][2]; /* pending pixels */ +int packsiz; /* actual packet size */ +int queuesiz = 0; /* number of pixels pending */ +extern double movepixel(); + main(argc, argv) /* interpolate pictures */ int argc; char *argv[]; { -#define check(olen,narg) if (argv[i][olen] || narg >= argc-i) goto badopt - extern double atof(); +#define check(ol,al) if (argv[an][ol] || \ + badarg(argc-an-1,argv+an+1,al)) \ + goto badopt int gotvfile = 0; + int doavg = -1; + int doblur = 0; char *zfile = NULL; - char *err; - int i, rval; + char *expcomp = NULL; + int i, an, rval; progname = argv[0]; - for (i = 1; i < argc && argv[i][0] == '-'; i++) { - rval = getviewopt(&ourview, argc-i, argv+i); + for (an = 1; an < argc && argv[an][0] == '-'; an++) { + rval = getviewopt(&ourview, argc-an, argv+an); if (rval >= 0) { - i += rval; + an += rval; continue; } - switch (argv[i][1]) { + switch (argv[an][1]) { + case 'e': /* exposure */ + check(2,"f"); + expcomp = argv[++an]; + break; case 't': /* threshold */ - check(2,1); - zeps = atof(argv[++i]); + check(2,"f"); + zeps = atof(argv[++an]); break; + case 'a': /* average */ + check(2,NULL); + doavg = 1; + break; + case 'B': /* blur views */ + check(2,NULL); + doblur = 1; + break; + case 'q': /* quick (no avg.) */ + check(2,NULL); + doavg = 0; + break; case 'n': /* dist. normalized? */ - check(2,0); + check(2,NULL); normdist = !normdist; break; case 'f': /* fill type */ - switch (argv[i][2]) { + switch (argv[an][2]) { case '0': /* none */ - check(3,0); - fill = 0; + check(3,NULL); + fillo = 0; break; case 'f': /* foreground */ - check(3,0); - fill = F_FORE; + check(3,NULL); + fillo = F_FORE; break; case 'b': /* background */ - check(3,0); - fill = F_BACK; + check(3,NULL); + fillo = F_BACK; break; case 'a': /* all */ - check(3,0); - fill = F_FORE|F_BACK; + check(3,NULL); + fillo = F_FORE|F_BACK; break; + case 's': /* sample */ + check(3,"i"); + fillsamp = atoi(argv[++an]); + break; case 'c': /* color */ - check(3,3); - deffill = backfill; - setcolr(backcolr, atof(argv[i+1]), - atof(argv[i+2]), atof(argv[i+3])); - i += 3; + check(3,"fff"); + fillfunc = backfill; + setcolor(backcolor, atof(argv[an+1]), + atof(argv[an+2]), atof(argv[an+3])); + setcolr(backcolr, colval(backcolor,RED), + colval(backcolor,GRN), + colval(backcolor,BLU)); + an += 3; break; case 'z': /* z value */ - check(3,1); - deffill = backfill; - backz = atof(argv[++i]); + check(3,"f"); + fillfunc = backfill; + backz = atof(argv[++an]); break; case 'r': /* rtrace */ - check(3,1); - deffill = rcalfill; - calstart(RTCOM, argv[++i]); + check(3,"s"); + fillfunc = rcalfill; + calstart(RTCOM, argv[++an]); break; default: goto badopt; } break; case 'z': /* z file */ - check(2,1); - zfile = argv[++i]; + check(2,"s"); + zfile = argv[++an]; break; case 'x': /* x resolution */ - check(2,1); - hresolu = atoi(argv[++i]); + check(2,"i"); + hresolu = atoi(argv[++an]); break; case 'y': /* y resolution */ - check(2,1); - vresolu = atoi(argv[++i]); + check(2,"i"); + vresolu = atoi(argv[++an]); break; case 'p': /* pixel aspect */ - check(2,1); - pixaspect = atof(argv[++i]); + if (argv[an][2] != 'a') + goto badopt; + check(3,"f"); + pixaspect = atof(argv[++an]); break; case 'v': /* view file */ - if (argv[i][2] != 'f') + if (argv[an][2] != 'f') goto badopt; - check(3,1); - gotvfile = viewfile(argv[++i], &ourview, 0, 0); - if (gotvfile < 0) { - perror(argv[i]); - exit(1); - } else if (gotvfile == 0) { + check(3,"s"); + gotvfile = viewfile(argv[++an], &ourview, 0, 0); + if (gotvfile < 0) + syserror(argv[an]); + else if (gotvfile == 0) { fprintf(stderr, "%s: bad view file\n", - argv[i]); + argv[an]); exit(1); } break; default: badopt: fprintf(stderr, "%s: command line error at '%s'\n", - progname, argv[i]); + progname, argv[an]); goto userr; } } /* check arguments */ - if ((argc-i)%2) + if ((argc-an)%2) goto userr; + if (fillsamp == 1) + fillo &= ~F_BACK; + if (doavg < 0) + doavg = (argc-an) > 2; + if (expcomp != NULL) + if (expcomp[0] == '+' | expcomp[0] == '-') { + expadj = atof(expcomp) + (expcomp[0]=='+' ? .5 : -.5); + if (doavg | doblur) + rexpadj = pow(2.0, atof(expcomp)); + else + rexpadj = pow(2.0, (double)expadj); + } else { + if (!isflt(expcomp)) + goto userr; + rexpadj = atof(expcomp); + expadj = log(rexpadj)/LOG2 + (rexpadj>1 ? .5 : -.5); + if (!(doavg | doblur)) + rexpadj = pow(2.0, (double)expadj); + } /* set view */ - if (err = setview(&ourview)) { - fprintf(stderr, "%s: %s\n", progname, err); + if (nextview(doblur ? stdin : (FILE *)NULL) == EOF) { + fprintf(stderr, "%s: no view on standard input!\n", + progname); exit(1); } normaspect(viewaspect(&ourview), &pixaspect, &hresolu, &vresolu); /* allocate frame */ - ourpict = (COLR *)malloc(hresolu*vresolu*sizeof(COLR)); - ourzbuf = (float *)calloc(hresolu*vresolu,sizeof(float)); - if (ourpict == NULL || ourzbuf == NULL) - syserror(); - /* get input */ - for ( ; i < argc; i += 2) - addpicture(argv[i], argv[i+1]); - /* fill in spaces */ - if (fill&F_BACK) - backpicture(); - else - fillpicture(); + if (doavg) { + ourspict = (COLOR *)bmalloc(hresolu*vresolu*sizeof(COLOR)); + ourweigh = (float *)bmalloc(hresolu*vresolu*sizeof(float)); + if (ourspict == NULL | ourweigh == NULL) + syserror(progname); + } else { + ourpict = (COLR *)bmalloc(hresolu*vresolu*sizeof(COLR)); + if (ourpict == NULL) + syserror(progname); + } + if (doblur) { + ourbpict = (COLOR *)bmalloc(hresolu*vresolu*sizeof(COLOR)); + if (ourbpict == NULL) + syserror(progname); + } + ourzbuf = (float *)bmalloc(hresolu*vresolu*sizeof(float)); + if (ourzbuf == NULL) + syserror(progname); + /* new header */ + newheader("RADIANCE", stdout); + /* run pictures */ + do { + bzero((char *)ourzbuf, hresolu*vresolu*sizeof(float)); + for (i = an; i < argc; i += 2) + addpicture(argv[i], argv[i+1]); + if (fillo&F_BACK) /* fill in spaces */ + backpicture(fillfunc, fillsamp); + else + fillpicture(fillfunc); + /* aft clipping */ + clipaft(); + } while (addblur() && nextview(stdin) != EOF); /* close calculation */ caldone(); /* add to header */ printargs(argc, argv, stdout); - if (gotvfile) { - printf(VIEWSTR); - fprintview(&ourview, stdout); - printf("\n"); + compavgview(); + if (doblur | gotvfile) { + fputs(VIEWSTR, stdout); + fprintview(&avgview, stdout); + putc('\n', stdout); } - if (pixaspect < .99 || pixaspect > 1.01) + if (pixaspect < .99 | pixaspect > 1.01) fputaspect(pixaspect, stdout); - if (ourexp > 0 && (ourexp < .995 || ourexp > 1.005)) + if (ourexp > 0) + ourexp *= rexpadj; + else + ourexp = rexpadj; + if (ourexp < .995 | ourexp > 1.005) fputexpos(ourexp, stdout); - printf("\n"); + if (strcmp(ourfmt, PICFMT)) /* print format if known */ + fputformat(ourfmt, stdout); + putc('\n', stdout); /* write picture */ writepicture(); /* write z file */ @@ -215,58 +321,129 @@ char *argv[]; exit(0); userr: fprintf(stderr, - "Usage: %s [view opts][-t eps][-z zout][-fT][-n] pfile zspec ..\n", + "Usage: %s [view opts][-t eps][-z zout][-e spec][-B][-a|-q][-fT][-n] pfile zspec ..\n", progname); exit(1); #undef check } +int headline(s) /* process header string */ char *s; { - static char *altname[] = {"rview","rpict","pinterp",VIEWSTR,NULL}; - register char **an; + char fmt[32]; - printf("\t%s", s); - + if (isheadid(s)) + return(0); + if (formatval(fmt, s)) { + if (globmatch(ourfmt, fmt)) { + wrongformat = 0; + strcpy(ourfmt, fmt); + } else + wrongformat = 1; + return(0); + } + if (nvavg < 2) { + putc('\t', stdout); + fputs(s, stdout); + } if (isexpos(s)) { theirexp *= exposval(s); + return(0); + } + if (isview(s) && sscanview(&theirview, s) > 0) + gotview++; + return(0); +} + + +nextview(fp) /* get and set next view */ +FILE *fp; +{ + char linebuf[256]; + char *err; + register int i; + + if (fp != NULL) { + do /* get new view */ + if (fgets(linebuf, sizeof(linebuf), fp) == NULL) + return(EOF); + while (!isview(linebuf) || !sscanview(&ourview, linebuf)); + } + /* set new view */ + if ((err = setview(&ourview)) != NULL) { + fprintf(stderr, "%s: %s\n", progname, err); + exit(1); + } + if (!nvavg) { /* first view */ + copystruct(&avgview, &ourview); + return(nvavg++); + } + /* add to average view */ + for (i = 0; i < 3; i++) { + avgview.vp[i] += ourview.vp[i]; + avgview.vdir[i] += ourview.vdir[i]; + avgview.vup[i] += ourview.vup[i]; + } + avgview.horiz += ourview.horiz; + avgview.vert += ourview.vert; + avgview.hoff += ourview.hoff; + avgview.voff += ourview.voff; + avgview.vfore += ourview.vfore; + avgview.vaft += ourview.vaft; + return(nvavg++); +} + + +compavgview() /* compute average view */ +{ + register int i; + double f; + + if (nvavg < 2) return; + f = 1.0/nvavg; + for (i = 0; i < 3; i++) { + avgview.vp[i] *= f; + avgview.vdir[i] *= f; + avgview.vup[i] *= f; } - for (an = altname; *an != NULL; an++) - if (!strncmp(*an, s, strlen(*an))) { - if (sscanview(&theirview, s+strlen(*an)) > 0) - gotview++; - break; - } + avgview.horiz *= f; + avgview.vert *= f; + avgview.hoff *= f; + avgview.voff *= f; + avgview.vfore *= f; + avgview.vaft *= f; + if (setview(&avgview) != NULL) /* in case of emergency... */ + copystruct(&avgview, &ourview); + pixaspect = viewaspect(&avgview) * hresolu / vresolu; } addpicture(pfile, zspec) /* add picture to output */ char *pfile, *zspec; { - extern double atof(); FILE *pfp; int zfd; char *err; COLR *scanin; float *zin; struct position *plast; + struct bound *xlim, ylim; int y; /* open picture file */ - if ((pfp = fopen(pfile, "r")) == NULL) { - perror(pfile); - exit(1); - } + if ((pfp = fopen(pfile, "r")) == NULL) + syserror(pfile); /* get header with exposure and view */ theirexp = 1.0; + copystruct(&theirview, &stdview); gotview = 0; - printf("%s:\n", pfile); - getheader(pfp, headline); - if (!gotview || fgetresolu(&thresolu, &tvresolu, pfp) - != (YMAJOR|YDECR)) { - fprintf(stderr, "%s: picture view error\n", pfile); + if (nvavg < 2) + printf("%s:\n", pfile); + getheader(pfp, headline, NULL); + if (wrongformat || !gotview || !fgetsresolu(&tresolu, pfp)) { + fprintf(stderr, "%s: picture format error\n", pfile); exit(1); } if (ourexp <= 0) @@ -278,38 +455,60 @@ char *pfile, *zspec; exit(1); } /* compute transformation */ - pixform(theirs2ours, &theirview, &ourview); - /* allocate scanlines */ - scanin = (COLR *)malloc(thresolu*sizeof(COLR)); - zin = (float *)malloc(thresolu*sizeof(float)); - plast = (struct position *)calloc(thresolu, sizeof(struct position)); - if (scanin == NULL || zin == NULL || plast == NULL) - syserror(); + hasmatrix = pixform(theirs2ours, &theirview, &ourview); /* get z specification or file */ + zin = (float *)malloc(scanlen(&tresolu)*sizeof(float)); + if (zin == NULL) + syserror(progname); if ((zfd = open(zspec, O_RDONLY)) == -1) { double zvalue; register int x; - if (!isfloat(zspec) || (zvalue = atof(zspec)) <= 0.0) { - perror(zspec); - exit(1); - } - for (x = 0; x < thresolu; x++) + if (!isflt(zspec) || (zvalue = atof(zspec)) <= 0.0) + syserror(zspec); + for (x = scanlen(&tresolu); x-- > 0; ) zin[x] = zvalue; } - /* load image */ - for (y = tvresolu-1; y >= 0; y--) { - if (freadcolrs(scanin, thresolu, pfp) < 0) { + /* compute transferrable perimeter */ + xlim = (struct bound *)malloc(numscans(&tresolu)*sizeof(struct bound)); + if (xlim == NULL) + syserror(progname); + if (!getperim(xlim, &ylim, zin, zfd)) { /* overlapping area? */ + free((char *)zin); + free((char *)xlim); + if (zfd != -1) + close(zfd); + fclose(pfp); + return; + } + /* allocate scanlines */ + scanin = (COLR *)malloc(scanlen(&tresolu)*sizeof(COLR)); + plast = (struct position *)calloc(scanlen(&tresolu), + sizeof(struct position)); + if (scanin == NULL | plast == NULL) + syserror(progname); + /* skip to starting point */ + for (y = 0; y < ylim.min; y++) + if (freadcolrs(scanin, scanlen(&tresolu), pfp) < 0) { fprintf(stderr, "%s: read error\n", pfile); exit(1); } - if (zfd != -1 && read(zfd,(char *)zin,thresolu*sizeof(float)) - < thresolu*sizeof(float)) { - fprintf(stderr, "%s: read error\n", zspec); + if (zfd != -1 && lseek(zfd, + (long)ylim.min*scanlen(&tresolu)*sizeof(float), 0) < 0) + syserror(zspec); + /* load image */ + for (y = ylim.min; y <= ylim.max; y++) { + if (freadcolrs(scanin, scanlen(&tresolu), pfp) < 0) { + fprintf(stderr, "%s: read error\n", pfile); exit(1); } - addscanline(y, scanin, zin, plast); + if (zfd != -1 && read(zfd, (char *)zin, + scanlen(&tresolu)*sizeof(float)) + < scanlen(&tresolu)*sizeof(float)) + syserror(zspec); + addscanline(xlim+y, y, scanin, zin, plast); } /* clean up */ + free((char *)xlim); free((char *)scanin); free((char *)zin); free((char *)plast); @@ -320,11 +519,15 @@ char *pfile, *zspec; pixform(xfmat, vw1, vw2) /* compute view1 to view2 matrix */ -register double xfmat[4][4]; +register MAT4 xfmat; register VIEW *vw1, *vw2; { double m4t[4][4]; + if (vw1->type != VT_PER & vw1->type != VT_PAR) + return(0); + if (vw2->type != VT_PER & vw2->type != VT_PAR) + return(0); setident4(xfmat); xfmat[0][0] = vw1->hvec[0]; xfmat[0][1] = vw1->hvec[1]; @@ -352,66 +555,50 @@ register VIEW *vw1, *vw2; m4t[2][2] = vw2->vdir[2]; m4t[3][2] = -DOT(vw2->vp,vw2->vdir); multmat4(xfmat, xfmat, m4t); + return(1); } -addscanline(y, pline, zline, lasty) /* add scanline to output */ +addscanline(xl, y, pline, zline, lasty) /* add scanline to output */ +struct bound *xl; int y; COLR *pline; float *zline; struct position *lasty; /* input/output */ { - extern double sqrt(); - double pos[3]; + FVECT pos; struct position lastx, newpos; + double wt; register int x; lastx.z = 0; - for (x = thresolu-1; x >= 0; x--) { - pos[0] = (x+.5)/thresolu + theirview.hoff - .5; - pos[1] = (y+.5)/tvresolu + theirview.voff - .5; + for (x = xl->max; x >= xl->min; x--) { + pix2loc(pos, &tresolu, x, y); pos[2] = zline[x]; - if (theirview.type == VT_PER) { - if (normdist) /* adjust for eye-ray distance */ - pos[2] /= sqrt( 1. - + pos[0]*pos[0]*theirview.hn2 - + pos[1]*pos[1]*theirview.vn2 ); - pos[0] *= pos[2]; - pos[1] *= pos[2]; - } - multp3(pos, pos, theirs2ours); - if (pos[2] <= 0) { + if ((wt = movepixel(pos)) <= FTINY) { lasty[x].z = lastx.z = 0; /* mark invalid */ continue; } - if (ourview.type == VT_PER) { - pos[0] /= pos[2]; - pos[1] /= pos[2]; - } - pos[0] += .5 - ourview.hoff; - pos[1] += .5 - ourview.voff; + /* add pixel to our image */ newpos.x = pos[0] * hresolu; newpos.y = pos[1] * vresolu; newpos.z = zline[x]; - /* add pixel to our image */ - if (pos[0] >= 0 && newpos.x < hresolu - && pos[1] >= 0 && newpos.y < vresolu) { - addpixel(&newpos, &lastx, &lasty[x], pline[x], pos[2]); - lasty[x].x = lastx.x = newpos.x; - lasty[x].y = lastx.y = newpos.y; - lasty[x].z = lastx.z = newpos.z; - } else - lasty[x].z = lastx.z = 0; /* mark invalid */ + addpixel(&newpos, &lastx, &lasty[x], pline[x], wt, pos[2]); + lasty[x].x = lastx.x = newpos.x; + lasty[x].y = lastx.y = newpos.y; + lasty[x].z = lastx.z = newpos.z; } } -addpixel(p0, p1, p2, pix, z) /* fill in pixel parallelogram */ +addpixel(p0, p1, p2, pix, w, z) /* fill in pixel parallelogram */ struct position *p0, *p1, *p2; COLR pix; +double w; double z; { double zt = 2.*zeps*p0->z; /* threshold */ + COLOR pval; /* converted+weighted pixel */ int s1x, s1y, s2x, s2y; /* step sizes */ int l1, l2, c1, c2; /* side lengths and counters */ int p1isy; /* p0p1 along y? */ @@ -419,18 +606,20 @@ double z; register int x, y; /* final position */ /* compute vector p0p1 */ - if (fill&F_FORE && ABS(p1->z-p0->z) <= zt) { + if (fillo&F_FORE && ABS(p1->z-p0->z) <= zt) { s1x = p1->x - p0->x; s1y = p1->y - p0->y; l1 = ABS(s1x); if (p1isy = (ABS(s1y) > l1)) l1 = ABS(s1y); + else if (l1 < 1) + l1 = 1; } else { l1 = s1x = s1y = 1; p1isy = -1; } /* compute vector p0p2 */ - if (fill&F_FORE && ABS(p2->z-p0->z) <= zt) { + if (fillo&F_FORE && ABS(p2->z-p0->z) <= zt) { s2x = p2->x - p0->x; s2y = p2->y - p0->y; if (p1isy == 1) @@ -440,35 +629,203 @@ double z; if (p1isy != 0 && ABS(s2x) > l2) l2 = ABS(s2x); } + if (l2 < 1) + l2 = 1; } else l2 = s2x = s2y = 1; /* fill the parallelogram */ + if (averaging) { + colr_color(pval, pix); + scalecolor(pval, w); + } for (c1 = l1; c1-- > 0; ) { x1 = p0->x + c1*s1x/l1; y1 = p0->y + c1*s1y/l1; for (c2 = l2; c2-- > 0; ) { x = x1 + c2*s2x/l2; + if (x < 0 | x >= hresolu) + continue; y = y1 + c2*s2y/l2; - if (zscan(y)[x] <= 0 || zscan(y)[x]-z + if (y < 0 | y >= vresolu) + continue; + if (averaging) { + if (zscan(y)[x] <= 0 || zscan(y)[x]-z > zeps*zscan(y)[x]) { - zscan(y)[x] = z; + copycolor(sscan(y)[x], pval); + wscan(y)[x] = w; + zscan(y)[x] = z; + } else if (z-zscan(y)[x] <= zeps*zscan(y)[x]) { + addcolor(sscan(y)[x], pval); + wscan(y)[x] += w; + } + } else if (zscan(y)[x] <= 0 || zscan(y)[x]-z + > zeps*zscan(y)[x]) { copycolr(pscan(y)[x], pix); + zscan(y)[x] = z; } } } } -backpicture() /* background fill algorithm */ +double +movepixel(pos) /* reposition image point */ +register FVECT pos; { + FVECT pt, tdir, odir; + double d; + register int i; + + if (pos[2] <= 0) /* empty pixel */ + return(0); + if (usematrix) { + pos[0] += theirview.hoff - .5; + pos[1] += theirview.voff - .5; + if (normdist & theirview.type == VT_PER) + d = sqrt(1. + pos[0]*pos[0]*theirview.hn2 + + pos[1]*pos[1]*theirview.vn2); + else + d = 1.; + pos[2] += d*theirview.vfore; + if (theirview.type == VT_PER) { + pos[2] /= d; + pos[0] *= pos[2]; + pos[1] *= pos[2]; + } + multp3(pos, pos, theirs2ours); + if (pos[2] <= ourview.vfore) + return(0); + if (ourview.type == VT_PER) { + pos[0] /= pos[2]; + pos[1] /= pos[2]; + } + pos[0] += .5 - ourview.hoff; + pos[1] += .5 - ourview.voff; + pos[2] -= ourview.vfore; + } else { + if (viewray(pt, tdir, &theirview, pos[0], pos[1]) < -FTINY) + return(0); + if (!normdist & theirview.type == VT_PER) /* adjust */ + pos[2] *= sqrt(1. + pos[0]*pos[0]*theirview.hn2 + + pos[1]*pos[1]*theirview.vn2); + pt[0] += tdir[0]*pos[2]; + pt[1] += tdir[1]*pos[2]; + pt[2] += tdir[2]*pos[2]; + viewloc(pos, &ourview, pt); + if (pos[2] <= 0) + return(0); + } + if (pos[0] < 0 | pos[0] >= 1-FTINY | pos[1] < 0 | pos[1] >= 1-FTINY) + return(0); + if (!averaging) + return(1); + if (ourview.type == VT_PAR) /* compute our direction */ + VCOPY(odir, ourview.vdir); + else + for (i = 0; i < 3; i++) + odir[i] = (pt[i] - ourview.vp[i])/pos[2]; + d = DOT(odir,tdir); /* compute pixel weight */ + if (d >= 1.-1./MAXWT/MAXWT) + return(MAXWT); /* clip to maximum weight */ + return(1./sqrt(1.-d)); +} + + +getperim(xl, yl, zline, zfd) /* compute overlapping image area */ +register struct bound *xl; +struct bound *yl; +float *zline; +int zfd; +{ + int step; + FVECT pos; + register int x, y; + /* set up step size */ + if (scanlen(&tresolu) < numscans(&tresolu)) + step = scanlen(&tresolu)/NSTEPS; + else + step = numscans(&tresolu)/NSTEPS; + if (step < MINSTEP) { /* not worth cropping? */ + yl->min = 0; + yl->max = numscans(&tresolu) - 1; + x = scanlen(&tresolu) - 1; + for (y = numscans(&tresolu); y--; ) { + xl[y].min = 0; + xl[y].max = x; + } + return(1); + } + yl->min = 32000; yl->max = 0; /* search for points on image */ + for (y = step - 1; y < numscans(&tresolu); y += step) { + if (zfd != -1) { + if (lseek(zfd, (long)y*scanlen(&tresolu)*sizeof(float), + 0) < 0) + syserror("lseek"); + if (read(zfd, (char *)zline, + scanlen(&tresolu)*sizeof(float)) + < scanlen(&tresolu)*sizeof(float)) + syserror("read"); + } + xl[y].min = 32000; xl[y].max = 0; /* x max */ + for (x = scanlen(&tresolu); (x -= step) > 0; ) { + pix2loc(pos, &tresolu, x, y); + pos[2] = zline[x]; + if (movepixel(pos) > FTINY) { + xl[y].max = x + step - 1; + xl[y].min = x - step + 1; /* x min */ + if (xl[y].min < 0) + xl[y].min = 0; + for (x = step - 1; x < xl[y].max; x += step) { + pix2loc(pos, &tresolu, x, y); + pos[2] = zline[x]; + if (movepixel(pos) > FTINY) { + xl[y].min = x - step + 1; + break; + } + } + if (y < yl->min) /* y limits */ + yl->min = y - step + 1; + yl->max = y + step - 1; + break; + } + } + /* fill in between */ + if (y < step) { + xl[y-1].min = xl[y].min; + xl[y-1].max = xl[y].max; + } else { + if (xl[y].min < xl[y-step].min) + xl[y-1].min = xl[y].min; + else + xl[y-1].min = xl[y-step].min; + if (xl[y].max > xl[y-step].max) + xl[y-1].max = xl[y].max; + else + xl[y-1].max = xl[y-step].max; + } + for (x = 2; x < step; x++) + copystruct(xl+y-x, xl+y-1); + } + if (yl->max >= numscans(&tresolu)) + yl->max = numscans(&tresolu) - 1; + y -= step; + for (x = numscans(&tresolu) - 1; x > y; x--) /* fill bottom rows */ + copystruct(xl+x, xl+y); + return(yl->max >= yl->min); +} + + +backpicture(fill, samp) /* background fill algorithm */ +int (*fill)(); +int samp; +{ int *yback, xback; int y; - COLR pfill; register int x, i; /* get back buffer */ yback = (int *)malloc(hresolu*sizeof(int)); if (yback == NULL) - syserror(); + syserror(progname); for (x = 0; x < hresolu; x++) yback[x] = -2; /* @@ -510,28 +867,53 @@ backpicture() /* background fill algorithm */ xback = x-1; } /* - * Check to see if we have no background for - * this pixel. If not, use background color. + * If we have no background for this pixel, + * use the given fill function. */ - if (xback < 0 && yback[x] < 0) { - (*deffill)(x,y); - continue; - } + if (xback < 0 && yback[x] < 0) + goto fillit; /* * Compare, and use the background that is * farther, unless one of them is next to us. + * If the background is too distant, call + * the fill function. */ if ( yback[x] < 0 || (xback >= 0 && ABS(x-xback) <= 1) || ( ABS(y-yback[x]) > 1 && zscan(yback[x])[x] < zscan(y)[xback] ) ) { - copycolr(pscan(y)[x],pscan(y)[xback]); + if (samp > 0 && ABS(x-xback) >= samp) + goto fillit; + if (averaging) { + copycolor(sscan(y)[x], + sscan(y)[xback]); + wscan(y)[x] = wscan(y)[xback]; + } else + copycolr(pscan(y)[x], + pscan(y)[xback]); zscan(y)[x] = zscan(y)[xback]; } else { - copycolr(pscan(y)[x],pscan(yback[x])[x]); + if (samp > 0 && ABS(y-yback[x]) > samp) + goto fillit; + if (averaging) { + copycolor(sscan(y)[x], + sscan(yback[x])[x]); + wscan(y)[x] = + wscan(yback[x])[x]; + } else + copycolr(pscan(y)[x], + pscan(yback[x])[x]); zscan(y)[x] = zscan(yback[x])[x]; } + continue; + fillit: + (*fill)(x,y); + if (fill == rcalfill) { /* use it */ + clearqueue(); + xback = x; + yback[x] = y; + } } else { /* full pixel */ yback[x] = -2; xback = -2; @@ -541,139 +923,218 @@ backpicture() /* background fill algorithm */ } -fillpicture() /* paint in empty pixels with default */ +fillpicture(fill) /* paint in empty pixels using fill */ +int (*fill)(); { register int x, y; for (y = 0; y < vresolu; y++) for (x = 0; x < hresolu; x++) if (zscan(y)[x] <= 0) - (*deffill)(x,y); + (*fill)(x,y); + if (fill == rcalfill) + clearqueue(); } -writepicture() /* write out picture */ +clipaft() /* perform aft clipping as indicated */ { + register int x, y; + int adjtest = ourview.type == VT_PER & zisnorm; + double tstdist; + double yzn2, vx; + + if (ourview.vaft <= FTINY) + return(0); + tstdist = ourview.vaft - ourview.vfore; + for (y = 0; y < vresolu; y++) { + if (adjtest) { /* adjust test */ + yzn2 = (y+.5)/vresolu + ourview.voff - .5; + yzn2 = 1. + yzn2*yzn2*ourview.vn2; + tstdist = (ourview.vaft - ourview.vfore)*sqrt(yzn2); + } + for (x = 0; x < hresolu; x++) + if (zscan(y)[x] > tstdist) { + if (adjtest) { + vx = (x+.5)/hresolu + ourview.hoff - .5; + if (zscan(y)[x] <= (ourview.vaft - + ourview.vfore) * + sqrt(vx*vx*ourview.hn2 + yzn2)) + continue; + } + if (averaging) + bzero(sscan(y)[x], sizeof(COLOR)); + else + bzero(pscan(y)[x], sizeof(COLR)); + zscan(y)[x] = 0.0; + } + } + return(1); +} + + +addblur() /* add to blurred picture */ +{ + COLOR cval; + double d; + register int i; + + if (!blurring) + return(0); + i = hresolu*vresolu; + if (nvavg < 2) + if (averaging) + while (i--) { + copycolor(ourbpict[i], ourspict[i]); + d = 1.0/ourweigh[i]; + scalecolor(ourbpict[i], d); + } + else + while (i--) + colr_color(ourbpict[i], ourpict[i]); + else + if (averaging) + while (i--) { + copycolor(cval, ourspict[i]); + d = 1.0/ourweigh[i]; + scalecolor(cval, d); + addcolor(ourbpict[i], cval); + } + else + while (i--) { + colr_color(cval, ourpict[i]); + addcolor(ourbpict[i], cval); + } + /* print view */ + printf("VIEW%d:", nvavg); + fprintview(&ourview, stdout); + putchar('\n'); + return(1); +} + + +writepicture() /* write out picture (alters buffer) */ +{ int y; + register int x; + double d; - fputresolu(YMAJOR|YDECR, hresolu, vresolu, stdout); + fprtresolu(hresolu, vresolu, stdout); for (y = vresolu-1; y >= 0; y--) - if (fwritecolrs(pscan(y), hresolu, stdout) < 0) - syserror(); + if (blurring) { + for (x = 0; x < hresolu; x++) { /* compute avg. */ + d = rexpadj/nvavg; + scalecolor(bscan(y)[x], d); + } + if (fwritescan(bscan(y), hresolu, stdout) < 0) + syserror(progname); + } else if (averaging) { + for (x = 0; x < hresolu; x++) { /* average pixels */ + d = rexpadj/wscan(y)[x]; + scalecolor(sscan(y)[x], d); + } + if (fwritescan(sscan(y), hresolu, stdout) < 0) + syserror(progname); + } else { + if (expadj) + shiftcolrs(pscan(y), hresolu, expadj); + if (fwritecolrs(pscan(y), hresolu, stdout) < 0) + syserror(progname); + } } -writedistance(fname) /* write out z file */ +writedistance(fname) /* write out z file (alters buffer) */ char *fname; { - extern double sqrt(); - int donorm = normdist && ourview.type == VT_PER; + int donorm = normdist & !zisnorm ? 1 : + ourview.type == VT_PER & !normdist & zisnorm ? -1 : 0; int fd; int y; - float *zout; - if ((fd = open(fname, O_WRONLY|O_CREAT|O_TRUNC, 0666)) == -1) { - perror(fname); - exit(1); - } - if (donorm - && (zout = (float *)malloc(hresolu*sizeof(float))) == NULL) - syserror(); + if ((fd = open(fname, O_WRONLY|O_CREAT|O_TRUNC, 0666)) == -1) + syserror(fname); for (y = vresolu-1; y >= 0; y--) { if (donorm) { - double vx, yzn2; + double vx, yzn2, d; register int x; - yzn2 = y - .5*(vresolu-1); + yzn2 = (y+.5)/vresolu + ourview.voff - .5; yzn2 = 1. + yzn2*yzn2*ourview.vn2; for (x = 0; x < hresolu; x++) { - vx = x - .5*(hresolu-1); - zout[x] = zscan(y)[x] - * sqrt(vx*vx*ourview.hn2 + yzn2); + vx = (x+.5)/hresolu + ourview.hoff - .5; + d = sqrt(vx*vx*ourview.hn2 + yzn2); + if (donorm > 0) + zscan(y)[x] *= d; + else + zscan(y)[x] /= d; } - } else - zout = zscan(y); - if (write(fd, (char *)zout, hresolu*sizeof(float)) - < hresolu*sizeof(float)) { - perror(fname); - exit(1); } + if (write(fd, (char *)zscan(y), hresolu*sizeof(float)) + < hresolu*sizeof(float)) + syserror(fname); } - if (donorm) - free((char *)zout); close(fd); } -isfloat(s) /* see if string is floating number */ -register char *s; -{ - for ( ; *s; s++) - if ((*s < '0' || *s > '9') && *s != '.' && *s != '-' - && *s != 'e' && *s != 'E' && *s != '+') - return(0); - return(1); -} - - backfill(x, y) /* fill pixel with background */ int x, y; { - register BYTE *dest = pscan(y)[x]; - - copycolr(dest, backcolr); + if (averaging) { + copycolor(sscan(y)[x], backcolor); + wscan(y)[x] = 1; + } else + copycolr(pscan(y)[x], backcolr); zscan(y)[x] = backz; } -calstart(prog, args) /* start fill calculation */ +calstart(prog, args) /* start fill calculation */ char *prog, *args; { char combuf[512]; - int p0[2], p1[2]; + char *argv[64]; + int rval; + register char **wp, *cp; if (childpid != -1) { fprintf(stderr, "%s: too many calculations\n", progname); exit(1); } - sprintf(combuf, prog, PACKSIZ, args); - if (pipe(p0) < 0 || pipe(p1) < 0) - syserror(); - if ((childpid = vfork()) == 0) { /* fork calculation */ - close(p0[1]); - close(p1[0]); - if (p0[0] != 0) { - dup2(p0[0], 0); - close(p0[0]); - } - if (p1[1] != 1) { - dup2(p1[1], 1); - close(p1[1]); - } - execl("/bin/sh", "sh", "-c", combuf, 0); - perror("/bin/sh"); - _exit(127); + strcpy(combuf, prog); + strcat(combuf, args); + cp = combuf; + wp = argv; + for ( ; ; ) { + while (isspace(*cp)) /* nullify spaces */ + *cp++ = '\0'; + if (!*cp) /* all done? */ + break; + *wp++ = cp; /* add argument to list */ + while (*++cp && !isspace(*cp)) + ; } - if (childpid == -1) - syserror(); - close(p0[0]); - close(p1[1]); - if ((psend = fdopen(p0[1], "w")) == NULL) - syserror(); - if ((precv = fdopen(p1[0], "r")) == NULL) - syserror(); + *wp = NULL; + /* start process */ + if ((rval = open_process(PDesc, argv)) < 0) + syserror(progname); + if (rval == 0) { + fprintf(stderr, "%s: command not found\n", argv[0]); + exit(1); + } + packsiz = rval/(6*sizeof(float)) - 1; + if (packsiz > PACKSIZ) + packsiz = PACKSIZ; queuesiz = 0; } -caldone() /* done with calculation */ +caldone() /* done with calculation */ { - int pid; - if (childpid == -1) return; - clearqueue(1); - while ((pid = wait(0)) != -1 && pid != childpid) - ; + clearqueue(); + close_process(PDesc); childpid = -1; } @@ -681,8 +1142,8 @@ caldone() /* done with calculation */ rcalfill(x, y) /* fill with ray-calculated pixel */ int x, y; { - if (queuesiz >= PACKSIZ) /* flush queue if needed */ - clearqueue(0); + if (queuesiz >= packsiz) /* flush queue if needed */ + clearqueue(); /* add position to queue */ queue[queuesiz][0] = x; queue[queuesiz][1] = y; @@ -690,46 +1151,64 @@ int x, y; } -clearqueue(done) /* process queue */ -int done; +clearqueue() /* process queue */ { FVECT orig, dir; - float fbuf[6]; + float fbuf[6*(PACKSIZ+1)]; + register float *fbp; register int i; + double vx, vy; + if (queuesiz == 0) + return; + fbp = fbuf; for (i = 0; i < queuesiz; i++) { viewray(orig, dir, &ourview, (queue[i][0]+.5)/hresolu, (queue[i][1]+.5)/vresolu); - fbuf[0] = orig[0]; fbuf[1] = orig[1]; fbuf[2] = orig[2]; - fbuf[3] = dir[0]; fbuf[4] = dir[1]; fbuf[5] = dir[2]; - fwrite((char *)fbuf, sizeof(float), 6, psend); + *fbp++ = orig[0]; *fbp++ = orig[1]; *fbp++ = orig[2]; + *fbp++ = dir[0]; *fbp++ = dir[1]; *fbp++ = dir[2]; } - if ((done ? fclose(psend) : fflush(psend)) == EOF) - syserror(); + /* mark end and get results */ + bzero((char *)fbp, 6*sizeof(float)); + if (process(PDesc, fbuf, fbuf, 4*sizeof(float)*(queuesiz+1), + 6*sizeof(float)*(queuesiz+1)) != + 4*sizeof(float)*(queuesiz+1)) { + fprintf(stderr, "%s: error reading from rtrace process\n", + progname); + exit(1); + } + fbp = fbuf; for (i = 0; i < queuesiz; i++) { - if (fread((char *)fbuf, sizeof(float), 4, precv) < 4) { - fprintf(stderr, "%s: read error in clearqueue\n", - progname); - exit(1); - } if (ourexp > 0 && ourexp != 1.0) { - fbuf[0] *= ourexp; - fbuf[1] *= ourexp; - fbuf[2] *= ourexp; + fbp[0] *= ourexp; + fbp[1] *= ourexp; + fbp[2] *= ourexp; } - setcolr(pscan(queue[i][1])[queue[i][0]], - fbuf[0], fbuf[1], fbuf[2]); - zscan(queue[i][1])[queue[i][0]] = fbuf[3]; + if (averaging) { + setcolor(sscan(queue[i][1])[queue[i][0]], + fbp[0], fbp[1], fbp[2]); + wscan(queue[i][1])[queue[i][0]] = 1; + } else + setcolr(pscan(queue[i][1])[queue[i][0]], + fbp[0], fbp[1], fbp[2]); + if (zisnorm) + zscan(queue[i][1])[queue[i][0]] = fbp[3]; + else { + vx = (queue[i][0]+.5)/hresolu + ourview.hoff - .5; + vy = (queue[i][1]+.5)/vresolu + ourview.voff - .5; + zscan(queue[i][1])[queue[i][0]] = fbp[3] / sqrt(1. + + vx*vx*ourview.hn2 + vy*vy*ourview.vn2); + } + fbp += 4; } - if (done) - fclose(precv); queuesiz = 0; } -syserror() /* report error and exit */ +syserror(s) /* report error and exit */ +char *s; { - perror(progname); + perror(s); exit(1); }