## Pothole Detection Based on Disparity Transformation and Road Surface Modeling

Rui Fan, Member, IEEE , Umar Ozgunalp, Member, IEEE , Brett Hosking, Member, IEEE , Ming Liu, Senior Member, IEEE , Ioannis Pitas, Fellow, IEEE

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% GLYPH&lt;147&gt; and the overall pixel-level accuracy is approximately 99 6% GLYPH&lt;147&gt; .

Index Terms -pothole detection, computer vision, road surface modeling, disparity map, golden section search, dynamic programming, surface normal.

## I. INTRODUCTION

R OAD 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,

- R. Fan is with the Robotics Institute, Department of Electronic and Computer Engineering, the Hong Kong University of Science and Technology, Hong Kong SAR, China (e-mail: rui.fan@ieee.org).
- U. Ozgunalp is with the Electrical and Electronic Engineering, Near East University, Nicosia, Cyprus (e-mail: umar.ozgunalp@neu.edu.tr).
- B. Hosking is with the National Oceanography Centre, Southampton, U.K. (e-mail: wilski@noc.ac.uk).
- M. Liu is with the Robotics Institute, Department of Electronic and Computer Engineering, the Hong Kong University of Science and Technology, Hong Kong SAR, China (e-mail: eelium@ust.hk).
- I. Pitas is with the School of Informatics, Aristotle University of Thessaloniki, Thessaloniki, Greece (e-mail: pitas@aiia.csd.auth.gr).

This work was supported by the National Natural Science Foundation of China, under grant No. U1713211, the Research Grant Council of Hong Kong SAR Government, China, under Projects No. 11210017 and No. 21202816, and the Shenzhen Science, Technology and Innovation Commission (SZSTI), under grant JCYJ20160428154842603, awarded to Prof. Ming Liu.

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.

## A. State of the Art in Road Pothole Detection

1) 2D Image Analysis-Based Pothole Detection Algorithms: There are typically four main steps used in 2D image analysisbased 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.

2) 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 [4], [12]-[16]. 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.

## 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 [4], [12]. 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 modelingbased 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.

## 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.

## 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.

## 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(a)) 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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-4d030126-e864-4f86-b17b-c6f6c16c08b2.jpg)

Fig. 1. The block diagram of the proposed pothole detection algorithm.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-9145383b-44fc-43d4-9277-38b4e1645c37.jpg)

(a)

Fig. 2. Disparity map when the roll angle does not equal zero: (a) stereo road images; (b) disparity map; (c) v-disparity map.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-b574e0e5-c14e-4c4b-b925-74f29bab41c9.jpg)

parabolas. The disparity map is shown in Fig. 2(b), and its corresponding v-disparity map is shown in Fig. 2(c). A vdisparity 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:

$$g ( u, v ) = c _ { 0 } + c _ { 1 } u + c _ { 2 } v + c _ { 3 } u ^ { 2 } + c _ { 4 } v ^ { 2 } + c _ { 5 } u v, \quad ( 1 ) \quad \text{this re} _ { c e }$$

where 𝑢 and 𝑣 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 𝑐 1 , 𝑐 3 and 𝑐 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

Fig. 3. Roll angle definition in a stereo vision system.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-e6f38092-d2eb-4f1d-909b-d6a5ead54abe.jpg)

on the v-disparity map can be assumed to be a parabola of the form [21]:

$$g ( v ) = \alpha _ { 0 } + \alpha _ { 1 } v + \alpha _ { 2 } v ^ { 2 }. \text{ \quad \ \ } ( 2 )$$

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 𝜃 (see Fig. 3) into the imaging process, where 𝑇 𝑐 and ℎ represent the baseline and the height of the stereo rig, respectively. 𝑜 𝒞 𝑙 and 𝑜 𝒞 𝑟 are the origins of the left and right camera coordinate systems, respectively. 𝑂 𝒲 is the origin of the world coordinate system. An example of the resulting disparity map is shown in Fig. 2(b), 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 𝜃 . 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.

1) 𝜃 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 [21], [28]-[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 𝑔 𝑢GLYPH&lt;150&gt; ' 𝑣 ' = 𝑐 0 , 𝑐 1 𝑢 , 𝑐 2 𝑣 . The roll angle was then calculated as arctan 'GLYPH&lt;0&gt; 𝑐 1 GLYPH&lt;157&gt; 𝑐 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 𝜶 = » 𝛼 GLYPH&lt;150&gt; 0 𝛼 GLYPH&lt;150&gt; 1 𝛼 2 … &gt; , storing the disparity projection model coefficients, can be estimated by solving a least squares problem as follows:

$$\alpha = \underset { \alpha } { \arg \min } \ E _ { 0 }, \text{ \quad \ \ } \text{ \quad \ } \text{(4),}$$

where and

$$e = d - V \alpha. \quad \text{ \quad \ \ } ( 5 ) \quad \text{ang}$$

The column vector 𝒅 = » 𝑑 GLYPH&lt;150&gt; 0 𝑑 GLYPH&lt;150&gt; 1 GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;150&gt; 𝑑 𝑛 … &gt; stores the disparity values. 𝑽 is a matrix of size ' 𝑛 , 1 ' GLYPH&lt;2&gt; 3 given as follows:

$$V = \begin{bmatrix} 1 & v _ { 0 } & v _ { 0 } ^ { 2 } \\ 1 & v _ { 1 } & v _ { 1 } ^ { 2 } \\ \vdots & \vdots & \vdots \\ 1 & v _ { n } & v _ { n } ^ { 2 } \end{bmatrix}, \quad & \text{the} \, \begin{matrix} & \text{us} \\ & \text{the} \\ & & ( 6 ) & \\ & \text{with} \\ & & \text{is ro} \end{matrix}$$

4

5

This optimization problem has a closed form solution:

$$\alpha = ( V ^ { \top } V ) ^ { - 1 } V ^ { \top } d. \quad \ \ \ ( 7 ) ^ { \quad \text{the} \, 1 }$$

The minimum energy 𝐸 0 min can also be obtained by combining (4), (5) and (7):

$$E _ { 0 _ { \min } } = d ^ { \top } d - d ^ { \top } V ( V ^ { \top } V ) ^ { - 1 } V ^ { \top } d. \quad \quad ( 8 ) \quad \text{step s}$$

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

To rotate the disparity map around a given angle 𝜃 , each set of original coordinates » 𝑢GLYPH&lt;150&gt; 𝑣 … &gt; is transformed to a set of new coordinates » 𝑥GLYPH&lt;150&gt; 𝑦 … &gt; as follows [31]:

$$x & = u \cos \theta + v \sin \theta, & \quad \text{(9)} \quad \text{row} \\ y & = v \cos \theta - u \sin \theta. & \quad \text{(9)} \quad \text{(and)}$$

$$E _ { 0 } = e ^ { \top } e, \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \quad \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Fig. 4. 𝐸 0 min function versus 𝜃 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-f9fbe1b8-61a7-4fb0-bbc2-9c246376e526.jpg)

(5) can now be rewritten as follows:

$$e ( \theta ) = d - Y ( \theta ) \alpha ( \theta ), \text{ \quad \ \ } ( 1 0 )$$

where 𝒀 ' 𝜃 ' is an ' 𝑛 , 1 ' GLYPH&lt;2&gt; 3 matrix:

$$Y ( \theta ) = \begin{bmatrix} 1 & y _ { 0 } ( \theta ) & y _ { 0 } ( \theta ) ^ { 2 } \\ 1 & y _ { 1 } ( \theta ) & y _ { 1 } ( \theta ) ^ { 2 } \\ \vdots & \vdots & \vdots \\ 1 & y _ { n } ( \theta ) & y _ { n } ( \theta ) ^ { 2 } \end{bmatrix}. \quad \quad ( 1 1 )$$

(7) is rewritten as follows:

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

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

$$E _ { 0 _ { \min } } ( \theta ) = d ^ { \top } d - d ^ { \top } Y ( \theta ) ( Y ( \theta ) ^ { \top } Y ( \theta ) ) ^ { - 1 } Y ( \theta ) ^ { \top } d. \quad ( 1 3 )$$

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 𝐸 0 min . 𝐸 0 min with respect to different 𝜃 is illustrated in Fig. 4. Giving a set of coordinates » 𝑢GLYPH&lt;150&gt; 𝑣 … &gt; , the new coordinate 𝑦 can be calculated using (9). The corresponding 𝐸 0 min can be computed from (13). Due to the fact that cos ' 𝜃 , 𝜋 ' = GLYPH&lt;0&gt; cos 𝜃 and sin ' 𝜃 , 𝜋 ' = GLYPH&lt;0&gt; sin 𝜃 , the disparity maps rotated around 𝜃 and 𝜃 , 𝜋 are symmetric with respect to the origin of the coordinate system. Namely, (13) outputs the same 𝐸 0 min no matter how the disparity map is rotated around 𝜃 or 𝜃 , 𝜋 . Therefore, we set the interval of 𝜃 to 'GLYPH&lt;0&gt; 𝜋 GLYPH&lt;157&gt; 2 GLYPH&lt;150&gt; 𝜋 GLYPH&lt;157&gt; 2 … . The estimation of 𝜃 is achieved by finding the position of the local minima between GLYPH&lt;0&gt; 𝜋 GLYPH&lt;157&gt; 2 and 𝜋 GLYPH&lt;157&gt; 2 .

However, finding the local minima is a computationally intensive task, because it involves performing the necessary calculations through the whole interval of 𝜃 . Furthermore, the step size 𝜀 𝜃 has to be set to a very small and practical value, in order to obtain an accurate value of 𝜃 [31]. Hence in this paper, GSS is utilized to reduce the searching times. The procedure of the proposed 𝜃 estimation algorithm is given in Algorithm 1, where 𝜅 = ' p 5 GLYPH&lt;0&gt; 1 'GLYPH&lt;157&gt; 2 is the golden section ratio [20].

The disparity map is then rotated around 𝜃 , as illustrated in Fig. 5(a). A y-disparity map (see Fig. 5(b)) 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.

Fig. 5. Elimination of the effects caused by the non-zero roll angle: (a) rotated disparity map; (b) the y-disparity map of Fig. 5(a).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-3aeaa8f8-4fd3-4eb5-aa8d-06b986b0d93c.jpg)

## Algorithm 1: 𝜃 estimation using GSS.

Data:

disparity map

Result: 𝜃

1 set 𝜃 1 and 𝜃 2 to GLYPH&lt;0&gt; 𝜋 GLYPH&lt;157&gt; 2 and 𝜋 GLYPH&lt;157&gt; 2 , respectively;

2 compute 𝐸 0 min ' 𝜃 1 ' using Eq. 13;

3 compute 𝐸 0 min ' 𝜃 2 ' using Eq. 13;

- 4 while 𝜃 2 GLYPH&lt;0&gt; 𝜃 1 ¡ 𝜀 𝜃 do
- 5 set 𝜃 3 and 𝜃 4 to 𝜅𝜃 1 , ' 1 GLYPH&lt;0&gt; 𝜅 ' 𝜃 2 and 𝜅𝜃 2 , ' 1 GLYPH&lt;0&gt; 𝜅 ' 𝜃 1 , respectively;
- 6 compute 𝐸 0 min ' 𝜃 3 ' using Eq. 13;
- 7 compute 𝐸 0 min ' 𝜃 4 ' using Eq. 13;

8

if

𝐸

0

min

'

𝜃

3

'

¡ 𝐸

0

min

'

𝜃

4

'

then

9

10

11

12

𝜃

else

𝜃

end

13

end

- 2) 𝜶 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:

$$E _ { 1 } ( d, y ) = - m ( d, y ) + \min _ { \tau } [ E _ { 1 } ( d + 1, y + \tau ) + \lambda \tau ] \quad \begin{array} { c } \text{more un} \\ \text{value o} \\ \text{determi} \end{array} \quad ( 1 4 ) \quad \text{(14)} \quad \text{of $\varepsilon_{\alpha} \, \mathfrak r$}$$

where 𝑚 𝑑GLYPH&lt;150&gt; ' 𝑦 ' represents the y-disparity value at the position » 𝑑GLYPH&lt;150&gt; 𝑦 … &gt; and 𝜆 is a positive smoothness term [21]. 𝐸 1 represents the energy of a possible target path in the y-disparity map. 𝜏 min is typically set to 0. 𝜏 max depends entirely on 𝜶 , and it is set to 10 in this paper. The target path 𝑀 = f' 𝑑 GLYPH&lt;150&gt; 𝑖 𝑦 𝑖 ' GLYPH&lt;150&gt; 𝑖 = 0 1 GLYPH&lt;150&gt; GLYPH&lt;150&gt; GLYPH&lt;149&gt; GLYPH&lt;149&gt; GLYPH&lt;149&gt; GLYPH&lt;150&gt; 𝑛 g can be found by minimizing the energy function in (14), where ' 𝑑 GLYPH&lt;150&gt; 𝑖 𝑦 𝑖 ' 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 𝜶 . The disparity map can, therefore, be transformed using 𝜃 and 𝜶 . However, the outliers in the target path may greatly affect the accuracy of 𝜶 estimation. We, therefore, use RANSAC to update the values in 𝜶 . The full list of procedures involved for 𝜶 estimation are detailed in Algorithm 2.

1

is

2

is replaced by

replaced by

𝜃

3

𝜃

4

;

;

## Algorithm 2: 𝜶 estimation using DP.

Input

:

𝑀

Output: 𝜶

1 create a 𝑡 GLYPH&lt;2&gt; ' 𝑠 , 3 ' matrix 𝑻 ;

- 2 for 𝑖 1 to 𝑡 do

3

randomly select

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

𝑝

pairs of coordinates from estimate

𝜶

using Eq. 3;

𝑻

'

𝑖GLYPH&lt;150&gt;

𝑠

,

1 :

for

𝑠

,

3

'

𝜶

&gt;

𝑗

1

to

𝑠

do set the tolerance to

𝜀

𝑗

𝜶

GLYPH&lt;157&gt;

2

GLYPH&lt;0&gt;

1

;

compute the number of inliers

𝑛

compute the number of outliers inlier

;

𝑛

outlier compute the ratio

𝜂

=

𝑛

inlier

GLYPH&lt;157&gt;

𝑛

outlier

𝑻

'

𝑖GLYPH&lt;150&gt;

𝑗

'

𝜂

;

end end

𝑗

=

do

0

;

𝑗

,

𝑗

1

;

find the highest 𝜂 in the 𝑗 th column of 𝑻 ;

while the highest

𝜂

in the

𝑗

th column of

𝑻

corresponds to more than one 𝜶 ;

- 19 find 𝜶 which corresponds to the highest 𝜂 in the 𝑗 th column of 𝑻 ;

RANSAC is iterated 𝑡 times. Selecting a higher 𝑡 raises the possibility of finding the best 𝜶 but also increases the processing time. In order to minimize the trade-off between speed and robustness, 𝑡 is set to 50 in this paper. In each iteration, we select 𝑝 pairs of coordinates » 𝑑GLYPH&lt;150&gt; 𝑦 … &gt; from the target path to estimate 𝜶 . For a smaller 𝑝 , there is less chance that any outliers will influence optimization. In this paper 𝑝 is set to 3 , which is the smallest possible value for determining 𝜶 . The ratio 𝜂 of inliers versus outliers can then be computed with respect to a given tolerance 𝜀 𝜶 . The best 𝜶 corresponds to the highest ratio 𝜂 . However, the selection of an appropriate 𝜀 𝜶 can be challenging, as it is possible that there could be more than one satisfying value for 𝜶 . Hence, in this paper, the value of 𝜀 𝜶 also changes 𝑠 times. In each iteration, the value of 𝜀 𝜶 reduces by half and 𝜂 is computed. The best 𝜶 can be determined by finding the highest ratio 𝜂 . In this paper, the values of 𝜀 𝜶 and 𝑠 are both set to 4 .

Finally, each disparity is transformed using the following equation:

$$\int _ { \text{here} } ^ {., \, n _ { \text{J} } } \quad \tilde { d } ( u, v, d, \theta ) = [ 1 \ y ( u, v, \theta ) \ y ( u, v, \theta ) ^ { 2 } ] \alpha ( \theta ) - d + \delta, \ \ ( 1 5 )$$

where 𝛿 is a constant value set to guarantee that all the transformed disparity values are non-negative. In this paper, we set 𝛿 to 30 . The transformed disparity map is shown in Fig. 6(a). 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.

;

;

;

𝑀

;

Fig. 6. Disparity transformation and undamaged road area extraction: (a) transformed disparity map; (b) extracted undamaged road areas.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-7ed877f4-3159-4a5f-a572-c93444e81e12.jpg)

## B. Undamaged Road Area Extraction

Next, we utilize Otsu's thresholding method to segment the transformed disparity map. The segmentation threshold 𝑇 o can be obtained by maximizing the inter-class variance 𝜎 2 o as follows [22]:

$$\sigma _ { o } ^ { 2 } ( T _ { o } ) = P$$

where

$$\intertext {$$

represent the probabilities of damaged and undamaged road areas, respectively. 𝑝 𝑖 ' ' is the probability of ˜ 𝑑 = 𝑖 . The average disparity values of the damaged and undamaged road areas are given by:

$$\cdot \\ \mu _ { 0 } ( T _ { o } ) = \frac$$

The segmentation result is shown in Fig. 6(b). 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.

## 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]:

$$X ^ { \mathcal { W } } = \frac { u T _ { c }$$

where 𝑓 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.

Fig. 7. Surface normal vectors mapping on a sphere.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-91abaa30-59b9-4434-b1d9-cfb5e3ccc999.jpg)

To model the disparity map, 𝒄 = » 𝑐 0 GLYPH&lt;150&gt; 𝑐 1 GLYPH&lt;150&gt; 𝑐 2 GLYPH&lt;150&gt; 𝑐 3 GLYPH&lt;150&gt; 𝑐 4 GLYPH&lt;150&gt; 𝑐 5 … &gt; , storing the quadratic surface model coefficients can be estimated as follows:

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

where

$$W = \begin{bmatrix} 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 } \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & u _ { n } & v _ { n } & u _ { n } ^ { 2 } & v _ { n } ^ { 2 } & u _ { n } v _ { n } \end{bmatrix}. \quad \begin{matrix} ( 2 1 ) \\ \end{matrix} \\. \end{bmatrix}.$$

4

5

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:

- GLYPH&lt;15&gt; it is located in one of the damaged road areas.
- GLYPH&lt;15&gt; its surface normal vector differs greatly from the optimal one.
- GLYPH&lt;15&gt; 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.

1) Optimal Normal Vector Estimation: For each point 𝒑 𝑖 = » 𝑢 GLYPH&lt;150&gt; 𝑖 𝑣 GLYPH&lt;150&gt; 𝑖 𝑑 𝑖 … &gt; in the undamaged road area, we would like to estimate a normal vector 𝒏 𝑖 = » 𝑛 𝑢𝑖 GLYPH&lt;150&gt; 𝑛 𝑣𝑖 GLYPH&lt;150&gt; 𝑛 𝑑𝑖 … &gt; from a set of 𝑘 points in its neighborhood 𝑸 𝑖 = » 𝒒 𝑖 1 GLYPH&lt;150&gt; 𝒒 𝑖 2 GLYPH&lt;150&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;150&gt; 𝒒 𝑖𝑘 … &gt; , where 𝒒 𝑖 𝑗 ≠ 𝒑 𝑖 . Here, we define the augmented neighbor matrix 𝑸 , 𝑖 which contains all neighbors and the point 𝒑 𝑖 itself as follows:

$$Q _ { i } ^ { + } = [ p _ { i }, \, Q _ { i } ^ { \top } ] ^ { \top }.$$

Existing normal vector estimation methods are generally classified into one of two categories: optimization-based and

Fig. 8. Disparity map modeling and pothole detection: (a) modeled disparity map, (b) detected potholes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-3da762b4-77a4-45df-94eb-47f65e5d2346.jpg)

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. 𝒏 𝑖 can be estimated as follows:

$$n _ { i } = \underset { n _ { i } } { \arg \min } \left \| \left [ Q _ { i } ^ { + } - \bar { Q } _ { i } ^ { + } \right ] n _ { i } \right \| _ { 2 }, \quad \quad ( 2 3 ) \quad \text{exec}$$

where

$$\bar { Q } _ { i } ^ { + } = \mathbf 1 _ { k + 1 } \left ( \frac { 1 } { k + 1 } Q _ { i } ^ { + \top } \mathbf 1 _ { k + 1 } \right ) ^ { \top }. \quad \quad \ ( 2 4 ) \, \stackrel { \text{case} } { \text{areas} } \, \\ \text{cotict} \,.$$

1 𝑚 represents an 𝑚 GLYPH&lt;2&gt; 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 ˆ 𝒏 can be determined by finding the position at which the projections distribute most intensively.

Since the projection of ˆ 𝒏 is also on the sphere, it can be written in spherical coordinates as follows:

$$\hat { n } = [ \sin \varphi _ { 1 } & \cos \varphi _ { 2 }, \ \sin \varphi _ { 1 } \sin \varphi _ { 2 }, \ \cos \varphi _ { 1 } ] ^ { \top } \\ s. t. \ \varphi _ { 1 } & \in [ 0, \pi ], \ \varphi _ { 2 } \in [ 0, 2 \pi ). \quad \text{ (25)} \quad \text{ equal}. \\ \quad \text{ and see}.$$

It can be estimated by minimizing 𝐸 2 :

$$E _ { 2 } = { \mathbf 1 } _ { n + 1 } ^ { \top } m, \quad \ \ \ \ ( 2 6 ) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where

$$m = [ - n _ { 0 } \cdot \hat { n }, \ - n _ { 1 } \cdot \hat { n }, \ \cdots, \ - n _ { n } \cdot \hat { n } ] ^ { \top }. \quad ( 2 7 ) \ \text{produ}$$

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

$$\tan \varphi _ { 1 } = \frac { 1 _ { n + 1 } ^ { \top } n _ { u } \cos \varphi _ { 2 } + 1 _ { n + 1 } ^ { \top } n _ { v } \sin \varphi _ { 2 } } { 1 _ { n + 1 } ^ { \top } n _ { d } }, \quad ( 2 8 ) \quad \text{higes} _ { \text{ } \text{Algori} }$$

$$\tan \varphi _ { 2 } = \frac { 1 _ { n + 1 } ^ { \top } n _ { v } } { 1 _ { n + 1 } ^ { \top } n _ { u } }, \quad \quad ( 2 9 ) \, \frac { \sin } { 8 ( a ) }$$

where 𝒏 𝑢 = » 𝑛 𝑢 0 GLYPH&lt;150&gt; 𝑛 𝑢 1 GLYPH&lt;150&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;150&gt; 𝑛 𝑢𝑛 … &gt; , 𝒏 𝑣 = » 𝑛 𝑣 0 GLYPH&lt;150&gt; 𝑛 𝑣 1 GLYPH&lt;150&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;150&gt; 𝑛 𝑣𝑛 … &gt; and 𝒏 𝑑 = » 𝑛 𝑑 0 GLYPH&lt;150&gt; 𝑛 𝑑 1 GLYPH&lt;150&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;150&gt; 𝑛 𝑑𝑛 … &gt; . Due to the fact that 𝜑 2 is between 0 and 2 𝜋 , (29) will have two solutions:

$$\varphi _ { 2 } = \arctan \frac { 1 _ { n + 1 } ^ { \top } n _ { v } } { 1 _ { n + 1 } ^ { \top } n _ { u } } + k \pi, \ k \in \{ 0, 1 \}. \quad \ ( 3 0 ) \quad \text{becau} \colon$$

Substituting each 𝜑 2 into (28) produces a value for 𝜑 1 . The two pairs of » 𝜑 GLYPH&lt;150&gt; 1 𝜑 2 … &gt; correspond to the maxima and minima of 𝐸 2 , respectively. By substituting each pair of » 𝜑 GLYPH&lt;150&gt; 1 𝜑 2 … &gt;

Fig. 9. The point clouds of the detected potholes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-0a71ae56-70ea-47b2-9603-d292056c5bf9.jpg)

into Eq. 26 and comparing the two obtained values, we can find the optimal normal vector. If the angle between ˆ 𝒏 and 𝒏 𝑖 exceeds a pre-set threshold 𝜀 𝒏 , 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, 𝜀 𝒏 is set to 𝜋 GLYPH&lt;157&gt; 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.

2) 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 𝑡 times. In each iteration, a subset of disparities are selected randomly to estimate 𝒄 . 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 𝑟 GLYPH&lt;2&gt; 𝑟 . As 𝑟 becomes smaller, more disparities will be used for surface fitting, which increases the computational complexity. In contrast, selecting a higher value for 𝑟 results in less computational complexity, but potentially increases noise sensitivity. In this paper, the value of 𝑟 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 𝜂 between the inliers and outliers are obtained. 𝒄 which corresponds to the highest 𝜂 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(a).

3) 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 𝜀 𝑑 . Before labeling different potholes using CCL, the connected components containing fewer than 𝑤 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 𝜀 𝑑 and 𝑤 is discussed in Section III-D. The detected potholes are shown in Fig. 8(b). Finally, the point clouds of the detected

Fig. 10. Experimental set-up for acquiring stereo road images.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-deaa2a3c-2a3a-40eb-aa8e-f10ed8f31b87.jpg)

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.

## A. Experimental Set-Up

In this work, we utilized a ZED stereo camera 1 to capture stereo road images. An example of the experimental setup 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 GLYPH&lt;2&gt; 1730 , 1030 GLYPH&lt;2&gt; 1720 , 1028 GLYPH&lt;2&gt; 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.

## 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 'GLYPH&lt;0&gt; 𝜋 2 GLYPH&lt;150&gt; 𝜋 2 … to find the local minima. Therefore, the computational complexity is O' 𝜋 𝜀 𝜃 ' . In our method, GSS reduces the interval size exponentially. As a result, the interval size then becomes 𝜅 𝑛 𝜋 after the 𝑛 -th iteration. Therefore, the proposed roll angle estimation algorithm reduces the computational complexity to O' log 𝜅 𝜀 𝜃 𝜋 ' . The proposed roll angle estimation algorithm needs 21 iterations to produce a roll angle, with an accuracy higher than 𝜋 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

1 https://www.stereolabs.com/

Fig. 11. Experimental results of the EISATS synthesized stereo dataset: (a) left stereo images (the areas in magenta are our manually selected road areas); (b) ground truth disparity maps; (c) transformed disparity maps.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-7b9bf6eb-6494-429d-9a25-1259bf0d849b.jpg)

utilized to estimate the roll angle ˆ 𝜃 . The absolute difference between the actual and estimated roll angles, i.e., Δ 𝜃 = j 𝜃 GLYPH&lt;0&gt; ˆ 𝜃 j , is computed for each frame. The average Δ 𝜃 is approximately 1 129 GLYPH&lt;149&gt; GLYPH&lt;2&gt; 10 GLYPH&lt;0&gt; 4 rad which is lower than 𝜋 18000 rad. Therefore, the proposed algorithm is capable of estimating the roll angle with high accuracy.

## 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 𝜎 𝑑 of the transformed disparity values as follows:

$$\sigma _ { d } = \sqrt { \frac { 1 } { m + 1 } \left \| \tilde { d } - \frac { \tilde { d } ^ { \top } \mathbf 1 _ { m + 1 } } { m + 1 } \right \| _ { 2 } }. \quad \quad ( 3 1 )$$

where ˜ 𝒅 = » ˜ 𝑑 GLYPH&lt;150&gt; 0 ˜ 𝑑 GLYPH&lt;150&gt; 1 GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;1&gt; GLYPH&lt;150&gt; ˜ 𝑑 𝑚 … &gt; is a column vector storing the transformed disparity values. The average 𝜎 𝑑 values of the two KITTI stereo datasets are provided in Table I.

TABLE I COMPARISON OF THE AVERAGE 𝜎 𝑑 VALUES BETWEEN THE TWO CASES.

| Dataset                            |   Case 1 |   Case 2 |
|------------------------------------|----------|----------|
| KITTI stereo 2012 training dataset |   0.4862 |   0.78   |
| KITTI stereo 2015 training dataset |   0.5506 |   0.9428 |

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

Fig. 12. Experimental results of the KITTI stereo datasets: (a) left stereo images (the areas in purple are the labeled road regions); (b) ground truth disparity maps; (c) transformed disparity maps.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-50b05d35-2b0f-46dc-92ab-1cc2ea17bd1f.jpg)

Fig. 13. The sum of Δ 𝑛 PD with respect to different 𝜀 𝑑 and 𝑤 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-9d7bcd2b-810f-4757-8b8e-8dadc881a4ac.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-ba2411d6-75b5-4650-80c0-d69e4925c691.jpg)

zero. From Table I, we can observe that the average 𝜎 𝑑 values of these two datasets reduce by approximately half when the effects caused by the non-zero roll angle are eliminated. The average 𝜎 𝑑 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.

## 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 𝑤 pixels and the disparity difference of each pixel exceeds 𝜀 𝑑 , it will be identified as a pothole. In our experiments, we utilize the brute-force search method to find the best values of 𝜀 𝑑 and 𝑤 . The search range for 𝜀 𝑑 and 𝑤 are set to » 3 0 GLYPH&lt;149&gt; GLYPH&lt;150&gt; 8 5 GLYPH&lt;149&gt; … and » 100 5000 GLYPH&lt;150&gt; … , respectively. The step sizes for searching 𝜀 𝑑 and 𝑤 are set to 0 1 GLYPH&lt;149&gt; and 100 , respectively.

For the first step, we go through the whole search range and record the number of detected potholes ˆ 𝑛 PD in each frame. The absolute difference Δ 𝑛 PD between each ˆ 𝑛 PD and the expected pothole number 𝑛 PD is then computed. The sum of Δ 𝑛 PD with respect to a pair of given 𝜀 𝑑 and 𝑤 can therefore be obtained, as illustrated in Fig. 13. In our experiments, the least sum of Δ 𝑛 PD is equal to one and it is achieved only when 𝜀 𝑑 = 6 2 GLYPH&lt;149&gt; and 𝜀 𝑑 = 3100 . The corresponding incorrect pothole detection

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-7918038e-2085-446f-a5ad-1c048c088b8f.jpg)

(a)

Fig. 14. Experimental result of incorrect pothole detection: (a) left stereo image; (b) transformed disparity map; (c) detection result; (d) ground truth.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-6169cedf-3900-41fb-b8bd-5baf43b02f14.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-866f63d9-43a6-4641-a373-ca5fb39add86.jpg)

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% GLYPH&lt;149&gt; and 84 8% GLYPH&lt;149&gt; , respectively. Compared to them, our proposed algorithm can detect potholes more accurately with a successful detection accuracy of 98 7% GLYPH&lt;149&gt; .

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

$$\text{precision} = \frac { n _ { \text{TP} } } { n _ { \text{TP} } + n _ { \text{FP} } },$$

$$\text{recall} = \frac { n _ { \text{TP} } } { n _ { \text{TP} } + n _ { \text{FN} } },$$

Fig. 15. Experimental results of pothole detection: (a) left stereo images; (b) transformed disparity maps; (c) pothole detection results obtained using the algorithm proposed in [15]; (d) pothole detection results obtained using the algorithm presented in [16]; (e) pothole detection results obtained using the proposed algorithm; (f) pothole ground truth.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00894v3-figure-8b335c40-d505-49ee-8693-9a3ba6954b69.jpg)

TABLE II COMPARISON OF SUCCESSFUL POTHOLE DETECTION ACCURACY.

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

$$F \text{-score} = 2 \times \frac { \text{precision} \times \text{recall} } { \text{precision} + \text{recall} }, \quad \text{ (34)} \quad \text{false} \quad \text{--end}.$$

false negative and true negative pixel numbers, respectively.

$$\text{accuracy} = \frac { n _ { \text{TP} } + n _ { \text{TN} } } { n _ { \text{TP} } + n _ { \text{TN} } + n _ { \text{FP} } + n _ { \text{FN} } }, \quad \quad ( 3 5 ) \quad \text{illustr} \\ \text{trithm}$$

where 𝑛 TP, 𝑛 FP, 𝑛 FN and 𝑛 TN are true positive, false positive,

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.

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 |
| Dataset 1 | [16]     |   0.4622 |      0.9976 |    0.6317 |     0.9936 |
| Dataset 1 | proposed |   0.499  |      0.9871 |    0.6629 |     0.994  |
| Dataset 2 | [15]     |   0.9754 |      0.9712 |    0.9733 |     0.9987 |
| Dataset 2 | [16]     |   0.8736 |      0.9907 |    0.9285 |     0.9968 |
| Dataset 2 | proposed |   0.9804 |      0.9797 |    0.98   |     0.9991 |
| Dataset 3 | [15]     |   0.6119 |      0.7714 |    0.6825 |     0.9948 |
| Dataset 3 | [16]     |   0.5339 |      0.992  |    0.6942 |     0.9957 |
| Dataset 3 | proposed |   0.5819 |      0.9829 |    0.731  |     0.9961 |
| Total     | [15]     |   0.7799 |      0.822  |    0.8004 |     0.9942 |
| Total     | [16]     |   0.6948 |      0.9921 |    0.8173 |     0.9954 |
| Total     | proposed |   0.7709 |      0.9815 |    0.8635 |     0.9964 |

## 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% GLYPH&lt;149&gt; and the pixel-level accuracy is approximately 99 6% GLYPH&lt;149&gt; .

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.

## REFERENCES

- [1] J. S. Miller and W. Y. Bellinger, 'Distress identification manual for the long-term pavement performance program,' United States. Federal Highway Administration. Office of Infrastructure Research and Development, Tech. Rep., 2014.
- [2] S. Mathavan, K. Kamal, and M. Rahman, 'A review of threedimensional imaging technologies for pavement distress detection and measurements,' IEEE Transactions on Intelligent Transportation Systems , vol. 16, no. 5, pp. 2353-2362, 2015.
- [3] T. Kim and S.-K. Ryu, 'Review and analysis of pothole detection methods,' Journal of Emerging Trends in Computing and Information Sciences , vol. 5, no. 8, pp. 603-608, 2014.
- [4] R. Fan, J. Jiao, J. Pan, H. Huang, S. Shen, and M. Liu, 'Real-time dense stereo embedded in a uav for road inspection,' in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops , 2019, pp. 0-0.
- [5] C. Koch, K. Georgieva, V. Kasireddy, B. Akinci, and P. Fieguth, 'A review on computer vision based defect detection and condition assessment of concrete and asphalt civil infrastructure,' Advanced Engineering Informatics , vol. 29, no. 2, pp. 196-210, 2015.
- [6] Y.-C. Tsai and A. Chatterjee, 'Pothole detection and classification using 3d technology and watershed method,' Journal of Computing in Civil Engineering , vol. 32, no. 2, p. 04017078, 2017.
- [7] I. Pitas, Digital image processing algorithms and applications . John Wiley &amp; Sons, 2000.
- [8] S. Li, C. Yuan, D. Liu, and H. Cai, 'Integrated processing of image and gpr data for automated pothole detection,' Journal of computing in civil engineering , vol. 30, no. 6, p. 04016015, 2016.
- [9] E. Buza, S. Omanovic, and A. Huseinovic, 'Pothole detection with image processing and spectral clustering,' in Proceedings of the 2nd International Conference on Information Technology and Computer Networks , 2013, pp. 48-53.
- [10] C. Koch and I. Brilakis, 'Pothole detection in asphalt pavement images,' Advanced Engineering Informatics , vol. 25, no. 3, pp. 507-515, 2011.
- [11] M. R. Jahanshahi, F. Jazizadeh, S. F. Masri, and B. Becerik-Gerber, 'Unsupervised approach for autonomous pavement-defect detection and quantification using an inexpensive depth sensor,' Journal of Computing in Civil Engineering , vol. 27, no. 6, pp. 743-754, 2012.
- [12] R. Fan, X. Ai, and N. Dahnoun, 'Road surface 3D reconstruction based on dense subpixel disparity map estimation,' IEEE Transactions on Image Processing , vol. PP, no. 99, p. 1, 2018.
- [13] C. Zhang and A. Elaksher, 'An unmanned aerial vehicle-based imaging system for 3d measurement of unpaved road surface distresses,' Computer-Aided Civil and Infrastructure Engineering , vol. 27, no. 2, pp. 118-129, 2012.
- [14] R. Fan, 'Real-time computer stereo vision for automotive applications,' Ph.D. dissertation, University of Bristol, 2018.
- [15] Z. Zhang, X. Ai, C. Chan, and N. Dahnoun, 'An efficient algorithm for pothole detection using stereo vision,' in Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on . IEEE, 2014, pp. 564-568.
- [16] A. Mikhailiuk and N. Dahnoun, 'Real-time pothole detection on tms320c6678 dsp,' in Imaging Systems and Techniques (IST), 2016 IEEE International Conference on . IEEE, 2016, pp. 123-128.
- [17] L. Cruz, L. Djalma, and V. Luiz, 'Kinect and rgbd images: Challenges and applications graphics,' in 2012 25th SIBGRAPI Conference on Patterns and Images Tutorials (SIBGRAPI-T) , 2012.
- [18] R. Hartley and A. Zisserman, Multiple view geometry in computer vision . Cambridge university press, 2003.
- [19] M. A. Fischler and R. C. Bolles, 'Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,' Communications of the ACM , vol. 24, no. 6, pp. 381-395, 1981.
- [20] P. Pedregal, Introduction to optimization . Springer Science &amp; Business Media, 2006, vol. 46.
- [21] U. Ozgunalp, R. Fan, X. Ai, and N. Dahnoun, 'Multiple lane detection algorithm based on novel dense vanishing point estimation,' IEEE Transactions on Intelligent Transportation Systems , vol. 18, no. 3, pp. 621-632, 2017.
- [22] N. Otsu, 'A threshold selection method from gray-level histograms,' IEEE transactions on systems, man, and cybernetics , vol. 9, no. 1, pp. 62-66, 1979.
- [23] M. B. Dillencourt, H. Samet, and M. Tamminen, 'A general approach to connected-component labeling for arbitrary image representations,' Journal of the ACM (JACM) , vol. 39, no. 2, pp. 253-280, 1992.
- [24] H. Eric Tseng, L. Xu, and D. Hrovat, 'Estimation of land vehicle roll and pitch angles,' Vehicle System Dynamics , vol. 45, no. 5, pp. 433-443, 2007.
- [25] J. Oh and S. B. Choi, 'Vehicle roll and pitch angle estimation using a cost-effective six-dimensional inertial measurement unit,' Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering , vol. 227, no. 4, pp. 577-590, 2013.
- [26] J. Ryu, E. J. Rossetter, and J. C. Gerdes, 'Vehicle sideslip and roll parameter estimation using gps,' in Proceedings of the AVEC International Symposium on Advanced Vehicle Control , 2002, pp. 373-380.

- [27] J. Ryu and J. C. Gerdes, 'Estimation of vehicle roll and road bank angle,' in American Control Conference, 2004. Proceedings of the 2004 , vol. 3. IEEE, 2004, pp. 2110-2115.
- [28] R. Labayrade and D. Aubert, 'A single framework for vehicle roll, pitch, yaw estimation and obstacles detection by stereovision,' in Intelligent Vehicles Symposium, 2003. Proceedings. IEEE . IEEE, 2003, pp. 31-36.
- [29] P. Skulimowski, M. Owczarek, and P. Strumiłło, 'Ground plane detection in 3d scenes for an arbitrary camera roll rotation through 'vdisparity' representation,' in 2017 Federated Conference on Computer Science and Information Systems (FedCSIS) . IEEE, 2017, pp. 669-674.
- [30] M. Evans, R. Fan, and N. Dahnoun, 'Iterative roll angle estimation from dense disparity map,' in Proc. 7th Mediterranean Conf. Embedded Computing (MECO) , Jun. 2018, pp. 1-4.
- [31] R. Fan, M. J. Bocus, and N. Dahnoun, 'A novel disparity transformation algorithm for road segmentation,' Information Processing Letters , vol. 140, pp. 18-24, 2018.
- [32] I. Haller and S. Nedevschi, 'Design of interpolation functions for subpixel-accuracy stereo-vision systems,' IEEE Transactions on image processing , vol. 21, no. 2, pp. 889-898, 2012.
- [33] K. Klasing, D. Althoff, D. Wollherr, and M. Buss, 'Comparison of surface normal estimation methods for range sensing applications,' in Proc. IEEE Int. Conf. Robotics and Automation , May 2009, pp. 32063211.
- [34] C. Wang, H. Tanahashi, H. Hirayu, Y. Niwa, and K. Yamamoto, 'Comparison of local plane fitting methods for range data,' in Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001 , vol. 1. IEEE, 2001, pp. I-I.
- [35] T. Vaudrey, C. Rabe, R. Klette, and J. Milburn, 'Differences between stereo and motion behaviour on synthetic and real-world stereo sequences,' in Image and Vision Computing New Zealand, 2008. IVCNZ 2008. 23rd International Conference . IEEE, 2008, pp. 1-6.
- [36] A. Wedel, C. Rabe, T. Vaudrey, T. Brox, U. Franke, and D. Cremers, 'Efficient dense scene flow from sparse or dense stereo data,' in European conference on computer vision . Springer, 2008, pp. 739751.
- [37] A. Geiger, P. Lenz, and R. Urtasun, 'Are we ready for autonomous driving? the kitti vision benchmark suite,' in Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on . IEEE, 2012, pp. 3354-3361.
- [38] M. Menze, C. Heipke, and A. Geiger, 'Joint 3d estimation of vehicles and scene flow,' ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences , vol. 2, p. 427, 2015.
Rui Fan Umar Ozgunalp Brett Hosking Ming Liu Ioannis Pitas
Pothole Detection Based on Disparity Transformation and Road Surface Modeling
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%.