INPUTS: xyz: an N*3 array of N points contained in cluster Vol : Volume of ellipsoid OUTPUTS: U: eigenvector (principal direction) r :radius of principal axis for Vol Ncor : Number of points inside the ellipsoid Copyright (C) 2011, ATR All Rights Reserved. License : New BSD License(see VBMEG_LICENSE.txt)
0001 function [Ncor,r,U,S] = vb_ellipsoid_check(xyz, Vol, R) 0002 % 0003 %INPUTS: 0004 % xyz: an N*3 array of N points contained in cluster 0005 % Vol : Volume of ellipsoid 0006 % 0007 %OUTPUTS: 0008 % U: eigenvector (principal direction) 0009 % r :radius of principal axis for Vol 0010 % Ncor : Number of points inside the ellipsoid 0011 % 0012 % Copyright (C) 2011, ATR All Rights Reserved. 0013 % License : New BSD License(see VBMEG_LICENSE.txt) 0014 0015 N = size(xyz,1); 0016 X0 = sum(xyz)/N; 0017 xyz = xyz - repmat(X0, [N 1]); 0018 0019 % covariance matrix 0020 XX = xyz'*xyz; 0021 0022 % eigenvector of covariance matrix 0023 [U, S] = eig(XX); 0024 0025 % eigenvalue of covariance matrix 0026 S = sqrt(diag(S)); 0027 0028 if exist('R','var'), 0029 R = sort(R); 0030 [S, ix] = sort(S); 0031 U = U(:,ix); 0032 S = R; 0033 end; 0034 0035 0036 % Vol = (4*pi/3)*R^3 = (4*pi/3)*r(1)*r(2)*r(3) 0037 R3 = Vol*(3/pi*4); 0038 0039 % radius of principal axis for Vol 0040 r = S*(R3/prod(S))^(1/3); 0041 0042 % transform to principal coordinate 0043 y = xyz*U; 0044 0045 % normalized radius 0046 dd = (y(:,1)/r(1)).^2 + (y(:,2)/r(2)).^2 + (y(:,3)/r(3)).^2 ; 0047 0048 % find inside of ellipsoid 0049 ix = find( dd <= 1 ); 0050 0051 % Number of points inside the ellipsoid 0052 Ncor = length(ix); 0053 0054