V_MINSPANE calculate minimum spanning tree using euclidean distance [p,s]=X
Inputs: x(n,d) d-dimensional data points, one per row
Outputs: p(n-1,1) indices of the parent of each node within the tree
(data point n is the root)
s(n-1,1) list of the edges in ascending order of euclidean
distance. Thus the shortest edge goes from node
s(1) to node p(s(1))
The minimum spanning tree (or shortest spanning tree ) defines a set of
n-1 links that interconnect n points (or nodes) with the minimum total
length. We represent these links in the form of a tree with node n as
the root. Each node (except node n) has a unique parent node that is
given in the output vector p; it is possible for several nodes to share
the same parent.
It can be useful to know which of the links are the longest, so the output
argument lists them in ascending order. s could be calculated directly by
calculaing the length, l, of each link as follows:
l=sqrt(sum((x(1:n-1,:) - x(p,:)).^2),2);
[v,s]=sort(l);0001 function [p,s]=v_minspane(x) 0002 %V_MINSPANE calculate minimum spanning tree using euclidean distance [p,s]=X 0003 % 0004 % Inputs: x(n,d) d-dimensional data points, one per row 0005 % 0006 % Outputs: p(n-1,1) indices of the parent of each node within the tree 0007 % (data point n is the root) 0008 % s(n-1,1) list of the edges in ascending order of euclidean 0009 % distance. Thus the shortest edge goes from node 0010 % s(1) to node p(s(1)) 0011 % 0012 % The minimum spanning tree (or shortest spanning tree ) defines a set of 0013 % n-1 links that interconnect n points (or nodes) with the minimum total 0014 % length. We represent these links in the form of a tree with node n as 0015 % the root. Each node (except node n) has a unique parent node that is 0016 % given in the output vector p; it is possible for several nodes to share 0017 % the same parent. 0018 % It can be useful to know which of the links are the longest, so the output 0019 % argument lists them in ascending order. s could be calculated directly by 0020 % calculaing the length, l, of each link as follows: 0021 % l=sqrt(sum((x(1:n-1,:) - x(p,:)).^2),2); 0022 % [v,s]=sort(l); 0023 0024 % Copyright (C) Mike Brookes 2000-2009 0025 % Version: $Id: v_minspane.m 10865 2018-09-21 17:22:45Z dmb $ 0026 % 0027 % VOICEBOX is a MATLAB toolbox for speech processing. 0028 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0029 % 0030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0031 % This program is free software; you can redistribute it and/or modify 0032 % it under the terms of the GNU General Public License as published by 0033 % the Free Software Foundation; either version 2 of the License, or 0034 % (at your option) any later version. 0035 % 0036 % This program is distributed in the hope that it will be useful, 0037 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0038 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0039 % GNU General Public License for more details. 0040 % 0041 % You can obtain a copy of the GNU General Public License from 0042 % http://www.gnu.org/copyleft/gpl.html or by writing to 0043 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0045 0046 [np,nd]=size(x); 0047 0048 % first do delauny tessalation to find feasible edges 0049 0050 t=delaunayn(x); % nd+1 vertices per polytope 0051 nt=size(t,1); 0052 cn=v_choosenk(nd+1,2); 0053 nf=size(cn,1); 0054 ee=zeros(nt,nf,2); % space for all the edges 0055 ee(:,:,1)=t(:,cn(:,1)); 0056 ee(:,:,2)=t(:,cn(:,2)); 0057 ee=reshape(ee,[],2); 0058 mk=ee(:,1)>ee(:,2); 0059 ee(mk,:)=ee(mk,2:-1:1); % make all edges in ascending order 0060 % ees=sparse(ee(:,1),ee(:,2),1); % remove duplicates 0061 [er,ec]=find(sparse(ee(:,1),ee(:,2),1)); % remove duplicates 0062 ee=[er,ec]; 0063 ne=size(ee,1); 0064 0065 % now apply Kruskal shortest spanning tree algorithm 0066 0067 sz=sum((x(ee(:,1),:)-x(ee(:,2),:)).^2,2); % length^2 of each edge 0068 [vz,mz]=sort(sz); 0069 ee=ee(mz,:); % sort edges into ascenging length order 0070 ts=ones(ne,1); % size of each component 0071 tp=zeros(np,1); % root node links 0072 ei=zeros(np-1,1); % index of shortest spanning tree edges 0073 k=0; 0074 for i=1:ne 0075 i1=ee(i,1); 0076 j=tp(i1); 0077 while j % find root node for ee(i,1) 0078 i1=j; 0079 j=tp(i1); 0080 end 0081 i2=ee(i,2); 0082 j=tp(i2); 0083 while j % find root node for ee(i,2) 0084 i2=j; 0085 j=tp(i2); 0086 end 0087 if i1~=i2 % if they are different, merge them 0088 k=k+1; 0089 ei(k)=i; % add to the shortest spanning tree 0090 if ts(i1)>ts(i2) 0091 tp(i2)=i1; % make i2 a sub tree of i1 0092 ts(i1)=ts(i1)+ts(i2); 0093 else 0094 tp(i1)=i2; % make i1 a sub tree of i2 0095 ts(i2)=ts(i1)+ts(i2); 0096 end 0097 end 0098 end 0099 ee=ee(ei,:); % refine the edges to include only those in the minimum spanning tree 0100 0101 % now arrange as a tree with point np as the head 0102 0103 eet=sparse(ee(:,1),ee(:,2),1:np-1,np,np); 0104 p=zeros(np-1,1); % points to parent node 0105 s=zeros(np-1,1); % sorted index 0106 chn=np; % start with the root node of the tree 0107 [rf,cf]=find(eet(:,chn)); % find any nodes that connect to it 0108 [rg,cg]=find(eet(chn,:)); 0109 while ~isempty(rf) || ~isempty(rg) 0110 p(rf)=chn(cf); % set the parents 0111 rcf=rf(:)+np*(chn(cf(:))-1); 0112 s(eet(rcf))=rf(:); 0113 eet(rcf)=0; % delete the edges linking them to their parents 0114 p(cg)=chn(rg); 0115 rcg=chn(rg(:))+np*(cg(:)-1); 0116 s(eet(rcg))=cg(:); 0117 eet(rcg)=0; 0118 chn=[rf(:); cg(:)]; % now search for their children 0119 [rf,cf]=find(eet(:,chn)); 0120 [rg,cg]=find(eet(chn,:)); 0121 end 0122