Simplification by constrained anti-aliasing in PostGIS

Question

I vectorized a raster data set, therefore, I'm facing a "steps" effect, or aliasing, due to the fact that the vectorization strictly follows the raster cells. So, the geometry is quite large (several dozen MB).

I'm wondering; what would be the best and most efficient function(s) to use within PostGIS to operate an anti-aliasing like simplification which follows these simple rules:

• the original shape should never be cut by its anti-aliased version,
• the anti-aliased version should be as close as possible to the original shape (i.e. the differences should be globally as minimal as possible)

The ideal would be to draw an outline joining all the "summits" of the steps, I would say.

The best I can come with today is this blue outline:

which is the result of such query:

SELECT
  ST_Buffer(
    ST_SimplifyPreserveTopology(
      ST_GeomFromWKB(geom, {SRS}),
      20),
  10)

But that means I have to be "aware" of the original raster cell size (20 meters) and the result is far from being perfect because this query sometimes cuts the original shape (in the areas depicted by the red ellipses).

A convex hull is not an option because the shape is globally blobby, at a lower frequency than the "steps", so it has some local concave curvatures portions.

I'm using PG14 with PostGIS 3.2 on Ubuntu 22.04.

Answer 1

Let's see if we can simulate @dr_jts Concave Polygon Hull algorithm within PostGIS - I haven't actually read the JTS algorithm, and have not tested this for 100% robustness, but it seems to work quite well...

---

The core idea is to check the (vector) direction of a given Point relative to the vector between its pre- and successor and the enclosing Polygon area, to determine it's concavity predicate. We do this by finding the determinant sign of the cross product between the preceding and succeeding Points, and the test Point; the sign then tells us if the so formed triangle is concave or convex with respect to the enclosed Polygon area.

First, let's create a utility function to calculate the cross product:

CREATE FUNCTION _ST_CrossProduct (
  pt GEOMETRY(POINT),
  v1 GEOMETRY(POINT),
  v2 GEOMETRY(POINT)
) RETURNS FLOAT LANGUAGE 'plpgsql' AS
  $BODY$
    DECLARE
        x1 FLOAT := ST_X($1); x2 FLOAT := ST_X($2); x3 FLOAT := ST_X($3);
        y1 FLOAT := ST_Y($1); y2 FLOAT := ST_Y($2); y3 FLOAT := ST_Y($3);
        
    BEGIN
        RETURN ( x3-x2 ) * ( y1-y2 ) - ( y3-y2 ) * ( x1-x2 );
    END;
  $BODY$
  IMMUTABLE STRICT PARALLEL SAFE
;

Then, we want to feed it the Point sequence of the Polygon boundary; this is an example of how to realize that most performantly:

SELECT  _ST_CrossProduct(
            COALESCE(
                ST_PointN(bdy, i-1),
                ST_PointN(bdy, -2)
            ),
            ST_PointN(bdy, i),
            ST_PointN(bdy, i+1)
        ) AS cross_product,
        i,
        ST_PointN(bdy, i) AS geom
FROM    ply,
        LATERAL ST_Boundary(ply.geom) AS bdy,
        LATERAL Generate_Series(1, ST_NPoints(bdy)) AS i
;

We use COALESCE to have the first Point in the sequence get compared to its respective pre- and successors.

Now, if the cross_product is

• negative, this Point is to the left of the vector between its pre- and successor
• positive, this Point is to the right of the vector between its pre- and successor
• 0, this Point is co-linear with its pre- and successor

It is important to put this in perspective of the vertex order of your Polygon - for this example I assume a clockwise order!

---

With the above in mind, run this to get an estimated Polygon Concave Hull for your Polygon:

SELECT ST_MakePolygon(ST_AddPoint(geom, ST_PointN(geom, 1))) AS geom
FROM   (
    SELECT ST_MakeLine(ST_PointN(bdy, i) ORDER BY i) AS geom
    FROM   ply,
           LATERAL ST_Boundary(ply.geom) AS bdy,
           LATERAL Generate_Series(1, ST_NPoints(bdy)) AS i
    WHERE  _ST_CrossProduct(
                COALESCE(
                    ST_PointN(bdy, i-1),
                    ST_PointN(bdy, -2)
                ),
                ST_PointN(bdy, i),
                ST_PointN(bdy, i+1)
            ) <= 0.0       -- <= 0 for an OUTER hull; >= 0 for an INNER hull
) q
;

You can ever so slightly go beyond 0.0 (e.g. <= 0.0001 or >= -0.0005) to include almost co-linear Points of the opposite concavity! The below image is with <= 0.0:

[Grey area]: reference Polygon | [Green line]: outer concave hull boundary | [Blue line]: inner concave hull boundary | [Black dots]: Polygon vertices

---

As it stands, this will only work with simple Polygon geometries, having no interior rings! For anything more sophisticated, this should be packed into a more efficient PL/pgSQL function.

If I find the time, I'll write up a function and add a link to the Github Gist.

Answer 2

One way to do this is to compute the Outer Concave Polygon Hull of the vectorized polygon. Here's an example of how it can eliminate the "jaggies" (with a suitable choice of control parameter value):

This has just been added to the JTS Topology Suite. Hopefully it will b be ported to GEOS and PostGIS as well sometime this year. In the meantime, perhaps you can find a way do the processing using JTS externally to the database.

UPDATE: This will be available in PostGIS 3.3 as ST_SimplifyPolygonHull.

Answer 3

So, if I understand the question correctly, run this simple customizable SQL-code:

WITH tbla AS (SELECT (ST_Dump(geom)).geom FROM <polygon>),
     tblb AS (SELECT ST_Boundary(geom) geom FROM tbla),
     tblc AS (SELECT (ST_DumpPoints(geom)).geom FROM tblb),
     tbld AS (SELECT ST_Buffer(ST_Buffer(geom, 0.01),-0.01) geom FROM tbla),
     tble AS (SELECT ST_MakeLine(a.geom) geom FROM tblc a WHERE NOT EXISTS (SELECT 1 FROM tbld b WHERE ST_Intersects(a.geom, b.geom)))
     SELECT ST_MakePolygon(ST_AddPoint(geom, ST_StartPoint(geom))) geom FROM tble

The input is the toothed polygon (EPSG:4326) shown in the figure below, with no holes in one instance.

Figure 1.

The output is a smoothed polygon by outer points, see figure 2 below:

Figure 2

Create a spatial SQL function

CREATE OR REPLACE FUNCTION ST_SmoothedToothedPolygonExternal(
geom GEOMETRY,
radius float)
RETURNS GEOMETRY AS
$BODY$
WITH tbla AS (SELECT (ST_Dump(geom)).geom),
     tblb AS (SELECT ST_Boundary(geom) geom FROM tbla),
     tblc AS (SELECT (ST_DumpPoints(geom)).geom FROM tblb),
     tbld AS (SELECT ST_Buffer(ST_Buffer(geom, radius),-radius) geom FROM tbla),
     tble AS (SELECT ST_MakeLine(a.geom) geom FROM tblc a WHERE NOT EXISTS (SELECT 1 FROM tbld b WHERE ST_Intersects(a.geom, b.geom)))
     SELECT ST_MakePolygon(ST_AddPoint(geom, ST_StartPoint(geom))) geom FROM tble
$BODY$
LANGUAGE SQL;

RUN

SELECT ST_SmoothedToothedPolygonExternal(geom, 0.01) geom FROM poly_grid_raster

OR

CREATE OR REPLACE FUNCTION ST_SmoothedToothedPolygonInternal(
geom GEOMETRY,
radius float)
RETURNS GEOMETRY AS
$BODY$
WITH tbla AS (SELECT (ST_Dump(geom)).geom),
     tblb AS (SELECT ST_Boundary(geom) geom FROM tbla),
     tblc AS (SELECT (ST_DumpPoints(geom)).geom FROM tblb),
     tbld AS (SELECT ST_Buffer(ST_Buffer(geom, radius),-radius) geom FROM tbla),
     tble AS (SELECT ST_MakeLine(a.geom) geom FROM tblc a WHERE EXISTS (SELECT 1 FROM tbld b WHERE ST_Intersects(a.geom, b.geom)))
     SELECT ST_MakePolygon(ST_AddPoint(geom, ST_StartPoint(geom))) geom FROM tble
$BODY$
LANGUAGE SQL;

RUN

SELECT ST_SmoothedToothedPolygonInternal(geom, 0.01) geom FROM poly_grid_raster

If necessary, adjust the value of radius ST_Buffer()...

Original spatial solutions :-)