点列をCatmull Rom曲線を通してBezier曲線列に変換する

前回，レイトレーシングでBezier曲線を描画する方法を紹介しました．しかし，Bezier曲線そのままでは意図した図形を描画することは大変です．今回は，特定の点列を通るCatmull Rom曲線という曲線をBezier曲線に変換する方法を紹介します．
当初，B-spline曲線からBezier曲線に変換する，Boehm's algorithmを使用しようかと考えたのですが，設定等が複雑なことと，点列を必ず通る精密な図形を描きたかったため，点列を必ず通るという特性を持つ．Catmull Rom曲線を使うことにしました．

scheme(Gauche)による実装はこちら
scheme-raytrace/points.scm at master · soma-arc/scheme-raytrace · GitHub
Processingによる実装も公開している人がいました．
Mike's Processing code snippet repository
この方法と，前回紹介したBezier曲線をレイトレする手法を用いて描いた曲線がこちらです．

各種曲線の概要はこのサイトが詳しいです．
デジタル・フロンティア-Digital Frontier | DF TALK | スプライン曲線の話
また，変換を理解するにあたって，このサイトが非常に勉強になりました．
t-pot 『3次曲線』
これに従い，順に曲線を見ていくことで，自然と変換の方法がわかります．

Ferguson/Coons曲線

基本となる曲線です．三次曲線を表すため，始点と終点，初速度，最終速度の4種のパラメータを用います．これを用いた曲線の計算式はt-potにすべて載っていますが，ここではあまり重要でないので，記述しません．Bezier曲線もCatmull Rom曲線も，この曲線を元に速度の求め方を考えていきます．

Bezier曲線

Ferguson/Coons曲線において，速度を直接決めるのは割と面倒ですので，指定された制御点から，速度を求めてやるようにします．これが，三次のBezier曲線です．を始点，終点としたとき，速度を以下のようにして求めます．


３という係数がBezier曲線の次数ですね．速度が求められれば，Ferguson/Coons曲線の式が使えます．

Catmull Rom曲線

Catmull Rom曲線では，点列すべてを通る曲線を定めることができます．点列のセグメントの前後合わせて4つの点を合わせてパラメータを決めます．ある点における速度を前後の移動量を足し合わせたものにというパラメータをかけたものとします．この値は通常，0.5が使われるため，基本的には，移動量の平均値をとることになります．


ここで注意が必要なのが，点列における始点，終点です．始点，終点では速度を求めるために必要な点がないため，違う速度を使ってあげる必要があります．t-potでは，最初と最後の速度の平均速度を使うとしています．僕は面倒だったので，実装する際には，あらかじめ点を挿入しておいてある，と仮定して無視しました．
パラメータに関しては，このpdfで知りました．
https://www.cs.cmu.edu/~462/projects/assn2/assn2/catmullRom.pdf

点列からBezier曲線の制御点を求める

さて，Catmull Rom曲線と同じ挙動をとるBezier曲線の制御点を求めてみましょう．Bezier曲線とCarmull Rom曲線のは等しくなることを利用して，方程式を解きます．ベジェ曲線の制御点をとすると点列の四点，を用いて以下のように求まります．


先に述べたように，は普通0.5です．これで，点列のあるセグメントに関するBezier曲線の制御点を求めることができました．これを点列全体に関して行えば，Bezier曲線列を得ることができます．

おわり

点列をCatmull Rom曲線を用いてBezier曲線に変換する方法を紹介しました．この方法と前回の記事の内容をあわせることで，点列で表現される図形をレイトレで描画することができます．しかし，点列が多い図形だと曲線の数が増えて，どうしても時間がかかってしまいます．BVHなどを用いて効率よくレイとの交差計算を行わなければ，とても使えたものではないでしょう．
Bezier曲線といえば，内分点をとっていく導出しか知りませんでしたが，今回見たような意味づけは知らなかったので勉強になりました．

レイトレでベジェ曲線を描画する(Ray tracing for curves primitive)

レイトレーシングで曲線を描画する方法をまとめます．とりあえず3つほど見つけたのですが，今回はこちらの方法を実装してみました．

• Koji Nakamaru and Yoshio Ono, RAY TRACING FOR CURVES PRIMITIVE, WSCG 2002.

ベジェ曲線を直接幅を持ったリボンのようなものとして描く方法で,髪の毛などの細いものを描画するのに向いているそうです.
レンダリング結果はこのような感じになりました．


青い球の位置が制御点です．太くしすぎたため，少しがたつきがでています．
アルゴリズム自体は割と単純で，擬似コードも載っているので実装は簡単だと思いますが，その意味を考えると結構つっかかるところが多かったため，メモしておきます．

scheme(Gauche)で書いたソースコードはこちら．このパストレーサでは構造体をベクタで表現しているため，コードはかなり汚いです…
scheme-raytrace/bezier.scm at master · soma-arc/scheme-raytrace · GitHub
ベジェ曲線の接線や分割に関してはこちらのサイトを参考に実装しました．
ベジエ曲線について - s.h’s page

この論文でも紹介されている他の方法はこちら

• レイマーチングを用いて，generalized cylinderとして描画する方法

Sphere Tracing: A Geometric Method for the Antialiased Ray Tracing of Implicit Surfaces | Computer Graphics Illinois

• 連続した球の移動として描く方法

Ray tracing objects defined by sweeping a sphere

アルゴリズムの流れ

• 論文のアルゴリズムはz軸が高さ(Z-UP)になっていますが，僕が今回使ったパストレーサーはy軸が高さ(Y-UP)となっているので，処理を行う前に簡単な変換を行っておきます.
• あらかじめ指定されたベジェ曲線の幅からベジェ曲線の分割の深さを推定
• 全体に座標変換をかけ，原点がレイの原点，レイの方向がz軸正の方向になるように変換します．
• ベジェ曲線を分割
• 一定の領域内に入っている分割されたベジェ曲線のセグメントの中で，レイがベジェ曲線の内側に入っているものを見つける．
• その中からレイに最も近いものを見つける

最大分割数の推定

ベジェ曲線を分割していくことで，接触を判定していきますが，分割する前に，その分割数を描画する曲線の幅から推定します．これにはWang’s methodと呼ばれるものを使います．この式にはというパラメータがありますが，論文ではという値を提案しています．これは髪などを描画する際に有効な値のようで，太い曲線を描画するためには分割数が小さくて，先の画像のように正しい帯にならないようです．太い曲線を描く際には注意が必要です．

ベジェ曲線を変換する

次に，ベジェ曲線の制御点に変換行列をかけて，レイが原点中心，z軸方向となるように座標系を変換します．以下のような変換行列が提示されています．

(Oはレイの原点，lはレイの方向，)
d=0の時はの符号に応じたx軸に関する90°回転とするようです．

まず，レイの原点まで平行移動をかけ，レイの向きがｚ軸の正の方向を向くように変換してやる必要があります．向きの変換には，z方向がレイの向きである，以下のような正規直交基底を用いることで，そのような変換をかけることができます．

正規直交基底の求め方は以下のようになります．
まず，レイの方向をz方向とすることから，z軸のベクトルが決まります．

次に，x方向を考えますが，この時に，y軸方向は０とすることで，z軸方向のベクトルをXZ平面において90°回転して正規化することで作る事ができます．

最後のy方向のベクトルは，二つのベクトルに直交するという条件

から，

が求められます．

分割

曲線の幅を持った矩形を考え，曲線の矩形領域がその矩形にかぶっていれば，ベジェ曲線を分割するという操作を最大分割数に到達するまで繰り返します．ベジェ曲線は凸包性があり，全体を矩形で囲えるため，それを利用することで簡単に領域を見ることができます．

ベジェ曲線の分割は，パラメータが0.5の部分で分割していきます．分割が終わると，そこにはベジェ曲線の列ができることになります．
実装としては再帰関数になります．scheme実装は多値を返すようにしました．
関数の引数はcがベジェ曲線のセグメント，v0, vnがベジェ曲線全体を見たときのパラメータの値，tがレイとの距離です．

(converge
 (lambda (depth c v0 vn t)
   (let* ((b (g:bounding-box c 0 0))
          (b-min (g:aabb-min b))
          (b-max (g:aabb-max b)))
     (cond ((or (>= (v:z b-min) t) (<= (v:z b-max) 0.000001)
                (>= (v:x b-min) width1) (<= (v:x b-max) (- width1))
                (>= (v:y b-min) width1) (<= (v:y b-max) (- width1)))
            (values #f #f))
           ((< depth 0)
            ;; ...
            )
           (else (let*-values
                  (((vm) (/ (+ v0 vn) 2))
                   ((cl cr) (bezier-split c 0.5))
                   ((hit-left? nt-left) (converge (- depth 1) cl v0 vm t))
                   ((hit-right? nt-right) (converge (- depth 1) cr vm vn t)))
                  (if (and hit-left? (< nt-left t))
                      (set! t nt-left))
                  (if (and hit-right? (< nt-right t))
                      (set! t nt-right))
                  (values (or hit-left? hit-right?) t)))))))

ここで気になったのが，擬似コード上では，

return
 converge(depth,   cl,   v0,   vm,   t)
 ||   converge(depth,   cr,   vm,   vn,   t);

と表現されている部分です．多くのプログラムは，||演算子で先に評価された値がtrueとなった場合に，後の演算は評価されないように最適化されるはずです．しかし，そうなるとレイと最も近い場所を得ることができません．擬似コードの場合，そのような最適化を考えないのが前提となるのでしょうか…そのため，今回の実装ではorを使いませんでした．

位置を見る

分割後は，各ベジェ曲線に対して，レイ（原点）が現在見ている曲線の内側に入っているかを考えます．ある点Qが，正しい領域にいるときに，以下の式を満たすそうです．

(dPは接線，Pは制御点)
今回考えるべきなのは，レイの位置なので，Qは原点の０となります．
また，カーブが急な場合はこのままではうまく判定できないため，となる場合に接線の向きを反転してやります.
このプロセスのコードを抜き出すを以下のようになります．一つ注意が必要なのは，接線との内積を取る時は，XY平面上でとり，Zの値は無視する必要があるということです．ベジェ曲線の制御点は三次元上の点として扱ってきたので，そのままの内積をとってしまいがちです．

(let* ((dir (v:diff (bezier-cp c 3) (bezier-cp c 0)))
       (dp0 (bezier-tan-vec c 0))
       (dp0 (if (< (dot2d dir dp0) 0)
                (v:scale dp0 -1)
                dp0))
       (dpn #f)
       (w #f)
       (v #f)
       (p #f))
  (if (< (dot2d dp0 (v:scale (bezier-cp c 0) -1)) 0)
      (values #f #f) ;; return
      (begin
       (set! dpn (bezier-tan-vec c 1))
       (set! dpn (if (< (dot2d dir dpn) 0)
                     (v:scale dpn -1)
                     dpn))
       (if (< (dot2d dpn (bezier-cp c (- n 1))) 0)
           (values #f #f) ;; return
           ;; ...
           ))))

レイの位置のパラメータを求める

次に，レイの位置のベジェ曲線のパラメータを求めます．この時，考えるベジェ曲線は直線として考えます．この直線に，原点から垂線を引き，その交点の位置からベジェ曲線のパラメータを求めます．

論文で出されているパラメータ，

をというベクトルで書きかえると

となり，wが，正射影ベクトルの長さで，をw:(1-w)で分割するときのパラメータであることがわかります．
あとは，分割前のベジェ曲線全体のパラメータである，から，wで分割する値を求め，ベジェ曲線上の点を計算します．

ここでもう一つ気になる点があり，直観的にvはベジェ曲線全体で考えたときのパラメータとなるはずなのですが，論文の擬似コードでは，現在見ているベジェ曲線のセグメントにパラメータを与えて評価しています．逆に，全体のベジェ曲線でパラメータを評価すると正しく動きませんでした．また，十分短い線分ならば，わざわざこのようなパラメータを求める必要がない気もします．ここだけは意図が理解できませんでした．

レイとの近さを比較する

最後に，計算された点が，幅をはみ出ていないか等を確かめ，レイとの近さ，すなわちzの値を比較して返します．

(begin
 (set! w (+ (* (v:x dir) (v:x dir))
            (* (v:y dir) (v:y dir))))
 (if (= w 0)
     (values #f #f)
     (begin
      (set! w (/ (+ (* (v:x (bezier-cp c 0)) (v:x dir))
                    (* (v:y (bezier-cp c 0)) (v:y dir)))
                 (- w)))
      (set! w (clamp w 0 1))
      (set! v (+ (* v0 (- 1 w)) (* vn w)))
      (set! p (bezier-point c v))
      (if (or (>= (+ (* (v:x p) (v:x p))
                     (* (v:y p) (v:y p))) width2)
              (<= (v:z p) 0.0001)
              (< t (v:z p)))
          (values #f #f)
          (values #t (v:z p))))))

法線について

法線ベクトルに関しては，特に論文中で言及がなかったので，レイの向かってきた方向の逆をそのまま返してやりました．円柱だと考えて返すこともできるでしょうけれど，細いものだとそんなに変わりはないかもしれません.

おわり

少々つっかかるところもありましたが，レイトレ―シングを用いてベジェ曲線を描画する方法を紹介しました．ベジェ曲線を直接扱うのは少し面倒ですが，曲線が描けるようになると表現の幅も広がると思います．次は点列からベジェ曲線列を求める方法をまとめます．

Djangoのtemplate engineをstandaloneで使う

DjangoはPythonのWebアプリケーションフレームワークです．そのテンプレートエンジンのみを使いたいときの設定方法を調べてみたところ新旧様々な情報があり，苦労しましたのでまとめておきます．しかし，後から気づいたのですが，jinja2というテンプレートエンジンがDjangoライクな文法を提供しているようで，テンプレートエンジンのみを使う場合はこちらの方が良いかもしれません…
ちなみに，DjangoのテンプレートエンジンはdjulaというCommon Lisp実装も存在しています．

環境はpython 2.7.10, django1.9.6です．Django自体はpipを用いて，

pip install django==1.9.6

でインストール．もしくは，現時点ではGitHubからリポジトリをクローンしてきても動きました．

import django
from django.conf import settings
from django.template import Context, Template, loader

settings.configure(TEMPLATE_DEBUG=True,
                   TEMPLATES=[{"BACKEND":'django.template.backends.django.DjangoTemplates',
                               "DIRS":"./"}])
django.setup()

t = Template("My name is {{ name }}")
c = Context({"name":"foo"})

print(t.render(c))

t = loader.get_template("index.tmpl")
c = Context({"users":
                 [{"id":0, "name": "hoge"},
                  {"id":1, "name": "foo"}]
             })

print(t.render(c))

index.tmpl
{% for user in users %}
   {{ user.id }} ... {{ user.name }}
{% endfor %}

テンプレートエンジンのみを使うためにはsettings.configureで設定を記述した後に，django.setup()を呼び出す必要があります．BACKENDが使用するテンプレートエンジンの設定で，DIRSが，get_templateで読み込まれるテンプレートを置くディレクトリになるはずです．
ContextにはJSONを読み込んでそのまま流し込むこともできます．ただし，真偽値(True, False)はJSONの値としてそのまま書けないので注意が必要です．

Orbit trapを考える

マンデルブロ集合などでよく使われるカラーリングのアルゴリズムにOrbit trapと呼ばれる手法があります．以前から試してみたかったのですが，ようやくある程度理解でき，このような作品をつくることができました．

Kleinyan Cat at Shadertoy
今考えるとたいした事はないアルゴリズムだったのですが，割と理解するのに時間がかかったのでまとめておきます．やはり，実際に手で書いて試してみるのが大事ですね…

Orbit trapとマンデルブロ集合

マンデルブロ集合のようなフラクタルは，Escape-time Fractalとも呼ばれ，ある式を繰り返し計算していき，その計算結果の収束，発散をみることで描かれるフラクタルです．
マンデルブロ集合の計算の仕方を簡単にまとめると，複素平面上に格子点を取ってその座標を複素数Cとし，以下の漸化式を繰り返し計算します．

が無限大に発散しないものがマンデルブロ集合となります．また，となると無限大に発散することがわかっています．
複素平面上の複素数の積は点の回転，拡大，縮小，複素数の和は点の平行移動になりますので，複素平面上に点の軌道が描かれることがわかります．
この軌道上にtrapとなる図形などを設定しておき，そこからの距離や角度を保持して，これを用いてカラーリングを行います．

点を用いたOrbit trap

まずは最も簡単な点を用いたOrbit trapを見てみます．

private double getDistance(Complex c,
                           Complex point,
                           int maxIteration)
{        
    double distance = 1e20;
    Complex z = new Complex(0, 0);
        
    for(int i=0; i<maxIteration; i++)
    {
        //Perform Mandelbrot iteration
        z = z.multiply(z);
        z = z.add(c);
               
        //Set new distance dist = min( dist, |z-point| )
        Complex zMinusPoint = new Complex(z);
        zMinusPoint = zMinusPoint.subtract(point);
            
        double zMinusPointModulus = zMinusPoint.magnitude();
        if(zMinusPointModulus < distance)
            distance = zMinusPointModulus;
    }
        
    return distance;
}

Orbit trap - Wikipedia, the free encyclopedia

レンダリング結果は載っていませんが，アルゴリズムの概要をWikipediaから引用しました．こうしてみると単純ですね．点を用いたOrbit trapで，あらかじめ設定した点と最も近い時の距離を用いてレンダリングを色を決定します．これとほぼ同じアルゴリズムで(-0.5, 2.0)の点でトラップし，distanceをそのままプロットしたのが以下の画像です．この画像はShadertoyのMandelbrot - orbit trapsを少し変更して生成したものです．

黒い点が連なっています．マンデルブロ集合では点の軌道がわかりにくく，点が連なる理由が見えてきません．個人的になれている接触ショットキー群でも考えてみたいと思います．

Orbit trapと接触ショットキー群

４つの円で構成される接触ショットキー群を考えます．この画像は，マンデルブロ集合と同じように，複素平面上に格子点を取り，その座標が４つの基本となる円盤に載っていた場合にその円に関する鏡映反転を行う，という操作を点がすべての円の外に出るまで繰り返し，その繰り返した回数に応じて色を塗った画像です．この点にも軌道ができますので，Orbit trapを行うことができます．
Kleinyan Catの191行目あたりのコメントアウトを外してCompileし直すと，次のような画像が得られます．

    vec3 hsv = hsv2rgb(vec3(0.04 * loopNum,1.0,1.0));
   
    col.xyz = mix( hsv, col.xyz, col.w );
    fragColor = col;
 
    if(trap1/1000. < 0.8){
    	fragColor = vec4(vec3(trap1/1000.), 1.0);
    }else{
    	fragColor = vec4( hsv, 1.0);
    }

これはt1c = vec2(-100, -100)と軌道の点が最も近い時の距離を保持し，スケーリングしてプロットしたものです．パラメータが動いても動かない点がありますが，これが最初に設定した点(-100, -100)です．
二つの円盤に点がかぶってしまっているのであまりよくない画像ですが，この図形の点の大体の軌道を考えてみるとこんな感じになると思います．

それぞれの点は反転を行うと，円の中心から点の方向に，円盤ひとつずつ降りていくような挙動を取ります．トラップする点，t1cを弄ってみるとより点列ができる理由がわかりやすいかもしれません．
ちなみにこの図形は二つの領域に円の極限（極限集合）で分けられており，この内側の点は内側に，外側の点は必ず外側に移るようになっています．つまり，内側にトラップする点を置いておくと，点列は内側にしか現れません．t1cを(0, 0)に置いてみるとよくわかります．

Orbit trapを用いたカラーリング

アルゴリズムは案外単純であることがわかりましたが，果たしてこれで面白いカラーリングができるのでしょうか…iq先生のMandelbrot - orbit trapsでは点との距離を見るトラップとaverage trap distance to lineというトラップを行い，それらをうまく変換しているようです．ズームや，経過時間に応じてトラップする点をずらすなど，色々とうまくやっているようで，使いこなすには経験とセンスが必要な気がします…定石などがあるのでしょうか．
今回はlineに関するtrapは試していませんが，64行目あたりの

col += 2.0*sqrt(c1*col1*col2);

を

col += c2;

にすると大体の挙動がわかります．

Bitmap Orbit traps

さて，次にtrapに画像を用いるBitmap Orbit trapを行ってみたいと思います．こちらは絵を移していくものなので，普通に色付けするよりは簡単に面白いものができるのではないかと思います．Bitmap Orbit trapの実装はFractal Nyancatが参考になります．Kleinyan Catのトラップの処理は以下のようになります．

if (col.w < 0.1)
	col = getNyanCatColor( pos - cPos[0] * scale);
if (col.w < 0.1)
	col = getNyanCatColor( pos - cPos[1] * scale);
if (col.w < 0.1)
	col = getNyanCatColor( pos - cPos[2] * scale);
if (col.w < 0.1)
	col = getNyanCatColor( pos - cPos[3] * scale);

cPos * scaleがトラップする点です．点の位置から引くことで画像の始点からの位置を求めています．画像は原点から四円の中心への直線の延長直線状に置いています．絵がない部分には，透明度が設定されているので，絵が取得されていれば画素値は取得されないようになっています．
getNyanCatColorという関数は，Fractal Nyancatで使われているものを整えました．

vec4 getNyanCatColor( vec2 p ){
	p += vec2(100, 100);//offset
	p = p / 200.;
	p = clamp(p,0.0,1.0);
	p.x = p.x*40.0/256.0;
	float fr = floor( mod( 20.0*iGlobalTime, 6.0 ) );
	p.x += fr*40.0/256.0;
	return texture2D( iChannel0, p);
}

トラップに使用する点の中心に画像を置きたいので，適当に移動した後にスケーリングします．ここら辺は実験的に決めましたが，ちゃんと考えれば一般的な解はあるはずです．
clampすることで，0以下は0に，1以上は1に丸まられ，画像から離れすぎた場所は角の画素が参照されることになります．こうすることで，パフォーマンスに影響がでる条件分岐を削る事もできます．参照しているテクスチャは6フレームの猫が収められているので，うまく参照してやります．トラップに用いる最初の猫の位置はこのような感じです．(フレームの外側に2体隠れていて，合計4体描かれているはずです．)

点の軌道上に，最初に置いた猫の画像があれば，そこから画素値を取得し，その色を取得する感じです．結果としてこのような画像が得られました．

移された猫は隣り合う円盤同士で鏡写しになり，歪みが入っているのがわかるでしょうか．これが円に関する鏡映変換の作用です．また，この結果はクライン群を生成する４つの生成元で一つの画像を写していくことと同じようです．以前描いた以下のようなものが大元となる一体の蝶を4つの生成元で移したものです．
GitHub - soma-arc/KleinianWalker: Draw limit set of Kleinian group and orbit of figure

生成元で直接移していくためには，生成元の組み合わせの木構造を幅優先で探索していかなければならないのですが，探索が深くなるにつれて，指数オーダーで計算量が増えていくので，細かく描こうとするほど時間がかかります．Bitmap Orbit trapが使えれば，そのような探索が必要でなくなり，高速な描画ができるということになります．

おわりに

今回は点と画像を用いたOrbit trapのみを試しましたが，他の種類のものも試す必要がありますね．また，Bitmap Orbit trapを用いることである画像を生成元で移していくような作品も高速に描画できることがわかりました．このようなテクニックはまだまだあるようで，一つ一つ実際に実装して理解していく必要がありそうです．日々精進あるのみですね．

Reference

1. Orbit trap - Wikipedia, the free encyclopedia
2. geometric orbit traps by Inigo Quilez
3. bitmap orbit traps by Inigo Quilez
4. Mandelbrot - orbit traps by iq
5. Fractal Nyancat by iq
6. Make Fractals of Images

WindowsでOCamlを学ぶ環境を作るメモ

前々からやろうと思っていたOCaml学ぶ環境を整えました．WindowsでOCamlはつらいという話をよく聞きますが，入門する程度には問題ないでしょう．
ちなみに読もうとしている本はプログラミング in OCamlです．
プログラミング in OCaml 〜関数型プログラミングの基礎からGUI構築まで〜 | Gihyo Digital Publishing … 技術評論社の電子書籍

僕の使っているEmacsはGnuPack11.00のものなのですが，windows側で設定されたパスをうまく引き継いでくれません．いくつかのパスを手で設定してやる必要がありましたので，メモしておきます．

処理系のインストール

ここでインストーラをダウンロードできます．インストール時にはcygwinをインストールするか聞かれますが，僕はインストールしませんでした．
The OCaml installer for Windows

Emacsの設定

tuareg-modeというのが良く使われているそうです．
GitHub - ocaml/tuareg: Emacs OCaml mode

インストールした後は以下の記事を参考に設定しました．
Emacsでocamlを書く設定 - 雑記

(add-to-list 'auto-mode-alist '("\\.ml[iylp]?" . tuareg-mode))
(autoload 'tuareg-mode "tuareg" "Major mode for editing OCaml code" t)
(autoload 'tuareg-run-ocaml "tuareg" "Run an inferior OCaml process." t)
(autoload 'ocamldebug "ocamldebug" "Run the OCaml debugger" t)

基本的な設定が済んだら.mlファイルを開いた状態でC-c C-s (tuareg-run-ocaml)でOCaml toplevel to runというプロンプトでコマンドを聞かれた後にREPLが起動します．
パスが通っていない場合はそのまま動かないので以下の設定のうちどちらかを加えます．

;;パスに加える．
(add-to-list 'exec-path "path/to/OCaml/bin/")
;;プロンプトの初期値を設定する．
(setq tuareg-interactive-program "/path/to/OCaml/bin/ocaml")

起動後にexception typetexp.errorといったエラーが出る場合は恐らくライブラリへのパスが通っていませんので環境変数を設定してやります．

(setenv "OCAMLLIB" "C:/tools/OCaml/lib")

一応環境が整いました．

TokyoDemoFest2016に参加しました

2月20-21日に行われた日本で唯一のデモパーティ，Tokyo Demo Festに参加してきました．GLSLのFragment Shaderでグラフィックスを表示するGLSL Graphics Compoに"Indra's Bubbles"という作品を出して二位をいただきました．

GLSL Graphics Compoで二位をいただきました.(左側はFMS_Catさん)日々精進あるのみですね.今日も一日ありがとうございます. #TokyoDemoFest pic.twitter.com/NF8MUQYAqk

— 蘇摩清良 (@soma_arc) 2016, 2月 21

作品の様子

#TokyoDemoFest pic.twitter.com/0dFBbDLVtA

— tokkyo (@tokkyo) 2016, 2月 21

Shadertoyに投稿したコード，動くものはこちら
Indra's Bubbles

デモとの出会い

元々，デモシーンにはレイマーチングを通じて出会いました．sphairahedronの高速な描画にレイマーチングを用いて物体同士のブール演算を用いて描画する方法が適していることがわかったため，これに用いるためにレイマーチングを学びました．次にMandelboxやPseudo-Kleinianの存在です．フラクタルアーティスト界隈でよく使われるそれらのフラクタル図形にはクライン群に見られるような特徴をみることができることに気付きました．
www.youtube.com
これらはクライン群関連の論文でも見たことがありませんでしたので調査を始め，結果としてデモで使われるようなテクニックを用いて本物のクライン群の高速描画アルゴリズムを指導教員と共に開発することができました．

作品について

中央に正八面体の頂点に配置した球で構成される群によるクライン群とその周囲にPseudo-Kleinianの式で構成されるフラクタル図形を配置したものです，クライン群とPseudo-Kleinianの一部は鏡面反射されるようにいじってあり，特に立方体の隅や特定の一面は全てが反射されています．鏡面反射する部位については，Orbit trapを行って計算していますが，Orbit trapそのものの理解が不十分であり，他の人の作例を見て，実験的につけたものです．
個人的には完全な球形とその次の形が良い感じに反射してくれてお気に入りです．
しかし，中央のクライン群の描画アルゴリズムありきの作品であり，レンダリング技術やエフェクト等がなっていません．特にフォグの実装が駄目なのか，低い解像度だとノイズが目立ちます．レンダリング技術を極め，より綺麗な映像にしたいところです．
中央のクライン群の可視化アルゴリズムについては，近いうちにまとめた文章を出す予定です．

イベントについて

初めてのデモパーティーでしたが，会場の雰囲気としては，ぎっちりとイベントが詰まっているわけではなく，適度に空白の時間があり，その間はパソコンに向かって作業する，もしくは他の人と交流するような感じになっていました．今回はセミナーをおこなう場所が分かれており，聞きたい人が聞きにいく感じになっていました．チョコレートや寿司，ピザ，ビールを振舞っている人もいて，パーティーという感じでしたね．ただ，イベント会場のArts Chiyodaの体育館は寒かったです．

おわり

どのデモも非常に面白く，素晴らしいリアルタイム生成のグラフィクス作品を見ることができました．このような作品があるのならば，自分の頭の中にある理想もどうにか作品にすることができるのではないかと，少しだけ希望を持つことができました．
来年もどうにか，クライン群，フラクタルのリアルタイムレンダリングを極めて作品を出したいと思います．また，Combined demo compoは音楽が必要なのでどうにかゴアトランスを作れるようになりたいところ…
それでは，Demoscenerの皆様の一年がより良きものでありますように．

Clojure (overtone + quil)によるコーディングとコードを音楽にするライブコーディング環境の試作

www.youtube.com
github.com

Parensymphonyはコーディングとコードが音楽になるライブコーディング環境の試作品です．ある授業の課題として作り，パフォーマンスをしました．動画を見ていただけるとわかるように，コーディングの操作に音がつきます．また，評価されたコードの構造からフレーズを生成していき，最終的に重なった音楽になります．動画中では，書いたコードで太鼓のパターン生成も行っています．つまり，コードが楽譜とプログラム二つの意味を持っていることになります．
ちなみに名称の由来はCommon LispのJavaScriptトランスレータライブラリであるParenscriptからいただきました．恐らくParenthesisとscriptからつけられた名称だと思います．

なぜつくったのか

ライブコーディングで音楽をする作品には色々な環境，作品があります．Clojure(overtone)を用いたものには以下のような作品，ユニットが存在しています．
Sam Aaron - Hacking Overtone - Live @ Arnolfini on Vimeo
Repl Electric - Live Coding
しかし，これらはLispコードを評価し，その結果として音楽や映像が出力されるものです．コーディング過程やコードそのものが映像，音楽表現になるものを見たことがありませんでした．僕はLispのコードを初めて来たときにこのプログラミング言語はなんと格好良いのだろうと思っておりましたので，Lispコードが主役のものを見てみたかったのです．丸括弧を入力するとジャーンと和音がなって格好良く括弧が現れるというものが僕が最初に夢想したものです．
また，Ｓ式は構文木を表しているわけですけれども，そこから音楽的な構造にならないかとも考えました．詳しくは分かりませんが，実際に楽譜を木構造としてみる音楽理論が存在するようです.（CiNii 論文 -  音楽理論GTTMに基づく木構造を用いたメロディ生成手法）
Lispのコードはprogramとしてもdataとしても解釈できます．そこでprogram = data = scoreを目指しました．

環境

今回はClojureで作りました．SuperCollider，Processingのラッパーであるovertoneとquilを使用しています．
Overtone - Collaborative Programmable Music
quil/quil · GitHub

ちなみにCommon Lispにもcl-colliderというSuperColliderのクライアントが存在するようですね．ただ，グラフィクス等を扱う良いライブラリがないのでClojureを使うことになりました．
byulparan/cl-collider · GitHub

実装

エディタ

エディタ部分は実装する知識，時間がなかったので，一つのフレームにつき，一つのＳ式を編集するという簡単なものです．ここに関しては，Emacsを用いてグローバルキーフックをかけるという方法も考えられましたが，ちょっと面倒そうなのと，やはり映像表現を自由に行いたかったので，quilの上に作りました．
ただ，結局コードを用いた映像表現に関しては文字のアニメーション表現等は経験がなく，今回は実装できませんでした．無理に実装してもゴタゴタしてつまらないだけでしょうし，各フレームが音に合わせてフラッシュするという単純なものにしました．

音

overtoneのチュートリアルで紹介されているplucked-stringのみを使用しています．フレーズはペンタトニックスケールを用いて構成されています．
括弧やスペース，タブを入力すると和音，その他のキーは単音が鳴ります．単音はフレーズが設定されており，入力するごとに順番に鳴り，スペースが入力されると別のフレーズに切り替わります．和音は完全にランダムです．本来はコード進行をつけようとしていたのですが，ペンタトニックスケールに合うような進行がよくわかりませんでした．
音が重なって聞こえにくくならないようにエディタの各フレームでは音の高さが変えられています．例えば，左上のフレームの単音はC2がルートのペンタトニックスケールで構成され，和音はその一つ上の高さになっています．

フレーズの生成

弦楽器でペンタトニックスケールは何でも合って凄いですね．
まず，スケールの山を作ります．

 (let [penta (scale :c3 :pentatonic)]
    (concat penta (reverse penta)))
-> (48 50 53 55 57 60 62 65 65 62 60 57 55 53 50 48)

そこからランダムに歯抜けにした無限シーケンスを作ります．

;;たとえば2つごとに抜きます．
(cycle (take-nth 2 notes))
-> (48 53 57 62 65 60 55 50 48 53 57 62 65 60 55 50  ...)

最初のいくつかを切り捨てて完成です．

;;例えば三番目以降を取得．
(nthrest phrase 3)
-> (62 65 60 55 50 48 53 57 62 65 60 55 50 48 53 57  ...)

途中にreverseを仕込んで反転させても面白いのですが，統一感がなくなったりしたのでやめました．このようにして生成されたフレーズがキー入力で鳴ったり，コードからの自動生成時に使用されます．

コードから音楽の生成

式の要素を順番に見ていき，文字列の長さ分だけフレーズを割り当てます．フレーズは各要素ごとにランダムに選びなおされます．要素の間には休符が入れられます．コードは以下のような感じ．

(defn gen-pattern [key-index n]
  (cond
    (or (number? n) (symbol? n)
        (keyword? n) (char? n)) (take (count (str n)) (get-phrase key-index))
    (empty? n) nil
    :else (concat (gen-pattern key-index (first n))
                  '(rest)
                  (gen-pattern key-index (rest n)))))

現状では要素の文字の長さとＳ式の構造によって休符をつけることのみしかしていません．もっとコードの情報をうまく使えると良いですね．

パフォーマンス

どのようなプログラムを書くか

環境自体が完成しても，パフォーマンスができなければ意味がありません．当初は競技プログラミングや4Clojureの簡単な課題を解くことも考えたのですが，音楽に全く関係ありませんし，4つコードでできることも限られています．また，Lisp(Clojure)についての事前知識がなければ何をやっているかわかりません．
今回はFizzBuzz(FibBuzz FibBuzz in Clojure · GitHub)から太鼓のフレーズを生成してみたところ，非常にいい感じになったのでこれを使いました．FizzBuzzならば，プログラミングを行っている人ならばわかる人も多いでしょう．
今回は非常にうまくはまってくれましたが，これ以外のものを考えるとなると難しいです．
overtoneを用いた作曲アプローチとしてleipzig(ctford/leipzig · GitHub)といった作曲ライブラリのようにノートのシーケンスを用いる方法があります．音の長さ，や高さをリストで与えることで演奏します．[1/3 1/3 1/3]のような感じで音符を与えるなどするのですが，このようにリストを書くようなプログラムは見ていてあまり面白くありません．
今回のように，なんらかのアルゴリズムを記述し，そのアルゴリズムを用いて生成された音楽に新たな要素を加えるというようなものが一番面白いと思います．書いたアルゴリズムがどんなに美しくてもそれが音楽として面白いかはまた別の話ですけれども…
竹内関数による音楽生成は非常に面白い例ですね．
竹内関数で音楽生成 - aikeの日記

タイピング

タイピングすることによって音楽ができるようにもなっていますが，自動演奏されるフレーズとの兼ね合いを考えてタイピングの速度を変える必要があります．あんまり早いと焦ってしまうので遅い方が良いのではないかと思います．

今後の課題

きちんとしたエディタを実装する

言わずもがな．Emacsキーバインド，Pareditあたりは必須です．

一般化

現在のParensymphonyは一つのパフォーマンスのために作られたものになっています．より一般化し，様々なパフォーマンスに対応できるようなフレームワークにすることが最終的な目標になるのだろうと思います．

コードを用いた映像表現

Lispコードをより格好良く見せる．

音楽知識の習得

言わずもがな．より良いフレーズの生成のために．

ゴアトランスを作りたい．

DTMで作れるようになってから．

シンセサイザーに関する知識の習得

自分の欲しい音が作れるようにならないと話になりません．サンプルを探してみてもなかなか欲しいものにはたどり着けませんでした．

非同期制御

お気づきだと思いますが，コードを評価するたびに．音楽が一度止まっています．ここら辺を改良するには，Clojureの非同期的な制御を行わなければならないのですが，僕が学習できていなかったのと，安定性を考えて一度止めて再度再生するようにしています．
かなり気になる部分ですので，改良は必須です．

感想

Clojure自体は初めてでしたが，楽しい言語でした．無限シーケンスが音楽的なことをするのに非常に便利でした．overtoneができるのも分かりますね．しかし，Clojureは関数型的な特性が強くて，そこら辺に手間取る事が多かったです．Common Lispで実装することも考えています．
音楽の知識，経験が足らず，大変でしたが，最終的にLispならではのものができたと思います．どのくらい時間がかけられるかわかりませんがより進化させたいとは思っています．なんにせよ，まずは知識不足を埋めることが一番の課題でしょうか．日々精進あるのみですね．