orientedBox3d

ORIENTEDBOX3D Object-oriented bounding box of a set of 3D points.

   OOBB = orientedBox3d(PTS)
   Returns the oriented bounding box of the collection of points in the
   N-by-3 array PTS. The result is given as:
   [XC YC ZC  L W H  PHI THETA PSI]
   where (XC,YC,ZC) corresponds to the center of the box, (L,W,H)
   corresponds to the length, width, and depth of the box, and (PHI,
   THETA, PSI) is the orientation of the box as Euler angles.

   [OOBB, ROT] = orientedBox3d(PTS)
   Also returns the rotation matrix of the point cloud, as a 3-by-3
   numeric array.

   Example
     [v, f] = sphereMesh;
     phi=-360+720*rand; theta=-360+720*rand; psi=-360+720*rand;
     angles = [phi, theta, psi];
     rotMat = eulerAnglesToRotation3d(angles);
     rotMat(1:3,4) = randi([-100,100],3,1);
     scale = [randi([7,9],1,1), randi([4,6],1,1), randi([1,3],1,1)];
     pts = transformPoint3d(bsxfun(@times, v, scale), rotMat);
     box3d = orientedBox3d(pts);
     figure; drawPoint3d(pts, '.'); 
     axis equal; xlabel('x'); ylabel('y'); zlabel('z');
     drawCuboid(box3d, 'FaceColor', 'none');

   See also
     meshes3d, drawCuboid, rotation3dToEulerAngles

0001 function [box3d, rotMat] = orientedBox3d(pts)
0002 %ORIENTEDBOX3D Object-oriented bounding box of a set of 3D points.
0003 %
0004 %   OOBB = orientedBox3d(PTS)
0005 %   Returns the oriented bounding box of the collection of points in the
0006 %   N-by-3 array PTS. The result is given as:
0007 %   [XC YC ZC  L W H  PHI THETA PSI]
0008 %   where (XC,YC,ZC) corresponds to the center of the box, (L,W,H)
0009 %   corresponds to the length, width, and depth of the box, and (PHI,
0010 %   THETA, PSI) is the orientation of the box as Euler angles.
0011 %
0012 %   [OOBB, ROT] = orientedBox3d(PTS)
0013 %   Also returns the rotation matrix of the point cloud, as a 3-by-3
0014 %   numeric array.
0015 %
0016 %   Example
0017 %     [v, f] = sphereMesh;
0018 %     phi=-360+720*rand; theta=-360+720*rand; psi=-360+720*rand;
0019 %     angles = [phi, theta, psi];
0020 %     rotMat = eulerAnglesToRotation3d(angles);
0021 %     rotMat(1:3,4) = randi([-100,100],3,1);
0022 %     scale = [randi([7,9],1,1), randi([4,6],1,1), randi([1,3],1,1)];
0023 %     pts = transformPoint3d(bsxfun(@times, v, scale), rotMat);
0024 %     box3d = orientedBox3d(pts);
0025 %     figure; drawPoint3d(pts, '.');
0026 %     axis equal; xlabel('x'); ylabel('y'); zlabel('z');
0027 %     drawCuboid(box3d, 'FaceColor', 'none');
0028 %
0029 %   See also
0030 %     meshes3d, drawCuboid, rotation3dToEulerAngles
0031 
0032 % ------
0033 % Author: David Legland
0034 % E-mail: david.legland@inrae.fr
0035 % Created: 2015-12-01, using Matlab 8.6.0.267246 (R2015b)
0036 % Copyright 2015-2024 INRA - Cepia Software Platform
0037 
0038 tri = convhulln(pts);
0039 nFaces = size(tri, 1);
0040 
0041 %% identify index of face with smallest width
0042 indMinBreadth = 0;
0043 minBreadth = Inf;
0044 for iFace = 1:nFaces
0045     faceInds = tri(iFace, :);
0046     plane = createPlane(pts(faceInds, :));
0047     
0048     breadth = max(abs(distancePointPlane(pts, plane)));
0049     
0050     if breadth < minBreadth
0051         minBreadth = breadth;
0052         indMinBreadth = iFace;
0053     end
0054 end
0055 
0056 % compute projection on reference plane
0057 refPlane = createPlane(pts(tri(indMinBreadth, :), :));
0058 pts2d = planePosition(projPointOnPlane(pts, refPlane), refPlane);
0059 
0060 % compute 2D OOBB for projected points
0061 box2d = orientedBox(pts2d);
0062 
0063 % extract reference points from planar OOBB: the center, and two direction
0064 % vectors
0065 center2d = box2d(1:2);
0066 L1 = box2d(3);
0067 L2 = box2d(4);
0068 markers2d = [0 0; L1/2 0; 0 L2/2];
0069 
0070 % orient reference points to 2d basis
0071 theta2d = box2d(5);
0072 rot = createRotation(deg2rad(theta2d));
0073 tra = createTranslation(center2d);
0074 transfo = tra * rot;
0075 markers2d = transformPoint(markers2d, transfo);
0076 
0077 % backprojection to 3D space
0078 markers3d = planePoint(refPlane, markers2d);
0079 
0080 % compute 3D vectors and center
0081 centerProj = markers3d(1,:);
0082 v1n = normalizeVector3d(markers3d(2,:) - centerProj);
0083 v2n = normalizeVector3d(markers3d(3,:) - centerProj);
0084 
0085 % compute rotation matrix and convert to Euler Angles
0086 v3n = crossProduct3d(v1n, v2n);
0087 rotMat = [v1n' v2n' v3n'];
0088 boxAngles = rotation3dToEulerAngles(rotMat);
0089 
0090 % create 3D object-oriented bounding box
0091 boxCenter3d = centerProj + v3n * minBreadth/2;
0092 box3d = [boxCenter3d L1 L2 minBreadth boxAngles];