Here's an approach that unions the primitives as rasters, meshes, and smooths.
Data from OP:
showmaze = Uncompress[FromCharacterCode @@ ImageData[Import["https://i.sstatic.net/XVJcP.png"], "Byte"]];
Primitives with an extended inlet and outlet:
prims = CapsuleShape @@@ Cases[showmaze, _Cylinder, ∞];
prims = prims /. {{5., 5., 5.} -> {5.5, 5., 5.}, {1., 1., 1.} -> {1., 0.5, 1.}};
Now let's rasterize the model. To ensure a quality rasterization, I'll rasterize each capsule separately and join them together:
ims = RegionImage[#, {{0.3`, 5.7`}, {0.3`, 5.7`}, {0.3`, 5.7`}}, RasterSize -> 100] & /@ prims;
im = ImageApply[Max, ims]
We'd like the complement of this image:
im = ImageTake[ColorNegate[im], {5, -5}, {5, -5}, {5, -5}]
Now mesh the data. Here I've clipped to show inside:
Show[bmr = ImageMesh[im, Method -> "DualMarchingCubes"], PlotRange -> {{0, 91}, {1, 92}, {0, 91}}]
This looks pretty good, and for applications like 3D printing it’s very much sufficient, but there are some artifacts we could smooth. One approach is to use GraphDiffusionFlow defined here, but I was unable to find parameters that smoothed out the caps nicely.
Instead I've gone ahead and implemented a version of the approach outlined here. The code for this is at the bottom of the post.
cube = smoothMeshRegion[bmr];
This looks a bit better, however the outer part of the cube has very soft edges:
{cube, Show[cube, PlotRange -> {{-1, 91}, {1, 93}, {-1, 91}}]}
We can fix this by clipping:
cube = BoundaryDiscretizeRegion[smoothed, {{1, 91}, {1, 91}, {1, 91}}];
{cube, Show[cube, PlotRange -> {{1, 90}, {2, 91}, {1, 90}}]}
Finally a wireframe view:
BoundaryMeshRegion[cube, MeshCellStyle -> {1 -> {Thin, Black}}, PlotTheme -> "Lines"]
Edit
I couldn't resist printing this maze
and solving it with food coloring.
Code Dump
(* https://pdfs.semanticscholar.org/c04a/52ad1287385b18464b61f190d1888bf95efd.pdf
Note: the paper suggests to work with the face centroids. I saw no discernible difference between using them or just the face vertices, so for speed and ease of implementation I'm going with the former. *)
Options[smoothMeshRegion] = {"FeatureVertices" -> Automatic, "LaplacianMatrixMethod" -> Automatic, "VertexPenalty" -> Automatic};
smoothMeshRegion[mr:(_MeshRegion|_BoundaryMeshRegion), OptionsPattern[]] :=
Block[{coords, cells, n, L, findices, Fdiag, μ, F, m, A, b, At, ncoords},
(* ------------ mesh data ------------ *)
coords = MeshCoordinates[mr];
cells = MeshCells[mr, 2, "Multicells" -> True];
n = Length[coords];
(* ------------ Laplacian matrix ------------ *)
Switch[OptionValue["LaplacianMatrixMethod"],
"UniformWeight", L = UniformWeightLaplacianMatrix[mr],
Automatic|"Contangent"|_, L = CotangentLaplacianMatrix[mr]
];
(* ------------ vertex penalty matrix ------------ *)
{findices, μ} = featureVertices[mr, OptionValue["FeatureVertices"]];
Fdiag = vertexPenalty[coords, OptionValue["VertexPenalty"]];
Fdiag = ReplacePart[Fdiag, Thread[findices -> μ]];
F = DiagonalMatrix[SparseArray[Fdiag]];
m = Length[F];
(* ------------ global matrix ------------ *)
A = Join[L, F];
(* ------------ right hand side ------------ *)
b = ConstantArray[0., {Length[A], 3}];
b[[n + 1 ;; n + m]] = Fdiag * coords;
(* ------------ solve the system ------------ *)
At = Transpose[A];
ncoords = Quiet[LinearSolve[At.A, At.b, Method -> "Pardiso"]];
(* for large enough μ, ensure the feature vertices are truly fixed *)
If[TrueQ[μ >= $μ],
ncoords[[findices]] = coords[[findices]]
];
(* ------------ construct mesh ------------ *)
Head[mr][ncoords, cells]
]
CotangentLaplacianMatrix[mr_] :=
Block[{n, cells, prims, p1, p2, cos, cots, inds, spopt, L},
n = MeshCellCount[mr, 0];
cells = MeshCells[mr, 2, "Multicells" -> True][[1, 1]];
prims = MeshPrimitives[mr, 2, "Multicells" -> True][[1, 1]];
p1 = prims - RotateRight[prims, {0, 1}];
p2 = -RotateRight[p1, {0, 1}];
cos = Total[p1 p2, {3}] Power[Total[p1^2, {3}]*Total[p2^2, {3}], -0.5];
cots = .5Flatten[cos*Power[1 - cos^2, -.5]];
inds = Transpose[{Flatten[cells], Flatten[RotateLeft[cells, {0, 1}]]}];
Internal`WithLocalSettings[
spopt = SystemOptions["SparseArrayOptions"];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
L = SparseArray[inds -> cots, {n, n}],
SetSystemOptions[spopt]
];
L = L + Transpose[L];
L = L - SparseArray[{Band[{1, 1}] -> Total[L, {2}]}];
L
]
UniformWeightLaplacianMatrix[mr_] :=
Block[{C00, W, II},
C00 = Unitize[#.Transpose[#]]&[mr["ConnectivityMatrix"[0, 2]]];
W = SparseArray[-Power[Map[Length, C00["MatrixColumns"]] - 1, -1.0]];
II = IdentityMatrix[Length[C00], SparseArray];
SparseArray[(C00 - II)W + II]
]
$μ = 5.0;
$λ = 0.5;
featureVertices[mr_, fv:Except[{_, _?Positive}]] := featureVertices[mr, {fv, $μ}];
featureVertices[_, {None, μ_}] := {{}, μ};
featureVertices[mr_, {Automatic, μ_}] := {Union @@ Region`InternalBoundaryEdges[mr][[All, 1]], μ}
featureVertices[_, {vinds_, μ_}] /; VectorQ[vinds, IntegerQ] && Min[vinds] <= 1 := {vinds, μ}
featureVertices[_, {_, μ_}] := {{}, μ};
vertexPenalty[coords_, λ:(Automatic|_?NumericQ)] := ConstantArray[If[TrueQ[NonNegative[λ]], λ, $λ], Length[coords]]
vertexPenalty[coords_, vpfunc_] := vpfunc[coords]
