【Unity】shaderでクォータニオンを使った回転を行う

## はじめに

shaderでクォータニオンを使った回転をします。クォータニオンを数学的に理解するのは難しいですが、使うだけなら簡単です。

本記事ではクォータニオンについての解説は行いません。少し調べるだけでとてもわかりやすい解説がたくさん見つかるからです。

クォータニオンについて知りたい方は以下の記事を読むことをおすすめします。
CGのための数学 09 クォータニオン
https://zenn.dev/mebiusbox/books/132b654aa02124/viewer/2966c7

## shaderでクォータニオンを使う

クォータニオンは軸と角度を決めて回転を行うのでジンバルロックが発生しなかったり、最適化次第で他の回転手法より処理が軽かったりします。

### 回転

```hlsl
float4 quaternion(float rad, float3 axis)
{
    return float4(normalize(axis) * sin(rad * 0.5), cos(rad * 0.5));
}

float3 rotateQuaternion(float rad, float3 axis, float3 pos)
{
    float4 q = quaternion(rad, axis);
    return (q.w*q.w - dot(q.xyz, q.xyz)) * pos + 2.0 * q.xyz * dot(q.xyz, pos) + 2 * q.w * cross(q.xyz, pos);
}

v2f vert (appdata v)
{
    v2f o;
    v.vertex.xyz = rotateQuaternion(_Rad, _Axis, v.vertex.xyz);
    o.vertex = UnityObjectToClipPos(v.vertex);
    o.uv = TRANSFORM_TEX(v.uv, _MainTex);
    return o;
}
```
_Rad = _Time.y , _Axis = (1.0, 0, 0) の時の動作
X軸で回転してますね
![quaternion.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/440297/1c376df8-a2a0-ba74-3311-37fa6d10ba69.gif)


### クォータニオンの合成

```hlsl
float4 quaternion(float rad, float3 axis)
{
    return float4(normalize(axis) * sin(rad * 0.5), cos(rad * 0.5));
}

float4 mulQuaternion(float4 q, float4 r)
{
    return float4( q.w*r.xyz + r.w*q.xyz + cross(q.xyz, r.xyz), q.w*r.w - dot(q.xyz,r.xyz));
}

float3 rotateQuaternion(float rad, float3 axis, float3 pos)
{
    float4 q = quaternion(rad, axis);
    return (q.w*q.w - dot(q.xyz, q.xyz)) * pos + 2.0 * q.xyz * dot(q.xyz, pos) + 2 * q.w * cross(q.xyz, pos);
}

v2f vert (appdata v)
{
    v2f o;
    float4 q1 = quaternion(_Rad1, _Axis1);
    float4 q2 = quaternion(_Rad2, _Axis2);
    float4 q  = mulQuaternion(q1, q2);
    v.vertex.xyz = rotateQuaternion(q, v.vertex.xyz);
    o.vertex = UnityObjectToClipPos(v.vertex);
    o.uv = TRANSFORM_TEX(v.uv, _MainTex);
    return o;
}
```

_Axis1(1.0, 0, 0) _Rad1(0.785) , _Axis2(0, 1.0, 0) _Rad2(0.785) の時の動作
X軸で45度、Y軸で45度のクォータニオンが合成されています

![mulQuaternion.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/440297/0202a9be-13b4-f6fe-abc9-b626796a8b47.png)



### クォータニオンの球面線形補間

```hlsl
float4 quaternion(float rad, float3 axis)
{
    return float4(normalize(axis) * sin(rad * 0.5), cos(rad * 0.5));
}

float4 sLerpQuaternion(float4 q, float4 r, float t)
{
    float qr = dot(q,r);
    float hs = 1.0 - qr *qr;
    if(hs <= 0.0){
        return q;
    }
    else{
        hs = sqrt(hs);
        if(abs(hs) < 0.0001){
            return q*0.5 + r*0.5;
        }else{
            float theta = acos(qr);
            return sin((1.0 - t)*theta)*q/sin(theta) + r*sin(t*theta)/sin(theta);
        }
    }
    return float4(0,0,0,1);
}

float3 rotateQuaternion(float rad, float3 axis, float3 pos)
{
    float4 q = quaternion(rad, axis);
    return (q.w*q.w - dot(q.xyz, q.xyz)) * pos + 2.0 * q.xyz * dot(q.xyz, pos) + 2 * q.w * cross(q.xyz, pos);
}

v2f vert (appdata v)
{
    v2f o;
    float4 q1 = quaternion(_Rad1, _Axis1);
    float4 q2 = quaternion(_Rad2, _Axis2);
    float4 q  = sLerpQuaternion(q1, q2, _SlerpPer);
    v.vertex.xyz = rotateQuaternion(q, v.vertex.xyz);
    o.vertex = UnityObjectToClipPos(v.vertex);
    o.uv = TRANSFORM_TEX(v.uv, _MainTex);
    return o;
}
```

lerp関数のように第三引数に0から1までの値を使うことで補間ができます

if文が多いので少し処理負荷が高いです

_Axis1(1.0, 0, 0) _Rad1(0.785) , _Axis2(0, 1.0, 0) _Rad2(0.785) ,  _SlerpPer(frac(_Time.y)) の時の動作
X軸で45度回転したときの姿勢とY軸で45度回転したときの姿勢の間を補完しています

![sLerpQuaternion.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/440297/b3973930-b7c3-4382-3887-4403507a9124.gif)


### オイラー角からクォータニオンへの変換

```hlsl
#define PI acos(-1.0)

float3 angleToRadian(float3 angle)
{
    return PI*angle/180.0;
}

// Unityの回転順はZXY
float4 eulerToQuaternion(float3 rad)
{
    rad = rad*0.5;
    return float4(cos(rad.x)*cos(rad.y)*cos(rad.z) + sin(rad.x)*sin(rad.y)*sin(rad.z),
                  sin(rad.x)*cos(rad.y)*cos(rad.z) + cos(rad.x)*sin(rad.y)*sin(rad.z),
                  cos(rad.x)*sin(rad.y)*cos(rad.z) - sin(rad.x)*cos(rad.y)*sin(rad.z),
                  cos(rad.x)*cos(rad.y)*sin(rad.z) - sin(rad.x)*sin(rad.y)*cos(rad.z));
}

float3 rotateQuaternion(float4 q, float3 pos)
{
    return (q.w*q.w - dot(q.xyz, q.xyz)) * pos + 2.0 * q.xyz * dot(q.xyz, pos) + 2 * q.w * cross(q.xyz, pos);
}

v2f vert (appdata v)
{
    v2f o;
    float4 q = eulerToQuaternion(angleToRadian(_Rotation.xyz));
    v.vertex.xyz = rotateQuaternion(q, v.vertex.xyz);
    o.vertex = UnityObjectToClipPos(v.vertex);
    o.uv = TRANSFORM_TEX(v.uv, _MainTex);
    return o;
}
```
オイラー角からクォータニオンに変換しているのでジンバルロックに注意する必要があります。
Unityの回転順はZXYなのでshader側でもZXYの順で回転するようにしています。

_Rotation(45, 45, 45) の時の動作

![eulerToQuaternion.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/440297/b486c685-779d-9065-75de-8b34ebc2a040.png)

## おまけ
### クォータニオンからオイラー角への変換
```hlsl
// クォータニオンをオイラー角(ラジアン)に変換する
float3 quaternionToEuler(float4 qua)
{
    qua.xyzw = qua.yzwx;
    float sx = -(2.0 * qua.y * qua.z - 2.0 * qua.x * qua.w);
    float unlocked = abs(sx) < 0.99999;
    float3 euler = 0;
    euler.x = asin(-(2.0 * qua.y * qua.z - 2.0 * qua.x * qua.w));
    euler.y = unlocked ?
        atan2((2.0 * qua.x * qua.z + 2.0 * qua.y * qua.w), (2.0 * qua.w * qua.w + 2.0 * qua.z * qua.z - 1.0)) :
        atan2(-(2.0 * qua.x * qua.z - 2.0 * qua.y * qua.w), 2.0 * qua.w * qua.w + 2.0 * qua.x * qua.x - 1.0);
    euler.z = unlocked ? 
        atan2((2.0 * qua.x * qua.y + 2.0 * qua.z * qua.w), (2.0 * qua.w * qua.w + 2.0 * qua.y * qua.y - 1.0)) : 0.0;
    return euler;
}
```

## 参考文献
CGのための数学 09 クォータニオン
https://zenn.dev/mebiusbox/books/132b654aa02124/viewer/2966c7
Conversion between quaternions and Euler angles
https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles