Как разместить ограничивающий эллипс вокруг набора 2D точек

Question

Учитывая набор 2d точек (в декартовой форме), мне нужно найти эллипс минимальной площади, такой, чтобы каждая точка в множестве лежала либо внутри, либо внутри эллипса.

Я нашел решение в виде псевдокода на этом сайте, но моя попытка реализовать решение в C++ была неудачной.

На следующем рисунке графически показано, как выглядит решение моей проблемы:

В своей попытке я использовал библиотеку Eigen для различных операций над матрицами.

//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;

Ошибка возникает в следующей строке:

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

и это гласит следующее:

    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)

Может кто-нибудь указать, почему эта ошибка происходит или даже лучше, дайте мне несколько советов о том, как подогнать эллипс к набору точек с помощью C++?

Answer 1

С помощью Eigen вы можете получить диагональный вектор из матрицы с .diagonal(); Вы можете рассматривать вектор как диагональную матрицу с .asDiagonal(); но вы не можете рассматривать плотную матрицу как диагональную матрицу. Так что эта строка должна быть

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