DICOM画像の再構成
任意の平面や曲面、3次元表示など、3次元の画像再構成は、放射線医療に欠かせないものとなっています。
その昔、フィルム時代、画像を取得後に、何らかの処理を行う「後処理（post-processing）」が自由にできなかった時代と比べて、デジタル画像が主流となった今はできるのが当然と思われるほどに画像再構成の存在感が増しています。
画像再構成を含むPost-processingは、DICOMデータオブジェクトでデジタルに画像が収集されない限り、実運用は不可能と言っていいでしょう。
CTスキャンのようなデジタル画像について考えてみましょう。
DICOMは、ピクセルデータだけでなく、距離、座標、患者を基準にした画像の面の方向を記録しています。
(0028, 0030) Pixel Spacing 属性：ボクセルの縦横サイズ（mm）と、画像の幅と高さがピクセル単位でわかっていれば、画像の物理サイズ ( Field Of View; FOV ) を求めることができます。例えば、Row, Columnが512、ピクセルのサイズが縦と横でぞれぞれ 0.4 mm のとき、FOVは縦横共に204.8mmとわかります。
(0020, 0032) Image Position (Patient) 属性：画像の左上隅（最初のピクセル）のx, y, z座標です。単位はmm。
(0020, 0037) Image Orientation (Patient) 属性：画像の行ベクトル（vr）と列ベクトル（vc）のx, y, z成分からなる方向余弦です。これらは、Image Position (Patient)属性で示された画像位置点が始点になります。画像の面の向きを定義するために利用されます。
(0018, 0088) Spacing between Slices 属性：連続するスライスの、そのスライイス方向における、ボクセル間の物理的な距離です。スライス方向の実際の距離の測定を可能にします。例えば、スライス方向がZ方向の場合、スライス画像間の間隔が3.0 mmであることが事前に分かれば、10スライス画像分の距離は、( 10 -1 ) × 3 = 27 mmとわかります。逆に、1番目と2番目のスライス画像で同じ位置にあるボクセルを考える場合、この2つのボクセルの間の距離がこの（0018, 0088）属性の値と同じになります。
その他：例えば、MRIの灌流画像を用いる解析では、時間軸で管理されたシリーズ画像から情報を処理する必要があるため、画像を取得した時刻が重要になったります。
この解説では最終的に画像平面をスライスしていきます。そのために必要な情報を早歩きで確認していきます。
DICOMの座標系
DICOM画像は現実世界を離散化した3次元空間から作成されます。この空間の軸は、x, y, z 軸で表現されます。DICOMではこの座標システムを、「Reference Coordinate System; RCS」と呼んでいます。
RCS は、DICOM Frame of Referenceにおける空間座標系である。
RCS は DICOM フレームの空間座標系であり、直交座標系における Image IE の原点、方向、空間スケールを決定するものである。
RCS は右手系の直交座標（right-handed Cartesian coordinate
）で、正の x 軸に沿った単位ベクトルと正の y 軸に沿った単位ベクトルのベクトル積は、正の z 軸に沿った単位ベクトルと等しくなる。
単位長さは1ミリメートルである。
一般的に、Image IEは、RCSの直交空間領域に対する画像サンプルの関係を指定する空間マッピングを含む。
DICOM Frame of Referenceは、一回の検査で複数のシリーズが作成された場合に、シリーズ間の空間的な位置を共有するための概念です。Image IEは画像そのものだと思ってください。
あとは、右手系？ってなりますよね（私だけでしょうか）。
これは、right-handed ruleというもので、右手の法則と同じ考え方です。どの指を z 軸にするかでパターンが分かれます。
（必ず x, y, z 軸は直交する）
一般に、DICOMでは、寝台（患者の体）方向を z 軸として考えます。そして、岸部〇伴が体操で使いそうな奇妙な右手の法則で軸を捉える必要があります（パターン１の形で、手首を後屈するのです）。
まず、CTスキャナのガントリと向き合って立ち、そして、自分の右手をまるでスパイダーマンが手首から蜘蛛の糸を出すようにガントリに向けます（掌が上）。そして、中指は寝台方向へ、人差し指は下方向、親指は右に向けます（ババアアア～ン！ゴゴゴゴゴゴという奇妙な描写が思い浮かびます）。
このようなイメージで、いいでしょう。このように患者を基準にして軸を考えます。この図でいうと、x 軸は患者の左側に向かって増加、y 軸は患者の後方に向かって増加、z 軸は患者の頭部に向かって増加します。
実は、この座標系は二足歩行（BIPED）な生き物を対象とするか、四足歩行（QUADUPED）動物を対象とするかで扱いに注意が必要になりますが、基本的にはこの座標系をもとに考えていくことになります（Anatomical Orientation Type (0010,2210)を参照してください）。
画像を表示する際には、患者の向きを考慮しなければならないことがあります。例えば、横向きで寝ている場合や、足がガントリ側に向いているなどです。患者の向きによって、画像を表示する際に、アプリケーションで画像を回転させたり、順序を入れ替えたりします。これは、本解説の範疇を超えるので、これ以上解説しません。ここでは、BIPED、かつ、頭部がガントリに入っていて、仰向けな状態、を対象として話を進めます。
画像の位置
DICOM画像は、画像の左上のピクセルが、RCS空間のどこに位置しているかを記録しています。これは、(0020, 0032) Image Position (Patient) 属性に保存されており、画像の左上隅（最初のピクセル）のx, y, z 座標になっています。単位はmmです。
例：
0020,0032 Image Position (Patient): -105.89111931994-104.02358965575-25.117109510232
方向余弦（Direction Cosine）
方向余弦（Direction Cosine）というのは、RCS基準軸からのベクトルの角度のことです。
(0020, 0037) Image Orientation (Patient)属性は、画像の面の方向を、方向余弦を、Row方向とColumn方向それぞれで記録しています。扱う角度の単位はラジアンです（度数法ではないです）。
例えば、一番単純な直交横断面を考えてみましょう。このとき、Image Orientation (Patient)属性は、次のようになります。
1 \ 0 \ 0 \ 0 \ 1 \ 0
これをRow方向成分、Column方向成分に分解すると次のようになります。
Row方向成分：[ 1 , 0 , 0 ]（最初の３つ）
Column方向成分：[ 0 , 1 , 0 ]（後ろの３つ）
これは、次の通りの意味です。
Row方向ベクトル：cos α = 1, cos β = 0, cos γ = 0
Column方向ベクトル：cos α = 0, cos β = 1, cos γ = 0
方向余弦はラジアンで計算されるので、度数法に直して意味を調べることができます。アークコサインで逆算して確かめます。アークコサインもラジアンで計算されるので、まず、アークコサインを計算した後に、その結果を角度に戻すと、次のようになります。
Row方向成分：
arc_cos α ( 1 ) to degree = 0
arc_cos β ( 0 ) to degree = 90
arc_cos γ ( 0 ) to degree = 90
Column方向成分：
arc_cos α ( 1 ) to degree = 90
arc_cos β ( 0 ) to degree = 0
arc_cos γ ( 0 ) to degree = 90
これで想像しやすくなりました。
Row方向は、RCS空間で、x 軸から0度、y 軸から90度、z 軸から90度なベクトルです。つまり、x 軸の方向と同じですね。
同様に、Column方向は、RCS空間で、x 軸から90度、y 軸から0度、z 軸から90度なベクトルです。つまり、y 軸の方向と同じです。
なので、この画像はXY平面を表している、ということが分かります。このように、画像の面の向きが計算できるようになっています。
3D空間内のボクセル位置の計算
画像の位置、ボクセルのサイズ（Pixel SpacingとSpacing between Slices）、面の向きを使って新しい任意の画像断面を再構成するにはどうすればよいかを見ていきます。
課題は大きく3つです。
新しく作る画像の画像の面の向き、左上隅の位置、再構成面のピクセルをそれぞれ取得することです。
画像の面の向きは、新しく作る画像のRow方向の方向余弦と、Column方向の方向余弦を求めることで達成できます。これは単純に、RCS空間で、作成したい画像の面がどのような角度になるかを決めるのみです。例えば、RowとColumnのベクトルの角度を度数法で取得し、ラジアンに変換して、方向余弦で計算するなどです。ここで、RowベクトルとColumnベクトルは必ず直交するようにしなければなりません。
次に、作成される画像の左上隅位置は、比較的簡単に計算できます。例えば、画像の始点となる左上隅ピクセルの位置が、元画像シリーズの画像空間のどこにあるかを知りたいとき、次の式から計算できます。
それぞれの意味は次の通りです。
P_x, y, z：画像平面におけるボクセル ( i, j ) の座標 ( mm )。
S_x, y, z：Image Position (Patient) の 3 つの値。RCSの原点からの位置( mm )。
X_x, y, z：Image Orientation (Patient)のRow方向の方向余弦。
Y_x, y, z：Image Orientation (Patient)のColumn方向の方向余弦。
i ：画像平面に対するColumnの位置。最初の列はインデックスは 0 。
Δi ：Pixel Spacing (0028,0030) のColumnサイズ( mm )。
j ：画像平面に対するRowの位置。最初の行のインデックスは 0 。
Δj ：Pixel Spacing (0028,0030)のRowサイズ( mm )。
これに従って行列計算（内積）をすれば、指定したピクセルの位置が取得できます。
最後の再構成面のピクセルは、ソースとなる元の画像シリーズから、新しく作成する画像ピクセルの位置に対応する画素値を計算することになります。
簡単な例で確認していきます。
実践のイメージ
今、手元に横断像の頭部MRシリーズがあり、このシリーズ画像を使って、異なる画像平面をスライス（リスライス）し、新しく画像を作成したいとします。頭部MRシリーズに含まれている画像の属性は次の通りだったと仮定します。
画像サイズ：512×512 マトリクス
画像枚数：135 スライス（0～134）
Pixel Spacing：y (row) = x (column) = 0.449219
Spacing between Slices = 1.2 mm （スライス厚として利用）
例として、冠状断面（CORONAL）を作ることにしましょう。
（体の頭尾方向は、寝台の向きと同じと考えて）
まずは、作成する画像平面を定義します。必要なものは、画像の向きを示すRow方向とColumn方向の方向余弦、新しく作る画像の左上のピクセル位置および再構成されたピクセルデータです。
今回は、ちょうど横断面画像の中央からCOR画像を作ることにします。
まず、COR断面のRowとColumnの方向余弦を決めます。
COR断面は、次のような方向余弦になります。
Row方向ベクトル = [ 1, 0, 0 ]
Column方向ベクトル = [ 0, 0, -1 ]
ここで、Column方向ベクトルのcos zがマイナスになっていますが、これは、RCSの z 軸に対して、マイナス z 軸方向という意味です。
次に、新しく作る画像の左上隅の位置を求めます。
頭部MRI画像は、鼻のあたりから頭頂まで、尾頭方向に画像が並んでいます。なので、ここでは、MR画像シリーズの最後の画像（頭のてっぺんの画像。0 から数えると、134番目のスライス）から、新しく作る画像の左上隅の位置を計算します。
134番の横断面画像の左上隅の位置は、
0020,0032 Image Position (Patient): -121.6217-116.967\88.78404
であったとします。
つまり、[ x : -121.6217 mm, y : -116.967 mm, z : 88.78404 mm ]
新しく作る画像の左上隅の位置を計算する方法は、先に述べた式を用いて、次のようにできます。
{ 1 × 0.449219, 0 × 0.449219, 0, -121.6217 }
行列A { 0 × 0.449219, 1 × 0.449219, 0, -116.967 }
{ 0 × 1.2 , 0 × 1.2 , 0, 88.78404 }
{ 0 , 0 , 0, 1 }
{ 0 } # 横断面画像中、X軸方向の左隅位置
行列B { 255 } # 横断面画像中、Y軸方向の中央の位置
{ 0 }
{ 1 }
ここで、方向余弦は、元画像の方向余弦であることに注意してください。また、各要素を { } でくくっていますが、実際には先の式のような行列と思っていただければ幸いです。
行列Aと行列Bの内積で、RCS空間内の位置[mm単位]が求まります。求めると、
import numpy as np
a = [
[1*0.449219, 0*0.449219, 0,-121.6217],
[0*0.449219, 1*0.449219, 0, -116.967],
[0*1.2,0*1.2, 0, 88.78404],
[0,0,0,1]
]
b = [[0],[255],[0], [1]]
a_ = np.array(a)
b_ = np.array(b)
new_ipp = np.dot(a_, b_)
#array([[-121.6217 ],
# [ -2.416155],
# [ 88.78404 ],
# [ 1. ]])この結果から、135枚目（1から数えて）の画像の、中央の位置の一番左のピクセル（再構成するCOR画像の左上隅）は、RCS空間上で、
[ -121.6217, -2.416155, 88.78404 ]であるとわかります。
再構成画像面の方向、左上隅の位置が分かったので、最後に画像ピクセルデータを構成します。
今回は、先の図に示したように、横断面画像のピクセルから取得すればよいということがわかります。よって、今回作成されるCOR画像のスライス厚は、ピクセルサイズである0.449219となります。（細かいことをいうと、横断面画像シリーズの枚数は135枚しかないので、CORの縦方向のサイズを135×(1.2/0.44929)した画像サイズへ変換します。）
以上の情報をDICOM属性に加えて、その他、MR Image IODに必要な属性も与えれば、無事に、DICOMのCOR画像が完成です。
しかし、この説明だけでは腑に落ちない方も多いでしょう。なぜならば、曲面や、任意の角度の場合、スライス厚を変えたい場合など、この説明では触れていないからです。これらのようなピクセルデータの計算は、実際には画像のピクセル補間法を用いて行われます。この説明は、非常に画像工学的で、DICOMの範疇を簡単に超えると思われますので、ここでは触れていません。
このような処理は、3Dレンダリングが可能な画像診断用ワークステーションで可能となっています。興味深いことに、3D レンダリング方法全体は、何年も前に、ゲーム業界で作られた技術です。これらの技術が放射線技術に役立っています。
ただ、どんなに優れた技術であっても、DICOMで扱われる「属性」情報がなければ、その処理は実現できなかったことでしょう。DICOM は、単純な斜め切り画像だけでなく、血管の走行に沿って断面を作るなど、精巧な画像再構築に必要とされる属性情報を提供しています。
Stay Visionary
この記事が気に入ったらサポートをしてみませんか？