💔

【平面分割】GISを作るデータ構造とアルゴリズム【計算幾何学】

に公開

概要

前回の記事では, GISで一般的に実装されている図形の結合、切り抜き、交差、、、といった「地図の重ね合わせ」の実装に入った. 特に, 「地図の重ね合わせ」処理のなかで行われている「線分交差」(図形同士の重ね合わせの交点)の処理を, 平面走査法というアルゴリズムで実装した.

今回からは平面走査法を駆使し,いよいよ「地図の重ね合わせ」処理本体を実装する.

図形の再定義

「地図の重ね合わせ」は, 図形の結合(和集合)、切り抜き(差集合)、交差(積集合)とあるように, 図形のブール演算とみなせる. 一般的な記号(True or False)のブール演算と違い, 図形のブール演算はどう行なっているのだろうか? こんなことをまずは想像してみてほしい.

つぎに,前回の記事にも載せたが, それぞれの「地図の重ね合わせ」処理のイメージをみてほしい.
**「重ねあわせ」の処理の入出力は次のようにすべて点集合で表現できる.**と述べた.

  • 入力: 図形P=[p_1,p_2, ..., p_5], 図形 Q=[q_1,q_2,...,q_6]
  • 出力:
    • 融合: 図形[p_1,p_2,A, q_2, ..., B, p_4,p_5]
    • 交差: 図形[q_1,A, p_3, B]
    • 差分: 図形[A, q_2, ..., q_6, B, p_3]

どの「地図の重ね合わせ」処理であっても, 入力時にはなかった交点A,Bが必要になる可能性があったことから,前回の記事では平面走査法による線分交差をまず紹介した.

今度の疑問はこうだ.
『「地図の重ね合わせ」処理の種類(融合、交差、差分)に合わせて, 出力を変える方法は?」「

人間であれば, 絵をみえるだけで, どこが融合、交差、差分の領域かわかるが, 計算機には自明ではない. 具体的には, 融合処理の結果として, 「なぜ頂点p_2を通って交点Aに差し掛かった段階で, 頂点p_3,q_1ではなく,頂点q_2へと向かうのか...」と考えないといけない.

ここで述べたいのは, 新たに「面」の定義が必要だということだ. また,ブール演算を実施するためには, 明示的に「面」の意味(ラベル)を識別できる必要がある.

これらを次に述べる 「二重辺リスト」と「境界サイクル」 を用いて改善する.

二重辺リスト

まずは, 図形を形成する辺の定義をしていこう.

図形は複数の辺がつながって構成されているため, 線分(辺)の集合とみなせることは想像に難くないだろう. ただここでは, 無向辺ではなく, 2本の有向辺で元の図形の一辺を表現する.
そしてそれらの2本の有向辺は互いに向きが反対となっている.

こうした2本の有向辺の複数集合として図形を表現したデータ構造を二重辺リストとよぶ.

上の図に二重辺リスト化した図形P,Qを示していく. いままでの無向の辺が, 向きが相互に反対の有向辺2つに置き換わっているのでわかりやすいだろう.

なぜ二重辺なのか,,,それは,図形の内側と外側を表現できるようにするためだ.
ここで図形の内側と外側を有向辺の向きに合わせて定義する.

  1. 図形の外側:図形を構成する有向辺を反時計回りに移動する
  2. 図形の内側:図形を構成する有向辺を時計回りに移動する

上図の通りに,図形の内側と外側を有向辺の向きに合わせて定義しておくとのちのち便利だ.

また, 交点周りなど, 複数の移動先がある場合の曲がり方のルールを決めておこう.
上図の右のように, 有向辺の終点を今度は始点とする有向辺が複数ある場合, もっとも時計回り側, つまり, 自身に対してもっとも右側に向かって移動する有向辺へと次の経路を進める[1]

以後,この曲がり方ルールを, 時計回りルールと呼ぶ.

二重辺リストと時計回りルールを利用すると, 図形の分割範囲を明確に定義できる!! 👍

方法は簡単で,

  1. どれでも好きな有向辺()を1つ選び, 時計回りルールで次の有向辺へ移動していく.元の有向辺()に戻ってきたら, 走査済みとして通った辺をラベルづけする
  2. 次に, 走査済みではない有向辺(*)をまたランダムに選び,1.と同じ時計回りルールで移動する.
  3. 1,2を繰り返し, すべての有向辺が走査済みとなったら終了

この方法を図形P,Qの重ね合わせ後の図に適用すると,計4つの巡回ルートが見つかる.
うち3つは分割後の個々の図形A,B,Cで,のこりは, 全図形を合わせたものになる.

なんと, 二重辺リストをうまくつかえば,図形という面をうまく区別できそうだ.

境界サイクル

めでたしめでたし...とはならず, まだ区別されたそれぞれの図形が何者かがわからない.つまり, 図形Aは, 元の図形P,Qのどのような重ね合わせ処理の結果で生まれたのか識別できていない.

この分割後の図形を識別する情報を,面分割後の二重辺を使って定義していこう.

上の図を例に考える. まず, 境界サイクルをすべてみつけよう. 境界サイクルは上述した時計回りルールで二重辺をたどっていった場合にできる一巡の辺集合である.

上の図だと, 境界サイクルは① ~ ⑥までできる. あと, それぞれの境界サイクルが図形の内側と外側のどちらを構成するかを分けて整理する. (辺の回転方向が反時計や時計回りかで判別)

  1. なにかの図形の外側: ②、③、⑥
  2. なにかの図形の内側: ①、④、⑤

次に上の図の青色部分の面をどのように定義するか考える. 見る感じ, 境界サイクル②、③、⑤を使えば、この面の境界を定義できそうだ、と気づくと思う.

逆にいえば, 境界サイクル②、③、⑤は, ある面に関連するものであると, 関係性を紐づける必要がある.

また「どうすれば💦」となる読者もいると思うが、ここでテクニック。外側境界サイクルに着目する.具体的には次のように各境界サイクルの関係性を分析する.

  1. すべての境界サイクルをノードとして登録(この時点ではまだエッジはない)
  2. まだ走査していない外側境界サイクルをひとつ選ぶ
  3. 選んだ外側境界サイクルを構成する頂点のうち, 最もx座標が小さいものを探す.
  4. 3で選んだ頂点からx座標の小さい方向へ平行線を伸ばす.
  5. 4の線が内側境界サイクルに当たった場合、この内側境界サイクルと,2の外側境界サイクルを示すノード間をエッジで結ぶ
  6. 4の線が外側境界サイクルに当たった場合、この外側境界サイクルと,2の外側境界サイクルを示すノード間をエッジで結ぶ
  7. すべての外側境界サイクルが走査済みになるまで2-6をくり返す

この手続きを図の図形に適用すると、4つの部分木(グループ)に分かれる. 一番外側の⑥を除くと, 残りのグループそれぞれを構成する境界サイクルで囲まれた面は、それぞれ重なりのない関係になっていることがわかる. これらの面はすべて,面分割後の独立した面(重なりのない面)であり, それらすべて1つまたは複数の境界サイクルで定義できる.

ブール演算

面分割後の独立した面(重なりのない面)であり, それらすべて1つまたは複数の境界サイクル, つまりは分割後の二重辺リストで定義できることがわかった.あともう少し.

最後に, この分割後の二重辺リストがもともとの図形のうちどこに所属しているものか、という情報を使おう.

わかりやすく冒頭の二つの図形(PとQ)の重ねあわせの例に戻る. いま, "A"の領域を構成する境界サイクルは, オリジナルの図形P,Qのそれぞれ内側と外側の辺で構成されている. 同様に, "B","C","D"もオリジナルの図形P,Qの内側か外側のセットの辺で構成されていることがわかる.

ここで, "内"を0,"外"を1として, 分割後の面を整理する.すると,図のような表ができあがる.[2]

結論をいえば, この表を使えば面分割のブール演算が導けてしまう.
ひとつ例をとると, P \cup Qは, (P=1) OR (Q=1)の関係なので, 表でP=1Q=1となっている面を集めると,P \cup Q = A + B + Cで表現できることができる.🥳

無事, 図形のブール演算に到達できた.

つまり...

ここまで長々と話したが, ポイントを押さえておこう.

  1. 二重辺リストを使って、面の外側・内側を定義.このとき, 各辺がオリジナルのどの図形を構成するものかを記録しておく.
  2. 平面走査法で交点を探しつつ, 1の二重辺リストの関係性や数を更新.
  3. 分割後の独立の面を, 2の分割処理後の二重辺リストを使って定義(境界サイクルとして)
  4. 3を通して、分割後の独立の面がオリジナルの図形との関係性を整理する
  5. 4の表を使ってブール演算を実施する.

データ構造とアルゴリズム

前節では, 例をつかいながら「どのようにしたら図形のブール演算が定義できそうか」を見てきた.
何よりも図形を二重辺リストで定義することがキーになることがわかっただろう.

二重辺リストには,

  1. 元の図形への所属関係も保存しなければいけない.
  2. 境界サイクルを定義できるよう各二重辺には前後関係が決められている.
  3. 最後に, 二重辺リストは面分割後で更新されている

という特徴をおさえる必要がある. これを, 面レコード、辺レコード(,頂点レコード) というデータ構造で機能をつくる.

イメージ

以下の表のような配列やリストをつくる.

講義資料より引用

Python実装例
二重辺リスト(面レコード、辺レコードなど)
def create_doubly_connected_list(polygons, polygon_label, line_label):
    '''
    与えられたポリゴンから二重辺リスト
    Input:
        polygons: 計算領域内の複数のポリゴン情報 [[[x,y],[x,y],...], [[x,y],[x,y],...]]
        polygon_label: 各ポリゴンのラベル 
    Output:
        vertex_xy : 全ポリゴンの座標情報 {"vid":[x,y],...}
        link_list : 片辺の頂点の接続関係 [[p,q],[q,p],...]
        v_adj_link: 頂点ごとの接続辺ID
        polygon_ps: ポリゴンごとに所属している点ID
        twin      : ある片辺の双子辺リスト
        lprev     : 接続-前リンク
        lnext     : 接続-後リンク
        plane_record: 面レコード
    '''
    vertex_num = 0               # 頂点数
    link_num   = 0               # リンクの数(両辺)
    
    vertex_xy  = {}              # 先頭はダミー
    link_list  = [[-1,-1]]       # 先頭はダミー
    v_adj_link = [-1]            # 頂点ごとの接続辺ID

    polygon_ps = [[-1]]          # ポリゴンごとの 所属している点ID
    twin       = [-1]            # 双子辺
    lprev      = [-1]            # 接続前リンク
    lnext      = [-1]            # 接続後リンク
    
    plane_label = []

    """
    頂点レコードの作成
    """
    for i, polygon_i in enumerate(polygons):
        size      = len(polygon_i) 
        _within_p = []
        for i in range(size):
            x,y   = polygon_i[i]
            index = in_vertex(x,y, vertex_xy)
            if index < 0:                                 # 未登録ならば
                v_adj_link.append(-1)                     # 頂点数に合わせてサイズを大きくしていく
                vertex_num += 1
                vertex_xy[vertex_num] = [x,y]
                index = vertex_num
            if i > 0:
                px,py  = polygon_i[i-1]
                _index = in_vertex(px,py, vertex_xy)     # 一つまえの点
                if _index < 0:
                    raise ValueError
                link_list.append([_index, index])         # リンクの追加
                link_list.append([index, _index])         # 反対側も
                link_num += 2
                twin.extend([link_num, link_num - 1])      # リンクを見つけるごとに双子辺を追加していく
                if v_adj_link[_index] == -1:              # とりあえず最初に見つかったリンクをつなげておく
                    v_adj_link[_index] = link_num - 1     # 片辺リスト
            _within_p.append(index)
        polygon_ps.append(_within_p)

    for i, value in enumerate(v_adj_link):
        if i > 0:
            if value == -1: # まだ未決定なものがあれば
                _index = in_adj_list(i, link_list)
                v_adj_link[i] = _index # 片辺リスト
                
    """
    片辺レコードの作成
    """
    lprev = [-1]*(len(twin))
    lnext = [-1]*(len(twin))
    for polygon_p in polygon_ps:  # ポリゴンごとにリストの接続関係を登録していく
        psize = len(polygon_p)
        for pi in polygon_p:
            if pi > 0:
                link_from_p = search_for_next_vertex_from_p(pi, link_list) # piを始点とするリンクを抽出
                for link_id, start_p in link_from_p.items(): # 時計回りと反時計回り
                    start_link_id   = link_id
                    current_link_id = link_id
                    next_p          = start_p
                    while True:
                        link_from_p2 = search_for_next_vertex_from_p(next_p, link_list) # next_pを始点とするリンクを抽出
                        next_link_id = -1
                        for link_id2, _p in link_from_p2.items():
                            if link_id2 == twin[current_link_id]:
                                continue
                            else:
                                next_link_id = link_id2
                                next_p = _p
                        lprev[next_link_id]          = current_link_id
                        lnext[current_link_id]       = next_link_id
                        lnext[twin[next_link_id]]    = twin[current_link_id]
                        lprev[twin[current_link_id]] = twin[next_link_id]

                        current_link_id = next_link_id
                        if current_link_id == start_link_id:
                            print("2重辺リストの作成完了")
                            break
                    break
                break  # 一度回せばよい
    print(link_num)
    """
    面レコードの作成
    """
    plane_record = [-1]*(link_num+1)
    polygon_record = [-1]*(link_num+1)

    plane_label = [[value + "外", value] for value in polygon_label] # ["外側成分", "内側成分"]

    for pi, polygon_i in enumerate(polygons):
        mp  = search_for_min_x(polygon_i) # ポリゴン内で最小の点(x座標最小, 等しい場合はy座標最小)
        mp  = polygon_ps[pi+1][mp]
        mp_start = [] 
        # mpを始点とする線分を二つ取り出す
        for i, value in enumerate(link_list):
            if i > 0:
                if mp == value[0]:
                    mp_start.append(value[1])
        px, py = vertex_xy[mp]
        qx, qy = vertex_xy[mp_start[0]]
        rx, ry = vertex_xy[mp_start[1]]
        theta1 = np.degrees(np.arctan2(qy-py,qx-px))
        theta2 = np.degrees(np.arctan2(ry-py,rx-px))

        if theta1 > theta2: # mp_start[0]へのリンクが外側境界
            link_id = search_for_link_id(mp, mp_start[0], link_list)
            while plane_record[link_id] == -1:
                plane_record[link_id] = plane_label[pi][0]
                l                     = (link_id-1) // 2
                polygon_record[link_id] = line_label[l]
                link_id               = lnext[link_id]
            link_id = search_for_link_id(mp, mp_start[1], link_list)
            while plane_record[link_id] == -1:
                plane_record[link_id] = plane_label[pi][1]
                l                     = (link_id-1) // 2
                polygon_record[link_id] = line_label[l]
                link_id               = lnext[link_id]
        else:
            link_id = search_for_link_id(mp, mp_start[1], link_list)
            while plane_record[link_id] == -1:
                plane_record[link_id] = plane_label[pi][0]
                print(link_id)
                l                     = (link_id-1) // 2
                polygon_record[link_id] = line_label[l]
                link_id               = lnext[link_id]
            link_id = search_for_link_id(mp, mp_start[0], link_list)
            while plane_record[link_id] == -1:
                plane_record[link_id] = plane_label[pi][1]
                l                     = (link_id-1) // 2
                polygon_record[link_id] = line_label[l]
                link_id               = lnext[link_id]
    print("-----------------------------")
    print("          頂点レコード        ")
    print("-----------------------------")
    for i in range(vertex_num):
        if i > 0:
            print(str(i) + " "*5 + str(vertex_xy[i][0]) + " "*3 + str(vertex_xy[i][1]) + " "*5 + str(v_adj_link[i]))
        else:
            print("PID" + " "*4 + " x " + " "* 3 + "y" + " "*5 + "Cij")
     
    print("-----------------------------")
    print("          片辺レコード        ")
    print("-----------------------------")
    for i, value in enumerate(link_list):
        if i > 0:
            print(str(i) + " "*7 + str(value[0]) + " "*5 + str(twin[i]) + " "*5 + str(plane_record[i]) +" "*5 + str(lnext[i]) + " "*5 + str(lprev[i]))
        else:
            print("LinkID" + " "*2 + "始点" +" "*3 +"双子辺" + " "*3 + "面" +" "*3 + "次" + " "*3 + "前")
    
    return vertex_xy, link_list, v_adj_link, polygon_ps, twin, lprev, lnext, plane_record, plane_label, polygon_record

片辺レコードについては, 線分の分割, つまり, 平面走査法により交点を見つけるたびに更新される.

講義資料より引用

ここでの処理の重要なポイントは,

  1. 交点を起点とした新しい片辺を追加すること. また, 既存の片辺の頂点ペアを更新すること.
  2. 片辺と双子辺の関係性を追加・更新すること.
  3. 片辺同士の前後関係を、時計回りルールで更新すること.

である.

最後に, 面分割後の片辺レコードを使って, 独立した面の特定・オリジナルの図形との関係性を整理する.

分割面のブール演算
def create_planar_label(vertex_xy,link_list, twin, lnext, lprev, left_link_event, plane_record):
    '''
    分割後の面に識別用のラベルを張り付ける処理
    Input:
        vertex_xy : 対象となる計算領域内のポリゴンの頂点座標 {"頂点ID": [x座標, y座標],...}
        link_list : 二重辺リスト [[始点p,終点q], [始点q, 終点q],....]
        twin      : 双子辺リスト
        lnext     : ある辺の巡回ルートのうち次に隣接する辺
        lprev     : ある辺の巡回ルートのうち前に隣接する辺
    Output;
        link_group_label: 各辺が所属する面のラベルリスト
        group_connect   : 面ラベルの接続関係
    '''
    vertex_visited_list = [] # 頂点ごとの走査完了リスト
    vertex_connect_list = {} # 頂点ごとの残りの接続辺のリスト
    link_visited_list   = [] # 片辺ごとの走査完了リスト
    link_group_label    = [] # リンクのグループ割り

    vertex_visited_list = [False for _ in range(len(vertex_xy)+1)] # 初期化(False)
    link_visited_list   = [False for _ in range(len(link_list))]   # 初期化(False)
    link_group_label    = [-1 for _ in range(len(link_list))]      # 初期化

    for pi in vertex_xy.keys():
        from_ps    = search_for_next_vertex_from_p(pi, link_list)
        from_links = list(from_ps)
        to_ps      = search_for_prev_vertex_to_p(pi, link_list)
        to_links   = list(to_ps)
        connect_link = from_links + to_links
        vertex_connect_list[pi] = connect_link
    
    '''
    ここから各閉じた面に個別のラベルを当てていく
    '''
    group_ID = 0
    group_connect = ["F0"] # 無限境界のラベルF0
    while True:
        sp = return_min_vertext(vertex_xy, vertex_visited_list) # 走査未完了の最小位置頂点を探す
        if sp > 0:
            # 外回りと内回りのリンクの向きを確定させる
            while sp in vertex_connect_list:                  # spにまだ接続リンクがあったら
                group_ID += 1
                rotation_links      = vertex_connect_list[sp] # spを接続とする片辺リスト
                link1               = rotation_links[0]       # とりあえず一つだしてみる(回してみる)
                loop_list           = [link1]
                link_visited_list[link1] = True               # 走査済み
                clink               = lnext[link1]            # 次のリンクに向かって回転開始
                while clink != link1:                         # 同じ片辺にかえって来ない限り
                    loop_list.append(clink)
                    link_visited_list[clink] = True
                    clink                    = lnext[clink]   # 次へ
                
                # 周回終わりには nlinkとplinkの角度を比較する
                # 点spから出ていくベクトルを軸として角度を計算
                npi,nqi             = link_list[clink]
                if npi == sp:                                 # このベクトルは頂点spから出て行くベクトル
                    plink           = lprev[clink]            # 前のベクトル
                    ppi, pqi        = link_list[plink]        #  p -> q
                    tvector         = np.array([vertex_xy[ppi][0] - vertex_xy[pqi][0], vertex_xy[ppi][1]-vertex_xy[pqi][1]], dtype = np.float64)
                    bvector         = np.array([vertex_xy[nqi][0] - vertex_xy[npi][0], vertex_xy[nqi][1]-vertex_xy[npi][1]], dtype = np.float64)
                elif nqi == sp:                               # このベクトルは頂点spに入っていくベクトル
                    nlink           = lnext[clink]            # 前のベクトル
                    ppi, pqi        = link_list[nlink]        #  p -> q
                    bvector         = np.array([vertex_xy[pqi][0] - vertex_xy[ppi][0], vertex_xy[pqi][1]-vertex_xy[ppi][1]], dtype = np.float64)
                    tvector         = np.array([vertex_xy[npi][0] - vertex_xy[nqi][0], vertex_xy[npi][1]-vertex_xy[nqi][1]], dtype = np.float64)
                
                # 走査したリストを排除しながらラベル付け        
                for link_id in loop_list:
                    link_group_label[link_id] = group_ID
                    for p in link_list[link_id]:
                        if link_id in vertex_connect_list[p]:
                            llist  = vertex_connect_list[p]
                            lindex = llist.index(link_id)
                            llist.pop(lindex)
                            vertex_connect_list[p] = llist 
                        if len(vertex_connect_list[p]) == 0:   # 空っぽになった
                            del vertex_connect_list[p]
                            vertex_visited_list[p] = True      
                 
                angle = vector_to_vector_rotation_angle(bvector, tvector) # bvector -> tvectorへの回転角を求める
                if angle < 0:
                    angle += 360
                if angle < 180:                               # この一連のサイクルは外側境界のサイクル
                    group_connect.append("F" + str(group_ID))
                else:                                         # この一連のサイクルは穴の境界 -> 外側境界との連結へ
                    if left_link_event[sp] == -1:             # この時は連結先は無限境界
                        group_connect.append(0)               # Index 0 つまり無限境界"F0"を示すように
                    else:
                        left_neighbour_link = left_link_event[sp] # 隣のリンクがある
                        group_label         = link_group_label[left_neighbour_link]
                        group_connect.append(group_label)     # 数字で接続
        else:
            break                                             # 走査完了
    print(link_group_label)
    print(group_connect)
    print("-----------------------------")
    print("          片辺レコード        ")
    print("-----------------------------")
    for i, value in enumerate(link_list):
        if i > 0:
            if i >= len(plane_record):
                print(str(i) + " "*7 + str(value[0]) + " "*5 + str(twin[i]) + " "*5 + "新規" + " "*5 + "F" + str(link_group_label[i]) + " "*5 + str(lnext[i]) + " "*5 + str(lprev[i]))
            else:
                print(str(i) + " "*7 + str(value[0]) + " "*5 + str(twin[i]) + " "*5 + str(plane_record[i]) + " "*5 + "F" + str(link_group_label[i]) + " "*5 + str(lnext[i]) + " "*5 + str(lprev[i]))
        else:
            print("LinkID" + " "*2 + "始点" +" "*3 +"双子辺" + " "*3 + "面" +" "*3 + "新ラベル"+ " "*3 + "次" + " "*3 + "前")
    return link_group_label, group_connect

def create_polygon_from_label(link_group_label, group_connect, filter_func, plane_label, plane_record, polygon_record, vertex_xy, link_list, lnext):
    '''
    新しい面ラベルと更新前のラベルの関係から 必要となるラベルをもつ面、および構成する頂点のリストを返す
    Input:
        link_group_label: 各辺が所属する面のラベルリスト
        group_connect   : 構成された新ラベルのリスト
        filter_func     : 抽出したいラベル(1:対象ポリゴン,0:対象外ポリゴン)
        plane_label     : 各ポリゴンの [外側成分, 内側成分]の多重リスト
        plane_record    : 各片辺の内側ポリゴンのラベルリスト
        polygon_record  : 各片辺の元の所属ポリゴンのリスト(内外の区別なし)
    Output:
        target_polygon  : 選択したラベルに一致した分割後のポリゴン [ [[x,y],...], [[x,y],...]
        lable_list      : target_polygonの並びに連動した(旧)ラベル [[label1,label2,...], [label1,label2,...]]     
    '''
    '''
    無限境界への接続関係を修正
    '''
    modify_label = []
    for plane_id in group_connect:
        if plane_id == "F0": # 無限境界側
            l = group_connect.index("F0") # "F0"の要素番号
            for ind, f in enumerate(group_connect):
                if ind != l:
                    if f == l:
                        modify_label.append(ind)
    """
    ここには 自動的に元の情報と新しい情報のラベルの対応付けの処理が必要となる
    """
    label_dict = {}
    for i, label in enumerate(link_group_label):
        if label in modify_label:
            label = 0 # 無限境界へとラベルを修正
        if i > 0:
            if i < len(plane_record):
                prev_label = plane_record[i]
                if label in label_dict:
                    a = label_dict[label]
                    if prev_label in a:
                        continue
                    else:
                        a.append(prev_label)
                        label_dict[label] = a
                else:
                    label_dict[label] = [prev_label]

    print(label_dict)
    """
    ラベルごとのポリゴンを抽出
    """
    # 検索用のビットとする(1はポリゴン内,0はポリゴン外, それぞれ該当するポリゴンとする)   
    filter_label = set([plane_label[i][choice] for i, choice in enumerate(filter_func)])
    print(plane_label)
    print("filter_label:",filter_label)
    target_label = []
    for label, value in label_dict.items():
        if set(value) == filter_label:
            print(label)
            target_label.append(label) 
    if len(target_label) == 0: # 見つからなかったら
        for i, fc in enumerate(filter_func): 
            if fc == 1:
                filter_label = set([plane_label[i][1]])
                print(filter_label)
        for label, value in label_dict.items():
            if set(value) == filter_label:
                target_label.append(label)

    target_polygon = []
    label_list     = []
    print("target_label:",target_label)
    if target_label == [0]:
        target_label = [1]
    print(link_group_label)
    for tlabel in target_label:
        a              = [] 
        l              = [] 
        st_link        = link_group_label.index(tlabel)
        print(tlabel, st_link)
        p,q            = link_list[st_link]
        a.append(vertex_xy[p])
        l.append(polygon_record[st_link])
        next_link = lnext[st_link]
        while p != q:
            a.append(vertex_xy[q])
            l.append(polygon_record[next_link])
            _, q      = link_list[next_link]
            next_link = lnext[next_link]
        a.append(vertex_xy[p])

        target_polygon.append(a)
        label_list.append(l)

    return target_polygon, label_list

実施例

最後に, 構築したプログラムを通して, 面分割のブール演算を試した結果を載せる.

演算の内容に合わせて処理後の図形も変化していることがわかる.検証図形は冒頭の図PとQをイメージしたものである.

まとめ

今回は, GISで最も多用されるであろう「地図の重ね合わせ」「面分割のブール演算」を説明した.
平面走査法+二重辺リストを駆使して、分割前・後の各面を定義することが重要であるとみてきた.

次回は, ボロノイ分割か三角形分割か、などを検討しているのでお楽しみに.

雑談

あまりにもタフすぎる🥹.
「地図の重ね合わせ」「面分割のブール演算」を説明する記事は「ない」と思っているので, すごく貴重な記事を残せたと自負している. が, 本プログラムはあくまでも個人例なので, もっと上手な実装があればそちらを参考にしてほしい.

前回同様, 要望があれば全体コードをGithubで公開予定である.

脚注
  1. 図では,A,B,Cへの選択肢のうち,Aを選ぶ ↩︎

  2. なんか見たことがある... ↩︎

Discussion