VTK-m  2.0
TriangleIntersections.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 #ifndef vtk_m_rendering_raytracing_TriangleIntersections_h
11 #define vtk_m_rendering_raytracing_TriangleIntersections_h
12 
13 #include <vtkm/Math.h>
14 
15 namespace vtkm
16 {
17 namespace rendering
18 {
19 namespace raytracing
20 {
21 
22 class Moller
23 {
24 public:
25  template <typename Precision>
26  VTKM_EXEC void IntersectTri(const vtkm::Vec<Precision, 3>& a,
27  const vtkm::Vec<Precision, 3>& b,
28  const vtkm::Vec<Precision, 3>& c,
29  const vtkm::Vec<Precision, 3>& dir,
30  Precision& distance,
31  Precision& u,
32  Precision& v,
33  const vtkm::Vec<Precision, 3>& origin) const
34  {
35  const vtkm::Float32 EPSILON2 = 0.0001f;
36 
37  vtkm::Vec<Precision, 3> e1 = b - a;
38  vtkm::Vec<Precision, 3> e2 = c - a;
39 
40  vtkm::Vec<Precision, 3> p;
41  p[0] = dir[1] * e2[2] - dir[2] * e2[1];
42  p[1] = dir[2] * e2[0] - dir[0] * e2[2];
43  p[2] = dir[0] * e2[1] - dir[1] * e2[0];
44  Precision dot = e1[0] * p[0] + e1[1] * p[1] + e1[2] * p[2];
45  if (dot != 0.f)
46  {
47  dot = 1.f / dot;
48  vtkm::Vec<Precision, 3> t;
49  t = origin - a;
50 
51  u = (t[0] * p[0] + t[1] * p[1] + t[2] * p[2]) * dot;
52  if (u >= (0.f - EPSILON2) && u <= (1.f + EPSILON2))
53  {
54 
55  vtkm::Vec<Precision, 3> q; // = t % e1;
56  q[0] = t[1] * e1[2] - t[2] * e1[1];
57  q[1] = t[2] * e1[0] - t[0] * e1[2];
58  q[2] = t[0] * e1[1] - t[1] * e1[0];
59  v = (dir[0] * q[0] + dir[1] * q[1] + dir[2] * q[2]) * dot;
60  if (v >= (0.f - EPSILON2) && v <= (1.f + EPSILON2) && !(u + v > 1.f))
61  {
62  distance = (e2[0] * q[0] + e2[1] * q[1] + e2[2] * q[2]) * dot;
63  }
64  }
65  }
66  }
67 
68 }; //Moller
69 
70 
71 
72 // TODO: optimization for sorting ray dims before this call.
73 // This is called multiple times and kz,kx, and ky are
74 // constant for the ray
75 
76 
77 class WaterTight
78 {
79 public:
80  template <typename Precision>
81  VTKM_EXEC inline void FindDir(const vtkm::Vec<Precision, 3>& dir,
82  vtkm::Vec<Precision, 3>& s,
83  vtkm::Vec<Int32, 3>& k) const
84  {
85  //Find max ray direction
86  k[2] = 0;
87  if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
88  {
89  if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
90  k[2] = 0;
91  else
92  k[2] = 2;
93  }
94  else
95  {
96  if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
97  k[2] = 1;
98  else
99  k[2] = 2;
100  }
101 
102  k[0] = k[2] + 1;
103  if (k[0] == 3)
104  k[0] = 0;
105  k[1] = k[0] + 1;
106  if (k[1] == 3)
107  k[1] = 0;
108 
109  if (dir[k[2]] < 0.f)
110  {
111  vtkm::Int32 temp = k[1];
112  k[1] = k[0];
113  k[0] = temp;
114  }
115 
116  s[0] = dir[k[0]] / dir[k[2]];
117  s[1] = dir[k[1]] / dir[k[2]];
118  s[2] = 1.f / dir[k[2]];
119  }
120 
121  template <typename Precision>
122  VTKM_EXEC_CONT inline void IntersectTri(const vtkm::Vec<Precision, 3>& a,
123  const vtkm::Vec<Precision, 3>& b,
124  const vtkm::Vec<Precision, 3>& c,
125  const vtkm::Vec<Precision, 3>& dir,
126  Precision& distance,
127  Precision& u,
128  Precision& v,
129  const vtkm::Vec<Precision, 3>& origin) const
130  {
131  vtkm::Vec<Int32, 3> k;
132  vtkm::Vec<Precision, 3> s;
133  //Find max ray direction
134  k[2] = 0;
135  if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
136  {
137  if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
138  k[2] = 0;
139  else
140  k[2] = 2;
141  }
142  else
143  {
144  if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
145  k[2] = 1;
146  else
147  k[2] = 2;
148  }
149 
150  k[0] = k[2] + 1;
151  if (k[0] == 3)
152  k[0] = 0;
153  k[1] = k[0] + 1;
154  if (k[1] == 3)
155  k[1] = 0;
156 
157  if (dir[k[2]] < 0.f)
158  {
159  vtkm::Int32 temp = k[1];
160  k[1] = k[0];
161  k[0] = temp;
162  }
163 
164  s[0] = dir[k[0]] / dir[k[2]];
165  s[1] = dir[k[1]] / dir[k[2]];
166  s[2] = 1.f / dir[k[2]];
167 
168  vtkm::Vec<Precision, 3> A, B, C;
169  A = a - origin;
170  B = b - origin;
171  C = c - origin;
172 
173  const Precision Ax = A[k[0]] - s[0] * A[k[2]];
174  const Precision Ay = A[k[1]] - s[1] * A[k[2]];
175  const Precision Bx = B[k[0]] - s[0] * B[k[2]];
176  const Precision By = B[k[1]] - s[1] * B[k[2]];
177  const Precision Cx = C[k[0]] - s[0] * C[k[2]];
178  const Precision Cy = C[k[1]] - s[1] * C[k[2]];
179 
180  //scaled barycentric coords
181  u = Cx * By - Cy * Bx;
182  v = Ax * Cy - Ay * Cx;
183  Precision w = Bx * Ay - By * Ax;
184  if (u == 0.f || v == 0.f || w == 0.f)
185  {
186  vtkm::Float64 CxBy = vtkm::Float64(Cx) * vtkm::Float64(By);
187  vtkm::Float64 CyBx = vtkm::Float64(Cy) * vtkm::Float64(Bx);
188  u = vtkm::Float32(CxBy - CyBx);
189 
190  vtkm::Float64 AxCy = vtkm::Float64(Ax) * vtkm::Float64(Cy);
191  vtkm::Float64 AyCx = vtkm::Float64(Ay) * vtkm::Float64(Cx);
192  v = vtkm::Float32(AxCy - AyCx);
193 
194  vtkm::Float64 BxAy = vtkm::Float64(Bx) * vtkm::Float64(Ay);
195  vtkm::Float64 ByAx = vtkm::Float64(By) * vtkm::Float64(Ax);
196  w = vtkm::Float32(BxAy - ByAx);
197  }
198  Precision low = vtkm::Min(u, vtkm::Min(v, w));
199  Precision high = vtkm::Max(u, vtkm::Max(v, w));
200 
201  Precision det = u + v + w;
202 
203  if (!((low < 0.) && (high > 0.)) && (det != 0.))
204  {
205  // Intersection is valid
206  const Precision Az = s[2] * A[k[2]];
207  const Precision Bz = s[2] * B[k[2]];
208  const Precision Cz = s[2] * C[k[2]];
209 
210  det = 1.f / det;
211 
212  u = u * det;
213  v = v * det;
214 
215  distance = (u * Az + v * Bz + w * det * Cz);
216  u = v;
217  v = w * det;
218  }
219  else
220  {
221  // Intersection is invalid
222  distance = -1.;
223  }
224  }
225 
226  template <typename Precision>
227  VTKM_EXEC inline void IntersectTriSn(const vtkm::Vec<Precision, 3>& a,
228  const vtkm::Vec<Precision, 3>& b,
229  const vtkm::Vec<Precision, 3>& c,
230  const vtkm::Vec<Precision, 3>& s,
231  const vtkm::Vec<Int32, 3>& k,
232  Precision& distance,
233  Precision& u,
234  Precision& v,
235  const vtkm::Vec<Precision, 3>& origin) const
236  {
237  vtkm::Vec<Precision, 3> A, B, C;
238  A = a - origin;
239  B = b - origin;
240  C = c - origin;
241 
242  const Precision Ax = A[k[0]] - s[0] * A[k[2]];
243  const Precision Ay = A[k[1]] - s[1] * A[k[2]];
244  const Precision Bx = B[k[0]] - s[0] * B[k[2]];
245  const Precision By = B[k[1]] - s[1] * B[k[2]];
246  const Precision Cx = C[k[0]] - s[0] * C[k[2]];
247  const Precision Cy = C[k[1]] - s[1] * C[k[2]];
248 
249  //scaled barycentric coords
250  u = Cx * By - Cy * Bx;
251  v = Ax * Cy - Ay * Cx;
252  Precision w = Bx * Ay - By * Ax;
253  if (u == 0.f || v == 0.f || w == 0.f)
254  {
255  vtkm::Float64 CxBy = vtkm::Float64(Cx) * vtkm::Float64(By);
256  vtkm::Float64 CyBx = vtkm::Float64(Cy) * vtkm::Float64(Bx);
257  u = vtkm::Float32(CxBy - CyBx);
258 
259  vtkm::Float64 AxCy = vtkm::Float64(Ax) * vtkm::Float64(Cy);
260  vtkm::Float64 AyCx = vtkm::Float64(Ay) * vtkm::Float64(Cx);
261  v = vtkm::Float32(AxCy - AyCx);
262 
263  vtkm::Float64 BxAy = vtkm::Float64(Bx) * vtkm::Float64(Ay);
264  vtkm::Float64 ByAx = vtkm::Float64(By) * vtkm::Float64(Ax);
265  w = vtkm::Float32(BxAy - ByAx);
266  }
267 
268  Precision low = vtkm::Min(u, vtkm::Min(v, w));
269  Precision high = vtkm::Max(u, vtkm::Max(v, w));
270 
271  Precision det = u + v + w;
272 
273  if (!((low < 0.) && (high > 0.)) && (det != 0.))
274  {
275  // Intersection is valid
276  const Precision Az = s[2] * A[k[2]];
277  const Precision Bz = s[2] * B[k[2]];
278  const Precision Cz = s[2] * C[k[2]];
279 
280  det = 1.f / det;
281 
282  u = u * det;
283  v = v * det;
284 
285  distance = (u * Az + v * Bz + w * det * Cz);
286  u = v;
287  v = w * det;
288  }
289  else
290  {
291  // Intersection is invalid
292  distance = -1.;
293  }
294  }
295 }; //WaterTight
296 
297 template <>
298 VTKM_EXEC inline void WaterTight::IntersectTri<vtkm::Float64>(const vtkm::Vec3f_64& a,
299  const vtkm::Vec3f_64& b,
300  const vtkm::Vec3f_64& c,
301  const vtkm::Vec3f_64& dir,
302  vtkm::Float64& distance,
303  vtkm::Float64& u,
304  vtkm::Float64& v,
305  const vtkm::Vec3f_64& origin) const
306 {
307  //Find max ray direction
308  int kz = 0;
309  if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
310  {
311  if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
312  kz = 0;
313  else
314  kz = 2;
315  }
316  else
317  {
318  if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
319  kz = 1;
320  else
321  kz = 2;
322  }
323 
324  vtkm::Int32 kx = kz + 1;
325  if (kx == 3)
326  kx = 0;
327  vtkm::Int32 ky = kx + 1;
328  if (ky == 3)
329  ky = 0;
330 
331  if (dir[kz] < 0.f)
332  {
333  vtkm::Int32 temp = ky;
334  ky = kx;
335  kx = temp;
336  }
337 
338  vtkm::Float64 Sx = dir[kx] / dir[kz];
339  vtkm::Float64 Sy = dir[ky] / dir[kz];
340  vtkm::Float64 Sz = 1. / dir[kz];
341 
342 
343 
344  vtkm::Vec3f_64 A, B, C;
345  A = a - origin;
346  B = b - origin;
347  C = c - origin;
348 
349  const vtkm::Float64 Ax = A[kx] - Sx * A[kz];
350  const vtkm::Float64 Ay = A[ky] - Sy * A[kz];
351  const vtkm::Float64 Bx = B[kx] - Sx * B[kz];
352  const vtkm::Float64 By = B[ky] - Sy * B[kz];
353  const vtkm::Float64 Cx = C[kx] - Sx * C[kz];
354  const vtkm::Float64 Cy = C[ky] - Sy * C[kz];
355 
356  //scaled barycentric coords
357  u = Cx * By - Cy * Bx;
358  v = Ax * Cy - Ay * Cx;
359 
360  vtkm::Float64 w = Bx * Ay - By * Ax;
361 
362  vtkm::Float64 low = vtkm::Min(u, vtkm::Min(v, w));
363  vtkm::Float64 high = vtkm::Max(u, vtkm::Max(v, w));
364 
365  vtkm::Float64 det = u + v + w;
366 
367  if (!((low < 0.) && (high > 0.)) && (det != 0.))
368  {
369  // Intersection is valid
370  const vtkm::Float64 Az = Sz * A[kz];
371  const vtkm::Float64 Bz = Sz * B[kz];
372  const vtkm::Float64 Cz = Sz * C[kz];
373 
374  det = 1. / det;
375 
376  u = u * det;
377  v = v * det;
378 
379  distance = (u * Az + v * Bz + w * det * Cz);
380  u = v;
381  v = w * det;
382  }
383  else
384  {
385  // Intersection is invalid
386  distance = -1.;
387  }
388 }
389 }
390 }
391 } //namespace vtkm::rendering::raytracing
392 #endif //vtk_m_rendering_raytracing_TriagnleIntersections_h
VTKM_EXEC
#define VTKM_EXEC
Definition: ExportMacros.h:51
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
VTKM_EXEC_CONT
#define VTKM_EXEC_CONT
Definition: ExportMacros.h:52
vtkm::rendering::raytracing::WaterTight::FindDir
VTKM_EXEC void FindDir(const vtkm::Vec< Precision, 3 > &dir, vtkm::Vec< Precision, 3 > &s, vtkm::Vec< Int32, 3 > &k) const
Definition: TriangleIntersections.h:81
vtkm::rendering::raytracing::Moller::IntersectTri
VTKM_EXEC void IntersectTri(const vtkm::Vec< Precision, 3 > &a, const vtkm::Vec< Precision, 3 > &b, const vtkm::Vec< Precision, 3 > &c, const vtkm::Vec< Precision, 3 > &dir, Precision &distance, Precision &u, Precision &v, const vtkm::Vec< Precision, 3 > &origin) const
Definition: TriangleIntersections.h:26
Math.h
vtkm::rendering::raytracing::Moller
Definition: TriangleIntersections.h:22
vtkm::rendering::raytracing::WaterTight::IntersectTriSn
VTKM_EXEC void IntersectTriSn(const vtkm::Vec< Precision, 3 > &a, const vtkm::Vec< Precision, 3 > &b, const vtkm::Vec< Precision, 3 > &c, const vtkm::Vec< Precision, 3 > &s, const vtkm::Vec< Int32, 3 > &k, Precision &distance, Precision &u, Precision &v, const vtkm::Vec< Precision, 3 > &origin) const
Definition: TriangleIntersections.h:227
vtkm::rendering::raytracing::WaterTight
Definition: TriangleIntersections.h:77
vtkm::Vec
A short fixed-length array.
Definition: Types.h:767
vtkm::Float32
float Float32
Definition: Types.h:154
vtkm::Int32
int32_t Int32
Definition: Types.h:160
vtkm::Float64
double Float64
Definition: Types.h:155
vtkm::rendering::raytracing::WaterTight::IntersectTri
VTKM_EXEC_CONT void IntersectTri(const vtkm::Vec< Precision, 3 > &a, const vtkm::Vec< Precision, 3 > &b, const vtkm::Vec< Precision, 3 > &c, const vtkm::Vec< Precision, 3 > &dir, Precision &distance, Precision &u, Precision &v, const vtkm::Vec< Precision, 3 > &origin) const
Definition: TriangleIntersections.h:122
Generated by   doxygen 1.8.17