Pothole Detection Based on Disparity Transformation and Road Surface Modeling

TL;DR

This paper introduces a robust and efficient computer vision-based method for pothole detection that combines disparity transformation, surface modeling, and outlier removal to achieve high accuracy in road maintenance tasks.

Contribution

The proposed approach uniquely integrates disparity transformation with 3D surface modeling and outlier reduction for improved pothole detection accuracy.

Findings

Detection accuracy of 98.7%

Pixel-level accuracy of 99.6%

Effective outlier handling with RANSAC

Abstract

Pothole detection is one of the most important tasks for road maintenance. Computer vision approaches are generally based on either 2D road image analysis or 3D road surface modeling. However, these two categories are always used independently. Furthermore, the pothole detection accuracy is still far from satisfactory. Therefore, in this paper, we present a robust pothole detection algorithm that is both accurate and computationally efficient. A dense disparity map is first transformed to better distinguish between damaged and undamaged road areas. To achieve greater disparity transformation efficiency, golden section search and dynamic programming are utilized to estimate the transformation parameters. Otsu's thresholding method is then used to extract potential undamaged road areas from the transformed disparity map. The disparities in the extracted areas are modeled by a quadratic surface using least squares fitting. To improve disparity map modeling robustness, the surface normal is also integrated into the surface modeling process. Furthermore, random sample consensus is utilized to reduce the effects caused by outliers. By comparing the difference between the actual and modeled disparity maps, the potholes can be detected accurately. Finally, the point clouds of the detected potholes are extracted from the reconstructed 3D road surface. The experimental results show that the successful detection accuracy of the proposed system is around 98.7% and the overall pixel-level accuracy is approximately 99.6%.

Tables

Table 1. TABLE I: Comparison of the Average σ d subscript 𝜎 𝑑 \sigma_{d} values between the two cases.

Dataset	Case 1	Case 2
KITTI stereo 2012 training dataset	0.4862	0.7800
KITTI stereo 2015 training dataset	0.5506	0.9428

Table 2. TABLE II: Comparison of Successful pothole Detection accuracy.

Dataset	Method	Total Potholes	Correct Detection	Incorrect Detection	Misdetection
Dataset 1	algorithm in [15]		11	11	0
	algorithm in [16]	22	22	0	0
	our algorithm		22	0	0
Dataset 2	algorithm in [15]		42	10	0
	algorithm in [16]	52	40	8	4
	our algorithm		51	1	0
Dataset 3	algorithm in [15]		5	0	0
	algorithm in [16]	5	5	0	0
	our algorithm		5	0	0
Total	algorithm in [15]		58	21	0
	algorithm in [16]	79	67	8	4
	our algorithm		78	1	0

Table 3. TABLE III: Comparison of Pixel-Level Precision, Recall, F-score and Accuracy.

Dataset	Method	recall	precision	F-score	accuracy
Dataset 1	[15]	0.5199	0.5427	0.5311	0.9892
	[16]	0.4622	0.9976	0.6317	0.9936
	proposed	0.4990	0.9871	0.6629	0.9940
Dataset 2	[15]	0.9754	0.9712	0.9733	0.9987
	[16]	0.8736	0.9907	0.9285	0.9968
	proposed	0.9804	0.9797	0.9800	0.9991
Dataset 3	[15]	0.6119	0.7714	0.6825	0.9948
	[16]	0.5339	0.9920	0.6942	0.9957
	proposed	0.5819	0.9829	0.7310	0.9961
Total	[15]	0.7799	0.8220	0.8004	0.9942
	[16]	0.6948	0.9921	0.8173	0.9954
	proposed	0.7709	0.9815	0.8635	0.9964

Equations

$$g(u,v)=c_0+c_{1}u+c_{2}v+c_3{u}^2+c_4{v}^2+c_{5}uv,$$

$$g(v)={\alpha }_0+{\alpha }_{1}v+{\alpha }_2{v}^2.$$

$$\begin{array}{r}\mathbf{\alpha }=\underset{\mathbf{\alpha }}{\mathrm{argmin}}{E}_0,\end{array}$$

$$E_0={\mathbf{e}}^{\top }\mathbf{e},$$

$$\mathbf{e}=\mathbf{d}-\mathbf{V}\mathbf{\alpha }.$$

$$\begin{array}{r}\mathbf{V}=\left[\begin{array}{ccc}1& v_0& {v_0}^2\\ 1& v_1& {v_1}^2\\ ⋮& ⋮& ⋮\\ 1& v_n& {v_n}^2\end{array}\right].\end{array}$$

$$\mathbf{\alpha }=({\mathbf{V}}^{\top }\mathbf{V}{)}^{-1}{\mathbf{V}}^{\top }\mathbf{d}.$$

$$E_{0_{\text{min}}}={\mathbf{d}}^{\top }\mathbf{d}-{\mathbf{d}}^{\top }\mathbf{V}({\mathbf{V}}^{\top }\mathbf{V}{)}^{-1}{\mathbf{V}}^{\top }\mathbf{d}.$$

$$\begin{array}{r}x=u\cos \theta +v\sin \theta ,\\ y=v\cos \theta -u\sin \theta .\end{array}$$

$$\mathbf{e}(\theta )=\mathbf{d}-\mathbf{Y}(\theta )\mathbf{\alpha }(\theta ),$$

$$\begin{array}{r}\mathbf{Y}(\theta )=\left[\begin{array}{ccc}1& y_0(\theta )& {y_0(\theta )}^2\\ 1& y_1(\theta )& {y_1(\theta )}^2\\ ⋮& ⋮& ⋮\\ 1& y_n(\theta )& {y_n(\theta )}^2\end{array}\right].\end{array}$$

$$\mathbf{\alpha }(\theta )=(\mathbf{Y}(\theta {)}^{\top }\mathbf{Y}(\theta ){)}^{-1}\mathbf{Y}(\theta {)}^{\top }\mathbf{d}.$$

$$E_{0_{\text{min}}}(\theta )={\mathbf{d}}^{\top }\mathbf{d}-{\mathbf{d}}^{\top }\mathbf{Y}(\theta )(\mathbf{Y}(\theta {)}^{\top }\mathbf{Y}(\theta ){)}^{-1}\mathbf{Y}(\theta {)}^{\top }\mathbf{d}.$$

$$\begin{array}{r}{E}_1(d,y)=-m(d,y)+\underset{\tau }{\min }[E_1(d+1,y+\tau )+\lambda \tau ]\\ \text{s.t.}\tau \in [{\tau }_{\text{min}},{\tau }_{\text{max}}],\end{array}$$

$$\tilde{d}(u,v,d,\theta )=[1y(u,v,\theta ){y(u,v,\theta )}^2]\mathbf{\alpha }(\theta )-d+\delta ,$$

$${\sigma }_{\text{o}}^2(T_{\text{o}})=P_0(T_{\text{o}})P_1(T_{\text{o}})[{\mu }_0(T_{\text{o}})-{\mu }_1(T_{\text{o}}){]}^2,$$

$$P_0(T_{\text{o}})=\sum _{i={\tilde{d}}_{\text{min}}}^{T_{\text{o}}-1}p(i),P_1(T_{\text{o}})=\sum _{i=T_{\text{o}}}^{{\tilde{d}}_{\text{max}}}p(i)$$

$$\begin{array}{r}{\mu }_0(T_{\text{o}})=\frac{1}{P_0(T_{\text{o}})}\sum _{i={\tilde{d}}_{\text{min}}}^{T_{\text{o}}-1}ip(i),\\ {\mu }_1(T_{\text{o}})=\frac{1}{P_1(T_{\text{o}})}\sum _{i=T_{\text{o}}}^{{\tilde{d}}_{\text{max}}}ip(i).\end{array}$$

$$X^{\mathcal{W}}=\frac{u{T}_c}{d},Y^{\mathcal{W}}=\frac{v{T}_c}{d},Z^{\mathcal{W}}=\frac{f{T}_c}{d},$$

$$\mathbf{c}=({\mathbf{W}}^{\top }\mathbf{W}{)}^{-1}{\mathbf{W}}^{\top }\mathbf{d},$$

$$\begin{array}{r}\mathbf{W}=\left[\begin{array}{cccccc}1& u_0& v_0& {u_0}^2& {v_0}^2& u_0{v}_0\\ 1& u_1& v_1& {u_1}^2& {v_1}^2& u_1{v}_1\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& u_n& v_n& {u_n}^2& {v_n}^2& u_n{v}_n\end{array}\right].\end{array}$$

$${\mathbf{Q}}_i^{+}=[{\mathbf{p}}_i,{{\mathbf{Q}}_i}^{\top }{]}^{\top }.$$

$$\boldsymbol{n}_{i}=\operatorname*{arg\,min}_{\boldsymbol{n}_{i}}\Big{|}\Big{|}\Big{[}\boldsymbol{Q}_{i}^{+}-\bar{\boldsymbol{Q}}_{i}^{+}\Big{]}\boldsymbol{n}_{i}\Big{|}\Big{|}_{2},$$

$$\bar{\boldsymbol{Q}}_{i}^{+}=\boldsymbol{1}_{k+1}{\big{(}\frac{1}{k+1}{\boldsymbol{Q}_{i}^{+}}^{\top}\boldsymbol{1}_{k+1}\big{)}}^{\top}.$$

$$\begin{array}{r}\hat{\mathbf{n}}=[\sin {\phi }_1\cos {\phi }_2,\sin {\phi }_1\sin {\phi }_2,\cos {\phi }_1{]}^{\top }\\ \text{s.t.}{\phi }_1\in [0,\pi ],{\phi }_2\in [0,2\pi ).\end{array}$$

$$E_2={1_{n+1}}^{\top }\mathbf{m},$$

$$\mathbf{m}=[-{\mathbf{n}}_0\cdot \hat{\mathbf{n}},-{\mathbf{n}}_1\cdot \hat{\mathbf{n}},\cdots ,-{\mathbf{n}}_n\cdot \hat{\mathbf{n}}{]}^{\top }.$$

$$\tan {\phi }_1=\frac{{1_{n+1}}^{\top }{\mathbf{n}}_u\cos {\phi }_2+{1_{n+1}}^{\top }{\mathbf{n}}_v\sin {\phi }_2}{{1_{n+1}}^{\top }{\mathbf{n}}_d},$$

$$\tan {\phi }_2=\frac{{1_{n+1}}^{\top }{\mathbf{n}}_v}{{1_{n+1}}^{\top }{\mathbf{n}}_u},$$

$${\phi }_2=\arctan \frac{{1_{n+1}}^{\top }{\mathbf{n}}_v}{{1_{n+1}}^{\top }{\mathbf{n}}_u}+k\pi ,k\in \{0,1\}.$$

Full text

I Introduction

ROAD potholes are considerably large structural failures on the road surface. They are caused by contraction and expansion of the road surface as rainwater permeates into the ground [1]. To ensure traffic safety, it is crucial and necessary to frequently inspect and repair road potholes [2]. Currently, potholes are regularly detected and reported by certified inspectors and structural engineers [3]. This task is, however, time-consuming and tedious [4]. Furthermore, the detection results are always subjective, because they depend entirely on personnel experience [5]. Therefore, automated pothole detection systems have been developed to recognize and localize potholes both efficiently and objectively.

Over the past decade, various technologies, such as active and passive sensing, have been utilized to acquire road data and aid personnel in detecting road potholes [5]. For example, Tsai and Chatterjee [6] mounted two laser scanners on a digital inspection vehicle (DIV) to collect 3D road surface data. These data were then processed using either semi or fully automatic methods for pothole detection. Such systems ensure personnel safety, but also reduce the need for manual intervention [6]. Furthermore, by comparing the road data collected over different periods, the traffic flow can be evaluated and the future road condition can be predicted [2]. The remainder of this section presents the state-of-the-art pothole detection algorithms and highlights the motivation, contributions and outline of this paper.

I-A State of the Art in Road Pothole Detection

I-A1 2D Image Analysis-Based Pothole Detection Algorithms

There are typically four main steps used in 2D image analysis-based pothole detection algorithms: a) image preprocessing; b) image segmentation; c) shape extraction; d) object recognition [5]. A color or gray-scale road image is first preprocessed, e.g., using morphological filters [7], to reduce image noise and enhance the pothole outline [5, 8]. The preprocessed road image is then segmented using histogram-based thresholding methods, such as Otsu’s [9] or the triangle [10] method. Otsu’s method minimizes the intra-class variance and performs better in terms of separating damaged and undamaged road areas [9]. The extracted region is then modeled by an ellipse [10]. Finally, the image texture within the ellipse is compared with the undamaged road area texture. If the former is coarser than the latter, the ellipse is considered to be a pothole [5].

However, both color and gray-scale image segmentation techniques are severely affected by various factors, most notably by poor illumination conditions [11]. Therefore, some authors proposed to perform segmentation on the depth maps, which has shown to achieve better performance when separating damaged and undamaged road areas [6, 11]. Furthermore, the shapes of actual potholes are always irregular, making the geometric and textural assumptions occasionally unreliable. Moreover, the pothole’s 3D spatial structure cannot always be explicitly illustrated in 2D road images [12]. Therefore, 3D road surface information is required to measure pothole volumes. In general, 3D road surface modeling-based pothole detection algorithms are more than capable of overcoming the aforementioned disadvantages.

I-A2 3D Road Surface Modeling-Based Pothole Detection Algorithms

The 3D road data used for pothole detection is generally provided by laser scanners [6], Microsoft Kinect sensors [11], or passive sensors [12, 13, 14, 15, 16, 4]. Laser scanners mounted on DIVs are typically used for accurate 3D road surface reconstruction. However, the purchase of laser scanning equipment and its long-term maintenance are still very expensive [3]. The Microsoft Kinect sensors were initially designed for indoor use. Therefore, they greatly suffer from infra-red saturation in direct sunlight [17]. For this reason, passive sensors, such as a single movable camera or multiple synchronized cameras, are more suitable for acquiring 3D road data and for pothole detection[12, 18]. For example, Zhang and Elaksher [13] mounted a single camera on an unmanned aerial vehicle (UAV) to reconstruct the road surface via Structure from Motion (SfM) [18]. A variety of stereo vision-based pothole detection methods have been developed as well [15, 16]. The 3D point cloud generated from a disparity map was interpolated into a quadratic surface using least squares fitting (LSF) [15]. The potholes were then recognized by comparing the difference between the 3D point cloud and the fitted quadratic surface. In [16], the surface modeling was performed on disparity maps instead of the point clouds, and random sample consensus (RANSAC) [19] was used to improve the pothole detection robustness.

I-B Motivation

Currently, laser scanning is still the main technology used to provide 3D road information for pothole detection, while other technologies, such as passive sensing, are under-utilized [2]. However, DIVs are not widely used, primarily because of the initial cost but also because their routine operation and long-term maintenance are still very costly [10]. Therefore, the trend is to equip DIVs with inexpensive, portable and durable sensors, such as digital cameras, for 3D road data acquisition. Stereo road image pairs can be used to calculate the disparity maps [14], which essentially represent the 3D road surface geometry. Recently, due to some major advances in computer stereo vision, road surface geometry can be reconstructed with a three-millimeter accuracy [12, 4]. Additionally, stereo cameras used for road data acquisition are inexpensive, portable and adaptable for different DIV types.

So far, comprehensive studies have been made in both 2D image analysis-based and 3D road surface modeling-based pothole detection. Unfortunately, these algorithms are always used independently [5]. Furthermore, pothole detection accuracy is still far from satisfactory [5]. Exploring effective approaches for disparity map preprocessing, by applying 2D image processing algorithms, is therefore also a popular area of research that requires more attention. Only the disparities in the potential undamaged road areas are then used for disparity map modeling.

Moreover, the surface normal vector is a very important descriptor, which is, however, rarely utilized in existing 3D road surface modeling-based pothole detection algorithms. In this paper, we improve disparity map modeling by eliminating the disparities whose surface normals differ significantly from the expected ones.

I-C Novel Contributions

In this paper, a robust stereo vision-based pothole detection system is introduced. The main contributions are: a) a novel disparity transformation algorithm; b) a robust disparity map modeling algorithm; c) three pothole detection datasets which have been made publicly available for research purposes. These datasets are also used in our experiments for assessing pothole detection accuracy.

Since the disparities in damaged road areas can severely affect the accuracy of disparity modeling, we first transform the disparity maps to better distinguish between damaged and undamaged road areas. To achieve greater processing efficiency, we use golden section search (GSS) [20] and dynamic programming (DP) [21] to estimate the transformation parameters. Otsu’s thresholding method [22] is then performed on the transformed disparity map to extract the undamaged road areas, where the disparities can be modeled by a quadratic surface using LSF. To improve the robustness of disparity map modeling, the surface normal information is also integrated into the modeling process. Furthermore, RANSAC is utilized to reduce the effects caused by any potential outliers. By comparing the difference between the actual and modeled disparity maps, the potholes can be detected effectively. Finally, different potholes are labeled using connected component labeling (CCL) [23] and their point clouds are extracted from the reconstructed 3D road surface.

I-D Paper Outline

The remainder of this paper is structured as follows: Section II details the proposed pothole detection algorithm. The experimental results for performance evaluation are illustrated in Section III. Finally, Section IV summarizes the paper and provides recommendations for future work.

II Pothole Detection Algorithm

The block diagram of the proposed pothole detection algorithm is illustrated in Fig. 1, where the algorithm consists of three main components: a) disparity transformation; b) undamaged road area extraction; c) disparity map modeling and pothole detection.

II-A Disparity Transformation

The input of this procedure is a dense disparity map having subpixel accuracy. Since the performance of disparity map modeling relies entirely on the disparity estimation accuracy, the dense disparity map was obtained from a stereo road image pair (see Fig. 2) through our disparity estimation algorithm [12], where the stereo matching search range propagates iteratively from the bottom of the image to the top, and the subpixel disparity map is refined by iteratively minimizing an energy function with respect to the interpolated correlation parabolas. The disparity map is shown in Fig. 2, and its corresponding v-disparity map is shown in Fig. 2. A v-disparity map can be created by computing the histogram of each horizontal row of the disparity map. The proposed pothole detection algorithm is based on the work presented in [16], where the disparities of the undamaged road surface are modeled by a quadratic surface as follows:

[TABLE]

where $u$ and $v$ are the horizontal and vertical disparity map coordinates, respectively. The origin of the coordinate system in (1) is at the center of the disparity map. Since in our experiments the stereo rig is mounted at a relatively low height, the curvature of the reconstructed road surface is not very high. This makes the values of $c_{1}$, $c_{3}$ and $c_{5}$ in (1) very close to zero, when the stereo rig is perfectly parallel to the horizontal road surface. In this case, the projection of the road disparities on the v-disparity map can be assumed to be a parabola of the form [21]:

[TABLE]

However, in practice, the stereo rig baseline is not always perfectly parallel to the horizontal road surface. This fact can introduce a non-zero roll angle $\theta$ (see Fig. 3) into the imaging process, where $T_{c}$ and $h$ represent the baseline and the height of the stereo rig, respectively. ${o}_{l}^{\mathscr{C}}$ and ${o}_{r}^{\mathscr{C}}$ are the origins of the left and right camera coordinate systems, respectively. ${O}^{\mathscr{W}}$ is the origin of the world coordinate system. An example of the resulting disparity map is shown in Fig. 2, where readers can clearly see that the disparity values change gradually in the horizontal direction, which makes the approach of representing the disparity projection using (2) somewhat problematic. In this regard, we first estimate the value of the roll angle. The effects caused by the non-zero roll angle are then eliminated by rotating the disparity map by $\theta$. Finally, the coefficients of the disparity projection model in (2) are estimated, and the disparity map is transformed to better distinguish between damaged and undamaged road areas.

II-A1 $\theta$ Estimation and Disparity Map Rotation

Over the past decade, considerable effort has been made to improve the roll angle estimation. The most commonly used device for this task is an inertial measurement unit (IMU). An IMU can measure the angular rate of a vehicle by analyzing the data acquired using different sensors, such as accelerometers, gyroscopes and magnetometers [24, 25]. In these approaches, the road bank angle is always assumed to be zero, and only the roll angle is considered in the estimation process. However, the estimation of both these two angles is always required in many real-world applications. Unfortunately, this cannot be realized with the use of only IMUs [26, 27].

In recent years, many authors have turned their focus towards estimating the roll angle from disparity maps [28, 29, 21, 30, 31]. For example, the road surface is assumed to be a horizontal ground plane, and an effective roll angle estimation algorithm was proposed based on v-disparity map analysis [28, 29]. In [21], the disparities in a selected small area were modeled by a plane $g(u,v)=c_{0}+c_{1}u+c_{2}v$. The roll angle was then calculated as $\arctan(-c_{1}/c_{2})$. However, finding a proper disparity map area for plane fitting is always challenging, because the selected area may contain an obstacle or a pothole, which can severely affect the fitting accuracy [30]. Furthermore, the above-mentioned algorithms are only suitable for planar road surfaces. Hence, in this subsection, we introduce a roll angle estimation algorithm, which can work effectively for both planar and non-planar road surfaces.

When the roll angle is equal to zero, the vector $\boldsymbol{\alpha}=[\alpha_{0},\alpha_{1},\alpha_{2}]^{\top}$, storing the disparity projection model coefficients, can be estimated by solving a least squares problem as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

The column vector $\boldsymbol{d}=[d_{0},\ d_{1},\ \cdots,\ d_{n}]^{\top}$ stores the disparity values. $\boldsymbol{V}$ is a matrix of size $(n+1)\times 3$ given as follows:

[TABLE]

This optimization problem has a closed form solution:

[TABLE]

The minimum energy $E_{0_{\text{min}}}$ can also be obtained by combining (4), (5) and (7):

[TABLE]

However, when the roll angle does not equal zero, the disparity distribution on each row becomes less compact (see Fig. 2). This greatly affects the accuracy of least squares fitting and produces a much higher $E_{0_{\text{min}}}$.

To rotate the disparity map around a given angle $\theta$, each set of original coordinates $[u,v]^{\top}$ is transformed to a set of new coordinates $[x,y]^{\top}$ as follows [31]:

[TABLE]

(5) can now be rewritten as follows:

[TABLE]

where $\boldsymbol{Y}(\theta)$ is an $(n+1)\times 3$ matrix:

[TABLE]

(7) is rewritten as follows:

[TABLE]

(4), (10) and (12) result in the following expression:

[TABLE]

Therefore, the main consideration of the proposed roll angle estimation algorithm is to rotate the disparity map at different angles, and find the angle which minimizes $E_{0_{\text{min}}}$. $E_{0_{\text{min}}}$ with respect to different $\theta$ is illustrated in Fig. 4. Giving a set of coordinates $[u,v]^{\top}$, the new coordinate $y$ can be calculated using (9). The corresponding $E_{0_{\text{min}}}$ can be computed from (13). Due to the fact that $\cos(\theta+\pi)=-\cos\theta$ and $\sin(\theta+\pi)=-\sin\theta$, the disparity maps rotated around $\theta$ and $\theta+\pi$ are symmetric with respect to the origin of the coordinate system. Namely, (13) outputs the same $E_{0_{\text{min}}}$ no matter how the disparity map is rotated around $\theta$ or $\theta+\pi$. Therefore, we set the interval of $\theta$ to $(-\pi/2,\pi/2]$. The estimation of $\theta$ is achieved by finding the position of the local minima between $-\pi/2$ and $\pi/2$.

However, finding the local minima is a computationally intensive task, because it involves performing the necessary calculations through the whole interval of $\theta$. Furthermore, the step size $\varepsilon_{\theta}$ has to be set to a very small and practical value, in order to obtain an accurate value of $\theta$ [31]. Hence in this paper, GSS is utilized to reduce the searching times. The procedure of the proposed $\theta$ estimation algorithm is given in Algorithm 1, where $\kappa=(\sqrt{5}-1)/2$ is the golden section ratio [20].

The disparity map is then rotated around $\theta$, as illustrated in Fig. 5. A y-disparity map (see Fig. 5) can be created by computing the histogram of each horizontal row of the rotated disparity map. We can observe that the disparity values on each row become more uniform. The evaluation of the proposed roll angle estimation algorithm will be discussed in Section III-B.

II-A2 $\boldsymbol{\alpha}$ Estimation and Disparity Transformation

In this subsection, we utilize DP [21] to extract the road disparity projection model from the y-disparity map. For the purpose of convenient notation, the projection model is also referred to as the target path. The energy of every possible solution is first computed as follows:

[TABLE]

where $m(d,y)$ represents the y-disparity value at the position $[d,y]^{\top}$ and $\lambda$ is a positive smoothness term [21]. $E_{1}$ represents the energy of a possible target path in the y-disparity map. $\tau_{\text{min}}$ is typically set to 0. $\tau_{\text{max}}$ depends entirely on $\boldsymbol{{\alpha}}$, and it is set to 10 in this paper. The target path $M=\{(d_{i},y_{i}),\ i=0,1,\dots,n\}$ can be found by minimizing the energy function in (14), where $(d_{i},y_{i})$ stores the horizontal and vertical coordinates of the target path, respectively.

By substituting the horizontal and vertical coordinates of the target path into (3), (4), (10) and (11), we can obtain $\boldsymbol{{\alpha}}$. The disparity map can, therefore, be transformed using $\theta$ and $\boldsymbol{\alpha}$. However, the outliers in the target path may greatly affect the accuracy of $\boldsymbol{{\alpha}}$ estimation. We, therefore, use RANSAC to update the values in $\boldsymbol{{\alpha}}$. The full list of procedures involved for $\boldsymbol{\alpha}$ estimation are detailed in Algorithm 2.

RANSAC is iterated $t$ times. Selecting a higher $t$ raises the possibility of finding the best $\boldsymbol{\alpha}$ but also increases the processing time. In order to minimize the trade-off between speed and robustness, $t$ is set to $50$ in this paper. In each iteration, we select $p$ pairs of coordinates $[d,y]^{\top}$ from the target path to estimate $\boldsymbol{{\alpha}}$. For a smaller $p$, there is less chance that any outliers will influence optimization. In this paper $p$ is set to $3$, which is the smallest possible value for determining $\boldsymbol{{\alpha}}$. The ratio $\eta$ of inliers versus outliers can then be computed with respect to a given tolerance $\varepsilon_{\boldsymbol{{\alpha}}}$. The best $\boldsymbol{{\alpha}}$ corresponds to the highest ratio $\eta$. However, the selection of an appropriate $\varepsilon_{\boldsymbol{{\alpha}}}$ can be challenging, as it is possible that there could be more than one satisfying value for $\boldsymbol{{\alpha}}$. Hence, in this paper, the value of $\varepsilon_{\boldsymbol{{\alpha}}}$ also changes $s$ times. In each iteration, the value of $\varepsilon_{\boldsymbol{{\alpha}}}$ reduces by half and $\eta$ is computed. The best $\boldsymbol{{\alpha}}$ can be determined by finding the highest ratio $\eta$. In this paper, the values of $\varepsilon_{\boldsymbol{{\alpha}}}$ and $s$ are both set to $4$.

Finally, each disparity is transformed using the following equation:

[TABLE]

where $\delta$ is a constant value set to guarantee that all the transformed disparity values are non-negative. In this paper, we set $\delta$ to $30$. The transformed disparity map is shown in Fig. 6. It can be clearly seen that the disparity values in the undamaged road areas become more uniform, while they differ significantly from those in the damaged areas (potholes). This makes the extraction of undamaged road regions much simpler.

II-B Undamaged Road Area Extraction

Next, we utilize Otsu’s thresholding method to segment the transformed disparity map. The segmentation threshold $T_{\text{o}}$ can be obtained by maximizing the inter-class variance ${\sigma}_{\text{o}}^{2}$ as follows [22]:

[TABLE]

where

[TABLE]

represent the probabilities of damaged and undamaged road areas, respectively. $p(i)$ is the probability of $\tilde{d}=i$. The average disparity values of the damaged and undamaged road areas are given by:

[TABLE]

The segmentation result is shown in Fig. 6. We can see that the undamaged road area is successfully extracted. However, Otsu’s thresholding method will always classify the disparities into two categories, even if the transformed disparity map does not contain any road damage. We, therefore, carry out disparity map modeling to ensure that the potholes are correctly detected.

II-C Disparity Map modeling and Pothole Detection

A common practice in 3D modeling-based pothole detection algorithms [15, 16] is to fit a quadratic surface to either a 3D point cloud or a 2D disparity map. In parallel axis stereo vision, the point cloud is generated from the disparity map as follows [21]:

[TABLE]

where $f$ is the camera focal length. A disparity error larger than one pixel may result in a non-negligible difference in the point cloud [32]. Therefore, disparity map modeling can avoid such errors generated from (19), producing greater accuracy compared to point cloud modeling.

To model the disparity map, $\boldsymbol{c}=[c_{0},c_{1},c_{2},c_{3},c_{4},c_{5}]^{\top}$, storing the quadratic surface model coefficients can be estimated as follows:

[TABLE]

where

[TABLE]

However, potential outliers can severely affect the accuracy of disparity map modeling and therefore need to be discarded beforehand. In this subsection, a disparity map point is determined as an outlier if it fulfills one of the following conditions:

• •

it is located in one of the damaged road areas.

• •

its surface normal vector differs greatly from the optimal one.

• •

its disparity value differs greatly from the one computed using Eq. 1.

In Section II-B, the undamaged road areas are successfully extracted and we only use the disparities in this area to model the disparity map. The rest of this subsection presents the approaches for determining the outliers which satisfy the last two conditions and the process of modeling the disparity map without using these outliers.

II-C1 Optimal Normal Vector Estimation

For each point $\boldsymbol{p}_{i}=[u_{i},\ v_{i},\ d_{i}]^{\top}$ in the undamaged road area, we would like to estimate a normal vector $\boldsymbol{n}_{i}=[n_{ui},\ n_{vi},\ n_{di}]^{\top}$ from a set of $k$ points in its neighborhood $\boldsymbol{Q}_{i}=[\boldsymbol{q}_{i1},\ \boldsymbol{q}_{i2},\ \cdots,\ \boldsymbol{q}_{ik}]^{\top}$, where $\boldsymbol{q}_{ij}\neq\boldsymbol{p}_{i}$. Here, we define the augmented neighbor matrix $\boldsymbol{Q}_{i}^{+}$ which contains all neighbors and the point $\boldsymbol{p}_{i}$ itself as follows:

[TABLE]

Existing normal vector estimation methods are generally classified into one of two categories: optimization-based and averaging-based [33]. Although the performance of normal vector estimation depends primarily on the application itself, PlanePCA [34], an optimization-based normal vector estimation method, has superior performance in terms of both speed and accuracy. Hence in this subsection, we utilize PlanePCA to estimate the normal vectors of the disparities. $\boldsymbol{n}_{i}$ can be estimated as follows:

[TABLE]

where

[TABLE]

$\boldsymbol{1}_{m}$ represents an $m\times 1$ vector of ones. Due to the fact that the normal vectors are normalized, they can be projected on a sphere, as shown in Fig. 7, where we can clearly see that the projections are distributed in a small area. Therefore, the optimal normal vector $\hat{\boldsymbol{n}}$ can be determined by finding the position at which the projections distribute most intensively.

Since the projection of $\hat{\boldsymbol{n}}$ is also on the sphere, it can be written in spherical coordinates as follows:

[TABLE]

It can be estimated by minimizing $E_{2}$:

[TABLE]

where

[TABLE]

By applying (25) and (27) to (26), the following expressions are derived:

[TABLE]

[TABLE]

where $\boldsymbol{n}_{u}=[n_{u0},n_{u1},\cdots,n_{un}]^{\top}$, $\boldsymbol{n}_{v}=[n_{v0},n_{v1},\cdots,n_{vn}]^{\top}$ and $\boldsymbol{n}_{d}=[n_{d0},n_{d1},\cdots,n_{dn}]^{\top}$. Due to the fact that $\varphi_{2}$ is between [math] and $2\pi$, (29) will have two solutions:

[TABLE]

Substituting each $\varphi_{2}$ into (28) produces a value for $\varphi_{1}$. The two pairs of $[\varphi_{1},\ \varphi_{2}]^{\top}$ correspond to the maxima and minima of $E_{2}$, respectively. By substituting each pair of $[\varphi_{1},\ \varphi_{2}]^{\top}$ into Eq. 26 and comparing the two obtained values, we can find the optimal normal vector. If the angle between $\hat{\boldsymbol{n}}$ and $\boldsymbol{n}_{i}$ exceeds a pre-set threshold $\varepsilon_{\boldsymbol{n}}$, the corresponding disparity will be considered as an outlier and will not be used for disparity map modeling. In our experiments, we assume that the second category of outliers account for $10\%$ of the undamaged road areas, and therefore, $\varepsilon_{\boldsymbol{n}}$ is set to $\pi/36$ rad. The outliers satisfying the first two conditions can then be successfully removed. The third category of outliers are removed along with the disparity map modeling.

II-C2 Disparity Map Modeling

To model the disparity map with more robustness, we use RANSAC to reduce the effects caused by the third category of outliers described in Section II-C. Here, RANSAC is iterated $t$ times. In each iteration, a subset of disparities are selected randomly to estimate $\boldsymbol{c}$. To ensure uniform distribution of the selected disparities, we equally divide the disparity map into a group of square blocks and select one disparity from each block. The disparity block size is $r\times r$. As $r$ becomes smaller, more disparities will be used for surface fitting, which increases the computational complexity. In contrast, selecting a higher value for $r$ results in less computational complexity, but potentially increases noise sensitivity. In this paper, the value of $r$ is set to $125$, which produces approximately $100$ square blocks for our disparity maps. In each iteration, the differences between the actual and fitted disparities are computed and the ratio $\eta$ between the inliers and outliers are obtained. $\boldsymbol{c}$ which corresponds to the highest $\eta$ is then selected as the desirable surface coefficients. Algorithm 2 presents more details on the least squares fitting using RANSAC. The modeled disparity map is shown in Fig. 8.

II-C3 Pothole Detection

The potholes can then be detected by finding the regions where the differences between the actual and modeled disparities are larger than a pre-set threshold $\varepsilon_{d}$. Before labeling different potholes using CCL, the connected components containing fewer than $w$ pixels are removed, because they severely affect the pothole labeling accuracy. Furthermore, the small holes in each connected component are filled, as they are considered to be noise. The selection of $\varepsilon_{d}$ and $w$ is discussed in Section III-D. The detected potholes are shown in Fig. 8. Finally, the point clouds of the detected potholes are extracted from the reconstructed 3D road surface. The corresponding results are shown in Fig. 9.

III Experimental Results

In this section, the performance of the proposed pothole detection algorithm is evaluated both qualitatively and quantitatively. The proposed algorithm was implemented in MATLAB on an Intel Core i7-8700K CPU using a single thread. The following subsections provide details on the experimental set-up and the evaluation of the proposed algorithm.

III-A Experimental Set-Up

In this work, we utilized a ZED stereo camera111https://www.stereolabs.com/ to capture stereo road images. An example of the experimental set-up is shown in Fig. 10. The stereo camera is calibrated manually using the stereo calibration toolbox in MATLAB R2018b. Using the above-mentioned experimental set-up, we created three datasets containing 67 pairs of stereo images. The image resolutions of dataset 1, 2 and 3 are $1028\times 1730$, $1030\times 1720$, $1028\times 1710$ pixels, respectively. The disparity maps are estimated using our previously published algorithm [12]. All datasets are publicly available and can be found at: ruirangerfan.com.

The following subsections analyze the accuracy of roll angle estimation, disparity transformation and pothole detection.

III-B Evaluation of Roll Angle Estimation

In this subsection, we first analyze the computational complexity of the proposed roll angle estimation algorithm. When estimating the roll angle without using GSS, we have to search through the whole interval of $(-\frac{\pi}{2},\frac{\pi}{2}]$ to find the local minima. Therefore, the computational complexity is $\mathcal{O}(\frac{\pi}{\varepsilon_{\theta}})$. In our method, GSS reduces the interval size exponentially. As a result, the interval size then becomes $\kappa^{n}\pi$ after the $n$-th iteration. Therefore, the proposed roll angle estimation algorithm reduces the computational complexity to $\mathcal{O}(\log_{\kappa}\frac{\varepsilon_{\theta}}{\pi})$. The proposed roll angle estimation algorithm needs 21 iterations to produce a roll angle, with an accuracy higher than $\frac{\pi}{18000}$ rad.

To evaluate the accuracy of the proposed roll angle estimation algorithm, we utilize a synthesized stereo dataset from EISATS [35, 36] where the roll angle is perfectly zero. Some experimental results are shown in Fig. 11. The road areas (see the magenta regions in the first row of Fig. 11) are manually selected and the disparities in these areas are utilized to estimate the roll angle $\hat{\theta}$. The absolute difference between the actual and estimated roll angles, i.e., $\Delta\theta=|\theta-\hat{\theta}|$, is computed for each frame. The average $\Delta\theta$ is approximately $1.129\times 10^{-4}$ rad which is lower than $\frac{\pi}{18000}$ rad. Therefore, the proposed algorithm is capable of estimating the roll angle with high accuracy.

III-C Evaluation of Disparity Transformation

Since the datasets we created only contain the ground truth of potholes, KITTI stereo datasets [37, 38] are utilized to quantify the performance of our proposed disparity transformation algorithm (the numbers of disparity maps in the KITTI stereo 2012 and 2015 datasets are 194 and 200, respectively). Some experimental results are shown in Fig. 12. Due to the fact that the proposed algorithm focuses entirely on the road surface, we manually selected a region of interest (see the purple areas in the first row) in each image to evaluate the performance of our algorithm. The corresponding transformed disparity maps are shown in the third row of Fig. 12, where readers can clearly see that the disparities in the road areas tend to have similar values. To quantify the accuracy of the transformed disparity maps, we compute the standard deviation $\sigma_{d}$ of the transformed disparity values as follows:

[TABLE]

where $\boldsymbol{\tilde{d}}=[\tilde{d}_{0},\ \tilde{d}_{1},\ \cdots,\ \tilde{d}_{m}]^{\top}$ is a column vector storing the transformed disparity values. The average $\sigma_{d}$ values of the two KITTI stereo datasets are provided in Table I.

In case 1, the effects caused by the non-zero roll angle are eliminated, while in case 2, the roll angle is assumed to be zero. From Table I, we can observe that the average $\sigma_{d}$ values of these two datasets reduce by approximately half when the effects caused by the non-zero roll angle are eliminated. The average $\sigma_{d}$ value of these two datasets is only 0.5188 pixels. Therefore, the proposed disparity transformation algorithm performs accurately and the transformed disparity values become very uniform. The runtime of disparity transformation is about 142 ms. In the next subsection, we will analyze the accuracy of pothole detection.

III-D Evaluation of Pothole Detection

In Section II-C, a set of randomly selected disparities are modeled as a quadratic surface. The potholes are detected by comparing the difference between the actual disparity map and the modeled quadratic surface. If a connected component contains more than $w$ pixels and the disparity difference of each pixel exceeds $\varepsilon_{d}$, it will be identified as a pothole. In our experiments, we utilize the brute-force search method to find the best values of $\varepsilon_{d}$ and $w$. The search range for $\varepsilon_{d}$ and $w$ are set to $[3.0,8.5]$ and $[100,5000]$, respectively. The step sizes for searching $\varepsilon_{d}$ and $w$ are set to $0.1$ and $100$, respectively.

For the first step, we go through the whole search range and record the number of detected potholes $\hat{n}_{\text{PD}}$ in each frame. The absolute difference $\Delta n_{\text{PD}}$ between each $\hat{n}_{\text{PD}}$ and the expected pothole number $n_{\text{PD}}$ is then computed. The sum of $\Delta n_{\text{PD}}$ with respect to a pair of given $\varepsilon_{d}$ and $w$ can therefore be obtained, as illustrated in Fig. 13. In our experiments, the least sum of $\Delta n_{\text{PD}}$ is equal to one and it is achieved only when $\varepsilon_{d}=6.2$ and $\varepsilon_{d}=3100$. The corresponding incorrect pothole detection result is shown in Fig. 14. Incorrect detection occurs when the middle of the first pothole subsides and the selected parameters cause the system to detect two potholes instead of one. Some examples of successful detection results are shown in the fifth row of Fig. 15, and the corresponding ground truth is shown in the sixth row.

We also compare our proposed algorithm with those produced in [15] and [16]. The pothole detection results obtained using the algorithms presented in [15] and [16] are shown in the third and forth rows of Fig. 15, respectively. The comparative pothole detection results are provided in Table II, where we can see that the successful detection accuracy achieved using [15] and [16] are $73.4\%$ and $84.8\%$, respectively. Compared to them, our proposed algorithm can detect potholes more accurately with a successful detection accuracy of $98.7\%$.

We also compare the proposed algorithm with [15] and [16] with respect to the pixel-level precision, recall, F-score and accuracy:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $n_{\text{TP}}$, $n_{\text{FP}}$, $n_{\text{FN}}$ and $n_{\text{TN}}$ are true positive, false positive, false negative and true negative pixel numbers, respectively.

The comparisons with respect to these four indicators are illustrated in Table III. It can be seen that our proposed algorithm outperforms [15] and [16], in terms of both pixel-level accuracy and F-score. It achieves an intermediate performance in terms of precision and recall. In addition, the precision and recall achieved using our proposed algorithm are very close to the highest values between [15] and [16]. Therefore, the proposed pothole detection algorithm performs both robustly and accurately.

IV Conclusion and Future Work

The main contributions of this paper are a novel disparity transformation algorithm and a disparity map modeling algorithm. Using our method, undamaged road areas are better distinguishable in the transformed disparity map and can be easily extracted using Otsu’s thresholding method. This greatly improves the robustness of disparity map modeling. To achieve greater processing efficiency, GSS and DP were utilized to estimate the transformation parameters. Furthermore, the disparities, whose normal vectors differ greatly from the optimal one, were also discarded in the process of disparity map modeling, which further improves the accuracy of the modeled disparity map. Finally, the potholes were detected by comparing the difference between the actual and modeled disparity maps. The point clouds of the detected potholes were then extracted from the reconstructed 3D road surface. In addition, we also created three datasets to contribute to stereo vision-based pothole detection research. The experimental results show that the overall successful detection accuracy of our proposed algorithm is around $98.7\%$ and the pixel-level accuracy is approximately $99.6\%$.

However, the parameters set for pothole detection cannot be applied to all cases. Therefore, we plan to train a deep neural network to detect potholes from the transformed disparity map. Furthermore, road surfaces cannot always be considered to be quadratic. Thus, we aim to design an algorithm to segment the reconstructed road surfaces into a group of localized planes prior to applying the proposed pothole detection algorithm.