Maximum volume inscribed ellipsoid in a polytope/set of points
Asked Answered
U

4

12

Later Edit: I uploaded here 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:

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 >= 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)

I have a set of points (n > 1000) in 3D space that describe a hollow ovoid like shape. What I would like is to fit an ellipsoid (3D) that is inside all of the points. I am looking for the maximum volume ellipsoid fitting inside the points.

I tried to adapt the code from Minimum Enclosing Ellipsoid (aka outer bounding ellipsoid)
by modifying the threshold err > tol, with my logic begin that all points should be smaller than < 1 given the ellipsoid equation. But no success.

I also tried the Loewner-John adaptation on mosek, but I didn't figure how to describe the intersection of a hyperplane with 3D polytope (the Ax <= b representation) so I can use it for the 3D case. So no success again.

inscribed ellipsoid

The code from the outer ellipsoid:

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 > 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)

Which gives me this result for the outer ellipsoid on my sample points: enter image description here

The main question: How do I fit an ellipsoid (3D) inside a cloud of 3D points using Python?

Is it possible to modify the algorithm for the outer ellipsoid to get the maximum inscribed (inner) ellipsoid?

I am looking for code in Python ideally.

Underlaid answered 17/5, 2020 at 21:31 Comment(13)
so what is the question? I'm not sure I see a question here.Wende
@Wende the question is. how do I modify the algorithm for the outer ellipsoid to get the maximum inscribed ellipsoid?Underlaid
Could you perhaps start with the containing ellipsoid and then shrink it? It's not clear that will give the optimal answer but it might give an answer. Strictly speaking I don't think the problem is precise: you need to say what ypu mean by 'inside'.Dish
@Dish I think the problem cannot be precisely defined as it is. What I am looking for is the maximum volume that could fit inside those points (without touching them). However, I think there is no computationally efficient ( or fast way) to prove that any inscribed volume is the maximum one. so, yes what you suggest makes sense for the moment. how should I go about shrinking the outer ellipsoid and how do I check it's inside the points?Underlaid
To clarify your question "ellipsoid is inside a point cloud", do you mean that the ellipsoid is contained within the convex hull of the point cloud? Or do you mean that the ellipsoid should satisfy two requirements: 1. the ellipsoid is contained within the convex hull of the point cloud. 2. none of the point in the point cloud is contained within the ellipsoid?Godderd
Not exactly the answer you are looking for but I have seen this article that solves a similar problem in LocalSolver (has Python API) localsolver.com/news.html?id=96 . Maybe it will help.Nucleotide
It'd be handy if you could edit your point extraction code to show an MWE of loading the data from the example set you uploaded.Herzog
@Herzog not sure I know what an MWE is, but I added and tested the full point extraction for the dataset I uploaded. I also updated the plot to reflect the points from the sample data.Underlaid
@Roxanne: It's a minimum, working example. I'm on the road now, but will try to take another look in the next few days.Herzog
Can any assumption be made on the input points?Blumenthal
@Blumenthal in theory yes, give me an example of assumption. the points come from a 3D segmentation of an organ.Underlaid
For example, can you assume that center of mass of the points will be in the desired ellipsoid?Blumenthal
@Blumenthal certainly, that's a very accurate assumption and how it should be according the medical theory behind the treatment. the points, inner and outer ellipsoid should have the same center of mass.Underlaid
G
9

Problem statement

Given a number of points v₁, v₂, ..., vₙ, find a large ellipsoid satisfying two constraints:

  1. The ellipsoid is in the convex hull ℋ = ConvexHull(v₁, v₂, ..., vₙ).
  2. None of the points v₁, v₂, ..., vₙ is within the ellipsoid.

I propose an iterative procedure to find a large ellipsoid satisfying these two constraints. In each iteration we need to solve a semidefinite programming problem. This iterative procedure is guaranteed to converge, however it not guaranteed to converge to the globally largest ellipsoid.

Approach

Find a single ellipsoid

The core of our iterative procedure is that in each iteration, we find an ellipsoid satisfying 3 conditions:

  • The ellipsoid is contained within ConvexHull(v₁, v₂, ..., vₙ) = {x | Ax<=b}.
  • For a set of points u₁, ... uₘ (where {v₁, v₂, ..., vₙ} ⊂ {u₁, ... uₘ}, namely the given point in the point clouds belongs to this set of points u₁, ... uₘ), the ellipsoid doesn't contain any point in u₁, ... uₘ. We call this set u₁, ... uₘ as "outside points".
  • For a set of points w₁,..., wₖ (where {w₁,..., wₖ} ∩ {v₁, v₂, ..., vₙ} = ∅, namely none of the point in v₁, v₂, ..., vₙ belongs to {w₁,..., wₖ}), the ellipsoid contains all of the points w₁,..., wₖ. We call this set w₁,..., wₖ as "inside points".

The intuitive idea is that the "inside points" w₁,..., wₖ indicate the volume of the ellipsoid. We will append new point to "inside points" so as to increase the ellipsoid volume.

To find such an ellipsoid through convex optimization, we parameterize the ellipsoid as

{x | xᵀPx + 2qᵀx  ≤ r}

and we will search for P, q, r.

The condition that the "outside points" u₁, ... uₘ are all outside of the ellipsoid is formulated as

uᵢᵀPuᵢ + 2qᵀuᵢ >= r ∀ i=1, ..., m

this is a linear constraint on P, q, r.

The condition that the "inside points" w₁,..., wₖ are all inside the ellipsoid is formulated as

wᵢᵀPwᵢ + 2qᵀwᵢ <= r ∀ i=1, ..., k

This is also a linear constraint on P, q, r.

We also impose the constraint

P is positive definite

P being positive definite, together with the constraint that there exists points wᵢ satisfying wᵢᵀPwᵢ + 2qᵀwᵢ <= r guarantees that the set {x | xᵀPx + 2qᵀx ≤ r} is an ellipsoid.

We also have the constraint that the ellipsoid is inside the convex hull ℋ={x | aᵢᵀx≤ bᵢ, i=1,...,l} (namely there are l halfspaces as the H-representation of ℋ). Using s-lemma, we know that a necessary and sufficient condition for the halfspace {x|aᵢᵀx≤ bᵢ} containing the ellipsoid is that

∃ λᵢ >= 0,
s.t [P            q -λᵢaᵢ/2]  is positive semidefinite.
    [(q-λᵢaᵢ/2)ᵀ     λᵢbᵢ-r]

Hence we can solve the following semidefinite programming problem to find the ellipsoid that contains all the "inside points", doesn't contain any "outside points", and is within the convex hull ℋ

find P, q, r, λ
s.t uᵢᵀPuᵢ + 2qᵀuᵢ >= r ∀ i=1, ..., m
    wᵢᵀPwᵢ + 2qᵀwᵢ <= r ∀ i=1, ..., k
    P is positive definite.
    λ >= 0,
    [P            q -λᵢaᵢ/2]  is positive semidefinite.
    [(q-λᵢaᵢ/2)ᵀ     λᵢbᵢ-r]

We call this P, q, r = find_ellipsoid(outside_points, inside_points, A, b).

The volume of this ellipsoid is proportional to (r + qᵀP⁻¹q)/power(det(P), 1/3).

Iterative procedure.

We initialize "outside points" as all the points v₁, v₂, ..., vₙ in the point cloud, and "inside points" as a single point w₁ in the convex hull ℋ. In each iteration, we use find_ellipsoid function in the previous sub-section to find the ellipsoid within ℋ that contains all "inside points" but doesn't contain any "outside points". Depending on the result of the SDP in find_ellipsoid, we do the following

  • If the SDP is feasible. We then compare the newly found ellipsoid with the largest ellipsoid found so far. If this new ellipsoid is larger, then accept it as the largest ellipsoid found so far.
  • If the SDP is infeasible, then we remove the last point in "inside points", add this point to "outside point".

In both cases, we then take a new sample point in the convex hull ℋ, add that sample point to "inside points", and then solve the SDP again.

The complete algorithm is as follows

  1. Initialize "outside points" to v₁, v₂, ..., vₙ, initialize "inside points" to a single random point in the convex hull ℋ.
  2. while iter < max_iterations:
  3. Solve the SDP P, q, r = find_ellipsoid(outside_points, inside_points, A, b).
  4. If SDP is feasible and volume(Ellipsoid(P, q, r)) > largest_volume, set P_best = P, q_best=q, r_best = r.
  5. If SDP is infeasible, pt = inside_points.pop_last(), outside_points.push_back(pt).
  6. Randomly sample a new point in ℋ, append the point to "inside points", iter += 1. Go to step 3.

Code

from scipy.spatial import ConvexHull, Delaunay
import scipy
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import dirichlet
from mpl_toolkits.mplot3d import Axes3D  # noqa


def get_hull(pts):
    dim = pts.shape[1]
    hull = ConvexHull(pts)
    A = hull.equations[:, 0:dim]
    b = hull.equations[:, dim]
    return A, -b, hull


def compute_ellipsoid_volume(P, q, r):
    """
    The volume of the ellipsoid xᵀPx + 2qᵀx ≤ r is proportional to
    power(r + qᵀP⁻¹q, dim/2)/sqrt(det(P))
    We return this number.
    """
    dim = P.shape[0]
    return np.power(r + q @ np.linalg.solve(P, q)), dim/2) / \
        np.sqrt(np.linalg.det(P))


def uniform_sample_from_convex_hull(deln, dim, n):
    """
    Uniformly sample n points in the convex hull Ax<=b
    This is copied from
    https://mcmap.net/q/1009000/-how-to-get-uniformly-distributed-points-in-convex-hull
    @param deln Delaunay of the convex hull.
    """
    vols = np.abs(np.linalg.det(deln[:, :dim, :] - deln[:, dim:, :]))\
        / np.math.factorial(dim)
    sample = np.random.choice(len(vols), size=n, p=vols / vols.sum())

    return np.einsum('ijk, ij -> ik', deln[sample],
                     dirichlet.rvs([1]*(dim + 1), size=n))


def centered_sample_from_convex_hull(pts):
    """
    Sample a random point z that is in the convex hull of the points
    v₁, ..., vₙ. z = (w₁v₁ + ... + wₙvₙ) / (w₁ + ... + wₙ) where wᵢ are all
    uniformly sampled from [0, 1]. Notice that by central limit theorem, the
    distribution of this sample is centered around the convex hull center, and
    also with small variance when the number of points are large.
    """
    num_pts = pts.shape[0]
    pts_weights = np.random.uniform(0, 1, num_pts)
    z = (pts_weights @ pts) / np.sum(pts_weights)
    return z


def find_ellipsoid(outside_pts, inside_pts, A, b):
    """
    For a given sets of points v₁, ..., vₙ, find the ellipsoid satisfying
    three constraints:
    1. The ellipsoid is within the convex hull of these points.
    2. The ellipsoid doesn't contain any of the points.
    3. The ellipsoid contains all the points in @p inside_pts
    This ellipsoid is parameterized as {x | xᵀPx + 2qᵀx ≤ r }.
    We find this ellipsoid by solving a semidefinite programming problem.
    @param outside_pts outside_pts[i, :] is the i'th point vᵢ. The point vᵢ
    must be outside of the ellipsoid.
    @param inside_pts inside_pts[i, :] is the i'th point that must be inside
    the ellipsoid.
    @param A, b The convex hull of v₁, ..., vₙ is Ax<=b
    @return (P, q, r, λ) P, q, r are the parameterization of this ellipsoid. λ
    is the slack variable used in constraining the ellipsoid inside the convex
    hull Ax <= b. If the problem is infeasible, then returns
    None, None, None, None
    """
    assert(isinstance(outside_pts, np.ndarray))
    (num_outside_pts, dim) = outside_pts.shape
    assert(isinstance(inside_pts, np.ndarray))
    assert(inside_pts.shape[1] == dim)
    num_inside_pts = inside_pts.shape[0]

    constraints = []
    P = cp.Variable((dim, dim), symmetric=True)
    q = cp.Variable(dim)
    r = cp.Variable()

    # Impose the constraint that v₁, ..., vₙ are all outside of the ellipsoid.
    for i in range(num_outside_pts):
        constraints.append(
            outside_pts[i, :] @ (P @ outside_pts[i, :]) +
            2 * q @ outside_pts[i, :] >= r)
    # P is strictly positive definite.
    epsilon = 1e-6
    constraints.append(P - epsilon * np.eye(dim) >> 0)

    # Add the constraint that the ellipsoid contains @p inside_pts.
    for i in range(num_inside_pts):
        constraints.append(
            inside_pts[i, :] @ (P @ inside_pts[i, :]) +
            2 * q @ inside_pts[i, :] <= r)

    # Now add the constraint that the ellipsoid is in the convex hull Ax<=b.
    # Using s-lemma, we know that the constraint is
    # ∃ λᵢ > 0,
    # s.t [P            q -λᵢaᵢ/2]  is positive semidefinite.
    #     [(q-λᵢaᵢ/2)ᵀ     λᵢbᵢ-r]
    num_faces = A.shape[0]
    lambda_var = cp.Variable(num_faces)
    constraints.append(lambda_var >= 0)
    Q = [None] * num_faces
    for i in range(num_faces):
        Q[i] = cp.Variable((dim+1, dim+1), PSD=True)
        constraints.append(Q[i][:dim, :dim] == P)
        constraints.append(Q[i][:dim, dim] == q - lambda_var[i] * A[i, :]/2)
        constraints.append(Q[i][-1, -1] == lambda_var[i] * b[i] - r)

    prob = cp.Problem(cp.Minimize(0), constraints)
    try:
        prob.solve(verbose=False)
    except cp.error.SolverError:
        return None, None, None, None

    if prob.status == 'optimal':
        P_val = P.value
        q_val = q.value
        r_val = r.value
        lambda_val = lambda_var.value
        return P_val, q_val, r_val, lambda_val
    else:
        return None, None, None, None


def draw_ellipsoid(P, q, r, outside_pts, inside_pts):
    """
    Draw an ellipsoid defined as {x | xᵀPx + 2qᵀx ≤ r }
    This ellipsoid is equivalent to
    |Lx + L⁻¹q| ≤ √(r + qᵀP⁻¹q)
    where L is the symmetric matrix satisfying L * L = P
    """
    fig = plt.figure()
    dim = P.shape[0]
    L = scipy.linalg.sqrtm(P)
    radius = np.sqrt(r + q@(np.linalg.solve(P, q)))
    if dim == 2:
        # first compute the points on the unit sphere
        theta = np.linspace(0, 2 * np.pi, 200)
        sphere_pts = np.vstack((np.cos(theta), np.sin(theta)))
        ellipsoid_pts = np.linalg.solve(
            L, radius * sphere_pts - (np.linalg.solve(L, q)).reshape((2, -1)))
        ax = fig.add_subplot(111)
        ax.plot(ellipsoid_pts[0, :], ellipsoid_pts[1, :], c='blue')

        ax.scatter(outside_pts[:, 0], outside_pts[:, 1], c='red')
        ax.scatter(inside_pts[:, 0], inside_pts[:, 1], s=20, c='green')
        ax.axis('equal')
        plt.show()
    if dim == 3:
        u = np.linspace(0, np.pi, 30)
        v = np.linspace(0, 2*np.pi, 30)

        sphere_pts_x = np.outer(np.sin(u), np.sin(v))
        sphere_pts_y = np.outer(np.sin(u), np.cos(v))
        sphere_pts_z = np.outer(np.cos(u), np.ones_like(v))
        sphere_pts = np.vstack((
            sphere_pts_x.reshape((1, -1)), sphere_pts_y.reshape((1, -1)),
            sphere_pts_z.reshape((1, -1))))
        ellipsoid_pts = np.linalg.solve(
            L, radius * sphere_pts - (np.linalg.solve(L, q)).reshape((3, -1)))
        ax = plt.axes(projection='3d')
        ellipsoid_pts_x = ellipsoid_pts[0, :].reshape(sphere_pts_x.shape)
        ellipsoid_pts_y = ellipsoid_pts[1, :].reshape(sphere_pts_y.shape)
        ellipsoid_pts_z = ellipsoid_pts[2, :].reshape(sphere_pts_z.shape)
        ax.plot_wireframe(ellipsoid_pts_x, ellipsoid_pts_y, ellipsoid_pts_z)
        ax.scatter(outside_pts[:, 0], outside_pts[:, 1], outside_pts[:, 2],
                   c='red')
        ax.scatter(inside_pts[:, 0], inside_pts[:, 1], inside_pts[:, 2], s=20,
                   c='green')
        ax.axis('equal')
        plt.show()


def find_large_ellipsoid(pts, max_iterations):
    """
    We find a large ellipsoid within the convex hull of @p pts but not
    containing any point in @p pts.
    The algorithm proceeds iteratively
    1. Start with outside_pts = pts, inside_pts = z where z is a random point
       in the convex hull of @p outside_pts.
    2. while num_iter < max_iterations
    3.   Solve an SDP to find an ellipsoid that is within the convex hull of
         @p pts, not containing any outside_pts, but contains all inside_pts.
    4.   If the SDP in the previous step is infeasible, then remove z from
         inside_pts, and append it to the outside_pts.
    5.   Randomly sample a point in the convex hull of @p pts, if this point is
         outside of the current ellipsoid, then append it to inside_pts.
    6.   num_iter += 1
    When the iterations limit is reached, we report the ellipsoid with the
    maximal volume.
    @param pts pts[i, :] is the i'th points that has to be outside of the
    ellipsoid.
    @param max_iterations The iterations limit.
    @return (P, q, r) The largest ellipsoid is parameterized as
    {x | xᵀPx + 2qᵀx ≤ r }
    """
    dim = pts.shape[1]
    A, b, hull = get_hull(pts)
    hull_vertices = pts[hull.vertices]
    deln = hull_vertices[Delaunay(hull_vertices).simplices]

    outside_pts = pts
    z = centered_sample_from_convex_hull(pts)
    inside_pts = z.reshape((1, -1))

    num_iter = 0
    max_ellipsoid_volume = -np.inf
    while num_iter < max_iterations:
        (P, q, r, lambda_val) = find_ellipsoid(outside_pts, inside_pts, A, b)
        if P is not None:
            volume = compute_ellipsoid_volume(P, q, r)
            if volume > max_ellipsoid_volume:
                max_ellipsoid_volume = volume
                P_best = P
                q_best = q
                r_best = r
            else:
                # Adding the last inside_pts doesn't increase the ellipsoid
                # volume, so remove it.
                inside_pts = inside_pts[:-1, :]
        else:
            outside_pts = np.vstack((outside_pts, inside_pts[-1, :]))
            inside_pts = inside_pts[:-1, :]

        # Now take a new sample that is outside of the ellipsoid.
        sample_pts = uniform_sample_from_convex_hull(deln, dim, 20)
        is_in_ellipsoid = np.sum(sample_pts.T*(P_best @ sample_pts.T), axis=0)\
            + 2 * sample_pts @ q_best <= r_best
        if np.all(is_in_ellipsoid):
            # All the sampled points are in the ellipsoid, the ellipsoid is
            # already large enough.
            return P_best, q_best, r_best
        else:
            inside_pts = np.vstack((
                inside_pts, sample_pts[np.where(~is_in_ellipsoid)[0][0], :]))
            num_iter += 1

    return P_best, q_best, r_best

if __name__ == "__main__":
    pts = np.array([[0., 0.], [0., 1.], [1., 1.], [1., 0.], [0.2, 0.4]])
    max_iterations = 10
    P, q, r = find_large_ellipsoid(pts, max_iterations)

I also put the code in the github repo

Results

Here is the result on a simple 2D example, say we want to find a large ellipsoid that doesn't contain the five red points in the figure below. Here is the result after the first iteration iteration1_result. The red points are the "outside points" (the initial outside points are v₁, v₂, ..., vₙ), the green point is the initial "inside points".

In the second iteration, the ellipsoid becomes

iteration2_result. By adding one more "inside point" (green dot), the ellipsoid gets larger.

This gif shows the animation of the 10 iteations.all_iterations_result

Godderd answered 20/5, 2020 at 5:34 Comment(7)
thanks, this is a very in-depth mathematical explanation. however what I struggle with is bringing the maths into actual code, since it seems I might need some very specific obscure libraries which are not your standard Python librariesUnderlaid
@Roxanne, I can write the python code for the approach. This approach should only use cvxpy (I later realized I don't need pycddlib). Before I write the code, could you please clarify my question: for all the points in the point cloud, do you want the ellipsoid not to contain any of the point? The approach above assumes that none of the point in the point cloud is within the ellipsoid.Godderd
yes, exactly. for all the points in the cloud, the inner ellipsoid should not to contain any of the points. it's just that my points have this really weird shape. see attached example of sample data points in the updated original questions.Underlaid
@Roxane, I just updated my approach. This new approach solves the problem in an iterative manner. Each iteration we need to solve a semidefinite programming problem (and I recommend using Mosek instead of SCS for this problem). Also I put the code both here and on the github repo.Godderd
@Roxanne, BTW, the proposed approach is not deterministic, as it need to sample "inside points". So you could run the procedure for multiple times, and pick the result with the largest volume.Godderd
I got this error: Traceback (most recent call last): File "......................py", line 251, in <module> P, q, r = find_large_ellipsoid(pts, max_iterations) File "....................py", line 235, in find_large_ellipsoid is_in_ellipsoid = np.sum(sample_pts.T*(P_best @ sample_pts.T), axis=0)\ UnboundLocalError: local variable 'P_best' referenced before assignmentAbolition
@seghier. From reading the code, P_best is unset either you have max_iterations <= 0, or volume returned from volume = compute_ellipsoid_volume(P, q, r) is NAN. Could you check the returned value of volume? BTW, we also have a github repo github.com/hongkai-dai/large_inscribed_ellipsoid, it would be easier to discuss bug by filing an issue on that repo.Godderd
H
8

Whether this answer works depends on how much noise is in your data. The idea is to first find the point cloud's convex hull and then find the largest ellipsoid which fits within this hull. If most of your points lie close to the surface of the ellipsoid they describe than this approximation won't be "too bad".

To do so, note that a convex hull can be described by a set of linear inequalities Ax<=b.

Note that the bounding ellipsoid can be described by E={Bx+d for ||x||_2<=1}, where B is a positive semi-definite matrix describing how and which directions the ellipsoid is stretched and d is a vector describing its offset.

Note that the volume of the ellipsoid is given by det(B^-1). If we were to try to maximize or minimize this determinant we would fail because that would give a non-convex problem. However, applying a log transform log(det(B^-1)) makes the problem convex again. The optimization program we're going to use doesn't allow matrix inverses, but it's easy to show that the foregoing is equivalent to -log(det(B)).

Finally, some bracing algebraic manipulation gives us the optimization problem:

minimize -log(det(B))
s.t.     ||B*A_i||_2 + a_i^T * d <= b_i, i = 1, ..., m
         B is PSD

We can solve this in Python using CVXPY as follows:

#!/usr/bin/env python3

from mpl_toolkits.mplot3d import axes3d
from scipy.spatial import ConvexHull
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
import sklearn.datasets

#From: https://mcmap.net/q/1009002/-how-to-generate-a-random-sample-of-points-from-a-3-d-ellipsoid-using-python
def random_point_ellipsoid(a,b,c,x0,y0,z0):
    """Generate a random point on an ellipsoid defined by a,b,c"""
    u = np.random.rand()
    v = np.random.rand()
    theta = u * 2.0 * np.pi
    phi = np.arccos(2.0 * v - 1.0)
    sinTheta = np.sin(theta);
    cosTheta = np.cos(theta);
    sinPhi = np.sin(phi);
    cosPhi = np.cos(phi);
    rx = a * sinPhi * cosTheta;
    ry = b * sinPhi * sinTheta;
    rz = c * cosPhi;
    return rx, ry, rz

def random_point_ellipse(W,d):
  # random angle
  alpha = 2 * np.pi * np.random.random()
  # vector on that angle
  pt = np.array([np.cos(alpha),np.sin(alpha)])
  # Ellipsoidize it
  return W@pt+d

def GetRandom(dims, Npts):
  if dims==2:
    W = sklearn.datasets.make_spd_matrix(2)
    d = np.array([2,3])
    points = np.array([random_point_ellipse(W,d) for i in range(Npts)])
  elif dims==3:
    points = np.array([random_point_ellipsoid(3,5,7,2,3,3) for i in range(Npts)])
  else:
    raise Exception("dims must be 2 or 3!")
  noise = np.random.multivariate_normal(mean=[0]*dims, cov=0.2*np.eye(dims), size=Npts)
  return points+noise

def GetHull(points):
  dim  = points.shape[1]
  hull = ConvexHull(points)
  A    = hull.equations[:,0:dim]
  b    = hull.equations[:,dim]
  return A, -b, hull #Negative moves b to the RHS of the inequality

def Plot(points, hull, B, d):
  fig = plt.figure()
  if points.shape[1]==2:
    ax = fig.add_subplot(111)
    ax.scatter(points[:,0], points[:,1])
    for simplex in hull.simplices:
      plt.plot(points[simplex, 0], points[simplex, 1], 'k-')
    display_points = np.array([random_point_ellipse([[1,0],[0,1]],[0,0]) for i in range(100)])
    display_points = display_points@B+d
    ax.scatter(display_points[:,0], display_points[:,1])
  elif points.shape[1]==3:
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(points[:,0], points[:,1], points[:,2])
    display_points = np.array([random_point_ellipsoid(1,1,1,0,0,0) for i in range(1000)])
    display_points = display_points@B+d
    ax.scatter(display_points[:,0], display_points[:,1], display_points[:,2])
  plt.show()

def FindMaximumVolumeInscribedEllipsoid(points):
  """Find the inscribed ellipsoid of maximum volume. Return its matrix-offset form."""
  dim = points.shape[1]
  A,b,hull = GetHull(points)

  B = cp.Variable((dim,dim), PSD=True) #Ellipsoid
  d = cp.Variable(dim)                 #Center

  constraints = [cp.norm(B@A[i],2)+A[i]@d<=b[i] for i in range(len(A))]
  prob = cp.Problem(cp.Minimize(-cp.log_det(B)), constraints)
  optval = prob.solve()
  if optval==np.inf:
    raise Exception("No solution possible!")
  print(f"Optimal value: {optval}") 

  Plot(points, hull, B.value, d.value)

  return B.value, d.value

FindMaximumVolumeInscribedEllipsoid(GetRandom(dims=2, Npts=100))
FindMaximumVolumeInscribedEllipsoid(GetRandom(dims=3, Npts=100))

Solutions are computed quickly.

Visually, this gives (for 2D):

2D maximum volume inscribed ellipsoid

Note that I've added a lot of noise to emphasize what's going on.

And for 3D:

3D maximum volume inscribed ellipsoid

Although the code above is written for two or three dimensions, you could easily adapt it for any number of dimensions, though the visualization will become more difficult.

If the convex hull is not good and you want some kind of "interior convex hull", that'll be harder: this hull isn't well-defined. However, you could use alpha shapes to try to find such a hull and then use the algorithm above to solve for it.

Note also that since we're using a convex polytope to bound the ellipse, rather than the points themselves, even if the points perfectly described an ellipsoid, we end up with an underestimated volume. We can visualize this, as below:

Circle inscribed in a square

If the vertices of the square are the points, then the square is their convex hull. The circle bounded by the hull is clearly smaller than the circle that would be bounded only by the points.

EDIT: TO get the volume, you need to convert pixel indices to the coordinate system of your DICOM image, like so (NOTE: I'm not sure if I've scaled the correct coordinates by the correct values, but you'll be able to figure this out given your knowledge of the data):

from mpl_toolkits.mplot3d import axes3d
from scipy.spatial import ConvexHull
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
import os
import sklearn.datasets
import SimpleITK as sitk
import code



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 >= 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(dcm_img):
    x_spacing, y_spacing, z_spacing = dcm_img.GetSpacing()
    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)

    surface_points = surface_points.astype(np.float64)

    surface_points[:,0] *= x_spacing/10
    surface_points[:,1] *= y_spacing/10
    surface_points[:,2] *= z_spacing/10

    return surface_points

def get_dcm_image(folder_path):
    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:
        raise Exception('Non-readable DICOM Data: ', folder_path)

    return dcm_img

def GetHull(points):
  dim  = points.shape[1]
  hull = ConvexHull(points)
  A    = hull.equations[:,0:dim]
  b    = hull.equations[:,dim]
  return A, -b, hull #Negative moves b to the RHS of the inequality

def FindMaximumVolumeInscribedEllipsoid(points):
  """Find the inscribed ellipsoid of maximum volume. Return its matrix-offset form."""
  dim = points.shape[1]
  A,b,hull = GetHull(points)

  B = cp.Variable((dim,dim), PSD=True) #Ellipsoid
  d = cp.Variable(dim)                 #Center

  constraints = [cp.norm(B@A[i],2)+A[i]@d<=b[i] for i in range(len(A))]
  prob = cp.Problem(cp.Minimize(-cp.log_det(B)), constraints)
  optval = prob.solve()
  if optval==np.inf:
    raise Exception("No solution possible!")
  print(f"Optimal value: {optval}") 

  return B.value, d.value

#From: https://mcmap.net/q/1009002/-how-to-generate-a-random-sample-of-points-from-a-3-d-ellipsoid-using-python
def random_point_ellipsoid(a,b,c,x0,y0,z0):
    """Generate a random point on an ellipsoid defined by a,b,c"""
    u = np.random.rand()
    v = np.random.rand()
    theta = u * 2.0 * np.pi
    phi = np.arccos(2.0 * v - 1.0)
    sinTheta = np.sin(theta);
    cosTheta = np.cos(theta);
    sinPhi = np.sin(phi);
    cosPhi = np.cos(phi);
    rx = a * sinPhi * cosTheta;
    ry = b * sinPhi * sinTheta;
    rz = c * cosPhi;
    return rx, ry, rz

def Plot(points, B, d):
  hull = ConvexHull(points)

  fig = plt.figure()
  ax = fig.add_subplot(111, projection='3d')
  ax.scatter(points[:,0], points[:,1], points[:,2], marker=".")
  display_points = np.array([random_point_ellipsoid(1,1,1,0,0,0) for i in range(1000)])
  display_points = display_points@B+d
  ax.scatter(display_points[:,0], display_points[:,1], display_points[:,2])
  plt.show()


folder_path = r"data"
dcm_img = get_dcm_image(folder_path)
points = get_surface_points(dcm_img)

B, d = FindMaximumVolumeInscribedEllipsoid(points)

Plot(points, B, d)

ball_vol = 4/3.0*np.pi*(1.0**3)

print("DCM vol: ", get_volume_ml(dcm_img))
print("Ellipsoid Volume: ", np.linalg.det(B) * ball_vol)

This gives

DCM vol:  16.2786318359375
Ellipsoid Volume:  11.947614772444393
Herzog answered 20/5, 2020 at 19:49 Comment(16)
thanks, this seems to be what I am looking for. I am well-aware that the ellipsoid volume might not represent the maximum possible volume, but that's okay for my requirements. correct me, if I am wrong - the radii of the inner ellipsoid are the diagonal values of the B matrix? B would represent the scaling matrix in a SVD of the matrix describing the ellipsoid?Underlaid
is the inner ellipsoid volume computed as vol = det(B) and is it always smaller than the volume described by the points?Underlaid
I believe the volume is det(B^-1). The eigenvalues of B are proportional to the axis lengths.Herzog
The ellipsoid I describe will always have less volume than the convex hull of the points. If your points all lie on a surface than the volume given by the ellipsoid I describe will always underestimate that volume. However, if your points don't form a well-defined surface (due to noise) than the points don't have a well-defined volume.Herzog
I am not sure about the volume being equal to det(B^-1). After I tested the algorithm on my points I get the outer ellipsoid volume = 62 (mL), det(B^-1)= 0.00019 and actual volume of the points 16 mL. My points describe the surface of a tumor object from a DICOM img and I computed the volume by counting all the pixels (nxnynz). So I find it hard to believe that the inner ellipsoid volume is so far underestimated. What do you think?Underlaid
reading Section 8.4.2 from Boyd & Vandenberghe "Convex Optimization" it says that maximum volume inscribed ellipsoid is proportional to det(B) web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf. in which case I get det(B)=5, but then again it looks way too small again compared to the outer ellipsoid and the actual volume of the pts.Underlaid
@Roxanne: The exact volume is sqrt(det(B)) * Vol(Ball(0, 1)). That is, the square root of the determinant multiplied by the volume of the 3D unit ball (link). Were you able to plot your results to see if the ellipsoid being generated looks reasonable given your data?Herzog
thanks. I plotted the results and it looks "reasonable" on my data, however a lot of points are outside. if you have time, you can check out the actual data I am testing on (see updated Question with uploaded data) and see what I mean.Underlaid
@Roxanne: Could you post how you are calculating the volume by counting pixels? Since this is an interior ellipsoid, you should expect that there will be points outside.Herzog
see updated Question, I added function def get_volume_ml(image). However, I compute the inner ellipsoid volume using the formula you found volume_inner_ellipsoid = np.sqrt(la.det(B) / 1000) * 4.19, where 4.19 ~= Vol(Ball(0,1)Underlaid
@Roxanne: See my updated answer. You need to convert the pixel indices into coordinates appropriate for your problem.Herzog
@Roxanne: You could also think about using an L1 or L2 best-fit ellipsoid, though I haven't had much luck with this yet numerically (see here).Herzog
@Roxanne: I wonder what the larger picture is with fitting the ellipsoid to the tumor volume? Get in touch if you'd like to discuss it. My email's on my website.Herzog
To use this code we need many modules, and visual studio and blas+lapack it is a nightmareAbolition
@seghier: The packages should all install pretty easily off of pip or conda and blas+lapack should have precompiled packages in most distro repos.Herzog
Reference: <web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf>.Sophy
D
0

One, perhaps the mathematically standard, way to represent (the surface of) an ellipsoid is that it is the set

{ X | (X-a)'*inv(C)*(X-a) = 1}
the solid ellipsoid is then 
{ X | (X-a)'*inv(C)*(X-a) <= 1}

Here C is a 3x3 symmetric positive definite matrix and a is the 'centre' of the ellipsoid.

We can make this a bit easier to deal with by using the cholesky decompoisition, that is finding a lower triangular matrix L so that

C = L*L'

and using the inverse M of L (L being lower triangular this is easy to compute). We than have that the solid ellipsoid is

{ X | (M*(X-a))'*(M*(X-a)) <= 1}
= { | ||M*(X-a))|| <= 1} where ||.|| is the euclidean 

norm

We have a bunch of points X[] and an ellipsoid (C,a) containing them, that is

for all i ||M*(X[i]-a)|| <= 1
i.e. for all i ||Y[i]|| <= 1 where Y[i] = M*(X[i]-a)

Now we want to transform the ellipsoid (ie change C and a) so that all the points are outside the transformed ellipsoid. We may as well transform M and a instead.

The very simplest thing to do would just be to scale M by a constant s, and leave a alone. This amounts to scaling all the Y[] and in this case its easy to see that the scale factor to use would be one over the minimum of the ||Y[i]||. This way all the points will be outside or on the transformed ellipsoid, and some will be on it, so the transformed ellipse is as big as possible.

In terms of D,a the new ellipse is then

D = (1/(s*s))*C

If this simple approach gives acceptable results, it's what I'd use.

Without moving the centre, the most general thing to do, I think, would be to change

M to N*M

with the constraints that N is upper triangular and has positive numbers on the diagonal. We require of N that

N*Y[i] >= 1 for all i

We need a criterion for selecting N. One would be that it should reduce the volume as little as possible, that is that the determinant (which for a lower triangular matrix is just the product of the diagonal elements) should be as small as possible, subject to the constraints.

There may well be packages that can do this sort of thing, but alas I don't know which ones (which should be taken more as an indication of my ignorance than as an indication that there aren't such packages).

Once N is found, the transformed C matrix is

D = L*inv(N)*inv(N')*L'

You could also change a. The details I leave to the interested reader...

Dish answered 18/5, 2020 at 17:5 Comment(0)
B
0

I think that if you can assume that the center of mass of the ellipsoids and your points are the same, you could just solve the equation for the ellipsoid passing through the closest or the farthest n points from the center of mass. I am not sure I will have time to beef up this answer but this approach should be quite simple to implement with standard Python tools, e.g.:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.center_of_mass.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html

and perhaps SymPy to solve the analytical equation.

Blumenthal answered 25/5, 2020 at 21:48 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.