--- ray/src/rt/srcdraw.c	1996/03/15 21:26:42	2.1
+++ ray/src/rt/srcdraw.c	2004/03/30 16:13:01	2.11
@@ -1,17 +1,16 @@
-/* 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.11 2004/03/30 16:13:01 schorsch Exp $";
 #endif
-
 /*
  * Draw small sources into image in case we missed them.
+ *
+ *  External symbols declared in ray.h
  */
 
-#include  "ray.h"
+#include "copyright.h"
 
+#include  "ray.h"
 #include  "view.h"
-
 #include  "source.h"
 
 
@@ -22,12 +21,37 @@ 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(RREAL p[2], double cv, int crit);
+static void clipregion(RREAL a[2], RREAL b[2], double cv, int crit, RREAL r[2]);
+static int hp_clip_poly(RREAL vl[][2], int nv, double cv, int crit,
+		RREAL vlo[][2]);
+static int box_clip_poly(RREAL vl[MAXVERT][2], int nv,
+		double xl, double xr, double yb, double ya, RREAL vlo[MAXVERT][2]);
+static double minw2(RREAL vl[][2], int nv, double ar2);
+static void convex_center(RREAL vl[][2], int nv, RREAL cv[2]);
+static double poly_area(RREAL vl[][2], int nv);
+static int convex_hull(RREAL vl[][2], int nv, RREAL vlo[][2]);
+static void spinsert(int sn, RREAL vl[][2], int nv);
+static int sourcepoly(int sn, RREAL sp[MAXVERT][2]);
+
+
 static int
-inregion(p, cv, crit)			/* check if vertex is in region */
-FLOAT	p[2];
-double	cv;
-int	crit;
+inregion(			/* check if vertex is in region */
+	RREAL	p[2],
+	double	cv,
+	int	crit
+)
 {
 	switch (crit) {
 	case CLIP_ABOVE:
@@ -43,12 +67,14 @@ int	crit;
 }
 
 
-static
-clipregion(a, b, cv, crit, r)		/* find intersection with boundary */
-register FLOAT	a[2], b[2];
-double	cv;
-int	crit;
-FLOAT	r[2];	/* return value */
+static void
+clipregion(		/* find intersection with boundary */
+	register RREAL	a[2],
+	register RREAL	b[2],
+	double	cv,
+	int	crit,
+	RREAL	r[2]	/* return value */
+)
 {
 	switch (crit) {
 	case CLIP_ABOVE:
@@ -66,14 +92,15 @@ FLOAT	r[2];	/* return value */
 
 
 static int
-hp_clip_poly(vl, nv, cv, crit, vlo)	/* clip polygon to half-plane */
-FLOAT	vl[][2];
-int	nv;
-double	cv;
-int	crit;
-FLOAT	vlo[][2];	/* return value */
+hp_clip_poly(	/* clip polygon to half-plane */
+	RREAL	vl[][2],
+	int	nv,
+	double	cv,
+	int	crit,
+	RREAL	vlo[][2]	/* return value */
+)
 {
-	FLOAT	*s, *p;
+	RREAL	*s, *p;
 	register int	j, nvo;
 
 	s = vl[nv-1];
@@ -93,13 +120,17 @@ 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];
-int	nv;
-double	xl, xr, yb, ya;
-FLOAT	vlo[MAXVERT][2];	/* return value */
+box_clip_poly(	/* clip polygon to box */
+	RREAL	vl[MAXVERT][2],
+	int	nv,
+	double	xl,
+	double	xr,
+	double	yb,
+	double	ya,
+	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);
@@ -112,13 +143,14 @@ FLOAT	vlo[MAXVERT][2];	/* return value */
 
 
 static double
-minw2(vl, nv, ar2)			/* compute square of minimum width */
-FLOAT	vl[][2];
-int	nv;
-double	ar2;
+minw2(			/* compute square of minimum width */
+	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;
@@ -144,11 +176,12 @@ double	ar2;
 }
 
 
-static
-convex_center(vl, nv, cv)		/* compute center of convex polygon */
-register FLOAT	vl[][2];
-int	nv;
-FLOAT	cv[2];		/* return value */
+static void
+convex_center(		/* compute center of convex polygon */
+	register RREAL	vl[][2],
+	int	nv,
+	RREAL	cv[2]		/* return value */
+)
 {
 	register int	i;
 					/* simple average (suboptimal) */
@@ -163,12 +196,13 @@ FLOAT	cv[2];		/* return value */
 
 
 static double
-poly_area(vl, nv)			/* compute area of polygon */
-register FLOAT	vl[][2];
-int	nv;
+poly_area(			/* compute area of polygon */
+	register RREAL	vl[][2],
+	int	nv
+)
 {
 	double	a;
-	FLOAT	v0[2], v1[2];
+	RREAL	v0[2], v1[2];
 	register int	i;
 
 	a = 0.;
@@ -185,13 +219,14 @@ int	nv;
 
 
 static int
-convex_hull(vl, nv, vlo)		/* compute polygon's convex hull */
-FLOAT	vl[][2];
-int	nv;
-FLOAT	vlo[][2];	/* return value */
+convex_hull(		/* compute polygon's convex hull */
+	RREAL	vl[][2],
+	int	nv,
+	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,44 +255,71 @@ FLOAT	vlo[][2];	/* return value */
 }
 
 
-int
-sourcepoly(vw, sn, sp)			/* compute image polygon for source */
-VIEW	*vw;
-int	sn;
-FLOAT	sp[MAXVERT][2];
+static void
+spinsert(			/* insert new source polygon */
+	int	sn,
+	RREAL	vl[][2],
+	int	nv
+)
 {
-	static char	cubeord[8][6] = {{1,3,2,6,4,5},{0,4,5,7,3,2},
+	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;
+}
+
+
+static int
+sourcepoly(			/* compute image polygon for source */
+	int	sn,
+	RREAL	sp[MAXVERT][2]
+)
+{
+	static short	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},
 					 {0,2,6,7,5,1},{0,4,6,7,3,1},
 					 {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 */
 		}
 		for (j = 0; j < 4; j++) {	/* four corners */
 			for (i = 0; i < 3; i++) {
 				ap[i] = s->sloc[i];
-				if (j==1|j==2) ap[i] += s->ss[SU][i];
+				if ((j==1)|(j==2)) ap[i] += s->ss[SU][i];
 				else ap[i] -= s->ss[SU][i];
-				if (j==2|j==3) ap[i] += s->ss[SV][i];
+				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 +328,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 +343,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 +352,112 @@ 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 */
-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 */
+			/* initialize by finding sources smaller than rad */
+extern void
+init_drawsources(
+	int	rad				/* source sample size */
+)
 {
-	int	sn;
-	FLOAT	spoly[MAXVERT][2], ppoly[MAXVERT][2];
-	int	nsv, npv;
-	int	xmin, xmax, ymin, ymax, x, y, i;
-	FLOAT	cxy[2];
-	double	pa;
-	RAY	sr;
+	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 (sn = 0; sn < nsources; sn++) {
+	for (i = nsources; i--; ) {
 					/* compute image polygon for source */
-		if (!(nsv = sourcepoly(vw, sn, spoly)))
+		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);
+	}
+}
+
+extern void			/* add sources smaller than rad to computed subimage */
+drawsources(
+	COLOR	*pic[],				/* subimage pixel value array */
+	float	*zbf[],				/* subimage distance array (opt.) */
+	int	x0,				/* origin and size of subimage */
+	int	xsiz,
+	int	y0,
+	int	ysiz
+)
+{
+	RREAL	spoly[MAXVERT][2], ppoly[MAXVERT][2];
+	int	nsv, npv;
+	int	xmin, xmax, ymin, ymax, x, y;
+	RREAL	cxy[2];
+	double	w;
+	RAY	sr;
+	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;
+				w = poly_area(ppoly, npv) * hres * vres;
+				if (w < .95) {		/* subpixel source */
+					rayorigin(&sr, NULL, SHADOW, 1.0);
+					sr.rsrc = sp->sn;
+				} else
+					rayorigin(&sr, NULL, PRIMARY, 1.0);
 				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.99*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.01))	/* source sample */
+					setcolor(pic[y-y0][x-x0], 0., 0., 0.);
+				scalecolor(sr.rcol, w);
+				scalecolor(pic[y-y0][x-x0], 1.-w);
 				addcolor(pic[y-y0][x-x0], sr.rcol);
 			}
 	}