 # Ordering points that are close around a non-convex polygon

Hello all.

I have a list of 2D points in random order (typically 100s).

These points correspond to the pixel coordinates of the boundary of an object. They are in ‘raster’ order. If I plot them I get this:

I am trying to order them so that if I plot a line between successive pairs, them I will have a closed polygon.

There is an algorithm that works well if the polygon is convex: You get a point inside the polygon, and order points by the angle of the line that join the center to them.

But as you can see, my contours are sometimes not convex. I have read a bit on the subject, and in summary as soon as the polygons are not convex, you cannot get a unique polygon with this algorithm. If I run it on this set of points using their centroid, I get this:

It’s not great…

So my question is:
There is no algorithm that returns a unique “good” polygon for non-convex polygons, but do you guys know an algorithm that can handle such polygons when we know that each points are close to each other (they are on the pixel grid and tangent)?

Right now we are recreating a pseudo image with these coordinates and filling the polygon before extracting the bounds of the BW patch obtained. But I am looking for an algorithm that operates directly on the XY coordinates.

4 Likes

Hey @tinevez,

can you share the shown example XY coordinate list? I think I have an algorithm available but would like to test first Cheers,
Robert

Here you are:

``````        1302        1078
1303        1074
1303        1080
1304        1072
1304        1081
1305        1070
1305        1081
1306        1068
1306        1085
1307        1067
1307        1085
1308        1066
1308        1085
1309        1066
1309        1085
1310        1065
1310        1085
1311        1066
1311        1085
1312        1066
1312        1085
1313        1066
1313        1085
1314        1067
1314        1086
1315        1067
1315        1086
1316        1067
1316        1087
1317        1067
1317        1087
1318        1067
1318        1088
1319        1066
1319        1089
1320        1066
1320        1090
1321        1065
1321        1090
1322        1065
1322        1089
1323        1064
1323        1089
1324        1063
1324        1089
1325        1063
1325        1089
1326        1063
1326        1089
1327        1063
1327        1089
1328        1063
1328        1090
1329        1063
1329        1090
1330        1063
1330        1091
1331        1064
1331        1092
1332        1066
1332        1093
1333        1068
1333        1094
1334        1072
1334        1100
1335        1074
1335        1101
1336        1085
1336        1102
1337        1088
1337        1102
1338        1090
1338        1102
1339        1093
1339        1102
1340        1095
1340        1102
1341        1096
1341        1102
1342        1098
1342        1101
``````
1 Like

Hmmm, there is a discrepancy between the image you showed first and the coordinate list. If I mark the pixels in the list white, it doesn’t fit 100%ly: Solving the connections here looks to me like a NP-hard problem, a special version of the travelling salesman problem. Tricky If your input looked indeed like this: You could connect all pixel coordinates, which are less than 2 pixels apart. This has some limitations, for example that connected pixels are not allowed in an L-shape as shown here:

But maybe it solves your initally mentioned problem. You can try it out with an adaped input-table:
jy_polygon2.csv (1.1 KB)

And this CLIJ macro:

``````run("Close All");
print("\\Clear");

open("C:/structure/data/jy_polygon2.csv");
Table.rename("jy_polygon2.csv", "Results");

// init GPU
run("CLIJ2 Macro Extensions", "cl_device=");
Ext.CLIJ2_clear();

// import coordinate table, and bring it in the right format
// (first row: X, second row Y coordinates)
pointlist_cols = "pointlist_cols";
Ext.CLIJ2_resultsTableToImage2D(pointlist_cols);
Ext.CLIJ2_transposeXY(pointlist_cols, pointlist);

// determine the distance from any point to any point
Ext.CLIJ2_generateDistanceMatrix(pointlist, pointlist, distance_matrix);

// draw a mesh; scaled to make it look nicer
scaling_factor = 10;
distance_threshold = 2;
Ext.CLIJ2_create2D(mesh, 500, 500, 32);
Ext.CLIJ2_multiplyImageAndScalar(pointlist, scaled_pointlist, scaling_factor);
Ext.CLIJ2_distanceMatrixToMesh(scaled_pointlist, distance_matrix, mesh, distance_threshold);
Ext.CLIJ2_pull(mesh);
``````

I know, changing the input data is a bit like cheating, but maybe it helps anyway.

Thanks for the riddle! You made my day Cheers,
Robert

True! My bad! I sent you the data after I reduced the number of points.
Do you care getting the actual points?
In the meantime H. Gluender pointed me to a paper I adapted so that I don’t have to create an image. Maybe we could compare?

1 Like

My macro also loads a table and needs no input image Sure! Curious!

First here is the complete list of points:

``````
1302   1078
1303   1074
1303   1075
1303   1076
1303   1077
1303   1079
1303   1080
1304   1072
1304   1073
1304   1081
1305   1070
1305   1071
1305   1081
1306   1068
1306   1069
1306   1082
1306   1083
1306   1084
1306   1085
1307   1067
1307   1085
1308   1066
1308   1085
1309   1066
1309   1085
1310   1065
1310   1085
1311   1066
1311   1085
1312   1066
1312   1085
1313   1066
1313   1085
1314   1067
1314   1086
1315   1067
1315   1086
1316   1067
1316   1087
1317   1067
1317   1087
1318   1067
1318   1088
1319   1066
1319   1089
1320   1066
1320   1090
1321   1065
1321   1090
1322   1065
1322   1089
1323   1064
1323   1089
1324   1063
1324   1089
1325   1063
1325   1089
1326   1063
1326   1089
1327   1063
1327   1089
1328   1063
1328   1090
1329   1063
1329   1090
1330   1063
1330   1091
1331   1064
1331   1065
1331   1092
1332   1066
1332   1067
1332   1093
1333   1068
1333   1069
1333   1070
1333   1071
1333   1094
1334   1072
1334   1073
1334   1095
1334   1100
1335   1074
1335   1075
1335   1076
1335   1077
1335   1078
1335   1079
1335   1080
1335   1081
1335   1082
1335   1083
1335   1084
1335   1096
1335   1097
1335   1098
1335   1099
1335   1101
1336   1085
1336   1086
1336   1087
1336   1102
1337   1088
1337   1089
1337   1102
1338   1090
1338   1091
1338   1092
1338   1102
1339   1093
1339   1094
1339   1102
1340   1095
1340   1102
1341   1096
1341   1097
1341   1102
1342   1098
1342   1099
1342   1100
1342   1101
``````

just in case.

1 Like

Second, the paper H. Gluender pointed me to is:

Chang F., Chen C.-J. and Lu C.-J. (2004)
“A linear-time component-labeling algorithm using contour tracing technique.”
Computer Vision and Image Understanding 93: 206-220. doi:10.1016/j.cviu.2003.09.002

The Java-code is available from Burger & Burge:
https://imagingbook.com/source/

In my case I did not want to create an image and directly work on the points themselves.
This changed the way the algo above is used. It turns into a simple routine where one look for the closest points (in Eucidean distance) and remove it from the list.

In MATLAB, my first implementation is the following:

``````function P2 = find_countour( P )
%FIND_COUNTOUR Sort x,y coordinates along a continuous contour.

P2 = NaN( size( P ) );
current = 1;

% Start point
ps = P( 1, : );
set_visited( 1 )
P2( current, : ) = ps;
current = current + 1;

% Loop.
done = size( P, 1 ) == 1;
while ~done

done = current == size( P, 1 );

id = find_next_point( ps );
ps = P( id, : );
set_visited( id )
P2( current, : ) = ps;
current = current + 1;

end

function set_visited( id )
P( id, : ) = NaN;
end

function id = find_next_point( p0 )
[ ~, id ] = pdist2( P, p0, 'Euclidean', 'Smallest', 1 );
end

end
``````

But it is not super fast.

To treat 3000 polygons of about 50 - 200 points it takes 45s.
That is 0.322 ms / point.

2 Likes

The clij macro above needs 116 ms - including loading the table from disc and showing the mesh on screen. 18 ms without that But yes, you need to reshape your input data - can you get rid of the L-shaped pixels before saving the points in the table?

So that’s the problem of reconstructing a surface (or a curve) form a point cloud ?

I’ve read a bit about this problem because Limeseg has to do this at some point. There’s a review here (https://hal.inria.fr/hal-01348404v2/document), and the term that comes to my mind is the so called ‘Poisson surface reconstruction’ where you turn this problem into a Poisson formulation (https://dl.acm.org/doi/abs/10.5555/1281957.1281965 ). It was too complicated for my use case, so I did not follow through. And I don’t know if there’s a java implementation of it. But maybe that’s something to look at ?

3 Likes

That is 0.1488 ms / point.
You beat me by a factor of 3.
Probably the factor would be greater if I gave you the 3000 contours to analyze at once.

My limitation is that I must stay in pure MATLAB 1 Like

How about CLATLAB?

No I think the difference is that in our case we know that we have points exactly on the boundary already.
Points cloud “stain” the whole volume.

Let me try. Happy about feedback!

Then maybe look at powercrust ? https://www.cs.ucdavis.edu/~amenta/pubs/sm.pdf

But yeah then you don’t take advantage of the fact that each points are close to each other, so it’s probably overkill.

2 Likes

Ok, can you help me convert the macro to MATLAB?
Apparently I miss plenty of arguments during the function calls.

1 Like

Ok I have a new version that makes just one call to `pdist` per polygon.
And I am down now to 0.0410 ms per points ( x7.8 speed-up ).
But that is the same approach:

``````function P2 = find_countour( P )
%FIND_COUNTOUR Sort x,y coordinates along a continuous contour.

P2 = NaN( size( P ) );
current = 1;

D = squareform( pdist ( P ) );
D( D == 0 ) = Inf;

% Start point
id = 1;
ps = P( id, : );
P2( current, : ) = ps;
current = current + 1;

% Loop.
done = size( P, 1 ) == 1;
while ~done

done = current == size( P, 1 );

prev_id = id;
id = find_next_point( id );
ps = P( id, : );
set_visited( prev_id )
P2( current, : ) = ps;
current = current + 1;

end

function set_visited( id )
D( id, : ) = Inf;
D( :, id ) = Inf;
end

function id = find_next_point( id )
[ ~, id ] = min( D( :, id ) );
end

end

``````
2 Likes

Which function? The documentation lists arguments for macro and matlab - and it should be very close to identical:

1 Like

Extra factor of 2 by manually inlining functions:

``````function P2 = find_countour( P )
%FIND_COUNTOUR Sort x,y coordinates along a continuous contour.

P2 = NaN( size( P ) );
current = 1;

D = squareform( pdist ( P ) );
D( D == 0 ) = Inf;

% Start point
id = 1;
ps = P( id, : );
P2( current, : ) = ps;
current = current + 1;

% Loop.
done = size( P, 1 ) == 1;
while ~done

done = current == size( P, 1 );

prev_id = id;
[ ~, id ] = min( D( :, id ) );

ps = P( id, : );
P2( current, : ) = ps;
current = current + 1;

D( prev_id, : ) = Inf;
D( :, prev_id ) = Inf;

end
end

``````
1 Like