# Three-dimensional ray tracing in python

I was looking for an algorithm for ray tracing in 3-D, and it took me some while to find one which can easily be programmed. Finally, I ran into this page, where a 3-D ray-tracing algorithm is clearly explained along with an example for programming. I decided to share my Python 2.7 implementation so that searchers either find this code and use it, or spend less time finding the reference to program their own code.

I wrote the algorithm in Python, slightly modified compared to the one in the link such that it is more suitable for checking whether a point is inside a (not necessarily convex) polygon.

```import numpy as np

# source: http://geomalgorithms.com/a06-_intersect-2.html

def ray_intersect_triangle(p0, p1, triangle):
# Tests if a ray starting at point p0, in the direction
# p1 - p0, will intersect with the triangle.
#
# arguments:
# p0, p1: numpy.ndarray, both with shape (3,) for x, y, z.
# triangle: numpy.ndarray, shaped (3,3), with each row
#     representing a vertex and three columns for x, y, z.
#
# returns:
#    0.0 if ray does not intersect triangle,
#    1.0 if it will intersect the triangle,
#    2.0 if starting point lies in the triangle.

v0, v1, v2 = triangle
u = v1 - v0
v = v2 - v0
normal = np.cross(u, v)

b = np.inner(normal, p1 - p0)
a = np.inner(normal, v0 - p0)

# Here is the main difference with the code in the link.
# Instead of returning if the ray is in the plane of the
# triangle, we set rI, the parameter at which the ray
# intersects the plane of the triangle, to zero so that
# we can later check if the starting point of the ray
# lies on the triangle. This is important for checking
# if a point is inside a polygon or not.

if (b == 0.0):
# ray is parallel to the plane
if a != 0.0:
# ray is outside but parallel to the plane
return 0
else:
# ray is parallel and lies in the plane
rI = 0.0
else:
rI = a / b

if rI < 0.0:
return 0

w = p0 + rI * (p1 - p0) - v0

denom = np.inner(u, v) * np.inner(u, v) - \
np.inner(u, u) * np.inner(v, v)

si = (np.inner(u, v) * np.inner(w, v) - \
np.inner(v, v) * np.inner(w, u)) / denom

if (si < 0.0) | (si > 1.0):
return 0

ti = (np.inner(u, v) * np.inner(w, u) - \
np.inner(u, u) * np.inner(w, v)) / denom

if (ti < 0.0) | (si + ti > 1.0):
return 0

if (rI == 0.0):
# point 0 lies ON the triangle. If checking for
# point inside polygon, return 2 so that the loop
# over triangles can stop, because it is on the
# polygon, thus inside.
return 2

return 1```