ViewVC Help

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

root/Development/ray/src/util/rsensor.c

Comparing ray/src/util/rsensor.c (file contents):
Revision 2.6 by greg, Sat Dec 13 00:44:05 2008 UTC vs.
Revision 2.24 by greg, Wed Apr 23 02:35:27 2025 UTC

#	Line 9 | Line 9 | static const char RCSid[] = "$Id$";
9		*/
10
11		#include "ray.h"
12	+	#include "platform.h"
13		#include "source.h"
14	+	#include "func.h"
15		#include "view.h"
16		#include "random.h"
17
18		#define DEGREE          (PI/180.)
19
20	<	#define MAXNT           180     /* maximum number of theta divisions */
20	>	#define MAXNT           181     /* maximum number of theta divisions */
21		#define MAXNP           360     /* maximum number of phi divisions */
22
23		extern char     *progname;      /* global argv[0] */
#	Line 26 | Line 28 | VIEW           ourview =  {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,0.,
28		1.,180.,180.,0.,0.,0.,0.,
29		{0.,0.,0.},{0.,0.,0.},0.,0.};
30
31	<	unsigned long   nsamps = 10000; /* desired number of initial samples */
30	<	unsigned long   nssamps = 9000; /* number of super-samples */
31	>	long            nsamps = 10000; /* desired number of initial samples */
32		int             ndsamps = 32;   /* number of direct samples */
33		int             nprocs = 1;     /* number of rendering processes */
34
35		float           *sensor = NULL; /* current sensor data */
36		int             sntp[2];        /* number of sensor theta and phi angles */
37		float           maxtheta;       /* maximum theta value for this sensor */
38	<	float           tvals[MAXNT+1]; /* theta values (1-D table of 1-cos(t)) */
39	<	float           *pvals = NULL;  /* phi values (2-D table in radians) */
38	>	float           tvals[MAXNT+1]; /* theta prob. values (1-D table of 1-cos(t)) */
39	>	float           *pvals = NULL;  /* phi prob. values (2-D table in radians) */
40		int             ntheta = 0;     /* polar angle divisions */
41		int             nphi = 0;       /* azimuthal angle divisions */
42		double          gscale = 1.;    /* global scaling value */
#	Line 59 | Line 60 | print_defaults()               /* print out default parameters */
60		over_options();
61		printf("-n %-9d\t\t\t# number of processes\n", nprocs);
62		printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
62	–	/* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
63		printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
64		printf("-vp %f %f %f\t# view point\n",
65		ourview.vp[0], ourview.vp[1], ourview.vp[2]);
#	Line 71 | Line 71 | print_defaults()               /* print out default parameters */
71		print_rdefaults();
72		}
73
74	+
75	+	void
76	+	quit(ec)                        /* make sure exit is called */
77	+	int     ec;
78	+	{
79	+	if (ray_pnprocs > 0)    /* close children if any */
80	+	ray_pclose(0);
81	+	else if (ray_pnprocs < 0)
82	+	_exit(ec);      /* avoid flush in child */
83	+	exit(ec);
84	+	}
85	+
86	+
87		int
88		main(
89		int     argc,
#	Line 89 | Line 102 | main(
102		directrelay = 3;
103		ambounce = 1;
104		maxdepth = -10;
105	+	/* initialize calcomp routines */
106	+	initfunc();
107		/* get options from command line */
108		for (i = 1; i < argc; i++) {
109		while ((rval = expandarg(&argc, &argv, i)) > 0)
#	Line 109 | Line 124 | main(
124		putchar('\n');
125		}
126		/* start process(es) */
127	<	ray_pinit(argv[argc-1], nprocs);
127	>	if (strcmp(argv[argc-1], "."))
128	>	ray_pinit(argv[argc-1], nprocs);
129		}
130		comp_sensor(argv[i]);   /* process a sensor file */
131		continue;
#	Line 117 | Line 133 | main(
133		if (argv[i][1] == 'r') {        /* sampling options */
134		if (argv[i][2] == 'd')
135		nsamps = atol(argv[++i]);
120	–	else if (argv[i][2] == 's')
121	–	nssamps = atol(argv[++i]);
136		else {
137		sprintf(errmsg, "bad option at '%s'", argv[i]);
138		error(USER, errmsg);
#	Line 183 | Line 197 | main(
197		}
198		i += rval;
199		}
200	<	if (!ray_pnprocs)
200	>	if (sensor == NULL)
201		error(USER, i<argc ? "missing sensor file" : "missing octree");
202		quit(0);
203		}
#	Line 195 | Line 209 | load_sensor(
209		char    *sfile
210		)
211		{
212	+	int     warnedneg;
213		char    linebuf[8192];
214	+	int     last_pos_val = 0;
215		int     nelem = 1000;
216		float   *sarr = (float *)malloc(sizeof(float)*nelem);
217		FILE    *fp;
#	Line 223 | Line 239 | load_sensor(
239		cp = fskip(cp);
240		if (cp == NULL)
241		break;
242	+	if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]+FTINY) {
243	+	sprintf(errmsg,
244	+	"Phi values not monotinically increasing in sensor file '%s'",
245	+	sfile);
246	+	error(USER, errmsg);
247	+	}
248		++ntp[1];
249		}
250	+	warnedneg = 0;
251		ntp[0] = 0;                             /* get thetas + data */
252		while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
253		++ntp[0];
#	Line 242 | Line 265 | load_sensor(
265		cp = fskip(cp);
266		if (cp == NULL)
267		break;
268	+	if (sarr[i] < .0) {
269	+	if (!warnedneg++) {
270	+	sprintf(errmsg,
271	+	"Negative value(s) in sensor file '%s' (ignored)\n", sfile);
272	+	error(WARNING, errmsg);
273	+	}
274	+	sarr[i] = .0;
275	+	} else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1))
276	+	last_pos_val = i;
277		++i;
278		}
279	<	if (i == ntp[0]*(ntp[1]+1))
279	>	if (i == ntp[0]*(ntp[1]+1))     /* empty line? */
280		break;
281	+	if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
282	+	sarr[(ntp[0]-1)*(ntp[1]+1)]) {
283	+	sprintf(errmsg,
284	+	"Theta values not monotinically increasing in sensor file '%s'",
285	+	sfile);
286	+	error(USER, errmsg);
287	+	}
288		if (i != (ntp[0]+1)*(ntp[1]+1)) {
289		sprintf(errmsg,
290		"bad column count near line %d in sensor file '%s'",
#	Line 253 | Line 292 | load_sensor(
292		error(USER, errmsg);
293		}
294		}
295	<	nelem = i;
295	>	/* truncate zero region */
296	>	ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1;
297	>	nelem = (ntp[0]+1)*(ntp[1]+1);
298		fclose(fp);
299		errmsg[0] = '\0';                       /* sanity checks */
300	<	if (ntp[0] <= 0)
301	<	sprintf(errmsg, "no data in sensor file '%s'", sfile);
300	>	if (!last_pos_val)
301	>	sprintf(errmsg, "no positive sensor values in file '%s'", sfile);
302		else if (fabs(sarr[ntp[1]+1]) > FTINY)
303		sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
304		sfile);
305		else if (fabs(sarr[1]) > FTINY)
306		sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
307		sfile);
308	<	else if (sarr[ntp[1]] <= FTINY)
308	>	else if (sarr[ntp[1]] < 270.-FTINY)
309		sprintf(errmsg,
310	<	"maximum phi must be positive in sensor file '%s'",
310	>	"maximum phi must be 270 or greater in sensor file '%s'",
311		sfile);
312	<	else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
312	>	else if (sarr[ntp[1]] >= 360.-FTINY)
313		sprintf(errmsg,
314	<	"maximum theta must be positive in sensor file '%s'",
314	>	"maximum phi must be less than 360 in sensor file '%s'",
315		sfile);
316		if (errmsg[0])
317		error(USER, errmsg);
#	Line 283 | Line 324 | init_ptable(
324		char    *sfile
325		)
326		{
327	<	int     samptot = nsamps;
327	>	long    samptot = nsamps;
328		float   *rowp, *rowp1;
329		double  rowsum[MAXNT], rowomega[MAXNT];
330		double  thdiv[MAXNT+1], phdiv[MAXNP+1];
#	Line 315 | Line 356 | init_ptable(
356		error(INTERNAL, errmsg);
357		}
358		/* compute boundary angles */
359	<	maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
359	>	maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
360	>	if (maxtheta > PI)
361	>	maxtheta = PI;
362		thdiv[0] = .0;
363		for (t = 1; t < sntp[0]; t++)
364		thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
365	<	thdiv[sntp[0]] = maxtheta*DEGREE;
366	<	phdiv[0] = .0;
365	>	thdiv[sntp[0]] = maxtheta;
366	>	phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
367		for (p = 1; p < sntp[1]; p++)
368		phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
369	<	phdiv[sntp[1]] = 2.*PI;
369	>	phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
370		/* size our table */
371	<	tsize = 1. - cos(maxtheta*DEGREE);
372	<	psize = PI*tsize/(maxtheta*DEGREE);
371	>	tsize = 1. - cos(maxtheta);
372	>	psize = PI*tsize/maxtheta;
373		if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
374		samptot = sntp[0]*sntp[1];
375	<	ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
375	>	ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1;
376		if (ntheta > MAXNT)
377		ntheta = MAXNT;
378		nphi = samptot/ntheta;
379	<	pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
379	>	pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
380		if (pvals == NULL)
381		error(SYSTEM, "out of memory in init_ptable()");
382		gscale = .0;                    /* compute our inverse table */
383		for (i = 0; i < sntp[0]; i++) {
384		rowp = &s_val(i,0);
385	<	rowsum[i] = 0.;
385	>	rowsum[i] = 1e-20;
386		for (j = 0; j < sntp[1]; j++)
387		rowsum[i] += *rowp++;
388		rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
389		rowomega[i] *= 2.*PI / (double)sntp[1];
390		gscale += rowsum[i] * rowomega[i];
391		}
392	+	if (gscale <= FTINY) {
393	+	sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
394	+	error(USER, errmsg);
395	+	}
396		for (i = 0; i < ntheta; i++) {
397		prob = (double)i / (double)ntheta;
398		for (t = 0; t < sntp[0]; t++)
#	Line 357 | Line 404 | init_ptable(
404		tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
405		frac*cos(thdiv[t+1]) );
406		/* offset b/c sensor values are centered */
407	<	if (t <= 0 || frac > 0.5)
407	>	if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) {
408		frac -= 0.5;
409	<	else if (t >= sntp[0]-1 || frac < 0.5) {
409	>	} else {
410		frac += 0.5;
411		--t;
412		}
413	<	pvals[i*(nphi+1)] = .0f;
413	>	pvals[i*(nphi+1)] = phdiv[0];
414		for (j = 1; j < nphi; j++) {
415		prob = (double)j / (double)nphi;
416		rowp = &s_val(t,0);
417		rowp1 = &s_val(t+1,0);
418	<	for (p = 0; p < sntp[1]; p++) {
418	>	for (p = 0; p < sntp[1]; p++)
419		if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
420		frac*rowp1[p]/rowsum[t+1]) <= .0)
421		break;
422	<	if (p >= sntp[1])
423	<	error(INTERNAL,
424	<	"code error 2 in init_ptable()");
378	<	frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
379	<	+ frac*rowp1[p]/rowsum[t+1]);
380	<	pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
381	<	frac1*phdiv[p+1];
422	>	if (p >= sntp[1]) {     /* should never happen? */
423	>	p = sntp[1] - 1;
424	>	prob = .5;
425		}
426	+	frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
427	+	+ frac*rowp1[p]/rowsum[t+1]);
428	+	pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
429	+	frac1*phdiv[p+1];
430		}
431	<	pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
431	>	pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
432		}
433	+	/* duplicate final row */
434	+	memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1),
435	+	sizeof(*pvals)*(nphi+1));
436		tvals[0] = .0f;
437		tvals[ntheta] = (float)tsize;
438		}
#	Line 429 | Line 479 | sens_val(
479		int     t, p;
480
481		dv[2] = DOT(dvec, ourview.vdir);
482	<	theta = (float)((1./DEGREE) * acos(dv[2]));
482	>	theta = acos(dv[2]);
483		if (theta >= maxtheta)
484		return(.0f);
485		dv[0] = DOT(dvec, ourview.hvec);
486		dv[1] = DOT(dvec, ourview.vvec);
487	<	phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
488	<	while (phi < .0f) phi += 360.f;
487	>	phi = atan2(-dv[0], dv[1]);
488	>	while (phi < .0f) phi += (float)(2.*PI);
489		t = (int)(theta/maxtheta * sntp[0]);
490	<	p = (int)(phi*(1./360.) * sntp[1]);
490	>	p = (int)(phi*(1./(2.*PI)) * sntp[1]);
491		/* hack for non-uniform sensor grid */
492	+	theta *= (float)(1./DEGREE);
493	+	phi *= (float)(1./DEGREE);
494		while (t+1 < sntp[0] && theta >= s_theta(t+1))
495		++t;
496		while (t-1 >= 0 && theta <= s_theta(t-1))
#	Line 450 | Line 502 | sens_val(
502		return(s_val(t,p));
503		}
504
505	+	/* Print origin and direction */
506	+	static void
507	+	print_ray(
508	+	FVECT rorg,
509	+	FVECT rdir
510	+	)
511	+	{
512	+	printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
513	+	rorg[0], rorg[1], rorg[2],
514	+	rdir[0], rdir[1], rdir[2]);
515	+	}
516	+
517		/* Compute sensor output */
518		static void
519		comp_sensor(
#	Line 460 | Line 524 | comp_sensor(
524		ndsamps > 0 ? 1 : 0;
525		char    *err;
526		int     nt, np;
527	<	COLOR   vsum;
527	>	SCOLOR  vsum;
528		RAY     rr;
529		double  sf;
530		int     i, j;
#	Line 474 | Line 538 | comp_sensor(
538		/* assign probability table */
539		init_ptable(sfile);
540		/* stratified MC sampling */
541	<	setcolor(vsum, .0f, .0f, .0f);
541	>	scolorblack(vsum);
542		nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
543		np = nsamps/nt;
544		sf = gscale/nsamps;
#	Line 484 | Line 548 | comp_sensor(
548		get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
549		if (ourview.vfore > FTINY)
550		VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
551	+	if (!ray_pnprocs) {
552	+	print_ray(rr.rorg, rr.rdir);
553	+	continue;
554	+	}
555		rr.rmax = .0;
556	<	rayorigin(&rr, PRIMARY, NULL, NULL);
557	<	scalecolor(rr.rcoef, sf);
556	>	rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
557	>	scalescolor(rr.rcoef, sf);
558		if (ray_pqueue(&rr) == 1)
559	<	addcolor(vsum, rr.rcol);
559	>	saddscolor(vsum, rr.rcol);
560		}
561		/* remaining rays pure MC */
562		for (i = nsamps - nt*np; i-- > 0; ) {
#	Line 496 | Line 564 | comp_sensor(
564		get_direc(rr.rdir, frandom(), frandom());
565		if (ourview.vfore > FTINY)
566		VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
567	+	if (!ray_pnprocs) {
568	+	print_ray(rr.rorg, rr.rdir);
569	+	continue;
570	+	}
571		rr.rmax = .0;
572	<	rayorigin(&rr, PRIMARY, NULL, NULL);
573	<	scalecolor(rr.rcoef, sf);
572	>	rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
573	>	scalescolor(rr.rcoef, sf);
574		if (ray_pqueue(&rr) == 1)
575	<	addcolor(vsum, rr.rcol);
575	>	saddscolor(vsum, rr.rcol);
576		}
577	+	if (!ray_pnprocs)                       /* just printing rays */
578	+	return;
579		/* scale partial result */
580	<	scalecolor(vsum, sf);
580	>	scalescolor(vsum, sf);
581		/* add direct component */
582		for (i = ndirs; i-- > 0; ) {
583		SRCINDEX        si;
#	Line 513 | Line 587 | comp_sensor(
587		if (sf <= FTINY)
588		continue;
589		sf *= si.dom/ndirs;
590	<	scalecolor(rr.rcoef, sf);
590	>	scalescolor(rr.rcoef, sf);
591		if (ray_pqueue(&rr) == 1) {
592	<	multcolor(rr.rcol, rr.rcoef);
593	<	addcolor(vsum, rr.rcol);
592	>	smultscolor(rr.rcol, rr.rcoef);
593	>	saddscolor(vsum, rr.rcol);
594		}
595		}
596		}
597	<	/* finish our calculation */
598	<	while (ray_presult(&rr, 0) > 0) {
599	<	multcolor(rr.rcol, rr.rcoef);
526	<	addcolor(vsum, rr.rcol);
597	>	while (ray_presult(&rr, 0) > 0) {       /* finish our calculation */
598	>	smultscolor(rr.rcol, rr.rcoef);
599	>	saddscolor(vsum, rr.rcol);
600		}
601	<	/* print our result */
602	<	printf("%.4e %.4e %.4e\n", colval(vsum,RED),
603	<	colval(vsum,GRN), colval(vsum,BLU));
601	>	for (i = 0; i < NCSAMP; i++)            /* print our result */
602	>	printf(" %.4e", vsum[i]);
603	>	fputc('\n', stdout);
604		}

Diff Legend

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