Minimum Volume Enclosing Ellipsoid

1 http://stackoverflow.com/questions/1768197/bounding-ellipse/1768440#1768440 (stack)

2 http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.116.7691 (document)

3 https://gist.github.com/Gabriel-p/4ddd31422a88e7cdf953 (python mvee code)

 

 

$ {(x_i -c)}^TE{(x_i -c)} \leq 1 \hspace{10mm} i=1,...,m $

$Vol(\varepsilon) = {v_0 \over \sqrt{det(E)}} = v_0 det(E^{-1})^{1 \over 2}$

 

document 를 참조하면 center form 을 이용하여 convex function 으로 변환 후 이를 이용하여

Largrange dual problem 을 가정하여 최적화를 진행하는 방식입니다. 세부 수식은 2 document 를 참조하셔야 됩니다.

 

 

import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

import cv2

def mvee(points, tol = 0.01):
    """
    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

img = cv2.imread('/home/roadcom/Downloads/partial_ellipse.jpeg', cv2.IMREAD_GRAYSCALE)

ret, th = cv2.threshold(img,0,255,cv2.THRESH_BINARY + cv2.THRESH_OTSU)

conts, _ = cv2.findContours(th,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_NONE)

# max bloc filter
max_blob_id = np.argmax([cv2.contourArea(p) for p in conts])

A, centroid = mvee(np.squeeze(conts[max_blob_id]))
U, D, V = la.svd(A)

# x, y radii.
rx, ry = 1./np.sqrt(D)
major, minor = max(rx,ry), min(rx,ry)
    
cxy = tuple(np.int0(centroid))
aixs = tuple(np.int0([major, minor]))

arcsin = -1. * np.rad2deg(np.arcsin(V[0][0]))
arccos = np.rad2deg(np.arccos(V[0][1]))

# Orientation angle (with respect to the x axis counterclockwise).
alpha = arccos if arcsin > 0. else -1. * arccos

rst_img = cv2.ellipse(img.copy(), cxy, aixs, alpha, 0, 360,(255,255,255),3)
B, G, R = cv2.split(rst_img) 
matplot_img = cv2.merge([R,G,B])

fig, axs = plt.subplots(1,2, figsize=(16,8))
axs[0].imshow(img, interpolation='quadric')
axs[0].set_title('Origin[B,G,R]')
axs[0].axis('off')

axs[1].imshow(matplot_img, interpolation='quadric')
axs[1].set_title('MVEE')
axs[1].axis('off')