直线和球的交点

贡献者： addis

预备知识　Matlab 的函数，一元二次方程

　　直线参数方程可以用几何矢量表示为 $ \boldsymbol{\mathbf{P}} _0 + \lambda \hat{\boldsymbol{\mathbf{v}}} $，其中 $ \hat{\boldsymbol{\mathbf{v}}} $ 是单位矢量。球心为 $ \boldsymbol{\mathbf{P}} _c$，半径为 $R$，那么直线和球相交的点满足

\begin{equation} ( \boldsymbol{\mathbf{P}} _0 + \lambda \hat{\boldsymbol{\mathbf{v}}} - \boldsymbol{\mathbf{P}} _c)^2 = R^2~. \end{equation}

其中左边的平方表示矢量和自身的内积，即矢量模长的平方。展开后得到关于 $\lambda$ 的二次方程

\begin{equation} \lambda^2 + 2 \hat{\boldsymbol{\mathbf{v}}} \boldsymbol\cdot ( \boldsymbol{\mathbf{P}} _0 - \boldsymbol{\mathbf{P}} _c) \lambda + ( \boldsymbol{\mathbf{P}} _0 - \boldsymbol{\mathbf{P}} _c)^2 - R^2 = 0~. \end{equation}

直接使用求根公式得

\begin{equation} \lambda = \frac{-b \pm \sqrt{b^2 - 4c}}{2}~. \end{equation}

其中

\begin{equation} b = 2 \boldsymbol{\mathbf{v}} \boldsymbol\cdot ( \boldsymbol{\mathbf{P}} _0 - \boldsymbol{\mathbf{P}} _c)~, \qquad c = ( \boldsymbol{\mathbf{P}} _0 - \boldsymbol{\mathbf{P}} _c)^2 - R^2~. \end{equation}

　　以下给出 Matlab 代码

未完成：应该支持多条直线而不是多个球……

代码 1：LSph2P.m

% 计算直线和球的焦点
% numel(L)=6, size(P)=[N,3], numel(R)=1
function [P1,P2] = LSph2P(L6, P, R)
L6 = L6(:)';
P0 = L6(1:3); v0 = L6(4:6)/norm(L6(4:6));
P0_P = P0-P;
b = 2*P0_P*v0';
c = dot(P0_P, P0_P, 2) - R^2;
d = b.^2 - 4*c;
if d < 0
    P1 = []; P2 = []; return;
end
d = sqrt(d);
P1 = 0.5*(-b -d).*v0 + P0;
if d == 0
    P2 = []; return;
end
P2 = 0.5*(-b +d).*v0 + P0;
end