4.5

4.5 | 5 ratings Rate this file 35 Downloads (last 30 days) File Size: 3.61 KB File ID: #35488

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

by

 

06 Mar 2012 (Updated )

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

| Watch this File

File Information
Description

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

Required Products MATLAB
MATLAB release MATLAB 7.9 (R2009b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (15)
23 Aug 2014 Frank Pennekamp

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

29 Sep 2013 Tristan Ursell

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!

29 Sep 2013 Tristan Ursell

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

28 Sep 2013 Zheng

Hi, I like your idea but every time I run your examples with your points2contour function, errors occur. Please help asap!

This is what showed up in the prompt window:
Undefined function or variable "bad_pts".

Error in points2contour (line 276)
temp1=find(~bad_pts);

13 Aug 2013 Tristan Ursell

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

23 Jul 2013 Tristan Ursell

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

23 Jul 2013 Roy

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!

10 Apr 2013 Jay Cheng

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!

15 Mar 2012 Tristan Ursell

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

13 Mar 2012 Catherine

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.

13 Mar 2012 Tristan Ursell

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.

13 Mar 2012 Catherine  
13 Mar 2012 Catherine

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.

12 Mar 2012 Tristan Ursell

@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.

12 Mar 2012 Catherine

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

Updates
08 Mar 2012

sped up distance for loop

15 Mar 2012

Added new option to exclude outlying points.

24 Jul 2013

added a (requested) indexing feature

30 Sep 2013

fixed small bug

Contact us