🐢

[F#]亀の甲割りの怠惰な実装

2021/05/04に公開

SVG

幾何計算

tech

私のかかわっている建築構造の世界では、「亀の甲割り」という用語があります。

具体的には、「あるポリゴンの中の任意の位置について、ポリゴンを構成する辺のうち最も近い辺に従属させるとして、ポリゴンを各辺の支配領域に区切る」といった問題になります。

計算としては、各辺に対する2等分線を引いていったりすることになるのですが、ポリゴンが複雑だと結構悩みます。

そこで今回は怠惰に、ポリゴンをメッシュ分割して各辺との距離を算出することで領域分けを行うアルゴリズムを作ってみました。

処理手順

以下のようなポリゴンを考えます。

これを包絡するようなメッシュを作ります。

メッシュの中心がポリゴン内に入っているもののみ対象とします。

それぞれのメッシュとポリゴンの各辺の距離を計算し、最も小さい距離となる辺と同じ色で着色して領域分けします。

結果として幾何学的な領域分割がやや力業で求まりました。

解説

1ステップずつ追っていけば大したものではないのですが、数学の知識をだいぶ忘れていたせいで苦戦してしまいました。

以下のアルゴリズムが書ければ簡単に実装できます。

ある点がポリゴンに内包されるかの判定
点と直線の距離

1については、Crossing Number Algorithm という方法を用いました。以下が参考になりました。

2はおそらく中学か高校の数学の問題になるので難しくはないのですが、私はほとんど記憶に残っておらず苦戦しました。
最終的には以下を参考にしました。

まずはあるメッシュがポリゴン内に含まれるかの判定です。
以下のように実装しました。

(**
    ポリゴンに対して指定の点が内包されるか
*)
let isInternal (points : (float * float) array) (x : float, y : float) = 
    // ポリゴンを構成する点から線分を取得する
    // 最後の線分はインデックス0と接続するため、
    // 剰余算を用いる
    let polyLines =
        [| 0 .. points.Length - 1 |]
        |> Array.map (fun i ->
            let p1 = points.[(i    ) % points.Length]
            let p2 = points.[(i + 1) % points.Length]
            (p1, p2)
        )

    polyLines
    // c点からの(1,0)方向の水平線と
    // ポリゴン各辺の交点を持つ辺のみ抽出する
    |> Array.filter (fun ((x1, y1), (x2, y2)) ->
        let (vx, vy) = (x2 - x1, y2 - y1) 

        // 水平線と交点を持つ=y座標が一致するということなので
        // そのときのベクトル倍率を計算する
        let t  = (y - y1) / vy
        let cx = x1 + t * vx
        let cy = y1 + t * vy

        let minY = [ y1; y2 ] |> List.min
        let maxY = [ y1; y2 ] |> List.max

        x <= cx &&              // 左側で交差する場合は除く
        minY <= y && y <= maxY  // 線分外で交差する場合は除く
    )
    // 交点が奇数なら内側
    |> Array.length
    |> fun count -> count % 2 = 1

点と直線の距離の計算は以下のように実装しました。


(**
    直線ABに対するC点を通る垂線の交点座標
*)
let verticalCrossPoint (pointA : (float * float), pointB : (float * float))
                       (pointC : float * float)  = 
    let (ax, ay) = pointA 
    let (bx, by) = pointB
    let (cx, cy) = pointC
    let (abx, aby) = ((bx - ax), (by - ay))
    let (acx, acy) = ((cx - ax), (cy - ay))
    let t = (abx * acx + aby * acy) / (abx ** 2.0 + aby ** 2.0)
    (ax + t * abx, ay + t * aby)

(**
    直線ABと点Cの距離
*)
let distance (pointA : (float * float), pointB : (float * float))
             (pointC : float * float) = 
    let (cx, cy) = pointC
    let (px, py) = verticalCrossPoint (pointA, pointB) pointC
    ((px - cx) ** 2.0 + (py - cy) ** 2.0) ** 0.5

あとはこれらを組み合わせ、内包されるメッシュに対して最も近い辺と同じ色をつけたのが最初に貼り付けた図になります。

図はSVGで以下を使って作成しました。

今度は正攻法に幾何計算を用いて領域分割にトライしたいと思います。

処理手順

解説

Discussion