如何围绕一组2D点拟合边界椭圆 [英] How to fit a bounding ellipse around a set of 2D points

问题描述

给定一组2d点（笛卡尔形式），我需要找到最小面积的椭圆，使得该集中的每个点都位于椭圆上或椭圆内。

我已经在此站点上以伪代码的形式

在我的尝试中，我使用了本征库，用于对矩阵进行各种运算。

  //拟合椭圆
的误差双精度= 0.2; 
 int n = 10； //点数
 int d = 2; //尺寸
 MatrixXd p = MatrixXd :: Random（d，n）; //填充具有随机点的矩阵
 
 MatrixXd q = p; 
 q.conservativeResize（p.rows（）+ 1，p.cols（））; 
 
 for（size_t i = 0; i< q.cols（）; i ++）
 {
 q（q.rows（）-1，i）= 1; 
} 
 
 int计数= 1； 
 double err = 1; 
 
 const double init_u = 1.0 /（double）n; 
 MatrixXd u = MatrixXd :: Constant（n，1，init_u）; 
 
 
而（err>公差）
 {
 MatrixXd Q_tr = q.transpose（）; 
 cout<<  1<<恩德尔
 MatrixXd X = q * u.asDiagonal（）* Q_tr; 
 cout<<  1a<<恩德尔
 MatrixXd M =（Q_tr * X.inverse（）* q）.asDiagonal（）; 
 cout<<  1b<<恩德尔
 
 
 
 int j_x，j_y; 
 double maximum = M.maxCoeff（& j_x，& j_y）; 
 double step_size =（最大值-d-1）/（（d + 1）*（最大值+ 1））; 
 
 MatrixXd new_u =（1-step_size）* u; 
 new_u（j_x，0）+ = step_size; 
 
 cout<<  2<<恩德尔
 
 //查找错误
 MatrixXd u_diff = new_u-u; 
 for（size_t i = 0; i< u_diff.rows（）; i ++）
 {
 for（size_t j = 0; j  u_diff（i，j）* = u_diff（i，j）; //将矩阵
的每个元素平方} 
 err = sqrt（u_diff.sum（））; 
 count ++; 
 u = new_u; 
} 
 
 cout<<  3<<恩德尔
 MatrixXd U = u.asDiagonal（）; 
 MatrixXd A =（1.0 /（double）d）*（p * U * p.transpose（）-（p * u）*（p * u）.transpose（））。inverse（）; 
 MatrixXd c = p * u; 
  

错误发生在以下行：

  MatrixXd M =（Q_tr * X.inverse（）* q）.asDiagonal（）; 
  

，其内容如下：

 运行：/usr/include/eigen3/Eigen/src/Core/DenseBase.h:261：void Eigen :: DenseBase< Derived> :: resize（Eigen :: Index，Eigen ::索引）[具有衍生的=本征::对角线本征::矩阵，0>; Eigen :: Index = long int]：断言`rows == this-> rows（）&& cols == this-> cols（）&&  DenseBase :: resize（）实际上不允许调整大小。’失败。 
已中止（核心已转储）
  

有人可以指出为什么发生此错误，甚至更好，给我一些有关如何使用C ++将椭圆拟合到一组点的建议？

解决方案

使用Eigen，您可以使用 .diagonal（）从矩阵获取对角矢量；您可以使用 .asDiagonal（）将向量视为对角矩阵；但是您不能将稠密矩阵视为对角矩阵。因此，该行应为

  MatrixXd M =（Q_tr * X.inverse（）* q）.diagonal（）; 
  

Given a set of 2d points (in Cartesian form), I need to find the minimum-area ellipse such that every point in the set lies either on or inside the ellipse.

I have found the solution in the form of pseudo-code on this site, but my attempt to implement the solution in C++ has been unsuccessful.

The following image illustrates graphically what the solution to my problem looks like:

In my attempt, I used the Eigen library for the various operations on matrices.

//The tolerance for error in fitting the ellipse
double tolerance = 0.2;
int n = 10; // number of points
int d = 2; // dimension
MatrixXd p = MatrixXd::Random(d,n); //Fill matrix with random points

MatrixXd q = p;
q.conservativeResize(p.rows() + 1, p.cols());

for(size_t i = 0; i < q.cols(); i++)
{
    q(q.rows() - 1, i) = 1;
}

int count = 1;
double err = 1;

const double init_u = 1.0 / (double) n;
MatrixXd u = MatrixXd::Constant(n, 1, init_u);


while(err > tolerance)
{
    MatrixXd Q_tr = q.transpose();
    cout << "1 " << endl;
    MatrixXd X = q * u.asDiagonal() * Q_tr;
    cout << "1a " << endl;
    MatrixXd M = (Q_tr * X.inverse() * q).asDiagonal();
    cout << "1b " << endl;



    int j_x, j_y;
    double maximum = M.maxCoeff(&j_x, &j_y);
    double step_size = (maximum - d - 1) / ((d + 1) * (maximum + 1));

    MatrixXd new_u = (1 - step_size) * u;
    new_u(j_x, 0) += step_size;

    cout << "2 " << endl;

    //Find err
    MatrixXd u_diff = new_u - u;
    for(size_t i = 0; i < u_diff.rows(); i++)
    {
        for(size_t j = 0; j < u_diff.cols(); j++)
            u_diff(i, j) *= u_diff(i, j); // Square each element of the matrix
    }
    err = sqrt(u_diff.sum());
    count++;
    u = new_u;
}

cout << "3 " << endl;
MatrixXd U = u.asDiagonal();
MatrixXd A = (1.0 / (double) d) * (p * U * p.transpose() - (p * u) * (p * u).transpose()).inverse();
MatrixXd c = p * u;

The error occurs on the following line:

MatrixXd M = (Q_tr * X.inverse() * q).asDiagonal();

and it reads as follows:

    run: /usr/include/eigen3/Eigen/src/Core/DenseBase.h:261: void Eigen::DenseBase<Derived>::resize(Eigen::Index, Eigen::Index) [with Derived = Eigen::Diagonal<Eigen::Matrix<double, -1, -1>, 0>; Eigen::Index = long int]: Assertion `rows == this->rows() && cols == this->cols() && "DenseBase::resize() does not actually allow to resize."' failed.
Aborted (core dumped)

Can someone please point out why this error is occurring or even better, give me some advice on how to fit an ellipse to a set of points using C++?

解决方案

With Eigen, you can get the diagonal vector from a matrix with .diagonal(); you can treat a vector as a diagonal matrix with .asDiagonal(); but you cannot treat a dense matrix as a diagonal matrix. So that line should be

MatrixXd M = (Q_tr * X.inverse() * q).diagonal(); 

这篇关于如何围绕一组2D点拟合边界椭圆的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！