Crack Propagation Velocity Determination by High-speed Camera Image Sequence Processing

The determination of crack propagation velocities can provide valuable information for a better understanding of damage processes of concrete. The spatio-temporal analysis of crack patterns developing at a speed of several hundred meters per second is a rather challenging task. In the paper, a photogrammetric procedure for the determination of crack propagation velocities in concrete specimens using high-speed camera image sequences is presented. A cascaded image sequence processing which starts with the computation of displacement vector fields for a dense pattern of points on the specimen’s surface between consecutive time steps of the image sequence chain has been developed. These surface points are triangulated into a mesh, and as representations of cracks, discontinuities in the displacement vector fields are found by a deformation analysis applied to all triangles of the mesh. Connected components of the deformed triangles are computed using region-growing techniques. Then, the crack tips are determined using the principal component analysis. The tips are tracked in the image sequence and the velocities between the time stamps of the images are derived. A major advantage of this method as compared to the established techniques is in the fact that it allows spatio-temporally resolved, full-field measurements rather than point-wise measurements. Furthermore, information on the crack width can be obtained simultaneously. To validate the experimentation, the authors processed image sequences of tests on four compact-tension specimens performed on a split-Hopkinson tension bar. The images were taken by a high-speed camera at a frame rate of 160,000 images per second. By applying the developed image sequence processing procedure to these datasets, crack propagation velocities of about 800 m/s were determined with a precision in the order of 50 m/s.


Introduction
To understand and describe the structural behavior of concrete, in-depth knowledge of the damage mechanisms is required. A crack formed in concrete is the result of a damage process, which starts with the formation and merging of micro-cracks. The evaluation of the resulting crack propagation velocity can provide valuable information on the process and rate of the damage. Such information is crucial to explaining the behavior of concrete at high loading rates. Several researchers have tried to explain the rate-dependent behavior of concrete based on the combined effects of high loading rate The image sequence processing chain, which is used for crack detection and crack measurement in the work presented here, bears the name "cascaded image analysis" in ref. [17]. It has been described in detail by Liebold and Maas [18][19][20] and is briefly summarized here. We limit ourselves to the analysis of specimens with planar surfaces as recorded by monocular image sequences.
As a first step, a grid of points is defined in the camera image of the first epoch under zero load. The grid points are tracked through the consecutive images of the sequence using a sub-pixel accuracy least-squares image matching technique [21] in order to obtain displacement fields for each epoch with regard to the first epoch. Figure 1 shows the displacement field from a bend test on a concrete beam with a thinned-out point set, wherein cracks cause discontinuities in the displacement vector field.
Materials 2020, 13, x FOR PEER REVIEW  3 of 15 concrete beam with a thinned-out point set, wherein cracks cause discontinuities in the displacement vector field. As an alternative to defining a regular grid of points to be tracked through the image sequence, an interest operator may be used to detect surface points with good image contrast, thus guaranteeing precise least-squares matching results. These points, which are not distributed in perfect regularity, are then triangulated into a mesh. The triangles of all images after the zero-load epoch are subsequently analyzed for deformations. Triangles with cracks penetrating them entirely show significant deformations, while the remaining triangles are unchanged. For that purpose, the relative translation vector ⃗ can be computed for each triangle based on the relative movement of its vertices [19,20]. Figure 2a schematically depicts a divided triangle used as the assumed model for the computation of the relative translation vector. The norm of this vector || ⃗ || can be used as a deformation quantity for brittle material. Figure  2b shows a color-coded map of this quantity. The norm is checked for exceeding a threshold value , that is on an order of magnitude of the precision of the displacement field, i.e., approximately a tenth of a pixel. Due to the development of some distributional noise in the deformation field, filter methods such as the bilateral filter technique can be applied [18]. In classifying crack regions, the connected components of triangles where || ⃗ || are determined by applying a region-growing technique including the hysteresis method with a second threshold, as performed by Canny [22]. Figure 3 shows an example of the connected components of deformed triangles from a load test with a concrete beam. As an alternative to defining a regular grid of points to be tracked through the image sequence, an interest operator may be used to detect surface points with good image contrast, thus guaranteeing precise least-squares matching results. These points, which are not distributed in perfect regularity, are then triangulated into a mesh. The triangles of all images after the zero-load epoch are subsequently analyzed for deformations. Triangles with cracks penetrating them entirely show significant deformations, while the remaining triangles are unchanged. For that purpose, the relative translation vector → t rel can be computed for each triangle based on the relative movement of its vertices [19,20]. Figure 2a schematically depicts a divided triangle used as the assumed model for the computation of the relative translation vector.
Materials 2020, 13, x FOR PEER REVIEW  3 of 15 concrete beam with a thinned-out point set, wherein cracks cause discontinuities in the displacement vector field. As an alternative to defining a regular grid of points to be tracked through the image sequence, an interest operator may be used to detect surface points with good image contrast, thus guaranteeing precise least-squares matching results. These points, which are not distributed in perfect regularity, are then triangulated into a mesh. The triangles of all images after the zero-load epoch are subsequently analyzed for deformations. Triangles with cracks penetrating them entirely show significant deformations, while the remaining triangles are unchanged. For that purpose, the relative translation vector ⃗ can be computed for each triangle based on the relative movement of its vertices [19,20]. Figure 2a schematically depicts a divided triangle used as the assumed model for the computation of the relative translation vector. The norm of this vector || ⃗ || can be used as a deformation quantity for brittle material. Figure  2b shows a color-coded map of this quantity. The norm is checked for exceeding a threshold value , that is on an order of magnitude of the precision of the displacement field, i.e., approximately a tenth of a pixel. Due to the development of some distributional noise in the deformation field, filter methods such as the bilateral filter technique can be applied [18]. In classifying crack regions, the connected components of triangles where || ⃗ || are determined by applying a region-growing technique including the hysteresis method with a second threshold, as performed by Canny [22]. Figure 3 shows an example of the connected components of deformed triangles from a load test with a concrete beam. The norm of this vector → t rel can be used as a deformation quantity for brittle material. Figure 2b shows a color-coded map of this quantity. The norm is checked for exceeding a threshold value δ, that is on an order of magnitude of the precision of the displacement field, i.e., approximately a tenth of a pixel. Due to the development of some distributional noise in the deformation field, filter methods such as the bilateral filter technique can be applied [18]. In classifying crack regions, the connected components of triangles where → t rel > δ are determined by applying a region-growing technique including the hysteresis method with a second threshold, as performed by Canny [22]. Figure 3 shows an example of the connected components of deformed triangles from a load test with a concrete beam.
Materials 2020, 13, x FOR PEER REVIEW 4 of 15 Figure 3. Connected components of crack triangles where || ⃗ || according to [19]. Furthermore, it is also possible to derive crack widths and crack opening vectors for deformed crack triangles with the help of the relative translation vector ⃗ as shown by Liebold and Maas [19,20]. Crack widths can be computed as the projection of the relative translation vector ⃗ onto the crack normal ⃗, which has to be determined first; see Figure 4a. For each deformed triangle, the crack normal can be estimated by a line fitted through neighboring deformed triangles; see Figure  4b.
(a) (b) Figure 4. (a) The crack width r is computed by the projection of ⃗ onto the normal ⃗; (b) crack normal estimation according to [19,20].
The result of the image sequence processing chain described is a complete crack pattern for each image of the sequence plus subpixel accuracy information on local crack width. If a crack propagates, the crack pattern on the specimen's surface and the measured crack widths will change. In the next step, these changes are analyzed to quantify the crack tip motion.
It is also possible to use and extend the algorithms presented here to measurements with stereo camera systems in order to analyze non-planar surfaces as shown by Liebold et al. [23,24].

Crack Tip Detection in Multi-temporal Crack Pattern Images
To estimate the propagation velocity of cracks in an image sequence, it is necessary to detect the crack tips in each image. Direct detection of the crack tip in the images is restricted by the subpixel width of a crack near its tip. As a more reliable alternative, the crack tip position and motion can be determined indirectly using the algorithms presented in the previous section. Accordingly, the displacement field and the relative translation vectors are computed for each triangle of the mesh. Crack candidates are found by thresholding where || ⃗ || , and connected components are determined with a region-growing technique using the hysteresis method with a second threshold (δ2 = 0.5 × δ1); see Figure 5.  t rel > δ according to [19]. Furthermore, it is also possible to derive crack widths and crack opening vectors for deformed crack triangles with the help of the relative translation vector → t rel as shown by Liebold and Maas [19,20].
Crack widths can be computed as the projection of the relative translation vector → t rel onto the crack normal → n , which has to be determined first; see Figure 4a. For each deformed triangle, the crack normal can be estimated by a line fitted through neighboring deformed triangles; see Figure 4b.  [19]. Furthermore, it is also possible to derive crack widths and crack opening vectors for deformed crack triangles with the help of the relative translation vector ⃗ as shown by Liebold and Maas [19,20]. Crack widths can be computed as the projection of the relative translation vector ⃗ onto the crack normal ⃗, which has to be determined first; see Figure 4a. For each deformed triangle, the crack normal can be estimated by a line fitted through neighboring deformed triangles; see Figure  4b.   [19,20].
The result of the image sequence processing chain described is a complete crack pattern for each image of the sequence plus subpixel accuracy information on local crack width. If a crack propagates, the crack pattern on the specimen's surface and the measured crack widths will change. In the next step, these changes are analyzed to quantify the crack tip motion.
It is also possible to use and extend the algorithms presented here to measurements with stereo camera systems in order to analyze non-planar surfaces as shown by Liebold et al. [23,24].

Crack Tip Detection in Multi-temporal Crack Pattern Images
To estimate the propagation velocity of cracks in an image sequence, it is necessary to detect the crack tips in each image. Direct detection of the crack tip in the images is restricted by the subpixel width of a crack near its tip. As a more reliable alternative, the crack tip position and motion can be determined indirectly using the algorithms presented in the previous section. Accordingly, the displacement field and the relative translation vectors are computed for each triangle of the mesh. Crack candidates are found by thresholding where || ⃗ || , and connected components are determined with a region-growing technique using the hysteresis method with a second threshold (δ2 = 0.5 × δ1); see Figure 5.  [19,20].
The result of the image sequence processing chain described is a complete crack pattern for each image of the sequence plus subpixel accuracy information on local crack width. If a crack propagates, the crack pattern on the specimen's surface and the measured crack widths will change. In the next step, these changes are analyzed to quantify the crack tip motion.
It is also possible to use and extend the algorithms presented here to measurements with stereo camera systems in order to analyze non-planar surfaces as shown by Liebold et al. [23,24].

Crack Tip Detection in Multi-Temporal Crack Pattern Images
To estimate the propagation velocity of cracks in an image sequence, it is necessary to detect the crack tips in each image. Direct detection of the crack tip in the images is restricted by the subpixel width of a crack near its tip. As a more reliable alternative, the crack tip position and motion can be determined indirectly using the algorithms presented in the previous section. Accordingly, Materials 2020, 13, 4415 5 of 15 the displacement field and the relative translation vectors are computed for each triangle of the mesh. Crack candidates are found by thresholding where → t rel > δ 1 , and connected components are determined with a region-growing technique using the hysteresis method with a second threshold (δ 2 = 0.5 × δ 1 ); see Figure 5. In the following, the crack tip must be found in the connected component of the deformed triangles. Two possible methods are presented here, an approach based on triangle neighborhood analysis and a principal component analysis approach.
The former method considers the neighborhood of each crack triangle. Figure 6 shows four possible cases in the neighborhood analysis of the crack candidates, shown as bold red triangles. In Figure 6 a-c, connected components of non-deformed triangles, i.e., darker blue triangles with edges marked in thicker, colored lines, in the set of the neighbor triangles are determined. The number of these components are counted. Figure 6a shows the standard case with two components, Figure 6b represents a crack junction with more than two connected components around the dark red crack candidate (sets of triangles with yellow or green edge marking). In Figure 6c, only one component of non-deformed neighbor triangles is determined, implying a possible crack tip. Figure 6d shows a crack candidate at the boundary of the mesh, which can be analyzed with the help of the halfedge data structure [25]. A boundary triangle has at least one halfedge with no opposite halfedge partner. Sometimes, the classification of crack tips can fail as shown in Figure 7a. This might also be the case if the thickness of the area of deformed triangles is larger than one triangle, which can be caused by filtering methods, the smearing of || ⃗ ||. Figure 7b depicts the connected component of crack triangles (blue) in an experiment. Mesh boundary triangles are colored yellow and triangles that are identified as crack tips are depicted in green. In the following, the crack tip must be found in the connected component of the deformed triangles. Two possible methods are presented here, an approach based on triangle neighborhood analysis and a principal component analysis approach.
The former method considers the neighborhood of each crack triangle. Figure 6 shows four possible cases in the neighborhood analysis of the crack candidates, shown as bold red triangles. In Figure 6a-c, connected components of non-deformed triangles, i.e., darker blue triangles with edges marked in thicker, colored lines, in the set of the neighbor triangles are determined. The number of these components are counted. Figure 6a shows the standard case with two components, Figure 6b represents a crack junction with more than two connected components around the dark red crack candidate (sets of triangles with yellow or green edge marking). In Figure 6c, only one component of non-deformed neighbor triangles is determined, implying a possible crack tip. Figure 6d shows a crack candidate at the boundary of the mesh, which can be analyzed with the help of the halfedge data structure [25]. A boundary triangle has at least one halfedge with no opposite halfedge partner. In the following, the crack tip must be found in the connected component of the deformed triangles. Two possible methods are presented here, an approach based on triangle neighborhood analysis and a principal component analysis approach.
The former method considers the neighborhood of each crack triangle. Figure 6 shows four possible cases in the neighborhood analysis of the crack candidates, shown as bold red triangles. In Figure 6 a-c, connected components of non-deformed triangles, i.e., darker blue triangles with edges marked in thicker, colored lines, in the set of the neighbor triangles are determined. The number of these components are counted. Figure 6a shows the standard case with two components, Figure 6b represents a crack junction with more than two connected components around the dark red crack candidate (sets of triangles with yellow or green edge marking). In Figure 6c, only one component of non-deformed neighbor triangles is determined, implying a possible crack tip. Figure 6d shows a crack candidate at the boundary of the mesh, which can be analyzed with the help of the halfedge data structure [25]. A boundary triangle has at least one halfedge with no opposite halfedge partner. Sometimes, the classification of crack tips can fail as shown in Figure 7a. This might also be the case if the thickness of the area of deformed triangles is larger than one triangle, which can be caused by filtering methods, the smearing of || ⃗ ||. Figure 7b depicts the connected component of crack triangles (blue) in an experiment. Mesh boundary triangles are colored yellow and triangles that are identified as crack tips are depicted in green.  (1) and (2) as well as Figure 8a.
where ⃗ = center of triangle i of the connected component, ⃗ = mean of the ⃗ values, n = number of triangles of the connected component, Σ = covariance matrix of the ⃗ values. After this, the covariance matrix can be decomposed using an eigenvalue decomposition: where V = rotation matrix with the eigenvectors, Λ = eigenvalue matrix (diagonal). The center points of the triangles of the connected components can then be transformed: The minimum and maximum values in the x and y directions belong to possible crack tip candidates. The dimension with the greater extension is chosen; see Equation (5) Figure 8c illustrates the extension of the transformed coordinates of the triangle centers as a rotated bounding box. If the crack tip candidate is a mesh boundary triangle, the candidate should be rejected.  (1) and (2) as well as Figure 8a.
where → x i = center of triangle i of the connected component,  The latter method in crack tip identification is the principal component analysis (PCA), as shown in Figure 8. The set of center points of the deformed triangles from the connected component are analyzed by computing the mean and the covariance matrix; see Equations (1) and (2) as well as Figure 8a.
where ⃗ = center of triangle i of the connected component, ⃗ = mean of the ⃗ values, n = number of triangles of the connected component, Σ = covariance matrix of the ⃗ values. After this, the covariance matrix can be decomposed using an eigenvalue decomposition: where V = rotation matrix with the eigenvectors, Λ = eigenvalue matrix (diagonal). The center points of the triangles of the connected components can then be transformed: The minimum and maximum values in the x and y directions belong to possible crack tip candidates. The dimension with the greater extension is chosen; see Equation (5) Figure 8c illustrates the extension of the transformed coordinates of the triangle centers as a rotated bounding box. If the crack tip candidate is a mesh boundary triangle, the candidate should be rejected. After this, the covariance matrix can be decomposed using an eigenvalue decomposition: where V = rotation matrix with the eigenvectors, Λ = eigenvalue matrix (diagonal).
The center points of the triangles of the connected components can then be transformed: The minimum and maximum values in the x and y directions belong to possible crack tip candidates. The dimension with the greater extension is chosen; see Equation (5) and Figure 8b. Figure 8c illustrates the extension of the transformed coordinates of the triangle centers as a rotated bounding box. If the crack tip candidate is a mesh boundary triangle, the candidate should be rejected.
Due to possible movements of the entire specimen between the images, the crack tip position of the previous image should be transformed into the current time step. The correction can be done by adding the mean vertex displacement changes to the previous crack tip triangle center.
The PCA method works well for cracks without branches, but may fail in cases of branching. During the experimentation, the authors could observe only one crack without branches. Here, the PCA method produced fewer outliers than the first presented method and consequently was used.

Crack Velocity Determination
Crack velocities can be obtained as derivatives of the crack lengths over time. This is often done by numerical differentiation with, for instance, finite differences, i.e., forward and central differences. As this may lead to high fluctuations due to uncertainties, the authors decided to compute regression lines. The least-squares method offers the advantage of providing standard deviations as internal error measures. The derivatives are computed at each data point, including the four nearest neighbors. The mathematic model for the regression is: where d = crack length, t = time, m = velocity, c = y-intercept. This can also be written in matrix notation: where A = design matrix, L = observation vector, x = vector of unknowns, v = residual vector. For the i th data point, the matrices are: Materials 2020, 13, 4415 8 of 15 The measurements are weighted using binomial weights (1, 4, 6, 4, 1). The parameter vector x can be obtained by minimizing the residuals and solving the normal equations: where P = weight matrix.
Then, the residual vector v can be calculated: The square of the a posteriori standard deviation σ 0 of the unit weight can be computed as follows: where n = number of observations, u = number of unknowns. The covariance matrix Σ xx of the parameters can be computed as follows: In addition to the actual velocity, the standard deviation of the velocity can be derived from Equation (15):

Velocity and Error Estimation of the Lacquer Barriers
The velocity between two lacquer barriers is: where ∆l = length between the conductive lacquer barriers, ∆t = time difference. It is not possible to compute a standard deviation using the least-squares method to evaluate the error of the velocity measurement. However, the error can be estimated using variance propagation. For uncorrelated measurements, the propagated variance is: To estimate the standard deviation σ m , the standard deviations of the measurements σ ∆l and σ ∆t are required. The thickness of a conductive lacquer barrier is about 3 mm and the expected velocity is about 800 m/s such that a crack needs about 4 µs to traverse it. It is not certain at which position the electrical contact is interrupted. Two barriers must be passed and the angle of impact is not necessarily 90 • . Based on these considerations, we assume roughly σ ∆t = 5 µs and σ ∆l = 3 mm. The standard deviation σ m increases with decreasing ∆t.

Material under Investigation
A normal-strength matrix containing cement, fly ash, and fine sand were used to produce the compact-tension specimens. Table 1 provides the mixture compositions. The finely grained matrix is used as the constituent matrix of strain hardening cement-based composites (SHCC), and its behavior under impact tensile loading has been investigated previously. Since the matrix was named M1 in the previous investigations [26,27], the same name is used in this study.

Dynamic Test Setup
The dynamic tests on compact-tension specimens were performed in a split-Hopkinson tension bar (SHTB). Figure 9 provides the schematic view of the setup and geometry of the specimens. The SHTB's input and transmitter bars, both with a diameter of 24 mm, are made of brass with a Young's modulus of 98 GPa. Adapters were used to fix the compact-tension specimens in the SHTB. The loading principle in the SHTB, based on the propagation of an elastic wave in the input bar, was used to achieve a high displacement rate in compact-tension specimens. A detailed description of the SHTB can be found in ref. [28]. Figure 10a shows the specimen in the test setup.

Dynamic Test Setup
The dynamic tests on compact-tension specimens were performed in a split-Hopkinson tension bar (SHTB). Figure 9 provides the schematic view of the setup and geometry of the specimens. The SHTB's input and transmitter bars, both with a diameter of 24 mm, are made of brass with a Young's modulus of 98 GPa. Adapters were used to fix the compact-tension specimens in the SHTB. The loading principle in the SHTB, based on the propagation of an elastic wave in the input bar, was used to achieve a high displacement rate in compact-tension specimens. A detailed description of the SHTB can be found in ref. [28]. Figure 10a shows the specimen in the test setup. A high-speed camera of type Photron Fastcam SA-X2 (produced by Photron Limited, Tokyo, Japan) was used to monitor the front surface of the specimen during the experiment. The frame rate of the camera was set to 160,000 fps at a resolution of 256 × 184 pixels which covered approximately a rectangle of 40 mm × 30 mm on the specimen's surface; see Figure 10b. The recorded frames were subsequently used for performing the image sequence processing chain as outlined above.

Experimental Results
For the detection of the deformed triangles, the algorithm described in Sections 3 and 4 was applied. Figure 11 shows color-coded maps of the relative translation vector lengths for each triangle for different time steps of one of the experiments. Then, the thresholding was done for the norm of the relative translation vector with different threshold values (0.08, 0.10, 0.12 px). The connected components and the crack tips were determined by means of the PCA method. Figure 12 shows an example with a threshold of 0.10 px.  A high-speed camera of type Photron Fastcam SA-X2 (produced by Photron Limited, Tokyo, Japan) was used to monitor the front surface of the specimen during the experiment. The frame rate of the camera was set to 160,000 fps at a resolution of 256 × 184 pixels which covered approximately a rectangle of 40 mm × 30 mm on the specimen's surface; see Figure 10b. The recorded frames were subsequently used for performing the image sequence processing chain as outlined above. Additionally, four barriers made of conductive lacquer were placed on the rear surface of the specimens for validation as presented in Figure 10c. The conductive barriers were each connected a voltage supply to the signal acquisition device. An oscilloscope of type Pico 4824 produced by Pico Technology, St Neots, UK, was used to record the voltage signals at 20 MHz. Upon propagation of the crack through the barriers, voltage signals were reduced to zero so that the time needed for the crack to propagate between barriers can be obtained. In total, four compact-tension specimens made of M1 matrix were tested in the SHTB.

Experimental Results
For the detection of the deformed triangles, the algorithm described in Sections 3 and 4 was applied. Figure 11 shows color-coded maps of the relative translation vector lengths for each triangle for different time steps of one of the experiments. Then, the thresholding was done for the norm of the relative translation vector with different threshold values δ (0.08, 0.10, 0.12 px). The connected components and the crack tips were determined by means of the PCA method. Figure 12 shows an example with a threshold of 0.10 px.

Experimental Results
For the detection of the deformed triangles, the algorithm described in Sections 3 and 4 was applied. Figure 11 shows color-coded maps of the relative translation vector lengths for each triangle for different time steps of one of the experiments. Then, the thresholding was done for the norm of the relative translation vector with different threshold values (0.08, 0.10, 0.12 px). The connected components and the crack tips were determined by means of the PCA method. Figure 12 shows an example with a threshold of 0.10 px.    Crack triangles (blue) with crack tips, i.e., green colored triangles and the yellow-colored mesh border triangles for the experiment presented in Figure 11.
The propagation length between the consecutive images, and accordingly its velocity, were determined automatically using the image processing chain described in the previous sections. The distance from the notch tip, that is the beginning of the crack, to the first crack tip was measured manually at the left border of the image since the notch tip was not included in the measuring field of the optical sensor. This measure is not necessary for the velocity determination but is useful for the comparison to the lacquer barriers in the crack length over the course of time. Figure 13 shows the crack propagation lengths as a function of time for the different thresholds δ for → t rel as well as for the first three conductive lacquer barriers for the four specimens tested. The last barrier was ignored due to crack branching which happened between the third and fourth barriers. In addition, the optical measurement range does not cover the last barrier; see also Figure 10b. As expected, the photogrammetric method with the threshold of δ = 0.08 px is the most sensitive, detecting crack tips earlier in the image sequence than if higher thresholds are used. The magenta-dotted line shows the crack length as a function of time for the conductive lacquer barriers. Comparing the methods, both techniques show monotonically increasing behavior of similar intensity. There is a slight mismatch in the time domain, which might be caused by synchronization errors but is not relevant to velocity determination. Crack triangles (blue) with crack tips, i.e., green colored triangles and the yellow-colored mesh border triangles for the experiment presented in Figure 11.
The propagation length between the consecutive images, and accordingly its velocity, were determined automatically using the image processing chain described in the previous sections. The distance from the notch tip, that is the beginning of the crack, to the first crack tip was measured manually at the left border of the image since the notch tip was not included in the measuring field of the optical sensor. This measure is not necessary for the velocity determination but is useful for the comparison to the lacquer barriers in the crack length over the course of time. Figure 13 shows the crack propagation lengths as a function of time for the different thresholds δ for ||⃗ || as well as for the first three conductive lacquer barriers for the four specimens tested. The last barrier was ignored due to crack branching which happened between the third and fourth barriers. In addition, the optical measurement range does not cover the last barrier; see also Figure 10b. As expected, the photogrammetric method with the threshold of = 0.08 px is the most sensitive, detecting crack tips earlier in the image sequence than if higher thresholds are used. The magenta-dotted line shows the crack length as a function of time for the conductive lacquer barriers. Comparing the methods, both techniques show monotonically increasing behavior of similar intensity. There is a slight mismatch in the time domain, which might be caused by synchronization errors but is not relevant to velocity determination. Next, the velocities were computed as derivatives of the crack lengths with respect to time. Figure 14 shows the results for the four specimens. The 5-point least-squares velocities obtained by the photogrammetric method for the different thresholds and the two velocities calculated between the conductive lacquer barriers are presented. In addition to the velocities, errors bars are plotted, Figure 13. Crack length over time plots using the photogrammetric method with different thresholds for deformed triangles (red, green, blue) and using the conductive lacquer barriers (magenta): (a) 1st sample; (b) 2nd sample; (c) 3rd sample; (d) 4th sample.
Next, the velocities were computed as derivatives of the crack lengths with respect to time. Figure 14 shows the results for the four specimens. The 5-point least-squares velocities obtained by the photogrammetric method for the different thresholds and the two velocities calculated between the conductive lacquer barriers are presented. In addition to the velocities, errors bars are plotted, showing the range of ± one standard deviation. The photogrammetric curves show some dependence on the threshold. This leads to an offset of the determined velocities on the time axis and to a certain fluctuation in the velocities as determined. In fact, as noted above, lower threshold values lead to an earlier detection of the crack tip. However, the resulting shift on the time axis is irrelevant to the velocity calculations. For the conductive lacquer barriers, only two velocity values can be derived. This is not sufficient for a comparison over the entire time span in which the photogrammetric method was used. However, the velocities are on the same order of magnitude. The standard deviations of the velocity values obtained with the photogrammetric method are in a range between 40 and 170 m/s, the standard deviations of the lacquer barrier based velocities reach several hundred m/s. Due to the relatively high standard deviations, the image sequence based instantaneous velocities can also be translated into average velocities using the least-squares method for all data points. Figure 15 shows the corresponding results: The average velocities obtained with the photogrammetric techniques vary between 650 and 900 m/s with standard deviations in the range of 30 to 90 m/s. The velocities between the conductive lacquer barriers are lower, in a range between 490 and 600 m/s, with higher standard deviations of up to 240 m/s. It is obvious that the velocities obtained from the high-speed camera data are consistently higher than the values obtained from the conductive lacquer barriers. However, it should be noted also that the standard deviations of the lacquer barrier measurements are much higher, and thus less reliable as only three data points are used in the computation. The number of conductive lacquer barriers in the measuring area of the camera could be increased in order to have more data points and to compare the methods, but this Due to the relatively high standard deviations, the image sequence based instantaneous velocities can also be translated into average velocities using the least-squares method for all data points. Figure 15 shows the corresponding results: The average velocities obtained with the photogrammetric techniques vary between 650 and 900 m/s with standard deviations in the range of 30 to 90 m/s. The velocities between the conductive lacquer barriers are lower, in a range between 490 and 600 m/s, with higher standard deviations of up to 240 m/s. It is obvious that the velocities obtained from the high-speed camera data are consistently higher than the values obtained from the conductive lacquer barriers. However, it should be noted also that the standard deviations of the lacquer barrier measurements are much higher, and thus less reliable as only three data points are used in the computation. The number of conductive lacquer barriers in the measuring area of the camera could be increased in order to have more data points and to compare the methods, but this may also lead to increasing measurement errors due to the small distances between the barriers. This indicates that the quality of the conductive lacquer barriers measurements cannot be considered as a reference for the results obtained from the photogrammetric image processing chain due to their lower quality.
Materials 2020, 13, x FOR PEER REVIEW 13 of 15 indicates that the quality of the conductive lacquer barriers measurements cannot be considered as a reference for the results obtained from the photogrammetric image processing chain due to their lower quality.

Conclusions
The article at hand presents a novel technique for determining crack propagation velocity on the basis of recorded high-speed camera image sequences. A precise, automatic, photogrammetric image sequence processing chain was developed, implemented, and validated. While using a high-speed camera means a relatively high instrumental effort, it provides the obvious advantage of offering spatio-temporally resolved, non-contact measurements rather than single-point measurements. The technique is able to detect cracks with a width on the order of 0.1 pixel and may thus be applied even in tracking fissures not visible to the human eye. The crack propagation velocities obtained from four laboratory experiments conducted in this study were on the order of 800 m/s.
For the experiments, a camera with a frame rate of 160,000 fps was used. The high frame rate comes with the drawback of an image format limited to 256 × 184 pixels. Future high-speed cameras will offer higher frame rates and larger image formats, which will be beneficial both for the number of cracks to be analyzed simultaneously and for the precision of the velocity values obtained by the method developed.  Figure 15. Least-squares velocities obtained using all data points as average velocities and the corresponding standard deviations (error bars: ±one standard deviation) for all specimens.

Conclusions
The article at hand presents a novel technique for determining crack propagation velocity on the basis of recorded high-speed camera image sequences. A precise, automatic, photogrammetric image sequence processing chain was developed, implemented, and validated. While using a high-speed camera means a relatively high instrumental effort, it provides the obvious advantage of offering spatio-temporally resolved, non-contact measurements rather than single-point measurements. The technique is able to detect cracks with a width on the order of 0.1 pixel and may thus be applied even in tracking fissures not visible to the human eye. The crack propagation velocities obtained from four laboratory experiments conducted in this study were on the order of 800 m/s.
For the experiments, a camera with a frame rate of 160,000 fps was used. The high frame rate comes with the drawback of an image format limited to 256 × 184 pixels. Future high-speed cameras will offer higher frame rates and larger image formats, which will be beneficial both for the number of cracks to be analyzed simultaneously and for the precision of the velocity values obtained by the method developed.