MorphoGraphX
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Geometry.hpp
Go to the documentation of this file.
1 #ifndef GEOMETRY_H
2 #define GEOMETRY_H
3 
10 #include <Config.hpp>
11 
12 #include <Util.hpp>
13 #include <Vector.hpp>
14 #include <Matrix.hpp>
15 #include <BoundingBox.hpp>
16 
17 #include <cmath>
18 #include <cuda/CudaGlobal.hpp>
19 #include <float.h>
20 #include <math.h>
21 
22 namespace mgx {
23 
30  using util::min;
31  using util::max;
32  using util::trim;
33  typedef unsigned char ubyte;
34  typedef unsigned int uint;
35  typedef unsigned short ushort;
36  typedef unsigned long ulong;
37 
38  typedef util::Vector<2, float> Point2f;
39  typedef util::Vector<3, float> Point3f;
40  typedef util::Vector<4, float> Point4f;
41  typedef util::Vector<5, float> Point5f;
42  typedef util::Vector<6, float> Point6f;
43  typedef util::Vector<12, float> Point12f;
44 
45  typedef util::Vector<2, double> Point2d;
46  typedef util::Vector<3, double> Point3d;
47  typedef util::Vector<4, double> Point4d;
48  typedef util::Vector<5, double> Point5d;
49  typedef util::Vector<6, double> Point6d;
50  typedef util::Vector<16, double> Point16d;
51 
52  typedef util::Vector<2, int> Point2i;
53  typedef util::Vector<3, int> Point3i;
54  typedef util::Vector<4, int> Point4i;
55  typedef util::Vector<5, int> Point5i;
56  typedef util::Vector<6, int> Point6i;
57 
58  typedef util::Vector<2, uint> Point2u;
59  typedef util::Vector<3, uint> Point3u;
60  typedef util::Vector<4, uint> Point4u;
61  typedef util::Vector<5, uint> Point5u;
62  typedef util::Vector<6, uint> Point6u;
63 
64  typedef util::Vector<2, size_t> Point2s;
65  typedef util::Vector<3, size_t> Point3s;
66  typedef util::Vector<4, size_t> Point4s;
67  typedef util::Vector<5, size_t> Point5s;
68  typedef util::Vector<6, size_t> Point6s;
69 
70  typedef util::Vector<3, ushort> Point3us;
71 
72  typedef util::Matrix<2, 2, float> Matrix2f;
73  typedef util::Matrix<3, 3, float> Matrix3f;
74  typedef util::Matrix<4, 4, float> Matrix4f;
75 
76  typedef util::Matrix<2, 2, double> Matrix2d;
77  typedef util::Matrix<3, 3, double> Matrix3d;
78  typedef util::Matrix<4, 4, double> Matrix4d;
79 
80  typedef util::BoundingBox<3, uint> BoundingBox3u;
81  typedef util::BoundingBox<3, int> BoundingBox3i;
82  typedef util::BoundingBox<3, float> BoundingBox3f;
83 
84  CU_HOST_DEVICE
85  Point3u toVoxelsCeil(const Point3f& p, const Point3f& step)
86  {
87  Point3u r(0, 0, 0);
88  if(p.x() > 0 and step.x() > 0)
89  r.x() = uint(ceil(p.x() / step.x()));
90  if(p.y() > 0 and step.y() > 0)
91  r.y() = uint(ceil(p.y() / step.y()));
92  if(p.z() > 0 and step.z() > 0)
93  r.z() = uint(ceil(p.z() / step.z()));
94  return r;
95  }
96 
97  // Multiple openGL 4D matrix by 3D point and return 3D point
98  template <typename T>
99  CU_HOST_DEVICE util::Vector<3, T> multMatrix4Point3(const util::Matrix<4, 4, T>& m, const util::Vector<3, T>& p)
100  {
101  util::Vector<4, T> t(p.x(), p.y(), p.z(), T(1.0));
102  t = m * t;
103  util::Vector<3, T> res(t.x() / t.t(), t.y() / t.t(), t.z() / t.t());
104  return res;
105  }
106 
107  template <typename T> T CU_HOST_DEVICE interpolate(const T a, const T b, const T s)
108  {
109  T t = trim(s, (T)0.0, (T)1.0);
110  return ((1.0 - t) * a + t * b);
111  }
112 
113  // CU_HOST_DEVICE
114  // Point3f getEulerAngles(Matrix4d &m)
115  //{
116  // Point3f r(0,0,0);
117  // if(fabs(m(1,0)) > .1)
118  // r.x() = atan2(m(0,2), m(2,2)) / M_PI * 180;
119  // else
120  // r.x() = atan2(-m(2,0), m(0,0)) / M_PI * 180;
121  // r.y() = atan2(-m(1,2), m(2,2)) / M_PI * 180;
122  // r.z() = asin(m(1,0)) / M_PI * 180;
123  // return(r);
124  //}
125 
126  // Plane-Line Intersection, 0 <= s <= 1, on line, s >= 0 ray intersect
127  CU_HOST_DEVICE
128  bool planeLineIntersect(Point3f p, Point3f n, Point3f u1, Point3f u2, float& s, Point3f& u)
129  {
130  n.normalize();
131  Point3f u1u2 = u2 - u1;
132  float t = n * u1u2;
133  // Check if parallel
134  if(t == 0)
135  return (false);
136  s = n * (p - u1) / t;
137  u = u1 + s * u1u2;
138  return (true);
139  }
140 
141  // Triangle area
142  CU_HOST_DEVICE
143  float triangleArea(Point3f a, Point3f b, Point3f c) {
144  return ((a - b).cross(a - c).norm() / 2.0);
145  }
146 
147  // Signed volume of a tetrahedron
148  CU_HOST_DEVICE
149  float signedTetraVolume(Point3f a, Point3f b, Point3f c)
150  {
151  return ((-c.x() * b.y() * a.z() + b.x() * c.y() * a.z() + c.x() * a.y() * b.z() - a.x() * c.y() * b.z()
152  - b.x() * a.y() * c.z() + a.x() * b.y() * c.z()) / 6.0f);
153  }
154 
155  // Get (a set of) basis vectors from a plane
156  CU_HOST_DEVICE
157  void getBasisFromPlane(const Point3f& nrml, Point3f& x, Point3f& y, Point3f& z)
158  {
159  Point3f n = nrml;
160  n.normalize();
161  x = Point3f(1, 0, 0);
162  y = Point3f(0, 1, 0);
163  z = Point3f(0, 0, 1);
164  if(z.cross(n).norm() != 0) {
165  y = z.cross(n);
166  y.normalize();
167  x = y.cross(n);
168  x.normalize();
169  } else
170  x = n * z * x;
171  z = n;
172  }
173 
174  // Intersect a ray with a 3D triangle
175  // Input: ray r1,r2 triangle t0,t1,t2
176  // Output: intp = intersection point (when it exists)
177  // Return: -1 = triangle is degenerate (a segment or point)
178  // 0 = disjoint (no intersect)
179  // 1 = intersect in unique point I1
180  // 2 = are in the same plane
181  CU_HOST_DEVICE
182  int rayTriangleIntersect(const Point3f& r0, const Point3f& r1, const Point3f& t0, const Point3f& t1, const Point3f& t2,
183  Point3f& intp)
184  {
185  // get triangle edge vectors and plane normal
186  Point3f u = t1 - t0;
187  Point3f v = t2 - t0;
188  Point3f n = u.cross(v);
189  if(n.x() == 0 && n.y() == 0 && n.z() == 0) // triangle is degenerate
190  return -1;
191 
192  Point3f dir = r1 - r0; // ray direction vector
193  Point3f w0 = r0 - t0;
194  float a = -n * w0;
195  float b = n * dir;
196  if(b == 0) { // ray is parallel to triangle plane
197  if(a == 0) // ray lies in triangle plane
198  return 2;
199  else
200  return 0; // ray disjoint from plane
201  }
202 
203  // get intersect point of ray with triangle plane
204  float r = a / b;
205  if(r < 0.0f) // ray goes away from triangle
206  return 0; // => no intersect
207  // for a segment, also test if (r > 1.0) => no intersect
208 
209  intp = r0 + r * dir; // intersect point of ray and plane
210 
211  // is intersect point inside the triangle?
212  float uu, uv, vv, wu, wv, d;
213  uu = u * u;
214  uv = u * v;
215  vv = v * v;
216  Point3f w = intp - t0;
217  wu = w * u;
218  wv = w * v;
219  d = uv * uv - uu * vv;
220 
221  // get and test parametric coords
222  float s, t;
223  s = (uv * wv - vv * wu) / d;
224  if(s < 0.0f || s > 1.0f) // outside
225  return 0;
226  t = (uv * wu - uu * wv) / d;
227  if(t < 0.0f || (s + t) > 1.0f) // outside
228  return 0;
229 
230  return 1; // in
231  }
232 
233  // Return distance from a point to a line (or segment)
234  CU_HOST_DEVICE
235  float distLinePoint(Point3f v1, Point3f v2, Point3f p, bool segment)
236  {
237  v2 -= v1;
238  float n = norm(v2);
239  // Points are the same
240  if(n == 0)
241  return norm(p - v1);
242 
243  Point3f vn = v2 / n;
244  p -= v1;
245  Point3f q = vn * (p * vn);
246  float dist = norm(p - q);
247  if(!segment)
248  return dist;
249 
250  // Projected point on line segment
251  float qn = p * vn;
252  if(qn >= 0 and qn <= n)
253  return dist;
254 
255  // Otherwise pick closest endpoint
256  dist = norm(p);
257  float d = norm(p - v2);
258  if(d < dist)
259  dist = d;
260 
261  return dist;
262  }
263 
264 }
265 #endif
CU_HOST_DEVICE Vector cross(const Vector &other) const
Compute the cross product as this x other.
Definition: Vector.hpp:641
CU_HOST_DEVICE Vector & normalize(void)
Normalize the vector.
Definition: Vector.hpp:548
Common definitions and utilities This file is shared by cuda, do not include headers that nvcc can't ...
Defines the Matrix class template This file is shared by cuda, do not include headers that nvcc can't...
Defines the Vector class template This file is shared by cuda, do not include headers that nvcc can't...