JP4105176B2 - 画像処理方法および画像処理プログラム - Google Patents

Description

本発明は、ボリュームレンダリングによる画像処理方法および画像処理プログラムに関する。

コンピュータを用いた画像処理技術の進展により人体の内部構造を直接観測することを可能にしたＣＴ（Computed Tomography）装置、ＭＲＩ（Magnetic Resonance Imaging）装置の出現は医療分野に革新をもたらし、生体の断層画像を用いた医療診断が広く行われている。さらに、近年では、断層画像だけでは分かりにくい複雑な人体内部の３次元構造を可視化する技術として、例えば、ＣＴ装置により得られた物体の３次元デジタルデータから３次元構造のイメージを直接描画するボリュームレンダリングが医療診断に使用されている。

また、ＣＴ装置から得られる画像データを扱う場合に、観察対象にＲＯＩ（関心領域：region of interest）を設定することが行われる。２次元的な領域を表す領域ＲＯＩは、通常、臨床的に関心のある部位を指定するために使用され、ＣＴ画像や核医学画像の上に正方形、円、任意形状の領域として設定される。ＲＯＩが設定されると、そのＲＯＩ内の画素の濃度（画素値）のヒストグラム分布、平均値、標準偏差、面積などが演算装置によって求められＣＴ画像とともに表示される。医師は、これらの画像を見ながら診断を行っている。

図１７は、３次元医療画像で用いられるヒストグラムの例を示す。３次元医療画像のヒストグラムは、ボクセル値（ＣＴ値など）ごとに、ボクセル数を集計したものである。ボクセル値は組織によって異なるため、ヒストグラムは物質の組成を表す。このため、ヒストグラムから、特定の組織の存在を確認したり、その体積や体積比を計算したりすることができる。

図１８は、従来の画像処理方法において、ヒストグラムを計算するフローチャートを示す。従来の画像処理方法では、まず、メモリ上でヒストグラム領域として、ボクセル値の取りうる範囲のVmin〜Vmaxの配列である頻度freqを確保する（ステップＳ１０１）。

次に、頻度freq[Vmin〜Vmax]を0で初期化し（ステップＳ１０２）、ボリュームデータVol(x,y,z)を構成する配列のx,y,z各要素の要素数x_max, y_max, z_max を取得する（ステップＳ１０３）。

次に、ボリュームデータを構成する配列のx,y,z要素の反復子i, j, k を確保する（ステップＳ１０４）。ボリュームを走査する３重ループの開始としてk = 0に設定し（ステップＳ１０５）、k < z_maxかどうか判断する（ステップＳ１０６）。

そして、k < z_max の場合（Ｙｅｓ）は、j = 0に設定し（ステップＳ１０７）、j < y_maxかどうか判断する（ステップＳ１０８）。そして、j < y_maxの場合（Ｙｅｓ）は、i = 0に設定し（ステップＳ１０９）、i < x_maxかどうか判断する（ステップＳ１１０）。

そして、i < x_max の場合（Ｙｅｓ）は、vox = Vol(i, j, k)により、ボリュームデータ内の各位置でのボクセル値を取得する（ステップＳ１１１）。次に、freq(vox) = freq(vox) +1により、ボクセル値に対応するヒストグラム値のカウントを１加算する（ステップＳ１１２）。

そして、i = i+1を計算し（ステップＳ１１３）、ステップＳ１１０に戻る。ステップＳ１１０において、i < x_maxでない場合（Ｎｏ）は、j = j+1を計算し（ステップＳ１１４）、ステップＳ１０８に戻る。ステップＳ１０８において、j < y_maxでない場合（Ｎｏ）は、k = k+1を計算し（ステップＳ１１５）、ステップＳ１０６に戻る。ステップＳ１０６において、k < z_maxでない場合（Ｎｏ）は、終了する（ステップＳ１１６）。

なお、２次元画像において、その画像データのヒストグラムを計算し、更にそのヒストグラムをユーザに操作させることによって画像データを変更するものがある（例えば、非特許文献１参照）。

一方、３Ｄ画像の関連技術として、３Ｄ画像上でＲＯＩを設定し、ＲＯＩ内の画素値の統計量を求め、ヒストグラムを表示するものがある（例えば、特許文献１参照）。また、ボリュームデータ全体のヒストグラムをもとに、カラーＬＵＴ(Look Up Table)関数を設定するものがある(例えば、特許文献２参照)。

特開平10−21362号公報 米国特許第6658080号明細書 Adobe社PhotoShopElements2.0取扱説明書

しかしながら、上記従来の画像処理方法にあっては、ヒストグラムを計算する対象となるボクセル群は、例えば、表示されている画像データ全体や３Ｄ画像上で設定したＲＯＩ内の画素値等、通常固定されている。このため、医師等が表示角度や拡大率を変化させてもヒストグラムは変化せず、必要とする情報が得られない場合がある。

本発明は、上記従来の事情に鑑みてなされたものであって、画像表示とヒストグラム表示とにより医療診断を行う場合に、ユーザの意図に応じて適切なヒストグラムを表示させることができる画像処理方法および画像処理プログラムを提供することを目的としている。

本発明の画像処理方法は、ボリュームレンダリングによる画像処理方法であって、第１のボリュームデータを用いてレンダリングによって生成した２次元画像に設定した領域に対応する前記第１のボリュームデータの領域に含まれるボクセル群のヒストグラムを動的に生成する。

上記構成によれば、例えばユーザが表示を変更した場合など、領域の設定を変更する度にヒストグラムが再計算される為、画像表示とヒストグラム表示とにより医療診断を行う場合に、ユーザの意図に応じて適切なヒストグラムを表示させることができる。

また、本発明の画像処理方法は、ボリュームレンダリングによる画像処理方法であって、第１のボリュームデータを用いてレンダリングによって生成した２次元画像に設定した領域に対応する前記第１のボリュームデータの領域に位置が関連する第２のボリュームデータの領域に含まれるボクセル群のヒストグラムを動的に生成する。

また、本発明の画像処理方法は、ボリュームレンダリングによる画像処理方法であって、第１のボリュームデータを用いてレンダリングによって生成した２次元画像に設定した領域に対応する前記第1のボリュームデータの領域に含まれるボクセル群のヒストグラム、及び、前記領域と位置が関連する第２のボリュームデータの領域に含まれるボクセル群のヒストグラムをそれぞれ動的に生成する。

また、本発明の画像処理方法は、前記第１もしくは前記第２のボリュームデータの前記領域は更にマスク領域に含まれる領域であるものである。また、本発明の画像処理方法は、前記ボクセル群のそれぞれのボクセルの不透明度に応じた重率を掛けたヒストグラムを生成する。また、本発明の画像処理方法は、複数のボリュームデータのヒストグラムをそれぞれ生成する。

また、本発明の画像処理方法は、前記レンダリングが、平行投影を用いるものである。また、本発明の画像処理方法は、前記レンダリングが、透視投影を用いるものである。また、本発明の画像処理方法は、前記レンダリングが、円筒投影を用いるものである。また、本発明の画像処理プログラムは、コンピュータに、本発明の画像処理方法を実行させるための画像処理プログラムである。

本発明にかかる画像処理方法および画像処理プログラムによれば、画像表示とヒストグラム表示とにより医療診断を行う場合に、ユーザの意図に応じて適切なヒストグラムを表示させることができる。

図１９は、本発明の一実施形態にかかる画像処理方法で使用されるコンピュータ断層撮影（ＣＴ）装置を概略的に示す。コンピュータ断層撮影装置は、被検体の組織等を可視化するものである。Ｘ線源１からは同図に鎖線で示す縁部ビームを有するピラミッド状のＸ線ビーム束２が放射される。Ｘ線ビーム束２は、例えば患者３である被検体を透過しＸ線検出器４に照射される。Ｘ線源１及びＸ線検出器４は、本実施形態の場合にはリング状のガントリー５に互いに対向配置されている。リング状のガントリー５は、このガントリーの中心点を通るシステム軸線６に対して、同図に示されていない保持装置に回転可能（矢印a参照）に支持されている。

患者３は、本実施形態の場合には、Ｘ線が透過するテーブル７上に寝ている。このテーブルは、図示されていない支持装置によりシステム軸線６に沿って移動可能（矢印ｂ参照）に支持されている。

従って、Ｘ線源１及びＸ線検出器４は、システム軸線６に対して回転可能でありかつシステム軸線６に沿って患者3に対して相対的に移動可能である測定システムを構成するので、患者３はシステム軸線６に関して種々の投影角及び種々の位置のもとで投射されることができる。その際に発生するＸ線検出器４の出力信号は、ボリュームデータ生成部１１１に供給され、ボリュームデータに変換される。

シーケンス走査の場合には患者３の層毎の走査が行なわれる。その際に、Ｘ線源１及びＸ線検出器４はシステム軸線６を中心に患者３の周りを回転し、Ｘ線源１及びＸ線検出器４を含む測定システムは患者３の２次元断層を走査するために多数の投影を撮影する。その際に取得された測定値から、走査された断層を表示する断層像が再構成される。相連続する断層の走査の間に、患者３はその都度システム軸線６に沿って移動される。この過程は全ての関心断層が捕捉されるまで繰り返される。

一方、スパイラル走査中は、Ｘ線源１及びＸ線検出器４を含む測定システムはシステム軸線６を中心に回転し、テーブル７は連続的に矢印ｂの方向に移動する。すなわち、Ｘ線源１及びＸ線検出器４を含む測定システムは、患者３に対して相対的に連続的にスパイラル軌道上を、患者３の関心領域が全部捕捉されるまで移動する。本実施形態の場合、同図に示されたコンピュータ断層撮影装置により、患者３の診断範囲における多数の相連続する断層信号がボリュームデータ生成部１１１に供給される。

ボリュームデータ生成部１１１で生成されたボリュームデータセットは、画像処理部１１８内のパス生成部１１２に導かれる。パス生成部１１２は、例えば観察対象の組織の中心線を表現するパスを生成し、検査のパスとして設定する。パス生成部１１２において生成されたパスは、投影画像生成部１１５に供給される。

一方、画像処理部１１８内の領域決定部１１４は、ヒストグラムを生成する対象領域を設定し、投影画像生成部１１５に供給する。なお、ヒストグラムを生成する対象領域は、後述する操作部１１３からの指示によりインタラクティブに変更可能である。

投影画像生成部１１５は、パス生成部１１２から供給されたパスに沿って仮想視点を移動させながら、領域決定部１１４から供給されたデータに従って仮想光線を放射し、観察対象組織の投影画像を生成する。投影画像生成部１１５で生成された投影画像は後処理部１１６に供給される。後処理部１１６は、投影画像とヒストグラムの合成表示、複数の投影画像の並列表示、複数の投影画像を順次表示するアニメーション表示、あるいは仮想内視鏡（ＶＥ）画像との同時表示などの処理を行う。後処理部１１６で処理された投影画像はディスプレイ１１７に供給され表示される。

また、操作部１１３は、キーボードやマウスなどからの操作信号に応じて、ヒストグラムを生成する対象領域の変更、投影画像の切り替え等の制御信号を生成し画像処理部１１８に供給する。これにより、ディスプレイ１１７に表示された投影画像を見ながら投影画像をインタラクティブに変更し、病巣を詳細に観察することができる。

図１は、本発明の実施形態にかかる画像処理方法において、動的なヒストグラム表示を行う場合の説明図である。本実施形態の画像処理方法は、ボリュームレンダリングによる画像処理方法であって、ボリュームデータ（画像データ）１１に対するレンダリングによって生成した画像内に領域を設定し、設定した領域内に投影されるボリュームデータのヒストグラムを生成するものである。例えば、図１（ａ）に示すように、設定した領域内に投影されるボリュームデータが対象ボクセル群１２の場合、対象ボクセル群１２に対するヒストグラム１３は、図１（ｃ）に示すものとなる。また、領域の設定を変更し、設定した領域内に投影されるボリュームデータを対象ボクセル群１４とした場合、対象ボクセル群１４に対するヒストグラム１５が動的に生成され、図１（ｄ）に示すものとなる。

このように本実施形態の画像処理方法では、ヒストグラム１３，１５の計算対象を、設定した領域内に投影されるボリュームデータ（対象ボクセル群１２，１４）に限定し動的に処理するため、領域の設定を変更するとヒストグラムもリアルタイムに変化する。これにより、医師等のユーザは、表示を変化させた場合にはヒストグラムが動的に生成されて表示されるため、ユーザの意図に応じた適切なヒストグラムを表示させることができ、効率的に診断を行うことができる。

次に、本実施形態の画像処理方法において、ヒストグラムの計算対象領域の指定方法を、（１）レンダリング領域、（２）マスク領域、（３）不透明領域、または上記の組み合わせとする場合について説明する。

（実施形態１）
図２は、本実施形態の画像処理方法において、レンダリング領域のヒストグラムを表示する場合を示す。図２（ａ）に示すように、レンダリング領域２２は、３次元画像データ２１のうち、レンダリングによって生成した画像内に設定した２次元領域である２次元画面２３に投影する領域である。すなわち、３次元画像データ２１を仮想光線２５で切り取った領域がレンダリング領域２２となる。この場合、レンダリング領域２２は、視点２４や２次元画面２３の大きさ（拡大率）を変える（２次元領域の設定を変える）ことにより変化する。

図２（ｂ），（ｃ），（ｄ）は、拡大率を変えた場合の画像表示と対応するヒストグラムを示す。同図に示すように、２次元画面２６の拡大率を変えると、それに対応して脂肪２８および骨２９を表わすヒストグラム２７が動的に生成されて表示される。これにより、医師等のユーザは、表示させる２次元画面２６を２次元領域として指定するだけで、２次元画面２６に対応するヒストグラム２７が動的に生成されて表示されるので、ユーザの意図に応じた適切なヒストグラムを見ることができ、円滑に診断を行うことができる。

（実施形態２）
図３は、本実施形態の画像処理方法において、マスク領域のヒストグラムを表示する場合を示す。図３（ａ）に示すように、マスク領域３２は、３次元画像データ３１のうち、３次元マスクで指定された領域である。各種演算によってマスク領域３２を編集すると、それに対応してヒストグラムが動的に生成される。

図３（ｂ），（ｃ），（ｄ）は、マスク領域３２を変えた場合の画像表示と対応するヒストグラムを示す。同図に示すように、マスク領域３２を変えると、それに対応して脂肪３６および骨３７を表わすヒストグラム３５が変化する。これにより、医師等は、マスク領域３２を指定するだけで、マスク領域３２に対応するヒストグラム３５が表示されるので、ヒストグラムの計算対象となる領域を指定する手間を省いて円滑に診断を行うことができる。

（実施形態３）
ボリュームレンダリングにおいてはボクセル値よりLUT(Look Up Table)関数を用いて、カラー値と不透明度を取得しレンダリングに用いることが行われている。特にこの時のLUT(Look Up Table)関数のうち不透明度を計算する関数をオパシティ関数と言う。
図４は、本実施形態の画像処理方法において、ボクセル値よりオパシティ関数を用いて計算した不透明度に応じたヒストグラムを表示する場合を示す。図４（ａ）に示すように、不透明領域４１は、オパシティ関数上の、不透明度がゼロでない領域（可視領域）である。不透明度は、オパシティ関数を変えると変化する。また、個々のボクセル値のヒストグラムに対する寄与度は不透明度に応じて重率をかけても良い。

図４（ｂ），（ｃ），（ｄ）は、オパシティ関数上の不透明領域４１を変えた場合の画像表示と対応するヒストグラムを示す。同図に示すように、オパシティ関数上の不透明領域４１を変えると、それに対応して脂肪４６および骨４７を表わすヒストグラム４５が変化する。これにより、医師等のユーザは、オパシティ関数上の不透明領域４１を指定するだけで、オパシティ関数上の不透明領域４１に対応するヒストグラム４５が動的に生成され表示されるので、ユーザの意図に応じた適切なヒストグラムを見て、円滑に診断を行うことができる。

図５は、本実施形態の画像処理方法において、画像内の背後に隠れている物体の情報を提示する例を示す。図５（ａ）に示すように、マスク領域５１に骨５２が表示され、ヒストグラム５３により骨５４があることを確認した場合に、図５（ｂ）に示すように、画面上のマウス操作でマスク領域５５から骨を削除したとしても、ヒストグラム５６により、画面上には表示されていない領域に骨５７がまだ残っていることがわかる。

そこで、図５（ｃ）に示すように、画面の表示位置およびマスク領域５８を変えて骨５９を探す。そして、ヒストグラム６０を見ながら残った骨６１を削除する。また、図５（ｄ）に示すヒストグラム６３により骨がすべて削除できたことを確認する。

このように、本実施形態の画像表示方法によれば、骨等を削除して組織を観察する場合に、領域を指定するだけでその領域に対応するヒストグラムが動的に生成され表示されるので、画像内の背後に隠れている物体の情報をリアルタイムで容易に取得することができる。

また、ＣＲＴ等の画面に表示されたヒストグラム上をマウスによりクリックすると、ヒストグラム上のクリック箇所の表現するボクセル値が存在するので、そのボクセル値のボリュームデータ上のボクセル値を持つボクセルを選択することが出来る。この場合、ヒストグラムの軸上のクリックすることによってもボクセル値を取得することが出来るので、そのボクセル値のボリュームデータ上のボクセル値を持つボクセルを選択することが出来る。選択されたボクセルが該当ボクセル値を持つボクセルのうちの一つをプログラムが選択しても、該当ボクセル値を持つボクセルの全てを選択しても良く、該当ボクセル値を持つボクセルを固まり毎にボクセル群として領域分けしその内の最大体積のボクセル群を選択しても良い。選択されたボクセルが直ちに画面に表示されるように画像を変更しても良いし、選択されたボクセルをボリュームデータから削除することもできるし、選択されたボクセルをマスク領域としても良い。また、ヒストグラム上で複数のボクセル値を同時に選択することも出来る。これによってユーザの関心のある箇所をボクセル値を通じて簡単に指摘することが出来る。

また、ボリュームデータは時系列情報を持っていても良く、この場合は画像がアニメーションするのにあわせて、ヒストグラムを動的に変更することが出来る。

また、互いに座標が関連づけられたボリュームデータが複数存在していても良い。複数ボリュームデータが存在する場合はボリュームレンダリングによって生成される画像は複数のボリュームの内の一つを用いて行われることもあるが、複数のボリュームデータを複合してボリュームレンダリングして画像を生成することもある。この時にヒストグラムは表示に用いられているボリュームデータに関わり合い無く、複数のボリュームデータから一つを選択してヒストグラムを表示しても良いし、複数のヒストグラムを表示しても良い。この時に表示に用いられていないボリュームデータの内のヒストグラム生成に用いるボクセルとは表示に用いられているボリュームデータのヒストグラムを計算するとしたら用いる３次元領域に対応する、ヒストグラム生成に用いるボリュームデータ内の３次元領域のボクセルである。これは、例えば、ＰＥＴ装置とＣＴ装置の双方からボリュームデータを取得して座標を関連づけて一方もしくは双方を合成して表示するような場合に効果的である。

また、本実施形態の画像処理方法は、仮想光線の投射が途中で中断した場合は中断した位置までの領域のヒストグラムを求めることが出来る。特に、レイキャスト法であっては仮想光線がボクセルを通過するのにあわせて光量が減衰する処理を行うが、光量が0になった箇所までのヒストグラムを表示することができる。このようにすれば半透明領域を表示しているときに半透明領域のうしろでかすんで見える物体をヒストグラムで表示しつつ不透明領域の背後はヒストグラムに反映させないことが出来る。

また、本実施形態の画像処理方法は、投射した仮想光線上のボクセルの最大値を取得して画像処理する方法であるＭＩＰ（Maximum Intensity Projection）に適用することもできる。ＭＩＰ法であっては光線の通過する領域に対する領域のヒストグラムを求めることも出来るし、最大値のヒストグラムを求めても良い。ＭＩＰは、ボリュームレンダリングの中では比較的簡単な計算で行うことができ、類似処理に最小値、平均値、加算値を取得する方法などがある。特に、最小値を取得するものをＭＩＮＩＰ（Minimum Intensity Projection）と言う。さらに、本実施形態の画像処理方法は、ＭＰＲのような断面に厚みをつけて切り出した上でＭＩＰ処理を行う厚み付きＭＩＰや厚み付きＭＩＮＩＰに適用することもできる。ようはボリュームデータを用いるボリュームレンダリング法であれば適用できる。

図６は、本実施形態の画像処理方法において、平行投影の場合に仮想光線の通過した点でのみヒストグラムを計算する実施例を示す。同図に示すように、２次元領域６２は、仮想光線６３がボリュームデータ６１を平行に通過する領域に対応し、その領域内のボクセルデータが投影される。このとき、本実施例では、ボリュームデータ６１のうち仮想光線６３が通過した点でのみヒストグラムを計算する。

本実施例によれば、ヒストグラムの計算が、レンダリング時の仮想光線６３の投射と同時に行えるので簡便であり、仮想光線６３の方向を指定するだけで、仮想光線６３が通過する領域（２次元領域）に対応するヒストグラムが動的に生成され表示される。このため、医師等のユーザは、ヒストグラムを計算する領域を指定する必要がなく、ユーザの意図に応じた適切なヒストグラムを見ながら、円滑に診断を行うことができる。なお、この場合、仮想光線６３が疎に投射される場合には誤差を伴う為、投射間隔を誤差が発生しない程度とするのが好ましい。

図７は、本実施形態の画像処理方法において、所定のボクセルがヒストグラムを計算する対象となるボクセルかどうか判断する方法を示す。同図に示すように、ボリュームデータ７１内の各ボクセル７３，７４が、ヒストグラムを計算する対象となるボクセルかどうかは、２次元領域７２内に投影されうるかどうかで判断する。すなわち、ボクセル７３は、２次元領域７２内に投影されるので、ヒストグラムを計算する対象となるボクセルであり、ボクセル７４は、２次元領域７２外に投影されるので、ヒストグラムを計算する対象となるボクセルではない。これにより、各ボクセルがヒストグラムの対象となるかどうかを直接的に判断することができる。

図８は、本実施形態の画像処理方法において、２次元領域を元に３次元領域を判断する関数を作成する場合を示す。この場合は、２次元領域８２を元に、ボリュームデータ８１内の３次元領域８３を判断する関数を作成する。この関数は視点などの変更により動的に変更され、これによってボクセルがヒストグラムを計算する対象となっているか否かを容易に判断できるようになる。本関数は２次元領域８２が矩形の場合は、矩形の４辺を含む直線がそれぞれが３次元空間に投影されることによって４つの平面を生成するので、ボリュームデータ８１内の領域は４つの平面によって区切られる領域内であり、容易に計算できる。本関数を用いてボリュームデータ８１の走査範囲を限定することによってヒストグラム計算ができる。

図９は、表示する画像と２次元領域は必ずしも一致しないことを示す。図９（ａ）は、表示する画像９１が、ヒストグラムを計算する２次元領域と一致する場合（表示画面全体を２次元領域として設定する場合）を示し、図９（ｂ）は、表示する画像９１が、ヒストグラムを計算する２次元領域９２と一致しない場合（表示画面の一部を２次元領域として設定する場合）を示す。同図に示すように、ボリュームデータの内、設定した２次元領域９２に投影されるボクセルのヒストグラムが計算されるので、視線方向（仮想光線の方向）に応じてヒストグラムを変更し、円滑に診断を行うことができる。

図１０は、透視投影法により画像を作成する場合に、所定のボクセルがヒストグラムを計算する対象となるボクセルかどうか判断する方法を示す。同図に示すように、透視投影法では、視点１０２から放射状に仮想光線１０６が放射される。

ボリュームデータ１０１内の各ボクセル１０３，１０４が、ヒストグラムを計算する対象となるボクセルかどうかは、各ボクセルが２次元領域１０５内に投影されうるかどうかで判断する。すなわち、ボクセル１０４は、２次元領域１０５内に投影されるので、ヒストグラムを計算する対象となるボクセルであり、ボクセル１０３は、２次元領域１０５外に投影されるので、ヒストグラムを計算する対象となるボクセルではない。これにより、各ボクセルがヒストグラムの対象となるかどうかを直接的に判断することができる。なお、円筒投影法の場合も同様に判断できる。

図１１は、仮想光線を基準にヒストグラムを求める処理の大枠のフローチャートを示す。この場合は、医師等の操作により表示画像が更新されると、メモリ上でヒストグラム領域として、ボクセル値の取りうる範囲のVmin〜Vmaxの配列である頻度freqを確保する（ステップＳ１１）。

次に、頻度freq[Vmin〜Vmax]を0で初期化し（ステップＳ１２）、描画範囲(x,y)を構成する配列のx,y各要素の要素数x_max, y_maxを取得する（ステップＳ１３）。次に、描画範囲を構成する配列のx,y要素の反復子i, j を確保する（ステップＳ１４）。画像を走査する２重ループの開始としてj = 0に設定する（ステップＳ１５）。

次に、j < y_maxかどうか判断し（ステップＳ１６）、j < y_maxの場合（Ｙｅｓ）は、i = 0に設定し（ステップＳ１７）、i < x_maxを判断する（ステップＳ１８）。そして、i < x_maxの場合（Ｙｅｓ）は、描画範囲内のP(i, j)にてボリュームデータに対して仮想光線を投射する（ステップＳ１９）。

次に、i = i+1を計算し（ステップＳ２０）、ステップＳ１８に戻る。ステップＳ１８において、i < x_maxでない場合（Ｎｏ）は、j = j+1を計算し（ステップＳ２１）、ステップＳ１６に戻る。ステップＳ１６において、j < y_maxでない場合（Ｎｏ）は終了する（ステップＳ２２）。なお、ステップ１９の詳細は図１２および図１３において説明する。

図１２は、仮想光線を基準にヒストグラムを求める処理の詳細を示す。これは、図１１のステップＳ１９の詳細なフローチャートである。この場合は、描画範囲内のP(i, j)に対応する投影開始点O（x,y,z）及びサンプリング間隔ΔS(x, y, z)を設定し（ステップＳ３１）、仮想光線現在位置X(x,y,z) = Oを設定する（ステップＳ３２）。

次に、vox = Vol(X(x,y,z))により、ボリュームデータ内の各位置でのボクセル値を取得する（ステップＳ３３）。そして、freq(vox) = freq(vox) +1により、ボクセル値に対応するヒストグラム値のカウントを１加算する（ステップＳ３４）。

次に、Xが終端位置まで来たかどうかを判断し（ステップＳ３５）、Xが終端位置でない場合（Ｎｏ）は、X(x,y,z) = X(x,y,z) + ΔS(x, y, z)により、仮想光線現在位置を前進させ（ステップＳ３６）、ステップＳ３３に戻る。一方、ステップＳ３５において、Xが終端位置まで来た場合（Ｙｅｓ）は、親ルーチン（図１１のステップＳ１９）に戻る（ステップＳ３７）。

図１３は、仮想光線を基準にヒストグラムを求める処理の詳細で、更にマスク領域とオパシティ関数の不透明度を考慮する場合を示す。これも、図１１のステップＳ１９の詳細なフローチャートである。この場合は、描画範囲内のP(i, j)に対応する投影開始点O（x,y,z）及びサンプリング間隔ΔS(x, y, z)を設定し（ステップＳ４１）、仮想光線現在位置X(x,y,z) = Oに設定する（ステップＳ４２）。

次に、vox = Vol(X(x,y,z))により、ボリュームデータ内の各位置でのボクセル値を取得し（ステップＳ４３）、op = Opacity(vox)により、voxに対応する不透明度を取得する（ステップＳ４４）。また、msk = Mask(X(x,y,z))により、位置Xに対応する不透明度を取得する（ステップＳ４５）。

次に、freq(vox) = freq(vox) + op*mskにより、ボクセル値に対応するヒストグラム値のカウントを各種不透明度に応じて加算し（ステップＳ４６）、Xが終端位置まで来たかどうかを判断する（ステップＳ４７）。そして、Xが終端位置でない場合（Ｎｏ）は、X(x,y,z) = X(x,y,z) + ΔS(x, y, z)により、仮想光線現在位置を前進させ（ステップＳ４８）、ステップＳ４３に戻る。一方、ステップＳ４７において、Xが終端位置まで来た場合（Ｙｅｓ）は、親ルーチン（図１１のステップＳ１９）に戻る（ステップＳ４９）。

図１４は、領域投射してヒストグラムを求める処理の大枠を示すフローチャートである。この場合は、医師等の操作により表示画像が更新されると、ヒストグラム領域として、ボクセル値の取りうる範囲のVmin〜Vmaxの配列により頻度freqを確保する（ステップＳ５１）。

次に、頻度freq[Vmin〜Vmax]を0で初期化し（ステップＳ５２）、ボリュームデータVol(x,y,z)を構成する配列のx,y,z各要素の要素数x_max, y_max, z_max を取得する（ステップＳ５３）。

次に、ボリュームデータ上の点から２次元面に投影する関数p(x2,y2)= Proj(x,y,z)の定義を行う（ステップＳ５４）。これは、例えば、平行投影ならp(x2, y2) = A(投影行列)(x,y,z) + B(オフセット)で関数は定義できる。

次に、２次元領域関数 flag = Area(x2,y2) {返り値flagは「領域外」か「領域内」}の定義を行う（ステップＳ５５）。これは、例えば、領域が矩形Rect(left,right,top,bottom)なら式：left<=x2 and x2<right and top<=y2 and y2<bottomが真なら「領域内」、偽なら「領域外」で関数は定義できる。次に、ヒストグラムを求めるためのループ（ステップＳ５６）となるが、これは図１５および図１６で詳細に説明する。

図１５は、領域投射してヒストグラムを求める処理において、ヒストグラムを求めるためのループを示すフローチャートである。この場合は、ボリュームデータを構成する配列のx,y,z要素の反復子i, j, k を確保し（ステップＳ６１）、ボリュームを走査する３重ループの開始としてk = 0に設定する（ステップＳ６２）。

次に、k < z_maxかどうか判断し（ステップＳ６３）、k < z_maxの場合（Ｙｅｓ）は、j = 0に設定し（ステップＳ６４）、j < y_maxかどうか判断する（ステップＳ６５）。そして、j < y_maxの場合（Ｙｅｓ）は、i = 0に設定し（ステップＳ６６）、i < x_maxかどうか判断する（ステップＳ６７）。

そして、i < x_maxの場合（Ｙｅｓ）は、関数p(x2,y2)= Proj(x,y,z)よりボリュームデータ内の位置(i,j,k)→２次元面上の位置Pos(x2, y2)を求める（ステップＳ６８）。そして、Area(Pos(x2,y2))が領域内かどうかを判断し（ステップＳ６９）、Area(Pos(x2,y2))が領域内の場合は、vox = Vol(i, j, k)により、ボリュームデータ内の各位置でのボクセル値を取得する（ステップＳ７０）。

次に、freq(vox) = freq(vox) +1により、ボクセル値に対応するヒストグラム値のカウントを１加算し（ステップＳ７１）、i = i+1を計算し（ステップＳ７２）、ステップＳ６７に戻る。一方、ステップＳ６９において、Area(Pos(x2,y2))が領域外の場合も、i = i+1を計算し（ステップＳ７２）、ステップＳ６７に戻る。

次に、ステップＳ６７において、i < x_max でない場合（Ｎｏ）は、j = j+1を計算し（ステップＳ７３）、ステップＳ６５に戻る。ステップＳ６５において、j < y_maxでない場合（Ｎｏ）は、k = k+1を計算し（ステップＳ７４）、ステップＳ６３に戻る。ステップＳ６３において、k < z_maxでない場合（Ｎｏ）は、終了する（ステップＳ７５）。

図１６は、領域投射してヒストグラムを求める処理のヒストグラムを求めるためのループにおいて、更にマスク領域とオパシティ関数の不透明度を考慮する場合を示すフローチャートである。この場合は、ボリュームデータを構成する配列のx,y,z要素の反復子i, j, k を確保し（ステップＳ８１）、ボリュームを走査する３重ループの開始としてk = 0に設定する（ステップＳ８２）。

次に、k < z_maxかどうか判断し（ステップＳ８３）、k < z_maxの場合（Ｙｅｓ）は、j = 0に設定し（ステップＳ８４）、j < y_maxかどうか判断する（ステップＳ８５）。そして、j < y_maxの場合（Ｙｅｓ）は、i = 0に設定し（ステップＳ８６）、i < x_maxかどうか判断する（ステップＳ８７）。

そして、i < x_maxの場合（Ｙｅｓ）は、関数p(x2,y2)= Proj(x,y,z)よりボリュームデータ内の位置(i,j,k)→２次元面上の位置Pos(x2, y2)を求める（ステップＳ８８）。そして、Area(Pos(x2,y2))が領域内かどうかを判断する（ステップＳ８９）。

そして、Area(Pos(x2,y2))が領域内の場合は、vox = Vol(i, j, k)により、ボリュームデータ内の各位置でのボクセル値を取得する（ステップＳ９０）。また、op = Opacity(vox)により、voxに対応する不透明度を取得する（ステップＳ９１）。また、msk = Mask(X(x,y,z))により、位置Xに対応する不透明度を取得する（ステップＳ９２）。

次に、freq(vox) = freq(vox) + op*mskにより、ボクセル値に対応するヒストグラム値のカウントを各種不透明度に応じて加算し（ステップＳ９３）、i = i+1を計算し（ステップＳ９４）、ステップＳ８７へ戻る。一方、ステップＳ８９において、Area(Pos(x2,y2))が領域外の場合も、i = i+1を計算し（ステップＳ９４）、ステップＳ８７へ戻る。

次に、ステップＳ８７において、i < x_max でない場合（Ｎｏ）は、j = j+1を計算し（ステップＳ９５）、ステップＳ８５に戻る。ステップＳ８５において、j < y_maxでない場合（Ｎｏ）は、k = k+1を計算し（ステップＳ９６）、ステップＳ８３に戻る。ステップＳ８３において、k < z_maxでない場合（Ｎｏ）は終了する（ステップＳ９７）。

以上説明したように、上記実施形態の画像処理方法によれば、画像表示とヒストグラム表示とにより医療診断を行う場合に、医師等のユーザの意図に応じて適切なヒストグラムを表示させることができる。また、対象領域の変更が対話的にヒストグラムに反映されるので、領域を指定する操作の補助として使用することができる。また、対象領域として、画像上に表示されていない領域を指定することもできるので、画像からは得られない情報をユーザに示すことができる。さらに、レンダリング領域のヒストグラムにおいては、画像内の背後に隠れている物体の情報を提示することができるようになる。

本発明の実施形態にかかる画像処理方法において、動的なヒストグラム表示を行う場合の説明図 本実施形態の画像処理方法において、レンダリング領域のヒストグラムを表示する場合を示す説明図 本実施形態の画像処理方法において、マスク領域のヒストグラムを表示する場合を示す説明図 本実施形態の画像処理方法において、不透明領域のヒストグラムを表示する場合を示す説明図 本実施形態の画像処理方法において、画像内の背後に隠れている物体の情報を提示する実施例 本実施形態の画像処理方法において、平行投影の場合に仮想光線の通過した点でのみヒストグラムを計算する実施例 本実施形態の画像処理方法において、所定のボクセルがヒストグラムを計算する対象となるボクセルかどうか判断する方法を示す説明図 本実施形態の画像処理方法において、２次元領域を元に３次元領域を判断する関数を作成する場合を示す説明図 表示する画像と２次元領域は必ずしも一致しないことを示す説明図 透視投影法により画像を作成する場合に、所定のボクセルがヒストグラムを計算する対象となるボクセルかどうか判断する方法を示す説明図 仮想光線を基準にヒストグラムを求める処理の大枠のフローチャート 仮想光線を基準にヒストグラムを求める処理の詳細を示すフローチャート 仮想光線を基準にヒストグラムを求める処理の詳細で、更にマスク領域とオパシティ関数の不透明度を考慮する場合を示すフローチャート 領域投射してヒストグラムを求める処理の大枠を示すフローチャート 領域投射してヒストグラムを求める処理において、ヒストグラムを求めるためのループを示すフローチャート 領域投射してヒストグラムを求める処理のヒストグラムを求めるためのループにおいて、更にマスク領域とオパシティ関数の不透明度を考慮する場合を示すフローチャート ３次元医療画像で用いられるヒストグラムの例 従来の画像処理方法においてヒストグラムを計算するフローチャート 本発明の一実施形態にかかる画像処理方法で使用されるコンピュータ断層撮影装置の概略ブロック図

符号の説明

１ Ｘ線源
２ Ｘ線ビーム束
３ 患者
４ Ｘ線検出器
５ ガントリー
６ システム軸線
７ テーブル
１１，２１，３１ 画像データ
１２，１４ 対象ボクセル群
１３，１５，２７，３５，４５，５３，５６，６０，６３ ヒストグラム
２２ レンダリング領域
２３，２６，３３，４２ 画面
２４，１０２ 視点
２５，６３，１０６ 仮想光線
２８，３６，４６ 脂肪
２９，３７，４７，５２，５４，５７，５９，６１ 骨
３２，５１，５５，５８，６２ マスク領域
４１ 不透明領域
６１，７１，８１，１０１ ボリュームデータ
６２，７２，８２，９２，１０５ ２次元領域
７３，７４ ボクセル
８３ ３次元領域
９１ 画像
１０３ 領域外ボクセル
１０４ 領域内ボクセル
１１１ ボリュームデータ生成部
１１２ パス生成部
１１３ 操作部
１１４ 領域決定部
１１５ 投影画像生成部
１１６ 後処理部
１１７ ディスプレイ
１１８ 画像処理部

Claims (10)

1. ボリュームレンダリングによる画像処理方法であって、
   第１のボリュームデータを用いてレンダリングによって生成され表示されている２次元画像を更新したときに、
   前記第１のボリュームデータを用いてレンダリングによって新たに生成した２次元画像に設定した領域に対応する前記第１のボリュームデータの領域に含まれるボクセル群を取得するステップと、
   取得した前記ボクセル群のヒストグラムを動的に生成するステップと、を有する画像処理方法。
2. ボリュームレンダリングによる画像処理方法であって、
   第１のボリュームデータを用いてレンダリングによって生成され表示されている２次元画像を更新したときに、
   前記第１のボリュームデータを用いてレンダリングによって新たに生成した２次元画像に設定した領域に対応する前記第１のボリュームデータの領域に位置が関連する第２のボリュームデータの領域に含まれるボクセル群を取得するステップと、
   取得した前記ボクセル群のヒストグラムを動的に生成するステップと、を有する画像処理方法。
3. ボリュームレンダリングによる画像処理方法であって、
   第１のボリュームデータを用いてレンダリングによって生成され表示されている２次元画像を更新したときに、
   前記第１のボリュームデータを用いてレンダリングによって新たに生成した２次元画像に設定した領域に対応する前記第１のボリュームデータの領域に含まれるボクセル群、及び、前記領域と位置が関連する第２のボリュームデータの領域に含まれるボクセル群を取得するステップと、
   取得した前記ボクセル群のヒストグラムをそれぞれ動的に生成するステップと、を有する画像処理方法。
4. 請求項１ないし３のいずれか一項記載の画像処理方法であって、
   前記第１もしくは前記第２のボリュームデータの前記領域は更にマスク領域に含まれる領域である画像処理方法。
5. 請求項１ないし３のいずれか一項記載の画像処理方法であって、
   前記ボクセル群のそれぞれのボクセルの不透明度に応じた重率を掛けたヒストグラムを生成する画像処理方法。
6. 請求項１ないし３のいずれか一項記載の画像処理方法であって、
   複数のボリュームデータのヒストグラムをそれぞれ生成する画像処理方法。
7. 請求項１ないし３のいずれか一項記載の画像処理方法であって、
   前記レンダリングは、平行投影を用いる画像処理方法。
8. 請求項１ないし３のいずれか一項記載の画像処理方法であって、
   前記レンダリングは、透視投影を用いる画像処理方法。
9. 請求項１ないし３のいずれか一項記載の画像処理方法であって、
   前記レンダリングは、円筒投影を用いる画像処理方法。
10. コンピュータに、請求項１ないし９のいずれか一項記載の各ステップを実行させるための画像処理プログラム。

Images