1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
"""
.. _moeller_ray_trace_example:
Visualize the Moeller-Trumbore Algorithm
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
This example demonstrates the Moeller-Trumbore intersection algorithm
using pyvista.
For additional details, please reference the following:
- `Moeller-Trumbore intersection algorithm <https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm>`_
- `Fast Minimum Storage Ray Triangle Intersectio <https://cadxfem.org/inf/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf>`_
First, define the ray triangle intersection method.
"""
from __future__ import annotations
import numpy as np
import pyvista as pv
def ray_triangle_intersection(ray_start, ray_vec, triangle):
"""Moeller-Trumbore intersection algorithm.
Parameters
----------
ray_start : np.ndarray
Length three numpy array representing start of point.
ray_vec : np.ndarray
Direction of the ray.
triangle : np.ndarray
``3 x 3`` numpy array containing the three vertices of a
triangle.
Returns
-------
bool
``True`` when there is an intersection.
tuple
Length three tuple containing the distance ``t``, and the
intersection in unit triangle ``u``, ``v`` coordinates. When
there is no intersection, these values will be:
``[np.nan, np.nan, np.nan]``
"""
# define a null intersection
null_inter = np.array([np.nan, np.nan, np.nan])
# break down triangle into the individual points
v1, v2, v3 = triangle
eps = 0.000001
# compute edges
edge1 = v2 - v1
edge2 = v3 - v1
pvec = np.cross(ray_vec, edge2)
det = edge1.dot(pvec)
if abs(det) < eps: # no intersection
return False, null_inter
inv_det = 1.0 / det
tvec = ray_start - v1
u = tvec.dot(pvec) * inv_det
if u < 0.0 or u > 1.0: # if not intersection
return False, null_inter
qvec = np.cross(tvec, edge1)
v = ray_vec.dot(qvec) * inv_det
if v < 0.0 or u + v > 1.0: # if not intersection
return False, null_inter
t = edge2.dot(qvec) * inv_det
if t < eps:
return False, null_inter
return True, np.array([t, u, v])
# %%
# Create a basic triangle within pyvista
points = np.array([[0, 0, 0], [0, 1, 0], [1, 0, 0]])
faces = np.array([3, 0, 1, 2])
tri = pv.PolyData(points, faces)
# cast a ray above pointed downwards
start = np.array([0.3, 0.25, 1])
direction = np.array([0, 0, -1])
# compute if the intersection exists
inter, tuv = ray_triangle_intersection(start, direction, points)
t, u, v = tuv
print('Intersected', inter)
print('t:', t)
print('u:', u)
print('v:', v)
# %%
# Plot the problem setup and the intersection
if inter:
# reconstruct intersection point in barycentric coordinates. See
# https://en.wikipedia.org/wiki/Barycentric_coordinate_system
a, b, c = (1 - u - v), u, v
point = tri.points[0] * a + tri.points[1] * b + tri.points[2] * c
pl = pv.Plotter()
pl.add_text(f'Intersected at ({point[0]:.3}, {point[0]:.3}, {point[0]:.3})', font_size=26)
pl.add_mesh(tri)
_ = pl.add_arrows(
np.array([start]),
np.array([direction]),
show_scalar_bar=False,
color='r',
style='wireframe',
)
pl.add_points(np.array([point]), point_size=20, render_points_as_spheres=True, color='b')
pl.add_point_labels(tri, [f'a = {1 - u - v:.3}', f'b = {u:.3}', f'c = {v:.3}'], font_size=40)
pl.show_bounds()
pl.camera_position = 'xy'
pl.show()
else: # no intersection
pl = pv.Plotter()
pl.add_text('No intersection')
_ = pl.add_arrows(
np.array([start]),
np.array([direction]),
show_scalar_bar=False,
color='r',
style='wireframe',
)
pl.add_mesh(tri)
pl.show_bounds()
pl.camera_position = 'xy'
pl.show()
|
