Message Boards Message Boards

Tilings and constraint programming

Posted 3 years ago


The goal of this post is to start from images like this example one :


and generate pictures like:




Mathematica at least 12.1 will be required since Mixed Integer programming is used.


Let's take the first image (black and white) as an example.

Let's assume we have a collection of 16 tiles:


The problem to solve is how to place the tiles on the picture so that the gray content of a tile is close to the gray content of the picture below it and at the same time the topological constraints are satisfied.

By topological constraints, I mean that the tiles must be compatible.

This is allowed:


This is forbidden because the dark horizontal band is continuing as a white horizontal band.


We are going to translate this problem into a set of equations on integer variables and with a linear cost function to optimize. The final problem will be solved with the LinearOptimization function from Mathematica 12.1.

Each pixel of the image is encoded by a vector eq1 because there are 16 different tiles in this example. The components of the vector can be either 0 or 1.

This is expressed as:

VectorLessEqual[{0, v}], VectorLessEqual[{v, 1}], v \[Element] Vectors[nbTiles, Integers]

For each pixel, only one tile can be used. We cannot put several tiles on a pixel but only choose one and only one.

If we use the constraint:


then we express that only one tile can be used. Indeed, since the component are integers and equal to 0 or 1, then the only way to satisfy this equation is that one of the components, and only one, is equal to one.

Expressing the topological constraints is similar.

We have equations like:


This equation is describing a relationship between pixel (x,y) and pixel (x+1,y).

The values on left and right side can either be 0 or 1 (at same time). When zero it means : none of the tiles is used. When one, it means one of the tile is used. So, the translation of the equation is:

If the tile 0,3 or 5 are used at pixel (x,y) then the tiles 3,6,9 or 11 must be used at pixel (x+1,y).

(I have not checked if it makes sense with the set of tiles I am using as example. The constraints for those tiles are probably different.).

To describe the topological constraint of the tiles, we have functions like:

rightSide[tileA[{a_,b_,c_}]] :={b,c};

This is giving a key describing the right side of the tile. Those keys are then used in associations to build the topological constraints. The key can be anything so if you want to add new tiles, you can just use the key you want to describe the sides of your tiles.

The generic function rightSide must be extended with new cases when new tiles are added.

Then, we need to express how good the tiles are approximating the original picture.

For this, an error function is created. It is a sum of terms:


It means that if the tile i is selected for pixel (x,y) then the approximation error is Subscript[f, k]

The function averageColor must be extended with new tiles. It returns the color content of a tile : a RGBColor. The code is using a color distance to compute the errors.

That's why the input picture is always converted to RGB and the alpha channel removed.

So, finally we have translated our problem into a set of integer constraints and with a linear cost function to optimize. It is a mixed integer programming problem which can be solved with LinearOptimization.

The tiles must be displayed. It is done by the function tileDraw and the tile is drawn in a square from (0,0) to (1,1) corners. The code is rasterizing those graphics into 50x50 pixel images.

I have had lots of problems with those pictures due to rounding errors ... probably due to the very old GPU on my very old computer. So I tuned the vector code assuming the final tile image is 50x50 pixels. Now the pictures are well aligned, there is no more one row or one column of wrong pixels on the boundary of the tiles.

But this may cause a problem on your configuration. So, if the tile pictures are not rendering correctly on your side, you'll need to tune my vector graphic code again.

If the picture is too big, solving the full mixed integer programming problem may take too long. But we can solve a sub-optimal problem. We divide the picture into sub-pictures and solve the problem on each sub-picture then we recombine the solutions. For it to work : we must add new constraints to express compatibility between the pictures.

For instance, the left side of a picture at (row,col) must be compatible with the right side of the picture at (row,col-1). So, the problem must be solved in a given order so that the constraints can be propagated from one sub-picture to the other.

For some tiles, the sub-optimal solution can be very good from an artistic point of view (the dark tiles below are working well). For other tiles (the smith tiles in the notebook) either the sub-optimal problem cannot always be solved (because the constraints coming from the previous pictures can't be satisfied) or the sub-optimal problem will look bad from time to time.

So this idea of using sub-picture is really dependent on the kind of tiles used. You need to experiment. But it is art after all.

Example of use

First, we get an example picture:

srcImage = ImageCrop[ExampleData[{"TestImage", "Girl3"}], {190, 270}]

The picture is resized, converted to RGB and any alpha channel removed.

imgToAnalyze = 
   RemoveAlphaChannel[ColorConvert[srcImage, "RGB"], Black]], {50, 

For the dark tiles (knots), we decide to only use gray levels. First color is the background of the tiles. Other colors are for the circles and vertical and horizontal bands.

tileData = 
   0.5, 0.5, 0.5], {RGBColor[0., 0., 0.], RGBColor[1., 1., 1.]}];

It gives a total of 16 tiles. The more tiles, the more difficult it is to solve the problem. 16 is ok on my old computer.

tileData["allTilesImg"] // Length

The problem is solved on 15x15 sub pictures.

solution = partitionSolve[tileData, imgToAnalyze, 15];

The final picture is generated from the tiles and the solution.

img = createPict[tileData, solution];

And you'll get:


Have fun ! I hope the vectorial code will not have to be tuned to generate the tile pictures (rounding errors).

The notebook is attached to the post.

16 Replies

Congratulations! Your post was highlighted on the Wolfram's official social media channels. Thank you for your contribution. We are looking forward to your future posts.

POSTED BY: Moderation Team

This is really cool! There is a chapter about something similar to this in a book called Opt Art that you might enjoy. The whole book is about creating art with linear optimization.

Thank you !

Indeed, I got the idea from this book ! But since I don't have access to powerful solvers like Gorubi, I had to improve a bit the initial idea to be able to work with bigger images.

I have coded a few other algorithms from the book and I'll post them to the forum at some point (just need to find time to clean them).

Posted 2 years ago

Looks like since V12.2 Gorubi can be accessed in WL, given a license for it.

enter image description here

POSTED BY: Greg Hurst

Very nice! Thanks for sharing!

POSTED BY: Sander Huisman

Thanks for sharing, this is excellent! I wanted to make a movie of how one would converge to a (suboptimal) solution using e.g. a stochastic hill-climber.

We first import our mystery image, resize and convert to Grayscale to use e.g. the dark tiles:

  Import[""], \
{50, Automatic}], "Grayscale"]

We use the functions from the OP to define the 16 tiles and topological constraints:

tileDataGrayscale = 
   0.5, 0.5, 0.5], {RGBColor[0., 0., 0.], RGBColor[1., 1., 1.]}];
nextHorizontallyAllowedGrayscale = 
   Join @@ Table[
     tileDataGrayscale["topologyH"][[i, 1, j]] -> 
      tileDataGrayscale["topologyH"][[i, 2]], {j, 
      Length[tileDataGrayscale["topologyH"][[1, 1]]]}, {i, 
nextVerticallyAllowedGrayscale = 
   Join @@ Table[
     tileDataGrayscale["topologyV"][[i, 1, j]] -> 
      tileDataGrayscale["topologyV"][[i, 2]], {j, 
      Length[tileDataGrayscale["topologyV"][[1, 1]]]}, {i, 

We define a cost function:

totalPixedlCostGrayscale[tileData_][target_, solution_] := 
   ColorDistance[GrayLevel[#1], tileData["tileColorValue"][[#2]], 
     DistanceFunction -> "CIE2000"] &, {target, solution}, 2], 2]

And make our starting guess:

topLeftGrayscale = RandomInteger[{1, 16}];
leftGrayscale = 
   topLeftGrayscale, 41];
solutionGrayscale = Transpose[
   Prepend[ConstantArray[1, {41, 50}], 
     topLeftGrayscale, 49]]], 1 -> leftGrayscale]];
Do[solutionGrayscale[[i, j]] = 
    nextVerticallyAllowedGrayscale[solutionGrayscale[[i - 1, j]]], 
     solutionGrayscale[[i, j - 1]]]]], {i, 2, 42}, {j, 2, 50}]
initialCostGrayscale = 
 costGrayscale = 
   ImageData[img], solutionGrayscale]
Show[createPict[tileDataGrayscale, solutionGrayscale], 
 ImageSize -> 300]

enter image description here

We write a simple mutation code to randomly flip some of these tiles (subject to the topological constraints). We accept the mutation for all cases where the cost function is less than the current cost function and according to the Metropolis algorithm for small positive cost function deltas (to hopefully avoid getting stuck in local minima):




RandomReal[]<Quiet[Exp[-normalizedCostDelta 5000]],

we initialize our climber and iterate:

countGrayscale = 0;
incrementalSolutionsGrayscale = {solutionGrayscale};
incrementalCostsGrayscale = {{countGrayscale, 
    Rescale[costGrayscale, {linearOptimizationCostGrayscale, 

  Rescale[costGrayscale, {linearOptimizationCostGrayscale, 
Do[mutateGrayscale[], 10^6]

Some 50,000 iterations later we get: enter image description here

For reference, the linear optimization solution is given by: enter image description here

Thanks ! It is awesome. I managed to make it work. But my computer being slow I do not yet have any cool result to share.

I find I keep coming back to this application every time I learn of a new optimization technique!

This time, it's projection onto (non-)convex sets, e.g. the "difference map" algorithm which was previously used to solve Sudoku puzzles (which is of-course another problem amenable to Linear Programming).

Here's an animation of the solver for a 30x30 image (note this seems to be as large as I could make it, the tiling projection is memory intensive..) using 24 grayscale Smith tiles:

enter image description here

A slightly longer walk-through of the projections I used can be found here, and I've attached the source-code. See the full discussion here:

Hope you find watching the algorithm trying to climb out of a local minimum as mesmerizing as I do!


this application is very musical, dependent upon scales and harmonics... here is where the visual and the musical begin to dance...

POSTED BY: Drew Lesso
Posted 3 years ago

Really cool stuff!

Just made this with your code:

enter image description here

and this one:

enter image description here

POSTED BY: Oliver Seipel

I am glad to see that the hack I did for the rounding problems in the vector code is also working on your configuration. The tiles are well aligning.

Thanks for sharing.

There is a related blog by Theodore Gray from 2008 that might be of interest. Different constraints, but similar idea in creating a mosaic.

POSTED BY: Daniel Lichtblau

With some changes to the notebook, the photo mosaic case can be supported:

  • The topological constraints must be disabled
  • A new constraint must be added to prevent a tile from being reused more than once on the whole picture.

If I can find some time, I’ll update the notebook to support this case too.

Amazing idea - and nicely done! It truly generates a piece of art! Thanks for sharing!

POSTED BY: Henrik Schachner

Thank you !

enter image description here -- you have earned Featured Contributor Badge enter image description here Your exceptional post has been selected for our editorial column Staff Picks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract