(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-10-31
(45)【発行日】2022-11-09
(54)【発明の名称】ボリュームレンダリング装置
(51)【国際特許分類】
G06T 15/08 20110101AFI20221101BHJP
【ＦＩ】
G06T15/08
【請求項の数】
16
(21)【出願番号】P 2018239791
(22)【出願日】2018-12-21
(65)【公開番号】P2020102004
(43)【公開日】2020-07-02
【審査請求日】2021-10-25
(73)【特許権者】
【識別番号】000002897
【氏名又は名称】大日本印刷株式会社
(74)【代理人】
【識別番号】100122529
【弁理士】
【氏名又は名称】藤枡 裕実
(74)【代理人】
【識別番号】100135954
【弁理士】
【氏名又は名称】深町 圭子
(74)【代理人】
【識別番号】100119057
【弁理士】
【氏名又は名称】伊藤 英生
(74)【代理人】
【識別番号】100131369
【弁理士】
【氏名又は名称】後藤 直樹
(74)【代理人】
【識別番号】100171859
【弁理士】
【氏名又は名称】立石 英之
(72)【発明者】
【氏名】茂出木 敏雄
【審査官】村松 貴士
(56)【参考文献】
【文献】特開２００６－１０７０９３（ＪＰ，Ａ）
【文献】特開２００９－１６０３０６（ＪＰ，Ａ）
【文献】特開平１１－１７５７４３（ＪＰ，Ａ）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｔ １５／０８
(57)【特許請求の範囲】
【請求項1】
所定の対象物について所定の間隔で撮影された複数の２次元の断層画像で構成されるボクセル画像に基づいて、信号値に対して色成分の値と不透明度を対応付けて定義されたカラーマップを参照してレンダリング画像を生成するためのボリュームレンダリング装置であって、 前記ボクセル画像と前記カラーマップを用いてレンダリング画像を生成するレンダリング手段を備え、
前記レンダリング手段は、前記レンダリング画像の各画素において、前記ボクセル画像に対して所定の座標変換を行いながら前記ボクセル画像の複数のボクセルを参照し、参照する各ボクセルに対して、前記カラーマップを参照しながら当該ボクセルの色成分および不透明度を取得するとともに、当該ボクセルの複数の近傍のボクセルの不透明度を取得するレイキャスティング手段と、
前記レイキャスティング手段により取得された当該ボクセルの複数の近傍のボクセルの不透明度を基に当該ボクセルにおける勾配ベクトルを算出するとともに、算出された勾配ベクトルに対する所定の光源ベクトルの正反射ベクトルを算出し、前記勾配ベクトルと前記光源ベクトルとの内積値を基に当該ボクセルの拡散反射光を算出し、前記正反射ベクトルと所定の視線ベクトルとの内積値を基に当該ボクセルの鏡面反射光を算出し、算出された拡散反射光と鏡面反射光を基に算出した陰影値を前記レイキャスティング手段により取得された当該ボクセルの色成分に乗算し、当該ボクセルの色成分を更新する陰影付加手段と、を有するボリュームレンダリング装置。
【請求項2】
前記所定の光源ベクトルは、あらかじめ定義された光源ベクトルに対して前記所定の座標変換における回転処理を行ったものであり、前記所定の視線ベクトルは、あらかじめ定義された視線ベクトルに対して前記所定の座標変換における回転処理を行ったものである請求項１に記載のボリュームレンダリング装置。
【請求項3】
前記回転処理は、所定の３次元座標系における回転を定義した３×３行列を用いて、前記あらかじめ定義された光源ベクトルまたは前記あらかじめ定義された視線ベクトルに対して回転処理を行う請求項２に記載のボリュームレンダリング装置。
【請求項4】
前記ボクセル画像のｘｙ軸方向の解像度をＲｘｙ、ｚ軸方向の解像度をＲｚとすると、前記陰影付加手段は、
当該ボクセルに対してｘ軸方向に＋１ボクセルおよび－１ボクセルずれたボクセルの不透明度の差分を前記勾配ベクトルのｘ軸方向成分とし、
当該ボクセルに対してｙ軸方向に＋１ボクセルおよび－１ボクセルずれたボクセルの不透明度の差分を前記勾配ベクトルのｙ軸方向成分とし、
当該ボクセルに対してｚ軸方向に＋１ボクセルおよび－１ボクセルずれたボクセルの不透明度の差分にＲｚ／Ｒｘｙを乗じた値をを前記勾配ベクトルのｚ軸方向成分としている請求項１から請求項３のいずれか一項に記載のボリュームレンダリング装置。
【請求項5】
前記ボクセル画像のｘｙ軸方向の解像度をＲｘｙ、ｚ軸方向の解像度をＲｚとすると、前記陰影付加手段は、
当該ボクセルに対してｘ軸方向に＋１ボクセルおよび－１ボクセルずれたボクセルの不透明度の差分を前記勾配ベクトルのｘ軸方向成分とし、
当該ボクセルに対してｙ軸方向に＋１ボクセルおよび－１ボクセルずれたボクセルの不透明度の差分を前記勾配ベクトルのｙ軸方向成分とし、
当該ボクセルに対して、ｚ軸方向に＋Ｒｚ／Ｒｘｙずれた座標の近傍の２つのボクセルの不透明度に所定の重みをかけて加算した不透明度と、ｚ軸方向に－Ｒｚ／Ｒｘｙずれた座標の近傍の２つのボクセルの不透明度に所定の重みをかけて加算した不透明度との差分を前記勾配ベクトルのｚ軸方向成分としている請求項１から請求項３のいずれか一項に記載のボリュームレンダリング装置。
【請求項6】
前記陰影付加手段は、あらかじめ所定の拡散反射係数，鏡面反射係数，鏡面光沢度，環境光を各々定義し、
前記正反射ベクトルは、前記勾配ベクトルとのなす角が前記光源ベクトルと前記勾配ベクトルとのなす角と同一になるように定義し、
前記勾配ベクトル、前記正反射ベクトル、前記光源ベクトル、前記視線ベクトルのいずれも単位ベクトルとするとき、
前記拡散反射光として、前記勾配ベクトルと前記光源ベクトルとの内積値に前記拡散反射係数を乗算した値とし、
前記鏡面反射光として、前記正反射ベクトルと前記視線ベクトルとの内積値に対して、前記鏡面光沢度で累乗した値に前記鏡面反射係数を乗算した値とし、
前記陰影値として、前記拡散反射光と前記鏡面反射光に前記環境光を加算した値としている請求項１から請求項５のいずれか一項に記載のボリュームレンダリング装置。
【請求項7】
前記レンダリング画像の生成領域をｘ軸方向にＮＢ（ＮＢ＞１）個の分割生成領域に分割し、前記レンダリング手段は、ＮＢ個の互いに異なる分割生成領域に対して、ＮＢ個の互いに異なるプロセッサで同時に実行させることにより、ＮＢ個の分割画像で構成されるレンダリング画像を生成するようにしている請求項１から請求項６のいずれか一項に記載のボリュームレンダリング装置。
【請求項8】
前記各断層画像の各画素のうち、前記カラ－マップを参照して得られた不透明度が０でない画素が存在する領域に３次元の各軸方向において外接する領域である不透明外接領域を算出する不透明外接領域算出手段と、を更に有し、前記レンダリング手段は、前記各断層画像のうち、前記不透明外接領域に対応する領域を有効領域として定義し、
前記レイキャスティング手段は、前記レンダリング画像の各画素において、前記有効領域と、視点に最も近い視線方向のベクトルとの交点を算出し、前記レンダリング画像の各画素において、前記算出される交点より前記ボクセル画像に対して視線方向に仮想光線を照射し、レイキャスティング法に基づいて所定の座標変換を行いながら前記ボクセル画像の複数のボクセルを参照する請求項１から請求項７のいずれか一項に記載のボリュームレンダリング装置。
【請求項9】
前記断層画像に対応して、０（表示しない）または１（表示する）の値をもつマスクデータが定義されている場合、前記不透明外接領域算出手段は、更に前記マスクデータを参照して、前記マスクデータの値が１で、かつ前記不透明度が０でない画素を特定するようにしている請求項８に記載のボリュームレンダリング装置。
【請求項10】
前記不透明外接領域算出手段は、所定の分割数をＮＡ（ＮＡ＞１）として、各前記断層画像をＮＡ個の断層画像群に分割し、ＮＡ個の各断層画像群に対して、前記不透明外接領域の一部を作成する処理を並行して実行する請求項８または請求項９に記載のボリュームレンダリング装置。
【請求項11】
２次元の前記各断層画像の１つ目の次元の最小の座標をＸｓｏ、１つ目の次元の最大の座標をＸｅｏ、２つ目の次元の最小の座標をＹｓｏ、２つ目の次元の最大の座標をＹｅｏ、第１番目の断層画像が配置されるｚ座標をＺｓｏ、第Ｍ番目の断層画像が配置されるｚ座標をＺｅｏとすると、Ｘｓｏ≦Ｘｓ＜Ｘｅ≦Ｘｅｏ、Ｙｓｏ≦Ｙｓ＜Ｙｅ≦Ｙｅｏ、Ｚｓｏ≦Ｚｓ＜Ｚｅ≦Ｚｅｏを満たす、６つの座標Ｘｓ，Ｘｅ，Ｙｓ，Ｙｅ，Ｚｓ，Ｚｅが定義されている場合、前記不透明外接領域算出手段は、Ｘｓ以上Ｘｅ以下、Ｙｓ以上Ｙｅ以下、Ｚｓ以上Ｚｅの範囲において前記不透明外接領域を算出する請求項８から請求項１０のいずれか一項に記載のボリュームレンダリング装置。
【請求項12】
前記レイキャスティング手段は、前記ボクセル画像を前記レンダリング画像に投影変換した座標系を視点座標系とすると、視点座標系において視線方向であるｚ軸の上限値を、前記レンダリング画像の画素（ｘ，ｙ）ごとに前記算出された交点に対応するＺｃに設定し、視点座標系のｚ軸上のＺｃから所定の下限値の範囲で仮想光線の照射を開始する起点座標Ｚｏを探索する起点座標探索手段を備えている請求項８から請求項１１のいずれか一項に記載のボリュームレンダリング装置。
【請求項13】
前記起点座標探索手段は、前記視点座標系において、前記レンダリング画像の各画素（ｘ，ｙ）におけるｚ座標を探索開始座標Ｚｓｔより前記下限値まで所定数ずつ変化させながら、前記視点座標系の３次元座標（ｘ，ｙ，ｚ）ごとに前記所定の座標変換を行って対応する前記ボクセル画像のボクセルのボクセル値を順次取得し、最初に見つかった不透明なボクセルの視点座標系におけるｚ座標を、起点座標Ｚｏとして探索するようにしている請求項１２に記載のボリュームレンダリング装置。
【請求項14】
前記レイキャスティング手段は、前記レンダリング画像の各画素（ｘ，ｙ）に対して、前記起点座標探索手段より探索された起点座標よりｚ軸に沿って下限値に向けて、所定の光強度をもつ仮想光線を照射する際、前記視点座標系の３次元座標（ｘ，ｙ，ｚ）ごとに前記所定の座標変換を行って、対応する前記ボクセル画像のボクセルのボクセル値として前記カラーマップを参照しながら不透明度を取得し、取得された不透明度が０でない場合、ボクセル値として前記カラーマップを参照しながら色成分（ＲＧＢ）を更に取得し、
前記陰影付加手段は、当該ボクセルの複数の近傍のボクセルの不透明度を基に算出した陰影値を取得した色成分に乗算することにより前記色成分を更新し、
前記レンダリング手段は、
当該ボクセルの不透明度を基に前記光強度を減衰させるとともに、当該ボクセルの不透明度及び前記更新された色成分および前記減衰させた光強度に基づいて累積輝度値を算出する処理を行い、算出された累積輝度値を基に、前記レンダリング画像の当該画素（ｘ，ｙ）に対応する画素値（ＲＧＢ）として与えるようにしている請求項１２に記載のボリュームレンダリング装置。
【請求項15】
所定の３次元座標系における回転を定義した３×３行列、ｘｙｚ各軸方向のオフセット値および微小オフセット値、ｘｙｚ各軸方向の拡大又は縮小倍率、ｚ軸方向の変倍率を含む前記所定の座標変換のパラメータを取得し、前記視点座標系の整数値の３次元座標（ｘ，ｙ，ｚ）を、前記パラメータに基づいて前記ボクセル画像の座標系に変換を行って、前記ボクセル画像の実数値のボクセル座標を算出し、
算出した実数値のボクセル座標の近傍の複数の整数値のボクセル座標に対応する前記ボクセル画像の複数のボクセルを特定し、
特定した複数のボクセルの整数値のボクセル座標が前記有効領域に含まれる場合は、前記複数のボクセルの信号値を前記カラーマップを参照して色成分または不透明度に変換し、変換された色成分または不透明度の各々に所定の重みを掛けながら加算した値を、ボクセル値として色成分または不透明度を算出するようにしている請求項１２から請求項１４のいずれか一項に記載のボリュームレンダリング装置。
【請求項16】
請求項１から請求項１５のいずれか一項に記載のボリュームレンダリング装置として、コンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本開示は、２次元の断層画像を複数枚用いて、３次元的に可視化するためのボリュームレンダリング技術に関する。
【背景技術】
【0002】
従来、ＣＴ、ＭＲＩ、ＰＥＴなど医療画像診断機器（モダリティ）により所定のスライス間隔で連続的に撮像され、ＤＩＣＯＭ形式等でＰＡＣＳ（Picture Archiving and Communication Systems）等の医療情報システムに保管されている複数の断層画像を基に、臓器等を３次元的に可視化することが行われている。
【0003】
一方、３Ｄコンピュータグラフィックス分野においては、所定の方向に仮想的な光源を設定し、陰影を付加して、より立体的に表現することが行われている（特許文献１～３参照）。特許文献１では、リアルな陰影を付与するため、ボクセルの不透明度の勾配を基に臓器領域の境界面を抽出し、境界面における法線ベクトルに基づいて境界面に陰影を付加する手法が提案されている。
【0004】
また、特許文献２では、各ボクセルの信号値の勾配からグラディエントを算出し、各ボクセルに直接陰影を付与する手法が提案されている。また、特許文献３では、各ボクセルの信号値の勾配からグラディエントを算出し、各ボクセルに直接陰影を付与する手法が提案されている。
【先行技術文献】
【特許文献】
【0005】
【文献】特許第５０６５７４０号公報
【文献】特許第３５４２７９９号公報
【文献】特許第２６２７６０７号公報
【発明の概要】
【発明が解決しようとする課題】
【0006】
しかしながら、特許文献１に記載の技術では、各ボクセルに直接陰影を付与する方法に比べ、陰影が境界面に限定して付与されるため、読影しやすい画像が得られるという利点があるが、計算負荷が大きく、鏡面反射をサポートできないという問題がある。また、特許文献２、３に記載の技術では、グラディエントの算出を不透明度で行わず、信号値で直接行っているため、カラーマップで陰影を制御できず、鏡面反射をサポートできないという問題がある。
【0007】
そこで、本開示は、複数の断層画像に対してカラーマップを適用して、レンダリング画像を生成する際、計算負荷を抑えつつ陰影を付与することが可能なボリュームレンダリング装置を提供することを課題とする。
【課題を解決するための手段】
【0008】
本開示は、上記課題を解決する手段を複数含んでいるが、その一例を挙げるならば、
所定の対象物について所定の間隔で撮影された複数の２次元の断層画像で構成されるボクセル画像に基づいて、信号値に対して色成分の値と不透明度を対応付けて定義されたカラーマップを参照してレンダリング画像を生成するためのボリュームレンダリング装置であって、
前記ボクセル画像と前記カラーマップを用いてレンダリング画像を生成するレンダリング手段と、を備え、
前記レンダリング手段は、前記レンダリング画像の各画素において、前記ボクセル画像に対して所定の座標変換を行いながら前記ボクセル画像の複数のボクセルを参照し、参照する各ボクセルに対して、前記カラーマップを参照しながら当該ボクセルの色成分および不透明度を取得するとともに、当該ボクセルの複数の近傍のボクセルの不透明度を取得するレイキャスティング手段と、
前記レイキャスティング手段により取得された当該ボクセルの複数の近傍のボクセルの不透明度を基に当該ボクセルにおける勾配ベクトルを算出するとともに、算出された勾配ベクトルに対する所定の光源ベクトルの正反射ベクトルを算出し、前記勾配ベクトルと前記光源ベクトルとの内積値を基に当該ボクセルの拡散反射光を算出し、前記正反射ベクトルと所定の視線ベクトルとの内積値を基に当該ボクセルの鏡面反射光を算出し、算出された拡散反射光と鏡面反射光を基に算出した陰影値を前記レイキャスティング手段により取得された当該ボクセルの色成分に乗算し、当該ボクセルの色成分を更新する陰影付加手段と、を有する。
所定の対象物について所定の間隔で撮影された複数の２次元の断層画像で構成されるボクセル画像に基づいて、信号値に対して色成分の値と不透明度を対応付けて定義されたカラーマップを参照してレンダリング画像を生成するためのボリュームレンダリング装置であって、
前記ボクセル画像と前記カラーマップを用いてレンダリング画像を生成するレンダリング手段と、を備え、
前記レンダリング手段は、前記レンダリング画像の各画素において、前記ボクセル画像に対して所定の座標変換を行いながら前記ボクセル画像の複数のボクセルを参照し、参照する各ボクセルに対して、前記カラーマップを参照しながら当該ボクセルの色成分および不透明度を取得するとともに、当該ボクセルの複数の近傍のボクセルの不透明度を取得するレイキャスティング手段と、
前記レイキャスティング手段により取得された当該ボクセルの複数の近傍のボクセルの不透明度を基に当該ボクセルにおける勾配ベクトルを算出するとともに、算出された勾配ベクトルに対する所定の光源ベクトルの正反射ベクトルを算出し、前記勾配ベクトルと前記光源ベクトルとの内積値を基に当該ボクセルの拡散反射光を算出し、前記正反射ベクトルと所定の視線ベクトルとの内積値を基に当該ボクセルの鏡面反射光を算出し、算出された拡散反射光と鏡面反射光を基に算出した陰影値を前記レイキャスティング手段により取得された当該ボクセルの色成分に乗算し、当該ボクセルの色成分を更新する陰影付加手段と、を有する。
【発明の効果】
【0009】
本開示によれば、複数の断層画像に対してカラーマップを適用して、レンダリング画像を生成する際、計算負荷を抑えつつ陰影を付与することが可能となる。
【図面の簡単な説明】
【0010】
【図1】複数の断層画像と、ボリュームレンダリング像として生成されるカラーのレンダリング画像を示す図である。
【図2】本開示の一実施形態に係るボリュームレンダリング装置１００のハードウェア構成図である。
【図3】本実施形態に係るボリュームレンダリング装置の構成を示す機能ブロック図である。
【図4】本実施形態に係るボリュームレンダリング装置の処理動作を示すフローチャートである。
【図5】本実施形態で用いるカラーマップを示す図である。
【図9】ＲＯＩと不透明外接領域を比較した図である。
【図11】座標変換の概念図である。
【図15】コンピュータによりボリュームレンダリング装置１００を実現する場合の並列化処理を示す図である。
【図16】コンピュータによりボリュームレンダリング装置１００を実現する場合の並列化処理のソフトウェアによる流れを示す図である。
【発明を実施するための形態】
【0011】
以下、本開示の好適な実施形態について図面を参照して詳細に説明する。
本開示は、複数の断層画像に基づいて、信号値に対して色成分の値と不透明度を対応付けて定義されたカラ－マップを参照してレンダリング画像を生成するものである。図１は、複数の断層画像と、レンダリング画像を示す図である。図１（ａ）は複数の断層画像、図１（ｂ）は１枚の断層画像、図１（ｃ）はレンダリング画像である。図１（ａ）（ｂ）に示す断層画像は、ＤＩＣＯＭ形式のものを示している。図１（ａ）は、３７０枚の断層画像群であり、図１（ｂ）は、そのうちの１枚を拡大したものである。後述する実施形態のボリュームレンダリング装置では、図１（ａ）に示したような複数の断層画像に対してレンダリング処理を実行して、図１（ｃ）に示したような１枚のレンダリング画像を生成する。図１（ｃ）に示したレンダリング画像は、５１２×５１２画素であって、１画素についてＲＧＢ各色８ビットの計２４ビットで表現したものである。本実施形態では、２次元の断層画像を、もう１つの別の次元となる方向に積層した状態を３次元のボクセル画像として扱う。３次元のボクセル画像における各ボクセルは、複数の２次元の断層画像における各画素に１対１で対応している。
本開示は、複数の断層画像に基づいて、信号値に対して色成分の値と不透明度を対応付けて定義されたカラ－マップを参照してレンダリング画像を生成するものである。図１は、複数の断層画像と、レンダリング画像を示す図である。図１（ａ）は複数の断層画像、図１（ｂ）は１枚の断層画像、図１（ｃ）はレンダリング画像である。図１（ａ）（ｂ）に示す断層画像は、ＤＩＣＯＭ形式のものを示している。図１（ａ）は、３７０枚の断層画像群であり、図１（ｂ）は、そのうちの１枚を拡大したものである。後述する実施形態のボリュームレンダリング装置では、図１（ａ）に示したような複数の断層画像に対してレンダリング処理を実行して、図１（ｃ）に示したような１枚のレンダリング画像を生成する。図１（ｃ）に示したレンダリング画像は、５１２×５１２画素であって、１画素についてＲＧＢ各色８ビットの計２４ビットで表現したものである。本実施形態では、２次元の断層画像を、もう１つの別の次元となる方向に積層した状態を３次元のボクセル画像として扱う。３次元のボクセル画像における各ボクセルは、複数の２次元の断層画像における各画素に１対１で対応している。
【0012】
＜１．装置構成＞
図２は、本開示の一実施形態に係るボリュームレンダリング装置１００のハードウェア構成図である。本実施形態に係るボリュームレンダリング装置１００は、汎用のコンピュータで実現することができ、図２に示すように、ＣＰＵ（Central Processing Unit）１と、コンピュータのメインメモリであるＲＡＭ（Random Access Memory）２と、ＣＰＵ１が実行するプログラムやデータを記憶するためのハードディスク、ＳＳＤ（Solid State Drive），フラッシュメモリ等の大容量の記憶装置３と、キーボード、マウス等の指示入力Ｉ／Ｆ（インターフェース）４と、データ記憶媒体等の外部装置とデータ通信するためのデータ入出力Ｉ／Ｆ（インターフェース）５と、液晶ディスプレイ等の表示デバイスである表示部６と、グラフィックスに特化した演算処理部であるＧＰＵ（Graphics Processing Unit）７と、表示部６に表示する画像を保持するフレームメモリ８と、を備え、互いにバスを介して接続されている。ＧＰＵ７による演算結果はフレームメモリ８に書き込まれるため、ＧＰＵ７とフレームメモリ８は、表示部６へのインタフェースを備えたビデオカードに搭載されて汎用のコンピュータにバス経由で装着されていることが多い。
図２は、本開示の一実施形態に係るボリュームレンダリング装置１００のハードウェア構成図である。本実施形態に係るボリュームレンダリング装置１００は、汎用のコンピュータで実現することができ、図２に示すように、ＣＰＵ（Central Processing Unit）１と、コンピュータのメインメモリであるＲＡＭ（Random Access Memory）２と、ＣＰＵ１が実行するプログラムやデータを記憶するためのハードディスク、ＳＳＤ（Solid State Drive），フラッシュメモリ等の大容量の記憶装置３と、キーボード、マウス等の指示入力Ｉ／Ｆ（インターフェース）４と、データ記憶媒体等の外部装置とデータ通信するためのデータ入出力Ｉ／Ｆ（インターフェース）５と、液晶ディスプレイ等の表示デバイスである表示部６と、グラフィックスに特化した演算処理部であるＧＰＵ（Graphics Processing Unit）７と、表示部６に表示する画像を保持するフレームメモリ８と、を備え、互いにバスを介して接続されている。ＧＰＵ７による演算結果はフレームメモリ８に書き込まれるため、ＧＰＵ７とフレームメモリ８は、表示部６へのインタフェースを備えたビデオカードに搭載されて汎用のコンピュータにバス経由で装着されていることが多い。
【0013】
本実施形態において、ＣＰＵ１は、マルチコアＣＰＵであり、複数のＣＰＵコアを有し、並列処理が可能となっている。図２の例では、ＲＡＭ２が１つだけ示されているが、ＣＰＵ１の各ＣＰＵコアが、１つのＲＡＭ２にアクセスするように構成されている。なお、ＣＰＵ１は複数であってもよい。またマルチコアＣＰＵは、各々のＣＰＵコアが論理的に複数の実行スレッドを受け入れるＣＰＵであってもよい。すなわち、ＣＰＵ１は、複数の実行スレッドそれぞれを受け入れるプロセッサを有している。
【0014】
図３は、本実施形態に係るボリュームレンダリング装置の構成を示す機能ブロック図である。図３において、１０は断層画像読込手段、２０はカラーマップ読込手段、３０はＲＯＩクリッピング設定手段、４０はマスク設定手段、５０は不透明外接領域算出手段、６０はレンダリング手段、６１はレイキャスティング手段、６２は起点座標探索手段、６３は陰影付加手段、７０はレンダリング画像出力手段である。
【0015】
断層画像読込手段１０は、ＣＴ、ＭＲＩ、ＰＥＴなどの医用画像診断機器により収集および蓄積されたＰＡＣＳから、記憶媒体や入力部を経由して、複数の断層画像を読み込む手段である。断層画像は、対象物に対して所定の間隔で撮影されて得られたものであり、各画素に信号値が付与された２次元の断層画像である。カラーマップ読込手段２０は、図示されていないカラーマップ作成手段により、あらかじめ作成されたカラーマップを記憶媒体や入力部から、カラーマップを読み込む手段である。ＲＯＩクリッピング設定手段３０は、読み込んだ複数の断層画像のうち、レンダリングする対象を、関心領域であるＲＯＩとして設定する手段である。マスク設定手段４０は、読み込んだ複数の断層画像のうち、表示させないマスク領域を設定する手段である。不透明外接領域算出手段５０は、信号値に対して色成分の値と不透明度を対応付けて定義されたカラーマップを参照して、複数の断層画像の各画素に不透明度を対応付けて、その不透明度に基づいて不透明外接領域を算出する手段である。詳しくは後述するが、不透明外接領域とは、不透明度が０でない画素が存在する範囲に３次元の各軸方向において外接する領域である。
【0016】
レンダリング手段６０は、レンダリング画像の生成領域を設定し、複数の断層画像で構成されるボクセル画像を投影した画素値を、カラーマップを参照しながら決定して、レンダリング画像を生成する手段である。レンダリング手段６０は、さらにレイキャスティング手段６１、陰影付加手段６３を備えている。レイキャスティング手段６１は、レンダリング画像の各画素において、後述する有効領域と、視点に最も近い視線方向のベクトルとの交点を算出し、レンダリング画像の各画素において、算出される交点よりボクセル画像に対して視線方向に仮想光線を照射し、レイキャスティング法に基づいて所定の座標変換を行いながらボクセル画像の複数のボクセルを参照する手段である。陰影付加手段６３は、参照されたボクセル画像の信号値を基にカラーマップを参照して変換された各ボクセルの（Ｒ，Ｇ，Ｂ）で構成される色成分に対して、そのボクセルの近傍のボクセルの信号値を基に同様にカラーマップを参照して変換された不透明度αを基に、そのボクセルにおける勾配ベクトルを算出し、算出した勾配ベクトルを用いて、陰影値を算出し、色成分（Ｒ，Ｇ，Ｂ）の各値に陰影値を乗算した値に置き換える手段である。さらに、レイキャスティング手段６１は、起点座標探索手段６２を備えている。起点座標探索手段６２は、レイキャスティング手段６１が画素値の演算を行う際、視点座標系において投影の起点となる座標を探索する手段である。
【0017】
断層画像読込手段１０およびカラーマップ読込手段２０は、ＣＰＵ１が補助しながら、主にデータ入出力Ｉ／Ｆ５において実現される。ＲＯＩクリッピング設定手段３０およびマスク設定手段４０は、ＣＰＵ１が補助しながら、主に指示入力Ｉ／Ｆ４において実現される。不透明外接領域算出手段５０、レイキャスティング手段６１、陰影付加手段６３、レンダリング手段６０、起点座標探索手段６２は、ＣＰＵ１が、記憶装置３に記憶されているプログラムを実行することにより実現される。レンダリング画像出力手段７０は、ＣＰＵ１が補助しながら、主にフレームメモリ８と表示部６において実現される。
【0018】
図３に示した各構成手段は、現実には図２に示したように、コンピュータおよびその周辺機器等のハードウェアに専用のプログラムを搭載することにより実現される。すなわち、コンピュータが、専用のプログラムに従って各手段の内容を実行することになる。なお、本明細書において、コンピュータとは、ＣＰＵ、ＧＰＵ等の演算処理部を有し、データ処理が可能な装置を意味し、パーソナルコンピュータなどの汎用コンピュータだけでなく、ＧＰＵを搭載するタブレットなどの携帯端末や様々な装置に組み込まれたコンピュータも含む。
【0019】
＜２．処理動作＞
次に、図２、図３に示したボリュームレンダリング装置の処理動作について説明する。図４は、本実施形態に係るボリュームレンダリング装置の処理動作を示すフローチャートである。まず、断層画像読込手段１０が、複数の断層画像を読み込む。複数の断層画像は、図１（ａ）に示したような態様のものである。本実施形態では、ＤＩＣＯＭ形式の断層画像群Ｄｏ（ｘ，ｙ，ｚ）（－３２７６８≦Ｄｏ（ｘ，ｙ，ｚ）≦３２７６７，０≦ｘ≦Ｓｘ－１，０≦ｙ≦Ｓｙ－１，０≦ｚ≦Ｓｚ－１; 解像度：Ｒｘｙ，Ｒｚ）を読み込むものとする。断層画像群Ｄｏ（ｘ，ｙ，ｚ）は１画素が１６ビットで記録されるため、－３２７６８～３２７６７の信号値Ｄｏ（ｘ，ｙ，ｚ）が記録されている。断層画像群Ｄｏ（ｘ，ｙ，ｚ）は、ｘ座標のサイズ（画素数）がＳｘ、ｙ座標のサイズ（画素数）がＳｙの断層画像が、Ｓｚ枚積層されたものである。なお、以後の説明では、断層画像の枚数ＳｚをＭとも表現する。図１（ｂ）における１枚の断層画像の左右方向がｘ軸、上下方向がｙ軸に対応し、図１（ａ）に示した複数の断層画像を当該断層画像のＤＩＣＯＭヘッダに記録されている所定のスライス間隔（隣接する断層画像のImage Positionタグ情報の差分）で重ねた場合の積層方向がｚ軸に対応する。
次に、図２、図３に示したボリュームレンダリング装置の処理動作について説明する。図４は、本実施形態に係るボリュームレンダリング装置の処理動作を示すフローチャートである。まず、断層画像読込手段１０が、複数の断層画像を読み込む。複数の断層画像は、図１（ａ）に示したような態様のものである。本実施形態では、ＤＩＣＯＭ形式の断層画像群Ｄｏ（ｘ，ｙ，ｚ）（－３２７６８≦Ｄｏ（ｘ，ｙ，ｚ）≦３２７６７，０≦ｘ≦Ｓｘ－１，０≦ｙ≦Ｓｙ－１，０≦ｚ≦Ｓｚ－１; 解像度：Ｒｘｙ，Ｒｚ）を読み込むものとする。断層画像群Ｄｏ（ｘ，ｙ，ｚ）は１画素が１６ビットで記録されるため、－３２７６８～３２７６７の信号値Ｄｏ（ｘ，ｙ，ｚ）が記録されている。断層画像群Ｄｏ（ｘ，ｙ，ｚ）は、ｘ座標のサイズ（画素数）がＳｘ、ｙ座標のサイズ（画素数）がＳｙの断層画像が、Ｓｚ枚積層されたものである。なお、以後の説明では、断層画像の枚数ＳｚをＭとも表現する。図１（ｂ）における１枚の断層画像の左右方向がｘ軸、上下方向がｙ軸に対応し、図１（ａ）に示した複数の断層画像を当該断層画像のＤＩＣＯＭヘッダに記録されている所定のスライス間隔（隣接する断層画像のImage Positionタグ情報の差分）で重ねた場合の積層方向がｚ軸に対応する。
【0020】
ＤＩＣＯＭ形式の複数の断層画像を読み込んだら、ＲＯＩクリッピング設定手段３０が、ＲＯＩクリッピング設定を行う（ステップＳ１０）。ＲＯＩとはＲｅｇｉｏｎ ｏｆ Ｉｎｔｅｒｅｓｔの略であり、本明細書では３次元空間における関心領域を意味する（学術的には、ＲＯＩは２次元領域を指し、３次元領域の場合、ＶＯＩ（Ｖｏｌｕｍｅ ｏｆ Ｉｎｔｅｒｅｓｔ）と称するが、医用画像業界では３次元領域でもＲＯＩが多用される）。ここでは、レンダリング画像の作成対象とする範囲を示す。ＲＯＩをクリッピング（切り取り）して設定することにより、３次元的に任意の位置で被写体を断裁したボリュームレンダリング像を生成することができ、体表や骨に隠れた臓器や臓器の内部を描出するのに用いられる。ステップＳ４０の処理において、所定の範囲にボクセル画像の参照範囲が限定されるように、ステップＳ１０においてＲＯＩを設定する。この設定は、ｘｙｚ各軸方向のそれぞれの座標範囲を数値で指定することにより行われる。数値の指定は、直接行ってもよいし、マウス等により画面上で位置を指示し、指示入力Ｉ／Ｆ４を介して数値に変換することにより行ってもよい。
【0021】
読み込んだ各断層画像のｘ軸方向の最小の座標をＸｓｏ、ｘ軸方向の最大の座標をＸｅｏ、ｙ軸方向の最小の座標をＹｓｏ、ｙ軸方向の最大の座標をＹｅｏ、第１スライス（先頭：第１番目）の断層画像に対応する座標をＺｓｏ、最終スライス（末尾：第Ｍ番目）の断層画像に対応する座標をＺｅｏとすると、Ｘｓｏ≦Ｘｓ＜Ｘｅ≦Ｘｅｏ、Ｙｓｏ≦Ｙｓ＜Ｙｅ≦Ｙｅｏ、Ｚｓｏ≦Ｚｓ＜Ｚｅ≦Ｚｅｏを満たす［Ｘｓ，Ｘｅ］、［Ｙｓ，Ｙｅ］、［Ｚｓ，Ｚｅ］の直方体で定義される範囲をＲＯＩと定義し処理対象として特定する。すなわち、レンダリング画像を生成する対象となるボクセル画像領域を指定する６つの座標Ｘｓ，Ｘｅ，Ｙｓ，Ｙｅ，Ｚｓ，Ｚｅが定義される。この６つの座標により、ｘｙｚ各軸に沿った辺をもつ直方体が定義され、その直方体の内部がクリッピングされたＲＯＩとなる。本実施形態では、後述のステップＳ３０において不透明外接領域を算出するが、この不透明外接領域もｘｙｚ各軸の所定の範囲で定義される直方体形状の領域としている。
【0022】
ステップＳ１０におけるＲＯＩクリッピング設定の結果、処理対象の画素数が（Ｘｅｏ－Ｘｓｏ＋１）×（Ｙｅｏ－Ｙｓｏ＋１）×（Ｚｅｏ－Ｚｓｏ＋１）画素から、（Ｘｅ－Ｘｓ＋１）×（Ｙｅ－Ｙｓ＋１）×（Ｚｅ－Ｚｓ＋１）画素に削減されることになる。ＲＯＩクリッピング設定手段３０によるステップＳ１０のＲＯＩクリッピング設定を行うことにより、観察対象の臓器を露出させたボリュームレンダリング像が得られるという効果に加え、処理対象の画素が減って処理負荷を抑え応答性を向上させることができるという二重のメリットがあるため、設定される場合が多いが、必須の処理ではなく省略することも可能である。尚、ステップＳ１０のＲＯＩクリッピング設定の段階では、ＲＯＩクリッピングを適用する範囲を定義するだけで、具体的なＲＯＩのクリッピング処理は、ステップＳ４０の処理で実行される。各々の処理において、ボクセル画像のＲＯＩ内に限定して演算が行われることにより処理負荷が軽減され、ステップＳ４０のレンダリング処理でクリッピング（断裁）されたボリュームレンダリング像が生成される。
【0023】
ＲＯＩクリッピング設定手段３０により、ＲＯＩクリッピング設定が行われたら、次に、マスク設定手段４０が、マスク設定を行う（ステップＳ２０）。マスク設定とは、断層画像の各画素を表示させるか否かの設定である。マスク設定の手法は、特に限定されないが、本実施形態では、ステップＳ１０において読み込んだ断層画像群をＭＰＲ画像（アキシャル・コロナル・サジタルの３画面）として画面表示するとともに、本明細書に記載の方法または公知の手法により作成されるボリュームレンダリング像を表示させ、これら２次元投影された複数の画像の中から最も指示がしやすい画像を選択し、選択した画像の２次元領域を指示しながら奥行方向に延長させた３次元領域を、マスクして表示させない範囲とマスクせずに表示させる領域を特定する。このマスク設定により、断層画像群に対応した３次元空間におけるマスクデータが得られる。マスクデータＭａｓｋ（ｘ，ｙ，ｚ）は、断層画像群Ｄｏ（ｘ，ｙ，ｚ）の各画素に対応して“０”か“１”かのいずれかの値が設定されたものである。このような３次元マスクを作成する処理は必須の処理に近いが、省略することも可能である。観察対象の臓器を描出する方法として、上述のＲＯＩクリッピングもあるが、３次元マスクに比べ操作が簡便であり、処理負荷も減少するため、ＲＯＩクリッピング処理も併用される頻度が多い。ＲＯＩクリッピング設定手段３０、マスク設定手段４０により表示させるべき特定の領域が設定されることになるが、両手段で設定されている場合は、いずれも表示すべき領域として設定された領域が有効な領域となり、いずれか一方が設定されている場合は、設定された領域が有効な領域となる。ＲＯＩクリッピング設定手段３０、マスク設定手段４０のどちらも用いられない場合は、読み込まれた断層画像群の全領域が有効な領域となる。
【0024】
次に、不透明外接領域算出手段５０が、不透明外接領域の算出処理を行う（ステップＳ３０）。不透明外接領域とは、不透明度が０でない画素が存在する範囲に３次元の各軸方向において外接する領域である。外接する領域としているのは、不透明度が０でない画素が存在する範囲に外接する所定の辺（線分）で囲まれる領域とすることにより、レンダリング処理時に、回転演算を行った場合に、視線と交差する範囲を演算で高速に算出できるようにするためである。そのため、不透明外接領域は不透明度が０でない画素が存在する範囲に外接する３次元座標軸に鉛直な直方体領域に限定される。（斜め直方体や曲面形状では外接面の算出に処理負荷がかかり、かつ視線との交点を線形な数式で解析的に算出できず所望の高速化効果を達成できない。）従って、断層画像の２つの座標軸の方向と、断層画像を積層した方向に沿った辺で規定される直方体の内部の領域とすることが好ましい。具体的な処理としては、カラーマップ読込手段２０から読み込んだカラーマップを用いて、クリッピング設定後の断層画像のＲＯＩに含まれる各画素の信号値に対して対応する不透明度を付与し、付与した不透明度に基づいて不透明外接領域を算出する。
【0025】
図５は、本実施形態で用いるカラーマップを示す図である。図５において、信号値は、断層画像の各画素に記録された信号値であり、Ｒ、Ｇ、Ｂは、それぞれ赤、緑、青の色成分であり、αはそのボクセルの不透明度である。図５は、断層画像の信号値が１６ビットで記録され、－３２７６８～３２７６７の値をとる場合のカラーマップである。なお、ＣＴ画像の場合は、水（water）成分を０として－２０４８から２０４８の範囲に正規化される場合が多い。図５に示すように、カラーマップには、断層画像の各画素の信号値に対して、Ｒ、Ｇ、Ｂの各色成分と不透明度αが対応付けられている。このカラーマップは、レンダリング手段６０によるレンダリング処理においても用いられる（後述：ステップＳ４０）。ステップＳ３０における不透明外接領域の算出処理においては、カラーマップのうち不透明度αのみを用いる。カラーマップを参照して得られた不透明度αを用いることにより、断層画像においてα＝０となる不透明な画素が存在する領域を不透明外接領域として算出することができる。この不透明外接領域に対応する領域を有効領域として定義する。「不透明外接領域に対応する領域」として、不透明外接領域をそのまま処理対象として有効な有効領域とすることもできるが、本実施形態では、後述するように、不透明外接領域を１画素分広げた領域を、不透明外接領域に対応する有効領域としている。
【0026】
図５に示したようなカラーマップを用いることにより、断層画像において所定の信号値を有する箇所を、任意の色成分、任意の不透明度で表現することができる。断層画像においては、臓器、器官等の部位に応じて信号値が異なるため、ボリュームレンダリング像において、カラーマップを定義することにより、信号値の相違に基づいて、人為的に、特定の臓器を色分けしたり、特定の器官を透明にして背後に隠れた臓器を露出させたりすることができる。特定の部位の不透明度αを最大値（α＝２５５）に設定すれば、その部位だけが明瞭に表示され、その部位の奥に位置する部位は描出されない。逆に特定の部位の不透明度αを最小値（α＝０）に設定すれば、その部位は不可視（透明）になり、その部位の奥に位置する部位が露出される。特定の部位の不透明度αを中間調（０＜α＜２５５）に設定すれば、その部位が半透明で表示され、その奥に位置する部位が透かし合成されて表示される。また、色成分Ｒ、Ｇ、Ｂとしては、目的とする部位を視認し易くするための任意の値を設定することができる。
【0027】
したがって、カラーマップは、観察対象の部位ごとに設定しておくことができ、同じ断層画像群に対しても、異なるカラーマップを使用することにより、異なるボリュームレンダリング像が得られることになる。すなわち、肺用のカラーマップを用いれば、肺の部分が目立ったボリュームレンダリング像が得られ、心臓用のカラーマップを用いれば、心臓の部分が目立ったボリュームレンダリング像が得られる。ただし、各臓器の信号値分布は互いにオーバーラップしていることもあるため、特定の臓器のみを描出するカラーマップを作成することはできないことが多い。例えば、心臓用のカラーマップでは肋骨・胸骨も同時に描出され、心臓はこれらの骨に隠れてしまうため、ＲＯＩクリッピングを設定して肋骨を仮想的に断裁するなどの工夫が必要になる。
【0028】
本実施形態では、指示入力Ｉ／Ｆ４を介した操作により、使用中のカラーマップに変更を加えたり、異なるカラーマップを読み込んだりすることが可能となっている。そして、ボリュームレンダリング装置は、変更されたカラーマップまたは新たに読み込まれたカラーマップを参照して、ボリュームレンダリング像として提示するためのレンダリング画像を生成する。本実施形態のボリュームレンダリング装置は、断層画像からカラーマップを参照してＲＧＢα形式のボクセル画像を生成せずに、レンダリング画像を生成する過程で複数の断層画像で構成されるボクセル画像から計算に必要とするボクセルに限定してカラーマップを参照しながら色成分および不透明度に変換する方法をとっているため、カラーマップの変更に対してＲＧＢα形式のボクセル構造体を再構築する必要が無く、リアルタイムにボリュームレンダリング像を更新できる。
【0029】
本実施形態では、ステップＳ３０における不透明外接領域の算出処理は、複数のＣＰＵコアにより並行して行われる。ここで、ステップＳ３０における不透明外接領域の算出処理の詳細について説明する。図６は、ステップＳ３０における不透明外接領域の算出処理を並行して行う場合を示すフローチャートである。不透明外接領域算出手段５０は、ステップＳ１０で読み込んだ断層画像群をＮＡ個（ＮＡ≧２の整数）に分割し、各ＣＰＵコアで並列処理を行う。（ｘ，ｙ）座標を備えた各断層画像が積層された方向がｚ軸方向となるため、ｚ座標の値に応じてＮＡ個に分割する。具体的には、各断層画像群の先頭の断層画像の座標をＺ１、末尾の断層画像の座標をＺ２とし、先頭座標Ｚ１、末尾座標Ｚ２に、読み込んだ断層画像を割り当てていく。例えば、１番目のＣＰＵコアには、Ｚ１＝Ｚｓ、Ｚ２＝Ｚｓ＋（Ｚｅ－Ｚｓ＋１）／ＮＡ－１が設定される。各断層画像群には、断層画像が順に配置されているが、一端を先頭、他端を末尾としている。
【0030】
次に、処理対象とする断層画像群に対する不透明外接領域を算出する（ステップＳ３３０）。ステップＳ３３０における処理は、ＮＡ個のＣＰＵコアにより並行して処理が行われる。ステップＳ３３０における処理の詳細については後述する。
【0031】
ステップＳ３３０における断層画像群に対する不透明外接領域の算出を終えたら、ＮＡ個の断層画像群からなる断層画像全体に対して、不透明外接領域の算出を行う（ステップＳ３６０）。ステップＳ３６０における処理の詳細についても後述する。
【0032】
本実施形態では、複数に分割した断層画像群を複数のＣＰＵコアで並行に処理する。なお、本実施形態では、ｚ＝Ｚ１を先頭、ｚ＝Ｚ２を末尾として説明したが、ｚ＝Ｚ１を末尾、ｚ＝Ｚ２を先頭としてもよい。結局は、読み込んだ全ての断層画像の積層構成を変更せず、一方から他方に向かって処理することができるように、先頭と末尾を設定すればよい。
【0033】
上述のステップＳ３３０の処理の詳細について説明する。図７は、ステップＳ３３０における断層画像群に対する不透明外接領域の算出処理を示すフローチャートである。ステップＳ３３０の処理は、分割されたＮＡ個の各断層画像群に対して行われる。ステップＳ３３０においては、各断層画像群の先頭の断層画像の座標をＺ１、末尾の断層画像の座標をＺ２とし、先頭座標Ｚ１、末尾座標Ｚ２に、読み込んだ断層画像を割り当てる。例えば、１番目のＣＰＵコアには、Ｚ１＝Ｚｓ、Ｚ２＝Ｚｓ＋（Ｚｅ－Ｚｓ＋１）／ＮＡ－１が設定される。各断層画像群には、断層画像が順に配置されているが、一端を先頭、他端を末尾としている。
【0034】
まず、初期設定を行う（ステップＳ３３１）。具体的には、断層画像を特定するｚと、断層画像の座標値であるｘ，ｙと、ｘの最小値、最大値をそれぞれ決定するための変数Ｘｍｉｎ（ｔｉｄ）、Ｘｍａｘ（ｔｉｄ）と、ｙの最小値、最大値をそれぞれ決定するための変数Ｙｍｉｎ（ｔｉｄ）、Ｙｍａｘ（ｔｉｄ）と、ｚの最小値、最大値をそれぞれ決定するための変数Ｚｍｉｎ（ｔｉｄ）、Ｚｍａｘ（ｔｉｄ）に初期値を与える。（ｔｉｄ）は、断層画像群を特定する識別情報であり、ＮＡ個の各断層画像群に対応する値が与えられる。このｔｉｄは、処理スレッドを特定する識別情報でもある。
【0035】
ｘ，ｙとしては、ステップＳ１０においてＲＯＩクリッピング設定が行われている場合、その最小値Ｘｓ、最小値Ｙｓをそれぞれ与える。また、Ｘｍｉｎ（ｔｉｄ）、Ｘｍａｘ（ｔｉｄ）、Ｙｍｉｎ（ｔｉｄ）、Ｙｍａｘ（ｔｉｄ）、Ｚｍｉｎ（ｔｉｄ）、Ｚｍａｘ（ｔｉｄ）としては、ステップＳ１０においてＲＯＩクリッピング設定が行われている場合、最小値を決定するための変数には、ＲＯＩクリッピング設定による最大値を与え、最大値を決定するための変数には、ＲＯＩクリッピング設定による最小値を与える。したがって、ＲＯＩクリッピング設定が行われている場合、以下の〔数式１〕に示すように初期値が設定される。
【0036】
〔数式１〕
Ｘｍｉｎ（ｔｉｄ）＝Ｘｅ
Ｘｍａｘ（ｔｉｄ）＝Ｘｓ
Ｙｍｉｎ（ｔｉｄ）＝Ｙｅ
Ｙｍａｘ（ｔｉｄ）＝Ｙｓ
Ｚｍｉｎ（ｔｉｄ）＝Ｚｅ
Ｚｍａｘ（ｔｉｄ）＝Ｚｓ
Ｘｍｉｎ（ｔｉｄ）＝Ｘｅ
Ｘｍａｘ（ｔｉｄ）＝Ｘｓ
Ｙｍｉｎ（ｔｉｄ）＝Ｙｅ
Ｙｍａｘ（ｔｉｄ）＝Ｙｓ
Ｚｍｉｎ（ｔｉｄ）＝Ｚｅ
Ｚｍａｘ（ｔｉｄ）＝Ｚｓ
【0037】
なお、ステップＳ１０においてＲＯＩクリッピング設定が省略されている場合には、［Ｘｓ，Ｘｅ］、［Ｙｓ，Ｙｅ］、［Ｚｓ，Ｚｅ］に代えて断層画像群の最小値、最大値である［Ｘｓｏ，Ｘｅｏ］、［Ｙｓｏ，Ｙｅｏ］、［Ｚｓｏ，Ｚｅｏ］が初期値として設定される。
【0038】
初期設定が行われたら、断層画像の対象画素の信号値により図５に示したカラーマップを参照し、不透明度αを取得する（ステップＳ３３２）。次に、取得した不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の少なくとも一方が“０”であるか否かを判定する（ステップＳ３３３）。なお、ステップＳ２０におけるマスク設定が省略されている場合には、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）が“０”であるか否かの判定は行わず、不透明度αのみが“０”であるか否かの判定を行う。
【0039】
ステップＳ３３３において、“Ｎｏ”すなわち不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でない場合には、不透明とすべき画素であると考えられるため、最小値を決定するための変数Ｘｍｉｎ（ｔｉｄ）、Ｙｍｉｎ（ｔｉｄ）、Ｚｍｉｎ（ｔｉｄ）、最大値を決定するための変数Ｘｍａｘ（ｔｉｄ）、Ｙｍａｘ（ｔｉｄ）、Ｚｍａｘ（ｔｉｄ）を変更する（ステップＳ３３４）。具体的には、以下の〔数式２〕に従った処理を実行することにより、変数の値を変更する。
【0040】
〔数式２〕
ｘ＜Ｘｍｉｎ（ｔｉｄ）の場合、Ｘｍｉｎ（ｔｉｄ）＝ｘ
ｘ＞Ｘｍａｘ（ｔｉｄ）の場合、Ｘｍａｘ（ｔｉｄ）＝ｘ
ｙ＜Ｙｍｉｎ（ｔｉｄ）の場合、Ｙｍｉｎ（ｔｉｄ）＝ｙ
ｙ＞Ｙｍａｘ（ｔｉｄ）の場合、Ｙｍａｘ（ｔｉｄ）＝ｙ
ｚ＜Ｚｍｉｎ（ｔｉｄ）の場合、Ｚｍｉｎ（ｔｉｄ）＝ｚ
ｚ＞Ｚｍａｘ（ｔｉｄ）の場合、Ｚｍａｘ（ｔｉｄ）＝ｚ
ｘ＜Ｘｍｉｎ（ｔｉｄ）の場合、Ｘｍｉｎ（ｔｉｄ）＝ｘ
ｘ＞Ｘｍａｘ（ｔｉｄ）の場合、Ｘｍａｘ（ｔｉｄ）＝ｘ
ｙ＜Ｙｍｉｎ（ｔｉｄ）の場合、Ｙｍｉｎ（ｔｉｄ）＝ｙ
ｙ＞Ｙｍａｘ（ｔｉｄ）の場合、Ｙｍａｘ（ｔｉｄ）＝ｙ
ｚ＜Ｚｍｉｎ（ｔｉｄ）の場合、Ｚｍｉｎ（ｔｉｄ）＝ｚ
ｚ＞Ｚｍａｘ（ｔｉｄ）の場合、Ｚｍａｘ（ｔｉｄ）＝ｚ
【0041】
〔数式２〕により、対象とする断層画像群（ｔｉｄ）における不透明外接領域を定義するｘｙｚ各軸方向の最小値または最大値の各値が、前述の不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でないという条件を満足する断層画像（ｚ）における画素（ｘ，ｙ）の３次元座標値（ｘ，ｙ，ｚ）により更新される。
【0042】
ステップＳ３３３において、すなわち不透明度α等が“０”であると判定された場合、“０”でないと判定されてステップＳ３３４により変数を変更された場合のいずれにおいても、対象とする画素に対する処理を終えたことになるので、次に、対象画素をｘ軸方向において変更する（ステップＳ３３５）。具体的には、ｘ座標の値をインクリメントする。ｘ座標の変更の結果、ｘ≦Ｘｅである場合は、ステップＳ３３３に戻って、次の画素に対する処理を行う。一方、ステップＳ３３５において、ｘ＞Ｘｅとなった場合は、対象画素をｙ軸方向において変更する（ステップＳ３３６）。具体的には、ｙ座標の値をインクリメントする。ｘ座標の変更の結果、ｙ≦Ｙｅである場合は、ステップＳ３３３に戻って、次の画素に対する処理を行う。このようにして、あるｚで特定される断層画像の全ての画素（ｘ，ｙ）に対して処理を終えたら、対象画素をｚ軸方向において変更する（ステップＳ３３７）。すなわち、対象とする断層画像を変更する。ステップＳ３３７において、Ｚ＞ｚ２となった場合は、１つの断層画像群における（Ｚ２－Ｚ１＋１）個の断層画像に対して処理を終えたことになるので、図７に示したステップＳ３３０の処理を終了する。
【0043】
１つの断層画像群における全ての断層画像、画素について処理を行うことにより、不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でない全ての断層画像、画素について、ステップＳ３３４における〔数式２〕に従った処理が実行されることになる。これにより、途中に不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）いずれかが“０”の画素があったとしても不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でない最も外側の断層画像、画素の位置を取得することができる。
【0044】
図７に示したステップＳ３３０の処理をｚ＝Ｚ１～ｚ＝Ｚ２まで行うことにより、分割された１つの断層画像群における不透明外接領域が算出される。そして、ステップＳ３６０において、ＮＡ個の全ての断層画像群を集めた断層画像群全体の不透明外接領域を算出する。ステップＳ３６０の処理の詳細について説明する。図８は、ステップＳ３６０における断層画像群全体の不透明外接領域の算出処理を示すフローチャートである。まず、初期設定を行う（ステップＳ３６１）。具体的には、断層画像群を特定する識別情報であるｔｉｄと、処理スレッド数を示す固定値ＮＡと、ｘの最小値、最大値をそれぞれ決定するための変数Ｘｍｉｎ、Ｘｍａｘと、ｙの最小値、最大値をそれぞれ決定するための変数Ｙｍｉｎ、Ｙｍａｘと、ｚの最小値、最大値をそれぞれ決定するための変数Ｚｍｉｎ、Ｚｍａｘに初期値を与える。
【0045】
次に、各分割画像群において算出された最小値Ｘｍｉｎ（ｔｉｄ）、Ｙｍｉｎ（ｔｉｄ）、Ｚｍｉｎ（ｔｉｄ）、最大値Ｘｍａｘ（ｔｉｄ）、Ｙｍａｘ（ｔｉｄ）、Ｚｍａｘ（ｔｉｄ）を用いて、最小値を決定するための変数Ｘｍｉｎ、Ｙｍｉｎ、Ｚｍｉｎ、最大値を決定するための変数Ｘｍａｘ、Ｙｍａｘ、Ｚｍａｘを変更する（ステップＳ３６２）。具体的には、以下の〔数式３〕に従った処理を実行することにより、変数の値を変更する。
【0046】
〔数式３〕
Ｘｍｉｎ（ｔｉｄ）＜Ｘｍｉｎの場合、Ｘｍｉｎ＝Ｘｍｉｎ（ｔｉｄ）
Ｘｍａｘ（ｔｉｄ）＞Ｘｍａｘの場合、Ｘｍａｘ＝Ｘｍａｘ（ｔｉｄ）
Ｙｍｉｎ（ｔｉｄ）＜Ｙｍｉｎの場合、Ｙｍｉｎ＝Ｙｍｉｎ（ｔｉｄ）
Ｙｍａｘ（ｔｉｄ）＞Ｙｍａｘの場合、Ｙｍａｘ＝Ｙｍａｘ（ｔｉｄ）
Ｚｍｉｎ（ｔｉｄ）＜Ｚｍｉｎの場合、Ｚｍｉｎ＝Ｚｍｉｎ（ｔｉｄ）
Ｚｍａｘ（ｔｉｄ）＞Ｚｍａｘの場合、Ｚｍａｘ＝Ｚｍａｘ（ｔｉｄ）
Ｘｍｉｎ（ｔｉｄ）＜Ｘｍｉｎの場合、Ｘｍｉｎ＝Ｘｍｉｎ（ｔｉｄ）
Ｘｍａｘ（ｔｉｄ）＞Ｘｍａｘの場合、Ｘｍａｘ＝Ｘｍａｘ（ｔｉｄ）
Ｙｍｉｎ（ｔｉｄ）＜Ｙｍｉｎの場合、Ｙｍｉｎ＝Ｙｍｉｎ（ｔｉｄ）
Ｙｍａｘ（ｔｉｄ）＞Ｙｍａｘの場合、Ｙｍａｘ＝Ｙｍａｘ（ｔｉｄ）
Ｚｍｉｎ（ｔｉｄ）＜Ｚｍｉｎの場合、Ｚｍｉｎ＝Ｚｍｉｎ（ｔｉｄ）
Ｚｍａｘ（ｔｉｄ）＞Ｚｍａｘの場合、Ｚｍａｘ＝Ｚｍａｘ（ｔｉｄ）
【0047】
〔数式３〕により、対象とする断層画像群（ｔｉｄ）における不透明外接領域のｘｙｚ各軸方向の最小値または最大値の各値が、全体の不透明外接領域を定義する範囲を超える場合、当該断層画像群（ｔｉｄ）における不透明外接領域のｘｙｚ各軸方向の最小値または最大値の値により、全体の不透明外接領域を定義する範囲が拡大更新される。
【0048】
次に、対象とする断層画像群を変更する（ステップＳ３６３）。具体的には、ｔｉｄの値をインクリメントする。ｔｉｄの値をインクリメントした結果、ｔｉｄがＮＡより小さい場合は、全ての断層画像群に対して処理を終えていないことを意味するので、ステップＳ３６２に戻って、次の断層画像群に対して最小値、最大値を決定するための変数を変更する処理を行う。ｔｉｄの値をインクリメントした結果、ｔｉｄがＮＡ以上である場合は、全ての断層画像群に対して処理を終えたことを意味するので、不透明外接領域の境界処理を行う（ステップＳ３６４）。ステップＳ３６４では、モアレ対策として、不透明度αが０の透明な画素で構成される境界面（x=Xisまたはx=Xieまたはy=Yisまたはy=Yieまたはz=Zisまたはz=Zie）を不透明外接領域の最外周となる外接直方体のエッジに１画素幅だけ付加する。具体的には、ステップＳ３６２を繰り返すことにより算出された領域をｘｙｚ各軸方向に１画素ずつ拡張する処理を行う。より具体的には、以下の〔数式４〕に従った処理を実行することにより、変数の値を変更する。
【0049】
〔数式４〕
Ｘｉｓ＝Ｘｍｉｎ－１とし、Ｘｉｓ＜Ｘｓの場合、Ｘｉｓ＝Ｘｓ
Ｘｉｅ＝Ｘｍａｘ＋１とし、Ｘｉｅ＞Ｘｅの場合、Ｘｉｅ＝Ｘｅ
Ｙｉｓ＝Ｙｍｉｎ－１とし、Ｙｉｓ＜Ｙｓの場合、Ｙｉｓ＝Ｙｓ
Ｙｉｅ＝Ｙｍａｘ＋１とし、Ｙｉｅ＞Ｙｅの場合、Ｙｉｅ＝Ｙｅ
Ｚｉｓ＝Ｚｍｉｎ－１とし、Ｚｉｓ＜Ｚｓの場合、Ｚｉｓ＝Ｚｓ
Ｚｉｅ＝Ｚｍａｘ＋１とし、Ｚｉｅ＞Ｚｅの場合、Ｚｉｅ＝Ｚｅ
Ｘｉｓ＝Ｘｍｉｎ－１とし、Ｘｉｓ＜Ｘｓの場合、Ｘｉｓ＝Ｘｓ
Ｘｉｅ＝Ｘｍａｘ＋１とし、Ｘｉｅ＞Ｘｅの場合、Ｘｉｅ＝Ｘｅ
Ｙｉｓ＝Ｙｍｉｎ－１とし、Ｙｉｓ＜Ｙｓの場合、Ｙｉｓ＝Ｙｓ
Ｙｉｅ＝Ｙｍａｘ＋１とし、Ｙｉｅ＞Ｙｅの場合、Ｙｉｅ＝Ｙｅ
Ｚｉｓ＝Ｚｍｉｎ－１とし、Ｚｉｓ＜Ｚｓの場合、Ｚｉｓ＝Ｚｓ
Ｚｉｅ＝Ｚｍａｘ＋１とし、Ｚｉｅ＞Ｚｅの場合、Ｚｉｅ＝Ｚｅ
【0050】
ステップＳ３６４においては、境界となる画素を１画素分だけ外側に移動させている。このようにして外側に１画素分広げた領域を処理対象として有効な有効領域とする。有効領域の最外側の画素は、不透明度αが“０”であるはずなので、有効領域の最外周の画素、いわゆる有効領域の境界面に位置する画素は、全て透明な画素となる。したがって、不透明外接領域算出手段５０は、各軸方向における最外周となる画素に対応する不透明度が０となるように有効領域を算出することになる。ただし、有効領域の境界面に位置する画素の色成分（ＲＧＢ値）は、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）が０またはＲＯＩの範囲外であるか否かにかかわらず、断層画像の信号値を基にカラーマップで定義されるＲＧＢ値をそのまま付与する。これにより、有効領域の内部と外部で、不透明度は不連続になるがＲＧＢ値は連続性が維持され、この有効領域に対応するボクセル画像に対してレンダリング処理を行った場合、有効領域の境界面に発生するモアレを抑制することが可能となる。
【0051】
ステップＳ３６４における処理の結果得られる［Ｘｉｓ，Ｘｉｅ］、［Ｙｉｓ，Ｙｉｅ］、［Ｚｉｓ，Ｚｉｅ］の直方体で定義される範囲を有効領域と定義し処理対象として特定する。なお、図８のステップＳ３６４においては、不透明外接領域の外側に１画素分広げた領域を有効領域とすることにより、有効領域の各軸方向における最外周となるボクセル（画素）に対応する不透明度が０となるようにした。しかし、これに限定されず、不透明外接領域で規定される領域をそのまま有効領域とし、最外周となるボクセルの不透明度が０でない場合であっても、最外周となるボクセルの不透明度を強制的に０に置き換えることにより、最外周となるボクセルに対応する不透明度が０となるようにしてもよい。
【0052】
上記のようにして、ステップＳ３０の処理により不透明外接領域が算出される。ここで、ステップＳ１０のＲＯＩクリッピング処理により得られるＲＯＩと、ステップＳ３０の処理により得られる不透明外接領域を比較してみる。図９は、ＲＯＩと不透明外接領域を比較した図である。このうち、図９（ａ）は、ＲＯＩ、図９（ｂ）は、不透明外接領域を示している。図９（ａ）（ｂ）において、楕円状の実線は、複数の断層画像による空間内における不透明画素（α＞０）の分布範囲を示している。ここで、不透明画素とは、カラーマップを参照した結果、不透明度αが０より大きい値をもつ画素であり、不透明画素の分布範囲とは、不透明画素が存在している最大範囲を示している。すなわち、不透明画素の分布範囲内には、不透明度αが０である透明画素も存在している。ここでは、図９（ａ）（ｂ）に示す不透明画素の分布範囲は同一のものを示している。
【0053】
図９（ａ）においては、実線の直方体がＲＯＩを示し、破線の直方体が不透明外接領域を示している。図９（ｂ）においては、実線の直方体が不透明外接領域を示している。上述のように、ＲＯＩは、ステップＳ１０におけるＲＯＩクリッピング設定において、ｘｙｚ各軸方向のそれぞれの座標範囲を数値で指定することにより行われる。また、不透明外接領域は、ステップＳ３０における不透明外接領域の算出において、不透明画素の分布範囲に外接する不透明外接領域を算出することにより行われる。通常、ステップＳ１０におけるＲＯＩクリッピング設定においては、不透明画素の分布範囲を含むように利用者が設定し、ステップＳ３０における不透明外接領域の算出においては、不透明画素の分布範囲を含む最小の不透明外接領域が算出される。そのため、図９（ａ）に実線の直方体と破線の直方体で示したように、不透明外接領域がＲＯＩに完全に包含され、不透明外接領域の方がＲＯＩよりも小さくなる。
【0054】
しかし、ＲＯＩは、ステップＳ１０におけるＲＯＩクリッピング設定において任意に設定できるため、必ずしも不透明画素の分布範囲全体を含むように設定されるとは限らない。その場合、有効領域は、ＲＯＩと不透明外接領域の両方に含まれる範囲内で定義される。また、マスク設定処理（Ｓ２０）においてマスクが定義されている場合、図９のα＞０の範囲として示されている楕円体のような任意の３次元形状で囲まれる範囲の外側がマスクされα＝０となるが、不透明画素の分布範囲とは独立して定義される。即ち、図９のα＞０の範囲として示されている楕円体のような分布が２種類存在し、これらの２種類の分布が重なったα＞０の分布範囲に外接するように不透明外接領域が設定されることになり、不透明外接領域は更に小さくなる。
【0055】
ここで、不透明外接領域の内部における画素の分布について説明する。図９（ｃ）は、図９（ｂ）に対し、ドーナツ型のマスクデータＭａｓｋ（ｘ，ｙ，ｚ）を追加したもので、図９（ｃ）の三重の楕円で挟まれた白色の領域がＭａｓｋ（ｘ，ｙ，ｚ）＝０となるマスク領域である。図９（ｃ）においては、不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でない画素を網掛けで示しており、不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）のいずれかが“０”の画素を網掛けせずに示している。図９（ｃ）に示すように、不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でない画素の分布が内側と外側の二重になっていたとしても、不透明度α、マスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でない画素の分布の最も外側を含むように、不透明外接領域は設定される。内側のマスクデータＭａｓｋ（ｘ，ｙ，ｚ）の両方が“０”でない画素の分布に対しては、不透明外接領域は設定されない。
【0056】
このようにして、不透明外接領域に対応する範囲に縮小した有効領域を定義してボリュームレンダリングを行っても、得られるレンダリング画像は、従来の断層画像全体（ボクセル画像）を有効領域としてボリュームレンダリングを行う場合と差異は無い。しかし、後述するレンダリング画像の各画素に対応して探索される仮想光線上の有効領域における起点座標を迅速に求めることが可能になる。
【0057】
ステップＳ３０における不透明外接領域算出処理を終了したら、レンダリング手段６０が、レンダリング処理を行う（ステップＳ４０）。図１０は、ステップＳ４０におけるレンダリング処理を各ＣＰＵコアで並行して行う処理を示すフローチャートである。
【0058】
レンダリング画像の生成は、設定されたサイズにより定まる生成領域を複数に分割し、得られた各分割生成領域について並行して処理を行って分割画像を生成することにより行う。レンダリング画像のサイズは、ｘ軸方向ｙ軸方向ともに同一であり、画素数Ｓｉｚｅ×Ｓｉｚｅに設定される。まず、レンダリング手段６０は、並列処理を行うための分割する分割生成領域を初期設定する（ステップＳ４１０）。本実施形態では、ｙ座標については分割せず、ｘ座標についてＮＢ個に分割する。具体的には、各分割生成領域のｘ座標の先頭座標をＸ１、末尾座標をＸ２とし、先頭座標Ｘ１、末尾座標Ｘ２に、設定されたサイズＳｉｚｅ×Ｓｉｚｅで特定される全体のレンダリング画像の分割生成領域を割り当てていく。初期値として、Ｘ１＝０、Ｘ２＝Ｓｉｚｅ／ＮＢ－１が設定される。
【0059】
次に、生成する座標変換パラメータの設定を行う（ステップＳ４２０）。レンダリング手段６０は、以後のステップにおいて座標変換を実施する際に共通に用いるパラメータである座標変換パラメータを設定する。座標変換パラメータの設定は、ボクセルに対して仮想光線を照射する際、ボクセル単位に逐次座標変換処理を円滑に行うことができるようにするために、視点を基準とする視点座標系における３次元座標を基にボクセル画像が定義されているボクセル座標系の３次元座標に変換し、ボクセル座標系において対応するボクセルの信号値を取得するための処理である。
【0060】
ここで、座標変換の概念図と座標変換パラメータの具体例については説明する。図１１は、座標変換の概念図である。レンダリングにより得られるレンダリング画像は視点から見た画像であるので、レンダリング画像の各画素（ｘ，ｙ）の基準となる視点座標系は必ずしもボクセル座標系とは一致していない。そこでレンダリング手段６０は、ボクセル座標系に定義されているボクセル画像を視点座標系に変換する必要がある。すなわち、本来は、図１１（ａ）に示すように、投影面の各点からボクセル画像に対して任意の方向に仮想光線（レイ）を発射し、仮想光線が交差するボクセル画像内の指定範囲における信号値を複数抽出し、抽出された複数のボクセルの信号値を基にレンダリング画像の画素値を算出する。しかし、計算を簡単にするため、ボクセル画像に座標変換を施し、図１１（ｂ）に示すように、投影面の各点よりｚ軸負方向に指定範囲における画素値を算出する。ただし、図１１（ｂ）に示すような座標変換されたボクセル画像は実際には別途作成されるものではなく、各ボクセルごとに逆方向に座標変換を行いながらボクセル画像の画素値を都度取得するようにしている。
【0061】
座標変換を実施する際には、様々な座標変換パラメータを用いる。視点座標系とボクセル座標系の位置関係が同じであれば、座標変換パラメータも同じであると考えられる。即ち、レンダリング手段６０において、視点座標系で参照される全てのボクセルに対して、ボクセル座標系に座標変換を実施するための座標変換パラメータは同一である。そこでレンダリング手段６０は、ステップＳ４２０において座標変換パラメータをあらかじめ算出し、以後のステップにおいてはその座標変換パラメータを共通の座標変換パラメータとして流用する。座標変換パラメータとしては例えば以下のようなものが挙げられる。
【0062】
・所定の３次元座標系における回転を定義した３×３行列（逆方向に投影変換が行われるため、指示入力Ｉ／Ｆ４により指示され定義された３×３行列の逆行列）
・ｘｙｚ軸方向のオフセット値：Ｘｏｆｆ、Ｙｏｆｆ、Ｚｏｆｆ（単位：ボクセル）
・ＲＯＩクリッピングで設定された領域：ｘ軸方向：Ｘｓ～Ｘｅ，ｙ軸方向：Ｙｓ～Ｙｅ，ｚ軸方向：Ｚｓ～Ｚｅ
・不透明外接領域：ｘ軸方向：Ｘｉｓ～Ｘｉｅ，ｙ軸方向：Ｙｉｓ～Ｙｉｅ，ｚ軸方向：Ｚｉｓ～Ｚｉｅ
・拡大又は縮小倍率：Ｓｃａｌｅ（ｘｙｚ各軸について同じ値を用いる。）
・ｚ軸方向変倍率：Ｓｃｚ＝Ｒｘｙ／Ｒｚ（ｘｙ軸方向の解像度Ｒｘｙとｚ軸方向の解像度Ｒｚの間の比率）
・座標変換サブサンプル・微小オフセット値：ｘ軸方向ｄｘ，ｙ軸方向ｄｙ，ｚ軸方向ｄｚ
・仮想光線のサブサンプリング倍率：Ｓｒａｙ（通常は１、値が大きいと粗くなる、１未満だと高精細になる）
・ｘｙｚ軸方向のオフセット値：Ｘｏｆｆ、Ｙｏｆｆ、Ｚｏｆｆ（単位：ボクセル）
・ＲＯＩクリッピングで設定された領域：ｘ軸方向：Ｘｓ～Ｘｅ，ｙ軸方向：Ｙｓ～Ｙｅ，ｚ軸方向：Ｚｓ～Ｚｅ
・不透明外接領域：ｘ軸方向：Ｘｉｓ～Ｘｉｅ，ｙ軸方向：Ｙｉｓ～Ｙｉｅ，ｚ軸方向：Ｚｉｓ～Ｚｉｅ
・拡大又は縮小倍率：Ｓｃａｌｅ（ｘｙｚ各軸について同じ値を用いる。）
・ｚ軸方向変倍率：Ｓｃｚ＝Ｒｘｙ／Ｒｚ（ｘｙ軸方向の解像度Ｒｘｙとｚ軸方向の解像度Ｒｚの間の比率）
・座標変換サブサンプル・微小オフセット値：ｘ軸方向ｄｘ，ｙ軸方向ｄｙ，ｚ軸方向ｄｚ
・仮想光線のサブサンプリング倍率：Ｓｒａｙ（通常は１、値が大きいと粗くなる、１未満だと高精細になる）
【0063】
座標変換パラメータが設定されたら、座標変換サブサンプルのオフセットを行う（ステップＳ４３０）。具体的には、ｄｘ＝ｄｙ＝ｄｚ＝０に設定する。また、ｌ＝０に設定する。次に、分割生成領域について分割画像の生成処理を行う（ステップＳ４４０）。上述のように、レンダリング画像の生成領域は、ｘ座標においてのみＮＢ個に分割されているため、分割画像のｙ座標の画素数はＳｉｚｅであり、分割画像のｘ座標の画素数は（Ｘ２－Ｘ１＋１）となる。したがって、ステップＳ４４０においては、Ｓｉｚｅ×（Ｘ２－Ｘ１＋１）画素の分割画像を生成することになる。このステップＳ４４０の処理については後述する。
【0064】
続いて、座標変換サブサンプルの更新を行う（ステップＳ４５０）。具体的には、l←l+1, dx←dx+1/L, dy←dy+1/L, dz←dz+1/L を実行する。座標変換サブサンプルを更新しながら、ステップＳ４４０の処理はＬ回実行されることになる。
【0065】
ステップＳ４４０の処理は、複数のＣＰＵコアにより並行して行われる。このため、最初の分割範囲について分割画像の生成処理を１つのＣＰＵコアに割り当てて実行処理を起動したら、処理終了を待たずに次の分割生成領域を設定および生成処理の起動を進める（ステップＳ４６０）。具体的には、Ｘ１、Ｘ２の値をそれぞれＳｉｚｅ／ＮＢだけ増加し、次の分割生成領域におけるｘ座標の先頭座標Ｘ１、末尾座標Ｘ２を新たに設定する。
【0066】
ステップＳ４６０において新たな分割生成領域を設定した後、新たな分割生成領域がレンダリング画像全体のサイズを超えている場合、すなわち、ここでは、Ｘ２＞Ｓｉｚｅ－１を満たす場合は、全ての分割生成領域について設定および分割画像生成の処理の起動を終えたとして、図１０に示すステップＳ４０の処理を終了する。Ｘ２≦Ｓｉｚｅ－１を満たす場合は、全ての分割生成領域について処理を終えていないため、ステップＳ４２０に戻って、他のＣＰＵコアに割り当て、分割画像の生成を起動する処理を繰り替えす。
【0067】
図１０のステップＳ４４０の分割画像生成処理について説明する。図１２は、ステップＳ４４０の分割画像生成処理の詳細を示すフローチャートである。ここで、生成するレンダリング画像の分割画像を定義しておく。レンダリング画像は、各２次元座標(x,y) (X1≦x≦X2, 0≦y≦Size-1)で特定される画素においてＲＧＢ各色８ビット計２４ビットの画素値を有する画像である。そして、このレンダリング画像に対して、まず背景色を設定する（ステップＳ６１０）。具体的には、８ビットで０～２５５の値をとる背景色bg(c)、c=0(R),1(G),2(B)を定義する。そして、対象とするレンダリング画像の画素をｘ＝Ｘ１、ｙ＝０とする。
【0068】
次に、仮想光線強度と累積輝度値を初期化する（ステップＳ６２０）。具体的には、仮想光線強度Ｔｒａｎｓ＝１．０、累積輝度値Ｅｎｅｒｇｙ（ｃ）＝０．０に初期化する。累積輝度値Ｅｎｅｒｇｙ（ｃ）は、Ｒ，Ｇ，Ｂにそれぞれ対応するｃ＝０，１，２ごとに設定される。
【0069】
次に、有効領域と仮想光線との交点ｚ座標を算出する（ステップＳ６３０）。具体的には、レイキャスティング処理の対象とするレンダリング画像上の対象画素（ｘ，ｙ）に対して、視点座標系におけるｚ＝Ｚｏ（＝Ｓｉｚｅ／Ｓｒａｙ－１）からｚ＝０の範囲で、視点ｚ＝Ｚｏに最も近い有効領域（［Ｘｉｓ，Ｘｉｅ］、［Ｙｉｓ，Ｙｉｅ］、［Ｚｉｓ，Ｚｉｅ］）との交点Ｚｃを算出する。
【0070】
交点Ｚｃの算出にあたり、まず、仮想光線ベクトルを算出する。具体的には、視点座標系の（ｘ，ｙ，Ｚｏ）および（ｘ，ｙ，０）に対して上記座標変換パラメータを用いて各々座標変換を行い、実数値のボクセル座標を各々（ｘ１，ｙ１，ｚ１）、（ｘ２，ｙ２，ｚ２）とし、仮想光線ベクトルの各要素ｖｘ＝ｘ２－ｘ１、ｖｙ＝ｙ２－ｙ１、ｖｚ＝ｚ２－ｚ１を算出する。
【0071】
次に、仮想光線ベクトルと有効領域の６面との交点を算出する。具体的には、以下の〔数式６〕に従った処理を実行する。
【0072】
〔数式６〕
ｔｘ＝ｔｙ＝ｔｚ＝１０に初期化
｜ｖｘ｜≧１の場合、ｔｘ１＝（Ｘｉｓ－ｘ１）／ｖｘ，ｔｘ２＝（Ｘｉｅ－ｘ１）／ｖｘを算出し、ｔｘ１、ｔｘ２のうち小さい方をｔｘに設定する。
｜ｖｙ｜≧１の場合、ｔｙ１＝（Ｙｉｓ－ｙ１）／ｖｙ，ｔｙ２＝（Ｙｉｅ－ｙ１）／ｖｙを算出し、ｔｙ１、ｔｙ２のうち小さい方をｔｙに設定する。
｜ｖｚ｜≧１の場合、ｔｚ１＝（Ｚｉｓ－ｚ１）／ｖｚ，ｔｚ２＝（Ｚｉｅ－ｚ１）／ｖｚを算出し、ｔｚ１、ｔｚ２のうち小さい方をｔｚに設定する。
ｔｘ＝ｔｙ＝ｔｚ＝１０に初期化
｜ｖｘ｜≧１の場合、ｔｘ１＝（Ｘｉｓ－ｘ１）／ｖｘ，ｔｘ２＝（Ｘｉｅ－ｘ１）／ｖｘを算出し、ｔｘ１、ｔｘ２のうち小さい方をｔｘに設定する。
｜ｖｙ｜≧１の場合、ｔｙ１＝（Ｙｉｓ－ｙ１）／ｖｙ，ｔｙ２＝（Ｙｉｅ－ｙ１）／ｖｙを算出し、ｔｙ１、ｔｙ２のうち小さい方をｔｙに設定する。
｜ｖｚ｜≧１の場合、ｔｚ１＝（Ｚｉｓ－ｚ１）／ｖｚ，ｔｚ２＝（Ｚｉｅ－ｚ１）／ｖｚを算出し、ｔｚ１、ｔｚ２のうち小さい方をｔｚに設定する。
【0073】
〔数式６〕に従った処理により、ｔｘ、ｔｙ、ｔｚが得られる。続いて、視点（ｚ＝Ｚｏ）に最も近い交点Ｚｃを決定する。具体的には、以下の〔数式７〕に従った処理を実行する。
【0074】
〔数式７〕
Ｚｃ＝－１に初期化し、以下の１）～３）を実行
１）ｔｘ≦１の場合、Ｙ＝ｖｙ・ｔｘ＋ｙ１，Ｚ＝ｖｚ・ｔｘ＋ｚ１を算出、Ｙｉｓ≦Ｙ≦ＹｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅならば、Ｚはボクセル座標系における交点のｚ座標であるから視点座標系における交点のｚ座標をＺｖとしてＺｖ＝Ｚｏ・（１－ｔｘ）を算出、Ｚｖ＞ＺｃならばＺｃ＝Ｚｖに置換
２）ｔｙ≦１の場合、Ｘ＝ｖｘ・ｔｙ＋ｘ１，Ｚ＝ｖｚ・ｔｙ＋ｚ１を算出、Ｘｉｓ≦Ｘ≦ＸｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅならば、Ｚｖ＝Ｚｏ・（１－ｔｙ）を算出、Ｚｖ＞ＺｃならばＺｃ＝Ｚｖに置換
３）ｔｚ≦１の場合、Ｘ＝ｖｘ・ｔｚ＋ｘ１，Ｙ＝ｖｙ・ｔｚ＋ｙ１算出、Ｘｉｓ≦Ｘ≦ＸｉｅかつＹｉｓ≦Ｙ≦Ｙｉｅならば、Ｚｖ＝Ｚｏ・（１－ｔｚ）を算出、Ｚｖ＞ＺｃならばＺｃ＝Ｚｖに置換
Ｚｃ＝－１に初期化し、以下の１）～３）を実行
１）ｔｘ≦１の場合、Ｙ＝ｖｙ・ｔｘ＋ｙ１，Ｚ＝ｖｚ・ｔｘ＋ｚ１を算出、Ｙｉｓ≦Ｙ≦ＹｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅならば、Ｚはボクセル座標系における交点のｚ座標であるから視点座標系における交点のｚ座標をＺｖとしてＺｖ＝Ｚｏ・（１－ｔｘ）を算出、Ｚｖ＞ＺｃならばＺｃ＝Ｚｖに置換
２）ｔｙ≦１の場合、Ｘ＝ｖｘ・ｔｙ＋ｘ１，Ｚ＝ｖｚ・ｔｙ＋ｚ１を算出、Ｘｉｓ≦Ｘ≦ＸｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅならば、Ｚｖ＝Ｚｏ・（１－ｔｙ）を算出、Ｚｖ＞ＺｃならばＺｃ＝Ｚｖに置換
３）ｔｚ≦１の場合、Ｘ＝ｖｘ・ｔｚ＋ｘ１，Ｙ＝ｖｙ・ｔｚ＋ｙ１算出、Ｘｉｓ≦Ｘ≦ＸｉｅかつＹｉｓ≦Ｙ≦Ｙｉｅならば、Ｚｖ＝Ｚｏ・（１－ｔｚ）を算出、Ｚｖ＞ＺｃならばＺｃ＝Ｚｖに置換
【0075】
〔数式７〕において、１）～３）は、複数実行される場合もあれば、１つも実行されない場合もある。〔数式７〕に従った処理を実行することにより交点Ｚｃが得られる。ただし、得られた結果が視点（ｚ＝Ｚｏ）より手前（Ｚｃ＞Ｚｏ）となった場合、Ｚｃ＝Ｚｏと補正し、視点の位置が交点となる。ステップＳ６３０においては、レイキャスティング手段６１は、レンダリング画像の各画素において、有効領域と、視点に最も近い視線方向のベクトルとの交点Ｚｃを算出している。
【0076】
このようにしてステップＳ６３０の処理により、有効領域と仮想光線との交点ｚ座標Ｚｃが算出されたら、次にボクセルの参照処理を実行する（ステップＳ６４０）。ボクセルの参照処理においては、参照したボクセルの信号値に対してカラーマップを参照しながら不透明度に変換して不透明なボクセルを探索し探索された不透明なボクセルの信号値に対して更にカラーマップを参照しながらＲＧＢ値に変換して、仮想光線強度Ｔｒａｎｓ、ＲＧＢ各色（ｃ＝０，１，２）の累積輝度値Ｅｎｅｒｇｙ（ｃ）を算出する。ステップＳ６４０の処理については後述する。
【0077】
仮想光線強度Ｔｒａｎｓ、ＲＧＢ各色（ｃ＝０，１，２）の累積輝度値Ｅｎｅｒｇｙ（ｃ）が算出されたら、レンダリング画像の画素（ｘ，ｙ）における画素値を決定する（ステップＳ６５０）。具体的には、以下の〔数式８〕に従った処理を実行する。
【0078】
〔数式８〕
Ｉｍａｇｅ（ｘ，ｙ，ｃ）←Ｉｍａｇｅ（ｘ，ｙ，ｃ）＋Ｋ・｛Ｅｎｅｒｇｙ（ｃ）＋Ｔｒａｎｓ・ｂｇ（ｃ）｝／Ｌ
Ｉｍａｇｅ（ｘ，ｙ，ｃ）←Ｉｍａｇｅ（ｘ，ｙ，ｃ）＋Ｋ・｛Ｅｎｅｒｇｙ（ｃ）＋Ｔｒａｎｓ・ｂｇ（ｃ）｝／Ｌ
【0079】
〔数式８〕において、Ｋは強度倍率であり、本実施形態では、Ｋ＝１．０である。このようにして、画素（ｘ，ｙ）における画素値が決定したら、レンダリング画像の対象画素をｘ軸方向に１つだけ移動させた（ステップＳ６６０）後、ステップＳ６２０～ステップＳ６５０の処理を繰り返して実行し、次の画素（ｘ，ｙ）を対象として画素値を算出する。ｘ軸方向の上限Ｘ２に達したら、レンダリング画像の対象画素をｘ軸方向の初期値に移動させる（ステップＳ６７０）。この際、ｘ＝Ｘ１とする。そして、ステップＳ６２０～ステップＳ６５０の処理を繰り返して実行し、次の画素（ｘ，ｙ）を対象として画素値を算出する。このようにして、Ｘ１≦ｘ≦Ｘ２、０≦ｙ≦Ｓｉｚｅ－１の範囲の全ての画素について画素値を算出したら、図１２に示すステップＳ４４０のレイキャスティング処理を終了する。
【0080】
ここで、ステップＳ６４０のボクセルの参照処理について説明する。図１３は、ステップＳ６４０のボクセルの参照処理の詳細を示すフローチャートである。まず、探索開始座標を設定する（ステップＳ６４１）。具体的には、探索開始のｚ座標Ｚｓｔ＝Ｚｃに設定する。次に、探索開始座標Ｚｓｔから先頭の不透明ボクセルのｚ座標を起点座標として探索する（ステップＳ６４２）。
【0081】
ステップＳ６４２の起点座標探索処理については、図１４を用いて説明する。図１４は、起点座標探索手段６２が各レンダリング画像の画素において起点座標を探索する処理を示すフローチャートである。起点座標探索手段６２は、各対象画素（ｘ，ｙ）について図１４に示したフローチャートに従った処理を実行する。
【0082】
起点座標探索手段６２は、起点座標を探索する対象画素（ｘ，ｙ）、視点座標系の３次元座標（ｘ，ｙ，ｚ）をセットする（ステップＳ９０１）。探索開始座標であるｚ座標の初期値はｚ＝Ｚｓｔ（上限値）とする。ｚ座標の下限値は０、上限値はＺｓｔである。上限値であるＺｓｔは視点のｚ座標である。
【0083】
起点座標探索手段６２は、３次元座標（ｘ，ｙ，ｚ）について、座標変換を実施して対応するボクセル画像のボクセルが、透明か否かを判定する（ステップＳ９０２）。座標変換は、視点座標系をボクセル座標系に変換する操作であり、ＧＵＩ側の変換操作とは逆になる。ＧＵＩ側ではＲＯＩによるクリッピング、スケーリング、ｚ軸方向変倍処理、回転、オフセット（ｘｙｚ軸方向同時）の順に行うものと仮定し、与えられた視点座標系の３次元座標（ｘ，ｙ，ｚ）（整数値）に対応するＤＩＣＯＭ形式の複数の断層画像Ｄｏ（ｘ，ｙ，ｚ）で構成されるボクセル画像のボクセルに対応する実数の座標（ｘｒ，ｙｒ，ｚｒ）を算出する。
【0084】
具体的には、まず、整数値である視点座標系の３次元座標（ｘ，ｙ，ｚ）を実数値（ｘｘ，ｙｙ，ｚｚ）に変換する。視点座標系は、Ｓｉｚｅ×Ｓｉｚｅ×Ｓｉｚｅ／Ｓｒａｙの直方体であり、Ｓｘ×Ｓｙ×Ｓｚのボクセル画像とは独立して定義される。ＤＩＣＯＭ形式の断層画像のサイズＳｘ×Ｓｙは、通常Ｓｘ＝Ｓｙ＝５１２であるため、Ｓｉｚｅ＝５１２とすることが多い。
【0085】
視点座標系のｚ軸方向は仮想光線サブサンプル（１／Ｓｒａｙ）倍（Ｓｒａｙ≦１）に拡大されていることを考慮して、最初に、以下の〔数式９〕に従った処理を実行してｘｙｚ軸方向同時にオフセット処理を行う。
【0086】
〔数式９〕
ｘｘ＝ｘ－Ｓｉｚｅ／２＋ｄｘ＋Ｘｏｆｆ
ｙｙ＝ｙ－Ｓｉｚｅ／２＋ｄｙ＋Ｙｏｆｆ
ｚｚ＝ｚ・Ｓｒａｙ－Ｓｉｚｅ／２＋ｄｚ＋Ｚｏｆｆ
ｘｘ＝ｘ－Ｓｉｚｅ／２＋ｄｘ＋Ｘｏｆｆ
ｙｙ＝ｙ－Ｓｉｚｅ／２＋ｄｙ＋Ｙｏｆｆ
ｚｚ＝ｚ・Ｓｒａｙ－Ｓｉｚｅ／２＋ｄｚ＋Ｚｏｆｆ
【0087】
続いて、生成された３×３回転逆行列を用いて、以下の〔数式１０〕に従った処理を実行して回転処理を行う。
【0088】
〔数式１０〕
ｘｘ←Ｒ１１・ｘｘ＋Ｒ１２・ｙｙ＋Ｒ１３・ｚｚ
ｙｙ←Ｒ２１・ｘｘ＋Ｒ２２・ｙｙ＋Ｒ２３・ｚｚ
ｚｚ←Ｒ３１・ｘｘ＋Ｒ３２・ｙｙ＋Ｒ３３・ｚｚ
ｘｘ←Ｒ１１・ｘｘ＋Ｒ１２・ｙｙ＋Ｒ１３・ｚｚ
ｙｙ←Ｒ２１・ｘｘ＋Ｒ２２・ｙｙ＋Ｒ２３・ｚｚ
ｚｚ←Ｒ３１・ｘｘ＋Ｒ３２・ｙｙ＋Ｒ３３・ｚｚ
【0089】
回転処理後の座標として（ｘｘ，ｙｙ，ｚｚ）が得られる。続いて、以下の〔数式１１〕に従った処理を実行して、スケーリング、ｚ軸方向変倍処理を同時に行う。
【0090】
〔数式１１〕
ｘｒ＝ｘｘ／Ｓｃａｌｅ＋Ｓｘ／２
ｙｒ＝ｙｙ／Ｓｃａｌｅ＋Ｓｙ／２
ｚｒ＝ｚｚ／（Ｓｃａｌｅ・Ｓｃｚ）＋Ｓｚ／２
ｘｒ＝ｘｘ／Ｓｃａｌｅ＋Ｓｘ／２
ｙｒ＝ｙｙ／Ｓｃａｌｅ＋Ｓｙ／２
ｚｒ＝ｚｚ／（Ｓｃａｌｅ・Ｓｃｚ）＋Ｓｚ／２
【0091】
このようにしてボクセル座標系の座標（ｘｒ，ｙｒ，ｚｒ）が算出される。この座標（ｘｒ，ｙｒ，ｚｒ）は実数値である。次に、起点座標探索手段６２は、算出されたボクセル座標系の座標（ｘｒ，ｙｒ，ｚｒ）を整数値化する。具体的には、以下の〔数式１２〕に従った処理を実行して、整数値化を行う。
【0092】
〔数式１２〕
ｘｉ＝ＩＮＴ［ｘｒ＋０．５］
ｙｉ＝ＩＮＴ［ｙｒ＋０．５］
ｚｉ＝ＩＮＴ［ｚｒ＋０．５］
ｘｉ＝ＩＮＴ［ｘｒ＋０．５］
ｙｉ＝ＩＮＴ［ｙｒ＋０．５］
ｚｉ＝ＩＮＴ［ｚｒ＋０．５］
【0093】
〔数式１２〕において、ＩＮＴは［］内の数値の小数点以下を切り捨てて整数化することを意味する。したがって、〔数式１２〕に従った処理により、実数値の座標（ｘｒ，ｙｒ，ｚｒ）の各値は、四捨五入されて、整数値の座標（ｘｉ，ｙｉ，ｚｉ）が得られる。
【0094】
このようにして、視点座標系の座標（ｘ，ｙ，ｚ）に対応するボクセル座標系の座標（ｘｉ，ｙｉ，ｚｉ）が得られたら、有効領域、マスクデータＭａｓｋおよび不透明度補正テーブルＳαを用いて、以下の〔数式１３〕に従った処理を実行して、ボクセル座標系の座標（ｘｉ，ｙｉ，ｚｉ）に対応するボクセルが、透明か、不透明か、有効領域外であるかを求める。不透明度補正テーブルＳα（ｘ，ｙ，ｚ）は、あらかじめ３次元空間に不透明度に対する補正倍率が定義された多次元伝達関数である。同関数は、内臓が存在する中心部の領域の補正倍率を最大に設定し、周辺に向けて補正倍率を所定の関数で減衰させ、骨領域より外側（体表面：脂肪・筋肉・皮膚層、体外：寝台、ヘッドレスト・固定治具）において補正倍率が０になるように定義される。これを用いて、ボリュームレンダリングを実行すると、骨領域より外側が透明になり、内臓領域のみが露出したボリュームレンダリング像を得ることができる。マスクデータＭａｓｋおよび不透明度補正テーブルＳαが定義されていない場合、Ｍａｓｋ（ｘ，ｙ，ｚ）＝Ｓα（ｘ，ｙ，ｚ）＝１とする。
【0095】
〔数式１３〕
ｘｉ＜Ｘｉｓまたはｘｉ＞Ｘｉｅまたはｙｉ＜Ｙｉｓまたはｙｉ＞Ｙｉｅまたはｚｉ＜Ｚｉｓまたはｚｉ＞Ｚｉｅの場合（すなわち有効領域の外の場合）：Ｖα＝－１（無効値）
上記以外の場合（すなわち有効領域の場合）：Ｖα＝Ｖｏ（ｘｉ，ｙｉ，ｚｉ）
ただし、Ｖｏ（ｘ，ｙ，ｚ）＝Ｃｍａｐ（Ｄｏ（ｘ，ｙ，ｚ），３）・Ｍａｓｋ（ｘ，ｙ，ｚ）・Ｓα（ｘ，ｙ，ｚ）と定義する。Ｃｍａｐ（ｄ，ｃ）はカラーマップで、ｄは断層画像の信号値、ｃは色および不透明度の成分でｃ＝０（Ｒ），１（Ｇ），２（Ｂ），３（α）である。
ｘｉ＜Ｘｉｓまたはｘｉ＞Ｘｉｅまたはｙｉ＜Ｙｉｓまたはｙｉ＞Ｙｉｅまたはｚｉ＜Ｚｉｓまたはｚｉ＞Ｚｉｅの場合（すなわち有効領域の外の場合）：Ｖα＝－１（無効値）
上記以外の場合（すなわち有効領域の場合）：Ｖα＝Ｖｏ（ｘｉ，ｙｉ，ｚｉ）
ただし、Ｖｏ（ｘ，ｙ，ｚ）＝Ｃｍａｐ（Ｄｏ（ｘ，ｙ，ｚ），３）・Ｍａｓｋ（ｘ，ｙ，ｚ）・Ｓα（ｘ，ｙ，ｚ）と定義する。Ｃｍａｐ（ｄ，ｃ）はカラーマップで、ｄは断層画像の信号値、ｃは色および不透明度の成分でｃ＝０（Ｒ），１（Ｇ），２（Ｂ），３（α）である。
【0096】
〔数式１３〕において、有効領域の外である場合は、Ｖαに－１を設定しているが、信号値が有効領域の外であって無効であることを示すことができる値であれば、どのような値であってもよい。
【0097】
Ｖα＝０であればステップＳ９０３へ進み、Ｖα＜０であればステップＳ９０４へ進み、Ｖα＞０であればステップＳ９０５へ進む。すなわち、ボクセルが透明であればステップＳ９０３へ進み、ボクセルが有効領域外であればステップＳ９０４へ進み、ボクセルが不透明であればステップＳ９０５へ進む。
【0098】
ステップＳ９０３において、起点座標探索手段６２は、ｚ軸方向に所定数のボクセルをスキップする。具体的には現在のｚ座標から所定数ｍ２（例えばｍ２＝４）を減算する。なお、ステップＳ７４０、Ｓ７６０で用いられる所定数ｍ１と比較すると、ｍ２＞ｍ１である。本実施形態ではｍ１＝１である。減算した結果、ｚ座標が０以上である場合はステップＳ９０２へ戻って同様の処理を繰り返す。ｚ座標が０未満になった場合は、ステップＳ９０４へ進む。
【0099】
ステップＳ９０４において、起点座標探索手段６２は、現在の対象画素（ｘ，ｙ）に対応するｚ座標に対して、起点座標が見つからなかった旨を示す値（ｚの本来の下限値が０である場合、例えば－１などの負値）をセットして起点座標の探索処理を終了する。
【0100】
ステップＳ９０５に進んだ場合、起点座標探索手段６２は、変数ｉを初期化した後（ステップＳ９０５）、ｍ１だけ加算する（ステップＳ９０６）。ｉがｍ２まで到達した場合、または現在のｚ座標にｉを加算すると視点ｚ座標である上限値Ｚｓｔに達する場合は、ステップＳ９０８に進み、ｚ座標に対して、現在のｚ座標にｍ１を加算する前のｉを加算した値をセットする。ｉがｍ２まで到達しておらず、かつ現在のｚ座標にｉを加算しても視点ｚ座標である上限値Ｚｓｔまで達しない場合は、ステップＳ９０７へ進む。
【0101】
ステップＳ９０７において、起点座標探索手段６２は、３次元座標（ｘ，ｙ，ｚ＋ｉ）について、座標変換を実施して対応するボクセルが透明か否かを判定する。判定手順はステップＳ９０２と同じである。ただし、有効領域外であるか否かの判定は行わない。Ｖα＞０であればステップＳ９０６へ戻ってｉにｍ１を加算して同様の処理を繰り返す。Ｖα＝０であればステップＳ９０８へ進む。
【0102】
ステップＳ９０８において、起点座標探索手段６２は、現在の対象画素（ｘ，ｙ）に対応して、ステップＳ９０７の１つ手前のｚ座標（ｚ＋ｉ－ｍ１）をセットして起点座標の探索を終了する。
【0103】
結局、起点座標探索手段６２は、図１４に示した処理に従い、視点座標系において、レンダリング画像の各画素（ｘ，ｙ）におけるｚ座標を探索開始座標Ｚｓｔより下限値０まで所定数ｍ２ずつ変化させながら、視点座標系の３次元座標（ｘ，ｙ，ｚ）ごとに所定の座標変換を行って対応するボクセル画像のボクセルのボクセル値を順次取得し、最初に見つかった不透明なボクセルの視点座標系におけるｚ座標を、起点座標として探索している。より具体的には、ｚ座標を所定の上限値より所定の下限値の範囲で前記所定数ｍ１より大きい所定数ｍ２（ｍ２＞ｍ１）で減算させて、各３次元座標に対して、座標変換を行ってボクセルが不透明か否かを判定し（Ｓ９０２）、所定数ｍ２の減算後のｚ座標におけるボクセルの信号値が不透明である場合（Ｓ９０５）は、減算後のｚ座標より所定数ｍ１を増加させて（Ｓ９０６）各３次元座標（ｘ，ｙ，ｚ）に対して、座標変換を行ってボクセルが不透明か否かを判定している（Ｓ９０７）。
【0104】
図１４に示した処理を実行して、起点座標ｚが探索されたら、図１３に示すステップＳ６４２の処理を終えたことになる。得られた起点座標ｚがｚ＜０である場合は、不透明なボクセルが存在しなかったことを意味するので、図１３に示したステップＳ６４０の処理を終了する。得られた起点座標ｚがｚ≧０である場合は、３次元座標（ｘ，ｙ，ｚ）を座標変換してボクセルの不透明度を取得する（ステップＳ６４３）。
【0105】
具体的には、ステップＳ９０２、Ｓ９０７と同様に〔数式９〕～〔数式１１〕に従った処理を実行して、与えられた視点座標系の３次元座標（ｘ，ｙ，ｚ）（整数値）に対応するＤＩＣＯＭ形式の複数の断層画像Ｄｏ（ｘ，ｙ，ｚ）の各画素に１対１で対応するボクセル画像のボクセルに対応する実数の座標（ｘｒ，ｙｒ，ｚｒ）を算出する。
【0106】
次に、算出されたボクセル座標系の座標（ｘｒ，ｙｒ，ｚｒ）を整数値化する。ステップＳ９０２、Ｓ９０７の処理とは異なり、まず、小数点以下を切り捨て整数化した座標を（ｘｉａ，ｙｉａ，ｚｉａ）とし、切り捨てた小数点以下の端数を（ｗｘ，ｗｙ，ｗｚ）とする。すなわち、ｘｒ＝ｘｉａ＋ｗｘ，ｙｒ＝ｙｉａ＋ｗｙ，ｚｒ＝ｚｉａ＋ｗｚとなる。
【0107】
次に、有効領域およびマスクデータＭａｓｋを用いて、以下の〔数式１４〕に従った処理を実行して、ボクセル座標系の座標（ｘｉａ，ｙｉａ，ｚｉａ）に対応するボクセルの不透明度Ｖαを求める。以下の〔数式１４〕においては、３通りに分けて処理が行われる。〔数式１４〕において、Ｖｏ（ｘ，ｙ，ｚ）＝Ｃｍａｐ（Ｄｏ（ｘ，ｙ，ｚ），３）・Ｍａｓｋ（ｘ，ｙ，ｚ）・Ｓα（ｘ，ｙ，ｚ）と定義する。
【0108】
〔数式１４〕
１）、ｘｉａ＜Ｘｉｓまたはｘｉａ＞Ｘｉｅまたはｙｉａ＜Ｙｉｓまたはｙｉａ＞Ｙｉｅまたはｚｉａ＜Ｚｉｓまたはｚｉａ＞Ｚｉｅ（有効領域外）の場合：Ｖα＝－１
２）上記１）を満たさず、ｘｉａ＋１＞Ｘｉｅまたはｙｉａ＋１＞Ｙｉｅまたはｚｉａ＋１＞Ｚｉｅの場合、補間処理を実行せずにＶα＝Ｖｏ（ｘｉａ，ｙｉａ，ｚｉａ）
３）上記１）２）を満たさない場合、以下の補間処理を実行
Ｖα＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ，ｚｉａ）
＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ，ｚｉａ）
＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ＋１，ｚｉａ）
＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ）
＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ，ｚｉａ＋１）
＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ，ｚｉａ＋１）
＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ＋１，ｚｉａ＋１）
＋ｗｚ・ｗｙ・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ＋１）
１）、ｘｉａ＜Ｘｉｓまたはｘｉａ＞Ｘｉｅまたはｙｉａ＜Ｙｉｓまたはｙｉａ＞Ｙｉｅまたはｚｉａ＜Ｚｉｓまたはｚｉａ＞Ｚｉｅ（有効領域外）の場合：Ｖα＝－１
２）上記１）を満たさず、ｘｉａ＋１＞Ｘｉｅまたはｙｉａ＋１＞Ｙｉｅまたはｚｉａ＋１＞Ｚｉｅの場合、補間処理を実行せずにＶα＝Ｖｏ（ｘｉａ，ｙｉａ，ｚｉａ）
３）上記１）２）を満たさない場合、以下の補間処理を実行
Ｖα＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ，ｚｉａ）
＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ，ｚｉａ）
＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ＋１，ｚｉａ）
＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ）
＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ，ｚｉａ＋１）
＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ，ｚｉａ＋１）
＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｖｏ（ｘｉａ，ｙｉａ＋１，ｚｉａ＋１）
＋ｗｚ・ｗｙ・ｗｘ・Ｖｏ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ＋１）
【0109】
〔数式１４〕の３）においては、補間処理を行う際に、実数値のボクセル座標（ｘｒ，ｙｒ，ｚｒ）の近傍に位置する整数値のボクセル座標（ｘｉａ，ｙｉａ，ｚｉａ）等を８個取得し、これら８個の座標値に対応する断層画像、マスクデータ、不透明度補正テーブルの各値を取得し、カラーマップを参照しながら断層画像より取得した各信号値を不透明度に変換した値に所定の重みｗｘ、ｗｙ、ｗｚ、（１－ｗｘ）、（１－ｗｙ）、（１－ｗｚ）を掛けながら加算した値を、３次元座標に対応するボクセルの不透明度Ｖαとして与えている。〔数式１４〕において、有効領域の外である場合、すなわち有効な領域の外である場合は、Ｖαに－１を設定しているが、無効であることを示すことができる値であれば、どのような値であってもよい。さらに、ステップＳ６４３においては、不透明度Ｖαを８ビットで取り得る最大値である２５５で除算して、最大値が１以下となるＡｌｐｈａ（＝Ｖα／２５５）を算出する。
【0110】
ステップＳ６４３における処理の結果、Ａｌｐｈａ＝０かつｚ＞０であれば、有効領域内で透明な画素を意味するので、ステップＳ６４７へ進む。Ａｌｐｈａ＞０であれば、不透明な画素であることを意味するので、ステップＳ６４４へ進む。Ａｌｐｈａ＜０、またはＡｌｐｈａ＝０かつｚ＝０であれば、有効領域外であるか、探索範囲の下限に達したことを意味するので、図１３に示したステップＳ６４０の処理を終了する。
【0111】
ステップＳ６４７において、レイキャスティング手段６１は、視点の座標を所定数ｍ１だけ減じる。この所定数ｍ１は整数である。上述のように、本実施形態では、ｍ１＝１であり、ｍ１＜ｍ２である。そして、ステップＳ６４２へ戻って同様の処理を繰り返す。
【0112】
上述のように、Ａｌｐｈａ＞０であれば、不透明な画素であることを意味するので、レイキャスティング手段６１は、ＲＧＢ値Ｖ（ｃ）を取得する（ステップＳ６４４）。
【0113】
具体的には、ステップＳ６４３と同様、〔数式９〕～〔数式１１〕に従った処理を実行して、与えられた視点座標系の３次元座標（ｘ，ｙ，ｚ）（整数値）に対応するＤＩＣＯＭ形式の複数の断層画像Ｄｏ（ｘ，ｙ，ｚ）で構成されるボクセル画像のボクセルに対応する実数の座標（ｘｒ，ｙｒ，ｚｒ）を算出し、小数点以下を切り捨て整数化した座標を（ｘｉａ，ｙｉａ，ｚｉａ）とし、切り捨てた小数点以下の端数を（ｗｘ，ｗｙ，ｗｚ）とする。すなわち、ｘｒ＝ｘｉａ＋ｗｘ，ｙｒ＝ｙｉａ＋ｗｙ，ｚｒ＝ｚｉａ＋ｗｚとなる。
【0114】
次に、有効領域を用いて、以下の〔数式１５〕に従った処理を実行して、ボクセル座標系の座標（ｘｉａ，ｙｉａ，ｚｉａ）に対応するボクセルＲＧＢ値Ｖ（ｃ）を求める。以下の〔数式１５〕においては、３通りに分けて処理が行われる。〔数式１５〕において、Ｖｃ（ｘ，ｙ，ｚ，ｃ）＝Ｃｍａｐ（Ｄｏ（ｘ，ｙ，ｚ），ｃ）（ｃ＝０（Ｒ），１（Ｇ），２（Ｂ））と定義する。尚、ＲＧＢ値Ｖ（ｃ）を算出する際。〔数式１４〕と異なり、マスクデータＭａｓｋおよび不透明度補正テーブルＳαは参照しない。
【0115】
〔数式１５〕
１）ｘｉａ＜Ｘｉｓまたはｘｉａ＞Ｘｉｅまたはｙｉａ＜Ｙｉｓまたはｙｉａ＞Ｙｉｅまたはｚｉａ＜Ｚｉｓまたはｚｉａ＞Ｚｉｅ（有効領域外）の場合：Ｖ（ｃ）＝０（ｃ＝０，１，２）
２）上記１）を満たさず、ｘｉａ＋１＞Ｘｉｅまたはｙｉａ＋１＞Ｙｉｅまたはｚｉａ＋１＞Ｚｉｅの場合、補間処理を実行せずにＶ（ｃ）＝Ｖｃ（ｘｉａ，ｙｉａ，ｚｉａ，ｃ）
３）上記１）２）を満たさない場合、以下の補間処理を実行
Ｖ（ｃ）＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ，ｚｉａ，ｃ）＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ，ｚｉａ，ｃ）
＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ＋１，ｚｉａ，ｃ）
＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ，ｃ）
＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ，ｚｉａ＋１，ｃ）
＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ，ｚｉａ＋１，ｃ）
＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ＋１，ｚｉａ＋１，ｃ）
＋ｗｚ・ｗｙ・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ＋１，ｃ）
１）ｘｉａ＜Ｘｉｓまたはｘｉａ＞Ｘｉｅまたはｙｉａ＜Ｙｉｓまたはｙｉａ＞Ｙｉｅまたはｚｉａ＜Ｚｉｓまたはｚｉａ＞Ｚｉｅ（有効領域外）の場合：Ｖ（ｃ）＝０（ｃ＝０，１，２）
２）上記１）を満たさず、ｘｉａ＋１＞Ｘｉｅまたはｙｉａ＋１＞Ｙｉｅまたはｚｉａ＋１＞Ｚｉｅの場合、補間処理を実行せずにＶ（ｃ）＝Ｖｃ（ｘｉａ，ｙｉａ，ｚｉａ，ｃ）
３）上記１）２）を満たさない場合、以下の補間処理を実行
Ｖ（ｃ）＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ，ｚｉａ，ｃ）＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ，ｚｉａ，ｃ）
＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ＋１，ｚｉａ，ｃ）
＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ，ｃ）
＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ，ｚｉａ＋１，ｃ）
＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ，ｚｉａ＋１，ｃ）
＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｖｃ（ｘｉａ，ｙｉａ＋１，ｚｉａ＋１，ｃ）
＋ｗｚ・ｗｙ・ｗｘ・Ｖｃ（ｘｉａ＋１，ｙｉａ＋１，ｚｉａ＋１，ｃ）
【0116】
〔数式１５〕の３）においては、補間処理を行う際に、実数値のボクセル座標（ｘｒ，ｙｒ，ｚｒ）の近傍に位置する整数値のボクセル座標（ｘｉａ，ｙｉａ，ｚｉａ）等を８個取得し、これら８個の座標値に対応する断層画像の各信号値を取得し、カラーマップを参照しながら断層画像より取得した各信号値をＲＧＢ値に変換した値に所定の重みｗｘ、ｗｙ、ｗｚ、（１－ｗｘ）、（１－ｗｙ）、（１－ｗｚ）を掛けながら加算した値を、３次元座標に対応するボクセルのＲＧＢ値Ｖ（ｃ）として与えている。〔数式１５〕の１）において、有効領域の外である場合は、Ｖ（ｃ）に０を設定している。
【0117】
本実施形態では、ステップＳ６４４においてＲＧＢ値Ｖ（ｃ）を取得する際、さらに陰影付加処理を行っている。ステップＳ６４４においては、陰影付加手段６３が、陰影付加処理を行う。これは、所定の照明環境に対応した陰影をボクセル画像に付加する処理であり、具体的にはＸｓ≦ｘ≦ＸｅかつＹｓ≦ｙ≦ＹｅかつＺｓ≦ｚ≦Ｚｅの範囲のボクセルの不透明度αを参照しながら、Ｘｓ＜ｘ＜ＸｅかつＹｓ＜ｙ＜ＹｅかつＺｓ＜ｚ＜Ｚｅの範囲のボクセルの色成分を更新する処理である。このために、ある平行光源と環境光成分を設定し、付加すべき陰影の計算を行う。具体的には、視線の向きを示す単位ベクトルとして視線ベクトル（Ｅｘ，Ｅｙ，Ｅｚ）、平行光源の向きを示す単位ベクトルとして光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）、環境光成分Ａｂ、拡散反射係数Ｄｆ（０以上の実数）、鏡面反射係数Ｓｐ（０以上の実数）および鏡面光沢度Ｓｈ（１以上の整数）を設定する。
【0118】
視線ベクトル（Ｅｘ，Ｅｙ，Ｅｚ）は通常ｚ軸負方向で（Ｅｘ，Ｅｙ，Ｅｚ）＝（０，０，－１）の固定値で、光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）、拡散反射係数Ｄｆ、鏡面反射係数Ｓｐ、鏡面光沢度Ｓｈ、環境光成分Ａｂの数値は、表現したい像の状況に応じて適宜設定することができるが、例えば、（Ｌｘ，Ｌｙ，Ｌｚ）＝（０．５７７３５，０．５７７３５，－０．５７７３５）、Ｄｆ＝１．２、Ｓｐ＝１．７、Ｓｈ＝８、Ａｂ＝０．２と設定する。Ｌｘ，Ｌｙ，Ｌｚは絶対値が１未満の実数値（負値を含む）、Ｄｆは０以上の実数値、Ｓｐは０以上の実数値、Ｓｈは１以上の整数値、Ａｂは０以上１以下の実数値である。
【0119】
まず、変換パラメータとして生成された３×３回転逆行列を用いて、以下の〔数式１６〕に従った処理を実行して、光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）、視線ベクトル（Ｅｘ，Ｅｙ，Ｅｚ）の回転処理を行い、ボクセル座標系におけるベクトルに変換する。
【0120】
〔数式１６〕
Ｌｘ←Ｒ１１・Ｌｘ＋Ｒ１２・Ｌｙ＋Ｒ１３・Ｌｚ
Ｌｙ←Ｒ２１・Ｌｘ＋Ｒ２２・Ｌｙ＋Ｒ２３・Ｌｚ
Ｌｚ←Ｒ３１・Ｌｘ＋Ｒ３２・Ｌｙ＋Ｒ３３・Ｌｚ
Ｅｘ←Ｒ１１・Ｅｘ＋Ｒ１２・Ｅｙ＋Ｒ１３・Ｅｚ
Ｅｙ←Ｒ２１・Ｅｘ＋Ｒ２２・Ｅｙ＋Ｒ２３・Ｅｚ
Ｅｚ←Ｒ３１・Ｅｘ＋Ｒ３２・Ｅｙ＋Ｒ３３・Ｅｚ
Ｌｘ←Ｒ１１・Ｌｘ＋Ｒ１２・Ｌｙ＋Ｒ１３・Ｌｚ
Ｌｙ←Ｒ２１・Ｌｘ＋Ｒ２２・Ｌｙ＋Ｒ２３・Ｌｚ
Ｌｚ←Ｒ３１・Ｌｘ＋Ｒ３２・Ｌｙ＋Ｒ３３・Ｌｚ
Ｅｘ←Ｒ１１・Ｅｘ＋Ｒ１２・Ｅｙ＋Ｒ１３・Ｅｚ
Ｅｙ←Ｒ２１・Ｅｘ＋Ｒ２２・Ｅｙ＋Ｒ２３・Ｅｚ
Ｅｚ←Ｒ３１・Ｅｘ＋Ｒ３２・Ｅｙ＋Ｒ３３・Ｅｚ
【0121】
回転処理後の座標としてボクセル座標系における光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）、視線ベクトル（Ｅｘ，Ｅｙ，Ｅｚ）が得られる。次に、陰影付加の対象とするボクセル（ｘｉａ，ｙｉａ，ｚｉａ）に対して、ボクセル画像のｘｙ軸方向の解像度をＲｘｙ（画素間隔の逆数でｍｍあたりの画素数で、ｘ軸方向とｙ軸方向は同一）、ボクセル画像のｚ軸方向の解像度をＲｚ（断層画像のスライス間隔の逆数で、ｍｍあたりの断層画像のスライス数）として、Ｒｚ／Ｒｘｙ＜１の場合、以下の〔数式１７〕に従った処理を実行して、実数値ｚ座標ｚｒ１、ｚｒ２を算出する。これは、ｚ軸方向の解像度がｘｙ軸方向の解像度に比べて低い場合、算出される勾配ベクトルおよび陰影値にｚ軸方向に沿って段差が発生するため、勾配ベクトルに対してｚ軸方向に平滑化を施すためである。
【0122】
〔数式１７〕
ｚｒ１＝（ｚｉａ＋１）・Ｒｚ／Ｒｘｙ
ｚｒ２＝（ｚｉａ－１）・Ｒｚ／Ｒｘｙ
ｚｒ１＝（ｚｉａ＋１）・Ｒｚ／Ｒｘｙ
ｚｒ２＝（ｚｉａ－１）・Ｒｚ／Ｒｘｙ
【0123】
次に、算出された実数値ｚ座標ｚｒ１、ｚｒ２を整数値化する。まず、小数点以下を切り捨て整数化した実数値ｚ座標をｚｉ１、ｚｉ２とし、切り捨てた小数点以下の端数成分をｗｚ１、ｗｚ２とする。すなわち、ｚｒ１＝ｚｉ１＋ｗｚ１，ｚｒ２＝ｚｉ２＋ｗｚ２となる。
【0124】
続いて、ボクセル（ｘ，ｙ，ｚ）における不透明度をＶα（ｘ，ｙ，ｚ）＝Ｃｍａｐ（Ｄｏ（ｘ，ｙ，ｚ），３）として、ボクセル（ｘｉａ，ｙｉａ，ｚｉａ）に対して、ｘｙ軸方向に各々正負２近傍、ｚ軸正負方向に各々２近傍からなる全８近傍ボクセルの不透明度であるＶα（ｘｉａ＋１，ｙｉａ，ｚｉａ）、Ｖα（ｘｉａ－１，ｙｉａ，ｚｉａ）、Ｖα（ｘｉａ，ｙｉａ＋１，ｚｉａ）、Ｖα（ｘｉａ，ｙｉａ－１，ｚｉａ）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ１）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ１＋１）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ２）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ２＋１）を取得する。
【0125】
Ｒｚ／Ｒｘｙ≧１の場合、〔数式１７〕は実行されないため、Ｖα（ｘｉａ，ｙｉａ，ｚｉ１）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ１＋１）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ２）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ２＋１）の代わりに、ｚ軸方向に正負２近傍のボクセルの不透明度Ｖα（ｘｉａ，ｙｉａ，ｚｉａ＋１）とＶα（ｘｉａ，ｙｉａ，ｚｉａ－１）を取得する。
【0126】
ただし、
ｘｉａ＋１＞ＸｉｅのときＶα（ｘｉａ＋１，ｙｉａ，ｚｉａ）＝０、
ｘｉａ－１＜ＸｉｓのときＶα（ｘｉａ－１，ｙｉａ，ｚｉａ）＝０、
ｙｉａ＋１＞ＹｉｅのときＶα（ｘｉａ，ｙｉａ＋１，ｚｉａ）＝０、
ｙｉａ－１＜ＹｉｓのときＶα（ｘｉａ，ｙｉａ－１，ｚｉａ）＝０、
Ｒｚ／Ｒｘｙ＜１の場合、
ｚｉ１＞ＺｉｅのときＶα（ｘｉａ，ｙｉａ，ｚｉ１）＝０、
ｚｉ１＋１＜ＺｉｅのときＶα（ｘｉａ，ｙｉａ，ｚｉ１＋１）＝０、
ｚｉ２＜ＺｉｓのときＶα（ｘｉａ，ｙｉａ，ｚｉ２）＝０、
ｚｉ２＋１＞ＺｉｓのときＶα（ｘｉａ，ｙｉａ，ｚｉ２＋１）＝０、
Ｒｚ／Ｒｘｙ≧１の場合、
ｚｉａ＋１＞ＺｉｅのときＶα（ｘｉａ，ｙｉａ＋１，ｚｉａ＋１）＝０、
ｚｉａ－１＜ＺｉｓのときＶα（ｘｉａ，ｙｉａ－１，ｚｉａ－１）＝０、
とする。
ｘｉａ＋１＞ＸｉｅのときＶα（ｘｉａ＋１，ｙｉａ，ｚｉａ）＝０、
ｘｉａ－１＜ＸｉｓのときＶα（ｘｉａ－１，ｙｉａ，ｚｉａ）＝０、
ｙｉａ＋１＞ＹｉｅのときＶα（ｘｉａ，ｙｉａ＋１，ｚｉａ）＝０、
ｙｉａ－１＜ＹｉｓのときＶα（ｘｉａ，ｙｉａ－１，ｚｉａ）＝０、
Ｒｚ／Ｒｘｙ＜１の場合、
ｚｉ１＞ＺｉｅのときＶα（ｘｉａ，ｙｉａ，ｚｉ１）＝０、
ｚｉ１＋１＜ＺｉｅのときＶα（ｘｉａ，ｙｉａ，ｚｉ１＋１）＝０、
ｚｉ２＜ＺｉｓのときＶα（ｘｉａ，ｙｉａ，ｚｉ２）＝０、
ｚｉ２＋１＞ＺｉｓのときＶα（ｘｉａ，ｙｉａ，ｚｉ２＋１）＝０、
Ｒｚ／Ｒｘｙ≧１の場合、
ｚｉａ＋１＞ＺｉｅのときＶα（ｘｉａ，ｙｉａ＋１，ｚｉａ＋１）＝０、
ｚｉａ－１＜ＺｉｓのときＶα（ｘｉａ，ｙｉａ－１，ｚｉａ－１）＝０、
とする。
【0127】
続いて、Ｒｚ／Ｒｘｙ＜１の場合、Ｖα（ｘｉａ，ｙｉａ，ｚｉ１）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ１＋１）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ２）、Ｖα（ｘｉａ，ｙｉａ，ｚｉ２＋１）の４近傍ボクセルに対して、端数成分ｗｚ１、ｗｚ２を用いて、以下の〔数式１８〕に従った処理を実行して、ｚ方向に補間し、ｚ軸方向補間成分Ｖ´ａ、Ｖ´ｂを算出する。
【0128】
〔数式１８〕
Ｖ´ａ＝（１－ｗｚ１）・Ｖα（ｘｉａ，ｙｉａ，ｚｉ１）＋ｗｚ１・Ｖα（ｘｉａ，ｙｉａ，ｚｉ１＋１）
Ｖ´ｂ＝（１－ｗｚ２）・Ｖα（ｘｉａ，ｙｉａ，ｚｉ２）＋ｗｚ２・Ｖα（ｘｉａ，ｙｉａ，ｚｉ２＋１）
Ｖ´ａ＝（１－ｗｚ１）・Ｖα（ｘｉａ，ｙｉａ，ｚｉ１）＋ｗｚ１・Ｖα（ｘｉａ，ｙｉａ，ｚｉ１＋１）
Ｖ´ｂ＝（１－ｗｚ２）・Ｖα（ｘｉａ，ｙｉａ，ｚｉ２）＋ｗｚ２・Ｖα（ｘｉａ，ｙｉａ，ｚｉ２＋１）
【0129】
そして、ボクセル画像のボクセル（ｘｉａ，ｙｉａ，ｚｉａ）における勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）を以下の〔数式１９〕に従った処理を実行することにより算出する。
【0130】
〔数式１９〕
Ｇｘ＝Ｖα（ｘｉａ＋１，ｙｉａ，ｚｉａ）－Ｖα（ｘｉａ－１，ｙｉａ，ｚｉａ）
Ｇｙ＝Ｖα（ｘｉａ，ｙｉａ＋１，ｚｉａ）－Ｖα（ｘｉａ，ｙｉａ－１，ｚｉａ）
Ｒｚ／Ｒｘｙ＜１の場合、Ｇｚ＝Ｖ´ａ－Ｖ´ｂ
Ｒｚ／Ｒｘｙ≧１の場合、
Ｇｚ＝Ｖα（ｘｉａ，ｙｉａ，ｚｉａ＋１）－Ｖα（ｘｉａ，ｙｉａ，ｚｉａ－１）
Ｇ＝（Ｇｘ2＋Ｇｙ2＋Ｇｚ2）1/2
Ｇｘ＝Ｖα（ｘｉａ＋１，ｙｉａ，ｚｉａ）－Ｖα（ｘｉａ－１，ｙｉａ，ｚｉａ）
Ｇｙ＝Ｖα（ｘｉａ，ｙｉａ＋１，ｚｉａ）－Ｖα（ｘｉａ，ｙｉａ－１，ｚｉａ）
Ｒｚ／Ｒｘｙ＜１の場合、Ｇｚ＝Ｖ´ａ－Ｖ´ｂ
Ｒｚ／Ｒｘｙ≧１の場合、
Ｇｚ＝Ｖα（ｘｉａ，ｙｉａ，ｚｉａ＋１）－Ｖα（ｘｉａ，ｙｉａ，ｚｉａ－１）
Ｇ＝（Ｇｘ2＋Ｇｙ2＋Ｇｚ2）1/2
【0131】
〔数式１９〕において、勾配ベクトルのｚ軸方向成分Ｇｚは、当該ボクセルに対して、ｚ軸方向に＋Ｒｚ／Ｒｘｙずれた座標の近傍の２つのボクセルの不透明度に所定の重みをかけて加算した不透明度と、ｚ軸方向に－Ｒｚ／Ｒｘｙずれた座標の近傍の２つのボクセルの不透明度に所定の重みをかけて加算した不透明度との差分としている。〔数式１９〕において、Ｇは勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）の大きさである。〔数式１９〕に示すように、陰影付加手段６３は、ボクセル画像の当該ボクセルのｘ軸方向、ｙ軸方向、ｚ軸方向の各々２近傍（Ｒｚ／Ｒｘｙ＜１の場合、ｚ軸方向は４近傍）のボクセルの不透明度の差分値を勾配ベクトルのｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分として算出し、あらかじめ定義されたｚ軸方向変倍率に基づいてｘ軸方向・ｙ軸方向成分またはｚ軸方向成分に補正を施した単位ベクトルを、当該ボクセルにおける勾配ベクトルとして算出する。さらに、Ｇ≧１の場合、勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）をＧで除算し、勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）を単位ベクトルに置き換える。Ｇ＜１の場合、ボクセル（ｘｉａ，ｙｉａ，ｚｉａ）において不透明度の勾配が無いため陰影値を算出できず、（Ｇｘ，Ｇｙ，Ｇｚ）＝（０，０，０）とし陰影係数Ｓとして環境光成分Ａｂを与える。
【0132】
勾配ベクトルが算出されたら、公知のＰｈｏｎｇシェーディングモデルにより陰影係数Ｓを算出する。具体的には、まず、鏡面反射係数Ｓｐ＞０の場合、単位ベクトルである勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）、光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）を用いて、以下の〔数式２０〕に従った処理を実行することにより、正反射ベクトル（Ｆｘ，Ｆｙ，Ｆｚ）を算出する。鏡面反射係数Ｓｐ＝０の場合、（Ｆｘ，Ｆｙ，Ｆｚ）＝（０，０，０）とする。
【0133】
〔数式２０〕
Ｆｘ＝２・Ｇｘ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）－Ｌｘ
Ｆｙ＝２・Ｇｙ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）－Ｌｙ
Ｆｚ＝２・Ｇｚ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）－Ｌｚ
Ｆｘ＝２・Ｇｘ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）－Ｌｘ
Ｆｙ＝２・Ｇｙ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）－Ｌｙ
Ｆｚ＝２・Ｇｚ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）－Ｌｚ
【0134】
〔数式２０〕においては、正反射ベクトルは、勾配ベクトルとのなす角が光源ベクトルと勾配ベクトルとのなす角と同一で、かつ光源ベクトルとのなす角が光源ベクトルと勾配ベクトルとのなす角の２倍になるように定義している。次に、勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）、光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）、正反射ベクトル（Ｆｘ，Ｆｙ，Ｆｚ）、視線ベクトル（Ｅｘ，Ｅｙ，Ｅｚ）を全て単位ベクトルとして用い、反射係数Ｄｆ（０以上の実数）、鏡面反射係数Ｓｐ（０以上の実数）、鏡面光沢度Ｓｈ（１以上の整数）、環境光成分Ａｂの各係数を用いて、以下の〔数式２１〕に従った処理を実行することにより、ボクセル（ｘｉａ，ｙｉａ，ｚｉａ）における陰影係数Ｓ（ｘｉａ，ｙｉａ，ｚｉａ）を算出する。
【0135】
〔数式２１〕
Ｓ（ｘｉａ，ｙｉａ，ｚｉａ）＝Ｄｆ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）＋Ｓｐ・（Ｆｘ・Ｅｘ＋Ｆｙ・Ｅｙ＋Ｆｚ・Ｅｚ）Sh＋Ａｂ
Ｓ（ｘｉａ，ｙｉａ，ｚｉａ）＝Ｄｆ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）＋Ｓｐ・（Ｆｘ・Ｅｘ＋Ｆｙ・Ｅｙ＋Ｆｚ・Ｅｚ）Sh＋Ａｂ
【0136】
〔数式２１〕における“Ｄｆ・（Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ）”は拡散反射光であり、 “Ｓｐ・（Ｆｘ・Ｅｘ＋Ｆｙ・Ｅｙ＋Ｆｚ・Ｅｚ）Sh” は鏡面反射光である。〔数式２１〕においては、算出された勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）と光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）との内積値を基にボクセルの拡散反射光を算出し、鏡面反射係数Ｓｐに正の値が設定されている場合、更に勾配ベクトル（Ｇｘ，Ｇｙ，Ｇｚ）に対する所定の光源ベクトル（Ｌｘ，Ｌｙ，Ｌｚ）の正反射ベクトル（Ｆｘ，Ｆｙ，Ｆｚ）を算出し、正反射ベクトル（Ｆｘ，Ｆｙ，Ｆｚ）と所定の視線ベクトル（Ｅｘ，Ｅｙ，Ｅｚ）との内積値を基にボクセルの鏡面反射光を算出し、算出された拡散反射光と鏡面反射光を基に算出した陰影値をレイキャスティング手段により取得された当該ボクセルの色成分に乗算し、当該ボクセルの色成分を下記のように更新する。また、鏡面反射光として、正反射ベクトル（Ｆｘ，Ｆｙ，Ｆｚ）と視線ベクトル（Ｅｘ，Ｅｙ，Ｅｚ）との内積値に対して、鏡面光沢度Ｓｈで累乗した値に鏡面反射係数Ｓｐを乗算した値としている。陰影係数Ｓ（ｘｉａ，ｙｉａ，ｚｉａ）が算出されたら、陰影係数Ｓ（ｘｉａ，ｙｉａ，ｚｉａ）を用いて、〔数式１５〕により算出されたボクセル（ｘｉａ，ｙｉａ，ｚｉａ）におけるＲＧＢ値Ｖ（ｃ）＝０（ｃ＝０，１，２）を、以下の〔数式２２〕に従った処理を実行することにより、補正する。すなわち、陰影を付加する。
【0137】
〔数式２２〕
Ｖ（ｃ）←Ｓ（ｘｉａ，ｙｉａ，ｚｉａ）・Ｖ（ｃ）
Ｖ（ｃ）←Ｓ（ｘｉａ，ｙｉａ，ｚｉａ）・Ｖ（ｃ）
【0138】
〔数式２２〕において、ｃ＝０，１，２であり、それぞれＲ，Ｇ，Ｂの色成分に対応する。すなわち、〔数式２２〕においては、色成分のみを補正する。結局、陰影付加手段６３は、〔数式１６〕～〔数式２２〕に従った処理を実行することにより、各ボクセルを含むｘｙｚ軸方向の近傍に位置する画素の不透明度αを参照して、各ボクセルの色成分の値を補正し、陰影付加を行っている。本実施形態では、好ましい例として、不透明度αを参照する画素として隣接する画素を用いたが、隣接する画素に限定されず、近傍の他の画素を参照するようにしてもよい。
【0139】
ステップＳ６４５において、陰影付加されたボクセルのＲＧＢ値Ｖ（ｃ）が得られたら、累積輝度値Ｅｎｅｒｇｙ（ｃ）と透過光強度Ｔｒａｎｓの更新を行う（ステップＳ６４６）。具体的には、以下の〔数式２３〕に従った処理を実行して、α値を補正した後、累積輝度値Ｅｎｅｒｇｙ（ｃ）と透過光強度Ｔｒａｎｓを更新する。
【0140】
〔数式２３〕
Ａｌｐｈａ←１－（１－Ａｌｐｈａ）1/Sray
Ｅｎｅｒｇｙ（ｃ）←Ｅｎｅｒｇｙ（ｃ）＋Ｔｒａｎｓ・Ａｌｐｈａ・Ｖ（ｃ）／２５５
Ｔｒａｎｓ←Ｔｒａｎｓ・（１．０－Ａｌｐｈａ）
Ａｌｐｈａ←１－（１－Ａｌｐｈａ）1/Sray
Ｅｎｅｒｇｙ（ｃ）←Ｅｎｅｒｇｙ（ｃ）＋Ｔｒａｎｓ・Ａｌｐｈａ・Ｖ（ｃ）／２５５
Ｔｒａｎｓ←Ｔｒａｎｓ・（１．０－Ａｌｐｈａ）
【0141】
累積輝度Ｅｎｅｒｇｙ（ｃ）、透過光強度Ｔｒａｎｓが更新されたら、Ｔｒａｎｓ＜０．００１である場合は、図１３に示したステップＳ６４０の処理を終了する。Ｔｒａｎｓ≧０．００１である場合は、ｚを所定数ｍ１だけ減じる（ステップＳ６４６）。本実施形態では、ｍ１＝１であり、ｍ１＜ｍ２である。そして、ｚが０以上であれば、ステップＳ６４３へ戻って同様の処理を繰り返す。ｚが０未満であれば、対象とするボクセルがなくなるので、図１３に示したステップＳ６４０の処理を終了する。
【0142】
〔数式９〕の第３式に示したように、ステップＳ９０２におけるオフセット処理の際、ｚは仮想光線のサブサンプリング倍率Ｓｒａｙと乗じられる。そのため、実質的には、ｚ軸方向においては、所定の参照間隔（ｍ１・Ｓｒａｙ）ごとにボクセルが参照されてボクセル画素値Ｖ（ｃ）が取得されることになる。したがって、（ｍ１・Ｓｒａｙ）＞１であれば、ｚ軸方向に対しては、全てのボクセルを参照して処理するのではなく、間引いたボクセルを参照して処理を行うことになる。そのため、ｚ軸方向に対して処理するボクセル数が少なくなり効率化が図られるが、ボリュームレンダリングの場合は画質劣化が目立ってくる（信号値でレンダリング像を生成する３Ｄ－ＭＩＰの場合は処理負荷がボリュームレンダリングより大きいため、（ｍ１・Ｓｒａｙ）＞１に設定し、通常はｍ１＝１、Ｓｒａｙ＝２に設定する）。本実施形態では、ｍ１＝１およびＳｒａｙ≦１であることが好ましく、通常はＳｒａｙ＝１に設定する。
【0143】
図１２に示したステップＳ６３０、Ｓ６４０においては、レイキャスティング手段６１は、ボクセル画像をレンダリング画像に投影変換した座標系を視点座標系とすると、視点座標系において視線方向であるｚ軸の上限値を、レンダリング画像の画素（ｘ，ｙ）ごとに算出された視線方向のベクトルと直方体との交点に対応するＺｃに設定し、視線方向であるｚ軸の下限値を、視点座標系のｚ軸方向の下限値である０に設定し、起点座標探索手段６２が、探索開始座標をＺｃとして視点座標系のｚ軸上のＺｃから０の範囲で仮想光線の照射を開始する起点座標Ｚｓｔを探索している。
【0144】
以上のようにして各画素の色成分を得る、図１０に示したステップＳ４４０の処理を、ＮＢ個の各ＣＰＵコアが、ＮＢ個に分割された各分割画像に対して実行する。そして、得られた分割画像を統合することにより、設定されたＳｉｚｅ×ＳｉｚｅのＲＧＢのカラー２次元投影画像であるレンダリング画像が生成される。なお、必ずしも設定されたＳｉｚｅ×Ｓｉｚｅの２次元投影画像に一度に統合する必要はなく、各分割画像を処理が完了次第、順次表示させながら、最終的に全ての分割画像を統合させるようにしてもよい。
【0145】
有効領域の境界面、すなわち有効領域の最外周のボクセルの不透明度αの値を０に設定するとともに色成分のＲＧＢ値には画素の信号値を基にカラーマップで定義される自然なＲＧＢ値を与えることにより、処理負荷を抑制しながら境界面に発生するモアレを抑制することができる。
【0146】
上述のように、本実施形態に係るボリュームレンダリング装置１００では、マルチコアＣＰＵであるＣＰＵ１がプログラムを実行することにより、不透明外接領域算出手段５０、レンダリング手段６０が実現される。このプログラムは、読み込んだ断層画像を複数（例えばＮＡまたはＮＢ個）の断層画像群に分割し、各ＣＰＵスレッドに各断層画像群を割り当てて、各ＣＰＵスレッドが不透明外接領域の一部を算出する処理、分割生成領域に画素値を割り当てる処理を並行して行う。ボリュームレンダリング装置１００を、汎用的なコンピュータであるパソコンで実現した場合について説明する。図１５は、コンピュータによりボリュームレンダリング装置１００を実現する場合の並列化処理を示す図である。図１５は不透明外接領域算出手段５０をＮＡ＝４個の断層画像群に分割して実行させる例を示しているが、レンダリング手段６０をＮＢ＝４個の分割生成領域に分割して実行させる場合も同様な構成になる。図１５に示すように、不透明外接領域算出手段５０、レンダリング手段６０が、アプリケーションプログラムにより実現される。例えばＮＡ＝４個の断層画像群に分割する場合、アプリケーションプログラムは、４個の各断層画像群に対する処理を４つの処理スレッド♯ｎ１～処理スレッド♯ｎ４に分割する。例えばＮＢ＝４個の分割生成領域に分割する場合、アプリケーションプログラムは、４個の各分割生成領域に対する処理を４つの処理スレッド♯ｎ１～処理スレッド♯ｎ４に分割する。そして、不透明外接領域算出手段５０、レンダリング手段６０（アプリケーションプログラム）が稼働するマルチタスク・オペレーティングシステムのジョブスケジューラが、複数のＣＰＵコアが有するＣＰＵスレッド＃１～ＣＰＵスレッド＃８に割り当てる。
【0147】
図１６は、コンピュータによりボリュームレンダリング装置１００を実現する場合の並列化処理のソフトウェアによる流れを示す図である。図１６は不透明外接領域算出手段５０をＮＡ＝４個の断層画像群に分割して実行させる例を示しているが、レンダリング手段６０をＮＢ＝４個の分割生成領域に分割して実行させる場合も同様な構成になる。図１６に示すように、不透明外接領域算出手段５０、レンダリング手段６０を実現するアプリケーションプログラムは、４個の各断層画像群に対する処理、４個の各分割生成領域に対する処理を４つの処理スレッド♯ｎ１～処理スレッド♯ｎ４に分割した後、それぞれ並列処理関数を起動する。そうすると、マルチタスク・オペレーティングシステムのジョブスケジューラが、４つのスレッド♯ｎ１～スレッド♯ｎ４を複数のＣＰＵコアが有するＣＰＵスレッド＃３、ＣＰＵスレッド＃５、ＣＰＵスレッド＃７、ＣＰＵスレッド＃２、にそれぞれ割り当てる。このようにして、オペレーティングシステムにより空いているＣＰＵコアのＣＰＵスレッドに、処理すべき処理スレッドが割り当てられるため、並列処理が可能となる。図１６の例では、アプリケーションプログラムが各ＣＰＵスレッドからの終了メッセージを待ち、全てのＣＰＵスレッドから終了メッセージが得られたら、次の処理に移行する。
【0148】
４コア構成のＣＰＵ(Intel CORE-i7)が実装されているWindows10/64bitsのパーソナルコンピュータで実測した結果、図１６の構成で512×512×370のボクセル画像を用いてレンダリング手段６０をＮＢ＝８個の分割生成領域に分割して陰影計算なし実行させる場合にはシングルスレッド実行時(512×512のフル解像度レンダリング画像のＣＰＵによる生成時間：約１．０５秒)の３．５倍弱程度の処理速度が得られた (同生成時間：約０．３秒)。しかし、一部のＣＰＵスレッドに処理待ちが生じた。そこで、本願提案の鏡面反射を含む陰影を付加し、８処理スレッドで並列処理を行ったところ、シングルスレッド実行時(同生成時間：約１．８秒)のほぼ４倍の処理性能が得られた(同生成時間：約0.45sec)。さらに、ＮＡ＞８に上げて実験した結果、処理性能は殆ど変化せず４倍で頭打ちであり、陰影計算を加えると絶対的な処理時間は増大するが、４コア構成による並列化効果は得られやすかった。尚、上記処理時間には不透明外接領域算出手段５０による処理時間も含まれており、レンダリング手段６０に比べ桁違いに処理時間が短いため、不透明外接領域算出手段５０による単独のマルチスレッド化による並列処理の効果は実測困難であった。
【0149】
＜３．断層画像の階調落とし＞
上記実施形態では、各画素が１６ビットの信号値を記録した断層画像をそのまま用いてレンダリング画像を作成するようにしたが、断層画像を８ビットに階調変換して階調低下画像を作成した後、レンダリング画像を作成するようにしてもよい。断層画像を８ビットに階調変換することにより、各画素の処理における負荷を削減することができ、ボリュームレンダリング像の高速な生成に寄与する。実測の結果、２割程度の削減効果が得られ、４コアＣＰＵ構成で８スレッドで陰影付きレンダリングを実行した場合、１６ビット信号値の断層画像を用いた場合、約０．４５ｓｅｃであったのに対し、８ビット信号値の階調低下断層画像を用いた場合、約０．３６ｓｅｃに短縮され、生成されるレンダリング画像の差異は殆ど認められなかった。
上記実施形態では、各画素が１６ビットの信号値を記録した断層画像をそのまま用いてレンダリング画像を作成するようにしたが、断層画像を８ビットに階調変換して階調低下画像を作成した後、レンダリング画像を作成するようにしてもよい。断層画像を８ビットに階調変換することにより、各画素の処理における負荷を削減することができ、ボリュームレンダリング像の高速な生成に寄与する。実測の結果、２割程度の削減効果が得られ、４コアＣＰＵ構成で８スレッドで陰影付きレンダリングを実行した場合、１６ビット信号値の断層画像を用いた場合、約０．４５ｓｅｃであったのに対し、８ビット信号値の階調低下断層画像を用いた場合、約０．３６ｓｅｃに短縮され、生成されるレンダリング画像の差異は殆ど認められなかった。
【0150】
この場合、断層画像階調変換手段を更に設け、断層画像読込手段１０が単一の断層画像を読み込むごとに、階調変換を行う。具体的には、一連のＤＩＣＯＭ形式の断層画像Ｄｏ（ｘ，ｙ，ｚ）（－３２７６８≦Ｄｏ（ｘ，ｙ，ｚ）≦３２７６７，０≦ｘ≦Ｓｘ－１，０≦ｙ≦Ｓｙ－１，０≦ｚ≦Ｓｚ－１; 解像度：Ｒｘｙ，Ｒｚ）の中で、Ｓｚ／２番目の中間の断層画像Ｄｏ（ｘ，ｙ，ｚ／２）を最初に読み込み、同断層画像における全ての画素の最小値をＤｍｉｎ、最大値をＤｍａｘとする。ＤｍｉｎとＤｍａｘは全ての断層画像より最小値と最大値を求める方法をとっても良いが、中間の断層画像だけで決定する方法をとると、大容量の１６ビットの原断層画像の全てのスライス分を読み込んでメモリ上に保持することなく、階調低下画像である階調圧縮断層画像Ｄ８（ｘ，ｙ，ｚ）だけをメモリ上に直接構築できる。続いて、以下の〔数式２４〕に従った処理を実行して、下限値Ｌｍｉｎ および上限値Ｌｍａｘを算出する。
【0151】
〔数式２４〕
下限値Ｌｍｉｎ＝(Ｄｍａｘ－Ｄｍｉｎ)・γ＋Ｄｍｉｎ
上限値Ｌｍａｘ＝(Ｄｍａｘ－Ｄｍｉｎ)・(１－γ)＋Ｄｍｉｎ
下限値Ｌｍｉｎ＝(Ｄｍａｘ－Ｄｍｉｎ)・γ＋Ｄｍｉｎ
上限値Ｌｍａｘ＝(Ｄｍａｘ－Ｄｍｉｎ)・(１－γ)＋Ｄｍｉｎ
【0152】
〔数式２４〕において、γは階調圧縮画像のコントラスト調整幅で、０に近いほどコントラストが増大するが輝度が小さくなる。γは０．３未満の実数値であるが、通常はγ＝０．１に設定する。レンダリング画像の輝度コントラストはカラーマップなど種々の設定で調整できるので、γは固定で良い。そして、下限値Ｌｍｉｎおよび上限値Ｌｍａｘの範囲を２５６段階に圧縮して階調圧縮断層画像Ｄ８（ｘ，ｙ，ｚ）を得る。具体的には、階調圧縮断層画像Ｄ８（ｘ，ｙ，ｚ）（０≦Ｄ８（ｘ，ｙ，ｚ）≦２５５，０≦ｘ≦Ｓｘ－１，０≦ｙ≦Ｓｙ－１，０≦ｚ≦Ｓｚ－１; 解像度：Ｒｘｙ，Ｒｚ）としたとき、以下の〔数式２５〕に従った処理を実行して、階調圧縮断層画像Ｄ８（ｘ，ｙ，ｚ） を算出する。
【0153】
〔数式２５〕
Ｄ８（ｘ，ｙ，ｚ）＝（Ｄｏ（ｘ，ｙ，ｚ）－Ｌｍｉｎ）・２５５／（Ｌｍａｘ－Ｌｍｉｎ）
Ｄ８（ｘ，ｙ，ｚ）＞２５５の場合：Ｄ８（ｘ，ｙ，ｚ）＝２５５、Ｄ８（ｘ，ｙ，ｚ）＜０の場合：Ｄ８（ｘ，ｙ，ｚ）＝０
Ｄ８（ｘ，ｙ，ｚ）＝（Ｄｏ（ｘ，ｙ，ｚ）－Ｌｍｉｎ）・２５５／（Ｌｍａｘ－Ｌｍｉｎ）
Ｄ８（ｘ，ｙ，ｚ）＞２５５の場合：Ｄ８（ｘ，ｙ，ｚ）＝２５５、Ｄ８（ｘ，ｙ，ｚ）＜０の場合：Ｄ８（ｘ，ｙ，ｚ）＝０
【0154】
〔数式２５〕の第２式に示すように、信号値が上限値Ｌｍａｘを超える場合は２５５、下限値Ｌｍｉｎを下回る場合は０として、下限値Ｌｍｉｎから上限値Ｌｍａｘの範囲を０から２５５に線形変換する。なお、階調圧縮断層画像を用いてボクセル画像を作成する場合、図５のカラーマップに代えて、階調低下画像用のカラーマップとして、８ビットの信号値用のカラーマップを用意しておく。具体的には、信号値０～２５５に対応付けてＲ，Ｇ，Ｂ、αが記録されたカラーマップを用いる。
【0155】
以上、本開示の好適な実施形態について説明したが、本開示は上記実施形態に限定されず、種々の変形が可能である。例えば、上記実施形態では、ＣＰＵ１として、複数のＣＰＵコアを有し、並列処理が可能なマルチコアＣＰＵを用いるようにしたが、１つのＣＰＵコアが１つの処理スレッドの処理のみを行うシングルコアＣＰＵを用いるようにしてもよい。
【0156】
また、上記実施形態では、モアレ対策を行うため、有効領域の最外周のボクセルの不透明度を０とするようにしたが、モアレが目立たない場合は、必ずしも有効領域の最外周のボクセルの不透明度を０とする必要はない。
【符号の説明】
【0157】
１・・・ＣＰＵ（Central Processing Unit）
２・・・ＲＡＭ（Random Access Memory）
３・・・記憶装置
４・・・指示入力Ｉ／Ｆ
５・・・データ入出力Ｉ／Ｆ
６・・・表示部
７・・・ＧＰＵ
８・・・フレームメモリ
１０・・・断層画像読込手段
２０・・・カラーマップ読込手段
３０・・・ＲＯＩクリッピング設定手段
４０・・・マスク設定手段
５０・・・不透明外接領域算出手段
６０・・・レンダリング手段
６１・・・レイキャスティング手段
６２・・・起点座標探索手段
６３・・・陰影付加手段
７０・・・レンダリング画像出力手段
１００・・・ボリュームレンダリング装置
２・・・ＲＡＭ（Random Access Memory）
３・・・記憶装置
４・・・指示入力Ｉ／Ｆ
５・・・データ入出力Ｉ／Ｆ
６・・・表示部
７・・・ＧＰＵ
８・・・フレームメモリ
１０・・・断層画像読込手段
２０・・・カラーマップ読込手段
３０・・・ＲＯＩクリッピング設定手段
４０・・・マスク設定手段
５０・・・不透明外接領域算出手段
６０・・・レンダリング手段
６１・・・レイキャスティング手段
６２・・・起点座標探索手段
６３・・・陰影付加手段
７０・・・レンダリング画像出力手段
１００・・・ボリュームレンダリング装置
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】