Real-Time Runway Detection for Infrared Aerial Image Using Synthetic Vision and an ROI Based Level Set Method

Abstract

We present a new method for real-time runway detection embedded in synthetic vision and an ROI (Region of Interest) based level set method. A virtual runway from synthetic vision provides a rough region of an infrared runway. A three-thresholding segmentation is proposed following Otsu’s binarization method to extract a runway subset from this region, which is used to construct an initial level set function. The virtual runway also gives a reference area of the actual runway in an infrared image, which helps us design a stopping criterion for the level set method. In order to meet the needs of real-time processing, the ROI based level set evolution framework is implemented in this paper. Experimental results show that the proposed algorithm is efficient and accurate.

1. Introduction

Low visibility is a primary cause of flight accidents and disruptions in flight scheduling [1,2]. Within a limited visual range and under low visibility conditions, pilots are likely to operate and potentially make errors in their judgement. The instrument landing system (ILS) [3,4] is widely used to assist pilots’ visual contact with a runway during landing, taxiing and takeoff. However, it is not suitable for all airports and runways due to its demanding installation criteria. Under these circumstances, enhanced flight vision system (EFVS) and synthetic vision system (SVS) [5] undoubtedly complement the limits of natural vision. An EFVS is a system that provides pilots with an image typically on a heads-up display (HUD) from multiple sensors, such as charge-coupled device (CCD) cameras, infrared cameras or radar. An SVS is a computer mediated reality system that uses 3D scene topography to provide pilots with a rendered image. This image simulates what pilots actually see in that position and pose and is derived from airborne global positioning system (GPS) and attitude-heading reference system (AHRS). In this paper, we call such a rendered image the “virtual image.” An EFVS and a SVS may be combined to create fused information to significantly improve pilots’ visual perception [6].

Runway intrusion is an important aviation safety concern [7,8]. Therefore, runway detection can help mitigate these risks. Runways are regarded as a general landmark in the course of landing and take-off. In published literature, runway detection methods are divided into two categories: line segments based [9,10,11,12] and feature learning based runway detection. In remote sensing images, straight line segments are used as shape primitives to extract runways. In more concrete terms, image edges are first detected, followed by line segments extraction, such as Hough transform, line segmentation detector [13,14] or EDLines [15]. Parallelity is imposed upon line segments to determine candidate runways. Geometric shape is the distinct trait that existing methods use. This approach is similar to validation in remote sensing images. However, under low visibility in infrared images, chances are that edge of runway cannot be extracted. Furthermore, in the phase of approaching and landing, runways cannot be distinguished from line segments, distracted by highway, buildings and fields. As a result, it is difficult to exclude false detections from runway candidates. Feature learning based airport detection algorithms have been proposed in [16,17,18,19,20,21], in which pairs of features were extracted. Then, the Adaboost or an unsupervised learning algorithm was used to identify runways. Convolutional neural networks [22,23] have a strong feature representation power, which is good for object detection. The feature learning method is time-consuming and can only coarsely localize runways, which is not suitable for real-time processing. Simultaneously, a vast sample set is necessary, which confines the feature learning method in applications of runway detection for aerial images under low visibilities.

EFVS can integrate data acquired by weather penetrating sensors to overcome the problems incurred under adverse weather conditions. In combination with line detector and multi-sensor data fusion, runways can be recognized [24,25]. SVS can generate virtual images which roughly match the corresponding IR (Infrared Radiation) images. In comparison with the EFVS based method, SVS is provided with prior knowledge for runway extraction. Thus, the location of a runway in a virtual image, namely the virtual runway, can be considered as prior knowledge. After segmentation, we are able to acquire a subset of the IR runway, which serves as an initial contour. Considering the level set method, introduced by Osher et al. [26], we are able to detect weak edges; this is taken into account in our real-time Runway Detection framework. Numerical implementations of a level set method are complicated and heavily iterative. Li et al. [27] introduced a new functional to constrain level set functions to approach a signed distance function. As a result, re-initialization can be eliminated, which reduces the computation load. Meanwhile, some researchers focused on the fast marching level set method [28,29], which uses only points close to the curve, namely narrow bands, at every time step to decrease the computational cost of the standard level set method. It greatly simplifies the iteration process.

Through repeated experiments on infrared runways, the touch zones lie so close to runways (see Figure 1) that it may be difficult to divide the two areas without an edge prior. Li’s method [27] defined an edge indicator function to make the dynamic curves move toward the object boundaries. Further discussion on this issue will be presented in Section 4. Li’s method makes it possible to differentiate touch zones and runways. Using a level set method to detect runways requires three issues to be addressed: initial contour, real-time processing and a stopping criterion. In this paper, a virtual runway induced by synthetic vision can provide an initial contour and a stopping criterion for the level set method. Incorporating synthetic vision into level set methods, we propose a real-time runway detection method. The main ideas include:

(1)

Locating a virtual runway using heuristic search to find out a subset of runways from infrared images. We generalize Otsu’s thresholding method to trichotomize the region overlapping a virtual runway, one of which is used as the initial contour in a level set method.

(2)

An ROI based level set method is proposed, which makes it possible to detect a runway in real time.

(3)

Shape of the virtual runway provides a rough area of the actual runway, which imposes a stopping strategy for evolution in the level set method.

2. Generation of Virtual Runway

Generation of the virtual runway is shown in Figure 2. First, high-definition texture maps are overlaid on to its digital elevation model (DEM) [30,31,32]. An architecture layout is used to create a high accuracy runway model with 3D professional modeling tools, like the Creator. Subsequently, terrain and runway local coordinate systems are transformed into a unified coordinate framework. Thereby, the terrain model and runway model are merged into an enhanced 3D terrain model. In addition, GPS offers longitude, latitude, and altitude to indicate the position of onboard cameras and AHRS offers heading, pitch and roll angle to indicate the pose of cameras. Finally, with the position and pose of onboard cameras three-dimensional engineering tools convert the six degrees of freedom data into a flat database format to generate a virtual image. The virtual image is intended to project objects from the enhanced 3D terrain model as they should appear in the correct perspective from a pilot’s position and orientation. The setup renders the enhanced 3D terrain model in response to incoming navigation sensor data. Since point coordinates of the 3D runway model, camera position and pose are known, the location of the virtual runway can be determined. The above-mentioned flow chart is demonstrated in Figure 2a. For a more visual explanation, procedures are also represented diagrammatically in Figure 2b, in which the region bounded by the quadrilateral in green is the virtual runway.

3. Real-Time Runway Detection Based on ROI Level Set Method

Li et al. [27] defined an improved level set method skipping re-initialization. The evolution equation of the level set function $\varphi$ proposed in [32] can be given by:

$$\left\{\begin{array}{c}\frac{\partial \varphi }{\partial t}=\mu \left[∆\varphi -\mathrm{div}\left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)\right]+\lambda \delta \left(\varphi \right)\mathrm{div}\left(g\frac{\nabla \varphi }{|\nabla \varphi |}\right)+\nu g\delta \left(\varphi \right),\hfill \\ {\varphi }^{n+1}={\varphi }^n+\tau \frac{\partial \varphi }{\partial t},n=0,1,\dots ,\hfill \end{array}\right.$$

(1)

where $\mu >0$ is a penalizing coefficient to prevent a level set function from deviating from a signed distance function, $\lambda >0$ and $\nu$ are parameters controlling the shape and size of the evolved curves, $\tau >0$ is an iteration step and g is an edge indicator function, given by:

$$g\triangleq \frac{1}{1+{\left|\nabla G_{\sigma }\ast I\right|}^2},$$

(2)

where $G_{\sigma }\ast I$ is the convolution of image I with the Gaussian kernel $G_{\sigma }$. The superscript n is the number of iterations, and ${\varphi }^0$ is the initial level set function. Empirically, if the initial contour defined by ${\varphi }^0$ is located outside the object of interest, v will take a positive value. Otherwise, a negative value is used to accelerate curve propagation towards actual contours of the object of interest. Using level set methods, there are three challenges in extracting runways from infrared images: (i) a proper initial contour for level set methods is necessary to detect a runway excluding touch zones and other connected regions; (ii) real-time processing is essential; and (iii) we need a stopping criterion to identify the exact runway. If a subset of the runway can be extracted, its boundary can naturally serve as an initial contour in a level set method. The histogram of the virtual runway bounded region (highlighted in green) in the infrared image is shown in Figure 3a, in which there are three or four peaks neglecting the effect of noise. Therefore, the traditional Otsu’s method cannot discriminate the touch zone (see the bottom part circle bounded in Figure 3b) from the runway. However, our method mentioned in Section 3.1 can properly partition the runway subset (in red) from the virtual runway bounded region.

3.1. Three-Thresholding Segmentation for Initial Contour

Runway structure in Figure 1 implies that the choice of initial level set function is a crucial factor leading to the accuracy and efficiency of runway detection. It is apparent that the virtual runway roughly outlines the region of IR runway. Specifically, the area where the virtual runway covers the IR image can be divided into three parts: background, touch zone and runway subregion according to the gray value difference. The brightest part is considered as a subset of the IR runway. Henceforth, the contour of such part serves as an initial contour in the level set method.

The classic Otsu method is practical and efficient for automatic thresholding [33,34]. In Otsu’s principle [35], the main goal is to maximize inter-class variance and minimize intra-class variance of classifications. In this paper, the corresponding region of the virtual runway in an IR image is regarded as a ROI, which can be represented as the convex hull of the vertices of the virtual runway $P_i(x,y),i=1,2,3,4$ with the mathematical formulation:

$$\mathrm{ROI}\triangleq \left\{{\displaystyle \sum _{i=1}^4}{\alpha }_i{P}_i|{\alpha }_i\ge 0,i=1,\dots ,4,{\displaystyle \sum _{i=1}^4}{\alpha }_i=1\right\}.$$

(3)

Let the pixels of the ROI be represented as $[g_{\min },\dots ,$ $i,\dots ,g_{\max }]$, where subscript min, max represents the minimum and maximum gray levels of the ROI, respectively. The number of pixels at level i is denoted by $n_i$ and the total number of pixels by $N={\displaystyle \sum _{i=g_{\min }}^{g_{\max }}}{n}_i$.

The gray level histogram is normalized as:

$$p_i=\frac{n_i}{N},p_i\ge 0,{\displaystyle \sum _{i=g_{\min }}^{g_{\max }}}{p}_i=1.$$

(4)

Assume that the three classes of the ROI are background denoted by $C_0$, touch zone denoted by $C_1$, and runway subset denoted by $C_2$. The pixels of the ROI can be trichotomized by two thresholds $k_1$ and $k_2$, with $g_{\min }\le k_1<k_2\le g_{\max }$. Motivated by the term “binary thresholding method” [35], the aforementioned process is named after the three-thresholding method. The probabilities of class occurrence are given by:

$$w_0=P_r\left(C_0\right)={\displaystyle \sum _{i=g_{\min }}^{k_1}}{p}_i=\omega \left(k_1\right),$$

(5)

$$w_1=P_r\left(C_1\right)={\displaystyle \sum _{i=k_1+1}^{k_2}}{p}_i=\omega \left(k_2\right)-\omega \left(k_1\right),$$

(6)

$$w_2=P_r\left(C_2\right)={\displaystyle \sum _{i=k_2+1}^{g_{\max }}}{p}_i=1-\omega \left(k_2\right).$$

(7)

Hence, the class mean levels are given by:

$$u_0={\displaystyle \sum _{i=g_{min}}^{k_1}}i{P}_r\left(i\right|C_0)={\displaystyle \sum _{i=g_{min}}^{k_1}}\frac{i{p}_i}{w_0}=\frac{\mu \left(k_1\right)}{\omega \left(k_1\right)},$$

(8)

where $\mu \left(k_1\right)={\sum }_{i=g_{min}}^{k_1}i{p}_i$ and

$$u_1={\displaystyle \sum _{i=k_1+1}^{k_2}}i{P}_r\left(i\right|C_1)={\displaystyle \sum _{i=k_1+1}^{k_2}}\frac{i{p}_i}{w_1}=\frac{\mu \left(k_2\right)-\mu \left(k_1\right)}{\omega \left(k_2\right)-\omega \left(k_1\right)},$$

(9)

$$u_2={\displaystyle \sum _{i=k_2+1}^{g_{\max }}}i{P}_r\left(i\right|C_2)={\displaystyle \sum _{i=k_2+1}^{g_{\max }}}\frac{i{p}_i}{w_2}=\frac{{\mu }_T-\mu \left(k_2\right)}{1-\omega \left(k_2\right)},$$

(10)

where

$${\mu }_T=\mu \left(g_{\max }\right)={\displaystyle \sum _{i=g_{\min }}^{g_{\max }}}i{p}_i.$$

(11)

Based on Equations (5)–(11), the following relations always hold for any choices of $k_1$ and $k_2$:

$${\displaystyle \sum _{i=0}^2}{w}_i=1,{\displaystyle \sum _{i=0}^2}{w}_i{u}_i={\mu }_T.$$

(12)

The class variances are represented as:

$${\sigma }_0^2={\displaystyle \sum _{i=g_{\min }}^{k_1}}{(i-u_0)}^2{P}_r\left(i\right|C_0)=\frac{{\sum }_{i=g_{\min }}^{k_1}{(i-u_0)}^2{p}_i}{w_0},$$

(13)

$${\sigma }_1^2={\displaystyle \sum _{i=k_1+1}^{k_2}}{(i-u_1)}^2{P}_r\left(i\right|C_1)=\frac{{\sum }_{i=k_1+1}^{k_2}{(i-u_1)}^2{p}_i}{w_1},$$

(14)

$${\sigma }_2^2={\displaystyle \sum _{i=k_2+1}^{g_{\max }}}{(i-u_2)}^2{P}_r\left(i\right|C_2)=\frac{{\sum }_{i=k_2+1}^{g_{\max }}{(i-u_2)}^2{p}_i}{w_2}.$$

(15)

In order to measure class separability, three class variances are introduced [36]. The within-class variance of levels is given by:

$${\sigma }_W^2=w_0{\sigma }_0^2+w_1{\sigma }_1^2+w_2{\sigma }_2^2.$$

(16)

The between-class variance of levels is given by:

$${\sigma }_B^2=w_0{(u_0-{\mu }_T)}^2+w_1{(u_1-{\mu }_T)}^2+w_2{(u_2-{\mu }_T)}^2.$$

(17)

The total variance of levels is given by:

$${\sigma }_T^2={\displaystyle \sum _{i=g_{\min }}^{g_{\max }}}{(i-{\mu }_T)}^2{p}_i.$$

(18)

Thus, the following discriminant criteria are commonly used:

$$\lambda =\frac{{\sigma }_B^2}{{\sigma }_W^2},\kappa =\frac{{\sigma }_T^2}{{\sigma }_W^2},\eta =\frac{{\sigma }_B^2}{{\sigma }_T^2}.$$

(19)

In order to discuss the relations of the indices mentioned above, we draw the following conclusion:

Theorem 1.

${\sigma }_T^2={\sigma }_W^2+{\sigma }_B^2,forany{k}_1,k_2.$

Proof. 

Substituting Equations (5)–(7), (13)–(15) into Equation (16) and expanding it, we obtain:

$$\begin{array}{cc}\hfill {\sigma }_W^2=& {\displaystyle {\sum }_{i=g_{\min }}^{k_1}}{(i-u_0)}^2{p}_i+{\displaystyle {\sum }_{i=k_1+1}^{k_2}}{(i-u_1)}^2{p}_i\hfill \\ & +{\displaystyle {\sum }_{i=k_2+1}^{g_{\max }}}{(i-u_2)}^2{p}_i\hfill \\ \\ \hfill =& {\displaystyle {\sum }_{i=g_{\min }}^{g_{\max }}}{i}^2{p}_i-\frac{{\mu }^2\left(k_1\right)}{\omega \left(k_1\right)}-\frac{{(\mu \left(k_2\right)-\mu \left(k_1\right))}^2}{\omega \left(k_2\right)-\omega \left(k_1\right)}-\frac{{({\mu }_T-\mu \left(k_2\right))}^2}{1-\omega \left(k_2\right)}.\hfill \end{array}$$

(20)

Similarly, substituting Equations (5)–(10) into Equation (17), we obtain:

$$\begin{array}{cc}\hfill {\sigma }_B^2=& \frac{{(\mu \left(k_1\right)-{\mu }_T\omega \left(k_1\right))}^2}{\omega \left(k_1\right)}+\frac{{\left[\mu \left(k_2\right)-\mu \left(k_1\right)-{\mu }_T(\omega \left(k_2\right)-\omega \left(k_1\right))\right]}^2}{\omega \left(k_2\right)-\omega \left(k_1\right)}\hfill \\ & +\frac{{\left[{\mu }_T-\mu \left(k_2\right)-{\mu }_T(1-\omega \left(k_2\right))\right]}^2}{1-\omega \left(k_2\right)}.\hfill \end{array}$$

(21)

Adding Equations (20) and (21) and simplifying, we get:

$${\sigma }_W^2+{\sigma }_B^2={\displaystyle \sum _{i=g_{\min }}^{g_{\max }}}{i}^2{p}_i-{\mu }_T^2.$$

(22)

On the other hand, expanding Equation (18) in connection with Equations (4) and (11), gives us:

$$\begin{array}{cc}\hfill {\sigma }_T^2=& {\displaystyle {\sum }_{i=g_{\min }}^{g_{\max }}}{i}^2{p}_i-2{\mu }_T{\displaystyle {\sum }_{i=g_{\min }}^{g_{\max }}}i{p}_i+{\mu }_T^2{\displaystyle {\sum }_{i=g_{\min }}^{g_{\max }}}{p}_i\hfill \\ \hfill =& {\displaystyle {\sum }_{i=g_{\min }}^{g_{\max }}}{i}^2{p}_i-{\mu }_T^2.\hfill \end{array}$$

(23)

Therefore, comparing the right-hand sides of Equations (22) and (23), we have ${\sigma }_T^2={\sigma }_W^2+{\sigma }_B^2$ as desired. ☐

Obviously, $\kappa$ and $\eta$ in Equation (19) can be recast in the form of $\lambda$:

$$\kappa =1+\lambda ;\eta =1-\frac{1}{\lambda +1}.$$

(24)

Therefore, the three discriminant criteria are equivalent to each other as follows:

$$max\lambda (k_1,k_2)\iff max\kappa (k_1,k_2)\iff max\eta (k_1,k_2).$$

Without loss of generality, the optimal value of $k_1$ and $k_2$ can be chosen as:

$$[k_1^{\ast },k_2^{\ast }]=argmax\left\{\lambda (k_1,k_2)\right\}.$$

(25)

Thus, the pixels of the ROI are trichomized into three classes:

$$g(x,y)\triangleq \left\{\begin{array}{cccc}255& f(x,y)\ge k_2,\hfill & \left(C_2\right),\hfill & \\ 128,& k_1\le f(x,y)<k_2,\hfill & \left(C_1\right),\hfill & (x,y)\in \mathrm{ROI},\hfill \\ 0& f(x,y)<k_1,\hfill & \left(C_0\right),\hfill & \end{array}\right.$$

(26)

where f is the grayscale intensity IR image and g the segmentation output.

We assume that the runway subset is the brightest part of the three classifications, namely the class $C_2$. Trace region boundaries in a binary image containing the $C_2$, and the boundary with the largest area is chosen as an initial contour in the level set method. Let the image domain be $\Omega \triangleq \left\{\right(x,y\left)\right|1\le x\le \mathrm{W},1\le y\le \mathrm{H}\}$, and the region bounded by the initial contour be ${\Omega }^0$. Thus, the initial level set function ${\varphi }^0$ is derived as:

$${\varphi }^0(x,y)\triangleq \left\{\begin{array}{cc}\hfill -\rho ,& (x,y)\in {\Omega }^0-\partial {\Omega }^0,\hfill \\ \hfill 0,& (x,y)\in \partial {\Omega }^0,\hfill \\ \hfill \rho ,& (x,y)\in \Omega -{\Omega }^0,\hfill \end{array}\right.$$

(27)

where $\rho >0$ is a predefined constant, and ∂ represents the boundary of the plane domain.

3.2. ROI Based Implementation for the Level Set Method

There is a product rule of divergence operator [37]: if $\phi$ is a scalar valued functional, and $\mathbf{F}$ is a vector field, then:

$$\mathrm{div}\left(\phi \mathbf{F}\right)=\nabla \phi ·\mathbf{F}+\phi \mathrm{div}\left(\mathbf{F}\right).$$

Applying the above rule to Equation (1), it will be recast as:

$$\left\{\begin{array}{cc}\hfill \frac{\partial \varphi }{\partial t}=& \mu \left[∆\varphi -\mathrm{div}\left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)\right]+\nu g{\delta }_{\epsilon }\left(\varphi \right),\hfill \\ & +\lambda {\delta }_{\epsilon }\left(\varphi \right)\left(\nabla g·\frac{\nabla \varphi }{|\nabla \varphi |}+g\mathrm{div}\left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)\right),\hfill \\ \hfill {\varphi }^{n+1}(x,y)=& {\varphi }^n(x,y)+\tau \frac{\partial \varphi }{\partial t},(x,y)\in \Omega ,n=0,1,\dots .\hfill \end{array}\right.$$

(28)

In numerical analysis, the $\delta \left(x\right)$ in Equation (1) is smoothed as:

$${\delta }_{\epsilon }\left(x\right)\triangleq \left\{\begin{array}{cc}\hfill 0,& \left|x\right|>\epsilon ,\hfill \\ \hfill \frac{1}{2\epsilon }[1+cos\left(\frac{\pi x}{\epsilon }\right)],& \left|x\right|\le \epsilon .\hfill \end{array}\right.$$

For ease of description, denote $∆\varphi \triangleq L$, $\frac{\nabla \varphi }{|\nabla \varphi |}\triangleq \mathbf{N}$, $\mathrm{div}\left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)=\mathrm{div}\left(\mathbf{N}\right)\triangleq K$. Note that, in Equation (28), the terms g and $\nabla g$ do not vary during the iterative procedure. However, the terms L, $\mathbf{N}$ and K are updated in every iteration. Ref. [38] observed that the front can be moved by updating the level set function at a small set of points in the neighborhood of the zero set instead of updating it at all the points on the grid. Hence, the narrow band based evolution method was conceived. Unfortunately, how to identify whether the points fall in the narrow band and how to retrieve the pre-computed values like g and $\nabla g$ promptly is a challenging problem. Considering this issue, we propose a rectangular ROI to implement level set methods without re-initialization. Denoting the ROI in the n-th iteration by ${\Omega }^n\subset \Omega$ and the discretized level set function ${\varphi }_{i,j}^n\triangleq \varphi (ih,jh,n\tau )$, where a grid in the domain $\Omega$ with spacing h and the time with spacing $\tau$, we can rewrite Equation (28) as follows:

$$\left\{\begin{array}{cc}\hfill {\varphi }_{i,j}^{n+1}=& {\varphi }_{i,j}^n+\tau \left[\mu (L_{i,j}^n-K_{i,j}^n)+\nu g_{i,j}{\delta }_{\epsilon }\left({\varphi }_{i,j}^n\right)\right.\hfill \\ & +\left.\lambda {\delta }_{\epsilon }\left({\varphi }_{i,j}^n\right)\left({\left(\nabla g\right)}_{i,j}·{\mathbf{N}}_{i,j}^n+g_{i,j}{K}_{i,j}^n\right)\right],\hfill \\ & (ih,jh)\in {\Omega }^n,\hfill \\ \hfill {\varphi }_{i,j}^0=& \varphi (ih,jh,0),(ih,jh)\in {\Omega }^0.\hfill \end{array}\right.$$

(29)

After n iterative evolutions, the region in which an object of interest lies is depicted by $A^n\triangleq \left\{(ih,jh)|{\varphi }_{i,j}^n\le 0\right\}$, whose minimal up-right bounding rectangle is given by:

$$B^n\triangleq \left\{(x,y)|x_{\min }^n\le x\le x_{\max }^n,y_{\min }^n\le y\le y_{\max }^n\right\}.$$

Assume that the current ROI can be expanded by $∆x$ and $∆y$ in the x- and y-directions, respectively. Hence, the ROI in the following iteration, i.e., ${\Omega }^{n+1}$ is determined by:

$${\Omega }^{n+1}\triangleq \left\{\begin{array}{cc}(x,y)|\hfill & max(x_{\min }^n-∆x,1)\le x\le min(x_{\max }^n+∆x,\mathrm{W}),\hfill \\ & max(y_{\min }^n-∆y,1)\le y\le min(y_{\max }^n+∆y,\mathrm{H}).\hfill \end{array}\right\}$$

3.3. Stopping Criteria

It is difficult to guarantee that when the computation stops the final shape is exactly the desired runway. We have designed two strategies. One is dependent on the area of the region $A^n$, once the area is close to the area of the virtual runway, the computation comes to a halt. The other one is to measure the variation of the last two level set functions. If all the points on the front come to a stop, it brings the computation to an end. If the values of the specific points of the last two level set functions satisfy:

$$\frac{{\sum }_{i,j}|{\varphi }_{i,j}^n-{\varphi }_{i,j}^{n-1}|}{|A^n|}\le \alpha \tau ,\mathrm{for}(ih,jh)\in A^n,$$

(30)

where $\alpha <1$ is a tolerance parameter related to time step $\tau$, it will stop the computations, and $A^n$ will be the desired final runway.

In this paper, we have a referenced area from the virtual runway. Thus, we choose the first strategy to stop. More specifically, once the current area of $A^n$ approximates the area for the virtual runway, iteration stops.

4. Analysis of Results and Comparisons

The algorithm was programmed using the Microsoft Visual Studio development platform (version 9.0, Redmond, WA, U.S.) with an Open CV (Computer Vision) library. How to select the related parameters is dependent on the processed objects. The authors [27] suggest $\tau \mu <\frac{1}{4}$ so as to maintain stable level set evolution. Commonly, p can take any positive value, $\sigma \in [1.0,4.0]$ and $\epsilon ,h,\delta$ are independent of images to some extent. Except for the parameters, considering the stability and time performance, through multiple trials, the relevant parameters for all the experiments in this paper are listed in

Table 1. Parameters related to experiments.

Image Dimension	$\mathit{\sigma }$	$\mathit{\rho }$	$\mathit{\epsilon }$	$\mathit{\mu }$	$\mathit{\nu }$	$\mathit{\lambda }$	$\mathit{\tau }$	h	∆x	∆y
720 × 576	1.5	3	1.5	0.06	−4.0	4.0	1.5	1	15	15

.

4.1. Analysis of Proposed Runway Detection

An overview of the proposed method for runway detection is shown in Figure 4. The original infrared image is shown in Figure 4a. Based on (2), quantifying the value from [0,1] to [0,255], the edge indicator g is presented in Figure 4b. Three classification segmentation via the proposed three-thresholding method is illustrated in Figure 4c, red parts of which are candidates for the runway. Followed by tracing boundaries of the region in red, an initial contour is outlined in Figure 4d, which will expand to the final runway. After 16 iterations, the curve is moved as shown in Figure 4e. In Figure 4d,e, the regions enclosed by the white rectangles are ROIs for the next iteration. When the level set function meets the stopping criterion, evolution halts and the final runway is highlighted in Figure 4f. The outcome is almost identical with human visual judgment.

Unfortunately, there is no ground truth to assess the error. In order to quantize the performance of our proposed method, we assume that the shape of a runway is quadrilateral and draw the ground truth manually for 100 continuous frames. The initial contour in the true runway is significant. Therefore, we define three indices below, which depict the error of the initial contour, the final contour and the area differece between ground truth and experimental results, respectively.

(1)

The error of the initial contour $e_0$. It is apparent that, if the initial contour lies in the range of the ground truth, it is appropriate for the proposed method. For the error in the initial contour, we employ a “one-vote negation evaluation strategy”, namely, we define the initial error $e_0=0$ if all the points of an initial contour lie within the ground truth; otherwise, $e_0=1$.

(2)

The error of the final contour $e_1$. The distance between the points of the final contour and the ground truth reflects the accuracy of the proposed method [39]. Thus, we define the final error $e_1$ as follows:

$$e_1=\frac{{\sum }_{i=1}^N\left|\mathrm{dist}(P_i,{\mathbf{C}}_{\mathbf{g}})\right|}{N}\mathrm{pel},$$

(31)

where the function $\mathrm{dist}(·)$ returns a signed distance between the point $P_i$ in the final contour and the nearest ground truth edge. Based on the premise that the manually drawn ground truth is sufficiently precise, it is obvious that the smaller the $e_1$ is, the more accurate the runway detection is.

(3)

The area differece $e_2$. Jaccard Similarity measure calculates [40] the similarity in the G ground truth area and the D detected area. In this measurement, the maximum value 1 means completly similar, while 0 means completely dissimilar. Naturally, we defined the third error index, given as below:

$$e_2=1-J=1-\frac{|G\bigcap D|}{|G\bigcup D|}.$$

(32)

Based on the error indices mentioned above, the diagram for error $e_0$, $e_1$ and $e_3$ is shown in Figure 5. The $e_0=0$ for all the tested images, $e_1$ takes values in the interval $[0.73,3.11]$ and $e_2$ in the interval $[0.16,0.51]$. As for the index $e_2$, the acreage of target runway is very small, a slight difference between the ground truth and tested area may lead to an obviously large error. However, if we replace $\frac{|G\bigcap D|}{\left|G\right|}$ with $\frac{|G\bigcap D|}{|G\bigcup D|}$ in Equation (32), the error drops to 0.15 on the average. This means that the runway extracted overlaps the 85% percent of the area ground truth covered. In order to demonstrate the error distinctly, we show a graphical representation of the runways detected and the corresponding ground truth, as shown in

Table 2. Comparison of results with ground truth.

Sequence Number	${\mathit{e}}_{\mathbf{1}}$ (pel)	Proposed Method	Ground Truth
100	3.11
84	2.62
97	2.57
16	2.57
75	2.47

, whose $e_1$ are the largest ones. Note that there is only a slight visual difference in the runways between the proposed method and the ground truth. Moreover, the error $e_1$ is not ideally zero for some cases. However, it is sufficient to improve a pilots’ visual perception and runways detected automatically can significantly improve a pilots’ spatial awareness. The relevant parameters are employed as shown in Table 1. The detection performance verifies that the parameters listed are image independent.

4.2. Comparison with Existing Methods

Compared to line segment based runway detection, Canny edge detector and EDLines have excellent accuracy and resistance toward false detections. In Figure 6, edge images by Canny’s method are listed in the first column. It is clear that this method succeeds in extracting the edge of runways, while it is intricate to determine runways from patches of line segments during post processing. In Figure 6b, the line segments of runway by the EDLines are missing. Although three line segments of runway are detected in Figure 6b, three out of 192 line segments is still a tough decision. On the contrary, the proposed method gives a good solution, as listed in the third column.

Compared to Li’s method [27], our ROI-based level set method has significantly better time performance. The average execution time of Li’s method and ours with the same computer configuration (CPU: Intel(R) Core(TM) i7-4790 (3.60 GHZ), RAM: 16 GB and OS: Windows 10 (64 bit)) is listed in

Table 3. Algorithm execution time (s/f).

Processing Time	Min	Max	Average
Li’s Method [27]	0.172	1.844	0.821
Proposed Method	0.015	0.063	0.035

. Our processing speed is 28f/s, which meets the demand for real-time processing.

In addition, we placed a subset of a runway as the initial curve (see Figure 7b) manually. Our final result in Figure 7c or f is similar to Li’s in Figure 7d after 100 iterations. However, it is difficult for Li’s method to design an appropriate stopping criterion when the curve evolves to the desired runway. It is obvious that the curve in Figure 7e has a bigger size than the exact runway after 200 iterations. In this application, a traditional level set with fixed maximal iterations is inadvisable. Simultaneously, it is worth pointing out that the virtual runway can provide clues about the initial contour and the stopping criterion leading to the accurate runway, which has practical significance in automatic processing.

Figure 8 demonstrates the comparison between Chan’s model [41,42] and our method. Due to lack of edge information, Chan’s model cannot separate the runway from the touch zone, as shown in Figure 8d. In this case, our method prevails in accuracy with the assistance of the edge indicator, as shown in Figure 8b. The final result of our method is shown in Figure 8c, which excludes other irrelevant regions apart from the runway.

As for time performance, it was reported that Balla-Arabe’s method [43] reaches 0.7 s per frame for $481\times 321$ plane images (CPU: AMD Athlon(TM) 5200 2.3 GHZ, RAM: 2 GB, Matlab 2010b, (MathWorks, Natick, MA, U.S.). Shi [29] proposed a real-time algorithm for the approximation of Level-Set-Based curve evolution and can segment the plane image in 0.0249 s (CPU: Intel 1.6 GHZ). Nevertheless, a circle outside the plane was plotted manually as an initial curve in Shi’s method, which is not feasible in batch processing.

5. Conclusions

We proposed an entirely new approach for runway detection from infrared aerial images incorporating synthetic vision. Synthetic vision is capable of generating virtual runways in reality. However, chances are that pilots are not able to confirm the authenticity of this information if errors occurred in either the navigation data or in the 3D scene database. Therefore, pilots are provided with credible visual cues by means of structures such as the runway, river, freeway, overpass and other landmarks explored under real circumstances [18,19,44]. A framework for runway detection on the basis of synthetic vision and an ROI level set method was presented. We first analyzed the structure of a runway and found that the runway adjoins a touch zone in the infrared image. In order to differentiate the two regions, edge information was used. Subsequently, a level set method without reinitialization was employed, which comprised an edge indicator. Then, synthetic vision provided an initial level set function and a stop condition. For the initial level set functions, we presented a three thresholding method to guarantee that the initial level set function was induced by the subset of IR runway. Furthermore, an ROI based evolution strategy was proposed to achieve real-time processing. Compared to traditional runway detection, the virtual runway proposed in this paper detects an approximate location of an IR runway; hence, it can output an exact runway without any follow-up recognition. Other methods need to carry out identification after runway extraction. In this paper, we assumed that the virtual runway can overlap with the IR runway. This assumption demands slightly reliable synthetic vision. If this prerequisite is not met, an incorrect sign of v in Equation (28) may make the level set function evolve in the wrong direction. This is a problem that we will focus on in future research, along with creating more reliable synthetic vision.