Orbit trapで遊ぼう
どうも、108Hassiumです。
この記事では、フラクタル図形の描画方法の一種である「orbit trap」というアルゴリズムを解説します。
点トラップ
Orbit trapは、数列$${z_{n+1}=f(z_n)}$$と「トラップ」と呼ばれる図形との間の距離の最小値を使った描画法です。
例えば、トラップとして原点を使って$${z^2+c}$$のマンデルブロ集合を描画するコードは以下のようになります。
void setup(){
size(2000,2000);
background(0);
noStroke();
double x,y,px,py,cx,cy,r;
boolean o;
for(int a=0;a<2000;a++){
for(int b=0;b<2000;b++){
cx=(double)a/500.0-2.0;
cy=(double)b/500.0-2.0;
x=0;
y=0;
r=100.0;
o=true;
for(int c=0;c<5000&&o;c++){
px=x;
py=y;
x=px*px-py*py+cx;
y=2.0*px*py+cy;
if(x*x+y*y<r){
r=x*x+y*y;
}
if(10000<x*x+y*y){
o=false;
}
}
fill(sh(r*10.0)*255);
rect(a,b,1,1);
}
}
}
float sh(double x){
return 1.0/((float)(x*x)+1.0);
}※言語はProcessing
描画結果はこちらです。
基本的な構造はマンデルブロ集合を普通に描画するコードと変わりませんが、大きく違うのは真ん中少し下あたりにあるif文です。
...
if(x*x+y*y<r){
r=x*x+y*y;
}
...この部分では$${z_n}$$と原点の距離(の2乗)がdouble型の変数rより小さいか判定し、小さければrの値を更新する、という計算をしています。
※実際に計算に使っているのは距離ではなく距離の2乗なのですが、面倒なので以降ずっと距離と呼ぶことにします。
これにより「$${z_n}$$と原点の距離の最小値」を求め、以下の部分で彩色しています。
...
fill(sh(r*10.0)*255);
rect(a,b,1,1);
...sh(x)という関数はx=0のとき1、x→∞のとき0に漸近していく関数なので、fill(sh(r*100.0)*255)によりrが0に近いほど明るく、0から離れると暗く彩色されます。
ちなみにsh(r*10.0)のrの係数は、小さくすると全体的に明るく描画され、大きくすると暗くなります。
トラップとなる点の位置を変えると、以下のようになります。
※y軸が下向きが正になるような座標系を使用しているため、上下逆に描画されています。
これは$${i}$$をトラップにしたもので、距離計算の部分を以下のように書き換えることで描画できます。
...
if(x*x+(y-1.0)*(y-1.0)<r){
r=x*x+(y-1.0)*(y-1.0);
}
...また、2個以上の点をトラップにすることもできます。
これは$${1}$$と$${1+i}$$をトラップにしたもので、以下のようにrの更新判定を2つ並べるだけで描画できます。
...
if((x-1.0)*(x-1.0)+(y-1.0)*(y-1.0)<r){
r=(x-1.0)*(x-1.0)+(y-1.0)*(y-1.0);
}
if((x-1.0)*(x-1.0)+y*y<r){
r=(x-1.0)*(x-1.0)+y*y;
}
...更に、rが更新されたときのcの値($${z_n}$$の$${n}$$)を記録しておくことで、模様に色を付けることができるようになります。
void setup(){
size(2000,2000);
background(0);
noStroke();
double x,y,px,py,cx,cy,r;
int t;
boolean o;
for(int a=0;a<2000;a++){
for(int b=0;b<2000;b++){
cx=(double)a/500.0-2.0;
cy=(double)b/500.0-2.0;
x=0;
y=0;
r=100.0;
t=0;
o=true;
for(int c=0;c<5000&&o;c++){
px=x;
py=y;
x=px*px-py*py+cx;
y=2.0*px*py+cy;
if(x*x+y*y<r){
r=x*x+y*y;
t=c;
}
if(10000<x*x+y*y){
o=false;
}
}
fill(cr(t*99)*sh(r*10.0),cr(t*99+85)*sh(r*10.0),cr(t*99+170)*sh(r*10.0));
rect(a,b,1,1);
}
}
}
float sh(double x){
return 1.0/((float)(x*x)+1.0);
}
float cr(float c){
return (c%256)*(256-(c%256))/65;
}
さて、ここまでの話は$${z^2+c}$$以外の関数のマンデルブロ集合やジュリア集合にもほぼそのまま応用できます。
ジュリア集合に使う場合、rの計算をする位置を$${z_{n+1}}$$を計算する前にずらしてください。
...
if(x*x+y*y<r){
r=x*x+y*y;
}
px=x;
py=y;
x=px*px-py*py+cx
y=2.0*px*py+cy
...
なお、点をトラップにしてジュリア集合を描画する場合、パラメータの組み合わせ次第では生理的嫌悪感を催させるような仕上がりになることがあります。
この記事の執筆中に、運悪くどこにも公開できないレベルのグロ画像が生成される組み合わせを発見してしまったので、ソースコードだけ載せておきます。
どうしても気になるという人は自己責任で実行してみてください。
void setup(){
size(2000,2000);
background(0);
noStroke();
double x,y,px,py,cx,cy,r;
boolean o;
for(int a=0;a<2000;a++){
for(int b=0;b<2000;b++){
x=(double)a/500.0-2.0;
y=(double)b/500.0-2.0;
cx=0.3;
cy=0.5;
r=100.0;
o=true;
for(int c=0;c<5000&&o;c++){
if(x*x+y*y<r){
r=x*x+y*y;
}
px=x;
py=y;
x=px*px-py*py+cx;
y=2.0*px*py+cy;
if(10000<x*x+y*y){
o=false;
}
}
fill(sh(r*100.0)*255);
rect(a,b,1,1);
}
}
}
float sh(double x){
return 1.0/((float)(x*x)+1.0);
}線トラップ
点だけでなく、線をトラップとして使うこともできます。
例えば、座標平面上の点$${(x,y)}$$と$${x}$$軸との間の距離の2乗は$${y^2}$$なので、$${x}$$軸(複素平面では実軸)をトラップとするorbit trapは以下のように実装できます。
void setup(){
size(2000,2000);
background(0);
noStroke();
double x,y,px,py,cx,cy,r;
boolean o;
for(int a=0;a<2000;a++){
for(int b=0;b<2000;b++){
cx=(double)a/500.0-2.0;
cy=(double)b/500.0-2.0;
x=0;
y=0;
r=100.0;
o=true;
for(int c=0;c<5000&&o;c++){
px=x;
py=y;
x=px*px-py*py+cx;
y=2.0*px*py+cy;
if(y*y<r){ //ここが違う
r=y*y;
}
if(10000<x*x+y*y){
o=false;
}
}
fill(sh(r*10.0)*255);
rect(a,b,1,1);
}
}
}
float sh(double x){
return 1.0/((float)(x*x)+1.0);
}実行結果はこうなります。
直線以外の図形に対して距離を計算するのは困難ですが、距離と似た性質を持つ「距離もどき」なら計算できます。
$${f(x,y)=0}$$という式で表せる曲線に対し、$${f(x,y)^2}$$の値は
$${f(x,y)=0}$$ならば$${f(x,y)^2=0}$$
$${f(x,y)\neq0}$$ならば$${0<f(x,y)^2}$$
・・・という性質を満たすので、距離の代わりとして機能することがあります。(複雑な式で表される図形だと上手くいかないことがあります)
例えば$${f(x,y)=y-x^2+1}$$とすると、$${y=x^2-1}$$で表される放物線をトラップにすることができます。
※y軸が下向きが正になるような座標系を使用しているため、上下逆に描画されています。
rを計算する部分は以下のようになります。
...
if((y-x*x+1.0)*(y-x*x+1.0)<r){
r=(y-x*x+1.0)*(y-x*x+1.0);
}
...他の図形も試してみます。
ジュリア集合でもきちんと動作します。
点トラップと同じように、複数の図形をトラップにすることもできます。
$${f(x,y)g(x,y)=0}$$という方程式で表される図形は$${f(x,y)=0}$$と$${g(x,y)=0}$$のグラフを重ね合わせたものになるのですが、どうやら「$${f(x,y)=0}$$と$${g(x,y)=0}$$をトラップにした場合」と「$${f(x,y)g(x,y)=0}$$をトラップにした場合」では微妙に異なった結果になるようです。
並べて比較してみると、線の交差の仕方が異なっています。
z^2+cのマンデルブロ集合
交差の仕方といえば、orbit trapで画像検索したときに出てくる画像には、以下のように線が交差せずに立体的に重なっているような見た目のものが見られます。
これはrの更新判定を以下のように書き換えることで描画できます。
...
if((x*x+y*y-1.0)*(x*x+y*y-1.0)<0.01){
r=(x*x+y*y-1.0)*(x*x+y*y-1.0);
o=false;
}
...また、rを更新した後のo=false;を削除すると、線の重なる順序が逆になります。
応用
応用です。
Orbit portal
Orbit portalは、不等式で表される領域をトラップとして使う手法です。(一応私が自力で発見した手法ですが、もしかしたら前例があって違う名前がついているかもしれません)
void setup(){
size(2000,2000);
background(0);
noStroke();
double x,y,px,py,cx,cy;
boolean o;
for(int a=0;a<2000;a++){
for(int b=0;b<2000;b++){
cx=(double)a/500.0-2.0;
cy=(double)b/500.0-2.0;
x=0;
y=0;
o=true;
for(int c=0;c<5000&&o;c++){
px=x;
py=y;
x=px*px-py*py+cx;
y=2.0*px*py+cy;
if(x*x+(y-3.0)*(y-3.0)<1.0){
cx=x*2.0;
cy=(y-3.0)*2.0;
x=0;
y=0;
}
if(10000<x*x+y*y){
fill(cr(c*7),cr(c*8),cr(c*9));
rect(a,b,1,1);
o=false;
}
}
}
}
}
float cr(float c){
return (c%256)*(256-(c%256))/65;
}$${z_n}$$がトラップの中に入ったら$${c}$$の値を再設定し、$${z_n}$$も0にリセットして計算し直す、という感じです。
ジュリア集合の場合は以下のようになります。
void setup(){
size(2000,2000);
background(0);
noStroke();
double x,y,px,py,cx,cy;
boolean o;
for(int a=0;a<2000;a++){
for(int b=0;b<2000;b++){
x=(double)a/500.0-2.0;
y=(double)b/500.0-2.0;
cx=-0.71;
cy=0.24;
o=true;
for(int c=0;c<5000&&o;c++){
px=x;
py=y;
x=px*px-py*py+cx;
y=2.0*px*py+cy;
if(x*x+(y-2.0)*(y-2.0)<1.0){
x=x*2.0;
y=(y-2.0)*2.0;
}
if(10000<x*x+y*y){
fill(cr(c*7),cr(c*8),cr(c*9));
rect(a,b,1,1);
o=false;
}
}
}
}
}
float cr(float c){
return (c%256)*(256-(c%256))/65;
}トラップに入った時の処理が変わっています。
マンデルブロ集合のときと同じように$${c}$$を書き換える処理に変更した場合、「マンデルブロ集合に片思いしまくるジュリア集合」みたいな見た目になります。
逆にマンデルブロ集合を描画する際に$${z_n}$$を更新する方の処理をすると、こんな感じになります。
Cellular coloring
Cellular coloringは、以下のような見た目の画像を描画できる手法です。
初出はおそらく以下のページです。
https://fractalforums.org/programming/11/cellular-coloring-of-mandelbrot-insides/3264
具体的なコードは載っていないのですが、断片的な説明の内容からおそらく以下のようなコードだろうと推測しました。
void setup(){
size(2000,2000);
background(0);
noStroke();
int t;
double x,y,px,py,cx,cy,r,pr;
boolean o;
for(int a=0;a<2000;a++){
for(int b=0;b<2000;b++){
cx=(double)a/500.0-2.0;
cy=(double)b/500.0-2.0;
x=0;
y=0;
pr=100.0;
r=100.0;
o=true;
t=0;
for(int c=0;c<5000&&o;c++){
px=x;
py=y;
x=px*px-py*py+cx;
y=2.0*px*py+cy;
if(x*x+y*y<r){
pr=r;
r=x*x+y*y;
t=c;
}else if(x*x+y*y<pr){
pr=x*x+y*y;
}
if(10000<x*x+y*y){
o=false;
}
}
fill(cr(t*99)*(1.0-sh(10.0*(1.0-r/pr))),cr(t*99+85)*(1.0-sh(10.0*(1.0-r/pr))),cr(t*99+170)*(1.0-sh(10.0*(1.0-r/pr))));
rect(a,b,1,1);
}
}
}
float sh(double x){
return 1.0/((float)(x*x)+1.0);
}
float cr(float c){
return (c%256)*(256-(c%256))/65;
}rが更新されたタイミングを記録するtのほか、$${z_n}$$とトラップの距離のうち2番目に小さいものを記録するprという変数が追加されています。
tを使って彩色することでタイル状の模様が現れるのはすでに見せた通りですが、どうやらタイルの境界線上ではr=prになり、それを利用することでフチの部分の色を変えられるようです。(詳しい仕組みは私にもよくわかりません)
トラップを複数個にしたり線をトラップにするといった複雑な手法とは相性が悪いようですが、他の関数のマンデルブロ集合やジュリア集合を描画するのは問題ないようです。
ジュリア集合の場合、周期サイクルへの収束のスピードによっては影が濃くなってかなり独特な見た目になります。
$${z^2+c}$$以外のマンデルブロ集合はこんな感じです。
追記
遊び足りない気がしたので続きを書きました。