複数の閉領域をマージして 結合後の領域の座標点を取得したい

Question

実現したいこと

(X, Y)の座標点を複数含んだ閉領域が存在している状況で、閉領域同士が接している際に、双方の領域をマージしたいです。
また、マージ後の領域について、外形を示す形状座標点を取得できるようにしたいです。

前提

Pythonで解析用の領域を管理させるシステムを開発中です。
ユーザーの都合上、同じ条件の領域を複数構築して、それを結合させることで大きな領域を定義させたいと考えています。(詳細は伏せさせてください)
当初、単純に線分の一致、ベクトル等を考慮して判断させようかと考えていたのですが、線分が接するけど両者の長さが異なる例などが多々あり、単純に処理ができない状況です。
プログラミング初学者につき、検討がつかないので質問させていただいた次第です。エラーではなく、計算のアルゴリズム等についてなのですが、アドバイスいただけますと幸いです。

例:対象としている閉領域の座標点(X, Y)

領域1: (10,0)(10,10)(20,10)(20,15)(25,15)(25,5)(30,5)(30,0)
領域2: (5,25)(25,25)(25,15)(20,15)(20,10)(15,10)(15,20)(5,20)

例:出力したいマージ後の形状座標点

(5,25)(25,25)(25,5)(30,5)(30,0)(10,0)(10,10)(15,10)(15,20)(5,20)

処理概要と出力イメージ

Comments

[winterboum]

ここでの質問は
(25,15)(25,5) と (25,25)(25,15) が接している部分だ、というのを「発見する」部分では不要で
発見後の処理のみですか？

[退会済みユーザー]

> 閉領域同士が接している
斜め同士はありなのか隣同士だけなのか等、図示イメージがあった方が回答が得やすいです。
後は画像処理系をPythonでやろうとするとOpenCVやScipyあたりが相性が良いと思います。イメージ的にはちょっと変な形（縦横だけなのかななめもありなのかで変わるはずです）のカーネルでモルフォロジー変換をすればいけそうな気がします。

[YuYH]

ご回答ありがとうございます。
処理内容としては「発見する」部分も悩みどころでして、長さは違うけど一直線上に存在する線分同士をどう判定させるのか等について悩んでいます。
やはりここは座標点の情報から計算させるのではなくて画像処理を活用するのが得策なのでしょうか

[TakaiY]

領域の線分は水平垂直のみですか？
例では領域は接しているもののみですが、領域の重なりを考慮する必要がありますか？

[winterboum]

TakaiYさんのコメントは「発見する」部分に関してですね？
ですと YuYHさんも一つ、
同じ座標かどうか の判断では 実数の微妙な差は気にせず 数値１==数値２ で比べて良いのか、「ある誤差範囲なら同じ」と見るべきなのか。微妙な差で 重なったりすき間があいたり する可能性が在るので。

[fana]

例を見た感じだと「接している（重なってはいない）」ならば，「一方の図形の頂点が他方の図形の辺上にある」っていう形になりますよね．
で，図形の定義をグラフだと見たとき，そういう頂点を全組み合わせチェックででもして探せば，グラフを１つにすることはできますよね．
で，「外周を辿る」とき，分岐点（「そういう頂点」のところ）ではどっちに行けばいいのか？っていうと，「今来た方向から見てどっち側の方向に辿れば良いのか」っていう話がありますよね．（右回りとか左回りとか）

…っていう感じでやれませんかね．結果として辿らなかった場所が「接した部分」になります．
（あー，でもこんなのだとマージ結果として穴ができる下側の例を扱えないか．辿らなかった部分グラフがループになってるか？とか，あるいはその中の頂点に「そういう」じゃないやつがいないか？　とかそういうチェックをして穴に対処する必要がありそう．）

[winterboum]

たどったところを消して新しい 道筋に登録していく、結果 消し残りが 接しているか穴。ということで後処理で出来ないかな

Answer 1

Shapelyのunary_unionが使えそうです。
結合後、直線区間中に連続する点は残りますが、それを除く方法はおそらく難しくなく、Shapelyの関数としても用意されているかもしれません。

追記

不要な点の削除はnormalizeしてsimplifyすればよい感じです。
参考：Simplify method can not simplify duplicated vertices of a ring #1638

import matplotlib.pyplot as plt
import shapely
from shapely.geometry import Polygon

# ポリゴンは反時計回りで定義
p1 = Polygon([(10,0),(10,10),(20,10),(20,15),(25,15),(25,5),(30,5),(30,0)][::-1])
p2 = Polygon([(5,25),(25,25),(25,15),(20,15),(20,10),(15,10),(15,20),(5,20)][::-1])
p3 = shapely.unary_union([p1,p2])

plt.plot(*p3.exterior.xy, marker="o")
plt.show()
print(p3)
# POLYGON ((25 15, 25 5, 30 5, 30 0, 10 0, 10 10, 15 10, 15 20, 5 20, 5 25, 25 25, 25 15))

# 不要な点を取り除く
p3 = p3.normalize().simplify(0.1)
plt.plot(*p3.exterior.xy, marker="o")
plt.show()
print(p3)
# POLYGON ((5 20, 5 25, 25 25, 25 5, 30 5, 30 0, 10 0, 10 10, 15 10, 15 20, 5 20))

# 重ならない
p1 = Polygon([(0,0),(1,0),(1,1),(0,1),(0,0)])
p2 = Polygon([(2,0),(3,0),(3,1),(2,1),(2,0)])
print(shapely.unary_union([p1,p2]))
# MULTIPOLYGON (((0 0, 1 0, 1 1, 0 1, 0 0)), ((2 0, 3 0, 3 1, 2 1, 2 0)))

# 内包
p1 = Polygon([(0,0),(3,0),(3,3),(0,3),(0,0)])
p2 = Polygon([(1,1),(2,1),(2,2),(1,2),(1,1)])
print(shapely.unary_union([p1,p2]))
# POLYGON ((0 0, 0 3, 3 3, 3 0, 0 0))

# 辺で接する
p1 = Polygon([(0,0),(1,0),(1,1),(0,1),(0,0)])
p2 = Polygon([(1,0),(2,0),(2,1),(1,1),(1,0)])
print(shapely.unary_union([p1,p2]))
# POLYGON ((1 0, 0 0, 0 1, 1 1, 2 1, 2 0, 1 0))

# 点で接する
p1 = Polygon([(0,0),(1,0),(1,1),(0,1),(0,0)])
p2 = Polygon([(1,1),(2,1),(2,2),(1,2),(1,1)])
print(shapely.unary_union([p1,p2]))
# MULTIPOLYGON (((1 1, 1 0, 0 0, 0 1, 1 1)), ((1 1, 1 2, 2 2, 2 1, 1 1)))

# 重複
p1 = Polygon([(0,0),(2,0),(2,2),(0,2),(0,0)])
p2 = Polygon([(0,1),(2,1),(2,3),(0,1)])
print(shapely.unary_union([p1,p2]))
# POLYGON ((2 1, 2 0, 0 0, 0 1, 0 2, 1 2, 2 3, 2 2, 2 1))

# 穴ができる
p1 = Polygon([(0,0),(3,0),(3,1),(0,1),(0,0)])
p2 = Polygon([(2,1),(3,1),(3,3),(0,3),(0,1),(1,1),(1,2),(2,2),(2,1)])
print(shapely.unary_union([p1,p2]))
# POLYGON ((3 1, 3 0, 0 0, 0 1, 0 3, 3 3, 3 1), (1 1, 2 1, 2 2, 1 2, 1 1))

Comments

[YuYH]

検討ありがとうございます。
Shapelyのライブラリを活用すれば特に問題なく領域のマージと座標点の取得まで行うことができました。

Answer 2

fourteenlengthさんから別解が入るかも、を期待します。

いまさらですが…。
しかもコメント欄で「モルフォロジー変換云々」とっぽいことを言っていたくせに、いざ実装する段階では普通に塗りつぶしで検知するという暴挙に出ております。

追記：
適当な座標データがないので試していませんが、ヒエラルキーの値を考慮してDrawContoursしてFindContoursをすれば、内側の空白（琵琶湖みたいな構造）も残せるはずです。

---

[10  0]
[10 10]
[14 10]
[15 11]
[15 19]
[14 20]
[ 5 20]
[ 5 25]
[25 25]
[25  6]
[26  5]
[29  5]
[29  0]

import cv2
import numpy as np

# 座標を定義
area_1 = np.asarray([[10, 0],[10,10],[20,10],[20,15],[25,15],[25, 5],[30, 5],[30, 0]])
area_2 = np.asarray([[ 5,25],[25,25],[25,15],[20,15],[20,10],[15,10],[15,20],[ 5,20]])

# ブランクイメージを作る
img = np.zeros((30,30))

# 塗りつぶす
img = cv2.drawContours(img ,[area_1], -1, 255, -1)
img = cv2.drawContours(img ,[area_2], -1, 255, -1)

# 座標を知りたい
# REF:http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_contours/py_contours_begin/py_contours_begin.html#id5
cnts, _ = cv2.findContours(img.astype(np.uint8),cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)

# 様子見
cnts = np.asarray(cnts).squeeze()
for cnt in cnts:
    print(cnt)

Answer 3

では「発見する」のは別の課題にさせていただいて(考えていなかったので）「結合させる」方について。
画像処理には知見がないので、画像処理 でこの問題をどの様に処理できるのか、はわかりません。fourteenlengthさんから別解が入るかも、を期待します。素人の脊椎反射では 「同じ色の領域を拾い出す」というような処理が使えるのかな？ と
で、
座標点の情報から計算 する方法です。
図をみると「接しているのは一つの線分だけ」ではないのですね。これらを一度に結合させるのではなく、一つずつ結合してください。

手順を５つに分けます
[1] 共通線分を抜き出す：今回はこれには触れない。連続している部分はひと塊に抜き出しておく。
[2]前処理。連続している部分を縮退させる
[3]１つ目の結合。これにより一つの領域になる。ただ切れ目が入っている
[4]一つ目の結合の切れ目を塗り消す
[5]残りの切れ目を塗りつぶす

番号付ですが、時計回りにつける と標準化してあるとします。
[1] は出来ているものとしますが、上の図の様に連続して接している場合は A1-A2 が B4-B3と A2-A3 が B3-B2と、、、とバラバラではなく A1-A2-A3-A4 が B4-B3-B2-B1 と接している と言うようにまとまってる様にしておく。
[2] 上の図だと A1-A2-A3-A4 が B4-B3-B2-B1 で接しています。
ここでは階段１stepですが、10段あって、凸凹にも斜面にもなってる、ことでしょう。
すると例えば A1-A2-A3-A4-A5-A6-A7-A8 が B8-B7-B6-B5-B4-B3-B2-B1 で接している。
このようなときに両端の共通だけ残して A1-A2-A7-A8 が B8-B7-B2-B1 で接している、と直します。
領域情報からも共通線分情報からも A3-A6、B3-B6情報は削除。
A2からA7(すなわち B7からB2)への斜面にしてしまう感じ。

[3]１つ目の共通部分の結合
領域AのA1-A2-A7-A8 が 領域BとB8-B7-B2-B1 で接している。このとき A1-A2とB8-B7で結合させます。
AからBへ移るのは
A1,A2 と B8の関係は(1) A1の外に在る。(2) A1-A2間にある。
(1)の場合 A1の次にB8へバック(A2と逆方向）しB7へ進む
(2)の場合 A1の次にB8へ進み(A2方向）B7へ進む
あら、同じですね。ので A0→A1→B8→B9 となります。

BからAに移るのは
(3) B7がA1,A2 間にある
(4) B7がA1,A2の外に在る
(3)の場合、B7からA1,A2方向へ進み、A2まで行き、A3へ向かう
(4)の場合、B7からA2に進み、A2まで行き、A3へ向かう
これも同じなので B2→B7→A2→A3 となります。

[4]一つ目の結合の切れ目を塗り消す
ここまでくると領域は一つですが、A8-A7-A2-B7-B2-B1 という切れ目が入っています。
このうち A7-A2、B7-B2は完全に重なってますから、領域情報からも共通線分情報からも削除します。切れ目は A8-A7-B2-B1 という折返しです（A7とB2は同じ場所）
(1)B1がA8-A7の外に在る。 B1→ B2(A7) へ行って → A8へ戻る。のですから B1→A8
(2)B1はA8-A7の内にある。 B1→ B2(A7) へ行って → A8へ戻る。のですから B1→A8
どちらも同じですから B0→B1→A8→A9

[5]残りの切れ目を塗りつぶす
下の図です。上をA領域、下をB領域とすると、Aから時計回りで最初の共通を処理したから、右が完結で左が残ってます。
この図では直線ですが、実際には段が在るでしょう。その段も前処理で一段になってます。
すると A1-A2-A7-A8 と B8-B7-B2-B1 の切れ目。この内 A2-A7 と B7-B2は共通ですから消す。すると A1-A2-B7-B8、A8-A7-B2-B1 という2つの折返しの切れ目となります。
この部分は[4]でわかったように A1-B8、A8-B1 になります。

頭の中で図を書いて番号振ってるので、typo在るでしょう、です。
ので、
質問するときは上の回答に従った 図を書いてそれ載せていただけたら。

追記
切れ目を消すという作業は 「隣り合った２つの座標が同じ位置だったら、ふたつまとめて削除」というやり方で、[4][5]行けそうかも。考えてください。