図1は、本実施形態の画像処理方法において仮想光線がボリュームデータを通過するときの仮想光線上のボクセル値(補間されたボクセル値を含む)の動き(プロファイルパターン)を表す。プロファイルパターンは仮想光線毎に定まり、仮想光線が通過する物体毎に特徴がある。ここでは、仮想光線が石灰化領域等の障害領域を通過する場合のボクセル値のプロファイルパターンを示す。図1(a)に示すように仮想光線が障害領域10の中心部分を通過する場合は、障害領域10に対応するボクセル値は大きく跳ね上がる。一方、図1(b)に示すように仮想光線が障害領域10の縁11をかすめる場合は、障害領域10に対応するボクセル値の上昇は限定的で比較的平坦である。
図2は、仮想光線が石灰化領域20の中心部分、その縁21および血流22を通過する場合のボクセル値のプロファイルパターンを示す。仮想光線が石灰化領域20の中心部分を通過する場合は、ボクセル値は高いピーク(最大値)を持ち、仮想光線が石灰化領域の縁21を通過する場合は、ボクセル値は低いピークを持つ。また、血流22に対応するボクセル値は、平坦な丘状になる。
図5は、本実施形態の画像処理方法において、MIP処理に勾配(グラディエント)を用いる場合の説明図である。本実施形態では、除外の計算をするのに、しきい値を二つ用意し、ボクセル値を値によって3つに区分する。すなわち、完全除外(十分にボクセル値の高い箇所は石灰化と見なし除去)、除外可能(ボクセル値のみでは判断がつかない領域)、非除去(通常組織)の3つに区分し、除外可能な領域のみを勾配を用いて判断する。
本実施例では仮想光線の方向ベクトルと仮想光線の通過するボクセルの勾配情報を用いることによって、石灰化領域そのものと石灰化領域の2次元的輪郭の両方を識別できる。ボクセルの勾配情報は対象ボクセルの周辺3×3×3ボクセル領域の差分を求めることによって求めることができる。 これに対して、従来の勾配による判断を行わない一定範囲の除去を行った場合は、ボクセル値の一定値以上が自動的に除去されるので、仮想光線上の最大値は必ず除去の判断のしきい値に固定されてしまう。それでは、石灰化領域の前後に存在する血流などが観察できない。本実施例では中間領域を設けて、その間では勾配を用いることによって中間領域が石灰化領域に属するのか否かの判断を行って望ましい箇所の最大値が得られる。
図17は、本実施形態の画像処理方法において、仮想光線に対応するボクセル値を除外するかどうかの判断方法の説明図である。この場合は、石灰化領域の中心部分を除外し、その縁だけを採用するための判断方法として、三次元的なボクセル値の勾配ベクトルGを用いる。
図17(a)に示すように、仮想光線が石灰化領域20の中心部分を通過する場合は、ボクセルの勾配ベクトルG30が仮想光線の進行方向と平行に近くなる。したがって、勾配ベクトルG30が仮想光線の進行方向と平行に近くなっている部分は、石灰化領域20の中心部分であると判断し、その部分を表示データから除外する。
一方、図17(b)に示すように、仮想光線が石灰化領域の縁21をかすめる場合には、ボクセルの勾配ベクトルG31が仮想光線の進行方向と垂直に近くなる。したがって、勾配ベクトルG31が仮想光線の進行方向と垂直に近くなっている部分は、石灰化領域の縁21であると判断し、その部分を表示データとして採用する。
この場合、仮想光線の進行方向ベクトルMと、ボクセル値の最大値の候補位置における勾配ベクトルGとの内積(M・G)を求め、その絶対値|M・G|を許容平行度しきい値TP(図24において詳述する)と比較することにより、表示データの除外/採用の判断を行うことができる。また、許容平行度しきい値TPは、観察者がGUIにより動的に変更することができ、表示された画像を見ながら石灰化領域の除去の程度を変更することができる。
図6は、本実施形態の画像処理方法において、除外可能の判断を行う場合の説明図である。除外可能の判断は、勾配ベクトルGと仮想光線ベクトルDを用いた評価関数fを用いて判断する。すなわち、
f(G, D) > 一定値
の条件式を満たせば除外と判断できる。例えば、
f(G, D) = (1-(1-G・D / |G|・|D|)^2 * h(|G|))
ここで、1-(1-G・D / |G|・|D|)は、勾配ベクトルGと仮想光線ベクトルDの交差角に関わる値であり、h(|G|) は、勾配ベクトルGの大きさを正規化したものである。
図7は、本実施形態の画像処理方法において、勾配(グラディエント)を用いたMIP処理(折り返し無し)のフローチャートである。これは、画面上の各ピクセルの計算の仕方であり、以下の計算を画像上の全ピクセル分行う。まず、投影位置より、投影開始点O(x,y,z)、およびサンプリング間隔ΔS(x,y,z)仮想光線方向D(x,y,z)を設定する(ステップS301)。
次に、最大値M=システム最小値、現在計算位置X(x,y,z)=Oに初期化する(ステップS302)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置、補完ボクセル値 Vを求める(ステップS303)。
次に、Vの値を判断し(ステップS304)、Vの値が除去可能な場合は、VのX位置の勾配Gを計算する(ステップS305)。そして、一定値Tと比較し(ステップS306)、f(G,D)>一定値Tでない場合(no)は、MとVを比較し(ステップS307)、M < V の場合(yes)は、M=V とする(ステップS308)。
次に、Xは終了位置まで来たかどうかを判断し(ステップS309)、Xが終了位置まで来ていない場合(no)は、X(x,y,z)にΔS(x,y,z)加算し、計算現在位置を前進させ(ステップS310)、ステップS303以降を繰り返す。一方、Xが終了位置まで来た場合(yes)は、最大値Mを計算ピクセルのピクセル値とする(ステップS311)。
図8は、本実施形態の画像処理方法において、MIP処理に勾配(グラディエント)および折り返しを用いる場合の説明図である。本実施形態の画像処理方法では、プロファイルパターンを「折り返し」と「勾配要素」を用いて改変する。それにより、最大値位置が変わるのでそれを利用する。すなわち、プロファイルパターン上のデータを置き換えることによって第2のプロファイルパターンを作成できる。このような処理を行うことによってより仮想光線上の障害物を除去するのに範囲を決定して除外するのと比較してより、なめらかな変化を画像に与えることができ、除外領域の境界がエイリアスとして描画されることを防ぐことができる。更に、画像処理のより柔軟性を与えることができる。
本実施形態では、第1のプロファイルパターンに処理を加えて第2のプロファイルパターンを作成し、第2のプロファイルパターン上での最大値を求める。従来のMIP処理では第1のプロファイルパターン上の最大値を求めるという処理を行っていたが、最大値は必ずしも観察したい箇所ではないので石灰化領域が画面に表示されると言う結果をもたらしていた。これに対して、置き換えを行うと最大値の位置がプロファイルパターン上で動くため、これを利用してより適切な値が最大値として採用される。
図9は、勾配要素によるプロファイルパターンの改変を示す説明図である。また、図10は、折り返しによるプロファイルパターンの改変を示す説明図である。プロファイルパターンの変換はプロファイルパターン上の各点を関数で変換することに等しい。すなわち、ボクセル値をV、折り返したボクセル値をflipV、仮想光線方向ベクトルをR、勾配をG = grad(V)、勾配大きさをL = |G|、仮想光線との角度をθ = arccos(G・R / |G|・|R|)として、上記を用いて、変換プロファイルパターンのボクセル値V2 = f(flipV, L, θ)を計算する。この場合、fはいかなる関数であっても良いが、f(flipV, L, θ) = flipV*LUT(L, θ) : lookup tableを例示する。
図11は、LUT関数例と石灰化領域の表示のされ方の説明図である。通常の石灰化領域表示では輪郭が不明瞭となるのに対して、本手法による石灰化領域の表示では輪郭が明瞭になる。また、勾配の方向と大きさの両方を考慮することによってより望ましい画像がえられる様子を示した。更に、LUT関数は説明のために2値のものを例示したが、LUT関数は多値であっても良い。そのようにすれば、よりなめらかな結果画像が得られる。
図12は、本実施形態の画像処理方法において、勾配(グラディエント)および折り返しを用いたMIP処理のフローチャートを示す。これは、画面上の各ピクセルの計算の仕方であり、以下の計算を画像上の全ピクセル分行う。まず、投影位置より、投影開始点O(x,y,z)、およびサンプリング間隔ΔS(x,y,z)、仮想光線方向D(x,y,z)を設定する(ステップS321)。
次に、最大値M=システム最小値、最大値M2=システム最小値、最大値位置VXを初期化し、現在計算位置X(x,y,z)=Oに初期化し、折り返ししきい値Tを設定する(ステップS322)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置、補完ボクセル値 Vを求める(ステップS323)。
次に、VのX位置の勾配Gを計算し(ステップS324)、flipV =折り返し関数(V, T)を計算し(ステップS325)、L = |G|を計算し(ステップS326)、θ= arccos(G・R / |G|・|R|)を計算し、(ステップS327)、V2 = f(flipV, L, θ)を計算する(ステップS328)。なお、プロファイルパターンはあらかじめ準備しないで動的に計算しても良い。
次に、M2とV2を比較し(ステップS329)、 M2 < V2の場合(yes)は、M=V、M2=V2とする(ステップS330)。そして、Xは終了位置まで来たかを判断し(ステップS331)、Xが終了位置まで来ていない場合(no)は、X(x,y,z)にΔS(x,y,z)加算し、計算現在位置を前進させ(ステップS332)、ステップS323以降の処理を繰り返す。一方、Xが終了位置まで来た場合(yes)は、最大値Mを計算ピクセルのピクセル値とする(ステップS333)。
図13は、本実施形態の画像処理方法において、石灰化領域の除外の判断として、プロファイルパターンの局所的な傾きの大きさを計算する処理をする。更に局所的な傾きの大きさを顕在化させるために第2のプロファイルパターンを作成する。本実施例では第2のプロファイルパターンは元となる仮想光線のプロファイルパターンの一定値以上を折り返しによって置き換えたものである。手順1では、例えば、ボクセル値の一定値(しきい値T:図22において詳述する)以上を折り返して置き換えデータである第2のプロファイルパターンを作成する。これは、ボクセル値の高い石灰化領域が描画されるのを防止し、かつ次の手順2で、ボクセル値の変化が大きい箇所を除外するのに用いるためであり、さらに、描画対象の輪郭を明らかにするためである。本処理の特徴として石灰化領域の境界を厳密に求める必要がないという点がある。また、勾配を用いる方法と比較して仮想光線上のデータのみで処理が完結するので計算が簡単で高速である。
これにより本実施形態の画像処理方法では、輪郭を描画しつつ、ボクセル値が高い石灰化領域が描画されることを防止することができる。また、最大値(しきい値T)の設定が容易となり再現性を高くすることができる。また、しきい値Tは、観察者がGUIにより動的に変更することができ、表示された画像を見ながら石灰化領域の除去の程度を変更することができる。
図3は、本実施形態の画像処理方法により取得するプロファイルパターン(1)を示す。プロファイルパターンは仮想光線上のボクセル値(補間されたボクセル値を含む)を表す。すなわち、図3(a)に示すように、血流22を通過する仮想光線には血流22に対応するボクセル値がプロファイルパターンに現れる。また、図3(b)に示すように、石灰化領域の縁21を通過する仮想光線には石灰化領域の縁21に対応するボクセル値がプロファイルパターンに現れる。また、図3(c)に示すように、石灰化領域20の中心部分を通過する仮想光線には、石灰化領域20ではなく血流22に対応するボクセル値がプロファイルパターンに現れる。
このため、表示画像に、図3(d)に示すように、石灰化領域20の中心部分を表示させずに、血流22と石灰化領域の縁(輪郭)21を表示させ、石灰化領域20の中心部分の領域には、その石灰化領域20の奥に存在する血流22を表示させることができる。これにより、石灰化領域20の大きさ(輪郭)と、石灰化領域20によって狭められた血管の血流22を同時に把握することができる。
図4は、本実施形態の画像処理方法により取得するプロファイルパターン(2)を示す。すなわち、図4(a)に示すように、血流22がほとんど存在しない場合は周辺の組織と同じボクセル値を用いる。また、図4(b)に示すように、石灰化領域の縁21を通過する仮想光線には石灰化領域の縁21に対応するボクセル値を用いる。また、図4(c)に示すように、石灰化領域20の中心部分を通過する仮想光線には、前後に血流22がないため、除外された範囲以外の最大値、例えば石灰化領域20の近辺に対応するボクセル値を用いる。これにより、観察者が、画像を見ながら石灰化領域20の大きさ(輪郭)を的確に把握することができる。
図14は、本実施形態の画像処理方法における手順2を示す。手順2では、変化の大きい箇所を描画対象に含めないようにするために、例えば、仮想光線上における二階微分の絶対値の大きさがしきい値より大きい箇所を除外する。これにより、石灰化領域の縁を描画して石灰化領域の大きさを把握できるようにするとともに、従来は描画できなかった石灰化領域の前後の血流を描画することができる。この処理を行うことによって手順1のような置き換え、折り返し処理を行っても結果的に図13のように折り返しを行ったしきい値が必ず最大値として得られるにとどまってしまうと言う問題を克服できる。
図15は、本実施形態の画像処理方法において、仮想光線に対応するボクセル値を除外するかどうかの判断方法(1)の説明図である。また、同図(a)は、仮想光線が石灰化領域の縁を通過する場合のボクセル値を示す。この場合は、例えば、ボクセル値の上下変動の大きさ、あるいは勾配の大きさを測り、例えば、許容変化しきい値TD(図23において詳述する)と比較することにより、仮想光線が通過した領域が石灰化領域の縁であることを検出し、そのボクセルは表示すべきデータ(採用データ)と判断する。
一方、図15(b)は、仮想光線が石灰化領域の中央部分を通過する場合のプロファイルパターンを示す。この場合は、ボクセル値の上下変動の大きさ、あるいは勾配の大きさを測ることにより、仮想光線が通過した領域が石灰化領域の中央部分であることを検出し、プロファイルパターン上の一定範囲をボクセルを表示すべきデータから除外する。これにより、石灰化領域の中央部分を描画せず、その縁すなわち輪郭だけを描画し、さらに石灰化領域の中央部分の前後に存在する血流を描画することができる。また、許容変化しきい値TDは、観察者がGUIにより動的に変更することができ、表示された画像を見ながら石灰化領域の除去の程度を変更することができる。
図16は、本実施形態の画像処理方法において、仮想光線に対応するボクセル値を除外するかどうかの判断方法(2)の説明図であり、仮想光線が石灰化領域の縁を通過する場合のボクセル値を示す。この場合には、折り返し点(スパイク)前後の勾配の大きさを参照し、折り返し点前後の勾配の片側だけが大きいので、この部分は石灰化領域の縁であると判断し、表示データとして採用する。簡単化のために、折り返し点前後のボクセル値を用いて二階差分値を取得してもよい。
p0 − 2*p1 + p2 > threshold
図18は、本実施形態の画像処理方法において、しきい値を用いた折り返し以外の変換方法を用いて石灰化領域の縁を検出する場合の説明図を示す。すなわち、ボクセル値をしきい値で折り返すのではなく、ボクセル値を所定のしきい値で丸め込み、丸め込まれたグラフの勾配(グラディエント)を計算し、その勾配の変化により表示データの除外/採用を判断する。
すなわち、図18(a)に示すように、所定のしきい値に丸め込まれたボクセル値の勾配が小さい場合、その領域は石灰化領域の縁であると判断して表示データとして採用する。一方、図18(b)に示すように、所定のしきい値に丸め込まれたボクセル値の勾配が大きい場合は、その領域は石灰化領域の中央部分であると判断して表示データから除外する。この方法によれば、しきい値での折り返し処理を行わずに、表示データの除外/採用を判断することができる。
図19は、本実施形態の画像処理方法において、石灰化領域の縁の表示データとして、折り返す前の元データを用いる場合の手順(3)を示す。本実施形態において、折り返しの基準となる最大値は、観察者がGUIにより動的に設定することができるが、仮想光線の通過位置によっては、図19(b)に示すように、石灰化領域の縁に対応するボクセル値も折り返しの対象になることがある。この場合は、図19(a)に示すように、折り返し処理を行う前の石灰化領域の縁に対応するボクセル値(元データ)を用いることにより、石灰化領域の縁を明瞭に描画することができる。
図20は、本実施形態の画像処理方法における望ましい実施例を示す。図20(a)は、仮想光線上に石灰化領域の中心部分、石灰化領域の縁および血流が存在する場合であり、最大値における折り返し処理を行ってボクセル値の変化の大きさを検出することにより、仮想光線が通過した領域が石灰化領域の中心部分か、あるいは石灰化領域の縁かを判断し、石灰化領域の中心部分を表示せずに石灰化領域の縁だけを表示することができる。
また、図20(b)は、仮想光線上に石灰化領域の中心部分および血流が存在する場合であり、最大値における折り返し処理を行ってボクセル値の変化の大きさを検出することにより、石灰化領域の中心部分を判断し、石灰化領域の中心部分を表示せずにその奥の血流を表示することができる。
図21は、通常のMIP画像(a)と本実施形態の画像処理方法によって作成したMIP画像(b)の例を示す。図21(a)に示す通常のMIP画像では、石灰化の存在によって石灰化の前後の血流がどのようになっているかわからないが、図21(b)に示す本実施形態のMIP画像では、石灰化は輪郭のみが表示され、前後に血流が残っているか否かが把握できるようになっている。
図22は、本実施形態の画像処理方法において、処理の全体像を示すフローチャートを示す。画像の各ピクセル値を求めるには、まず、図13に示したように、注目する組織、例えば血液のボクセル値よりいくぶん大きいしきい値Tを決定する(ステップS11)。
次に、仮想光線を投射し(ステップS12)、仮想光線の上のボクセル値を配列A1(元の(第1の)プロファイルパターン)として取得する(ステップS13)。次に、配列A1のしきい値T以上の値をしきい値Tで折り返した配列A2(第2のプロファイルパターン)を作成し(ステップS14)、配列A2上で一部データ、例えば、石灰化領域の中心部分に相当する折り返えされたデータを除外する(ステップS15)。
次に、配列A2上で最大値となる値M1を求め(ステップS16)、配列A1上で値M1に対応する値M2(図19(a)の元データ)を求める(ステップS17)。そして、値M2をこの仮想光線に対するピクセル値とする(ステップS18)。
本実施形態では、ステップS14において、バッファを作成して置き換え配列A2を求め、続くステップS15において、配列A2上の一部データを除外しているが、GUIによる観察者からの指示に応答し、仮想光線の進行にあわせて動的に置き換え配列A2を求め、一部データを除外することもできる。これにより、観察者が、注目する対象や観察する方向を変えながら対象を詳細に観察することができる。
図23は、本実施形態の画像処理方法において、ボクセルデータの変化を用いてデータ除外の判断を行う場合のフローチャートを示す。まず、図15に示したように、ボクセルデータの変化(グラフの傾斜)の大きさを判断するための許容変化しきい値TDを決定する(ステップS21)。
次に、配列A2上を走査開始し(ステップS22)、配列A2上の位置Pを用いて走査する(ステップS23)。次に、配列A2上の位置P−1の値V0、配列A2上の位置Pの値V1、配列A2上の位置P+1の値V2をそれぞれ取得する(ステップS24)。
次に、変化D0 = |V1-V0|、変化D1 = |V2-V1|をそれぞれ取得し(ステップS25)、変化D0および変化D1と許容変化しきい値TDの大きさを比較する(ステップS26)。そして、変化D0および変化D1が許容変化しきい値TDより大きい場合(yes)は、位置P上の点を除外する(ステップS27)。
一方、変化D0および変化D1が許容変化しきい値TDより大きくない場合(no)は、位置Pが配列の末尾まで来たかどうかを判断し(ステップS28)、位置Pが配列の末尾でない場合(no)は、位置Pを+1前進させ(ステップS30)、ステップS24以降の処理を繰り返す。一方、ステップS28において、位置Pが配列の末尾まで来た場合(yes)は、配列A2上の走査を終了する(ステップS29)。
これにより、石灰化領域の中央部分を描画せず、その縁すなわち輪郭だけを描画し、さらに石灰化領域の中央部分の前後に存在する血流を描画することができる。また、許容変化しきい値TDは、観察者がGUIにより動的に変更することができ、表示された画像を見ながら石灰化領域の除去の程度を変更することができる。
図24は、本実施形態の画像処理方法において、勾配を用いてデータ除外を判断する場合のフローチャートを示す。この場合は、図17に示したように、まず、仮想光線の進行方向ベクトルMを取得し(ステップS41)、許容平行度しきい値TPを決定する(ステップS42)。
次に、配列A2上を走査開始し(ステップS43)、配列A2上の位置Pを用いて走査する(ステップS44)。また、配列A2上の位置Pに対応するボリュームデータ上のボクセルより勾配ベクトルGを求める(ステップS45)。
次に、仮想光線の進行方向ベクトルM とボリュームデータの勾配ベクトルGとの内積の絶対値と許容平行度しきい値TPを比較し(ステップS46)、進行方向ベクトルM と勾配ベクトルGの内積の絶対値が許容平行度しきい値TPより大きい場合(yes)は、位置P上の点を除外する(ステップS47)。
一方、進行方向ベクトルM と勾配ベクトルGの内積の絶対値が許容平行度しきい値TPより大きくない場合(no)は、位置Pが配列の末尾まで来たかどうかを判断し(ステップS48)、位置Pが配列の末尾でない場合(no)は、位置Pを+1前進させ(ステップS50)、ステップS45以降の処理を繰り返す。一方、位置Pが配列の末尾まで来た場合(yes)は、配列A2上の走査を終了する(ステップS49)。
これにより、勾配ベクトルGが仮想光線の進行方向ベクトルMと平行に近くなっている部分は、石灰化領域の中心部分であると判断して表示データから除外し、一方、勾配ベクトルGが仮想光線の進行方向ベクトルMと垂直に近くなっている部分は、石灰化領域の縁であると判断してその部分を表示データとして採用することができる。また、許容平行度しきい値TPは、観察者がGUIにより動的に変更することができ、表示された画像を見ながら石灰化領域の除去の程度を変更することができる。
さて、レイキャスト法であっては各ボクセルに不透明度が設定可能である。そこで、仮想光線上の範囲を除外するのではなく、仮想光線上のデータに不透明度を割り当てることができる。
これによれば、例えば、ボクセル値の勾配の急な場所を硬い場所、ボクセル値の勾配の緩やか場所を柔らかい場所、と言った表現が行えるようになる。また、これによって勾配の急な場所を選択して表示しないことが可能になる。さらに、勾配の急な場所(石灰化の境界面)を検出することができるようになるので、マスク処理などを伴わずに石灰化部分の除去された画像を提示することができる。
図25は、障害領域の透明度を変更し、傷害領域の中心部分を表示させないようにする場合(レイキャスト法に対する応用1)の説明図を示す。描画対象のボクセル値の透明度を変更することによっても、障害となる領域を除外することができる。すなわち、描画対象の除外の度合いαを計算し、除外の度合いαを透明度と関連付ける。
この場合は、折り返したデータの仮想光線の方向ベクトルVと勾配ベクトルGの内積を用いて、除外の度合いα=|G*V|を計算し、方向ベクトルVと勾配ベクトルGが平行で除外の度合いαが高い部分は、石灰化領域の中心部分なので、そのボクセル値の透明度を高くして石灰化領域を表示しないようにする。これによれば、折り返し計算等の処理を行うことなく、石灰化領域の中心部分を表示せずに、石灰化領域の縁および石灰化領域の前後の血流を表示することが可能である。
図26は、ボクセル値の勾配(グラディエント)を用いて透明度を変更する場合(レイキャスト法に対する応用)のフローチャートを示す。この処理は画面上の各ピクセルの計算であり、以下の計算を画像上の全ピクセル分行う。まず、投影位置より、投影開始点O(x,y,z)、およびサンプリング間隔ΔS(x,y,z)を設定する(ステップS101)。
次に、反射光Eを「0」、残存光Iを「1」、現在計算位置X(x,y,z)を「O」に初期化する(ステップS102)。そして、現在計算位置X(x,y,z)の周辺のボクセルデータより現在計算位置X、補完ボクセル値Vを求める(ステップS103)。また、補完ボクセル値Vに対応する不透明度αを得る(ステップS104)。この場合、α1=f(V),α2=1-G*(X-O),α=α1* (1-α2)により、位置Xの除外度合いを計算する(ステップS112)。
次に、補完ボクセル値Vに対応するカラー値Cを得る(ステップS105)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置および勾配Gを求める。また、光線方向X-Oと勾配Gよりシェーディング係数βを求める(ステップS106)。
次に、X(x,y,z)位置の減衰光D(D=I*α)及び部分反射光F(F=β*D*C)を計算する(ステップS107)。そして、反射光Eと残存光Iを更新し(I=I-D,E=E+F)、現在計算位置を進行させる(X=X+ΔS)(ステップS108)。
次に、現在計算位置Xが終了位置まで来たかどうか、または残存光Iが「0」になったかどうかを判断し(ステップS109)、現在計算位置Xが終了位置でなく、残存光Iが「0」でない場合(no)は、X(x,y,z)にΔS(x,y,z)を加算し、現在計算位置を前進させ(ステップS110)、ステップS103以降の処理を繰り返す。
一方、現在計算位置Xが終了位置まで来たか、または残存光Iが「0」になった場合(yes)は、反射光Eを計算ピクセルのピクセル値として計算を終了する(ステップS111)。このように、従来はレイキャスト法であっては反射光量を計算するためにボクセルの勾配と光源方向を掛け合わせたシェーディング係数を用いたのだが、勾配の物理的意味を離れて勾配と仮想光線方向を掛け合わせることによって石灰化領域とその2次元的境界面が検出可能である。更に検出された情報によってボクセルの不透明度を変更することによって、石灰化領域の中心部分を表示せずに、石灰化領域の縁および石灰化領域の前後の血流を表示することが可能である。
図27は、透明度を変更する場合(レイキャスト法に対する応用2)の説明図を示す。この場合は、折り返したボクセル値の変化から除外の度合い(α)を求める。すなわち、除去の判定部分において、折り返したデータに対して、例えば、ボクセル値が急変して折り返している箇所をα=1とし、ボクセル値が急変して折り返していない箇所をα=0.5とし、それ以外の箇所をα=0とする。そして、α=1の領域を障害領域として表示せず、α=0.5の領域を障害領域の輪郭部分として表示する。このように、折り返したボクセル値の変化から除外の度合い(α)を求めることにより、石灰化領域の中心部分を表示せずに、石灰化領域の縁および石灰化領域の前後の血流を表示することが可能である。
図28は、ボクセル値の変化を用いて透明度を変更する場合(レイキャスト法に対する応用)のフローチャートを示す。この処理は画面上の各ピクセルの計算であり、以下の計算を画像上の全ピクセル分行う。まず、投影位置より、投影開始点O(x,y,z)、およびサンプリング間隔ΔS(x,y,z)を設定する(ステップS121)。
次に、反射光Eを「0」、残存光Iを「1」、現在計算位置X(x,y,z)を「O」に初期化する(ステップS122)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置、補完ボクセル値Vを求める(ステップS123)。また、補完ボクセル値Vと「除外」を考慮した不透明度αを得る(ステップS124)。この場合、α1 = f( V),α2 = g(X),α = α1 * (1-α2)により位置Xの除外度合いを計算する(ステップS132)。
次に、補完ボクセル値Vに対応するカラー値Cを得る(ステップS125)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置、および勾配Gを求める。また、光線方向X-Oと勾配Gよりシェーディング係数βを求める(ステップS126)。
次に、X(x,y,z)位置の減衰光D(D=I*α)及び部分反射光F (F=β*D*C)を計算する(ステップS127)。そして、反射光Eと残存光Iを更新し(I=I-D,E=E+F)、現在計算位置を進行させる(X=X+ΔS)(ステップS128)。
次に、Xが終了位置まで来たかどうか、または残存光Iが「0」になったかどうかを判断し(ステップS129)、Xが終了位置でなく、残存光Iが「0」でない場合(no)は、X(x,y,z)にΔS(x,y,z)を加算して現在計算位置を前進させ(ステップS130)、ステップS123以降の処理を繰り返す。
一方、Xが終了位置まで来たか、または残存光Iが「0」になった場合(yes)は、反射光Eを計算ピクセルのピクセル値として計算を終了する(ステップS131)。このように、ボクセル値の変化を用いて透明度を変更することにより、石灰化領域の中心部分を表示せずに、石灰化領域の縁および石灰化領域の前後の血流を表示することが可能である。
図29は、勾配処理と折り返し処理を行うMIPフローチャートを示す。この処理は画面上の各ピクセルの計算であり、以下の計算を画像上の全ピクセル分行う。まず、投影位置より、投影開始点O(x,y,z)、およびサンプリング間隔ΔS(x,y,z)を設定する(ステップS141)。
次に、最大値Mをシステム最小値に、現在計算位置X(x,y,z)を「O」に初期化し(ステップS142)、X(x,y,z)位置周辺のボクセルデータよりX位置、補完ボクセル値Vを求める(ステップS143)。次に、X位置周辺のボクセルデータVS1を取得し(ステップS144)、折り返しデータVS2をVS1から計算する(ステップS145)。
次に、折り返しデータVS2のX位置の勾配Gを計算し(ステップS146)、光線方向X-Oと勾配Gの内積Iを求める(ステップS147)。そして、一定値Tに対して -T<I<T の条件を満たすかどうかを判断し(ステップS148)、その条件を満たす場合(yes)は、最大値M と補完ボクセル値Vを比較し(ステップS149)、最大値M が補完ボクセル値Vより小さい場合(yes)は、最大値Mを補完ボクセル値Vとする(ステップS150)。
次に、Xが終了位置まで来たかどうかを判断し(ステップS151)、Xが終了位置まで来ていない場合(no)は、X(x,y,z)にΔS(x,y,z)を加算して現在計算位置を前進させ(ステップS152)、ステップS143以降の処理を繰り返す。また、Xが終了位置まで来た場合(yes)は、最大値Mを計算ピクセルのピクセル値として終了する(ステップS153)。このように、勾配処理と折り返し処理を組み合わせることにより、石灰化領域の中心部分を表示せずに、石灰化領域の縁および石灰化領域の前後の血流を表示することが可能である。
また、本実施形態の画像処理方法は、GPU(Graphic Processing Unit)により行うことができる。GPUは、汎用のCPUと比較して特に画像処理に特化した設計がなされている演算処理装置で、通常CPUとは別個にコンピュータに搭載される。
また、本実施形態の画像処理方法は、ボリュームレンダリングの計算を所定の画像領域、ボリュームの領域等で分割し、後で重ね合わせることができるので、並列処理やネットワーク分散処理、専用プロッセッサ、或いはそれらの複合により行うことができる。
以上、本発明の実施形態について説明したが、本発明は上記の実施形態に限定されるものではなく、以下の態様も可能である。
(その1)
例えば、観察対象付近のコントラストの強い塊をAとし、投影方法およびその方向に応じて定まるベクトルをBとする。ここで、Bは平行投影の場合は見る方向、VE(Virtual Endoscope)の場合は放射状ベクトル、あるいは見る方向と直行するベクトルである。
この場合に、レイキャストによるレンダリング処理において、Aが観察対象を観察する妨げになることを防ぐ目的で、レイキャストの各ステップにおいて動的にAの外殻部が特に選択的に強調される機構を持ち、強調された結果とBに応じて、(1)当該ステップをスキップする、または(2)元のデータ値を他の値と置換して処理をする、または(3)光の減衰をシミュレートするレンダリング法の場合に、減衰量を変化させる、のいずれかの処理を行う。
(その2)
観察対象付近のコントラストの強い塊をAとし、観察したい対象物体をBとし、投影方法およびその方向に応じて定まるベクトルをCとする。ここで、Cは平行投影なら見る方向、VEなら放射状ベクトル、または見る方向と直行するベクトルである。
この場合に、レイキャストによるレンダリング処理において、AがBの観察の妨げになっている状況において、(1)Aと重ならない部分ではBの画像を損ねることなく、(2)Aと重なる部分ではAの特徴を簡潔に示すAの画像の一部を残して除去あるいは半透明にして、(3)Bの全体像をAに妨げられず観察できるようにするという目的で、レイキャストの各ステップにおいて動的にAの外殻部が特に選択的に強調される機構を持ち、強調された結果とCに応じて、(1)当該ステップをスキップする、または(2)元のデータ値を他の値と置換して処理をする、または(3)光の減衰をシミュレートするレンダリング法の場合、減衰量を変化させる、のいずれかの処理を行う。
(その3)
例えば大腸のエアー像(半透明で輪郭だけ見え易くしたVR像)にも応用できる。エアー像の場合、「A=B」であり「Aの特徴を簡潔に示す部分=Bの観察したい部分」となる。
この場合、コントラストの強い塊をAとし、観察したい対象物体をB(BがAと同一である場合もある)とし、投影方法およびその方向に応じて定まるベクトルをCとする。Cは平行投影なら見る方向、VEなら放射状ベクトル、あるいは見る方向と直行するベクトルである。
レイキャストによるレンダリング処理において、AがBの観察の妨げになっている状況において、(1)AとBが重ならない部分ではBの観察に対して悪影響を及ぼすような変化をBの画像に与えること無く、(2)AとBが重なる部分ではAの特徴を簡潔に示すAの画像の一部を残して除去あるいは半透明にして、(3)Bの全体像をAに妨げられず観察できるようにするという目的で、レイキャストの各ステップにおいて動的にAの外殻部が特に選択的に強調される機構を持ち、強調された結果とCに応じて、(1)当該ステップをスキップする、または(2)元のデータ値を他の値と置換して処理をする、または(3)光の減衰をシミュレートするレンダリング法の場合、減衰量を変化させる、のいずれかの処理を行う。
(その4)
コントラストの強い塊をAとし、観察したい対象物体をB(BがAと同一である場合もある)とし、投影方法およびその方向に応じて定まるベクトルをCとする。Cは平行投影なら見る方向、VEなら放射状ベクトル、あるいは見る方向と直行するベクトルである。
この場合、レイキャストによる画像投影方法において、AがBの観察の妨げになっている状況において、(1)AとBが重ならない部分ではBの観察に対して悪影響を及ぼすような変化をBの画像に与えること無く、(2)AとBが重なる部分ではAの特徴を簡潔に示すAの画像の一部を残して除去あるいは半透明にして、(3)Bの全体像をAに妨げられず観察できるようにするという目的で、レイキャストの各ステップにおいて動的にAの外殻部が特に選択的に強調される機構を持ち、強調された結果とCが、あらかじめ定められた条件に適合する場合は、(1)当該ステップをスキップ、または(2)元のデータ値を他の値と置換してから通常処理、または(3)通常処理の前段階として強調結果とCによる変換を付加した処理、または(4)光の減衰をシミュレートするレンダリング法の場合、減衰量を変化させる、のいずれかの処理を通常処理の代用として行う。
(その5)
以上は、Aの輪郭を残すことを前提としているが、Aの輪郭を残さない態様も可能である。すなわち、コントラストの強い塊をAとし、観察したい対象物体をB(BがAと同一である場合もある)とし、投影方法およびその方向に応じて定まるベクトルをCとする。ここで、Cは平行投影なら見る方向、VEなら放射状ベクトル、あるいは見る方向と直行するベクトルである。
この場合、レイキャストによる画像投影方法において、AがBの観察の妨げになっている状況において、(1)AとBが重ならない部分ではBの観察に対して悪影響を及ぼすような変化をBの画像に与えること無く、(2)AとBが重なる部分では、(2.1)Aの特徴を簡潔に示すAの投影像の一部以外、あるいは(2.2)Aの投影像の全てを除去あるいは半透明にして、(3)Bの全体像をAに妨げられず観察できるようにするという目的で、レイキャストの各ステップにおいて動的にAの外殻部が特に選択的に強調される機構を持ち、レイキャストの各ステップが対象とする位置とその近傍における「元データ値/強調された結果/C」が、あらかじめ定められた条件に適合する場合は、(1)当該ステップをスキップ、(2)元のデータ値を他の値と置換してから通常処理、(3)通常処理の前段階として強調結果とCによる変換を付加した処理、(4)光の減衰をシミュレートするレンダリング法の場合、減衰量を変化させる、のいずれかの処理を通常処理の代用として行う。
(その6)
本発明の画像処理方法は、(A)ボリュームレンダリングにおいて、(B)仮想光線の上の各サンプル点において対象サンプル点の値を用いるか否か(除外)判断するステップを保有する。(C)目的は障害領域の輪郭のみを描画することである。また、(D)除外の判断にはボリュームデータから作成したデータを用いるとより望ましい(これは必須要素ではない)。
本発明の特徴として、従来はボリュームを改変したり、マスクボリュームを作成することによって障害領域を除去してきたが、本発明においては仮想光線投射時に動的に障害領域を計算するので投影方向を加味した画像(視点方向から見ての輪郭の描出)が作成できる。その為に、視線方向を動かしたときに輪郭表現が視線に追従し望ましい画像が容易に得られる。
本発明の画像処理方法は、(A)ボリュームレンダリング全般で使用できるが、特にMIP法、レイキャスト法において有効に機能する。ボリュームレンダリングでは、MIP法に代表される(A1)仮想光線上のデータから一点のサンプル点を選択して画素を決定するMIP法、MinIP法等の手法と、レイキャスト法に代表される(A2)仮想光線上の複数のサンプル点を用いて画素を決定するレイキャスト法、RaySUM法、平均値法等の手法がある。
本発明はいずれの方法であっても適用が可能であり、A1であってはBのステップを経て「仮想光線上のデータ」からBで選択されたサンプル点を除外することで実装できる(MIP法での実装例)。一方、A2であってはBのステップを経て「仮想光線上の複数の点」からBで選択されたサンプル点を除外することで実装できる。
また、レイキャスト法等のA2の方法であっては、上記のようにA,B,Cをそのまま実装することもできるが、必ずしもサンプル点を除外する必要はなく、画素決定に対する寄与度を下げることもできる。
例えば、レイキャスト法では仮想光線が通過するサンプル点に対してボクセル値を元に不透明度αを計算する。ここで、Bで選択されたサンプル点に対して不透明度αの値を操作することができる。
(B)仮想光線の上の各サンプル点において対象サンプル点の値を用いるか否か(除外)判断するステップには、上述した(B1)勾配を用いた方法(注:折り返したデータでは勾配は全サンプル点において計算しても、一部だけ(しきい値近傍など)計算しても良い)、および(B2)変化を用いる方法のほかに、仮想光線上のデータを(a)周波数解析して高周波成分を除外する。(b)データの分散を計算して、分散の一定値以上のサンプル点を除外する。(c)雑音除去フィルターの使用等の方法によって除外の判断ができる。また、除去は「除く」「除かない」の二値ではなく「除去の度合い」として多値で求めることもできる。
(C)本発明の目的は、障害領域の輪郭のみを描画することであるが、そのほか、輪郭に限らず、データの二次元的特徴を表現することができる。また、輪郭に関わりなく動的に除外を計算することによって医療画像における新しい画像を作成することも可能である。
(D)また、除外の判断にはボリュームデータから作成したデータを用いるとより望ましい(これは必須要素ではない)。この場合、(D1)折り返しデータを用いる。(D2)しきい値以上を切ったデータを用いる(勾配を用いるの場合)。(D3)勾配を用いる場合、折り返しを全く行わず、(1)ボクセル値がしきい値よりだいぶ大きい時は除去、(2)ボクセル値がしきい値近傍の時は勾配を計算して除去を判断、(3)ボクセル値がしきい値よりだいぶ小さい時は非除去としても同等の効果が得られる。
(D4)マスクデータを用いる。すなわち、マスクデータに対して勾配を計算しても同等の効果が得られる。(D5)ボリュームデータから作成したデータは、(a)仮想光線投射時に計算してもいいし、(b)あらかじめ第二のボリュームデータとして作成しておいても良い。
上記各実施形態では、石灰化領域の識別を行ったが障害領域であれば何でも良い。例えばステント、医療クリップ、医療コイルなどの人体に挿入した医療器具の観察には特に有効である。また、血流や臓器の領域でもよい。また、造影剤によって信号を強調された造影領域であっても良い。造影領域は低濃度造影による造影領域を含む。また、障害領域は高信号領域に限らず、例えば気泡のような低信号領域の2次元的輪郭を表示するのにも用いることができる。また、脂肪のような中間的な信号強度をもつ領域の2次元的輪郭を表示するのにも用いることができる。
上記各実施形態では、障害領域は塊状領域であったが障害領域であればいかなる毛上であっても良い。例えば、図30のように血流が複雑に交錯する画像であれば血流の2次元的輪郭を表示することによって血流を正確に認識することが可能となる。また、腸の造影であれば、腸の2次元的輪郭を表示することによって複雑に折りたたまれている腸壁の位置関係や形状がよくわかる。
上記各実施形態では、本発明による画像を単独で表示したが本発明による画像はその他の画像と組み合わせて表示しても良い。例えば、本発明の画像と従来画像を並べて表示したり、重ね合わせて表示したりしても良い。この場合、同じ画角、平行投影法、透視投影法、円筒投影法による画像、仮想光線上のデータは一括して取得しても良いし、計算の必要に応じて逐次取得してもかまわない。
上記各実施形態では、画像全体を同一の計算手法で作成していたが、画像の一部にのみ本手法による画像処理を適用しても良い。例えば、マウスなどのポインティングデバイスの周辺に対してのみ本手法による画像処理を行っても良い。このようにすれば、従来手法による画像処理と本手法による画像処理を比較しつつ観察することが容易になる。
上記各実施形態では、画像全体を同一の計算手法で作成していたが、画像の一部にのみ本手法による画像処理を適用しても良い。特に本手法であっては障害領域の存在しない箇所では従来手法と同一の画像を得ることができるので、あらかじめ、従来の方法で画像を作成し、特に障害領域が含まれている領域を自動的に求めることによって、処理速度を上げることができる。特にMIP法を用いた実施例に関しては、仮想光線上の最大値が一定値以上の箇所でのみ本手法を適用することができる。
上記各実施形態では、勾配情報を対象ボクセルの周辺3×3×3ボクセル領域の差分を求めることによって求めたが、勾配情報であれば求め方は問わない。例えば対象ボクセルではなく、仮想光線の通過する点の値を周辺ボクセルの補間によって求める場合は対象ボクセルの周辺3×3×3ボクセルもそれぞれ補間によって求めて良い。また、3×3×3ボクセルでなくとも2×2×2領域であっても良い。また、仮想光線方向にボクセルを正規化した上で勾配を求めて良い。
上記各実施形態では、仮想光線上のボクセルの画像に対する寄与度が直感的に判断することが難しいので、ボリュームデータに含まれるボクセルの画像に対する寄与度を別個に可視化した画像を表示することができる。例えば、仮想光線上のボクセルの寄与度を表現したグラフを表示することができる。また、複数の仮想光線上に関して一括して表示することもできる。例えば、画像上の線に対応する仮想光線群で構成される面を用いたCPR(Curved Multi Planer Reconstruction)画像を作成し、CPR画像上のそれぞれの点に関わる仮想光線上のボクセルの寄与度を重ね合わせてマップすることができる。
上記各実施形態では、仮想光線上のボクセルの画像計算において用いられたか除外されたかを直感的に判断することが難しいので、ボリュームデータに含まれるボクセルが画像作成において用いられたか除外されたかを別個に可視化した画像を表示することができる。例えば、仮想光線上のボクセルが画像作成において用いられたか除外されたかグラフを表示することができる。また、複数の仮想光線上に関して一括して表示することもできる。例えば、画像上の線に対応する仮想光線群で構成される面を用いたCPR(Curved Multi Planer Reconstruction)画像を作成し、CPR画像上のそれぞれの点に関わる仮想光線上のボクセルが画像作成において用いられたか除外されたかを重ね合わせてマップすることができる。
上記各実施形態では、仮想光線上のプロファイルパターンがどのように置き換えられたかを直感的に判断することが難しいので、仮想光線上の置き換え後のプロファイルパターンを表示しても良い。例えば、ポインティングデバイスの指し示す点に対応する仮想光線上の置き換え後プロファイルパターンを画像に重ね合わせ表示することや別ウィンドウ上に表示することができる。また、元のプロファイルパターンも重ねて表示することもできる。
上記各実施形態では、特に仮想光線上の最大値などの単一の点のみを用いた画像処理方法であっては単一の点の取得された3次元上の位置を直感的に判断することが難しいので、単一の点の取得された3次元上の位置を可視化することができる。例えば、単一の点の取得された3次元上の位置を深度情報としてマップした画像が作成できる。また、マップ画像を本実施方法で作成した画像や従来の計算法を用いて作成した画像と重ね合わせたり並べて表示したりできる。
上記各実施形態では、単一の演算装置を用いて計算したが複数の演算装置を用いて並列処理しても良い。例えば、本実施方法であっては仮想光線毎に独立して計算が行えるので、仮想光線毎に並列して計算が行える。
上記各実施形態では、単一の中央演算処理装置を用いて画像処理を行ったが、画像処理はプログラマブルシェーダを装備したGPU(Graphics Processing Unit)を用いても良い。また、その他の演算処理装置を用いても良い。
上記各実施形態では、画像処理に用いるパラメータをあらかじめ定めてあったが、ユーザーの操作にあわせて動的に変更されるようにしても良い。例えば、折り返し処理のしきい値を画面上に設定したスライダを操作することによって変化させることができる。
上記各実施形態では、仮想光線の方向ベクトルを用いたがそれ以外の方向ベクトルであってもかまわない。例えば仮想光線の方向と斜めに交差する方向ベクトルを用いることによってMIP法であっても影に相当する表現ができるようになる。また、例えば血管の中心線の方向ベクトルを用いることで血管の手前の領域と血管の奥の領域とで異なった表現ができるようになる。
上記各実施形態では、画素値の決定に当たって仮想光線上のデータの1以上の点の値を用い、1以上の点の位置関係が相互交換可能なボリュームレンダリング法として、MIP法を用いたがそれ以外の方法であってもかまわない。例えば、最小値を用いるMinIP法(Minimum Intensity Projection)を用いても良い。また、例えば、1以上の点の平均値を用いる平均値法を用いても良い。また、例えば、1以上の点の合計値を用いるRay Sum法を用いても良い。また、例えば、仮想光線上の上位10の値を取得しそれらの平均値を用いるTop 10 MIP法であっても良い。
上記各実施形態では、ボリューム内の各位置の勾配を用いたが、ボリューム内の仮想光線の通過点より動かした位置の勾配を用いてもかまわない。このようにすることによって2次元的輪郭の描出される位置が偏差し、2次元的輪郭の存在を観察者に示唆しつつ境界領域の対してより精密な観察が行えるようになる。
上記各実施形態では、ボリューム内の各位置の勾配を用いたが、勾配情報は画素の決定委に用いるのと他のボリュームを用いても良い。例えば、時系列的に異なる他のボリュームデータを用いることによって、2次元的輪郭の移動を描画できるようになる。また、例えば、事前に作成したマスクボリュームを用いると、輪郭抽出アルゴリズム或いは観察者によって作成された輪郭情報を描画できるようになる。
上記各実施形態では、静止画を作成したが、動画を作成しても良い、また、観察者の操作にあわせて表示される動画を動的に作成しても良い。例えば、観察対象の周囲を視点が旋回する動画を作成すれば、観察対象の2次元的輪郭をより精密に観察できるようになる。また、例えば、画像処理に関わるパラメータを変化させる動画を用いた場合も2次元的輪郭をより精密に観察できるようになる。
上記各実施形態では、観察対象の2次元的輪郭を表示したが、表示するのは2次元的輪郭にとどまらない。例えば、仮想光線の方向ベクトルと勾配ベクトルとの角度が小さい箇所を除去するのではなく角度が大きい箇所を除去すれば観察対象の中央部分が強調される。また、例えば、仮想光線の方向ベクトルと勾配ベクトルとの外積ベクトルの方向によって識別すれば2次元的輪郭のうち特定の方向を向いているものを表示することができる。
上記各実施形態では、レイキャスト法を用いた場合に不透明度を設定することをしたが、不透明度の限らず画像に対する寄与度を設定できる。例えば、平均値法において寄与度を設定し、平均値に代わって加重平均を用いて画像が作成できる。
上記各実施形態では、石灰化領域の2次元的輪郭を表示したが、表示するのは石灰化領域にとどまらない。例えばステント、医療クリップ、医療コイルなどの人体に挿入した医療器具の観察には特に有効である。また、高信号領域に限らず、例えば気泡のような低信号領域の2次元的輪郭を表示するのにも用いることができる。
上記各実施形態では、CT装置より得られたボリュームデータを用いたが、ボリュームデータであれば作成手段は問わない。例えば、MRI装置(Magnetic Resonance Imaging),PET装置(Positron Emission Tomography)より得られたボリュームデータを用いることができる。また、ボリュームデータにフィルタなどを作用させて改変したボリュームデータや、複数のボリュームデータを合成したボリュームデータを用いても良い。
上記各実施形態では、第1のベクトル情報と第2のベクトル情報の角度を計算したが、ここで角度は負の値であっても直角より大きくても良い。このようにすれば手間側の壁のみを表示するようなことができるようになる。