检查点是否在多维椭球体内

Question

检查点是否在多维椭球体内

3

我有一个问题，需要确定某些点是否在多维椭球内。真正的问题是由于可能的变换（旋转、平移）引起的，我不知道如何将它们应用到我的计算中。

在我的工作中，我会有一组点，需要找到最小体积包围椭球。然后，我将检查其他一组点，并确定它们中的哪些属于找到的椭球。

我决定使用这里的代码。

function [] = Checkup()
points  = [[ 0.53135758, -0.25818091, -0.32382715] 
    [ 0.58368177, -0.3286576,  -0.23854156,] 
    [ 0.18741533,  0.03066228, -0.94294771] 
    [ 0.65685862, -0.09220681, -0.60347573]
    [ 0.63137604, -0.22978685, -0.27479238]
    [ 0.59683195, -0.15111101, -0.40536606]
    [ 0.68646128,  0.0046802,  -0.68407367]
    [ 0.62311759,  0.0101013,  -0.75863324]];
P = points\'; % <- remove the \ symbol here
dimen = length(P(:,1));

% c - vector with ellipse centers
[A, c] = MinVolEllipse(P, 0.01);
[~, Q, V] = svd(A);
radiuses = 1:dimen;

% Calculate radiuses
for i = 1:dimen
    radiuses(1, i) = 1 / sqrt(Q(i,i));
end

% Check if points lie within ellipse and print Ok for every point inside ellipsoid
for i = 1:length(P(1,:)) % length(P(1,:)) is number of points
    value = 0;
    for j = 1:dimen
        %adding ((p_i - c_i) / (r_i))^2 value
        value = value + ( ( (P(j,i) - c(j, 1) )^2 ) / (radiuses(1, j)^2));

    end
    if value <= 1
        disp('Ok')
    end
end

当前这段代码并不打印“Ok”文本（但它应该打印8次）。

根据此链接所述：“V是给出椭球方向的旋转矩阵”，我认为这是我需要用于我的代码工作的最后一个元素。

我的方法

好的，我对我的代码进行了一些更改。现在我尝试做这样的事情：

假设我有椭球的中心和半径。

对于每个点：

更新位置如下：point = point-center -> 换句话说，将其平移到椭球在点(0, 0,... 0)处的中心。
旋转它，就好像椭球与所有轴平行：point = invert(V) * point
使用等式检查点是否在椭球内：(point.x / radius.x)^2 + ... + (point.i / radius.i)^2，其中i是维数
如果等式结果小于等于1，则点位于椭圆内。

这种方法好吗？我的意思是 - 我不知道是否允许反转V矩阵并将其用作之前旋转的反转...

适当的解决方案

看起来我的提议方法很好，但还有一种更好的方法：

function [result] = Ellipse()
% Returns vector consisting of 10 entries
% that represent error ratio for every test set 

result = 0:9;

% Run tests for all point sets
for pointSet = 0:9
    [Ptraining, Ptest] = LoadPointSet(pointSet);
    count = 0;

    % c - vector with ellipse centers
    [A, c] = MinVolEllipse(Ptraining, 0.00001);

    % Check if points lie within ellipse
    for i = 1:length(Ptest(1,:)) % length(P(1,:)) is number of points
        pp = Ptest(:, i) - c;
        value = pp' * A * pp;    
        if value > 1
            count = count + 1;
        end
    end

    % Get the error ratio
    index = pointSet + 1;
    result(index) = count / length(Ptest(1,:));
end

end

- Waszker

细节在哪里？你使用了什么方程？它是轴对齐的椭球体还是非轴对齐的？你一直在尝试调试什么问题？我会使用视觉检查来验证正确性，例如对一组随机点进行3D绘图，并根据测试椭球体的谓词函数设置它们的颜色...这样你就可以直接看到问题所在...另外，你可以将椭球体绘图添加进去，以查看谓词是否与应该相符。 - Spektre

你不明白我的问题的哪一部分？我试图确定某个点是否属于找到的椭球体。在我提供的Matlab代码中，我为3D空间中的8个点找到了最小体积包围椭球体。可在描述中早先提供的链接中找到可视化结果。 - Waszker

问题不在于任务没有清晰地陈述...问题在于您没有描述您是如何解决这个问题以及您当前的方法有什么问题。如果您按照我提出的方式进行可视化，您将会看到谓词认为什么是椭球体内部，从而更容易推断出问题所在... - Spektre

现在怎么样？我更新了我的方法部分。 - Waszker

1

没有任何理由关闭这个问题。对于那些不懂数学的人，请离开。你们不需要关闭那些我们懂得回答的问题。 - David Hammen

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- David Hammen · Accepted Answer

现在这段代码没有打印出“Ok”文本（但它应该打印8次）。

首先，不应该期望近似算法如MinVolEllipse函数能够给出完全精确的最小体积包围椭圆。您向MinVolEllipse提供了一个公差，因此您应该预期该函数的结果存在一些误差。

更重要的是，在这里您不需要使用奇异值分解（但请参见下面的详细信息）。椭球面的方程式为(x-c)^TA(x-c)=1。您只需检查每个点是否满足(x-c)^TA(x-c)小于一（加上一些公差）即可：

tol_mvee = 0.01;
tol_dist = 0.1;

[A, c] = MinVolEllipse(P, tol_mvee);

% Check if points lie within ellipse and print Ok for every point inside ellipsoid
for i = 1:length(P(1,:)) % length(P(1,:)) is number of points
    x   = P(i,:);
    x_c = x-c;
    d   = dot(x_c, A*x_c);
    if d < 1+tol_dist
        disp('Ok');
    end
end

你的第一个点明显在椭球体内，其余七个点都在真实最小容积包络椭球的边界上或非常靠近边界。算法会遇到一些麻烦，还会遇到一个问题，即所得到的椭球形状明显不是球形的（请参见条件数下面的内容）。

奇异值分解可能对告诉你来自MinVolEllipse的A矩阵是否存在病态很有用，这在这个特定问题中有些情况。您的矩阵的条件数超过200，这意味着您最好不要期望得到比1e-6更小的公差的结果。