От треугольника Серпинского до «снежинок». Реализовано на Python

Здравствуй, читатель. Это моя первая статья на Habr. В этой статье я поделюсь своим опытом построения треугольника Серпинского и расскажу, к чему привели дальнейшие эксперименты с фракталами подобного типа. Заранее прошу сильно не бить, я не являюсь специалистом в области фракталов, а код был написан для себя и особо не причесывался.

Но начнем сначала... Это было мое первое знакомство с фракталами, и начальной задачей было построение треугольника Серпинского. Строится данный фрактал согласно нескольким правилам:

• задаются координаты вершин треугольника;

• задаются координаты первой точки фрактала. Тут имеется несколько вариантов: либо выбирается любая вершина треугольника, либо задается точка со случайными координатами, но внутри треугольника;

• каждая следующая точка треугольника Серпинского находится на отрезке, соединяющем любую имеющуюся точку фрактала со случайной вершиной треугольника. В простом случае данная точка лежит на середине описанного отрезка.

В конечном счете должно получится подобное изображение.

Код в Python имеет данный вид. В представленном варианте треугольник Серпинского состоит из 10000 точек.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root

#Вершины треугольника
A=[0,0]
B=[0.5,np.sqrt(3)/2]
C=[1,0]
ABC=[A,B,C]

#Случайная точка внутри треугольника
fx=np.random.randint(1,100)/100
if fx<=0.5:
    fy=np.random.randint(1,fx*np.tan(np.pi/3)*100)/100
else:
    fy=np.random.randint(1,(1-fx)*np.tan(np.pi/3)*100)/100
p=[[fx,fy]]

#Генерирую остальные точки согласно правилу: каждая точка располагается на середине
#отрезка, соединяющего любую имеющуюся точку с любой вершиной треугольника.

for _ in range(9999):
    U=ABC[np.random.randint(0,len(ABC))]
    G=p[np.random.randint(0,len(p))]
    dx=np.absolute(U[0]-G[0])/2
    dy=np.absolute(U[1]-G[1])/2
    if U[0]>G[0]:
        nx=G[0]+dx
    else:
        nx=G[0]-dx
    if U[1]>G[1]:
        ny=G[1]+dy
    else:
        ny=G[1]-dy
    p.append([nx,ny])


plt.figure(figsize = (10,10))
plt.xlim(0,1)
plt.ylim(0,1)
for i in range(len(p)):
    plt.plot(p[i][0], p[i][1], 'bo')

Результатом выполнения данного кода является данное изображение.

Необходимый результат был достигнут! Однако стало интересно, что будет, если построить фрактал аналогичного типа для квадрата. Код особо ничем не отличается. Задаются четыре вершины, и в этот раз в качестве первой точки фрактала была взята случайная вершина квадрата.

Хм, получается какое-то случайное распределение. Было решено изменить расположение новой точки на отрезке между существующей точкой и вершиной. Для этого в коде требуется изменить делитель при нахождении dx и dy. Логично предположить, что, если взять делитель равный 1, то все точки соберутся в вершинах квадрата, так что изменяя делитель от 1 до 2, я остановился на наиболее приглянувшемся изображении при делителе равном 1,7.

A=[0,0]
B=[0,1]
C=[1,1]
D=[1,0]
ABCD=[A,B,C,D]

p=[[A[0],A[1]]]

for _ in range(10000):
    U=ABCD[np.random.randint(0,len(ABCD))]
    G=p[np.random.randint(0,len(p))]
    dx=np.absolute(U[0]-G[0])/1.7
    dy=np.absolute(U[1]-G[1])/1.7
    if U[0]>G[0]:
        nx=G[0]+dx
    else:
        nx=G[0]-dx
    if U[1]>G[1]:
        ny=G[1]+dy
    else:
        ny=G[1]-dy
    p.append([nx,ny])


plt.figure(figsize = (10,10))
for i in range(len(p)):
    plt.plot(p[i][0], p[i][1], 'ro')

Результат интересный, так что логичным продолжением было построение фрактала для пятиугольника. Код приводить не вижу смыла, он аналогичен предыдущему.

Выглядит все это красиво, но для фигур с большим количеством углом писать отдельный код уже не хотелось, и было решено создать функцию, которая строила бы фракталы подобного типа для любого количества углов (больше 3). Привожу код данной функции.

def frac(N=3, X=10, Y=10, R=5, n=10000, d=2):
    
    s=[]
    for i in range(0,N):
        a = 2*np.pi*i/N
        s.append([R*np.cos(a)+X, R*np.sin(a)+Y])
        
    p=[[s[0][0],s[0][1]]]
    
    for _ in range(n-1):
        U=s[np.random.randint(0,len(s))]
        G=p[np.random.randint(0,len(p))]
        dx=np.absolute(U[0]-G[0])/d
        dy=np.absolute(U[1]-G[1])/d
        if U[0]>G[0]:
            nx=G[0]+dx
        else:
            nx=G[0]-dx
        if U[1]>G[1]:
            ny=G[1]+dy
        else:
            ny=G[1]-dy
        p.append([nx,ny])

    plt.figure(figsize = (10,10))
    for i in range(len(p)):
        plt.plot(p[i][0], p[i][1], 'ko', ms =3)
    for i in range(len(s)):
        plt.plot(s[i][0], s[i][1], 'ro', ms =6)

Изначально количество углов (N) задано 3, X (10) и Y (10) задают центр окружности, R (5) - радиус окружности, необходимой для нахождения вершин фигуры, число точек (n) изначально равно 10000, а делитель (d) равен 2. Все эти данные являются необязательными параметрами функции. Вот примеры работы функции для некоторых фигур.

На этом на данный момент мои эксперименты заканчиваются. Надеюсь было интересно, а может быть кому-то даже полезно. Спасибо за внимание!