Connect Randomly Ordered 2D Points into a Minimal Nearest-Neighbor Closed Contour

points2contour
Tristan Ursell
Sept 2013

[Xout,Yout]=points2contour(Xin,Yin,P,direction)
[Xout,Yout]=points2contour(Xin,Yin,P,direction,dlim)
[Xout,Yout,orphans]=points2contour(Xin,Yin,P,direction,dlim)
[Xout,Yout,orphans,indout]=points2contour(Xin,Yin,P,direction,dlim)

Given any list of 2D points (Xin,Yin), construct a singly connected nearest-neighbor path in either the 'cw' or 'ccw' directions. The code has been written to handle square and hexagon grid points, as well as any non-grid arrangement of points.

'P' sets the point to begin looking for the contour from the original ordering of (Xin,Yin), and 'direction' sets the direction of the contour, with options 'cw' and 'ccw', specifying clockwise and counter-clockwise, respectively.

The optional input parameter 'dlim' sets a distance limit, if the distance between a point and all other points is greater than or equal to 'dlim', the point is left out of the contour.

The optional output 'orphans' gives the indices of the original (Xin,Yin) points that were not included in the contour.

The optional output 'indout' is the order of indices that produces Xin(indout)=Xout and Yin(indout)=Yout.

There are many (Inf) situations where there is no unique mapping of points into a connected contour -- e.g. any time there are more than 2 nearest neighbor points, or in situations where the nearest neighbor matrix is non-symmetric. Picking a different P will result in a different contour. Likewise, in cases where one point is far from its neighbors, it may be orphaned, and only connected into the path at the end, giving strange results.

The input points can be of any numerical class.

Note that this will *not* necessarily form the shortest path between all the points -- that is the NP-Hard Traveling Salesman Problem, for which there is no deterministic solution. This will, however, find the shortest path for points with a symmetric nearest neighbor matrix.

see also: bwtraceboundary

Example 1: continuous points
N=200;
P=1;
theta=linspace(0,2*pi*(1-1/N),N);
[~,I]=sort(rand(1,N));
R=2+sin(5*theta(I))/3;

Xin=R.*cos(theta(I));
Yin=R.*sin(theta(I));

[Xout,Yout]=points2contour(Xin,Yin,P,'cw');

figure;
hold on
plot(Xin,Yin,'b-')
plot(Xout,Yout,'r-','Linewidth',2)
plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
plot(Xout(1),Yout(1),'g.','Markersize',15)
plot(Xout(end),Yout(end),'r.','Markersize',15)
xlabel('X')
ylabel('Y')
axis equal tight
title(['Black = original points, Blue = original ordering, Red = new ordering, Green = starting points'])
box on

Example 2: square grid
P=1;

Xin=[1,2,3,4,4,4,4,3,2,1,1,1];
Yin=[0,0,0,0,1,2,3,3,2,2,1,0];

[Xout,Yout]=points2contour(Xin,Yin,P,'cw');

figure;
hold on
plot(Xin,Yin,'b-')
plot(Xout,Yout,'r-','Linewidth',2)
plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
plot(Xout(1),Yout(1),'g.','Markersize',15)
plot(Xout(end),Yout(end),'r.','Markersize',15)
xlabel('X')
ylabel('Y')
axis equal tight
box on

Example 3: continuous points, pathological case
N=200;
P=1;
theta=linspace(0,2*pi*(1-1/N),N);
[~,I]=sort(rand(1,N));
R=2+sin(5*theta(I))/3;

Xin=(1+rand(1,N)/2).*R.*cos(theta(I));
Yin=(1+rand(1,N)/2).*R.*sin(theta(I));

[Xout,Yout]=points2contour(Xin,Yin,P,'cw');

figure;
hold on
plot(Xin,Yin,'b-')
plot(Xout,Yout,'r-','Linewidth',2)
plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
plot(Xout(1),Yout(1),'g.','Markersize',15)
plot(Xout(end),Yout(end),'r.','Markersize',15)
xlabel('X')
ylabel('Y')
axis equal tight
title(['Black = original points, Blue = original ordering, Red = new ordering, Green = starting points'])
box on

Example 4: continuous points, distance limit applied
N=200;
P=1;
theta=linspace(0,2*pi*(1-1/N),N);
[~,I]=sort(rand(1,N));
R=2+sin(5*theta(I))/3;
R(2)=5; %the outlier

Xin=(1+rand(1,N)/16).*R.*cos(theta(I));
Yin=(1+rand(1,N)/16).*R.*sin(theta(I));

[Xout,Yout,orphans,indout]=points2contour(Xin,Yin,P,'cw',1);

figure;
hold on
plot(Xin,Yin,'b-')
plot(Xin(orphans),Yin(orphans),'kx')
plot(Xin(indout),Yin(indout),'r-','Linewidth',2)
plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
plot(Xout(1),Yout(1),'g.','Markersize',15)
plot(Xout(end),Yout(end),'r.','Markersize',15)
xlabel('X')
ylabel('Y')
axis equal tight
title(['Black = original points, Blue = original ordering, Red = new ordering, Green = starting points'])
box on