📜 ⬆️ ⬇️

Hedgehog in fractal fog

This article is the last of my fractal series. In Habrastier "Drawing pictures with the Hilbert Curve," it was told about a kitten named Gav, in Habrastier "Kosh on a complex plane" - about the flow of fractals into the horizon, in Habrastier " Fractals Night" about the escape time algorithm. In this article we will talk about a hedgehog in the fog and, of course, about the cat.




')
Anyone who has ever experienced a bit of fractals has seen beautiful pictures of Julia sets, which are defined by a square polynomial; but it is interesting to solve the inverse problem: let there be some picture, and we go over it to come up with such a polynomial, which will give some approximation of the original picture. A picture with a hedgehog above demonstrates this idea.

So, consider the cat K.



We need to come up with a polynomial f such that inside a cat the sequence f ( z ), f ( f ( z )), f ( f ( f ( z ))), ... is bounded, and outside the cat it tends to infinity. For a long time I did not think about this question, but immediately looked on the Internet and found a wonderful article by Kathryn A. Lindsey , “Shapes of polynomial Julia sets” , 2013, arXiv: 1209.0143 . This article proves that for a “good” cat and for a given accuracy δ, such a polynomial can be invented, and the proof is constructive.

Consider this design. Let φ 'be a conformal mapping that the exterior of the unit circle translates into the exterior of a cat. Let φ ( z ) = φ '((1 + ε) z ) be the original map corrected to some small ε.



Further, let ω ( z ) = c - n ( z - r 0 ) ( z - r 1 ) ⋅… ⋅ ( z - r n −1 ), where c is the coefficient at z in the expansion of the map φ in the Laurent series, and r k = φ (e ik / n ) are images of roots of unit nth degree. It is proved that the polynomial f ( z ) = z (ω ( z ) + 1) is the desired one for sufficiently large n .

Therefore, to solve our problem, it is necessary to learn how to generate the corresponding conformal mapping. As a result of a small search on the Internet, the zipper package was found. This package allows you to find a conformal mapping of the interior of the unit circle to the area bounded by the polyline. Although this package was written twenty years ago in an ancient dialect , it was not difficult to assemble it and use it.

A piece of ancient dialect
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 

Coat the zipper with python glue
 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 

To use this package, you need to build a broken line in the image. I didn’t do some kind of automatic solution, but loaded the image into Inkscape editor and circled it, and parsing the SVG format is simple. This is convenient in that it allows you to experiment with a polyline.

Parsim paths in SVG file
 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 



Next, we need a mapping from the exterior to the exterior, and the package finds the mapping from the inside to the inside. Here you just need to pair the generated map with the inverse. That is, φ '( z ) = 1 / (1 / z ), where ψ is the map that generates the packet. And it is necessary to feed the already inverted broken line to the input of the packet.

We mate
 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) 



In the definition of the polynomial f , the coefficient c is also involved. For an approximate calculation of its value, we use the following technique. Let be

f ( z ) = c z + a + a 1 / z + a 2 / z 2 + a 3 / z 3 + ...

Consider some finite, but rather large, part of this series. Substitute in a series of values ​​of the roots of units of n -th degree, where n is greater than the number of members of this piece.

f (e ik / n ) = c e ik / n + a + a 1 e −2π ik / n + a 2 e −4π ik / n + a 3 e − ik / n + ..., k = 0, ..., n - 1.

Now multiply each line by e −2π ik / n and add all the lines. Since the sum of all roots of a unit is 0, only n c will remain on the right. Therefore, we can put

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

Putting it all together and check what happened. The figure below shows the “thermal” graph of the logarithm of the module f ( z ) (the redder, the greater the value). As you can see, inside the cat the values ​​of the polynomial are small, and outside the cat the polynomial increases, so it should be. Notice also how the values ​​of f (e ik / n ) are distributed (green dots), because of this effect I had to draw a cat whose legs are far enough apart from each other.



Now we simply use some kind of algorithm to draw the Julia set and get a cat whose boundary is fractal.

For example, the escape time algorithm
 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) 



If our image consists of several areas A 1 , A 2 , ..., A q , then in the above article it is proposed to use instead of a polynomial a rational map f ( z ) defined as follows. For each region Ar , we write the product ω r ( z ) as above. Then the desired rational function f ( z ) is defined by the formula f ( z ) = z / (1 / ω 1 ( z ) +… + 1 / ω q ( z )).

For example, suppose we have a hedgehog in the fog E.



We will circle the hedgehog and the hill (the lower border of the hill is outside the picture).



Hedgehog in a vacuum.



A hill without a hedgehog.



Together.



Enlarged belly of a hedgehog with fleas, hedgehogs.



As always, the source code can be found on a githaba .

I repeat the link to the article with the help of which everything turned out: Kathryn A. Lindsey , “Shapes of polynomial Julia sets” , 2013, arXiv: 1209.0143 . I note, the article is easy to read, so you can pick up the evidence for a cup of tea.

I hope it was fun.

Failed doubles





Source: https://habr.com/ru/post/230009/


All Articles