この記事は FOSS4G Advent Calendar 2023 21日目の記事でございます。
はじめに
大学時代に、PostGIS の ST_HausdorffDistance を少し触っていたので、復習がてらどういった関数なのかまとめたいと思います。
２つのジオメトリ間の距離を表現するものとして様々な計算手法があります。ハウスドルフ距離もその1つで、ユークリッド幾何学に代表される2点間の距離ではなく、2つの部分集合に対する距離を示します。私も数学はあまり深く学んでいないため、ざっくり「２つの集合に対する距離」という認識でいます。
さて、PostgreSQLの空間情報の拡張機能であるPostGISにおいても、2つのジオメトリに対してハウスドルフ距離を計測する関数が備わっています。それが ST_HausdorffDistance です。
ハウスドルフ距離
ハウスドルフ距離は Wikipedia を見ると少し難しそうな数式が書いてありますが、原理はそこまで難しくありません。
文章で簡潔に説明すると、ジオメトリA（線形）とジオメトリB（線形）があったとして「ジオメトリAを構成する各点からジオメトリBへの最短距離を算出したうちの最大距離」が、ハウスドルフ距離です。
図で表現してみます。例えば、以下の画像のように、複数の点で構成されている２つの線分があったとします。
まず、各点から線分に対する最短距離を算出します。以下の画像では、線分Bの各点から、線分Aに対する最小距離を算出します（青い線）、その中で一番最大の距離がハウスドルフ距離となります（赤い線）。
ST_Distance は何を返す？
ST_Distance は、線形同士であっても最短距離をそのまま返します。実際に、以下の画像の線分データに対して ST_Distance を実行してみます。
SELECT
ST_Distance(geomA, geomB)
FROM (
SELECT
'LINESTRING(0 0, 1 1, 2 0, 3 1, 4 0)'::geometry AS geomA,
'LINESTRING(0 4, 1 3, 2 5, 3 3, 4 4)'::geometry AS geomB
) as t;
これを実行すると、結果は２と返ってきます。
st_distance
-------------
2
なぜこの値が返ってくるのかというと、線分Aの (1,1) と、線分Bの (1,3) が最短距離となり、その距離が2であるからです。
ST_HausdorffDistance
続いて、PostGISで ST_HausdorffDistance を使ってみます。先程と同じく、以下の画像の線分データに対して ST_HausdorffDistance を実行してみます。
SELECT
ST_HausdorffDistance(geomA, geomB)
FROM (
SELECT
'LINESTRING(0 0, 1 1, 2 0, 3 1, 4 0)'::geometry AS geomA,
'LINESTRING(0 4, 1 3, 2 5, 3 3, 4 4)'::geometry AS geomB
) as t;
これを実行すると、結果は 4.123105625617661 と返ってきます。
st_hausdorffdistance
----------------------
4.123105625617661
なぜこの値が返ってくるのかというと、先程の原理の説明のように線分Aの (2,5) と、線分Bの (3,1) の距離がハウスドルフ距離となり、その距離が 4.123105625617661 であるからです。
各点の最短距離を出してみる
実際にハウスドルフ距離の原理に従っているかどうか調べるため、線分Bの各点から線分Aに対する最短距離を PostGIS で算出してみます。
以下は、線分Bの各点から線分Aに対する最短距離を算出するSQL文です。地道に各点からの最短距離を ST_Distance を使って算出します。
SELECT
ST_Distance(geomA, pointA),
ST_Distance(geomA, pointB),
ST_Distance(geomA, pointC),
ST_Distance(geomA, pointD),
ST_Distance(geomA, pointE)
FROM (
SELECT
'LINESTRING(0 0, 1 1, 2 0, 3 1, 4 0)'::geometry AS geomA,
'POINT(0 4)'::geometry AS pointA,
'POINT(1 3)'::geometry AS pointB,
'POINT(2 5)'::geometry AS pointC,
'POINT(3 3)'::geometry AS pointD,
'POINT(4 4)'::geometry AS pointE
) as t;
上記のSQL文をPostgreSQLのデータベースに打つと、結果は以下のようになります。
st_distance | st_distance | st_distance | st_distance | st_distance
--------------------+-------------+-------------------+-------------+--------------------
3.1622776601683795 | 2 | 4.123105625617661 | 2 | 3.1622776601683795
実際にハウスドルフ距離として算出されている値は、それぞれの最短距離の最大値であることが、上記の結果からわかります。
ST_Distance と ST_HausdorffDistance の比較
以下の画像で、各関数の距離の算出方法を見てみます。 ST_Distance は赤色の線の距離を算出しています。２つの線分に対する最短距離です。 ST_HausdorffDistance は緑色の線の距離を算出しています。
補足事項
PostGIS の公式ドキュメントの ST_HausdorffDistance の説明 には、以下のNoteが貼られています。
このアルゴリズムは標準的なハウスドルフ距離と等価ではありません。しかし、使用可能な場面の大部分で正しくなる近似計算がなされています。重要なものに、それぞれが概ね平行で概ね等しい長さのラインストリングがあります。これはラインのマッチングに使える基準です。
上記の詳しい検証まではできていない＆言葉の意味を完全に理解はできていません。文面をそのまま読むと、PostGISでは標準とは少し違うアルゴリズムでハウスドルフ距離の算出をしているようですが、値が正確になるような近似計算が施されているとのことです。
また、似たようなジオメトリ間の距離算出の関数として ST_FrechetDistane (フレシェ距離) があります。これは線の位置と点の並び順を考慮して距離を算出します。