【Python初心者向け】医療画像（MRI/CT/PET等）を読み込んで任意に回転させた画像を出力する

こんにちは、OguLinks代表の小倉（@jun_ogulinks）です。

コメディカル向けにPythonを使ったプログラミングでこんなことできるよ的な情報を発信しています。

プログラミング未経験者・初心者にも「へー、そうなんだぁ」と思ってもらえるように、あまり突っ込み過ぎずにお話していきます。

今回のテーマはこちら！

1. Pythonで医療用画像（3次元）を読み込む方法

2. 3次元の画像を回転する方法

3. 回転後の画像を出力・保存する方法

では、さっそく解説です。

ちょうどこの記事を書くきっかけになったのは、GoogleがOpenAIのChatGPTに対抗して公開したBardの精度がどれ程なのか知りたかったからです。

実際、この記事のタイトルにある内容をBardにお願いしてみた結果がこちらです。

初心者の方には呪文にしか見えませんが、今は気にしないでください。

NiBabelというのは、Pythonで医療画像を読み込むのに一般的に使われるパッケージです。

パッケージというのはPythonで0から自分でコードを書かなくても、特定の機能がまとめてあるもののことを言います。

この自動生成されたコードはちゃんと動くのでしょうか？

試すとこのようになります。

( ,,`･ω･´)ﾝﾝﾝ？

なんかおかしいですね。

roll方向に90度回転して欲しいところ、なぜか180度も回転しちゃってます。

該当するコードはこの部分（右側にグレーでコメントされている行）。

import nibabel as nib

# 入力画像の読み込み
img = nib.load('input.nii.gz')

# 画像のロール方向に90度回転
affine = img.affine.copy()
affine[:, 2] = -affine[:, 2]                  # ← この部分が回転処理を表しています
img = nib.Nifti1Image(img.get_data(), affine)

# 画像の保存
nib.save(img, 'output.nii.gz')

回転・平行移動はアフィン変換によって行われます。

読み込んだ画像のヘッダー情報からアフィン行列を抜き出してきて、それに処理を加えるということをここではしているのですが、これが間違っています。

Bardで指示を出すときにプロンプトをもっと工夫すれば解決できるかもしれませんが、今回自動生成されたコードでは正直、使いにくいですね。

なので、このコードをベースに自分ならこうやって書く！というのを紹介します（初心者の方は　(　´_ゝ｀)ﾌｰﾝ　くらいでOK）。 

# 必要なパッケージの読み込み
import nibabel as nib
import numpy as np
import math
from nilearn import plotting
import matplotlib.pyplot as plt

# 回転処理用の独自関数を定義
def rotate(theta_x, theta_y, theta_z, original_affine):
    theta_x = theta_x
    theta_y = theta_y
    theta_z = theta_z

    # x軸を中心とした回転の座標回転行列
    r_mat_x = np.array([[1, 0, 0, 0],
                        [0, math.cos(math.radians(theta_x)), -math.sin(math.radians(theta_x)), 0],
                        [0, math.sin(math.radians(theta_x)), math.cos(math.radians(theta_x)), 0],
                        [0, 0, 0, 1]])

    # y軸を中心とした回転の座標回転行列
    r_mat_y = np.array([[math.cos(math.radians(theta_y)), 0, math.sin(math.radians(theta_y)), 0],
                        [0, 1, 0, 0],
                        [-math.sin(math.radians(theta_y)), 0, math.cos(math.radians(theta_y)), 0],
                        [0, 0, 0, 1]])

    # z軸を中心とした回転の座標回転行列
    r_mat_z = np.array([[math.cos(math.radians(theta_z)), -math.sin(math.radians(theta_z)), 0, 0],
                        [math.sin(math.radians(theta_z)), math.cos(math.radians(theta_z)), 0, 0],
                        [0, 0, 1, 0],
                        [0, 0, 0, 1]])

    # 回転後のpointsを格納するための配列を生成
    tmpMat = np.array([[0,0,0,0],
                       [0,0,0,0],
                       [0,0,0,0],
                       [0,0,0,1]])

    # x軸回転後
    if theta_x!=0:
        tmpMat = np.dot(r_mat_x, original_affine)

    # y軸回転後
    if theta_y!=0:
        if theta_x!=0:
            tmpMat = np.dot(r_mat_y, tmpMat)
        else:
            tmpMat = np.dot(r_mat_y, original_affine)

    # z軸回転
    if theta_z!=0:
        if theta_x!=0 | theta_y!=0:
            tmpMat = np.dot(r_mat_z, tmpMat)
        else:
            tmpMat = np.dot(r_mat_z, original_affine)

    return tmpMat

#------------------------------------
# Bardで自動生成されたコードを少し修正

# 入力画像の読み込み
img = nib.load('input.nii')

# 読み込みデータからヘッダーのコピーを作製
header = img.header.copy()

# 画像のロール方向に90度回転
affine = img.affine.copy()
new_affine = rotate(0, 90, 0, affine)               # ← この部分がメインの変更点
convert_img = nib.Nifti1Image(img.get_fdata(), new_affine, header)

# 画像の保存
nib.save(convert_img, 'output.nii')


#------------------------------------
# オリジナル（回転前の）画像を描画
before_img = nib.load('input.nii')

fig = plt.figure(figsize=(10,5))
plotting.plot_anat(before_img, figure=fig)


#------------------------------------
# 回転後の画像を描画
after_img = nib.load('output.nii')

fig = plt.figure(figsize=(10,5))
plotting.plot_anat(after_img, figure=fig)

いきなり長いコードでビビるかもしれませんが、これは使いやすくするために必要なものになります。

というわけで、roll方向に90度回転できているか確認します。

できていますね、問題なさそうです。

自動生成されたコードでは回転処理のところをアフィン行列の一部にマイナスをかけていましたが、それだと「どのように回転しているのか」イメージしずらいですよね？

なので、私の作ったコードではx,y,z軸のそれぞれを何度回転するか指定できるように改造してあります。

この3行目の部分ですね。

# 画像のロール方向に90度回転
affine = img.affine.copy()
new_affine = rotate(0, 90, 0, affine)               # ← この部分がメインの変更点
convert_img = nib.Nifti1Image(img.get_fdata(), new_affine, header)

rotate(x軸の回転角度、y軸の回転角度、z軸の回転角度、元画像のアフィン行列) という形式になっています。

rotate関数という独自の関数を定義することでこのようなことができるようになります（この関数を定義している分、コードが長くなっているのですね）。

もちろんこれには3次元のアフィン行列の変換方法についての知識も必要になってきますが、それは必要になったときに調べれば良いです。

正直、私も最初はそんなこと知りませんでした。

「3次元画像の回転ってどうやるんだろう？」と思ってネットで公式を調べて、それをPythonに実装しただけです。

なので、初心者の方でも調べる力があればこういったことにも挑戦できます。

そもそも、コメディカルで3次元画像の回転について知っているのはごくわずかだと思います。

私は臨床検査技師の資格を持っていますが、そんなことは一切授業で出てきません。

ちなみに、ここで紹介したコードはどんなときに役に立つか、最後に説明して終わりにしたいと思います。

データ分析をする際に、使用する機種ごとに一定のずれが生じることがあります。

ヒトの場合はほとんどないかもしれませんが、動物用の装置の場合は機種によって意味わからないくらい座標が変わることがあります。

そのため、回転・平行移動がある程度一定になっている場合は、3D Slicerのようなビュワーで表示する前に大まかな位置調整をこのコードですることができます。

ニッチな用途かもしれませんが、私は仕事の都合上これすごい使っています。

Python初心者の場合、ここに書かれている情報だけで同じことができるかと言われると難しいですが、ひとまずこんなことがPythonでできるんだと思っていただけたのであれば幸いです。