Given three distinct points on a sphere, find the unique round circle they live in

Question

Say you have three (distinct) points on the unit sphere in Euclidean space

$$p_1, p_2, p_3 \in S^n = \{ x \in \mathbb R^{n+1} : |x| = 1 \}$$

I'd like to find, as efficiently and robustly as possible, a description of the unique round circle in $S^n$ that contains the three points. By round circle I mean the intersection of an affine $2$-dimensional subspace of $\mathbb R^{n+1}$ with $S^n$.

I'm mostly interested in the $n=3$ case, although the $n>3$ cases are of some interest to me as well.

Off the top of my head the most sensible way to accomplish this would be:

1. Find the smallest-norm convex-combination of the $p_i$'s, i.e. solve $$\min \left\|\sum_i \alpha_i p_i\right\|, \hskip 1cm \sum_i \alpha_i = 1$$ which is a calculus problem.

2. Replace the $p_i$'s by $p_i - q$ where $q$ is the above norm-minimizer.

3. Compute the orthogonal complement of $span(p_1, p_2, p_3)$.

If $q_1, q_2, \cdots, q_{n-2}$ is a basis for the orthogonal complement then that would give a system of equations describing the circle.

$$q_i \cdot (x-p_j) = 0$$

for all $i,j$ (these equations would technically be independent of j).

The nice thing about this setup is it's just linear algebra. One problem with this solution is step (2) -- if the norm-minimizer results in a very small but non-zero norm there might be numerical instabilities. In the application I have in mind, there will be (potentially) billions of such computations and these kinds of instabilities will be difficult to avoid.

The $n=2$ case has a rather cute, stable and efficient solution which (off the top of my head) I don't see how to generalize

$$Det \pmatrix{ x & y & z & 1 \cr \cdot & \cdot & \cdot & 1 \cr \cdot & \cdot & \cdot & 1 \cr \cdot & \cdot & \cdot & 1} = 0$$

where you plug the three points $p_1, p_2, p_3$ into the dotted rows.

If there was a more general solution of this kind, that would be wonderful as it solves the stability issue. On top of that, it's a simple closed-form solution and my application could use a solution that is easily differentiated. I don't need that, but it would be useful.

edit: I need an answer that gives the "universal bundle" description of the circle, i.e. an equation of the plane the circle lives in. You could think of this as a point in the Grassmannian $G_{n+1,2}$ together with a vector in the orthogonal complement of the 2-dimensional subspace. This is because I need ready access to the Hausdorff distance function (minimum distance) from points in the sphere to the circle. i.e. a parametrization of the circle is not enough.

edit 2: in my comment below I refer to the matrix formulation of the Grassmannian. In this formulation, the space of 2-dimensional subspaces of $\mathbb R^n$ is the space of $n\times n$ matrices $A$ such that: $A^t = A, A^2=A, \text{ and } tr(A) = 2$. From this perspective the matrix $A$ represents the orthogonal projection map onto its image, which is a $2$-dimensional subspace.

Answer 1

Let $(P)$ be the plane through $p_1,p_2,p_3$.

$(P)$ can be described as the set of points

$$p_{u,v}:=p_1+u(p_2-p_1)+v(p_3-p_1), \ \ \ \ \ \text{for} \ u,v \in \mathbb{R}.$$

Let $\Omega$ be the orthogonal projection of the origin $O$ onto $(P)$ (obtained for example by minimizing the norm of $p_{u,v}$).

$\Omega$ is naturally the center of the circumscribed circle to triangle $p_1,p_2,p_3$.

It suffices then to orthonormalize the (vector) basis $\{(p_2-p_1),(p_3-p_1)\}$ into the (vector) basis $\{q_1,q_2\}$ to be able to describe the circle in the following way:

$$p=\Omega+r \cos(t) q_1 + r \sin(t)q_2, \ \ \ t \in [0,2 \pi]$$

with $r>0$ defined by equality $\vec{O\Omega}^2+r^2=1.$

Remark: this works in any dimension.

Answer 2

A student helped me realize there are a few terrific answers to my question. The flashiest, although not exactly what I want, is via the Cayley-Menger formula. This is the formula that gives the $n$-dimensional content of the convex hull of $(n+1)$-points in Euclidean space, as a determinant.

If you take the convex hull of $(n+1)$ points in $\mathbb R^m$, say the points $p_0, p_1, \cdots, p_n$ and let $d_{ij}$ be the distance between the $i$-th and $j$-th point, $d_{ij} = |p_i-p_j|$ then the $n$-dimensional content of the convex hull of these points is given by:

$$\mu_n (Hull(p_0,p_1,\cdots,p_n)) = \sqrt{\frac{(-1)^{n+1}}{(n!)^22^n}\begin{vmatrix} 0 & d_{01}^2 & d_{02}^2 & \cdots & d_{0n}^2 & 1 \\ d_{01}^2 & 0 & d_{12}^2 & \cdots & d_{1n}^2 & 1 \\ d_{02}^2 & d_{12}^2 & 0 & \cdots & d_{2n}^2 & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ d_{0n}^2 & d_{1n}^2 & d_{2n}^2 & \cdots & 0 & 1 \\ 1 & 1 & 1 & \cdots & 1 & 0 \end{vmatrix}\ }.$$

To relate this to my question, given three distinct points $p_1,p_2,p_3 \in S^3$ then the equation for $q \in S^3$ to be on the unique round circle containing the points $p_1,p_2,p_3$ is:

$$\mu_3(Hull(q,p_1,p_2,p_3)) = 0.$$

The problem with this is four points on a round circle is co-dimension two in the space of all 4-tuples of points in $S^3$, so $0$ can't be a regular value. I'd really like the manifold to be the pre-image of a regular value. But the normal bundle gives the clue, as 4-tuples that are on a round circle have a normal bundle induced from $TS^2$, the tangent bundle.

There is a map

$$f : C_4(S^3) \to S^2 \times S^2$$

where $C_4(S^3)$ is the space of distinct $4$-tuples of points in $S^3$, such that the diagonal $\Delta S^2 \subset S^2 \times S^2$ is a regular value, and $f^{-1}(\Delta S^2)$ is the subspace of $4$-tuples of points on a round circle. $\Delta S^2 = \{(p,p) : p \in S^2 \}$.

The way you define the map is you use the quaternionic structure of $S^3$ to write $C_4(S^3) \simeq S^3 \times C_3(\mathbb R^3)$. Then the 4-tuples of points on a round circle corresponds to the triples of points on a line in $\mathbb R^3$ (cross $S^3$). So the function is:

$$g : C_3(\mathbb R^3) \to S^2 \times S^2$$

given by $$g(p_1,p_2,p_3) = \left(\frac{p_2-p_1}{|p_2-p_1|}, \frac{p_3-p_2}{|p_3-p_2|}\right)$$

There's two further variations of this function for the three different cyclic orderings of the four points. This is the kind of function I can use for my problem -- I need to do a Newton-method approximation.