Producing a random Wang Tile tiling image more efficiently

I'm following along with this SIGGRAPH 2006 paper Recursive Wang Tiles for Real-Time Blue Noise - there's a video here too. Eventually I want to try to produce the blue noise results in the paper, and experiment with texture synthesis, but I'm just getting started using Wang Tiles at the moment.

Here I'm trying to turn the 8-tile tile set that appears in the video and produce a random valid tiling. I'm using a ResourceFunction["BacktrackSearch"] and an ArrayFilter to check all the tiles match edges, or match to a zero tile / padding that has not yet been filled in.

It's a bit slow. Is there a more efficient way to generate a tiling?

Remove["Global`*"]

(* tile num to colours 1,2,3,4 *)
tileset = <|0 -> {0, 0, 0, 0}, 1 -> {1, 2, 3, 4}, 2 -> {1, 4, 3, 2},
3 -> {3, 4, 3, 4}, 4 -> {1, 4, 1, 4}, 5 -> {3, 4, 1, 2},
6 -> {3, 2, 1, 4}, 7 -> {3, 2, 3, 2}, 8 -> {1, 2, 1, 2}|>;

(* given a tile id, produce an image *)
tileToImage[num_] :=
With[{ts = tileset[num],
cols = {
RGBColor[0.9999785876752382, 0.475513385976072, 0.6629426677719442],
RGBColor[0.9999991551385031, 0.6620017739781776, 0.9994575868819303],
RGBColor[0.7921405246863186, 0.9455041087128487, 0.3831022964631722],
RGBColor[0.5908239188023992, 0.7286211433501548, 0.0005627139280530491]}
},
Rasterize[
Graphics[
MapIndexed[{#1,
Rotate[Triangle[],
45° - 90° *(#2[] - 1), {0, 0}]} &,
cols[[ts]]], ImageMargins -> 0, ImageSize -> Full,
PlotRange -> {{-1, 1}, {-1, 1}}/Sqrt], RasterSize -> {32, 32}]]

(* match a piece with center or zero - zero is just padding and always matches *)
moz[x_, a_] := x == a || x == 0 || a == 0;
(* take the 3x3 block b and match the n,s,e,w quadrants with corresponding opposite sides on neighbours *)
matching[b_] := With[{c = tileset[b[[2, 2]]]
, top = tileset[b[[1, 2]]]
, right = tileset[b[[2, 3]]]
, bottom = tileset[b[[3, 2]]]
, left = tileset[b[[2, 1]]]},
moz[c[], top[]] &&
moz[c[], right[]] &&
moz[c[], bottom[]] &&
moz[c[], left[]]
]

(* Reshape the tile list into a grid. A choice of tiles is valid if for every element its immediate neighbours all match or match with zero *)
valid[tiles_, dimensions_] :=
With[{grid = ArrayReshape[tiles, dimensions]},
Count[ArrayFilter[matching, grid, 1, Padding -> 0], False, 2] == 0]

result = Module[{dimensions = {8, 8}, problemsize},
problemsize = Times @@ dimensions;
ArrayReshape[
(* find a list of tiles we can organize into a dim x dim array *)

ResourceFunction["BacktrackSearch"][
Table[RandomSample@DeleteCases[Keys[tileset], 0], problemsize],
valid[#, dimensions] &,
valid[#, dimensions] && Length[#] == problemsize &
], dimensions]
];

ImageAssemble[Map[tileToImage, result, {2}]] Update:

This works a bit better by precomputing all the tile images and not doing the ArrayFilter on the whole array, but instead only matching up the last tile added with the tile above and the tile to the left. I can now get 16x16 tilings in a more reasonable time frame, though it's still a bit slow for some reason:

(** ^^^ tileset and tileToImage ^^^ **)

(* match a piece with center or zero - zero is just padding and always matches *)
moz[x_, a_] := x == a || x == 0 || a == 0;

(* match the top and left edges up *)
matching[thisone_, above_, left_] :=
With[{t = tileset[thisone], a = tileset[above], l = tileset[left]},
moz[t[], a[]] && moz[t[], l[]]]

lastValid[tiles_, dims_] := With[{
width = Last@dims,
offset = Length@tiles},
matching[tiles[[-1]],
If[offset > width, tiles[[offset - width]], 0],
If[Mod[offset, width] == 1, 0, tiles[[offset - 1]]]
]]

result = Module[{dimensions = {16, 16}, problemsize},
problemsize = Times @@ dimensions;
ArrayReshape[
(* find a list of tiles we can organize into a dim x dim array *)

ResourceFunction["BacktrackSearch"][
Table[RandomSample@DeleteCases[Keys[tileset], 0], problemsize],
lastValid[#, dimensions] &,
Length[#] == problemsize &
], dimensions]
];

tileImages = AssociationMap[tileToImage, Keys@tileset];
ImageAssemble[Map[tileImages[#] &, result, {2}]]
• There is a Wolfram Demonstrations project example here - it's also a bit slow and the source has no comments so it's hard to follow. Dec 1 '21 at 14:03

Haven't looked much to check what ResourceFunction["BacktrackSearch"] does, but this appears similar in spirit to OP's edit (i.e. using horizontal/vertical topology constraints to select the next tile). This is inspired by this excellent community post on tiling constraints and my comment there.

Pre-define the topology constraints (manually in this case, since it was easy - I suppose it can be automated fairly easily during the tileset information):

topologyConstraints["right"]=<|1->{2,5,7,8},2->{1,3,4,6},3->{1,3,4,6},4->{1,3,4,6},5->{1,3,4,6},6->{2,5,7,8},7->{2,5,7,8},8->{2,5,7,8}|>;
topologyConstraints["bottom"]=<|1->{3,5,6,7},2->{3,5,6,7},3->{3,5,6,7},4->{1,2,4,8},5->{1,2,4,8},6->{1,2,4,8},7->{3,5,6,7},8->{1,2,4,8}|>;

We then start by randomly selecting a starting tile on the top-left, and construct the first row and column from that:

{dimx, dimy} = {16, 16};
topLeft = RandomInteger[{1, 8}];
left = NestList[RandomChoice@*topologyConstraints["bottom"], topLeft,dimy - 1];
solution = Transpose[
ReplacePart[Transpose[
Prepend[ConstantArray[0, {dimy - 1, dimx}],
NestList[RandomChoice@*topologyConstraints["right"], topLeft,
dimx - 1]]], 1 -> left]];

Finally, we start from the (2,2) tile and make our way down by randomly selecting a tile at the intersection of the available horizontal/vertical topological constraints:

Do[solution[[i, j]] =
RandomChoice[
Intersection[topologyConstraints["bottom"][solution[[i - 1, j]]],
topologyConstraints["right"][solution[[i, j - 1]]]]], {i, 2,dimy}, {j, 2, dimx}];

ImageAssemble[Map[tileToImage, solution, {2}]] Note: to speed-up the visualization, it suffices to memoize the rastering function, i.e. add tileToImage[num_] := tileToImage[num] = ... to the definition.