|
|
| (One intermediate revision by the same user not shown) |
| Line 1: |
Line 1: |
| | 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. Demonstration [http://www.star-control.com/forum/index.php/topic,1660.msg35209.html#msg35209 here]. | | 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. Demonstration [http://www.star-control.com/forum/index.php/topic,1660.msg35209.html#msg35209 here]. |
| | ==Minimally Connected Opaque Friend Graph== | | ==Minimally Connected Opaque Friend Graph== |
| − | The code below for Boruvka's algorithm also contains an option to create a Minimally Connected Opaque Friend Graph (MCOFG). This is a graph I invented for the purpose of adding links to a completed MST to get closed loops, which were present in the original Star Control 1 starmaps. Actually, the MCOFG has another very interesting use in my real life work, but that's another story. | + | The matlab/octave code below for Boruvka's algorithm also contains an option to create a Minimally Connected Opaque Friend Graph (MCOFG). This is a graph I invented for the purpose of adding links to a completed MST to get closed loops, which were present in the original Star Control 1 starmaps. Actually, the MCOFG has another very interesting use in my real life work, but that's another story. |
| | | | |
| | The basic properties of the graph are | | The basic properties of the graph are |
| Line 19: |
Line 19: |
| | <br>((r_2)/(2r_1))*((\pi/2)/(\theta))*(1/f) | | <br>((r_2)/(2r_1))*((\pi/2)/(\theta))*(1/f) |
| | <br>to vary the number of links. Higher values of f lead to more connections, more friends. The locality factor l influences the radius within which points are considered for friendship, which is by default the length of the longest local MST link. l = 2, say, would double this. | | <br>to vary the number of links. Higher values of f lead to more connections, more friends. The locality factor l influences the radius within which points are considered for friendship, which is by default the length of the longest local MST link. l = 2, say, would double this. |
| − | ==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=MCOFG=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
| |
| − |
| |
| − | MCOFG=0; % to create Minimally Connected Opaque Friend Graph (superset of MST)
| |
| − | friendliness=1; % default value is one
| |
| − | localityfactor=2; % default value is one
| |
| − |
| |
| − | % detect input properties and set up index
| |
| − | sizeinputpositions=size(inputpositions);
| |
| − | lengthinputpositions=sizeinputpositions(1);
| |
| − | dimensions=sizeinputpositions(2);
| |
| − |
| |
| − | index=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;
| |
| − | if xextent==0
| |
| − | xextent=eps;
| |
| − | end
| |
| − |
| |
| − | miny=min(inputpositionys);
| |
| − | maxy=max(inputpositionys);
| |
| − | yextent=miny*maxy;
| |
| − | if yextent==0
| |
| − | yextent=eps;
| |
| − | end
| |
| − |
| |
| − | if dimensions>2
| |
| − | minz=min(inputpositionzs);
| |
| − | maxz=max(inputpositionzs);
| |
| − | zextent=minz*maxz;
| |
| − | if zextent==0
| |
| − | zextent=eps;
| |
| − | end
| |
| − | 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(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
| |
| − |
| |
| − | if MCOFG==1
| |
| − | sizebranches=size(branches);
| |
| − | lengthbranches=sizebranches(1);
| |
| − | branches=[branches,repmat(1,[lengthbranches 1])]; % designate existing links as MST links (type 1)
| |
| − |
| |
| − | localmaxs=[];
| |
| − | for i=linspace(1,lengthinputpositions,lengthinputpositions)
| |
| − | MSTlengths=branches(:,1);
| |
| − | branchlist=((branches(:,2))==i)+((branches(:,3))==i);
| |
| − | localMSTlengths=MSTlengths(find(branchlist==1));
| |
| − | localmax=max(localMSTlengths); % get length of longest adjacent MST link
| |
| − | localmaxs=[localmaxs;localmax];
| |
| − | end
| |
| − | localmaxs=localityfactor*localmaxs;
| |
| − | for i=linspace(1,lengthinputpositions,lengthinputpositions)
| |
| − | localmax=localmaxs(i);
| |
| − |
| |
| − | neighbours=[branches(:,2),branches(:,3)];
| |
| − | branchlist=((branches(:,2))==i)+((branches(:,3))==i);
| |
| − |
| |
| − | neighbours=neighbours(find(branchlist==1),:);
| |
| − | neighbours=[neighbours(:,1);neighbours(:,2)];
| |
| − |
| |
| − | localcentrex=inputpositionxs(i);
| |
| − | localcentrey=inputpositionys(i);
| |
| − | if dimensions>2
| |
| − | localcentrez=inputpositionzs(i);
| |
| − | end
| |
| − |
| |
| − | % retrieve all points within radius localmax
| |
| − | boxminx=localcentrex-localmax;
| |
| − | boxmaxx=localcentrex+localmax;
| |
| − | boxminy=localcentrey-localmax;
| |
| − | boxmaxy=localcentrey+localmax;
| |
| − | if dimensions>2
| |
| − | boxminz=localcentrez-localmax;
| |
| − | boxmaxz=localcentrez+localmax;
| |
| − | end
| |
| − |
| |
| − | if dimensions==2
| |
| − | closemask=(boxminx<=inputpositionxs).*(inputpositionxs<=boxmaxx).*(boxminy<=inputpositionys).*(inputpositionys<=boxmaxy);
| |
| − | end
| |
| − | if dimensions==3
| |
| − | closemask=(boxminx<=inputpositionxs).*(inputpositionxs<=boxmaxx).*(boxminy<=inputpositionys).*(inputpositionys<=boxmaxy).*(boxminz<=inputpositionzs).*(inputpositionzs<=boxmaxz);
| |
| − | end
| |
| − |
| |
| − | closedistances=[];
| |
| − | localindices=[];
| |
| − |
| |
| − | xcomponents=[];
| |
| − | ycomponents=[];
| |
| − | if dimensions>2
| |
| − | zcomponents=[];
| |
| − | end
| |
| − |
| |
| − | for j=transpose(find(closemask==1))
| |
| − | otherpointx=inputpositionxs(j);
| |
| − | otherpointy=inputpositionys(j);
| |
| − | if dimensions>2
| |
| − | otherpointz=inputpositionzs(j);
| |
| − | end
| |
| − |
| |
| − | xcomponent=otherpointx-localcentrex;
| |
| − | ycomponent=otherpointy-localcentrey;
| |
| − | if dimensions>2
| |
| − | zcomponent=otherpointz-localcentrez;
| |
| − | end
| |
| − |
| |
| − | if dimensions==2
| |
| − | closedistance=((xcomponent).^2+(ycomponent).^2).^(1/2);
| |
| − | end
| |
| − | if dimensions==3
| |
| − | closedistance=((xcomponent).^2+(ycomponent).^2+(zcomponent).^2).^(1/2);
| |
| − | end
| |
| − |
| |
| − | if closedistance<=localmax
| |
| − | closedistances=[closedistances;closedistance];
| |
| − | localindices=[localindices;j];
| |
| − |
| |
| − | xcomponents=[xcomponents;xcomponent];
| |
| − | ycomponents=[ycomponents;ycomponent];
| |
| − | if dimensions>2
| |
| − | zcomponents=[zcomponents;zcomponent];
| |
| − | end
| |
| − | end
| |
| − | end
| |
| − |
| |
| − | % remove local centre point from prospective linkees
| |
| − | notcentre=find(closedistances>(min(closedistances)));
| |
| − |
| |
| − | closedistances=closedistances(notcentre);
| |
| − | localindices=localindices(notcentre);
| |
| − | xcomponents=xcomponents(notcentre);
| |
| − | ycomponents=ycomponents(notcentre);
| |
| − | if dimensions>2
| |
| − | zcomponents=zcomponents(notcentre);
| |
| − | end
| |
| − |
| |
| − | sizeclosedistances=size(closedistances);
| |
| − | lengthclosedistances=sizeclosedistances(1);
| |
| − | for j=linspace(1,lengthclosedistances,lengthclosedistances)
| |
| − | closedistance=closedistances(j);
| |
| − | localindex=localindices(j);
| |
| − |
| |
| − | xcomponent=xcomponents(j);
| |
| − | ycomponent=ycomponents(j);
| |
| − | if dimensions>2
| |
| − | zcomponent=zcomponents(j);
| |
| − | end
| |
| − |
| |
| − | if sum(localindex==neighbours)==0 % if there's already a link to this other point, or it's the same as our local centre point, leave it alone
| |
| − | if sum(closedistances<closedistance)>0 % if there are actually points that are closer, apart from the local centre point
| |
| − | inclusionfactor=1; % ``inclusion" within a modified radius . . . sort of
| |
| − | for k=transpose(find(closedistances<closedistance));
| |
| − | if dimensions==2
| |
| − | dotproduct=(xcomponent*xcomponents(k))+(ycomponent*ycomponents(k));
| |
| − | end
| |
| − | if dimensions==3
| |
| − | dotproduct=(xcomponent*xcomponents(k))+(ycomponent*ycomponents(k))+(zcomponent*zcomponents(k));
| |
| − | end
| |
| − |
| |
| − | angle=acos(dotproduct/(closedistance*closedistances(k)));
| |
| − | inclusionfactor=inclusionfactor*((closedistance/2)/(closedistances(k)))*((pi/2)/angle)/friendliness; % I think multiplying the inclusion factor is better than some additive scheme, since this way one point can block a connection even if there are many other points working against this effect; also, dividing by the friendliness at each step might be best to make this work the same way at all scales. That the presence of neighbours can ``create" connections on the opposite side of the centre point is also intuitively acceptable, in my view
| |
| − | end
| |
| − | if inclusionfactor<=1
| |
| − | branches=[branches;[closedistance,i,localindex,2]];
| |
| − | end
| |
| − | end
| |
| − | end
| |
| − | end
| |
| − | 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== | + | ==Efficiency Method== |
| − | Also downloadable [http://www.physics.usyd.edu.au/~asmith/files/sc1/MST_JarnikPrim.m here]
| + | In discussion on the forum Death 999 contributed a method that has become our preferred method for closed loop generation - I'll insert a description when I get a moment. It is implemented in both versions of the code available below. |
| − | % minimal spannng tree in matlab
| + | |
| − | % Jarnik/Prim's algorithm - runs pretty slowly for any more than 50 input points
| + | ==Code== |
| − | % see http://www.star-control.com/forum/index.php/topic,1660.msg33938.html#msg33938
| + | When prototyping, I wrote a [http://www.physics.usyd.edu.au/~asmith/files/sc1/MSTclosure.m matlab/octave version available here], and have recently done a [http://www.physics.usyd.edu.au/~asmith/files/sc1/MSTclosure.c C conversion available here]. |
| − | % 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
| |