$\begingroup$

======= 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.

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