12
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Szabolcs Nov 22 '18 at 16:21
  • $\begingroup$ Something like this: i.stack.imgur.com/ICbYL.png This is not yet answer-quality, but it is still useful to share. $\endgroup$ – Szabolcs Nov 22 '18 at 16:26
  • $\begingroup$ This is not completely unrelated... mathematica.stackexchange.com/q/9928/1089 $\endgroup$ – chris Nov 27 '18 at 6:07
16
$\begingroup$

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]

enter image description here

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

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

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

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

{im3d, medialaxis}

enter image description here

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

enter image description here


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.

$\endgroup$
6
$\begingroup$

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}

enter image description here

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

enter image description here

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

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

enter image description here

Other images:

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

$\endgroup$
  • $\begingroup$ Great work, thank you +1 $\endgroup$ – UDB Dec 3 '18 at 15:45

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.