4次元行列と透視行列

主に3DCG関係で自分がよく使用する4次元行列関係の数式の自分用メモです。

#表記

M = \begin{pmatrix}
M_{00} & M_{01} & M_{02} & M_{03}\\
M_{10} & M_{11} & M_{12} & M_{13}\\
M_{20} & M_{21} & M_{22} & M_{23}\\
M_{30} & M_{31} & M_{32} & M_{33}
\end{pmatrix}

#行列式
4次元行列の行列式は少し複雑です。
疑似コードで表現します。

function determinat(M)
  d := 0
  for i=0 to 3
    j := (i == 0) ? 1 : 0
    k := (i <= 1) ? 2 : 1
    l := (i <= 2) ? 3 : 2
    a := M[1,j] * (M[2,k]*M[3,l] - M[3,k]*M[2,l])
    a += M[1,k] * (M[2,l]*M[3,j] - M[3,l]*M[2,j])
    a += M[1,l] * (M[2,j]*M[3,k] - M[3,j]*M[2,k])
    if (i mod 2) == 0
      d += M[0,i] * a
    else
      d -= M[0,i] * a
    endif
  next
  return d
end

#逆行列
4次元行列の逆行列も少し複雑です。
疑似コードで表現します。

function inverse(M)
  R := ( 1, 0, 0, 0,
         0, 1, 0, 0,
         0, 0, 1, 0,
         0, 0, 0, 1 )
  d := 1.0 / determinat(M)
  for i=0 to 3
    j0 := (i == 0) ? 1 : 0
    j1 := (i <= 1) ? 2 : 1
    j2 := (i <= 2) ? 3 : 2
    for j=0 to 3
      i0 := (j == 0) ? 1 : 0
      i1 := (j <= 1) ? 2 : 1
      i2 := (j <= 2) ? 3 : 2
      a := M[i0,j0] * (M[i1,j1]*M[i2,j2] - M[i2,j1]*M[i1,j2])
      a += M[i0,j1] * (M[i1,j2]*M[i2,j0] - M[i2,j2]*M[i1,j0])
      a += M[i0,j2] * (M[i1,j0]*M[i2,j1] - M[i2,j0]*M[i1,j1])
      if ((i + j) mod 2) == 0
        R[i,j] := a * d
      else
        R[i,j] := -a * d
      endif
    next
  next
  return R
end

#透視変換（z遠方制限なし）
ほとんど使用しないですが、ｚが無限遠までクリップされない視錐台
z = zs/2 ～ ∞
z座標の手前側の制限については注意が必要。手前側の制限についてはｚ制限付きの一般的な透視変換式から任意の値を導出可能

width : 画面幅
height : 画面高
x0 : 画面の中心ｘ座標
y0 : 画面の中心ｙ座標
zs : 投影面のｚ座標（透視座標と等倍になるｚ）

\begin{align}
ws &= 0.5 \cdot width\\
hs &= 0.5 \cdot height\\
xs &= zs / ws\\
ys &= -zs / hs\\
xo &= (x0 - ws) / ws\\
yo &= -(y0 - hs) / hs\\
Proj &= \begin{pmatrix}
xs & 0 & xo & 0 \\
0 & ys & yo & 0 \\
0 & 0 & 1 & -zs \\
0 & 0 & 1 & 0
\end{pmatrix}
\end{align}

※俺々エンジンの座標系から OpenGL への変換なので注意
　俺々エンジンはｙ座標は下に延びていて上下反転している。

#透視変換
一般的な視錐台。
但し、俺々エンジンの座標系のためｙ座標が反転することに注意。

width : 画面幅
height : 画面高
x0 : 画面の中心ｘ座標
y0 : 画面の中心ｙ座標
zs : 投影面のｚ座標（透視座標と等倍になるｚ）
near : 視錐台の手前側ｚ座標
far : 視錐台の奥側ｚ座標

\begin{align}
ws &= 0.5 \cdot width\\
hs &= 0.5 \cdot height\\
xs &= zs / ws\\
ys &= -zs / hs\\
xo &= (x0 - ws) / ws\\
yo &= -(y0 - hs) / hs\\
zn &= near\\
zf &= 2 \cdot far\\
Proj &= \begin{pmatrix}
xs & 0 & xo & 0 \\
0 & ys & yo & 0 \\
0 & 0 & \frac{zf+zn}{zf-zn} & -\frac{2 \cdot zf \cdot zn}{zf - zn} \\
0 & 0 & 1 & 0
\end{pmatrix}
\end{align}