ViewVC Help

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

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

Comparing ray/src/cv/bsdfinterp.c (file contents):
Revision 2.2 by greg, Sat Oct 20 07:02:00 2012 UTC vs.
Revision 2.3 by greg, Tue Oct 23 05:10:42 2012 UTC

#	Line 13 | Line 13 | static const char RCSid[] = "$Id$";
13		#include <string.h>
14		#include <math.h>
15		#include "bsdfrep.h"
16	–	/* migration edges drawn in raster fashion */
17	–	MIGRATION               *mig_grid[GRIDRES][GRIDRES];
16
19	–	#ifdef DEBUG
20	–	#include "random.h"
21	–	#include "bmpfile.h"
22	–	/* Hash pointer to byte value (must return 0 for NULL) */
23	–	static int
24	–	byte_hash(const void *p)
25	–	{
26	–	size_t  h = (size_t)p;
27	–	h ^= (size_t)p >> 8;
28	–	h ^= (size_t)p >> 16;
29	–	h ^= (size_t)p >> 24;
30	–	return(h & 0xff);
31	–	}
32	–	/* Write out BMP image showing edges */
33	–	static void
34	–	write_edge_image(const char *fname)
35	–	{
36	–	BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
37	–	BMPWriter       *wtr;
38	–	int             i, j;
39	–
40	–	fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
41	–	hdr->compr = BI_RLE8;
42	–	for (i = 256; --i; ) {                  /* assign random color map */
43	–	hdr->palette[i].r = random() & 0xff;
44	–	hdr->palette[i].g = random() & 0xff;
45	–	hdr->palette[i].b = random() & 0xff;
46	–	/* reject dark colors */
47	–	i += (hdr->palette[i].r + hdr->palette[i].g +
48	–	hdr->palette[i].b < 128);
49	–	}
50	–	hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
51	–	/* open output */
52	–	wtr = BMPopenOutputFile(fname, hdr);
53	–	if (wtr == NULL) {
54	–	free(hdr);
55	–	return;
56	–	}
57	–	for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
58	–	for (j = 0; j < GRIDRES; j++)
59	–	wtr->scanline[j] = byte_hash(mig_grid[i][j]);
60	–	if (BMPwriteScanline(wtr) != BIR_OK)
61	–	break;
62	–	}
63	–	BMPcloseOutput(wtr);                    /* close & clean up */
64	–	}
65	–	#endif
66	–
67	–	/* Draw edge list into mig_grid array */
68	–	void
69	–	draw_edges(void)
70	–	{
71	–	int             nnull = 0, ntot = 0;
72	–	MIGRATION       *ej;
73	–	int             p0[2], p1[2];
74	–
75	–	memset(mig_grid, 0, sizeof(mig_grid));
76	–	for (ej = mig_list; ej != NULL; ej = ej->next) {
77	–	++ntot;
78	–	pos_from_vec(p0, ej->rbfv[0]->invec);
79	–	pos_from_vec(p1, ej->rbfv[1]->invec);
80	–	if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
81	–	++nnull;
82	–	mig_grid[p0[0]][p0[1]] = ej;
83	–	continue;
84	–	}
85	–	if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
86	–	const int       xstep = 2*(p1[0] > p0[0]) - 1;
87	–	const double    ystep = (double)((p1[1]-p0[1])*xstep) /
88	–	(double)(p1[0]-p0[0]);
89	–	int             x;
90	–	double          y;
91	–	for (x = p0[0], y = p0[1]+.5; x != p1[0];
92	–	x += xstep, y += ystep)
93	–	mig_grid[x][(int)y] = ej;
94	–	mig_grid[x][(int)y] = ej;
95	–	} else {
96	–	const int       ystep = 2*(p1[1] > p0[1]) - 1;
97	–	const double    xstep = (double)((p1[0]-p0[0])*ystep) /
98	–	(double)(p1[1]-p0[1]);
99	–	int             y;
100	–	double          x;
101	–	for (y = p0[1], x = p0[0]+.5; y != p1[1];
102	–	y += ystep, x += xstep)
103	–	mig_grid[(int)x][y] = ej;
104	–	mig_grid[(int)x][y] = ej;
105	–	}
106	–	}
107	–	if (nnull)
108	–	fprintf(stderr, "Warning: %d of %d edges are null\n",
109	–	nnull, ntot);
110	–	#ifdef DEBUG
111	–	write_edge_image("bsdf_edges.bmp");
112	–	#endif
113	–	}
114	–
115	–	/* Identify enclosing triangle for this position (flood fill raster check) */
116	–	static int
117	–	identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
118	–	int px, int py)
119	–	{
120	–	const int       btest = 1<<(py&07);
121	–
122	–	if (vmap[px][py>>3] & btest)            /* already visited here? */
123	–	return(1);
124	–	/* else mark it */
125	–	vmap[px][py>>3] |= btest;
126	–
127	–	if (mig_grid[px][py] != NULL) {         /* are we on an edge? */
128	–	int     i;
129	–	for (i = 0; i < 3; i++) {
130	–	if (miga[i] == mig_grid[px][py])
131	–	return(1);
132	–	if (miga[i] != NULL)
133	–	continue;
134	–	miga[i] = mig_grid[px][py];
135	–	return(1);
136	–	}
137	–	return(0);                      /* outside triangle! */
138	–	}
139	–	/* check neighbors (flood) */
140	–	if (px > 0 && !identify_tri(miga, vmap, px-1, py))
141	–	return(0);
142	–	if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
143	–	return(0);
144	–	if (py > 0 && !identify_tri(miga, vmap, px, py-1))
145	–	return(0);
146	–	if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
147	–	return(0);
148	–	return(1);                              /* this neighborhood done */
149	–	}
150	–
17		/* Insert vertex in ordered list */
18		static void
19		insert_vert(RBFNODE **vlist, RBFNODE *v)
#	Line 219 | Line 85 | order_triangle(MIGRATION *miga[3])
85		return(1);
86		}
87
88	+	/* Determine if we are close enough to the given edge */
89	+	static int
90	+	on_edge(const MIGRATION *ej, const FVECT ivec)
91	+	{
92	+	double  cos_a = DOT(ej->rbfv[0]->invec, ivec);
93	+	double  cos_b = DOT(ej->rbfv[1]->invec, ivec);
94	+	double  cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec);
95	+	double  cos_aplusb = cos_a*cos_b -
96	+	sqrt((1.-cos_a*cos_a)*(1.-cos_b*cos_b));
97	+
98	+	return(cos_aplusb - cos_c < .01);
99	+	}
100	+
101	+	/* Determine if we are inside the given triangle */
102	+	static int
103	+	in_tri(const RBFNODE *v1, const RBFNODE *v2, const RBFNODE *v3, const FVECT p)
104	+	{
105	+	FVECT   vc;
106	+	int     sgn1, sgn2, sgn3;
107	+	/* signed volume test */
108	+	VCROSS(vc, v1->invec, v2->invec);
109	+	sgn1 = (DOT(p, vc) > 0);
110	+	VCROSS(vc, v2->invec, v3->invec);
111	+	sgn2 = (DOT(p, vc) > 0);
112	+	if (sgn1 != sgn2)
113	+	return(0);
114	+	VCROSS(vc, v3->invec, v1->invec);
115	+	sgn3 = (DOT(p, vc) > 0);
116	+	return(sgn2 == sgn3);
117	+	}
118	+
119	+	/* Compute intersection with the given position over remaining mesh */
120	+	static int
121	+	in_mesh(MIGRATION *miga[3], unsigned char *emap, int nedges,
122	+	const FVECT ivec, MIGRATION *mig)
123	+	{
124	+	MIGRATION       *ej1, *ej2;
125	+	RBFNODE         *tv;
126	+	int             ejndx;
127	+	/* check visitation record */
128	+	if (mig->rbfv[0]->ord > mig->rbfv[1]->ord)
129	+	ejndx = mig->rbfv[1]->ord + (nedges-1)*mig->rbfv[0]->ord;
130	+	else
131	+	ejndx = mig->rbfv[0]->ord + (nedges-1)*mig->rbfv[1]->ord;
132	+	if (emap[ejndx>>3] & 1<<(ejndx&07))     /* tested already? */
133	+	return(0);
134	+	emap[ejndx>>3] |= 1<<(ejndx&07);        /* else mark & test it */
135	+	if (on_edge(mig, ivec)) {
136	+	miga[0] = mig;                  /* close enough to edge */
137	+	return(1);
138	+	}
139	+	/* do triangles either side */
140	+	for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
141	+	ej1 = nextedge(mig->rbfv[0],ej1)) {
142	+	if (ej1 == mig)
143	+	continue;
144	+	tv = opp_rbf(mig->rbfv[0],ej1);
145	+	for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
146	+	if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
147	+	if (in_mesh(miga, emap, nedges, ivec, ej1))
148	+	return(1);
149	+	if (in_mesh(miga, emap, nedges, ivec, ej2))
150	+	return(1);
151	+	if (in_tri(mig->rbfv[0], mig->rbfv[1],
152	+	tv, ivec)) {
153	+	miga[0] = mig;
154	+	miga[1] = ej1;
155	+	miga[2] = ej2;
156	+	return(1);
157	+	}
158	+	}
159	+	}
160	+	return(0);
161	+	}
162	+
163		/* Find edge(s) for interpolating the given vector, applying symmetry */
164		int
165		get_interp(MIGRATION *miga[3], FVECT invec)
#	Line 242 | Line 183 | get_interp(MIGRATION *miga[3], FVECT invec)
183		return(-1);                     /* outside range! */
184		}
185		{                                       /* else use triangle mesh */
186	<	const int       sym = use_symmetry(invec);
187	<	unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
188	<	int             pstart[2];
189	<	RBFNODE         *vother;
190	<	MIGRATION       *ej;
191	<	int             i;
192	<
193	<	pos_from_vec(pstart, invec);
194	<	memset(floodmap, 0, sizeof(floodmap));
195	<	/* call flooding function */
196	<	if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
197	<	return(-1);             /* outside mesh */
198	<	if ((miga[0] == NULL) | (miga[2] == NULL))
199	<	return(-1);             /* should never happen */
200	<	if (miga[1] == NULL)
201	<	return(sym);            /* on edge */
202	<	/* verify triangle */
203	<	if (!order_triangle(miga)) {
186	>	int             sym = use_symmetry(invec);
187	>	int             nedges = 0;
188	>	MIGRATION       *mep;
189	>	unsigned char   *emap;
190	>	/* clear visitation map */
191	>	for (mep = mig_list; mep != NULL; mep = mep->next)
192	>	++nedges;
193	>	emap = (unsigned char *)calloc((nedges*(nedges-1) + 7)>>3, 1);
194	>	if (emap == NULL) {
195	>	fprintf(stderr, "%s: Out of memory in get_interp()\n",
196	>	progname);
197	>	exit(1);
198	>	}
199	>	/* identify intersection  */
200	>	if (!in_mesh(miga, emap, nedges, invec, mig_list))
201	>	sym = -1;               /* outside mesh */
202	>	else if (miga[1] != NULL &&
203	>	(miga[2] == NULL || !order_triangle(miga))) {
204		#ifdef DEBUG
205		fputs("Munged triangle in get_interp()\n", stderr);
206		#endif
207	<	vother = NULL;          /* find triangle from edge */
267	<	for (i = 3; i--; ) {
268	<	RBFNODE     *tpair[2];
269	<	if (get_triangles(tpair, miga[i]) &&
270	<	(vother = tpair[ is_rev_tri(
271	<	miga[i]->rbfv[0]->invec,
272	<	miga[i]->rbfv[1]->invec,
273	<	invec) ]) != NULL)
274	<	break;
275	<	}
276	<	if (vother == NULL) {   /* couldn't find 3rd vertex */
277	<	#ifdef DEBUG
278	<	fputs("No triangle in get_interp()\n", stderr);
279	<	#endif
280	<	return(-1);
281	<	}
282	<	/* reassign other two edges */
283	<	for (ej = vother->ejl; ej != NULL;
284	<	ej = nextedge(vother,ej)) {
285	<	RBFNODE *vorig = opp_rbf(vother,ej);
286	<	if (vorig == miga[i]->rbfv[0])
287	<	miga[(i+1)%3] = ej;
288	<	else if (vorig == miga[i]->rbfv[1])
289	<	miga[(i+2)%3] = ej;
290	<	}
291	<	if (!order_triangle(miga)) {
292	<	#ifdef DEBUG
293	<	fputs("Bad triangle in get_interp()\n", stderr);
294	<	#endif
295	<	return(-1);
296	<	}
207	>	sym = -1;
208		}
209	+	free(emap);
210		return(sym);                    /* return in standard order */
211		}
212		}

Diff Legend

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