画像と周波数領域と1/fゆらぎとsaliency

Posted at 2022-12-02

はじめに

この記事は、LIGHTzアドベントカレンダー2022の2日目の記事です。

@KTasknこと、たすく、です。
LIGHTzは業務に支障がでない範囲で自由な働き方ができる風土が整っています。
私は現在、そんな弊社で働きながら大学院のの修士１回生(２年目(!?))として、データサイエンス・コンピュータビジョンの研究に励んでいます。

というのは，去年のアドベントカレンダーの書き出しをそのままコピってきたのですが，

今年のアドベントカレンダーでは，そのコンピュータビジョンを研究・勉強するうえで，おもしろかった内容を投稿していきたいと思っています．

画像のフーリエ変換

みなさん，フーリエ変換はしたことありますか？毎日していますか？
関数を周波数領域に変換するアレです．フーリエ変換だけで教科書一冊ある奥の深そうな話なのでその話はすみませんが割愛します．
ウィキペディアの周波数領域の記事に記載の記事がとてもわかりやすい．

いろいろな関数は，単純な波の足し合わせで表現できるんですね．

上の例だと，１次元な値なんですが，実は２次元の画像でも同じことができるというのが今回のお話でして，
下の画像は私が昨夜ぽちぽちとドット絵エディタで打ったMegamanのドット絵になります．これだって実は波の足し合わせで表せるんですよ．

二次元の波ってなんじゃという話なんですが，上の画像を構成する波が下の画像ら（抜粋）．(幅，高さ)=(32,32)の画像なら32x32=1024の波の足し合わせで表現されます．

PythonとOpenCVを使うと下記のコードで確認することができます．

# jupyter notebook
from PIL import Image
import numpy as np
import cv2
import matplotlib.pyplot as plt

# 画像を開いて，グレースケールにする
a_image = Image.open("megaman.png").convert("L")
# 画像を表示
display(a_image)

# 高速フーリエ変換する
fimage = np.fft.fft2(a_image)

# 周波数領域を描画
display(Image.fromarray(np.uint8(np.abs(fimage))))

images = []
for i in range(32):
    for j in range(32):
        # 周波数ごとに波形をとりだす
        # 0で初期化した周波数行列を作成する
        d_fimage = np.zeros((32, 32), dtype=complex)
        # 周波数帯 i, j をとりだし，行列に埋め込む
        d_fimage[j, i] = fimage[j, i]
        # 周波数領域 -> 時間領域(画像)に変換する
        d_ifimage = np.fft.ifft2(d_fimage)
        # 表示する
        display(Image.fromarray(np.uint8(d_ifimage)))
        # あとで足し合わせるためにリストに保存しておく
        images.append(d_ifimage)

# リストに保存した波の画像を足し合わせて表示する
display(Image.fromarray(np.uint8(np.abs(np.sum(images, axis=0)))))

# 時間領域に再変換
ifimage = np.fft.ifft2(fimage).real
display(Image.fromarray(np.uint8(ifimage)))

風景と被写体と，周波数領域と1/fゆらぎ

そんな周波数領域なのですが，不思議なことに自然な写真だと，周波数とパワーが綺麗に反比例することが知られており，世の中で$\frac{1}{f}$なんていわれる法則に従います．下の画像だと，オレンジ色の線が実際の周波数に対するパワーの値で，青い線が$\frac{1}{f}$のガイドになります．かなり線にのっかっているのがわかります．