边界椭圆

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

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

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

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

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

由于算法以矩阵forms返回椭圆的方程,

矩阵formshttp://mathurl.com/yz7flxe.png

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

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

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

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

像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, ie 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(ab) = 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;