前回、ガボールフィルタを実装して「特定の向き・特定の細かさの模様に反応するフィルタ」を実際に確認した。あのフィルタは、ガウス関数とcos波を組み合わせたものだった。触っているうちに、素朴な疑問が浮かんだ。cos波を局所的に当てるのがガボールフィルタなら、画像全体をそのままcos波(とsin波)に分解したら何が見えるのか。 それがフーリエ変換だ。

というわけで、今回は画像の2次元フーリエ変換を実際にPythonで実装し、スペクトルを覗いてみる。

フーリエ変換は何をしているのか

2次元離散フーリエ変換(DFT)は、画像 $f(x, y)$ を、いろいろな向き・周波数のsin波・cos波の足し合わせとして書き直す変換だ。

$$F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) \, e^{-2\pi i \left( \frac{ux}{M} + \frac{vy}{N} \right)}$$

式だけ見ると身構えるが、やっていることはガボールフィルタと大差ない。$(u, v)$ が「向きと周波数」を指定するパラメータで、$F(u, v)$ はその向き・周波数の波が元画像にどれだけ含まれているかを表す複素数だ。ガボールフィルタが「特定の1つの向き・周波数だけ」を局所的に取り出していたのに対し、フーリエ変換は「画像全体を、あらゆる向き・周波数の波の足し合わせ」として一気に分解する。

$F(u, v)$ は複素数なので、絶対値(振幅)と偏角(位相)の2つの情報を持つ。振幅がその波の強さ、位相がその波の位置(ずれ)を表す。この振幅と位相、後半で実際に入れ替えてみる。

人工的な縞模様で試す

まず、ガボール記事でも使った人工的な縞模様で挙動を確認する。numpy.fft.fft2 で2次元DFTを計算し、fftshift で低周波成分を中心に寄せ、log(1 + |F|) でスケールを圧縮して可視化する。

import numpy as np

def make_stripe(size, angle_deg, wavelength):
    x, y = np.meshgrid(np.arange(size), np.arange(size))
    theta = np.deg2rad(angle_deg)
    xr = x * np.cos(theta) + y * np.sin(theta)
    return np.sin(2 * np.pi * xr / wavelength)

def spectrum(img):
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    return np.log(1 + np.abs(fshift))

「垂直・細かい」「垂直・粗い」「斜め45度」の3パターンを用意し、それぞれのスペクトルを並べた。

垂直の縞2つは、想定通り中心から左右に伸びる軸上に点が現れ、縞が細かいほど点が中心から離れる(=高周波)。ここまでは予想通りだった。

意外だったのは斜め45度の縞だ。点ではなく、十字型に薄く伸びる模様になった。最初はバグかと思ったが、原因は境界条件にある。DFTは「画像が上下左右に無限に繰り返しタイル状に並んでいる」という前提で計算される。垂直・水平の縞は画像の一辺と綺麗に噛み合うため境界で不連続が起きにくいが、45度の縞は画像の端で模様が食い違い、そこに不連続な境目ができる。この不連続がスペクトル上でノイズのように広い周波数帯に漏れ出し、十字型の模様(スペクトル漏れ、spectral leakage)として見えていた。理論の教科書には出てこない、実際に手を動かして初めて気づいた挙動だった。

実写真のスペクトルを見る

人工的な縞ではなく、実写真(定番のテスト画像camera)でも同じことをやってみる。

from skimage import data

camera = data.camera().astype(np.float64)
mag_camera = spectrum(camera)

こちらにも、さっきの人工縞と同じ十字型の筋がくっきり出ている。理由は同じで、写真の四辺(画像の境界)がそのまま不連続点になり、水平・垂直方向にスペクトル漏れが発生するためだ。実は写真のフーリエスペクトルを見るとほぼ必ずこの十字が現れる。「画像の端」という、内容とは無関係な要因がスペクトルに影響しているという事実は、最初にガボールフィルタを触ったときには意識していなかった点だ。

低周波・高周波だけを取り出す

スペクトル上で円形のマスクをかけて一部の周波数だけを残し、逆変換で画像に戻してみる。中心(低周波)だけ残せばぼかし、逆に中心以外(高周波)だけ残せば輪郭が浮き出るはずだ。

f = np.fft.fft2(camera)
fshift = np.fft.fftshift(f)

rows, cols = camera.shape
crow, ccol = rows // 2, cols // 2
Y, X = np.ogrid[:rows, :cols]
dist = np.sqrt((X - ccol) ** 2 + (Y - crow) ** 2)

radius = 30
low_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * (dist <= radius))))
high_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * (dist > radius))))

想定通り、低周波だけ残すとぼかし画像に、高周波だけ残すと輪郭線画像になった。ただし低周波側の画像をよく見ると、うっすらと波打つような模様(リンギング)が乗っている。これは円形マスクで周波数を「スパッと」切り落としたことによる副作用で、理想的な急峻なフィルタ(いわゆるブリックウォールフィルタ)をかけると、逆変換した画像に波紋状のアーティファクトが出るという、信号処理の教科書でよく説明される現象を実際に自分の画像で確認できた形になる。実務でぼかし処理をするときにガウシアンブラーのような滑らかなフィルタが好まれる理由が、手を動かすと腑に落ちた。

振幅と位相、画像の「形」を決めているのはどっちか

ここが今回一番面白かった実験。フーリエ変換した複素数 $F(u,v)$ を振幅 $|F|$ と位相 $\angle F$ に分解し、2枚の別々の写真の振幅と位相を入れ替えて逆変換したら何が起きるか試した。

photo_a = data.camera().astype(np.float64)
photo_b = transform.resize(
    color.rgb2gray(data.astronaut()), photo_a.shape, anti_aliasing=True
) * 255

この2枚のフーリエ変換から、振幅と位相を入れ替えて合成する。

Fa, Fb = np.fft.fft2(photo_a), np.fft.fft2(photo_b)
amp_a, phase_a = np.abs(Fa), np.angle(Fa)
amp_b, phase_b = np.abs(Fb), np.angle(Fb)

mix1 = amp_a * np.exp(1j * phase_b)  # Aの振幅 + Bの位相
mix2 = amp_b * np.exp(1j * phase_a)  # Bの振幅 + Aの位相

結果は予想以上にはっきりしていた。「Aの振幅+Bの位相」は明るさの分布こそAに近いのに、シルエットは完全にB(宇宙飛行士)。逆に「Bの振幅+Aの位相」は、Aの構図(カメラマン)がそのまま浮かび上がる。つまり、画像の「形・輪郭・構図」を決めているのはほぼ位相の方で、振幅はコントラストや質感程度の情報しか持っていない。

これは事前に文献で知識としては見たことがあったが、自分の手元の画像で実際にやってみると想像以上に振幅側の情報量のなさに驚いた。ガボールフィルタの記事で「画像から特徴を抽出するとは何か」と考えていたが、少なくとも古典的なフーリエの枠組みでは「形」の情報のほとんどは位相に閉じ込められているらしい。

まとめ

ガボールフィルタが「特定の向き・周波数を局所的に覗き見るフィルタ」だとすれば、フーリエ変換は「画像全体を、あらゆる向き・周波数の波の重ね合わせとして一気に展開する」操作だった。実際に手を動かしてみて印象に残ったのは、教科書的な低周波=ぼかし・高周波=輪郭という話よりも、画像の境界が生む十字型のスペクトル漏れと、位相が画像の「形」のほとんどを担っているという2点だった。どちらも、式を眺めているだけでは実感として残らなかったと思う。

• 実践記
• 機械学習
• 画像処理

関連記事

• 今さらガボールフィルタを触ってみた。画像の「向き」を見るとはどういうことか
• 二重降下を実験で再現した。補間閾値でテスト誤差が76.8倍に跳ね上がり、そこから64分の1まで下がったが、シンプルな次数5モデルには最後まで届かなかった
• 「破局的忘却」を実際に起こしてみたら、Task A精度は0.9834から0.0000まで消えた。10%のリハーサルでも一度は9.9%まで沈んでから戻ってきた