Difference between revisions of "User:Zeracles/Starmap/Minimal Spanning Tree"

From Ultronomicon
Jump to navigation Jump to search
m (renamed target and added line)
(Boruvka's algorithm is faster, at least the way I've done it)
Line 1: Line 1:
−
=Jarnik/Prim's Algorithm=
+
I originally coded up [[wikipedia:Prim's_algorithm|Jarnik/Prim's algorithm]]. In the process of optimising it, I unwittingly reinvented [[wikipedia:Bor%C5%AFvka%27s_algorithm|Boruvka's algorithm]], which handles large numbers of input points better, at least the way I've done it.
 +
==Boruvka's Algorithm==
 +
Also downloadable [http://www.physics.usyd.edu.au/~asmith/files/sc1/MST_Boruvka.m here]
 +
% minimal spannng tree in matlab
 +
% Boruvka's algorithm (independently invented by yours truly)
 +
% see http://www.star-control.com/forum/index.php/topic,1660
 +
% by Anthony Smith (Zeracles)
 +
% astroant@hotmail.com
 +
 +
% if only after an MST, leave separation=pruning=0
 +
 +
load('inputpositions')
 +
 +
localcount=8;
 +
 +
separation=0;
 +
separationlength=0.03;
 +
 +
pruning=0;
 +
prunesteps=1; % how many times to have a go at the tree
 +
prunelevel=10; % length of branch to cut each time
 +
 +
% detect input properties and set up index
 +
sizeinputpositions=size(inputpositions);
 +
lengthinputpositions=sizeinputpositions(1);
 +
dimensions=sizeinputpositions(2);
 +
 +
index=transpose(linspace(1,lengthinputpositions,lengthinputpositions));
 +
 +
inputpositionxs=inputpositions(:,1);
 +
inputpositionys=inputpositions(:,2);
 +
if dimensions>2
 +
inputpositionzs=inputpositions(:,3);
 +
end
 +
 +
% initialise search parameters
 +
minx=min(inputpositionxs);
 +
maxx=max(inputpositionxs);
 +
xextent=minx*maxx;
 +
 +
miny=min(inputpositionys);
 +
maxy=max(inputpositionys);
 +
yextent=miny*maxy;
 +
 +
if dimensions>2
 +
minz=min(inputpositionzs);
 +
maxz=max(inputpositionzs);
 +
zextent=minz*maxz;
 +
end
 +
 +
if dimensions==2
 +
measure=xextent*yextent;
 +
localmeasure=measure/lengthinputpositions;
 +
searchincrement=((localcount*localmeasure)/pi)^(1/2);
 +
end
 +
if dimensions==3
 +
measure=xextent*yextent*zextent;
 +
localmeasure=measure/lengthinputpositions;
 +
searchincrement=((localcount*localmeasure)/(4/3)/pi)^(1/3);
 +
end
 +
 +
treenumber=index;
 +
numbertrees=length(unique(treenumber));
 +
branches=[];
 +
searchfactor=1; % setting this to zero can lead to infinite loop if this is scaled multiplicatively
 +
searchfactors=zeros(lengthinputpositions,1);
 +
 +
while numbertrees>1 % each input point is treated as a separate tree, and we join them together until we only have one tree
 +
searchfactor=2*searchfactor; % calculate interpoint distances for ``close" points, with ``close" defined by this variable
 +
while min(searchfactors)<searchfactor % searchfactors stores the searchfactors searched for each existing tree
 +
tree=treenumber(min(find(searchfactors<(searchfactor)))); % choose a tree to search which has been searched with the smallest distance so far
 +
treesearched=0;
 +
while treesearched==0 % while we're not yet done looking around the nodes of this tree for links
 +
maskvector=treenumber==tree;
 +
treeindices=find(maskvector);
 +
lengthtreepositions=sum(maskvector);
 +
 +
remainingpositionxs=inputpositionxs(find(maskvector==0));
 +
remainingpositionys=inputpositionys(find(maskvector==0));
 +
if dimensions>2
 +
remainingpositionzs=inputpositionzs(find(maskvector==0));
 +
end
 +
 +
remainingpositionindices=index(find(maskvector==0));
 +
 +
mindistances=[];
 +
for i=linspace(1,lengthtreepositions,lengthtreepositions) % search around each of the tree positions for remainingpositions
 +
% retrieve coordinates of a tree position
 +
treepositionx=inputpositions(treeindices(i),1);
 +
treepositiony=inputpositions(treeindices(i),2);
 +
if dimensions>2
 +
treepositionz=inputpositions(treeindices(i),3);
 +
end
 +
treeindex=treeindices(i);
 +
 +
% search for remainingpositions
 +
boxradius=searchfactor*searchincrement;
 +
 +
boxminx=treepositionx-boxradius;
 +
boxmaxx=treepositionx+boxradius;
 +
boxminy=treepositiony-boxradius;
 +
boxmaxy=treepositiony+boxradius;
 +
if dimensions>2
 +
boxminz=treepositionz-boxradius;
 +
boxmaxz=treepositionz+boxradius;
 +
end
 +
 +
if dimensions==2
 +
closerremainingmask=(boxminx<=remainingpositionxs).*(remainingpositionxs<=boxmaxx).*(boxminy<=remainingpositionys).*(remainingpositionys<=boxmaxy);
 +
end
 +
if dimensions==3
 +
closerremainingmask=(boxminx<=remainingpositionxs).*(remainingpositionxs<=boxmaxx).*(boxminy<=remainingpositionys).*(remainingpositionys<=boxmaxy).*(boxminz<=remainingpositionzs).*(remainingpositionzs<=boxmaxz);
 +
end
 +
 +
mindistancecandidates=[];
 +
mindistancecandidateindices=[];
 +
if sum(closerremainingmask)>0 % if there are remainingpositions inside the box, calculate the actual distances to them
 +
for k=transpose(find(closerremainingmask==1));
 +
branchcandidatex=remainingpositionxs(k);
 +
branchcandidatey=remainingpositionys(k);
 +
if dimensions>2
 +
branchcandidatez=remainingpositionzs(k);
 +
end
 +
branchcandidateindex=remainingpositionindices(k); % set up to remember the index of each prospective link destination
 +
 +
if dimensions==2
 +
mindistancecandidate=((branchcandidatex-treepositionx).^2+(branchcandidatey-treepositiony).^2).^(1/2);
 +
end
 +
if dimensions==3
 +
mindistancecandidate=((branchcandidatex-treepositionx).^2+(branchcandidatey-treepositiony).^2+(branchcandidatez-treepositionz).^2).^(1/2);
 +
end
 +
if mindistancecandidate<=boxradius
 +
mindistancecandidates=[mindistancecandidates;mindistancecandidate];
 +
mindistancecandidateindices=[mindistancecandidateindices;branchcandidateindex];
 +
end
 +
end
 +
end
 +
if numel(mindistancecandidates)>0 % minimum distance from this position on the tree
 +
mindistance=[min(mindistancecandidates),treeindex,mindistancecandidateindices(find(mindistancecandidates==min(mindistancecandidates)))]; % minimum distances carry the indices of the prospective link destination
 +
mindistances=[mindistances;mindistance];
 +
end
 +
end
 +
if numel(mindistances)==0
 +
searchfactors(maskvector)=searchfactor;
 +
treesearched=1;
 +
end
 +
if numel(mindistances)>0
 +
branchdistance=min(mindistances(:,1));
 +
branchindex=find(branchdistance==mindistances(:,1));
 +
numberbranchcandidates=numel(branchindex);
 +
if numberbranchcandidates>1 % if there's a tie for minimum length
 +
branchindex=branchindex(ceil(numberbranchcandidates*rand(1)));
 +
end
 +
 +
% add this branch to the tree
 +
branch=mindistances(branchindex,:);
 +
branches=[branches;branch];
 +
 +
% join the two trees together
 +
targettree=treenumber(branch(3)); % need to update the treenumber variable, which remembers which points belong to which trees
 +
treenumber(find(treenumber==targettree))=tree;
 +
numbertrees=numbertrees-1; % so that we are one step closer to finishing
 +
end
 +
end
 +
end
 +
end
 +
 +
% do custom stuff with this tree
 +
% separation - remove all branches above a certain length
 +
if separation==1
 +
separationmask=(branches(:,1))<=separationlength;
 +
branches=branches(find(separationmask),:);
 +
end
 +
 +
% pruning - remove wandering branches
 +
if pruning==1
 +
degrees=[];
 +
for i=linspace(1,lengthinputpositions,lengthinputpositions)
 +
degree=(sum((branches(:,2))==i))+(sum((branches(:,3))==i));
 +
degrees=[degrees;degree];
 +
end
 +
 +
while prunesteps>0
 +
protectednodes=index(find(degrees>2));
 +
 +
tempprunelevel=prunelevel;
 +
while tempprunelevel>0
 +
ends=index(find(degrees==1)); % find branch end nodes
 +
for i=transpose(ends)
 +
if sum(i==protectednodes)==0
 +
if (degrees(i))>0
 +
deletemask=((branches(:,2))==i)+((branches(:,3))==i);
 +
cutbranch=branches(find(deletemask),:);
 +
 +
branches=branches(find(deletemask==0),:);
 +
degrees(cutbranch(2))=(degrees(cutbranch(2)))-1;
 +
degrees(cutbranch(3))=(degrees(cutbranch(3)))-1;
 +
end
 +
end
 +
end
 +
tempprunelevel=tempprunelevel-1;
 +
end
 +
prunesteps=prunesteps-1;
 +
end
 +
end
 +
 +
save('branches','branches')
 +
 +
% for plotting
 +
sizebranches=size(branches);
 +
numberofbranches=sizebranches(1);
 +
 +
figure
 +
hold on
 +
box on
 +
if dimensions==2
 +
plot(inputpositions(:,1),inputpositions(:,2),'k.')
 +
end
 +
if dimensions==3
 +
plot3(inputpositions(:,1),inputpositions(:,2),inputpositions(:,3),'k.')
 +
end
 +
 +
for i=linspace(1,numberofbranches,numberofbranches)
 +
if dimensions==2
 +
startindex=branches(i,2);
 +
finishindex=branches(i,3);
 +
 +
startx=inputpositions(startindex,1);
 +
starty=inputpositions(startindex,2);
 +
 +
finishx=inputpositions(finishindex,1);
 +
finishy=inputpositions(finishindex,2);
 +
 +
line([startx;finishx],[starty;finishy],'Color','k')
 +
end
 +
if dimensions==3
 +
startindex=branches(i,2);
 +
finishindex=branches(i,3);
 +
 +
startx=inputpositions(startindex,1);
 +
starty=inputpositions(startindex,2);
 +
startz=inputpositions(startindex,3);
 +
 +
finishx=inputpositions(finishindex,1);
 +
finishy=inputpositions(finishindex,2);
 +
finishz=inputpositions(finishindex,3);
 +
 +
line([startx;finishx],[starty;finishy],[startz;finishz],'Color','k')
 +
end
 +
end
 +
==Jarnik/Prim's Algorithm==
 
Also downloadable [http://www.physics.usyd.edu.au/~asmith/files/sc1/MST_JarnikPrim.m here]
 
Also downloadable [http://www.physics.usyd.edu.au/~asmith/files/sc1/MST_JarnikPrim.m here]
 
  % minimal spannng tree in matlab
 
  % minimal spannng tree in matlab

Revision as of 14:49, 14 February 2010

I originally coded up Jarnik/Prim's algorithm. In the process of optimising it, I unwittingly reinvented Boruvka's algorithm, which handles large numbers of input points better, at least the way I've done it.

Boruvka's Algorithm

Also downloadable here 

%	minimal spannng tree in matlab
%		Boruvka's algorithm (independently invented by yours truly)
%	see http://www.star-control.com/forum/index.php/topic,1660
%	by Anthony Smith (Zeracles)
%	astroant@hotmail.com

%	if only after an MST, leave separation=pruning=0

load('inputpositions')

localcount=8;

separation=0;
separationlength=0.03;

pruning=0;
prunesteps=1; % how many times to have a go at the tree
prunelevel=10; % length of branch to cut each time

%	detect input properties and set up index
sizeinputpositions=size(inputpositions);
lengthinputpositions=sizeinputpositions(1);
dimensions=sizeinputpositions(2);

index=transpose(linspace(1,lengthinputpositions,lengthinputpositions));

inputpositionxs=inputpositions(:,1);
inputpositionys=inputpositions(:,2);
if dimensions>2
	inputpositionzs=inputpositions(:,3);
end

%	initialise search parameters
minx=min(inputpositionxs);
maxx=max(inputpositionxs);
xextent=minx*maxx;

miny=min(inputpositionys);
maxy=max(inputpositionys);
yextent=miny*maxy;

if dimensions>2
	minz=min(inputpositionzs);
	maxz=max(inputpositionzs);
	zextent=minz*maxz;
end

if dimensions==2
	measure=xextent*yextent;
	localmeasure=measure/lengthinputpositions;
	searchincrement=((localcount*localmeasure)/pi)^(1/2);
end
if dimensions==3
	measure=xextent*yextent*zextent;
	localmeasure=measure/lengthinputpositions;
	searchincrement=((localcount*localmeasure)/(4/3)/pi)^(1/3);
end

treenumber=index;
numbertrees=length(unique(treenumber));
branches=[];
searchfactor=1; % setting this to zero can lead to infinite loop if this is scaled multiplicatively
searchfactors=zeros(lengthinputpositions,1);

while numbertrees>1 % each input point is treated as a separate tree, and we join them together until we only have one tree
	searchfactor=2*searchfactor; % calculate interpoint distances for ``close" points, with ``close" defined by this variable
	while min(searchfactors)<searchfactor % searchfactors stores the searchfactors searched for each existing tree
		tree=treenumber(min(find(searchfactors<(searchfactor)))); % choose a tree to search which has been searched with the smallest distance so far
		treesearched=0;
		while treesearched==0 % while we're not yet done looking around the nodes of this tree for links
			maskvector=treenumber==tree;
			treeindices=find(maskvector);
			lengthtreepositions=sum(maskvector);
			
			remainingpositionxs=inputpositionxs(find(maskvector==0));
			remainingpositionys=inputpositionys(find(maskvector==0));
			if dimensions>2
				remainingpositionzs=inputpositionzs(find(maskvector==0));
			end
			
			remainingpositionindices=index(find(maskvector==0));
			
			mindistances=[];
			for i=linspace(1,lengthtreepositions,lengthtreepositions) % search around each of the tree positions for remainingpositions
				%	retrieve coordinates of a tree position
				treepositionx=inputpositions(treeindices(i),1);
				treepositiony=inputpositions(treeindices(i),2);
				if dimensions>2
					treepositionz=inputpositions(treeindices(i),3);
				end
				treeindex=treeindices(i);
				
				%	search for remainingpositions
				boxradius=searchfactor*searchincrement;
				
				boxminx=treepositionx-boxradius;
				boxmaxx=treepositionx+boxradius;
				boxminy=treepositiony-boxradius;
				boxmaxy=treepositiony+boxradius;
				if dimensions>2
					boxminz=treepositionz-boxradius;
					boxmaxz=treepositionz+boxradius;
				end
				
				if dimensions==2
					closerremainingmask=(boxminx<=remainingpositionxs).*(remainingpositionxs<=boxmaxx).*(boxminy<=remainingpositionys).*(remainingpositionys<=boxmaxy);
				end
				if dimensions==3
					closerremainingmask=(boxminx<=remainingpositionxs).*(remainingpositionxs<=boxmaxx).*(boxminy<=remainingpositionys).*(remainingpositionys<=boxmaxy).*(boxminz<=remainingpositionzs).*(remainingpositionzs<=boxmaxz);
				end
				
				mindistancecandidates=[];
				mindistancecandidateindices=[];
				if sum(closerremainingmask)>0 % if there are remainingpositions inside the box, calculate the actual distances to them
					for k=transpose(find(closerremainingmask==1));
						branchcandidatex=remainingpositionxs(k);
						branchcandidatey=remainingpositionys(k);
						if dimensions>2
							branchcandidatez=remainingpositionzs(k);
						end
						branchcandidateindex=remainingpositionindices(k); % set up to remember the index of each prospective link destination
						
						if dimensions==2
							mindistancecandidate=((branchcandidatex-treepositionx).^2+(branchcandidatey-treepositiony).^2).^(1/2);
						end
						if dimensions==3
							mindistancecandidate=((branchcandidatex-treepositionx).^2+(branchcandidatey-treepositiony).^2+(branchcandidatez-treepositionz).^2).^(1/2);
						end
						if mindistancecandidate<=boxradius
							mindistancecandidates=[mindistancecandidates;mindistancecandidate];
							mindistancecandidateindices=[mindistancecandidateindices;branchcandidateindex];
						end
					end
				end
				if numel(mindistancecandidates)>0 % minimum distance from this position on the tree
					mindistance=[min(mindistancecandidates),treeindex,mindistancecandidateindices(find(mindistancecandidates==min(mindistancecandidates)))]; % minimum distances carry the indices of the prospective link destination
					mindistances=[mindistances;mindistance];
				end
			end
			if numel(mindistances)==0
				searchfactors(maskvector)=searchfactor;
				treesearched=1;
			end
			if numel(mindistances)>0
				branchdistance=min(mindistances(:,1));
				branchindex=find(branchdistance==mindistances(:,1));
				numberbranchcandidates=numel(branchindex);
				if numberbranchcandidates>1 % if there's a tie for minimum length
					branchindex=branchindex(ceil(numberbranchcandidates*rand(1)));
				end
				
				%	add this branch to the tree
				branch=mindistances(branchindex,:);
				branches=[branches;branch];
				
				%	join the two trees together
				targettree=treenumber(branch(3)); % need to update the treenumber variable, which remembers which points belong to which trees
				treenumber(find(treenumber==targettree))=tree;
				numbertrees=numbertrees-1; % so that we are one step closer to finishing
			end
		end
	end
end

%	do custom stuff with this tree
%		separation - remove all branches above a certain length
if separation==1
	separationmask=(branches(:,1))<=separationlength;
	branches=branches(find(separationmask),:);
end

%		pruning - remove wandering branches
if pruning==1
	degrees=[];
	for i=linspace(1,lengthinputpositions,lengthinputpositions)
		degree=(sum((branches(:,2))==i))+(sum((branches(:,3))==i));
		degrees=[degrees;degree];
	end
	
	while prunesteps>0
		protectednodes=index(find(degrees>2));
		
		tempprunelevel=prunelevel;
		while tempprunelevel>0
			ends=index(find(degrees==1)); % find branch end nodes
			for i=transpose(ends)
				if sum(i==protectednodes)==0
					if (degrees(i))>0
						deletemask=((branches(:,2))==i)+((branches(:,3))==i);
						cutbranch=branches(find(deletemask),:);
						
						branches=branches(find(deletemask==0),:);
						degrees(cutbranch(2))=(degrees(cutbranch(2)))-1;
						degrees(cutbranch(3))=(degrees(cutbranch(3)))-1;
					end
				end
			end
			tempprunelevel=tempprunelevel-1;
		end
		prunesteps=prunesteps-1;
	end
end

save('branches','branches')

%	for plotting
sizebranches=size(branches);
numberofbranches=sizebranches(1);

figure
hold on
box on
if dimensions==2
	plot(inputpositions(:,1),inputpositions(:,2),'k.')
end
if dimensions==3
	plot3(inputpositions(:,1),inputpositions(:,2),inputpositions(:,3),'k.')
end

for i=linspace(1,numberofbranches,numberofbranches)
	if dimensions==2
		startindex=branches(i,2);
		finishindex=branches(i,3);
		
		startx=inputpositions(startindex,1);
		starty=inputpositions(startindex,2);
		
		finishx=inputpositions(finishindex,1);
		finishy=inputpositions(finishindex,2);
		
		line([startx;finishx],[starty;finishy],'Color','k')
	end
	if dimensions==3
		startindex=branches(i,2);
		finishindex=branches(i,3);
		
		startx=inputpositions(startindex,1);
		starty=inputpositions(startindex,2);
		startz=inputpositions(startindex,3);
		
		finishx=inputpositions(finishindex,1);
		finishy=inputpositions(finishindex,2);
		finishz=inputpositions(finishindex,3);
		
		line([startx;finishx],[starty;finishy],[startz;finishz],'Color','k')
	end
end

Jarnik/Prim's Algorithm

Also downloadable here 

%	minimal spannng tree in matlab
%		Jarnik/Prim's algorithm - runs pretty slowly for any more than 50 input points
%	see http://www.star-control.com/forum/index.php/topic,1660.msg33938.html#msg33938
%	by Anthony Smith (Zeracles)
%	astroant@hotmail.com
load('inputpositions')

remainingpositions=inputpositions;

%	detect input properties and set up index
sizeremainingpositions=size(remainingpositions);
lengthremainingpositions=sizeremainingpositions(1);
dimensions=sizeremainingpositions(2);

index=transpose(linspace(1,lengthremainingpositions,lengthremainingpositions));

remainingpositionxs=remainingpositions(:,1);
remainingpositionys=remainingpositions(:,2);
if dimensions>2
	remainingpositionzs=remainingpositions(:,3);
end

%	initialise search parameters
minx=min(remainingpositionxs);
maxx=max(remainingpositionxs);
xextent=minx*maxx;

miny=min(remainingpositionys);
maxy=max(remainingpositionys);
yextent=miny*maxy;

if dimensions>2
	minz=min(remainingpositionzs);
	maxz=max(remainingpositionzs);
	zextent=minz*maxz;
end

if dimensions==2
	measure=xextent*yextent;
	localmeasure=measure/lengthremainingpositions;
	searchincrement=((10*localmeasure)/pi)^(1/2);
end
if dimensions==3
	measure=xextent*yextent*zextent;
	localmeasure=measure/lengthremainingpositions;
	searchincrement=((10*localmeasure)/(4/3)/pi)^(1/3);
end

%	initialise with a seed point
seedindex=1;
maskvector=zeros(lengthremainingpositions,1);
maskvector(seedindex)=1;

remainingpositionxs=remainingpositionxs(find(maskvector==0));
remainingpositionys=remainingpositionys(find(maskvector==0));
if dimensions>2
	remainingpositionzs=remainingpositionzs(find(maskvector==0));
end
remainingpositionindices=index(find(maskvector==0));

lengthremainingpositions=lengthremainingpositions-1;

treeindices=index(seedindex,:);
lengthtreepositions=1;

branches=[];

while lengthremainingpositions>0 % keep going until nothing remains unlinked to the tree (the ``remainingpositions")
	mindistances=[];
	for i=linspace(1,lengthtreepositions,lengthtreepositions) % search around each of the tree positions for remainingpositions
		%	retrieve coordinates of a tree position
		treepositionx=inputpositions(treeindices(i),1);
		treepositiony=inputpositions(treeindices(i),2);
		if dimensions>2
			treepositionz=inputpositions(treeindices(i),3);
		end
		treeindex=treeindices(i);
		
		occupiedbox=0;
		searchfactor=1;
		while occupiedbox==0 % search for remainingpositions inside a box that grows until some are found
			boxradius=searchfactor*searchincrement;
			
			boxminx=treepositionx-boxradius;
			boxmaxx=treepositionx+boxradius;
			boxminy=treepositiony-boxradius;
			boxmaxy=treepositiony+boxradius;
			if dimensions>2
				boxminz=treepositionz-boxradius;
				boxmaxz=treepositionz+boxradius;
			end
			
			if dimensions==2
				closerremainingmask=(boxminx<=remainingpositionxs).*(remainingpositionxs<=boxmaxx).*(boxminy<=remainingpositionys).*(remainingpositionys<=boxmaxy);
			end
			if dimensions==3
				closerremainingmask=(boxminx<=remainingpositionxs).*(remainingpositionxs<=boxmaxx).*(boxminy<=remainingpositionys).*(remainingpositionys<=boxmaxy).*(boxminz<=remainingpositionzs).*(remainingpositionzs<=boxmaxz);
			end
			
			mindistancecandidates=[];
			mindistancecandidateindices=[];
			if sum(closerremainingmask)>0 % if there are remainingpositions inside the box, calculate the actual distances to them
				for k=find(closerremainingmask==1)
					branchcandidatex=remainingpositionxs(k);
					branchcandidatey=remainingpositionys(k);
					if dimensions>2
						branchcandidatez=remainingpositionzs(k);
					end
					branchcandidateindex=remainingpositionindices(k); % set up to remember the index of each prospective link destination
					
					if dimensions==2
						mindistancecandidate=((branchcandidatex-treepositionx).^2+(branchcandidatey-treepositiony).^2).^(1/2);
					end
					if dimensions==3
						mindistancecandidate=((branchcandidatex-treepositionx).^2+(branchcandidatey-treepositiony).^2+(branchcandidatez-treepositionz).^2).^(1/2);
					end
					if mindistancecandidate<=boxradius
						mindistancecandidates=[mindistancecandidates;mindistancecandidate];
						mindistancecandidateindices=[mindistancecandidateindices;branchcandidateindex];
					end
				end
			end
			if numel(mindistancecandidates)>0
				mindistance=[min(mindistancecandidates),treeindex,mindistancecandidateindices(find(mindistancecandidates==min(mindistancecandidates)))]; % minimum distances carry the indices of the prospective link destination
				mindistances=[mindistances;mindistance];
				
				occupiedbox=1;
			end
			searchfactor=searchfactor+1; % make the search box bigger
		end
	end
	branchdistance=min(mindistances(:,1));
	branchindex=find(branchdistance==mindistances(:,1));
	numberbranchcandidates=numel(branchindex);
	if numberbranchcandidates>1 % if there's a tie for minimum length
		branchindex=branchindex(ceil(numberbranchcandidates*rand(1)));
	end
	
	%	add this branch to the tree
	branch=mindistances(branchindex,:);
	branches=[branches;branch];
	
	%	add that location to the tree
	treeindices=[treeindices;branch(3)];
	lengthtreepositions=lengthtreepositions+1;
	
	%	remove that location from the remainingpositions
	remainingmask=ones(lengthremainingpositions,1);
	removed=find(remainingpositionindices==branch(3));
	remainingmask(removed)=0;
	
	remainingpositionindices=remainingpositionindices(find(remainingmask));
	remainingpositionxs=remainingpositionxs(find(remainingmask));
	remainingpositionys=remainingpositionys(find(remainingmask));
	if dimensions>2
		remainingpositionzs=remainingpositionzs(find(remainingmask));
	end

	lengthremainingpositions=lengthremainingpositions-1;
end

save('branches','branches')

%	for plotting
sizebranches=size(branches);
numberofbranches=sizebranches(1);

figure
hold on
box on
plot(inputpositions(:,1),inputpositions(:,2),'k.')

for i=linspace(1,numberofbranches,numberofbranches)
	startindex=branches(i,2);
	finishindex=branches(i,3);
	
	startx=inputpositions(startindex,1);
	starty=inputpositions(startindex,2);
	
	finishx=inputpositions(finishindex,1);
	finishy=inputpositions(finishindex,2);
	
	line([startx;finishx],[starty;finishy],'Color','k')
end
Retrieved from "https://wiki.uqm.stack.nl/index.php?title=User:Zeracles/Starmap/Minimal_Spanning_Tree&oldid=23873"