Matching an image with a fractal, advice needed

advicefractalsnumerical methods

Introduction

I'm trying to find a fractal function that spells out a name on the complex plane. I'm doing this by numerically 'training' a program to match the desired image / name.

I don't expect any answers regarding the coding, but any help is appreciated and I'm very open to critique. It is written in C++, do keep in mind that it's just a working prototype and that I'm pretty new to C++. (the github link is in the last section)

Objective

Given a region $E$ in the complex plane. Find a complex fractal function $f:\mathbb{C}\to\mathbb{C}$ that when the escape time algorithm is applied it is convergent in all points of $E$ and divergent outside of $E$. I'm not concerned here with making the fractal look good as I have plenty of programs that can draw fractals and make them look appealing, I'd just like to be able to find a complex function $f$ that converges in $E$ to spell out the name. Generally it also doesn't matter where this region $E$ is or how much area it covers, as long as it spells out the name.

I'm not hoping that the fractal perfectly matches the given region $E$, but rather that it matches it to a point where you can easily read it and it of course looks appealing. Ideally I'm hoping for it to match the letters but have mandelbrot like shapes on the edges.

My Method

fractal function

I guess a form of the function $f(z;\boldsymbol\theta)$, here $\boldsymbol\theta$ denotes the parameters this fractal function takes. So far in practice I have only used function series, here are a few examples.
$$
\sum_{n=1}^N\theta_nz^n+c\\
\sum_{n=1}^N\exp{(in\theta_nz)}+c
$$

Where $c\in\mathbb C$ is the coordinate (just like in the mandelbrot $f(z)=z^2+c$), which I leave out of the notation of $f$ for simplicity.

In these examples there are $N$ parameters, but I have used other forms so the amount of parameters may vary. Expected is that if the amount of parameters goes up, the fractal can take more complicated shapes. There is of course a balance between execution time and getting the desired degree of complexity.

I have also used complex parameters, but after lots of trial and error I just started to use two real parameters for each term inside the sum as these forms of $f$ gave very diverse fractals which made me think these were more suited to match the region $E$. The form of $f$ I have been using predominantly is the following:
$$
\begin{align}
\tag{1}
\sum_{n=1}^N \theta_{n,1}(\exp(i\theta_{n,2}z)+c)
\end{align}
$$

Where $\theta_{n,1}$ and $\theta_{n,2}$ are real parameters. This means that for this example there are $2N$ parameters.

To give some insight into the coding process and to show I'm not making any mistakes, equation (1) can be translated into code like the following:

for (int n=1; n<N; n++)
{
    double A = exp(-theta_2[n] * y);
    double a = theta_2[n] * x;
    //temporary variables so x and y don't change during the summation
    tx += theta_1[n] * (A * cos(a) + c); 
    ty += theta_1[n] * (A * sin(a) + c);    
}
x = tx;
y = ty;

cost function

I then define a cost function $C_f(\boldsymbol{\theta})$ which takes all the parameters as input and outputs a number between $0$ and $1$ using the fractal function $f$. If this cost is $0$ the current fractal function $f(z;\boldsymbol{\theta})$, with its current configuration of parameters, will converge exactly in every point of $E$ and diverge in every point outside of it, hence the ideal solution. Consequently if the cost is $1$ it implies that it diverges inside of every point of $E$ and converges outside of $E$. As discussed earlier I'm not looking for a fractal function and parameter configuration , $f(z;\boldsymbol\theta)$, that perfectly matches the region $E$, so the cost should generally be very low but not $0$.

In practice I have found that a lot of the parameter configurations result in completely divergent fractals. Hence I have added a weight to the cost function so that it prefers convergent pixels over divergent ones. This also means that the cost can never be $1$, but this doesn't matter as we don't care about the least optimal solution.

Let $v_c\in[0,1]$ be the convergence value of $f(z;\boldsymbol\theta)$ in the point $c$, being $0$ if the fractal converges and $1$ if it diverges. The value $v_c$ is obtained using a simple escape time algorithm. The cost function $C_f$ is then defined as:
$$
C_f(\boldsymbol\theta)=
\begin{cases}
\tag 2
\sqrt{\dfrac{1}{w}\sum\limits_c(v_c-d_c)^2}:v_c < d_c\\
\sqrt{\sum\limits_c(v_c-d_c)^2}:v_ \ge d_c
\end{cases}
$$

Where the sum is over all the coordinates of the relevant pixels, for now we can assume it is just a simple rectangle in the complex plane where we want the name to be spelled out.

$d_c$ is the desired value of the fractal in the point $c$. So if the name, spelled out on the complex plane, crosses over the point $c$ then $d=0$. Keep in mind that the value $v_c$ can be any number between $0$ and $1$. So the desired pixels $d_c$ are also chosen to be between $0$ and $1$, which just corresponds to the brightness (black=fully convergent, white=fully divergent) of the image spelling out the name.

The cost function is also divided by a weight $w$ for values lower than the desired $d_c$. This is to put a preference on convergent pixels as discussed earlier.

We can also clearly see that the cost function returns $0$ if all iterated values equal the desired values.

descending the gradient of the cost function

So to achieve my goal of finding a fractal that converges inside $E$, I try to minimize the cost function by changing the parameters $\boldsymbol\theta$.

I do this by numerically calculating the gradient $\nabla_{\boldsymbol\theta}C_f$ and taking a small step in the opposite direction. The elements of this gradient are nothing more than:
$$
\tag 3
\dfrac{\partial C_f(\boldsymbol\theta)}{\partial\theta_n}=\dfrac{C_f(\theta_1,..,\theta_n+\Delta\theta_n,..)-C_f(\theta_1,..,\theta_n,..)}{\Delta\theta_n}
$$

I then take a small step along $-\nabla_{\boldsymbol\theta}C_f$, this just means to subtract $s\frac{\partial C_f(\boldsymbol\theta)}{\partial \theta_n}$ from the parameters $\theta_n$. Where $s$ is the step size, I have tried all kinds of step sizes and I'm still not really sure if there is a most efficient way of choosing one.

I then calculate the cost again with those new parameters and compare it to the previous cost. Here two things can generally happen:

  • Any small step taken in the direction of the negative gradient does not go downhill, or in other words the new cost is higher than the previous cost.
  • There is a step that does go downhill, so the cost decreases after taking the step.

If the latter happens, nothing needs to be done besides repeat the process until the cost is essentially $0$ and the fractal matches the given region $E$.

Now it almost always does not go downhill. I'm not really sure why this happens, but my suspicion is that the parameter space is very chaotic and possesses a lot of local minima. It also seems to have a harder time going downhill if there are more parameters, which I can't really explain. I have also observed that the partial derivatives seem to be very inconsistent. They usually grow very big for small $\Delta \theta_n$ and generally decrease as $\Delta\theta_n$ increases ($\Delta\theta_n\in[10^{-6},1]$). Which leads me to believe my method of calculating the partial derivatives isn't very accurate at all.

There are two things I can do here if the step taken does not go downhill.

  • Choose a new random point in the parameter space, by this I mean picking a new set of random parameters $\boldsymbol\theta$ and starting the process over again.
  • Taking a small step in the negative gradient anyway and continuing the process.

Neither of these really get me to the desired solution of finding the correct parameters $\boldsymbol\theta$ so that $f(z;\boldsymbol\theta)$ matches $E$.

The Question

The main question really is:

Is my current method a good approach? Is there a better way of doing this? Or is there anything that makes the current method flawed?

Note that I'm not expecting any answers on the programming point of view, but any input is greatly appreciated!

Continuing with the current method. I have programmed a working prototype in C++, but I lack the mathematical knowledge to find a good form of the function $f$, so:

What is a good form of the function $f(z;\boldsymbol\theta)?$

I'm continuously changing the parameters $\boldsymbol\theta$ to lower the cost, but I'm not even sure that a parameter configuration which has a cost very close to zero exists for the current form of $f$.

I wasn't even sure if I should post this question on here, it feels kind of vague what I'm asking but I've really hit a dead end. I don't know where to go from here and would really appreciate any ideas, tips or feedback anyone has to offer. Also my apologies if this post is too long, I wanted to make sure everything was super clear to avoid any confusion.

The Code

Again I don't expect any answers regarding the coding, but if you'd like to check it out please do so on https://github.com/catmousedog/FractalDrawer.git.

The code does do one thing I haven't fully explained in the section cost function, but that isn't all too necessary. It doesn't iterate the fractal function $f$ or the cost function $C_f$ over a simple rectangle in the complex plane, but rather over a set of predetermined relevant pixels. It's easiest to understand this by looking at the following picture.

This picture represents the desired set of pixels $d$. The fractal function $f$ should ideally converge if the pixel is black ($d=0$) and diverge if it is white ($d=1$), these two make up for all the relevant pixels. Hence the red pixels are irrelevant, they are ignored in all calculations.

Adding the concept of relevant pixels was naively done by me thinking it would increase the chance of finding a function which fits the desired region $E$ as there are no constraints on the function outside of the relevant pixels. Whether this is a good or bad idea, I don't really know.

Best Answer

Shapes can be approximated by Julia sets of rational functions with a large number of roots [1]:

we are [...] interested in factored rationals of the form

$$ R(q) = e^C \frac{ \prod (q - t_i)^{T_i} }{ \prod (q - b_i)^{B_i} } $$

Property 1: If an iterate $q$ is sufficiently close to a top root $t_i$, then it is inside the filled Julia set $J$.

Property 2: If an iterate $q$ is sufficiently close to a bottom root $b_i$, then it is outside of $J$.

Property 3: The $T_i$ term expands and contracts the attractive region of $t_i$. The $B_i$ does the same for repulsion with $b_i$.

The paper uses quaternions for 3D shapes, but presumably the same techniques can apply in 2D with complex numbers, given that the paper references [2] which is about doing something similar for Jordan curves with polynomials, and combinations of Jordan curves with rational functions:

figure 3

Both papers are complementary, [2] has lots of proofs while [1] has practical tips too.

[1]: "Quaternion Julia Set Shape Optimization", Theodore Kim, Eurographics Symposium on Geometry Processing, Volume 34 (2015), Number 5, https://doi.org/10.1111/cgf.12705 online https://tkim.graphics/JULIA/

[2]: "Shapes of Polynomial Julia Sets", Kathryn A. Lindsey, Ergodic Theory and Dynamical Systems, Volume 35, Issue 6, September 2015, pp. 1913 - 1924, DOI: https://doi.org/10.1017/etds.2014.8 preprint https://arxiv.org/abs/1209.0143