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.

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

https://reference.wolfram.com/language/ref/method/Gurobi.html

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:

img=ColorConvert[
 ImageResize[
  Import["https://www.beyonddream.com/images/product/23892024.jpg"], \
{50, Automatic}], "Grayscale"]

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

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

We define a cost function:

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

And make our starting guess:

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

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):

mutateGrayscale[]:=Block[{randomRow=RandomInteger[{1,42}],randomCol=RandomInteger[{1,50}],mutation=solutionGrayscale,mutationCost,allowable,normalizedCostDelta},
Do[
allowable=Which[
i==1&&j==1,RandomInteger[{1,16}],
i==1&&j>1,RandomChoice[nextHorizontallyAllowedGrayscale[mutation[[i,j-1]]]],
j==1&&i>1,RandomChoice[nextVerticallyAllowedGrayscale[mutation[[i-1,j]]]],
i>1&&j>1,Intersection[nextVerticallyAllowedGrayscale[mutation[[i-1,j]]],nextHorizontallyAllowedGrayscale[mutation[[i,j-1]]]]];

If[i==randomRow&&j==randomCol,
If[Length[allowable]>1,mutation[[i,j]]=RandomChoice[DeleteCases[allowable,mutation[[i,j]]]],Break[]],

If[MemberQ[allowable,mutation[[i,j]]],Continue[],
mutation[[i,j]]=RandomChoice[allowable]]],
{i,randomRow,42},{j,randomCol,50}];

countGrayscale=countGrayscale+1;
mutationCost=totalPixedlCostGrayscale[tileDataGrayscale][ImageData[ColorConvert[obama,"Grayscale"]],mutation];
normalizedCostDelta=Rescale[mutationCost,{linearOptimizationCostGrayscale,initialCostGrayscale}]-Rescale[costGrayscale,{linearOptimizationCostGrayscale,initialCostGrayscale}];
If[
RandomReal[]<Quiet[Exp[-normalizedCostDelta 5000]],
AppendTo[incrementalCostsGrayscale,{countGrayscale,Rescale[costGrayscale=mutationCost,{linearOptimizationCostGrayscale,initialCostGrayscale}]}];AppendTo[incrementalSolutionsGrayscale,solutionGrayscale=mutation];
]]

we initialize our climber and iterate:

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

Dynamic[{countGrayscale, 
  Rescale[costGrayscale, {linearOptimizationCostGrayscale, 
    initialCostGrayscale}]}]
Do[mutateGrayscale[], 10^6]

Some 50,000 iterations later we get:

For reference, the linear optimization solution is given by:

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:

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: https://community.wolfram.com/groups/-/m/t/2823307

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

Attachments:

difference-map-t...nb

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

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

-- you have earned Featured Contributor Badge Your exceptional post has been selected for our editorial column Staff Picks http://wolfr.am/StaffPicks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!