fromPoints method

OBB fromPoints(List<Vector3> points)

Computes the minimum enclosing OBB for the given set of points. The method is an implementation of {@link http://gamma.cs.unc.edu/users/gottschalk/main.pdf Collision Queries using Oriented Bounding Boxes} by Stefan Gottschalk. According to the dissertation, the quality of the fitting process varies from the respective input. This method uses the best approach by computing the covariance matrix based on the triangles of the convex hull (chapter 3.4.3).

However, the implementation is susceptible to {@link https://en.wikipedia.org/wiki/Regular_polygon regular polygons} like cubes or spheres. For such shapes, it's recommended to verify the quality of the produced OBB. Consider to use an AABB or bounding sphere if the result is not satisfying.

Implementation

/// of the produced OBB. Consider to use an AABB or bounding sphere if the result
/// is not satisfying.
OBB fromPoints(List<Vector3> points ) {
	final convexHull = ConvexHull().fromPoints( points );

	// 1. iterate over all faces of the convex hull and triangulate
	final faces = convexHull.faces;
	final List<HalfEdge> edges = [];
	final List<double> triangles = [];

	for ( int i = 0, il = faces.length; i < il; i ++ ) {
		final face = faces[ i ];
		HalfEdge? edge = face.edge;

		edges.length = 0;

		// gather edges
		do{
			edges.add( edge! );
			edge = edge.next;
		} while ( edge != face.edge );

		// triangulate
		final triangleCount = ( edges.length - 2 );

		for ( int j = 1, jl = triangleCount; j <= jl; j ++ ) {
			final v1 = edges[ 0 ].vertex;
			final v2 = edges[ j + 0 ].vertex;
			final v3 = edges[ j + 1 ].vertex;

			triangles.addAll([ v1.x, v1.y, v1.z ]);
			triangles.addAll( [v2.x, v2.y, v2.z] );
			triangles.addAll( [v3.x, v3.y, v3.z] );
		}
	}

	// 2. build covariance matrix

	final p = Vector3();
	final q = Vector3();
	final r = Vector3();

	final qp = Vector3();
	final rp = Vector3();

	final v = Vector3();

	final mean = Vector3();
	final weightedMean = Vector3();
	double areaSum = 0;

	double cxx, cxy, cxz, cyy, cyz, czz;
	cxx = cxy = cxz = cyy = cyz = czz = 0;

	for ( int i = 0, l = triangles.length; i < l; i += 9 ) {

		p.fromArray( triangles, i );
		q.fromArray( triangles, i + 3 );
		r.fromArray( triangles, i + 6 );

		mean.set( 0, 0, 0 );
		mean.add( p ).add( q ).add( r ).divideScalar( 3 );

		qp.subVectors( q, p );
		rp.subVectors( r, p );

		final area = v.crossVectors( qp, rp ).length / 2; // .length() represents the frobenius norm here
		weightedMean.add( v.copy( mean ).multiplyScalar( area ) );

		areaSum += area;

		cxx += ( 9.0 * mean.x * mean.x + p.x * p.x + q.x * q.x + r.x * r.x ) * ( area / 12 );
		cxy += ( 9.0 * mean.x * mean.y + p.x * p.y + q.x * q.y + r.x * r.y ) * ( area / 12 );
		cxz += ( 9.0 * mean.x * mean.z + p.x * p.z + q.x * q.z + r.x * r.z ) * ( area / 12 );
		cyy += ( 9.0 * mean.y * mean.y + p.y * p.y + q.y * q.y + r.y * r.y ) * ( area / 12 );
		cyz += ( 9.0 * mean.y * mean.z + p.y * p.z + q.y * q.z + r.y * r.z ) * ( area / 12 );
		czz += ( 9.0 * mean.z * mean.z + p.z * p.z + q.z * q.z + r.z * r.z ) * ( area / 12 );
	}

	weightedMean.divideScalar( areaSum );

	cxx /= areaSum;
	cxy /= areaSum;
	cxz /= areaSum;
	cyy /= areaSum;
	cyz /= areaSum;
	czz /= areaSum;

	cxx -= weightedMean.x * weightedMean.x;
	cxy -= weightedMean.x * weightedMean.y;
	cxz -= weightedMean.x * weightedMean.z;
	cyy -= weightedMean.y * weightedMean.y;
	cyz -= weightedMean.y * weightedMean.z;
	czz -= weightedMean.z * weightedMean.z;

	final covarianceMatrix = Matrix3();

	covarianceMatrix.elements[ 0 ] = cxx;
	covarianceMatrix.elements[ 1 ] = cxy;
	covarianceMatrix.elements[ 2 ] = cxz;
	covarianceMatrix.elements[ 3 ] = cxy;
	covarianceMatrix.elements[ 4 ] = cyy;
	covarianceMatrix.elements[ 5 ] = cyz;
	covarianceMatrix.elements[ 6 ] = cxz;
	covarianceMatrix.elements[ 7 ] = cyz;
	covarianceMatrix.elements[ 8 ] = czz;

	// 3. compute rotation, center and half sizes

	covarianceMatrix.eigenDecomposition( eigenDecomposition );

	final unitary = eigenDecomposition['unitary']!;

	final v1 = Vector3();
	final v2 = Vector3();
	final v3 = Vector3();

	unitary.extractBasis( v1, v2, v3 );

	double u1 = - double.infinity;
	double u2 = - double.infinity;
	double u3 = - double.infinity;
	double l1 = double.infinity;
	double l2 = double.infinity;
	double l3 = double.infinity;

	for ( int i = 0, l = points.length; i < l; i ++ ) {
		final p = points[ i ];

		u1 = math.max( v1.dot( p ), u1 );
		u2 = math.max( v2.dot( p ), u2 );
		u3 = math.max( v3.dot( p ), u3 );

		l1 = math.min( v1.dot( p ), l1 );
		l2 = math.min( v2.dot( p ), l2 );
		l3 = math.min( v3.dot( p ), l3 );
	}

	v1.multiplyScalar( 0.5 * ( l1 + u1 ) );
	v2.multiplyScalar( 0.5 * ( l2 + u2 ) );
	v3.multiplyScalar( 0.5 * ( l3 + u3 ) );

	// center
	center.add( v1 ).add( v2 ).add( v3 );

	halfSizes.x = u1 - l1;
	halfSizes.y = u2 - l2;
	halfSizes.z = u3 - l3;

	// halfSizes
	halfSizes.multiplyScalar( 0.5 );

	// rotation
	rotation.copy( unitary );

	return this;
}