Frisvad の方法# Frisvad の方法を確認する Frisvad [2012] では, 接空間の z z z n ⃗ \vec{n} n ⊥ ω ⃗ {}_{\perp}\vec{\omega} ⊥ ω ω ⃗ \vec{\omega} ω s \boldsymbol{s} s t \boldsymbol{t} t q ^ \hat{q} q ^
q ^ = ( q v , q w ) = ( s × t 2 ( 1 + ⟨ s , t ⟩ ) , 2 ( 1 + ⟨ s , t ⟩ ) 2 ) (1) \hat{q} =(q_{v}, q_{w}) =\left(\frac{\boldsymbol{s} \times \boldsymbol{t}}{\sqrt{2(1 + \langle \boldsymbol{s}, \boldsymbol{t} \rangle)}}, \frac{\sqrt{2(1 + \langle \boldsymbol{s}, \boldsymbol{t} \rangle)}}{2} \right) \tag{1} q ^ = ( q v , q w ) = ( 2 ( 1 + ⟨ s , t ⟩) s × t , 2 2 ( 1 + ⟨ s , t ⟩) ) ( 1 ) z z z ( 0 , 0 , 1 ) T (0,0,1)^{\mathsf{T}} ( 0 , 0 , 1 ) T n ⃗ = ( n x , n y , n z ) T \vec{n}=(n_{x}, n_{y}, n_{z})^{\mathsf{T}} n = ( n x , n y , n z ) T ( 1 ) (1) ( 1 )
q ^ = ( ( − n y , n x , 0 ) T 2 ( 1 + n z ) , 2 ( 1 + n z ) 2 ) (2) \hat{q} =\left(\frac{(-n_{y}, n_{x}, 0)^{\mathsf{T}}}{\sqrt{2(1 + n_{z})}}, \frac{\sqrt{2(1 + n_{z})}}{2} \right) \tag{2} q ^ = ( 2 ( 1 + n z ) ( − n y , n x , 0 ) T , 2 2 ( 1 + n z ) ) ( 2 ) 代数的な操作の後, n ⃗ \vec{n} n
ω ⃗ = ( q ^ ( ⊥ ω ⃗ , 0 ) q ^ ∗ ) v = ( x ( 1 − n x 2 1 + n z − n x n y 1 + n z − n x ) + y ( − n x n y 1 + n z 1 − n y 2 1 + n z − n y ) + z ( n x n y n z ) ) \begin{align*} \vec{\omega} &=(\hat{q} ({}_{\perp}\vec{\omega}, 0) \hat{q}^{*})_{v}\\ &= \begin{pmatrix} x \begin{pmatrix} 1 - \frac{n^{2}_{x}}{1 + n_{z}}\\ - \frac{n_{x} n_{y}}{1 + n_{z}}\\ - n_{x} \end{pmatrix} + y \begin{pmatrix} - \frac{n_{x} n_{y}}{1 + n_{z}}\\ 1 - \frac{n^{2}_{y}}{1 + n_{z}}\\ - n_{y} \end{pmatrix} + z \begin{pmatrix} n_{x}\\ n_{y}\\ n_{z} \end{pmatrix} \end{pmatrix} \tag{3} \end{align*} ω = ( q ^ ( ⊥ ω , 0 ) q ^ ∗ ) v = x 1 − 1 + n z n x 2 − 1 + n z n x n y − n x + y − 1 + n z n x n y 1 − 1 + n z n y 2 − n y + z n x n y n z ( 3 ) この方程式は, 1 + n z = 0 1 + n_{z}=0 1 + n z = 0 n ⃗ = ( 0 , 0 , − 1 ) T \vec{n}=(0, 0, -1)^{\mathsf{T}} n = ( 0 , 0 , − 1 ) T ω ⃗ = ( − y , − x , − z ) T \vec{\omega}=(-y, -x , -z)^{\mathsf{T}} ω = ( − y , − x , − z ) T
z z z n ⃗ \vec{n} n ( 3 ) (3) ( 3 )
b ⃗ 1 = ( 1 − n x 2 1 + n z , − n x n y 1 + n z , − n x ) T b ⃗ 2 = ( − n x n y 1 + n z , 1 − n y 2 1 + n z , − n y ) T n ⃗ = ( n x , n y , n z ) T \begin{align*} \vec{b}_{1} &= \left( 1 - \frac{n^{2}_{x}}{1 + n_{z}}, - \frac{n_{x} n_{y}}{1 + n_{z}}, - n_{x} \right)^{\mathsf{T}} \tag{4}\\ \vec{b}_{2} &= \left( - \frac{n_{x} n_{y}}{1 + n_{z}}, 1 - \frac{n^{2}_{y}}{1 + n_{z}}, - n_{y} \right)^{\mathsf{T}} \tag{5}\\ \vec{n} &= \left( n_{x}, n_{y}, n_{z} \right)^{\mathsf{T}} \tag{6}\\ \end{align*} b 1 b 2 n = ( 1 − 1 + n z n x 2 , − 1 + n z n x n y , − n x ) T = ( − 1 + n z n x n y , 1 − 1 + n z n y 2 , − n y ) T = ( n x , n y , n z ) T ( 4 ) ( 5 ) ( 6 ) これらの式は, n ⃗ \vec{n} n z z z
NOTE 特異点への対処について: 1 + n z = 0 1 + n_{z}=0 1 + n z = 0 b ⃗ 1 = ( 0 , − 1 , 0 ) T \vec{b}_{1}=(0, -1, 0)^{\mathsf{T}} b 1 = ( 0 , − 1 , 0 ) T b ⃗ 2 = ( − 1 , 0 , 0 ) T \vec{b}_{2}=(-1, 0, 0)^{\mathsf{T}} b 2 = ( − 1 , 0 , 0 ) T
ω ⃗ = ( b ⃗ 1 , b ⃗ 2 , n ⃗ ) ⊥ ω ⃗ = ( 0 − 1 0 − 1 0 0 0 0 − 1 ) ( x y z ) = ( − y − x − z ) \vec{\omega} =(\vec{b}_{1}, \vec{b}_{2}, \vec{n}) {}_{\perp}\vec{\omega} = \begin{pmatrix} 0 & -1 & 0\\ -1 & 0 & 0\\ 0 & 0 & -1 \end{pmatrix} \begin{pmatrix} x\\ y\\ z \end{pmatrix} = \begin{pmatrix} -y\\ -x\\ -z \end{pmatrix} ω = ( b 1 , b 2 , n ) ⊥ ω = 0 − 1 0 − 1 0 0 0 0 − 1 x y z = − y − x − z 実装は次のようになる.
// 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.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); } Issues# Max [2017] によると, 単精度浮動小数点数の場合に特異点付近で許容できない誤差が発生することが報告されている. また, 最大誤差を最小化することにも触れており,
As t increases away from -1, the maximum errors for vectors n with n.z ≥ t will decrease, and those errors for n.z < t will increase, so there should be an optimum trade-off value of t where they are equal. When considering only the normality conditions, the maximum normality error was minimized at 0.00222844 when t = -0.9999805689, in which case the maximum orthogonality error was 0.00622611. When considering only the orthogonality conditions, the maximum error was minimized to the value 0.004901 when t = -0.99998796, in which case the maximum normality condition error was 0.00180161. This same value of t minimized the maximum errors when both the orthogonality and normality conditions were considered together. So if these errors of about a half of one percent are acceptable, this threshold of -0.99998796 would work in Frisvad’s method.
閾値を次のように変更する.
... if (n.z < - 0.99998796 f ) // Handle the singularity ... 確度の比較 Pros and Cons# Pros外積と正規化がないため, 既存の方法を凌駕する実行速度を得る 実質的に分岐がない 四元数を基礎としているため, 接ベクトルの方向に一貫性がある (Frisvad [2014]) Cons単精度浮動小数点数の場合に, 特異点付近で許容できない誤差が発生する Appendix 1: 最短円弧の四元数について# Melax [2000] を参考に最短円弧の四元数を導出する.
安定した式の導出# 回転を表す単位四元数は, 次のように表される.
q = e θ 2 u ^ = cos ( θ 2 ) + sin ( θ 2 ) u ^ q =e^{\frac{\theta}{2} \hat{\boldsymbol{u}}} =\cos{\left(\frac{\theta}{2}\right)} + \sin{\left(\frac{\theta}{2}\right)} \hat{\boldsymbol{u}} q = e 2 θ u ^ = cos ( 2 θ ) + sin ( 2 θ ) u ^ 2 つのベクトル a ^ , b ^ ∈ R 3 \hat{\boldsymbol{a}}, \hat{\boldsymbol{b}} \in \mathbb{R}^{3} a ^ , b ^ ∈ R 3 ∥ a ^ ∥ = ∥ b ^ ∥ = 1 \|\hat{\boldsymbol{a}}\|=\|\hat{\boldsymbol{b}}\|=1 ∥ a ^ ∥ = ∥ b ^ ∥ = 1 d = ⟨ a ^ , b ^ ⟩ d=\langle \hat{\boldsymbol{a}}, \hat{\boldsymbol{b}} \rangle d = ⟨ a ^ , b ^ ⟩ c = a ^ × b ^ \boldsymbol{c}=\hat{\boldsymbol{a}}\times \hat{\boldsymbol{b}} c = a ^ × b ^ q q q q v = ( q x , q y , q z ) q_{v}=(q_{x}, q_{y}, q_{z}) q v = ( q x , q y , q z ) θ \theta θ sin ( θ 2 ) \sin{\left(\frac{\theta}{2} \right)} sin ( 2 θ ) c \boldsymbol{c} c θ \theta θ sin ( θ ) \sin{(\theta)} sin ( θ )
q v = sin ( θ 2 ) u ^ = sin ( θ 2 ) sin θ c q_{v} =\sin{\left(\frac{\theta}{2} \right)} \hat{\boldsymbol{u}} =\frac{\sin{\left(\frac{\theta}{2} \right)}}{\sin{\theta}} \boldsymbol{c} q v = sin ( 2 θ ) u ^ = sin θ sin ( 2 θ ) c sin ( θ 2 ) sin θ \frac{\sin{\left(\frac{\theta}{2} \right)}}{\sin{\theta}} s i n θ s i n ( 2 θ ) sin ( θ 2 ) \sin{\left(\frac{\theta}{2} \right)} sin ( 2 θ )
sin ( θ 2 ) = 1 − cos θ 2 \sin{\left(\frac{\theta}{2} \right)} =\sqrt{\frac{1 - \cos{\theta}}{2}} sin ( 2 θ ) = 2 1 − cos θ また, 次の三角関数の恒等式を用いる.
sin 2 θ + cos 2 θ = 1 \sin^{2}{\theta} + \cos^{2}{\theta}=1 sin 2 θ + cos 2 θ = 1 すると, 次のようになる.
sin ( θ 2 ) sin θ = 1 − cos θ 2 1 − cos 2 θ \frac{\sin{\left(\frac{\theta}{2} \right)}}{\sin{\theta}} =\frac{\sqrt{\frac{1 - \cos{\theta}}{2}}}{\sqrt{1 - \cos^{2}{\theta}}} sin θ sin ( 2 θ ) = 1 − cos 2 θ 2 1 − c o s θ d d d
1 − d 2 1 − d 2 = 1 − d 2 ( 1 − d 2 ) = 1 − d 2 ( 1 + d ) ( 1 − d ) = 1 2 ( 1 + d ) \frac{\sqrt{\frac{1 - d}{2}}}{\sqrt{1 - d^{2}}} =\sqrt{\frac{1 - d}{2(1 - d^{2})}} =\sqrt{\frac{1 - d}{2(1 + d)(1 - d)}} =\frac{1}{\sqrt{2(1 + d)}} 1 − d 2 2 1 − d = 2 ( 1 − d 2 ) 1 − d = 2 ( 1 + d ) ( 1 − d ) 1 − d = 2 ( 1 + d ) 1 それらを代入するとベクトル部は次のようになる.
q v = c 2 ( 1 + d ) q_{v} =\frac{\boldsymbol{c}}{\sqrt{2(1 + d)}} q v = 2 ( 1 + d ) c スカラー部は, 三角関数の半角公式を用いると簡単になる.
q s = cos ( θ 2 ) = 1 + cos θ 2 = 1 + d 2 q_{s} =\cos{\left(\frac{\theta}{2} \right)} =\sqrt{\frac{1 + \cos{\theta}}{2}} =\sqrt{\frac{1 + d}{2}} q s = cos ( 2 θ ) = 2 1 + cos θ = 2 1 + d 平方根(コード上の sqrt)の呼び出しを減らすため, 平方根内に 2 2 \frac{2}{2} 2 2
q s = 1 + d 2 = 2 2 1 + d 2 = 2 ( 1 + d ) 2 q_{s} =\sqrt{\frac{1 + d}{2}} =\sqrt{\frac{2}{2} \frac{1 + d}{2}} =\frac{\sqrt{2(1 + d)}}{2} q s = 2 1 + d = 2 2 2 1 + d = 2 2 ( 1 + d ) 残存する不安定性条件# もし a ^ = − b ^ \hat{\boldsymbol{a}}=-\hat{\boldsymbol{b}} a ^ = − b ^ d = − 1 d=-1 d = − 1
実装例# vec4 rotation_arc (in vec3 a , in vec3 b ){ vec3 c = cross (a, b); float d = dot (a, b); float s = sqrt ( 2. * ( 1. + d)); return vec4 (c / s, s / 2. ); } 