ViewVC Help

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

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

Comparing ray/src/gen/mksource.c (file contents):
Revision 2.3 by schorsch, Mon Sep 19 12:32:12 2005 UTC vs.
Revision 2.13 by greg, Tue Jun 3 21:31:51 2025 UTC

#	Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7
8		#include "ray.h"
9		#include "random.h"
10	+	#include "paths.h"
11		#include "resolu.h"
12
13		#define NTRUNKBR        4               /* number of branches at trunk */
#	Line 17 | Line 18 | static const char RCSid[] = "$Id$";
18		/* Data structure for geodesic samples */
19
20		typedef struct tritree {
21	<	FVECT           gdv[3];         /* spherical triangle vertex direc. */
22	<	FVECT           sd;             /* sample direction if leaf */
21	>	int32           gdv[3];         /* spherical triangle vertex direc. */
22	>	int32           sd;             /* sample direction if leaf */
23		struct tritree  *kid;           /* 4 children if branch node */
24		COLR            val;            /* sampled color value */
25		} TRITREE;
26
27		typedef struct lostlight {
28		struct lostlight        *next;  /* next in list */
29	<	FVECT           sd;             /* lost source direction */
29	>	int32           sd;             /* lost source direction */
30		COLOR           intens;         /* output times solid angle */
31		} LOSTLIGHT;
32
32	–	char    *progname;
33	–
33		FVECT   scene_cent;             /* center of octree cube */
34		RREAL   scene_rad;              /* radius to get outside cube from center */
35
#	Line 56 | Line 55 | vol_sign(const FVECT v1, const FVECT v2, const FVECT v
55
56		/* Is the given direction contained within the specified spherical triangle? */
57		int
58	<	intriv(FVECT tri[3], const FVECT sdir)
58	>	intriv(const int32 trid[3], const FVECT sdir)
59		{
60		int     sv[3];
61	+	FVECT   tri[3];
62
63	+	decodedir(tri[0], trid[0]);
64	+	decodedir(tri[1], trid[1]);
65	+	decodedir(tri[2], trid[2]);
66		sv[0] = vol_sign(sdir, tri[0], tri[1]);
67		sv[1] = vol_sign(sdir, tri[1], tri[2]);
68		sv[2] = vol_sign(sdir, tri[2], tri[0]);
#	Line 91 | Line 94 | leafsample(TRITREE *leaf)
94		{
95		RAY     myray;
96		RREAL   wt[3];
97	+	FVECT   sdir;
98		int     i, j;
99		/* random point on triangle */
100		i = random() % 3;
#	Line 98 | Line 102 | leafsample(TRITREE *leaf)
102		j = random() & 1;
103		wt[(i+2-j)%3] = 1. - wt[i] -
104		(wt[(i+1+j)%3] = (1.-wt[i])*frandom());
105	<	leaf->sd[0] = leaf->sd[1] = leaf->sd[2] = .0;
106	<	for (i = 0; i < 3; i++)
107	<	VSUM(leaf->sd, leaf->sd, leaf->gdv[i], wt[i]);
108	<	normalize(leaf->sd);                    /* record sample direction */
105	>	sdir[0] = sdir[1] = sdir[2] = .0;
106	>	for (i = 0; i < 3; i++) {
107	>	FVECT   vt;
108	>	decodedir(vt, leaf->gdv[i]);
109	>	VSUM(sdir, sdir, vt, wt[i]);
110	>	}
111	>	normalize(sdir);                        /* record sample direction */
112	>	leaf->sd = encodedir(sdir);
113		/* evaluate at inf. */
114	<	VSUM(myray.rorg, scene_cent, leaf->sd, scene_rad);
115	<	VCOPY(myray.rdir, leaf->sd);
114	>	VSUM(myray.rorg, scene_cent, sdir, scene_rad);
115	>	VCOPY(myray.rdir, sdir);
116		myray.rmax = 0.;
117		ray_trace(&myray);
118	<	setcolr(leaf->val, colval(myray.rcol,RED),
111	<	colval(myray.rcol,GRN),
112	<	colval(myray.rcol,BLU));
118	>	scolor_colr(leaf->val, myray.rcol);
119		}
120
121		/* Initialize a branch node contained in the given spherical triangle */
122		void
123	<	subdivide(TRITREE *branch, FVECT dv[3])
123	>	subdivide(TRITREE *branch, const int32 dv[3])
124		{
125	<	FVECT   sdv[3];
125	>	FVECT   dvv[3], sdv[3];
126	>	int32   sd[3];
127		int     i;
128
129	<	for (i = 0; i < 3; i++)         /* copy spherical triangle */
130	<	VCOPY(branch->gdv[i], dv[i]);
129	>	for (i = 0; i < 3; i++) {       /* copy spherical triangle */
130	>	branch->gdv[i] = dv[i];
131	>	decodedir(dvv[i], dv[i]);
132	>	}
133		for (i = 0; i < 3; i++) {       /* create new vertices */
134		int     j = (i+1)%3;
135	<	VADD(sdv[i], dv[i], dv[j]);
135	>	VADD(sdv[i], dvv[i], dvv[j]);
136		normalize(sdv[i]);
137	+	sd[i] = encodedir(sdv[i]);
138		}
139		/* allocate leaves */
140		branch->kid = (TRITREE *)calloc(4, sizeof(TRITREE));
141		if (branch->kid == NULL)
142		error(SYSTEM, "out of memory in subdivide()");
143		/* assign subtriangle directions */
144	<	VCOPY(branch->kid[0].gdv[0], dv[0]);
145	<	VCOPY(branch->kid[0].gdv[1], sdv[0]);
146	<	VCOPY(branch->kid[0].gdv[2], sdv[2]);
147	<	VCOPY(branch->kid[1].gdv[0], sdv[0]);
148	<	VCOPY(branch->kid[1].gdv[1], dv[1]);
149	<	VCOPY(branch->kid[1].gdv[2], sdv[1]);
150	<	VCOPY(branch->kid[2].gdv[0], sdv[1]);
151	<	VCOPY(branch->kid[2].gdv[1], dv[2]);
152	<	VCOPY(branch->kid[2].gdv[2], sdv[2]);
153	<	VCOPY(branch->kid[3].gdv[0], sdv[0]);
154	<	VCOPY(branch->kid[3].gdv[1], sdv[1]);
155	<	VCOPY(branch->kid[3].gdv[2], sdv[2]);
144	>	branch->kid[0].gdv[0] = dv[0];
145	>	branch->kid[0].gdv[1] = sd[0];
146	>	branch->kid[0].gdv[2] = sd[2];
147	>	branch->kid[1].gdv[0] = sd[0];
148	>	branch->kid[1].gdv[1] = dv[1];
149	>	branch->kid[1].gdv[2] = sd[1];
150	>	branch->kid[2].gdv[0] = sd[1];
151	>	branch->kid[2].gdv[1] = dv[2];
152	>	branch->kid[2].gdv[2] = sd[2];
153	>	branch->kid[3].gdv[0] = sd[0];
154	>	branch->kid[3].gdv[1] = sd[1];
155	>	branch->kid[3].gdv[2] = sd[2];
156		}
157
158		/* Recursively subdivide the given node to the specified quadtree depth */
#	Line 155 | Line 165 | branchsample(TRITREE *node, int depth)
165		return;
166		if (isleaf(node)) {                     /* subdivide leaf node */
167		TRITREE branch, *moved_leaf;
168	+	FVECT   sdir;
169		subdivide(&branch, node->gdv);
170	<	moved_leaf = findleaf(&branch, node->sd);
170	>	decodedir(sdir, node->sd);
171	>	moved_leaf = findleaf(&branch, sdir);
172		if (moved_leaf != NULL) {       /* bequeath old sample */
173	<	VCOPY(moved_leaf->sd, node->sd);
173	>	moved_leaf->sd = node->sd;
174		copycolr(moved_leaf->val, node->val);
175		}
176		for (i = 0; i < 4; i++)         /* compute new samples */
#	Line 172 | Line 184 | branchsample(TRITREE *node, int depth)
184
185		/* Sample sphere using triangular geodesic mesh */
186		TRITREE *
187	<	geosample(int nsamps)
187	>	geosample(long nsamps)
188		{
189		int     depth;
190		TRITREE *tree;
#	Line 181 | Line 193 | geosample(int nsamps)
193		/* figure out depth */
194		if ((nsamps -= 4) < 0)
195		error(USER, "minimum number of samples is 4");
196	<	nsamps = nsamps*3/NTRUNKBR;     /* round up */
196	>	nsamps = nsamps*(NTRUNKBR-1)/NTRUNKBR;  /* round up */
197		for (depth = 0; nsamps > 1; depth++)
198		nsamps >>= 2;
199		/* make base tetrahedron */
#	Line 193 | Line 205 | geosample(int nsamps)
205		trunk[1][1] = sqrt(1. - trunk[1][2]*trunk[1][2]);
206		spinvector(trunk[2], trunk[1], trunk[0], 2.*PI/3.);
207		spinvector(trunk[3], trunk[1], trunk[0], 4.*PI/3.);
208	<	VCOPY(tree->gdv[0], trunk[0]);
197	<	VCOPY(tree->gdv[1], trunk[0]);
198	<	VCOPY(tree->gdv[2], trunk[0]);
208	>	tree->gdv[0] = tree->gdv[1] = tree->gdv[2] = encodedir(trunk[0]);
209		tree->kid = (TRITREE *)calloc(NTRUNKBR, sizeof(TRITREE));
210		if (tree->kid == NULL) goto memerr;
211		/* grow our tree from trunk */
212		for (i = 0; i < NTRUNKBR; i++) {
213		for (j = 0; j < 3; j++)         /* XXX works for tetra only */
214	<	VCOPY(tree->kid[i].gdv[j], trunk[(i+j)%NTRUNKVERT]);
214	>	tree->kid[i].gdv[j] = encodedir(trunk[(i+j)%NTRUNKVERT]);
215		leafsample(&tree->kid[i]);
216		branchsample(&tree->kid[i], depth);
217		}
#	Line 229 | Line 239 | get_ehisto(const TRITREE *node, long exphisto[256])
239		double
240		get_threshold(const TRITREE *tree)
241		{
242	+	long    samptotal = 0;
243		long    exphisto[256];
233	–	long    samptotal;
244		int     i;
245		/* compute sample histogram */
246	<	memset((void *)exphisto, 0, sizeof(exphisto));
246	>	memset(exphisto, 0, sizeof(exphisto));
247		for (i = 0; i < NTRUNKBR; i++)
248		get_ehisto(&tree->kid[i], exphisto);
249		/* use 98th percentile */
#	Line 269 | Line 279 | findemax(TRITREE *node, int *expp)
279
280		/* Compute solid angle of spherical triangle (approx.) */
281		double
282	<	tri_omegav(FVECT v[3])
282	>	tri_omegav(const int32 vd[3])
283		{
284	<	FVECT   e1, e2, vcross;
285	<
284	>	FVECT   v[3], e1, e2, vcross;
285	>
286	>	decodedir(v[0], vd[0]);
287	>	decodedir(v[1], vd[1]);
288	>	decodedir(v[2], vd[2]);
289		VSUB(e1, v[1], v[0]);
290		VSUB(e2, v[2], v[1]);
291		fcross(vcross, e1, e2);
292		return(.5*VLEN(vcross));
293		}
294
295	<	/* Sum intensity times direction for non-zero leaves */
295	>	/* Sum intensity times direction for above-threshold perimiter within radius */
296		void
297		vector_sum(FVECT vsum, TRITREE *node,
298	<	const FVECT cent, double mincos, int ethresh)
298	>	FVECT cent, double maxr2, int ethresh)
299		{
300		if (isleaf(node)) {
301		double  intens;
302	+	FVECT   sdir;
303		if (node->val[EXP] < ethresh)
304	<	return;
305	<	if (DOT(node->sd,cent) < mincos)
306	<	return;
304	>	return;                         /* below threshold */
305	>	if (fdir2diff(node->sd,cent) > maxr2)
306	>	return;                         /* too far away */
307		intens = colrval(node->val,GRN) * tri_omegav(node->gdv);
308	<	VSUM(vsum, vsum, node->sd, intens);
308	>	decodedir(sdir, node->sd);
309	>	VSUM(vsum, vsum, sdir, intens);
310		return;
311		}
312	<	if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
313	<	DOT(node->gdv[0],cent) > mincos &&
314	<	DOT(node->gdv[1],cent) > mincos &&
315	<	DOT(node->gdv[2],cent) > mincos)
316	<	return;
317	<	vector_sum(vsum, &node->kid[0], cent, mincos, ethresh);
318	<	vector_sum(vsum, &node->kid[1], cent, mincos, ethresh);
319	<	vector_sum(vsum, &node->kid[2], cent, mincos, ethresh);
320	<	vector_sum(vsum, &node->kid[3], cent, mincos, ethresh);
312	>	if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
313	>	fdir2diff(node->gdv[0],cent) < maxr2 &&
314	>	fdir2diff(node->gdv[1],cent) < maxr2 &&
315	>	fdir2diff(node->gdv[2],cent) < maxr2)
316	>	return;                                 /* containing node */
317	>	vector_sum(vsum, &node->kid[0], cent, maxr2, ethresh);
318	>	vector_sum(vsum, &node->kid[1], cent, maxr2, ethresh);
319	>	vector_sum(vsum, &node->kid[2], cent, maxr2, ethresh);
320	>	vector_sum(vsum, &node->kid[3], cent, maxr2, ethresh);
321		}
322
323		/* Claim source contributions within the given solid angle */
324		void
325	<	claimlight(COLOR intens, TRITREE *node, const FVECT cent, double mincos)
325	>	claimlight(COLOR intens, TRITREE *node, FVECT cent, double maxr2)
326		{
327		int     remaining;
328		int     i;
329		if (isleaf(node)) {             /* claim contribution */
330		COLOR   contrib;
331		if (node->val[EXP] <= 0)
332	<	return;
333	<	if (DOT(node->sd,cent) < mincos)
334	<	return;
332	>	return;         /* already claimed */
333	>	if (fdir2diff(node->sd,cent) > maxr2)
334	>	return;         /* too far away */
335		colr_color(contrib, node->val);
336		scalecolor(contrib, tri_omegav(node->gdv));
337		addcolor(intens, contrib);
338		copycolr(node->val, blkclr);
339		return;
340		}
341	<	if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
342	<	DOT(node->gdv[0],cent) > mincos &&
343	<	DOT(node->gdv[1],cent) > mincos &&
344	<	DOT(node->gdv[2],cent) > mincos)
345	<	return;
341	>	if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
342	>	fdir2diff(node->gdv[0],cent) < maxr2 &&
343	>	fdir2diff(node->gdv[1],cent) < maxr2 &&
344	>	fdir2diff(node->gdv[2],cent) < maxr2)
345	>	return;                 /* previously claimed node */
346		remaining = 0;                  /* recurse on children */
347		for (i = 0; i < 4; i++) {
348	<	claimlight(intens, &node->kid[i], cent, mincos);
348	>	claimlight(intens, &node->kid[i], cent, maxr2);
349		if (!isleaf(&node->kid[i]) || node->kid[i].val[EXP] != 0)
350		++remaining;
351		}
352		if (remaining)
353		return;
354		/* consolidate empties */
355	<	free((void *)node->kid); node->kid = NULL;
355	>	free(node->kid); node->kid = NULL;
356		copycolr(node->val, blkclr);
357	<	VCOPY(node->sd, node->gdv[0]);  /* doesn't really matter */
357	>	node->sd = node->gdv[0];        /* doesn't really matter */
358		}
359
360		/* Add lost light contribution to the given list */
361		void
362	<	add2lost(LOSTLIGHT **llp, COLOR intens, const FVECT cent)
362	>	add2lost(LOSTLIGHT **llp, COLOR intens, FVECT cent)
363		{
364		LOSTLIGHT       *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
365
366		if (newll == NULL)
367		return;
368		copycolor(newll->intens, intens);
369	<	VCOPY(newll->sd, cent);
369	>	newll->sd = encodedir(cent);
370		newll->next = *llp;
371		*llp = newll;
372		}
373
374		/* Check lost light list for contributions */
375		void
376	<	getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cent, double omega)
376	>	getlost(LOSTLIGHT **llp, COLOR intens, FVECT cent, double omega)
377		{
378	<	const double    mincos = 1. - omega/(2.*PI);
378	>	const double    maxr2 = omega/PI;
379		LOSTLIGHT       lhead, *lastp, *thisp;
380
381		lhead.next = *llp;
382		lastp = &lhead;
383		while ((thisp = lastp->next) != NULL)
384	<	if (DOT(thisp->sd,cent) >= mincos) {
384	>	if (fdir2diff(thisp->sd,cent) <= maxr2) {
385		LOSTLIGHT       *mynext = thisp->next;
386		addcolor(intens, thisp->intens);
387	<	free((void *)thisp);
387	>	free(thisp);
388		lastp->next = mynext;
389		} else
390		lastp = thisp;
#	Line 380 | Line 395 | getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cen
395		void
396		mksources(TRITREE *samptree, double thresh, double maxang)
397		{
398	+	#define MAXITER         100
399		const int       ethresh = (int)(log(thresh)/log(2.) + (COLXS+.5));
400		const double    maxomega = 2.*PI*(1. - cos(PI/180./2.*maxang));
401		const double    minintens = .05*thresh*maxomega;
402	+	int             niter = MAXITER;
403		int             nsrcs = 0;
404		LOSTLIGHT       *lostlightlist = NULL;
405		int             emax;
#	Line 408 | Line 425 | mksources(TRITREE *samptree, double thresh, double max
425		*/
426		if (thresh <= FTINY)
427		return;
428	<	for ( ; ; ) {
428	>	while (niter--) {
429		emax = ethresh;         /* find brightest unclaimed */
430		startleaf = NULL;
431		for (i = 0; i < NTRUNKBR; i++) {
#	Line 419 | Line 436 | mksources(TRITREE *samptree, double thresh, double max
436		if (startleaf == NULL)
437		break;
438		/* claim it */
439	<	VCOPY(curcent, startleaf->sd);
439	>	decodedir(curcent, startleaf->sd);
440		curomega = tri_omegav(startleaf->gdv);
441		currad = sqrt(curomega/PI);
442		growstep = 3.*currad;
#	Line 435 | Line 452 | mksources(TRITREE *samptree, double thresh, double max
452		vsum[0] = vsum[1] = vsum[2] = .0;
453		for (i = 0; i < NTRUNKBR; i++)
454		vector_sum(vsum, &samptree->kid[i],
455	<	curcent, cos(currad+growstep),
455	>	curcent, 2.-2.*cos(currad+growstep),
456		thisethresh);
457		if (normalize(vsum) == .0)
458		break;
459	<	movedist = acos(DOT(vsum,curcent));
459	>	movedist = Acos(DOT(vsum,curcent));
460		if (movedist > growstep) {
461		VSUB(vsum, vsum, curcent);
462		movedist = growstep/VLEN(vsum);
#	Line 451 | Line 468 | mksources(TRITREE *samptree, double thresh, double max
468		curomega = 2.*PI*(1. - cos(currad));
469		for (i = 0; i < NTRUNKBR; i++)
470		claimlight(curintens, &samptree->kid[i],
471	<	curcent, cos(currad));
471	>	curcent, 2.-2.*cos(currad));
472		} while (curomega < maxomega &&
473		bright(curintens)/curomega > thisthresh);
474		if (bright(curintens) < minintens) {
#	Line 471 | Line 488 | mksources(TRITREE *samptree, double thresh, double max
488		printf("0\n0\n4 %f %f %f %f\n",
489		curcent[0], curcent[1], curcent[2],
490		2.*180./PI*currad);
491	+	niter = MAXITER;
492		}
493	+	#undef MAXITER
494		}
495
496		int
#	Line 482 | Line 501 | main(int argc, char *argv[])
501		TRITREE *samptree;
502		double  thresh = 0;
503		int     i;
504	<
505	<	progname = argv[0];
504	>
505	>	fixargv0(argv[0]);              /* sets global progname */
506	>
507		for (i = 1; i < argc && argv[i][0] == '-'; i++)
508		switch (argv[i][1]) {
509		case 'd':               /* number of samples */
#	Line 536 | Line 556 | userr:
556		}
557
558		void
559	<	eputs(char  *s)
559	>	eputs(const char  *s)
560		{
561		static int  midline = 0;
562
#	Line 554 | Line 574 | eputs(char  *s)
574		}
575
576		void
577	<	wputs(char *s)
577	>	wputs(const char *s)
578		{
579		/* no warnings */
580		}

Diff Legend

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