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$?

  • $\begingroup$ Are $Q$ and $R$ completely arbitrary? (for instance, $Q$ does not necessarily lie on the first sphere, or anything?) Also, for the arc $AB$, I take for granted that it's always the shorter one? What if $A$ and $B$ are exactly opposite on the first sphere (centered in $C$), that is, $A,C,B$ are aligned, in that order? Last, $A$, $B$ and $Q$ are given by $x,y,z$ coordinates? (notice that you must then ensure that the radii are equal, i.e. $AC=BC$, up to floating point accuracy) $\endgroup$ – Jean-Claude Arbaut Sep 22 '16 at 14:04

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.

  • $\begingroup$ Great explanation, but one small question: what is the value T? $\endgroup$ – 5-to-9 Sep 22 '16 at 16:25
  • $\begingroup$ See 8th equation. $\endgroup$ – John Hughes Sep 22 '16 at 16:33
  • $\begingroup$ Also...I was assuming u, v, w were unit vecs, so I need to slightly edit...but won't do so until tomorrow at the soonest. $\endgroup$ – John Hughes Sep 22 '16 at 16:35
  • $\begingroup$ I've edited, but may have introduced errors in the process. I suspect that $h$ should be $r_2 / r_1$ rather than the other way around, because that'd eliminate any trace of $r_1^4$, which seems as if it shouldn't be there. $\endgroup$ – John Hughes Sep 23 '16 at 1:03
  • $\begingroup$ "each of u and w have the same length, the radius of the first sphere divided by the radius of the second" well, at least that part makes no sense, as you now have defined their length as $ ||w|| = ||u|| = 1$. So at this point I ask myself what $h$ actually is beside an arbitrary ratio. $\endgroup$ – 5-to-9 Sep 23 '16 at 13:14

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.