Arc Circle Intersection in Euclidean 3D

Question

While programming on a GPS based game I came upon the following problem for 3D euclidean space:

Assume we have an arc (given by a sphere center coordinate C and two euclidean coordinates $A$ and $B$ on that sphere, with $AC = AB$ and $A !=B$) and another sphere (given by an arbitrary center coordinate $Q$ and radius $R > 0$). All points given in cartesian coordinates $(x, y, z)$. We will exclude the edge cases where $A, B, C$ are aligned.

Now we search for intersect points $X$ (trivial: $0 <= |X| <= 2$) between the arc and the sphere.

Is there a simple formula with acceptable numerical stability (we have GPS related precision limits anyway) for efficiently computing those two (or fewer) coordinates based on the given 5 parameters $A,B,C,Q,R$?

Answer 1

I thought about it, and the second answer I proposed earlier is not all that bad. I'm going to treat the points $A$, $B$, $C$, $Q$, as coordinate triples, so that I can add and subtract them as vectors. I'm also going to write $r_2$ for $R$, the radius of the second sphere, and $r_1$ for the radius of the first. Here goes: $$ \newcommand{\a}{\mathbf a} \newcommand{\b}{\mathbf b} \newcommand{\c}{\mathbf c} \newcommand{\u}{\mathbf u} \newcommand{\v}{\mathbf v} \newcommand{\w}{\mathbf w} $$ \begin{align} r_1 &= \|A - C\| \\ r_2 & = R\\ \a &= (A - Q)/r_2 \\ \b &= (B - Q)/r_2\\ \c &= (C - Q)/r_2. \end{align}

In doing this, we've changed to a coordinate system in which the second sphere is at the origin, and has radius 1; we'll convert back at the end. Let \begin{align} \u &= \frac{r_2}{r_1}(\a - \c) \\ \v &= \frac{r_2}{r_1}(\b - \c) \\ \v' &= \v - (\u \cdot \v) \u \\ \w & = \v' / \| \v' \| \end{align}

Note that $\|u\| = \|v \| = \|w\| = 1$. And $\u \cdot \w = 0$.

Now the arc from $\u$ to $\w$ passes through $\v$ at some angle $T$. What angle? \begin{align} T &= \arccos ({\v \cdot \u}) \end{align} Let \begin{align} \gamma(t) = \c + r_1\cos(t) \u + r_1\sin(t) \w \end{align} for $0 \le t \le T$.

Then $\gamma$ parameterizes the arc from $\a$ to $\b$ on the (translated and scaled) first sphere. We seek points where \begin{align} \| \gamma(t) \|^2 &= 1. \end{align} because those are the ones that lie on the second sphere (in the new coordinate system).

Writing this out, we have \begin{align} 1 &= \| \gamma(t) \|^2 \\ &= \gamma(t) \cdot \gamma(t)\\ &= (\c + r_1\cos(t) \u + r_1\sin(t) \w) \cdot (\c + r_1\cos(t) \u + r_1\sin(t) \w) \\ &= (\c\cdot \c) + 2r_1\cos(t) (\u \cdot \c) + 2 r_1\sin(t) (\w \cdot \c) + r_1^2\cos^2(t) (\u \cdot \u) + 2r_1^2\cos(t)\sin(t) (\u \cdot \w) + r_1^2\sin^2(t) (\w \cdot \w) \end{align} Fortunately, $\u$ and $\w$ are orthogonal, and each of $\u$ and $\w$ is a unit vector. So we have \begin{align} 1 &= \| \gamma(t) \|^2 \\ &= (\c\cdot \c) + 2r_1\cos(t) (\u \cdot \c) + 2 r_1\sin(t) (\w \cdot \c) + r_1^2\cos^2(t) + r_1^2\sin^2(t) \\ &= (\c\cdot \c) + 2r_1\cos(t) (\u \cdot \c) + 2 r_1 \sin(t) (\w \cdot \c) + r_1^2 \\ \end{align}

Now we write $s = \tan(\frac{t}{2})$, which turns out to mean that $\sin t = \frac{2s}{1+s^2}$ and $\cos t = \frac{1-s^2}{1+s^2}$. So our equation becomes \begin{align} 1 &= (\c\cdot \c) + 2r_1\frac{1-s^2}{1+s^2} (\u \cdot \c) + 2 r_1\frac{2s}{1+s^2} (\w \cdot \c) + r_1^2 \\ \end{align}

Multiplying through by $1+s^2$ gives \begin{align} 1+s^2 &= (1+s^2)(\c\cdot \c) + 2r_1(1-s^2)(\u \cdot \c) + 4r_1s (\w \cdot \c) + r_1^2(1+s^2) \\ 0 &= -1 -s^2 + (1+s^2)(\c\cdot \c) + 2r_1(1-s^2)(\u \cdot \c) + 4r_1s (\w \cdot \c) + r_1^2(1+s^2) \\ 0 &= s^2 (\|c\|^2 - 1 - 2 r_1\u \cdot \c + r_1^2) + s 4r_1 (\w \cdot \c) + 2r_1(\u \cdot \c) + (\c \cdot \c) - 1 + r_1^2 \end{align} which is a quadratic $As^2 + Bs + C = 0$, where the coefficients are \begin{align} A &= (\|c\|^2 - 1 - 2 r_1u \cdot c + r_1^2)\\ B &= 4r_1(w \cdot c)\\ C &= 2r_1(u \cdot c) + (c \cdot c) - 1 + r_1^2 \end{align}

(My apologies here for reusing the names $A, B, C$.) This quadratic can be solved with the quadratic formula to get solutions $s_1, s_2$. If both solutions have an imaginary part (which happens when $B^2 - 4AC < 0$), then the arc does not intersect the sphere. Otherwise, they are real, but may be identical. Let's see how to translate them back into an overall solution to the problem.

At this point, now that we've solved the quadratic, "A, B, C" go back to their original meanings as points in 3-space, OK?

From $s_i$ ($i = 1, 2$ henceforth), we can compute $$ t_i = \arctan(2 s_i). $$ If either value of $t_i$ is outside the interval $[0, T]$, we ignore it -- it does not correspond to a point of intersection of the arc and the sphere.

From this, we can compute $\cos t_i$ and $\sin t_i$, and thus compute $$ \gamma_i = \c + r_1\cos t_i \u + r_1\sin t_i \w $$ and then, going back to the original coordinate system, compute $$ X_i = r_2\gamma_i + Q $$ and the points $X_i$ are the intersections you seek.

Answer 2

I don;t know about a formula, but this should be relatively easy:

First, let's call the arc AB by the name J; I assume that A and B are not antipodal (for then J is not unique) or equal (a trivial case), and that by "the arc from A to B" you mean "the shorted of the two pieces of the unique great circle that passes through A and B.

In that case, there are at most two points of intersection.

Look at the function on the first sphere defined by $f(X) = $ distance from $X$ to $Q$. On the arc $J$, this is either monotone or has a single extremum.

You can then use Golden Section search on $f$ to find the extremum (or to discover that it's monotone). If it's monotone, then bisection will find a point (if any) for which $f(X) = R$. If it's unimodal, Golden Section Search will find a point $D$ on the arc $AB$ such that the function is monotone on $AD$ and on $DB$, and you can apply bisection to each half to find a point for which $f(X) = R$.

Alternatively, you can parameterize the arc from $A$ to $B$ with a path $\gamma$, subtract $Q$ from everything so that $Q$ is at the origin, and divide everything by $R$ so that the radius of the second sphere is $1$. Then you're looking for values $t$ with $\| \gamma(t) \|^2 = 1$. If you use some sort of constant-speed parameterization, this should turn into a quadratic that you can solve with the quadratic formula, but it's not likely to be pretty.