[ViewVC] Diff of: radiance/ray/src/cv/pabopto2xml.c

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.10 by greg, Tue Sep 18 02:50:13 2012 UTC vs.
Revision 2.15 by greg, Sun Sep 23 16:45:20 2012 UTC

#	Line 88 \| Line 88 \| *static MIGRATION mig_grid[GRIDRES][GRIDRES];**
88
89		#define round(v) (int)((v) + .5 - ((v) < -.5))
90
91	+	char *progname;
92	+
93	+	#ifdef DEBUG /* percentage to cull (<0 to turn off) */
94	+	int pctcull = -1;
95	+	#else
96	+	int pctcull = 90;
97	+	#endif
98	+	/* sampling order (set by data density) */
99	+	int samp_order = 0;
100	+
101		/* Compute volume associated with Gaussian lobe */
102		static double
103		rbf_volume(const RBFVAL *rbfp)
#	Line 136 \| Line 146 \| *eval_rbfrep(const RBFNODE rp, const FVECT outvec)**
146		double sig2;
147		int n;
148
149	+	if (rp == NULL)
150	+	return(.0);
151		rbfp = rp->rbfa;
152		for (n = rp->nrbf; n--; rbfp++) {
153		ovec_from_pos(odir, rbfp->gx, rbfp->gy);
#	Line 152 \| Line 164 \| static void
164		insert_dsf(RBFNODE *newrbf)
165		{
166		RBFNODE rbf, rbf_last;
167	<
167	>	/* check for redundant meas. */
168	>	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
169	>	if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
170	>	fputs("Duplicate incident measurement (ignored)\n", stderr);
171	>	free(newrbf);
172	>	return;
173	>	}
174		/* keep in ascending theta order */
175		for (rbf_last = NULL, rbf = dsf_list;
176		single_plane_incident & (rbf != NULL);
#	Line 173 \| Line 191 \| static RBFNODE *
191		make_rbfrep(void)
192		{
193		int niter = 16;
194	+	int minrad = ANG2R(pow(2., 1.-samp_order));
195		double lastVar, thisVar = 100.;
196		int nn;
197		RBFNODE *newnode;
#	Line 204 \| Line 223 \| make_rbfrep(void)
223		newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
224		newnode->rbfa[nn].gx = i;
225		newnode->rbfa[nn].gy = j;
226	+	if (newnode->rbfa[nn].crad < minrad)
227	+	minrad = newnode->rbfa[nn].crad;
228		++nn;
229		}
230		/* iterate to improve interpolation accuracy */
231		do {
232	<	double dsum = .0, dsum2 = .0;
232	>	double dsum = 0, dsum2 = 0;
233		nn = 0;
234		for (i = 0; i < GRIDRES; i++)
235		for (j = 0; j < GRIDRES; j++)
#	Line 236 \| Line 257 \| make_rbfrep(void)
257		newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
258
259		insert_dsf(newnode);
260	+	/* adjust sampling resolution */
261	+	samp_order = log(2./R2ANG(minrad))/M_LN2 + .5;
262	+
263		return(newnode);
264		}
265
#	Line 544 \| Line 568 \| *migration_step(MIGRATION mig, double src_rem, double*
568		{
569		static double *src_cost = NULL;
570		int n_alloc = 0;
571	<	const double maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
571	>	const double maxamt = .1; /* 2./(mtx_nrows(mig)mtx_ncols(mig)); /
572		double amt = 0;
573		struct {
574		int s, d; /* source and destination */
#	Line 606 \| Line 630 \| *migration_step(MIGRATION mig, double src_rem, double*
630		return(best.amt);
631		}
632
633	+	#ifdef DEBUG
634	+	static char *
635	+	thetaphi(const FVECT v)
636	+	{
637	+	static char buf[128];
638	+	double theta, phi;
639	+
640	+	theta = 180./M_PI*acos(v[2]);
641	+	phi = 180./M_PI*atan2(v[1],v[0]);
642	+	sprintf(buf, "(%.0f,%.0f)", theta, phi);
643	+
644	+	return(buf);
645	+	}
646	+	#endif
647	+
648		/* Compute (and insert) migration along directed edge */
649		static MIGRATION *
650	<	make_migration(RBFNODE from_rbf, RBFNODE to_rbf)
650	>	make_migration(RBFNODE from_rbf, RBFNODE to_rbf, int creat_only)
651		{
652		const double end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
653	<	float *pmtx = price_routes(from_rbf, to_rbf);
654	<	MIGRATION newmig = (MIGRATION )malloc(sizeof(MIGRATION) +
655	<	sizeof(float) *
617	<	(from_rbf->nrbf*to_rbf->nrbf - 1));
618	<	double src_rem = (double )malloc(sizeof(double)*from_rbf->nrbf);
619	<	double dst_rem = (double )malloc(sizeof(double)*to_rbf->nrbf);
653	>	float *pmtx;
654	>	MIGRATION *newmig;
655	>	double src_rem, dst_rem;
656		double total_rem = 1.;
657		int i;
658	<
658	>	/* check if exists already */
659	>	for (newmig = from_rbf->ejl; newmig != NULL;
660	>	newmig = nextedge(from_rbf,newmig))
661	>	if (newmig->rbfv[1] == to_rbf)
662	>	return(creat_only ? (MIGRATION *)NULL : newmig);
663	>	/* else allocate */
664	>	pmtx = price_routes(from_rbf, to_rbf);
665	>	newmig = (MIGRATION )malloc(sizeof(MIGRATION) + sizeof(float)
666	>	(from_rbf->nrbf*to_rbf->nrbf - 1));
667	>	src_rem = (double )malloc(sizeof(double)from_rbf->nrbf);
668	>	dst_rem = (double )malloc(sizeof(double)to_rbf->nrbf);
669		if ((newmig == NULL) \| (src_rem == NULL) \| (dst_rem == NULL)) {
670		fputs("Out of memory in make_migration()\n", stderr);
671		exit(1);
672		}
673		#ifdef DEBUG
674		{
675	<	double theta, phi;
676	<	theta = acos(from_rbf->invec[2])*(180./M_PI);
677	<	phi = atan2(from_rbf->invec[1],from_rbf->invec[0])*(180./M_PI);
632	<	fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
633	<	round(theta), round(phi));
634	<	theta = acos(to_rbf->invec[2])*(180./M_PI);
635	<	phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
636	<	fprintf(stderr, "(%d,%d)\n", round(theta), round(phi));
675	>	fprintf(stderr, "Building path from (theta,phi) %s ",
676	>	thetaphi(from_rbf->invec));
677	>	fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
678		}
679		#endif
680		newmig->next = NULL;
#	Line 647 \| Line 688 \| *make_migration(RBFNODE from_rbf, RBFNODE to_rbf)*
688		for (i = to_rbf->nrbf; i--; )
689		dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
690		/* move a bit at a time */
691	<	while (total_rem > end_thresh)
691	>	while (total_rem > end_thresh) {
692		total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
693	+	#ifdef DEBUG
694	+	/* fputc('.', stderr); */
695	+	fprintf(stderr, "\n%.9f remaining...", total_rem);
696	+	#endif
697	+	}
698	+	#ifdef DEBUG
699	+	fputs("done.\n", stderr);
700	+	#endif
701
702		free(pmtx); /* free working arrays */
703		free(src_rem);
#	Line 700 \| Line 749 \| *get_triangles(RBFNODE rbfv[2], const MIGRATION mig)*
749		RBFNODE *tv;
750
751		rbfv[0] = rbfv[1] = NULL;
752	+	if (mig == NULL)
753	+	return(0);
754		for (ej = mig->rbfv[0]->ejl; ej != NULL;
755		ej = nextedge(mig->rbfv[0],ej)) {
756		if (ej == mig)
#	Line 716 \| Line 767 \| *get_triangles(RBFNODE rbfv[2], const MIGRATION mig)*
767		return((rbfv[0] != NULL) + (rbfv[1] != NULL));
768		}
769
770	+	/* Check if prospective vertex would create overlapping triangle */
771	+	static int
772	+	overlaps_tri(const RBFNODE bv0, const RBFNODE bv1, const RBFNODE *pv)
773	+	{
774	+	const MIGRATION *ej;
775	+	RBFNODE *vother[2];
776	+	int im_rev;
777	+	/* find shared edge in mesh */
778	+	for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
779	+	const RBFNODE *tv = opp_rbf(pv,ej);
780	+	if (tv == bv0) {
781	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
782	+	ej->rbfv[1]->invec, bv1->invec);
783	+	break;
784	+	}
785	+	if (tv == bv1) {
786	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
787	+	ej->rbfv[1]->invec, bv0->invec);
788	+	break;
789	+	}
790	+	}
791	+	if (!get_triangles(vother, ej)) /* triangle on same side? */
792	+	return(0);
793	+	return(vother[im_rev] != NULL);
794	+	}
795	+
796		/* Find context hull vertex to complete triangle (oriented call) */
797		static RBFNODE *
798		find_chull_vert(const RBFNODE rbf0, const RBFNODE rbf1)
799		{
800	<	FVECT vmid, vor;
800	>	FVECT vmid, vejn, vp;
801		RBFNODE rbf, rbfbest = NULL;
802	<	double dprod2, bestdprod2 = 0.5;
802	>	double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
803
804	+	VSUB(vejn, rbf1->invec, rbf0->invec);
805		VADD(vmid, rbf0->invec, rbf1->invec);
806	<	if (normalize(vmid) == 0)
806	>	if (normalize(vejn) == 0 \|\| normalize(vmid) == 0)
807		return(NULL);
808		/* XXX exhaustive search */
809	+	/* Find triangle with minimum rotation from perpendicular */
810		for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
811		if ((rbf == rbf0) \| (rbf == rbf1))
812		continue;
813	<	tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
814	<	dprod2 = DOT(vor, vmid);
736	<	if (dprod2 <= FTINY)
813	>	tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
814	>	if (DOT(vp, vmid) <= FTINY)
815		continue; /* wrong orientation */
816	<	dprod2 *= dprod2 / DOT(vor,vor);
817	<	if (dprod2 > bestdprod2) { /* more convex than prev? */
818	<	rbfbest = rbf;
819	<	bestdprod2 = dprod2;
820	<	}
816	>	area2 = .25*DOT(vp,vp);
817	>	VSUB(vp, rbf->invec, rbf0->invec);
818	>	dprod = -DOT(vp, vejn);
819	>	VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
820	>	dprod = DOT(vp, vmid) / VLEN(vp);
821	>	if (dprod <= bestdprod + FTINY(1 - 2(area2 < bestarea2)))
822	>	continue; /* found better already */
823	>	if (overlaps_tri(rbf0, rbf1, rbf))
824	>	continue; /* overlaps another triangle */
825	>	rbfbest = rbf;
826	>	bestdprod = dprod; /* new one to beat */
827	>	bestarea2 = area2;
828		}
829	<	return(rbf);
829	>	return(rbfbest);
830		}
831
832		/* Create new migration edge and grow mesh recursively around it */
833		static void
834	<	mesh_from_edge(RBFNODE rbf0, RBFNODE rbf1)
834	>	mesh_from_edge(MIGRATION *edge)
835		{
836	<	MIGRATION *newej;
836	>	MIGRATION ej0, ej1;
837		RBFNODE *tvert[2];
838	<
839	<	if (rbf0 < rbf1) /* avoid migration loops */
840	<	newej = make_migration(rbf0, rbf1);
756	<	else
757	<	newej = make_migration(rbf1, rbf0);
838	>
839	>	if (edge == NULL)
840	>	return;
841		/* triangle on either side? */
842	<	get_triangles(tvert, newej);
843	<	if (tvert[0] == NULL) { /* recurse on new right edge */
844	<	tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
842	>	get_triangles(tvert, edge);
843	>	if (tvert[0] == NULL) { /* grow mesh on right */
844	>	tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
845		if (tvert[0] != NULL) {
846	<	mesh_from_edge(rbf0, tvert[0]);
847	<	mesh_from_edge(rbf1, tvert[0]);
846	>	if (tvert[0] > edge->rbfv[0])
847	>	ej0 = make_migration(edge->rbfv[0], tvert[0], 1);
848	>	else
849	>	ej0 = make_migration(tvert[0], edge->rbfv[0], 1);
850	>	if (tvert[0] > edge->rbfv[1])
851	>	ej1 = make_migration(edge->rbfv[1], tvert[0], 1);
852	>	else
853	>	ej1 = make_migration(tvert[0], edge->rbfv[1], 1);
854	>	mesh_from_edge(ej0);
855	>	mesh_from_edge(ej1);
856		}
857	<	}
858	<	if (tvert[1] == NULL) { /* recurse on new left edge */
768	<	tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
857	>	} else if (tvert[1] == NULL) { /* grow mesh on left */
858	>	tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
859		if (tvert[1] != NULL) {
860	<	mesh_from_edge(rbf0, tvert[1]);
861	<	mesh_from_edge(rbf1, tvert[1]);
860	>	if (tvert[1] > edge->rbfv[0])
861	>	ej0 = make_migration(edge->rbfv[0], tvert[1], 1);
862	>	else
863	>	ej0 = make_migration(tvert[1], edge->rbfv[0], 1);
864	>	if (tvert[1] > edge->rbfv[1])
865	>	ej1 = make_migration(edge->rbfv[1], tvert[1], 1);
866	>	else
867	>	ej1 = make_migration(tvert[1], edge->rbfv[1], 1);
868	>	mesh_from_edge(ej0);
869	>	mesh_from_edge(ej1);
870		}
871		}
872		}
873
874	+	#ifdef DEBUG
875	+	#include "random.h"
876	+	#include "bmpfile.h"
877	+	/* Hash pointer to byte value (must return 0 for NULL) */
878	+	static int
879	+	byte_hash(const void *p)
880	+	{
881	+	size_t h = (size_t)p;
882	+	h ^= (size_t)p >> 8;
883	+	h ^= (size_t)p >> 16;
884	+	h ^= (size_t)p >> 24;
885	+	return(h & 0xff);
886	+	}
887	+	/* Write out BMP image showing edges */
888	+	static void
889	+	write_edge_image(const char *fname)
890	+	{
891	+	BMPHeader *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
892	+	BMPWriter *wtr;
893	+	int i, j;
894	+
895	+	fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
896	+	hdr->compr = BI_RLE8;
897	+	for (i = 256; --i; ) { /* assign random color map */
898	+	hdr->palette[i].r = random() & 0xff;
899	+	hdr->palette[i].g = random() & 0xff;
900	+	hdr->palette[i].b = random() & 0xff;
901	+	/* reject dark colors */
902	+	i += (hdr->palette[i].r + hdr->palette[i].g +
903	+	hdr->palette[i].b < 128);
904	+	}
905	+	hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
906	+	/* open output */
907	+	wtr = BMPopenOutputFile(fname, hdr);
908	+	if (wtr == NULL) {
909	+	free(hdr);
910	+	return;
911	+	}
912	+	for (i = 0; i < GRIDRES; i++) { /* write scanlines */
913	+	for (j = 0; j < GRIDRES; j++)
914	+	wtr->scanline[j] = byte_hash(mig_grid[i][j]);
915	+	if (BMPwriteScanline(wtr) != BIR_OK)
916	+	break;
917	+	}
918	+	BMPcloseOutput(wtr); /* close & clean up */
919	+	}
920	+	#endif
921	+
922		/* Draw edge list into mig_grid array */
923		static void
924		draw_edges()
#	Line 816 \| Line 962 \| draw_edges()
962		if (nnull)
963		fprintf(stderr, "Warning: %d of %d edges are null\n",
964		nnull, ntot);
965	+	#ifdef DEBUG
966	+	write_edge_image("bsdf_edges.bmp");
967	+	#endif
968		}
969
970		/* Build our triangle mesh from recorded RBFs */
#	Line 823 \| Line 972 \| static void
972		build_mesh()
973		{
974		double best2 = M_PI*M_PI;
975	<	RBFNODE rbf, rbf_near = NULL;
975	>	RBFNODE *shrt_edj[2];
976	>	RBFNODE rbf0, rbf1;
977		/* check if isotropic */
978		if (single_plane_incident) {
979	<	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
980	<	if (rbf->next != NULL)
981	<	make_migration(rbf, rbf->next);
979	>	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
980	>	if (rbf0->next != NULL)
981	>	make_migration(rbf0, rbf0->next, 1);
982		return;
983		}
984	<	/* find RBF nearest to head */
985	<	if (dsf_list == NULL)
986	<	return;
987	<	for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
988	<	double dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
984	>	/* start w/ shortest edge */
985	>	shrt_edj[0] = shrt_edj[1] = NULL;
986	>	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
987	>	for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
988	>	double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
989		if (dist2 < best2) {
990	<	rbf_near = rbf;
990	>	shrt_edj[0] = rbf0;
991	>	shrt_edj[1] = rbf1;
992		best2 = dist2;
993		}
994		}
995	<	if (rbf_near == NULL) {
996	<	fputs("Cannot find nearest point for first edge\n", stderr);
995	>	if (shrt_edj[0] == NULL) {
996	>	fputs("Cannot find shortest edge\n", stderr);
997		exit(1);
998		}
999		/* build mesh from this edge */
1000	<	mesh_from_edge(dsf_list, rbf_near);
1000	>	if (shrt_edj[0] < shrt_edj[1])
1001	>	mesh_from_edge(make_migration(shrt_edj[0], shrt_edj[1], 0));
1002	>	else
1003	>	mesh_from_edge(make_migration(shrt_edj[1], shrt_edj[0], 0));
1004		/* draw edge list into grid */
1005		draw_edges();
1006		}
#	Line 870 \| Line 1024 \| *identify_tri(MIGRATION miga[3], unsigned char vmap[GR**
1024		return(1);
1025		if (miga[i] != NULL)
1026		continue;
873	–	while (i > 0 && miga[i-1] > mig_grid[px][py]) {
874	–	miga[i] = miga[i-1];
875	–	--i; /* order vertices by pointer */
876	–	}
1027		miga[i] = mig_grid[px][py];
1028		return(1);
1029		}
#	Line 891 \| Line 1041 \| *identify_tri(MIGRATION miga[3], unsigned char vmap[GR**
1041		return(1); /* this neighborhood done */
1042		}
1043
1044	+	/* Insert vertex in ordered list */
1045	+	static void
1046	+	insert_vert(RBFNODE *vlist, RBFNODE v)
1047	+	{
1048	+	int i, j;
1049	+
1050	+	for (i = 0; vlist[i] != NULL; i++) {
1051	+	if (v == vlist[i])
1052	+	return;
1053	+	if (v < vlist[i])
1054	+	break;
1055	+	}
1056	+	for (j = i; vlist[j] != NULL; j++)
1057	+	;
1058	+	while (j > i) {
1059	+	vlist[j] = vlist[j-1];
1060	+	--j;
1061	+	}
1062	+	vlist[i] = v;
1063	+	}
1064	+
1065	+	/* Sort triangle edges in standard order */
1066	+	static int
1067	+	order_triangle(MIGRATION *miga[3])
1068	+	{
1069	+	RBFNODE *vert[7];
1070	+	MIGRATION *ord[3];
1071	+	int i;
1072	+	/* order vertices, first */
1073	+	memset(vert, 0, sizeof(vert));
1074	+	for (i = 3; i--; ) {
1075	+	if (miga[i] == NULL)
1076	+	return(0);
1077	+	insert_vert(vert, miga[i]->rbfv[0]);
1078	+	insert_vert(vert, miga[i]->rbfv[1]);
1079	+	}
1080	+	/* should be just 3 vertices */
1081	+	if ((vert[3] == NULL) \| (vert[4] != NULL))
1082	+	return(0);
1083	+	/* identify edge 0 */
1084	+	for (i = 3; i--; )
1085	+	if (miga[i]->rbfv[0] == vert[0] &&
1086	+	miga[i]->rbfv[1] == vert[1]) {
1087	+	ord[0] = miga[i];
1088	+	break;
1089	+	}
1090	+	if (i < 0)
1091	+	return(0);
1092	+	/* identify edge 1 */
1093	+	for (i = 3; i--; )
1094	+	if (miga[i]->rbfv[0] == vert[1] &&
1095	+	miga[i]->rbfv[1] == vert[2]) {
1096	+	ord[1] = miga[i];
1097	+	break;
1098	+	}
1099	+	if (i < 0)
1100	+	return(0);
1101	+	/* identify edge 2 */
1102	+	for (i = 3; i--; )
1103	+	if (miga[i]->rbfv[0] == vert[0] &&
1104	+	miga[i]->rbfv[1] == vert[2]) {
1105	+	ord[2] = miga[i];
1106	+	break;
1107	+	}
1108	+	if (i < 0)
1109	+	return(0);
1110	+	/* reassign order */
1111	+	miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1112	+	return(1);
1113	+	}
1114	+
1115		/* Find edge(s) for interpolating the given incident vector */
1116		static int
1117		get_interp(MIGRATION *miga[3], const FVECT invec)
#	Line 913 \| Line 1134 \| *get_interp(MIGRATION miga[3], const FVECT invec)**
1134		}
1135		return(0); /* outside range! */
1136		}
1137	<	{ /* else use triagnle mesh */
1137	>	{ /* else use triangle mesh */
1138		unsigned char floodmap[GRIDRES][(GRIDRES+7)/8];
1139		int pstart[2];
1140	+	RBFNODE *vother;
1141	+	MIGRATION *ej;
1142	+	int i;
1143
1144		pos_from_vec(pstart, invec);
1145		memset(floodmap, 0, sizeof(floodmap));
#	Line 926 \| Line 1150 \| *get_interp(MIGRATION miga[3], const FVECT invec)**
1150		return(0); /* should never happen */
1151		if (miga[1] == NULL)
1152		return(1); /* on edge */
1153	<	return(3); /* else in triangle */
1153	>	/* verify triangle */
1154	>	if (!order_triangle(miga)) {
1155	>	#ifdef DEBUG
1156	>	fputs("Munged triangle in get_interp()\n", stderr);
1157	>	#endif
1158	>	vother = NULL; /* find triangle from edge */
1159	>	for (i = 3; i--; ) {
1160	>	RBFNODE *tpair[2];
1161	>	if (get_triangles(tpair, miga[i]) &&
1162	>	(vother = tpair[ is_rev_tri(
1163	>	miga[i]->rbfv[0]->invec,
1164	>	miga[i]->rbfv[1]->invec,
1165	>	invec) ]) != NULL)
1166	>	break;
1167	>	}
1168	>	if (vother == NULL) { /* couldn't find 3rd vertex */
1169	>	#ifdef DEBUG
1170	>	fputs("No triangle in get_interp()\n", stderr);
1171	>	#endif
1172	>	return(0);
1173	>	}
1174	>	/* reassign other two edges */
1175	>	for (ej = vother->ejl; ej != NULL;
1176	>	ej = nextedge(vother,ej)) {
1177	>	RBFNODE *vorig = opp_rbf(vother,ej);
1178	>	if (vorig == miga[i]->rbfv[0])
1179	>	miga[(i+1)%3] = ej;
1180	>	else if (vorig == miga[i]->rbfv[1])
1181	>	miga[(i+2)%3] = ej;
1182	>	}
1183	>	if (!order_triangle(miga)) {
1184	>	#ifdef DEBUG
1185	>	fputs("Bad triangle in get_interp()\n", stderr);
1186	>	#endif
1187	>	return(0);
1188	>	}
1189	>	}
1190	>	return(3); /* return in standard order */
1191		}
1192		}
1193
#	Line 961 \| Line 1222 \| *e_advect_rbf(const MIGRATION mig, const FVECT invec)**
1222		for (i = 0; i < mtx_nrows(mig); i++)
1223		for (j = 0; j < mtx_ncols(mig); j++)
1224		n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
1225	<
1225	>	#ifdef DEBUG
1226	>	fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
1227	>	mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
1228	>	#endif
1229		rbf = (RBFNODE )malloc(sizeof(RBFNODE) + sizeof(RBFVAL)(n-1));
1230		if (rbf == NULL)
1231		goto memerr;
#	Line 1008 \| Line 1272 \| advect_rbf(const FVECT invec)
1272		{
1273		MIGRATION *miga[3];
1274		RBFNODE *rbf;
1275	<	int n, i, j;
1275	>	float mbfact, mcfact;
1276	>	int n, i, j, k;
1277	>	FVECT v0, v1, v2;
1278		double s, t;
1279
1280		if (!get_interp(miga, invec)) /* can't interpolate? */
1281		return(NULL);
1282	<	if (miga[1] == NULL) /* along edge? */
1282	>	if (miga[1] == NULL) /* advect along edge? */
1283		return(e_advect_rbf(miga[0], invec));
1284	+	#ifdef DEBUG
1285	+	if (miga[0]->rbfv[0] != miga[2]->rbfv[0] \|
1286	+	miga[0]->rbfv[1] != miga[1]->rbfv[0] \|
1287	+	miga[1]->rbfv[1] != miga[2]->rbfv[1]) {
1288	+	fputs("Triangle vertex screw-up!\n", stderr);
1289	+	exit(1);
1290	+	}
1291	+	#endif
1292		/* figure out position */
1293	<
1293	>	fcross(v0, miga[2]->rbfv[0]->invec, miga[2]->rbfv[1]->invec);
1294	>	normalize(v0);
1295	>	fcross(v2, miga[1]->rbfv[0]->invec, miga[1]->rbfv[1]->invec);
1296	>	normalize(v2);
1297	>	fcross(v1, invec, miga[1]->rbfv[1]->invec);
1298	>	normalize(v1);
1299	>	s = acos(DOT(v0,v1)) / acos(DOT(v0,v2));
1300	>	geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
1301	>	s, GEOD_REL);
1302	>	t = acos(DOT(v1,invec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
1303	>	n = 0; /* count migrating particles */
1304	>	for (i = 0; i < mtx_nrows(miga[0]); i++)
1305	>	for (j = 0; j < mtx_ncols(miga[0]); j++)
1306	>	for (k = (miga[0]->mtx[mtx_ndx(miga[0],i,j)] > FTINY) *
1307	>	mtx_ncols(miga[2]); k--; )
1308	>	n += (miga[2]->mtx[mtx_ndx(miga[2],i,k)] > FTINY &&
1309	>	miga[1]->mtx[mtx_ndx(miga[1],j,k)] > FTINY);
1310	>	#ifdef DEBUG
1311	>	fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
1312	>	miga[0]->rbfv[0]->nrbf, miga[0]->rbfv[1]->nrbf,
1313	>	miga[2]->rbfv[1]->nrbf, n);
1314	>	#endif
1315		rbf = (RBFNODE )malloc(sizeof(RBFNODE) + sizeof(RBFVAL)(n-1));
1316		if (rbf == NULL) {
1317		fputs("Out of memory in advect_rbf()\n", stderr);
#	Line 1024 \| Line 1319 \| advect_rbf(const FVECT invec)
1319		}
1320		rbf->next = NULL; rbf->ejl = NULL;
1321		VCOPY(rbf->invec, invec);
1027	–	rbf->vtotal = 0;
1322		rbf->nrbf = n;
1323	<
1323	>	n = 0; /* compute RBF lobes */
1324	>	mbfact = s * miga[0]->rbfv[1]->vtotal/miga[0]->rbfv[0]->vtotal *
1325	>	(1.-t + t*miga[1]->rbfv[1]->vtotal/miga[1]->rbfv[0]->vtotal);
1326	>	mcfact = (1.-s) *
1327	>	(1.-t + t*miga[2]->rbfv[1]->vtotal/miga[2]->rbfv[0]->vtotal);
1328	>	for (i = 0; i < mtx_nrows(miga[0]); i++) {
1329	>	const RBFVAL *rbf0i = &miga[0]->rbfv[0]->rbfa[i];
1330	>	const float w0i = rbf0i->peak;
1331	>	const double rad0i = R2ANG(rbf0i->crad);
1332	>	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1333	>	for (j = 0; j < mtx_ncols(miga[0]); j++) {
1334	>	const float ma = miga[0]->mtx[mtx_ndx(miga[0],i,j)];
1335	>	const RBFVAL *rbf1j;
1336	>	double rad1j, srad2;
1337	>	if (ma <= FTINY)
1338	>	continue;
1339	>	rbf1j = &miga[0]->rbfv[1]->rbfa[j];
1340	>	rad1j = R2ANG(rbf1j->crad);
1341	>	srad2 = (1.-s)(1.-t)rad0irad0i + s(1.-t)rad1jrad1j;
1342	>	ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
1343	>	geodesic(v1, v0, v1, s, GEOD_REL);
1344	>	for (k = 0; k < mtx_ncols(miga[2]); k++) {
1345	>	float mb = miga[1]->mtx[mtx_ndx(miga[1],j,k)];
1346	>	float mc = miga[2]->mtx[mtx_ndx(miga[2],i,k)];
1347	>	const RBFVAL *rbf2k;
1348	>	double rad2k;
1349	>	FVECT vout;
1350	>	int pos[2];
1351	>	if ((mb <= FTINY) \| (mc <= FTINY))
1352	>	continue;
1353	>	rbf2k = &miga[2]->rbfv[1]->rbfa[k];
1354	>	rbf->rbfa[n].peak = w0i * ma * (mbmbfact + mcmcfact);
1355	>	rad2k = R2ANG(rbf2k->crad);
1356	>	rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + trad2krad2k));
1357	>	ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
1358	>	geodesic(vout, v1, v2, t, GEOD_REL);
1359	>	pos_from_vec(pos, vout);
1360	>	rbf->rbfa[n].gx = pos[0];
1361	>	rbf->rbfa[n].gy = pos[1];
1362	>	++n;
1363	>	}
1364	>	}
1365	>	}
1366	>	rbf->vtotal = miga[0]->rbfv[0]->vtotal * (mbfact + mcfact);
1367		return(rbf);
1368		}
1369
1370	+	/* Interpolate and output isotropic BSDF data */
1371	+	static void
1372	+	interp_isotropic()
1373	+	{
1374	+	const int sqres = 1<<samp_order;
1375	+	FILE *ofp = NULL;
1376	+	char cmd[128];
1377	+	int ix, ox, oy;
1378	+	FVECT ivec, ovec;
1379	+	double bsdf;
1380	+	#if DEBUG
1381	+	fprintf(stderr, "Writing isotropic order %d ", samp_order);
1382	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1383	+	else fputs("raw data\n", stderr);
1384	+	#endif
1385	+	if (pctcull >= 0) { /* begin output */
1386	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1387	+	pctcull, samp_order);
1388	+	fflush(stdout);
1389	+	ofp = popen(cmd, "w");
1390	+	if (ofp == NULL) {
1391	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1392	+	exit(1);
1393	+	}
1394	+	} else
1395	+	fputs("{\n", stdout);
1396	+	/* run through directions */
1397	+	for (ix = 0; ix < sqres/2; ix++) {
1398	+	RBFNODE *rbf;
1399	+	SDsquare2disk(ivec, (ix+.5)/sqres, .5);
1400	+	ivec[2] = input_orient *
1401	+	sqrt(1. - ivec[0]ivec[0] - ivec[1]ivec[1]);
1402	+	rbf = advect_rbf(ivec);
1403	+	for (ox = 0; ox < sqres; ox++)
1404	+	for (oy = 0; oy < sqres; oy++) {
1405	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1406	+	ovec[2] = output_orient *
1407	+	sqrt(1. - ovec[0]ovec[0] - ovec[1]ovec[1]);
1408	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1409	+	if (pctcull >= 0)
1410	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1411	+	else
1412	+	printf("\t%.3e\n", bsdf);
1413	+	}
1414	+	free(rbf);
1415	+	}
1416	+	if (pctcull >= 0) { /* finish output */
1417	+	if (pclose(ofp)) {
1418	+	fprintf(stderr, "Error running '%s'\n", cmd);
1419	+	exit(1);
1420	+	}
1421	+	} else {
1422	+	for (ix = sqressqressqres/2; ix--; )
1423	+	fputs("\t0\n", stdout);
1424	+	fputs("}\n", stdout);
1425	+	}
1426	+	}
1427	+
1428	+	/* Interpolate and output anisotropic BSDF data */
1429	+	static void
1430	+	interp_anisotropic()
1431	+	{
1432	+	const int sqres = 1<<samp_order;
1433	+	FILE *ofp = NULL;
1434	+	char cmd[128];
1435	+	int ix, iy, ox, oy;
1436	+	FVECT ivec, ovec;
1437	+	double bsdf;
1438	+	#if DEBUG
1439	+	fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1440	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1441	+	else fputs("raw data\n", stderr);
1442	+	#endif
1443	+	if (pctcull >= 0) { /* begin output */
1444	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1445	+	pctcull, samp_order);
1446	+	fflush(stdout);
1447	+	ofp = popen(cmd, "w");
1448	+	if (ofp == NULL) {
1449	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1450	+	exit(1);
1451	+	}
1452	+	} else
1453	+	fputs("{\n", stdout);
1454	+	/* run through directions */
1455	+	for (ix = 0; ix < sqres; ix++)
1456	+	for (iy = 0; iy < sqres; iy++) {
1457	+	RBFNODE *rbf;
1458	+	SDsquare2disk(ivec, (ix+.5)/sqres, (iy+.5)/sqres);
1459	+	ivec[2] = input_orient *
1460	+	sqrt(1. - ivec[0]ivec[0] - ivec[1]ivec[1]);
1461	+	rbf = advect_rbf(ivec);
1462	+	for (ox = 0; ox < sqres; ox++)
1463	+	for (oy = 0; oy < sqres; oy++) {
1464	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1465	+	ovec[2] = output_orient *
1466	+	sqrt(1. - ovec[0]ovec[0] - ovec[1]ovec[1]);
1467	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1468	+	if (pctcull >= 0)
1469	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1470	+	else
1471	+	printf("\t%.3e\n", bsdf);
1472	+	}
1473	+	free(rbf);
1474	+	}
1475	+	if (pctcull >= 0) { /* finish output */
1476	+	if (pclose(ofp)) {
1477	+	fprintf(stderr, "Error running '%s'\n", cmd);
1478	+	exit(1);
1479	+	}
1480	+	} else
1481	+	fputs("}\n", stdout);
1482	+	}
1483	+
1484		#if 1
1485	+	/* Read in BSDF files and interpolate as tensor tree representation */
1486	+	int
1487	+	main(int argc, char *argv[])
1488	+	{
1489	+	RBFNODE *rbf;
1490	+	double bsdf;
1491	+	int i;
1492	+
1493	+	progname = argv[0];
1494	+	if (argc > 2 && !strcmp(argv[1], "-t")) {
1495	+	pctcull = atoi(argv[2]);
1496	+	argv += 2; argc -= 2;
1497	+	}
1498	+	if (argc < 3) {
1499	+	fprintf(stderr,
1500	+	"Usage: %s [-t pctcull] meas1.dat meas2.dat .. > bsdf.xml\n",
1501	+	progname);
1502	+	return(1);
1503	+	}
1504	+	for (i = 1; i < argc; i++) { /* compile measurements */
1505	+	if (!load_pabopto_meas(argv[i]))
1506	+	return(1);
1507	+	compute_radii();
1508	+	cull_values();
1509	+	make_rbfrep();
1510	+	}
1511	+	build_mesh(); /* create interpolation */
1512	+	/* xml_prologue(); /* start XML output */
1513	+	if (single_plane_incident) /* resample dist. */
1514	+	interp_isotropic();
1515	+	else
1516	+	interp_anisotropic();
1517	+	/* xml_epilogue(); /* finish XML output */
1518	+	return(0);
1519	+	}
1520	+	#else
1521		/* Test main produces a Radiance model from the given input file */
1522		int
1523		main(int argc, char *argv[])

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing ray/src/cv/pabopto2xml.c (file contents): Revision 2.10 by greg, Tue Sep 18 02:50:13 2012 UTC vs. Revision 2.15 by greg, Sun Sep 23 16:45:20 2012 UTC

Diff Legend

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.10 by greg, Tue Sep 18 02:50:13 2012 UTC vs.
Revision 2.15 by greg, Sun Sep 23 16:45:20 2012 UTC