--- ray/src/gen/mkillum2.c	2007/09/18 19:51:07	2.20
+++ ray/src/gen/mkillum2.c	2007/09/21 05:53:21	2.21
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: mkillum2.c,v 2.20 2007/09/18 19:51:07 greg Exp $";
+static const char	RCSid[] = "$Id: mkillum2.c,v 2.21 2007/09/21 05:53:21 greg Exp $";
 #endif
 /*
  * Routines to do the actual calculation for mkillum
@@ -13,16 +13,11 @@ static const char	RCSid[] = "$Id: mkillum2.c,v 2.20 20
 #include  "random.h"
 
 
-static void mkaxes(FVECT u, FVECT v, FVECT n);
-static void rounddir(FVECT dv, double alt, double azi);
-static void flatdir(FVECT dv, double alt, double azi);
+COLORV *	distarr = NULL;		/* distribution array */
+int		distsiz = 0;
 
 
-static COLORV *	distarr = NULL;		/* distribution array */
-static int	distsiz = 0;
-
-
-static void
+void
 newdist(			/* allocate & clear distribution array */
 	int siz
 )
@@ -35,17 +30,18 @@ newdist(			/* allocate & clear distribution array */
 		return;
 	}
 	if (distsiz < siz) {
-		free((void *)distarr);
-		distarr = (COLORV *)malloc(sizeof(COLORV)*3*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(COLORV)*3*siz);
+	memset(distarr, '\0', sizeof(COLOR)*siz);
 }
 
 
-static int
+int
 process_ray(RAY *r, int rv)
 {
 	COLORV	*colp;
@@ -62,7 +58,7 @@ process_ray(RAY *r, int rv)
 }
 
 
-static void
+void
 raysamp(	/* queue a ray sample */
 	int  ndx,
 	FVECT  org,
@@ -84,7 +80,7 @@ raysamp(	/* queue a ray sample */
 }
 
 
-static void
+void
 rayclean()			/* finish all pending rays */
 {
 	RAY	myRay;
@@ -94,6 +90,60 @@ rayclean()			/* finish all pending rays */
 }
 
 
+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,
@@ -117,14 +167,15 @@ my_face(		/* make an illum face */
 )
 {
 #define MAXMISS		(5*n*il->nsamps)
-	int  dim[3];
-	int  n, nalt, nazi, h;
+	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);
@@ -133,14 +184,22 @@ my_face(		/* make an illum face */
 		return(my_default(ob, il, nm));
 	}
 				/* set up sampling */
-	if (il->sampdens <= 0)
-		nalt = nazi = 1;
-	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 = nalt*nazi;
+	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
+		n = nazi*nalt;
 	newdist(n);
 				/* take first edge >= sqrt(area) */
 	for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
@@ -170,17 +229,22 @@ my_face(		/* make an illum face */
 	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 */
-		    h = ilhash(dim, 3) + i;
-		    multisamp(sp, 2, urand(h));
-		    r1 = (dim[1] + sp[0])/nalt;
-		    r2 = (dim[2] + 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];
+		    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];
+		    }
 					/* random location */
 		    do {
 			multisamp(sp, 2, urand(h+4862+nmisses));
@@ -194,17 +258,22 @@ my_face(		/* make an illum face */
 			objerror(ob, WARNING, "bad aspect");
 			rayclean();
 			freeface(ob);
-			free((void *)distarr);
 			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] += .0001*fa->norm[j];
+			org[j] += r1*fa->norm[j];
 					/* send sample */
-		    raysamp(dim[1]*nazi+dim[2], org, dir);
+		    raysamp(dim[1], org, dir);
 		}
 	rayclean();
+	if (il->sd != NULL)	/* run distribution through BSDF */
+		redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm);
 				/* write out the face and its distribution */
-	if (average(il, distarr, nalt*nazi)) {
+	if (average(il, distarr, n)) {
 		if (il->sampdens > 0)
 			flatout(il, distarr, nalt, nazi, u, v, fa->norm);
 		illumout(il, ob);
@@ -241,6 +310,8 @@ my_sphere(	/* make an illum sphere */
 		nalt = sqrt(2./PI*n) + .5;
 		nazi = PI/2.*nalt + .5;
 	}
+	if (il->sd != NULL)
+		objerror(ob, WARNING, "BSDF ignored");
 	n = nalt*nazi;
 	newdist(n);
 	dim[0] = random();
@@ -271,7 +342,7 @@ my_sphere(	/* make an illum sphere */
 		}
 	rayclean();
 				/* write out the sphere and its distribution */
-	if (average(il, distarr, nalt*nazi)) {
+	if (average(il, distarr, n)) {
 		if (il->sampdens > 0)
 			roundout(il, distarr, nalt, nazi);
 		else
@@ -291,12 +362,14 @@ my_ring(		/* make an illum ring */
 	char  *nm
 )
 {
-	int  dim[3];
-	int  n, nalt, nazi;
-	double  sp[4], 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);
@@ -308,37 +381,57 @@ my_ring(		/* make an illum ring */
 		nalt = sqrt(n/PI) + .5;
 		nazi = PI*nalt + .5;
 	}
-	n = nalt*nazi;
+	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
+		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++) {
 					/* next sample point */
-		    multisamp(sp, 4, urand(ilhash(dim,3)+i));
+		    h = ilhash(dim,2) + i;
 					/* random direction */
-		    r1 = (dim[1] + sp[0])/nalt;
-		    r2 = (dim[2] + sp[1] - .5)/nazi;
-		    flatdir(dn, r1, r2);
-		    for (j = 0; j < 3; j++)
+		    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 */
+		    multisamp(sp, 2, urand(h+8371));
 		    r3 = sqrt(CO_R0(co)*CO_R0(co) +
-			    sp[2]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
-		    r2 = 2.*PI*sp[3];
+			    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] + r2*v[j] +
-						.0001*co->ad[j];
+						r3*co->ad[j];
 					/* send sample */
-		    raysamp(dim[1]*nazi+dim[2], org, dir);
+		    raysamp(dim[1], org, dir);
 		}
 	rayclean();
+	if (il->sd != NULL)	/* run distribution through BSDF */
+		redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
 				/* write out the ring and its distribution */
-	if (average(il, distarr, nalt*nazi)) {
+	if (average(il, distarr, n)) {
 		if (il->sampdens > 0)
 			flatout(il, distarr, nalt, nazi, u, v, co->ad);
 		illumout(il, ob);
@@ -347,58 +440,4 @@ my_ring(		/* make an illum ring */
 				/* clean up */
 	freecone(ob);
 	return(1);
-}
-
-
-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);
 }