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

I hope it is useful for those who found this page!