Helecho de Barnsley

El helecho de Barnsley usa cuatro transformaciones afines. La fórmula de una transformación es la siguiente:

${\displaystyle f(x,y)={\begin{bmatrix}a&b\\c&d\end{bmatrix}}{\begin{bmatrix}x\\y\end{bmatrix}}+{\begin{bmatrix}e\\f\end{bmatrix}}}$

Barnsley muestra el código "IFS" para su fractal de helecho "Black Spleenwort" como una matriz de valores que se muestran en una tabla.^[3] En la tabla, las columnas de la "a" a la "f" son los coeficientes de la ecuación y "p" representa el factor de probabilidad.

Estos corresponden a las siguientes transformaciones:

${\displaystyle f_{1}(x,y)={\begin{bmatrix}\ 0.00&\ 0.00\ \\0.00&\ 0.16\end{bmatrix}}{\begin{bmatrix}\ x\\y\end{bmatrix}}}$

${\displaystyle f_{2}(x,y)={\begin{bmatrix}\ 0.85&\ 0.04\ \\-0.04&\ 0.85\end{bmatrix}}{\begin{bmatrix}\ x\\y\end{bmatrix}}+{\begin{bmatrix}\ 0.00\\1.60\end{bmatrix}}}$

${\displaystyle f_{3}(x,y)={\begin{bmatrix}\ 0.20&\ -0.26\ \\0.23&\ 0.22\end{bmatrix}}{\begin{bmatrix}\ x\\y\end{bmatrix}}+{\begin{bmatrix}\ 0.00\\1.60\end{bmatrix}}}$

${\displaystyle f_{4}(x,y)={\begin{bmatrix}\ -0.15&\ 0.28\ \\0.26&\ 0.24\end{bmatrix}}{\begin{bmatrix}\ x\\y\end{bmatrix}}+{\begin{bmatrix}\ 0.00\\0.44\end{bmatrix}}}$

Aunque, en teoría, el helecho de Barnsley podría trazarse a mano con un bolígrafo y papel cuadriculado, el número de iteraciones necesarias asciende a decenas de miles, lo que hace que el uso de una computadora sea prácticamente obligatorio. Muchos modelos informáticos diferentes del helecho de Barnsley son populares entre los matemáticos contemporáneos. Siempre que las matemáticas se programen correctamente utilizando la matriz de constantes de Barnsley, se producirá la misma forma de helecho.

El primer punto dibujado está en el origen (x₀ = 0, y₀ = 0) y a continuación los nuevos puntos se calculan iterativamente aplicando aleatoriamente una de las siguientes cuatro transformaciones de coordenadas:^[4][5]

ƒ₁

x_{n + 1} = 0

y_{n + 1} = 0.16 y_n.

Esta transformación de coordenadas se elige el 1% del tiempo y simplemente asigna cualquier punto a un punto en el primer segmento de línea en la base del tallo. Esta parte de la figura es la primera que se completa durante el transcurso de las iteraciones.

ƒ₂

x_{n + 1} = 0.85 x_n + 0.04 y_n

y_{n + 1} = −0.04 x_n + 0.85 y_n + 1.6.

Esta transformación de coordenadas se elige el 85% de las veces y calcula cualquier punto dentro del foliolo representado por el triángulo rojo a un punto dentro del folliolo opuesto, más pequeño, representado por el triángulo azul en la figura.

ƒ₃

x_{n + 1} = 0.2 x_n - 0.26 y_n

y_{n + 1} = 0.23 x_n + 0.22 y_n + 1.6.

Esta transformación de coordenadas se elige el 7% del tiempo y calcula cualquier punto dentro del foliolo representado por el triángulo azul a un punto dentro del triángulo correspondiente alterno a través del tallo (lo voltea).

ƒ₄

x_{n + 1} = −0.15 x_n + 0.28 y_n

y_{n + 1} = 0.26 x_n + 0.24 y_n + 0.44.

Esta transformación de coordenadas se elige el 7% de las veces y calcula cualquier punto dentro del foliolo representado por el triángulo azul a un punto dentro del triángulo correspondiente alterno a través del tallo (sin voltearlo).

La primera transformación de coordenadas dibuja el tallo. La segunda genera copias sucesivas del tallo y las frondas inferiores para formar el helecho completo. La tercera dibuja la fronda inferior a la izquierda. El cuarto dibuja la fronda inferior a la derecha. La naturaleza recursiva del SIF garantiza que el conjunto sea una réplica más grande de cada fronda. Debe tenerse en cuenta que el helecho completo está dentro del rango −2,1820 < x < 2.6558 y 0 ≤ y < 9.9983.

Jugando con los coeficientes, es posible crear variedades de helechos mutantes. En su artículo sobre fractales V-variables, Barnsley llama a este rasgo un "superfractal".^[2]

Sin embargo, un experimentador ha elaborado una tabla de coeficientes para producir otro helecho de aspecto notablemente natural, parecido a los helechos Cyclosorus o Thelypteridaceae. Estos valores paramétricos son:^[6][7]

Puede usarse el código siguiente para dibujar el helecho con un ordenador:

import turtle
import random

pen = turtle.Turtle()
pen.speed(0)
pen.color("green")
pen.penup()

x = 0
y = 0
for n in range(11000):
    pen.goto(65 * x, 37 * y - 252)  # scale the fern to fit nicely inside the window
    pen.pendown()
    pen.dot(3)
    pen.penup()
    r = random.random()
    if r < 0.01:
        x, y =  0.00 * x + 0.00 * y,  0.00 * x + 0.16 * y + 0.00
    elif r < 0.86:
        x, y =  0.85 * x + 0.04 * y, -0.04 * x + 0.85 * y + 1.60
    elif r < 0.93:
        x, y =  0.20 * x - 0.26 * y,  0.23 * x + 0.22 * y + 1.60
    else:
        x, y = -0.15 * x + 0.28 * y,  0.26 * x + 0.24 * y + 0.44

# Barnsley's Fern

# create function of the probability and the current point
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
}

# how many reps determines how detailed the fern will be
reps <- 10000

# create a vector with probability values, and a matrix to store coordinates
p <- runif(reps)

# initialise a point at the origin
coords <- c(0, 0)

# compute Fractal Coordinates
m <- Reduce(fractal_fern2, p, accumulate = T, init = coords)
m <- t(do.call(cbind, m))

# Create plot
plot(m, type = "p", cex = 0.1, col = "darkgreen",
     xlim = c(-3, 3), ylim = c(0, 10), 
     xlab = NA, ylab = NA, axes = FALSE)

/*
    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();
  }
}

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();
  }
}

<canvas id="canvas" height="700" width="700">
</canvas>

<script>
    let canvas;
    let canvasContext;

    let x = 0, y = 0;

    window.onload = function () {
        canvas = document.getElementById("canvas");
        canvasContext = canvas.getContext('2d');

        canvasContext.fillStyle = "black";
        canvasContext.fillRect(0, 0, canvas.width, canvas.height);

        setInterval(() => {
            // Update 20 times every frame
            for (let i = 0; i < 20; i++)
                update();
                
        }, 1000/250); // 250 frames per second
    };

    function update() {

        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, 1, "green");

        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>

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)

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

/* results table */
declare @fern table (Fun int, X float, Y float, Seq int identity(1,1) primary key, DateAdded datetime default getdate())

declare @i int = 1	/* iterations */
declare @fun int	/* random function */

declare @x float = 0	/* initialise X = 0 */
declare @y float = 0	/* initialise Y = 0 */
declare @rand float
insert into @fern (Fun, X, Y) values (0,0,0)	/* set starting point */

while @i < 5000 /* how many points? */
begin

	set @rand = rand()

	select @Fun = case	/* get random function to use -- @fun = f1 = 1%, f2 = 85%, f3 = 7%, f4 = 7% */
		when @rand <= 0.01 then 1
		when @rand <= 0.86 then 2
		when @rand <= 0.93 then 3
		when @rand <= 1 then 4
	end

	select top 1 @X = X, @Y = Y from @fern order by Seq desc /* get previous point */

	insert into @fern(Fun, X, Y)	/* transform using four different function expressions  */
	select @fun,
		case @fun
			when 1 then 0
			when 2 then 0.85*@x+0.04*@y
			when 3 then 0.2*@x-0.26*@y
			when 4 then -0.15*@x + 0.28*@y
		end X,
		case @fun
			when 1 then 0.16*@y
			when 2 then -0.04*@x + 0.85*@y + 1.6
			when 3 then 0.23*@x + 0.22*@y + 1.6
			when 4 then 0.26*@x + 0.24*@y + 0.44
		end Y

	set @i=@i+1
end

select top 5000 *,geography::Point(Y, X, 4326) from @fern 
order by newid()

N = 1000000;

xy = [0; 0];
fern = zeros(N, 2);

f_1 = [0 0; 0 0.16];
f_2 = [0.85 0.04; -0.04 0.85];
f_3 = [0.2 -0.26; 0.23 0.22];
f_4 = [-0.15 0.28; 0.26 0.24];

P = randsample(1:4, N, true, [0.01 0.85 0.07 0.07]);
for i = 2:N
    p = P(i - 1);
    if p == 1 % Stem
        xy = f_1 * xy;
    elseif p == 2 % Sub-leaflets
        xy = f_2 * xy + [0; 1.6];
    elseif p == 3 % Left leaflet
        xy = f_3 * xy + [0; 1.6];
    else % Right leaflet
        xy = f_4 * xy + [0; 0.44];
    end
    
    fern(i, 1) = xy(1);
    fern(i, 2) = xy(2);
end
clearvars -except N fern % R2008a+

% Plotting the fern

%{
% Better detail, slower performance
c = linspace(0, 0.35, N - 1); c(end + 1) = 1;
colormap(summer(N));
set(gcf, 'Color', 'k', 'position', [10, 50, 800, 600]);
scatter(fern(:, 1), fern(:, 2), 0.1, c, 'o');
set(gca, 'Color', 'k');
%}

%
% Less detail, better performance
c = linspace(1, 0.2, N - 1); c(end + 1) = 0;
colormap(summer(N));
set(gcf, 'Color', 'k', 'position', [10, 50, 800, 600]);
scatter(fern(:, 1), fern(:, 2), 0.1, c, '.');
set(gca, 'Color', 'k');
%}

BarnsleyFern[iterations_] := (
  it = iterations; 
  noElement = 1;
  BarnsleyfernTransformation(*obtained from the affine transformation matrix for simpler
coding*)= {
    {0, +0.16 y},
    {0.85 x + 0.04 y, 1.6 - 0.04 x + 0.85 y},
    {0.2 x - 0.26 y, 1.6 + 0.23 x + 0.22 y},
    {-0.15 x + 0.28 y, 0.44 + 0.26 x + 0.24 y}};
   (*list of coordinates*)
  point = {{0, 0}};
  While[it > 0 && noElement <= iterations,
   AppendTo[point,
    BarnsleyfernTransformation[[(*determines which transformation \
applies to the previous point*)
       RandomInteger[{1, 4}]]] /. {x -> Part[point, noElement, 1], 
      y -> Part[point, noElement, 2] }];
   (*update var for the while loop*)
   it -= 1; noElement += 1];
  ListPlot[point, AspectRatio -> Automatic, PlotRange -> Full])