ImageJのフーリエ変換を雑にマクロで紐解いてみる

私は画像のフーリエ変換をするときにImageJを使用することがあります．いちいちコードをかかずとも，マウスを数回クリックするだけで複雑な計算をしてくれるImageJはすばらしいソフトウェアだと思います．プログラムのコードを書くことは好きですが，画像一枚をフーリエ変換するようなときはGUIで行ったほうが効率的です．

ImageJのデフォルトの設定でフーリエ変換を行うと画像が一枚作成されるのですが，僕はこれが何なのか知りませんでした．よくない使い方ですが今まで雰囲気でなんとなくフーリエ変換を行っていました．以前の記事の最初に示したなんちゃってk-spaceもImageJのFFTを使用して作成しました．

radmodel.hatenablog.com

知らないままにしておくのはまずいと思い，ImageJのデフォルト設定のFFTについて勉強しました．その結果，デフォルトのFFT後に現れる画像はパワースペクトルであり，画像としては不正確だが見た目を整えて人間が見やすい形にした画像であることがわかりました．表示されていないところで正確なデータが保持されており，計算をする際にはそちらの値を使用するそうです．そのためデフォルトのFFT後に表示される画像を保存したとしても，他の解析には使えないということがわかりました．

この記事ではImageJのFFTのデフォルト設定で行われている操作をマクロを使って雑に検証していきます．雑にというのはプログラムのソースコードを読めばわかるのに，~~僕はJavaが分からないので~~公式ドキュメントを参考にして自分なりに検証したからです．そのため計算が正確でないかもしれないし，どこか間違っている可能性もあります．その点をご了承ください．

ImageJのFFTはFHTがベースとなっている

そもそもImageJのFFTは純粋な高速フーリエ変換ではないようです．ImageJ日本語情報のFFTの項目を見てみると以下の記載があります．

この機能は、NIHImage の派生ソフト ImageFFT の作者である Arlo Reeves が製作した2次元の高速ハートレー変換（Fast Hartley Transform; FHT）を実装したものがベースになっている。

ImageJのFFTは高速ハートレー変換をベースにしているそうです．高速ハートレー変換なんてはじめて聞きました．

高速ハートレー変換（Fast Hartley Transform; FHT）

では高速ハートレー変換とはなんなのでしょう．日本語で以下の文献を見つけました．

www.jstage.jst.go.jp

この文献の緒言から高速ハートレー変換について簡単にまとめてみると

フーリエ変換と極めて類似した変換である
実数の変換である
順・逆変換が等しい

といった特徴があるようです．（離散化）ハートレー変換，ハートレー逆変換を式で確認してみると

f:id:radmodel:20201011044812p:plain

f:id:radmodel:20201011044710p:plain

です．ここでcas(x)は

f:id:radmodel:20201011044438p:plain

らしいです．

一方，（離散化）フーリエ変換，フーリエ逆変換は以下になります．

f:id:radmodel:20201011045158p:plain

f:id:radmodel:20201011050147p:plain

たしかに上で挙げたような特徴が式からうかがえます．~~よくわかりませんが，~~ImageJのFFTはハートレー変換がベースとなっているようです．

ImageJのFFT後の画像はパワースペクトル

公式ドキュメントを読み進めていくと，ImageJのFFT後の画像はパワースペクトルであるとの記載があります．さらに，このパワースペクトルは対数変換されており，8bitになっているようです．これらを理解するために，デフォルトのFFTで表示される画像（パワースペクトル）をマクロで作成していきます．

複素フーリエ変換

実数の2次元画像に対してフーリエ変換を行うと，複素数の2次元画像に変換されます．元画像f(n)をフーリエ変換して算出されるF(k)の右辺にはiが存在しているため，元画像f(n)が実数であっても複素数に変換されるということです．上のフーリエ変換の式のf(n)が画像でいうところの各画素値にあたり，F(k)がフーリエ変換後の各画素値になります．

複素数はそのままだと扱いづらいですが，実部と虚部に分けることで実数平面の2つの画像として表示することができます．まずは画像をフーリエ変換し，実部と虚部に分けて表示させます．ImageJを開きctrl+shift+Nを同時に入力，その後表示されたText Windowに以下のマクロを入力し，ctrl+Rで実行します．

run("T1 Head (2.4M, 16-bits)");
setSlice(64);
run("FFT Options...", "complex do");

このようにすることでImageJのサンプル画像（T1WI head sagital）を表示し，そのフーリエ変換後の画像を実部と虚部に分けて表示することができます．

f:id:radmodel:20201016003055p:plain

f:id:radmodel:20201016005043p:plain

ここでImageJのFFTはFHTベースと書きましたが，なぜ複素数の実部と虚部の値が出て来るのかが不思議だったのですが，wikipediaを見てみますと正規化の係数として1/2^1/2を用いる時のフーリエ変換とハートレー変換は以下の式で表されるそうです．

f:id:radmodel:20201017141117j:plain

つまりハートレー変換はフーリエ変換に変換できるようです．

パワースペクトル

次にこれらの実部と虚部の画像から理論上のパワースペクトルを作成します．パワースペクトルはF(k)の絶対値の二乗であらわされます．つまり実部と虚部の二乗和になります．マクロは以下の通りです．

selectWindow("Complex of t1-head.tif");
run("Square", "stack");
run("Stack to Images");
imageCalculator("Add create 32-bit", "Real","Imaginary");

なおここで作成した理論上のパワースペクトルはFFT OptionsのなかのRaw Power Spectrumを選択すれば作成可能です．手動で作成したものと，Raw Power Spectrumで作戦される画像は一致していました．

f:id:radmodel:20201016013506p:plain

対数変換と8bit画像への変換

ImageJのデフォルトのFFTではこの理論上のパワースペクトルの対数を取り，最小値が1、最大値が254になるようになんらかの処理を行なっています．そして32bit画像であったパワースペクトルを8bit画像に変換しているようです．

この処理が公式ドキュメントを読むだけではいまいち理解できませんでした．仕方ないので読めないソースコードをなんとなく雰囲気で読み，こちらの質問も参考にしてみました．

//対数変換
selectWindow("Result of Real");
run("Log");
//(1～254に)線形変換
//この画像の最大値が30.413、最小値が6.18
run("Macro...", "code=v=(v-6.18)/(30.413-6.18)*(254-1)+1"); 
//8bitに変換
run("Conversions...", " ");
run("8-bit");

f:id:radmodel:20201016013830p:plain

デフォルトのFFT後のパワースペクトルと比較

デフォルトのFFTで作成された画像とヒストグラムを表示し，比較してみます．

//ImageJのデフォルト設定のFFTを実行
selectWindow("t1-head.tif");
run("FFT Options...", "fft do");
//それぞれのヒストグラムを作成
selectWindow("Result of Real");
run("Histogram", "bins=256 x_min=0 x_max=255");
selectWindow("FFT of t1-head.tif");
run("Histogram", "bins=256 x_min=0 x_max=255");

f:id:radmodel:20201016013940p:plain

f:id:radmodel:20201016014146p:plain

ヒストグラムは左側が手動，右側がデフォルトのパワースペクトルのものです．デフォルトのフーリエ変換後の画像は手動で計算した画像と類似しており，ヒストグラムの形状も同様に類似していますが，値は同じになりませんでした．どこかで計算間違いをしているからだと思いますが，原因はわかりません．誰か教えてください．

まとめ

デフォルトのパワースペクトルを手動で（なんとなく）計算できました．ImageJのデフォルトのFFTがいろいろな処理をしていることがわかりました．

なぜImageJのFFTがこのような仕様なのかを自分なりに考察してみました．まず，FFTをパワースペクトルで表示することについてですが，パワースペクトルが最も適した形式であったからだと思いました．そもそもフーリエ変換は周波数解析をするために行われるので，各周波数の成分の大きさを観察するためにはパワースペクトルで表示するのが良いはずです．また，複素数を2次元画像として表現するのが困難なことも理由の一つだと思います．

ちなみにパワースペクトルがパワーと呼ばれる所以を調べてみたところ，以下の説明を見つけました．なんとなくわかった気になれます．

物理学のいろいろな問題で振動のエネルギーが1/2(定数)×(振幅)²で表されます．この定数は問題によって，物体の質量であったりコンデンサーの静電容量（キャパシタンス）であったりコイルの誘導係数（インダクタンス）であったりします．ですから，振幅の二乗をエネルギーとみなすのは自然です．しかし，ここでは物理学でいうエネルギーというより，形式的に考えた”エネルギーのようなもの”と考えてください．（金谷健一（2003）『これなら分かる応用数学教室最小二乗法からウェーブレットまで』共立出版，p.99）

その他の処理についてですが，対数を取るのは小さな値を底上げする（観察しやすくする）ためであり，線形変換をして8bitにするのは人間の目が8bit程度しか認識できないからそうしているのだと思います．

つまり，ImageJのデフォルトのFFTはわざわざ丁寧に人間が目で見やすいような形に調整してくれているのだと思います．

まとまりがないですが，最後に今回使ったImageJのマクロを以下にまとめて終わりにしたいと思います．

//サンプル画像を読み込み
run("T1 Head (2.4M, 16-bits)");
setSlice(64);

//(複素)フーリエ変換
run("FFT Options...", "complex do");

//パワースペクトルの算出
selectWindow("Complex of t1-head.tif");
run("Square", "stack");
run("Stack to Images");
imageCalculator("Add create 32-bit", "Real","Imaginary");

//ImageJのデフォルト設定のFFTから算出される画像に寄せる
selectWindow("Result of Real");
//対数変換
run("Log");
//(1～254に)線形変換
//この画像の最大値が30.413、最小値が6.18
run("Macro...", "code=v=(v-6.18)/(30.413-6.18)*(254-1)+1");
//8bitに変換
run("Conversions...", " ");
run("8-bit");

//デフォルトのパワースペクトルと作成したパワースペクトルのヒストグラムを比較
//ImageJのデフォルト設定のFFTを実行
selectWindow("t1-head.tif");
run("FFT Options...", "fft do");
//それぞれのヒストグラムを作成
selectWindow("Result of Real");
run("Histogram", "bins=256 x_min=0 x_max=255");
selectWindow("FFT of t1-head.tif");
run("Histogram", "bins=256 x_min=0 x_max=255");
````

放射線技術×データサイエンス

診療放射線技師がPythonとRとImageJを使って色々します