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