Improved Cuckoo Search with Luus-Jakoola Heuristics for the IFS Inverse Problem of Binary Self-Similar Fractal Images

. This paper addresses the following problem: how to reconstruct a given binary self-similar fractal image through iterated functions systems. This means to obtain an iterated function system (IFS) whose attractor is a good approximation of the input image. This problem is known to be a very diﬃcult multivariate nonlinear continuous optimization problem. To tackle this issue, this paper introduces a new hybrid method comprised of a modiﬁcation of the original cuckoo search method for global optimization called improved cuckoo search (ICS) along with the Luus-Jakoola heuristics for local search. This hybrid methodology is applied to three fractal examples with 3, 4, and 26 contractive functions. Our experimental results show that the method performs very well and provides visually satisfactory solutions for the instances in our benchmark. The numerical values of the similarity index used in this work also show that the results are not optimal yet, suggesting that the method might arguably be further improved.


Introduction
Fractals are one of the most challenging and intriguing mathematical shapes ever defined.Basically, they are geometric figures created by repeating a simple process over and over so that it yields a self-similar pattern across different scales.Interestingly, in the case of fractals, the scale factor of this replicating pattern is not an integer number, but a real one.Such a number is called the fractal dimension and it is usually larger than the topological dimension of the fractal [3,9].Fractals have become ubiquitous objects in popular culture, particularly since the 80s of last century, owing to the technological advances in hardware and software and the widespread availability of personal computers.They are also very popular in science due to their ability to describe many growing patterns and natural structures commonly found in real-life objects: branches of trees, river networks, coastlines, mountain ranges, and so on.Furthermore, fractals have found remarkable applications in computer graphics, scientific visualization, image processing, dynamical systems, telecommunications, medicine, biology, arts, and many other fields [1,3,9,11,12].
There are several methods described in the literature to obtain fractal images.They include the Brownian motion, escape-time fractals, finite subdivision rules, L-systems, strange attractors of dynamical systems, and many others [1,3].One of the most popular methods is the Iterated Function Systems (IFS), originally conceived by J.E. Hutchinson [13] and popularized by M. Barnsley in [1].Roughly, an IFS consists of a finite system of contractive maps on a complete metric space.It can be proved that the Hutchinson operator over the set of all compact subsets of this space has a unique non-empty compact fixed set for the induced Hausdorff metric, called the attractor of the IFS.The graphical representation of this attractor is (at least approximately) a self-similar fractal image.Conversely, each self-similar fractal image can be represented by an IFS.Obtaining the parameters of such IFS is called the IFS inverse problem.Basically, it consists of solving an image reconstruction problem: given a fractal image, compute the IFS whose attractor approximates the input image accurately.
This IFS inverse problem has shown to be extremely difficult.In fact, the general case is still unsolved and only partial solutions have been reached so far.In this paper we propose a new approach to address this problem for the case of binary fractal images.Our methodology consists of the hybridization of a bio-inspired metaheuristics based on the cuckoo search algorithm for global optimization and a local search procedure.In particular, we consider a modification of the original cuckoo search method called the improved cuckoo search (ICS), which is based on the idea of allowing the method parameters to change over the generations [16].This method is hybridized with the Luus-Jakoola (LJ) heuristics, a search method aimed at improving the local search step to refine the quality of the solution.
The structure of this paper is as follows: Section 2 introduces the basic concepts and definitions about the iterated function systems and the IFS inverse problem.Then, Section 3 describes the original and the improved cuckoo search algorithms.Our proposed method is described in detail in Section 4, while the experimental results are briefly discussed in Section 5.The paper closes with the main conclusions and some ideas about future work in the field.

Iterated Function Systems
An Iterated Function System (IFS) is a finite set Øφ i Ù i 1,...,η of contractive maps φ i : Ω Ω defined on a complete metric space M ÔΩ, ΨÕ, where Ω R n and Ψ is a distance on Ω.We refer to the IFS as W ØΩ; φ i , . . ., φ η Ù.For visualization purposes, in this paper we consider that the metric space ÔΩ, ΨÕ is R 2 along with the Euclidean distance d 2 , which is a complete metric space.Note, however, that our method can be applied to any other complete metric space of any dimension without further modifications.In this two-dimensional case, the affine transformations φ κ are of the form: 1.In fact, µ κ detÔΘ κ Õ 1 meaning that φ κ shrinks distances between points.Let us now define a transformation called the Hutchinson operator, Υ , on the compact subsets of Ω, HÔΩÕ, by: with B È HÔΩÕ.If all the φ κ are contractions, Υ is also a contraction in HÔΩÕ with the induced Hausdorff metric [1,13].Then, according to the fixed point theorem, Υ has a unique fixed point, Υ ÔAÕ A, called the attractor of the IFS.
Consider a set of probabilities P Øω 1 , . . ., ω η Ù, with η κ 1 ω κ 1.There exists an efficient method, known as probabilistic algorithm, for the generation of the attractor of an IFS.It follows from the result ØΞ j Ù j A provided that Ξ 0 È Ω, where: Ξ j φ κ ÔΞ j¡1 Õ with probability ω κ 0, see [2].Picking an initial point Ξ 0 , one of the mappings in the set Øφ i , . . ., φ η Ù is chosen at random using the weights Øω 1 , . . ., ω η Ù.The selected map is then applied to generate a new point, and the same process is repeated again with the new point and so on.As a result of this stochastic iterative process, we obtain a sequence of points that converges to the fractal as the number of points increases [1,10].

The Collage Theorem
The Collage Theorem basically says that any digital image I can be approximated through an IFS [1,10].In particular, it states that given a non-empty compact subset I È HÔΩÕ, the Hausdorff metric HÔ., .Õ, a non-negative real threshold value ǫ 0, and an IFS W ØΩ; φ i , . . ., φ η Ù on Ω with contractivity factor 0 s 1 (the maximum of the contractivity factors

The IFS Inverse Problem
Suppose that we are given an initial fractal image F .The Collage Theorem says that it is possible to obtain an IFS W whose attractor has a graphical representation F that approximates F accurately according to a similarity function S, which measures the graphical distance between F and F [1].Note that once W is computed, F Υ ÔI Õ for any (not necessarily fractal) initial image I .Mathematically, this means that we have to solve the optimization problem: minimize The problem ( 3) is a continuous constrained optimization problem, because all free variables in ØΘ κ , Σ κ , ω κ Ù i are real-valued and must satisfy the condition that the corresponding functions φ κ have to be contractive.It is also a multimodal problem, since there can be several global or local minima of the similarity function.The problem is so difficult that only partial solutions have been reported so far, but the general problem still remains unsolved to a large extent.In this paper we address this problem by using a hybrid approach based on the cuckoo search algorithm described in next paragraphs.
3 The Cuckoo Search Algorithms

Original Cuckoo Search (CS)
The Cuckoo search (CS) is a powerful metaheuristic algorithm originally proposed by Yang and Deb in 2009 [18].Since then, it has been successfully applied to difficult optimization problems [5,14,17,19].The algorithm is inspired by the obligate interspecific brood-parasitism of some cuckoo species that lay their eggs in the nests of host birds of other species to escape from the parental investment in raising their offspring and minimize the risk of egg loss to other species, as the cuckoos can distributed their eggs amongst a number of different nests.This interesting and surprising breeding behavioral pattern is the metaphor of the cuckoo search metaheuristic approach for solving optimization problems.In the cuckoo search algorithm, the eggs in the nest are interpreted as a pool of candidate solutions of an optimization problem while the cuckoo egg represents a new coming solution.The ultimate goal of the method is to use these new (and potentially better) solutions associated with the parasitic cuckoo eggs to replace the current solution associated with the eggs in the nest.This replacement, carried out iteratively, will eventually lead to a very good solution of the problem.
In addition to this representation scheme, the CS algorithm is also based on three idealized rules [18,19]: 1.Each cuckoo lays one egg at a time, and dumps it in a randomly chosen nest; 2. The best nests with high quality of eggs (solutions) will be carried over to the next generations; 3. The number of available host nests is fixed, and a host can discover an alien egg with a probability p a È Ö0, 1×.In this case, the host bird can either throw the egg away or abandon the nest so as to build a completely new nest in a new location.For simplicity, this assumption can be approximated by a fraction p a of the n nests being replaced by new nests (with new random solutions at new locations).The basic steps of the CS algorithm are summarized in the pseudocode shown in Table 1.Basically, the CS algorithm starts with an initial population of n host nests and it is performed iteratively.The initial values of the jth component of the ith nest are determined by the expression x j i Ô0Õ rand.Ôup j i ¡low j i Õ low j i , where up j i and low j i represent the upper and lower bounds of that jth component, respectively, and rand represents a standard uniform random number on the interval Ô0, 1Õ.With this choice, the initial values are within the search space domain.These boundary conditions are also controlled in each iteration step.
For each iteration t, a cuckoo egg i is selected randomly and new solutions x t 1 i are generated by using the Lévy flight.According to the original creators of the method, the strategy of using Lévy flights is preferred over other simple random walks because it leads to better overall performance of the CS.The general equation for the Lévy flight is given by: where t indicates the number of the current generation, and α 0 indicates the step size, which should be related to the scale of the particular problem under study.The symbol is used in Eq. ( 4) to indicate the entry-wise multiplication.Note that Eq. ( 4) is essentially a Markov chain, since next location at generation t 1 only depends on the current location at generation t and a transition probability, given by the first and second terms of Eq. ( 4), respectively.This transition probability is modulated by the Lévy distribution as: which has an infinite variance with an infinite mean.From the computational standpoint, the generation of random numbers with Lévy flights is comprised of two steps: firstly, a random direction according to a uniform distribution is chosen; then, the generation of steps following the chosen Lévy distribution is carried out.The authors suggested to use the Mantegna's algorithm for symmetric distributions (see [19] for details), which computes the factor: where Γ denotes the Gamma function and β 3 2 in the original implementation by Yang and Deb [19].This factor is used in Mantegna's algorithm to compute the step length ς as: , where u and v follow the normal distribution of zero mean and deviation σ 2 u and σ 2 v , respectively, where σ u obeys the Lévy distribution given by Eq. ( 6) and σ v 1.Then, the stepsize ζ is computed as ζ 0.01 ς Ôx ¡ x best Õ.Finally, x is modified as: x x ζ.∆ where ∆ is a random vector of the dimension of the solution x and that follows the normal distribution N Ô0, 1Õ.The CS method then evaluates the fitness of the new solution and compares it with the current one.In case the new solution brings better fitness, it replaces the current one.On the other hand, a fraction of the worse nests (according to the fitness) are abandoned and replaced by new solutions so as to increase the exploration of the search space looking for more promising solutions.The rate of replacement is given by the probability p a , a parameter of the model that has to be tuned for better performance.Moreover, for each iteration step, all current solutions are ranked according to their fitness and the best solution reached so far is stored as the vector x best .

Improved Cuckoo Search (ICS)
The improved cuckoo search (ICS) method was proposed in [16] to enhance the performance of the original CS.It is based on the idea of allowing its parameters p a and α to change over the generations, as opposed to their fixed values in the original CS.In ICS the parameter p a is modified as: where the subscripts M and m are used to indicate the maximum and minimum values of the parameter respectively, and Λ indicates the total number of iterations.According to Eq. ( 7), the parameter p a is now decreased linearly with the number of iterations from a maximum value p aM until a minimum one, p am .At early iterations, its value is high to enforce the diversity of solutions in the algorithm.This diversity is decreasing over the time so as to intensify the search using the best candidates of the population in final iterations for a better finetuning of the solutions.The parameter α, also assumed constant in the CS, is primarily used to promote exploration of the search space.Therefore, it makes sense to modify it dynamically starting from a high value, α M , to perform a extensive exploration and gradually reducing it until a low value, α m , to promote exploitation and eventually homing into the optimum, in a rather similar way to the inertia weight in PSO.Consequently, it is modified as: 4 Proposed Approach

Hybrid ICS with Luus-Jakoola Heuristics
To address the IFS inverse problem described in Section 2.2, we propose a new hybrid scheme for proper balance between exploration and exploitation.Firstly, we consider the improved cuckoo search method for global optimization described in Section 3.2.This modified scheme is improved by its hybridization with a local search procedure.In particular, we apply the Luus-Jaakola (LJ) method, a gradient-free heuristics firstly proposed in [15] to solve nonlinear programming problems.LJ starts with an initialization step, where random uniform values are chosen within the search space by computing the upper and lower bounds for each dimension.Then, a random uniform value in-between is sampled for each component.This value is added to the current position of the potential solution to generate a new candidate solution, which replaces the current one only when the fitness improves; otherwise, the sampling space is multiplicatively decreased by a factor (usually of value 95%, but other values can also be used).In practice, we found that it is better to consider a self-adaptive size for this factor, with the effect of speeding up the convergence to the steady state.This process is repeated iteratively.With each iteration, the neighborhood of the point decreases, so the procedure eventually collapses to a point.

Application to the IFS inverse problem
Given a 2D self-similar binary fractal image I comprised of η functions φ κ , we apply our hybrid method to solve the IFS inverse problem.We consider an initial population of χ individuals ØC i Ù i 1,...,χ , where each individual is a collection of η real-valued vectors C i κ of the free variables of Eq. ( 1), as: These individuals are initialized with uniform random values in Ö¡1, 1× for the variables in Θ κ and Σ κ , and in Ö0, 1× for the ω i κ , such that η κ 1 ω i κ 1.After this initialization step, we compute the contractive factors µ κ and reinitialize all functions φ κ with µ κ 1 to ensure that only contractive functions are included in the initial population.Before applying our method, we also need to define a suitable fitness function.In this paper we use the Hamming distance: the fractal images are stored as bitmap images of 0s and 1s for a given resolution defined by a mesh size parameter, m s .Then, we divide the number of different values between the original and the reconstructed matrices by the total number of boxes in the image.This yields the normalized similarity error rate index between both images, denoted by S , the fitness function used in this work.

Parameter Tuning
It is well-known that the parameter tuning of metaheuristic methods is troublesome and problem-dependent.The cuckoo search is specially advantageous in this regard, as it depends on only two parameters: the population size, χ, set to χ 100, and the probability p a , calculated taking: p aM 0.5 and p am 0.005 in Eq. (7).Moreover, the method is executed for Λ iterations.In our simulations, we found that Λ 1500 is enough to reach convergence in all cases.Our hybrid method also requires to define suitable values for α M and α m in Eq. ( 8).Similar to [16], in this work they are set to 0.5 and 0.01, respectively.Finally, our method requires to define the mesh size, ν, set to ν 80 in this work.

Experimental Results
Our method has been applied to several examples of fractals.Only three of them are included here because of limitations of space: rotated triangle, Crystal, and the AIAI-2018 fractal (especially created for this paper).They are shown in Figs. 1 to 3, respectively.The input fractal images are displayed twice in each figure, one in red and another with one color for each individual contractive function.As the reader can see, the three fractals are comprised of 3, 4, and 26 functions, respectively.The application of our method to these input images returns the IFS minimizing the S index, which are then used to render the reconstructed fractal images, shown in blue and with different colors as above.
As the reader can see, although the matching is not optimal, our method captures the underlying structure of the fractal images with good visual quality.This is a remarkable result because our initial population is totally random, meaning that their corresponding images at early generations are all totally different from the given image.Our method is able to obtain a final image that replicates well the input image and is also a self-similar fractal for all instances  number of contractive functions η, number of free variables of the optimization problem, and best and mean value of the S index for 20 independent executions of the method.As the reader can see, the similarity index for the best execution yields errors ranging from 20% to 30%, depending on the example.Although the reconstructed images are perfectly recognizable in all cases, this similarity error increases with the complexity of the model, with the third ex-ample exhibiting the larger value.This result is very reasonable because this example requires to solve an optimization problem with many more variables than the other two examples.In addition, the numerical values of the similarity index used in this work indicate that the results are not optimal yet, suggesting that the method might arguably be further improved.
All computations in this paper have been performed on a 2.6 GHz.Intel Core i7 processor with 16 GB. of RAM.The source code has been implemented by the authors in the native programming language of the popular scientific program Matlab version 2015a and using the numerical libraries for fractals in [4,[6][7][8].Regarding the CPU times, they depend on the complexity of the model and its number of contractive functions.In general, we noticed that the method is slow and time-consuming.For illustration, each single execution takes about 25 45 minutes for the first two examples, and more than 1 hour for the third one.

Conclusions and Future Work
In this paper we introduced a new hybrid method to solve the following optimization problem: given any binary self-similar fractal image, the goal is to determine an IFS whose attractor is a good approximation of this input image.This problem, called the IFS inverse problem, is known to be a very difficult multivariate nonlinear continuous optimization problem.In this new hybrid method, a modification of the original cuckoo search method for global optimization called improved cuckoo search (ICS) is coupled with the Luus-Jakoola heuristics for local search.This hybrid methodology is applied to three fractal examples: rotated triangle, crystal, and the AIAI-2018 fractal.The method replicates the input images very well, yielding visually satisfactory results in all cases.
The numerical results show however that our final solutions are not optimal yet, so the method might be improved in several ways.On one hand, we want to modify our fitness function to obtain a better measure of the quality of the reconstructed fractal.On the other hand, we would like to obtain automatically the optimal value of the number of contractive functions for the IFS.We also wish to extend our results to the cases of images that are neither binary nor self-similar.Reducing our CPU times is also one of our future goals in the field.

Fig. 3 .
Fig. 3. Example of the AIAI-2018 fractal (from top to bottom): original image in red, original image with different colors for each contractive function, reconstructed image in blue, and reconstructed image with different colors for each contractive function.

Table 1 .
[18,19]search algorithm via Lévy flights as originally proposed in[18,19].ÔxÕ, x Ôx1, . . ., xDÕ T Generate initial population of n host nests xi Ôi 1, 2, . . ., nÕ while Ôt M axGenerationÕ or (stop criterion) Get a cuckoo (say, i) randomly by Lévy flights Evaluate its fitness Fi Choose a nest among n (say, j) randomly if (Fi FjÕ Replace j by the new solution end A fraction (pa) of worse nests are abandoned and new ones are built via Lévy flights Keep the best solutions (or nests with quality solutions) Rank the solutions and find the current best end while Postprocess results and visualization end

Table 2 .
Numerical results of the best and mean values of the S index for the examples in Figs. 1 to 3 with the proposed method.

Table 2 .
For each fractal example (in rows), the table shows (in columns):