Cannot understand the approach of this code

Question

I wanted to find the properties of a circle from 3 points that are on its circumference, I understand how the problem should be solved, but the provided code after the explanation, does not match the approach, the code is as below:

def define_circle(p1, p2, p3):
    """
    Returns the center and radius of the circle passing the given 3 points.
    In case the 3 points form a line, returns (None, infinity).
    """
    temp = p2[0] * p2[0] + p2[1] * p2[1]
    bc = (p1[0] * p1[0] + p1[1] * p1[1] - temp) / 2
    cd = (temp - p3[0] * p3[0] - p3[1] * p3[1]) / 2
    det = (p1[0] - p2[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p2[1])
    
    if abs(det) < 1.0e-6:
        return (None, np.inf)
    
    # Center of circle
    cx = (bc*(p2[1] - p3[1]) - cd*(p1[1] - p2[1])) / det
    cy = ((p1[0] - p2[0]) * cd - (p2[0] - p3[0]) * bc) / det
    
    radius = np.sqrt((cx - p1[0])**2 + (cy - p1[1])**2)
    return ((cx, cy), radius)

also the link for the explanation is: web archive

Can someone explain to me how this code is working?

Answer 1

Let $u,v,w$ be the three points. Let $u=(a,\alpha)$, $v=(b,\beta)$, and let $w=(c,\gamma)$. Let $(d,\delta)$ be the center of the circle. Consider the matrix equation $$2\begin{pmatrix} a-b & \alpha - \beta \\ b-c & \beta-\gamma\end{pmatrix} \begin{pmatrix}d \\ \delta\end{pmatrix}.$$ The first component of the resulting vector is \begin{align*} 2(a-b)d + 2(\alpha-\beta)\delta & = 2ad +2\alpha \delta - 2bd - 2\beta \delta \\ & = (a^2+2ad+d^2)+(\alpha^2+2\alpha \delta+\delta^2) -a^2-\alpha^2 \\ & - (b^2+2bd+d^2) -(\beta^2 +2\beta\delta + \delta^2) + b^2+\beta^2 \\ & = [(a-d)^2+(\alpha-\delta)^2] - (a^2+\alpha^2) \\ & - [(b-d)^2+(\beta-\delta)^2]+b^2+\beta^2 \\ & = b^2+\beta^2 - (a^2+\alpha^2).\end{align*} The terms in the square boxes $[\ldots]$ cancel because they're each equal to the squared radius of the circle. The first is the squared distance from $(a,\alpha)$ to the center (which is $r^2$), and the second is the squared distance from $(b,\beta)$ to the center. This is $2\times$ the $bc$ term from the code.

A similar calculation yields that the second entry of the matrix product is $2\times$ the $cd$ term. In other words, $$2\begin{pmatrix} a-b & \alpha - \beta \\ b-c & \beta-\gamma\end{pmatrix} \begin{pmatrix}d \\ \delta\end{pmatrix}=\begin{pmatrix} \|(b,\beta)\|^2 - \|(a,\alpha)\|^2 \\ \|(b,\beta)\|^2 - \|(c,\gamma)\|^2\end{pmatrix}.$$ So we multiple both sides of this equation by the inverse matrix of $2\begin{pmatrix} a-b & \alpha - \beta \\ b-c & \beta-\gamma\end{pmatrix}$, which is what the code is doing to get the center.

If the matrix is not invertible, it means the three points lie on a line, because the determinant is a multiple of the signed area of the parallelogram formed by $(a,\alpha)-(b,\beta)$ and $(b,\beta)-(c,\gamma)$, which is why it returns no center and infinite radius when $\det$=0.

After we know a point on the circle (actually, three) and the center, the radius is straightforward.

EDIT: For where the first multiplication came from, we want to know the center of the circle, $(d,\delta)$. We write down $$(a-d)^2+(\alpha-\delta)^2=r^2,\tag{$1$}$$ $$(b-d)^2+(\beta-\delta)^2=r^2,\tag{$2$}$$ $$(c-d)^2+(\gamma-\delta)^2=r^2.\tag{$3$}$$ Subtract the second equation from the first and all $d^2,\delta^2$, and $r^2$ terms cancel. You're left with an equation involving $a,b,\alpha,\beta$ (all known) and $d,\delta$ (unknown). Subtract the third equation from the second and something similar happens. The two resulting equations when I subtract $(2)$ from $(1)$ and $(3)$ from $(2)$ are exactly what the first matrix equation is saying (although there's a little rewriting to get it into the same form).