shaderで内積、外積の使いどころ

#はじめに
内積はA⋅B、外積はA×Bと表記する。この演算子の形から`dot()`,`cross()`と関数名がきたようだ。`cross()`に関しては性質も表していて、良いネーミングだと、いつも思ってます。
数学の知識は余り無いので、ほぼ感覚で使いこなす説明になると思います。
#２Dでの内積、外積
glslには２Dの外積がない。理由はよくわからない。なので自前で用意しなければならない。

```glsl
float cross(vec2 a, vec2 b) 
{
  return a.x * b.y - a.y * b.x;
}
```
以後のスクリプトには、この`cross()`を使っていきます。参考までに`dot()`を関数で書いてみます。

```glsl
float dot(vec2 a, vec2 b) 
{
  return a.x * b.x + a.y * b.y;
}
```
内積、外積の数式を見ていきます。

``` 
A⋅B = |A||B|cos(Θ)
A×B = |A||B|sin(Θ)
```
以上の情報だけで、話を進めていきます。
shaderでグラフィックをするなら距離関数は避けては通れません。そして、これが色々な事を教えてくれます。
なので`cross()`の性質をみるために距離関数を書いてみます。

```glsl
vec2 p = (gl_FragCoord.xy * 2.0 - resolution) / resolution.y;
float a = radians(30.0);
vec2 v = vec2(cos(a), sin(a));
float d = abs(cross(p, v));
gl_FragColor = vec4(d, d, d, 1.0);
```
![ダウンロード (2).png](https://qiita-image-store.s3.amazonaws.com/0/153571/47bf8685-f42f-c060-8028-0d0bfabec6e9.png)



これは、原点を通るベクトルvの直線の距離関数になります。
では

```glsl
float d = abs(cross(p, v));
```
を

```glsl
float d = abs(dot(p, v));
```
に変えると

![ダウンロード (1).png](https://qiita-image-store.s3.amazonaws.com/0/153571/14e48772-478d-6d18-1043-a12b8b534b40.png)

先程の直線と直交する直線の距離関数になります。
改めて先程の`cross()`と`dot()`の関数を良くみてください。

```glsl
float cross(vec2 a, vec2 b) 
{
  b = vec2(-b.y, b.x); // 直交するベクトルに変換
  return dot(a,b);
}
```
こういう構造になってます。

#２D座標の回転
この性質を使うと座標を回転させられます。

```glsl
float angle = radians(30.0);
vec2 v = vec2(cos(a), sin(a));
p = vec2(dot(p,v), cross(p,v));
```
`cross()`を使わずに書くと

```glsl
float angle = radians(30.0);
vec2 v1 = vec2(cos(a), sin(a));
vec2 v2 = vec2(-sin(a), cos(a)); // v1の直交ベクトル
p = vec2(dot(p,v1), dot(p,v2));
```
これを見ると２Dの回転行列の仕組みがみえてきます。

```glsl
float angle = radians(30.0);
vec2 v1 = vec2(cos(a), sin(a));
vec2 v2 = vec2(-sin(a), cos(a)); // v1の直交ベクトル
p = mat2(v1, v2) * p;
//p = mat2(cos(a), sin(a), -sin(a), cos(a)) * p;
```
Don't think! Feel. by Bruce Lee

# dot()の使いどころ
`dot()`の用法として、長さの比較に使われたりしてます。
ベクトルの長さの取得は`length()`を使ってますが、別の書き方をすると

```glsl
sqrt(dot(v,v))
つまり
sqrt(v.x * v.x + v.y * v.y) // 三平方の定理
```
単純に比較するだけなら`sqrt()`は無駄な処理になります。`dot(v,v)`だけで比較します。
iq氏のdistans functionsでも使ってます。

```glsl
float dot2( in vec3 v ) { return dot(v,v); }
float udTriangle( vec3 p, vec3 a, vec3 b, vec3 c )
{
    vec3 ba = b - a; vec3 pa = p - a;
    vec3 cb = c - b; vec3 pb = p - b;
    vec3 ac = a - c; vec3 pc = p - c;
    vec3 nor = cross( ba, ac );

    return sqrt(
    (sign(dot(cross(ba,nor),pa)) +
     sign(dot(cross(cb,nor),pb)) +
     sign(dot(cross(ac,nor),pc))<2.0)
     ?
     min( min(
     dot2(ba*clamp(dot(ba,pa)/dot2(ba),0.0,1.0)-pa),
     dot2(cb*clamp(dot(cb,pb)/dot2(cb),0.0,1.0)-pb) ),
     dot2(ac*clamp(dot(ac,pc)/dot2(ac),0.0,1.0)-pc) )
     :
     dot(nor,pa)*dot(nor,pa)/dot2(nor) );
}
```
同じくiq氏のdistans functionsの中の

```glsl
float sdCapsule( vec3 p, vec3 a, vec3 b, float r )
{
    vec3 pa = p - a, ba = b - a;
    float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
    return length( pa - ba*h ) - r;
}
```
これの

```glsl
dot(pa,ba)/dot(ba,ba)
```
この部分。数式にすれば

```
|A||B|cos(theta) / |A||A|cos(0)
|B|cos(theta) / |A|
```
これで比率を取り出し、パラメトリック関数のt値にして位置情報を所得する手法。

内積の式を見てもらえばわかるけど、2つのベクトルの角度を取得する時に

```glsl
acos(dot(normalize(a),normalize(b)))
```
と`dot()`を利用できます。

他は`fract(sin())`でノイズを作ったりする時に、`vec2`とか`vec3`とかを適当な`float`にしたい時に使ったりしています。

# ３Dのcross()について
`dot()`については２Dの延長の考え方でいけるけど、`cross()`は別物です。とりあえず、イメージで説明します。アナログ時計の時針と分針を２つのベクトルとすれば、２つのベクトルの`cross()`は軸棒になります。軸棒と時針は直交してます。軸棒と分針は直交してます。これだけのルールですね。この簡単なルールが威力絶大で、前後左右上下のわからない３D空間においては、`cross()`だけが頼りです。`cross()`を理解してからは、日常に重力と言う一本のベクトルが常にある、ありがたみがわかりました。それは置いておいて`cross()`の性質を探ってみます。既にある関数ですが、`cross()`の関数を書いてみます。

```glsl
vec3 cross(vec3 a, vec3 b) 
{
  return vec3(
    a.y * b.z - b.y * a.z,
    a.z * b.x - b.z * a.x,
    a.x * b.y - b.x * a.y);
}
```
こんな中身をしてます。
２Dの`cross()`を、この`vec3`の`cross()`を使って書いてみます。

```glsl
float cross(vec2 a, vec2 b) 
{
  return cross(vec3(a,0),vec3(b,0)).z;
}
```
これを書くことによって見えてくるものがあります。
`cross()`に入力する２つのベクトルがxy平面上にあるなら、`cross()`の出力ベクトルにはｚ要素しか存在しません。２つのベクトルの角度が変わると、`cross()`の出力ベクトルのスカラーは変化します。ここで[いつか読んだtwitterを思い出しました。](https://twitter.com/gam0022/status/846405574283616256)
`acos(dot(normalize(a),normalize(b)))`で2つのベクトルの角度が取得出来るので、試してはないけど`cross()`のスカラーを`asin()`に入力すれは角度が取得できる気がします。でも`acos(dot())`の方がいいかな。
２つのベクトルが重なると`cross()`のスカラーは０になります。普段オイラー角を使って回転させないのでジンバルロックを経験したことが無いけど、たぶんこれが原因かなと思ってます。やはり１つのベクトルだと直交するベクトルは無限個になっちゃうから、２つのベクトルが必要なのですね。なんか少し哲学が入りそうです。

#３D座標の回転
`cross()`は３D座標の回転時に自分の姿勢をkeepすることが、主な使われ方だと思います。ベクトルを使ってオブジェクトの向きのコントロールとかにも使えます。今回はraymarchigで使えるLookAt行列を作る方法の説明をします。
upというベクトルが必要なのは`cross()`には２つのベクトルが必要だからです。upには、とりあえず重力方向の`vec3(0,1,0)`を使ってください。

```glsl
mat3 lookat(vec3 eye, vec3 target, vec3 up)
{
	vec3 w = normalize(target-eye);
	vec3 u = normalize(cross(w,up));
	vec3 v = normalize(cross(u,w)); //w,uは直交してるのでnormalize()無しでも良いと思うけど、一応。
    return mat3(u,v,w);
}
```
すでに２Dの回転行列を書いておいたから、雰囲気はわかると思います。


```glsl
vec3 w = normalize(target-eye);
vec3 u = normalize(cross(w,up));
vec3 v = normalize(cross(u,w));
p = p.x * u + p.y * v + p.z * w;
```
こう使って回転させる事もできます。
視点から目標への正規ベクトルがwになります。wと補助のベクトルupとのcross()で画面水平方向の正規ベクトルuを取得。upはwと直交していないので、改めてu,wの`cross()`を取って画面垂直方向の正規ベクトルvを取得。それを`mat3`に入力。という手順です。雰囲気で言えばuvwはxyzにあたるって感じですかね。

####ここに差し込む形で追記
```glsl
vec3 w = normalize(target-eye);
vec3 u = normalize(cross(w,up));
vec3 v = normalize(cross(u,w));
p = vec3(dot(p,u), dot(p,v), dot(p,w));
```
今まで勘違いして、これも同等と言っていたが違いました。
`inverse(mat3(u,v,w))*p`と同等でした。
2Dの方も直さないといけないけれど、説明が絡んでるので、その内に直します。使うには問題ないです。３Dは理解しておかないとヒドイ目に会うかな。
(2020/10/06)


Don't think! Feel. by Bruce Lee

上で書いた理由でwとupは違うベクトルでなければいけません。upは必ず`vec3(0,1,0)`でなくては良いので、自由に変えれます。真下を見たいならupの`vec3(0,1,0)`と重なってしまうので`vec3(0,0,1)`とかにします。upを`sin()`とかを使って揺らすのもありですけど、３D酔いに注意。

#最後に
ここまでくれば、iq氏のdistans functionsの平面の距離関数が読み解けるんじゃないですかね。
法線と座標のdot()で距離ってことですね。

```glsl
float sdPlane( vec3 p, vec4 n )
{
  // n must be normalized
  return dot(p,n.xyz) + n.w;
}
```
#追記
ここの応用です。
あるベクトルvのCylinderの距離関数。座標をp、ベクトルをv、通過点をc、半径をtすれば

```glsl
float deCylinder(vec3 p, vec3 v, vec3 c, float t)
{
	return length(cross(p - c, normalize(v))) - t;
}
```
これもイケます。
正規化されたベクトル vと座標 pの距離を求める需要もあるみたいですね。通過点をcとした場合の距離は

```glsl
    float d = length(cross(p - c, v));
```
同じ事ですけど、一応書いておきます。距離関数をブラックボックスで扱っていると、こういう事が意外とスルーされちゃったりします。
***

法線nの平面にベクトルaを投影させるには、

```glsl
cross(cross(n,a),n)
```
もしくは、

```glsl
a-dot(n,a)*n
```
任意軸の折りたたみは

```glsl
p = p - 2.0 * min(0.0, dot(p, n)) * n // pは座標、nは法線
```
なので、この投影の式に似ている。回転に関する関数を作るのにdot()は、色々な場所で登場する気がしてます。

***
投影の方法を眺めていたら、任意軸での回転を思いつきました。
任意軸で回転することは、任意軸に垂直な平面に投影して、その面で回転して、元の高さに戻す事。これを関数にしてみました。
平面の回転の関数は、

```glsl
vec2 rotate(vec2 p, float theta)
{
    return p * cos(theta) + vec2(-p.y, p.x) * sin(theta);
}
```
こう書けます。座標(ベクトル)と`cos(theta)`を乗算したものと、直交座標(ベクトル)と`sin(theta)`を乗算したものを加算します。pを`vec3(1,0)`、pの直交座標を`vec3(0,1)`と考えれば、見たことのある`vec2(cos(theta),sin(theta))`と同じになるのでイメージしやすいと思います。投影面に対して、これを使います。
投影したベクトルに`cos(theta)`を乗算して、投影したベクトルに投影面上で直交するベクトルに`sin(theta)`を乗算して、それらを加算。回転軸方向に高さを戻す。

```glsl
vec3 rotate(vec3 p,vec3 axis,float theta)
{
    axis=normalize(axis);
    vec3 v = cross(axis,p), u = cross(v, axis);
    return u * cos(theta) + v * sin(theta) + axis * dot(p, axis);   
}
```
単体で使うなら、面倒な行列やクォータニオンから解放されます。

x軸、y軸、z軸で回すのを2Dでする方法があるのだけど、意外と浸透していない気がするので、ちょっと書いておきます。上の2Dの`rotate()`を使っています。別に`mat2(cos(a),sin(a),-sin(a),cos(a))`とかで回しても問題ないです。好みの問題です。

```glsl
p.yz = rotate(p.yz, 1.0); // x軸回し
p.zx = rotate(p.zx, 1.0); // y軸回し
p.xy = rotate(p.xy, 1.0); // z軸回し
```
この方法だと、xyzそれぞれの軸用の回転関数を用意しなくて良いので、お気楽です。

***
[iq氏の距離関数](http://iquilezles.org/www/articles/distfunctions/distfunctions.htm)の中にあるVertical Capsuleを垂直方向限定ではなく、任意方向にするのに`dot()`が使えます。

```glsl
float sdVerticalCapsule( vec3 p, float h, float r )
{
    p.y -= clamp( p.y, 0.0, h );
    return length( p ) - r;
}
```
これが元の距離関数。任意方向をvec3 v とした場合、(vはnormalize()必須)

```glsl
float deCapsule( vec3 p, vec3 v, float h, float r )
{
    p -= clamp(dot(p,v), 0.0, h) * v;
    return length( p ) - r;
}    
```
となります。`clamp()`を使って空間を間引く感じですかね。面白い使い方です。
***
どうも任意方向のベクトルの操作みたいなのがdot()の主な使われ方のような気がしてきました。任意軸の折りたたみでもdot()を使ってるので、ちょっと解説してみます。
一番単純な折りたたみのabs()を使います。

```glsl
p.x = abs(p.x);
```
こんな感じで使います。この折りたたみを違う書き方をします。

```glsl
p.x = p.x-2.0*min(0.0,p.x); 
```
ちょっと見に理解しづらいけど、正の数、負の数と分けてみれば理解できると思います。min()を使ってるのでp.xが正の数なら 0 に変換されます。つまりp.xが正の数の時は、何もしない。負の数の時は

```glsl
p.x = p.x - 2.0 * p.x;
つまり
p.x = -p.x;
```
なので、負の数は正の数に変換されて、`abs()`と同等になります。
後は上と同じように任意方向に対応するように改造。

```glsl
p = p - 2.0 * min(0.0, dot(p, n)) * n // pは座標、nは法線
```
復習を兼ねて、もう一度書くと、`dot(p,n)`でスカラーが取得できて、単位ベクトルn に乗算してるだけ。
Don't think! Feel. by Bruce Lee
これには、`min()`が使われてます。これは非常にラッキーな事でsmooth minimum の関数(世間じゃあsmin()とか言う名前で関数名を付けている奴)をmin()の代わりに使ってやると面白い絵になります。
[参考：任意方向の折りたたみ](http://blog.hvidtfeldts.net/index.php/2011/08/distance-estimated-3d-fractals-iii-folding-space/)

***
直交行列ってものもあるみたいですね。こいつを使っての回転もできます。単位行列は `mat2(1,0,0,1)` 。なんとなく虚数が回転と関係があるって感じさせる `-1`だね。 

```glsl
p = vec2(dot(p,v),dot(p,mat2(0,1,-1,0)*v));
```
***
ベクトルをどうこう言うスキルとは、ちょっと離れちゃうけど、ユニークなseedを作るのに使えるスキル。
例えば

```glsl
p = mod(q, 4.0) - 2.0;
```
と折りたたむ時に、各マス目にユニークなseedを付けたい事があります。そんな時の一つのアイディアとして

```glsl
vec3 u = floor(p / 4.0);
float seed = u.x * 1.0 + u.y * 10.0 + u.z * 100.0;
```
という方法があります。これを`dot()`を使って書き換えると

```glsl
float seed = dot(u, vec3(1.0, 10.0, 100.0));
```
まあ、これだけのことですけど、知っていると色々と応用が効きます。 `vec3(1.0, 10.0, 100.0)` はxyzの各要素が同じ数字じゃなきゃ問題無いと思います。 `fract(sin())` でnoiseを作っても、数字の組み合わせによっては変な偏りもでる時があるけど、そこは適当に調整してます。uが`vec3(0)`に成るのを嫌うなら、適当な数字を加算しても良いと思います。
`vec3`のseedが欲しければ、

```glsl
vec3 seed = vec3(
    dot(u, vec3(1.0, 10.0, 100.0),
    dot(u, vec3(123.0, 456.0, 789.0),
    dot(u, vec3(0.123, 0.456, 0.789))
```
で、このseedを `fract(sin())` に入れてやるか、数字の組み合わせによっては `fract()` だけで足りるかも？
***
オブジェクトを、ある点に向かせるのに前述のlookAtでは、任意軸回転になるため動きとして無理がでてくる。そこでオイラー角を使ってみた。もっとシンプルにしたかったけど、難しくて、これが限界です。とりあえず関数として書いておきます。pが座標、dがz軸が法線となる面が向く点。

```glsl
vec3 pointAt(vec3 p, vec3 d)
{
	vec3 u,v,w;
	w = normalize(cross(cross(vec3(0,1,0),d),vec3(0,1,0)));
	u = normalize(cross(w,vec3(0,1,0)));
	v = vec3(0,1,0);//normalize(cross(u,w));
	p = mat3(u,v,w)*p;
	w = normalize(vec3(0, d.y,length(d.xz)));
	u = vec3(1,0,0);//normalize(cross(vec3(0,1,0),w));
	v = normalize(cross(w,u));
	return mat3(u,v,w)*p;
}
```
短く詰めると

```glsl
vec3 pointAt(vec3 p, vec3 d)
{
	vec3 a = normalize(cross(cross(vec3(0,1,0),d),vec3(0,1,0)));
	vec3 b = normalize(vec3(0, d.y,length(d.xz)));
	return mat3(vec3(1,0,0),normalize(cross(b,vec3(1,0,0))),b)*
			mat3(normalize(cross(a,vec3(0,1,0))),vec3(0,1,0),a)*p;
}
```
２つある行列の乗算する順番を入れ替えると、任意軸回転と同じ動きになります。その辺りがシンプルにするヒントかもしてません。
今は文章での説明は、できそうもありません。もう少し理解が進んだら追記しておきます。

#####追記の追記
乗算する順番を入れ替えるがヒントでした。結論から言うと逆行列を使うが正解でした。相変わらず説明が出来ません。経験則です。
改めて関数を書いてみます。pが座標、dがz軸が法線となる面が向く点、upが補助ベクトル。

``` glsl
vec3 pointAt(vec3 p, vec3 d, vec3 up)
{
	vec3 u = normalize(cross(d,up));
	vec3 v = normalize(cross(u,d)); //u,dは直交してるのでnormalize()無しでも良いと思うけど、一応。
	return p * mat3(u,v,d);
}
```
この関数はオブジェクトがマウスポインターの方を向くとかに使えるのじゃないですかね。
以前の回りくどい関数は消去しようと思ったけど、なにかのヒントになる事を期待して消さないでおきます。
`p * mat3(u,v,d)` は通常の乗算する順番と逆になってます。これで逆行列を乗算したと同等になります。
参考までに他の逆行列の使い方を書いておきます。

```glsl
    p = transpose(mat3(u,v,w)) * p; // 試した事ないけどwebGL1.0でも使えるかな。
    p = inverse(mat3(u,v,w)) * p; // #version 300 es より使えます。
    p = p * mat3(u,v,w);
    p *= mat3(u,v,w);
```
`transpose()`は回転行列の時は逆行列として使えると、どこかで読みました。理由は、わかりません。
追記 :　回転行列の転置が逆行列となる理由が書いてあるtwitterが出てきました。
https://twitter.com/gam0022/status/1153643512644063234
この部分のサンプルを書きました。
https://gaz.hateblo.jp/entry/2019/05/12/164845

***
ここまで、わかってくると、任意軸回転のmat3()を導き出したくなります。そのうちにチャレンジしてみます。
***
追記を書いてきて見えてきたのは、オブジェクトの向きの操作は、cross()を多用し3直角のベクトルを特定しmat3に入力して、それに今の座標を乗算すればいい。解れば単純な話でした。なら、わざと3直角を崩してエフェクトってのもコントロールしやすくなりそうです。
右手系、左手系の回転とか、あるみたいですけど、その辺りは適当に書いているのでよくわかりません。`cross()`の引数に順番を入れ替えるだけなので、適宜調整してください。うかっかりすれば上下反対に関数があるかもしれませんが、それも調整してください。
グラフィックプログラミングの良いところは、正誤が素人でも解る事。大した数学の能力も無いけど、なんとかなるものですね。