(57)【特許請求の範囲】
【請求項１】
標的質量が受ける放射線量を計算する方法であって、該方法は、
前記の標的質量の第１部位を囲むボクセル容積を入力する工程と、
前記の各ボクセル容積内の標的質量の密度に応じたボクセル質量密度を入力する工程と、
前記の標的質量の第２部位を囲み、前記の第１部位とオーバラップするドーセル容積を定義する工程と、
前記のドーセル容積とボクセル容積との間の共通容積を算出する工程と、
増加ドーセル質量を得るため、前記の共通容積を対応する前記ボクセル質量密度で乗算する工程と、
ドーセル質量を得るため、前記ドーセル容積に応じた増加ドーセル質量を加算する工程とを備えるとともに、
前記各工程が終了した後に、
前記の標的質量におけるエネルギー堆積の位置を特定する工程と、
前記のドーセル容積内の標的質量が受ける放射線量を算出する工程を行うことを特徴とする放射線量の計算方法。
【請求項２】
前記の共通容積の算出工程が、
前記のドーセル容積がボクセル容積に対する周期的反復パターンを形成できるよう、ドーセル容積の大きさを縮小して全体数を増加させる工程と、
前記の周期的反復パターンの第１周期とボクセル容積との間の周期的共通容積を算出する工程と、
前記の周期的反復パターンの次の周期内の、次の共通容積の組を算出するため、前記の周期的共通容積を再利用する工程とを備えた、請求項１記載の放射線の計算方法。
【請求項３】
前記の特定工程が、
前記の標的質量から離脱する、入射非荷電粒子のトラッキングを終了する工程と、
前記の標的質量と相互作用するような、入射非荷電粒子のトラッキングを終了する工程と、
入射非荷電粒子の相互作用の位置を特定する工程と、
その位置でのエネルギー堆積を特定する工程とを備えた、請求項１記載の放射線量の計算方法。
【請求項４】
前記の特定工程が、
前記の標的質量から離脱する、前記のドーセル容積内で形成された娘核非荷電粒子のトラッキングを終了する工程と、
前記の標的質量と相互作用するような、前記のドーセル容積内で形成された娘核非荷電粒子のトラッキングを終了する工程と、
娘核非荷電粒子の相互作用の位置を特定する工程と、
その位置でのエネルギー堆積を特定する工程とを備えた、請求項１記載の放射線量の計算方法。
【請求項５】
前記の特定工程が、
前記の標的質量から離脱する、前記のドーセル容積外で形成され、かつ、それに向かって移動する娘核非荷電粒子のトラッキングを終了する工程と、
前記の標的質量と相互作用するような、前記のドーセル容積外で形成され、かつ、それに向かって移動する娘核非荷電粒子のトラッキングを終了する工程と、
娘核非荷電粒子の相互作用の位置を特定する工程と、
その位置でのエネルギー堆積を特定する工程とを備えた、請求項１記載の放射線量の計算方法。
【請求項６】
前記の特定工程が、
分散減少法を活用できるような、前記のドーセル容積外で形成され、かつ、それから離反する娘核非荷電粒子をトラッキングする工程を備えた、請求項１記載の放射線量の計算方法。
【請求項７】
前記の特定工程が、
ロシアンルーレット法を活用できるような、前記のドーセル容積外で形成され、かつ、それから離反する娘核非荷電粒子をトラッキングする工程を備えた、請求項１記載の放射線量の計算方法。
【請求項８】
前記の特定工程が、
前記のドーセル容積を通過する荷電粒子のトラッキングをする工程と、
前記の各荷電粒子の経路に沿った１組のエネルギー堆積の位置を定義する工程と、
前記の１組のエネルギー堆積の位置にエネルギー堆積を割り当てる工程とを備えた、請求項１記載の放射線量の計算方法。
【請求項９】
前記の特定工程が、
前記のドーセル容積外の荷電粒子の位置を特定する工程と、
前記の荷電粒子に関連するエネルギー値の組である第１の組を算出する工程と、
前記の荷電粒子がドーセル容積に到達するのに要するエネルギー値の組である第２の組を算出する工程と、
前記の第１の組の荷電粒子に対応するエネルギー値が、前記の第２の組の荷電粒子に対応するエネルギー値よりも小さい場合には、前記荷電粒子中の各荷電粒子のトラッキングを終了する工程とを備えた、請求項１記載の放射線量の計算方法。
【請求項１０】
前記の特定工程が、
前記のドーセル容積外の荷電粒子の位置を特定する工程と、
前記の荷電粒子に関連するエネルギー値の組である第１の組を算出する工程と、
前記の荷電粒子がドーセル容積に到達するのに要するエネルギー値の組である第２の組を算出する工程と、
前記の第１の組の荷電粒子に対応するエネルギー値が、前記の第２の組の荷電粒子に対応するエネルギー値よりも小さくない場合には、前記荷電粒子中の各荷電粒子をトラッキングする工程とを備えた、請求項１記載の放射線量の計算方法。
【請求項１１】
前記の特定工程が、
前記のドーセル容積とオーバラップしない１組のボクセル容積を特定する工程と、
各荷電粒子が、前記の１組のボクセル容積の各ボクセル容積からドーセル容積へ移動するのに要するエネルギー値を算出する工程と、
前記のエネルギー値をボクセル容積に割り当てる工程とを備えた、請求項１記載の放射線量の計算方法。
【請求項１２】
標的質量が受ける放射線量を計算するシステムであって、該システムは、
前記の標的質量の第１部位を囲むボクセル容積を入力する手段と、
前記の各ボクセル容積内の標的質量の密度に応じたボクセル質量密度を入力する手段と、
前記の標的質量の第２部位を囲み、前記の第１部位とオーバラップするドーセル容積を定義する手段と、
前記のドーセル容積とボクセル容積との間の共通容積を算出する手段と、
増加ドーセル質量を得るため、前記の共通容積を対応する前記ボクセル質量密度で乗算する手段と、
ドーセル質量を得るため、前記ドーセル容積に応じた増加ドーセル質量を加算する手段とを備えるとともに、
前記各手段に加え、さらに、
前記の標的質量におけるエネルギー堆積の位置を特定する手段と、
前記のドーセル容積内の標的質量が受ける放射線量を算出する手段とを備えたことを特徴とする放射線量の計算システム。
【請求項１３】
前記の共通容積を算出する手段が、
前記のドーセル容積がボクセル容積に対する周期的反復パターンを形成できるよう、ドーセル容積の大きさを縮小して全体数を増加させる手段と、
前記の周期的反復パターンの第１周期とボクセル容積との間の周期的共通容積を算出する手段と、
前記の周期的反復パターンの次の周期内の、次の共通容積の組を算出するため、前記の周期的共通容積を再利用する手段とを備えた、請求項１２記載の放射線の計算システム。
【請求項１４】
前記の特定手段が、
前記の標的質量から離脱する、入射非荷電粒子のトラッキングを終了する手段と、
前記の標的質量と相互作用するような、入射非荷電粒子のトラッキングを終了する手段と、
入射非荷電粒子の相互作用の位置を特定する手段と、
その位置でのエネルギー堆積を特定する手段とを備えた、請求項１２記載の放射線量の計算システム。
【請求項１５】
前記の特定手段が、
前記の標的質量から離脱する、前記のドーセル容積内で形成された娘核非荷電粒子のトラッキングを終了する手段と、
前記の標的質量と相互作用するような、前記のドーセル容積内で形成された娘核非荷電粒子のトラッキングを終了する手段と、
娘核非荷電粒子の相互作用の位置を特定する手段と、
その位置でのエネルギー堆積を特定する手段とを備えた、請求項１２記載の放射線量の計算システム。
【請求項１６】
前記の特定手段が、
前記の標的質量から離脱する、前記のドーセル容積外で形成され、かつ、それに向かって移動する娘核非荷電粒子のトラッキングを終了する手段と、
前記の標的質量と相互作用するような、前記のドーセル容積外で形成され、かつ、それに向かって移動する娘核非荷電粒子のトラッキングを終了する手段と、
娘核非荷電粒子の相互作用の位置を特定する手段と、
その位置でのエネルギー堆積を特定する手段とを備えた、請求項１２記載の放射線量の計算システム。
【請求項１７】
前記の特定手段が、
分散減少法を活用できるような、前記のドーセル容積外で形成され、かつ、それから離反する娘核非荷電粒子をトラッキングする手段を備えた、請求項1２記載の放射線量の計算システム。
【請求項１８】
前記の特定手段が、
ロシアンルーレット法を活用できるような、前記のドーセル容積外で形成され、かつ、それから離反する娘核非荷電粒子をトラッキングする手段を備えた、請求項1２記載の放射線量の計算システム。
【請求項１９】
前記の特定手段が、
前記のドーセル容積を通過する荷電粒子のトラッキングをする手段と、
前記の各荷電粒子の経路に沿った１組のエネルギー堆積の位置を定義する手段と、
前記の１組のエネルギー堆積の位置にエネルギー堆積を割り当てる手段とを備えた、請求項１２記載の放射線量の計算システム。
【請求項２０】
前記の特定手段が、
前記のドーセル容積外の荷電粒子の位置を特定する手段と、
前記の荷電粒子に関連するエネルギー値の組である第１の組を算出する手段と、
前記の荷電粒子がドーセル容積に到達するのに要するエネルギー値の組である第２の組を算出する手段と、
前記の第１の組の荷電粒子に対応するエネルギー値が、前記の第２の組の荷電粒子に対応するエネルギー値よりも小さい場合には、前記荷電粒子中の各荷電粒子のトラッキングを終了する手段とを備えた、請求項１２記載の放射線量の計算システム。
【請求項２１】
前記の特定手段が、
前記のドーセル容積外の荷電粒子の位置を特定する手段と、
前記の荷電粒子に関連するエネルギー値の組である第１の組を算出する手段と、
前記の荷電粒子がドーセル容積に到達するのに要するエネルギー値の組である第２の組を算出する手段と、
前記の第１の組の荷電粒子に対応するエネルギー値が、前記の第２の組の荷電粒子に対応するエネルギー値よりも小さくない場合には、前記荷電粒子中の各荷電粒子をトラッキングする手段とを備えた、請求項１２記載の放射線量の計算システム。
【請求項２２】
前記の特定手段が、
前記のドーセル容積とオーバラップしない１組のボクセル容積を特定する手段と、
各荷電粒子が、前記の１組のボクセル容積の各ボクセル容積からドーセル容積へ移動するのに要するエネルギー値を算出する手段と、
前記のエネルギー値をボクセル容積に割り当てる手段とを備えた、請求項１２記載の放射線量の計算システム。
【請求項２３】
コンピュータのプログラムコードを収録するコンピュータ利用可能媒体であって、
標的質量の第１部位を囲むボクセル容積を入力する工程と、
前記の各ボクセル容積内の標的質量の密度に応じたボクセル質量密度を入力する工程と、
前記の標的質量の第２部位を囲み、前記の第１部位とオーバラップするドーセル容積を定義する工程と、
前記のドーセル容積とボクセル容積との間の共通容積を算出する工程と、
増加ドーセル質量を得るため、前記の共通容積を対応する前記ボクセル質量密度で乗算する工程と、
ドーセル質量を得るため、前記ドーセル容積に応じた増加ドーセル質量を加算する工程とを備えるとともに、
前記各工程が終了した後に、
前記の標的質量におけるエネルギー堆積の位置を特定する工程と、
前記のドーセル容積内の標的質量が受ける放射線量を算出する工程を実行することにより、標的質量が受ける放射線量をコンピュータに計算させるための、コンピュータのプログラムコードを収録するコンピュータ利用可能媒体。
【請求項２４】
前記の共通容積の算出工程が、
前記のドーセル容積がボクセル容積に対する周期的反復パターンを形成できるよう、ドーセル容積の大きさを縮小して全体数を増加させる工程と、
前記の周期的反復パターンの第1周期とボクセル容積との間の周期的共通容積を算出する工程と、
前記の周期的反復パターンの次の周期内の次の共通容積の組を算出するため、前記の周期的共通容積を再利用する工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項２５】
前記の特定工程が、
前記の標的質量から離脱する、入射非荷電粒子のトラッキングを終了する工程と、
前記の標的質量と相互作用するような、入射非荷電粒子のトラッキングを終了する工程と、
入射非荷電粒子の相互作用の位置を特定する工程と、
その位置でのエネルギー堆積を特定する工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項２６】
前記の特定工程が、
前記の標的質量から離脱する、前記のドーセル容積内で形成された娘核非荷電粒子のトラッキングを終了する工程と、
前記の標的質量と相互作用するような、前記のドーセル容積内で形成された娘核非荷電粒子のトラッキングを終了する工程と、
娘核非荷電粒子の相互作用の位置を特定する工程と、
その位置でのエネルギー堆積を特定する工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項２７】
前記の特定工程が、
前記の標的質量から離脱する、前記のドーセル容積外で形成され、かつ、それに向かって移動する娘核非荷電粒子のトラッキングを終了する工程と、
前記の標的質量と相互作用するような、前記のドーセル容積外で形成され、かつ、それに向かって移動する娘核非荷電粒子のトラッキングを終了する工程と、
娘核非荷電粒子の相互作用の位置を特定する工程と、
その位置でのエネルギー堆積を特定する工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項２８】
前記の特定工程が、
分散減少法を活用できるような、前記のドーセル容積外で形成され、かつ、それから離反する娘核非荷電粒子をトラッキングする工程を備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項２９】
前記の特定工程が、
ロシアンルーレット法を活用できるような、前記のドーセル容積外で形成され、かつ、それから離反する娘核非荷電粒子をトラッキングする工程を備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項３０】
前記の特定工程が、
前記のドーセル容積を通過する荷電粒子のトラッキングをする工程と、
前記の各荷電粒子の経路に沿った１組のエネルギー堆積の位置を定義する工程と、
前記の１組のエネルギー堆積の位置にエネルギー堆積を割り当てる工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項３１】
前記の特定工程が、
前記のドーセル容積外の荷電粒子の位置を特定する工程と、
前記の荷電粒子に関連するエネルギー値の組である第１の組を算出する工程と、
前記の荷電粒子がドーセル容積に到達するのに要するエネルギー値の組である第２の組を算出する工程と、
前記の第１の組の荷電粒子に対応するエネルギー値が、前記の第２の組の荷電粒子に対応するエネルギー値よりも小さい場合には、前記荷電粒子中の各荷電粒子のトラッキングを終了する工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項３２】
前記の特定工程が、
前記のドーセル容積外の荷電粒子の位置を特定する工程と、
前記の荷電粒子に関連するエネルギー値の組である第１の組を算出する工程と、
前記の荷電粒子がドーセル容積に到達するのに要するエネルギー値の組である第２の組を算出する工程と、
前記の第１の組の荷電粒子に対応するエネルギー値が、前記の第２の組の荷電粒子に対応するエネルギー値よりも小さくない場合には、前記荷電粒子中の各荷電粒子をトラッキングする工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項３３】
前記の特定工程が、
前記のドーセル容積とオーバラップしない１組のボクセル容積を特定する工程と、
各荷電粒子が、前記の１組のボクセル容積の各ボクセル容積からドーセル容積へ移動するのに要するエネルギー値を算出する工程と、
前記のエネルギー値をボクセル容積に割り当てる工程とを備えた、請求項２３記載のコンピュータ利用可能媒体。
【請求項３４】
標的質量が受ける放射線量を計算するシステムであって、
前記の標的質量の第１部位を囲むボクセル容積を入力し、前記の標的質量の第２部位を囲み、前記の第１部位とオーバラップするドーセル容積を定義し、前記のドーセル容積とボクセル容積との間の共通容積を算出するための共通容積算出モジュールと、
増加ドーセル質量を得るため、前記の共通容積を、ボクセル質量密度で乗算し、かつ、ドーセル質量を得るため、ドーセル容積に応じた増加ドーセル質量を加算するための、前記の各ボクセル容積内の標的質量の密度に応じたボクセル質量密度を入力するためのドーセル質量モジュールと、
前記の標的質量におけるエネルギー堆積の位置を特定するための、前記の共通容積算出モジュールに接続された放射線トランスポートモジュールと、
前記のドーセル容積内の標的質量が受ける放射線量を算出するための、前記の共通容積算出モジュールと前記の放射線トランスポートモジュールと前記ドーセル質量モジュールとに接続された線量計算モジュールとを備えた放射線量の計算システム。
【発明の詳細な説明】
【０００１】
合衆国政府は、ローレンスリバモア国立研究所の研究に関する合衆国エネルギー省とカリフォルニア大学間で結ばれた契約番号Ｗ−７４０５−ＥＮＧ−４８により、本発明における権利を有する。
【０００２】
この米国特許出願の優先日は、バーグストロム他による1998年１１月１１日出願の、先行の米国仮特許出願番号60／108092の「モンテカルロ法に基づく粒子トランスポートグリッド（Monte  Carlo  Based  Particle  Transport  Grid）のサブ容積（Sub-Volumes）内の放射線量計算のためのシステムと方法」により確定している。
【０００３】
本出願は、チャンドラー他による１９９６年３月５日出願の、参考同時係属米国特許出願番号08／610,917、「全粒子モンテカルロトランスポート（All  Particle  Monte  Carlo  Transport）を使用した、放射線治療線量の計算」に関連し、一部包含するものである。この関連出願は、通例としてカリフォルニア大学理事会に譲渡されている。
【０００４】
（技術分野）
本発明は、放射線量計算のシステムと方法に関するものであり、特に、モンテカルロ法に基づいた粒子トランスポートグリッドのサブ容積内における放射線量の計算方法に関する。
【０００５】
（技術背景）
人間に対する有効的な放射線療法とは、患者の病状診断、治療計画、療法投与という３つの関連する要素から成り立つ。過去数十年の間に、患者診断や療法投与の技術は大きく改善されてきているが、現状において、（ＣＴやＭＲＩ走査法などの）撮影や治療の最新装置を、放射線治療分野で最大限に利用するには問題がある。それゆえ、それらの装置利用の実現には、高い見識と正確な治療計画が必要である。
【０００６】
治療計画において重要な点は、患者へ照射する放射線量を的確に判断することにある。放射線単位と測定法に関する国際委員会（ＩＣＲＵ）の報告書「高エネルギーフォトンとエレクトロンによる、外部ビーム放射線治療手順におけるコンピュータの使用」、ＩＣＲＵレポート４２（１９８７）では、低線量の勾配域（low  dose  gradient  regions）における線量計算の目標精度を２％と推奨している。しかしながら、現状での不均一媒質（heterogeneous  media）においての計算エラーは３％〜１０％以上であることが、Ｍ．Ｅ．マスターソン他による「外部フォトン線量計算の検証における学会内部実験」、国際腫瘍学生理学物理学連携学会誌、21  37‐58（１９９１）、Ｋ．アヤンガー他による、「特別設計の不均質ファントムを使った３次元線量計算アルゴリズムの実験的検証」、医学物理学学会誌、21(2)325‐329（１９９３）、Ｃ．ハークマンス他による「頭部と頚部におけるフォトン線量計算でのペンシルビームの制限」、放射線療法腫瘍学学会誌、39  74‐80（１９９５）、Ｔ．クノースらによる、「肺組織におけるフォトン線量計算でのペンシルビームの制限」、物理学医学生物学学会誌、40  1411‐1420（１９９５）で、述べられている。
【０００７】
畳み込み（convolution）計算と重ね合わせ（superposition）計算のアルゴリズムは、従来方法に比べるとかなり正確であるが、計算に必要な消費時間が大幅に増加する。畳み込み計算と重ね合わせ計算のアルゴリズムについては、A．ボーイアとE．モックによる、「畳み込み計算を採用したフォトン線量分配モデル」医学物理学学会誌、12  169‐177（１９８５）、Ｔ．R．マッキー他による、「１５MV  Ｘ線の線量を計算するための畳み込み計算法」医学物理学学会誌、12  188‐169（１９８５）、Ｒ．モハン他による「フォトンのための差動ペンシルビーム線量計算モデル」、医学物理学学会誌、13  64‐73（１９８６）、Ａ．アネスヨ他による「高エネルギーフォトンビームによる治療計画のための点分散関数（point  spread  functions）の計算と応用」、編纂腫瘍学会誌26  49‐56（１９８７）、Ｒ．モハンとＣ．チュイによる、「３次元治療計画のための、不規則形状領域の線量分布計算における高速フーリエ変換法の利用」医学物理学学会誌、14  70‐77（１９８７）、Ｐ．Ｅ．メカルフ他による、「重ね合わせ計算による多色高エネルギーフォトンビームのモデル化」、オーストリア物理化学医学会誌、12  138‐148（１９８９）、Ｐ．Ｗ．ホーバン他による、「肺における１０ＭＶのフォトンの重ね合せ計算」、オーストリア物理化学医学会誌13  81‐92（１９９０）に、述べてある。
【０００８】
モンテカルロ式線量計算技法は、しかしながら、線量計算法の好適な手段として広く受け入れられている。モンテカルロ技法は、ほとんど全ての環境において正確に線量を算出することができる。関連記事として、Ｒ．モーハンによる「放射線治療計画のための線量計算方法」（Ｔ．Ｍ．ジェンキンス、Ｗ．Ｒ．ネルソン、Ａ．リンディ編、「エレクトロンとフォトンのモンテカルロトランスポート」549‐471（１９８８）収録）、Ｄ．Ｗ．Ｏ．ロジャーズとＡ．Ｆ．ビェラジューによる「イオン化放射線量法のためのエレクトロンとフォトンのトランスポートにおけるモンテカルロ技法」（イオン化放射線線量計測、第3巻、Ｋ．Ｒ．カセ、Ｂ．Ｅ．ビヨルガード、Ｒ．Ｈ．アティックス編、アカデミックニューヨーク（１９９０）、pp．427‐539収録）、Ｔ．Ｒ．マキーによる「放射線療法におけるモンテカルロ方法の応用例」（『イオン化放射能の線量法』第3巻、Ｋ．Ｒ．カセ、Ｂ．Ｅ．ビヨルガード、Ｒ．Ｈ．アティクス編、アカデミックニューヨーク（１９９０）pp．541‐620収録）、があげられる。モンテカルロ技法は、また、食物照射のモデル化や理解に対する手段にも応用することも可能である。
【０００９】
残念ながら、単一治療計画にモンテカルロ計算法を適用するには、数日あるいは数週間要する場合がある。Ａ．Ｅ．ナーウンの「フォトンとエレクトロンのモンテカルロ法の概要」（Ｔ．Ｍ．ジェンキンス、Ｗ．Ｒ．ネルソン、Ａ．リンディ編、エレクトロンとフォトンのモンテカルロトランスポート、３-２０（１９８８）収録）を参照してほしい。モンテカルロ法のアルゴリズムは、現状では専門的なユーザによるセットアップが必要であり、その実行速度はかなり低速である。さらに、モンテカルロ法のアルゴリズムは、ＣＴやＭＲＩの走査全体における線量計算において利用可能である。このため、ＣＴやＭＲＩ走査の小規模における線量計算が必要な場合は、汎用モンテカルロ法アルゴリズムでは必要以上の情報が提供されるため、コストや計算時間が増える。
【００１０】
上記で述べた問題点の解決するには、従来技術における問題を克服するような放射線量計算法のシステムと方法が必要である。
【００１１】
（発明の開示）
本発明は、モンテカルロ粒子トランスポートグリッドのサブ容積内における放射線量の計算のシステムおよび方法に関する。本発明は、エネルギー値つまり線量の計算を必要とする多くの課題に適用可能である。本発明の方法においては、最初に、標的質量（target  mass）の第１部位を囲むボクセル容積（voxel  volumes）を入力する。次の第２工程では、標的質量の第２部位を囲み第１部位とオーバラップするドーセル容積（dosel  volumes）を定義する(define)。第３工程で、ドーセル容積とボクセル容積間の共通容積（common  volumes）を算出する。第４工程では、標的質量内でのエネルギー堆積（deposits）の位置を特定する(identify)。そして、第５工程で、ドーセル容積内における標的質量が受ける放射線量を算出するのである。
【００１２】
本発明のシステムにおいては、共通容積算出モジュールにより、標的質量の第１部位を囲むボクセル容積が入力され、各ボクセル容積内の標的質量の密度に対応したボクセル質量密度が入力され、標的質量の第２部位を囲み第１部位とオーバラップするドーセル容積が定義され、ドーセル容積とボクセル容積の間の共通容積が算出されるのである。ドーセル質量モジュールでは、増加ドーセル質量を計算するために共通容積を対応するボクセル質量密度で乗算し、ドーセル質量を計算するためドーセル容積に対応した前記増加ドーセル質量を加算する。放射線トランスポートモジュールは、標的質量内のエネルギー堆積の位置を特定する。そして、前記の共通容積算出モジュールと放射線トランスポートモジュールに接続された線量計算モジュールで、ドーセル容積内の標的質量が受ける放射線量が算出されるのである。
【００１３】
本発明のシステムや方法を使えば、放射線量計算のための汎用モンテカルロ法アルゴリズムを、ＣＴやＭＲＩ走査より小さなサブ容積における線量計算を実行する専用アルゴリズムに効率良く変換できる。ドーセル容積は、ユーザが所望する位置や解像度でのドーセル計算を可能にする。またドーセル容積を利用すれば、患者の（人体を横切る）横方向では精度が高いが、頭部から足部への方向においてはやや粗い精度で、一定でないようなＣＴやＭＲＩ走査でのスライス間隔に比べて、線量計算精度が高くなる。
【００１４】
ドーセル容積のサイズによっては、計算値はモンテカルロ法のエネルギー堆積値の低値部分に収束する。加えて、本発明は分散減少法を使って、ＣＴやＭＲＩ走査よりも小さな線量記録領域をユーザが定義できるため、現在の分析方法よりもかなり迅速に粒子をトラッキングすることができる。
本発明の上記の特徴やその他の特徴は、下記の詳細な説明、図面、請求の範囲などを参照すれば当業者には明白となろう。
【００１５】
（発明を実施するための最良の形態）
図１は、本発明によるモンテカルロ法に基づく粒子トランスポートグリッドのサブ容積内の放射線量の計算のためのシステムの図表１００である。放射線量とは、（図示しない）標的質量内の単位質量あたりのエネルギー量として定義される、周知の用語である。標的質量が受ける放射線量の計算は、３つの基本的な段階で行われる。まず第１の段階では、放射線量の計算の対象たる標的質量内のサブ容積を決定する。第２に、標的質量へ放射される放射線量をモニターして、サブ容積内に投下されるエネルギー量を記録する。第３に、サブ容積内に投下されたエネルギー量をサブ容積内での標的質量単位で割って、サブ容積単位の放射線量を算出するのである。
【００１６】
前記の図表１００には、標的質量とオーバラップするサブ容積の集合体１０２が２次元で図示されている。標的質量は、放射線量を受けるガン患者の腫瘍患部などの構造体である。図１では、簡略化のため前記のサブ容積が２次元で示されているが、実際の標的質量は３次元であるため、本発明では３次元で説明する。
【００１７】
前記のサブ容積の集合体１０２は、ボクセル（すなわち容積要素）グリッド１０４で構成されている。ボクセルやボクセルグリッドは、当業者にはよく知られた用語である。ボクセルグリッド１０４は、粒子トランスポートグリッドと呼ぶこともできる。ボクセルグリッド１０４は、ＣＴ法やＭＲＩ法で標的質量を走査する場合に使われる。ＣＴ法やＭＲＩ法では、ボクセルグリッド１０４内の直方体をボクセル１０６と呼ぶ。各ボクセル１０６は、ＣＴ法やＭＲＩ法で設定可能な一定の質量密度と原子構成を有する。
【００１８】
現状では、各毎ボクセルの放射線量を、モンテカルロ法に基づく放射線量計算システムによって算出する。モンテカルロ計算法では、計算電力（computing  power）が大量に必要になるため、特にユーザが小規模範囲における放射線量情報を必要とする場合、各毎ボクセルの放射線量を算出するのに過剰な時間がかかる。
【００１９】
これに対して本発明では、放射線量情報が必要な標的質量における被算出範囲をカバーするドーセルグリッド１０８を、ユーザが定義することができる。一般的に、前記のドーセルグリッド１０８で定義される標的質量の大きさは、前記のボクセルグリッド１０４で定義されたものよりも小さくなる。しかしながらこのドーセルグリッド１０８を拡大することは可能である。ドーセルグリッド１０８を縮小した場合は、モンテカルロ計算に必要な時間と電力消費は大幅に減少する。また、前記のドーセルグリッド１０８は前記のボクセルグリッド１０４から独立しているため、ユーザはモンテカルロ計算の対象範囲を指定する(specify)のに十分な柔軟性を有することになる。また、前記のドーセルグリッドとボクセルグリッド間の間隔の割合を、ユーザの操作によって変更できる。このため、ＣＴ法あるいはＭＲＩ法がひとたび開始されると、ボクセルグリッド１０４の解像度が固定されてしまうのに対して、ドーセルグリッド１０８によって、ユーザは多様な解像度の線量計算を実行することができる。
【００２０】
前記のドーセルグリッド１０８は、ドーセルの集合体１１０から構成される。各ドーセルは、ドーセルグリッド点１１１を中心とするドーセル容積をもつ。このドーセル１１０は、球体、立方体、あるいは四面体のなどのいかなる形状であっても構わない。また、ドーセル１１０の互いの重なり度合いも多様であってよい。この好適実施例では、ドーセル１１０は一定の大きさと間隔をもつ。前記の各ドーセル１１０の半径は、ユニットセル１１３に対して対角にあるドーセル１１０がドーセルグリッド１０８に接触するよう定義される。そのため、３次元で見た場合、ユニットセル１１３は立方体であることが好ましい。ドーセルグリッド１０８内の複数のドーセルは重なり合うため、連続するエネルギー堆積がオーバーサンプルされ、収束時間は短縮される。
【００２１】
前記のドーセル容積が定義されると、ＣＴ法かＭＲＩ法によって算出した各ドーセルが重なり合うボクセルの質量密度データをもとに、各ドーセルの質量が計算される。閉形式（A  closed  form）の解決法を、以下に示す。粒子トランスポートビーム１１２は、前記の標的質量に向けて放射線を投下する。この粒子トランスポートビーム１１２は、フォトン、ニュートロン、プロトン、α粒子、エレクトロン、ポジトロンなどの数種類の亜原子粒子で構成されていても構わない。前記の粒子トランスポートビーム１１２が標的質量と相互作用し、ビーム１１２の入射粒子が１１４の方向に散乱し、エレクトロン１１６あるいはポジトロンを放出する。本発明では、モンテカルロ法アルゴリズムを使用して、入射分散された粒子、エレクトロン、ポジトロンのうちの、前記のドーセルグリッド１０８で定義された標的質量の対象範囲を最も通過し易い粒子を選択的にトラッキングできるよう分散減少法（variance  reductions  techniques）を採用している。標的質量の対象範囲を通過しにくい、入射分散された粒子、エレクトロン、ポジトロンはトラッキングされない。この選択的トラッキングによって、計算の必要電力と時間の大幅削減が可能になる。ドーセルグリッド１０８の設定、ドーセル質量計算、放射線トランスポート（すなわち粒子トラッキング）、線量計算の詳細を以下に述べる。
【００２２】
図２には、１組の入力パラメータ２００のブロック図を示す。この入力パラメータ２００は、粒子ビーム特性２０２、粒子照射設定２０４、ボクセルグリッド設定２０６、標的質量特性２０８、ドーセルグリッド設定２１０から構成される。この粒子ビーム特性２０２には、前記の粒子ビーム１１２の粒子の種類と、粒子ビーム１１２を定量化するような周知の多様な原子や原子核データが含まれている。この粒子照射設定２０４には、粒子ビームアラインメントと、全ビーム数についてのパラメータが含まれる。このボクセルグリッド設定(set-up)２０６は、前記のボクセルグリッド１０４の配置と間隔を指定する。この標的質量特性２０８では、前記の各ボクセル１０６内の標的質量の質量密度や原子構成などを、ＣＴ法あるいはＭＲＩ法で得たデータから算出する。本発明を説明する目的で前記のドーセル１１０は球体にしてあるが、ドーセル１１０を立方体にすると、オーバラップする部分の容積計算の数学的複雑さが大幅に軽減可能であることが当業者には理解できよう。
【００２３】
前記のドーセルグリッド設定２１０は、前記のドーセル１１０をドーセルグリッド１０８内に幾何学的に配置するための様々な細部を指定する。つまり、各ドーセル１１０の半径、標的質量内のドーセル１１０の中心の始点（ｘ,ｙ,ｚ）、その始点を中心とするドーセル１１０の配列（ix,jy,kz）である。記号i,j,kは、デカルト座標系の各球体の位置を指定する変数である。この配列は、始点（1x,1y,1z）を中心とする単一球体、一列状の複数球体（ix,1y,1z）、平面状の複数球体（ix,jy,1z）あるいは、複数球体の容積（ix,jy,kz）のいずれでもよい。このようにして、前記のドーセルグリッド１０８は、標的質量内のいかなる所望の容積もカバーできる。
【００２４】
図３には、本発明を実施するためのコンピュータ３０２のブロック図を示す。このコンピュータ３０２は、本発明の多様な工程を実行する一連のソフトウェアモジュールを備える。このコンピュータ３０２には、一般的な単一プロセッサをもつパソコンを使用しても構わないが、もっとパワフルな装置で実施することも可能であることは、当業者にとっては明白であろう。コンピュータ３０２は、コンピュータ内の処理ユニットのデータの受信、伝送、出力を制御するコンピュータプログラム命令が格納される内部メモリーを備えており、この内部メモリーは、揮発性メモリーと不揮発性メモリーの両方を備える。この内部メモリーを、コンパクトディスク、フラッシュメモリーカード、磁気ディスクドライブ、ダイナミックランダムアクセスメモリーなどのコンピュータ記録媒体で補足可能であることは、当業者には理解できるであろう。
【００２５】
前記のコンピュータ３０２が前記の粒子トランスポートビーム１１２から投下された放射線量を算出する前に、ユーザは前記のドーセルグリッド１０８の解像度（resolution）と物理的次元を指定する必要がある。高解像度における計算では、統計的不確定性を決めるのに時間がかかるが、ドーセルの空間的分布のより詳細な情報を得ることができる。低解像度における計算では、同等の統計的不確定性（statistical  uncertainty）をより速く算出できる。
【００２６】
前記のコンピュータ３０２に内蔵された一般的な容積計算モジュール３０４は、ライン３０６で入力パラメータ２００を受信し、ドーセル集合体１１０内の各ドーセルと、そのドーセルとオーバラップするボクセルのサブ集合１０６との間の共通容積を算出する。一般的に前記のドーセルグリッド１０８には多数のドーセル１１０が含まれ、ボクセルグリッド１０４も多数のボクセル１０６をもつ傾向があるため、一連の共通容積を計算するには、かなりの計算時間が必要となるであろう。
【００２７】
計算速度をさらに上げるために、本発明では、球体ドーセルと直方体ボクセル間の共通容積を算出する閉形式解法を開示している。球体ドーセルと直方体ボクセル間の共通容積の計算法である「閉形式解法」という表現は、ここで始めて提示するものなので、その定理を展開しておく。
・まず、半径（ｒ）の球体を考慮する。
・前記の球体の中心を原点とする（x,y,z）座標系を定義(define)し、そのx,y,zを
前記の球形の半径で正規化する（つまりx,y,zをｒで割る）
・前記の球体の１/８象限にのみオーバラップするボクセルを考慮する。[注意：このボクセルは複数のサブボクセルにスライスされ、各サブボクセルを球体の1つの１/８象限に配置する]
・前記のボクセルを、前記の球体のx,y,z全てが正値である１/８象限に配置する。[閉形式解法をその他の１/８象限へ同様に拡張する。]
・「コーナ容積（corner  volume）」V3（x,y,z）、つまり、点（x,y,z）において、以下によって限定される容積を定義(define)する。
・（x,y,z）を通るｘｙ面
・（x,y,z）を通るｘｚ面
・（x,y,z）を通るｙｚ面
・前記の球体の表面
この定義に合う８つの容積のうち、（x,y,z）が位置する球体の１/８象限内に全部が収まる容積を選択する。
・原点に最も近接するボクセルのコーナを（ｘ０、ｙ０、ｚ０）とし、そのボクセルの各面の幅をそれぞれdx,dy,dzとすると、x,y,zはｒで正規化されているので、前記のボクセル幅値はdx*ｒとなる。
・次に、前記のボクセルと前記の球体間の共通容積の閉形式解法式は、代入して下記式（１）のようになる。
【数５】
0003730514.tif 000002
【００２８】
この式の最初の項は、前記の球体とボクセルの共通容積に、その他の容積を加えたものであり、後の項は、それら容積を差し引いたものである。
【００２９】
前記の閉形式解法の手順を以下に示す。
・まず、Ｖ３（ｘ，ｙ，ｚ）は、下記の式（２）で表す。
【数６】
0003730514.tif 000003
ただし、
【数７】
0003730514.tif 000004
かつ、
【数８】
0003730514.tif 000005
【００３０】
上記の式（２）は、以下のように求める。
・まず、
【数９】
0003730514.tif 000006
・  V2（x,y）は以下で定義される容積の１／２である。
・（x,y,０）を通るxz平面
・（x,y,０）を通るyz平面
・前記の球体の表面
これら４つの容積から、（x,y）を含む１/４象限に位置する容積を選択する。
球体の対称性から、V3（x,y,z）とV2（x,y）は、各引数の置換において不変値であることがわかる。V3（x,y,z）をV2（x,y）、V2（y,ｚ）などで表記することも可能であるから、V2（x,y）の式が導き出される。
・点（x,y,z）を通るｘｙ、ｘｚ、ｙｚの平面を描く。これらの平面は、点（x,y,z）を含む球体の１/８象限を、８つのサブ１/８象限に分割する。+ｘ、-ｙ、+ｚ方向に沿って位置する前記のサブ１/８象限の容積をV+-+で表す。V2（x,y）は、前記の１/８象限の一部であり、点（x,y,z）よりも正方向のｘとｙの値をもつ。V2（x,０）は、点（x,y,z）よりも正方向のｘ値をもつ１/８象限の一部である。
・代入から、
【数１０】
0003730514.tif 000007
・上記の式を合わせると、
【数１１】
0003730514.tif 000008
ただし、Voctantは、球体の１/８象限の総容積である。
【００３１】
前記の式（７）の左側部分に、Ｖ+++とＶ---以外のサブ１/８象限が含まれる。しかし、Ｖ+++＝Ｖ３（x，y，z）およびＶ---＝ｒ３ｘｙｚである。
・そこで、式（７）は以下のように書き替えることができる。
【数１２】
0003730514.tif 000009
【００３２】
上記の式（２）の手順を完結するには、V2（x,y）の式を得るための積分が必要である。先ず単位球体上の容積を計算し、その結果をｒ３で乗算する。
・図７と図８を参照してほしい。図７は、３次元球体のグラフである。図８は、任意値ｙにおける図７の３次元球体の２次元スライスを示すグラフである。
・図７と図８から、以下の式が成り立つ。
【数１３】
0003730514.tif 000010
【００３３】
Ａ（ｘ０、ｙ）は2番目の図の、パイ形状区域から三角形状を差し引いたものであることが判る。
・次に、ｙからｔへ積分変数を変更する。
ただし、
【数１４】
0003730514.tif 000011
・「Mathematica」（ウルフラム研究所、イリノイ州、キャンペーンユバナ）のソフトを使って、式１０の不定積分値ｈ（ｘ０，ｔ）を算出する。
・その結果を簡略化すると、以下のようになる。
【数１５】
0003730514.tif 000012
【００３４】
前記のMathematicaには、積分定数が含まれている。式（１１）の導関数から、それが式（１０）の積分値であることが判明する。
・式（１０）からの積分値を限定すると、下記のように簡略化できる。
【数１６】
0003730514.tif 000013
ゆえに、式（１３）が成り立つ。
【数１７】
0003730514.tif 000014
・式（１２）と（１３）を式（８）に代入した結果が、式（２）である。
【００３５】
同等のドーセルグリッド容積をカバーできるように、単位セル１１３（つまり、ドーセルグリッドの間隔）を縮小し、ドーセルの総数を増加することで、共通容積の計算速度をさらに高めることが可能である。さらに多くの数のドーセルでボクセルグリッド１０４を構成する場合、ドーセルとボクセルの組み合わせパターンがドーセルグリッド１０８内で周期的に反復されることになる。この周期性の結果、１周期間におけるドーセルとボクセルの間の共通容積1組を算出するだけでよい。また、この特性から、最初の１組の共通容積を、その後に続く全ての共通容積に再使用することもできる。
【００３６】
ドーセルグリッドの間隔は、次のように縮小されることが好ましい。
ｍDｘｖ  ＝ｎDｘｄ        （１４）
ただし、
・ｎとｍは、それぞれ整数である。
・Dｘｄは、ドーセルグリッド１０８の間隔である。
・Dｘｖは、ｘあるいはｙのボクセルグリッド１０４の間隔であり、
このｘとｙとは、単一CTあるいはMRIスライスのグリッド軸である。
式（１４）で使用されるｘあるいはｙの座標は、ドーセルグリッド１０８の領域内に多数のボクセルグリッド点を有する座標であり、
「長方向」と分類されている。ボクセルグリッド点とは、ボクセル１０６の各コーナである。
・ｍ、ｎ、Dｘｄは、（ｎDｘｄ）が前記の長方向におけるドーセルグリッドの幅のおよそ１／ｑなるよう選択される。共通容積1組を迅速に計算する
ことと、単位セル１１３内の分数変化を小さくすることの適当なバランスを取るには、「ｑ」を５に設定すればよい。
【００３７】
前記の式（１４）が真のとき、ボクセルグリッド１０４とドーセルグリッド１０８とは相同期性（co-synchronous）をもつ。相同期性の状態では、ボクセル１０６とドーセル１１０の共通容積は、ドーセルグリッド１０８において周期ｎ、ボクセルグリッド１０４において周期ｍ、という周期性をもつため、第１組のｎ個のドーセルで計算するだけでよい。さらに、Dｘｖ＝Dｙｖの場合、ボクセルグリッド１０４とドーセルグリッド１０８とが、ｘとｙの両方向において相同期性をもつので、ボクセル１０６とドーセル１１０の共通容積を求めるには、両軸の第１組のｎ個のドーセルを計算するだけでよい。このことによりｑ＝５での計算速度は、２０倍ほど向上する。DｘｖがDｙｖと等しくない場合、ボクセルグリッド１０４とドーセルグリッド１０８の相同期性は、長方向においてのみ存在する。Dｙｄが２０％以上変化すると、ユーザが入力したDｙｄの元の値は保存され、全部のドーセル１１０の共通容積が算出される。
【００３８】
前記の閉形式解法を使って、ドーセルグリッドの間隔を縮小しておけば、共通容積計算モジュール３０４で、各ドーセル内に位置するボクセルグリッド１０４の各グリッド点のＶ３（x,y,z）値を算出できる。ドーセルの１つ以上の１/８象限内に存在するボクセル１０６の部分は、式（２）を利用するためには、個別のボクセル値として扱わねばならない。このことは、ボクセルグリッド１０４上にはない下記のそれぞれの点におけるＶ３値を計算することを意味する。
・ボクセルの辺とｘ＝０、ｙ＝０、ｚ＝０の平面との交点。
（ｘ,ｙ,ｚ）＝（０,０,０）はドーセルの中心である。
・ボクセルの平面とｘ軸、ｙ軸、ｚ軸との交点。
・前記のドーセルの中心点。
【００３９】
これら新規の点がボクセルグリッド１０４内の点と同じ場合は、Ｖ３は両方の組の点にて算出できる。これらのボクセルグリッド以外の点は、共通容積計算に影響しないゼロ容積のボクセルを増加させ、不必要な計算を行わせる結果となる。それゆえ、共通容積計算モジュール３０４はＶ３の値を使い、式（１）でドーセルの共通容積とドーセルとオーバラップする各ボクセルを計算する。
【００４０】
共通容積計算モジュール３０４で算出した共通容積の集合体を使って、コンピュータ３０２内のドーセル質量モジュール３０８は、共通容積集合体内の共通容積を、共通容積にオーバラップする各ボクセル１０６の質量密度で乗算して、増加ドーセル質量を得る。ドーセル質量モジュール３０８は、入力パラメータ３０６から各ボクセル１０６の質量密度を入力する。そして、共通容積の集合体内の各共通容積の増加ドーセル質量が、計算されるのである。次に、ドーセル質量モジュール３０８では、全部の増加ドーセル質量を合計して１ドーセル分のドーセル質量を算出する。ドーセル集合体１１０内の各ドーセル毎に、ドーセル質量を計算する。
【００４１】
エネルギー勾配グリッドの計算モジュール３０９は、後で説明する荷電粒子トラッキングの準備工程において、エネルギー勾配グリッドを算出する。エネルギー勾配グリッドは、１つのドーセル１１０に到達して浸入できるよう、ドーセルグリッド１０８外に位置するような、荷電粒子に要するエネルギー量を示す1組のエネルギー値で構成されている。ドーセルに近接の荷電粒子は、ドーセルから遠く離れた荷電粒子よりも必要なドーセル浸入エネルギーが少なくすむ。その結果、エネルギー勾配グリッド計算モジュール３０９が算出したエネルギー値は、勾配をもつ。ドーセルグリッド１０８内に全体が収まっていない各ボクセルについても、エネルギー値を計算する。
【００４２】
つまり、本発明における荷電粒子は、その運動エネルギーＴで決定されるような残留の連続減速近似（ＣＳＤＡ）の移動範囲がそのエネルギー堆積が周知なほど小さい場合には、トラッキングの必要がないのである。そのため、計算速度を上げることができる。エネルギー勾配グリッド計算モジュール３０９により、各ボクセルに最小トラッキングエネルギー値Ｔｍｉｎが充てられる。その後の、後述の荷電粒子トランスポート計算工程中に、荷電粒子のＴ値を、荷電粒子を含むボクセルのＴｍｉｎ値と比較する。Ｔ値がそれ以下の場合には、荷電粒子をトラッキングせず、そのエネルギーは標的質量内に堆積されたもの、つまり、娘核粒子に付与されたものとみなして、その娘核粒子を離散的にトラッキングするのである。
【００４３】
前記のＴｍｉｎ値を求めるさい、最初に、正確な平行六面体状のボクセルグリッドを決める。ボクセルグリッド内のボクセルも、正確な平行六面体である。ボクセル内の密度や材質特性は、一定である。放射線トランスポートはボクセルグリッド全体で算出するが、線量はドーセルグリッド内にのみ格納される。Ｔｍｉｎ値の算出手順は、荷電粒子がドーセルグリッドの内側あるいは外側にあるかよって異なる。
【００４４】
ドーセルグリッドの内側にある場合、対応するＣＳＤＡ範囲がボクセルの最小サイズの０．３倍になるよう、各ボクセルにＴｍｉｎ値を割り当てる。完全にドーセルグリッド外部にある第1シェルのボクセルにも、同様にＴｍｉｎ値を割り当てる。ただし、Ｔｍｉｎ値は、液体ボクセルよりも気体ボクセルのほうが小さい。そのエネルギー値にて荷電粒子トラッキングを停止した後に発生する線量の空間誤差は、標的質量についての知識の空間的解決策と比べて、小さくなければならない。
【００４５】
ドーセルグリッド外においては、Ｔｍｉｎ値が、可能な経路全部で、荷電粒子がドーセルグリッド外の第１シェル内のボクセルに到達するのに要する最小運動エネルギー値であるため、Ｔ＜Ｔｍｉｎの荷電粒子はドーセルグリッドには達しない。Ｔ値の距離による変動値は、ＣＳＤＡ近似法を使って求められる。
【００４６】
前記のＴｍｉｎ値の解決法とは、最小エネルギー損失の経路が、低密度域への最短幾何学経路となり、ドーセルグリッドまで曲線経路に沿って低密度値をたどることを確認することである。Ｔｍｉｎ値は、ドーセルグリッドへの経路を実際に最初から最後までたどる荷電粒子の分数近似値に影響されることはない。Ｔｍｉｎ値は、運動エネルギー値がＴ０となるよう、周知の経路をたどる荷電粒子のＴ値を求めるような単純な問題を最初に考慮することにより計算できる。図４に示すような経路に沿った後ろ方向の移動では、荷電粒子は以下ように仮定する。
・運動エネルギーＴ１で材質＃１に入り、曲線経路Ｄｘ１の所定距離を移動した後、Ｔ０でそこから出る。
・運動エネルギーＴ２で材質＃２に入り、曲線経路Ｄｘ２の所定距離を移動した後、Ｔ１でそこから出る。
・同様。
まず、以下のように仮定する。
【数１８】
0003730514.tif 000015
次に、以下のように等式を連続して解いて、Ｔ１，Ｔ２，Ｔ３…を求める。
【数１９】
0003730514.tif 000016
【００４７】
これにより、各材質の境界での等式が解決できる。より迅速な計算処理は、下記のように可能となる。まず、ボクセルグリッド内の材質全部の各Ｔ値の最小阻止能をＳｕ（Ｔ）とする。
Ｓu（Ｔ）＝Min（Ｓ1（Ｔ），Ｓ2（Ｔ），Ｓ3（Ｔ），…）  （３）
上記の等式（１）のＲｕ（Ｔ）を、対応ＣＳＤＡ範囲とする。ボクセルグリッド内の各ボクセルの阻止能をＳｕ（Ｔ）で置き換えると、等式（２）は以下のようになる。
【数２０】
0003730514.tif 000017
つまり、
【数２１】
0003730514.tif 000018
Ｔｉを式（２）の解とし、Ｔｉ’を式（５）の対応解とする。Ｔｉ’＜＝Ｔｉである。それゆえ、式（２）に換えて式（５）を使うと、荷電粒子が必要量より低いエネルギー量でトラッキングされても、Ｔｍｉｎの計算速度を早くすることができる。
【００４８】
さて、Ｔｏ＝Ｔｍｉｎの式（５）にしたがって第１シェル内の各ボクセルのＴ値を最小にするような、ドーセルグリッド外の所定点からドーセルグリッド外の第１シェルのボクセルのまでの経路を求める問題を考えてみる。その課題は、Ｒｕ（Ｔ）がＴの単純増加関数であるため、下記の式（６）で示す経路長Ｐを最小にすることと等価である。
【数２２】
0003730514.tif 000019
ここで、ｎはＴｍｉｎが算出されるボクセルである。
【００４９】
問題を扱い安くするため、以下のような近似条件を使う。
１）Ｐは、各ボクセルの中心点にて算出する。
２）ｘ、ｙ、ｚ軸に平行な方向の中心点間の経路だけを考慮する。
これらの近似条件は、以下の項目により補正できる。
１）Ｐから、短経路がボクセルに入る軸に沿ってボクセル長さの１/２を減算する。これにより、ドーセルグリッドに最近接のボクセルの端部で、ボクセルのＴｍｉｎ値が決定する。
２）Ｐｍｉｎに１/ｓｑｒｔ（３）を乗算する。これにより、許容されない対角線経路を補正できる。
【００５０】
前記の短経路は、ドーセルグリッドからドーセルグリッドを囲む連続複数シェルのボクセルまで解答値を外向きに延ばすことにより求められる。前述したように、第１シェルボクセル内においては、式（６）のＴ０値は対応するＣＳＤＡ範囲がボクセルの最小サイズの０．３倍となるよう決められる。第２シェルボクセルでは、Ｐｍｉｎの初期値が、第１シェルボクセルから外側へまっすぐ経路を伸ばして設定される。そしてプログラムにより、第２シェルボクセルの内側を掃引して短経路を求めるのである。
【００５１】
前記のような掃引法が、図５の一次元図を説明するのに最も簡単である。ここで、第２シェルボクセル内の線分に沿ったボクセルインデクスをｊとし、Ｐｊの実線３３０は、第１シェルボクセルから外側へ経路を延ばして、等式（６）から得られる経路長である。
【数２３】
0003730514.tif 000020
ただし、Ｄｘはボクセルの幅である。これにより、図の右側に点線３３２が形成される。その点線の各ボクセルの荷電粒子の経路は短くて、第１シェルボクセルへとまっすぐ下降する代わりに、第２シェルボクセルの左へ移動してから第１シェルボクセルまで下方へ伸びている。図の左側の点線３２４は、下記の式（８）のごとく、ｊを減少させる方向に掃引すれば形成できる。
【数２４】
0003730514.tif 000021
それゆえに、解の短経路は、掃引の後、現在のボクセルの前のボクセルが現在のボクセルに繋がってその経路を短くできるかどうかの判断により、低点３３６から伸ばすことができる。
【００５２】
二次元の低点からの延長は、+ｘ方向および+ｙ方向を見据えつつ、+ｘ方向および+ｙ方向へのラスタ走査により行える。その結果、（+ｘ，+ｙ）１/４象限内での低点を大きくできる。低点を全方向に大きくするためには、特に、（+ｘ，+ｙ）、（+ｘ，−ｙ）、（−ｘ，+ｙ）、（−ｘ，−ｙ）の４つの方向にラスタ走査する必要がある。
【００５３】
ゆえに、短経路を求める一般的な方法は、下記のようにまとめられる。短経路を求める手順は、以下の工程からなる。
１）対応するＣＳＤＡ範囲が第１シェルボクセル内の各ボクセルの最小サイズの０．３倍となるよう、ドーセルグリッドおよびドーセルグリッド外の第１シェルボクセルのＴｍｉｎ値を算出する。
２）いくつかの表面が無くなるような、第１シェルボクセル内の各ボクセルから第２シェルボクセル内の隣接ボクセルまで直線の経路を求める。
３）以下の工程を実行して、下記の低点が第２シェルボクセルのコーナにて形成できるよう第２シェルボクセルをラスタ走査する。
Ａ）Ｚ平面のｚの最小値を選択する。
Ｂ）図６の経路３３８に沿って第２シェルボクセルを走査する。
Ｃ）経路３３８を２回まわる。
Ｄ）経路３３８を走査している間、第１のボクセルがｘ方向かｙ方向のい
ずれに存在するか、また、ｚ方向に存在する第２のボクセルが現在のボ
クセルへ繋ぐことにより経路長を減少させるかどうかを調べる。
Ｅ）その次に高い値にｚを設定する。
Ｆ）全部のｚ値に対して上記の工程ＡからＥまでを繰り返す。
４）前記の工程３は、右回りネジのネジ軸下方へネジ溝に沿った移動方向に存在するボクセルを調べつつ、右回りネジのネジ溝に沿って正のｚ方向に移動するのと同じことである。次に、移動方向を調べつつ、以下のような走査を行う。
Ａ）右回りネジに沿って負のｚ方向に移動する。
Ｂ）左回りネジに沿って正負のｚ方向に移動する。
Ｃ）ｘ方向やｙ方向についても同じようなネジ走査を行う。これで、第２
シェルボクセルの走査が完了する。
５）前記の工程２から４までを繰り返して、各追加シェルのボクセルの経路長を求める。
６）最も外側のシェルのボクセルを処理した後、工程３と４を繰り返しつつ、シェルを越えてドーセル領域へと戻る。この内方向の走査では、３つ全部の移動方向（現在シェルの２方向、現在シェルに垂直な内方向）に対する検査を行う。これで、エネルギー勾配グリッド計算が完了する。
【００５４】
次に、コンピュータ３０２内の放射線トランスポートモジュール３１０では、共通容積計算モジュール３０４からのドーセルとボクセルの位置、間隔、共通容積などの情報が入力される。また、放射線トランスポートモジュール３１０は、エネルギー勾配グリッド計算モジュール３０９からエネルギー勾配グリッドを受け取る。さらに、放射線トランスポートモジュール３１０には、ライン３１２の粒子トランスポートビームデータも入力される。粒子トランスポートビームデータには、位置、エネルギー量、統計学的重量、粒子の種類、方向、照射源で生成された粒子が標的に到達して相互作用する分布確率などの、粒子トランスポートビーム１１２内の粒子の情報が含まれている。本発明におけるドーセルグリッド１０８はボクセルグリッド１０４と分離されているため、ドーセルグリッド１０８内の、あるいは、それに浸入しそうな粒子のみが放射線トランスポートモジュール３１０でトラッキングされる。それゆえ、計算能力は、ドーセルグリッド１０８外の、あるいは、それに侵入しそうでない粒子のトラッキングには向けられない。
【００５５】
なお、フォトンやニュートロンなどの非荷電粒子、および、電子、ポジトロン、プロトン、アルファ粒子などの荷電粒子は、異なる手順により放射線トランスポートモジュール３１０でトラッキングされる。
【００５６】
非荷電粒子は、放射線トランスポートモジュール３１０の相互作用毎処理によりトラッキングされ、明確な座標（ｘ，ｙ，ｚ）は、各相互作用に関連するエネルギー堆積量に対して特定される。全部の一次非荷電粒子（つまり、粒子トランスポートビーム１１２の一部）のトラッキングは、標的質量と相互作用するまで、あるいは、ボクセルグリッド１０４から離脱する(leave)まで続けられる。しかしながら、二次非荷電粒子（つまり、娘核（daughter  particle）とも呼ばれるような、粒子トランスポートビーム１１２と標的質量との衝突によって生じた粒子）は、分散減少法(variance reduction)で処理される。
【００５７】
二次非荷電粒子がドーセルグリッド１０８内で生成された場合には、一次非荷電粒子で説明した方法でトラッキングされる。また、二次非荷電粒子がドーセルグリッド１０８外で生成された場合でも、ドーセルグリッド１０８へ向かう場合には、同様に、一次非荷電粒子の方法でトラッキングされる。
【００５８】
けれども、二次非荷電粒子がドーセルグリッド１０８外で生成されて、ドーセルグリッド１０８から離反する方向へ向かう場合には、ロシアンルーレット法や粒子スプリット法などの周知の分散減少法でそれらの二次非荷電粒子を処理する。ロシアンルーレット法や粒子スプリット法については、Ｌ．カーターとＥ．Ｄ．キャシュウェルらの「モンテカルロ法による粒子トランスポートシミュレーション」、ＥＲＤＡクリチカルレビュシリーズ、１９７５年刊、１６ページ、および、アイバン．ラックスとラズロー．コブリンガーらの「モンテカルロ粒子トランスポート法：ニュートロンとフォトン計算」、ＣＲＣ出版、ボカラトン、１９９１年刊、５８ページに記載されている。
【００５９】
もし、ロシアンルーレット法を行った後でもまだ二次非荷電粒子が存在している場合は、それらを一次非荷電粒子の場合と同様の方法でトラッキングする。しかしながら、粒子スプリット法は、二次非荷電粒子がドーセルグリッド１０８に再浸入したときに、統計的変動（statistical  fluctuations）を減らすため使われる。
【００６０】
前記の放射線トランスポートモジュール３１０により、ドーセルグリッド１０８を通過する荷電粒子全部がトラッキングされる。ドーセルグリッド１０８外の荷電粒子については、放射線トランスポートモジュール３１０で分散減少法の範囲阻止を行う。範囲阻止（range  rejection）は、荷電粒子をドーセルグリッド１０８内にトランスポートできるエネルギー（範囲）量をもたない粒子をトラッキングしない（阻止する）方法である。範囲阻止法については、ジェンキンスらによる「電子とフォトンのモンテカルロトランスポート」、プレナム出版、ニューヨーク、1988年刊に詳しく記述されている。
【００６１】
範囲阻止を始めるにあたって、放射線トランスポートモジュール３１０は、ライン３１２に入力された粒子ビームトランスポートデータを使って、荷電粒子に関連するエネルギー値（Ｅ０）を算出する。放射線トランスポートモジュール３１０は、エネルギー勾配グリッド計算モジュール３０９で予め作成されたエネルギー勾配グリッドにアクセスして、荷電粒子がドーセルグリッド１０８に浸入するために要するエネルギー量（Ｅ１）を求める。Ｅ０がＥ１より小さい場合には、荷電粒子はドーセルグリッド１０８に浸入するために要するエネルギー量をもたないので、放射線トランスポートモジュール３１０は荷電粒子のトラッキングを停止する。Ｅ０がＥ１より大きいか同じである場合には、放射線トランスポートモジュール３１０による荷電粒子のトラッキングが継続される。範囲阻止法によるエネルギー勾配グリッドの利用を、ここでは、勾配荷電粒子トラッキングと呼ぶ。
【００６２】
トラッキングされる荷電粒子について、放射線トランスポートモジュール３１０では、荷電粒子を、連続的かつ離散的な相互作用からエネルギーを堆積するかのようにモデル化する。そのモデルを使って、放射線トランスポートモジュール３１０は、ボクセルグリッド１０４に達するまでの荷電粒子の経路をトラッキングする。荷電粒子は、経路中に所定のエネルギー量を堆積させる。
次に、放射線トランスポートモジュール３１０は、経路に沿った数カ所のランダム地点（ｘ，ｙ，ｚ）にエネルギー値を割り当てる。荷電粒子の経路のランダム地点に割り当てられたエネルギー量の和は、元の所定エネルギー量と同じである。ランダム地点のいくつかは、ドーセルグリッド１０８内に位置し、それ以外はその外に位置しても構わない。反対に、荷電粒子が所定位置（ｘ，ｙ，ｚ）における離散的な相互作用でエネルギーを堆積する場合、放射線トランスポートモジュール３１０は所定位置に堆積エネルギー量を割り当てる。
【００６３】
前記の放射線トランスポートモジュール３１０からの出力は、さまざまな地点（ｘ，ｙ，ｚ）における多量のエネルギー堆積量から成る。このことは、ボクセルあたりのエネルギー堆積量を生成するたけの標準的なモンテカルロ放射線トランスポートコード法とは、まったく逆である。
【００６４】
前記のコンピュータ３０２内の線量計算モジュール３１４では、ドーセル質量モジュール３０８からの各ドーセル１１０のドーセル質量と、放射線トランスポートモジュール３１０からの位置（ｘ，ｙ，ｚ）のエネルギー堆積量とが入力される。始めに、線量計算モジュール３１４はエネルギー堆積の位置（ｘ，ｙ，ｚ）を特定する。続いて、線量計算モジュール３１４は、そのオーバラップするドーセル１１０が位置情報（ｘ，ｙ，ｚ）を含むかを特定する。位置（ｘ，ｙ，ｚ）に関連するエネルギー堆積量を、（ｘ，ｙ，ｚ）を含む各ドーセルの質量で割って、各ドーセルの増加線量を求める。線量計算モジュール３１４により、各エネルギー堆積工程が実行されて、各ドーセルの総線量が決まる。線量計算モジュール３１４が、各ドーセル１１０の１組の線量値をライン３１６に出力する。
【００６５】
図９は、モンテカルロ法を基にした粒子トランスポートグリッドのサブ容積内の放射線量の計算方法のフローチャートである。計算方法の最初の工程６０２では、共通容積計算モジュール３０４とエネルギー勾配グリッド計算モジュール３０９に、粒子ビーム特性、粒子照射設定、ボクセルグリッド、多様な標的質量特性、ドーセルグリッドなどが入力される。工程６０４では、エネルギー勾配グリッド計算モジュール３０９により、荷電粒子トラッキングのためのエネルギー勾配グリッドが作成される。工程６０４は、図１０に拡大図示されている。次の工程６０６では、共通容積計算モジュール３０４により、ドーセルグリッド１０８内の各球形ドーセルとボクセルグリッド１０４内の各平行六面体ボクセル間の共通容積の集合体が算出される。ドーセル質量モジュール３０８では、各共通容積が、対応するボクセル質量密度により乗算されて、増加ドーセル質量値が求められ、続く工程６０８で、対応するドーセルの総ドーセル質量に加算されるのである。さらに、工程６１０で、放射線トランスポートモジュール３１０により、標的質量を通過する非荷電粒子の選択的トラッキングが行われる。工程６１０は、図１１に拡大図示されている。そして工程６１２において、放射線トランスポートモジュール３１０は、標的質量を通過する荷電粒子を選択的にトラッキングする。工程６１２も、図１２に拡大図示されている。それから、放射線トランスポートモジュール３１０により、工程６１４で、前述のように、標的質量と離散的に相互作用する非荷電粒子や荷電粒子のエネルギー堆積量と位置とが特定される。工程６１６では、放射線トランスポートモジュール３１０により、同じく前述のように、各連続ランダム位置間の間隔に比例するエネルギー量が堆積される、標的質量と連続的に相互作用する荷電粒子のランダム位置が作成される。次の工程６１８では、線量計算モジュール３１４が、各ドーセルに対して、ドーセル内に堆積されたエネルギー量を加算し、ドーセル質量で割り、ドーセルに付与される放射線量を算出するのである。工程６１８の後、放射線量計算の手順は終了する。
【００６６】
図１０は、エネルギー勾配グリッドの計算方法（図９の工程６０４）のフローチャートである。最初の工程７０２では、エネルギー勾配グリット計算モジュール３０９により、ドーセルグリッド１０８内に完全に収まらないボクセルの集合体が特定される。次の工程７０４では、エネルギー勾配グリット計算モジュール３０９は、荷電粒子をそのボクセル集合体内のボクセルからドーセルグリッドまで移動させるにはどれくらいのエネルギー量が必要かを算出する。工程７０６では、エネルギー勾配グリット計算モジュール３０９が、その算出エネルギー値をボクセルに割り当てる。続いて、工程７０８で、エネルギー勾配グリット計算モジュール３０９により、ボクセル集合体内の残りのボクセルの対するエネルギー値の計算と割り当ての処理が繰り返される。そして工程７０８の後、エネルギー勾配グリット計算の手順が終了する。
【００６７】
図１１は、非荷電粒子のトラッキング方法（図９の工程６１０）のフローチャートである。最初の工程８０２では、放射線トランスポートモジュール３１０により、入力（一次）非荷電粒子が標的質量と相互作用するか、あるいは、それから離脱する(leave)までそれら非荷電粒子はトラッキングされる。次の工程８０４では、放射線トランスポートモジュール３１０が、ドーセルグリッド１０８内で形成された娘核（二次）非荷電粒子を、標的質量と相互作用するか、あるいは、それから離脱するまでトラッキングする。工程８０６では、放射線トランスポートモジュール３１０により、ドーセルグリッド１０８外で形成されその方向に向かう娘核非荷電粒子が、標的質量と相互作用するか、あるいは、それから離脱するまでトラッキングされる。続く工程８０８では、放射線トランスポートモジュール３１０により、ロシアンルーレット法や粒子スプリット法などの分散減少法を利用する場合、ドーセルグリッド１０８外で形成され離間する方向に向かう娘核非荷電粒子がトラッキングされる。そして工程８０８の後、非荷電粒子トラッキング工程は終了する。
【００６８】
図１２は、荷電粒子のトラッキング方法（図９の工程６１２）のフローチャートである。最初の工程９０２では、放射線トランスポートモジュール３１０により、ドーセルグリッド１０８を通過する全部の荷電粒子がトラッキングされる。次の工程９０４では、放射線トランスポートモジュール３１０が、ドーセルグリッド外の荷電粒子について、荷電粒子の位置を特定する。工程９０６では、放射線トランスポートモジュール３１０により、荷電粒子に関連するエネルギー量（Ｅ０）が算出される。次の工程９０８では、放射線トランスポートモジュール３１０は、エネルギー勾配グリッドにアクセスして、荷電粒子が前記の位置から移動してドーセルグリッド１０８に浸入するために要するエネルギー量（Ｅ１）を求める。工程９１０では、放射線トランスポートモジュール３１０により、Ｅ０がＥ１より小さい場合には、荷電粒子のトラッキングが停止される。続く工程９１２で、放射線トランスポートモジュール３１０は、Ｅ０がＥ１よりも大きいか同じである場合には、荷電粒子のトラッキングを継続する。そして工程９１２の後、荷電粒子トラッキング工程を終了する。
【００６９】
以上、本発明の好適実施例を説明したが、それに関するさまざまな変更が可能なことが当業者には理解できよう。本発明による好適実施例の変更や修正は、付随する発明の請求範囲内に制限されるものである。
【図面の簡単な説明】
【図１】  本発明による、モンテカルロ法に基づいた粒子トランスポートグリッドのサブ容積内における放射線量計算システムを示した図である。
【図２】１組の入力パラメータを示すブロック図である。
【図３】本発明を実施するコンピュータシステムのブロック図である。
【図４】荷電粒子の公知の経路を示す図である。
【図５】ある次元掃引方法を示す図である。
【図６】ボクセルの第２シェルを走査する経路を示す図である。
【図７】３次元球体を示すグラフである。
【図８】任意値ｙにおける図４の３次元球体の２次元スライスを示すグラフである。
【図９】モンテカルロ法に基づいた粒子トランスポートグリッドのサブ容積内における、放射線量計算の手順を示すフローチャートである。
【図１０】エネルギー勾配グリッド計算の手順を示すフローチャートである。
【図１１】非電荷粒子トラッキングの手順を示すフローチャートである。
【図１２】電荷粒子トラッキングの手順を示すフローチャートである。
【図１】
0003730514.tif 000022
【図２】
0003730514.tif 000023
【図３】
0003730514.tif 000024
【図４】
0003730514.tif 000025
【図５】
0003730514.tif 000026
【図６】
0003730514.tif 000027
【図７】
0003730514.tif 000028
【図８】
0003730514.tif 000029
【図９】
0003730514.tif 000030
【図１０】
0003730514.tif 000031
【図１１】
0003730514.tif 000032
【図１２】
0003730514.tif 000033