Strategy for Crack Width Measurement of Multiple Crack Patterns in Civil Engineering Material Testing Using a Monocular Image Sequence Analysis

An image sequence analysis procedure is developed to quantitatively analyze complex multiple crack patterns in tension tests of fiber-reinforced composite specimens. Planar textured surfaces of such specimens can be observed with a monocular image sequence using a camera of suitable spatial and temporal resolution. Due to the narrow crack paths, a dense high-precision displacement vector field is computed applying least-squares image matching techniques. Some uniformly distributed matching points are triangulated into a mesh. To measure deformations, principal strains and crack widths are computed for each face. Stretched triangles presumably containing one or multiple cracks are subdivided into three new triangles to densify the mesh in critical regions. The subdivision is repeated for some iterations. The crack width computation of the triangles requires at least three vertices and its displacements. Due to the dense displacement vector field, there are more points available. In this paper, an algorithm for the crack width computation in a least-squares fit is presented.


Review of Related Work
Photogrammetry has a high potential in deformation measurement in civil engineering material testing due to its high accuracy and high resolution. Classical methods such as inductive displacement transducers, inclinometers, or strain gauges offer only point-wise measurements. During the last few years, several contributions about photogrammetry in material testing were published. Early photogrammetric applications in the field of deformation measurement dealt with artificial targets that were tracked in an image sequence of monocular or multi-ocular camera systems (Maas 1998;Fraser and Riedel 2000;Whiteman et al. 2002;Benning et al. 2004;Hampel and Maas 2003; Barazzetti and Scaioni 2010). For example, Barazzetti and Scaioni (2009) presented a 2D image-based method for crack analysis using a digital camera, an orientation frame, and a pair of signalized supports. Other authors used digital image correlation techniques for tracking natural points or artificial random patterns on the surface of the probe (Maas and Hampel 2006;Hampel and Maas 2009;Koschitzki et al. 2011;Maas 2016, 2018). There are further publications considering image processing techniques (edge detection) to extract cracks as the Route-Finder algorithm as well as the Fly-Fisher algorithm (Dare et al. 2002;Niemeier et al. 2008;Detchev et al. 2012). In addition, Detchev et al. (2016) also used a multi-camera system to observe beams in a cyclic load test. There, the amplitudes and offsets of a sinusoidal function for each coordinate as a function of time were determined in a least-squares fit. Concerning crack movements,  pointed out that the crack opening vector has got three components. The first one is the crack width that is normal to the crack, the second one is parallel to the crack course, and the third one is perpendicular to the first two components (perpendicular to the surface).  also refer to Irwin (1958) describing several theoretical modes of fracture. Not all of them can be captured correctly by monocular image observations. In case of movements perpendicular to the surface (out-of-plane movements), 3D systems like stereo cameras have to be used. Görtz (2004) presented an algorithm to compute crack widths in rectangle elements based on the displacement vectors of the vertices. He considered two components: parallel and perpendicular to the crack course. However, global rotations between the reference and the subsequent epochs are not considered.
Furthermore, there exist several companies offering stereo systems and software using digital image correlation (DIC) techniques (e.g., GOM ARAMIS from GOM or VIC-3D from Correlated Solutions, Inc.). These commercial software packages analyze image sequences, compute displacement fields, and visualize principal strains. However, actual metric crack widths have to be measured manually by clicking points to define distances that should be observed. In the work of Liebold and Maas (2018), three approaches for automatic crack width computation in triangle meshes are presented. These methods also consider rotations between the epochs of monocular image sequences and represent an extension to the approach of Görtz (2004). In this publication, the work of Liebold and Maas (2018) is continued.
One of the algorithms is extended and applied to analyze multiple crack patterns rather than single cracks. The extension includes the subdivision of the mesh in critical areas and a least-squares refinement for the crack width computation. Multiple crack structures appear in tension tests of specimens consisting of special composites (here, SHCC: Strain-Hardening Cement-based Composites). This multiple cracking ensures the ductility behaviour that is intended for special applications in civil engineering and the measurement of crack widths is an important issue (Curosu et al. 2017). Figure 1 shows an example of a multiple crack pattern of a deformed SHCC specimen.

Image Analysis-Basic Algorithm
The basic algorithm of crack detection on the basis of analyzing discontinuities in deformation vector fields determined by least-squares matching has been described in detail by Liebold and Maas (2018). An overview is also given here: a tension test is performed on an SHCC specimen with a planar surface. This experiment is observed by a monocular camera system whose optical axis is perpendicular to the surface. During the whole experiment, the relative orientation between camera and surface must not change and deformations must only appear in the plane to be observed. In addition, the surface must show a suitable natural or artificial pattern, such that there is enough texture for the matching process. During the tension test, an image sequence is recorded whose first image under zero load is defined as the reference image. In this image, a set of points is defined. These vertices can be arranged in a regular grid or can be computed by an interest operator; for instance, the Harris operator (Harris and Stephens 1988). In the subsequent images, the displacements of the points are determined with sub-pixel accuracy by least-squares matching (Ackermann 1984;Förstner 1984;Grün 1985Grün , 2012. In each epoch, the vertices are triangulated into a mesh using the Delaunay algorithm, see Fig. 2. Subsequently, the triangles are analyzed for changes. Deformations can be detected by the computation of principal strains (Appendix A) as shown by Liebold and Maas (2016). A visualization of these values is depicted in Fig. 3.
The principal strains only show deformed regions. However, there is no information about the actual crack widths yet.

Limitations of Monocular Observations
Monocular observations of crack patterns require some assumptions. If these assumptions are not fulfilled, systematic errors can influence the measurements: • There are alignment errors if the optical axis is not perfectly perpendicular to the surface. • Moreover, errors due to relative movements between camera and specimen during the experiment are possible. • In addition, there are perspective errors due to out-ofplane movements on the surface of the specimen. • If the measured values are transformed from image to object space, scaling errors can appear.
• There are further errors due to lens distortion which are minimized by camera calibration for instance using the Brown parameters (Brown 1971).
The projective errors of the first and the second point could be minimized using an orientation frame with at least four targets with known coordinates (similar to approach of Barazzetti and Scaioni 2009). The inner geometry of the frame has to stay constant during the load test and the frame has to be attached to the surface, such that it stays undeformed but parallel to it. The measured image coordinates can then be corrected using a projective transformation to the frame system. However, this is not applied in the experiments of this publication.

Differences of Principal Strains and Crack Widths
This subsection explains why we use another model for the deformation measurement (crack width algorithm). Several free and commercial software packages as well as Liebold and Maas (2016) use the computation of principal strains to detect deformed areas as it is shown for triangle meshes in Appendix A. Another way for deformation detection is the determination of crack widths in triangle meshes is presented by Liebold and Maas (2018), a short overview is given in Appendix B. An important difference between both quantities is the underlying model. For principal strains, it is assumed that the surface element is deformed in a nonrigid way. Mathematically, the model is based on an affine transformation (see Appendix A). On the other hand, the crack width computation assumes that the triangle is split into two parts and one of these parts has experienced a relative rigid movement, see also Appendix B, Figs. 4 and 5 for the effect in the mesh. Another difference is that the principal strain is a unit-less quantity, whereas the crack width is a metric quantity. The principal strain can be interpreted as a ratio of lengths, while  Differences of models. a Reference triangle, b model of nonrigid deformation which is assumed for principal strains, and c splitting the triangle into two parts at the crack front and relative rigid movement of the upper part (assumed for crack width computation) the crack width represents the distance of the perpendicular movement between the parts left and right of the crack. In the following, some additional considerations to the work of Liebold and Maas (2018) are presented. First, the effect of the triangle size on principle strains is considered. The principle strains are varying according to the size, because they represent a relative change. The following example illustrates this effect, see Fig. 6. There are two triangles; the bigger one has twice the size of the smaller one. Due to a crack opening with the width r, the triangles are deformed.
We consider the simple case that the heights h i of the triangles are parallel to the crack normal. Then, the principal strain for the smaller triangle can be computed as follows: it can be obtained by the stretch ratio d 1 of the heights of the subsequent and the reference triangle: Or rather, the Cauchy strain e 1 is: The strain d 2 of the second triangle with twice the size of the first triangle ( h 2 = 2 ⋅ h 1 ) and the same crack width r is: (2) e 1 = d 1 − 1 = r h 1 . (3) And the Cauchy strain e 2 is: The Cauchy strain is halved if the triangle edge sizes are doubled in this example: In particular, this effect leads to two problems: first, if there are different triangle sizes in a mesh, then the principal strain values are not comparable. Second, it is difficult to define thresholds for critical strains if the sizes differ. Figure 7 depicts the latter case. In the figure, some matches fail (red points), because the crack runs through the corresponding patches (grey squares). Therefore, some triangles are larger in the crack region. The computation of crack widths avoids these problems, as the crack width is an absolute quantity and not a ratio as the strain. Therefore, crack widths are used for the analysis in this paper. Nevertheless, the calculation of crack widths also has got some limits. An important intermediate result is the relative translation vector. Different relative movements within a triangle mesh could lead to the same relative translation vector of the upper point (single point on one side of the crack) in the triangle, see Fig. 8. The cases (a) and (b) (4) e 2 = d 2 − 1 = r 2 ⋅ h 1 .
(5) e 1 = 2 ⋅ e 2 . illustrate that it is not possible to say where the crack runs through the triangle. It is somewhere between the baseline and the upper point. Moreover, it is not unambiguous which direction the crack has got if only the displacements of the triangle vertices are analyzed; compare (a) and (c). Further information is required. With the help of the crack normal � ⃗ n , the direction of the crack is given and the crack width can be computed by projecting the relative translation vector ⃗t rel onto this normal vector � ⃗ n . Otherwise, only an upper limit for the crack width can be computed with the length of the relative translation vector. Appendix C shows the computation of the crack normals as shown by Liebold and Maas (2018). As one can see from Fig. 8d, it is possible that the upper part has an additional relative rotation. If only the displacements of the triangle vertices are given, the relative rotation cannot be derived. In the model, it is assumed that there is no relative rotation.
Although the whole SHCC specimen has a ductile behaviour, the concrete between cracks can be considered as a stiff material, especially if there are multiple cracks, such that the model with two rigid parts inside the triangle fits better than the affine model using the principle strains. For other, more ductile materials as steel, the affine model should be used, especially if cracks do not appear.

Image Analysis for Multiple Crack Patterns
In this section, a strategy for the deformation analysis of multiple crack patterns is presented. The single steps are explained in the following subsections. Figure 9 depicts the steps of the presented algorithm in a flowchart. The algorithm is designed hierarchically to analyze critical regions by a denser mesh. The workflow begins with the tracking of a dense point grid due to the narrow crack paths. In the second step, a uniformly distributed selection of these points is triangulated into a mesh to avoid small triangles due to inaccuracies in computing crack widths. Next, the mesh is densified in regions where cracks appear. Then, the relative translation vectors are computed including points inside the triangles from the dense point grid. At the end, crack normals and widths are calculated for deformed triangles.

Dense Displacement Field
Due to the narrow crack structures of the fiber-reinforced probes, a dense regular grid of points is used, such that the patches have a big overlap. These points are tracked using least-squares matching (LSM) in the following epochs. If cracks cross matching patches, the assumption of a linear patch deformation in LSM may not be justified, and either the standard deviation of the shifts becomes large or the algorithm does not converge and fails, see Fig. 10a. The first case is typical for thin cracks and large patch sizes. In this approach, it is considered to be failed if the standard deviations of the shifts exceed a threshold. To increase the success rate, the patches of the matching points where LSM fails can be adapted, such that the crack does not cross the patch anymore. LSM is repeated with different predefined patch sizes, but in several cases, LSM fails for all these patch sizes. As shown in Fig. 10b, LSM succeeds for the blue patches now, but it still fails for the red patches. In Fig. 10, possible overlap of the patches is not shown. Successful matching requires suitable texture within the patch in both coordinate directions, and, thus, requires a certain minimum patch size.

Triangle Mesh Creation and Analysis
The crack width analysis for small triangles is inaccurate. In addition, when using larger triangles, areas not affected by cracks can be detected more efficiently. Therefore, a coarse set of points is used for the triangulation. These points can be obtained as follows: all successful matching points from the dense grid are inserted in a regular grid with a defined grid size. Then, the nearest neighbours to the centers of the grid cells are used to define the coarse set of points, see also Fig. 11a. If the cell is empty (no successful matching points inside), it is not considered. After this, the coarse set of points (red encircled points in Fig. 11) is triangulated into a mesh using the Delaunay algorithm, see Fig. 11b. The blue points are the centers of the grid cells. Points where matching fails have a grey colour. The triangle edges of the thinned-out mesh are depicted in green. Figure 12 shows a thinned-out triangle mesh of a regular coarse set of points from an experiment.
Next, for each triangle, the norm of the relative translation vector ||⃗t rel || is computed using the algorithm shown in Appendix B.

Densification of the Mesh
Subsequently, the resolution of the mesh is increased in regions where cracks appear. To achieve the densification, triangles with lengths of the relative translation vectors (Appendix B) larger than a user-defined threshold ( ||⃗t rel || > ) are split into three parts, see Fig. 13. The threshold should be in the same order of magnitude as the precision of least-squares matching [sub-pixel precision, in bad cases 0.1 px (Grün 2012)]. In the experiments presented in this approach, a threshold of 0.075 px is used. The triangle is split as follows: the new vertex for the mesh is in the set of successful matching points inside the triangle (green points in Fig. 13) and it is also the nearest neighbour to the triangle center (red point in Fig. 13). A triangle will not be subdivided if the new triangles would have too small edge lengths. The minimal possible edge length is set to 5 px. This subdivision procedure is repeated for some iterations (here, five iterations).
The effect of densification is visualized in Fig. 14.

Crack Width Computation in a Dense Displacement Field
This section deals with the refinement of the relative translation vectors, which is the central innovation of the algorithm presented in this paper.

Model
First, the displacement vectors of the matching points are computed as follows: To model the crack movement, the displacement vectors are partitioned into two clusters (set M 1 and set M 2 ). Each set belongs to one side of the crack, see Fig. 15. A method for clustering is described in the following subsection. The mathematical description can be done as follows: The points of set M 1 are transformed using a rigid transformation with the parameters ⃗t and due to some movement in the planar space as a consequence of the application of a force to the specimen. For the second set ( M 2 ), a relative translation vector ⃗t rel is added to the transformation due to the crack opening, see Eq. (7). The presented model is a simplification of the model used for extended finite-element method (Moes et al. 1999). For each triangle, there is a separate set of parameters: where � ⃗ p j includes the coordinates of the subsequent epoch, � ⃗ p ref,j includes the coordinates of the reference epoch, ⃗t is the translation vector, is the rotation matrix, and ⃗t rel is the relative translation vector.
The rotation matrix is parameterized with two parameters (c and s, see Eq. 8). They are linear in the observation equations, but necessitate a constraint: where c is cos , first rotation parameter, s is sin , second rotation parameter, and is the rotation angle.
The first step for the computation is the cluster analysis to know which point belongs to which set. The next subsection concentrates on this.

Clustering
For the classification, the parameters from Appendix B ( , ⃗t , ⃗t rel ) derived from the three vertices of the triangle can be used. According to Eq. (7), it is tested for each point inside the triangle if the upper part or the lower part of the equation leads to a smaller deviation i .
Equation 9 shows how to decide whether the point with index k is assigned to the set M 1 or to the set M 2 : Figure 16 shows an example from experimental data. One triangle is shown, and the red and blue points belong to the two sets of the clustering result. The crack running through the triangle is visible and separates the two sets.

Alternative Method for Clustering
There are also alternatives to partition sets. Several methods are known from machine learning. In our case, the unsupervised k-means algorithm (Lloyd 1982) with k = 2 can be used. The input data for k-means are the displacement vectors of the matching points. Due to the relative translation of set M 2 , the displacement vectors should differ from the first set. Initial centroids for the clusters can be obtained with k-means++ (Arthur and Vassilvitskii 2007). Figure 17a depicts the displacement vectors between the epochs if there is only a translation and a relative shift. Figure 17b shows the two clusters of their coordinates. The scattering is only caused by random errors in the matching process.
If there are global rotations, the displacement vectors will show according systematic effects. Figure 18a, b depict an example. The scattering of the vector clusters that is visible is only caused by the global rotation, not by random errors from matching. This could lead to wrong classifications if the relative translations are smaller than the systematic scattering. In most of our experiments, global rotations were very small such that the effect can be ignored in the clustering with k-means.
If there are significant rotations, the k-means method can fail. In that case, the first presented algorithm should be preferred for clustering.

Least-Squares Adjustment
The dense displacement field offers the possibility to consider more than the three triangle vertices for computing the relative translation vector. Therefore, the model is overdetermined and the parameters can be calculated in a leastsquares fit. The observation equations are derived from Eq. 7: where � ⃗ v j is the residual vector of point j. The number of points n is equal to the cardinality (#) of the union of sets M 1 and M 2 and is also the sum of the individual cardinalities of the sets.
The coordinates of the subsequent epoch are considered as observations. The translation vector ⃗t , the rotation matrix of the rigid body transformation and the relative translation vector ⃗t rel are unknowns. Next, the observation equations have to be linearized. Only rotation parameters are linearized applying Taylor's theorem because they are nonlinear in the constraint equation: where c 0 and s 0 are the initial values of c and s, dc and ds are the corrections to c and s.
If there are small global rotations, the initial value of c 0 can be set to 1 and s 0 can be set to 0. The linearized observation equations can be expressed as: where ⃗ l j is the reduced observation for � ⃗ p j , 0 is the initial rotation matrix, is the matrix with corrections to . The parameters are collected in the vector of unknowns : and it can be decomposed in the vector of initial parameters 0 and the vector of corrections to the unknowns .
The initial parameter vector 0 is = t x t y c s t rel,x t rel,y T (15) = 0 + . The vector of corrections to the unknowns is In matrix notation, the linearized observation equations can be written as: where is the reduced observation vector, is the residual vector, is the Jacobian matrix, is the condition matrix, is the vector of inconsistencies.
(17) = t x t y dc ds t rel,x t rel,y T .
(18) + = ⋅ subject to ⋅ = Due to the constraint, the Gauss-Markov model is extended with the method of the Lagrange multipliers .
The solution of this system can be obtained with the extended normal equations: where = 1 2 ⋅ , vector of Lagrangian multipliers.
where n 2 is the number of points in set M 2 .
The upper part of the normal matrix T ⋅ can be computed directly (see Eq. 25) such that it is not necessary to compute the Jacobian matrix and the observation vector in order to be more efficient. The right hand side vector T ⋅ is expressed in Eq. 26. The parameter vector is obtained by solving the extended normal equations (Eq. 24). The translation parameters ⃗t and ⃗t rel can be found directly in the vector. The rotation parameters c and s have to be corrected: The process of the computing of the normal equations with the new -matrix should be repeated until the absolute corrections to the unknowns dc and ds fall below a threshold. Considering the relative translation vector ⃗t rel , it can be

The Jacobian matrix is
where n 1 is the number of points in set M 1 .

The reduced observation vector is
The condition matrix is The vector of inconsistencies can be written as: decomposed in a part perpendicular to the crack and another part parallel to the crack. The crack width r is computed by the scalar projection of the relative translation vector ⃗t rel onto the crack normal � ⃗ n , see Appendix B (crack width computation in triangles), Appendix C (crack normal determination), Fig. 32b and Eq. 28.
The algorithm for the computation of the absolute value of the relative translation vector ||⃗t rel || can also be applied on triangles on 3D surfaces. The 3D coordinates of the vertices of the triangle can be transformed to a local 2D system in the reference and in the subsequent epoch. Then, the 2D algorithm can be applied. Global translation and rotation have to be discarded due to the transformation to 2D.

Experimental Results
In this section, some results of a quasi-static tension test of a SHCC specimen are shown. The loading force is increased stepwise that leads to increasing multiple cracking. The resolution of the observing camera is 5184 × 2912 px. The width of the probe is 4 cm and the length between the clamps is approximately 10 cm. The following parameters are used for the geometric analysis: • The size of the grid cells for the definition of the coarse subset of points is set to 150 px. • The number of iterations for the densification of the mesh is set to 5. • A triangle is subdivided if ||⃗t rel || > = 0.075 px and if the side lengths of the new triangles are greater than or equal to 5 px. • The 1st method of Sect. 2.4.2 is used for clustering in the least-squares method. • If ||⃗t rel || > = 0.15 px, the triangle is considered as crack candidate (Appendices B and C).
The crack widths r (or rather ||⃗t rel || for triangles with ||⃗t rel || ≤ ) can be visualized in a colour-coded map. The triangle crack widths are depicted for four epochs of an experiment in Fig. 19. The increasing multiple cracking is shown and this behaviour is typical for SHCC. It ensures the ductility of the material. Furthermore, it is possible to create 3D visualizations where triangles are transformed to prisms whose heights correspond to the crack widths in addition to the colour code, see Fig. 20.
As already described in Sect. 1.4, principal strains depend on triangle size, whereas the crack width values should not change. In Fig. 21, the influence of different triangle sizes on the principal strains and on the crack widths is shown. In the right images, the edge lengths of the triangles are divided in half, and as expected, the strains appear much greater in the right image. Considering the crack widths, only small differences caused by discretization can be seen. As already mentioned, the presented results in this section show a typical behaviour of a tension test with a SHCC specimen. The comparison to other measuring methods is very difficult: Strain gauges cannot be fixed to the probe, because they would influence the development of the multiple cracking. Two inductive displacement transducers are used to measure the entire extension between the measuring area between the clamps, but they only give single values.

Detailed Test of the Algorithm on the Basis of Synthetic Image Data
The photogrammetric crack pattern analysis procedure as described above delivers results with rather high internal precision figures and rather high spatial resolution. As a consequence, it is almost impossible to provide independent reference measurements which could serve for an external accuracy test. Therefore, we decided to use a synthetic data set with synthetic images containing defined deformations as a basis to test the developed algorithms. Herein, precision, accuracy, and reliability are analyzed. Accuracy describes the deviations between the measurements and the true values including systematic errors, whereas precision shows how measurements differ from each other due to random errors. In this paper, reliability describes the robustness of the algorithm, and therefore, the ratio of outliers is used for the evaluation. The given shifts are compared to the measured vectors in two ways: the first one is the analysis of the 2D relative translation vectors and the second one is the 1D analysis of the computed crack widths. The following parameters are used for the geometric analysis: • The size of the grid cells for the definition of the coarse subset of points is set to 75 px. • The number of iterations for the densification of the mesh is set to 5. • A triangle is subdivided if ||⃗t rel || > = 0.075 px and if the side lengths of the new triangles are greater than or equal to 5 px. • The 1st method of Sect. 2.4.2 is used for clustering in the least-squares method. • If ||⃗t rel || > = 0.15 px, the triangle is considered as crack candidate (Appendices B and C).

Generation of the Images
To get a reference image (undeformed state), a random pattern is generated onto a grey (almost white) background. The resolution of the image is set to 300 × 2000 px. Because of sharp edges, the image is blurred with a Gaussian smoothing filter, see Fig. 22a. For the deformed state, a rotation with an angle of = 2 • (rotation about the image center x c ; y c T ) and a translation of ⃗t = 3 px; −1 px T was simulated for the whole image. For the right part of the image, an additional relative

Measurements
For all generated images, the procedure of Sect. 2 is applied. For each triangle, the relative translation vector and the crack width are computed and the data are analyzed below. The test is done with the three-point algorithm (labeled with 3p, Appendix B) where only the three vertices of the triangle are used and it is done with the least-squares adjustment algorithm (labeled with ls, Sect. 2.4) including also the matching points not belonging to the mesh inside the triangle. Figure 23 shows the result of the algorithm for the crack width computation of one of the cases. The deformed triangles are detected correctly.

Statistics of the Relative Translation Vectors
First, the relative translation vectors are considered. The relative shifts are measured in the coordinate system of the deformed state. To compare the given values used in the (29) For the left side of the crack: x c y c For the right side of the crack: where i is the index of the crack triangle. The mean of the n relative translation vectors for all crack triangles is: where n is the number of crack triangles/observations. The empirical covariance matrix is computed as follows: where s 2 x is the variance of the x values of ⃗t rel , s 2 y is the variance of the y values of ⃗t rel , and s xy is the covariance of x and y of ⃗t rel . To find multivariate outliers, the Mahalanobis distance can be used if the data is normally distributed. The squared Mahalanobis distance MD 2 is distributed according to the 2 distribution: where b is the dimension (here b = 2). According to the 3-rule in the one-dimensional case, the confidence level = 1 − is set to 99.73% . A data vector is considered as outlier if the following condition for the squared Mahalanobis distance is fulfilled: where is the significance level.
The outlier test is done in an iterative process. The mean vector and the empirical covariance matrix are computed in each iteration and only the data vector with the highest Mahalanobis distance is rejected if Eq. (35) is fulfilled. The process is repeated until the maximum of the squared Mahalanobis distances is below the critical value from the 2 distribution. To evaluate the precision, an eigenvalue decomposition of the empirical covariance is conducted: where is the eigenvector matrix; is the eigenvalue matrix (diagonal). s 1 and s 2 are the square roots of the eigenvalues and can be used as a quantity for the precision measurement. Figure 24 shows the scatter plots for the smallest shift distance of v = 0.2 px computed with the three-point algorithm. Table 1 depicts the relative translation vectors in the reference coordinate system and the corresponding vectors in the deformed state. On the left side of Fig. 24, some outliers are visible. In the center plots, the outliers are removed and the confidence ellipses with a confidence level of 95 % (red dashed ellipses) and 99.73 % (red dotted ellipses) are plotted. On the right side, the histograms of the crack widths and the kernel density estimation (magenta dashed line) are depicted. There, the red dotted vertical lines mark the reference values. Figure 25 depicts the plots for the simulated shift v = 0.2 px computed with the least-squares algorithm. In Figs. 26 and 27, there are the plots for the simulated shift of v = 0.6 px.
To evaluate the accuracy, the empirical covariance matrix * is computed using the given expected vector � ⃗ used for the image generation: where s * 2 x is the variance of the x values of ⃗t rel , s * 2 y is the variance of the y values of ⃗t rel , and s * xy is the covariance of x and y of ⃗t rel . The principal standard deviations can be obtained by an eigenvalue decomposition of the empirical covariance matrix * : where * is the eigenvector matrix; * is the eigenvalue matrix (diagonal). s * 1 and s * 2 can be used as a quantity for the accuracy measurement. In addition, the accuracy can also be estimated by the following three quantities: The distance from the mean to the reference vector is: The mean distance to reference vector is: The maximum distance to reference vector is: To conclude the results, the precision and accuracy measurements are plotted in diagrams, see Fig. 28. The plots for The reliability is evaluated with the help of the outlier ratio. The Gaussian noise added to the images leads to noise in the displacements. Sometimes, the wrong base edge is determined by Eq. 51 in Appendix B due to noise effects. In such cases, outliers appear. Figure 29 depicts the outlier ratios of all the experiments. The mean outlier ratio is approximately 2%. The highest outlier ratio appears in the experiment with ⃗t rel = 0.05 px; 1.50 px T . Again, the experiment with ⃗t rel = 0.01 px; 0.20 px T shows a special behaviour, too. It has one of the highest outlier ratios and it has a high number of crack triangles that were not detected (false negatives). In this case, the small relative shift is almost perpendicular to one triangle edge crossed by the crack, such that its length rarely changes (smaller than the precision of LSM). Because of that, there are two edges that have almost constant side lengths. In such cases, it is difficult for the algorithm to decide which edge is the baseline for the crack width computation. If the wrong side is chosen, the computed relative translation vector is an outlier.
In conclusion, the precision and accuracy of the relative translation vector components is below 0.06 px. Under suboptimal conditions, the accuracy may reach values up to 0.2 px. The ratio of outliers as a measure of reliability is on average 2%, but the maximum value is about 7%.

Statistics of the Crack Widths
After computing the crack width by projecting the relative translation vector onto the crack normal (Eq. 28), these values can also be analyzed statistically. Because of the projection, the crack normal � ⃗ n is a further quantity that influences the scattering. It is to be expected that the accuracy of the crack normal � ⃗ n is not that high due to discretization errors (line fit with the center points of the crack triangle and its neighbours, see Appendix C), such that the accuracy of the crack widths should be worse than the accuracy of the relative translation vectors. Precision and accuracy of the relative translation vectors for all reference crack widths for the three-point (3p) and the leastsquares algorithm (ls). Left: precision using principle standard deviations from the empirical covariance matrix; center: accuracy using principle standard deviations from the covariance matrix with the given reference shift as average; right: other quantities for accuracy measurement, distance from the mean to the reference relative shift, mean distance from the single to the reference relative translation vector, and maximum distance to the reference relative shift vector The expected crack widths r are the x components of the relative translation vector ⃗t rel,ref used for image generation: The mean of the measured crack widths r m is: where r i is the crack width of the ith crack triangle; n is the number of crack triangles.
The standard deviation of the crack widths s r is used as a quantity for precision: The standard deviation s * r with the reference crack width r as expected value is used as a quantity for the accuracy: is an exception, because the mean and the expected value differ ( s r,y < s * r,y ). The expected value r,y is zero, because there are only movements along the crack, whereas the mean is greater than zero, because all computed crack widths are positive (absolute values), see also the histograms in the second rows of Figs. 24,25,26,and 27. The crack widths are not normally distributed. For the shifts in y as well as in x and y, the precision and accuracy values are higher in case of higher shifts v. This behaviour is expected due to the fluctuation in the crack normal estimation. Higher relative shifts should lead to higher deviations.
In summary, the precision and accuracy of the crack widths is better than 0.10 px for crack widths below 2 px. The accuracy and the precision depends on the crack width itself. Crack widths up to 5 px can lead to accuracy values up to 0.25 px.

Further Remarks
In real monocular experiments, there are further systematic errors. Some of the points are already listed in Sect. 1.3. In addition to these points, in case of measurements with inductive displacement transducers as comparing method, clamps can cause occlusions.
However, the error analysis in this section only considers the computation of displacements as well as the geometric analysis and is separated from the other effects mentioned in Sect. 1.3 and at the begin of this subsection, although the results may be too optimistic.

Conclusion and Outlook
This paper presents a strategy to detect cracks and compute their widths in multiple crack structures. The method is based on the deformation analysis of triangle meshes. The reliability, the precision, and the accuracy are checked using different simulated images with known deformations. Crack widths can be determined with an accuracy of better than 0.1 px in most cases, with lesser accuracy up to 0.25 px under suboptimal conditions. Further work should concentrate on model extensions using angles obtained by leastsquares matching. Another extension could be the rotation of the upper part of the triangle. In addition, the accuracy can be determined with other reference measurements. A further interesting issue would be the application of the algorithm on triangle meshes of 3D surfaces.
In this appendix, a short overview of the algorithm according to the approach of Liebold and Maas (2018) is given. Between the reference and the subsequent epoch, a rigid movement for the three vertices of the triangle is assumed, and in case of a crack running through the triangle, an additional relative translation is added to one of the vertices. Figure 31 shows an example of the movement of the vertices of a crack triangle.
(48) 2 = ⋅ = ⋅ T = ⋅ T . First, the rigid transformation concerning all vertices has to be separated from the relative translation concerning only one vertex. For that purpose, the minimal absolute difference of the edge lengths between the reference and the subsequent epoch is determined and the corresponding edge is considered as constant base edge. The indices of the vertices corresponding to the edge with the minimal absolute distance change i min and j min are the baseline indices b 1 and b 2 : The remaining index of the triangle belongs to the upper point. The coordinates of the two vertices of the base edges in the reference and subsequent epoch are used to compute the parameters of the rigid body transformation ( ⃗t and ): where ⃗t the translation vector; is the rotation matrix.
The formula is extended with the relative translation vector ⃗t rel for the upper vertex � ⃗ p up : Then, the upper vertex of the reference epoch is transformed in the subsequent epoch using the parameters ⃗t and . According to Eq. (53), the relative translation is computed by the difference of the upper point in the subsequent epoch and transformed reference point: Fig. 32 a Transformed reference and subsequent triangle according to Liebold and Maas (2018); b components of the relative translation vector Fig. 33 Mechanical shear forces; red: extended triangles; blue: triangles with || rel || ≤ ; direction of the relative translation vector rel is not parallel to the crack normal anymore. This figure is according to Liebold and Maas (2018) left side and right side of the crack is composed of a perpendicular part and a parallel part, see Fig. 32b. The presented algorithm is based on the approach of Liebold and Maas (2018). In case of shear effects, the crack normal can be determined as follows: • First of all, the approximate crack width is computed with Eq. (56) for each triangle. If the absolute value of the relative shift is greater than a threshold ( ||⃗t rel || > ), the triangle is considered as a critical candidate. • Then, for each critical candidate, all critical triangles of the second-order neighbourhood (neighbours and neighbours of neighbours) are determined. Triangles are considered as neighbours if they have at least one common vertex. • After this, a line fit of the center points of the critical triangles is computed. • The perpendicular on this line is an approximation for the crack normal � ⃗ n (Fig. 33).
In Fig. 34, the line fit is visualized.