LoginSignup
1
0

More than 1 year has passed since last update.

@aa_debdeb(Atsushi Asakura)

【画像処理】NumpyでCannyエッジ検出

Last updated at Posted at 2021-12-26

NumpyでCannyエッジ検出を実装してみます。

画像読み込み

画像を読み込み、カラー画像をグレースケール画像に変換します。

original_image = plt.imread(image_name)
if np.issubdtype(original_image.dtype, np.integer):
  original_image = original_image / np.iinfo(original_image.dtype).max
gray_image = 0.2116 * original_image[:,:,0] + 0.7152 * original_image[:,:,1] + 0.0722 * original_image[:,:,2]
plt.imshow(gray_image, cmap='gray')

GrayImage.png

畳み込み関数

畳み込みを行う関数convolve2dを定義しておきます。詳細は以下の記事を参照してください。
【画像処理】Numpyで空間フィルタリング(畳み込み演算) - Qiita

def _convolve2d(image, kernel):
  shape = (image.shape[0] - kernel.shape[0] + 1, image.shape[1] - kernel.shape[1] + 1) + kernel.shape
  strides = image.strides * 2
  strided_image = np.lib.stride_tricks.as_strided(image, shape, strides)
  return np.einsum('kl,ijkl->ij', kernel, strided_image)

def _convolve2d_multichannel(image, kernel):
  convolved_image = np.empty((image.shape[0] - kernel.shape[0] + 1, image.shape[1] - kernel.shape[1] + 1, image.shape[2]))
  for i in range(image.shape[2]):
    convolved_image[:,:,i] = _convolve2d(image[:,:,i], kernel)
  return convolved_image

def _pad_singlechannel_image(image, kernel_shape, boundary):
  return np.pad(image, ((int(kernel_shape[0] / 2),), (int(kernel_shape[1] / 2),)), boundary)

def _pad_multichannel_image(image, kernel_shape, boundary):
  return  np.pad(image, ((int(kernel_shape[0] / 2),), (int(kernel_shape[1] / 2),), (0,)), boundary)

def convolve2d(image, kernel, boundary='edge'):
  if image.ndim == 2:
    pad_image = _pad_singlechannel_image(image, kernel.shape, boundary) if boundary is not None else image
    return _convolve2d(pad_image, kernel)
  elif image.ndim == 3:
    pad_image = _pad_multichannel_image(image, kernel.shape, boundary) if boundary is not None else image
    return _convolve2d_multichannel(pad_image, kernel)

画像の平滑化(ノイズ除去)

ガウシアンフィルターで画像を平滑化して、ノイズを除去します。

def create_gaussian_kernel(size=(3, 3), sigma=1):
  center = ((size[0] - 1) / 2, (size[1] - 1) / 2)
  sigma2 = 2 * sigma * sigma
  kernel = np.fromfunction(lambda y, x: np.exp(-((x - center[1]) ** 2 + (y - center[0]) ** 2) / sigma2), size)
  kernel = kernel / np.sum(kernel)
  return kernel 

gaussian_kernel = create_gaussian_kernel()
gaussian_image = convolve2d(gray_image, gaussian_kernel)
plt.imshow(gaussian_image, cmap='gray')

GaussianImage.png

画像の微分

Sobelフィルターで画像の微分を行い、勾配を求めます。

まず横方向の微分です。

sobel_x_kernel = np.array([[-1, 0, 1],
                           [-2, 0, 2],
                           [-1, 0, 1]])
sobel_x_image = convolve2d(gaussian_image, sobel_x_kernel) / 8
value_range = max(abs(sobel_x_image.min()), abs(sobel_x_image.max()))
plt.imshow(sobel_x_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()

SobelXImage.png

次に縦方向の微分です。

sobel_y_kernel = np.array(([[1, 2, 1],
                            [0, 0, 0],
                            [-1, -2, -1]]))
sobel_y_image = convolve2d(gaussian_image, sobel_y_kernel) / 8
value_range = max(abs(sobel_y_image.min()), abs(sobel_y_image.max()))
plt.imshow(sobel_y_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()

SobelYImage.png

各方向の微分から勾配の強度と方向を求めます。

sobel_intensity_image = np.sqrt(sobel_x_image ** 2 + sobel_y_image ** 2)
plt.imshow(sobel_intensity_image, cmap='gray')
plt.colorbar()

SobelIntensityImage.png

sobel_direction_image = np.arctan(np.divide(sobel_y_image, sobel_x_image, out=np.zeros_like(sobel_y_image), where=sobel_x_image!=0))
plt.imshow(sobel_direction_image, cmap='hsv')
plt.colorbar()

SobelDirectionImage.png

勾配方向の非極値抑制

細線化のために、先ほど求めた勾配強度のうち勾配方向で最大値となるものだけを残します。

def _suppress_nonmaximum(array):
  direction = array[9]
  if direction < 0:
    direction += np.pi
  if direction < 0.125 * np.pi:
    v0, v1, v2 = (array[3], array[4], array[5])
  elif direction < 0.375 * np.pi:
    v0, v1, v2 = (array[2], array[4], array[6])
  elif direction < 0.625 * np.pi:
    v0, v1, v2 = (array[1], array[4], array[7])
  elif direction < 0.875 * np.pi:
    v0, v1, v2 = (array[0], array[4], array[8])
  else:
    v0, v1, v2 = (array[3], array[4], array[5])
  return v1 if v1 >= v0 and v1 >= v2 else 0

def suppress_nonmaximum(intensity_image, direction_image):
  pad_image = np.pad(intensity_image, 1, 'edge')
  shape = intensity_image.shape + (3, 3)
  strides = pad_image.strides * 2
  strided_image = np.lib.stride_tricks.as_strided(pad_image, shape, strides).reshape(shape[0], shape[1], 9)
  return np.apply_along_axis(_suppress_nonmaximum, 2, np.dstack((strided_image, direction_image)))

nonmaximum_suppression_image = suppress_nonmaximum(sobel_intensity_image, sobel_direction_image)
plt.imshow(nonmaximum_suppression_image, cmap='gray')
plt.colorbar()

SuppressNonmaximumImage.png

ヒステリシス閾値処理

非極値抑制した画像からエッジ画像を作成します。まず、値が閾値$T_{high}$よりも大きい画素をエッジとみなします。さらに値が閾値$T_{low}$よりも大きく、かつエッジが隣接している画素もエッジとみなすことで連結したエッジを検出します。

def histeresis_thresholding(image, threshold):
  edge_image = np.where(image >= threshold[1], 1.0, 0.0)

  while True:
    pad_edge_image = np.pad(edge_image, 1, 'edge')
    shape = image.shape + (3, 3)
    strides = pad_edge_image.strides * 2
    neighbour_edge_image = np.lib.stride_tricks.as_strided(pad_edge_image, shape, strides)
    edge_candidate_index = np.where((image >= threshold[0]) & (edge_image != 1.0))

    changed = False
    for index in list(zip(*edge_candidate_index)):
      if np.any(neighbour_edge_image[index] == 1.0):
        edge_image[index] = 1.0
        changed = True
    if not changed:
      break;

  return edge_image

canny_image = histeresis_thresholding(nonmaximum_suppression_image, (0.03, 0.1))
plt.imshow(canny_image, cmap='gray')

HisteresisThresholdingImage.png

まとめ

以上の処理をまとめて、Cannyエッジ検出を行う関数を作成すると次のようになります。

def canny(image, gaussian_size=(3, 3), gaussian_sigma=1, threshold=(0.02, 0.05)):
  # Noise reduction by gaussian filter
  gaussian_kernel = create_gaussian_kernel(size=gaussian_size, sigma=gaussian_sigma)
  gaussian_image = convolve2d(gray_image, gaussian_kernel)

  # Gradient by sobel filter
  sobel_x_kernel = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
  sobel_x_image = convolve2d(gaussian_image, sobel_x_kernel) / 8
  sobel_y_kernel = np.array(([[1, 2, 1], [0, 0, 0], [-1, -2, -1]]))
  sobel_y_image = convolve2d(gaussian_image, sobel_y_kernel) / 8
  sobel_intensity_image = np.sqrt(sobel_x_image ** 2 + sobel_y_image ** 2)
  sobel_direction_image = np.arctan(np.divide(sobel_y_image, sobel_x_image, out=np.zeros_like(sobel_y_image), where=sobel_x_image!=0))

  # Non-maximum suppression
  nonmaximum_suppression_image = suppress_nonmaximum(sobel_intensity_image, sobel_direction_image)

  # Histresis thresholding
  return histeresis_thresholding(nonmaximum_suppression_image, threshold)

パラメータを変えてCannyエッジ検出を試してみます。

canny_image1 = canny(gray_image)
plt.imshow(canny_image1, cmap='gray')

CannyImage1.png

canny_image2 = canny(gray_image, gaussian_size=(5, 5), gaussian_sigma=0.2, threshold=(0.01, 0.04))
plt.imshow(canny_image2, cmap='gray')

CannyImage2.png

canny_image3 = canny(gray_image, gaussian_size=(3, 3), gaussian_sigma=1.2, threshold=(0.06, 0.08))
plt.imshow(canny_image3, cmap='gray')

CannyImage3.png

今回実装したコードはGoogle Colaboratoryに置いてあります。

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0