やってみる

アウトプットすべく己を導くためのブログ。その試行錯誤すらたれ流す。

SQLite3学習 Geopoly(2次元ベクタ画像の生成)

 基本図形くらいしかできない。

成果物

情報源

前提

Geopoly

Geopolyモジュールは、GeoJSON表記(RFC-7946)を使用して2次元ポリゴンを記述するR-Tree拡張の代替インターフェースです。Geopolyには、あるポリゴンが別のポリゴンに含まれているか、別のポリゴンと重なっていることを検出し、ポリゴンで囲まれた領域を計算し、ポリゴンを線形変換し、ポリゴンをSVGとしてレンダリングするなどの機能が含まれます 。

  • 単純なポリゴンを記述する頂点のJSON配列形式
  • [[X,Y],[X,Y],...]
  • 反時計回り

 たとえば以下はX軸上にある二等辺三角形

[[0,0]、[1,0]、[0.5,1]、[0,0]]

R-Tree拡張との違い

項目 Geopoly R-Tree拡張
次元 2 1〜5
多角形 長方形

 細かい形で重なり判定したいならGeopolyを使うといい。

 R-Treeは前に触った。

create,insert

 geopolyを使うには仮想テーブルを作る必要がある。これには_shape列が自動的に含まれる。

CREATE VIRTUAL TABLE mygeo0 USING geopoly();

 _shape列に頂点データを入れる。

INSERT INTO mygeo0 VALUES('[[0,0],[1,0],[0.5,1],[0,0]]');
.headers on
SELECT * FROM mygeo0;
_shape

 出力結果はバイナリのようだ。

補助列

CREATE VIRTUAL TABLE newtab USING geopoly(a,b,c);

 a,b,cは任意で作れる列。もちろん名前も数も自由。検索するためのメタデータにでもすればいいと思う。なにせ最小構成はIDすらない。

select

CREATE VIRTUAL TABLE mygeo1 USING geopoly(name);
INSERT INTO mygeo1(name,_shape) \
  VALUES('二等辺三角形','[[0,0],[1,0],[0.5,1],[0,0]]');
INSERT INTO mygeo1(name,_shape) \
  VALUES('正方形','[[10,10],[11,10],[11,11],[10,11],[10,10]]');

geopoly_overlap()

 _shapeが指定領域の一部にでも重なるならヒットする。

SELECT * FROM mygeo1 
WHERE geopoly_overlap(_shape, '[[0,0],[0.5,0],[0.25,0.5],[0,0]]');
_shape|name
|二等辺三角形
SELECT * FROM mygeo1 
WHERE geopoly_overlap(_shape, '[[10,10],[10.5,10],[10.25,10.5],[10,10]]');
_shape|name
|正方形
SELECT * FROM mygeo1 
WHERE geopoly_overlap(_shape, 
  '[[100,100],[100.5,100],[100.25,100.5],[100,100]]');



geopoly_within()

 _shapeが指定領域に完全に含まれるならヒットする。

SELECT * FROM mygeo1 
WHERE geopoly_within(_shape, '[[0,0],[0.5,0],[0.25,0.5],[0,0]]');



 同じか、それ以上だとヒットする。

SELECT * FROM mygeo1 
WHERE geopoly_within(_shape, '[[0,0],[1,0],[0.5,1],[0,0]]');
_shape|name
|二等辺三角形
SELECT * FROM mygeo1 
WHERE geopoly_within(_shape, '[[0,0],[2,0],[1,2],[0,0]]');
_shape|name
|二等辺三角形
SELECT * FROM mygeo1 
WHERE geopoly_within(_shape, '[[0,0],[1.1,0],[0.55,1.1],[0,0]]');
_shape|name
|二等辺三角形

 以下はヒットしない。Yが最も高い頂点は中央0.55ではなく0.6である。そのせいで_shapeの一部がはみ出して完全には含まれなくなってしまうため。

SELECT * FROM mygeo1 
WHERE geopoly_within(_shape, '[[0,0],[1.1,0],[0.6,1.1],[0,0]]');



使ってみる

CREATE VIRTUAL TABLE newtab USING geopoly(a,b,c);
INSERT INTO newtab(_shape) 
  VALUES('[[0,0],[1,0],[0.5,1],[0,0]]');
SELECT * FROM newtab 
  WHERE geopoly_overlap(_shape, $query_polygon);

関数一覧

関数 説明
geopoly_overlap(P1,P2) 戻り値0: P1P2に重なり無し
戻り値 非0:重なり有り
geopoly_within(P1,P2) 戻り値0: P1P2に含まれない
戻り値 非0: P1P2に含まれる
geopoly_area(P) ポリゴンPで囲まれた面積を返す
geopoly_blob(P) ポリゴンPのバイナリエンコーディングをBLOB形式で返す
geopoly_json(P) ポリゴンPのGeoJSON表現をTEXT文字列として返す
geopoly_svg(P,...) ポリゴンPSVG表現であるTEXT文字列を返す
geopoly_bbox(P) ポリゴンPを完全に囲む最小の長方形である新しいポリゴンを返す
geopoly_group_bbox(P) 集約中にみられるすべてのポリゴンPを囲む最小の長方形を返す
geopoly_contains_point(P,X,Y) 戻り値0: 座標X,YP内にない
戻り値 非0: 座標X,YP内にある
geopoly_xform(P,A,B,C,D,E,F) ポリゴンPのアフィン変換。
x1 = A * x0 + B * y0 + E
y1 = C * x0 + D * y0 + F
geopoly_regular(X,Y,R,N) 基本図形を返す。
Rは外接半径。Nは角数。3<=N<=1000N<3ならNULLを返す。1000<Nなら1000とする
geopoly_ccw(J) 反時計回り回転でポリゴンを返す

 なお、上記関数の引数PがポリゴンでないならNULLを返す。

geopoly_area(P)

 ポリゴンPで囲まれた面積を返す。

select geopoly_area('[[0,0],[1,0],[0.5,1],[0,0]]');
0.5

 テーブルの列を渡してみる。

CREATE VIRTUAL TABLE mygeo1 USING geopoly(name);
INSERT INTO mygeo1(name,_shape) 
  VALUES('二等辺三角形','[[0,0],[1,0],[0.5,1],[0,0]]');
INSERT INTO mygeo1(name,_shape) 
  VALUES('正方形','[[10,10],[11,10],[11,11],[10,11],[10,10]]');
select geopoly_area(_shape) from mygeo1 
where name='二等辺三角形';
0.5

 正方形でやってみる。

select geopoly_area(_shape) from mygeo1 
where name='正方形';
1.0

 面積の計算式。

三角形の面積=(底辺×高さ)÷2=(1*1)/2=0.5
四角形の面積=幅×高さ=1*1=1

 翻訳したら「領域」と出たが「面積」だと思う。

geopoly_blob(P)

 ポリゴンPのバイナリエンコーディングをBLOB形式で返す。

select geopoly_blob('[[0,0],[1,0],[0.5,1],[0,0]]');



 もしやselectして取得される値か。

geopoly_json(P)

 ポリゴンPのGeoJSON表現をTEXT文字列として返す。

select geopoly_json('[[0,0],[1,0],[0.5,1],[0,0]]');
[[0.0,0.0],[1.0,0.0],[0.5,1.0],[0.0,0.0]]

 小数値になった。

 バイナリ値からでも使えた。

select geopoly_json(geopoly_blob('[[0,0],[1,0],[0.5,1],[0,0]]'));
[[0.0,0.0],[1.0,0.0],[0.5,1.0],[0.0,0.0]]

geopoly_svg(P)

 ポリゴンPSVG表現であるTEXT文字列を返す。

select geopoly_svg('[[0,0],[1,0],[0.5,1],[0,0]]');
<polyline points='0,0 1,0 0.5,1 0,0'></polyline>

 BLOB形式からでも可能。

select geopoly_svg(geopoly_blob('[[0,0],[1,0],[0.5,1],[0,0]]'));
<polyline points='0,0 1,0 0.5,1 0,0'></polyline>

 任意の属性も付与できる。

select geopoly_svg('[[0,0],[1,0],[0.5,1],[0,0]]', 'class="mycls"', 'style="fill:blue"');
<polyline points='0,0 1,0 0.5,1 0,0' class="mycls" style="fill:blue"></polyline>

 ファイル出力してみる。

sqlite3 :memory: \
"SELECT '<svg width=\"480\" height=\"320\">';" \
"SELECT geopoly_svg(
  '[[100,100],[200,100],[150,200],[100,100]]', 
  'style=\"fill:blue\"');" \
"SELECT '</svg>';" > geopoly_svg.0.svg

 アフィン変換を利用して100倍にしつつ50pixcel平行移動。

sqlite3 :memory: \
"SELECT '<svg width=\"480\" height=\"320\">';" \
"SELECT geopoly_svg(
  geopoly_xform(
    '[[0,0],[1,0],[0.5,1],[0,0]]', 100, 0, 0, 100, 50, 50
  ), 'style=\"fill:blue\"');" \
"SELECT '</svg>';" > geopoly_svg.1.svg

geopoly_bbox(P)

 ポリゴンPを完全に囲む最小の長方形である新しいポリゴンを返す。

select geopoly_bbox('[[0,0],[1,0],[0.5,1],[0,0]]');



 バイナリ形式で返ってくるらしい。bboxってbinary_boxか?

select geopoly_json(geopoly_bbox('[[0,0],[1,0],[0.5,1],[0,0]]'));
[[0.0,0.0],[1.0,0.0],[1.0,1.0],[0.0,1.0],[0.0,0.0]]

 自作の三角形とそれを完全に囲む最小矩形を重ねて表示する。

geopoly_bbox.3.sql

.param set @polygon '[[0,0],[1,0],[0.5,1],[0,0]]'
SELECT '<svg width="300" height="300">';
WITH t1(svg, color) AS (VALUES
  ((select @polygon), 'red'),
  ((select geopoly_bbox(@polygon)), 'green')
)
SELECT
  geopoly_svg(
    geopoly_xform(svg, 100, 0, 0, 100, 200, 200), 
    printf('style="fill:%s;stroke:%s;stroke-width:2"',color), 
    printf('fill-opacity="0.5"')
  )
  FROM t1;
SELECT '</svg>';
sqlite3 :memory: < geopoly_bbox.3.sql > geopoly_bbox.3.svg
  • 座標を100倍にした
  • 200平行移動した
  • 半透明にして重なりが見えるようにした

geopoly_group_bbox(P)

 集約中にみられるすべてのポリゴンPを囲む最小の長方形を返す。

select geopoly_group_bbox(_shape) 
from (
  select '[[0,0],[1,0],[0.5,1],[0,0]]' as _shape 
  union 
  select '[[10,10],[11,10],[11,11],[10,11],[10,10]]'
);



 バイナリで出た。GeoJson形式にする。

select 
  geopoly_json(geopoly_group_bbox(_shape)) 
from (
  select '[[0,0],[1,0],[0.5,1],[0,0]]' as _shape 
  union 
  select '[[10,10],[11,10],[11,11],[10,11],[10,10]]'
);
[[0.0,0.0],[11.0,0.0],[11.0,11.0],[0.0,11.0],[0.0,0.0]]

 どちらも含む矩形領域が出た。

geopoly_contains_point(P,X,Y)

 座標X,YP内にないなら戻り値0。座標X,YP内にあるなら戻り値は非0

 集約中にみられるすべてのポリゴンPを囲む最小の長方形を返す。

 以下は含まれる。

select geopoly_contains_point(
  '[[0,0],[1,0],[0.5,1],[0,0]]',0.5,0.5);
2

 1じゃないのか。ふつう非0といえば1だと思うのだが。中点だから特別なのかと思ったがY0.2にしても同じだった。

 以下は含まれない。

select geopoly_contains_point(
  '[[0,0],[1,0],[0.5,1],[0,0]]',1,1);
0

geopoly_xform(P,A,B,C,D,E,F)

 ポリゴンPのアフィン変換。以下のように計算する。

x1 = A * x0 + B * y0 + E
y1 = C * x0 + D * y0 + F

平行移動

 変形せず平行移動するだけなら以下。

.param set @polygon '[[0,0],[1,0],[0.5,1],[0,0]]'
.param set @DX 5
.param set @DY 7
select geopoly_xform(@polygon, 1, 0, 0, 1, @DX, @DY);



 バイナリで出たのでGeoJsonに変換。

.param set @polygon '[[0,0],[1,0],[0.5,1],[0,0]]'
.param set @DX 5
.param set @DY 7
select geopoly_json(
  geopoly_xform(@polygon, 1, 0, 0, 1, @DX, @DY));
[[5.0,7.0],[6.0,7.0],[5.5,8.0],[5.0,7.0]]

 OK。

回転

 座標[0,0]を中心にRラジアンで回転する。

.param set @polygon '[[0,0],[1,0],[0.5,1],[0,0]]'
.param set @DX 5
.param set @DY 7
.param set @R (PI()/2)
select geopoly_json(
  geopoly_xform(@polygon, cos(@R), sin(@R), -sin(@R), cos(@R), 0.5, 0.5));
[[0.0,0.0],[6.12323e-17,-1.0],[1.0,-0.5],[0.0,0.0]]

 もし以下のように「関数がない」と怒られたら前回のでビルド&ロードする。

Error: no such function: sin

 ラジアンについてはこちらを参考にした。

ラジアン
30° π/6
45° π/4
60° π/3
90° π/2

 @R=(PI()/2)=π/2=90°のはず。つまり[0,0]軸を中心にして90°回転させたつもり。頂点データだけだとイメージしづらい……。単純に右へ90°傾けたかったのだが、なんか違うっぽい。数学わからん。

geopoly_regular(X,Y,R,N)

 基本図形を返す。

 Rは外接半径。Nは角数。3<=N<=1000N<3ならNULLを返す。1000<Nなら1000とする

 XYを中心とし、Rの外接半径をもつ、側面の角がN個の、単純・正・等辺・等角多角形を返す。

 以下、外接半径1の正三角形のつもり

select geopoly_json(geopoly_regular(0, 0, 1, 3));
[[1.00007,0.0],[-0.499953,0.866088],[-0.499953,-0.866088],[1.00007,0.0]]

 3〜20角までの基本図形が並んだSVGを出力する。

geopoly_regular.0.sql

SELECT '<svg width="600" height="300">';
WITH t1(x,y,n,color) AS (VALUES
   (100,100,3,'red'),
   (200,100,4,'orange'),
   (300,100,5,'green'),
   (400,100,6,'blue'),
   (500,100,7,'purple'),
   (100,200,8,'red'),
   (200,200,10,'orange'),
   (300,200,12,'green'),
   (400,200,16,'blue'),
   (500,200,20,'purple')
)
SELECT
   geopoly_svg(geopoly_regular(x,y,40,n),
        printf('style="fill:none;stroke:%s;stroke-width:2"',color))
   || printf(' <text x="%d" y="%d" alignment-baseline="central" text-anchor="middle">%d</text>',x,y+6,n)
  FROM t1;
SELECT '</svg>';
sqlite3 :memory: < geopoly_regular.0.sql > geopoly_regular.0.svg

3 4 5 6 7 8 10 12 16 20

geopoly_ccw(J)

 反時計回り回転でポリゴンを返す。

.param set @polygon '[[0,0],[1,0],[0.5,1],[0,0]]'
select geopoly_json(geopoly_ccw(@polygon));
[[0.0,0.0],[1.0,0.0],[0.5,1.0],[0.0,0.0]]

 元から反時計回りなので変化なし。

 では時計回りなら?

.param set @polygon_r '[[1,0],[0,0],[0.5,1],[1,0]]'
select geopoly_json(geopoly_ccw(@polygon_r));
[[1.0,0.0],[0.5,1.0],[0.0,0.0],[1.0,0.0]]

 反時計回りに変換された。

 GeoJsonにおいては反時計回りが正しいようなので、この関数は便利そう。

対象環境

$ uname -a
Linux raspberrypi 4.19.42-v7+ #1218 SMP Tue May 14 00:48:17 BST 2019 armv7l GNU/Linux

前回まで