pythonを使って線の数を数えたい

Question

実現したいこと

下記画像で赤線に塗っている部分のオブジェクトを輪郭検出したい。

前提

画像中の線の数を数えるプログラムを作りたいのですが、見切れているためか検出されない線があります。
どういった手法なら以下の画像中に存在する見切れた線も検出することができるでしょうか？

※7/21追記

みなさんのコメントのおかげでモルフォロジー変換でノイズを飛ばすことができました！

モルフォロジー変換でノイズとして消すことのできなかった部分は領域の大きさで除外しようと思っていたのですが
下記のコードで領域の大きさを視覚化したところどうも赤線部分のオブジェクトは極端に領域の大きさが小さくとらえられているようです。
なにかコードに不備があるのでしょうか？
上記の画像から太線のみを検出できる方法があれば教えていただきたいです。

from google.colab.patches import cv2_imshow
import cv2
import numpy as np

# 画像を読み込む
image = cv2.imread('/content/28896_12765.png')

# 画像をグレースケールに変換
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

# Cannyエッジ検出を適用
threshold1 = 100
threshold2 = 650
edges = cv2.Canny(gray, threshold1, threshold2)

# 輪郭を検出
contours, hierarchy = cv2.findContours(edges, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

# 面積が一定以上の輪郭のみを抽出
filtered_contours = [cnt for cnt in contours if cv2.contourArea(cnt) > 2]

# 線を描画したコピーを作成
image_lines = image.copy()
cv2.drawContours(image_lines, filtered_contours, -1, (0, 255, 0), 2)

# 各輪郭の面積を描画
for cnt in filtered_contours:
    # 輪郭の面積を取得
    area = cv2.contourArea(cnt)

    # 輪郭の中心座標を取得
    M = cv2.moments(cnt)
    cx = int(M["m10"] / M["m00"])
    cy = int(M["m01"] / M["m00"])

    # テキストを描画
    cv2.putText(image_lines, f"S: {area}", (cx - 50, cy), cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0, 0, 255), 2)

# 画像を表示
cv2_imshow(image_lines)

# 検出された輪郭の数を出力
line_count = len(filtered_contours)
print("線の数:", line_count)

発生している問題・エラーメッセージ

エラー自体は出ていないのですが、上記の画像のように上手く検出できない線があります。
上から1/4ほどのところで画像を横切っている太い線や画像左下部の3つ並んだ線、画像右下から伸びる少し穴の開いた線なども検出できるようにしたいです。

from google.colab.patches import cv2_imshow
import cv2
import numpy as np

# 画像を読み込む
image = cv2.imread('/content/gazo1.png')

# 画像をグレースケールに変換
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

# Cannyエッジ検出を適用
threshold1 = 50
threshold2 = 150
edges = cv2.Canny(gray, threshold1, threshold2)

# 輪郭を検出
contours, hierarchy = cv2.findContours(edges, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)

# 輪郭の面積が40以上のもののみ抽出
filtered_contours = [cnt for cnt in contours if cv2.contourArea(cnt) > 40]

# 線を描画したコピーを作成
image_lines = image.copy()
cv2.drawContours(image_lines, filtered_contours, -1, (0, 255, 0), 2)

# 画像を表示
cv2_imshow(image_lines)

# 検出された輪郭の数を出力
line_count = len(filtered_contours)
print("線の数:", line_count)

試したこと

しきい値やthreshold値の変更

Comments

[y_waiwai]

このままではコードが読みづらいので、質問を編集し、</>(コードの挿入)ボタンを押し、出てくる’’’の枠の中にコードを貼り付けてください

[退会済みユーザー]

どの辺までがノイズで、どの辺までが線か（MSPAINTでの塗りつぶし程度でもいいので）表示できますか？
ぱっと思いつくだけで方法が２、３あり、どれがいいかが変わってくる話のためです。

[misox]

ご指摘ありがとうございます。
使い方を理解できていませんでした。
コードの貼り付け方を修正しました。

画像は赤線で塗った太い部分が検出対象です。
それ以外の細かい点や細い線は全てノイズです。

[yominet]

ノイズ消しでよく使われる手法
グレースケール画像に対し
１．erodeで、細かいをもの消す
２．dilateで、残ったものの太さを戻す

◆ついでに
赤色にぬった基準はなんでしょうか？
たとえば、太さとか長さとか？

[misox]

ありがとうございます。
なるほど、検出の前にノイズ消しの処理をした方がいいのですね。

赤色にぬった基準は長さ、というよりはまとまり具合でしょうか
これは堤防の検出を目的として、標高しきい値が一定範囲の部分を白色で表示させているものです。
赤色に塗っている部分が堤防で他の部分はたまたま、しきい値が条件に合致してしまったノイズです。

ですので赤色にぬった基準は「一定以上の長さでピクセルがまとまっている部分」です。

[fana]

話が２値化済みの画像からの出発になっている場合，まぁ「ノイズを消しましょう」とかそういう話になるのだろうと思うが，本件で扱いたい話というのはそういう事柄で良いのか否か？

何と言うか，その絵が「やりたいことに対しての良い２値化結果」とは見えないのだが…
そもそもやりたいことに対して「行った２値化自体がイマイチ」という話も多分にあり得るのではないかと個人的には思える．

現状が【「とりあえず」標高データに閾値処理をしてみた】程度の段階（確固たる何かが無いような状況）なのであれば

・処理の流れの中でいつかは「２値化」するにしても，もっとマシな結果を得る方法は無いのか？
・そもそも「標高が一定以上のところ」という単純な話「だけ」していれば良いのか？

みたいなところをまずは考えるべきでは…？　とか思ったりする．
（例えば，堤防は確かに周囲よりは高いのだろうけど，高さだけでなく高さの変化量には特徴は無いのか？　とか）
簡単に「ノイズ」という言葉で切り捨てているけども，「それらは何なの？」っていうのを検討すべきかと．
やるべきことは　それらと堤防との間の差とは何か？　というのを見つけることなのではあるまいか．
（「現在の２値画像での領域形状(?)から十分に判別できる」みたいな 理屈/検討結果/etc があっての話なのであればそもそも質問していないだろうし？）

[misox]

おっしゃる通りだと思います。
自分でもまだ検出に足る画像データではないと思っています。
しかし、検出にどういった手法があるのかも知らずに闇雲にデータを集めても有効な情報を集められるかわからなかったので一度今あるデータで検出を試み、その上でなぜこの赤線部分は上手く検出できないのかをここで質問させて頂きました。
高さの変化量ですが「標高しきい値」を同じ意味の言葉だと思って使っていました。申し訳ありません。

[fana]

変化量というのは，「なだらかに変化するのか急激に変わるのか」みたいな話です．
（高台があったとしてその側面が坂になってるのか壁になってるのか，というような）

それはそれとして，本件の画像について「なぜうまく検出できないのか」を知るには，
Cannyの結果や，findContours の結果等，各ステップの結果を可視化（表示）して見るとよいのではないかと．
ある処理段階であなたの想定とは異なる結果を生じていれば，そこが原因でしょう．

[jbpb0]

元画像も質問にアップしたら、質問のコードでの処理内容とは別のアイデアを書いてくれる人がいるかも

> # 画像をグレースケールに変換

をしてるのなら、おそらく元画像はカラーですよね
色の情報が使えないかな？？ とか思ったり
(いろんな画像で、それぞれ違うのかもしれませんが)

[fana]

実際に（←もちろん仮想的に）水を流してみたらわかったりしないのかな？
Wikipediaで「堤防」見たら何やらいろいろと種類があるみたいで，判別には相応の知見が必要そうだが…

[退会済みユーザー]

> これは堤防の検出を目的として、標高しきい値が一定範囲の部分を白色で表示させているものです。
だとすれば、ドローンから空撮した写真、もしくはLiDARか何かで撮った距離情報、もしかしたら超高解像度の衛星写真の端くれだったりしますか？
高さ情報以外にRGB情報もあるならさばきようが変わるとも思います。

[misox]

みなさんありがとうございます。
やはりこの画像だけでは検出は難しいのでしょうか？

ちなみにこの画像は国土地理院の標高タイル
https://maps.gsi.go.jp/development/demtile.html
を利用して得た標高情報を可視化させたものです。
実際に質問のコードで使っているデータ：https://cyberjapandata.gsi.go.jp/xyz/dem5a/15/28896/12765.txt
表示の参考にしたサイト：https://www.kkaneko.jp/db/map/kokudo.html
です

[misox]

RGB情報は同じく国土地理院の発行している地図データから得られるかもしれませんが複数ファイルの情報を組み合わせることもできるのでしょうか？
https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/15/28896/12765.jpg
https://cyberjapandata.gsi.go.jp/xyz/std/15/28896/12765.png

[jbpb0]

> 表示の参考にしたサイト

https://www.kkaneko.jp/db/map/kokudo.html
の
「標高が高いところが赤，低いところが青． になるように、自動的に色がつく．」
のやり方でカラー画像を作って保存し、それを質問のコード中でグレースケールに変換してる、ということでしょうか？

もしそうなら、カラー画像を作るのではなく、
https://www.kkaneko.jp/db/map/kokudo.html
のコードの
plt.imshow(Z)
を
plt.imshow(Z, cmap = "gray")
と変えて、直接グレースケール画像を作る方がいいかも

カラー画像を経由したら、元データ(csvファイル)の数値の大小と、グレースケールデータの明暗が一致しないかもしれないので

あるいは、画像を経由しないで、csvファイルを読んだ元データをnumpy配列にして、質問のコードの「gray」として使ったら、どうでしょうか
(整数に変換する必要があるかも)

【追記】
画像を経由する場合ですが、カラーではない単純な二次元データをmatplotlibのimshowで表示すると、たしかデータの最小値〜最大値の範囲に合わせて勝手に正規化されて表示されるので、それを使って傾斜量を計算しても、別の場所のデータ間で値の大小が一致しないと思います

勝手に正規化されず、全ての場所のデータで傾斜量の大小が一致するように、
https://www.kkaneko.jp/db/map/kokudo.html
のコードの「plt.imshow(...」に、「vmin」と「vmax」を追加するといいと思います

[misox]

なるほど！自分の見やすさしか考えていませんでしたが確かに直接グレースケール画像を作った方がいいですね。
「vmin」「vmax」も知りませんでした、ありがとうございます。

[jbpb0]

> 複数ファイルの情報を組み合わせることもできるのでしょうか？

複数の画像間の、実際の地理的な座標の関係が分かる(A画像のある画素と同じ場所が、B画像ではどこの画素なのかが分かる)のであれば、できますよ

それぞれの画像から「堤防らしさ」をそれぞれ算出して、それぞれの「堤防らしさ」の同じ場所(実際の地理的な座標)の値を組み合わせて、判定します
組み合わせ方は、たとえば
・どちらも閾値を超えた場所だけ抽出する
・掛け算する
・足し算する
とか、いろいろ考えられます

https://cyberjapandata.gsi.go.jp/xyz/std/15/28896/12765.png
を見ると、川は水色っぽい色で表示されてますよね

質問者さんが抽出したい「堤防」の定義次第ですが、近くに川や海がある場合に限定されるのなら、近くに川や海の色があるかどうかも「堤防らしさ」の算出に使えば、誤判定が減らせるような気がします

[退会済みユーザー]

> 複数ファイルの情報を組み合わせることもできるのでしょうか？
こういうのをマスク処理と言います。あとは、国土地理院の写真で高さ情報と、RGBの同期されたデータがあるので、あとはどこまで気合を入れて解決したいか、です。

ガチの地図データばりの精度を目指すのであれば、MRCNNだとかpix2pixとかで探すと、海外勢で「地図データから建物を抽出してみた」系の記事があるハズです。
多少の誤差は気にしない、であれば、モルフォロジー変換でカスを飛ばしてその後にscipyのスケルトンだとか分水嶺だとかを駆使してっぽい部分を抽出してからハフ検出でうまく行くハズです。

あとは、どの方向性でどのくらいの精度でやりたいかを質問に加えると、誰かの引き出しにしっかりハマって回答をもらいやすくなります。

[misox]

なるほど、ありがとうございます！
ゆくゆくはガチの地図データばりの精度も目指してみたいのですが
今回は多少の誤差は気にしない程度の検出を考えています。
scipyというものがあったのですね！まさにこういった処理をしたいと考えていました。
モルフォロジー変換のノイズ除去はできたのですが明線化のやり方が分からないのでよければ教えていただけないでしょうか？

Answer 1

閾値調整はきちんとしていませんが雰囲気は出ていると思います。

生	スケルトン	モルフォロジー変換	ハフ検出

import numpy as np
import cv2
from skimage.morphology import skeletonize

img = cv2.imread('test.png')

# グレースケールがいい
if len(img.shape) ==3:
    img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

# 二値化したい
img[img<128] = 0
img[img!=0] = 255

# さばきたいイメージ
cv2.imshow('img',img )

# スケルトン化した後扱いやすいようにuint8にする
img_skeleton = (skeletonize(img)*255).astype(np.uint8)

# とりあえず骨にはなった
cv2.imshow('img_skeleton',img_skeleton )
cv2.imwrite('img_skeleton.png',img_skeleton)


bool_bypass = False
 
if bool_bypass:
    img_bold = img_skeleton
else:
    # 多分線が細すぎるので太くする
    # http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_morphological_ops/py_morphological_ops.html
    kernel = np.ones((3,3),np.uint8)
    img_bold = cv2.dilate(img_skeleton,kernel,iterations = 2)
    img_bold = cv2.morphologyEx(img_bold, cv2.MORPH_CLOSE, kernel,iterations = 5)
    cv2.imshow('img_bold',img_bold )

# ハフ検出したい
# http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_houghlines/py_houghlines.html
rho = 1
theta = np.pi/180
min_vote = 0
minLineLength = 10
maxLineGap = 10
lines = cv2.HoughLinesP(img_bold,rho,theta ,min_vote ,minLineLength,maxLineGap)

if lines is None:
    print('線不検出')
    cv2.waitKey(1000)
else:
    print('検出した線:',len(lines))
    # 描画用にRGBに戻す
    img_bold = cv2.cvtColor(img_bold,cv2.COLOR_GRAY2BGR)
     
    for line in lines:
        x1,y1,x2,y2 = line[0]
        print('デバグ座標群:',x1,y1,x2,y2)
        img_bold = cv2.line(img_bold ,(x1,y1),(x2,y2),(0,255,0),3)

    # cv2.imwrite('houghlines.jpg',img_bold )
    cv2.imshow('img_hough',img_bold )
    cv2.imwrite('houghlines.png',img_bold )

    cv2.waitKey(0)

Comments

[退会済みユーザー]

モルフォロジー変換でごみとりはしたほうがいいですがスケルトンはなくてもよかったかもしれません。
特に、HoughLinesPは「途切れていても適当に補って検出」「長さが短ければ無視」など相当柔軟に対応するので、もしかしたら上の画像でいう二値化しただけの「生」を食わせるのが一番いい、もあり得ます。

この辺は検討の余地ありです。あとは何とかなりそうですか:)

[misox]

はい！ありがとうございます。
無事ある程度の検出はできるようになりました。
あとは自分でいろいろ微調整していこうと思います。