JP2010501100A - 準モンテカルロ法を使用するマルコフ連鎖の同時シミュレーション - Google Patents

準モンテカルロ法を使用するマルコフ連鎖の同時シミュレーション Download PDF

Info

Publication number
JP2010501100A
JP2010501100A JP2009524776A JP2009524776A JP2010501100A JP 2010501100 A JP2010501100 A JP 2010501100A JP 2009524776 A JP2009524776 A JP 2009524776A JP 2009524776 A JP2009524776 A JP 2009524776A JP 2010501100 A JP2010501100 A JP 2010501100A
Authority
JP
Japan
Prior art keywords
graphics system
computer graphics
simulating
ray
sorting
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2009524776A
Other languages
English (en)
Other versions
JP4947394B2 (ja
Inventor
アレクサンドリア ケラー，
カーステン ヴェヒター，
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nvidia ARC GmbH
Original Assignee
Mental Images GmbH
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Mental Images GmbH filed Critical Mental Images GmbH
Publication of JP2010501100A publication Critical patent/JP2010501100A/ja
Application granted granted Critical
Publication of JP4947394B2 publication Critical patent/JP4947394B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/00—3D [Three Dimensional] image rendering
    • G06T15/06—Ray-tracing

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Graphics (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Generation (AREA)
  • Image Analysis (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

コンピュータグラフィックスシステムがマルコフ連鎖（したがって光子の軌道など）をより効率的にシミュレートすることを可能にするように動作可能な、方法、システム、装置、及びコンピュータソフトウェア／コンピュータコード製品は、準モンテカルロ法を使用してマルコフ連鎖をシミュレートするステップ、及び／又はシミュレートするための手段を含み、マルコフ連鎖をシミュレートするステップは、状態をソートするステップを含み、ソートするステップは、近接性ソーティングを含む。
【選択図】 図９

Description

関連出願の相互参照
本特許出願（整理番号ＭＥＮＴ−１０１−Ｂ）は、２００６年８月１５日に出願された米国仮特許出願第６０／８２２，４１７号の優先権の利益を主張する。本出願は、２００６年６月２３日に出願された米国特許出願第１１／４７４，５１７号（整理番号ＭＥＮＴ−１０１−ＵＳ）の一部継続出願でもあり、同出願は、２００５年６月２３日に出願された米国仮特許出願第６０／６９３，２３１号（整理番号ＭＥＮＴ−１０１−ＰＲ）の利益を主張し、また２００２年１１月１９日に出願された米国特許出願第１０／２９９，５７４号（整理番号ＭＥＮＴ−０７５）の一部継続出願でもある。
米国特許出願第１０／２９９，５７４号（整理番号ＭＥＮＴ−０７５）は、２００１年６月１９日に出願された米国特許出願第０９／８８４，８６１号（整理番号ＭＥＮＴ−０６１）の一部継続出願であり、同出願は、２００１年２月１日に出願された米国仮特許出願第６０／２６５，９３４号、及び２０００年６月１９日に出願された米国仮特許出願第６０／２１２，２８６号の優先権の利益を主張する。
整理番号ＭＥＮＴ−０６１、ＭＥＮＴ−０７５、ＭＥＮＴ−１０１−ＰＲ、及びＭＥＮＴ−１０１−ＵＳを含むが、それらに限定されない、上に列挙された特許出願、並びにそれらの仮特許出願の各々は、その全体が本明細書に記載されているかのように、その全体が参照により本明細書に組み込まれる。
コンピュータプログラムリストの参照
本出願には、以下にソースコードリストが添付される。そのソースコードリストは、本明細書では「コンピュータプログラムリスト」と呼ばれ、「１．１．１」の形式の３桁の参照番号によって識別されるセクションに構成される。
発明の分野
本発明は、一般に、動画及び他の応用例のためなどの、デジタルコンピューティングシステムにおける、またデジタルコンピューティングシステムによる、画像レンダリングのための方法及びシステムに関し、詳細には、リアルタイム高精度レイトレーシングのための方法、システム、デバイス、及びコンピュータソフトウェアに関する。
発明の背景
「レイトレーシング」という用語は、光源をカメラと結びつけるすべての光路を識別し、これらの寄与を総計することによって、フォトリアリスティックな画像を合成するための技法を指すものである。シミュレーションは、可視性を決定するために視線に沿って光線をトレースし、照度を決定するために光源から光線をトレースする。
レイトレーシングは、動画及び他の応用例における主流になっている。しかし、現在のレイトレーシング技法は、数値的問題、動的シーンを処理する能力の限界、アクセラレーションデータ構造（ａｃｃｅｌｅｒａｔｉｏｎ ｄａｔａ ｓｔｒｕｃｔｕｒｅ）の設定の遅さ、及び大きなメモリフットプリントを含む、多くの知られた限界及び弱点により影響をこうむる。したがって、現在のレイトレーシング技法は、風になびく森や人の髪の毛など、全体的に動きのあるシーンを効率的に処理するための能力を欠いている。現在のレイトレーシングシステムの限界の克服は、例えば、映画制作におけるより高品質なモーションブラー（ｍｏｔｉｏｎ ｂｌｕｒ）のレンダリングも可能にする。
レイトレーシングシステムの性能を改善するための現在の試みは、多くの理由で不十分である。例えば、現在のリアルタイムレイトレーシングシステムは一般に、そのアクセラレーション構造として、軸平行バイナリ空間分割（ａｘｉｓ−ａｌｉｇｎｅｄ ｂｉｎａｒｙ ｓｐａｃｅ ｐａｒｔｉｔｉｏｎ）に基づいた３Ｄツリーを使用する。これらのシステムの主な重点は、静的シーンをレンダリングすることにあるので、これらのシステムは一般に、全体的に動きのあるシーンに関して必要なデータ構造を構築するのに必要とされる、かなり長い設定時間に対処することができない。この線に沿って、ある製造業者は、効率的な３Ｄツリーを構築し、ツリーを横断するのに必要とされる時間を短縮できる技法を開発することによって、リアルタイムレイトレーシングを改良した。しかし、レイトレーシングされる物体の数が増加するにつれて、システムの予想メモリ要件が２次的に増大することが示され得る。
別の製造業者は、システム性能を改善するために境界ボリューム階層（ｂｏｕｎｄｉｎｇ ｖｏｌｕｍｅ ｈｉｅｒａｒｃｈｙ）を使用する、レイトレーシング集積回路を設計した。しかし、そのアーキテクチャの性能は、あまりに多くのインコヒーレントな２次光線がトレースされる場合、損なわれることが判明した。
加えて、フィールドプログラマブルゲートアレイ（ＦＰＧＡ）を使用して３Ｄツリー横断技法を実施することによって、システム性能を改善する試みも行われた。これらのシステムの処理速度の主な増加は、コヒーレントな光線の束をトレースし、高速ハードワイヤード計算を実行するＦＰＧＡの能力を活用することによって獲得される。
アクセラレーション構造の構築は、いまだハードウェアでは実施されていない。ＦＰＧＡ実装は一般に、精度を抑えた浮動小数点法（ｆｌｏａｔｉｎｇ ｐｏｉｎｔ ｔｅｃｈｎｉｑｕｅ）を使用する。
多くの異なる技法が、マルコフ連鎖（Ｍａｒｋｏｖ ｃｈａｉｎ）を含む光子軌道をシミュレートするために使用できる。マルコフ連鎖をシミュレートする際の準モンテカルロ法（ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ ｔｅｃｈｎｉｑｕｅ）の使用が、例えば、Ｌ’Ｅｃｕｙｅｒ他、Ｒａｎｄｏｍｉｚｅｄ Ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ Ｓｉｍｕｌａｔｉｏｎ ｏｆ Ｍａｒｋｏｖ Ｃｈａｉｎｓ ｗｉｔｈ ａｎ Ｏｒｄｅｒｅｄ Ｓｔａｔｅ Ｓｐａｃｅ（２００４）（「Ｌ’Ｅｃｕｙｅｒ」）において説明されており、同文献は、その全体があたかも説明されているかのように、参照により本明細書に組み込まれる。
Ｌ’Ｅｃｕｙｅｒは、全順序（離散又は連続）状態空間を有するマルコフ連鎖の各ステップにおける状態分布を推定するための、ランダム化準モンテカルロ法（ｒａｎｄｏｍｉｚｅｄ ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ ｍｅｔｈｏｄ）の使用について説明している。連鎖内のステップの数は、ランダムで無制限とすることができる。この方法は、特に、各ステップにおいて状態依存コストが支払われる場合の、何らかのランダムな停止時間までの予想トータルコストの、低分散不偏推定器（ｌｏｗ−ｖａｒｉａｎｃｅ ｕｎｂｉａｓｅｄ ｅｓｔｉｍａｔｏｒ）を取得するために使用することができる。当該論文は、数値的な説明を提供しており、その説明では、標準モンテカルロと比べて、分散低減が大きい。
Ｌ’Ｅｃｕｙｅｒは、全順序状態空間を有する離散時間及び離散状態マルコフ連鎖について、固定数のステップにわたって過渡的な量（ｔｒａｎｓｉｅｎｔ ｍｅａｓｕｒｅ）を推定するための決定論的準モンテカルロ法（ｄｅｔｅｒｍｉｎｉｓｔｉｃ ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ ｍｅｔｈｏｄ）が、以下の論文において提案され、研究されたことを述べており、以下の論文も、その全体が参照により本明細書に組み込まれる。
Ｌｅｃｏｔ他、「Ｑｕａｓｉ−Ｒａｎｄｏｍ Ｗａｌｋ Ｍｅｔｈｏｄｓ」、Ｍｏｎｔｅ Ｃａｒｌｏ ａｎｄ Ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ Ｍｅｔｈｏｄｓ ２０００、ｐｐ．６３〜８５、Ｓｐｒｉｎｇｅｒ−Ｖｅｒｌａｇ（Ｂｅｒｌｉｎ、２００２）
Ｌｅｃｏｔ他、「Ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ Ｍｅｔｈｏｄｓ ｆｏｒ Ｅｓｔｉｍａｔｉｎｇ Ｔｒａｎｓｉｅｎｔ Ｍｅａｓｕｒｅｓ ｏｆ Ｄｉｓｃｒｅｔｅ Ｔｉｍｅ Ｍａｒｋｏｖ Ｃｈａｉｎｓ」、Ｍｏｎｔｅ Ｃａｒｌｏ ａｎｄ Ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ Ｍｅｔｈｏｄｓ ２００２、ｐｐ．３２９〜４３、Ｓｐｒｉｎｇｅｒ−Ｖｅｒｌａｇ（Ｂｅｒｌｉｎ、２００４）（「Ｌｅｃｏｔ ２００４」）
Ｌ’Ｅｃｕｙｅｒは、基数２の（０，２）−列を使用して、連鎖のｎ＝２ｋ個のコピーを同数のステップにわたって並列にシミュレートする技法について説明している。連鎖のステップｊにおいて、当該技法は、ｎ個のコピーをそれらの状態に従って再度順序づけ、シミュレーションを行うために一様乱数の代わりに（０，２）−列のｎｊからｎｊ＋ｎ−１までの要素を利用することによって、ｎ個のコピーについて遷移（即ち次の状態）をシミュレートする。当該技法は、連鎖の各遷移のシミュレーションが単一の一様確率変数（ｕｎｉｆｏｒｍ ｒａｎｄｏｍ ｖａｒｉａｔｅ）を必要とすることを仮定する。Ｌｅｃｏｔ ２００４において、マルコフ連鎖の遷移行列の構造についてのある条件下で、正しい値への収束が証明された。
Ｌ’Ｅｃｕｙｅｒは、上記の技法を、ランダムで無制限な数のステップを有し、連続状態空間を有するマルコフ連鎖に一般化することを試みており（そのような技法は特に再生シミュレーションを扱うことを可能にすると主張されている）、それに対しては、マルコフ連鎖の１つのステップにおいて次の状態を生成するのに必要とされる一様確率変数の数ｄは、１より大きくなり得る。
しかし、光子軌道などをシミュレートするために、様々な研究者がマルコフ連鎖の使用を前提としているが、マルコフ連鎖、例えば光子軌道をシミュレートするために、準モンテカルロ法を使用するこれまでの技法は、ランダムサンプリングよりわずかに良い程度の改善を提供したにすぎなかった。しかし、以下に説明される本発明によるＱＭＣ法は、解演算子（ｓｏｌｕｔｉｏｎ ｏｐｅｒａｔｏｒ）の本質的な低次元構造を活用することができ、はるかに効率的な技法をもたらす。
したがって、本発明以前には、マルコフ連鎖を同時にシミュレートするための、より効率的な方法、システム、装置、デバイス、及びソフトウェアプログラム／コード製品に対する必要がまだ存在していた。
したがって、準モンテカルロ法及びソーティング戦略（ｓｏｒｔｉｎｇ ｓｔｒａｔｅｇｙ）を使用する、マルコフ連鎖をシミュレートするための、改良されたより効率的な方法、システム、装置、デバイス、及びソフトウェアプログラム／コード製品を提供するのが望ましい。
本発明は、画像内のピクセルのピクセル値を生成するためのコンピュータグラフィックスシステムにおける使用に適した、方法、システム、装置、及びコンピュータソフトウェア／コンピュータコード製品を提供し、ピクセル値は、シーン内の点を、シミュレートされるカメラの画像平面上に記録されるように表し、コンピュータグラフィックスシステムは、選択された方向に沿ってピクセルからシーンに向かって放たれた少なくとも１つの光線をシミュレートするステップを含む、選択されたレイトレーシング法を使用して、画像のピクセル値を生成するように構成され、レイトレーシング法は、光線とシーン内の物体の表面との交わりを計算するステップと、シーン内の物体を照明する光線の軌道をシミュレートするステップとをさらに含み、軌道をシミュレートするステップは、マルコフ連鎖を使用する。
本発明の一態様では、本発明による方法、システム、装置、及びコンピュータソフトウェア／コンピュータコード製品は、準モンテカルロ法を使用してマルコフ連鎖をシミュレートするステップ（及び／又はシミュレートするための手段）を含み、マルコフ連鎖をシミュレートするステップは、状態をソートするステップを含み、ソートするステップは、近接性ソーティング（ｐｒｏｘｉｍｉｔｙ ｓｏｒｔｉｎｇ）を含む。
本発明のさらなる態様は、複数のマルコフ連鎖を同時にシミュレートするステップを含み、複数のマルコフ連鎖を同時にシミュレートするステップは、Ｈａｌｔｏｎ列、Ｈａｌｔｏｎ列の変形、格子列、又は（ｔ，ｓ）−列のいずれかを使用して選択され得る、準モンテカルロ点を利用するステップを含む。複数のマルコフ連鎖を同時にシミュレートするステップは、任意の数の連鎖を利用することができる。
本発明のさらなる態様では、複数のマルコフ連鎖を同時にシミュレートするステップは、軌道分割（ｔｒａｊｅｃｔｏｒｙ ｓｐｌｉｔｔｉｎｇ）などによって、追加的な軌道をオンザフライで追加するステップを含む。シミュレートするステップは、粒子吸収をシミュレートする技法を利用するステップも含むことができ、ロシアンルーレット法（Ｒｕｓｓｉａｎ Ｒｏｕｌｅｔｔｅ ｔｅｃｈｎｉｑｕｅ）を使用するステップを含むことができる。
本発明の一態様では、ｓ次元問題について、複数のマルコフ連鎖を同時にシミュレートするステップは、以下の技法を含む。
ｎ：＝０
準モンテカルロ点ｘｌを使用して、Ｘ０，ｌを初期化する。
ループ：
− 適切な順序を使用して、状態ベクトルをソートする。
− ｎ：＝ｎ＋１
− 後続の準モンテカルロ点ｘｌを使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる。
本発明は、状態空間における近接性によって状態を順序づけるステップ、及びランダム化の技法も含むことができる。任意の選択されたランダム化を施された任意の選択された点集合又は点列が、本発明の方法と併せて利用できる。
別の態様では、本発明は、以下の技法をさらに含む。
ｎ：＝０、点列ｘｌを選択する。
ランダム化をインスタンス化する。
ｒａｎｄ（ｘｌ）を使用して、Ｘ０，ｌを初期化する。
ループ
・ 何らかの順序を使用して、状態ベクトルをソートする。
・ ｎ：＝ｎ＋１
・ ランダム化をインスタンス化する。
・ ｒａｎｄ（ｘｌ）を使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる。
ランダム化技法は、各時間ステップｎにおいて準モンテカルロ点をランダム化するステップを含むことができ、ランダム化するステップは、より詳細には、後続のサンプルがそこから取り出される、基数ｂ＝２の（ｔ，ｓ）−列を選択するステップと、サンプルにｓ次元ランダムベクトルによるＸＯＲ演算を施すステップであって、ランダムベクトルが、各遷移ステップの後に生成される、ステップを含む。
本発明の別の態様は、明るさ、又は状態のノルムによってソートするステップと、空間階層を使用してソートするステップと、少なくとも１つの空間充填曲線（ｓｐａｃｅ ｆｉｌｌｉｎｇ ｃｕｒｖｅ）によって列挙されたバケットを使用してソートするステップとを含む。空間階層を使用してソートするステップは、バイナリ階層を利用するステップを含むことができ、空間階層は、ＢＳＰツリー、ｋＤツリー、若しくはＢＳＰツリー構造の他の軸平行サブセットのいずれか、又は境界ボリューム階層、又は規則的ボクセル（ｖｏｘｅｌ）を含むことができる。バイナリ階層は、選択された発見的手法によって選択された平面を使用して再帰的に細分を行い、その後、階層のリーフを近接性の順序で列挙するために、階層を順序正しく横断することによって、構築することができる。より効率的な横断は、ツリーを左平衡化（ｌｅｆｔ−ｂａｌａｎｃｅ）させ、ツリーを配列データ構造に保存することによって提供することができる。空間階層を使用してソートするステップは、バケットソーティング及び選択された空間充填曲線を利用するステップと、本発明の一態様では、規則的ボクセルを使用するステップとを含むこともできる。
本発明の別の態様では、シミュレーション及び／又はソーティング処理は、
軸平行境界ボックス（ａｘｉｓ−ａｌｉｇｎｅｄ ｂｏｕｎｄｉｎｇ ｂｏｘ）によって、レンダリングされる物体の境界を定めるステップと、
物体の選択された点を通過する分割平面を適用することによって、境界ボックスを、各々がリーフで終了する、連続的な左ブランチ及び右ブランチに再帰的又は反復的に分割し、それによって、複数のリーフを生成するステップと、
リーフをそれぞれの明るさに従って階層的に順序づけるステップと、
選択された数のバケットに分割された、選択されたサイズの行列上で、バケットソートを実行するステップと、
行列内をたどるために選択された空間充填曲線を使用するステップであって、個々のバケット内において、より小さな空間充填曲線が、行列内の個々のセル内をたどるために使用される、ステップと
を含む。
本発明のシミュレーションプロセスは、積分方程式を解くための構成用に適合することができる。シミュレーションは、リアルに見える画像を合成するための光輸送（ｌｉｇｈｔ ｔｒａｎｓｐｏｒｔ）シミュレーションを可能にするように動作可能である。本発明のこの態様によれば、光輸送シミュレーションの基礎をなす積分方程式は、経路空間をサンプリングすることがマルコフ連鎖をシミュレートすることに対応するように、経路積分として定式化し直すことができ、その際、経路は、レイトレーシング及び散乱事象によって確立され、初期分布は、光源又はセンサの放出特性によって決定され、遷移確率は、画像内の表面についての双方向反射分布関数（ｂｉ−ｄｉｒｅｃｔｉｏｎａｌ ｒｅｆｌｅｃｔａｎｃｅ ｄｉｓｔｒｉｂｕｔｉｏｎ ｆｕｎｃｔｉｏｎ）及び双方向表面下散乱分布関数（ｂｉ−ｄｉｒｅｃｔｉｏｎａｌ ｓｕｂｓｕｒｆａｃｅ ｓｃａｔｔｅｒｉｎｇ ｄｉｓｔｒｉｂｕｔｉｏｎ ｆｕｎｃｔｉｏｎ）によって決定される。透過効果は、経路積分、初期分布、及び遷移確率を利用することによって、処理することができる。本発明の関連態様では、経路積分は、高次元準モンテカルロ点を利用することによって、又は低次元準モンテカルロ点をパディングすることによって、解くことができ、光輸送シミュレーションの基礎をなす積分方程式は、Ｆｒｅｄｈｏｌｍ積分方程式又はＶｏｌｔｅｒｒａ積分方程式と見なされる。
さらなる詳細、実施形態、実践、及び例が、添付の図面と併せて読まれる以下の詳細な説明において説明され、又は図面と併せて詳細な説明を読むことから、当業者に容易に明らかとなるであろう。
本発明が配備され得る、従来のデジタル処理システムの概略図である。 本発明が配備され得る、従来のパーソナルコンピュータ又は同様のコンピューティング装置の概略図である。 本発明の第１の態様による、全体的な方法を示す図である。 自己交差の問題を図説する、レイトレーシング手順の図である。 本発明のさらなる態様による、アクセラレーションデータ構造として使用される分割された軸平行境界ボックスの立面図による図である。 Ｌ平面及びＲ平面を用いた境界ボックスの分割を示す、図５に示された軸平行境界ボックスの等角図による図である。 Ｌ平面及びＲ平面を用いた境界ボックスの分割を示す、図５に示された軸平行境界ボックスの等角図による図である。 Ｌ平面及びＲ平面を用いた境界ボックスの分割を示す、図５に示された軸平行境界ボックスの等角図による図である。 本発明のさらなる態様による、レイトレーシング法のフローチャートである。 本発明のさらなる態様による、レイトレーシング法のフローチャートである。 シフトΦＳ（ｎ）を図説する図である。 例示的なＨａｌｔｏｎ列をその階層構造特性とともに示す図である。 双方向経路トレーシングによってサンプリング輸送経路空間を示す図である。 モンテカルロ法を使用するシミュレーションと配列ランダム化準モンテカルロ法を使用するシミュレーションの間の差を示した、レンダリングされた画像の図である。 モンテカルロ法を使用するシミュレーションと配列ランダム化準モンテカルロ法を使用するシミュレーションの間の差を示した、レンダリングされた画像の図である。 モンテカルロ法を使用するシミュレーションと配列ランダム化準モンテカルロ法を使用するシミュレーションの間の差を示した、レンダリングされた画像の図である。 モンテカルロ法を使用するシミュレーションと配列ランダム化準モンテカルロ法を使用するシミュレーションの間の差を示した、レンダリングされた画像の図である。 光源上の開始点からの光子軌道を示す図である。 光源上の開始点からの光子軌道を示す図である。 光源上の開始点からの光子軌道を示す図である。 光源上の開始点からの光子軌道を示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 状態を空間階層のリーフに分類し、それによって、階層を順序正しく横断することによって、近接性の順序を定義するステップを示す図である。 テストシーン「Ｓｈｉｒｌｅｙ ６」の図である。 マスタ解に対する平均ＲＭＳ誤差を示すグラフである。 大域的分散を示すグラフである。 ピクセルベースの分散を示すグラフである。 テストシーン「Ｉｎｖｉｓｉｂｌｅ Ｄａｔｅ」の図である。 マスタ解に対する平均ＲＭＳ誤差を示すグラフである。 大域的分散を示すグラフである。 ピクセルベースの分散を示すグラフである。 モンテカルロを使用してレンダリングされたテストシーン「Ｉｎｖｉｓｉｂｌｅ Ｄａｔｅ」の視覚的比較の図である。 ランダム化準モンテカルロを使用してレンダリングされたテストシーン「Ｉｎｖｉｓｉｂｌｅ Ｄａｔｅ」の視覚的比較の図である。 ソーティングを伴うランダム化準モンテカルロを使用してレンダリングされたテストシーン「Ｉｎｖｉｓｉｂｌｅ Ｄａｔｅ」の視覚的比較の図である。 迷宮テストシーンの概略図である。 各技法を使用した際のカメラによって受け取られた総放射の平均量を示すグラフである。 グラフの拡大部分の図である。 マスタ解からの偏差を表示する表である。 本発明の説明された態様による、技法及び副技法のフローチャートである。 本発明の説明された態様による、技法及び副技法のフローチャートである。 本発明の説明された態様による、技法及び副技法のフローチャートである。 本発明の説明された態様による、技法及び副技法のフローチャートである。 本発明の説明された態様による、技法及び副技法のフローチャートである。 本発明の説明された態様による、技法及び副技法のフローチャートである。 本発明の説明された態様による、技法及び副技法のフローチャートである。
本発明は、レイトレーシングのため、及び高速レイトレーシング用に必要とされるアクセラレーションデータ構造の効率的な構築のための、改良された技法を提供する。以下の説明は、これらの技法による方法、構造、及びシステムについて説明する。
説明される方法及びシステムが、Ｍｉｃｒｏｓｏｆｔ Ｗｉｎｄｏｗｓ、Ｌｉｎｕｘ、又はＵｎｉｘなどの従来のオペレーティングシステムに従って動作する、又は従来のオペレーティングシステムをエミュレートする、スタンドアロン構成の、又はネットワーク上に分散する、パーソナルコンピュータ（ＰＣ）又は等価なデバイスなどの従来のコンピュータ装置を使用して、ソフトウェア、ハードウェア、又はソフトウェアとハードウェアの組合せで実施され得ることは、当業者であれば理解されよう。以下に説明され、特許請求の範囲で列挙される、様々な処理手段及び計算手段は、したがって、適切に構成されたデジタル処理デバイス又はデバイスのネットワークの、ソフトウェア及び／又はハードウェア要素で実施されてよい。処理は、順次又は並列に実行されてよく、専用又は再構成可能なハードウェアを使用して実施されてよい。
本発明による方法、デバイス、又はソフトウェア製品は、スタンドアロンであるか、ネットワーク接続されているか、ポータブルであるか、又は固定式であるかを問わない、図１（例えばネットワークシステム１００）に例として示されたような、広範な従来のコンピューティングデバイス及びシステムのいずれかにおいて動作することができ、従来のコンピューティングデバイス及びシステムは、従来のＰＣ１０２、ラップトップ１０４、ハンドヘルド若しくはモバイルコンピュータ１０６を含み、又はインターネット若しくは他のネットワーク１０８を介して接続され、インターネット若しくは他のネットワーク１０８は、同様にサーバ１１０及びストレージ１１２を含むことができる
従来のコンピュータソフトウェア及びハードウェア実施に沿って、本発明に従って構成されたソフトウェアアプリケーションは、例えば、図２に示されるようなＰＣ１０２内で動作することができ、ＰＣ１０２において、プログラム命令は、ＣＤ−ＲＯＭ１１６、磁気ディスク、又は他のストレージ１２０から読み出され、ＣＰＵ１１８による実行のためにＲＡＭ１１４にロードされることができる。データは、従来のキーボード、スキャナ、マウス、又は他の要素１０３を含む、任意の知られたデバイス又は手段を介して、システムに入力することができる。
本発明の本説明は、２つの大きなセクションに分割される。
Ｉ．リアルタイム高精度レイトレーシング
ＩＩ．マルコフ連鎖及び準モンテカルロ法を使用してシミュレートされる光線のためのソーティング戦略
Ｉ．リアルタイム高精度レイトレーシング
図３は、本発明の一態様による全体的な方法２００を示す図である。この方法は、コンピュータグラフィックスシステムに関連して実施され、システムでは、画像内の各ピクセルについて、ピクセル値が生成される。生成された各ピクセル値は、シーン内の点を、シミュレートされるカメラの画像平面上に記録されるように表す。コンピュータグラフィックスシステムは、選択されたレイトレーシング法を使用して、画像のピクセル値を生成するように構成される。選択されたレイトレーシング法は、選択された方向に沿ってピクセルからシーンに向かって放たれた少なくとも１つの光線を含む光線ツリー（ｒａｙ ｔｒｅｅ）の使用を含み、光線とシーン内の物体（及び／又は物体の表面）との交わりの計算をさらに含む。
図３の方法２００において、光線とシーン内の表面との交わりを計算するために、境界ボリューム階層が使用される。ステップ２０１において、シーンの境界ボックスが計算される。ステップ２０２において、所定の終了基準が満たされているかどうかが判定される。満たされていない場合、ステップ２０３において、軸平行境界ボックスが精緻化される。プロセスは、終了基準が満たされるまで、再帰的に継続される。本発明の一態様によれば、終了基準は、境界ボックス座標が光線／表面の交点の浮動小数点表現から解像度（ｒｅｓｏｌｕｔｉｏｎ）の１単位分だけ異なるという条件として定義される。しかし、本発明の範囲は、他の終了基準にまで及んでいる。
アクセラレーション構造としての境界ボリューム階層の使用は、多くの理由で有利である。境界ボリューム階層のためのメモリ要件は、レイトレーシングされる物体の数に線形に制限され得る。また、以下に説明されるように、境界ボリューム階層は、３Ｄツリーよりもはるかに効率的に構築することができ、このことが、全体的に動きのあるシーンに必要とされる償却解析（ａｍｏｒｔｉｚｅｄ ａｎａｌｙｓｉｓ）にとって、境界ボリューム階層を非常に適したものにする。
以下の説明は、レイトレーシング法におけるある種の問題と、それらの問題に対処する本発明の特定の態様について詳細に説明する。図４は、「自己交差（ｓｅｌｆ−ｉｎｔｅｒｓｅｃｔｉｏｎ）」問題を図説する図である。図４は、画像表面３０２と、観測点３０４と、光源３０６とを含む、レイトレーシング手順３００を示している。表面の画像を合成するため、一連の計算が、観測点３０４と表面３０２の間に延在する光線を見出すために実行される。図４は、１つのそのような光線３０８を示している。理想的には、その場合、光線３０８と表面３０２との正確な交点３１０が計算される。
しかし、コンピュータにおける浮動小数点算術計算のせいで、計算された光線／表面の交点３１２が実際の交点３１０とは異なることが時にあり得る。さらに、図４に図説されるように、計算された点３１２が表面３０２の「誤った」側に見出されることがあり得る。その場合、計算された光線／表面の交点３１２から光源３０６まで延在する第２の光線３１４を見出すために計算が実行されたとき、これらの計算は、第２の光線３１４が、光源３０６まで直接的に延在する代わりに、第２の交点３１６において表面３０２に打ち当たることを示し、したがって、画像誤差をもたらす。
自己交差問題に対する１つの知られた解決策は、表面３０２から安全な距離のところから各光線３０８を出発させることである。この安全な距離は一般に、グローバル浮動小数点（ｇｌｏｂａｌ ｆｌｏａｔｉｎｇ ｐｏｉｎｔ）εとして表現される。しかし、グローバル浮動小数点εの決定は、画像が合成される、シーン及びシーン自体の中の特定の位置に大きく依存する。
本発明の一態様は、より高精度の代替案を提供する。計算された光線／表面の交点３１２に到達した後、実際の交点３１０により近い交点を再計算するために、計算された点３１２及び光線３０８の方向が使用される。交点のこの再計算は、精度を高める反復として、レイトレーシング法に組み込まれる。反復的に計算される交点が表面３０２の「誤った」側に存在すると判明した場合、その交点は、表面３０２の「正しい」側に移動される。反復的に計算される交点は、表面の法線に沿って、又は法線の最長成分によって決定される軸に沿って移動することができる。グローバル浮動小数点εを使用する代わりに、小数点は浮動小数点の仮数の最終ビットの方へ整数εだけ移動される。
説明された手順は、倍精度の計算を回避し、浮動小数点数の指数によって決定される浮動小数点数の位取り（ｓｃａｌｅ）に暗黙に適応するという利点を有する。したがって、この実施では、すべての２次光線は、εオフセットを不要にするこれらの修正された点から直接出発する。したがって、交わり計算の最中、有効な光線区間は、何らかのオフセットの代わりに、０から開始すると仮定することができる。
仮数の整数表現の修正は、どの点がどの側に存在するかを決定するために三角形と平面とを交わらせる場合、数値問題も回避する。
凸結合の凸閉包性を活用して、光線と自由曲面との交わりは、光線の源に最も近い交点を含む軸平行境界ボックスを精緻化することによって見出すことができる。この精緻化は、浮動小数点数の解像度に達するまで、即ち、境界ボックス座標が浮動小数点表現から解像度の１単位分だけ異なるまで、継続することができる。その場合、自己交差問題は、境界ボックスの中心で表面法線に最も近い境界ボックスの角を選択することによって回避される。この角点は、その後、第２の光線を出発させるために使用される。この「光線物体交わりテスト」は、非常に効率的であり、自己交差問題の回避から利益を得る。
アクセラレーションデータ構造を構築した後、三角形が面内で変形される。新しい表現は、特別な苦労なしに交わりテストが縮退三角形（ｄｅｇｅｎｅｒａｔｅ ｔｒｉａｎｇｌｅ）を扱えるように、縮退三角形を符号化する。新しい表現は、もちろん、縮退三角形がグラフィックスパイプラインに入ることを単に防止することも可能である。コンピュータプログラムリストのセクション２．２．１及び２．２．２は、本発明の本態様によるソースコードのリストを示している。
テストは最初に、光線と三角形の平面との交わりを決定し、その後、光線上の有効区間］０，ｒｅｓｕｌｔ．ｔｆａｒ］の外部の交わりを除外する。これは、ただ１回の整数テストによって達成される。＋０が有効区間から除外されることに留意されたい。これは、非正規化浮動小数点数がサポートされない場合に重要である。この最初の決定が成功した場合、交わりの重心座標を計算することによって、テストは先に進む。完全包含テストを実行するために、やはり整数テストだけが、即ち、より詳細には２つのビットのテストだけが必要とされることに留意されたい。したがって、ブランチの数は最小である。この効率的なテストを可能にするため、三角形の辺及び法線は、変形ステップにおいて、適切にスケーリングされる。
テストの精度は、誤った又は失われた光線交わりを回避するのに十分である。しかし、横断の最中、堅牢な交わりテストのために、三角形を拡大するのが適切である状況が生じることがある。これは、三角形を変形する前に行うことができる。三角形はその法線の最長成分によって識別される軸に沿って射影されるので、この射影ケースは、保存されなければならない。これは、アクセラレーションデータ構造のリーフノードにおいて、カウンタによって達成される。三角形参照（ｔｒｉａｎｇｌｅ ｒｅｆｅｒｅｎｃｅ）は、射影ケースによってソートされ、リーフは、各クラス内の三角形の数のためのバイトを含む。
本発明のさらなる態様は、レイトレーシング用のアクセラレーションデータ構造を構築するための改良された手法を提供する。多くの異なる最適化に従った従来のソフトウェア実施と比較して、本明細書で説明される手法は、優れたレイトレーシング性能をもたらす著しく平たいツリーを生みだす。
分割平面の候補は、分割される軸平行境界ボックスの内部の三角形頂点の座標によって与えられる。これは、実際には境界ボックスの外部にあるが、境界ボックスによって定義される３つの区間の１つに含まれる座標を少なくとも１つ有する、頂点を含むことに留意されたい。これらの候補から、現在の軸平行境界ボックスの最長辺の中央に最も近い平面が選択される。さらなる最適化は、表面法線の最長成分が潜在的な分割平面の法線と一致する三角形の候補のみを選択する。分割平面を三角形頂点を通るように置くことは、分割平面によって分割される三角形の数を暗黙的に減少させるので、この手順は、はるかに平たいツリーを生みだす。加えて、表面が緊密に近似され、空の空間が最大化される。三角形の数が指定された閾値よりも多くなり、分割平面の候補がもはや存在しない場合、ボックスは、その最長辺に沿って中央で分割される。これは、例えば、長対角線物体の使用を含む、他の手法の非効率性を回避する。
どの三角形が階層内のノードの左側及び右側チャイルドに属するかを決定する再帰的手順は一般に、広範囲のブックキーピングと、メモリ割り当てとを必要とした。例外ケースにだけうまくいかない、はるかに簡略な手法が存在する。レイトレーシングされる物体への参照の配列が、２つだけ割り当てられる。第１の配列は、物体参照（ｏｂｊｅｃｔ ｒｅｆｅｒｅｎｃｅ）を用いて初期化される。再帰的な空間分割の最中、左側の要素のスタックは、配列の先頭から成長し、一方、右側に分類される要素は、配列の末尾から中央に向かって成長するスタック上に維持される。分割平面と交わっている要素、即ち、左側でも右側でもある要素を速やかに回復できるように、第２の配列が、それらのスタックを維持する。したがって、バックトラッキングは、効率的で簡単である。
表面積ヒューリスティック（ｓｕｒｆａｃｅ ａｒｅａ ｈｅｕｒｉｓｔｉｃ）を使用することによってツリーのブランチを刈り込む代わりに、ツリーの深さは、固定された深さから開始するバイナリ空間分割を近似的に左平衡化することによって刈り込まれる。網羅的な実験によって観測されたように、グローバルな固定深さパラメータが、非常に多様なシーンにわたって指定され得る。これは、ある回数のバイナリ空間分割の後、空間において相対的に平らな接続された要素が通常は残っていることを観測することによって理解できる。コンピュータプログラムリストのセクション２．３．１は、本発明のこの態様によるソースコードのリストを示している。
境界ボリューム階層を使用して、レイトレーシングされる各物体は、正確に１度だけ参照される。その結果、３Ｄツリーとは対照的に、階層の横断の最中に、物体と光線との複数の交わりを防止するために、メールボックスメカニズム（ｍａｉｌｂｏｘ ｍｅｃｈａｎｉｓｍ）は必要とされない。これは、システム性能の観点から見て大きな利点であり、共用メモリシステム上での実施をはるかに簡単にする。第２の重要な結果は、境界ボリューム階層のツリーでは、レイトレーシングされる物体の総数よりも多い内部ノードが存在できないことである。したがって、アクセラレーションデータ構造のメモリフットプリントは、構築前に物体の数に線形に制限され得る。そのような事前の制限は、３Ｄツリーの構築では利用可能でなく、３Ｄツリーでは、メモリ計算量は、レイトレーシングされる物体の数につれて２次的に増大することが予想される。
したがって、現在の３Ｄツリーレイトレーシング法よりも著しく高速であり、レイトレーシングされる物体の数につれて、メモリ要件が２次的に予想されるのではなく線形に増大する、境界ボリューム階層の新しい概念が、今から説明される。境界ボリューム階層が３Ｄツリーよりも性能的にまさることを可能にする中核概念は、境界ボリューム自体に焦点を絞る代わりに、空間がどのように分割され得るかに焦点を絞ることである。
３Ｄツリーにおいては、境界ボックスは、単一の平面によって分割される。本発明の本態様によれば、２つの軸平行境界ボックスを定義するために、２つの平行な平面が使用される。図５は、主要なデータ構造４００を示す図である。
図５は、軸平行境界ボックス４０２を立面図で示している。軸平行であり、互いに平行である、Ｌ平面４０２及びＲ平面４０４は、境界ボックス４０２を左軸平行境界ボックスと右軸平行境界ボックスとに分割するために使用される。左境界ボックスは、元の境界ボックス４０２の左壁４０６からＬ平面４０２まで広がる。右境界ボックスは、Ｒ平面４０４から元の境界ボックス４０２の右壁４０８まで広がる。したがって、左境界ボックスと右境界ボックスとは、互いに重なり合うことがある。光線４１２の横断は、光線４１０の有効な区間［Ｎ，Ｆ］４１２に関するＬ平面４０２及びＲ平面４０４の交わりの位置によって決定される。
図５のデータ構造４００では、Ｌ平面４０２及びＲ平面４０４は、境界ボックス４０２内に含まれる空間ではなく、元の境界ボックス４０２内に含まれる物体の集合を分割するように、互いに対して位置づけられる。３Ｄツリー分割とは対照的に、２つの平面をもつことが、２つの平面の間の空の空間を最大化する可能性を提供する。その結果、シーンの境界が、はるかに高速に近似され得る。
図６〜図８は、データ構造４００をさらに説明する、一連の３次元図である。図６は、境界ボックス４０２の図を示している。説明の目的で、境界ボックス４０２内に仮想物体が、抽象的な円４１６として示されている。図７及び図８に示されるように、境界ボックス４０２を左境界ボックス４０２ａと右境界ボックス４０２ｂとに分割するために、次にＬ平面４０４及びＲ平面４０６が使用される。Ｌ平面及びＲ平面は、それらの間の空の空間が最大化されるように選択される。各仮想物体４１６は、最終的には左境界ボックス４０２ａ又は右境界ボックス４０２ｂ内に入る。図８の下図に示されるように、仮想物体４１６は、「左」物体４１６ａと「右」物体４１６ｂとに分割される。結果の境界ボックス４０２ａ、４０２ｂの各々は、それ自体が分割され、それ以降も、終了基準が満たされるまで同様である。
図９は、説明される方法５００のフローチャートである。ステップ５０１において、シーンの境界ボックスが計算される。ステップ５０２において、軸平行境界ボックスを、重なり合うことがある左軸平行境界ボックスと右軸平行境界ボックスとに分割するために、Ｌ平面及びＲ平面が使用される。ステップ５０３において、元の軸平行境界ボックスに含まれる仮想物体の集合を、左物体の集合と右物体の集合とに分割するために、左境界ボックス及び右境界ボックスが使用される。ステップ５０４において、左物体及び右物体が、終了基準が満たされるまで、再帰的に処理される。
以前の実施において使用された１つの分割パラメータの代わりに、２つの分割パラメータが、ノード内に保存される。ノードの数はレイトレーシングされる物体の数によって線形的に制限されるので、すべてのノードの配列は、１度だけ割り当てられ得る。したがって、構築中の３Ｄツリーの高コストのメモリ管理が不要になる。
構築技法は、３Ｄツリー構築用の類似物（ａｎａｌｏｇ）よりはるかに簡単であり、再帰的方法で、又は反復バージョン及びスタックを使用することによって、容易に実施される。物体のリスト及び軸平行境界ボックスを与えると、Ｌ平面及びＲ平面が決定され、物体の集合がしかるべく決定される。その後、左物体及び右物体が、何らかの終了基準が満たされるまで、再帰的に処理される。内部ノードの数は制限されているので、ただ１つの物体が残されているときに、終了を当てにしても安全である。
分割が、ｘ軸、ｙ軸、及びｚ軸に垂直な平面に沿って、物体を分類することだけに依存しており、それは非常に効率的で、数値的に絶対安定であることに留意されたい。３Ｄツリーとは対照的に、物体と分割平面との正確な交わりは計算される必要がなく、それは、数値的に堅牢な方法で達成するには、より大きなコストがかかり困難である。頂点及び辺沿いでの三角形の喪失など、３Ｄツリーの数値問題は、境界ボリューム階層の構築前に三角形を拡大することによって回避され得る。また、３Ｄツリーでは、重なり合う物体は、左軸平行境界ボックスと右軸平行境界ボックスの両方に分類されなければならず、それによって、ツリーの予想される２次成長を引き起こす。
Ｌ平面及びＲ平面を、したがって実際のツリーレイアウトを決定するために、様々な技法が使用されてよい。図６〜図８に戻ると、１つの技法は、上述の３Ｄツリー構築技法を使用して平面Ｍ４１８を決定し、新しい軸平行境界ボックスの結果のＬ平面及びＲ平面の重なりが、提案された分割平面Ｍ４１８と最も小さく重なり合うように、物体を分割することである。結果のツリーは、対応する３Ｄツリーと非常に類似しているが、空間の代わりに、物体集合が分割されるので、結果のツリーは、はるかに平たい。別の手法は、チャイルドボックスの重なりが最小になり、空の空間が可能な限り最大化されるような方法で、Ｒ平面及びＬ平面を選択することである。ある物体にとっては、軸平行境界ボックスが非効率的であることに留意されたい。そのような状況の一例は、軸平行境界ボックスの対角線上の、半径の小さい長い円柱である。
図１０は、本発明のこの態様による、方法６００のフローチャートである。ステップ６０１において、シーンの境界ボックスが計算される。ステップ６０２において、分割平面Ｍを決定するために、３Ｄツリー構築が実行される。ステップ６０３において、軸平行境界ボックスを、分割平面Ｍと最も小さく重なり合う左軸平行境界ボックスと右軸平行境界ボックスとに分割するために、平行なＬ平面及びＲ平面が使用される。ステップ６０４において、元の軸平行境界ボックスに含まれる仮想物体の集合を、左物体の集合と右物体の集合とに分割するために、左境界ボックス及び右境界ボックスが使用される。ステップ６０５において、左物体及び右物体が、終了基準が満たされるまで、再帰的に処理される。図１０で説明される方法６００は、図３で説明された方法２００と同様に、３Ｄツリー構築、リアルタイム処理、バケットソーティング、自己交差などに関する技法を含む、本明細書で説明される他の技法と組み合わされてよいことに留意されたい。
３Ｄツリーの場合、空間細分は、物体周囲の空間の空の部分を切り取るように継続される。説明される境界ボリューム階層の場合、そのような物体をより小さな物体に分割することが、同様の挙動をもたらす。メモリ要件の予測可能性を維持するため、最大境界ボックスサイズが定義される。最大境界ボックスサイズを超える大きさをもつすべての物体は、要件を満たすように、より小さい部分に分割される。最大許容サイズは、すべての物体のうちの最小の大きさを求めてデータセットをスキャンすることによって見出すことができる。
本明細書で説明されるデータ構造は、高速３Ｄツリー横断の原理を境界ボリューム階層に移植することを可能にする。（１）左チャイルドのみ、（２）右チャイルドのみ、（３）左チャイルドの後、右チャイルド、（４）右チャイルドの後、左チャイルド、又は（５）光線が分割平面の間（即ち空の空間）にある、という、横断のケースは同様である。説明される技法における１つのノードは２つの平行な平面によって分割されるので、ボックスをどのような順序で横断するかは、光線の方向によって決定される。図６は、上述の技法を含むソースコードリストを説明している。
これまでの境界ボリューム階層技法は、チャイルドノード、又はヒープデータ構造の更新など必要とされる付加的な作業をどのような順序で横断するかを効率的に決定できなかった。加えて、全境界ボリュームが、ロードされ、光線に対してテストされなければならなかったが、本手法は、２つの平面距離を必要とするだけである。しかし、ソフトウェアで２つの平面に対して光線をチェックすることは、よりコスト高であるように思える。横断は、３Ｄツリーにおけるボトルネックであり、もう少し多くの計算をここで行うことが、メモリアクセスの待ち時間をより良く隠蔽する。加えて、境界ボリューム階層ツリーは、同じ性能の対応する３Ｄツリーよりもはるかに小さくなる傾向にある。
本明細書では新しい境界ボリューム階層が説明されるが、３Ｄツリーの横断と強い関連が存在する。Ｌ＝Ｒと設定すると、従来のバイナリ空間分割が獲得され、横断技法は、３Ｄツリーの横断技法へと崩落する。
説明される境界ボリューム階層は、自由曲面を細分することによって、光線自由曲面交わりを効率的に見つけるためにも適用することができる。そのようにすることは、自由曲面と凸閉包性との交わり、及び実際の浮動小数点計算に応じて浮動小数点精度まで効率的に計算される細分技法を可能にする。細分ステップは、例えば、多項式曲面、有理曲面、及び近似細分化曲面（ａｐｐｒｏｘｉｍａｔｉｎｇ ｓｕｂｄｉｖｉｓｉｏｎ ｓｕｒｆａｃｅ）について実行される。空間内の各軸について、重なり合う可能性のある境界ボックスが、上で説明されたように決定される。バイナリ細分の場合、新しいメッシュの新しい境界ボックスについてのＬボックスの交わり及びＲＬボックスの交わり。今は、ボックスの空間的順序が知られているので、上述の横断が効率的に実行できる。境界ボリュームの階層を事前計算する代わりに、それはオンザフライで計算することができる。この手順は、自由曲面について効率的であり、アクセラレーションデータ構造のためのメモリを節約することを可能にし、それは、バックトラッキングによって横断されなければならない境界ボリュームの小さなスタックによって置き換えられる。細分は、光線表面交わりが浮動小数点精度の点又は小さなサイズの区間へと崩落した境界ボリュームに入るまで、継続される。コンピュータプログラムリストのセクション２．１．１は、本発明のこの態様によるコードリストを示している。
レイトレーシングにおけるアクセラレーションデータ構造としての規則的グリッドの使用は簡単であるが、効率は、空間適応性の欠如、及び多くの空のグリッドセルのその後の横断により影響をこうむる。階層的な規則的グリッドは、状況を改善できるが、境界ボリューム階層及び３Ｄツリーと比較すると依然として劣っている。しかし、アクセラレーションデータ構造の構築速度を改善するために、規則的グリッドが使用できる。アクセラレーションデータ構造を構築するための技法は、クイックソーティングと類似しており、○（ｎｌｏｇｎ）で動作することが予想される。線形時間で動作するバケットソーティングを適用することによって、改善が獲得され得る。したがって、物体の軸平行境界ボックスは、ｎｘ×ｎｙ×ｎｚ個の軸平行ボックスに分割される。その後、各物体は、１つの選択点によって、これらのボックスの１つに正確に分類され、例えば、重心又は各三角形の第１の頂点が使用できる。その後、各グリッドセル内の物体の実際の軸平行境界ボックスが決定される。これらの軸平行境界ボックスは、ボックスが分割平面の１つと交わらない限り、それらが含む物体の代わりに使用される。交わる場合は、ボックスが開けられ、代わりに、ボックス内の物体が直接使用される。この手順は、多くの比較及びメモリアクセスを省き、構築技法のオーダの定数を目覚しく改善し、再帰的に適用することもできる。上記の技法は、物体のストリームを処理することによって実現され得るので、特にハードウェア実装にとって魅力的である。
アクセラレーションデータ構造は、要求時に、即ち、物体を有する特定の軸平行境界ボックスを光線が横断したときに、構築することができる。その場合、一方では、アクセラレーションデータ構造は、光線の射さない空間の領域では、精緻になることはなく、キャッシュは、決して触れられることのないデータによって汚染されない。他方、精緻化の後、光線と交わる可能性のある物体は、すでにキャッシュ内に存在する。
上記の説明から、本発明が、レイトレーシングにおける長く知られていた問題に対処し、改善された精度、全体的速度、及びアクセラレーションデータ構造のメモリフットプリントを改善した、レイトレーシングのための技法を提供することが分かる。数値精度の改善は、他の数体系への転換ばかりでなく、例えば、ＡＲＴレイトレーシングチップのハードウェアで使用される対数体系への転換も行う。プロセッサ又は専用ハードウェアへのＩＥＥＥ浮動小数点標準の具体的な実装は、性能に深刻な影響を及ぼし得ることに留意されたい。例えば、Ｐｅｎｔｉｕｍ ４チップ上で、非正規化数は、１００分の１以下に性能を低下させ得る。上で説明されたように、本発明の実施は、これらの例外を回避する。本明細書で説明される境界ボリューム階層の概念は、境界ボリューム階層をリアルタイムレイトレーシングにとって適したものにする。償却解析においては、説明された技法は、当技術分野のこれまでの状況よりも性能が優れており、したがって、制作セッティング（ｐｒｏｄｕｃｔｉｏｎ ｓｅｔｔｉｎｇ）などにおいて、例えば、全体的に動きのあるシーンでモーションブラーを計算するための技法など、より正確な技法の使用を可能にする。説明された境界ボリューム階層が、３Ｄツリー及び他の技法と比較した場合、特にハードウェア実装において、また巨大シーンについて、著しい利点を有することは、上記の説明から明らかである。償却解析においては、説明された境界ボリューム階層は、現在の３Ｄツリーよりも少なくとも２倍、性能が優れている。加えて、メモリフットプリントは、事前に決定することができ、物体の数に対して線形である。
ＩＩ．マルコフ連鎖及び準モンテカルロ法を使用してシミュレートされる光線のためのソーティング戦略
以下に説明される本発明の態様は、光線表面交わりを計算するためにレイトレーシング法を使用する、画像内のピクセルのピクセル値を計算するためにコンピュータ又は他のデジタルプロセッサが使用されるコンピュータグラフィックス及び他の応用例における使用に適した、光輸送の効率的なフォトリアリスティックシミュレーションのための方法、技法、及びシステムを提供する。
特に、本発明は、光子など光のパケット又はユニットによって照明されるシーン内の物体の高速処理を可能にするソーティング戦略を適用することによって、光線表面交わりテストをより効率的にするための方法を含む、従来のレイトレーシング法に対する改良を提供し、光のパケット又はユニットの軌道は、マルコフ連鎖を使用してシミュレートされ、マルコフ連鎖は、適切な準モンテカルロ（ＱＭＣ）法を使用してシミュレートされる。ソーティング戦略は、近接性ソーティングを含み、本発明は、有利には、そのような近接性ソーティングを使用する低次元のＱＭＣ法を活用する。
説明される技法及びシステムは、従来のシステムでは非常な困難を提示した課題であった、リアルタイムに近い時間で積分方程式を解くことのために使用することができる。
これまでは、マルコフ連鎖、例えば光子の軌道をシミュレートするための準モンテカルロ法の使用は、ランダムサンプリングよりわずかに良い程度の改善を提供したにすぎなかった。しかし、本発明によるＱＭＣ法は、解演算子の本質的な低次元構造を活用することができ、はるかに効率的な技法をもたらす。
本発明のこの態様の説明は、以下のように構成される。
１．マルコフ連鎖及び準モンテカルロ法
２．マルコフ連鎖の同時シミュレーション
２．１ １次元での決定論的技法の分析
２．１．１ 技法の簡略化
２．１．２ いつ、なぜ、それが機能するか
２．２ ｓ次元での簡略化技法
２．３ ランダム化
２．４ 配列ランダム化準モンテカルロ
３．光輸送シミュレーションへの応用
３．１ Ｆｒｅｄｈｏｌｍ積分及びＶｏｌｔｅｒｒａ積分
４．ソーティング技法及びシステム
４．１ 状態のノルム
４．２ 空間階層
４．３ バケットソーティング及び空間充填曲線
５．結果
６．結論
７．一般的技法
１．マルコフ連鎖及び準モンテカルロ法
数学的に、マルコフ連鎖は、システム状態の離散時間確率系列である。マルコフ連鎖は、現在の状態の知識が与えられれば、以降の状態の確率を予想する上で、以前の状態は無関係である、「記憶いらずの」プロセスである。マルコフ連鎖は、多くの異なる現象をモデル化し、分析するために使用することができる。例えば、英語は、１つの単語が別の単語の後に続く相対頻度の計算に基づいて、マルコフ連鎖を使用してモデル化された。これらの遷移確率を使用して、英語の文を生成することが可能である。例えば、待ち行列、ブラウン運動、粒子輸送の分析、及び本説明に特に関係が深い光輸送シミュレーションを含む、他の多くの進展プロセスも、マルコフ連鎖を使用して説明することができる。
より具体的には、マルコフ連鎖は、シミュレート画像のシーン内の物体を照明する光子の軌道をシミュレートするために使用できる。光輸送シミュレーションは、複雑な領域及び不連続性と、高次元積分への還元とによって特徴づけられる。加えて、光輸送シミュレーションは、同一演算子の反復適用による低次元構造によって特徴づけられる。
マルコフ連鎖自体も、準モンテカルロ（ＱＭＣ）法を使用して、シミュレートされ得る。ＱＭＣ法の背後にある一般的な考えは、何らかの方法でランダムに生成された数の調整を試みる代わりに、選択区間の非常に一様な被覆（ｃｏｖｅｒａｇｅ）を提供するために、知られた決定論的系列から開始するというものである。ランダム化ＱＭＣでは、数はランダムな様相を有するが、規則的な方法で選択区間を被覆するという所望の特性をグループとして保持するような方法で、系列はランダム化される。
ＱＭＣを使用してマルコフ連鎖をシミュレートするために利用できる多くの手法が存在する。１つの手法は、高次元低食い違い量点集合（ｈｉｇｈ−ｄｉｍｅｎｓｉｏｎａｌ，ｌｏｗ−ｄｉｓｃｒｅｐａｎｃｙ ｐｏｉｎｔ ｓｅｔ）を適用することである。しかし、この手法は、中次元に対してだけ著しい利点を提供する。別の手法は、低次元ランダム化点集合のパディングである。この第２の手法は、実施するのがより簡単であり、ほぼ良好な結果を達成する。第３の手法は、状態のソーティングとパディングとを組み合わせる。この第３の手法は、ＱＭＣの能力が目に見える、性能の上昇をもたらし、決定論的技法及びランダム化技法の両方と併せて使用することができる。この手法の実施は、相対的に容易である。
プロセスの特性は、複数のマルコフ連鎖の寄与を平均することによって推定することができる。各軌道を個々にシミュレートする代わりに、相関サンプルを使用してマルコフ連鎖を同時にシミュレートし、各遷移ステップの後に１群の状態をソートすることで、収束を顕著に改善できることが見出された。
状態空間上で順序を与えた場合、これらの手法は、１群のマルコフ連鎖をソートすることによって改善され得る。以下に説明される本発明の態様によれば、決定論的手法の使用は、説明される技法において簡略化をもたらし、説明されたソーティングがいつ、なぜ、増大した性能をもたらすかに関する有益な徴候を提供する。以下にさらに説明されるように、説明されるソーティング戦略の一部又は全部は、有利には、光輸送シミュレーションの効率を増大させるために使用することができる。
２．マルコフ連鎖の同時シミュレーション
効率を増大させるため、マルコフ連鎖は、同時にシミュレートされ得る。準モンテカルロ法を用いたマルコフ連鎖の同時シミュレーションを中間ソーティングステップによって改善するというアイデアは、ボルツマン方程式を扱い、後に熱方程式を扱った一連の論文において、最初にＬｅｃｏｔによって紹介された。その後、このアイデアは、グリッド上で熱方程式を解くために、Ｍｏｒｏｋｏｆｆ及びＣａｆｌｉｓｃｈによって使用及び洗練され、最近、当該技法のランダム化バージョンと、希少事象シミュレーションのための分割とを組み込むために、Ｌ’Ｅｃｕｙｅｒ、Ｔｕｆｆｉｎ、Ｄｅｍｅｒｓ他によって拡張された。独立した研究が、上記の手法に関連した方法で、コンピュータグラフィックスの分野において行われてきた。
以下の説明では、上記の方式にまさる改善を実施する技法が説明される。これらの改良技法の背景は、Ｈ．Ｎｉｅｄｅｒｒｅｉｔｅｒ編、Ｍｏｎｔｅ Ｃａｒｌｏ ａｎｄ Ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ Ｍｅｔｈｏｄｓ ｉｎ Ｓｃｉｅｎｔｉｆｉｃ Ｃｏｍｐｕｔｉｎｇ ２００２、ｐｐ．３２９〜４４、Ｓｐｒｉｎｇｅｒ Ｖｅｒｌａｇ（２００４）に記載の、Ｌｅｃｏｔ他、「Ｑｕａｓｉ−Ｍｏｎｔｅ Ｃａｒｌｏ Ｍｅｔｈｏｄｓ ｆｏｒ Ｅｓｔｉｍａｔｉｎｇ Ｔｒａｎｓｉｅｎｔ Ｍｅａｓｕｒｅｓ ｏｆ Ｄｉｓｃｒｅｔｅ Ｔｉｍｅ Ｍａｒｋｏｖ Ｃｈａｉｎｓ」において説明されている技法によって部分的に提供され、同文献は、その全体があたかも説明されているかのように、参照により本明細書に組み込まれる。
以下の説明は、説明される技法が、いつ、なぜ、従来の手法よりも優れたものとなるかに関する、概念的フレームワークを提供する。説明される技法の有利な実施も説明される。
２．１ １次元での決定論的技法の分析
１次元問題についての決定論的技法が今から説明される。Ｌｅｃｏｔにおいて提示された技法は、初期分布
Figure 2010501100

及び遷移行列
Figure 2010501100

を有する離散状態空間Ｅ上で、マルコフ連鎖を同時にシミュレートする。基数ｂの（ｔ，２）−列
Figure 2010501100

を使用して、Ｎ＝ｂｍ個の連鎖が、並列にシミュレートされ、ここで、Ｘｎ，ｌは、時間ステップｎにおける連鎖ｌの状態であり、
Figure 2010501100

は、ゼロを含む自然数の集合である。さらに、当該技法は、ｍ＞ｔについて、Ｎ＝ｂｍ個の後続要素ｘi，１が、基数ｂの（０，ｍ，１）−ネットを形成することを要求する。Ｌｅｃｏｔに示されるように、以下の基準が満たされる場合、当該技法は収束する。
Figure 2010501100
Ｘｎ，ｌは、時間ステップｎにおける連鎖ｌの状態である。したがって、シミュレーションは、以下のように実行され得る。
ｎ：＝０
０≦ｌ＜Ｎについて、ｘｌ，２を使用して、
Figure 2010501100

を初期化する。
ループ
何らかの順序を使用して、状態ベクトルをソートする。
ｎ：＝ｎ＋１
０≦ｌ＜Ｎについて、
Figure 2010501100

を使用して、
Figure 2010501100

を計算することによって、連鎖を継続させる。
Ｎ＝ｂｍ個の連鎖の場合、ｍ＞ｔについて、ｂｍ個の後続する要素ｘｉ，１が、（０，ｍ，１）−ネットを形成するように、基数ｂの（ｔ，２）−列 ｘｉが選択される。
基数ｂ＝２のＳｏｂｏｌ’（０，２）−列（φ２（ｌ），φＳ（ｌ））が、以下の要件を満たすことが分かる。
ｍ＝３＞ｔ＝０について、Ｎ＝２３＝８となり、
Figure 2010501100

について、
Figure 2010501100

となる。
すべての置換は同一である。
２．１．１ 技法の簡略化
一様分布を仮定した場合、安定なソーティングは、状態の順序を変化させないことが分かる。実際、収束条件は無変化であり続けるので、すべての置換は省略することができる。したがって、上述の技法は、簡略化することができる。
基数ｂ＝２の（０，２）−列であり、収束条件を維持するために必要とされる仮定を満たす、Ｓｏｂｏｌ’列
Figure 2010501100

について検討する。元のＳｏｂｏｌ’列は、Ｉ．Ｓｏｂｏｌ’、「Ｏｎ ｔｈｅ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ Ｐｏｉｎｔｓ ｉｎ ａ Ｃｕｂｅ ａｎｄ ｔｈｅ Ａｐｐｒｏｘｉｍａｔｅ Ｅｖａｌｕａｔｉｏｎ ｏｆ Ｉｎｔｅｇｒａｌｓ」、Ｚｈ．Ｖｙｃｈｉｓｌ．Ｍａｔ．ｉＭａｔ．Ｆｉｚ．、７（４）：７８４−８０２（１９６７）において定義されており、同文献は、その全体があたかも説明されているかのように、参照により本明細書に組み込まれる。
Figure 2010501100

についての簡単なコード例は、Ｔ．Ｋｏｌｌｉｇ及びＡ．Ｋｅｌｌｅｒ、「Ｅｆｆｉｃｉｅｎｔ Ｍｕｌｔｉｄｉｍｅｎｓｉｏｎａｌ Ｓａｍｐｌｉｎｇ」、Ｃｏｍｐｕｔｅｒ Ｇｒａｐｈｉｃｓ Ｆｏｒｕｍ、２１（３）：５５７〜５６３（２００２年９月）において提供されており、同文献は、その全体があたかも説明されているかのように、参照により本明細書に組み込まれる。
シミュレーション自体は、時間ステップｎ＝０で開始し、μの実現のため、ｘｌ，２を使用して、０≦ｌ≦Ｎについて、状態
Figure 2010501100

を初期化する。その後、当該技法は、状態をソートすること（以下のセクション３において詳細に説明される）によって続行され、Ｐに従った遷移を実現するため、０≦ｌ≦Ｎについて、
Figure 2010501100

を使用して、
Figure 2010501100

を計算することによって、連鎖を継続させる。遷移の次の状態を選択するためのインデックス
Figure 2010501100

は、実際には、基数２のｖａｎ ｄｅｒ Ｃｏｒｐｕｔ列Φ２を使用し、これは、（０，１）−列であり、したがって、（０，ｍ，１）−ネットの系列である。例えば、ｍ＝３＞ｔ＝０を選択すると、Ｎ＝２３＝８となり、
Figure 2010501100

について、
Figure 2010501100

となる。したがって、異なる時間ステップｎの間に使用されるすべてのインデックスは、実際には、すべてのｍについて同一である。
一様確率
Figure 2010501100

を仮定すると、収束定理（ｃｏｎｖｅｒｇｅｎｃｅ ｔｈｅｏｒｅｍ）が依然として適用されるが、より重要な安定なソーティングは、状態の順序を変化させない。したがって、実際には、インデックス置換は、収束条件に触れることなく、恒等置換（ｉｄｅｎｔｉｔｙ）として選択できることになる。同じことが、初期状態Ｘ０，ｌを選択するために適用され、それが、（０，２）−列 ｘｌの第２の要素だけを使用する、以下の簡略化された、しかし等価の技法をもたらす。
・ ｎ：＝０
・ ０≦ｌ＜２ｍについて、ｘｌ，２を使用して、Ｘ０，ｌを初期化する。
・ ループ
− 適切な順序を使用して、状態ベクトルをソートする。
− ｎ：＝ｎ＋１
− ０≦ｌ＜２ｍについて、ΦＳ（ｌ＋ｎ・２ｍ）を使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる。
２．１．２ いつ、なぜ、それが機能するか
多くの応用例において観測された、当該方式の改善された収束は、Ｐに従ったＸｎ，ｌの遷移を実現するために使用されるサンプル
ΦＳ（ｌ＋ｎ・２ｍ）
の構造によって引き起こされることが分かる。これは、暗黙の階層構造を露呈させる、根基逆関数（ｒａｄｉｃａｌ ｉｎｖｅｒｓｅ）の分解
Figure 2010501100

によって理解することができ、ΦＳ（ｌ）は、状態数ｌに依存する間隔１／（２ｍ）を有するオフセットであり、一方、シフトΦＳ（ｎ）は、時間ステップｎにおけるすべての区間について同一である。図１１は、シフトΦＳ（ｎ）を図説する２次元配列７００を示している。
ここで、ｌ＝０，．．．，２ｍ−１についてのΦＳ（ｌ）の全体は、実際には、（０，ｍ，１）−ネットでなければならず、したがって、サンプルの等距離集合でなければならないので、低次元設定は、サンプルがシフトされた格子又は階層化されたサンプルであるという、誤解を招く解釈をもたらすことがある。しかし、良好な性能は、ΦＳ（ｌ）が（ｔ，ｓ）−列であり、したがって、ｔ≦ｍ’≦ｍである任意のｍ’についての（ｔ，ｍ’，ｓ）−ネットの系列であるという特性から生じる。これは、状態空間において同様であり、したがって、ソーティング後の順序によって連なるｂｍ’個の状態が、（ｔ，ｍ’，ｓ）−ネットによって遷移をサンプリングすることを意味し、それが、良好な離散密度近似を保証する。最大の改善は、２ｍ個の連鎖すべてが同じ状態にある場合に獲得される。連鎖の状態が状態空間において離れていればいるだけ、性能改善は小さくなる。
２．２ ｓ次元での簡略化技法
（ｔ，ｍ，ｓ）−ネットの系列である、基底ｂの（ｔ，ｓ）−列を使用して、当該方式は、ｓ次元でも機能する。ソーティング後に状態が同様になるマルコフ連鎖は、低食い違い量サンプルによる遷移確率のサンプリングを保証される。ｓ次元での簡略化技法は、今では以下のような様相を呈する。
・ ｎ：＝０
・ 準モンテカルロ点ｘｌを使用して、Ｘ０，ｌを初期化する。
・ ループ
− 適切な順序を使用して、状態ベクトルをソートする。
− ｎ：＝ｎ＋１
− （ｔ，ｓ）−列からの後続のサンプルｘｌを使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる。
いくつかのシミュレーションは、ある局所的な微妙な効果を獲得するために、軌道分割を必要とする。これは、より複雑な手法を使用して、すでに対処されていたが、実際には、分割される状態について、単に（ｔ，ｓ）−列からより多くのサンプルを取り出すことによって、より簡単な方法で達成することができる。
これは、正確にｂｍ個の連鎖を同時にシミュレートする必要すらないという事実の結果である。最大性能利得を可能にするためには、（ｔ，ｓ）−列から後続のサンプルを取り出すこと、及び後続の（ｔ，ｍ，ｓ）−ネット内の点の数ｂｍを最小化することだけが重要である。しかし、（ｔ，ｓ）−列の選択は、（０，ｓ）−列はｂ≧ｓの場合にのみ存在し、ｍ＞ｔであるという条件によって制限される。
Ｈａｌｔｏｎ列又はそのスクランブルされた変形などの他の根基逆変換（ｒａｄｉｃａｌ ｉｎｖｅｒｓｉｏｎ）ベースの点集合は、（ｔ，ｓ）−列と同様の特性を満たし、同様の性能利得をもたらすことに留意されたい。図１２は、Ｈａｌｔｏｎ列、点ｘｌ＝（φ２（ｌ），φ３（ｌ））を図説するグラフ７２０を示している。示されたＨａｌｔｏｎ列では、根基逆変換は、最大の「穴」に落ち込む後続の点を有する、階層構造を暗示している。この特性は、任意の開始点にとって好適であり、（ｔ，ｓ）−列も、同じ特性を有する。
２．３ ランダム化
より高次元における決定論的技法の収束についての一般的証明は現在存在していないが、当該技法は、各時間ステップｎにおいて準モンテカルロ点を新たにランダム化することによって不偏になる。実際には、これはパディングされた複製のサンプリングの場合であるので、不偏性についての議論はより簡単になる。しかし、ランダム化は、特別な注意を払うに値する。
最も効率的な実施は、基底ｂ＝２の（ｔ，ｓ）−列を選択し、それから後続のサンプルを取り出し、そのサンプルにｓ次元ランダムベクトルによるＸＯＲ演算を施すことから成る。このランダムベクトルは、各遷移ステップの後、新たに取り出される。しかし、ランダムなスクランブリングは、点が列挙される順序を変化させるので、（ｔ，ｍ，ｓ）−ネットの系列の局所的特性も、変化させられる。
この結果は、これまでの手法を使用した場合に見られる効果のいくつかについての説明として利用することができる。配列（Ｒ）ＱＭＣシミュレーションで使用されるＳｏｂｏｌ点及びＫｏｒｏｂｏｖ点は、分散低減の１０倍までは、その変換された等価物（Ｓｏｂｏｌに対するグレイコード、Ｋｏｒｏｂｏｖに対するＢａｋｅｒ変換）よりも悪い。これについての説明は、点の構造において見出される。Ｓｏｂｏｌ列から抽出される（ｔ，ｍ，ｓ）−ネットの系列は、局所的には、そのグレイコード変形よりも悪い。同じことが、Ｋｏｒｏｂｏｖ格子とそれを変換した変形についても言える。
２．４ 配列ランダム化準モンテカルロ
以下は、ｓ次元問題のためのランダム化技法である。
ｎ：＝０、点列ｘｌを選択する。
ランダム化をインスタンス化する。
０≦ｌ＜Ｎについて、ｒａｎｄ（ｘｌ）を使用して、Ｘ０，ｌを初期化する。
ループ
・ 適切な順序を使用して、状態ベクトルをソートする。
・ ｎ：＝ｎ＋１
・ ランダム化をインスタンス化する。
・ ０≦ｌ＜Ｎについて、ｒａｎｄ（ｘｌ）を使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる。
当該技法は、多くの興味深い特性を有する。第１に、当該技法は、ｘｌの表にまとめられた１つの集合（ｏｎｅ ｔａｂｕｌａｔｅｄ ｓｅｔ）について、及びｘｌの連続ストリームについて不偏である。第２に、ストリームを使用する場合、分割及び吸収は本質的である。これは、ロシアンルーレット方式において真実であり、ｘｌからより多くの連続サンプルを取ることによって分割する場合にも真実である。
３．光輸送シミュレーション
数値的な証拠のため、上述の技法が、リアリスティックな画像を合成するための光輸送シミュレーションに適用される。基礎をなす積分方程式は、経路積分として定式化し直すことができる。
図１３は、双方向経路トレーシングによってサンプリング輸送経路空間７４０を示す図である。目７４２からシーン物体７４６までの軌道と光源７４４からシーン物体７４６までの軌道とは、マルコフ連鎖によって生成され、輸送される光の量を決定するために接続される。図１３に示されるサンプリング経路空間は、マルコフ連鎖をシミュレートすることに対応し、経路は、レイトレーシング及び散乱事象によって確立される。初期分布は、光源７４４の放出特性によって決定され、遷移確率は、衝突表面上における双方向反射分布関数によって与えられる。
経路積分を解くため、少なくとも２つの基本戦略が存在し、それらは、高次元低食い違い量点を使用するか、又は低次元低食い違い量点をパディングする。後者の手法は、上述の発見に適合しており、高次元事象は、マルコフ連鎖の後続する遷移として構成される。計測してみると、マルコフ連鎖シミュレーションのために高次元系列を使用することに対する差はないと言えるほど小さいが、パディングされた系列の使用は、計算的により高速であり、よりわずかな実施上の労力しか必要としない。それは、レンダリング業界の実務者にとってもより簡単である。
加えて、（ｔ，ｓ）−列又はＨａｌｔｏｎ列及びそのスクランブルされた変形の階層構造特性は、上で説明されたように、次元が小さい場合にはるかに良好になるので、低次元手法は、はるかに良好な結果を可能にする。
図１４Ａ〜Ｂ及び図１５Ａ〜Ｂは、モンテカルロ法を使用するシミュレーションと配列ランダム化準モンテカルロ（Ａ−ＲＱＭＣ）法を使用するシミュレーションの間の差を示した、レンダリングされた画像７６０、７７０、７８０、７９０である。モンテカルロ法を使用するシミュレーションは、図１４Ａ及び図１５Ａの画像７６０及び７８０に示されている。配列ランダム化準モンテカルロ法を使用するシミュレーションは、図１４Ｂ及び図１５Ｂの画像７７０及び７９０に示されている。各レンダリング画像の上部の矢印７６２、７７２、７８２、７９２は、放出位置及び方向を示している。図１４Ａ〜Ｂ及び図１５Ａ〜Ｂから、これらの例では、配列ランダム化ＱＭＣ法を使用してシミュレートされた光路は、モンテカルロ法を使用してシミュレートされた光路よりも、はるかに一様に分布していることが分かる。
３．１ Ｆｒｅｄｈｏｌｍ積分及びＶｏｌｔｅｒｒａ積分
光輸送の基礎をなす積分方程式は、Ｆｒｅｄｈｏｌｍ積分方程式又はＶｏｌｔｅｒｒａ積分方程式と見なすことができるが、そのことは実施にとって重要となる。
Figure 2010501100

は、点ｘにおいて方向ωの方に表面から反射された放射であり、積分演算子の定義域
Figure 2010501100

は、ｘにおける法線ｎ（ｘ）と軸平行な半球である。
図１６Ａ〜Ｄは、光源上の複数の開始点８０４、８１４、８２４、８３４から開始させた光子軌道８０２、８１２、８２２、８３２を示した、一連の図８００、８１０、８２０、８３０である。光線をトレースした後、表面８０６、８１６、８２６、８３６に衝突したとき、散乱の方向を決定して経路を継続するために、双方向反射分布関数がサンプリングされる。
ｆｒは、光表面特性（ｏｐｔｉｃａｌ ｓｕｒｆａｃｅ ｐｒｏｐｅｒｔｙ）を記述する双方向反射分布関数であり、Ｌは、入射放射である。この積分演算子の使用は、積分定義域がｘに依存するような、第２種のＶｏｌｔｅｒｒａ積分方程式
Figure 2010501100

をもたらし、他方で、ｘとは独立に単位球Ｓ２の全方向にわたって積分が実行されるような、第２種のＦｒｅｄｈｏｌｍ積分方程式をもたらす。
全方向を生成し、表面法線ｎ（ｘ）に対して負のスカラ積を有する方向を排除するという後者の手法の使用は、計算的に魅力があるが、第１の手法は、局所的な表面フレーム（ｓｕｒｆａｃｅ ｆｒａｍｅ）に変換されなければならない、半球における方向の生成を必要とし、これはよりコスト高になる。単位球からＳ２又は
Figure 2010501100

への写像は、当技術分野において知られている。
全方向を生成するためにはるかに重要な議論は、これから説明される手法に関連する。ソーティングによって、例えば角にある、異なる表面法線を有する２つのすぐ近くの表面点が、散乱方向を生成するために、（ｔ，ｓ）−列の後続のサンプルを使用することが起こり得る。全方向の生成は、今は良好に機能するが、後続のサンプルを使用する、２つの異なる局所的フレームにおける方向の生成は、低食い違い量性を打ち壊す。これらの不連続性は、はっきりと目に見えるようになる。しかし、全方向の使用は、ｆｒによるインポータンス（ｉｍｐｏｒｔａｎｃｅ）のサンプリング、又はコサイン項を可能にせず、それはしばしば不都合であり、さらなる研究に値する。
４．ソーティング技法及びシステム
上述の技法の効率を増大させるために、ソーティングを使用することが可能である。本明細書で説明されるように、（ｔ，ｓ）−列は、基底ｂの（ｔ，ｍ，ｓ）−ネットの系列であり、Ｈａｌｔｏｎ列及びその変形に類似している。類似の状態Ｘｎ，ｌは、連続する点を使用する。したがって、状態は、状態空間における近接性によって順序づけられるべきである。状態空間において状態がより近いほど、状態空間はより一様に探索される。
状態を互いに可能な限り近くするため、状態は、近隣の状態の距離の和が最小になるような順序で列挙されなければならない。これは、実際には、巡回セールスマン問題に関係するが、この問題では、与えられた都市の集合と、１つの都市から別の都市まで旅をするコストに対して、各都市を正確に１回訪れた後、出発した都市に戻る、最も安価な巡回旅行が求められる。この問題はすでに、１８世紀初頭にオイラーによって研究され、残念なことに、ＮＰ困難であることを示すことができる。
われわれの問題は、ルートの最後の状態から最初の状態に戻る必要がないことを除いて、非常に類似している。近似ではあるが最適に近い、巡回セールスマン問題の解を効率的に計算するために、いくつかの技法がすでに利用可能である。しかし、これらの技法の実行時間は、距離行列の計算だけでもＯ（Ｎ２）の計算量を示すので、シミュレーションでは許容可能ではないが、古典的なモンテカルロ法のＯ（Ｎ）の計算量に可能な限り近く、技法を維持することが望ましい。
以下では、高次元状態空間用の高速ソーティングを達成するための、本発明の態様による多くの順序が説明される。
４．１ 状態のノルム
クイックソートの平均計算量は、Ｏ（ＮｌｏｇＮ）であるが、あるシナリオの場合、基数ソート（ｒａｄｉｘ ｓｏｒｔ）など、Ｏ（Ｎ）の技法すら存在し、しかし、それは追加の一時メモリを必要とする。これらの１次元ソーティング技法を使用するため、多次元状態は、１次元まで引き下げられなければならない。多くの選択肢の中で、しばしば、何らかのノルム
Figure 2010501100

が、状態空間上で順序を定義するために使用される。しかし、類似のノルムが、必ずしも状態空間における近接性を示すわけではない。これについての簡単な例は、輸送シミュレーションにおける、空間内で遠く離れて配置された粒子の類似のエネルギーである。
４．２ 空間階層
多次元ソーティングを可能にする第２の可能性は、状態上で順序を定義するための、空間階層の使用である。この目的のために効率的なデータ構造は、ＢＳＰツリー、その特殊な軸平行サブセット、ｋＤツリー、又は境界ボリューム階層である。そのようなバイナリ階層の構築は簡単である。空間は、何らかの発見的方法によって選択された平面を使用して、再帰的に細分される。構築は、平均でＯ（ＮｌｏｇＮ）で実行される。階層の順序正しい横断が、近接性の順序でリーフを列挙する。この横断は、ツリーが左平衡化され、その結果、配列内に保存できる場合、自明なものとなる。
例えばレイトレーシングを加速するために、空間階層がとにかく使用されなければならない場合、階層のための付加的な構築時間は存在しない。その場合、階層のリーフに結びつけられたリンクリストとして、粒子が保存される。
図１７Ａ〜Ｇ、図１８Ａ〜Ｂ、及び図１９Ａ〜Ｃは、空間階層のリーフへの状態の分類が、階層を順序正しく横断することによって、近接性の順序を定義することを示した、一連の図である。
図１７Ａは、光源８５０と、レンダリングされる物体８６０とを含む、例示的なシーン８４０を示している。物体８６０は、五頂点の星形をしている。図１７Ｂでは、物体８６０は、軸平行境界ボックス８７０によって境界を定められる。図１７Ｃ〜Ｇでは、境界ボックス８７０が、再帰的にブランチに細分される。各ブランチは、リーフで終了する。
図１７Ｃでは、境界ボックス８７０が、物体８６０の第１の頂点８６１を通過する分割平面８７１によって、第１世代の左ブランチＬと第１世代の右ブランチＲとに分割される。
図１７Ｄでは、左ブランチＬが、物体８６０の第２の頂点８６２を通過する第２の分割平面８７２によって分割され、第２世代の左ブランチＬＬと第２世代の右ブランチＬＲとを生じさせる。
図１７Ｅでは、第３の分割平面８７３が、物体８６０の第３の頂点８６３を通過し、ブランチＬＲを第３世代の左ブランチＬＲＬと右ブランチＬＲＲとに分割する。
図１７Ｆでは、第４及び第５の分割平面８７４及び８７５が、物体８６０の第４及び第５の頂点８６４及び８６５を通過する。図１７Ｆに示されるように、結果は、星形８６０の５つの点の各々に対して１つの、５つのリーフＬＲＬ、ＬＲＲ、ＲＲ、ＬＬＬ、ＬＬＲ、及びＲＬである。
図１７Ｇは、各リーフ内の個々の点の照明を示している。図１８Ａでは、各リーフは、その明るさを示すために陰影づけされている。図１８Ｂでは、リーフは、それぞれの明るさに従って、階層的に順序づけられる。「最も暗い」から「最も明るい」まで、リーフは、以下のように、ＬＲＬ、ＬＲＲ、ＬＬＲ、ＲＲ、ＲＬの順に並べられる。
残念ながら、この順序の品質は、シミュレーション用に使用される空間階層の品質によって強く決定され、そのことは、階層内のリーフの数が連鎖Ｎの数よりもはるかに少ない場合、これは複数の状態が同じリーフにマッピングされることを引き起こすので、特に問題となる。
４．３ バケットソーティング及び空間充填曲線
線形時間の計算量を保証するため、バケットソーティングが使用できる。セクション２．１で略述された簡単な技法のｓ次元への拡張では、多次元の状態が、状態の第１の次元によってバケットに分類され、その後、各バケットの状態が、第２の次元に従ってバケットに分類され、以降も同様に行われた。この手順は、うまく機能するが、状態空間内で近い状態がソーティング手順によって引き離され得るという問題を有する。加えて、各次元の階層構造が使用されなければならず、それが、同時にシミュレートされるマルコフ連鎖の数において、次元の呪い（ｃｕｒｓｅ ｏｆ ｄｉｍｅｎｓｉｏｎ）を引き起こす。
したがって、状態空間は、バケットとして役立つ、等しいボリュームセル即ち「ボクセル」になる。各状態のバケットは、ボクセルグリッドの解像度に従って、状態座標を切り捨てることによって見出される。これは、連続状態空間にも、離散状態空間にも適用されることに留意されたい。近接性によるボクセルの列挙は、状態に所望の順序をもたらし、ボクセルの数に対する線形時間で実行できる。
近接性によってボクセルを列挙する順序は、ペアノ曲線、ヒルベルト曲線、又はＨインデクシング（Ｈ−ｉｎｄｅｘｉｎｇ）などの空間充填曲線によって与えられる。これらの曲線は、すべてのボクセルが正確に１度だけ訪れられ、全体の経路長が相対的に短いことを保証する。われわれ自身のシミュレーションにおける場合である、大きな形状に関連する問題にとって、これは、巡回セールスマン問題に対する高速でメモリ効率のよい近似解を生成するための、数少ない可能性の１つにすらなり得る。しかし、これらの曲線は、むしろ評価するのにコストがかかり、効率的であるために表にまとめられる必要があり、又はより高い次元に対しては利用可能でない。
幸い、ルベーグ曲線又はモートン順序（Ｍｏｒｔｏｎ ｏｒｄｅｒ）としても知られるＺ曲線は、これらの問題を回避する。多次元空間においてバケットの整数座標を与えると、その１次元Ｚ曲線インデックスが、図１９Ａ〜Ｃに示されるように、座標値をビット単位にインタリーブすることによって、簡単に計算される。
図１９Ａ〜Ｃは、２×２（図１９Ａ）、４×４（図１９Ｂ）、及び１６×１６（図１９Ｃ）のバケットについて、２次元のＺ曲線８９１、８９３、８９５を示した、一連の図８９０、８９２、８９４を示している。左上を原点（０，０）とすると、×によって印を付けられた点は、整数座標（３，４）を有し、それは、２進法の（０１１，１００）２に対応する。そのバイナリＺ曲線インデックス１００１０１２は、バイナリ座標をビット単位にインタリーブすることによって計算される。
これは、任意の次元及び問題サイズについて、計算が容易で高速であり、追加的なメモリを必要としない。その結果は、全体的な文脈では、例えばヒルベルト曲線ほど良好でないことが分かった。しかし、平均的な場合の分割品質、及び平均的な場合／最悪の場合の対数インデックス範囲は、同程度に良好である。問題は、数値的実験用に使用された、図２０Ａに示されたＳｈｉｒｌｅｙの「Ｓｃｅｎｅ ６」のような対称性の高いシーンにおいて生じる。（ｘ，ｙ）平面に平行な壁上の状態は、非常に良好にソートされるが、（ｙ，ｚ）平面に平行な他の２つの壁上に配置された状態は、曲線によって交互に訪れられ、それが、いくつかのシナリオでは、相関アーチファクトをもたらし得る。
５．数値的な証拠
上述の説明に続いて、不偏誤差推定を得るために、Ｃｒａｎｌｅｙ−Ｐａｔｔｅｒｓｏｎ回転によってランダム化された、Ｆａｕｒｅによる置換を行った、Ｈａｌｔｏｎ列が使用された。ソーティングについては、Ｚ曲線順序が、特定のセッティングに対して最良の結果をもたらし、以下の実験のために使用された。ランダム化の省略は結果の精度に顕著な影響をもたないことが数値的に検証されたことにさらに留意されたい。数値的な実験では、マルコフ連鎖をシミュレートするための４つの手法が比較された。
ＭＣ：メルセンヌツイスタ（Ｍｅｒｓｅｎｎｅ Ｔｗｉｓｔｅｒ）によって生成された一様乱数が、古典的モンテカルロサンプリング用に使用された。
ＲＱＭＣ：Ｃｒａｎｌｅｙ−Ｐａｔｔｅｒｓｏｎ回転によってランダム化された、Ｆａｕｒｅによる置換を行った、高次元Ｈａｌｔｏｎ列が使用され、要素のペアが、２次元放出及び散乱事象をサンプリングするために使用された。
ｌｏ−ｄｉｍ ＲＱＭＣＳ：Ｃｒａｎｌｅｙ−Ｐａｔｔｅｒｓｏｎ回転によってランダム化された、２次元Ｈａｌｔｏｎ列が使用された。Ｚ曲線が、バケットソーティングされた状態を列挙するために使用された。
ｈｉ−ｄｉｍ ＲＱＭＣＳ：Ｃｒａｎｌｅｙ−Ｐａｔｔｅｒｓｏｎ回転によってランダム化された、Ｆａｕｒｅによる置換を行った、高次元Ｈａｌｔｏｎ列が使用された。Ｚ曲線が、バケットソーティングされた状態を列挙するために使用された。
第１の実験では、経路積分を計算するために、強い全体的な照明技法が使用された。図２０Ａは、相対的に簡単な光輸送設定である、テストシーン「Ｓｈｉｒｌｅｙ ６」を示し、図２１Ａは、より複雑な光輸送設定である、テストシーン「Ｉｎｖｉｓｉｂｌｅ Ｄａｔｅ」を示している。図２０Ｂ及び図２１Ｂは、マスタ解に対する平均ＲＭＳ誤差を示す、グラフ９０１及び９１１であり、図２０Ｃ及び図２１Ｃは、（画像全体にわたって平均化された）大域的分散を示す、グラフ９０２及び９１２であり、図２０Ｄ及び図２１Ｄは、ピクセルベースの分散を示す、グラフ９０３及び９１３である。数値は、マルコフ連鎖の数を変化させた、１０回の独立実行を平均することによって獲得された。測定された数値は、簡単なテストシーンについてのみ説得力がある。複雑なケースでさえ、独立実行の回数が少なすぎるため、モンテカルロサンプリングを上回る性能利得は測定できない。より多くの実験は、過大な実行時間のせいで可能でなかった。
しかし、視覚誤差は、図２２Ａ〜Ｃに見られるような、劇的に異なったストーリを語り、非常に困難なセッティングにおいてさえ、新しい技法の紛れもない優位性が明らかになる。図２２Ａ〜Ｃは、シミュレーション用に３００個の連鎖を使用した、テストシーン「Ｉｎｖｉｓｉｂｌｅ Ｄａｔｅ」についての視覚的比較を示している。図２２Ａは、モンテカルロを使用してレンダリングされたテストシーンを示し、図２２Ｂは、ランダム化準モンテカルロを使用してレンダリングされたテストシーンを示し、図２２Ｃは、ソーティングを伴うランダム化準モンテカルロを使用してレンダリングされたテストシーンを示している。
図２２Ａ〜Ｃでは、唯一の光源は、隣の部屋の天井に配置されているので、見えていない。光子のより良い分布のため、ランダム化準モンテカルロ（ＲＱＭＣ）（図２２Ｂ）は、減少した陰影アーチファクトによって分かるように、モンテカルロ（ＭＣ）（図２２Ａ）よりも視覚的に性能がよい。（バケットソートのために２５６３個のボクセルを使用した）ソーティングを伴うＲＱＭＣ（ＲＱＭＣＳ）（図２２Ｃ）は、より多くの光子が第２の部屋に入り込み、ドアの裏側さえも非常によく照明されるので、さらに優れている。
このケースは例外でなく、多くのテストケースについて観測することができる。標準的な誤差測定は視覚的品質にとって適切な誤差測定でなく、このことは、コンピュータグラフィックスにおける、知られてはいるが未解決の問題であることだけを強調する。
図２３Ａ〜Ｄは、非常に困難な光輸送問題についての測定を示しており、光子は光源から直接トレースされ、最終経路セグメントはカメラに接続された（双方向経路トレーシングの１つの技法）。図２３Ａは、説明のために床と天井が取り除かれた、迷宮テストシーン９４０の概略図を示している。カメラは下方の部屋の中に配置され、光源は接続トンネルの他の終端に配置される。図２３Ｂに示されたグラフ９５０は、「箱ひげ」図であり、各技法を使用した５０回の独立した実現について、１０４８５７６のシミュレートされた光経路の、カメラによって受け取られた総放射の平均量（×によって印される）を示している。図２３Ｃに示されたグラフ９６０は、図２３Ｂに示されたグラフの関心のある部分を拡大したものである。図２３Ｄに示された表９７０は、５０回の独立した実現にわたって平均された、カメラ及び各ピクセル（２５６×２５６のピクセル解像度）によって受け取られた総放射のＬ１及びＬ２ノルム（分散）を使用した、マスタ解からの偏差を最終的に表示している。この困難なセッティングにおいて、新しい方法（ＲＱＭＣＳ）は、紛れもなく優れている。
上記の測定とは反対に、ただ１つの同時にシミュレートされるマルコフ連鎖が検討される。今は、十分な回数の実験が、計算的に実現可能であり、新しい技法の優位性が、明瞭に見えるものとなった。
６．結論
上述のシステム及び技法は、従来技術にまさる多くの利点を有することが分かる。これらは、以下のもの、即ち、直観的には系列特性に焦点が絞られた技法の簡略化、改善された収束、線形の計算量、及び光輸送シミュレーションへの適用を含む。上述のシステム及び技法は、技法を積分として定式化することに一歩近づき、その点において、一貫性を示すことが相対的に容易である。
発明者らは、首尾よく、マルコフ連鎖を同時にシミュレートするための技法を簡略化し、いつ、なぜ、状態のソーティングが収束を改善できるかについての直観を提供した。加えて、技法はもはや次元の呪いによって制限されず、シミュレーションは時間につれて変化し得る遷移確率Ｐ≡Ｐｎをまさに使用できるので、斉次マルコフ連鎖（ｈｏｍｏｇｅｎｅｏｕｓ Ｍａｒｋｏｖ ｃｈａｉｎ）への制限は存在しない。したがって、本明細書で説明された技法は、非斉次マルコフ連鎖（ｉｎｈｏｍｏｇｅｎｅｏｕｓ Ｍａｒｋｏｖ ｃｈａｉｎ）をシミュレートするために、即ち、Ｐ≡Ｐｎが時間ステップｎに依存する場合に、使用することができる。
実験は、必ずしもすべての（ｔ，ｓ）−列又は根基逆変換ベースの点列が等しく良好ではないことを示した。これは、さらなる特性解析に値する。
本明細書で説明された技法は、ランク１格子の系列が適用できる場合、さらに簡単になる。しかし、これまでの構築は、改善された性能にとって必要とされる、（ｔ，ｓ）−列の特性を欠いている。本発明のさらなる一態様によれば、所望の結果を達成するために、適切なランク１格子の系列を構築することが可能なこともある。
発明者らは全方向を使用しているので、即ち、球の積（ｐｒｏｄｕｃｔ ｏｆ ｓｐｈｅｒｅｓ）にわたって積分するので、その方向で最近の研究との関連性を確立することも興味深い。
７．一般的技法
図２４〜図３０は、本発明の上述の態様による、多くの技法及び副技法のフローチャートである。個々のボックス及びボックスのグループを含む、以下の技法及び副技法の一部又は全部は、本明細書で説明された本発明を実施するために、望みどおりに組み合わされてよいことに留意されたい。
図２４は、本発明による、第１の一般的技法１０００のフローチャートである。技法１０００は以下を含む。
ボックス１００１：準モンテカルロ法を使用して、マルコフ連鎖をシミュレートする。
ボックス１００２：状態をソートする。
（近接性ソーティングを含む、状態によるソーティングを含むことができる）
ボックス１００３：複数のマルコフ連鎖を同時にシミュレートする。
（準モンテカルロ点を利用できる。Ｈａｌｔｏｎ列、Ｈａｌｔｏｎ列の変形、格子列、又は（ｔ，ｓ）−列のいずれかを使用して、準モンテカルロ点を選択できる。任意の数の連鎖を利用できる。軌道分割によるなど、軌道をオンザフライで追加するステップを含むことができる。ロシアンルーレット法又は粒子吸収をシミュレートする他の技法を利用できる）
図２５は、ｓ次元問題について複数のマルコフ連鎖を同時にシミュレートするための、本発明による、さらなる技法１１００のフローチャートである。技法１１００は以下を含む。
ボックス１１０１：ｎ：＝０
（ｎを０に等しく設定する）
ボックス１１０２：準モンテカルロ点ｘｌを使用して、Ｘ０，ｌを初期化する。
ボックス１１０３：適切な順序を使用して、状態ベクトルをソートする。
ボックス１１０４：ｎ：＝ｎ＋１
（ｎを１だけ増やす）
ボックス１１０５：次の準モンテカルロ点ｘｌを使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる。
（任意の選択されたランダム化を施された任意の選択された点集合又は点列、及び状態空間における近接性による状態の順序づけも利用できる）
ボックス１１０３〜１１０５は、ループ内で実行される。
図２６は、追加技法１２００のフローチャートであり、その１つ又は複数は、図２５に示された技法１１００と併せて実行されてよい。これらの技法１２００は以下を含む。
ボックス１２０１：明るさ、又は状態のノルムによってソートする。
ボックス１２０２：空間階層を使用してソートする。
（バイナリ階層、又はＢＳＰツリー、ｋＤツリー、若しくはＢＳＰツリー構造の他の軸平行サブセットのいずれか、又は境界ボリューム階層、又は規則的ボクセルを利用するステップを含むことができる。バイナリ階層は、選択された発見的手法によって選択された平面を使用して再帰的に細分を行い、その後、階層のリーフを近接性の順序で列挙するために、階層を順序正しく横断することによって構築することができる。階層の効率的な横断を可能にするために、ツリーを左平衡化させるステップと、ツリーを配列データ構造に保存するステップも含むことができる）
ボックス１２０３：少なくとも１つの空間充填曲線によって列挙されたバケットを使用してソートする。
図２７は、本発明による、さらなる技法１３００のフローチャートである。技法１３００は以下を含む。
ボックス１３０１：ｎ：＝０、点列ｘｌを選択する。
ボックス１３０２：ランダム化をインスタンス化する。
ボックス１３０３：ｒａｎｄ（ｘｌ）を使用して、Ｘ０，ｌを初期化する。
ボックス１３０４：何らかの順序を使用して、状態ベクトルをソートする。
ボックス１３０５：ｎ：＝ｎ＋１
ボックス１３０６：ランダム化をインスタンス化する。
ボックス１３０７：ｒａｎｄ（ｘｌ）を使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる。
ボックス１３０４〜１３０７は、ループ内で実行される。
図２８は、図２７の技法１３００と併せて実行されてよい、追加技法１４００のフローチャートである。追加技法１４００は、各時間ステップｎにおいて準モンテカルロ点をランダム化し、以下を含む。
ボックス１４０１：後続のサンプルがそこから取り出される、基数ｂ＝２の（ｔ，ｓ）−列を選択する。
ボックス１４０２：サンプルにｓ次元ランダムベクトルによるＸＯＲ演算を施し、ランダムベクトルは、各遷移ステップの後に生成される。
図２９は、本発明による、さらなる技法１５００のフローチャートである。技法１５００は以下を含む。
ボックス１５０１：軸平行境界ボックスによって、レンダリングされる物体の境界を定める。
ボックス１５０２：物体の選択された点を通過する分割平面を適用することによって、境界ボックスを、各々がリーフで終了する、連続的な左ブランチ及び右ブランチに再帰的又は反復的に分割し、それによって、複数のリーフを生成する。
ボックス１５０３：リーフをそれぞれの明るさに従って階層的に順序づける。
ボックス１５０４：選択された数のバケットに分割された、選択されたサイズの行列上で、バケットソートを実行する。
ボックス１５０５：行列内をたどるために選択された空間充填曲線を使用し、個々のバケット内において、より小さな空間充填曲線が、行列内の個々のセル内をたどるために使用される。
図３０は、本発明による、さらなる技法１６００のフローチャートである。技法１６００は、リアルに見える画像を合成するための光輸送シミュレーションを可能にするシミュレーションを構成し、以下を含む。
ボックス１６０１：光輸送シミュレーションの基礎をなす積分方程式を、経路空間をサンプリングすることがマルコフ連鎖をシミュレートすることに対応するように、経路積分として定式化し直し、その際、経路は、レイトレーシング及び散乱事象によって確立される。
ボックス１６０２：初期分布を、光源又はセンサの放出特性によって決定し、遷移確率を、画像内の表面についての双方向反射分布関数及び双方向表面下散乱分布関数によって決定する。
ボックス１６０３：透過効果を、経路積分、初期分布、及び遷移確率を利用することによって処理する。
ボックス１６０４：経路積分を、任意の種類のランダム化を施して、又は施さずに、高次元準モンテカルロ点を利用することによって、又は低次元準モンテカルロ点をパディングすることによって解く。
（光輸送シミュレーションの基礎をなす積分方程式を、Ｆｒｅｄｈｏｌｍ積分方程式又はＶｏｌｔｅｒｒａ積分方程式と見なす）
上記の説明は、当業者が本発明を実施することを可能にする詳細を含むが、説明が本質的に例示的なものであること、これらの教示から利益を得る当業者には本発明の多くの修正及び変形が明らかであることを理解されたい。したがって、本明細書の発明が、本明細書に添付された特許請求の範囲によってのみ確定されること、特許請求の範囲が、従来技術によって許される限り広く解釈されることが意図されている。
コンピュータプログラムリスト

CODE LISTING 2.2.1

void Triangle::Transform()
{
Point *p = (Point *)this;
Vector n3d;
Vector n_abs = n3d =(p[1]-p[0])|(p[2]-p[0]);
// search largest component forprojection (0=x,1=y,2=z)
uintCast(n_abs.dx) &= 0x7FFFFFFF;
uintCast(n_abs.dy) &= 0x7FFFFFFF;
uintCast(n_abs.dz) &= 0x7FFFFFFF;
// Degenerated Triangles must be handled(set edge-signs)
if(!((n_abs.dx+n_abs.dy+n_abs.dz) >DEGEN_TRI_EPSILON))
//(!(...) > EPS) to handle NaN’s
{
d = p[0].x;
p0.u = -p[0].y;
p0.v = -p[0].z;
n.u=n.v = 0.0f;
e[0].u = e[1].u = e[0].v = e[1].v = 1.0f;
return;
}
U32 axis = 2;
if(n_abs.dx > n_abs.dy)
{
if(n_abs.dx > n_abs.dz)
axis = 0;
}
else if(n_abs.dy > n_abs.dz)
axis = 1;
Point p03d = p[0];
Point p13d = p[1];
Point p23d = p[2];
float t_inv = 2.0f/n3d[axis];
e[0].u =(p23d[PlusOneMod3[axis]]-p03d[PlusOneMod3[axis]])*t_inv;
e[0].v =(p23d[PlusOneMod3[axis+1]]-p03d[PlusOneMod3[axis+1]])*t_inv;
e[1].u =(p13d[PlusOneMod3[axis]]-p03d[PlusOneMod3[axis]])*t_inv;
e[1].v =(p13d[PlusOneMod3[axis+1]]-p03d[PlusOneMod3[axis+1]])*t_inv;
t_inv *= 0.5f;
n.u = n3d[PlusOneMod3[axis]] *t_inv;
n.v = n3d[PlusOneMod3[axis+1]]*t_inv;
p0.u = -p03d[PlusOneMod3[axis]];
p0.v = -p03d[PlusOneMod3[axis+1]];
d = p03d[axis] +n.u*p03d[PlusOneMod3[axis]] + n.v*p03d[PlusOneMod3[axis+1]];
}


CODE LISTING 2.2.2

U32 *idx = pointer_to_face_indices;
U32 ofs = projection_case;
for(U32 ii = num_triData; ii;ii--,idx++)
{
float t = (triData[*idx].d -ray.from[ofs]
-triData[*idx].n.u*ray.from[PlusOneMod3[ofs]]
-triData[*idx].n.v*ray.from[PlusOneMod3[ofs+1]])
/ (ray.d[ofs] +triData[*idx].n.u*ray.d[PlusOneMod3[ofs]]
+ triData[*idx].n.v*ray.d[PlusOneMod3[ofs+1]]);
if(uintCast(t)-1 >uintCast(result.tfar)) //-1 for +0.0f
continue;
float h1 = t*ray.d[PlusOneMod3[ofs]] +ray.from[PlusOneMod3[ofs]]
+ triData[*idx].p0.u;
float h2 = t*ray.d[PlusOneMod3[ofs+1]] +ray.from[PlusOneMod3[ofs+1]]
+ triData[*idx].p0.v;
float u = h1*triData[*idx].e[0].v -h2*triData[*idx].e[0].u;
float v = h2*triData[*idx].e[1].u -h1*triData[*idx].e[1].v;
float uv = u+v;
if((uintCast(u) | uintCast(v) |uintCast(uv)) > 0x40000000)
continue;
result.tfar = t;
result.tri_index = *idx;
}


CODE LISTING 2.3.1

Point *p = (Point*)&triData[tri_index];
int boxMinIdx, boxMaxIdx;
// boxMinIdx and boxMaxIdx index thesmallest and largest vertex of the triangle
// in the component dir[0] of the splitplane
if(p[0][dir[0]] < p[1][dir[0]])
{
if(p[2][dir[0]] < p[0][dir[0]])
{
boxMinIdx = 2;
boxMaxIdx = 1;
}
else
{
boxMinIdx = 0;
boxMaxIdx = p[2][dir[0]] <p[1][dir[0]] ? 1 : 2;
}
}
else
{
if(p[2][dir[0]] < p[1][dir[0]])
{
boxMinIdx = 2;
boxMaxIdx = 0;
}
else
{
boxMinIdx = 1;
boxMaxIdx = p[2][dir[0]] <p[0][dir[0]] ? 0 : 2;
}
}
/* If the triangle is in the split planeor completely on one side of the split plane
is decided without any numerical errors,i.e. at the precision the triangle is
entered to the rendering system. Usingepsilons here is wrong and not necessary.
*/
if((p[boxMinIdx][dir[0]] == split)&& (p[boxMaxIdx][dir[0]] == split)) // in split plane ?
{
on_splitItems++;
if(split < middle_split) // put tosmaller volume
left_num_divItems++;
else
{
unsorted_border--;
U32 t = itemsList[unsorted_border];
right_num_divItems--;
itemsList[right_num_divItems] =itemsList[left_num_divItems];
itemsList[left_num_divItems] = t;
}
}
else if(p[boxMaxIdx][dir[0]] <=split) // triangle completely left ?
left_num_divItems++;
else if(p[boxMinIdx][dir[0]] >=split) // triangle completely right ?
{
unsorted_border--;
U32 t = itemsList[unsorted_border];
right_num_divItems--;
itemsList[right_num_divItems] =itemsList[left_num_divItems];
itemsList[left_num_divItems] = t;
}
else
// and now detailed decision, trianglemust intersect split plane ...
{
/* In the sequel we determine whether atriangle should go left and/or right, where we already know that it mustintersect the split plane in a line segment.
All computations are ordered so that themore precise computations are done first. Scalar products and cross productsare evaluated last.
In some situations it may be necessaryto expand the bounding box by
an epsilon. This, however, will blow upthe required memory by large amounts.
If such a situation is encountered, itmay be better to analyze it numerically in order not to use any epsilons...
Arriving here we know thatp[boxMaxIdx][dir[0]] < split < p[boxMaxIdx][dir[0]]
and that p[boxMidIdx][dir[0]] \in[p[boxMaxIdx][dir[0]], p[boxMaxIdx][dir[0]]].
We also know, that the triangle has anon-empty intersection with the current
voxel. The triangle also cannot lie inthe split plane, and its vertices cannot
lie on one side only.
*/
int boxMidIdx = 3 - boxMaxIdx -boxMinIdx; // missing index, found by 3 = 0 + 1 + 2
/* We now determine the vertex that isalone on one side of the split plane.
Depending on whether the lonely vertexis on the left or right side,
we have to later swap the decision,whether the
triangle should be going to the left orright.
*/
int Alone = (split <p[boxMidIdx][dir[0]]) ? boxMinIdx : boxMaxIdx;
int NotAlone = 3 - Alone - boxMidIdx;
// == (split < p[boxMidIdx][dir[0]])? boxMaxIdx : boxMinIdx;
// since sum of idx = 3 = 0 + 1 + 2
float dist = split - p[Alone][dir[0]];
U32 swapLR = uintCast(dist)>>31;// == p[Alone][dir[0]] > split;
/* Now the line segments connecting thelonely vertex with the remaining two verteces
are intersected with the split plane. a1and a2 are the intersection points.
The special case "if(p[boxMidIdx][dir[0]]== split)" [yields a x / x, which could
be optimized] does not help at all sinceit only can happen as often as the highest
valence of a vertex of the mesh is...
*/
float at = dist / (p[boxMidIdx][dir[0]]- p[Alone][dir[0]]);
float at2 = dist / (p[NotAlone][dir[0]]- p[Alone][dir[0]]);
float a1x = (p[boxMidIdx][dir[1]] -p[Alone][dir[1]]) * at;
float a1y = (p[boxMidIdx][dir[2]] -p[Alone][dir[2]]) * at;
float a2x = (p[NotAlone][dir[1]] -p[Alone][dir[1]]) * at2;
float a2y = (p[NotAlone][dir[2]] - p[Alone][dir[2]])* at2;
// n is a vector normal to the line ofintersection a1a2 of the triangle
// and the split plane
float nx = a2y - a1y;
float ny = a2x - a1x;
// The signs indicate the quadrant ofthe vector normal to the intersection line
U32 nxs = uintCast(nx)>>31; // == (nx< 0.0f)
U32 nys = uintCast(ny)>>31; // ==(ny < 0.0f)
/* Numerical precision: Due tocancellation, floats of approximately same exponent
should be subtracted first, beforeadding something of a different order of
magnitude. All brackets in the sequelare ESSENTIAL for numerical precision.
Change them and you will see more errorsin the BSP...
pMin is the lonely point in thecoordinate system with the origin at
bBox.bMinMax[0]
pMax is the lonely point in thecoordinate system with the origin at
bBox.bMinMax[1]
*/
float pMinx = p[Alone][dir[1]] -bBox.bMinMax[0][dir[1]];
float pMiny = p[Alone][dir[2]] -bBox.bMinMax[0][dir[2]];
float pMaxx = p[Alone][dir[1]] -bBox.bMinMax[1][dir[1]];
float pMaxy = p[Alone][dir[2]] -bBox.bMinMax[1][dir[2]];
// Determine coordinates of the boundingbox, however, with respect to p + a1 being the origin.
float boxx[2];
float boxy[2];
boxx[0] = (pMinx + a1x) * nx;
boxy[0] = (pMiny + a1y) * ny;
boxx[1] = (pMaxx + a1x) * nx;
boxy[1] = (pMaxy + a1y) * ny;
/* Test, whether line of intersection ofthe triangle and the split plane passes by the
bounding box. This is done by indexingthe coordinates of the bounding box by the
quadrant of the vector normal to theline of intersection. In fact this is
the nifty implementation of the 3d testintroduced by in the book with Haines:
"Real-Time Rendering"
By the indexing the vertices areselected, which are farthest from the line.
Note that the triangle CANNOT completelypass the current voxel, since it must have
a nonempty intersection with it.
*/
U32 resultS;
if(pMinx + MAX(a1x,a2x) < 0.0f) //line segment of intersection a1a2 left of box
resultS = uintCast(pMinx)>>31;
else if(pMiny + MAX(a1y,a2y) < 0.0f)// line segment of intersection a1a2 below box
resultS = uintCast(pMiny)>>31;
else if(pMaxx + MIN(a1x,a2x) > 0.0f)// line segment of intersection a1a2 right of box
resultS = (pMaxx > 0.0f);
else if(pMaxy + MIN(a1y,a2y) > 0.0f)// line segment of intersection a1a2 above box
resultS = (pMaxy > 0.0f);
else if(boxx[1^nxs] > boxy[nys])
// line passes beyond bbox ? =>triangle can only be on one side
resultS = (a1y*a2x > a1x*a2y);
// sign of cross product a1 x a2 ischecked to determine side
else if(boxx[nxs] < boxy[1^nys])
resultS = (a1y*a2x < a1x*a2y);
else
// Ok, now the triangle must be bothleft and right
{
stackList[currStackItems++] =itemsList[left_num_divItems];
unsorted_border--;
itemsList[left_num_divItems] =itemsList[unsorted_border];
continue;
}
if(swapLR != /*^*/ resultS)
{
unsorted_border--;
U32 t = itemsList[unsorted_border];
right_num_divItems--;
itemsList[right_num_divItems] =itemsList[left_num_divItems];
itemsList[left_num_divItems] = t;
}
else
left_num_divItems++;
}


CODE LISTING 2.4.1

Intersection Boundary::Intersect(Ray&ray) //ray.tfar is changed!
{
// Optimized inverse calculation (saves2 of 3 divisions)
float inv_tmp =(ray.d.dx*ray.d.dy)*ray.d.dz;
if((uintCast(inv_tmp)&0x7FFFFFFF)> uintCast(DIV_EPSILON))
{
inv_tmp = 1.0f / inv_tmp;
ray.inv_d.dx =(ray.d.dy*ray.d.dz)*inv_tmp;
ray.inv_d.dy = (ray.d.dx*ray.d.dz)*inv_tmp;
ray.inv_d.dz =(ray.d.dx*ray.d.dy)*inv_tmp;
}
else
{
ray.inv_d.dx =((uintCast(ray.d.dx)&0x7FFFFFFF) > uintCast(DIV_EPSILON)) ?
(1.0f / ray.d.dx) :INVDIR_LUT[uintCast(ray.d.dx) >> 31];
ray.inv_d.dy =((uintCast(ray.d.dy)&0x7FFFFFFF) > uintCast(DIV_EPSILON)) ?
(1.0f / ray.d.dy) :INVDIR_LUT[uintCast(ray.d.dy) >> 31];
ray.inv_d.dz =((uintCast(ray.d.dz)&0x7FFFFFFF) > uintCast(DIV_EPSILON)) ?
(1.0f / ray.d.dz) :INVDIR_LUT[uintCast(ray.d.dz) >> 31];
}
Intersection result;
result.tfar = ray.tfar;
result.tri_index = -1;
//
//BBox-Check
//
float tnear = 0.0f;
worldBBox.Clip(ray,tnear);
if(uintCast(ray.tfar) == 0x7F7843B0)//ray.tfar==3.3e38f //!!
return(result);
//
U32 current_bspStack = 1; //wegen emptystack case == 0
U32 node = 0;
//
//BSP-Traversal
//
const U32 whatnode[3] ={(uintCast(ray.inv_d.dx)>>27) & sizeof(BSPNODELEAF),
(uintCast(ray.inv_d.dy)>>27) &sizeof(BSPNODELEAF),
(uintCast(ray.inv_d.dz)>>27) &sizeof(BSPNODELEAF)};
U32 bspStackNode[128];
float bspStackFar[128];
float bspStackNear[128];
bspStackNear[0] = -3.4e38f; // sentinel
do
{
//Ist Node ein Leaf (type<0) oder nur neVerzweigung (type>=0)
while(((BSPNODELEAF&)bspNodes[node]).type >= 0)
{
//Split-Dimension (x|y|z)
U32 proj =((BSPNODELEAF&)bspNodes[node]).type & 3;
float distl =(((BSPNODELEAF&)bspNodes[node]).splitlr[whatnode[proj]>>4]
- ray.from[proj])*ray.inv_d[proj];
float distr =(((BSPNODELEAF&)bspNodes[node]).splitlr[(whatnode[proj]>>4)^1]
- ray.from[proj])*ray.inv_d[proj];
node =(((BSPNODELEAF&)bspNodes[node]).type - proj) | whatnode[proj];
//type & 0xFFFFFFF0
if(tnear <= distl)
{
if(ray.tfar >= distr)
{
bspStackNear[current_bspStack] =MAX(tnear,distr);
bspStackNode[current_bspStack] =node^sizeof(BSPNODELEAF);
bspStackFar[current_bspStack] =ray.tfar;
current_bspStack++;
}
ray.tfar = MIN(ray.tfar,distl);
}
else
if(ray.tfar >= distr)
{
tnear = MAX(tnear,distr);
node ^= sizeof(BSPNODELEAF);

}
else
goto stackPop;
}
//
//Faces-Intersect
... code omitted ...
//
//
//Hit gefunden?
//
do //!! NEEDS bspStackNear[0] = -3.4e38f;!!!!!
{
stackPop:
current_bspStack--;
tnear = bspStackNear[current_bspStack];
}while(result.tfar < tnear);
if(current_bspStack == 0)
return(result);
node = bspStackNode[current_bspStack];
ray.tfar =bspStackFar[current_bspStack];
} while (true);
}

Claims (37)

  1. 画像内のピクセルのピクセル値を生成するためのコンピュータグラフィックスシステムであり、
    前記ピクセル値が、シーン内の点を、シミュレートされるカメラの画像平面上に記録されるように表し、
    前記コンピュータグラフィックスシステムが、選択された方向に沿って前記ピクセルからシーン内に向かう少なくとも１つの光線をシミュレートするステップを含む選択されたレイトレーシング法を使用して、画像の前記ピクセル値を生成するように構成され、
    前記レイトレーシング法が、光線と前記シーン内の物体の表面との交わりを計算するステップと、前記シーン内の物体を照明する光線の軌道をシミュレートするステップとをさらに含み、軌道をシミュレートする前記ステップが、マルコフ連鎖を使用する、前記コンピュータグラフィックスシステムの改良であって、
    準モンテカルロ法を使用して前記マルコフ連鎖をシミュレートするステップを含み、
    マルコフ連鎖をシミュレートする前記ステップが、状態をソートするステップを含み、
    ソートする前記ステップが、近接性ソーティングを含む、コンピュータグラフィックスシステムの改良。
  2. 複数のマルコフ連鎖を同時にシミュレートするステップを含む、請求項１に記載のコンピュータグラフィックスシステムの改良。
  3. 複数のマルコフ連鎖を同時にシミュレートする前記ステップが、準モンテカルロ点を利用するステップを含む、請求項２に記載のコンピュータグラフィックスシステムの改良。
  4. 前記準モンテカルロ点が、Ｈａｌｔｏｎ列、Ｈａｌｔｏｎ列の変形、格子列、又は（ｔ，ｓ）−列のいずれかを使用して選択される、請求項３に記載のコンピュータグラフィックスシステムの改良。
  5. 複数のマルコフ連鎖を同時にシミュレートする前記ステップが、任意の数の連鎖を利用する、請求項２に記載のコンピュータグラフィックスシステムの改良。
  6. 複数のマルコフ連鎖を同時にシミュレートする前記ステップが、追加的な軌道をオンザフライで追加するステップをさらに含む、請求項２に記載のコンピュータグラフィックスシステムの改良。
  7. 複数のマルコフ連鎖を同時にシミュレートする前記ステップが、軌道分割をさらに含む、請求項２に記載のコンピュータグラフィックスシステムの改良。
  8. 複数のマルコフ連鎖を同時にシミュレートする前記ステップが、粒子吸収をシミュレートする技法を利用するステップを含む、請求項２に記載のコンピュータグラフィックスシステムの改良。
  9. 複数のマルコフ連鎖を同時にシミュレートする前記ステップが、ロシアンルーレット法を使用するステップをさらに含む、請求項８に記載のコンピュータグラフィックスシステムの改良。
  10. ｓ次元問題について、複数のマルコフ連鎖を同時にシミュレートする前記ステップが、
    ｎ：＝０
    準モンテカルロ点ｘｌを使用して、Ｘ０，ｌを初期化する
    ループ：
    − 適切な順序を使用して、状態ベクトルをソートする
    − ｎ：＝ｎ＋１
    − 次の準モンテカルロ点ｘｌを使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる
    というステップを含む、請求項２に記載のコンピュータグラフィックスシステムの改良。
  11. 状態空間における近接性によって状態を順序づけるステップを含む、請求項１０に記載のコンピュータグラフィックスシステムの改良。
  12. ランダム化を含む、請求項１１に記載のコンピュータグラフィックスシステムの改良。
  13. 任意の選択されたランダム化を施された任意の選択された点集合又は点列を利用するステップを含む、請求項１２に記載のコンピュータグラフィックスシステムの改良。
  14. ランダム化と、以下の技法、即ち、
    ｎ：＝０、点列ｘｉを選択する
    ランダム化をインスタンス化する
    ｒａｎｄ（ｘｌ）を使用して、Ｘ０，ｌを初期化する
    ループ
    ・ 何らかの順序を使用して、状態ベクトルをソートする
    ・ ｎ：＝ｎ＋１
    ・ ランダム化をインスタンス化する
    ・ ｒａｎｄ（ｘｌ）を使用して、Ｘｎ，ｌを計算することによって、連鎖を継続させる
    を適用するステップとを含む、請求項１１に記載のコンピュータグラフィックスシステムの改良。
  15. 各時間ステップｎにおいて前記準モンテカルロ点をランダム化するステップを含み、
    ランダム化する前記ステップが、
    後続のサンプルがそこから取り出される、基数ｂ＝２の（ｔ，ｓ）−列を選択するステップと、
    前記サンプルにｓ次元ランダムベクトルによるＸＯＲ演算を施すステップであって、前記ランダムベクトルが、各遷移ステップの後に生成される、ステップと
    を含む、請求項１１に記載のコンピュータグラフィックスシステムの改良。
  16. 選択された数の追加的な軌道をオンザフライで追加するステップを含む、請求項１５に記載のコンピュータグラフィックスシステムの改良。
  17. 軌道分割を含む、請求項１６に記載のコンピュータグラフィックスシステムの改良。
  18. 粒子吸収をシミュレートする技法を利用するステップを含む、請求項１７に記載のコンピュータグラフィックスシステムの改良。
  19. ロシアンルーレット法を使用するステップを含む、請求項１８に記載のコンピュータグラフィックスシステムの改良。
  20. 明るさ、又は状態のノルムによってソートするステップと、
    空間階層を使用してソートするステップと、
    少なくとも１つの空間充填曲線によって列挙されたバケットを使用してソートするステップと
    を含む、請求項２に記載のコンピュータグラフィックスシステムの改良。
  21. 空間階層を使用してソートするステップが、バイナリ階層を利用するステップを含む、請求項２０に記載のコンピュータグラフィックスシステムの改良。
  22. 前記空間階層が、ＢＳＰツリー、ｋＤツリー、若しくはＢＳＰツリー構造の他の軸平行サブセットのいずれか、又は境界ボリューム階層、又は規則的ボクセルを含む、請求項２１に記載のコンピュータグラフィックスシステムの改良。
  23. バイナリ階層が、選択された発見的手法によって選択された平面を使用して再帰的に細分を行い、
    その後、前記階層のリーフを近接性の順序で列挙するために、前記階層を順序正しく横断することによって、構築される、請求項２０に記載のコンピュータグラフィックスシステムの改良。
  24. 前記ツリーを左平衡化させ、前記ツリーを配列データ構造に保存して、前記階層の効率的な横断を可能にするステップを含む、請求項２３に記載のコンピュータグラフィックスシステムの改良。
  25. 前記空間階層がリーフを含み、
    前記空間階層を使用してソートするステップが、バケットソーティング及び選択された空間充填曲線を利用するステップをさらに含む、請求項２４に記載のコンピュータグラフィックスシステムの改良。
  26. 前記空間階層を使用してソートするステップが、規則的ボクセルを使用するステップをさらに含む、請求項２４に記載のコンピュータグラフィックスシステムの改良。
  27. 軸平行境界ボックスによって、レンダリングされる物体の境界を定めるステップと、
    前記物体の選択された点を通過する分割平面を適用することによって、前記境界ボックスを、各々がリーフで終了する、連続的な左ブランチ及び右ブランチに再帰的又は反復的に分割し、それによって、複数のリーフを生成するステップと、
    前記リーフをそれぞれの明るさに従って階層的に順序づけるステップと、
    選択された数のバケットに分割された、選択されたサイズの行列上で、バケットソートを実行するステップと、
    前記行列内をたどるために選択された空間充填曲線を使用するステップであって、個々のバケット内において、より小さな空間充填曲線が、前記行列内の個々のセル内をたどるために使用されるステップと
    を含む、請求項２５に記載のコンピュータグラフィックスシステムの改良。
  28. 前記シミュレーションが積分方程式を解くための構成用に動作可能である、請求項２に記載のコンピュータグラフィックスシステムの改良。
  29. 前記シミュレーションが、リアルに見える画像を合成するための光移送シミュレーションを可能にするように動作可能である、請求項２８に記載のコンピュータグラフィックスシステムの改良。
  30. 光移送シミュレーションの基礎をなす積分方程式が、経路空間をサンプリングすることがマルコフ連鎖をシミュレートすることに対応するように、経路積分として定式化し直され、その際、経路が、レイトレーシング及び散乱事象によって確立され、
    初期分布が、光源又はセンサの放出特性によって決定され、遷移確率が、前記画像内の表面についての双方向反射分布関数及び双方向表面下散乱分布関数によって決定される、請求項２９に記載のコンピュータグラフィックスシステムの改良。
  31. 前記経路積分、初期分布、及び遷移確率を利用することによって、透過効果を処理するステップを含む、請求項３０に記載のコンピュータグラフィックスシステムの改良。
  32. 前記経路積分が、高次元準モンテカルロ点を利用することによって、又は低次元準モンテカルロ点をパディングすることによって、解かれる、請求項３０に記載のコンピュータグラフィックスシステムの改良。
  33. 光移送シミュレーションの基礎をなす前記積分方程式が、Ｆｒｅｄｈｏｌｍ積分方程式又はＶｏｌｔｅｒｒａ積分方程式と見なされる、請求項３２に記載のコンピュータグラフィックスシステムの改良。
  34. 画像内のピクセルのピクセル値を生成するためのコンピュータグラフィックスシステムであり、
    前記ピクセル値が、シーン内の点を、シミュレートされるカメラの画像平面上に記録されるように表し、
    前記コンピュータグラフィックスシステムが、選択された方向に沿って前記ピクセルからシーン内に向かう少なくとも１つの光線をシミュレートするステップを含む選択されたレイトレーシング法を使用して、画像の前記ピクセル値を生成するように構成され、
    前記レイトレーシング法が、光線と前記シーン内の物体の表面との交わりを計算するステップと、前記シーン内の物体を照明する光線の軌道をシミュレートするステップとをさらに含み、
    軌道をシミュレートする前記ステップが、マルコフ連鎖を使用する、前記コンピュータグラフィックスシステムにおける方法であって、
    準モンテカルロ法を使用して前記マルコフ連鎖をシミュレートするステップを含み、
    マルコフ連鎖をシミュレートする前記ステップが、状態をソートするステップを含み、
    ソートする前記ステップが、近接性ソーティングを含む、方法。
  35. 画像内のピクセルのピクセル値を生成するためのコンピュータグラフィックスシステムであり、
    前記ピクセル値が、シーン内の点を、シミュレートされるカメラの画像平面上に記録されるように表し、
    前記コンピュータグラフィックスシステムが、選択された方向に沿って前記ピクセルからシーン内に向かう少なくとも１つの光線をシミュレートするステップを含む選択されたレイトレーシング法を使用して、画像の前記ピクセル値を生成するように構成され、
    前記レイトレーシング法が、光線と前記シーン内の物体の表面との交わりを計算するステップと、前記シーン内の物体を照明する光線の軌道をシミュレートするステップとをさらに含み、
    軌道をシミュレートする前記ステップが、マルコフ連鎖を使用する、前記コンピュータグラフィックスシステムにおけるシステムであって、
    準モンテカルロ法を使用して前記マルコフ連鎖をシミュレートするための手段を備え、
    マルコフ連鎖をシミュレートするための前記手段が、状態をソートするための手段を備え、
    ソートするための前記手段が、近接性ソーティングのための手段を備える、システム。
  36. 画像内のピクセルのピクセル値を生成するためのコンピュータグラフィックスシステムであり、
    前記ピクセル値が、シーン内の点を、シミュレートされるカメラの画像平面上に記録されるように表し、
    前記コンピュータグラフィックスシステムが、選択された方向に沿って前記ピクセルからシーン内に向かう少なくとも１つの光線をシミュレートするステップを含む選択されたレイトレーシング法を使用して、画像の前記ピクセル値を生成するように構成され、
    前記レイトレーシング法が、光線と前記シーン内の物体の表面との交わりを計算するステップと、前記シーン内の物体を照明する光線の軌道をシミュレートするステップとをさらに含み、
    軌道をシミュレートする前記ステップが、マルコフ連鎖を使用し、
    コンピュータプログラム製品が、コンピュータ可読媒体上に保存されたコンピュータ実行可能プログラムコードを含む、前記コンピュータグラフィックスシステムにおける前記コンピュータ実行可能プログラムコードであって
    準モンテカルロ法を使用して前記マルコフ連鎖をシミュレートするためのコンピュータコード手段を含み、
    マルコフ連鎖をシミュレートするための前記コンピュータコード手段が、状態をソートするためのコンピュータコード手段を含み、
    ソートするための前記コンピュータコード手段が、近接性ソーティングのためのコンピュータコード手段を含む、コンピュータ実行可能プログラムコード。
  37. 非斉次マルコフ連鎖をシミュレートするステップを含む、請求項１９に記載のコンピュータグラフィックスシステムの改良。
JP2009524776A 2006-08-15 2007-08-15 準モンテカルロ法を使用するマルコフ連鎖の同時シミュレーション Expired - Fee Related JP4947394B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US82241706P 2006-08-15 2006-08-15
US60/822,417 2006-08-15
PCT/US2007/075966 WO2008022173A2 (en) 2006-08-15 2007-08-15 Simultaneous simulation of markov chains using quasi-monte carlo techniques

Publications (2)

Publication Number Publication Date
JP2010501100A true JP2010501100A (ja) 2010-01-14
JP4947394B2 JP4947394B2 (ja) 2012-06-06

Family

ID=39083085

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2009524776A Expired - Fee Related JP4947394B2 (ja) 2006-08-15 2007-08-15 準モンテカルロ法を使用するマルコフ連鎖の同時シミュレーション

Country Status (4)

Country Link
EP (1) EP2052366A2 (ja)
JP (1) JP4947394B2 (ja)
CA (1) CA2660190A1 (ja)
WO (1) WO2008022173A2 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011243029A (ja) * 2010-05-19 2011-12-01 International Business Maschines Corporation シミュレーション・システム、方法及びプログラム
JP2013037691A (ja) * 2011-08-04 2013-02-21 Nvidia Corp 加速構造を構築するためのシステム、方法、及びコンピュータプログラムプロダクト

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8131770B2 (en) * 2009-01-30 2012-03-06 Nvidia Corporation System, method, and computer program product for importance sampling of partitioned domains
US10235601B1 (en) 2017-09-07 2019-03-19 7D Labs, Inc. Method for image analysis
US11334762B1 (en) 2017-09-07 2022-05-17 Aurora Operations, Inc. Method for image analysis

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005505810A (ja) * 2001-06-07 2005-02-24 メンタル イメージズ ゲーエムベーハー サンプル点生成のためのリカーシブローテーション（ｒｅｃｕｒｓｉｖｅｒｏｔａｔｉｏｎｓ）を用いた決定論的手法（ｓｔｒｉｃｔｌｙｄｅｔｅｒｍｉｎｉｓｔｉｃｍｅｔｈｏｄｏｌｏｇｉｅｓ）を用いた画像レンダリング

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1264281A4 (en) * 2000-02-25 2007-07-11 Univ New York State Res Found ARRANGEMENT AND METHOD FOR PROCESSING AND PLAYING A VOLUME
US7184071B2 (en) * 2002-08-23 2007-02-27 University Of Maryland Method of three-dimensional object reconstruction from a video sequence using a generic model

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005505810A (ja) * 2001-06-07 2005-02-24 メンタル イメージズ ゲーエムベーハー サンプル点生成のためのリカーシブローテーション（ｒｅｃｕｒｓｉｖｅｒｏｔａｔｉｏｎｓ）を用いた決定論的手法（ｓｔｒｉｃｔｌｙｄｅｔｅｒｍｉｎｉｓｔｉｃｍｅｔｈｏｄｏｌｏｇｉｅｓ）を用いた画像レンダリング
JP2005506607A (ja) * 2001-06-07 2005-03-03 メンタル イメージス， ゲーエムベーハー グローバルイルミネーションを評価するためのロシアンルーレット法を用いた描画

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011243029A (ja) * 2010-05-19 2011-12-01 International Business Maschines Corporation シミュレーション・システム、方法及びプログラム
US8942962B2 (en) 2010-05-19 2015-01-27 International Business Machines Corporation Method and system for measured value simulation
JP2013037691A (ja) * 2011-08-04 2013-02-21 Nvidia Corp 加速構造を構築するためのシステム、方法、及びコンピュータプログラムプロダクト

Also Published As

Publication number Publication date
CA2660190A1 (en) 2008-02-21
JP4947394B2 (ja) 2012-06-06
WO2008022173A2 (en) 2008-02-21
EP2052366A2 (en) 2009-04-29
WO2008022173A3 (en) 2008-12-04

Similar Documents

Publication Publication Date Title
US7773088B2 (en) Simultaneous simulation of markov chains using quasi-monte carlo techniques
US8411088B2 (en) Accelerated ray tracing
US7952583B2 (en) Quasi-monte carlo light transport simulation by efficient ray tracing
US7499053B2 (en) Real-time precision ray tracing
US7659894B2 (en) Terminating spatial partition hierarchies by a priori bounding memory
JP4858795B2 (ja) 瞬時光線追跡
Jönsson et al. Historygrams: Enabling interactive global illumination in direct volume rendering using photon mapping
WO2009044282A2 (en) Quasi-monte carlo light transport simulation by efficient ray tracing
US20060066616A1 (en) Diffuse photon map decomposition for parallelization of global illumination algorithm
Wang et al. GPU-based out-of-core many-lights rendering
JP4947394B2 (ja) 準モンテカルロ法を使用するマルコフ連鎖の同時シミュレーション
EP1899896A2 (en) Real-time precision ray tracing
Shareef et al. Rapid previewing via volume-based solid modeling
Bruckner Dynamic Visibility‐Driven Molecular Surfaces
Krayer et al. Generating signed distance fields on the GPU with ray maps
Morrical et al. Quick clusters: A GPU-parallel partitioning for efficient path tracing of unstructured volumetric grids
Cook et al. Image-space visibility ordering for cell projection volume rendering of unstructured data
Ali et al. Compressed facade displacement maps
Figueiredo et al. An efficient collision detection algorithm for point cloud models
Byrne et al. Applications of the VOLA format for 3D data knowledge discovery
Wächter et al. Efficient simultaneous simulation of Markov chains
Berk et al. Direct volume rendering of unstructured grids
Gao et al. A GPU algorithm for convex hull
Ye et al. 3D Model Occlusion Culling Optimization Method Based on WebGPU Computing Pipeline.
Tuft System for Collision Detection Between Deformable Models Built on Axis Aligned Bounding Boxes and GPU Based Culling

Legal Events

Date Code Title Description
A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20111025

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20120113

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

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20120207

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20120222

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20150316

Year of fee payment: 3

R150 Certificate of patent or registration of utility model

Free format text: JAPANESE INTERMEDIATE CODE: R150

LAPS Cancellation because of no payment of annual fees