フーリエフィルタリングによるノイズ除去 - 周波数領域での画像処理
フーリエ変換と画像の周波数表現
フーリエ変換 (Fourier Transform) は、画像を空間領域から周波数領域に変換する数学的操作です。空間領域では各ピクセルの輝度値として表現される画像が、周波数領域では様々な周波数の正弦波の重ね合わせとして表現されます。この変換により、画像の周波数成分を直接操作してノイズ除去やエッジ強調が可能になります。
2D 離散フーリエ変換 (DFT):
M×N の画像 f(x,y) の 2D DFT は F(u,v) = ΣΣ f(x,y) × exp(-j2π(ux/M + vy/N)) で定義されます。F(u,v) は周波数 (u,v) における複素数値で、振幅 |F(u,v)| と位相 arg(F(u,v)) を持ちます。振幅スペクトルは各周波数成分の強さを、位相スペクトルは各成分の位置情報を表します。
周波数の意味:
画像における低周波成分は、緩やかな輝度変化 (大きな構造、全体的な明るさ) に対応します。高周波成分は、急激な輝度変化 (エッジ、テクスチャ、ノイズ) に対応します。DC 成分 F(0,0) は画像の平均輝度を表します。スペクトル画像では中心が低周波、周辺が高周波を示します。
FFT による高速計算:
DFT の直接計算は O(N⁴) ですが、FFT (Fast Fourier Transform) アルゴリズムにより O(N² log N) に削減されます。NumPy の np.fft.fft2() や OpenCV の cv2.dft() が FFT を実装しています。1,024 × 1,024画像の FFT は現代の CPU で数ミリ秒で完了します。
ローパスフィルタ - ノイズ除去の基本
ローパスフィルタ (Low-Pass Filter) は、低周波成分を通過させ高周波成分を減衰させるフィルタです。ノイズは一般に高周波成分として現れるため、ローパスフィルタによりノイズを除去できます。ただし、エッジも高周波成分を含むため、ノイズ除去とエッジ保持のトレードオフが存在します。
理想ローパスフィルタ:
カットオフ周波数 D0 以下の成分をすべて通過させ、D0 以上をすべて遮断する理想的なフィルタです。H(u,v) = 1 (D(u,v) <= D0), H(u,v) = 0 (D(u,v) > D0)。ただし、急峻な遮断によりリンギング (Gibbs 現象) が発生し、画像にリング状のアーティファクトが生じます。実用的にはほとんど使用されません。
Butterworth ローパスフィルタ:
H(u,v) = 1 / (1 + (D(u,v)/D0)^(2n)) で定義される滑らかな遮断特性を持つフィルタです。次数 n が大きいほど遮断が急峻になり、n=1 で最も滑らかな遮断、n→∞ で理想フィルタに近づきます。n=2 が実用的なバランスとして広く使用されています。リンギングが少なく、自然なぼかし効果が得られます。
ガウシアンローパスフィルタ:
H(u,v) = exp(-D(u,v)² / (2D0²)) で定義されるガウス関数型のフィルタです。リンギングが全く発生しない最も滑らかな遮断特性を持ちます。空間領域でのガウシアンぼかし (cv2.GaussianBlur) と等価で、σ = M/(2πD0) の関係があります。ノイズ除去の品質は Butterworth と同程度ですが、リンギングフリーな点が利点です。
カットオフ周波数の選択:
カットオフ周波数 D0 の選択がフィルタリング結果を大きく左右します。D0 が小さいほどノイズ除去効果は高いですが、画像のぼけも大きくなります。一般的には、パワースペクトルの 90-95% を含む周波数を D0 とする方法や、ノイズのスペクトル特性から D0 を決定する方法が使用されます。
ハイパスフィルタとバンドパスフィルタ
ハイパスフィルタは高周波成分を通過させ低周波を減衰させるフィルタで、エッジ検出やシャープニングに使用されます。バンドパスフィルタは特定の周波数帯域のみを通過させ、特定のテクスチャや構造の抽出に使用されます。
ハイパスフィルタ:
ハイパスフィルタは H_HP(u,v) = 1 - H_LP(u,v) で定義されます (H_LP はローパスフィルタ)。画像のエッジやテクスチャを強調し、全体的な明るさ (DC 成分) を除去します。結果画像は暗くなるため、元画像に加算する「アンシャープマスク」として使用するのが一般的です。enhanced = original + k × highpass (k はシャープニング強度) の形式です。
バンドパスフィルタ:
H_BP(u,v) = H_LP(D_H) - H_LP(D_L) で定義され、D_L から D_H の周波数帯域のみを通過させます。特定のスケールのテクスチャを抽出したり、特定の周波数帯のノイズを分離したりするのに使用されます。医療画像での特定構造の強調や、リモートセンシングでの地物分類に応用されています。
バンドリジェクトフィルタ:
H_BR(u,v) = 1 - H_BP(u,v) で定義され、特定の周波数帯域を除去します。周期的なノイズ (電源ハムによる縞模様など) の除去に特に有効です。ノイズの周波数が既知の場合、その周波数帯域のみを除去することで、画像の他の情報を保持したままノイズを除去できます。
ホモモルフィックフィルタリング:
画像を照明成分 (低周波) と反射成分 (高周波) に分離し、それぞれを独立に処理する手法です。対数変換後にフーリエ変換を適用し、低周波を抑制・高周波を強調するフィルタを適用します。不均一な照明条件下で撮影された画像のコントラスト改善に効果的で、顔認識の前処理などに使用されます。
ノッチフィルタ - 周期ノイズの除去
ノッチフィルタ (Notch Filter) は、周波数スペクトル上の特定の点 (ノッチ) を除去するフィルタです。周期的なノイズ (モアレ、電磁干渉による縞模様、スキャン時のバンディングなど) の除去に特化しています。
周期ノイズの特徴:
周期的なノイズは、周波数スペクトル上で明確なピーク (スパイク) として現れます。例えば、水平方向の縞模様は垂直軸上の特定周波数にピークを生じ、斜め方向の縞模様は対応する角度の位置にピークを生じます。これらのピークを特定し、除去することで周期ノイズを効果的に除去できます。
ノッチフィルタの設計:
ノッチフィルタは、スペクトル上の特定の点 (u0, v0) とその対称点 (-u0, -v0) を中心とした小さな領域を除去します。フーリエ変換の対称性により、実画像のスペクトルは原点対称なので、ノッチは常にペアで設定します。ガウシアン型のノッチ H(u,v) = 1 - exp(-D1²/(2σ²)) × exp(-D2²/(2σ²)) が滑らかな除去を実現します (D1, D2 はノッチ中心からの距離)。
ノッチ位置の特定:
周期ノイズのノッチ位置は、パワースペクトル (|F(u,v)|²) を可視化して手動で特定するか、自動検出アルゴリズムで特定します。自動検出では、DC 成分を除外した後にスペクトルのピークを検出します。ピーク検出には局所最大値の探索や、閾値ベースの検出が使用されます。
実践例 - スキャン画像のモアレ除去:
印刷物をスキャンした画像には、印刷の網点パターンとスキャナの解像度の干渉によるモアレが発生します。スペクトルを観察すると、網点の周波数に対応する規則的なピークが確認できます。これらのピークにノッチフィルタを適用することで、モアレを除去しつつ画像の内容を保持できます。ノッチの幅 (σ) は、ピークの広がりに合わせて調整します。
ウィーナーフィルタと逆フィルタリング
ウィーナーフィルタ (Wiener Filter) は、ノイズと信号の統計的特性を考慮した最適なフィルタリング手法です。単純なローパスフィルタと異なり、信号対ノイズ比 (SNR) に基づいて周波数ごとに最適な減衰量を決定します。
逆フィルタリングの問題:
劣化画像 G(u,v) = H(u,v)F(u,v) + N(u,v) から元画像 F を復元する最も単純な方法は、逆フィルタ F_hat = G/H です。しかし、H(u,v) が 0 に近い周波数ではノイズ N/H が増幅され、結果が破綻します。これが逆フィルタリングの根本的な問題です。
ウィーナーフィルタの定式化:
ウィーナーフィルタは、復元画像と元画像の平均二乗誤差を最小化する最適フィルタです。W(u,v) = H*(u,v) / (|H(u,v)|² + S_n(u,v)/S_f(u,v)) で定義されます。H* は H の複素共役、S_n はノイズのパワースペクトル、S_f は信号のパワースペクトルです。SNR が高い周波数では逆フィルタに近づき、SNR が低い周波数では減衰が強くなります。
正則化パラメータ:
実際にはノイズと信号のパワースペクトルが未知のことが多いため、S_n/S_f を定数 K で近似します。W(u,v) = H*(u,v) / (|H(u,v)|² + K) の形式で、K はノイズレベルに応じて調整します。K が大きいほどノイズ抑制が強く (ぼけが大きく)、K が小さいほど復元が鋭く (ノイズが残る) なります。
応用例 - モーションブラーの除去:
カメラの手ぶれによるモーションブラーは、周波数領域では特定方向の sinc 関数として表現されます。ブラーの方向と長さを推定し、対応する劣化関数 H を構築してウィーナーフィルタを適用することで、ぶれた画像を復元できます。ブラーパラメータの推定には、スペクトルのゼロ交差パターンの分析やケプストラム解析が使用されます。
実装ガイド - Python でのフーリエフィルタリング
NumPy と OpenCV を使用したフーリエフィルタリングの実装を解説します。FFT の計算からフィルタ設計、逆変換まで、実践的なパイプラインを構築します。
FFT と逆 FFT の基本:
NumPy での FFT は F = np.fft.fft2(image) で計算し、F_shift = np.fft.fftshift(F) で低周波を中心に移動します。フィルタ適用後は F_ishift = np.fft.ifftshift(F_filtered) で元に戻し、result = np.abs(np.fft.ifft2(F_ishift)) で空間領域に逆変換します。
ガウシアンローパスフィルタの実装:
フィルタマスクの生成は、画像中心からの距離マップを計算し、ガウス関数を適用します。rows, cols = image.shape; crow, ccol = rows//2, cols//2; D = np.sqrt((np.arange(rows)[:,None]-crow)**2 + (np.arange(cols)[None,:]-ccol)**2); H = np.exp(-D**2/(2*D0**2)) で生成し、F_filtered = F_shift * H で適用します。
ノッチフィルタの実装:
スペクトルを可視化して周期ノイズのピーク位置を特定し、その位置にガウシアンノッチを配置します。magnitude = 20*np.log(np.abs(F_shift)+1) でスペクトルを可視化し、ピーク座標 (u0, v0) を特定します。ノッチフィルタは H = 1 - np.exp(-D1**2/(2*sigma**2)) * np.exp(-D2**2/(2*sigma**2)) で生成します (D1, D2 はノッチ中心からの距離)。
パフォーマンスの考慮:
大きな画像では FFT の計算時間が問題になることがあります。画像サイズを 2 のべき乗にパディングすると FFT が最も効率的に動作します。optimal_size = cv2.getOptimalDFTSize(n) で最適サイズを取得できます。また、実数画像の FFT は対称性を利用した np.fft.rfft2() で計算量を半減できます。GPU 加速には CuPy の cupyx.scipy.fft.fft2() が利用可能で、CPU の 10-50 倍の高速化が期待できます。