In my current project I’m doing a fair amount of geometry, and one small problem I needed to solve a while back was finding whether a point is inside a cylinder.
The accepted answer for this on math.stackexchange.com wasn’t ideal — part of it was very over-complicated, and also didn’t work under some circumstances. So I contributed my own answer. In this post, in addition to the maths, I’ll give an implementation in Python.
Method
We can solve this problem by constructing the cylinder negatively:
Start with an infinite space
Throw out everything that isn't within the cylinder.
This is a classic mathematician’s approach, but it works great here, and it also works pretty well for any simply-connected solid object with only straight or concave surfaces, depending on how complex those surfaces are.
First some definitions:
Our cylinder is defined by two points, A and B, and a radius R
The point we want to test is P
The vectors from the origin to points A, B and P are \(\boldsymbol{r}_A\), \(\boldsymbol{r}_B\) and \(\boldsymbol{r}_P\) respectively.
We start with the infinite space, that is we assume all points are within the cylinder until we show they aren’t.
Then we construct 3 cuts to exclude certain values of \(\boldsymbol{r}_P\).
First, a cylindrical cut of radius R about an infinite line that goes through A and B. (This was taken from John Alexiou's answer in the link above):
-
The vector from point A to point B is:
\begin{equation*} \boldsymbol{e} = \boldsymbol{r}_B-\boldsymbol{r}_A \end{equation*}This defines the direction of the line through A and B.
-
The distance from any point P at vector \(\boldsymbol{r}_P\) to the line is:
\begin{equation*} d = \frac{\| \boldsymbol{e}\times\left(\boldsymbol{r}_{P}-\boldsymbol{r}_{A}\right) \|}{\|\boldsymbol{e}\|} \end{equation*}This is based on finding the distance of a point to a line, and using point A as an arbitrary point on the line. We could equally have used B.
-
We then simply exclude all points with \(d > R\).
This can be optimised slightly by squaring both sides of the comparison to avoid two square root operations.
Second, a planar cut that throws away the space above the "top" of the cylinder, which I'm calling A.
The plane is defined by any point on it (A will do), and any normal pointing out of the cylinder, \(-\boldsymbol{e}\) will do (i.e. in the opposite direction to \(\boldsymbol{e}\) as defined above).
-
We can see which side a point is of this plane as per Relation between a point and a plane:
The point is "above" the plane (in the direction of the normal) if:
\begin{equation*} (\boldsymbol{r}_P - \boldsymbol{r}_A) \cdot -\boldsymbol{e} > 0 \end{equation*} We exclude points which match the above.
Third, a planar cut which throws away the space below the bottom of the cylinder B.
-
This is the same as the previous step, but with the other end of the cylinder and the normal vector in the other direction, so the condition is:
\begin{equation*} (\boldsymbol{r}_P - \boldsymbol{r}_B) \cdot \boldsymbol{e} > 0 \end{equation*}
Python implementation
Below is a minimal implementation with zero dependencies outside the standard lib, in which there are just enough classes, with just enough methods, to express the algorithm neatly, following the above steps exactly.
In a real implementation:
You might separate out a
Pointclass as being semantically different fromVecYou could move some functions to be methods (in my real implementation, I have a
Cylinder.contains_point()method, for example)You should probably use
@dataclass(frozen=True)– immutable objects are a good default, I didn’t use them here because there aren’t needed and I’m focusing on clarity of the code.-
Conversely, if performance is more of a consideration, and you don’t have a more general need for classes like
VecandCylinder:you might use a more efficient representation, such as a tuple or list, or a numpy array especially if you wanted bulk operations.
you could use more generic dot product functions etc. from numpy.
you might inline more of the algorithm into a single function.
For clarity, I also have not implemented the optimisation mentioned in which you can avoid doing some square root operations.
from __future__ import annotations import math from dataclasses import dataclass # Implementation of https://math.stackexchange.com/questions/3518495/check-if-a-general-point-is-inside-a-given-cylinder # See the accompanying blog post http://lukeplant.me.uk/blog/posts/check-if-a-point-is-in-a-cylinder-geometry-and-code/ # -- Main algorithm -- def cylinder_contains_point(cylinder: Cylinder, point: Vec) -> bool: # First condition: distance from axis cylinder_direction: Vec = cylinder.end - cylinder.start point_distance_from_axis: float = abs( cross_product( cylinder_direction, (point - cylinder.start), ) ) / abs(cylinder_direction) if point_distance_from_axis > cylinder.radius: return False # Second condition: point must lie below the top plane. # Third condition: point must lie above the bottom plane # We construct planes with normals pointing out of the cylinder at both # ends, and exclude points that are outside ("above") either plane. start_plane = Plane(cylinder.start, -cylinder_direction) if point_is_above_plane(point, start_plane): return False end_plane = Plane(cylinder.end, cylinder_direction) if point_is_above_plane(point, end_plane): return False return True # -- Supporting classes and functions -- @dataclass class Vec: """ A Vector in 3 dimensions, also used to represent points in space """ x: float y: float z: float def __add__(self, other: Vec) -> Vec: return Vec(self.x + other.x, self.y + other.y, self.z + other.z) def __sub__(self, other: Vec) -> Vec: return self + (-other) def __neg__(self) -> Vec: return -1 * self def __mul__(self, scalar: float) -> Vec: return Vec(self.x * scalar, self.y * scalar, self.z * scalar) def __rmul__(self, scalar: float) -> Vec: return self * scalar def __abs__(self) -> float: return math.sqrt(self.x**2 + self.y**2 + self.z**2) @dataclass class Plane: """ A plane defined by a point on the plane, `origin`, and a `normal` vector to the plane. """ origin: Vec normal: Vec @dataclass class Cylinder: """ A closed cylinder defined by start and end points along the center line and a radius """ start: Vec end: Vec radius: float def cross_product(a: Vec, b: Vec) -> Vec: return Vec( a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x, ) def dot_product(a: Vec, b: Vec) -> float: return a.x * b.x + a.y * b.y + a.z * b.z def point_is_above_plane(point: Vec, plane: Plane) -> bool: """ Returns True if `point` is above the plane — that is on the side of the plane which is in the direction of the plane `normal`. """ # See https://math.stackexchange.com/a/2998886/78071 return dot_product((point - plane.origin), plane.normal) > 0 # -- Tests -- def test_cylinder_contains_point(): # Test cases constructed with help of Geogebra - https://www.geogebra.org/calculator/tnc3arfm cylinder = Cylinder(start=Vec(1, 0, 0), end=Vec(6.196, 3, 0), radius=0.5) # In the Z plane: assert cylinder_contains_point(cylinder, Vec(1.02, 0, 0)) assert not cylinder_contains_point(cylinder, Vec(0.98, 0, 0)) # outside bottom plane assert cylinder_contains_point(cylinder, Vec(0.8, 0.4, 0)) assert not cylinder_contains_point(cylinder, Vec(0.8, 0.5, 0)) # too far from center assert not cylinder_contains_point(cylinder, Vec(0.8, 0.3, 0)) # outside bottom plane assert cylinder_contains_point(cylinder, Vec(1.4, -0.3, 0)) assert not cylinder_contains_point(cylinder, Vec(1.4, -0.4, 0)) # too far from center assert cylinder_contains_point(cylinder, Vec(6.2, 2.8, 0)) assert not cylinder_contains_point(cylinder, Vec(6.2, 2.2, 0)) # too far from center assert not cylinder_contains_point(cylinder, Vec(6.2, 3.2, 0)) # outside top plane # Away from Z plane assert cylinder_contains_point(cylinder, Vec(1.02, 0, 0.2)) assert not cylinder_contains_point(cylinder, Vec(1.02, 0, 1)) # too far from center assert not cylinder_contains_point(cylinder, Vec(0.8, 0.3, 2)) # too far from center, and outside bottom plane