Neighborhoods: Experimenting with Cyclic Cellular Automata

On candy stripe legs the Spiderman comes, softly through the shadow of the evening sun (Lullaby, The Cure)

Cellular automata are an inmense source of artistical images. Today, I experimented with Cyclic automata, which are ruled with these simple rules:

  • Create a grid of cells.
  • Give a state to each cell randomly; states a numbers between 0 and M-1 (you choose the value of M previously).
  • For each cell, count how many of its neighbouring cells have their state value exactly 1 unit greater than the cell’s state.
  • If the resulting number is greater than a certain thresold (that you also choose previously) increment the state of the cell by 1; if cell state reaches value of M, then you have to put 0 (in other words, you have to add 1 modulus M).
  • Repeat two previous steps a number of times.

A key concept in this algorithm is defining which is the neighborhood of a cell. There are two of them quite famous: Moore and Von Neumann neighborhoods, but you can define your own ones. Once you decide to stop iterating, you can give color to each cell according its final state, and you will obtain images like this one:

If you have a look to the code, you will see the next parameters:

  • neighborhood: the pattern of neighbouring cells
  • range: the depth of neighborhood
  • states: maximum number of states allowed (the M of the algorithm)
  • thresold
  • iter: number of iterations
  • width and height of the grid

Apart from Moore and Von Neumann, I implemented some other neighborhoods. This chart shows some of them. In columns you can find the folowing: M (Moore), N (Von Neumann), Mr (Moore remote), Nr (Von Neumann remote), Cr (Cross), S1 (S-Shape #1), Bl (Blade), C2 (Corners #2) and TM (Tick Mark). In rows, you can find different ranges for each neighborhood, from 1 to 4:

You will find more neighborhood in the code of the experiment, which is here. There are infinite combinations of the previous parameters, and each one results in a different image:

I used again COLOULovers palettes As always, I encourage you to experiment with the code, create your own neighborhoods, and see how they work. Happy New Year 2021!


Spinning on that dizzy edge (Just Like Heaven, The Cure)

This post talks about a generative system called Physarum model, which simulates the evolution of a colony of extremely simple organisms that, under certain environmental conditions, result into complex behaviors. Apart from the scientific interest of the topic, this model produce impressive images like this one, that I call The Death of a Red Dwarf:

You can find a clear explanation of how a physarum model works in this post, by Sage Jenson. A much deeper explanation can be found in this paper by Jeff Jones, from the University of the West of England. Briefly, a physarum model evolves a set of particles (agents), making them move over a surface. Agents turn towards locations with higher concentrations of a pheromone trail. Once they move, they make a deposition of pheromone as well. These are the steps of a single iteration of the model:

  • Sensor stage: Each agent looks to three positions of the trail map (left, front and right) according a certain sensor angle
  • Motor stage: then it moves to the place with the higher concentration with some rules to deal with ties
  • Deposition stage: once in the new location, the agent deposites a certain amount of pheromone.
  • Diffuse stage: the pheromones diffusing over the surface to blur the trail array.
  • Decay stage: this make to decay the concentration of pheromone on the surface.

Sage Jenson explains the process with this illustrative diagram:

My implementation is a bit different from Jones’ one. The main difference is that I do not apply the diffuse stage after deposition: I prefer a high defined picture instead blurriyng it. I also play with the initial arrangement of agents (location and heading angle) as well with the initial configuration of environment. For example, In The Death of a Red Dwarf, agents start from a circle and the environment is initialized in a dense disc. You can find the details in the code. There you will see that the system is governed by the following parameters:

  • Front and left sensor angles from forward position
  • Agent rotation angle
  • Sensor offset distance
  • Step size (how far agent moves per step)
  • Chemoattractant deposition per step
  • Trail-map chemoattractant diffusion decay factor

In adition to them (specific of the physarum model), you also can change others like colors, noise of angles (parameter amount of jitter function), number of agents and iterations as well as the initial arrangement of the environment and the location of agents. I invite you to do it. You will discover many abstractions: butterfly wings, planets, nets, explosions, supernovas … I have spent may hours playing with it. Some examples:

You can find the code here. Please, let me know if you do something interesting with it. Share your artworks with me in Twitter or drop me an email (you can find my address here).


Moça do corpo dourado
Do Sol de Ipanema
O seu balançado
É mais que um poema
(Garota de Ipanema, João Gilberto)

Sometimes I think about the reasons why I spend so many time doing experiments and writing my discoveries in a blog. Even although the main reason to start this blog was some kind of vanity, today I have pretty clear why I still keep writing it: to keep my mind tuned. I really enjoy looking for ideas, learning new algorithms, figuring out the way to translate them into code and trying to discover new territories going a step further. I cannot imagine my life without coding. Many good times in the last years have been in front of my laptop listening music and drinking a beer. In these strange times, confined at house, coding has became in something more important. It keeps me ahead from the sad news and moves my mind to places where everything is quiet, friendly and perfect. Blogging is my therapy, my mindfulness.

This post is inspired in this post from Softology, an amazing blog I recommend you to read. In it, you can find a description of the stepping stone cellular automaton as well as a appealing collection of images generated using this technique. I modified the original algorithm described in the post to create images like these, which remind me a watercolor painting:

I begin with a 400 x 400 null matrix. After that, I choose a number of random pixels that will act as centers of circles. Around them I substitute the initial zeros by numbers drawned from a normal distribution which mean depends on the distance of pixels to the center. The next step is to apply the stepping stone algorithm. For each pixel, I substitute its value by a weighted average of itself and the value of some of its neighbors, choosen randomly. I always mix values of the pixels. The original algorithm, as described in the Softology’s blog, performs these mixings randomly. Another difference is that I mix values intead interchanging them, as the original algorithm does. Once I repeat this process a number of times, I pick a nice palette from COLOURLovers and turn values of pixels into colors with ggplot:

The code is here. Let me know if you do something interesting with it. Turning numbers into bright colors: I cannot imagine a better way to spend some hours in these shadowy times.


Ooh, see the fire is sweepin’
Our very street today
Burns like a red coal carpet
Mad bull lost its way
(Gimme Shelter, The Rolling Stones)

After following this easy tutorial, you will be able to create tiled images from a photograph. You may want to use your own portrait or some other as I did. I use geom_tile: one of my preferred geometries of ggplot, which was the one I used in some other experiments like space invaders or Newton’s fractals. I used original photos from some of the most terrific monsters of the cinema: Frankenstein, Dracula and The Mummy. I love how rough squares create textures and sense of depth. This is Frankenstein after the transformation:

The process is quite simple:

  • Load the image and convert it to grayscale. I use imager, a very useful and easy to use package for image processing with R.
  • Reduce the resolution of the image as well as its dimension. Each new (big) pixel is summarized with the mean of the grayscale values of the pixels inside it.
  • Divide these average values into a number of groups using cut function.
  • Represent pixels with ggplot, using geom_tile. There are two important parameters: size and color of lines of each tile. Both of them depend on the value of the group previously calculated in which falls the tile. The graph is composed layer by layer depending on these gropus.

This is Dracula after being tiled:

And this is The Mummy:

These are the stunning original images (I love them all):

You can find the code here: let me know if you do something with it.

Reaction Diffusion

Sin patria ni banderas, ahora vivo a mi manera; y es que me siento extranjero fuera de tus agujeros (Tercer movimiento: Lo de dentro, Extremoduro)

The technique I experimented with in this post is an endless source to obtain amazing images. It is called reaction-diffusion and simulates the evolution of a system where several substances interact chemically transforming into each other (reaction) and spreading out over a surface in space (diffusion). In my case there are just two substances in a 2D space and the evolution of system is simulated using the Gray-Scott algorithm which is ruled by several parameters that, once determined, can produce images like this one:

This article by Karl Sims, is a very clear and comprehensive explanation of the Gray-Scott algorithm. Briefly, the Gray-Scott model simulates the evolution of two substances, A and B in a two dimensional grid. In each cell of the grid, the substance A is added at a given feed rate f. Then, both substances react following this rule: two particles of B convert a particle of A into a particle of B. To prevent overpopulation, B particles are killed at a given kill rate k. In the diffusion phase, both substances spread out accross the cells of the grid at a given diffusion rates Da and Db. The diffusion is performed using a 2D Lapacian operator with a 3x3 convolution matrix L.

In the article you can find the equations that rule the system, which depend on the parameters described in the previous paragraph. To obtain all the images of this post, I maintained most of them always constant and just changed the following:

  • Feed and kill rates (f and k respectively)
  • The initial proportion of both substances A and B. I always started with A=0 in each cell and B=1 in some (small) amount of them selected randomly or according to some geometrical considerations (and inner circle or square, for example). I let you to try with your own ideas.

Sometimes the system converges and remains stable after a number of iterations. For example, these images are obtained iterating 5000 times:

Before converging you can obtain also nice patterns in between:

The variety of patterns is amazing: tessellations, lattices, caleidoscopic … some more examples:

I used again Rcpp to iterate efficiently but now I tried RcppArmadillo, a C++ library for linear algebra and scientific computing because it contains a class called cube, which is a 3D matrix that fits perfectly into a 2D grid where the 3rd dimension is the amount of particles A and B.

I like to share my code because I think it may serve as a source of inspiration to someone. You can find the code here with just one of the many configurations I tried to generate the previous images. It may serve you as a good starting point to explore you own ones. Let me know if you reach interesting places.

Happy New Year 2020.

The Chaos Game: an experiment about fractals, recursivity and creative coding

Mathematics, rightly viewed, possesses not only truth, but supreme beauty (Bertrand Russell)

You have a pentagon defined by its five vertex. Now, follow these steps:

  • Step 0: take a point inside the pentagon (it can be its center if you want to do it easy). Keep this point in a safe place.
  • Step 1: choose a vertex randomly and take the midpoint between both of them (the vertex and the original point). Keep also this new point. Repeat Step 1 one more time.
  • Step 2: compare the last two vertex that you have chosen. If they are the same, choose another with this condition: if it’s not a neighbor of the last vertex you chose, keep it. If it is a neighbor, choose another vertex randomly until you choose a not-neighbor one. Then, take the midpoint between the last point you obtained and this new vertex. Keep also this new point.
  • Step 3: Repeat Step 2 a number of times and after that, do a plot with the set of points that you obtained.

If you repeat these steps 10 milion times, you will obtain this stunning image:

I love the incredible ability of maths to create beauty. More concretely, I love the fact of how repeating extremely simple operations can bring you to unexpected places. Would you expect that the image created with the initial naive algorithm would be that? I wouldn’t. Even knowing the result I cannot imagine how those simple steps can produce it.

The image generated by all the points repeat itself at different scales. This characteristic, called self-similarity, is property of fractals and make them extremely attractive. Step 2 is the key one to define the shape of the image. Apart of comparing two previous vertex as it’s defined in the algorithm above, I implemented two other versions:

  • one version where the currently chosen vertex cannot be the same as the previously chosen vertex.
  • another one where the currently chosen vertex cannot neighbor the previously chosen vertex if the three previously chosen vertices are the same (note that this implementation is the same as the original but comparing with three previous vertex instead two).

These images are the result of applying the three versions of the algorithm to a square, a pentagon, a hexagon and a heptagon (a row for each polygon and a column for each algorithm):

From a technical point of view I used Rcppto generate the set of points. Since each iteration depends on the previous one, the loop cannot easily vectorised and C++ is a perfect option to avoid the bottleneck if you use another technique to iterate. In this case, instead of writing the C++ directly inside the R file with cppFunction(), I used a stand-alone C++ file called chaos_funcs.cpp to write the C++ code that I load into R using sourceCpp().

Some days ago, I gave a tutorial at the coding club of the University Carlos III in Madrid where we worked with the integration of C++ and R to create beautiful images of strange attractors. The tutorial and the code we developed is here. You can also find the code of this experiment here. Enjoy!

Colonizing Franky

Y otra vez me arranco despacito, al sentir que nada necesito (Locura transitoria, Extremoduro)

One of my favorite sites in the Internet is algorithmic botany . It’s always a source of inspiration for me. I recently discovered there the space colonization algorithm, concretely in this paper. Originally, the algorithm was developed to simulate leaf venation patterns as well as the branching structure of trees and it works by simulating the competition for space between growing veins (or branches). Given a initial set of attractor points (3.000 points in my case), and a initial node (also a point located randomly inside the picture) the algorithm performs the next steps iteratively:

  • measure distances between attractors and nodes
  • assign the closest node to each attractor
  • keep just those pairs (node, attractor) which distance is between a minimum and maximum
  • normalize and obtain the unit average vectors for each node
  • create a new set of nodes using previous vectors and a predefined longitude

Once again, I used this image of Boris Karloff as Frankenstein’s monster to experiment with the algorithm. I have to say that I coded my own version of the algorithm which maybe is not the canonical one inspired by the previous paper as well as by this tutorial by the coding challenge. This is the resulting drawing using my preferred parametrization. As you will see in the code, I overimpose several layers of colonizations. I love the resulting images:

Play with the parameters to create your own images. If you improve my code, do not hesitate to do a pull request in Github if you want.

Clustering Frankenstein

Necesito para estar sentado, un arbolito en este descampado (Desarraigo, Extremoduro)

From time to time I come back to experiment with this stunning photograph of Boris Karloff as Frankenstein’s monster. I have done several of them previously: from decomposing it into Voronoi regions, to draw it as a single line portrait using an algorithm to solve the travelling salesman problem. I also used this last technique to do a pencil portrait of the image. Today I will use a machine learning algorithm to reinterpret the monster once again. Concretely, I will use hierarchical clustering to do drawings like this one:

The idea is simple: once loaded the photograph, the first step is to binarize it into a black and white image using thresold function of imager package. After that, a random sample of black points is taken. Here comes the clustering algorithm, which starts measuring the euclidean distance between each pair of points. Then, a hierarchical clustering is done so I can reproduce how points are gathered walking through the resulting dendrogram of the previous clustering, starting from the maximum number of clusters (each cluster is an individual point) and ending with the minimum one (just one cluster with the whole sample). The next image shows and example of this process for a sample of 25 points. The left plot shows the population of points and the right one the way that points are connected once the dendrogram is analyzed following the steps described before:

Applying this technique to a big amount of points (between 2.000 and 5.000) result in very interesting drawings. To make process faster, I used map function from purrr package. To render the graph I use ggplot function with geom_curve. This geometry draws a curve between two points, named (x, y) and (xend, yend) respectively. Among others, there are two important parameters to control its shape: curvature (negative values produce left-hand curves, positive values produce right-hand curves, and zero produces a straight line) and angle (values less than 90 skew the curve towards the start point and values greater than 90 skew the curve towards the end point). Playing with this paramaters, as well as with the sample size, you can generate a wide variety of drawings (note that here only appear the segments, since now I removed the points of their extremes):

You can find the code of this experiment here. If you do something interesting with it, please let me know. Thanks a lot for reading my post.


For a lonely soul, you’re having such a nice time (Nothing in my way, Keane)

In my previous post, I created the P2 Penrose tessellation according to the instructions of this post. Now it’s time to create the P3 tessellation following the same technique I described already. This is the image of the P3 tessellation:

Note that all tiles are rhombuses. I recognize that I like the P3 more than the P2, I do not really why. What about you? Here you have the code to play with it if you want.

Kites and Darts: the Penrose Tiling

Agarrada a mis costillas le cuelgan las piernas (Godzilla, Leiva)

Penrose tilings are amazing. Apart of the inner beauty of tesselations, they have two interesting properties: they are non-periodic (they lack any translational symmetry) and self-similar (any finite region appears an infinite number of times in the tiling). Both characteristics make them a kind of chaotical as well as ordered mathematical object that make them really appealing.

In this experiment I create Penrose tilings. Concretely, the P2 tiling, according to this article from Simon Tatham, that I will follow and provides a perfect explanation of how these tessellations can be constructed. The code is available here, and you can use it to create Penrose tilings like this one:

I will not explain in depth how to build the P2 tiling, since the article I mentioned before does it perfectly. Instead of that, I will give some highlights of the process together with a brief explanation of the code involved in it.

Everything has to do with triangles. Concretely, everything has to do with two types of triangles. To differenciate them I name their sides with numbers. The first triangle has labels 1, 2 and 3 and the other one has labels 1, 2 and 4. Two triangles of type 123 forms a kyte like this:

On the other hand, two triangles of type 124 forms a dart like this:

Actually, kites and darts don’t contain their inner segments so both of them are polygons of 4 sides. The building of a Penrose tiling is an iterative process that begins with 5 kites (i.e. 10 triangles of type 123) gathered like this:

You can start with many other patterns but this one will result in a round shape tiling and I like it. I build the tiling by subdividing triangles as Simon describes in his article. A triangle of type 123 is subdivided into three triangles: two of type 123 and one of type 124. The following image shows a 123 triangle (left) and the result after its division (right):

A triangle of type 124 is subdivided into two triangles: one of type 123 and one of type 124. The following image shows a 124 triangle (left) and the result after its division (right):

In each iteration, all triangles are subdivided according its type. After 5 iterations, the resulting pattern is like this:

To make calculations easier I arranged the data frame following a segment structure, in which the sides of triangles are defined by two coordinates: (x, y) and (xend, yend). The bad side of it is that I have rounding problems after making some iterations. It makes the points that would be the same differs slightly because they come from different triangles. I fix it using a hierarchical clustering and substituting points by its centroids after cutting up the dendogram using a very low thresold. Once this problem is solved I can remove the inner segments of all kites and darts, which are segments of type 3 or 4. Apart of removing them, I join the xx triengles to form 4-sides polygons. All these tasks are done with the function Arrange_df (remember that the code is here). This is the result:

This pattern is quite similar to its previous one but now the data frame is ready to be arranged as a polygon using the function Create_Polygon. At least, I calculate the area of each polygon with the Shoelace formula to create a columns called area which I use to fill polygons with two nice colors.

I hope that these explanations will help you to understand and improve the code as well as to invite you to create your own Penrose tilings.