--- ray/src/gen/mkillum2.c	1991/07/23 15:42:42	1.1
+++ ray/src/gen/mkillum2.c	2011/08/16 18:09:53	2.37
@@ -1,50 +1,820 @@
-/* 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.37 2011/08/16 18:09:53 greg Exp $";
 #endif
-
 /*
- * Routines to do the actual calcultion and output for mkillum
+ * Routines to do the actual calculation for mkillum
  */
 
-#include  "mkillum.h"
+#include <string.h>
 
+#include  "mkillum.h"
 #include  "face.h"
-
 #include  "cone.h"
+#include  "source.h"
+#include  "paths.h"
 
+#ifndef NBSDFSAMPS
+#define NBSDFSAMPS	256		/* BSDF resampling count */
+#endif
 
-printobj(mod, obj)		/* print out an object */
-char  *mod;
-register OBJREC  *obj;
+COLORV *	distarr = NULL;		/* distribution array */
+int		distsiz = 0;
+COLORV *	direct_discount = NULL;	/* amount to take off direct */
+
+
+void
+newdist(			/* allocate & clear distribution array */
+	int siz
+)
 {
+	if (siz <= 0) {
+		if (distsiz > 0)
+			free(distarr);
+		distarr = NULL;
+		distsiz = 0;
+		return;
+	}
+	if (distsiz < siz) {
+		if (distsiz > 0)
+			free(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);
+}
+
+
+static void
+new_discount()			/* allocate space for direct contrib. record */
+{
+	if (distsiz <= 0)
+		return;
+	direct_discount = (COLORV *)calloc(distsiz, sizeof(COLOR));
+	if (direct_discount == NULL)
+		error(SYSTEM, "out of memory in new_discount");
+}
+
+
+static void
+done_discount()			/* clear off direct contrib. record */
+{
+	if (direct_discount == NULL)
+		return;
+	free(direct_discount);
+	direct_discount = NULL;
+}
+
+
+int
+process_ray(			/* process a ray result or report error */
+	RAY *r,
+	int rv
+)
+{
+	COLORV	*colp;
+
+	if (rv == 0)			/* no result ready */
+		return(0);
+	if (rv < 0)
+		error(USER, "ray tracing process died");
+	if (r->rno >= distsiz)
+		error(INTERNAL, "bad returned index in process_ray");
+	multcolor(r->rcol, r->rcoef);	/* in case it's a source ray */
+	colp = &distarr[r->rno * 3];
+	addcolor(colp, r->rcol);
+	if (r->rsrc >= 0 &&		/* remember source contrib. */
+			direct_discount != NULL) {
+		colp = &direct_discount[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|SPECULAR, NULL, NULL);
+	myRay.rno = ndx;
+					/* queue ray, check result */
+	process_ray(&myRay, ray_pqueue(&myRay));
+}
+
+
+void
+srcsamps(			/* sample sources from this surface position */
+	struct illum_args *il,
+	FVECT org,
+	FVECT nrm,
+	MAT4 ixfm
+)
+{
+	int  nalt=1, nazi=1;
+	SRCINDEX  si;
+	RAY  sr;
+	FVECT	v;
+	double  d;
+	int  i, j;
+						/* get sampling density */
+	if (il->sd == NULL && il->sampdens > 0) {
+		i = PI * il->sampdens;
+		nalt = sqrt(i/PI) + .5;
+		nazi = PI*nalt + .5;
+	}
+	initsrcindex(&si);			/* loop over (sub)sources */
+	for ( ; ; ) {
+		VCOPY(sr.rorg, org);		/* pick side to shoot from */
+		if (il->sd != NULL) {
+			int  sn = si.sn;
+			if (si.sp+1 >= si.np) ++sn;
+			if (sn >= nsources) break;
+			if (source[sn].sflags & SDISTANT)
+				d = DOT(source[sn].sloc, nrm);
+			else {
+				VSUB(v, source[sn].sloc, org);
+				d = DOT(v, nrm);
+			}
+		} else
+			d = 1.0;		/* only transmission */
+		if (d < 0.0)
+			d = -1.0001*il->thick - 5.*FTINY;
+		else
+			d = 5.*FTINY;
+		VSUM(sr.rorg, sr.rorg, nrm, d);
+		samplendx++;			/* increment sample counter */
+		if (!srcray(&sr, NULL, &si))
+			break;			/* end of sources */
+						/* index direction */
+		if (ixfm != NULL)
+			multv3(v, sr.rdir, ixfm);
+		else
+			VCOPY(v, sr.rdir);
+		if (il->sd != NULL) {
+			i = getBSDF_incndx(il->sd, v);
+			if (i < 0)
+				continue;	/* must not be important */
+			sr.rno = i;
+			d = 1.0/getBSDF_incohm(il->sd, i);
+		} else {
+			if (v[2] >= -FTINY)
+				continue;	/* only sample transmission */
+			v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2];
+			sr.rno = flatindex(v, nalt, nazi);
+			d = nalt*nazi*(1./PI) * v[2];
+		}
+		d *= si.dom;			/* solid angle correction */
+		scalecolor(sr.rcoef, d);
+		process_ray(&sr, ray_pqueue(&sr));
+	}
+}
+
+
+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;
 
-	printf("\n%s %s %s", mod, ofun[obj->otype].funame, obj->oname);
-	printf("\n%d", obj->oargs.nsargs);
-	for (i = 0; i < obj->oargs.nsargs; i++)
-		printf(" %s", obj->oargs.sarg[i]);
-#ifdef  IARGS
-	printf("\n%d", obj->oargs.niargs);
-	for (i = 0; i < obj->oargs.niargs; i++)
-		printf(" %d", obj->oargs.iarg[i]);
-#else
-	printf("\n0");
-#endif
-	printf("\n%d", obj->oargs.nfargs);
-	for (i = 0; i < obj->oargs.nfargs; i++) {
-		if (i%3 == 0)
-			putchar('\n');
-		printf(" %18.12g", obj->oargs.farg[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);
+}
+
+
+void
+flatdir(		/* compute uniform hemispherical direction */
+	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
+flatindex(		/* compute index for hemispherical direction */
+	FVECT	dv,
+	int	nalt,
+	int	nazi
+)
+{
+	double	d;
+	int	i, j;
+	
+	d = 1.0 - dv[2]*dv[2];
+	i = d*nalt;
+	d = atan2(dv[1], dv[0]) * (0.5/PI);
+	if (d < 0.0) d += 1.0;
+	j = d*nazi + 0.5;
+	if (j >= nazi) j = 0;
+	return(i*nazi + j);
+}
+
+
+int
+printgeom(		/* print out detailed geometry for BSDF */
+	struct BSDF_data *sd,
+	char *xfrot,
+	FVECT ctr,
+	double s1,
+	double s2
+)
+{
+        static char	mgftemp[] = TEMPLATE;
+	char		cmdbuf[64];
+	FILE		*fp;
+	double		sca;
+
+	if (sd == NULL || sd->mgf == NULL)
+		return(0);
+	if (sd->dim[0] <= FTINY || sd->dim[1] <= FTINY)
+		return(0);
+	if ((s1 > s2) ^ (sd->dim[0] > sd->dim[1])) {
+		sca = s1; s1 = s2; s2 = sca;
 	}
-	putchar('\n');
+	s1 /= sd->dim[0];
+	s2 /= sd->dim[1];
+	sca = s1 > s2 ? s1 : s2;
+	strcpy(mgftemp, TEMPLATE);
+	if ((fp = fopen(mktemp(mgftemp), "w")) == NULL)
+		error(SYSTEM, "cannot create temporary file for MGF");
+					/* prepend our transform */
+	fprintf(fp, "xf%s -s %.5f -t %.5g %.5g %.5g\n",
+			xfrot, sca, ctr[0], ctr[1], ctr[2]);
+					/* output given MGF description */
+	fputs(sd->mgf, fp);
+	fputs("\nxf\n", fp);
+	if (fclose(fp) == EOF)
+		error(SYSTEM, "error writing MGF temporary file");
+					/* execute mgf2rad to convert MGF */
+	strcpy(cmdbuf, "mgf2rad ");
+	strcpy(cmdbuf+8, mgftemp);
+	fflush(stdout);
+	system(cmdbuf);
+	unlink(mgftemp);		/* clean up */
+	return(1);
 }
 
 
-mkillum(ob, il, rt)		/* make an illum object */
-OBJREC  *ob;
-struct illum_args  *il;
-struct rtproc  *rt;
+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);
+	printobj(il->altmat, ob);
+	return(1);
+}
+
+
+int
+my_face(		/* make an illum face */
+	OBJREC  *ob,
+	struct illum_args  *il,
+	char  *nm
+)
+{
+	int  dim[2];
+	int  n, nalt, nazi, alti;
+	double  sp[2], r1, r2;
+	int  h;
+	FVECT  dn, org, dir;
+	FVECT  u, v;
+	double  ur[2], vr[2];
+	MAT4  xfm;
+	char  xfrot[64];
+	int  nallow;
+	FACE  *fa;
+	int  i, j;
+				/* get/check arguments */
+	fa = getface(ob);
+	if (fa->area == 0.0) {
+		freeface(ob);
+		return(my_default(ob, il, nm));
+	}
+				/* set up sampling */
+	if (il->sd != NULL) {
+		if (!getBSDF_xfm(xfm, fa->norm, il->udir, xfrot)) {
+			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++) {
+		r1 = DOT(VERTEX(fa,i),u);
+		if (r1 < ur[0]) ur[0] = r1;
+		if (r1 > ur[1]) ur[1] = r1;
+		r2 = DOT(VERTEX(fa,i),v);
+		if (r2 < vr[0]) vr[0] = r2;
+		if (r2 > vr[1]) vr[1] = r2;
+	}
+				/* output detailed geometry? */
+	if (!(il->flags & IL_LIGHT) && il->sd != NULL && il->sd->mgf != NULL &&
+			il->thick <= FTINY) {
+		for (j = 3; j--; )
+			org[j] = .5*(ur[0]+ur[1])*u[j] +
+					.5*(vr[0]+vr[1])*v[j] +
+					fa->offset*fa->norm[j];
+		printgeom(il->sd, xfrot, org, ur[1]-ur[0], vr[1]-vr[0]);
+	}
+	dim[0] = random();
+				/* sample polygon */
+	nallow = 5*n*il->nsamps;
+	for (dim[1] = 0; dim[1] < n; dim[1]++)
+		for (i = 0; i < il->nsamps; i++) {
+					/* 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 {
+			multisamp(sp, 2, urand(h+4862+nallow));
+			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) && nallow-- > 0);
+		    if (nallow < 0) {
+			objerror(ob, WARNING, "bad aspect");
+			rayclean();
+			freeface(ob);
+			return(my_default(ob, il, nm));
+		    }
+		    if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY)
+			r1 = -1.0001*il->thick - 5.*FTINY;
+		    else
+			r1 = 5.*FTINY;
+		    for (j = 0; j < 3; j++)
+			org[j] += r1*fa->norm[j];
+					/* send sample */
+		    raysamp(dim[1], org, dir);
+		}
+				/* add in direct component? */
+	if (il->flags & IL_LIGHT || il->sd != NULL) {
+		MAT4	ixfm;
+		if (il->sd == NULL) {
+			for (i = 3; i--; ) {
+				ixfm[i][0] = u[i];
+				ixfm[i][1] = v[i];
+				ixfm[i][2] = fa->norm[i];
+				ixfm[i][3] = 0.;
+			}
+			ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
+			ixfm[3][3] = 1.;
+		} else {
+			if (!invmat4(ixfm, xfm))
+				objerror(ob, INTERNAL,
+					"cannot invert BSDF transform");
+			if (!(il->flags & IL_LIGHT))
+				new_discount();
+		}
+		dim[0] = random();
+		nallow = 10*il->nsamps;
+		for (i = 0; i < il->nsamps; i++) {
+					/* randomize location */
+		    h = dim[0] + samplendx++;
+		    do {
+			multisamp(sp, 2, urand(h+nallow));
+			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) && nallow-- > 0);
+		    if (nallow < 0) {
+			objerror(ob, WARNING, "bad aspect");
+			rayclean();
+			freeface(ob);
+			return(my_default(ob, il, nm));
+		    }
+					/* sample source rays */
+		    srcsamps(il, org, fa->norm, ixfm);
+		}
+	}
+				/* wait for all rays to finish */
+	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);
+		done_discount();
+		if (!il->sampdens)
+			il->sampdens = nalt*nazi/PI + .999;
+	}
+				/* 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);
+	return(0);
+}
+
+
+int
+my_sphere(	/* make an illum sphere */
+	register OBJREC  *ob,
+	struct illum_args  *il,
+	char  *nm
+)
+{
+	int  dim[3];
+	int  n, nalt, nazi;
+	double  sp[4], r1, r2, r3;
+	FVECT  org, dir;
+	FVECT  u, v;
+	register int  i, j;
+				/* check arguments */
+	if (ob->oargs.nfargs != 4)
+		objerror(ob, USER, "bad # of arguments");
+				/* set up sampling */
+	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;
+	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]++)
+		for (i = 0; i < il->nsamps; i++) {
+					/* next sample point */
+		    multisamp(sp, 4, urand(ilhash(dim,3)+i));
+					/* random direction */
+		    r1 = (dim[1] + sp[0])/nalt;
+		    r2 = (dim[2] + sp[1] - .5)/nazi;
+		    rounddir(dir, r1, r2);
+					/* random location */
+		    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);
+		    for (j = 0; j < 3; j++) {
+			org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
+					r3*dir[j];
+			dir[j] = -dir[j];
+		    }
+					/* send sample */
+		    raysamp(dim[1]*nazi+dim[2], org, dir);
+		}
+				/* wait for all rays to finish */
+	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 */
+	return(1);
+}
+
+
+int
+my_ring(		/* make an illum ring */
+	OBJREC  *ob,
+	struct illum_args  *il,
+	char  *nm
+)
+{
+	int  dim[2];
+	int  n, nalt, nazi, alti;
+	double  sp[2], r1, r2, r3;
+	int  h;
+	FVECT  dn, org, dir;
+	FVECT  u, v;
+	MAT4  xfm;
+	CONE  *co;
+	int  i, j;
+				/* get/check arguments */
+	co = getcone(ob, 0);
+				/* set up sampling */
+	if (il->sd != NULL) {
+		if (!getBSDF_xfm(xfm, co->ad, il->udir, NULL)) {
+			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] < n; dim[1]++)
+		for (i = 0; i < il->nsamps; i++) {
+					/* 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];
+		    }
+					/* randomize location */
+		    multisamp(sp, 2, urand(h+8371));
+		    r3 = sqrt(CO_R0(co)*CO_R0(co) +
+			    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 - 5.*FTINY;
+		    else
+			r3 = 5.*FTINY;
+		    for (j = 0; j < 3; j++)
+			org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
+						r3*co->ad[j];
+					/* send sample */
+		    raysamp(dim[1], org, dir);
+		}
+				/* add in direct component? */
+	if (il->flags & IL_LIGHT || il->sd != NULL) {
+		MAT4	ixfm;
+		if (il->sd == NULL) {
+			for (i = 3; i--; ) {
+				ixfm[i][0] = u[i];
+				ixfm[i][1] = v[i];
+				ixfm[i][2] = co->ad[i];
+				ixfm[i][3] = 0.;
+			}
+			ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
+			ixfm[3][3] = 1.;
+		} else {
+			if (!invmat4(ixfm, xfm))
+				objerror(ob, INTERNAL,
+					"cannot invert BSDF transform");
+			if (!(il->flags & IL_LIGHT))
+				new_discount();
+		}
+		dim[0] = random();
+		for (i = 0; i < il->nsamps; i++) {
+					/* randomize location */
+		    h = dim[0] + samplendx++;
+		    multisamp(sp, 2, urand(h));
+		    r3 = sqrt(CO_R0(co)*CO_R0(co) +
+			    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);
+		    for (j = 0; j < 3; j++)
+			org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j];
+					/* sample source rays */
+		    srcsamps(il, org, co->ad, ixfm);
+		}
+	}
+				/* wait for all rays to finish */
+	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);
+		done_discount();
+		if (!il->sampdens)
+			il->sampdens = nalt*nazi/PI + .999;
+	}
+				/* 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);
+	return(1);
+}
+
+
+void
+redistribute(		/* pass distarr ray sums through BSDF */
+	struct BSDF_data *b,
+	int nalt,
+	int nazi,
+	FVECT u,
+	FVECT v,
+	FVECT w,
+	MAT4 xm
+)
+{
+	int	nout = 0;
+	MAT4	mymat, inmat;
+	COLORV	*idist;
+	COLORV	*cp;
+	FVECT	dv;
+	double	wt;
+	int	i, j, k, c, o;
+	COLOR	col, cinc;
+					/* copy incoming distribution */
+	if (b->ninc > distsiz)
+		error(INTERNAL, "error 1 in redistribute");
+	idist = (COLORV *)malloc(sizeof(COLOR)*b->ninc);
+	if (idist == NULL)
+		error(SYSTEM, "out of memory in redistribute");
+	memcpy(idist, distarr, sizeof(COLOR)*b->ninc);
+					/* compose direction transform */
+	for (i = 3; i--; ) {
+		mymat[i][0] = u[i];
+		mymat[i][1] = v[i];
+		mymat[i][2] = w[i];
+		mymat[i][3] = 0.;
+	}
+	mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
+	mymat[3][3] = 1.;
+	if (xm != NULL)
+		multmat4(mymat, xm, mymat);
+	for (i = 3; i--; ) {		/* make sure it's normalized */
+		wt = 1./sqrt(	mymat[0][i]*mymat[0][i] +
+				mymat[1][i]*mymat[1][i] +
+				mymat[2][i]*mymat[2][i]	);
+		for (j = 3; j--; )
+			mymat[j][i] *= wt;
+	}
+	if (!invmat4(inmat, mymat))	/* need inverse as well */
+		error(INTERNAL, "cannot invert BSDF transform");
+	newdist(nalt*nazi);		/* resample distribution */
+	for (i = b->ninc; i--; ) {
+		int	direct_out = -1;
+		COLOR	cdir;
+		getBSDF_incvec(dv, b, i);	/* compute incident irrad. */
+		multv3(dv, dv, mymat);
+		if (dv[2] < 0.0) {
+			dv[0] = -dv[0]; dv[1] = -dv[1]; dv[2] = -dv[2];
+			direct_out += (direct_discount != NULL);
+		}
+		wt = getBSDF_incohm(b, i);
+		wt *= dv[2];			/* solid_angle*cosine(theta) */
+		cp = &idist[3*i];
+		copycolor(cinc, cp);
+		scalecolor(cinc, wt);
+		if (!direct_out) {		/* discount direct contr. */
+			cp = &direct_discount[3*i];
+			copycolor(cdir, cp);
+			scalecolor(cdir, -wt);
+			if (b->nout != b->ninc)
+				direct_out = flatindex(dv, nalt, nazi);
+			else
+				direct_out = i;	/* assumes dist. mirroring */
+		}
+		for (k = nalt; k--; )		/* loop over distribution */
+		  for (j = nazi; j--; ) {
+		    int	rstart = random();
+		    for (c = NBSDFSAMPS; c--; ) {
+			double  sp[2];
+			multisamp(sp, 2, urand(rstart+c));
+			flatdir(dv, (k + sp[0])/nalt,
+					(j + .5 - sp[1])/nazi);
+			multv3(dv, dv, inmat);
+						/* evaluate BSDF @ outgoing */
+			o = getBSDF_outndx(b, dv);
+			if (o < 0) {
+				nout++;
+				continue;
+			}
+			wt = BSDF_value(b, i, o) * (1./NBSDFSAMPS);
+			copycolor(col, cinc);
+			if (b->nout != b->ninc)
+				o = k*nazi + j;
+			if (o == direct_out)
+				addcolor(col, cdir);	/* minus direct */
+			scalecolor(col, wt);
+			cp = &distarr[3*(k*nazi + j)];
+			addcolor(cp, col);	/* sum into distribution */
+		    }
+		  }
+	}
+	free(idist);			/* free temp space */
+	if (nout) {
+		sprintf(errmsg, "missing %.1f%% of BSDF directions",
+				100.*nout/(b->ninc*nalt*nazi*NBSDFSAMPS));
+		error(WARNING, errmsg);
+	}
 }