Can anyone provide an example SQL for smoothing linestrings from postgis table using Bezier curves or Iterative Averaging (SIA) algorithm?
I created a small, naive script which converts input LineStrings to CompoundCurves based on some heuristics.
What it does:
- Cuts down sharp corners to create a visually more appealing results than the original data.
- Uses plpgsql. No additional extensions required.
- Accepts an optional "smoothing factor" between 0 and 100 besides a geometry.
What it doesn't do:
- Processes MultiLineStrings. For any other geometry type, it simply returns the input.
- Uses Z and M values. It simply drops them. Use this only for 2D cartographic purposes.
- Creates mathematically correct results. The results are far from correct, and might be even visually unaesthetic in some cases (e.g. sharp corners). I didn't test it thoroughly. Always review the results!
- Runs fast. I'm sure it can be rewritten to a far more optimal form.
- Does real smoothing. There are much better algorithms (e.g. Chaiken or the ones mentioned in the question) to use for real smoothing. This answer targets people like me searching for a pure PostGIS approach automatically creating some kind of curved lines from real data.
CREATE OR REPLACE FUNCTION CreateCurve(geom geometry, percent int DEFAULT 40) RETURNS geometry AS $$ DECLARE result text; p0 geometry; p1 geometry; p2 geometry; intp geometry; tempp geometry; geomtype text := ST_GeometryType(geom); factor double precision := percent::double precision / 200; i integer; BEGIN IF percent < 0 OR percent > 100 THEN RAISE EXCEPTION 'Smoothing factor must be between 0 and 100'; END IF; IF geomtype != 'ST_LineString' OR factor = 0 THEN RETURN geom; END IF; result := 'COMPOUNDCURVE(('; p0 := ST_PointN(geom, 1); IF ST_NPoints(geom) = 2 THEN p1:= ST_PointN(geom, 2); result := result || ST_X(p0) || ' ' || ST_Y(p0) || ',' || ST_X(p1) || ' ' || ST_Y(p1) || '))'; ELSE FOR i IN 2..(ST_NPoints(geom) - 1) LOOP p1 := ST_PointN(geom, i); p2 := ST_PointN(geom, i + 1); result := result || ST_X(p0) || ' ' || ST_Y(p0) || ','; tempp := ST_Line_Interpolate_Point(ST_MakeLine(p1, p0), factor); p0 := ST_Line_Interpolate_Point(ST_MakeLine(p1, p2), factor); intp := ST_Line_Interpolate_Point( ST_MakeLine( ST_Line_Interpolate_Point(ST_MakeLine(p0, p1), 0.5), ST_Line_Interpolate_Point(ST_MakeLine(tempp, p1), 0.5) ), 0.5); result := result || ST_X(tempp) || ' ' || ST_Y(tempp) || '),CIRCULARSTRING(' || ST_X(tempp) || ' ' || ST_Y(tempp) || ',' || ST_X(intp) || ' ' || ST_Y(intp) || ',' || ST_X(p0) || ' ' || ST_Y(p0) || '),('; END LOOP; result := result || ST_X(p0) || ' ' || ST_Y(p0) || ',' || ST_X(p2) || ' ' || ST_Y(p2) || '))'; END IF; RETURN ST_SetSRID(result::geometry, ST_SRID(geom)); END; $$ LANGUAGE 'plpgsql' IMMUTABLE;
As it returns curves in a geometry type, if you would like to use it in a GIS such as QGIS, you have to wrap it into PostGIS functions converting them around. The intended use syntax is:
SELECT ST_AsText(ST_CurveToLine(CreateCurve(geom))) AS geom FROM linestringtable;
This is still an open issue in PostGIS (and other GIS tools) as stated in the book "PostGIS in Action" in chapter 2.2.6 "Curved geometries".
Here are some references to algorithms and code:
- http://www.antigrain.com/research/bezier_interpolation/index.html (also used by LineSmoother function in JTS!)