2014/09/03
PostGIS >  3次ベジェ曲線の近似
実質的に昨日（SVG のパスを取り込む）の続きだが、他で使うかもしれないので標題を分けた。実行環境は 2014/08/13 の実機を参照。3次ベジェ曲線を定義する四点の座標から、細かい直線で近似したジオメトリを返すストアド関数を作った。定義文は後掲。下が使用例。


上の例で近似したベジェを Adobe Illustrator で描くと ↓ のようになる（ただし縦座標の向きは PostGIS と逆）。ストアドの第一引数に p1、c1、c2、p2 の順に XY 座標を羅列して渡す。第二引数は近似する細かさのパラメータ。詳細は後述。


↓ ストアドの定義文。全体の流れは、まず p1-c1、c1-c2、c2-p2 という三つの直線を作り（WITH 句の c ブロック）最大の長さを出す（d ブロック）。この最大長と第二引数で渡されたパラメータから、三つの直線の分割位置を計算（e ブロック）。後は3次ベジェの定義どおりに曲線の軌跡となる点を作り（f〜h ブロック）、最後に一つの線ジオメトリにまとめて返す。
CREATE OR REPLACE FUNCTION "201409"."03_bezier3"(float[], float)
RETURNS geometry(LineString) LANGUAGE sql IMMUTABLE AS $BODY$

-- 第一引数（配列, 8要素）起点XY, 起点の制御点XY, 終点の制御点XY, 終点XY
-- 第二引数 区間の細かさ（最小区間長に近い数値。小さいほど精度大）

WITH a (lab, geom) AS (
	VALUES ('p1', ST_Point($1[1], $1[2]))	-- 起点XY
		, ('c1', ST_Point($1[3], $1[4]))-- 起点の制御点XY
		, ('c2', ST_Point($1[5], $1[6]))-- 終点の制御点XY
		, ('p2', ST_Point($1[7], $1[8]))-- 終点XY
), c (lab, geom) AS (
	-- 三つの直線を作成
	SELECT text 'm1', ST_MakeLine(geom) FROM a WHERE lab IN ('p1', 'c1')
	UNION ALL
	SELECT 'm2', ST_MakeLine(geom) FROM a WHERE lab IN('c1', 'c2')
	UNION ALL
	SELECT 'm3', ST_MakeLine(geom) FROM a WHERE lab IN('c2', 'p2')
), d AS (
	SELECT max(ST_Length(geom)) FROM c
), e AS (
	-- 区間の分割位置を準備
	SELECT generate_series(0, floor(max / $2) :: int) * $2 / max loc
	FROM d
	UNION SELECT 1
), f AS (
	SELECT lab, loc, ST_LineInterpolatePoint(geom, loc) geom FROM c, e
), g AS (
	SELECT loc, ST_LineInterpolatePoint(ST_MakeLine(geom), loc) geom
	FROM f WHERE lab <> 'm3' GROUP BY loc
	UNION ALL
	SELECT loc, ST_LineInterpolatePoint(ST_MakeLine(geom), loc)
	FROM f WHERE lab <> 'm1' GROUP BY loc
), h AS (
	SELECT loc, ST_LineInterpolatePoint(ST_MakeLine(geom), loc) geom
	FROM g
	GROUP BY loc ORDER BY loc
)
SELECT ST_MakeLine(geom) FROM h ;
$BODY$ ;

-- Example
SELECT ST_AsText("201409"."03_bezier3"(
	ARRAY[0, 25, 25, 0, 35, 85, 50, 25], 1)) ;

ベジェの定義、ストアドで使った主な関数の詳細は ↓ を参照。今回、特に役立った PostGIS の関数は ST_LineInterpolatePoint。線ジオメトリと 0〜1 の位置から線上の内挿点を得られる。

■ Wikipedia : ベジェ曲線
http://ja.wikipedia.org/wiki/%E3%83%99%E3%82%B8%E3%82%A7%E6%9B%B2%E7%B7%9A

■ Rui's Blog : 中学生でもわかるベジェ曲線
http://ruiueyama.tumblr.com/post/11197882224

■ PostGIS 2.1 Manual : ST_MakeLine, ST_LineInterpolatePoint
http://postgis.net/docs/manual-2.1/ST_MakeLine.html
http://postgis.net/docs/manual-2.1/ST_Line_Interpolate_Point.html

冒頭の使用例の結果を R で可視化。↓ 第二引数の細かさによる違いも見てみた。
rm(list=ls())
library(RPostgreSQL)
library(maptools)
library(rgeos)
con = dbConnect(PostgreSQL()
	, dbname='xxx', port=xxx, user='xxx', password='xxx')
df1 = dbGetQuery(con, '
	SELECT ST_AsText("201409"."03_bezier3"(
		ARRAY[0, 25, 25, 0, 35, 85, 50, 25]
		, unnest(ARRAY[2, 5])))
')
invisible(dbDisconnect(con))
if(is.null(dev.list())) windows(height=(300-46)/96, width=(444-8)/96)
par(las=1, mgp=c(25, 5, 0)/10, plt=c(10, 95, 15, 95)/100, tck=-0.02)
nrw = nrow(df1)
col = c(1, 2)
for(i in 1:nrw){
	plot(add=ifelse(i==1, F, T)
		, axes=ifelse(i==1, T, F)
		, readWKT(df1[i,1]), col=col[i])
}
legend(bty='n', col=col, leg=c(2, 5), lty=1
	, title='fineness', x='topleft')



先ほどの Illustrator 画面より縦に圧縮されているが、形状は同じで大丈夫そう。細かさパラメータが小さいほど各直線の長さも小さく高精度になる。このスケールだと 2 程度で十分。

↓ 出力された近似線の点の数を ST_NumPoints 関数で確認する例。細かさパラメータを決める時の参考になるかも。
WITH a AS (
	SELECT "201409"."03_bezier3"(
		ARRAY[0, 25, 25, 0, 35, 85, 50, 25]
		, unnest(ARRAY[0.1, 1, 5])) geom
)
SELECT ST_NumPoints(geom), ST_AsText(geoM) FROM a ;


これで一応 SVG の3次ベジェ曲線を PostGIS に取り込む準備ができた。昨日の続きに使いながら、ストアドでおかしい所があったら修正する。
<< SVG のパスを取り込む（2）
SVG のパスを取り込む（1） >>