JP2009069930A - 粒子法シミュレーションのためのスライスデータ構造、およびスライスデータ構造を利用した粒子法シミュレーションのｇｐｕへの実装方法 - Google Patents

Abstract

【課題】粒子を用いた流体計算手法に基づく粒子法シミュレーションで近傍探索に用いる計算格子の用い方につきメモリの使用効率を良好にすることができる粒子法シミュレーションのためのスライスデータ構造、およびＧＰＵへの実装方法を提供する。
【解決手段】ＣＰＵやＧＰＵ１２を用いた粒子法シミュレーションのためのスライスデータ構造は、計算空間３１についてのデータ構造であって、空間は多数のボクセル３２で構成される３次元計算空間であり、Ｙ軸に直交する複数のスライス３１ａを形成して、多数のボクセルを複数の２次元スライスで分割し、複数の２次元スライスの各々で、粒子が存するボクセルの存在範囲について最大と最小のボクセルの各々の開始座標を計算し、存在範囲を矩形の形状で囲むバウンディングボックスを決定し、複数の２次元スライス領域の各々のバウンディングボックスに含まれるボクセルについてメモリを確保する。
【選択図】図４

Description

本発明は、粒子を用いた流体計算手法に基づく粒子法シミュレーションで、近傍探索に用いる計算格子の用い方についてメモリ効率を良好にしかつ近傍粒子探索を効率化する粒子法シミュレーションのためのスライスデータ構造、およびスライスデータ構造を利用した粒子法シミュレーションのＧＰＵへの実装方法に関する。

「粒子法シミュレーション」とは、対象とする物体を多数の粒子の集まり（集合）としてとらえ、多数の粒子の個々の挙動に注目して計算を行うシミュレーション手法である。コンピュータでの計算上、上記粒子には、通常、直径が一定である「球」が用いられる。物体の形状を球の集合体によって表現するとき、当該球の大きさ（直径）の決め方に依存して物体の形状の近似度または空間解像度が変化する。また、球の直径を変化させ、空間解像度を変化させることにより、コンピュータによるシミュレーション計算における計算精度と計算速度を制御することが可能となる。

近年、上記の粒子法シミュレーションを利用して流体等の挙動計算が行われ、コンピュータの表示装置の画面上で流動等の流体シミュレーションの研究が行われている（例えば非特許文献１を参照）。粒子法シミュレーション技術の研究は、コンピュータグラフィックス（ＣＧ）の基礎となる技術として重要である。

粒子法を用いた流体シミュレーションの手法にはＭＰＳ法（Moving Particle Semi-implicit Method）とＳＰＨ（Smoothed Particle Hydrodynamics）がある。ＭＰＳ法は連続の式から導出するポワソン方程式を解いて非圧縮流れを計算するが、ＳＰＨは非圧縮条件を課さず、疑似非圧縮流れを計算する。そのためＭＰＳ法はＳＰＨに比較して計算コストが高い。このため、コンピュータグラフィックスの分野ではＳＰＨの手法が主流となっている。

上記の粒子法シミュレーションの研究では、従来、主にＣＰＵを利用して計算アルゴリズムを作る方法は多く行われていた。しかし、近年では、ＧＰＵ（Graphics Processing Unit：一般的にストリーミングプロセッサ）を利用した粒子法シミュレーションの研究も行われつつある。ＧＰＵは、グラフィックスの処理に特化するために、内部に複数のプロセッサが並列に演算処理を行うという構成になっている。このため、頂点の座標変換等の処理であれば、ＧＰＵはＣＰＵに比較して極めて高速に処理することができる。またＧＰＵは、シェーダと呼ばれるプログラムを書くことによって処理を制御できるため、グラフィックス以外の様々な処理に用いることもできる。そこで、粒子法シミュレーションに係る計算もＧＰＵで高速化することも試みられている。例えば本発明者は、ＳＰＨのすべてのシミュレーションをＧＰＵを用いて行う手法を開発した（例えば非特許文献２，３を参照）。
田中正幸等、「粒子法を用いた剛体計算手法の開発とコンピュータ・グラフィックスへの適用」、日本機械学会、第１９回計算力学講演会講演論文集、７０１−７０２頁、２００６ T. Harada, S. Koshizuka, and Y. Kawaguchi Real-time fluid simulation coupled with cloth. In Proc. of Theory and Practice of Computer Graphics, 2007. T. Harada, S. Koshizuka, and Y. Kawaguchi Smoothed particle hydrodynamics on gpus. In Proc. of Computer Graphics International, pages 63-70, 2007.

本発明者が提案した上記の手法では、粒子法シミュレーションにおける近傍粒子の探索に関して空間上に固定された格子を用いたため、メモリ効率が悪く、計算領域について制限が生じるという問題があった。

例えば、粒子を用いる流体計算手法は、格子を用いる流体計算手法に比して、一般的に計算格子を用いないため、細かい飛沫等も容易に計算でき、実際の流体挙動を表現できるという利点がある。他方、粒子を用いる流体計算手法は、計算粒子同士が接続関係を持たず、動的に粒子の配置を変更するため、或る座標での物理量を計算するためには各タイムステップで、当該座標の近傍に存在する粒子を探索する処理が必要となる。

近傍粒子探索を効率的に行うためには、計算領域に格子（計算格子）を配置し、或る粒子の番号を、当該粒子が入っているボクセル（またはバケットもしくはグリッド）に格納することで効率化を図ることができる。一般的には計算領域を内包するように計算格子が用意される。

上記において、計算格子を固定化して用意すると、流体の計算領域内の粒子分布が変化するときには、粒子が存在しないボクセルが多くなり、メモリの使用効率が悪くなる。シミュレーションプログラムが使用できるメモリの大きさは限りがあるため、計算格子に用いるメモリを少なくすればするほど、流体を表現する粒子に用いることができるメモリが増え、より大規模な流体シミュレーションを行うことができる。

以上の問題点から、粒子法シミュレーションの近傍粒子探索に用いる計算格子のメモリ効率を高めることが望まれている。

本発明の目的は、上記の課題に鑑み、粒子を用いた流体計算手法に基づく粒子法シミュレーションで、近傍粒子探索に用いる計算格子の用い方について、メモリの使用効率を良好にし、かつ近傍粒子探索を効率化することができる粒子法シミュレーションのためのスライスデータ構造を提供することにある。

また本発明の他の目的は、メモリの使用効率を良好にすることができる上記のスライスデータ構造を容易にＧＰＵに実装化することができるスライスデータ構造を利用した粒子法シミュレーションのＧＰＵへの実装方法を提供することにある。

本発明は、上記の目的を達成するために、次のように構成される。

請求項１に係る本発明の粒子法シミュレーションのためのスライスデータ構造は、ＣＰＵやＧＰＵ等を用いて粒子法シミュレーションの計算で用いられる空間についてのデータ構造であって、空間は多数のボクセルで構成される３次元計算空間であり、３次元計算空間で任意の１つの軸に直交する複数の２次元スライス領域（スライス）を形成することにより、多数のボクセル（格子）を複数の２次元スライス領域で分割し、複数の２次元スライス領域の各々で、粒子が存するボクセルの存在範囲について最大と最小のボクセルの各々の開始座標を計算することにより、存在範囲を矩形の形状で囲む矩形範囲（バウンディングボックス）を決定し、複数の２次元スライス領域の各々の矩形範囲内に含まれるボクセルについてメモリを確保することを特徴とする。

請求項２に係る本発明の粒子法シミュレーションのためのスライスデータ構造は、上記の構成において、粒子が存在するボクセルについて確保されたメモリの格納要素には当該粒子を特定する値（例えば粒子番号等）が格納されることを特徴とする。

請求項３に係る本発明の粒子法シミュレーションのためのスライスデータ構造は、上記の構成において、複数の２次元スライス領域の各々の矩形範囲内における或る点（ｘ，ｙ，ｚ）が存在するボクセルのボクセル番号ｖ（ｘ，ｙ，ｚ）の計算は、各々の２次元スライス領域の矩形範囲を定義する値、スライス番号を決定する値、先頭番号を用いて行われることを特徴とする。

請求項４に係る本発明の粒子法シミュレーションのためのスライスデータ構造は、上記の構成において、粒子法シミュレーションの計算は近傍粒子探索の計算であり、近傍粒子探索の計算を効率化するために３次元計算空間が２次元スライス領域により分割されることを特徴とする。

請求項５に係る本発明の粒子法シミュレーションのためのスライスデータ構造は、上記の構成において、粒子法シミュレーションの計算について３次元計算空間を２次元スライス領域により分割する計算はＧＰＵを用いて行われることを特徴とする。

請求項６に係る本発明のスライスデータ構造を利用した粒子法シミュレーションのＧＰＵへの実装方法は、粒子法シミュレーションの近傍粒子探索の計算を効率化するためのデータ構造を構築し、そのデータ構造を用いて近傍粒子探索をＧＰＵ上で行うため必要なデータをテクスチャとしてビデオメモリに格納してＧＰＵへ実装する方法であって、多数のボクセルで構成される３次元計算空間で任意の１つの軸に直交する複数の２次元スライス領域を形成して多数のボクセルを複数の２次元スライス領域で分割し、複数の２次元スライス領域の各々で、粒子が存するボクセルの存在範囲について最大と最小のボクセルの各々の開始座標を計算することにより、存在範囲を矩形の形状で囲む矩形範囲を決定し、さらに、少なくとも１次元のテクスチャを用意し、複数の２次元スライス領域の各々を定義するに最低限必要な４つの値をテクスチャのＲＧＢＡチャンネルに格納し、２次元テクスチャを用意し、粒子の粒子番号を前記２次元テクスチャに格納することを特徴とする。

請求項７に係る本発明のスライスデータ構造を利用した粒子法シミュレーションのＧＰＵへの実装方法は、最低限の４つの値は、任意の２次元スライス領域（ｉ）における先頭のボクセルの番号ｐ_ｉ、最小のボクセルの開始座標ｂｘ_{ｉ，ｍｉｎ}とｂｚ_{ｉ，ｍｉｎ}、所定の軸方向のボクセル数ｎｘ_ｉであることを特徴とする。

本発明によれば、固定格子を用いず、３次元の計算空間を１つの軸方向に直交する２次元スライス領域に分割して、所定の矩形範囲すなわちバウンディングボックスと考え方を利用してメモリに格納するボクセルを制限したため、ＧＰＵ等のビデオメモリを効率よく利用することができ、ＧＰＵ等を用いた粒子シミュレーションの近傍粒子探索に用いるボクセル（格子）のメモリ効率を良好にすることができ、近傍粒子探索を効率化することができる。

以下に、本発明の好適な実施形態（実施例）を添付図面に基づいて説明する。

図１は、ＧＰＵを用いた粒子法シミュレーションで、近傍粒子探索の計算でスライスデータ構造を用い、このスライスデータ構造を利用して構築されるシミュレーション方法が実施されるコンピュータシステムの構成を示す。

コンピュータシステム１０は、例えばＰＣ（Personal Computer）を利用して構成され、演算処理部としての通常のコンピュータの演算処理を行うＣＰＵ１１と、画像処理に特化された演算処理を行うＧＰＵ１２を備えている。またコンピュータシステム１０は、記憶装置としてメインメモリ１３Ａとビデオメモリ（ＶＲＡＭ）１３Ｂを備え、周辺装置として少なくとも入力装置１４および表示装置１５を備えている。ＣＰＵ１１とＧＰＵ１２とメインメモリ１３Ａと入力装置１４と表示装置１５はバス１６を介して相互に接続されている。入力装置１４は入力インタフェース１７を介してバス１６に接続され、表示装置１５は出力インタフェース１８を介してバス１６に接続されている。またビデオメモリ１３Ｂは、ＧＰＵ１２に直接に接続されている。ビデオメモリ１３Ｂに対してはＧＰＵ１２を通してデータの行き来が行われる。

ＣＰＵ（中央演算処理装置）１１は、例えばＰＣ等に装備される通常的な演算処理部である。またＧＰＵ１２は「Graphics Processing Unit（グラフィックス演算プロセッサ）」であり、画像演算処理を実行する。ＣＰＵには例えば「Core X6800」が使用される。

ＧＰＵ１２は、複数のプロセッサを内蔵して成り、並列計算機として用いられる。ＧＰＵ１２としては、例えばＮＶＩＤＩＡ社製の「GeForce 8800GTX」が使用される。なお、ＧＰＵ１２で使用される製品はこれに限定されず、類似した構造・機能を有する他のストリーミングプロセッサを用いることができる。このＧＰＵ１２は、メインメモリ１３Ａに格納された近傍粒子探索プログラム１３Ａ−１を実行し、粒子法シミュレーションにおける近傍粒子探索演算を実施する。近傍粒子探索プログラム１３Ａ−１は、近傍粒子探索に用いるデータ構造の構築を実行する。本発明に係る特徴的な構造を有するデータ構造の構築方法については後で詳細に説明する。

またＧＰＵ１２は、メインメモリ１３Ａに格納された通常的な画像処理プログラム１３Ａ−２を実行し、近傍粒子探索の計算で得られかつテクスチャに格納された画像データを用いて画像処理を実施し、作成した画像データを用いて表示装置１５の画面に描画（表示のためのレンダリング）する。表示装置１５の画面には、粒子法シミュレーションの計算に基づいて、ＣＧ画像が表示される。

本実施形態の以下の説明では、例えば、粒子法シミュレーションに基づく流体シミュレーションの計算例が説明される。この流体シミュレーションの演算では、主にＧＰＵ１２の特徴的な内部構造および機能を利用して、近傍粒子探索に用いるためのデータ構造の構築手法が実行される。本実施形態で説明されるデータ構造を、後述するごとく、「スライスデータ構造（またはスライスベースデータ構造」と呼ぶこととする。なおＣＰＵを利用して実行することも可能である。

次に、図２を参照して上記のＧＰＵ１２内の要部の基本構造を説明する。図２はＧＰＵ１２の内部構造を模式的に示すブロック図である。ＧＰＵ２１は、並列的に配置された複数のプロセッサ２１を備え、これらのプロセッサによる並列計算処理に基づいて、表示（描画）のためのレンダリングが実行される。特に本実施形態に係る近傍粒子探索に用いるデータ構造の構築（ボクセル（またはバケット）の構築）の方法によれば、１つのボクセルにに４個の粒子が存在するとき、４個の粒子の粒子番号を１つのピクセルのＲＧＢＡチャンネルの各々に書き込む。ＧＰＵ１２における当該ボクセルの構築の処理によって画像演算処理の高速化が実現される。

ＧＰＵ１２は、通常、前段側で並列的に配置された例えば８個のバーテックス（頂点）シェーダ（ＶＳ）２２と、後段側で並列的に配置された例えば２４個のフラグメントシェーダ（ＦＳ）２３とを有している。複数のバーテックス（頂点）プロセッサ（２１Ａ）の各々で実行されるプログラムがバーテックスシェーダ２２であり、フラグメントプロセッサ（２１Ｂ）の各々で実行されるのがフラグメントシェーダ２３である。

上記において、バーテックスプロセッサは、並列的に配置されたプロセッサであり、多数の頂点の座標変換（入力された頂点座標を、レンダリングされる画像の空間での座標に変換すること）を一斉に処理することに特化されたプロセッサである。ＧＰＵ２１は、頂点の座標変換を、ＣＰＵに比較して高速に計算することができる。バーテックスプロセッサの動作を指定するものがバーテックスシェーダである。

さらにＧＰＵ１２には、深度テスト機能部２４、カラーマスク機能部２５、ステンシルテスト機能部２６が備えられている。またブロック２７はシェーダ指令分配部である。

ＧＰＵ１２を用いた粒子法シミュレーションで使用される本実施形態に係るデータ構造としての「スライスデータ構造」の基本的な概念について説明する。

粒子法に基づく流体計算では、各タイムステップで実行される近傍粒子探索のための計算に用いる空間（または領域）を設定する。この空間を「計算空間」と呼ぶ。この計算空間は、図３に示すごとく、互いに直交するＸ軸、Ｙ軸、Ｚ軸の３軸で定義される３次元の空間３１である。計算空間３１においては多数の格子３２が用意される。多数の格子３２は計算空間３１を形成し、各格子３２は一辺が「ｄ」の小さい立方体である。本実施形態の説明では、当該格子３２を「ボクセル」と呼ぶことにする。なお、「ボクセル」については、他に「バケット」または「グリッド」という呼び方もある。

上記の多数の格子すなわちボクセル３２は無限の広がりを持っており、計算中には計算空間３１上でその配置は変更されない。本実施形態による計算方法では、計算空間３１を形成する上記の多数のボクセル３２のうち、後述のごとく、必要なものだけビデオメモリ１３Ｂ上に保持する。

なお従来の計算手法のごとく固定の格子を用いる場合には、計算空間内において計算領域のためのバウンディングボックスを定め、当該バウンディングボックスの内部のボクセルのために粒子の有無に拘わらずメモリを確保するようにしていた。つまり、粒子が存在しないボクセルについてまでメモリを確保するようにしていた。これに対して、本実施形態の計算手法では、原則的に、粒子が存在しているボクセルにメモリを確保するようにしている。

本実施形態の計算手法では、まず、計算空間３１に関して１つの軸に直交する複数のスライス３１ａを形成し、計算空間３１に用意された多数のボクセル３２を複数のスライス３１ａのそれぞれに含まれるように分割する。各スライス３１ａの次元は、計算空間３１の次元（３次元）よりも１次元低いものとなり、２次元である。各スライス３１ａの厚みに相当する寸法はボクセル３２の一辺（ｄ）にほぼ等しい。すなわち、スライス３１ａが直交する１つの軸の方向には１ボクセル分だけを有することになる。この「１つの軸」としては、例えばＸ軸、Ｙ軸、Ｚ軸などの任意の１つの軸を選択することができる。例えば「１つの軸」としてＹ軸を選択すると、スライス３１ａはＸＺ軸で決まるＸＺ平面方向に広がる２次元空間となる。以下では、３次元の計算空間３１をＹ軸に垂直な２次元のスライス３１ａに分割した場合の例を説明する。

計算空間３１を複数のスライス３１ａに分割した後に、それぞれのスライス３１ａでＸＺ軸方向に関してのバウンディングボックスを定義する。当該バウンディングボックスは、１つのスライス３１ａ内において粒子が存在する領域を矩形の領域によって可能限り狭く確定するための境界で区切られる範囲を意味している。このスライス３１ａにおいてバウンディングボックスを定義することは、ＸＺ軸方向における最大と最小のボクセルの開始座標ｂｘ_{ｉ，ｍａｘ}、ｂｘ_{ｉ，ｍｉｎ}、ｂｚ_{ｉ，ｍａｘ}、ｂｚ_{ｉ，ｍｉｎ}を計算することを意味する。ここで「ｉ」はスライスを特定する番号を意味する。ＸＺ軸方向における最大と最小のボクセルの開始座標ｂｘ_{ｉ，ｍａｘ}、ｂｘ_{ｉ，ｍｉｎ}、ｂｚ_{ｉ，ｍａｘ}、ｂｚ_{ｉ，ｍｉｎ}を計算することにより、当該１つのスライス３１ａのバウンディングボックスに含まれるボクセルの数、換言するとスライス３１ａに含まれるＸＺ軸方向のボクセルの数ｎｘ_ｉ，ｎｚ_ｉは下記（数１）の（１）式と（２）式で求められる。

上記のｎｘ_ｉ，ｎｚ_ｉを用いることにより番号ｉのスライス３１ａでのボクセルの数ｎ_ｉは、ｎ_ｉ＝ｎｘ_ｉ×ｎｚ_ｉと求めることができる。本実施形態の計算手法では、各スライス内のバウンディングボックスの内部の領域のみのボクセルに関してメモリを確保する。

ここで、計算空間３１はｎ個のスライス｛Ｓ_０，Ｓ_１，Ｓ_２，…，Ｓ_ｎ−１｝に分解されるとする。最初のスライスＳｉでの先頭のボクセルの番号ｐ_ｉはスライスＳ_０からＳ_ｉ−１にあるボクセルの数の和として下記の（数２）の（３）式を計算する。

図４に、従来の計算手法による固定格子の構造を用いる場合のバウンディングボックスの例（Ａ）と、本実施形態の計算手法によるボクセルを決めるためのバウンディングボックスの例（Ｂ）とを示す。図４の（Ａ）と（Ｂ）では、図示されるごとくＸ軸、Ｙ軸、Ｚ軸が定義されている。図４の（Ａ）と（Ｂ）は、Ｚ軸方向に対して直交する１つのＸＹ軸方向の１つの正方形の面（ＸＹ平面）で示している。図４の（Ａ）と（Ｂ）において、正方形４１で定められた範囲内の小さい正方形４１ａは１つのボクセルを示し、さらに円４２はボクセル内に存在している粒子を示している。従って、正方形４１の範囲内には粒子４２が存在するボクセルと粒子４２が存在しないボクセルとが含まれている。

図４の（Ａ）と（Ｂ）でＹ軸方向に直交する番号４３で示される層状の部分が前述したスライス３１ａに対応している。

なお、図４において上記の番号４１ａ，４２，４３はそれぞれ一部のみを指しており、その他の同じものには番号を付するのを省略している。

従来の計算手法によれば、固定格子の構造であるため、図４の（Ａ）に示されるごとく粒子の有無に関係なく正方形の範囲４１（太い実線で示す）に含まれるすべてのボクセルについてメモリが確保されるように構築されるものであった。すなわち、粒子４２の存在しない多数のボクセルについてまでメモリを確保する必要があった。これに対して、本実施形態の計算手法では、各スライスで定義されたバウンディングボックスが合成されて作られた領域４４（太い実線で示す）、すなわち可能な限り粒子４２が存在するボクセルだけを含むように狭く作られた領域４４に含まれたボクセルについてメモリを確保するように構築される。

計算空間３１に関して上記スライスを利用して前述のごとくバウンディングボックスを形成することによって、或る点（ｘ，ｙ，ｚ）が存在するボクセルの番号は次のように計算される。

まずＹ軸方向の最小のボクセルの開始座標ｂｙ_ｍｉｎを用いて、その点が含まれるスライスの番号ｉを下記の（数３）に示す（４）式で計算する。

そして点（ｘ，ｙ，ｚ）のボクセル番号はｖ（ｘ，ｙ，ｚ）は、下記の（数４）で示す（５）式で計算される。

こうして上記の点（ｘ，ｙ，ｚ）が存在するボクセルのボクセル番号が計算される。式（５）を参照すると、或る点が存在するボクセルの番号を計算するために必要な値は、各スライスのバウンディングボックスを定義する値ｂｘ_{ｉ，ｍｉｎ}、ｂｚ_{ｉ，ｍｉｎ}、ｎｘ_ｉ、スライス番号を決定する値ｂｙ_ｍｉｎ、各スライスの先頭番号ｐ_ｉ、そしてボクセルの一辺の長さｄであることが分かる。スライスを利用した上記のデータ構造を用いると、図４の（Ｂ）に示されるように、メモリを効率よく用いてボクセル（格子）を作ることができる。

以上によって、本実施形態による計算手法で特徴的なデータ構造である「スライスデータ構造」が作られる。上記のスライスデータ構造を実際に構築するためには、ｎ個のスライス｛Ｓ_０，Ｓ_１，Ｓ_２，…，Ｓ_ｎ−１｝の各々のバウンディングボックスを決定し、決定されたバウンディングボックスから各スライスの先頭番号を計算する。そして、ボクセル番号ｖ（ｘ，ｙ，ｚ）を計算するための上記４つの値（ｐ_ｉ、ｂｘ_{ｉ，ｍｉｎ}、ｂｚ_{ｉ，ｍｉｎ}、ｎｘ_ｉ）を求めた後、これらの値をボクセルに書き込んでいく。これらの４つの値はバウンディングボックスを定義するための最低限の値であり、これらの値のみに限定されるものではない。

次に、上記スライスデータ構造に基づく計算手法がＧＰＵ１２を用いて実行できるようにＧＰＵ１２上に実装する仕組みについて説明する。

ＧＰＵ１２を用いた本実施形態の計算手法の実装について図５に示す。図５に示されるごとく、ビデオメモリ１３Ｂでのデータ構造の設定の段階（Ｓ１１）、計算空間３１を分割する多数のスライス３１ａの各々でのバウンディングボックスの決定の段階（Ｓ１２）、各スライス３１ａでの粒子の先頭番号の計算の段階（Ｓ１３）、値の格納の段階（Ｓ１４）、および値の読み出しの段階（Ｓ１５）の５つの段階によって構成される。

なお上記の「スライスデータ構造」については、「スライスデータ構造の構築」という観点、および「スライスデータ構造の利用」という観点で把握することができる。「スライスデータ構造の構築」という観点によれば、「スライスデータ構造の構築」は、上記の５つの段階Ｓ１１〜Ｓ１５のうち、「データ構造の設定」、「バウンディングボックスの決定」、「先頭番号の計算」、および「値の格納」によって構成される。また「スライスデータ構造の利用」という観点によれば、上記の５つの段階Ｓ１１〜Ｓ１５のうち、「値の読み出し」により構成される。

以下、５つの段階Ｓ１１〜Ｓ１５での処理内容を分説する。

データ構造の設定の段階Ｓ１１：
最初に、図６の処理フローを参照して、ＧＰＵ１２上での実装を可能にするためのデータ構造を説明する。ＧＰＵ１２上で計算を行うためには、計算空間におけるスライス、ボクセル、および粒子に関するデータをテクスチャ１３Ｂ−１としてビデオメモリ１３Ｂに格納することが必要である。

前述の通り計算空間３１は複数のスライス３１ａに分割されるので、各スライス３１ａに必要とされる値、すなわち前述の４つの値ｐ_ｉ、ｂｘ_{ｉ，ｍｉｎ}、ｂｚ_{ｉ，ｍｉｎ}、ｎｘ_ｉを格納するために１次元のテクスチャを用意する（ステップＳ２１）。複数のスライス３１ａの数が１次元テクスチャの大きさの最大値を超える場合には、２次元のテクスチャとすることもできる。各スライス３１ａに必要な４つの値ｐ_ｉ、ｂｘ_{ｉ，ｍｉｎ}、ｂｚ_{ｉ，ｍｉｎ}、ｎｘ_ｉ は１次元テクスチャの１つのＲＧＢＡチャンネルに格納される（ステップＳ２２）。ＲＧＢＡチャンネルを利用して格納するため、通常、１枚のテクスチャで十分に対応することができる。

また、格子すなわちボクセルに保持する値を格納するために２次元テクスチャを用意する（ステップＳ２３）。本実施形態の例において、ボクセルに保持する値は粒子番号である。２次元テクスチャを利用してボクセルに保持する粒子番号を格納する（ステップＳ２４）。粒子番号を格納する２次元テクスチャを「インデックスプールテクスチャ」と呼ぶこととする。本実施形態の例では、１個のボクセルに１テクセルを割り当てるようにしている。この場合において、２次元テクスチャのテクセル数が、インデックスプールテクスチャに格納できるボクセルの最大数となる。

バウンディングボックスの決定の段階Ｓ１２：
次に、図７の処理フローを参照してバウンディングボックスの決定の手順を説明する。スライス３１ａでのバウンディングボックスを決定するためには、まず各スライスでの最大と最小のボクセルの開始座標を計算する（ステップＳ３１）。この計算の処理では、粒子のＹ座標値から当該粒子が属するスライスの番号を計算し、さらに各スライスに存在する粒子の座標の最大値と最小値を計算することにより行われる。ＧＰＵ１２を用いているので、当該ＧＰＵ１２では、粒子１個につき１個の頂点を割り当て、用意した前述の１次元テクスチャに書き込む（ステップＳ３２）。

ＧＰＵ１２のバーテックスシェーダでは、粒子のＹ座標値から上記の式（４）を用いてスライス番号（ｉ）を計算し、すなわち、上記した１次元テクスチャに書き込む位置を計算し（ステップＳ３３）、計算した当該位置に頂点を書き込む（ステップＳ３４）。そしてＧＰＵ１２のフラグメントシェーダでは、Ｘ座標値およびＺ座標値を「色」として出力する（ステップＳ３５）。なお、最大値と最小値の選択は、ＧＰＵ１２のアルファブレンディングの機能を用いて行うことが可能である。

先頭番号の計算の段階Ｓ１３：
通常、複数のスライス３１ａの各々でのＸＺ座標の最大値と最小値を用いることでそのスライス上のボクセルの数を計算することが可能である。先頭番号の計算は、各スライスにおいて式（３）を計算しなければならない。この処理を、スライスｉについてｎ_０からｎ_ｉ−１までの値を読み出して計算することになる。スライスの総数がｍであるとすると、従来の計算手法によれば、（ｍ（ｍ＋１））／２回のメモリアクセスが必要になり、このためメモリ転送速度により処理の速度が制限されることになる。

そこで、本実施形態の計算手法では、ＧＰＵ１２のバーテックスシェーダを利用した機能で実行される。

すなわち、スライスｉに存在するボクセルの数ｎｉを読み出して（ｉ＋１）からｍ番目のスライスの先頭番号配列に加算していく。これによれば、ｍ回のメモリアクセスで先頭番号を計算することができる。この処理はいわゆる「Scattering Operation」であり、ＧＰＵ１２を用いる場合にはバーテックスシェーダを用いて行うことができる。この処理は、スライスがｍ個あるとすると、要素数ｍの１次元配列に書き込むという処理である。つまり、ＧＰＵ１２を用いる場合には、ｍ個のテクセルから成る１次元テクスチャに書き込む処理になる。スライスｉに存在するボクセル数ｎｉを（ｉ＋１）番目からｍ番目のスライスの先頭番号配列に書き込む処理は、（ｉ＋１）番目のテクセルからｍ番目のテクセルまでの線分を描画することで行うことができる。

上記の処理は、ＧＰＵ１２の機能によれば、フラグメントシェーダの機能によりスライス内のボクセルの数を色として出力する。さらにＧＰＵ１２のアルファブレンディングの機能により加算処理を行うことにより式（３）を評価し、各スライスの先頭番号を計算することが可能となる。

図８は、上記のスライスの「先頭番号の計算」の考え方を図解したものである。図８の（Ａ）は、前述した図４の（Ｂ）の内容と同じである。図８の（Ｂ）は、（Ａ）で示した下から上まで例えば１０層のスライスの各々におけるバウンディングボックス内に含まれるボクセルの数を対応する箇所に数値で示している。図８の（Ｃ）は、各スライスのボクセルの数が加算される状態をイメージ的に図解したものである。図８の（Ｄ）は下側から上側に向かってスライスの総数に関する累積値を示すものである。

値の格納の段階Ｓ１４：
上記のごとくして計算したスライスの先頭番号とバウンディングボックスを定義する値を用いて粒子番号をインデックスプールテクスチャに格納していく。

値の読み出しの段階Ｓ１５：
図９に示した処理フローを参照する。粒子のＹ座標値からスライス番号を計算し（ステップＳ４１）、当該スライスの先頭番号ｐ_ｉとスライスでの最小値ｂｘ_{ｉ，ｍｉｎ}、ｂｚ_{ｉ，ｍｉｎ}、スライスでのＸ軸方向のボクセル数ｎｘ_ｉを読み出す（ステップＳ４２）。これらの値を用いて上記の式（５）を計算し、その粒子に割り当てられているボクセルのインデックスを計算する（ステップＳ４３）。計算されたインデックスに対応するテクセルへの値の書き込みは、粒子１個に頂点を１個を割り当て、そのテクセルの粒子番号を読み出すことで行われる（ステップＳ４４）。

或る座標が存在するボクセルに格納されている値の読み出しは、その座標から値の格納時と同様にボクセルのインデックス、すなわちメモリのアドレスを計算することによって行われる。

スライスデータ構造を利用した本実施形態に係る計算手法は、粒子番号を格納する格子データ構造の構築と、粒子番号の格納までの処理をＧＰＵ１２で行うことができ、それにより粒子法シミュレーションのすべての工程をＧＰＵ１２上で行うことができる。

本実施形態に係る計算手法は、固定格子を用いず、計算領域の制限を受けないので、広い計算空間であっても容易に計算することができ、粒子が存在しないボクセルを排除することによりメモリを効率よく使用することができる。さらに、粒子の個数を大きくしても通常的な使用されるメモリを用いて粒子法シミュレーションを実行することができる。

本発明は、ＧＰＵ等を用いた粒子シミュレーションの計算でビデオメモリの効率を良好にするのに利用される。

本発明に係る粒子法シミュレーションが実行されるコンピュータシステムの基本的な構成を示すシステム構成図である。 コンピュータシステムに実装されたＧＰＵの内部の基本的構成を示すブロック図である。 ３次元の計算空間を図示する斜視図である。 従来の固定格子（Ａ）と本実施形態に係る格子（ボクセル）の決め方（Ｂ）とを対比して説明するための図である。 スライスデータ構造による計算手法をＧＰＵに実装するための手順を示すフローチャートである。 データ構造の設定の段階についてのより具体的な手順を示すフローチャートである。 バウンディングボックスの決定の段階についてのより具体的な手順を示すフローチャートである。 スライスの「先頭番号の計算」の考え方を図解したものである。 値の読み出しの段階についてのより具体的な手順を示すフローチャートである。

符号の説明

１０ コンピュータシステム
１１ ＣＰＵ
１２ ＧＰＵ
１３Ａ メインメモリ
１３Ａ−１ 近傍粒子探索プログラム
１３Ａ−２ 画像処理プログラム
１３Ｂ ビデオメモリ
１３Ｂ−１ テクスチャ

Claims (7)

1. 粒子法シミュレーションの計算で用いられる空間についてのデータ構造であって、
   前記空間は多数のボクセルで構成される３次元計算空間であり、
   前記３次元計算空間で任意の１つの軸に直交する複数の２次元スライス領域を形成することにより、前記多数のボクセルを前記複数の２次元スライス領域で分割し、
   前記複数の２次元スライス領域の各々で、粒子が存するボクセルの存在範囲について最大と最小のボクセルの各々の開始座標を計算することにより、前記存在範囲を矩形の形状で囲む矩形範囲を決定し、
   前記複数の２次元スライス領域の各々の前記矩形範囲内に含まれる前記ボクセルについてメモリを確保することを特徴とする粒子法シミュレーションのためのスライスデータ構造。
2. 前記粒子が存在する前記ボクセルについて確保された前記メモリの格納要素には前記粒子を特定する値が格納されることを特徴とする請求項１記載の粒子法シミュレーションのためのスライスデータ構造。
3. 前記複数の２次元スライス領域の各々の前記矩形範囲内における或る点（ｘ，ｙ，ｚ）が存在する前記ボクセルのボクセル番号ｖ（ｘ，ｙ，ｚ）の計算は、各々の前記２次元スライス領域の前記矩形範囲を定義する値、スライス番号を決定する値、先頭番号を用いて行われることを特徴とする請求項１記載の粒子法シミュレーションのためのスライスデータ構造。
4. 前記粒子法シミュレーションの前記計算は近傍粒子探索の計算であり、前記近傍粒子探索の計算を効率化するために前記３次元計算空間が前記２次元スライス領域により分割されることを特徴とする請求項１記載の粒子法シミュレーションのためのスライスデータ構造。
5. 前記粒子法シミュレーションの前記計算について前記３次元計算空間を前記２次元スライス領域により分割する計算はＧＰＵを用いて行われることを特徴とする請求項１〜４のいずれか１項に記載の粒子法シミュレーションのためのスライスデータ構造。
6. 粒子法シミュレーションの近傍粒子探索の計算を効率化するためのデータ構造を構築し、そのデータ構造を用いて近傍粒子探索をＧＰＵ上で行うため必要なデータをテクスチャとしてビデオメモリに格納してＧＰＵへ実装する方法であって、
   多数のボクセルで構成される３次元計算空間で任意の１つの軸に直交する複数の２次元スライス領域を形成して前記多数のボクセルを前記複数の２次元スライス領域で分割し、
   前記複数の２次元スライス領域の各々で、粒子が存するボクセルの存在範囲について最大と最小のボクセルの各々の開始座標を計算することにより、前記存在範囲を矩形の形状で囲む矩形範囲を決定し、
   少なくとも１次元のテクスチャを用意し、
   前記複数の２次元スライス領域の各々を定義するに最低限必要な４つの値を前記テクスチャのＲＧＢＡチャンネルに格納し、
   ２次元テクスチャを用意し、
   前記粒子の粒子番号を前記２次元テクスチャに格納する、
   ことを特徴とするスライスデータ構造を利用した粒子法シミュレーションのＧＰＵへの実装方法。
7. 最低限の前記４つの値は、任意の２次元スライス領域（ｉ）における先頭のボクセルの番号ｐ_ｉ、最小のボクセルの開始座標ｂｘ_{ｉ，ｍｉｎ}とｂｚ_{ｉ，ｍｉｎ}、所定の軸方向のボクセル数ｎｘ_ｉであることを特徴とする請求項３記載のスライスデータ構造を利用した粒子法シミュレーションのＧＰＵへの実装方法。

Images