Ray Sphere Intersection Brute Force

In this notebook is the general solution to a ray sphere interesection worked out in what might be considered the obvious manner. It is worth understanding as an illustration of the general approach. It is NOT the best way to compute ray sphere intersections as will be highlighted in the notebook to follow.

Ross Beveridge, October 9, 2018

latex.matrix_delimiters("[", "]")
import sympy
var('r, s');
Cv = vector(SR, 3, var('cx, cy, cz'));
Lv = vector(SR, 3, var('lx, ly, lz'));
Uv = vector(SR, 3, var('ux, uy, uz'));
Pv = vector(SR, 3, var('x, y, z'));
Tv = vector(SR, 3, var('tx, ty, tz'));
pretty_print("Sphere center: C = ", Cv);
pretty_print("Ray start:     L = ", Lv);
pretty_print("Ray direction: U = ", Uv);
pretty_print("Base to Center: T = ", Cv - Lv);

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Sphere|\phantom{\verb!x!}\verb|center:|\phantom{\verb!x!}\verb|C|\phantom{\verb!x!}\verb|=| \left(\mathit{cx},\,\mathit{cy},\,\mathit{cz}\right)$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Ray|\phantom{\verb!x!}\verb|start:|\phantom{\verb!xxxxx!}\verb|L|\phantom{\verb!x!}\verb|=| \left(\mathit{lx},\,\mathit{ly},\,\mathit{lz}\right)$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Ray|\phantom{\verb!x!}\verb|direction:|\phantom{\verb!x!}\verb|U|\phantom{\verb!x!}\verb|=| \left(\mathit{ux},\,\mathit{uy},\,\mathit{uz}\right)$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Base|\phantom{\verb!x!}\verb|to|\phantom{\verb!x!}\verb|Center:|\phantom{\verb!x!}\verb|T|\phantom{\verb!x!}\verb|=| \left(\mathit{cx} - \mathit{lx},\,\mathit{cy} - \mathit{ly},\,\mathit{cz} - \mathit{lz}\right)$$

The implicit form that defines all points on the sphere of radius r is shown below as is the parametric form for the ray. A note about presentation, the ray is being shown as a vector that is written out as a row vector. This makes no substantive difference, but the pretty print version of the ray may not at first match your expectation.

eq1 = ((Pv - Cv)*(Pv - Cv) - r^2 == 0)
ray = Lv + s * Uv; 
pretty_print(eq1)
pretty_print(LatexExpr("{R(s)\:=\:}"),ray)

$$\newcommand{\Bold}[1]{\mathbf{#1}}{\left(\mathit{cx} - x\right)}^{2} + {\left(\mathit{cy} - y\right)}^{2} + {\left(\mathit{cz} - z\right)}^{2} - r^{2} = 0$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}{R(s)\:=\:} \left(s \mathit{ux} + \mathit{lx},\,s \mathit{uy} + \mathit{ly},\,s \mathit{uz} + \mathit{lz}\right)$$

There is an obvious way to take these two equations and create a new single equation in a single unknown, s, the ray parameter. Namely, substitute the X, Y, and Z components from the ray into the sphere equation.

eq2 = (eq1.lhs()(x=ray[0],y=ray[1],z=ray[2]) == 0)
pretty_print(eq2)

$$\newcommand{\Bold}[1]{\mathbf{#1}}{\left(s \mathit{ux} - \mathit{cx} + \mathit{lx}\right)}^{2} + {\left(s \mathit{uy} - \mathit{cy} + \mathit{ly}\right)}^{2} + {\left(s \mathit{uz} - \mathit{cz} + \mathit{lz}\right)}^{2} - r^{2} = 0$$

Now use the T vector introduced above that goes from the sphere center to the base of the ray.

eq2s = sympy.sympify(eq2.lhs());
eq2ss = SR(eq2s.subs([(cx-lx,tx),(cy-ly,ty),(cz-lz,tz)]));
eq3 = (eq2ss == 0); 
pretty_print(eq3)

$$\newcommand{\Bold}[1]{\mathbf{#1}}{\left(s \mathit{ux} - \mathit{tx}\right)}^{2} + {\left(s \mathit{uy} - \mathit{ty}\right)}^{2} + {\left(s \mathit{uz} - \mathit{tz}\right)}^{2} - r^{2} = 0$$

Expand the terms and collect relative to the variable s.

eq4  = eq3.expand(); 
eq4l = eq4.lhs().collect(s);
eq4  = (eq4l == 0); 
pretty_print(eq4);

$$\newcommand{\Bold}[1]{\mathbf{#1}}{\left(\mathit{ux}^{2} + \mathit{uy}^{2} + \mathit{uz}^{2}\right)} s^{2} - r^{2} - 2 \, {\left(\mathit{tx} \mathit{ux} + \mathit{ty} \mathit{uy} + \mathit{tz} \mathit{uz}\right)} s + \mathit{tx}^{2} + \mathit{ty}^{2} + \mathit{tz}^{2} = 0$$

Now a key constraint enters into play, the direction vector U is unit length so the s squared term's lead constant equals 1

eq4s  = sympy.sympify(eq4.lhs());
eq4l = SR(eq4s.subs(ux^2 + uy^2 + uz^2,1));
eq4l = eq4l.collect(s);
eq5 = (eq4l == 0); 
pretty_print(eq5);

$$\newcommand{\Bold}[1]{\mathbf{#1}}-r^{2} - 2 \, {\left(\mathit{tx} \mathit{ux} + \mathit{ty} \mathit{uy} + \mathit{tz} \mathit{uz}\right)} s + s^{2} + \mathit{tx}^{2} + \mathit{ty}^{2} + \mathit{tz}^{2} = 0$$

At this point we have a simple quadratic equation to solve in order to find the two possible values of s. That there are two values in general makes geometric sense, they are the points where the ray enters and leaves the sphere. Also, should the values of s be imaginary that is indicating a geometric configuration where the ray does not intersect the sphere.

res = solve(eq5,s);
show(res[0]); show(res[1]);

$$\newcommand{\Bold}[1]{\mathbf{#1}}s = \mathit{tx} \mathit{ux} + \mathit{ty} \mathit{uy} + \mathit{tz} \mathit{uz} - \sqrt{\mathit{tx}^{2} \mathit{ux}^{2} + 2 \, \mathit{tx} \mathit{ty} \mathit{ux} \mathit{uy} + \mathit{ty}^{2} \mathit{uy}^{2} + \mathit{tz}^{2} \mathit{uz}^{2} + r^{2} - \mathit{tx}^{2} - \mathit{ty}^{2} - \mathit{tz}^{2} + 2 \, {\left(\mathit{tx} \mathit{tz} \mathit{ux} + \mathit{ty} \mathit{tz} \mathit{uy}\right)} \mathit{uz}}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}s = \mathit{tx} \mathit{ux} + \mathit{ty} \mathit{uy} + \mathit{tz} \mathit{uz} + \sqrt{\mathit{tx}^{2} \mathit{ux}^{2} + 2 \, \mathit{tx} \mathit{ty} \mathit{ux} \mathit{uy} + \mathit{ty}^{2} \mathit{uy}^{2} + \mathit{tz}^{2} \mathit{uz}^{2} + r^{2} - \mathit{tx}^{2} - \mathit{ty}^{2} - \mathit{tz}^{2} + 2 \, {\left(\mathit{tx} \mathit{tz} \mathit{ux} + \mathit{ty} \mathit{tz} \mathit{uy}\right)} \mathit{uz}}$$

In closing, understanding how this general problem was solved by substituting the X, Y and Z elments of the array into the equation for the sphere is useful. This approach generalizes to other geometric objects defined by implicit forms, and do note we used the implicit form to define our sphere.

Equally important, while you may want to code this solution as a debugging check for you faster algorithm, you will want to use the faster ray sphere intersection being introduced next.