7

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:

anti-aliasing in PostGIS

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.

3 Answers 3

4

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):

enter image description here

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.

4

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:

Polygon Concave Hull [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.

4
  • 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?
    – geozelot
    Apr 26 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.
    – dr_jts
    Apr 26 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.
    – geozelot
    Apr 26 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!
    – dr_jts
    Apr 26 at 17:23
1

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.

enter image description here

Figure 1.

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

enter image description here

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 :-)

2
  • 1
    Interesting. The original outer corner points don't intersect the reversed buffer mainly because of the line segment approximation to the first buffer's curve around those corners?
    – hgb
    Apr 28 at 3:28
  • @hgb, All right, the first key point is to create a double buffer that allows us to select/unselect external/internal points and create smoothed polygonal objects, the second key point is a good behavior of the aggregate function ST_MakeLine()... :-)... Apr 29 at 19:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.