はじめに

数式を暗記していた自分が、もう少し概念的に理解できるようにがんばった記録を（自己満で）共有させていただきます。間違いやミスにつきましては、修正していきますので、お気づきの点はコメントください。

画像は波でできている

画像はさまざまな波の画像からできています（言い換えると、波の画像から再構成できます）。
波の画像は、基底関数 ( basis function ) と呼ばれています。

波の計算

山は山側、谷は谷側、山と谷は差し引き、波の計算は簡単です。

フーリエ変換（Fourier Transform）

先ほどの例のように、いろいろな波の足し合わせで画像のピクセル値（輝度）が表現できます。

波を分解したり、元に戻したりできる？→Yes。フーリエ変換でできます。

一次元フーリエ変換

一つの波がどのような波の組み合わせでできているかを計算するためには、周波数に分解するという作業を行います。

連続した世界

F(\omega)=\int_{-\infty}^{\infty}f(x)e^{-j2\pi\omega x}dx, \omega\in[-\infty,\infty]

離散化した世界

F(u)=\sum_{n=0}^{N-1}f(x)e^{-j\frac{2\pi ux}{N}}, [0 \leq u \leq N-1]

インテグラルは連続変数の範囲での合計を表し、シグマは離散変数の合計を表します。離散化した世界でいうf(x)は一次元の数値の配列です。

複素数としての扱い

フーリエ変換の結果は複素数になります。

e^{-j\frac{2\pi ux}{N}}=cos\frac{2\pi ux}{N}-j\frac{2\pi ux}{N}

ここで、jは虚数（二乗すると-1になる数の概念）です。
cos成分は実部（虚数を含まない）、sin成分は虚部と捉えることができます。これはオイラーの公式に従って導かれます。

e^{jx}=cosx+jsinx

e^{-jx}=cosx-jsinx

フーリエ変換の概念イメージ

実部はそのままの値でコンピュータ上で実数として扱えます。
フーリエ変換した実部を正弦波として表示した例を示します。
（波の成分ごとに分かれていることがわかります。この画像は256×256サイズのため、本当は256の波がありますが、ここではナイキスト周波数から16番目までのみを表示しています。）

from numpy.fft import fft

# 256*256 胸部画像の行データを利用する
x = c_row
#フーリエ変換を実施
freq = fft(x)
#結果を絶対値で取得（結果が複素数で返ってくるため）
freq_abs = np.abs(freq) # fft result
#グラフにして、左右でシンメトリーになることを確認。
#ナイキスト周波数の位置（N/2）以降が虚像（pretense）
#特殊な操作であるが，虚像を分かりやすくするために配列の前半と後半を入れ替える
swap = np.zeros(256)
swap[:128] = freq_abs[128:256]
swap[128] = freq_abs[0]
swap[129:] = freq_abs[1:128]
#グラフ用のデータへ
data = {'freq':swap,
        'sampling size':np.linspace(0,256,256)
        }
df = pd.DataFrame(data=data)
#シンメトリーをグラフへ
fig = px.line(df,x='sampling size',y='freq',title='fft result with the pretense at row128')
fig.show()
#実際の計算にはナイキスト周波数の位置（N/2）以降の虚像部分は不要。
#ナイキスト周波数の位置（N/2）までをグラフへ
#この例ではN=256のため、その位置は128。
data = {'freq':freq_abs[0:128],
        'sampling size':np.linspace(0,128,128)
        }
df = pd.DataFrame(data=data)
fig = px.line(df,x='sampling size',y='freq',title='fft result without the pretense at row128')
fig.show()

# 3D 可視化
#再び胸部の128行目のピクセル配列を利用
x = c_row
#ナイキスト周波数までを取得
freq = fft(x)[0:128]
freq_max = np.max(np.abs(freq))

N = len(freq) # データ数
n = np.arange(N) # 1きざみの配列をN個

#各周波数をsin波の信号に変換
wave_arr = []
name_arr = []
location_arr = []
for i in range(N):
  norm = np.abs(freq[i])/freq_max # 正規化
  freq_sin_wave = np.sin(np.abs(freq[i])*2*np.pi*(n/N))*norm
  name = ['freq'+str(i)] * N
  if i < 16 : # 16番目までの低周波成分を可視化
    wave_arr.extend(freq_sin_wave)
    name_arr.extend(name)
    location_arr.extend(n)
#グラフデータへ
data = {'amplitude':wave_arr,
        'sampling location':location_arr,
        'freq':name_arr,
        }
df = pd.DataFrame(data=data)
fig = px.line_3d(df, x="sampling location", y="freq", z="amplitude", color='freq')
fig.show()

二次元離散フーリエ変換

一次元フーリエ変換を行・列にわけて行うことで、二次元上にあるピクセルの輝度を周波数空間に変換します。
2D-DFT, 2 Dimensional-Discrete Fourier Transform

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

これは次のようにも書けます。
このように変換すると、行と列で計算していることがわかりやすいでしょう。

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

基底関数

expの部分は、基底関数として表現することができます。

入力画像のサイズ、周波数空間の座標によって、基底関数は決定されます。それぞれ、入力画像と同じ16×16のサイズになります（M×N）。
縞模様が細かいところは高周波成分、大きいところは低周波成分を表現します。

実部

虚部

# 基底関数をつくる
import numpy as np
import matplotlib.pyplot as plt

'''
img_size : square size of image f(x,y)
u,v : spatial space indice
'''
def basis_function(img_size=256, u=0, v=0):
  N = img_size
  x = np.linspace(0, N-1, N)
  y = np.linspace(0, N-1, N)
  x_, y_ = np.meshgrid(x, y)
  # print(x_.shape)
  # print(y_.shape)
  # print(x_)
  # print(y_)
  bf = np.exp(-1j*2*np.pi*(u*x_/N+v*y_/N))
  if u == 0 and v == 0:
    bf = np.round(bf)
  real = np.real(bf)#実部を取得
  imag = np.imag(bf)#虚部を取得
  return real, imag

# r,i = basis_function(img_size=3, u=1, v=1)
# plt.imshow(r, cmap='gray')
# plt.show()
# '''
# 虚部はjを外した値で描画される
# a = np.array([1+2j, 3+4j, 5+6j])
# a.imag
# array([2.,  4.,  6.])
# '''
# plt.imshow(i, cmap='gray')
# plt.show()

# basis functionの確認
size = 16
bf_arr_real = np.zeros((size*size,size,size))
bf_arr_imag = np.zeros((size*size,size,size))
ind = 0
for col in range(size):
  for row in range(size):
    re,imag = basis_function(img_size=size, u=row, v=col)
    bf_arr_real[ind] = re
    bf_arr_imag[ind] = imag
    ind += 1

# real part
_, axs = plt.subplots(size, size, figsize=(7, 7))
axs = axs.flatten()
for img, ax in zip(bf_arr_real, axs):
  # ax.set_title(str(h)+","+str(w))
  ax.set_axis_off()
  ax.imshow(img,cmap='gray')
#plt.tight_layout()
plt.show()

# imaginary part
_, axs = plt.subplots(size, size, figsize=(7, 7))
axs = axs.flatten()
for img, ax in zip(bf_arr_imag, axs):
  # ax.set_title(str(h)+","+str(w))
  ax.set_axis_off()
  ax.imshow(img,cmap='gray')
#plt.tight_layout()
plt.show()

#2次元離散フーリエ逆変換

周波数空間から元の画像を復元できます。
2D-iDFT, 2 Dimensional-Inverse Discrete Fourier Transform

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

2D-DFTと2D-iDFT

以上が2D-DFTと2D-iDFTです。

import cv2
import numpy as np
from numpy.fft import *
import matplotlib.pyplot as plt

H = 16
W = 16

chest = cv2.imread('path to ChestXray256.png',cv2.IMREAD_GRAYSCALE).astype(float)
# サイズが大きいので16×16へダウンサンプル
chest = cv2.resize(chest,(W,H),cv2.INTER_LANCZOS4)
# リサイズした画像を確認
plt.subplot(1,3,1)
plt.title("original")
plt.tight_layout()
plt.imshow(chest,cmap='gray')

# 2次元フーリエ変換
spectrum = np.fft.fft2(chest)
#spectrum = np.fft.fftshift(spectrum)# swap
# print(spectrum.dtype) # complex128（複素数）
real = np.real(spectrum)#実部を取得
imag = np.imag(spectrum)#虚部を取得
# 2DFFT結果（周波数空間画像）を表示
plt.subplot(1,3,2)
plt.tight_layout()
plt.imshow(real)
plt.title("real")
plt.subplot(1,3,3)
plt.tight_layout()
plt.imshow(imag)
plt.title("imaginary")

# 2次元フーリエ逆変換は関数でシンプルにコーディングできる。
invert = ifft2(spectrum)
# 逆変換結果を表示
plt.imshow(np.abs(invert), cmap='gray')#絶対値で取得する
plt.title("i-2D-FFT")
plt.show()#逆変換の表示

スワップの操作

高周波成分と低周波成分とを分けて見やすくする、あるいは、周波数フィルター（後述）目的で、計算結果はスワップされることがあります。（基底関数の並びがスワップされていたら、周波数空間のデータもスワップする必要があります。）

References

https://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm
C言語による画像再構成入門-フーリエ変換の基礎と応用　第4章
DIP Lecture 7: The 2D Discrete Fourier Transform, https://www.youtube.com/watch?v=NbQY1x8H6QQ