3D Thinning algorithm in Mathematica
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];
padding = Abs[Subtract @@@ minmax];
vpts = Join[pts, (minmax[[All, 1]] - padding)IdentityMatrix[3], (minmax[[All, 2]] + padding)IdentityMatrix[3]];
dm = DelaunayMesh[vpts];
prims = MeshPrimitives[dm, 3, "Multicells" -> True][[1, 1]];
vnodes = First@*Circumsphere /@ prims;
conn = dm["ConnectivityMatrix"[3, 2]];
adj = conn.Transpose[conn];
vlines = UpperTriangularize[adj, 1]["NonzeroPositions"];
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]];
Thread[Polygon[faces] /. Dispatch[Thread[Range[Length[vnodes]] -> vnodes]]]
]
(* 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];
im3dc = ImagePad[im3d, -bd];
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[
ImagePad[#, -BorderDimensions[#, 0]]&[RegionImage[RegionUnion[
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}]
Addendum
Based on papers of Amenta and Attali, a better approximation can be found, but we need to be able to find intersections of spheres efficiently.
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}},
Module[{cb = nb, adj, pp},
cb[[14]] = 0;
adj = a27;
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}];
adj = adj[[pp, pp]];
adj = adj + Transpose[adj];
adj = adj.adj;
If[Min[adj] > 0,
True,
adj = adj.adj;
Min[adj] > 0
]
],
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},
idata = ImageData[ImagePad[im, 1]];
ii = Replace[iter, Except[_Integer?NonNegative] -> 10^9, {0}];
res = morphologicalThinning3DC[idata, ii];
ImagePad[Image3D[res, "Bit"], -1]
]
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]]]
