title: Maximum volume inscribed ellipsoid in a polytope/set of points
body: <p><strong>Later Edit</strong>: I uploaded <a href="https://send.firefox.com/download/5897807ea7bddc1e/#H4Zc7WIqCtkCCE3WBxz0og" rel="nofollow noreferrer">here</a> a sample of my original data. It's actually a segmentation image in the DICOM format. The volume of this structure as it is it's ~ 16 mL, so I assume the inner ellipsoid volume should be smaller than that. to extract the points from the DICOM image I used the following code:</p>
  
  <pre><code>import os
  import numpy as np
  import SimpleITK as sitk
  
  
  def get_volume_ml(image):
      x_spacing, y_spacing, z_spacing = image.GetSpacing()
      image_nda = sitk.GetArrayFromImage(image)
      imageSegm_nda_NonZero = image_nda.nonzero()
      num_voxels = len(list(zip(imageSegm_nda_NonZero[0],
                                imageSegm_nda_NonZero[1],
                                imageSegm_nda_NonZero[2])))
      if 0 &gt;= num_voxels:
          print('The mask image does not seem to contain an object.')
          return None
      volume_object_ml = (num_voxels * x_spacing * y_spacing * z_spacing) / 1000
      return volume_object_ml
  
  
  def get_surface_points(folder_path):
      """
      :param folder_path: path to folder where DICOM images are stored
      :return: surface points of the DICOM object
      """
      # DICOM Series
      reader = sitk.ImageSeriesReader()
      dicom_names = reader.GetGDCMSeriesFileNames(os.path.normpath(folder_path))
      reader.SetFileNames(dicom_names)
      reader.MetaDataDictionaryArrayUpdateOn()
      reader.LoadPrivateTagsOn()
      try:
          dcm_img = reader.Execute()
      except Exception:
          print('Non-readable DICOM Data: ', folder_path)
          return None
      volume_obj = get_volume_ml(dcm_img)
      print('The volume of the object in mL:', volume_obj)
      contour = sitk.LabelContour(dcm_img, fullyConnected=False)
      contours = sitk.GetArrayFromImage(contour)
      vertices_locations = contours.nonzero()
  
      vertices_unravel = list(zip(vertices_locations[0], vertices_locations[1], vertices_locations[2]))
      vertices_list = [list(vertices_unravel[i]) for i in range(0, len(vertices_unravel))]
      surface_points = np.array(vertices_list)
  
      return surface_points
  
  folder_path = r"C:\Users\etc\TTT [13]\20160415 114441\Series 052 [CT - Abdomen WT 1 0 I31f 3]"
  points = get_surface_points(folder_path)
  </code></pre>
  
  <p>I have a set of points (n > 1000) in 3D space that describe a hollow ovoid like shape. <strong>What I would like is to fit an ellipsoid (3D) that is inside all of the points.</strong>  I am looking for the maximum volume ellipsoid fitting inside the points.</p>
  
  <p>I tried to adapt the code from <a href="https://stackoverflow.com/questions/14016898/port-matlab-bounding-ellipsoid-code-to-python">Minimum Enclosing Ellipsoid</a> (aka outer bounding ellipsoid)<br>
   by modifying the threshold <code>err &gt; tol</code>, with my logic begin that all points should be smaller than &lt; 1 given the ellipsoid equation. But no success.</p>
  
  <p>I also tried the  Loewner-John adaptation on <a href="https://docs.mosek.com/9.2/pythonfusion/case-studies-ellipsoids.html" rel="nofollow noreferrer">mosek</a>, but I didn't figure how to describe the intersection of a hyperplane with 3D polytope (the Ax &lt;= b representation) so I can use it for the 3D case. So no success again.</p>
  
  <p><a href="https://i.stack.imgur.com/eU2S2.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/eU2S2.png" alt="inscribed ellipsoid"></a></p>
  
  <p>The code from the outer ellipsoid:</p>
  
  <pre><code>import numpy as np
  import numpy.linalg as la
  import matplotlib.pyplot as plt
  from mpl_toolkits.mplot3d import Axes3D
  
  pi = np.pi
  sin = np.sin
  cos = np.cos
  
  def plot_ellipsoid(A, centroid, color, ax):
  """
  
  :param A: matrix
  :param centroid: center
  :param color: color
  :param ax: axis
  :return:
  """
  centroid = np.asarray(centroid)
  A = np.asarray(A)
  U, D, V = la.svd(A)
  rx, ry, rz = 1. / np.sqrt(D)
  u, v = np.mgrid[0:2 * np.pi:20j, -np.pi / 2:np.pi / 2:10j]
  x = rx * np.cos(u) * np.cos(v)
  y = ry * np.sin(u) * np.cos(v)
  z = rz * np.sin(v)
  E = np.dstack((x, y, z))
  E = np.dot(E, V) + centroid
  x, y, z = np.rollaxis(E, axis=-1)
  ax.plot_wireframe(x, y, z, cstride=1, rstride=1, color=color, alpha=0.2)
  ax.set_zlabel('Z-Axis')
  ax.set_ylabel('Y-Axis')
  ax.set_xlabel('X-Axis')
  
  def mvee(points, tol = 0.001):
      """
      Finds the ellipse equation in "center form"
      (x-c).T * A * (x-c) = 1
      """
      N, d = points.shape
      Q = np.column_stack((points, np.ones(N))).T
      err = tol+1.0
      u = np.ones(N)/N
      while err &gt; tol:
          # assert u.sum() == 1 # invariant
          X = np.dot(np.dot(Q, np.diag(u)), Q.T)
          M = np.diag(np.dot(np.dot(Q.T, la.inv(X)), Q))
          jdx = np.argmax(M)
          step_size = (M[jdx]-d-1.0)/((d+1)*(M[jdx]-1.0))
          new_u = (1-step_size)*u
          new_u[jdx] += step_size
          err = la.norm(new_u-u)
          u = new_u
      c = np.dot(u,points)        
      A = la.inv(np.dot(np.dot(points.T, np.diag(u)), points)
                 - np.multiply.outer(c,c))/d
      return A, c
  
  folder_path = r"" # path to a DICOM img folder
  points = get_surface_points(folder_path) # or some random pts 
  
  A, centroid = mvee(points)    
  U, D, V = la.svd(A)    
  rx_outer, ry_outer, rz_outer = 1./np.sqrt(D)
  # PLOT
  fig = plt.figure()
  ax1 = fig.add_subplot(111, projection='3d')
  ax1.scatter(points[:, 0], points[:, 1], points[:, 2], c='blue')
  plot_ellipsoid(A, centroid, 'green', ax1)
  </code></pre>
  
  <p>Which gives me this result for the outer ellipsoid on my sample points:
      <a href="https://i.stack.imgur.com/odiXl.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/odiXl.png" alt="enter image description here"></a></p>
  
  <p><strong>The main question: How do I fit an ellipsoid (3D) inside a cloud of 3D points using Python?</strong>  </p>
  
  <p>Is it possible to modify the algorithm for the outer ellipsoid to get the maximum inscribed (inner) ellipsoid?</p>
  
  <p>I am looking for code in <code>Python</code> ideally.</p>