/* Copyright (c) 1994 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * Interpolate and extrapolate pictures with different view parameters.
 *
 *	Greg Ward	09Dec89
 */

#include "standard.h"

#include <ctype.h>

#include "view.h"

#include "color.h"

#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 averaging	(ourweigh != NULL)

#define MAXWT		1000.		/* maximum pixel weight (averaging) */

#define F_FORE		1		/* fill foreground */
#define F_BACK		2		/* fill background */

#define PACKSIZ		256		/* max. calculation packet size */

#define RTCOM		"rtrace -h- -ovl -fff "

#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 */
double	pixaspect = 1.0;		/* pixel aspect ratio */

double	zeps = .02;			/* allowed z epsilon */

COLR	*ourpict;			/* output picture (COLR's) */
COLOR	*ourspict;			/* output pixel sums (averaging) */
float	*ourweigh = NULL;		/* output pixel weights (averaging) */
float	*ourzbuf;			/* corresponding z-buffer */

char	*progname;

int	fillo = F_FORE|F_BACK;		/* selected fill options */
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;			/* original picture exposure */
int	expadj = 0;			/* exposure adjustment (f-stops) */
double	rexpadj = 1;			/* real exposure adjustment */

VIEW	theirview = STDVIEW;		/* input view */
int	gotview;			/* got input view? */
int	wrongformat = 0;		/* input in another format? */
RESOLU	tresolu;			/* input resolution */
double	theirexp;			/* input picture exposure */
MAT4	theirs2ours;			/* transformation matrix */
int	hasmatrix = 0;			/* has transformation matrix */

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;			/* number of pixels pending */

extern double	movepixel();


main(argc, argv)			/* interpolate pictures */
int	argc;
char	*argv[];
{
#define  check(ol,al)		if (argv[i][ol] || \
				badarg(argc-i-1,argv+i+1,al)) \
				goto badopt
	int	gotvfile = 0;
	int	doavg = -1;
	char	*zfile = NULL;
	char	*expcomp = NULL;
	char	*err;
	int	i, rval;

	progname = argv[0];

	for (i = 1; i < argc && argv[i][0] == '-'; i++) {
		rval = getviewopt(&ourview, argc-i, argv+i);
		if (rval >= 0) {
			i += rval;
			continue;
		}
		switch (argv[i][1]) {
		case 'e':				/* exposure */
			check(2,"f");
			expcomp = argv[++i];
			break;
		case 't':				/* threshold */
			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;
			break;
		case 'f':				/* fill type */
			switch (argv[i][2]) {
			case '0':				/* none */
				check(3,NULL);
				fillo = 0;
				break;
			case 'f':				/* foreground */
				check(3,NULL);
				fillo = F_FORE;
				break;
			case 'b':				/* background */
				check(3,NULL);
				fillo = F_BACK;
				break;
			case 'a':				/* all */
				check(3,NULL);
				fillo = F_FORE|F_BACK;
				break;
			case 's':				/* sample */
				check(3,"i");
				fillsamp = atoi(argv[++i]);
				break;
			case 'c':				/* color */
				check(3,"fff");
				fillfunc = backfill;
				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 */
				check(3,"f");
				fillfunc = backfill;
				backz = atof(argv[++i]);
				break;
			case 'r':				/* rtrace */
				check(3,"s");
				fillfunc = rcalfill;
				calstart(RTCOM, argv[++i]);
				break;
			default:
				goto badopt;
			}
			break;
		case 'z':				/* z file */
			check(2,"s");
			zfile = argv[++i];
			break;
		case 'x':				/* x resolution */
			check(2,"i");
			hresolu = atoi(argv[++i]);
			break;
		case 'y':				/* y resolution */
			check(2,"i");
			vresolu = atoi(argv[++i]);
			break;
		case 'p':				/* pixel aspect */
			if (argv[i][2] != 'a')
				goto badopt;
			check(3,"f");
			pixaspect = atof(argv[++i]);
			break;
		case 'v':				/* view file */
			if (argv[i][2] != 'f')
				goto badopt;
			check(3,"s");
			gotvfile = viewfile(argv[++i], &ourview, 0, 0);
			if (gotvfile < 0)
				syserror(argv[i]);
			else if (gotvfile == 0) {
				fprintf(stderr, "%s: bad view file\n",
						argv[i]);
				exit(1);
			}
			break;
		default:
		badopt:
			fprintf(stderr, "%s: command line error at '%s'\n",
					progname, argv[i]);
			goto userr;
		}
	}
						/* check arguments */
	if ((argc-i)%2)
		goto userr;
	if (fillsamp == 1)
		fillo &= ~F_BACK;
	if (doavg < 0)
		doavg = (argc-i) > 2;
	if (expcomp != NULL)
		if (expcomp[0] == '+' | expcomp[0] == '-') {
			expadj = atof(expcomp) + (expcomp[0]=='+' ? .5 : -.5);
			if (doavg)
				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)
				rexpadj = pow(2.0, (double)expadj);
		}
						/* set view */
	if ((err = setview(&ourview)) != NULL) {
		fprintf(stderr, "%s: %s\n", progname, err);
		exit(1);
	}
	normaspect(viewaspect(&ourview), &pixaspect, &hresolu, &vresolu);
						/* allocate frame */
	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 (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]);
							/* fill in spaces */
	if (fillo&F_BACK)
		backpicture(fillfunc, fillsamp);
	else
		fillpicture(fillfunc);
							/* close calculation */
	caldone();
							/* aft clipping */
	clipaft();
							/* add to header */
	printargs(argc, argv, stdout);
	if (gotvfile) {
		fputs(VIEWSTR, stdout);
		fprintview(&ourview, stdout);
		putc('\n', stdout);
	}
	if (pixaspect < .99 | pixaspect > 1.01)
		fputaspect(pixaspect, stdout);
	if (ourexp > 0)
		ourexp *= rexpadj;
	else
		ourexp = rexpadj;
	if (ourexp < .995 | ourexp > 1.005)
		fputexpos(ourexp, stdout);
	fputformat(COLRFMT, stdout);
	putc('\n', stdout);
							/* write picture */
	writepicture();
							/* write z file */
	if (zfile != NULL)
		writedistance(zfile);

	exit(0);
userr:
	fprintf(stderr,
	"Usage: %s [view opts][-t eps][-z zout][-e spec][-a|-q][-fT][-n] pfile zspec ..\n",
			progname);
	exit(1);
#undef check
}


headline(s)				/* process header string */
char	*s;
{
	char	fmt[32];

	if (isheadid(s))
		return;
	if (formatval(fmt, s)) {
		wrongformat = strcmp(fmt, COLRFMT);
		return;
	}
	putc('\t', stdout);
	fputs(s, stdout);

	if (isexpos(s)) {
		theirexp *= exposval(s);
		return;
	}
	if (isview(s) && sscanview(&theirview, s) > 0)
		gotview++;
}


addpicture(pfile, zspec)		/* add picture to output */
char	*pfile, *zspec;
{
	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)
		syserror(pfile);
					/* get header with exposure and view */
	theirexp = 1.0;
	gotview = 0;
	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)
		ourexp = theirexp;
	else if (ABS(theirexp-ourexp) > .01*ourexp)
		fprintf(stderr, "%s: different exposure (warning)\n", pfile);
	if (err = setview(&theirview)) {
		fprintf(stderr, "%s: %s\n", pfile, err);
		exit(1);
	}
					/* compute transformation */
	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 (!isflt(zspec) || (zvalue = atof(zspec)) <= 0.0)
			syserror(zspec);
		for (x = scanlen(&tresolu); x-- > 0; )
			zin[x] = zvalue;
	}
					/* 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 && 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);
		}
		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);
	fclose(pfp);
	if (zfd != -1)
		close(zfd);
}


pixform(xfmat, vw1, vw2)		/* compute view1 to view2 matrix */
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];
	xfmat[0][2] = vw1->hvec[2];
	xfmat[1][0] = vw1->vvec[0];
	xfmat[1][1] = vw1->vvec[1];
	xfmat[1][2] = vw1->vvec[2];
	xfmat[2][0] = vw1->vdir[0];
	xfmat[2][1] = vw1->vdir[1];
	xfmat[2][2] = vw1->vdir[2];
	xfmat[3][0] = vw1->vp[0];
	xfmat[3][1] = vw1->vp[1];
	xfmat[3][2] = vw1->vp[2];
	setident4(m4t);
	m4t[0][0] = vw2->hvec[0]/vw2->hn2;
	m4t[1][0] = vw2->hvec[1]/vw2->hn2;
	m4t[2][0] = vw2->hvec[2]/vw2->hn2;
	m4t[3][0] = -DOT(vw2->vp,vw2->hvec)/vw2->hn2;
	m4t[0][1] = vw2->vvec[0]/vw2->vn2;
	m4t[1][1] = vw2->vvec[1]/vw2->vn2;
	m4t[2][1] = vw2->vvec[2]/vw2->vn2;
	m4t[3][1] = -DOT(vw2->vp,vw2->vvec)/vw2->vn2;
	m4t[0][2] = vw2->vdir[0];
	m4t[1][2] = vw2->vdir[1];
	m4t[2][2] = vw2->vdir[2];
	m4t[3][2] = -DOT(vw2->vp,vw2->vdir);
	multmat4(xfmat, xfmat, m4t);
	return(1);
}


addscanline(xl, y, pline, zline, lasty)	/* add scanline to output */
struct bound	*xl;
int	y;
COLR	*pline;
float	*zline;
struct position	*lasty;		/* input/output */
{
	FVECT	pos;
	struct position	lastx, newpos;
	double	wt;
	register int	x;

	lastx.z = 0;
	for (x = xl->max; x >= xl->min; x--) {
		pix2loc(pos, &tresolu, x, y);
		pos[2] = zline[x];
		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];
		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, 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? */
	int	x1, y1;				/* p1 position */
	register int	x, y;			/* final position */

					/* compute vector p0p1 */
	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 (fillo&F_FORE && ABS(p2->z-p0->z) <= zt) {
		s2x = p2->x - p0->x;
		s2y = p2->y - p0->y;
		if (p1isy == 1)
			l2 = ABS(s2x);
		else {
			l2 = ABS(s2y);
			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 (y < 0 || y >= vresolu)
				continue;
			if (averaging) {
				if (zscan(y)[x] <= 0 || zscan(y)[x]-z
						> zeps*zscan(y)[x]) {
					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 */
register FVECT	pos;
{
	FVECT	pt, tdir, odir;
	double	d;
	register int	i;
	
	if (pos[2] <= 0)		/* empty pixel */
		return(0);
	if (!averaging && hasmatrix) {
		pos[0] += theirview.hoff - .5;
		pos[1] += theirview.voff - .5;
		if (theirview.type == VT_PER) {
			if (normdist)			/* adjust 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(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);
		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;
	register int	x, i;
							/* get back buffer */
	yback = (int *)malloc(hresolu*sizeof(int));
	if (yback == NULL)
		syserror(progname);
	for (x = 0; x < hresolu; x++)
		yback[x] = -2;
	/*
	 * Xback and yback are the pixel locations of suitable
	 * background values in each direction.
	 * A value of -2 means unassigned, and -1 means
	 * that there is no suitable background in this direction.
	 */
							/* fill image */
	for (y = 0; y < vresolu; y++) {
		xback = -2;
		for (x = 0; x < hresolu; x++)
			if (zscan(y)[x] <= 0) {		/* empty pixel */
				/*
				 * First, find background from above or below.
				 * (farthest assigned pixel)
				 */
				if (yback[x] == -2) {
					for (i = y+1; i < vresolu; i++)
						if (zscan(i)[x] > 0)
							break;
					if (i < vresolu
				&& (y <= 0 || zscan(y-1)[x] < zscan(i)[x]))
						yback[x] = i;
					else
						yback[x] = y-1;
				}
				/*
				 * Next, find background from left or right.
				 */
				if (xback == -2) {
					for (i = x+1; i < hresolu; i++)
						if (zscan(y)[i] > 0)
							break;
					if (i < hresolu
				&& (x <= 0 || zscan(y)[x-1] < zscan(y)[i]))
						xback = i;
					else
						xback = x-1;
				}
				/*
				 * If we have no background for this pixel,
				 * use the given fill function.
				 */
				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] ) ) {
					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 {
					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;
			}
	}
	free((char *)yback);
}


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)
				(*fill)(x,y);
}


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 (alters buffer) */
{
	int	y;
	register int	x;
	double	d;

	fprtresolu(hresolu, vresolu, stdout);
	for (y = vresolu-1; y >= 0; y--)
		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 */
char	*fname;
{
	int	donorm = normdist && ourview.type == VT_PER &&
			!averaging && hasmatrix;
	int	fd;
	int	y;
	float	*zout;

	if ((fd = open(fname, O_WRONLY|O_CREAT|O_TRUNC, 0666)) == -1)
		syserror(fname);
	if (donorm
	&& (zout = (float *)malloc(hresolu*sizeof(float))) == NULL)
		syserror(progname);
	for (y = vresolu-1; y >= 0; y--) {
		if (donorm) {
			double	vx, yzn2;
			register int	x;
			yzn2 = (y+.5)/vresolu + ourview.voff - .5;
			yzn2 = 1. + yzn2*yzn2*ourview.vn2;
			for (x = 0; x < hresolu; x++) {
				vx = (x+.5)/hresolu + ourview.hoff - .5;
				zout[x] = zscan(y)[x]
					* sqrt(vx*vx*ourview.hn2 + yzn2);
			}
		} else
			zout = zscan(y);
		if (write(fd, (char *)zout, hresolu*sizeof(float))
				< hresolu*sizeof(float))
			syserror(fname);
	}
	if (donorm)
		free((char *)zout);
	close(fd);
}


backfill(x, y)				/* fill pixel with background */
int	x, y;
{
	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 */
char	*prog, *args;
{
	char	combuf[512];
	char	*argv[64];
	int	rval;
	register char	**wp, *cp;

	if (childpid != -1) {
		fprintf(stderr, "%s: too many calculations\n", progname);
		exit(1);
	}
	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))
			;
	}
	*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 */
{
	if (childpid == -1)
		return;
	clearqueue();
	close_process(PDesc);
	childpid = -1;
}


rcalfill(x, y)				/* fill with ray-calculated pixel */
int	x, y;
{
	if (queuesiz >= packsiz)	/* flush queue if needed */
		clearqueue();
					/* add position to queue */
	queue[queuesiz][0] = x;
	queue[queuesiz][1] = y;
	queuesiz++;
}


clearqueue()				/* process queue */
{
	FVECT	orig, dir;
	float	fbuf[6*(PACKSIZ+1)];
	register float	*fbp;
	register int	i;

	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);
		*fbp++ = orig[0]; *fbp++ = orig[1]; *fbp++ = orig[2];
		*fbp++ = dir[0]; *fbp++ = dir[1]; *fbp++ = dir[2];
	}
					/* mark end and get results */
	bzero((char *)fbp, 6*sizeof(float));
	if (process(PDesc, fbuf, fbuf, 4*sizeof(float)*queuesiz,
			6*sizeof(float)*(queuesiz+1)) !=
			4*sizeof(float)*queuesiz) {
		fprintf(stderr, "%s: error reading from rtrace process\n",
				progname);
		exit(1);
	}
	fbp = fbuf;
	for (i = 0; i < queuesiz; i++) {
		if (ourexp > 0 && ourexp != 1.0) {
			fbp[0] *= ourexp;
			fbp[1] *= ourexp;
			fbp[2] *= ourexp;
		}
		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;
	}
	queuesiz = 0;
}


syserror(s)			/* report error and exit */
char	*s;
{
	perror(s);
	exit(1);
}