How to check if a 3D point is inside a cylinder

Question

Given two 3d points and another list of 3d points, I want to check which one is inside the cylinder defined as the 3d line between the two points with a radius of r. I've implemented a numeric solution for that, which isn't accurate and too slow:

def point_in_cylinder(pt1, pt2, points, r, N=100):
    dist = np.linalg.norm(pt1 - pt2)
    ori = (pt2 - pt1) / dist 
    line = np.array([pt1 + ori*t for t in np.linspace(0, dist, N)])
    dists = np.min(cdist(line, points), 0)
    return np.where(dists <= r)[0]

I'm sure there is a better solution for that...

***** Edit *****

I speed up this function a little bit by replacing the listcomp (where the line is declared) with matrices multiplication:

line = (pt1.reshape(3, 1) + elc_ori.reshape(3, 1) @ np.linspace(0, dist, N).reshape(1, N)).T

Answer 1

To my understanding, you are generating a very large list of uniformly spaced reference points on the cylinder's axis, then checking if the minimum distance from these points to the test point is less than the radius of the cylinder.

This method is slow because each of these tests is O(N), when O(1) can be achieved (see later), but crucially it is also inaccurate because the actual test volume does not fill the entire cylinder.

The diagram below illustrates why:

^{Please forgive the bad quality.}

The white space between the cylindrical surface and the sphere leads to false negatives. To reduce this error you would need to increase N, which would in turn make the algorithm slower.

Even with (hypothetically) an infinite list of reference points covering the entire axis, the test region would still converge to a capsule, not the entire cylinder.

---

Precise O(1) method

Assuming the cylinder is defined by (p1,p2,r), check that the test point q -

1. Lies between the planes of the two circular facets of the cylinder:

   

2. Lies inside the curved surface of the cylinder:

   

---

numpy (pseudo-)code

def points_in_cylinder(pt1, pt2, r, q):
    vec = pt2 - pt1
    const = r * np.linalg.norm(vec)
    return np.where(np.dot(q - pt1, vec) >= 0 and np.dot(q - pt2, vec) <= 0 \ 
           and np.linalg.norm(np.cross(q - pt1, vec)) <= const)