Subversion Repositories gelsvn

Rev

Rev 330 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
290 jrf 1 
#include <cfloat>
2 
#include "CGLA/statistics.h"
3 
#include "CGLA/eigensolution.h"
4 
#include "CGLA/Mat4x4f.h"
5 
#include "AABox.h"
6 
#include "OBox.h"
291 jrf 7 
#include "Triangle.h"
290 jrf 8 
 
9 
using namespace std;
10 
using namespace CGLA;
11 
 
12 
namespace
13 
{
14 
 
291 jrf 15 
	Mat3x3f compute_rotation(const vector<Geometry::Triangle>& invec)
290 jrf 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 
			{
291 jrf 25 
				const Geometry::Triangle& tri = invec[i];
290 jrf 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 
 
291 jrf 85 
namespace Geometry
86 
{
87 
 
290 jrf 88 
bool OBox::intersect(const CGLA::Vec3f& p, const CGLA::Vec3f& d) const
89 
{
90 
	Vec3f pr = R * p;
91 
	Vec3f dr = R * d;
92 
	return aabox.intersect(pr, dr);
93 
}
94 
 
95 
 
96 
OBox OBox::box_triangle(const Triangle& t)
97 
{
98 
	Vec3f e0 = t.get_v1()-t.get_v0();
99 
	Vec3f e1 = t.get_v2()-t.get_v1();
100 
	Vec3f e2 = t.get_v0()-t.get_v2();
101 
 
102 
	Vec3f X,Y,Z;
103 
	if(sqr_length(e0) > sqr_length(e1))
104 
		{
105 
			if(sqr_length(e0) > sqr_length(e2))
106 
				{
107 
					X = normalize(e0);
108 
					Y = normalize(e1 - X * dot(X, e1));
109 
				}
110 
			else
111 
				{
112 
					X = normalize(e2);
113 
					Y = normalize(e0 - X * dot(X, e0));
114 
				}
115 
		}
116 
	else
117 
		{
118 
			if(sqr_length(e1) > sqr_length(e2))
119 
				{
120 
					X = normalize(e1);
121 
					Y = normalize(e2 - X * dot(X, e2));
122 
				}
123 
			else
124 
				{
125 
					X = normalize(e2);
126 
					Y = normalize(e0 - X * dot(X, e0));
127 
				}
128 
		}
129 
	Z = cross(X,Y);
130 
 
131 
	const Mat3x3f Rot(X,Y,Z);
132 
 
133 
	Vec3f p0 = Rot * t.get_v0();
134 
	Vec3f p1 = Rot * t.get_v1();
135 
	Vec3f p2 = Rot * t.get_v2();
136 
	Vec3f pmin = v_min(p0, v_min(p1, p2));
137 
	Vec3f pmax = v_max(p0, v_max(p1, p2));
138 
 
139 
	Vec3f centre_close = v_max(pmin, v_min(pmax, Rot * t.get_centre()));
140 
	return OBox(Rot, AABox(pmin, pmax, centre_close));
141 
}
142 
 
143 
 
144 
OBox OBox::box_and_split(const std::vector<Triangle>& invec,
145 
													 std::vector<Triangle>& lvec,
146 
													 std::vector<Triangle>& rvec)
147 
{
148 
	// Obtain the rotation matrix for the OBB
149 
	const Mat3x3f Rot = compute_rotation(invec);
150 
	const int N_tri = invec.size();
151 
	const int N_pts = 3*N_tri;
152 
 
153 
	// Compute the rotated set of points and the extents of the point aligned
154 
	// BBox.
155 
	vector<Vec3f> pts(N_pts);
156 
	Vec3f tri_pmin(FLT_MAX), tri_pmax(-FLT_MAX);
157 
	for(int i=0;i<N_tri;++i)
158 
		{
159 
			const Triangle& tri = invec[i];
160 
 
161 
			int offs = 3*i;
162 
			pts[offs  ] = Rot*tri.get_v0();
163 
			pts[offs+1] = Rot*tri.get_v1();
164 
			pts[offs+2] = Rot*tri.get_v2();
165 
 
166 
			for(int j=0;j<3;++j)
167 
				{
168 
					tri_pmin = v_min(pts[offs+j], tri_pmin);
169 
					tri_pmax = v_max(pts[offs+j], tri_pmax);
170 
				}
171 
		}
172 
 
173 
	// Find the point closest to the centre.
174 
	const Vec3f centre = tri_pmin + 0.5f*(tri_pmax-tri_pmin);
175 
	Vec3f centre_close;
176 
	float min_dist = FLT_MAX;
177 
	for(int i=0;i<N_pts;++i)
178 
		{
179 
			Vec3f v = pts[i];
180 
			float sl = sqr_length(centre-v);
181 
			if(sl < min_dist)
182 
				{
183 
					min_dist = sl;
184 
					centre_close = v;
185 
				}
186 
		}
187 
 
188 
	// Partition the triangles
189 
	const float thresh = centre[0];
190 
	for(int i=0;i<N_tri;++i)
191 
		{
192 
			Vec3f p = Rot * invec[i].get_centre();
193 
			if( p[0] > thresh)
194 
				rvec.push_back(invec[i]);
195 
			else
196 
				lvec.push_back(invec[i]);
197 
		}
198 
 
199 
	// If all triangles landed in one box, split them naively.
200 
	if(lvec.empty() || rvec.empty())
201 
		{
202 
			lvec.clear();
203 
			lvec.insert(lvec.end(),
204 
									invec.begin(),
205 
									invec.begin()+N_tri/2);
206 
			rvec.clear();
207 
			rvec.insert(rvec.end(),
208 
									invec.begin()+N_tri/2,
209 
									invec.end());
210 
		}
211 
 
212 
	return OBox(Rot, AABox(tri_pmin, tri_pmax, centre_close));
213 
}
214 
 
330 jab 215 
}