JP7275669B2 - マスク生成装置、３次元再構成像生成装置、マスク生成方法、及びプログラム - Google Patents

Description

本開示は、マスク生成装置、３次元再構成像生成装置、マスク生成方法、及びプログラムに関する。

医用画像診断において、特定の人体組織の観察に適した３次元再構成像等を得たい場面がある。例えば、胸部や頭部にある内臓や脳等の所定の人体組織の観察を行いたい場合に、内臓や脳等は骨に囲まれており、骨領域はむしろ診断の妨げになる。ＣＴ画像から生成される３次元再構成像等では一般に骨が鮮明に描写され、内臓や血管は隠れてしまうことがあるため、可視化にあたっては切断を行うなど工夫が必要となることがある。

たとえば、領域拡張法（リージョングローイング法）を用いて非観察対象としたい骨領域を抽出した３次元マスクを作成し、３次元マスクを参照しながらマスク処理により非観察対象としたい骨領域を除去したボリュームレンダリング像を生成する方法が開示されている（特許文献２、３参照）。領域拡張法とは、非観察対象領域の画素を指定し、その画素を開始点（拡張開始点）として、近傍画素を次々と抽出する方法である。

特開２００９－２４０５６９号公報 特許４０８７５１７号公報 特許５２５７９５８号公報

しかしながら、領域拡張法を用いた３次元マスクの作成は、ユーザによる複数の拡張開始点の指示が必須で、拡張の打ち切り段階もユーザが指示する必要があるため、自動化が難しく、ユーザごとに結果にバラつきが発生するという問題がある。また領域拡張法に限らず、３次元マスクの作成は、ユーザのスキルや経験が必要であり、その作成に時間や手間がかかる。

本開示の実施形態は、このような点を鑑みてなされたものであり、所望の領域を可視化するためのマスクデータを容易に作成することが可能な、マスク生成装置等を提供することを目的とする。

本開示の一実施形態によると、複数の断層画像を取得する画像取得手段と、断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段と、断層画像毎に、断層画像の各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、を備えるマスク生成装置が提供される。
本開示の一実施形態によると、複数の断層画像を取得する画像取得手段と、断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、断層画像毎に、断層画像の各画素が前記幾何形状の内部か外部のいずれに含まれるかを定義したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、を備え、前記被写体領域算出手段は、前記被写体領域として、信号値が所定の閾値より大きい画素が存在する領域の最外側を囲う矩形領域である第１の矩形領域を算出する、マスク生成装置が提供される。
本開示の一実施形態によると、複数の断層画像を取得する画像取得手段と、断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、断層画像毎に、断層画像の各画素が前記幾何形状の内部か外部のいずれに含まれるかを定義したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、を備え、前記被写体領域算出手段は、前記被写体領域として、信号値が所定の閾値より大きい画素が存在する領域の最外側を囲う矩形領域である第１の矩形領域の内部に定義される第２の矩形領域を算出する、マスク生成装置が提供される。
本開示の一実施形態によると、複数の断層画像を取得する画像取得手段と、断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、断層画像毎に、断層画像の各画素が前記幾何形状の内部か外部のいずれに含まれるかを定義したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、断層画像毎に、前記被写体領域算出手段により算出された被写体領域の縦方向のサイズと横方向のサイズの比率に応じた補正方法により該被写体領域を補正する被写体領域補正手段と、を備えるマスク生成装置が提供される。

本開示の一実施形態によると、複数の断層画像を取得する画像取得手段と、断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段と、断層画像毎に、断層画像の各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、前記断層画像の各画素をボクセル空間に配置して構成される３次元ボクセルと前記マスクデータに基づいて３次元再構成像を生成する３次元再構成像生成手段と、を備える３次元再構成像生成装置が提供される。

本開示の一実施形態によると、コンピュータが、複数の断層画像を取得する画像取得ステップと、断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出ステップと、断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出ステップと、断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出ステップと、断層画像毎に、断層画像の各画素に対して重み値算出ステップにより算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成ステップと、を実行するマスク生成方法が提供される。

本開示の一実施形態によると、コンピュータを、複数の断層画像を取得する画像取得手段、断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段、断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段、断層画像毎に、断層画像の各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段、として機能させるプログラムが提供される。

本開示により、所望の領域を可視化するためのマスクデータを容易に作成することができる。

マスク生成装置１の機能構成を示す図 マスク生成装置１のハードウェア構成を示す図 マスク生成装置１の動作を示すフローチャート 描画重みテーブル作成処理の流れを示すフローチャート 矩形領域４１、矩形領域４２、被写体候補領域５０を説明する図 被写体領域算出処理の流れを示すフローチャート 被写体領域算出処理の概要を示す図 第２重心の算出について説明する図 第３重心の算出について説明する図 算出された矩形領域４２を示す図 （ａ）領域内部の信号値が比較的大きい場合の矩形領域４２、（ｂ）領域内部の信号値が比較的小さい場合の矩形領域４２ 被写体候補領域５０との面積比率に基づいて矩形領域４２を補正する内容を示す図 面積比率に基づく補正を上胸部に適用した場合を示す図 矩形領域４２の縦横比率に基づいて矩形領域４２を補正（縮小補正）する内容を示す図 被写体領域補正処理の流れを示すフローチャート 重み値Ｓｖ（ｘ、ｙ、ｚ）の２次元分布パターンを示す図 マスクデータ生成処理の流れを示すフローチャート マスクデータ生成処理の流れを示すフローチャート 描画重みテーブル作成処理およびマスクデータ作成処理を並行して実行する概要を示す図 描画重みテーブル作成処理およびマスクデータ作成処理を並行して実行させる仕組みを示す図 ３次元再構成像生成装置２が実行する再構成像生成処理の概要を示す図 ３次元再構成像生成装置２の機能構成を示す図 ３次元再構成像生成装置２のハードウェア構成を示す図 ３次元再構成像生成装置２の動作を示すフローチャート ボリュームレンダリング画像を生成する処理の全体の流れを示すフローチャート 有効ボクセル領域Ｖｒを示す図 ボリュームレンダリング画像を生成する具体的な処理の流れを示すフローチャート レイキャスティング処理の流れを示すフローチャート 有効ボクセル領域Ｖｒと視線ベクトルとの交点を算出する様子を示す図 有効ボクセル領域Ｖｒと視線ベクトルとの交点を算出する様子を示す図 起点座標探索処理の流れを示すフローチャート ボクセルα値取得処理（補間なし）の流れを示すフローチャート ボクセルα値取得処理（補間あり）の流れを示すフローチャート ＭＩＰ画像を生成する処理の全体の流れを示すフローチャート 有効ボクセル領域Ｖｒを示す図 ＭＩＰ画像を生成する具体的な処理の流れを示すフローチャート レイキャスティング処理の流れを示すフローチャート 起点座標探索処理の流れを示すフローチャート ボクセル信号値取得処理（補間なし）の流れを示すフローチャート ボクセル信号値取得処理（補間あり）の流れを示すフローチャート ＭＰＲ画像を生成する処理の流れを示すフローチャート ＭＰＲ画像（体軸断面、冠状断面、矢状断面）について説明する図 胸腹部のボリュームレンダリング画像の表示例 胸腹部のＭＩＰ画像の表示例 胸腹部のＭＰＲ画像（体軸断面、冠状断面、矢状断面）の表示例 胸腹部のボリュームレンダリング画像の比較表示例

以下図面に基づいて、本開示の実施の形態を詳細に説明する。

（１．マスク生成装置１の機能構成）
図１は、マスク生成装置１の機能構成を示す図である。図に示すように、マスク生成装置１は、画像取得部２１、描画重みテーブル作成部２２、マスクデータ生成部２３を備える。

画像取得部２１は、ＣＴ装置により撮影された断層画像群Ｄｏを取得する。断層画像群Ｄｏは、被検体（人体）を頭尾軸に沿って所定の間隔で連続的に撮影した複数の断層画像（ＣＴ画像）である。各断層画像はＤＩＣＯＭ形式の２次元の画像データである。ＤＩＣＯＭ形式は、１ファイルにヘッダ部と画像データ部を含む医療画像で一般的に用いられる画像フォーマットであり、画像撮影時のパラメータや診断情報を保存しておくことができる。

１つの断層画像は、例えば、５１２×５１２ピクセルの画像（ＣＴ画像）である。断層画像の各画素には、信号値ｖが付与されており、ＣＴ画像の場合、信号値ｖはＣＴ値である。本実施の形態では、信号値ｖ（ＣＴ値）は１６ビット（－３２７６８≦ｖ≦３２７６７）のデータとする（但し、信号値ｖのビット数は特に限定されない）。

ＣＴ値は、水を基準として表現した組織のＸ線減弱係数であり、ＣＴ値により組織や病変の種類等が判断できるようになっている(単位はＨＵ（Hounsfield Unit））。ＣＴ値は、水と空気のＸ線減弱係数で標準化されており、水のＣＴ値は０、空気のＣＴ値を－１０００である。また、脂肪のＣＴ値は－１２０～－１００程度であり、通常組織のＣＴ値は０～１２０程度であり、骨のＣＴ値は１０００前後を示す。

断層画像群Ｄｏは、ＸＹの２次元データである断層画像をＺ軸方向に積層したものであり、ＸＹＺの座標軸で規定されるボクセル空間Ｒに配置可能なボクセルデータ（３次元ボクセル）として表現可能である。例えば、断層画像群Ｄｏ（３次元ボクセル）は、以下のように定義される。なお本開示において、Ｘ軸は人体の左右軸、Ｙ軸は人体の背腹軸（前後軸）、Ｚ軸は人体の頭尾軸（上下軸）に相当するものとする。ＸＹ軸は、例えば、ＣＴ装置での撮影時に設定される。また、本開示において、Ｘ軸方向を横方向、Ｙ軸方向を縦方向と呼ぶ場合がある。

（式１）
－３２７６８≦Ｄｏ（ｘ、ｙ、ｚ）≦３２７６７
０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１、０≦ｚ≦Ｓｚ－１
解像度：Ｒｘｙ、解像度Ｒｚ

（式１）において、Ｓｘは断層画像の横方向の画素数（ボクセル数）、Ｓｙは断層画像の縦方向の画素数（ボクセル数）、Ｓｚはスライス枚数（ボクセル数）を表し、スキャン方向に所定の間隔で配置される。Ｒｘｙは横方向及び縦方向の断層画像の解像度であり、横方向と縦方向の解像度は通常同一で（ＤＩＣＯＭ規格上は異なった解像度を指定できるが、異なった解像度でスキャンするモダリティは現存しない）、画素の間隔（Ｗｘｙ）の逆数、すなわち単位距離あたりの画素数を示す。Ｒｚはスライスの解像度であり、断層画像のスライス間隔（Ｗｚ）（例えば、０．５ｍｍや１ｍｍ）の逆数、すなわち単位距離あたりのスライス枚数を表す。

断層画像群Ｄｏの信号値は、必要に応じて８ビットに階調圧縮される場合がある。８ビットに階調圧縮された断層画像群Ｄｏ（ｘ、ｙ、ｚ）は以下のように定義される。

（式２）
０≦Ｄｏ（ｘ、ｙ、ｚ）≦２５５
０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１、０≦ｚ≦Ｓｚ－１

描画重みテーブル作成部２２は、断層画像Ｄｏ（ｚ）毎に、断層画像Ｄｏ（ｚ）の各画素（ｘ、ｙ、ｚ）を描画対象とする重み値を格納した描画重みテーブルＳｖを作成する。描画重みテーブルＳｖは以下のように定義される。

（式３）
０≦Ｓｖ（ｘ、ｙ、ｚ）≦１
０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１、０≦ｚ≦Ｓｚ－１

描画重みテーブル作成部２２は、図１に示すように、被写体領域算出部３１、被写体領域補正部３２、被写体形状算出部３３、重み値算出部３４から構成される。

被写体領域算出部３１は、断層画像Ｄｏ（ｚ）毎に、断層画像Ｄｏ（ｚ）の各画素の信号値に基づいて被写体領域を算出する。

被写体領域補正部３２は、胸部ＣＴ画像の場合に、被写体領域算出部３１により算出された被写体領域を補正する。

被写体形状算出部３３は、断層画像Ｄｏ（ｚ）毎に、被写体領域算出部３１により算出された被写体領域または被写体領域補正部３２により補正された被写体領域の外郭を所定の幾何形状（被写体形状）で近似し、当該幾何形状を規定するパラメータ（幾何パラメータＰ（ｚ））を算出する。

重み値算出部３４は、断層画像Ｄｏ（ｚ）毎に、被写体形状算出部３３により算出された幾何パラメータＰ（ｚ）に基づいて、断層画像Ｄｏ（ｚ）の各画素（ｘ、ｙ、ｚ）を描画対象とする重み値Ｓｖ（ｘ、ｙ、ｚ）を算出し、重み値Ｓｖ（ｘ、ｙ、ｚ）を格納した３次元の描画重みテーブルＳｖを作成する。

なお、重み値Ｓｖ（ｘ、ｙ、ｚ）は３次元のデータ配列であり、そのままデータテーブル（描画重みテーブル）として扱えるため、重み値と描画重みテーブルは実質的に同一である。このため、重み値と描画重みテーブルは同一の符号「Ｓｖ」で表す。

マスクデータ生成部２３は、描画重みテーブル作成部２２により作成された描画重みテーブルＳｖに基づいて、断層画像群Ｄｏの各画素（ｘ、ｙ、ｚ）の重み値を所定の閾値で二値化したマスク値を算出し、各画素（ｘ、ｙ、ｚ）に対応するマスク値を保持する３次元のマスクデータＭａｓｋを生成する。マスクデータＭａｓｋは以下のように定義される。

（式４）
Ｍａｓｋ（ｘ、ｙ、ｚ）＝ ０（非描画）または１（描画）
０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１、０≦ｚ≦Ｓｚ－１

（２．マスク生成装置１のハードウェア構成）
図２は、本実施の形態におけるマスク生成装置１のハードウェア構成を示すブロック図である。図に示すように、マスク生成装置１は、制御部１１、記憶部１２、メディア入出力部１３、通信制御部１４、入力部１５、表示部１６、周辺機器Ｉ／Ｆ部１７等が、バス１８を介して接続される汎用のコンピュータで実現される。但し、これに限ることなく、用途、目的に応じて様々な構成を採ることが可能である。

制御部１１は、単一プロセッサ、或いはマルチコアＣＰＵで構成されるプロセッサ又はプロセッサが複数個で構成されるマルチプロセッサ（以降、いずれも単に「ＣＰＵ」と表記し、ＣＰＵが有するコアをＣＰＵコアと表記する）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フレームメモリ(Frame Memory)等によって構成される。ＣＰＵは、記憶部１２、ＲＯＭ、記録媒体等に格納されるレンダリングプログラムをＲＡＭ、フレームメモリ上のワークメモリ領域に呼び出して実行し、バス１８を介して接続された各装置を駆動制御し、マスク生成装置１が行う後述する処理を実現する。なお、本開示では、複数のＣＰＵコアを用いた処理を後述するが、複数のＣＰＵコアには物理的に複数のＣＰＵコア（区別して説明する場合には、物理コアと表記する）だけでなく、論理的に複数のＣＰＵコア（区別して説明する場合には、論理コアと表記する）も含まれる。

ＲＯＭは、不揮発性メモリであり、コンピュータのブートプログラムやＢＩＯＳ等のプログラム、データ等を恒久的に保持している。ＲＡＭ、フレームメモリは、揮発性メモリであり、記憶部１２、ＲＯＭ、記録媒体等からロードしたプログラム、データ等を一時的に保持するとともに、制御部１１が各種処理を行う為に使用するワークエリアを備える。

記憶部１２は、ＨＤＤ（Hard Disk Drive）等であり、制御部１１が実行するプログラム、プログラム実行に必要なデータ、ＯＳ（Operating System）等が格納される。プログラムに関しては、ＯＳに相当する制御プログラムや、後述する処理をコンピュータに実行させるためのアプリケーションプログラムが格納されている。これらの各プログラムコードは、制御部１１により必要に応じて読み出されてＲＡＭ、フレームメモリに移され、ＣＰＵに読み出されて各種の手段として実行される。

メディア入出力部１３（ドライブ装置）は、データの入出力を行い、例えば、ＣＤドライブ（－ＲＯＭ、－Ｒ、－ＲＷ等）、ＤＶＤドライブ（－ＲＯＭ、－Ｒ、－ＲＷ等）等のメディア入出力装置を有する。通信制御部１４は、通信制御装置、通信ポート等を有し、コンピュータとネットワーク間の通信を媒介する通信インタフェースであり、ネットワークを介して、他のコンピュータ間との通信制御を行う。ネットワークは、有線、無線を問わない。

入力部１５は、データの入力を行い、例えば、キーボード、マウス等のポインティングデバイス、テンキー等の入力装置を有する。入力部１５を介して、コンピュータに対して、操作指示、動作指示、データ入力等を行うことができる。
表示部１６は、液晶パネル等のディスプレイ装置、ディスプレイ装置と連携してコンピュータのビデオ機能を実現するための論理回路等（ビデオアダプタ等）を有する。なお、入力部１５及び表示部１６は、タッチパネルディスプレイのように、一体となっていてもよい。

周辺機器Ｉ／Ｆ（Interface）部１７は、コンピュータに周辺機器を接続させるためのポートであり、周辺機器Ｉ／Ｆ部１７を介してコンピュータは周辺機器とのデータの送受信を行う。周辺機器Ｉ／Ｆ部１７は、ＵＳＢ（Universal Serial Bus）やＩＥＥＥ１３９４やＲＳ－２３２Ｃ等によって構成されており、通常複数の周辺機器Ｉ／Ｆを有する。周辺機器との接続形態は有線、無線を問わない。バス１８は、各装置間の制御信号、データ信号等の授受を媒介する経路である。

（３．マスク生成装置１の動作）
図３のフローチャートを参照しながら、マスク生成装置１の動作について説明する。

まず、マスク生成装置１の制御部１１は、ＣＴ装置により撮影された断層画像群Ｄｏを取得する（ステップＳ１）。

続いて、制御部１１は、断層画像Ｄｏ（ｚ）毎に、断層画像Ｄｏ（ｚ）の各画素（ｘ、ｙ、ｚ）を描画対象とする重み値Ｓｖ（ｘ、ｙ、ｚ）を格納した描画重みテーブルＳｖを作成する（ステップＳ２）。描画重みテーブルＳｖを作成する処理については後述する（「４．描画重みテーブルの作成」参照）。

そして、制御部１１は、ステップＳ２において作成された描画重みテーブルＳｖを参照して、断層画像群Ｄｏの各画素（ｘ、ｙ、ｚ）の重み値Ｓｖ（ｘ、ｙ、ｚ）を所定の閾値で二値化したマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を算出し、各画素（ｘ、ｙ、ｚ）に対応するマスク値を保持する３次元のマスクデータＭａｓｋを生成する（ステップＳ３）。マスクデータＭａｓｋを生成する処理については後述する（「５．マスクデータの生成」参照）。

（４．描画重みテーブルの作成）
図４～図１６を参照しながら、図３のステップＳ２において実行される描画重みテーブルＳｖを作成する処理を説明する。
図４は、描画重みテーブルＳｖを作成する処理の流れを示すフローチャートである。

＜被写体領域算出処理＞
制御部１１は、まず、断層画像Ｄｏ（ｚ）毎に、断層画像Ｄｏ（ｚ）の各画素の信号値に基づいて、被写体領域として、矩形領域４１（第１の矩形領域）または矩形領域４２（第２の矩形領域）を算出する（ステップＳ２１）。

図５に示すように、矩形領域４１（第１の矩形領域）とは、被写体候補領域５０の最外側を囲う矩形領域である。被写体候補領域５０とは、断層画像Ｄｏ（ｚ）のうち信号値が所定の閾値（例えば、空気でない信号値以上）より大きい画素領域である。なお、被写体領域として矩形領域４１を算出する方法は、本発明者が特願２０１８－１０３７６７に既に提案しているため、ここでは、主に、被写体領域として矩形領域４２（第２の矩形領域）を算出する方法を説明する。

矩形領域４２（第２の矩形領域）とは、図５に示すように、矩形領域４１の内部に定義される矩形領域である。矩形領域４１は、信号値が所定の閾値より大きい全ての画素を含むため、図５に示すように、人体５０ａ以外の寝台５０ｂやヘッドレストなどの不要領域が含まれる場合がある。そこで、本開示では、被写体候補領域５０の信号値分布を考慮して、矩形領域４１の内部に定義される矩形領域４２を算出することで、寝台やヘッドレストなどの不要領域が含まれないように被写体領域を抽出する。

図６～図１４を参照しながら、図４のステップ２１において実行される被写体領域（矩形領域４２）を算出する処理を具体的に説明する。
図６は、被写体領域（矩形領域４２）を算出する処理の流れを示すフローチャートである。なお、断層画像Ｄｏ（ｚ）の各画素の信号値が１６ビットの場合は、あらかじめ、以下のように８ビットの信号値に変換しておく。

（式５）
Ｄ８（ｘ、ｙ、ｚ）＝｛Ｄｏ（ｘ、ｙ、ｚ）＋７００｝・２５５／１７２４
ただし、Ｄ８（ｘ、ｙ、ｚ）＜０の場合はＤ８（ｘ、ｙ、ｚ）＝０、Ｄ８（ｘ、ｙ、ｚ）＞２５５の場合はＤ８（ｘ、ｙ、ｚ）＝２５５とする。

まず、制御部１１は、被写体候補領域５０の画素の信号値で重み付けした重心座標を、矩形領域４２の中心座標として算出する（ステップＳ３１、図７（ａ））。

具体的には、制御部１１は、算出対象の矩形領域４２の中心座標を（Ｃ（ｘ）、Ｃｙ（ｚ））、重みの総和をＷｓとし、Ｃｘ（ｚ）＝Ｃｙ（ｚ）＝Ｗｓ＝０に初期化し、断層画像Ｄ８（ｘ、ｙ、ｚ）の各画素（ｘ、ｙ）（０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１）に対して、Ｄ８（ｘ、ｙ、ｚ）＞Ｖｍｉｎ（例えば、Ｖｍｉｎ＝１０）を満たす場合、被写体候補領域５０の画素であるとし、以下のようにＣｘ（ｚ）、Ｃｙ（ｚ）、Ｗｓを更新する。

（式６）
ｗ＝｛Ｄ８（ｘ、ｙ、ｚ）－Ｖｍｉｎ｝^２として、
Ｃｘ（ｚ）←Ｃｘ（ｚ）＋ｘ・ｗ
Ｃｙ（ｚ）←Ｃｙ（ｚ）＋ｙ・ｗ
Ｗｓ←Ｗｓ＋ｗ

制御部１１は、０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１の範囲のＳｘ×Ｓｙ個の画素に対して上記処理を実行したのち、Ｗｓ＞１の場合、

（式７）
Ｃｘ（ｚ）＝Ｃｘ（ｚ）／Ｗｓ
Ｃｙ（ｚ）＝Ｃｙ（ｚ）／Ｗｓ

とする。
以上のように、被写体候補領域５０の画素の信号値で重み付けした重心座標を求め、この重心座標を矩形領域４２の中心座標（Ｃｘ（ｚ）、Ｃｙ（ｚ））として算出する。

続いて、制御部１１は、算出された中心座標（Ｃｘ（ｚ）、Ｃｙ（ｚ））を起点に上下左右の４方向に矩形領域４１（第１の矩形領域）の範囲で信号値が所定の閾値より大きい画素（被写体候補領域５０の画素）の信号値で重み付けした重心（以下、「第２重心」と呼ぶ）を算出する（ステップＳ３２、図７（ｂ））。

すなわち、図８（ａ）（ｂ）に示すように、制御部１１は、中心座標（Ｃｘ（ｚ）、Ｃｙ（ｚ））（被写体候補領域５０の重心）を起点に被写体候補領域５０をＸ軸方向に２分割し、分割された各被写体候補領域５０における画素の信号値で重み付けした重心（第２重心）のＸ座標Ｇｘ１、Ｇｘ２（Ｇｘ１＜Ｇｘ２）を算出する。

同様に、図８（ｃ）（ｄ）に示すように、制御部１１は、中心座標（Ｃｘ（ｚ）、Ｃｙ（ｚ））（被写体候補領域５０の重心）を起点に被写体候補領域５０をＹ軸方向に２分割し、分割された各被写体候補領域５０における画素の信号値で重み付けした重心（第２重心）のＹ座標Ｇｙ１、Ｇｙ２（Ｇｙ１＜Ｇｙ２）を算出する。

第２重心を求める具体的な処理は以下の通りである。
制御部１１は、まず、各第２重心座標に対する重みの総和をＷｘ１、Ｗｘ２、Ｗｙ１、Ｗｙ２とし、Ｇｘ１＝Ｇｘ２＝Ｇｙ１＝Ｇｙ２＝Ｗｘ１＝Ｗｘ２＝Ｗｙ１＝Ｗｙ２＝０に初期化する。

そして、制御部１１は、断層画像Ｄ８（ｘ、ｙ、ｚ）の各画素（ｘ、ｙ）（０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１）に対して、Ｄ８（ｘ、ｙ、ｚ）＞Ｖｍｉｎ（例えば、Ｖｍｉｎ＝１０）を満たす場合、被写体候補領域５０の画素であるとし、以下のようにＧｘ１、Ｇｘ２、Ｇｙ１、Ｇｙ２、Ｗｘ１、Ｗｘ２、Ｗｙ１、Ｗｙ２を更新する。

（式８）
ｗ＝｛Ｄ８（ｘ、ｙ、ｚ）－Ｖｍｉｎ｝^２として、
ｘ＜Ｃｘ（ｚ）の場合、Ｇｘ１←Ｇｘ１＋ｘ・ｗ、Ｗｘ１←Ｗｘ１＋ｗ
ｘ≧Ｃｘ（ｚ）の場合、Ｇｘ２←Ｇｘ２＋ｘ・ｗ、Ｗｘ２←Ｗｘ２＋ｗ
ｙ＜Ｃｙ（ｚ）の場合、Ｇｙ１←Ｇｙ１＋ｙ・ｗ、Ｗｙ１←Ｗｙ１＋ｗ
ｙ≧Ｃｙ（ｚ）の場合、Ｇｙ２←Ｇｙ２＋ｙ・ｗ、Ｗｙ２←Ｗｙ２＋ｗ

制御部１１は、０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１の範囲のＳｘ×Ｓｙ個の画素に対して上記処理を実行したのち、

（式９）
Ｇｘ１＝Ｇｘ１／Ｗｘ１
Ｇｘ２＝Ｇｘ２／Ｗｘ２
Ｇｙ１＝Ｇｙ１／Ｗｙ１
Ｇｙ２＝Ｇｙ２／Ｗｙ２

とし、分割された各被写体候補領域５０における第２重心座標Ｇｘ１、Ｇｘ２（Ｇｘ１＜Ｇｘ２）、Ｇｙ１、Ｇｙ２（Ｇｙ１＜Ｇｙ２）を求める。

続いて、制御部１１は、ステップＳ３２において算出された第２重心座標を起点に更に上下左右の４方向に矩形領域４１（第１の矩形領域）の範囲で信号値が所定の閾値より大きい画素（被写体候補領域５０の画素）の信号値で重み付けした重心（以下、「第３重心」と呼ぶ）を算出する（ステップＳ３３、図７（ｃ））。

すなわち、図９（ａ）（ｂ）に示すように、制御部１１は、第２重心Ｇｘ１、Ｇｘ２を起点に被写体候補領域５０をＸ軸方向に３分割し、左右端の２つの被写体候補領域５０における画素の信号値で重み付けした重心（第３重心）のＸ座標Ｇｔｘ１、Ｇｔｘ２（Ｇｔｘ１＜Ｇｔｘ２）を算出する。

同様に、図９（ｃ）（ｄ）に示すように、制御部１１は、第２重心Ｇｙ１、Ｇｙ２を起点に被写体候補領域５０をＹ軸方向に３分割し、上下端の２つの被写体候補領域５０における画素の信号値で重み付けした重心（第３重心）のＹ座標Ｇｔｙ１、Ｇｔｙ２（Ｇｔｙ１＜Ｇｔｙ２）を算出する。

第３重心を求める具体的な処理は以下の通りである。
制御部１１は、まず、各第３重心座標に対する重みの総和をＷｘ１、Ｗｘ２、Ｗｙ１、Ｗｙ２とし、Ｇｔｘ１＝Ｇｔｘ２＝Ｇｔｙ１＝Ｇｔｙ２＝Ｗｘ１＝Ｗｘ２＝Ｗｙ１＝Ｗｙ２＝０に初期化する。

そして、制御部１１は、断層画像Ｄ８（ｘ、ｙ、ｚ）の各画素（ｘ、ｙ）（０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１）に対して、Ｄ８（ｘ、ｙ、ｚ）＞Ｖｍｉｎ（例えば、Ｖｍｉｎ＝１０）を満たす場合、被写体候補領域５０の画素であるとし、以下のようにＧｔｘ１、Ｇｔｘ２、Ｇｔｙ１、Ｇｔｙ２、Ｗｘ１、Ｗｘ２、Ｗｙ１、Ｗｙ２を更新する。

（式１０）
ｗ＝｛Ｄ８（ｘ、ｙ、ｚ）－Ｖｍｉｎ｝^２として、
ｘ＜Ｇｘ１の場合、Ｇｔｘ１←Ｇｔｘ１＋ｘ・ｗ、Ｗｘ１←Ｗｘ１＋ｗ
ｘ＞Ｇｘ２の場合、Ｇｔｘ２←Ｇｔｘ２＋ｘ・ｗ、Ｗｘ２←Ｗｘ２＋ｗ
ｙ＜Ｇｙ１の場合、Ｇｔｙ１←Ｇｔｙ１＋ｙ・ｗ、Ｗｙ１←Ｗｙ１＋ｗ
ｙ＞Ｇｙ２の場合、Ｇｔｙ２←Ｇｔｙ２＋ｙ・ｗ、Ｗｙ２←Ｗｙ２＋ｗ

制御部１１は、０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１の範囲のＳｘ×Ｓｙ個の画素に対して上記処理を実行したのち、

（式１１）
Ｇｔｘ１＝Ｇｔｘ１／Ｗｘ１
Ｇｔｘ２＝Ｇｔｘ２／Ｗｘ２
Ｇｔｙ１＝Ｇｔｙ１／Ｗｙ１
Ｇｔｙ２＝Ｇｔｙ２／Ｗｙ２

とし、分割された各被写体候補領域５０における第３重心座標Ｇｔｘ１、Ｇｔｘ２（Ｇｔｘ１＜Ｇｔｘ２）、Ｇｔｙ１、Ｇｔｙ２（Ｇｔｙ１＜Ｇｔｙ２）を求める。

そして、制御部１１は、ステップＳ３３において算出された第３重心座標に基づいて、以下のように、矩形領域４２の横方向のサイズＷ（ｚ）、縦方向のサイズＨ（ｚ）を算出する（ステップＳ３４、図７（ｄ））。

（式１２）
横方向のサイズ：Ｗ（ｚ）＝Ｇｔｘ２－Ｇｔｘ１
縦方向のサイズ：Ｈ（ｚ）＝Ｇｔｙ２－Ｇｔｙ１

また、制御部１１は、以下のように、中心座標（Ｃｘ（ｚ）、Ｃｙ（ｚ））を起点として横方向のサイズと縦方向のサイズを２種類算出する。

（式１３）
第１の横方向のサイズ：Ｗ１（ｚ）＝Ｃｘ（ｚ）－Ｇｔｘ１
第２の横方向のサイズ：Ｗ２（ｚ）＝Ｇｔｘ２－Ｃｘ（ｚ）
第１の縦方向のサイズ：Ｈ１（ｚ）＝Ｃｙ（ｚ）－Ｇｔｙ１
第２の縦方向のサイズ：Ｈ２（ｚ）＝Ｇｔｙ２－Ｃｙ（ｚ）

図１０は、算出された矩形領域４２を示す図である。図に示すように、矩形領域４２の横方向のサイズＷ（ｚ）は、右方向の第３重心Ｇｔｘ２から左方向の第３重心Ｇｔｘ１までの距離に相当し、矩形領域４２の縦方向のサイズＨ（ｚ）は、下方向の第３重心Ｇｔｙ２から上方向の第３重心Ｇｔｙ１までの距離に相当する。

また、第１の横方向のサイズＷ１（ｚ）は、左方向の第３重心Ｇｔｘ１から矩形領域４２の中心座標までの距離に相当し、第２の横方向のサイズＷ２（ｚ）は、右方向の第３重心Ｇｔｘ２から矩形領域４２の中心座標までの距離に相当する。第１の縦方向のサイズＨ１（ｚ）は、上方向の第３重心Ｇｔｙ１から矩形領域４２の中心座標までの距離に相当し、第２の縦方向のサイズＨ２（ｚ）は、下方向の第３重心Ｇｔｙ２から矩形領域４２の中心座標までの距離に相当する。

以上の処理により、被写体領域として、中心座標（Ｃｘ（ｚ）、Ｃｙ（ｚ））、横方向のサイズＷ（ｚ）（第１の横方向のサイズＷ１（ｚ）、第２の横方向のサイズＷ２（ｚ））、縦方向のサイズＨ（ｚ）（第１の縦方向のサイズＨ１（ｚ）、第２の縦方向のサイズＨ２（ｚ））により規定される矩形領域４２が算出される。

なお、ステップＳ３４において、制御部１１は、ステップＳ３２において算出された第２重心に基づいて矩形領域４２の横方向のサイズＷ（ｚ）（第１の横方向のサイズＷ１（ｚ）、第２の横方向のサイズＷ２（ｚ））と縦方向のサイズＨ（ｚ）（第１の縦方向のサイズＨ１（ｚ）、第２の縦方向のサイズＨ２（ｚ））を算出してもよい。

また、制御部１１は、ステップＳ３３において算出された第３重心を起点に更に上下左右の４方向に第４重心を算出してもよい。この場合、制御部１１は、第４重心に基づいて矩形領域４２の横方向のサイズＷ（ｚ）（第１の横方向のサイズＷ１（ｚ）、第２の横方向のサイズＷ２（ｚ））と縦方向のサイズＨ（ｚ）（第１の縦方向のサイズＨ１（ｚ）、第２の縦方向のサイズＨ２（ｚ））を算出する。

＜被写体領域補正処理＞
ところで、胸部ＣＴの場合、矩形領域４２の縦横サイズを中心点から上下左右方向の第２重心や第３重心で算出する方法をとると、骨領域内の臓器の信号値分布によっては、矩形領域４２が骨領域を除去し骨領域内の臓器等を可視化するのに適したサイズとならない場合がある。

例えば、図１１（ａ）のように領域内部の信号値が比較的大きい場合（例えば、腹部肝臓領域の場合）、算出される重心（第２重心、第３重心）が中心に寄りやすくなるため、矩形領域４２が小さめに算出される。すなわち、矩形領域４２が骨領域より内側に定義される場合がある。このような矩形領域４２に基づいてマスクデータを作成すると、本体描画対象である内臓領域が除去される場合がある。

また、図１１（ｂ）のように領域内部の信号値が比較的小さい場合（例えば、胸部肺領域の場合）、算出される重心が外側に寄りやすくなるため、矩形領域４２が大きめに算出される。すなわち、矩形領域４２が骨領域より外側に定義される場合がある。このような矩形領域４２に基づいてマスクデータを作成すると、非描画対象である骨領域が精度よく除去できない場合がある。

そこで、上記問題を解決するため、図１２に示すように、算出された矩形領域４２の横方向および縦方向のサイズに対して、被写体候補領域５０と矩形領域４２との面積比を倍率として乗算することで、矩形領域４２のサイズを補正する。

これにより、図１１（ａ）のように領域内部の信号値が比較的大きい場合（例えば、腹部肝臓領域の場合）には矩形領域４２が拡大補正され、図１１（ｂ）のように領域内部の信号値が比較的小さい場合（例えば、胸部肺領域の場合）には矩形領域４２が縮小補正され、矩形領域４２が骨領域を除去し骨領域内の臓器等を可視化するのに適したサイズに補正される。

ただし、上記補正を施すと元来不要な上胸部の矩形領域４２（被写体領域）が拡大補正されてしまうという別の問題が生じる。
図１３は、上記補正を上胸部に適用した場合を示す。図に示すように、上胸部は骨領域に囲まれた形態ではなく、高信号の骨領域が単独で分布している。そのため、矩形領域４２に対する被写体候補領域５０の面積比率が腹部並みに大きくなってしまい、腹部と同様に拡大補正されてしまう。しかしながら、上胸部は非描画対象である骨領域が大半を占めるため、上胸部領域において矩形領域４２（被写体領域）が拡大されることは望ましくない。

ところで、上胸部の骨領域は横方向に伸びているため、矩形領域４２が扁平になる特徴がある。そこで、本開示では、図１４に示すように、矩形領域４２の縦横比率（縦方向のサイズＨ（ｚ）／横方向のサイズＷ（ｚ））を算出し、この値が所定の閾値（例えば、０．６）より小さければ上胸部であると識別し、矩形領域４２の横方向および縦方向のサイズに対して、矩形領域４２と被写体候補領域５０との面積比率を乗算する代わりに、縦横比率を乗算して補正（縮小補正）する。一方、矩形領域４２の縦横比率が所定の閾値（例えば、０．６）以上であれば胸部または腹部であると判断し、矩形領域４２の横方向および縦方向のサイズに対して、矩形領域４２と被写体候補領域５０との面積比率を乗算して補正する。

このように、本開示では、矩形領域４２の縦横比率（縦方向のサイズ／横方向のサイズ）に基づいて被写体領域の部位（上胸部／胸部または腹部）を判断し、部位に応じた補正方法により矩形領域４２を補正する。

図１５は、ステップＳ２２の被写体領域の補正処理の流れを示すフローチャートである。
制御部１１は、まず、以下のように、矩形領域４２の縦横比率Ｒｖｈを算出する（ステップＳ４１）。

（式１４）
縦横比率Ｒｖｈ＝Ｈ（ｚ）／Ｗ（ｚ）

縦横比率に対する閾値をＳｒｖｈ（例えば、０．６２）とし、Ｒｖｈ＜Ｓｒｖｈの場合（ステップＳ４２；Ｎｏ）、制御部１１は、上胸部であると判断し、以下のように矩形領域４２を補正（縮小補正）する（ステップＳ４３、図１４の補正）。

（式１５）
横方向のサイズ：Ｗ（ｚ）＝Ｗ（ｚ）・Ｒｖｈ
縦方向のサイズ：Ｈ（ｚ）＝Ｈ（ｚ）・Ｒｖｈ

（式１６）
第１の横方向のサイズ：Ｗ１（ｚ）＝Ｗ１（ｚ）・Ｒｖｈ
第２の横方向のサイズ：Ｗ２（ｚ）＝Ｗ２（ｚ）・Ｒｖｈ
第１の縦方向のサイズ：Ｈ１（ｚ）＝Ｈ１（ｚ）・Ｒｖｈ
第２の縦方向のサイズ：Ｈ２（ｚ）＝Ｈ２（ｚ）・Ｒｖｈ

一方、Ｒｖｈ≧Ｓｒｖｈの場合（ステップＳ４２；Ｙｅｓ）、制御部１１は、
胸部または腹部であると判断し、まず、断層画像Ｄ８（ｘ、ｙ、ｚ）の各画素（ｘ、ｙ）（０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１）に対して、Ｄ８（ｘ、ｙ、ｚ）＞Ｖｍｉｎ（例えば、Ｖｍｉｎ＝１０）を満たす画素数をカウントし、カウント値Ｃｏｂｊを被写体候補領域５０の面積として算出する（ステップＳ４４）。

続いて、制御部１１は、矩形領域４２との面積比率Ｒｏｂｊを

（式１７）
Ｒｏｂｊ＝Ｃｏｂｊ／（Ｈ（ｚ）・Ｗ（ｚ））

により算出し、横方向、縦方向の倍率Ｒｏｘ、ＲｏｙをＲｏｘ＝Ｒｏｙ＝Ｒｏｂｊとする（ステップＳ４５）。

変倍率の上限をＭｏｘ、Ｍｏｙに設定し、Ｒｏｘ＞Ｍｏｘの場合はＲｏｘ＝Ｍｏｘ、Ｒｏｙ＞Ｍｏｙの場合はＲｏｙ＝Ｍｏｙとし、制御部１１は、以下のように矩形領域４２を補正する（ステップＳ４６、図１２の補正）。

（式１８）
横方向のサイズ：Ｗ（ｚ）＝Ｗ（ｚ）・Ｒｏｘ
縦方向のサイズ：Ｈ（ｚ）＝Ｈ（ｚ）・Ｒｏｙ

（式１９）
第１の横方向のサイズ：Ｗ１（ｚ）＝Ｗ１（ｚ）・Ｒｏｘ
第２の横方向のサイズ：Ｗ２（ｚ）＝Ｗ２（ｚ）・Ｒｏｘ
第１の縦方向のサイズ：Ｈ１（ｚ）＝Ｈ１（ｚ）・Ｒｏｙ
第２の縦方向のサイズ：Ｈ２（ｚ）＝Ｈ２（ｚ）・Ｒｏｙ

＜被写体形状算出処理＞
続いて、制御部１１は、断層画像Ｄｏ（ｚ）毎に、ステップＳ２１において算出された被写体領域（矩形領域４２）、またはステップＳ２２において補正された被写体領域（矩形領域４２）の外郭を所定の幾何形状（被写体形状）で近似し、当該幾何形状を規定するパラメータ（幾何パラメータＰ（ｚ））を算出する（ステップＳ２３）。

本開示では、以下の各領域（４象限）において異なる径が設定できる特殊な楕円（以下、「変形楕円」とも呼ぶ）で近似するものとし、幾何パラメータＰ（ｚ）として変形楕円を構成する各領域（４象限）の楕円を規定するパラメータ（楕円パラメータＰ１（ｚ）～Ｐ４（ｚ）））を算出する。
（領域１）ｘ＜Ｃｘ（ｚ）かつｙ＜Ｃｙ（ｚ）
（領域２）ｘ≧Ｃｘ（ｚ）かつｙ＜Ｃｙ（ｚ）
（領域３）ｘ＜Ｃｘ（ｚ）かつｙ≧Ｃｙ（ｚ）
（領域４）ｘ≧Ｃｘ（ｚ）かつｙ≧Ｃｙ（ｚ）

（式２０）
（領域１における楕円パラメータＰ１（ｚ））
楕円の中心座標Ｃｘ（ｚ）、Ｃｙ（ｚ）
楕円の横方向の２分割サイズｒｘ１＝Ｗ１（ｚ）・ａ１（横径）
楕円の縦方向の２分割サイズｒｙ１＝Ｈ１（ｚ）・ｂ１（縦径）
（領域２における楕円パラメータＰ２（ｚ））
楕円の中心座標Ｃｘ（ｚ）、Ｃｙ（ｚ）
楕円の横方向の２分割サイズｒｘ２＝Ｗ２（ｚ）・ａ２（横径）
楕円の縦方向の２分割サイズｒｙ１＝Ｈ１（ｚ）・ｂ１（縦径）
（領域３における楕円パラメータＰ３（ｚ））
楕円の中心座標Ｃｘ（ｚ）、Ｃｙ（ｚ）
楕円の横方向の２分割サイズｒｘ１＝Ｗ１（ｚ）・ａ１（横径）
楕円の縦方向の２分割サイズｒｙ２＝Ｈ２（ｚ）・ｂ２（縦径）
（領域４における楕円パラメータＰ４（ｚ））
楕円の中心座標Ｃｘ（ｚ）、Ｃｙ（ｚ）
楕円の横方向の２分割サイズｒｘ２＝Ｗ２（ｚ）・ａ２（横径）
楕円の縦方向の２分割サイズｒｙ２＝Ｈ２（ｚ）・ｂ２（縦径）

ａ１、ａ２、ｂ１、ｂ２は、楕円の縦横径Ｗ１（ｚ）、Ｗ２（ｚ）、Ｈ１（ｚ）、Ｈ２（ｚ）の各々に対応する補正係数（実数値で補正しない場合は１．０）であり、左方向（ａ１、Ｗ１（ｚ））、右方向（ａ２、Ｗ２（ｚ））、正面方向（ｂ１、Ｈ１（ｚ））、背面方向（ｂ２、Ｈ２（ｚ））の４方向別に設定する。

例えば、以下のように設定される。
・胸部ＣＴの場合：ａ１＝ａ２＝０．９５、ｂ１＝１．０５、ｂ２＝０．５５
・頭部ＣＴの場合：ａ１＝ａ２＝１．１５、ｂ１＝１．２、ｂ２＝１．２５

＜重み値算出処理＞
そして、制御部１１は、断層画像Ｄｏ（ｚ）毎に、ステップＳ２３において算出された変形楕円の幾何パラメータＰ（ｚ）（楕円パラメータＰ１（ｚ）～Ｐ４（ｚ））に基づいて各画素（ｘ、ｙ、ｚ）を描画対象とする重み値Ｓｖ（ｘ、ｙ、ｚ）を算出し、描画重みテーブルＳｖを作成する（ステップＳ２４）。
例えば、ｚ番目の断層画像Ｄｏ（ｚ）の各画素（ｘ、ｙ、ｚ）の重み値Ｓｖ（ｘ、ｙ、ｚ）は、以下のように算出される。

（式２１）
Ｓｖ（ｘ、ｙ、ｚ）＝Ｓｖｘｙ（ｘ、ｙ、ｚ）・Ｓｖｚ（ｚ）
（０≦Ｓｖ（ｘ、ｙ、ｚ）、Ｓｖｘｙ（ｘ、ｙ、ｚ）、Ｓｖｚ（ｚ）≦１）
Ｓｖｘｙ＝１－ｋ・ｒ（ｘ、ｙ、ｚ）
（領域１）ｘ＜Ｃｘ（ｚ）かつｙ＜Ｃｙ（ｚ）の場合：
ｒ（ｘ、ｙ、ｚ）＝｛（ｘ－Ｃｘ（ｚ））／ｒｘ１｝^２＋｛（ｙ－Ｃｙ（ｚ））／ｒｙ１｝^２
（領域２）ｘ≧Ｃｘ（ｚ）かつｙ＜Ｃｙ（ｚ）の場合：
ｒ（ｘ、ｙ、ｚ）＝｛（ｘ－Ｃｘ（ｚ））／ｒｘ２｝^２＋｛（ｙ－Ｃｙ（ｚ））／ｒｙ１｝^２
（領域３）ｘ＜Ｃｘ（ｚ）かつｙ≧Ｃｙ（ｚ）の場合：
ｒ（ｘ、ｙ、ｚ）＝｛（ｘ－Ｃｘ（ｚ））／ｒｘ１｝^２＋｛（ｙ－Ｃｙ（ｚ））／ｒｙ２｝^２
（領域４）ｘ≧Ｃｘ（ｚ）かつｙ≧Ｃｙ（ｚ）の場合：
ｒ（ｘ、ｙ、ｚ）＝｛（ｘ－Ｃｘ（ｚ））／ｒｘ２｝^２＋｛（ｙ－Ｃｙ（ｚ））／ｒｙ２｝^２
Ｓｖｚ（ｚ）＝１．０－ｋｚ・｛（ｚ－Ｓｚ／２）／ｒｚ｝^２
ｒｚ＝Ｓｚ・ｃ／２

ｋはＸＹ平面、ｋｚはＺ方向の各々減衰係数で、胸部ＣＴの場合はｋ＝１．０、ｋｚ＝０．０、頭部ＣＴの場合はｋ＝ｋｚ＝１．０に設定する。

ｃは頭部ＣＴの場合に設定されるＺ方向幅Ｓｚ／２に対応する補正係数（実数値で補正しない場合は１．０）であり、例えば、ｃ＝０．９である。

図１６は、式２１により算出される重み値Ｓｖ（ｘ、ｙ、ｚ）の２次元分布パターンを示す図である。図に示すように、変形楕円の中心（Ｃｘ（ｚ）、Ｃｙ（ｚ））の重み値Ｓｖは１（最大）となり、変形楕円内部は中心から遠ざかるにつれ重み値Ｓｖが小さくなる。具体的には、変形楕円内部は、変形楕円を構成する各楕円の横径、縦径の比率および中心座標が同一となる各同心楕円上の重み値Ｓｖが、各同心楕円の横径、縦径が大きいほど小さくなる（径の２乗に反比例して小さくなる）。変形楕円外部の重み値Ｓｖは一律に０とする。

また、頭部ＣＴの場合は、Ｓｖｚ（ｚ）によって、断層画像Ｄｏ（ｚ）のスライス位置が中央から末端に位置するにつれ（Ｓｚ／２とスライス順位ｚとの差が大きくなるにつれ）、重み値Ｓｖ（ｘ、ｙ、ｚ）が一律に減衰する（変形楕円の中心であってもＳｖ（ｘ、ｙ、ｚ）は１未満の値となる）。具体的には、断層画像のスライス順位をｚとすると、（ｚ－Ｓｚ／２）（Ｓｚ／２）の２乗に反比例して、重み値Ｓｖ（ｘ、ｙ、ｚ）が一律に減衰する

制御部１１は、全ての断層画像Ｄｏ（ｚ）について、重み値Ｓｖ（ｘ、ｙ、ｚ）（０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１、０≦ｚ≦Ｓｚ－１）を算出し、算出された重み値Ｓｖ（ｘ、ｙ、ｚ）を格納する３次元の描画重みテーブルＳｖを作成する。

以上、図４～図１６を参照しながら、描画重みテーブルＳｖを作成する処理について説明した。

（５．マスクデータの生成）
図１７のフローチャートを参照しながら、図３のステップＳ３において実行されるマスクデータＭａｓｋを生成する処理を説明する。
図１７の処理は全ての画素（ｘ、ｙ、ｚ）に対して実行される。

制御部１１は、各画素（ｘ、ｙ、ｚ）の重み値Ｓｖ（ｘ、ｙ、ｚ）を描画重みテーブルＳｖから取得し、重み値Ｓｖ（ｘ、ｙ、ｚ）・２５５＞閾値Ｔｈ（例えば、Ｔｈ＝１０）を満たす画素（ｘ、ｙ、ｚ）については（ステップＳ５１；Ｙｅｓ）、当該画素（ｘ、ｙ、ｚ）に対応するマスク値を１．０に設定し（Ｍａｓｋ（ｘ、ｙ、ｚ）＝１．０、ステップＳ５２）、重み値Ｓｖ（ｘ、ｙ、ｚ）・２５５＞１０を満たさない画素（ｘ、ｙ、ｚ）については（ステップＳ５１；Ｎｏ）、当該画素（ｘ、ｙ、ｚ）に対応するマスク値を０に設定する（Ｍａｓｋ（ｘ、ｙ、ｚ）＝０、ステップＳ５３）。以上のように、各画素（ｘ、ｙ、ｚ）の重み値Ｓｖ（ｘ、ｙ、ｚ）を所定の閾値で二値化し、各画素（ｘ、ｙ、ｚ）のマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を保持するマスクデータＭａｓｋを生成する。

図１８は、マスクデータＭａｓｋを生成する処理の別例である。
図１８の処理は全ての画素（ｘ、ｙ、ｚ）に対して実行される。
制御部１１は、まず、各画素の８ビットの信号値Ｄ８（ｘ、ｙ、ｚ）に描画重みテーブルＳｖから取得した重み値Ｓｖ（ｘ、ｙ、ｚ）を乗算して信号値を補正する（Ｄ８’ （ｘ、ｙ、ｚ）＝Ｄ８（ｘ、ｙ、ｚ）・Ｓｖ（ｘ、ｙ、ｚ）、ステップＳ６１）。

そして、制御部１１は、補正された信号値Ｄ８’（ｘ、ｙ、ｚ）がＤ８’（ｘ、ｙ、ｚ）＞閾値Ｔｈ（例えば、Ｔｈ＝１０）を満たす画素（ｘ、ｙ、ｚ）については（ステップＳ６２；Ｙｅｓ）、当該画素（ｘ、ｙ、ｚ）に対応するマスク値を１．０に設定し（Ｍａｓｋ（ｘ、ｙ、ｚ）＝１．０、ステップＳ６３）、Ｄ８’（ｘ、ｙ、ｚ）＞閾値Ｔｈを満たさない画素（ｘ、ｙ、ｚ）については（ステップＳ６２；Ｎｏ）、当該画素（ｘ、ｙ、ｚ）に対応するマスク値を０に設定する（Ｍａｓｋ（ｘ、ｙ、ｚ）＝０、ステップＳ６４）。

（６．描画重みテーブル作成処理とマスクデータ生成処理の並行化）
前述した図３のステップＳ２およびステップＳ３の処理（描画重みテーブル作成処理およびマスクデータ作成処理）は、図１９に示すように、断層画像群Ｄｏの処理領域をＺ軸方向に等分割（図の例では４分割）したうえで、各処理領域（Ｄｏ１～Ｄｏ４）ごとに並行して実行することが望ましい。これにより、高速にマスクデータＭａｓｋを生成することができる。

図２０は、描画重みテーブル作成処理およびマスクデータ作成処理を並行して実行させる仕組みの例を示す。図２０の例では、マスク生成装置１（コンピュータ）に搭載されたマザーボード（図下）は、４つのＣＰＵ（物理コア）を備え、かつ、各ＣＰＵが２つの論理コアを有する構成（マルチコアＣＰＵ）とする。すなわち、マスク生成装置１は、実質８論理コア（スレッド＃１～＃８）を備えている。はじめに、マスク生成装置１に格納されているマスク生成プログラムが起動すると、Windowsマルチタスク・オペレーティングシステムのジョブスケジューラによりＣＰＵコア＃１のスレッド＃１に割り当てられて実行される。続いて、マスク生成プログラムは、処理領域Ｄｏ１に係る描画重みテーブル作成処理およびマスクデータ作成処理（図の並行処理スレッド＃１）、処理領域Ｄｏ２に係る描画重みテーブル作成処理およびマスクデータ作成処理（図の並行処理スレッド＃２）、処理領域Ｄｏ３に係る描画重みテーブル作成処理およびマスクデータ作成処理（図の並行処理スレッド＃３）、処理領域Ｄｏ４に係る描画重みテーブル作成処理およびマスクデータ作成処理（図の並行処理スレッド＃４）の各スレッドを起動し、Windowsマルチタスク・オペレーティングシステムのジョブスケジューラに登録する。

Windowsマルチタスク・オペレーティングシステムでは、コアの空き状況等によって、登録された各並行処理スレッド＃１～＃４を、８論理コア（スレッド＃１～＃８）のいずれかに割り当て、各並行処理スレッド＃１～＃４を並行して実行する。図２０の例では、並行処理スレッド＃１、並行処理スレッド＃２および並行処理スレッド＃３は、ジョブスケジューラにより物理コアが空いているＣＰＵコア＃２のスレッド＃３、ＣＰＵコア＃３のスレッド＃５、ＣＰＵコア＃４のスレッド＃７に各々割り当てられる。残る並行処理スレッド＃４は、物理コアは空いていないがＣＰＵコア＃１の論理スレッド＃２に割り当てられている。ＣＰＵコア＃１の論理スレッド＃１ではレンダリングプログラム本体が使用中であるが、マスク生成プログラム本体は４つの並行処理スレッドの起動を行って、これらの終了待ち状態になっているため、ＣＰＵコア＃１の論理スレッド＃２に割り当てられている並行処理スレッド＃４が主に実行することになる。

以上、マスク生成装置１の動作について説明した。マスク生成装置１は、ＣＴ装置により撮影された断層画像を取得し（図３のステップＳ１）、断層画像Ｄｏ（ｚ）毎に各画素（ｘ、ｙ）の信号値に基づいて被写体領域（矩形領域４２）を算出し（図４のステップＳ２１、（式７）（式１２）（式１３）参照）、算出された被写体領域（矩形領域４２）の外郭を所定の幾何形状（変形楕円）で近似し、当該幾何形状（変形楕円）を規定する幾何パラメータＰ（ｚ）（４象限の楕円パラメータＰ１（ｚ）～Ｐ４（ｚ））を算出する（図４のステップＳ２３、（式２０）参照）。続いて、断層画像Ｄｏ（ｚ）毎に、幾何パラメータＰ（ｚ）に基づいて断層画像の各画素（ｘ、ｙ）を描画対象とする重み値Ｓｖ（ｘ、ｙ、ｚ）を算出する（図４のステップＳ２４、（式２１）参照）。そして、各画素（ｘ、ｙ、ｚ）の重み値Ｓｖ（ｘ、ｙ、ｚ）を所定の閾値で二値化したマスク値を算出し、各画素（ｘ、ｙ、ｚ）のマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を保持するマスクデータＭａｓｋを生成する（図３のステップＳ３）。
これにより、被写体領域の外郭を近似した幾何形状に基づいて所望の領域を可視化するためのマスクデータを容易に作成することが可能となる。

また、被写体候補領域５０の信号値分布を考慮して被写体領域（矩形領域４２）を算出することで（図６参照）、人体以外の寝台やヘッドレストなどの不要領域を不可視とするマスクデータを生成することができる。

また、胸部ＣＴにおいては、被写体領域（矩形領域４２）の縦横比率に基づいて部位（上胸部／胸部または腹部）を判断し、部位に応じた補正方法により被写体領域（矩形領域４２）を補正することで（図１５参照）、被写体領域を高精度に算出することができる。
さらに、被写体領域（矩形領域４２）の外郭を、胸骨側と椎骨側で楕円の径を変えることが可能な変形楕円で近似することで（図１６参照）、肋骨、胸骨を含む骨領域が的確に除去され、肺、心臓、血管、肝臓、腎臓などの診断上重要な内臓領域等が良好に可視化されるマスクデータを生成することができる。

なお、本開示では各画素（ｘ、ｙ、ｚ）を描画対象とする重み値Ｓｖ（ｘ、ｙ、ｚ）を格納する描画重みテーブルＳｖを作成することとしたが、描画重みテーブルＳｖの作成は必須ではない。すなわち、図４のステップＳ２４（（式２１）参照）で算出された重み値Ｓｖ（ｘ、ｙ、ｚ）に対して都度、図１７または図１８の処理によりマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を計算するようにしてもよい。この場合、描画重みテーブルＳｖのデータ領域をメモリ上に確保する必要がない。

更に、（式２１）で示される重み値Ｓｖ（ｘ、ｙ、ｚ）の算出自体も必須ではない。図１７に示す実施例の場合、Ｓｖ（ｘ、ｙ、ｚ）・２５５＞Ｔｈを満たす場合、マスク値Ｍａｓｋ（ｘ、ｙ、ｚ）＝１となる。胸部ＣＴにおいては、Ｓｖｘｙ（ｘ、ｙ、ｚ）＝１－ｋ・ｒ（ｘ、ｙ、ｚ）＞Ｔｈ／２５５を満たす場合、マスク値Ｍａｓｋ（ｘ、ｙ、ｚ）＝１となる。これは、ボクセル座標（ｘ、ｙ、ｚ）が図１６に示される複数の同心の変形楕円のいずれかの単一の変形楕円の内部に含まれる場合、マスク値Ｍａｓｋ（ｘ、ｙ、ｚ）＝１となると解釈できる。即ち、胸部ＣＴにおいては、マスク値Ｍａｓｋ（ｘ、ｙ、ｚ）は図１６に示されるいずれかの単一の変形楕円の内部（マスク値：１）か外部（マスク値：０）かを定義したものであり、ボクセル座標（ｘ、ｙ、ｚ）の単一の変形楕円に対する内外判定処理によりマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を算出することが可能である。頭部ＣＴにおいても、Ｓｖ（ｘ、ｙ、ｚ）＝Ｓｖｘｙ（ｘ、ｙ、ｚ）・Ｓｖｚ（ｚ）で定義される単一の変形楕円体に対する内外判定処理によりマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を算出することは可能である。ただし、特に頭部ＣＴの場合、内外判定処理によりマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を算出する方法より、（式２１）に基づいて描画重みテーブルＳｖを作成し、Ｔｈで二値化してマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を算出する前述の方法の方が計算負荷は少ない。

また、本開示では、図４のステップＳ２３において、被写体領域（矩形領域４２）の外郭を変形楕円（図１６参照）で近似する場合を説明したが、近似する幾何形状は変形楕円に限定されない。例えば、本発明者が特願２０１８－１０３７６７において既に提案しているように、被写体領域（矩形領域４２）の外郭を単一の楕円で近似してもよいし、特願２０１８－１２４１２５において既に提案しているように、被写体領域（矩形領域４２）の外郭を複数の楕円で近似してもよい。複数の楕円で近似する場合も、ボクセル座標（ｘ、ｙ、ｚ）の複数の楕円に対する内外判定処理を行い、判定結果のＡＮＤまたはＯＲ演算を行うことによりマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を算出することが可能である。

また、図４のステップＳ２１において、被写体領域として矩形領域４２（第２の矩形領域、図５参照）を算出したが、矩形領域４１（第１の矩形領域、図５参照）を算出してもよい。この場合、図４のステップＳ２２における被写体領域の補正は行わず、図４のステップＳ２３では、例えば、被写体領域の外郭を単一の楕円（既提案の特願２０１８－１０３７６７参照）や複数の楕円（既提案の特願２０１８－１２４１２５参照）で近似する。

［マスクデータＭａｓｋの適用例］
以降、生成されたマスクデータＭａｓｋの適用例として、マスクデータＭａｓｋを用いて３次元再構成像を生成する処理を説明する。

（７．３次元再構成像生成処理の概要）
図２１は、３次元再構成像を生成する３次元再構成像生成装置２が実行する再構成像生成処理の概要を示す図である。図に示すように、３次元再構成像生成装置２は、複数の断層画像（図の例は、５１２×５１２ピクセルの３７０枚の胸部ＣＴ画像）に基づいてマスクデータＭａｓｋを用いた再構成像生成処理を実行し、所定の視点から被検体を観察したボリュームレンダリング画像、ＭＩＰ（Maximum Intensity Projection）画像、ＭＰＲ（Multi-Planar
Reconstruction）画像などの３次元再構成像を生成する。

（８．３次元再構成像生成装置２の機能構成）
図２２は、３次元再構成像生成装置２の機能構成を示す図である。図に示すように、３次元再構成像生成装置２は、画像取得部２１、描画重みテーブル作成部２２、マスクデータ生成部２３、クリッピング領域設定部２４、有効ボクセル領域算出部２５、３次元再構成像生成部２６、３次元再構成像表示部２７を備える。画像取得部２１、描画重みテーブル作成部２２（被写体領域算出部３１、被写体領域補正部３２、被写体形状算出部３３、重み値算出部３４）、マスクデータ生成部２３は、図１のマスク生成装置１の機能構成と同様であるため、説明を省略する。

クリッピング領域設定部２４は、被写体のクリッピング領域ＲＯＩ（関心領域）を設定する。本開示では、クリッピングＲＯＩは直方体とし、以下のように、Ｘ方向ＲＯＩ（Ｘ方向の開区間）、Ｙ方向ＲＯＩ（Ｙ方向の開区間）、Ｚ方向ＲＯＩ（Ｚ方向の開区間）によって定義する。

（式２２）
Ｘ軸方向ＲＯＩ：Ｘｓ－Ｘｅ
Ｙ軸方向ＲＯＩ：Ｙｓ－Ｙｅ
Ｚ軸方向ＲＯＩ：Ｚｓ－Ｚｅ
尚、クリッピング領域を設定しない場合は、Ｘｓ＝０、Ｘｅ＝Ｓｘ－１、Ｙｓ＝０、Ｙｅ＝Ｓｙ－１、Ｚｓ＝０、Ｚｅ＝Ｓｚ－１となる。

有効ボクセル領域算出部２５は、ボクセル値が描画対象のボクセルを規定するように予め定められた条件を満たすボクセル（以下、「有効ボクセル」と表記する場合がある）を全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。有効ボクセル領域Ｖｒは、以下のように、Ｘ軸方向ＲＯＩ、Ｙ軸方向ＲＯＩ、Ｚ軸方向ＲＯＩによって定義される。ここで「ボクセル値」とは、各ボクセルに対応づけられる数値情報であり、モダリティで測定された信号値（ＣＴ値）や、信号値をカラーマップＣｍａｐに通して得られる不透明度や色値、などを含む。

（式２３）
Ｘ軸方向ＲＯＩ：Ｘｉｓ－Ｘｉｅ
Ｙ軸方向ＲＯＩ：Ｙｉｓ－Ｙｉｅ
Ｚ軸方向ＲＯＩ：Ｚｉｓ－Ｚｉｅ

有効ボクセル領域Ｖｒを算出することによって、ボリュームレンダリング画像やＭＩＰ画像を生成する処理（レイキャスティング処理）が、有効ボクセル領域Ｖｒ内のボクセルのみを参照して実行されるため、計算負荷を大幅に抑えることができ、境界が直方体のためレイキャスティング処理における起点座標を直線（仮想光線）との交点に基づいて迅速に探索できる。

３次元再構成像生成部２６は、マスクデータＭａｓｋを参照しながら、断層画像群Ｄｏから３次元再構成像を生成する。３次元再構成像生成部２６は、図２２に示すように、ボリュームレンダリング画像生成部５１、ＭＩＰ画像生成部５２、ＭＰＲ画像生成部５３から構成される。

ボリュームレンダリング画像生成部５１は、後述するレイキャスティング処理において、マスクデータＭａｓｋに基づいて断層画像群Ｄｏの各ボクセルの不透明度に対するマスク値（透明／不透明）を設定したうえで（後述する図３３のステップＳ１４４）、カラーマップ設定部５１０により設定されたカラーマップＣｍａｐを断層画像群Ｄｏの各ボクセルに適用し（後述する図３３のステップＳ１４５）、ボリュームレンダリング画像ＩｍａｇｅＶＲを生成する。

カラーマップＣｍａｐは、信号値ｖと色値（具体的にはＲＧＢ値）及び不透明度(α値)との対応関係を定義したものであり、信号値ｖを２４ビットの色（ＲＧＢ値）及び８ビットの不透明度（α値）に変換する関数（実体的には２次元のデータテーブル）として表現される。例えば、１６ビットの断層画像群Ｄｏ（（式１）参照）に適用されるカラーマップＣｍａｐは次のように定義される。

（式２４）
０≦Ｃｍａｐ（ｖ、ｎ）≦２５５
－３２７６８≦ｖ≦３２７６７
ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）、３（α）

また、８ビットの断層画像群Ｄｏ（（式２）参照）に適用されるカラーマップＣｍａｐ（ｖ、ｎ）は次のように定義される。

（式２５）
０≦Ｃｍａｐ（ｖ、ｎ）≦２５５
０≦ｖ≦２５５
ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）、３（α）

（式２４）、（式２５）に示すカラーマップＣｍａｐ（ｖ、ｎ）（０≦ｎ≦３）のうち、特にカラーマップＣｍａｐ（ｖ、ｎ）（０≦ｎ≦２）は、信号値ｖを色値（ＲＧＢ値）に変換する関数に相当する。

また、（式２４）、（式２５）に示すカラーマップＣｍａｐ（ｖ、ｎ）（０≦ｎ≦３）のうち、特にカラーマップＣｍａｐ（ｖ、３）は、信号値ｖを不透明度に変換する関数に相当し、一般的に「オパシティカーブ」と呼ばれる。

カラーマップＣｍａｐ（ｖ、３）（オパシティカーブ）は、例えば、通常組織（１６ビットの信号値ｖ＝０～１２０程度）の不透明度が２５５に設定（不透明度＝２５５は不透明（光がすべて反射）であることを示す）され、骨領域（１６ビットの信号値ｖ＝１０００前後）の不透明度が１～２５４の中間値に設定（不透明度＝１～２５４は半透明（入射光の一部が反射され、その他は透過）であることを示す）され、その他の組織等の部位を不透明度が０に設定（不透明度＝０は透明（入射光の全てが透過）であることを示す）される。

生成されるボリュームレンダリング画像ＩｍａｇｅＶＲは、フルカラー（２４ビット）の３次元再構成像（Ｓｉｚｅ×Ｓｉｚｅの画像）であり、以下のように定義される。

（式２６）
０≦ＩｍａｇｅＶＲ（ｘ、ｙ、ｎ）≦２５５
０≦ｘ≦Ｓｉｚｅ－１、０≦ｙ≦Ｓｉｚｅ－１
ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）

ＭＩＰ画像生成部５２は、後述するレイキャスティング処理において、マスクデータＭａｓｋに基づいて断層画像群Ｄｏの各ボクセルの信号値に対するマスク値を設定したうえで（後述する図４０のステップＳ２２４）、断層画像群Ｄｏの各ボクセルに信号値を取得し（後述する図４０のステップＳ２２５）、ＭＩＰ画像ＩｍａｇｅＭＩＰを生成する。

生成されるＭＩＰ画像ＩｍａｇｅＭＩＰは、モノクロ（１６ビットまたは８ビット）の３次元再構成像（Ｓｉｚｅ×Ｓｉｚｅの画像）であり、以下のように定義される。ＭＩＰ画像ＩｍａｇｅＭＩＰを生成する段階では、階調は元のＤＩＣＯＭ画像の信号値のままモノクロ（１６ビット）で計算する方法が一般的であり、ＩｍａｇｅＭＩＰ（ｘ、ｙ）は－３２７６８≦ＩｍａｇｅＭＩＰ（ｘ、ｙ）≦３２７６７の値をもつ場合があるが、表示される段階では以下のようにモノクロ（８ビット）に変換される。

（式２７）
０≦ＩｍａｇｅＭＩＰ（ｘ、ｙ）≦２５５
０≦ｘ≦Ｓｉｚｅ－１、０≦ｙ≦Ｓｉｚｅ－１

ＭＰＲ画像生成部５３は、断層画像群Ｄｏに基づいて任意断面のＭＰＲ画像ＩｍａｇｅＭＰＲを生成するともに、ＭＰＲ画像ＩｍａｇｅＭＰＲにマスクデータＭａｓｋに基づいて生成されたマスク画像を合成する。

３次元再構成像表示部２７は、３次元再構成像生成部２６（ボリュームレンダリング画像生成部５１、ＭＩＰ画像生成部５２、ＭＰＲ画像生成部５３）により生成された３次元再構成像を表示画面（表示部１６）に表示する。

（９．３次元再構成像生成装置２のハードウェア構成）
図２３は、３次元再構成像生成装置２のハードウェア構成を示す図である。図に示すように、３次元再構成像生成装置２は、制御部１１、記憶部１２、メディア入出力部１３、通信制御部１４、入力部１５、表示部１６、周辺機器Ｉ／Ｆ部１７等が、バス１８を介して接続される汎用のコンピュータで実現される。各構成は図２と同様であるため、記載を省略する。

（１０．３次元再構成像生成装置２の動作）
図２４のフローチャートを参照しながら、３次元再構成像生成装置２の全体の動作について説明する。

まず、３次元再構成像生成装置２の制御部１１は、ＣＴ装置により撮影された１６ビットの断層画像群Ｄｏ（（式１）参照）を取得する（ステップ７１）。この際、制御部１１は、必要に応じて、１６ビットの断層画像群Ｄｏを８ビットの断層画像群Ｄｏ（（式２）参照）に階調圧縮してもよい。

（１）断層画像群ＤｏのＳｚ／２番目の中間スライスＤｏ（ｚ／２）における全ての画素の最小値Ｄｍｉｎ、最大値Ｄｍａｘを算出する。なおＤｍｉｎとＤｍａｘは全てのスライスより最小値と最大値を求める方法をとっても良いが、本実施の形態のように中間スライスだけで決定する方法をとることで、大容量の１６ビットの断層画像群Ｄｏをメモリ上に保持する必要がなくなる。

（２）下限値Ｌｍｉｎ＝（Ｄｍａｘ－Ｄｍｉｎ）・γ＋Ｄｍｉｎ、上限値Ｌｍａｘ＝（Ｄｍａｘ－Ｄｍｉｎ）・（１－γ）＋Ｄｍｉｎを設定する。ここで、γは階調圧縮画像のコントラスト調整幅で、０に近いほどコントラストが増大する（但し、輝度が小さくなる）。通常はγ＝０．１に設定する。

（３）断層画像群Ｄｏを以下の計算式で８ビットの断層画像群Ｄｏ（（式２）参照）に変換する（下限値Ｌｍｉｎ～上限値Ｌｍａｘの範囲で２５６段階に圧縮する）。

（式２８）
Ｄｏ（ｘ、ｙ、ｚ）
＝（Ｄｏ（ｘ、ｙ、ｚ）－Ｌｍｉｎ）・２５５／（Ｌｍａｘ－Ｌｍｉｎ）
但し、Ｄｏ（ｘ、ｙ、ｚ）＞２５５の場合はＤｏ（ｘ、ｙ、ｚ）＝２５５、Ｄｏ（ｘ、ｙ、ｚ）＜０の場合はＤｏ（ｘ、ｙ、ｚ）＝０に飽和させる。

このように階調圧縮をすることで、断層画像群を保持するためのメモリ容量を半分に抑えることができる。たとえ信号値の階調が１６ビットあっても、カラーマップＣｍａｐにより、変換される色値（ＲＧＢ値）及び不透明度（α値）の階調はディスプレイの階調に合わせて８ビットに制限されるため、階調圧縮に伴う画質劣化は殆ど生じない。

続いて、制御部１１は、断層画像Ｄｏ（ｚ）毎に、断層画像Ｄｏ（ｚ）の各画素（ｘ、ｙ、ｚ）を描画対象とする重み値Ｓｖ（ｘ、ｙ、ｚ）を格納した描画重みテーブルＳｖを作成し（ステップＳ７２）、作成された描画重みテーブルＳｖを参照して、断層画像群Ｄｏの各画素（ｘ、ｙ、ｚ）の重み値Ｓｖ（ｘ、ｙ、ｚ）を所定の閾値で二値化したマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を算出し、各画素（ｘ、ｙ、ｚ）に対応するマスク値Ｍａｓｋ（ｘ、ｙ、ｚ）を保持する３次元のマスクデータＭａｓｋを生成する（ステップＳ７３）。ステップＳ７２、ステップＳ７３の処理は、図３のステップＳ２、ステップＳ３と同様である。

そして、制御部１１は、ステップＳ７３において生成されたマスクデータＭａｓｋを参照しながら、断層画像群Ｄｏから３次元再構成像を生成し（ステップＳ７４）、生成した３次元再構成像を表示する（ステップＳ７５）。

以降、ステップＳ７４において、ボリュームレンダリング画像、ＭＩＰ画像、およびＭＰＲ画像を生成する例を順に説明する。

（１１．ボリュームレンダリング画像の生成）
図２５～図３３を参照しながら、ボリュームレンダリング画像を生成する処理を説明する。

図２５は、ボリュームレンダリング画像を生成する処理の全体の流れを示すフローチャートである。
制御部１１は、ユーザからカラーマップＣｍａｐの設定を受け付ける（ステップＳ８１）。ここで、断層画像群Ｄｏの信号値が１６ビットの場合（（式１）参照）は、１６ビット対応のカラーマップＣｍａｐ（（式２４）参照）を設定し、断層画像群Ｄｏの信号値が８ビットに圧縮されている場合（（式２）参照）は、８ビット対応のカラーマップＣｍａｐ（（式２５）参照）を設定する。

カラーマップＣｍａｐの設定方法は特に限定されず、例えば、制御部１１は、予め用意されている複数のカラーマップＣｍａｐの中から、断層画像群Ｄｏに適用するカラーマップＣｍａｐをユーザに選択させることで設定することができる。また、制御部１１は、前回レンダリング時に使用したカラーマップＣｍａｐを記憶部１２から読込み、断層画像群Ｄｏに適用するカラーマップＣｍａｐとして設定しても良い。

また、制御部１１は、カラーマップ設定画面（不図示）においてユーザが作成したカラーマップＣｍａｐを設定しても良い。この場合、制御部１１は、カラーマップ設定画面において、代表的な信号値ｖ（例えば、１０箇所程度の特徴的な信号値）と当該信号値ｖに対する色値（ＲＧＢ値）及び不透明度（α値）をユーザに設定させ、ユーザが設定した代表的な信号値ｖと色値及び不透明度との対応関係に基づいて所定範囲の信号値ｖを色値及び不透明度に変換するカラーマップＣｍａｐを自動生成する。

続いて、制御部１１は、必要に応じて、ユーザからクリッピング領域ＲＯＩ（（式２２）参照）の設定を受け付ける（ステップＳ８２）。

続いて、制御部１１は、有効ボクセル領域Ｖｒを算出する（ステップＳ８３）。具体的には、図２６に示すように、制御部１１は、不透明度が０でない値（α＞０）をもつボクセル（不透明ボクセル）を有効ボクセルとし、これらの有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。クリッピング領域ＲＯＩやマスクデータＭａｓｋが設定されている場合は、制御部１１は、クリッピングＲＯＩやマスクデータＭａｓｋにより規定された描画範囲内において、不透明度が０でない値（α＞０）をもつボクセル（不透明ボクセル）を有効ボクセルとし、これらの有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。

そして、制御部１１は、ステップＳ７１（図２４）において取得した断層画像群Ｄｏ（３次元ボクセル）、ステップＳ７３（図２４）において生成されたマスクデータＭａｓｋ、ステップＳ８１（図２５）において設定されたカラーマップＣｍａｐ、ステップＳ８２（図２５）において設定されたクリッピング領域ＲＯＩ、ステップＳ８３（図２５）において算出された有効ボクセル領域Ｖｒを参照して、ボリュームレンダリング画像ＩｍａｇｅＶＲを生成する（ステップＳ８３）。

特に本開示では、制御部１１は、レイキャスティング処理において、マスクデータＭａｓｋに基づいて断層画像群Ｄｏの各ボクセルの不透明度に対するマスク値を設定したうえで（後述する図３３のステップＳ１４４）、断層画像群Ｄｏの各ボクセルにカラーマップＣｍａｐ（オパシティカーブ）を適用し（後述する図３３のステップＳ１４５）、ボリュームレンダリング画像ＩｍａｇｅＶＲを生成する。これにより、被写体形状（幾何形状）に基づいて生成したマスクデータＭａｓｋの描画・非描画情報が反映されたボリュームレンダリング画像が生成される。

図２７は、ステップＳ８４において実行される、ボリュームレンダリング画像ＩｍａｇｅＶＲを生成する具体的な処理の流れを示す。
まず、制御部１１は、座標変換パラメータを設定する（ステップＳ９１）。本開示では、後述するレイキャスティング処理（ステップＳ９３、図２８）で、視点座標系の各３次元座標毎にボクセルを参照するため、逐次座標変換を行う方法をとる。そのため、制御部１１は、以下のように、各座標変換処理に共通する座標変換パラメータを事前に設定・算出しておく。

＜座標変換パラメータ＞
回転パラメータ行列Ｒ：
Ｒ＝［Ｒ１１ Ｒ１２ Ｒ１３；
Ｒ２１ Ｒ２２ Ｒ２３；
Ｒ３１ Ｒ３２ Ｒ３３］
（ボクセル座標系から視点座標系への座標変換を行うための３×３の行列の逆行列、ＧＵＩ側はボクセル座標系から視点座標系への座標変換を指示するが、レンダリング側は視点座標系からボクセル座標系に座標変換を行う。）
ＸＹＺ軸方向のオフセット：
Ｘｏｆｆ、Ｙｏｆｆ、Ｚｏｆｆ（２次元画面上で指定するため、通常Ｚｏｆｆ＝０）
クリッピング領域：
Ｘ軸方向ＲＯＩ：Ｘｓ－Ｘｅ
Ｙ軸方向ＲＯＩ：Ｙｓ－Ｙｅ
Ｚ軸方向ＲＯＩ：Ｚｓ－Ｚｅ
有効ボクセル領域（外接直方体領域）：
Ｘ軸方向ＲＯＩ：Ｘｉｓ－Ｘｉｅ
Ｙ軸方向ＲＯＩ：Ｙｉｓ－Ｙｉｅ
Ｚ軸方向ＲＯＩ：Ｚｉｓ―Ｚｉｅ
拡大縮小倍率Ｓｃａｌｅ（ＸＹＺ軸方向で同一）
Ｘ軸方向変倍率：Ｓｃｘ、Ｙ軸方向変倍率：Ｓｃｙ、Ｚ軸方向変倍率Ｓｃｚ
座標変換サブサンプル・オフセット：Ｘ軸方向ｄｘ、Ｙ軸方向ｄｙ、Ｚ軸方向ｄｚ
生成するボリュームレンダリング画像のサイズ：Ｓｉｚｅ（ＸＹ軸方向で同一）
仮想光線のサブサンプリング倍率：Ｓｒａｙ（通常は１で、値が１より大きいと粗くなり、１未満だと高精細になる。ボリュームレンダリング画像を生成する場合はＳｒａｙ＝１、ＭＩＰ画像を生成する場合はＳｒａｙ＝２などに設定する。）

制御部１１は、上記した回転パラメータ行列Ｒに対して、初期値として単位行列を設定する。すなわち、Ｒ１１＝１、Ｒ１２＝０、Ｒ１３＝０、Ｒ２１＝０、Ｒ２２＝１、Ｒ２３＝０、Ｒ３１＝０、Ｒ３２＝０、Ｒ３３＝１と設定する。
そして、ＧＵＩの指示に従い、Ｘ軸中心回転Ｒｘ、Ｙ軸中心回転Ｒｙ、Ｚ軸中心回転Ｒｚ（角度単位：ラジアン）のいずれかを逐次指定し、以下のように、各々回転行列Ａを生成して回転パラメータ行列Ｒに右から乗算して、回転パラメータ行列Ｒを更新する。これにより、ＧＵＩの指示により生成されるボクセル座標系から視点座標系への回転行列の逆行列が算出される。

回転行列Ａを
Ａ＝［Ａ１１ Ａ１２ Ａ１３；
Ａ２１ Ａ２２ Ａ２３；
Ａ３１ Ａ３２ Ａ３３］
とすると、
Ｘ軸中心回転Ｒｘの場合の回転行列Ａの各要素は、
Ａ１１＝１、Ａ１２＝０、Ａ１３＝０
Ａ２１＝０、Ａ２２＝ｃｏｓＲｘ、Ａ２３＝ｓｉｎＲｘ
Ａ３１＝０、Ａ３２＝ｓｉｎＲｘ、Ａ３３＝ｃｏｓＲｘ
Ｙ軸中心回転Ｒｙの場合の回転行列Ａの各要素は、
Ａ１１＝ｃｏｓＲｙ、Ａ１２＝０、Ａ１３＝ｓｉｎＲｙ
Ａ２１＝０、Ａ２２＝１、Ａ２３＝０
Ａ３１＝－ｓｉｎＲｙ、Ａ３２＝０、Ａ３３＝ｃｏｓＲｙ
Ｚ軸中心回転Ｒｚの場合の回転行列Ａの各要素は、
Ａ１１＝ｃｏｓＲｚ、Ａ１２＝ｓｉｎＲｚ、Ａ１３＝０
Ａ２１＝－ｓｉｎＲｚ、Ａ２２＝ｃｏｓＲｚ、Ａ２３＝０
Ａ３１＝０、Ａ３２＝０、Ａ３３＝１
となる。
回転パラメータ行列Ｒは、Ｒ←Ｒ×Ａと更新される。

以上で定義された座標変換パラメータに基づく座標変換処理は、レイキャスティング処理（ステップＳ９３、図２８）の各処理の中で逐次実行される。

続いて、制御部１１は、処理回数ｌ＝０、サブサンプル・オフセット値をｄｘ＝ｄｙ＝ｄｚ＝０に初期化する（ステップＳ９２）。
そして、制御部１１は、レイキャスティング処理（各座標（ｘ、ｙ）毎に色値を算出する処理）を実行する（ステップＳ９３）。

＜レイキャスティング処理＞
図２８のフローチャートを参照して、レイキャスティング処理について説明する。
制御部１１は、まず、生成する２４ビット（ＲＧＢ）のボリュームレンダリング画像ＩｍａｇｅＶＲ（ｘ、ｙ、ｎ）の初期値を全て０に設定する（ＩｍａｇｅＶＲ（ｘ、ｙ、ｎ）＝０、ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ））。そして、サブサンプル回数Ｌとして、各２次元座標（ｘ、ｙ）（０≦ｘ≦Ｓｉｚｅ－１、０≦ｙ≦Ｓｉｚｅ－１）に対して、以下の処理を実行する。

制御部１１は、背景色ｂｇ（ｎ）（ｎ＝０（Ｒ）、１（Ｇ）、２（Ｂ）、０≦ｂｇ（ｎ）≦２５５）を設定し、ｘ＝０、ｙ＝０とする（ステップＳ１０１）。

続いて、制御部１１は、仮想光強度Ｔｒａｎｓ＝１．０、累積輝度値Ｅｎｅｒｇｙ（ｎ）＝０．０（０≦ｎ≦２）に初期化する（ステップＳ１０２）。

続いて、制御部１１は、有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する（ステップＳ１０３）。有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する処理については後述する（「１１－１．交点算出処理」参照）。

続いて、制御部１１は、ＺｓｔをＺｓｔ＝Ｚｃに設定し（ステップＳ１０４）、Ｚｓｔから先頭の有効ボクセル（不透明ボクセル）のＺ座標ｚ（仮想光線の照射を開始する起点座標）を探索する（ステップＳ１０５）。すなわち、ステップＳ１０３において算出された交点から視線方向に向かってレイキャスティング計算を開始する起点座標を探索する。起点座標ｚを探索する処理については後述する（「１１－２．起点座標探索処理」参照）。

ｚ＜０の場合、ステップＳ１１１へ移行する。

一方、ｚ≧０の場合、制御部１１は、３次元座標（ｘ、ｙ、ｚ）を座標変換してボクセルα値（Ｖα）を取得し、Ａｌｐｈａ＝Ｖα／２５５を求める（ステップＳ１０６）。ボクセルα値（Ｖα）を取得する処理については後述する（「１１－４．ボクセルα値、ＲＧＢ値の取得（補間あり）」参照）。

一方、ステップＳ１０６において算出したＡｌｐｈａがＡｌｐｈａ＜０（有効ボクセル領域Ｖｒを通過）、またはＡｌｐｈａ＝０かつｚ＝０の場合、ステップＳ１１１へ移行し、画素（ｘ、ｙ）におけるＲＧＢ画素値を決定する。これにより、有効ボクセル領域Ｖｒを逸脱（通過）している場合（Ａｌｐｈａ＜０）は、その先に描画対象のボクセルは存在しないため、レイキャスティング処理を早期に打ち切り、冗長な処理を省略する。

一方、ステップＳ１０６において算出したＡｌｐｈａがＡｌｐｈａ＞０の場合、制御部１１は、３次元座標（ｘ、ｙ、ｚ）を座標変換してボクセルＲＧＢ値Ｖｃ（ｎ）を取得する（ステップＳ１０８）。ボクセルＲＧＢ値（Ｖｃ（ｎ））を取得する処理については後述する（「１１－４．ボクセルα値、ＲＧＢ値の取得（補間あり）」参照）。

続いて、制御部１１は、仮想光線サブサンプルに基づくα値を、
Ａｌｐｈａ’＝１－（１－Ａｌｐｈａ）^{１／Ｓｒａｙ}
と補正し、
累積輝度を、
Ｅｎｅｒｇｙ（ｎ）＝Ｅｎｅｒｇｙ（ｎ）＋Ｔｒａｎｓ／Ａｌｐｈａ・Ｖｃ（ｎ）／２５５
透過光強度を、
Ｔｒａｎｓ＝Ｔｒａｎｓ・（１．０－Ａｌｐｈａ）
と更新する（ステップＳ１０９）。

Ｔｒａｎｓ＜０．００１の場合、ステップＳ１１１へ移行し、画素（ｘ、ｙ）におけるＲＧＢ画素値を決定する。Ｔｒａｎｓ≧０．００１の場合、制御部１１は、ｚ＝ｚ－１に更新し（ステップＳ１１０）、ｚ≧０の場合、ステップＳ１０６に戻り、ｚ＜０の場合、ステップＳ１１１へ移行し、画素（ｘ、ｙ）におけるＲＧＢ画素値を決定する。

ステップＳ１１１において、制御部１１は、画素（ｘ、ｙ）におけるＲＧＢ画素値を次のように決定する。

（式２９）
ＩｍａｇｅＶＲ（ｘ、ｙ、ｎ）＝ＩｍａｇｅＶＲ（ｘ、ｙ、ｎ）＋ｋ・Ｅｎｅｒｇｙ（ｎ）・Ｌｉｇｈｔ（ｎ）／Ｌ
ここで、ｋは強度倍率であり、初期値はｋ＝１．０に設定されている。

続いて、制御部１１は、ｘ←ｘ＋１に更新し（ステップＳ１１２）、ｘ＜Ｓｉｚｅの場合、ステップＳ１０２に戻り、次の画素ｘのＲＧＢ値を決定する処理（ステップＳ１０２～Ｓ１１１）を実行する。ｘ≧Ｓｉｚｅの場合、ｙ←ｙ＋１に更新し（ステップＳ１１３）、ｙ≦Ｓｉｚｅの場合、ｘ＝０としたうえで、ステップＳ１０２に戻り、ｙ行目の画素（ｘ、ｙ）のＲＧＢ値を決定する処理（ステップＳ１０２～Ｓ１１１）を繰り返す。ｙ＞Ｓｉｚｅの場合、制御部１１は、処理を終了する。
すなわち、全画素のＲＧＢ値が算出されるまで（ステップＳ１１２；ｘ≧Ｓｉｚｅ、かつ、ステップＳ１１３；ｙ＞Ｓｉｚｅ）、ステップＳ１０２～Ｓ１１１の処理を繰り返す。

図２７のフローチャートに戻る。
制御部１１は、処理回数をｌ←ｌ＋１に更新し、サブサンプル・オフセット値をｄｘ←ｄｘ＋１／Ｌ、ｄｙ←ｄｙ＋１／Ｌ、ｄｚ←ｄｚ＋１／Ｌに更新する（ステップＳ９４）。制御部１１は、処理回数ｌがｌ＞Ｌ－１を満たすまで（ステップＳ９４；ｌ＞Ｌ－１）、ステップＳ９３のレイキャスティング処理を繰り返す。レイキャスティング処理が終了すると、制御部１１は、生成したボリュームレンダリング画像ＩｍａｇｅＶＲを出力する（ステップＳ９５）。

（１１－１．交点算出処理）
図２８のステップＳ１０３において実行される有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する処理について説明する。

制御部１１は、以下の処理手順で、レイキャスティング処理の算出対象画素（ｘ、ｙ）に対して、視点座標系において、ｚ＝Ｚｏ（＝Ｓｉｚｅ／Ｓｒａｙ－１）からｚ＝０の範囲で、視点ｚ＝Ｚｏに最も近い有効ボクセル領域Ｖｒ（［Ｘｉｓ、Ｘｉｅ］、［Ｙｉｓ、Ｙｉｅ］、［Ｚｉｓ、Ｚｉｅ］）との交点Ｚｃを算出する。

（１）視線ベクトル（仮想光線ベクトル）の算出
制御部１１は、図２９に示すように、視点座標系における視点座標（ｘ、ｙ、Ｚｏ）および下限座標（ｘ、ｙ、０）に対して各々座標変換を行い、ボクセル座標系における視点座標（ｘ１、ｙ１、ｚ１）および下限座標（ｘ２、ｙ２、ｚ２）を求め、視線ベクトル（ｖｘ、ｖｙ、ｖｚ）の各要素を次のように算出する。

ｖｘ＝ｘ２－ｘ１
ｖｙ＝ｙ２－ｙ１
ｖｚ＝ｚ２－ｚ１

（２）視線ベクトルと有効ボクセル領域Ｖｒとの交点の算出
制御部１１は、ｔｘ＝ｔｙ＝ｔｚ＝１０と初期化し、
１）｜ｖｘ｜≧１の場合、ｔｘ１＝（Ｘｉｓ－ｘ１）／ｖｘ、ｔｘ２＝（Ｘｉｅ－ｘ１）／ｖｘを算出し、いずれか小さい方をｔｘに設定する。
２）｜ｖｙ｜≧１の場合、ｔｙ１＝（Ｙｉｓ－ｙ１）／ｖｙ、ｔｙ２＝（Ｙｉｅ－ｙ１）／ｖｙを算出し、いずれか小さい方をｔｙに設定する。
３）｜ｖｚ｜≧１の場合、ｔｚ１＝（Ｚｉｓ－ｚ１）／ｖｚ、ｔｚ２＝（Ｚｉｅ－ｚ１）／ｖｚを算出し、いずれか小さい方をｔｚに設定する。

これにより、図２９（右図）に示すように、ボクセル座標系における交点座標が求まる。具体的には、交点座標は（ｖｘ・ｔ＋ｘ１，ｖｙ・ｔ＋ｙ１，ｖｚ・ｔ＋ｚ１）で与えられ（ｔは（２）で算出されたｔｘ、ｔｙ、ｔｚのいずれかの値）、図２９の例ではｚ＝Ｚｉｓおよびｚ＝Ｚｉｅの２面に交点が存在し、ｚ＝Ｚｉｅにおける交点座標は（ｖｘ・ｔｚ＋ｘ１，ｖｙ・ｔｚ＋ｙ１，ｖｚ・ｔｚ＋ｚ１）と算出できる。交点座標は、視線ベクトルと有効ボクセル領域Ｖｒ（直方体）を構成する６面との交点（通常、交点は２点以上存在する）のうち、視点に最も近い交点である。

（３）視点座標系における交点Ｚｃの決定
制御部１１は、視点座標系における交点座標ＺｃをＺｃ＝－１と初期化し、
１）ｔｘ≦１の場合
Ｙ＝ｖｙ・ｔｘ＋ｙ１、Ｚ＝ｖｚ・ｔｘ＋ｚ１を算出する。
Ｙｉｓ≦Ｙ≦ＹｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅの場合（ボクセル座標系における交点が有効ボクセル領域Ｖｒ内に存在する場合）、対応する視点座標系における交点Ｚ２＝Ｚｏ・（１－ｔｘ）を算出し、Ｚ２＞ＺｃならばＺｃ＝Ｚ２に置換する。
２）ｔｙ≦１の場合
Ｘ＝ｖｘ・ｔｙ＋ｘ１、Ｚ＝ｖｚ・ｔｙ＋ｚ１を算出する。
Ｘｉｓ≦Ｘ≦ＸｉｅかつＺｉｓ≦Ｚ≦Ｚｉｅの場合（ボクセル座標系における交点が有効ボクセル領域Ｖｒ内に存在する場合）、対応する視点座標系における交点Ｚ２＝Ｚｏ・（１－ｔｙ）を算出し、Ｚ２＞ＺｃならばＺｃ＝Ｚ２に置換する。
３）ｔｚ≦１の場合
Ｘ＝ｖｘ・ｔｚ＋ｘ１、Ｙ＝ｖｙ・ｔｚ＋ｙ１を算出する。
Ｘｉｓ≦Ｘ≦ＸｉｅかつＹｉｓ≦Ｙ≦Ｙｉｅの場合（ボクセル座標系における交点が有効ボクセル領域Ｖｒ内に存在する場合）、対応する視点座標系における交点Ｚ２＝Ｚｏ・（１－ｔｚ）を算出し、Ｚ２＞ＺｃならばＺｃ＝Ｚ２に置換する。

これにより、図２９（左図）に示すように、視点座標系における視点に最も近い交点座標（ｘ、ｙ、Ｚｃ）が求まる。

（４）視点より手前に算出される交点の補正
ｔｘ、ｔｙ、ｔｚは負値をとる場合があり、その場合は交点が視点より手前になる。制御部１１は、Ｚｃ＞Ｚｏの場合、Ｚｃ＝Ｚｏに補正する。すなわち、図３０のように、視点座標系における交点座標が視点より手前（不可視領域）に算出された場合、当該交点座標を視点位置に補正する。

（１１－２．起点座標探索処理）
図３１のフローチャートを参照して、図２８のステップＳ１０５において実行される、起点座標ｚを探索する処理について説明する。

まず、制御部１１は、ｚをｚ＝Ｚｓｔ（＝交点座標Ｚｃ）に初期化し、探索対象画素（ｘ、ｙ）を入力する（ステップＳ１２１）。続いて、制御部１１は、３次元座標（ｘ、ｙ、ｚ）に対して座標変換を行い、ボクセルのα値を取得する（ステップＳ１２２）。ボクセルのα値を取得する処理は後述する（「１１－３．ボクセルα値を取得（補間なし）」参照）。

ステップＳ１２２において取得したボクセルのα値がα＝０の場合（ボクセルが透明の場合）、制御部１１は、ｚ←ｚ－ｍ（例えばｍ＝４）に更新する（ステップＳ１２３）。ｍは探索時のＺ軸方向のスキップ幅であり、ｍが大きいほど高速に探索できるが、スキップ幅が大きすぎると起点座標の探索精度が悪化する場合がある。

ステップＳ１２３において更新したｚがｚ＜０の場合（ボクセル空間外の場合）、制御部１１は、ｚ＝－１（起点座標が存在しないことを示す値）を返し（ステップＳ１２４）、処理を終了する。ｚ≧０の場合、制御部１１は、ステップＳ１２２に戻る。

ステップＳ１２２において取得したボクセルのα値がα＜０の場合（有効ボクセル領域Ｖｒを通過した場合）は、制御部１１は、ｚ＝－１（起点座標が存在しないことを示す値）を返し（ステップＳ１２４）、処理を終了する。

ステップＳ１２２において取得したボクセルのα値がα＞０の場合（ボクセルが不透明の場合）、以下のように、Ｚ軸の＋方向（視線逆方向）に戻って正確な起点座標を探索する。

制御部１１は、まず、ｉをｉ＝０に初期化し（ステップＳ１２５）、ｉ←ｉ＋１に更新する（ステップＳ１２６）。

ｉ≧ｍまたはｚ＋ｉ≧Ｚｓｔの場合、制御部１１は、ｚ＝ｚ＋ｉ－１を返し（ステップＳ１２８）、処理を終了する。

ｉ＜ｍかつｚ＋ｉ＜Ｚｓｔの場合、３次元座標（ｘ、ｙ、ｚ＋ｉ）に対して座標変換を行い、ボクセルのα値を取得する（ステップＳ１２７）。ボクセルのα値を取得する処理は後述する（「１１－３．ボクセルα値を取得（補間なし）」参照）。

ステップＳ１２７において取得したボクセルのα値がα＝０の場合（ボクセルが透明の場合）、制御部１１は、ｚ＝ｚ＋ｉ－１を返し（ステップＳ１２８）、処理を終了する。

ステップＳ１２７において取得したボクセルのα値がα＞０の場合（ボクセルが不透明の場合）、ステップＳ１２６に戻り、ステップＳ１２６～Ｓ１２７の処理を繰り返す。

以上のように、先頭の不透明ボクセルのＺ座標ｚ（仮想光線の照射を開始する起点座標）を所定のステップ幅で探索することで（ステップＳ１２２、Ｓ１２３）、起点座標を高速に特定できる。また、不透明ボクセルが探索された場合、当該ボクセルの位置とステップ幅だけ視線逆方向（Ｚ軸の＋方向）に戻った位置との間に存在する可能性がある不透明ボクセルを再探索することで（ステップＳ１２６、Ｓ１２７）、起点座標を正確に特定することができる。

（１１－３．ボクセルα値を取得（補間なし））
図３２のフローチャートを参照して、起点座標探索処理（図３１）のステップＳ１２２およびステップＳ１２７において実行される、ボクセルのα値を取得する処理について説明する。

まず、制御部１１は、与えられた視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対して座標変換を行い、対応するボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を算出する（ステップＳ１３１）。ここで座標変換とは、視点座標系をボクセル座標系に変換する処理であり、ＧＵＩ側の変換処理とは逆になる。ＧＵＩ側では関心領域ＲＯＩによるクリッピング、スケーリング、Ｚ軸方向変倍処理、オフセット（ＸＹＺ軸方向同時）、回転、透視変換の順に行うものと仮定し、制御部１１は、与えられた視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を以下のように算出する。

まず、制御部１１は、視点座標系の座標値（ｘ、ｙ、ｚ）（整数値）を次のように実数値（ｘｘ、ｙｙ、ｚｚ）に変換する。視点座標系はＳｉｚｅ（Ｘ軸方向のサイズ）×Ｓｉｚｅ（Ｙ軸方向のサイズ）×Ｓｉｚｅ／Ｓｒａｙ（Ｚ軸方向のサイズ）の直方体で、Ｓｘ×Ｓｙ×Ｓｚの断層画像群Ｄｏとは独立して定義される。通常Ｓｘ＝Ｓｙ＝５１２であればＳｉｚｅ＝５１２となる。視点座標系のＺ軸方向は仮想光線サブサンプル１／Ｓｒａｙ倍に伸縮されていることを考慮して、はじめに以下のようにオフセット処理（ＸＹＺ軸方向同時）を行う。

（式３０）
ｘｘ＝ｘ－Ｓｉｚｅ／２＋ｄｘ＋Ｘｏｆｆ
ｙｙ＝ｙ－Ｓｉｚｅ／２＋ｄｙ＋Ｙｏｆｆ
ｚｚ＝ｚ・Ｓｒａｙ－Ｓｉｚｅ／２＋ｄｚ＋Ｚｏｆｆ

続いて、制御部１１は、以下のように、回転パラメータ行列Ｒを用いて回転処理（座標変換）を行う。

（式３１）
ｘｘ’＝Ｒ１１・ｘｘ＋Ｒ１２・ｙｙ＋Ｒ１３・ｚｚ
ｙｙ’＝Ｒ２１・ｘｘ＋Ｒ２２・ｙｙ＋Ｒ２３・ｚｚ
ｚｚ’＝Ｒ３１・ｘｘ＋Ｒ３２・ｙｙ＋Ｒ３３・ｚｚ
回転処理後の（ｘｘ’、ｙｙ’、ｚｚ’）を（ｘｘ、ｙｙ、ｚｚ）とする。

そして、制御部１１は、スケーリング、ＸＹＺ軸方向変倍処理を同時に行い、以下のように、ボクセルの座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）を算出する。

（式３２）
ｘｒ＝ｘｘ／Ｓｃａｌｅ／Ｓｃｘ＋Ｓｘ／２
ｙｒ＝ｙｙ／Ｓｃａｌｅ／Ｓｃｙ＋Ｓｙ／２
ｚｒ＝ｚｚ／Ｓｃａｌｅ／Ｓｃｚ＋Ｓｚ／２

続いて、制御部１１は、算出したボクセルの座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）に対して、各値に０．５を加算して、以下のように、小数点以下を切り捨て整数化した座標値（ｘｉ、ｙｉ、ｚｉ）（四捨五入した整数値）に変換する（ステップＳ１３２）。

（式３３）
ｘｉ＝ＩＮＴ［ｘｒ＋０．５］
ｙｉ＝ＩＮＴ［ｙｒ＋０．５］
ｚｉ＝ＩＮＴ［ｚｒ＋０．５］

そして、制御部１１は、有効ボクセル領域Ｖｒ、マスクデータＭａｓｋを考慮して、以下のようにボクセルα値を取得する（ステップＳ１３３）。

（式３４）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外の場合）
α＝－１
２）上記１）を満たさない場合
α＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ、ｚｉ）、３）・Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）

（１１－４．ボクセルα値、ＲＧＢ値の取得（補間あり））
レイキャスティング処理（図２８）のステップＳ１０６、Ｓ１０８において実行される、ボクセルα値（Ｖα）、ボクセルＲＧＢ値を取得する処理を説明する。

（１１－４－１．ボクセルα値の取得）
図３３は、ボクセルα値を取得する処理を示すフローチャートである。

まず、制御部１１は、図３２のステップ１３１と同様の方法で、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対して座標変換を行い、ボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を算出する（ステップＳ１４１）。

続いて、制御部１１は、算出したボクセルの座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）を、小数点以下を切り捨て整数化した座標値（ｘｉ、ｙｉ、ｚｉ）（整数値）に変換する（ステップＳ１４２）。ここで、切り捨てた小数点以下の端数を（ｗｘ、ｗｙ、ｗｚ）（０≦ｗｘ、ｗｙ、ｗｚ＜１）とする（すなわち、ｘｒ＝ｘｉ＋ｗｘ、ｙｒ＝ｙｉ＋ｗｙ、ｚｒ＝ｚｉ＋ｗｚ）。

そして、制御部１１は、有効ボクセル領域Ｖｒを考慮して、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応する補間対象のボクセルの信号値を断層画像群Ｄｏに基づいて以下のように抽出する（ステップＳ１４３）。

（式３５）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｄ０００＝－３２７６９（無効の値）
２）上記１）の条件を満たさない場合において、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（中央ボクセルを抽出）
Ｄ０００＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
３）上記１）２）の条件を満たさない場合（８近傍ボクセルの信号値を抽出）
Ｄ０００＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
Ｄ１００＝Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ）
Ｄ０１０＝Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ）
Ｄ１１０＝Ｄｏ（ｘｉ＋１、ｙｉ＋１、ｚｉ）
Ｄ００１＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ＋１）
Ｄ１０１＝Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ＋１）
Ｄ０１１＝Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ＋１）
Ｄ１１１＝Ｄｏ（ｘｉ＋１、ｙｉ＋１、ｚｉ＋１）

続いて、制御部１１は、マスクデータＭａｓｋに基づいて、抽出した補間対象の各ボクセルのα値（不透明度）に対するマスク値Ｓ０００～Ｓ１１１（０または１）を以下のように設定する（ステップＳ１４４）。すなわち、各ボクセルのα値（不透明度）を計算に用いるか否かを設定する。

（式３６）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｓ０００＝０
２）上記１）の条件を満たさない場合において、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（中央ボクセルのマスク値を設定）
Ｓ０００＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）
３）上記１）２）の条件を満たさない場合（８近傍ボクセルのマスク値を設定）
Ｓ０００＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）
Ｓ１００＝Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ）
Ｓ０１０＝Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ）
Ｓ１１０＝Ｍａｓｋ（ｘｉ＋１、ｙｉ＋１、ｚｉ）
Ｓ００１＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ＋１）
Ｓ１０１＝Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ＋１）
Ｓ０１１＝Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ＋１）
Ｓ１１１＝Ｍａｓｋ（ｘｉ＋１、ｙｉ＋１、ｚｉ＋１）

そして、制御部１１は、抽出した補間対象ボクセルの信号値に基づいて、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセルのα値Ｖαを以下のように取得する（ステップＳ１４５）。

（式３７）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｖα＝－１
２）上記１）の条件を満たさず、ｘｉ＝Ｘｉｓ、ｘｉ＝Ｘｉｅ、ｙｉ＝Ｙｉｓ、ｙｉ＝Ｙｉｅ、ｚｉ＝Ｚｉｓ、又はｚｉ＝Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒの境界面）
Ｖα＝０
３）上記１）２）の条件を満たさず、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（補間しない）
Ｖα＝Ｃｍａｐ（Ｄ０００、３）・Ｓ０００
４）上記１）２）３）の条件を満たさない場合（補間する）
Ｖα＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｃｍａｐ（Ｄ０００、３）・Ｓ０００＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１００、３）・Ｓ１００＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｃｍａｐ（Ｄ０１０、３）・Ｓ０１０＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１０、３）・Ｓ１１０＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｃｍａｐ（Ｄ００１、３）・Ｓ００１＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１０１、３）・Ｓ１０１＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｃｍａｐ（Ｄ０１１、３）・Ｓ０１１＋ｗｚ・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１１、３）・Ｓ１１１

（１１－４－２．ボクセルＲＧＢ値の取得）
制御部１１は、ステップＳ１４３において抽出した補間対象ボクセルの信号値に基づいて、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセルＲＧＢ値（Ｖｃ（ｎ））（０≦ｎ≦２）を以下のように取得する。なお、ＲＧＢ値を取得する場合には、Ｓ０００、Ｓ１００、・・・、Ｓ１１１を乗算するマスク処理を行わない。これにより、マスク境界面でのモアレの発生を抑制する。

（式３８）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｖｃ（ｎ）＝０ （０≦ｎ≦２）
２）上記１）の条件を満たさず、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（補間しない）
Ｖｃ（ｎ）＝Ｃｍａｐ（Ｄ０００、ｎ） （０≦ｎ≦２）
３）上記１）２）の条件を満たさない場合（補間する）
Ｖｃ（ｎ）＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｃｍａｐ（Ｄ０００、ｎ）＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１００、ｎ）＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｃｍａｐ（Ｄ０１０、ｎ）＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１０、ｎ）＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｃｍａｐ（Ｄ００１、ｎ）＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｃｍａｐ（Ｄ１０１、ｎ）＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｃｍａｐ（Ｄ０１１、ｎ）＋ｗｚ・ｗｙ・ｗｘ・Ｃｍａｐ（Ｄ１１１、ｎ） （０≦ｎ≦２）

（１１－５．陰影計算）
制御部１１は、必要に応じて、次のように陰影計算を行ってもよい。
まず、光源ベクトル（Ｌｘ、Ｌｙ、Ｌｚ）（単位ベクトル）を設定する。例えば、（Ｌｘ、Ｌｙ、Ｌｚ）＝（０．５７７３５、０．５７７３５、０．５７７３５）と設定する。また、環境光成分Ａｂ（０≦Ａｂ≦１、例えばＡｂ＝０．２）を設定する。

続いて、制御部１１は、次のように、座標変換により算出された３次元座標（ｘｉ、ｙｉ、ｚｉ）における６近傍ボクセルの不透明度Ｖ１００、Ｖ２００、Ｖ０１０、Ｖ０２０、Ｖ００１、Ｖ００２を取得する。但し、ｘｉ＋１＞ＸｉｅのときＶ１００＝０、ｘｉ－１＜ＸｉｓのときＶ２００＝０、ｙｉ＋１＞ＹｉｅのときＶ０１０＝０、ｙｉ－１＜ＹｉｓのときＶ０２０＝０、ｚｉ＋１＜ＺｉｓのときＶ００１＝０、ｚｉ－１＜ＺｉｓのときＶ００２＝０とする。すなわち、隣接するボクセルがクリッピング範囲の場合は不透明度を０として扱う。

（式３９）
Ｖ１００＝Ｃｍａｐ（Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ）、３）・Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ）
Ｖ２００＝Ｃｍａｐ（Ｄｏ（ｘｉ－１、ｙｉ、ｚｉ）、３）・Ｍａｓｋ（ｘｉ－１、ｙｉ、ｚｉ）
Ｖ０１０＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ）、３）・Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ）
Ｖ０２０＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ－１、ｚｉ）、３）・Ｍａｓｋ（ｘｉ、ｙｉ－１、ｚｉ）
Ｖ００１＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ、ｚｉ＋１）、３）・Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ＋１）
Ｖ００２＝Ｃｍａｐ（Ｄｏ（ｘｉ、ｙｉ、ｚｉ－１）、３）・Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ－１）

続いて、制御部１１は、座標変換により算出された３次元座標（ｘｉ、ｙｉ、ｚｉ）における勾配ベクトル（Ｇｘ、Ｇｙ、Ｇｚ）を、以下の式で算出する。

（式４０）
Ｇｘ＝（Ｖ１００－Ｖ２００）・Ｓｃｘ
Ｇｙ＝（Ｖ０１０－Ｖ０２０）・Ｓｃｙ
Ｇｚ＝（Ｖ００１－Ｖ００２）・Ｓｃｚ
Ｇ＝｛Ｇｘ^２＋Ｇｙ^２＋Ｇｚ^２｝^１／２

続いて、Ｇ≧１の場合、制御部１１は、輝度値（陰影値）Ｓ（０≦Ｓ≦１）として拡散反射成分を算出し、

（式４１）
Ｓ＝（１－Ａｂ）｜Ｇｘ・Ｌｘ＋Ｇｙ・Ｌｙ＋Ｇｚ・Ｌｚ｜／Ｇ＋Ａｂ

を与える。Ｇ＜１の場合（αが変化しない場合）、制御部１１は、輝度値（陰影値）Ｓとして、Ｓ＝０を与える。

そして、制御部１１は、

（式４２）
Ｖｃ’（ｎ）＝Ｓ・Ｖｃ（ｎ）

とし、算出されたＲＧＢ値Ｖｃ（ｎ）（０≦ｎ≦２）の成分を改変する。

以上、図２５～図３３を参照しながら、ボリュームレンダリング画像を生成する処理を説明した。マスクデータＭａｓｋに基づいて断層画像群Ｄｏの各ボクセルの不透明度に対するマスク値を設定したうえで、断層画像群Ｄｏの各ボクセルにカラーマップＣｍａｐ（オパシティカーブ）を適用して、ボクセルα値を取得（補間計算）することで、マスクデータＭａｓｋの描画・非描画情報が反映されたボリュームレンダリング画像を生成することができる。

（１２．ＭＩＰ画像の生成）
次に、図３４～図４０を参照しながら、ＭＩＰ画像を生成する処理を説明する。

図３４は、ＭＩＰ画像を生成する処理の全体の流れを示すフローチャートである。

まず、制御部１１は、必要に応じて、ユーザからクリッピング領域ＲＯＩ（（式２２）参照）の設定を受け付ける（ステップＳ１６１）。

続いて、制御部１１は、有効ボクセル領域Ｖｒを算出する（ステップＳ１６２）。具体的には、図３５に示すように、制御部１１は、信号値が所定の閾値以上（例えば、空気でない信号値以上）であるボクセルを有効ボクセルとし、これらの有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。クリッピング領域ＲＯＩやマスクデータＭａｓｋが設定されている場合は、制御部１１は、クリッピングＲＯＩやマスクデータＭａｓｋにより規定された描画範囲内において、信号値が所定の閾値以上（例えば、空気でない信号値以上）であるボクセルを有効ボクセルとし、これらの有効ボクセルを全て含み、かつ有効ボクセルに外接する直方体を有効ボクセル領域Ｖｒとして算出する。

そして、制御部１１は、ステップＳ７１（図２４）において取得した断層画像群Ｄｏ（３次元ボクセル）、ステップＳ７３（図２４）において生成されたマスクデータＭａｓｋ、ステップＳ１６１（図３４）において設定されたクリッピング領域ＲＯＩ、ステップＳ１６２（図３４）において算出された有効ボクセル領域Ｖｒを参照して、ＭＩＰ画像ＩｍａｇｅＭＩＰを生成する（ステップＳ１６３）。

特に本開示では、制御部１１は、レイキャスティング処理において、マスクデータＭａｓｋに基づいて断層画像群Ｄｏの各ボクセルの信号値に対するマスク値を設定したうえで（後述する図４０のステップＳ２２４）、断層画像群Ｄｏの各ボクセルの信号値を取得（補間計算）し（後述する図４０のステップＳ２２５）、ＭＩＰ画像ＩｍａｇｅＭＩＰを生成する。これにより、被写体形状（幾何形状）に基づいて生成したマスクデータＭａｓｋの描画・非描画情報が反映されたＭＩＰ画像が生成される。

図３６は、ステップＳ１６３において実行される、ＭＩＰ画像ＩｍａｇｅＭＩＰを生成する具体的な処理の流れを示す。
まず、制御部１１は、座標変換パラメータを設定する（ステップＳ１７１）。座標変換パラメータの設定は、前記した図２７のステップＳ９１と同様である。

続いて、制御部１１は、レイキャスティング処理（各座標（ｘ、ｙ）毎に代表信号値を算出する処理）を実行する（ステップＳ１７２）。

＜レイキャスティング処理＞
図３７のフローチャートを参照して、レイキャスティング処理について説明する。生成するＭＩＰ画像ＩｍａｇｅＭＩＰ（ｘ、ｙ）は断層画像群Ｄｏの階調に基づいてモノクロ８ビットまたは１６ビットの値をもつが、図２４のステップＳ７１において断層画像群Ｄｏに対する階調圧縮処理を施さない場合が多く、その場合はモノクロ１６ビットになる。以下、モノクロ（１６ビット）のＭＩＰ画像ＩｍａｇｅＭＩＰ（ｘ、ｙ）を生成する場合について説明する。

制御部１１は、まず、生成するモノクロ（１６ビット）のＭＩＰ画像ＩｍａｇｅＭＩＰ（ｘ、ｙ）の初期値を全て０に設定する（ＩｍａｇｅＭＩＰ（ｘ、ｙ）＝０）。また、代表信号値算出モードｍｏｄｅ（０：ＭＩＰ、１：ＭｉｎＩＰ、２：ＲａｙＳｕｍ）を設定し、各２次元座標（ｘ、ｙ）（０≦ｘ≦Ｓｉｚｅ－１、０≦ｙ≦Ｓｉｚｅ－１）に対して、以下の処理を実行する。

制御部１１は、ｘ＝０、ｙ＝０と設定し（ステップＳ１８１）、有効ボクセル領域Ｖｒと仮想光線との交点のＺ座標（Ｚｃ）を算出する（ステップＳ１８２）。交点のＺ座標を算出する処理は、前記した図２８のステップＳ１０３の処理と同様である（「１１－１．交点算出処理」参照）。

続いて、制御部１１は、代表信号値算出モードｍｏｄｅに応じて、代表信号値Ｖｍとカウンタｃｎｔを以下のように初期化し、Ｚｓｔ＝Ｚｃとする（ステップＳ１８３）。

（式４３）
ｍｏｄｅ＝０（ＭＩＰ）：Ｖｍ＝－３２７６８（最小値を設定）
ｍｏｄｅ＝１（ＭｉｎＩＰ）：Ｖｍ＝３２７６７（最大値を設定）
ｍｏｄｅ＝２（ＲａｙＳｕｍ）：Ｖｍ＝０、ｃｎｔ＝０

続いて、制御部１１は、Ｚｓｔから後続する先頭の有効ボクセル（信号値が所定の閾値以上のボクセル）のＺ座標ｚ（起点座標）を探索する（ステップＳ１８４）。すなわち、ステップＳ１８２において算出された交点から視線方向に向かってレイキャスティング計算を開始する起点座標を探索する。起点座標を探索する処理については後述する（「１２－１．起点座標探索処理」、図３８参照）。

ｚ＜０の場合、ステップＳ１８９へ移行する。

一方、ｚ≧０の場合、制御部１１は、３次元座標（ｘ、ｙ、ｚ）を座標変換してボクセル信号値Ｖｓを取得する（ステップＳ１８５）。ボクセル信号値Ｖｓを取得する処理については後述する（「１２－３．ボクセル信号値Ｖｓを取得（補間あり）」参照）。

ステップＳ１８５において取得したボクセル信号値ＶｓがＶｓ＝－３２７６９の場合（有効ボクセル領域Ｖｒを通過した場合）、ステップＳ１８９へ移行する。これにより、有効ボクセル領域Ｖｒを通過している場合（Ｖｓ＝－３２７６９）は、その先に描画対象のボクセルは存在しないため、レイキャスティング処理を早期に打ち切り、冗長な処理を省略する。

ステップＳ１８５において取得したボクセル信号値ＶｓがＶｓ＜－３２７６９の場合（無効値の場合）、Ｚｓｔ＝ｚ－１としたうえで（ステップＳ１８６）、ステップＳ１８４に戻り、起点座標を再度探索する。

ステップＳ１８５において取得したボクセル信号値ＶｓがＶｓ≧－３２７６８の場合（有効値の場合）、制御部１１は、代表信号値算出モードｍｏｄｅに応じて、以下のように代表信号値Ｖｍを算出する（ステップＳ１８７）。

（式４４）
ｍｏｄｅ＝０（ＭＩＰ ｍｏｄｅ）：
Ｖｓ＞ＶｍならばＶｍ＝Ｖｓ
ｍｏｄｅ＝１（ＭｉｎＩＰ ｍｏｄｅ）：
Ｖｓ＜ＶｍならばＶｍ＝Ｖｓ
ｍｏｄｅ＝２（ＲａｙＳｕｍ ｍｏｄｅ）：
Ｖｍ←Ｖｍ＋Ｖｓ、ｃｎｔ←ｃｎｔ＋１

ｚ＝ｚ－１に更新し（ステップＳ１８８）、ｚ≧０の場合、ステップＳ１８５に戻り、ｚ＜０の場合、ステップＳ１８９へ移行する。

ステップ１８９において、ｍｏｄｅ＝０、１の場合（ステップＳ１８９；ｍｏｄｅ＜２）には、制御部１１は、ステップＳ１８７で得られた代表信号値Ｖｍを画素値として記録する（ＩｍａｇｅＭＩＰ（ｘ、ｙ）＝Ｖｍ、ステップＳ１９１）。一方、ｍｏｄｅ＝２（ＲａｙＳｕｍ ｍｏｄｅ）の場合（ステップＳ１８９；ｍｏｄｅ＝２）には、制御部１１は、Ｖｍ←Ｖｍ／ｃｎｔのように信号値の平均を計算したうえで（ステップＳ１９０）、代表信号値Ｖｍを画素値として記録する（ＩｍａｇｅＭＩＰ（ｘ、ｙ）＝Ｖｍ、ステップＳ１９１）。

続いて、制御部１１は、ｘ←ｘ＋１に更新し（ステップＳ１９２）、ｘ＜Ｓｉｚｅの場合、ステップＳ１８２に戻り、次の画素ｘの代表信号値Ｖｍを求める処理（ステップＳ１８２～Ｓ１９１）を実行する。ｘ≧Ｓｉｚｅの場合、ｙ←ｙ＋１に更新し（ステップＳ１９３）、ｙ≦Ｓｉｚｅの場合、ｘ＝０としたうえで、ステップＳ１８２に戻り、ｙ行目の画素（ｘ、ｙ）の代表信号値Ｖｍを求める処理（ステップＳ１８２～Ｓ１９１）を繰り返す。ｙ＞Ｓｉｚｅの場合、制御部１１は、処理を終了する。
すなわち、画素領域の全画素の代表信号値Ｖｍが得られるまで（ステップＳ１９２；ｘ≧Ｓｉｚｅ、かつ、ステップＳ１９３；ｙ＞Ｓｉｚｅ）、ステップＳ１８２～Ｓ１９１の処理を繰り返す。

図３６のフローチャートに戻る。
制御部１１は、ステップＳ１７２において生成したＭＩＰ画像ＩｍａｇｅＭＩＰを出力する（ステップＳ１７３）。

（１２－１．起点座標探索処理）
図３８のフローチャートを参照して、図３７のステップＳ１８４において実行される、起点座標を探索する処理について説明する。

まず、制御部１１は、ｚをｚ＝Ｚｓｔと初期化し、探索対象画素（ｘ、ｙ）を入力する（ステップＳ２０１）。続いて、制御部１１は、３次元座標（ｘ、ｙ、ｚ）に対して座標変換を行い、ボクセル信号値Ｖｓを取得する（ステップＳ２０２）。ボクセル信号値Ｖｓを取得する処理は後述する（「１２－２．ボクセル信号値Ｖｓを取得（補間なし）」参照）。

ステップＳ２０２において取得したボクセル信号値ＶｓがＶｓ＝－３２７６９の場合（有効ボクセル領域Ｖrを通過した場合）、制御部１１は、ｚ＝－１（起点座標が存在しないことを示す値）を返し（ステップＳ２０４）、処理を終了する。

一方、ステップＳ２０２において取得したボクセル信号値Ｖｓ＜－３２７６９の場合（無効値の場合）は、制御部１１は、ｚ←ｚ－１に更新する（ステップＳ２０３）。

ステップＳ２０３において更新したｚがｚ＜０の場合、制御部１１は、ｚ＝－１（起点座標が存在しないことを示す値）を返し（ステップＳ２０４）、処理を終了する。ｚ≧０の場合、制御部１１は、ステップＳ２０２に戻る。

一方、ステップＳ２０２において取得したボクセル信号値ＶｓがＶｓ≧－３２７６８の場合（有効ボクセルの場合）、制御部１１は、ｚを起点座標として出力する。

（１２－２．ボクセル信号値Ｖｓを取得（補間なし））
図３９のフローチャートを参照して、起点座標探索処理（図３８）のステップＳ２０２において実行される、ボクセル信号値Ｖｓを取得する処理について説明する。

まず、制御部１１は、図３２のステップＳ１３１と同様の方法で、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対して座標変換を行い、ボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を算出し（ステップＳ２１１）、図３２のステップＳ１３２と同様の方法で、算出したボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）を整数の座標値（ｘｉ、ｙｉ、ｚｉ）に変換する（ステップＳ２１２）。

そして、制御部１１は、有効ボクセル領域Ｖｒ、マスクデータＭａｓｋを考慮して、以下のようにボクセル信号値Ｖｓを取得する（ステップＳ２１３）。

（式４５）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外の場合）
Ｖｓ＝－３２７６９
２）上記１）を満たさず、Ｍａｓｋ（ｘ、ｙ、ｚ）＝０の場合（無効値）
Ｖｓ＝－９９９９９
３）上記１）２）を満たさない場合（有効値）
Ｖｓ＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）

（１２－３．ボクセル信号値Ｖｓの取得（補間あり））
図４０のフローチャートを参照して、レイキャスティング処理（図３７）のステップＳ１８５において実行される、ボクセル信号値Ｖｓを取得する処理を説明する。

まず、制御部１１は、図３３のステップＳ１４１と同様の方法で、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対して座標変換を行い、ボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）を算出し（ステップＳ２２１）、図３３のステップＳ１４２と同様の方法で、算出したボクセルの実数の座標値（ｘｒ、ｙｒ、ｚｒ）（実数値）を、小数点以下を切り捨て整数化した座標値（ｘｉ、ｙｉ、ｚｉ）に変換する（ステップＳ２２２）。切り捨てた小数点以下の端数を（ｗｘ、ｗｙ、ｗｚ）（０≦ｗｘ、ｗｙ、ｗｚ＜１）とする（すなわち、ｘｒ＝ｘｉ＋ｗｘ、ｙｒ＝ｙｉ＋ｗｙ、ｚｒ＝ｚｉ＋ｗｚ）。

続いて、制御部１１は、以下のように、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応する補間対象のボクセルの信号値Ｄ０００、Ｄ１００、Ｄ０１０、Ｄ００１、Ｄ１０１、Ｄ０１１、Ｄ１１１を抽出する（ステップＳ２２３）。

（式４６）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｄ０００＝－３２７６９（無効の値）
２）上記１）の条件を満たさない場合において、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（中央ボクセルを抽出）
Ｄ０００＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
３）上記１）２）の条件を満たさない場合（８近傍ボクセルの信号値を抽出）
Ｄ０００＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ）
Ｄ１００＝Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ）
Ｄ０１０＝Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ）
Ｄ１１０＝Ｄｏ（ｘｉ＋１、ｙｉ＋１、ｚｉ）
Ｄ００１＝Ｄｏ（ｘｉ、ｙｉ、ｚｉ＋１）
Ｄ１０１＝Ｄｏ（ｘｉ＋１、ｙｉ、ｚｉ＋１）
Ｄ０１１＝Ｄｏ（ｘｉ、ｙｉ＋１、ｚｉ＋１）
Ｄ１１１＝Ｄｏ（ｘｉ＋１、ｙｉ＋１、ｚｉ＋１）

続いて、制御部１１は、マスクデータＭａｓｋに基づいて、抽出した補間対象の各ボクセルの信号値に対するマスク値Ｓ０００～Ｓ１１１（０または１）を以下のように設定する（ステップＳ２２４）。すなわち、各ボクセルの信号値を計算に用いるか否かを設定する。

（式４７）
１）１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｓ０００＝０
２）上記１）の条件を満たさない場合において、ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、又はｚｉ＋１＞Ｚｉｅのいずれかを満たす場合（中央ボクセルのマスク値を設定）
Ｓ０００＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）
３）上記１）２）の条件を満たさない場合（８近傍ボクセルのマスク値を設定）
Ｓ０００＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ）
Ｓ１００＝Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ）
Ｓ０１０＝Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ）
Ｓ１１０＝Ｍａｓｋ（ｘｉ＋１、ｙｉ＋１、ｚｉ）
Ｓ００１＝Ｍａｓｋ（ｘｉ、ｙｉ、ｚｉ＋１）
Ｓ１０１＝Ｍａｓｋ（ｘｉ＋１、ｙｉ、ｚｉ＋１）
Ｓ０１１＝Ｍａｓｋ（ｘｉ、ｙｉ＋１、ｚｉ＋１）
Ｓ１１１＝Ｍａｓｋ（ｘｉ＋１、ｙｉ＋１、ｚｉ＋１）

そして、制御部１１は、抽出した補間対象ボクセルの信号値に基づいて、視点座標系の３次元座標値（ｘ、ｙ、ｚ）（整数値）に対応するボクセル信号値Ｖｓを以下のように取得する（ステップＳ２２５）。

（式４８）
１）ｘｉ＜Ｘｉｓ、ｘｉ＞Ｘｉｅ、ｙｉ＜Ｙｉｓ、ｙｉ＞Ｙｉｅ、ｚｉ＜Ｚｉｓ、又はｚｉ＞Ｚｉｅのいずれかを満たす場合（有効ボクセル領域Ｖｒ外）
Ｖｓ＝－３２７６９
２）上記１）の条件を満たさず、Ｓ０００＝０の場合（無効値）
Ｖｓ＝－９９９９９
３）上記１）２）を満たさない場合において、
ｘｉ＋１＞Ｘｉｅ、ｙｉ＋１＞Ｙｉｅ、またはｚｉ＋１＞Ｚｉｅのいずれかを満たすか、或いはＳ０００、Ｓ１００、Ｓ０１０、Ｓ１１０、Ｓ００１、Ｓ１０１、Ｓ０１１、Ｓ１１１のいずれかが０の場合（補間しない）
Ｖｓ＝Ｄ０００
４）上記１）２）３）の条件を満たさない場合（補間する）
Ｖｓ＝（１－ｗｚ）（１－ｗｙ）（１－ｗｘ）・Ｄ０００＋（１－ｗｚ）（１－ｗｙ）・ｗｘ・Ｄ１００＋（１－ｗｚ）・ｗｙ・（１－ｗｘ）・Ｄ０１０＋（１－ｗｚ）・ｗｙ・ｗｘ・Ｄ１１０＋ｗｚ・（１－ｗｙ）（１－ｗｘ）・Ｄ００１＋ｗｚ・（１－ｗｙ）・ｗｘ・Ｄ１００＋ｗｚ・ｗｙ・（１－ｗｘ）・Ｄ０１１＋ｗｚ・ｗｙ・ｗｘ・Ｄ１１１

以上、モノクロ（１６ビット）のＭＩＰ画像ＩｍａｇｅＭＩＰ（ｘ、ｙ）を生成する場合について説明したが、制御部１１は、生成されたＭＩＰ画像ＩｍａｇｅＭＩＰ（ｘ、ｙ）を８ビットに変換したうえで、表示する。具体的には、ＤＩＣＯＭ規格のＷｉｎｄｏｗ Ｗｉｄｔｈ（ウィンドウ幅：ＷＷ）、Ｗｉｎｄｏｗ Ｌｅｖｅｌ（ウィンドウレベル：ＷＬ）を設定し、上限値Ｌｍａｘ＝ＷＬ＋ＷＷ／２、下限値Ｌｍｉｎ＝ＷＬ－ＷＷ／２に変換し、以下のようにＭＩＰ画像ＩｍａｇｅＭＩＰ（ｘ、ｙ）を８ビットに変換する。

ＩｍａｇｅＭＩＰ‘（ｘ、ｙ）＝（ＩｍａｇｅＭＩＰ（ｘ、ｙ）－Ｌｍｉｎ）・２５５／（Ｌｍａｘ－Ｌｍｉｎ）
ＩｍａｇｅＭＩＰ‘（ｘ、ｙ）＜０の場合、ＩｍａｇｅＭＩＰ‘（ｘ、ｙ）＝０とする。
ＩｍａｇｅＭＩＰ‘（ｘ、ｙ）＞２５５の場合、ＩｍａｇｅＭＩＰ‘（ｘ、ｙ）＝２５５とする。

以上、図３４～図４０を参照しながら、ＭＩＰ画像を生成する処理を説明した。マスクデータＭａｓｋに基づいて断層画像群Ｄｏの各ボクセルの信号値に対するマスク値を設定したうえで（図４０のステップＳ２２４）、断層画像群Ｄｏの各ボクセルの信号値を取得（補間計算）することで（図４０のステップＳ２２５）、マスクデータＭａｓｋの描画・非描画情報が反映されたＭＩＰ画像を生成することができる。

（１３．ＭＰＲ画像の生成）
次に、図４１、図４２を参照しながら、ＭＰＲ画像を生成する処理を説明する。

図４１は、ＭＰＲ画像を生成する処理の流れを示すフローチャートである。

まず、断層画像群Ｄｏの信号値が１６ビットの場合（（式１）参照）には、制御部１１は、断層画像群Ｄｏを８ビットに変換する（ステップＳ２４１）。具体的には、制御部１１は、断層画像群Ｄｏ（－３２７６８≦Ｄｏ（ｘ、ｙ、ｚ）≦３２７６７、０≦ｘ≦Ｓｘ－１、０≦ｙ≦Ｓｙ－１、０≦ｚ≦Ｓｚ－１；３方向の変倍率（Ｓｃｘ、Ｓｃｙ、Ｓｃｚ））に対して、ＤＩＣＯＭ規格のＷｉｎｄｏｗ Ｗｉｄｔｈ（ウィンドウ幅：ＷＷ）、Ｗｉｎｄｏｗ Ｌｅｖｅｌ（ウィンドウレベル：ＷＬ）を設定し、上限値Ｌｍａｘ＝ＷＬ＋ＷＷ／２、下限値Ｌｍｉｎ＝ＷＬ－ＷＷ／２に変換し、以下のように断層画像群Ｄｏを８ビットに変換する。

（式４９）
Ｄｏ（ｘ、ｙ、ｚ）＝（Ｄｏ（ｘ、ｙ、ｚ）－Ｌｍｉｎ）・２５５／（Ｌｍａｘ－Ｌｍｉｎ）
Ｄｏ（ｘ、ｙ、ｚ）＜０の場合、Ｄｏ（ｘ、ｙ、ｚ）＝０とする。
Ｄｏ（ｘ、ｙ、ｚ）＞２５５の場合、Ｄｏ（ｘ、ｙ、ｚ）＝２５５とする。

続いて、制御部１１は、８ビットの断層画像群Ｄｏに基づいて、任意断面のＭＰＲ画像ＩｍａｇｅＭＰＲを生成する（ステップＳ２４２）。例えば、図４２に示すように、ＸＹ平面に平行な体軸断面（アキシャル（Ａｘｉａｌ）断面）の場合、指定されたスライス位置ｚにおけるＸＹ平面上のボクセルに対して（Ｓｃｘ、Ｓｃｙ）変倍を加えながら各画素（ｘ、ｙ）を取得して２次元画像を再構成する。ＸＺ平面に平行な冠状断面（コロナル（Ｃｏｒｏｎａｌ）断面）の場合、指定されたスライス位置ｙにおけるＸＺ平面上のボクセルに対して（Ｓｃｘ、Ｓｃｚ）変倍を加えながら各画素（ｘ、ｚ）を取得して２次元画像を再構成する。ＺＹ平面に平行な矢状断面（サジタル（Ｓａｇｉｔｔａｌ）断面）の場合、指定されたスライス位置ｘにおけるＺＹ平面上のボクセルに対して（Ｓｃｚ、Ｓｃｙ）変倍を加えながら各画素（ｚ、ｙ）を取得して２次元画像を再構成する。このようにして３方向のＭＰＲ画像を生成する。

また、ＸＹＺ空間を斜めに切断した斜断面（オブリーク（Ｏｂｌｉｑｕｅ）断面）のＭＰＲ画像を生成する場合は、３方向のＭＰＲ画像のいずれかを２次元上で回転させながら３次元空間における回転を指示し、生成される回転行列に基づいて前記８ビットの断層画像群Ｄｏに対して３次元空間で回転を加えた上で、各２次元画像の再構成を行う。更に、スラブ厚を指定して体軸断面（アキシャル（Ａｘｉａｌ）断面）のＭＰＲ像を生成する場合は、指定されたスライス位置ｚの近傍で指定された厚み（スライス枚数）だけスライス位置を微小に移動させて、ＸＹ平面上のボクセルに対して（Ｓｃｘ、Ｓｃｙ）変倍を加えながら各画素（ｘ、ｙ）を複数通り取得して、各画素が複数個の信号値の代表信号値で構成される２次元画像を再構成する。複数個の信号値から代表信号値を算出する方法として、前述のＭＩＰ像を生成する場合と同様に、最大値（ＭＩＰ）、最小値（ＭｉｎＩＰ）、平均値（ＲａｙＳｕｍ）のいずれかを設定できる。

そして、制御部１１は、ステップＳ２４２において生成されたＭＰＲ画像ＩｍａｇｅＭＰＲにマスクデータＭａｓｋに基づいて生成されたマスク画像を合成（重畳）して出力する（ステップＳ２４３）。例えば、制御部１１は、ＭＰＲ画像ＩｍａｇｅＭＰＲに合成（重畳）するマスク画像として、Ｍａｓｋ（ｘ、ｙ、ｚ）＝０（非描画対象）の画素またはＭａｓｋ（ｘ、ｙ、ｚ）＝１（描画対象）の画素を所定の色で塗りつぶした画像、或いはＭａｓｋ（ｘ、ｙ、ｚ）＝０（非描画対象）の画素とＭａｓｋ（ｘ、ｙ、ｚ）＝１（描画対象）の画素を異なる色で色分けした画像を生成し、ＭＰＲ画像に合成する。

以上、図４１、図４２を参照しながら、ＭＰＲ画像を生成する処理を説明した。なお、上記の例では、ＭＰＲ画像にマスクデータＭａｓｋに基づいて生成したマスク画像を合成するようにしたが、ボリュームレンダリング画像やＭＩＰ画像の場合と同様に、マスクデータＭａｓｋにより規定される描画対象の画素のみを表示したＭＰＲ画像を生成してもよい。

以上、マスクデータＭａｓｋの適用例として、マスクデータＭａｓｋを用いて３次元再構成像（ボリュームレンダリング画像／ＭIＰ画像／ＭＰＲ画像）を生成する例を説明した。

（１４．３次元再構成像の表示例）
最後に、本開示の３次元再構成像生成装置２により生成された３次元再構成像（ボリュームレンダリング画像、ＭＩＰ画像、ＭＰＲ画像）の表示例を示す。

図４３は、本開示の３次元再構成像生成装置２により生成された胸腹部のボリュームレンダリング画像の表示例を示す図である。
図４３（ａ）は正面から観察したボリュームレンダリング画像であり、図４３（ｂ）は左側面から観察したボリュームレンダリング画像である。
肋骨、胸骨を含む全ての骨領域が十分に除去され、肺・心臓・血管・肝臓・腎臓などの診断上重要な内臓領域等が鮮明に可視化されていることが分かる。

図４４は、本開示の３次元再構成像生成装置２により生成された胸腹部のＭＩＰ画像の表示例を示す図である。
図４４（ａ）は正面から観察したＭＩＰ画像であり、図４５（ｂ）は左側面から観察したＭＩＰ画像である。
図４３のボリュームレンダリング画像と同様に、肋骨、胸骨を含む全ての骨領域が十分に除去され、肺・心臓・血管・肝臓・腎臓などの診断上重要な内臓領域等が鮮明に可視化されていることが確認できる。

図４５は、本開示の３次元再構成像生成装置２により生成された胸腹部のＭＰＲ画像の表示例を示す図である。図に示すように、Ａｘｉａｌ（体軸断面）、Ｃｏｒｏｎａｌ（冠状断面）、Ｓａｇｉｔｔａｌ（矢状断面）の各断面についてマスク画像が合成（重畳）されたＭＰＲ画像が表示される。図のマスク画像は、非描画対象を「赤色」（図はグレースケール画像のため「グレー」で表示されている）で塗りつぶした画像である。このように、マスク画像がＭＰＲ画像に合成表示されるので、ユーザは断層画像中の描画・非描画領域を容易に把握することができる。

比較例として、図４６に、本発明者が既に提案している特願２０１８－１２４１２５における不透明度の制御手法を用いて生成した胸腹部のボリュームレンダリング画像の表示例を示しておく。図４４（ａ）は正面から観察したボリュームレンダリング画像であり、図４４（ｂ）は左側面から観察したボリュームレンダリング画像である。特願２０１８－１２４１２５は、被写体領域の外郭を幾何形状で近似し、当該幾何形状を規定するパラメータから各画素の補正倍率を求め、ボリュームレンダリング（レイキャスティング処理）の過程で、各画素の信号値をオパシティカーブに通して得られる各画素の不透明度に対して前記補正倍率を乗算することで、不透明度を制御する方法である。

本開示の方法により得られるボリュームレンダリング画像（図４３）のほうが、骨領域が十分に除去され、骨領域の内側の内臓領域等が鮮明に可視化されていることが確認できる。
特願２０１８－１２４１２５の方法は、ユーザが設定するオパシティカーブに補正を施すため、オパシティカーブが変更されると、ＣＴ画像の変更と同様に骨除去の形態が変わり、パラメータ調整をやり直す必要があった。また、オパシティカーブが変更されると、観察対象の臓器の不透明度に補正倍率が乗算され、ユーザの意図とは異なって臓器が半透明（不鮮明）になることがあり、このように半透明になるとレイキャスティング法では処理負担の増大につながるという問題があった。

これに対し、本開示の方法では、ボリュームレンダリングの過程で不透明度を制御する方法をとらず、重み値（特願２０１８－１２４１２５の補正倍率に相当する値）に基づいて直接マスクデータを生成し、一般的なマスク参照によるレンダリング処理を実行するので、オパシティカーブが変更されることによるパラメータ調整の必要はなくなり、また、ユーザの意図とは異なって臓器が半透明（不鮮明）になることがなくなり、さらに、ボリュームレンダリングの過程で補正倍率を乗算する必要がないため処理負担も軽減される。

また、特願２０１８－１２４１２５の方法は、カラーマップ（オパシティカーブ）に補正を施すため、カラーマップを使用しないＭＩＰ画像やＭＰＲ画像には適用できなかった。一方、本開示の方法では、前記したように、ボリュームレンダリングの過程で不透明度を制御する方法をとらず、重み値（特願２０１８－１２４１２５の補正倍率に相当する値）に基づいて直接マスクデータを生成し、一般的なマスク参照による３次元再構成像生成処理を実行するため、カラーマップを使用しないＭＩＰ画像やＭＰＲ画像の生成処理にも適用可能となる。

また、特願２０１８－１２４１２５の方法では、信号値が所定の閾値より大きい全ての画素領域に基づいて被写体領域（図５の矩形領域４１）を算出していたため、被写体領域に寝台・ヘッドレストなどが含まれる場合があった。これに対し、本開示の方法では、信号値分布を考慮して被写体領域（図５の矩形領域４２）を算出するので、寝台・ヘッドレストなどを含まない領域を被写体領域として高精度に抽出することができる。

また特願２０１８－１２４１２５の方法では、３つの楕円で被写体領域の外郭を近似するが、椎骨を除去するための中央の楕円の径が小さめに設定されるため、近傍の大動脈などが同時に削られやすいという問題があった。また、設定パラメータが増え、パラメータ調整に手間がかかり、種々のＣＴ画像に適合させるのが困難な場合があった。これに対し、本開示の方法では、被写体領域の外郭を、胸骨側と椎骨側で楕円の径を変えることが可能な１つの変形楕円で近似することで、設定パラメータ数が抑えられるとともに、肋骨、胸骨を含む全ての骨領域を自動除去し、肺、心臓、血管、肝臓、腎臓などの診断上重要な内臓領域等を高精度に自動抽出することが可能となる。

以上、添付図面を参照しながら、本開示に係る好適な実施形態について説明したが、本開示はかかる例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本開示の技術的範囲に属するものと了解される。

なお、断層画像から被写体領域を算出する処理を実行する被写体領域算出装置を構成することもできる。この被写体領域算出装置は、図１の画像取得部２１、被写体領域算出部３１、被写体領域補正部３２を備え、図２と同様のハードウェア構成で実現される。そして、被写体領域算出装置の制御部が、図３のステップＳ１、図４のステップＳ２１、ステップＳ２２の処理を実行することで、断層画像から被写体領域を算出する。

１・・・・・・マスク作成装置
２・・・・・・３次元再構成像生成装置
２１・・・・・画像取得部
２２・・・・・描画重みテーブル作成部
２３・・・・・マスクデータ生成部
２４・・・・・クリッピング領域設定部
２５・・・・・有効ボクセル領域算出部
２６・・・・・３次元再構成像生成部
２７・・・・・３次元再構成像表示部
３１・・・・・被写体領域算出部
３２・・・・・被写体領域補正部
３３・・・・・被写体形状算出部
３４・・・・・重み値算出部
４１・・・・・矩形領域（第１の矩形領域）
４２・・・・・矩形領域（第２の矩形領域）
５０・・・・・被写体候補領域
５１・・・・・ボリュームレンダリング画像生成部
５２・・・・・ＭＩＰ画像生成部
５３・・・・・ＭＰＲ画像生成部
５１０・・・・カラーマップ設定部

Claims (23)

1. 複数の断層画像を取得する画像取得手段と、
   断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、
   断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、
   断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段と、
   断層画像毎に、断層画像の各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、
   を備えるマスク生成装置。
2. 前記重み値算出手段は、前記重み値を、前記幾何形状の中心が最大となり、前記幾何形状の外部が最小となり、前記幾何形状の内部は前記中心から遠ざかるにつれ小さくなるように算出する、
   請求項１に記載のマスク生成装置。
3. 前記重み値算出手段は、断層画像のスライス位置が中央から末端に位置するほど当該断層画像に対応する前記重み値が小さくなるように補正する、
   請求項１に記載のマスク生成装置。
4. 前記マスクデータ生成手段は、断層画像毎に、前記重み値により重み付けされた各画素の信号値を所定の閾値で二値化したマスク値を算出し、各画素のマスク値を保持するマスクデータを生成する、
   請求項１に記載のマスク生成装置。
5. 複数の断層画像を取得する画像取得手段と、
   断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、
   断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、
   断層画像毎に、断層画像の各画素が前記幾何形状の内部か外部のいずれに含まれるかを定義したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、を備え、
   前記被写体領域算出手段は、前記被写体領域として、信号値が所定の閾値より大きい画素が存在する領域の最外側を囲う矩形領域である第１の矩形領域を算出する、
   マスク生成装置。
6. 複数の断層画像を取得する画像取得手段と、
   断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、
   断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、
   断層画像毎に、断層画像の各画素が前記幾何形状の内部か外部のいずれに含まれるかを定義したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、を備え、
   前記被写体領域算出手段は、前記被写体領域として、信号値が所定の閾値より大きい画素が存在する領域の最外側を囲う矩形領域である第１の矩形領域の内部に定義される第２の矩形領域を算出する、
   マスク生成装置。
7. 前記被写体領域算出手段は、前記第２の矩形領域の中心座標が、前記第１の矩形領域内の信号値が前記所定の閾値より大きい画素の信号値で重み付けした重心座標となるように、前記第２の矩形領域を算出する、
   請求項６に記載のマスク生成装置。
8. 前記被写体領域算出手段は、
   前記第２の矩形領域の中心座標を起点に上下左右の４方向に前記第１の矩形領域の範囲で信号値が前記所定の閾値より大きい画素の信号値で重み付けした重心座標である第２重心を４箇所算出し、
   算出された左方向の第２重心と右方向の第２重心に基づいて前記第２の矩形領域の横方向のサイズを算出し、
   算出された上方向の第２重心と下方向の第２重心に基づいて前記第２の矩形領域の縦方向のサイズを算出する、
   請求項７に記載のマスク生成装置。
9. 前記被写体領域算出手段は、
   前記算出した４個の第２重心を起点に更に上下左右の４方向に前記第１の矩形領域の範囲で信号値が前記所定の閾値より大きい画素の信号値で重み付けした重心座標である第３重心を４箇所算出し、
   算出された左方向の第３重心と右方向の第３重心に基づいて前記第２の矩形領域の横方向のサイズを算出し、
   算出された上方向の第３重心と下方向の第３重心に基づいて前記第２の矩形領域の縦方向のサイズを算出する、
   請求項８に記載のマスク生成装置。
10. 前記被写体領域算出手段は、
    算出された左方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第１の横方向サイズを算出し、
    算出された右方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第２の横方向サイズを算出し、
    算出された上方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第１の縦方向サイズを算出し、
    算出された下方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第２の縦方向サイズを算出する、
    請求項８または請求項９に記載のマスク生成装置。
11. 前記断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段と、を更に備え、
    前記マスクデータ生成手段は、前記断層画像毎に、各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成し、
    前記幾何形状は、中心座標、横方向のサイズ、および縦方向のサイズによって規定される幾何形状であり、
    前記被写体形状算出手段は、
    前記幾何形状の中心座標として、前記第１または前記第２の矩形領域の中心座標を算出し、
    前記幾何形状の横方向のサイズとして、前記第１または前記第２の矩形領域の横方向のサイズを算出し、
    前記幾何形状の縦方向のサイズとして、前記第１または前記第２の矩形領域の縦方向のサイズを算出し、
    前記重み値算出手段は、算出した前記幾何形状の中心座標、横方向のサイズ、および縦方向のサイズに基づいて、前記重み値を算出する、
    請求項５から請求項９のいずれかに記載のマスク生成装置。
12. 前記被写体形状算出手段は、算出した前記幾何形状の横方向のサイズおよび縦方向のサイズに所定の横方向のサイズに対する倍率および縦方向のサイズに対する倍率を乗算して前記横方向のサイズおよび前記縦方向のサイズを補正する、
    請求項１１に記載のマスク生成装置。
13. 前記断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段と、を更に備え、
    前記マスクデータ生成手段は、前記断層画像毎に、各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成し、
    前記幾何形状は、中心座標、中心座標から左端までの距離である第１の横方向のサイズ、中心座標から右端までの距離である第２の横方向のサイズ、中心座標から上端までの距離である第１の縦方向のサイズ、および中心座標から下端までの距離である第２の縦方向のサイズによって規定される幾何形状であり、
    前記被写体形状算出手段は、
    前記幾何形状の中心座標として、前記第２の矩形領域の中心座標を算出し、
    前記幾何形状の第１の横方向のサイズとして、前記第２の矩形領域の第１の横方向のサイズを算出し、
    前記幾何形状の第２の横方向のサイズとして、前記第２の矩形領域の第２の横方向のサイズを算出し、
    前記幾何形状の第１の縦方向のサイズとして、前記第２の矩形領域の第１の縦方向のサイズを算出し、
    前記幾何形状の第２の縦方向のサイズとして、前記第２の矩形領域の第２の縦方向のサイズを算出し、
    前記重み値算出手段は、算出した前記幾何形状の中心座標、第１の横方向サイズ、第２の横方向サイズ、第１の縦方向サイズ、および第２の縦方向サイズに基づいて、前記重み値を算出する、
    請求項１０に記載のマスク生成装置。
14. 前記被写体形状算出手段は、
    算出した前記幾何形状の第１の横方向のサイズおよび第２の横方向のサイズに互いに異なる横方向のサイズに対する倍率を乗算して前記第１の横方向のサイズおよび前記第２の横方向のサイズを補正し、
    算出した前記幾何形状の第１の縦方向のサイズおよび第２の縦方向のサイズに互いに異なる縦方向のサイズに対する倍率を乗算して前記第１の縦方向のサイズおよび前記第２の縦方向のサイズを補正する、
    請求項１３に記載のマスク生成装置。
15. 複数の断層画像を取得する画像取得手段と、
    断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、
    断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、
    断層画像毎に、断層画像の各画素が前記幾何形状の内部か外部のいずれに含まれるかを定義したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、
    断層画像毎に、前記被写体領域算出手段により算出された被写体領域の縦方向のサイズと横方向のサイズの比率に応じた補正方法により該被写体領域を補正する被写体領域補正手段と、
    を備えるマスク生成装置。
16. 前記被写体領域算出手段は、前記被写体領域として、信号値が所定の閾値より大きい画素が存在する領域の最外側を囲う矩形領域である第１の矩形領域の内部に定義される第２の矩形領域を算出し、
    前記被写体領域補正手段は、前記被写体領域算出手段により算出された前記第２の矩形領域の縦方向のサイズを横方向のサイズで除算した縦横比率を算出し、
    算出された前記縦横比率が所定の閾値より小さい場合、前記被写体領域算出手段により算出された前記第２の矩形領域の横方向のサイズおよび縦方向のサイズの各々に前記縦横比率を乗算して補正し、
    算出された前記縦横比率が前記所定の閾値以上の場合、前記第１の矩形領域の内部で信号値が前記所定の閾値より大きい画素の個数を前記第２の矩形領域の面積で除算した面積比率を算出し、前記被写体領域算出手段により算出された前記第２の矩形領域の横方向のサイズおよび縦方向のサイズの各々に前記面積比率を乗算して補正する、
    請求項１５に記載のマスク生成装置。
17. 前記被写体領域算出手段は、前記被写体領域として、信号値が所定の閾値より大きい画素が存在する領域の最外側を囲う矩形領域である第１の矩形領域の内部に定義される第２の矩形領域を算出し、
    前記被写体領域算出手段は、前記第２の矩形領域の中心座標が、前記第１の矩形領域内の信号値が前記所定の閾値より大きい画素の信号値で重み付けした重心座標となるように、前記第２の矩形領域を算出し、
    前記被写体領域算出手段は、
    前記第２の矩形領域の中心座標を起点に上下左右の４方向に前記第１の矩形領域の範囲で信号値が前記所定の閾値より大きい画素の信号値で重み付けした重心座標である第２重心を４箇所算出し、
    算出された左方向の第２重心と右方向の第２重心に基づいて前記第２の矩形領域の横方向のサイズを算出し、
    算出された上方向の第２重心と下方向の第２重心に基づいて前記第２の矩形領域の縦方向のサイズを算出し、
    前記被写体領域算出手段は、
    前記算出した４個の第２重心を起点に更に上下左右の４方向に前記第１の矩形領域の範囲で信号値が前記所定の閾値より大きい画素の信号値で重み付けした重心座標である第３重心を４箇所算出し、
    算出された左方向の第３重心と右方向の第３重心に基づいて前記第２の矩形領域の横方向のサイズを算出し、
    算出された上方向の第３重心と下方向の第３重心に基づいて前記第２の矩形領域の縦方向のサイズを算出し、
    前記被写体領域算出手段は、
    算出された左方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第１の横方向サイズを算出し、
    算出された右方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第２の横方向サイズを算出し、
    算出された上方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第１の縦方向サイズを算出し、
    算出された下方向の第２重心または第３重心から前記第２の矩形領域の中心座標までの距離である第２の縦方向サイズを算出し、
    前記被写体領域補正手段は、前記被写体領域算出手段により算出された前記第２の矩形領域の縦方向のサイズを横方向のサイズで除算した縦横比率を算出し、
    算出された前記縦横比率が所定の閾値より小さい場合、前記被写体形状算出手段により算出された前記第２の矩形領域の第１の横方向のサイズ、第２の横方向のサイズ、第１の縦方向のサイズ、および第２の縦方向のサイズの各々に前記縦横比率を乗算して補正し、
    算出された前記縦横比率が前記所定の閾値以上の場合、前記第１の矩形領域の内部で信号値が前記所定の閾値より大きい画素の個数を前記第２の矩形領域の面積で除算した面積比率を算出し、前記被写体領域算出手段により算出された前記第２の矩形領域の第１の横方向のサイズ、第２の横方向のサイズ、第１の縦方向のサイズ、および第２の縦方向のサイズの各々に前記面積比率を乗算して補正する、
    請求項１５に記載のマスク生成装置。
18. 複数の断層画像を取得する画像取得手段と、
    断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段と、
    断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段と、
    断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段と、
    断層画像毎に、断層画像の各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段と、
    前記断層画像の各画素をボクセル空間に配置して構成される３次元ボクセルと前記マスクデータに基づいて３次元再構成像を生成する３次元再構成像生成手段と、
    を備える３次元再構成像生成装置。
19. 前記３次元再構成像生成手段は、前記マスクデータに基づいて各ボクセルの不透明度を計算に用いるか否かを設定する、
    請求項１８に記載の３次元再構成像生成装置。
20. 前記３次元再構成像生成手段は、前記マスクデータに基づいて各ボクセルの信号値を計算に用いるか否かを設定する、
    請求項１８に記載の３次元再構成像生成装置。
21. 前記３次元再構成像生成手段は、前記３次元再構成像として、前記断層画像に基づいて生成された３次元再構成像に前記マスクデータに基づいて生成されたマスク画像を合成した画像を生成する、
    請求項１８に記載の３次元再構成像生成装置。
22. コンピュータが、
    複数の断層画像を取得する画像取得ステップと、
    断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出ステップと、
    断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出ステップと、
    断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出ステップと、
    断層画像毎に、断層画像の各画素に対して重み値算出ステップにより算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成ステップと、
    を実行するマスク生成方法。
23. コンピュータを、
    複数の断層画像を取得する画像取得手段、
    断層画像毎に、各画素の信号値に基づいて被写体領域を算出する被写体領域算出手段、
    断層画像毎に、算出された被写体領域の外郭を所定の幾何形状で近似し、当該幾何形状のパラメータを算出する被写体形状算出手段、
    断層画像毎に、前記幾何形状のパラメータに基づいて断層画像の各画素を描画対象とする重み値を算出する重み値算出手段、
    断層画像毎に、断層画像の各画素に対して重み値算出手段により算出された重み値を所定の閾値で二値化したマスク値を算出し、各画素に対するマスク値を保持するマスクデータを生成するマスクデータ生成手段、
    として機能させるプログラム。

Images