ViewVC Help

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

root/Development/ray/src/rt/srcsamp.c

Comparing ray/src/rt/srcsamp.c (file contents):
Revision 2.2 by greg, Mon Nov 25 09:52:49 1991 UTC vs.
Revision 2.23 by greg, Fri Dec 13 00:50:55 2024 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1991 Regents of the University of California */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ LBL";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		* Source sampling routines
6	+	*
7	+	*  External symbols declared in source.h
8		*/
9
10	+	#include "copyright.h"
11	+
12		#include  "ray.h"
13
14		#include  "source.h"
15
16		#include  "random.h"
17
18	+	#ifdef SSKIPOPT
19	+	/* The following table is used for skipping sources */
20	+	static uby8     *srcskipflags = NULL;           /* source inclusion lookup */
21	+	static int      ssf_count = 0;                  /* number of flag entries */
22	+	static int      ssf_max = 0;                    /* current array size */
23	+	static uby8     *ssf_noskip = NULL;             /* set of zero flags */
24
25	+	uby8            *ssf_select = NULL;             /* sources we may skip */
26	+
27	+	/* Find/allocate source skip flag entry (free all if NULL) */
28	+	int
29	+	sskip_rsi(uby8 *flags)
30	+	{
31	+	uby8    *flp;
32	+	int     i;
33	+
34	+	if (flags == NULL) {            /* means clear all */
35	+	efree(srcskipflags); srcskipflags = NULL;
36	+	ssf_count = ssf_max = 0;
37	+	sskip_free(ssf_noskip);
38	+	sskip_free(ssf_select);
39	+	return(0);
40	+	}
41	+	if (ssf_noskip == NULL)         /* first call? */
42	+	ssf_noskip = sskip_new();
43	+
44	+	if (sskip_eq(flags, ssf_noskip))
45	+	return(-1);             /* nothing to skip */
46	+	/* search recent entries */
47	+	flp = srcskipflags + ssf_count*SSKIPFLSIZ;
48	+	for (i = ssf_count; i-- > 0; )
49	+	if (sskip_eq(flp -= SSKIPFLSIZ, flags))
50	+	return(-2-i);   /* found it! */
51	+	/* else tack on new entry */
52	+	if (ssf_count >= ssf_max) {     /* need more space? */
53	+	ssf_max = ssf_count + (ssf_count>>2) + 64;
54	+	if (ssf_max <= ssf_count &&
55	+	(ssf_max = ssf_count+1024) <= ssf_count)
56	+	error(SYSTEM, "out of space in sskip_rsi()");
57	+
58	+	srcskipflags = (uby8 *)erealloc(srcskipflags,
59	+	ssf_max*SSKIPFLSIZ);
60	+	}
61	+	sskip_cpy(srcskipflags + ssf_count*SSKIPFLSIZ, flags);
62	+
63	+	return(-2 - ssf_count++);       /* return index (< -1) */
64	+	}
65	+
66	+	/* Get skip flags associated with RAY rsrc index (or NULL) */
67	+	uby8 *
68	+	sskip_flags(int rsi)
69	+	{
70	+	if (rsi >= -1)
71	+	return(ssf_noskip);
72	+
73	+	if ((rsi = -2 - rsi) >= ssf_count)
74	+	error(CONSISTENCY, "bad index to sskip_flags()");
75	+
76	+	return(srcskipflags + rsi*SSKIPFLSIZ);
77	+	}
78	+
79	+	/* OR in a second set of flags into a first */
80	+	void
81	+	sskip_addflags(uby8 *dfl, const uby8 *sfl)
82	+	{
83	+	int     nb = SSKIPFLSIZ;
84	+
85	+	while (nb--)
86	+	*dfl++ |= *sfl++;
87	+	}
88	+	#endif
89	+
90	+	int
91	+	srcskip(                        /* pre-emptive test for source to skip */
92	+	int  sn,
93	+	RAY  *r
94	+	)
95	+	{
96	+	SRCREC  *sp = source + sn;
97	+
98	+	if (sp->sflags & SSKIP)
99	+	return(1);
100	+	#ifdef SSKIPOPT
101	+	if (r->rsrc < -1 &&     /* ray has custom skip flags? */
102	+	sskip_chk(sskip_flags(r->rsrc), sn))
103	+	return(1);
104	+	#endif
105	+	if ((sp->sflags & (SPROX|SDISTANT)) != SPROX)
106	+	return(0);
107	+
108	+	return(dist2(r->rorg, sp->sloc) >
109	+	(sp->sl.prox + sp->srad)*(sp->sl.prox + sp->srad));
110	+	}
111	+
112		double
113	<	nextssamp(r, si)                /* compute sample for source, rtn. distance */
114	<	register RAY  *r;               /* origin is read, direction is set */
115	<	register SRCINDEX  *si;         /* source index (modified to current) */
113	>	nextssamp(                      /* compute sample for source, rtn. distance */
114	>	RAY  *r,                /* origin is read, direction is set */
115	>	SRCINDEX  *si           /* source index (modified to current) */
116	>	)
117		{
118		int  cent[3], size[3], parr[2];
119	<	FVECT  vpos;
119	>	SRCREC  *srcp;
120	>	double  vpos[3];
121		double  d;
122	<	register int  i;
122	>	int  i;
123		nextsample:
124		while (++si->sp >= si->np) {    /* get next sample */
125		if (++si->sn >= nsources)
126		return(0.0);    /* no more */
127	<	if (source[si->sn].sflags & SSKIP)
127	>	if (srcskip(si->sn, r))
128		si->np = 0;
129		else if (srcsizerat <= FTINY)
130		nopart(si, r);
#	Line 47 | Line 143 | nextsample:
143		if (!skipparts(cent, size, parr, si->spt))
144		error(CONSISTENCY, "bad source partition in nextssamp");
145		/* compute sample */
146	+	srcp = source + si->sn;
147		if (dstrsrc > FTINY) {                  /* jitter sample */
148		dimlist[ndims] = si->sn + 8831;
149		dimlist[ndims+1] = si->sp + 3109;
150		d = urand(ilhash(dimlist,ndims+2)+samplendx);
151	<	if (source[si->sn].sflags & SFLAT) {
151	>	if (srcp->sflags & SFLAT) {
152		multisamp(vpos, 2, d);
153	<	vpos[2] = 0.5;
153	>	vpos[SW] = 0.5;
154		} else
155		multisamp(vpos, 3, d);
156		for (i = 0; i < 3; i++)
157		vpos[i] = dstrsrc * (1. - 2.*vpos[i]) *
158	<	(double)size[i]/MAXSPART;
158	>	(double)size[i]*(1.0/MAXSPART);
159		} else
160		vpos[0] = vpos[1] = vpos[2] = 0.0;
161
162	<	for (i = 0; i < 3; i++)
163	<	vpos[i] += (double)cent[i]/MAXSPART;
162	>	VSUM(vpos, vpos, cent, 1.0/MAXSPART);
163	>	/* avoid circular aiming failures */
164	>	if ((srcp->sflags & SCIR) && (si->np > 1) | (dstrsrc > 0.7)) {
165	>	FVECT   trim;
166	>	if (srcp->sflags & (SFLAT|SDISTANT)) {
167	>	d = 1.12837917;         /* correct setflatss() */
168	>	trim[SU] = d*sqrt(1.0 - 0.5*vpos[SV]*vpos[SV]);
169	>	trim[SV] = d*sqrt(1.0 - 0.5*vpos[SU]*vpos[SU]);
170	>	trim[SW] = 0.0;
171	>	} else {
172	>	trim[SW] = trim[SU] = vpos[SU]*vpos[SU];
173	>	d = vpos[SV]*vpos[SV];
174	>	if (d > trim[SW]) trim[SW] = d;
175	>	trim[SU] += d;
176	>	d = vpos[SW]*vpos[SW];
177	>	if (d > trim[SW]) trim[SW] = d;
178	>	trim[SU] += d;
179	>	if (trim[SU] > FTINY*FTINY) {
180	>	d = 1.0/0.7236; /* correct sphsetsrc() */
181	>	trim[SW] = trim[SV] = trim[SU] =
182	>	d*sqrt(trim[SW]/trim[SU]);
183	>	} else
184	>	trim[SW] = trim[SV] = trim[SU] = 0.0;
185	>	}
186	>	for (i = 0; i < 3; i++)
187	>	vpos[i] *= trim[i];
188	>	}
189		/* compute direction */
190		for (i = 0; i < 3; i++)
191	<	r->rdir[i] = source[si->sn].sloc[i] +
192	<	vpos[SU]*source[si->sn].ss[SU][i] +
193	<	vpos[SV]*source[si->sn].ss[SV][i] +
194	<	vpos[SW]*source[si->sn].ss[SW][i];
191	>	r->rdir[i] = srcp->sloc[i] +
192	>	vpos[SU]*srcp->ss[SU][i] +
193	>	vpos[SV]*srcp->ss[SV][i] +
194	>	vpos[SW]*srcp->ss[SW][i];
195
196	<	if (!(source[si->sn].sflags & SDISTANT))
197	<	for (i = 0; i < 3; i++)
76	<	r->rdir[i] -= r->rorg[i];
196	>	if (!(srcp->sflags & SDISTANT))
197	>	VSUB(r->rdir, r->rdir, r->rorg);
198		/* compute distance */
199		if ((d = normalize(r->rdir)) == 0.0)
200		goto nextsample;                /* at source! */
201
202		/* compute sample size */
203	<	si->dom  = source[si->sn].ss2;
204	<	if (source[si->sn].sflags & SFLAT) {
205	<	si->dom *= sflatform(si->sn, r->rdir);
206	<	si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART);
207	<	} else if (source[si->sn].sflags & SCYL) {
208	<	si->dom *= scylform(si->sn, r->rdir);
88	<	si->dom *= (double)size[SU]/MAXSPART;
203	>	if (srcp->sflags & SFLAT) {
204	>	si->dom = sflatform(si->sn, r->rdir);
205	>	si->dom *= size[SU]*size[SV]*(1.0/MAXSPART/MAXSPART);
206	>	} else if (srcp->sflags & SCYL) {
207	>	si->dom = scylform(si->sn, r->rdir);
208	>	si->dom *= size[SU]*(1.0/MAXSPART);
209		} else {
210	<	si->dom *= (double)(size[SU]*size[SV]*size[SW]) /
211	<	(MAXSPART*MAXSPART*MAXSPART) ;
210	>	si->dom = size[SU]*size[SV]*(double)size[SW] *
211	>	(1.0/MAXSPART/MAXSPART/MAXSPART) ;
212		}
213	<	if (source[si->sn].sflags & SDISTANT)
213	>	if (srcp->sflags & SDISTANT) {
214	>	si->dom *= srcp->ss2;
215		return(FHUGE);
216	<	if (si->dom <= FTINY)
216	>	}
217	>	if (si->dom <= 1e-4)
218		goto nextsample;                /* behind source? */
219	<	si->dom /= d*d;
219	>	si->dom *= srcp->ss2/(d*d);
220		return(d);              /* sample OK, return distance */
221		}
222
223
224	<	skipparts(ct, sz, pp, pt)               /* skip to requested partition */
225	<	int  ct[3], sz[3];              /* center and size of partition (returned) */
226	<	register int  pp[2];            /* current index, number to skip (modified) */
227	<	unsigned char  *pt;             /* partition array */
224	>	int
225	>	skipparts(                      /* skip to requested partition */
226	>	int  ct[3],
227	>	int  sz[3],             /* center and size of partition (returned) */
228	>	int  pp[2],             /* current index, number to skip (modified) */
229	>	unsigned char  *pt      /* partition array */
230	>	)
231		{
232	<	register int  p;
232	>	int  p;
233		/* check this partition */
234		p = spart(pt, pp[0]);
235		pp[0]++;
236	<	if (p == S0)                    /* leaf partition */
236	>	if (p == S0) {                  /* leaf partition */
237		if (pp[1]) {
238		pp[1]--;
239		return(0);      /* not there yet */
240		} else
241		return(1);      /* we've arrived */
242	+	}
243		/* else check lower */
244		sz[p] >>= 1;
245		ct[p] -= sz[p];
#	Line 130 | Line 256 | unsigned char  *pt;            /* partition array */
256		}
257
258
259	<	nopart(si, r)                   /* single source partition */
260	<	register SRCINDEX  *si;
261	<	RAY  *r;
259	>	void
260	>	nopart(                         /* single source partition */
261	>	SRCINDEX  *si,
262	>	RAY  *r
263	>	)
264		{
265		clrpart(si->spt);
266		setpart(si->spt, 0, S0);
#	Line 140 | Line 268 | RAY  *r;
268		}
269
270
143	–	cylpart(si, r)                  /* partition a cylinder */
144	–	SRCINDEX  *si;
145	–	register RAY  *r;
146	–	{
147	–	double  dist2, safedist2, dist2cent, rad2;
148	–	FVECT  v;
149	–	register SRCREC  *sp;
150	–	int  pi;
151	–	/* first check point location */
152	–	clrpart(si->spt);
153	–	sp = source + si->sn;
154	–	rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
155	–	v[0] = r->rorg[0] - sp->sloc[0];
156	–	v[1] = r->rorg[1] - sp->sloc[1];
157	–	v[2] = r->rorg[2] - sp->sloc[2];
158	–	dist2 = DOT(v,sp->ss[SU]);
159	–	safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
160	–	dist2 *= dist2 / safedist2;
161	–	dist2cent = DOT(v,v);
162	–	dist2 = dist2cent - dist2;
163	–	if (dist2 <= rad2) {            /* point inside extended cylinder */
164	–	si->np = 0;
165	–	return;
166	–	}
167	–	safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
168	–	if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
169	–	dist2cent >= safedist2) {       /* or too far */
170	–	setpart(si->spt, 0, S0);
171	–	si->np = 1;
172	–	return;
173	–	}
174	–	pi = 0;
175	–	si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
176	–	sp->sloc, sp->ss[SU], safedist2);
177	–	}
178	–
179	–
271		static int
272	<	cyl_partit(ro, pt, pi, mp, cent, axis, d2)      /* slice a cylinder */
273	<	FVECT  ro;
274	<	unsigned char  *pt;
275	<	register int  *pi;
276	<	int  mp;
277	<	FVECT  cent, axis;
278	<	double  d2;
272	>	cyl_partit(                             /* slice a cylinder */
273	>	FVECT  ro,
274	>	unsigned char  *pt,
275	>	int  *pi,
276	>	int  mp,
277	>	FVECT  cent,
278	>	FVECT  axis,
279	>	double  d2
280	>	)
281		{
282		FVECT  newct, newax;
283		int  npl, npu;
#	Line 216 | Line 309 | double  d2;
309		}
310
311
312	<	flatpart(si, r)                         /* partition a flat source */
313	<	register SRCINDEX  *si;
314	<	register RAY  *r;
312	>	void
313	>	cylpart(                        /* partition a cylinder */
314	>	SRCINDEX  *si,
315	>	RAY  *r
316	>	)
317		{
318	<	register FLOAT  *vp;
318	>	double  dist2, safedist2, dist2cent, rad2;
319		FVECT  v;
320	<	double  du2, dv2;
320	>	SRCREC  *sp;
321		int  pi;
322	<
322	>	/* first check point location */
323		clrpart(si->spt);
324	<	vp = source[si->sn].sloc;
325	<	v[0] = r->rorg[0] - vp[0];
326	<	v[1] = r->rorg[1] - vp[1];
327	<	v[2] = r->rorg[2] - vp[2];
328	<	vp = source[si->sn].snorm;
329	<	if (DOT(v,vp) <= FTINY) {       /* behind source */
324	>	sp = source + si->sn;
325	>	rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
326	>	v[0] = r->rorg[0] - sp->sloc[0];
327	>	v[1] = r->rorg[1] - sp->sloc[1];
328	>	v[2] = r->rorg[2] - sp->sloc[2];
329	>	dist2 = DOT(v,sp->ss[SU]);
330	>	safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
331	>	dist2 *= dist2 / safedist2;
332	>	dist2cent = DOT(v,v);
333	>	dist2 = dist2cent - dist2;
334	>	if (dist2 <= rad2) {            /* point inside extended cylinder */
335		si->np = 0;
336		return;
337		}
338	<	dv2 = 2.*r->rweight/srcsizerat;
339	<	dv2 *= dv2;
340	<	vp = source[si->sn].ss[SU];
341	<	du2 = dv2 * DOT(vp,vp);
342	<	vp = source[si->sn].ss[SV];
343	<	dv2 *= DOT(vp,vp);
338	>	safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
339	>	if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
340	>	dist2cent >= safedist2) {       /* or too far */
341	>	setpart(si->spt, 0, S0);
342	>	si->np = 1;
343	>	return;
344	>	}
345		pi = 0;
346	<	si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
347	<	source[si->sn].sloc,
247	<	source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
346	>	si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
347	>	sp->sloc, sp->ss[SU], safedist2);
348		}
349
350
351		static int
352	<	flt_partit(ro, pt, pi, mp, cent, u, v, du2, dv2)        /* partition flatty */
353	<	FVECT  ro;
354	<	unsigned char  *pt;
355	<	register int  *pi;
356	<	int  mp;
357	<	FVECT  cent, u, v;
358	<	double  du2, dv2;
352	>	flt_partit(                             /* partition flatty */
353	>	FVECT  ro,
354	>	unsigned char  *pt,
355	>	int  *pi,
356	>	int  mp,
357	>	FVECT  cent,
358	>	FVECT  u,
359	>	FVECT  v,
360	>	double  du2,
361	>	double  dv2
362	>	)
363		{
364		double  d2;
365		FVECT  newct, newax;
#	Line 299 | Line 403 | double  du2, dv2;
403		}
404
405
406	+	void
407	+	flatpart(                               /* partition a flat source */
408	+	SRCINDEX  *si,
409	+	RAY  *r
410	+	)
411	+	{
412	+	RREAL  *vp;
413	+	FVECT  v;
414	+	double  du2, dv2;
415	+	int  pi;
416	+
417	+	clrpart(si->spt);
418	+	vp = source[si->sn].sloc;
419	+	v[0] = r->rorg[0] - vp[0];
420	+	v[1] = r->rorg[1] - vp[1];
421	+	v[2] = r->rorg[2] - vp[2];
422	+	vp = source[si->sn].snorm;
423	+	if (DOT(v,vp) <= 0.) {          /* behind source */
424	+	si->np = 0;
425	+	return;
426	+	}
427	+	dv2 = 2.*r->rweight/srcsizerat;
428	+	dv2 *= dv2;
429	+	vp = source[si->sn].ss[SU];
430	+	du2 = dv2 * DOT(vp,vp);
431	+	vp = source[si->sn].ss[SV];
432	+	dv2 *= DOT(vp,vp);
433	+	pi = 0;
434	+	si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
435	+	source[si->sn].sloc,
436	+	source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
437	+	}
438	+
439	+
440		double
441	<	scylform(sn, dir)               /* compute cosine for cylinder's projection */
442	<	int  sn;
443	<	register FVECT  dir;            /* assume normalized */
441	>	scylform(                       /* compute cosine for cylinder's projection */
442	>	int  sn,
443	>	FVECT  dir              /* assume normalized */
444	>	)
445		{
446	<	register FLOAT  *dv;
446	>	RREAL  *dv;
447		double  d;
448
449		dv = source[sn].ss[SU];

Diff Legend

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