PostGIS の ST_HausdorffDistance とは何か？


この記事は [FOSS4G Advent Calendar 2023](https://qiita.com/advent-calendar/2023/foss4g) 21日目の記事でございます。

## はじめに

大学時代に、PostGIS の ST_HausdorffDistance を少し触っていたので、復習がてらどういった関数なのかまとめたいと思います。

２つのジオメトリ間の距離を表現するものとして様々な計算手法があります。ハウスドルフ距離もその1つで、ユークリッド幾何学に代表される2点間の距離ではなく、2つの部分集合に対する距離を示します。私も数学はあまり深く学んでいないため、ざっくり「２つの集合に対する距離」という認識でいます。

さて、PostgreSQLの空間情報の拡張機能であるPostGISにおいても、2つのジオメトリに対してハウスドルフ距離を計測する関数が備わっています。それが [ST_HausdorffDistance](https://postgis.net/docs/manual-dev/postgis-ja.html#ST_HausdorffDistance) です。

## ハウスドルフ距離

ハウスドルフ距離は [Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%8F%E3%82%A6%E3%82%B9%E3%83%89%E3%83%AB%E3%83%95%E8%B7%9D%E9%9B%A2) を見ると少し難しそうな数式が書いてありますが、原理はそこまで難しくありません。

文章で簡潔に説明すると、ジオメトリA（線形）とジオメトリB（線形）があったとして「ジオメトリAを構成する各点からジオメトリBへの最短距離を算出したうちの最大距離」が、ハウスドルフ距離です。

図で表現してみます。例えば、以下の画像のように、複数の点で構成されている２つの線分があったとします。

![hausdorff_blank.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/93706/3065c50b-4f23-1696-2a9a-4d43234d21d1.png)

まず、各点から線分に対する最短距離を算出します。以下の画像では、線分Bの各点から、線分Aに対する最小距離を算出します（青い線）、その中で一番最大の距離がハウスドルフ距離となります（赤い線）。

![hausdorff.drawio (3).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/93706/c72cfdc2-d134-17a9-2805-b57491e6a386.png)

## ST_Distance は何を返す？

`ST_Distance` は、線形同士であっても最短距離をそのまま返します。実際に、以下の画像の線分データに対して　`ST_Distance`　を実行してみます。

![hausdorff_blank.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/93706/3065c50b-4f23-1696-2a9a-4d43234d21d1.png)

```sql
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;

```

これを実行すると、結果は２と返ってきます。

```zsh
 st_distance 
-------------
           2
```

なぜこの値が返ってくるのかというと、線分Aの (1,1) と、線分Bの　(1,3) が最短距離となり、その距離が2であるからです。

## ST_HausdorffDistance

続いて、PostGISで `ST_HausdorffDistance` を使ってみます。先程と同じく、以下の画像の線分データに対して　`ST_HausdorffDistance`　を実行してみます。

![hausdorff_blank.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/93706/3065c50b-4f23-1696-2a9a-4d43234d21d1.png)

```sql
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 と返ってきます。

```zsh
 st_hausdorffdistance 
----------------------
    4.123105625617661
```

なぜこの値が返ってくるのかというと、先程の原理の説明のように線分Aの (2,5) と、線分Bの　(3,1) の距離がハウスドルフ距離となり、その距離が 4.123105625617661 であるからです。

### 各点の最短距離を出してみる

実際にハウスドルフ距離の原理に従っているかどうか調べるため、線分Bの各点から線分Aに対する最短距離を PostGIS で算出してみます。

以下は、線分Bの各点から線分Aに対する最短距離を算出するSQL文です。地道に各点からの最短距離を `ST_Distance` を使って算出します。

```sql
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のデータベースに打つと、結果は以下のようになります。

```zsh
    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` は緑色の線の距離を算出しています。

![hausdorff.drawio.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/93706/a518df0c-5455-3638-c09c-a673cc54f915.png)

## 補足事項

[PostGIS の公式ドキュメントの ST_HausdorffDistance の説明](https://postgis.net/docs/manual-dev/postgis-ja.html#ST_HausdorffDistance) には、以下のNoteが貼られています。

> このアルゴリズムは標準的なハウスドルフ距離と等価では*ありません*。しかし、使用可能な場面の大部分で正しくなる近似計算がなされています。重要なものに、それぞれが概ね平行で概ね等しい長さのラインストリングがあります。これはラインのマッチングに使える基準です。

上記の詳しい検証まではできていない＆言葉の意味を完全に理解はできていません。文面をそのまま読むと、PostGISでは標準とは少し違うアルゴリズムでハウスドルフ距離の算出をしているようですが、値が正確になるような近似計算が施されているとのことです。

また、似たようなジオメトリ間の距離算出の関数として [ST_FrechetDistane (フレシェ距離)](https://postgis.net/docs/manual-dev/postgis-ja.html#ST_FrechetDistance) があります。これは線の位置と点の並び順を考慮して距離を算出します。

## 参考文献

* [大人になってからの再学習 ハウスドルフ距離 (Hausdorff distance)](https://zellij.hatenablog.com/entry/20111206/p1)