【UnityShader】ロドリゲスの回転公式【2】#82

前回の成果

ロドリゲスの回転公式のベクトル表現を理解した。

soramamenatan.hatenablog.com

今回やること

ロドリゲスの回転公式のベクトル表現を行列で表現します。

正射影ベクトル

まずは、 $r \cos\theta + (1 - \cos\theta)(r \cdot n)n + (n \times r) \sin\theta$ の $(r \cdot n)n$ を行列表現します。
そのために、ベクトル射影の行列表現というものを使用します。
公式は以下となります。

正射影ベクトルの公式ベクトル $\displaystyle \vec{b}$ を $\displaystyle \vec{a}$ が定める直線に正射影したベクトルは、
$\displaystyle \frac{\vec{a} \cdot \vec{b}}{|\vec{a}|^2} \vec{a}$

では証明していこうと思います。

証明

以下の画像を用いて証明します。

正射影ベクトルの公式の証明と使い方 | 高校数学の美しい物語：より引用

$\vec{v}$ を求める

f:id:soramamenatan:20201128162236p:plain

まず、求めたいベクトルを $\vec{v}$ と置きます。
$\vec{v}$ と $\vec{a}$ は平行ですので、実数 $k$ を利用すると

$\vec{v}$ を求める $\displaystyle \vec{v} = k \vec{a} \\ |\vec{v}| = \pm k |\vec{a}|$

とすることができます。

$\vec{a}$ と $\vec{b}$ の内積

f:id:soramamenatan:20201128162718p:plain

次に内積について考えます。

$|\vec{v}| = \pm k |\vec{a}|$ ですので、以下となります。

$\vec{a} \cdot \vec{b}$ を求める $\displaystyle \begin{align} \vec{a} \cdot \vec{b} &= |\vec{a}||\vec{b}| \cos \theta \\ &= \pm |\vec{a}||\vec{v}| \\ \end{align}$

実数 $k$ を求める

$|\vec{v}|$ と $\vec{a} \cdot \vec{b}$ を求めることができたので、 $k$ を求めます。

$k$ を求める $\displaystyle \begin{align} \vec{a} \cdot \vec{b} &= \pm |\vec{a}||\vec{v}| \\ |\vec{v}| &= \pm \frac{\vec{a} \cdot \vec{b}}{|\vec{a}|} \\ \end{align}$
$\displaystyle |\vec{v}| = \pm k |\vec{a}|$ より、
$\displaystyle \begin{align} \pm k |\vec{a}| &= \pm \frac{\vec{a} \cdot \vec{b}}{|\vec{a}|} \\ k &= \pm \frac{\vec{a} \cdot \vec{b}}{|\vec{a}|^2} \end{align}$

代入する

これで $k$ を求めることができました。
$\vec{v} = k \vec{a}$ なので以下となり、証明完了です。

証明 $\displaystyle \begin{align} \vec{v} &= k \vec{a} \\ &= \frac{\vec{a} \cdot \vec{b}}{|\vec{a}|^2} \vec{a} \end{align}$

正射影ベクトルの行列表現

正射影ベクトルの原理が理解できたので、これを行列で表します。

ベクトル射影の表現行列：より引用

ベクトル $r$ のベクトル $a$ への射影

ベクトル $r$ のベクトル $a$ への射影は以下となります。

ベクトル $r$ のベクトル $a$ への射影 $\displaystyle \begin{align} r_{a} &= \frac{r \cdot a}{a^2}a \\ &= \frac{r_{1}a_{1} + r_{2}a_{2} + r_{3}a_{3}}{a^2}a \\ &= \frac{1}{a^2} \begin{pmatrix} (r_{1}a_{1} + r_{2}a_{2} + r_{3}a_{3})a_{1} \\ (r_{1}a_{1} + r_{2}a_{2} + r_{3}a_{3})a_{2} \\ (r_{1}a_{1} + r_{2}a_{2} + r_{3}a_{3})a_{3} \end{pmatrix} \\ &= \frac{1}{a^2} \begin{pmatrix} r_{1}a^2_{1} & r_{2}a_{1}a_{2} & r_{3}a_{1}a_{3} \\ r_{1}a_{1}a_{2} & r_{2}a^2_{2} & r_{3}a_{2}a_{3} \\ r_{1}a_{1}a_{3} & r_{2}a_{2}a_{3} & r_{3}a^2_{3} \end{pmatrix} \\ &= \frac{1}{a^2} \begin{pmatrix} a^2_{1} & a_{1}a_{2} & a_{1}a_{3} \\ a_{1}a_{2} & a^2_{2} & a_{2}a_{3} \\ a_{1}a_{3} & a_{2}a_{3} & a^2_{3} \end{pmatrix} \begin{pmatrix} r_{1} \\ r_{2} \\ r_{3} \end{pmatrix} \\ &= Ar \end{align}$

正射影ベクトルの証明や、行列表現が長かったため忘れがちですが、
これにより、 $(r \cdot n)n$ を行列表現することができます。

$(r \cdot n)n$ の行列表現 $\displaystyle \begin{align} \frac{r \cdot a}{a^2}a &= \frac{1}{a^2} \begin{pmatrix} a^2_{1} & a_{1}a_{2} & a_{1}a_{3} \\ a_{1}a_{2} & a^2_{2} & a_{2}a_{3} \\ a_{1}a_{3} & a_{2}a_{3} & a^2_{3} \end{pmatrix} \begin{pmatrix} r_{1} \\ r_{2} \\ r_{3} \end{pmatrix} \\ \end{align}$
より、
$\displaystyle \begin{align} (r \cdot n)n &= \begin{pmatrix} n^2_{1} & n_{1}a_{2} & n_{1}n_{3} \\ n_{1}n_{2} & n^2_{2} & n_{2}n_{3} \\ n_{1}n_{3} & n_{2}a_{3} & n^2_{3} \end{pmatrix} r \end{align}$

また、この行列は転置行列ですので以下のような表し方もできます。

転置行列の表し方 $\displaystyle \begin{align} (r \cdot n)n &= \begin{pmatrix} n^2_{1} & n_{1}a_{2} & n_{1}n_{3} \\ n_{1}n_{2} & n^2_{2} & n_{2}n_{3} \\ n_{1}n_{3} & n_{2}a_{3} & n^2_{3} \end{pmatrix} r \\ &= \begin{pmatrix} n_{1} \\ n_{2} \\ n_{3} \end{pmatrix} \begin{pmatrix} n_{1} & n_{2} & n_{3} \end{pmatrix} r \\ =& nn^{\mathrm{T}}r \end{align}$

ベクトル積の表現行列

$r \cos\theta + (1 - \cos\theta)(r \cdot n)n + (n \times r) \sin\theta$ の $(r \cdot n)n$ を行列表現しました。
最後に $n \times r$ を行列表現します。

$n \times r$ の行列表現

ベクトル同士の乗算ですので、外積を表しています。

$n \times r$ の行列表現 $\displaystyle \begin{align} n \times r &= \begin{pmatrix} n_{2}r_{3} - n_{3}r_{2} \\ n_{3}r_{1} - n_{1}r_{3} \\ n_{1}r_{2} - n_{2}r_{1} \end{pmatrix} \\ &= \begin{pmatrix} 0 & -n_{3} & n_{2} \\ n_{3} & 0 & -n_{1} \\ -n_{2} & n_{1} & 0 \end{pmatrix} \begin{pmatrix} r_{1} \\ r_{2} \\ r_{3} \end{pmatrix} \\ \end{align}$

また、この行列は歪対称行列ですので以下で表すこともできます。

歪対称行列の表し方 $\displaystyle \begin{align} n \times r &= \begin{pmatrix} 0 & -n_{3} & n_{2} \\ n_{3} & 0 & -n_{1} \\ -n_{2} & n_{1} & 0 \end{pmatrix} \begin{pmatrix} r_{1} \\ r_{2} \\ r_{3} \end{pmatrix} \\ &= \left\lbrack n \right\rbrack_{\times} r \end{align}$

これで $r \cos\theta + (1 - \cos\theta)(r \cdot n)n + (n \times r) \sin\theta$ を行列表現できます。

$r \cos\theta + (1 - \cos\theta)(r \cdot n)n + (n \times r) \sin\theta$ の行列表現

長くなりましたが、これでロドリゲスの回転行列を表すことができます。

画像をもう一度出してから行列に置き換えていきます。

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

位置ベクトル $r$ を角 $\theta$ だけ回転させる回転行列を $R_{n}(\theta)$ とします。
そうすると、回転後の位置ベクトル $r'$ は以下で表されます。

$r'$ $\displaystyle r' = R_{n}(\theta)r$

次に、ロドリゲスの回転公式のベクトル表現より以下が定義されます。

ロドリゲスの回転公式のベクトル表現 $\displaystyle r' = r \cos\theta + (1 - \cos\theta)(r \cdot n)n + (n \times r) \sin\theta$

そして、ベクトル射影の行列表現より以下となります。

ベクトル射影の行列表現 $\displaystyle \begin{align} (r \cdot n)n &= nn^{\mathrm{T}}r \\ &= \begin{pmatrix} n^2_{1} & n_{1}a_{2} & n_{1}n_{3} \\ n_{1}n_{2} & n^2_{2} & n_{2}n_{3} \\ n_{1}n_{3} & n_{2}a_{3} & n^2_{3} \end{pmatrix} r \\ \end{align}$

更に、ベクトル積の行列表現により以下になります。

ベクトル積の行列表現 $\displaystyle \begin{align} n \times r &= \left\lbrack n \right\rbrack_{\times} r \\ &= \begin{pmatrix} 0 & -n_{3} & n_{2} \\ n_{3} & 0 & -n_{1} \\ -n_{2} & n_{1} & 0 \end{pmatrix} r \end{align}$

ベクトル射影の行列表現と、ベクトル積の行列表現により、ロドリゲスの回転公式のベクトル表現は以下に置き換えることができます。
$I$ は単位行列です。

ロドリゲスの回転公式の置き換え $\displaystyle \begin{align} r' =( \cos \theta I + (1 - \cos \theta) nn^{\mathrm{T}} + \sin \theta \left\lbrack n \right\rbrack_{\times} ) r \end{align}$

よって、以下でロドリゲスの回転公式を行列表現できます。

ロドリゲスの回転公式の行列表現 $\displaystyle \begin{align} R_{n}(\theta) &= \cos \theta I + (1 - \cos \theta) nn^{\mathrm{T}} + \sin \theta \left\lbrack n \right\rbrack_{\times} \\ &= \cos \theta \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} + (1 - \cos \theta) \begin{pmatrix} n^2_{1} & n_{1}a_{2} & n_{1}n_{3} \\ n_{1}n_{2} & n^2_{2} & n_{2}n_{3} \\ n_{1}n_{3} & n_{2}a_{3} & n^2_{3} \end{pmatrix} + \sin \theta \begin{pmatrix} 0 & -n_{3} & n_{2} \\ n_{3} & 0 & -n_{1} \\ -n_{2} & n_{1} & 0 \end{pmatrix} \\ &= \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} \end{align}$

これでロドリゲスの回転公式の行列表現が求めることができました。

結果

Thetaを変更することによって、Axisの軸で回転すれば成功です。

f:id:soramamenatan:20201128192114g:plain

inspector

f:id:soramamenatan:20201128192116p: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);
                // 0行列を回避
                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
        }
    }
}