Zenn
Open7

SceneKitのスニペット・メモ [MSL/Swift]

metal
SceneKit
poteChpoteCh

前に投稿した記事 で投影テクスチャマッピングのことを書いていたときに、MSLには逆行列を作る関数がないなと思って調べたので、そのコードを残しておきます。
解説を足したり他のコードスニペットもまとめたりしてから記事にしようかと考えていたけど忘れそうなので備忘録も兼ねて。

追記：コメントが増えて長くなってきたので簡易目次を入れます。

環境

  • Swift Playgrounds 4.5.4
  • Mac Catalyst 15.2 （iOS 15.2・macOS 12.1 互換）
  • Metal Shading Language 3.2
  • Xcode 16.2 (16C5032a)

ドキュメント

poteChpoteCh

MSLで逆行列を作る関数

メモ

OpenGL ES 3.0には inverse 関数がある。

mat2 inverse(mat2 m);
mat3 inverse(mat3 m);
mat4 inverse(mat4 m);

MSLには determinant 関数と transpose 関数はあるが inverse 相当の関数がない。
逆行列の計算は極力避けるようにするか、SceneKitで SCNMatrix4Invert 関数や simd_float4x4.inverse プロパティ（2x2、3x3などにもある）が使えるので、そちら側であらかじめ準備した逆行列の変数をシェーダーへ渡す流れを想定していそうだが、シェーダー側で計算したいケースも出てくるかもしれないのでコードを書いてみる。

float2x2

\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}^{-1} = \frac{1}{ a_{11}a_{22} - a_{12}a_{21} } \times \begin{bmatrix} a_{22} & -a_{12} \\ -a_{21} & a_{11} \end{bmatrix} 
float2x2 inverse(float2x2 mA) {
    return 1.0 / determinant(mA) * float2x2(
         mA[1].y, -mA[0].y,
        -mA[1].x,  mA[0].x
    );
}

float3x3

\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}^{-1} = \frac{1}{ \begin{split} a_{11}a_{22}a_{33} + a_{12}a_{23}a_{31} + a_{13}a_{21}a_{32} \\ - a_{13}a_{22}a_{31} - a_{12}a_{21}a_{33} - a_{11}a_{23}a_{32} \end{split} } \\[1.0em] \times \begin{bmatrix} a_{22}a_{33} - a_{23}a_{32} & -(a_{12}a_{33} - a_{13}a_{32}) & a_{12}a_{23} - a_{13}a_{22} \\ -(a_{21}a_{33} - a_{23}a_{31}) & a_{11}a_{33} - a_{13}a_{31} & -(a_{11}a_{23} - a_{13}a_{21}) \\ a_{21}a_{32} - a_{22}a_{31} & -(a_{11}a_{32} - a_{12}a_{31}) & a_{11}a_{22} - a_{12}a_{21} \end{bmatrix} 
float3x3 inverse(float3x3 mA) {
    return 1.0 / determinant(mA) * float3x3(
          mA[1].y * mA[2].z - mA[1].z * mA[2].y,   -(mA[0].y * mA[2].z - mA[0].z * mA[2].y),    mA[0].y * mA[1].z - mA[0].z * mA[1].y,
        -(mA[1].x * mA[2].z - mA[1].z * mA[2].x),    mA[0].x * mA[2].z - mA[0].z * mA[2].x,   -(mA[0].x * mA[1].z - mA[0].z * mA[1].x),
          mA[1].x * mA[2].y - mA[1].y * mA[2].x,   -(mA[0].x * mA[2].y - mA[0].y * mA[2].x),    mA[0].x * mA[1].y - mA[0].y * mA[1].x
    );
}

float4x4

\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}^{-1} = \frac{1}{ \begin{split} a_{11}( a_{22}a_{33}a_{44} + a_{23}a_{34}a_{42} + a_{24}a_{32}a_{43} - a_{24}a_{33}a_{42} - a_{23}a_{32}a_{44} - a_{22}a_{34}a_{43} ) \\ + a_{12}( a_{21}a_{34}a_{43} + a_{23}a_{31}a_{44} + a_{24}a_{33}a_{41} - a_{24}a_{31}a_{43} - a_{23}a_{34}a_{41} - a_{21}a_{33}a_{44} ) \\ + a_{13}( a_{21}a_{32}a_{44} + a_{22}a_{34}a_{41} + a_{24}a_{31}a_{42} - a_{24}a_{32}a_{41} - a_{22}a_{31}a_{44} - a_{21}a_{34}a_{42} ) \\ + a_{14}( a_{21}a_{33}a_{42} + a_{22}a_{31}a_{43} + a_{23}a_{32}a_{41} - a_{23}a_{31}a_{42} - a_{22}a_{33}a_{41} - a_{21}a_{32}a_{43} ) \end{split} } \\[1.0em] \times \begin{bmatrix} a_{22}a_{33}a_{44} + a_{23}a_{34}a_{42} + a_{24}a_{32}a_{43} - a_{24}a_{33}a_{42} - a_{23}a_{32}a_{44} - a_{22}a_{34}a_{43} & -a_{12}a_{33}a_{44} - a_{13}a_{34}a_{42} - a_{14}a_{32}a_{43} + a_{14}a_{33}a_{42} + a_{13}a_{32}a_{44} + a_{12}a_{34}a_{43} & a_{12}a_{23}a_{44} + a_{13}a_{24}a_{42} + a_{14}a_{22}a_{43} - a_{14}a_{23}a_{42} - a_{13}a_{22}a_{44} - a_{12}a_{24}a_{43} & -a_{12}a_{23}a_{34} - a_{13}a_{24}a_{32} - a_{14}a_{22}a_{33} + a_{14}a_{23}a_{32} + a_{13}a_{22}a_{34} + a_{12}a_{24}a_{33} \\ -a_{21}a_{33}a_{44} - a_{23}a_{34}a_{41} - a_{24}a_{31}a_{43} + a_{24}a_{33}a_{41} + a_{23}a_{31}a_{44} + a_{21}a_{34}a_{43} & a_{11}a_{33}a_{44} + a_{13}a_{34}a_{41} + a_{14}a_{31}a_{43} - a_{14}a_{33}a_{41} - a_{13}a_{31}a_{44} - a_{11}a_{34}a_{43} & -a_{11}a_{23}a_{44} - a_{13}a_{24}a_{41} - a_{14}a_{21}a_{43} + a_{14}a_{23}a_{41} + a_{13}a_{21}a_{44} + a_{11}a_{24}a_{43} & a_{11}a_{23}a_{34} + a_{13}a_{24}a_{31} + a_{14}a_{21}a_{33} - a_{14}a_{23}a_{31} - a_{13}a_{21}a_{34} - a_{11}a_{24}a_{33} \\ a_{21}a_{32}a_{44} + a_{22}a_{34}a_{41} + a_{24}a_{31}a_{42} - a_{24}a_{32}a_{41} - a_{22}a_{31}a_{44} - a_{21}a_{34}a_{42} & -a_{11}a_{32}a_{44} - a_{12}a_{34}a_{41} - a_{14}a_{31}a_{42} + a_{14}a_{32}a_{41} + a_{12}a_{31}a_{44} + a_{11}a_{34}a_{42} & a_{11}a_{22}a_{44} + a_{12}a_{24}a_{41} + a_{14}a_{21}a_{42} - a_{14}a_{22}a_{41} - a_{12}a_{21}a_{44} - a_{11}a_{24}a_{42} & -a_{11}a_{22}a_{34} - a_{12}a_{24}a_{31} - a_{14}a_{21}a_{32} + a_{14}a_{22}a_{31} + a_{12}a_{21}a_{34} + a_{11}a_{24}a_{32} \\ -a_{21}a_{32}a_{43} - a_{22}a_{33}a_{41} - a_{23}a_{31}a_{42} + a_{23}a_{32}a_{41} + a_{22}a_{31}a_{43} + a_{21}a_{33}a_{42} & a_{11}a_{32}a_{43} + a_{12}a_{33}a_{41} + a_{13}a_{31}a_{42} - a_{13}a_{32}a_{41} - a_{12}a_{31}a_{43} - a_{11}a_{33}a_{42} & -a_{11}a_{22}a_{43} - a_{12}a_{23}a_{41} - a_{13}a_{21}a_{42} + a_{13}a_{22}a_{41} + a_{12}a_{21}a_{43} + a_{11}a_{23}a_{42} & a_{11}a_{22}a_{33} + a_{12}a_{23}a_{31} + a_{13}a_{21}a_{32} - a_{13}a_{22}a_{31} - a_{12}a_{21}a_{33} - a_{11}a_{23}a_{32} \end{bmatrix} 
float4x4 inverse(float4x4 mA) {
    return 1.0 / determinant(mA) * float4x4(
         mA[1].y * mA[2].z * mA[3].w + mA[1].z * mA[2].w * mA[3].y + mA[1].w * mA[2].y * mA[3].z - mA[1].w * mA[2].z * mA[3].y - mA[1].z * mA[2].y * mA[3].w - mA[1].y * mA[2].w * mA[3].z,
        -mA[0].y * mA[2].z * mA[3].w - mA[0].z * mA[2].w * mA[3].y - mA[0].w * mA[2].y * mA[3].z + mA[0].w * mA[2].z * mA[3].y + mA[0].z * mA[2].y * mA[3].w + mA[0].y * mA[2].w * mA[3].z,
         mA[0].y * mA[1].z * mA[3].w + mA[0].z * mA[1].w * mA[3].y + mA[0].w * mA[1].y * mA[3].z - mA[0].w * mA[1].z * mA[3].y - mA[0].z * mA[1].y * mA[3].w - mA[0].y * mA[1].w * mA[3].z,
        -mA[0].y * mA[1].z * mA[2].w - mA[0].z * mA[1].w * mA[2].y - mA[0].w * mA[1].y * mA[2].z + mA[0].w * mA[1].z * mA[2].y + mA[0].z * mA[1].y * mA[2].w + mA[0].y * mA[1].w * mA[2].z,

        -mA[1].x * mA[2].z * mA[3].w - mA[1].z * mA[2].w * mA[3].x - mA[1].w * mA[2].x * mA[3].z + mA[1].w * mA[2].z * mA[3].x + mA[1].z * mA[2].x * mA[3].w + mA[1].x * mA[2].w * mA[3].z,
         mA[0].x * mA[2].z * mA[3].w + mA[0].z * mA[2].w * mA[3].x + mA[0].w * mA[2].x * mA[3].z - mA[0].w * mA[2].z * mA[3].x - mA[0].z * mA[2].x * mA[3].w - mA[0].x * mA[2].w * mA[3].z,
        -mA[0].x * mA[1].z * mA[3].w - mA[0].z * mA[1].w * mA[3].x - mA[0].w * mA[1].x * mA[3].z + mA[0].w * mA[1].z * mA[3].x + mA[0].z * mA[1].x * mA[3].w + mA[0].x * mA[1].w * mA[3].z,
         mA[0].x * mA[1].z * mA[2].w + mA[0].z * mA[1].w * mA[2].x + mA[0].w * mA[1].x * mA[2].z - mA[0].w * mA[1].z * mA[2].x - mA[0].z * mA[1].x * mA[2].w - mA[0].x * mA[1].w * mA[2].z,

         mA[1].x * mA[2].y * mA[3].w + mA[1].y * mA[2].w * mA[3].x + mA[1].w * mA[2].x * mA[3].y - mA[1].w * mA[2].y * mA[3].x - mA[1].y * mA[2].x * mA[3].w - mA[1].x * mA[2].w * mA[3].y,
        -mA[0].x * mA[2].y * mA[3].w - mA[0].y * mA[2].w * mA[3].x - mA[0].w * mA[2].x * mA[3].y + mA[0].w * mA[2].y * mA[3].x + mA[0].y * mA[2].x * mA[3].w + mA[0].x * mA[2].w * mA[3].y,
         mA[0].x * mA[1].y * mA[3].w + mA[0].y * mA[1].w * mA[3].x + mA[0].w * mA[1].x * mA[3].y - mA[0].w * mA[1].y * mA[3].x - mA[0].y * mA[1].x * mA[3].w - mA[0].x * mA[1].w * mA[3].y,
        -mA[0].x * mA[1].y * mA[2].w - mA[0].y * mA[1].w * mA[2].x - mA[0].w * mA[1].x * mA[2].y + mA[0].w * mA[1].y * mA[2].x + mA[0].y * mA[1].x * mA[2].w + mA[0].x * mA[1].w * mA[2].y,

        -mA[1].x * mA[2].y * mA[3].z - mA[1].y * mA[2].z * mA[3].x - mA[1].z * mA[2].x * mA[3].y + mA[1].z * mA[2].y * mA[3].x + mA[1].y * mA[2].x * mA[3].z + mA[1].x * mA[2].z * mA[3].y,
         mA[0].x * mA[2].y * mA[3].z + mA[0].y * mA[2].z * mA[3].x + mA[0].z * mA[2].x * mA[3].y - mA[0].z * mA[2].y * mA[3].x - mA[0].y * mA[2].x * mA[3].z - mA[0].x * mA[2].z * mA[3].y,
        -mA[0].x * mA[1].y * mA[3].z - mA[0].y * mA[1].z * mA[3].x - mA[0].z * mA[1].x * mA[3].y + mA[0].z * mA[1].y * mA[3].x + mA[0].y * mA[1].x * mA[3].z + mA[0].x * mA[1].z * mA[3].y,
         mA[0].x * mA[1].y * mA[2].z + mA[0].y * mA[1].z * mA[2].x + mA[0].z * mA[1].x * mA[2].y - mA[0].z * mA[1].y * mA[2].x - mA[0].y * mA[1].x * mA[2].z - mA[0].x * mA[1].z * mA[2].y
    );
}

余因子行列の計算式を愚直にコードへおとしただけなので改善の余地はある。

参考サイト

poteChpoteCh

float4x4の逆行列のコードをすこし整理した版

float4x4 inverse(float4x4 mA) {
    float idet = 1.0 / determinant(mA);
    float4x4 c;
    c[0].x =  determinant(float3x3(mA[1].yzw, mA[2].yzw, mA[3].yzw));
    c[0].y = -determinant(float3x3(mA[0].yzw, mA[2].yzw, mA[3].yzw));
    c[0].z =  determinant(float3x3(mA[0].yzw, mA[1].yzw, mA[3].yzw));
    c[0].w = -determinant(float3x3(mA[0].yzw, mA[1].yzw, mA[2].yzw));
    c[1].x = -determinant(float3x3(mA[1].xzw, mA[2].xzw, mA[3].xzw));
    c[1].y =  determinant(float3x3(mA[0].xzw, mA[2].xzw, mA[3].xzw));
    c[1].z = -determinant(float3x3(mA[0].xzw, mA[1].xzw, mA[3].xzw));
    c[1].w =  determinant(float3x3(mA[0].xzw, mA[1].xzw, mA[2].xzw));
    c[2].x =  determinant(float3x3(mA[1].xyw, mA[2].xyw, mA[3].xyw));
    c[2].y = -determinant(float3x3(mA[0].xyw, mA[2].xyw, mA[3].xyw));
    c[2].z =  determinant(float3x3(mA[0].xyw, mA[1].xyw, mA[3].xyw));
    c[2].w = -determinant(float3x3(mA[0].xyw, mA[1].xyw, mA[2].xyw));
    c[3].x = -determinant(float3x3(mA[1].xyz, mA[2].xyz, mA[3].xyz));
    c[3].y =  determinant(float3x3(mA[0].xyz, mA[2].xyz, mA[3].xyz));
    c[3].z = -determinant(float3x3(mA[0].xyz, mA[1].xyz, mA[3].xyz));
    c[3].w =  determinant(float3x3(mA[0].xyz, mA[1].xyz, mA[2].xyz));
    return idet * c;
}

誤差があるので計算結果は整理前のコードと完全に同じではない。

float3x3の逆行列のコードを同様に整理した版

float3x3 inverse(float3x3 mA) {
    float idet = 1.0 / determinant(mA);
    float3x3 c;
    c[0].x =  determinant(float2x2(mA[1].yz, mA[2].yz));
    c[0].y = -determinant(float2x2(mA[0].yz, mA[2].yz));
    c[0].z =  determinant(float2x2(mA[0].yz, mA[1].yz));
    c[1].x = -determinant(float2x2(mA[1].xz, mA[2].xz));
    c[1].y =  determinant(float2x2(mA[0].xz, mA[2].xz));
    c[1].z = -determinant(float2x2(mA[0].xz, mA[1].xz));
    c[2].x =  determinant(float2x2(mA[1].xy, mA[2].xy));
    c[2].y = -determinant(float2x2(mA[0].xy, mA[2].xy));
    c[2].z =  determinant(float2x2(mA[0].xy, mA[1].xy));
    return idet * c;
}

文字数＆行数はむしろ増えるが、ぱっと見でのわかりやすさはこちらに分がありそう。

poteChpoteCh

MSLでの変換行列

ジオメトリシェーダー向けのfloat4x4の変換行列を作る関数

移動 (Translate)

float4x4 translate(float3 _translate) {
    return float4x4(
        1.0, 0.0, 0.0, 0.0,
        0.0, 1.0, 0.0, 0.0,
        0.0, 0.0, 1.0, 0.0,
        _translate.x, _translate.y, _translate.z, 1.0
    );
}

SCNNodeの simdPosition プロパティや position プロパティと同様の変形ができる。

拡大縮小 (Scale)

float4x4 scale(float3 _scale) {
    return float4x4(
        _scale.x, 0.0,      0.0,      0.0,
        0.0,      _scale.y, 0.0,      0.0,
        0.0,      0.0,      _scale.z, 0.0,
        0.0,      0.0,      0.0,      1.0
    );
}

SCNNodeの simdScale プロパティや scale プロパティと同様の変形ができる。

回転 (Rotate)

オイラー角による回転

// X軸を中心に回転（ピッチ）
float4x4 rotate_x(float _angle) {
    return float4x4(
        1.0,  0.0,         0.0,         0.0,
        0.0,  cos(_angle), sin(_angle), 0.0,
        0.0, -sin(_angle), cos(_angle), 0.0,
        0.0,  0.0,         0.0,         1.0
    );
}

// Y軸を中心に回転（ヨー）
float4x4 rotate_y(float _angle) {
    return float4x4(
        cos(_angle), 0.0, -sin(_angle), 0.0,
        0.0,         1.0,  0.0,         0.0,
        sin(_angle), 0.0,  cos(_angle), 0.0,
        0.0,         0.0,  0.0,         1.0
    );
}

// Z軸を中心に回転（ロール）
float4x4 rotate_z(float _angle) {
    return float4x4(
         cos(_angle), sin(_angle), 0.0, 0.0,
        -sin(_angle), cos(_angle), 0.0, 0.0,
         0.0,         0.0,         1.0, 0.0,
         0.0,         0.0,         0.0, 1.0
    );
}

// オイラー角による回転
float4x4 rotate_euler_angles(float3 _angle) {
    return rotate_z(_angle.z) * rotate_y(_angle.y) * rotate_x(_angle.x);
}

SCNNodeの simdEulerAngles プロパティや eulerAngles プロパティと同様の変形ができる。

回転軸と角度の指定による回転

float4x4 rotate_axis(float3 _axis, float _angle) {
    float3 axis = normalize(_axis);
    float cos_angle = cos(_angle);
    float3 inv_cos_angle_axis = (1.0 - cos_angle) * axis;
    float3 sin_angle_axis = sin(_angle) * axis;

    return float4x4(
        inv_cos_angle_axis.x * axis.x + cos_angle,
        inv_cos_angle_axis.x * axis.y + sin_angle_axis.z,
        inv_cos_angle_axis.x * axis.z - sin_angle_axis.y,
        0.0,

        inv_cos_angle_axis.y * axis.x - sin_angle_axis.z,
        inv_cos_angle_axis.y * axis.y + cos_angle,
        inv_cos_angle_axis.y * axis.z + sin_angle_axis.x,
        0.0,

        inv_cos_angle_axis.z * axis.x + sin_angle_axis.y,
        inv_cos_angle_axis.z * axis.y - sin_angle_axis.x,
        inv_cos_angle_axis.z * axis.z + cos_angle,
        0.0,

        0.0, 0.0, 0.0, 1.0
    );
}

SCNNodeの simdRotation プロパティや rotation プロパティと同様の変形ができる。

クォータニオンによる回転

float4x4 quat_to_mat4(float4 _q) {
    return float4x4(
        1.0 - 2.0 * (_q.y * _q.y + _q.z * _q.z),
        2.0 * (_q.x * _q.y + _q.z * _q.w),
        2.0 * (_q.x * _q.z - _q.y * _q.w),
        0.0,

        2.0 * (_q.y * _q.x - _q.z * _q.w),
        1.0 - 2.0 * (_q.x * _q.x + _q.z * _q.z),
        2.0 * (_q.y * _q.z + _q.x * _q.w),
        0.0,

        2.0 * (_q.z * _q.x + _q.y * _q.w),
        2.0 * (_q.z * _q.y - _q.x * _q.w),
        1.0 - 2.0 * (_q.x * _q.x + _q.y * _q.y),
        0.0,

        0.0, 0.0, 0.0, 1.0
    );
}

// 虚数と実数を指定
float4x4 rotate_quat(float _ix, float _iy, float _iz, float _r) {
    return quat_to_mat4(float4(_ix, _iy, _iz, _r));
}

// 角度と回転軸を指定
float4x4 rotate_quat(float _angle, float3 _axis) {
    // SCNNodeのsimdOrientationプロパティと挙動を揃える目的で
    // この関数内では引数で渡された回転軸の値を正規化していない
    float half_angle = _angle * 0.5;
    float4 q = float4(_axis * sin(half_angle), cos(half_angle));
    return quat_to_mat4(q);
}

SCNNodeの simdOrientation プロパティや orientation プロパティと同様の変形ができる。

メモ

SceneKitではローカル座標やビュー座標は 右手座標系 、ビューポート座標やテクスチャ座標は 左手座標系 になっている。
ジオメトリシェーダーで rotate_z(45.0 / 180.0 * M_PI_F) * _geometry.position の変形を行った場合の回転の向きは、Z軸のプラス側（手前側）の位置から原点の方向を見て 反時計回り に45度の回転になる。

剪断 (せん断 / Shear / Skew)

float4x4 shear(float2 xy, float2 yz, float2 zx) {
    return float4x4(
        1.0,   xy[1], zx[0], 0.0,
        xy[0], 1.0,   yz[1], 0.0,
        zx[1], yz[0], 1.0,   0.0,
        0.0,   0.0,   0.0,   1.0
    );
}

参考サイト

poteChpoteCh

MSLでの法線の再計算

ジオメトリシェーダーで頂点座標を動かしたあとに法線を再計算するコード例

// 法線の再計算用に頂点の近傍点を用意する
float delta = 0.01;
float3 temp_tangent = _geometry.tangent.xyz;
float3 temp_binormal = normalize(cross(_geometry.tangent.xyz, _geometry.normal)) * _geometry.tangent.w;
temp_tangent = _geometry.position.xyz + temp_tangent * delta;
temp_binormal = _geometry.position.xyz + temp_binormal * delta;
// ※ _geometry.tangent.w は通常では -1 になっている
// テクスチャ座標のXあるいはYの方向が逆向きに
// マッピングされている場合 +1 になる (SCNCylinderの上面・底面などが該当)

// 頂点の座標を編集する
float4x4 transform = rotate_z(0.25 * 3.1415927) * scale(float3(1.5, 0.5, 1.0));
_geometry.position = transform * _geometry.position;
// ※ rotate_zとscaleの関数については先に投稿したコメントを参照

// 近傍点に対して_geometry.positionと同じように座標の編集を行う
temp_tangent = (transform * float4(temp_tangent, 1.0)).xyz;
temp_binormal = (transform * float4(temp_binormal, 1.0)).xyz;

// 編集した頂点と近傍点の座標から接空間を再計算する
temp_tangent -= _geometry.position.xyz;
temp_binormal -= _geometry.position.xyz;

// 接空間から頂点の接線と法線を再計算する
_geometry.tangent.xyz = normalize(temp_tangent);
_geometry.normal = normalize(cross(temp_binormal, temp_tangent)) * _geometry.tangent.w;

メモ

binormal（bitangent）を計算するところとnormalを再計算するところでは、

// binormalの計算
float3 temp_binormal = normalize(cross(_geometry.tangent.xyz, _geometry.normal)) * _geometry.tangent.w;

// normalの再計算
_geometry.normal = normalize(cross(temp_binormal, temp_tangent)) * _geometry.tangent.w;

としている。
このクロス積の引数の順序を入れ替えて * _geometry.tangent.w も省いた、

// binormalの計算
float3 temp_binormal = normalize(cross(_geometry.normal, _geometry.tangent.xyz));

// normalの再計算
_geometry.normal = normalize(cross(temp_tangent, temp_binormal));

の組み合わせにしても最終的に得られる法線は同じ。
Unityでの法線の再計算に関する記事では後者の順序の書きかたも見かける。ただしSceneKitのジオメトリシェーダーでこのように書くと、binormalの向きが逆方向になってしまうケースが存在する。

別のコード例

ジオメトリの変形が変換行列で処理できるものなら、先のコメントに記載した逆行列を作る関数を利用して法線を再計算する方法もある。

// 変換行列を作る --- (1)
float4x4 transform = rotate_z(0.25 * 3.1415927) * scale(float3(1.5, 0.5, 1.0));
// 変換行列の逆行列の転置行列（逆転置行列）を作る --- (2)
float3x3 transform_it = transpose(
    inverse(
        float3x3(transform[0].xyz, transform[1].xyz, transform[2].xyz)
    )
);

// (1) の変換行列で頂点の座標を編集する
_geometry.position = transform * _geometry.position;
// (2) の変換行列で法線を再計算する
_geometry.normal = normalize(transform_it * _geometry.normal);
poteChpoteCh

MSLでの座標系の変換

座標の変換

変換元 変換先 コード
ローカル座標系 ワールド座標系 scn_node.modelTransform * _geometry.position
ワールド座標系 ローカル座標系 scn_node.inverseModelTransform * world_position
ローカル座標系 ビュー座標系 scn_node.modelViewTransform * _geometry.position
ビュー座標系 ローカル座標系 scn_node.inverseModelViewTransform * float4(_surface.position, 1.0)
ローカル座標系 クリップ座標系 scn_node.modelViewProjectionTransform * _geometry.position
クリップ座標系 ローカル座標系 scn_node.inverseModelViewProjectionTransform * projection_position
ワールド座標系 ビュー座標系 scn_frame.viewTransform * world_position
ビュー座標系 ワールド座標系 scn_frame.inverseViewTransform * float4(_surface.position, 1.0)
ビュー座標系 クリップ座標系 scn_frame.projectionTransform * float4(_surface.position, 1.0)
クリップ座標系 ビュー座標系 scn_frame.inverseProjectionTransform * projection_position
ワールド座標系 クリップ座標系 scn_frame.viewProjectionTransform * world_position
クリップ座標系 ワールド座標系 scn_frame.inverseViewProjectionTransform * projection_position
※ドキュメントにはこの変換行列の記載なし
ビュー座標系
（座標と法線）		 キューブマップ
テクスチャ座標		 (scn_frame.viewToCubeTransform * float4(reflect(_surface.position, _surface.normal), 0.0)).xyz

法線の変換

変換元 変換先 コード
ローカル座標系 ビュー座標系 normalize((scn_node.normalTransform * float4(_geometry.normal, 0.0)).xyz)
ビュー座標系 ローカル座標系 normalize((inverse(scn_node.normalTransform) * float4(_surface.normal, 0.0)).xyz)
※inverse関数については 先に投稿したコメント を参照
ローカル座標系 ワールド座標系 normalize((transpose(scn_node.inverseModelTransform) * float4(_geometry.normal, 0.0)).xyz);
ワールド座標系 ローカル座標系 normalize((transpose(scn_node.modelTransform) * float4(world_normal, 0.0)).xyz);
ワールド座標系 ビュー座標系 normalize((transpose(scn_frame.inverseViewTransform) * float4(world_normal, 0.0)).xyz);
※ビュー変換行列が直交行列なら
normalize((scn_frame.viewTransform * float4(world_normal, 0.0)).xyz);
でも変換できる
ビュー座標系 ワールド座標系 normalize((transpose(scn_frame.viewTransform) * float4(_surface.normal, 0.0)).xyz);
※ビュー変換行列が直交行列なら
normalize((scn_frame.inverseViewTransform * float4(_surface.normal, 0.0)).xyz);
でも変換できる

メモ１

上の表ではコードを短くまとめるために４×４の変換行列のままで法線の計算をしているが、変換行列の左上３×３の成分だけを使って、

normalize(
    inverse(
        float3x3(
            scn_node.normalTransform[0].xyz,
            scn_node.normalTransform[1].xyz,
            scn_node.normalTransform[2].xyz
        )
    ) * _surface.normal
);

のように書いたほうが無駄な計算を減らせる。

メモ２

基本的には座標の変換行列の逆転置行列との乗算で法線を変換できる。
変換行列が直交行列ならば、逆行列 ＝ 転置行列 なので、逆行列の転置行列 ＝ 転置行列の転置行列 ＝ もとの変換行列 になり、（逆転置ではない）座標の変換行列と法線との乗算でも変換ができる。

ビュー変換行列は直交行列になっていることが多いが、カメラの設定によっては直交行列にならないケースも存在する。そういったケースでは、上の表のコードで ワールド座標系 <--> ビュー座標系 の法線の変換ができないことがある（カメラノードに対して、回転＋アスペクト比が変わるような拡大縮小の変形を設定していたりすると、この問題が起こる場合がある）。
その状況下でも scn_node.normalTransform の変換行列を使って ローカル座標系 <--> ビュー座標系 の法線の変換をするぶんには問題ないため、ビュー座標系での法線からワールド座標系での法線を得たいときは、いったんローカル座標系への変換を経由させる ビュー座標系 -> ローカル座標系 -> ワールド座標系 の流れで法線を計算すると簡単で誤差も少なくてよいと思われる。

参考ドキュメント

poteChpoteCh

MSLでの行列のデータ順と変換行列の合成

行列のデータ順

MSLの行列では、格納するデータの順序が 列優先（column-major order） になっている。

float3x2 であれば ３列２行 の行列を表す。

\small \left[ \begin{array}{c|c|c} a_{1} & b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \end{array} \right] \hspace{1.2em}\large = \small\hspace{1.0em} \begin{array}{l} \texttt{float3x2 m = float3x2(} \\ \hspace{2.0em} \def\arraystretch{1.2} \begin{array}{l} \mathtt{a_{1},} & \mathtt{a_{2},} \\ \hline \mathtt{b_{1},} & \mathtt{b_{2},} \\ \hline \mathtt{c_{1},} & \mathtt{c_{2}} \end{array} \\ \texttt{);} \end{array} \hspace{1.0em}\large = \small\hspace{1.0em} \def\arraystretch{1.2} \begin{array}{l} \texttt{float3x2 m;} \\ \texttt{m[0] = float2(} \mathtt{ a_{1},\; a_{2} } \texttt{);} \\ \texttt{m[1] = float2(} \mathtt{ b_{1},\; b_{2} } \texttt{);} \\ \texttt{m[2] = float2(} \mathtt{ c_{1},\; c_{2} } \texttt{);} \\ \vphantom{ } \end{array}

行列とベクトルの計算において、たとえば 行列 m と ベクトル v をそれぞれ仮に

float3x3 m = float3x3(
    a1, a2, a3,
    b1, b2, b3,
    c1, c2, c3
);

float3 v = float3(x, y, z);

としたとき、以下の６種類の計算結果はどれも同じになる。

float3 out1 = m * v;

float3 out2;
out2.x = (a1 * x) + (b1 * y) + (c1 * z);
out2.y = (a2 * x) + (b2 * y) + (c2 * z);
out2.z = (a3 * x) + (b3 * y) + (c3 * z);

float3 out3;
out3.x = (m[0][0] * v.x) + (m[1][0] * v.y) + (m[2][0] * v.z);
out3.y = (m[0][1] * v.x) + (m[1][1] * v.y) + (m[2][1] * v.z);
out3.z = (m[0][2] * v.x) + (m[1][2] * v.y) + (m[2][2] * v.z);

float3 out4;
out4.x = dot( float3(m[0][0], m[1][0], m[2][0]), v );
out4.y = dot( float3(m[0][1], m[1][1], m[2][1]), v );
out4.z = dot( float3(m[0][2], m[1][2], m[2][2]), v );

float3 out5 = float3(0.0);
out5 += float3(a1, a2, a3) * x;
out5 += float3(b1, b2, b3) * y;
out5 += float3(c1, c2, c3) * z;

float3 out6 = float3(0.0);
out6 += m[0] * v.x;
out6 += m[1] * v.y;
out6 += m[2] * v.z;

変換行列の合成

MSLで変換行列を合成する際は、あとから適用したい変換 × さきに適用したい変換 の順序で行列を掛け算する。
変換対象のベクトルは変換行列の後ろから掛ける。

同次座標 p を 拡大縮小 -> 回転 -> 移動 の順で変形させるなら、

// 座標
float4 p = float4(x, y, z, 1.0);
// 拡大縮小の変換行列
float4x4 scale_matrix = scale(_scale.xyz);
// 回転の変換行列
float4x4 rotate_matrix = rotate_axis(_axis.xyz, _angle);
// 移動の変換行列
float4x4 translate_matrix = translate(_translate.xyz);

// 移動 * 回転 * 拡大縮小 * 座標
p = translate_matrix * rotate_matrix * scale_matrix * p;

のように掛けあわせる。
※変換行列を作る関数については 先に投稿したコメント を参照。

なお、この例に関しては次のように書いて計算量を減らすこともできる。

float4 p = float4(x, y, z, 1.0);
float4x4 transform_matrix = rotate_axis(_axis.xyz, _angle);
transform_matrix[0].xyz *= _scale.x;
transform_matrix[1].xyz *= _scale.y;
transform_matrix[2].xyz *= _scale.z;
transform_matrix[3].xyz = _translate.xyz;

p = transform_matrix * p;