2019年5月27日 更新

Extracting Straight Lines〜画像から境界線を効率よく求める方法〜

本記事では,IEEE Transactions on Pattern Analysis and Machine Intelligence に投稿された,「Extracting Straight Lines」の論文を説明し,実装を行います.1986年に投稿された少し古めの論文ですが,アイディアは非常に面白いです.

1,172 view お気に入り 0

Line Support

次の段階として,このGradient の向きが同じ方向を向いている領域を探します. 境界線付近のGradient は,境界線に垂直な方向に揃っているはずです. 向きの揃っている領域,というのがLine Support です.

LSD: a Line Segment Detecto...

LSD: a Line Segment Detector [2] からの引用.

LSDのアルゴリズムでは,このLine Support 領域がある閾値よりも大きい場合,それを境界として認識します.

本記事では,Line Support 領域を求め,Line Support 領域を連結していくことにより,エッジを検出します.

実装

import numpy as np
import matplotlib.pyplot as plt
import cv2
import os
import scipy as sp
from scipy import signal
import skimage
import sys
pi = np.pi
01.py

画像の作成

Is = np.zeros([50,80])
shp = np.shape(Is)

for i in range(shp[0]):
    for j in range(int(shp[1]/2)):
        Is[i,j] = 125

def addGaussianNoise(src):
    [row,col] = np.shape(src)
    mean = 0
    var = 5.5
    sigma = 135
    gauss = np.random.normal(mean,sigma,(row,col))
    gauss = gauss.reshape(row,col)
    noisy = src + gauss

    return noisy

Is = Is + addGaussianNoise(Is)


## Normalize

Is = np.round(255*(Is - np.min(Is))/(np.max(Is) - np.min(Is)))
02.py

Gradient の計算

kernel_y = np.array([[-1, -1],
                   [1, 1]])
kernel_x = np.array([[-1, 1],
                   [-1, 1]])


Is_grady = cv2.filter2D(src=Is, kernel=kernel_y, ddepth=-1)
Is_gradx = cv2.filter2D(src=Is, kernel=kernel_x, ddepth=-1)


import matplotlib.patches as pat



plt.figure()
ax = plt.subplot(1,1,1)
plt.quiver(Is_gradx,Is_grady,angles="xy",headwidth=3,color="#444444")
rec1 = pat.Rectangle(xy = (35, 15), width = 10, height = 10, 
                     angle = 0, color = "blue", alpha = 0.5)
ax.add_patch(rec1)
ax.set_aspect('equal')


plt.figure()
ax2 = plt.subplot(1,1,1)
plt.quiver(Is_gradx[15:25,35:45],Is_grady[15:25,35:45],angles="xy",headwidth=3,color="#444444")
ax2.set_aspect('equal')
03.py

Gradient Orientation を16方向へ分類

Is_grad_norm = np.sqrt(Is_gradx**2 + Is_grady**2)

## Parameter1
TOL_GRAD = 100

Is_grad_theta = np.ones(shp)*pi/2
Is_grad_theta[Is_gradx != 0] = np.arctan(Is_grady[Is_gradx != 0]/Is_gradx[Is_gradx != 0])
Is_grad_theta = Is_grad_theta * (Is_grad_norm > Thre_grad)

## 
theta = pi/8
theta_all = np.arange(-pi,pi,theta)

Is_grad_cls = -pi*np.ones(shp)
Is_grad_theta_k = np.ones(shp) * sys.maxsize
ones = np.ones(shp)
for theta_i in theta_all:
    abs_grad_thetai = np.abs(Is_grad_theta - theta_i)
    grad_cls =  abs_grad_thetai < Is_grad_theta_k
    Is_grad_cls[grad_cls] = theta_i*ones[grad_cls]
    Is_grad_theta_k[grad_cls] = abs_grad_thetai[grad_cls]


tg = Is_grad_norm > Thre_grad
Is_grad_cls = Is_grad_cls * tg
04.py

Line Support & 境界検出

TOL_ANG = pi/6
theta_ch = [-pi,-pi/2,0,pi/2]

TOL_PIX = np.min(shp)*0.3

for theta_i in theta_ch:
    Is_grad_region = np.logical_and(Is_grad_cls > (theta_i - TOL_ANG),Is_grad_cls < (theta_i + TOL_ANG))
    Is_grad_region = (Is_grad_region * tg + 0).astype("uint8")
    
    plt.figure()
    plt.imshow(Is,cmap="gray")
    plt.title(str(theta_i))
    
    if theta_i in [-pi,0]:
        Is_grad_region_TOL = np.sum(Is_grad_region,axis=0) > TOL_PIX
        if any(Is_grad_region_TOL) == True:
            Is_grad_region_TOL_ind = np.where(Is_grad_region_TOL)[0]
            for Is_grad_region_TOL_i in Is_grad_region_TOL_ind:
                x1 = Is_grad_region_TOL_i
                x2 = Is_grad_region_TOL_i
                y1 = np.min(np.where([Is_grad_region[:,Is_grad_region_TOL_i]][0]))
                y2 = np.max(np.where([Is_grad_region[:,Is_grad_region_TOL_i]][0]))
                plt.plot([x1,x2],[y1,y2],'--r')
    
    if theta_i in [-pi/2,pi/2]:
        Is_grad_region_TOL = np.sum(Is_grad_region,axis=1) > TOL_PIX
        if any(Is_grad_region_TOL) == True:
            Is_grad_region_TOL_ind = np.where(Is_grad_region_TOL)[0]
            for Is_grad_region_TOL_i in Is_grad_region_TOL_ind:
                x1 = Is_grad_region_TOL_i
                x2 = Is_grad_region_TOL_i
                y1 = np.min(np.where([Is_grad_region[Is_grad_region_TOL_i,:]][0]))
                y2 = np.max(np.where([Is_grad_region[Is_grad_region_TOL_i,:]][0]))
                plt.plot([x1,x2],[y1,y2],'--k')
        

05.py
Output

Output

33 件

関連する記事 こんな記事も人気です♪

ディジタル画像処理~pythonによる空間フィルタリングpart1~ 

ディジタル画像処理~pythonによる空間フィルタリングpart1~ 

ディジタル画像処理を解説します.今回は,代表的な空間フィルタリングをpythonで実行してみました。
亀谷 桃子 | 2,132 view
Morphology (モルフォロジー) 変換の実装 ~ Python + OpenCV ~

Morphology (モルフォロジー) 変換の実装 ~ Python + OpenCV ~

画像処理の一つ,モルフォロジー変換をPython と OpenCVのライブラリを用いて実装し,それを2値画像に対して適用します.
三好 裕之 | 2,362 view
等角写像による画像の変換〜Schwarz-Christoffel 変換〜part 2

等角写像による画像の変換〜Schwarz-Christoffel 変換〜part 2

前回の記事「等角写像による画像の変換〜Schwarz-Christoffel 変換〜part 1」の続きです. 実際に実装をして,写像を確かめてみます.
三好 裕之 | 1,819 view
スパースモデリング ~少量データから画像を復元~

スパースモデリング ~少量データから画像を復元~

matlab のコードをもとに,スパースモデリングとは何かについて説明します.今回はスパースモデリングの説明です.
三好 裕之 | 13,124 view
点像分布関数の話

点像分布関数の話

光学系を理解するうえで前提として知っておいた方が良い内容
Yoshiyuki Arai | 1,255 view

この記事のキーワード

この記事のキュレーター

三好 裕之 三好 裕之