はじめに
CT画像のように断層像を積み重ねたような医用画像（ボリュームデータ）は，画素値を格納する3D配列に加えて，その配列のどの向きを正とするか という方向に関する情報を持ちます．
画素値を格納する配列の座標と，身体の向きを基準にして見る座標との関係は非常にややこしいため，備忘録としてここにまとめます．
なお，本投稿で使用している画像はすべて下記パブリックデータからのものです：
- https://www.cancerimagingarchive.net/collection/hcc-tace-seg/
- https://github.com/wasserth/TotalSegmentator?tab=Apache-2.0-1-ov-file
コンテンツ
1. 医用画像に関する3つの座標系
2. Image transformation
3. 補足：3DSlicerに表示される変換行列について
4. おわりに
5. ご参考
医用画像に関する3つの座標系
- World coordinate system
- Anatomical coordinate system
- Image coordinate system
https://slicer.readthedocs.io/en/latest/user_guide/coordinate_systems.html
1. World coordinate system
スキャナから見たときの座標系．
2. Anatomical coordinate system
日本語で患者座標系と言います．身体の向きについて，下記３面と6方向を定義します：
- Coronal面：地面に垂直な面．右手（R）⇔左手（L）方向のスライス．
- Sagittal面：地面に垂直な面．前（A）⇔後（P）方法のスライス．
- Axial面：地面に平行な面．つま先（I）⇔頭（S）方向のスライス．
（I,S,R,L,A,Pはそれぞれ，inferior, superior, right, left, anterior, posteriorの略）
どの向きを正とするかを表すときは，向かう方向のアルファベットを取って表します．
例えば，LPSなら
- 右→左（L）を正
- 前→後（P）を正
- つま先→頭（S）を正
とすることを意味します．ややこしいのは，どの向きを正とするかがデータやアプリによって異なる点です．DICOMやITKではLPSが採用されており，このLPSがスタンダードの位置づけですが，3DSlicerではRASが採用されています（詳細は後述しますが左右前後の座標の符号が反転することになります）．
この投稿では1ボリュームデータの基準点をoriginと言います．例えば，SimpleITKというpythonライブラリでoriginの座標を確認することができます．
import SimpleITK as sitk
img = sitk.ReadImage(r"3 C-A-P - acquisitionNumber 1.mhd")
print(img.GetOrigin(), "/mm")
(-224.80000299999992, -200.0, -374.9999999999999) /mm
データがmhd/raw形式なら，mhdをメモ帳などで開くと，"Offset"という項目で確認できます：
このボリュームがもしLPSなら，ボリューム内originの位置は下図のようになります：
3. Image coordinate system
ここではボリューム座標と呼ぶことにします．画素値が格納されている3D配列の座標（自然数）です．(0, 0, 0)がorigin点に一致します．
Image transformation
次にややこしいのが，ボリューム座標が必ずしもLPSになっているとは限らないということです．そこで，ボリューム座標から患者座標に変換する行列が必要になります．患者座標を$\vec{x}$，ボリューム座標を$(i, j, k)'$，変換行列を$A$，平行移動（origin）を$\vec{t}$としたとき，患者座標→ボリューム座標変換は
\vec{x} = A(i, j, k)' + \vec{t}
となり，成分ごとに書くと
\begin{pmatrix}
x \\
y \\
z \\
\end{pmatrix}
=
\begin{pmatrix}
A_{11}\Delta{i} & A_{12}\Delta{j} & A_{13}\Delta{k}\\
A_{21}\Delta{i} & A_{22}\Delta{j} & A_{23}\Delta{k}\\
A_{31}\Delta{i} & A_{32}\Delta{j} & A_{33}\Delta{k}\\
\end{pmatrix}
\begin{pmatrix}
i \\
j \\
k \\
\end{pmatrix}
+
\begin{pmatrix}
t_1 \\
t_2 \\
t_3 \\
\end{pmatrix}
となります1．ここで，$\Delta{i}, \Delta{j}, \Delta{k}$は解像度（spacing）を表します．
この式で出てくる$A_{11}, A_{12}, ... A_{33}$と$\Delta{i}, \Delta{j}, \Delta{k}$も，下記の通りSimpleITKで確認することができます：
import SimpleITK as sitk
img = sitk.ReadImage(r"3 C-A-P - acquisitionNumber 1.mhd")
print(img.GetDirection())
print(img.GetSpacing())
(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
(0.78125, 0.78125, 5.0)
データがmhd/raw形式なら，mhdをメモ帳などで開くと，"TransformMatrix"と"ElementSpacing"という項目で確認できます：
この場合は，変換行列に符号の反転がないため，ボリューム座標の向きと患者座標の向きが一致していてLPSということになります．
ここからは実際に3DSlicerを使って，ボリューム座標を変換した値が患者座標に一致するか確認してみます．
3DSlicerのDataモジュールを開くと，左下のData Proveでカーソルが置かれている患者座標とボリューム座標を確認することができます：
ここで注意しなければならないのは，Data Proveに出ている患者座標の見方です．今，(L 73.4, P 5.3, I 305.0)と出ていて符号がすべて正であることは，変換行列と全く関係がありません． カーソルを移動してみるとわかりますが，他の位置に行けばLPはひっくり返ってRAになります：
つまり，このData Proveで表示される患者座標は，必ず正を取るようにマイナスの場合はアルファベットを反転します．
今回，もとの画像は変換行列からLPSということがわかっているため，(L 73.4, P 5.3, I 305.0)と表示されている患者座標を(L 73.4, P 5.3, S -305.0)とみなして，ボリューム座標(382, 263, 14)→患者座標(L 73.4, P 5.3, S -305.0)の変換の検算をしてみます．
\begin{align}
\begin{pmatrix}
x \\
y \\
z \\
\end{pmatrix}
&=
\begin{pmatrix}
A_{11}\Delta{i} & A_{12}\Delta{j} & A_{13}\Delta{k}\\
A_{21}\Delta{i} & A_{22}\Delta{j} & A_{23}\Delta{k}\\
A_{31}\Delta{i} & A_{32}\Delta{j} & A_{33}\Delta{k}\\
\end{pmatrix}
\begin{pmatrix}
i \\
j \\
k \\
\end{pmatrix}
+
\begin{pmatrix}
t_1 \\
t_2 \\
t_3 \\
\end{pmatrix}\\
&=
\begin{pmatrix}
1\times 0.78125 & 0 & 0 \\
0 & 1\times 0.78125 & 0 \\
0 & 0 & 1\times 5.0 \\
\end{pmatrix}
\begin{pmatrix}
382 \\
263 \\
14 \\
\end{pmatrix}
+
\begin{pmatrix}
-224.80000299999992 \\
-200.0 \\
-374.9999999999999 \\
\end{pmatrix}
\end{align}
したがって，
\begin{align}
&L: 1 \times 0.78125\times 382 -224.80000299999992 = 73.637498\\
&P: 1\times 0.78125\times 263 -200.0 = 5.46875\\
&S: 1\times 5.0\times 14 -374.9999999999999 = -304.9999999999999
\end{align}
となり，多少誤差がありますがほぼ一致していて，考え方が正しいことが確認できました．
RASのデータについても確認してみます．
origin, 変換行列，spacingの情報を確認します：
import SimpleITK as sitk
img = sitk.ReadImage(r"ct.nii.gz")
print(img.GetOrigin(), "/mm")
print(img.GetDirection())
print(img.GetSpacing())
(157.68359375, 0.18359375, -869.0) /mm
(-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0)
(1.5, 1.5, 1.5)
変換行列から，ボリューム座標の正の向きと患者座標の正の向きが反転することがわかります．
originは下図のようになります：
この点について検算してみます：
- ボリューム座標：(165, 98, 25)
- 患者座標：(R 90.4, A 146.4, I 831.5) (=(L -90.4, P -146.4, S -831.5))
\begin{align}
&L: -1\times 1.5\times 165 +157.68359375 = -89.816407\\
&P: -1\times 1.5\times 98 +0.18359375 = -146.816406\\
&S: 1\times 1.5\times 25 -869.0 = -831.5
\end{align}
となり，確認できました．
補足：3DSlicerに表示される変換行列について
3DSlicerの表記はRAS基準になっているようです．例えば，DataモジュールのNode informationにある"Origin"はLPの符号が反転していますし，"IJKtoRASDirections"とあるように，名前の通りRASへの変換行列になっていて，その場合はもとが"IJKtoLPS"であるため変換行列のLPの符号が反転しています：
また，Segmentation editorからSegmentation geometryを確認すると，こちらもRAS基準で書かれていることがわかります：
おわりに
間違い等ありましたらコメントいただけると幸いです．
ご参考
- https://slicer.readthedocs.io/en/latest/user_guide/coordinate_systems.html
- https://medium.com/geekculture/coordinate-systems-in-medical-data-science-a-rant-90394f60b27
- https://nipy.org/nibabel/dicom/dicom_orientation.html?source=post_page-----90394f60b27--------------------------------
- https://theaisummer.com/medical-image-coordinates/?source=post_page-----90394f60b27--------------------------------#introduction-to-dicom-for-machine-learning-engineers
- https://www.mathworks.com/help/medical-imaging/ug/medical-image-coordinate-systems.html#
-
ご参考-3にあるように，4x4の行列にまとめて書くこともできます ↩