# 3D Thinning algorithm in Mathematica

as far as I understood Thinning[] only works on 2D images thus far. Is there a library/piece of code that performs 3D Thinning in Mathematica?

Max

edit: I want lines

• You should clarify though what sort of result you want: lines or sheets. In 3D, both are possible. I am wondering if one could somehow post-process the output of DistanceTransform, which does work in 3D. Nov 22, 2018 at 16:21
• Something like this: i.stack.imgur.com/ICbYL.png This is not yet answer-quality, but it is still useful to share. Nov 22, 2018 at 16:26
• This is not completely unrelated... mathematica.stackexchange.com/q/9928/1089 Nov 27, 2018 at 6:07

One approach is to apply a topology preserving morphological thinning method. I'll implement an algorithm by Lee. The main idea is to look at every 3x3x3 cube in the image and to decide which pixels are valid to delete. For the most part I stuck to his approach, except I did not implement the function Octree_Labeling. It's wildly recursive and Compile does not allow for recursion. Instead I took a transitive closure approach.

The compiled routines are quite messy looking since a lot of things are baked into precomputed lookup tables:

With[{
inds = {25,26,16,17,22,23,13,27,24,18,15,26,23,17,19,22,10,13,20,23,11,21,24,20,23,12,15,11,7,16,8,17,4,13,5,9,8,18,17,6,5,15,1,10,4,13,2,11,5,3,2,12,11,6,5,15},
δG6 = Riffle[{1,-1,-1,1,-3,-1,-1,1,-1,1,1,-1,3,1,1,-1,-3,-1,3,1,1,-1,3,1,-1,1,1,-1,3,1,1,-1,-3,3,-1,1,1,3,-1,1,-1,1,1,-1,3,1,1,-1,1,3,3,1,5,3,3,1,-1,1,1,-1,3,1,1,-1,-7,-1,-1,1,-3,-1,-1,1,-1,1,1,-1,3,1,1,-1,-3,-1,3,1,1,-1,3,1,-1,1,1,-1,3,1,1,-1,-3,3,-1,1,1,3,-1,1,-1,1,1,-1,3,1,1,-1,1,3,3,1,5,3,3,1,-1,1,1,-1,3,1,1,-1},ConstantArray[0,128]]
},

preservesEulerCharacteristicQ = Compile[{{nb, _Integer, 1}},
Module[{n},
n = Partition[nb[[inds]], 7].{128,64,32,16,8,4,2};
Total[δG6[[n + 1]]] == 0
],
RuntimeOptions -> "Speed",
CompilationTarget -> "C"
];

];

With[{a27 = ConstantArray[0, {27, 27}], tups = Tuples[{-1,0,1}, 3], dig = IntegerDigits[Range[0, 26], 3, 3] + 1, rng = Range[27]},
With[{c1 = Select[Transpose[Transpose[tups] + #], 0 < Min[#] <= Max[#] < 4&].{9,3,1}-12& /@ dig},
With[{c = Table[Select[c1[[i]], i <= # <= 27 && i != 14&], {i, 27}]},
With[{conn1=c[[1]],conn2=c[[2]],conn3=c[[3]],conn4=c[[4]],conn5=c[[5]],conn6=c[[6]],conn7=c[[7]],conn8=c[[8]],conn9=c[[9]],conn10=c[[10]],conn11=c[[11]],conn12=c[[12]],conn13=c[[13]],conn15=c[[15]],conn16=c[[16]],conn17=c[[17]],conn18=c[[18]],conn19=c[[19]],conn20=c[[20]],conn21=c[[21]],conn22=c[[22]],conn23=c[[23]],conn24=c[[24]],conn25=c[[25]],conn26=c[[26]],conn27=c[[27]]},

simplePointQ = Compile[{{nb, _Integer, 1}},
cb[[14]] = 0;

If[Total[cb] > 23, Return[True]];

If[cb[[1]] == 1, Do[If[cb[[j]] == 1, adj[[1, j]] = 1], {j, conn1}]];
If[cb[[2]] == 1, Do[If[cb[[j]] == 1, adj[[2, j]] = 1], {j, conn2}]];
If[cb[[3]] == 1, Do[If[cb[[j]] == 1, adj[[3, j]] = 1], {j, conn3}]];
If[cb[[4]] == 1, Do[If[cb[[j]] == 1, adj[[4, j]] = 1], {j, conn4}]];
If[cb[[5]] == 1, Do[If[cb[[j]] == 1, adj[[5, j]] = 1], {j, conn5}]];
If[cb[[6]] == 1, Do[If[cb[[j]] == 1, adj[[6, j]] = 1], {j, conn6}]];
If[cb[[7]] == 1, Do[If[cb[[j]] == 1, adj[[7, j]] = 1], {j, conn7}]];
If[cb[[8]] == 1, Do[If[cb[[j]] == 1, adj[[8, j]] = 1], {j, conn8}]];
If[cb[[9]] == 1, Do[If[cb[[j]] == 1, adj[[9, j]] = 1], {j, conn9}]];
If[cb[[10]] == 1, Do[If[cb[[j]] == 1, adj[[10, j]] = 1], {j, conn10}]];
If[cb[[11]] == 1, Do[If[cb[[j]] == 1, adj[[11, j]] = 1], {j, conn11}]];
If[cb[[12]] == 1, Do[If[cb[[j]] == 1, adj[[12, j]] = 1], {j, conn12}]];
If[cb[[13]] == 1, Do[If[cb[[j]] == 1, adj[[13, j]] = 1], {j, conn13}]];
If[cb[[15]] == 1, Do[If[cb[[j]] == 1, adj[[15, j]] = 1], {j, conn15}]];
If[cb[[16]] == 1, Do[If[cb[[j]] == 1, adj[[16, j]] = 1], {j, conn16}]];
If[cb[[17]] == 1, Do[If[cb[[j]] == 1, adj[[17, j]] = 1], {j, conn17}]];
If[cb[[18]] == 1, Do[If[cb[[j]] == 1, adj[[18, j]] = 1], {j, conn18}]];
If[cb[[19]] == 1, Do[If[cb[[j]] == 1, adj[[19, j]] = 1], {j, conn19}]];
If[cb[[20]] == 1, Do[If[cb[[j]] == 1, adj[[20, j]] = 1], {j, conn20}]];
If[cb[[21]] == 1, Do[If[cb[[j]] == 1, adj[[21, j]] = 1], {j, conn21}]];
If[cb[[22]] == 1, Do[If[cb[[j]] == 1, adj[[22, j]] = 1], {j, conn22}]];
If[cb[[23]] == 1, Do[If[cb[[j]] == 1, adj[[23, j]] = 1], {j, conn23}]];
If[cb[[24]] == 1, Do[If[cb[[j]] == 1, adj[[24, j]] = 1], {j, conn24}]];
If[cb[[25]] == 1, Do[If[cb[[j]] == 1, adj[[25, j]] = 1], {j, conn25}]];
If[cb[[26]] == 1, Do[If[cb[[j]] == 1, adj[[26, j]] = 1], {j, conn26}]];
If[cb[[27]] == 1, Do[If[cb[[j]] == 1, adj[[27, j]] = 1], {j, conn27}]];

pp = Complement[rng cb, {0}];

True,
]
],
RuntimeOptions -> "Speed",
CompilationTarget -> "C"
];

]]]];

morphologicalThinning3DC = Compile[{{oim, _Integer, 3}, {iter, _Integer}},
Block[{cnt = 0, X, Y, Z, im = oim, im2 = oim, samefaces = 0, nb, bag, nface, yo, sbp, il, jl, kl},
{X, Y, Z} = Dimensions[im];
While[cnt++ < iter && samefaces < 6,
samefaces = 0;
Do[
nface = 1;
Do[
If[im[[i,j,k]] == 0, Continue[]];

If[!Or[
f == 1 && im[[i,j-1,k]] == 0,
f == 2 && im[[i,j+1,k]] == 0,
f == 3 && im[[i+1,j,k]] == 0,
f == 4 && im[[i-1,j,k]] == 0,
f == 5 && im[[i,j,k+1]] == 0,
f == 6 && im[[i,j,k-1]] == 0
], Continue[]];

nb = Flatten[im[[i-1;;i+1, j-1;;j+1, k-1;;k+1]]];

If[Total[nb] != 2 && preservesEulerCharacteristicQ[nb] && simplePointQ[nb],
If[!simplePointQ[Flatten[im2[[i-1;;i+1, j-1;;j+1, k-1;;k+1]]]],
im2[[i, j, k]] = 1,
im2[[i, j, k]] = 0;
nface = 0;
];
],

{k, 2, Z-1},
{j, 2, Y-1},
{i, 2, X-1}
];
im = im2;
samefaces += nface,
{f, 1, 6}
];
];

im
],
RuntimeOptions -> "Speed",
CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True},
CompilationTarget -> "C"
];


The main routine:

MorphologicalThinning3D[im_, 0] := im

MorphologicalThinning3D[im_Image3D, iter_:Infinity] /; ImageType[im] === "Bit" :=
Block[{idata, ii, res},
ii = Replace[iter, Except[_Integer?NonNegative] -> 10^9, {0}];

res = morphologicalThinning3DC[idata, ii];

]


Some examples:

flat = Import["https://i.stack.imgur.com/heX24.png"];
im = Image3D[First[ImagePartition[flat, {56, 76}]], "Bit"];
thin = MorphologicalThinning3D[im]; // AbsoluteTiming

{1.07601, Null}

{im, thin}


Show[ImageMesh[thin, BaseStyle -> Red], ImageMesh[im, BaseStyle -> Opacity[.3]]]


The second argument represents how many iterations of thinning to perform. The default is Infinity:

frames = Table[MorphologicalThinning3D[im, i], {i, 0, 8}];


Other images:

im = FillingTransform[Closing[ImageCrop[Binarize[ExampleData[{"TestImage3D", "Orbits"}]]], 1]];
thin = MorphologicalThinning3D[im]; // AbsoluteTiming

{3.88302, Null}

{im, thin}


Show[ImageMesh[thin, BaseStyle -> Red], ImageMesh[im, BaseStyle -> Opacity[.3]]]


im = ImagePad[ImageCrop[Binarize[ExampleData[{"TestImage3D", "CTengine"}]]], 1];
thin = MorphologicalThinning3D[im]; // AbsoluteTiming

{33.3316, Null}

{im, thin}


Show[ImageMesh[thin, BaseStyle -> Red], ImageMesh[im, BaseStyle -> Opacity[.3]]]


• Great work, thank you +1
– UDB
Dec 3, 2018 at 15:45

I'm currently working on a fast/scalable solution to question 18893 and it's at a place where I think I can provide an answer here. Though I will admit there's a lot of room for improvement.

In this answer I'll refer to

• the medial surface of a 3D object $$O$$ as the surface comprising the points equidistant from at least 2 points on $$\partial O$$
• the medial axis of a 3D object as the singular (or non-manifold) edges of the medial surface.

The main idea of this answer is that for a sufficiently dense sampling of points on the boundary of a volume, the vertices, edges, and faces of the Voronoi diagram gives an approximation to the medial surface. Since OP is interested in lines, I will use this to find the medial axis.

— Since the medial surface is geometrically unstable, this answer might be a dead end. —

First here's the code to find Voronoi faces for a set of points:

VoronoiFaces3D[pts_] :=
Block[{minmax, padding, vpts, dm, prims, vnodes, conn, adj, vlines, mr1d, linepairs, normals, planar, faces},
minmax = MinMax /@ Transpose[pts];

dm = DelaunayMesh[vpts];
prims = MeshPrimitives[dm, 3, "Multicells" -> True][[1, 1]];
vnodes = First@*Circumsphere /@ prims;

conn = dm["ConnectivityMatrix"[3, 2]];

mr1d = Quiet @ MeshRegion[vnodes, Line[vlines]];

conn = mr1d["ConnectivityMatrix"[0, 1]];
adj = Unitize[conn.Transpose[conn]] - IdentityMatrix[Length[conn], SparseArray];
linepairs = Catenate[Insert[#1, 2] /@ Subsets[{##2}, {2}]& @@@ adj["MatrixColumns"]];

normals = lpNormals[vnodes, linepairs];
planar = GatherByList[linepairs, Round[normals, 10^-8.]];

(* this is might not work for adjacent coplanar faces... *)
faces = FindFundamentalCycles[Apply[UndirectedEdge, Catenate[Union @@@ validEdges[planar, Range[Length[planar]]]], {1}]][[All, All, 1, 1]];

]

(* function introduced below *)
GatherByList[list_, representatives_] := Module[{func},
func /: Map[func, _] := representatives;
GatherBy[list, func]
]

lpNormals = Compile[{{coords, _Real, 2}, {lp, _Integer, 1}},
Module[{cross, sgn},
cross = Cross[coords[[lp[[2]]]] - coords[[lp[[1]]]], coords[[lp[[3]]]] - coords[[lp[[1]]]]];
cross = cross/Sqrt[cross.cross];
sgn = Sign[First[cross]];
If[sgn == 0, sgn = Sign[cross[[2]]]];
If[sgn == 0, sgn = Sign[cross[[3]]]];

sgn*cross
],
RuntimeAttributes -> {Listable},
CompilationTarget -> "C",
Parallelization -> True,
RuntimeOptions -> "Speed"
];

validEdges = Compile[{{inds, _Integer, 1}, {k, _Integer}},
Module[{i1, i2, i3},
i1 = inds[[1]];
i2 = inds[[2]];
i3 = inds[[3]];
Which[
i1 > i2 && i3 > i2, {{{i2, k}, {i1, k}}, {{i2, k}, {i3, k}}},
i2 > i1 && i3 > i2, {{{i1, k}, {i2, k}}, {{i2, k}, {i3, k}}},
i1 > i2 && i2 > i3, {{{i2, k}, {i1, k}}, {{i3, k}, {i2, k}}},
True, {{{i1, k}, {i2, k}}, {{i3, k}, {i2, k}}}
]
],
RuntimeAttributes -> {Listable},
CompilationTarget -> "C",
Parallelization -> True,
RuntimeOptions -> "Speed"
];


GatherByList introduced here.

Since I have yet to answer 18893 (with this code), here's an example:

SeedRandom[10000];
pts = RandomReal[{0, 1}, {1000, 3}];

faces = VoronoiFaces3D[pts];
Graphics3D[MapIndexed[{Lighter[ColorData[112] @@ #2, .3], #1} &, faces], PlotRange -> {{0, 1}, {0, 1}, {0, 1}}, Lighting -> {{"Ambient", White}}, Boxed -> False]


The following function takes in a binary Image3D object, converts it into a BoundaryMeshRegion, approximates the medial surface, finds the approximate medial axis, then converts back into an image. To approximate the medial surface, I simply take the Voronoi faces completely contained in the object.

MedialAxisThinning3D[im3d_Image3D, testpts_:1000, smallsize_:Automatic, seed_:10000] /; ImageType[im3d] === "Bit" :=
Block[{bd, im3dc, bmr, pts, faces, rmf, ms, se, ma, medialaxis},
bd = BorderDimensions[im3d, 0];
bmr = ImageMesh[im3dc, Method -> "DualMarchingCubes"];

BlockRandom[
SeedRandom[seed];
pts = RandomPoint[RegionBoundary[bmr], testpts];
];
faces = VoronoiFaces3D[pts];
rmf = RegionMember[bmr];

ms = DiscretizeGraphics[Polygon[Select[faces[[All, 1]], And @@ rmf[#]&]]];
se = Pick[Range[MeshCellCount[ms, 1]], UnitStep[2 - Differences[ms["ConnectivityMatrix"[1, 2]]["RowPointers"]]], 0];
ma = MeshRegion[MeshCoordinates[ms], MeshCells[ms, {1, se}]];

medialaxis = ImageResize[
ma,
Point[Tuples[RegionBounds[bmr]]]
]]],
ImageDimensions[im3dc]
];

ImagePad[DeleteSmallComponents[Binarize[medialaxis], Replace[smallsize, Automatic -> Sequence[], {0}], CornerNeighbors -> False], bd]
]


We can see there are some heuristic parameters. I haven't found a way to automatically determine appropriate values.

Note too that this routine can be slow. I think there's a lot of room for improvement.

Some examples:

im3d = ImageCrop[Binarize[FillingTransform[Closing[ExampleData[{"TestImage3D", "Orbits"}], 1]]]];

medialaxis = MedialAxisThinning3D[im3d];

{im3d, medialaxis}


bd = ImagePad[EdgeDetect[ImagePad[im3d, 1]], -1];
pp = PixelValuePositions[bd, 1];
towhite = Pick[pp, UnitStep[Transpose[pp][[1]] - 40], 1];
slice = im3d*Dilation[ReplacePixelValue[bd, towhite -> 0], 1];
{ColorCombine[{bd, medialaxis, .7 medialaxis}], ColorCombine[{slice, medialaxis, .7 medialaxis}]}


im3d = ImageCrop[Binarize[ExampleData[{"TestImage3D", "CTengine"}]]];

(* this is slow... *)
medialaxis = MedialAxisThinning3D[im3d, 10000, 40];

{im3d, medialaxis}


ColorCombine[{EdgeDetect[im3d], medialaxis, .7 medialaxis}]