边界椭圆

Question

我已经获得了图形模块的分配,其中一部分是计算一组任意形状的最小边界椭圆.椭圆不必是轴对齐的.

这是在使用AWT形状的java(euch)中工作,因此我可以使用所有工具形状来检查对象的包含/交集.

Answer 1

您正在寻找最小体积封闭椭球体,或在您的2D情况下,寻找最小面积.该优化问题是凸的并且可以有效地解决.查看我所包含的链接中的MATLAB代码 - 实现是微不足道的,并不需要比矩阵求逆更复杂的东西.

任何对数学感兴趣的人都应该阅读本文档.

此外,绘制椭圆也很简单 - 这可以在这里找到,但是你需要一个MATLAB特定的函数来生成椭圆上的点.

但由于算法以矩阵形式返回椭圆的方程,

矩阵形式http://mathurl.com/yz7flxe.png

您可以使用此代码查看如何将方程式转换为规范形式,

规范http://mathurl.com/y86tlbw.png

使用奇异值分解(SVD).然后使用规范形式绘制椭圆非常容易.

这是MATLAB代码在一组10个随机2D点(蓝色)上的结果.

像PCA这样的其他方法并不能保证从分解中获得的椭圆(本征 /奇异值)将是最小边界椭圆,因为椭圆外的点是方差的指示.

编辑:

因此,如果有人阅读该文档,有两种方法可以在2D中进行此操作:这是最优算法的伪代码 - 在文档中清楚地解释了次优算法:

最优算法:

Input: A 2x10 matrix P storing 10 2D points 
       and tolerance = tolerance for error.
Output: The equation of the ellipse in the matrix form, 
        i.e. a 2x2 matrix A and a 2x1 vector C representing 
        the center of the ellipse.

// Dimension of the points
d = 2;   
// Number of points
N = 10;  

// Add a row of 1s to the 2xN matrix P - so Q is 3xN now.
Q = [P;ones(1,N)]  

// Initialize
count = 1;
err = 1;
//u is an Nx1 vector where each element is 1/N
u = (1/N) * ones(N,1)       

// Khachiyan Algorithm
while err > tolerance
{
    // Matrix multiplication: 
    // diag(u) : if u is a vector, places the elements of u 
    // in the diagonal of an NxN matrix of zeros
    X = Q*diag(u)*Q'; // Q' - transpose of Q    

    // inv(X) returns the matrix inverse of X
    // diag(M) when M is a matrix returns the diagonal vector of M
    M = diag(Q' * inv(X) * Q); // Q' - transpose of Q  

    // Find the value and location of the maximum element in the vector M
    maximum = max(M);
    j = find_maximum_value_location(M);

    // Calculate the step size for the ascent
    step_size = (maximum - d -1)/((d+1)*(maximum-1));

    // Calculate the new_u:
    // Take the vector u, and multiply all the elements in it by (1-step_size)
    new_u = (1 - step_size)*u ;

    // Increment the jth element of new_u by step_size
    new_u(j) = new_u(j) + step_size;

    // Store the error by taking finding the square root of the SSD 
    // between new_u and u
    // The SSD or sum-of-square-differences, takes two vectors 
    // of the same size, creates a new vector by finding the 
    // difference between corresponding elements, squaring 
    // each difference and adding them all together. 

    // So if the vectors were: a = [1 2 3] and b = [5 4 6], then:
    // SSD = (1-5)^2 + (2-4)^2 + (3-6)^2;
    // And the norm(a-b) = sqrt(SSD);
    err = norm(new_u - u);

    // Increment count and replace u
    count = count + 1;
    u = new_u;
}

// Put the elements of the vector u into the diagonal of a matrix
// U with the rest of the elements as 0
U = diag(u);

// Compute the A-matrix
A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' );

// And the center,
c = P * u;