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;
     rotMat = eulerAnglesToRotation3d(30, 20, 10);
     pts = transformPoint3d(bsxfun(@times, v, [5 3 1]), rotMat);
     box3d = orientedBox3d(pts);
     figure; drawPoint3d(pts, '.'); hold on;
     axis equal; axis([-6 6 -6 6 -5 5]);
     h = drawCuboid(box3d);
     set(h, 'facecolor', 'none');

   See also
     meshes3d, drawCuboid, rotation3dToEulerAngles

0001 function [box3d, rotMat] = orientedBox3d(pts)
0002 % 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 %     rotMat = eulerAnglesToRotation3d(30, 20, 10);
0019 %     pts = transformPoint3d(bsxfun(@times, v, [5 3 1]), rotMat);
0020 %     box3d = orientedBox3d(pts);
0021 %     figure; drawPoint3d(pts, '.'); hold on;
0022 %     axis equal; axis([-6 6 -6 6 -5 5]);
0023 %     h = drawCuboid(box3d);
0024 %     set(h, 'facecolor', 'none');
0025 %
0026 %   See also
0027 %     meshes3d, drawCuboid, rotation3dToEulerAngles
0028  
0029 % ------
0030 % Author: David Legland
0031 % e-mail: david.legland@inrae.fr
0032 % Created: 2015-12-01,    using Matlab 8.6.0.267246 (R2015b)
0033 % Copyright 2015 INRA - Cepia Software Platform.
0034 
0035 tri = convhulln(pts);
0036 nFaces = size(tri, 1);
0037 
0038 %% identify index of face with smallest width
0039 indMinBreadth = 0;
0040 minBreadth = Inf;
0041 for iFace = 1:nFaces
0042     faceInds = tri(iFace, :);
0043     plane = createPlane(pts(faceInds, :));
0044     
0045     breadth = max(abs(distancePointPlane(pts, plane)));
0046     
0047     if breadth < minBreadth
0048         minBreadth = breadth;
0049         indMinBreadth = iFace;
0050     end
0051 end
0052 
0053 % compute projection on reference plane
0054 refPlane = createPlane(pts(tri(indMinBreadth, :), :));
0055 pts2d = planePosition(projPointOnPlane(pts, refPlane), refPlane);
0056 
0057 % compute 2D OOBB for projected points
0058 box2d = orientedBox(pts2d);
0059 
0060 % extract reference points from planar OOBB: the center, and two direction
0061 % vectors
0062 center2d = box2d(1:2);
0063 L1 = box2d(3);
0064 L2 = box2d(4);
0065 markers2d = [0 0; L1/2 0; 0 L2/2];
0066 
0067 % orient reference points to 2d basis
0068 theta2d = box2d(5);
0069 rot = createRotation(deg2rad(theta2d));
0070 tra = createTranslation(center2d);
0071 transfo = tra * rot;
0072 markers2d = transformPoint(markers2d, transfo);
0073 
0074 % backprojection to 3D space
0075 markers3d = planePoint(refPlane, markers2d);
0076 
0077 % compute 3D vectors and center
0078 centerProj = markers3d(1,:);
0079 v1n = normalizeVector3d(markers3d(2,:) - centerProj);
0080 v2n = normalizeVector3d(markers3d(3,:) - centerProj);
0081 
0082 % compute rotation matrix and convert to Euler Angles
0083 v3n = crossProduct3d(v1n, v2n);
0084 rotMat = [v1n' v2n' v3n'];
0085 boxAngles = rotation3dToEulerAngles(rotMat);
0086 
0087 % create 3D object-oriented bounding box
0088 boxCenter3d = centerProj + v3n * minBreadth/2;
0089 box3d = [boxCenter3d L1 L2 minBreadth boxAngles];