2023-4-19 计算地平线剔除点如图所示，点O为单位球球心，点C为地形瓦片包围球球心，视线是从H->P，我们发现无

在上一节地平线剔除中，实现了通过地平线剔除一个点，这个点的坐标是当作已知的，并且满足这个点被剔除时，整个物体必然被全部剔除，现在我们称该点为地平线剔除点，那么在实际情况中，遇到复杂物体，比如地形瓦片，地平线剔除点是在包围球球心还是包围球表面上或者其他位置？本节专门来讨论这个问题

如图所示，点O为单位球球心，点C为地形瓦片包围球球心，视线是从H->P，我们要求无论相机从哪个方向转过来，在看到地形瓦片之前必然先看到点P,因此P点即为所要计算的地平线剔除点

算法

我们取V为物体上任意一点，可以计算点V所对应的点P，即 HV与OC的交点循环计算瓦片上每一个点，与OC必然计算出多个交点，我们取最远一个，即 $\color{orange}|\vec{OP}|$ 最大的一个

现在明确已知点O,V,需计算P
我们设 $\angle{POV}为 \angle{\alpha}$
设 $\angle{VOH}为\angle{\beta}$
则可以标出其余角，如图所示: 首先根据正弦定理，即对边/对角=2R,R为外接圆半径

\frac{a}{\sin{A}}=\frac{b}{\sin{B}}=\frac{c}{\sin{C}}=2R\\ \frac{|\vec{OP}|}{\sin{(90+\beta)}}=\frac{|\vec{OV}|}{\sin{(90-(\alpha+\beta))}}

根据公式 $\sin{(90+\theta)}=\cos{\theta}$ 得:

\frac{|\vec{OP}|}{\cos{\beta}}=\frac{|\vec{OV}|}{\cos{(\alpha+\beta)}}\\ |\vec{OP}|=\frac{|\vec{OV}|*\cos{\beta}}{\cos{(\alpha+\beta)}}\\ =\frac{|\vec{OV}|*\frac{|\vec{OH}|}{|\vec{OV}|}}{\cos{(\alpha+\beta)}}\\ =\frac{|\vec{OH}|}{\cos{(\alpha+\beta)}}\\ =\frac{1}{\cos{(\alpha+\beta)}}\\ =\frac{1}{\cos{\alpha}*\cos{\beta}-\sin{\alpha}*\sin{\beta}}

现在需要计算 $\cos{\alpha},\cos{\beta},\sin{\alpha},\sin{\beta}$

\sin{\beta}=\frac{|\vec{HV}|}{|\vec{OV}|}\\ =\frac{\sqrt{|\vec{OV}|^2-|\vec{OH}|^2}}{|\vec{OV}|}\\ =\frac{\sqrt{|\vec{OV}|^2-1}}{|\vec{OV}|}\\ \cos{\beta}=\frac{|\vec{OH}|}{|\vec{OV}|}\\ \cos{\alpha}=\hat{OV}\cdot\hat{OP}\\ =\hat{OV}\cdot\hat{OC}\\

由 $|\vec{a}\times\vec{b}|=|\vec{a}||\vec{b}|*\sin{\theta}$ 得

|\hat{OP}\times\vec{OV}|=|\vec{OV}|*\sin{\alpha}\\ 则\sin{\alpha}=\frac{|\hat{OP}\times\vec{OV}|}{|\vec{OV}|}

综上：

\vec{OP}=\hat{OC}*|\vec{OP}|

代码

cesium中实现如下

function computeMagnitude(ellipsoid, position, scaledSpaceDirectionToPoint) {
    var scaledSpacePosition = ellipsoid.transformPositionToScaledSpace(position);
    var magnitudeSquared = scaledSpacePosition.magnitudeSquared();
    var magnitude = Math.sqrt(magnitudeSquared);
    var direction = scaledSpacePosition.divideByScalar(magnitude);

    // For the purpose of this computation, points below the ellipsoid
    // are considered to be on it instead.
    magnitudeSquared = Math.max(1.0, magnitudeSquared);
    magnitude = Math.max(1.0, magnitude);

    var cosAlpha = direction.dot(scaledSpaceDirectionToPoint);
    var sinAlpha = direction.cross(scaledSpaceDirectionToPoint).magnitude();
    var cosBeta = 1.0 / magnitude;
    var sinBeta = Math.sqrt(magnitudeSquared - 1.0) * cosBeta;

    return 1.0 / (cosAlpha * cosBeta - sinAlpha * sinBeta);
}