正規直交基底の作り方について、改めて勉強しました


先日Youtubeを見ながら、とあるライブコーディングを写経していたのですが、説明が無かった以下の関数がありました。動画のコメント欄に、説明が欠落した旨が記されていたのですが、よくわからなかったので少し調べてみました。この関数は、とある単位ベクトルから、そのベクトルを一つの軸とした、正規直交基底を構築するために、他の二つの互いに直交したベクトルを算出するためのコードです。(Github: mitsuba-rei/rt より引用)

std::tuple tangentSpace(const V& n)
{
  const double s = std::copysign(1, n.z);
  const double a = -1 / (s + n.z);
  const double b = n.x*n.y*a;
  return {
  V(1 + s * n.x*n.x*a,s*b,-s * n.x),
  V(b,s + n.y*n.y*a,-n.y)
  };
}

せっかくなので、自分のわかる範囲から少し順を追って、基底座標系の構築を見てみたいと思います。

Gram Schmidtの正規直交化法

まずはGram Schmidtの正規直交化法を用いて構築する関数です。以下のコードは、[Frisvad 2012]からの引用です。

void naive ( const Vec3f & n , Vec3f & b1 , Vec3f & b2 )
{
  // If n is near the x-axis , use the y- axis . Otherwise use the x- axis .
  if(n.x > 0.9 f ) b1 = Vec3f (0.0 f , 1.0 f , 0.0 f );
  else b1 = Vec3f (1.0 f , 0.0 f , 0.0 f );

  b1 -= n* dot (b1 , n ); // Make b1 orthogonal to n
  b1 *= rsqrt ( dot (b1 , b1 )); // Normalize b1
  b2 = cross (n , b1 ); // Construct b2 using a cross product
}

こちらはベクトルの内積、外積が分かれば、直感的理解できるものとなっています。
Gram Schmidtの正規直交化法を用いて、線形独立な一つの軸を直交化し、正規化することで、基底ベクトルとします。もう一つの基底ベクトルは、外積を用いて算出します。

Hughes and Mollerの手法

次に、Building an Orthonormal Basis from a Unit Vector[Hughes and Moller 99]の手法です。下記のコードは、[Frisvad 2012]からの引用です。

void hughes_moeller ( const Vec3f & n , Vec3f & b1 , Vec3f & b2 )
{
  // Choose a vector orthogonal to n as the direction of b2.
  if( fabs (n.x ) > fabs (n. z )) b2 = Vec3f (-n.y , n .x , 0.0 f );
  else b2 = Vec3f (0.0 f , -n.z , n.y );

  b2 *= rsqrt ( dot (b2 , b2 )); // Normalize b2
  b1 = cross (b2 , n ); // Construct b1 using a cross product
}

こちらも直感的にに理解できると思います。ベクトルnを、xy平面、yz平面いずれか、投影したベクトルが大きくなるほうに投影し、それを90度回転させることで、直交したベクトルを算出します。それを正規化し、一つの基底ベクトルとして、もう一つの軸を外積を用いて算出します。Gram Schmidtの正規直交化法を用いた場合よりは高速ですが、いずれの場合も、ベクトルの正規化のために、rsqrt()の呼び出しが必要となっています。

Frisvadの手法

x軸(1,0,0),y軸(0,1,0),z軸(0,0,1)の基底座標があったとして、その中の一つの軸、例えばz軸を、ベクトルnの方向に向ける回転が計算できれば、それは、nを一つの軸とした、正規直交基底を計算できることと同義です。
任意のベクトルtを任意のベクトルsの向きに回転させるのに使用するquaternionは、以下の式で算出されます。かっこの中の前半が、3要素でs,tに直交したベクトル成分の虚部となり、後半が1要素の実部になります。

\begin{aligned}  \hat{q} = (q_v, q_w) = (\frac{1}{\sqrt{2(1+s \cdot t)}}(s \times t),\frac{\sqrt{2(1+s \cdot t)}}{2})  \end{aligned}

\hat{q} は単位quaternionなので、逆の回転は共役となります。
\begin{aligned}  \hat{q}^{-1} = (-q_v, q_w)  \end{aligned}

これらを用いて、まず、z軸(0,0,1)をベクトルnに回転させるのに使用するquaternionは、以下のように計算されます。
\begin{aligned}  \hat{q}= (\frac{(-n_y, n_x, 0)}{\sqrt{2(1+n_z)}}, \frac{1}{2}\sqrt{2(1+n_z)})  \end{aligned}

このquaternionで表現される回転を、x軸(1, 0, 0), y軸(0,1,0)にそれぞれ適用すれば、正規直交基底のベクトルが算出できるという訳です。実部成分の回転を打ち消しつつ回転させるため、左からqを掛け、右からqの共役を掛けます。

\begin{aligned}  b1 = \hat{q}(1, 0, 0, 0) \hat{q}^{-1} \\  b2 = \hat{q}(0, 1, 0, 0) \hat{q}^{-1}  \end{aligned}

このquaternionの計算は、回転される軸の成分にゼロが多いため、非常に簡単になります。
\begin{aligned}  b1 = (1 - n^2_x / (1+n_z), -n_x n_y / (1+n_z), -n_x) \\  b2 = (-n_x n_y / (1+n_z), 1 - n^2_y / (1+n_z), -n_y)  \end{aligned}

これをプログラムにすると以下のようになります。単位ベクトルに回転をかける演算なので、正規化の処理が必要ありません。rsqrt()を呼び出さないため、高速に実行できます。
また、式から明らかなとおり、この式には特異点があり、nが(0,0,-1)の時に、ゼロ除算となってしまいます。また、その周囲では、扱う数値が、極小、極大化する傾向があるので、浮動小数点の演算では精度が著しく低下します。そのため、それを避けるための条件式が必要となります。

void frisvad ( const Vec3f & n , Vec3f & b1 , Vec3f & b2 )
{
  if(n.z < -0.9999999 f) // Handle the singularity
  {
    b1 = Vec3f ( 0.0 f , -1.0 f , 0.0 f );
    b2 = Vec3f ( -1.0 f , 0.0 f , 0.0 f );
    return ;
  }

  const float a = 1.0 f /(1.0 f + n.z );
  const float b = -n.x*n .y*a ;
  b1 = Vec3f (1.0 f - n .x*n. x*a , b , -n .x );
  b2 = Vec3f (b , 1.0 f - n .y*n. y*a , -n .y );
}

Duff et al.,の手法について

Frisvadの手法は、コストの高いrsqrt()の呼び出しが無く、非常に早く計算できる優れた方法ですが、nが(0,0,-1)近傍の際に発生する誤差、直交性及びに、誤った系(左手系、右手系)になるケースに関する指摘がBuilding an Orthonormal Basis, Revisited[Duff et al., 2017]でなされています。この論文で示される解決法はシンプルで、ベクトルnのz成分が負の時は、nを(0,0,1)に回転するためのquaternionを求める代わりに、nを(0,0,-1)に回転するためのquaternionを求め、その結果算出された基底ベクトルを反転することで、Frisvadの手法にある特異点を回避しています。

void revisedONB(const Vec3f &n, Vec3f &b1, Vec3f &b2)
{
  if(n.z<0.){
    const float a = 1.0f / (1.0f - n.z);
    const float b = n.x * n.y * a;
    b1 = Vec3f(1.0f - n.x * n.x * a, -b, n.x);
    b2 = Vec3f(b, n.y * n.y*a - 1.0f, -n.y);
  } else{
    const float a = 1.0f / (1.0f + n.z);
    const float b = -n.x * n.y * a;
    b1 = Vec3f(1.0f - n.x * n.x * a, b, -n.x);
    b2 = Vec3f(b, 1.0f - n.y * n.y * a, -n.y);
  }
}

しかし、このコードは、ランダムなベクトルnを用いたテストでは、分岐の方向が偏らないため投機実行がうまく動作しないせいか、Frisvadのオリジナルコードのほぼ2倍ほどの処理時間となってしまったそうです。しかし、条件分岐を避けて以下のように書くことで、Frisvadのオリジナルのコードとほぼ同等の処理時間で、基底座標系の構築ができるそうです。

void branchlessONB(const Vec3f &n, Vec3f &b1, Vec3f &b2)
{
  float sign = copysignf(1.0f, n.z);
  const float a = -1.0f / (sign + n.z);
  const float b = n.x * n.y * a;
  b1 = Vec3f(1.0f + sign * n.x * n.x * a, sign * b, -sign * n.x);
  b2 = Vec3f(b, sign + n.y * n.y * a, -n.y);
}

これは、最初に示した、mituba-reiが使用している関数と等価です。

参照

Building an Orthonormal Basis from a 3D Unit Vector Without Normalization[Frisvad, 2012]
http://orbit.dtu.dk/files/126824972/onb_frisvad_jgt2012_v2.pdf

Building an Orthonormal Basis, Revisited[Duff et al., 2017]
https://graphics.pixar.com/library/OrthonormalB/paper.pdf

CG技術系バーチャルYoutuber、レイトレーシングしてみた[三葉レイのCG技術チャンネル]
https://youtu.be/4XeJEDuhyPs
Github: mitsuba-rei/rt
https://github.com/mitsuba-rei/rt

広告