geometry - How is ST_PointOnSurface calculated?

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

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.