MRI画像の画素値がしたがう確率分布をPythonで理解する

以下の有名な論文があります．

pubmed.ncbi.nlm.nih.gov

この論文によればMRI画像は

• 正規分布
• Rician (Rice) 分布
• Rayleigh分布

の３つのいずれかの確率分布にしたがうとされています．簡単に言えば，SNRが良いピクセルは正規分布，悪いピクセルはRice分布，信号がない箇所（バックグラウンド）はRayleigh分布にしたがうようです．

今回の記事では，MRI画像にノイズを付加した画像を作成し，上述した確率分布が得られるのかを検証しました．MRI画像はエムアールアイシミュレーションズさんのBlochSolverで作成した画像を用いています．

www.blochsolver.org

また，今回検証に使用した画像は，以下の以前の記事で作成したファントムデータのT2WIを使用しています．

radmodel.hatenablog.com

ノイズを付加した画像の作成

MRI画像はなんらかのノイズの影響を受けているため，確率分布にしたがってピクセル値は変動していると考えられています．ノイズの原因はさまざまですが，今回とりあげるのは熱雑音です．以下にwikipediaの熱雑音の説明を載せています．

熱雑音（ねつざつおん、英: thermal noise）は、抵抗体内の自由電子の不規則な熱振動（ブラウン運動[1]）によって生じる雑音のことをいう。

MRIの場合，熱雑音はコイルなどから発生していると考えられます．また，熱雑音は正規分布に従うと仮定されることが多いようです．

以前の記事で取り上げたように，MRIは直角位相検波という手法でreal信号とimaginary信号を取得しています．そのため，real信号とimaginary信号に対して，正規分布にしたがうノイズを付加することで，実際のMRI画像と同様のノイズを付加したことを模擬できます．以下，そのコードです．

import numpy as np 
import matplotlib.pyplot as plt 
from PIL import Image

path = 'result/2020-11-26/2D_SpinEcho-15-complex.flt'
matrix = 256

# numpyでバイナリデータを読みこむ
def bin2kspc(path, matrix):
    raw = np.fromfile(path, dtype='f', sep='')
    real = raw[::2].reshape(matrix, matrix)
    imag = raw[1::2].reshape(matrix, matrix)
    return real, imag

# real imageとimaginary imageに正規分布にしたがうノイズを付与
def add_noise(real, imag, sd):
    shape = np.shape(real)
    noise_r = np.random.normal(0, sd, shape)
    noise_i = np.random.normal(0, sd, shape)
    real += noise_r
    imag += noise_i
    return real, imag

# realとimaginaryを複素数としてひとまとめにする
def make_comp(real, imag):
    real = real.astype(np.complex64)
    imag = imag.astype(np.complex64)
    imag *= 1j
    kspace = real + imag
    return kspace

# 逆フーリエ変換
def ifft(kspace):
    ifft = np.fft.ifft2(kspace)
    ifft = np.fft.ifftshift(ifft)
    real_ifft = ifft.real
    imag_ifft = ifft.imag
    abs_ifft = np.abs(ifft)
    return real_ifft, imag_ifft, abs_ifft

def main(sd):
    real, imag = bin2kspc(path, matrix)
    n_real, n_imag = add_noise(real, imag, sd)
    kspace = make_comp(n_real, n_imag)
    real, imag, magn = ifft(kspace)
    return real, imag, magn

# sd=10とsd=50のノイズを付加したで画像を作成
real_sd10, imag_sd10, magn_sd10 = main(sd=10)
real_sd50, imag_sd50, magn_sd50 = main(sd=50)

# データを表示
fig = plt.figure(figsize=(7,4))

plt.axis('off')

ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

ax1.imshow(magn_sd10, cmap='gray',vmin=0,vmax=1)
ax1.set_title('Magnitude (sd=10)')
ax1.axis("off")
ax2.imshow(magn_sd50, cmap='gray',vmin=0,vmax=1)
ax2.set_title('Magnitude (sd=50)')
ax2.axis("off")

plt.savefig('magnitude.png')

k-space上に標準偏差が10と50の正規分布にしたがうノイズを加えたMRI画像を作成しました．sd=50の画像はかなりノイジーに見えます．

大量のデータの作成

次に同じsdの正規分布にしたがうノイズを発生させ，そのノイズを付加した画像を複数枚作成します．こうすることで，擬似的に複数回ファントムを撮像した状態を模擬できます．今回知りたいのはMRIの画素値がどのような確率分布にしたがうのかということなので，シミュレーションによりノイズが付加された画像を大量に作成します．

画像を作成した後，すべての画像上の左上（バックグラウンド）と画像中の3×3の丸がある中の2行目のすべての丸のなかの画素値をひとつ取得します．つまり４つのROI (1pixel) を設定します．下の画像のA,B,C,Dの位置です．

import pandas as pd
from tqdm import tqdm

images = {'sd10':[],'sd50':[]}
d = {'magn_sd10_A':[], 'magn_sd10_B':[], 'magn_sd10_C':[], 'magn_sd10_D':[],
     'magn_sd50_A':[], 'magn_sd50_B':[], 'magn_sd50_C':[], 'magn_sd50_D':[]}

for j in tqdm(range(5000)):
    images['sd10'].append(main(sd=10)[2])
    images['sd50'].append(main(sd=50)[2])

for i in tqdm(range(5000)):
    # A = [0,0], B = [80,140], C = [130,140], D = [180,140]
    d['magn_sd10_A'].append(images['sd10'][i][0,0])
    d['magn_sd10_B'].append(images['sd10'][i][80,140])
    d['magn_sd10_C'].append(images['sd10'][i][130,140])
    d['magn_sd10_D'].append(images['sd10'][i][180,140])
    d['magn_sd50_A'].append(images['sd50'][i][0,0])
    d['magn_sd50_B'].append(images['sd50'][i][80,140])
    d['magn_sd50_C'].append(images['sd50'][i][130,140])
    d['magn_sd50_D'].append(images['sd50'][i][180,140])

df = pd.DataFrame.from_dict(d)

100%|██████████| 5000/5000 [00:47<00:00, 104.65it/s]
100%|██████████| 5000/5000 [00:00<00:00, 29689.48it/s]

5000回撮像したファントム画像のA,B,C,Dの位置の画素値のデータを取得しました．

次にfitterというPythonのライブラリを使って，それぞれのピクセルがどの確率分布にしたがうのか検証します．fitterは最尤法を用いて，ピクセルの集合のヒストグラムから当てはまりの良い確率分布を探し出すことができるライブラリです．

from fitter import Fitter

result = {'magn_sd10_A':[], 'magn_sd10_B':[], 'magn_sd10_C':[], 'magn_sd10_D':[],
          'magn_sd50_A':[], 'magn_sd50_B':[], 'magn_sd50_C':[], 'magn_sd50_D':[]}

for i in result.keys():
    result[i] = Fitter(df[i], distributions=['norm', 'rice', 'rayleigh'])

以上でsd=10とsd=50のノイズを加えた画像のA,B,C,Dの位置のピクセルがしたがう確率分布の結果を計算できました．

それではMRI画像がどのような確率分布にしたがうのかを見ていきます．まずはAの位置（バックグラウンド）のsd=10とsd=50のバックグラウンドのピクセルがしたがう確率分布を見ていきます．

# sd=10
result['magn_sd10_A'].fit()
result['magn_sd10_A'].summary()

	sumsquare_error	aic	bic	kl_div
rayleigh	70.486597	-137.259822	-21291.818686	inf
rice	70.968561	-130.667313	-21249.229462	inf
norm	283.554135	-23.668768	-14331.916283	inf

# sd=50
result['magn_sd50_A'].fit()
result['magn_sd50_A'].summary()

	sumsquare_error	aic	bic	kl_div
rice	2.401464	129.260283	-38180.021089	inf
rayleigh	2.406282	124.249910	-38178.517287	inf
norm	10.665122	209.537389	-30734.037491	inf

表は左列から二乗和誤差，AIC，BIC，カルバック・ライブラー情報量を示しています．どの指標も値が小さいほど良いことを示します．グラフはヒストグラムがデータの分布で，曲線はデータを確率分布にフィッティングさせたときに求めたパラメータで作成された確率分布を示しています．

sd=10ではRayleigh分布の二乗和誤差がもっとも低値でしたが，sd=50ではRice分布の二乗和誤差がもっとも低値を示しました．よく見てみると，sd=50ではAICとBICはRayleigh分布の方が低値を示しているので，実質的にどちらのsdの画像もRayleigh分布にしたがっていると考えることができそうです．

次にBの位置を見ていきます．

# sd=10
result['magn_sd10_B'].fit()
result['magn_sd10_B'].summary()

	sumsquare_error	aic	bic	kl_div
norm	20.220786	-7.820432	-27535.376053	inf
rice	20.230288	-5.692187	-27524.509975	inf
rayleigh	770.452108	-203.635577	-9334.044093	inf

# sd=50
result['magn_sd50_B'].fit()
result['magn_sd50_B'].summary()

	sumsquare_error	aic	bic	kl_div
norm	1.450872	235.523199	-40708.108874	inf
rice	1.451824	228.592473	-40696.311037	inf
rayleigh	16.732950	115.898379	-28482.032280	inf

こちらは表の値とヒストグラムをみた感じでは，どちらも正規分布にしたがっていると考えて良さそうです．

次にCの位置です．

# sd=10
result['magn_sd10_C'].fit()
result['magn_sd10_C'].summary()

	sumsquare_error	aic	bic	kl_div
rayleigh	70.767290	-194.946403	-21271.947146	inf
rice	70.774084	-192.916473	-21262.949908	inf
norm	274.531438	-109.723563	-14493.602664	inf

# sd=50
result['magn_sd50_C'].fit()
result['magn_sd50_C'].summary()

	sumsquare_error	aic	bic	kl_div
rayleigh	3.147475	139.233666	-36835.929129	inf
rice	3.147479	141.245974	-36827.405657	inf
norm	13.302968	226.980239	-29628.995801	inf

Cの位置では，どちらのsdの画像もRayleigh分布もしくはRice分布にしたがっているように見えます．なお，DはBほぼ同様の結果になったので省略します．

まとめ

バックグラウンド領域はRayleigh分布に，信号発生領域では正規分布もしくはRice分布にしたがっていることが確認できました．最初にあげた論文では，バックグラウンド領域はRayleigh分布にしたがうとありましたので，その通りになっていることが確認できました．

気になるのは信号発生領域はどのようなときに正規分布もしくはRice分布にしたがうのかということです．論文中ではSNR<2のピクセルはRice分布，SNR>2のピクセルは正規分布にしたがうとありました．そこで今回のBとCのSNRを算出してみます．

def SNR(d):
    return np.mean(d)/np.std(d)

print(f"SNR of B (sd=10): {SNR(df['magn_sd10_B'])}\nSNR of B (sd=50): {SNR(df['magn_sd50_B'])}\nSNR of C (sd=10): {SNR(df['magn_sd10_C'])}\nSNR of C (sd=50): {SNR(df['magn_sd50_C'])}")

SNR of B (sd=10): 12.391050217474877
SNR of B (sd=50): 2.8108878515892477
SNR of C (sd=10): 1.9194703847218806
SNR of C (sd=50): 1.9089855205754447

BはどちらのsdにおいてもSNR>2でしたが，CはどちらのsdにおいてもSNR<2でした．論文と整合性がとれてることが確認できました．

参考にした記事

rmizutaa.hatenablog.com