检查点是否位于多维椭球内

Question

我无法确定某些点是否位于多维椭球内。 真正的问题是由潜在的转换（旋转、平移）引起的，我不知道如何将它们应用到我的计算中。

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

我决定使用来自here的代码：

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 次）。

根据this：“V 是旋转矩阵，为您提供椭圆体的方向” 我认为这是我的代码需要使用的最后一个元素。

我的方法

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

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

对于每一点：

1. 如下更新位置：点 = 点 - 中心 -> 换句话说，将其转换为椭圆体的中心在点 (0, 0, ... 0)

2. 像椭球平行于所有轴一样旋转它：point = invert(V) * point

3. 使用等式检查点是否位于椭圆体内：(point.x / radius.x)^2 + ... + (point.i / radius.i)^2 其中 i 是维数

4. 如果方程给出的结果

这种方法好吗？我的意思是 - 我不知道是否允许我反转 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

Answer 1

现在这段代码不能打印 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 的结果。