PRINCIPALAXES Principal axes of a set of ND points. [CENTER, ROTMAT] = principalAxes(PTS) [CENTER, ROTMAT, SCALES] = principalAxes(PTS) Computes the principal axes of a set of points given in a N-by-ND array and returns the result in two or three outputs: CENTER is the centroid of the points, as a 1-by-ND row vector ROTMAT represents the orientation of the point cloud, as a ND-by-ND rotation matrix SCALES is the scaling factor along each dimension, as a 1-by-ND row vector. Example pts = randn(100, 2); pts = transformPoint(pts, createScaling(5, 2)); pts = transformPoint(pts, createRotation(pi/6)); pts = transformPoint(pts, createTranslation(3, 4)); [center, rotMat] = principalAxes(pts); See also equivalentEllipse, equivalentEllipsoid, principalAxesTransform
0001 function varargout = principalAxes(points) 0002 %PRINCIPALAXES Principal axes of a set of ND points. 0003 % 0004 % [CENTER, ROTMAT] = principalAxes(PTS) 0005 % [CENTER, ROTMAT, SCALES] = principalAxes(PTS) 0006 % Computes the principal axes of a set of points given in a N-by-ND array 0007 % and returns the result in two or three outputs: 0008 % CENTER is the centroid of the points, as a 1-by-ND row vector 0009 % ROTMAT represents the orientation of the point cloud, as a ND-by-ND 0010 % rotation matrix 0011 % SCALES is the scaling factor along each dimension, as a 1-by-ND row 0012 % vector. 0013 % 0014 % Example 0015 % pts = randn(100, 2); 0016 % pts = transformPoint(pts, createScaling(5, 2)); 0017 % pts = transformPoint(pts, createRotation(pi/6)); 0018 % pts = transformPoint(pts, createTranslation(3, 4)); 0019 % [center, rotMat] = principalAxes(pts); 0020 % 0021 % See also 0022 % equivalentEllipse, equivalentEllipsoid, principalAxesTransform 0023 % 0024 0025 % ------ 0026 % Author: David Legland 0027 % e-mail: david.legland@inrae.fr 0028 % Created: 2019-08-12, using Matlab 9.6.0.1072779 (R2019a) 0029 % Copyright 2019 INRAE - Cepia Software Platform. 0030 0031 % number and dimension of points 0032 n = size(points, 1); 0033 nd = size(points, 2); 0034 0035 % compute centroid 0036 center = mean(points); 0037 0038 % compute the covariance matrix 0039 covPts = cov(points) / n; 0040 0041 % perform a principal component analysis to extract principal axes 0042 [rotMat, S] = svd(covPts); 0043 0044 % extract length of each semi axis 0045 radii = sqrt(diag(S) * n); 0046 0047 % sort axes from greater to lower 0048 [radii, ind] = sort(radii, 'descend'); 0049 radii = radii'; 0050 0051 % format U to ensure first axis points to positive x direction 0052 rotMat = rotMat(ind, :); 0053 if rotMat(1,1) < 0 && nd > 2 0054 rotMat = -rotMat; 0055 % keep matrix determinant positive 0056 rotMat(:,3) = -rotMat(:,3); 0057 end 0058 0059 % format output 0060 if nargout == 2 0061 varargout = {center, rotMat}; 0062 elseif nargout == 3 0063 varargout = {center, rotMat, radii}; 0064 end