UPTEC F 16044
Examensarbete 30 hp Juni 2016
Runway detection in LWIR video
Real time image processing and presentation of sensor data
Erasmus Cedernaes
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Abstract
Runway detection in LWIR video
Erasmus Cedernaes
Runway detection in long wavelength infrared (LWIR) video could potentially increase the number of successful landings by increasing the situational awareness of pilots and verifying a correct approach.
A method for detecting runways in LWIR video was therefore proposed and evaluated for robustness, speed and FPGA acceleration.
The proposed algorithm improves the detection probability by making assumptions of the runway appearance during approach, as well as by using a modified Hough line transform and a symmetric search of peaks in the accumulator that is returned by the Hough line transform.
A video chain was implemented on a Xilinx ZC702 Development card with input and output via HDMI through an expansion card. The video frames were buffered to RAM, and the detection algorithm ran on the CPU, which however did not meet the real-time requirement. Strategies were proposed that would improve the processing speed by either
acceleration in hardware or algorithmic changes.
ISSN: 1401-5757, UPTEC F 16044
Examinator: Tomas Nyberg
Ämnesgranskare: Anders Hast
Handledare: Joakim Lindén
Populärvetenskaplig sammanfattning
För att kunna landa ett flygplan på en landningsbana behöver en pilot antingen visuell kontakt med landningsbanan eller ett styrsystem som sköter landningen automatiskt. I många fall finns inte autopilot tillgänglig för piloten eftersom det kräver dyra investeringar för ägaren av landningsbanan. Då måste piloten själv kunna bekräfta att han kan se landningsbanan i tillräckligt god tid för att kunna avbryta en landning. I disigt väder kan dessa marginaler vara knappa eller obefintliga, och piloten blir tvungen att ha ett externt system att förlita sig på. Därför har detta examensarbete undersökt huruvida värmekameror monterade på flygplanet skulle kunna hjälpa till vid landningar, så att piloten har ett stödsystem även om landningsbanan inte är utrustad med något. Eftersom värmekameror både kan se bättre i mörker och genom dis än vanliga kameror sågs dessa som bra kandidater att undersöka.
Videosekvenser i det långvågiga IR-spektrat (det vill säga värmestrålning) studerades för att hitta detaljer som kunde användas för att kunna urskilja landningsbanor. Land- ningsbanan bildar en form av parallelltrapets i kameraobjektivet under inflygningen; ko- rtsidorna på landningsbanan bildar två parallella, horisontella linjer, medan långsidorna bildar två vertikala linjer, svagt vinklade ut från lodlinjen med ungefär lika men motsatt vinkel. När planet närmar sig landningsbanan blir vinkelskillnaden större mellan långsi- dorna på landningsbanan, och långsidorna ser större ut i bilden. En algoritm skrevs som tog fram linjer i bilderna som stämde in på detta mönster.
För att kunna använda algoritmen i realtid behövdes någon form av hårdvaruaccelerering.
Därför utreddes det om de olika delarna av algoritmen kunde köras på en FPGA (Field Programmable Gate Array). Xilinx Zynq sågs som en lämplig plattform för ändamålet eftersom den kombinerar en CPU med en FPGA i samma förpackning. Flera block togs fram som omvandlade inkommande video, via HDMI-kabel, till lämpligt format så att det kunde sparas till arbetsminnet, där den sedan kunde läsas och manipuleras.
Den slutgiltiga implementationen av algoritmen kördes i CPU-delen på plattformen men uppnåde ej målet att köras i realtid. Flera förslag lades därför fram för hur olika delar av algoritmen skulle kunna överföras till FPGA-delen för att uppnå en ökad prestanda, och även hur algoritmen skulle kunna simplifieras: främst genom att använda information från tidigare bilder.
Studiens slutsats är att IR-kameror kan användas som ett kompletterande hjälpmedel vid
landningar i dålig sikt.
Acknowledgements
I would like to express my sincere gratitude toward all of those who have been by my side during the development of this project. Many have shown great interest in the project and given me interesting ideas to explore, often without much effort from my side.
Special thanks are directed toward my family: Alice, the love of my life; my parents; my two brothers; Mischa, Saatchi and Tootsie; Mormor and Morfar.
I would like to thank the people at Saab: my supervisor Joakim Lindén for the immense dedication to the project, Per Emilsson for pushing me to finish on time, Johan Zandén for providing IR sensor data, all of the proofreaders for helping me improve the quality of the report, and the rest at Avionics Systems who share interesting anecdotes and brighten up each day.
I am grateful to my supervisor at Uppsala University, Anders Hast, for listening to my
ideas and providing feedback to keep me at a good pace.
Abbreviations
AMBA Advanced Microcontroller Bus Architecture
Open standard created by ARM for interconnections between func- tional blocks.
AXI Advanced eXtensible Interface
Part of the AMBA specification meant for high performance data transfers.
CPU Central Processing Unit DMA Direct Memory Access
Access to memory without need of processor engagement.
FIFO First In First Out
A type of memory buffer. The first data that is written will be the first to be read (queue analogy: the first person in line will be first served, and so on).
FPGA Field Programmable Gate Array A reconfigurable integrated circuit.
HDMI High-Definition Multimedia Interface
A standard for protocols, signals, electrical interfaces and mechanical requirement of transfer of uncompressed video.
HT Hough Transform
Algorithm for finding parametrized shapes in images.
IC (Monolithic) Integrated Circuit A set of electronic circuits on a single chip.
IP
block/core
Intellectual Property block/core
A licensed functional block for programmable logic designs.
LWIR Long Wavelength InfraRed
Electromagnetic radiation that is emitted by body-temperature black bodies.
OS Operating System
PL Programmable Logic, such as an FPGA PS Processing System, i.e. the CPU
SoC System on a Chip
An IC that integrates all of the components of a computer, or a similar system, on a single chip.
VHDL VHSIC (Very High Speed Integrated Circuit) Hardware Description Language
A language for describing digital electronics architecture and be-
haviour. Commonly used as a programming language for FPGAs.
Useful words
Approach The final stage in a landing before touching ground.
Bare-metal application
A standalone application, i.e. runs without an OS, and thus has complete and direct access to the hardware.
ChipScope A Xilinx tool that is used for recording and viewing internal FPGA signals. Can be useful for debugging.
False positive (negative)
Erroneously classifying something as true (false) when it actually is false (true).
Real-time This suggests that a system can process data at least as fast as it is received (throughput is equal to the input rate). It can also put limitations on the delay (commonly called lag) between the start of a data transmission and when the data has been processed.
True positive (negative)
Successfully classifying something as true (false).
Xilinx A supplier of programmable logic devices and inventor of the FPGA.
[25]
Contents
1 Introduction 1
1.1 Goal . . . . 1
1.2 Delimitation . . . . 1
1.3 Research questions . . . . 1
1.4 Related research . . . . 1
1.4.1 Pin point landing . . . . 1
1.4.2 Lane assist . . . . 2
1.4.3 Safe landing site detection . . . . 2
1.4.4 Runway incursion detection . . . . 2
2 Theory 4 2.1 Rectification from trapezoidal to square coordinate system . . . . 4
2.2 2D camera projection . . . . 4
2.3 Hough transform (HT) . . . . 4
2.3.1 Hough line transform . . . . 5
2.4 Bresenham’s line algorithm . . . . 5
2.5 Integral image . . . . 6
2.6 Spatial correlation and convolution . . . . 6
2.7 Sobel operator . . . . 7
2.8 Canny edge detection . . . . 7
3 Hypotheses - Algorithms 8 3.1 Edge detection and HT after perspective correction . . . . 8
3.2 Edge detection with a filtered HT . . . . 8
4 Results - Algorithm prototyping 9 4.1 Edge detection and HT after perspective correction . . . . 9
4.1.1 Perspective in-correction . . . . 11
4.1.2 LWIR-runways . . . . 11
4.2 Edge detection with a filtered HT . . . . 11
4.2.1 Peak search . . . . 13
4.2.2 Hough space traversal . . . . 15
4.2.3 Angle filter . . . . 16
4.2.4 Attitude estimation . . . . 18
5 Hardware 20 5.1 Xilinx Zynq-7000 . . . . 20
5.1.1 PS-PL communication . . . . 21
5.1.2 Protocol summary . . . . 21
5.1.3 Operating System (OS) . . . . 22
5.2 Xilinx ZC702 Development Board . . . . 22
5.3 Vivado Design Suite . . . . 23
5.3.1 Vivado IP Integrator (IPI) . . . . 23
5.3.2 Xilinx Software Development Kit (SDK) . . . . 23
5.3.3 Vivado High Level Synthesis (HLS) . . . . 23
6 Result - Implementation of video stream in hardware 24
6.1 FMC Imageon HDMI in/out card . . . . 24
6.1.1 YCbCr . . . . 24
6.1.2 Sync words . . . . 24
6.1.3 Blanking . . . . 25
6.1.4 Video dimensions . . . . 25
6.2 Simple video output generator . . . . 26
6.2.1 Multiple clock sources and metastability . . . . 26
6.3 Video to AXI-Stream . . . . 26
6.3.1 AXI-Stream for video . . . . 27
6.4 FIFO for AXI-Streamed video . . . . 27
6.5 Video Direct Memory Access (VDMA) Stream to Memory Map (S2MM) . 27 6.5.1 First implementation . . . . 28
6.5.2 Second implementation: Increasing the Master S2MM clock frequency 29 6.5.3 Third implementation: Connecting a FIFO for AXI-Streamed Video 29 6.6 Video Direct Memory Access (VDMA) Memory Map to Stream (MM2S) . 30 6.6.1 Video flow control . . . . 30
6.6.2 Video Output Generator . . . . 30
6.7 Video contrast scaling . . . . 30
7 Result - Software implementation of Algorithm on development card 32 7.1 Contrast adjust . . . . 32
7.2 Horizontal Sobel . . . . 33
7.3 Magnitude filter . . . . 33
7.4 Weighted multiplication . . . . 33
7.5 Hough transform . . . . 34
7.5.1 Suggested improvement for speed . . . . 34
7.6 Peak search . . . . 34
7.7 Plot lines . . . . 34
7.8 Profiling . . . . 35
7.8.1 NEON SIMD instructions . . . . 35
7.8.2 Dual core processing . . . . 36
8 Further development - Hardware implementation of algorithm on development card 37 8.1 Contrast adjust . . . . 37
8.2 Horizontal Sobel . . . . 38
8.3 Magnitude filter . . . . 39
8.4 Weighted multiplication . . . . 39
8.5 Extract nonzero . . . . 39
8.6 Possibility of HT in FPGA . . . . 40
8.6.1 HT with buffered nonzero pixels . . . . 40
9 Discussion 41 9.1 Real-time performance . . . . 41
9.2 LWIR for runway detection . . . . 41
9.3 Approach behaviour . . . . 41
10 Conclusions 43
References 44
1 Introduction
Avionics Systems is a business within the global aerospace and security company Saab, which develop and manufacture qualified electronics, software and mechanics for aircrafts, helicopters and other demanding applications. Most of their projects are conducted in an international environment and they deliver, among other, products to two of the world’s largest aircraft manufacturers; Airbus and Boeing.
In the world of avionics, the requisites are strict regarding safety and robustness of the technical systems. Using a variety of electro-optical sensors, such as visual cameras, infrared cameras and radar can help the pilot when visibility conditions are not optimal.
A problem that arises is how to best extract and present the information from one or more such sensors.
1.1 Goal
The goal is to develop a method for detecting runways in LWIR-video. The suggested method shall then be evaluated for hardware acceleration and robustness, such that it can run in real-time during approach (the final stage of landing before touching ground). The final presentation shall show the runway location and possibly a positional estimation of the aircraft relative to the runway.
1.2 Delimitation
This report has limited its application to standard runway approaches: the runway is in front of the aircraft; the aircraft is almost levelled with the horizon; and the approach angle is about 3
◦(normal approach angle).
1.3 Research questions
This report will reach its goal by also trying to answer the following questions:
• Is a combination of FPGA and CPU advantageous for image processing?
• Are LWIR sensors suitable for detecting runways?
1.4 Related research
Focus was put on research topics that bring up important sides of the posed problem:
accuracy, security and ability to accelerate in hardware. Some of them are described in the following sections.
1.4.1 Pin point landing
Kaufman et. al. [10] developed a positioning method for pin point landings on stellar objects without atmosphere, such as the moon. The method identifies shadows on the stel- lar object, formed by craters or mountains, and then matches them with a pre-calculated multi-resolution binary shadow map. This resulted in a high precision estimation of the lunar lander’s position as long as the pre-calculated images did not differ too much in size from the actual live footage.
Mourikis et. al. [13] developed VISINAV, an algorithm for precision planetary landing. It
integrates absolute positional measurements via matching of Mapped Landmarks (land-
marks that are extracted together with their global position before entry) as well as a
relative positional measurement with Opportunistic Features (significant points in live
footage that can be tracked between consecutive images). The algorithm was tested with
a sounding rocket test and attained a positional touchdown error of 6.4 m. The VISINAV
algorithm has been partially implemented in an FPGA by Morfopolous et. al. [12], but did however not achieve real-time execution.
Both these algorithms can achieve great absolute accuracy at the cost of high computa- tional complexity and the requirement of image data before entry. This is very common in space exploration since excursions can be precisely planned, both in space and time.
However, these algorithms cannot handle the case when the exact position and appearance of the landing area is unknown before entry/approach.
1.4.2 Lane assist
Jung et. al. [9] developed an efficient lane departure detection algorithm for low com- puting power systems. It uses the integral image to calculate directional features in the images, similar to Haar-features [21]. Furthermore, the line candidates can be verified if the mounting position of the camera is known.
This problem is very similar to finding a runway in images as both a car lane and an airport runway have line markings with similar characteristics. The main difference be- tween posed problems is that the relative position and apparent direction of the car lane markings are very static as the car is always on the road, at the same height and trav- elling almost always in the same direction relative to the road. Therefore, some of the generalisations that make this algorithm as fast and accurate cannot be made for runway markings, as their relative position and direction can vary to a much higher degree.
1.4.3 Safe landing site detection
Shen et al. [16] developed an algorithm for detecting safe landing sites through ocular assessment. The algorithm is separated into four blocks of image processing that are applied to an image that covers a large field of view. The first block detects the horizon.
This is done by first applying the Canny edge detector and then Hough transforming the results. A selection ofthe largest peaks from the Hough transform are considered and each is valued by its distance to strong edges in its neighbourhood. The final horizon candidate is then adjusted by nudging it into the direction of the strongest edges in its neighbourhood. The second block assesses the roughness of the ground (the area below the horizon) by first separating it into square bins of equal size (whose size are adjusted for the height of the aircraft). Each bin then is assigned a cumulative hazard strength, which is the number of pixels within a bin whose intensity exceeds a prespecified safeness threshold. The third block separates the roughness measure by the means of a K-means clustering method and groups areas by a region-growing method. The last block then measures the dimensions of each region by transforming from the image space into the real space, using the attitude and altitude of the aircraft.
This method is clearly very computationally demanding as it consists of many complex steps. The importance of detecting the horizon is emphasized as it narrows down the search region for landing sites.
1.4.4 Runway incursion detection
Vegula et. al. [20] developed an algorithm for detecting runway incursions. This is done
by first detecting the horizon via a grey level threshold per column in the image. The
runway search area is then limited to a triangular area with its vantage point horizontally
centred in the image, vertically placed at the horizon, with the other two corners placed
in the lower corners of the image. Then, by applying the Canny edge detector, edges are
extracted, and sent onward to a Hough line transformation, in which the two strongest
peaks correspond to the runway edges. Another method is then proposed for tracking the runway when detected. This is done by means of cropping out a part of the image around the detected runway and using it for template matching in the following frames.
This method suffers from the limitation of the runway location in the image. As the search
area is very narrow, it thus allows very little space for manoeuvring, such as rolling or
yawing, which means the runway can potentially go outside of the search area, even
though it is still visible in the image. Also, the assumption that the two strongest peaks
in the Hough transform correspond to the runway is a simplification that will work for
some runways, which however is not always true if there are other long lines in the image,
such as approach lights in front of the runway.
2 Theory
The following sections provide a background with theory of image processing algorithms and operations that are used in the proposed algorithms. Understanding this section is important to be able to understand the details of this report, but is not necessary to able to understand the general context.
Other theory was studied as well, but for various reasons did not end up being used.
The most thoroughly studied example is the Retinex algorithm, which is used for image improvement in low contrast imagery [14]. It did not show up in the final implementation since none of the example footage was taken during inclement weather. Other examples are feature detection and extraction algorithms, which were not used since they heavily rely on beforehand taken footage of the scene that can be related to in position and attitude.
2.1 Rectification from trapezoidal to square coordinate system
Consider a normalized coordinate system (i.e. 0 ≤ x ≤ 1, 0 ≤ y ≤ 1) with corners in p
i, i = 1, 2, 3, 4 that should be linearly transformed such that the corners map onto the points p
0i. A point p = (x, y) is then translated to the point
p
0(x, y) = p
1+ x(p
2− p
1) + y(p
3− p
1) + xy(p
1+ p
4− p
2− p
3). (1) This can easily be verified by inserting the corners p
i∈ {p
1, p
2, p
3, p
4} = {(0, 0), (1, 0), (0, 1), (1, 1)}, which yields p
0i.
2.2 2D camera projection
Consider if we want to create a 2D projection of a scene using a pinhole model. This can be achieved by three simple transformations [1, pp. 190-193, 256-259]:
1. Translate the coordinate system such that the camera is placed in the origin. This is equivalent to subtracting each point in the coordinate system with the coordinates of the camera.
2. Rotate the coordinate system such the camera becomes aligned with a unit vector.
This is equivalent to multiplying each coordinate with the inverse rotational matrix of the camera. A rotational matrix is orthogonal by construction and therefore equal to its transpose, that is M
−1= M
Tfor an orthogonal matrix M .
3. Project the world coordinates to the sensor by multiplying by the ratio F/z, where F is the focal length (a parameter which is determined by the camera in use).
2.3 Hough transform (HT)
A parameterized shape can be found by applying the HT to an image. An infinite amount of shapes can go through each pixel of an image, however, each of them must satisfy the equation of the shape. The computational attractiveness of the HT comes from the fact that it uses a discretised accumulator with one dimension for each degree of freedom for the description of the shape. Each point in the accumulator then corresponds to a shape with specific parameter values, discretised to a chosen resolution. The value of a point in the accumulator is the number of pixels which fulfills that specific parametrization, which makes it easy to sort out which parametrization is most prominent in a figure.
The computational complexity of the HT is linear with respect to the number of pixels
considered since each pixel generates a fixed amount of parametrization candidates.
2.3.1 Hough line transform
Lines can be found in images by using the line version of the HT. The parametrization then specifically becomes
ρ = x cos θ + y sin θ, −D
−≤ ρ ≤ D
+, −π/2 ≤ π ≤ π/2 (2) where ρ is the distance from the origin, x and y, the position of a pixel in the image, θ the angle of the line and D
−, D
+the bounds for the perpendicular distance from the origin to any other pixel in the image. [15, pp. 330-336] [6, pp. 775-760]
Using this parameterization, each point in the accumulator corresponds to a line in dis- play space (Figure 1a), and each point in the image corresponds to a sine curve in the accumulator (Figure 1b).
θ y ρ
x
(a)
ρ
θ
(b)
Figure 1: HT point correspondence (a) A point in the accumulator corresponds to a straight line in the image (b) A point in the image corresponds to a sine curve in the accumulator
The Hough line transform is robust against white noise since there needs to be structure to generate many votes for a parameter pair. Furthermore, a single picture can gener- ate many (often millions of) votes, which means that even “short” lines can get many votes.
When line candidates have been chosen, the pixels that correspond to each line segment can be found using Bresenham’s line algoritm, described in the following section.
2.4 Bresenham’s line algorithm
Consider a line that goes between two points p
1= (x
1, y
1) and p
2= (x
2, y
2) with dif- ferences dx = x
2− x
1and dy = y
2− y
1. If the differences are positive and the slope is limited (dx ≥ 0, dy ≥ 0 and dy ≤ dx), then an eight-connected line can be drawn between p
1and p
2with a simple formulation of Bresenham’s line algoritm. The integer version is described in Listing 1.
Listing 1: Pseudo code for Bresenham’s line algorithm
P l o t B r e s e n h a m L i n e( p1 , p2 ) { dx ← p2 . x - p1 . x
dy ← p2 . y - p1 . y D ← 0
x ← p1 . x y ← p1 . y
w h i l e ( x ≤ p2 . x ) { if (2* D ≥ dx ) {
D ← D + dy - dx y ← y + 1
} e l s e { D ← D + dy }
p l o t P i x e l( x , y ) x ← x + 1 }
}
The assumption dy ≤ dx implies that each column needs only to have one pixel drawn in order to guarantee a connected line segment. Extending the algorithm to slopes outside of that interval can be solved by switching places between the x- and/or y-coordinates of the start and/or end-points, which then mirrors the coordinates into the valid octant. [1, pp. 355-357] [22]
2.5 Integral image
An integral image is defined as the cumulative sum of the rows and columns of an image.
It can be described by the following recursion formula I
0,0int= I
0,0I
(x+1,y)int= I
(x,y)int+ I
(x+1,y)I
(x,y+1)int= I
(x,y)int+ I
(x,y+1),
(3)
where I
intis the integral image and I is the original image. [21]
2.6 Spatial correlation and convolution
Gonzalez and Woods clearly differentiate between the two terms correlation and con- volution in [6, pp. 168-172], which can otherwise be a source of confusion when used imprecisely: Assume that a filter κ is of size n × m and, n = 2a + 1, m = 2b + 1, with a and b as positive integers, and the image I has been padded appropriately (with zeroes for example). The spatial correlation between κ and I is then defined as
κ(x, y) I(x, y) =
a
X
s=−a b
X
t=−b
κ(s, t)I(x + s, y + t), (4)
while the convolution is defined as κ(x, y) ? I(x, y) =
a
X
s=−a b
X
t=−b
κ(s, t)I(x − s, y − t). (5)
Thus, the spatial correlation and convolution are related, such that, flipping either the
filter or the image along each dimension turns a correlation into a convolution and vice
versa.
2.7 Sobel operator
Edges can be enhanced in images by spatially correlating them with a filter kernel that approximates the derivative. One of the most common edge enhancers is the Sobel op- erator, which approximates the gradient in the image by first calculating approximations to both the horizontal and vertical derivatives with the kernels:
κ
x=
−1 0 1
−2 0 2
−1 0 1
, κ
y=
−1 −2 −1
0 0 0
1 2 1
. (6)
The magnitude and direction of the edges can finally be calculated with
|G| =
q(κ
xI)
2+ (κ
yI)
2and arg(G) = atan2(κ
yI, κ
xI), (7) where G can be seen as the gradient of image I (note the use of the correlation oper- ator rather than the convolution). One can look upon the two Sobel kernels as being combinations of an averaging operator in one direction and a central difference operator perpendicular to that direction since the two Sobel kernels can be written as
κ
x=
1 2 1
h
−1 0 1
i, κ
y=
−1 0 1
h
1 2 1
i. (8)
The averaging operator smooths out the result to reduce jitter perpendicular to the di- rection of the derivative. [6, pp. 728-731]
2.8 Canny edge detection
The Canny edge detector algorithm is a popular and efficient algorithm for finding edges in images. It is based on three objectives: Low error rate, good localization and single edge point response. The algorithm is carried out as follows [15, pp. 326-329] [6, pp.
741-747]:
1. Smooth the image I with a Gaussian filter kernel.
2. Apply a simple edge enhancer (such as Sobel, Prewitt, Roberts).
3. Use the horizontal and vertical edge enhancement results to calculate magnitude and direction of edges with Equation (7).
4. Suppress non-maxima to remove edges that are weaker than ones in its neighbour- hood. Set magnitude to zero for pixels that have a neighbour with similar direction (within a tolerance of 8-connectedness) but larger magnitude.
5. Detect edges: Start an edge-follow for pixels whose magnitude is larger than an
upper threshold T
highand label each traversed pixel as an edge. Keep on following
the edge via eight-connected neighbours until the magnitude reaches below the lower
threshold T
low.
3 Hypotheses - Algorithms
Two different methods were tested for finding runway candidates in images. Those two hypotheses are presented in this section along with motivations as why either one would work.
3.1 Edge detection and HT after perspective correction
This algorithm needs an attitude (pitch, roll and yaw) and position approximation as input to calculate the position of the runway in the image using its 2D projection (Section 2.2). Then an approximate region of interest is cropped out from the image and rectified into square coordinates (section 2.1), thus simulating a top-down view. This would turn the runway into a rectangle irrespective of view angle. To be able to find the rectangle corresponding to the runway, two further image processing steps are taken. First, the rectified image is filtered to detect edges. Then, those edges are Hough transformed, to be able to find long lines. The final step is to search the HT accumulator for potential runway candidates. Since the runway has been rectified, parallel lines will be located in the same θ-column in the HT accumulator, which results in a simple search of two simultaneous peaks for each θ-column.
This method is mainly chosen as a candidate since it significantly reduces the number of runway candidates given by the HT by reducing the peak combinations to search for.
Also, the perspective correction can produce a rectified image with a resolution that can be decoupled from the input image resolution by resampling at specific resolutions. This fact can be of significant importance when deciding what to accelerate in hardware, as programmable logic cannot handle dynamic memory allocations. Further, by rectifying the runway, the perspective distortion is removed and the texture of the runway (such as markings and lights) can then be used in the verification process, which has been proven to be very effective by Jackson et. al. [8].
3.2 Edge detection with a filtered HT
The image is horizontally Sobel filtered to find vertical lines. A reasonable amount of pixels can then be passed on for HT by using adaptive thresholding. By keeping the sign (positive/negative) of the filtered image, lines found by the HT can be matched with similar, but opposite signed line segments. Thus, the left and right side of the runway can more easily be combined.
The runway’s movement in the images during landing can be utilised to simplify the problem. When the runway is far away, it will be located in the upper section of the image (near the horizon) and the sides of the runway will seem almost vertical, and parallel, in the image. When the airplane gets closer to the runway, the sides of the runway will become elongated (stretched further down in the image) and at the same time more angled; forming a vantage point in the upper part of the image. Therefore, points on the runway edges that are close to the top of the image can contribute to both vertical lines and angled lines, however, points close to the bottom cannot contribute to vertical lines since the runway will only be visible in those parts of the images when viewed from closeby.
This method was chosen as a candidate as it consists of several computationally sim-
ple operations (except for the HT), which are almost all possible to, accelerate in an
FPGA.
4 Results - Algorithm prototyping
The following sections describe the process of evaluating the algorithm candidates. The potential upsides and downsides are presented in the end of the corresponding section for each algorithm candidate.
4.1 Edge detection and HT after perspective correction
The prototype algorithm described in this section corresponds to the hypothesis described in Section 3.1.
Figure 2 illustrates two cut-outs from an image taken from Google maps. The left figure shows a rectangular cut-out and the right the corresponding perspective corrected cut- out (taken from the dashed red line in the left one). The right image was perspective corrected by mapping each pixel to its corresponding position in the left image (Section 2.1). Pixels that are not mapped exactly onto other are bilinearly interpolated between the closest neighbours.
One can clearly see that the perspective corrected cut-out turns all lines that head toward the vanishing point into parallel lines. However, image artefacts arise during perspective correction, especially for details further from the camera (i.e. closer to the top of Figure 2). This is due to the bilinear upscaling being sampled at a nonlinear grid.
(a) (b)
Figure 2: Cutouts (a) Rectangular cutout (b) Its perspective correction (from the red trapezoidal area)
Figure 3 shows the Canny edge detection of corresponding subfigures in Figure 2. The
black pixels are considered as edges.
(a) (b)
Figure 3: Canny edge detection results (a) Canny edge detection of Figure 2a (b) Canny edge detection of Figure 2b
Figure 4 shows the accumulators that are returned when applying the HT to Figure 3, with black corresponding to a high value and white a low. The five largest peaks in each HT accumulator are marked with squares and the true peaks (ground truth) corresponding to the left and right side of the runway are marked with black crosses. As one can see, the lines corresponding to the runway sides are not among the five largest before perspective correction, but are among the five largest after perspective correction. Also, as expected, the highest peaks are located at very similar angles after perspective correction since many lines are directed toward the vanishing point.
−50 0 50
−200 0 200 400 600 800
θ [
◦]
ρ [pixels]
(a)
−50 0 50
−200
−100 0 100 200 300
θ [
◦]
ρ [pixels]
(b)
Figure 4: HT results with five highest peaks marked with red squares (a) HT of Figure 3a (b) HT of Figure 3b
Figure 5 shows the detected lines that correspond to the five highest peaks in Figure
4 respectively. Perspective correction makes it easier for the HT to find lines that are
parallel and better separated in image space, and therefore narrows down the search area
further than a simple rectangular cut-out. This makes it very easy to find rectangular
portions of the image, since right angles on the ground also become perpendicular in the
rectified image. Also, it is easy to distinguish the markings on the runway, which in turn
can be used to verify where the runway is located as Jackson et. al. have shown [8].
(a) (b)
Figure 5: Resulting lines (a) Lines extracted from Figure 4a (b) Lines extracted from Figure 4b
4.1.1 Perspective in-correction
A significant problem that might appear is that errors in the attitude estimation propagate to produce an incorrect perspective correction. This is because an incorrect estimation of the camera angle leads to an incorrect estimate of the position of the vanishing point.
This results in straight lines becoming curved after perspective correction, more so closer to the vanishing point. This effect can be reduced by cropping away a part of the top of the perspective corrected image, since the bottom part of the image is less affected.
This does, however, not eliminate the problem, and introduces other issues, such as lines becoming shorter and thus giving less votes in Hough space.
4.1.2 LWIR-runways
Real LWIR-footage was acquired at a later stage of this study, which also introduced another, before unknown, deficit with this method: In the LWIR-spectrum, the runway markings become very faint, and almost impossible to distinguish except when viewed from a very short distance. Therefore, the runway markings are not suitable for confirming a runway candidate when using an LWIR sensor.
As a consequence, this method was abandoned in favour of another method, described in the next section, which has no requirement of an attitude estimation, but can improve if one is added.
4.2 Edge detection with a filtered HT
The prototype algorithm described in this section corresponds to the hypothesis described in Section 3.2.
The input IR-image is first contrast enhanced to make sure that small objects in the image with very high/low values become levelled with the rest of the image. This must be done to make sure that the approach lights on and beside the runway does not blow out the image, thus underexposing the rest of the image. Matlab has a built-in function called imadjust, which does this.
The contrast enhanced image is then horizontally Sobel filtered, which emphasizes the
vertical lines/edges in the image. This is advantageous compared to a full Sobel filter
since the vertical lines of the runway are larger than the horizontal lines in the image
projection. If we assume that the lines in the rest of the image are uniformly rotated,
then, by only filtering out the vertical lines we would provide a better “signal to noise ratio” since a larger portion of the vertical lines correspond the the runway than the portion of horizontal lines. Also, homogeneously coloured areas become bounded by a positive and a negative flank of the horizontally Sobel filtered image, which makes it easier to determine if two lines correspond to eachother.
A 3 × 3 kernel is used for the Sobel filter, i.e. κ
x=
−1 0 1
−2 0 2
−1 0 1
, since it reduces the
vertical noise compared to the 1 × 3 kernel κ
x=
h−1 0 1
i, however, at the cost of a larger kernel. Implementing this filtering in the PL would involve two line buffers, with the added delay of at least one row and one pixel.
Figure 6 shows a sample image from an IR-camera with the runway in the upper middle of the image. 7 shows the result when horizontally Sobel filtering that image.
Figure 6: Input image Figure 7: Horizontal Sobel filtering
Only a small percentage of the strongest edges of the Sobel filtered image are then kept by setting a threshold for the magnitude. This has two advantages: First, the amount of stray pixels is reduced, which decreases the amount of noise. Second, the number of pixels processed by the HT are reduced, which shortens the computational time. By tweaking the threshold, a balance can be achieved between false positives and positional accuracy of the runway. For example, Figure 8 shows the pixels that correspond to the 2% strongest edges of Figure 7.
Using the threshold alone is not sufficient for finding the sides of the runway since they are short compared to other lines in the image. Therefore, some extra data needs to be added to the pixels before we can find the runway. One solution is to weight each pixel in the HT according to some principle that would be advantageous for the runway lines.
In our case, the runway is generally very bright and has large derivatives at the vertical edges. Thus, we can weight the luminance in the IR-spectrum with the Sobel image and use this as an input to the HT.
Figure 9 shows the weighted input to the HT, which, in fact, is purely a pixel wise
multiplication between Figure 6, Figure 7 and Figure 8.
Figure 8: Sobel magnitude threshold Figure 9: Weighted input to HT
Figure 10a shows the Hough line transformation of Figure 9, with the most negative (trough), and most positive peak (ridge) marked. The corresponding lines are shown in Figure 10b, with matching colour for each peak. As can be seen, the approach lights in front of the runway give a large contribution in the HT accumulator and thus contribute to the largest peaks. A peak search method is proposed in the next section to make it easier to distinguish the runway from other long lines during approach.
−60 −40 −20 0 20 40 60
100
200
300
400
500
θ [
◦]
ρ [pixels]
Ridge Trough
(a)
Ridge Trough
(b)
Figure 10: Original (a) HT with peaks marked (b) Corresponding lines
4.2.1 Peak search
It is not necessarily the strongest peaks in the HT accumulator that correspond to the runway. Therefore, some clever tricks are needed to increase the chances of finding the runway.
We can assume that the runway gives peaks in the Sobel filtered image of opposite signs on opposing sides if the runway has a similar colour on the left/right side in the IR-spectrum and is surrounded by similar terrain on both sides. In our case, the runway was bright in the IR-spectrum with, relatively, dark terrain surrounding it.
When the aircraft is levelled with the horizon, the projection of the runway will become
symmetric along the vertical axis, with the left side of the runway at equal, but opposite,
angle from the vertical axis as the right side of the runway. Thus, to be able to find the
runway, we should search for two concurrent peaks in the HT accumulator of opposite signs, which are located at equal angle from the vertical axis. The pseudo code for the search method used is shown in Listing 2.
Listing 2: Symmetric peak search
S y m m e t r i c P e a k S e a r c h( H , θ _ s y m m e t r y , m a x A n g l e ) {
% M a k e v e c t o r s w i t h the m a x i m u m / m i n i m u m v a l u e in H for e a c h a n g l e Hp ← m a x i m u m V a l u e F o r E a c h A n g l e ( H );
Hn ← m i n i m u m V a l u e F o r E a c h A n g l e ( H );
% I n i t i a l i s e p e a k . a n g l e ← 0;
p e a k . v a l u e ← 0;
% L o o p for e a c h a n g l e a r o u n d the s y m m e t r y for ( dθ in 0 to m a x A n g l e ) {
c u r r e n t V a l u e ← Hp (θ _ s y m m e t r y + dθ) * Hn (θ _ s y m m e t r y - dθ );
% K e e p l a r g e s t n e g a t i v e p e a k if ( c u r r e n t V a l u e < p e a k V a l u e ) {
p e a k . a n g l e ← dθ;
p e a k . v a l u e ← c u r r e n t V a l u e ; }
}
r e t u r n p e a k ; }
The algorithm searches for a maximum and minimum pair on both sides of the symmetry angle (which is the angle between the upward direction in the image and the true vertical direction). The max/min pair is rated by their product, and the angle corresponding to the most negative product is returned. The product is used since it is only large when both the positive and negative angles have peaks but of opposite signs, which prevents a very large line from being returned if there is no corresponding mirrored line. The algorithm can also be extended by allow an angle tolerance for the symmetry angle, and do several searches slightly shifted to the input symmetry angle, θ_symmetry.
Figure 11 shows the resulting HT and lines which are extracted when using a symmetric
peak search. Appearently, the true combination of (ρ, θ) is not strong enough to be found
in the HT accumulator. A method is proposed in the following section for giving the
runway a larger impact in the HT accumulator.
−60 −40 −20 0 20 40 60 100
200
300
400
500
θ [
◦]
ρ [pixels]
Ridge Trough
(a)
Ridge Trough
(b)
Figure 11: Regular HT with symmetric peak search (a) HT with peaks marked (b) Cor- responding lines
4.2.2 Hough space traversal
We can use the fact that the lines travel downward in the image when the runway gets closer, as well as becoming more angled from the vertical line. Thus when performing the HT, we can exclude lines of specific angles for different parts of the image. For example, vertical lines are expected in the upper part of the image, since the runway sides will have an almost vertical image projection when looked upon from far away. However, vertical lines are not expected further down in the image since we only expect the runway to be visible there when it is very close. The edges were therefore tracked to provide training data for exclusion of angles versus position.
The runway edges were tracked by hand for the two provided image sequences. The angle of the edges could then be plotted versus the vertical position over the course of time.
These probability densities, with normalised sum per angle are shown in Figure 12a and
Figure 12b, which represent the results from image sequence 1 and 2 respectively. It is
clear that only an interval of y-positions correspond to each angle, almost symmetrically
around θ = 0
◦. We can therefore narrow down the search range further by only considering
viable combinations of y and θ, thus reducing noise as well as processing load.
−60 −40 −20 0 20 40 60 0
100 200 300 400
θ [
◦]
y [pixels]
(a)
−60 −40 −20 0 20 40 60
0
50
100
θ [
◦]
y [pixels]
(b)
Figure 12: Probability density for the image sequences (a) Image sequence 1 (b) Image sequence 2
Consider an aircraft flying levelled, straight toward a rectangular runway. In this case, the image projection of the runway will be completely symmetrical along the vertical. If the plane rolls to any side, the symmetry angle of the runway will shift an equal, but mirrored angle in the image projection. But, by knowing the vertical direction when flying, we can correct for this shift in real time and keep the symmetry.
Figure 13a and Figure 13b show the results when correcting Figure 12a and Figure 12b for symmetry, thus reducing the effect roll has on the image projection. This is done by, for each frame, subtracting the mean angle of the runway sides so that the angles of the runway sides are symmetric around θ = 0
◦at all times.
−60 −40 −20 0 20 40 60
0 100 200 300 400
θ [
◦]
y [pixels]
(a)
−60 −40 −20 0 20 40 60
0
50
100
θ [
◦]
y [pixels]
(b)
Figure 13: Probability density for the image sequences with angle corrected for symmetry (a) Image sequence 1 (b) Image sequence 2
4.2.3 Angle filter
These probability densities now show a predictable pattern. With a very rough approxi-
mation, we can mask out these possible areas with partially linear polygons using simple
mathematical expressions and thus reduce the amount of data to process as well as nar- row down the possible candidates in Hough space. For example, setting a bound on the vertical position in the image, y, versus the angle relative to the symmetry line, θ,
y
min(θ) ≤ y(θ) ≤ y
max(θ), (9)
where
y
min(θ) = A · |θ| + B, y
max(θ) = C · |θ| + D, (10) (|| denotes absolute value) gives us a simple rule for which angles to include in our compu- tation. For image sequence 1, we can heuristically, set the parameters to A = 1.5, B = 10, C = 3.0 and D = 100. This gives a 87.6% coverage of all the probable values and 65.4%
less pixels to process. We also get an exclusion of 83.9% of the improbable values but only exclude 20.0% of the probable values. This means that we improve the signal to noise ratio drastically (considering pixels not corresponding to a runway as noise). Figure 14a and Figure 14b show y
min(red line) and y
max(orange line) drawn into the probability density maps. The upper and lower bounds have been scaled for image sequence 2 (Figure 14b) such that the proportions between the image size and bounds are the same (since the image sequences provided are of different sizes).
−60 −40 −20 0 20 40 60
0 100 200 300 400
θ [
◦]
y [pixels]
(a)
−60 −40 −20 0 20 40 60
0
50
100
θ [
◦]
y [pixels]
(b)
Figure 14: Probability density for the image sequences with angle corrected for symmetry.
The upper and lower bounds for the filtering are displayed as a red and yellow line respectively (a) Image sequence 1 (b) Image sequence 2
The number of pixels to process can be reduced further, and the coverage can be improved by searching for “better” combinations of parameters (this can quite easily be done with an optimization routine by formulating an objective function that defines what is “better”).
Be that as it may, even if these measures improve, they might only improve with regards to the data set in use (so called overfitting). Thus, if the runway will traverse the image projection differently, then using the parameters that were optimal for this data set might give worse signal to noise ratio. Therefore, it can be well motivated to use heuristically picked parameters when not enough training data is available; as in this case. We choose to leave the parameters at these values so we do not exclude too much data before the decision taking stage.
Figure 15 shows the results from the proposed HT together with the symmetric peak
search. Figure 15a shows the HT accumulator and the strongest positive and negative peak. Figure 15b shows the lines which correspond to those peaks.
−60 −40 −20 0 20 40 60
100
200
300
400
500
θ [
◦]
ρ [pixels]
Ridge Trough
(a)
Ridge Trough
(b)
Figure 15: Filtered HT with symmetric peak search (a) HT accumulator with peaks marked (b) Corresponding lines
4.2.4 Attitude estimation
The aircraft is, however, seldom levelled. Since landing is a manoeuvre that need to be very precise, the aircraft will most likely need to make several adjustments of yaw, pitch and roll before touching ground. Therefore, the upward direction in the image will often differ from the true vertical direction.
To improve the results, one should try to estimate the attitude of the camera. The instantaneous rotational velocity can be measured with a gyroscope so that one can estimate the rotational velocity of the aircraft. However, one must have an estimation of the true downward direction to get an absolute angle, as well as to counteract any noise from the gyroscope. An accelerometer can do this by approximating the downward direction with the direction of the gravitational force. However, the measured acceleration is affected by the aircraft’s manoeuvres, which distort the true downward direction. To minimize the effects of noise, a sophisticated observer such as the Kalman filter should be implemented to utilise the combination of gyroscope and accelerometer signals. This can typically be implemented in a small separate system/sensor to relieve the image processing system from this task, which then can connect via an interface to read the estimation when needed.
An estimation of the attitude can also be extracted with image processing, possibly with a combination of different approaches, such as: SLAM (Simultaneous Localization And Mapping) which would be able to map the ground, such that the attitude and position can be calculated; Horizon detection, which can approximate the roll and pitch of the aircraft in relation to the horizon; and optical flow, which can approximate the rotational velocity in all three directions, however, without reference point.
Despite the many approaches with image processing, none of them can reliably approx-
imate the attitude and position of the aircraft unless the image is adequately rich in
detail and within visual range of the ground (which certainly cannot be guaranteed in
all weather conditions). Furthermore, image processing is computationally intense, which
would result in a much higher power consumption and much dedicated processing power
compared to a sensor combination of gyroscopes and accelerometers.
5 Hardware
A large part of the project was to implement the algorithm on a platform, on which it should run in real-time. Thus, a platform had to be chosen as well as the architecture for it. The following sections describe the chosen hardware as well as things which were needed to be done before starting the development.
5.1 Xilinx Zynq-7000
The Xilinx Zynq-7000 is a monolithic SoC with a combined dual core ARM Cortex A9 Processor and an Artix/Kintex 7-series FPGA [26]. It was chosen as the main core of the platform since the combination of an FPGA and CPU opens up new possibilities regarding algorithmic construction and processing power. The PL can potentially alleviate the PS by accelerating heavy tasks in “soft” (i.e. programmable) hardware blocks. The alleged benefits of replacing a separated FPGA and CPU with a Zynq SoC are many, according to Xilinx [24] [27]. The most important ones are:
• Efficiency as SoC
Since the PS and PL are manufactured on a single chip, data paths become short between them. This reduces the capacitance of the signals paths, which lowers the energy needed to toggle signals and also increases the maximum operational frequency. All in all, this gives a potential of a high performance per watt.
• Efficiency through connectivity
The PS and PL are highly integrated with multiple pathways to transfer data in between them. Extra pathways also allow the PL to have direct access to the on chip memory of the PS, as well as the possibility to employ the 64-bit ACP (Access Coherency Port) port, which is connected directly to the PS’s snoop control unit to maintain cache coherency in memory.
• Security
Since the PS and PL are placed on the same chip, signals between the two are hidden from the outside and cannot be probed in a simple fashion. Additionally, the Zynq offers user authentication through RSA, encryption with AES-256 and automatic
“zeroing” on tamper detect to keep its information safe after configuration.
Figure 16 shows a simple hierarchical sketch of the Zynq architecture (together with the
main memory), how the CPU, FPGA and memory are connected together with an AXI
interconnect. The three marked data paths GP, ACP and HP are explained in Section
5.1.1.
RAM etc.
Peripherals
CPU
AXI Interconnect
FPGA
Zynq SoC
GP
HP ACP
Figure 16: Important features of the Zynq architecture
5.1.1 PS-PL communication
There are several paths to send data back and forth between the PS and PL in the Zynq. The most common paths are through the General Purpose ports (GP), the Access Coherency Port (ACP) and the High Performance Port (HP). The main properties of them are (source: Table 22-8 [28]):
• General purpose Port (GP)
Used for medium throughput transfers where the CPU is master. It requires CPU attention and is therefore suitable for lighter tasks, such as control functions or peripheral access.
• Access Coherency Port (ACP)
Used for high throughput transfers where the PL is master. It has a very low latency thanks to the optional cache coherency. Writing over the cache can however cause cache thrashing. This port is suitable for transfers of smaller datasets which require high performance.
• High performance Port (HP)
Used for high throughput transfers where the PL is master. HP ports offers the, in total, highest throughput since there are several of them. This port is suitable for transfers of large data sets which need high throughput, such as video streams.
5.1.2 Protocol summary
Xilinx has adapted the AXI protocol as a standard for their IP cores. AXI is an open standard developed by ARM. Here follows shows a summary of the AXI protocols used by Xilinx as well as some of their properties. [29]
• AXI-Lite
A memory mapped interface which has a low throughput with high latency. Used for reading and writing control registers.
• AXI-Full
A memory mapped interface which has a high throughput. Used for reading and
writing larger chunks of data.
• AXI-Stream
A streaming interface which has a high throughput with low latency. Used for unidirectional data transfers, such as real-time processing.
5.1.3 Operating System (OS)
There are several alternatives regarding the general setup of the processing system when writing software for the Zynq. Here follows a summary of the advantages/disadvantages of the main OS alternatives (Bare-metal is synonymous to running the program without an OS): [30]
• Bare-metal
Running without an OS is light-weight. There is only one process running on the CPU, hence, all resources will be uniquely dedicated to that process. This results in very reliable performance, with low variation in execution time and low latency.
• Linux
An OS adds some latency and strain on the CPU since it manages the whole processing system: memory (dynamic memory allocation); which processes to run (threads); and peripherals (drivers). This makes it easier to develop applications since many tasks can be solved at a higher level of abstraction.
5.2 Xilinx ZC702 Development Board
The hardware and software implementations in this report are all run on a Xilinx ZC702 development board. The board is built to simplify prototyping with the Xilinx Zynq SoC (described in Section 5.1), as many peripherals are included. Here follows a summary of the most important peripherals/properties of the board, for this project: [31]
• Zynq-Z7020
Middle-end Zynq SoC
• USB to UART
Serial communication, makes terminal input/output possible. Useful for debugging
• USB to JTAG
Used for configuring the FPGA/programming the CPU as well as hardware/software debugging
• DDR3 RAM
Memory for both FPGA and CPU
• FMC connectors
Attaching additional peripherals such as the FMC Imageon card for additional HDMI input and output ports
• SD card reader
Non-volatile storage of programs, for example an OS
• Ethernet
Connecting to LAN or internet
The board itself has however only a single HDMI-out port, which is why an Imageon
extension card was attached, which though an FMC connector supplies an extra HDMI-
input and HDMI-output. This enables us to both receive a video stream on one cable,
then manipulate it and transmit the result on the other cable. Figure 17 shows the basic
setup of the ZC702 development board, with the FMC Imageon card supplying an input and output video stream to the ZC702 board.
FMC Imageon Card
ZC702 board
3 R D D
ZYNQ
