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

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

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)))

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')

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

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')