File Exchange

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

version 1.5 (3.61 KB) by

Connects randomly ordered 2D points into a minimal nearest neighbor contour.

Updated

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.

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

Almog

Bing Liu

### Bing Liu (view profile)

I had not carefully test "the bug" I mentioned in the previous post, If I find time, I will check it out.

Tristan Ursell

### Tristan Ursell (view profile)

Hi Bing -- the 'bad case' you show below is simply a case where the order you originally present to the function is not the order of shortest distances between neighboring points. The point one chooses to start from can have a significant effect on the outcome of the identified contour. For instance, if you start with P=6, you'll get "what you think you should get" for the data set you give -- the function is working fine. The general problem of trying to find the shortest path through N points is NP hard and is, in general, not computationally tractable. This is a nearest neighbor approximation.

Bing Liu

### Bing Liu (view profile)

Xin = [123.7302214;120.7407074;117.7265833;114.2966447;110.9905859;108.3515822;54.5009171;053.24921186;52.60716221;50.69945203;50.89511149;53.89977212;60.19004599;67.18463273;132.0701569;132.7761246;127.0586296];

Yin = [296.0433838;301.0230798;305.2253093;308.6279995;311.5592887;314.9352201;317.9818372;310.329296;303.0774656;296.0433801;288.707404;281.9942383;277.3895304;274.8017206;273.1972414;281.3838842;289.9496757];

Actually the input ordered as a contour, after Point2Contour function, it is crossed.

Bing Liu

### Bing Liu (view profile)

@Tristan Ursell thank you very much for this function, I am no sure is there a bug for this line of code.

indout=circshift(temp1(indout0),[P,0]);

should it be

indout=circshift(temp1(indout0),[1 - P,0]);

because I am a newbie in this area, so it is just a guess...

thank you so much for your wonderful work again.

Tristan Ursell

### Tristan Ursell (view profile)

@Frank -- That's an interesting idea, but I worry that it's too dependent on which point you choose as the first point. In other words, for the same vector, different starting points would give different results (although the same is true for the algorithm in general). Probably the right way to do something like that is calculate a distance matrix for all points and then do a bit of graph theory to group the points by mutual distance. Give it a try :)

Frank Pennekamp

### Frank Pennekamp (view profile)

Hi,
This is a great function you've made here! But I was wondering if it is possible to find multiple closed contours in one X and Y vector?

maybe using an if statement when the starting point is closer than any other data point it closes the contour and starts a new one on the next point of your vector?

I think this would be a very nice addition to your program for when your working with big spatial datasets.

Cheers

Tristan Ursell

### Tristan Ursell (view profile)

Zheng -- found a small bug that was responsible. New version uploaded, and awaiting approval from Mathworks (or just add "bad_pts=zeros(size(Xin));") at line 235. Thanks!

Tristan Ursell

### Tristan Ursell (view profile)

Hi Zheng, I presume this was with the last example? I'll take a look.

Zheng

### Zheng (view profile)

This is what showed up in the prompt window:

Error in points2contour (line 276)

Tristan Ursell

### Tristan Ursell (view profile)

@Roy, I included the requested feature late last month.

Tristan Ursell

### Tristan Ursell (view profile)

@Roy -- good suggestion, I will try to implement that soon.

Roy

### Roy (view profile)

This looks great, but I have a feature request:

The built-in command 'sort' has the optional syntax

[y,i]=sort(x)

where y is a sorted version of x, and x(i)=y

This is useful if you have several data vectors that all need to be re-ordered along with x.

I think adding this feature to the program would be pretty simple, but that you could do it much faster than I could.

Thanks!

Jay Cheng

### Jay Cheng (view profile)

Awesome work! It works on my data set while built in poly2cw does not.

my x=[121
121
121
121
121
121
121
121
121
121
121
121
121
121
121
121
121
122
122
122
122
122
122
122
122
122
123
123
123
123
123
124
124
124
125
125
125
125
125
125
125
125
126
126
126
127
127
127
128
128
128
128
129
129
129
130
130
130
130
131
131
131
131
132
132
133
133
134
134
135
135
136
136
137
137
138
138
139
139
140
140
141
141
142
142
143
143
144
144
145
145
145
145
146
146
147
147
148
148
149
149
150
150
151
151
152
152
153
153
154
154
155
155
156
156
157
157
158
158
159
159
160
160
161
161
162
162
163
163
164
164
165
165
166
166
167
167
168
168
169
169
170
170
171
171
171
171
172
172
173
173
174
174
175
175
175
175
175
175
176
176
176
176
177
177
178
178
179
179
180
180
181
181
182
182
183
183
184
184
184
184
185
185
186
186
187
187
188
188
189
189
190
190
190
190
190
191
191
191
191
191
191
191
191
191
192
192
192
192
192
192
192
192
192
192
192
192
192
193
193
193
193
193
193
193
];
y=[121
121
121
121
121
121
121
121
121
121
121
121
121
121
121
121
121
122
122
122
122
122
122
122
122
122
123
123
123
123
123
124
124
124
125
125
125
125
125
125
125
125
126
126
126
127
127
127
128
128
128
128
129
129
129
130
130
130
130
131
131
131
131
132
132
133
133
134
134
135
135
136
136
137
137
138
138
139
139
140
140
141
141
142
142
143
143
144
144
145
145
145
145
146
146
147
147
148
148
149
149
150
150
151
151
152
152
153
153
154
154
155
155
156
156
157
157
158
158
159
159
160
160
161
161
162
162
163
163
164
164
165
165
166
166
167
167
168
168
169
169
170
170
171
171
171
171
172
172
173
173
174
174
175
175
175
175
175
175
176
176
176
176
177
177
178
178
179
179
180
180
181
181
182
182
183
183
184
184
184
184
185
185
186
186
187
187
188
188
189
189
190
190
190
190
190
191
191
191
191
191
191
191
191
191
192
192
192
192
192
192
192
192
192
192
192
192
192
193
193
193
193
193
193
193
];

Maybe I misunderstood poly2cw. But what ever, your points2contour can sort my coordinates into correct order! Many thanks!

Tristan Ursell

### Tristan Ursell (view profile)

Here's a partial fix of the problem you mention Catherine. There is no general solution to your problem without solving the TSP -- which is NP hard :) There is however another level of fix that would require time to write, where you sample N-points ahead in the contour and search exhaustively for the shortest local N-path. That would likely get rid of most crossings and

Catherine

### Catherine (view profile)

Yes, that is interesting, in theory. In practice, this unexpected part of my project is a very annoying bottleneck. ;) Luckily, mine is a very rough shoreline (taken off a grid), and so is not really scale-independent.

Tristan Ursell

### Tristan Ursell (view profile)

Hi Catherine --

I'll be away for a few days, but yes, I was going to add a feature like that, and will try to get to it this week. Interesting that you bring up shorelines, which are fractals, and fractals, having no absolute scale do not have a definitely shortest contour. I realize that for your example its discrete and hence there is some shortest path, but it's still really hard to connect such structures because they tend to have this scale free property.

Catherine

Catherine

### Catherine (view profile)

Tristan, do you have any hints for dealing with messy dat sets? I am tracing points along a shoreline, and your program gets it almost right, just a few places where it joins points across large distances that should not be connected. Something like 'if distance between points>X, redo' in some way or other.

Tristan Ursell

### Tristan Ursell (view profile)

@Catherine, glad to hear it's working well for you :) Please consider leaving a good rating so that others are more likely to try it.

Catherine

### Catherine (view profile)

This works great, and is very fast (on my admittedly small dataset). Thanks! I've been looking for something like this.