Message Boards Message Boards

GROUPS:

3D Thinning algorithm in Wolfram Language

Posted 14 days ago
202 Views
|
1 Reply
|
4 Total Likes
|

Open in Cloud | See Original | Download to Desktop via Attachments Below

enter image description here

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.

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.

Part II

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

Attachments:

enter image description here - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract