--- ray/src/px/pinterp.c	1989/12/12 15:53:10	1.3
+++ ray/src/px/pinterp.c	1990/01/03 12:52:30	1.9
@@ -21,7 +21,7 @@ static char SCCSid[] = "$SunId$ LBL";
 
 VIEW	ourview = STDVIEW(512);		/* desired view */
 
-double	zeps = 0.001;			/* allowed z epsilon */
+double	zeps = .02;			/* allowed z epsilon */
 
 COLR	*ourpict;			/* output picture */
 float	*ourzbuf;			/* corresponding z-buffer */
@@ -31,7 +31,9 @@ char	*progname;
 VIEW	theirview = STDVIEW(512);	/* input view */
 int	gotview;			/* got input view? */
 
+double	theirs2ours[4][4];		/* transformation matrix */
 
+
 main(argc, argv)			/* interpolate pictures */
 int	argc;
 char	*argv[];
@@ -49,6 +51,14 @@ char	*argv[];
 			check(2,1);
 			zeps = atof(argv[++i]);
 			break;
+		case 'x':				/* x resolution */
+			check(2,1);
+			ourview.hresolu = atoi(argv[++i]);
+			break;
+		case 'y':				/* y resolution */
+			check(2,1);
+			ourview.vresolu = atoi(argv[++i]);
+			break;
 		case 'v':				/* view */
 			switch (argv[i][2]) {
 			case 't':				/* type */
@@ -163,10 +173,10 @@ addpicture(pfile, zfile)		/* add picture to output */
 char	*pfile, *zfile;
 {
 	FILE	*pfp, *zfp;
-	COLR	*scanin;
-	float	*zin;
 	char	*err;
-	int	xres, yres;
+	COLR	*scanin;
+	float	*zin, *zlast;
+	int	*plast;
 	int	y;
 					/* open input files */
 	if ((pfp = fopen(pfile, "r")) == NULL) {
@@ -181,75 +191,170 @@ char	*pfile, *zfile;
 	printf("%s:\n", pfile);
 	gotview = 0;
 	getheader(pfp, headline);
-	if (!gotview || fgetresolu(&xres, &yres, pfp) != (YMAJOR|YDECR)) {
+	if (!gotview || fgetresolu(&theirview.hresolu, &theirview.vresolu, pfp)
+			!= (YMAJOR|YDECR)) {
 		fprintf(stderr, "%s: picture view error\n", pfile);
 		exit(1);
 	}
-	theirview.hresolu = xres;
-	theirview.vresolu = yres;
 	if (err = setview(&theirview)) {
 		fprintf(stderr, "%s: %s\n", pfile, err);
 		exit(1);
 	}
+					/* compute transformation */
+	pixform(theirs2ours, &theirview, &ourview);
 					/* allocate scanlines */
-	scanin = (COLR *)malloc(xres*sizeof(COLR));
-	zin = (float *)malloc(xres*sizeof(float));
-	if (scanin == NULL || zin == NULL) {
+	scanin = (COLR *)malloc(theirview.hresolu*sizeof(COLR));
+	zin = (float *)malloc(theirview.hresolu*sizeof(float));
+	zlast = (float *)calloc(theirview.hresolu, sizeof(float));
+	plast = (int *)calloc(theirview.hresolu, sizeof(int));
+	if (scanin == NULL || zin == NULL || zlast == NULL || plast == NULL) {
 		perror(progname);
 		exit(1);
 	}
 					/* load image */
-	for (y = yres-1; y >= 0; y--) {
-		if (freadcolrs(scanin, xres, pfp) < 0) {
+	for (y = theirview.vresolu-1; y >= 0; y--) {
+		if (freadcolrs(scanin, theirview.hresolu, pfp) < 0) {
 			fprintf(stderr, "%s: read error\n", pfile);
 			exit(1);
 		}
-		if (fread(zin, sizeof(float), xres, zfp) < xres) {
+		if (fread(zin, sizeof(float), theirview.hresolu, zfp)
+				< theirview.hresolu) {
 			fprintf(stderr, "%s: read error\n", zfile);
 			exit(1);
 		}
-		addscanline(y, scanin, zin);
+		addscanline(y, scanin, zin, plast, zlast);
 	}
 					/* clean up */
 	free((char *)scanin);
 	free((char *)zin);
+	free((char *)plast);
+	free((char *)zlast);
 	fclose(pfp);
 	fclose(zfp);
 }
 
 
-addscanline(y, pline, zline)		/* add scanline to output */
+pixform(xfmat, vw1, vw2)		/* compute view1 to view2 matrix */
+register double	xfmat[4][4];
+register VIEW	*vw1, *vw2;
+{
+	double	m4t[4][4];
+
+	setident4(xfmat);
+	xfmat[0][0] = vw1->vhinc[0];
+	xfmat[0][1] = vw1->vhinc[1];
+	xfmat[0][2] = vw1->vhinc[2];
+	xfmat[1][0] = vw1->vvinc[0];
+	xfmat[1][1] = vw1->vvinc[1];
+	xfmat[1][2] = vw1->vvinc[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->vhinc[0]/vw2->vhn2;
+	m4t[1][0] = vw2->vhinc[1]/vw2->vhn2;
+	m4t[2][0] = vw2->vhinc[2]/vw2->vhn2;
+	m4t[3][0] = -DOT(vw2->vp,vw2->vhinc)/vw2->vhn2;
+	m4t[0][1] = vw2->vvinc[0]/vw2->vvn2;
+	m4t[1][1] = vw2->vvinc[1]/vw2->vvn2;
+	m4t[2][1] = vw2->vvinc[2]/vw2->vvn2;
+	m4t[3][1] = -DOT(vw2->vp,vw2->vvinc)/vw2->vvn2;
+	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);
+}
+
+
+addscanline(y, pline, zline, lasty, lastyz)	/* add scanline to output */
 int	y;
 COLR	*pline;
 float	*zline;
+int	*lasty;			/* input/output */
+float	*lastyz;		/* input/output */
 {
-	FVECT	p, dir;
-	double	xnew, ynew, znew;
+	extern double	sqrt(), fabs();
+	double	pos[3];
+	int	lastx = 0;
+	double	lastxz = 0;
+	double  zt;
+	int	xpos, ypos;
 	register int	x;
-	register int	xpos, ypos;
 
-	for (x = 0; x < theirview.hresolu; x++) {
-		rayview(p, dir, &theirview, x+.5, y+.5);
-		p[0] += zline[x]*dir[0];
-		p[1] += zline[x]*dir[1];
-		p[2] += zline[x]*dir[2];
-		pixelview(&xnew, &ynew, &znew, &ourview, p);
-		if (znew <= 0.0 || xnew < 0 || xnew > ourview.hresolu
-				|| ynew < 0 || ynew > ourview.vresolu)
+	for (x = theirview.hresolu-1; x >= 0; x--) {
+		pos[0] = x - .5*(theirview.hresolu-1);
+		pos[1] = y - .5*(theirview.vresolu-1);
+		pos[2] = zline[x];
+		if (theirview.type == VT_PER) {
+			/*
+			 * The following (single) statement can go
+			 * if z is along the view direction rather
+			 * than an eye ray.
+			 */
+			pos[2] /= sqrt( 1.
+					+ pos[0]*pos[0]*theirview.vhn2
+					+ pos[1]*pos[1]*theirview.vvn2 );
+			pos[0] *= pos[2];
+			pos[1] *= pos[2];
+		}
+		multp3(pos, pos, theirs2ours);
+		if (pos[2] <= 0)
 			continue;
-					/* check current value at position */
-		xpos = xnew;
-		ypos = ynew;
-		if (zscan(ypos)[xpos] <= 0.0
-				|| zscan(ypos)[xpos] - znew
-					> zeps*zscan(ypos)[xpos]) {
-			zscan(ypos)[xpos] = znew;
-			copycolr(pscan(ypos)[xpos], pline[x]);
+		if (ourview.type == VT_PER) {
+			pos[0] /= pos[2];
+			pos[1] /= pos[2];
 		}
+		pos[0] += .5*ourview.hresolu;
+		pos[1] += .5*ourview.vresolu;
+		if (pos[0] < 0 || (xpos = pos[0]) >= ourview.hresolu
+			|| pos[1] < 0 || (ypos = pos[1]) >= ourview.vresolu)
+			continue;
+					/* add pixel to our image */
+		zt = 2.*zeps*zline[x];
+		addpixel(xpos, ypos,
+			(fabs(zline[x]-lastxz) <= zt) ? lastx - xpos : 1,
+			(fabs(zline[x]-lastyz[x]) <= zt) ? lasty[x] - ypos : 1,
+			pline[x], pos[2]);
+		lastx = xpos;
+		lasty[x] = ypos;
+		lastxz = lastyz[x] = zline[x];
 	}
 }
 
 
+addpixel(xstart, ystart, width, height, pix, z)	/* fill in area for pixel */
+int	xstart, ystart;
+int	width, height;
+COLR	pix;
+double	z;
+{
+	register int	x, y;
+					/* make width and height positive */
+	if (width < 0) {
+		width = -width;
+		xstart = xstart-width+1;
+	} else if (width == 0)
+		width = 1;
+	if (height < 0) {
+		height = -height;
+		ystart = ystart-height+1;
+	} else if (height == 0)
+		height = 1;
+					/* fill pixel(s) within rectangle */
+	for (y = ystart; y < ystart+height; y++)
+		for (x = xstart; x < xstart+width; x++)
+			if (zscan(y)[x] <= 0
+					|| zscan(y)[x]-z > zeps*zscan(y)[x]) {
+				zscan(y)[x] = z;
+				copycolr(pscan(y)[x], pix);
+			}
+}
+
+
 fillpicture()				/* fill in empty spaces */
 {
 	int	*yback, xback;
@@ -263,47 +368,62 @@ fillpicture()				/* fill in empty spaces */
 		return;
 	}
 	for (x = 0; x < ourview.hresolu; x++)
-		yback[x] = -1;
+		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 < ourview.vresolu; y++)
+	for (y = 0; y < ourview.vresolu; y++) {
+		xback = -2;
 		for (x = 0; x < ourview.hresolu; x++)
-			if (zscan(y)[x] <= 0.0) {	/* found hole */
-				xback = x-1;
-				do {			/* find boundary */
-					if (yback[x] < 0) {
-						for (i = y+1;
-							i < ourview.vresolu;
-								i++)
-							if (zscan(i)[x] > 0.0)
-								break;
-						if (i < ourview.vresolu
+			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 < ourview.vresolu; i++)
+						if (zscan(i)[x] > 0)
+							break;
+					if (i < ourview.vresolu
 				&& (y <= 0 || zscan(y-1)[x] < zscan(i)[x]))
-							yback[x] = i;
-						else
-							yback[x] = y-1;
-					}
-				} while (++x < ourview.hresolu
-						&& zscan(y)[x] <= 0.0);
-				i = xback;		/* pick background */
-				if (x < ourview.hresolu
-			&& (i < 0 || zscan(y)[i] < zscan(y)[x]))
-					xback = x;
-							/* fill hole */
-				if (xback < 0) {
-					while (++i < x)
-						if (yback[i] >= 0)
-				copycolr(pscan(y)[i],pscan(yback[i])[i]);
-				} else {
-					while (++i < x)
-						if (yback[i] < 0
-				|| ABS(i-xback) <= 1 || (ABS(y-yback[i]) > 1
-				&& zscan(yback[i])[i] < zscan(y)[xback]))
-					copycolr(pscan(y)[i],pscan(y)[xback]);
-						else
-					copycolr(pscan(y)[i],pscan(yback[i])[i]);
+						yback[x] = i;
+					else
+						yback[x] = y-1;
 				}
-		} else
-			yback[x] = -1;			/* clear boundary */
+				/*
+				 * Next, find background from left or right.
+				 */
+				if (xback == -2) {
+					for (i = x+1; x < ourview.hresolu; i++)
+						if (zscan(y)[i] > 0)
+							break;
+					if (i < ourview.hresolu
+				&& (x <= 0 || zscan(y)[x-1] < zscan(y)[i]))
+						xback = i;
+					else
+						xback = x-1;
+				}
+				if (xback < 0 && yback[x] < 0)
+					continue;	/* no background */
+				/*
+				 * Compare, and use the background that is
+				 * farther, unless one of them is next to us.
+				 */
+				if (yback[x] < 0 || ABS(x-xback) <= 1
+					|| ( ABS(y-yback[x]) > 1
+				&& zscan(yback[x])[x] < zscan(y)[xback] ))
+					copycolr(pscan(y)[x],pscan(y)[xback]);
+				else
+					copycolr(pscan(y)[x],pscan(yback[x])[x]);
+			} else {				/* full pixel */
+				yback[x] = -2;
+				xback = -2;
+			}
+	}
 	free((char *)yback);
 }