Phase unwrapping by N-connected TRW-S algorithm for InSAR images

(1)

by

Mehrnaz Movahed

B.Sc., University of Tehran, 2016

A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of

MASTER OF APPLIED SCIENCE

in the Department of Electrical and Computer Engineering

c

Mehrnaz Movahed, 2020 University of Victoria

(2)

Phase Unwrapping by N-Connected TRW-S Algorithm for InSAR Images by Mehrnaz Movahed B.Sc., University of Tehran, 2016 Supervisory Committee Dr. A. Baniasadi, Supervisor

(Department of Electrical and Computer Engineering)

Dr. N. J. Dimopoulos, Departmental Member

(3)

ABSTRACT

Synthetic aperture radar (SAR) data has been an interesting subject to study and investigate for researchers. One of the applications of SAR data is obtaining inter-ferometry synthetic aperture radar (InSAR) images [1]. However, these images as they are, are not useful for research purposes. Similar to other applications involving images, the presence of noise can distort a significant amount or all of the desired information. Various processes are applied to these images to extract the desired information [2]. One of the pieces of information that is extracted is phase image. Because of the periodic nature of phase, phase images only contain values between (-π,π]. Thus, resulting in an ambiguous phase image. Phase unwrapping is applied to recover the true value of the phase from the phase image [3].

Many algorithms have been proposed for phase unwrapping. Some of these al-gorithms, such as belief propagation (BP) and graph cuts, are energy minimization algorithms [4]. That is, they convert the phase unwrapping into energy minimization and find a lower bound on the energy function. In the process, these algorithms find the value of the real phase in the phase images. However, these algorithms suf-fer from various drawbacks such as limited number of applicable energy functions, getting stuck in a loop, or relatively low accuracy [5]. To overcome some of these drawbacks, sequential tree re-weighted algorithm (TRWS) was proposed. TRWS has better but still relatively low accuracy.

In this thesis, we improved the accuracy of TRWS. To achieve better accuracy we use more dense graphs (graphs with more edges) than what was originally proposed in TRWS. To compare the performance of our proposal with TRWS in terms of accuracy we tested our algorithm on synthetic images similar to actual InSAR images. The main reason for not using actual InSAR images is because there is no true value to compare the simulation results with. Other reasons include relatively huge image sizes and insufficient resources required to process them. Simulation results show that our proposal has more accuracy than TRWS. This shows that our proposal has a better performance than TRWS for phase unwrapping.

(4)

List of Tables

Table 4.1 The mean square values for Fig. 4.2 . . . 35 Table 4.2 The mean square values for Fig. 4.3 . . . 35 Table 4.3 The average of MSE of executing the algorithm for different

neigh-borhoods in 21 images with no noise. . . 35 Table 4.4 The MSE values of Fig. 4.4 . . . 36 Table 4.5 The MSE values of Fig. 4.5 . . . 36 Table 4.6 The average MSE of executing the algorithm for different

neigh-borhoods in 21 images for different noise magnitudes. . . 39 Table 4.7 The standard deviation of executing the algorithm for different

neighborhoods in 21 images for different noise magnitudes. . . . 39 Table 4.8 The confidence interval for different neighborhoods with different

noise levels. . . 40 Table 4.9 The average MSE of executing the algorithm for different

neigh-borhoods in 21 images with high density of fringes. . . 41 Table 4.10The overhead time of using OpenGM for each neighborhood. . . 43

(7)

List of Figures

Figure 1.1 Volcano Kilauea, Hawaii. Taken on October 4, 1994 and

pub-lished on May 1, 1999. Courtesy NASA/JPL-Caltech . . . 2

Figure 1.2 Izmit, Turkey 1999 earthquake interferogram . Courtesy NASA/JPL-Caltech . . . 4

Figure 1.3 Orbital effect for two satellites on two different orbital paths [14] 5 Figure 1.4 An everywhere smooth prior function [27]. . . 8

Figure 1.5 Piecewise constant prior function [27]. . . 8

Figure 1.6 Piecewise smooth prior function [27]. . . 9

Figure 2.1 All the cliques in a 3 × 3 grid . . . 14

Figure 3.1 A 4-connected neighborhood. The black circle is the main pixel and the white ones are the neighbors. . . 22

Figure 4.1 Examples of N-connected neighborhood . . . 29

Figure 4.2 First example of PU on an image with no noise. . . 33

Figure 4.3 Second example of PU on an image with no noise. . . 34

Figure 4.4 First example of PU on an noisy image. . . 37

Figure 4.5 Second example of PU on an noisy image. . . 38

Figure 4.6 Comparing images with different fringe density. . . 40

Figure 4.7 MSE/time diagrams of different noise levels for different neigh-borhoods . . . 42

(8)

Introduction

Interferometric synthetic aperture radar (InSAR) is a technique that uses radars in mapping the changes on the earth’s surface [1]. In this technique two or more synthetic aperture radar (SAR) images captured by remote sensing satellites are used. Various information are extracted from these images such as phase which is the interest of InSAR technique. By comparing the images and detecting the changes of phase, InSAR enables us to identify the changes of distance between the satellite and the earths surface. Any change in this distance shows a surface deformation which can be results of natural events like earthquakes and volcanoes, subsidence and uplift induced by nature incident or human intervention, glacial motions, etc. InSAR technique has several applications like generating digital elevation map (DEM), natural hazards studies like volcanoes and earthquakes, monitoring structural stability, infrastructure, and building monitoring, monitoring landscape features, etc [6]– [9]. In this chapter, we provide a brief introduction to SAR images, interferometry, and various processes used to extract useful information from them.

1.1 Synthetic aperture radar

Synthetic aperture radars (SAR) illuminate a target on the surface of the earth by producing a narrow beam of electromagnetic waves successively within the range of microwave frequency. They generate amplitude and phase images based on the reflected signal. The atmospheric absorbance of electromagnetic waves within the microwave frequency range is low. Therefore, the satellite can take images in the presence of any atmosphere media such as clouds, fog or smoke and at any time of

(9)

day or night.

The phase and amplitude images show different features of the captured terrain. The amplitude image is used to measure the roughness and the slope of the surface. It shows the reflectivity of the target terrain, i.e. how much of the electromagnetic waves are sent back to the satellite [6]. For example, amplitude images of water are shown as dark spots. This is because water is a great reflector and reflects the signal away from the receivers. Thus, nothing is recorded on the amplitude image. Fig. 1.1 shows an amplitude image of the volcano Kilauea, Hawaii.

Figure 1.1: Volcano Kilauea, Hawaii. Taken on October 4, 1994 and published on May 1, 1999. Courtesy NASA/JPL-Caltech

The phase image shows the distance between the satellite and the terrain. There are different scatterers on the surface in different distances from the radar with dif-ferent slant ranges. Each cause a delay during the two-way travel of the signal from and back to the satellite. Since, the transmitted signal is purely sinusoidal, the phase difference, denoted by φ, between the sent and received signal is equivalent to the delay. Thus, φ is proportional to the two-way travel distance of the signal, denoted by 2R through the following,

φ = 2π λ 2R =

4π λ R,

where λ is the wavelength of the signal [10]. Because of the periodic nature of the signal, two travel distances that are an integer multiple of 2π radians far from each

(10)

other, have the same phase. In other words, the measured phase only shows the fraction of the distance that is less than one wavelength. There are different factors that affect the phase of the received signal [10]. For example, one factor is the huge difference between a signal wavelength and the dimension of a target on the earths surface known as a resolution cell. The dimension of the resolution cell is a few meters and the wavelength is just a few centimeters (e.g. for European Remote Sensing satellite λ is ≈ 5.6 cm). Another factor is numerous scatterers in a resolution cell. The satellite receives the signal after multiple reflections between scatterers. In addition, the structure and properties that the scatterers are made of effect its reflectivity and as a result change the phase. Because of these reasons and many more, the phases in a phase image appear to be random with no correlation. Hence, a single SAR phase image is not very useful in visualizing the topography of the target terrain. Having more than one image and comparing them, may help in removing some of these factors.

1.2 Interferometry

As previously mentioned, isolating the factors that result in a lack of correlation in a phase image will lead to a meaningful result. The methods of isolating these factors is known as interferometry. In interferometry, the difference of two phase images from the same terrain, known as an interferogram, is produced. The phase difference is measured in radians and is shown in the form of colored contours (fringes) covering a single 2π cycle. The color of each contour represents a certain phase within this cycle. Thus, the phase of interferograms is referred to as wrapped. In order to find the true phase, phase unwrapping is required. Fig. 1.2 shows an interferogram of Izmit, Turkey, created by two images taken on August 13, 1999 and September 17, 1999. This interferogram shows the terrain deformation caused by an earthquake on August 17, 1999. Each contour represents a displacement on the surface of the earth. Note that in interferograms the phase difference is shown in the form of repeated contours. In this example (Fig. 1.2) the contours show 28 mm vertical motion (towards the satellite) or about 70 mm horizontal motion. The red line in Fig. 1.2 represents the location of fault breaks.

There are two methods for obtaining the images required to find an interferogram image. The first method is across-track interferometry. In this method, two images of a terrain are taken at the same time but at slightly different positions. The second

(11)

Figure 1.2: Izmit, Turkey 1999 earthquake interferogram . Courtesy NASA/JPL-Caltech

method is along-track interferometry. In this method, two images of a terrain are taken at two different times but at the same position. Across-track interferometry is used in obtaining surface topography. Along-track interferometry is used to measure surface deformation at different time intervals. Along-track interferometry is used in generating DEMs.

Factors that contribute to producing erroneous interferograms include:

• Orbital effect : To find a useful interferogram in along-track interferometry, the position of the second satellite, taking the second image of the terrain, should be as close as possible to the same spatial position of the first satellite. But this rarely happens. There is always a slight difference in position of the first and second satellite. This is caused by using two different satellites orbiting in two different orbiting paths or by a slight change in a satellites orbital path around the earth. This is known as orbital effect. Fig. 1.3 shows the orbital effect for two satellites orbiting on two different paths. This effect can be modeled and removed from the interferogram [13].

• Thermal effect : When images are obtained at different times, the temperature of the sensors and the weather conditions are probably different in two image acquisitions. This is known as thermal effect and produces a phase noise in the interferogram. This effect can be modeled and removed from the interferogram [13].

(12)

Figure 1.3: Orbital effect for two satellites on two different orbital paths [14]

surface that have occurred between acquisitions such as soil moisture, changes in vegetation, the roughness of the surface, etc. This effect like other previously mentioned effects can be modeled and removed [13].

1.3 Phase unwrapping

A phase of a signal is usually measured within the interval (-π,π]. If the phase is not within this interval it is wrapped by removing integer multiples of 2π until it falls within (-π,π]. Phase unwrapping (PU) is the technique of recovering the 2π integer multiples of the wrapped phase.This technique has various applications such as InSAR [11], 3D surface measurement [15], digital holography [16], and magnetic resonance imaging (MRI) [17]. In the past two decades, the concept of PU has been widely studied and many algorithms have been developed for this task [18]. In this thesis PU algorithms on two dimensional (2D) images are analyzed. There are two classes of 2D phase unwrapping methods:

• Local PU methods : Local PU methods are based on path following paradigms. That is they start from a certain point (pixel) on the image and follow a pre-defined path. The rate of success depends on which path has been chosen. If the starting pixel is not surrounded by many fringes, the difference between the obtained phase and the true phase (known as unwrapping error) is mini-mized [19]. Different approaches have been proposed to find a better path for local PU methods [20]– [24]. Since there are multiple paths to be chosen within an image, local PU methods may not give a unique result [25]. Based on the

(13)

image and the density of fringes, there can be more than one solution. This becomes a concern when the images contain noise such as InSAR data. That is because noise is an effective factor in choosing a path.

• Global PU methods : In global methods the phase unwrapping problem is formulated in a generalized minimum-norm sense [18]. These methods take the whole image into account and are independent of a path The approach in global PU methods is to minimize a cost function in the form of Lp _{norm function [18].}

The goal is to find the absolute phase of the image by finding the minimum of the Lp _{norm of the differential value of wrapped phase and unwrapped phase}

derivatives. The general form of Lp norm minimization function is shown below [18]: min    M −2 X i=0 N −1 X j=0

|φi+1,j − φi,j− ∆xi,j| p + M −1 X i=0 N −2 X j=0 |φi,j+1− φi,j− ∆ y i,j| p    , (1.1)

where φi,j = ψi,j+ 2ki,jπ, (i, j) is the coordination of a pixel, ψi,j is the wrapped

phase, φi,j is the unwrapped phase, and ki,j is the desired integer. In Eq. 1.1

∆x_i,j and ∆y_i,j are wrapped phase differences that can be computed as below [18]: ∆x_i,j = Wψi+1,j − ψi,j , i = 0, ..., M − 2, j = 0, ..., N − 1

∆x_i,j = 0 otherwise, and

∆y_i,j = Wψi,j+1− ψi,j , i = 0, ..., M − 1, j = 0, ..., N − 2

∆y_i,j = 0 otherwise,

W is a wrapping operator that wraps its argument into the range of (-π,π].

1.4 Energy minimization

One of the problems in early vision (low-level computer vision) is estimating a local quantity (known as a label) of a pixel in data combined with noise [26]. In stereo and motion problem, and image restoration this quantity represents disparity and intensity, respectively. This quantity tends to change smoothly piecewise, i.e. there

(14)

is no sudden or rapid change from one pixel to the next. However, in some cases for example over the object boundaries, there is a significant change in quantity from one pixel to an adjacent one. Let L and P denote a set of integers (refer to labels) and a set of pixels, respectively. The solution to this problem is to find a labeling system such as f that assigns a label, denoted by fp, to each pixel p ∈ P such that the

changes between a pixel and its neighbors stay smooth. The outcome of the labeling system should be consistent with the input, i.e. should reflect the actual changes in the measured data.

Pixel labeling problems can be formulated as energy minimization problems [26]. An energy function has two terms; one term maintains consistency between the result and the observed data and the second term maintains smooth changes between a pixel and its neighbors. The framework of formalizing a labeling system f into an energy minimization problem is shown in [26]. For self sufficiency of this thesis this framework is provided below. Let E(f ) be an energy function. Hence

E (f ) = Esmooth(f ) + Edata(f ) .

In an energy minimization problem the goal is to find the minimum of E(f ). This im-plies minimizing Edata(f ) and Esmooth(f ). Minimizing Edata(f ) means minimizing the

inconsistency between the result and the observed data. Minimizing Esmooth(f ) means

minimizing the sudden changes that reduce the smoothness of the result. Edata(f ) is

in the form shown below [26]

Edata(f ) =

X

p∈P

Dp fp , (1.2)

where Dp fp is a measurement of how suitable a label fp is for pixel p based on the

observed data. A general form Esmooth(f ) is [27]:

Esmooth(f ) =

X

{p,q}∈P

Vpq fp, fq , (1.3)

where Vpq(fp, fq) measures the smoothness between a pixel p and its neighbor q with

labels fp and fq, respectively. The outcome of the function Vpq(fp, fq) is the cost

based on the difference of labels p and q.

In visual problems having large discontinuities in the image is common. If the smoothness function enforces smoothness throughout the entire image it does not

(15)

let a large discontinuity show itself in the result. This type of smoothness function is called an everywhere smooth prior function [27]. Fig. 1.4 shows an example of an everywhere smooth prior function. This figure shows that when |fp − fq| (the

horizontal axis) is small the cost (the vertical axis) is a small value. However for larger values of |fp− fq|the cost becomes larger without any limit.

Figure 1.4: An everywhere smooth prior function [27].

Discontinuity-preserving energy functions are functions that are bounded on the cost value. There are two types of discontinuity-preserving energy functions. The first are known as piecewise constant prior function. Fig. 1.5 shows an example of this type. It can be seen that if |fp − fq| is zero, the value of cost function is zero.

Otherwise, the cost is a constant value. This type of function is more suitable for images with large discontinuities [27]. The second type of discontinuity-preserving

Figure 1.5: Piecewise constant prior function [27].

(16)

for small values of |fp − fq| is small. For larger values of |fp − fq| up to a certain

point the function is linear. For values beyond that certain point, the cost function is constant. Fig. 1.6 shows a piecewise smooth function. This function has shown better results than the other two types [27].

Figure 1.6: Piecewise smooth prior function [27].

As previously mentioned, one of the applications of the energy minimization tech-nique is the pixel labeling problem. PU is a pixel labeling problem and can be solved in the energy minimization framework [31]. This is a generalization of the global PU method. The goal of global phase unwrapping methods is to minimize a cost function in the form of an Lp norm function (Eq. 1.1) in order to find suitable labels for the pixels and thus obtain their true phase. An energy minimization framework with piecewise smooth prior smoothness function is used to reach this goal. In the next chapters, we will provide more details on the mathematical aspect of this problem and its solution.

1.5 Related works

Phase unwrapping is an easy task when there is no noise or aliasing problem. In reality, InSAR data sets always contain a large amount of noise, and even the un-noisy areas can be affected by the contaminated areas. That will result in inaccurate PU output with a high error percentage. Therefore, this process becomes more com-plicated and finding a suitable method is a challenge. There have been numerous researches on finding new algorithms for the PU problem or improving the previous ones from different aspects like efficiency, accuracy, memory usage in the case of large grid phase unwrapping, etc.

(17)

In [28], the authors propose a new discrete energy minimization framework for the PU problem. The objective function is a first order Markov random field where the corresponding energy has only two-way interactions. The energy minimization is carried out via a two iterative scheme similar to [29] and [21]. In this method each step is projected onto a graph min-cut/max-flow problem, mainly relying on [30], [27], and [26]. When the clique potentials that are chosen are convex, a low order polynomial algorithm is devised (specifically, the classical Lp _{norm when p ≥ 1).}

On the other hand, for nonconvex clique potentials, the proposed problem is NP-hard ( [30] and [27]). In this setting, neither the iterative scheme, nor the proposed projecting onto a graph min-cut/max-flow method is possible. Therefore, the authors propose a discontinuity preserving algorithm. This is a modified version of the original projecting onto a graph min-cut/max-flow algorithm based on two main ideas: each binary problem has a wide configuration space and applying major minimization concepts to the energy function.

In [32] the authors propose another algorithm based on energy minimization to solve the PU problem. In this algorithm they only consider the energy of the four close neighbor pixels. The algorithm counts the energy difference of these pixels and then it calculates the probability to accept or reject the new labels for the pixels. The same approach is used by considering the four pixels in the third level (the neighbors that are three pixels far from the main pixel in every direction) in [31].

1.6 Contributions and outline

In this chapter, we provided a brief introduction to InSAR technique and its appli-cations. One of these applications is generating DEM. This is a map that shows the changes on a surface of the earth during a time period. One of the challenges in determining the changes on a surface is identifying the true phase. This is because the phase part of interferograms which shows changes in distance, is wrapped. PU is a method for obtaining the true phase from the wrapped one. We explained how PU can be solved in the framework of energy minimization. That is, by defining appropriate energy functions and applying minimization technique on them to solve the pixel labeling problem. Since PU is a pixel labeling problem, hence we can use the aforementioned energy functions to obtain the unwrapped phase.

In the second chapter, we will provide the background knowledge to understand energy minimization technique. Specifically, in this thesis tree re-weighted message

(18)

passing (TRWS) algorithm is of interest. In this chapter, we discuss concepts such as markov random field (MRF) and maximum a posteriori (MAP) estimation. MRF is a means of modeling the interactions of pixels in space. But since these interactions can not be obtained directly, MAP is used to estimate them. By the end of this chapter, it will be easier to understand how the PU problem is formulated into an energy minimization problem.

In Chapter 3, we will present the TRWS algorithm. This algorithm is used in solving the energy minimization problem. This algorithm works by setting a lower bound equal to zero and increasing it in each iteration until the lower bound meets the minimum of the energy function. The advantage of this algorithm over other algorithms is that the lower bound always increases (i.e. never decreases) and thus guarantees the lower bound will meet the minimum of the energy function. In addi-tion, TRWS uses half the memory that other minimization algorithms such as belief propagation use. However, TRWS is not completely robust to noise. This means that in the presence of noise from the sources that were mentioned in this chapter, the minimum of the energy function found by TRWS may not be the correct answer. The advantages that we mentioned are the reasons that we chose this algorithm to be improved in this thesis.

In Chapter 4, we will present an enhancement to minimize the effect of noise on TRWS. This improvement is based on utilizing N-connected neighborhood rather than 4-connected neighborhood that is used in TRWS. To determine the phase of a pixel, TRWS uses the information from that pixel and 4 surrounding pixels. This is referred to a 4-connected neighborhood. In this chapter, to determine the phase of a pixel, we use the information of N pixels surrounding it. Thus, we extend 4-connected neighborhood to N-connected neighborhood. Simulation results and calculating the mean square error between N-connected TRWS and 4-connected TRWS will show that this approach results in higher robustness against noise. Finally, in Chapter 5, we will provide a summary of this thesis and other possible approaches to further improve N-connected TRWS.

(19)

Chapter 2 Preliminaries

The focus of this thesis is improving TRWS for solving the PU problem. As previously mentioned, TRWS is a message passing algorithm, similar to belief propagation, used in energy minimization. An energy function has two terms. The first term is known as a unary term (i.e. data energy function). The second term is a pairwise interaction term (i.e. smoothness energy function). This energy function is derived from a Markov random field (MRF) [33]. An MRF is used to obtain a configuration of pixels and labels by applying maximum a-posteriori (MAP) estimation [27]. In this chapter, we provide the preliminaries required to understand the energy minimization framework and the TRWS algorithm.

2.1 Message passing algorithms

Message passing algorithms are used in solving various problems including energy minimization, optimization, constraint satisfaction, and inference [34]. There are different classes of message passing algorithms. The two most well known of them are belief propagation (BP) and divide and conquer (DC). These algorithms are based on sending messages between nodes. There are different graphical models that message passing algorithms use to define and visualize the problem (such as factor graphs, MRF, etc) [34].

In energy minimization problem, there is a set of random variables that we call labels. Every node has a number of possible labels. Therefore, there are many possible configurations for the nodes of a graph. However, there is one configuration that minimizes the energy function that has been defined for the problem. Message

(20)

passing algorithms propagate messages between the nodes and pass the energy or cost of each node to the next one. The algorithms update the messages once they have passed them through the nodes. By passing and updating, the message passing algorithms try to decrease the total energy of the graph and find a minimum.

Some of the graphical models (e.g factor graph) includes loops. In energy mini-mization problems, the presence of loops sometimes causes uncertainty for an algo-rithm to terminate the process. An MRF model is an undirected graph that makes it possible to follow a determined direction in order to reach convergence and termi-nate the algorithm. It is a well known graphical model for solving computer vision problems. The TRWS algorithm is based on the MRF model. This algorithm does not include loops and pass the messages through trees which are undirected. That is why TRWS guarantees convergence. In the next section, we will explain MRF and we will provide details about how we derive the energy function from an MRF model.

2.2 Markov random fields

Suppose P is a group of pixels in a graph G with a set of labels L . Let Np denote

the pixels that are connected to p ∈ P. Np is known as neighborhood of p. For any

p, q ∈ P, Np and Nq have the following properties:

(a) p /∈ Np, q /∈ Nq,

(b) if p ∈ Nq then q ∈ Np.

Suppose F = F1, F2, . . . , FN is a set of random variables on P where n = |P| is

the cardinality of P. Each random variable Fi, 1 6 i 6 n, takes a value from

L = f1, f2, · · ·. Let f = fp|p ∈ P be a configuration on F and F is the set of all

configurations. In this thesis, P Fi = fp which is the probability of random variable

Fi being fpwill be shown as P fp. For F to be an MRF it must satisfy the following:

(i) P (f ) > 0 ∀f ∈ F , (ii) P fp|{fq}q∈P\p

= P fp|{fq}q∈Np .

The second criteria means that the assigned label for each pixel depends only on the labels of its neighbors.

There are two methods to define an MRF. The first method is by their joint probability distribution and the second method is by the local conditional probability.

(21)

In the latter method, there are non-trivial consistency constraints that make this method difficult [27]. Hence, the first method is the most commonly used. The Hammersley-Clifford theorem provides a method to characterize an MRF [35]. This theorem shows that MRFs and Gibbs random fields are equivalent. To define a Gibbs random fields, it is necessary to define a clique. A subset c of P is a clique if any member of c is a neighbor of all other members. C denotes the set of all the cliques. Fig. 2.1 shows all the available cliques in a 3 × 3 grid. A Gibbs random field can be

Figure 2.1: All the cliques in a 3 × 3 grid defined as follows [27] P (f ) = Z−1· exp −X c∈C Vc(f ) , (2.1)

where Z is a constant with the purpose of normalizing and {Vc(f )} are clique potential

functions. If a clique contains more than two nodes, then Vc(f ) = 0. If there are two

nodes, p and q, in a clique with labels fp and fq, respectively, then Vc(f ) is defined

as:

Vc(f ) = V{p,q} fp, fq , (2.2)

where V{p,q} fp, fq is the same as Section 1.4. Therefore, equation 2.1 is:

P (f ) = Z−1· exp − X {p,q}∈N V{p,q} fp, fq , (2.3)

(22)

thus defining the probability of a configuration in MRF. Since there are many possible configurations, it is an NP-hard problem to find the most probable f ∈ F on F . Therefore, we need to use likelihood of f based on the observed data to form a posterior probability. Maximizing the posterior probability (MAP estimation) gives us the optimum answer.

2.3 Maximum a posteriori estimation

The posterior probability is defined as the probability of configuration f occurring given that an observation d has happened. Using Bayes theorem

P f |d = P d|f P (f )

P (d) . (2.4)

MAP estimation is finding the configuration f∗ that maximizes P f |d. Since d is the observed data and P (d) is known, then:

f∗ = arg max

f ∈F

P d|f P (f ) (2.5)

A suitable model of P d|f will be ( [27]): P d|f

= Y

p∈P

P dp|fp , (2.6)

where dp is the observed data of pixel p. Suppose the following as P dp|fp definition

( [27]): P dp|fp = Cp· exp −Dp fp fp ∈ L, (2.7)

where Cp is a constant with the purpose of normalizing and Dp fp is as defined in

Section 1.4. From equations 2.6 and 2.7 the following can be concluded:

P d|f ∝ exp −X p∈P Dp fp . (2.8)

From equations 2.3 and 2.8, the definition f∗ will be:

f∗ = arg max f ∈F exp − X {p,q}∈N V{p,q} fp, fq − X p∈P Dp fp . (2.9)

(23)

Maximizing the expression in Eq. 2.9 is equivalent to minimizing the exponent. Hence, if we define the exponent as E (f ) (i.e. an energy function),

E (f ) = X {p,q}∈N V{p,q} fp, fq + X p∈P Dp fp , (2.10)

then finding a maximum of a posteriori probability and thus an optimum configuration is equivalent to an energy minimization problem.

2.4 Summary

In this chapter we provided the background concepts that help in understanding the energy minimization technique and the TRWS algorithm. We explained how the TRWS algorithm uses message passing method to solve an energy minimization problem. Message passing algorithms use graphical models to model the problem and solve it. We presented MRF which is a graphical model that is very well known for modeling computer vision problems. The configuration that MRF provides can not be obtained easily. Thus, we explained MAP estimation to obtain an approximation of the configuration. At the end of this chapter, we showed how the pixel labeling problem becomes an energy minimization problem.

In the next chapter we will present the TRWS algorithm. We will provide the details of this algorithm and explain how we use this algorithm to solve the PU problem.

(24)

Chapter 3 Sequential tree re-weighted

message passing algorithm

Sequential tree re-weighted message passing algorithm (TRWS), is an improved ver-sion of tree re-weighted message passing algorithm (TRW) and is an energy minimiza-tion algorithm [5]. As previously menminimiza-tioned in Chapter 1, the advantage of TRWS over TRW is that the lower bound always increases (i.e. never decreases). Thus, it guarantees the lower bound will meet the minimum of the energy function [5]. Hence, unlike TRW, TRWS always converges. TRWS shows promising results for pixel labeling problems in an MRF setting [4]. As previously mentioned, PU is a pixel labeling problem. In this thesis TRWS is exploited to solve this problem for InSAR applications. In this chapter we provide a summary of the TRWS algorithm.

3.1 Introduction

As shown in Chapter 2, many early vision problems can be presented in the form of an energy minimization problem where the energy function is in the form of Eq. 2.10. Finding the minimum of an energy function is an NP-hard problem [5]. Thus, the focus has been on approximate minimization algorithms. Graph cuts algorithm is one of the best in this area [5]. In many vision applications the results of graph cuts algorithm are considered the most accurate amongst other energy minimization algorithms [5]. The disadvantage of graph cuts algorithm is that it is only applicable to a limited number of energy functions [36].

(25)

func-tion in the form of Eq. 2.10 max-product BP is applicable. However, it does suffer from some drawbacks [5]. One drawback is that in some cases it gets stuck in a loop and never converges [5]. Another drawback is that for some functions where both max product BP and graph cuts are applicable, max product BP returns a result greater than that of graph cuts algorithm.

In [37], a new energy minimization algorithm called TRW is introduced. The idea of TRW is maximizing a concave lower bound. When a certain condition, known as the tree agreement, is satisfied the algorithm returns the global minimum of the energy function. Similar to max product BP, TRW is applicable to any energy function in the form of Eq. 2.10. One advantage of TRW over max product BP and graph cuts is it returns lower result for an energy function. However, it suffers from some drawbacks. One of the drawbacks is it can get stuck in a loop and not converge. Another drawback is that TRW sometimes decreases the lower bound resulting in unnecessary execution of the algorithm or not finding the minimum of the energy function [5].

TRWS is a modified version of TRW that guarantees no decrease in lower bound on the energy function. This is achieved by defining a condition known as weak tree agreement (WTA). The algorithm always converges to a result that satisfies this condition [5]. In the rest of this chapter, we explain the TRWS algorithm.

3.1.1 Notations

In order to explain how TRWS works, it is required to introduce some notations. These notations are taken from [5]. Let G = (V, E ) denote the graph of pixels obtained from an image where V is the set of nodes and E is the set of edges. The energy function in Eq. 2.10 can be rewritten as

E x | θ = X s∈V θs(xs) + X (s,t)∈E θst(xs, xt) , (3.1)

where s and t are two nodes of the graph and xs and xt are the labels of these two

nodes. θs(·) is a unary data energy (same as Dp) and θst(·, ·) is a pairwise interaction

energy (same as V{p,q}). The vector x is a vector of length n = |V| where each

element is the label of the n nodes of G. For a node s, the variable xs ∈ Xs is a

possible label where Xs is the set of all possible labels for node s. Therefore, x∈ X

(26)

If we denote θs(xs = j) and θst(xs= j, xt= k) by θs;jand θst;jk, respectively, then

θ = {θα|α ∈ I} where

I = {(s; j)} ∪ {(st; jk)}.

It is shown in [5] that an energy function can be written as the inner product of two vectors dependent on θ and x. Hence, θ and x should be in the same space. This is why in [5], the mapping φ : X → Rd is introduced. With this mapping, the energy function in Eq. 3.1 can be written as

E x | θ = hθ, φ (x)i = X

α∈I

θαφα(x) ,

where φα :X → R consists of the following functions:

φs;j(x) = [xs = j]

φst;jk(x) = [xs = j, xt = k],

and [.] is one if its argument is true, and zero otherwise. Other functions that are required in explaining the TRWS algorithm are the three functions that represent the minimum of energy under certain conditions. These functions are:

Φ (θ) = minx∈X E x|θ Φs;j(θ) = minx∈X ,xs=j E x|θ Φst;jk(θ) = minx∈X ,xs=j,xt=k E x|θ

Φs;j is called min-marginal for node s and Φst;jk is the min-marginal for edge (s, t ).

As apparent from the name of the algorithm, TRWS operates on a set of trees. Let T be the set of all trees in graph G and T denote a tree in this set. Let ρT _be

the probability that T is chosen from T . Note that ρT 6= 0 for all T ∈ T and every edge of G belongs to at least one tree in T . Suppose,

IT _{= {(s; j) |s ∈ V}T_{} ∪ {(st; jk) | (s, t) ∈ E}T_},

(27)

energy parameter associated to T [5]. Thus, AT _{can be defined as [5]:}

AT _{= {θ}T

∈ Rd_|θT

α = 0 ∀ α ∈ I\I T_}.

Suppose θ is the collection of all energy parameters of all trees in T , then A can be defined as [5]:

A = {θ ∈ Rd×|T |_{| θ}T _{∈ A}T _{f or all} _{T ∈ T }.}

Using the functions and notations that are defined the lower bound of the energy function, Φρ, is as follows [5]: Φρ(θ) = X T ρTΦθT = X T ρT min x∈Xhθ T_{, φ (x)i.}

Hence, an energy minimization problem using TRWS is equivalent to maximizing Φρ

which is a concave function of θ [5], that is: max

θ∈A,P

Tρ

T_θT_{= ¯}_θ Φρ(θ) . (3.2)

Note that the obtained lower bound of the energy function, Φρ, gives the optimal

value for vector ¯θ =P

T ρ

T_θT _[5].

3.2 The TRWS algorithm

In [5] three different versions of the algorithm are described. One of them is shown to be the efficient implementation of TRWS. This is the algorithm that will be described in this thesis. The reason this algorithm is more efficient compared to the other two is that it requires less memory.

Before the first step of the TRWS algorithm, the following needs to be done [5]: • Assigning a number to each node on the graph from {1, 2, · · · , |V|} to have an

order for nodes. This order is denoted by i (.).

• Selecting trees from T based on i (.). These trees should be monotonic and cover all edges of graph G. A graph G and its trees are called monotonic if an ordering i (u) , u ∈ V, of nodes exists such that for consecutive nodes uT

1, · · · , uTn(T ) of any tree T ∈ T , i uT 1 , · · · , i uT n(T ) is a monotonic sequence.

(28)

• Choosing a probability distribution,ρ, for T such that ρT _{> 0 for every T ∈ T}

and P

T ρ

T _{= 1.}

Now we can run the algorithm. Note that the input energy function is shown as ¯θ. Suppose, for nodes s, t ∈ V, the message of directed edge (s, t) ∈ E is shown by Mst.

This message is a vector and its number of elements is equal to |Xt| (i.e. the number

of possible labels for destination node). The steps of the TRWS algorithm are as follows [5]:

1. Initialize all the messages on edges of G to zero.

2. Set Ebound = 0. For nodes s ∈ V in the order of increasing i (s) do the following:

– Compute ˆθs = ¯θs +

P

(u,s)∈EMus. Normalize vector ˆθs using:

δ := min

j

ˆ

θs;j, θˆs;j := ˆθs;j− δ, Ebound := Ebound+ δ.

– For each edge (s, t) ∈ E (i (s) < i (t)), message Mst should be updated and

normalized using in the following: Mst;k:= min j { γstθˆs;j− Mts;j + ¯θst;jk} δ := min k Mst;k, Mst;k = Mst;k− δ, Ebound = Ebound+ δ.

3. Set i (u) := |V| + 1 − i (u).

4. Check the stopping criterion. If it is satisfied the algorithm is done and otherwise back to step 2.

The coefficient γst used in step 2 is defined as γst = ρst/ρs where ρst and ρs are the

probabilities of a random chosen tree T having edge (s, t) and node s, respectively. It means that γst is the probability that a randomly chosen tree has edge (s, t) given that

it contains node s. For further clarification, it is better to explain how the monotonic trees are constructed. In [5], it is described that the monotonic trees are chosen in a way that for the first node s, an edge (u, s) such that u has lower ordering than s (i (u) < i (s)) can not be found. A similar statement holds for the last node. For the last node t, there is no edge (t, v) such that v has a higher ordering than t. For every tree that is constructed, the edges forming that tree are removed from G. In the end,

(29)

there should be no edges left. Thus, it can be concluded that every edge is covered by only one tree. The number of trees that node s appears in is [5]

ns= max{| (u, s) ∈ E : i (u) < i (s) |, | (s, v) ∈ E : i (v) > i (s) |}.

Thus, γst = 1/ns. For example, consider a 4-connected neighborhood shown in Fig.

3.1, in a neighborhood for every node s, there are two nodes with higher ordering. Hence, ns= 2 and γst = 0.5.

Figure 3.1: A 4-connected neighborhood. The black circle is the main pixel and the white ones are the neighbors.

The TRWS algorithm requires half the memory space that BP does. In BP, for every edge, both forward and reverse messages are stored. However, in TRWS there is no need to store the reverse messages because they are already updated when they need to be used [5]. Hence, the same memory space can be used to keep either Mst

or Mts. After edge (s, t), is processed in step 2 of the algorithm, when the current

node is s, we replace Mts with Mst.

3.2.1 Weak tree agreement

The TRWS algorithm does not specify a stopping criterion [5]. In this section, we describe WTA as a condition that needs to be satisfied for TRWS to stop running. Some methods for checking whether WTA has been satisfied is provided.

WTA guarantees that TRWS will never get stuck in a loop. To define WTA a con-sistent set needs to be defined. Suppose OPTT θT_{is a set of optimal configurations}

with the parameter θT _{. Let OPT}T_{(θ) be the set {OPT}T _θT_{| T ∈ T }. This set is}

the set of all optimal configurations for vectors θT and OPTT (θ) ⊆ Xn [5]. For two sets S, ˜S ∈ Xn, S ⊆ ˜S if ST ⊆ ˜ST for every tree T. Now we can define a consistent set. A set S ∈ Xn is consistent if it satisfies the following three conditions [5]:

(30)

• For every tree T, ST _{6= ∅.}

• If node s is on tree T and T0_{, for every configuration x ∈ S}T _{there exists a}

configuration x0 ∈ ST0

such that xs = x0s.

• If edge (s, t) is in trees T and T0_{, for every configuration x ∈ S}T _{there exists a}

configuration x0 ∈ ST0 _{such that x}

s = x0s, xt= x0t.

At this point we can define WTA condition: vector θ = {θT_{} ∈ A satisfies the WTA}

if there exists a consistent S ⊆ OPT (θ).

Based on the definition of WTA, one method of determining if it is satisfied or not is to add more minimum configurations to S until it becomes consistent or there are no more configurations remaining. This is a time consuming method [5]. A simpler method is to check the lower bound (Ebound). If there is no increase for a certain

number of iterations (e.g. 10 iterations) in Ebound, the WTA condition is satisfied.

Thus, the algorithm can be terminated. Another method is to predefine the number of iterations and not let the algorithm run beyond them. Obviously choosing each method results in different speed and accuracy.

3.3 Using TRWS in phase unwrapping

As previously mentioned, the focus of this thesis is on the PU problem. Since PU is a pixel labeling problem, the TRWS algorithm can be used to solve it. In this section, the equations of Chapters 1 and 2 are matched with the steps of the TRWS algorithm. Recall from Chapter 1, in a PU problem the input is an ambiguous phase image where the phase of a pixel is in (−π, +π]. In other words, it is wrapped and does not show the true phase. If φs and ψs are the unwrapped and wrapped phase of

pixel s, respectively,

φs = ψs+ 2ksπ,

where ks is an integer known as the label. In Eq. 3.1 the label for pixel s is denoted

as xs.

The energy that we try to minimize in the TRWS algorithm consists of the smooth-ness energy (V{p,q} fp, fq) and the data energy (Dp fp) in Eq. 2.10. The

(31)

smooth-ness energy for a N × M image is derived from Eq. 1.1 with p = 1: Esmoothness= M −2 X i=0 N −1 X j=0

φi+1,j − φi,j − ∆xi,j

+ M −1 X i=0 N −2 X j=0

φi,j+1− φi,j − ∆yi,j

, where    ∆x

i,j = Wψi+1,j − ψi,j , i = 0, ..., M − 2, j = 0, ..., N − 1,

∆x i,j = 0, otherwise, and   

∆y_i,j = Wψi,j+1− ψi,j , i = 0, ..., M − 1, j = 0, ..., N − 2,

∆y_i,j = 0, otherwise.

This is the smoothness energy for a 4-connected neighborhood because the only neigh-bors that matter are one pixel far from the current pixel. Now if we replace φi,j with

ψi,j + 2ki,jπ: Esmoothness= M −2 X i=0 N −1 X j=0

ψi+1,j + 2πki+1,j − ψi,j+ 2πki,j − ∆xi,j

+ M −1 X i=0 N −2 X j=0

ψi,j+1+ 2πki,j+1 − ψi,j + 2πki,j − ∆yi,j

.

Assume that Ni+1,j = ψi+1,j − ψi,j− ∆xi,j and Ni,j+1 = ψi,j+1− ψi,j − ∆yi,j. Since 2π

is a constant we can divide by 2π to obtain:

Esmoothness= M −2 X i=0 N −1 X j=0

ki+1,j− ki,j + Ni+1,j/2π + M −1 X i=0 N −2 X j=0

ki,j+1− ki,j+ Ni,j+1/2π .

(3.3) As previously mentioned in Section 3.2, the second step of the TRWS algorithm involves ¯θst;jk in updating the messages. Let s = (i, j) and t = (i + 1, j), then

¯

θst;jk = Ni+1,j and is called a shift. An array of the shift for all pixels and

differ-ent combinations of possible labels can be computed before starting the algorithm. The goal of minimizing the smoothness energy function is to minimize the difference between the labels of adjacent pixels and consequently minimizing the difference be-tween the phase of adjacent pixels after unwrapping. However, we can minimize the

(32)

smoothness energy with different sets of labels as long as the difference between every two labels stays the same. Hence, there are many possible combinations and it is nearly impossible to find one answer. This is where the data energy makes solving the problem easier. The data energy function is the summation of a set of labels. Thus, minimizing the total energy function which includes the smoothness and the data energy functions, puts a limit to the number of sets of labels as possible answers. Obtaining the data energy is easier than the smoothness energy. The data energy is

Edata =

X

s∈V

ks− βs, (3.4)

where βs is the label for pixel s that we desire. Sometimes we have a prior knowledge

of the terrain that the phase images are captured from, or we want the phase difference between pixels in the final result to be within a specific range. In this case, we can use βs to manipulate the final result. In this project, we conduct our experiments

(the details are provided in Chapter 4) on synthetic interferograms. These images are random (they do not correspond to a specific area) and we do not expect the results to be in any specific range. Hence we consider βs = 0 for all pixels.

Based on Section 3.2 in every iteration of TRWS, the messages are updated and a new lower bound for the energy function is calculated. This means that in each iteration a new configuration and set of labels is derived. Hence in each iteration, new values for smoothness and data energy are obtained. The optimum point is where the lower bound meets the total energy. In that case, we terminate the algorithm and the set of labels that has caused this condition is the optimum answer. However, because of the time limitation that we may have, we define a total number of iterations where we terminate the algorithm thereafter even if we have not found the optimum answer.

3.4 Summary

In this chapter, we explained the TRWS algorithm. We compared the TRWS algo-rithm with BP, graph cuts, and TRW algoalgo-rithms and pointed out the advantages of using the TRWS algorithm over others. We provided some notations in order to explain the TRWS algorithm steps. After that, we explained WTA as a condition that needs to be satisfied for the algorithm to terminate. To determine if WTA has been satisfied, we provided some methods. Finally, in the last section, we reviewed the PU problem and how we can use the TRWS algorithm to solve it.

(33)

In the next chapter, we will explain our contribution in this thesis. We will provide the details of extending a 4-connected neighborhood to an N-connected neighborhood. We will explain how this change affects the robustness of the results in the presence of noise. We will show the results of our experiments in the forms of images and numerical results and compare them with the results of a 4-connected TRWS. We show that our model can improve the accuracy of the result and in some cases the improvement is more than 98%.

(34)

Chapter 4 Experiments and results

As explained in Chapter 3, the TRWS algorithm uses a 4-connected neighborhood (Fig. 3.1) to pass messages through a graph (or an image) and update them. This implies that the smoothness energy function in TRWS (Section 1.3), prevents sudden changes and maintains the smoothness in a small group of pixels. The presence of noise in an image can cause false fringes and sudden changes. In this thesis we propose a method to improve the robustness of the TRWS algorithm. This proposition can be used in solving the PU problem in noisy images of InSAR and enhancing the accuracy of the results. Improving the accuracy of the TRWS algorithm for solving PU problems, allows us to gain a more accurate understanding of the terrain or the surface that is being captured. Our proposition, which we will refer to as N-connected neighborhood, is to use an extended neighborhood for TRWS. In this chapter, we will provide the details of our contribution. We will show the results of our experiments and compare the results with a 4-connected TRWS.

4.1 N-connected neighborhood

As previously mentioned in Section 1.2, there are different sources of noise (e.g. orbital effect) that affect interferograms. Some of these sources can be modeled and removed. However, the wrapped phase images that need unwrapping still contain noise. This noise affects the accuracy of the results.

We use energy minimization algorithms for solving the PU problem. Finding the minimum of the energy function in the PU problem is an NP-hard problem [5]. Hence, the algorithm solves the problem by finding an approximation of the answer.

(35)

Having an approximated answer is another reason that affects the accuracy. In this thesis, our goal is to enhance the accuracy of the results of solving the PU problem for InSAR interferograms. TRWS is an energy minimization algorithm. As previously mentioned in Chapter 3, using TRWS has some advantages over similar algorithms such as BP, graph cuts, and TRW [5]. TRWS gives results that are closer to the global minimum compared to other energy minimization algorithms [4]. Hence, we focus on this algorithm.

In Chapter 2 we explained how an energy function consists of data and smoothness energy functions. Data energy function shows the consistency between the result and the observed data. Minimizing this function means minimizing the inconsistency between each pixel and the observed data. The smoothness energy function shows the smoothness between each pixel and its neighbors. Minimizing the smoothness function is equivalent to minimizing sudden changes (i.e. smooth transition from a pixel to its adjacent one).

In the TRWS algorithm, the smoothness energy measures the smoothness in a 4-connected neighborhood (Fig. 3.1) for each pixel. In interferograms, the changes in phase of a small neighborhood of pixels are very small. Hence, the changes in a small neighborhood are smooth. If large changes are observed in an interferogram, it usually means that noise is present. To reduce the effect of noise in the final unwrapped image, it is required to enhance the robustness of the algorithm that is used. Robustness means that the algorithm is resistant to noise.

One method for increasing robustness is to force smoothness in a neighborhood larger than 4 pixels. Our contribution in this thesis is to introduce a new neigh-borhood for improving the accuracy of the result. We refer to this neighneigh-borhood as N-connected neighborhood. This is because by just changing one factor we can extend the neighborhood. N-connected neighborhoods are considered in the shape of squares. Therefore, we can start with one pixel at the center and then add neighbors on each side. Thus, the neighborhoods can be 3 × 3, 5 × 5, 7 × 7, etc. In a 3 × 3 neighborhood, the number of neighbors will be 8 as shown in Fig. 4.1a. In Fig. 4.1b and 4.1c, 5 × 5 and 7 × 7 neighborhoods are shown, respectively.

Although N can be any integer, we recommend to choose N based on the density of the fringes in the image. For example if the image contains a peak of a mountain or edge of a building, N should be a small integer to ensure that the sudden changes (peak of the mountain and edge of the building) are not ignored (smoothed out). In addition, the size of the input image is important for choosing N. Choosing a large

(36)

(a) 8-connected (b) 24-connected (c) 48-connected

Figure 4.1: Examples of N-connected neighborhood

neighborhood (compared to the size of the image) will not improve the accuracy and it is only time and memory consuming.

Furthermore, the choice of N is directly related to the time complexity of the PU problem. Having a large neighborhood means processing more edges for each pixel and thus takes more time. In addition to that, having a larger neighborhood means needing more memory because we need to store longer vectors. Therefore, there is a trade-off between time/memory and the accuracy of the results.

4.2 N-connected TRWS

There are some changes that need to be made to convert a 4-connected TRWS to an N-connected one. As we mentioned in Section 4.1, the size of the neighborhood is effective for the smoothness energy function. Eq. 3.3 is the smoothness energy function for solving the PU problem in InSAR application. This function includes two summations. Each summation belongs to a neighbor that has a higher ordering than the main pixel (in Eq. 3.3 the main pixel is denoted by (i, j)). In a 4-connected neighborhood, there are two neighbors with a higher ordering than the main pixel which is half of the total number of neighbors. It is the same for N-connected neigh-borhoods. For 8, 24, and 48- connected neighborhoods, the numbers of neighbors with higher ordering than the main pixel are 4, 12, and 24, respectively. Hence, the new form of the smoothness energy is

Esmoothness= M −2 X i=0 N −1 X j=0 n/2 X c=1

ki+c,j− ki,j + Ni+c,j/2π+ M −1 X i=0 N −2 X j=0 n/2 X c=1

ki,j+c− ki,j+ Ni,j+c/2π ,

(37)

where n is the total number of neighbors. Since the data energy function (Eq. 3.4) does not depend on neighbors, it does not need any change for the N-connected algorithm.

In addition to the energy functions, we need to extend the arrays of message and shift. For any neighbor and for each possible label in a defined neighborhood, we have a message and a shift. Hence, we need to define the size of these arrays based on the same factor that was mentioned earlier in Section 4.1 for switching between different sizes of N-connected neighborhoods. For this factor, we used the number of pixels in each side of a square from the center. In Fig. 4.1, the radius of 8, 24, and 48 is 1, 2, and 3, respectively.

In Section 3.2, we described the steps of the TRWS algorithm. The second step included two parts. First one was computing ˆθs and the second part was updating

and normalizing Mst;k. In both parts, we had to repeat the computations for every

neighbor. Thus, to convert a 4-connected TRWS to an N-connected TRWS we need to extend the loops in the code that implement the second step, therefore we include all the new neighbors.

4.3 Experiments

• Data set:

To conduct our experiments, we decided to use synthetic images. This decision was mainly made because of lack of ground truth for real interferograms. Hence, we were not able to compare the accuracy of the results of different neighbor-hoods. In addition, the size of a real interferogram is another reason for this decision. The size of real interferograms is very large (for instance they can be 10000 × 10000 pixels). Thus, processing data with this size and repeating the experiments for different values of N requires a significant amount of time and memory.

To produce a synthetic data set, we utilized a generator (developed by 3vGeo-matics). The generated synthetic images are a simple proxy in lieu of real data. They are not meant to be particularly accurate to real interferograms, but they have some similar properties. For instance we can introduce noise through the generator that is similar to the noise structure we can find in real data. We can change the magnitude of noise to observe the effects on the results. By using synthetic images to test different algorithms, we can find out which method

(38)

fails in solving the PU problem. The images are generated by using a simple procedural noise methodology, in particular, perlin noise [38]. The procedural noise is used in different ways to generate the expected structure and noise that may realistically exist in some types of real data. Besides the magnitude of noise, we determine the size of images and number of layers. The latter vari-able is the initial number of images that the generator produces. Afterwards, the generator subtracts any two images from each other. The results of sub-tracting two same layers in different positions are counted the same. Hence, if the number of layers is equal to 7, then the number of produced images is

7×6

2 = 21. The produced images are similar to interferograms after PU (our

ground truth images). Further, by having a wrapper function, we produce the wrapped interferograms. Thus, we have an input data set and a ground truth data set.

• Size of images and neighborhoods:

The size of images we conducted our experiments on is 500 × 500 pixels. Besides the 4-connected TRWS, we tested the first three sizes of N-connected TRWS which are 8, 24, and 48-connected. The reason we did not consider larger neighborhoods is that for the size of our input images, a larger neighborhood does not give us a significant improvement and it is only time consuming to process.

• Numerical comparison:

It is not always easy to visually compare the results from different neighbor-hoods. We need a numerical comparison to have a more thorough conclusion. For this purpose, we calculate the mean square error (MSE). The MSE value is found by the following

M SE = 1 n n X i=1 φi− φ gt i 2 ,

where n is the number of pixels, and φ and φgt _{are unwrapped phase values of}

the result and the ground truth, respectively. • Parameters:

For each experiment, we show two examples and then repeat the experiment with 21 images for different neighborhoods to compare the average MSEs. There

(39)

are some parameters that we can change in order to produce different data sets such as noise, spatial parameter, and size of the images. Spatial parameter determines how close the fringes are to each other. We keep these parameters the same in order to produce 21 images to obtain an average MSE for each neighborhood. The following table contains the values of the aforementioned parameters in three experiments that we have conducted (more details in Section 4.4):

Exp1 Exp2 Exp3

Size 500 × 500 500 × 500 500 × 500

Noise 0 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 1.0

Spatial 6 (3vGeomatics default) 6 10

• Machine features:

For conducting the experiments, we used a personal laptop with following fea-tures:

CPU Intel Core i7-5500U CPU @ 2.40GHz × 4R

RAM 8GB OS Linux Ubuntu 16.04 L1i Cache 32K L1d Cache 32K L2 Cache 256K L3 Cache 4096K

4.4 Results and comparison

In the first experiment, the input data (wrapped images) do not contain any noise. It is clear that the data always has noise in real applications. Initially, we wanted to find out how well the TRWS algorithm works with N-connected neighborhoods in perfect conditions. Figures 4.2 and 4.3 are examples of the data with no noise. The color bar on the right side of each image shows the value of the phase at each pixel. In the wrapped image (Fig. 4.2a), the phase is between (−π, π] indicated by the color bar beside it. Figures 4.2d, 4.2e, and 4.2f show the results for the 8, 24, and 48-connected neighborhoods, respectively. Note that these figures are the same

(40)

(a) Wrapped image (b) The ground truth

(c) 4-connected result (d) 8-connected result

(e) 24-connected result (f) 48-connected result

(41)

(42)

neighborhood 4-conn 8-conn 24-conn 48-conn

MSE 0.7021 0 0 0

Table 4.1: The mean square values for Fig. 4.2 neighborhood 4-conn 8-conn 24-conn 48-conn

MSE 1.3617 0 0 0

Table 4.2: The mean square values for Fig. 4.3

as the ground truth (Fig. 4.2b). This is not the case for Fig. 4.2c which is the result of the 4-connected neighborhood. In Fig. 4.2c, note that the upper left side of the image shows inconsistency with the ground truth. Fig. 4.3 illustrates the same approaches as Fig. 4.2. Note that in Fig. 4.3c there is inconsistency in the bottom of the image compared to the ground truth. The MSE values for all N-connected neighborhoods in Figures 4.2 and 4.3 are shown in Tables 4.1 and 4.2, respectively. Note that the error of the proposed method is zero compared to the 4-connected neighborhood. To further examine using N-connected neighborhoods in the TRWS algorithm we executed the algorithm for different neighborhoods in 21 images and obtained the average MSE value for each neighborhood. These values are shown in Table 4.3. Based on Table 4.3, N-connected neighborhoods shows greater accuracy than 4-connected.

neighborhood 4-conn 8-conn 24-conn 48-conn

avg MSE 0.2541 0 0 0

Table 4.3: The average of MSE of executing the algorithm for different neighborhoods in 21 images with no noise.

The next experiment we conduct is on noisy data. As it was mentioned before in Section 4.3, the generator uses perlin noise function to produce the synthetic images. To add noise, we simply increase the magnitude of noise. The noise magnitude is the factor of the noise function in the generator. When this factor is zero it means that the noise function is not included in producing the images. By increasing the noise magnitude, we are increasing the strength of the noise function in producing the images by the generator. For the data set without any noise (e.g. Fig. 4.2) the magnitude of noise was zero. As we previously stated, the synthetic images are not supposed to demonstrate exact features of real interferograms. However, the

(43)

added noise will affect the data the same way as observed noise in real interferograms although we can not correspond the noise to a specific source of noise. Figures 4.5 and 4.6 are examples for the second experiment. Figures 4.4a and 4.4b are wrapped image and the ground truth, respectively. The effect of noise is visible in the wrapped image. The wrapped image in Fig. 4.4a is not as smooth as the no noise image in Fig. 4.2a. Figures 4.4c, 4.4d, 4.4e, and 4.4f are the results of 4, 8, 24, and 48-connected respectively. Fig. 4.5 illustrates the same approaches as Fig. 4.4. Tables 4.4 and 4.5 show the MSE values of different neighborhoods of Fig. 4.4 and 4.5, respectively.

neighborhood 4-conn 8-conn 24-conn 48-conn MSE 15.6235 0.1511 0.1020 0.0871

Table 4.4: The MSE values of Fig. 4.4

neighborhood 4-conn 8-conn 24-conn 48-conn MSE 2.7538 0.1034 0.0697 0.0634

Table 4.5: The MSE values of Fig. 4.5

Based on these tables, we can find a significant difference between the MSE values of 4-connected and the smallest N-connected TRWS (8-connected). For example in Fig. 4.4, there is a 99% improvement in 8-connected compared to 4-connected.

We decided to gradually increase the noise magnitude to observe the effect of the noise on accuracy. For each level of noise we executed the algorithm for 21 images for different neighborhoods to obtain average MSEs. The different noise magnitudes that we tested the algorithms with are 0.5, 1, 1.5, 2, 2.5, and 3. The noise magnitude in Figures 4.4 and 4.5 is 1. We did not increase the magnitude beyond 3 because the accuracy of the results is very low and the running time very large. Table 4.6, shows the average MSE results of executing the algorithm for different neighborhoods in 21 images for different levels of noise. As the results show, by increasing the noise magnitude, there is an increase in MSE value of each neighborhood (i.e. a decrease in the accuracy level). However, N-connected algorithms (4, 24, and 48) decrease the MSE value significantly. For instance, when the magnitude is 1, there is %98 improvement by comparing the results of 4 and 48-connected neighborhoods. Table 4.7 shows the standard deviation for each neighborhood and noise magnitude. Since the noise is generated by a random function, there is no guarantee that running the

(44)

(45)

(46)

neighborhood 4-conn 8-conn 24-conn 48-conn magnitude = 0.5 0.1097 0.0001 0.00003 0.00001 magnitude = 1.0 6.1768 0.0972 0.0873 0.0607 magnitude = 1.5 7.8298 0.7459 0.5932 0.5189 magnitude = 2.0 9.1679 2.042 1.663 1.396 magnitude = 2.5 10.5941 3.6047 2.9925 2.4596 magnitude = 3.0 13.3317 5.074 4.8351 3.948

Table 4.6: The average MSE of executing the algorithm for different neighborhoods in 21 images for different noise magnitudes.

neighborhood 4-conn 8-conn 24-conn 48-conn magnitude = 0.5 0.4068 0.0001 0.00008 0.00004 magnitude = 1.0 2.7004 0.0077 0.0032 0.003 magnitude = 1.5 4.9192 0.04746 0.0158 0.0142 magnitude = 2.0 7.2339 0.1615 0.14 0.635 magnitude = 2.5 9.3592 1.0661 0.9258 0.7388 magnitude = 3.0 12.1087 2.0925 1.7208 1.0871

Table 4.7: The standard deviation of executing the algorithm for different neighbor-hoods in 21 images for different noise magnitudes.

experiment for another 21 times, will result in the same average MSE obtained in Table 4.6. Hence, we can define a confidence interval which is an interval where the average MSE obtained from running the experiments for several times will fall within it. To calculate a confidence interval we need to determine a confidence level. This is a level of certainty (expressed in percentage) which indicates that the average MSE will be in the obtained confidence interval after running the experiment for an arbitrary number of times. We calculated the confidence interval for the data in Table 4.6. The confidence level is 95%. Table 4.8 shows the result of this calculation.

As previously mentioned, the density of fringes in the image is effective in choosing the size of the neighborhood. If there are sudden changes in a small neighborhood, then a smaller N can show more accurate results than a greater one. By changing the spatial factor of the generator we can produce images with different densities of fringes. Fig. 4.6a is an example of an image with a higher density of fringes.

(47)

neighborhood 4-conn 8-conn 24-conn 48-conn magnitude = 0.5 0.1097 ± 0.174 0.0001 ± 0 0.00003 ± 0 0.00001 ± 0 magnitude = 1.0 6.1768 ± 1.155 0.0972 ± 0.0033 0.0873 ± 0.0014 0.0607 ± 0.0013 magnitude = 1.5 7.8298 ± 2.104 0.7459 ± 0.0203 0.5932 ± 0.0068 0.5189 ± 0.0061 magnitude = 2.0 9.1679 ± 3.094 2.042 ± 0.0691 1.663 ± 0.0599 1.396 ± 0.0272 magnitude = 2.5 10.5941 ± 4.003 3.6047 ± 0.456 2.9925 ± 0.396 2.4596 ± 0.316 magnitude = 3.0 13.3317 ± 5.179 5.074 ± 0.895 4.8351 ± 0.736 3.948 ± 0.465 Table 4.8: The confidence interval for different neighborhoods with different noise levels.

(a) An example of an image with high density of fringes.

.

(b) An example of an image with low density of fringes.

Figure 4.6: Comparing images with different fringe density.

We do PU on 21 images similar to Fig. 4.7.a to show that in images with high density of fringes, smaller connected TRWS works better than 4-connected and larger N-connected TRWS. The noise magnitude of these 21 images is 1 and the spatial factor stays the same for all images. Table 4.9 shows the average MSEs of PU on these images. Based on Table 4.9, greater value of N, shows less accuracy and greater MSE. However, 8-connected TRWS, has better accuracy than 4-connected. These results prove that our method, even for images with sudden changes in small areas, shows better accuracy than the 4-connected TRWS. Therefore, it is important to analyze the input data before implementing the algorithm to choose the most suitable size of neighborhood.

Phase unwrapping by N-connected TRW-S algorithm for InSAR images

Contents

List of Tables

List of Figures

Introduction

1.1

Synthetic aperture radar

1.2

Interferometry

1.3

Phase unwrapping

1.4

Energy minimization

1.5

Related works

1.6

Contributions and outline

Chapter 2

Preliminaries

2.1

Message passing algorithms

2.2

Markov random fields

2.3

Maximum a posteriori estimation

2.4

Summary

Chapter 3

Sequential tree re-weighted

message passing algorithm

3.1

Introduction

3.1.1

Notations

3.2

The TRWS algorithm

3.2.1

Weak tree agreement

3.3

Using TRWS in phase unwrapping

3.4

Summary

Chapter 4

Experiments and results

4.1

N-connected neighborhood

4.2

N-connected TRWS

4.3

Experiments

4.4

Results and comparison