こんにちは。traP Advent Calendar 2016 22日目担当のArkです。
今回はクリエイティブコーディングについて書きます。

1. クリエイティブコーディングとは

みなさん、クリエイティブコーディングを知っていますか？


僕は知りません。


・・・。


ちょうどいい言葉ないかなーとググってみたらこのワードが出てきたのでタイトルに使ってみました。
wikipediaによれば、

Creative coding is a type of computer programming in which the goal is to create something expressive instead of something functional.

とのこと。

要は「プログラミングによってクリエイティブな作品を創ろう」というものだと思います（たぶん）。
今年は クリエイティブコーディング Advent Calendar 2016 もやってるみたいですし、ブームなんでしょうか。

面白そうなので、なにか創ってみましょう！

2. この記事の内容

今回はこんなのを作ってみます。（「▷」ボタンを押すと動きます）
円と星が交互に変形するという単純なデモです。

順を追ってこのアニメーションの作り方を説明します。

3. 開発環境

今回は、Shadertoyを使います。（※注意: ページが重たいのでスマートフォンなどで開かないでください。）
ShaderとかShadertoyってなんぞという方は、sobayaが書いたこちらの記事を参照ください。
Shaderというのは†お絵かきツール†です。

4. 円を描く

// 中心vec2(0.0, 0.0), 半径rの円からの距離を返す
float circle(vec2 p, float r) {
    return length(p) - r;
}

vec3 calc(vec2 p) {
    float d = circle(p, 0.4);
    return vec3(pow(clamp(1.0-d, 0.0, 1.0), 5.0));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord ) {
	vec2 uv = (2.0*fragCoord.xy-iResolution.xy) / min(iResolution.x, iResolution.y);
	fragColor = vec4(calc(uv),1.0);
}

これで円が描けます。float circle(vec2 p, float r)は、半径$r$の円と点$p$との距離を返します。

図に表すとこんな感じ。
このように、与えられた点に対し、その点との距離を返す関数をdistance functionと言います。（そのまんまですね。）
distance functionは、数学で定義される距離関数とは違って負の値も許されます。
このcircle関数は、円の内側だと負の値が返されるようになります。

円はできたのですが、せっかくなので楕円まで拡張させちゃいましょう。

// 中心がvec2(0.0, 0.0)でr.x, r.yをそれぞれの軸の半径とする楕円との距離を返す
float ellipse(vec2 p, vec2 r) {
    return (length(p/r) - 1.0) * min(r.x, r.y);
}

5. 星を描く

次は星（五芒星の内側を塗りつぶしたもの）を描きます。・・・が、星のdistance functionをつくるのは一筋縄では行きません。
段階的に構成していきます。
（もっといい構成方法があるかもしれないですが、とりあえず今回はこの方法でやっていきます。）

5.1. 線分のdistance functionをつくる

float crs(vec2 v1, vec2 v2) {
    return v1.x*v2.y - v1.y*v2.x;
}

// 2点v1, v2を端点に持つ線分との距離を返す
float line(vec2 p, vec2 v1, vec2 v2) {
    p  -= v1;
    vec2 v = v2-v1;
    float t = dot(p, normalize(v));
    if (t<0.0) {
        return length(p);
    } else if (t>length(v)) {
        return length(p-v);
    } else {
        return abs(crs(p, normalize(v)));
    }
}

これが線分のdistance functionです。
細かい説明は省きますが$p$の位置を3つの場合に分けて距離を計算しています。

5.2. 五芒星のdistance functionをつくる

次は五芒星です。

• $\left(r\cos (0\cdot \frac{2\pi }{5}),r\sin (0\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (2\cdot \frac{2\pi }{5}),r\sin (2\cdot \frac{2\pi }{5})\right)$
• $\left(r\cos (1\cdot \frac{2\pi }{5}),r\sin (1\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (3\cdot \frac{2\pi }{5}),r\sin (3\cdot \frac{2\pi }{5})\right)$
• $\left(r\cos (2\cdot \frac{2\pi }{5}),r\sin (2\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (4\cdot \frac{2\pi }{5}),r\sin (4\cdot \frac{2\pi }{5})\right)$
• $\left(r\cos (3\cdot \frac{2\pi }{5}),r\sin (3\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (0\cdot \frac{2\pi }{5}),r\sin (0\cdot \frac{2\pi }{5})\right)$
• $\left(r\cos (4\cdot \frac{2\pi }{5}),r\sin (4\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (1\cdot \frac{2\pi }{5}),r\sin (1\cdot \frac{2\pi }{5})\right)$

五芒星は、これらの点を結んだ5本の線分で構成されます。
つまり、この5本の線分との距離のうち、最も小さい距離が五芒星との距離です。

線分との距離は[5.1.]でつくったline関数がそのまま使えます。
よってコードにすると次のようにかけます。

#define PI 3.14159265359
#define INF 1e10

// 中心vec2(0.0, 0.0)で半径rの五芒星との距離を返す
float pentagram(vec2 p, float r) {
    float d = INF;
    for (int i=0; i<5; i++) {
        float rad1 = 2.0*PI*float(i)/5.0;
        float rad2 = 2.0*PI*float(i+2)/5.0;
        vec2 v1 = vec2(cos(rad1), sin(rad1)) * r;
        vec2 v2 = vec2(cos(rad2), sin(rad2)) * r;
        d = min(d, line(p, v1, v2));
    }
    return d;
}

5.3. 三角形の内外判定をする

星のdistance functionをつくるには、内側に$p$が来たときに距離が負の値（or $0$）にならないといけません。
星の内外判定をするための準備段階として、三角形の内外判定をする関数をつくります。

さて、三角形の内外判定ですが、衝突判定でよく使われるアルゴリズムとして外積を使ったものがあります。
それを今回は使いましょう。

float crs(vec2 v1, vec2 v2) {
    return v1.x*v2.y - v1.y*v2.x;
}

// 3点v1, v2, v3を頂点とする三角形の内部にpはあるか？
bool innerTriangle(vec2 p, vec2 v1, vec2 v2, vec2 v3) {
    float c1 = crs(v2-v1, p-v1);
    float c2 = crs(v3-v2, p-v2);
    float c3 = crs(v1-v3, p-v3);
    return (c1>0.0&&c2>0.0&&c3>0.0) || (c1<0.0&&c2<0.0&&c3<0.0);
}

アルゴリズムを知りたい人は
点と三角形の当たり判定( 内外判定 )
にわかりやすい説明がありましたのでそちらを見てください。

5.4. 星のdistance functionをつくる

• $\left(0,0\right)$と$\left(r\cos (0\cdot \frac{2\pi }{5}),r\sin (0\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (2\cdot \frac{2\pi }{5}),r\sin (2\cdot \frac{2\pi }{5})\right)$
• $\left(0,0\right)$と$\left(r\cos (1\cdot \frac{2\pi }{5}),r\sin (1\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (3\cdot \frac{2\pi }{5}),r\sin (3\cdot \frac{2\pi }{5})\right)$
• $\left(0,0\right)$と$\left(r\cos (2\cdot \frac{2\pi }{5}),r\sin (2\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (4\cdot \frac{2\pi }{5}),r\sin (4\cdot \frac{2\pi }{5})\right)$
• $\left(0,0\right)$と$\left(r\cos (3\cdot \frac{2\pi }{5}),r\sin (3\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (0\cdot \frac{2\pi }{5}),r\sin (0\cdot \frac{2\pi }{5})\right)$
• $\left(0,0\right)$と$\left(r\cos (4\cdot \frac{2\pi }{5}),r\sin (4\cdot \frac{2\pi }{5})\right)$と$\left(r\cos (1\cdot \frac{2\pi }{5}),r\sin (1\cdot \frac{2\pi }{5})\right)$

これらの3点を頂点に持つ三角形の内側であれば距離を負にすればよいです。（具体的には$-1$をかける。）

つまり、コードにするとこうなります。

float starPolygonFive(vec2 p, float r) {
    float d = INF;
    for (int i=0; i<5; i++) {
        float rad1 = 2.0*PI*float(i)/5.0;
        float rad2 = 2.0*PI*float(i+2)/5.0;
        vec2 v1 = vec2(cos(rad1), sin(rad1)) * r;
        vec2 v2 = vec2(cos(rad2), sin(rad2)) * r;
        bool flg = innerTriangle(p, vec2(0.0), v1, v2);
        d = min(d, line(p, v1, v2) * (flg?-1.0:1.0));
    }
    return d;
}

これで星のdistance functionができました！やったね！

ついでにn芒星まで拡張させちゃいましょう。

#define ITER_MAX 10000

float starPolygon(vec2 p, int n, int m, float r) {
    float d = INF;
    for (int i=0; i<ITER_MAX; i++) {
        if (i >= n) break;
        
        float rad1 = 2.0*PI*float(i)/float(n);
        float rad2 = 2.0*PI*float(i+m)/float(n);
        vec2 v1 = vec2(cos(rad1), sin(rad1)) * r;
        vec2 v2 = vec2(cos(rad2), sin(rad2)) * r;
        bool flg = innerTriangle(p, vec2(0.0), v1, v2);
        d = min(d, line(p, v1, v2) * (flg?-1.0:1.0));
    }
    return d;
}

5行目でfor (int i=0; i<n; i++) {とかきたいところですが、GLSLではループ時のカウンタ変数を変数と比較することができないので、ITER_MAXという適当な定数を使っています。

6. 2つの図形間での変形

やっと円と星の描画が終わりました。次は、これらの間で変形をしていきます。

どうやって変形させるか？
せっかくdistance functionを用いたのでそれを利用するよりほかはないでしょう。

min関数を使います。

例えばこんなのは、

vec3 calc(vec2 p) {
    float d1 = ellipse(p, vec2(0.3));
    float d2 = starPolygon(p, 5, 2, 0.6);
    float d = min(d1, d2);
    return vec3(pow(clamp(1.0-d, 0.0, 1.0), 5.0));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord ) {
	vec2 uv = (2.0*fragCoord.xy-iResolution.xy) / min(iResolution.x, iResolution.y);
	fragColor = vec4(calc(uv),1.0);
}

これでOKです。4行目のfloat d = min(d1, d2);がポイントです。

あとは、星の半径を時間に応じて動かすだけです。

vec3 calc(vec2 p) {
    float r = (sin(iGlobalTime*PI/2.0)+1.0)*0.4;
    float d1 = ellipse(p, vec2(0.3));
    float d2 = starPolygon(p, 5, 2, r);
    float d = min(d1, d2);
    return vec3(pow(clamp(1.0-d, 0.0, 1.0), 5.0));
}

calc関数をこのように書き換えてみましょう。
とりあえず、円と星との間での変形には成功しました！

7. smooth minimum

変形はできたけど、できたモーションはなんか微妙・・・。
不連続感というか刺々しいというか見てて面白みがないです。

原因はmin関数を使ってるから。
この関数のせいで、2つのdistance functionの†接合部分†が滑らかではなくなってるのです。

じゃあどうするか？

smooth minimumを使いましょう。

普通のmin関数の定義は

$\mathrm{m}\mathrm{i}\mathrm{n}(a,b)=\{\begin{array}{ll}a& (\mathrm{i}\mathrm{f}a\le b)\\ b& (\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e})\end{array}$

それに対しsmooth minimumを用いたsmin関数の定義は

$\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}(a,b,k)=-\frac{\log \left({\mathrm{e}}^{-ka}+{\mathrm{e}}^{-kb}\right)}{k}$

となります。（正しくはexponential smooth minimumです。）

• $f(x)=\cos (10x)$
• $g(x)=\log x$
• 赤色 $\cdots y=f(x)$
• 橙色 $\cdots y=g(x)$
• 緑色 $\cdots y=\mathrm{m}\mathrm{i}\mathrm{n}(f(x),g(x))$
• 青色 $\cdots y=\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}(f(x),g(x),2)$
• 紫色 $\cdots y=\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}(f(x),g(x),10)$

グラフで見てみると$\mathrm{m}\mathrm{i}\mathrm{n}$に対し、$\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}$だと2つの関数が滑らかに†接続†されていることがわかります。

それではさっそくsmooth minimumを使ってみましょう。

#define EPS 1e-10

float smin(float a, float b, float k) {
    float res = exp(-k*a) + exp(-k*b);
    return -log(res)/k;
}

vec3 calc(vec2 p) {
    float r = (sin(iGlobalTime*PI/2.0)+1.0)*0.4+EPS;
    float d1 = ellipse(p, vec2(0.2));
    float d2 = starPolygon(p, 5, 2, r);
    float d = smin(d1, d2, 3.5);
    return vec3(pow(clamp(1.0-d, 0.0, 1.0), 5.0));
}

これで滑らかに変形するようになりました。

smooth minimumについてもっと詳しく知りたい方はIQさん（神みたいな存在の人）の記事がありますのでこちらを見てください。

8. 仕上げ

これで所期の目的は達成しました。
実際に動くコードはこちらにあります。

あとは、適当に着色とか変形とかしてポってやるとこんなのもできます。

色々遊んでみましょう！

9. 補足情報

今回用いたdistance functionは、レイマーチングというものによく使われます。
レイマーチングも面白いので、興味ある方は是非やってみてください。
三次元版のdistance functionについてはこちらに大量の例があるので、見てみると参考になるかと思います。

また、今回はアニメーションに三角関数を用いました。「周期性を出したいから三角関数を使う」は間違ってないですが、動き方としては単調で面白くないです。Easing functionという便利なものがあるので使ってみると表現の幅が広がるでしょう。（もちろん三角関数を使ったほうがいい場面もあります）

ところでこの記事に関連するものとして、GPU で暖を取りたい人のための GLSL Advent Calendar 2016というAdCがあります。GLSLやGPUを使った創作に興味がある人にとっては特に面白いんじゃないでしょうか。
（自分もなにか記事を書きたいなあ・・・でも、ここに書けるほどのネタもないしなあ、と悩んでいる間にあっという間に埋まってました＞＜。来年こそは参戦したいし、技量ももっとつけたい）

10. おわりに

今回の内容はここまでです。ここまで読んでくださった方、ありがとうございました。

もっと色々書きたいこともあるので、またいつか続編を書きたいと思います。
そのときも是非読んでいただけるとうれしいです。

以上、「クリエイティブコーディングをしよう！【その１】GLSLで遊ぶ」でした！
明日のtraP Advent Calendar 2016はDavid、gotoh、satorikuの3名がお送りいたします。