oriented bounding box for points

oobb.py

#!/usr/bin/env python3

# https://stackoverflow.com/questions/32892932/create-the-oriented-bounding-box-obb-with-python-and-numpy
def oobb(pts):
    '''
    given: 
    	pts, 2d or 3d
    return:
    	corners, center of the oobb
    '''
    # covariance matrix 
    cot_mat = np.cov(pts, rowvar = False, bias = 1) 
    # eigen vector
    v, vect = np.linalg.eig( cot_mat )
    tvect = np.transpose(vect)
    #use the inverse of the eigenvectors as a rotation matrix and
    #rotate the points so they align with the x and y axes
    rpts = np.dot(pts, np.linalg.inv(tvect))  
    # min and max aligned with the rotated axes
    minb = np.min(rpts, axis = 0)
    maxb = np.max(rpts, axis = 0)
    diff = (maxb - minb) * 0.5
    # the center is just half way between the min and max xyz
    center = minb + diff
    # 2d points
    if points.shape[-1] == 2:
        x, y = diff
        corners = np.array([
        center+[-x,-y],
        center+[ x,-y],
        center+[ x, y],
        center+[-x, y],
                ])
    # 3d points
    else:
        x, y, z = diff
        corners = np.array([
        center+[-x,-y,-z],
        center+[ x,-y,-z],
        center+[ x, y,-z],
        center+[-x, y,-z],
        center+[-x,-y, z],
        center+[ x,-y, z],
        center+[ x, y, z],
        center+[-x, y, z],
                ])
    #use the the eigenvectors as a rotation matrix and
    #rotate the corners and the centerback
    corners = np.dot(corners,tvect)
    center = np.dot(center,tvect)
    return corners, center


# to test the function in 2d
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
n = 20
theta = np.linspace(0, 2*np.pi, n)
points = np.hstack( (  ( 0.25*np.cos(theta) +  0.02 * np.random.rand(n) ).reshape(-1,1),
                       ( 0.25*np.sin(theta) +  0.02 * np.random.rand(n) ).reshape(-1,1)))
points = points + np.array((0.5,0.5))

corners, center = oobb(points)

plt.plot(points[:,0], points[:,1])
plt.scatter(corners[:,0], corners[:,1])
plt.scatter(center[0], center[1])
plt.axis('scaled')
plt.show()