Subversion Repositories gelsvn

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
290 jrf 1 
#include <cfloat>
2 
#include <GL/gl.h>
3 
#include "CGLA/statistics.h"
4 
#include "CGLA/eigensolution.h"
5 
#include "CGLA/Mat4x4f.h"
6 
#include "AABox.h"
7 
#include "OBox.h"
8 
 
9 
using namespace std;
10 
using namespace CGLA;
11 
 
12 
namespace
13 
{
14 
 
15 
	Mat3x3f compute_rotation(const vector<Triangle>& invec)
16 
	{
17 
		const int N_tri = invec.size();
18 
 
19 
		Mat3x3f C;
20 
		float a_H = 0;
21 
		Vec3f m_H(0);
22 
 
23 
		for(int i=0;i<N_tri;++i)
24 
			{
25 
				const Triangle& tri = invec[i];
26 
 
27 
				float a_k = tri.area();
28 
				a_H += a_k;
29 
 
30 
				Vec3f m_k = tri.centre();
31 
				m_H += a_k * m_k;
32 
 
33 
				Vec3f p_k = tri.get_v0();
34 
				Vec3f q_k = tri.get_v1();
35 
				Vec3f r_k = tri.get_v2();
36 
 
37 
				Mat3x3f M,P,Q,R;
38 
				outer_product(m_k,m_k, M);
39 
				outer_product(p_k,p_k, P);
40 
				outer_product(q_k,q_k, Q);
41 
				outer_product(r_k,r_k, R);
42 
 
43 
				C += (a_k/12.0f) * (9*M+P+Q+R);
44 
			}
45 
		m_H /= a_H;
46 
		C /= a_H;
47 
		Mat3x3f M;
48 
		outer_product(m_H, m_H, M);
49 
		C -= M;
50 
 
51 
		Mat3x3f Q,L;
52 
		const int n_eig = power_eigensolution(C,Q,L,2);
53 
 
54 
		Vec3f X = normalize(Q[0]);
55 
		Vec3f Y = normalize(Q[1]);
56 
		Vec3f Z;
57 
 
58 
		float xy_ortho = fabs(dot(X,Y));
59 
		if(n_eig == 2 && xy_ortho < 0.3)
60 
			{
61 
				if(xy_ortho > CGLA::TINY)
62 
					Y = normalize(Y-X*dot(X,Y));
63 
				Z = normalize(cross(X,Y));
64 
			}
65 
		else if(n_eig==1)
66 
			{
67 
				Y = Vec3f(X[2],X[0],X[1]);
68 
				Y = normalize(Y-X*dot(X,Y));
69 
				Z = normalize(cross(X,Y));
70 
			}
71 
		else
72 
			{
73 
				X=Vec3f(1,0,0);
74 
				Y=Vec3f(0,1,0);
75 
				Z=Vec3f(0,0,1);
76 
			}
77 
		return Mat3x3f(X,Y,Z);
78 
 
79 
	}
80 
 
81 
 
82 
 
83 
}
84 
 
85 
bool OBox::intersect(const CGLA::Vec3f& p, const CGLA::Vec3f& d) const
86 
{
87 
	Vec3f pr = R * p;
88 
	Vec3f dr = R * d;
89 
	return aabox.intersect(pr, dr);
90 
}
91 
 
92 
void OBox::minmax_sq_dist(const CGLA::Vec3f& p,
93 
													float& dmin, float& dmax) const
94 
{
95 
	Vec3f pr = R * p;
96 
	aabox.minmax_sq_dist(pr, dmin, dmax);
97 
}
98 
 
99 
void OBox::gl_draw() const
100 
{
101 
	Mat4x4f m = identity_Mat4x4f();
102 
	copy_matrix(R, m);
103 
	glPushMatrix();
104 
	glMultMatrixf(m.get());
105 
	aabox.gl_draw();
106 
	glPopMatrix();
107 
}
108 
 
109 
OBox OBox::box_triangle(const Triangle& t)
110 
{
111 
	Vec3f e0 = t.get_v1()-t.get_v0();
112 
	Vec3f e1 = t.get_v2()-t.get_v1();
113 
	Vec3f e2 = t.get_v0()-t.get_v2();
114 
 
115 
	Vec3f X,Y,Z;
116 
	if(sqr_length(e0) > sqr_length(e1))
117 
		{
118 
			if(sqr_length(e0) > sqr_length(e2))
119 
				{
120 
					X = normalize(e0);
121 
					Y = normalize(e1 - X * dot(X, e1));
122 
				}
123 
			else
124 
				{
125 
					X = normalize(e2);
126 
					Y = normalize(e0 - X * dot(X, e0));
127 
				}
128 
		}
129 
	else
130 
		{
131 
			if(sqr_length(e1) > sqr_length(e2))
132 
				{
133 
					X = normalize(e1);
134 
					Y = normalize(e2 - X * dot(X, e2));
135 
				}
136 
			else
137 
				{
138 
					X = normalize(e2);
139 
					Y = normalize(e0 - X * dot(X, e0));
140 
				}
141 
		}
142 
	Z = cross(X,Y);
143 
 
144 
	const Mat3x3f Rot(X,Y,Z);
145 
 
146 
	Vec3f p0 = Rot * t.get_v0();
147 
	Vec3f p1 = Rot * t.get_v1();
148 
	Vec3f p2 = Rot * t.get_v2();
149 
	Vec3f pmin = v_min(p0, v_min(p1, p2));
150 
	Vec3f pmax = v_max(p0, v_max(p1, p2));
151 
 
152 
	Vec3f centre_close = v_max(pmin, v_min(pmax, Rot * t.get_centre()));
153 
	return OBox(Rot, AABox(pmin, pmax, centre_close));
154 
}
155 
 
156 
 
157 
OBox OBox::box_and_split(const std::vector<Triangle>& invec,
158 
													 std::vector<Triangle>& lvec,
159 
													 std::vector<Triangle>& rvec)
160 
{
161 
	// Obtain the rotation matrix for the OBB
162 
	const Mat3x3f Rot = compute_rotation(invec);
163 
	const int N_tri = invec.size();
164 
	const int N_pts = 3*N_tri;
165 
 
166 
	// Compute the rotated set of points and the extents of the point aligned
167 
	// BBox.
168 
	vector<Vec3f> pts(N_pts);
169 
	Vec3f tri_pmin(FLT_MAX), tri_pmax(-FLT_MAX);
170 
	for(int i=0;i<N_tri;++i)
171 
		{
172 
			const Triangle& tri = invec[i];
173 
 
174 
			int offs = 3*i;
175 
			pts[offs  ] = Rot*tri.get_v0();
176 
			pts[offs+1] = Rot*tri.get_v1();
177 
			pts[offs+2] = Rot*tri.get_v2();
178 
 
179 
			for(int j=0;j<3;++j)
180 
				{
181 
					tri_pmin = v_min(pts[offs+j], tri_pmin);
182 
					tri_pmax = v_max(pts[offs+j], tri_pmax);
183 
				}
184 
		}
185 
 
186 
	// Find the point closest to the centre.
187 
	const Vec3f centre = tri_pmin + 0.5f*(tri_pmax-tri_pmin);
188 
	Vec3f centre_close;
189 
	float min_dist = FLT_MAX;
190 
	for(int i=0;i<N_pts;++i)
191 
		{
192 
			Vec3f v = pts[i];
193 
			float sl = sqr_length(centre-v);
194 
			if(sl < min_dist)
195 
				{
196 
					min_dist = sl;
197 
					centre_close = v;
198 
				}
199 
		}
200 
 
201 
	// Partition the triangles
202 
	const float thresh = centre[0];
203 
	for(int i=0;i<N_tri;++i)
204 
		{
205 
			Vec3f p = Rot * invec[i].get_centre();
206 
			if( p[0] > thresh)
207 
				rvec.push_back(invec[i]);
208 
			else
209 
				lvec.push_back(invec[i]);
210 
		}
211 
 
212 
	// If all triangles landed in one box, split them naively.
213 
	if(lvec.empty() || rvec.empty())
214 
		{
215 
			lvec.clear();
216 
			lvec.insert(lvec.end(),
217 
									invec.begin(),
218 
									invec.begin()+N_tri/2);
219 
			rvec.clear();
220 
			rvec.insert(rvec.end(),
221 
									invec.begin()+N_tri/2,
222 
									invec.end());
223 
		}
224 
 
225 
	return OBox(Rot, AABox(tri_pmin, tri_pmax, centre_close));
226 
}
227