How is ST_PointOnSurface calculated?

Question

The PostGIS documentation states that ST_PointOnSurface returns "a POINT guaranteed to lie on the surface". It seems like this function could be trivially implemented to give results that satisfy the documentation but provide little real-world utility, though I'm certain that PostGIS provides a non-trivial implementation.

This introduction to PostGIS provides a nice comparison and contrast of ST_Centroid with ST_PointOnSurface and says that "[ST_PointOnSurface] is substantially more computationally expensive than the centroid operation".

Is there a more thorough explanation of how ST_PointOnSurface is calculated? I've been using ST_Centroid, but have encountered some edge cases in my data where the centroid is outside the geometry. I believe that ST_PointOnSurface is the correct substitute, but the function name and documentation leave room for uncertainty.

Further, is the computational expense of ST_PointOnSurface incurred even if the centroid does lie within the geometry already?

Answer 1

Based on a few experiments, I think ST_PointOnSurface() works roughly like this, if the geometry is a polygon:

1. Trace an east-west ray, lying half-way between the northern and southern extents of the polygon.
2. Find the longest segment of the ray that intersects the polygon.
3. Return the point that is half-way along said segment.

That may not make sense, so here's a sketch of a polygon with a ray dividing it into a northern and southern parts:

             _
            / \             <-- northern extent
           /   \
          /     \
         /       \
        /         \      __
       /           \    /  \
      /_ _ _ P _ _ _\  / _ _\  P = point-on-surface
     /               \/      \
    /                         \
   /            C              \   C = centroid
  /                             \
 /                              /
/______________________________/  <-- southern extent

Thus, ST_PointOnSurface() and ST_Centroid() are usually different points, even on convex polygons.

The only reason for the "surface" in the name, I think, is that if the geometry has 3D lines then the result will be simply be one of the vertexes.

I would agree that more explanation (and better naming) would have been useful and hope a GEOS programmer might shed some more light on the matter.

Answer 2

The Algorithm mentioned by Martin F is not guaranteed to work for every polygon. I am not sure if this is the algorithm used by said function or not, however consider the example below:

In this example, the longest line segment that the bisector crosses is between the points (-0.147,2.5) and (3.207,2.5), and the midpoint of this segment is the point (1.53,2.5), which is certainly not on the surface of the polygon.

Edit: I do realise this should've been a reply to Mark F's answer, however I didn't know how to include or reference an image in a reply.