ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/gen/mkillum2.c

Comparing ray/src/gen/mkillum2.c (file contents):
Revision 2.22 by greg, Sun Sep 30 19:05:26 2007 UTC vs.
Revision 2.37 by greg, Tue Aug 16 18:09:53 2011 UTC

#	Line 10 | Line 10 | static const char      RCSid[] = "$Id$";
10		#include  "mkillum.h"
11		#include  "face.h"
12		#include  "cone.h"
13	<	#include  "random.h"
13	>	#include  "source.h"
14	>	#include  "paths.h"
15
16	+	#ifndef NBSDFSAMPS
17	+	#define NBSDFSAMPS      256             /* BSDF resampling count */
18	+	#endif
19
20		COLORV *        distarr = NULL;         /* distribution array */
21		int             distsiz = 0;
22	+	COLORV *        direct_discount = NULL; /* amount to take off direct */
23
24
25		void
#	Line 22 | Line 27 | newdist(                       /* allocate & clear distribution array */
27		int siz
28		)
29		{
30	<	if (siz == 0) {
30	>	if (siz <= 0) {
31		if (distsiz > 0)
32	<	free((void *)distarr);
32	>	free(distarr);
33		distarr = NULL;
34		distsiz = 0;
35		return;
36		}
37		if (distsiz < siz) {
38		if (distsiz > 0)
39	<	free((void *)distarr);
39	>	free(distarr);
40		distarr = (COLORV *)malloc(sizeof(COLOR)*siz);
41		if (distarr == NULL)
42		error(SYSTEM, "out of memory in newdist");
#	Line 41 | Line 46 | newdist(                       /* allocate & clear distribution array */
46		}
47
48
49	+	static void
50	+	new_discount()                  /* allocate space for direct contrib. record */
51	+	{
52	+	if (distsiz <= 0)
53	+	return;
54	+	direct_discount = (COLORV *)calloc(distsiz, sizeof(COLOR));
55	+	if (direct_discount == NULL)
56	+	error(SYSTEM, "out of memory in new_discount");
57	+	}
58	+
59	+
60	+	static void
61	+	done_discount()                 /* clear off direct contrib. record */
62	+	{
63	+	if (direct_discount == NULL)
64	+	return;
65	+	free(direct_discount);
66	+	direct_discount = NULL;
67	+	}
68	+
69	+
70		int
71	<	process_ray(RAY *r, int rv)
71	>	process_ray(                    /* process a ray result or report error */
72	>	RAY *r,
73	>	int rv
74	>	)
75		{
76		COLORV  *colp;
77
78	<	if (rv == 0)
78	>	if (rv == 0)                    /* no result ready */
79		return(0);
80		if (rv < 0)
81		error(USER, "ray tracing process died");
82		if (r->rno >= distsiz)
83		error(INTERNAL, "bad returned index in process_ray");
84	+	multcolor(r->rcol, r->rcoef);   /* in case it's a source ray */
85		colp = &distarr[r->rno * 3];
86		addcolor(colp, r->rcol);
87	+	if (r->rsrc >= 0 &&             /* remember source contrib. */
88	+	direct_discount != NULL) {
89	+	colp = &direct_discount[r->rno * 3];
90	+	addcolor(colp, r->rcol);
91	+	}
92		return(1);
93		}
94
95
96		void
97	<	raysamp(        /* queue a ray sample */
97	>	raysamp(                        /* queue a ray sample */
98		int  ndx,
99		FVECT  org,
100		FVECT  dir
#	Line 73 | Line 108 | raysamp(       /* queue a ray sample */
108		VCOPY(myRay.rorg, org);
109		VCOPY(myRay.rdir, dir);
110		myRay.rmax = .0;
111	<	rayorigin(&myRay, PRIMARY, NULL, NULL);
111	>	rayorigin(&myRay, PRIMARY|SPECULAR, NULL, NULL);
112		myRay.rno = ndx;
113		/* queue ray, check result */
114		process_ray(&myRay, ray_pqueue(&myRay));
#	Line 81 | Line 116 | raysamp(       /* queue a ray sample */
116
117
118		void
119	+	srcsamps(                       /* sample sources from this surface position */
120	+	struct illum_args *il,
121	+	FVECT org,
122	+	FVECT nrm,
123	+	MAT4 ixfm
124	+	)
125	+	{
126	+	int  nalt=1, nazi=1;
127	+	SRCINDEX  si;
128	+	RAY  sr;
129	+	FVECT   v;
130	+	double  d;
131	+	int  i, j;
132	+	/* get sampling density */
133	+	if (il->sd == NULL && il->sampdens > 0) {
134	+	i = PI * il->sampdens;
135	+	nalt = sqrt(i/PI) + .5;
136	+	nazi = PI*nalt + .5;
137	+	}
138	+	initsrcindex(&si);                      /* loop over (sub)sources */
139	+	for ( ; ; ) {
140	+	VCOPY(sr.rorg, org);            /* pick side to shoot from */
141	+	if (il->sd != NULL) {
142	+	int  sn = si.sn;
143	+	if (si.sp+1 >= si.np) ++sn;
144	+	if (sn >= nsources) break;
145	+	if (source[sn].sflags & SDISTANT)
146	+	d = DOT(source[sn].sloc, nrm);
147	+	else {
148	+	VSUB(v, source[sn].sloc, org);
149	+	d = DOT(v, nrm);
150	+	}
151	+	} else
152	+	d = 1.0;                /* only transmission */
153	+	if (d < 0.0)
154	+	d = -1.0001*il->thick - 5.*FTINY;
155	+	else
156	+	d = 5.*FTINY;
157	+	VSUM(sr.rorg, sr.rorg, nrm, d);
158	+	samplendx++;                    /* increment sample counter */
159	+	if (!srcray(&sr, NULL, &si))
160	+	break;                  /* end of sources */
161	+	/* index direction */
162	+	if (ixfm != NULL)
163	+	multv3(v, sr.rdir, ixfm);
164	+	else
165	+	VCOPY(v, sr.rdir);
166	+	if (il->sd != NULL) {
167	+	i = getBSDF_incndx(il->sd, v);
168	+	if (i < 0)
169	+	continue;       /* must not be important */
170	+	sr.rno = i;
171	+	d = 1.0/getBSDF_incohm(il->sd, i);
172	+	} else {
173	+	if (v[2] >= -FTINY)
174	+	continue;       /* only sample transmission */
175	+	v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2];
176	+	sr.rno = flatindex(v, nalt, nazi);
177	+	d = nalt*nazi*(1./PI) * v[2];
178	+	}
179	+	d *= si.dom;                    /* solid angle correction */
180	+	scalecolor(sr.rcoef, d);
181	+	process_ray(&sr, ray_pqueue(&sr));
182	+	}
183	+	}
184	+
185	+
186	+	void
187		rayclean()                      /* finish all pending rays */
188		{
189		RAY     myRay;
#	Line 127 | Line 230 | rounddir(              /* compute uniform spherical direction */
230		}
231
232
233	<	static void
233	>	void
234		flatdir(                /* compute uniform hemispherical direction */
235	<	register FVECT  dv,
235	>	FVECT  dv,
236		double  alt,
237		double  azi
238		)
#	Line 145 | Line 248 | flatdir(               /* compute uniform hemispherical direction *
248
249
250		int
251	+	flatindex(              /* compute index for hemispherical direction */
252	+	FVECT   dv,
253	+	int     nalt,
254	+	int     nazi
255	+	)
256	+	{
257	+	double  d;
258	+	int     i, j;
259	+
260	+	d = 1.0 - dv[2]*dv[2];
261	+	i = d*nalt;
262	+	d = atan2(dv[1], dv[0]) * (0.5/PI);
263	+	if (d < 0.0) d += 1.0;
264	+	j = d*nazi + 0.5;
265	+	if (j >= nazi) j = 0;
266	+	return(i*nazi + j);
267	+	}
268	+
269	+
270	+	int
271	+	printgeom(              /* print out detailed geometry for BSDF */
272	+	struct BSDF_data *sd,
273	+	char *xfrot,
274	+	FVECT ctr,
275	+	double s1,
276	+	double s2
277	+	)
278	+	{
279	+	static char     mgftemp[] = TEMPLATE;
280	+	char            cmdbuf[64];
281	+	FILE            *fp;
282	+	double          sca;
283	+
284	+	if (sd == NULL || sd->mgf == NULL)
285	+	return(0);
286	+	if (sd->dim[0] <= FTINY || sd->dim[1] <= FTINY)
287	+	return(0);
288	+	if ((s1 > s2) ^ (sd->dim[0] > sd->dim[1])) {
289	+	sca = s1; s1 = s2; s2 = sca;
290	+	}
291	+	s1 /= sd->dim[0];
292	+	s2 /= sd->dim[1];
293	+	sca = s1 > s2 ? s1 : s2;
294	+	strcpy(mgftemp, TEMPLATE);
295	+	if ((fp = fopen(mktemp(mgftemp), "w")) == NULL)
296	+	error(SYSTEM, "cannot create temporary file for MGF");
297	+	/* prepend our transform */
298	+	fprintf(fp, "xf%s -s %.5f -t %.5g %.5g %.5g\n",
299	+	xfrot, sca, ctr[0], ctr[1], ctr[2]);
300	+	/* output given MGF description */
301	+	fputs(sd->mgf, fp);
302	+	fputs("\nxf\n", fp);
303	+	if (fclose(fp) == EOF)
304	+	error(SYSTEM, "error writing MGF temporary file");
305	+	/* execute mgf2rad to convert MGF */
306	+	strcpy(cmdbuf, "mgf2rad ");
307	+	strcpy(cmdbuf+8, mgftemp);
308	+	fflush(stdout);
309	+	system(cmdbuf);
310	+	unlink(mgftemp);                /* clean up */
311	+	return(1);
312	+	}
313	+
314	+
315	+	int
316		my_default(     /* default illum action */
317		OBJREC  *ob,
318		struct illum_args  *il,
#	Line 166 | Line 334 | my_face(               /* make an illum face */
334		char  *nm
335		)
336		{
169	–	#define MAXMISS         (5*n*il->nsamps)
337		int  dim[2];
338	<	int  n, nalt, nazi, h, alti;
338	>	int  n, nalt, nazi, alti;
339		double  sp[2], r1, r2;
340	+	int  h;
341		FVECT  dn, org, dir;
342		FVECT  u, v;
343		double  ur[2], vr[2];
344		MAT4  xfm;
345	<	int  nmisses;
345	>	char  xfrot[64];
346	>	int  nallow;
347		FACE  *fa;
348	<	register int  i, j;
348	>	int  i, j;
349		/* get/check arguments */
350		fa = getface(ob);
351		if (fa->area == 0.0) {
#	Line 185 | Line 354 | my_face(               /* make an illum face */
354		}
355		/* set up sampling */
356		if (il->sd != NULL) {
357	<	if (!getBSDF_xfm(xfm, fa->norm, il->udir)) {
357	>	if (!getBSDF_xfm(xfm, fa->norm, il->udir, xfrot)) {
358		objerror(ob, WARNING, "illegal up direction");
359		freeface(ob);
360		return(my_default(ob, il, nm));
#	Line 227 | Line 396 | my_face(               /* make an illum face */
396		if (r2 < vr[0]) vr[0] = r2;
397		if (r2 > vr[1]) vr[1] = r2;
398		}
399	+	/* output detailed geometry? */
400	+	if (!(il->flags & IL_LIGHT) && il->sd != NULL && il->sd->mgf != NULL &&
401	+	il->thick <= FTINY) {
402	+	for (j = 3; j--; )
403	+	org[j] = .5*(ur[0]+ur[1])*u[j] +
404	+	.5*(vr[0]+vr[1])*v[j] +
405	+	fa->offset*fa->norm[j];
406	+	printgeom(il->sd, xfrot, org, ur[1]-ur[0], vr[1]-vr[0]);
407	+	}
408		dim[0] = random();
409		/* sample polygon */
410	<	nmisses = 0;
410	>	nallow = 5*n*il->nsamps;
411		for (dim[1] = 0; dim[1] < n; dim[1]++)
412		for (i = 0; i < il->nsamps; i++) {
413	<	/* random direction */
413	>	/* randomize direction */
414		h = ilhash(dim, 2) + i;
415		if (il->sd != NULL) {
416		r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
#	Line 246 | Line 424 | my_face(               /* make an illum face */
424		dir[j] = -dn[0]*u[j] - dn[1]*v[j] -
425		dn[2]*fa->norm[j];
426		}
427	<	/* random location */
427	>	/* randomize location */
428		do {
429	<	multisamp(sp, 2, urand(h+4862+nmisses));
429	>	multisamp(sp, 2, urand(h+4862+nallow));
430		r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
431		r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
432		for (j = 0; j < 3; j++)
433		org[j] = r1*u[j] + r2*v[j]
434		+ fa->offset*fa->norm[j];
435	<	} while (!inface(org, fa) && nmisses++ < MAXMISS);
436	<	if (nmisses > MAXMISS) {
435	>	} while (!inface(org, fa) && nallow-- > 0);
436	>	if (nallow < 0) {
437		objerror(ob, WARNING, "bad aspect");
438		rayclean();
439		freeface(ob);
440		return(my_default(ob, il, nm));
441		}
442		if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY)
443	<	r1 = -1.0001*il->thick - .0001;
443	>	r1 = -1.0001*il->thick - 5.*FTINY;
444		else
445	<	r1 = .0001;
445	>	r1 = 5.*FTINY;
446		for (j = 0; j < 3; j++)
447		org[j] += r1*fa->norm[j];
448		/* send sample */
449		raysamp(dim[1], org, dir);
450		}
451	+	/* add in direct component? */
452	+	if (il->flags & IL_LIGHT || il->sd != NULL) {
453	+	MAT4    ixfm;
454	+	if (il->sd == NULL) {
455	+	for (i = 3; i--; ) {
456	+	ixfm[i][0] = u[i];
457	+	ixfm[i][1] = v[i];
458	+	ixfm[i][2] = fa->norm[i];
459	+	ixfm[i][3] = 0.;
460	+	}
461	+	ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
462	+	ixfm[3][3] = 1.;
463	+	} else {
464	+	if (!invmat4(ixfm, xfm))
465	+	objerror(ob, INTERNAL,
466	+	"cannot invert BSDF transform");
467	+	if (!(il->flags & IL_LIGHT))
468	+	new_discount();
469	+	}
470	+	dim[0] = random();
471	+	nallow = 10*il->nsamps;
472	+	for (i = 0; i < il->nsamps; i++) {
473	+	/* randomize location */
474	+	h = dim[0] + samplendx++;
475	+	do {
476	+	multisamp(sp, 2, urand(h+nallow));
477	+	r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
478	+	r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
479	+	for (j = 0; j < 3; j++)
480	+	org[j] = r1*u[j] + r2*v[j]
481	+	+ fa->offset*fa->norm[j];
482	+	} while (!inface(org, fa) && nallow-- > 0);
483	+	if (nallow < 0) {
484	+	objerror(ob, WARNING, "bad aspect");
485	+	rayclean();
486	+	freeface(ob);
487	+	return(my_default(ob, il, nm));
488	+	}
489	+	/* sample source rays */
490	+	srcsamps(il, org, fa->norm, ixfm);
491	+	}
492	+	}
493	+	/* wait for all rays to finish */
494		rayclean();
495		if (il->sd != NULL) {   /* run distribution through BSDF */
496		nalt = sqrt(il->sd->nout/PI) + .5;
497		nazi = PI*nalt + .5;
498		redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm);
499	+	done_discount();
500	+	if (!il->sampdens)
501	+	il->sampdens = nalt*nazi/PI + .999;
502		}
503		/* write out the face and its distribution */
504		if (average(il, distarr, n)) {
#	Line 286 | Line 510 | my_face(               /* make an illum face */
510		/* clean up */
511		freeface(ob);
512		return(0);
289	–	#undef MAXMISS
513		}
514
515
#	Line 344 | Line 567 | my_sphere(     /* make an illum sphere */
567		/* send sample */
568		raysamp(dim[1]*nazi+dim[2], org, dir);
569		}
570	+	/* wait for all rays to finish */
571		rayclean();
572		/* write out the sphere and its distribution */
573		if (average(il, distarr, n)) {
#	Line 374 | Line 598 | my_ring(               /* make an illum ring */
598		FVECT  u, v;
599		MAT4  xfm;
600		CONE  *co;
601	<	register int  i, j;
601	>	int  i, j;
602		/* get/check arguments */
603		co = getcone(ob, 0);
604		/* set up sampling */
605		if (il->sd != NULL) {
606	<	if (!getBSDF_xfm(xfm, co->ad, il->udir)) {
606	>	if (!getBSDF_xfm(xfm, co->ad, il->udir, NULL)) {
607		objerror(ob, WARNING, "illegal up direction");
608		freecone(ob);
609		return(my_default(ob, il, nm));
#	Line 403 | Line 627 | my_ring(               /* make an illum ring */
627		for (i = 0; i < il->nsamps; i++) {
628		/* next sample point */
629		h = ilhash(dim,2) + i;
630	<	/* random direction */
630	>	/* randomize direction */
631		if (il->sd != NULL) {
632		r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
633		} else {
#	Line 413 | Line 637 | my_ring(               /* make an illum ring */
637		r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
638		flatdir(dn, r1, r2);
639		for (j = 0; j < 3; j++)
640	<	dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
640	>	dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
641		}
642	<	/* random location */
642	>	/* randomize location */
643		multisamp(sp, 2, urand(h+8371));
644		r3 = sqrt(CO_R0(co)*CO_R0(co) +
645		sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
#	Line 423 | Line 647 | my_ring(               /* make an illum ring */
647		r1 = r3*cos(r2);
648		r2 = r3*sin(r2);
649		if (il->sd != NULL && DOT(dir, co->ad) < -FTINY)
650	<	r3 = -1.0001*il->thick - .0001;
650	>	r3 = -1.0001*il->thick - 5.*FTINY;
651		else
652	<	r3 = .0001;
652	>	r3 = 5.*FTINY;
653		for (j = 0; j < 3; j++)
654		org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
655		r3*co->ad[j];
656		/* send sample */
657		raysamp(dim[1], org, dir);
658		}
659	+	/* add in direct component? */
660	+	if (il->flags & IL_LIGHT || il->sd != NULL) {
661	+	MAT4    ixfm;
662	+	if (il->sd == NULL) {
663	+	for (i = 3; i--; ) {
664	+	ixfm[i][0] = u[i];
665	+	ixfm[i][1] = v[i];
666	+	ixfm[i][2] = co->ad[i];
667	+	ixfm[i][3] = 0.;
668	+	}
669	+	ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
670	+	ixfm[3][3] = 1.;
671	+	} else {
672	+	if (!invmat4(ixfm, xfm))
673	+	objerror(ob, INTERNAL,
674	+	"cannot invert BSDF transform");
675	+	if (!(il->flags & IL_LIGHT))
676	+	new_discount();
677	+	}
678	+	dim[0] = random();
679	+	for (i = 0; i < il->nsamps; i++) {
680	+	/* randomize location */
681	+	h = dim[0] + samplendx++;
682	+	multisamp(sp, 2, urand(h));
683	+	r3 = sqrt(CO_R0(co)*CO_R0(co) +
684	+	sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
685	+	r2 = 2.*PI*sp[1];
686	+	r1 = r3*cos(r2);
687	+	r2 = r3*sin(r2);
688	+	for (j = 0; j < 3; j++)
689	+	org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j];
690	+	/* sample source rays */
691	+	srcsamps(il, org, co->ad, ixfm);
692	+	}
693	+	}
694	+	/* wait for all rays to finish */
695		rayclean();
696		if (il->sd != NULL) {   /* run distribution through BSDF */
697		nalt = sqrt(il->sd->nout/PI) + .5;
698		nazi = PI*nalt + .5;
699		redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
700	+	done_discount();
701	+	if (!il->sampdens)
702	+	il->sampdens = nalt*nazi/PI + .999;
703		}
704		/* write out the ring and its distribution */
705		if (average(il, distarr, n)) {
#	Line 448 | Line 711 | my_ring(               /* make an illum ring */
711		/* clean up */
712		freecone(ob);
713		return(1);
714	+	}
715	+
716	+
717	+	void
718	+	redistribute(           /* pass distarr ray sums through BSDF */
719	+	struct BSDF_data *b,
720	+	int nalt,
721	+	int nazi,
722	+	FVECT u,
723	+	FVECT v,
724	+	FVECT w,
725	+	MAT4 xm
726	+	)
727	+	{
728	+	int     nout = 0;
729	+	MAT4    mymat, inmat;
730	+	COLORV  *idist;
731	+	COLORV  *cp;
732	+	FVECT   dv;
733	+	double  wt;
734	+	int     i, j, k, c, o;
735	+	COLOR   col, cinc;
736	+	/* copy incoming distribution */
737	+	if (b->ninc > distsiz)
738	+	error(INTERNAL, "error 1 in redistribute");
739	+	idist = (COLORV *)malloc(sizeof(COLOR)*b->ninc);
740	+	if (idist == NULL)
741	+	error(SYSTEM, "out of memory in redistribute");
742	+	memcpy(idist, distarr, sizeof(COLOR)*b->ninc);
743	+	/* compose direction transform */
744	+	for (i = 3; i--; ) {
745	+	mymat[i][0] = u[i];
746	+	mymat[i][1] = v[i];
747	+	mymat[i][2] = w[i];
748	+	mymat[i][3] = 0.;
749	+	}
750	+	mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
751	+	mymat[3][3] = 1.;
752	+	if (xm != NULL)
753	+	multmat4(mymat, xm, mymat);
754	+	for (i = 3; i--; ) {            /* make sure it's normalized */
755	+	wt = 1./sqrt(   mymat[0][i]*mymat[0][i] +
756	+	mymat[1][i]*mymat[1][i] +
757	+	mymat[2][i]*mymat[2][i] );
758	+	for (j = 3; j--; )
759	+	mymat[j][i] *= wt;
760	+	}
761	+	if (!invmat4(inmat, mymat))     /* need inverse as well */
762	+	error(INTERNAL, "cannot invert BSDF transform");
763	+	newdist(nalt*nazi);             /* resample distribution */
764	+	for (i = b->ninc; i--; ) {
765	+	int     direct_out = -1;
766	+	COLOR   cdir;
767	+	getBSDF_incvec(dv, b, i);       /* compute incident irrad. */
768	+	multv3(dv, dv, mymat);
769	+	if (dv[2] < 0.0) {
770	+	dv[0] = -dv[0]; dv[1] = -dv[1]; dv[2] = -dv[2];
771	+	direct_out += (direct_discount != NULL);
772	+	}
773	+	wt = getBSDF_incohm(b, i);
774	+	wt *= dv[2];                    /* solid_angle*cosine(theta) */
775	+	cp = &idist[3*i];
776	+	copycolor(cinc, cp);
777	+	scalecolor(cinc, wt);
778	+	if (!direct_out) {              /* discount direct contr. */
779	+	cp = &direct_discount[3*i];
780	+	copycolor(cdir, cp);
781	+	scalecolor(cdir, -wt);
782	+	if (b->nout != b->ninc)
783	+	direct_out = flatindex(dv, nalt, nazi);
784	+	else
785	+	direct_out = i; /* assumes dist. mirroring */
786	+	}
787	+	for (k = nalt; k--; )           /* loop over distribution */
788	+	for (j = nazi; j--; ) {
789	+	int rstart = random();
790	+	for (c = NBSDFSAMPS; c--; ) {
791	+	double  sp[2];
792	+	multisamp(sp, 2, urand(rstart+c));
793	+	flatdir(dv, (k + sp[0])/nalt,
794	+	(j + .5 - sp[1])/nazi);
795	+	multv3(dv, dv, inmat);
796	+	/* evaluate BSDF @ outgoing */
797	+	o = getBSDF_outndx(b, dv);
798	+	if (o < 0) {
799	+	nout++;
800	+	continue;
801	+	}
802	+	wt = BSDF_value(b, i, o) * (1./NBSDFSAMPS);
803	+	copycolor(col, cinc);
804	+	if (b->nout != b->ninc)
805	+	o = k*nazi + j;
806	+	if (o == direct_out)
807	+	addcolor(col, cdir);    /* minus direct */
808	+	scalecolor(col, wt);
809	+	cp = &distarr[3*(k*nazi + j)];
810	+	addcolor(cp, col);      /* sum into distribution */
811	+	}
812	+	}
813	+	}
814	+	free(idist);                    /* free temp space */
815	+	if (nout) {
816	+	sprintf(errmsg, "missing %.1f%% of BSDF directions",
817	+	100.*nout/(b->ninc*nalt*nazi*NBSDFSAMPS));
818	+	error(WARNING, errmsg);
819	+	}
820		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)