Distance from circle to line segment: the intuition or geometry behind this algorithm

algorithmscollision detectiongeometryvectors

I am using this algorithm to compute the distance between a circle and a line segment. I do not understand the part of the algorithm that involves dot products (and that kind of bothers me, I want to understand what I'm doing). To quote it in case the site ever goes down:

  // get length of the line
  float distX = x1 - x2;
  float distY = y1 - y2;
  float len = sqrt( (distX*distX) + (distY*distY) );

  // get dot product of the line and circle
  float dot = ( ((cx-x1)*(x2-x1)) + ((cy-y1)*(y2-y1)) ) / pow(len,2);

  // find the closest point on the line
  float closestX = x1 + (dot * (x2-x1));
  float closestY = y1 + (dot * (y2-y1));

  // is this point actually on the line segment?
  // if so keep going, but if not, return false
  boolean onSegment = linePoint(x1,y1,x2,y2, closestX,closestY);
  if (!onSegment) return false;

Where x1, y1, x2, y2 are the coordinates defining the line segment, and cx, cy and r describe the circle's coordinates and radius.

Thank you. I'm not great at math (but not terrible) so emphasis on the intuition.

(side note: if there's an easier way to do this, let me know!)

Best Answer

There are two things needed here (after it has already been tested that the endpoints are outside the circle):

  • The distance $d$ from the center $C$ of the circle to the line, and
  • The location of the point $Q$ on the line which is closest to $C$.

If $d > R$, then the line does not intersect the circle and we are done (if the line has no intersection, then neither does the line segment).

If $d \le R$, then the line does intersect the circle. But that doesn't mean the line segment intersects. It could be on a portion of the line away from the circle. For this we need $Q$, the point on the line closest to $C$. If $Q$ is between $P_1$ and $P_2$, then the segment intersects the circle.

If $P_1$ is between $Q$ and $P_2$, or vice versa, then the segment does not intersect the circle. We've already concluded that $P_1$ and $P_2$ are outside the circle. If there were a point of intersection between the segment and the circle, since $Q$ is inside the circle, everything between the intersection and $Q$ would also be inside the circle, including one of the two endpoints that we know are outside - a contradiction. So no such point of intersection can exist.

So the algorithm finds the coordinates of $Q$, calculates the distance $d$ to $C$, and checks to see if $Q$ is on the line segment $\overline{P_1P_2}$, all in a way that woefully inefficient. I've already pointed out the most glaring of the inefficiencies in a comment.

The key to finding $Q$ is that the line $\overline{CQ}$ is perpendicular to $\overline{P_1P_2}$. That is, it is the projection of $C$ onto the line. For reasons that 3Blue1Brown can explain much better than I can, the dot product of the vector $\vec{P_1C}$ and the vector $\vec{P_1P_2}$ is (the length of the projected vector $\vec{P_1Q}$) multiplied by (the length $L$ of $\vec{P_1P_2}$).

The value dot in your algorithm is not the dot product itself, but the dot product divided by $L^2$ ($L$ is len in the algorithm). Dividing by $L$ once gives the length of $\vec{P_1Q}$, that is, the distance from $P_1$ to $Q$. Dividing by $L$ a second time gives the ratio of the distance $P_1Q$ to the distance $P_1P_2$: $$\mathtt{dot} = \dfrac{P_1Q}{P_1P_2}$$ Once they know this ratio, they find $Q$ by the vector equation $$\vec Q = \vec P_1 + (\mathtt{dot})(\vec P_2 - \vec P_1)$$ which is the equation for the line through $P_1$ and $P_2$.

That explains their algorithm.


Now if you want a good algorithm, there are better ways to do it.

The direction of $\vec{CQ}$ is perpendicular to $\vec{P_1P_2}$. Because we are in a plane, there is only one direction (up to a change in sign) perpendicular to any nonzero vector. If $\vec{P_1P_2} = \begin{bmatrix} p_x\\p_y\end{bmatrix}$, then we are looking for a vector $\begin{bmatrix} n_x\\n_y\end{bmatrix}$ whose dot product with $\begin{bmatrix} p_x\\p_y\end{bmatrix}$ is $0$: $$n_xp_x + n_yp_y = 0.$$ An easy trick for accomplishing this is just to set $n_x = p_y, n_y = -p_x$. This gives a vector $\vec n$ either parallel to $\vec{CQ}$ or anti-parallel (i.e., parallel to $\vec{QC}$). The length of $\vec n$ is the same as $\vec{P_1P_2}$, which is $L$.

Now instead of taking the dot product of $\vec{P_1C}$ with $\vec{P_1P_2}$, we take the dot product of $\vec{CP_1}$ with $\vec n$. Again per the interpretation of the dot product discussed in the 3Blue1Brown video, that dot product is going to be the length of the projection of $\vec{CP_1}$ onto the line of $\vec n$ times the length of $\vec n$, or its opposite. But the line of $\vec n$ is $\overline{CQ}$, so the projection is going to be $\vec{CQ}$ (because $\angle CQP_1$ is a right angle), and its length is just $d$, while the length of $\vec n$ is $L$: $$\vec{CP_1} \cdot \vec n = \pm dL$$ so $$d^2 = \frac{(\vec{CP_1} \cdot \vec n)^2}{L^2}$$ The full line intersects the circle if and only if $d\le R$, or equivalently $d^2 \le R^2$ (no need to take square roots in this version).

Thus you do not need to know $Q$ in order to find $d$. However, you still need to know where $Q$ falls on the line through $P_1, P_2$. But all you need for that is dot. Another major inefficiency in their algorithm was calling a "linePoint" function which assuredly was written to test arbitrary points to see if they are on the line through $P_1, P_2$, to test if $Q$ is between $P_1$ and $P_2$. We already know that $Q$ is on the line, so most of the testing in that function is not necessary. And very likely, it tests for lying between $P_1$ and $P_2$ by recalculating the value of dot. For points $Q$ on the line, $0 \le \mathtt{dot} \le 1$ when $Q$ is between $P_1$ and $P_2$. (It is $< 0$ when $Q$ is on the other side of $P_1$ and $>1$ when $Q$ is on the other side of $P_2$, this follows from dot being the ratio of $P_1Q$ to $L$, where $P_1Q$ is negative when $Q$ is outside of $P_1$.)

This gives a nicer algorithm:

boolean IsCollisionCircleLine(float cx, float cy, float r, float x1, float y1, float x2, float y2)
    // Change coordinate origin to be at the circle center
    x1 = x1 - cx
    y1 = y1 - cy
    x2 = x2 - cx
    y2 = y2 - cy

    // Square the radius to avoid needing any square roots
    float r2 = r * r
    
    // Check if endpoints are in the circle
    if (x1 * x1 + y1 * y1 <= r2) return true
    if (x2 * x2 + y2 * y2 <= r2) return true
   
    // Calculate the line segment length squared
    float len2 = (x1 - x2) * (x1 - x2) + (y1 - y2)*(y1 - y2)

    // find the perpendicular vector to the line segment
    float nx = y2 - y1
    float ny = x1 - x2

    // find the distance squared from center to line, times len2
    float dist2 = nx * x1 + ny * y1
    dist2 = dist2 * dist2

    // Check if full line intersects the circle. If not, return false
    if (dist2 > len2 * r2) return false

    // Calculate the distance from (x1,y1) to the point of closest approach) times the segment length
    float index = (x1*(x1 - x2) + y1*(y1 - y2)) 

    // Check if point of closest approach is inside the segment.
    if (index < 0) return false
    if (index > len2) return false
    return true
Related Question