SUB-PIXEL ACCURACY CRACK WIDTH DETERMINATION ON CONCRETE BEAMS IN LOAD TESTS BY TRIANGLE MESH GEOMETRY ANALYSIS

This paper deals with the determination of crack widths of concrete beams during load tests from monocular image sequences. The procedure starts in a reference image of the probe with suitable surface texture under zero load, where a large number of points is defined by an interest operator. Then a triangulated irregular network is established to connect the points. Image sequences are recorded during load tests with the load increasing continuously or stepwise, or at intermittently changing load. The vertices of the triangles are tracked through the consecutive images of the sequence with sub-pixel accuracy by least squares matching. All triangles are then analyzed for changes by principal strain calculation. For each triangle showing significant strain, a crack width is computed by a thorough geometric analysis of the relative movement of the vertices.


INTRODUCTION
For the examination of the behavior of concrete structures, civil engineers conduct load tests on concrete beams.For the understanding of the evolution of cracks during the process, the automatic measurement of quantities such as the number of cracks, crack localization and crack widths is important.Different measuring systems are used, for instance strain gauges, inclinometers, inductive displacement transducers or acoustic emission sensors.In addition to these sensors, photogrammetric methods are applied because they offer their high spatial resolution and a high accuracy.Several publications deal with photogrammetry in civil engineering material testing.(Whiteman et al., 2002) and (Fraser and Riedel, 2000) measured vertical displacements of targets placed on a line on the specimens surface with multi-ocular camera systems.There are also methods that are not based on image comparison to reference images.For instance, (Dare et al., 2002) computed polygons along the crack using the fly-fisher algorithm and the route-finder algorithm.Furthermore, they presented a method for crack width measurement based on the analysis of profiles perpendicular to the polygons.(Hampel and Maas, 2003) and (Benning et al., 2004) used multi-ocular camera systems for displacement measurement in image sequences of planar plates with a grid of targets.The advantage of discrete targets is the high accuracy of the displacement that could be achieved.But due to the distance between the targets, there is a poor crack location resolution.(Maas and Hampel, 2006) and (Hampel and Maas, 2009) used least squares matching (LSM) to determine dense image point shifts and compared it with target grids.Crack widths were estimated by the analysis of profiles in x and y direction of the image coordinate systems.(Benning et al., 2004) and (Lange, 2009) presented an algorithm to compute crack widths for each square of a grid based on the method of (Görtz, 2004), where the direction of the crack was considered.(Koschitzki et al., 2011) computed interest points in a zero load image (reference image) and tracked them with LSM in an monocular image sequence.The points were meshed to a triangulated irregular network (TIN), and the ratio of their areas to the reference were visualized.The surface of the concrete specimen had to be textured for matching.(Barazzetti and Scaioni, 2010) presented three image-based methods for displacement measurement in civil engineering material testing.They applied the Wallis filter on natural texture for contrast enhancement.The FAST interest operator (Rosten et al., 2010) was used to get points that were tracked with LSM and cross correlation techniques.(Detchev et al., 2013) used a multi-camera and projector configuration to measure deformations at loaded concrete beams.(Fedele et al., 2013) and (Fedele et al., 2014) combined digital image correlation with finite element methods to obtain dense displacement fields.
The approach presented in this paper continues the work described in (Liebold and Maas, 2016) and (Koschitzki et al., 2011).A short overview of this approach is given in the following chapter.The sections after this concentrate on crack width determination in triangle meshes.The experimental data of quasi-static load tests was obtained together with the Institute for Concrete Structures at Leibniz Universität Hannover and the Institute for Concrete Structures at Dresden University of Technology.Comparisons between other sensors were done by (Schacht, 2014).

IMAGE ANALYSIS FOR DEFORMATION MEASUREMENT
A load test is conducted on a concrete specimen with a planar surface applying a force F to the beam.The optical axis of the camera should be perpendicular to the side face of the beam to be able to analyze the deformations in 2D, see Figure 1.
In the approach presented here, the area of the surface of the beam was about 2.4 m x 0.3 m.For the photogrammetric setup, a Nikon D300 camera with a focal length of 20 mm was used and the frame rate was set to 0.5 fps.The observed region of interest was about 3600 x 1000 px in the image space or 1.1 m x 0.3 m in the object space and is according to one half of the beam.The distance between camera and object was about 1 m.
Due to the low contrast texture of concrete, the surface to be observed has to be prepared with a suitable artificial pattern to guar-Figure 1. Experimental setup for the load test of a concrete beam antee a reliable image template matching.During the experiment, an image sequence is recorded.The probe will develop a crack pattern with increasing load.The first image of the sequence is taken as reference image.In this image under zero load, points are defined on a regular grid or by applying an interest operator, for instance the Harris operator (Harris and Stephens, 1988).
The points in the subsequent images are tracked by Least Squares Matching ( (Ackermann, 1984); (Förstner, 1984); (Grün, 1985)).Initial values for the shifts can be obtained with normalized crosscorrelation.The points are triangulated into a mesh using the Delaunay algorithm (Figure 2).Each triangle is tested for changes in each time epoch by computing principal strains.

Figure 2. Triangle mesh of interest points
In order to calculate the principal strains for each triangle, an affine transformation is computed with given coordinate pairs (x, ỹ) and (x, y) of the triangle vertices: The deformation tensor F is filled with affine parameters (a12, a13, a22, a23).The translation (a11, a21) is discarded.The matrix F can be decomposed in a symmetric and rotation matrix (polar decomposition): where F = deformation tensor R = rotation matrix (orthogonal) V = left strain tensor (symmetric) The left Cauchy-Green deformation tensor V 2 is the product of the deformation matrix and its transpose: In text next step, an eigenvalue decomposition of left Cauchy-Green deformation tensor is conducted.
where C = eigenvector matrix (orthogonal matrix) Λ = eigenvalue matrix (diagonal matrix) The principal strains are the eigenvalues of V .That is why, the square root of the eigenvalues of V 2 are computed.The corresponding eigenvectors (columns of matrix C) define the directions of the strains.
The square root of the larger eigenvalue is the principal strain s that is used for the next steps, see equation 5.The principal strain is a dimensionless quantity that describes a ratio of a distance to its reference.Values larger than 1 stand for an extension.The rotation matrix is discarded because it has no influence to the strain.
The principal strains of each triangle can be visualized in colorcoded map, see Figure 3.If a crack runs through a triangle, the triangle will become extended and will thus have a larger principal strain.Because of the typical accuracy of the matching algorithm (in extrem cases 0.10 px and even better, see (Grün, 2012)), some noise occurs in the principal strain values.To get rid of these noise effects, a bilateral filter can be applied (Liebold and Maas, 2016).The extended triangles with principal strains larger than a threshold can be merged to a region and labeled, see Figure 4.Each region represents one crack.

CRACK WIDTHS IN TRIANGLE MESHES
The principal strain represents a ratio of the current distance to its reference.Whereas crack widths are absolute (metric) values, the strain depends on the size of a triangle in the images.Hence, it is easier to interpret crack widths.Therefore, metric crack widths shall be derived in a next step.Figure 5 shows a triangle that is separated because a crack runs through it.n describes the direction of the crack normal.
As a first idea, it seems obvious that the crack normal is parallel to the principal strain direction.Based on this assumption, the following two sections show two ways how to compute the crack width in a triangle.Later on, shifts along the crack direction are also taken into consideration that leads to systematic errors in the assumption.

Crack widths with differences of the distances in crack normal direction
If a crack runs through a triangle, there will be one vertex on one side and two vertices on the other side of the crack.For the first presented algorithm, the vertex on the first side of the crack is kept fixed.Because of that, the other two vertices are shifted, see Figure 6.A line is defined by p1 and the crack normal and is intersected with the edge s12.The difference of the distances between the intersection point f and p1 and its corresponding value in the reference state is the crack width, see Figure 7.
First, the edge vectors in the deformed triangle are calculated: It is also done for the reference state: The vector of the crack normal n should be normalized: The following equation for the deformed state can be set up with simple vector algebra, see Figure 7.
where d = distance between p3 and f v = second unknown, factor of the edge s23 Equation 9 can be reordered and can be considered as a linear system: Cramer's rule (Cramer, 1750) is used for the computation of the distance d between p1 and the edge s23 in the deformed state: The factor v for the edge s23 can also be calculated with Cramer's rule: The length factor v for the edge s23 should be constant for the reference and the deformed state: Hence, the intersection point of the crack normal and the edge s23 can be determined as follows: The distance between p 1ref and the edge s 23ref in reference state is computed in the next step: The crack width r corresponds to the difference between the distances in the deformed and the reference state: 3.2 Crack widths with differences of the altitudes of the triangles In the following algorithm, one edge ( s12) of the triangle is considered as constant.Because of that, the third point (p3) is shifted, see Figure 8.As one can see in Figure 9, the crack width can be computed with the help of the heights in the triangles and the angle between the crack normal and altitude: Figure 9. Crack width, the constant edge s 12ref is shifted to edge s12 At first, the perpendicular foot f is calculated with the projection of edge s13 on egde s12: The same is also be done for the reference state: The altitude vector in the deformed ( h) and the reference ( h ref ) triangle is the difference between p3 and the foot point f : In the next step, the difference of the lengths of the altitude vectors is computed: The cosine of the angle α between the crack normal and altitude can be determined with the dot product of the altitude vector and the crack normal in the deformed state: The crack width r is then computed with the trigonometric function:

Crack width considering shift effects along the crack
Additional considerations are required if there are mechanical shear forces along the crack (Figure 10).The shifts in this direction cause a change in the principal strain directions obtained from the triangles, in a way, that the principal strain direction is not parallel to the crack normal anymore.Therefore, the crack normal has to be obtained in another way.
Figure 10.Triangle mesh, red: extended triangles; blue: triangles without strains; shear forces affect along crack direction in opposite direction; the principal strain vector s is not parallel to the crack normal n anymore

Determination of the crack normal:
The crack normal in the deformed state under the presence of shear effects can be obtained as follows: At first, only triangles with principal strains larger than a threshold are considered as crack candidates.For each triangle in this crack region, the second order neighbor triangles are determined.Triangles are considered as neighbors if they share at least one vertex.Second order neighbors are neighbors of neighbors.If the triangles of the second neighborhood also belong to the crack region, their geometric centers are used for a weighted Principal Component Analysis (PCA, (Pearson, 1901)).The direction of the second principal component corresponds to the crack normal direction.This method can be seen as a local line fit, see also Figure 11.Each neighbor triangle center j, that is used for the PCA of triangle i, gets a Gaussian weight ω ij that penalizes greater distances: where mi = geometric center of triangle i mj = geometric center of a neighborhood triangle j σ = distance where the weight gets 1 e e = Euler's number The weights ωij are normed for each neighbor triangle j: In the next step, the weighted barycenter b is computed: Then, the weighted covariance matrix Z can be computed: Next, an eigenvalue decomposition is applied on the covariance matrix Z: where C = eigenvector matrix (orthogonal matrix) Λ = eigenvalue matrix (diagonal matrix) The second principal component is the column of eigenvector matrix C that corresponds to the smallest eigenvalue.It is regarded as the direction of the crack normal n.

Crack widths and translation:
The two methods for the crack width determination shown in 3.1 and 3.2 will not work correctly if there are significant shifts along the crack direction.Shifts along the crack direction can be incorporated by using vector algebra, see Figure 12 and Figure 13.First, the coordinates of the triangle in the reference state have to be transformed in a way, that the base edges of the reference and the deformed state have the same orientation.Therefore, p 1ref should have the same coordinates as p1.Because of that, p 1ref has to be transformed: The edge s 12ref should have the same direction as s12: The x-coordinate dx 13ref,t of the edge s 12ref,t can be determined by projecting the edge s 13ref onto the edge s 12ref : The x direction ex,t is set parallel to s 12ref,t : The y component is computed with the help of the Pythagorean theorem: The y direction ey,t is set perpendicular to s 12ref,t : The transformed edge s 13ref,t is composed of its x and y components: The third vertex p 3ref,t is computed by adding the edge vector s 13ref,t to the first vertex: As one can derive from Figure 13, the crack width r is calculated as follows: Optionally, the translation t along the crack can be determined by vector subtraction:

EXPERIMENTAL RESULTS
First, the algorithms are numbered due to an easier description: number of algorithm corresponding section 1 3.1 2 3.2 3 3.3 Table 1.Numbering of algorithms For each extended triangle, a crack width can be computed and displayed in a color-coded map, see Figure 15.The visualizations of the three algorithms look very similar.In the center of the right crack, some differences of the crack widths can be seen.Applying algorithm 3, translations in crack direction could also be computed and the length of the vector could be visualized.
There exist some shifts of up to 1.5 px in crack direction in the center of the right crack.
One could expect that algorithm 1 and algorithm 2 have identical results.To prove this, the residuals between them are computed (equation 38) and analyzed, see Table 2 and Figure 14.
Considering the average and median, the values are nearly the same because the values are almost zero.But looking at the histogram in Figure 14, there is a small tendency to the negative side.r algorithm1 seems to be a little bit larger due to the onesided histogram, the negative median and average.The RM S12 is 0.11 px. 12 = r algorithm2 − r algorithm1 To verify the results, some points were set manually next to the crack in a way that the line between these points is perpendicular to the crack direction.In Figure 16, these points are white.The colored points between the white were set on the crack in order to find the corresponding triangle.
This can only be done in the deformed image where the crack is visible.The shift of the points to the reference image is determined with LSM.With the help of this, the crack width can be determined as the difference of the distances between these points in the deformed and the reference state (equation 39).The corresponding triangle can be found by a nearest neighbor query of the center points of the triangles in the deformed state.
where r manual = crack width of the manual method d ref erence = distance between the points (reference) d def ormed = distance between the points (deformed) As this procedure comes with its own error budget, it rather has the character of a plausibility check.A procedure for providing rigorous reference measurements to verify the accuracy potential of the crack width determination procedure yet has to be developed.
Table 3 shows some statistical analysis concerning the residuals.
The RMS is about 0.3 px.The corresponding histograms are plotted in Figure 17.For algorithm 3, two residuals are computed.4. The residuals of the three algorithms are very similar.Overall, the obtained precision is a about five times worse than the LSM accuracy (0.05 px).This can partly be explained by variance propagation and some uncertainties of the normals concerning algorithm 3.An independent measuring system with a higher precision would be required for a better analysis of the accuracy. RMS

CONCLUSION AND OUTLOOK
An algorithm for crack detection in image sequences was introduced in this paper.The method is based on the analysis of triangle changes with principal strains.Different algorithms for crack width computation based on a thorough geometric analysis of deformed triangles were presented.First results are visualized in color-coded maps.The crack widths can be determined with subpixel accuracy.Future work should concentrate on further tests and evaluations including other independent measuring methods with a higher precision.In addition, the method could extended towards handling 3D meshes.
1) where x, ỹ = coordinates of the subsequent epoch x, y = coordinates of the reference epoch aij = affine parameters

Figure 3 .
Figure 3. Color-coded visualization of the principal strains on the right half of a concrete beam under load

Figure 4 .
Figure 4. Labeled regions with triangles with principal strains being larger than a threshold are merged together.

Figure 5 .
Figure 5. Separated triangle caused by a crack

Figure 6 .
Figure 6.(a): crack through the reference triangle; (b): new triangle in the deformed state; (c): both triangles, vertex p 1ref in reference and deformed state have the same position

Figure 7 .
Figure 7. Reference and deformed triangle, vertex p 1ref in reference and deformed state have the same position

Figure 8 .
Figure 8. (a): crack through the reference triangle; (b): new triangle in the deformed state; (c): both triangles, where the edges s12 and s 12ref in reference and deformed state are identical

Figure 12 .
Figure 12. (a): crack through the reference triangle with a shift in crack direction; (b): new triangle in the deformed state; (c): both triangles, where the edges s12 and s 12ref in reference and deformed state are identical Figure 14.Histograms of the residuals between algorithm 1 and algorithm 2

Figure 16 .
Figure 15.(a): crack width visualization with algorithm 1; (b): crack width visualization with algorithm 2; (c): crack width visualization with algorithm 3; (d): visualization of the absolute values of translations of algorithm 3

Figure 17 .
Figure 17.Histograms of the residuals to the manual crack width

Table 3 .
Statistical values of the residuals to the manual measurements in px

Table 4 .
Statistical values of the residuals to the manual measurements in px without outliers