Little script to compute the minimal oriented bounding box of a given mesh. Takes inspiration from "Fast oriented bounding box optimization on the rotation group SO(3, R)" from Chang et al. 2011, and adapts it to the python ecosystem.

oriented_bounding_box.py

#!/usr/bin/env python3

"""
Little script to compute the minimal oriented bounding box (OBB)
of a given mesh.

Takes inspiration from 
"Fast oriented bounding box optimization on the rotation group SO(3, R)"
from Chang et al. 2011, and adapts it to the python ecosystem.
"""

import argparse
import time

import numpy as np
import trimesh
from scipy.optimize import minimize, differential_evolution
from scipy.spatial.transform import Rotation


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("model", help="Path to 3D-model")
    parser.add_argument(
        "--rotate",
        action="store_true",
        help="Whether to rotate the model such, that its Axis-Aligned BB becomes minimal. Overwrites 'model' file.",
    )
    args = parser.parse_args()

    mesh = trimesh.load(args.model, force="mesh")
    volume_initial = mesh.bounding_box.volume
    volume_obb_trimesh = mesh.bounding_box_oriented.volume

    T_aabb2obb = np.eye(4)
    T_aabb2obb[:3, :3] = compute_obb(mesh)

    mesh.apply_transform(T_aabb2obb)
    volume_obb = mesh.bounding_box.volume

    print(
        f"BBox Volume: initial={volume_initial:.3f} obbtrimesh={volume_obb_trimesh:.3f} obbthis={volume_obb:.3f}"
    )
    if args.rotate:
        mesh.export(args.model)


def compute_obb(mesh):
    dims_orig = np.max(mesh.convex_hull.vertices, axis=0) - np.min(
        mesh.convex_hull.vertices, axis=0
    )
    vol_orig = np.product(dims_orig)

    def f(q_wxyz):
        q_wxyz /= np.linalg.norm(q_wxyz)
        verts_new = Rotation.from_quat(q_wxyz).apply(mesh.convex_hull.vertices)

        E_vol = (
            np.product(np.max(verts_new, axis=0) - np.min(verts_new, axis=0)) / vol_orig
        )
        # try to rotate as little as possible, by including rotation magnitude
        E_rot = Rotation.from_quat(q_wxyz).magnitude() / np.pi

        return (0.5 * (E_vol**2 + E_rot**2)) ** (1 / 2)

    # The default polish uses L-BFGS-B, which in my tests was worse than polishing with
    # Nelder-Mead, 60% of the time. Using L-BFGS-B after Nelder-Mead improved results
    # further, 80% of the time. However, here we're entering diminishing returns territory
    t0 = time.perf_counter()
    rde = differential_evolution(f, [(-1, 1), (-1, 1), (-1, 1), (-1, 1)], polish=False)
    t1 = time.perf_counter()
    res_nm = minimize(f, rde.x, method="Nelder-Mead")
    t2 = time.perf_counter()
    res_comb = minimize(f, res_nm.x, method="L-BFGS-B")
    t3 = time.perf_counter()
    final_wxyz = res_comb.x

    dur_de_ms = (t1 - t0) * 1000
    dur_nm_ms = (t2 - t1) * 1000
    dur_bf_ms = (t3 - t2) * 1000
    print(
        f"Timings: genetic={dur_de_ms:.0f}ms nelderm={dur_nm_ms:.0f}ms bfgs={dur_bf_ms:.0f}ms"
    )

    R_aabb2obb = Rotation.from_quat(final_wxyz).as_matrix()
    return R_aabb2obb


if __name__ == "__main__":
    main()

$ for file in *.ply ; do echo $file ; ./obb.py $file ; done
suzanne2012.ply
Timings: genetic=202ms nelderm=7ms bfgs=16ms
BBox Volume: initial=8.106 obbtrimesh=6.331 obbthis=6.277

suzanne507.ply
Timings: genetic=103ms nelderm=7ms bfgs=12ms
BBox Volume: initial=9.168 obbtrimesh=7.309 obbthis=7.162

suzanne7958.ply
Timings: genetic=196ms nelderm=13ms bfgs=50ms
BBox Volume: initial=8.126 obbtrimesh=6.317 obbthis=6.305

tricky_cube.ply
Timings: genetic=151ms nelderm=5ms bfgs=6ms
BBox Volume: initial=0.016 obbtrimesh=0.016 obbthis=0.016

weird_thing.ply
Timings: genetic=150ms nelderm=7ms bfgs=14ms
BBox Volume: initial=277.871 obbtrimesh=62.007 obbthis=53.626