Crack modeling via minimum-weight surfaces in 3d Voronoi diagrams

As the number one building material, concrete is of fundamental importance in civil engineering. Understanding its failure mechanisms is essential for designing sustainable buildings and infrastructure. Micro-computed tomography (μCT) is a well-established tool for virtually assessing crack initiation and propagation in concrete. The reconstructed 3d images can be examined via techniques from the fields of classical image processing and machine learning. Ground truths are a prerequisite for an objective evaluation of crack segmentation methods. Furthermore, they are necessary for training machine learning models. However, manual annotation of large 3d concrete images is not feasible. To tackle the problem of data scarcity, the image pairs of cracked concrete and corresponding ground truth can be synthesized. In this work we propose a novel approach to stochastically model crack structures via Voronoi diagrams. The method is based on minimum-weight surfaces, an extension of shortest paths to 3d. Within a dedicated image processing pipeline, the surfaces are then discretized and embedded into real μCT images of concrete. The method is flexible and fast, such that a variety of different crack structures can be generated in a short amount of time.


Introduction
The shortest path problem (SPP) is formulated on a planar (directed) graph consisting of vertices and weighted arcs.Given a start and an end vertex, one is interested in finding a path of minimal weight connecting these vertices.Shortest paths represent a useful tool in 2d image processing in several contexts including image quilting [19] and image segmentation.In particular, shortest paths have been used to segment linear, 1d structures such as cracks in 2d images [2].Minimum-weight surfaces were introduced in [21] and [15] as an extension of shortest paths to 3d.The minimum-weight surface problem (MSP) is formulated on a cellular complex consisting of vertices, arcs, weighted facets and cells.Here, the goal is to find a set of facets of minimal weight that is bounded by an input cycle.The voxel lattice in a 3d image can be interpreted as a cellular complex.In this context, minimum-weight surfaces have been used to segment flat structures in medical 3d images [15].The segmentation of cracks in images of concrete is a broad field of research.A large variety of segmentation methods have been studied for 2d and 3d images.Comprehensive overviews and comparison studies in 2d and 3d can be found in [9] and [4], respectively.In image segmentation, ground truths are necessary for an objective output evaluation.Furthermore, they are prerequisite for training machine learning models such as convolutional neural networks [8] and random forests [7].In 2d, annotated data is abundantly available, for example SDNET2018 [11].For 3d images, manual annotation is not feasible due to the large amount of data.Therefore, 3d ground truth images are scarce and one has to rely on artificial data.Previously, the generation of synthetic cracks has been realized via fractional Brownian surfaces [1,4].However, the Hurst index -a measure of roughness -is the only parameter that can be used to control the shape of the surface.Therefore, this model is rather limited.Voronoi diagrams have been used previously for simulating crack propagation via finite element methods [13,16].In this work, we propose a novel method to synthesize crack structures in 3d images using Voronoi diagrams generated from random point processes.The method includes two aspects: First, we use minimum-weight surfaces as a tool to model the macrostructure of cracks.Bounded 3d Voronoi diagrams serve as the underlying cell complexes.This leaves several degrees of freedom such as the choice of the generating point processes, bounding cycles and the facet weights.Second, the computed surfaces are discretized to 3d binary images.The rough microstructure of cracks is modeled via a second Voronoi diagram on a finer scale.Then, the cracks are embedded into patches of real computed tomography (CT) images of concrete.This paper is structured as follows.In Section 2, we outline the concept of shortest paths.Their extension to minimum-weight surfaces is described in Section 3. Section 4 comprises our crack modeling pipeline using Voronoi diagrams.It includes a description of our approach for macro-and microstructure modeling.Section 5 serves as conclusion and outlook to possible future research.

Shortest paths
Let G = (V, A) be a directed graph consisting of a set of vertices V and directed arcs A ⊆ V × V .We use the notation α(a) for the start vertex and ω(a) for the end vertex of an arc a ∈ A. Furthermore, let c : A → R >0 be a function assigning a non-negative weight to every arc in A. A path in G is a finite sequence of vertices and arcs, It is called a cycle (or closed contour) if no arc and no vertex is included more than once except for v 0 = v k .The weight of a path P is given as c(P ) = k−1 i=0 c(a i ).Given two vertices s, t ∈ V , the SPP is looking for a path of minimum weight from s to t.Many algorithms for solving the SPP exist, for example Dijkstra's [10] or Bellman's and Ford's algorithm [5].SPPs can also be formulated as binary integer programs [17,22].Note that this approach is less efficient than the ones described in [10] or [5].However, it gives a good intuition for an analogue approach to compute minimum-weight surfaces which we describe in Section 3. Let s be the start-and t the end vertex of the path and let x be a vector of binary variables x i ∈ {0, 1} assigned to each arc a i .The SPP can then be formulated as subject to Bx = p (2) The vertex-arc incidence matrix B and the vector p in constraints (2) are given as The constraints (2) are flow-conservation constraints.They ensure that every vertex that is part of the path is incident to exactly one incoming and one outgoing arc, except for the start-and end vertex.This results in the fact that any feasible solution must contain a path from s to t and, by the assumption that all costs are strictly positive, this ensures that the optimal solution is in fact a path without repeated vertices.The variables x i are binary by constraint ( 3) and indicate whether arc a i belongs to the path (x i = 1) or not (x i = 0).Finally, the objective in ( 1) is to minimize the weight over all paths from s to t.Note that, in case of negative arc costs, the problem of finding a shortest path without repeated vertices becomes NP-hard (which can be seen by a reduction from the Hamiltonian Path Problem [12]).Thus, additional constraints become necessary and, hence, we assume that all arc costs are strictly positive.

Minimum-weight surfaces
Minimum-weight surfaces have been presented in [21] and [15] as an extension of shortest paths to 3d.For our purposes, let K = (V, A, F, C) be a cellular complex consisting of a set of vertices Further, let w : F → R >0 be a function assigning a non-negative weight to every facet in F .Note that (V, A) defines a (directed) graph.Given a cycle H on (V, A), the MSP is looking for a connected set of facets in K of mimimum weight that is bounded by H.
To this end, arc directions may be assigned arbitrarily.Every facet is considered twice, once in clockwise and once in counterclockwise orientation.The orientation of H must be chosen to be either clockwise or counterclockwise.If the direction of arc a coincides with the direction of its counterpart in an incident facet f (or cycle H), we call a and f (or a and H) coherent.If it does not coincide, we call them anti-coherent.
The MSP can be formulated as a binary integer program analogously to the one in Section 2. It is given as subject to Dy = q (5) The arc-facet incidence matrix D and the vector q in constraints ( 5) are given as 1 if a j and f i are incident and coherent, −1 if a j and f i are incident and anti-coherent, 0 else and Similarly to the flow-conservation constraints in (2), constraints (5) ensure that we obtain a connected set of facets that is bounded by H.The variables y i are binary by constraint ( 6) and indicate whether facet f i belongs to the surface (y i = 1) or not (y i = 0).Finally, the objective in ( 4) is to minimize the weight over all possible surfaces that are bounded by H.
We always assume positive costs to ensure that the output is indeed a connected surface.

Crack modeling
Minimum-weight surfaces have been used for 3d image segmentation with the 3d grid being the underlying cell complex [15].Reversing this approach, minimum-weight surfaces can also be used to model surface-like, connected structures such as cracks on a macroscopic level.
Here, we focus on minimum-weight surfaces in bounded 3d Voronoi diagrams.Our goal is to develop an approach for modeling 2d crack structures to generate semi-synthetic 3d images of cracked concrete.Using this approach in combination with various point process models allows us to control the geometry of the resulting crack structure.
The most striking geometric characteristics of cracks in 3d CT images of concrete have been identified and discussed in [18].In particular, we observe that 1. crack widths are varying and cracks may appear on multiple scales, 2. crack surfaces are not totally smooth but rather rough due to the granularity of the concrete's cement matrix, 3. when propagating through concrete, cracks may branch.An example is given in Fig. 1.These observations are used in the discretization procedure for modeling the cracks' microstructure.
The C i are called the cells of the Voronoi diagram.Note that, if S is not bounded, the Voronoi diagram contains cells of infinite size.In practice, it is often convenient to only consider bounded cells.Therefore, we restrict our attention to the bounded Voronoi diagram given by W ∩ Q = {C 1 ∩ Q, . . ., C n ∩ Q} for some bounded region Q ⊂ S. Note that this operation yields additional vertices, arcs and facets on the boundary of Q that belong to the bounded Voronoi diagram.

Minimum-weight Voronoi surfaces
Bounded Voronoi diagrams in 3d can be considered as a cellular complex.Therefore, given a cycle on the arcs of the cell complex induced by a 3d Voronoi diagram and weighted facets, we are able to compute a minimum-weight surface by solving the optimization problem given in Section 3. As a result, we obtain a connected set of facets that we call a minimum-weight Voronoi surface.

Crack generation
We propose the following method to simulate crack structures via minimum-weight Voronoi surfaces.
2. Compute a random point pattern R ⊂ Q as a realization of some point process model.

3.
Compute the Voronoi diagram generated by R, bounded by Q.
4. Define functions c, w assigning a non-negative weight to each of the arcs and facets, respectively.
5. Choose a vertex on each of the four vertical edges of Q. Denote them by u 1 , u 2 , u 3 , u 4 .Compute shortest paths from u 1 to u 2 , u 2 to u 3 , u 3 to u 4 and u 4 to u 1 , via Dijkstra's algorithm, only using arcs that lie on the boundary of Q. Denote the paths by P 1 , P 2 , P 3 , P 4 .Then, H = ∪ 4 i=1 P i is a cycle on the boundary of Q.The approach above leaves us several degrees of freedom.Facet shape and variability can be controlled by choice of the underlying point process model, while the intensity of the generator process influences the mean facet size.Additionally, the size of the bounding cuboid, the input cycle and the weighting functions c and w can be varied.Minimum-weight surfaces for the same realization of a Poisson point process but different choices of input cycles are given in Fig. 3.Note that the input cycle is not restricted to lie on the boundary of a cube but may be chosen arbitrarily.Moreover, Fig. 4 shows minimum-weight surfaces in Voronoi diagrams generated by Poisson point processes, Matérn cluster processes and regular processes obtained by a force-biased sphere packing [3,6] with a volume fraction of 60%.Arcs and facets were weighted by their lengths and areas, respectively.For the cluster process, we observe a bimodal distribution of the facet areas.The surfaces resulting from the regular model

Discretization
The method described in Section 4.3 outputs a set of vertices of convex facets.In this section, we describe a method to transfer this representation to a discrete image.Let I denote a 3d label image and J a 3d binary image, both of size The discretization procedure is given as follows.
1. Compute a minimum-weight Voronoi surface in a cuboid of size 2. Discretize the Voronoi diagram.For every voxel (p, q, r) do: Set I(p, q, r) = l if (p, q, r) is contained in cell l.
3. Discretize the minimum-weight surface.For every two neighboring voxels (p, q, r), (p , q , r ) (with respect to the 26-neighborhood) do: Set J(p, q, r) = 1 if I(p, q, r) = j and I(p , q , r ) = k for generators j, k whose cells share a facet that is part of the minimum-weight surface.Output J.
The procedure is visualized in Fig. 5.It yields a binary image whose foreground is a piecewise planar structure of constant width.However, real crack structures usually are far more complex as has been pointed out in Section 4. In the following, we propose techniques to account for these observations.

Adaptive dilation
Our goal is to model cracks of varying thickness.We propose a procedure to dilate the foreground of image J: Every x-slice of J is dilated separately and iteratively.We choose a quadratic structuring element of size 2 × 2. The number of iterations depends on a random walk with Bernoulli-distributed increments and index set {0, 1, . . ., d 1 }.The increments are either 1 with probability p or 0 with probability 1 − p.Thus, the crack thickness can be controlled via parameter p.The procedure is visualized in Fig. 6 for different choices of p.Note that the random walk can be substituted by any suitable stochastic process, for example to produce decreasing crack widths.

Microstructure modeling
In order to model the rough microstructure on the boundary of cracks, we compute a second (Poisson-) Voronoi diagram with a higher intensity than the one used for the computation of the minimum-weight surface.Then, for every foreground voxel in the dilated crack image J, we identify the Voronoi cell it is contained in.The whole cell is then discretized with voxel value 1 according to the approach described in Section 4.4.The procedure is visualized in Fig. 7. Afterwards, we apply a median filter to the resulting image.

Crack branching
Crack branches emerge when a crack splits into two or more cracks.Often, the thickness of these branches lies in a range of 1-2 voxels.Branching cracks can be modeled by combining two minimum-weight surfaces obtained from different cycles on the cuboid.If the underlying set of generators for the Voronoi diagram is identical, the surfaces may share several facets, see Fig. 8.

Crack embedding
An approach for embedding crack structures in real 3d concrete images has been proposed in [4].We extract image patches of the same size as the ground truth images from the real CT images.These patches are multiplied voxelwise with the inverse ground truth images.This leads to crack voxels having grayvalue 0 while the background does not change.Cracks and air pores both consist of air.Therefore, they should possess the same grayvalue distribution.We assume these grayvalues to be i.i.d.normally distributed.Mean and standard deviation are estimated via sample mean and sample standard deviation of the empirical distribution of air pore grayvalues.Then the crack voxels are simulated according to that distribution.To smooth the transition between background and crack, we apply a Gaussian filter to crack voxels and all voxels in their 26-neighborhood.The final image together with its ground truth is given in Fig. 8.A corresponding 3d rendering is given in Fig. 9.

Conclusion
In this work, we have presented a novel method to generate artificial crack images.It includes the generation of a macrostructure via minimum-weight surfaces and a discretization procedure for generating its microstructure.
The shape and size of the output can be controlled by several parameters.Thus it allows for the generation of a wide range of surface structures.Our next steps include the generation of a full semi-synthetic data set for training machine learning models and evaluating segmentation methods on multiscale crack images.Furthermore, the model may be extended to include the grayvalue information of real concrete images.We can assume that cracks, when propagating through concrete, take the path of least resistance.Certain parts of the concrete mixture are less prone to cracking than other parts.In particular, this holds true for parts with a higher density.Therefore, facet weights may be derived from the mean voxel grayvalue in their vicinity.

Figure 1 :
Figure 1: 2d slice of a 3d CT image of cracked concrete of size 1000 × 550 × 880 voxels with a voxel edge length of 22.7µm.The image stems from a normal-strength concrete sample that was exposed to a tensile test.Sample: Department of Civil Engineering, University of Kaiserslautern, Imaging: Fraunhofer ITWM, Kaiserslautern

6 .
Compute a minimum-weight surface bounded by P by solving the integer program from Section 3. The concept is illustrated in Fig. 2. It shows a minimum-weight surface in a Voronoi diagram bounded by Q = [0, 1]×[0, 1]×[0, 1].Its generators are a realization of a Poisson point process in Q of intensity 500.Arcs and facets are both assigned unit weights.

Figure 2 :
Figure 2: Minimum-weight surface generation.Left: Poisson-Voronoi diagram bounded by a cube, middle: cycle on the boundary of the Voronoi diagram, right: minimum-weight surface.

Figure 3 :
Figure 3: Examples of minimum-weight surfaces in the same Voronoi diagram but with different choices of input cycles.

Figure 5 :
Figure 5: Slice view visualizing the discretization procedure.Left: Label image obtained from discretizing the Voronoi diagram.Red voxels indicate facets that are part of the minimum-weight surface.Right: Output binary image.

Figure 6 :
Figure 6: Adaptive dilation applied to the input image with p = 0.02 (left) and p = 0.05 (right).

Figure 7 :
Figure 7: Microstructure generation.Left: Discretized Poisson-Voronoi diagram and input crack (red).Right: Set of Voronoi cells (white) that are intersected by the input crack.

Figure 8 :
Figure 8: Crack embedding.Left: Ground truth image after applying a median filter, right: synthesized crack image.

Figure 9 :
Figure 9: 3d volume rendering of the simulated crack structure.