--- ray/src/px/pinterp.c 1993/04/01 11:21:18 2.10 +++ ray/src/px/pinterp.c 1994/12/23 22:35:12 2.19 @@ -1,4 +1,4 @@ -/* Copyright (c) 1993 Regents of the University of California */ +/* Copyright (c) 1994 Regents of the University of California */ #ifndef lint static char SCCSid[] = "$SunId$ LBL"; @@ -21,8 +21,13 @@ static char SCCSid[] = "$SunId$ LBL"; #include "resolu.h" #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 averaging (ourweigh != NULL) +#define MAXWT 100. /* maximum pixel weight (averaging) */ + #define F_FORE 1 /* fill foreground */ #define F_BACK 2 /* fill background */ @@ -34,6 +39,11 @@ static char SCCSid[] = "$SunId$ LBL"; 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,7 +51,9 @@ 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 (-a option) */ +float *ourweigh = NULL; /* output pixel weights (-a option) */ float *ourzbuf; /* corresponding z-buffer */ char *progname; @@ -51,6 +63,7 @@ int fillsamp = 0; /* sample separation (0 == inf) */ extern int backfill(), rcalfill(); /* fill functions */ 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 */ @@ -60,7 +73,7 @@ int gotview; /* got input view? */ 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 PDesc[3] = {-1,-1,-1}; /* rtrace process descriptor */ @@ -69,7 +82,9 @@ unsigned short queue[PACKSIZ][2]; /* pending pixels */ int packsiz; /* actual packet size */ int queuesiz; /* number of pixels pending */ +extern double movepixel(); + main(argc, argv) /* interpolate pictures */ int argc; char *argv[]; @@ -78,6 +93,7 @@ char *argv[]; badarg(argc-i-1,argv+i+1,al)) \ goto badopt int gotvfile = 0; + int doavg = -1; char *zfile = NULL; char *err; int i, rval; @@ -95,6 +111,14 @@ char *argv[]; check(2,"f"); zeps = atof(argv[++i]); break; + case 'a': /* average */ + check(2,NULL); + doavg = 1; + break; + case 'q': /* quick (no avg.) */ + check(2,NULL); + doavg = 0; + break; case 'n': /* dist. normalized? */ check(2,NULL); normdist = !normdist; @@ -124,8 +148,11 @@ char *argv[]; case 'c': /* color */ check(3,"fff"); fillfunc = backfill; - setcolr(backcolr, atof(argv[i+1]), + setcolor(backcolor, atof(argv[i+1]), atof(argv[i+2]), atof(argv[i+3])); + setcolr(backcolr, colval(backcolor,RED), + colval(backcolor,GRN), + colval(backcolor,BLU)); i += 3; break; case 'z': /* z value */ @@ -186,17 +213,30 @@ char *argv[]; if (fillsamp == 1) fillo &= ~F_BACK; /* set view */ - if (err = setview(&ourview)) { + if ((err = setview(&ourview)) != NULL) { fprintf(stderr, "%s: %s\n", progname, err); exit(1); } normaspect(viewaspect(&ourview), &pixaspect, &hresolu, &vresolu); /* allocate frame */ - ourpict = (COLR *)bmalloc(hresolu*vresolu*sizeof(COLR)); + if (doavg < 0) + doavg = (argc-i) > 2; + 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); + } ourzbuf = (float *)bmalloc(hresolu*vresolu*sizeof(float)); - if (ourpict == NULL || ourzbuf == NULL) + if (ourzbuf == NULL) syserror(progname); bzero((char *)ourzbuf, hresolu*vresolu*sizeof(float)); + /* new header */ + newheader("RADIANCE", stdout); /* get input */ for ( ; i < argc; i += 2) addpicture(argv[i], argv[i+1]); @@ -207,6 +247,8 @@ char *argv[]; fillpicture(fillfunc); /* close calculation */ caldone(); + /* aft clipping */ + clipaft(); /* add to header */ printargs(argc, argv, stdout); if (gotvfile) { @@ -241,8 +283,9 @@ char *s; { char fmt[32]; - if (isformat(s)) { - formatval(fmt, s); + if (isheadid(s)) + return; + if (formatval(fmt, s)) { wrongformat = strcmp(fmt, COLRFMT); return; } @@ -267,6 +310,7 @@ char *pfile, *zspec; COLR *scanin; float *zin; struct position *plast; + struct bound *xlim, ylim; int y; /* open picture file */ if ((pfp = fopen(pfile, "r")) == NULL) @@ -290,38 +334,59 @@ char *pfile, *zspec; } /* compute transformation */ hasmatrix = pixform(theirs2ours, &theirview, &ourview); - /* allocate scanlines */ - scanin = (COLR *)malloc(scanlen(&tresolu)*sizeof(COLR)); + /* get z specification or file */ zin = (float *)malloc(scanlen(&tresolu)*sizeof(float)); - plast = (struct position *)calloc(scanlen(&tresolu), - sizeof(struct position)); - if (scanin == NULL || zin == NULL || plast == NULL) + if (zin == NULL) syserror(progname); - /* get z specification or file */ if ((zfd = open(zspec, O_RDONLY)) == -1) { double zvalue; register int x; - if (!isfloat(zspec) || (zvalue = atof(zspec)) <= 0.0) + if (!isflt(zspec) || (zvalue = atof(zspec)) <= 0.0) syserror(zspec); - exit(1); for (x = scanlen(&tresolu); x-- > 0; ) zin[x] = zvalue; } - /* load image */ - for (y = 0; y < numscans(&tresolu); y++) { + /* 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, - scanlen(&tresolu)*sizeof(float)) - < scanlen(&tresolu)*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); @@ -332,7 +397,7 @@ 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]; @@ -372,7 +437,8 @@ register VIEW *vw1, *vw2; } -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; @@ -380,38 +446,37 @@ struct position *lasty; /* input/output */ { FVECT pos; struct position lastx, newpos; + double wt; register int x; lastx.z = 0; - for (x = scanlen(&tresolu); x-- > 0; ) { + for (x = xl->max; x >= xl->min; x--) { pix2loc(pos, &tresolu, x, y); pos[2] = zline[x]; - if (movepixel(pos) < 0) { + if ((wt = movepixel(pos)) <= FTINY) { lasty[x].z = lastx.z = 0; /* mark invalid */ continue; } + /* 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? */ @@ -425,6 +490,8 @@ double z; l1 = ABS(s1x); if (p1isy = (ABS(s1y) > l1)) l1 = ABS(s1y); + if (l1 < 1) + l1 = 1; } else { l1 = s1x = s1y = 1; p1isy = -1; @@ -440,9 +507,15 @@ 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; @@ -453,54 +526,165 @@ double z; y = y1 + c2*s2y/l2; if (y < 0 || y >= vresolu) continue; - if (zscan(y)[x] <= 0 || zscan(y)[x]-z + 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; } } } } +double movepixel(pos) /* reposition image point */ FVECT pos; { - FVECT pt, direc; + double d0, d1; + FVECT pt, tdir, odir; + register int i; if (pos[2] <= 0) /* empty pixel */ - return(-1); - if (hasmatrix) { + return(0); + if (normdist && theirview.type == VT_PER) { /* adjust distance */ + d0 = pos[0] + theirview.hoff - .5; + d1 = pos[1] + theirview.voff - .5; + pos[2] /= sqrt(1. + d0*d0*theirview.hn2 + d1*d1*theirview.vn2); + } + if (!averaging && hasmatrix) { pos[0] += theirview.hoff - .5; pos[1] += theirview.voff - .5; 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) - return(-1); + 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; + } else { + if (viewray(pt, tdir, &theirview, pos[0], pos[1]) < -FTINY) + return(0); + 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]; + d0 = DOT(odir,tdir); /* compute pixel weight */ + if (d0 <= FTINY) + return(0); /* relative angle >= 90 degrees */ + if (d0 >= 1.-1./MAXWT/MAXWT) + return(MAXWT); /* clip to maximum weight */ + return(1./sqrt(1.-d0)); +} + + +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); } - if (viewray(pt, direc, &theirview, pos[0], pos[1]) < 0) - return(-1); - pt[0] += direc[0]*pos[2]; - pt[1] += direc[1]*pos[2]; - pt[2] += direc[2]*pos[2]; - viewloc(pos, &ourview, pt); - if (pos[2] <= 0) - return(-1); - return(0); + 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; + for (x = numscans(&tresolu) - 1; x > y; x--) /* fill bottom rows */ + copystruct(xl+x, xl+y); + return(yl->max >= yl->min); } @@ -574,12 +758,25 @@ int samp; < zscan(y)[xback] ) ) { if (samp > 0 && ABS(x-xback) >= samp) goto fillit; - copycolr(pscan(y)[x],pscan(y)[xback]); + 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 { if (samp > 0 && ABS(y-yback[x]) > samp) goto fillit; - copycolr(pscan(y)[x],pscan(yback[x])[x]); + 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; @@ -611,13 +808,56 @@ int (*fill)(); } +clipaft() /* perform aft clipping as indicated */ +{ + register int x, y; + double tstdist; + double yzn2, vx; + + if (ourview.vaft <= FTINY) + return; + tstdist = ourview.vaft - ourview.vfore; + for (y = 0; y < vresolu; y++) { + if (ourview.type == VT_PER) { /* adjust distance */ + 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 (ourview.type == VT_PER) { + 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; + } + } +} + + writepicture() /* write out picture */ { int y; + register int x; + double d; fprtresolu(hresolu, vresolu, stdout); for (y = vresolu-1; y >= 0; y--) - if (fwritecolrs(pscan(y), hresolu, stdout) < 0) + if (averaging) { + for (x = 0; x < hresolu; x++) { /* average pixels */ + d = 1./wscan(y)[x]; + scalecolor(sscan(y)[x], d); + } + if (fwritescan(sscan(y), hresolu, stdout) < 0) + syserror(progname); + } else if (fwritecolrs(pscan(y), hresolu, stdout) < 0) syserror(progname); } @@ -658,23 +898,14 @@ char *fname; } -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; } @@ -696,14 +927,14 @@ char *prog, *args; cp = combuf; wp = argv; for ( ; ; ) { - while (isspace(*cp)) cp++; - if (!*cp) break; - *wp++ = cp; - while (!isspace(*cp)) - if (!*cp++) goto done; - *cp++ = '\0'; + while (isspace(*cp)) /* nullify spaces */ + *cp++ = '\0'; + if (!*cp) /* all done? */ + break; + *wp++ = cp; /* add argument to list */ + while (*++cp && !isspace(*cp)) + ; } -done: *wp = NULL; /* start process */ if ((rval = open_process(PDesc, argv)) < 0) @@ -774,8 +1005,13 @@ clearqueue() /* process queue */ fbp[1] *= ourexp; fbp[2] *= ourexp; } - setcolr(pscan(queue[i][1])[queue[i][0]], - fbp[0], fbp[1], fbp[2]); + 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]); zscan(queue[i][1])[queue[i][0]] = fbp[3]; fbp += 4; }