ViewVC Help

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

root/Development/ray/src/cv/bsdfrep.c

Comparing ray/src/cv/bsdfrep.c (file contents):
Revision 2.21 by greg, Fri Mar 7 21:21:02 2014 UTC vs.
Revision 2.38 by greg, Sun May 18 01:46:05 2025 UTC

#	Line 9 | Line 9 | static const char RCSid[] = "$Id$";
9
10		#define _USE_MATH_DEFINES
11		#include <stdlib.h>
12	–	#include <string.h>
12		#include <math.h>
13		#include "rtio.h"
14		#include "resolu.h"
15		#include "bsdfrep.h"
16	+	#include "random.h"
17		/* name and manufacturer if known */
18		char                    bsdf_name[256];
19		char                    bsdf_manuf[256];
#	Line 29 | Line 29 | int                    single_plane_incident = -1;
29		int                     input_orient = 0;
30		int                     output_orient = 0;
31
32	+	/* represented color space */
33	+	RBColor                 rbf_colorimetry = RBCunknown;
34	+
35	+	const char              *RBCident[] = {
36	+	"CIE-Y", "CIE-XYZ", "Spectral", "Unknown"
37	+	};
38	+
39		/* BSDF histogram */
40		unsigned long           bsdf_hist[HISTLEN];
41
42		/* BSDF value for boundary regions */
43		double                  bsdf_min = 0;
44	+	double                  bsdf_spec_val = 0;
45	+	double                  bsdf_spec_rad = 0;
46
47		/* processed incident DSF measurements */
48		RBFNODE                 *dsf_list = NULL;
#	Line 44 | Line 53 | MIGRATION              *mig_list = NULL;
53		/* current input direction */
54		double                  theta_in_deg, phi_in_deg;
55
56	+	/* header line sharing callback */
57	+	int                     (*sir_headshare)(char *s) = NULL;
58	+
59		/* Register new input direction */
60		int
61		new_input_direction(double new_theta, double new_phi)
62		{
51	–	if (!input_orient)              /* check input orientation */
52	–	input_orient = 1 - 2*(new_theta > 90.);
53	–	else if (input_orient > 0 ^ new_theta < 90.) {
54	–	fprintf(stderr,
55	–	"%s: Cannot handle input angles on both sides of surface\n",
56	–	progname);
57	–	return(0);
58	–	}
63		/* normalize angle ranges */
64		while (new_theta < -180.)
65		new_theta += 360.;
#	Line 65 | Line 69 | new_input_direction(double new_theta, double new_phi)
69		new_theta = -new_theta;
70		new_phi += 180.;
71		}
68	–	if ((theta_in_deg = new_theta) < 1.0)
69	–	return(1);              /* don't rely on phi near normal */
72		while (new_phi < 0)
73		new_phi += 360.;
74		while (new_phi >= 360.)
75		new_phi -= 360.;
76	+	/* check input orientation */
77	+	if (!input_orient)
78	+	input_orient = 1 - 2*(new_theta > 90.);
79	+	else if (input_orient > 0 ^ new_theta < 90.) {
80	+	fprintf(stderr,
81	+	"%s: Cannot handle input angles on both sides of surface\n",
82	+	progname);
83	+	return(0);
84	+	}
85	+	if ((theta_in_deg = new_theta) < 1.0)
86	+	return(1);              /* don't rely on phi near normal */
87		if (single_plane_incident > 0)  /* check input coverage */
88		single_plane_incident = (round(new_phi) == round(phi_in_deg));
89		else if (single_plane_incident < 0)
#	Line 91 | Line 104 | new_input_direction(double new_theta, double new_phi)
104		int
105		use_symmetry(FVECT vec)
106		{
107	<	const double    phi = get_phi360(vec);
107	>	double  phi = get_phi360(vec);
108	>	/* because of -0. issue */
109	>	while (phi >= 360.) phi -= 360.;
110	>	while (phi < 0.) phi += 360.;
111
112		switch (inp_coverage) {
113		case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
#	Line 198 | Line 214 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
214		int                     pos[2];
215		int                     n;
216
217	<	for (n = ((-.01 > phi) | (phi > .01))*rbf->nrbf; n-- > 0; ) {
217	>	for (n = (cos(phi) < 1.-FTINY)*rbf->nrbf; n-- > 0; ) {
218		ovec_from_pos(outvec, rbf->rbfa[n].gx, rbf->rbfa[n].gy);
219		spinvector(outvec, outvec, vnorm, phi);
220		pos_from_vec(pos, outvec);
#	Line 209 | Line 225 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
225		}
226
227		/* Compute outgoing vector from grid position */
228	+	#if 1
229		void
230		ovec_from_pos(FVECT vec, int xpos, int ypos)
231	<	{
232	<	double  uv[2];
231	>	{                               /* precomputed table version */
232	>	static int      qsiz = 0;
233	>	static float    (*q_uv)[2] = NULL;
234	>
235	>	if (vec == NULL) {      /* just free table? */
236	>	if (q_uv) free(q_uv);
237	>	qsiz = 0;
238	>	return;
239	>	}
240	>	if (qsiz != grid_res>>1) {
241	>	int     x, y;   /* (re)make positive quadrant table */
242	>	RREAL   uv[2];
243	>	double  r;
244	>	if (q_uv) free(q_uv);
245	>	qsiz = grid_res>>1;
246	>	q_uv = (float (*)[2])malloc(sizeof(float)*2*qsiz*qsiz);
247	>	for (y = qsiz; y--; )
248	>	for (x = qsiz; x--; ) {
249	>	square2disk(uv, 0.5 + (x+.5)/grid_res,
250	>	0.5 + (y+.5)/grid_res);
251	>	/* uniform hemispherical projection */
252	>	r = sqrt(2. - uv[0]*uv[0] - uv[1]*uv[1]);
253	>	q_uv[qsiz*y + x][0] = (float)(r*uv[0]);
254	>	q_uv[qsiz*y + x][1] = (float)(r*uv[1]);
255	>	}
256	>	}
257	>	/* put in positive quadrant */
258	>	if (xpos >= qsiz) { xpos -= qsiz; vec[0] = 1.; }
259	>	else { xpos = qsiz-1 - xpos; vec[0] = -1.; }
260	>	if (ypos >= qsiz) { ypos -= qsiz; vec[1] = 1.; }
261	>	else { ypos = qsiz-1 - ypos; vec[1] = -1.; }
262	>
263	>	vec[0] *= (RREAL)q_uv[qsiz*ypos + xpos][0];
264	>	vec[1] *= (RREAL)q_uv[qsiz*ypos + xpos][1];
265	>	vec[2] = output_orient*sqrt(1. - vec[0]*vec[0] - vec[1]*vec[1]);
266	>	}
267	>	#else
268	>	void
269	>	ovec_from_pos(FVECT vec, int xpos, int ypos)
270	>	{                               /* table-free version */
271	>	RREAL   uv[2];
272		double  r2;
273	<
274	<	SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
273	>
274	>	if (vec == NULL)
275	>	return;
276	>
277	>	square2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
278		/* uniform hemispherical projection */
279		r2 = uv[0]*uv[0] + uv[1]*uv[1];
280		vec[0] = vec[1] = sqrt(2. - r2);
#	Line 223 | Line 282 | ovec_from_pos(FVECT vec, int xpos, int ypos)
282		vec[1] *= uv[1];
283		vec[2] = output_orient*(1. - r2);
284		}
285	+	#endif
286
287		/* Compute grid position from normalized input/output vector */
288		void
289		pos_from_vec(int pos[2], const FVECT vec)
290		{
291	<	double  sq[2];          /* uniform hemispherical projection */
291	>	RREAL   sq[2];          /* uniform hemispherical projection */
292		double  norm = 1./sqrt(1. + fabs(vec[2]));
293
294	<	SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
294	>	disk2square(sq, vec[0]*norm, vec[1]*norm);
295
296		pos[0] = (int)(sq[0]*grid_res);
297		pos[1] = (int)(sq[1]*grid_res);
#	Line 257 | Line 317 | rbf_volume(const RBFVAL *rbfp)
317		return(integ);
318		}
319
320	<	/* Evaluate RBF for DSF at the given normalized outgoing direction */
321	<	double
322	<	eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
320	>	/* Evaluate BSDF at the given normalized outgoing direction in color */
321	>	SDError
322	>	eval_rbfcol(SDValue *sv, const RBFNODE *rp, const FVECT outvec)
323		{
324		const double    rfact2 = (38./M_PI/M_PI)*(grid_res*grid_res);
265	–	double          minval = bsdf_min*output_orient*outvec[2];
325		int             pos[2];
326		double          res = 0;
327	+	double          usum = 0, vsum = 0;
328		const RBFVAL    *rbfp;
329		FVECT           odir;
330		double          rad2;
331		int             n;
332	+	/* assign default value */
333	+	sv->spec = c_dfcolor;
334	+	sv->cieY = bsdf_min;
335		/* check for wrong side */
336	<	if (outvec[2] > 0 ^ output_orient > 0)
337	<	return(.0);
338	<	/* use minimum if no information avail. */
339	<	if (rp == NULL)
340	<	return(minval);
336	>	if (outvec[2] > 0 ^ output_orient > 0) {
337	>	strcpy(SDerrorDetail, "Wrong-side scattering query");
338	>	return(SDEargument);
339	>	}
340	>	if (rp == NULL)         /* return minimum if no information avail. */
341	>	return(SDEnone);
342		/* optimization for fast lobe culling */
343		pos_from_vec(pos, outvec);
344		/* sum radial basis function */
#	Line 282 | Line 346 | eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
346		for (n = rp->nrbf; n--; rbfp++) {
347		int     d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
348		(pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
349	+	double  val;
350		rad2 = R2ANG(rbfp->crad);
351		rad2 *= rad2;
352		if (d2 > rad2*rfact2)
353		continue;
354		ovec_from_pos(odir, rbfp->gx, rbfp->gy);
355	<	res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
355	>	val = rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
356	>	if (rbf_colorimetry == RBCtristimulus) {
357	>	usum += val * (rbfp->chroma & 0xff);
358	>	vsum += val * (rbfp->chroma>>8 & 0xff);
359	>	}
360	>	res += val;
361		}
362	<	if (res < minval)       /* never return less than minval */
363	<	return(minval);
364	<	return(res);
362	>	sv->cieY = res / COSF(outvec[2]);
363	>	if (sv->cieY < bsdf_min) {      /* never return less than bsdf_min */
364	>	sv->cieY = bsdf_min;
365	>	} else if (rbf_colorimetry == RBCtristimulus) {
366	>	C_CHROMA        cres = (int)(usum/res + frandom());
367	>	cres |= (int)(vsum/res + frandom()) << 8;
368	>	c_decodeChroma(&sv->spec, cres);
369	>	}
370	>	return(SDEnone);
371		}
372
373	+	/* Evaluate BSDF at the given normalized outgoing direction in Y */
374	+	double
375	+	eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
376	+	{
377	+	SDValue sv;
378	+
379	+	if (eval_rbfcol(&sv, rp, outvec) == SDEnone)
380	+	return(sv.cieY);
381	+
382	+	return(0.0);
383	+	}
384	+
385		/* Insert a new directional scattering function in our global list */
386		int
387		insert_dsf(RBFNODE *newrbf)
#	Line 304 | Line 392 | insert_dsf(RBFNODE *newrbf)
392		for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
393		if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
394		fprintf(stderr,
395	<	"%s: Duplicate incident measurement (ignored)\n",
396	<	progname);
395	>	"%s: Duplicate incident measurement ignored at (%.1f,%.1f)\n",
396	>	progname, get_theta180(newrbf->invec),
397	>	get_phi360(newrbf->invec));
398		free(newrbf);
399		return(-1);
400		}
#	Line 393 | Line 482 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
482		return((rbfv[0] != NULL) + (rbfv[1] != NULL));
483		}
484
485	+	/* Return single-lobe specular RBF for the given incident direction */
486	+	RBFNODE *
487	+	def_rbf_spec(const FVECT invec)
488	+	{
489	+	RBFNODE         *rbf;
490	+	FVECT           ovec;
491	+	int             pos[2];
492	+
493	+	if (input_orient > 0 ^ invec[2] > 0)    /* wrong side? */
494	+	return(NULL);
495	+	if ((bsdf_spec_val <= bsdf_min) | (bsdf_spec_rad <= 0))
496	+	return(NULL);                   /* nothing set */
497	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE));
498	+	if (rbf == NULL)
499	+	return(NULL);
500	+	ovec[0] = -invec[0];
501	+	ovec[1] = -invec[1];
502	+	ovec[2] = invec[2]*(2*(input_orient==output_orient) - 1);
503	+	pos_from_vec(pos, ovec);
504	+	rbf->ord = 0;
505	+	rbf->next = NULL;
506	+	rbf->ejl = NULL;
507	+	VCOPY(rbf->invec, invec);
508	+	rbf->nrbf = 1;
509	+	rbf->rbfa[0].peak = bsdf_spec_val * COSF(ovec[2]);
510	+	rbf->rbfa[0].chroma = c_dfchroma;
511	+	rbf->rbfa[0].crad = ANG2R(bsdf_spec_rad);
512	+	rbf->rbfa[0].gx = pos[0];
513	+	rbf->rbfa[0].gy = pos[1];
514	+	rbf->vtotal = rbf_volume(rbf->rbfa);
515	+	return(rbf);
516	+	}
517	+
518		/* Advect and allocate new RBF along edge (internal call) */
519		RBFNODE *
520		e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
#	Line 403 | Line 525 | e_advect_rbf(const MIGRATION *mig, const FVECT invec,
525		double          t, full_dist;
526		/* get relative position */
527		t = Acos(DOT(invec, mig->rbfv[0]->invec));
528	<	if (t < M_PI/grid_res) {                /* near first DSF */
528	>	if (t <= .001) {                        /* near first DSF */
529		n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
530		rbf = (RBFNODE *)malloc(n);
531		if (rbf == NULL)
#	Line 413 | Line 535 | e_advect_rbf(const MIGRATION *mig, const FVECT invec,
535		return(rbf);
536		}
537		full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
538	<	if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
538	>	if (t >= full_dist-.001) {              /* near second DSF */
539		n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
540		rbf = (RBFNODE *)malloc(n);
541		if (rbf == NULL)
#	Line 449 | Line 571 | tryagain:
571		const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
572		const float         peak0 = rbf0i->peak;
573		const double        rad0 = R2ANG(rbf0i->crad);
574	+	C_COLOR             cc0;
575		FVECT               v0;
576		float               mv;
577		ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
578	+	c_decodeChroma(&cc0, rbf0i->chroma);
579		for (j = 0; j < mtx_ncols(mig); j++)
580		if ((mv = mtx_coef(mig,i,j)) > cthresh) {
581		const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
#	Line 462 | Line 586 | tryagain:
586		rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
587		rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
588		rad0*rad0/rad2;
589	+	if (rbf_colorimetry == RBCtristimulus) {
590	+	C_COLOR cres;
591	+	c_decodeChroma(&cres, rbf1j->chroma);
592	+	c_cmix(&cres, 1.-t, &cc0, t, &cres);
593	+	rbf->rbfa[n].chroma = c_encodeChroma(&cres);
594	+	} else
595	+	rbf->rbfa[n].chroma = c_dfchroma;
596		rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
597		ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
598		geodesic(v, v0, v, t, GEOD_REL);
#	Line 498 | Line 629 | clear_bsdf_rep(void)
629		inp_coverage = 0;
630		single_plane_incident = -1;
631		input_orient = output_orient = 0;
632	+	rbf_colorimetry = RBCunknown;
633		grid_res = GRIDRES;
634	+	memset(bsdf_hist, 0, sizeof(bsdf_hist));
635	+	bsdf_min = 0;
636	+	bsdf_spec_val = 0;
637	+	bsdf_spec_rad = 0;
638		}
639
640		/* Write our BSDF mesh interpolant out to the given binary stream */
#	Line 515 | Line 651 | save_bsdf_rep(FILE *ofp)
651		fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf);
652		fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
653		fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
654	+	fprintf(ofp, "COLORIMETRY=%s\n", RBCident[rbf_colorimetry]);
655		fprintf(ofp, "GRIDRES=%d\n", grid_res);
656		fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
657	+	if ((bsdf_spec_val > bsdf_min) & (bsdf_spec_rad > 0))
658	+	fprintf(ofp, "BSDFSPEC= %f %f\n", bsdf_spec_val, bsdf_spec_rad);
659		fputformat(BSDFREP_FMT, ofp);
660		fputc('\n', ofp);
661	+	putint(BSDFREP_MAGIC, 2, ofp);
662		/* write each DSF */
663		for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
664		putint(rbf->ord, 4, ofp);
#	Line 529 | Line 669 | save_bsdf_rep(FILE *ofp)
669		putint(rbf->nrbf, 4, ofp);
670		for (i = 0; i < rbf->nrbf; i++) {
671		putflt(rbf->rbfa[i].peak, ofp);
672	+	putint(rbf->rbfa[i].chroma, 2, ofp);
673		putint(rbf->rbfa[i].crad, 2, ofp);
674	<	putint(rbf->rbfa[i].gx, 1, ofp);
675	<	putint(rbf->rbfa[i].gy, 1, ofp);
674	>	putint(rbf->rbfa[i].gx, 2, ofp);
675	>	putint(rbf->rbfa[i].gy, 2, ofp);
676		}
677		}
678		putint(-1, 4, ofp);             /* terminator */
#	Line 567 | Line 708 | save_bsdf_rep(FILE *ofp)
708		static int
709		headline(char *s, void *p)
710		{
711	<	char    fmt[32];
711	>	char    fmt[MAXFMTLEN];
712	>	int     i;
713
714	+	if (isheadid(s))
715	+	return(0);
716		if (!strncmp(s, "NAME=", 5)) {
717		strcpy(bsdf_name, s+5);
718		bsdf_name[strlen(bsdf_name)-1] = '\0';
719	+	return(1);
720		}
721		if (!strncmp(s, "MANUFACT=", 9)) {
722		strcpy(bsdf_manuf, s+9);
723		bsdf_manuf[strlen(bsdf_manuf)-1] = '\0';
724	+	return(1);
725		}
726		if (!strncmp(s, "SYMMETRY=", 9)) {
727		inp_coverage = atoi(s+9);
728		single_plane_incident = !inp_coverage;
729	<	return(0);
729	>	return(1);
730		}
731		if (!strncmp(s, "IO_SIDES=", 9)) {
732		sscanf(s+9, "%d %d", &input_orient, &output_orient);
733	<	return(0);
733	>	return(1);
734		}
735	+	if (!strncmp(s, "COLORIMETRY=", 12)) {
736	+	fmt[0] = '\0';
737	+	sscanf(s+12, "%s", fmt);
738	+	for (i = RBCunknown; i >= 0; i--)
739	+	if (!strcmp(fmt, RBCident[i]))
740	+	break;
741	+	if (i < 0)
742	+	return(-1);
743	+	rbf_colorimetry = i;
744	+	return(1);
745	+	}
746		if (!strncmp(s, "GRIDRES=", 8)) {
747		sscanf(s+8, "%d", &grid_res);
748	<	return(0);
748	>	return(1);
749		}
750		if (!strncmp(s, "BSDFMIN=", 8)) {
751		sscanf(s+8, "%lf", &bsdf_min);
752	<	return(0);
752	>	return(1);
753		}
754	<	if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
755	<	return(-1);
754	>	if (!strncmp(s, "BSDFSPEC=", 9)) {
755	>	sscanf(s+9, "%lf %lf", &bsdf_spec_val, &bsdf_spec_rad);
756	>	return(1);
757	>	}
758	>	if (formatval(fmt, s))
759	>	return (strcmp(fmt, BSDFREP_FMT) ? -1 : 0);
760	>	if (sir_headshare != NULL)
761	>	return ((*sir_headshare)(s));
762		return(0);
763		}
764
#	Line 611 | Line 774 | load_bsdf_rep(FILE *ifp)
774		if (ifp == NULL)
775		return(0);
776		if (getheader(ifp, headline, NULL) < 0 || (single_plane_incident < 0) |
777	<	!input_orient | !output_orient) {
777	>	!input_orient | !output_orient |
778	>	(grid_res < 16) | (grid_res > 0xffff)) {
779		fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
780		progname);
781		return(0);
782		}
783	+	if (getint(2, ifp) != BSDFREP_MAGIC) {
784	+	fprintf(stderr, "%s: bad magic number for BSDF interpolant\n",
785	+	progname);
786	+	return(0);
787	+	}
788		memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */
789		while ((rbfh.ord = getint(4, ifp)) >= 0) {
790		RBFNODE         *newrbf;
#	Line 636 | Line 805 | load_bsdf_rep(FILE *ifp)
805		*newrbf = rbfh;
806		for (i = 0; i < rbfh.nrbf; i++) {
807		newrbf->rbfa[i].peak = getflt(ifp);
808	+	newrbf->rbfa[i].chroma = getint(2, ifp) & 0xffff;
809		newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
810	<	newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
811	<	newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
810	>	newrbf->rbfa[i].gx = getint(2, ifp) & 0xffff;
811	>	newrbf->rbfa[i].gy = getint(2, ifp) & 0xffff;
812		}
813		if (feof(ifp))
814		goto badEOF;

Diff Legend

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