How can I calculate a jigsaw puzzle cut path?

Question

I want to generate a path to cut an arbitrary shape into a set of jigsaw puzzle pieces.

1. All pieces must be unique to preclude placing a piece in the wrong spot.
2. Pieces must be interlocking such that each piece is held by adjacent pieces.
3. It must be possible to generate different paths (sets) for a specific shape that are not merely rotations or reflections of the first.

Since jigsaw puzzles are mass produced there must be a known solution. However, my interest is not in simply implementing this solution, rather I am curious to know how one might approach this task, aided by Mathematica. I am also interested in alternatives to the standard puzzle piece such as these (even if pieces are not 100% interlocking):

Answer 1

======= Update =========

Great question! It inspired this Wolfram Blog article and includes most of the code below plus some apps and fractal layouts like this:

I think it make sense to keep the older code blow for archival and historic purposes.

======= Older implementation =========

Excellent motivating creativity question. This is a bit big for a comment, so here are a few thoughts.

• There is an obvious relation to tiling problems which are well represented in the Wolfram Demonstration Project.

• One approach could be to morph some non-interlocking tilling pieces to have interlocking parts (if interlocking is important)

• Perhaps some twist to a puzzle can come from

  • a-periodic tilling
  • large piece set tilling

• Not all tilings are appropriate, because of, for example, presence of gaps. So here are some candidates:

  • Robinson Tiling
  • Penrose Tiles
  • Pinwheel Tiling
  • Ammann Tiles
  • Pentagon Tilings
  • Voronoi Image
  • Neat Alternating Tilings of the Plane
  • Nowhere-Neat Tilings of the Plane, Part 2

==== Voronoi practical implementation ====

Let's start from writing the following function:

bsc[p1_, p2_] := 
 With[{rc = RandomChoice[{-1, 1}], d = EuclideanDistance[p1, p2], 
   pm = (p1 + p2)/2, dp = (p2 - p1)/5},
  If[d < .1,
   Line[{p1, p2}],
   BSplineCurve[{p1, pm, pm - dp + rc {1, -1} Reverse[dp], 
     pm + dp + rc {1, -1} Reverse[dp], pm, p2}, 
    SplineWeights -> {1, 5, 5, 5, 5, 1}/22]
   ]]

which will morph a long enoug line into a line with a "tongue". It will put the tongue in a random direction for more random generation of puzzle pieces. And some comparison function that will show what points we are adding to wrap BSplineCurve around.

f[p1_, p2_] := 
 With[{d = EuclideanDistance[p1, p2], pm = (p1 + p2)/2, 
   dp = (p2 - p1)/5},
  If[d < .1,
   Line[{p1, p2}],
   Line[{p1, pm, pm - dp + {1, -1} Reverse[dp], 
     pm + dp + {1, -1} Reverse[dp], pm, p2}]
   ]]

Here is a Manipulate to test it out:

Manipulate[
 Graphics[{f @@ pt, {Red, Thick, bsc @@ pt}}, ImageSize -> {300, 300},
   Axes -> True, Frame -> True, AspectRatio -> 1, 
  PlotRange -> {{0, 1}, {0, 1}}], {{pt, {{0, 0}, {1, 1}}}, Locator}, 
 FrameMargins -> 0]

Now this will create a simple Voronoi diagram:

gr = ListDensityPlot[RandomReal[{}, {35, 3}], InterpolationOrder -> 0,
   Mesh -> All, Frame -> False]

And this will extract lines out of it and replace long enoug lines with our tongues function:

Graphics[bsc @@@ 
  Union[Sort /@ 
    Flatten[Partition[#, 2, 1] & /@ 
      Map[gr[[1, 1, #]] &, 
       Flatten[Cases[gr, Polygon[_], Infinity][[All, 1]], 1]], 1]]]

This can be superimposed on an image or simply colorized (with added outer frame):

MorphologicalComponents[
  Binarize@Graphics[{Thick, 
     bsc @@@ Union[
       Sort /@ Flatten[
         Partition[#, 2, 1] & /@ 
          Map[gr[[1, 1, #]] &, 
           Flatten[Cases[gr, Polygon[_], Infinity][[All, 1]], 1]], 
         1]]}, Frame -> True, FrameTicks -> False, 
    PlotRangePadding -> 0, FrameStyle -> Thick]] // Colorize

You have to execute code a few times to find best random colorization - it is based on random tongue orientation.

=== Yet another way - Hilbert & Moore curves====

I slightly modified this Demonstration and cut resulting curves with grid lines:

LSystem[axiom_, rules_List, n_Integer?NonNegative, False] := 
  LSystem[axiom, rules, n, False] = 
   Nest[StringReplace[#, rules] &, axiom, n];
LSystem[axiom_, rules_List, n_Integer?NonNegative, True] := 
  LSystem[axiom, rules, n, True] = 
   NestList[StringReplace[#, rules] &, axiom, n];
LSeed["Hilbert curve"] = "+RF-LFL-FR+"; 
LSeed["Moore curve"] = "+LFL+F+LFL-";
LRules["Hilbert curve"] = {"L" -> "+RF-LFL-FR+", 
  "R" -> "-LF+RFR+FL-"}; 
LRules["Moore curve"] = {"L" -> "-RF+LFL+FR-", "R" -> "+LF-RFR-FL+"};
LPoints[lstring_String] := 
  LPoints[lstring] = Map[First, Split[Chop[Map[First,
       FoldList[Function[{pta, c},
         Switch[c,
          "+", {pta[[1]], pta[[2]] + 90 Degree},
          "-", {pta[[1]], pta[[2]] - 90 Degree},

          "F", {{pta[[1]][[1]] + Cos[pta[[2]]], 
            pta[[1]][[2]] + Sin[pta[[2]]]}, pta[[2]]},
          _, pta]],
        {{0, 0}, 0.},
        Characters[lstring]]]]]];
LPoints[lstring_String, level_Integer?NonNegative] := 
  Map[(2^(5 - level)*# + ((2^(5 - level) - 1)/2)) &, LPoints[lstring]];
LLine[lstring_String, 
   level_Integer?NonNegative] := {AbsoluteThickness[2*(5 - level)], 
   BSplineCurve[LPoints[lstring, level]]};
Manipulate[
 MorphologicalComponents[
   Binarize@Graphics[If[showPreviousLevels == False,
      LLine[
       LSystem[LSeed[whichcurve], LRules[whichcurve], n - 1, 
        showPreviousLevels], n], 
      MapIndexed[LLine[#1, First[#2]] &, 
       LSystem[LSeed[whichcurve], LRules[whichcurve], n - 1, True]]], 
     GridLinesStyle -> Directive[Thick, Black], FrameStyle -> Thick, 
     GridLines -> {None, Table[x, {x, 4, 30, 4}]}, Frame -> True, 
     FrameTicks -> False, PlotRangePadding -> 0]] // Colorize
 , {{n, 4}, 
  ControlType -> None}, {{whichcurve, "Hilbert curve", 
   ""}, {"Hilbert curve", "Moore curve"}, 
  ControlType -> RadioButtonBar}, {{showPreviousLevels, False}, 
  ControlType -> None}, SaveDefinitions -> True]

Answer 2

The first part of the problem is partitioning a shape into smaller parts of a roughly equal area. Then we can add little "tongues" on the pieces to make them interlock.

One idea for partitioning is using either a Delaunay triangulation of a set of points (for triangular pieces) or a Voronoi tessellation (for many-sided polygons).

Let's take for example this method of generating a set of random points with a minimum distance. I modified the algorithm a little to squeeze as many points into a region as possible:

canvas = Image@ConstantArray[0, {100, 100}];
distance = 15;

{img, {pts}} = Reap[
   NestWhile[
    ImageCompose[#, SetAlphaChannel[#, #] &@Image@DiskMatrix[distance],
      Sow@RandomChoice@Position[
         Transpose@ImageData[#, DataReversed -> True], 0.]] &,
    canvas,
    Count[Flatten@ImageData@Binarize[#], 0] > 0 &]];

Then the Delaunay triangulation looks like this:

<<ComputationalGeometry`
PlanarGraphPlot[pts, LabelPoints -> False]

(For a better result we'd need to add points from the edges.)

The Voronoi tessellation looks like this:

Show[DiagramPlot[pts], PlotRange -> 100 {{0, 1}, {0, 1}}]

It can also be shown with ListDensityPlot:

showTiles[pts_] := 
 ListDensityPlot[ArrayPad[pts, {{0, 0}, {0, 1}}], Mesh -> All, 
  InterpolationOrder -> 0, Frame -> False, ColorFunction -> (White &)]

showTiles[pts]

Another possibility is to start with a regular grid of points (e.g. a hexagonal grid),

hex = Join @@ Table[{x, Sqrt[3] y}, {x, 0, 4}, {y, 0, 2}];
pts = Join[hex, TranslationTransform[{1/2, Sqrt[3]/2}] /@ hex];

showTiles[pts]

and distort it randomly:

pts = RandomReal[0.1 {-1, 1}, Dimensions[pts]] + pts;
showTiles[pts]

Just for some fun, we can actually create the tiles and shuffle them a bit. Who wants to add some code to make them draggable and rotatable?

Graphics[
 {EdgeForm[Black],
  Texture@ExampleData[{"TestImage", "Sailboat"}], 

  GeometricTransformation[#, 
     Composition[
      TranslationTransform@RandomReal[0.1 {-1, 1}, 2], 
      RotationTransform[RandomReal[{-Pi, Pi}], Mean@First[#]]]] & /@ 
   Cases[Normal@showTiles[Rescale[pts]], 
    Polygon[p_, ___] :> Polygon[p, VertexTextureCoordinates -> p], 
    Infinity]}
 ]