バーンズリーのシダ

バーンズリーのシダ（Barnsley fern）は、イギリスの数学者にちなんで名付けられたフラクタルで、その著書『"Fractals Everywhere』"で初めて記述された。チャセンシダの1種であるミヤマカワリシダ（）を模したものである。
歴史.
シダは、自己相似図形の代表例の1つである。すなわち、数学的に生成されたパターンであり、そのパターンは任意に拡大または縮小しても同一の形状が現れる自己相似性を示す。シェルピンスキーの三角形のように、バーンズリーのシダは、コンピューターで数式を繰り返し使用することで、視覚的に美しい図形を構築する方法を示している。バーンズリーの1988年の本『"Fractals Everywhere"』は、 ジョージア工科大学数学部の"フラクタル幾何学"と呼ばれるコースにおける学部生と大学院生に対する講義内容に基づいている。この本の出版後、"フラクタル測定理論"と呼ばれる第二のコースが設置された。バーンズリーの作品は、数学モデルで自然を模倣しようとするグラフィックアーティストのインスピレーションの源となっている。
バーンズリーによって開発されたシダのコードは、フラクタルを作成するための反復関数系（IFS：Iterated Function System）の例である。これはに従う。彼はフラクタルを使用して、科学技術におけるさまざまな現象、とりわけ植物構造、についてモデル化している。 —マイケル・バーンズリー"ら" 
シダの生成.
バーンズリーのシダは、4つのアフィン変換を使用する。変換の式は次のとおりである。
あるいは
バーンズリーは、前述した著書で、シダを生成する"IFS"コードを、以下の表のように、値の行列の形で示した。下表では、列「a」から「f」は前記の等式f(x,y)の係数を示しており、「p」はそれぞれの値の組に対する確率因子を表している。
上記の表は、次の4つの変換に対応している。
コンピューターによる生成.
バーンズリーのシダは、理論的にはペンとグラフ用紙を使用して手でプロットできるが、必要な反復回数は数万回に達するので、実際にはコンピューターを使用する必要がある。バーンズリーのシダの多くのさまざまな計算モデルは、現代の数学者に人気がある。バーンズリーの係数行列を使用し、数式が正しくプログラムされている限り、同じ形状をもったシダが生成される。
最初に描画される点は原点（ "x" 0 = 0、 "y" 0 = 0）にあり、次の4つの座標変換ƒ1, ƒ2, ƒ3, ƒ4のいずれかを設定した発生確率pでランダムに適用することにより、新しい点を反復的に計算する 。すなわち、どの座標変換を使用するか確率に応じた乱数範囲を設定しておき、コンピューターで乱数を発生させ、乱数の値に応じてどの座標変換を使用するのか決め、変換計算を反復的に実行する。これは、1988年にバーンズリーがランダム・アルゴリズムと名付けた初めての例であり、現在ではカオスゲーム、確率論的反復関数系とも呼ばれている。
"ƒ"1 
座標変換ƒ1は1％の確率(p=0.01)で実行される。この変換は、任意の点を、シダの茎基部の最初の線分にマッピングするだけである。フラクタル図のこの部分は、反復の過程で最初に完了する部分である。 
"ƒ"2
この座標変換ƒ2は、85％の確率(p=0.85)で実行され、図の赤い三角形で表される葉の中の任意のポイントを、青い三角形で表される反対側の小さな葉の中の点にマッピングする。
"ƒ"3
この座標変換ƒ3は、7％の確率(p=0.07)で実行され、青い三角形で表される葉（または羽片）の中の任意の点を、茎を交差し交互に対応する三角形内にマッピングする（反転する）。
"ƒ"4
この座標変換ƒ4は、7％の確率(p=0.07)で実行され、青い三角形で表される葉（または羽片）の中の任意の点を、茎を交差し交互に対応する三角形内にマッピングする（反転せずに）。
最初の座標変換は茎を描画する。2番目は茎と下葉を連続コピーして、完全なシダを生成する。3番目は左下葉を描き、4番目は、右下の葉を描く。IFSの再帰的な性質により、全体が各葉の部分の大きな複製となっている。シダ全体が-2.1820 < x < 2.6558 および 0 ≤ "y" <   9.9983　の範囲内にあることに注意。
様々なシダの変種.
係数をいろいろと変えることで、様々なシダの品種、変種を作成することが可能である。V変数フラクタルに関する彼の論文において、バーンズリーはこの特性を"スーパー"フラクタルと呼んでいる。
ある実験者は係数の表を作成して、'または'シダに類似した、非常に自然に見える別のシダを作成した（下表に示す） 。
構文例.
以下の各プログラム言語の構文を使用して、シダを自分で描くことができる。
Python.
import turtle
import random
pen = turtle.Turtle()
pen.speed(20)
pen.color('green')
pen.penup()
x = 0
y = 0
for n in range(100000):
 pen.goto(85*x, 57*y-275) # 57 is to scale the fern and -275 is to start the drawing from the bottom.
 pen.pendown()
 pen.dot()
 pen.penup()
 r = random.random() # to get probability
 r = r * 100
 xn = x
 yn = y
 if r < 1: # elif ladder based on the probability
 x = 0
 y = 0.16 * yn
 elif r < 86:
 x = 0.85 * xn + 0.04 * yn
 y = -0.04 * xn + 0.85 * yn + 1.6
 elif r < 93:
 x = 0.20 * xn - 0.26 * yn
 y = 0.23 * xn + 0.22 * yn + 1.6
 else:
 x = -0.15 * xn + 0.28 * yn
 y = 0.26 * xn + 0.24 * yn + 0.44
R.
fractal_fern2 <- function(x, p){
 if (p <= 0.01) {
 m <- matrix(c(0, 0, 0, .16), 2, 2)
 f <- c(0, 0)
 } else if (p <= 0.86) {
 m <- matrix(c(.85, -.04, .04, .85), 2, 2)
 f <- c(0, 1.6)
 } else if (p <= 0.93) {
 m <- matrix(c(.2, .23, -.26, .22), 2, 2)
 f <- c(0, 1.6)
 } else {
 m <- matrix(c(-.15, .26, .28, .24), 2, 2)
 f <- c(0, .44)
 m %*% x + f
reps <- 10000
p <- runif(reps)
coords <- c(0, 0)
m <- Reduce(fractal_fern2, p, accumulate = T, init = coords)
m <- t(do.call(cbind, m))
plot(m, type = "p", cex = 0.1, col = "darkgreen",
 xlim = c(-3, 3), ylim = c(0, 10), 
 xlab = NA, ylab = NA, axes = FALSE)
Processing.
 Barnsley Fern for Processing 3.4
// declaring variables x and y
float x, y;
// creating canvas
void setup() {
 size(600, 600);
 background(255);
/* setting stroke, mapping canvas and then
 plotting the points */
void drawPoint() {
 stroke(34, 139, 34);
 strokeWeight(1);
 float px = map(x, -2.1820, 2.6558, 0, width);
 float py = map(y, 0, 9.9983, height, 0);
 point(px, py);
/* algorithm for calculating value of (n+1)th
 term of x and y based on the transformation
 matrices */
void nextPoint() {
 float nextX, nextY;
 float r = random(1);
 if (r < 0.01) {
 nextX = 0;
 nextY = 0.16 * y;
 } else if (r < 0.86) {
 nextX = 0.85 * x + 0.04 * y;
 nextY = -0.04 * x + 0.85 * y + 1.6;
 } else if (r < 0.93) {
 nextX = 0.20 * x - 0.26 * y;
 nextY = 0.23 * x + 0.22 * y + 1.6;
 } else {
 nextX = -0.15 * x + 0.28 * y;
 nextY = 0.26 * x + 0.24 * y + 0.44;
 x = nextX;
 y = nextY;
/* iterate the plotting and calculation
 functions over a loop */
void draw() {
 for (int i = 0; i < 100; i++) {
 drawPoint();
 nextPoint();
P5.JS.
let x = 0;
let y = 0;
function setup() {
 createCanvas(600, 600);
 background(0);
//range −2.1820 < x < 2.6558 and 0 ≤ y < 9.9983.
function drawPoint() {
 stroke(255);
 strokeWeight(1);
 let px = map(x, -2.1820, 2.6558, 0, width);
 let py = map(y, 0, 9.9983, height, 0);
 point(px, py);
function nextPoint() {
 let nextX;
 let nextY;
 let r = random(1);
 if (r < 0.01) {
 //1
 nextX = 0;
 nextY = 0.16 * y;
 } else if (r < 0.86) {
 //2
 nextX = 0.85 * x + 0.04 * y;
 nextY = -0.04 * x + 0.85 * y + 1.60;
 } else if (r < 0.93) {
 //3
 nextX = 0.20 * x + -0.26 * y;
 nextY = 0.23 * x + 0.22 * y + 1.60;
 } else {
 //4
 nextX = -0.15 * x + 0.28 * y;
 nextY = 0.26 * x + 0.24 * y + 0.44;
 x = nextX;
 y = nextY;
function draw() {
 for (let i = 0; i < 1000; i++) {
 drawPoint();
 nextPoint();
} 
JavaScript (HTML5).
<canvas id="canvas" height="700" width="700">
</canvas>
<script>
 let canvas;
 let canvasContext;
 let move;
 let fps = 250;
 let pointColor = "green";
 let backgroundColor = "black";
 let pointRadius = 1;
 let x = 0;
 let y = 0;
 window.onload = function () {
 canvas = document.getElementById("canvas");
 canvasContext = canvas.getContext('2d');
 canvasContext.fillStyle = backgroundColor;
 canvasContext.fillRect(0, 0, canvas.width, canvas.height);
 move = setInterval(doIt, 1000/fps);
 let doIt = function () {
 for (let i = 0; i < 20; i++)
 moveAll()
 function moveAll() {
 let nextX, nextY;
 let r = Math.random();
 if (r < 0.01) {
 nextX = 0;
 nextY = 0.16 * y;
 } else if (r < 0.86) {
 nextX = 0.85 * x + 0.04 * y;
 nextY = -0.04 * x + 0.85 * y + 1.6;
 } else if (r < 0.93) {
 nextX = 0.20 * x - 0.26 * y;
 nextY = 0.23 * x + 0.22 * y + 1.6;
 } else {
 nextX = -0.15 * x + 0.28 * y;
 nextY = 0.26 * x + 0.24 * y + 0.44;
 // Scaling and positioning
 let plotX = canvas.width * (x + 3) / 6;
 let plotY = canvas.height - canvas.height * ((y + 2) / 14);
 drawFilledCircle(plotX, plotY, pointRadius, pointColor);
 x = nextX;
 y = nextY;
 const drawFilledCircle = (centerX, centerY, radius, color) => {
 canvasContext.beginPath();
 canvasContext.fillStyle = color;
 canvasContext.arc(centerX, centerY, radius, 0, 2 * Math.PI, true);
 canvasContext.fill();
</script>
QBasic.
SCREEN 12
WINDOW (-5, 0)-(5, 10)
RANDOMIZE TIMER
COLOR 10
DO
 SELECT CASE RND
 CASE IS < .01
 nextX = 0
 nextY = .16 * y
 CASE .01 TO .08
 nextX = .2 * x - .26 * y
 nextY = .23 * x + .22 * y + 1.6
 CASE .08 TO .15
 nextX = -.15 * x + .28 * y
 nextY = .26 * x + .24 * y + .44
 CASE ELSE
 nextX = .85 * x + .04 * y
 nextY = -.04 * x + .85 * y + 1.6
 END SELECT
 x = nextX
 y = nextY
 PSET (x, y)
LOOP UNTIL INKEY$ = CHR$(27)
VBA（CorelDraw）.
Sub Barnsley()
Dim iEnd As Long
Dim i As Long
Dim x As Double
Dim y As Double
Dim nextX As Double
Dim nextY As Double
Dim sShapeArray() As Shape
Dim dSize As Double
Dim sColor As String
 dSize = 0.01 ' Size of the dots
 sColor = "0,0,100" ' RGB color of dots, value range 0 to 255
 iEnd = 5000 ' Number of iterations
 ReDim sShapeArray(iEnd)
 ' In Corel, each object drawn requires a variable name of its own
 Randomize ' Initialize the Rnd function
 For i = 0 To iEnd ' Iterate ...
 Select Case Rnd
 Case Is < 0.01
 ' f1 = Draw stem
 nextX = 0
 nextY = 0.16 * y
 Case 0.01 To 0.08
 ' f3
 nextX = 0.2 * x - 0.26 * y
 nextY = 0.23 * x + 0.22 * y + 1.6
 Case 0.08 To 0.15
 ' f4
 nextX = -0.15 * x + 0.28 * y
 nextY = 0.26 * x + 0.24 * y + 0.44
 Case Else
 ' f2
 nextX = 0.85 * x + 0.04 * y
 nextY = -0.04 * x + 0.85 * y + 1.6
 End Select
 x = nextX
 y = nextY
 Set sShapeArray(i) = ActiveLayer.CreateEllipse2(x + 2.5, y + 0.5, dSize)
 sShapeArray(i).Style.StringAssign "{""fill"":{""primaryColor"":""RGB255,USER," & sColor & ",100,00000000-0000-0000-0000-000000000000"",""secondaryColor"":""RGB255,USER,255,255,255,100,00000000-0000-0000-0000-000000000000"",""type"":""1"",""fillName"":null},""outline"":{""width"":""0"",""color"":""RGB255,USER,0,0,0,100,00000000-0000-0000-0000-000000000000""},""transparency"":{}}"
 DoEvents
 Next
 End Sub
Amola.
addpackage("Forms.dll")
set("x", 0)
set("y", 0)
set("width", 600)
set("height", 600)
method setup()
 createCanvas(width, height)
 rect(0, 0, 600, 600, color(0, 0, 0))
end
method drawPoint()
 set("curX", div(mult(width, add(x, 3)), 6))
 set("curY", sub(height, mult(height, div(add(y, 2), 14))))
 set("size", 1)
 //log(curX)
 //log(curY)
 rect(round(curX - size / 2), round(curY - size / 2), round(curX + size / 2), round(curY + size / 2), color(34, 139, 34))
end
method nextPoint()
 set("nextX", 0)
 set("nextY", 0)
 set("random", random(0, 100))
 if(random < 1)
 set("nextX", 0)
 set("nextY", 0.16 * y)
 end
 else
 if(random < 86)
 set("nextX", 0.85 * x + 0.04 * y)
 set("nextY", -0.04 * x + 0.85 * y + 1.6)
 end
 else
 if(random < 93)
 set("nextX", 0.2 * x - 0.26 * y)
 set("nextY", 0.23 * x + 0.22 * y + 1.6)
 end
 else
 set("nextX", -0.15 * x + 0.28 * y)
 set("nextY", 0.26 * x + 0.24 * y + 0.44)
 end
 end
 end
 set("x", nextX)
 set("y", nextY)
end
setup()
while(true)
 drawPoint()
 nextPoint()
end
表計算ソフトによる計算例.
プログラム言語の構文と同じ計算を、以下のように表計算ソフトの関数を使用することでも実現できる。
なお、B8,C8,D8のセルには以下のような複数条件判定の関数(ネスティング参照)を入力する。
最終8行目をオートフィルで適当な行数だけコピーし、XY散布図とするとバーンズリーのシダのフラクタル図が得られる。各変換式ƒの係数a,b,c,d,e,fと確率pは任意に変更可能である。
各列は以下のような計算を行っている。