Їжачок у фрактальному тумані

Ця стаття — остання з серії моїх хабрастатей про фрактали. У хабрастатье «Малюємо за допомогою кривої Гільберта» розповідалося про кошеня по імені Гав, в хабрастатье «Кіш на комплексній площині» — про перетікання фракталами в горизонт, в хабрастатье «Ніч фракталів» — про алгоритм часу утікання. У цій статті піде мова про їжачка в тумані і, звичайно ж, про кота.





Кожен, хто хоч трохи стикався з фракталами, бачив красиві картинки множин Жюлиа, які визначаються квадратним многочленом; але цікаво вирішувати зворотну задачу: нехай є деяка картинка, а ми ходимо по ній придумати такий многочлен, який дасть деяке наближення вихідної картинки. Картинка з їжачком вище демонструє цю ідею.

Отже, розглянемо кота K.



Нам потрібно придумати такий многочлен f, щоб усередині кота послідовність f(z), f(f(z)), f(f(f(z))),… була обмеженою, а поза кота — прагнула до нескінченності. Довго над цим питанням я не думав, а відразу ж пошукав в інтернеті і знайшов чудову статтю Kathryn A. Lindsey, «Shapes of поліноміальні Julia sets», 2013, arXiv:1209.0143. У цій статті доводиться, що для «хорошого» кота і для заданої точності δ такий многочлен придумати можна, причому доказ конструктивно.

Розглянемо цю конструкцію. Нехай φ' — конформное відображення, яке зовнішність одиничної окружності переводить в зовнішність кота. Нехай φ(z) = φ'((1 + ε)z) — подправленное на деякий малий ε початкове відображення.



Далі, нехай ω(z) = cn(z r0)(z r1) ⋅… ⋅ (z rn-1), де c — коефіцієнт при z в розкладанні відображення φ в ряд Лорана, а rk = φ(eik/n) — образи коренів з одиниці n-ой ступеня. Доводиться, що многочлен f(z) = z(ω(z) + 1) буде шуканим при досить великому n.

Отже, для рішення нашої задачі необхідно навчитися генерувати відповідне конформное відображення. В результаті невеликого пошуку в інтернеті, був знайдений пакет zipper. Цей пакет дозволяє знаходити конформное відображення нутрощі одиничного кола на область, обмежену ламаної. Хоча цей пакет написав двадцять років тому на стародавньому говіркою, зібрати його і скористатися ним не склало праці.

Шматочок на стародавньому говіркою
call invers
write(4,*)z1,z2,z3,zrot1,zto0,zto1,angler,zrot2
do 981 j=4,n-2,2
981 write(4,999)a(j),b(j),c(j)
zm=dcmplx(0.d0,0.d0)
ierr=0
do 982 j=1,n
if(cdabs(zm-z(j)).lt.1.d-16)then
ierr=1
jm=j-1
write(*,*)' WARNING: prevertices',j,' and',jm,' are equal'
endif
zm=z(j)
x=dreal(z(j))
y=dimag(z(j))
982 write(3,999)x,y
if(ierr.eq.1)then

Обмазуємо zipper клеєм на пітоні
def calc_phi(points):
with open("init.dat", "w") as output_file:
for point in points:
output_file.write(str(point.real) + " " + str(point.imag) + "\n")
output_file.write("\n0.0 0.0\n")
os.system("echo \"init.dat\n200\npoly.dat\" | ./polygon")
os.system("./zipper")
os.system("rm init.dat")

def transform_points(points):
with open("fftpts.dat", "w") as output_file:
for point in points:
output_file.write(str(point.real) + " " + str(point.imag) + "\n")
os.system("echo \"0\nfftpts.dat\nfftpts.img\n0\" | ./forward")
transformed_points = []
with open("fftpts.img", "r") as input_file:
for line in input_file.readlines():
x, y = float(line[0:25]), float(line[25:])
transformed_points.append(complex(x, y))
os.system("rm fftpts.dat")
os.system("rm fftpts.img")
return transformed_points

Щоб використовувати цей пакет, потрібно по зображенню побудувати ламану. Я не став робити якогось автоматичного рішення, а завантажив зображення в редактор Inkscape і обвів його, а розпарсити SVG формат простіше простого. Це зручно тим, що дозволяє експериментувати з ламаною.

Парсим шляху в SVG-файл
def read_points_from_svg(file_name, path_n):
with open(file_name, "r") as input_file:
content = input_file.read()
soup = BeautifulSoup.BeautifulSoup(content)
path = soup.findAll("path")[path_n]
data = path.get("d").split(" ")
x, y = 0, 0
is_move_to = False
is_relative = False
points = []
for d in data:
if d == "m":
is_move_to = True
is_relative = True
elif d == "M":
is_move_to = True
is_relative = False
elif d == "z":
pass
elif d == "Z":
pass
elif d == "l":
is_move_to = False
is_relative = True
elif d == "L":
is_move_to = False
is_relative = False
else:
dx, dy = d.split(",")
dx = float(dx)
dy = float(dy)
if is_move_to:
x = dx
y = dy
is_move_to = False
else:
if is_relative:
x += dx
y += dy
else:
x = dx
y = dy
points.append(complex(x, y))
return points



Далі, нам потрібно відображення зовнішності в зовнішність, а пакет знаходить відображення нутрощі в середину. Тут потрібно просто сполучити генерується відображення інверсією. Тобто φ'(z) = 1 / ψ(1/z), де ψ — відображення, яке генерує пакет. А на вхід пакета треба подавати вже обернену ламану.

Сопрягаем
def invert_points_1(points, epsilon):
inverted_points = []
for point in points:
inverted_points.append(1 / ((1 + epsilon)*point))
return inverted_points

def invert_points_2(points, shift):
inverted_points = []
for point in points:
inverted_points.append(1 / point + shift)
return inverted_points


if __name__ == "__main__":
...
inverted_points = invert_points_1(points, epsilon)
transformed_points = transform_points(inverted_points)
inverted_transformed_points = invert_points_2(transformed_points, shift)



У визначенні многочлена f ще бере участь коефіцієнт c. Для наближеного обчислення його значення скористаємося наступним прийомом. Нехай

f(z) = cz + a + a1/z + a2/z2 + a3/z3 +…

Розглянемо деяку кінцеву, але досить велику, частину цього ряду. Підставимо в ряд значення коренів з одиниці n-ой ступеня, де n більше числа членів цього шматка.

f(eik/n) = ceik/n + a + a1e-2πik/n + a2e-4πik/n + a3e-6πik/n + ..., k = 0, ..., n — 1.

Тепер помножимо кожний рядок на e-2πik/n і складемо всі рядки. Так як сума всіх коренів з одиниці дорівнює 0, то справа залишиться тільки nc. Тому можна покласти

c = (f(1) ⋅ 1 + f(ei/n)e-2πi/n +… + f(ei(n — 1)/n)e-2πi(n — 1)/n)) / n.

Зберемо всі разом і перевіримо, що вийшло. На малюнку нижче наведений «теплової» графік логарифма модуля f(z) (ніж красно, тим більше значення). Як видно, всередині кота значення многочлена маленькі, а поза кота многочлен зростає, так і повинно бути. Зауважте також, як розподіляються значення f(eik/n) (зелені крапки), з-за цього ефекту довелося малювати кота, у якого ноги знаходяться досить далеко один від одного.



Тепер просто-напросто застосовуємо який-небудь алгоритм для малювання безлічі Жюлиа і отримуємо кота, межа якого фрактальна.

Наприклад, алгоритм часу утікання
def get_value(points, z, radius):
result = 1
for point in points:
result *= (z - point) / radius
return z * (1 + result)


def get_radius(points):
n = len(points)
result = 0
for k in range(n):
result += points[k] * cmath.exp(-2j*math.pi*k/n)
return result / n


def draw_fractal(image, points, scale, shift, bound, max_num_iter, radius):
width, height = image.size 
draw = ImageDraw.Draw(image)
for y in range(0, height):
for x in range(0, width):
z = (complex(x, y) - shift) * scale
n = 0
while abs(z) < bound and n < max_num_iter:
z = get_value(points, z, radius)
n += 1
if n < max_num_iter:
color = ((100 * n) % 255, 128 + (50 * n) % 255, 128 + (75 * n) % 255)
draw.point((x, y, color)



Якщо ж у нас зображення складається з декількох областей A1, A2, ..., Aq, то у вищезгаданій статті пропонується використовувати замість многочлена раціональне відображення f(z), визначене наступним чином. Для кожної області Ar, запишемо твір ωr(z) як зроблено вище. Тоді шукана раціональна функція f(z) визначається формулою f(z) = z / (1/ ω1(z) +… + 1 / ωq(z)).

Наприклад, нехай у нас є їжачок в тумані Е.



Обведемо їжачка і пагорб (нижня межа пагорба знаходиться за межами картинки).



Їжачок у вакуумі.



Пагорб без їжачка.



Всі разом.



Збільшене черевце їжачка з блохами-їжаками.



Як завжди, вихідний код можна знайти на гітхабі.

Ще раз повторю посилання на статтю, за допомогою якої все вийшло: Kathryn A. Lindsey, «Shapes of поліноміальні Julia sets», 2013, arXiv:1209.0143. Зазначу, стаття легко читається, так що можете поразбирать докази за чашкою чаю.

Сподіваюся було весело.

Невдалі дублі






Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.