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 nearestneighbor path in either the 'cw' or 'ccw' directions. The code has been written to handle square and hexagon grid points, as well as any nongrid 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 counterclockwise, 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 nonsymmetric. 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 NPHard 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*(11/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:end1),Yout(2:end1),'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:end1),Yout(2:end1),'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*(11/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:end1),Yout(2:end1),'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*(11/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:end1),Yout(2:end1),'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
1.5.0.0  fixed small bug 

1.4.0.0  added a (requested) indexing feature 

1.3.0.0  Added new option to exclude outlying points. 

1.1.0.0  sped up distance for loop 
Create scripts with code, output, and formatted text in a single executable document.
Ting Wang (view profile)
Fabricio Castro (view profile)
Fabricio Castro (view profile)
The same is achieved using the following short function wrote by myself. He's the code:
function [x,y]=points2contour(X,Y)
coord = cat(2,X,Y);
coord = unique(coord,'rows');
x = coord(:,1);
y = coord(:,2);
N=length(x);
x_=x(1,1); y_=y(1,1);
C = [];
coord_=coord;
for i=1:N1
id=any(coord_(:,1:2)~=[x_,y_],2);
coord_=coord_(id,1:2);
IDX = knnsearch(coord_,cat(2,x_,y_));
temp_=coord_(IDX,:);
x_=temp_(1,1);
y_=temp_(1,2);
C = [C;coord_(IDX,:)];
end
X=C(:,1); x=cat(1,x(1,1),X);
Y=C(:,2); y=cat(1,y(1,1),Y);
end
The above function it's more modest than yours, once it sort the random scattered points just in one direction, but it works quite good.
Christian Krogh (view profile)
Almog (view profile)
Bing Liu (view profile)
@Tristan Ursell thanks for your reply, I understand it now :)
I had not carefully test "the bug" I mentioned in the previous post, If I find time, I will check it out.
Thanks for your reply and wonderful work again. :)
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 (view profile)
bad case for this function:
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 (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 (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 (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 (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 (view profile)
Hi Zheng, I presume this was with the last example? I'll take a look.
Zheng (view profile)
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);
Tristan Ursell (view profile)
@Roy, I included the requested feature late last month.
Tristan Ursell (view profile)
@Roy  good suggestion, I will try to implement that soon.
Roy (view profile)
This looks great, but I have a feature request:
The builtin 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 reordered 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 (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 (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 Npoints ahead in the contour and search exhaustively for the shortest local Npath. That would likely get rid of most crossings and
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 scaleindependent.
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 (view profile)
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 (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 (view profile)
This works great, and is very fast (on my admittedly small dataset). Thanks! I've been looking for something like this.