matlab - 在 Matlab 中设置线性编程代码以获得 4D 或更高维度的顶点

Question

假设我有一个 4D 空间中给出的点云，它可以设置为任意密集。这些点将位于多面体的表面上，此时多面体是未知物体。 (只是为了提供一些无用的可视化，这里是我给出的点云的相当任意的投影 - 即绘制 100000x4 数组的前 3 列):

我的目标是获得这个多面体的顶点，而我为获得这些顶点所做的众多尝试之一就是计算点云的凸包。这里的问题是，生成的顶点数量与指定的容差不同，并且没有先验的方法来检查要使用哪个顶点。我还注意到，这里使用 qhull 算法并不是最佳的，因为所有给定的点实际上都位于船体上。现在在我的former question我接受的答案建议设置面部损失函数并使用梯度下降。这似乎是一个不错的方法，但当时我还无法实现它。我认为这是一个线性规划问题，不幸的是我对此一无所知。我扫描了一些在线资源，似乎在该设置中经常遇到找到这样的极值点(这就是凸优化的含义吗？)，这就是为什么我希望提出关于 SO 的后续问题是有意义的.

我也希望如果我将建议的方法粘贴到此处不会违反礼仪:

Let xi in R4, i = 1, ..., m, be the PCA-reduced data points.

Let F = (a, b) be a face, where a is in R4 with a • a = 1 and b is in R.

We define the face loss function L as follows, where λ+, λ- > 0 are parameters you choose. λ+ should be a very small positive number. λ- should be a very large positive number.

L(F) = sumi(λ+ • max(0, a • xi + b) - λ- • min(0, a • xi + b))

We want to find minimal faces F with respect to the loss function L. The small set of all minimal faces will describe your polytope. You can solve for minimal faces by randomly initializing F and then performing gradient descent using the partial derivatives ∂L / ∂aj, j = 1, 2, 3, 4, and ∂L / ∂b. At each step of gradient descent, constrain a • a to be 1 by normalizing.

∂L / ∂aj = sumi(λ+ • xj • [a • xi + b > 0] - λ- • xj • [a • xi + b < 0]) for j = 1, 2, 3, 4

∂L / ∂b = sumi(λ+ • [a • xi + b > 0] - λ- • [a • xi + b < 0])

Note Iverson brackets: [P] = 1 if P is true and [P] = 0 if P is false.

有人可以提供一些建议，我应该如何将其放入代码中吗？我想我可以使用 Matlab 的 linprog 来实现这一点吗？例如，我不太确定随机启动一张脸意味着什么，以及此时 a 和 b 是否未知。

另外，这是一个 link到数据集之一。

Answer 1

以下使用 RANSAC -将 2D 平面拟合到 nD 点云上的方法。然后可以通过这些平面的交集找到顶点。

算法思想:

RANSAC 的总体思想非常简单。

1. 选择称为“假设内点”的随机数据点
2. 根据假设的内点生成模型
3. 找到其他适合模型的点，称之为共识集
4. 重复步骤 1-3，直到找到具有所需量级共识集的模型。
5. 根据最大共识集改进模型。

在我们的示例中，我们可以按以下方式使用它:

1. 选择三个假设的内点。
2. 这三个点所跨越的平面就是我们的模型。
3. 计算整个数据集的点到平面的正交距离。选择预定义容差内的这些点。
4. 重复步骤 1-3 并选择具有最大共识集的平面。
5. (我们现在可以找到适合平面的最小二乘法，但在本例中似乎没有必要。)
6. 保存找到的平面并从我们需要搜索的点中删除其共识集。
7. 重复上述步骤，直到找到所有平面。
8. 使平面相交以获得顶点。

实现此功能所需的代码:

我们需要一些代码来实现这一点，这超出了本答案的范围:

• 计算数据点与模型平面的距离:distancePointsAffineSpace .
• 相交平面以获得边/顶点: affineSpaceIntersection .
• 从最终集中删除重复点:uniquetol 。 (从版本 R2015a 开始，此功能是内置的。)

实现:

有了这个 RANSAC 部分就非常简单了:

%% // Match planes to dataset X:
%  // Choose 3 Points randomly. Generate plane. Find points within tol.
pointsWithinTolOf = @(Points,tol,Space) ...
                      distancePointsAffineSpace(Points, Space)<tol;
availablePoints = 1:size(X,1);
[foundPlane, pointsOfPlane] = deal(cell(0));
for i = 1:maxNumPlanes
    disp(['Plane #',num2str(i)]);
    bestNumInliers = 0;
    for j = 1:numIterations
        randomPointsIdxs = availablePoints(randperm(numel(availablePoints),3));
        currentPlane = X(randomPointsIdxs,:);
        inliers = find(pointsWithinTolOf(X, PointPlaneDistTol, currentPlane));
        numInliers = numel(inliers);
        if numInliers > bestNumInliers
            bestCurrentPlane = currentPlane;
            bestInliers = inliers;
            bestNumInliers = numInliers;
        end
    end
    foundPlane{i} = bestCurrentPlane;
    pointsOfPlane{i} = bestInliers;
    availablePoints = setdiff(availablePoints, bestInliers);
    if isempty(availablePoints)
        break;
    end
end
numPlanes = numel(foundPlane);

其余的代码相当长，所以我只写下基本思想并链接到代码。

找到模型飞机后，我们:

• 通过与所有模型平面相交并检查相交线上是否有数据点来查找每个平面的边缘。
• 通过与每个平面的边相交来计算其顶点。

您可以下载此 here 的完整代码.

可视化:

您似乎正在处理由三个挤压三角形和三个挤压四边形组成的 3D 莫比乌斯带的边界。这是该 strip 的 4D 旋转投影为 3D(投影到 2D 屏幕上)。