Apollonian Gasketを描画してみる
皆様、ここの所かなり暑くなってきていますが、いかがお過ごしでしょうか。このところ、お庭のマルベリーの収穫が最盛期を迎えていまして、早朝の5時頃から収穫作業です。夜、帰ってきたら下処理をしてジャムを作ります。今年はこちらもかなり豊作です。しかし、この収穫作業はなかなか大変だったりします。よく見ないと毛虫を掴んでたりしますからねw
さて、本題です。こういった好事家好みのネタはまっとうな方々は手を出さないと言う事がよく分かりました。
誰もマトモにやってないんです。どこ検索してもAIに本格的なフラクタル描かせようなんてヤツはいません。こんな事考えるのは、世間的にはただのモノ好きと言う事かも知れません。
まずは出力結果をご覧いただきましょう。
細部を拡大してみます。
み〜んなくっついているのが良くわかります。
しかし、フラクタルって、魅力的です。エレガントです。このガスケットなんかの整然と隙間に収まって行く円にハアハアしてしまいそうです。まあ、こんな世迷い言はスルーして頂いてかまいません。
今回のポイントは複素平面での処理を行っている点だったりします。
該当部分を抜き出してみます。
// complex numbers
function C(r,i){
this.r = r; //real
this.i = i; //imaginally
}
これが複素数オブジェクト用。複素平面と言っても実際は直交座標系ですので、x,yを実数部と虚数部にあてはめて処理する事になります。しかしながら虚数部にはi(虚数)がくっつきますので単純に計算して良いものではありません。で、専用の演算ファンクションが必要となります。
虚数と言うのは2乗すると-1になる変なヤツですから…
//compute complex value
function Cconj(c){ //conjunction
return new C(c.r, -c.i);
}
function Cadd(c,d){ //add
return new C(c.r+d.r,c.i+d.i);
}
function Csub(c,d){ //sub
return new C(c.r-d.r,c.i-d.i);
}
function Cmul(c,d){ //multi
return new C(c.r*d.r-c.i*d.i,c.r*d.i+c.i*d.r);
}
function Csqrt(c){ //square root
return new C(
Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i))*Math.cos(Math.atan2(c.i,c.r)/2),
Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i))*Math.sin(Math.atan2(c.i,c.r)/2)
);
}
以上が複素数用四則演算+αファンクション群です。引数は全て複素数となります。
続いてcircleオブジェクトも見ておきましょう。
function circle(x, y, rad) {
this.x = x;
this.y = y;
this.r = rad;
this.draw = function() {
var cir = app.activeDocument.pathItems.ellipse (
(this.x+this.r)*scl+ofst, (this.y-this.r)*scl+ofst,
this.r*2*scl, this.r*2*scl);
cir.filled = false;
cir.stroked = true;
cir.strokeColor = clr;
cir.strokeWidth = 0.25;
return cir;
}
}
座標と半径をプロパティに持つオブジェクトです。実際に円を描画するためにはメソッドを呼び出す必要があります。このファンクションはグローバルで定義される、オフセット値と倍率を参照しますので注意してください。
もうひとつ、重要なのがこのファンクションです。
function tng(a, b, c, flg) { //a,b,c are circle objects
var k1 = 1 / a.r;
var z1 = new C(a.x, a.y);
var kz1 = Cmul(new C(k1, 0), z1);
var k2 = 1 / b.r;
var z2 = new C(b.x, b.y);
var kz2 = Cmul(new C(k2, 0), z2);
var k3 = 1 / c.r;
var z3 = new C(c.x, c.y);
var kz3 = Cmul(new C(k3, 0), z3);
var k4p = k1 + k2 + k3 + 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
var k4m = k1 + k2 + k3 – 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
var kz4p = Cadd(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
var kz4m = Csub(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
var k4,kz4,k4b,kz4b;
var cs = new Array();
if (k4p>k4m){
k4 = k4p;
kz4 = kz4p;
k4b = k4m;
kz4b = kz4m;
} else {
k4 = k4m;
kz4 = kz4m;
k4b = k4p;
kz4b = kz4p;
}
var cc = new circle(kz4.r/k4,kz4.i/k4,Math.abs(1/k4));
var dx = a.x – cc.x
var dy = a.y – cc.y
var dr = a.r + cc.r
if (Math.abs(dx*dx+dy*dy-dr*dr)>0.0001) {
cc = new circle(kz4b.r/k4,kz4b.i/k4,Math.abs(1/k4));
}
cs.push(cc);
if (flg) {
cc = new circle(kz4b.r/k4b,kz4b.i/k4b,Math.abs(1/k4b));
cs.push(cc);
}
return cs;
}
3つの接する円を投入すると、それらにくっつく円を返却します。デカルトの定理やら三角関数を利用したジオメトリカルな演算が中心です。
メインルーチンでは初期単位円を生成した後、ランダムに円を生成し、続いてループ内で再帰的に前述のファンクションを呼び出しながら円を描画しています。
以下が全体です。
//Draw Apolonian Gasket…//more information in… http://en.wikipedia.org/wiki/Descartes%27_theorem
var clr = new CMYKColor;
clr.cyan = 0;
clr.magenta = 0;
clr.yellow = 0;
clr.black = 100;
var ofst = 200;
var scl =100;
var minsize = 0.001;
main ();
//tangent circles
function tng(a, b, c, flg) { //a,b,c are circle objects
var k1 = 1 / a.r;
var z1 = new C(a.x, a.y);
var kz1 = Cmul(new C(k1, 0), z1);
var k2 = 1 / b.r;
var z2 = new C(b.x, b.y);
var kz2 = Cmul(new C(k2, 0), z2);
var k3 = 1 / c.r;
var z3 = new C(c.x, c.y);
var kz3 = Cmul(new C(k3, 0), z3);
var k4p = k1 + k2 + k3 + 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
var k4m = k1 + k2 + k3 – 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
var kz4p = Cadd(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
var kz4m = Csub(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
var k4,kz4,k4b,kz4b;
var cs = new Array();
if (k4p>k4m){
k4 = k4p;
kz4 = kz4p;
k4b = k4m;
kz4b = kz4m;
} else {
k4 = k4m;
kz4 = kz4m;
k4b = k4p;
kz4b = kz4p;
}
var cc = new circle(kz4.r/k4,kz4.i/k4,Math.abs(1/k4));
var dx = a.x – cc.x
var dy = a.y – cc.y
var dr = a.r + cc.r
if (Math.abs(dx*dx+dy*dy-dr*dr)>0.0001) {
cc = new circle(kz4b.r/k4,kz4b.i/k4,Math.abs(1/k4));
}
cs.push(cc);
if (flg) {
cc = new circle(kz4b.r/k4b,kz4b.i/k4b,Math.abs(1/k4b));
cs.push(cc);
}
return cs;
}
function main() {
var dc = app.documents.add();
//draw primally circle
var b = new circle(0, 0, -1);
b.draw ();
//generate 2nd and 3rd circles //.54 pi/6
var tr = 1-Math.random()*0.5;
var pa = Math.PI/2 – Math.asin(Math.random()*(1-tr)/tr);
var px = tr * Math.cos(pa);
var py = tr * Math.sin(pa);
var pr = 1 – tr;
var qr = (1 – pr) * (1 – Math.cos(pa + Math.PI/2))
/ (1 + pr – (1 – pr) * Math.cos(pa + Math.PI/2));
var qx = 0;
var qy = qr – 1;
var p = new circle(px, py, pr);
var q = new circle(qx, qy, qr);
p.draw();
q.draw();
var crcl = new Array();
var cs = tng(b, p, q, true);
crcl.push([b, p, cs[0]]);
crcl.push([b, cs[0], q]);
crcl.push([cs[0], p, q]);
crcl.push([b, p, cs[1]]);
crcl.push([b, cs[1], q]);
crcl.push([cs[1], p, q]);
cs[0].draw ();
cs[1].draw();
for (var c=0;c<crcl.length;c++){
var nc = tng(crcl[c][0], crcl[c][1], crcl[c][2])[0];
if (nc.r>minsize){
crcl.push([nc, crcl[c][1], crcl[c][2]]);
crcl.push([crcl[c][0], nc, crcl[c][2]]);
crcl.push([crcl[c][0], crcl[c][1], nc]);
nc.draw ();
}
}
}
function circle(x, y, rad) {
this.x = x;
this.y = y;
this.r = rad;
this.draw = function() {
var cir = app.activeDocument.pathItems.ellipse (
(this.x+this.r)*scl+ofst, (this.y-this.r)*scl+ofst,
this.r*2*scl, this.r*2*scl);
cir.filled = false;
cir.stroked = true;
cir.strokeColor = clr;
cir.strokeWidth = 0.25;
return cir;
}
}
// complex numbers
function C(r,i){
this.r = r; //real
this.i = i; //imaginally
}
//compute complex value
function Cconj(c){//conjunction
return new C(c.r, -c.i);
}
function Cadd(c,d){ //add
return new C(c.r+d.r,c.i+d.i);
}
function Csub(c,d){ //sub
return new C(c.r-d.r,c.i-d.i);
}
function Cmul(c,d){ //multi
return new C(c.r*d.r-c.i*d.i,c.r*d.i+c.i*d.r);
}
function Csqrt(c){ //square ro
ot
return new C(
Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i))*Math.cos(Math.atan2(c.i,c.r)/2),
Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i))*Math.sin(Math.atan2(c.i,c.r)/2)
);
}
こんな感じです。短いコードですが、ステップインで回すと演算の度に飛びまくるので把握しにくいのが難点。