正規直交基底の構築, 再考

2017-09-07
#math

“Building an Orthonormal Basis, Revisited“という、1つの単位ベクトルから正規直交基底を求める方法について書かれた論文を読んでみた。

Revisitedというだけあって、Frisvadが提案した方法だと誤差が大きくなるケースがあるのを修正する、という内容。 そちらも知らなかったのであたってみた。

ナイーブな方法

1つの単位ベクトル から3次元の正規直交基底を求めるナイーブな方法は、平行しない適当なベクトルを選び、直交成分を取り出して正規化したものを とし、また外積を取ったものを として求める：

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.9f) b1 = Vec3f(0.0f, 1.0f, 0.0f);
  else            b1 = Vec3f(1.0f, 0.0f, 0.0f);
  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
}

Hughes-Möller法

適当なベクトルとして、 のx, y, z成分の中で絶対値が一番小さい成分を0として、残り2つを入れ替え、片方の符号を反転したものを選ぶと に直交する、ということを利用してナイーブな方法よりちょっとお得：

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.0f);
  else                       b2 = Vec3f(0.0f, -n.z, n.y);
  b2 *= rsqrt(dot(b2, b2));
  b1 = cross(b2, n);
}

Frisvad法

クォータニオンによる任意軸の回転を利用して、 を回転させた基底を計算する。 正規化をしなくてすむのでお得：

// Listing 2. New way of finding an orthonormal basis from a unit 3D vector.
void frisvad(const Vec3f& n, Vec3f& b1, Vec3f& b2) {
  if (n.z < -0.9999999f) {  // Handle the singularity
    b1 = Vec3f(0.0f, -1.0f, 0.0f);
    b2 = Vec3f(-1.0f, 0.0f, 0.0f);
    return;
  }
  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);
}

Frisvad法の問題点

が-1に近づくと誤差がひどくなり、また行列式が負になることもある。

修正方法

z成分が-1に近づくと問題が起きるが、そうでなければ問題ない。 じゃあということでzが負の場合には逆にから任意軸回転させてやればよい：

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);
  }
}

しかしこの変更によってFrisvad法より2倍も遅くなってしまう！ これはFrisvad法の場合には条件分岐がほとんど成り立たないので分岐予測がほぼ当たるのに対し、修正版では半々になってしまうため。 これを低減するためにcopysignfを使い分岐をなくす：

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);
}

これにより5%~12%遅い程度で済む。

参照

関連記事

新着記事