Go to Qiita Advent Calendar 2024 Top
LoginSignup
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@Glaesses

C言語で回転する立方体をつくる

Posted at

はじめに

まずは見てみましょう．

一応，この動画を最初から見ればどのような流れでコーディングをしているのかわかるのですが，めちゃ高速なので理解に苦しみます．したがって，実際に手順を踏んでコーディングし，その背後にある数学についても解説します．

前提知識

今回は3次元の回転する立体を2次元に落とし込んで表示させます．したがって，3次元での回転を考える必要があります．ここで便利なのが，回転行列です．

簡単のため，2次元の回転を考えることにします．点$(x,y)$を原点中心に角度$\theta$だけ回転した後の点の座標を求めてみましょう．回転には極座標を用いると簡単です．点$(x,y)$の動径を$r$，偏角を$\alpha$とすると，
$$ x = r\cos\alpha ,\ y = r\sin\alpha$$
となります．これを$\theta$回転した後の点$(x',y')$は，
$$ x' = r\cos(\alpha+\theta) ,\ y' = r\sin(\alpha+\theta)$$
となります．加法定理を用いると，

\begin{align}
x'  &= r\cos(\alpha+\theta)\\
    &= r(\cos\alpha\cos\theta-\sin\alpha\sin\theta)\\
    &= x\cos\theta-y\sin\theta\\
y'  &= r\sin(\alpha+\theta)\\
    &= r(\sin\alpha\cos\theta+\cos\alpha\sin\theta)\\
    &= y\cos\theta+x\sin\theta
\end{align}

ここで行列を用いると，

\begin{align}
\left(\begin{array}{c}
  x' \\
  y' \\
\end{array}\right)
&=
\left(\begin{array}{cc}
  \cos\theta & -\sin\theta\\
  \sin\theta & \cos\theta\\
\end{array}\right)
\left(\begin{array}{c}
  x \\
  y \\
\end{array}\right)\\
&= R(\theta)\left(
\begin{array}{c}
  x \\
  y \\
\end{array}\right)
\end{align}

と表せます．この$R(\theta)$が2次元における回転行列です．

この3次元ver. も存在し，$x$軸，$y$軸，$z$軸周りの回転を表す回転行列は次にようになります．

\begin{align}
R_x(\theta) = 
\left(\begin{array}{ccc}
  1 & 0 & 0 \\
  0 & \cos\theta & -\sin\theta\\
  0 & \sin\theta & \cos\theta\\
\end{array}\right)\\
R_y(\theta) = 
\left(\begin{array}{ccc}
  \cos\theta & 0 & \sin\theta \\
  0 & 1 & 0 \\
  -\sin\theta & 0 & \cos\theta\\
\end{array}\right)\\
R_z(\theta) = 
\left(\begin{array}{ccc}
  \cos\theta & -\sin\theta & 0\\
  \sin\theta & \cos\theta & 0\\
  0 & 0 & 1
\end{array}\right)
\end{align}

ここでの回転の方向は，$R_x$は$y$軸を$z$軸に向ける方向，$R_y$は$z$軸を$x$軸に向ける方向，$R_z$は$x$軸を$y$軸に向ける方向です．（回転行列 - Wikipedia）

勝手に自分でプログラミングする際はこれで十分なんですが，今回はちょっと異なります．というのも，動画では点を回転させるのではなく，点は固定して軸（座標・視点）を回転させているためです．詳しい説明は省きますが，“座標を$\phi$回転させた状態の点の位置”は，“もとの状態の座標で点を$-\phi$回転させた位置”と等しいことが分かれば，簡単です．すなわち，$\sin A$を$\sin(-A) = -\sin A$，$\cos A$を$\cos(-A) = \cos A$に変換すればいいんです．したがって，

\begin{align}
R'_x(\theta) = 
\left(\begin{array}{ccc}
  1 & 0 & 0 \\
  0 & \cos\theta & \sin\theta\\
  0 & -\sin\theta & \cos\theta\\
\end{array}\right)\\
R'_y(\theta) = 
\left(\begin{array}{ccc}
  \cos\theta & 0 & -\sin\theta \\
  0 & 1 & 0 \\
  \sin\theta & 0 & \cos\theta\\
\end{array}\right)\\
R'_z(\theta) = 
\left(\begin{array}{ccc}
  \cos\theta & \sin\theta & 0\\
  -\sin\theta & \cos\theta & 0\\
  0 & 0 & 1
\end{array}\right)
\end{align}

これを掛けていきます！

補足 列ベクトル$(x,y,z)$に右から$R_x, R_y, R_z$をかけていっても，座標を回転できます．動画ではこれを用いています．

以上のことが分かれば，もうコード書けます！

実際に書いてみよう！

1. 回転後の座標を求める

行列の計算は面倒くさいので，Linear Algebra CalculatorというWebツールを用います．

点$(i,j,k)$を$x$軸回りに$A$，$y$軸回りに$B$，$z$軸回りに$C$だけ回転させた点の座標を求めます．

結果はこのようになりました．（動画と同じになるように項順を入れ替えています）

\begin{align}
\left(\begin{array}{ccc}
  \cos C & \sin C & 0 \\ 
  -\sin C & \cos C & 0 \\
  0 & 0 & 1 
\end{array}\right)
\left(\begin{array}{ccc}
  \cos B & 0 & -\sin B \\
  0 & 1 & 0 \\
  \sin B & 0 & \cos B
\end{array}\right)
\left(\begin{array}{ccc}
  1 & 0 & 0 \\
  0 & \cos A & \sin A \\
  0 & -\sin A & \cos A
\end{array}\right)
\left(\begin{array}{c}
  i\\ j\\ k
\end{array}\right)\\
=
\left(\begin{array}{c}
  j \sin A \sin B \cos C - k \cos A \sin B \cos C + j \cos A \sin C + k \sin A \sin C + i \cos B \cos C \\
  j \cos A \cos C + k \sin A \cos C - j \sin A \sin B \sin C + k \cos A \sin B \sin C - i \cos B \sin C \\
  k \cos A \cos B - j \sin A \cos B + i \sin B 
\end{array}\right)
\end{align}

以上より，

float calculateX(int i, int j, int k) {
  return j * sin(A) * sin(B) * cos(C) - k * cos(A) * sin(B) * cos(C) +
         j * cos(A) * sin(C) + k * sin(A) * sin(C) + i * cos(B) * cos(C);
}

float calculateY(int i, int j, int k) {
  return j * cos(A) * cos(C) + k * sin(A) * cos(C) -
         j * sin(A) * sin(B) * sin(C) + k * cos(A) * sin(B) * sin(C) -
         i * cos(B) * sin(C);
}

float calculateZ(int i, int j, int k) {
  return k * cos(A) * cos(B) - j * sin(A) * cos(B) + i * sin(B);
}

というようにコーディングできます．$A, B, C$は回転角度，$i, j, k$は点の座標です．

2. 3Dを2Dに

次に，3次元上にある点を2次元に落とし込みます．動画では，立方体の$z$座標にdistanceFromCam（初期値100）を足し，$xy$平面に映し出しています．求める式は次の通りです．$height, width$は投影する画面の縦横サイズで，$cubeWidth$は立方体の1辺の半分の長さです．

\begin{align}
xp &= \frac{width}{2} - 2 \times cubeWidth + \frac{K1 \times x \times 2}{z}\\
yp &= \frac{height}{2} + \frac{K1 \times y}{z}
\end{align}

両式とも，1項目の$\frac{width}{2}, \frac{height}{2}$は，投影画面の中心を起点に持ってくるためです．上式の2項目$- 2 \times cubeWidth$も位置調整です．すこし左に寄せています．最後の項における$K1$は表示倍率です．大きく表示したければこの値を大きく，小さくしたければこの値を小さくすればできます．遠い点ほど小さく表示させたいので，$z$で割って（動画では$ozz = 1/z$を掛けて）います．$x$座標にだけ$2$を掛けているのは，表示する文字が高さより幅が短いためです．幅を2倍して，正方形に近づくようにしています．

立方体ですから，手前側にある面だけが表示されてほしいです．そのためにzBufferという配列に$z$座標の逆数を記録しておきます．もし新たな点の表示位置が重なれば，zBufferの値と比較し，新たな点の方が大き（i.e. $z$座標が小さ，近）ければ更新します．

3. 描写！

考え方は以上の通りです．これらを駆使して，描写していきます．

変数は以下の通りです．

float A, B, C;                  // 座標の回転角

float cubeWidth = 20;           // 立方体の1辺の長さの半分
int width = 160, height = 44;   // screenの縦と横の長さ
float zBuffer[160 * 44];        // z座標の逆数を記録するバッファ
char buffer[160 * 44];          // 表示する文字を記録するバッファ
int backgroundASCIICode = ' ';  // 背景の文字
int distanceFromCam = 100;      // カメラ（視点・screen）から立方体までの距離（z座標）
float K1 = 40;                  // 表示倍率（拡大係数）

float incrementSpeed = 0.6;

float x, y, z;
float ooz;
int xp, yp;
int idx;

まず，screenをきれいにします．

printf("\x1b[2J");

こうすることで，全面がまっさらな状態になります．

次に，配列の初期化です．

memset(buffer, backgroundASCIICode, width * height);
memset(zBuffer, 0, width * height * 4);

zBufferの方はchar型なので4倍してあります．値は0，すなわち無限遠として初期化します．

次に，立方体の各点の計算です．

for (float cubeX = - cubeWidth; cubeX < cubeWidth; cubeX += incrementSpeed) {
  for (float cubeY = - cubeWidth; cubeY < cubeWidth; cubeY += incrementSpeed) {
    calculateForSurface(cubeX, cubeY, -cubeWidth, '.');
    calculateForSurface(cubeWidth, cubeY, cubeX, '$');
    calculateForSurface(-cubeWidth, cubeY, -cubeX, '~');
    calculateForSurface(-cubeX, cubeY, cubeWidth, '#');
    calculateForSurface(cubeX, -cubeWidth, -cubeY, ';');
    calculateForSurface(cubeX, cubeWidth, cubeY, '+');
  }
}

void calculateForSurface(float cubeX, float cubeY, float cubeZ, int ch) {
  x = calculateX(cubeX, cubeY, cubeZ);
  y = calculateY(cubeX, cubeY, cubeZ);
  z = calculateZ(cubeX, cubeY, cubeZ) + distanceFromCam;

  ooz = 1 / z;

  xp = (int)(width / 2 - 2 * cubeWidth + K1 * ooz * x * 2);
  yp = (int)(height / 2 + K1 * ooz * y);

  idx = xp + yp * width;
  if (idx >= 0 && idx < width * height) {
    if (ooz > zBuffer[idx]) {
      zBuffer[idx] = ooz;
      buffer[idx] = ch;
    }
  }
}

各面を.$~#;+で塗分けます．$z$軸を垂直上方向とすると，底面が.，上面が#という具合です．idxは，配列が1次元なのでそのインデックスを求めています．

最後に実際に標準出力に表示させます．

printf("\x1b[H");
for (int k = 0; k < width * height; k++) {
  putchar(k % width ? buffer[k] : 10);
}
A += 0.005;   // 回転
B += 0.005;   // 回転
usleep(1000); // 描画速度調整

1行目でカーソルを左上に移動します．上書きするような感じです．screenの幅分表示したら改行（ASCIIで10）しています．

参考

https://youtu.be/p09i_hoFdd0
https://rikei-tawamure.com/entry/2019/11/04/184049
https://ja.wikipedia.org/wiki/%E5%9B%9E%E8%BB%A2%E8%A1%8C%E5%88%97
https://www.symbolab.com/solver/linear-algebra-calculator
https://www.mm2d.net/main/legacy/c/c-06.html
https://github.com/servetgulnaroglu/cube.c

1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?