Intersection of ray with circle in 3D subject to endpoint and orientation constrains

Imagine a circle of radius R in 3D space with a line l running threw it's center C in a direction perpendicular to the plane of the circle. Basically, like the axel of a wheel.

From a given point P that is not on the circle or on l, a ray extends to intersect both l and the circle. What would be the equations used to find the intersection points the ray make with the circle and l? You are given the coordinates of C and P, the radius R and the orientation of l.

I am trying to model looking from a point P onto a wheel-axis shape and find from point of view P the point of the edge of the circle that would appear to intersect with it's axis. Of course it doesn't but it is how this 3d structure would appear in a 2d image if a camera was situated at point P.

You should have simply stated your question in terms of orthographic projection from the get go. Try the following algorithm, I hope it works (if, of course, I understand correctly what your goal is).

This time you assume you are given the vector $\vec{p}$ instead of the point $P$ that determines the direction of all rays illuminating the circle. Last time, all rays were emanating from the point $P$. This time, all rays are parallel to $\vec{p}$. Everything else, including the notations, stays the same. The theoretical arguments are almost word for word the same, except this time the ray $s$ does not pass through point $P$ but is parallel to $\vec{p}$ instead.

The algorithm is a bit simpler this time, depending what points you really need to find. I will find all of them: $Q_1, Q_2, R_1, R_2$ but you can choose which ones you really need.

1. Calculate $$\cos(\theta) = \frac{\big( \vec{v} \cdot \vec{p}\big)}{|\vec{v}||\vec{p}|}$$ where $\theta = \angle\,(\vec{v},\vec{p})$ is angle between first vector $\vec{v}$ and second vector $\vec{p}$, following that order.

2. If $\cos(\theta) < 0$ then set $$\vec{v} := -\vec{v}$$ $$\cos(\theta) := - \cos(\theta)$$ Else, keep $\vec{v}$ and $\cos(\theta)$ the same. This way $\vec{v}$ and $\vec{p}$ are in the same half-space with respect to the plane $\beta$ determined by the circle (and so $\beta$ is orthogonal to $\vec{v}$).

Observe that since $Q_1R_1$ is parallel to $\vec{p}$ $$\angle \, Q_1R_1C = \angle (\vec{v}, \vec{p}) = \theta$$

1. Calculate $$\vec{u} = \frac{ \vec{v} \times \big(\vec{v}\times\vec{p} \big)}{|\vec{v} \times \big(\vec{v}\times\vec{p}\big)|}$$ This vector is calculated so that it points from $C$ to $Q_1$ which is the point behind $l$.

2. Calculate $$\vec{CQ_1} = r \vec{u} \,\,\, \text{ and } \,\,\, \vec{CQ_2} = - r \vec{u}$$ which leads to $$\vec{OQ_1} =\vec{OC} + r \vec{u} \,\,\, \text{ and } \,\,\, \vec{OQ_2} = \vec{OC} - r \vec{u}$$

If you look at triangle $CQ_1R_1$, you will see that it is right-angled triangle with angle $\angle \, Q_1R_1C = \theta$ and $|CQ_1| = r$ so $$|CQ_1| = r \, \cot(\theta) = \frac{r \,\cos(\theta)}{\sqrt{1 - \cos^2(\theta)}}$$ Thus

1. Calculate $$\vec{CR_1} = \frac{r \,\cos(\theta)}{\sqrt{1 - \cos^2(\theta)}} \, \frac{\vec{v}}{|\vec{v}|}\,\,\,\,\, \text{ and } \,\,\,\,\, \vec{CR_2} = - \, \frac{r \,\cos(\theta)}{\sqrt{1 - \cos^2(\theta)}} \, \frac{\vec{v}}{|\vec{v}|}$$ which lead to $$\vec{OR_1} = \vec{OC} + \frac{r \,\cos(\theta)}{\sqrt{1 - \cos^2(\theta)}} \, \frac{\vec{v}}{|\vec{v}|}\,\,\,\,\, \text{ and } \,\,\,\,\, \vec{OR_2} = \vec{OC} - \, \frac{r \,\cos(\theta)}{\sqrt{1 - \cos^2(\theta)}} \, \frac{\vec{v}}{|\vec{v}|}$$
• Thank you for all your efforts. Really nice work! The CAD system I use can switch from orthographic to perspective projections. So both of your articles will help and you did not waste time on the perspective approach at all. – Mark Apr 30 '17 at 7:35
• The apparent intersections are mapping onto k and l as they should. They overlap perfectly in the view direction and when I change view, they lay on the line and circle as they should. But I must have done something wrong because in step 3, the vector from C to Q1 should point "behind" l. Mine points "in front" of l, closer to the camera. – Mark Apr 30 '17 at 18:09
• @Mark Oh wait, maybe there is nothing wrong, maybe you and I interpret the vector $\vec{p}$ differently, that's why. Does the vector $\vec{p}$ you use point from you to the point $C$ or from point $C$ to you? I have assumed it points from $C$ to you, which may be the source of the discrepancy. You may first invert it by $\vec{p} : = - \vec{p}$ to match my interpretation, or in my interpretation multiply the vector $\vec{p}$ by $-1$ everywhere. – Futurologist Apr 30 '17 at 19:12
• Yup that is it. The vector points from me to point C (i.e. out from the camera). Great, thanks, have a good weekend ... I am going dancing. LOL – Mark Apr 30 '17 at 19:54

No need to rotate anything, it is a matter of very simple geometry, which if you follow through, gives you a very simple explicit algorithm for computing the points you need.

Assume the orientation of the lines $l$ is given by a vector $\vec{v}$. Then the circle, call it $k$, has given center $C$ and radius $r$ and lies in the plane $\beta$ passing through $C$ and orthogonal to $\vec{v}$. Denote by $s$ the ray through point $P$ that intersects both line $l$ and the circle $s$. Then, since $s$ intersects $l$, the two together determines a plane $\alpha$, which is orthogonal to the plane of the circle $\beta$ and is transverse to the circle itself. The ray $s$ lies in this plane $\alpha$ and at the same time intersects $k$ so the ray $s$ passes through the point of intersection of $k$ and $\alpha$. So all you have to do is find the intersection points of the plane $\alpha$ and the circle $k$ (technically, you have two solutions of your problem). It becomes even simpler when you notice that the planes $\alpha$ and $\beta$ (the one of the circle $k$) intersect at a common line $l_C$ that passes thorough the circle center $C$, that is $l_C = \alpha \cap \beta$. Therefore the intersection points between $k$ and $\alpha$ are in fact the intersection points of $l_C$ and $k$. In other words, the two points you are looking for are exactly $l_C \cap k$. Which, by the way, are the two points on $l_C$ at a distance $r$ from point $C \in l_C$ (on either side of $C$ on $l_C$).

All of the above observations prompt the following algorithm:

1. If the dot product $\big(\vec{v} \cdot \vec{CP}\big) < 0$ then set $\vec{v} := -\vec{v}$. This way we make sure both vectors $\vec{v}$ and $\vec{CP}$ are in the same half space with respect to the plane $\beta$ defined by the circle (recall $\beta$ is orthogonal to $\vec{v}$).

2. Define vector $$\vec{n} = \vec{v} \times \vec{CP}$$ (cross product) which is orthogonal to $\alpha$, and thus orthogonal to $l_C$.

3. Then define vector $$\vec{w} = \vec{v} \times \vec{n} = \vec{v} \times\big(\vec{v} \times \vec{CP}\big)$$ and normalize it to $$\vec{u} = \frac{\vec{w}}{|\vec{w}|} = \frac{ \vec{v} \times\big(\vec{v} \times \vec{CP}\big)}{| \vec{v} \times\big(\vec{v} \times \vec{CP}\big)|}$$ Vector $\vec{u}$ is parallel to line $l_C$ because $l_C$ is orthogonal to both vectors $\vec{v}$ and $\vec{n}$, and vector $\vec{u}$ is also orthogonal to both of them (cross product of the two).

4. If point $O$ is the origin of the coordinate system, a point $X$ lies on the line $l_C$ if and only if $$\vec{OX}= \vec{OC} + t \, \vec{u}$$

5. The two intersection points $Q_1$ and $Q_2$ you are looking for are $$\vec{OQ_1}= \vec{OC} + r \, \vec{u}$$ $$\vec{OQ_2}= \vec{OC} - r \, \vec{u}$$

If I am not wrong, according to the way I have deliberately defined the relative location of the vectors, point $Q_1$ should be "behind" the line $l$ and point $Q_2$ "in front" when looking from point $P$. For further reference I will use the notations $R_1 = PQ_1 \cap l$ and $R_2 = PQ_2 \cap l$.

1. For $i=1,2$ calculate the vectors $$\vec{Q_iP} = \vec{CP} - (-1)^{i} \, r \, \vec{u} \,\,\, \text{ and } \,\,\, |Q_iP| = \sqrt{\big(\vec{Q_iP} \cdot \vec{Q_iP} \big)}$$ i.e. the latter is a dot product and then square root. The former equality holds because $$\vec{Q_iP} = \vec{OP} - \vec{OQ_i} = \vec{CP} - \vec{CQ_i} = \vec{CP} - (-1)^i \, r \, \vec{u}$$

2. For $i=1,2$ calculate $\cos(\alpha_i) = \cos\big( \angle \, CQ_iP\big)$ by calculating the dot products $$\cos(\alpha_i) = (-1)^i\, \frac{\big(\, \vec{u} \cdot\vec{Q_iP} \,\big)}{|{Q_iP}|}$$

3. For $i=1,2$ calculate $$\vec{OR_i} = \vec{OQ_i} + \left(\,\,\frac{r}{|Q_iP| \,\cos(\alpha_i)}\,\right) \, \vec{Q_iP}$$ where by construction $R_1$ is in between $Q_1$ and $P$ (i.e. $Q_1$ is behind the line $l$ when looking from point $P$) while $Q_2$ is outside the straight segment formed by $Q_2$ and $P$.

If you manage to write a computer implementation of this algorithm, let me know if it works or not. If not, I will look up what needs to be corrected.

• Thanks for help on finding the points on the edge of the circle that correspond to the "apparent" intersections. I tried it on my 3d CAD program and it worked great. Now i would like to find where these apparent points correspond to points on the axis. Also, find whether the portion of the circle around the apparent intersections is "in front" or "behind" the axis. Thanks again. – Mark Apr 28 '17 at 0:21
• @Mark I have added extra steps in the algorithm. Read my post again and let me know if it works. – Futurologist Apr 29 '17 at 4:12
• Thank you again! I will look it over this weekend and let you know next week. One thing is to distinguish between orthographic projections and perspective ones. I am working with orthographic ones but it is a minor issues to use parallel lines to the camera view direction. – Mark Apr 29 '17 at 5:38
• @Mark I added a post with an algorithm for the case of orthographic projection, i.e. with parallel rays instead of rays emanating from a common point $P$. Check it out and let me know if it helps. – Futurologist Apr 30 '17 at 3:49

Rotate everything so that the disk lies in the $xy$-plane with its center at the origin, and the line $L$ coincides with the $z$-axis. One more rotation about the $z$-azis places $p$ to project onto the $x$-axis. Now it is easy to connect $p$ through $L$ to the circle. Then reverse all rotations.

• If i have the points on the circle that go with the apparent intersections then i need the associated points on the axis. Also, i need to know whether the circle at an apparent intersection is closer or further than the axis from the camera. I thought a line parallel to the view direction but through an apparent intersection would intersect the axis. Distances from the camera of the point on the circle and axis could then be compared. However, there maybe a more clever/elegant solution. – Mark Apr 27 '17 at 23:19
• @Mark: "i need the associated points on the axis." Easily achieved. Just similar triangles. But I cannot follow the remainder of your remarks. I think we are viewing the situation rather differently. – Joseph O'Rourke Apr 27 '17 at 23:28
• Maybe it is better to describe what i want to achieve. Basically. I want to illustrate either the cicle "in front" of the axis at an apparent intersection or "behind" it by removing some of the axis or some of the circle at these points. I am using a 3d cad program so when i change perspectives, i can have a script that redraws the axis and circle, cutting out small portions to leave white space around these apparent intersections. Hope i made sense. LOL – Mark Apr 27 '17 at 23:44
• Sorry, @Mark, I don't think I can help you. Your needs are rather specific to your goals. – Joseph O'Rourke Apr 27 '17 at 23:53