基于Matlab实现空间数据的主成分局部均值聚类.zip

% function YY = PCLM_estimate( X, d,p,b, n,c,K, P ) % % usage: YY = PCLM_estimate( X, 3,1,2, 30,1,6, 1 ); % PCLM = Principal Components Local Mean spatial clustering % % Performs Principal Curve estimation of 2D,3D point clouds, on the basis % of Local Means (LM) and principal components (PC) of local Cov Matrices. % It is a "blurred method" and needs only a FEW iterations to cluster; % however, standard criteria of convergence are not applicable. % % INPUT: % X = Data matrix Nx4: [x, y, z, w], % (Longitude, Latitude, Depth, Magnitude), % (Add columns as [X, ones(N,1)], if necessary) % % d = 2, 3, Space Dimension of the estimates (2D,3D) ... d=2, % p <= d-1, dimension of the Projection ... p=1, % (when d=3, p=1 estimates Surfaces /sheets, % p=2 estimates Curves /snakes ), % b <= d, dimension of the Neighbors' search ... b=d, % >= 2, (when d=3, b=2 provides vertical ordering), % and 3d estimates look like sections % n = 10, 20 ... number of Nearest neighbors searched ... n=30, % c = 0, 1, 2: global, local, mixed Covariance matrix ... c=1, % (use c=0,2 when the curves have the same direction) % K = maximum number of iterations ... K=10, % (choose a mild value as the method is blurring) % % P = 1, 0, yes/no plot of the results, P=1, % % OUTPUT: % YY = a structure array with K estimates Y (Nxd) % % USAGE: % X = [X, ones(length(X),2)]; % if size(X,2)<3; % YY = PCLM_estimate( X, 2,1,2, 30,1,6, 1 ); % YY = PCLM_estimate( X, 3,2,3, 30,1,6, 1 ); % YY = PCLM_estimate( X, 3,1,2, 30,1,6, 1 ); % % ----------------------------------------------------------- % % Author: Carlo Grillenzoni (c) (carlog@iuav.it) 2015 % IUAV: Institute of Architecture, University of Venice % % Reference: Smoothing Three-Dimensional Manifold Data, % with Application to Tectonic Fault Detection, % Mathematical Geosciences, 48, 2016 % https://link.springer.com/article/10.1007/s11004-015-9630-x % N = size(X,1); Xx = X(:,1:d); Y=Xx; Yy=Y*0; w = X(:,end); m0 = (w'*Xx)/sum(w); S0=(repmat(w,1,d).*Xx)'*(Xx/sum(w))-(m0'*m0); for k=1:K if P==1, k end [I, D]= knn_search( Y(:,1:b), Y(:,1:b), n); Yo=Y; for i=1:N yi = Y(i,:); Yi = [ yi; Y(I(i,:),:) ]; wi = [ w(i); w(I(i,:)) ]; mi = (wi'*Yi)/sum(wi); % Si =(repmat(wi,1,d).*Yi)'*(Yi/sum(wi))-(mi'*mi); if c==0, S=S0; elseif c==1, S=Si; else, S=(S0+Si)/2; end % [V, L] = eig(S); [ls, is] = sort(diag(L)); % ,'ascend'); Vs = V(:,is); v=Vs(:,1:p); u=Vs(:,p+1:d); % v=V(:,1:p); u=V(:,p+1:d); % no sort Yy(i,:) = yi*(u*u') + mi*(v*v') ; end % i Y = Yy; YY{k} = Y; end % k if P==1, if d==2, if K>3, figure for k=1:4 h=k; % round(k*(K/4)); subplot(2,2,k); plot(X(:,1),X(:,2),'.b') hold on; Yk=YY{h}; plot(Yk(:,1),Yk(:,2),'.r') axis tight; title(sprintf('%d-th iteration',h)) end end % K figure; plot(X(:,1),X(:,2),'xb') hold on; plot(Y(:,1),Y(:,2),'.r') axis tight; title(sprintf('%d-th iteration',K)) end % d if d==3, figure subplot(221); plot(X(:,1),X(:,2),'.b') % 'xb','MarkerSize',3) hold on; plot(Y(:,1),Y(:,2),'.r'); axis tight; ax=axis; title('2D view') subplot(222); plot(Y(:,1),Y(:,2),'.k') axis(ax); grid on; title('Estimate') subplot(223); plot3(Y(:,1),Y(:,2),Y(:,3),'.r') hold on; plot3(X(:,1),X(:,2),X(:,3),'.b') box on; grid on; axis tight; ; view(-40,60) title('3D view') subplot(224); plot3(Y(:,1),Y(:,2),Y(:,3),'.k') box on; grid on; axis tight; view(-40,60) title('Estimate') figure; plot3(Y(:,1),Y(:,2),Y(:,3),'.r') hold on; plot3(X(:,1),X(:,2),X(:,3),'.b') box on; grid on; axis tight; ; view(-40,60) title('3D view') end % d end % P % pro memoria % PCLM_estimate( [X, ones(1300,1)], 2,1,2, 35,1,9, 1 ); % simulated % PCLM_estimate( XX(1:1500,:), 2,1,2, 20,2,9, 1 ); % california % PCLM_estimate( Xx(1:1500,:), 2,1,2, 40,1,9, 1 ); % minnesota