--- ray/src/gen/mkillum2.c	1991/07/25 10:59:07	1.5
+++ ray/src/gen/mkillum2.c	2007/11/05 23:40:26	2.23
@@ -1,69 +1,222 @@
-/* Copyright (c) 1991 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: mkillum2.c,v 2.23 2007/11/05 23:40:26 greg Exp $";
 #endif
-
 /*
  * Routines to do the actual calculation for mkillum
  */
 
-#include  "mkillum.h"
+#include <string.h>
 
+#include  "mkillum.h"
 #include  "face.h"
-
 #include  "cone.h"
-
 #include  "random.h"
 
 
-o_default(ob, il, rt, nm)	/* default illum action */
-OBJREC  *ob;
-struct illum_args  *il;
-struct rtproc  *rt;
-char  *nm;
+COLORV *	distarr = NULL;		/* distribution array */
+int		distsiz = 0;
+
+
+void
+newdist(			/* allocate & clear distribution array */
+	int siz
+)
 {
+	if (siz == 0) {
+		if (distsiz > 0)
+			free((void *)distarr);
+		distarr = NULL;
+		distsiz = 0;
+		return;
+	}
+	if (distsiz < siz) {
+		if (distsiz > 0)
+			free((void *)distarr);
+		distarr = (COLORV *)malloc(sizeof(COLOR)*siz);
+		if (distarr == NULL)
+			error(SYSTEM, "out of memory in newdist");
+		distsiz = siz;
+	}
+	memset(distarr, '\0', sizeof(COLOR)*siz);
+}
+
+
+int
+process_ray(RAY *r, int rv)
+{
+	COLORV	*colp;
+
+	if (rv == 0)
+		return(0);
+	if (rv < 0)
+		error(USER, "ray tracing process died");
+	if (r->rno >= distsiz)
+		error(INTERNAL, "bad returned index in process_ray");
+	colp = &distarr[r->rno * 3];
+	addcolor(colp, r->rcol);
+	return(1);
+}
+
+
+void
+raysamp(	/* queue a ray sample */
+	int  ndx,
+	FVECT  org,
+	FVECT  dir
+)
+{
+	RAY	myRay;
+	int	rv;
+
+	if ((ndx < 0) | (ndx >= distsiz))
+		error(INTERNAL, "bad index in raysamp");
+	VCOPY(myRay.rorg, org);
+	VCOPY(myRay.rdir, dir);
+	myRay.rmax = .0;
+	rayorigin(&myRay, PRIMARY, NULL, NULL);
+	myRay.rno = ndx;
+					/* queue ray, check result */
+	process_ray(&myRay, ray_pqueue(&myRay));
+}
+
+
+void
+rayclean()			/* finish all pending rays */
+{
+	RAY	myRay;
+
+	while (process_ray(&myRay, ray_presult(&myRay, 0)))
+		;
+}
+
+
+static void
+mkaxes(			/* compute u and v to go with n */
+	FVECT  u,
+	FVECT  v,
+	FVECT  n
+)
+{
+	register int  i;
+
+	v[0] = v[1] = v[2] = 0.0;
+	for (i = 0; i < 3; i++)
+		if (n[i] < 0.6 && n[i] > -0.6)
+			break;
+	v[i] = 1.0;
+	fcross(u, v, n);
+	normalize(u);
+	fcross(v, n, u);
+}
+
+
+static void
+rounddir(		/* compute uniform spherical direction */
+	register FVECT  dv,
+	double  alt,
+	double  azi
+)
+{
+	double  d1, d2;
+
+	dv[2] = 1. - 2.*alt;
+	d1 = sqrt(1. - dv[2]*dv[2]);
+	d2 = 2.*PI * azi;
+	dv[0] = d1*cos(d2);
+	dv[1] = d1*sin(d2);
+}
+
+
+static void
+flatdir(		/* compute uniform hemispherical direction */
+	register FVECT  dv,
+	double  alt,
+	double  azi
+)
+{
+	double  d1, d2;
+
+	d1 = sqrt(alt);
+	d2 = 2.*PI * azi;
+	dv[0] = d1*cos(d2);
+	dv[1] = d1*sin(d2);
+	dv[2] = sqrt(1. - alt);
+}
+
+
+int
+my_default(	/* default illum action */
+	OBJREC  *ob,
+	struct illum_args  *il,
+	char  *nm
+)
+{
 	sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
 			nm, ofun[ob->otype].funame, ob->oname);
 	error(WARNING, errmsg);
-	if (!(il->flags & IL_LIGHT))
-		printobj(il->altmat, ob);
+	printobj(il->altmat, ob);
+	return(1);
 }
 
 
-o_face(ob, il, rt, nm)		/* make an illum face */
-OBJREC  *ob;
-struct illum_args  *il;
-struct rtproc  *rt;
-char  *nm;
+int
+my_face(		/* make an illum face */
+	OBJREC  *ob,
+	struct illum_args  *il,
+	char  *nm
+)
 {
 #define MAXMISS		(5*n*il->nsamps)
-	int  dim[4];
-	int  n, nalt, nazi;
-	float  *distarr;
-	double  r1, r2;
+	int  dim[2];
+	int  n, nalt, nazi, h, alti;
+	double  sp[2], r1, r2;
 	FVECT  dn, org, dir;
 	FVECT  u, v;
 	double  ur[2], vr[2];
+	MAT4  xfm;
 	int  nmisses;
-	register FACE  *fa;
+	FACE  *fa;
 	register int  i, j;
 				/* get/check arguments */
 	fa = getface(ob);
 	if (fa->area == 0.0) {
 		freeface(ob);
-		o_default(ob, il, rt, nm);
-		return;
+		return(my_default(ob, il, nm));
 	}
 				/* set up sampling */
-	n = PI * il->sampdens;
-	nalt = sqrt(n/PI) + .5;
-	nazi = PI*nalt + .5;
-	n = nalt*nazi;
-	distarr = (float *)calloc(n, 3*sizeof(float));
-	if (distarr == NULL)
-		error(SYSTEM, "out of memory in o_face");
-	mkaxes(u, v, fa->norm);
+	if (il->sd != NULL) {
+		if (!getBSDF_xfm(xfm, fa->norm, il->udir)) {
+			objerror(ob, WARNING, "illegal up direction");
+			freeface(ob);
+			return(my_default(ob, il, nm));
+		}
+		n = il->sd->ninc;
+	} else {
+		if (il->sampdens <= 0) {
+			nalt = nazi = 1;	/* diffuse assumption */
+		} else {
+			n = PI * il->sampdens;
+			nalt = sqrt(n/PI) + .5;
+			nazi = PI*nalt + .5;
+		}
+		n = nazi*nalt;
+	}
+	newdist(n);
+				/* take first edge >= sqrt(area) */
+	for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
+		u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
+		u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
+		u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
+		if ((r1 = DOT(u,u)) >= fa->area-FTINY)
+			break;
+	}
+	if (i < fa->nv) {	/* got one! -- let's align our axes */
+		r2 = 1.0/sqrt(r1);
+		u[0] *= r2; u[1] *= r2; u[2] *= r2;
+		fcross(v, fa->norm, u);
+	} else			/* oh well, we'll just have to wing it */
+		mkaxes(u, v, fa->norm);
+				/* now, find limits in (u,v) coordinates */
 	ur[0] = vr[0] = FHUGE;
 	ur[1] = vr[1] = -FHUGE;
 	for (i = 0; i < fa->nv; i++) {
@@ -77,63 +230,76 @@ char  *nm;
 	dim[0] = random();
 				/* sample polygon */
 	nmisses = 0;
-	for (dim[1] = 0; dim[1] < nalt; dim[1]++)
-	    for (dim[2] = 0; dim[2] < nazi; dim[2]++)
+	for (dim[1] = 0; dim[1] < n; dim[1]++)
 		for (i = 0; i < il->nsamps; i++) {
-					/* random direction */
-		    dim[3] = 1;
-		    r1 = (dim[1]+urand(urind(ilhash(dim,4),i)))/nalt;
-		    dim[3] = 2;
-		    r2 = (dim[2]+urand(urind(ilhash(dim,4),i)))/nazi;
-		    flatdir(dn, r1, r2);
-		    for (j = 0; j < 3; j++)
-			dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*fa->norm[j];
-					/* random location */
+					/* randomize direction */
+		    h = ilhash(dim, 2) + i;
+		    if (il->sd != NULL) {
+			r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
+		    } else {
+			multisamp(sp, 2, urand(h));
+			alti = dim[1]/nazi;
+			r1 = (alti + sp[0])/nalt;
+			r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
+			flatdir(dn, r1, r2);
+			for (j = 0; j < 3; j++)
+			    dir[j] = -dn[0]*u[j] - dn[1]*v[j] -
+						dn[2]*fa->norm[j];
+		    }
+					/* randomize location */
 		    do {
-			dim[3] = 3;
-			r1 = ur[0] + (ur[1]-ur[0]) *
-					urand(urind(ilhash(dim,4),i+nmisses));
-			dim[3] = 4;
-			r2 = vr[0] + (vr[1]-vr[0]) *
-					urand(urind(ilhash(dim,4),i+nmisses));
+			multisamp(sp, 2, urand(h+4862+nmisses));
+			r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
+			r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
 			for (j = 0; j < 3; j++)
 			    org[j] = r1*u[j] + r2*v[j]
 					+ fa->offset*fa->norm[j];
 		    } while (!inface(org, fa) && nmisses++ < MAXMISS);
 		    if (nmisses > MAXMISS) {
 			objerror(ob, WARNING, "bad aspect");
-			rt->nrays = 0;
+			rayclean();
 			freeface(ob);
-			free((char *)distarr);
-			o_default(ob, il, rt, nm);
-			return;
+			return(my_default(ob, il, nm));
 		    }
+		    if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY)
+			r1 = -1.0001*il->thick - .0001;
+		    else
+			r1 = .0001;
 		    for (j = 0; j < 3; j++)
-			org[j] += .001*fa->norm[j];
+			org[j] += r1*fa->norm[j];
 					/* send sample */
-		    raysamp(distarr+dim[1]*nazi+dim[2], org, dir, rt);
+		    raysamp(dim[1], org, dir);
 		}
-	rayflush(rt);
-				/* write out the face w/ distribution */
-	flatout(il, distarr, nalt, nazi, u, v, fa->norm);
-	illumout(il, ob);
+	rayclean();
+	if (il->sd != NULL) {	/* run distribution through BSDF */
+		nalt = sqrt(il->sd->nout/PI) + .5;
+		nazi = PI*nalt + .5;
+		redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm);
+	}
+				/* write out the face and its distribution */
+	if (average(il, distarr, n)) {
+		if (il->sampdens > 0)
+			flatout(il, distarr, nalt, nazi, u, v, fa->norm);
+		illumout(il, ob);
+	} else
+		printobj(il->altmat, ob);
 				/* clean up */
 	freeface(ob);
-	free((char *)distarr);
+	return(0);
 #undef MAXMISS
 }
 
 
-o_sphere(ob, il, rt, nm)	/* make an illum sphere */
-register OBJREC  *ob;
-struct illum_args  *il;
-struct rtproc  *rt;
-char  *nm;
+int
+my_sphere(	/* make an illum sphere */
+	register OBJREC  *ob,
+	struct illum_args  *il,
+	char  *nm
+)
 {
-	int  dim[4];
+	int  dim[3];
 	int  n, nalt, nazi;
-	float  *distarr;
-	double  r1, r2, r3;
+	double  sp[4], r1, r2, r3;
 	FVECT  org, dir;
 	FVECT  u, v;
 	register int  i, j;
@@ -141,35 +307,32 @@ char  *nm;
 	if (ob->oargs.nfargs != 4)
 		objerror(ob, USER, "bad # of arguments");
 				/* set up sampling */
-	n = 4.*PI * il->sampdens;
-	nalt = sqrt(n/PI) + .5;
-	nazi = PI*nalt + .5;
+	if (il->sampdens <= 0)
+		nalt = nazi = 1;
+	else {
+		n = 4.*PI * il->sampdens;
+		nalt = sqrt(2./PI*n) + .5;
+		nazi = PI/2.*nalt + .5;
+	}
+	if (il->sd != NULL)
+		objerror(ob, WARNING, "BSDF ignored");
 	n = nalt*nazi;
-	distarr = (float *)calloc(n, 3*sizeof(float));
-	if (distarr == NULL)
-		error(SYSTEM, "out of memory in o_sphere");
+	newdist(n);
 	dim[0] = random();
 				/* sample sphere */
 	for (dim[1] = 0; dim[1] < nalt; dim[1]++)
-	    for (dim[2] = 0; dim[2] < nazi; dim[2]++) {
-		if (il->nsamps > 2 && nazi > 20) {
-		    rounddir(dir, (dim[1]+.5)/nalt, (dim[2]+.5)/nazi);
-		    mkaxes(u, v, dir);
-		}
+	    for (dim[2] = 0; dim[2] < nazi; dim[2]++)
 		for (i = 0; i < il->nsamps; i++) {
+					/* next sample point */
+		    multisamp(sp, 4, urand(ilhash(dim,3)+i));
 					/* random direction */
-		    dim[3] = 1;
-		    r1 = (dim[1]+urand(urind(ilhash(dim,4),i)))/nalt;
-		    dim[3] = 2;
-		    r2 = (dim[2]+urand(urind(ilhash(dim,4),i)))/nazi;
+		    r1 = (dim[1] + sp[0])/nalt;
+		    r2 = (dim[2] + sp[1] - .5)/nazi;
 		    rounddir(dir, r1, r2);
 					/* random location */
-		    if (il->nsamps <= 2 || nazi <= 20)
-			mkaxes(u, v, dir);		/* yuck! */
-		    dim[3] = 3;
-		    r3 = sqrt(urand(urind(ilhash(dim,4),i)));
-		    dim[3] = 4;
-		    r2 = 2.*PI*urand(urind(ilhash(dim,4),i));
+		    mkaxes(u, v, dir);		/* yuck! */
+		    r3 = sqrt(sp[2]);
+		    r2 = 2.*PI*sp[3];
 		    r1 = r3*ob->oargs.farg[3]*cos(r2);
 		    r2 = r3*ob->oargs.farg[3]*sin(r2);
 		    r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
@@ -179,161 +342,110 @@ char  *nm;
 			dir[j] = -dir[j];
 		    }
 					/* send sample */
-		    raysamp(distarr+dim[1]*nazi+dim[2], org, dir, rt);
+		    raysamp(dim[1]*nazi+dim[2], org, dir);
 		}
-	    }
-	rayflush(rt);
-				/* write out the sphere w/ distribution */
-	roundout(il, distarr, nalt, nazi);
-	illumout(il, ob);
+	rayclean();
+				/* write out the sphere and its distribution */
+	if (average(il, distarr, n)) {
+		if (il->sampdens > 0)
+			roundout(il, distarr, nalt, nazi);
+		else
+			objerror(ob, WARNING, "diffuse distribution");
+		illumout(il, ob);
+	} else
+		printobj(il->altmat, ob);
 				/* clean up */
-	free((char *)distarr);
+	return(1);
 }
 
 
-o_ring(ob, il, rt, nm)		/* make an illum ring */
-OBJREC  *ob;
-struct illum_args  *il;
-struct rtproc  *rt;
-char  *nm;
+int
+my_ring(		/* make an illum ring */
+	OBJREC  *ob,
+	struct illum_args  *il,
+	char  *nm
+)
 {
-	int  dim[4];
-	int  n, nalt, nazi;
-	float  *distarr;
-	double  r1, r2, r3;
+	int  dim[2];
+	int  n, nalt, nazi, alti;
+	double  sp[2], r1, r2, r3;
+	int  h;
 	FVECT  dn, org, dir;
 	FVECT  u, v;
-	register CONE  *co;
+	MAT4  xfm;
+	CONE  *co;
 	register int  i, j;
 				/* get/check arguments */
 	co = getcone(ob, 0);
 				/* set up sampling */
-	n = PI * il->sampdens;
-	nalt = sqrt(n/PI) + .5;
-	nazi = PI*nalt + .5;
-	n = nalt*nazi;
-	distarr = (float *)calloc(n, 3*sizeof(float));
-	if (distarr == NULL)
-		error(SYSTEM, "out of memory in o_ring");
+	if (il->sd != NULL) {
+		if (!getBSDF_xfm(xfm, co->ad, il->udir)) {
+			objerror(ob, WARNING, "illegal up direction");
+			freecone(ob);
+			return(my_default(ob, il, nm));
+		}
+		n = il->sd->ninc;
+	} else {
+		if (il->sampdens <= 0) {
+			nalt = nazi = 1;	/* diffuse assumption */
+		} else {
+			n = PI * il->sampdens;
+			nalt = sqrt(n/PI) + .5;
+			nazi = PI*nalt + .5;
+		}
+		n = nazi*nalt;
+	}
+	newdist(n);
 	mkaxes(u, v, co->ad);
 	dim[0] = random();
 				/* sample disk */
-	for (dim[1] = 0; dim[1] < nalt; dim[1]++)
-	    for (dim[2] = 0; dim[2] < nazi; dim[2]++)
+	for (dim[1] = 0; dim[1] < n; dim[1]++)
 		for (i = 0; i < il->nsamps; i++) {
-					/* random direction */
-		    dim[3] = 1;
-		    r1 = (dim[1]+urand(urind(ilhash(dim,4),i)))/nalt;
-		    dim[3] = 2;
-		    r2 = (dim[2]+urand(urind(ilhash(dim,4),i)))/nalt;
-		    flatdir(dn, r1, r2);
-		    for (j = 0; j < 3; j++)
+					/* next sample point */
+		    h = ilhash(dim,2) + i;
+					/* randomize direction */
+		    if (il->sd != NULL) {
+			r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
+		    } else {
+			multisamp(sp, 2, urand(h));
+			alti = dim[1]/nazi;
+			r1 = (alti + sp[0])/nalt;
+			r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
+			flatdir(dn, r1, r2);
+			for (j = 0; j < 3; j++)
 			dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
-					/* random location */
-		    dim[3] = 3;
+		    }
+					/* randomize location */
+		    multisamp(sp, 2, urand(h+8371));
 		    r3 = sqrt(CO_R0(co)*CO_R0(co) +
-				urand(urind(ilhash(dim,4),i))*
-				(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
-		    dim[3] = 4;
-		    r2 = 2.*PI*urand(urind(ilhash(dim,4),i));
+			    sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
+		    r2 = 2.*PI*sp[1];
 		    r1 = r3*cos(r2);
 		    r2 = r3*sin(r2);
+		    if (il->sd != NULL && DOT(dir, co->ad) < -FTINY)
+			r3 = -1.0001*il->thick - .0001;
+		    else
+			r3 = .0001;
 		    for (j = 0; j < 3; j++)
-			org[j] = CO_P0(co)[j] + r1*u[j] + r1*v[j] +
-					.001*co->ad[j];
-
+			org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
+						r3*co->ad[j];
 					/* send sample */
-		    raysamp(distarr+dim[1]*nazi+dim[2], org, dir, rt);
+		    raysamp(dim[1], org, dir);
 		}
-	rayflush(rt);
-				/* write out the ring w/ distribution */
-	flatout(il, distarr, nalt, nazi, u, v, co->ad);
-	illumout(il, ob);
+	rayclean();
+	if (il->sd != NULL) {	/* run distribution through BSDF */
+		nalt = sqrt(il->sd->nout/PI) + .5;
+		nazi = PI*nalt + .5;
+		redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
+	}
+				/* write out the ring and its distribution */
+	if (average(il, distarr, n)) {
+		if (il->sampdens > 0)
+			flatout(il, distarr, nalt, nazi, u, v, co->ad);
+		illumout(il, ob);
+	} else
+		printobj(il->altmat, ob);
 				/* clean up */
 	freecone(ob);
-	free((char *)distarr);
-}
-
-
-raysamp(res, org, dir, rt)	/* compute a ray sample */
-float  res[3];
-FVECT  org, dir;
-register struct rtproc  *rt;
-{
-	register float  *fp;
-
-	if (rt->nrays == rt->bsiz)
-		rayflush(rt);
-	rt->dest[rt->nrays] = res;
-	fp = rt->buf + 6*rt->nrays++;
-	*fp++ = org[0]; *fp++ = org[1]; *fp++ = org[2];
-	*fp++ = dir[0]; *fp++ = dir[1]; *fp = dir[2];
-}
-
-
-rayflush(rt)			/* flush buffered rays */
-register struct rtproc  *rt;
-{
-	register int  i;
-
-	if (rt->nrays <= 0)
-		return;
-	i = 6*rt->nrays + 3;
-	rt->buf[i++] = 0.; rt->buf[i++] = 0.; rt->buf[i] = 0.;
-	if ( process(rt->pd, (char *)rt->buf, (char *)rt->buf,
-			3*sizeof(float)*rt->nrays,
-			6*sizeof(float)*(rt->nrays+1)) <
-			3*sizeof(float)*rt->nrays )
-		error(SYSTEM, "error reading from rtrace process");
-	i = rt->nrays;
-	while (i--) {
-		rt->dest[i][0] += rt->buf[3*i];
-		rt->dest[i][1] += rt->buf[3*i+1];
-		rt->dest[i][2] += rt->buf[3*i+2];
-	}
-	rt->nrays = 0;
-}
-
-
-mkaxes(u, v, n)			/* compute u and v to go with n */
-FVECT  u, v, n;
-{
-	register int  i;
-
-	v[0] = v[1] = v[2] = 0.0;
-	for (i = 0; i < 3; i++)
-		if (n[i] < 0.6 && n[i] > -0.6)
-			break;
-	v[i] = 1.0;
-	fcross(u, v, n);
-	normalize(u);
-	fcross(v, n, u);
-}
-
-
-rounddir(dv, alt, azi)		/* compute uniform spherical direction */
-register FVECT  dv;
-double  alt, azi;
-{
-	double  d1, d2;
-
-	dv[2] = 1. - 2.*alt;
-	d1 = sqrt(1. - dv[2]*dv[2]);
-	d2 = 2.*PI * azi;
-	dv[0] = d1*cos(d2);
-	dv[1] = d1*sin(d2);
-}
-
-
-flatdir(dv, alt, azi)		/* compute uniform hemispherical direction */
-register FVECT  dv;
-double  alt, azi;
-{
-	double  d1, d2;
-
-	d1 = sqrt(alt);
-	d2 = 2.*PI * azi;
-	dv[0] = d1*cos(d2);
-	dv[1] = d1*sin(d2);
-	dv[2] = sqrt(1. - d1*d1);
+	return(1);
 }