フーリエ変換。学生の頃からよく聞く単語であって、なんとなくでしか理解していなかった私ですが今回、MRIの画像のフーリエ変換をして分かったことを記事にします。
imageJでフーリエ変換
それでは、まずimageJで上の画像をフーリエ変換してみます。
imageJのツールバーに画像をドラックドロップしてツールバーの「Process」⇒「FFT」⇒「FFT」で以下のフーリエ変換した画像を取得することができます。
今回は、画像をフーリエ変換し再度元の画像に戻す。ここまでを目標としました。
Pythonでやってみた
まずはpydicomで画像を読み込み、ピクセルデータを取得します。
そのピクセルデータをnumpyを使ってフーリエ変換して表示してみます。
numpyでフーリエ変換は
numpy.fft.fft〇(配列)
〇の部分は２次元のデータをフーリエ変換する場合は２を入れます。3次元であれば3を入力。（1次元のデータの場合は省略できます。）
まずは、ここまでコードにしてみました。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
import pydicom import numpy as np import matplotlib.pyplot as plt img = pydicom.dcmread('IMG9') pix_arr = img.pixel_array # 変数pix_arrにピクセルデータを格納 fft_arr = np.fft.fft2(pix_arr) # 変数fft_arrにpix_arrをフーリエ変換したデータを格納 fig = plt.figure() ax1 = fig.add_subplot(1,2,1) ax2 = fig.add_subplot(1,2,2) ax1.imshow(pix_arr, cmap='gray') ax2.imshow(fft_arr, cmap='gray') plt.show() |
上記コードを実行すると・・・・エラーが
TypeError: Image data of dtype complex128 cannot be converted to float
データ型のエラー？「complex128」って？？？初めて見る型でしたが調べてみると複素数型らしいです。確かに、学校で習いました。フーリエ変換を行うと実数部と虚数部に分かれると。
フーリエ変換したデータを表示してみます。上記コードの9行目に
|
1 |
print(fft_arr) |
を追加して表示してみました。
確かに、複素数を現す「j」が付いたものが表示されています。
よく教科書でみる形ですね。
画像部分は実数の部分だという記憶があったので、フーリエ変換後のデータから実数部分だけを取り出して画像表示を試してみます。
|
1 |
ax2.imshow(fft_arr, cmap='gray') |
上記コードの fft_arr の部分を fft_arr.real に変換して再度挑戦。変換後は以下になります。
|
1 |
ax2.imshow(fft_arr.real, cmap='gray') |
あれ、右側の画像、真っ黒です・・・データがおかしいのかなと思ったのですが、右側の画像上でカーソルを動かすと右下の値が変化しているのでデータはありそうです。
データ型の確認をしてみました。
|
1 |
print(fft_arr) |
の次に以下のコードを足してデータの型を確認してみます。
データ型は 「float64」 の型でした。DICOM画像のピクセルデータは基本的に16ビット画像ですが、フーリエ変換後の実数部のデータは浮動小数点型64ビットとなっています。という事は、閾値を設定してあげれば表示されるかもしれない。
|
1 |
ax2.imshow(fft_arr.real, cmap='gray') |
の部分を
|
1 |
ax2.imshow(fft_arr.real, cmap='gray', vmax =50000) |
と変更して再度チャレンジ。結果は。。。
今度は真っ白になってしまいました。でも、やっぱりカーソルを合わせると値が右下に表示されているんでデータはあるのだと思うのですが・・・・
そういえば、フーリエ変換後のデータはfloat64なのでデータの幅が非常に広いのかなと思い、下限値も設定してみました。
|
1 |
ax2.imshow(fft_arr.real, cmap='gray',vmin = -50000, vmax =50000) |
なにやらうっすら見えてきましたが、何かが違う・・・・
なんでなのか悩んでいると、そうそう、MRI画像は０以下のデータは持たず、絶対値表示されているという事を思い出しました。
なので
|
1 |
ax2.imshow(fft_arr.real, cmap='gray', vmax =50000) |
の部分を
|
1 |
ax2.imshow(np.abs(fft_arr.real), cmap='gray', vmax =50000) |
に変更し絶対値表示に変更してみると
やっと、画像が確認できましたが・・・・・imageJでやったフーリエ変換の画像と違う・・・・・。
先ほどやったimageJのFFT画像では中心部が白くて周りに行くほど白くなる画像のはずですが・・・
調べてみると、位相シフトが必要らしいです。
numpy.fft.fftshift（配列）
で位相シフトができるようです。
|
1 2 |
fft_arr = np.fft.fft2(pix_arr) # 変数fft_arrにpix_arrをフーリエ変換したデータを格納 |
に部分を
|
1 2 3 |
fft_arr = np.fft.fft2(pix_arr) # 変数fft_arrにpix_arrをフーリエ変換したデータを格納 fft_arr_shift = np.fft.fftshift(fft_arr) |
と付け加えて、位相シフトさせます。そして
|
1 |
ax2.imshow(np.abs(fft_arr.real), cmap='gray', vmax =50000) |
の部分で位相シフトした配列に変更します。
|
1 |
ax2.imshow(np.abs(fft_arr_shift.real), cmap='gray', vmax =50000) |
おお！！なんかそれらしい画像が得られました！！再度、imageJでのフーリエ変換の画像を確認してみましょう。
画像のコントラストは違っていますが、それらしいものになりました！！ただ、これがきちんと正しいのか確認をしてみてみないことには終われません。
なので、再度今の工程を戻っていきます。フーリエ変換でできた画像を再度、位相シフトし、逆フーリエ変換して表示してみます。逆フーリエ変換は
numpy.fft.ifft〇(配列)
でできます。コードを以下に示します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
import pydicom import numpy as np import matplotlib.pyplot as plt img = pydicom.dcmread('IMG9') pix_arr = img.pixel_array # 変数pix_arrにピクセルデータを格納 fft_arr = np.fft.fft2(pix_arr) # 変数fft_arrにpix_arrをフーリエ変換したデータを格納 fft_arr_shift = np.fft.fftshift(fft_arr) print(fft_arr) print(fft_arr.real.dtype) fft_rev = np.fft.fftshift(fft_arr_shift) # 位相シフトを戻す fft_rev = np.fft.ifft2(fft_rev) #逆フーリエ変換 fig = plt.figure() ax1 = fig.add_subplot(1,3,1) ax2 = fig.add_subplot(1,3,2) ax3 = fig.add_subplot(1,3,3) ax1.imshow(pix_arr, cmap='gray') ax2.imshow(np.abs(fft_arr_shift.real), cmap='gray', vmax =10000) ax3.imshow(np.abs(fft_rev.real), cmap='gray') plt.show() |
結果は
一番右側の画像が、戻した画像になります。無事に元の画像に戻りました。
分かったこと
numpyを使ってピクセルデータをフーリエ変換した結果は実数と複素数の結果として返される。
フーリエ変換後の実数部分のデータ型はfloat64型でデータ範囲が広域になるので範囲を絞る必要がある。
フーリエ変換の結果を表示すると、低周波性部分は四隅にあり、位相シフトすることにより中心に低周波成分、四隅に高周波成分を配置することができる。（MRIの本とかに書いてあるsin波、cos波の位相シフトとはこのことを言っているのかな？）
最後に
今回は、フーリエ変換をやってみました。
いままでよくわからず、とりあえず苦手克服のためにやってみました。もしかしたら間違っているかもしれませんが、お気づきの際にはコメント等いただけると幸いです。
おつかれさまでした。
