--- ray/src/rt/srcdraw.c	1996/03/15 21:26:42	2.1
+++ ray/src/rt/srcdraw.c	2003/06/26 00:58:10	2.7
@@ -1,13 +1,14 @@
-/* Copyright (c) 1996 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: srcdraw.c,v 2.7 2003/06/26 00:58:10 schorsch Exp $";
 #endif
-
 /*
  * Draw small sources into image in case we missed them.
+ *
+ *  External symbols declared in ray.h
  */
 
+#include "copyright.h"
+
 #include  "ray.h"
 
 #include  "view.h"
@@ -22,10 +23,21 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #define MAXVERT		10
 
+typedef struct splist {
+	struct splist	*next;			/* next source in list */
+	int	sn;				/* source number */
+	short	nv;				/* number of vertices */
+	RREAL	vl[3][2];			/* vertex array (last) */
+} SPLIST;				/* source polygon list */
 
+extern VIEW	ourview;		/* our view parameters */
+extern int	hres, vres;		/* our image resolution */
+static SPLIST	*sphead = NULL;		/* our list of source polys */
+
+
 static int
 inregion(p, cv, crit)			/* check if vertex is in region */
-FLOAT	p[2];
+RREAL	p[2];
 double	cv;
 int	crit;
 {
@@ -45,10 +57,10 @@ int	crit;
 
 static
 clipregion(a, b, cv, crit, r)		/* find intersection with boundary */
-register FLOAT	a[2], b[2];
+register RREAL	a[2], b[2];
 double	cv;
 int	crit;
-FLOAT	r[2];	/* return value */
+RREAL	r[2];	/* return value */
 {
 	switch (crit) {
 	case CLIP_ABOVE:
@@ -67,13 +79,13 @@ FLOAT	r[2];	/* return value */
 
 static int
 hp_clip_poly(vl, nv, cv, crit, vlo)	/* clip polygon to half-plane */
-FLOAT	vl[][2];
+RREAL	vl[][2];
 int	nv;
 double	cv;
 int	crit;
-FLOAT	vlo[][2];	/* return value */
+RREAL	vlo[][2];	/* return value */
 {
-	FLOAT	*s, *p;
+	RREAL	*s, *p;
 	register int	j, nvo;
 
 	s = vl[nv-1];
@@ -94,12 +106,12 @@ FLOAT	vlo[][2];	/* return value */
 
 static int
 box_clip_poly(vl, nv, xl, xr, yb, ya, vlo)	/* clip polygon to box */
-FLOAT	vl[MAXVERT][2];
+RREAL	vl[MAXVERT][2];
 int	nv;
 double	xl, xr, yb, ya;
-FLOAT	vlo[MAXVERT][2];	/* return value */
+RREAL	vlo[MAXVERT][2];	/* return value */
 {
-	FLOAT	vlt[MAXVERT][2];
+	RREAL	vlt[MAXVERT][2];
 	int	nvt, nvo;
 
 	nvt = hp_clip_poly(vl, nv, yb, CLIP_BELOW, vlt);
@@ -113,12 +125,12 @@ FLOAT	vlo[MAXVERT][2];	/* return value */
 
 static double
 minw2(vl, nv, ar2)			/* compute square of minimum width */
-FLOAT	vl[][2];
+RREAL	vl[][2];
 int	nv;
 double	ar2;
 {
 	double	d2, w2, w2min, w2max;
-	register FLOAT	*p0, *p1, *p2;
+	register RREAL	*p0, *p1, *p2;
 	int	i, j;
 				/* find minimum for all widths */
 	w2min = FHUGE;
@@ -146,9 +158,9 @@ double	ar2;
 
 static
 convex_center(vl, nv, cv)		/* compute center of convex polygon */
-register FLOAT	vl[][2];
+register RREAL	vl[][2];
 int	nv;
-FLOAT	cv[2];		/* return value */
+RREAL	cv[2];		/* return value */
 {
 	register int	i;
 					/* simple average (suboptimal) */
@@ -164,11 +176,11 @@ FLOAT	cv[2];		/* return value */
 
 static double
 poly_area(vl, nv)			/* compute area of polygon */
-register FLOAT	vl[][2];
+register RREAL	vl[][2];
 int	nv;
 {
 	double	a;
-	FLOAT	v0[2], v1[2];
+	RREAL	v0[2], v1[2];
 	register int	i;
 
 	a = 0.;
@@ -186,12 +198,12 @@ int	nv;
 
 static int
 convex_hull(vl, nv, vlo)		/* compute polygon's convex hull */
-FLOAT	vl[][2];
+RREAL	vl[][2];
 int	nv;
-FLOAT	vlo[][2];	/* return value */
+RREAL	vlo[][2];	/* return value */
 {
 	int	nvo, nvt;
-	FLOAT	vlt[MAXVERT][2];
+	RREAL	vlt[MAXVERT][2];
 	double	voa, vta;
 	register int	i, j;
 					/* start with original polygon */
@@ -220,11 +232,36 @@ FLOAT	vlo[][2];	/* return value */
 }
 
 
+static
+spinsert(sn, vl, nv)			/* insert new source polygon */
+int	sn;
+RREAL	vl[][2];
+int	nv;
+{
+	register SPLIST	*spn;
+	register int	i;
+
+	if (nv < 3)
+		return;
+	if (nv > 3)
+		spn = (SPLIST *)malloc(sizeof(SPLIST)+sizeof(RREAL)*2*(nv-3));
+	else
+		spn = (SPLIST *)malloc(sizeof(SPLIST));
+	if (spn == NULL)
+		error(SYSTEM, "out of memory in spinsert");
+	spn->sn = sn;
+	for (i = spn->nv = nv; i--; ) {
+		spn->vl[i][0] = vl[i][0]; spn->vl[i][1] = vl[i][1];
+	}
+	spn->next = sphead;		/* push onto global list */
+	sphead = spn;
+}
+
+
 int
-sourcepoly(vw, sn, sp)			/* compute image polygon for source */
-VIEW	*vw;
+sourcepoly(sn, sp)			/* compute image polygon for source */
 int	sn;
-FLOAT	sp[MAXVERT][2];
+RREAL	sp[MAXVERT][2];
 {
 	static char	cubeord[8][6] = {{1,3,2,6,4,5},{0,4,5,7,3,2},
 					 {0,1,3,7,6,4},{0,1,5,7,6,2},
@@ -232,16 +269,16 @@ FLOAT	sp[MAXVERT][2];
 					 {0,2,3,7,5,4},{1,5,4,6,2,3}};
 	register SRCREC	*s = source + sn;
 	FVECT	ap, ip;
-	FLOAT	pt[6][2];
+	RREAL	pt[6][2];
 	int	dir;
 	register int	i, j;
 
 	if (s->sflags & (SDISTANT|SFLAT)) {
-		if (s->sflags & SDISTANT && vw->type == VT_PAR)
+		if (s->sflags & SDISTANT && ourview.type == VT_PAR)
 			return(0);		/* all or nothing case */
 		if (s->sflags & SFLAT) {
 			for (i = 0; i < 3; i++)
-				ap[i] = s->sloc[i] - vw->vp[i];
+				ap[i] = s->sloc[i] - ourview.vp[i];
 			if (DOT(ap, s->snorm) >= 0.)
 				return(0);	/* source faces away */
 		}
@@ -253,11 +290,11 @@ FLOAT	sp[MAXVERT][2];
 				if (j==2|j==3) ap[i] += s->ss[SV][i];
 				else ap[i] -= s->ss[SV][i];
 				if (s->sflags & SDISTANT) {
-					ap[i] *= 1. + vw->vfore;
-					ap[i] += vw->vp[i];
+					ap[i] *= 1. + ourview.vfore;
+					ap[i] += ourview.vp[i];
 				}
 			}
-			viewloc(ip, vw, ap);		/* find image point */
+			viewloc(ip, &ourview, ap);	/* find image point */
 			if (ip[2] <= 0.)
 				return(0);		/* in front of view */
 			sp[j][0] = ip[0]; sp[j][1] = ip[1];
@@ -266,7 +303,7 @@ FLOAT	sp[MAXVERT][2];
 	}
 					/* identify furthest corner */
 	for (i = 0; i < 3; i++)
-		ap[i] = s->sloc[i] - vw->vp[i];
+		ap[i] = s->sloc[i] - ourview.vp[i];
 	dir =	(DOT(ap,s->ss[SU])>0.) |
 		(DOT(ap,s->ss[SV])>0.)<<1 |
 		(DOT(ap,s->ss[SW])>0.)<<2 ;
@@ -281,7 +318,7 @@ FLOAT	sp[MAXVERT][2];
 			if (cubeord[dir][j] & 4) ap[i] += s->ss[SW][i];
 			else ap[i] -= s->ss[SW][i];
 		}
-		viewloc(ip, vw, ap);		/* find image point */
+		viewloc(ip, &ourview, ap);	/* find image point */
 		if (ip[2] <= 0.)
 			return(0);		/* in front of view */
 		pt[j][0] = ip[0]; pt[j][1] = ip[1];
@@ -290,77 +327,103 @@ FLOAT	sp[MAXVERT][2];
 }
 
 
-			/* add sources smaller than rad to computed subimage */
-drawsources(vw, xr, yr, pic, zbf, x0, xsiz, y0, ysiz, rad)
-VIEW	*vw;				/* full image view */
-int	xr, yr;				/* full image dimensions */
+			/* initialize by finding sources smaller than rad */
+init_drawsources(rad)
+int	rad;				/* source sample size */
+{
+	RREAL	spoly[MAXVERT][2];
+	int	nsv;
+	register SPLIST	*sp;
+	register int	i;
+					/* free old source list if one */
+	for (sp = sphead; sp != NULL; sp = sphead) {
+		sphead = sp->next;
+		free((void *)sp);
+	}
+					/* loop through all sources */
+	for (i = nsources; i--; ) {
+					/* compute image polygon for source */
+		if (!(nsv = sourcepoly(i, spoly)))
+			continue;
+					/* clip to image boundaries */
+		if (!(nsv = box_clip_poly(spoly, nsv, 0., 1., 0., 1., spoly)))
+			continue;
+					/* big enough for standard sampling? */
+		if (minw2(spoly, nsv, ourview.vn2/ourview.hn2) >
+				(double)rad*rad/hres/hres)
+			continue;
+					/* OK, add to our list */
+		spinsert(i, spoly, nsv);
+	}
+}
+
+void			/* add sources smaller than rad to computed subimage */
+drawsources(pic, zbf, x0, xsiz, y0, ysiz)
 COLOR	*pic[];				/* subimage pixel value array */
 float	*zbf[];				/* subimage distance array (opt.) */
 int	x0, xsiz, y0, ysiz;		/* origin and size of subimage */
-int	rad;				/* source sample size */
 {
-	int	sn;
-	FLOAT	spoly[MAXVERT][2], ppoly[MAXVERT][2];
+	RREAL	spoly[MAXVERT][2], ppoly[MAXVERT][2];
 	int	nsv, npv;
-	int	xmin, xmax, ymin, ymax, x, y, i;
-	FLOAT	cxy[2];
-	double	pa;
+	int	xmin, xmax, ymin, ymax, x, y;
+	RREAL	cxy[2];
+	double	w;
 	RAY	sr;
-					/* loop through all sources */
-	for (sn = 0; sn < nsources; sn++) {
-					/* compute image polygon for source */
-		if (!(nsv = sourcepoly(vw, sn, spoly)))
-			continue;
+	register SPLIST	*sp;
+	register int	i;
+					/* check each source in our list */
+	for (sp = sphead; sp != NULL; sp = sp->next) {
 					/* clip source poly to subimage */
-		nsv = box_clip_poly(spoly, nsv,
-				(double)x0/xr, (double)(x0+xsiz)/xr,
-				(double)y0/yr, (double)(y0+ysiz)/yr, spoly);
+		nsv = box_clip_poly(sp->vl, sp->nv,
+				(double)x0/hres, (double)(x0+xsiz)/hres,
+				(double)y0/vres, (double)(y0+ysiz)/vres, spoly);
 		if (!nsv)
 			continue;
-					/* is it big so sampling found it? */
-		if (minw2(spoly, nsv, vw->vn2/vw->hn2) > (double)rad*rad/xr/xr)
-			continue;
 					/* find common subimage (BBox) */
 		xmin = x0 + xsiz; xmax = x0;
 		ymin = y0 + ysiz; ymax = y0;
 		for (i = 0; i < nsv; i++) {
-			if ((double)xmin/xr > spoly[i][0])
-				xmin = spoly[i][0]*xr + FTINY;
-			if ((double)xmax/xr < spoly[i][0])
-				xmax = spoly[i][0]*xr - FTINY;
-			if ((double)ymin/yr > spoly[i][1])
-				ymin = spoly[i][1]*yr + FTINY;
-			if ((double)ymax/yr < spoly[i][1])
-				ymax = spoly[i][1]*yr - FTINY;
+			if ((double)xmin/hres > spoly[i][0])
+				xmin = spoly[i][0]*hres + FTINY;
+			if ((double)xmax/hres < spoly[i][0])
+				xmax = spoly[i][0]*hres - FTINY;
+			if ((double)ymin/vres > spoly[i][1])
+				ymin = spoly[i][1]*vres + FTINY;
+			if ((double)ymax/vres < spoly[i][1])
+				ymax = spoly[i][1]*vres - FTINY;
 		}
 					/* evaluate each pixel in BBox */
 		for (y = ymin; y <= ymax; y++)
 			for (x = xmin; x <= xmax; x++) {
 							/* subarea for pixel */
 				npv = box_clip_poly(spoly, nsv,
-						(double)x/xr, (x+1.)/xr,
-						(double)y/yr, (y+1.)/yr, ppoly);
+						(double)x/hres, (x+1.)/hres,
+						(double)y/vres, (y+1.)/vres,
+						ppoly);
 				if (!npv)
 					continue;	/* no overlap */
 				convex_center(ppoly, npv, cxy);
-				if ((sr.rmax = viewray(sr.rorg, sr.rdir, vw,
-						cxy[0], cxy[1])) < -FTINY)
+				if ((sr.rmax = viewray(sr.rorg,sr.rdir,&ourview,
+						cxy[0],cxy[1])) < -FTINY)
 					continue;	/* not in view */
-				if (source[sn].sflags & SSPOT &&
-						spotout(sr, source[sn].sl.s))
+				if (source[sp->sn].sflags & SSPOT &&
+						spotout(&sr, source[sp->sn].sl.s))
 					continue;	/* outside spot */
 				rayorigin(&sr, NULL, SHADOW, 1.0);
-				sr.rsrc = sn;
+				sr.rsrc = sp->sn;
 				rayvalue(&sr);		/* compute value */
 				if (bright(sr.rcol) <= FTINY)
 					continue;	/* missed/blocked */
 							/* modify pixel */
 				if (zbf[y-y0] != NULL &&
-						sr.rt < zbf[y-y0][x-x0])
+						sr.rt < 0.999*zbf[y-y0][x-x0])
 					zbf[y-y0][x-x0] = sr.rt;
-				pa = poly_area(ppoly, npv);
-				scalecolor(sr.rcol, pa*xr*yr);
-				scalecolor(pic[y-y0][x-x0], (1.-pa*xr*yr));
+				else if (!bigdiff(sr.rcol, pic[y-y0][x-x0],
+						0.001))	/* source sample */
+					setcolor(pic[y-y0][x-x0], 0., 0., 0.);
+				w = poly_area(ppoly, npv) * hres * vres;
+				scalecolor(sr.rcol, w);
+				scalecolor(pic[y-y0][x-x0], 1.-w);
 				addcolor(pic[y-y0][x-x0], sr.rcol);
 			}
 	}