使用GJK距离子算法找到简单形最接近原点的点

Question

我正在尝试实现Gilbert-Johnson-Keerthi距离算法（GJK），但我在“距离子算法”中遇到了问题，这也称为“约翰逊算法”，该算法用于确定简单形式上最接近原点的点。我得到了错误的结果，但是我找不到代码中的任何错误，所以问题必须在我的算法解释上。

在约翰逊算法中（如Gino van den Bergen的书《交互式3D环境中的碰撞检测》中所述），最靠近原点的简单形式的仿射包的点由以下公式给出：X = {yi : i ∈ Ix}。

Δi^X值是按X的基数递增的顺序递归确定的：

...并且Δ^X由以下给出：

对于二维，我使用以下方法找到距离原点最近的点：

Point ClosestPointToOrigin(Simplex simplex)
{
    float dx = 0;
    for (int i = 0; i < simplex.size(); ++i)
        dx += dj(simplex, i);

    Point closest_point(0,0);
    for (int i = 0; i < simplex.size(); ++i)
        closest_point += dj(simplex, i) / dx * simplex[i];

    return closest_point;
}

其中Δi的值由以下确定：

float dj(Simplex simplex, int j)
{
    if (j == 0)
    {
        return 1;
    }
    else
    {
        float d = 0;

        for (int i = 0; i < j; ++i)
            d += dj(simplex, i) * (simplex[0] - simplex[j]).dotProduct(simplex[i]);

        return d;
    }
}

对于一个二维单体 X = {y1, y2}，其中y1 = (1,1)，y2 = (1，-1)，上面的代码返回(1.0，-0.333333)，而最接近的点实际上是(1, 0)。

我一定做错了什么，但我无法弄清楚是什么。

- Adrian Lopez

---

也许我错了，但据我所知，在维度n中，单形体应该有n+1个顶点，所以在你的例子中缺少一个顶点，这可能会导致问题？ - Holt

虽然n维单纯形确实有n+1个顶点，但GJK在“单纯形”中包括所有可在所选坐标空间中表示的n-单纯形。例如，在ℝ³中，单纯形可以是四面体（3-单纯形）、三角形（2-单纯形）、线段（1-单纯形）或点（0-单纯形）。因此，上面的线段示例是可以的。尽管如此，感谢您的回复。 - Adrian Lopez

你能发布原始的数学方程吗？因为你提供的代码中有很多未定义的东西，比如 I_x 和 I_X，或者 Delta_i（你定义了 Delta_i^X 但没有定义 Delta_i）。 - Holt

在完整的GJK算法中，I_x代表来自单纯形Y = {y_i : i ∈ I_y}的索引子集I_y。在这种情况下，目标是测试所有开放子集X和找到其中一个包含在aff（X）中的点，该点也在conv（Y）中且最靠近原点。 "Delta_i"和"I_X"是打字错误（应为"Delta_i^X"和"I_x"），我已经进行了更正。 - Adrian Lopez

Answer 1

你的错误在于 dj 函数，可能是你对 dxi 方程的理解有误或者没有书写正确。

我会尝试解释一下，如果你有任何不理解的地方，请随时评论（我将编写伪 Python 代码，但应该很容易理解）。

假设我有以下的单纯形表：

S  = Simplex({
    1: Point (1, 1)  # y1
    2: Point (1,-1)  # y2
})

我可以立即计算出2个Delta值：

然后，我可以计算出另外2个Delta值：

希望现在你能看到你的错误：Δ值是基于索引的，因此对于每个维度为n的Simplex X，您都有n个Δ值。你的一个错误是假设你可以计算出Δ^X₀和Δ^X_i而不考虑X的内容，这是错误的。

现在是最后一个Δ：

请注意：

一旦你到了这里：

这是一个用Python编写的代码，如果你不理解它，我会尝试用你理解的语言来写一个：

import numpy


class Point(numpy.ndarray):
    def __new__(cls, x, y):
        return numpy.asarray([x, y]).astype(float).view(cls)

    def __str__(self):
        return repr(self)

    def __repr__(self):
        return "Point ({}, {})".format(self.x, self.y)

    x = property(fget=lambda s: s[0])
    y = property(fget=lambda s: s[1])


class Simplex(dict):
    def __init__(self, points):
        super(Simplex, self).__init__(enumerate(points))

    def __str__(self):
        return repr(self)

    def __repr__(self):
        return "Simplex <" + dict.__repr__(self) + ">"


def closest_point(s):
    dx = sum(dj(s, i) for i in range(len(s)))
    return sum(dj(s, i) / dx * v for i, v in s.items())


def dj(s, j):
    if len(s) == 0 or (len(s) == 1 and j not in s):
        print(s, j)
        raise ValueError()
    if len(s) == 1:
        return 1
    ts = s.copy()
    yj = s[j]
    del ts[j]
    return sum(
        dj(ts, i) * (ts[list(ts.keys())[0]] - yj).dot(v) 
        for i, v in ts.items()
    )


S = Simplex([Point(1, 1), Point(1, -1)])

print(closest_point(S))

- Holt