JP7205189B2 - レンダリング装置、レンダリング方法、及びプログラム - Google Patents

レンダリング装置、レンダリング方法、及びプログラム Download PDF

Info

Publication number
JP7205189B2
JP7205189B2 JP2018218847A JP2018218847A JP7205189B2 JP 7205189 B2 JP7205189 B2 JP 7205189B2 JP 2018218847 A JP2018218847 A JP 2018218847A JP 2018218847 A JP2018218847 A JP 2018218847A JP 7205189 B2 JP7205189 B2 JP 7205189B2
Authority
JP
Japan
Prior art keywords
voxel
effective
calculated
voxels
pixel
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2018218847A
Other languages
English (en)
Other versions
JP2020086803A (ja
Inventor
敏雄 茂出木
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Dai Nippon Printing Co Ltd
Original Assignee
Dai Nippon Printing Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Dai Nippon Printing Co Ltd filed Critical Dai Nippon Printing Co Ltd
Priority to JP2018218847A priority Critical patent/JP7205189B2/ja
Publication of JP2020086803A publication Critical patent/JP2020086803A/ja
Application granted granted Critical
Publication of JP7205189B2 publication Critical patent/JP7205189B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Generation (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Description

本開示は、レンダリング装置、レンダリング方法、及びプログラムに関する。
医用画像診断においては、Ｘ線ＣＴ（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がｚ＜０の場合（ボクセル空間外の場合）、制御部１１は、ｚ＝－１（起点座標が存在しないことを示す値）を返し（図２５のステップＳ１５４）、処理を終了する。ｚ≧０の場合、制御部１１は、ステップＳ１５２に戻る。
一方、ステップＳ１５２において取得したボクセル信号値ＶｓがＶｓ≧－３２７６８の場合（有効ボクセルの場合）、制御部１１は、ｚを起点座標として出力する。
（１２－２．ボクセル信号値Ｖｓを取得（補間なし））
起点座標探索処理（図２５）のステップＳ１５２において実行される、座標変換を行いボクセル信号値Ｖｓを取得する処理について説明する。
（１２－２－１．座標変換）
まず、前記した（１１－３－１．座標変換（実数の座標値を取得））および（１１－３－２．座標変換（整数の座標値に変換））と同様の方法で、３次元座標値（ｘ、ｙ、ｚ）に対応するボクセルの整数の座標値（ｘｉ、ｙｉ、ｚｉ）を算出する。
（１２－２－２．ボクセル信号値Ｖｓを取得）
そして、制御部１１は、有効ボクセル領域Ｖｒ、マスクデータＭａｓｋ、または補正テーブルＳαが定義されている場合、これらを考慮して、以下のようにボクセル信号値Ｖｓを取得する。
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外の場合）
Ｖｓ＝－３２７６９
２）上記１）を満たさず、Ｍａｓｋ（ｘ、ｙ、ｚ）＝０の場合（無効値）
Ｖｓ＝－９９９９９
３）上記１）２）を満たさない場合（有効値）
Ｖｓ＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
（１２－３．ボクセル信号値Ｖｓの取得（補間あり））
（１２－３－１．座標変換）
まず、前記した（１１－３－１．座標変換（実数の座標値を取得））と同様の方法で、３次元座標値（ｘ、ｙ、ｚ）に対応するボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を算出する。
（１２－３－２．ボクセル信号値の抽出）
続いて、前記した（１０－２．ボクセル信号値の抽出）と同様の方法で、補間対象のボクセル信号値Ｄ０００、Ｄ１００、Ｄ０１０、Ｄ００１、Ｄ１０１、Ｄ０１１、Ｄ１１１を抽出し、α値スケールＳ０００、Ｓ１００、Ｓ０１０、Ｓ１１０、Ｓ００１、Ｓ１０１、Ｓ０１１、Ｓ１１１を設定する。
（１２－３－３．ボクセル信号値Ｖｓの取得）
そして、制御部１１は、抽出した補間対象ボクセルの信号値に基づいて、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセル信号値Ｖｓを以下のように決定する。
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｖｓ＝－３２７６９
２）上記１）の条件を満たさず、Ｓ０００＝０の場合（無効値）
Ｖｓ＝－９９９９９
３）上記１）２）を満たさない場合において、
ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、またはｚｉ＋１＞Ｚｉｅのいずれかを満たすか、或いはＳ０００、Ｓ１００、Ｓ０１０、Ｓ１１０、Ｓ００１、Ｓ１０１、Ｓ０１１、Ｓ１１１のいずれかが０の場合（補間しない）
Ｖｓ＝Ｄ０００
４）上記１）２）３）の条件を満たさない場合（補間する）
Ｖｓ＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｄ０００＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｄ１００＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｄ０１０＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｄ１１０＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｄ００１＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｄ１００＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｄ０１１＋ｗｚ・ｗｙ・ｗｘ・Ｄ１１１
以上、レンダリング処理として、レイキャスティング法に基づいてＭＩＰ画像を生成する処理について説明した。
以上、本開示の実施の形態について説明した。本開示の実施の形態によれば、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する。そして、視線方向のベクトルと有効ボクセル領域との交点を幾何計算により算出し、算出された交点から視線方向に向かってレイキャスティング計算を開始する起点座標を探索する。これにより、有効ボクセル領域の外側における冗長な探索処理を削減することができ、レンダリング処理の高速化が実現できる。例えば、有効ボクセル領域との交点が存在しない背景部の画素は計算対象から外され、冗長な処理が大幅に削減される。また、背景部以外の画素についても、交点から起点座標の探索を開始すればよいため、投影面から交点までの範囲の探索を省略できる。本開示の方法は、レイキャスティング計算を開始する起点座標を探索するにあたり、投影面から有効ボクセルとの距離の分布を示すデプスマップを作成したりする必要がないので、レンダリング処理の高速化を簡便な手段で実現することができる。
また、起点座標の探索範囲が絞り込まれるので、従来に比べて起点座標を探索する際のスキップ幅を小さくすることができる。これにより、従来に比べて起点座標をより正確に求めることができ、精度の高いレンダリング画像を得ることができる。
図２６～図２８は、スキップ幅を変えた場合のレンダリング画像の画像精度を比較した図である。図２６はスキップ幅＝１６、図２７はスキップ幅＝８、図２８はスキップ幅＝４とした場合のレンダリング画像を示す。（ａ）はフルカラー（２４ビット）のボリュームレンダリング画像、（ｂ）はモノクロ（８ビット）のボリュームレンダリング画像とＭＩＰ画像を示す。なお、図２８（ｂ）のＭＩＰ画像のみスキップ幅を１としている。
図２６、図２７のように、スキップ幅をある程度大きくするとレンダリング処理を高速に実行できるが、有効ボクセルの見落としが多くなるため、画素抜け（図２６（ａ）、図２７（ａ）参照）や画像エッジ部の波形歪み（図２６（ｂ）、図２７（ｂ）参照）が発生する。これに対し、図２８のようにスキップ幅を小さくすると、有効ボクセルの見落としが少なくなるため、画素抜け（図２８（ａ）参照）や画像エッジ部の波形歪み（図２８（ｂ）参照））の発生が抑えられる。
しかしながら、従来（図３（ａ）参照）のように、投影面から起点座標を探索する方法でスキップ幅を小さくすると、計算負荷が大きくなり、現実的な処理時間でレンダリング処理を実行することが困難となる。この点、本開示（図３（ｂ）参照）のように、視線ベクトルと有効ボクセル領域との交点を算出し、この交点から起点座標を探索する方法をとることで、図２８のようにスキップ幅を小さく設定しても現実的な処理時間でレンダリング処理を実行することができる。すなわち、本開示の方法によれば、レンダリング処理を高速かつ高精度に実行することが可能となる。
さらに本開示では、レイキャスティング計算の途上において、有効ボクセル領域を逸脱（通過）した場合（図１９のステップＳ９６；Ａｌｐｈａ＜０、図２４のステップＳ１３５；Ｖｓ＝－３２７６９）、その先に描画対象のボクセルは存在しないため、レイキャスティング処理を早期に打ち切るようにする。これにより、有効ボクセル領域の外側におけるレイキャスティング計算を省略できる。
また本開示では、断層画像群Ｄｏから構成されるボクセル空間を分割し、分割したボクセル空間ごとに有効ボクセル領域を並行して算出し、算出した各有効ボクセル領域を統合する方法をとることで高速に有効ボクセル領域を求めることができる。また、レンダリングの投影面を複数の画素領域に分割し、各画素領域におけるレンダリング処理を並行して実行させることで、高速にレンダリング画像を得ることができる。例えば、並行化のスレッド数をＮ（例えばＮ＝４）とし、レンダリング装置１で使用可能なＣＰＵコアがＮ個以上存在し、Ｎ個のスレッドに異なるＣＰＵコアが割り当てられれば、Ｎ倍の処理速度で有効ボクセル領域算出処理およびレンダリング処理を実行することが可能となる。
また本開示では、上記した並行化処理を１台のコンピュータ（レンダリング装置１）で完結させることで、リモートデクストップ等の構成で複数のユーザが使用しても、サーバーコンピュータ間の通信が発生しないので、通信経路上でデータが喪失（データロス）したり、ネットワーク遅延等に起因するレンダリング処理の遅延が生じない。また、機密性の高い断層画像群Ｄｏ（ＤＩＣＯＭ画像）が複数のクライアントコンピュータにコピーされるようなことがないため、セキュリティ対策や個人情報の取り扱いが容易となる。
なお、上記した並行化処理は必須ではない。この場合、図１０のステップＳ３において、制御部１１は、ボクセル空間を分割することなく、ボクセル空間全体における有効ボクセル領域を直接算出する。また、図１０のステップＳ４において、制御部１１は、投影面を分割することなく、投影面全体を対象としてレンダリング処理（レイキャスティング処理）を実行する。
以上、添付図面を参照しながら、本開示に係るレンダリング装置１等の好適な実施形態について説明したが、本開示はかかる例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本開示の技術的範囲に属するものと了解される。
１ ：レンダリング装置
Ｄｏ ：断層画像群
Ｃｍａｐ ：カラーマップ
ＲＯＩ ：クリッピング領域
Ｍａｓｋ ：マスクデータ
Ｓα ：補正テーブル
Ｉｍａｇｅ ：レンダリング画像(ボリュームレンダリング画像、ＭＩＰ画像)
Ｖｒ ：有効ボクセル領域

Claims (31)

  1. 複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
    を備え、
    クリッピング領域が設定されている場合には、
    前記有効ボクセル領域算出手段は、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記クリッピング領域に含まれるボクセルを、有効ボクセルとする、
    レンダリング装置。
  2. 複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
    を備え、
    マスク領域が設定されている場合には、
    前記有効ボクセル領域算出手段は、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記マスク領域に含まれるボクセルまたは前記マスク領域に含まれないボクセルを、有効ボクセルとする、
    レンダリング装置。
  3. 複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
    を備え、
    前記有効ボクセル領域算出手段は、
    前記複数の断層画像を複数の断層画像群に分割する手段と、
    分割した断層画像群ごとに、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たす有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する手段と、
    算出した各有効ボクセル領域に基づいてボクセル空間全体における有効ボクセル領域を算出する手段と、を備え、
    複数のＣＰＵコアを用いて、各断層画像群に対応する有効ボクセル領域を並行して算出する、
    レンダリング装置。
  4. 複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
    を備え、
    前記有効ボクセル領域算出手段は、
    信号値が所定の閾値以上のボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
    前記レンダリング手段は、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上において所定の閾値以上の信号値をもつボクセルを抽出し、抽出したボクセルの信号値に基づいて代表信号値を算出し、レンダリング画像の画素値とする、
    レンダリング装置。
  5. 前記レンダリング手段は、
    前記代表信号値として、抽出したボクセルの信号値の最大値、最小値または平均値のいずれかを算出する、
    請求項４に記載のレンダリング装置。
  6. 複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
    を備え、
    前記有効ボクセル領域算出手段は、
    不透明度が０でないボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
    前記レンダリング手段は、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルの色値および不透明度に基づいて、各画素の画素値を算出する、
    レンダリング装置。
  7. 複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
    を備え、
    前記レンダリング手段は、
    視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
    仮想光線上のボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
    ボクセル座標系における座標に基づいて近傍の複数のボクセルを特定し、
    特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、前記複数のボクセルのボクセル値に基づいて、前記画素値を算出する、
    レンダリング装置。
  8. 前記レンダリング手段は、
    特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、特定した複数のボクセルの各信号値に所定の重みを付けて加算し、前記画素値を算出する、
    請求項７に記載のレンダリング装置。
  9. 前記レンダリング手段は、
    特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、特定した複数のボクセルの各色値および各不透明度にそれぞれ所定の重みを付けて加算し、前記画素値を算出する、
    請求項７に記載のレンダリング装置。
  10. 複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段と、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段と、
    を備え、
    前記レンダリング手段は、
    レンダリング画像の画素ごとに、算出された視点に近い交点を探索開始点として、視線方向に向かって最も近い有効ボクセルの座標を、仮想光線の照射を開始する起点座標として探索する起点座標探索手段と、を更に備える、
    レンダリング装置。
  11. 前記起点座標探索手段は、
    算出された視点から遠い交点までの範囲を探索する、
    請求項１０に記載のレンダリング装置。
  12. 前記起点座標探索手段は、
    レンダリング画像の画素ごとに、視線方向に向かって所定のステップ幅で当該画素に対応するボクセルのボクセル値を順次取得しながら、有効ボクセルを探索する、
    請求項１０または請求項１１に記載のレンダリング装置。
  13. 前記起点座標探索手段により有効ボクセルが探索されなかった場合、予め決められた所定値を返す、
    請求項１２に記載のレンダリング装置。
  14. 前記起点座標探索手段は、
    有効ボクセルが探索された場合、当該有効ボクセルの位置と前記ステップ幅だけ視線方向と逆の方向に戻った位置との間に存在する有効ボクセルを探索する、
    請求項１２に記載のレンダリング装置。
  15. 前記起点座標探索手段は、
    視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
    ボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
    ボクセル座標系における座標が有効ボクセル領域に含まれる場合は、当該ボクセルの信号値をボクセル値として取得する、
    請求項１２から請求項１４のいずれかに記載のレンダリング装置。
  16. 前記起点座標探索手段は、
    視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
    ボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
    ボクセル座標系における座標が有効ボクセル領域に含まれる場合は、当該ボクセルの不透明度をボクセル値として取得する、
    請求項１２から請求項１４のいずれかに記載のレンダリング装置。
  17. 前記レンダリング手段は、
    レンダリング画像の投影面を複数の画素領域に分割する手段と、
    分割した画素領域ごとに、当該画素領域の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を求めるレンダリング処理を実行する手段と、
    算出した各画素領域のレンダリング結果に基づいてレンダリング画像を生成する手段と、を備え、
    複数のＣＰＵコアを用いて、各画素領域に対する各レンダリング処理を並行して実行する、
    請求項１から請求項１６のいずれかに記載のレンダリング装置。
  18. コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
    を含み、
    クリッピング領域が設定されている場合には、
    前記有効ボクセル領域算出ステップは、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記クリッピング領域に含まれるボクセルを、有効ボクセルとする、
    レンダリング方法。
  19. コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
    を含み、
    マスク領域が設定されている場合には、
    前記有効ボクセル領域算出ステップは、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記マスク領域に含まれるボクセルまたは前記マスク領域に含まれないボクセルを、有効ボクセルとする、
    レンダリング方法。
  20. コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
    を含み、
    前記有効ボクセル領域算出ステップは、
    前記複数の断層画像を複数の断層画像群に分割するステップと、
    分割した断層画像群ごとに、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たす有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出するステップと、
    算出した各有効ボクセル領域に基づいてボクセル空間全体における有効ボクセル領域を算出するステップと、を含み、
    複数のＣＰＵコアを用いて、各断層画像群に対応する有効ボクセル領域を並行して算出する、
    レンダリング方法。
  21. コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
    を含み、
    前記有効ボクセル領域算出ステップは、
    信号値が所定の閾値以上のボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
    前記レンダリングステップは、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上において所定の閾値以上の信号値をもつボクセルを抽出し、抽出したボクセルの信号値に基づいて代表信号値を算出し、レンダリング画像の画素値とする、
    レンダリング方法。
  22. コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
    を含み、
    前記有効ボクセル領域算出ステップは、
    不透明度が０でないボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
    前記レンダリングステップは、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルの色値および不透明度に基づいて、各画素の画素値を算出する、
    レンダリング方法。
  23. コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
    を含み、
    前記レンダリングステップは、
    視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
    仮想光線上のボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
    ボクセル座標系における座標に基づいて近傍の複数のボクセルを特定し、
    特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、前記複数のボクセルのボクセル値に基づいて、前記画素値を算出する、
    レンダリング方法。
  24. コンピュータが、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング方法であって、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出ステップと、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリングステップと、
    を含み、
    前記レンダリングステップは、
    レンダリング画像の画素ごとに、算出された視点に近い交点を探索開始点として、視線方向に向かって最も近い有効ボクセルの座標を、仮想光線の照射を開始する起点座標として探索する起点座標探索ステップと、を更に含む、
    レンダリング方法。
  25. コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
    前記コンピュータを、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
    として機能させるプログラムであって、
    クリッピング領域が設定されている場合には、
    前記有効ボクセル領域算出手段は、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記クリッピング領域に含まれるボクセルを、有効ボクセルとする、
    プログラム。
  26. コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
    前記コンピュータを、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
    として機能させるプログラムであって、
    マスク領域が設定されている場合には、
    前記有効ボクセル領域算出手段は、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たし、かつ前記マスク領域に含まれるボクセルまたは前記マスク領域に含まれないボクセルを、有効ボクセルとする、
    プログラム。
  27. コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
    前記コンピュータを、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
    として機能させるプログラムであって、
    前記有効ボクセル領域算出手段は、
    前記複数の断層画像を複数の断層画像群に分割する手段と、
    分割した断層画像群ごとに、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たす有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する手段と、
    算出した各有効ボクセル領域に基づいてボクセル空間全体における有効ボクセル領域を算出する手段と、を備え、
    複数のＣＰＵコアを用いて、各断層画像群に対応する有効ボクセル領域を並行して算出する、
    プログラム。
  28. コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
    前記コンピュータを、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
    として機能させるプログラムであって、
    前記有効ボクセル領域算出手段は、
    信号値が所定の閾値以上のボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
    前記レンダリング手段は、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上において所定の閾値以上の信号値をもつボクセルを抽出し、抽出したボクセルの信号値に基づいて代表信号値を算出し、レンダリング画像の画素値とする、
    プログラム。
  29. コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
    前記コンピュータを、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
    として機能させるプログラムであって、
    前記有効ボクセル領域算出手段は、
    不透明度が０でないボクセルを有効ボクセルとし、前記有効ボクセル領域を算出し、
    前記レンダリング手段は、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルの色値および不透明度に基づいて、各画素の画素値を算出する、
    プログラム。
  30. コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
    前記コンピュータを、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
    として機能させるプログラムであって、
    前記レンダリング手段は、
    視点座標系における回転を定義した回転行列、ＸＹＺ軸各方向のオフセット値、ＸＹＺ軸方向の拡大又は縮小倍率、Ｚ軸方向の変倍率を含む座標変換パラメータを取得し、
    仮想光線上のボクセルの視点座標系における座標を、前記座標変換パラメータに基づいて、ボクセル座標系における座標に変換し、
    ボクセル座標系における座標に基づいて近傍の複数のボクセルを特定し、
    特定した複数のボクセルのボクセル座標系における座標が前記有効ボクセル領域に含まれる場合は、前記複数のボクセルのボクセル値に基づいて、前記画素値を算出する、
    プログラム。
  31. コンピュータを、複数の断層画像から構成される３次元ボクセルに基づいてレンダリング画像を生成するレンダリング装置として機能させるプログラムであって、
    前記コンピュータを、
    ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセルを有効ボクセルとし、全ての有効ボクセルを含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域として算出する有効ボクセル領域算出手段、
    レンダリング画像の各画素に対して、視線方向のベクトルと前記有効ボクセル領域との交点を算出し、算出された視点に近い交点から視線方向に向かって仮想光線を照射し、仮想光線上のボクセルのボクセル値に基づいて、各画素の画素値を算出するレンダリング手段、
    として機能させるプログラムであって、
    前記レンダリング手段は、
    レンダリング画像の画素ごとに、算出された視点に近い交点を探索開始点として、視線方向に向かって最も近い有効ボクセルの座標を、仮想光線の照射を開始する起点座標として探索する起点座標探索手段と、を更に備える、
    プログラム。
JP2018218847A 2018-11-22 2018-11-22 レンダリング装置、レンダリング方法、及びプログラム Active JP7205189B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2018218847A JP7205189B2 (ja) 2018-11-22 2018-11-22 レンダリング装置、レンダリング方法、及びプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2018218847A JP7205189B2 (ja) 2018-11-22 2018-11-22 レンダリング装置、レンダリング方法、及びプログラム

Publications (2)

Publication Number Publication Date
JP2020086803A JP2020086803A (ja) 2020-06-04
JP7205189B2 true JP7205189B2 (ja) 2023-01-17

Family

ID=70908208

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2018218847A Active JP7205189B2 (ja) 2018-11-22 2018-11-22 レンダリング装置、レンダリング方法、及びプログラム

Country Status (1)

Country Link
JP (1) JP7205189B2 (ja)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002306480A (ja) 2001-04-12 2002-10-22 Toshiba Corp 画像処理装置及びその方法
JP6418344B1 (ja) 2018-02-23 2018-11-07 大日本印刷株式会社 コンピュータプログラム、画像処理装置及び画像処理方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH08147496A (ja) * 1994-11-21 1996-06-07 Hitachi Medical Corp 回転像の作成方法および装置
JP3527796B2 (ja) * 1995-06-29 2004-05-17 株式会社日立製作所 高速３次元画像生成装置および方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002306480A (ja) 2001-04-12 2002-10-22 Toshiba Corp 画像処理装置及びその方法
JP6418344B1 (ja) 2018-02-23 2018-11-07 大日本印刷株式会社 コンピュータプログラム、画像処理装置及び画像処理方法

Also Published As

Publication number Publication date
JP2020086803A (ja) 2020-06-04

Similar Documents

Publication Publication Date Title
US10546415B2 (en) Point cloud proxy for physically-based volume rendering
US10893262B2 (en) Lightfield rendering based on depths from physically-based volume rendering
US10692272B2 (en) System and method for removing voxel image data from being rendered according to a cutting region
US20050143654A1 (en) Systems and methods for segmented volume rendering using a programmable graphics pipeline
JP4588736B2 (ja) 画像処理方法および装置並びにプログラム
JP5225999B2 (ja) 結合された強度の投影
US20050237336A1 (en) Method and system for multi-object volumetric data visualization
JP2006502508A (ja) ３ｄモデリングシステム
US8149237B2 (en) Information processing apparatus and program
JP2016135252A (ja) 医用画像処理装置及び医用画像診断装置
JP2006055213A (ja) 画像処理装置、及びプログラム
JP2017076397A (ja) メッシュのボクセル化
JP4885042B2 (ja) 画像処理方法および装置ならびにプログラム
US20220343589A1 (en) System and method for image processing
JP7275669B2 (ja) マスク生成装置、３次元再構成像生成装置、マスク生成方法、及びプログラム
JP2006000127A (ja) 画像処理方法および装置並びにプログラム
JP6544472B1 (ja) レンダリング装置、レンダリング方法、及びプログラム
JP4376944B2 (ja) 中間画像生成方法および装置ならびにプログラム
JP7247577B2 (ja) ３次元再構成像表示装置、３次元再構成像表示方法、プログラム、及び画像生成方法
JP7205189B2 (ja) レンダリング装置、レンダリング方法、及びプログラム
JP5065740B2 (ja) 画像処理方法および装置並びにプログラム
JP7003635B2 (ja) コンピュータプログラム、画像処理装置及び画像処理方法
JP2023527438A (ja) リアルタイム深度マップを用いたジオメトリ認識拡張現実効果
JP7167699B2 (ja) ボリュームレンダリング装置
JP2006000338A (ja) 画像処理方法および装置並びにプログラム

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20210928

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20220622

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20220628

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20220825

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

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20221129

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20221212

R150 Certificate of patent or registration of utility model

Ref document number: 7205189

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150