Introduction
Although there is much to like in this code, there are a few things that are not really optimal. For instance, randomly combining rotation matrices and hoping you'll eventually get all necessary rotations just doesn't feel good. Also, randomly trying the pieces is not really efficient.
I made a version directly aimed at solving the SOMA puzzle that Anton Antonov referred to above, but which can be easily changed into the original one. This version should be more efficient, and uses a fast method to find suitable locations for every piece you want to fit next in the fitting process.
My version uses a straightforward and classical backtracking procedure using a recursive algorithm where new pieces are added until either the last piece has been successfully placed or placing a new piece is not possible anymore. In that case, the algorithm goes back to the last stage which still had untried alternatives and tries again from there.
The pieces
pieces =
{
{{0, 0, 0}, {1, 0, 0}, {0, 1, 0}},
{{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 2, 0}},
{{0, 0, 0}, {0, 1, 0}, {0, 2, 0}, {1, 1, 0}},
{{0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {1, 2, 0}},
{{0, 0, 0}, {1, 0, 0}, {0, 0, 1}, {0, 1, 1}},
{{0, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 0, 1}},
{{0, 0, 0}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1}}
};
Visualization defs
makeCubes := Cuboid[# - {1, 1, 1}/2, # + {1, 1, 1}/2] &
showPieces := Graphics3D[{#},
PlotRange -> {{-1/2, 2 + 1/2}, {-1/2, 2 + 1/2}, {-1/2, 2 + 1/2}},
Boxed -> True, ImageSize -> Tiny
] &
colors = {RGBColor[0.5, 0.5, 0.6], RGBColor[0.8, 0.14, 0.17],
RGBColor[0.85, 0.81, 0.16], RGBColor[0.24, 0.27, 0.87],
RGBColor[0.21, 0.63, 0.24], RGBColor[0.71, 0.44, 0.72],
RGBColor[0.98, 0.53, 0.14]};
A convenience function to translate all parts of a piece to a new position:
translateAll[positions_List, shift_List] := (shift + #) & /@ positions;
Display the SOMA pieces
MapIndexed[
Graphics3D[Prepend[#1, colors[[#2[[1]]]]], Boxed -> False] &,
Map[makeCubes, pieces, {2}]
] // Row

Rotation handling
Four steps of 90 degrees rotations in each of the orthogonal directions are combined, then duplicate results are removed. Total number of rotations is 24. Additionally, below a convenience function to generate all rotational variations of a given piece. Note that it translates each piece to the origin and duplicates (due to symmetries in the pieces) are removed (the Sort
in the code makes this work correctly).
xr = NestList[RotationMatrix[?/2, {1, 0, 0}].# &, DiagonalMatrix[{1, 1, 1}], 3];
yr = NestList[RotationMatrix[?/2, {0, 1, 0}].# &, DiagonalMatrix[{1, 1, 1}], 3];
zr = NestList[RotationMatrix[?/2, {0, 0, 1}].# &, DiagonalMatrix[{1, 1, 1}], 3];
rots = Flatten[Outer[Dot, xr, yr, zr, 1], 2] // DeleteDuplicates;
doRotations[pointList_] :=
Module[{rotRes},
rotRes = Outer[Dot, rots, pointList, 1];
Function[{coordList}, (# - Min /@ Transpose[coordList]) & /@ coordList // Sort] /@ rotRes // DeleteDuplicates
]
Rotational variations of each of the pieces
rotatedPieces = doRotations /@ pieces;
Showing them:
Framed@*Multicolumn /@ Map[showPieces, Map[makeCubes, rotatedPieces, {3}], {2}] // Multicolumn

Changing object definition
Currently, the pieces are defined by the individual coordinates of their building blocks. However, I will need a different representation further on, where the objects are defined in a 3D matrix, with entries indicating either empty or filled space. The following routines will take care of the conversion:
fillObject[object_, fillPositions_List, filling_Integer] :=
Module[{ob = object},
Do[With[{in = Sequence @@ (i + 1)}, ob[[in]] = filling], {i, fillPositions}]; ob
]
coordinatesToObject[coordList_] :=
Module[ {range},
range = Max /@ Transpose[coordList] + 1;
fillObject[ConstantArray[0, range], coordList, 1]
]
We can now convert the SOMA pieces to this representation:
piecesAsObjects = Map[coordinatesToObject, rotatedPieces, {2}];
Visualizing those with Image3D
:
Map[Image3D, piecesAsObjects, {2}] // Grid

Finding a solution
With everything setup above, finding a solution is now relatively easy. As mentioned in the introduction it will be a backtracking, recursive algorithm.
The core of the algorithm is the use of ListCorrelate
(which was the reason to move to a different object representation). We use it to move a piece over the available volume (also defined in terms of a matrix indicating free and occupied space). The result will be a matrix that will contain a number equal to the size of the piece on all locations where that piece could be placed. We then use Position
to find all those locations. The locations will then be tried successively in finding the solution.
fitPieces[fitVolume_, pieceNr_, sol_] :=
Module[{lookup, possiblePositions, newVolume, placedPiece, newSol},
Do[
lookup = piecesAsObjects[[pieceNr, rotNr]];
possiblePositions = translateAll[Position[Quiet@ListCorrelate[lookup, fitVolume], Total[lookup, 3]], -{1, 1, 1}];
Do[
placedPiece = translateAll[rotatedPieces[[pieceNr, rotNr]], pos];
newSol = Append[sol, {pieceNr, rotNr, pos}];
If[pieceNr == Length[pieces],
Throw[newSol];
(* else *),
newVolume = fillObject[fitVolume, placedPiece, 0];
fitPieces[newVolume, pieceNr + 1, newSol];
],
{pos, possiblePositions}
],
{rotNr, Length@rotatedPieces[[pieceNr]]}
]
]
Note that I used Throw
to escape from the function as soon as a solution is found. A Catch
around the initial call to fitPieces
will be necessary. If you want to find all solutions, replace Throw
with Sow
and Catch
with Reap
. The attached notebook has a demonstration of this at the bottom.
Testing
A cube:
fillSpace = {({
{1, 1, 1},
{1, 1, 1},
{1, 1, 1}
}), ({
{1, 1, 1},
{1, 1, 1},
{1, 1, 1}
}), ({
{1, 1, 1},
{1, 1, 1},
{1, 1, 1}
})}
res = Catch[fitPieces[fillSpace, 1, {}]]
{{1, 1, {0, 0, 1}}, {2, 1, {0, 0, 0}}, {3, 1, {0, 0, 2}}, {4, 10,
{0, 2, 0}}, {5, 7, {1, 0, 0}}, {6, 3, {1, 1, 1}}, {7, 2, {1, 0,
1}}}}
Interpretation of each part: piece number, rotation variant and position in solution. Visualization routine:
showResult[res_] :=
Graphics3D[
MapIndexed[
Prepend[#1, colors[[#2[[1]]]]] &,
Map[makeCubes, (translateAll[rotatedPieces[[#1, #2]], #3] & @@@ res), {2}]
],
Boxed -> False,
SphericalRegion -> True,
Lighting -> "Neutral"
]
showResult[res]

More examples from the SOMA chart
426
fillSpace = {({
{1, 1, 1},
{1, 1, 1},
{1, 1, 1},
{1, 1, 1},
{1, 1, 1},
{1, 1, 1}
}), ({
{1, 1, 0},
{1, 1, 0},
{1, 1, 0},
{1, 1, 0},
{0, 0, 0},
{0, 0, 0}
}), ({
{1, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0}
})};
res = Catch[fitPieces[fillSpace, 1, {}]];
showResult[res]
It took 6 seconds to find this one.
427
fillSpace = {({
{1, 1, 1},
{1, 1, 0},
{1, 0, 0}
}), ({
{1, 1, 1},
{1, 0, 0},
{1, 0, 0}
}), ({
{1, 1, 1},
{1, 1, 0},
{1, 0, 0}
}), ({
{1, 1, 0},
{1, 1, 0},
{0, 0, 0}
}), ({
{1, 1, 0},
{1, 0, 0},
{0, 0, 0}
}), ({
{1, 1, 0},
{1, 0, 0},
{0, 0, 0}
})};
res = Catch[fitPieces[fillSpace, 1, {}]];
showResult[res]

Found in 1.2 seconds.
445
This one is nasty as the picture is misleading, apparently showing 28 cubes instead of 27. One solution is:

Found in 0.1 seconds.
Attachments: