Subversion Repositories gelsvn

Rev

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

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