はじめに
内積はA⋅B、外積はA×Bと表記する。この演算子の形からdot(),cross()と関数名がきたようだ。cross()に関しては性質も表していて、良いネーミングだと、いつも思ってます。
数学の知識は余り無いので、ほぼ感覚で使いこなす説明になると思います。
２Dでの内積、外積
glslには２Dの外積がない。理由はよくわからない。なので自前で用意しなければならない。
float cross(vec2 a, vec2 b)
{
return a.x * b.y - a.y * b.x;
}
以後のスクリプトには、このcross()を使っていきます。参考までにdot()を関数で書いてみます。
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()の性質をみるために距離関数を書いてみます。
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);
これは、原点を通るベクトルvの直線の距離関数になります。
では
float d = abs(cross(p, v));
を
float d = abs(dot(p, v));
に変えると
先程の直線と直交する直線の距離関数になります。
改めて先程のcross()とdot()の関数を良くみてください。
float cross(vec2 a, vec2 b)
{
b = vec2(-b.y, b.x); // 直交するベクトルに変換
return dot(a,b);
}
こういう構造になってます。
２D座標の回転
この性質を使うと座標を回転させられます。
float angle = radians(30.0);
vec2 v = vec2(cos(a), sin(a));
p = vec2(dot(p,v), cross(p,v));
cross()を使わずに書くと
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の回転行列の仕組みがみえてきます。
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()を使ってますが、別の書き方をすると
sqrt(dot(v,v))
つまり
sqrt(v.x * v.x + v.y * v.y) // 三平方の定理
単純に比較するだけならsqrt()は無駄な処理になります。dot(v,v)だけで比較します。
iq氏のdistans functionsでも使ってます。
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の中の
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;
}
これの
dot(pa,ba)/dot(ba,ba)
この部分。数式にすれば
|A||B|cos(theta) / |A||A|cos(0)
|B|cos(theta) / |A|
これで比率を取り出し、パラメトリック関数のt値にして位置情報を所得する手法。
内積の式を見てもらえばわかるけど、2つのベクトルの角度を取得する時に
acos(dot(normalize(a),normalize(b)))
とdot()を利用できます。
他はfract(sin())でノイズを作ったりする時に、vec2とかvec3とかを適当なfloatにしたい時に使ったりしています。
３Dのcross()について
dot()については２Dの延長の考え方でいけるけど、cross()は別物です。とりあえず、イメージで説明します。アナログ時計の時針と分針を２つのベクトルとすれば、２つのベクトルのcross()は軸棒になります。軸棒と時針は直交してます。軸棒と分針は直交してます。これだけのルールですね。この簡単なルールが威力絶大で、前後左右上下のわからない３D空間においては、cross()だけが頼りです。cross()を理解してからは、日常に重力と言う一本のベクトルが常にある、ありがたみがわかりました。それは置いておいてcross()の性質を探ってみます。既にある関数ですが、cross()の関数を書いてみます。
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()を使って書いてみます。
float cross(vec2 a, vec2 b)
{
return cross(vec3(a,0),vec3(b,0)).z;
}
これを書くことによって見えてくるものがあります。
cross()に入力する２つのベクトルがxy平面上にあるなら、cross()の出力ベクトルにはｚ要素しか存在しません。２つのベクトルの角度が変わると、cross()の出力ベクトルのスカラーは変化します。ここでいつか読んだtwitterを思い出しました。
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)を使ってください。
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の回転行列を書いておいたから、雰囲気はわかると思います。
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にあたるって感じですかね。
ここに差し込む形で追記
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()で距離ってことですね。
float sdPlane( vec3 p, vec4 n )
{
// n must be normalized
return dot(p,n.xyz) + n.w;
}
追記
ここの応用です。
あるベクトルvのCylinderの距離関数。座標をp、ベクトルをv、通過点をc、半径をtすれば
float deCylinder(vec3 p, vec3 v, vec3 c, float t)
{
return length(cross(p - c, normalize(v))) - t;
}
これもイケます。
正規化されたベクトル vと座標 pの距離を求める需要もあるみたいですね。通過点をcとした場合の距離は
float d = length(cross(p - c, v));
同じ事ですけど、一応書いておきます。距離関数をブラックボックスで扱っていると、こういう事が意外とスルーされちゃったりします。
法線nの平面にベクトルaを投影させるには、
cross(cross(n,a),n)
もしくは、
a-dot(n,a)*n
任意軸の折りたたみは
p = p - 2.0 * min(0.0, dot(p, n)) * n // pは座標、nは法線
なので、この投影の式に似ている。回転に関する関数を作るのにdot()は、色々な場所で登場する気がしてます。
投影の方法を眺めていたら、任意軸での回転を思いつきました。
任意軸で回転することは、任意軸に垂直な平面に投影して、その面で回転して、元の高さに戻す事。これを関数にしてみました。
平面の回転の関数は、
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)を乗算して、それらを加算。回転軸方向に高さを戻す。
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))とかで回しても問題ないです。好みの問題です。
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氏の距離関数の中にあるVertical Capsuleを垂直方向限定ではなく、任意方向にするのにdot()が使えます。
float sdVerticalCapsule( vec3 p, float h, float r )
{
p.y -= clamp( p.y, 0.0, h );
return length( p ) - r;
}
これが元の距離関数。任意方向をvec3 v とした場合、(vはnormalize()必須)
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()を使います。
p.x = abs(p.x);
こんな感じで使います。この折りたたみを違う書き方をします。
p.x = p.x-2.0*min(0.0,p.x);
ちょっと見に理解しづらいけど、正の数、負の数と分けてみれば理解できると思います。min()を使ってるのでp.xが正の数なら 0 に変換されます。つまりp.xが正の数の時は、何もしない。負の数の時は
p.x = p.x - 2.0 * p.x;
つまり
p.x = -p.x;
なので、負の数は正の数に変換されて、abs()と同等になります。
後は上と同じように任意方向に対応するように改造。
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()の代わりに使ってやると面白い絵になります。
参考：任意方向の折りたたみ
直交行列ってものもあるみたいですね。こいつを使っての回転もできます。単位行列は mat2(1,0,0,1) 。なんとなく虚数が回転と関係があるって感じさせる -1だね。
p = vec2(dot(p,v),dot(p,mat2(0,1,-1,0)*v));
ベクトルをどうこう言うスキルとは、ちょっと離れちゃうけど、ユニークなseedを作るのに使えるスキル。
例えば
p = mod(q, 4.0) - 2.0;
と折りたたむ時に、各マス目にユニークなseedを付けたい事があります。そんな時の一つのアイディアとして
vec3 u = floor(p / 4.0);
float seed = u.x * 1.0 + u.y * 10.0 + u.z * 100.0;
という方法があります。これをdot()を使って書き換えると
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が欲しければ、
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軸が法線となる面が向く点。
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;
}
短く詰めると
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が補助ベクトル。
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) は通常の乗算する順番と逆になってます。これで逆行列を乗算したと同等になります。
参考までに他の逆行列の使い方を書いておきます。
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()の引数に順番を入れ替えるだけなので、適宜調整してください。うかっかりすれば上下反対に関数があるかもしれませんが、それも調整してください。
グラフィックプログラミングの良いところは、正誤が素人でも解る事。大した数学の能力も無いけど、なんとかなるものですね。