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 :wink:

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:

image

Solving the connections here looks to me like a NP-hard problem, a special version of the travelling salesman problem. Tricky :wink:

If your input looked indeed like this:
image

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 :slight_smile:

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 :wink:

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 :slight_smile:

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 :slight_smile:

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. :slight_smile:

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