Geometric Methods for Low Level Image Processing ...

3 downloads 0 Views 11MB Size Report
Gösta H. Granlund. Preview. The multi-resolution approach has been ...... [44] Johansson M. The Hilbert transform. Master's thesis, Vaxjo University. [45] Weiss J.
Geometric Methods for Low Level Image Processing M Sc (Med Phys) E Ulises Moya-S´ anchez

Electrical and Computer Sciences Deparment CINVESTAV Unidad Guadalajara, Jalisco, M´exico [email protected]

Abstract

The visual system is the most widely studied sensory modality. Some physiological results claim that lines (local even signal) and edges (local odd signal) are the basic features in the processing of the mammalian visual system. Therefore, it is easy to understand that the local feature detection plays an important role in human perception and in image processing and computer vision. In this work we address the topic of image processing for computer vision applications. Image processing is a wide field, therefore, we focus on low-level and mid-level image processing, that corresponds to the first layers of a bottom-up computer vision system. The main aim of this report is to propose a geometric approach, based on atomic function filter bank and local phase, in the framework of geometric algebra. In this regard, we show the advantages to carry out computations of local phase using the quaternion atomic function and using the atomic Riesz transform, instead of the classical global-phase approach. In addition, we explain how an atomic function, formulated in the quaternionic algebra framework, can be used as a building block to perform multiple analytical operations commonly used in image processing, such as low-pass filter, derivatives, lines or edge detection, Hilbert and Riesz transform, steering filters, multi-resolution analysis and symmetry feature detection. One advantage of this approach is that filters, entities, and operations are defined in geometric algebra by using hypercomplex numbers. The experimental part shows that this is promising approach to low and mid-level image processing using the geometric algebra framework.

i

Resumen El sistema visual de los mamiferos es el sistema visual m´as ampliamente estudiado. Algunos resultados fisiol´ogicos afirman que las l´ıneas y los bordes son los elementos b´asicos en el procesamiento visual de los mam´ıferos. Por lo tanto, es f´acil de entender que la detecci´on de caracter´ısticas locales en las im´agenes juegan un papel importante en la percepci´on humana y en el procesamiento de im´agenes y visi´on computacional. En este trabajo se aborda el tema de procesamiento de im´agenes para aplicaciones de visi´on computacional. El procesamiento de im´agenes es un campo muy amplio, por lo tanto, nos enfocamos en el procesamiento de im´agenes de bajo y medio nivel , que corresponde a las primeras capas de un sistema de visi´on computacional. El objetivo principal de este trabajo es el de proponer un modelo geom´etrico, para el procesamiento de bajo y medio nivel usando el ´algebra de geom´etrica, basado en una funci´on at´omica y la fase local. En este sentido, se muestran las ventajas de llevar a cabo los c´alculos de la fase local usando la funci´on at´omica con las transformadas de Hilbert y de Riesz. Esta funci´on se puede utilizar como un bloque de construcci´on para llevar a cabo m´ ultiples operaciones anal´ıticas utilizadas en el procesamiento de im´agenes, tales como filtro de paso bajas, derivadas, extracci´on de l´ıneas o detecci´on de bordes, an´alisis multi-resoluci´on y detecci´on de elementos caracter´ısticos de simetr´ıa. Una ventaja de este enfoque es que los filtros, las entidades, y las operaciones se definen en ´algebra geom´etrica mediante el uso de n´ umeros hipercomplejos. La parte experimental muestra que este modelo es prometedor para el procesamiento de bajo y medio nivel de im´agenes.

Contents

1 Introduction

1

1.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Main Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3

Main Parts of our Approach . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.4

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5

Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.6

Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.7

Outline of this Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2 Atomic Function Theory

9

2.1

Atomic Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.2

Atomic Functions Properties . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2.1

Atomic Windows and Wavelets . . . . . . . . . . . . . . . . . . . . .

11

Mother Atomic Function up(x) . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.3.1

Remarkable Properties of up(x) . . . . . . . . . . . . . . . . . . . . .

12

2.3.2

Multidimensional up Functions . . . . . . . . . . . . . . . . . . . . .

15

Radial Atomic Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.4.1

Radial Atomic Function P lop . . . . . . . . . . . . . . . . . . . . . .

17

2.5

The derivative of up, dup(x) . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2.6

2D derivatives of up(x) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

2.3

2.4

iii

iv

Contents

3 Geometric Algebra and Quaternion Algebra

21

3.1

Why do we use Geometric Algebra? . . . . . . . . . . . . . . . . . . . . . . .

22

3.2

Geometric Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

3.2.1

Conformal Model of Geometric Algebra . . . . . . . . . . . . . . . . .

25

3.2.2

Geometric Transforms in Exponential Form . . . . . . . . . . . . . .

25

3.2.3

2D Symmetries and Point Groups and GA . . . . . . . . . . . . . . .

27

3.3

Quaternion Algebra

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.3.1

Complex Fourier Transform . . . . . . . . . . . . . . . . . . . . . . .

30

3.3.2

Quaternion Fourier Transform . . . . . . . . . . . . . . . . . . . . . .

31

3.3.3

Quaternion Atomic Function Qup(x, y) . . . . . . . . . . . . . . . . .

33

4 Phase Information

35

4.1

Phase Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.2

Fourier Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.3

Instantaneous Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.3.1

Hilbert Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.3.2

Analytic Signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

4.4

Hilbert Transform Using AF . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

4.5

Local Phase Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.6

2D Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

4.6.1

Quaternion Analytic Function . . . . . . . . . . . . . . . . . . . . . .

45

Monogenic Signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

4.7.1

Riesz Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

4.8

Riesz Transform Using AF . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

4.9

Symmetry and Local Phase . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

4.7

5 Multi-resolution and Steerable Filters

55

5.1

Multi-resolution

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

5.2

Quaternion Atomic Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

5.3

Steerable Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

Contents 5.4

v Riesz Atomic Multiscale . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 Experimental Analysis

58 61

6.1

Quaternion Atomic Function Qup(x, y) . . . . . . . . . . . . . . . . . . . . .

61

6.2

Derivatives of up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

6.3

Phase Computation

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

6.3.1

Quaternionic Phase Computation . . . . . . . . . . . . . . . . . . . .

64

6.3.2

Local Phase and Orientation with Atomic Riesz Transform . . . . . .

67

6.4

Multi-resolution and Steerable Filters

. . . . . . . . . . . . . . . . . . . . .

70

6.4.1

Multi-scale and Steerable Filters Using Qup(x, y) . . . . . . . . . . .

70

6.4.2

Atomic Riesz Multi-scale . . . . . . . . . . . . . . . . . . . . . . . . .

70

6.5

Symmetry Feature Extraction . . . . . . . . . . . . . . . . . . . . . . . . . .

71

6.6

Linking to High Levels in Image Processing . . . . . . . . . . . . . . . . . . .

74

7 Conclusions

79

A Apendix 1: Geometric Algebra

83

A.1 Geometric Algebra: An Outline . . . . . . . . . . . . . . . . . . . . . . . . .

83

A.2 Conformal Geometric Algebra . . . . . . . . . . . . . . . . . . . . . . . . . .

84

A.2.1 The point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

A.2.2 Spheres and Planes . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

A.2.3 Circles and Lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85

A.3 Rigid Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85

Bibliography

89

Chapter

1

Introduction It is not surprising at all that almost all blind mathematicians are geometers. The spatial intuition that sighted people have is based on the image of the world that is projected on their retinas; thus it is a two (and not three) dimensional image that is analysed in the brain of a sighted person. A blind persons spatial intuition on the other hand, is primarily the result of tile and operational experience. It is also deeper - in the literal as well as the metaphorical sense. [...] Prof. Alexei Sossinsky

1.1

Introduction

The visual system is the most widely studied sensory modality [1, 2]. The human vision system has been studied in different ways, such as light perception, motion and pattern recognition [3, 2]. Some physiological results claim that lines (local even signal) and edges (local odd signal) are the basic features in a mammalian visual system [4, 5]. Therefore, it is easy to understand that the local features play an important role in human perception, image processing and computer vision implementation. In this work we address the topic of image processing for computer vision applications. Image processing is an extensive field, therefore we focus on low-level and mid-level image processing, that corresponds to the first layers of a bottom-up computer vision system [1]. The computer vision layers are listed below. • High level (Image measurements in, image description out): Image understanding. • Middle level (Image in, data out): Image analysis, descriptors, measurements. • Low level image processing (image in, image out): feature detection, lines, edges, etc. 1

2

Chapter 1. Introduction

Computer vision is considered as the way in which computers can see, by which, we mean, what objects the system can detect, recognize and locate and, in general, understand the scenes [1]. In visual scenes, one uses primitives to model the world and combine them, with operations or certain transformations such as rigid body motions, or projective imaging. Each of those elements of geometry and their relationships need to be modeled consistently in proper algebra for geometry [6, 7]. According to [6, 7, 8] the linear algebra falls in two ways: in its representational power due to the low level of its primitives and in its lack of universal geometric operations and transformation. As a consequence, present-day software for geometry contains many tricks and ad hoc extensions, such as overloading the use of vectors, or invoking quaternions [6, 7, 8]. The various formalisms which are used require increasingly more involved and specialized conversions, degrading performance and flexibility, leading to implementation errors. To hide these problems, one may encapsulate algebra elements in objects, but this causes an extra layer of software. One of the main objectives of this work is to propose a new approach for low-level image processing and analysis capable of extracting low and mid-level geometric information, in which it is possible to link the computer vision levels using geometric algebra GA. In this regard, a new geometric approach was proposed. That approach is based in 3 parts: • Atomic functions as building blocks in which it is possible to compute and to perform multiple analytical operations commonly used in image processing as low-pass, bandpass, derivatives, Hilbert and Riesz transform and multiscale. • Geometric algebra and quaternion algebra which include complex and hypercomplex numbers. These entities allow us to carry out low-level image computations as local phase and local orientation. Consequently, this information can be then used to define a geometric object in an image. • A multi-scale approach to extract different features, sizes or orientations. Complex and hypercomplex numbers which are contained in GA play a substantial role in image processing, especially to filter the images and to obtain local features. In the case of 1D signals, the analytic signal is typically computed using a complex number and its Hilbert transform. This complex representation allows us to compute the phase and magnitude of the signal in space domain. For high-dimension signals, the Riesz transform and partial Hilbert transforms have been used as extensions of the analytic signal to compute the local phase and local orientation. However, Riesz and Hilbert transforms are global transforms. In order to select a local feature it is necessary to select a window of the size of the locality. In this work we use the atomic filter to compute the local phase and orientation. The phase information can be used to detect low level geometric characteristics such as

1.2. Main Objective

3

lines or edges using the local decomposition of the structure according to its symmetries [2, 3, 9]. Moreover, the local phase information can be invariant to illumination and to rotations [10, 11, 12]. In order to show some applications in low level image processing, we find lines and oriented edges using the quaternion atomic function its derivatives and the local phase computed by the quaternion analytic signal and the monogenic signal. The link between the low level image processing and higher levels is established by 2D symmetries. In this work we define some geometrics objects in terms of the centroid, corner and main axis in the plane of the image as well as the normal vector in order to formulate GA constrains. Hence, a practical example using a geometric neural Network is shown.

1.2

Main Objective

One of the main objectives of this work is to propose a new approach for low-level image processing and analysis capable of extracting low and mid-level geometric information in which it is possible to link the computer vision levels using geometric algebra GA. The ideal characteristics of this approach are:

1. It should be preferably linear due to linear systems or models, which are typically simpler in terms of features and properties. 2. It must be capable of enhancing images, as well as representing and extracting low-level geometric information, such as edges and lines. 3. It should be simple and compact.

Each of those ideal characteristics and their relationships were straightforward modeled using Geometric Algebra. We chose e a geometrical model, because geometric algebra allows us to model geometric primitives, operations and transformations such as rigid body motions, or projective geometry in the same mathematical framework. Besides, our approach resembles to the mammalian visual system; we use the same input lines and edges in a compact representation as in a mammalian visual system. Additionally, in this work we propose a new way to use an atomic function up(x) as a kernel to compute the Hilbert and Riesz transforms.

4

Chapter 1. Introduction

1.3

Main Parts of our Approach

In order to achieve the above mentioned objectives, we have proposed a new approach. This Approach is divided into three parts, as presented in Figure 1.1

Figure 1.1: Main parts of our proposed approach.

1. Compact support filters (windows) based on solutions of differential functional equations, called Atomic Functions (AF ). By definition, an Atomic Function AF is a compactly supported and infinitely differentiable solution of functional differential equations. We use the most important AF , the up(x) function [13] because it has some useful properties, such as good local properties in space-domain, the fact that a derivative of any order can be expressed easily, and the approximation to other functions is relatively simple [13]. Additionally, the Hilbert and the Riesz transform was developed in terms of its derivatives. Moreover, radial atomic functions were used to generate a radial band-pass in a multi-scale approach. 2. The second part is the GA and the quaternion framework. We use that to compute to the local phase and orientation information in the images [10, 11]. The phase information shows the geometric structure of the image, as the local decomposition of the image (2D signal) according to its symmetries [10] and it has also been shown by Oppenheim et al [14] that the phase information is more important for the perceptual impression of an image than the magnitude. Additionally we use the conformal geometry in order to show how it is possible to define geometric objects by using local features. Finally, some transformation was used in a geometric neural network in order to encode the pattern angle. The intrinsic dimension is an important characteristic of the image because its possible to define the feature neighborhood (see Figure 1.2). 3. The third part of our approach is a filter bank in a multi-resolution approach and steerable filters (see Fig 1.3). In this work we propose a multi-scale processing with steerable filters based on a quaternionic atomic function Qup(x) and on an Atomic Riesz for multi-scale to get different features, sizes and different orientation information preserving an illumination invariance and rotation invariance.

1.4. Overview

5 0iD

1iD

2iD

Figure 1.2: Intrinsic dimension. From left to right 0iD, 1iD and 2iD.

Figure 1.3: A multi-scale and steerable filter bank of up(x, y).

1.4

Overview

Our approach will be applied by the geometric algebra group at CINVESTAV-Guadalajara. In this group, a mobile robotic and a medical robotic system are being developed. The Figure 1.4 shows the work scenario of the medical robotic system with stereo-endoscopic cameras and a US probe. This system can be helpful in assisting the surgeon when performing laparoscopic procedures.

1.5

Related Work

This thesis is primary based on the quaternion analytic signal proposed by T. B¨ ulow [11] and the monogenic signal developed by M. Felsberg [10]. Both approaches propose a generalization of the analytic signal into 2D dimensions to compute the local phase information for local feature extraction. These approaches were based on hypercomplex and GA formulations and they claim that complex algebra is not enough to formulate the 2D analytic signal [10]. The most important differences between these two approaches, is that the quaternionic analytic functions use the Hilbert transform and the partial Hilbert transform, whereas the monogenic signal use the Riesz transform and its invariant to rotations. The practical implementation of these approaches use the quaternionic Gabor and Log-Gabor filters in a multi-scale approach respectively. It has been shown that the monogenic signal is the most general approach for performing the 2D extension of the analytic signal [10, 15, 16].

6

Chapter 1. Introduction

Figure 1.4: The main application of our approach will be in stereo endoscopic images. The main difference between these two approaches and our approach, is that, we propose a another filter with good local property, and our approach is useful for both, low and mid-level image processing. Instead of using a Gabor or a Log-Gabor filter, we use a compact support window (with good local properties) the AF up(x, y), which is easy to derivative with dyadic operations and obtain the Hilbert and Riesz transform. A local compact support window allows an accurate local phase estimation with lower computational load than the Gabor or Log-Gabor; due to the uncertainty of this filters [10]. Moreover, our approach links the phase information and the geometric feature extraction with image analysis and symmetry, which finally helps in the geometric objects definitions taking into account its transformations.

1.6

Contribution

Our main contributions are listed below. 1. Novel filters and operations were proposed in this work. The atomic function up(x) was formulated in the quaternionic algebra framework. Additionally, we use these windows as building blocks to perform multiple analytical operations commonly used in image processing such as low-pass filter, derivatives, local phase processing (by computing the

1.7. Outline of this Thesis

7

local Hilbert and Riesz transforms), steering quaternionic filters and multi-resolution analysis. A radial atomic function was used to obtain a radial bandpass in order to implement an atomic monogenic signal. 2. In this work we are presenting the way in which the phase and its symmetry response can be used to obtain geometric features and we present one example on how the geometric features can be connected with high level processing by using the local symmetry.

1.7

Outline of this Thesis

The following list presents the outline as follows: • The main subject of chapter 2 is devoted to the atomic functions. Not only by highlighting the main characteristics of the AF , but also by introducing some advantages of these functions. Hence, we focus on the mother atomic function up(x) and its radial representation. This chapter also introduces the basis for the atomic functions theory and its properties. • In chapter 3 we answer two main questions: Why is geometric algebra used? and What is geometric algebra?. Additionally, in this chapter the reader will find basic definitions for quaternion algebra, the quaternion Fourier transform and the quaternion atomic function. The last part of this chapter is devoted to the basis of geometric algebra and its applications to model geometric objects. • In Chapter 4 the local phase concept is explained. As we will see, the phase information offers a natural model to disentangle the information according to their even and odd symmetries. Additionally, in this chapter, we will introduce the Fourier phase, the instantaneous phase, and the local phase. Each of this phases are explained and their main differences are highlighted. As we mentioned in previous sections, we are interested in the local phase, therefore, we devote the main part of this chapter to it. We explain some approaches to compute the local phase such as the quaternion analytic signal and the monogenic signal. In this regard, some properties of the Hilbert and the Riesz transform are highlighted. • Chapter 5 is devoted to generate a multi-resolution and steerable filters bank for both quaternion atomic functions and Riesz atomic transform. In this chapter we show how a radial compact support window can be performed by using a difference of radial atomic function.. • In Chapter 6 some applications are developed using the theory of the chapters 2, 3, 4 and 5. • Finally, the conclusions are presented in chapter 7.

Chapter

2

Atomic Function Theory If, in some cataclysm, all of scientific knowledge were to be destroyed, and only one sentence passed on to the next generation of creatures, what statement would contain the most information in the fewest words? I believe it is the atomic hypothesis that All things are made of atoms-little particles that move around in perpetual motion, attracting each other when they are a little distance apart, but repelling upon being squeezed into one another. In that one sentence, you will see, there is an enormous amount of information about the world, if just a little imagination and thinking are applied. Richard Feynman

Preview This chapter introduces the basis for the atomic functions theory and its properties. The atomic functions (AF ) are one of three main parts of our approach. The use of AF s has been important contributions to mathematical analysis, signal processing, pattern recognition, image processing, numerical methods, and so forth. To put it simply, we use them as windows or filters. By definition, an Atomic Functions AF are compactly supported and infinitely differentiable solutions of functional differential equations. We especially focus on the mother atomic function up(x) and its derivatives dup(x) In the AF class, the function up(x) is the simplest and at the same time, the most useful primitive function to generate other kinds of atomic functions. In this chapter the main properties of the up(x) as well as the multidimensional atomic functions and its derivatives are presented. 9

10

Chapter 2. Atomic Function Theory

2.1

Atomic Functions

The Atomic Functions1 (AF ) were first developed in the 1970s, jointly by V. L and V. A Rvachev [17, 18]. Since then, the use of AF s has generated very important contributions in mathematical analysis, signal processing, pattern recognition, image processing, numerical methods, and so forth [19, 17, 18, 20, 9, 21, 22, 23].

Definition 2.1.1 By definition, the atomic functions are compactly supported and infinitely differentiable solutions of functional differential equations (see Eq 2.1) with a shifted argument [13]. The atomic functions have the property in which the nth-order derivative of an AF is a linear combination of the same function (or its derivatives) with a shift in the argument i.e. Lf (x) = λ

m X

c(k)f (ax − b(k)),

(2.1)

k=1

where, |a| > 1, b, c, mλ ∈ N L = constant coefficients.

dn dxn

n−1

d + a1 dx n−1 + ... + an is a linear differential operator with

A functional is widely used in classical mechanics to compute the Lagragian. A functional 0 differential equation is a differential equation in which the derivative y (t) of an unknown function y has a value at t which is related to y by a function of other function, i.e. a function that takes functions and returns functions. Some functions (Atomic Functions) which can solve the Eq 2.1 are listed below: 1 up(x) = 2π 1 f upN (x) = 2π

Z



iνx

e −∞

Z



k=1 iνx

e −∞

∞ Y sin(ν2−k )

ν2−k

dν,

∞ sin(ν2) N Y sin(ν2−k ) ( ) dν, ν2 ν2−k k=1

∞ X 1 (a − 1)πk ha (x) = (a − 1)( + sinc[ ] cos[(a − 1)πkx]), 2π k=1 ak Z ∞ ∞ Y ¡ sin ν(n + 1)−k ¢n 1 Ξ(x) = eiνx dν. −k 2π −∞ ν(n + 1) k=1 1

Also known as R-functions

(2.2) (2.3) (2.4) (2.5)

2.2. Atomic Functions Properties

2.2

11

Atomic Functions Properties

In this work we mostly use the up(x) function (Eq 2.2). Many useful properties of the AF have been reported in [21, 22, 13, 19, 17, 18] and in this work we take advantage of some properties such as: 1. There are explicit equations for the values of the moments and Fourier transform; hence, several operations can be computed easily. 2. The AF s are compactly supported functions in the space domain. Therefore, we can obtain good local characteristics in space domain. 3. Any polynomial can be represented by means of their translations thus, we can approximate any polynomial function. 4. Since derivatives of any order can be represented in terms of simple shifts, we can easily represent any derivative operator. 5. AF s are infinitely differentiable (C ∞ ), as a result, the AF and their Fourier transforms are rapidly decreasing functions. Their Fourier transforms decrease on the real axis faster than any power function.

2.2.1

Atomic Windows and Wavelets

The AF s has been used as filters or windows [13, 18]. The AF s windows were compared with classic ones by means of parameters such as: the equivalent noise bandwidth, the 50% overlapping region correlation, the parasitic modulation amplitude, the maximum conversion losses (in decibels), the maximum side lobe level (in decibels), the asymptotic decay rate of the side lobes (in decibels per octave), the window width at the six-decibel level, the coherent gain, etc. All atomic windows exceed classic ones in terms of the asymptotic decay rate [13, 18]. Additionally, real and complex wavelets based on AF theory have been proposed by [13, 18, 24, 21, 22]. AF s are well localized due to their compactness and their Fourier transforms decrease on the real axis, faster than any power function [13]. In fact, the Fourier transform of an AF corresponds to infinite products of periodic functions; therefore the atomic functions themselves generate a wavelet system. The idea of Atomic Wavelets (AW ) consists of translations and dilations of some mother functions as other wavelets approach. The AF wavelets have been compared with other wavelets, such as Daubechies, and a wavelet based on AF has shown a similar or a better performance [21]. As we will see in chapter 3 a new hypercomplex Atomic function was proposed the Qup based on quaternion Fourier transform. Moreover, a quaternion atomic multi-scale and steerable filters are presented in chapter 5.

12

2.3

Chapter 2. Atomic Function Theory

Mother Atomic Function up(x)

In the AF class, the function up(x) is the simplest and at the same time, the most useful primitive function to generate other kinds of atomic functions [13]. The Atomic Function up(x) is generated by infinite product of sinc functions. The function up(x) has the following representation in terms of the Fourier transform [25, 17, 18]:

1 up(x) = 2π 1 = 2π

Z Z

∞ ∞ Y

sin(ν2−k ) iνx e dν, ν2−k −∞ k=1 ∞

iνx up(ν)e ˆ dν

(2.6) (2.7)

−∞

Figure 2.1 shows five convolutions in space domain of rectangle functions. As can be seen from this figure, only with 5 convolutions is enough for obtain a up(x) function. In the following Figure 2.2 the plots of up(x) and its Fourier transform up(ν) ˆ are shown. According to Figure 2.2 the up(x) has good local property due to its compact support in space domain. 1.2

1

0.8

0.6

0.4

0.2

0

−0.2

0

20

40

60

80

100

120

140

160

180

200

x

Figure 2.1: Convolutions of square functions in space domain, only with 5 convolutions is enough for obtain a up(x) function.

2.3.1

Remarkable Properties of up(x)

From previous works, the main properties of the up(x) are listed below [25, 17, 18, 13]: • The function up(x) has compact support [−1, 1] in space domain.

2.3. Mother Atomic Function up(x)

13 F(up(x)) 1

0.8

0.8

|F(up(x))|

I(a.u)

up(x) 1

0.6

0.4

0.2

0.6

0.4

0.2

0 −1

−0.5

0

0.5

0 −200

1

Space domian (x)

−100

0

100

200

Frequency domain (ν)

Figure 2.2: Left. Atomic function up(x). Right. Fourier Transform up(ν). ˆ • up(x) is even, up(−x) = up(x) R1 • Its maximum value is up(0) = 1 and −1 up(x)dx = 1. • n-order derivatives are defined by:

(2.8)

n

d

(n)

up(x) = 2

n(n+1)/2

2 X

δk up(2n x + 2n + 1 − 2k),

(2.9)

k=1

where δ2k = −δk , δ2k−1 = δk , δ1 = 1. • Translations with smaller steps yield polynomials (xn ) of any degree, i.e. ∞ X

ck up(x − (k2n )) ≡ xn , ck , x ∈ R.

(2.10)

k=−∞

Figure 2.3 introduces an example of how the up(x) can be used to represent any polynomial by means of translation.The upper curve is approximate by up functions. • We can express the up(x) functions in terms of shifts up(x) = 1 − up(1 − x). • Moments of even order of up(x). Z an =

1

xn up(x)dx,

(2.11)

(2.12)

−1 n

a2n where a0 = 1 and a2n+1 = 0.

(2n)! X a2n−2k = 2n , 2 − 1 k=1 (2n − 2k)!(2k + 1)!

(2.13)

14

Chapter 2. Atomic Function Theory

Aproximation using up(x) 0.2

I(a.u)

0.15

0.1

0.05

0

0

20

40

60

80

100

Space domain (x)

Figure 2.3: Approximation by means of translations of up(x). • Moments of odd order of up(x) Z

1

bn =

xn up(x)dx,

(2.14)

0

b2n+1

¶ µ n+1 X 2k 1 a2n−2k+2 , = 2n + 2k (n + 1)22n+1 k=1

(2.15)

where b2n = a2n . • The Fourier expansion of up(x) in even harmonics is up(x) = 0.5 +

∞ X

up(πk)cos(π(2k ˆ − 1)x)

(2.16)

k=−∞

and its possible to express the functional differential equation (and its solutions) such as [26]: a2 (ha (ax + 1) − ha (ax − 1)) , 2 n m 1 X Y sin(kπ(a − 1)2−i ) ha (x, m, n) = + 2 k=1 i=1 kπ(a − 1)2−i cos(πkx), ha =

(2.17) (2.18)

1 1 where a > 1, and its support is [− a−1 , a−1 ]. In the following Figure 2.4 we use Equation 2.18 with i = 5, a = 2 and different k = 1, 2, 3.

2.3. Mother Atomic Function up(x)

15

1.2 k=1 k=2 k=3

1 0.8 0.6 0.4 0.2 0 −0.2

0

10

20

30

40

50

Figure 2.4: Equation 2.18 with i = 5, a = 2 and different k = 1, 2, 3.

2.3.2

Multidimensional up Functions

Digital processing of multidimensional signals is one of the most promising applications of AF s [13]. A natural extension of the up(x) to the case of many variables is based on the usual outer product of 1D up(x) [23]. The simplest case of 2D atomic function is computed by a simple outer product as follows up(x, y) = up(x)up(y).

(2.19)

According to [13] as in the 1D signals, due to the infinite smoothness of the AF s, the side lobes of these windows are characterized by infinite asymptotic decay rate. Figure 2.5 depicts the square function and its Fourier transform (Sinc(x)), which is equivalent to using k = 1 in Eq 2.2. It is shown in Figure 2.6 the Eq 2.19 using k = 10, 100 in Eq 2.2. It can be seen in the Figures 2.5 and 2.6 that they don’t have a radial symmetry in the space domain.

Figure 2.5: Square function in space domain (left) and in frequency domain (right). k = 10.

16

Chapter 2. Atomic Function Theory

Figure 2.6: up(x, y) in space domain (left) and frequency domain using k = 100 see Eq 2.2 Figure 2.7 shows 2D up(x, y) functions in space and frequency domain. In this Figure, the 2D up(x, y) shows a compact support in space domain.

Figure 2.7: Top: Two views of up(x, y). Bottom: Fourier Transform up(ν, ˆ υ).

2.4

Radial Atomic Functions

Other 2D atomic functions were defined in [13] such as ³p ´ Radial(up) = up x2 + y 2 , √ √ x+y 3 x−y 3 Hexagonal(up) = up(x)up( )up( ). 2 2

(2.20) (2.21)

2.4. Radial Atomic Functions

17

Unfortunately, the tensor product and rotation don’t allow us to synthesize two-dimensional filters with arbitrary support domains especially for concave domains, for example; cross or star shape (there are other Atomic functions in order to approximate the compact support) [13]. In this work we use the Fourier series approximation (see Eq 2.18) in order to build an atomic radial function i.e.

where a > 1, r =

p

a2 ha = (ha (ar + 1) − ha (ar − 1)) , 2 n m 1 X Y sin(kπ(a − 1)2−i ) ha (r, m, n) = + 2 k=1 i=1 kπ(a − 1)2−i cos(πkr)

(2.22) (2.23)

1 1 , a−1 ]. x2 + y 2 and its support is [− a−1

p Figure 2.8: up( x2 + y 2 ) in space domain and in frequency domain using Fourier series approximation. As we will see in the chapter 4 it is important to use a radial function in order to compute the Riesz transform.

2.4.1

Radial Atomic Function P lop

Another atomic radial function, P lop(r), was introduced by [27]. The P lop(r) function was defined as a radial infinite differentiable function with compact support (see Fig 2.9), in frequency domain i.e.

P lop(ν, ˆ υ) =

∞ X ∞ Y h=0 k=0

[−(u2 + ν 2 )]k 32k(h+1) [(k + 1)!]2

and its solution of the following functional-differential equation [27] Z 2 ∇ P lop(x, y) = λ P lop[3(x − ξ1 ), 3(y − ξ2 )]ds + µP lop(3x, 3y), ∂S

(2.24)

(2.25)

18

Chapter 2. Atomic Function Theory 2

2

∂ ∂ 2 2 5 where ∇2 = ∆ = e1 ∂x 2 +e2 ∂y 2 , ξ1 +ξ2 = 4/9, µ = −4πλ/3 and λ = 3 /4π. Figure 2.9 depicts the P lop in Fourier domain in comparison with the radial up function. As it can be seen in the last figure a radial atomic function offers an invariant rotation response in frequency and space domain.

√ Figure 2.9: Left, up( ˆ ν 2 + υ 2 ); right, P lop(ν, ˆ υ).

2.5

The derivative of up, dup(x)

There are some mask operators like Sobel, Prewitt, Kirsh, which are used to extract edges in images. A common drawback of these operations is that it is impossible to ensure in their performance the required characteristics in their performance over a wide range of the working band of the processed signals and retune to adapt to the signal parameters [19]. This means that adaptation of the differential operator to the behavior of the input signal by broadening or narrowing its band is desirable, in order to assure maximum signal-to-noise ratio [19]. This problem is reduced to the synthesis of infinitely differentiable finite functions with small wide-band which are used for constructing the weighting windows [25, 18, 19]. One of the most effective solution is obtained with the help of the atomic functions [19]. The AF can be used in two ways: the construction of a window in frequency region to obtain the required improvement in properties of the pulse characteristic, or direct synthesis based on Eq 2.1 [19]. Therefore the function up(x) satisfies the Eq. 2.1,in this way:

dup(x) = 2up(2x + 1) − 2up(2x − 1),

(2.26)

2.6. 2D derivatives of up(x)

19

and if we compute the Fourier Transform of Eq. 2.26 we obtain ¡ ¢ iνF [up(2x)] = eiν − e−iν F [up(2x)]), F (dup(x)) = 2i sin(ν)F (up(2x)).

(2.27) (2.28)

By differentiating Eq. 2.1 term by term we obtain n

(n)

d

up(x) = 2

n(n+1)/2

2 X

δk up(2n x + 2n + 1 − 2k),

(2.29)

k=1

where δ2k = −δk , δ2k−1 = δk , δ1 = 1. The function dup(x) provides a good window in the spatial frequency regions because, the side lobes have been completely eliminated [19].

2.6

2D derivatives of up(x)

Similar to Eq 2.19 we can get a 2D expression of each derivative:

dup(x, y)x = dup(x)up(y) dup(x, y)y = up(x)dup(y) dup(x, y)x,y = dup(x)dup(y)

(2.30) (2.31) (2.32)

Figure 2.10 presents two graphics of dup(x, y)y , dup(x, y)x,y in space domain. It is possible to implement an approximation of derivatives with the difference of two functions. The figure 2.11 shows a difference of Gauss and a difference of radial rup. In this work we use a difference of radial atomic function in order to implement an approximation of the monogenic signal.

20

Chapter 2. Atomic Function Theory

Figure 2.10: Left dup(x, y)x,y , rigth dup(x, y)y .

0.25 0.2

0.08

0.15

0.06

0.1 0.04 0.05 0.02 0 0 −0.05 −0.02 100 −0.04 0

20

40

50 60

80

100

0

−0.1 0 50 100

0

20

40

60

80

Figure 2.11: Left, difference of radial up and on the right, difference of Gauss.

100

Chapter

3

Geometric Algebra and Quaternion Algebra As long as algebra and geometry have been separated, their progress has been slow and their uses limited; but when these two sciences have been united, they have lent each mutual forces, and have marched together towards perfection. J. L. Lagrange

Preview The main aim of this chapter is to provide an answer for a fundamental question: Why do we use Geometric Algebra? In order to answer this question, we claim that GA offers a complete mathematical approach. It has primitives to model the geometric elements, and it combines operations or transformations with those elements, such as rotations, reflections or translations in order to apply rigid body motions, or projective geometry. In fact, we use the GA framework to extend the up(x, y) in an even-sub algebra, which is isomorphic to quaternion algebra; that allows us to obtain new filters with different symmetries and permits us to compute the local phases of the images and link the low-level image processing with the upper layers of computer vision systems. Additionally, in this chapter the reader will find the basis of geometric algebra and its applications to model geometric objects. The last part of this chapter is devoted to basic definition of quaternion algebra, quaternion Fourier transform and quaternion atomic function. 21

22

Chapter 3. Geometric Algebra and Quaternion Algebra

3.1

Why do we use Geometric Algebra?

Geometric algebra GA is one area of Clifford algebra; its name was introduced1 by David Hestenes in order to emphasize the geometric interpretation [29]. GA allows describing and solving geometrical problems in a compact manner [28, 7, 29, 30, 8]. By using GA, the modeling (descriptions and operations) of geometric objects can be done simply, and at the same time better [31, 32]. But that simplicity is not always totally clear at the beginning, as it performs basic operations in a different manner and notation [7, 31, 32, 8]. In order to clarify this point, we enumerate the main advantages of GA that become relevant in our work as follows 1. GA is a compact and efficient language of geometric entities [33, 28]. In others words, we can combine geometric transformations (rotations, reflections, etc) with the construction and intersection of geometric entities in a single framework. 2. The GA allows direct computing with subspaces, which simplifies many geometrical operations, as subspaces form an algebra with well-defined products that have direct geometric significance. Some researchers claim that geometric algebra is a better choice than (vector calculus) a language that always reduces everything to vectors (which are just 1D subspaces) [34, 8]. 3. The classical approaches make use of a mixture of two kinds of entities; complex valued entities and matrices. In that formulation it is harder to switch between these two representations [7, 29, 8]. 4. According to B¨ ulow and Felsberg, it has been shown that the complex algebra is not enough to formulate and derive a 2D phase (generalization of 2D analytic signal). 5. In this work, we use a quaternion algebra H, which is isomorphic to the even geometric sub-algebra of the 3D euclidian geometric algebra (see section 3.2). The quaternions is a popular method to represent rotations in a 3D space instead of matrices [29, 30], because their components are easier to interpret (rotation axis direction and angle). 6. The rotations can have some numerical advantages. For instance, two rotations can be represented by quaternions (16 multiplications and 12 additions) or by 3 × 3 matrices (27 multiplications and 18 additions) [29]. 7. The scaling-rotations operations can be written as only one operation using quaternions [10]. 1

Suggested by Clifford himself [28]

3.2. Geometric Algebra

23

8. Many research groups have used the GA in order to represent and solve a variety of problems in different areas including the group of CINVESTAV Guadalajara. These groups use the GA framework because the complexity of modeling and the computational overhead can be substantially reduced. As mentioned above, many applications of GA have been performed in computer vision applications. However, these applications use conventional methods in order to obtain the input information (points, lines, corners, edges). Therefore, we propose a new approach in low-level image processing and analysis capable of extracting low and mid-level geometric information in which it is possible to link all computer vision levels using GA.

3.2

Geometric Algebra

The GA is a linear space of dimension 2n , n = p + q + r, and (p, q, r, ∈ N ). The notation Gp,q,r represents the p,q,r vector basis that squares to +1, -1, 0 respectively. We use the R3 (G3,0,0 ) which has the following elements (blades) basis:  

 

1 , e1 , e2 , e3 , e23 , e31 , e12 , e123 ≡ I3  |{z} | {z } | {z } | {z } scalar

vectors

bivectors

(3.1)

trivectors

In Eq. (3.1) the scalar is a blade of zero degree, vectors are of first degree, bivectors are blades of second degree and pseudo scalar I3 is a blade of third degree. The Clifford product (of two vectors) is defined as the sum of the inner product and the wedge product: ab = a · b + a ∧ b

(3.2)

The wedge product is anticonmutative a ∧ b = −b ∧ a. The scalar product of the basis vectors is given by the Kroneckers delta (δ) ei · ej = δi,j .

(3.3)

The outer product or wedge product is given by (if i 6= j)

ei ∧ ej = eij .

(3.4)

The psedoscalar I3 is used to convert a geometric entity from its standard form to its dual form and vice versa.

24

Chapter 3. Geometric Algebra and Quaternion Algebra

1 e1 e2 e3 e23 e31 e12 I3

1 1 e1 e2 e3 e23 e31 e12 I3

e1 e1 1 −e12 e31 I3 e3 −e2 e23

e2 e2 e12 1 −e23 −e3 I3 e1 e31

e3 e3 −e31 e23 1 e2 −e1 I3 e12

e23 e23 I3 e3 −e2 -1 e12 −e31 −e1

e31 e31 −e3 I3 e1 −e12 −1 e23 −e2

e12 e12 e2 −e1 I3 e31 −e23 −1 −e3

I3 I3 e23 e31 e12 −e1 −e2 −e3 −1

Table 3.1: Geometric product of the basis.

I3 = e1 e2 e3 = e123

(3.5)

All products between the basis are presented in the Table 3.1 By definition a multivector is a sum of different blades. Let be Ar y Bs multivectors of grade r and s respectively, then the Clifford product of two multivectors can be expressed by Ar Bs = hAr Bs ir+s + hAr Bs ir+s−2 + ... + hAr Bs i|r−s| ,

(3.6)

Ar · Bs = hAr Bs i|r−s| ,

(3.7)

Ar ∧ Bs = hAr Bs ir+s .

(3.8)

where by definition

The GA of R3 can be used to represent the group of rotations and scale transformations by spinors [10]. The main idea of the spinors is to generalize a rotation using a unit complex number. These representation combine the most compact representation, with the numerically fastest evaluation of the quaternion product. Let R be a multivector with unit magnitude [10, 35] R = a + be23 + ce23 + de12 , R = e

−θ n 2

,

(3.9) (3.10)

where a2 + b2 + c2 + d2 = 1, n is a unit bivector and θ represents the rotation angle. The spinors form an even subalgebra which is isomorphic to quaternion algebra H. In the next section we focus on the quaternion algebra because it can be expresses as subalgebra of the 3D euclidean geometric algebra and at this moment it is enough to describe the 2D phases of the image.

3.2. Geometric Algebra

3.2.1

25

Conformal Model of Geometric Algebra

The conformal model is commonly applied in computer vision and robotics [36, 7, 8], since the conformal model of the GA extends the 3D space in order to express lines spheres, planes or circles as elements of the algebra and include, translations, dilations or inversion. According to [11], and our previous results, we claim that the conformal model of GA is the appropriate approach in which it is possible to model all the levels of image processing. In the following list the elements of Conformal Geometric Algebra CGA in each levels of image processing are shown: 1. The low-level image processing operations, can be computed by using the bivector part of CGA, which is isomorphic to H, in order to detect the points of interest, lines or edges, at different orientations or sizes [11, 10, 20, 9, 37]. In fact, the quaternion analytic signal and the monogenic signal are used to compute the local phase (see chapter 4). 2. In order to face the mid-level image processing, the even and odd symmetry relations in the phase and orientation can be used, which is also encoded in the bivector part. 3. The high-level of image processing is connected by the local symmetry. The symmetry concept is given by two main parts: the object definition and the transformation of the object and many geometric objects and transforms were defined in the CGA [7, 36, 6, 8]. As is shown in Table 3.2 there are some geometric entities in CGA. In this case the input (points, directions, etc) for those entities were taken from the previous levels. Additionally in the Table 3.3 a few geometric transforms in CGA are shown. Inner Product Null Space (IPNS) and Outer Product Null Space (OPNS) are the type of representations that the geometric entities belong to.

3.2.2

Geometric Transforms in Exponential Form

In this section we introduce the use of bivector in order to represent 3D rotations. A rotor e = 1 and it is defined as follows R is an even element of G3 in which RR µ ¶ θ R = exp − l , (3.11) 2 + The components of this rotor R ∈ G3 are equal to the rotor in G4,1 . In this case a unit bivector is used l ∈ hG3 i2 ⊆ G4,1 , which represents a normalized line as a rotation axis, with a rotation angle θ. The equation (3.11) can be rewritten using the Euler representation as

26

Chapter 3. Geometric Algebra and Quaternion Algebra

Entity Sphere Point Plane

Line

Circle

Par of points

Representation Grade 1 2 2 s = p + 2 (p − ρ )e + e0 1 1 2 x = x + 2 x e + e0 1 P = NIE − de 1 N = (a − b) ∧ (a − c) d = (a ∧ b ∧ c)IE L = P1 ∧ P2 2 = rIE + eMIE r=a−b M=a∧b z = s1 ∧ s1 2 Pz = z · e, Lz ∗ = z ∧ e z2 pz = Pz ∨ Lz , ρ = (e∧z) 2 P P = s1 ∧ s2 ∧ s3 3

Representation Dual Grade ∗ s =a∧b∧c∧d 4 1 2 ∗ x = (−Ex − 2 x e + e0 )IE 4 P∗ = e ∧ a ∧ b ∧ c 4 L∗ = e ∧ a ∧ b

3

z∗ = a ∧ b ∧ c

3

P P ∗ = a ∧ b, X ∗ = e ∧ x

2

Table 3.2: List of conformal entities in IPNS and OPNS. Transform Inversion Reflection Translation Rotation

g(x) ∈ Rn ρ2 x−p

Versor in Gn+1,1 s = p − 12 ρ2 e

+p

σ(x) ³

x−p ρ

−nxn + 2xδ

P = n + δe

1

x−t

T = 1 + 12 te

1

R = R + e(c × R)

1

R(x − c)R−1 + c

´2

Table 3.3: Conformal transforms and its representation in the versor form. µ

¶ θ R = exp − l 2 µ ¶ µ ¶ θ θ = cos − sen l 2 2

(3.12)

and the equation can be defined the equation (3.12) in terms of bivector basis as R = rc − rs l = rc − rs1 e23 − rs2 e31 − rs3 e12

(3.13)

e = rc + rs l = rc + rs1 e23 + rs2 e31 + rs3 e12 , R

(3.14)

e and for the R

3.2. Geometric Algebra

27

where rc and rl ∈ R. The rotation of an geometric entity can be carry out as follows e P 0 = RP R,

(3.15)

where P is a plane and P 0 is the rotated plane. Moreover, the product of two rotors is a new rotor e = (R2 R1 )P (R f1 R f2 ), P 0 = RP R (3.16) The traslation of a geometric entity regarding t ∈ hG3 i1 is made using T ∈ G4,1 , µ T = exp T = (1 +

et 2

¶ ,

et ), 2

(3.17) (3.18)

The rotation of a geometric object is defined by e P 0 = TP T,

(3.19)

where P is a plane and P 0 is the traslation of the plane. If you apply a rotor and a traslation a motor is obtained (M) M = TR, (3.20) as example, the motor under P can be described as f = TRP R e T, e P 0 = MP M

3.2.3

(3.21)

2D Symmetries and Point Groups and GA

As we mentioned before, in order to express the symmetry, we need an object definition and a transform definition [38]. The 2D symmetries can be defined in 2 groups [39]: 1. The reflection (Cn ) cyclic n-fold rotations ( 2π ) and the reflection P (see Table 3.3). n 2. Dn dihedral n-fold reflections generated by a Cn element. In this work we define some geometrics objects in terms of the centroid Cen, corner Cor and main axis M a in the plane of the image P .

28

Chapter 3. Geometric Algebra and Quaternion Algebra

P ∗ = e ∧ a ∧ b ∧ c, 1 Cen = cen + cen2 e + e0 , 2 1 Cori = cor + cor2 e + e0 , 2 L∗ = e ∧ a ∧ b.

(3.22) (3.23) (3.24) (3.25)

The euclidian elements, cen, cor, etc. of the conformal points Cen, Cor are extracted by the phase computation, as reader will see in the chapter 6. Whereas the symmetric operations defined in GA are given by the equations 3.15, 3.15, 3.17, 3.18 and 3.20. By using centroids as rotation reference and the corners and its rotations or reflections it is possible to define a regular polygon as in the work [39]. In order to know how close a corner is extracted we use the following expressions fi = 0, Cori+1 ∧ Mi Cori M q fi |, Dif f erenceof corners = 2|ci+1 · Mi ci M

(3.26) (3.27)

as example, the equation of an equilateral triangle is obtained by the centroid and a corner. Then, the rotation in this case is ∆Θ = 120. It is easy to compute the equations of the following polygons [39]

Figure 3.1: Geometric polygons.

3.3

Quaternion Algebra

T. B¨ ulow and M. Felsberg [11, 10], claims that, the complex algebra is not enough to formulate the n−dimensional generalization of the analytic signal. Therefore, they used the quaternion + + algebra H or the G3,0,0 to develop their works. The sub-algebra G3,0,0 (bivector basis) is isomorphic to H which is an associative non-commutative four-dimensional algebra, it consists of one real scalar and three imaginary elements, q = a + bi + cj + dk q = a + be23 + ce13 + de21 ,

(3.28) (3.29)

3.3. Quaternion Algebra

29

where a, b, c, d ² R and i → e23 , j → e13 , k → e21 . The units i, j obey the relations i2 = j 2 = ijk = −1, ij = k.

(3.30)

Quaternion algebra is geometrically inspired and the imaginary components can be described in terms of the basis of R3 space [35]. The real part of q is noted by Re(q) = a and the pure part is P u(q) = bi + cj + dk [35]. The quaternions can be used to represent rotations (as bivector) in R3 and R4 and translations (as vector) in R3 [35]. One of the main differences between GA and quaternion is that, GA is able to naturally distinguish objects and operations, and in GA the rotations can be applied for a different 3D space [10, 7, 8]. Another important characteristic of the H is the phase concept. There are different representations of the polar form of q. In this work we use the following representation [11]: q = |q|eiφ ekψ ejθ

(3.31)



where |q| = q q¯, where q¯ is the conjugate of q = a − bi − cj − dk and the angles (φ, θ, ψ) represent three quaternionic phases [11]. The Eq 3.31 is based on a 3D rotation (R = Rx (2φ)Rz (2ψ)Ry (2θ)) which is an isomorphism of a unit sphere S 3 to the special orthogonal group in 3D SO(3) [11, 35]. The rotation matrix M can be factorized into three rotations R = Rx (2φ), Rz (2ψ), Ry (2θ) in the coordinate axes as follows. M(q) = M(q1 )M(q2 )M(q3 ) q1 = eiφ , q2 = ejθ , q3 = ekψ

(3.32) (3.33)

 a2 + b2 − c2 − d2 2(bc − ad) 2(bd + ac)  2(bc + ad) a2 − b2 + c2 − d2 2(cd − ab) M(q) =  2 2 2 2 2(bd − ac) 2(cd + ab) a −b −c +d

(3.34)



   R=  

cos(2ψ) cos(2θ) − sin(2ψ) cos(2ψ) sin(2θ) cos(2φ) sin(2ψ) cos(2θ) cos(2φ) cos(2ψ) cos(2φ) sin(2ψ) sin(2θ) + sin(2φ) sin(2θ) − sin(2φ) cos(2θ) sin(2φ) sin(2ψ) cos(2θ) sin(2φ) cos(2ψ) sin(2φ) sin(2ψ) sin(2θ) − cos(2φ) sin(2θ) + cos(2φ) cos(2θ)

     

(3.35)

The quaternion algebra contains 3 complex subfields with orthogonal imaginary units [11]. The phases are expressed by the following rules [11]:

ψ=−

arcsin(2(bc − ad)) 2

(3.36)

30

Chapter 3. Geometric Algebra and Quaternion Algebra

• If ψ ∈] − π4 , π4 [, then φ =

¯ )) argi (q Tj (q 2

and θ =

• If ψ = ± π4 , then select either φ=0 and θ =

¯ )q ) argj (Ti (q . 2

¯ )q ) argj (Tk (q 2

or θ=0 and φ =

¯ )) argi (q Tk (q . 2

• If eiφ ekψ ejθ = −q and φ ≥0, then φ → φ − π. • If eiφ ekψ ejθ = −q and φ 0;  1 0 if ν = 0; sign(ν) =  −1 if ν < 0.

(4.7)

4.3. Instantaneous Phase

39

The π/2 phase interpretation of the Hilbert transform is easy to see in the Eq 4.7 via the imaginary value i π −i = e−i 2 if ν > 0; π i = ei 2 if ν < 0;

Hilbert Transform Properties The most important properties of the Hilbert Transform are listed below [11, 45, 43]. 1. Linearity. The Hilbert transform of f (x) = c1 f1 (x) + c2 f2 (x) if the Hilbert transform of f1 (x) y f2 (x) is given by fH (f (x)) = c1 fH (f (x1 )) + c2 fH (f (x2 )).

(4.8)

2. |fH (f (0))| = 0, thats means the direct current DC component is equal to zero. 2 3. fH [fH (f (x))] = fH (f (x)) = −f (x).

4. An interesting property occurs when multiple Hilbert transform are used, we can get the inverse Hilbert transform 3 fH (f (x))fH (f (x)) = f (x), −1 3 fH (f (x)) = fH (f (x)).

(4.9) (4.10)

5. The Hilbert transform and the original function are orthogonal over the infinite interval. 6. If f (x) ∈ R, the fH (f (x)) ∈ R. 7. The Hilbert transform is anti-symmetric. 8. The Hilbert transform of the convolution is given by fH [f (x) ∗ g(x)] = fH [f (x)] ∗ g(x) = fH [f (x)] ∗ [g(x)].

(4.11)

9. The Hilbert transform of the derivative of the function is equivalent to the derivate of the Hilbert transform of a function d 0 f (x)fH fH (f (x)). ← → dx

(4.12)

40

4.3.2

Chapter 4. Phase Information

Analytic Signal

For 1D signals, the phase computation is carried out using of the analytical signal which is 1 based on the Hilbert transform (fH (x) = πx ) and is given by fA (f (x)) = f (x) + ifH (x), fA (f (x)) = |A|eiθ ,

(4.13) (4.14)

p where |A| = f (x)2 + fH (x)2 and θ = arctan( ffH(x) ) allow us to extract the magnitude and (x) phase independently. As a example of 1D the analytical signal, the figure 4.2 shows the original signal sin(x) and its Hilbert transform; the phase θ and the plot of the real part real(fA (f (x))) = sin(x) versus the imaginary part imag(fA (f (x))) = fH (sin(x)). The circle in figure 4.2 represents the angle in the complex plane.

Amp

2

sin(x) \F_H(sin_x)

1

0

−1 −15

−10

−5

0

5

10

15

Angle

π fractions

4

phase

2 0 −2 −4 −15

−10

−5

0

5

10

15

Angle 2

Angle(rad)

Imag

1 0 −1 −2 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Real

Figure 4.2: From top to down. The original signal sin(x) (dash-line), and its Hilbert transform. The instantaneous phase and the relation between the real (original signal) and the imaginary part (the Hilbert transform of the original signal). The figure 4.3 highlights the relation of some parts of the signal with the phase. On the top of the figure, the original signal and the odd part (left) and the even part (right) may be

4.3. Instantaneous Phase

41

seen. The down-plots show a complex plane, a real part (x − axis) and an imaginary part (y − axis). An edge (odd symmetry) with negative slope has a phase range (0, π), whereas a positive slope edge has (π, 2π). The even part of the signal is related with the line, if the line has positive values it will have a phase range of ( 3π , π ) meanwhile if the lines have negative 4 2 values, the phase range is ( π2 , 3π ). 4 1

sin(x)

sin(x)

Amp

Amp

1

0

−1

0

4

0

−1

8

0

Angle

7

1.5

Imag

Imag

1.5

0

−1.5 −1

3.5

Angle

0

Real

1

0

−1.5 −1

0

1

Real

Figure 4.3: On the top of the figure the original signal and the odd part (left) and the even part (right) may be seen. The down-plots show a complex plane, a real part (x − axis) and an imaginary part (y − axis) which are related to the even part and the odd part of the signal. In image processing, the instantaneous phase components express different image primitives such as lines or edges [2]. The Figure 4.4 compare the global and the instantaneous phase (using the Hilbert transform) of a image. The Figure 4.4 highlights the main difference between the global phase and the instantaneous phase of a image. First, the global phase is more sensitive to the change of brightness than the instantaneous phase. Second, the instantaneous phase response is related to the odd and even responses;, this is in relation to the lines and edges.

Figure 4.4: From left to right. Original image, global magnitude, global phase, instantaneous magnitude and instantaneous phase.

42

4.4

Chapter 4. Phase Information

Hilbert Transform Using AF

Theorem 4.4.1 The Hilbert transform and the partials Hilbert transforms, can be expressed in terms of the up and its derivative as follows, ¶ 1 gHi f (x, y) = f (x, y) ∗ dup(x, y)x ∗ − log(|x|) , π µ ¶ 1 gHj f (x, y) = f (x, y) ∗ dup(x, y)y ∗ − log(|y|) , π µ ¶ 1 gHk f (x, y) = f (x, y) ∗ dup(x, y)xy ∗ 2 log(|x|) log(|y|) . π µ

(4.15) (4.16) (4.17)

Prof. The Hilbert transform in terms of the convolution is defined as

fHi (x) = f (x) ∗

1 , πx

(4.18)

using the association and the distribution properties of the convolution operation f ∗ (g ∗ h) = (f ∗ g) ∗ h ∇ (f ∗ g) = ∇f ∗ g = f ∗ ∇g,

(4.19) (4.20)

∂ ∂ + e2 ∂y allows us rewrite the Hilbert transform in terms of where f , g ∈