# Image-Based Framework for Concrete Surface Crack Monitoring and Quantification

Authors: ZhiQiang Chen; Tara C. Hutchinson
Journal: Advances in Civil Engineering
 (2010)

---

## Abstract

In the engineering community, nondestructive imaging has been widely used for damage identification by capturing anomalies on the surface or inside of structural elements. In this paper,
we focus on one of the most common damage types observed in civil engineering, namely,
concrete surface cracks. To identify this type of damage, we propose an image-based framework,
whereby optical cameras provide the source images. The framework involves several
advanced image processing methods, including: (i) the determination of damage occurrence
using time-series images, (ii) the localization of damage at each image frame, and (iii) the
geometric quantification of damage. Challenges that may arise when images are obtained in the
laboratory or field environment are addressed. Two application examples are provided to
demonstrate the use and effectiveness of the proposed approach.


---

## Body

## 1. Introduction

In the communities of structural health monitoring (SHM) and nondestructive evaluation (NDE), development of an automated structural damage identification solution has been a key objective. Most SHM methodologies rely on sensing of one-dimensional vibration signals, which are used to extract global modal features as signatures of structural integrity by using system identification-based or statistical pattern recognition-based methods [1, 2]; therefore, the awareness of structural damage is usually based on classifying these indirect signatures. Different from these SHM methodologies, NDE methods, especially when 2- or higher-dimensional imaging methods are employed, are able to provide a direct characterization of local structural damage. Depending on the nature of images, spatial extent and spectral variation of damage usually manifest themselves in the captured images. Subsequently, one can employ photogrammetric or image analysis methods to extract these damage characteristics quantitatively.A variety of imaging technologies have been developed in the NDE community in an effort to detect local structural damage to civil/mechanical systems or components. Widely used imaging devices include infrared thermography, microwave imaging, acoustic imaging, X-ray imaging, and other radiography-based methods [3]. Common types of structural damage encountered in practice include external or internal cracks, voids, delamination, and ablation, to name a few. To detect these different types of damage, the imaging device must often be customized to capture the characteristics associated with the physical damage; hence the resulting equipment may be expensive. In civil engineering, however, optical imaging by means of commercially inexpensive cameras has become an important, economically feasible utility for recording structural damage. Its use has been extensive in the laboratory setting; however, use in the field has also recently become commonplace due to the low cost of high-resolution optical cameras and associated support equipment. To expedite and automate damage recognition using the captured digital image, computational image processing tools have been attempted to detect or track structural damage in static images or dynamic image frames [4–8].In this paper, we consider one of the most commonly observed types of damage in civil engineering, namely, concrete surface cracks. Concrete is the most widely used man-made material for constructed systems, and surface cracks are a reality due to its low tensile strength and environmental sensitivity. As a result, it is well known that a primary damage modality of concrete structures is that of surface cracks. In engineering practice, the identification of concrete surface cracks is an indispensable task in the process of conducting damage- or failure- prevention structural upgrading. Very narrow to hairline cracks, however, as may occur due to the natural process of shrinkage for example, are usually not significant enough to change the global structural integrity of the element hence they are not identifiable using SHM technologies. However, the appearance of these cracks and their propagation in concrete structures may indicate potential inception of degradation of structural integrity or associated softening of the system. In practice, visual inspection of the concrete surface is a common approach. For example, a considerable amount of human effort is devoted to periodic visual inspection of concrete bridges in the United States [9].In Figure1, two groups of manually prepared image mosaics are shown. In Figure 1(a), the example surfaces display normal concrete finish whereas in Figure 1(b) concrete surfaces with evident cracks are shown. Generally, the surface objects in optical images display distinct features in terms of their spectral intensities and their topographical patterns. More specifically, the normal concrete surface generally can be viewed as a type of texture with lighter intensities and visually identifiable spotted air pockets. Cracks, on the other hand, are usually darker with approximately uniform spectral intensities. In addition cracks in images can be geometrically described as elongated “narrow” 2D objects with closed boundaries. From an image analysis point of view, the concrete surface in images provides rich and diverse patterns for exploring different analysis methods.Example concrete surface images: (a) concrete surfaces with normal finish; (b) concrete surfaces with evident cracks.
(a)(b)


### 1.1. Problem Definition and Previous Work

Detection of cracks with images of concrete surfaces has been explored by many researchers using a variety of methods. However, a comprehensive image-based damage identification framework is still lacking in the literature. Such a framework needs to provide a complete solution to the final goal, which is damage quantification. Herein, we define damage identification as a three-level probing process, which includes(1)
temporal monitoring of the inception or propagation of damage using time-series images (the occurrence determination problem);(2)
detection of the boundary of damage at each image frame (the localization problem);(3)
quantitative, geometric characterization of individual damage objects (the quantification problem).If only static images instead of multiple images in a temporal sequence are available, identification reduces to the problems of damage boundary localization and damage quantification (Levels 2 and 3). In previous work, efforts have stressed more on boundary detection [4, 5, 10–12]. Among these efforts, the simplest solution is to determine the intensity gradients across the image, and boundaries between different image objects are extracted based on local variations of gradient magnitudes and orientations. This detection process in image processing is termed gradient-based edge detection. Common gradient-based edge detectors include the Roberts, Sobel, Prewitt, and Canny edge detector [13]. In these detectors, the approximation of intensity gradients is based on the concept of convolution with a set of directional derivative masks at a specified scale. To deal with the multiresolution nature of edges, researchers have proposed an automatic scale-selection approach [14]; however, this approach has not been investigated in the crack detection literature. Another approach to overcom the multiresolution problem is to use wavelet-based image transform. Wavelets are special basis functions, localized in both spatial and frequency domains. Short-duration wavelets can isolate fine variations of image intensity while long-duration wavelets can isolate coarse variations. Image structures, such as edges, at different scales are obtained by computing wavelet coefficients at different resolutions [15]. Subsequently, the obtained wavelet coefficients are thresholded to determine edges at different scales.


### 1.2. Scope of This Work

Regardless of the previous methods adopted, a major difficulty lies in the postprocessing steps. When geometric quantities of detected cracks are to be extracted, the processing usually becomes a tedious step with many manual interventions. This is due to the lack of a geometric model for the cracks. On the other hand, there is no previous work that has addressed an image-based monitoring problem in the context of concrete surface crack identification. In what follows, first, we discuss the implementation challenges when raw laboratory or field images are used. Subsequently, a unified framework is proposed, which involves (i) a manifold-distance computing algorithm for image-based damage monitoring and (ii) a level-set active contour model-based method for damage boundary detection and geometric quantification. To demonstrate the use of the proposed framework, two application examples are presented. One example uses a field image, and the other uses time-series images taken in a controlled laboratory environment. The later provides direct comparison with analog measurements for validation of the analysis results.


## 1.1. Problem Definition and Previous Work

Detection of cracks with images of concrete surfaces has been explored by many researchers using a variety of methods. However, a comprehensive image-based damage identification framework is still lacking in the literature. Such a framework needs to provide a complete solution to the final goal, which is damage quantification. Herein, we define damage identification as a three-level probing process, which includes(1)
temporal monitoring of the inception or propagation of damage using time-series images (the occurrence determination problem);(2)
detection of the boundary of damage at each image frame (the localization problem);(3)
quantitative, geometric characterization of individual damage objects (the quantification problem).If only static images instead of multiple images in a temporal sequence are available, identification reduces to the problems of damage boundary localization and damage quantification (Levels 2 and 3). In previous work, efforts have stressed more on boundary detection [4, 5, 10–12]. Among these efforts, the simplest solution is to determine the intensity gradients across the image, and boundaries between different image objects are extracted based on local variations of gradient magnitudes and orientations. This detection process in image processing is termed gradient-based edge detection. Common gradient-based edge detectors include the Roberts, Sobel, Prewitt, and Canny edge detector [13]. In these detectors, the approximation of intensity gradients is based on the concept of convolution with a set of directional derivative masks at a specified scale. To deal with the multiresolution nature of edges, researchers have proposed an automatic scale-selection approach [14]; however, this approach has not been investigated in the crack detection literature. Another approach to overcom the multiresolution problem is to use wavelet-based image transform. Wavelets are special basis functions, localized in both spatial and frequency domains. Short-duration wavelets can isolate fine variations of image intensity while long-duration wavelets can isolate coarse variations. Image structures, such as edges, at different scales are obtained by computing wavelet coefficients at different resolutions [15]. Subsequently, the obtained wavelet coefficients are thresholded to determine edges at different scales.


## 1.2. Scope of This Work

Regardless of the previous methods adopted, a major difficulty lies in the postprocessing steps. When geometric quantities of detected cracks are to be extracted, the processing usually becomes a tedious step with many manual interventions. This is due to the lack of a geometric model for the cracks. On the other hand, there is no previous work that has addressed an image-based monitoring problem in the context of concrete surface crack identification. In what follows, first, we discuss the implementation challenges when raw laboratory or field images are used. Subsequently, a unified framework is proposed, which involves (i) a manifold-distance computing algorithm for image-based damage monitoring and (ii) a level-set active contour model-based method for damage boundary detection and geometric quantification. To demonstrate the use of the proposed framework, two application examples are presented. One example uses a field image, and the other uses time-series images taken in a controlled laboratory environment. The later provides direct comparison with analog measurements for validation of the analysis results.


## 2. Practical Challenges

Although image samples shown in Figure1 are valuable in the context of building library of damage patterns and testing new crack detection methods, significant practical difficulties exist when raw images captured in the laboratory or field are considered. Considering Figure 2, an image taken in the field of a shear crack at a bridge abutment, although it is visually easy to identify the crack in the image, many noisy artifacts are found, which prevent from realizing an automated, algorithm-based identification. These noisy artifacts include, for example, decoloration, occlusion, and solar illumination-related photometric variabilities (e.g., shadows). In addition, the relative large size of image compared to the cracked subregion impedes a computationally efficient implementation.Figure 2
A field image taken under the over-crossing bridge of Wailua river, Kapaa, Hawaii, which shows a long crack on the bridge abutment.Even in a controlled environment, for example, laboratory testing, although test configuration is adjustable, the optimal one is often not for imaging. Hence, use of images for a quantitative damage identification faces significant challenges as well. Figures3(a) and 3(b) show two images, which display the same concrete wall specimen that has been subject to biaxial loading at two different cyclic displacement peaks [16]. From Figure 3(a) to Figure 3(b), it can be seen that additional cracks were developed due to the increased displacement demand. One may attempt to use a boundary detection method to locate the cracks in either of the two images. However, it is not difficult to reason that the obtained boundaries would include many falsealarms (e.g., those associated with the diagonal displacement transducers), thereby preventing from an automated damage quantification.A laboratory biaxial test specimen under cyclic loading (a) at drift ratio ofγ=0.21% and (b) γ=0.41% (images courtesy of Soppe and Hutchinson, [16]).
(a)(b)


## 3. Proposed Framework

A straightforward approach to addressing the aforementioned adverse factors may be achieved by performing nonsupervised image segmentation. However, in most image segmentation methods, the number of image modalities must be known a priori [17]. Therefore, for a field image as shown in Figure 2 where the number of different texture regions is arbitrary, segmentation approach is not practical.This difficulty may be solved if continuous imaging through routine inspection is conducted, or more ideally, permanent camera sensors are installed to monitor the surface change of a critical structural component. Therefore, multitemporal images for the same structural surface can be obtained. This situation admits a digital change detection-based solution [18]. The applicability of a change detection method is based on a reasonable assumption that two temporally adjacent images may only differ in illumination-induced photometric distortion (e.g., solar or laboratory illumination variations), structural motion or camera motion-induced geometric distortion, and most importantly, structural damage. Meanwhile, most of the image artifacts in the scene, including decoloration and other texture objects, remain in the image sequence; hence, simple differencing operations can theoretically remove their influence if the image distortions are handled appropriately. To deal with geometric or photometric distortions, relatively simple models can be used; if these models or model parameters are estimated, a measure that is sensitive to the occurrence of structural damage by distortion can be subsequently extracted (i.e., distortion-invariant measure). As shown in the later formulation, this is the basic idea underlying the proposed manifold distance-based structural damage monitoring.From the reasoning above, one can see that an appropriate modeling for geometric or photometric distortions becomes the critical path in realizing a distortion-invariant damage detection. For mitigating illumination-induced photometric distortion, a common technique is normalization of spectral intensities in small image regions through subdividing the image domain to many small subdomains. For geometric distortion, the prevailing adverse factor that complicates this concept is the motion of structural objects. For a structure under loading, any point within the structure potentially undergoes some spatial motion across time. To extract a robust feature representing crack damage from a time series of images, modeling the motion of structural objects in the images is a critical issue.With the background information introduced in the above, we propose an image-based structural damage monitoring and quantification framework. As illustrated in Figure4, the proposed framework has three operational components. First, image acquisition is conducted. The success of this element relies upon reliable assembly of hardware and careful design of acquisition software, which is not in the scope of this paper. It should be noted that computing is necessary in this step, which usually includes camera calibration and image correction. When a 3D scene is projected into an image plane, geometric distortion is associated with the projection determined by the intrinsic and extrinsic parameters of the digital camera [17]. In practice, a process known as camera calibration has to be conducted. Once calibration parameters are estimated, acquired images can be easily corrected to remove this distortion. In the following, we describe the subsequent two key elements to the framework, which are temporal image-based damage monitoring and geometric image-based quantification of damage.Figure 4
Proposed image-based structural damage monitoring and quantification framework.


## 4. Manifold-Distance Concrete CrackMonitoring Using Time-Series Images

We approach our image-based damage monitoring problem as a temporal change detection problem, wherein structural damage is treated as a type of local anomaly with its amount measurable between temporal image frames. Due to the aforementioned adverse effects of structural motion in images, a seemly straightforward approach—image differencing—is not appropriate since it requires strict spatial alignment between the temporal image pair. A simple remedy is to first extract the underlying motion parameters, then with the estimated motion parameters, the temporal image frames can be coregistered, such that image differencing can be performed, and motion-invariant features useful for representing concrete cracks can be defined. This approach was proposed in our previous work in [19]. In this section, we will briefly present this combined motion estimation and feature extraction procedure.


### 4.1. Assumptions about Structural Motion in Images

Theoretically, the displacement field of the structural member is governed by the material's stress-strain relation and the geometric boundary condition of the structural member, which may be in a form or a combination of bending, axial, shear deformation, and rigid-body motion. However, the motion of the structural members cannot be fully captured without loss of accuracy due to the limit of image resolution. Suppose that a resolution of1mm/pixel may suffice for the purpose of monitoring crack damage for concrete structures. With this resolution and a common image size of 480×640, the field of view (FoV) on the surface of the target member is approximately a 480×640mm2 rectangular zone. For normal-weight concrete under compression, the ultimate compressive strain is about 0.003 [20]. At this limit state, the maximum relative deformation due to compression within this FoV is about 1.92mm. Since the original image domain is usually divided into multiple subdomains (e.g., 10×10 subdomains) in order to minimize photometric distortions, the maximum within a subdomain is much less than 1mm. A direct consequence is that the relative deformation captured in an image subdomain can be reasonably assumed to be negligible. It is noted that the above reasoning is based on the maximum compressive strain; however, cracks are due to local extensive tensile strain, which is much smaller. This strengthens that local deformation is generally not captured by optical cameras.We also limit the study in this paper to consider the structural motion plane only approximately parallel with the image plane, resulting in a linearly scaled planar motion in the image domain. Therefore, we assume that the structural motion considered is rigid-body motion only, which includes the two orthogonal translations and one in-plane rotation.


### 4.2. Formulation of Manifold Distance Measure

To present the formulation of both motion estimation and damage feature extraction using a time series of images, we use two notations to describe an image acquired at timet. First, u(x,t) is used to denote a continuous image function or an intensity value at the location x and the time t. Second, if a digital image is treated as a vector, it is denoted as u(t) with dimension of N indicating the number of pixels. For a sequence of image vectors across time, they can be denoted by {u(t0),u(t1),…,u(ti),…}.Given a current intensity value at timeti by u(x,ti) and its temporally previous value by u(x,ti-j)(j≥1), the current intensity value u(x,ti) can be obtained by evolving u(x,ti-j) through a spatial transformation of the coordinates x=(x,y) if there is no photometric transformation. Without loss of generality, we use ut(x) to replace u(x,ti) and u(x) to replace u(x,ti-j) for notational simplicity. Therefore, the above transformation is written as(1)ut(x)=u[ϕ(x,θ)],
where the function ϕ:ℝ2×ℝM→ℝ2 defines a generic transformation model ϕ(x,θ)=[ϕ(x,θ),ψ(x,θ)]T parameterized by θ (an M×1 vector quantity with M≪N). This expression in the visual tracking literature is traditionally called the image constancy assumption [21]. The parameterized transformation model ϕ(x,θ) is differentiable with respect to both x and θ. Additionally, an initialization is usually adopted for θ=0, which is u[ϕ(x,0)]=u(x). If the assumed transformation holds globally in the spatial domain 𝒳, then any image can be fully determined by another image and the spatial transformation model ϕ(x,θ) between them. Formally, images like these can be treated as a lower-dimensional manifold, denoted by u(θ), embedded in the original high-dimensional observation space ℝN. In this case, the intrinsic dimension of such a manifold space is the number of motion parameters M, which has six parameters for an affine transformation model or three parameters for a rigid-body motion model. Note that the intrinsic dimension M is generally far less than the dimension of the image observation space N.The geometric structure of an image manifold is usually augmented and complicated by some lighting-induced photometric perturbations and the development of cracks across the image sequence. Technically, the photometric perturbation can be modeled by some parametric intensity transformation [22]. As a result, more intrinsic dimensions are added into the previous M dimensions. In this paper, we accept the fact that lighting is time varying across the image sequence however, we assume that it is statistically stationary within a small image domain at one acquisition time. Therefore, this type of perturbation can be removed through a simple intensity standardization (i.e., subtracting the sample mean intensity from the raw image) to remove constant variation. Cracks may also be viewed as a sort of geometric deformation occurring at some intrinsic dimensions. However, the irregular topological shape of cracks in the image domain indicates that they have highly nonlinear curvature locally in both observation and manifold spaces. Therefore, it is hard to find an explicit parametric model defined globally in the image domain to represent the cracks.Leaving the potential crack development not modeled at a pixelx, we treat it as a residue expressed by |u[ϕ(x,θ*)]-ut(x)|2 and denote the total residue over the spatial domain by r(θ*,t)=(∫𝒳dx|u[ϕ(x,θ*)]-ut(x)|2)1/2, where θ* contains the ground truth values of motion parameters that fully accounts for the spatial transformation occurring in the image frames from u to ut. Ideally, if there are no cracks developing at time t (more exactly, within a duration from ti-j to ti), but only spatial transformation, the total residue will yield a value of zero. If crack damage does occur, the total residue should be larger than zero, and the same amount of cracks corresponds to the same value of the total residue yet with possibly different motion parameters θ*s. Therefore, r(θ*,t) is theoretically an ideal feature of crack occurrence that is invariant to the underlying motion of the structural component.If one replacesθ* with an unknown motion parameter vector θ in the previous definition of the total residue, a parameterized total residue is yielded, r(θ,t), which is essentially a varying Euclidean distance (ED) between a parameterized image vector with a base at u and the reference image vector ut. On the other hand, if the image vectors are viewed in the manifold space, then this distance metric quantifies the distance between the manifold u(θ) and the reference point ut. It is reasoned that given some potential crack occurrence left unmodeled and a certain amount of motion determined by θ*, the quantity r(θ*,t) is therefore the minimum distance between u(θ) and ut, which is called manifold distance (MD) in the literature [23, 24]. In this paper, we define an MD measure at time t given u and ut:(2)MD(t)=r(θ*,t)≡inf{r(θ,t)θ∣∈ℝM}.
In particular, if θ* is equal to 0, that is, there is no spatial motion, MD(t) degrades to the constant Euclidean distance between u and ut, which is ED(t)=(∫𝒳dx|u(x)-ut(x)|2)1/2. Therefore, the seeking of a motion-invariant crack damage feature entails a proper estimate of θ*, which leads to the following optimization problem:(3)minimizer2(θ,t)=∫𝒳dx|u[ϕ(x,θ)]-ut(x)|2,(4)subjecttou[ϕ(x,0)]=u(x).


### 4.3. Estimation of Motion Parameters

The objective functionr2(θ,t) in (3) unlikely owns a convex surface due to the highly nonlinear nature of an image function in terms of both x and θ. To address this, one may approximate the functional u[ϕ(x,θ)] in the manifold space at θ=0 by carrying out a Taylor series expansion by considering the parameterized model ϕ(x,θ). Suppose that the translation vector is denoted by (θ1,θ2)T, and the rotation angle is θ3, which constitutes the motion parameter vector θ=(θ1,θ2,θ3)T. The motion function is therefore(5)ϕ(x,θ)=(cosθ3sinθ3-sinθ3cosθ3){xy}+{θ1θ2}={θ1+xcosθ3+ysinθ3θ2-xsinθ3+ycosθ3}={ϕ(x,θ)ψ(x,θ)}.
With some steps of vector calculus computation, one arrives at an approximate form of (3):(6)r2¯(θ,t)=C0+C1θ1+C2θ2+C3θ3+C4θ12+C5θ22+C6θ32+C7θ1θ2+C8θ1θ3+⋯+C9θ2θ3+C10θ1θ32+C11θ2θ32+C12θ33+C13θ34,
where the constant coefficients C0~C13 are integrations over the image domain 𝒳. Detailed derivation of these constant coefficients may be found in [19]. The approximate squared residue in (6) is a fourth-order polynomial; hence it provides a tractable basis for studying the solution of this optimization problem defined in (3) and (4). Although analytical forms of the solution are available, they are extremely lengthy. A more common approach in practice is to use the standard Newton method for the numerical solutions, which is expressed an iterative procedure:(7)θ(k+1)=θ(k)-H-1(θ(k),t)g(θ(k),t),
where H(·) and g(·) are the Hessian matrix and the gradient vector of (6) with respect to θ, respectively, and their analytical forms are found in [19]. Upon convergence of the iteration, the MD measurement is computed by(8)MD(t)≃[r2¯(θ̂,t)]1/2,
where θ̂ is the estimated motion parameter vector from (7).If one is given a time series of images{u(t0),u(t1),u(t2),…,…}, there are two ways to compute a sequence of time-varying MD measurements. First, one can compute an MD measurement for any temporally adjacent image frames, that is, MD(ti) between u(ti-1) and u(ti). In this paper, MD measurements resulting from this method can be used to monitor the inception of concrete cracks. Another method is that upon cracks having appeared at time ti, the subsequent MD measurements at time tj are all computed relative to the time ti, that is, MD(tj) is obtained based on u(ti) and u(tj) where ti is fixed. This approach can be used to monitor the continuous growth of cracks.


### 4.4. Damage Localization

The MD measurement computed at an arbitrary time in essence offers alump sum indication of occurrence for the concrete cracks in an image domain. However, if the image domain is relatively large, it does not provide an indication regarding the location of cracks. In this paper, we localize concrete cracks in images by dividing the 2D image domain into multiple subdomains before computing the MD measurements. Another advantage of this treatment is that it can minimize the disturbance of photometric distortions. The criterion for the subdivision is that the size of the subdomain should be much larger than the potential motion change between two temporal image patches upon which an MD measurement is computed. In practice, one needs a priori knowledge about the magnitude of the potential motion change in order to make such subdivision.Suppose that the image domain is divided intoK×L subdomains, therefore, a matrix (K×L) of MD measurements can be obtained at an arbitrary acquisition time. With the MD measurements at different times, one can perform a nonsupervised outlier analysis to determine which MD measurement statistically guarantees an instance of crack occurrence. The outlier analysis, in its simplest form, is implemented by assuming that the underlying features follow a normal distribution (for more rigorous treatment about modeling the distribution of outliers, one may refer to [25]). The subdomains that are determined statistically to have crack occurrence hence provide a description of the initial location of damage.


## 4.1. Assumptions about Structural Motion in Images

Theoretically, the displacement field of the structural member is governed by the material's stress-strain relation and the geometric boundary condition of the structural member, which may be in a form or a combination of bending, axial, shear deformation, and rigid-body motion. However, the motion of the structural members cannot be fully captured without loss of accuracy due to the limit of image resolution. Suppose that a resolution of1mm/pixel may suffice for the purpose of monitoring crack damage for concrete structures. With this resolution and a common image size of 480×640, the field of view (FoV) on the surface of the target member is approximately a 480×640mm2 rectangular zone. For normal-weight concrete under compression, the ultimate compressive strain is about 0.003 [20]. At this limit state, the maximum relative deformation due to compression within this FoV is about 1.92mm. Since the original image domain is usually divided into multiple subdomains (e.g., 10×10 subdomains) in order to minimize photometric distortions, the maximum within a subdomain is much less than 1mm. A direct consequence is that the relative deformation captured in an image subdomain can be reasonably assumed to be negligible. It is noted that the above reasoning is based on the maximum compressive strain; however, cracks are due to local extensive tensile strain, which is much smaller. This strengthens that local deformation is generally not captured by optical cameras.We also limit the study in this paper to consider the structural motion plane only approximately parallel with the image plane, resulting in a linearly scaled planar motion in the image domain. Therefore, we assume that the structural motion considered is rigid-body motion only, which includes the two orthogonal translations and one in-plane rotation.


## 4.2. Formulation of Manifold Distance Measure

To present the formulation of both motion estimation and damage feature extraction using a time series of images, we use two notations to describe an image acquired at timet. First, u(x,t) is used to denote a continuous image function or an intensity value at the location x and the time t. Second, if a digital image is treated as a vector, it is denoted as u(t) with dimension of N indicating the number of pixels. For a sequence of image vectors across time, they can be denoted by {u(t0),u(t1),…,u(ti),…}.Given a current intensity value at timeti by u(x,ti) and its temporally previous value by u(x,ti-j)(j≥1), the current intensity value u(x,ti) can be obtained by evolving u(x,ti-j) through a spatial transformation of the coordinates x=(x,y) if there is no photometric transformation. Without loss of generality, we use ut(x) to replace u(x,ti) and u(x) to replace u(x,ti-j) for notational simplicity. Therefore, the above transformation is written as(1)ut(x)=u[ϕ(x,θ)],
where the function ϕ:ℝ2×ℝM→ℝ2 defines a generic transformation model ϕ(x,θ)=[ϕ(x,θ),ψ(x,θ)]T parameterized by θ (an M×1 vector quantity with M≪N). This expression in the visual tracking literature is traditionally called the image constancy assumption [21]. The parameterized transformation model ϕ(x,θ) is differentiable with respect to both x and θ. Additionally, an initialization is usually adopted for θ=0, which is u[ϕ(x,0)]=u(x). If the assumed transformation holds globally in the spatial domain 𝒳, then any image can be fully determined by another image and the spatial transformation model ϕ(x,θ) between them. Formally, images like these can be treated as a lower-dimensional manifold, denoted by u(θ), embedded in the original high-dimensional observation space ℝN. In this case, the intrinsic dimension of such a manifold space is the number of motion parameters M, which has six parameters for an affine transformation model or three parameters for a rigid-body motion model. Note that the intrinsic dimension M is generally far less than the dimension of the image observation space N.The geometric structure of an image manifold is usually augmented and complicated by some lighting-induced photometric perturbations and the development of cracks across the image sequence. Technically, the photometric perturbation can be modeled by some parametric intensity transformation [22]. As a result, more intrinsic dimensions are added into the previous M dimensions. In this paper, we accept the fact that lighting is time varying across the image sequence however, we assume that it is statistically stationary within a small image domain at one acquisition time. Therefore, this type of perturbation can be removed through a simple intensity standardization (i.e., subtracting the sample mean intensity from the raw image) to remove constant variation. Cracks may also be viewed as a sort of geometric deformation occurring at some intrinsic dimensions. However, the irregular topological shape of cracks in the image domain indicates that they have highly nonlinear curvature locally in both observation and manifold spaces. Therefore, it is hard to find an explicit parametric model defined globally in the image domain to represent the cracks.Leaving the potential crack development not modeled at a pixelx, we treat it as a residue expressed by |u[ϕ(x,θ*)]-ut(x)|2 and denote the total residue over the spatial domain by r(θ*,t)=(∫𝒳dx|u[ϕ(x,θ*)]-ut(x)|2)1/2, where θ* contains the ground truth values of motion parameters that fully accounts for the spatial transformation occurring in the image frames from u to ut. Ideally, if there are no cracks developing at time t (more exactly, within a duration from ti-j to ti), but only spatial transformation, the total residue will yield a value of zero. If crack damage does occur, the total residue should be larger than zero, and the same amount of cracks corresponds to the same value of the total residue yet with possibly different motion parameters θ*s. Therefore, r(θ*,t) is theoretically an ideal feature of crack occurrence that is invariant to the underlying motion of the structural component.If one replacesθ* with an unknown motion parameter vector θ in the previous definition of the total residue, a parameterized total residue is yielded, r(θ,t), which is essentially a varying Euclidean distance (ED) between a parameterized image vector with a base at u and the reference image vector ut. On the other hand, if the image vectors are viewed in the manifold space, then this distance metric quantifies the distance between the manifold u(θ) and the reference point ut. It is reasoned that given some potential crack occurrence left unmodeled and a certain amount of motion determined by θ*, the quantity r(θ*,t) is therefore the minimum distance between u(θ) and ut, which is called manifold distance (MD) in the literature [23, 24]. In this paper, we define an MD measure at time t given u and ut:(2)MD(t)=r(θ*,t)≡inf{r(θ,t)θ∣∈ℝM}.
In particular, if θ* is equal to 0, that is, there is no spatial motion, MD(t) degrades to the constant Euclidean distance between u and ut, which is ED(t)=(∫𝒳dx|u(x)-ut(x)|2)1/2. Therefore, the seeking of a motion-invariant crack damage feature entails a proper estimate of θ*, which leads to the following optimization problem:(3)minimizer2(θ,t)=∫𝒳dx|u[ϕ(x,θ)]-ut(x)|2,(4)subjecttou[ϕ(x,0)]=u(x).


## 4.3. Estimation of Motion Parameters

The objective functionr2(θ,t) in (3) unlikely owns a convex surface due to the highly nonlinear nature of an image function in terms of both x and θ. To address this, one may approximate the functional u[ϕ(x,θ)] in the manifold space at θ=0 by carrying out a Taylor series expansion by considering the parameterized model ϕ(x,θ). Suppose that the translation vector is denoted by (θ1,θ2)T, and the rotation angle is θ3, which constitutes the motion parameter vector θ=(θ1,θ2,θ3)T. The motion function is therefore(5)ϕ(x,θ)=(cosθ3sinθ3-sinθ3cosθ3){xy}+{θ1θ2}={θ1+xcosθ3+ysinθ3θ2-xsinθ3+ycosθ3}={ϕ(x,θ)ψ(x,θ)}.
With some steps of vector calculus computation, one arrives at an approximate form of (3):(6)r2¯(θ,t)=C0+C1θ1+C2θ2+C3θ3+C4θ12+C5θ22+C6θ32+C7θ1θ2+C8θ1θ3+⋯+C9θ2θ3+C10θ1θ32+C11θ2θ32+C12θ33+C13θ34,
where the constant coefficients C0~C13 are integrations over the image domain 𝒳. Detailed derivation of these constant coefficients may be found in [19]. The approximate squared residue in (6) is a fourth-order polynomial; hence it provides a tractable basis for studying the solution of this optimization problem defined in (3) and (4). Although analytical forms of the solution are available, they are extremely lengthy. A more common approach in practice is to use the standard Newton method for the numerical solutions, which is expressed an iterative procedure:(7)θ(k+1)=θ(k)-H-1(θ(k),t)g(θ(k),t),
where H(·) and g(·) are the Hessian matrix and the gradient vector of (6) with respect to θ, respectively, and their analytical forms are found in [19]. Upon convergence of the iteration, the MD measurement is computed by(8)MD(t)≃[r2¯(θ̂,t)]1/2,
where θ̂ is the estimated motion parameter vector from (7).If one is given a time series of images{u(t0),u(t1),u(t2),…,…}, there are two ways to compute a sequence of time-varying MD measurements. First, one can compute an MD measurement for any temporally adjacent image frames, that is, MD(ti) between u(ti-1) and u(ti). In this paper, MD measurements resulting from this method can be used to monitor the inception of concrete cracks. Another method is that upon cracks having appeared at time ti, the subsequent MD measurements at time tj are all computed relative to the time ti, that is, MD(tj) is obtained based on u(ti) and u(tj) where ti is fixed. This approach can be used to monitor the continuous growth of cracks.


## 4.4. Damage Localization

The MD measurement computed at an arbitrary time in essence offers alump sum indication of occurrence for the concrete cracks in an image domain. However, if the image domain is relatively large, it does not provide an indication regarding the location of cracks. In this paper, we localize concrete cracks in images by dividing the 2D image domain into multiple subdomains before computing the MD measurements. Another advantage of this treatment is that it can minimize the disturbance of photometric distortions. The criterion for the subdivision is that the size of the subdomain should be much larger than the potential motion change between two temporal image patches upon which an MD measurement is computed. In practice, one needs a priori knowledge about the magnitude of the potential motion change in order to make such subdivision.Suppose that the image domain is divided intoK×L subdomains, therefore, a matrix (K×L) of MD measurements can be obtained at an arbitrary acquisition time. With the MD measurements at different times, one can perform a nonsupervised outlier analysis to determine which MD measurement statistically guarantees an instance of crack occurrence. The outlier analysis, in its simplest form, is implemented by assuming that the underlying features follow a normal distribution (for more rigorous treatment about modeling the distribution of outliers, one may refer to [25]). The subdomains that are determined statistically to have crack occurrence hence provide a description of the initial location of damage.


## 5. Level-Set Representation of Concrete Cracks and Quantification

The image-processing procedure presented in the previous section offers a means of monitoring the inception or propagation of structural damage; meanwhile, a simple method for approximately locating the potential cracks is given. In the following, a method for geometric boundary detection and detailed damage quantification of concrete surface cracks is presented, in which we proposed a level-set active contour model-based approach for this purpose.


### 5.1. Active Contour Models for Boundary Detection

The basic idea in active contour models is to evolve a curve subject to constraints in a given image, and eventually the curve converges to the closed boundary of objects in the image. Since the invention of the classical active contour model [26], many geometric versions have been proposed, such as the Geodesic active contour model [27]. In the expressions of the resulting partial differential equations in these contour models, a stopping function is used to stop the curve evolution, which is usually defined based on image gradients.In cases of image-based detection of concrete cracks, another cue is that concrete surface images can be approximated by piecewise, multiphase cartoon images. For a cracked concrete surface image, the simplest form contains two phases, the normal background or the cracked areas. This implication inspires us to use the Mumford-Shah segmentation model [28], which is further extended by Chan and Vese in a level-set framework [29]. A level-set method is a general technique for evolving curves (surfaces) that undergo complex topographical changes such as merging and pinching [30, 31]. It has been proven efficient in recent years in problems ranging from tracking, modeling, and simulating motion of dynamic surfaces in the fields of graphics, image analysis, fluid dynamics, fracture mechanics, and optimal topology control.The Chan-Vese level-set active contour model can detect object boundaries without relying on image gradients, eliminating the need to define a gradient-based stopping function. In addition, the most appealing aspect of the Chan-Vese model is that it has a level-set representation, which offers the following benefits: (i) automatic handling of topographical merging and pinching of evolving contours, (ii) detection of interior contours, (iii) initial contours can be anywhere in images, and (iv) its level-set representation offers the capability to extract geometric quantities of concrete damage within a formal mathematical framework.


### 5.2. Level-Set Representation

Suppose that an original concrete surface image is denoted byu0(x,y) defined in a 2D domain Ω; an evolving curve attempting to locate objects of interest in the image is denoted by Γ; one further denotes the open sub-domains inside and outside of Γ by Ω+ and Ω-, respectively. In the spirit of level-set methods, the curve Γ is embedded as the zero level set of a Lipschitz continuous function ϕ:Ω→ℝ, denoted by {ϕ(x,y)=0}. Correspondingly, Ω+ is written as {ϕ(x,y)>0} and Ω- as {ϕ(x,y)<0}. Let H(·) be the Heaviside function. The energy Chan-Vese functional is(9)E[c1,c2,ϕ(x,y,t)∣u0(x,y)]=∫Ω|u0(x,y)-c1|2H(ϕ)dxdy+∫Ω|u0(x,y)-c2|2[1-H(ϕ)]dxdy+μ∫Ω|∇H(ϕ)|dxdy,
where c1 and c2 are the average intensities across the inside (Ω+) and the outside (Ω-) regions, μ is the associated length scale, and t is a pseudo time variable. Generally, the parameter μ controls the smoothness of the curve; therefore μ|Γ| is a regularizer in the above functional. The time-dependent (a pseudocomputational time will be introduced in evolving Γ) constants c1 and c2 vary as the curve evolves in u0(x,y).By minimizing (9) with respect to c1, c2, and ϕ, one arrives at the following Euler-Lagrange equations:(10)∂ϕ∂t=δ(ϕ)[μ∇·(∇ϕ|∇ϕ|)-|u0(x,y)-c1|2+|u0(x,y)-c2|2(∇ϕ|∇ϕ|)],(11)c1(t)=∫Ωu0(x,y)H(ϕ)dxdy∫ΩH(ϕ)dxdy,(12)c2(t)=∫Ωu0(x,y)[1-H(ϕ)]dxdy∫Ω[1-H(ϕ)]dxdy,
where δ is the Dirac delta function defined by δ(ϕ)=(d/dϕ)H(ϕ). With the initial contour denoted by ϕ(x,y,0)=ϕ0(x,y), one can see that (10)–(12) construct a recursive procedure when solving for the implicit level-set function ϕ.


### 5.3. Extraction of Geometric Properties

In this paper, the numerical schemes for solving (10) are not included; one may refer to [32] for more details. To summarize, by iteratively integrating (10) numerically, the level-set function ϕ(x,y,t) will eventually converge at a time, denoted by te; hence, we obtain the final level-set function ϕe=ϕ(x,y,te). Many geometric quantities of detected cracks, such as perimeter (total length of closed boundary), area, curvature, and normal direction of crack boundaries, can be computed directly with ϕe. However, we believe that the area, perimeter, and width are of primary engineering interest.


#### 5.3.1. Area and Perimeter of Cracks

In the Chan-Vese model, the computation of areas and perimeters of detected objects has been embedded in (10)–(12): the area A=∫ΩH(ϕe)dxdy, and the perimeter P=∫Ω|∇H(ϕe)|dxdy. However, for the purpose of computing areas and perimeters, a more “smeared-out” approximation will offer more numerical accuracy, which is(13)Hθ(ϕ)={1ifϕ>θ,0ifϕ<-θ,12[1+ϕθ+1πsin(πϕθ)]if|ϕ|≤θ,
where the parameter θ is selected comparable to the used space step. In the work of Osher and Sethian [31], θ=1.5h is recommended.


#### 5.3.2. Extraction of Crack Width

Another important geometric quantity of engineering interest is crack width. However, unlike the computation of areas or perimeters for concrete cracks, the extraction of crack width is a nontrivial problem. First, the concept ofwidth is empirical. Recall that we adopt a nonrigorous way to describe cracks, namely, cracks in images are spatially narrow, elongated objects with strong contrast with the background, where the qualitative terms “narrow” or “elongated” are related to “width” without explicit definitions. In engineering practice, first, one needs to manually identify two end points for a single crack or multiple end points for a crack with a complex pattern; then by starting any of these end points, widths can be measured along the cracks continuously. To compute the crack widths analytically, a rigorous definition for crack width is necessary.Considering a simple case (Figure5), where a crack, with its closed boundary denoted by Γ, has two end points A:(xa,ya) and B:(xb,yb), therefore, the boundary of this crack can be split into two open boundaries, Γ1 and Γ2, which both start at A and end at B. The width w at a pixel (x1,y1) along one side of the crack boundaries, say, Γ1, is the shortest Euclidean distance (L2 distance) from (x1,y1) to a point (x2,y2) in the other side of boundaries Γ2. For cracks with complex patterns, for example, a continuous crack with “T”, “+” shapes, or other multiple angle splits, more than two end points have to be specified. One can break these cracks into several segments, wherein each segment has two end points, hence two sides of boundaries. The Euclidean distance-based definition for widths associated with these crack segments can be defined similarly.Figure 5
Definition of crack width: the distance from(x1,y1) to (x2,y2) is the shortest compared with others, for example, from (x1,y1) to (x2′,y2′); the former defines the crack width w at (x1,y1).


#### 5.3.3. Automated Extraction Method

The preceding definition of crack width may give rise to an exact solution to compute widths with known end points for a given crack. However, as the topographical patterns of cracks become more complex, the manual intervention for specifying the possible end points will be more tedious. Herein, we propose an approximate method for width extraction. This approximate method relies on the extraction of centerlines of detected cracks as well as the construction of a so-called signed distance function (SDF). An SDF, denoted byφ(x,y), defines the distance values from (x,y) to the evolving zero level-set Γ(t)={ϕ(x,y,t)=0}, which are positive inside of Γ(t) and negative outside of Γ(t), and satisfies |∇φ|=1 everywhere. The steady state of φ(x,y) is obtained through solving(14)∂φ∂τ+sgn(ϕ)(|∇φ|-1)=0,
where sgn(ϕ)=2H(ϕ)-1 is a signum function, and τ is another pseudotime variable added in the dependents of φ. With the initial condition φ(x,y,0)=ϕ(x,y,t), the zero level set Γ(t) will be preserved during evolving φ(x,y,τ). For centerline extraction, although a level set-based extraction of centerlines is possible [33], we employ a simple morphological skeleton operator that can preserve the connectivity of the original topographical patterns [32]. Theoretically, for objects with general topographical shapes, morphological skeltonization does not usually produce the ideal centerlines; nonetheless, for narrow, elongated 2D objects, such as cracks, morphological skeltonization is sufficiently accurate.The proposed width extraction method is summarized in the following:(1)
fromHθ(ϕe)>0.5, a binary image b(x,y) is obtained, wherein the detected cracks are marked by unity;(2)
perform skeletonization onb(x,y), and a new binary image c(x,y) is produced wherein the centerlines are unity;(3)
solve for an SDFψ(x,y) using (14) with ψ(x,y,0)=sgn(ϕ)=2Hθ(ϕe)-1;(4)
the shortest distances from the centerlines to the zero level set are obtained by a simple Boolean operation, which is denoted byd(x,y)=ψ(x,y), forall(x,y)∈{(x,y)∣c(x,y)=1}, the widths along cracks are approximated by w(x,y)≃2d(x,y).In this approximate method, there is no restriction to the complexity of topographical patterns of cracks as well as no need to manually specify the end points of cracks. Regarding the accuracy of this method, the approximate widths will be slightly greater than the exact values for a simple crack similar to Figure5.


## 5.1. Active Contour Models for Boundary Detection

The basic idea in active contour models is to evolve a curve subject to constraints in a given image, and eventually the curve converges to the closed boundary of objects in the image. Since the invention of the classical active contour model [26], many geometric versions have been proposed, such as the Geodesic active contour model [27]. In the expressions of the resulting partial differential equations in these contour models, a stopping function is used to stop the curve evolution, which is usually defined based on image gradients.In cases of image-based detection of concrete cracks, another cue is that concrete surface images can be approximated by piecewise, multiphase cartoon images. For a cracked concrete surface image, the simplest form contains two phases, the normal background or the cracked areas. This implication inspires us to use the Mumford-Shah segmentation model [28], which is further extended by Chan and Vese in a level-set framework [29]. A level-set method is a general technique for evolving curves (surfaces) that undergo complex topographical changes such as merging and pinching [30, 31]. It has been proven efficient in recent years in problems ranging from tracking, modeling, and simulating motion of dynamic surfaces in the fields of graphics, image analysis, fluid dynamics, fracture mechanics, and optimal topology control.The Chan-Vese level-set active contour model can detect object boundaries without relying on image gradients, eliminating the need to define a gradient-based stopping function. In addition, the most appealing aspect of the Chan-Vese model is that it has a level-set representation, which offers the following benefits: (i) automatic handling of topographical merging and pinching of evolving contours, (ii) detection of interior contours, (iii) initial contours can be anywhere in images, and (iv) its level-set representation offers the capability to extract geometric quantities of concrete damage within a formal mathematical framework.


## 5.2. Level-Set Representation

Suppose that an original concrete surface image is denoted byu0(x,y) defined in a 2D domain Ω; an evolving curve attempting to locate objects of interest in the image is denoted by Γ; one further denotes the open sub-domains inside and outside of Γ by Ω+ and Ω-, respectively. In the spirit of level-set methods, the curve Γ is embedded as the zero level set of a Lipschitz continuous function ϕ:Ω→ℝ, denoted by {ϕ(x,y)=0}. Correspondingly, Ω+ is written as {ϕ(x,y)>0} and Ω- as {ϕ(x,y)<0}. Let H(·) be the Heaviside function. The energy Chan-Vese functional is(9)E[c1,c2,ϕ(x,y,t)∣u0(x,y)]=∫Ω|u0(x,y)-c1|2H(ϕ)dxdy+∫Ω|u0(x,y)-c2|2[1-H(ϕ)]dxdy+μ∫Ω|∇H(ϕ)|dxdy,
where c1 and c2 are the average intensities across the inside (Ω+) and the outside (Ω-) regions, μ is the associated length scale, and t is a pseudo time variable. Generally, the parameter μ controls the smoothness of the curve; therefore μ|Γ| is a regularizer in the above functional. The time-dependent (a pseudocomputational time will be introduced in evolving Γ) constants c1 and c2 vary as the curve evolves in u0(x,y).By minimizing (9) with respect to c1, c2, and ϕ, one arrives at the following Euler-Lagrange equations:(10)∂ϕ∂t=δ(ϕ)[μ∇·(∇ϕ|∇ϕ|)-|u0(x,y)-c1|2+|u0(x,y)-c2|2(∇ϕ|∇ϕ|)],(11)c1(t)=∫Ωu0(x,y)H(ϕ)dxdy∫ΩH(ϕ)dxdy,(12)c2(t)=∫Ωu0(x,y)[1-H(ϕ)]dxdy∫Ω[1-H(ϕ)]dxdy,
where δ is the Dirac delta function defined by δ(ϕ)=(d/dϕ)H(ϕ). With the initial contour denoted by ϕ(x,y,0)=ϕ0(x,y), one can see that (10)–(12) construct a recursive procedure when solving for the implicit level-set function ϕ.


## 5.3. Extraction of Geometric Properties

In this paper, the numerical schemes for solving (10) are not included; one may refer to [32] for more details. To summarize, by iteratively integrating (10) numerically, the level-set function ϕ(x,y,t) will eventually converge at a time, denoted by te; hence, we obtain the final level-set function ϕe=ϕ(x,y,te). Many geometric quantities of detected cracks, such as perimeter (total length of closed boundary), area, curvature, and normal direction of crack boundaries, can be computed directly with ϕe. However, we believe that the area, perimeter, and width are of primary engineering interest.


### 5.3.1. Area and Perimeter of Cracks

In the Chan-Vese model, the computation of areas and perimeters of detected objects has been embedded in (10)–(12): the area A=∫ΩH(ϕe)dxdy, and the perimeter P=∫Ω|∇H(ϕe)|dxdy. However, for the purpose of computing areas and perimeters, a more “smeared-out” approximation will offer more numerical accuracy, which is(13)Hθ(ϕ)={1ifϕ>θ,0ifϕ<-θ,12[1+ϕθ+1πsin(πϕθ)]if|ϕ|≤θ,
where the parameter θ is selected comparable to the used space step. In the work of Osher and Sethian [31], θ=1.5h is recommended.


### 5.3.2. Extraction of Crack Width

Another important geometric quantity of engineering interest is crack width. However, unlike the computation of areas or perimeters for concrete cracks, the extraction of crack width is a nontrivial problem. First, the concept ofwidth is empirical. Recall that we adopt a nonrigorous way to describe cracks, namely, cracks in images are spatially narrow, elongated objects with strong contrast with the background, where the qualitative terms “narrow” or “elongated” are related to “width” without explicit definitions. In engineering practice, first, one needs to manually identify two end points for a single crack or multiple end points for a crack with a complex pattern; then by starting any of these end points, widths can be measured along the cracks continuously. To compute the crack widths analytically, a rigorous definition for crack width is necessary.Considering a simple case (Figure5), where a crack, with its closed boundary denoted by Γ, has two end points A:(xa,ya) and B:(xb,yb), therefore, the boundary of this crack can be split into two open boundaries, Γ1 and Γ2, which both start at A and end at B. The width w at a pixel (x1,y1) along one side of the crack boundaries, say, Γ1, is the shortest Euclidean distance (L2 distance) from (x1,y1) to a point (x2,y2) in the other side of boundaries Γ2. For cracks with complex patterns, for example, a continuous crack with “T”, “+” shapes, or other multiple angle splits, more than two end points have to be specified. One can break these cracks into several segments, wherein each segment has two end points, hence two sides of boundaries. The Euclidean distance-based definition for widths associated with these crack segments can be defined similarly.Figure 5
Definition of crack width: the distance from(x1,y1) to (x2,y2) is the shortest compared with others, for example, from (x1,y1) to (x2′,y2′); the former defines the crack width w at (x1,y1).


### 5.3.3. Automated Extraction Method

The preceding definition of crack width may give rise to an exact solution to compute widths with known end points for a given crack. However, as the topographical patterns of cracks become more complex, the manual intervention for specifying the possible end points will be more tedious. Herein, we propose an approximate method for width extraction. This approximate method relies on the extraction of centerlines of detected cracks as well as the construction of a so-called signed distance function (SDF). An SDF, denoted byφ(x,y), defines the distance values from (x,y) to the evolving zero level-set Γ(t)={ϕ(x,y,t)=0}, which are positive inside of Γ(t) and negative outside of Γ(t), and satisfies |∇φ|=1 everywhere. The steady state of φ(x,y) is obtained through solving(14)∂φ∂τ+sgn(ϕ)(|∇φ|-1)=0,
where sgn(ϕ)=2H(ϕ)-1 is a signum function, and τ is another pseudotime variable added in the dependents of φ. With the initial condition φ(x,y,0)=ϕ(x,y,t), the zero level set Γ(t) will be preserved during evolving φ(x,y,τ). For centerline extraction, although a level set-based extraction of centerlines is possible [33], we employ a simple morphological skeleton operator that can preserve the connectivity of the original topographical patterns [32]. Theoretically, for objects with general topographical shapes, morphological skeltonization does not usually produce the ideal centerlines; nonetheless, for narrow, elongated 2D objects, such as cracks, morphological skeltonization is sufficiently accurate.The proposed width extraction method is summarized in the following:(1)
fromHθ(ϕe)>0.5, a binary image b(x,y) is obtained, wherein the detected cracks are marked by unity;(2)
perform skeletonization onb(x,y), and a new binary image c(x,y) is produced wherein the centerlines are unity;(3)
solve for an SDFψ(x,y) using (14) with ψ(x,y,0)=sgn(ϕ)=2Hθ(ϕe)-1;(4)
the shortest distances from the centerlines to the zero level set are obtained by a simple Boolean operation, which is denoted byd(x,y)=ψ(x,y), forall(x,y)∈{(x,y)∣c(x,y)=1}, the widths along cracks are approximated by w(x,y)≃2d(x,y).In this approximate method, there is no restriction to the complexity of topographical patterns of cracks as well as no need to manually specify the end points of cracks. Regarding the accuracy of this method, the approximate widths will be slightly greater than the exact values for a simple crack similar to Figure5.


## 5.3.1. Area and Perimeter of Cracks

In the Chan-Vese model, the computation of areas and perimeters of detected objects has been embedded in (10)–(12): the area A=∫ΩH(ϕe)dxdy, and the perimeter P=∫Ω|∇H(ϕe)|dxdy. However, for the purpose of computing areas and perimeters, a more “smeared-out” approximation will offer more numerical accuracy, which is(13)Hθ(ϕ)={1ifϕ>θ,0ifϕ<-θ,12[1+ϕθ+1πsin(πϕθ)]if|ϕ|≤θ,
where the parameter θ is selected comparable to the used space step. In the work of Osher and Sethian [31], θ=1.5h is recommended.


## 5.3.2. Extraction of Crack Width

Another important geometric quantity of engineering interest is crack width. However, unlike the computation of areas or perimeters for concrete cracks, the extraction of crack width is a nontrivial problem. First, the concept ofwidth is empirical. Recall that we adopt a nonrigorous way to describe cracks, namely, cracks in images are spatially narrow, elongated objects with strong contrast with the background, where the qualitative terms “narrow” or “elongated” are related to “width” without explicit definitions. In engineering practice, first, one needs to manually identify two end points for a single crack or multiple end points for a crack with a complex pattern; then by starting any of these end points, widths can be measured along the cracks continuously. To compute the crack widths analytically, a rigorous definition for crack width is necessary.Considering a simple case (Figure5), where a crack, with its closed boundary denoted by Γ, has two end points A:(xa,ya) and B:(xb,yb), therefore, the boundary of this crack can be split into two open boundaries, Γ1 and Γ2, which both start at A and end at B. The width w at a pixel (x1,y1) along one side of the crack boundaries, say, Γ1, is the shortest Euclidean distance (L2 distance) from (x1,y1) to a point (x2,y2) in the other side of boundaries Γ2. For cracks with complex patterns, for example, a continuous crack with “T”, “+” shapes, or other multiple angle splits, more than two end points have to be specified. One can break these cracks into several segments, wherein each segment has two end points, hence two sides of boundaries. The Euclidean distance-based definition for widths associated with these crack segments can be defined similarly.Figure 5
Definition of crack width: the distance from(x1,y1) to (x2,y2) is the shortest compared with others, for example, from (x1,y1) to (x2′,y2′); the former defines the crack width w at (x1,y1).


## 5.3.3. Automated Extraction Method

The preceding definition of crack width may give rise to an exact solution to compute widths with known end points for a given crack. However, as the topographical patterns of cracks become more complex, the manual intervention for specifying the possible end points will be more tedious. Herein, we propose an approximate method for width extraction. This approximate method relies on the extraction of centerlines of detected cracks as well as the construction of a so-called signed distance function (SDF). An SDF, denoted byφ(x,y), defines the distance values from (x,y) to the evolving zero level-set Γ(t)={ϕ(x,y,t)=0}, which are positive inside of Γ(t) and negative outside of Γ(t), and satisfies |∇φ|=1 everywhere. The steady state of φ(x,y) is obtained through solving(14)∂φ∂τ+sgn(ϕ)(|∇φ|-1)=0,
where sgn(ϕ)=2H(ϕ)-1 is a signum function, and τ is another pseudotime variable added in the dependents of φ. With the initial condition φ(x,y,0)=ϕ(x,y,t), the zero level set Γ(t) will be preserved during evolving φ(x,y,τ). For centerline extraction, although a level set-based extraction of centerlines is possible [33], we employ a simple morphological skeleton operator that can preserve the connectivity of the original topographical patterns [32]. Theoretically, for objects with general topographical shapes, morphological skeltonization does not usually produce the ideal centerlines; nonetheless, for narrow, elongated 2D objects, such as cracks, morphological skeltonization is sufficiently accurate.The proposed width extraction method is summarized in the following:(1)
fromHθ(ϕe)>0.5, a binary image b(x,y) is obtained, wherein the detected cracks are marked by unity;(2)
perform skeletonization onb(x,y), and a new binary image c(x,y) is produced wherein the centerlines are unity;(3)
solve for an SDFψ(x,y) using (14) with ψ(x,y,0)=sgn(ϕ)=2Hθ(ϕe)-1;(4)
the shortest distances from the centerlines to the zero level set are obtained by a simple Boolean operation, which is denoted byd(x,y)=ψ(x,y), forall(x,y)∈{(x,y)∣c(x,y)=1}, the widths along cracks are approximated by w(x,y)≃2d(x,y).In this approximate method, there is no restriction to the complexity of topographical patterns of cracks as well as no need to manually specify the end points of cracks. Regarding the accuracy of this method, the approximate widths will be slightly greater than the exact values for a simple crack similar to Figure5.


## 6. Application I: Field Image

In routine field inspection of concrete structures, a large amount of digital images are usually collected. Well-organized image archives as well as text-based documents provide an informative means of reporting the structural integrity of civil infrastructure objects. However, these field images are mostly used as a visual media only—readers or analysts read the images visually and may subjectively correlate visual identification with other engineering knowledge. To obtain more quantitative information, an effective image analysis method is needed. In what follows, we present a comparative study to demonstrate that the proposed framework provides a quantitative image-mining solution to this problem, which is unparalleled compared to other non-contour-based image analysis methods.Figure6 displays a field image of the surface of a concrete region that has clearly observed cracks. One first notices that the surface textures vary spatially showing different modalities (e.g., lighter versus darker texture). Second, there are a large number of small air pockets on the surface. Finally, cracks are not fully connected; rather discontinuities are observed between crack segments. The presence of these artifacts results in significant adverse effects and prevents from realizing an automatic identification of concrete cracks.Figure 6
Example image taken from a concrete component in the field. Note the non-crack image artifacts annotated. The subregion in the bottom-right bounding box is used to illustrate the proposed level-set detection approach.We first attempt to use a traditional approach. In Figure7, two crack detection results are shown using the classical Canny edge detector [34]. The two results are presented using different analysis parameters. It is seen that the two results are dramatically different, and it is very difficult to differentiate the boundaries of the concrete cracks from the noisy edges in the background. Although some statistical methods can mitigate this parameter selection problem (e.g., [35]), the critical difficulty that prohibits subsequent geometric quantification is that the detected edges are merely nonclosed segments rather than closed contours. This difficulty exists similarly if other non-model-based edge detectors (e.g., a Sobel detector or wavelet transform-based detectors) are used.Crack detection using the Canny edge detector (as implemented in the Matlab Image Processing Toobox, [8]), with different parameters: (a) threshold parameter = 0.1, Gaussian smoothing parameter = 1.0; (2) threshold parameter = 0.2, Gaussian smoothing parameter = 1.5.
(a)(b)The proposed level-set active contour approach can avoid this difficulty. Before applying this approach, we note that the purpose of crack identification in this situation is to obtain subsequent quantitative characteristics of the cracks rather than to solely determine the presence of the cracks (i.e., the analyst knows that there are cracks in the image). Therefore, we assume that the analyst interactively selects a region of interest (ROI) in the image domain then seeks a quantitative characterization. In Figure8(a), such an ROI is shown, which contains multiple intersecting cracks. Starting with the initialization in Figure 8(b), we apply the level-set based detection, which iteratively converges to the final detected image contours as shown in Figures 8(c), 8(d), and 8(e). In the final result of Figure 8(e), one observes that small contours that come from air pockets are detected as well. These small non-crack boundaries can be easily removed by a labeling and area thresholding process, and the final crack contours are shown in Figure 8(f).Crack detection using the level set-based active contour method: (a) the selected ROI, also shown in Figure6; (b) the initialized contours; (c) and (d) intermediate detection results; (e) the converged result; (f) the cleaned detection result with contours for the air pockets removed.
(a)(b)(c)(d)(e)(f)Different from the detection results shown in Figure7, each detected crack segment in Figure 8(f) is featured with a closed contour (and an underlying mathematical level-set function). This property determines that the quantitative extraction of crack width and other geometric properties can be readily produced using the proposed methods in the previous section. We demonstrate this application thoroughly in the next section.


## 7. Application II: Laboratory Images

To demonstrate a complete process of concrete surface crack monitoring and damage quantification using the proposed framework, a laboratory beam is loaded and monitored using an array of cameras. As shown in Figure9(a), the experiment involves a simply supported concrete beam loaded at mid-span controlled with linearly increasing displacements. A 2-point arrangement of loading is provided resulting in a region of constant moment. The beam span is 60inches (1524mm), and its section size is 5inches×7inches (127mm×177.8mm). It was designed to fail in a shear mode using 4000psi (27.6MPa) concrete with no shear-resisting stirrups. The front view of the beam was monitored during loading by four digital cameras. The resolution of the image frames in the recorded video streams is 0.89mm/pixel, and the frequency of camera capture was set at 40frames/sec.(a) Schematic of test set-up for a simple supported beam captured by four cameras (C1–C4 as seen in (b)); the fields of views (FOV1–FOV4) of the four cameras are overlapped and symmetric with the centerline; (b) the four cameras (C1–C4) in the real set-up; (c) the captured scene before loading by the camera C1 and the region of interest.
(a)(b)(c)To illustrate the use of image-based crack damage monitoring and quantification, we select a short sequence of image frames generated from Camera C1, which were continuously captured within a1.725-second duration, resulting in 70 image frames. In Figure 10, we illustrate six frames that are representative of the temporal development of crack occurrence. Note that the images shown here are cropped from the raw image frames, resulting in a smaller image domain, which is 240×268.Concrete surface at different acquisition times: (a) at 0 second (the first frame with dividing grids of the image domain); (b) at 0.625 second (the26th frame); (c) at 0.65 second (the 27th frame at which the first crack appeared); (d) at 0.80 second (the 33rd frame); (e) at 1.20 seconds (the 45th frame); (f) at 1.85 second (the 70th frame). In (a), the 10×10 subdivisions for the image domain are overlaid with the image.
(a)(b)(c)(d)(e)(f)


### 7.1. Crack Damage Monitoring

#### 7.1.1. Computation of Time-Series Damage Measures

In Figure10(a), we show that the image domain is divided into 10×10 subdomains; each of the subdomains has 24×26 pixels except the subdomains in the far right column which are of size 24×34. For each subdomain, we compute the time-varying MD measurements (as well as the motion estimates). Two methods of computing MD measurements are conducted in parallel. First, we compute the MD measurements using temporally adjacent image patches. In Figure 11(a), MD measurements resulting from this method are plotted against time. For completeness, we illustrate the extracted progressive motion quantities across time in Figure 11(b), which are the average motion quantities based on the center 2×2 image patches. The second method of computing MD measurements, which uses a fixed image frame at a certain time as the base, will not be started until crack damage is detected.(a) Manifold-distance measurements across time; (b) the corresponding motion of the centerline of the beam.
(a)(b)From Figure11(a), it can been seen that all MD measurements computed before t=0.650seconds from the 100 image patches are in a narrow band (about 0.1~0.2). At t=0.650seconds, peak values of the MD measurements from some subdomains emerge, which are a strong indication of crack damage occurrence. The evolution of extracted motion quantities in the center of the beam agrees with the time of this crack event (Figure 10(c)). However, we see that these motion quantities, if used as features of crack occurrence, are not as discriminative as the MD measures across time.


### 7.2. Outlier Analysis and Initial Damage Localization

The above conclusion regarding when the crack event occurs is based on subjective observation. In a practical structural damage monitoring context, an outlier analysis may be employed to automate this process. An outlier analysis, in its simplest form, is implemented by assuming that the underlying features follow a normal distribution (for more rigorous treatment about modeling the distribution of outliers and other approaches, one may refer to [25]). Herein, we take advantage of the fact that there is no crack in the second image frame (at 0.025seconds), therefore the first group of MD measures at this moment are used to extract the sample mean and the sample standard deviation of MD measures, denoted by μm and σm. At the next image frame (at 0.050seconds), the MD measures MD[i,j] (1≤i,j≤10) are used to test the hypothesis that if Ho:|MD[i,j]-μm|<λ,σm holds for a fixed threshold λ. Any MD[i,j]s that fail the test Ho are treated as outliers indicating crack damage; otherwise, they are normal MD measures, and are used to update the sample mean μm and the sample standard deviation σm. In addition, as the hypothesis test goes, at each time, a binary damage map for the divided image domain is generated progressively, where the subdomain is set unity if it is determined as a damaged region. In Figure 12, this initial damage localization is illustrated.Distance measurements before 0.650 seconds: (a) manifold distance; (b) crude damage localization att = 0.650 seconds, where the axis numbers mark the location of different subdomains shown in Figure 10(a).
(a)(b)


### 7.3. Crack Damage Boundary Detection and Quantification

In Figure13, images showing the cracked regions at later times are shown. These are automatically cropped based on the initial localization information presented in Figure 12(b). In the following numerical studies, we generally choose the parameters as follows: space step h=1, time step Δt=0.1 for (10) and Δt=0.25 for (14), and ϵ=h=1 and θ=1.5h=1.5 for the approximate Heaviside functions. The length scale μ in principle should vary relative to the scale of the target objects, and empirical investigation needs to be conducted to determine the suitable range of this parameter. In this paper, we use the average of these bounds, approximately μ=0.008·2552, throughout this work.Cropped images (72×160 pixels) of concrete cracks at different times, from (a) to (c), denoted by T1,T2, and T3. By synchronizing with the loading time, T1=30.5seconds, T2=34.75seconds, and T3=39.75seconds (these time stamps are different from the relative ones used in Figure 10.
(a)(b)(c)


#### 7.3.1. Damage Detection for Regions of Interest

Recall that the Chan-Vese active model in this paper is derived for two-phase images, namely, assuming that the images contain only two objects, the normal background and the potential structural damage. This limitation in practice should be considered when selecting the regions of interest in structural components. Thus, we only consider a small area in the left shear zone of the beam as shown in Figure9(c). Also, since a large number of image frames were recorded, we subsampled and selected three representative image frames in this work. Figures 13(a)–13(c) show the three images captured at different times for the same target region as shown in Figure 9(c). In Figures 14(a)–14(c), the boundaries of the cracks are successfully detected by applying the Chan-Vese active contour model. It is notable that in Figure 14(a), an interior noncracked area appears, which visually can be seen in Figure 13(a).(a)–(c) Boundaries of detected cracks at different times; (d)–(f) extracted centerlines of detected cracks by morphological skeltonization.
(a)
T1(b)
T2(c)
T3(d)
T1(e)
T2(f)
T3


#### 7.3.2. Extraction of Geometric Quantities: Areas and Perimeters

By visually inspecting the detected cracks in Figures14(a)–14(c), one can observe that their topographical patterns change significantly with time (T1→T3) due to increasing mid-span displacements. We first compute the areas and perimeters of these cracks using the methods discussed previously. In Table 1, these results are reported as well as their corresponding true values considering the image resolution. The true values are determined by manually summing the pixels. It is interesting to note that the areas of the cracks increase; however, the perimeters of the cracks with increasing mid-span displacements do not necessarily increase. This is due to the higher topographical variations of concrete cracks at lower displacement levels.Table 1
Extracted areas and perimeters of cracks (numbers in parentheses are true values of extracted quantities: for areas, the multiplier is0.892mm2; for perimeter, the multiplier is 0.89mm).
Geometric quantitiesCrack images atT1T2T3Areas79312401698(628mm2)(982mm2)(1345mm2)Perimeters308279289(274mm)(248mm)(257mm)


#### 7.3.3. Extraction of Crack Widths

In accordance with the approximate procedure for extraction of crack width, we first perform skeltonization over the binary detection results at timeT1,T2, and T3. Figures 14(d)–14(f) illustrate the centerlines of these multitemporal cracks. By solving (14) for the SDFs and conducting Boolean operations to extract distance values at the centerlines, one obtains the approximate widths along cracks, denoted by wT1(x,y), wT2(x,y), and wT3(x,y). Table 2 summarizes the statistics of these crack widths. Since it is difficult to plot extracted crack widths along the crack in a 2D domain, for width values in w(x,y), we plot the pairs (x,maxy[w(x,y)]), that is, the x-coordinate versus the maximum width in y direction given x. These plots are shown in Figures 15(a)–15(c).Table 2
Crack width statistics for detected cracks at different times: the top two rows are from the level set-based approximate method, and the bottom two rows are from the explicit searching method (numbers in parentheses are true values of crack width considering the image resolution).
Crack width statisticsCrack images atT1T2T3ApproximateMaximum11.8 (10.4mm)11.9  (10.6mm)16.7 (14.8mm)Average4.3 (3.8mm)7.5 (6.7mm)9.9 (8.8mm)ExactMaximum12.6 (11.2mm)11.1 (9.9mm)16.0 (14.2mm)Average4.3 (3.8mm)7.3 (6.5mm)9.8 (8.7mm)Plots of the extracted crack widths at timesT1,T2, and T3 (a)–(c) using the proposed level set-based approximate method and (d)–(f) using the explicit exact method.
(a)
T1(b)
T2(c)
T3(d)
T1(e)
T2(f)
T3To evaluate the accuracy of the above approximate crack widths, we use an explicit searching method to extract the exact crack widths: first, we explicitly start traveling along one side of the two boundaries (e.g., theΓ1 in Figure 5) from one end point, which can be realized by searching the 8 connected neighborhoods repeatedly; then at each pixel (x,y), the shortest distance to the other side of the two boundaries (Γ2) is computed by minimizing the distance values from (x,y) to all the pixels in Γ2. This explicit searching method yields the exact values of crack width along a simple crack. Certainly this method suffers from the limitations that manual intervention is needed to provide the end points, and further generalization to complex crack patterns is difficult. Nonetheless, it provides the baseline for comparison with the approximate method. The bottom two rows of Table 2 summarize the width statistics of the obtained exact crack widths, and in Figures 15(d)–15(f) the crack widths are again plotted along the x-axis.By comparing the statistics results from the two extraction methods in Table2, one can see that the basic statistics of extracted crack widths at the three different times are considerably close. One can further compare the plots of extracted width along x-axis in Figure 15; visually, at each time of T1,T2, or T3, the obtained widths from the two methods agree very well. Quantitatively, we compute the root-mean-square-error (RMSE) over the results from the exact method and from the approximate method: at T1, RMSE=2.2mm at T2, RMSE=0.7mm and at T3, RMSE=1.1mm. At time T1, slightly higher error is produced because in the exact method, the small interior noncracked region within the detected crack is ignored, therefore, higher widths are yielded around this region while in using the level set-based approximate method, the interior noncracked region is automatically excluded in computing the width.


#### 7.3.4. Correlation with Displacement and Load Amplitudes

In this laboratory experiment, the linearly increasing displacements at the mid-span and associated load amplitudes from the actuator were recorded at40Hz. Based on the synchronized computer times of the image-acquisition workstation and the actuator-controller workstation, each captured image frame can be mapped to a simultaneous mid-span displacement and load amplitude. In Figure 16(a), the linearly-increased mid-span displacement versus time is plotted, with the times when the crack images used previously were captured articulated. In Figure 16(b), the load-displacement curve is provided, wherein there simultaneous values of average crack width at times T1-T3 are annotated.(a) Displacements versus time at mid-span; (b) measured load amplitude versus displacement at the mid-span, and the extracted average crack widths atT1,T2, and T3.
(a)(b)These correlation plots illustrate that cracks in a critical region of a structural component can be captured continuously with a proper camera array setup, and if available, this data can be synchronized with other measurements. It is noted that this continuous sampling of realistic 3D objects, as made possible with the 2D images, is not possible with any other current measurement methods that are discontinuous in space. With an offline analysis proposed in this paper, the positions of cracks and the associated geometric quantities can be tracked continuously in time. Hence, one can correlate the crack characteristics to other quantities of interest, such as load amplitude, displacement, or other measured response data.


## 7.1. Crack Damage Monitoring

### 7.1.1. Computation of Time-Series Damage Measures

In Figure10(a), we show that the image domain is divided into 10×10 subdomains; each of the subdomains has 24×26 pixels except the subdomains in the far right column which are of size 24×34. For each subdomain, we compute the time-varying MD measurements (as well as the motion estimates). Two methods of computing MD measurements are conducted in parallel. First, we compute the MD measurements using temporally adjacent image patches. In Figure 11(a), MD measurements resulting from this method are plotted against time. For completeness, we illustrate the extracted progressive motion quantities across time in Figure 11(b), which are the average motion quantities based on the center 2×2 image patches. The second method of computing MD measurements, which uses a fixed image frame at a certain time as the base, will not be started until crack damage is detected.(a) Manifold-distance measurements across time; (b) the corresponding motion of the centerline of the beam.
(a)(b)From Figure11(a), it can been seen that all MD measurements computed before t=0.650seconds from the 100 image patches are in a narrow band (about 0.1~0.2). At t=0.650seconds, peak values of the MD measurements from some subdomains emerge, which are a strong indication of crack damage occurrence. The evolution of extracted motion quantities in the center of the beam agrees with the time of this crack event (Figure 10(c)). However, we see that these motion quantities, if used as features of crack occurrence, are not as discriminative as the MD measures across time.


## 7.1.1. Computation of Time-Series Damage Measures

In Figure10(a), we show that the image domain is divided into 10×10 subdomains; each of the subdomains has 24×26 pixels except the subdomains in the far right column which are of size 24×34. For each subdomain, we compute the time-varying MD measurements (as well as the motion estimates). Two methods of computing MD measurements are conducted in parallel. First, we compute the MD measurements using temporally adjacent image patches. In Figure 11(a), MD measurements resulting from this method are plotted against time. For completeness, we illustrate the extracted progressive motion quantities across time in Figure 11(b), which are the average motion quantities based on the center 2×2 image patches. The second method of computing MD measurements, which uses a fixed image frame at a certain time as the base, will not be started until crack damage is detected.(a) Manifold-distance measurements across time; (b) the corresponding motion of the centerline of the beam.
(a)(b)From Figure11(a), it can been seen that all MD measurements computed before t=0.650seconds from the 100 image patches are in a narrow band (about 0.1~0.2). At t=0.650seconds, peak values of the MD measurements from some subdomains emerge, which are a strong indication of crack damage occurrence. The evolution of extracted motion quantities in the center of the beam agrees with the time of this crack event (Figure 10(c)). However, we see that these motion quantities, if used as features of crack occurrence, are not as discriminative as the MD measures across time.


## 7.2. Outlier Analysis and Initial Damage Localization

The above conclusion regarding when the crack event occurs is based on subjective observation. In a practical structural damage monitoring context, an outlier analysis may be employed to automate this process. An outlier analysis, in its simplest form, is implemented by assuming that the underlying features follow a normal distribution (for more rigorous treatment about modeling the distribution of outliers and other approaches, one may refer to [25]). Herein, we take advantage of the fact that there is no crack in the second image frame (at 0.025seconds), therefore the first group of MD measures at this moment are used to extract the sample mean and the sample standard deviation of MD measures, denoted by μm and σm. At the next image frame (at 0.050seconds), the MD measures MD[i,j] (1≤i,j≤10) are used to test the hypothesis that if Ho:|MD[i,j]-μm|<λ,σm holds for a fixed threshold λ. Any MD[i,j]s that fail the test Ho are treated as outliers indicating crack damage; otherwise, they are normal MD measures, and are used to update the sample mean μm and the sample standard deviation σm. In addition, as the hypothesis test goes, at each time, a binary damage map for the divided image domain is generated progressively, where the subdomain is set unity if it is determined as a damaged region. In Figure 12, this initial damage localization is illustrated.Distance measurements before 0.650 seconds: (a) manifold distance; (b) crude damage localization att = 0.650 seconds, where the axis numbers mark the location of different subdomains shown in Figure 10(a).
(a)(b)


## 7.3. Crack Damage Boundary Detection and Quantification

In Figure13, images showing the cracked regions at later times are shown. These are automatically cropped based on the initial localization information presented in Figure 12(b). In the following numerical studies, we generally choose the parameters as follows: space step h=1, time step Δt=0.1 for (10) and Δt=0.25 for (14), and ϵ=h=1 and θ=1.5h=1.5 for the approximate Heaviside functions. The length scale μ in principle should vary relative to the scale of the target objects, and empirical investigation needs to be conducted to determine the suitable range of this parameter. In this paper, we use the average of these bounds, approximately μ=0.008·2552, throughout this work.Cropped images (72×160 pixels) of concrete cracks at different times, from (a) to (c), denoted by T1,T2, and T3. By synchronizing with the loading time, T1=30.5seconds, T2=34.75seconds, and T3=39.75seconds (these time stamps are different from the relative ones used in Figure 10.
(a)(b)(c)


### 7.3.1. Damage Detection for Regions of Interest

Recall that the Chan-Vese active model in this paper is derived for two-phase images, namely, assuming that the images contain only two objects, the normal background and the potential structural damage. This limitation in practice should be considered when selecting the regions of interest in structural components. Thus, we only consider a small area in the left shear zone of the beam as shown in Figure9(c). Also, since a large number of image frames were recorded, we subsampled and selected three representative image frames in this work. Figures 13(a)–13(c) show the three images captured at different times for the same target region as shown in Figure 9(c). In Figures 14(a)–14(c), the boundaries of the cracks are successfully detected by applying the Chan-Vese active contour model. It is notable that in Figure 14(a), an interior noncracked area appears, which visually can be seen in Figure 13(a).(a)–(c) Boundaries of detected cracks at different times; (d)–(f) extracted centerlines of detected cracks by morphological skeltonization.
(a)
T1(b)
T2(c)
T3(d)
T1(e)
T2(f)
T3


### 7.3.2. Extraction of Geometric Quantities: Areas and Perimeters

By visually inspecting the detected cracks in Figures14(a)–14(c), one can observe that their topographical patterns change significantly with time (T1→T3) due to increasing mid-span displacements. We first compute the areas and perimeters of these cracks using the methods discussed previously. In Table 1, these results are reported as well as their corresponding true values considering the image resolution. The true values are determined by manually summing the pixels. It is interesting to note that the areas of the cracks increase; however, the perimeters of the cracks with increasing mid-span displacements do not necessarily increase. This is due to the higher topographical variations of concrete cracks at lower displacement levels.Table 1
Extracted areas and perimeters of cracks (numbers in parentheses are true values of extracted quantities: for areas, the multiplier is0.892mm2; for perimeter, the multiplier is 0.89mm).
Geometric quantitiesCrack images atT1T2T3Areas79312401698(628mm2)(982mm2)(1345mm2)Perimeters308279289(274mm)(248mm)(257mm)


### 7.3.3. Extraction of Crack Widths

In accordance with the approximate procedure for extraction of crack width, we first perform skeltonization over the binary detection results at timeT1,T2, and T3. Figures 14(d)–14(f) illustrate the centerlines of these multitemporal cracks. By solving (14) for the SDFs and conducting Boolean operations to extract distance values at the centerlines, one obtains the approximate widths along cracks, denoted by wT1(x,y), wT2(x,y), and wT3(x,y). Table 2 summarizes the statistics of these crack widths. Since it is difficult to plot extracted crack widths along the crack in a 2D domain, for width values in w(x,y), we plot the pairs (x,maxy[w(x,y)]), that is, the x-coordinate versus the maximum width in y direction given x. These plots are shown in Figures 15(a)–15(c).Table 2
Crack width statistics for detected cracks at different times: the top two rows are from the level set-based approximate method, and the bottom two rows are from the explicit searching method (numbers in parentheses are true values of crack width considering the image resolution).
Crack width statisticsCrack images atT1T2T3ApproximateMaximum11.8 (10.4mm)11.9  (10.6mm)16.7 (14.8mm)Average4.3 (3.8mm)7.5 (6.7mm)9.9 (8.8mm)ExactMaximum12.6 (11.2mm)11.1 (9.9mm)16.0 (14.2mm)Average4.3 (3.8mm)7.3 (6.5mm)9.8 (8.7mm)Plots of the extracted crack widths at timesT1,T2, and T3 (a)–(c) using the proposed level set-based approximate method and (d)–(f) using the explicit exact method.
(a)
T1(b)
T2(c)
T3(d)
T1(e)
T2(f)
T3To evaluate the accuracy of the above approximate crack widths, we use an explicit searching method to extract the exact crack widths: first, we explicitly start traveling along one side of the two boundaries (e.g., theΓ1 in Figure 5) from one end point, which can be realized by searching the 8 connected neighborhoods repeatedly; then at each pixel (x,y), the shortest distance to the other side of the two boundaries (Γ2) is computed by minimizing the distance values from (x,y) to all the pixels in Γ2. This explicit searching method yields the exact values of crack width along a simple crack. Certainly this method suffers from the limitations that manual intervention is needed to provide the end points, and further generalization to complex crack patterns is difficult. Nonetheless, it provides the baseline for comparison with the approximate method. The bottom two rows of Table 2 summarize the width statistics of the obtained exact crack widths, and in Figures 15(d)–15(f) the crack widths are again plotted along the x-axis.By comparing the statistics results from the two extraction methods in Table2, one can see that the basic statistics of extracted crack widths at the three different times are considerably close. One can further compare the plots of extracted width along x-axis in Figure 15; visually, at each time of T1,T2, or T3, the obtained widths from the two methods agree very well. Quantitatively, we compute the root-mean-square-error (RMSE) over the results from the exact method and from the approximate method: at T1, RMSE=2.2mm at T2, RMSE=0.7mm and at T3, RMSE=1.1mm. At time T1, slightly higher error is produced because in the exact method, the small interior noncracked region within the detected crack is ignored, therefore, higher widths are yielded around this region while in using the level set-based approximate method, the interior noncracked region is automatically excluded in computing the width.


### 7.3.4. Correlation with Displacement and Load Amplitudes

In this laboratory experiment, the linearly increasing displacements at the mid-span and associated load amplitudes from the actuator were recorded at40Hz. Based on the synchronized computer times of the image-acquisition workstation and the actuator-controller workstation, each captured image frame can be mapped to a simultaneous mid-span displacement and load amplitude. In Figure 16(a), the linearly-increased mid-span displacement versus time is plotted, with the times when the crack images used previously were captured articulated. In Figure 16(b), the load-displacement curve is provided, wherein there simultaneous values of average crack width at times T1-T3 are annotated.(a) Displacements versus time at mid-span; (b) measured load amplitude versus displacement at the mid-span, and the extracted average crack widths atT1,T2, and T3.
(a)(b)These correlation plots illustrate that cracks in a critical region of a structural component can be captured continuously with a proper camera array setup, and if available, this data can be synchronized with other measurements. It is noted that this continuous sampling of realistic 3D objects, as made possible with the 2D images, is not possible with any other current measurement methods that are discontinuous in space. With an offline analysis proposed in this paper, the positions of cracks and the associated geometric quantities can be tracked continuously in time. Hence, one can correlate the crack characteristics to other quantities of interest, such as load amplitude, displacement, or other measured response data.


## 7.3.1. Damage Detection for Regions of Interest

Recall that the Chan-Vese active model in this paper is derived for two-phase images, namely, assuming that the images contain only two objects, the normal background and the potential structural damage. This limitation in practice should be considered when selecting the regions of interest in structural components. Thus, we only consider a small area in the left shear zone of the beam as shown in Figure9(c). Also, since a large number of image frames were recorded, we subsampled and selected three representative image frames in this work. Figures 13(a)–13(c) show the three images captured at different times for the same target region as shown in Figure 9(c). In Figures 14(a)–14(c), the boundaries of the cracks are successfully detected by applying the Chan-Vese active contour model. It is notable that in Figure 14(a), an interior noncracked area appears, which visually can be seen in Figure 13(a).(a)–(c) Boundaries of detected cracks at different times; (d)–(f) extracted centerlines of detected cracks by morphological skeltonization.
(a)
T1(b)
T2(c)
T3(d)
T1(e)
T2(f)
T3


## 7.3.2. Extraction of Geometric Quantities: Areas and Perimeters

By visually inspecting the detected cracks in Figures14(a)–14(c), one can observe that their topographical patterns change significantly with time (T1→T3) due to increasing mid-span displacements. We first compute the areas and perimeters of these cracks using the methods discussed previously. In Table 1, these results are reported as well as their corresponding true values considering the image resolution. The true values are determined by manually summing the pixels. It is interesting to note that the areas of the cracks increase; however, the perimeters of the cracks with increasing mid-span displacements do not necessarily increase. This is due to the higher topographical variations of concrete cracks at lower displacement levels.Table 1
Extracted areas and perimeters of cracks (numbers in parentheses are true values of extracted quantities: for areas, the multiplier is0.892mm2; for perimeter, the multiplier is 0.89mm).
Geometric quantitiesCrack images atT1T2T3Areas79312401698(628mm2)(982mm2)(1345mm2)Perimeters308279289(274mm)(248mm)(257mm)


## 7.3.3. Extraction of Crack Widths

In accordance with the approximate procedure for extraction of crack width, we first perform skeltonization over the binary detection results at timeT1,T2, and T3. Figures 14(d)–14(f) illustrate the centerlines of these multitemporal cracks. By solving (14) for the SDFs and conducting Boolean operations to extract distance values at the centerlines, one obtains the approximate widths along cracks, denoted by wT1(x,y), wT2(x,y), and wT3(x,y). Table 2 summarizes the statistics of these crack widths. Since it is difficult to plot extracted crack widths along the crack in a 2D domain, for width values in w(x,y), we plot the pairs (x,maxy[w(x,y)]), that is, the x-coordinate versus the maximum width in y direction given x. These plots are shown in Figures 15(a)–15(c).Table 2
Crack width statistics for detected cracks at different times: the top two rows are from the level set-based approximate method, and the bottom two rows are from the explicit searching method (numbers in parentheses are true values of crack width considering the image resolution).
Crack width statisticsCrack images atT1T2T3ApproximateMaximum11.8 (10.4mm)11.9  (10.6mm)16.7 (14.8mm)Average4.3 (3.8mm)7.5 (6.7mm)9.9 (8.8mm)ExactMaximum12.6 (11.2mm)11.1 (9.9mm)16.0 (14.2mm)Average4.3 (3.8mm)7.3 (6.5mm)9.8 (8.7mm)Plots of the extracted crack widths at timesT1,T2, and T3 (a)–(c) using the proposed level set-based approximate method and (d)–(f) using the explicit exact method.
(a)
T1(b)
T2(c)
T3(d)
T1(e)
T2(f)
T3To evaluate the accuracy of the above approximate crack widths, we use an explicit searching method to extract the exact crack widths: first, we explicitly start traveling along one side of the two boundaries (e.g., theΓ1 in Figure 5) from one end point, which can be realized by searching the 8 connected neighborhoods repeatedly; then at each pixel (x,y), the shortest distance to the other side of the two boundaries (Γ2) is computed by minimizing the distance values from (x,y) to all the pixels in Γ2. This explicit searching method yields the exact values of crack width along a simple crack. Certainly this method suffers from the limitations that manual intervention is needed to provide the end points, and further generalization to complex crack patterns is difficult. Nonetheless, it provides the baseline for comparison with the approximate method. The bottom two rows of Table 2 summarize the width statistics of the obtained exact crack widths, and in Figures 15(d)–15(f) the crack widths are again plotted along the x-axis.By comparing the statistics results from the two extraction methods in Table2, one can see that the basic statistics of extracted crack widths at the three different times are considerably close. One can further compare the plots of extracted width along x-axis in Figure 15; visually, at each time of T1,T2, or T3, the obtained widths from the two methods agree very well. Quantitatively, we compute the root-mean-square-error (RMSE) over the results from the exact method and from the approximate method: at T1, RMSE=2.2mm at T2, RMSE=0.7mm and at T3, RMSE=1.1mm. At time T1, slightly higher error is produced because in the exact method, the small interior noncracked region within the detected crack is ignored, therefore, higher widths are yielded around this region while in using the level set-based approximate method, the interior noncracked region is automatically excluded in computing the width.


## 7.3.4. Correlation with Displacement and Load Amplitudes

In this laboratory experiment, the linearly increasing displacements at the mid-span and associated load amplitudes from the actuator were recorded at40Hz. Based on the synchronized computer times of the image-acquisition workstation and the actuator-controller workstation, each captured image frame can be mapped to a simultaneous mid-span displacement and load amplitude. In Figure 16(a), the linearly-increased mid-span displacement versus time is plotted, with the times when the crack images used previously were captured articulated. In Figure 16(b), the load-displacement curve is provided, wherein there simultaneous values of average crack width at times T1-T3 are annotated.(a) Displacements versus time at mid-span; (b) measured load amplitude versus displacement at the mid-span, and the extracted average crack widths atT1,T2, and T3.
(a)(b)These correlation plots illustrate that cracks in a critical region of a structural component can be captured continuously with a proper camera array setup, and if available, this data can be synchronized with other measurements. It is noted that this continuous sampling of realistic 3D objects, as made possible with the 2D images, is not possible with any other current measurement methods that are discontinuous in space. With an offline analysis proposed in this paper, the positions of cracks and the associated geometric quantities can be tracked continuously in time. Hence, one can correlate the crack characteristics to other quantities of interest, such as load amplitude, displacement, or other measured response data.


## 8. Concerns Prior to Practical Application

### 8.1. Desired Resolution

In the previous two sections, the proposed framework and the component algorithms are applied in two different practical situations. In either of these example applications, if engineering properties are to be extracted, the resolution of the images must be carefully considered. In the field, with a large amount of (potentially manually captured) images, it may not be feasible to reliably configure the camera capture hardware to accomplish consistent and reasonable resolution. One of the easiest solutions is to place a small pattern, for example, drawn onto paper, within the field of view. Using this reference pattern, a simple calculation based on the known size of the pattern and its size in the image domain may give an approximate estimate of the image resolution (e.g., in terms of length unit per pixel).The laboratory application example demonstrates a complete application of the proposed framework. In this experiment, the image resolution is determined to be 0.89 mm/pixel. It is noted that narrower, hairline cracks, which may have widths as small as 0.1 mm cannot be captured with this set-up. This is primarily a limitation of the hardware set-up. Alternative camera type and placement relative to the field of view of interest can resolve this issue. If the width of the cracks is still less than the image resolution, a detection method with subpixel accuracy may be necessary. The level-set active contour method as presented in this paper has not been calibrated to solve this problem, as it internally relies on the calculation of the mean intensity values within contoured image areas. Multiple overlapping fields of view have been used with success in other applications to achieve subpixel accuracy; however, this would be a subject of future study using the currently presented methods.


### 8.2. Image Noise and Size

A practical problem in any image-based analysis situation is that noisy image artifacts may exist in the domain. In the case of concrete surface images, these artifacts include large decolored, roughed, or spalled concrete surface areas. In this situation, prior to turning to other image segmentation or object recognition methods, we propose that a subdomain approach may be used first. Such an approach may be in the form of dividing the image domain into multiple subdomains or in the form of identifying user-selected ROIs. It is noted that the influence of noisy artifacts will appear in the first application example if a different ROI is selected (Figure6). Another practical issue, with a related solution, is that of image size. If images with very large sizes are encountered, computational cost may increase significantly; meanwhile robustness to noisy artifacts likely decreases. In this situation, the subdomain approach is an effective solution as well.


### 8.3. Computational Cost

The algorithms in the proposed framework, including the manifold-based monitoring algorithm and the level set-based crack detection and quantification algorithms, comprise a large number of basic vector or matrix operations with the complexity order atO(N) or O(N2). As a result, the associated computational cost is considerably larger than the traditional simple filtering-based methods. Computational cost is a concern if the framework is needed for online monitoring and quantification of rapidly captured image sequences. For offline processing of single images, however, the computation time is acceptable based on our implementation of the proposed algorithms in a Matlab environment [36]. If the framework is used as an online monitoring solution, fine tuning of the image capture rate and the image processing time is subject to further investigation. A plausible solution may be to subsample the captured images (i.e., only process a fraction of the images). It is noted that most of the time consumed from an operational point of view is spent tuning the internal parameters. This issue is expected to be mitigated when more empirical knowledge is obtained.


## 8.1. Desired Resolution

In the previous two sections, the proposed framework and the component algorithms are applied in two different practical situations. In either of these example applications, if engineering properties are to be extracted, the resolution of the images must be carefully considered. In the field, with a large amount of (potentially manually captured) images, it may not be feasible to reliably configure the camera capture hardware to accomplish consistent and reasonable resolution. One of the easiest solutions is to place a small pattern, for example, drawn onto paper, within the field of view. Using this reference pattern, a simple calculation based on the known size of the pattern and its size in the image domain may give an approximate estimate of the image resolution (e.g., in terms of length unit per pixel).The laboratory application example demonstrates a complete application of the proposed framework. In this experiment, the image resolution is determined to be 0.89 mm/pixel. It is noted that narrower, hairline cracks, which may have widths as small as 0.1 mm cannot be captured with this set-up. This is primarily a limitation of the hardware set-up. Alternative camera type and placement relative to the field of view of interest can resolve this issue. If the width of the cracks is still less than the image resolution, a detection method with subpixel accuracy may be necessary. The level-set active contour method as presented in this paper has not been calibrated to solve this problem, as it internally relies on the calculation of the mean intensity values within contoured image areas. Multiple overlapping fields of view have been used with success in other applications to achieve subpixel accuracy; however, this would be a subject of future study using the currently presented methods.


## 8.2. Image Noise and Size

A practical problem in any image-based analysis situation is that noisy image artifacts may exist in the domain. In the case of concrete surface images, these artifacts include large decolored, roughed, or spalled concrete surface areas. In this situation, prior to turning to other image segmentation or object recognition methods, we propose that a subdomain approach may be used first. Such an approach may be in the form of dividing the image domain into multiple subdomains or in the form of identifying user-selected ROIs. It is noted that the influence of noisy artifacts will appear in the first application example if a different ROI is selected (Figure6). Another practical issue, with a related solution, is that of image size. If images with very large sizes are encountered, computational cost may increase significantly; meanwhile robustness to noisy artifacts likely decreases. In this situation, the subdomain approach is an effective solution as well.


## 8.3. Computational Cost

The algorithms in the proposed framework, including the manifold-based monitoring algorithm and the level set-based crack detection and quantification algorithms, comprise a large number of basic vector or matrix operations with the complexity order atO(N) or O(N2). As a result, the associated computational cost is considerably larger than the traditional simple filtering-based methods. Computational cost is a concern if the framework is needed for online monitoring and quantification of rapidly captured image sequences. For offline processing of single images, however, the computation time is acceptable based on our implementation of the proposed algorithms in a Matlab environment [36]. If the framework is used as an online monitoring solution, fine tuning of the image capture rate and the image processing time is subject to further investigation. A plausible solution may be to subsample the captured images (i.e., only process a fraction of the images). It is noted that most of the time consumed from an operational point of view is spent tuning the internal parameters. This issue is expected to be mitigated when more empirical knowledge is obtained.


## 9. Conclusions

In this paper, a novel image-based framework is proposed to monitor and quantify one of the most common types of structural damage in civil engineering, namely, concrete surface cracks. In this case, optical camera images are used. Using illustrative examples, we first discuss the technological challenges in implementing any image-based damage identification approach using field or laboratory images. These include various instances of photometric and geometric distortions. For the latter, structural motion in the images forms one of the major challenges.By means of integrating two advanced image modeling algorithms, the proposed framework provides a full solution to the common damage identification objectives encountered in practice. First, through the novel use of manifold-distance computation based on multi-temporal images, motion-invariant features of crack occurrence are extracted. Second, a level set-based active contour algorithm is applied to realize a model-based boundary detection for concrete surface cracks. Based on the obtained level-set crack boundaries, an approximate method for determining continuous crack widths is developed. Two example applications are provided to demonstrate use of the proposed framework. Results show the success in monitoring temporally and quantifying geometrically the concrete cracks.



---

*Source: 215295-2010-06-07.xml*