NumpyでLoG(Laplacian of Gaussian)フィルターを実装して、ゼロ交差によるエッジ検出を試してみます。
LoGフィルターとは名前の通りガウス関数のラプラシアン(二次微分)で、以下の式で表されます。
LoG(x,y)=\frac{x^2+y^2-2\sigma^2}{2\pi\sigma^6}exp(-\frac{x^2+y^2}{2\sigma^2})
まず、使用する画像を読み込んで、カラー画像をグレースケール画像に変換します。
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')
画像の畳み込み演算を行う関数を定義します。詳細は以下の記事を参照してください。
【画像処理】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)
LoGフィルターのカーネルを作成して、先ほどのグレースケール画像に畳み込みます。
def create_log_kernel(size=(13, 13), sigma=2):
center = ((size[0] - 1) / 2, (size[1] - 1) / 2)
sigma2 = sigma * sigma
sigma6 = sigma2 * sigma2 * sigma2
sigma2_2 = sigma2 * 2.0
sigma2_2_inv = 1.0 / sigma2_2
sigma6_2pi_inv = 1.0 / (sigma6 * 2.0 * np.pi)
def calc_weight(y, x):
sqDist = (x - center[1]) ** 2 + (y - center[0]) ** 2
return (sqDist - sigma2_2) * sigma6_2pi_inv * np.exp(-sqDist * sigma2_2_inv)
return np.fromfunction(calc_weight, size)
log_kernel = create_log_kernel()
log_image = convolve2d(gray_image, log_kernel)
value_range = max(abs(log_image.min()), abs(log_image.max()))
plt.imshow(log_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
LoGフィルターを畳み込んだ画像からゼロ交差を求めて、エッジを検出します。ここではある画素について左右または上下の画素の正負が異なるときに、その画素がゼロ交差とみなしています。
def detect_zero_crossing(image):
pad_image = np.pad(image, 1, 'edge')
shape = 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(
lambda array: 1 if array[3] * array[5] < 0 or array[2] * array[6] < 0 else 0,
2, strided_image)
zero_crossing_image = detect_zero_crossing(log_image)
plt.imshow(zero_crossing_image, cmap='gray')
ゼロ交差までできたので、いくつかパラメータを変更して試してみます。基本的には、sigmaが大きくなるほど、sizeも大きくする必要があります。
log_kernel1 = create_log_kernel(size=(7, 7), sigma=1)
log_image1 = convolve2d(gray_image, log_kernel1)
value_range = max(abs(log_image1.min()), abs(log_image1.max()))
plt.imshow(log_image1, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
zero_crossing_image1 = detect_zero_crossing(log_image1)
plt.imshow(zero_crossing_image1, cmap='gray')
log_kernel2 = create_log_kernel(size=(31, 31), sigma=4)
log_image2 = convolve2d(gray_image, log_kernel2)
value_range = max(abs(log_image2.min()), abs(log_image2.max()))
plt.imshow(log_image2, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
zero_crossing_image2 = detect_zero_crossing(log_image2)
plt.imshow(zero_crossing_image2, cmap='gray')
log_kernel3 = create_log_kernel(size=(41, 41), sigma=7)
log_image3 = convolve2d(gray_image, log_kernel3)
value_range = max(abs(log_image3.min()), abs(log_image3.max()))
plt.imshow(log_image3, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
zero_crossing_image3 = detect_zero_crossing(log_image3)
plt.imshow(zero_crossing_image3, cmap='gray')
実装したコードはGoogle Colabratoryに置いておきました。