【UnityShader】ロドリゲスの回転公式【1】#81

前回の成果

回転行列を理解した

soramamenatan.hatenablog.com

今回やること

前回の回転行列をもっと簡単に解きます。

事前準備

Scene上にCubeを配置します。
Cubeに今回制作したマテリアルをアタッチします。
自分は回転を見やすくするためにUnityのデフォルトのDafault-Checker-GrayをTextureにアタッチしています。

f:id:soramamenatan:20201115165516p:plain

ソースコード

Shader "Unlit/RodriguesRotation" {
    Properties {
        _MainTex ("Texture", 2D) = "white" {}
        _Axis ("Axis", Vector) = (0, 0, 0, 0)
        _Theta ("Theta", float) = 0
    }
    SubShader {
        Tags { "RenderType"="Opaque" }
        LOD 100

        Pass {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag

            #include "UnityCG.cginc"

            struct appdata {
                float4 vertex : POSITION;
                float2 uv : TEXCOORD0;
            };

            struct v2f {
                float2 uv : TEXCOORD0;
                float4 vertex : SV_POSITION;
            };

            sampler2D _MainTex;
            float4 _MainTex_ST;
            fixed4 _Axis;
            float _Theta;

            fixed3 rotate(fixed3 pos, fixed3 axis, float theta) {
                fixed3 axisN = normalize(axis);
                if (length(axisN) < 0.001) {
                    return pos;
                }
                half radTheta = radians(theta);
                fixed sinTheta = sin(radTheta);
                fixed cosTheta = cos(radTheta);
                fixed a = 1.0 - cosTheta;

                fixed3x3 m = fixed3x3 (
                    axisN.x * axisN.x * a + cosTheta,
                    axisN.y * axisN.x * a + axisN.z * sinTheta,
                    axisN.z * axisN.x * a - axisN.y * sinTheta,
                    axisN.x * axisN.y * a - axisN.z * sinTheta,
                    axisN.y * axisN.y * a + cosTheta,
                    axisN.z * axisN.y * a + axisN.x * sinTheta,
                    axisN.x * axisN.z * a + axisN.y * sinTheta,
                    axisN.y * axisN.z * a - axisN.x * sinTheta,
                    axisN.z * axisN.z * a + cosTheta
                );
                return mul(m, pos);
            }

            v2f vert (appdata v) {
                v2f o;
                v.vertex.xyz = rotate(v.vertex.xyz, _Axis.xyz, _Theta);
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.uv = TRANSFORM_TEX(v.uv, _MainTex);
                return o;
            }

            fixed4 frag (v2f i) : SV_Target {
                fixed4 col = tex2D(_MainTex, i.uv);
                return col;
            }
            ENDCG
        }
    }
}

前回の回転行列とShaderの中身は大きく変わっていません。
変化させているのは、回転用の関数であるrotateになります。

回転行列の比較

前回のものも今回のものも、結果としてはオブジェクトを回転させるものです。

前回の回転関数

float4 rotate3D(float4 pos) {
    float radX = radians(_ThetaX);
    float radY = radians(_ThetaY);
    float radZ = radians(_ThetaZ);

    float sinX = sin(radX);
    float cosX = cos(radX);

    float4x4 rotX = float4x4(
        float4(1, 0, 0, 0),
        float4(0, cosX, sinX, 0),
        float4(0, -sinX, cosX, 0),
        float4(0, 0, 0, 1)
    );

    float sinY = sin(radY);
    float cosY = cos(radY);

    float4x4 rotY = float4x4(
        float4(cosY, 0, -sinY, 0),
        float4(0, 1, 0, 0),
        float4(sinY, 0, cosY, 0),
        float4(0, 0, 0, 1)
    );

    float sinZ = sin(radZ);
    float cosZ = cos(radZ);

    float4x4 rotZ = float4x4(
        float4(cosZ, sinZ, 0, 0),
        float4(-sinZ, cosZ, 0, 0),
        float4(0, 0, 1, 0),
        float4(0, 0, 0, 1)
    );

    pos = mul(pos, rotX);
    pos = mul(pos, rotY);
    pos = mul(pos, rotZ);
    return pos;
}

今回の回転関数

fixed3 rotate(fixed3 pos, fixed3 axis, float theta) {
    fixed3 axisN = normalize(axis);
    if (length(axisN) < 0.001) {
        return pos;
    }
    half radTheta = radians(theta);
    fixed sinTheta = sin(radTheta);
    fixed cosTheta = cos(radTheta);
    fixed a = 1.0 - cosTheta;

    fixed3x3 m = fixed3x3 (
        axisN.x * axisN.x * a + cosTheta,
        axisN.y * axisN.x * a + axisN.z * sinTheta,
        axisN.z * axisN.x * a - axisN.y * sinTheta,
        axisN.x * axisN.y * a - axisN.z * sinTheta,
        axisN.y * axisN.y * a + cosTheta,
        axisN.z * axisN.y * a + axisN.x * sinTheta,
        axisN.x * axisN.z * a + axisN.y * sinTheta,
        axisN.y * axisN.z * a - axisN.x * sinTheta,
        axisN.z * axisN.z * a + cosTheta
    );
    return mul(m, pos);
}

比較してみると、前回の回転行列は各軸に対して毎回行列の計算をしており少し冗長に感じます。
それに比べ、今回のものは1回で行列の計算を済ませています。

ロドリゲスの回転公式

今回のShaderでは、ロドリゲスの回転公式と言われるものを使用しています。
説明は以下になります。

ベクトル空間において、与えられた回転軸に対して回転を行うための効率的なアルゴリズムを指す。

ロドリゲスの回転公式 - Wikipedia：より引用

そして、式は以下になります。

ロドリゲスの回転公式の行列表現 3次元空間において、原点 $O$ を通る回転軸(軸方向の単位ベクトルを $n$ とする)の周りに、
位置ベクトル $r$ を角 $\theta$ だけ回転させる回転行列を $R_{n}(\theta)$ とする。
そうすると、回転後の位置ベクトル $r'$ は
$r' = R_{n}(\theta)r$
となる。
直交座標系において、回転軸方向の単位ベクトルを $n = (n_{1}, n_{2}, n_{3})$ とする。
$\displaystyle R_{n}(\theta) = \begin{pmatrix} n^{2}_{1}(1 - \cos\theta) + \cos\theta & n_{1}n_{2}(1 - \cos\theta) - n_{3}\sin\theta & n_{1}n_{3}(1 - \cos\theta) + n_{2}\sin\theta \\ n_{1}n_{2}(1 - \cos\theta) + n_{3}\sin\theta & n^{2}_{2}(1 - \cos\theta) + \cos\theta & n_{2}n_{3}(1 - \cos\theta) - n_{1}\sin\theta \\ n_{1}n_{3}(1 - \cos\theta) - n_{2}\sin\theta & n_{2}n_{3}(1 - \cos\theta) + n_{1}\sin\theta & n^{2}_{3}(1 - \cos\theta) + cos\theta \end{pmatrix}$

ロドリゲスのの回転公式の表現行列：より引用

この式をソースコードに起こしたものが今回の回転回数となります。

なぜこの式で求めることが出来るのかを自分なりに理解したので、解説していきます。

ロドリゲスの回転公式のベクトル表現

式が長いので混乱しやすいですが、今回求めたいものは

$r' = R_{n}(\theta)r$
になります。
$r'$ は以下の式で求めることができます。

$r'$ の式 $\displaystyle r' = r \cos\theta + (1 - \cos\theta)(r \cdot n)n + (n \times r) \sin\theta$

まずは $r'$ を求める式をひとつずつ理解していきます。

$r$ を求める

f:id:soramamenatan:20201127122838p:plain

まずは $r$ を分解します。
ベクトルは垂直成分と平行成分の和に置き換えることができます。
ですので、以下の式で表すことができます。

$r$ の置き換え $\displaystyle r = \overrightarrow{OO'} + \overrightarrow{O'P}$

$\overrightarrow{OO'}$ を求める

f:id:soramamenatan:20201127124110p:plain

平行成分である $\overrightarrow{OO'}$ を置き換えます。

単位ベクトルと内積を取ると、単位ベクトルの方向の値が取れる性質があります。
原理までは割愛しますが、詳しくはこちらの画像と画像のサイト様が分かりやすかったです。

単位ベクトルと内積

ベクトル内積の意味をイメージで学ぶ。射影とは？なす角とは？ | ばたぱら：より引用

この性質を利用すると、 $\overrightarrow{OO'}$ は以下で表すことができます。

$\overrightarrow{OO'}$ の置き換え $\displaystyle \overrightarrow{OO'} = (r \cdot n)n$

$\overrightarrow{O'P}$ を求める

f:id:soramamenatan:20201127183503p:plain

垂直成分の $\overrightarrow{O'P}$ についても置き換えます。
まず、 $\overrightarrow{O'P}$ をベクトルの引き算を用いて2つのベクトルに分解します。

ベクトルの引き算

ベクトルの足し算・引き算：より引用

すると、以下のように表すことができます。

$\displaystyle \overrightarrow{O'P} = \overrightarrow{OP} - \overrightarrow{OO'}$

$\overrightarrow{OP}$ は $r$ 、 $\overrightarrow{OO'}$ は平行成分で求めたとおり $(r \cdot n)n$ と置き換えることができます。
ですので、最終的に $\overrightarrow{O'P}$ は以下で表すことができます。

$\overrightarrow{O'P}$ の置き換え $\displaystyle \begin{align} \overrightarrow{O'P} &= \overrightarrow{OP} - \overrightarrow{OO'} \\ &= r - (r \cdot n)n \end{align}$

$\overrightarrow{O'Q}$ を求める

f:id:soramamenatan:20201127185009p:plain

$r$ の平行成分である $\overrightarrow{OO'}$ と垂直成分である $\overrightarrow{O'P}$ に垂直なベクトルである、 $\overrightarrow{O'Q}$ を求めます。

$\overrightarrow{O'Q}$ は $n$ と $\overrightarrow{O'P}$ の外積から求めることが出来るので、

$\displaystyle \overrightarrow{O'Q} = n \times \overrightarrow{O'P}$

となります。

そして、 $\overrightarrow{O'P}$ は $\overrightarrow{OP} - \overrightarrow{OO'}$ で表せると垂直成分を求める時に出したので以下となります。

$\displaystyle \begin{align} \overrightarrow{O'Q} &= n \times \overrightarrow{O'P} \\ &= n \times (\overrightarrow{OP} - \overrightarrow{OO'}) \\ &= n \times \overrightarrow{OP} - n \times \overrightarrow{OO'} \end{align}$

$n$ と $\overrightarrow{OO'}$ は同じ向きですので、外積を求めると0になります。
ですので、最終的に以下となります。

$\overrightarrow{O'Q}$ の置き換え $\displaystyle \begin{align} \overrightarrow{O'Q} &= n \times \overrightarrow{O'P} \\ &= n \times (\overrightarrow{OP} - \overrightarrow{OO'}) \\ &= n \times \overrightarrow{OP} - n \times \overrightarrow{OO'} \\ &= n \times \overrightarrow{OP} - 0 \\ &= n \times \overrightarrow{OP} \\ &= n \times r \end{align}$

$\overrightarrow{O'P'}$ を求める

f:id:soramamenatan:20201127192426p:plain

$r'$ の垂直成分である $\overrightarrow{O'P'}$ を求めます。

点 $O'$ , $P$ , $P'$ , $Q$ は同じ円周上の点ですので、2次元の回転が計算に利用できます。

$\displaystyle \overrightarrow{O'P'} = \overrightarrow{O'P} \cos\theta + \overrightarrow{O'Q} \sin\theta$

そして、 $\overrightarrow{O'P}$ は $r - (r \cdot n )n$ 、 $\overrightarrow{O'Q}$ は $n \times r$ に置き換えられるので以下となります。

$\overrightarrow{O'P'}$ を求める $\displaystyle \begin{align} \overrightarrow{O'P'} &= \overrightarrow{O'P} \cos\theta + \overrightarrow{O'Q} \sin\theta \\ &= (r - (r \cdot n) n) \cos\theta + (n \times r)\sin\theta \end{align}$

$r'$ を求める

いよいよ最後です。
$r'$ は以下で表すことができます。

$\displaystyle \begin{align} r' &= \overrightarrow{OP'} \\ &= \overrightarrow{OO'} + \overrightarrow{O'P'} \end{align}$

そして、 $\overrightarrow{OO'}$ と $\overrightarrow{O'P'}$ は以下で表せます。

$\displaystyle \overrightarrow{OO'} = (r \cdot n) n \\ \overrightarrow{O'P'} = (r - (r \cdot n) n) \cos\theta + (n \times r)\sin\theta$

ですので、求めたかった $r'$ は以下で表せます。

$r'$ を求める $\displaystyle \begin{align} r' &= \overrightarrow{OP'} \\ &= \overrightarrow{OO'} + \overrightarrow{O'P'} \\ &= (r \cdot n) n + (r - (r \cdot n) n) \cos\theta + (n \times r)\sin\theta \\ &= r \cos\theta + (1 - \cos\theta)(r \cdot n)n + (n \times r) \sin\theta \end{align}$