# 3D Thinning algorithm in Wolfram Language

**URL:** <https://community.wolfram.com/t/3d-thinning-algorithm-in-wolfram-language/13370>\
**Category:** Staff Picks\
**Tags:** geometry, graphics-and-visualization, graphs-and-networks, optimization, 3d-printing, staff-pick, wolfram-repositories\
**Created:** [November 28, 2018, 6:08pm UTC](https://community.wolfram.com/t/3d-thinning-algorithm-in-wolfram-language/13370 "2018-11-28T18:08:29Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![ghurst](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/ghurst/32/20267_2.png) [@ghurst](https://community.wolfram.com/u/ghurst)\
**Post date:** [November 28, 2018, 6:08pm UTC](https://community.wolfram.com/t/3d-thinning-algorithm-in-wolfram-language/13370/1 "2018-11-28T18:08:29Z")

</div>

**[Open in Cloud](https://www.wolframcloud.com/objects/wolfram-community/3D-Thinning-algorithm-in-Wolfram-Language-by-Chip-Hurst) | [See Original](https://mathematica.stackexchange.com/questions/186503/3d-thinning-algorithm-in-mathematica) | Download to Desktop via Attachments Below**

 ![chip.jpg](https://global.discourse-cdn.com/wolfram/original/2X/7/7f4f3aea9eda5d3d04611abf3626dcfd458f269e.jpeg)

Here I will explain how to perform 3D Thinning with Wolfram Language. Though I will admit there’s a lot of room for improvement.

In this post 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](https://mathematica.stackexchange.com/questions/67466/while-partitioning-the-elements-in-a-list-using-gatherby-can-i-correspondingly/158637#158637).

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]

```

[![enter image description here](https://i.stack.imgur.com/zctlv.png)](https://i.stack.imgur.com/zctlv.png)

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}

```

[![enter image description here](https://i.stack.imgur.com/QfvrZ.png)](https://i.stack.imgur.com/QfvrZ.png)

```
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}]}

```

[![enter image description here](https://i.stack.imgur.com/Tw0gf.png)](https://i.stack.imgur.com/Tw0gf.png)

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

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

{im3d, medialaxis}

```

[![enter image description here](https://i.stack.imgur.com/OEjn3.png)](https://i.stack.imgur.com/OEjn3.png)

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

```

[![enter image description here](https://i.stack.imgur.com/iXLeF.png)](https://i.stack.imgur.com/iXLeF.png)

* * *

## Addendum

Based on papers of [Amenta](https://www.cs.jhu.edu/~misha/Fall05/Papers/amenta00A.pdf) and [Attali](https://www.sciencedirect.com/science/article/pii/S1077314297905361), a better approximation can be found, but we need to be able to find intersections of spheres efficiently.

# Part II

One approach is to apply a topology preserving morphological thinning method. I’ll implement an algorithm by [Lee](https://www.sciencedirect.com/science/article/pii/S104996528471042X). 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}

```

[![enter image description here](https://i.stack.imgur.com/VL512.png)](https://i.stack.imgur.com/VL512.png)

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

```

[![enter image description here](https://i.stack.imgur.com/siYrc.png)](https://i.stack.imgur.com/siYrc.png)

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

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

```

 ![chip.jpg](https://global.discourse-cdn.com/wolfram/original/2X/7/7f4f3aea9eda5d3d04611abf3626dcfd458f269e.jpeg)

Other images:

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

```

> ```
> {3.88302, Null}
> 
> ```

```
{im, thin}

```

[![enter image description here](https://i.stack.imgur.com/0BZD9.png)](https://i.stack.imgur.com/0BZD9.png)

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

```

[![enter image description here](https://i.stack.imgur.com/DvWlN.png)](https://i.stack.imgur.com/DvWlN.png)

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

```

> ```
> {33.3316, Null}
> 
> ```

```
{im, thin}

```

[![enter image description here](https://i.stack.imgur.com/uUFAP.png)](https://i.stack.imgur.com/uUFAP.png)

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

```

[![enter image description here](https://i.stack.imgur.com/8zfel.png)](https://i.stack.imgur.com/8zfel.png)

[3D-Thinning-algorithm-in-Wolfram-Language-by-Chip-Hurst.nb](https://community.wolfram.com/uploads/short-url/4BCZ1B55E93e8bruSdvmUjVd6dN.nb) (2.74 MB)

---

<div class="post-metadata">

**Author:** ![wolframcommunity](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/wolframcommunity/32/8335_2.png) [@wolframcommunity](https://community.wolfram.com/u/wolframcommunity)\
**Post date:** [November 30, 2018, 12:30pm UTC](https://community.wolfram.com/t/3d-thinning-algorithm-in-wolfram-language/13370/2 "2018-11-30T12:30:20Z")

</div>

![2015-11-02_10-07-11.png](https://global.discourse-cdn.com/wolfram/original/3X/5/0/50e59c9ead62f729e7130d3ea1482555a75602b6.png)

- Congratulations! This post is now a [Staff Pick](http://community.wolfram.com/content?curTag=staff+picks) as distinguished by a badge on [your profile](http://community.wolfram.com/web/ghurst)! Thank you, keep it coming, and consider contributing your work to the [The Notebook Archive](http://notebookarchive.org)!