MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
octree_tools.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#ifndef MATH_OCTREETOOLS_HEADER
11#define MATH_OCTREETOOLS_HEADER
12
13#include <limits>
14
15#include "math/defines.h"
16#include "math/vector.h"
17#include "math/matrix.h"
18#include "math/matrix_tools.h"
19#include "math/functions.h"
20
23
30template <typename T>
31bool
33 math::Vector<T,3> const& boxhalfsize, math::Vector<T,3> const& a,
34 math::Vector<T,3> const& b, math::Vector<T,3> const& c);
35
40template <typename T>
41bool
43 math::Vector<T,3> const& pos,
44 math::Vector<T,3> const& boxhalfsize);
45
50template <typename T>
51bool
52ray_box_overlap (math::Vector<T,3> const& origin, math::Vector<T,3> const& dir,
53 math::Vector<T,3> const& box_min, math::Vector<T,3> const& box_max);
54
61template <typename T>
62T
64 math::Vector<T,3> const& dir, math::Vector<T,3> const& a,
65 math::Vector<T,3> const& b, math::Vector<T,3> const& c,
66 float* bary = 0);
67
72template <typename T, int N>
73bool
74box_box_overlap (math::Vector<T,N> const& b1min, math::Vector<T,N> const& b1max,
75 math::Vector<T,N> const& b2min, math::Vector<T,N> const& b2max);
76
81template <typename T>
84 math::Vector<T,3> const& p2, math::Vector<T,3> const& d2);
85
90template <typename T, int N>
91bool
93 math::Vector<T, N> const& aabb_min, math::Vector<T, N> const& aabb_max);
94
95/* ------------------------ Implementation ------------------------ */
96
97#define AXISTEST_X01(a, b, fa, fb) { \
98 T p0 = a * v[0][1] - b * v[0][2]; \
99 T p2 = a * v[2][1] - b * v[2][2]; \
100 T min = std::min(p0, p2); T max = std::max(p0, p2); \
101 T rad = fa * boxhalfsize[1] + fb * boxhalfsize[2]; \
102 if (min > rad || max < -rad) return false; }
103
104#define AXISTEST_X2(a, b, fa, fb) { \
105 T p0 = a * v[0][1] - b * v[0][2]; \
106 T p1 = a * v[1][1] - b * v[1][2]; \
107 T min = std::min(p0, p1); T max = std::max(p0, p1); \
108 T rad = fa * boxhalfsize[1] + fb * boxhalfsize[2]; \
109 if (min > rad || max < -rad) return false; }
110
111#define AXISTEST_Y02(a, b, fa, fb) { \
112 T p0 = -a * v[0][0] + b * v[0][2]; \
113 T p2 = -a * v[2][0] + b * v[2][2]; \
114 T min = std::min(p0, p2); T max = std::max(p0, p2); \
115 T rad = fa * boxhalfsize[0] + fb * boxhalfsize[2]; \
116 if (min > rad || max < -rad) return false; }
117
118#define AXISTEST_Y1(a, b, fa, fb) { \
119 T p0 = -a * v[0][0] + b * v[0][2]; \
120 T p1 = -a * v[1][0] + b * v[1][2]; \
121 T min = std::min(p0, p1); T max = std::max(p0, p1); \
122 T rad = fa * boxhalfsize[0] + fb * boxhalfsize[2]; \
123 if (min > rad || max < -rad) return false; }
124
125#define AXISTEST_Z12(a, b, fa, fb) { \
126 T p1 = a * v[1][0] - b * v[1][1]; \
127 T p2 = a * v[2][0] - b * v[2][1]; \
128 T min = std::min(p1, p2); T max = std::max(p1, p2); \
129 T rad = fa * boxhalfsize[0] + fb * boxhalfsize[1]; \
130 if (min > rad || max < -rad) return false; }
131
132#define AXISTEST_Z0(a, b, fa, fb) { \
133 T p0 = a * v[0][0] - b * v[0][1]; \
134 T p1 = a * v[1][0] - b * v[1][1]; \
135 T min = std::min(p0, p1); T max = std::max(p0, p1); \
136 T rad = fa * boxhalfsize[0] + fb * boxhalfsize[1]; \
137 if (min > rad || max < -rad) return false; }
138
139template <typename T>
140bool
142 math::Vector<T,3> const& boxhalfsize, math::Vector<T,3> const& a,
143 math::Vector<T,3> const& b, math::Vector<T,3> const& c)
144{
145 math::Vector<T,3> v[3] = { a - boxcenter, b - boxcenter, c - boxcenter };
146 math::Vector<T,3> e[3] = { v[1] - v[0], v[2] - v[1], v[0] - v[2] };
147
148 {
149 T abs[3] = { std::abs(e[0][0]), std::abs(e[0][1]), std::abs(e[0][2]) };
150 AXISTEST_X01(e[0][2], e[0][1], abs[2], abs[1])
151 AXISTEST_Y02(e[0][2], e[0][0], abs[2], abs[0])
152 AXISTEST_Z12(e[0][1], e[0][0], abs[1], abs[0])
153 }
154
155 {
156 T abs[3] = { std::abs(e[1][0]), std::abs(e[1][1]), std::abs(e[1][2]) };
157 AXISTEST_X01(e[1][2], e[1][1], abs[2], abs[1])
158 AXISTEST_Y02(e[1][2], e[1][0], abs[2], abs[0])
159 AXISTEST_Z0 (e[1][1], e[1][0], abs[1], abs[0])
160 }
161
162 {
163 T abs[3] = { std::abs(e[2][0]), std::abs(e[2][1]), std::abs(e[2][2]) };
164 AXISTEST_X2 (e[2][2], e[2][1], abs[2], abs[1])
165 AXISTEST_Y1 (e[2][2], e[2][0], abs[2], abs[0])
166 AXISTEST_Z12(e[2][1], e[2][0], abs[1], abs[0])
167 }
168
169 for (int i = 0; i < 3; ++i)
170 {
171 T min = math::min(v[0][i], v[1][i], v[2][i]);
172 T max = math::max(v[0][i], v[1][i], v[2][i]);
173 if (min > boxhalfsize[i] || max < -boxhalfsize[i])
174 return false;
175 }
176
177 math::Vector<T,3> normal(e[0].cross(e[1])); // Not normalized
178 if (!plane_box_overlap(normal, v[0], boxhalfsize))
179 return false;
180
181 return true;
182}
183
184/* ---------------------------------------------------------------- */
185
186template <typename T>
187bool
189 math::Vector<T,3> const& pos,
190 math::Vector<T,3> const& boxhalfsize)
191{
192 math::Vector<T,3> vmin, vmax;
193 for (int q = 0; q < 3; ++q)
194 {
195 if (normal[q] > (T)0)
196 {
197 vmin[q] = -boxhalfsize[q] - pos[q];
198 vmax[q] = boxhalfsize[q] - pos[q];
199 }
200 else
201 {
202 vmin[q] = boxhalfsize[q] - pos[q];
203 vmax[q] = -boxhalfsize[q] - pos[q];
204 }
205 }
206
207 if (normal.dot(vmin) > (T)0)
208 return false;
209 if (normal.dot(vmax) >= (T)0)
210 return true;
211
212 return false;
213}
214
215/* ---------------------------------------------------------------- */
216
217/*
218 * Ray-box intersection using IEEE numerical properties to ensure that the
219 * test is both robust and efficient, as described in:
220 *
221 * Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley
222 * "An Efficient and Robust Ray-Box Intersection Algorithm"
223 * Journal of graphics tools, 10(1):49-54, 2005
224 *
225 */
226
227template <typename T>
228bool
230 math::Vector<T,3> const& box_min, math::Vector<T,3> const& box_max)
231{
232 math::Vector<T,3> box[2] = { box_min, box_max };
233 T idir[3] = { (T)1 / dir[0], (T)1 / dir[1], (T)1 / dir[2] };
234 bool sign[3] = { idir[0] < (T)0, idir[1] < (T)0, idir[2] < (T)0 };
235
236 float tmin = (box[sign[0]][0] - origin[0]) * idir[0];
237 float tmax = (box[1 - sign[0]][0] - origin[0]) * idir[0];
238 float tymin = (box[sign[1]][1] - origin[1]) * idir[1];
239 float tymax = (box[1 - sign[1]][1] - origin[1]) * idir[1];
240
241 if (tmin > tymax || tymin > tmax)
242 return false;
243 if (tymin > tmin)
244 tmin = tymin;
245 if (tymax < tmax)
246 tmax = tymax;
247
248 float tzmin = (box[sign[2]][2] - origin[2]) * idir[2];
249 float tzmax = (box[1 - sign[2]][2] - origin[2]) * idir[2];
250
251 if (tmin > tzmax || tzmin > tmax)
252 return false;
253 if (tzmin > tmin)
254 tmin = tzmin;
255 if (tzmax < tmax)
256 tmax = tzmax;
257
258 T const t0 = (T)0;
259 T const t1 = (T)1/(T)0;
260 return tmin < t1 && tmax > t0;
261}
262
263/* ---------------------------------------------------------------- */
264
265template <typename T>
266T
268 math::Vector<T,3> const& dir, math::Vector<T,3> const& a,
269 math::Vector<T,3> const& b, math::Vector<T,3> const& c,
270 float* bary)
271{
272 /* find vectors for two edges sharing vert0 */
273 math::Vector<T,3> edge1 = b - a;
274 math::Vector<T,3> edge2 = c - a;
275
276 /* begin calculating determinant - also used to calculate U parameter */
277 math::Vector<T,3> pvec = dir.cross(edge2);
278
279 /* if determinant is near zero, ray lies in plane of triangle */
280 T det = edge1.dot(pvec);
281 if (MATH_FLOAT_EQ(det, T(0)))
282 return T(0);
283 T inv_det = T(1) / det;
284
285 /* calculate distance from vert0 to ray origin */
286 math::Vector<T,3> tvec = origin - a;
287
288 /* calculate U parameter and test bounds */
289 T u = tvec.dot(pvec) * inv_det;
290 if (u < T(0) || u > T(1))
291 return (T)0;
292
293 /* prepare to test V parameter */
294 math::Vector<T,3> qvec = tvec.cross(edge1);
295
296 /* calculate V parameter and test bounds */
297 T v = dir.dot(qvec) * inv_det;
298 if (v < T(0) || u + v > T(1))
299 return T(0);
300
301 /* calculate t, ray intersects triangle */
302 T t = edge2.dot(qvec) * inv_det;
303
304 /* Since "0" is used to indicate failure, return minimum value. */
305 if (t == T(0))
306 return std::numeric_limits<T>::min();
307
308 if (bary != 0)
309 {
310 bary[0] = u;
311 bary[1] = v;
312 }
313
314 return t;
315}
316
317/* ---------------------------------------------------------------- */
318
319template <typename T, int N>
320bool
322 math::Vector<T,N> const& b2min, math::Vector<T,N> const& b2max)
323{
324 for (int i = 0; i < N; ++i)
325 if (b1min[i] > b2max[i] || b1max[i] < b2min[i])
326 return false;
327 return true;
328}
329
330/* ---------------------------------------------------------------- */
331
332template <typename T>
335 math::Vector<T,3> const& p2, math::Vector<T,3> const& d2)
336{
337 math::Vector<T,3> dx = d1.cross(d2);
339 m1[0] = p2[0] - p1[0]; m1[1] = d2[0]; m1[2] = dx[0];
340 m1[3] = p2[1] - p1[1]; m1[4] = d2[1]; m1[5] = dx[1];
341 m1[6] = p2[2] - p1[2]; m1[7] = d2[2]; m1[8] = dx[2];
342
344 m2[0] = p2[0] - p1[0]; m2[1] = d1[0]; m2[2] = dx[0];
345 m2[3] = p2[1] - p1[1]; m2[4] = d1[1]; m2[5] = dx[1];
346 m2[6] = p2[2] - p1[2]; m2[7] = d1[2]; m2[8] = dx[2];
347
350 return ret;
351}
352
353/* ---------------------------------------------------------------- */
354
355template <typename T, int N>
356bool
358 math::Vector<T, N> const& aabb_min, math::Vector<T, N> const& aabb_max)
359{
360 for (int i = 0; i < N; ++i)
361 if (point[i] < aabb_min[i] || point[i] > aabb_max[i])
362 return false;
363 return true;
364}
365
368
369#endif /* MATH_OCTREETOOLS_HEADER */
Matrix class for arbitrary dimensions and types.
Definition matrix.h:54
Vector class for arbitrary dimensions and types.
Definition vector.h:87
T square_norm(void) const
Computes the squared norm of the vector (much cheaper).
Definition vector.h:441
Vector< T, N > cross(Vector< T, N > const &other) const
Cross product between this and another vector.
Definition vector.h:549
T dot(Vector< T, N > const &other) const
Dot (or scalar) product between self and another vector.
Definition vector.h:542
#define MATH_GEOM_NAMESPACE_BEGIN
Definition defines.h:21
#define MATH_NAMESPACE_BEGIN
Definition defines.h:15
#define MATH_NAMESPACE_END
Definition defines.h:16
#define MATH_GEOM_NAMESPACE_END
Definition defines.h:22
#define MATH_FLOAT_EQ(x, v)
Definition defines.h:98
bool plane_box_overlap(math::Vector< T, 3 > const &normal, math::Vector< T, 3 > const &pos, math::Vector< T, 3 > const &boxhalfsize)
Returns true if the given plane (in Hesse normal form) and a box in the origin given by 'boxhalfsizes...
bool box_box_overlap(math::Vector< T, N > const &b1min, math::Vector< T, N > const &b1max, math::Vector< T, N > const &b2min, math::Vector< T, N > const &b2max)
Returns true if the given boxes overlap, false otherwise.
bool triangle_box_overlap(math::Vector< T, 3 > const &boxcenter, math::Vector< T, 3 > const &boxhalfsize, math::Vector< T, 3 > const &a, math::Vector< T, 3 > const &b, math::Vector< T, 3 > const &c)
Returns true if the given triangle and box overlap, false otherwise.
bool ray_box_overlap(math::Vector< T, 3 > const &origin, math::Vector< T, 3 > const &dir, math::Vector< T, 3 > const &box_min, math::Vector< T, 3 > const &box_max)
Returns true if the given ray intersects with the given axis-aligned bounding box,...
math::Vector< T, 2 > ray_ray_intersect(math::Vector< T, 3 > const &p1, math::Vector< T, 3 > const &d1, math::Vector< T, 3 > const &p2, math::Vector< T, 3 > const &d2)
Intersects the given rays with each other.
T ray_triangle_intersect(math::Vector< T, 3 > const &origin, math::Vector< T, 3 > const &dir, math::Vector< T, 3 > const &a, math::Vector< T, 3 > const &b, math::Vector< T, 3 > const &c, float *bary=0)
Intersects the given ray with the given triangle and returns the 't' parameter of the intersection po...
bool point_box_overlap(math::Vector< T, N > const &point, math::Vector< T, N > const &aabb_min, math::Vector< T, N > const &aabb_max)
Returns true if the given point overlaps with the axis-aligned box.
T const & max(T const &a, T const &b, T const &c)
Returns the maximum value of three arguments.
Definition functions.h:227
T const & min(T const &a, T const &b, T const &c)
Returns the minimum value of three arguments.
Definition functions.h:219
T matrix_determinant(Matrix< T, N, N > const &mat)
Calculates the determinant of the given matrix.
#define AXISTEST_Y02(a, b, fa, fb)
#define AXISTEST_Z12(a, b, fa, fb)
#define AXISTEST_Y1(a, b, fa, fb)
#define AXISTEST_X2(a, b, fa, fb)
#define AXISTEST_X01(a, b, fa, fb)
#define AXISTEST_Z0(a, b, fa, fb)