Stanford University/CS131/Homework 1の最後の部分をやる。これでやっと長かった宿題１とお別れできる。宿題１でこれだから、まじで先が思いやられる。

# Setup
import numpy as np
import matplotlib.pyplot as plt
from time import time
from skimage import io

from __future__ import print_function

%matplotlib inline
plt.rcParams['figure.figsize'] = 18.0, 14.0 # set default size of plots
plt.rcParams["font.size"] = "18"
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'

Part 3: Separable Filters¶

Theory¶

Consider a $M_1\times{N_1}$ image $I$ and a $M_2\times{N_2}$ filter $F$. A filter $F$ is separable if it can be written as a product of two 1D filters: $F=F_1F_2$.
$M_1\times{N_1}$ image $I$、$M_2\times{N_2}$ filter $F$と考え、filterが2個の1次元フィルタの積として $F=F_1F_2$になる場合、$F$は可分である。
For example,
$$F=
\begin{bmatrix}
1 & -1 \\
1 & -1
\end{bmatrix}
$$
can be written as a matrix product of
以下の行列積として書き表せる。
$$F_1=
\begin{bmatrix}
1 \\
1
\end{bmatrix},
F_2=
\begin{bmatrix}
1 & -1
\end{bmatrix}
$$
Therefore $F$ is a separable filter.
従って、$F$は可分フィルタである。
Prove that for any separable filter $F=F_1F_2$,
全ての可分フィルタ$F=F_1F_2$に対し以下の式が成立することを証明せよ。
$$I*F=(I*F_1)*F_2$$

Since F is separable, so
Fは可分なので
$$F[i,j] = F_1[i]F_2[j]$$
以下のようになる
\begin{align}
(I*F)[m,n] &= \sum_{i=-\infty}^\infty\sum_{j=-\infty}^\infty I[i,j]\cdot F[m-i,n-j]\\
&= \sum_{i=-\infty}^\infty\sum_{j=-\infty}^\infty I[i,j]\cdot F_1[m-i]\cdot F_2[n-j]\\
&= \sum_{j=-\infty}^\infty F_2[n-j]\sum_{i=-\infty}^\infty I[i,j]\cdot F_1[m-i]\\
&= \sum_{j=-\infty}^\infty F_2[n-j]\cdot(I*F_1)\\
&= (I*F_1)*F_2
\end{align}

Complexity comparison¶

(i) How many multiplications do you need to do a direct 2D convolution (i.e. $I*F$?)
(ii) How many multiplications do you need to do 1D convolutions on rows and columns (i.e. $(I*F_1)*F_2$)
(iii) Use Big-O notation to argue which one is more efficient in general: direct 2D convolution or two successive 1D convolutions?

(i)2次元直接畳み込み(例えば$I*F$)に必要な乗数は？
(ii)行と列に対する1次元畳み込み(例えば$(I*F_1)*F_2$)に必要な乗数は？
(iii)2次元直接畳み込みと、2連続の1次元畳み込みのどちらがより効率的であるかを、Big-O表記法を使って論ぜよ。

(i) $H_i*W_i*H_k*W_k$
(ii) $H_i*W_i*H_k+W_k*H_i*W_i$
(iii) direct 2D convolution is $O(n^3)$, two successive 1D convolutions is $O(n^2)$
(iii)2次元直接畳み込みは$O(n^3)$で、連続1次元畳み込みは$O(n^2)$

Now, we will empirically compare the running time of a separable 2D convolution and its equivalent two 1D convolutions. Gaussian kernel, widely used for blurring images, is one example of a separable filter. Run the code below to see its effect.
次に、可分2D畳み込みと、それに対応する連続1次元畳み込みの実行時間を実験的に比較する。ぼかし画像に広く使われているガウス核は、可分フィルタの1つの例だ。下のコードを実行して効果を確かめよ。

from filters import conv_nested
# Load image
img = io.imread('dog.jpg', as_gray=True)

# 5x5 Gaussian blur
kernel = np.array(
[
    [1,4,6,4,1],
    [4,16,24,16,4],
    [6,24,36,24,6],
    [4,16,24,16,4],
    [1,4,6,4,1]
])

t0 = time()
out = conv_nested(img, kernel)
t1 = time()
t_normal = t1 - t0

# Plot original image
plt.subplot(1,2,1)
plt.imshow(img)
plt.title('Original')
plt.axis('off')

# Plot convolved image
plt.subplot(1,2,2)
plt.imshow(out)
plt.title('Blurred')
plt.axis('off')

plt.show()

In the below code cell, define the two 1D arrays (k1 and k2) whose product is equal to the Gaussian kernel.
下のコードセルに、積がガウス核に等しい2つの1次元配列(k1とk2)を定義せよ。

# The kernel can be written as outer product of two 1D filters
k1 = None  # shape (5, 1)
k2 = None  # shape (1, 5)

### YOUR CODE HERE
k1 = np.array([1,4,6,4,1]).reshape((5,1))
k2 = np.array([1,4,6,4,1]).reshape((1,5))
### END YOUR CODE

# Check if kernel is product of k1 and k2
if not  np.all(k1 * k2 == kernel):
    print('k1 * k2 is not equal to kernel')
    
assert k1.shape == (5, 1), "k1 should have shape (5, 1)"
assert k2.shape == (1, 5), "k2 should have shape (1, 5)"

We now apply the two versions of convolution to the same image, and compare their running time. Note that the outputs of the two convolutions must be the same.
次に、同じ画像に2つのバージョンの畳み込みを適用して実行時間を比較する。2つの畳み込みの出力は同じでなければならないことに留意する。

# Perform two convolutions using k1 and k2
t0 = time()
out_separable = conv_nested(img, k1)
out_separable = conv_nested(out_separable, k2)
t1 = time()
t_separable = t1 - t0

# Plot normal convolution image
plt.subplot(1,2,1)
plt.imshow(out)
plt.title('Normal convolution')
plt.axis('off')

# Plot separable convolution image
plt.subplot(1,2,2)
plt.imshow(out_separable)
plt.title('Separable convolution')
plt.axis('off')

plt.show()

print("Normal convolution: took %f seconds." % (t_normal))
print("Separable convolution: took %f seconds." % (t_separable))

Normal convolution: took 6.983799 seconds.
Separable convolution: took 2.992070 seconds.

# Check if the two outputs are equal
assert np.max(out_separable - out) < 1e-10

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-45-97f743632d95> in <module>()
      1 # Check if the two outputs are equal
----> 2 assert np.max(out_separable1 - out) < 1e-10

AssertionError:

何故か2つの出力は等しくならなかった。conv_nestedの実装に問題があるのかもしれないので、一応調べてみることにした。

def conv_nested(image, kernel):
    """A naive implementation of convolution filter.
    This is a naive implementation of convolution using 4 nested for-loops.
    This function computes convolution of an image with a kernel and outputs
    the result that has the same shape as the input image.
    Args:
        image: numpy array of shape (Hi, Wi)
        kernel: numpy array of shape (Hk, Wk)
    Returns:
        out: numpy array of shape (Hi, Wi)
    """
    Hi, Wi = image.shape
    Hk, Wk = kernel.shape
    out = np.zeros((Hi, Wi))
    ### YOUR CODE HERE
    # for-loopsを4個作る
    for m in range(Hi): # hight(高さ)※-1がないとエラー
        for n in range(Wi): # width(横幅)※-1がないとエラー
            for i in range(Hk): # Hkernel(Kernel高さ)
                for j in range(Wk): # Wkernel(kernel幅)
                    # image=f, kernel=h ※+1がないとエラー
                    out[m,n]+=image[m-i][n-j]*kernel[i,j]
    ### END YOUR CODE
    return out

# Load image
img = io.imread('dog.jpg', as_gray=True)

# 5x5 Gaussian blur
kernel = np.array(
[
    [1,4,6,4,1],
    [4,16,24,16,4],
    [6,24,36,24,6],
    [4,16,24,16,4],
    [1,4,6,4,1]
])

t0 = time()
out = conv_nested(img, kernel)
t1 = time()
t_normal = t1 - t0

# Plot original image
plt.subplot(1,2,1)
plt.imshow(img)
plt.title('Original')
plt.axis('off')

# Plot convolved image
plt.subplot(1,2,2)
plt.imshow(out)
plt.title('Blurred')
plt.axis('off')

(-0.5, 331.5, 299.5, -0.5)

# The kernel can be written as outer product of two 1D filters
k1 = None  # shape (5, 1)
k2 = None  # shape (1, 5)

### YOUR CODE HERE
k1 = np.array([1,4,6,4,1]).reshape((5,1))
k2 = np.array([1,4,6,4,1]).reshape((1,5))
### END YOUR CODE

# Check if kernel is product of k1 and k2
if not  np.all(k1 * k2 == kernel):
    print('k1 * k2 is not equal to kernel')
    
assert k1.shape == (5, 1), "k1 should have shape (5, 1)"
assert k2.shape == (1, 5), "k2 should have shape (1, 5)"

# Perform two convolutions using k1 and k2
t0 = time()
out_separable = conv_nested(img, k1)
out_separable = conv_nested(out_separable, k2)
t1 = time()
t_separable = t1 - t0
# Plot normal convolution image
plt.subplot(1,2,1)
plt.imshow(out)
plt.title('Normal convolution')
plt.axis('off')
# Plot separable convolution image
plt.subplot(1,2,2)
plt.imshow(out_separable)
plt.title('Separable convolution')
plt.axis('off')
plt.show()
print("Normal convolution: took %f seconds." % (t_normal))
print("Separable convolution: took %f seconds." % (t_separable))

Normal convolution: took 1.828289 seconds.
Separable convolution: took 0.871541 seconds.

# Check if the two outputs are equal
assert np.max(out_separable - out) < 1e-10

やはりconv_nestedの実装に誤りがあったようだ。チョロっとコードをいじっただけで実行速度が3.82倍も向上している。しかし、HW1でこれだと、freshmanでCS131を履修しなければならないとなると、かなりしんどいだろうなぁと思われる。アメリカの大学の宿題はとにかく半端ない。例えば、18クレジット(4クレジット✕3、3クレジット✕2)を選択した場合、宿題の量はまじでやばい。外国人留学生の場合、原則として1セメスター12クレジットに規定されていたが、俺は貧乏で一刻も早く卒業しないといけなかったので18クレジットを選択して人生詰んだ。金さえあればと今でも悔やまれる。