Computing the horizon occlusion point 计算地平线遮挡面最佳视点

https://cesium.com/blog/2013/05/09/computing-the-horizon-occlusion-point/

你已经厌倦了地平线筛选了吗?太好了,我也不!
上次，我们解释了什么是地平线剔除，并展示了一个非常有效的方法来测试通过椭球体计算一个被遮挡面遮挡的点。然而，我们想要测试遮挡的对象并不是那么简单。特别是，我们希望能够测试通过椭球体遮挡面遮挡地形瓦片的情况。但地形瓦片是由数千个顶点组成的复杂对象。
Deron Ohlarik在之前的一篇博客中提到了这一点，他解释说，对于任何任意的几何图形，我们都可以计算出一个点的位置，我们称之为地平线遮挡点，它与几何图形有着特殊的关系。不管观察者从哪个方向接近几何图形，在同一时间或在几何图形的任何部分变得可见之前，这个点对观察者来说都是可见的。这正是我们所需要的!然而，如何计算这样一个点的许多细节留给读者作为练习。此外，还不清楚这种方法是否可以推广到椭球体而不是球体。本博客旨在填补这两个空白。
再一次，。
让我们看看我们的情况。和之前一样，我们将所有坐标转换到椭球空间坐标系下，方法是将每个分量(X、Y和Z)乘以沿该轴的椭球半径的倒数。

在该图中，地球显示为蓝色，地形块显示为棕色。在尺度空间中，地球是一个单位球。地形瓦块包围球的中心显示为点 C。边界球不是缩放空间中的球，但这与我们无关，因为我们只使用它的中心。
首先，我们做出一个任意的决定，我们的地平线遮挡点将位于这条中心线|$\vec{O}C$|上，这条中心线是从地球中心到地形瓦片的包围球中心的向量)。我们只需要计算它沿着这个向量的距离。
点V是地形块中的一个顶点。从点V的角度来看，点H是地平线上的一点。从点V的角度来看，有无穷多个水平点，它们在单位球面上形成一个圆，但是只有两个水平点通过与中心线相交形成一个矢量。其中一个显示为实线|$\vec{H}P$| 。另一个以虚线形式连接到点V。在虚线上，与中心线的交点出现在点V之前，因此它将比其他交点更接近椭球的中心，我们不需要关心它。如果点V是地形块中的唯一顶点，那么该图中的点P就是我们的地平线遮挡点。对于多个顶点，我们对每个顶点重复计算点P，然后选择离椭球面最远的一个。
那么我们如何为给定的地形平铺顶点计算点P呢?我们把图中各个角标出来。

上图中角α和β,我们可以表达另一个角度而言,仅仅使用三角形内角和是180度。
根据sin定律:

然后 我们应用二角和差公式：cos(α+β)=cos(α)cos(β)−sin(α)sin(β)

我们已经知道cos(β)， 然后用勾股定理计算sin(β)：

我们用点乘计算 cos(α)

我们可以用叉乘计算sin(α):

最后计算OP距离 |$\vec{O}P$|

注意： 【^OP】 是单位向量，我们可以根据地形瓦片中心点C到椭球中心点O计算出来 。

Cesium 中代码 Source\Core\EllipsoidalOccluder.js：

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

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

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

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