Juliaで二次元フーリエ変換と画像処理(ローパス、ハイパスフィルタ)

はじめに
２次元フーリエ変換
ローパスハイパスフィルタの実装
1. 関連

はじめに

JuliaのFFTについて、公式documentや日本語サイトがあまりなかったので、2次元フーリエ変換をする方法を備忘録として残しておこうと思う。

2Dフーリエ変換は画像処理分野や様々な分野でよく使われている。中でもレーザー光学の分野ではフーリエ光学と呼ばれる分野もある。特によく使われるところではFraunhofer回折積分やレンズのフーリエ変換作用などがある。ここではそれらの計算に使うためのツールとして基本的な２次元フーリエ変換をJuliaで実装しようと思う。

２次元フーリエ変換

using Pkg
Pkg.add("FFTW")

FFTのためにFFTWというパッケージを導入する。基本的なドキュメントはAbstractFFTs(https://juliamath.github.io/AbstractFFTs.jl/stable/api/)というサイトにある。基本的なFFTでよく使う機能はフーリエ変換である「fft」フーリエ変換後の座標を調整する「fftshift」逆フーリエ変換「ifft」の３つである。これらの諸機能については以下を参考にして下さい。

fft(A [, dims])
基本的な使い方はAにFFTしたいデータを入れればいいです、ここでAは複数次元でも勝手にFFTしてくれます。dimsには特定の軸だけFFTしたいときに使います。

以下、ドキュメントの和訳

配列 A の多次元高速フーリエ変換（FFT）を実行します。オプションの dims 引数は、変換を行う次元の部分集合（整数、範囲、タプル、または配列など）を指定します。変換対象の次元のサイズが小さな素因数の積であると、最も効率的に処理されます（Base.nextprod を参照）。さらに効率を高めたい場合は、plan_fft() も参照してください。

fftshift(x, [dim])
基本的な使い方はFFTしたあとに使います。XにはFFT後のデータを入れます。なぜ使うかは上記過去記事を読んでください。これを使った後に逆フーリエ変換するにはifftshiftをする必要があります。

以下、ドキュメントの和訳

周期信号 x に対して、指定された次元 dim に沿って円形シフト（circular shift）を行います。これにより、インデックス1に中心がある信号を、サイズ N に対して N÷2+1 を中心とする形にシフトします（ここで N はその次元のサイズです）。

この操作は ifftshift を使って元に戻すことができます。N が偶数の場合、これは配列の前半と後半を入れ替えることと等しくなり、そのため fftshift と ifftshift は同じ動作になります。

dim を指定しなかった場合、すべての次元に対してシフトが行われます。

fftshift は新しい配列に結果を出力します。あらかじめ用意した配列に出力したい場合は、fftshift! を使用してください。

ifft(A [, dims])
多次元の逆高速フーリエ変換（逆FFT）を実行します。基本的な使い方はfftと同じです。

ガウス関数はフーリエ変換したあともガウス関数であることが知られています。そのため、今回ガウス関数を使ってフーリエ変換の様子を確認してみます。

using FFTW
using PyPlot
plt = PyPlot

# 2D Gaussian function
gaussian2d(x, y, sigma) = exp(-(x^2 + y^2) / (2 * sigma^2)) / (2 * π * sigma^2)

# Create x, y meshgrid
max_size = 15
n = 100
x = range(-max_size, max_size, length=n)
y = range(-max_size, max_size, length=n)

# Generate 2D grid
X = repeat(reshape(x, 1, :), n, 1)
Y = repeat(reshape(y, :, 1), 1, n)

# Evaluate 2D Gaussian on grid
sigma = 1.0
gauss_mesh = gaussian2d.(X, Y, sigma)

# Plot original 2D Gaussian
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(gauss_mesh, extent=(-max_size, max_size, -max_size, max_size), cmap="hot")
plt.title("2D Gaussian Function")   
plt.xlabel("x")
plt.ylabel("y")
plt.colorbar()


# FFT and plot spectrum
FFT_gauss = fftshift((fft(gauss_mesh)))
FFT_gauss = abs.(FFT_gauss).^2  # Compute the power spectrum

plt.subplot(1, 2, 2)
plt.imshow(FFT_gauss, extent=(-max_size, max_size, -max_size, max_size), cmap="hot")
plt.title(L"2D Gaussian Spectrum (|FFT|^2)")
plt.xlabel("kx")
plt.ylabel("ky")
plt.colorbar()
display(plt.gcf())

無事、FFTされていることが確認できました。

ローパスハイパスフィルタの実装

ついでにいぜんやっていたローパスハイパスフィルタも実装しました。以下のコードをJuliaに直したものになります。

using FFTW
using Images
using PyPlot
using TestImages
plt = PyPlot

function lowpassfilter(shift_fft_img, pas)
    h, w = size(shift_fft_img)
    cy, cx = div(h, 2), div(w, 2)
    rh, rw = round(Int, pas * cy), round(Int, pas * cx)
    
    fdst = zeros(ComplexF64, h, w)
    fdst[cy - rh + 1 : cy + rh, cx - rw + 1 : cx + rw] .= shift_fft_img[cy - rh + 1 : cy + rh, cx - rw + 1 : cx + rw]
    
    fds = ifftshift(fdst)
    ds = real(ifft(fds))
    ds_img = Gray.(clamp.(ds, 0, 1))
    
    return fdst, ds_img
end

function highpassfilter(shift_fft_img, pas)
    h, w = size(shift_fft_img)
    cy, cx = div(h, 2), div(w, 2)
    rh, rw = round(Int, pas * cy), round(Int, pas * cx)
    
    filtered = copy(shift_fft_img)
    filtered[cy - rh + 1 : cy + rh, cx - rw + 1 : cx + rw] .= 0
    
    fds = ifftshift(filtered)
    ds = real(ifft(fds))
    ds_img = Gray.(clamp.(ds, 0, 1))
    
    return filtered, ds_img
end

function plot_4(img,am,fdst,ds)
    fig = plt.figure(figsize=(18, 9))
    ax1 = fig.add_subplot(1,4,1)
    ax2 = fig.add_subplot(1,4,2)
    ax3 = fig.add_subplot(1,4,3)
    ax4 = fig.add_subplot(1,4,4)
    ax1.imshow(Float64.(img), cmap = "gray")
    ax2.imshow(am, cmap = "gray")
    ax3.imshow(20 .* log.(abs.(fdst) .+ 1e-5), cmap = "gray")
    ax4.imshow(Float64.(ds), cmap = "gray")
    display(plt.gcf())
end

img = testimage("cameraman")
imgf = channelview(float.(img))  # ← Gray{N0f8} → Float32 → Array{Float32,2}
f = fftshift(fft(imgf))
amplitude = log.(1 .+ abs.(f))  # 振幅スペクトルの可視化用

# 低域フィルタ適用
filtered_fft, filtered_img = lowpassfilter(f, 0.2)

# 結果のプロット（関数名修正も必要）
plot_4(img, amplitude, filtered_fft, filtered_img)