JP2006238938A - 画像処理方法 - Google Patents

画像処理方法 Download PDF

Info

Publication number
JP2006238938A
JP2006238938A JP2005054863A JP2005054863A JP2006238938A JP 2006238938 A JP2006238938 A JP 2006238938A JP 2005054863 A JP2005054863 A JP 2005054863A JP 2005054863 A JP2005054863 A JP 2005054863A JP 2006238938 A JP2006238938 A JP 2006238938A
Authority
JP
Japan
Prior art keywords
image processing
processing method
value
image
virtual ray
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2005054863A
Other languages
English (en)
Other versions
JP4212564B2 (ja
Inventor
Kazuhiko Matsumoto
和彦 松本
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ziosoft Inc
Original Assignee
Ziosoft Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ziosoft Inc filed Critical Ziosoft Inc
Priority to JP2005054863A priority Critical patent/JP4212564B2/ja
Priority to US11/346,058 priority patent/US7853057B2/en
Priority to GB0602433A priority patent/GB2423684B/en
Priority to GB0617870A priority patent/GB2427811B/en
Priority to GB0617871A priority patent/GB2427990B/en
Priority to EP06004002A priority patent/EP1696386A3/en
Priority to DE102006009255A priority patent/DE102006009255B4/de
Publication of JP2006238938A publication Critical patent/JP2006238938A/ja
Application granted granted Critical
Publication of JP4212564B2 publication Critical patent/JP4212564B2/ja
Priority to US12/469,009 priority patent/US20090251464A1/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • G06T15/08Volume rendering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • G06T15/06Ray-tracing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/50Clinical applications
    • A61B6/507Clinical applications involving determination of haemodynamic parameters, e.g. perfusion CT

Abstract

【課題】MIP処理等の医療画像の計算時に、石灰化領域等の障害領域を除去しつつ、障害領域の輪郭を動的に判断して描画することができる画像処理方法を提供する。
【解決手段】注目する組織、例えば血液のボクセル値よりいくぶん大きいしきい値Tを決定する(ステップS11)。次に、仮想光線を投射し(ステップS12)、仮想光線の上のボクセル値を配列A1(元配列)として取得する(ステップS13)。次に、配列A1のしきい値T以上の値をしきい値Tで折り返した配列A2(置き換え配列)を作成し(ステップS14)、配列A2上で一部データ、例えば、石灰化領域の中心部分に相当する折り返えされたデータを除外する(ステップS15)。次に、配列A2上で最大値となる値M1を求め(ステップS16)、配列A1上で値M1に対応する値M2を求める(ステップS17)。そして、値M2をこの仮想光線に対するピクセル値とする(ステップS18)。
【選択図】 図22

Description

本発明は、ボリュームレンダリングによる画像処理方法に関する。
従来、コンピューター断層撮影(CT:computed tomography)装置、核磁気共鳴映像(MRI:magnetic resonance imaging)装置等で得られる3次元画像を、所望の方向に投影して投影画像を得ることが行われている。このような投影画像を得るための処理としてはボリュームレンダリングが広く用いられている。ボリュームレンダリングには例えば、投影方向について最大画素値を抽出して投影するMIP(Maximum Intensity Projection)処理や最小画素値を抽出して投影するMinIP(Minimum Intensity Projection)処理、投影方向に仮想光線を投射し物体からの反射光を計算するレイキャスト法等が知られている。
図31は、MIP処理の説明図であり、描画対象のボクセル値に対応する3Dデータと、表示データとして選択される最大値との関係を示す。MIP処理では、図の矢印で示す投影線上の3Dデータの最大値を表示データとするため、図31(a),(b),(c),(d)では、それぞれ3Dデータの最大値である4,8,8,8が表示データとなる。
図32は、レイキャスト(Raycast)画像とMIP画像を示す。図32(a)に示すレイキャスト画像は、ボユームレンダリング画像の一種であり、仮想光線上の複数のボクセルからの反射光を蓄積することによって画素を決定する。このため輪郭の描画に優れ写実的な画像が得られる。また、仮想光線がボクセルデータ間を通過する場合はボクセルデータそのものではなく、ボクセルデータを補間した情報を元に計算を行うこともある。
一方、図32(b)に示すMIP画像は、前述のように、仮想光線上から単一のボクセルを選択してその情報を用いて画素を決定し、ボクセル値をそのまま描画するので客観性に優れ、計算が高速である。このため、血管の描画にはMIP画像を用いることが多い。なお、補間値を用いる場合は、複数のボクセルを参照することになるが仮想光線上の一点の情報のみを用いるという点では変わらない。ただ、ボクセル値に特徴の無い臓器の描画が困難になる場合がある。
図33は、MIP画像において、血管内に石灰化領域50が付着し血流52に障害が生じた箇所の状況を説明するための図である。また、同図の(a)、(b)は、血管の同じ箇所を90度異なる方向から観察した場合を示す。
図33(a)に示すMIP画像では、血管内における高CT値の石灰化領域50の大きさは把握できるが、石灰化領域50で塞がれている狭窄箇所51の血流52を正確に測定できない場合がある。また、図33(b)に示すMIP画像では、石灰化領域50が障害になって血流52の観察が難しくなり、実際には石灰化領域50の奥あるいは手前に血流52があっても、その血流52を観察することができない。
図34は、血管内に高CT値の石灰化領域が存在する箇所における仮想光線上のボクセル値の変化を示す図である。仮想光線上において石灰化領域のボクセル値は、値が大きくかつ鋭角状のピークを持った値となる。一方、血流のボクセル値は、値が小さくなだらかな形状の値となる。
したがって、MIP画像では、仮想光線上のボクセル値の最大値を表示データとするため、石灰化領域が存在する血管を観察する場合に、ボクセル値が大きい石灰化領域が表示され、石灰化領域の奥あるいは手前に存在する血流を表示することができない。
図35は、MIP画像において石灰化領域の奥あるいは手前に存在する血流を観察する場合における従来の解決法を説明するための図である。従来は、同図に示すように、石灰化領域のCT値をなんらかの値に置き換えた置き換えボリュームデータ(例えば空気のボリュームデータ)を作成し、石灰化領域に対応するボクセル値を下げて血流を表示できるようにしていた。或いは、石灰化領域に相当する領域を描画対象から除くことによってもほぼ同様の効果が得られていた。しかし、上記二つの手法であってはあらかじめ石灰化領域を特定する領域抽出処理を行う必要があった。
すなわち、従来の解決法にあっては、ボリュームレンダリングの前段階で、しきい値あるいは他のアルゴリズムによる領域抽出法を用いて石灰化領域の検出を行い、領域抽出の結果を用いて、ボリュームデータを改変(石灰化領域を除去)したり、あるいはマスクデータ(マスクボリュームを用いて非描画領域を指定)を作成し、血流を表示できるようにしていた。
図36は、従来のレイキャスト法における画面上の各ピクセルの計算を示すフローチャートである。従来のレイキャスト法では、以下の計算を画像上の全ピクセル分行う。まず、投影位置より、投影開始点O(x,y,z)、およびサンプリング間隔ΔS(x,y,z)を設定する(ステップS201)。
次に、反射光Eを「0」、残存光Iを「1」、現在計算位置X(x,y,z)を「O」に初期化する(ステップS202)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置、補完ボクセル値Vを求める(ステップS203)。また、補完ボクセル値Vに対応する不透明度αを得る(ステップS204)。この場合、α=f( V)なる関数をあらかじめ準備しておく(ステップS212)。
次に、補完ボクセル値Vに対応するカラー値Cを得る(ステップS205)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置およびグラディントGを求め、光線方向X-OとグラディントGよりシェーディング係数βを求める(ステップS206)。
次に、X(x,y,z)位置の減衰光D(D=I*α)及び部分反射光F(F=β*D*C)を計算する(ステップS207)。そして、反射光Eおよび残存光Iを更新し(I=I-D,E=E+F)、現在計算位置を進行させる(X=X+ΔS)(ステップS208)。
次に、Xが終了位置まで来たかどうか、または残存光Iが「0」になったかどうかを判断し(ステップS209)、Xが終了位置でなく、残存光Iが「0」でない場合(no)は、X(x,y,z)にΔS(x,y,z)を加算し、現在計算位置を前進させ(ステップS210)、ステップS203以降の処理を繰り返す。一方、Xが終了位置まで来たか、または残存光Iが「0」になった場合(yes)は、反射光Eを計算ピクセルのピクセル値として計算を終了する(ステップS211)。
図37は、従来のMIP処理において、画面上の各ピクセルを計算するフローチャートを示す。従来のMIP処理では以下の計算を画像上の全ピクセル分行う。まず、投影位置より、投影開始点O(x,y,z)、およびサンプリング間隔ΔS(x,y,z)を設定する(ステップS221)。
次に、最大値Mをシステム最小値に、また現在計算位置X(x,y,z)を「O」に初期化する(ステップS222)。そして、X(x,y,z)位置周辺のボクセルデータよりX位置、補完ボクセル値Vを求める(ステップS223)。
次に、最大値Mと補完ボクセル値Vを比較し(ステップS224)、最大値Mが補完ボクセル値Vより小さい場合(yes)は、最大値Mを補完ボクセル値Vとする(ステップS225)。そして、現在計算位置Xが終了位置まで来たかどうかを判断し(ステップS226)、現在計算位置Xが終了位置まで来ていない場合(no)は、X(x,y,z)にΔS(x,y,z)を加算して現在計算位置を前進させ(ステップS227)、ステップS223以降の処理を繰り返す。一方、現在計算位置Xが終了位置まで来た場合(yes)は、最大値Mを計算ピクセルのピクセル値とする(ステップS228)。
また、特許文献1にあっては障害領域を含まない第2ボリュームデータを作成し、第2ボリュームデータにおける最大値を求め、最大値の位置のに対応する元のボリュームデータの位置における値を描画に用いている。
Medical diagnostic method for the two-dimensional imaging of structures US-6205350, U.S. Philips Corporation
しかしながら、上記従来の方法にあっては、ボリュームデータの置き換えにより石灰化領域等の障害領域が除去され、障害領域そのものの情報が完全に失われてしまう。また、厳密に障害領域のみを排除することは難しく、血流を正確に描画することができない。また、抽出領域はボクセル単位で設定されるので領域の境界にエイリアスが生じ画像の劣化を伴う。また、マスク情報や第2ボリュームデータの保持はシステムによけいなメモリ負荷をかけ、及び、ボリュームデータを改変した場合は元のデータとの比較が困難になる。また、個々の障害領域の抽出に時間がかかり観察者の主観要素が大きい。特に観察者の主観要素に左右されることによって観察者毎の再現性が低く客観的な診断情報としての普遍性に欠くことにつながっていた。その為に診断の現場では用いることが難しく現実にはあまり利用されていないという問題がある。
図38は、従来のMIP画像における問題点を説明するための図である。従来の方法にあっては、同図の(a),(b)に示すように、石灰化領域61の前後の血流60を観察するために石灰化領域61を除去するが、その場合に、血流60が存在する部分62も削除してしまう。また、従来の方法では、石灰化領域61が全く表示されなくなるので、病変部の判断が難しくなるとともに、必要な領域も削除されることが多いので信頼性が低くなってしまう。
この場合は、図38(c)に示すように、石灰化領域を除去した跡の部分63の情報が必要であり、特に、石灰化領域の輪郭部分64の情報が必要である。すなわち、石灰化領域の中身は表示せずに石灰化領域の輪郭のみが表示されると診断に有効である。
なお、石灰化領域は3次元領域であるのでその境界面は3次元空間上の面を構成する。そのために従来のマスク適用やボリュームの改変で石灰化を表現しようとした場合は画像の個々の画素は画素を構成する仮想光線と3次元的境界面との交わりを表現するので2次元的な輪郭を表現することができない。一方、画像を見ながら診察する場合は、画像計算時に石灰化領域の2次元的輪郭部分および石灰化領域の周辺、特に手前及び奥の情報が必要である。2次元的輪郭部分は、三次元的に仮想光線が石灰化領域の縁をかすめる箇所のみを描画できれば診断に有効である。
本発明は、上記従来の事情に鑑みてなされたものであって、MIP処理等の医療画像の計算時に、石灰化領域等の障害領域を除去しつつ、障害領域の2次元的輪郭を動的に判断して描画することができる画像処理方法を提供することを目的としている。
本発明の画像処理方法は、ボリュームレンダリングによる画像処理方法であって、第1のベクトル情報と第2のベクトル情報とを用いて、仮想光線上の、位置関係が相互交換可能な1以上の点の値を選択し、選択した1以上の点の値を用いて画素値を決定する。従来はMIP法などの反射光を計算しないボリュームレンダリング法ではベクトル情報を用いることがなかったのに対し、上記構成によれば、ベクトル情報を用いることによって、例えば、仮想光線方向に依存した2次元的効果(2次元的輪郭)を画像に付加することができる。これにより、例えば、2次元的輪郭を利用して障害領域の形状情報とその周辺の情報を一度に観察することが可能になる。
本発明の画像処理方法の処理対象はボリュームレンダリングである。また、画素値に対して「ボクセル」ではなくて「点」としたのは、実際にはボクセル情報を補間した値を用いることが多いためである。また、1以上の点の位置関係が相互交換可能としたのは、MIPのみを対照するのであれば「1の点」を用いても良いが、例えば「Top10MIP」と言う上位10のデータの平均値を表示する方法があり、Top10MIPは複数の点の情報を用いるので「1の点」の定義から外れるためである。また、MIP法では奥行き方向にボリュームデータを反転しても同じ画像がえられる。MIP法などでは仮想光線上の値の位置関係を用いないので不透明度と言う値が定義できない。これは、平均値法であっても同様で、平均値を構成する各値を入れ替えても結果には影響しない。一方、加重平均はボクセルの画素に対する寄与度として表現できるが、これも加重ごと入れ替えれば各値を入れ替えても結果には影響しない。よって、寄与度は多値であってもかまわない。このようなボリュームレンダリング法にあっては光線のシミュレートとは異なる方法で画像を作成するので従来から勾配要素を計算式に組み込むことは考えられていなかった。これに対して、レイキャスト法では不透明度の変更という実施形態になるが、この場合は仮想光線の光量減衰を扱うのでボリュームデータの前後を反転すると異なる画像がえられる。また、第1のベクトル情報は、主には仮想光線の方向ベクトルであるが、それ以外の方向ベクトルを用いるとそれに対応する画像が得られる。一方、第2のベクトル情報は、主にはボクセルの勾配(補間された勾配を含む)であるが、動き情報などのボクセルに結び付けられるベクトル情報があった場合はそれらを使うことができる。
また、本発明の画像処理方法において、前記第1のベクトル情報は、仮想光線方向ベクトルである。また、本発明の画像処理方法において、前記第2のベクトル情報は、勾配情報である。また、本発明の画像処理方法において、選択する点の数は、1である。
また、本発明の画像処理方法は、1以上の点の値を選択するのに仮想光線上のデータを置き換えたデータを用いる。また、本発明の画像処理方法において、前記置き換えデータの値は、元データ値をしきい値で折り返したものである。
また、本発明の画像処理方法は、1以上の点の値を選択するのに第2のベクトル情報の大きさを用いる。また、本発明の画像処理方法は、1以上の点の値を選択するのに第1のベクトル情報と第2のベクトル情報の間の角度関係を用いる。また、本発明の画像処理方法は、描画対象に含まれる物体の2次元的輪郭を表示させる。
また、本発明の画像処理方法は、ボリュームレンダリングによる画像処理方法であって、画素値の決定に当たって仮想光線上の1以上の点の値を利用し、仮想光線上の前記1以上の点の値の画素値を決定する処理への寄与度を変更する。従来のボリュームレンダリング法では、画素値を決定する処理への寄与度はマスク情報などを通して与えられていた為、寄与度は3次元的に決定され2次元的効果を付加することができなかったのに対し、上記構成によれば、仮想光線上で寄与度を判断するので2次元的効果を画像に付加することができる。
また、本発明の画像処理方法において、前記寄与度は0である。また、本発明の画像処理方法は、前記1以上の点の寄与度を決定するのに仮想光線上のデータを置き換えたデータを用いる。また、本発明の画像処理方法は、前記置き換えデータの値は元データ値をしきい値で折り返したものである。
また、本発明の画像処理方法は、前記1以上の点の寄与度を決定するのに、点の位置に対応するボリューム上の勾配ベクトルと仮想光線の進行方向ベクトルとを更に用いるものである。また、本発明の画像処理方法は、前記1以上の点の寄与度を決定するのに、仮想光線上のボクセル値の変化を用いるものである。
また、本発明の画像処理方法は、描画対象に含まれる物体の2次元的輪郭を表示させる。また、本発明の画像処理方法は、描画対象に含まれる物体を非表示とする。
また、本発明の画像処理方法は、画像を並べて表示、重ねて表示、または差分を表示させる。また、本発明の画像処理方法は、ユーザの指定した箇所にのみ適応する。また、本発明の画像処理方法は、画面上に設けた小窓内にのみ適用する。
また、本発明の画像処理方法は、輪郭の表現を連続的に変更する。また、本発明の画像処理方法は、並列処理を行う。また、本発明の画像処理方法は、GPUにより処理を行う。また、本発明の画像処理方法は、パラメータが変更可能なGUIにより処理を行う。
また、本発明の画像処理方法は、MIP、MinIP、Raysum、または平均値法により処理を行う。また、本発明の画像処理方法は、IP、MinIP、Raysum、平均値法、またはレイキャストにより処理を行う。また、本発明の画像処理方法は、距離で判断し、1Dグラディエントを行う。また、本発明の画像処理方法は、仮想光線を含む断面画像に選択された点を表示させる。
また、本発明の画像処理装置は、画素値の決定に当たって仮想光線上のデータの1以上の点の値を用い、1以上の点の位置関係が相互交換可能なボリュームレンダリング画像を表示する画像処理装置であって、描画対象に含まれる物体の2次元的輪郭を表示させるものである。
本発明の画像処理方法によれば、仮想光線の投射時に、仮想光線上のボクセルデータにおける一部のボクセルデータ、例えば障害領域の2次元的中心部分を非表示データとし、障害領域の2次元的輪郭を表示データとすることにより、診断対象を的確に表示することができる。
図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のベクトル情報の角度を計算したが、ここで角度は負の値であっても直角より大きくても良い。このようにすれば手間側の壁のみを表示するようなことができるようになる。
本発明は、MIP処理等の医療画像の計算時に、石灰化領域等の障害領域を除去しつつ、障害領域の輪郭を動的に判断して描画することができる画像処理方法に利用可能である。
本実施形態の画像処理方法において利用するボクセル値の特徴を説明するための図 仮想光線が石灰化領域20の中心部分、その縁21および血流22を通過する場合のボクセル値の変化を示す図 本実施形態の画像処理方法により取得するプロファイルパターン(1)を示す図 本実施形態の画像処理方法により取得するプロファイルパターン(2)を示す図 本実施形態の画像処理方法において、MIP処理に勾配(グラディエント)を用いる場合の説明図 本実施形態の画像処理方法において、除外可能の判断を行う場合の説明図 本実施形態の画像処理方法において、勾配(グラディエント)を用いたMIP処理(折り返し無し)のフローチャート 本実施形態の画像処理方法において、MIP処理に勾配(グラディエント)および折り返しを用いる場合の説明図 勾配要素によるプロファイルパターンの改変を示す説明図 折り返しによるプロファイルパターンの改変を示す説明図 LUT関数例と石灰化領域の表示のされ方の説明図 本実施形態の画像処理方法において、勾配(グラディエント)および折り返しを用いたMIP処理のフローチャート 本実施形態の画像処理方法において、石灰化領域の中心部分を削除してその前後の血流を表示し、かつ石灰化領域の縁を表示する手順1を示す図 本実施形態の画像処理方法における手順2を示す図 本実施形態の画像処理方法において、仮想光線に対応するボクセル値を除外するかどうかの判断方法(1)の説明図 本実施形態の画像処理方法において、仮想光線に対応するボクセル値を除外するかどうかの判断方法(2)の説明図 本実施形態の画像処理方法において、仮想光線に対応するボクセル値を除外するかどうかの判断方法(3)の説明図 本実施形態の画像処理方法において、しきい値を用いた折り返し以外の変換方法を用いて石灰化領域の縁を検出する場合の説明図 本実施形態の画像処理方法において、石灰化領域の縁のデータとして折り返す前の元データを用いる場合の手順(3)を示す図 本実施形態の画像処理方法における望ましい実施例を示す図 通常のMIP画像(a)と本実施形態の画像処理方法によって作成したMIP画像(b)の例を示す図 本実施形態の画像処理方法において、画像の各ピクセル値を求める処理を示すフローチャート 本実施形態の画像処理方法において、データ変化を用いてデータ除外の判断を行う場合のフローチャート 本実施形態の画像処理方法において、勾配を用いてデータ除外を判断する場合のフローチャート 透明度を変更する場合(レイキャスト法に対する応用1)の説明図 勾配を用いて透明度を変更する場合(レイキャスト法に対する応用)のフローチャート 透明度を変更する場合(レイキャスト法に対する応用2)の説明図 変化を用いて透明度を変更する場合(レイキャスト法に対する応用)のフローチャート 勾配処理と折り返し処理を行うMIPフローチャート 交錯する血流を可視化する場合の説明図 ボクセル値および3Dデータに対するMIP処理の説明図 レイキャスト(Raycast)画像とMIP画像の例 MIP画像において、血管内に石灰化領域50が付着し血流52に障害が生じた箇所の状況を説明するための図 血管内に高CT値の石灰化領域が存在する箇所における仮想光線上のボクセル値の変化を示す図 MIP画像において石灰化領域の奥あるいは手前に存在する血流を観察する場合における従来の解決法を説明するための図 従来のレイキャスト法における画面上の各ピクセルの計算を示すフローチャート 従来のMIP処理において画面上の各ピクセルを計算するフローチャート 従来のMIP画像における問題点を説明するための図
符号の説明
10 障害領域
11 障害領域の縁
20 石灰化領域
21 石灰化領域の縁
22 血流

Claims (29)

  1. ボリュームレンダリングによる画像処理方法であって、
    第1のベクトル情報と第2のベクトル情報とを用いて、仮想光線上の、位置関係が相互交換可能な1以上の点の値を選択し、
    選択した1以上の点の値を用いて画素値を決定する画像処理方法。
  2. 請求項1記載の画像処理方法であって、
    前記第1のベクトル情報は、仮想光線方向ベクトルである画像処理方法。
  3. 請求項1記載の画像処理方法であって、
    前記第2のベクトル情報は、勾配情報である画像処理方法。
  4. 請求項1記載の画像処理方法であって、
    選択する点の数は1である画像処理方法。
  5. 請求項1記載の画像処理方法であって、
    1以上の点の値を選択するのに仮想光線上のデータを置き換えたデータを用いる画像処理方法。
  6. 請求項5記載の画像処理方法であって、
    前記置き換えデータの値は、元データ値をしきい値で折り返したものである画像処理方法。
  7. 請求項1記載の画像処理方法であって、
    1以上の点の値を選択するのに第2のベクトル情報の大きさを用いる画像処理方法。
  8. 請求項1記載の画像処理方法であって、
    1以上の点の値を選択するのに第1のベクトル情報と第2のベクトル情報の間の角度関係を用いる画像処理方法。
  9. 請求項1記載の画像処理方法であって、
    描画対象に含まれる物体の2次元的輪郭を表示させる画像処理方法。
  10. ボリュームレンダリングによる画像処理方法であって、
    画素値の決定に当たって仮想光線上の1以上の点の値を利用し、
    仮想光線上の前記1以上の点の値の画素値を決定する処理への寄与度を変更する画像処理方法。
  11. 請求項10記載の画像処理方法であって、
    前記寄与度は0である画像処理方法。
  12. 請求項10記載の画像処理方法であって、
    前記1以上の点の寄与度を決定するのに仮想光線上のデータを置き換えたデータを用いる画像処理方法。
  13. 請求項12記載の画像処理方法であって、
    前記置き換えデータの値は元データ値をしきい値で折り返したものである画像処理方法。
  14. 請求項10記載の画像処理方法であって、
    前記1以上の点の寄与度を決定するのに、点の位置に対応するボリューム上の勾配ベクトルと仮想光線の進行方向ベクトルとを更に用いる画像処理方法。
  15. 請求項10記載の画像処理方法であって、
    前記1以上の点の寄与度を決定するのに、仮想光線上のボクセル値の変化を用いる画像処理方法。
  16. 請求項10記載の画像処理方法であって、
    描画対象に含まれる物体の2次元的輪郭を表示させる画像処理方法。
  17. 請求項10記載の画像処理方法であって、
    描画対象に含まれる物体を非表示とする画像処理方法。
  18. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    画像を並べて表示、重ねて表示、または差分を表示させる画像処理方法。
  19. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    ユーザの指定した箇所にのみ適応する画像処理方法。
  20. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    画面上に設けた小窓内にのみ適用する画像処理方法。
  21. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    輪郭の表現を連続的に変更する画像処理方法。
  22. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    並列処理を行う画像処理方法。
  23. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    GPUにより処理を行う画像処理方法。
  24. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    パラメータが変更可能なGUIにより処理を行う画像処理方法。
  25. 請求項1記載の画像処理方法であって、
    MIP、MinIP、Raysum、または平均値法により処理を行う画像処理方法。
  26. 請求項10記載の画像処理方法であって、
    MIP、MinIP、Raysum、平均値法、またはレイキャストにより処理を行う画像処理方法。
  27. 請求項1ないし17のいずれか一項記載の画像処理方法であって、
    仮想光線を含む断面画像に前記1以上の選択された点を表示させる画像処理方法。
  28. コンピュータに、請求項1ないし27のいずれか一項記載の各ステップを実行させるための画像処理プログラム。
  29. 画素値の決定に当たって仮想光線上のデータの1以上の点の値を用い、1以上の点の位置関係が相互交換可能なボリュームレンダリング画像を表示する画像処理装置であって、描画対象に含まれる物体の2次元的輪郭を表示させる画像処理装置。
JP2005054863A 2005-02-28 2005-02-28 画像処理方法および画像処理プログラム Active JP4212564B2 (ja)

Priority Applications (8)

Application Number Priority Date Filing Date Title
JP2005054863A JP4212564B2 (ja) 2005-02-28 2005-02-28 画像処理方法および画像処理プログラム
US11/346,058 US7853057B2 (en) 2005-02-28 2006-02-02 Image processing method and device using a virtual ray and vector information
GB0617870A GB2427811B (en) 2005-02-28 2006-02-07 Image processing method and image processing device
GB0617871A GB2427990B (en) 2005-02-28 2006-02-07 Image processing method and image processing device
GB0602433A GB2423684B (en) 2005-02-28 2006-02-07 Image processing method and image processing device
EP06004002A EP1696386A3 (en) 2005-02-28 2006-02-28 Method and device for volume rendering
DE102006009255A DE102006009255B4 (de) 2005-02-28 2006-02-28 Bildverarbeitungsverfahren und Bildverarbeitungsgerät
US12/469,009 US20090251464A1 (en) 2005-02-28 2009-05-20 Image processing method and image processing device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2005054863A JP4212564B2 (ja) 2005-02-28 2005-02-28 画像処理方法および画像処理プログラム

Related Child Applications (2)

Application Number Title Priority Date Filing Date
JP2008125016A Division JP2008200524A (ja) 2008-05-12 2008-05-12 画像処理方法
JP2008125015A Division JP2008246215A (ja) 2008-05-12 2008-05-12 画像処理方法

Publications (2)

Publication Number Publication Date
JP2006238938A true JP2006238938A (ja) 2006-09-14
JP4212564B2 JP4212564B2 (ja) 2009-01-21

Family

ID=36119639

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2005054863A Active JP4212564B2 (ja) 2005-02-28 2005-02-28 画像処理方法および画像処理プログラム

Country Status (5)

Country Link
US (2) US7853057B2 (ja)
EP (1) EP1696386A3 (ja)
JP (1) JP4212564B2 (ja)
DE (1) DE102006009255B4 (ja)
GB (3) GB2427990B (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008161688A (ja) * 2007-01-02 2008-07-17 General Electric Co <Ge> 自動冠動脈カルシウム検出およびラベリングのシステム
WO2012144266A1 (ja) * 2011-04-20 2012-10-26 株式会社日立メディコ 医用画像表示装置、医用画像表示方法

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090204338A1 (en) * 2008-02-13 2009-08-13 Nordic Bioscience A/S Method of deriving a quantitative measure of the instability of calcific deposits of a blood vessel
EP2524652B1 (en) * 2010-01-15 2014-11-05 Hitachi Medical Corporation Ultrasonic diagnostic device and ultrasonic image display method
JP2012016575A (ja) * 2010-06-07 2012-01-26 Toshiba Corp 画像処理装置及び医用画像診断装置
JP5576782B2 (ja) * 2010-12-16 2014-08-20 オリンパス株式会社 画像処理装置、画像処理方法、及び画像処理プログラム
US8867806B2 (en) * 2011-08-01 2014-10-21 Impac Medical Systems, Inc. Method and apparatus for correction of errors in surfaces
KR101859413B1 (ko) * 2011-12-26 2018-05-21 삼성디스플레이 주식회사 표시 모듈 및 이를 포함하는 표시 장치
KR102240564B1 (ko) * 2014-07-29 2021-04-15 삼성전자주식회사 영상 렌더링 장치 및 방법
US9659405B2 (en) * 2015-04-01 2017-05-23 Toshiba Medical Systems Corporation Image processing method and apparatus
US10223813B2 (en) * 2015-08-13 2019-03-05 InstaRecon Method and system for reprojection and backprojection for tomography reconstruction
CN111684495A (zh) * 2017-12-22 2020-09-18 奇跃公司 用于在混合现实系统中管理和显示虚拟内容的方法和系统
IL301443A (en) 2018-02-22 2023-05-01 Magic Leap Inc A browser for mixed reality systems
IL301281A (en) 2018-02-22 2023-05-01 Magic Leap Inc Shoot an object with physical manipulation
US10991133B2 (en) * 2019-01-25 2021-04-27 Siemens Healthcare Gmbh Volume rendering from three-dimensional medical data using quantum computing
JP7127212B2 (ja) * 2019-04-25 2022-08-29 富士フイルム株式会社 画像の向き設定装置、方法およびプログラム

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS62164174A (ja) * 1986-01-16 1987-07-20 Hitachi Ltd 三次元表示方式
JPH03147082A (ja) * 1989-11-01 1991-06-24 Hitachi Ltd 3次元データ表示方法
JPH10171976A (ja) * 1996-10-07 1998-06-26 Ge Yokogawa Medical Syst Ltd 画像処理方法及び画像処理装置
JPH11283018A (ja) * 1998-03-30 1999-10-15 Shimadzu Corp 医用画像処理装置
JP2002312809A (ja) * 2001-04-12 2002-10-25 Ziosoft Inc 3次元画像表示方法
JP2003153894A (ja) * 2001-11-26 2003-05-27 Ziosoft Inc 3次元画像処理方法、装置およびプログラム
JP2004141514A (ja) * 2002-10-28 2004-05-20 Toshiba Corp 画像処理装置及び超音波診断装置

Family Cites Families (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE3903838A1 (de) * 1988-02-09 1989-08-17 Toshiba Kawasaki Kk Verfahren und einrichtung zum darstellen dreidimensionaler bilder
US5204627A (en) * 1991-03-14 1993-04-20 Wisconsin Alumni Research Foundation Adaptive NMR angiographic reprojection method
US5368033A (en) * 1993-04-20 1994-11-29 North American Philips Corporation Magnetic resonance angiography method and apparatus employing an integration projection
JP3570576B2 (ja) * 1995-06-19 2004-09-29 株式会社日立製作所 マルチモダリティに対応した3次元画像合成表示装置
US6072497A (en) * 1997-05-30 2000-06-06 Hewlett-Packard Company Volumetric pre-clipping method that guarantees minimal number of sample points through a volume
DE19806728A1 (de) * 1998-02-18 1999-08-19 Philips Patentverwaltung Verfahren zur zweidimensionalen Abbildung von Strukturen für die medizinische Diagnostik
JP4421016B2 (ja) * 1999-07-01 2010-02-24 東芝医用システムエンジニアリング株式会社 医用画像処理装置
US6654012B1 (en) 1999-10-01 2003-11-25 Terarecon, Inc. Early ray termination in a parallel pipelined volume rendering system
DE10012174B4 (de) 2000-03-13 2005-02-17 Siemens Ag Abbildungsverfahren für einen Volumendatensatz
DE10012172A1 (de) 2000-03-13 2001-09-27 Siemens Ag Abbildungsverfahren für einen Volumendatensatz
AU2001256992A1 (en) * 2000-04-07 2001-10-23 Stephen R. Aylward Systems and methods for tubular object processing
US7022073B2 (en) * 2003-04-02 2006-04-04 Siemens Medical Solutions Usa, Inc. Border detection for medical imaging
JP2005054863A (ja) 2003-08-01 2005-03-03 Showa Corp 油圧緩衝器

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS62164174A (ja) * 1986-01-16 1987-07-20 Hitachi Ltd 三次元表示方式
JPH03147082A (ja) * 1989-11-01 1991-06-24 Hitachi Ltd 3次元データ表示方法
JPH10171976A (ja) * 1996-10-07 1998-06-26 Ge Yokogawa Medical Syst Ltd 画像処理方法及び画像処理装置
JPH11283018A (ja) * 1998-03-30 1999-10-15 Shimadzu Corp 医用画像処理装置
JP2002312809A (ja) * 2001-04-12 2002-10-25 Ziosoft Inc 3次元画像表示方法
JP2003153894A (ja) * 2001-11-26 2003-05-27 Ziosoft Inc 3次元画像処理方法、装置およびプログラム
JP2004141514A (ja) * 2002-10-28 2004-05-20 Toshiba Corp 画像処理装置及び超音波診断装置

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008161688A (ja) * 2007-01-02 2008-07-17 General Electric Co <Ge> 自動冠動脈カルシウム検出およびラベリングのシステム
WO2012144266A1 (ja) * 2011-04-20 2012-10-26 株式会社日立メディコ 医用画像表示装置、医用画像表示方法
CN103491877A (zh) * 2011-04-20 2014-01-01 株式会社日立医疗器械 医用图像显示装置、医用图像显示方法
US9135702B2 (en) 2011-04-20 2015-09-15 Hitachi Medical Corporation Image display device for medical applications, image display method for medical applications
JP5980776B2 (ja) * 2011-04-20 2016-08-31 株式会社日立製作所 医用画像表示装置、医用画像表示方法

Also Published As

Publication number Publication date
GB0617871D0 (en) 2006-10-18
EP1696386A2 (en) 2006-08-30
GB2427990B (en) 2007-09-19
US20060193510A1 (en) 2006-08-31
US7853057B2 (en) 2010-12-14
GB2423684B (en) 2007-11-07
GB2427990A (en) 2007-01-10
GB0602433D0 (en) 2006-03-22
GB0617870D0 (en) 2006-10-18
JP4212564B2 (ja) 2009-01-21
EP1696386A3 (en) 2011-09-07
DE102006009255A1 (de) 2006-09-07
DE102006009255B4 (de) 2007-10-11
GB2427811B (en) 2007-10-31
GB2427811A (en) 2007-01-03
US20090251464A1 (en) 2009-10-08
GB2423684A (en) 2006-08-30

Similar Documents

Publication Publication Date Title
JP4212564B2 (ja) 画像処理方法および画像処理プログラム
US8150111B2 (en) Methods, systems, and computer program products for processing three-dimensional image data to render an image from a viewpoint within or beyond an occluding region of the image data
CN107924580B (zh) 在医学成像中的表面体积混合模块的可视化
US9563978B2 (en) Image generation apparatus, method, and medium with image generation program recorded thereon
JP5415068B2 (ja) カーブした細長い構造の切断面の視覚化
US20080117203A1 (en) Methods and Apparatus for Visualizing Data
EP2856428B1 (en) Segmentation highlighter
JP2006075602A (ja) 血管構造の3d画像データセットからなるプラーク沈着の可視化方法
US20110082667A1 (en) System and method for view-dependent anatomic surface visualization
JP2006167287A (ja) 血管狭窄率解析システム
JP2010528750A (ja) 管状構造の検査
JP2012511380A (ja) 血管分析
JP2006055213A (ja) 画像処理装置、及びプログラム
JP2008289767A (ja) 画像処理方法及び画像処理プログラム
JP2006326312A (ja) 多分岐された血管および該血管の周囲環境を単一のイメージに同時に投影する方法
JP4885042B2 (ja) 画像処理方法および装置ならびにプログラム
JP2015515296A (ja) 対象物の画像情報の提供
JP2005523053A (ja) 対象物表面の折り畳まれた解剖学的部分の視覚化のための医用ビューイングシステム及び画像処理方法
JP2008200524A (ja) 画像処理方法
JP2008246215A (ja) 画像処理方法
Bartz et al. Visualization and exploration of segmented anatomic structures
JP2019010382A (ja) 画像位置合わせ装置、方法およびプログラム
Abrahams et al. Detecting blindspots in colonoscopy by modelling curvature
CN111724388B (zh) 医学图像数据的可视化
US20220343586A1 (en) Method and system for optimizing distance estimation

Legal Events

Date Code Title Description
RD04 Notification of resignation of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7424

Effective date: 20071129

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20080122

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20080312

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20080512

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20080716

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20080905

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20081008

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20081028

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20111107

Year of fee payment: 3

R150 Certificate of patent or registration of utility model

Ref document number: 4212564

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20121107

Year of fee payment: 4

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20121107

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20131107

Year of fee payment: 5

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250