Stupid Spherical Harmonics (SH) Tricks

要約

この論文は同名のGDC 2008の講義の手引である。これは、球面調和関数《Spherical Harmonics》(SH)の簡潔な概要を提供し、インタラクティブグラフィクスにおいて使用可能ないくつかの方法と発生するかもしれない問題を議論する。特に、「SHを用いて効率的にライティングモデルを求める方法」、「“リンギング《ringing》”とは何か、そして、それについてできることは何か」、「SH積の効率的な求め方とそれが使われる所」、といった問題に焦点を当てる。最新バージョンはhttps://www.ppsloan.org/publications/で入手できる。

はじめに

｜ラプラス方程式《Laplace’s equation》[@2]の解である調和関数《harmonic functions》は様々な分野で広く用いられる。球面調和関数は球に制限したときの解である¹。これは物理学において熱方程式《heat equation》 (時間上の温度変化のモデル化[@5; @25])、重力場や電場[@9]のようなポテンシャルの問題を解くために使われてきた。量子化学や量子力学では原子の電子配置《electron configuration》や量子角運動量《quantum angular momentum》をモデル化するためにも用いられてきた[@16; @51]。グラフィクスに近い所では、散乱現象をモデル化するために使われてきた。コンピュータグラフィックスにおいて、これは広く使われてきた。初期には、ボリューメトリック散乱効果[@18]、大局的な影を持たないマイクロファセットBRDFのための環境反射[@6]、ディフューズでないオフライン光輸送シミュレーション[@40]、BRDF表現[@53]、画像の再ライティング[@28]、制御可能なライティングを伴う画像ベースレンダリング[@54; @55]、光源発光のモデル化[@8]、といった所で使われた。時代が下ると、異方性散乱[@50]やコンピュータビジョン[@3]における研究が増えてくる。

本論の中心はインタラクティブグラフィクスに関連する技術にある。ゲームにおいて広く使われてきた最初の論文は効率的に球面調和関数を用いて放射照度環境マップを表現することで対処する[@35]ものであり、遠方の照明の下にあるディフューズオブジェクトのインタラクティブレンダリングを可能にする。これは同じ制約を持つ制限されたBRDFの種類を扱うために拡張された[@36]。Precomputed Radiance Transfer(PRT)[@41; @20; @24]はライティング環境への静的なオブジェクトやシーンの反応をモデル化するものであり、しばしばSHを用いて表現され、ソフトシャドウやディフューズと単純な｜艶やかな《glossy》マテリアルとの相互反射のような複雑な大域照明効果を含む。これは、より一般的なBRDFモデルを扱うため[@20; @23; @42]、表面下散乱を組み込むため[@42]、機械学習による圧縮技術を通じてレンダリング効率を大幅に改善するため[@42]、表面の詳細のような”局所的”なテクスチャをモデル化する様々な技術を扱うため[@43; @44; @45]に拡張された。SHは遠方のライティング環境による単一散乱をモデル化するために使われてきた[@49]。他には、遠方のライティング仮定を｜持ち上げる《lift》ために勾配を用いるとき[@1]、相互反射のサポート[@46; @33]を含めた動的オブジェクトに対処するいくつかの技術において[@56; @37]、一般的なBRDFモデルを持つオブジェクトのシャドウイングをモデル化するための可視性の表現として[@12]、変形可能なオブジェクトからの影をモデル化するためにスケーリング操作を用いるとき[@52]、屈折のパラメータ化として[@11]、法線マップによるLevel of Detailの問題に対処するための技術として[@15]用いられる。

より実践的な論文には、PRTの実装の詳細をカバーするもの[@13]、エンジンにこれらの技術を統合する方法[@30]、放射照度ボリュームに対してSHと勾配を使う方法[@31; @32]、投影に関する実践的な問題とSH係数を効率的に量子化する方法[@21]、解析的なスカイライトモデル[@34]をSHに投影して、モデルのパラメータの関数としてSHライトプロブを求めるために｜大局的多項式フィッティング《global polynomial fit》を用いる良い論文[@14]が含まれる。また、SHを用いて半球上に定義されるよりロバストな投影関数の計算技術《numerical techniques》も調査されてきた[@22]。

リアルタイムグラフィクスでの多くは、可視性、ライティング、反射率として、単に球面関数を表現する便利な方法として用いられる。他に、ウェーブレット[@39]、キューブマップ上のウェーブレット[@27]、球面放射基底関数《spherical radial basis functions》[@9]、など[@26]の使用可能な基底関数が多く存在するが、球面調和関数はこの文書で述べられるであろういくつかの素晴らしい特性を持つ。これらのその他の基底関数がより適切となるシナリオが存在することを強調することは重要である。

球面調和関数は幾分手強そうに見えるかもしれないが、実のところ素直である。これは単位円上のFourier基底に類似した球上のソレであり、数値的に計算し易い。信号処理において用いられるFourier基底のように、発生し得る”リンギング”アーティファクトを最小限にするため、(ビデオゲームでは常に行われるであろう)級数の打ち切りを行うときには気を付けなければならない。本論は、球面調和関数を用いて効率的に光を計算および表現する方法、SH表現から従来の光を取り出す方法、“リンギング”の説明とその影響を最小限にする軽減技術、球面調和関数を用いる関数の積の調査を述べる。そこで述べられるのは、これが有用である点と、最適化する価値がある特殊ケースについてである。

背景

放射照度環境マップ

放射照度《irradiance》環境《environment》マップ《map》はライトプロブとクランプされたコサイン関数との畳み込みによって生成される。これは放射輝度を表すために$\pi$で除算して正規化されるべきである。この畳み込みはSHを用いて効率的に行われることができ[@35]、同様に、直接的にSHから効率的にレンダリングされるのに十分な正確さを持つ。次数3のSHはこのカーネルを近似するのに良い働きをするが、HDR光源を用いているならば、次数5の使用を検討したくなるだろう(次数4のZHの係数はゼロなので、そのバンドは省略できる)。下図はクランプされたコサイン関数と次数3のSH近似の画像であり、赤の曲線がSH近似を示す。左側の図は$\theta$の関数としてのプロットであり、右側の図は関数の絶対値でスケールした極座標プロットである。

（画像）

下図は次数5の投影(青)も含めたプロットである。

（画像）

次数3のSH近似は$\theta = 0$(北極)で$\frac{1}{16}$だけ過剰推定し、南極に大きさ$\frac{1}{16}$を持つ｜不要なローブ《spurious lobe》を持つ。17の値を反射する方向光源は、法線がその方向を指すとき、反対方向に向かって1の値を反射するだろう(0を反射しているべきなのに)。次数5の近似は、同様の方向を指す法線で31を反射する方向光源では-1を反射するであろう、負のローブを持つ。これらの近似は正確であるが、特に非常に明るい光源を用いている場合、誤差を引き起こし得る。

付録A10は放射照度環境マップを効率的に求めるためのシェーダおよびCPUコードを含む。

ライティングモデル

SHにはライディングを表現する方法が様々存在する。最も簡単なものとして単にキューブマップから投影する方法があるが、計算が高価でなく、アーティストへ公開するものとして有用である可能性を秘めた解析モデルも存在する。最近の論文[@14]では実用的なスカイライトモデル[@34]をSHに投影し、そのモデルのパラメータ空間上でSH係数の大局的多項式フィッティングを行う。

float f[], s[];
float fWtSum = 0;
Foreach(cube map face)
    Foreach(texel)
        float fTmp = 1 + u^2 + v^2;
        float fWt = 4 / (sqrt(fTmp) * fTmp);
        EvalSHBasis(texel, s);
        f += t(texel) * fWt * s; // ベクトル
        fWtSum += fWt;
f *= 4 * Pi / fWtSum; // 球の面積

リンギングは、Gibbs現象とも呼ばれ、信号処理における一般的な問題である。不連続性《discontinuity》を持つ信号が(連続的な関数のみを表現できる)有限のFourier基底に投影されるとき、オーバーシュート《overshoot》やアンダーシュート《undershoot》がその不連続性の周囲に発生する。不連続性を持たない関数でも投影が｜切り詰められる《truncated》場合に同様の挙動を示す可能性がある。我々は、ライティングモデルを調査するとき、そして、放射照度環境マップの表現(クランプしたコサイン関数の投影)において、これらの問題をすでに確認していた。同様の問題は表面設計《surface design》において発生する。ここでは、一連の幾何的制約を満たそうとするとき、望ましくない振動《oscilations》が発生する可能性がある。これらの問題には2つの一般的な解法が存在する。

1. 切り詰めた投影係数にシグマファクタを用いて｜窓を掛ける《windowing》方法。これは信号処理において最も一般的な解法であり、球面調和関数で自明に使用できる[@8; @41; @37]。
2. 標準の最小二乗誤差ではなく、変分関数《variational function》のいくつかの形式を最小化する方法(例えば、曲率の大きさを最小化する)。これは｜コンピュータを利用した《computer aided》幾何的設計《geometric design》において一般的に行われるが、球面調和関数を用いて効率的に行われることもできる[@38]。

SH積

SHを用いて表される2つの関数の積のSH表現を計算することはしばしば有用である。例となるシナリオには以下がある。

1. 大きな飛行物体(可視性×光)、または、シーンの単純な可視性モデル(大きな建物など)に基づいてスカイライトモデルに穴を開ける。
2. 可視性関数を乗算する。これは動的で近似的な大域照明を行うときに発生する。
3. SHライトプロブをスケールまたは修正する。0から1までのある定数の乗算は、例えば、雲を近似するのに使われることができる。

周波数領域における積の計算は非常に複雑であり、要するに、2つのSHベクトルと”三重積テンソル”の掛け算になる。このコードは効率的に生成でき[@47]、本論では説明されないだろう。とはいえ、言及する価値のある特殊な場合がいくつかある。

おわりに

球面調和関数は、ゲームでは特にライティングで、極めて有用なツールである。本論が、これの使い方とこれを使う時に当たり得る課題の和らげ方、に光明を投じることを願うばかりである。本論で述べたアイデアを拡張する方法がいくつかある。(窓関数で用いる)窓掛け係数は負(または正)のローブの大きさ、または、ゼロになるべき方向を指すときに反射した放射輝度の大きさを最小化するために解くことができる。光を抽出するとき、非線形最適化に基づくより厳格な技術は使用できるかもしれず[@29; @44]、より一般的なライティングモデル(例えば、光源からの円錐の角度を含む)が抽出できるかもしれない。フィッティング処理に窓を掛けてみることの調査をする価値もあるかもしれない。まず、関数の滑らかなバージョンにフィットさせて、窓掛けの量を｜減らしていく《dial back》ことで、より良い極小値へ効果的に舵取りを行う。内容依存の窓掛けのアイデアは、特にPRTのような技術と統合するとき、具体化する必要がある。

球面調和関数で色々やっているときに非常に貴重であると判明したツールのひとつは記号的計算プログラムである。私はMaple(https://www.maplesoft.com/)を用いたが、Mathematica(http://www.wolfram.com/)のような他のプログラムも同様に機能するだろう。DirectX SDK(http://msdn.microsoft.com/directx)はPRTと放射照度環境マップの両方を使うサンプルを付属した、計算《evaluation》、回転、積、いくつかの解析的ライティングモデルのための関数を持つ。

謝辞

作業の多くは、Microsoft Researchのグラフィクスグループに配属されたときに始まり、他の協力者と、特に、長期に渡る実りの多い論文シリーズを執筆するJohn Snyderと共に行われた。また、Hao Chen(Bungie)、Arn Arndt、James Grieves、Clint Hanson(当時EA)、Loren McQuade(Blizzard在籍中)、Dan Baker(MSとFiraxis)、Alex Evans(当時Lionhead)、Tom Forsyth(Muckyfoot)、Willem de Boer(Muckyfoot)、Alex Manchor Ko(Naughty Dog)、Naty Hoffman(SCEA)、chris Oat(ATI)等、幾人かのゲーム開発者と球面調和関数の議論をさせて頂いた。Jason Sandlin、Ben Luna、Jon SteedはMicrosoftにてサンプルやテストコードに関してやってもらった。そして、SHやその他の話題についてGDAlgorithmsメーリングリストにて議論に参加できたことを嬉しく思う。

References

（訳注：原著を参照してください）

漸化式《recurrence relations》[@48]はSH基底関数の多項式の形式を効率的に求めるのに使用できる。基底関数の定式を再掲する。

単位球上の点を$(x, y, z)$とすると、(Zのみに依存し、$\sin\theta^m$で割られる)Legendre陪多項式はこれらの再帰性《recurrences》を用いて求めることができる($P_0^0 = 1$)。

ここで、$l$は内部ループで増加し、$m$は外部ループで増加する。

三角関数《trigonometric》の加法定理《addition formula》は(球上の座標$(x,y)$から計算できるように$\sin\theta^m$をかけられた)$\phi$依存の項を求めるのに使用できる。ここで、$S(0) = 0, C(0) = 1$である。

元の式では、$\sin|m|\phi$を$S(m)$に、$\cos m\phi$を$C(m)$に置き換えている。一般に、共通の部分式を自然と因数に分解するため、これらの式に基づくコードを生成することはより効率的である。

SH基底関数の多項式の形式は以下にリスト化される。$L$はバンド番号、$M$は基底関数である。Mapleが$L$と$M$の順をごちゃまぜにすることに注意する…

ラジアン角$a$に対する球光源があるとすると、最初の6つのバンドの記号的な積分は以下のようになる。

ラジアン角$a$に対する円錐があるとすると、その光源は北極で強度1を持ち、角度$a$で0に減衰する。6次では、この関数は単精度を用いて約8°以下の角度を求めるべきではない。平滑化関数の微分は北極と$a$で0となる。最初の6つのバンドは以下となる。

方向$d$にある方向光と元のライティング環境との間の近似誤差を最小化する強度$s$を計算できる。

ここで、$L_e$はライティング環境のSH表現であり、$L_d$は方向$d$にあるライティングモデルのSH表現である¹⁴。この解は以下となる。

アンビエント光を追加したいならば、誤差関数を最小化する必要がある。

ライティングに畳み込みを吸収させ、各変数に関して微分すると以下のようになる。

そして、最小値を見つけるために2つの式がゼロ¹⁵に等しいとして解く。$\hat{L}_a(s)$はスカラの関数であり¹⁶、SHの直交性により、その積分はすべて単純なSHベクトルの内積となる。これは以下の式となる。

ここで、

$c$と$a$について解くと、

同じ結果に達するより単純な方法が存在することが判明している。環境と光の両方がDC項を含まない方向光の強度を解き、スケールされた方向光を用いるときに環境のDC項を再構築するアンビエント項を計算する。これは最終的に以下となる。

上記の式において、内積はDC項を無視し、$\hat{L}_e《0》$は環境光のDC項である。2つや3つの光に対して同様の技術を用いるだろう(その場合、光の強度と同時に$a$を扱うと式は｜より汚く《much nastier》なる)。

複数の光は同様のやり方で行われることができ、まず、すべての光ベクトルからDCを取り除き、2つの光で、以下の誤差関数をもたらす。

その後、係数のそれぞれに関して微分して0に対して解く。以下の式で強度を得る。

ここで、

$A$が方向に独立して一定であることは指摘しておく価値がある(これは次数とベクトルが畳み込まれるかどうかに依存する)。そのアンビエント項は以下のようになる。

3つの光では、その強度は以下となる。

ここで、$C = L_{d0} \circ L_{d2}, D = L_{d1} \circ L_{d2}, E = 2BCD + A(A^2 - B^2 - C^2 - D^2)$である。

この行列は対称であり、このアンビエント係数は以下となる。

まず、正規直交基底関数に対して、実際に最小二乗投影が基底関数に対して積分することで達成されることを示そう。以下の関数を最小化する係数ベクトル$c$を見つけたいとする。

これは各係数に関して微分することで行われることができ、第1階微分を0に対して解く¹⁷。

基底関数は正規直交であるので、$\int y_i(s) y_j(s) ds = \delta_{ij}$であることが知られ、この事実を利用して0に対して解くと、以下を得る。

こうして、直接積分が最小二乗の結果をもたらすようにする。

そして、球上で積分された二乗Laplacianに基づくペナルティを含む誤差汎関数を最小化する係数ベクトル$g$を導く。上の式により、純粋な二乗誤差を最小化する係数ベクトル$c$は分かっている。indexingに使われる関数$h$を導入する。

そして、その誤差関数は以下となる。

微分すると、以下を得る。

0に対して解くと、以下を得る。

これはNewtons法を用いて極めて簡単に行われることができる。二乗Laplaciandは以下となる。

$L$の値($L_l = l^2 (l+1)^2$)の配列は静的であり、$B$の値の配列は$B_l = \sum_{m=-1}^l (f_l^m)^2$で計算できる。Newtons法は初期の推定値を取り(この問題では0でうまく機能する)、以下の漸化式を用いて洗練する。

ここで、$f$は根を探索する関数であり、$f'$は一階微分である。$f_l^m$を$\lambda$を含む$c_l^m$で置き換え、因数分解すると、以下となる。

ここで、$\Delta^2$は元の二乗Laplacianである。反復回数が最大に達するか、逐次近似の絶対値がしきい値を下回るまで反復する(実践では1e-6が上手く機能するように思える)。これは次数《degree》$(n-1)^2$の(すべての分母の積で乗算する)多項式として表現できる。ここで、$n$はその次数《order》である。これは2次かそれより下の次数でのみ有用であるだろう。この場合、根の閉形式の解が計算できたとしても、反復による解よりはるかに速いかどうかは明らかでない。

以下の2つは両方とも半球上で定義される関数である。ひとつ目は支配的な地面を持つシーンで役立ち得る定数である。ふたつ目はZでクランプしたコサインである。

// Zにおける半球での6次のライティングの乗算から6次のSH係数を生成する
void HemiMult(float* R, float* L) {
    R[0] = 0.433012702f * L[2] - 0.1653594569f * L[12] + 0.1036445247f * L[30] + 0.5f * L[0];
    R[1] = 0.5f * L[1] + 0.4192627457f * L[5] - 0.1711632992f * L[19];
    R[2] = -0.05412658775f * L[20] + 0.5f * L[2] + 0.2420614591f * L[6] + 0.433012702f * L[0];
    R[3] = 0.5f * L[3] - 0.1711632992f * L[21] + 0.4192627457f * L[7];
    R[4] = 0.5f * L[4] - 0.1713860232f * L[28] + 0.4133986423f * L[10];
    R[5] = -0.06477782793f * L[29] + 0.4192627457f * L[1] + 0.5f * L[5] + 0.2614562582f * L[11];
    R[6] = 0.5f * L[6] - 0.1448476267f * L[30] + 0.2420614591f * L[2] + 0.3697549864f * L[12];
    R[7] = -0.06477782793f * L[31] + 0.4192627457f * L[3] + 0.5f * L[7] + 0.2614562582f * L[13];
    R[8] = 0.5f * L[8] - 0.1713860232f * L[32] + 0.4133986423f * L[14];
    R[9] = 0.5f * L[9] + 0.41015625f * L[17];
    R[10] = 0.2685102947f * L[18] + 0.5f * L[10] + 0.4133986423f * L[4];
    R[11] = 0.360244363f * L[19] + 0.5f * L[11] + 0.2614562582f * L[5];
    R[12] = 0.2790440836f * L[20] + 0.5f * L[12] - 0.1653594569f * L[0] + 0.3697549864f * L[6];
    R[13] = 0.5f * L[13] + 0.360244363f * L[21] + 0.2614562582f * L[7];
    R[14] = 0.2685102947f * L[22] + 0.4133986423f * L[8] + 0.5f * L[14];
    R[15] = 0.41015625f * L[23] + 0.5f * L[15];
    R[16] = 0.408100316f * L[26] + 0.5f * L[16];
    R[17] = 0.2720668773f * L[27] + 0.5f * L[17] + 0.41015625f * L[9];
    R[18] = 0.2685102947f * L[10] + 0.35621916f * L[28] + 0.5f * L[18];
    R[19] = 0.2856107251f * L[29] - 0.1711632992f * L[1] + 0.5f * L[19] + 0.360244363f * L[11];
    R[20] = 0.2790440836f * L[12] + 0.5f * L[20] + 0.3498002708f * L[30] - 0.05412658775f * L[2];
    R[21] = 0.2856107251f * L[31] + 0.360244363f * L[13] + 0.5f * L[21] -
 0.1711632992f * L[3];
    R[22] = 0.5f * L[22] + 0.35621916f * L[32] + 0.2685102947f * L[14];
    R[23] = 0.2720668773f * L[33] + 0.5f * L[23] + 0.41015625f * L[15];
    R[24] = 0.5f * L[24] + 0.408100316f * L[34];
    R[25] = 0.5f * L[25];
    R[26] = 0.408100316f * L[16] + 0.5f * L[26];
    R[27] = 0.5f * L[27] + 0.2720668773f * L[17];
    R[28] = 0.5f * L[28] - 0.1713860232f * L[4] + 0.35621916f * L[18];
    R[29] = -0.06477782793f * L[5] + 0.5f * L[29] + 0.2856107251f * L[19];
    R[30] = 0.3498002708f * L[20] + 0.1036445247f * L[0] - 0.1448476267f * L[6] + 0.5f * L[30];
    R[31] = -0.06477782793f * L[7] + 0.2856107251f * L[21] + 0.5f * L[31];
    R[32] = 0.5f * L[32] - 0.1713860232f * L[8] + 0.35621916f * L[22];
    R[33] = 0.5f * L[33] + 0.2720668773f * L[23];
    R[34] = 0.5f * L[34] + 0.408100316f * L[24];
    R[35] = 0.5f * L[35];
}

// Zにおけるクランプしたコサインでの6次のライティングの乗算から6次のSH係数を生成する
void CosMult(float* R, float* L) {
    // 整形用の定数
    const float T1 = 0.09547032698f;
    const float T2 = 0.1169267933f;
    const float T3 = 0.2581988897f;
    const float T4 = 0.2886751347f;
    const float T5 = 0.2390457218f;
    const float T6 = 0.2535462764f;
    const float T7 = 0.2182178903f;
    const float T8 = 0.1083940384f;
    const float T9 = 0.2519763153f;
    const float TA = 0.2439750183f;
    const float TB = 0.3115234375f;
    const float TC = 0.2512594538f;
    const float TD = 0.31640625f;
    R[0] = 0.25f * L[0] - 0.03125f * L[20] + T4 * L[2] + 0.1397542486f * L[6];
    R[1] = 0.2236067977f * L[5] + T2* L[11] - 0.02896952533f * L[29] + 0.1875f * L[1];
    R[2] = 0.375f * L[2] + T4 * L[0] + T3 * L[6] - 0.01495979856f * L[30] + T1 * L[12];
    R[3] = 0.2236067977f * L[7] - 0.02896952533f * L[31] + T2* L[13] + 0.1875f * L[3];
    R[4] = 0.101487352f * L[18] + 0.1889822365f * L[10] + 0.15625f * L[4];
    R[5] = 0.2236067977f * L[1] + T5 * L[11] + 0.3125f * L[5] + 0.09568319309f * L[19];
    R[6] = T6 * L[12] + T3 * L[2] + 0.3125f * L[6] + 0.113550327f * L[20] + 0.1397542486f * L[0];
    R[7] = T5 * L[13] + 0.09568319309f * L[21] + 0.2236067977f * L[3] + 0.3125f * L[7];
    R[8] = 0.15625f * L[8] + 0.101487352f * L[22] + 0.1889822365f * L[14];
    R[9] = 0.09068895910f * L[27] + 0.1666666667f * L[17] + 0.13671875f * L[9];
    R[10] = T7 * L[18] + 0.2734375f * L[10] + 0.09068895910f * L[28] + 0.1889822365f * L[4];
    R[11] = T8 * L[29] + TA * L[19] + T5 * L[5] + 0.30078125f * L[11] + T2* L[1];
    R[12] = T9 * L[20] + 0.328125f * L[12] + T6 * L[6] + 0.1028316139f * L[30] + T1 * L[2];
    R[13] = 0.30078125f * L[13] + T8 * L[31] + TA * L[21] + T5 * L[7] + T2* L[3];
    R[14] = 0.2734375f * L[14] + 0.09068895910f * L[32] + 0.1889822365f * L[8] + T7 * L[22];
    R[15] = 0.0906889591f * L[33] + 0.13671875f * L[15] + 0.1666666667f * L[23];
    R[16] = 0.1507556723f * L[26] + 0.123046875f * L[16];
    R[17] = 0.24609375f * L[17] + 0.1666666667f * L[9] + 0.201007563f * L[27];
    R[18] = 0.101487352f * L[4] + T7 * L[10] + 0.2302830933f * L[28] + 0.28125f * L[18];
    R[19] = TD * L[19] + 0.09568319309f * L[5] + TA * L[11] + 0.2461829819f * L[29];
    R[20] = T9 * L[12] + TC * L[30] + 0.113550327f * L[6] - 0.03125f * L[0] + TD   * L[20];
    R[21] = 0.2461829819f * L[31] + TA * L[13] + 0.9568319309e-1f * L[7] + TD * L[21];
    R[22] = 0.101487352f * L[8] + 0.28125f * L[22] + T7 * L[14] + 0.2302830933f * L[32];
    R[23] = 0.201007563f * L[33] + 0.1666666667f * L[15] + 0.24609375f * L[23];
    R[24] = 0.1507556723f * L[34] + 0.123046875f * L[24];
    R[25] = 0.1127929688f * L[25];
    R[26] = 0.2255859375f * L[26] + 0.1507556723f * L[16];
    R[27] = 0.9068895910e-1f * L[9] + 0.201007563f * L[17] + 0.2631835938f * L[27];
    R[28] = 0.2302830933f * L[18] + 0.9068895910e-1f * L[10] + 0.30078125f * L[28];
    R[29] = 0.1083940385f * L[11] + 0.2461829819f * L[19] - 0.02896952533f * L[1] + TB * L[29];
    R[30] = 0.322265625f * L[30] + TC * L[20] - 0.01495979856f * L[2] + 0.1028316139f * L[12];
    R[31] = 0.2461829819f * L[21] - 0.02896952533f * L[3] + TB * L[31] + 0.1083940385f * L[13];
    R[32] = 0.09068895910f * L[14] + 0.2302830933f * L[22] + 0.30078125f * L[32];
    R[33] = 0.09068895910f * L[15] + 0.201007563f * L[23] + 0.2631835938f * L[33];
    R[34] = 0.2255859375f * L[34] + 0.1507556723f * L[24];
    R[35] = 0.1127929688f * L[35];
}