(57)【要約】粒子および波動のシミュレーションを行うための本システムおよび方法は、核スペクトルおよび全てのスペクトル放射輸送、量子粒子輸送、プラズマ輸送および帯電粒子輸送を伴う計算に関して有用である。本発明は、一般的な３次元問題に埋め込むための正確な変数を生成するメカニズムを提供し、一連の単純単一衝突相互作用有限要素を拡張して、複合多重衝突有限要素を構築することができる手段について説明する。【選択図】図１７
2007531566.tif 000002
【特許請求の範囲】
【請求項１】
媒体を透過する粒子の輸送をアルゴリズムによってシミュレートする方法であって、
（ａ）１組の粒子および環境の初期条件を設定する工程と、
（ｂ）実際の対象物または系からボクセルの計算グリッドシステムを生成する工程と、
（ｃ）計算アルゴリズムを用いて粒子分布の複数の放射線集合を設定する工程と、
（ｄ）放射線集合および適切な統合カーネルを用いて輸送乗数を決定する工程と、
（ｅ）ボクセル相互作用タリー内および／またはボクセルタリーサーフィス上における複数の離散粒子分布を適用することによってシミュレートされた粒子輸送を開始する工程と、
（ｆ）前記輸送乗数を適用して、離散粒子タリーを複数の第１ボクセルから複数の第２ボクセルへ輸送する工程と、
（ｇ）前記粒子分布の放射線集合が前記ボクセルのグリッドシステムを通って順次輸送される間、所定の限界に到達するまで、ボクセル内で粒子タリーを継続する工程、
（ｈ）コンピュータの記録場所内から粒子相互作用タリーをコンパイルし、相互作用モデルを適用することによって、ある期間または世代における散乱、状態および累積相互作用を決定する工程と、
（ｉ）相互作用応答率および／または世代固有値が実質的に収束するまで工程（ｆ）〜工程（ｈ）を繰り返す工程と、
（ｊ）シミュレートされた粒子輸送を示す出力を計算によって生成する工程とを備えていることを特徴とする方法。
【請求項２】
媒質を通過する粒子の輸送をアルゴリズムによってシミュレートする方法であって、
（ａ）１組の粒子および環境の初期条件を設定する工程と、
（ｂ）実際の対象物または系からボクセルの計算グリッドシステムを生成する工程と、
（ｃ）計算アルゴリズムを用いて粒子分布の複数の放射線集合を設定する工程と、
（ｄ）放射線集合および適切な統合カーネルを用いて輸送乗数を決定する工程と、
（ｅ）ボクセル相互作用タリー内および／またはボクセルタリーサーフィス上における複数の離散粒子分布を適用することによってシミュレートされた粒子輸送を開始する工程と、
（ｆ）前記輸送乗数を適用して、離散粒子タリーを複数の第１ボクセルから複数の第２ボクセルへ輸送する工程と、
（ｇ）コンピュータの記録場所内から粒子相互作用タリーをコンパイルし、相互作用モデルを適用することによって、ある期間または世代における散乱、状態および累積相互作用を決定する工程と、
（ｈ）相互作用応答率および／または世代固有値が実質的に収束するまで工程（ｆ）〜工程（ｈ）を繰り返す工程と、
（ｉ）シミュレートされた粒子輸送を示す出力を計算によって生成する工程とを備えていることを特徴とする方法。
【請求項３】
請求項１又は２の方法において、
工程（ｇ）の前に、１組の乗数に関連付けられた第１のボクセルタリー位置がゼロにされることを特徴とする。
【請求項４】
請求項１又は２の方法において、
前記粒子分布の複数の放射線集合を設定するためのアルゴリズム計算は、モンテカルロ技術を用いて行われることを特徴とする。
【請求項５】
請求項１又は２の方法は、
核放射線輸送をモデリングするために用いる複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項６】
請求項１又は２の方法は、
電磁粒子輸送をモデリングするために用いる複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項７】
請求項６の方法において、
前記電磁粒子輸送は、赤外線波、光波、ＵＶ波および電波のうちの１つまたはそれらの組み合わせを含むことを特徴とする方法。
【請求項８】
請求項１又は２の方法は、
放射熱移動をモデリングするために用いる複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項９】
粒子輸送を計算するために、離散粒子タリーを示すコンピュータメモリロケーションをリンクさせる方法であって、
（ａ）ボクセルボリューム、サーフィスおよび関数係数に離散粒子タリーを割り当てる工程と、
（ｂ）参照離散粒子タリーから輸送された粒子の一部を表す略正確な乗数を、複数の適切な離散位相状態変数を有する隣接する離散粒子タリー位置に割り当てる工程と、
（ｃ）参照タリー位置からの輸送乗数に関連する隣接する離散粒子タリー位置に対してポインタの規則的な指定を行う工程と、
（ｄ）参照粒子タリーの系およびそれに対応する乗数をスイープして、離散粒子を隣接するボクセルボリューム、サーフィスおよび関数係数に輸送する工程とを備えていることを特徴とする方法。
【請求項１０】
請求項９の方法において、
前記複数の離散位相空間／状態変数は、核放射線輸送をモデリングするために使用されることを特徴とする。
【請求項１１】
請求項９の方法において、
前記複数の離散位相空間／状態変数は、電磁粒子輸送をモデリングするために使用されることを特徴とする。
【請求項１２】
請求項１１の方法において、
前記電磁粒子輸送は、赤外線波、光波、ＵＶ波および電波のうちの１つまたはそれらの組み合わせを含むことを特徴とする。
【請求項１３】
請求項９の方法において、
前記複数の離散位相空間／状態変数は、放射熱移動をモデリングするために使用されることを特徴とする。
【請求項１４】
請求項９の方法において、
前記複数の離散位相空間／状態変数は、音波をモデリングするために使用されることを特徴とする。
【請求項１５】
参照ボクセルボリュームまたはサーフィスに関係する隣接ボクセルのローカルシステムを横断する１つ以上の固有の経路を含む放射線集合を用いて幾何学的粒子経路を特定する方法であって、前記放射線集合は、
（ａ）ボクセルのシステムを通る放射線集合探査経路の固有のディスクリプタと、
（ｂ）局所ボクセル群から隣接する局所ボクセル群の放射線集合へ放出される放射線集合をリンクさせる仕様と、
（ｃ）参照ボクセルボリュームまたはサーフィスに関係する隣接ボクセルの不変ローカルシステムを横断する、１つ以上の可変長の代表放射線と、
（ｄ）指定された新生の粒子分布に関係する特定の放射線を横断する粒子の一部と、
（ｅ）横断したボクセル内の放射線に関連付けられた長さとを含むことを特徴とする方法。
【請求項１６】
参照離散粒子タリーから適切な離散位相状態の隣接する離散粒子タリー位置へ輸送された粒子の一部を表す正確または略正確な乗数を計算する方法であって、
（ａ）ボクセル内の物質組成を特定する工程と、
（ｂ）一般システム内の局所ボクセル群を特定する工程と、
（ｃ）粒子と物質との相互作用によるボクセル内の粒子減衰を表す適切な統合カーネルを使用する工程と、
（ｄ）参照ボクセルボリュームまたはサーフィスに関係する隣接ボクセルのローカルシステムを横断する１つ以上の固有の経路を含む放射線集合を生成する工程と、
（ｅ）放射線集合データと共に統合カーネルを適用して、参照離散粒子の、隣接する離散粒子ボクセルボリューム、局所ボクセル群サーフィスおよび関数係数への集合輸送を表す固有の乗数を得る工程とを備えていることを特徴とする方法。
【請求項１７】
請求項１６の方法は、
核放射輸送をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項１８】
請求項１６の方法は、
電磁粒子輸送をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項１９】
請求項１８の方法において、
前記電磁粒子輸送は、赤外線波、光波、ＵＶ波および電波のうちの１つまたはそれらの組み合わせであることを特徴とする。
【請求項２０】
請求項１６の方法は、
放射熱移動をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項２１】
請求項１６の方法は、
音波をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項２２】
粒子輸送の計算方法であって、
（ａ）物理系を表すボクセルのグリッドシステムを構築する工程と、
（ｂ）参照ボクセルボリュームまたはサーフィスに関係する隣接ボクセルのローカルシステムを横断する１つ以上の固有の経路を含む放射線集合を用いて幾何学的粒子経路を特定する技術を用いて輸送乗数を構築する工程と、
（ｃ）粒子輸送を計算するための初期条件を設定する工程と、
（ｄ）離散粒子タリーを表すコンピュータメモリロケーションをリンクさせて、コンピュータメモリ適用輸送乗数の１スイープにおいて粒子を輸送する技術を用いて、インパルシブスイープによる関数係数、サーフィスおよびボリューム離散粒子タリー位置への粒子輸送を決定する工程と、
（ｅ）相互作用モデルおよび得られたボクセルボリュームへの蓄積粒子衝突タリーを用いて粒子の相互作用を計算する工程と、
（ｆ）収束基準に基づいて粒子輸送を終了する工程と、
（ｇ）結果を記憶する工程とを備えていることを特徴とする方法。
【請求項２３】
請求項２２の方法は、
核放射輸送をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項２４】
請求項２２の方法は、
電磁粒子輸送をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項２５】
請求項２４の方法は、
前記電磁粒子輸送は、赤外線波、光波、ＵＶ波および電波のうちの１つまたはそれらの組み合わせであることを特徴とする。
【請求項２６】
請求項２２の方法は、
放射熱移動をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項２７】
請求項２２の方法は、
音波をモデリングするために使用される複数の離散位相空間変数をさらに含むことを特徴とする。
【請求項２８】
請求項２３の方法は、
医療分野での放射線治療計画策定に使用することを特徴とする。
【請求項２９】
請求項２３の方法は、
強度変調放射線療法の３次元的治療計画策定に使用することを特徴とする。
【請求項３０】
請求項１５の方法は、
モンテカルロ法を用いて放射線集合を予め計算することをさらに含むことを特徴とする。
【請求項３１】
請求項３０の方法は、
標準的な幾何学的グリッドシステムに使用することを特徴とする。
【請求項３２】
請求項１５の方法は、
ボクセルシステムを通過する放射線経路を記述するハッシュエンコードされたデータを用いて放射線集合経路を特定することをさらに含むことを特徴とする。
【請求項３３】
請求項１５の方法は、
ニ分木法を用いて放射線集合経路を特定することをさらに含むことを特徴とする。
【請求項３４】
請求項１６の方法は、
モンテカルロ法を用いて放射線集合を予め計算すること、および
割当用ハッシュテーブル内の離散粒子メモリを離散粒子乗数に関連付けるメモリ参照ハッシュポインタを使用することをさらに含むことを特徴とする。
【請求項３５】
請求項１６の方法は、
（ａ）参照ボクセルについてサブサーフィスおよびボリュームを関連付ける工程と、
（ｂ）前記サブサーフィスによって表される離散粒子の一部について、初期粒子重量を決定する工程と、
（ｃ）延長された放射線エンドポイントを用いて隣接するＬＶＧ表面輸送を決定する工程と、
（ｄ）スタートポイントおよびエンドポイントに基づいて放射線集合を決定し、前記放射線集合を参照離散粒子に関連付けて、ＬＶＧ表面から放出される粒子についての放射線集合位相タリービンを生成するシステムを使用する工程とをさらに備えることを特徴とする。
【請求項３６】
請求項３５の方法は、
不規則な幾何学的グリッドシステムに使用することを特徴とする。
【請求項３７】
請求項３５の方法は、
ボクセルボリューム粒子源および相互作用タリーから離散粒子を発出するボリュームについて角度放射線重量を計算することを含むことを特徴とする。
【請求項３８】
請求項２２の方法は、
最適化モジュールを組み込んで、最適な粒子分布を決定する工程をさらに備え、この工程は、
（ａ）輸送計算結果から設計最適化特性を決定する工程と、
（ｂ）試験的最適化のための初期粒子分布を決定する工程とを含むことを特徴とする。
【請求項３９】
ボクセル内での粒子相互作用を計算する方法であって、
（ａ）ボクセルボリューム内での衝突確率を計算する工程と、
（ｂ）相互作用に関連した物理パラメータを計算する工程と、
（ｃ）関数係数、ボリュームまたはサーフィス離散粒子タリー分布を相互作用から計算する工程とを備えていることを特徴とする方法。
【請求項４０】
請求項２２の方法は、
距離モーメントを用いてサーフィス放射粒子分布およびボリューム放射粒子分布を微分する工程をさらに備えていることを特徴とする。
【請求項４１】
請求項２２の方法は、
相互作用モデルを計算するために用いられることを特徴とする。
【請求項４２】
請求項４１の方法は、
（ａ）より粒度の高い相互作用モデルを用いてサブボリューム相互作用を表す工程と、
（ｂ）サーフィス初期粒子分布を用いて相互作用モデル応答を決定する工程と、
（ｃ）ボイドを用いて、サーフィス放射粒子分布およびボリューム放射粒子分布を微分するために使用される適用可能な距離モーメントに依存した適切な放射線集合の割当てを実現する工程とをさらに備えていることを特徴とする。
【請求項４３】
請求項４２の方法は、
前記より粒度の高い相互作用モデルはボクセル内の粒子相互作用を計算することによって構成される方法であって、
（ａ）ボクセルボリューム内での衝突確率を計算する工程と、
（ｂ）相互作用に関連した物理パラメータを計算する工程と、
（ｃ）ボリュームまたはサーフィス離散粒子タリー分布を相互作用から計算する工程と、
（ｄ）高次の粒子散乱モデリングのための角度分布を表す関数係数を計算する工程とを備えていることを特徴とする。
【請求項４４】
請求項１６の方法は、
無限放射線相互作用の直接解析的解および幾何学系全体における輸送積分ならびに複数の状態変数をさらに含み、
（ａ）数値解または解析解を用い、離散位相空間における一定の離散粒子分布を仮定する工程、および／または
（ｂ）角度離散位相空間間隔における適切な角度分布を説明する工程をさらに備えていることを特徴とする。
【請求項４５】
請求項４４の方法は、
前記状態変数のうちの１つが、第１衝突相互作用モーメントについての立体角であることを特徴とする。
【請求項４６】
請求項９の方法は、
不変の離散位相空間値を用いて、隣接する離散粒子タリー乗数を、参照タリーから前記系内の全ての離散粒子へと拡張することをさらに含むことを特徴とする。
【請求項４７】
請求項１５の方法は、
不変の離散位相空間値を用いて、隣接する離散粒子タリー乗数を、参照タリーから前記系内の全ての離散粒子へと拡張することをさらに含むことを特徴とする。
【請求項４８】
請求項９の方法は、
アナログコンピュータシステムまたはアナログデジタルハイブリッドコンピュータシステムにおいて、粒子タリーを直接布線することをさらに含むことを特徴とする。
【請求項４９】
請求項１６の方法は、
パターンマッチングアルゴリズムを用いて輸送乗数の計算処理を高速化することをさらに含むことを特徴とする。
【請求項５０】
請求項３５の方法は、
パターンマッチングアルゴリズムを用いて輸送乗数の計算処理を高速化することをさらに含むことを特徴とする。
【請求項５１】
請求項４１の方法は、
系に入射する代表放射線外部ビームの明示的表現モデリングをさらに含み、
（ａ）前記系内へ流入する粒子の代表放射線モデリング、
（ｂ）前記相互作用モデルを用いて初期の散乱放射線源を生成すること、および
（ｃ）関数係数を用いて初期の散乱放射源分布を生成すること
をさらに備えていることを特徴とする。
【請求項５２】
請求項１５の方法は、
モンテカルロ法を用いて純粋な幾何学的特性を計算することをさらに含むことを特徴とする。
【請求項５３】
請求項１６の方法は、
ソフトウェアアプリケーションにおいて、代替の統合カーネルを用いて粒子輸送をモデリングすることをさらに含むことを特徴とする。
【請求項５４】
請求項１５の方法は、
統合カーネルが明示的に適用された放射線集合幾何学特性について、モンテカルロ事前計算法を用いて上限および下限を設定することをさらに含むことを特徴とする。
【請求項５５】
請求項１６の方法は、
前記統合カーネルが明示的に適用された放射線集合幾何学特性について、モンテカルロ事前計算法を用いて上限および下限を設定することをさらに含むことを特徴とする。
【請求項５６】
請求項３９の方法は、
衝突確率法を用いて相互作用モデルの中核が形成され、前記モデルは衝撃初期値の性質を許容するように改変可能であることを特徴とする。
【請求項５７】
請求項９の方法は、
特定の期間内における粒子輸送解を用いて、過渡システムをモデリングすることを特徴とする。
【請求項５８】
請求項２２の方法は、
特定の期間内における粒子輸送解を用いて、過渡システムをモデリングすることを特徴とする。
【請求項５９】
請求項２２の方法は、
初期衝撃に関連した絶対的システム境界収束を用いることをさらに含むことを特徴とする。
【請求項６０】
請求項５９の方法は、
非核分裂性の非時間的固有値問題に用いることを特徴とする。
【請求項６１】
請求項２３の方法は、
前記相互作用モデル内の核分裂パラメータおよびシステム生成固有値を用いることをさらに含むことを特徴とする。
【請求項６２】
請求項２３の方法は、
前記相互作用モデル内の核分裂パラメータおよびシステム生成固有値を用いることと、前記相互作用モデル内の生成モーメントを用いることとをさらに含むことを特徴とする。
【請求項６３】
請求項２３の方法は、
関数係数デポジション法を用いて、角度的且つ空間的な意味で、サーフィスまたはサブサーフィスタリー分布を空間的にマッピングすることをさらに含むことを特徴とする。
【請求項６４】
請求項２３の方法は、
関数係数デポジション法を用いて、ＬＶＧ境界についてのデータ圧縮ツールとして機能することをさらに含むことを特徴とする。
【請求項６５】
レイトレーシング結果を乗数フィールドに蓄積する方法であって、
（ａ）サーフィス上またはボリューム内の点を選択する工程と、
（ｂ）放射線に関連付けられた重量を決定する工程と、
（ｃ）前記選択された点から第１のボクセルを通過するレイトレーシングを行う工程と、
（ｄ）前記ボクセルのターミナルに関連付けられたポインタを特定する工程と、
（ｅ）前記ボクセル内で減衰した粒子の数を決定する工程と、
（ｆ）前記ポインタを、前記ボクセル内の前記減衰した粒子のターミナルに関連付けられた乗数に関係付ける工程と、
（ｇ）放射線に関係付けられた重量、ボクセル内における減衰、およびターミナルに関係付けられた重量を用いて、特定された乗数に結果を追加する工程と、
（ｈ）前記ボクセル内の前記放射線に適した全てのターミナルポインタが特定されて処理されるまで、工程（ｄ）〜工程（ｇ）を継続する工程と、
（ｉ）前記ボクセルから続く粒子の数を減少させる工程と、
（ｊ）問題エッジまたはターミナルサーフィスに到達するまで、隣接するボクセルへと継続し、ポインタを特定し、乗数を蓄積する工程（工程ｄ〜工程ｈ）と、
（ｋ）前記サーフィス上のターミナルに関係付けられたポインタを特定する工程と、
（ｌ）放射線に関連付けられた重量、減衰された粒子の一部、およびターミナルに関連付けられた重量を用いて、ポインタに関連付けられた乗数に結果を追加する工程と、
（ｍ）前記放射線に適した全てのターミナルが処理されるまで工程（ｋ）〜工程（ｍ）を継続する工程とを繰り返し行うことを特徴とする方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、媒質内を伝播する分子および／または波動の輸送をシミュレートする方法に関する。より具体的には、本発明は、特定の媒質内で伝播し、相互に作用する核放射線および電磁放射線の両方を、計算アルゴリズムおよび方法体系を用いてシミュレートする方法に関する。
【背景技術】
【０００２】
本発明は、多様な用途に利用可能であり、先行技術における数々の問題に対応するものであるが、医療分野での利用が本発明の記載に最も関連しており、先行技術を克服するものである。１世紀以上もの間、様々な種類の癌性腫瘍の治療を目的として、医師たちは患者に電離放射線照射を施してきた。今日までに、腫瘍学分野における放射線療法は大いに進歩している。その技術的進歩のひとつとして、患者への照射前に放射線強度を変調させる技法がある。放射線を変調させると、そのビームは均一でなくなり、線量分布形状が不規則になる。この不規則な線量分布の有利な点は、ビームの最大線量または強度を正確に腫瘍の標的体積の位置に限定して合わせることと、同時に腫瘍周辺の健康な組織に対するビームの線量または強度を最小にすることができることである。
【０００３】
現在、この手法は強度変調放射線治療と呼ばれているが、ここではＩＭＲＴと呼ぶことにする。ビーム強度の変調は慣例的な手法であり、何十年もの間、腫瘍放射線療法の分野において利用されている。最近では、マルチリーフコリメータ（ＭＬＣ）が３次元原体放射線治療（３ＤＣＲＴ）およびＩＭＲＴの両方を実施する機能を有するため、これを用いてビーム強度を変調するのが一般的になってきている。
【０００４】
しかしながら、強度変調ビームが複雑であること、および／またはＭＬＣで生成されたそれらのフィールド形状が変則的であることから、マニュアル計算により患者体内での深部線量分布を測定することは、不可能ではないにしても極めて困難である。このことから、計算アルゴリズムの使用が必要とされている。コンピュータを用いることにより、物質内を伝播し、物質と相互作用する電離放射線を、アルゴリズムによりシミュレートすることができる。したがって、このような計算アルゴリズムにより、面倒なマニュアル計算を行わなくても必要な深部線量分布を得ることができる。
【０００５】
長年にわたり、３次元での粒子輸送の測定にはモンテカルロ法または確率的方法が用いられてきた。このような粒子または波動シミュレーションのための計算方法体系は多くの用途に使用されている。例えば、強度変調放射線治療の３次元的治療計画策定においては、モンテカルロ計算だけでなく、従来の放射線ビルドアップ係数法に基づく半経験的な方法も利用される。また、組織ファントム全体での散乱線量を計算するために、高速フーリエ変換が重ね合わせを用いて反復的に用いられる。原子力工学の分野では、原子炉の炉心設計から核医学防護といった広範な用途のために、α、βおよびγの中性子およびプロトン放射を核輸送計算でシミュレートしている。また、放射線検出器や電子レンジ等の多くの機器の設計にも輸送計算は有益である。
【０００６】
計算方法が利用されるその他の分野として、燃焼器およびボイラ設計において有用な放射熱伝達がある。例えば、「離散的伝達方法」では、境界で接する反射面同士の間の放射線を用いて、検査ボリュームのグリッドシステムを横断する熱放射をモデリングする。この放射線に対応する特定の長さに様々の統合カーネルが適用され、システムを通過するビームの強度が測定される。放射線の強度は、固定された問題境界の検査ボリュームに放射線の反射を伴って出入りする時に計算される。
【０００７】
「不変量埋込み法」では、不変量解の大きなフレーム構造への埋め込みと、有限要素解とノード節点計算とのカップリングが行われる。しかしながら、この方法は、１次元を超える拡張には比較的不正確な、機能性粒子密度の連続性の概念を保持している（とはいえ、不変量を埋め込む有限要素を決定するための離散的輸送方法を使用すると、さらに大きな問題に放射線解を拡張するための実際的な手段が得られる）。
【０００８】
既存の方法は多くの計算用途に適しているものの、それらはいくつもの欠点を有している。１つの重要な欠点は、極めて大きな放射線集合を、一般化されたフィールドで分布ソースからの放射分布の計算に用いることができる簡単なスカラー乗数に換算するためのメカニズムを与えることができないことである。もう１つの欠点は、公知の方法体系では、放射線のローカルシステムがどのように一組の検査ボリュームから別の組の検査ボリュームに拡張されるのかが説明されていないということである。先行技術による方法では、一般化されたフィールドを用いて、大きなシステムに埋め込まれる不変ボクセル群を形成するためのメカニズムが含まれていない。また、既存の方法体系のさらに別の欠点は、一連の単純単一衝突相互作用有限要素（a series of single collision interaction finite elements）を拡張して、複合多重衝突有限要素（complex multi-collision finite element）を構築することができる手段を提供することができないことである。最後に、先行技術による方法は、全ての形状の放射粒子（電磁放射から量子的粒子、中性子から音波）のための一般化スキームが欠如している。したがって、先行技術におけるこれらの欠点を克服できる計算方法体系が要求されている。
【０００９】
上記のような欠点を有することから、先行技術による計算方法は、輸送シミュレーションの解法にはなり得るものの、満足できる速度および精度で粒子輸送をシミュレートすることは不可能である。モンテカルロ法のような典型的な方法では、特定の条件下であれば十分な精度のシミュレーションを行うことができるが、放射線腫瘍学クリニックで使用するのに適した時間内でシミュレーションはできない。また、電子機器設計等のその他の用途においてはあまりにも実験的で、費用も時間もかかってしまう。あるいは、その他のシミュレーション法では必要な速度が得られるが、最も不可欠な精度が犠牲になってしまう。このように、当該分野においては、満足できるシミュレーション速度および精度で放射線輸送シミュレーションを行うことができる計算方法体系が熱望されている。
【発明の開示】
【００１０】
本発明は、粒子および波動のシミュレーションを行うためのシステムおよび計算方法体系に関する。本発明は、核スペクトルおよび全てのスペクトル（電波スペクトルからガンマ線スペクトルまで）の電磁放射輸送だけでなく、量子粒子輸送、プラズマ輸送、および帯電粒子輸送の計算に関して特に有用である。本発明はさらに振動／音波計算および放射熱伝達にも有用である。
【００１１】
本発明は、直接解法で解かれた離散的有限要素における離散粒子の値を優先して、粒子密度の関数表現(functional representations)を減少させるという点において、先行技術からの概念的な脱却を意味する。本発明の有利な点によれば、一般的な３次元（３Ｄ）の問題に埋め込むための正確な不変量（R. Belman, G. m. Wing， An Instruction to Invariant Imbedding, R. Belman, SIAM (1992) ISBN 0-89871-304-8を参照）を生成するメカニズムが提供される。本発明はまた、一連の単純単一衝突相互作用有限要素（a series of simple single collision interaction finite elements）を拡張して複合多重衝突有限要素を構築することができる手段について記載する。
【００１２】
この結果、本明細書に記載する本発明によれば、先行技術に比べて少なくとも１０倍〜１０００倍も高速な計算方法が提供される。また、この結果によれば、直接解法と比較して、放射線効果集約軽減ベンチマーク問題（ray effect intensive mitigation benchmark problems）の顕著な改善が見られる。以下はデジタルコンピュータを用いて得られた結果である。しかしながら本発明では、高速粒子輸送計算という特殊な分野において、アナログ方式またはデジタル−アナログハイブリッド方式を用いることも考慮に入れている。例えば、本発明において、強度変調放射線治療の３次元的治療計画策定に高速アナログ制御システムを用いることができる。これにより、患者の位置データと併用して外部ビーム機器を使用する際に、粒子輸送計算結果がリアルタイムに提供される。
【００１３】
上記したように、ここに記載する実施例は、主に医療用途で放射線分布計算を行うための計算装置に関する。特に、癌の放射線治療に用いられる、人体組織中でのガンマ線およびＸ線粒子輸送を計算するための強度変調放射線治療の３次元的治療計画策定（以下「ＩＭＲＴ  ３ＤＲＴＰ」とする）に関する実施例について言及する。しかしながら、上記実施例は、その他の分野における本発明の応用についても説明する。その分野とは、医療用放射線治療計画策定シミュレーション、電子機器設計（例えば電子レンジおよびコンピュータチップ）、レーダ（電磁放射モデリング）、ソナー（音波シミュレーション）、送電機器設計、光学装置類、放射熱伝達、核遮蔽および原子炉シミュレーション、宇宙物理学、ボイラおよび燃焼器シミュレーション、電気通信（例えば携帯電話用電波塔の最適位置決定）および機械学を含むが、これらに限定されるわけではない。本発明の適合性は、部分的にはその特異なアルゴリズム機構部分に基づくものであり、（Ｃのファンクションポインタまたはオブジェクト指向コンピュータ言語に関して）置換可能な統合カーネルおよび相互作用モデルを用いることによって各種用途のモデリングを行うことが可能である。以下に定義するように、統合カーネルは図１０のＢ５、図１１のＢ３および図１２のＢ３で使用する。相互作用モデルは図７のブロック５に明示的に示される。
【００１４】
−定義−
以下の用語は、本発明の理解を助けるために規定された定義である。これらの用語は、本発明の属する分野において従来から用いられている意味と矛盾するものではない。
【００１５】
粒子
粒子は二重性を有する１つのエネルギーパケットとして定義される。エネルギーパケットは、離散的有限オブジェクトであってもよいし、特定の波長または周波数を有する波動であってもよい。エネルギーの二重性という概念は新しいものではなく、一般的な科学分野では慣例的に受け入れられている。
【００１６】
この定義に当てはまる粒子の例としては、γ線、Ｘ線等の電磁放射線、中性子、アルファ粒子、陽子および電子等のその他の核子が挙げられるが、これらに限定されるわけではない。△ｘ−Ｂｏｌｄ  △ｘは、変数ｘのｘからｘ＋△ｘへの有限間隔を示す。
【００１７】
表面に関して、△ｓは有限領域の離散的表面を示す（弧に対応した長さの単位ではない）。太字の△ｘは本発明に関するものであるが、強調表示していない△ｘは従来の数値増分を意味する。
【００１８】
離散粒子
本発明の重要な構成要素として、離散粒子は、多数の粒子が離散的有限状態変数（すなわち離散的位相空間）内に表現された実数タリー（real number tally）である。
【００１９】
ＩＭＲＴ  ３ＤＲＴＰにおいて、離散粒子はＰ［△ｓ、△Ω、△ｔ］と定義される。Ｐは有限領域△ｓの表面に関連し、立体角集合△Ωの範囲内および時間△ｔの有限間隔の範囲内を横断する粒子の数である。しかし、数値計算装置としては、それは１点にあるのではなく、有限間隔を横切る定数値にある。
【００２０】
上記の離散粒子状態変数△ｓ、△Ωおよび△ｔは、ＩＭＲＴ  ３ＤＲＴＰの放射線輸送をシミュレートするのに好ましい状態変数であるが、本発明はこれらの変数のみに限定されない。その他の各種の状態変数は、離散エネルギー群、電荷、量子スピンまたは放射線集合（以下に定義する）を含むが、これらに限定されない。
【００２１】
ボクセル
ボクセルは、三次元におけるピクセルと考えることができる。ピクセルは２次元平面上での画像に関連するのに対し、ボクセルは対象物の３次元画像に対応する。各ボクセルは、空間面と、その中に制限的に配置された均質な物質組成の両方を表現する。例えば６つの表面を有する６面体のような、ボクセルとしての規則的な形状を視覚化するのはそれほど難しいことではないが、各ボクセルが均一形状である必要はない。例えば、３ＤＲＴＰ  ＩＭＲＴでは、ボクセルは骨の一部として成形されてもよく、臓器の１または複数の部分を含んでいてもよい。このボクセルという用語は、ＩＭＲＴのコンピュータ解法に関する先行技術の記載と共通のものである。
【００２２】
グリッド
本発明において、グリッドは、それぞれ多数の面を有する隣接ボクセルの有限体系として考えられる。例えば、簡単な概念上の例は、立方体を形成する３×３×３のボクセルであり、各ボクセルはそれぞれ６つの面を有する。グリッドシステムは、３次元空間で物質の物理系を表現するために用いられる。ここで強調すべきは、ボクセルおよびグリッドの形状は、規則的であっても変則的であってもよいということである。例えば、３ＤＲＴＰ  ＩＭＲＴでは、ボクセルのグリッドシステムは、これに限定はされないが、例えば臓器および骨等を含むクリティカル領域に適合するように、変則的な形状となっている。
【００２３】
放射線集合
放射線集合は、１からｎまでの複数の代表放射線であり、ｎは有限整数である。これらの放射線は、サーフィスまたはボリュームから放出されてグリッドを横断する全ての放射線の平均相対部分とボクセル経路長を表す。１つの放射線集合のうちのそれぞれの放射線は、固有の部分および固有のボクセル経路長を有する。しかし、各放射線は、同じ放射線集合に含まれるその他の放射線と同じ経路でボクセルグリッドを横断する。
【００２４】
局所ボクセル群
局所ボクセル群（ＬＶＧ）は、グリッドのサブセット、若しくは、特定のボクセルサーフィスまたはボクセルボリュームから参照されたボクセルの集団または集合である。このボクセルサーフィスまたはボクセルボリュームから、放射線集合が局所ボクセルを通って放出され、局所ボクセル群の外部表面に到達する。特定のサーフィスに関して以外は、局所ボクセル群は重複しており、固定された有限要素ではないため、参照されたボクセルサーフィスは重要である。
【００２５】
相互作用モデル
相互作用モデル（ＩＭ）は、ボクセル内で介在物質と衝突したことで一度タリーされた（tallied）粒子の分布計算において適用される方法体系である。相互作用モデルは、粒子の相互作用、反射または散乱ならびに状態変化をモデル化することが要求される。相互作用モジュールは、ボクセルサーフィスから放出される離散粒子の分布を、全ての関連物質および状態変数を考慮することによって決定する。相互作用モジュールはまた、離散的立体角群△Ω内の放射線集合についての粒子分布を決定する。
【００２６】
相互作用タリー（Interaction Tallies）
相互作用タリーの概念は、相互作用モデルの概念と密接に関連している。相互作用タリーは、△ｓ、△Ｅ、△Ω、△ｔおよび△R等の粒子の入射状態の値、若しくは、これらに限定はされないが、量子スピン（量子力学）、力（機械分野）または温度（放射熱移動）といったその他の状態変数を具現化する高度なタリーと比較して、ボクセルボリューム内で生じる相互作用の単一同種計数と同程度に単純である。相互作用タリーは関数係数にも適用できる。例えば、本発明を利用して、角粒子分布の平面または球面の調和的表現に関連付けられた関数係数を直接的に取り込むことができる。
【００２７】
均質化
この用語は、粒子輸送計算に使用される、体積または相互作用の観点から平均された物質特性を指す。
【００２８】
統合カーネル
粒子輸送に関して、統合カーネルとは、粒子に対応した直線的で非衝突の強度カーネルを指す。放射粒子のほとんどの形状に適切な統合カーネルは単純にｅ-Σ△rである。Σは粒子−物質間相互作用の全マクロ断面積、△ｒは特定のボクセル内での横断距離である。これが非衝突の透過粒子を表しているのに対し、１−ｅ-Σ△rは相互作用成分を表している。相互作用モデルが修正Ｐ1等方散乱である場合、Σは全断面に対する輸送断面を表す。
【００２９】
輸送乗数
この用語は、図６および図１０に記載の数値であり、ソース位置（サーフィスまたはボリューム）からターミナル位置（サーフィス、ボリュームまたは関数係数）に輸送された離散粒子のうちの非衝突分の最終合計を示す。
【００３０】
ターミナル
ターミナルは、粒子または減衰粒子の一部の末端が蓄積される位置または関数係数を指す。１つのボクセルおよびサーフィスにつき複数のターミナルがあってもよい。ポインタは、１つのターミナル乗数、またはいくつかのターミナルに対応した複数の乗数を示してもよい。
【００３１】
ターミナル面
ターミナル面は、粒子のタリーを行った局所ボクセル群との境界面である。
【発明を実施するための最良の形態】
【００３２】
モンテカルロ法は、それぞれの離散粒子ヒストリーを正確に追跡し、何億もの（１０億とまではいかないにせよ）正確な粒子ヒストリーを用いて確率的結果を構築する（E. Cashwell & C. Everett, The Practical Manual on The Monte Carlo Method for Random Walk Problems, Pergamon Press (1959)）。本発明は、コンピュータメモリに記憶した離散粒子を、何百万もの異なる位相空間粒子計数値を本質的に表す、精密な位相空間において定義することによって、このプロセスを反転させる。つまり、本発明は、正確かつ効率的に、分布させた位相空間離散粒子のローカルサーフィス、関数係数およびボリュームに対する輸送を計算し、この計算結果を各サーフィス、関数係数およびボリュームに適した乗算フィールドへと還元する。この計算に伴う近似は、離散粒子自体の増分が完全に一定であるという仮定である。その後、「正確な」計算を用いて、ローカル領域において一般化されたフィールドの輸送乗数を決定して連続性を形成し、一般化領域へ拡張する。
【００３３】
サーフィスおよびボリュームから放出される離散粒子は、ＬＶＧ近傍のサーフィス、関数係数およびボリュームに、乗数とボクセルポインタの対を介して直接「布線（wired）」され、それによって粒子輸送のほぼ正確な局所解が得られる（離散粒子そのものの不変性が前提である）。ＬＶＧは、粒子数の寄与率（particle count contribution）をローカル参照ボクセルボリュームまたはサーフィスから外部ボクセルボリュームおよびサーフィスに換算する正確な局所解を与える。このことにより、以前の１次元を分析する不変量埋込み法（imbedded invariants methods that break down past one dimension）とは対照的に、３次元的問題に取り組むために必要な精度が得られる。粒子がローカルシステムに適当に分布することから、ＬＶＧ乗算領域は放射線効果を著しく減少させる。ＬＶＧに分布していない粒子は減衰され、ＬＶＧ表面から放出される。表面から放出した粒子の分布は、表面界面での角度分布を測定するために、直接的タリープロセスまたは最小二乗法での関数デポジション(function deposition)によって明示的に追跡される。同様に、例えば異方性Ｐ5等の高位の粒子散乱を可能にするために、関数係数タリーを複合相互作用モデルとともに使用してもよい。
【００３４】
放射線集合の概念を離散粒子の位相空間に導入することで、ＬＶＧ／ＬＶＧ表面界面で正確な粒子輸送が行われ、精度はさらに向上する。そして、コンピュータメモリシステムを掃引することにより、輸送乗数／ポインタシステムを用いて、離散粒子を抽象的概念上、サーフィスからサーフィス、サーフィスからボリューム、ボリュームからサーフィス、またボリュームからボリュームへと理論的に移動させることが可能になる（図７、ブロック４）。使用可能な相互作用モデル（例えば簡単な核粒子単一エネルギー等方性散乱モデル、多群モデル、異方性散乱モデル）を多数用いて、ボリューム内での粒子間相互作用を調節し、位相空間を適当に調節することにより、粒子輸送の掃引を継続する（図７、ブロック４、５および６）。本発明それ自体を利用して、適当な相互作用モデルを作成することが可能である。
【００３５】
多重衝突輸送から相互作用をデカップリングすることにより、ボクセルを通過する放射線経路に沿った正確かつ直接的、局所的な解析解（exact direct local analytic solutions along ray paths through voxels）が適用可能である。相互作用モデルは、さらなる離散粒子輸送掃引のための新しいボクセル離散粒子ソースを作成する。したがって、適当な離散粒子の正確な輸送解を使用することにより、多くの位相空間粒子を高精度で使用することができる。単一乗数はＬＶＧにおいて使用され、直接的な非確率的結果が極めて高速で得られる。
【００３６】
この方法と、従来のグリーン関数による方法［R. D. Lawrence and J. J. Dorning, A Nodal Green's Function Method for Multidimensional Diffusion Calculations, Nuclear Science and Engineering 76, 218-231(1980)］との根本的な相違点は、グリーン関数は、全時間領域および全散乱相互作用モーメントについて境界値問題の解を与えるため、境界条件に制限されるという点である。１次元または多次元的なグリーン関数応答を解くか否かに関わらず、グリーン関数は、２つの点の間に経時的に発生する全てのイベントを表現する。本発明では、１次元での第１飛行衝突解は、時間とは無関係な真空境界を有する輸送解である。その結果、近似は存在しない。つまり、これは真の輸送解である。外部反復モデルおよび分離相互作用モデルは、散乱相互作用および反応率／過渡変数を計上する。グリーン関数による方法は、散乱を含む全マトリクス応答を作成し、その結果、改良拡散理論等の近似化が行われる。本発明とグリーン関数との共通点は、位置結合である。しかしながら本発明は、別の時間および散乱反復において、散乱要素を選別することによってはるかに高い精度をもたらすことができる。
【００３７】
離散粒子の定義を利用することで、本発明のシステムは、全ての先行技術から完全に差別化される。また、メモリポインタ／乗数対を形成するために放射線集合データを縮小することは、本発明に特有のことである。また、正確なＬＶＧ−ＬＶＧ粒子輸送を拡張するための放射線集合△Rを使用することもまた同様である。
【００３８】
−図面の説明−
以下の図面は、先行技術の課題を解決するために利用することができる本発明およびその使用方法を順次説明するものである。以下に各図面を説明する。
【００３９】
図１
本図は特定の放射線のある粒子が横断する有限インクリメンタル表面を示す簡略図である。ボリュームからサーフィスまでのレイトレーシングにより、最小インクリメンタルサーフィス上に粒子が均等に付着する。サーフィスから放出される放射線は、サーフィスの中心から進行し、情報の一部を失うと考えられる。または、本発明の計算前の変形例においては、放射線の表面分布はサーフィス全体にわたって均一であってもよい。放射線の粒子の角度表示については、入射する全ての放射線集合について、全ての角度が明示的（explicit）であってもよい。あるいは、ボリュームからサーフィスへ移動する粒子を、最小二乗誤差マトリクスから構成される関数係数に挿入して、サーフィスデータ圧縮を行ってもよい（図６の記載を参照）。
【００４０】
本発明を洗練するための別の方法として、単一のサーフィス、サブサーフィスまたは複数のサーフィスにおける詳細な空間粒子分布を決定するために関数係数を使用してもよい。これらの方法は図６の記載でより詳細に説明する。
【００４１】
図２
本図はボクセルおよびグリッドシステムの概念を図示したものである。ボクセルの形状は規則的であってもよいし、図示するように変則的であってもよい。ボクセルは、均質物質であってもよいし、より大きなグリッドシステム内に存在する不均質サブボクセルと境界を接するサーフィスであってもよい。不均質ボクセルは、本発明によって形成され、グリッドシステムに埋め込まれてもよい。この場合には、埋め込まれたボクセル周辺のサーフィスは、そのボクセルと残りのグリッドシステムとの界面として機能する。また、埋め込まれたボクセルは、同様または異なる角度の放射線集合分布を有するサブボクセルで構成されていてもよい（図６の記載を参照）。サブボクセルまたは単一ボクセルをこのように用いる時、このサブシステムが上記の局所ボクセル群ＬＶＧとみなされる。
【００４２】
図３
本図は２Ｄ平面上のサーフィス点から放出される代表放射線を示す。各放射線は、ソース表面に対応するカップリング乗数を用いてトレースされる。図示した例の場合、各放射線トレースについて、適当な統合カーネルが適用されて、サーフィス−位置乗数(surface to position multipliers)が累積される。
【００４３】
図４
本図では、同一角度とサーフィス状態値△ｓおよび△Ωとを有する複数の放射線が、ボクセルシステムを横断する様々の経路に続くサブサーフィスから放出される。これらの経路を組み合わせて、サーフィスから周辺の各ノードまでの単一乗数を形成するが、多くの放射線はサーフィスから放出される。複数の放射線と単一の乗数とを組み合わせることにより、精度を犠牲にすることなく、処理時間が大幅に向上する。通常のボクセルシステムを横断する放射線集合を事前計算することによって、通常システムでの処理速度を向上させてもよい（図７参照）。
【００４４】
図５
本図は局所ボクセル群を形成するボクセル群の境界を示す。この２Ｄ図面においては、内側のボクセル群が外側のボクセル群から完全に分離していることが好ましい。この場合、ＬＶＧ内からトレーシングしてカップリングする放射線は境界表面で終了し、ＬＶＧの外側からトレーシングする放射線は境界で終了する。この分離により、個々のボクセルまたはボクセル群を任意の解像度に変更するための実際的なメカニズムが得られる。メモリが十分にあれば、放射線集合内または外で安定した放射線集合の角度依存性を維持することが好ましい。または、ＬＶＧ境界の内外からの異なる角度放射線集合分布をマッピングしてもよい。さらに別の改良点として、直接関数係数デポジション(direct function coefficient deposition)を用いて一般化された△Ω変換メカニズムを設けてもよい（下記参照）。
【００４５】
図６
本図は、本発明によるポインタ／乗数対およびポインタ／乗数群データの関係を示す。本発明では、ポインタ／乗数表のポインタは、遠隔のボクセル相互作用スコアターミナル、離散表面ターミナルまたは関数係数ターミナル群のいずれかを指す。本発明の好ましい実施形態としては、本図に記載のような参照ハッシュテーブルを使用して、特定の放射線集合をグリッド位置に基づくポインタで参照してもよい。この特殊なタスクを実行するために各種の方法があり、様々のデータ構造を使用することができる。本図には、特定の遠隔ポインタに対応する放射線集合インデックスを含む、好ましい小型のメモリフットプリントデータ構造が記載されている。この好ましい実施形態の値は、放射線集合インデックスを参照する時に、最低限の処理工程で遠隔のボクセル相互作用点、表面および関数を迅速に変更することができる値である。
【００４６】
例えば、特定の位置相互作用点を変更する際は、その点に関連する全てのポインタが特定され、その点を通過する全ての位置からの放射線集合が再計算されて、物質特性変化に影響を与える。放射線集合をターミナルポインタに関連付けるハッシング技術を用いることにより、物質変更処理が高速化する。
【００４７】
本発明のポインタ輸送乗数表は、ボクセル内またはボクセル表面の特定の点から、レイトレーシング法を用いて設定することができる。このようなレイトレーシング法は、２つの形態をとることができ、設定された角度群またはポイント・ツー・ポイント方式によって一般化される。事前計算放射線集合スキーム（図７参照）を用いてこれを行うことも可能である。
【００４８】
ポイント・ツー・ポイント方式および事前計算方式のいずれの場合でも、ボクセル内および表面での均等な点分布を、適当な角度分布Ｐ（Ω）と組み合わせて用いることにより、放射線集合乗数を累積する。図７の説明に記載されているように、この実施例において、角度重量が逆算される。
【００４９】
ボクセルの相互作用点からの前方向レイトレーシングでは、各ボクセル内の相互作用を表す複数の点を用いて輸送乗数を累積する。それぞれの点には、点で表されるボリュームを示す空間重量と位置との両方が与えられる。位置および重量を公正に選択することにより、各点に関連する数学演算が最小になる。
【００５０】
各放射線トレースについて、放射線の固有の方向を確定し、特定の放射線によって表される立体角によって、放射線をさらに重み付けする。一般的には、これは極性の方位球座標で分解されるｗ△Ω＝ｓｉｎ（θ）△θ△φ／４πとして検出される。ｗ△Ωは離散的角度重量を表す。ボクセルシステム通過時の減衰を表す、問題に適した統合カーネルが適用されて、乗数の重みがさらに減らされる。
【００５１】
全ての場合において、輸送乗数は、ボクセルシステム内をトレースするそれぞれの放射線または放射線集合について累積される（図３）。累積のために離散表面および単一ボクセル相互作用点の両方を重み付けする際には、各放射線トレース経路に累積した直接的な単一衝突統合カーネルを用いる。
【００５２】
本発明の好ましい実施形態として、図６の輸送／乗数システムは通常、高速アクセスのためのポインタに関するハッシュテーブルと共に形成される。乗数は放射線集合の関数および開始点として単純に累積される。
【００５３】
関数係数デポジションは、さまざまな様態で形成し得、２つの使用モードを有する本発明の好ましい実施形態である。係数に寄与する各放射線集合について直接的単一衝突統合カーネルおよび適切な関数重量の両方を重み付けすることによって、関数係数への各デポジションを累積する。関数係数デポジションは、事前計算オプションで予め決められていてもよいし、重みは通常、ある任意のグリッドシステムについて１回だけ計算される。
【００５４】
関数係数を累積するための好ましい方法では、係数に寄与する各放射線集合の角度に対応した曲率マトリクスの逆数である最小二乗誤差マトリクス(P. R. Bevington, Data Reduction and Error Analysis for The Physical Sciences, 153頁、McGraw Hill Book Company, (1969) Library of Congress Catalogue number 69-16492参照）を構成する。数値演算では、直行関数が与えられると、最小二乗法により、特定の独立したパラメータ標本点（例えばインデックスされた放射線集合の方向パラメータ）に関する係数マトリクスを作成できることがよく知られている。これは、特定の標本関数点を、特定の係数にデポジットするのに適した重み付けを示している。例えば、公知の球面調和関数に基づく関数の場合を考える。
【００５５】
【数１】
2007531566.tif 000003
【００５６】
式中、Plm（ｃｏｓθ）は対応するルジャンドルの多項式であり、上記多項式のｉは虚数を表す。離散ビンの関数としての粒子タリーを表す角度関数の構成はｆ（θ，φ）＝ΣlΣmＣlmｘＹlm（θ，φ）であり、形式的にｌの合計は０から∞、ｍの合計は−１から＋１である。
【００５７】
事前に知られているように、予想される全ての標本点に依存しないパラメータ（角度または方向余弦の離散放射線集合の値）、係数マトリクスＣlmを簡便形Ｃjに線形化することによって、最小二乗加重マトリクスが構成される。ｉ個の放射線集合に対する以下の係数加重マトリクスが構成される。
【００５８】
【数２】
2007531566.tif 000004
【００５９】
式中、ｗｉｊはｊの各放射線集合方向ｉに適した重みであり、ｘは線形化された係数であって、各関数はｉ点で１からｎまで評価される。輸送スイープにおける関数評価が常に既知点（known point）で行われているため、対称二乗係数マトリクス（symmetric square coefficient matrix）を減じるｋの合計が可能である。得られた加重マトリクスは、関数係数ポインタ／輸送乗数ターミナルＣjに累積された各放射線について輸送乗数を修正するために使用する。これらの線形化係数は、Ｃlmに対応するので、関数ｆ（θ，φ）は、最小二乗フィッティング精度（least squares fitting accuracy）で再構築される。重みｗｉｊは、実際のｆ（θ，φ）の値には依存しないが、むしろ線形化された関数係数に関連するｎ個の輸送乗数に寄与する特定のそれぞれの放射線に関連する公知の標本点θi，φiに依存する。関数は、輸送問題スイープで十分に構築される。
【００６０】
【数３】
2007531566.tif 000005
【００６１】
【数４】
2007531566.tif 000006
【００６２】
式中、ｖは特定のターミナルに寄与する各ボクセル、ｇivは係数が累積されたときのボクセルｖからのタリーへの寄与率である。全ての関連するｎ個の角度の合計は、セットアップ計算（図７、ブロック２または２Ａ）での輸送乗数プロセスの一部として計算された。ボクセルボリューム位置での各ボクセルｖとｔｖ初期タリースコアの合計を、輸送スイープ中に行う（図７、ブロック４）。このように、ある位置で散乱またはソースタリー、およびΣiｇiv*ｗijで表される計算された輸送乗数が与えられると、係数は、演算およびｎ（ｉとの合計）個の明示的角度群から要求される明示的タリー角度をより小さい数の係数に圧縮する（例えば、極端な放射線効果問題の場合、一方の側に４０００個の角度があるが、Ｐ5サーフィス調和近似の場合、係数は３６個のみである）。あるいは、ｎ個の表面から係数、またはその他の対象の関数状態値への圧縮であってもよい。
【００６３】
データ圧縮に関数を用いる時には、それは形式的には離散タリーの関数であることを忘れてはならない。実際の輸送乗数は、明示的な細分化された離散タリー放射線集合構造を用いて表面から進行する。関数は放射線集合角度システムを表面境界で翻訳するか、またはデータ圧縮を行うだけである。
【００６４】
例えば、体球関数の場合には、ユーティリティと第１の使用モードが明らかである。修正Ｐ1散乱にのみ適しているボクセル内の単一相互作用タリーにデポジットするよりも、球面調和関数形式は与えられたいずれの相互作用ボクセルにおける高位の異方性Ｐn散乱を可能にする。この種の球面調和関数を高位の散乱計算に使用することは慣例的である。高位の二重微分散乱データはこのような表現にすぐに適用可能な形式で表れる。
【００６５】
この実施形態を用いる第２の使用モードは、表面についてのデータ圧縮スキームである。これにより、ボクセルから表面への角度デポジション分布を正確に表現するために、何千もの個別の放射線集合アキュムレータをＬＶＧ表面に設ける必要が軽減される。個別に追跡された膨大な数のｉ線群を境界に蓄積しておくよりも、多数の異なるｉ線群角度を集約するＣj係数だけがあればよい。このように、Ｃj係数にデポジションを圧縮する際、通常は、正確な表面形状を構築するための乗算が少なくてすみ、ポインタ乗数対を格納するメモリロケーションが少なくてすむ。
【００６６】
データ圧縮の他にも、この関数形式は、角度群のうちの１つの群を別の表面界面で翻訳することを可能にする。ＬＶＧまたは非均質ボクセルを一般的システムから完全に分離するＬＶＧ表面は、異なる角度離散化スキームを使用することができる。例えば、外側グリッドシステムにポイント・ツー・ポイントシステムを使用し、レイトレーシング法または事前計算モデルを分離したＬＶＧの内側に使用してもよい。また、角度群の順序は異なるが、同様の離散角度計算システムを用いてもよい。例えば、システムの外側に何百もの放射線集合角度を有し、何千もの離散角度を分離ＬＶＧの内側に使用してもよい。関数形式は、機能的補間法（functional interpolation）による境界をまたいだ翻訳を可能にする。あらゆる場合において、関数はそれぞれの放射線集合の立体角方向について、角度タリーを再計算するために用いられる。１組の関数係数を内側から外側へ用いて、別の組の関数係数をＬＶＧサーフィス境界に沿って外側から用いられる。
【００６７】
これは関数係数の使用法の好ましい形式のひとつであるが、これはその他の関数でポジション手法を除外するものではない。上記の方法では、いずれの直行関数を用いてもよく、球面調和関数に関連する表面調和関数や、一般的な直行多項式を用いた特定の実験が行われる。
【００６８】
また、事前計算されたベジエ点とともにＢスプラインを用いてデータ圧縮および翻訳を行うこともできる。データ圧縮精度向上のためにウェーブレットを用いてもよい(Y. Nievergelt Wavelets Made Easy Birkhauser (1999) ISBN 0-817604061参照）。係数の作成は困難であるが、そのアプローチは、係数系に輸送乗数をデポジットする方法のうちのひとつである。このような方法は、タリー関数再構築に関連する係数を縮小するという点において有利である。しかしながら、データ整合性や歴史的・理論的理由から、体球関数または球面調和関数形式が好ましい。これは、整合性のある直行関数で、散乱を簡単に解明することを可能にする二重微分散乱断面積に用いる形式と一致する。
【００６９】
図７
図７は、本発明のアルゴリズムを説明するための一般的なフロー図である。各プロセスまたはフローブロックについて以下に説明する。
【００７０】
Ａ．物理的システムデータベース（図７、ブロック０）
粒子輸送のモデル化には、輸送媒体に関連した物質および幾何学的レイアウトをある程度指定することが要求される。物理的システムのインプットは通常、データベースから得られる。
【００７１】
Ｂ．グリッド構造（図７、ブロック１−１Ａ）
例えば、ボクセルのグリッドシステムについて考える。グリッドは、各ボクセルにおいて、モデル化の対象となる物理的システムを物質組成と重ね合わせる。物理的システムのグラフィックイメージをボクセルに変換することにより、グリッドシステムが形成される。または、ボクセルグリッドシステムの特定のインプットをシステムに入力してもよい（図７、１Ａ）。グリッド構造は物質組成の変化を考慮しなければならない。最も基本的なレベルでは、１個のボクセルが複数の媒体を包含する時は、均質化の何らかの形式を適用しなければならない。これは、単純なボリューム加重スキームから、フラックスボリューム加重スキームにわたって異なっていてもよい。このような方法は当該技術分野において周知である。本発明は、非均質サブボクセルを含む複合ボクセルを形成するために利用することができる。
【００７２】
ＩＭＲＴ  ３ＤＲＴＰにおいて、物理的システムとは、組織および骨だけでなく、金属製の挿入タブおよび人口器官をも備えた３Ｄの人体モデルである。グリッドシステムは、身体画像と、例えばγ線均質化されたマクロ断面積等の、それに含まれる物質特性とを重ね合わせる。このグリッドシステムは、市販の３Ｄグラフィックパッケージ情報を適当なグリッドシステムに変換したものであってもよい。好ましくは、このグリッドシステムは、ボクセル物質を均質化する必要を最小にすることによって精度を最大にするために、骨および組織に関して変則的な形状となっている。
【００７３】
Ｃ．放射線集合ＬＶＧ構成‐概説（図７、ブロック２Ａ−５）
インライン式（図７、２）または事前計算方式（図７、２Ａおよび２Ａ１）で放射線集合を形成することができる。インライン式の放射線集合／ＬＶＧ構成は、変則グリッドに最も適している。ＬＶＧ構成に用いる事前計算された放射線集合の幾何学的特性は、規則的なシステムの高速計算に適している。本発明の好ましい実施形態においては、放射線集合ＬＶＧ構成の全ての形態がコンピュータ内のオプションとして含まれる。事前計算による放射線集合の形状は、アナログコントロールシステムにおいて好ましい。
【００７４】
放射線集合モデリングのコアとなるものは、単一のボクセル内相互作用を各相互作用スイープに使用することである。各放射線に関し、衝突および非衝突粒子密度を用いて、相互作用する粒子の数を測定する。そして、相互作用モデル（図７、５）に従うか、またはＬＶＧ表面境界の横断を続けて、さらに粒子輸送を継続する（図７、４）。単一衝突相互作用モデリングの確率過程で放射線集合を生成してもよい（Cashwellら、supra）。また、適当な立体角ドメインに対する粒子分布の直接的集積化を用いて、適当な放射線集合の幾何学的係数を直接的に計算することもできる。この係数は、事前計算の一部として解析的、半解析的または確率的に導かれる（図７、２Ａ１）。次いで、これらを、特定の放射線、横断したボクセル内の放射線の個々の長さおよび単一相互作用粒子輸送カーネルに関連する、適当な離散角群△Ωの頻度とともに使用して、ＬＶＧボリューム相互作用および表面乗数を計算する（図７、２Ａ）。
【００７５】
あるいは、離散輸送法に類似した方法を用いて、参照ボクセルサーフィスまたはボリュームから放出されてＬＶＧを横断する放射線集合を測定することも可能である（Lockwoodら、supra； Cumber, supra）。本発明の文脈において、このような方法との大きな違いは、ボリュームから放出される放射線が質量中心だけから表現されてはいないことである。より重要には、ＬＶＧを横断する第１の飛行相互作用だけを考慮し、別の相互作用モデルを用いて散乱を処理することが好ましい。本発明は、反復的方法とともに前方向への直接的方法において放射線輸送を考慮することで、粒子散乱に対処する。この違いは、著しい精度の向上に寄与している。ＤＴＭに類似のサーフェス・ツー・サーフェス輸送カップリング（surface-to-surface transport coupling）の代わりに、ポイント・ツー・ポイント法を使用してもよいが、予め決められた放射線集合角度群で純粋なレイトレーシング法を用いると、良好な結果が得られ、より大規模なシステムへ不変量ボクセル群を簡単に埋め込むことができる。さらに、本発明によれば、図６の説明に記載したアプローチのメリットを生かしつつ、一般的な関数係数へ直接的に伝達を行うことができる。最後に、本発明は、１Ｄの放射線集合の結果を輸送乗数に累積することにより、計算性能を大幅に向上させる。
【００７６】
適当な粒子輸送方程式の標準的な直接解析解（standard analytic direct solutions）を、離散粒子ＬＶＧ輸送乗数のインライン式（図７、２）計算に用いてもよい。好ましい方法としては、極めて大型のサンプルサイズの確率的方法を、多数の表面を有する大型の正多面体グリッドシステムに使用して、ＬＶＧに放射線集合の長さを重ね合わせることが挙げられる。変則的なグリッドについては、参照ボクセルボリュームまたはボクセルサーフィス上に、有限表面△ｓおよびサブ表面の中心として、多数の点を設定する方法が好ましい。その後ポイント・ツー・ポイント法で放射線集合を計算する。
【００７７】
ポイント・ツー・ポイント法では、特にボクセルボリューム内のＩＭおよび表面のコサインに関する適当な角度分布が、参照ボクセルサーフィスまたはボリュームからの全ての点集合によって形成される立体角に基づき、参照点からＬＶＧ表面点に移動する放射線の対応部分に因数分解される。
【００７８】
前方向レイトレーシング法では、十分に多数の立体角離散群を使用することにより、表面Ｐ（Ω）の分布仮定（distribution assumption）の必要性が低くなる。
【００７９】
Ｄ．事前計算放射線集合（pre-computed ray set）（図７、ブロック２Ａ１）
事前計算放射線集合オプションについて、頻度およびボクセル長さを計算するためのプロセスを以下の方程式で説明する。サブ表面、サブ角度群で粒子放射線集合に関連する頻度は以下のようにして与えられる。
【００８０】
【数５】
2007531566.tif 000007
【００８１】
式中、ｐ（Ω）は固体群△Ωj内の適用可能なサーフィスを流れる粒子に適した粒子分布である。異なるレベルまたは衝突モーメントのｐ（Ω）が可能であることに留意されたい。近傍の粒子ソース分布ｐ（Ω）は、対象としている相互作用モデルに適している（例えばフラットなボクセルボリューム分布、等方性）。減衰分布ソースから比較的離れたところで、ｐ（Ω）は△Ｓに対するコサインの正規分布を表現する。△Ｓiおよび△Ωjは、適用可能な放射線集合△Rに適した、放射線集合境界が規定されたサーフィス(ray set bounded surfaces)および角度集合を表す。ｆi,j（△Ｓi，△Ωj）は、特定の放射線集合を横断する粒子の適当なインデックス付き頻度を表す。
【００８２】
ここで、ボクセルｌ内の放射線集合△Rの平均放射線長さを以下のように考える。
【００８３】
【数６】
2007531566.tif 000008
【００８４】
式中、△ｒl,i,j（△Ｓi，△Ωj）はボクセルおよび経路に依存する平均放射線長さである。
【００８５】
これらの積分の解は、２次元的埋め尽くしであっても、多くの場合自明でない。これらは多くの場合、△Ｓ、△Ωおよび△R間の複雑な従属性である。したがって、事前計算放射線集合の頻度および長さを解くための最も一般的なアプローチは、モンテカルロ法である。これは確率過程に立ち返ってはいるが、計算はオフラインであり、減衰を伴わないので、物質特性は重要ではないということに留意しなければならない。得られた幾何学的特性を後で物質特性と組み合わせることで、輸送乗数が決定する（図１０）。したがって、何十億もの粒子線が衝突せずに形成され、特性が平均化されて上記積分が効率よく解かれるため、この結果は極めて正確である。この結果は包括的でもあり、物質に関係なくモデル化されたグリッド形状を有する特定のｐ（Ω）に応用可能である。
【００８６】
図８
図８は、図７のブロック２Ａ１に説明したように、本発明において放射線集合を計算するために使用される簡単なモンテカルロアルゴリズムを示す図である。幾何学的な局所ボクセル群がセットアップされる（図８、Ｂ１）。放射線に対応した幾何学的特性のみが生成されるので、ボクセルに関連付けられた物質特性はない。特定のボクセル経路の放射線集合に関連する適当な多重放射線を決定するために、事前計算を行う（図８、Ｂ２−Ｂ５）。
【００８７】
参照ボクセルサーフィスまたはボリュームに適したｐ（Ω）分布を用いてソース粒子を生成する（図８、Ｂ２）。最も簡便な方法としては、カノニカルｕ、ｖおよびｗの方向余弦をデカルト３次元座標に作成する。局所ボクセルボリュームソースｐ（Ω）の分布を用いることは、各放射線集合ライブラリに特異的であり、いくつかの分布は実際の計算に使用してもよい。このような集合は、距離モーメントを表す。第１の距離モーメントはソース（ＩＭソースを含む）からＬＶＧ表面境界まで、第２の距離モーメントはＬＶＧ表面から遠隔表面までというようにである。したがって、図８の出力のいくつかは図７の処理で使用してもよい。
【００８８】
ＩＭＲＴ  ３ＤＲＴＰだけでなく、その他の粒子輸送について、事前計算モードでの本発明の好ましい実施形態では、２組の分布モーメント手段が用いられ、ボクセルボリュームからのソース粒子および散乱粒子と、それらの特定の分布とが１つの組とされる。ＬＶＧ表面を通過して放出される粒子は角度群△Ω内で分布したコサインであり、もう１組のモーメントを形成する。
【００８９】
しかしながら、図６の記載ですでに説明したように、ＬＶＧ境界の一方の側でノードに事前計算セットを使用し、境界のもう一方の側でレイトレーシングまたはポイント・ツー・ポイント放射線集合等の全く異なるスキームを使用してもよい。この種の手段はｐ（Ω）に関する仮定が不完全な場合に好ましい。
【００９０】
ソース粒子は、次に特定のボクセル経路に沿ってＬＶＧ境界に投射される（図８、Ｂ３）。粒子ヒストリーは経路を示す放射線集合に応じて記憶され、限定経路長が決定されて、該放射線集合についての代表放射線が適切に展開される（図８、Ｂ４）。事前計算最小タリーが達成されたのち（図８、Ｂ５）、適当なボクセル経路に依存した、放射線集合を表す放射線長さの範囲が形成される（図８、Ｂ６）。簡単な例では、事前計算での最大放射線長さおよび最小放射線長さの間で、均等なビン分布（bin distribution）に基づき、スコアリングビン（scoring bins）が設定される。放射線集合スコアと、スコアリングビンとを関連付けることは、特定の粒子線集合に関連する特有のハッシュストリングを構築することによって、最も効率的に行うことができる。また、このように定義された放射線集合は、明示的角度群△Ωだけでなく、局所参照システムと、最も重要である表面参照ＬＶＧ近傍システムとの両方について投射されたＬＶＧ境界を示す値を含んでもよい。後者の値は表面でＬＶＧをカップリングするのに用いられる。
【００９１】
ソース粒子、粒子投射およびスコアリング工程が繰り返される。（図８、Ｂ７−Ｂ９）、このとき、事前計算において展開された放射線集合内の代表放射線が、より詳細にスコアリングされる。さらに別の特異な放射線集合の経路が見つかった時は、該集合を構成する複数の放射線を除いてこれらの経路をスコアリングし、後のより詳細な放射線計算のために保存しておく。このことは、実際的な問題として、１０9の事前計算粒子ヒストリーであっても、角度群が多数あり且つ３Ｄ形状が大きい場合にもあてはまることである。
【００９２】
タリー信号または所定のカウントを受けて（図８、Ｂ１０）、事前計算された放射線集合の特性はデータベース、コンピュータファイルまたはその他の記憶媒体に保存される（図８、Ｂ１１）。この工程は、後に最粒度・高解像度の計算を行うために再開してもよい。これにより、放射線集合の幾何学的特性について、極めて高い統計精度が得られる（図８、Ｂ１２）。
【００９３】
図９
図９は、表面から放出されるフラットソースノードに関する図８の工程における、特定の放射線集合の出力を表す。この出力は、本発明のプロトタイプコードを用いて生成した。図９のライン出力を以下に記載する。
【００９４】
ライン１：放射線集合の数
サイズ、角度群および形状等の仕様に依存する数ダースから数十、または数十万の放射線集合があってもよい。これは６3で１２９２７番目の放射線集合であり、８個の角度群システムである。１０は横断したボクセルの数を表す。［０］［０］［０］はｕ、ｖおよびｗのコサイン角度群に関する角度群を表す。
【００９５】
ライン２−４
ライン２は、この特定の放射線集合を区別するためのハッシュコードを反映している。ハッシュコード（または２分木）を用いることにより、放射線集合を追跡するための効率的なメカニズムが得られる。方向余弦は、［０］［０］［０］であり、その後に３３６の２つのエントリが続く。この値は、特異的なＬＶＧ出口サーフィス点座標を表す。ローカル出口サーフィスが、隣接するＬＶＧシステムの基準として使用される場合には、この値は隣接するＬＶＧ出口サーフィス座標と偶然にも同じである。後に続く値は、ボクセルおよび立方体ボクセルの終了面インデックスである１〜６を表す数値の対である（なお、通常の立方体ボクセルは２４、５４またはそれ以上の面を有する）。ライン３および４は、便宜上ＬＶＧ表面の出口を表す。
【００９６】
ライン５
ライン５は、平均経路長（ユニット１の立方体ボクセル(unit 1 cubic voxel)に基づく）を粒子（非使用）の数で割った値を表す。これに放射線集合のうちの最も長い放射線と最も短い放射線が続き、最後にｐ（Ω）からサンプリングされた放射線集合の頻度が続く。
【００９７】
ライン６−１５
ライン６−１５は、詳細な経路を示す。これらのラインでは、[Ｉ][Ｊ][ｋ]座標における平均代表放射線のライン毎の各ボクセル、出現したサイド、各ボクセル内で通過された平均長が順次示されている。
【００９８】
ライン１６
ライン１６は、集合に３つの放射線が存在することを指定する。３パネルインテグレーションと同様の概念である。
【００９９】
ライン１７
ライン１７は、第１パネル、パネル０への情報と、それに続くＬＶＧ境界までの平均放射線長さ、最大長さおよび最小長さを与える。
【０１００】
ライン１８
ライン１８は、上記からの継続であり、放射線集合のうちの放射線の相対頻度を示し、使用される上位長さの間隔がそれに続く。
【０１０１】
ライン１９−２８
ライン１９−２８は、パネルを表す放射線長さを経路に与える。
【０１０２】
ライン２９−４０および４１−５２
ライン２９−４０および４１−５２は、他の２つのパネルである１と２のそれぞれのためのパネル０に上記のシーケンスを繰り返す。このことにより、３つのパネルを表す放射線集合についての情報が完成する。
【０１０３】
モンテカルロ法は一般的に複雑な形状について好ましい方法であるが、図９の事前計算プロセスと同様のデータを得ることができるのであれば、放射線集合頻度およびボクセル横断長さに関する上記の幾何学的積分を解くためのその他の方法を使用することもできる。
【０１０４】
図１０
事前計算による物質特定ＬＶＧの構成（pre-compute material specific LVG construction）（図７、ブロック２Ａ）
図９に示す平均放射線長さの出力が与えられた状態で、図７の一般的なグリッドシステムについての工程１または１Ａに示す物質特性を用いて、特定のＬＶＧ乗数を構成することができる。図１０にこの工程のブロック図を示す。
【０１０５】
参照ボクセル位置を処理するための第１の工程（図１０、Ｂ１）は、伝達アキュムレータの割り当てまたは再割り当てである（図１０、Ｂ２）。アキュムレータは各放射線集合△Rおよび粒子−物質間の相互作用の影響を受ける全ての明示状態の粒子の関数であり、例えばエネルギー△Ｅである。
【０１０６】
次の工程では、参照ボクセル位置のためのターミナルメモリを割り当てる（図１０、Ｂ３）。ターミナルは、表面上の離散粒子またはボクセルボリューム内の相互作用タリーのいずれか、および参照位置からターミナル位置までの輸送乗数を表すボクセルＬＶＧメモリロケーションへのポインタとして定義される。ターミナルは、△Rだけでなく、その他の物質状態変数の関数であってもよい。この割り当てに続いて、ターミナル離散粒子メモリポインタをキーとするハッシュテーブルを作成し、ＬＶＧのターミナル位置を参照することによって乗数値を累積する（図１０、Ｂ４）。
【０１０７】
最後に、参照ボクセル位置から各放射線集合のグリッド位置システムをウォークスルーして、参照ターミナルからＬＶＧターミナルまでの離散粒子乗数を計算して累積する（図１０、Ｂ５）。ハッシュテーブルは、ウォークスルーした位置の乗数を累積するための適当なビンを迅速に指定するのを支援する。伝達アキュムレータは、統合カーネルで各ボクセルへ伝達した値をタリーするために使用される。統合カーネルを放射線集合長さ△ｒl,i,j（△Ｓi，△Ωj）に適用し、放射線を横切る粒子の一部ｆi,j（△Ｓi，△Ωj）に設定された伝達アキュムレータから始めて、ＩＭで用いるボクセルボリュームについての衝突値およびＬＶＧ表面についての伝達値を累積する。参照位置から最終位置までの最終合計が輸送乗数である。
【０１０８】
例えば、ＩＭＲＴでのガンマ線減衰に関して、ボクセル表面の輸送乗数はｆ*Σｅ-△r*μである。合計（Σ）はサーフィスへの経路上の各ボクセルについて行う。変数ｆは上記で定義したように、特定の放射線を横断する部分、△ｒは各ボクセルの経路長さ、μはボクセル物質に応じた適当な減衰率である。例えばエネルギー等の、全ての物質位相空間状態についての乗数が存在する。表面のＬＶＧ点については、放射線集合または角度群を、精度向上のために特定の乗数と共に使用してもよい。位置、放射線集合、角度群および物質状態のそれぞれについての結果は、ポインタ−乗数対の線形配列であり、ポインタはコンピュータメモリ内のＬＶＧ粒子数アドレスを参照する（図６）。ポインタおよび乗数に関するその他の蓄積スキームは、例えば２つの同期配列が考えられる。これらの乗数が設定されていれば、参照離散粒子のＬＶＧ配列をスイープするだけで、粒子がＬＶＧ近傍に輸送される。
【０１０９】
本発明の一実施形態においては、メモリ使用量の観点から、ＬＶＧターミナル位置のスイープに対応するように正しく配列された乗数のベクトルを有することが望ましい。ポインタ・乗数対に替えて、乗数を、このスイープで使用したポインタ演算に対応する単一ベクターに配置してもよい。これは境界を含まない内側のＬＶＧについて最も簡単に実行することができる。
【０１１０】
事前計算した幾何学的特性のＬＶＧ構成プロセスを、グリッドシステムの全ての点において実行する必要がない、ということに留意することが重要である。グリッド内での物質インデックスのパターンマッチングを適用して、同一乗数を使用することができるシステムを特定することができ、単純なポインタ演算を適用して、ＬＶＧ配列値をその他の同一の物質位置に変換することができる。（J. Karkkainen & E. Ukkonen， Two- and Higher-Dimensional Pattern Matching in Optimal Expected Time, 29 SIAM J. COMPUT., 571-589 (1999)、効率的な多次元的パターンマッチングアルゴリズム等）。事前計算したＬＶＧ幾何学的システムを適用可能な規則的グリッドシステムについては、グリッド内の物質組成が高度に分化していない限り、高速計算のためのパターンマッチングのいくつかの形状が好ましい。同様に、パターンマッチングは、規則的なグリッドを使用する場合には、ＩＭＲＴ  ３ＤＲＴＰでの使用に好ましい本発明の実施形態である。
【０１１１】
図１１
Ａ．インライン式放射線集合ＬＶＧの構成（図７、ブロック２）
単一ＬＶＧ全体において、変則グリッドシステム内の参照位置から離散粒子輸送乗数を計算するためには、離散輸送法の放射線レイアウトに類似したポイント・ツー・ポイント方法体系を用いる方法が好ましい。インライン式ＬＶＧ構成の目的は、離散粒子輸送乗数の放射線集合に基づいた計算を直接的にアセンブルすることである。ソース離散粒子の空間的分布は一定であるため、また、その他の離散粒子空間への経路が分かっているため、単位離散ソースに基づいて輸送乗数を計算する際に、いずれの従来の直接的な放射線輸送方法を使用することができる（例えばBellmanら、An Introduction to Invariant Imbedding, SIAM (1992)；A. Shimizu, Development of Angular Eigenvalue Method for Radiation Transport Problems, 37 J. Nuclear Science and Technology, 15-25 (2000)； Olveyら、Accuracy Comparisons for Variational R, t and t-l Response Matrix Formulations, 14 Annals of Nuclear Energy, 203-209 (1987); Sternickら、The Theory & Practice of Intensity Modulated Radiation Therapy, Advanced Medical Publishing, 37-49 (1997)）。これらの構成を一度に使用することで、放射線効果を緩和するために多くの角度群を使用しなければならなくなる。従来のレイトレーシング法はこのようなアプローチにおいて使用することができ、またそれを必要とする時期がある。しかしながら、角度群よりもむしろ放射線集合の概念を維持して、任意または変則的なグリッドシステムにもＬＶＧ表面をリンクすることが好ましい。ポイント・ツー・ポイント法または従来のレイトレーシング法を用いるか否かに関係なく、アルゴリズムは実際的には同じであり、図１１に示されている。
【０１１２】
第１の工程は、メモリを表面に割り当てることである。ポイント・ツー・ポイント法においては、表面中心に対応させて角度分布を計算しなければならない。レイトレーシング法を用いる時は、角度系が予め定義される。いずれの場合においても、ボリューム−サーフィスカップリングのための最初の重み付けを決定するために、放射線によって表される相対的な立体各領域を使用する。表面−表面カップリングは、放射線集合ごとに行われ、関数係数を終止点として用いなければ、特別な重み付けを必要としない。ブロックＢ１以下については、レイトレーシングおよびポイント・ツー・ポイント法は同じである。
【０１１３】
図１１のブロック２において、ボリュームからサーフィスへのレイトレーシングおよびポインタ乗数対の累積が行われる。図１１のＢ３に示す典型的なレイトレーシングルーチンは、適当な統合カーネルを使用するこの工程で用いられる。ボリュームをサーフィスにカップリングしたのち、図１１のブロック４に示す工程で、サーフィスをボリュームにカップリングし、サーフィスをサーフィスにカップリングする。サーフィスにおいて、各放射線集合は、それぞれ個別に追跡されて、更なる伝達のために粒度の高い詳細が提供される。このことは、メモリが許す場合に好ましい。メモリが制約される時は、データを圧縮してメモリ要求を縮小するために、関数係数を表面終了点として用いてもよい。しかしながら、境界で正確な角度情報が失われるので、精度が低下してしまう。
【０１１４】
Ｂ．最初の離散粒子入力（図７、ブロック３）
ここでは、ボクセル境界において、適当な状態変数で指定された離散粒子の初期条件が存在すると考える。離散粒子カウントは、それを放出する全てのボクセル面に及ぶ。同時にまたは代替として、初期条件は、ソース粒子カウントとして与えられる、ボクセルボリュームから放出されるソース粒子の数として得ることもできる。しかしながら、この値はソースボクセル表面での最初の離散粒子カウントに変換することができる。ＩＭＲＴ  ３ＤＲＴＰにおいて、組織内で散乱放射をモデリングするための初期条件は、最初の直接的放射線計算（図７、Ｂ９Ａ）から進む。相互作用率と、場合によっては角度情報とは、ボクセル内での１回目の放射線衝突のために記録される。この値は、最初のガンマ線離散粒子カウントを構成するために使用する。また、本発明は、表面離散粒子カウントおよびコサイン距離モーメントとともに用いることもできる。しかしながら、典型的な放射線法を用いてシステムに入射する実際の放射線を明示的にモデリングする方が好ましい。本発明の相互作用モデルは、さらなる粒子輸送のソースとなる散乱放射をモデリングするために用いることができる。図６で説明した関数係数デポジションは、高位の散乱相互作用の処理に用いることができる。
【０１１５】
Ｃ．離散粒子輸送スイープ（図７、ブロック４）
離散粒子に関する初期ソース条件が指定されると、ブロック３またはブロック９Ａのいずれかにおいて、輸送スイープを行う。粒子スイープ中、粒子はボクセル離散粒子タリーへ輸送されて相互作用モデル計算が行われ、さらにＬＶＧ表面境界へ輸送される。
【０１１６】
これは、ゼロ以外のカウントを基準とする各離散粒子位置に対する簡単なスイープである。それぞれの基準において、ＬＶＧターミナル−乗数対の直線系をスイープし、各乗数を参照離散粒子カウントに当てはめて、端数の（fractional）粒子カウントをターミナル離散粒子ポインタ位置に累積する。
【０１１７】
この計算は、最も相互作用しない粒子がシステムから一掃された内部収束の状態になるまで行われる。離散時間状態エポック（discrete time state epoch）△ｔを無視することができない過渡的な問題については、異なる内部収束方法が必要となる場合がある。また、このような方法は、ＩＭに対応した計算コストによっては好ましい場合がある。しかしながら、輸送スイープは参照ＬＶＧについて行われ、内部収束は生じない。この場合、粒子がＬＶＧ境界に輸送された直後にＩＭを適用する。ＩＭＲＴ  ３ＤＲＴＰでは、局所ＬＶＧターミナルをスイープし、ＩＭを１つの反復スイープシステムとして用いてさらに散乱を計算することが好ましい。
【０１１８】
図１２
Ａ．相互作用モデル（図７、ブロック５）
相互作用モデルは、適当な状態変数とともにターミナル離散粒子タリーを受け、モデルタイプの選択に応じて、ボクセル表面に、またはボクセル内部から、新しい離散粒子を作成する。本発明は、比較的大きなボクセルのための相互作用モデルを作成することができる。これは、システムの好ましい例である。この方法体系を、有効な相互作用モデルを作成するための簡単な衝突確率アプローチによって以下に説明する。以前から、複合衝突確率法（complex collision probability methods）が用いられている。Marleauらは、中性子計算にこの方法を使用する例を挙げている（Analytic Reductions for Transmission and Leakage Probabilities in Finite Tubes and Hexahedra,104 Nucl. Sci. & Eng．， 209-216 (1990); A Transport Method for  Treating  Three Dimensional Lattices of Heterogeneous Cells, 101 Nucl. Sci. & Eng., 217-225 (1989)。これらは、本発明に関連して、一般化された粒子輸送に容易に適応可能である。
【０１１９】
角度粒子分布を表す球面調和関数を作成するために、上記の関数係数デポジションを用いることができる。これを高位の二重微分断面積と交替で使用することにより、詳細な角度散乱情報を計算することができる。
【０１２０】
相互作用モデルの目的は、衝突後の粒子配置を計算することである。これには、衝突後に相互作用する一次粒子の衝突パラメータ、状態変化および新しい状態変数を有する二次粒子の発生が含まれる。光学装置類、レーダ、ソナーおよび放射熱伝達には、一次粒子の衝突後パラメータの比較的簡単な計算が必要である。核放射線プロセスでは、例えばコンプトン散乱からの反跳電子または核分裂プロセスからの別の中性子といった二次的な放射を、相互作用モジュールによって発生させる必要がある。ＩＭＲＴ  ３ＤＲＴＰでは、ガンマ相互作用によって、低γエネルギーでの光電吸収が生じ、光電子が生成される。低エネルギー散乱光子および二次的な反跳電子の両方を用いたコンプトン散乱は、中位のγエネルギーでモデル化される。１．０２ＭｅＶ以上の極めて高いγエネルギーでは、対の生成プロセスをモデリングすることもできる。
【０１２１】
クリティカルシステムで核分裂を伴うプロセスについては、それぞれの生成無限増倍率を、全システムの固有値に乗じてローカルソース強度を決定するための特別なモデルが要求される。１つの漏れないパラメータ（non-leakage parameter）を各ボクセルに使用することができるとすると、これは図示するように自明なモデルである。
【０１２２】
Ｂ．放射粒子に関する簡単な衝突確率相互作用モデル
上記のように、相互作用モデルに簡単な衝突パラメータによるアプローチを用いることができる。このアプローチでは、マクロ断面積の比を利用して、ボクセル内で衝突する粒子の配置を決定する。このようなモデルを適用するためには、ボクセルボリュームにつき平均１つ以下の衝突が必要である。放射輸送の場合は、理想的には、状態平均自由工程１／Σ＞ｄcを限定することによってこの基準を満たすことができる。Σは最小実体状態の合計または輸送マクロ断面積であり、ｄcはボクセルを横断する最大経路長さである。しかし、実際の実験によれば、このような基準は有意に緩和することが可能である（図１５、１Ａｉｉの散乱結果および説明を参照）。
【０１２３】
この簡単なモデルでは、漏れない確率が計算され、それぞれのボクセルから放出される全ての散乱粒子および二次粒子に適用される。この確率は、フラットな分布を想定し、統合カーネル減衰で抜け出る粒子を計算することによって得られる。文献に記載されるその他の方法を使用することもまた可能である。
【０１２４】
このモデルにおいて収束を高速化するために、ボクセル内で相互作用する粒子の連続する世代に、この同様の漏れない確率を適用する。図１２にこのモデルの簡単なフロー図を示す。このモデルへの入力は、図７のブロック４に記載されたものと、図７のブロック９Ａからの推定される初期条件入力を含む。グリッドの初期化は、相互作用モデルに有用なボクセルについての物質パラメータの事前計算を包含していてもよい。これは、それぞれの状態（例えばエネルギー）、粒子吸収フラクション、状態輸送フラクションおよび伝達フラクションについての漏れない確率を含む。
【０１２５】
漏れない確率と伝達とに関して、２つの取りうるアプローチがある。第１のアプローチでは、全ての散乱粒子が内側および外側に向かってボクセル表面に出現すると仮定する。好ましいアプローチは、漏れない確率ｐnlを用いて、ボクセル内でのさらなる内部散乱を計上することである。ここで、簡単なエネルギー輸送モデルにおいて、状態エネルギー群ｇからｇ’へと散乱する粒子を考える。Σｒgg'を用いて群ｇからｇ’への散乱／除去マクロ断面積を表し、Σｔgを用いてボクセル物質の全マクロ断面積を表す。散乱する粒子（したがって衝突としてタリーされる）は、その最初の衝突について、群ｇ’への散乱の確率Σｒgg'／Σｔgを有する。しかしながら、粒子は群内で散乱する確率Σｒgg／Σｔgを有する。次の衝突については、
次のｇ’へのボクセル内散乱移動Σｒgg／Σｔg*ｐnl*Σｒgg'／Σｔg、それに続いて、後続の各世代については、Σｒgg／Σｔg2*ｐｎｌ2*Σｒgg'／Σｔg等と続く。
【０１２６】
ｐnlが１．０より小さく、また全ての輸送確率が１．０より小さいので、このようなボクセルにおいて群ｇからｇ’へ移動した全ての世代の粒子の合計数は、次の式で与えられるということを容易に実証することができる。
【０１２７】
【数７】
2007531566.tif 000009
【０１２８】
時間固有値の問題については、ｐnlを固有値問題λで乗算する必要がある。核分裂の問題については、群依存無限乗数係数（ｋ?）が除去粒子に対する生成粒子の比を表す。サブクリティカルボクセルについて上記したように、これは散乱モーメントを置換または増加させることができる。核分裂を反映させるための様々のスキームが利用可能である。
【０１２９】
あるいは、異なる値のｐnlが異なる衝突モーメントに適するように、散乱または核分裂生成アプローチをｐnlに利用することもできる。しかしながら、１／Σ＞ｄｃの衝突基準は、散乱伝達モーメントをボクセルボリューム相互作用モデルで使用する必要を制限するためには十分である。
【０１３０】
第１のインラインプロセスを続行して、ボクセル位置を設定する（図１２、Ｂ１）。その後、所定の状態タリーが得られる（図１２、Ｂ２）。相互作用タリーは、上記したように、全ての離散粒子状態値の関数であってもよい。それぞれの状態離散粒子について、相互作用パラメータが適用される。最初のボクセル伝達プロセスは、離散粒子相互作用タリーが状態値と共に与えられると、散乱プロセスを考慮する（図１２、Ｂ３）。Ｐn散乱では、相互作用タリーによってスコアリングされた入射放射線の角度群を、ボクセル離散粒子応答において、適当な散乱ｐ（Ω）を用いて考慮する。その他の離散表面を介したボクセルの離散粒子伝達を決定する際には、タリー入射面△Ｓiもまた考慮される。しかしながら、本発明の好ましい実施形態、つまり等方性散乱を用いる簡単で効率的な単一衝突ボクセルモデルでは、最初の相互作用モーメントのみがこの細部（detail）で処理される。これは、多重散乱の複雑性および後続する散乱の統計的な数が比較的小さいことによる。
【０１３１】
高位の散乱相互作用においては、関数係数によるアプローチを使用して、角粒子分布を決定する。これは、複雑なボクセル内散乱を計算するために使用することもできる。しかし、このような複雑性は、ほとんどのモデリングタスクに必要ではないため、このような複雑な演算を扱うよりも、より小さいボクセルサイズを用いる方が望ましい。
【０１３２】
その後、必要に応じて、ボクセル伝達から、例えば吸収、散乱およびエネルギー付与（線量）等の最終的なボクセルボリュームタリーに進行する（図１２、Ｂ４）。このプロセスに続いて、その他の全ての状態について、ボクセルボリュームタリーに対するボクセル内散乱寄与率を決定する（図１２、Ｂ５）。その後、ボクセル内散乱離散粒子の二次輸送を計算する（図１２、Ｂ６）。最後に、全ての粒子状態および位置を必要に応じてループしてＩＭスイープを計算する。
【０１３３】
３ＤＲＴＰ  ＩＭＲＴでは、散乱を利用する上記モデルは、簡単で好ましい放射モデルである。しかしながら、このモデルは、ここに記載するように、より大きいボクセルボリュームでの埋め込み不変量を決定するためのその他の計算にも使用することができる。
【０１３４】
Ｃ．その他の粒子衝突確率相互モデル
使用方法にもよるが、当該分野において公知である様々な方法によって、ボクセルボリューム相互作用を事前計算することができる。しかしながら、後続のボクセル離散粒子伝達では、これを明示的な粒子状態値として用いる場合には、適当な放射線集合状態においてボクセルを離れる粒子を正確に配置する必要がある。
【０１３５】
Ｄ．本発明に基づく相互作用モデルの作成
本発明の好ましい実施形態において、ボクセルボリューム相互作用モデリングに使用する埋め込み不変量を作成するために、発明そのものを使用することができる。したがって、小さいボクセルに限定されることがなく、上記で定義した基準１／Σ＞ｄcが満たされる。ＬＶＧの放射線集合特性を用いて、この演算をオフラインで事前計算として行うか、あるいは物質およびグリッド特性が設定された後にインラインで行う。いずれの場合においても、ボクセル群に初期境界条件を設定するだけでよい。しかしながらこのモデルでは、大きいボクセルを作成するために組み合わされた全てのボクセルは参照ＬＶＧ内になければならない。
【０１３６】
グリッドシステムをさらに分割し、衝突モーメント応答を明示的に追跡するグリッドシステム全域での離散粒子輸送を解く。サーフィス−ボリューム情報はＬＶＧをより大きなシステムに埋め込む際に使用するために記憶しておく。いくつかの表面を組み合わせて大きなボクセル入射面および出口表面を作成してもよい。離散粒子密度はこれらの表面に均等に分布されて外部入射粒子への適当な統合ボクセルシステム応答を形成する。サブ表面の中心からレイトレーシングまたはポイント・ツー・ポイント法を利用して、埋め込まれたシステム内の輸送を測定してもよい。あるいは、事前計算法を用いてもよい。
【０１３７】
入射および放出に空のボクセルを用いて、放射線集合の初期条件および放射線集合に基づくスコアリングを可能にすることもできる。空所は、放射線集合の距離モーメントのグループ分けを行う。このモードでは、ＬＶＧでのモデリングに必要な、明示的角度群の数を減らすためにこれらを用いることができる。入力側では、遠隔モーメントＬＶＧ放射線集合に対応する粒子流動を生じさせる。出力側では、これらは終了する放射線修後を計算するために使用される。全ての崩壊していない状態のグループについての初期離散粒子のグループ分けを明示的にモデル化する必要がある（図７、ブロックＢ９Ａ参照）。
【０１３８】
核分裂を伴う問題については、衝突時間の関数としての明示的システム応答を使用する必要がある。少なくとも１つまたは通常はいくつかの明示的な発生応答が測定されたのち、収束した無限システム応答を測定してもよい。どのような場合においても、核分裂システムでのボクセル群はローカルクリティカルシステムを形成する（Ｂｅｌｌｍａｎ，ｓｕｐｒａ）。最後に、グループ分けされたボクセルのデータは後の事前計算の一部としての計算または終了計算のいずれかに使用するために記憶しておく。
【０１３９】
Ｅ．収束（図７、ブロック６）
本発明では、非核分裂性粒子輸送の境界値問題とは対照的に、衝撃性初期値問題として粒子輸送を解く。核分裂性粒子輸送については、世代間固有値（generational Eigenvalue）は、世代間核分裂変化に基づいて計算することができる。固有値を、関連する反応率基準と組み合わせると、許容できる解への収束が成立する。
【０１４０】
ボクセルボリュームを基に残留相互作用タリーの全相互作用タリーに対する比を計算するとともに、初期離散粒子入力に関連する絶対的なシステム内粒子カウントを、収束を測定するために用いることができる。収束に続いて、スイープの一部としてシステムから外へ輸送されなかった残留散乱または二次ＩＭ粒子を反映させるために、結果の繰り込みを行う。これは、ボクセルの不変量応答を測定してより広範な問題の解に取り込むような問題においては特に重要である。
【０１４１】
Ｆ．結果のデータベース保存（図７、ブロック７）
計算結果は、従来のリレーショナルデータベースまたはオブジェクトデータベースに保存することが好ましい。放射線集合データも同様に保存することができる。これは本発明による好ましい保存方法である。
【０１４２】
Ｇ．最適化／設計エンジン（図７、ブロック８Ａ）
上記のように、本発明は、自動設計および／または治療計画策定の最適化に理想的である。ブロック８Ａに本発明の結果に基づく設計の最適化を示すが、正確な仕様は本発明の範囲外である。
【０１４３】
ＩＭＲＴ  ３ＤＲＴＰでは、対象となる腫瘍への線量を最大にして、健康な組織への線量を最小にするために、本発明の保存結果を利用し、シミュレートされた適当な外部放射線を発生させる。本発明を用いて調査のための初期計算を行うことにより、腫瘍に有意に寄与しない放射線を迅速に除外することができる。また、本発明は対象外の健康な組織に対する散乱放射寄与度をモデル化するのに理想的である。本発明の計算結果を用いて最適なビーム形状および粒子強度を選択するために、市販の最適化エンジンを利用してもよい。
【０１４４】
Ｈ．初期粒子分布（図７、ブロック９Ａ）
上記したように、最適化エンジンの結果によって試験粒子の初期分布が特定される。放射線ペンシルビームの場合に好ましいモデル化の手順では、散乱放射に対応する初期粒子分布を測定するために、放射線を直接モデル化して相互作用モデルを使用する（図７、ブロック５）。このような初期の放射線ビームはグリッドシステムの認識が必須である。これは、ある特別な場合のＬＶＧの再構築を必然的に伴う。ペンシルビームは代表的な放射線を用いてモデル化してもよい。ＩＭＲＴ  ３ＤＲＴＰでは、散乱放射線効果を計算するために本発明を用いた直接的な放射線ビーム計算が、対象組織へ直接的に放射線をドーズする際に好ましい。
【０１４５】
図１３
本図は、事前計算オプションを用いた規則的な形状のための予備問題セットアップを示す。この問題に使用するプロトタイプは、図７のブロック１Ａ、２ＡＩ、２Ａ、３、４、６および７を用いた。図１４にこの問題セットアップのための特定の結果を示す。本発明のＩＭ解釈を実証するその他のプロトタイプは以下の図面に示す。
【０１４６】
図１３の予備問題では、粒子は、オフコサインのソース分布で以て、影をつけたダクトの手前側入口を通過する。さもなければ完全にシールドされているサイドを通過する特定のソース分布は、それ自体が減衰を含まない１０×１０×１０ｃｍ3の隣接するソースボクセル全体に亘って均一に分布した等方性粒子（isotropic particles）のソース分布を有していた。ダクトシステムは、高さと幅がそれぞれ１０ｃｍの断面を有し、６０×６０×６０ｃｍ3の系に拡張している。放射線効果誤差を最大にするために、その他のソース粒子は境界から流動せず、散乱はモデル化されない。全ての面の境界条件は完全真空である。０．１ｃｍ-1の標準的な全粒子減衰断面を用いる。これは問題軸のうち３つに関して粒子の反射がないことを除き、国際ベンチマーク問題と同様である。
【０１４７】
この問題セットアップでの選択された結果の比較を図１４に示す。
【０１４８】
図１４
予備的平面相互作用率の結果(preliminary planer interaction rate results)。この図は、高度粒子計数モンテカルロ法（上位値）と比較した、本発明の結果（中央値）を示す。さらに、各１０3セルについて、相対差分を百分率で示している（低セル百分率値）。結果はソース平面から４０〜５０ｃｍ上方の平面で表す。０．０相互作用の１つのセルが空のダクトを示すのに対し、影をつけたセルはその平面についての最大誤差を示す。典型的には、直接解は、他の先行技術の直接法の場合、それらの距離において、２０％から４０％も減少している。信号は、極限放射線効果ストリーミング問題(extreme ray-effect streaming problem)について、これらの距離において、数万のオリジナルソースである。また、本発明は、この問題またはその他の問題について、モンテカルロ法よりも１０００倍以上の速度で複数のソース分布を処理するので、設計問題および３ＤＲＴＰ  ＩＭＲＴにとって理想的であることに留意されたい。
【０１４９】
使用したモンテカルロコードは本発明者が構築したものであり、同一のコンピュータで基本的なケースの計算に用いると、本発明よりも１０００倍遅いことがわかった。図１５、１６および１７に、インライン式放射線オプションを用いた標準国際参照ベンチマークに対する結果を示す。
【０１５０】
図１５
図１５に標準国際ベンチマーク問題に対する本発明の現在のプロトタイプ結果を示す。このプロトタイプは、簡単な相互作用モデルを含み、様々な直平行六面体ボクセルを使用する。図７のブロック１Ａ、２、４、５および６で十分に使用したアルゴリズムを使用した。
【０１５１】
参照問題および結果の比較は、雑誌"Progress in Nuclear Energy", Volume39 number 2 ISSN 0149-1970 (2001)の特集号に公表された"3D Radiation Transport Benchmarks for Simple Geometries with Void Region"から転載する。このベンチマークからモデル化した問題は、問題番号１である。
【０１５２】
この問題は、１００×１００×１００ｃｍ3の中心空隙を有する２００×２００×２００ｃｍ3の暗色吸収物質の立方体からなる。空隙の中心には２０×２０×２０ｃｍ3の暗色物質からなるソースが存在する。この問題は２つのモードで解決される。空隙領域の全マクロ断面積は１０-4ｃｍ-1であるのに対し、暗色吸収体の断面積は０．１ｃｍ-1である。この問題は、直接法では極めて難解である。なぜなら物質が暗色で、散乱がほとんどまたは全くなく、問題サイズが使用する断面積には大きいからである。
【０１５３】
この問題は２つのモードで解決される。第１のモードである問題１Ａｉでは、両方の領域が純粋な吸収体である。第２のモードである問題１Ａｉｉでは、両方の領域が５０％散乱を生じるので、空隙領域および暗色領域の両方において、吸収断面積および散乱断面積はそれぞれ０．５×１０-4ｃｍ-1、０．０５ｃｍ-1である。中心ブロックでのソース率は、均等に１粒子ｃｍ-3−ｓ-1である。１つの参照軸を比較のために与える。座標系は各方向に−１００ｃｍから＋１００ｃｍまで延びている。
【０１５４】
本発明と比較して、ＴＯＲＴ、ＡＲＤＲＡおよびＥＶＥＮＴ等の核輸送コードを考慮する。ＰｅｎｎＴｒａｎ等のその他のコードは正確な数値を公表していないが、雑誌に記載の図形から、いかなる場合でも本発明の結果の方が優れているように思われる。正確な解析フラックスまたはＧＶＭＰモンテカルロコード（ＭＣＮＰの変形）のいずれかを比較に用いた。モンテカルロコードは１Ａｉｉの結果を得るのに３７８，０００秒かかっていた（図１７参照）。
【０１５５】
本発明は、フラックスを直接計算しないので、２×２×２ｃｍ3の小さなノードサイズを用いてフラックス率を再構築した。これは、他のコードによって計算されたポイントフラックスと比較して、ノード平均フラックスがレポートされる際に、本発明にとって、誤差のさらなる原因となる。
【０１５６】
図１５の散乱について、本発明では、システム全体が完全にカップリングされた。散乱問題は、ノード軸について効果的な真空境界条件がある場合、より少ない数のノードのモデリングを要求した。散乱問題は、問題１において、全てのノードの完全なモデリングを要求した。ノードサイズは、最小の２×２×２ｃｍ3のノードから、測定軸から離れた２０×２０×２０ｃｍ3のノードへと変化した。このことは、他のコードに用いる方法と矛盾しなかった。これらの特定の結果に対してレイトレーシングを用いた。これは、９９７８個の立体角のモデリングを要求した。
【０１５７】
本発明は、特別のセットアップを用い且つ実行モードが分離された状態で実行された。このセットアップは、任意のソース分布が与えられた場合、実行時間が元の時間のごく一部になるように完了された。図１６の結果は、参照ベンチマーク問題を２つに分割するＬＶＧアプローチと共に本発明を説明している。図１７は、機械時間（machine time）の比較を提供する。
【０１５８】
図１６
この図面は、問題１ＡｉにおけるＬＶＧサーフィスカットの影響を示す。この問題が散乱を有さない場合、サーフィスが、モデリングに際して特に問題となる。選択されたサーフィスは、ボイド／吸収体界面から５０ｃｍの位置であった。３つの異なるサーフィス結果を示す。第１のサーフィス結果は、サーフィスカットを提示していない。第２のサーフィス結果は、各サイド（プロトタイプへの入力）につき４つのサブサーフィスを提示しており、放射線集合は該サーフィスを介して明示的に追跡される。サーフィスカットは、１つのカットサイドにつき、１つの６次のサーフィス調和関数計数フィットを用いた。カットした場合についてカット後の結果を示す。なお、カット前の結果は同一である。
【０１５９】
ＬＶＧサーフィスおよび４つのサブサーフィスのみの場合、妥当な結果が得られた。図６で説明した技術を用いて決定した５７個の係数を６次のサーフィス調和関数に適用すると、同じく良好な一致（agreement）が得られる。追加のサブサーフィスを用いて結果をさらに改善し得る。サーフィス調和関数は以下の通りであった。
【０１６０】
【数８】
2007531566.tif 000010
【０１６１】
上記式中、ｍの総和は１から６までの値であり、ｎの総和は１からｍまでの値であり、Ｐm（μ）はルジャンドル多項式であり、Ｐmn（μ）は関連するルジャンドル多項式である。該サーフィスに対して垂直なコサインがμによって与えられる。φは方位角である。係数ａ、ｂ、ｃは、図６に示した方法に従って線形化され、フィットされた。
【０１６２】
図１７
この図面は、図１５および図１６についてのタイミング結果の比較を示す。本発明を、高価ではないＰＣプロセッサ上で実行した。本発明の装置のクロック速度は他の場合よりも速かったため、タイミング比較が難しかった。本発明のセットアップ時間は、任意のソース分布最適化問題（３Ｄ ＩＭＲＴにおける問題など）について、一回のコスト（one-time cost）である。そういうわけで、この計算時間は、システム全体をカップリングするために１度使用される。このカップリングに続いて、実行時間が提示され、厳密に収束した散乱問題についても、この実行時間は、セットアップ時間よりもかなり速い。本発明についてのセットアップ時間も、比較した直説法のセットアップ時間よりも優れており、モンテカルロ法よりもかなり速い。
【０１６３】
このデータ並びに図１４の結果は、本発明を用いることにより、速度が１０００倍向上することを示している。
【０１６４】
上で説明した本発明は、本発明の好適な実施形態を参照して具体的に示し且つ説明したが、本発明の範囲内において形態および詳細についてのさまざまな改変が可能であり、本発明の趣旨および範囲から逸れることなくさまざまなアプリケーションを用い得ることを当業者は理解する。
【図面の簡単な説明】
【０１６５】
【図１】図１は有限領域△ｓのサーフィスを通過する実際の１個の粒子を示す簡略図である。
【図２】図２は湾曲した矩形領域を有する３×３×３配列のボクセルからコンパイルされた本発明のグリッドの一例を示す図である。
【図３】図３は基準面から２次元グリッドシステムを横断する放射線のサブセットを示す図である。
【図４】図４は有限領域のサーフィスから放出されて複数のボクセルを横断し、特定の放射線集合の立体角群△Ωを占有する放射線の集合を示す図である。
【図５】図５は局所ボクセル群（ＬＶＧ）の基準面およびそれを含むグリッドシステムを示す。
【図６】図６は輸送乗数およびポインタを格納するために使用されるタリーメモリ装置のデータ図である。
【図７】図７は好ましい論理プロセス経路および代替経路とともに本発明の好ましい実施形態を示すブロックボディ図である。図中、各ブロックは媒質中の粒子輸送をシミュレートするプロセスにおける工程またはコンピュータソフトウェア指令を表す。
【図８】図８は本発明において放射線集合の計算的シミュレーションを行うための論理的工程を説明するためのブロックボディ図である。
【図９】図９は図８のブロックボディ図に説明した工程を用いて得られた特定の放射線集合のシミュレーション出力を示す。
【図１０】図１０は本発明のＬＶＧ乗数を決定するために必要な工程を順次説明するための別のブロックボディ図である。
【図１１】図１１はインライン放射線集合に基づくＬＶＧ離散粒子輸送乗数（図７のブロック２実施形態）を示す図である。図１１はプロセスおよび／またはオブジェクトブロックを示す。
【図１２】図１２は相互作用モデルで粒子をシミュレートするための工程を示す簡単なフロー図である。
【図１３】図１３は図１４の課題レイアウトのサンプルを示す図である。
【図１４】図１４は本発明および従来技術によってシミュレートされた放射線輸送の出力結果を示す表である。両システムの精度は同程度であったが、本発明のシミュレーション計算速度は約１０００倍であった。
【図１５】図１５は"3D Radiation Transport Benchmark Problems and Results for Simple Geometries with Void Region", ISSN 0149-1970 "Progress in Nuclear Energy"に定められた国際標準ベンチマークの課題１Ａにおける放射線輸送のシミュレーション結果を比較して示す。本発明とその他の計算コードとを比較した結果、本発明では精度の向上が見られた。
【図１６】図１６は、タリーデポジションサーフィス（tally deposition surface）の背後で結果を正確に予測するための、問題１Ａの結果と本発明の能力を説明するための切断面との比較を示す。これは埋め込まれた問題をシミュレートする本発明の機能を説明する。
【図１７】図１７は図１５および図１６に示す問題１Ａの結果の解に関する時間とマシンとの比較を示す。この時間比較は、本発明の速度が１０００倍であるという利点を表している。
【図１】
2007531566.tif 000011
【図２】
2007531566.tif 000012
【図３】
2007531566.tif 000013
【図４】
2007531566.tif 000014
【図５】
2007531566.tif 000015
【図６】
2007531566.tif 000016
【図７】
2007531566.tif 000017
【図８】
2007531566.tif 000018
【図９】
2007531566.tif 000019
【図１０】
2007531566.tif 000020
【図１１】
2007531566.tif 000021
【図１２】
2007531566.tif 000022
【図１３】
2007531566.tif 000023
【図１４】
2007531566.tif 000024
【図１５】
2007531566.tif 000025
【図１６】
2007531566.tif 000026
【図１７】
2007531566.tif 000027