在 3D 空间中查找 3D 点和定向椭圆之间的距离（C++）答案

【问题标题】：Find the distance between a 3D point and an Orientated Ellipse in 3D space (C++)在 3D 空间中查找 3D 点和定向椭圆之间的距离（C++）
【发布时间】：2016-07-21 03:56:48
【问题描述】：

为了给这个问题一些背景知识，我正在创建一个游戏，它需要知道一个物体的“轨道”是否在另一个轨道的容差范围内。为了展示这一点，我使用目标轨道绘制了一个具有给定半径（公差）的圆环形状，现在我需要检查椭圆是否在该圆环内。

我在数学/堆栈交换的方程式中迷失了方向，因此寻求更具体的解决方案。为澄清起见，这是一张带有圆环和轨道（红线）的游戏图像。很简单，我想检查那个红色轨道是否在那个圆环形状内。

我认为我需要做的是，在其中一个轨道上的 World-Space 中绘制四个点（很容易做到）。然后我需要计算该点与其他轨道椭圆之间的最短距离。这是困难的部分。有几个例子可以找到点到椭圆的最短距离，但都是二维的，很难理解。

如果该距离小于所有四个点的容差，则认为这等于轨道位于目标圆环内。

为简单起见，所有这些轨道的原点始终位于世界原点 (0, 0, 0) - 我的坐标系是 Z-Up。每个轨道都有一系列定义它的参数（轨道元素）。

【问题讨论】：

如果你知道如何在 2D 中找到点到椭圆的距离，那么在 3D 中它很简单：将点投影到椭圆的平面上，找到投影点在该平面上的距离。现在从椭圆上的目标点，你有一条易于测量的直线到原点。干杯！
椭圆是椭圆还是圆？圆环是圆形的还是椭圆形的？他们有相同的中心吗？ [通常情况下，椭圆轨迹不在行星中心。]
嗨@YvesDaoust - 这两个对象确实是椭圆（不是圆形），然后共享相同的原点。
谢谢@Cheersandhth.-Alf，我最初的尝试确实使用了投影（这很容易做到）——但它找到了我努力的最短距离。知道是否有一些简单的代码示例吗？
@JamesBaxter：我没有这方面的经验。但我认为，为了生成一个对您有用的代码示例，一些关于您如何表示轨道、轨道参数的信息会有所帮助。

标签： c++ math 3d unreal-engine4 orbital-mechanics

【解决方案1】：

这里简单的做法：

将每个轨道采样到一组N 点。

让来自第一个轨道的点为A，来自第二个轨道的点为B。
```
const int N=36;
float A[N][3],B[N][3];
```
找到 2 个最近的点

所以d=|A[i]-B[i]| 是最小的。如果d 小于或等于您的边距/阈值，则轨道彼此太近。
速度与准确性

除非您对 #2 使用一些高级方法，否则它的计算将是 O(N^2)，这有点吓人。 N 越大，结果的准确性越高，但计算时间也越长。有一些方法可以解决这两种情况。例如：
1. 第一个小样本N
2. 当找到最近的点时再次采样两个轨道
  
  但仅在这些点附近（N 更高）。
3. 您可以通过循环 #2 来递归提高准确度，直到达到所需的精度
4. 测试d 椭圆是否彼此太近

【讨论】：

谢谢，这是一个很好的建议。所以对我来说关键是速度，其中有几个需要同时排列，并且每帧都需要检查轨道（每秒 30-60 次）。我将再次尝试上面建议的投影（将一个投影到另一个上，然后尝试在目标上的 4 个点处做最短距离椭圆）
好吧，我已经尝试过了，但不幸的是它太慢了。尤其是与每帧 15 个以上的轨道相比时。要比较的点数开始偏离轨道。对于手机游戏来说，太多了。让我感到困惑的是，在任何地方都没有这样的例子。当然，必须有一种数学方法来计算 3D 空间中点和椭圆之间的最短距离。

【解决方案2】：

我想我可能有一个新的解决方案。

在当前轨道（椭圆）上绘制四个点。
将这些点投影到目标轨道（环面）的平面上。
使用目标轨道倾角作为平面法线，计算每个（归一化）点与近点自变量之间的角度在目标轨道上。
将此角度用作平均异常值，并计算等效偏心异常值。
使用这些偏心异常来绘制目标轨道上的四个点 - 这应该是距离另一个轨道最近的点。
检查这些点之间的距离。

这里的困难在于计算角度并将其转换为另一个轨道上的异常。不过，这应该比递归函数更准确、更快。当我尝试这个时会更新。

编辑：

是的，这行得通！

    // The Four Locations we will use for the checks
TArray<FVector> CurrentOrbit_CheckPositions;
TArray<FVector> TargetOrbit_ProjectedPositions;
CurrentOrbit_CheckPositions.SetNum(4);
TargetOrbit_ProjectedPositions.SetNum(4);

// We first work out the plane of the target orbit.
const FVector Target_LANVector = FVector::ForwardVector.RotateAngleAxis(TargetOrbit.LongitudeAscendingNode, FVector::UpVector); // Vector pointing to Longitude of Ascending Node
const FVector Target_INCVector = FVector::UpVector.RotateAngleAxis(TargetOrbit.Inclination, Target_LANVector);                  // Vector pointing up the inclination axis (orbit normal)
const FVector Target_AOPVector = Target_LANVector.RotateAngleAxis(TargetOrbit.ArgumentOfPeriapsis, Target_INCVector);           // Vector pointing towards the periapse (closest approach)

// Geometric plane of the orbit, using the inclination vector as the normal.
const FPlane ProjectionPlane = FPlane(Target_INCVector, 0.f);   // Plane of the orbit. We only need the 'normal', and the plane origin is the Earths core (periapse focal point)

// Plot four points on the current orbit, using an equally-divided eccentric anomaly.
const float ECCAngle = PI / 2.f;
for (int32 i = 0; i < 4; i++)
{
    // Plot the point, then project it onto the plane
    CurrentOrbit_CheckPositions[i] = PosFromEccAnomaly(i * ECCAngle, CurrentOrbit);
    CurrentOrbit_CheckPositions[i] = FVector::PointPlaneProject(CurrentOrbit_CheckPositions[i], ProjectionPlane);

    // TODO: Distance from the plane is the 'Depth'. If the Depth is > Acceptance Radius, we are outside the torus and can early-out here

    // Normalize the point to find it's direction in world-space (origin in our case is always 0,0,0)
    const FVector PositionDirectionWS = CurrentOrbit_CheckPositions[i].GetSafeNormal();

    // Using the Inclination as the comparison plane - find the angle between the direction of this vector, and the Argument of Periapse vector of the Target orbit
    // TODO: we can probably compute this angle once, using the Periapse vectors from each orbit, and just multiply it by the Index 'I'
    float Angle = FMath::Acos(FVector::DotProduct(PositionDirectionWS, Target_AOPVector));

    // Compute the 'Sign' of the Angle (-180.f - 180.f), using the Cross Product
    const FVector Cross = FVector::CrossProduct(PositionDirectionWS, Target_AOPVector);
    if (FVector::DotProduct(Cross, Target_INCVector) > 0)
    {
        Angle = -Angle;
    }

    // Using the angle directly will give us the position at th eccentric anomaly. We want to take advantage of the Mean Anomaly, and use it as the ecc anomaly
    // We can use this to plot a point on the target orbit, as if it was the eccentric anomaly.
    Angle = Angle - TargetOrbit.Eccentricity * FMathD::Sin(Angle);
    TargetOrbit_ProjectedPositions[i] = PosFromEccAnomaly(Angle, TargetOrbit);}

我希望 cmets 描述它是如何工作的。经过几个月的摸索终于解决了。谢谢大家！

【讨论】：