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

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.12 by greg, Thu Sep 20 01:23:36 2012 UTC vs.
Revision 2.16 by greg, Sun Oct 14 22:31:20 2012 UTC

#	Line 8 \| Line 8 \| static const char RCSid[] = "$Id$";
8		* G.Ward
9		*/
10
11	+	#ifndef _WIN32
12	+	#include <unistd.h>
13	+	#include <sys/wait.h>
14	+	#include <sys/mman.h>
15	+	#endif
16		#define _USE_MATH_DEFINES
17		#include <stdio.h>
18		#include <stdlib.h>
#	Line 89 \| Line 94 \| *static MIGRATION mig_grid[GRIDRES][GRIDRES];**
94		#define round(v) (int)((v) + .5 - ((v) < -.5))
95
96		char *progname;
97	<	/* percentage to cull (<0 to turn off) */
97	>
98	>	#ifdef DEBUG /* percentage to cull (<0 to turn off) */
99	>	int pctcull = -1;
100	>	#else
101		int pctcull = 90;
102	<	/* sampling order */
102	>	#endif
103	>	/* number of processes to run */
104	>	int nprocs = 1;
105	>	/* number of children (-1 in child) */
106	>	int nchild = 0;
107	>
108	>	/* sampling order (set by data density) */
109		int samp_order = 0;
110
111		/* Compute volume associated with Gaussian lobe */
#	Line 160 \| Line 174 \| static void
174		insert_dsf(RBFNODE *newrbf)
175		{
176		RBFNODE rbf, rbf_last;
177	<
177	>	/* check for redundant meas. */
178	>	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
179	>	if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
180	>	fputs("Duplicate incident measurement (ignored)\n", stderr);
181	>	free(newrbf);
182	>	return;
183	>	}
184		/* keep in ascending theta order */
185		for (rbf_last = NULL, rbf = dsf_list;
186		single_plane_incident & (rbf != NULL);
#	Line 219 \| Line 239 \| make_rbfrep(void)
239		}
240		/* iterate to improve interpolation accuracy */
241		do {
242	<	double dsum = .0, dsum2 = .0;
242	>	double dsum = 0, dsum2 = 0;
243		nn = 0;
244		for (i = 0; i < GRIDRES; i++)
245		for (j = 0; j < GRIDRES; j++)
#	Line 248 \| Line 268 \| make_rbfrep(void)
268
269		insert_dsf(newnode);
270		/* adjust sampling resolution */
271	<	samp_order = log(2./R2ANG(minrad))/log(2.) + .5;
271	>	samp_order = log(2./R2ANG(minrad))/M_LN2 + .5;
272
273		return(newnode);
274		}
#	Line 558 \| Line 578 \| *migration_step(MIGRATION mig, double src_rem, double*
578		{
579		static double *src_cost = NULL;
580		int n_alloc = 0;
581	<	const double maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
581	>	const double maxamt = .1; /* 2./(mtx_nrows(mig)mtx_ncols(mig)); /
582		double amt = 0;
583		struct {
584		int s, d; /* source and destination */
#	Line 620 \| Line 640 \| *migration_step(MIGRATION mig, double src_rem, double*
640		return(best.amt);
641		}
642
643	<	/* Compute (and insert) migration along directed edge */
643	>	#ifdef DEBUG
644	>	static char *
645	>	thetaphi(const FVECT v)
646	>	{
647	>	static char buf[128];
648	>	double theta, phi;
649	>
650	>	theta = 180./M_PI*acos(v[2]);
651	>	phi = 180./M_PI*atan2(v[1],v[0]);
652	>	sprintf(buf, "(%.0f,%.0f)", theta, phi);
653	>
654	>	return(buf);
655	>	}
656	>	#endif
657	>
658	>	/* Create a new migration holder (sharing memory for multiprocessing) */
659		static MIGRATION *
660	<	make_migration(RBFNODE from_rbf, RBFNODE to_rbf)
660	>	new_migration(RBFNODE from_rbf, RBFNODE to_rbf)
661		{
662	+	size_t memlen = sizeof(MIGRATION) +
663	+	sizeof(float)(from_rbf->nrbfto_rbf->nrbf - 1);
664	+	MIGRATION *newmig;
665	+	#ifdef _WIN32
666	+	newmig = (MIGRATION *)malloc(memlen);
667	+	#else
668	+	if (nprocs <= 1) { /* single process? */
669	+	newmig = (MIGRATION *)malloc(memlen);
670	+	} else { /* else need to share memory */
671	+	newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ\|PROT_WRITE,
672	+	MAP_ANON\|MAP_SHARED, -1, 0);
673	+	if ((void *)newmig == MAP_FAILED)
674	+	newmig = NULL;
675	+	}
676	+	#endif
677	+	if (newmig == NULL) {
678	+	fprintf(stderr, "%s: cannot allocate new migration\n", progname);
679	+	exit(1);
680	+	}
681	+	newmig->rbfv[0] = from_rbf;
682	+	newmig->rbfv[1] = to_rbf;
683	+	/* insert in edge lists */
684	+	newmig->enxt[0] = from_rbf->ejl;
685	+	from_rbf->ejl = newmig;
686	+	newmig->enxt[1] = to_rbf->ejl;
687	+	to_rbf->ejl = newmig;
688	+	newmig->next = mig_list; /* push onto global list */
689	+	return(mig_list = newmig);
690	+	}
691	+
692	+	#ifdef _WIN32
693	+	#define await_children(n) (void)(n)
694	+	#define run_subprocess() 0
695	+	#define end_subprocess() (void)0
696	+	#else
697	+
698	+	/* Wait for the specified number of child processes to complete */
699	+	static void
700	+	await_children(int n)
701	+	{
702	+	if (n > nchild)
703	+	n = nchild;
704	+	while (n-- > 0) {
705	+	int status;
706	+	if (wait(&status) < 0) {
707	+	fprintf(stderr, "%s: missing child(ren)!\n", progname);
708	+	nchild = 0;
709	+	break;
710	+	}
711	+	--nchild;
712	+	if (status) {
713	+	if ((status = WEXITSTATUS(status)))
714	+	exit(status);
715	+	fprintf(stderr, "%s: subprocess died\n", progname);
716	+	exit(1);
717	+	}
718	+	}
719	+	}
720	+
721	+	/* Start child process if multiprocessing selected */
722	+	static pid_t
723	+	run_subprocess(void)
724	+	{
725	+	int status;
726	+	pid_t pid;
727	+
728	+	if (nprocs <= 1) /* any children requested? */
729	+	return(0);
730	+	await_children(nchild + 1 - nprocs); /* free up child process */
731	+	if ((pid = fork())) {
732	+	if (pid < 0) {
733	+	fprintf(stderr, "%s: cannot fork subprocess\n",
734	+	progname);
735	+	exit(1);
736	+	}
737	+	++nchild; /* subprocess started */
738	+	return(pid);
739	+	}
740	+	nchild = -1;
741	+	return(0); /* put child to work */
742	+	}
743	+
744	+	/* If we are in subprocess, call exit */
745	+	#define end_subprocess() if (nchild < 0) _exit(0); else
746	+
747	+	#endif /* ! _WIN32 */
748	+
749	+	/* Compute and insert migration along directed edge (may fork child) */
750	+	static MIGRATION *
751	+	create_migration(RBFNODE from_rbf, RBFNODE to_rbf)
752	+	{
753		const double end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
754	<	float *pmtx = price_routes(from_rbf, to_rbf);
755	<	MIGRATION newmig = (MIGRATION )malloc(sizeof(MIGRATION) +
756	<	sizeof(float) *
631	<	(from_rbf->nrbf*to_rbf->nrbf - 1));
632	<	double src_rem = (double )malloc(sizeof(double)*from_rbf->nrbf);
633	<	double dst_rem = (double )malloc(sizeof(double)*to_rbf->nrbf);
754	>	float *pmtx;
755	>	MIGRATION *newmig;
756	>	double src_rem, dst_rem;
757		double total_rem = 1.;
758		int i;
759	<
760	<	if ((newmig == NULL) \| (src_rem == NULL) \| (dst_rem == NULL)) {
761	<	fputs("Out of memory in make_migration()\n", stderr);
759	>	/* check if exists already */
760	>	for (newmig = from_rbf->ejl; newmig != NULL;
761	>	newmig = nextedge(from_rbf,newmig))
762	>	if (newmig->rbfv[1] == to_rbf)
763	>	return(NULL);
764	>	/* else allocate */
765	>	newmig = new_migration(from_rbf, to_rbf);
766	>	if (run_subprocess())
767	>	return(newmig); /* child continues */
768	>	pmtx = price_routes(from_rbf, to_rbf);
769	>	src_rem = (double )malloc(sizeof(double)from_rbf->nrbf);
770	>	dst_rem = (double )malloc(sizeof(double)to_rbf->nrbf);
771	>	if ((src_rem == NULL) \| (dst_rem == NULL)) {
772	>	fputs("Out of memory in create_migration()\n", stderr);
773		exit(1);
774		}
775		#ifdef DEBUG
776	<	{
777	<	double theta, phi;
778	<	theta = acos(from_rbf->invec[2])*(180./M_PI);
779	<	phi = atan2(from_rbf->invec[1],from_rbf->invec[0])*(180./M_PI);
646	<	fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
647	<	round(theta), round(phi));
648	<	theta = acos(to_rbf->invec[2])*(180./M_PI);
649	<	phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
650	<	fprintf(stderr, "(%d,%d)\n", round(theta), round(phi));
651	<	}
776	>	fprintf(stderr, "Building path from (theta,phi) %s ",
777	>	thetaphi(from_rbf->invec));
778	>	fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
779	>	/* if (nchild) */ fputc('\n', stderr);
780		#endif
653	–	newmig->next = NULL;
654	–	newmig->rbfv[0] = from_rbf;
655	–	newmig->rbfv[1] = to_rbf;
656	–	newmig->enxt[0] = newmig->enxt[1] = NULL;
657	–	memset(newmig->mtx, 0, sizeof(float)from_rbf->nrbfto_rbf->nrbf);
781		/* starting quantities */
782	+	memset(newmig->mtx, 0, sizeof(float)from_rbf->nrbfto_rbf->nrbf);
783		for (i = from_rbf->nrbf; i--; )
784		src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
785		for (i = to_rbf->nrbf; i--; )
786		dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
787		/* move a bit at a time */
788	<	while (total_rem > end_thresh)
788	>	while (total_rem > end_thresh) {
789		total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
790	+	#ifdef DEBUG
791	+	if (!nchild)
792	+	/* fputc('.', stderr); */
793	+	fprintf(stderr, "%.9f remaining...\r", total_rem);
794	+	#endif
795	+	}
796	+	#ifdef DEBUG
797	+	if (!nchild) fputs("done.\n", stderr);
798	+	#endif
799
800		free(pmtx); /* free working arrays */
801		free(src_rem);
#	Line 675 \| Line 808 \| *make_migration(RBFNODE from_rbf, RBFNODE to_rbf)*
808		for (j = to_rbf->nrbf; j--; )
809		newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
810		}
811	<	/* insert in edge lists */
812	<	newmig->enxt[0] = from_rbf->ejl;
680	<	from_rbf->ejl = newmig;
681	<	newmig->enxt[1] = to_rbf->ejl;
682	<	to_rbf->ejl = newmig;
683	<	newmig->next = mig_list;
684	<	return(mig_list = newmig);
811	>	end_subprocess(); /* exit here if subprocess */
812	>	return(newmig);
813		}
814
815		/* Get triangle surface orientation (unnormalized) */
#	Line 714 \| Line 842 \| *get_triangles(RBFNODE rbfv[2], const MIGRATION mig)*
842		RBFNODE *tv;
843
844		rbfv[0] = rbfv[1] = NULL;
845	+	if (mig == NULL)
846	+	return(0);
847		for (ej = mig->rbfv[0]->ejl; ej != NULL;
848		ej = nextedge(mig->rbfv[0],ej)) {
849		if (ej == mig)
#	Line 730 \| Line 860 \| *get_triangles(RBFNODE rbfv[2], const MIGRATION mig)*
860		return((rbfv[0] != NULL) + (rbfv[1] != NULL));
861		}
862
863	+	/* Check if prospective vertex would create overlapping triangle */
864	+	static int
865	+	overlaps_tri(const RBFNODE bv0, const RBFNODE bv1, const RBFNODE *pv)
866	+	{
867	+	const MIGRATION *ej;
868	+	RBFNODE *vother[2];
869	+	int im_rev;
870	+	/* find shared edge in mesh */
871	+	for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
872	+	const RBFNODE *tv = opp_rbf(pv,ej);
873	+	if (tv == bv0) {
874	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
875	+	ej->rbfv[1]->invec, bv1->invec);
876	+	break;
877	+	}
878	+	if (tv == bv1) {
879	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
880	+	ej->rbfv[1]->invec, bv0->invec);
881	+	break;
882	+	}
883	+	}
884	+	if (!get_triangles(vother, ej)) /* triangle on same side? */
885	+	return(0);
886	+	return(vother[im_rev] != NULL);
887	+	}
888	+
889		/* Find context hull vertex to complete triangle (oriented call) */
890		static RBFNODE *
891		find_chull_vert(const RBFNODE rbf0, const RBFNODE rbf1)
892		{
893	<	FVECT vmid, vor;
893	>	FVECT vmid, vejn, vp;
894		RBFNODE rbf, rbfbest = NULL;
895	<	double dprod2, bestdprod2 = 0.5;
895	>	double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
896
897	+	VSUB(vejn, rbf1->invec, rbf0->invec);
898		VADD(vmid, rbf0->invec, rbf1->invec);
899	<	if (normalize(vmid) == 0)
899	>	if (normalize(vejn) == 0 \|\| normalize(vmid) == 0)
900		return(NULL);
901		/* XXX exhaustive search */
902	+	/* Find triangle with minimum rotation from perpendicular */
903		for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
904		if ((rbf == rbf0) \| (rbf == rbf1))
905		continue;
906	<	tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
907	<	dprod2 = DOT(vor, vmid);
750	<	if (dprod2 <= FTINY)
906	>	tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
907	>	if (DOT(vp, vmid) <= FTINY)
908		continue; /* wrong orientation */
909	<	dprod2 *= dprod2 / DOT(vor,vor);
910	<	if (dprod2 > bestdprod2) { /* more convex than prev? */
911	<	rbfbest = rbf;
912	<	bestdprod2 = dprod2;
913	<	}
909	>	area2 = .25*DOT(vp,vp);
910	>	VSUB(vp, rbf->invec, rbf0->invec);
911	>	dprod = -DOT(vp, vejn);
912	>	VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
913	>	dprod = DOT(vp, vmid) / VLEN(vp);
914	>	if (dprod <= bestdprod + FTINY(1 - 2(area2 < bestarea2)))
915	>	continue; /* found better already */
916	>	if (overlaps_tri(rbf0, rbf1, rbf))
917	>	continue; /* overlaps another triangle */
918	>	rbfbest = rbf;
919	>	bestdprod = dprod; /* new one to beat */
920	>	bestarea2 = area2;
921		}
922	<	return(rbf);
922	>	return(rbfbest);
923		}
924
925		/* Create new migration edge and grow mesh recursively around it */
926		static void
927	<	mesh_from_edge(RBFNODE rbf0, RBFNODE rbf1)
927	>	mesh_from_edge(MIGRATION *edge)
928		{
929	<	MIGRATION *newej;
929	>	MIGRATION ej0, ej1;
930		RBFNODE *tvert[2];
931	<
932	<	if (rbf0 < rbf1) /* avoid migration loops */
933	<	newej = make_migration(rbf0, rbf1);
770	<	else
771	<	newej = make_migration(rbf1, rbf0);
931	>
932	>	if (edge == NULL)
933	>	return;
934		/* triangle on either side? */
935	<	get_triangles(tvert, newej);
936	<	if (tvert[0] == NULL) { /* recurse on new right edge */
937	<	tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
935	>	get_triangles(tvert, edge);
936	>	if (tvert[0] == NULL) { /* grow mesh on right */
937	>	tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
938		if (tvert[0] != NULL) {
939	<	mesh_from_edge(rbf0, tvert[0]);
940	<	mesh_from_edge(rbf1, tvert[0]);
939	>	if (tvert[0] > edge->rbfv[0])
940	>	ej0 = create_migration(edge->rbfv[0], tvert[0]);
941	>	else
942	>	ej0 = create_migration(tvert[0], edge->rbfv[0]);
943	>	if (tvert[0] > edge->rbfv[1])
944	>	ej1 = create_migration(edge->rbfv[1], tvert[0]);
945	>	else
946	>	ej1 = create_migration(tvert[0], edge->rbfv[1]);
947	>	mesh_from_edge(ej0);
948	>	mesh_from_edge(ej1);
949		}
950	<	}
951	<	if (tvert[1] == NULL) { /* recurse on new left edge */
782	<	tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
950	>	} else if (tvert[1] == NULL) { /* grow mesh on left */
951	>	tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
952		if (tvert[1] != NULL) {
953	<	mesh_from_edge(rbf0, tvert[1]);
954	<	mesh_from_edge(rbf1, tvert[1]);
953	>	if (tvert[1] > edge->rbfv[0])
954	>	ej0 = create_migration(edge->rbfv[0], tvert[1]);
955	>	else
956	>	ej0 = create_migration(tvert[1], edge->rbfv[0]);
957	>	if (tvert[1] > edge->rbfv[1])
958	>	ej1 = create_migration(edge->rbfv[1], tvert[1]);
959	>	else
960	>	ej1 = create_migration(tvert[1], edge->rbfv[1]);
961	>	mesh_from_edge(ej0);
962	>	mesh_from_edge(ej1);
963		}
964		}
965		}
966
967	+	#ifdef DEBUG
968	+	#include "random.h"
969	+	#include "bmpfile.h"
970	+	/* Hash pointer to byte value (must return 0 for NULL) */
971	+	static int
972	+	byte_hash(const void *p)
973	+	{
974	+	size_t h = (size_t)p;
975	+	h ^= (size_t)p >> 8;
976	+	h ^= (size_t)p >> 16;
977	+	h ^= (size_t)p >> 24;
978	+	return(h & 0xff);
979	+	}
980	+	/* Write out BMP image showing edges */
981	+	static void
982	+	write_edge_image(const char *fname)
983	+	{
984	+	BMPHeader *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
985	+	BMPWriter *wtr;
986	+	int i, j;
987	+
988	+	fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
989	+	hdr->compr = BI_RLE8;
990	+	for (i = 256; --i; ) { /* assign random color map */
991	+	hdr->palette[i].r = random() & 0xff;
992	+	hdr->palette[i].g = random() & 0xff;
993	+	hdr->palette[i].b = random() & 0xff;
994	+	/* reject dark colors */
995	+	i += (hdr->palette[i].r + hdr->palette[i].g +
996	+	hdr->palette[i].b < 128);
997	+	}
998	+	hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
999	+	/* open output */
1000	+	wtr = BMPopenOutputFile(fname, hdr);
1001	+	if (wtr == NULL) {
1002	+	free(hdr);
1003	+	return;
1004	+	}
1005	+	for (i = 0; i < GRIDRES; i++) { /* write scanlines */
1006	+	for (j = 0; j < GRIDRES; j++)
1007	+	wtr->scanline[j] = byte_hash(mig_grid[i][j]);
1008	+	if (BMPwriteScanline(wtr) != BIR_OK)
1009	+	break;
1010	+	}
1011	+	BMPcloseOutput(wtr); /* close & clean up */
1012	+	}
1013	+	#endif
1014	+
1015		/* Draw edge list into mig_grid array */
1016		static void
1017		draw_edges()
#	Line 830 \| Line 1055 \| draw_edges()
1055		if (nnull)
1056		fprintf(stderr, "Warning: %d of %d edges are null\n",
1057		nnull, ntot);
1058	+	#ifdef DEBUG
1059	+	write_edge_image("bsdf_edges.bmp");
1060	+	#endif
1061		}
1062
1063		/* Build our triangle mesh from recorded RBFs */
#	Line 837 \| Line 1065 \| static void
1065		build_mesh()
1066		{
1067		double best2 = M_PI*M_PI;
1068	<	RBFNODE rbf, rbf_near = NULL;
1068	>	RBFNODE *shrt_edj[2];
1069	>	RBFNODE rbf0, rbf1;
1070		/* check if isotropic */
1071		if (single_plane_incident) {
1072	<	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
1073	<	if (rbf->next != NULL)
1074	<	make_migration(rbf, rbf->next);
1072	>	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
1073	>	if (rbf0->next != NULL)
1074	>	create_migration(rbf0, rbf0->next);
1075	>	await_children(nchild);
1076		return;
1077		}
1078	<	/* find RBF nearest to head */
1079	<	if (dsf_list == NULL)
1080	<	return;
1081	<	for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
1082	<	double dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
1078	>	/* start w/ shortest edge */
1079	>	shrt_edj[0] = shrt_edj[1] = NULL;
1080	>	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
1081	>	for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
1082	>	double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
1083		if (dist2 < best2) {
1084	<	rbf_near = rbf;
1084	>	shrt_edj[0] = rbf0;
1085	>	shrt_edj[1] = rbf1;
1086		best2 = dist2;
1087		}
1088		}
1089	<	if (rbf_near == NULL) {
1090	<	fputs("Cannot find nearest point for first edge\n", stderr);
1089	>	if (shrt_edj[0] == NULL) {
1090	>	fputs("Cannot find shortest edge\n", stderr);
1091		exit(1);
1092		}
1093		/* build mesh from this edge */
1094	<	mesh_from_edge(dsf_list, rbf_near);
1094	>	if (shrt_edj[0] < shrt_edj[1])
1095	>	mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
1096	>	else
1097	>	mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
1098	>	/* complete migrations */
1099	>	await_children(nchild);
1100		/* draw edge list into grid */
1101		draw_edges();
1102		}
#	Line 901 \| Line 1137 \| *identify_tri(MIGRATION miga[3], unsigned char vmap[GR**
1137		return(1); /* this neighborhood done */
1138		}
1139
1140	+	/* Insert vertex in ordered list */
1141	+	static void
1142	+	insert_vert(RBFNODE *vlist, RBFNODE v)
1143	+	{
1144	+	int i, j;
1145	+
1146	+	for (i = 0; vlist[i] != NULL; i++) {
1147	+	if (v == vlist[i])
1148	+	return;
1149	+	if (v < vlist[i])
1150	+	break;
1151	+	}
1152	+	for (j = i; vlist[j] != NULL; j++)
1153	+	;
1154	+	while (j > i) {
1155	+	vlist[j] = vlist[j-1];
1156	+	--j;
1157	+	}
1158	+	vlist[i] = v;
1159	+	}
1160	+
1161	+	/* Sort triangle edges in standard order */
1162	+	static int
1163	+	order_triangle(MIGRATION *miga[3])
1164	+	{
1165	+	RBFNODE *vert[7];
1166	+	MIGRATION *ord[3];
1167	+	int i;
1168	+	/* order vertices, first */
1169	+	memset(vert, 0, sizeof(vert));
1170	+	for (i = 3; i--; ) {
1171	+	if (miga[i] == NULL)
1172	+	return(0);
1173	+	insert_vert(vert, miga[i]->rbfv[0]);
1174	+	insert_vert(vert, miga[i]->rbfv[1]);
1175	+	}
1176	+	/* should be just 3 vertices */
1177	+	if ((vert[3] == NULL) \| (vert[4] != NULL))
1178	+	return(0);
1179	+	/* identify edge 0 */
1180	+	for (i = 3; i--; )
1181	+	if (miga[i]->rbfv[0] == vert[0] &&
1182	+	miga[i]->rbfv[1] == vert[1]) {
1183	+	ord[0] = miga[i];
1184	+	break;
1185	+	}
1186	+	if (i < 0)
1187	+	return(0);
1188	+	/* identify edge 1 */
1189	+	for (i = 3; i--; )
1190	+	if (miga[i]->rbfv[0] == vert[1] &&
1191	+	miga[i]->rbfv[1] == vert[2]) {
1192	+	ord[1] = miga[i];
1193	+	break;
1194	+	}
1195	+	if (i < 0)
1196	+	return(0);
1197	+	/* identify edge 2 */
1198	+	for (i = 3; i--; )
1199	+	if (miga[i]->rbfv[0] == vert[0] &&
1200	+	miga[i]->rbfv[1] == vert[2]) {
1201	+	ord[2] = miga[i];
1202	+	break;
1203	+	}
1204	+	if (i < 0)
1205	+	return(0);
1206	+	/* reassign order */
1207	+	miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1208	+	return(1);
1209	+	}
1210	+
1211		/* Find edge(s) for interpolating the given incident vector */
1212		static int
1213		get_interp(MIGRATION *miga[3], const FVECT invec)
#	Line 926 \| Line 1233 \| *get_interp(MIGRATION miga[3], const FVECT invec)**
1233		{ /* else use triangle mesh */
1234		unsigned char floodmap[GRIDRES][(GRIDRES+7)/8];
1235		int pstart[2];
1236	+	RBFNODE *vother;
1237	+	MIGRATION *ej;
1238	+	int i;
1239
1240		pos_from_vec(pstart, invec);
1241		memset(floodmap, 0, sizeof(floodmap));
#	Line 936 \| Line 1246 \| *get_interp(MIGRATION miga[3], const FVECT invec)**
1246		return(0); /* should never happen */
1247		if (miga[1] == NULL)
1248		return(1); /* on edge */
1249	<	return(3); /* else in triangle */
1249	>	/* verify triangle */
1250	>	if (!order_triangle(miga)) {
1251	>	#ifdef DEBUG
1252	>	fputs("Munged triangle in get_interp()\n", stderr);
1253	>	#endif
1254	>	vother = NULL; /* find triangle from edge */
1255	>	for (i = 3; i--; ) {
1256	>	RBFNODE *tpair[2];
1257	>	if (get_triangles(tpair, miga[i]) &&
1258	>	(vother = tpair[ is_rev_tri(
1259	>	miga[i]->rbfv[0]->invec,
1260	>	miga[i]->rbfv[1]->invec,
1261	>	invec) ]) != NULL)
1262	>	break;
1263	>	}
1264	>	if (vother == NULL) { /* couldn't find 3rd vertex */
1265	>	#ifdef DEBUG
1266	>	fputs("No triangle in get_interp()\n", stderr);
1267	>	#endif
1268	>	return(0);
1269	>	}
1270	>	/* reassign other two edges */
1271	>	for (ej = vother->ejl; ej != NULL;
1272	>	ej = nextedge(vother,ej)) {
1273	>	RBFNODE *vorig = opp_rbf(vother,ej);
1274	>	if (vorig == miga[i]->rbfv[0])
1275	>	miga[(i+1)%3] = ej;
1276	>	else if (vorig == miga[i]->rbfv[1])
1277	>	miga[(i+2)%3] = ej;
1278	>	}
1279	>	if (!order_triangle(miga)) {
1280	>	#ifdef DEBUG
1281	>	fputs("Bad triangle in get_interp()\n", stderr);
1282	>	#endif
1283	>	return(0);
1284	>	}
1285	>	}
1286	>	return(3); /* return in standard order */
1287		}
1288		}
1289
#	Line 1015 \| Line 1362 \| memerr:
1362		return(NULL); /* pro forma return */
1363		}
1364
1018	–	/* Insert vertex in ordered list */
1019	–	static void
1020	–	insert_vert(RBFNODE *vlist, RBFNODE v)
1021	–	{
1022	–	int i, j;
1023	–
1024	–	for (i = 0; vlist[i] != NULL; i++) {
1025	–	if (v == vlist[i])
1026	–	return;
1027	–	if (v < vlist[i])
1028	–	break;
1029	–	}
1030	–	for (j = i; vlist[j] != NULL; j++)
1031	–	;
1032	–	while (j > i) {
1033	–	vlist[j] = vlist[j-1];
1034	–	--j;
1035	–	}
1036	–	vlist[i] = v;
1037	–	}
1038	–
1039	–	/* Sort triangle edges in standard order */
1040	–	static void
1041	–	order_triangle(MIGRATION *miga[3])
1042	–	{
1043	–	RBFNODE *vert[4];
1044	–	MIGRATION *ord[3];
1045	–	int i;
1046	–	/* order vertices, first */
1047	–	memset(vert, 0, sizeof(vert));
1048	–	for (i = 0; i < 3; i++) {
1049	–	insert_vert(vert, miga[i]->rbfv[0]);
1050	–	insert_vert(vert, miga[i]->rbfv[1]);
1051	–	}
1052	–	/* identify edge 0 */
1053	–	for (i = 0; i < 3; i++)
1054	–	if (miga[i]->rbfv[0] == vert[0] &&
1055	–	miga[i]->rbfv[1] == vert[1]) {
1056	–	ord[0] = miga[i];
1057	–	break;
1058	–	}
1059	–	/* identify edge 1 */
1060	–	for (i = 0; i < 3; i++)
1061	–	if (miga[i]->rbfv[0] == vert[1] &&
1062	–	miga[i]->rbfv[1] == vert[2]) {
1063	–	ord[1] = miga[i];
1064	–	break;
1065	–	}
1066	–	/* identify edge 2 */
1067	–	for (i = 0; i < 3; i++)
1068	–	if (miga[i]->rbfv[0] == vert[0] &&
1069	–	miga[i]->rbfv[1] == vert[2]) {
1070	–	ord[2] = miga[i];
1071	–	break;
1072	–	}
1073	–	miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1074	–	}
1075	–
1365		/* Partially advect between recorded incident angles and allocate new RBF */
1366		static RBFNODE *
1367		advect_rbf(const FVECT invec)
#	Line 1086 \| Line 1375 \| advect_rbf(const FVECT invec)
1375
1376		if (!get_interp(miga, invec)) /* can't interpolate? */
1377		return(NULL);
1378	<	if (miga[1] == NULL) /* along edge? */
1378	>	if (miga[1] == NULL) /* advect along edge? */
1379		return(e_advect_rbf(miga[0], invec));
1091	–	/* put in standard order */
1092	–	order_triangle(miga);
1380		#ifdef DEBUG
1381		if (miga[0]->rbfv[0] != miga[2]->rbfv[0] \|
1382		miga[0]->rbfv[1] != miga[1]->rbfv[0] \|
#	Line 1186 \| Line 1473 \| interp_isotropic()
1473		int ix, ox, oy;
1474		FVECT ivec, ovec;
1475		double bsdf;
1476	<
1476	>	#if DEBUG
1477	>	fprintf(stderr, "Writing isotropic order %d ", samp_order);
1478	>	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1479	>	else fputs("raw data\n", stderr);
1480	>	#endif
1481		if (pctcull >= 0) { /* begin output */
1482		sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1483		pctcull, samp_order);
1484		fflush(stdout);
1485		ofp = popen(cmd, "w");
1486		if (ofp == NULL) {
1487	<	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1487	>	fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
1488	>	progname);
1489		exit(1);
1490		}
1491		} else
#	Line 1220 \| Line 1512 \| interp_isotropic()
1512		}
1513		if (pctcull >= 0) { /* finish output */
1514		if (pclose(ofp)) {
1515	<	fprintf(stderr, "Error running '%s'\n", cmd);
1515	>	fprintf(stderr, "%s: error running '%s'\n",
1516	>	progname, cmd);
1517		exit(1);
1518		}
1519		} else {
#	Line 1240 \| Line 1533 \| interp_anisotropic()
1533		int ix, iy, ox, oy;
1534		FVECT ivec, ovec;
1535		double bsdf;
1536	<
1536	>	#if DEBUG
1537	>	fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1538	>	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1539	>	else fputs("raw data\n", stderr);
1540	>	#endif
1541		if (pctcull >= 0) { /* begin output */
1542		sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1543		pctcull, samp_order);
1544		fflush(stdout);
1545		ofp = popen(cmd, "w");
1546		if (ofp == NULL) {
1547	<	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1547	>	fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
1548	>	progname);
1549		exit(1);
1550		}
1551		} else
#	Line 1275 \| Line 1573 \| interp_anisotropic()
1573		}
1574		if (pctcull >= 0) { /* finish output */
1575		if (pclose(ofp)) {
1576	<	fprintf(stderr, "Error running '%s'\n", cmd);
1576	>	fprintf(stderr, "%s: error running '%s'\n",
1577	>	progname, cmd);
1578		exit(1);
1579		}
1580		} else
#	Line 1291 \| Line 1590 \| *main(int argc, char argv[])**
1590		double bsdf;
1591		int i;
1592
1593	<	progname = argv[0];
1594	<	if (argc > 2 && !strcmp(argv[1], "-t")) {
1595	<	pctcull = atoi(argv[2]);
1593	>	progname = argv[0]; /* get options */
1594	>	while (argc > 2 && argv[1][0] == '-') {
1595	>	switch (argv[1][1]) {
1596	>	case 'n':
1597	>	nprocs = atoi(argv[2]);
1598	>	break;
1599	>	case 't':
1600	>	pctcull = atoi(argv[2]);
1601	>	break;
1602	>	default:
1603	>	goto userr;
1604	>	}
1605		argv += 2; argc -= 2;
1606		}
1607	<	if (argc < 3) {
1608	<	fprintf(stderr,
1609	<	"Usage: %s [-t pctcull] meas1.dat meas2.dat .. > bsdf.xml\n",
1607	>	if (argc < 3)
1608	>	goto userr;
1609	>	#ifdef _WIN32
1610	>	if (nprocs > 1) {
1611	>	fprintf(stderr, "%s: multiprocessing not supported\n",
1612		progname);
1613		return(1);
1614		}
1615	+	#endif
1616		for (i = 1; i < argc; i++) { /* compile measurements */
1617		if (!load_pabopto_meas(argv[i]))
1618		return(1);
#	Line 1317 \| Line 1628 \| *main(int argc, char argv[])**
1628		interp_anisotropic();
1629		/* xml_epilogue(); /* finish XML output */
1630		return(0);
1631	+	userr:
1632	+	fprintf(stderr,
1633	+	"Usage: %s [-n nprocs][-t pctcull] meas1.dat meas2.dat .. > bsdf.xml\n",
1634	+	progname);
1635	+	return(1);
1636		}
1637		#else
1638		/* Test main produces a Radiance model from the given input file */

Diff Legend

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

Comparing ray/src/cv/pabopto2xml.c (file contents): Revision 2.12 by greg, Thu Sep 20 01:23:36 2012 UTC vs. Revision 2.16 by greg, Sun Oct 14 22:31:20 2012 UTC

Diff Legend

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.12 by greg, Thu Sep 20 01:23:36 2012 UTC vs.
Revision 2.16 by greg, Sun Oct 14 22:31:20 2012 UTC