Better “centerpoint” than centroid of a concave polygon

I'm using the centroid of polygons to attach a marker in a map application. This works definitely fine for convex polygons and quite good for many concave polygons.

However, some polygons (banana, donut) obviously don't produce the desired result: The centroid is in these cases outside the polygons area.

Does anybody know a better approach to find a suitable point within any polygons area (which may contain holes!) to attach a marker?

Thank you!

[Matt J]

Didn't we already solve this problem last month?

https://www.mathworks.com/matlabcentral/answers/510329-how-to-find-a-figure-s-centroid-inside-a-polyshape

[David Franco]

Hi Matt,

Yes, we did it! But for my 399 polyshapes it's consuming too much time and I have to do this many times...

(each polyshape has more then 30,000 points!)

[Image Analyst]

So, how much time exactly? Hours?

[David Franco]

I tried today and it took about 8 hours to show me the result... And many numbers are still outside the borders...

[Image Analyst]

It would be faster to convert it to an image, get the centroids, and then translate those back to your analytical graph. If it were an image, it would take seconds. For each polyshape, convert it to an image with poly2mask(). Then just find the centroid with regionprops(), find the boundaries with bwboundaries, then see if the centroid is inside the region, and if it's not, reassign it to the nearest boundary point. Pretty easy. Attach your polyshapes with a .mat file if you still need help. I haven't worked with them before but I'm assuming it's easy to get each set of boundary coordinates in an N-by-2 list of (x,y) coordinates, right?

[David Franco]

Thanks Image Analyst! The polyshapes are here: Google Drive

[Image Analyst]

I tried to download it and it locked up my browser. I'm not going to try again. Just attach it here if you can instead of on a third party website.

[David Franco]

The file has 13MB. The maximum size allowed here is 5 MB.

[Image Analyst]

Can you trim it down to a subset of shapes?

[David Franco]

Done!

Answer 1

[darova]

Here is the simplest solution for this task

x = rand(4,1);
y = rand(4,1);
x = [x; x(1)];
y = [y; y(1)];

ii = 2;   % corner to place point
n = 3;    % how far inside from corner

x0 = (x(ii-1)+2*n*x(ii)+x(ii+1))/(n+1)/2;
y0 = (y(ii-1)+2*n*y(ii)+y(ii+1))/(n+1)/2;
if ~inpolygon(x0,y0,x,y)
    x0 = x(ii) - (x0-x(ii));
    y0 = y(ii) - (y0-y(ii));
end
plot([x(ii) x0],[y(ii) y0],'o-r')
line(x,y)
axis equal

Answer 2

[Chad Greene]

I just ran into this problem when trying to place a text label in the middle of a crescent-shaped ice shelf. The centroid or the mean or median of the coordinates of the ice shelf polygon are all outside the bounds of the ice shelf. Here's the best solution I could come up with:

% Convert the outline to a polyshape:
P = polyshape(x,y); 

% And get the delaunay triangulation of the polygon:  
T = triangulation(P);

% Now find the center points of all the triangles:
[C,r] = circumcenter(T);

% Get the index of the centerpoint that has the largest radius from the boundary:
[~,ind] = max(r); 

% These center coordinates are in the center of the fattest part of the polygon:  
xc = C(ind,1); 
yc = C(ind,2); 

Answer 3

[Image Analyst]

Well if you compute the distance transform of the shapes you'll get a "spine" that runs along the midline of the shape. I don't know if any point along there is better than any other point. Maybe you can just pick the point half way from one end to the other, if there even ARE endpoints. The only way I know how to do it is with digital images, not analytically with a set of (x,y) vertex points. And it requires the Image Processing Toolbox so you can use bwskel() or bwdist().

[David Franco]

Thanks, but I am working with shapefiles and I need to find the skeleton of a polygon...

Answer 4

[Chad Greene]

I added a polycenter function to do exactly what you need in the Climate Data Toolbox for Matlab.

If you have a shapefile, the sytax is:

S = shaperead('myshapefile.shp'); 
[xc,yc] = polycenter(S); 

Then xc,yc are the centerpoints of any entries in the shapefile structure S.

[David Franco]

Thank you @Chad Greene!

[Sim]

Hi @Chad Greene, cool tool!

However, I got the following error:

>> [xc,yc] = polycenter(S);
Reference to non-existent field 'Lon'.

Error in polycenter (line 116)
         x = S(k).Lon;     

Any idea on how to solve it ?

[Chad Greene]

@Sim Check the field names of the shapefile and then use x,y coordinates as inputs, like this

[xc,yc] = polycenter(x,y);

which would be the same as

[lonc,latc] = polycenter(lon,lat);

if your input coordinates are geo.