Jak dopasować elipsę wokół zestawu punktów 2D

Question

Biorąc pod uwagę zestaw punktów 2d (w formie kartezjańskiej), muszę znaleźć elipsę o minimalnym obszarze, tak aby każdy punkt w kłamstwie znajdował się na elipsie lub wewnątrz niej .

Mam found the solution w postaci pseudokodu na tej stronie, ale moja próba wdrożenia rozwiązania w C++ zakończyła się niepowodzeniem.

Poniższy rysunek ilustruje graficznie co rozwiązanie mojego problemu wygląda następująco:

W mojej próbie, użyłem biblioteki Eigen dla różnych operacji na macierzach.

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

Błąd występuje w następującym wierszu:

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

i brzmi następująco:

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) 

Może ktoś proszę wskazać dlaczego ten błąd występuje, a nawet lepiej, dać mi kilka porad, jak dopasować elipsę do zestawu punktów za pomocą C++?

Answer 1

W Eigen można uzyskać wektor diagonalny z macierzy z .diagonal(); możesz traktować wektor jako diagonalną macierz z .asDiagonal(); ale nie można traktować gęstej matrycy jako matrycy diagonalnej. Tak więc linia powinna być

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