# Simplification by constrained anti-aliasing in PostGIS

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.

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.

• OP let me know how well this does on real data. @dr_jts since I have not tested for robustness - any immediate issues you can see? Commented Apr 26, 2022 at 13:12
• Clever! There might be an issue if the "raster" polygon boundary loops back on itself and touches, but that can probably be taken care of by `ST_MakeValid`. Commented Apr 26, 2022 at 14:48
• @dr_jts yeah I mentioned that it will only work with simple Polygons - as in valid. I'll go through your JTS algorithm later, see how you handle some of the ambiguities. Commented Apr 26, 2022 at 15:12
• The JTS algorithm is doing a fair bit more work to prevent invalidity occurring (but it's more general-purpose as well). My guess is that it is too complex to consider porting it to SQL - but feel free to try! Commented Apr 26, 2022 at 17:23

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`.

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,
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),
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,
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),
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
``````