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!)
My take on smoothing - Smoothing Iterative Averaging (SIA) + ST_SimplifyPreserveTopology using PostGIS:
Keep in mind for good results: ss < si, sl = 0 if ss > 1
sl is my addition to algorithm to not change longer segments as I use it for smoothing out GPS accuracy issues in line tracks. After SIA I also use St_SimplifyPreserveTopology to remove tiny nuances and downsize number of vertices. I don't care about 3D and have only Line Strings. Shouldn't be hard to adapt if you use multiline strings and 3D.
update table_name set geom = ST_SimplifyPreserveTopology(ST_SIA(geom, 1, 3, 20), 1);
-- Smoothing Iterative Averaging (SIA) -- Parameters: -- Smoothing Sensitivity (ss) - positive integer -- Smoothing Iteration (si) - positive integer -- Smoothing Limit (sl) - 0+ numeric - maximum distance between neighbouring vertices, 0 = do not apply limit create or replace function ST_SIA(geom geometry, ss int default 1, si int default 5, sl numeric default 0) returns geometry as $$ declare p_prev geometry; pi geometry; p_next geometry; new_line geometry; s integer; i integer; i_prev integer; i_next integer; p_count integer; begin if ss < 1 or si < 1 or sl < 0 then raise exception 'Smoothing sensitivity and smoothing iteration must be positive integer! Smoothing Limit 0+!'; end if; if ST_GeometryType(geom) != 'ST_LineString' or ST_NPoints(geom) <= ss*2+1 then return geom; end if; new_line := geom; p_count := st_npoints(geom); for s in 1..si loop for i in 2..p_count-1 loop i_prev := i-ss; i_next := i+ss; if i <= ss then i_prev := 1; end if; if i+ss > p_count then i_next := p_count; end if; p_prev := ST_PointN(geom, i_prev); pi := ST_PointN(geom, i); p_next := ST_PointN(geom, i_next); if sl = 0 or ST_Length(ST_MakeLine(p_prev, p_next)) < sl then new_line := ST_SetPoint(new_line, i-1, ST_Centroid(ST_MakeLine(array[p_prev,pi,p_next]))); end if; end loop; geom := new_line; end loop; return geom; end; $$ language plpgsql;