title: How to fit a bounding ellipse around a set of 2D points
body: <p>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.</p>
  
  <p>I have <a href="https://stackoverflow.com/a/1768440/4184412">found the solution</a> in the form of pseudo-code on this site, but my attempt to implement the solution in C++ has been unsuccessful.</p>
  
  <p>The following image illustrates graphically what the solution to my problem looks like:</p>
  
  <p><a href="https://i.stack.imgur.com/qv5t2.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qv5t2.png" alt="a set of points with a bounding ellipse"></a></p>
  
  <p>In my attempt, I used the <a href="http://eigen.tuxfamily.org/index.php?title=Main_Page" rel="nofollow noreferrer">Eigen</a> library for the various operations on matrices.</p>
  
  <pre><code>//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 &lt; 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 &gt; tolerance)
  {
      MatrixXd Q_tr = q.transpose();
      cout &lt;&lt; "1 " &lt;&lt; endl;
      MatrixXd X = q * u.asDiagonal() * Q_tr;
      cout &lt;&lt; "1a " &lt;&lt; endl;
      MatrixXd M = (Q_tr * X.inverse() * q).asDiagonal();
      cout &lt;&lt; "1b " &lt;&lt; endl;
  
  
  
      int j_x, j_y;
      double maximum = M.maxCoeff(&amp;j_x, &amp;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 &lt;&lt; "2 " &lt;&lt; endl;
  
      //Find err
      MatrixXd u_diff = new_u - u;
      for(size_t i = 0; i &lt; u_diff.rows(); i++)
      {
          for(size_t j = 0; j &lt; 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 &lt;&lt; "3 " &lt;&lt; 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;
  </code></pre>
  
  <p>The error occurs on the following line:</p>
  
  <pre><code>MatrixXd M = (Q_tr * X.inverse() * q).asDiagonal();
  </code></pre>
  
  <p>and it reads as follows:</p>
  
  <pre><code>    run: /usr/include/eigen3/Eigen/src/Core/DenseBase.h:261: void Eigen::DenseBase&lt;Derived&gt;::resize(Eigen::Index, Eigen::Index) [with Derived = Eigen::Diagonal&lt;Eigen::Matrix&lt;double, -1, -1&gt;, 0&gt;; Eigen::Index = long int]: Assertion `rows == this-&gt;rows() &amp;&amp; cols == this-&gt;cols() &amp;&amp; "DenseBase::resize() does not actually allow to resize."' failed.
  Aborted (core dumped)
  </code></pre>
  
  <p>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++?</p>