Equivalent ellipsoid of a set of 3D points. ELL = equivalentEllipsoid(PTS) Compute the equivalent ellipsoid of the set of points PTS. The result is an ellipsoid defined by: ELL = [XC YC ZC A B C PHI THETA PSI] where [XC YC ZY] is the center, [A B C] are the lengths of the semi-axes (in decreasing order), and [PHI THETA PSI] are Euler angles representing the ellipsoid orientation, in degrees. Example pts = randn(300, 3); pts = transformPoint3d(pts, createScaling3d([6 4 2])); pts = transformPoint3d(pts, createRotationOx(pi/6)); pts = transformPoint3d(pts, createRotationOy(pi/4)); pts = transformPoint3d(pts, createRotationOz(pi/3)); pts = transformPoint3d(pts, createTranslation3d([5 4 3])); elli = equivalentEllipsoid(pts); figure; drawPoint3d(pts); axis equal; hold on; drawEllipsoid(elli, ... 'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3); See also spheres, drawEllipsoid, equivalentEllipse, principalAxes principalAxesTransform, rotation3dToEulerAngles
0001 function ell = equivalentEllipsoid(points) 0002 % Equivalent ellipsoid of a set of 3D points. 0003 % 0004 % ELL = equivalentEllipsoid(PTS) 0005 % Compute the equivalent ellipsoid of the set of points PTS. The result 0006 % is an ellipsoid defined by: 0007 % ELL = [XC YC ZC A B C PHI THETA PSI] 0008 % where [XC YC ZY] is the center, [A B C] are the lengths of the 0009 % semi-axes (in decreasing order), and [PHI THETA PSI] are Euler angles 0010 % representing the ellipsoid orientation, in degrees. 0011 % 0012 % Example 0013 % pts = randn(300, 3); 0014 % pts = transformPoint3d(pts, createScaling3d([6 4 2])); 0015 % pts = transformPoint3d(pts, createRotationOx(pi/6)); 0016 % pts = transformPoint3d(pts, createRotationOy(pi/4)); 0017 % pts = transformPoint3d(pts, createRotationOz(pi/3)); 0018 % pts = transformPoint3d(pts, createTranslation3d([5 4 3])); 0019 % elli = equivalentEllipsoid(pts); 0020 % figure; drawPoint3d(pts); axis equal; 0021 % hold on; drawEllipsoid(elli, ... 0022 % 'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3); 0023 % 0024 % See also 0025 % spheres, drawEllipsoid, equivalentEllipse, principalAxes 0026 % principalAxesTransform, rotation3dToEulerAngles 0027 0028 % ------ 0029 % Author: David Legland 0030 % e-mail: david.legland@inrae.fr 0031 % Created: 2011-03-12, using Matlab 7.9.0.529 (R2009b) 0032 % Copyright 2011 INRA - Cepia Software Platform. 0033 0034 % number of points 0035 n = size(points, 1); 0036 0037 % compute centroid 0038 center = mean(points); 0039 0040 % compute the covariance matrix 0041 covPts = cov(points)/n; 0042 0043 % perform a principal component analysis with 3 variables, 0044 % to extract equivalent axes 0045 [U, S] = svd(covPts); 0046 0047 % extract length of each semi axis 0048 radii = sqrt(5) * sqrt(diag(S)*n)'; 0049 0050 % sort axes from greater to lower 0051 [radii, ind] = sort(radii, 'descend'); 0052 0053 % format U to ensure first axis points to positive x direction 0054 U = U(ind, :); 0055 if U(1,1) < 0 0056 U = -U; 0057 % keep matrix determinant positive 0058 U(:,3) = -U(:,3); 0059 end 0060 0061 % convert axes rotation matrix to Euler angles 0062 angles = rotation3dToEulerAngles(U); 0063 0064 % concatenate result to form an ellipsoid object 0065 ell = [center, radii, angles];