User:Zeracles/Starmap/Minimal Spanning Tree

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