はじめに

こちらの参考書の第二章、Fourier and Wavelet Transformsの個人的備忘として投稿。
本稿は2.7節のTwo-Dimensional Transforms and Image Processingについて記述する。

概要

前回記事までは、1次元の信号に対してフーリエ変換とウェーブレット変換を考えてきたが、これらは簡単により高次元空間へと拡張できる。例えば、2次元変換や3次元変換などがそれに勿論該当するのだが、これらは画像処理や圧縮の分野でよく使われる手法となっている。

画像データに対する2次元フーリエ変換

画像データ$\boldsymbol{X}\in\mathbb{R}^{n\times m}$に対する2次元フーリエ変換は、まず全ての行に対し、1次元フーリエ変換を施した後に、中間行列の全列に対して1次元フーリエ変換を実施することで実現できる。

この一連の変換を以下の図に例として示す。まずは適当な画像を読み込んで、グレースケールとして描画

# Import libraries
from matplotlib.image import imread
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams.update({'font.size': 18})

# Read data & check to display
fig = imread('../data/banta.jpg')
fig = np.mean(fig, -1)
plt.imshow(fig)
plt.show()

これは星野リゾートが運営する、沖縄県読谷村の西海岸に位置するバンタカフェで撮影した1枚である。
続いて、この行列に対して行方向のフーリエ変換を適用

# Compute row-wise FFT
row_fft_shift = np.zeros_like(fig, dtype='complex_')
row_fft = np.zeros_like(fig, dtype='complex_')

for j in range(fig.shape[0]):
    row_fft_shift[j,:] = np.fft.fftshift(np.fft.fft(fig[j,:]))
    row_fft[j,:] = np.fft.fft(fig[j,:])

plt.imshow(np.log(np.abs(row_fft_shift)))
plt.show()

さらに、この中間行列に対して列方向のフーリエ変換を適用

# Compute column-wise FFT
col_fft = np.zeros_like(row_fft)
for j in range(row_fft.shape[1]):
    col_fft[:,j] = np.fft.fft(row_fft[:,j])

plt.imshow(np.fft.fftshift(np.log(np.abs(col_fft))))
plt.show()

このような手順を踏まなくとも、numpyにはnp.fft.fft2というメソッドが用意されており、一発でこれと同じ結果を得る事が出来る。

# Much more efficient to use fft2
fft2 = np.fft.fft2(fig)

plt.imshow(np.fft.fftshift(np.log(np.abs(fft2))))
plt.show()

フーリエ変換後の可視化画像はログスケールとしているが、要はここで小さい係数（=ノイズと見なせる周波数帯）を落とすことで、画質を損なうことなく圧縮をすることが出来る。

# Sort by magnitude
fft2sort = np.sort(np.abs(fft2.reshape(-1)))

# Zero out all small coefficients and inverse transform
for keep in [0.05, 0.01, 0.002]:
    # Find small indices
    thresh = fft2sort[int(np.floor((1-keep) * len(fft2sort)))]
    ind = np.abs(fft2) > thresh
    
    # Threshold small indices
    tlow = fft2 * ind
    
    # Compressed image
    low = np.fft.ifft2(tlow).real
    
    # Show image
    plt.imshow(low)
    plt.title('Compressed image: keep = ' + str(keep * 100) + '% of FFT')
    
    # Add text as a legend in the upper right corner
    plt.text(0.95, 0.05,
             f'Length of fft2sort: {len(fft2sort)}\nThreshold: {thresh:.0f}\n#Transmitted: {np.count_nonzero(ind)}',
             horizontalalignment='right',
             verticalalignment='bottom',
             transform=plt.gca().transAxes,
             bbox=dict(facecolor='white', alpha=0.6))
    plt.show()

出力結果は以下の通り。

このようにFFTを用いれば特定の周波数帯域を分離することが容易であるので、同じ要領でノイズ除去や信号のフィルタリングなどにも利用される。
では、ガウシアンノイズを加えた画像に対するノイズキャンセリングを実際に行ってみる。高周波モードにおいてノイズが特に顕著であることが観察されるため、低周波数を含む所定の半径外のフーリエ係数をゼロにする。

まずは読み込んだ画像に適当なガウシアンノイズを加える。

# Import libraries
from matplotlib.image import imread
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams.update({'font.size': 12})

# Read data & add some noise
fig = imread('../data/banta.jpg')
fig = np.mean(fig, -1)
noise = fig + 200*np.random.randn(*fig.shape).astype('uint8')
plt.imshow(noise)
plt.show()

次にこのノイズデータに対して2次元フーリエ変換を実施

# 2-dim FFT
fft2noise = np.fft.fft2(noise)
shiftnoise = np.fft.fftshift(fft2noise)
F = np.log(np.abs(shiftnoise) + 1)
plt.imshow(F)
plt.show()

ここで、真ん中の高周波数モードを残すように適当な半径を決めて、その円外にあるモードを落とす処理を行う。

# Create mesh grid
nx, ny = fft2noise.shape
X, Y = np.meshgrid(
    np.arange(-ny/2 + 1, ny/2 + 1),
    np.arange(-nx/2 + 1, nx/2 + 1)
)

# Noise filtering
R2 = np.power(X, 2) + np.power(Y, 2)
ind = R2 < 150**2
shiftnoisefilt = shiftnoise * ind
Ffilt = np.log(np.abs(shiftnoisefilt)+1)

# Reconstruct the original image
filt = np.fft.ifftshift(shiftnoisefilt)
filt = np.fft.ifft2(filt).real

plt.imshow(filt)
plt.show()
plt.imshow(Ffilt)
plt.show()

画像データに対する2次元ウェーブレット変換

FFTと同じく、離散ウェーブレット変換（DWT）も画像処理/圧縮によく用いられる手法である。画像処理において、ウェーブレット変換は「時間」に相当する軸として、画像の水平軸（x軸）および垂直軸（y軸）を使用する。したがって、時間に相当する概念は、画像の空間的な位置（ピクセルの位置）に対応することとなる。

以下では画像データに対する2階層のウェーブレット変換を行ってみる。母ウェーブレットとしてHaarウェーブレット$\psi(t)$を与える。

\begin{align}
\psi(t)
=
\begin{cases}
1\quad&\text{for }0\leq t<1/2, \\
-1\quad&\text{for }1/2\leq t<1, \\
0\quad&\text{otherwise.}
\end{cases}
\end{align}

# Import libraries
from matplotlib.image import imread
import matplotlib.pyplot as plt
import numpy as np
import pywt
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams.update({'font.size': 12})

# Read data
fig = imread('../data/banta.jpg')
fig = np.mean(fig, -1)

# Wavelet decomposition (2 level)
n = 2
coeffs = pywt.wavedec2(fig, wavelet='db1', level=n)

# normalize each coefficient array
coeffs[0] /= np.abs(coeffs[0]).max()
for detail_level in range(n):
    coeffs[detail_level + 1] = [d/np.abs(d).max()
                                for d in coeffs[detail_level + 1]]
arr, coeff_slices = pywt.coeffs_to_array(coeffs)

plt.imshow(arr, cmap='gray', vmin=-0.1, vmax=0.1)
plt.show()

かなり見えにくいが、左上のブロックは低周波成分を示し、元の画像の大まかな形状や低周波成分を表す。その他のブロック($3\times2$階層$=6$ブロック)は高周波成分であり、エッジやテクスチャなどを表す。各ブロックは水平方向、垂直方向、および対角方向の詳細係数に対応する。

最後にウェーブレット変換の階層表現を用いた画像圧縮をデモする。

# Read data
fig = imread('../data/banta.jpg')
fig = np.mean(fig, -1)

# Wavelet decomposition (4 level)
n = 8
w = 'db1'
coeffs = pywt.wavedec2(fig, wavelet=w, level=n)

# Wavelet Compression
coeff_arr, coeff_slices = pywt.coeffs_to_array(coeffs)
Csort = np.sort(np.abs(coeff_arr.reshape(-1)))

for keep in [0.05, 0.01, 0.002]:
    # Find small indices
    thresh = Csort[int(np.floor((1-keep) * len(Csort)))]
    ind = np.abs(coeff_arr) > thresh
    
    # Threshold small indices
    Cfilt = coeff_arr * ind

    # Compressed image
    coeffs_filt = pywt.array_to_coeffs(
        Cfilt, coeff_slices, output_format='wavedec2')
    
    # Show image
    Arecon = pywt.waverec2(coeffs_filt, wavelet=w)
    plt.imshow(Arecon.astype('uint8'))
    plt.title('Compressed image: keep = ' + str(keep * 100) + '% of DWT')

    # Plot reconstruction
    plt.text(0.95, 0.05,
            f'Length of Csort: {len(Csort)}\nThreshold: {thresh:.0f}\n#Transmitted: {np.count_nonzero(ind)}',
            horizontalalignment='right',
            verticalalignment='bottom',
            transform=plt.gca().transAxes,
            bbox=dict(facecolor='white', alpha=0.6))
    plt.show()

DWT係数の0.2%しか残さないような攻めた圧縮でも、画像の粗い特徴は保持される。（勿論、分解する階層の数にも依存）例えば、帯域幅が制限される状況でデータを送信する際、DWT情報の多くが切り捨てられてもデータの最も重要な特徴を転送することが出来る。

おわりに

画像圧縮/処理する際にフーリエ変換とウェーブレット変換、どういうケースでどちらを用いるべきかがまだ把握しきれていないが、

エッジ検出やテクスチャに特化したい→ウェーブレット変換
特定の周波数成分を詳細に解析したい→フーリエ変換
だという理解だ。このあたりは必要に応じて調査したいと思う。