--- ray/src/common/image.c	2018/01/24 03:34:15	2.41
+++ ray/src/common/image.c	2022/01/13 00:26:09	2.54
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: image.c,v 2.41 2018/01/24 03:34:15 greg Exp $";
+static const char	RCSid[] = "$Id: image.c,v 2.54 2022/01/13 00:26:09 greg Exp $";
 #endif
 /*
  *  image.c - routines for image generation.
@@ -15,11 +15,6 @@ static const char	RCSid[] = "$Id: image.c,v 2.41 2018/
 #include  "paths.h"
 #include  "view.h"
 
-
-#define  FEQ(x,y)	(fabs((x)-(y)) <= FTINY)
-#define  VEQ(v,w)	(FEQ((v)[0],(w)[0]) && FEQ((v)[1],(w)[1]) \
-				&& FEQ((v)[2],(w)[2]))
-
 VIEW  stdview = STDVIEW;		/* default view parameters */
 
 static gethfunc gethview;
@@ -68,8 +63,8 @@ VIEW  *v
 			return(ill_horiz);
 		if (v->vert >= 180.0-FTINY)
 			return(ill_vert);
-		v->hn2 = 2.0 * tan(v->horiz*(PI/180.0/2.0));
-		v->vn2 = 2.0 * tan(v->vert*(PI/180.0/2.0));
+		v->hn2 = 2.0 * tan(v->horiz*(PI/360.));
+		v->vn2 = 2.0 * tan(v->vert*(PI/360.));
 		break;
 	case VT_CYL:				/* cylindrical panorama */
 		if (v->horiz > 360.0+FTINY)
@@ -77,7 +72,7 @@ VIEW  *v
 		if (v->vert >= 180.0-FTINY)
 			return(ill_vert);
 		v->hn2 = v->horiz * (PI/180.0);
-		v->vn2 = 2.0 * tan(v->vert*(PI/180.0/2.0));
+		v->vn2 = 2.0 * tan(v->vert*(PI/360.));
 		break;
 	case VT_ANG:				/* angular fisheye */
 		if (v->horiz > 360.0+FTINY)
@@ -92,18 +87,18 @@ VIEW  *v
 			return(ill_horiz);
 		if (v->vert > 180.0+FTINY)
 			return(ill_vert);
-		v->hn2 = 2.0 * sin(v->horiz*(PI/180.0/2.0));
-		v->vn2 = 2.0 * sin(v->vert*(PI/180.0/2.0));
+		v->hn2 = 2.0 * sin(v->horiz*(PI/360.));
+		v->vn2 = 2.0 * sin(v->vert*(PI/360.));
 		break;
 	case VT_PLS:				/* planispheric fisheye */
 		if (v->horiz >= 360.0-FTINY)
 			return(ill_horiz);
 		if (v->vert >= 360.0-FTINY)
 			return(ill_vert);
-		v->hn2 = 2.*sin(v->horiz*(PI/180.0/2.0)) /
-				(1.0 + cos(v->horiz*(PI/180.0/2.0)));
-		v->vn2 = 2.*sin(v->vert*(PI/180.0/2.0)) /
-				(1.0 + cos(v->vert*(PI/180.0/2.0)));
+		v->hn2 = 2.*sin(v->horiz*(PI/360.)) /
+				(1.0 + cos(v->horiz*(PI/360.)));
+		v->vn2 = 2.*sin(v->vert*(PI/360.)) /
+				(1.0 + cos(v->vert*(PI/360.)));
 		break;
 	default:
 		return("unknown view type");
@@ -125,6 +120,79 @@ VIEW  *v
 }
 
 
+char *
+cropview(			/* crop a view to the indicated bounds */
+VIEW *v,
+double x0,
+double y0,
+double x1,
+double y1
+)
+{
+	static char	ill_hemi[] = "illegal crop for hemispherical view";
+	double	d;
+					/* order crop extrema */
+	if (x0 > x1) { d=x0; x0=x1; x1=d; }
+	if (y0 > y1) { d=y0; y0=y1; y1=d; }
+
+	if ((x1-x0 <= FTINY) | (y1-y0 <= FTINY))
+		return("zero crop area");
+
+	d = x1 - x0;			/* adjust horizontal size? */
+	if (!FABSEQ(d, 1.))
+		switch (v->type) {
+		case VT_PER:
+			v->horiz = 360./PI*atan( d*tan(PI/360.*v->horiz) );
+			break;
+		case VT_PAR:
+		case VT_ANG:
+		case VT_CYL:
+			v->horiz *= d;
+			break;
+		case VT_HEM:
+			d *= sin(PI/360.*v->horiz);
+			if (d > 1.)
+				return(ill_hemi);
+			v->horiz = 360./PI*asin( d );
+			break;
+		case VT_PLS:
+			d *= sin(PI/360.*v->horiz) /
+					(1. + cos(PI/360.*v->horiz));
+			v->horiz = 360./PI*acos( (1. - d*d) / (1. + d*d) );
+			break;
+		}
+
+	d = y1 - y0;			/* adjust vertical size? */
+	if (!FABSEQ(d, 1.))
+		switch (v->type) {
+		case VT_PER:
+		case VT_CYL:
+			v->vert = 360./PI*atan( d*tan(PI/360.*v->vert) );
+			break;
+		case VT_PAR:
+		case VT_ANG:
+			v->vert *= d;
+			break;
+		case VT_HEM:
+			d *= sin(PI/360.*v->vert);
+			if (d > 1.)
+				return(ill_hemi);
+			v->vert = 360./PI*asin( d );
+			break;
+		case VT_PLS:
+			d *= sin(PI/360.*v->vert) /
+					(1. + cos(PI/360.*v->vert));
+			v->vert = 360./PI*acos( (1. - d*d) / (1. + d*d) );
+			break;
+		}
+					/* adjust offsets */
+	v->hoff = ((x0 + x1)*.5 - .5 + v->hoff) / (x1 - x0);
+	v->voff = ((y0 + y1)*.5 - .5 + v->voff) / (y1 - y0);
+
+	return(setview(v));		/* final error checks & set-up */
+}
+
+
 void
 normaspect(				/* fix pixel aspect or resolution */
 double  va,			/* view aspect ratio */
@@ -231,8 +299,9 @@ viewloc(			/* find image location for point */
 FVECT  ip,
 VIEW  *v,
 FVECT  p
-)
+)	/* Use VL_* flags to interpret return value */
 {
+	int	rflags = VL_GOOD;
 	double  d, d2;
 	FVECT  disp;
 
@@ -244,17 +313,19 @@ FVECT  p
 		break;
 	case VT_PER:			/* perspective view */
 		d = DOT(disp,v->vdir);
+		if ((v->vaft > FTINY) & (d >= v->vaft))
+			rflags |= VL_BEYOND;
 		ip[2] = VLEN(disp);
-		if (d < 0.0) {		/* fold pyramid */
+		if (d < -FTINY) {	/* fold pyramid */
 			ip[2] = -ip[2];
 			d = -d;
-		}
-		if (d <= FTINY)
-			return(0);	/* at infinite edge */
+		} else if (d <= FTINY)
+			return(VL_BAD);	/* at infinite edge */
 		d = 1.0/d;
 		disp[0] *= d;
 		disp[1] *= d;
 		disp[2] *= d;
+		if (ip[2] < 0.0) d = -d;
 		ip[2] *= (1.0 - v->vfore*d);
 		break;
 	case VT_HEM:			/* hemispherical fisheye */
@@ -269,45 +340,57 @@ FVECT  p
 		d = DOT(disp,v->hvec);
 		d2 = DOT(disp,v->vdir);
 		ip[0] = 180.0/PI * atan2(d,d2) / v->horiz + 0.5 - v->hoff;
-		d = d*d + d2*d2;
-		if (d <= FTINY*FTINY)
-			return(0);	/* at pole */
-		d = 1.0/sqrt(d);
+		d2 = d*d + d2*d2;
+		if (d2 <= FTINY*FTINY)
+			return(VL_BAD);	/* at pole */
+		if ((v->vaft > FTINY) & (d2 >= v->vaft*v->vaft))
+			rflags |= VL_BEYOND;
+		d = 1.0/sqrt(d2);
 		ip[1] = DOT(disp,v->vvec)*d/v->vn2 + 0.5 - v->voff;
 		ip[2] = VLEN(disp);
 		ip[2] *= (1.0 - v->vfore*d);
-		return(1);
+		goto gotall;
 	case VT_ANG:			/* angular fisheye */
 		ip[0] = 0.5 - v->hoff;
 		ip[1] = 0.5 - v->voff;
 		ip[2] = normalize(disp) - v->vfore;
 		d = DOT(disp,v->vdir);
 		if (d >= 1.0-FTINY)
-			return(1);
+			goto gotall;
 		if (d <= -(1.0-FTINY)) {
 			ip[0] += 180.0/v->horiz;
-			return(1);
+			goto gotall;
 		}
 		d = (180.0/PI)*acos(d) / sqrt(1.0 - d*d);
 		ip[0] += DOT(disp,v->hvec)*d/v->horiz;
 		ip[1] += DOT(disp,v->vvec)*d/v->vert;
-		return(1);
+		goto gotall;
 	case VT_PLS:			/* planispheric fisheye */
 		ip[0] = 0.5 - v->hoff;
 		ip[1] = 0.5 - v->voff;
 		ip[2] = normalize(disp) - v->vfore;
 		d = DOT(disp,v->vdir);
 		if (d >= 1.0-FTINY)
-			return(1);
+			goto gotall;
 		if (d <= -(1.0-FTINY))
-			return(0);
+			return(VL_BAD);
 		ip[0] += DOT(disp,v->hvec)/((1. + d)*sqrt(v->hn2));
 		ip[1] += DOT(disp,v->vvec)/((1. + d)*sqrt(v->vn2));
-		return(1);
+		goto gotall;
+	default:
+		return(VL_BAD);
 	}
 	ip[0] = DOT(disp,v->hvec)/v->hn2 + 0.5 - v->hoff;
 	ip[1] = DOT(disp,v->vvec)/v->vn2 + 0.5 - v->voff;
-	return(1);
+gotall:					/* add appropriate return flags */
+	if (ip[2] <= 0.0)
+		rflags |= VL_BEHIND;
+	else if ((v->type != VT_PER) & (v->type != VT_CYL))
+		rflags |= VL_BEYOND*((v->vaft > FTINY) &
+					(ip[2] >= v->vaft - v->vfore));
+	rflags |= VL_OUTSIDE*((0.0 >= ip[0]) | (ip[0] >= 1.0) |
+				(0.0 >= ip[1]) | (ip[1] >= 1.0));
+	return(rflags);
 }
 
 
@@ -371,14 +454,14 @@ int  ac,
 char  *av[]
 )
 {
-#define check(c,l)	if ((av[0][c]&&av[0][c]!=' ') || \
+#define check(c,l)	if ((av[0][c]&&!isspace(av[0][c])) || \
 			badarg(ac-1,av+1,l)) return(-1)
 
 	if (ac <= 0 || av[0][0] != '-' || av[0][1] != 'v')
 		return(-1);
 	switch (av[0][2]) {
 	case 't':			/* type */
-		if (!av[0][3] || av[0][3]==' ')
+		if (!av[0][3] || isspace(av[0][3]))
 			return(-1);
 		check(4,"");
 		v->type = av[0][3];
@@ -499,44 +582,44 @@ VIEW  *vp
 		sprintf(cp, " -vt%c", vp->type);
 		cp += strlen(cp);
 	}
-	if (!VEQ(vp->vp,stdview.vp)) {
+	if (!VABSEQ(vp->vp,stdview.vp)) {
 		sprintf(cp, " -vp %.6g %.6g %.6g",
 				vp->vp[0], vp->vp[1], vp->vp[2]);
 		cp += strlen(cp);
 	}
-	if (!FEQ(vp->vdist,stdview.vdist) || !VEQ(vp->vdir,stdview.vdir)) {
+	if (!FABSEQ(vp->vdist,stdview.vdist) || !VABSEQ(vp->vdir,stdview.vdir)) {
 		sprintf(cp, " -vd %.6g %.6g %.6g",
 				vp->vdir[0]*vp->vdist,
 				vp->vdir[1]*vp->vdist,
 				vp->vdir[2]*vp->vdist);
 		cp += strlen(cp);
 	}
-	if (!VEQ(vp->vup,stdview.vup)) {
+	if (!VABSEQ(vp->vup,stdview.vup)) {
 		sprintf(cp, " -vu %.6g %.6g %.6g",
 				vp->vup[0], vp->vup[1], vp->vup[2]);
 		cp += strlen(cp);
 	}
-	if (!FEQ(vp->horiz,stdview.horiz)) {
+	if (!FABSEQ(vp->horiz,stdview.horiz)) {
 		sprintf(cp, " -vh %.6g", vp->horiz);
 		cp += strlen(cp);
 	}
-	if (!FEQ(vp->vert,stdview.vert)) {
+	if (!FABSEQ(vp->vert,stdview.vert)) {
 		sprintf(cp, " -vv %.6g", vp->vert);
 		cp += strlen(cp);
 	}
-	if (!FEQ(vp->vfore,stdview.vfore)) {
+	if (!FABSEQ(vp->vfore,stdview.vfore)) {
 		sprintf(cp, " -vo %.6g", vp->vfore);
 		cp += strlen(cp);
 	}
-	if (!FEQ(vp->vaft,stdview.vaft)) {
+	if (!FABSEQ(vp->vaft,stdview.vaft)) {
 		sprintf(cp, " -va %.6g", vp->vaft);
 		cp += strlen(cp);
 	}
-	if (!FEQ(vp->hoff,stdview.hoff)) {
+	if (!FABSEQ(vp->hoff,stdview.hoff)) {
 		sprintf(cp, " -vs %.6g", vp->hoff);
 		cp += strlen(cp);
 	}
-	if (!FEQ(vp->voff,stdview.voff)) {
+	if (!FABSEQ(vp->voff,stdview.voff)) {
 		sprintf(cp, " -vl %.6g", vp->voff);
 		cp += strlen(cp);
 	}
@@ -563,7 +646,7 @@ char  *s
 	}
 					/* skip leading path */
 	cp = s;
-	while (*cp && *cp != ' ')
+	while (*cp && !isspace(*cp))
 		cp++;
 	while (cp > s && !ISDIRSEP(cp[-1]))
 		cp--;
@@ -615,7 +698,8 @@ RESOLU  *rp
 	if (rp != NULL && !fgetsresolu(rp, fp))
 		mvs.ok = 0;
 
-	fclose(fp);
+	if (fp != stdin)
+		fclose(fp);
 
 	return(mvs.ok);
 }