JP3542779B2

JP3542779B2 - 3D image display method

Info

Publication number: JP3542779B2
Application number: JP2001114133A
Authority: JP
Inventors: 和彦松本
Original assignee: Ziosoft Inc
Current assignee: Ziosoft Inc
Priority date: 2001-04-12
Filing date: 2001-04-12
Publication date: 2004-07-14
Anticipated expiration: 2021-04-12
Also published as: JP2002312809A

Description

【０００１】
【発明の属する技術分野】
本発明はレイキャスティングによりボリュームレンダリングを行う３次元画像表示方法、装置およびプログラムに関し、特に、レイキャスティングにおいて使用されるボクセルデータのグラディエントの算出方法に関する。
【０００２】
【従来の技術】
コンピュータを用いた画像処理技術の進展により人体の内部構造を直接観測することを可能にしたＣＴ（ＣｏｍｐｕｔｅｄＴｏｍｏｇｒａｐｈｙ）の出現は医療分野に革新をもたらした技術であり、生体の断層画像を用いた医療診断が広く行われている。さらに近年は、断層画像だけではわかり難い物体内部の３次元構造を可視化する技術として、ＣＴにより得られる物体の３次元ディジタルデータから３次元構造のイメージを直接描画するボリュームレンダリングが注目されている。
【０００３】
ボリュームレンダリングの優れた手法としてレイキャスティングが知られている。レイキャスティングは、仮想始点から物体に対して仮想光線（レイ）を照射し、物体内部からの仮想反射光の画像を仮想投影面に形成することにより、物体内部の３次元構造を透視するイメージ画像を形成する手法である。レイキャスティングについては、例えば、「新世代３次元ＣＴ診断」（１９９５年１１月１日、株式会社南江堂発行）に基本的な理論が述べられている。
【０００４】
以下にレイキャスティングの要点を説明する。物体の３次元領域の構成単位となる微小単位領域をボクセルと称し、ボクセルの濃度値等の特性を表わす固有のデータをボクセル値と称する。物体全体はボクセル値の３次元配列であるボクセルデータで表現される。通常、ＣＴにより得られる２次元の断層画像データを断層面に垂直な方向に沿って積層し、必要な補間を行うことにより３次元配列のボクセルデータが得られる。
【０００５】
仮想始点から物体に対して照射された仮想光線に対する仮想反射光は、ボクセル値に対して人為的に設定される不透明度に応じて生ずるものとする。さらに、仮想的な表面を捕捉するためにボクセルデータのグラディエントすなわち法線ベクトルを求め、仮想光線と法線ベクトルのなす角の余弦から陰影付けのシェーディング係数を計算する。仮想反射光は、ボクセルに照射される仮想光線の強度にボクセルの不透明度とシェーディング係数を乗じて算出される。
【０００６】
図４はレイキャスティングによる投影画像形成の概念を説明する図である。図４において、４１０は仮想投影面、４２０はボクセルデータで表現される３次元物体であり、４２１、４２２、４２３は、仮想投影面４１０上の座標点Ｏ（ｕ，ｖ）から３次元物体４２０に照射される仮想光線に沿って、仮想反射光を算出するために一定の間隔で刻まれたレイ到達点に位置するボクセルＶ（ｎ−１）、Ｖｎ、Ｖ（ｎ＋１）である。一般には仮想光線の方向およびレイ到達点を刻む間隔は、ボクセルデータの３次元配列の軸方向やボクセル間隔と一致しないので、補間計算により各レイ到達点におけるボクセルの特性値を求める必要がある。
【０００７】
仮想光線上のｎ番目のレイ到達点に位置するボクセルＶｎについて、そのボクセル値に対して与えられた不透明度をαｎ、ボクセルデータのグラディエントから得られたシェーディング係数をβｎとする。ここで、ボクセルＶｎに入射する仮想光線がボクセルＶｎの不透明度により減衰するものとして、その減衰分を減衰光とし、減衰光の分だけ減衰した光を残存光とする。ボクセルＶｎを通過した残存光が次のボクセルＶ（ｎ＋１）に入射する仮想光線となる。また、減衰光の分がボクセルＶｎの反射光を与えるものであると考える。ボクセルＶ（ｎ−１）からの残存光をＩ（ｎ−１）、ボクセルＶｎの残存光をＩｎ、ボクセルＶｎの減衰光をＤｎ、ボクセルＶｎによる部分反射光をＦｎとすると、＊を乗算記号として、以下の式を得る。
【０００８】
Ｄｎ＝αｎ＊Ｉ（ｎ−１）
Ｆｎ＝βｎ＊Ｄｎ＝βｎ＊αｎ＊Ｉ（ｎ−１）
Ｉｎ＝Ｉ（ｎ−１）−Ｄｎ＝（１−αｎ）＊Ｉ（ｎ−１）
【０００９】
仮想投影面４１０上の座標（ｕ，ｖ）に投影される仮想反射光Ｅ（ｕ，ｖ）は、レイ到達点ごとに算出される上記部分反射光を積算したものである。なお、この手法は物体内部の３次元構造の仮想的な透視画像を形成することが目的であり、実世界の物理現象に関する計算ではないので、ボクセルからの反射光に対する減衰は考えない。Ｏ（ｕ，ｖ）から照射する仮想光線の強度をＩ（０）とし、表記法として、Σ（ｉ＝１；ｎ）Ａｉを数列Ａ１からＡｎまでのｎ項の和記号、Π（ｉ＝１；ｎ）Ａｉを数列Ａ１からＡｎまでのｎ項の積記号とすると、次の式を得る。
Ｅ（ｕ，ｖ）＝Σ（ｉ＝１；終端条件）Ｆｉ＝Ｉ（０）＊（Σ（ｉ＝１；終端条件）βｉ＊αｉ＊（Π（ｊ＝１；ｉ−１）（１−αｊ）））
【００１０】
ここで、終端条件は仮想光線が物体を通り抜けるか、残存光が０になったときである。このようにして、仮想投影面４１０上のすべての座標（ｕ，ｖ）についてＥ（ｕ，ｖ）を算出することにより仮想的な３次元イメージの透視画像が形成される。
【００１１】
仮想投影面上の座標点から照射される仮想光線は仮想投影面に垂直な平行光線となるので、この方法は平行投影法と称されている。この方法は外側から見た物体内部の透視画像が得られるので直観的に理解し易いが、例えば医療分野では、臓器の内腔面を外側から観察することになるという欠点がある。これに対して、空間の任意の始点から放射状に仮想光線を照射して透視画像を形成する透視投影法がある。例えば、始点を人体の血管内部に置くことにより、血管内腔面を表面にして表示することが可能になる。
【００１２】
図５は透視投影法によるレイキャスティングの概念を説明する図である。図５において、５１０は空間におかれた任意の始点Ｏ（Ｏｘ，Ｏｙ，Ｏｚ）、５２０はボクセルデータで表現される３次元物体、５３０は仮想投影面である。始点５１０から照射される仮想光線に沿って仮想反射光を算出するためのレイ到達点が一定の間隔で刻まれる。５２１は始点Ｏ（Ｏｘ，Ｏｙ，Ｏｚ）からレイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）に至るレイベクトル、５２２はレイ到達点を刻む間隔を表わすステップベクトルΔＳ＝（ΔＳｘ，ΔＳｙ，ΔＳｚ）、５３１は仮想投影面５３０上の座標（ｕ，ｖ）に投影された仮想反射光Ｅ（ｕ，ｖ）である。ここで、Ｒｘ、Ｒｙ、Ｒｚ、ΔＳｘ、ΔＳｙ、ΔＳｚは（ｕ，ｖ）の関数になる。
【００１３】
図５では仮想投影面５３０が仮想光線の到達方向にあるが、始点５１０を仮想的な眼球の水晶体と考えると、この眼球の網膜５１２上の座標（ｕ，ｖ）に仮想反射光５１１がＥ’（ｕ，ｖ）として投影され、これを始点５１０に対して対称方向に拡大表示したものが仮想投影面５３０上の仮想反射光Ｅ（ｕ，ｖ）であると考えることができる。
【００１４】
ここで、仮想光線上のｎ番目のレイ到達点Ｒｎに位置するボクセルＶｎ（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ）について、そのボクセル値に対して与えられた不透明度をα（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ）、ボクセルデータのグラディエントから得られたシェーディング係数をβ（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ）、ボクセルＶ（ｎ−１）からの残存光をＩ（ｎ−１）、ボクセルＶｎの残存光をＩｎ、ボクセルＶｎの減衰光をＤｎ、ボクセルＶｎによる部分反射光をＦｎとすると、平行投影法と同様に、以下の式を得る。
【００１５】
Ｄｎ＝α（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ）＊Ｉ（ｎ−１）
Ｆｎ＝β（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ）＊Ｄｎ＝β（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ）＊α（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ）＊Ｉ（ｎ−１）
Ｉｎ＝Ｉ（ｎ−１）−Ｄｎ＝（１−α（Ｒｘｎ，Ｒｙｎ，Ｒｚｎ））＊Ｉ（ｎ−１）
【００１６】
平行投影法と同様に、仮想投影面５３０上の座標（ｕ，ｖ）に投影される仮想反射光Ｅ（ｕ，ｖ）を求めると、前述の表記法を用いて次の式を得る。
【００１７】
Ｅ（ｕ，ｖ）＝Σ（ｉ＝１；終端条件）Ｆｉ＝Ｉ（０）＊（Σ（ｉ＝１；終端条件）β（Ｒｘｉ，Ｒｙｉ，Ｒｚｉ）＊α（Ｒｘｉ，Ｒｙｉ，Ｒｚｉ）＊（Π（ｊ＝１；ｉ−１）（１−α（Ｒｘｊ，Ｒｙｊ，Ｒｚｊ））））
【００１８】
終端条件は仮想光線が物体を通り抜けるか、残存光が０になったときである。このようにして、仮想投影面５３０上のすべての座標（ｕ，ｖ）についてＥ（ｕ，ｖ）を算出することにより仮想的な３次元イメージの透視画像が形成される。
【００１９】
以上説明したような原理に基づいて行なわれる従来のレイキャスティングによる３次元画像表示方法の処理フローチャートを図６に示す。以降の説明では、投影法について図５を参照するが、平行投影法と透視投影法とを区別せずに一般化して述べる。また、図６の説明において、処理フローチャートの各ステップを（ステップの符号）で表記する。
【００２０】
図６において、（６０１）は、使用するボクセルデータを格子点（Ｘ，Ｙ，Ｚ）についてＶ（Ｘ，Ｙ，Ｚ）と表現することを示す。（６０２）では、各格子点（Ｘ，Ｙ，Ｚ）のグラディエントＧ（Ｘ，Ｙ，Ｘ）＝（Ｇｘ，Ｇｙ，Ｇｚ）を、格子点（Ｘ，Ｙ，Ｚ）の近傍のボクセルデータから計算してメモリ等に格納する。例えば６近傍であれば、グラディエントＧ（Ｘ，Ｙ，Ｚ）のベクトル成分（成分比）は以下の式になる。
【００２１】
Ｇｘ（Ｘ，Ｙ，Ｘ）＝Ｖ（（Ｘ＋１），Ｙ，Ｚ）−Ｖ（（Ｘ−１），Ｙ，Ｚ）Ｇｙ（Ｘ，Ｙ，Ｘ）＝Ｖ（Ｘ，（Ｙ＋１），Ｚ）−Ｖ（Ｘ，（Ｙ−１），Ｚ）Ｇｚ（Ｘ，Ｙ，Ｘ）＝Ｖ（Ｘ，Ｙ，（Ｚ＋１））−Ｖ（Ｘ，Ｙ，（Ｚ−１））
【００２２】
（６０３）では、投影面の各座標（ｕ，ｖ）への反射光Ｅ（ｕ，ｖ）を求めるために、投影面のピクセルサイズをＷ＊Ｈとして、以降の処理フローに対して、０≦ｕ＜Ｗ、０≦ｖ＜Ｈの繰り返し制御を行う。
【００２３】
（６０４）では、始点座標ＯおよびステップベクトルΔＳ＝（ΔＳｘ，ΔＳｙ，ΔＳｚ）を決定する。始点Ｏは、平行投影法では投影面のピクセル座標（ｕ，ｖ）にあり、透視投影法では空間に任意に設定した固定点である。ステップベクトルΔＳはレイ到達点を刻む間隔を長さに持つベクトルであり、平行投影法では光線方向に一定の向きを持ち、透視投影法では仮想光線方向に応じて向きが変化する。
【００２４】
（６０５）では、投影面の座標（ｕ，ｖ）の現在の値に対して反射光Ｅ（ｕ，ｖ）を求める計算の繰り返し制御を行うために、レイ到達点を始点Ｏに設定し、反射光と残存光について、Ｅ＝０、Ｉ＝Ｉ（０）に初期化する。
【００２５】
（６０６）では、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）の補間ボクセル値Ｖ（Ｒｘ，Ｒｙ，Ｒｚ）を、レイ到達点Ｒの周囲の格子点のＭ個のボクセルデータから算出する。通常はレイ到達点Ｒを囲む立方格子点８個のボクセルデータから、ｘ方向、ｙ方向、ｚ方向それぞれの補間計算を行う。求めた補間ボクセル値から、あらかじめ設定した変換関数により、レイ到達点Ｒの不透明度α（Ｒｘ，Ｒｙ，Ｒｚ）を得る。
【００２６】
次に、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）におけるシェーディング係数β（Ｒｘ，Ｒｙ，Ｒｚ）を算出する。従来のレイキャスティングにおいては、レイ到達点Ｒの周囲の格子点のグラディエントからそれぞれの格子点のシェーディング係数を求め、この格子点のシェーディング係数から補間計算によりレイ到達点Ｒにおける補間シェーディング係数を算出していた。
【００２７】
すなわち、（６０７）では、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）の周囲の格子点のＭ個のグラディエントＧｍ（ｍ＝１，２，．．，Ｍ）、通常はレイ到達点Ｒを囲む立方格子点８個のグラディエントを取得し、それぞれのグラディエントと始点Ｏからレイ到達点Ｒに至るレイベクトルのなす角度θｍ（ｍ＝１，２，．．，Ｍ）からそれぞれの格子点におけるシェーディング係数βｍ＝ＡＢＳ（ｃｏｓθｍ）を得る。ここでＡＢＳ（ｃｏｓθｍ）はｃｏｓθｍの絶対値である。次に、βｍから補間計算によりレイ到達点Ｒにおける補間シェーディング係数β（Ｒｘ，Ｒｙ，Ｒｚ）を算出する。
【００２８】
（６０８）では、不透明度αとシェーディング係数βから、レイ到達点Ｒにおける減衰光Ｄおよび部分反射光Ｆを求める。これらを用いて、（６０９）では、残存光Ｉおよび反射光Ｅを更新し、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）をステップベクトルΔＳだけ進める。（６１０）では、繰り返し制御のための終端条件をチェックする。終端条件が満たされると、（６１１）で、以上の計算から得られたＥの値を座標（ｕ，ｖ）のピクセル値とし、（６０３）に戻る。
【００２９】
以上のようにして（６０３）から（６１１）の一連の処理を繰り返し、投影面のＷ＊Ｈ個のピクセル座標（ｕ，ｖ）すべてについて反射光Ｅ（ｕ，ｖ）を計算することにより、始点Ｏから観察した物体内部の３次元構造を透視するイメージ画像を投影面上に描画することができる。
【００３０】
【発明が解決しようとする課題】
レイキャスティングによりボリュームレンダリングを行う３次元画像表示は、ボクセルデータから直接物体構造を描画することができるため、特に、人体において、骨、内臓等の組織が複雑に入り組んでいる場合にあっても、これらをほぼ明確に分離描画することができ、人体の内部組織の空間的な位置関係を直観的に効率よく把握できるという利点がある。
【００３１】
一方で、ボリュームレンダリングには、３次元空間に連続的に存在する物理量を離散的に測定して得られた有限個のボクセルデータで表現することに基づく原理的な精度の問題がある。精度に影響を与える要因として、測定ノイズやサンプリング間隔等のデータ取得に起因するものと、そのデータに対する計算法に起因するものがある。前者に対しては、空間フィルタリング等により画像のコントラストを改善する方法が提案されている。
【００３２】
計算法に関しては、レイキャスティングではグラディエントから得られるシェーディング係数が３次元画像の画質に直接影響するために、微妙な３次元形状が抽出できるようなグラディエントの計算方法がシェーディング手法として研究されている。前述した格子点（Ｘ，Ｙ，Ｚ）のグラディエントＧ（Ｘ，Ｙ，Ｘ）の算出法は、優れた立体感のある画像が得られるとされるグレーレベルグラディエント法によるものであり、図６の説明では最低限の６近傍の例を示したが、例えば２６近傍のボクセルデータを使用することにより精度を上げることができる。
【００３３】
しかしながら、図６の処理フローチャートに示した従来の計算法においては、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）のシェーディング係数β（Ｒｘ，Ｒｙ，Ｒｚ）として、レイ到達点Ｒの周囲の格子点（離散的な位置）におけるグラディエントを用いて求めたそれぞれの格子点のシェーディング係数βｍから補間計算により算出した補間シェーディング係数を使用するため、シェーディング係数βｍの計算点とレイ到達点Ｒとの間にずれがあり、これが画像に好ましくないアーティファクトを生ずる等の３次元画像の画質劣化を招く原因となる。また、従来の計算法においては、格子点のグラディエントを予め計算してＨＤＤまたはメモリに格納しているため、画像データの読み込み時にグラディエントを事前計算するために時間がかかるとともに、常に一定量の記憶領域が占有されるという問題もある。
【００３４】
本発明は上記事情に鑑みてなされたもので、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）に計算中心を置いたグラディエントの算出を可能にする計算法を提供することにより、レイ到達点Ｒのシェーディング係数β（Ｒｘ，Ｒｙ，Ｒｚ）を、レイ到達点Ｒの周囲の格子点（離散的な位置）におけるグラディエントに基づいて算出するのではなく、前記レイ到達点Ｒに計算中心を置いたグラディエントから算出することを可能にし、レイキャスティングにより描画される３次元画像の画質を向上させる３次元画像表示方法、装置およびプログラムを提供することを目的とする。
【００３５】
【課題を解決するための手段】
本発明の請求項１に係る３次元画像表示方法は、前述したレイキャスティングの手法によりボリュームレンダリングを行う３次元画像表示方法において、レイ到達点における複数個の近傍座標点を決める座標差分を任意に指定して前記近傍座標点を決定するステップと、前記近傍座標点それぞれの周囲のボクセルのボクセル値から前記近傍座標点それぞれの補間ボクセル値を算出するステップと、前記補間ボクセル値を用いてレイ到達点におけるグラディエントを算出するステップとを含むものである。また、請求項３に係る３次元画像表示プログラムは、ボクセルデータの各ボクセル値を変換して得られる各ボクセルの不透明度と前記ボクセルデータのグラディエントから得られるシェーディング係数とを用い、仮想始点から発する仮想光線に沿ってあらかじめ設定された間隔で刻まれたレイ到達点ごとに算出される部分反射光を積算して仮想反射光を得ることにより仮想投影面に画像を形成するレイキャスティングによりボリュームレンダリングを行う３次元画像表示プログラムにおいて、コンピュータを、前記レイ到達点における複数個の近傍座標点を決める座標差分を任意に指定して前記近傍座標点を決定する手段、前記近傍座標点それぞれの周囲のボクセルのボクセル値から前記近傍座標点それぞれの補間ボクセル値を算出する手段、前記補間ボクセル値を用いて前記レイ到達点におけるグラディエントを算出する手段、として機能させるものである。
【００３６】
請求項１および請求項４に係る発明によれば、レイ到達点におけるグラディエントの算出に必要な複数個の近傍座標点のボクセル値を求め、次に、その近傍座標点ボクセル値からレイ到達点におけるグラディエントを算出することで、レイ到達点のシェーディング係数をレイ到達点の周囲の格子点（離散的な位置）におけるグラディエントに基づいて算出する従来の計算法に比べて、レイ到達点とグラディエント計算中心が一致しており、レイキャスティングにより描画される３次元画像の画質を向上させることができる。また、格子点でのグラディエントを事前に計算して保持しておく必要がないため、メモリ資源を有効に活用することができる。
【００３７】
本発明の請求項２に係る３次元画像表示方法は、請求項１に記載の３次元画像表示方法において、レイ到達点の座標に対する前記近傍座標点の座標差分を任意に指定することを可能にするものである。また、請求項５に係る３次元画像表示プログラムは、請求項４に記載の３次元画像表示プログラムにおいて、前記レイ到達点の座標に対する前記近傍座標点の座標差分が任意に指定されるものである。
【００３８】
請求項２および請求項５に係る発明によれば、レイ到達点の近傍座標点を決定する座標差分の任意指定を可能にすることにより、ボクセルデータの性質に応じた指定や描画結果を観測しながらの可変的な指定が可能になり、描画する対象に応じて近傍座標点を最適に決定することができ、描画される３次元画像の画質をさらに向上させることができる。
【００３９】
【発明の実施の形態】
以下、本発明の実施の形態について図面を参照して説明する。図１は本発明の一実施の形態に係るレイキャスティングによる３次元画像表示方法を示す処理フローチャートである。図１の説明において、処理フローチャートの各ステップを（ステップの符号）で表記する。
【００４０】
図１において、（１０１）は、使用するボクセルデータを格子点（Ｘ，Ｙ，Ｚ）についてＶ（Ｘ，Ｙ，Ｚ）と表現することを示す。ここで、図６に示した従来の処理フローチャートの場合と異なり、各格子点（Ｘ，Ｙ，Ｚ）のグラディエントＧ（Ｘ，Ｙ，Ｘ）の算出は行なわない。
【００４１】
（１０２）では、投影面の各座標（ｕ，ｖ）への反射光Ｅ（ｕ，ｖ）を求めるために、投影面のピクセルサイズをＷ＊Ｈとして、以降の処理フローに対して、０≦ｕ＜Ｗ、０≦ｖ＜Ｈの繰り返し制御を行う。
【００４２】
（１０３）では、始点座標ＯおよびステップベクトルΔＳ＝（ΔＳｘ，ΔＳｙ，ΔＳｚ）を決定する。始点Ｏは、平行投影法では投影面のピクセル座標（ｕ，ｖ）にあり、透視投影法では空間に任意に設定した固定点である。ステップベクトルΔＳはレイ到達点を刻む間隔を長さに持つベクトルであり、平行投影法では仮想光線方向に一定の向きを持ち、透視投影法では仮想光線方向に応じて向きが変化する。
【００４３】
（１０４）では、投影面の座標（ｕ，ｖ）の現在の値に対して反射光Ｅ（ｕ，ｖ）を求める計算の繰り返し制御を行うために、レイ到達点を始点Ｏに設定し、反射光と残存光について、Ｅ＝０、Ｉ＝Ｉ（０）に初期化する。
【００４４】
（１０５）では、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）の補間ボクセル値Ｖ（Ｒｘ，Ｒｙ，Ｒｚ）を、レイ到達点Ｒの周囲のＭ個のボクセルデータから算出する。通常はレイ到達点Ｒを囲む立方格子点のボクセルデータから、ｘ方向、ｙ方向、ｚ方向それぞれの補間計算を行う。本実施形態では、立法格子点８個のボクセルデータから補間計算を行うものとする。求めた補間ボクセル値から、あらかじめ設定した変換関数により、レイ到達点Ｒの不透明度α（Ｒｘ，Ｒｙ，Ｒｚ）を得る。
【００４５】
前記補間計算は、例えば、レイ到達点Ｒを含むｚ方向に垂直な平面と８個のボクセルデータの格子点からなる立方体との交差線が描く正方形の４個の頂点の補間ボクセル値を、前記８個のボクセルデータからｚ方向の補間により求め、次に、レイ到達点Ｒを含むｘ方向の直線と前記正方形との交点の２点の補間ボクセル値を、前記４個の正方形の頂点の補間ボクセル値からｙ方向の補間により求め、最後に目的の補間ボクセル値を、この２点の補間ボクセル値からｘ方向の補間により求めることができる。
【００４６】
（１０６）では、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）のグラディエントを求めるために、グラディエントを算出するのに必要なＮ個の近傍座標点ＲＡｎ（ＲＡｘｎ，ＲＡｙｎ，ＲＡｚｎ）を決定する。近傍座標点はｘ方向、ｙ方向、ｚ方向それぞれに必要なので、最小限６近傍とするが、目的に応じて多近傍を選ぶことができ、例えば２６近傍等とすることもできる。
【００４７】
（１０７）では、各近傍座標点ＲＡｎ（ＲＡｘｎ，ＲＡｙｎ，ＲＡｚｎ）の補間ボクセル値ＶＡｎ（ＲＡｘｎ，ＲＡｙｎ，ＲＡｚｎ）を、各近傍座標点ＲＡｎの周囲のそれぞれＰ個のボクセルデータから算出する。通常は各近傍座標点ＲＡｎそれぞれを囲む立方格子点８個のボクセルデータから、ｘ方向、ｙ方向、ｚ方向それぞれの補間計算を行う。
【００４８】
（１０８）では、Ｎ個の近傍座標点ＲＡｎの補間ボクセル値ＶＡｎから、レイ到達点ＲのグラディエントＧ（Ｒｘ，Ｒｙ，Ｒｚ）を算出する。得られたグラディエントと始点Ｏからレイ到達点Ｒに至るレイベクトルとのなす角度θからレイ到達点Ｒのシェーディング係数β（Ｒｘ，Ｒｙ，Ｒｚ）＝ＡＢＳ（ｃｏｓθ）を得る。
【００４９】
（１０９）では、不透明度αとシェーディング係数βから、レイ到達点Ｒにおける減衰光Ｄおよび部分反射光Ｆを求める。これらを用いて、（１１０）では、残存光Ｉおよび反射光Ｅを更新し、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）をステップベクトルΔＳだけ進める。（１１１）では、繰り返し制御のための終端条件をチェックする。終端条件が満たされると、（１１２）で、以上の計算から得られたＥの値を座標（ｕ，ｖ）のピクセル値とし、（１０２）に戻る。
【００５０】
以上のようにして（１０２）から（１１２）の一連の処理を繰り返し、投影面のＷ＊Ｈ個のピクセル座標（ｕ，ｖ）すべてについて反射光Ｅ（ｕ，ｖ）を計算することにより、始点Ｏから観察した物体内部の３次元構造を透視するイメージ画像を投影面上に描画することができる。
【００５１】
図２はレイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）におけるシェーディング係数β（Ｒｘ，Ｒｙ，Ｒｚ）の算出方法について、従来のレイキャスティングによる方法と、本発明に係るレイキャスティングによる方法とを比較して説明する図である。図２において、２１０はレイ到達点Ｒであり、２２０は仮想光線を示し、図は便宜的に２次元で表現している。以下の説明も便宜的に２次元で説明する。
【００５２】
図２の（ａ）は、従来のレイキャスティングによるシェーディング係数の算出方法を示すものである。レイ到達点Ｒの周囲の格子点Ｒ１、Ｒ２、Ｒ３、Ｒ４におけるそれぞれのグラディエントと仮想光線２２０のなす角度から各格子点におけるシェーディング係数β１、β２、β３、β４を求める。ここで用いる格子点Ｒ１、Ｒ２、Ｒ３、Ｒ４のグラディエントは、図６の（６０１）で全格子点について算出しておいたものである。次に、各格子点のシェーディング係数β１、β２、β３、β４から補間計算によりレイ到達点Ｒの補間シェーディング係数βを算出する。
【００５３】
図２の（ｂ）は、本発明に係るレイキャスティングによるグラディエント算出方法を示すものである。図２の（ｂ）において、ＲＡ１、ＲＡ２、ＲＡ３、ＲＡ４はレイ到達点Ｒの近傍座標点であり、Ｖ１１、Ｖ１２、Ｖ１３、Ｖ１４は近傍座標点の１つであるＲＡ１の周囲の格子点におけるボクセルデータを示している。
【００５４】
グラディエントを算出するには、まず、それぞれの近傍座標点ＲＡｎの補間ボクセル値ＶＡｎを、それぞれの近傍座標点ＲＡｎを囲むボクセルデータから補間計算により求める。例えば、近傍座標点ＲＡ１の補間ボクセル値を、近傍座標点ＲＡ１を囲むボクセルデータＶ１１、Ｖ１２、Ｖ１３、Ｖ１４から補間計算する。次に、近傍座標点の補間ボクセル値ＶＡｎからレイ到達点ＲのグラディエントＧを求める。図２の（ｂ）の例では、（ＶＡ１−ＶＡ３）と（ＶＡ２−ＶＡ４）を直交成分に持つベクトルがグラディエントを与える。
【００５５】
このようにして、レイ到達点Ｒ（Ｒｘ，Ｒｙ，Ｒｚ）のグラディエントＧ（Ｒｘ，Ｒｙ，Ｒｚ）を近傍座標点の補間ボクセル値ＶＡｎから算出することで、レイ到達点Ｒとグラディエントの計算中心が一致するため、このグラディエントを用いてシェーディング係数β（Ｒｘ，Ｒｙ，Ｒｚ）を求めることにより、レイ到達点Ｒの周囲の格子点のグラディエント基づいて補間計算により補間シェーディング係数を求める従来の方法に比べて、レイキャスティングにより描画される３次元画像の画質を向上させることができる。
【００５６】
図３は、本発明に係るレイキャスティングによるグラディエント算出方法におけるレイ到達点の近傍座標点の決定方法を説明する図である。図３において、３１０はレイ到達点Ｒであり、３２０は仮想光線を示し、図は便宜的に２次元で表現している。
【００５７】
近傍座標点ＲＡｎのレイ到達点Ｒからの距離、すなわちｘ方向、ｙ方向、ｚ方向のいずれかの方向の座標差分は、図３の（ａ）に示すように、ボクセルデータの格子間距離、すなわちボクセル間の単位距離にするのが最も簡明である。これに対して、図３の（ｂ）の例では、近傍座標点ＲＡｎのレイ到達点Ｒからの距離をボクセルデータの格子間距離の１／２にしている。図示例では、いずれも差分計算の方向をボクセルの配列方向と同一にしているが、必ずしもボクセルの配列方向と同一にする必要はない。すなわち、任意に指定可能な座標差分とは、レイ到達点Ｒからの距離および方向を考慮して決定される。
【００５８】
この近傍座標点の決め方は描画される３次元画像の画質に影響を与えるが、その最適値はボクセルデータの性質によって異なる。本発明の請求項２に係るレイキャスティングによるグラディエント算出方法は、このレイ到達点の近傍座標点を決める座標差分を任意に指定することを可能にするものである。近傍座標点の位置をボクセルデータの性質に応じて指定することにより、あるいは描画結果を観測しながら可変的に指定することにより、描画する対象に応じて近傍座標点を最適に決定することができる。
【００５９】
【発明の効果】
以上説明したように、本発明によれば、レイ到達点のシェーディング係数をレイ到達点の周囲の格子点（離散的な位置）におけるグラディエントに基づいて算出せずに、レイ到達点におけるグラディエントの算出に必要な複数個の近傍座標点のボクセル値を求め、その近傍座標点ボクセル値からグラディエントを算出することでレイ到達点と計算中心が一致するグラディエントが得られ、このグラディエントを用いてシェーディング係数を求めることにより、レイキャスティングにより描画される３次元画像の画質を向上させることができる。さらに、格子点でのグラディエントを事前に計算して保持しておく必要がないため、メモリ資源を有効に活用することができる。
【００６０】
さらに本発明によれば、レイ到達点の近傍座標点を決定する座標差分の任意指定を可能にすることにより、ボクセルデータの性質に応じた指定や描画結果を観測しながらの可変的な指定が可能になり、描画する対象に応じて近傍座標点を最適に決定することができ、描画される３次元画像の画質をさらに向上させることができる。
【図面の簡単な説明】
【図１】本発明の一実施の形態に係るレイキャスティングによる３次元画像表示方法を示す処理フローチャートである。
【図２】レイ到達点におけるシェーディング係数の算出方法について、従来の方法と本発明に係る方法とを比較して説明する図である。
【図３】本発明に一実施の形態に係るレイキャスティングによる３次元画像表示方法において、グラディエントを算出するためのレイ到達点の近傍座標点の決定方法を説明する図である。
【図４】レイキャスティングによる投影画像形成の概念を説明する図である。
【図５】透視投影法によるレイキャスティングの概念を説明する図である。
【図６】従来のレイキャスティングによる３次元画像表示方法を示す処理フローチャートである。
【符号の説明】
２１０、３１０レイ到達点
２２０、３２０仮想光線
４１０、５３０仮想投影面
４２０、５２０ボクセルデータで表現される３次元物体
４２１、４２２、４２３レイ到達点に位置するボクセル
５１０始点
５１１、５３１仮想反射光
５１２仮想網膜
５２１レイベクトル
５２２ステップベクトル[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a three-dimensional image display method for performing volume rendering by ray casting., Device and programIn particular, the present invention relates to a method for calculating a gradient of voxel data used in ray casting.
[0002]
[Prior art]
The advent of CT (Computed Tomography), which has made it possible to directly observe the internal structure of the human body with the progress of image processing technology using a computer, is a technology that has brought innovation to the medical field, and medical treatment using tomographic images of a living body Diagnosis is widespread. Furthermore, in recent years, as a technique for visualizing a three-dimensional structure inside an object that is difficult to understand only with a tomographic image, volume rendering for directly drawing an image of the three-dimensional structure from three-dimensional digital data of the object obtained by CT has attracted attention.
[0003]
Ray casting is known as an excellent technique for volume rendering. In ray casting, a virtual ray (ray) is irradiated to an object from a virtual starting point, and an image of virtual reflected light from the inside of the object is formed on a virtual projection plane, so that a three-dimensional structure inside the object is seen through. Is a method of forming The basic theory of ray casting is described, for example, in "New Generation Three-Dimensional CT Diagnosis" (November 1, 1995, published by Nankodo Co., Ltd.).
[0004]
The main points of ray casting will be described below. A small unit area that is a constituent unit of the three-dimensional area of the object is called a voxel, and unique data representing characteristics such as a density value of the voxel is called a voxel value. The entire object is represented by voxel data, which is a three-dimensional array of voxel values. Usually, two-dimensional tomographic image data obtained by CT is stacked in a direction perpendicular to the tomographic plane, and necessary interpolation is performed to obtain three-dimensional voxel data.
[0005]
It is assumed that the virtual reflected light corresponding to the virtual light beam emitted from the virtual start point to the object is generated according to the opacity artificially set for the voxel value. Further, a gradient of the voxel data, that is, a normal vector is obtained to capture a virtual surface, and a shading coefficient for shading is calculated from a cosine of an angle formed by the virtual ray and the normal vector. The virtual reflected light is calculated by multiplying the intensity of the virtual light beam applied to the voxel by the opacity of the voxel and the shading coefficient.
[0006]
FIG. 4 is a diagram for explaining the concept of projection image formation by ray casting. In FIG. 4, reference numeral 410 denotes a virtual projection plane, 420 denotes a three-dimensional object represented by voxel data, and 421, 422, and 423 denote three-dimensional objects 420 from a coordinate point O (u, v) on the virtual projection plane 410. Are the voxels V (n-1), Vn, V (n + 1) located at the ray arrival points carved at regular intervals in order to calculate the virtual reflected light along the virtual light rays irradiated on the virtual rays. In general, the direction of the virtual ray and the interval at which the ray arrival point is carved do not match the axial direction of the three-dimensional array of voxel data or the voxel interval. Therefore, it is necessary to obtain the voxel characteristic value at each ray arrival point by interpolation calculation.
[0007]
With respect to the voxel Vn located at the n-th ray arrival point on the virtual ray, the opacity given to the voxel value is αn, and the shading coefficient obtained from the gradient of the voxel data is βn. Here, assuming that the virtual light beam incident on the voxel Vn is attenuated by the opacity of the voxel Vn, the amount of attenuation is defined as attenuated light, and the light attenuated by the attenuated light is defined as residual light. The residual light that has passed through the voxel Vn becomes a virtual ray incident on the next voxel V (n + 1). It is also assumed that the attenuated light gives the reflected light of the voxel Vn. If the remaining light from voxel V (n-1) is I (n-1), the remaining light of voxel Vn is In, the attenuated light of voxel Vn is Dn, and the partially reflected light by voxel Vn is Fn, a multiplication symbol * The following equation is obtained.
[0008]
Dn = αn * I (n-1)
Fn = βn * Dn = βn * αn * I (n-1)
In = I (n-1) -Dn = (1-αn) * I (n-1)
[0009]
The virtual reflected light E (u, v) projected on the coordinates (u, v) on the virtual projection plane 410 is obtained by integrating the partial reflected light calculated for each ray arrival point. The purpose of this method is to form a virtual perspective image of the three-dimensional structure inside the object, and it is not a calculation relating to physical phenomena in the real world. Therefore, attenuation of light reflected from voxels is not considered. The intensity of a virtual ray emitted from O (u, v) is defined as I (0). As a notation, ； (i = 1; n) Ai is a sum symbol of n terms from the sequence A1 to An, and Π (i = 1; n) If Ai is the product symbol of n terms from the sequences A1 to An, the following equation is obtained.
E (u, v) = Σ (i = 1; terminal condition) Fi = I (0) * (Σ (i = 1; terminal condition) βi * αi * (Π (j = 1; i−1) (1 -Αj)))
[0010]
Here, the termination condition is when the virtual ray passes through the object or the residual light becomes zero. In this way, by calculating E (u, v) for all coordinates (u, v) on the virtual projection plane 410, a virtual perspective image of a three-dimensional image is formed.
[0011]
This method is called a parallel projection method because a virtual ray emitted from a coordinate point on the virtual projection plane becomes a parallel ray perpendicular to the virtual projection plane. This method is easy to intuitively understand because a fluoroscopic image of the inside of the object viewed from the outside is obtained. However, for example, in the medical field, there is a disadvantage that the lumen surface of the organ is observed from the outside. On the other hand, there is a perspective projection method in which a virtual image is formed by irradiating a virtual ray radially from an arbitrary start point in space. For example, by placing the starting point inside a blood vessel of a human body, it is possible to display the blood vessel lumen surface as a surface.
[0012]
FIG. 5 is a diagram for explaining the concept of ray casting by the perspective projection method. In FIG. 5, reference numeral 510 denotes an arbitrary start point O (Ox, Oy, Oz) placed in a space; 520, a three-dimensional object represented by voxel data; and 530, a virtual projection plane. Ray arrival points for calculating virtual reflected light along the virtual light rays emitted from the starting point 510 are carved at regular intervals. 521 is a ray vector from the starting point O (Ox, Oy, Oz) to the ray arrival point R (Rx, Ry, Rz), 522 is a step vector ΔS = (ΔSx, ΔSy, ΔSz) representing an interval at which the ray arrival point is carved. Reference numeral 531 denotes virtual reflected light E (u, v) projected on the coordinates (u, v) on the virtual projection plane 530. Here, Rx, Ry, Rz, ΔSx, ΔSy, and ΔSz are functions of (u, v).
[0013]
In FIG. 5, the virtual projection plane 530 is in the arrival direction of the virtual ray. However, when the starting point 510 is considered to be the lens of the virtual eyeball, the virtual reflected light 511 is at the coordinates (u, v) on the retina 512 of this eyeball. It can be considered that the virtual reflected light E (u, v) on the virtual projection plane 530 is projected as' (u, v) and enlarged and displayed symmetrically with respect to the start point 510.
[0014]
Here, for the voxel Vn (Rxn, Ryn, Rzn) located at the n-th ray arrival point Rn on the virtual ray, the opacity given to the voxel value is α (Rxn, Ryn, Rzn), and the voxel The shading coefficient obtained from the data gradient is β (Rxn, Ryn, Rzn), the residual light from voxel V (n−1) is I (n−1), the residual light of voxel Vn is In, and the attenuation of voxel Vn is In. Assuming that the light is Dn and the partially reflected light by the voxel Vn is Fn, the following equation is obtained as in the parallel projection method.
[0015]
Dn = α (Rxn, Ryn, Rzn) * I (n−1)
Fn = β (Rxn, Ryn, Rzn) * Dn = β (Rxn, Ryn, Rzn) * α (Rxn, Ryn, Rzn) * I (n−1)
In = I (n−1) −Dn = (1−α (Rxn, Ryn, Rzn)) * I (n−1)
[0016]
Similarly to the parallel projection method, when the virtual reflected light E (u, v) projected on the coordinates (u, v) on the virtual projection plane 530 is obtained, the following expression is obtained using the above-described notation.
[0017]
E (u, v) = Σ (i = 1; terminal condition) Fi = I (0) * (Σ (i = 1; terminal condition) β (Rxi, Ryi, Rzi) * α (Rxi, Ryi, Rzi) * (Π (j = 1; i-1) (1-α (Rxj, Ryj, Rzj))))
[0018]
The termination condition is when the virtual ray passes through the object or the residual light becomes zero. In this way, by calculating E (u, v) for all coordinates (u, v) on the virtual projection plane 530, a virtual perspective image of a three-dimensional image is formed.
[0019]
FIG. 6 shows a processing flowchart of a conventional three-dimensional image display method by ray casting performed based on the principle as described above. In the following description, the projection method will be described with reference to FIG. 5. However, the parallel projection method and the perspective projection method will be generalized without distinction. In addition, in the description of FIG. 6, each step of the processing flowchart is denoted by (step reference).
[0020]
In FIG. 6, (601) indicates that the voxel data to be used is expressed as V (X, Y, Z) with respect to the grid point (X, Y, Z). In (602), the gradient G (X, Y, X) = (Gx, Gy, Gz) of each grid point (X, Y, Z) is calculated from voxel data near the grid point (X, Y, Z). Calculate and store in memory. For example, if it is near 6, the vector component (component ratio) of the gradient G (X, Y, Z) is represented by the following equation.
[0021]
Gx (X, Y, X) = V ((X + 1), Y, Z) -V ((X-1), Y, Z) Gy (X, Y, X) = V (X, (Y + 1), Z ) -V (X, (Y-1), Z) Gz (X, Y, X) = V (X, Y, (Z + 1))-V (X, Y, (Z-1))
[0022]
In (603), the pixel size of the projection plane is set to W * H in order to obtain the reflected light E (u, v) for each coordinate (u, v) of the projection plane, and 0 is set for the subsequent processing flow. ≦ u <W and 0 ≦ v <H are repeatedly controlled.
[0023]
In (604), the starting point coordinates O and the step vector ΔS = (ΔSx, ΔSy, ΔSz) are determined. The starting point O is located at pixel coordinates (u, v) on the projection plane in the parallel projection method, and is a fixed point arbitrarily set in space in the perspective projection method. The step vector ΔS is a vector having a length at which the interval between the ray arrival points is ticked. The parallel vector has a fixed direction in the ray direction in the parallel projection method, and the direction changes according to the virtual ray direction in the perspective projection method.
[0024]
In (605), the ray reaching point is set to the starting point O in order to perform the repetitive control of the calculation for obtaining the reflected light E (u, v) for the current value of the coordinates (u, v) of the projection plane, The reflected light and the residual light are initialized to E = 0 and I = I (0).
[0025]
In (606), the interpolation voxel value V (Rx, Ry, Rz) of the ray arrival point R (Rx, Ry, Rz) is calculated from the M voxel data of the lattice points around the ray arrival point R. Normally, interpolation calculation in each of the x, y, and z directions is performed from voxel data of eight cubic lattice points surrounding the ray arrival point R. From the obtained interpolated voxel values, the opacity α (Rx, Ry, Rz) of the ray arrival point R is obtained by a preset conversion function.
[0026]
Next, the shading coefficient β (Rx, Ry, Rz) at the ray arrival point R (Rx, Ry, Rz) is calculated. In the conventional ray casting, a shading coefficient of each grid point is obtained from a gradient of grid points around the ray arrival point R, and an interpolation shading coefficient at the ray arrival point R is calculated from the shading coefficient of the grid point by interpolation calculation. I was
[0027]
That is, in (607), the M gradients Gm (m = 1, 2,..., M) of the grid points around the ray arrival point R (Rx, Ry, Rz), usually surround the ray arrival point R 8 gradients of cubic lattice points are obtained, and the shading coefficient at each lattice point is obtained from the angle θm (m = 1, 2,..., M) formed by each gradient and the ray vector from the starting point O to the ray arrival point R. βm = ABS (cos θm) is obtained. Here, ABS (cos θm) is the absolute value of cos θm. Next, an interpolation shading coefficient β (Rx, Ry, Rz) at the ray arrival point R is calculated from βm by interpolation calculation.
[0028]
In (608), the attenuation light D and the partially reflected light F at the ray arrival point R are obtained from the opacity α and the shading coefficient β. Using these, in (609), the residual light I and the reflected light E are updated, and the ray arrival point R (Rx, Ry, Rz) is advanced by the step vector ΔS. At (610), the terminal condition for repetitive control is checked. When the termination condition is satisfied, the value of E obtained from the above calculation is set as the pixel value of the coordinates (u, v) in (611), and the process returns to (603).
[0029]
As described above, the series of processing from (603) to (611) is repeated, and the reflected light E (u, v) is calculated for all the W * H pixel coordinates (u, v) on the projection plane. An image image that sees through the three-dimensional structure inside the object observed from the starting point O can be drawn on the projection plane.
[0030]
[Problems to be solved by the invention]
Since the three-dimensional image display that performs volume rendering by ray casting can directly draw an object structure from voxel data, especially in the human body, even when tissues such as bones and internal organs are complicated, These can be separated and drawn almost clearly, and there is an advantage that the spatial positional relationship of the internal tissue of the human body can be intuitively and efficiently grasped.
[0031]
On the other hand, volume rendering has a problem of fundamental accuracy based on representing a finite number of voxel data obtained by discretely measuring physical quantities continuously present in a three-dimensional space. Factors that affect accuracy include those caused by data acquisition such as measurement noise and sampling intervals, and those caused by the calculation method for the data. For the former, a method of improving the contrast of an image by spatial filtering or the like has been proposed.
[0032]
Regarding the calculation method, since the shading coefficient obtained from the gradient directly affects the image quality of a three-dimensional image in ray casting, a gradient calculation method capable of extracting a delicate three-dimensional shape has been studied as a shading method. The above-described method of calculating the gradient G (X, Y, X) of the grid points (X, Y, Z) is based on the gray level gradient method, which is considered to provide an image with an excellent three-dimensional effect. In the above description, an example of a minimum of 6 neighborhoods is shown, but the accuracy can be improved by using voxel data of 26 neighborhoods, for example.
[0033]
However, in the conventional calculation method shown in the processing flowchart of FIG. 6, the shading coefficient β (Rx, Ry, Rz) of the ray arrival point R (Rx, Ry, Rz) is set as a grid point around the ray arrival point R. Since the interpolation shading coefficient calculated by the interpolation calculation from the shading coefficient βm of each grid point obtained by using the gradient at (discrete position) is used, a point between the calculation point of the shading coefficient βm and the ray reaching point R is used. There is a shift, which causes deterioration of the image quality of the three-dimensional image such as generation of undesired artifacts in the image. Further, in the conventional calculation method, since the gradient of the grid point is calculated in advance and stored in the HDD or the memory, it takes time to pre-calculate the gradient when reading the image data, and a fixed amount of storage is always required. There is also a problem that the area is occupied.
[0034]
The present invention has been made in view of the above circumstances, and provides a calculation method capable of calculating a gradient centered on a ray arrival point R (Rx, Ry, Rz) by calculating the gradient of the ray arrival point R. Instead of calculating the shading coefficient β (Rx, Ry, Rz) based on the gradient at grid points (discrete positions) around the ray arrival point R, a gradient centered on the ray arrival point R is calculated. Three-dimensional image display method for improving the image quality of a three-dimensional image drawn by ray casting, Device and programThe purpose is to provide.
[0035]
[Means for Solving the Problems]
The three-dimensional image display method according to claim 1 of the present invention is the three-dimensional image display method for performing volume rendering by the above-described ray casting method, wherein a plurality of neighboring coordinate points at a ray arrival point are determined.Arbitrarily specify the coordinate difference to be determined and determine the neighboring coordinate pointsAnd calculating an interpolated voxel value of each of the neighboring coordinate points from voxel values of voxels around each of the neighboring coordinate points, and calculating a gradient at a ray reaching point using the interpolated voxel values. Including. Claims3The three-dimensional image display program according to the above, using the opacity of each voxel obtained by converting each voxel value of the voxel data and a shading coefficient obtained from the gradient of the voxel data, along a virtual ray emitted from a virtual start point A three-dimensional image display that performs volume rendering by ray casting in which an image is formed on a virtual projection plane by integrating partial reflected light calculated for each ray arrival point carved at a preset interval to obtain virtual reflected light. In the program, a computer is used to calculate a plurality of nearby coordinate points at the ray arrival pointArbitrarily specify the coordinate difference to be determined and determine the neighboring coordinate pointsMeans for calculating an interpolated voxel value of each of the neighboring coordinate points from voxel values of voxels around each of the neighboring coordinate points, and means for calculating a gradient at the ray reaching point using the interpolated voxel values. It is to let.
[0036]
Claim 1And claim 4ToSuch inventionAccording to, the voxel values of a plurality of neighboring coordinate points required for calculating the gradient at the ray arrival point are obtained, and then the gradient at the ray arrival point is calculated from the voxel values of the neighboring coordinate points, thereby obtaining the ray arrival point. Compared to the conventional calculation method that calculates the shading coefficient of the ray based on the gradient at the grid points (discrete positions) around the ray arrival point, the ray arrival point and the gradient calculation center match, and drawing by ray casting The image quality of the three-dimensional image to be obtained can be improved. Further, since it is not necessary to calculate and hold the gradient at the grid point in advance, the memory resources can be effectively used.
[0037]
According to a three-dimensional image display method according to a second aspect of the present invention, in the three-dimensional image display method according to the first aspect, it is possible to arbitrarily specify a coordinate difference between the coordinates of the ray arrival point and the neighboring coordinate points. Is what you do.A three-dimensional image display program according to a fifth aspect is the three-dimensional image display program according to the fourth aspect, wherein a coordinate difference between the coordinates of the ray arrival point and the neighboring coordinate points is arbitrarily specified. .
[0038]
Claim 2And claim 5ToSuch inventionAccording to the above, it is possible to specify a coordinate difference for determining a coordinate point in the vicinity of a ray arrival point, thereby enabling specification according to the properties of voxel data and variable specification while observing a drawing result. Thus, the neighboring coordinate points can be optimally determined according to the object to be rendered, and the image quality of the rendered three-dimensional image can be further improved.
[0039]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings. FIG. 1 is a flowchart illustrating a method for displaying a three-dimensional image by ray casting according to an embodiment of the present invention. In the description of FIG. 1, each step of the processing flowchart is denoted by (step reference).
[0040]
In FIG. 1, (101) indicates that voxel data to be used is expressed as V (X, Y, Z) with respect to a grid point (X, Y, Z). Here, unlike the case of the conventional processing flowchart shown in FIG. 6, the calculation of the gradient G (X, Y, X) of each grid point (X, Y, Z) is not performed.
[0041]
In (102), in order to obtain the reflected light E (u, v) at each coordinate (u, v) of the projection plane, the pixel size of the projection plane is set to W * H, and 0 is applied to the subsequent processing flow. ≦ u <W and 0 ≦ v <H are repeatedly controlled.
[0042]
In (103), the starting point coordinates O and the step vector ΔS = (ΔSx, ΔSy, ΔSz) are determined. The starting point O is located at pixel coordinates (u, v) on the projection plane in the parallel projection method, and is a fixed point arbitrarily set in space in the perspective projection method. The step vector ΔS is a vector having a length in which the interval between the ray arrival points is ticked. The parallel vector has a fixed direction in the virtual ray direction in the parallel projection method, and the direction changes according to the virtual ray direction in the perspective projection method.
[0043]
In (104), the ray arrival point is set to the start point O in order to perform the repetitive control of the calculation for obtaining the reflected light E (u, v) for the current value of the coordinates (u, v) of the projection plane, The reflected light and the residual light are initialized to E = 0 and I = I (0).
[0044]
In (105), the interpolation voxel value V (Rx, Ry, Rz) of the ray arrival point R (Rx, Ry, Rz) is calculated from M voxel data around the ray arrival point R. Normally, interpolation calculations in the x, y, and z directions are performed from voxel data of a cubic lattice point surrounding the ray arrival point R. In the present embodiment, interpolation calculation is performed from voxel data of eight cubic lattice points. From the obtained interpolated voxel values, the opacity α (Rx, Ry, Rz) of the ray arrival point R is obtained by a preset conversion function.
[0045]
The interpolation calculation, for example, the interpolation voxel values of the four vertices of a square drawn by the intersection line of a plane perpendicular to the z direction including the ray arrival point R and a cube consisting of grid points of eight voxel data, Interpolated in the z-direction from the eight voxel data, and then interpolated voxel values of two points at the intersection of the square with the straight line in the x-direction including the ray arrival point R and the vertices of the four squares From the voxel values, it can be obtained by interpolation in the y direction, and finally, the target interpolation voxel value can be obtained from these two interpolation voxel values by interpolation in the x direction.
[0046]
In (106), in order to obtain the gradient of the ray arrival point R (Rx, Ry, Rz), N neighboring coordinate points RAn (RAxn, RAyn, RAzn) necessary for calculating the gradient are determined. Neighboring coordinate points are required in each of the x, y, and z directions, so the minimum is six neighbors. However, multiple neighbors can be selected according to the purpose, for example, 26 neighbors.
[0047]
In (107), an interpolated voxel value VAn (RAxn, RAyn, RAzn) of each neighboring coordinate point RAn (RAxn, RAyn, RAzn) is calculated from P voxel data around each neighboring coordinate point RAn. Normally, interpolation calculation in each of the x, y, and z directions is performed from voxel data of eight cubic lattice points surrounding each of the neighboring coordinate points RAn.
[0048]
In (108), the gradient G (Rx, Ry, Rz) of the ray arrival point R is calculated from the interpolation voxel values VAn of the N neighboring coordinate points RAn. The shading coefficient β (Rx, Ry, Rz) = ABS (cos θ) of the ray reaching point R is obtained from the angle θ between the obtained gradient and the ray vector from the starting point O to the ray reaching point R.
[0049]
In (109), the attenuation light D and the partially reflected light F at the ray arrival point R are obtained from the opacity α and the shading coefficient β. Using these, in (110), the residual light I and the reflected light E are updated, and the ray arrival point R (Rx, Ry, Rz) is advanced by the step vector ΔS. At (111), the terminal condition for repetitive control is checked. When the termination condition is satisfied, the value of E obtained from the above calculation is set as the pixel value of the coordinates (u, v) in (112), and the process returns to (102).
[0050]
As described above, the series of processing from (102) to (112) is repeated to calculate the reflected light E (u, v) for all the W * H pixel coordinates (u, v) on the projection plane, An image image that sees through the three-dimensional structure inside the object observed from the starting point O can be drawn on the projection plane.
[0051]
FIG. 2 shows a method of calculating a shading coefficient β (Rx, Ry, Rz) at a ray arrival point R (Rx, Ry, Rz) by comparing a conventional ray casting method with a ray casting method according to the present invention. FIG. In FIG. 2, 210 is a ray reaching point R, 220 is a virtual ray, and the figure is expressed in two dimensions for convenience. The following description is also made in two dimensions for convenience.
[0052]
FIG. 2A shows a method of calculating a shading coefficient by conventional ray casting. Shading coefficients β1, β2, β3, and β4 at each grid point are obtained from angles formed by the virtual rays 220 and the respective gradients at grid points R1, R2, R3, and R4 around the ray arrival point R. The gradients of the grid points R1, R2, R3, R4 used here are calculated for all grid points in (601) of FIG. Next, the interpolation shading coefficient β of the ray reaching point R is calculated by interpolation calculation from the shading coefficients β1, β2, β3, and β4 of each grid point.
[0053]
FIG. 2B shows a gradient calculation method by ray casting according to the present invention. In FIG. 2B, RA1, RA2, RA3, and RA4 are coordinate points near the ray arrival point R, and V11, V12, V13, and V14 are grid points around RA1, which is one of the nearby coordinate points. 14 shows voxel data.
[0054]
In order to calculate the gradient, first, an interpolation voxel value VAn of each neighboring coordinate point RAn is obtained by an interpolation calculation from voxel data surrounding each neighboring coordinate point RAn. For example, the interpolation voxel value of the neighboring coordinate point RA1 is interpolated from the voxel data V11, V12, V13, and V14 surrounding the neighboring coordinate point RA1. Next, a gradient G of the ray reaching point R is obtained from the interpolation voxel value VAn of the neighboring coordinate point. In the example of FIG. 2B, a vector having (VA1-VA3) and (VA2-VA4) as orthogonal components gives a gradient.
[0055]
By calculating the gradient G (Rx, Ry, Rz) of the ray arrival point R (Rx, Ry, Rz) from the interpolated voxel value VAn of the neighboring coordinate point in this way, the calculation center of the ray arrival point R and the gradient is calculated. Since the shading coefficient β (Rx, Ry, Rz) is obtained using this gradient, the conventional method of obtaining the interpolation shading coefficient by interpolation calculation based on the gradient of the grid points around the ray arrival point R is obtained. In comparison, the image quality of a three-dimensional image drawn by ray casting can be improved.
[0056]
FIG. 3 is a diagram illustrating a method of determining coordinate points near a ray arrival point in the gradient calculation method by ray casting according to the present invention. In FIG. 3, reference numeral 310 denotes a ray reaching point R, 320 denotes a virtual ray, and the figure is expressed in two dimensions for convenience.
[0057]
As shown in FIG. 3A, the distance between the neighboring coordinate point RAn and the ray arrival point R, that is, the coordinate difference in any of the x direction, the y direction, and the z direction, is the distance between the grids of the voxel data. That is, it is easiest to set the unit distance between voxels. On the other hand, in the example of FIG. 3B, the distance of the neighboring coordinate point RAn from the ray arrival point R is set to の of the inter-grid distance of the voxel data. In the illustrated example, the direction of the difference calculation is the same as the arrangement direction of the voxels in each case, but it is not necessarily required to be the same as the arrangement direction of the voxels. That is, the coordinate difference that can be arbitrarily specified is determined in consideration of the distance and the direction from the ray arrival point R.
[0058]
The method of determining the neighboring coordinate points affects the image quality of the three-dimensional image to be drawn. A gradient calculating method by ray casting according to a second aspect of the present invention enables a coordinate difference for determining a coordinate point near a ray arrival point to be arbitrarily specified. By specifying the position of the nearby coordinate point according to the properties of the voxel data, or variably specifying the position while observing the drawing result, the nearby coordinate point can be optimally determined according to the drawing target. .
[0059]
【The invention's effect】
As described above, according to the present invention, the gradient at the ray arrival point is calculated without calculating the shading coefficient of the ray arrival point based on the gradient at grid points (discrete positions) around the ray arrival point. The voxel values of a plurality of neighboring coordinate points required for are obtained, and a gradient is calculated from the voxel values of the neighboring coordinate points to obtain a gradient in which the ray arrival point and the calculation center match, and the shading coefficient is calculated using this gradient. By obtaining, the image quality of the three-dimensional image drawn by ray casting can be improved. Further, since it is not necessary to calculate and hold the gradient at the grid point in advance, the memory resources can be effectively used.
[0060]
Furthermore, according to the present invention, by allowing arbitrary designation of a coordinate difference for determining a coordinate point in the vicinity of a ray arrival point, it is possible to perform designation according to the properties of voxel data and variable designation while observing a drawing result. This makes it possible to optimally determine the neighboring coordinate points according to the object to be drawn, and to further improve the image quality of the three-dimensional image to be drawn.
[Brief description of the drawings]
FIG. 1 is a processing flowchart showing a three-dimensional image display method by ray casting according to one embodiment of the present invention.
FIG. 2 is a diagram for explaining a method of calculating a shading coefficient at a ray arrival point by comparing a conventional method and a method according to the present invention.
FIG. 3 is a diagram illustrating a method for determining coordinate points near a ray arrival point for calculating a gradient in a three-dimensional image display method using ray casting according to an embodiment of the present invention.
FIG. 4 is a diagram illustrating the concept of projection image formation by ray casting.
FIG. 5 is a diagram illustrating the concept of ray casting by the perspective projection method.
FIG. 6 is a processing flowchart showing a conventional three-dimensional image display method by ray casting.
[Explanation of symbols]
210, 310 Ray arrival point
220, 320 virtual rays
410, 530 Virtual projection plane
420, 520 3D object represented by voxel data
421, 422, 423 Voxel located at ray reaching point
510 starting point
511,531 Virtually reflected light
512 virtual retina
521 ray vector
522 step vector

Claims

ボクセルデータの各ボクセル値を変換して得られる各ボクセルの不透明度と前記ボクセルデータのグラディエントから得られるシェーディング係数とを用いてレイキャスティングによりボリュームレンダリングを行う３次元画像表示方法において、
前記レイキャスティングは仮想始点から発する仮想光線に沿ってあらかじめ設定された間隔で刻まれたレイ到達点ごとに算出される部分反射光を積算して仮想反射光を得ることにより仮想投影面に画像を形成する手法であって、前記レイ到達点における複数個の近傍座標点を決める座標差分を任意に指定して前記近傍座標点を決定するステップと、前記近傍座標点それぞれの周囲のボクセルのボクセル値から前記近傍座標点それぞれの補間ボクセル値を算出するステップと、前記補間ボクセル値を用いて前記レイ到達点におけるグラディエントを算出するステップと、を含むことを特徴とする３次元画像表示方法。A three-dimensional image display method for performing volume rendering by ray casting using opacity of each voxel obtained by converting each voxel value of voxel data and a shading coefficient obtained from a gradient of the voxel data,
The ray casting integrates the partial reflected light calculated for each ray arrival point carved at a preset interval along a virtual ray emitted from a virtual starting point to obtain a virtual reflected light, thereby forming an image on a virtual projection surface. Arbitrarily specifying a coordinate difference for determining a plurality of nearby coordinate points at the ray arrival point to determine the nearby coordinate points, and voxel values of voxels around each of the nearby coordinate points. A step of calculating an interpolated voxel value of each of the neighboring coordinate points from the following, and a step of calculating a gradient at the ray arrival point using the interpolated voxel value.

請求項１に記載の３次元画像表示方法を実施するための３次元画像表示装置。A three-dimensional image display device for performing the three-dimensional image display method according to claim 1.

ボクセルデータの各ボクセル値を変換して得られる各ボクセルの不透明度と前記ボクセルデータのグラディエントから得られるシェーディング係数とを用い、仮想始点から発する仮想光線に沿ってあらかじめ設定された間隔で刻まれたレイ到達点ごとに算出される部分反射光を積算して仮想反射光を得ることにより仮想投影面に画像を形成するレイキャスティングによりボリュームレンダリングを行う３次元画像表示プログラムにおいて、コンピュータを、Using the opacity of each voxel obtained by converting each voxel value of the voxel data and the shading coefficient obtained from the gradient of the voxel data, it is engraved at a preset interval along a virtual ray emitted from a virtual start point In a three-dimensional image display program for performing volume rendering by ray casting in which an image is formed on a virtual projection plane by integrating the partial reflected light calculated for each ray arrival point to obtain a virtual reflected light,
前記レイ到達点における複数個の近傍座標点を決める座標差分を任意に指定して前記近傍座標点を決定する手段、前記近傍座標点それぞれの周囲のボクセルのボクセル値から前記近傍座標点それぞれの補間ボクセル値を算出する手段、前記補間ボクセル値を用いて前記レイ到達点におけるグラディエントを算出する手段、として機能させることを特徴とする３次元画像表示プログラム。 Means for arbitrarily specifying a coordinate difference for determining a plurality of neighboring coordinate points at the ray arrival point to determine the neighboring coordinate points, interpolation of each of the neighboring coordinate points from voxel values of voxels around each of the neighboring coordinate points A three-dimensional image display program functioning as means for calculating a voxel value and means for calculating a gradient at the ray reaching point using the interpolated voxel values.