(57)【要約】      （修正有）【課題】複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置を提供する。【解決手段】ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、当該画素の画素値を算出するレンダリング手段と、を備える。【選択図】図１０
2020086803.tif 000002
【特許請求の範囲】
【請求項１】
複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
を備えるレンダリング装置。
【請求項２】
クリッピング領域が設定されている場合には、
前記有効ボクセル領域算出手段は、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記クリッピング領域に含まれるボクセルを、有効ボクセルとする、
請求項１に記載のレンダリング装置。
【請求項３】
マスク領域が設定されている場合には、
前記有効ボクセル領域算出手段は、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記マスク領域に含まれるボクセルまたは前記マスク領域に含まれないボクセルを、有効ボクセルとする、
請求項１または請求項２に記載のレンダリング装置。
【請求項４】
前記有効ボクセル領域算出手段は、
前記複数の断層画像を複数の断層画像群に分割する手段と、
分割した断層画像群ごとに、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たす有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する手段と、
算出した各有効ボクセル領域に基づいてボクセル空間全体における有効ボクセル領域を算出する手段と、を備え、
複数のＣＰＵコアを用いて、各断層画像群に対応する有効ボクセル領域を並行して算出する、
請求項１から請求項３のいずれかに記載のレンダリング装置。
【請求項５】
前記レンダリング手段は、
レンダリング画像の投影面を複数の画素領域に分割する手段と、
分割した画素領域ごとに、当該画素領域の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を求めるレンダリング処理を実行する手段と、
算出した各画素領域のレンダリング結果に基づいてレンダリング画像を生成する手段と、を備え、
複数のＣＰＵコアを用いて、各画素領域に対する各レンダリング処理を並行して実行する、
請求項１から請求項４のいずれかに記載のレンダリング装置。
【請求項６】
前記レンダリング手段は、
レンダリング画像の画素ごとに、算出された視点に近い交点を探索開始点として、視線方向に向かって最も近い有効ボクセルの座標を、仮想光線の照射を開始する起点座標として探索する起点座標探索手段と、を更に備える、
請求項１から請求項５のいずれかに記載のレンダリング装置。
【請求項７】
前記起点座標探索手段は、
算出された視点から遠い交点までの範囲を探索する、
請求項６に記載のレンダリング装置。
【請求項８】
前記起点座標探索手段は、
レンダリング画像の画素ごとに、視線方向に向かって所定のステップ幅で当該画素に対応するボクセルのボクセル値を順次取得しながら、有効ボクセルを探索する、
請求項６または請求項７に記載のレンダリング装置。
【請求項９】
前記起点座標探索手段により有効ボクセルが探索されなかった場合、予め決められた所定値を返す、
請求項８に記載のレンダリング装置。
【請求項１０】
前記起点座標探索手段は、
有効ボクセルが探索された場合、当該有効ボクセルの位置と前記ステップ幅だけ視線方向と逆の方向に戻った位置との間に存在する有効ボクセルを探索する、
請求項８に記載のレンダリング装置。
【請求項１１】
前記有効ボクセル領域算出手段は、
信号値が所定の閾値以上のボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
前記レンダリング手段は、
レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上において所定の閾値以上の信号値をもつボクセルを抽出し、抽出したボクセルの信号値に基づいて代表信号値を算出し、レンダリング画像の画素値とする、
請求項１から請求項１０のいずれかに記載のレンダリング装置。
【請求項１２】
前記レンダリング手段は、
前記代表信号値として、抽出したボクセルの信号値の最大値、最小値または平均値のいずれかを算出する、
請求項１１に記載のレンダリング装置。
【請求項１３】
前記有効ボクセル領域算出手段は、
不透明度が０でないボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
前記レンダリング手段は、
レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルの色値および不透明度に基づいて、各画素の画素値を算出する、
請求項１から請求項１０のいずれかに記載のレンダリング装置。
【請求項１４】
前記レンダリング手段は、
視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
仮想光線上のボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
ボクセル座標系における座標に基づいて近傍の複数のボクセルを特定し、
特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、前記複数のボクセルのボクセル値に基づいて、前記画素値を算出する、
請求項１から請求項１３のいずれかに記載のレンダリング装置。
【請求項１５】
前記レンダリング手段は、
特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、特定した複数のボクセルの各信号値に所定の重みを付けて加算し、前記画素値を算出する、
請求項１４に記載のレンダリング装置。
【請求項１６】
前記レンダリング手段は、
特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、特定した複数のボクセルの各色値および各不透明度にそれぞれ所定の重みを付けて加算し、前記画素値を算出する、
請求項１４に記載のレンダリング装置。
【請求項１７】
前記起点座標探索手段は、
視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
ボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
ボクセル座標系における座標が有効ボクセル領域に含まれる場合は、当該ボクセルの信号値をボクセル値として取得する、
請求項８から請求項１０のいずれかに記載のレンダリング装置。
【請求項１８】
前記起点座標探索手段は、
視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
ボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
ボクセル座標系における座標が有効ボクセル領域に含まれる場合は、当該ボクセルの不透明度をボクセル値として取得する、
請求項８から請求項１０のいずれかに記載のレンダリング装置。
【請求項１９】
コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
を含むレンダリング方法。
【請求項２０】
コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
前記コンピュータを、
ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
として機能させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本開示は、レンダリング装置、レンダリング方法、及びプログラムに関する。
【背景技術】
【０００２】
医用画像診断においては、Ｘ線ＣＴ（Computed 
Tomography）装置やＭＲＩ（Magnetic  Resonance 
Imaging）から得られる断層画像群（３次元のボリュームデータ）に対して、レイキャスティング法によるレンダリング処理を実行し、ボリュームレンダリング画像やＭＩＰ（Maximum Intensity Projection）画像などの２次元のレンダリング画像を生成・描画する技術が普及している。
【０００３】
レイキャスティング法によるレンダリング処理は、大容量の断層画像群を処理する必要があるため計算量が膨大となる。特に、レイキャスティング法において冗長な処理として、レイキャスティング計算を開始する起点（仮想光線を照射する起点）を探索する処理がある。この処理は、レンダリングの投影面から視線方向に向かってボクセルをしらみつぶしに確認しながら、視点に最も近い有効なボクセル（例えば、ボリュームレンダリングの場合は不透明なボクセル）を探索する必要がある。
【０００４】
このようにボクセルをしらみつぶしに探索する方法の最大の難点は、描画対象が存在しない背景部にあり、視線方向の全てのボクセルを確認しないと、背景部の画素であることを決定できない。すなわち、背景部の面積が広いほど多くの処理時間がかかってしまう。
【０００５】
これに対し、特許文献１では、レンダリングの投影面とボクセルデータ内の描画対象（３次元オブジェクト）との距離の分布を示すデプスマップを生成し、デプスマップを参照して、各仮想光線ごとに対応する距離以上の地点からレイキャスティング計算を行う方法が開示されている。デプスマップを参照することで、仮想光線上に描画対象が全く存在していない状態を認識することができるため、描画対象が存在しない背景部での無駄な計算が回避できる。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００６−０９９４２２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
ところで、特許文献１のようにデプスマップを生成する方法では、レンダリング結果の精度がデプスマップの生成精度に大きく依存する。このため、高精度なレンダリング結果を得るためには、デプスマップの生成精度も高める必要がある。しかしながら、デプスマップの生成精度を高めるためには、描画対象の３次元形状を詳細に反映したオクトリーデータ（ボクセルデータをブロック分割したデータ）を作成する必要があるため、演算量が膨大となり、レンダリング処理の高速化効果が十分に得られないといった問題があった。
【０００８】
本願は、上記した問題に鑑みてなされたものであり、デプスマップを生成することなく、レンダリング処理の高速化を実現する、レンダリング装置等を提供することを目的とする。
【課題を解決するための手段】
【０００９】
本開示の一実施形態によると、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、を備えるレンダリング装置が提供される。
【００１０】
また本開示の一実施形態によると、コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、を含むレンダリング方法が提供される。
【００１１】
また本開示の一実施形態によると、コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、前記コンピュータを、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、として機能させるプログラムが提供される。
【発明の効果】
【００１２】
本開示によれば、デプスマップを生成することなく、レンダリング処理の高速化を実現する。
【図面の簡単な説明】
【００１３】
【図１】ボリュームレンダリング処理の概要を示す図
【図２】レンダリング装置１のハードウェア構成を示す図
【図３】（ａ）従来のレイキャスティング法を表す概念図、（ｂ）新規のレイキャスティング法を表す概念図
【図４】ボリュームレンダリング画像を生成する場合における有効ボクセル領域Ｖｒを示す図
【図５】ＭＩＰ画像を生成する場合における有効ボクセル領域Ｖｒを示す図
【図６】ボクセル空間Ｒの分割例を示す図
【図７】有効ボクセル領域算出処理を並行して実行させる仕組みの例を示す図
【図８】レンダリングの投影面の分割例を示す図
【図９】レンダリングの投影面の分割例を示す図
【図１０】レンダリング装置１の全体の動作を示すフローチャート
【図１１】有効ボクセル領域Ｖｒを算出する処理の全体の流れを示すフローチャート
【図１２】ボリュームレンダリング画像を生成する場合において、各有効ボクセル領域Ｖｒ（ｉ）を算出する処理を示すフローチャート
【図１３】ＭＩＰ画像を生成する場合において、各有効ボクセル領域Ｖｒ（ｉ）を算出する処理を示すフローチャート
【図１４】有効ボクセル領域算出処理を並行して実行させる処理を示すフローチャート
【図１５】各ボクセル空間の各有効ボクセル領域に基づいて、ボクセル空間Ｒ全体における有効ボクセル領域Ｖｒを算出する処理を示すフローチャート
【図１６】レンダリング処理の全体の流れを示すフローチャート
【図１７】レンダリング処理を並行して実行させる処理を示す図
【図１８】レイキャスティング法に基づいてボリュームレンダリング画像を生成する処理の全体の流れを示すフローチャート
【図１９】レイキャスティング処理を示すフローチャート
【図２０】有効ボクセル領域Ｖｒと仮想光線との交点を算出する様子を示す図
【図２１】有効ボクセル領域Ｖｒと仮想光線との交点を算出する様子を示す図
【図２２】起点座標探索処理を示すフローチャート
【図２３】レイキャスティング法によるＭＩＰ画像を生成する処理の全体の流れを示すフローチャート
【図２４】レイキャスティング処理を示すフローチャート
【図２５】起点座標探索処理を示すフローチャート
【図２６】スキップ幅を変えた場合のレンダリング画像の画像精度を比較した図
【図２７】スキップ幅を変えた場合のレンダリング画像の画像精度を比較した図
【図２８】スキップ幅を変えた場合のレンダリング画像の画像精度を比較した図
【発明を実施するための形態】
【００１４】
以下図面に基づいて、本開示の実施の形態を詳細に説明する。本実施の形態では、Ｘ線ＣＴ装置により撮影されたＣＴ画像に基づいてレンダリング画像を生成する場合について説明するが、本開示はＣＴ画像以外の医用画像（ＭＲＩ画像やＰＥＴ画像等）に基づいてレンダリング画像を生成する場合にも適用可能である。
【００１５】
（１．レンダリング処理の概要）
図１は、レンダリング装置１が実行するレンダリング処理の概要を示す図である。図に示すように、レンダリング装置１は、複数の断層画像（図の例では、５１２×５１２ピクセルの３７０枚の胸部ＣＴ画像）に基づいてレイキャスティング法に基づくレンダリング処理を実行し、所定の視点から被検体を観察したレンダリング画像（ボリュームレンダリング画像やＭＩＰ（Maximum Intensity Projection）画像）を生成する。
【００１６】
（２．レンダリング装置１のハードウェア構成）
図２は、本実施の形態におけるレンダリング装置１のハードウェア構成を示すブロック図である。図２に示すように、レンダリング装置１は、制御部１１、記憶部１２、メディア入出力部１３、通信制御部１４、入力部１５、表示部１６、周辺機器Ｉ／Ｆ部１７等が、バス１８を介して接続される汎用のコンピュータで実現される。但し、これに限ることなく、用途、目的に応じて様々な構成を採ることが可能である。
【００１７】
制御部１１は、単一プロセッサ、或いはマルチコアＣＰＵで構成されるプロセッサ又はプロセッサが複数個で構成されるマルチプロセッサ（以降、いずれも単に「ＣＰＵ」と表記し、ＣＰＵが有するコアをＣＰＵコアと表記する）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フレームメモリ(Frame Memory)等によって構成される。ＣＰＵは、記憶部１２、ＲＯＭ、記録媒体等に格納されるレンダリングプログラムをＲＡＭ、フレームメモリ上のワークメモリ領域に呼び出して実行し、バス１８を介して接続された各装置を駆動制御し、レンダリング装置１が行う後述する処理を実現する。なお、本開示では、複数のＣＰＵコアを用いた処理を後述するが、複数のＣＰＵコアは物理的に複数のＣＰＵコアだけでなく、論理的に複数のＣＰＵコアも含まれる。
【００１８】
ＲＯＭは、不揮発性メモリであり、コンピュータのブートプログラムやＢＩＯＳ等のプログラム、データ等を恒久的に保持している。ＲＡＭ、フレームメモリは、揮発性メモリであり、記憶部１２、ＲＯＭ、記録媒体等からロードしたプログラム、データ等を一時的に保持するとともに、制御部１１が各種処理を行う為に使用するワークエリアを備える。
【００１９】
記憶部１２は、ＨＤＤ（Hard Disk Drive）等であり、制御部１１が実行するプログラム、プログラム実行に必要なデータ、ＯＳ（Operating System）等が格納される。プログラムに関しては、ＯＳに相当する制御プログラムや、後述する処理をコンピュータに実行させるためのアプリケーションプログラムが格納されている。これらの各プログラムコードは、制御部１１により必要に応じて読み出されてＲＡＭ、フレームメモリに移され、ＣＰＵに読み出されて各種の手段として実行される。
【００２０】
メディア入出力部１３（ドライブ装置）は、データの入出力を行い、例えば、ＣＤドライブ（−ＲＯＭ、−Ｒ、−ＲＷ等）、ＤＶＤドライブ（−ＲＯＭ、−Ｒ、−ＲＷ等）等のメディア入出力装置を有する。通信制御部１４は、通信制御装置、通信ポート等を有し、コンピュータとネットワーク間の通信を媒介する通信インタフェースであり、ネットワークを介して、他のコンピュータ間との通信制御を行う。ネットワークは、有線、無線を問わない。
【００２１】
入力部１５は、データの入力を行い、例えば、キーボード、マウス等のポインティングデバイス、テンキー等の入力装置を有する。入力部１５を介して、コンピュータに対して、操作指示、動作指示、データ入力等を行うことができる。表示部１６は、液晶パネル等のディスプレイ装置、ディスプレイ装置と連携してコンピュータのビデオ機能を実現するための論理回路等（ビデオアダプタ等）を有する。なお、入力部１５及び表示部１６は、タッチパネルディスプレイのように、一体となっていてもよい。
【００２２】
周辺機器Ｉ／Ｆ（Interface）部１７は、コンピュータに周辺機器を接続させるためのポートであり、周辺機器Ｉ／Ｆ部１７を介してコンピュータは周辺機器とのデータの送受信を行う。周辺機器Ｉ／Ｆ部１７は、ＵＳＢ（Universal Serial Bus）やＩＥＥＥ１３９４やＲＳ−２３２Ｃ等によって構成されており、通常複数の周辺機器Ｉ／Ｆを有する。周辺機器との接続形態は有線、無線を問わない。バス１８は、各装置間の制御信号、データ信号等の授受を媒介する経路である。
【００２３】
（３．データ内容）
次に、レンダリング装置１が扱うデータ内容について説明する。
【００２４】
＜断層画像群Ｄｏ＞
断層画像群Ｄｏは、被検体（人体）を頭尾軸に沿って所定のスライス間隔で連続的に撮影した複数の断層画像（ＣＴ画像）である。各断層画像はＤＩＣＯＭ形式の２次元の画像データである。ＤＩＣＯＭ形式は、１ファイルにヘッダ部と画像データ部を含む医療画像で一般的に用いられる画像フォーマットであり、画像撮影時のパラメータや診断情報を保存しておくことができる。
【００２５】
１つの断層画像は、例えば、５１２×５１２ピクセルの画像（ＣＴ画像）である。断層画像の各画素には、信号値ｖが付与されており、ＣＴ画像の場合、信号値ｖはＣＴ値である。本実施の形態では、信号値ｖ（ＣＴ値）は１６ビット（−３２７６８≦ｖ≦３２７６７）のデータとする（但し、信号値ｖのビット数は特に限定されない）。
【００２６】
ＣＴ値は、水を基準として表現した組織のＸ線減弱係数であり、ＣＴ値により組織や病変の種類等が判断できるようになっている(単位はＨＵ（Hounsfield Unit））。ＣＴ値は、水と空気のＸ線減弱係数で標準化されており、水のＣＴ値は０、空気のＣＴ値を−１０００である。また、脂肪のＣＴ値は−１２０〜−１００程度であり、通常組織のＣＴ値は０〜１２０程度であり、骨のＣＴ値は１０００前後を示す。
【００２７】
断層画像群Ｄｏは、ＸＹの２次元データである断層画像をＺ軸方向に積層したものであり、ＸＹＺの３次元データ（ボクセルデータ）として表現可能である。例えば断層画像群Ｄｏは、以下のように、ＸＹＺの３次元データ（Ｄｏ（ｘ、ｙ、ｚ））として定義される。尚、本開示において、Ｘ軸は人体の左右軸、Ｙ軸は人体の背腹軸（前後軸）、Ｚ軸は人体の頭尾軸（上下軸）に相当するものとする。ＸＹ軸は、例えば、ＣＴ装置での撮影時に設定される。また、本開示において、Ｘ軸方向を横方向、Ｙ軸方向を縦方向と呼ぶ場合がある。
【００２８】
（式１）
−３２７６８≦Ｄｏ（ｘ、ｙ、ｚ）≦３２７６７
０≦ｘ≦Ｓｘ−１、０≦ｙ≦Ｓｙ−１、０≦ｚ≦Ｓｚ−１
解像度：Ｒｘｙ、解像度Ｒｚ
【００２９】
（式１）において、Ｓｘは断層画像の横方向（Ｘ軸方向）の画素数、Ｓｙは断層画像の縦方向（Ｙ軸方向）の画素数、Ｓｚはスライス枚数（断層画像の枚数）を表す。ＲｘｙはＸ軸方向及びＹ軸方向の断層画像の解像度であり、画素の間隔の逆数、すなわち単位距離あたりの画素数を示す。Ｒｚはスライスの解像度であり、断層画像のスライス間隔（例えば、０．５ｍｍや１ｍｍ）の逆数、すなわち単位距離あたりのスライス枚数を表す。
【００３０】
後述するように（図１０のステップＳ２参照）、断層画像群Ｄｏの信号値は、必要に応じて８ビットに階調圧縮される場合がある。８ビットに階調圧縮された断層画像Ｄｏ（ｘ、ｙ、ｚ）は以下のように定義される。
【００３１】
（式２）
０≦Ｄｏ（ｘ、ｙ、ｚ）≦２５５      
０≦ｘ≦Ｓｘ−１、０≦ｙ≦Ｓｙ−１、０≦ｚ≦Ｓｚ−１
解像度：Ｒｘｙ、解像度Ｒｚ
【００３２】
＜カラーマップＣｍａｐ＞
カラーマップＣｍａｐは、信号値ｖと色値（具体的にはＲＧＢ値）及び不透明度(α値)との対応関係を定義するものであり、信号値ｖを２４ビットの色（ＲＧＢ値）及び８ビットの不透明度（α値）に変換する関数（実体的には２次元のデータテーブル）として表現可能である。例えば、１６ビットの断層画像群Ｄｏ（（式１）参照）に適用されるカラーマップＣｍａｐは次のように定義される。
【００３３】
（式３）
０≦Ｃｍａｐ（ｖ、ｎ）≦２５５
−３２７６８≦ｖ≦３２７６７
ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）、３（α）
【００３４】
また、８ビットの断層画像群Ｄｏ（（式２）参照）に適用されるカラーマップＣｍａｐ（ｖ、ｎ）は次のように定義される。
【００３５】
（式４）
０≦Ｃｍａｐ（ｖ、ｎ）≦２５５
０≦ｖ≦２５５
ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）、３（α）
【００３６】
（式３）、（式４）に示すカラーマップＣｍａｐ（ｖ、ｎ）（０≦ｎ≦３）のうち、特にカラーマップＣｍａｐ（ｖ、ｎ）（０≦ｎ≦２）は、信号値ｖを色値（ＲＧＢ値）に変換する関数に相当する。
【００３７】
また、（式３）、（式４）に示すカラーマップＣｍａｐ（ｖ、ｎ）（０≦ｎ≦３）のうち、特にカラーマップＣｍａｐ（ｖ、３）は、信号値ｖを不透明度に変換する関数に相当し、一般的に「オパシティカーブ」と呼ばれる。
【００３８】
カラーマップＣｍａｐ（ｖ、３）（オパシティカーブ）は、例えば、通常組織（１６ビットの信号値ｖ＝０〜１２０程度）の不透明度が２５５に設定（不透明度＝２５５は不透明（光がすべて反射）であることを示す）され、骨領域（１６ビットの信号値ｖ＝１０００前後）の不透明度が１〜２５４の中間値に設定（不透明度＝１〜２５４は半透明（入射光の一部が反射され、その他は透過）であることを示す）され、その他の組織等の部位を不透明度が０に設定（不透明度＝０は透明（入射光の全てが透過）であることを示す）される。
【００３９】
＜クリッピング領域ＲＯＩ＞
クリッピング領域ＲＯＩ（Region of Interest）は、断層画像群Ｄｏにおける関心領域であり、本実施の形態では、クリッピング領域ＲＯＩは直方体で設定される。クリッピング領域ＲＯＩは、以下のように、Ｘ軸方向ＲＯＩ、Ｙ軸方向ＲＯＩ、Ｚ軸方向ＲＯＩによって定義される。
【００４０】
（式５）
Ｘ軸方向ＲＯＩ：Ｘｓ−Ｘｅ
Ｙ軸方向ＲＯＩ：Ｙｓ−Ｙｅ
Ｚ軸方向ＲＯＩ：Ｚｓ−Ｚｅ
尚、クリッピング領域を設定しない場合は、Ｘｓ＝０、Ｘｅ＝Ｓｘ−１、Ｙｓ＝０、Ｙｅ＝Ｓｙ−１、Ｚｓ＝０、Ｚｅ＝Ｓｚ−１となる。
【００４１】
＜マスクデータＭａｓｋ＞
マスクデータＭａｓｋは、断層画像群Ｄｏの各画素（ｘ、ｙ、ｚ）の可視／不可視の情報を保持する３次元データであり、以下のように定義される。
【００４２】
（式６）
Ｍａｓｋ（ｘ、ｙ、ｚ）＝  ０（不可視）または１（可視）
０≦ｘ≦Ｓｘ−１、０≦ｙ≦Ｓｙ−１、０≦ｚ≦Ｓｚ−１
解像度：Ｒｘｙ、解像度Ｒｚ
【００４３】
例えば、断層画像群Ｄｏに円柱のマスク領域を設定し、マスク領域内を可視化対象とする場合、円柱内部の画素に対応するマスクデータＭａｓｋ（ｘ、ｙ、ｚ）には「１」が設定され、円柱外部の画素に対応するマスクデータＭａｓｋ（ｘ、ｙ、ｚ）には「０」が設定される。なお、マスク領域外を可視化対象としてもよい。この場合、円柱内部の画素に対応するマスクデータＭａｓｋ（ｘ、ｙ、ｚ）には「０」が設定され、円柱外部の画素に対応するマスクデータＭａｓｋ（ｘ、ｙ、ｚ）には「１」が設定される。
【００４４】
＜補正テーブルＳα＞
補正テーブルＳαは、断層画像群Ｄｏの各画素（ｘ、ｙ、ｚ）に対応する不透明度の補正倍率を格納するデータテーブルであり、以下のように定義される。補正テーブルＳαによって、オパシティカーブで規定される信号値に応じた不透明度を各画素毎に調整することができる。
【００４５】
（式７）
０≦Ｓα（ｘ、ｙ、ｚ）≦１
０≦ｘ≦Ｓｘ−１、０≦ｙ≦Ｓｙ−１、０≦ｚ≦Ｓｚ−１
【００４６】
＜ボリュームレンダリング画像Ｉｍａｇｅ＞
ボリュームレンダリング画像Ｉｍａｇｅは、レイキャスティング法に基づいてボリュームレンダリング画像を生成する処理（図１８参照）によって得られるフルカラー（２４ビット）の２次元画像（Ｓｉｚｅ×Ｓｉｚｅの画像）であり、以下のように定義される。
【００４７】
（式８）
０≦Ｉｍａｇｅ（ｘ、ｙ、ｎ）≦２５５
０≦ｘ≦Ｓｉｚｅ−１、０≦ｙ≦Ｓｉｚｅ−１
ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）
【００４８】
＜ＭＩＰ画像Ｉｍａｇｅ＞
ＭＩＰ画像Ｉｍａｇｅは、レイキャスティング法に基づいてＭＩＰ画像を生成する処理（図２３参照）によって得られるモノクロ（８ビット）の２次元画像（Ｓｉｚｅ×Ｓｉｚｅの画像）であり、以下のように定義される。ＭＩＰ画像Ｉｍａｇｅを生成する段階では、階調は元のＤＩＣＯＭ画像の信号値のままモノクロ（１６ビット）で計算する方法が一般的であり、Ｉｍａｇｅ（ｘ、ｙ）は−３２７６８≦Ｉｍａｇｅ（ｘ、ｙ）≦３２７６７の値をもつ場合があるが、表示される段階では以下のようにモノクロ（８ビット）に変換される。
【００４９】
（式９）
０≦Ｉｍａｇｅ（ｘ、ｙ）≦２５５
０≦ｘ≦Ｓｉｚｅ−１、０≦ｙ≦Ｓｉｚｅ−１
【００５０】
＜有効ボクセル領域Ｖｒ＞
有効ボクセル領域Ｖｒは、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセル（以下、「有効ボクセル」とも表記する）を全て含み、かつ有効ボクセルに外接する直方体領域であり、以下のように、Ｘ軸方向ＲＯＩ、Ｙ軸方向ＲＯＩ、Ｚ軸方向ＲＯＩによって定義される。ここで「ボクセル値」とは、各ボクセルに対応づけられる数値情報であり、モダリティで測定された信号値（例えばＣＴ撮影の場合はＣＴ値）や、信号値をカラーマップＣｍａｐに通して得られる不透明度や色値、などを含む。
【００５１】
（式１０）
Ｘ軸方向ＲＯＩ：Ｘｉｓ−Ｘｉｅ
Ｙ軸方向ＲＯＩ：Ｙｉｓ−Ｙｉｅ
Ｚ軸方向ＲＯＩ：Ｚｉｓ−Ｚｉｅ
【００５２】
（４．従来のレイキャスティング法）
図３（ａ）は、従来のレイキャスティング法を表す概念図である。レイキャスティング法では、投影面（レンダリング画像）の各画素から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出する。この際、図に示すように、レイキャスティング計算を開始する起点（視点座標系において仮想光線を照射する起点）の座標（起点座標）を探索する必要がある。
【００５３】
従来、図３（ａ）に示すように、断層画像群Ｄｏから構成されるボクセル空間Ｒにおいて、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルである有効ボクセル（例えば、不透明度が０でない値をもつボクセル）が見つかるまで、各画素から視線方向に向かってしらみつぶしにボクセルを確認していた。ボクセルを確認するにあたっては、視点座標系からボクセル座標系へ逐一座標変換を行う必要があり、また、ボリュームレンダリングではカラーマップＣｍａｐを逐一参照する必要があるため、それなりの計算負荷を伴う。
【００５４】
この方法の最大の難点は、図３（ａ）の画素１のように、有効ボクセルが存在しない背景部にあり、この場合、視線方向の全てのボクセル値を確認しないと、背景部の画素であることを決定できない。このため、背景部の面積が広いほど、多くの処理時間がかかってしまう。
【００５５】
（５．新規のレイキャスティング方法）
図３（ｂ）は、本開示における新規のレイキャスティング法を表す概念図である。本開示では、図３（ｂ）に示すように、ボクセル空間Ｒにおいて、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルである有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。
【００５６】
具体的には、ボリュームレンダリング画像を生成する場合には、図４に示すように、不透明度が０でない値（α＞０）をもつボクセル（不透明ボクセル）を有効ボクセルとし、有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出することが望ましい。
【００５７】
また、ＭＩＰ画像を生成する場合には、図５に示すように、信号値が所定の閾値以上（例えば、空気でない信号値以上）であるボクセルを有効ボクセルとし、有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出することが望ましい。
【００５８】
そして、図３（ｂ）に示すように、投影面（レンダリング画像）の各画素から伸びる視線方向のベクトル（仮想光線）と、有効ボクセル領域Ｖｒとの交点座標を幾何計算により算出し、交点座標から視線方向に向かって起点座標の探索を開始する。すなわち、従来のように投影面から起点座標の探索を開始する必要がないため、起点座標を高速に探索することが可能となる。
【００５９】
例えば、図３の画素１については、従来方法（図３（ａ）参照）では、探索範囲Ａの範囲のボクセルを探索する必要があったが、本開示の方法（図３（ｂ）参照）では、画素１については交点座標が存在しないため起点座標の探索を実行しない。すなわち、有効ボクセルが存在しない背景部において無駄な探索が一切行われない。
【００６０】
また画素２については、従来方法（図３（ａ）参照）では、探索範囲Ｂの範囲のボクセルを探索する必要があったが、本開示の方法（図３（ｂ）参照）では、交点座標から探索を開始するため、探索範囲Ｃの範囲だけの探索で済み（探索範囲Ｂ−探索範囲Ｃの範囲の探索を省略できる）、背景部以外においても起点座標の探索が高速化される。
【００６１】
なお、クリッピング領域ＲＯＩ（（式５）参照）が設定されている場合は、クリッピング領域ＲＯＩにより規定された描画範囲に含まれ、かつ、図４または図５の有効ボクセル領域に含まれるボクセル（不透明度が０でないボクセルまたは信号値が所定の閾値以上のボクセル）を有効ボクセルとする。そして、有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。
【００６２】
また、マスク領域（マスクデータＭａｓｋ、（式６）参照）が設定されている場合は、マスクにより規定された描画範囲（Ｍａｓｋ＝１のボクセル範囲）に含まれ、かつ、図４または図５の有効ボクセル領域に含まれるボクセル（不透明度が０でないボクセルまたは信号値が所定の閾値以上のボクセル）を有効ボクセルとする。そして、有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。
【００６３】
また、クリッピング領域ＲＯＩ（（式５）参照）およびマスク領域（マスクデータＭａｓｋ、（式６）参照）が設定されている場合は、クリッピング領域ＲＯＩにより規定された描画範囲およびマスクにより規定された描画範囲（Ｍａｓｋ＝１のボクセル範囲）に含まれ、かつ、図４または図５の有効ボクセル領域に含まれるボクセル（不透明度が０でないボクセルまたは信号値が所定の閾値以上のボクセル）を有効ボクセルとする。そして、有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。
【００６４】
（６．有効ボクセル領域算出処理の並行化）
また本開示では、複数のＣＰＵコアを用いて、有効ボクセル領域Ｖｒを高速に算出する。具体的には、断層画像群Ｄｏから構成されるボクセル空間ＲをＺ軸方向に等分割し、分割した各ボクセル空間に対してそれぞれ有効ボクセル領域を算出する処理（有効ボクセル領域算出処理）を実行する。この際、レンダリング装置１は、各ボクセル空間に対する各有効ボクセル領域算出処理を、レンダリング装置１が備える複数のＣＰＵコアに別々に割り当てることで、各有効ボクセル領域算出処理を並行して実行させる。そして、算出された各有効ボクセル領域を統合してボクセル空間Ｒ全体における有効ボクセル領域Ｖｒを算出する。このように、ボクセル空間Ｒを分割し、分割したボクセル空間ごとに有効ボクセル領域を並行して算出し、算出した各有効ボクセル領域を統合する方法をとることで高速に有効ボクセル領域Ｖｒを求める。
【００６５】
図６は、ボクセル空間Ｒの分割例を示す。図は、断層画像群Ｄｏから構成されるボクセル空間ＲをＺ軸方向に４分割（等分割）した例である。図に示すように、ボクセル空間Ｒが、４つのボクセル空間Ｒ（０）、Ｒ（１）、Ｒ（２）、Ｒ（３）に分割される。各ボクセル空間のＺ軸方向の上限座標、下限座標をＺ１、Ｚ２とすると、ボクセル空間Ｒ（０）はＺ１＝ＺｓからＺ２＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／４−１の範囲の断層画像群Ｄｏから構成され、ボクセル空間Ｒ（１）はＺ１＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／４からＺ２＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／２−１の範囲の断層画像群Ｄｏから構成され、ボクセル空間Ｒ（２）はＺ１＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／２からＺ２＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）・３／４−１の範囲の断層画像群Ｄｏから構成され、ボクセル空間Ｒ（３）はＺ１＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）・３／４からＺ２＝Ｚｅまでの範囲の断層画像群Ｄｏから構成される。
【００６６】
図７は、有効ボクセル領域算出処理を並行して実行させる仕組みの例を示している。図７の例では、レンダリング装置１（コンピュータ）に搭載されたマザーボード（図下）は、４つのＣＰＵ（物理コア）を備え、かつ、各ＣＰＵが２つの論理コアを有する構成（マルチコアＣＰＵ）とする。すなわち、レンダリング装置１は、実質８コア（スレッド＃１〜＃８）を備えている。はじめに、レンダリング装置１に格納されているレンダリングプログラムが起動すると、Windowsマルチタスク・オペレーティングシステムのジョブスケジューラによりＣＰＵコア＃１のスレッド＃１に割り当てられて実行される。続いて、レンダリングプログラムは、ボクセル空間Ｒ（０）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃１）、ボクセル空間Ｒ（１）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃２）、ボクセル空間Ｒ（２）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃３）、ボクセル空間Ｒ（３）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃４）の各スレッドを起動し、Windowsマルチタスク・オペレーティングシステムのジョブスケジューラに登録する。
【００６７】
Windowsマルチタスク・オペレーティングシステムでは、コアの空き状況等によって、登録された各並行処理スレッド＃１〜＃４を、８コア（スレッド＃１〜＃８）のいずれかに割り当て、各並行処理スレッド＃１〜＃４を並行して実行する。図７の例では、並行処理スレッド＃１、並行処理スレッド＃２および並行処理スレッド＃３は、ジョブスケジューラにより物理コアが空いているＣＰＵコア＃２のスレッド＃３、ＣＰＵコア＃３のスレッド＃５、ＣＰＵコア＃４のスレッド＃７に各々割り当てられる。残る並行処理スレッド＃４は、物理コアは空いていないがＣＰＵコア＃１の論理スレッド＃２に割り当てられている。ＣＰＵコア＃１の論理スレッド＃１ではレンダリングプログラム本体が使用中であるが、レンダリングプログラム本体は４つの並行処理スレッドの起動を行って、これらの終了待ち状態になっているため、ＣＰＵコア＃１の論理スレッド＃２に割り当てられている並行処理スレッド＃４が主に実行することになる。
【００６８】
なお本開示では、１台のコンピュータ（レンダリング装置１）で並行化処理を行うものとする。従って、例えばグリッドコンピューティングのように複数台のコンピュータを用いて並行化処理（分散処理）を行うことはしない。ただし、前記並行化処理を実行する１台のコンピュータをサーバーとして、コンピュータのＧＵＩ操作系のみをネットワーク接続された複数のクライアントコンピュータで実行させる、いわゆるリモートデスクトップ（マイクロソフト社の登録商標）等の構成は本開示に含む。リモートデスクトップにおけるサーバーを含め、並行化処理を１台のコンピュータで完結させることで、たとえ複数のユーザが同時に使用しても、サーバーコンピュータ間の通信が発生しないので、通信経路上でデータが喪失（データロス）したり、ネットワーク遅延等に起因する処理の遅延が生じない。また、機密性の高い断層画像群Ｄｏ（ＤＩＣＯＭ画像）が複数のクライアントコンピュータにコピーされるようなことがないため、セキュリティ対策や個人情報の取り扱いが容易となる。
【００６９】
（７．レンダリング処理の並行化）
また本開示では、複数のＣＰＵコアを用いて、レンダリング画像の投影面の画素値を計算するレンダリング処理を高速に実行する。具体的には、レンダリング画像の投影面を複数の画素領域に分割したうえで、各画素領域に対してそれぞれレンダリング処理を実行する。この際、レンダリング装置１は、各レンダリング処理を、レンダリング装置１が備える複数のＣＰＵコアに別々に割り当てることで、各レンダリング処理を並行して実行させる。レンダリング処理の並行化についても、図７と同様にWindowsマルチタスク・オペレーティングシステムの仕組みを用いて実現される。
【００７０】
図８は、レンダリングの投影面の分割例を示す。図は、Ｓｉｚｅ×Ｓｉｚｅの投影面ＰをＹ軸方向に４分割（等分割）した例である。図に示すように、投影面Ｐが、４つの画素領域Ｐ（０）、Ｐ（１）、Ｐ（２）、Ｐ（３）に分割される。各画素領域のＹ軸方向の上限座標、下限座標をＹ１、Ｙ２とすると、画素領域Ｐ（０）はＹ１＝０からＹ２＝Ｓｉｚｅ／４−１の範囲の画素から構成され、画素領域Ｐ（１）はＹ１＝Ｓｉｚｅ／４からＹ２＝Ｓｉｚｅ／２−１の範囲の画素から構成され、画素領域Ｐ（２）はＹ１＝Ｓｉｚｅ／２からＹ２＝Ｓｉｚｅ・３／４−１の範囲の画素から構成され、画素領域Ｐ（３）はＹ１＝Ｓｉｚｅ・３／４からＹ２＝Ｓｉｚｅ−１までの範囲の画素から構成される。
【００７１】
図８のように、画像を等分割することにより各ＣＰＵコアの処理時間のバランスをとり、効率のよい並行処理を行うことができる。なお、投影面の分割は図の例に限定されない。例えば、投影面をＸ軸方向に等分割してもよい（被写体のボリュームをＸ軸中心に反時計方向に９０度回転させた正面アングルでレンダリングする場合には、むしろＹ軸方向に分割する場合に比べメモリアクセスの効率が向上する）。
【００７２】
また、図９のように、投影面Ｐを２×２の画素ブロックに分割し、各ブロックの左上の画素から構成される領域を画素領域Ｐ（０）とし、各ブロックの右上の画素から構成される領域を画素領域Ｐ（１）とし、各ブロックの左下の画素から構成される領域を画素領域Ｐ（２）とし、各ブロックの右下の画素から構成される領域を処理領域Ｐ（３）としてもよい。なお、画素ブロックは２×２である必要はなく、３×３や４×４などでもよい。
【００７３】
図９の分割方法では、少なくとも１つの処理領域のレンダリング結果だけでも投影面の全体像を可視化できるメリットがある。例えば、ある画素領域のレンダリング処理が遅延した場合であっても、遅延した処理の結果を待たずに他の処理結果に基づきレンダリング結果（全体像に対して均質に間引きされた画像）を表示させ、投影面の全体像の概観を確認することができる。
【００７４】
レンダリング処理の並行化は、図７と同様の仕組みにより実現される。例えば、図７において、並列処理スレッド＃１が画素領域Ｐ（０）に係るレンダリング処理、並列処理スレッド＃２が画素領域Ｐ（１）に係るレンダリング処理、並列処理スレッド＃３が画素領域Ｐ（２）に係るレンダリング処理、並列処理スレッド＃４が画素領域Ｐ（３）に係るレンダリング処理、などに相当する。
【００７５】
（８．レンダリング装置１の動作）
次に、図１０〜図２５を参照しながら、レンダリング装置１の動作について説明する。
図１０は、レンダリング装置１の全体の動作を示すフローチャートである。
【００７６】
まず、レンダリング装置１の制御部１１は、断層画像群Ｄｏ（（式１）参照）、カラーマップＣｍａｐ（（式３）参照）、クリッピング領域ＲＯＩ（（式５）参照）、マスクデータＭａｓｋ（（式６）参照）、補正テーブルＳα（（式７）参照）を取得する（図１０のステップＳ１）。
【００７７】
続いて、制御部１１は、必要に応じて、断層画像群Ｄｏに対して階調圧縮処理を施す（図１０のステップＳ２）。具体的には、制御部１１は、以下の手順で１６ビットの断層画像群Ｄｏ（（式１）参照）を８ビットの断層画像群Ｄｏ（（式２）参照）に変換する。
【００７８】
（１）断層画像群ＤｏのＳｚ／２番目の中間スライスＤｏ（ｚ／２）における全ての画素の最小値Ｄｍｉｎ、最大値Ｄｍａｘを算出する。なおＤｍｉｎとＤｍａｘは全てのスライスより最小値と最大値を求める方法をとっても良いが、本実施の形態のように中間スライスだけで決定する方法をとることで、大容量の１６ビットの断層画像群Ｄｏをメモリ上に保持する必要がなくなる。
（２）下限値Ｌｍｉｎ＝（Ｄｍａｘ−Ｄｍｉｎ）・γ＋Ｄｍｉｎ、上限値Ｌｍａｘ＝（Ｄｍａｘ−Ｄｍｉｎ）・（１−γ）＋Ｄｍｉｎを設定する。ここで、γは階調圧縮画像のコントラスト調整幅で、０に近いほどコントラストが増大する（但し、輝度が小さくなる）。通常はγ＝０．１に設定する。
【００７９】
（３）断層画像群Ｄｏを以下の計算式で８ビットの断層画像群Ｄｏ（（式２）参照）に変換する（下限値Ｌｍｉｎ〜上限値Ｌｍａｘの範囲で２５６段階に圧縮する）。
Ｄｏ（ｘ、ｙ、ｚ）
＝（Ｄｏ（ｘ、ｙ、ｚ）−Ｌｍｉｎ）・２５５／（Ｌｍａｘ−Ｌｍｉｎ）
但し、Ｄｏ（ｘ、ｙ、ｚ）＞２５５の場合はＤｏ（ｘ、ｙ、ｚ）＝２５５、Ｄｏ（ｘ、ｙ、ｚ）＜０の場合はＤｏ（ｘ、ｙ、ｚ）＝０に飽和させる。
【００８０】
このように階調圧縮をすることで、断層画像群を保持するためのメモリ容量を半分に抑えることができる。たとえ信号値の階調が１６ビットあっても、カラーマップＣｍａｐにより、変換される色値（ＲＧＢ値）及び不透明度（α値）の階調はディスプレイの階調に合わせて８ビットに制限されるため、階調圧縮に伴う画質劣化は殆ど生じない。
【００８１】
続いて、制御部１１は、 有効ボクセル領域Ｖｒを算出する（図１０のステップＳ３）。有効ボクセル領域Ｖｒを算出する処理については後述する（「９．有効ボクセル領域算出処理」参照）。
【００８２】
続いて、制御部１１は、レンダリング処理を実行し、レンダリング画像Ｉｍａｇｅ（ボリュームレンダリング画像またはＭＩＰ画像）を生成する（図１０のステップＳ４）。この際、制御部１１は、ステップＳ３において算出した有効ボクセル領域Ｖｒを考慮してレンダリング処理を高速に実行する。レンダリング処理については後述する（「１０．レンダリング処理」「１１．ボリュームレンダリング画像の生成」「１２．ＭＩＰ画像の生成」参照）。
【００８３】
そして、制御部１１は、ステップＳ４において生成したレンダリング画像Ｉｍａｇｅ（ボリュームレンダリング画像またはＭＩＰ画像）を表示部１６に表示する（図１０のステップＳ５）。
【００８４】
（９．有効ボクセル領域算出処理）
図１１〜図１５を参照して、図１０のステップＳ３において実行される有効ボクセル領域Ｖｒを算出する処理について説明する。
図１１は、有効ボクセル領域Ｖｒを算出する処理の全体の流れを示すフローチャートである。
制御部１１は、まず、ボクセル空間ＲをＺ軸方向に等分割する（図１１のステップＳ１１）。ここでは、図６に示したように、断層画像群Ｄｏから構成されるボクセル空間ＲをＺ軸方向に４分割するものとする。すなわち、ボクセル空間Ｒ全体が、ボクセル空間Ｒ（０）、Ｒ（１）、Ｒ（２）、Ｒ（３）に分割される。
【００８５】
そして、制御部１１は、分割したボクセル空間Ｒ（ｉ）（ｉ＝０〜３）に対してそれぞれ有効ボクセル領域Ｖｒ（ｉ）（ｉ＝０〜３）を個別に算出する（図１１のステップＳ１２）。
【００８６】
図１２、図１３は、図１１のステップＳ１２において実行される各有効ボクセル領域Ｖｒ（ｉ）（ｉ＝０〜３）を算出する処理を示すフローチャートであり、図１２はボリュームレンダリング画像を生成する場合における処理、図１３はＭＩＰ画像を生成する場合における処理を示す。
【００８７】
図１２、図１３の処理は、各ボクセル空間Ｒｉ（ｉ＝０〜３）に対して並行して実行され、各ボクセル空間Ｒ（ｉ）における各有効ボクセル領域Ｖｒ（ｉ）（Ｘ軸方向ＲＯＩ、Ｙ軸方向ＲＯＩ、Ｚ軸方向ＲＯＩ）を定義するパラメータとして、各ボクセル空間Ｒ（ｉ）内の有効ボクセルのＸＹＺ軸方向の最小座標（Ｘｍｉｎ（ｉ）、Ｙｍｉｎ（ｉ）、Ｚｍｉｎ（ｉ））および最大座標（Ｘｍａｘ（ｉ）、Ｙｍａｘ（ｉ）、Ｚｍａｘ（ｉ））が算出される。
【００８８】
まず、図１２を参照して、ボリュームレンダリング画像を生成する場合における有効ボクセル領域Ｖｒ（ｉ）を算出する処理を説明する。
【００８９】
制御部１１は、まず、ボクセル（ｘ、ｙ、ｚ）の座標をｘ＝Ｘｓ、ｙ＝Ｙｓ、ｚ＝Ｚ１に初期化し、ボクセル空間Ｒ（ｉ）における有効ボクセルのＸＹＺ軸方向の最小座標（Ｘｍｉｎ（ｉ）、Ｙｍｉｎ（ｉ）、Ｚｍｉｎ（ｉ））および最大座標（Ｘｍａｘ（ｉ）、Ｙｍａｘ（ｉ）、Ｚｍａｘ（ｉ））を、Ｘｍｉｎ（ｉ）＝Ｘｅ、Ｘｍａｘ（ｉ）＝Ｘｓ、Ｙｍｉｎ（ｉ）＝Ｙｅ、Ｙｍａｘ（ｉ）＝Ｙｓ、Ｚｍｉｎ（ｉ）＝Ｚｅ、Ｚｍａｘ（ｉ）＝Ｚｓに初期化する（図１２のステップＳ３１）。
【００９０】
続いて、制御部１１は、ボクセル（ｘ、ｙ、ｚ）の不透明度Ｖαを、カラーマップＣｍａｐを参照して取得する（図１２のステップＳ３２）。
【００９１】
Ｖα＝Ｃｍａｐ（Ｄｏ（ｘ、ｙ、ｚ）、３）
【００９２】
ここでＶα＝０およびＭａｓｋ（ｘ、ｙ、ｚ）＝０のいずれも満たさない場合（有効ボクセルの場合）（図１２のステップＳ３３；Ｎｏ）、制御部１１は、次のように、ボクセル空間Ｒ（ｉ）における有効ボクセルのＸＹＺ軸方向の最小座標および最大座標を更新し（図１２のステップＳ３４）、ステップＳ３５へ移行する。
【００９３】
ｘ＜Ｘｍｉｎ（ｉ）の場合：  Ｘｍｉｎ（ｉ）＝ｘ
ｘ＞Ｘｍａｘ（ｉ）の場合：  Ｘｍａｘ（ｉ）＝ｘ
ｙ＜Ｙｍｉｎ（ｉ）の場合：  Ｙｍｉｎ（ｉ）＝ｙ
ｙ＞Ｙｍａｘ（ｉ）の場合：  Ｙｍａｘ（ｉ）＝ｙ
ｚ＜Ｚｍｉｎ（ｉ）の場合：  Ｚｍｉｎ（ｉ）＝ｚ
ｚ＞Ｚｍａｘ（ｉ）の場合：  Ｚｍａｘ（ｉ）＝ｚ
【００９４】
一方、Ｖα＝０またはＭａｓｋ（ｘ、ｙ、ｚ）＝０の場合（有効ボクセルでない場合）（図１２のステップＳ３３；Ｙｅｓ）、制御部１１は、ボクセル空間Ｒ（ｉ）における有効ボクセルのＸＹＺ軸方向の最小座標および最大座標を更新することなく、ステップＳ３５へ移行する。
【００９５】
続いて、制御部１１は、ボクセル（ｘ、ｙ、ｚ）のｘ座標をｘ←ｘ＋１に更新する（図１２のステップＳ３５）。ｘ≦Ｘｅの場合はステップＳ３２に戻り、ステップＳ３２〜Ｓ３４の処理を繰り返す。一方、ｘ＞Ｘｅの場合はボクセル（ｘ、ｙ、ｚ）のｙ座標をｙ←ｙ＋１に更新する（図１２のステップＳ３６）。ｙ≦Ｙｅの場合はステップＳ３２に戻り、ステップＳ３２〜Ｓ３４の処理を繰り返す。一方、ｙ＞Ｙｅの場合はボクセル（ｘ、ｙ、ｚ）のｚ座標を更新する（図１２のステップＳ３７）。ｚ≦Ｚ２の場合はステップＳ３２に戻り、ステップＳ３２〜Ｓ３４の処理を繰り返す。一方、ｚ＞Ｚ２の場合は処理を終了する。
すなわち、制御部１１は、ボクセル空間Ｒ（ｉ）内のボクセルをクリッピング領域ＲＯＩ内の範囲で走査し、ボクセル空間Ｒ（ｉ）内における有効ボクセル（不透明ボクセル）のＸＹＺ軸方向の最小座標および最大座標を決定する。
【００９６】
次に、図１３を参照して、ＭＩＰ画像を生成する場合における有効ボクセル領域Ｖｒ（ｉ）を算出する処理を説明する。
【００９７】
制御部１１は、まず、ボクセル（ｘ、ｙ、ｚ）の座標をｘ＝Ｘｓ、ｙ＝Ｙｓ、ｚ＝Ｚ１に初期化し、ボクセル空間Ｒ（ｉ）における有効ボクセルのＸＹＺ軸方向の最小座標（Ｘｍｉｎ（ｉ）、Ｙｍｉｎ（ｉ）、Ｚｍｉｎ（ｉ））および最大座標（Ｘｍａｘ（ｉ）、Ｙｍａｘ（ｉ）、Ｚｍａｘ（ｉ））を、Ｘｍｉｎ（ｉ）＝Ｘｅ、Ｘｍａｘ（ｉ）＝Ｘｓ、Ｙｍｉｎ（ｉ）＝Ｙｅ、Ｙｍａｘ（ｉ）＝Ｙｓ、Ｚｍｉｎ（ｉ）＝Ｚｅ、Ｚｍａｘ（ｉ）＝Ｚｓに初期化する（図１３のステップＳ４１）。
【００９８】
続いて、制御部１１は、次のように、ボクセル（ｘ、ｙ、ｚ）の８ビットの信号値Ｖ８を取得する（図１３のステップＳ４２）。
【００９９】
断層画像群Ｄｏの信号値が１６ビットの場合：
Ｖ８＝（Ｄｏ（ｘ、ｙ、ｚ）−Ｌｍｉｎ）・２５５／（Ｌｍａｘ−Ｌｍｉｎ）
断層画像群Ｄｏの信号値が８ビットの場合：
Ｖ８＝Ｄｏ（ｘ、ｙ、ｚ）
【０１００】
ここでＶ８＝０（信号値が所定の閾値未満であることを示す）およびＭａｓｋ（ｘ、ｙ、ｚ）＝０のいずれも満たさない場合（即ち、有効ボクセルの場合）（図１３のステップＳ４３；Ｎｏ）、制御部１１は、次のように、ボクセル空間Ｒ（ｉ）における有効ボクセルのＸＹＺ軸方向の最小座標および最大座標を更新し（図１３のステップＳ４４）、ステップＳ４５へ移行する。
【０１０１】
ｘ＜Ｘｍｉｎ（ｉ）の場合：  Ｘｍｉｎ（ｉ）＝ｘ
ｘ＞Ｘｍａｘ（ｉ）の場合：  Ｘｍａｘ（ｉ）＝ｘ
ｙ＜Ｙｍｉｎ（ｉ）の場合：  Ｙｍｉｎ（ｉ）＝ｙ
ｙ＞Ｙｍａｘ（ｉ）の場合：  Ｙｍａｘ（ｉ）＝ｙ
ｚ＜Ｚｍｉｎ（ｉ）の場合：  Ｚｍｉｎ（ｉ）＝ｚ
ｚ＞Ｚｍａｘ（ｉ）の場合：  Ｚｍａｘ（ｉ）＝ｚ
【０１０２】
一方、Ｖ８＝０またはＭａｓｋ（ｘ、ｙ、ｚ）＝０の場合（即ち、有効ボクセルでない場合）（図１３のステップＳ４３；Ｙｅｓ）、制御部１１は、ボクセル空間Ｒ（ｉ）における有効ボクセルのＸＹＺ軸方向の最小座標および最大座標を更新することなく、ステップＳ４５へ移行する。
【０１０３】
続いて、制御部１１は、ボクセル（ｘ、ｙ、ｚ）のｘ座標をｘ←ｘ＋１に更新する（図１３のステップＳ４５）。ｘ≦Ｘｅの場合はステップＳ４２に戻り、ステップＳ４２〜Ｓ４４の処理を繰り返す。一方、ｘ＞Ｘｅの場合はボクセル（ｘ、ｙ、ｚ）のｙ座標をｙ←ｙ＋１に更新する（図１３のステップＳ４６）。ｙ≦Ｙｅの場合はステップＳ４２に戻り、ステップＳ４２〜Ｓ３４の処理を繰り返す。一方、ｙ＞Ｙｅの場合はボクセル（ｘ、ｙ、ｚ）のｚ座標を更新する（図１３のステップＳ４７）。ｚ≦Ｚ２の場合はステップＳ４２に戻り、ステップＳ４２〜Ｓ４４の処理を繰り返す。一方、ｚ＞Ｚ２の場合は処理を終了する。
すなわち、制御部１１は、ボクセル空間Ｒ（ｉ）内のボクセルをクリッピング領域ＲＯＩ内の範囲で走査し、ボクセル空間Ｒ（ｉ）内における有効ボクセル（信号値が所定の閾値以上のボクセル）のＸＹＺ軸方向の最小座標および最大座標を決定する。
【０１０４】
図１２または図１３の処理を各ボクセル空間Ｒ（ｉ）（ｉ＝０〜３）ごとに実行することで、各ボクセル空間Ｒ（ｉ）における各有効ボクセル領域Ｖｒ（ｉ）（Ｘ軸方向ＲＯＩ、Ｙ軸方向ＲＯＩ、Ｚ軸方向ＲＯＩ）を定義するパラメータとして、各ボクセル空間Ｒ（ｉ）内の有効ボクセルのＸＹＺ軸方向の最小座標（Ｘｍｉｎ（ｉ）、Ｙｍｉｎ（ｉ）、Ｚｍｉｎ（ｉ））（ｉ＝０〜３）および最大座標（Ｘｍａｘ（ｉ）、Ｙｍａｘ（ｉ）、Ｚｍａｘ（ｉ））（ｉ＝０〜３）が算出される。
【０１０５】
ここで、前記したように、制御部１１は、各ボクセル空間Ｒ（ｉ）（ｉ＝０〜３）に対して実行される各有効ボクセル領域算出処理（図１１のステップＳ１２（図１２、図１３））を、レンダリング装置１が備える複数のＣＰＵコアに別々に割り当てることで、各有効ボクセル領域Ｖｒ（ｉ）（ｉ＝０〜３）を並行して算出する。
【０１０６】
図１４は、各ボクセル空間Ｒ（ｉ）（ｉ＝０〜３）に対する各有効ボクセル領域算出処理を並行して実行させる処理を示す。図に示すように、制御部１１（レンダリングプログラム）が、ボクセル空間Ｒ（０）（Ｚ１＝Ｚｓ〜Ｚ２＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／４−１）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃１）、ボクセル空間Ｒ（１）（Ｚ１＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／４〜Ｚ２＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／２−１）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃２）、ボクセル空間Ｒ（２）（Ｚ１＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）／２〜Ｚ２＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）・３／４−１）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃３）、ボクセル空間Ｒ（３）（Ｚ１＝Ｚｓ＋（Ｚｅ−Ｚｓ＋１）・３／４〜Ｚ２＝Ｚｅ）に係る有効ボクセル領域算出処理（図の並行処理スレッド＃４）の各並行処理スレッドを起動し（図１４のステップＳ２１）、各並行処理スレッドからの終了メッセージを待つ（図１４のステップＳ２２）。
【０１０７】
各並行処理スレッド＃１〜＃４は、図７において述べたように、Windowsマルチタスク・オペレーティングシステムのジョブスケジューラに基づいて並行して実行される。図の例では、ボクセル空間Ｒ（０）に係る並行処理スレッド＃１がＣＰＵコア＃２（スレッド＃３）で実行され、ボクセル空間Ｒ（１）に係る並行処理スレッド＃２がＣＰＵコア＃３（スレッド＃５）で実行され、ボクセル空間Ｒ（２）に係る並行処理スレッド＃３がＣＰＵコア＃４（スレッド＃７）で実行され、ボクセル空間Ｒ（３）に係る並行処理スレッド＃４がＣＰＵコア＃１（スレッド＃２）で実行される（図１４のステップＳ２３）。
【０１０８】
制御部１１は、各並行処理スレッド＃１〜＃４から処理の終了メッセージを受信し、全ての並行処理スレッドから終了メッセージを受信すると、並行処理を終了する（図１４のステップＳ２４）。
以上のように、複数のＣＰＵコアを用いることで、各ボクセル空間Ｒ（ｉ）（ｉ＝０〜３）に対する各有効ボクセル領域Ｖｒ（ｉ）（ｉ＝０〜３）が並行して算出される。
【０１０９】
図１１のフローチャートに戻る。
制御部１１は、ステップＳ１２（図１２、図１３）において算出した各ボクセル空間Ｒ（ｉ）（ｉ＝０〜３）の各有効ボクセル領域Ｖｒ（ｉ）（ｉ＝０〜３）に基づいて、ボクセル空間Ｒ全体における有効ボクセル領域Ｖｒを算出する（図１１のステップＳ１３）。
【０１１０】
図１５は、ボクセル空間Ｒ全体における有効ボクセル領域Ｖｒを算出する処理を示すフローチャートである。
制御部１１は、まず、ｉ＝０に初期化し、ボクセル空間Ｒの有効ボクセルのＸＹＺ軸方向の最小座標（Ｘｍｉｎ、Ｙｍｉｎ、Ｚｍｉｎ）および最大座標（Ｘｍａｘ、Ｙｍａｘ、Ｚｍａｘ）を、Ｘｍｉｎ＝Ｘｅ、Ｘｍａｘ＝Ｘｓ、Ｙｍｉｎ＝Ｙｅ、Ｙｍａｘ＝Ｙｓ、Ｚｍｉｎ＝Ｚｅ、Ｚｍａｘ＝Ｚｓに初期化する（図１５のステップＳ５１）。
【０１１１】
続いて、制御部１１は、次のように、ボクセル空間Ｒ内の有効ボクセルのＸＹＺ軸方向の最小座標および最大座標を更新する（図１５のステップＳ５２）。
【０１１２】
Ｘｍｉｎ（ｉ）＜Ｘｍｉｎの場合：  Ｘｍｉｎ＝Ｘｍｉｎ（ｉ）
Ｘｍａｘ（ｉ）＞Ｘｍａｘの場合：  Ｘｍａｘ＝Ｘｍａｘ（ｉ）
Ｙｍｉｎ（ｉ）＜Ｙｍｉｎの場合：  Ｙｍｉｎ＝Ｙｍｉｎ（ｉ）
Ｙｍａｘ（ｉ）＞Ｙｍａｘの場合：  Ｙｍａｘ＝Ｙｍａｘ（ｉ）
Ｚｍｉｎ（ｉ）＜Ｚｍｉｎの場合：  Ｚｍｉｎ＝Ｚｍｉｎ（ｉ）
Ｚｍａｘ（ｉ）＞Ｚｍａｘの場合：  Ｚｍａｘ＝Ｚｍａｘ（ｉ）
【０１１３】
続いて、制御部１１は、ｉ←ｉ＋１に更新し（図１５のステップＳ５３）、ｉ≧ボクセル分割数（＝４）となるまで、ステップＳ５２、Ｓ５３の処理を繰り返す。これにより、最小座標（Ｘｍｉｎ、Ｙｍｉｎ、Ｚｍｉｎ）には、各ボクセル空間Ｒ（ｉ）（ｉ＝０〜３）における有効ボクセルの最小座標（Ｘｍｉｎ（ｉ）、Ｙｍｉｎ（ｉ）、Ｚｍｉｎ（ｉ））（ｉ＝０〜３）のうち、最も小さい座標値が割り当てられ、最大座標（Ｘｍａｘ、Ｙｍａｘ、Ｚｍａｘ）には、各ボクセル空間Ｒ（ｉ）（ｉ＝０〜３）における最大座標（Ｘｍａｘ（ｉ）、Ｙｍａｘ（ｉ）、Ｚｍａｘ（ｉ））（ｉ＝０〜３）のうち、最も大きい座標値が割り当てられる。
【０１１４】
続いて、制御部１１は、モアレ対策として、以下のように、α値が０の透明ボクセルで構成される境界面（ｘ＝Ｘｉｓ、ｘ＝Ｘｉｅ、ｙ＝Ｙｉｓ、ｙ＝Ｙｉｅ、ｚ＝Ｚｉｓ、またはｚ＝Ｚｉｅのいずれかを満たす）を有効ボクセル領域Ｖｒ（外接直方体）のエッジに１ボクセル幅だけ付加する補正を行う（図１５のステップＳ５４）。
【０１１５】
Ｘｉｓ＝Ｘｍｉｎ−１とし、Ｘｉｓ＜Ｘｓの場合、Ｘｉｓ＝Ｘｓ
Ｘｉｅ＝Ｘｍａｘ＋１とし、Ｘｉｅ＞Ｘｅの場合、Ｘｉｅ＝Ｘｅ
Ｙｉｓ＝Ｙｍｉｎ−１とし、Ｙｉｓ＜Ｙｓの場合、Ｙｉｓ＝Ｙｓ
Ｙｉｅ＝Ｙｍａｘ＋１とし、Ｙｉｅ＞Ｙｅの場合、Ｙｉｅ＝Ｙｅ
Ｚｉｓ＝Ｚｍｉｎ−１とし、Ｚｉｓ＜Ｚｓの場合、Ｚｉｓ＝Ｚｓ
Ｚｉｅ＝Ｚｍａｘ＋１とし、Ｚｉｅ＞Ｚｅの場合、Ｚｉｅ＝Ｚｅ
【０１１６】
なお本補正は必須ではなく、補正をしない場合には、制御部１１は、Ｘｉｓ＝Ｘｍｉｎ、Ｘｉｅ＝Ｘｍａｘ、・・・のように、ステップＳ５２で求めた座標をそのまま設定する。また、ＭＩＰ画像を生成する場合は、モアレ対策が不要なため、本補正は行わない。
【０１１７】
以上の処理により、ボクセル空間Ｒ全体における有効ボクセル領域Ｖｒ（Ｘ軸方向ＲＯＩ、Ｙ軸方向ＲＯＩ、Ｚ軸方向ＲＯＩ）（（式１０）参照）を定義するパラメータとして、ボクセル空間Ｒ内の有効ボクセルのＸＹＺ軸方向の最小座標（Ｘｉｓ、Ｙｉｓ、Ｚｉｓ）および最大座標（Ｘｉｅ、Ｙｉｅ、Ｚｉｅ）が算出される。
【０１１８】
（１０．レンダリング処理）
図１６は、図１０のステップＳ４において実行されるレンダリング処理の全体の流れを示すフローチャートである。
制御部１１は、まず、レンダリングの投影面を複数の画素領域に分割する（図１６のステップＳ６１）。ここでは、図３に示したように、レンダリングの投影面をＹ軸方向に４分割（等分割）するものとする。すなわち、レンダリングの投影面が、４つの画素領域Ｐ（０）、Ｐ（１）、Ｐ（２）、Ｐ（３）に分割される。
【０１１９】
そして、制御部１１は、分割した各画素領域Ｐ（０）〜Ｐ（３）に対してそれぞれレンダリング処理を実行する（図１６のステップＳ６２）。ここで、前記したように、各画素領域Ｐ（０）〜Ｐ（３）に対して実行されるレンダリング処理を、レンダリング装置１が備える複数のＣＰＵコアに別々に割り当てることで、各レンダリング処理を並行して実行する。
【０１２０】
図１７は、レンダリング処理を並行して実行させる処理を示す図である。図に示すように、制御部１１（レンダリングプログラム）が、画素領域Ｐ（０）に係るレンダリング処理（図の並行処理スレッド＃１）、画素領域Ｐ（１）に係るレンダリング処理（図の並行処理スレッド＃２）、画素領域Ｐ（２）に係るレンダリング処理（図の並行処理スレッド＃３）、画素領域Ｐ（３）に係るレンダリング処理（図の並行処理スレッド＃４）の各並行処理スレッドを起動し（図１７のステップＳ７１）、各並行処理スレッドからの終了メッセージを待つ（図１７のステップＳ７２）。
【０１２１】
各並行処理スレッド＃１〜＃４は、Windowsマルチタスク・オペレーティングシステムのジョブスケジューラに基づいて並行して実行される。図の例では、画素領域Ｐ（０）に係る並行処理スレッド＃１がＣＰＵコア＃２（スレッド＃３）で実行され、画素領域Ｐ（１）に係る並行処理スレッド＃２がＣＰＵコア＃３（スレッド＃５）で実行され、画素領域Ｐ（２）に係る並行処理スレッド＃３がＣＰＵコア＃４（スレッド＃７）で実行され、画素領域Ｐ（３）に係る並行処理スレッド＃４がＣＰＵコア＃１（スレッド＃２）で実行される（図１７のステップＳ７３）。
【０１２２】
制御部１１は、各並行処理スレッド＃１〜＃４から処理の終了メッセージを受信し、全ての並行処理スレッドから終了メッセージを受信すると、並行処理を終了する（図１７のステップＳ７４）。なお、前述したように、図９のような分割方法を採用した場合には、図１７のステップＳ７４において、制御部１１は、全ての並行処理スレッドの終了を待たなくてもよい。例えば、個々の並行処理スレッドの処理の終了メッセージを受信する度に、後述の図１６のステップＳ６３において、制御部１１は、終了メッセージを受信した並行処理スレッドの結果（一部の画素領域のレンダリング結果）に基づいてレンダリング画像を部分的に生成してもよい。そうすると、一部の並行処理スレッドを実行させるためのＣＰＵコアの割り当てが遅延し、当該並行処理スレッドが終了していなくても、少なくとも１つの並行処理スレッドが終了していれば、間引きされた状態でレンダリング画像の全体像を確認することができる。
【０１２３】
図１６のフローチャートに戻る。
制御部１１は、ステップＳ６２において得られた各画素領域Ｐ（０）〜Ｐ（３）のレンダリング画像を結合して、全体のレンダリング画像Ｉｍａｇｅ（ボリュームレンダリング画像またはＭＩＰ画像）を生成する（図１６のステップＳ６３）。
【０１２４】
以降、図１６のステップＳ６２において実行されるレンダリング処理として、レイキャスティング法に基づいてボリュームレダンリグ画像を生成する処理例と、レイキャスティング法に基づいてＭＩＰ画像を生成する処理例を順に説明する。
【０１２５】
（１１．ボリュームレンダリング画像の生成）
図１８は、レイキャスティング法に基づいてボリュームレンダリング画像を生成する処理の全体の流れを示すフローチャートである。この処理は、前記したように、投影面を分割した各画素領域Ｐ（０）〜Ｐ（３）に対して並行して実行される。
【０１２６】
まず、制御部１１は、座標変換パラメータを設定する（図１８のステップＳ８１）。本開示では、後述するレイキャスティング処理（図１８のステップＳ８３、図１９）で、視点座標系の各３次元座標毎にボクセルを参照するため、逐次座標変換を行う方法をとる。そのため、制御部１１は、以下のように、各座標変換処理に共通する座標変換パラメータを事前に設定・算出しておく。
【０１２７】
＜座標変換パラメータ＞
回転パラメータ行列Ｒ：
Ｒ＝［Ｒ１１  Ｒ１２  Ｒ１３；
Ｒ２１  Ｒ２２  Ｒ２３；
Ｒ３１  Ｒ３２  Ｒ３３］
（ボクセル座標系から視点座標系への座標変換を行うための３×３の行列の逆行列、ＧＵＩ側はボクセル座標系から視点座標系への座標変換を指示するが、レンダリング側は視点座標系からボクセル座標系に座標変換を行う。）
ＸＹＺ軸方向のオフセット：
Ｘｏｆｆ、Ｙｏｆｆ、Ｚｏｆｆ（２次元画面上で指定するため、通常Ｚｏｆｆ＝０）
クリッピング領域：
Ｘ軸方向ＲＯＩ：Ｘｓ−Ｘｅ
Ｙ軸方向ＲＯＩ：Ｙｓ−Ｙｅ
Ｚ軸方向ＲＯＩ：Ｚｓ−Ｚｅ
有効ボクセル領域（外接直方体領域）：
Ｘ軸方向ＲＯＩ：Ｘｉｓ−Ｘｉｅ
Ｙ軸方向ＲＯＩ：Ｙｉｓ−Ｙｉｅ
Ｚ軸方向ＲＯＩ：Ｚｉｓ―Ｚｉｅ
拡大縮小倍率Ｓｃａｌｅ（ＸＹＺ軸方向で同一）
Ｚ軸方向変倍率Ｓｃｚ＝Ｒｘｙ／Ｒｚ
座標変換サブサンプル・オフセット：Ｘ軸方向ｄｘ、Ｙ軸方向ｄｙ、Ｚ軸方向ｄｚ
生成するボリュームレンダリング画像のサイズ：Ｓｉｚｅ（ＸＹ軸方向で同一）
仮想光線のサブサンプリング倍率：Ｓｒａｙ（通常は１で、値が１より大きいと粗くなり、１未満だと高精細になる。ボリュームレンダリング画像を生成する場合はＳｒａｙ＝１、ＭＩＰ画像を生成する場合はＳｒａｙ＝２などに設定する。）
【０１２８】
制御部１１は、上記した回転パラメータ行列Ｒに対して、初期値として単位行列を設定する。すなわち、Ｒ１１＝１、Ｒ１２＝０、Ｒ１３＝０、Ｒ２１＝０、Ｒ２２＝１、Ｒ２３＝０、Ｒ３１＝０、Ｒ３２＝０、Ｒ３３＝１と設定する。
そして、ＧＵＩの指示に従い、Ｘ軸中心回転Ｒｘ、Ｙ軸中心回転Ｒｙ、Ｚ軸中心回転Ｒｚ（角度単位：ラジアン）のいずれかを逐次指定し、以下のように、各々回転行列Ａを生成して回転パラメータ行列Ｒに右から乗算して、回転パラメータ行列Ｒを更新する。これにより、ＧＵＩの指示により生成されるボクセル座標系から視点座標系への回転行列の逆行列が算出される。
【０１２９】
回転行列Ａを
Ａ＝［Ａ１１  Ａ１２  Ａ１３；
Ａ２１  Ａ２２  Ａ２３；
Ａ３１  Ａ３２  Ａ３３］
とすると、
Ｘ軸中心回転Ｒｘの場合の回転行列Ａの各要素は、
Ａ１１＝１、Ａ１２＝０、Ａ１３＝０
Ａ２１＝０、Ａ２２＝ｃｏｓＲｘ、Ａ２３＝ｓｉｎＲｘ
Ａ３１＝０、Ａ３２＝ｓｉｎＲｘ、Ａ３３＝ｃｏｓＲｘ
Ｙ軸中心回転Ｒｙの場合の回転行列Ａの各要素は、
Ａ１１＝ｃｏｓＲｙ、Ａ１２＝０、Ａ１３＝ｓｉｎＲｙ
Ａ２１＝０、Ａ２２＝１、Ａ２３＝０
Ａ３１＝−ｓｉｎＲｙ、Ａ３２＝０、Ａ３３＝ｃｏｓＲｙ
Ｚ軸中心回転Ｒｚの場合の回転行列Ａの各要素は、
Ａ１１＝ｃｏｓＲｚ、Ａ１２＝ｓｉｎＲｚ、Ａ１３＝０
Ａ２１＝−ｓｉｎＲｚ、Ａ２２＝ｃｏｓＲｚ、Ａ２３＝０
Ａ３１＝０、Ａ３２＝０、Ａ３３＝１
となる。
回転パラメータ行列Ｒは、Ｒ←Ｒ×Ａと更新される。
【０１３０】
以上で定義された座標変換パラメータに基づく座標変換処理は、レイキャスティング処理（図１８のステップＳ８３、図１９）の各処理の中で逐次実行される。
【０１３１】
続いて、制御部１１は、処理回数ｌ＝０、サブサンプル・オフセット値をｄｘ＝ｄｙ＝ｄｚ＝０に初期化する（図１８のステップＳ８２）。
そして、制御部１１は、Ｓｉｚｅ×（Ｙ２−Ｙ１＋１）サイズ（投影面をＹ軸方向に４分割（等分割）したサイズ）の画素領域に対して、レイキャスティング処理（画素領域の各座標（ｘ、ｙ）毎に色値を算出する処理）を実行する（図１８のステップＳ８３）。
【０１３２】
図１９のフローチャートを参照して、このレイキャスティング処理について説明する。
制御部１１は、まず、生成する２４ビット（ＲＧＢ）のレンダリング画像Ｉｍａｇｅ（ｘ、ｙ、ｎ）の初期値を全て０に設定する（Ｉｍａｇｅ（ｘ、ｙ、ｎ）＝０、ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ））。そして、サブサンプル回数Ｌとして、各２次元座標（ｘ、ｙ）（０≦ｘ≦Ｓｉｚｅ−１、０≦ｙ≦Ｓｉｚｅ−１）に対して、以下の処理を実行する。
【０１３３】
制御部１１は、背景色ｂｇ（ｎ）（ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）、０≦ｂｇ（ｎ）≦２５５）を設定し、ｘ＝０、ｙ＝Ｙ１とする（図１９のステップＳ９１）。
【０１３４】
続いて、制御部１１は、仮想光強度Ｔｒａｎｓ＝１．０、累積輝度値Ｅｎｅｒｇｙ（ｎ）＝０．０（０≦ｎ≦２）に初期化する（図１９のステップＳ９２）。
【０１３５】
続いて、制御部１１は、有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する（図１９のステップＳ９３）。有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する処理については後述する（「１１−１．交点算出処理」参照）
【０１３６】
続いて、制御部１１は、ＺｓｔをＺｓｔ＝Ｚｃに設定し（図１９のステップＳ９４）、Ｚｓｔから先頭の有効ボクセル（不透明ボクセル）のＺ座標ｚ（仮想光線の照射を開始する起点座標）を探索する（図１９のステップＳ９５）。すなわち、ステップＳ９３において算出された交点から視線方向に向かってレイキャスティング計算を開始する起点座標を探索する。起点座標ｚを探索する処理については後述する（「１１−２．起点座標探索処理」、図２２参照）。
【０１３７】
ｚ＜０の場合、ステップＳ１０１へ移行する。
【０１３８】
一方、ｚ≧０の場合、制御部１１は、３次元座標（ｘ、ｙ、ｚ）を座標変換してボクセルα値（Ｖα）を取得し、Ａｌｐｈａ＝Ｖα／２５５を求める（図１９のステップＳ９６）。座標変換を行い、ボクセルα値（Ｖα）を取得する処理については後述する（「１１−４．ボクセルα値、ＲＧＢ値の取得（補間あり）」参照）。
【０１３９】
ステップＳ９６において算出したＡｌｐｈａがＡｌｐｈａ＝０かつｚ＝０の場合、Ｚｓｔ＝ｚ−１としたうえで（ステップＳ９６）、ステップＳ９５に戻り、不透明ボクセルのＺ座標ｚを再度探索する。
【０１４０】
一方、ステップＳ９６において算出したＡｌｐｈａがＡｌｐｈａ＜０（有効ボクセル領域Ｖｒを通過）、またはＡｌｐｈａ＝０かつｚ＝０の場合、ステップＳ１０１へ移行し、画素（ｘ、ｙ）におけるＲＧＢ画素値を決定する。これにより、有効ボクセル領域Ｖｒを逸脱（通過）している場合（Ａｌｐｈａ＜０）は、その先に描画対象のボクセルは存在しないため、レイキャスティング処理を早期に打ち切り、冗長な処理を省略する。
【０１４１】
一方、ステップＳ９６において算出したＡｌｐｈａがＡｌｐｈａ＞０の場合、制御部１１は、３次元座標（ｘ、ｙ、ｚ）を座標変換してボクセルＲＧＢ値Ｖｃ（ｎ）を取得する（図１９のステップＳ９８）。座標変換を行い、ボクセルＲＧＢ値（Ｖｃ（ｎ））を取得する処理については後述する（「１１−４．ボクセルα値、ＲＧＢ値の取得（補間あり）」参照）。
【０１４２】
続いて、制御部１１は、仮想光線サブサンプルに基づくα値を、
Ａｌｐｈａ’＝１−（１−Ａｌｐｈａ）１／Ｓｒａｙ
と補正し、
累積輝度を、
Ｅｎｅｒｇｙ（ｎ）＝Ｅｎｅｒｇｙ（ｎ）＋Ｔｒａｎｓ／Ａｌｐｈａ・Ｖｃ（ｎ）／２５５
透過光強度を、
Ｔｒａｎｓ＝Ｔｒａｎｓ・（１．０−Ａｌｐｈａ）
と更新する（図１９のステップＳ９９）。
【０１４３】
Ｔｒａｎｓ＜０．００１の場合、ステップＳ１０１へ移行し、画素（ｘ、ｙ）におけるＲＧＢ画素値を決定する。Ｔｒａｎｓ≧０．００１の場合、制御部１１は、ｚ＝ｚ−１に更新し（図１９のステップＳ１００）、ｚ≧０の場合、ステップＳ９６に戻り、ｚ＜０の場合、ステップＳ１０１へ移行し、画素（ｘ、ｙ）におけるＲＧＢ画素値を決定する。
【０１４４】
ステップＳ１０１において、制御部１１は、画素（ｘ、ｙ）におけるＲＧＢ画素値を次のように決定する。
【０１４５】
Ｉｍａｇｅ（ｘ、ｙ、ｎ）＝Ｉｍａｇｅ（ｘ、ｙ、ｎ）＋ｋ・Ｅｎｅｒｇｙ（ｎ）・Ｌｉｇｈｔ（ｎ）／Ｌ
ここで、ｋは強度倍率であり、初期値はｋ＝１．０に設定されている。
【０１４６】
続いて、制御部１１は、ｘ←ｘ＋１に更新し（図１９のステップＳ１０２）、ｘ＜Ｓｉｚｅの場合、ステップＳ９２に戻り、次の画素ｘのＲＧＢ値を決定する処理（ステップＳ９２〜Ｓ１０１）を実行する。ｘ≧Ｓｉｚｅの場合、ｙ←ｙ＋１に更新し（図１９のステップＳ１０３）、ｙ≦Ｙ２の場合、ｘ＝０としたうえで、ステップＳ９２に戻り、ｙ行目の画素（ｘ、ｙ）のＲＧＢ値を決定する処理（ステップＳ９２〜Ｓ１０１）を繰り返す。ｙ＞Ｙ２の場合、制御部１１は、処理を終了する。
すなわち、画素領域の全画素のＲＧＢ値が算出されるまで（ステップＳ１０２；ｘ≧Ｓｉｚｅ、かつ、ステップＳ１０３；ｙ＞Ｙ２）、ステップＳ９２〜Ｓ１０１の処理を繰り返す。
【０１４７】
図１８のフローチャートに戻る。
制御部１１は、処理回数をｌ←ｌ＋１に更新し、サブサンプル・オフセット値をｄｘ←ｄｘ＋１／Ｌ、ｄｙ←ｄｙ＋１／Ｌ、ｄｚ←ｄｚ＋１／Ｌに更新する（図ＬのステップＳ８４）。制御部１１は、処理回数ｌがｌ＞Ｌ−１を満たすまで（ステップＳ８４；ｌ＞Ｌ−１）、ステップＳ８３のレイキャスティング処理を繰り返す。レイキャスティング処理が終了すると、制御部１１は、画素領域に対応するレンダリング画像Ｉｍａｇｅ（ｘ、ｙ、ｎ）（ｎ≦０≦２）を出力する（図１８のステップＳ８５）。
【０１４８】
（１１−１．交点算出処理）
図１９のステップＳ９３において実行される有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する処理について説明する。
【０１４９】
制御部１１は、以下の処理手順で、レイキャスティング処理の算出対象画素（ｘ、ｙ）に対して、視点座標系において、ｚ＝Ｚｏ（＝Ｓｉｚｅ／Ｓｒａｙ−１）からｚ＝０の範囲で、視点ｚ＝Ｚｏに最も近い有効ボクセル領域Ｖｒ（［Ｘｉｓ、Ｘｉｅ］、［Ｙｉｓ、Ｙｉｅ］、［Ｚｉｓ、Ｚｉｅ］）との交点Ｚｃを算出する。
【０１５０】
（１）視線ベクトル（仮想光線ベクトル）の算出
制御部１１は、図２０に示すように、視点座標系における視点座標（ｘ、ｙ、Ｚｏ）および下限座標（ｘ、ｙ、０）に対して各々座標変換を行い、ボクセル座標系における視点座標（ｘ１、ｙ１、ｚ１）および下限座標（ｘ２、ｙ２、ｚ２）を求め、視線ベクトル（ｖｘ、ｖｙ、ｖｚ）の各要素を次のように算出する。
【０１５１】
ｖｘ＝ｘ２−ｘ１
ｖｙ＝ｙ２−ｙ１
ｖｚ＝ｚ２−ｚ１
【０１５２】
（２）視線ベクトルと有効ボクセル領域Ｖｒとの交点の算出
制御部１１は、ｔｘ＝ｔｙ＝ｔｚ＝１０と初期化し、
１）｜ｖｘ｜≧１の場合、ｔｘ１＝（Ｘｉｓ−ｘ１）／ｖｘ、ｔｘ２＝（Ｘｉｅ−ｘ１）／ｖｘを算出し、いずれか小さい方をｔｘに設定する。
２）｜ｖｙ｜≧１の場合、ｔｙ１＝（Ｙｉｓ−ｙ１）／ｖｙ、ｔｙ２＝（Ｙｉｅ−ｙ１）／ｖｙを算出し、いずれか小さい方をｔｙに設定する。
３）｜ｖｚ｜≧１の場合、ｔｚ１＝（Ｚｉｓ−ｚ１）／ｖｚ、ｔｚ２＝（Ｚｉｅ−ｚ１）／ｖｚを算出し、いずれか小さい方をｔｚに設定する。
【０１５３】
これにより、図２０（右図）に示すように、ボクセル座標系における交点座標が求まる。具体的には、交点座標は（ｖｘ・ｔ＋ｘ１，ｖｙ・ｔ＋ｙ１，ｖｚ・ｔ＋ｚ１）で与えられ（ｔは（２）で算出されたｔｘ、ｔｙ、ｔｚのいずれかの値）、図２０の例ではｚ＝Ｚｉｓおよびｚ＝Ｚｉｅの２面に交点が存在し、ｚ＝Ｚｉｅにおける交点座標は（ｖｘ・ｔｚ＋ｘ１，ｖｙ・ｔｚ＋ｙ１，ｖｚ・ｔｚ＋ｚ１）と算出できる。交点座標は、視線ベクトルと有効ボクセル領域Ｖｒ（直方体）を構成する６面との交点（通常、交点は２点以上存在する）のうち、視点に最も近い交点である。
【０１５４】
（３）視点座標系における交点Ｚｃの決定
制御部１１は、視点座標系における交点座標ＺｃをＺｃ＝−１と初期化し、
１）ｔｘ≦１の場合
Ｙ＝ｖｙ・ｔｘ＋ｙ１、Ｚ＝ｖｚ・ｔｘ＋ｚ１を算出する。
Ｙｉｓ≦Ｙ≦ＹｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅの場合（ボクセル座標系における交点が有効ボクセル領域Ｖｒ内に存在する場合）、対応する視点座標系における交点Ｚ２＝Ｚｏ・（１−ｔｘ）を算出し、Ｚ２＞ＺｃならばＺｃ＝Ｚ２に置換する。
２）ｔｙ≦１の場合
Ｘ＝ｖｘ・ｔｙ＋ｘ１、Ｚ＝ｖｚ・ｔｙ＋ｚ１を算出する。
Ｘｉｓ≦Ｘ≦ＸｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅの場合（ボクセル座標系における交点が有効ボクセル領域Ｖｒ内に存在する場合）、対応する視点座標系における交点Ｚ２＝Ｚｏ・（１−ｔｙ）を算出し、Ｚ２＞ＺｃならばＺｃ＝Ｚ２に置換する。
３）ｔｚ≦１の場合
Ｘ＝ｖｘ・ｔｚ＋ｘ１、Ｙ＝ｖｙ・ｔｚ＋ｙ１を算出する。
Ｘｉｓ≦Ｘ≦ＸｉｅかつＹｉｓ≦Ｙ≦Ｙｉｅの場合（ボクセル座標系における交点が有効ボクセル領域Ｖｒ内に存在する場合）、対応する視点座標系における交点Ｚ２＝Ｚｏ・（１−ｔｚ）を算出し、Ｚ２＞ＺｃならばＺｃ＝Ｚ２に置換する。
【０１５５】
これにより、図２０（左図）に示すように、視点座標系における視点に最も近い交点座標（ｘ、ｙ、Ｚｃ）が求まる。
【０１５６】
（４）視点より手前に算出される交点の補正
ｔｘ、ｔｙ、ｔｚは負値をとる場合があり、その場合は交点が視点より手前になる。制御部１１は、Ｚｃ＞Ｚｏの場合、Ｚｃ＝Ｚｏに補正する。すなわち、図２１のように、視点座標系における交点座標が視点より手前（不可視領域）に算出された場合、当該交点座標を視点位置に補正する。
【０１５７】
（１１−２．起点座標探索処理）
図２２のフローチャートを参照して、図１９のステップＳ９５において実行される、起点座標ｚを探索する処理について説明する。
【０１５８】
まず、制御部１１は、ｚをｚ＝Ｚｓｔ（＝交点座標Ｚｃ）に初期化し、探索対象画素（ｘ、ｙ）を入力する（図２２のステップＳ１１１）。続いて、制御部１１は、３次元座標（ｘ、ｙ、ｚ）に対して座標変換を行い、ボクセルのα値を取得する（図２２のステップＳ１１２）。座標変換を行い、ボクセルのα値を取得する処理は後述する（「１１−３．ボクセルα値を取得（補間なし）」参照）。
【０１５９】
ステップＳ１１２において取得したボクセルのα値がα＝０の場合（ボクセルが透明の場合）、制御部１１は、ｚ←ｚ−ｍ（例えばｍ＝４）に更新する（図２２のステップＳ１１３）。ｍは探索時のＺ軸方向のスキップ幅であり、ｍが大きいほど高速に探索できるが、スキップ幅が大きすぎると起点座標の探索精度が悪化する場合がある。
【０１６０】
ステップＳ１１３において更新したｚがｚ＜０の場合（ボクセル空間外の場合）、制御部１１は、ｚ＝−１（起点座標が存在しないことを示す値）を返し（図２２のステップＳ１１４）、処理を終了する。ｚ≧０の場合、制御部１１は、ステップＳ１１２に戻る。
【０１６１】
ステップＳ１１２において取得したボクセルのα値がα＜０の場合（有効ボクセル領域Ｖｒを通過した場合）は、制御部１１は、ｚ＝−１（起点座標が存在しないことを示す値）を返し（図２２のステップＳ１１４）、処理を終了する。
【０１６２】
ステップＳ１１２において取得したボクセルのα値がα＞０の場合（ボクセルが不透明の場合）、以下のように、Ｚ軸の＋方向（視線逆方向）に戻って正確な起点座標を探索する。
【０１６３】
制御部１１は、まず、ｉをｉ＝０に初期化し（図２２のステップＳ１１５）、ｉ←ｉ＋１に更新する（図２２のステップＳ１１６）。
【０１６４】
ｉ≧ｍまたはｚ＋ｉ≧Ｚｓｔの場合、制御部１１は、ｚ＝ｚ＋ｉ−１を返し（図２２のステップＳ１１８）、処理を終了する。
【０１６５】
ｉ＜ｍかつｚ＋ｉ＜Ｚｓｔの場合、３次元座標（ｘ、ｙ、ｚ＋ｉ）に対して座標変換を行い、ボクセルのα値を取得する（図２２のステップＳ１１７）。座標変換を行い、ボクセルのα値を取得する処理は後述する（「１１−３．ボクセルα値を取得（補間なし）」参照）。
【０１６６】
ステップＳ１１７において取得したボクセルのα値がα＝０の場合（ボクセルが透明の場合）、制御部１１は、ｚ＝ｚ＋ｉ−１を返し（図２２のステップＳ１１８）、処理を終了する。
【０１６７】
ステップＳ１１７において取得したボクセルのα値がα＞０の場合（ボクセルが不透明の場合）、ステップＳ１１６に戻り、ステップＳ１１６〜Ｓ１１７の処理を繰り返す。
【０１６８】
以上のように、先頭の不透明ボクセルのＺ座標ｚ（仮想光線の照射を開始する起点座標）を所定のステップ幅で探索することで（ステップＳ１１２、Ｓ１１３）、起点座標を高速に特定できる。また、不透明ボクセルが探索された場合、当該ボクセルの位置とステップ幅だけ視線逆方向（Ｚ軸の＋方向）に戻った位置との間に存在する可能性がある不透明ボクセルを再探索することで（ステップＳ１１５〜Ｓ１１７）、起点座標を正確に特定することができる。
【０１６９】
（１１−３．ボクセルα値を取得（補間なし））
起点座標探索処理（図２２）のステップＳ１１２およびステップＳ１１７において実行される、座標変換を行いボクセルのα値を取得する処理について説明する。
【０１７０】
（１１−３−１．座標変換（実数の座標値を取得））
まず座標変換について説明する。
座標変換は、視点座標系をボクセル座標系に変換する処理であり、ＧＵＩ側の変換処理とは逆になる。ＧＵＩ側では関心領域ＲＯＩによるクリッピング、スケーリング、Ｚ軸方向変倍処理、オフセット（ＸＹＺ軸方向同時）、回転、透視変換の順に行うものと仮定し、制御部１１は、与えられた視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を以下のように算出する。
【０１７１】
制御部１１は、視点座標系の座標値（ｘ、ｙ、ｚ）（整数値）を次のように実数値（ｘｘ、ｙｙ、ｚｚ）に変換する。視点座標系はＳｉｚｅ（Ｘ軸方向のサイズ）×Ｓｉｚｅ（Ｙ軸方向のサイズ）×Ｓｉｚｅ／Ｓｒａｙ（Ｚ軸方向のサイズ）の直方体で、Ｓｘ×Ｓｙ×Ｓｚの断層画像群Ｄｏとは独立して定義される。通常Ｓｘ＝Ｓｙ＝５１２であればＳｉｚｅ＝５１２となる。視点座標系のＺ軸方向は仮想光線サブサンプル１／Ｓｒａｙ倍に伸縮されていることを考慮して、はじめに以下のようにオフセット処理（ＸＹＺ軸方向同時）を行う。
【０１７２】
ｘｘ＝ｘ−Ｓｉｚｅ／２＋ｄｘ＋Ｘｏｆｆ
ｙｙ＝ｙ−Ｓｉｚｅ／２＋ｄｙ＋Ｙｏｆｆ
ｚｚ＝ｚ・Ｓｒａｙ−Ｓｉｚｅ／２＋ｄｚ＋Ｚｏｆｆ
【０１７３】
続いて、制御部１１は、以下のように、回転パラメータ行列Ｒを用いて回転処理を行う。
【０１７４】
ｘｘ’＝Ｒ１１・ｘｘ＋Ｒ１２・ｙｙ＋Ｒ１３・ｚｚ
ｙｙ’＝Ｒ２１・ｘｘ＋Ｒ２２・ｙｙ＋Ｒ２３・ｚｚ
ｚｚ’＝Ｒ３１・ｘｘ＋Ｒ３２・ｙｙ＋Ｒ３３・ｚｚ
回転処理後の（ｘｘ’、ｙｙ’、ｚｚ’）を（ｘｘ、ｙｙ、ｚｚ）とする。
【０１７５】
制御部１１は、スケーリング、Ｚ軸方向変倍処理を同時に行い、以下のように、ボクセルの座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）を取得する。
【０１７６】
ｘｒ＝ｘｘ／Ｓｃａｌｅ＋Ｓｘ／２
ｙｒ＝ｙｙ／Ｓｃａｌｅ＋Ｓｙ／２
ｚｒ＝ｚｚ／Ｓｃａｌｅ／Ｓｃｚ＋Ｓｚ／２
【０１７７】
（１１−３−２．座標変換（整数の座標値に変換））
続いて、制御部１１は、算出したボクセルの座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）に対して、各値に０．５を加算して、以下のように、小数点以下を切り捨て整数化した座標値を（ｘｉ、ｙｉ、ｚｉ）（四捨五入した整数値）を求める。
【０１７８】
ｘｉ＝ＩＮＴ［ｘｒ＋０．５］
ｙｉ＝ＩＮＴ［ｙｒ＋０．５］
ｚｉ＝ＩＮＴ［ｚｒ＋０．５］
【０１７９】
（１１−３−３．ボクセルα値を取得）
そして、制御部１１は、有効ボクセル領域Ｖｒ、マスクデータＭａｓｋ、または補正テーブルＳαが定義されている場合、これらを考慮して、以下のようにボクセルα値を取得する。
【０１８０】
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外の場合）
α＝−１
２）上記１）を満たさない場合
α＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ、ｚｉ）、３）・Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）・Ｓα（ｘｉ、ｙｉ、ｚｉ）
【０１８１】
（１１−４．ボクセルα値、ＲＧＢ値の取得（補間あり））
図１９のステップＳ９６、Ｓ９８において実行される、３次元座標値（ｘ、ｙ、ｚ）に対して座標変換を行い、ボクセルα値（Ｖα）、ボクセルＲＧＢ値を取得する処理について説明する。
【０１８２】
（１１−４−１．座標変換）
まず、前記した（１１−３−１．座標変換（実数の座標値を取得））と同様の方法で、３次元座標値（ｘ、ｙ、ｚ）に対応するボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を算出する。
【０１８３】
続いて、制御部１１は、算出したボクセルの座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）に対して、小数点以下を切り捨て整数化した座標値を（ｘｉ、ｙｉ、ｚｉ）（整数値）とし、切り捨てた小数点以下の端数を（ｗｘ、ｗｙ、ｗｚ）（０≦ｗｘ、ｗｙ、ｗｚ＜１）とする（すなわち、ｘｒ＝ｘｉ＋ｗｘ、ｙｒ＝ｙｉ＋ｗｙ、ｚｒ＝ｚｉ＋ｗｚ）。
【０１８４】
（１１−４−２．ボクセル信号値の抽出）
そして、制御部１１は、有効ボクセル領域Ｖｒを考慮して、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセルの信号値を断層画像群Ｄｏに基づいて以下のように抽出する。
【０１８５】
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｄ０００＝−３２７６９（無効の値）
２）上記１）の条件を満たさない場合において、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（補間しない）
Ｄ０００＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
【０１８６】
３）上記１）２）の条件を満たさない場合（８近傍ボクセルの信号値を抽出）
Ｄ０００＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
Ｄ１００＝Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ）
Ｄ０１０＝Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ）
Ｄ１１０＝Ｄｏ（ｘｉ＋１、ｙｉ＋１、ｚｉ）
Ｄ００１＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ＋１）
Ｄ１０１＝Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ＋１）
Ｄ０１１＝Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ＋１）
Ｄ１１１＝Ｄｏ（ｘｉ＋１、ｙｉ＋１、ｚｉ＋１）
【０１８７】
ここで、マスクデータＭａｓｋ、補正テーブルＳαが定義されている場合は、制御部１１は、抽出した補間対象ボクセルの信号値に対応するα値スケールＳ０００〜Ｓ１１１（０≦Ｓ０００、・・・、Ｓ１１１≦１）を以下のように設定する。なお、マスクデータＭａｓｋ、補正テーブルＳαを使用しない場合は、Ｓ０００＝Ｓ１００＝Ｓ０１０＝Ｓ１１０＝Ｓ００１＝Ｓ１０１＝Ｓ０１１＝Ｓ１１１＝１とする。
【０１８８】
１）１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｓ０００＝０
２）上記１）の条件を満たさない場合において、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（補間しない）
Ｓ０００＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）・Ｓα（ｘｉ、ｙｉ、ｚｉ）
【０１８９】
３）上記１）２）の条件を満たさない場合（８近傍ボクセルの信号値を抽出）
Ｓ０００＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）・Ｓα（ｘｉ、ｙｉ、ｚｉ）
Ｓ１００＝Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ）・Ｓα（ｘｉ＋１、ｙｉ、ｚｉ）
Ｓ０１０＝Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ）・Ｓα（ｘｉ、ｙｉ＋１、ｚｉ）
Ｓ１１０＝Ｍａｓｋ（ｘｉ＋１、ｙｉ＋１、ｚｉ）・Ｓα（ｘｉ＋１、ｙｉ＋１、ｚｉ）
Ｓ００１＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ＋１）・Ｓα（ｘｉ、ｙｉ、ｚｉ＋１）
Ｓ１０１＝Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ＋１）・Ｓα（ｘｉ＋１、ｙｉ、ｚｉ＋１）
Ｓ０１１＝Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ＋１）・Ｓα（ｘｉ、ｙｉ＋１、ｚｉ＋１）
Ｓ１１１＝Ｍａｓｋ（ｘｉ＋１、ｙｉ＋１、ｚｉ＋１）・Ｓα（ｘｉ＋１、ｙｉ＋１、ｚｉ＋１）
【０１９０】
（１１−４−３．ボクセルα値の取得）
そして、制御部１１は、抽出した補間対象ボクセルの信号値に基づいて、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセルのα値Ｖαを以下のように決定する。
【０１９１】
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｖα＝−１
２）上記１）の条件を満たさず、ｘｉ＝Ｘｉｓ、ｘｉ＝Ｘｉｅ、ｙｉ＝Ｙｉｓ、ｙｉ＝Ｙｉｅ、ｚｉ＝Ｚｉｓ、又はｚｉ＝Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒの境界面）
Ｖα＝０
【０１９２】
３）上記１）２）の条件を満たさず、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（補間しない）
Ｖα＝Ｃｍａｐ（Ｄ０００、３）・Ｓ０００
【０１９３】
４）上記１）２）３）の条件を満たさない場合（補間する）
Ｖα＝（１−ｗｚ）（１−ｗｙ）（１−ｗｘ）・Ｃｍａｐ（Ｄ０００、３）・Ｓ０００＋（１−ｗｚ）（１−ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１００、３）・Ｓ１００＋（１−ｗｚ）・ｗｙ・（１−ｗｘ）・Ｃｍａｐ（Ｄ０１０、３）・Ｓ０１０＋（１−ｗｚ）・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１０、３）・Ｓ１１０＋ｗｚ・（１−ｗｙ）（１−ｗｘ）・Ｃｍａｐ（Ｄ００１、３）・Ｓ００１＋ｗｚ・（１−ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１０１、３）・Ｓ１０１＋ｗｚ・ｗｙ・（１−ｗｘ）・Ｃｍａｐ（Ｄ０１１、３）・Ｓ０１１＋ｗｚ・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１１、３）・Ｓ１１１
【０１９４】
（１１−４−４．ボクセルＲＧＢ値の取得）
また、制御部１１は、抽出した補間対象ボクセルの信号値に基づいて、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセルのＲＧＢ値（Ｖｃ（ｎ））（０≦ｎ≦２）を以下のように決定する。なお、ＲＧＢ値を取得する場合には、Ｓ０００、Ｓ１００、・・・、Ｓ１１１を乗算するマスク処理を行わない。これにより、マスク境界面でのモアレの発生を抑制する。
【０１９５】
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｖｃ（ｎ）＝０  （０≦ｎ≦２）
２）上記１）の条件を満たさず、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（補間しない）
Ｖｃ（ｎ）＝Ｃｍａｐ（Ｄ０００、ｎ）  （０≦ｎ≦２）
【０１９６】
３）上記１）２）の条件を満たさない場合（補間する）
Ｖｃ（ｎ）＝（１−ｗｚ）（１−ｗｙ）（１−ｗｘ）・Ｃｍａｐ（Ｄ０００、ｎ）＋（１−ｗｚ）（１−ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１００、ｎ）＋（１−ｗｚ）・ｗｙ・（１−ｗｘ）・Ｃｍａｐ（Ｄ０１０、ｎ）＋（１−ｗｚ）・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１０、ｎ）＋ｗｚ・（１−ｗｙ）（１−ｗｘ）・Ｃｍａｐ（Ｄ００１、ｎ）＋ｗｚ・（１−ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１０１、ｎ）＋ｗｚ・ｗｙ・（１−ｗｘ）・Ｃｍａｐ（Ｄ０１１、ｎ）＋ｗｚ・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１１、ｎ）  （０≦ｎ≦２）
【０１９７】
（１１−５．陰影計算）
制御部１１は、必要に応じて、次のように陰影計算を行ってもよい。
まず、光源ベクトル（Ｌｘ、Ｌｙ、Ｌｚ）（単位ベクトル）を設定する。例えば、（Ｌｘ、Ｌｙ、Ｌｚ）＝（０．５７７３５、０．５７７３５、０．５７７３５）と設定する。また、環境光成分Ａｂ（０≦Ａｂ≦１、例えばＡｂ＝０．２）を設定する。
【０１９８】
続いて、制御部１１は、次のように、座標変換により算出された３次元座標（ｘｉ、ｙｉ、ｚｉ）における６近傍ボクセルの不透明度Ｖ１００、Ｖ２００、Ｖ０１０、Ｖ０２０、Ｖ００１、Ｖ００２を取得する。但し、ｘｉ＋１＞ＸｉｅのときＶ１００＝０、ｘｉ−１＜ＸｉｓのときＶ２００＝０、ｙｉ＋１＞ＹｉｅのときＶ０１０＝０、ｙｉ−１＜ＹｉｓのときＶ０２０＝０、ｚｉ＋１＜ＺｉｓのときＶ００１＝０、ｚｉ−１＜ＺｉｓのときＶ００２＝０とする。すなわち、隣接するボクセルがクリッピング範囲の場合は不透明度を０として扱う。
【０１９９】
Ｖ１００＝Ｃｍａｐ（Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ）、３）・Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ）・Ｓα（ｘｉ＋１、ｙｉ、ｚｉ）
Ｖ２００＝Ｃｍａｐ（Ｄｏ（ｘｉ−１、ｙｉ、ｚｉ）、３）・Ｍａｓｋ（ｘｉ−１、ｙｉ、ｚｉ）・Ｓα（ｘｉ−１、ｙｉ、ｚｉ）
Ｖ０１０＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ）、３）・Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ）・Ｓα（ｘｉ、ｙｉ＋１、ｚｉ）
Ｖ０２０＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ−１、ｚｉ）、３）・Ｍａｓｋ（ｘｉ、ｙｉ−１、ｚｉ）・Ｓα（ｘｉ、ｙｉ−１、ｚｉ）
Ｖ００１＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ、ｚｉ＋１）、３）・Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ＋１）・Ｓα（ｘｉ、ｙｉ、ｚｉ＋１）
Ｖ００２＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ、ｚｉ−１）、３）・Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ−１）・Ｓα（ｘｉ、ｙｉ、ｚｉ−１）
【０２００】
上記のように、各不透明度の要素にＭａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ）・Ｓα（ｘｉ＋１、ｙｉ、ｚｉ）を乗算すると、カット面の反射に豊かな陰影が付くがモアレが生じやすい。その場合、Ｖ１００＝Ｃｍａｐ（Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ）、３）のようにＭａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ）・Ｓα（ｘｉ＋１、ｙｉ、ｚｉ）の項を乗算しなければモアレを抑圧できるが、カット面の光沢が無くなりザラザラになる弊害がある。
【０２０１】
続いて、制御部１１は、座標変換により算出された３次元座標（ｘｉ、ｙｉ、ｚｉ）における勾配ベクトル（Ｇｘ、Ｇｙ、Ｇｚ）を、以下の式で算出する。
【０２０２】
Ｇｘ＝（Ｖ１００−Ｖ２００）・Ｒｘｙ／Ｒｚ
Ｇｙ＝（Ｖ０１０−Ｖ０２０）・Ｒｘｙ／Ｒｚ
Ｇｚ＝Ｖ００１−Ｖ００２
Ｇ＝｛Ｇｘ２＋Ｇｙ２＋Ｇｚ２｝１／２
【０２０３】
続いて、Ｇ≧１の場合、制御部１１は、輝度値（陰影値）Ｓ（０≦Ｓ≦１）は、拡散反射成分を算出し、
Ｓ＝（１−Ａｂ）｜Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ｜／Ｇ＋Ａｂ
を与える。Ｇ＜１の場合（αが変化しない場合）、制御部１１は、輝度値（陰影値）Ｓとして、Ｓ＝０を与える。
【０２０４】
そして、制御部１１は、
Ｖｃ’（ｎ）＝Ｓ・Ｖｃ（ｎ）
とし、算出されたＲＧＢ値Ｖｃ（ｎ）（０≦ｎ≦２）の成分を改変する。
【０２０５】
以上、レイキャスティング法に基づいてボリュームレンダリング画像を生成する処理について説明した。
【０２０６】
（１２．ＭＩＰ画像の生成）
次に、レンダリング処理として、レイキャスティング法に基づいてＭＩＰ画像を生成する処理を説明する。
【０２０７】
図２３は、レイキャスティング法によるＭＩＰ画像を生成する処理の全体の流れを示すフローチャートである。この処理は、前記したように、投影面を分割した各画素領域Ｐ（０）〜Ｐ（３）に対して並行して実行される。
【０２０８】
まず、制御部１１は、座標変換パラメータを設定する（図２３のステップＳ１２１）。座標変換パラメータの設定は、前記した図１８のステップＳ８１と同様である。
【０２０９】
続いて、制御部１１は、Ｓｉｚｅ×（Ｙ２−Ｙ１＋１）サイズ（投影面をＹ軸方向に４分割（等分割）したサイズ）の画素領域に対して、レイキャスティング処理（画素領域の各座標（ｘ、ｙ）毎に代表信号値を算出する処理）を実行する（図２３のステップＳ１２２）。
【０２１０】
図２４のフローチャートを参照して、このレイキャスティング処理について説明する。生成するＭＩＰ画像Ｉｍａｇｅ（ｘ、ｙ）は断層画像群Ｄｏの階調に基づいてモノクロ８ビットまたは１６ビットの値をもつが、図１０のステップＳ２に示される断層画像群Ｄｏに対する階調圧縮処理を施さない場合が多く、その場合はモノクロ１６ビットになる。以下、モノクロ（１６ビット）のＭＩＰ画像Ｉｍａｇｅ（ｘ、ｙ）を生成する場合について説明する。
制御部１１は、まず、生成するモノクロ（１６ビット）のＭＩＰ画像Ｉｍａｇｅ（ｘ、ｙ）の初期値を全て０に設定する（Ｉｍａｇｅ（ｘ、ｙ）＝０）。また、代表信号値算出モードｍｏｄｅ（０：ＭＩＰ、１：ＭｉｎＩＰ、２：ＲａｙＳｕｍ）を設定し、各２次元座標（ｘ、ｙ）（０≦ｘ≦Ｓｉｚｅ−１、０≦ｙ≦Ｓｉｚｅ−１）に対して、以下の処理を実行する。
【０２１１】
制御部１１は、ｘ＝０、ｙ＝Ｙ１と設定し（図２４のステップＳ１３１）、有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する（図２４のステップＳ１３２）。交点のＺ座標を算出する方法は、前記した図１９のステップＳ９３（「１１−１．交点算出処理」参照）と同様である。
【０２１２】
続いて、制御部１１は、代表信号値算出モードｍｏｄｅに応じて、代表信号値Ｖｍとカウンタｃｎｔを以下のように初期化し、Ｚｓｔ＝Ｚｃとする（図２４のステップＳ１３３）。
【０２１３】
ｍｏｄｅ＝０（ＭＩＰ）：Ｖｍ＝−３２７６８（最小値を設定）
ｍｏｄｅ＝１（ＭｉｎＩＰ）：Ｖｍ＝３２７６７（最大値を設定）
ｍｏｄｅ＝２（ＲａｙＳｕｍ）：Ｖｍ＝０、ｃｎｔ＝０
【０２１４】
続いて、制御部１１は、Ｚｓｔから後続する先頭の有効ボクセル（信号値が所定の閾値以上のボクセル）のＺ座標ｚ（起点座標）を探索する（図２４のステップＳ１３４）。すなわち、ステップＳ１３２において算出された交点から視線方向に向かってレイキャスティング計算を開始する起点座標を探索する。起点座標を探索する処理については後述する（「１２−１．起点座標探索処理」、図２５参照）。
【０２１５】
ｚ＜０の場合、ステップＳ１３９へ移行する。
【０２１６】
一方、ｚ≧０の場合、制御部１１は、３次元座標（ｘ、ｙ、ｚ）を座標変換してボクセル信号値Ｖｓを取得する（図２４のステップＳ１３５）。座標変換を行い、ボクセル信号値Ｖｓを取得する処理については後述する（「１２−３．ボクセル信号値Ｖｓを取得（補間あり）」参照）。
【０２１７】
ステップＳ１３５において取得したボクセル信号値ＶｓがＶｓ＝−３２７６９の場合（有効ボクセル領域Ｖｒを通過した場合）、ステップＳ１３９へ移行する。これにより、有効ボクセル領域Ｖｒを通過している場合（Ｖｓ＝−３２７６９）は、その先に描画対象のボクセルは存在しないため、レイキャスティング処理を早期に打ち切り、冗長な処理を省略する。
【０２１８】
ステップＳ１３５において取得したボクセル信号値ＶｓがＶｓ＜−３２７６９の場合（無効値の場合）、Ｚｓｔ＝ｚ−１としたうえで（図２４のステップＳ１３６）、ステップＳ１３４に戻り、起点座標を再度探索する。
【０２１９】
ステップＳ１３５において取得したボクセル信号値ＶｓがＶｓ≧−３２７６８の場合（有効値の場合）、制御部１１は、代表信号値算出モードｍｏｄｅに応じて、以下のように代表信号値Ｖｍを算出する（図２４のステップＳ１３７）。
【０２２０】
ｍｏｄｅ＝０（ＭＩＰ  ｍｏｄｅ）：
Ｖｓ＞ＶｍならばＶｍ＝Ｖｓ
ｍｏｄｅ＝１（ＭｉｎＩＰ  ｍｏｄｅ）：
Ｖｓ＜ＶｍならばＶｍ＝Ｖｓ
ｍｏｄｅ＝２（ＲａｙＳｕｍ  ｍｏｄｅ）：
Ｖｍ←Ｖｍ＋Ｖｓ、ｃｎｔ←ｃｎｔ＋１
【０２２１】
ｚ＝ｚ−１に更新し（図２４のステップＳ１３８）、ｚ≧０の場合、ステップＳ１３５に戻り、ｚ＜０の場合、ステップＳ１３９へ移行する。
【０２２２】
ステップ１３９において、ｍｏｄｅ＝０、１の場合（図２４のステップＳ１３９；ｍｏｄｅ＜２）には、制御部１１は、ステップＳ１３７で得られた代表信号値Ｖｍを画素値として記録する（Ｉｍａｇｅ（ｘ、ｙ）＝Ｖｍ、図２４のステップＳ１４１）。一方、ｍｏｄｅ＝２（ＲａｙＳｕｍ  ｍｏｄｅ）の場合（図２４のステップＳ１３９；ｍｏｄｅ＝２）には、制御部１１は、Ｖｍ←Ｖｍ／ｃｎｔのように信号値の平均を計算したうえで（図２４のステップＳ１４０）、代表信号値Ｖｍを画素値として記録する（Ｉｍａｇｅ（ｘ、ｙ）＝Ｖｍ、図２４のステップＳ１４１）。
【０２２３】
続いて、制御部１１は、ｘ←ｘ＋１に更新し（図２４のステップＳ１４２）、ｘ＜Ｓｉｚｅの場合、ステップＳ１３２に戻り、次の画素ｘの代表信号値Ｖｍを求める処理（ステップＳ１３２〜Ｓ１４１）を実行する。ｘ≧Ｓｉｚｅの場合、ｙ←ｙ＋１に更新し（図２４のステップＳ１４３）、ｙ≦Ｙ２の場合、ｘ＝０としたうえで、ステップＳ１３２に戻り、ｙ行目の画素（ｘ、ｙ）の代表信号値Ｖｍを求める処理（ステップＳ１３２〜Ｓ１４１）を繰り返す。ｙ＞Ｙ２の場合、制御部１１は、処理を終了する。
すなわち、画素領域の全画素の代表信号値Ｖｍが得られるまで（ステップＳ１４２；ｘ≧Ｓｉｚｅ、かつ、ステップＳ１４３；ｙ＞Ｙ２）、ステップＳ１３２〜Ｓ１４１の処理を繰り返す。
【０２２４】
図２３のフローチャートに戻る。
制御部１１は、画素領域に対応するＭＩＰ画像Ｉｍａｇｅ（ｘ、ｙ）を出力する（図２３のステップＳ１２３）。
【０２２５】
（１２−１．起点座標探索処理）
図２５のフローチャートを参照して、図２４のステップＳ１３４において実行される、起点座標を探索する処理について説明する。
【０２２６】
まず、制御部１１は、ｚをｚ＝Ｚｓｔと初期化し、探索対象画素（ｘ、ｙ）を入力する（図２５のステップＳ１５１）。続いて、制御部１１は、３次元座標（ｘ、ｙ、ｚ）に対して座標変換を行い、ボクセル信号値Ｖｓを取得する（図２５のステップＳ１５２）。座標変換を行い、ボクセル信号値Ｖｓを取得する処理は後述する（「１２−２．ボクセル信号値Ｖｓを取得（補間なし）」参照）。
【０２２７】
ステップＳ１５２において取得したボクセル信号値ＶｓがＶｓ＝−３２７６９の場合（有効ボクセル領域Ｖrを通過した場合）、制御部１１は、ｚ＝−１（起点座標が存在しないことを示す値）を返し（図２５のステップＳ１５４）、処理を終了する。
【０２２８】
一方、ステップＳ１５２において取得したボクセル信号値Ｖｓ＜−３２７６９の場合（無効値の場合）は、制御部１１は、ｚ←ｚ−１に更新する（図２５のステップＳ１５３）。
【０２２９】
ステップＳ１５３において更新したzがｚ＜０の場合（ボクセル空間外の場合）、制御部１１は、ｚ＝−１（起点座標が存在しないことを示す値）を返し（図２５のステップＳ１５４）、処理を終了する。ｚ≧０の場合、制御部１１は、ステップＳ１５２に戻る。
【０２３０】
一方、ステップＳ１５２において取得したボクセル信号値ＶｓがＶｓ≧−３２７６８の場合（有効ボクセルの場合）、制御部１１は、ｚを起点座標として出力する。
【０２３１】
（１２−２．ボクセル信号値Ｖｓを取得（補間なし））
起点座標探索処理（図２５）のステップＳ１５２において実行される、座標変換を行いボクセル信号値Ｖｓを取得する処理について説明する。
【０２３２】
（１２−２−１．座標変換）
まず、前記した（１１−３−１．座標変換（実数の座標値を取得））および（１１−３−２．座標変換（整数の座標値に変換））と同様の方法で、３次元座標値（ｘ、ｙ、ｚ）に対応するボクセルの整数の座標値（ｘｉ、ｙｉ、ｚｉ）を算出する。
【０２３３】
（１２−２−２．ボクセル信号値Ｖｓを取得）
そして、制御部１１は、有効ボクセル領域Ｖｒ、マスクデータＭａｓｋ、または補正テーブルＳαが定義されている場合、これらを考慮して、以下のようにボクセル信号値Ｖｓを取得する。
【０２３４】
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外の場合）
Ｖｓ＝−３２７６９
２）上記１）を満たさず、Ｍａｓｋ（ｘ、ｙ、ｚ）＝０の場合（無効値）
Ｖｓ＝−９９９９９
３）上記１）２）を満たさない場合（有効値）
Ｖｓ＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
【０２３５】
（１２−３．ボクセル信号値Ｖｓの取得（補間あり））
（１２−３−１．座標変換）
まず、前記した（１１−３−１．座標変換（実数の座標値を取得））と同様の方法で、３次元座標値（ｘ、ｙ、ｚ）に対応するボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を算出する。
【０２３６】
（１２−３−２．ボクセル信号値の抽出）
続いて、前記した（１０−２．ボクセル信号値の抽出）と同様の方法で、補間対象のボクセル信号値Ｄ０００、Ｄ１００、Ｄ０１０、Ｄ００１、Ｄ１０１、Ｄ０１１、Ｄ１１１を抽出し、α値スケールＳ０００、Ｓ１００、Ｓ０１０、Ｓ１１０、Ｓ００１、Ｓ１０１、Ｓ０１１、Ｓ１１１を設定する。
【０２３７】
（１２−３−３．ボクセル信号値Ｖｓの取得）
そして、制御部１１は、抽出した補間対象ボクセルの信号値に基づいて、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセル信号値Ｖｓを以下のように決定する。
【０２３８】
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｖｓ＝−３２７６９
２）上記１）の条件を満たさず、Ｓ０００＝０の場合（無効値）
Ｖｓ＝−９９９９９
３）上記１）２）を満たさない場合において、
ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、またはｚｉ＋１＞Ｚｉｅのいずれかを満たすか、或いはＳ０００、Ｓ１００、Ｓ０１０、Ｓ１１０、Ｓ００１、Ｓ１０１、Ｓ０１１、Ｓ１１１のいずれかが０の場合（補間しない）
Ｖｓ＝Ｄ０００
【０２３９】
４）上記１）２）３）の条件を満たさない場合（補間する）
Ｖｓ＝（１−ｗｚ）（１−ｗｙ）（１−ｗｘ）・Ｄ０００＋（１−ｗｚ）（１−ｗｙ）・ｗｘ・Ｄ１００＋（１−ｗｚ）・ｗｙ・（１−ｗｘ）・Ｄ０１０＋（１−ｗｚ）・ｗｙ・ｗｘ・Ｄ１１０＋ｗｚ・（１−ｗｙ）（１−ｗｘ）・Ｄ００１＋ｗｚ・（１−ｗｙ）・ｗｘ・Ｄ１００＋ｗｚ・ｗｙ・（１−ｗｘ）・Ｄ０１１＋ｗｚ・ｗｙ・ｗｘ・Ｄ１１１
【０２４０】
以上、レンダリング処理として、レイキャスティング法に基づいてＭＩＰ画像を生成する処理について説明した。
【０２４１】
以上、本開示の実施の形態について説明した。本開示の実施の形態によれば、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する。そして、視線方向のベクトルと有効ボクセル領域との交点を幾何計算により算出し、算出された交点から視線方向に向かってレイキャスティング計算を開始する起点座標を探索する。これにより、有効ボクセル領域の外側における冗長な探索処理を削減することができ、レンダリング処理の高速化が実現できる。例えば、有効ボクセル領域との交点が存在しない背景部の画素は計算対象から外され、冗長な処理が大幅に削減される。また、背景部以外の画素についても、交点から起点座標の探索を開始すればよいため、投影面から交点までの範囲の探索を省略できる。本開示の方法は、レイキャスティング計算を開始する起点座標を探索するにあたり、投影面から有効ボクセルとの距離の分布を示すデプスマップを作成したりする必要がないので、レンダリング処理の高速化を簡便な手段で実現することができる。
【０２４２】
また、起点座標の探索範囲が絞り込まれるので、従来に比べて起点座標を探索する際のスキップ幅を小さくすることができる。これにより、従来に比べて起点座標をより正確に求めることができ、精度の高いレンダリング画像を得ることができる。
【０２４３】
図２６〜図２８は、スキップ幅を変えた場合のレンダリング画像の画像精度を比較した図である。図２６はスキップ幅＝１６、図２７はスキップ幅＝８、図２８はスキップ幅＝４とした場合のレンダリング画像を示す。（ａ）はフルカラー（２４ビット）のボリュームレンダリング画像、（ｂ）はモノクロ（８ビット）のボリュームレンダリング画像とＭＩＰ画像を示す。なお、図２８（ｂ）のＭＩＰ画像のみスキップ幅を１としている。
【０２４４】
図２６、図２７のように、スキップ幅をある程度大きくするとレンダリング処理を高速に実行できるが、有効ボクセルの見落としが多くなるため、画素抜け（図２６（ａ）、図２７（ａ）参照）や画像エッジ部の波形歪み（図２６（ｂ）、図２７（ｂ）参照）が発生する。これに対し、図２８のようにスキップ幅を小さくすると、有効ボクセルの見落としが少なくなるため、画素抜け（図２８（ａ）参照）や画像エッジ部の波形歪み（図２８（ｂ）参照））の発生が抑えられる。
【０２４５】
しかしながら、従来（図３（ａ）参照）のように、投影面から起点座標を探索する方法でスキップ幅を小さくすると、計算負荷が大きくなり、現実的な処理時間でレンダリング処理を実行することが困難となる。この点、本開示（図３（ｂ）参照）のように、視線ベクトルと有効ボクセル領域との交点を算出し、この交点から起点座標を探索する方法をとることで、図２８のようにスキップ幅を小さく設定しても現実的な処理時間でレンダリング処理を実行することができる。すなわち、本開示の方法によれば、レンダリング処理を高速かつ高精度に実行することが可能となる。
【０２４６】
さらに本開示では、レイキャスティング計算の途上において、有効ボクセル領域を逸脱（通過）した場合（図１９のステップＳ９６；Ａｌｐｈａ＜０、図２４のステップＳ１３５；Ｖｓ＝−３２７６９）、その先に描画対象のボクセルは存在しないため、レイキャスティング処理を早期に打ち切るようにする。これにより、有効ボクセル領域の外側におけるレイキャスティング計算を省略できる。  
【０２４７】
また本開示では、断層画像群Ｄｏから構成されるボクセル空間を分割し、分割したボクセル空間ごとに有効ボクセル領域を並行して算出し、算出した各有効ボクセル領域を統合する方法をとることで高速に有効ボクセル領域を求めることができる。また、レンダリングの投影面を複数の画素領域に分割し、各画素領域におけるレンダリング処理を並行して実行させることで、高速にレンダリング画像を得ることができる。例えば、並行化のスレッド数をＮ（例えばＮ＝４）とし、レンダリング装置１で使用可能なＣＰＵコアがＮ個以上存在し、Ｎ個のスレッドに異なるＣＰＵコアが割り当てられれば、Ｎ倍の処理速度で有効ボクセル領域算出処理およびレンダリング処理を実行することが可能となる。
【０２４８】
また本開示では、上記した並行化処理を１台のコンピュータ（レンダリング装置１）で完結させることで、リモートデクストップ等の構成で複数のユーザが使用しても、サーバーコンピュータ間の通信が発生しないので、通信経路上でデータが喪失（データロス）したり、ネットワーク遅延等に起因するレンダリング処理の遅延が生じない。また、機密性の高い断層画像群Ｄｏ（ＤＩＣＯＭ画像）が複数のクライアントコンピュータにコピーされるようなことがないため、セキュリティ対策や個人情報の取り扱いが容易となる。
【０２４９】
なお、上記した並行化処理は必須ではない。この場合、図１０のステップＳ３において、制御部１１は、ボクセル空間を分割することなく、ボクセル空間全体における有効ボクセル領域を直接算出する。また、図１０のステップＳ４において、制御部１１は、投影面を分割することなく、投影面全体を対象としてレンダリング処理（レイキャスティング処理）を実行する。
【０２５０】
以上、添付図面を参照しながら、本開示に係るレンダリング装置１等の好適な実施形態について説明したが、本開示はかかる例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本開示の技術的範囲に属するものと了解される。
【符号の説明】
【０２５１】
１          ：レンダリング装置
Ｄｏ        ：断層画像群
Ｃｍａｐ    ：カラーマップ
ＲＯＩ      ：クリッピング領域
Ｍａｓｋ    ：マスクデータ
Ｓα        ：補正テーブル
Ｉｍａｇｅ  ：レンダリング画像(ボリュームレンダリング画像、ＭＩＰ画像)
Ｖｒ        ：有効ボクセル領域
【図１】
2020086803.tif 000003
【図２】
2020086803.tif 000004
【図３】
2020086803.tif 000005
【図４】
2020086803.tif 000006
【図５】
2020086803.tif 000007
【図６】
2020086803.tif 000008
【図７】
2020086803.tif 000009
【図８】
2020086803.tif 000010
【図９】
2020086803.tif 000011
【図１０】
2020086803.tif 000012
【図１１】
2020086803.tif 000013
【図１２】
2020086803.tif 000014
【図１３】
2020086803.tif 000015
【図１４】
2020086803.tif 000016
【図１５】
2020086803.tif 000017
【図１６】
2020086803.tif 000018
【図１７】
2020086803.tif 000019
【図１８】
2020086803.tif 000020
【図１９】
2020086803.tif 000021
【図２０】
2020086803.tif 000022
【図２１】
2020086803.tif 000023
【図２２】
2020086803.tif 000024
【図２３】
2020086803.tif 000025
【図２４】
2020086803.tif 000026
【図２５】
2020086803.tif 000027
【図２６】
2020086803.tif 000028
【図２７】
2020086803.tif 000029
【図２８】
2020086803.tif 000030