ゲーム作りとかCGに関わる数学③(中級)

# ゲーム作りとかCGに関わる数学(初歩)③
この辺から、理系の高校～大学教養課程くらいのレベルの話になってきます。ちょ～っと難しいかもしれないけど頑張りましょう。

## 微分積分
ここは実は微分積分をやってくださいというリクエストがあったので入れてみました。
人によっては「え？ゲームプログラミングに微分積分とか使うの？」という現役のゲームプログラマの方もいるかもしれませんが、実は必要なんですよ。
高校でやる試験問題を解くみたいな事はしませんが、概念の理解は必要になります。
そもそも速度加速度にかんしても微分や積分の概念ですしね。

### 極限について
極限…この辺も、厳密な話をしだすと割と大変な話になっちゃいますので、一般的な高校で教えるような極限の話をしますね。極限と言えば、lim(limitの略ね)なんですがこれは
```math
\displaylines{
\lim_{x\to 0} x\\
\lim_{h\to 0} f(x+h)
}
```
などのように書いて、limの下に書いてある変数を右側の数字に「限りなく近づける」というものです。この「限りなく近づける」というのがまた曲者で

「え？『限りなく』って何？何そのふわっとした説明！？はっきりしてよ！」

と混乱してしまった高校生も少なくないんじゃないかと思います。だから僕は言います「h→0と書いてあったらh=0です」と。厳密な数学クラスタの方からはお叱りを受けそうなんですが、そういう方は一旦この記事をそっ閉じしてください。スルースキルって大事ですよ(笑)

はい、ということで、x→0やh→0と書いてある部分のxやhを0にしてみます。
そうすると

```math
\displaylines{
\lim_{x\to 0} x=0\\
\lim_{h\to 0} f(x+h)=f(x)
}
```
当然その変数は0になり、なくなっちゃいますね？これでいいのだ。極限というのはそう言う事なのです。x→0はx=0って言う事なんだね！

じゃあなんでlimなんてクッソややこしい概念があるの？
というと数学にもプログラムには共通するルールとして「０で割らないでね？」というのがあります。０除算の禁止です。
また「無限大」なんて数値もないですよね。「無量大数」とは違うんです。あれはまぁ一応「値」が定義されていますから(定義が色々あるようですが)

「無量大数」https://ja.wikipedia.org/wiki/%E7%84%A1%E9%87%8F%E5%A4%A7%E6%95%B0

しかし、計算の過程においては0で割ったり、値を∞という事にしたい事があります。
要は禁じ手の計算を大目に見てもらうために仕方なく0や∞を用います。
例えばx→0であれば計算式を変形してなんとかかんとか0除算を回避します。回避が成功したらすみやかにxに0を代入します。

例えば$$g(x)=\frac{f(x)}{x}$$という式があるとします。この式のxを0にしたいのですが、そうすると$$g(0)=\frac{f(0)}{0}$$
となってしまい、0除算となってしまいます。これはマズい。ということで、「ほんとうは0なんだけど、0じゃない事にして計算」するために極限を使います

```math
\lim_{x\to 0} \frac{f(x)}{x}
```

x「うるうる、ぼく0じゃないよ🥺。だからゆるしてね」
数学神「ゆるす」

となります。こうやって数学神が大目に見てる間に分母の0を消しにかかります。例えば$$f(x)=\frac{x^3+2x^2+6x}{x}$$だとします。そうすると約分によって$$f(x)=x^2+2x+6$$になり、分母の0が消失します。こうなればしめたもので、何食わぬ顔で
x「ごめんね、ぼくは0だったんだ😇」と打ち明けます。
そうすると
$$f(0)=0^2+2×0+6=6$$
となります。

無限大のほうは、別に除算しても怒られませんし、$$\lim_{x\to \infty} f(x)$$で$$f(x)=x^2+2x+1$$だったとしても$$\lim_{x\to \infty} x^2+2x+1=\infty$$
となり、無限大は特に値を表さないため、この場合も「xが取り扱えないほど大きい数ならf(x)も取り扱えないほど大きい数」という事で、あまり意味ないんですが無限大の場合は「数列の和」やら「数列の積」として使われることが多いですね。
たとえば$$\sum\limits_{n=1}^\infty n = \infty$$
であれば、nを無限大まで足すわけですから当然結果も無限大(取り扱えない数)になります。
ところがこれがこういう式だと結果は1に収束してしまいます。

$$\sum\limits_{n=1}^\infty \frac{1}{2^n} = 1$$

永遠に足し続けるわけですから無限大になりそうなものなのですが、図にしてみると…
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/30494060-4319-7d44-cb68-c0141aa9cb28.png)
まずn=1のとき1/2です。次は1/4を足すわけですが、これは見方を変えると「1までの距離の半分」とも言えるわけで、結果的に1/2+1/4=3/4になります。
これを繰り返すと、１までの残りの距離の半分をずーっと足していくのですが「１までの残りの距離」の半分であるため１ににじり寄りつつも、１の手前で永遠に足踏みすることとなります。結果的にこれを無限に繰り返せば１になってしまうというわけです。

「アキレスと亀」のパラドックスみたいな感じですね。アキレスが亀に追いついても亀はその1/10だけ進むので、永遠に追いつけないってやつ。
(※あれも「アキレスは亀のいる地点までしか進まない」っていう暗黙の制約があるからで、そんなもん「追い越していいよ」って言うた時点でパラドックスではなくなると思うんですが…)

ともかく、無限大に足していってもこういう特定の値になってしまうのを「収束する」と言います。最初みたいに合計もまた無限大になるものは「発散する」と言います。

まぁ無限大の方はこんなもんでいいでしょう。次に微分の話をしていきますがこれはh→0を用いる話です。

### 微分について
そもそも微分、微分いうとりますけれども、微分って何でしょう？

微分というのは、特定の式 y=f(x) のグラフ上のとある点$$x=p_x$$における「接線の傾き」を計算するためのものです。さらに専門的な学習をしてる人にとっては別の解釈があるかもしれませんが、ひとまず「接線の傾き」という事にしてください。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/472a08b0-30f7-fd23-5afa-b408ccca080b.png)

「何で**傾き**なんて求めるん？」という疑問が起こるかもしれませんが、その疑問はひとまず置いておいて、ともかく「微分＝その地点での傾きを求めるためのもの」だと思っておいてください。

分かりやすい例で言うと、走ってる車の「とある瞬間の速度」を求めるようなときに「とある点の傾き」を使うと思ってください。車って常に60Km/hで走ってるわけじゃなくて、信号前では減速するし信号のない直進道路では60Km/hで走ったりします。

そういうアップダウンをするものですが、瞬間的な速度を求める必要があったりしますよね？そういう時に使うのが微分なんです。

では、この「傾き」はどうやって求めるのでしょうか？これを求める時に使う概念が「極限」なのです。

ちなみに、微分と言うと…
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/1b8c5fe0-151e-885d-4859-c79db3c1c705.png)
このように微分前の乗数が、微分後の係数にかけられて、乗数が１つずつ下がる…というような事を高校の頃に習った人がいるかもしれません。
でもゲームプログラミングにおいては、これを覚えること自体はそんなに重要じゃありません。ただ、この後に解説する事をイメージしやすくするために上のルールで

```math
\displaylines{
f(x)=x^3+2x^2+4x+1 ― ①
とすると\\
f'(x)=\frac{df(x)}{dx}=3x^2+4x+4
}
```
となる事をとりあえず頭に入れておいてください。で、微分というのは先ほども説明したように、クッソ微小区間の変化率です。
```math
f'(x)=\frac{df(x)}{dx}=\lim_{h\to 0}\frac{f(x+h)-(x)}{x+h-x}=\lim_{h\to 0}\frac{f(x+h)-f(x)}{h}
```
となります。ここで先ほどの①の式を当てはめて、さらに展開すると
```math
\displaylines{
\lim_{h\to 0}\frac{f(x+h)-f(x)}{h}=\\
\lim_{h\to 0}\frac{(x+h)^3+2(x+h)^2+4(x+h)+1-(x^3+2x^2+4x+1)}{h}=\\
\lim_{h\to 0}\frac{x^3+3x^2h+3xh^2+h^3+2x^2+4xh+2h^2+4x+4h+1-(x^3+2x^2+4x+1)}{h}=\\
\lim_{h\to 0}\frac{3x^2h+3xh^2+h^3+4xh+2h^2+4h}{h}=\\
\lim_{h\to 0}3x^2+3xh+h^2+4x+2h+4
}
```
となります。ここで、先ほどの説明通り分母のhがなくなってしまったので安心してh=0にしてみると残るのは$$3x^2+4x+4$$となり、所謂「公式」で計算したものと同じになります。ちなみにこれが$$(x^n)'$$に対しても常に成り立つのかというと成り立ちます。一応導出を書いておくと
```math
\displaylines{
\lim_{h\to 0}\frac{(x+h)^n-x^n}{h}=\\
\lim_{h\to 0}\frac{x^n+nx^{n-1}h+...+nxh^{n-1}+h^n-x^n}{h}=\\
\lim_{h\to 0}\frac{nx^{n-1}h+...+nxh^{n-1}+h^n}{h}=\\
\lim_{h\to 0}nx^{n-1}+...+nxh^{n-2}+h^{n-1}=\\
nx^{n-1}
}
```
となるわけです。1行目から2行目がよく分からない人は「二項定理」で検索してみてください。中学か高校の時にちょっとだけやった覚えがある法則が出てくるはずです。
[二項定理](https://rikeilabo.com/commentary-binomial-theorem)

微分についてもうちょっとだけ続きますが、忘れないように改めて書いておくと、プログラミングにおける微分について大事なのはこういう公式を覚える事じゃないです。とにかく「瞬間的な微小変化を求める」「接線を求める」事が目的だという事は忘れないでください。

次に合成関数の微分についてちょっとだけやってみましょう。あんまり必要じゃないんですが、抽象化の訓練だと思ってお付き合いください。
今までf(x)ばかりに登場してもらいましたが、兄弟のg(x)くんにも登場してもらいましょう。ここでは適当に

```math
\displaylines{
f(x)=x^2+3x+3\\
g(x)=2x+5
}
```
とでもしましょうか。まずわかりやすい例として$$(f(x)g(x))'$$について考えてみましょう。同行する前に最初に愚直に計算しちゃいます。
```math
\displaylines{
f(x)g(x)=2x^3+11x^2+21x+15\\
(f(x)g(x))'=6x^2+22x+21
}
```
はい、この結果を覚えておいてください。ここでまた極限による計算をしますが、具体的な関数ではなく抽象化されたf(x)g(x)のまま取り扱います。
```math
\lim_{h \to 0}\frac{f(x+h)g(x+h)-f(x)g(x)}{h}
```
ちょっとここで「ズルい」事やります。数学ではよくやるやつなんですが「特定の形にしたいために、敢えて一旦は形を複雑にしてしまう」って奴です。有理化とかの時によくやったアレです。収支があってれば問題にしないアレです。
具体的には真ん中にf(x+h)g(x)をねじ込んでみます。そうすると
```math
\displaylines{
\lim_{h \to 0}\frac{f(x+h)g(x+h)-f(x+h)g(x)+f(x+h)g(x)-f(x)g(x)}{h}\\
\lim_{h \to 0}\frac{f(x+h)(g(x+h)-g(x))+g(x)(f(x+h)-f(x))}{h}\\
\lim_{h \to 0}\frac{f(x+h)(g(x+h)-g(x))}{h}+\frac{g(x)(f(x+h)-f(x))}{h}
}
```
はい、この形にした時点で、前半分はg(x)の微分の形に、後ろ半分はf(x)の微分の形になっているのが分かると思います。なので
```math
\displaylines{
\lim_{h \to 0}\frac{f(x+h)(g(x+h)-g(x))}{h}+\frac{g(x)(f(x+h)-f(x))}{h}\\
\lim_{h \to 0}f(x+h)g'(x)+g(x)f'(x)\\
既に分母のhは亡くなっているため、h=0としてよいため\\
=f(x)g'(x)+f'(x)g(x)
}
```
となります。先ほどの具体例に当てはめてみましょう。
```math
\displaylines{
f'(x)=x^2+3x+3=2x+3\\
g'(x)=2x+5=2\\
f(x)g'(x)=2x^2+6x+6\\
f'(x)g(x)=4x^2+16x+15\\
(f(x)g(x))'=f(x)g'(x)+f'(x)g(x)=6x^2+22x+21
}
```
となっており、ひとまずは間違っていなさそうだという事が分かりますね。もうちょっとだけ頑張りましょう。
次は$$(\frac{1}{f(x)})'$$です。これもひとまず極限の式に当てはめてみます。

```math
\displaylines{
\lim_{h \to 0}\frac{\frac{1}{f(x+h)}-\frac{1}{f(x)}}{h}\\
\lim_{h \to 0}\frac{1}{h}\frac{f(x)-f(x+h)}{f(x+h)f(x)}\\
\lim_{h \to 0}\frac{f(x)-f(x+h)}{h}\frac{1}{f(x+h)f(x)}\\
\lim_{h \to 0}-\frac{f(x+h)-f(x)}{h}\frac{1}{f(x+h)f(x)}
}
```
う～ん、なんともややこしいですが左側の式がf(x)の微分の式になっていることがわかりますね？
そうすると
```math
\displaylines{
\lim_{h \to 0}-\frac{f(x+h)-f(x)}{h}\frac{1}{f(x+h)f(x)}\\
=\lim_{h \to 0}-f'(x)\frac{1}{f(x+h)f(x)}\\
=-f'(x)\frac{1}{f^2(x)}\\
=\frac{-f'(x)}{f^2(x)}\\
}
```
となります。ちなみにf(x)/g(x)やf(g(x))などは、この辺の考え方の応用でどうにかなるので、あとは各自で調べてみてください。
ここで慣れておいて欲しいのは関数を関数のまま扱う抽象思考です。

じゃあ実際ゲームプログラミングにおいてどのような面で役に立っているのかというと、とりあえず僕の知ってる限りでは、
+ ピクセルシェーダにおいて微小変化から法線(勾配)を求める事ができる
+ ベジェ曲線上の座標から媒介変数tの近似を求める事ができる(ニュートン法)

くらいですね。フーリエ変換の説明ができればその辺でも微分方程式とか使うのですが、ちょっとこの段階では早すぎるため上の二つについて解説します。

#### ピクセルシェーダの中で勾配を求める
ここで使うのは、高校の頃に習った微分とはちょっと違くて、**偏微分**とよばれるものなんですが、まぁ考え方は同じで、微小変化を何かに役立てようというものです。
さて、ピクセルシェーダ内で使用する微小変化というのは「テクスチャ」として渡された情報のX方向Y方向の変化量です。
また、その時の微小変化の単位はその時のテクスチャの解像度に依存しますので、渡されたテクスチャの1ピクセル分がhにあたると考えます。もちろん「無限小」に比べればはるかに大きな範囲になりますが、テクスチャにおいては変化の単位が1ピクセルより小さいことはないため、これで良しとします。

で、ピクセルシェーダの中でP(x,y,z)のそれぞれの座標が渡されるとします。これはワールド・ビュー行列までが計算された後だと考えます。
また、偏微分を求める関数としては、hlslの場合ddx,ddy関数(glslの場合はdFdx,dFdy関数)がありますので、それを使います。
```hlsl
float4 main(PInput input):SV_Target{
    float3 N = cross(ddx(input.pos),ddy(input.pos));
    return float4(normalize(N),1);
}
```
という感じで、もし法線情報が与えられてなくても、座標のxyの変化率から法線を求める事ができます。

#### ニュートン法で媒介変数tの近似値を求める
ニュートン法(ニュートン・ラフソン法とも)というのは、微分の考えによって、y=f(x)において「yはわかっているがxがわからない」そんな時に、微分によってなんとかかんとかxを求めようというものです。

どんな時に使うかって？例えば、色んなアニメーションツールにおいて、特定のパラメータの時間変化をベジェ曲線で設定することがあります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/ffef0619-90ad-13e3-edf2-43a7e5c43eb9.png)

このとき、xが時間方向で、yが実際のパラメータの変化を表しています。もしyが移動量だとするなら、この図の場合だと、最初はゆっくり、中間で急に早くなって、終盤直前ではまたゆっくりというような変化になります。

なお、ベジェ曲線というのは媒介変数tを用いてxとyが曲線的に変化するようにしたものです。図のベジェだと三次ベジェであるため、x=B(t)とy=B(t)の式はそれぞれ

```math
\displaylines{
x=B_x(t)=(1-t)^3x_0+3(1-t)^2tx_1+3(1-t)t^2x_2+t^3x_3\\
y=B_y(t)=(1-t)^3y_0+3(1-t)^2ty_1+3(1-t)t^2y_2+t^3y_3
}
```
と、まぁ、どっかで見たような式になっていますが、x₀, x₁, x₂, x₃はそれぞれコントロールポイントの座標値です。ともかくtが決まればxが決まる式になっています。

で、このtが先にわかっていれば問題ないんですが、「xを与えてyを知りたい」という要望の時に困っちゃいます。なぜならこの式においてはxとyの間には何の関連性もないからです。

それを行うためには、x=B(t)の逆関数を求めて、xからtを求め、そこで得られたtからyを求めるしかないわけです。ところが逆関数を作るのは大変です。

そこで、ニュートン法を用いてx座標からtの近似値を求め、tさえわかれば、あとはyを求めるのは容易であるということです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/c6fe601e-3cc1-5ae5-ea7f-d81d77014960.png)


まず、このようなグラフの交点を求めたいとします。横軸をtで縦軸をxとします。グラフがt軸と交点を持っていますが、これはx=0固定する必要はなく、必要に応じて定数を加減算すればこの形にできますので、ひとまずこの交点のtを求めたいと考えてください。

ともかく図においてはx=0は分かってるため、この時のtを求めればいいだけです。でもわかりません。グラフの形だけは分かっています。

そこで「確定で分かっている座標」をまず取ってきます。x=B(t)はわかっているため、適当な座標のt₀,B(t₀)は得られます。また、B'(t)の式も分りますね。

```math
\displaylines{
B_x'(t)=-(1 -3t+3t^2-t^3)'x_0+3(t-2t^2+t^3)'x_1+3(t^2-t^3)'x_2+(t^3)'x_3\\
=(-3+6t-3t^2)x_0 + (3-12t+9t^2)x_1 +(6t-9t^2)x_2+3t^2x_3
}
```
ともかく、これが微分ですが、前にも説明した通り、微分というのは接線の傾きを求めるもので、このB'(t)は傾きです。じゃあその適当な点P₀(t₀,B(t₀))における接線の式は

```math
x-B_x(t_0)=B_x'(t_0)(t-t_0)
```
で、この接線がt軸(x=0)との交点を持つとします。そうしたうえでt=の式に変換してみます

```math
\displaylines{
0-B_x(t_0)=B_x'(t_0)(t-t_0)\\
\frac{-B_x(t_0)}{B_x'(t_0)}=(t-t_0)\\
t=t_0-\frac{B_x(t_0)}{B_x'(t_0)}
}

```
となります。要はこれで、接線がt軸と交差するときのtの値が求まりました。で、なんでこんなけったいな事をやったのかというと、これを繰り返すことで、下の図のように求めたい真の値に近づくからです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/2dab5644-ba87-f0fa-4c69-6569c46a8f4d.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/e9e12285-1869-b037-0fd6-f85e0bc55f18.png)

まぁ、あくまでも近似値なので、場合によっては下図のようにオーバーしてしまう事もありますが

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275040/510bff50-3526-7374-5399-7b5cbf01eaf8.png)
何度か繰り返していれば、接線の性質上求めたい交点に近づいていくのが分かります。

これをプログラムに直すと
```C++
float t = x;//ひとまず適当なtをおく(ここでは=xとする)
float k0 = 1 + 3 * a.x - 3 * b.x;//t^3の係数
float k1 = 3 * b.x - 6 * a.x;//t^2の係数
float k2 = 3 * a.x;//tの係数
//誤差の範囲内かどうかに使用する定数
constexpr float epsilon = 0.0005f;
//ニュートン法
for (int i = 0; i < n; ++i) {
	//f(t)求めまーす
	auto ft = k0 * t*t*t + k1 * t*t + k2 * t - x;
	//もし結果が0に近い(誤差の範囲内)なら打ち切り
	if (ft <= epsilon && ft >= -epsilon)break;
	//f'(t)求めまーす
	auto fdt = 3 * k0*t*t + 2 * k1*t + k2;//微分結果
	if (fdt == 0)break;
	t -= ft / fdt;
}
```

でも、ここでの目的はx座標からy座標を求める事だったため、得られたtを元のベジェ式に入れなおしてyを得ます。

```C++
///ニュートン法
float GetYFromXOnBezierNT(float x,  Position2 a, Position2 b, unsigned int n = 8) {
	if (a.x == a.y&&b.x == b.y)return x;//計算不要
	float t = x;
	float k0 = 1 + 3 * a.x - 3 * b.x;//t^3の係数
	float k1 = 3 * b.x - 6 * a.x;//t^2の係数
	float k2 = 3 * a.x;//tの係数
	//誤差の範囲内かどうかに使用する定数
	constexpr float epsilon = 0.0005f;
	//ニュートン法
	for (int i = 0; i < n; ++i) {
		//f(t)求めまーす
		auto ft = k0 * t*t*t + k1 * t*t + k2 * t - x;
		//もし結果が0に近い(誤差の範囲内)なら打ち切り
		if (ft <= epsilon && ft >= -epsilon)break;
		//f'(t)求めまーす
		auto fdt = 3 * k0*t*t + 2 * k1*t + k2;//微分結果
		if (fdt == 0)break;
		t -= ft / fdt;
	}
	//既に求めたいtは求めているのでyを計算する
	auto r = 1 - t;
	return t * t*t + 3 * t*t*r*b.y + 3 * t*r*r*a.y;
}
```
はい、このように微分を応用してベジェのｘ座標からy座標を求めることができました。
実例映像を添付しますね。

https://x.com/CTsuchinoko/status/1802654422557679796

一応今回のベジェに参考になるソースコードは
https://github.com/boxerprogrammer/Mathmatics/tree/master/Bezier_3/BezierTime
に上げていますので、試してみてください

今回は積分までやる予定でしたが、ここまでで十分に長くなってしまったため、一旦ここで区切ります。
お疲れさまでした。