Edge Detection Methods Based on Modified Differential Phase Congruency of Monogenic Signal

Yan Yang School of Mathematics (Zhuhai), Sun Yat-Sen University,
Zhuhai, 519080, China.
Email: mathyy@sina.com Kit Ian Kou Department of Mathematics, University of Macau, Macao (Via Hong Kong).
Email: kikou@umac.mo Cuiming Zou Department of Mathematics, University of Macau, Macao (Via Hong Kong).
Email: zoucuiming2006@163.com

Abstract

Monogenic signal is regarded as a generalization of analytic signal from one dimensional to higher dimensional space, which has been received considerable attention in the literature. It is defined by an original signal with its isotropic Hilbert transform (the combination of Riesz transform). Like the analytic signal, monogenic signal can be written in the polar form. Then it provides the signal features representation, such as the local attenuation and the local phase vector. The aim of the paper is twofold: first, to analyze the relationship between the local phase vector and the local attenuation in the higher dimensional spaces. Secondly, a study on image edge detection using modified differential phase congruency is presented. Comparisons with competing methods on real-world images consistently show the superiority of the proposed methods.

keywords:

Hilbert transform; Phase space, Poisson operator
AMS Mathematical Subject Classification: 44A15; 70G10; 35105

^†^†journal: Multidimensional Systems and Signal Processing

1 Introduction

In the scale-space literature, there are a lot of papers discussing Gaussian scale-space as the only linear scale-space [3, 14, 15]. The Gaussian scale space is obtained as the solution of the heat equation. In [11], M. Felsberg and G. Sommer proposed a new linear scale-space which is generated by the Poisson kernel, it is the so-called Poisson scale-space in two-dimensional (2D) spaces.

The Poisson scale-space is obtained by Poisson filtering (the convolution of the original signal and the Poisson kernel). The harmonic conjugate (the convolution of the original signal and the conjugate Poisson kernel) yields the corresponding figure flow at all scales. The Poisson scale-space and its corresponding figure flow form the monogenic scale-space [11]. In mathematics, monogenic scale-space is the Hardy space in the upper half complex plane. The boundary value of a monogenic function in the upper half space is the monogenic signal. The monogenic scale signal gives deeper insight to low level image processing.

Monogenic signal is regarded as a generalization of analytic signal from one dimensional space to higher dimensional case, which is first studied by M. Felsberg and G. Sommer in 2001 [10]. It is defined by an original signal with its Riesz transform in higher dimensions. Under certain assumptions, monogenic function can be representation in the polar form and then it provides the signal features, such as the local attenuation and the local phase vector [10, 21, 20, 8]. In [21], we first defined the scalar-valued phase derivative (local frequency) of a multivariate signal in higher dimensions. Then we studied the applications in signal processing [20, 19, 24, 25, 26].

Monogenic signals at any scale $s>0$ form monogenic scale-space. The representation of monogenic scale-space is just a monogenic function in the upper half-space. Therefore, considering the monogenic scale space with scale $s$ instead of monogenic signals, provides us more analysis tools. In the monogenic scale-space, the important features in image processing, such as local phase-vector, and local attenuation (the log of local amplitude) involving through scale $s$ are given in [11]. The relationship between the local attenuation and the local phase in the intrinsically 1D cases are derived in [11]. However, the problem is open if the signal is not intrinsically 1D signal.

The contributions of this paper are summarized as follows.

1.

We give the solution of the problem: if the higher dimensional signal is not intrinsically 1D signal, we derive the relationship between the local phase-vector and and local attenuation.
2.

We proposed the local attenuation (LA) method for edge detection operator. We establish the theoretical and experiment results on the newly methods.
3.

We proposed the modified differential phase congruency (MDPC) method for the edge detection operator. We establish the theoretical and experiment results on the newly methods.
4.

We show that in higher dimensional space, the instantaneous frequency in higher dimensional spaces defined by is equal to the minus of the scale derivative of the local attenuation.

We show that the zero points of the differential phase congruency is not equal to the extrema of the local attenuation. The nonzero extra term

-{\rm Vec}\left[\left(\underline{D}\frac{\underline{v}}{|\underline{v}|}\right)\frac{\underline{v}}{|\underline{v}|}\right]\sin^{2}\theta+(\sin\theta\cos\theta-\theta){\partial\over\partial s}\left({\underline{v}\over|\underline{x}|}\right)

appears in high dimensional cases.

The rest of this paper is organized as follows. In order to make it self-contained, Section 2 gives a brief introduction to some general definitions and basic properties of Hardy space, analytic signal, Clifford algebra, monogenic signal and monogenic scale space. In Section 3 we derive the relationship between the local phase-vector and and local attenuation. Various edge detection methods are provided in Section 4. Finally, experiment results are drawn in Section 5.

2 Preliminaries

In the present section, we begin by reviewing some definitions and basic properties of analytic signal and Hardy space [5, 12, 13].

2.1 Analytic Signal and Hardy Space

Definition 2.1 (Analytic Signal).

For a square integrable real-valued function $f$ , the complex-valued signal $f_{A}$ whose imaginary part is the Hilbert transform of its real part is called the analytic signal. That is,

f_{A}(x):=f(x)+{\bf i}\mathcal{H}(f)(x),

where $\mathcal{H}(f)(x)$ is the Hilbert transform (HT) of $f$ defined by

\displaystyle\mathcal{H}[f](x):=\frac{1}{\pi}{\rm p.v.}\int_{-\infty}^{+\infty}\frac{f(s)}{x-s}ds=\frac{1}{\pi}\lim\limits_{\varepsilon\to 0^{+}}\int_{\varepsilon\leq|x-s|}\frac{f(s)}{x-s}ds,

provided this integral exists as a principal value ( ${\rm p.v.}$ means the Cauchy principle value).

Due to its definition, the real $u$ and imaginary parts $v$ of analytic signal $f_{A}=u+{\bf i}v$ form the Hilbert transform pairs

\displaystyle\mathcal{H}[u]=v.

(1)

To proceed the properties of analytic signal, we introduce the notion of Hardy space [12, 13], we will notice that the class of analytic signals is the class of boundary values of Hardy space functions.

Definition 2.2 (Hardy Space).

The Hardy space $H^{2}({\bf C}^{+})$ is the class of analytic functions $f$ on the upper half complex plane ${\bf C}^{+}:=\{x+{\bf i}y\,|\,x\in{\bf R},y>0\}$ which satisfies the growth condition

\left(\int_{-\infty}^{\infty}|f(x+{\bf i}y)|^{2}dx\right)^{1/2}<\infty,

for all $y>0$ .

Important properties of Hardy functions are given by Titchmarsh’s Theorem [16].

Theorem 2.1 (Titchmarsh’s Theorem).

Let $g:=u+{\bf i}v\in H^{2}({\bf C}^{+})$ . Then the following two assertions are equivalent:

1.

The Hardy function $g$ has no negative-frequency components. That is,

$g(z)={1\over 2\pi}\int_{0}^{\infty}G(\omega)e^{{\bf i}\omega z}d\omega,$

where $G(\omega):=\int_{{\bf R}}g(x)e^{-{\bf i}\omega x}dx$ is the Fourier transform of $g$ .
2.

The real and imaginary parts verify the formulas:

$u(x+{\bf i}y)=u\ast P_{y}(x)=\int_{{\bf R}}P_{y}(x-t)u(t)dt,$

and

$v(x+{\bf i}y)=v\ast Q_{y}(x)=\int_{{\bf R}}Q_{y}(x-t)v(t)dt,$

for all $y>0$ , where $P_{y}(x)={1\over\pi}{y\over x^{2}+y^{2}}$ and $Q_{y}(x)={1\over\pi}{x\over x^{2}+y^{2}}$ are the Poisson and conjugate Poisson kernel in ${\bf C}^{+}$ .

In this way, an analytic signal $f_{A}=f+{\bf i}\mathcal{H}[f]$ represents the boundary values of Hardy function $u+{\bf i}v$ in the upper half plane ${\bf C}^{+}$ [12]. That is,

f(x)=\lim_{y\rightarrow 0}u(x+{\bf i}y)

and

\mathcal{H}[f](x)=\lim_{y\rightarrow 0}v(x+{\bf i}y).

Starting from this concept we are going to study the higher dimensional generalization on Clifford algebra.

2.2 Clifford Algebra

For all what follows we will work in ${\it Cl}_{0,m}$ , the real Clifford algebra. Most of the basic knowledge and notations in relation to Clifford algebra are referred to [2, 7]. Let ${\bf e}_{1},...,{\bf e}_{m}$ be basic elements satisfying ${\bf e}_{i}{\bf e}_{j}+{\bf e}_{j}{\bf e}_{i}=-2\delta_{ij}$ , where $\delta_{ij}=1$ if $i=j,$ and $\delta_{ij}=0$ otherwise, $i,j=1,2,\cdots,m.$ The Clifford algebra ${\it Cl}_{0,m}$ is the associative algebra over the real field ${\bf R}$ . A general element in ${\it Cl}_{0,m}$ , therefore, is of the form $x=\sum_{S}x_{S}{\bf e}_{S},x_{s}\in{{\bf R}}$ , where ${\bf e}_{S}={\bf e}_{i_{1}}{\bf e}_{i_{2}}\cdots{\bf e}_{i_{l}},$ and $S$ runs over all the ordered subsets of $\{1,2,\cdots,m\},$ namely $S=\{1\leq i_{1}<i_{2}<\cdots<i_{l}\leq m\},\quad 1\leq l\leq m.$

Let

\ {\bf R}^{m}=\{\underline{x}\;|\;\underline{x}=x_{1}{\bf e}_{1}+\cdots+x_{m}{\bf e}_{m},x_{j}\in{\bf R},j=1,2,\cdots,m\}

be identical with the usual Euclidean space and an element in ${\bf R}^{m}$ is called a vector. Moreover, let

{\bf R}_{1}^{m}=\{x\;|\;x=x_{0}+\underline{x},x_{0}\in{\bf R},\underline{x}\in{\bf R}^{m}\}

be the para-vector space and an element in ${\bf R}_{1}^{m}$ is called a para-vector. The multiplication of two para-vectors $x_{0}+\underline{x}=\sum_{j=0}^{m}x_{j}{\bf e}_{j}$ and $y_{0}+\underline{y}=\sum_{j=0}^{m}y_{j}{\bf e}_{j}$ is given by $(x_{0}+\underline{x})(y_{0}+\underline{y})=(x_{0}y_{0}+\underline{x}\cdot\underline{y})+(x_{0}\underline{y}+y_{0}\underline{x})+(\underline{x}\wedge\underline{y})$ with $\underline{x}\cdot\underline{y}=-<\underline{x},\underline{y}>=-\sum_{j=1}^{m}x_{j}y_{j}$ and $\underline{x}\wedge\underline{y}=\sum_{i<j}{\bf e}_{ij}(x_{i}y_{j}-x_{j}y_{i}).$ Clearly, we have

\displaystyle\underline{x}\cdot\underline{y}=\frac{1}{2}(\underline{x}\underline{y}+\underline{y}\underline{x})

(2)

and

\underline{x}\wedge\underline{y}=\frac{1}{2}(\underline{x}\underline{y}-\underline{y}\underline{x}).

There are three parts of $(x_{0}+\underline{x})(y_{0}+\underline{y})$ . We denote them as follows

1.

the scalar part: $x_{0}y_{0}+\underline{x}\cdot\underline{y}={\rm{Sc}}[(x_{0}+\underline{x})(y_{0}+\underline{y})]$ ,
2.

the vector part : $x_{0}\underline{y}+y_{0}\underline{x}={\rm{Vec}}[(x_{0}+\underline{x})(y_{0}+\underline{y})]$ ,
3.

the bi-vector part : $\underline{x}\wedge\underline{y}={\rm{Bi}}[(x_{0}+\underline{x})(y_{0}+\underline{y})]$ .

In particular, we have $\underline{x}^{2}=-<\underline{x},\underline{x}>=-|\underline{x}|^{2}=-\sum_{j=1}^{m}x_{j}^{2},\mbox{ for }\underline{x}\in{\bf R}^{m}$ .

The conjugation and reversion of ${\bf e}_{S}$ are defined by $\overline{\bf e}_{S}=\overline{\bf e}_{il}\cdots\overline{\bf e}_{i1}$ and ${\overline{\bf e}_{j}=-{\bf e}_{j}}$ , respectively. Therefore, the Clifford conjugate of a para-vector $x_{0}+\underline{x}$ is $\overline{x_{0}+\underline{x}}=x_{0}-\underline{x}$ . It is easy to verify that $0\not=x_{0}+\underline{x}\in{\bf R}_{1}^{m}$ implies

(x_{0}+\underline{x})^{-1}:=\frac{\overline{x_{0}+\underline{x}}}{|x_{0}+\underline{x}|^{2}}.

The open ball with center $0$ and radius $r$ in ${\bf R}_{1}^{m}$ is denoted by $B(0,r)$ and the unit sphere in ${\bf R}_{1}^{m}$ is denoted by $S^{m}$ .

The natural inner product between $x$ and $y$ in ${\it Cl}_{0,m}$ is defined by $<x,y>:=\sum_{S}x_{S}\overline{y_{S}},$ where $x:=\sum_{S}x_{S}{\bf e}_{S},x_{S}\in{\bf R}$ and $y:=\sum_{S}y_{S}{\bf e}_{S},y_{S}\in{\bf R}$ . The norm associated with this inner product is defined by $|x|=<x,x>^{1\over 2}=(\sum_{S}|x_{S}|^{2})^{{1\over 2}}.$

Let $\Omega$ be an open subset of ${\bf R}_{1}^{m}$ with a piecewise smooth boundary. We say that function $f$ defined on $\Omega$ such that $f(x_{0}+\underline{x})=\sum_{S}f_{S}(x_{0}+\underline{x}){\bf e}_{S}$ is a Clifford-valued function or, briefly, a ${\it Cl}_{0,m}$ -valued function, where $f_{S}$ are real-valued functions defined in $\Omega$ .

A possibility to generalize complex analytic is offered by following the Riemann approach, which is introduced by means of the generalized Cauchy-Riemann operator $\frac{\partial}{\partial{x_{0}}}+\underline{D}$ , where $\underline{D}={\partial\over\partial x_{1}}{\bf e}_{1}+\cdots+{\partial\over\partial x_{m}}{\bf e}_{m}$ is the Dirac operator. Nullsolutions to this operator provide us with the class of the so-called monogenic functions.

Definition 2.3.

(Monogenic Function) An ${\it Cl}_{0,m}$ -valued function $f$ is called left (resp. right) monogenic in $\Omega$ if $\left(\frac{\partial}{\partial{x_{0}}}+\underline{D}\right)f=0$ (resp. $f\left(\frac{\partial}{\partial{x_{0}}}+\underline{D}\right)=0$ ) in $\Omega$ .

In the following, let

E(x_{0}+\underline{x})={\overline{x_{0}+\underline{x}}\over|x_{0}+\underline{x}|^{m+1}}

be the Cauchy kernel defined in ${\bf R}_{1}^{m}\setminus\{0\}$ . It is easy to verify that $E(x_{0}+\underline{x})$ is a monogenic function in ${\bf R}_{1}^{m}\setminus\{0\}$ [2, 7].

Remark 2.1.

1.

For a ${\it Cl}_{0,m}$ -valued function defined on an open subset of ${\bf R}^{m}$ , we apply the Dirac operator $\underline{D}$ for the monogenic function.
2.

Throughout the paper, and unless otherwise stated, we only use left ${\it Cl}_{0,m}$ -valued monogenic functions that, for simplicity, we call monoginic. Nevertheless, all results accomplished to left ${\it Cl}_{0,m}$ -valued monogenic functions can be easily adapted to right ${\it Cl}_{0,m}$ -valued monogenic functions.

We further introduce the right linear Hilbert space of integrable and square integrable ${\it Cl}_{0,m}$ -valued functions in $\Omega\subset{\bf R}^{m}$ that we denote by $L^{1}{(\Omega,{\it Cl}_{0,m})}$ and $L^{2}{(\Omega,{\it Cl}_{0,m})}$ , respectively. If $f\in L^{1}{({\bf R}^{m},{\it Cl}_{0,m})}$ , the Fourier transform of $f$ is defined by

\displaystyle\hat{f}(\underline{\xi})=\int_{{\bf R}^{m}}e^{-{\bf i}<\underline{x},\underline{\xi}>}f(\underline{x})d\underline{x},

(3)

if in addition, $\hat{f}\in L^{1}{({\bf R}^{m},{\it Cl}_{0,m})}$ , then function $f$ can be recovered by the inverse Fourier transform

f(\underline{x})=\frac{1}{(2\pi)^{m}}\int_{{\bf R}^{m}}e^{{\bf i}<\underline{x},\underline{\xi}>}\hat{f}(\underline{\xi})d\underline{\xi}.

The well-known Plancherel Theorem for Fourier transform of $f$ and $g\in L^{2}({\bf R}^{m},{\it Cl}_{0,m})$ holds

\displaystyle\int_{{\bf R}^{m}}f(\underline{x})g(\underline{x})d\underline{x}=\frac{1}{(2\pi)^{m}}\int_{{\bf R}^{m}}\hat{f}(\underline{\xi})\overline{\hat{g}(\underline{\xi})}d\underline{\xi}.

In a recent paper [10], the authors defined the notion of the monogenic signal. It is regarded as an extension of the notion of the analytic signal to multidimensional signals.

2.3 Monogenic Signal and Monogenic Scale Space

The monogenic signal was defined by an original signal and its ”isotropic Hilbert transform” in the higher dimensional spaces (a combination of the Riesz transforms).

Definition 2.4 (Monogenic Signal).

For $f\in L^{2}({\bf R}^{m},{\it Cl}_{0,m})$ , the monogenic signal $f_{M}\in L^{2}({\bf R}^{m},{\it Cl}_{0,m})$ is defined by

f_{M}(\underline{x}):=f(\underline{x})+H[f](\underline{x}),

where $H[f]$ is the isotropic Hilbert transform of $f$ defined by

$\displaystyle H[f](\underline{x})$	$\displaystyle:=$	$\displaystyle p.v.\frac{1}{\omega_{m}}\int_{{\bf R}^{m}}\frac{\overline{\underline{x}-\underline{t}}}{\|\underline{x}-\underline{t}\|^{m+1}}f(\underline{t})d\underline{t}$
	$\displaystyle=$	$\displaystyle\lim_{\epsilon\rightarrow 0^{+}}\frac{1}{\omega_{m}}\int_{\|\underline{x}-\underline{t}\|>\epsilon}\frac{\overline{\underline{x}-\underline{t}}}{\|\underline{x}-\underline{t}\|^{m+1}}f(\underline{t})d\underline{t}$
	$\displaystyle=$	$\displaystyle-\sum_{j=1}^{m}R_{j}(f)(\underline{x}){\bf e}_{j}.$

Furthermore,

R_{j}(f)(\underline{x}):=\lim_{\epsilon\rightarrow 0^{+}}\frac{1}{\omega_{m}}\int_{|\underline{x}-\underline{t}|>\epsilon}\frac{x_{j}-t_{j}}{|\underline{x}-\underline{t}|^{m+1}}f(\underline{t})d\underline{t},

is the jth-Reisz transform of $f$ and $\omega_{m}=\frac{2\pi^{\frac{m+1}{2}}}{\Gamma(\frac{m+1}{2})}$ is the area of the unit sphere $S^{m}$ in ${\bf R}^{m}_{1}$ .

Remark 2.2.

If $f(\underline{x})$ is real-valued, then by Definition 2.4, $H[f](\underline{x})$ is vector-valued.

Let us now generalize the notion of Hardy space to multidimensional space.

Definition 2.5 (Monogenic Scale Space).

The monogenic scale space $M^{2}({\bf R}_{1}^{m,+})$ is the class of monogenic functions $f^{+}(\underline{x},s)$ defined on half space

{\bf R}_{1}^{m,+}=\{x\;|\;x=(\underline{x},s),\underline{x}\in{\bf R}^{m},s>0\},

which satisfies the growth condition

\left(\int_{{\bf R}^{m}}|f^{+}(\underline{x},s)|^{2}d\underline{x}\right)^{1/2}<\infty,

for all scale $s>0$ .

Like in the complex case, a monogenic signal is the boundary value of the monogenic scale function in the half space ${\bf R}_{1}^{m,+}$ [4]. Some basic properties of the Monogenic scale space $M^{2}({\bf R}_{1}^{m,+})$ in the half space are summarized as follows. For the proof of Theorem 2.2 we refer the reader to [4] and [17].

Theorem 2.2.

Suppose $f^{+}(\underline{x},s):=u(\underline{x},s)+v(\underline{x},s)\in M^{2}({\bf R}_{1}^{m,+})$ . Then the following two assertions are equivalent:

The inverse Fourier transform of $f^{+}$ vanishes for $s<0$ . That is, the ${\it Cl}_{0,m}$ -valued function $f^{+}(\underline{x},s)$ has the form

\displaystyle f^{+}(\underline{x},s)=\frac{1}{(2\pi)^{m}}\int_{{\bf R}^{m}}e^{+}(s+\underline{x},\underline{t})\hat{f}(\underline{t})d\underline{t}

in the half space $s>0$ , where

e^{+}(s+\underline{x},\underline{t})=e^{-s|\underline{t}|}e^{{\bf i}<\underline{x},\underline{t}>}(1+{\bf i}\frac{\underline{t}}{|\underline{t}|}),

is monogenic in ${{\bf R}}_{1}^{m}$ and $\hat{f}$ is the Fourier transform of $f$ given by (3).

The functions $u$ and $v$ are constructed by the Poisson and the conjugate Poisson integrals, respectively. That is,

\displaystyle u(\underline{x},s)=u*P_{s}(\underline{x})=\frac{1}{\omega_{m}}\int_{{\bf R}^{m}}\frac{s}{|s+(\underline{x}-\underline{t})|^{m+1}}u(\underline{t})d\underline{t}

(4)

and

\displaystyle\underline{v}(\underline{x},s)=v*Q_{s}(\underline{x})=\frac{1}{\omega_{m}}\int_{{\bf R}^{m}}\frac{\overline{\underline{x}-\underline{t}}}{|s+(\underline{x}-\underline{t})|^{m+1}}v(\underline{t})d\underline{t},

(5)

where $P_{s}(\underline{x}):=\frac{1}{\omega_{m}}\frac{s}{|s+\underline{x}|^{m+1}}$ and $Q_{s}(\underline{x}):=\frac{1}{\omega_{m}}\frac{\overline{\underline{x}}}{|s+\underline{x}|^{m+1}}$ are the Poisson and the conjugate Poisson kernel in ${\bf R}_{1}^{m,+}$ , respectively.

3 Local Attenuation and Local Phase Vector

Note that it is possible to write the monogenic scale function $f\in M^{2}({\bf R}^{m,+}_{1})$ in polar coordinate. Let us review the local feature [21] as follows.

Definition 3.1 (Local Features Representation I).

Suppose $f:=u+\underline{v}\in M^{2}({\bf R}^{m,+}_{1})$ has the polar form

\displaystyle f(\underline{x},s)=A(f)(\underline{x},s)e^{\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}\theta(\underline{x},s)},

(6)

then

\displaystyle A(f)(\underline{x},s):=|f(\underline{x},s)|=\sqrt{u(\underline{x},s)^{2}+|\underline{v}(\underline{x},s)|^{2}}

(7)

is is called the local amplitude.

\displaystyle\theta(\underline{x},s):=\arctan\left(\frac{|\underline{v}(\underline{x},s)|}{u(\underline{x},s)}\right)

(8)

is called the phase angle that is between $0$ and $\pi$ .

\displaystyle\underline{r}(\underline{x},s):=\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}\theta(\underline{x},s),

(9)

is called the local phase vector. ${\rm{Sc}}\left[(\underline{D}\theta(\underline{x},s))\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}\right]$ is called the directional phase derivative and $e^{\underline{r}(\underline{x},s)}$ is called the phase direction. The phase derivative or instantaneous frequency is defined by

\displaystyle{\rm{Sc}}\left[\left({\underline{D}f(\underline{x},s)}\right)\left({f(\underline{x},s)}\right)^{-1}\right].

(10)

Building on the ideas of [11], we can have the alternative form.

Definition 3.2 (Local Features Representation II).

For nontrivial function $f:=u+\underline{v}\in M^{2}({\bf R}^{m,+}_{1})$ , the local amplitude is nonzero. We can rewrite (6) as

\displaystyle f(\underline{x},s)=e^{a(\underline{x},s)+\underline{r}(\underline{x},s)},

(11)

where

\displaystyle a(\underline{x},s):=\ln A(f)(\underline{x},s)=\frac{1}{2}\ln(u(\underline{x},s)^{2}+|\underline{v}(\underline{x},s)|^{2})

(12)

is called the local attenuation.

Remark 3.1.

In one-dimensional case, $\frac{\underline{v}(x,s)}{|\underline{v}(x,s)|}={\bf i},$ therefore the local phase vector $\underline{r}(x,s)={\bf i}\theta(x,s)$ .

Suppose $f(x,s):=u(x,s)+{\bf i}v(x,s)\in H^{2}({{\bf C}^{+}})$ has the form (11). That is, $f(x,s)=e^{a(x,s)+{\bf i}\theta(x,s)}$ has no zeros and isolated singularities in the half plane ${\bf C}^{+}$ , then the local attenuation $a(x,s)=\frac{1}{2}\ln(u^{2}+v^{2})$ and the local phase $\theta(x,s)=\arctan\left(\frac{v}{u}\right)$ are related by the Cauchy-Riemann system. The reason is that the composition of analytic function is analytic. If $f(x,s)=u(x,s)+{\bf i}v(x,s)$ is analytic and has no zeros and isolated singularities in the half plane ${\bf C}^{+}$ , then $a(x,s)+{\bf i}\theta(x,s)$ is also analytic in ${\bf C}^{+}$ .

Using the Cauchy-Riemann system for $a(x,s)+{\bf i}\theta(x,s)$ , we have

\displaystyle\frac{\partial a}{\partial s}+\frac{\partial\theta}{\partial x}=0,

\displaystyle\frac{\partial a}{\partial x}-\frac{\partial\theta}{\partial s}=0

From the above system, we notice that:

1.

The instantaneous frequency $\frac{\partial\theta}{\partial x}$ can be obtained by the minus of the scale derivative of the local attenuation ${\partial a\over\partial s}$ .
2.

The zero points of the scale derivative of the local phase ${\partial\theta\over\partial s}$ is given by the extrema of the local attenuation.

Building on the ideas of 1D, the authors [11] considered the intrinsically 1D monogenic signals.

Definition 3.3.

If $f(x,s):=u(x,s)+{\bf i}v(x,s)\in H^{2}({\bf C}^{+})$ has no zeros and isolated singularities in the half plane ${\bf C}^{+}$ , then the intrinsically 1D monogenic signal is defined by

$\displaystyle f(<\underline{x},\underline{n}>,s)$	$\displaystyle=$	$\displaystyle u(<\underline{x},\underline{n}>,s)+\overline{\underline{n}}v(<\underline{x},\underline{n}>,s)$	(13)
	$\displaystyle=$	$\displaystyle u(<\underline{x},\underline{n}>,s)+\underline{v}(<\underline{x},\underline{n}>,s$
	$\displaystyle=$	$\displaystyle e^{a(<\underline{x},\underline{n}>,s)+\underline{r}(<\underline{x},\underline{n}>,s)},$

where $\underline{x},\underline{n}\in{\bf R}^{m}$ and $\underline{n}$ is a fixed unit vector. The local attenuation is given by $a(<\underline{x},\underline{n}>,s)=\frac{1}{2}\ln(u^{2}+v^{2})$ and the local phase vector is given by $\underline{r}(<\underline{x},\underline{n}>,s)=\overline{n}\arctan\left(\frac{v}{u}\right)$ for intrinsically 1D signal.

Felsberg et al. [11] proved that for the intrinsically 1D signals, the local attenuation $a(<\underline{x},\underline{n}>,s)$ and the local phase-vector $\underline{r}(<\underline{x},\underline{n}>,s)$ are related by the Hilbert transform pairs (1). Moreover, by the analyticity, using the generalized Cauchy-Riemann operator ${\partial\over\partial s}+\underline{D}$ on $a(<\underline{x},\underline{n}>,s)+\underline{r}(<\underline{x},\underline{n}>,s)$ , we have

\displaystyle\frac{\partial a}{\partial s}+\underline{D}(\underline{r})=0,

(14)

\displaystyle\underline{D}a+\frac{\partial\underline{r}}{\partial s}=0.

(15)

In [11], the local frequency of the intrinsically 1D signal and the differential phase congruency are defined by $\underline{D}(\underline{r})$ and $\frac{\partial\underline{r}}{\partial s}$ , receptively. From system (14) and (15), we notice that:

1.

The local frequency in an intrinsically 1D signal $\underline{D}(\underline{r})$ can also be obtained by the minus of the scale derivative of the local attenuation ${\partial a\over\partial s}$ .
2.

The zero points of the differential phase congruency ${\partial r\over\partial s}$ is given by the extrema of the local attenuation.

Remark 3.2.

In the recent paper [21], the instantaneous frequency of $f:=u+\underline{v}=e^{a+\underline{r}}$ is given by (10)

			$\displaystyle{\rm{Sc}}\left[({\underline{D}f(\underline{x},s)})({f(\underline{x},s)})^{-1}\right]$		(16)
		$\displaystyle=$	$\displaystyle{\rm{Sc}}\left[\left(\underline{D}\frac{\underline{v}}{\|\underline{v}\|}\right)\sin\theta(\underline{x},s)\cos\theta(\underline{x},s)\right]+{\rm{Sc}}\left[\left(\underline{D}\theta(\underline{x},s)\right)\frac{\underline{v}}{\|\underline{v}\|}\right].$		(16)

In particular, if $\frac{\underline{v}}{|\underline{v}|}$ is a constant, the first term in (16) vanishes, then the instantaneous frequency is ${\rm{Sc}}\left[\left(\underline{D}\theta(\underline{x},s)\right)\frac{\underline{v}}{|\underline{v}|}\right]$ . It coincides with the local frequency defined in [11]. That is, when $\frac{\underline{v}}{|\underline{v}|}=\overline{\underline{n}}$ is a constant, the local frequency $\underline{D}(\underline{r})$ is given by $(\underline{D}\theta(<\underline{x},\underline{n}>,s))\overline{\underline{n}}$ .

Remark 3.3.

1.

In Clifford analysis [7], we notice that if $f(x,s)=u(x,s)+{\bf i}v(x,s)\in H^{2}({\bf C}^{+})$ , then for fixed unit vector $\underline{n}\in{\bf R}^{m}$ , the function

$f(<\underline{x},\underline{n}>,s)=u(<\underline{x},\underline{n}>,s)+\overline{\underline{n}}v(<\underline{x},\underline{n}>,s),$

is monogenic in ${\bf C}^{+}$ . It is called monogenic plane wave.
2.

Clearly, if $f(x,s)=e^{a(x,s)+i\theta(x,s)}\in H^{2}({\bf C}^{+})$ has no zeros and isolated singularities in ${\bf C}^{+}$ , then $a(x,s)+{\bf i}\theta(x,s)$ is also analytic in ${\bf C}^{+}$ . Consequently, the function $a(<\underline{x},\underline{n}>,s)+\overline{\underline{n}}\theta(<\underline{x},\underline{n}>,s)=a(<\underline{x},\underline{n}>,s)+\underline{r}(<\underline{x},\underline{n}>,s)$ is monogenic in ${\bf R}^{m,+}_{1}$ .

Problem 3.1.

What is the situation in the higher dimension if the signal is not intrinsically 1D signal?

The solution was not considered in [11] and [9]. While in higher dimension, the situation is more complicated. The theory does not hold in general. In fact, if $f(\underline{x},s)=u(\underline{x},s)+\underline{v}(\underline{x},s)=e^{a(\underline{x},s)+\underline{r}(\underline{x},s)}$ is monogenic in the half space ${\bf R}_{1}^{m,+}$ , $a(\underline{x},s)+\underline{r}(\underline{x},s)$ is not monogenic in general. Therefore, the local attenuation $a$ and the local phase vector $r$ are not related by the generalized Cauchy-Riemann system in higher dimensions. Let us now look at an example to illustrate the topic discussed above.

Example 3.1.

Let $f(\underline{x},s)=\frac{s}{|s+\underline{x}|^{m+1}}+\frac{\overline{\underline{x}}}{|s+\underline{x}|^{m+1}}=E(s+\underline{x})$ be the Cauchy kernel in ${\bf R}_{1}^{m}\setminus\{0\}$ , which is monogenic in ${\bf R}_{1}^{m}\setminus\{0\}$ . Then, by straightforward computations, we have

a(\underline{x},s)+\underline{r}(\underline{x},s)=-\frac{m}{2}\ln(s^{2}+|\underline{x}|^{2})+\frac{\overline{\underline{x}}}{|\underline{x}|}\arctan\left(\frac{|\underline{x}|}{s}\right).

Then we apply the generalized Cauchy-Riemann operator ${\partial\over\partial s}+\underline{D}$ on it, we have

			$\displaystyle\left(\frac{\partial}{\partial s}+\underline{D}\right)\left[-\frac{m}{2}\ln(s^{2}+\|\underline{x}\|^{2})+\frac{\overline{\underline{x}}}{\|\underline{x}\|}\arctan\left(\frac{\|\underline{x}\|}{s}\right)\right]$
		$\displaystyle=$	$\displaystyle\frac{(1-m)(s+\underline{x})}{s^{2}+\|\underline{x}\|^{2}}+\frac{m-1}{\|\underline{x}\|}\arctan\left(\frac{\|\underline{x}\|}{s}\right)\neq 0.$

Therefore, $a(\underline{x},s)+\underline{r}(\underline{x},s)$ is not monogenic.

Let us now describe the solution for Problem 3.1, Theorem 3.1 gives the relationship between the local phase vector $r$ and the local attenuation $a$ in higher dimensional spaces.

Theorem 3.1.

Let $f(\underline{x},s)=u(\underline{x},s)+\underline{v}(\underline{x},s)=e^{a(\underline{x},s)+\underline{r}(\underline{x},s)}\in M^{2}({\bf R}_{1}^{m,+})$ , where $a(\underline{x},s)$ and $\underline{r}(\underline{x},s)$ are the local attenuation and the local phase-vector defined by (12) and (9), respectively. If $f$ has no zeros and isolated singularities in the half space ${\bf R}_{1}^{m,+}$ . Then we have

\displaystyle\frac{\partial a}{\partial s}+{\rm Sc}[(\underline{D}e^{\underline{r}})e^{-\underline{r}}]=0,

(17)

\displaystyle\frac{\partial\underline{r}}{\partial s}+\underline{D}a-{\rm Vec}\left[\left(\underline{D}\frac{\underline{v}}{|\underline{v}|}\right)\frac{\underline{v}}{|\underline{v}|}\right]\sin^{2}\theta+(\sin\theta\cos\theta-\theta){\partial\over\partial s}\left({\underline{v}\over|\underline{x}|}\right)=0.

(18)

In particular, if $\frac{\underline{v}}{|\underline{v}|}$ is independent of $s$ , that is ${\partial\over\partial s}\left({\underline{v}\over|\underline{x}|}\right)=0$ , then we have the following corollary.

Corollary 3.1.

Let $f(\underline{x},s)=u(\underline{x},s)+\underline{v}(\underline{x},s)=e^{a(\underline{x},s)+\underline{r}(\underline{x},s)}\in M^{2}({\bf R}_{1}^{m,+})$ , where $a$ and $\underline{r}$ are the local attenuation and local phase-vector defined by (12) and (9), respectively. If $f$ has no zeros and isolated singularities in the half space ${\bf R}_{1}^{m,+}$ and the local orientation $\frac{\underline{v}}{|\underline{v}|}$ does not change through scale $s$ , then we have

\displaystyle\frac{\partial\underline{r}}{\partial s}+\underline{D}a-{\rm Vec}[(\underline{D}\frac{\underline{v}}{|\underline{v}|})\frac{\underline{v}}{|\underline{v}|}]\sin^{2}\theta=0.

(19)

Combining (17), (18) and (16), we conclude that

Theorem 3.2.

[Instantaneous Frequency]

1.

The instantaneous frequency in higher dimensional spaces defined by (10) is equal to the minus of the scale derivative of the local attenuation ${\partial a\over\partial s}$ .
2.

The zero points of the differential phase congruency ${\partial r\over\partial s}$ is not equal to the extrema of the local attenuation.

Remark 3.4.

By Theorem 3.2, we notice that, like in one dimensional case, the phase derivative in higher dimensions can also be given by the minus of the scale derivative of the local attenuation. However, the zero points of the phase congruency is not equal to the extrema of the local attenuation in high dimensional case. The nonzero extra term

-{\rm Vec}\left[\left(\underline{D}\frac{\underline{v}}{|\underline{v}|}\right)\frac{\underline{v}}{|\underline{v}|}\right]\sin^{2}\theta+(\sin\theta\cos\theta-\theta){\partial\over\partial s}\left({\underline{v}\over|\underline{x}|}\right)

appears in high dimensional cases.

We have divided the proof of Theorem 3.1 into a series of lemmas.

Lemma 3.1.

\displaystyle{\rm{Sc}}\left[(\frac{\partial}{\partial s}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}\right]=0.

(20)

Proof: By the generalized Euler formula $e^{\underline{r}(\underline{x},s)}=e^{\frac{\underline{v}}{|\underline{v}|}\theta}=\cos\theta+\frac{\underline{v}}{|\underline{v}|}\sin\theta$ , we have

	$\displaystyle\left(\frac{\partial}{\partial s}e^{\underline{r}(\underline{x},s)}\right)e^{-\underline{r}(\underline{x},s)}$	(21)
$\displaystyle=$	$\displaystyle\frac{\partial}{\partial s}\left(\cos\theta+\frac{\underline{v}}{\|\underline{v}\|}\sin\theta\right)\left(\cos\theta-\frac{\underline{v}}{\|\underline{v}\|}\sin\theta\right)$
$\displaystyle=$	$\displaystyle\left(-\sin\theta\frac{\partial\theta}{\partial s}+\frac{\partial\frac{\underline{v}}{\|\underline{v}\|}}{\partial s}\sin\theta+\frac{\underline{v}}{\|\underline{v}\|}\cos\theta\frac{\partial\theta}{\partial s}\right)\left(\cos\theta-\frac{\underline{v}}{\|\underline{v}\|}\sin\theta\right)$
$\displaystyle=$	$\displaystyle\frac{\underline{v}}{\|\underline{v}\|}\frac{\partial\theta}{\partial s}+\sin\theta\cos\theta\frac{\partial\frac{\underline{v}}{\|\underline{v}\|}}{\partial s}-\sin^{2}\theta\frac{\partial\frac{\underline{v}}{\|\underline{v}\|}}{\partial s}\frac{\underline{v}}{\|\underline{v}\|}.$

Clearly, the scalar part of $\left(\frac{\partial}{\partial s}e^{\underline{r}(\underline{x},s)}\right)e^{-\underline{r}(\underline{x},s)}$ is decided by the third part of equation (21). Let us now prove the following

{\rm{Sc}}\left[(\frac{\partial\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}}{\partial s})\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}\right]=0.

Denote $\underline{I}(\underline{x},s):=\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}$ , we have $\underline{I}^{2}(\underline{x},s)=-1$ . Then $\frac{\partial[{\underline{I}^{2}(\underline{x},s)}]}{\partial s}=0$ . By equation (2), we have

	$\displaystyle\frac{\partial[{\underline{I}^{2}(\underline{x},s)}]}{\partial s}$	$\displaystyle=$	$\displaystyle\frac{\partial[{\underline{I}(\underline{x},s)}]}{\partial s}\underline{I}(\underline{x},s)+\underline{I}(\underline{x},s)\frac{\partial[{\underline{I}(\underline{x},s)}]}{\partial s}$
		$\displaystyle=$	$\displaystyle 2{\rm{Sc}}[\frac{\partial{\underline{I}(\underline{x},s)}}{\partial s}\underline{I}(\underline{x},s)]=0.$

Thus, we obtain the desired result.

Lemma 3.2.

\displaystyle{\rm{Vec}}\left[(\frac{\partial}{\partial s}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}\right]=(\sin\theta\cos\theta-\theta)\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}+\frac{\partial\underline{r}}{\partial s}.

(22)

\displaystyle{\rm{Vec}}\left[(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}\right]=-\sin^{2}\theta{\rm{Vec}}\left[(\underline{D}\frac{\underline{v}}{|\underline{v}|})\frac{\underline{v}}{|\underline{v}|}\right].

(23)

Proof: From (21), we know that the vector part of $(\frac{\partial}{\partial s}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}$ is decided by $\frac{\underline{v}}{|\underline{v}|}\frac{\partial\theta}{\partial s}+\sin\theta\cos\theta\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}$ . Since $\underline{r}=\frac{\underline{v}}{|\underline{v}|}\theta$ , we have

\frac{\partial\underline{r}}{\partial s}=\frac{\partial\theta}{\partial s}\frac{\underline{v}}{|\underline{v}|}+\theta\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}.

Therefore, we obtain equation (22).

To prove equation (23), by direct calculation, we have

	$\displaystyle\left(\underline{D}e^{\underline{r}(\underline{x},s)}\right)e^{-\underline{r}(\underline{x},s)}$	(24)
$\displaystyle=$	$\displaystyle\underline{D}\left(\cos\theta+\frac{\underline{v}}{\|\underline{v}\|}\sin\theta\right)\left(\cos\theta-\frac{\underline{r}}{\|\underline{r}\|}\sin\theta\right)$
$\displaystyle=$	$\displaystyle\left[-\sin\theta(\underline{D}\theta)+(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})\sin\theta+\cos\theta(\underline{D}\theta)\frac{\underline{v}}{\|\underline{v}\|}\right]\left[\cos\theta-\frac{\underline{v}}{\|\underline{v}\|}\sin\theta\right]$
$\displaystyle=$	$\displaystyle\frac{\underline{v}}{\|\underline{v}\|}(\underline{D}\theta)+\sin\theta\cos\theta(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})-\sin^{2}\theta(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})\frac{\underline{v}}{\|\underline{v}\|}.$

The fist part and the second part of equation (24) are scalar and bi-vector, respectively. Therefore the vector part of $(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}$ is decided by the third part of equation (24). Thus we obtain (23).

We can now prove Theorem 3.1.

Proof of Theorem 3.1: Since $f(\underline{x},s)=e^{a(\underline{x},s)+\underline{r}(\underline{x},s)}\in M^{2}({\bf R}_{1}^{m,+})$ , we have

\left(\frac{\partial}{\partial s}+\underline{D}\right)e^{a(\underline{x},s)+\underline{r}(\underline{x},s)}=0.

By straightforward computation, we have

e^{a(\underline{x},s)}\frac{\partial a(\underline{x},s)}{\partial s}e^{\underline{r}(\underline{x},s)}+e^{a(\underline{x},s)}\frac{\partial e^{\underline{r}(\underline{x},s)}}{\partial s}+e^{a(\underline{x},s)}[\underline{D}a(\underline{x},s)]e^{\underline{r}(\underline{x},s)}+e^{a(\underline{x},s)}(\underline{D}e^{\underline{r}(\underline{x},s)})=0.

That is

\displaystyle\frac{\partial a(\underline{x},s)}{\partial s}+\frac{\partial e^{\underline{r}(\underline{x},s)}}{\partial s}e^{-\underline{r}(\underline{x},s)}+\underline{D}a(\underline{x},s)+(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}=0.

(25)

Therefore, the scalar part of (25) is zero. By combining Lemma 3.1, we have

			$\displaystyle{\rm Sc}\left[\frac{\partial a(\underline{x},s)}{\partial s}+\frac{\partial e^{\underline{r}(\underline{x},s)}}{\partial s}e^{-\underline{r}(\underline{x},s)}+\underline{D}a(\underline{x},s)+(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}\right]$
		$\displaystyle=$	$\displaystyle\frac{\partial a(\underline{x},s)}{\partial s}+{\rm Sc}[\frac{\partial e^{\underline{r}(\underline{x},s)}}{\partial s}e^{-\underline{r}(\underline{x},s)}]+{\rm Sc}[(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}]$
		$\displaystyle=$	$\displaystyle\frac{\partial a(\underline{x},s)}{\partial s}+{\rm Sc}[(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}]=0.$

Therefore, we get the desired result (17).

The vector part of Eq. (25) is also zero. By using Lemma 3.2, we obtain

			$\displaystyle{\rm Vec}\left[\frac{\partial a(\underline{x},s)}{\partial s}+\frac{\partial e^{\underline{r}(\underline{x},s)}}{\partial s}e^{-\underline{r}(\underline{x},s)}+\underline{D}a(\underline{x},s)+(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}\right]$
		$\displaystyle=$	$\displaystyle{\rm Vec}\left[\frac{\partial e^{\underline{r}(\underline{x},s)}}{\partial s}e^{-\underline{r}(\underline{x},s)}\right]+\underline{D}a(\underline{x},s)+{\rm Vec}\left[(\underline{D}e^{\underline{r}(\underline{x},s)})e^{-\underline{r}(\underline{x},s)}\right]$
		$\displaystyle=$	$\displaystyle\frac{\partial\underline{r}}{\partial s}+\underline{D}a+(\sin\theta\cos\theta-\theta)\frac{\partial\frac{\underline{v}}{\|\underline{v}\|}}{\partial s}-{\rm{Vec}}\left[(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})\frac{\underline{v}}{\|\underline{v}\|}\right]\sin^{2}\theta=0.$

This completes the proof.

If $f\in M^{2}({\bf R}_{1}^{m,+})$ has the axial form

f(\underline{x},s)=u(\rho,s)+\frac{\overline{\underline{x}}}{|\underline{x}|}v(\rho,s),\mbox{ }\rho=|\underline{x}|.

Then, in this case, the local orientation $\frac{\underline{v}}{|\underline{v}|}=\frac{\overline{\underline{x}}}{|\underline{x}|}$ does not change through the scale $s$ . By polar coordinate, $f(\underline{x},s)=e^{a(\rho,s)+\frac{\overline{\underline{x}}}{|\underline{x}|}\theta(\rho,s)}$ , using Theorem 3.1, we have the following corollary.

Corollary 3.2.

Let $f(\underline{x},s)=u(\rho,s)+\frac{\overline{\underline{x}}}{|\underline{x}|}v(\rho,s)=e^{a(\rho,s)+\frac{\overline{\underline{x}}}{|\underline{x}|}\theta(\rho,s)}\in M^{2}({\bf R}_{1}^{m,+})$ . Then we have

\displaystyle-\frac{\partial a}{\partial s}=\frac{\partial\theta}{\partial\rho}+\frac{m-1}{\rho}\sin\theta\cos\theta,

(28)

\displaystyle\frac{\partial\theta}{\partial s}=\frac{\partial a}{\partial\rho}+\frac{m-1}{\rho}\sin^{2}\theta.

(29)

It is easy to see that when $m=1$ , the above system ((28) and (29)) is just the Cauchy-Riemann system in the one dimensional case.

4 Edge Detection Methods

Edge detection by means of quadrature filters has two ways: either by detecting local maxima of the local amplitude or by detecting points of stationary phase in scale-space. In this section, we begin by reviewing the differential phase congruency method [11].

4.1 Differential Phase Congruency Methods

Method 4.1 (DPC).

For intrinsically 1D monogenic signal $f\in H^{2}({\bf C}^{+})$ given by (13), if $f$ has no zero and isolated singularities in the half plane ${\bf C}^{+}$ . Then the differential phase congruency (DPC) has the following formula

\displaystyle\frac{\partial\underline{r}_{in1}(\underline{x},s)}{\partial s}=\frac{u(\underline{x},s)\frac{\partial\underline{v}(\underline{x},s)}{\partial s}-\underline{v}(\underline{x},s)\frac{\partial u(\underline{x},s)}{\partial s}}{u(\underline{x},s)^{2}+|\underline{v}(\underline{x},s)|^{2}}=0,

(30)

where $\underline{r}_{in1}(\underline{x},s):=\underline{r}(<\underline{x},\underline{n}>,s)$ .

By (15), we notice that formula (30) can also be obtained by the $-\underline{D}a$ . However, the zero points of the differential phase congruency is not given by the extrema of the local attenuation in higher dimension.

4.2 Proposed Methods

Let’s introduce the local attenuation (LA) method for monogenic signals.

Method 4.2 (LA).

For $f\in M^{2}({\bf R}_{1}^{m,+})$ has no zeros and isolated singularities in the half space ${\bf R}^{m,+}_{1}$ , the local attenuation has the formula

\displaystyle\underline{D}a(\underline{x},s)=\frac{u(\underline{x},s)\underline{D}[u(\underline{x},s)]+|\underline{v}(\underline{x},s)|\underline{D}[|\underline{v}(\underline{x},s)|]}{u^{2}(\underline{x},s)+|\underline{v}(\underline{x},s)|^{2}}.

(31)

Applying Dirac operator $\underline{D}$ on the local attenuation $a$ , by direct computation on (7), formula (31) follows. Using (15), we know that for intrinsically 1D signals, the zero points of the differential phase congruency is given by the extrema of the local attenuation. Notice that formula (31) is equivalent to (30) for the intrinsically 1D monogenic signal.

Our second method is the modified differential phase congruency (MDPC) method. To proceed, we need the following technical lemma.

Lemma 4.1.

\displaystyle\frac{\partial\underline{r}(\underline{x},s)}{\partial s}=\left(\theta-\sin\theta\cos\theta\right)\frac{\partial\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}}{\partial s}+\frac{u(\underline{x},s)\frac{\partial\underline{v}(\underline{x},s)}{\partial s}-\underline{v}(\underline{x},s)\frac{\partial u(\underline{x},s)}{\partial s}}{u^{2}(\underline{x},s)+|\underline{v}(\underline{x},s)|^{2}}.

(32)

Proof: By Eq. (9), we have

\displaystyle\frac{\partial\underline{r}(\underline{x},s)}{\partial s}

\displaystyle=

\displaystyle\frac{\partial}{\partial s}(\frac{\underline{v}}{|\underline{v}|}\theta)=\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}\theta+\frac{\underline{v}}{|\underline{v}|}\frac{\partial}{\partial s}\theta.

(33)

By straightforward computation, we have

\displaystyle\frac{\partial}{\partial s}\theta=\frac{\partial}{\partial s}(\arctan\left(\frac{|\underline{v}|}{u}\right))=\frac{\frac{\partial|\underline{v}|}{\partial s}u-|\underline{v}|\frac{\partial u}{\partial s}}{u^{2}+|\underline{v}|^{2}}.

(34)

Then,

\displaystyle\frac{\underline{v}}{|\underline{v}|}\frac{\partial}{\partial s}(\arctan\left(\frac{|\underline{v}|}{u}\right))=\frac{\frac{\underline{v}}{|\underline{v}|}\frac{\partial|\underline{v}|}{\partial s}u-\underline{v}\frac{\partial u}{\partial s}}{u^{2}+|\underline{v}|^{2}}.

(35)

Using the equation

\displaystyle\frac{\partial\underline{v}}{\partial s}=\frac{\partial}{\partial s}(\frac{\underline{v}}{|\underline{v}|}|\underline{v}|)=|\underline{v}|\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}+\frac{\underline{v}}{|\underline{v}|}\frac{\partial|\underline{v}|}{\partial s},

we obtain

\displaystyle\frac{\underline{v}}{|\underline{v}|}\frac{\partial|\underline{v}|}{\partial s}=\frac{\partial\underline{v}}{\partial s}-|\underline{v}|\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}.

(36)

Applying Eq. (36) to Eq. (35), we have

\displaystyle\frac{\underline{v}}{|\underline{v}|}\frac{\partial}{\partial s}(\arctan\left(\frac{|\underline{v}|}{u}\right))=\frac{u\frac{\partial\underline{v}}{\partial s}-\underline{v}\frac{\partial u}{\partial s}}{u^{2}+|\underline{v}|^{2}}-\frac{u|\underline{v}|}{u^{2}+|\underline{v}|^{2}}\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}.

(37)

Combining Eq. (33) and Eq. (37), we obtain Eq. (32).

Remark 4.1.

Note that equation (30) is a special case of (32). The reason is in the intrinsically 1D neighborhood, the local orientation $\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}=\overline{\underline{n}}$ is a constant. So $\frac{\partial\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}}{\partial s}=0$ . In fact, formula (30) always holds if the local orientation is independent of $s$ .

Let us now define the points of modified differential phase congruency.

Definition 4.1.

Let $\underline{r}(\underline{x},s)$ be the local phase vector, given by (9), of function $f\in M^{2}({\bf R}_{1}^{m,+}$ . Points where

\frac{\partial\underline{r}(\underline{x},s)}{\partial s}-{\rm Vec}[(\underline{D}\frac{\underline{v}}{|\underline{v}|})\frac{\underline{v}}{|\underline{v}|}]\sin^{2}\theta+(\sin\theta\cos\theta-\theta)\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}=0

are called points of modified differential phase congruency (MDPC).

Remark 4.2.

From Theorem 3.1 we know that in any higher dimensional cases, edge detection by means of local amplitude maxima is equivalent to edge detection by modified differential phase congruency.

Using Eq. (32), we can now proposed our second method, the so-called modified differential phase congruency (MDPC) method.

Method 4.3 (MDPC).

For $f\in M^{2}({\bf R}_{1}^{m,+})$ has no zeros and isolated singularities in the half space ${\bf R}^{m,+}_{1}$ , the MDPC has the formula

			$\displaystyle\frac{\partial\underline{r}(\underline{x},s)}{\partial s}-{\rm Vec}[(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})\frac{\underline{v}}{\|\underline{v}\|}]\sin^{2}\theta+\left(\sin\theta\cos\theta-\theta\right)\frac{\partial\frac{\underline{v}}{\|\underline{v}\|}}{\partial s}$		(38)
		$\displaystyle=$	$\displaystyle\frac{u(\underline{x},s)\frac{\partial\underline{v}(\underline{x},s)}{\partial s}-\underline{v}(\underline{x},s)\frac{\partial u(\underline{x},s)}{\partial s}}{u^{2}(\underline{x},s)+\|\underline{v}(\underline{x},s)\|^{2}}-{\rm Vec}\left[(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})\frac{\underline{v}}{\|\underline{v}\|}\right]\sin^{2}\theta.$		(38)

Finally, we introduce a mixed method by combining local attenuation and modified differential phase congruency (LA+MDPC) for edge detection.

Method 4.4 (LA+MDPC).

For $f\in M^{2}({\bf R}_{1}^{m,+})$ has no zeros and isolated singularities in the half space ${\bf R}^{m,+}_{1}$ , the MDPC has the formula

			$\displaystyle\frac{\partial\underline{r}}{\partial s}-\underline{D}a-{\rm Vec}[(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})\frac{\underline{v}}{\|\underline{v}\|}]\sin^{2}\theta+(\sin\theta\cos\theta-\theta)\frac{\partial\frac{\underline{v}}{\|\underline{v}\|}}{\partial s}$		(39)
		$\displaystyle=$	$\displaystyle\frac{u(\underline{x},s)\frac{\partial\underline{v}(\underline{x},s)}{\partial s}-\underline{v}(\underline{x},s)\frac{\partial u(\underline{x},s)}{\partial s}}{u^{2}(\underline{x},s)+\|\underline{v}(\underline{x},s)\|^{2}}-\underline{D}a-{\rm Vec}[(\underline{D}\frac{\underline{v}}{\|\underline{v}\|})\frac{\underline{v}}{\|\underline{v}\|}]\sin^{2}\theta$		(39)

5 Experiments

In this section, we begin by showing the details of our proposed methods. Two classical edge detection methods, such as Canny and Sobel edge detectors, will be compared with our algorithms. The Canny edge detector will begin by applying Gaussian filter to the test images. Then Canny edge detector computes the gradients on the images. For the Sobel edge detector, we only apply its gradients to the original test images. For the DPC and our proposed methods, we first apply the Poisson filter to the test images, then we compute apply their formulas to the images.

By comparing with the classical methods, phased based methods may show more detail on image.

Refer to caption — Figure 1: Original images

5.1 Algorithms

Let us now give the details of the phase based algorithms. They are divided by the following steps.

Step 1.

Input image $f(\underline{x})$ . For simplicity, the color image is converted to the gray image.

Step 2.

Poisson filtering: $u(\underline{x},s)=f*P_{s}(\underline{x})$ and and $\underline{v}(\underline{x},s)=f*Q_{s}(\underline{x})$ for a fixed scale $s>0$ . We will discuss how to choose $s$ in Section LABEL:5.1. We consider $s$ in $0.1$ , $0.5$ , $1.0$ and $5.0$ . Moreover, we choose $s=0.5$ for all test images to compare with different methods.

Step 3.

Compute the local attenuation $a(\underline{x},s)=\frac{1}{2}\ln(u(\underline{x},s)^{2}+|\underline{v}(\underline{x},s)|^{2})$ and local phase vector $\underline{r}(\underline{x},s)=\frac{\underline{v}(\underline{x},s)}{|\underline{v}(\underline{x},s)|}\theta(\underline{x},s)$ , where the phase angle is given by $\theta(\underline{x},s)=\arctan\left(\frac{|\underline{v}(\underline{x},s)|}{u(\underline{x},s)}\right)$ .

Step 4.

Compute gradients by different methods to get the gradient maps.

The differential phase congruency (DPC) method: compute $\frac{\partial\underline{r}_{in1}(\underline{x},s)}{\partial s}$ by the formula

\frac{u(\underline{x},s)\frac{\partial\underline{v}(\underline{x},s)}{\partial s}-\underline{v}(\underline{x},s)\frac{\partial u(\underline{x},s)}{\partial s}}{u(\underline{x},s)^{2}+|\underline{v}(\underline{x},s)|^{2}}.

The local amplitude (LA) method: compute $\underline{D}a(\underline{x},s)$ , where $\underline{D}$ is the sum for the derivatives of image in vertical and horizontal directions. By theoretical analysis, $\underline{D}a(\underline{x},s)$ can be computed by

\frac{u(\underline{x},s)\underline{D}[u(\underline{x},s)]+|\underline{v}(\underline{x},s)|\underline{D}[|\underline{v}(\underline{x},s)|]}{u^{2}(\underline{x},s)+|\underline{v}(\underline{x},s)|^{2}}.

The modified differential phase congruency (MDPC) method: compute $\frac{\partial\underline{r}(\underline{x},s)}{\partial s}-{\rm Vec}[(\underline{D}\frac{\underline{v}}{|\underline{v}|})\frac{\underline{v}}{|\underline{v}|}]\sin^{2}\theta+\left(\sin\theta\cos\theta-\theta\right)\frac{\partial\frac{\underline{v}}{|\underline{v}|}}{\partial s}$ , which equals to

\frac{u(\underline{x},s)\frac{\partial\underline{v}(\underline{x},s)}{\partial s}-\underline{v}(\underline{x},s)\frac{\partial u(\underline{x},s)}{\partial s}}{u^{2}(\underline{x},s)+|\underline{v}(\underline{x},s)|^{2}}-{\rm Vec}\left[(\underline{D}\frac{\underline{v}}{|\underline{v}|})\frac{\underline{v}}{|\underline{v}|}\right]\sin^{2}\theta.

The mixed method by using local attenuation and modified differential phase congruency (LA+MDCP): compute

\frac{u(\underline{x},s)\frac{\partial\underline{v}(\underline{x},s)}{\partial s}-\underline{v}(\underline{x},s)\frac{\partial u(\underline{x},s)}{\partial s}}{u^{2}(\underline{x},s)+|\underline{v}(\underline{x},s)|^{2}}-\underline{D}a-{\rm Vec}[(\underline{D}\frac{\underline{v}}{|\underline{v}|})\frac{\underline{v}}{|\underline{v}|}]\sin^{2}\theta.

Step 5.

Applying Non-maximum suppress to these gradient maps, which is the same as for the Canny edge detector. After non-maximum suppression, the edge will become thinner [22, 23]. For a fair comparison, all the six methods aforementioned utilize the non-maximum suppression method with the same parameters. Concretely, we choose the radius $r=1.5$ and the lower and upper threshold values are $1.0$ and $3.5$ , respectively.

5.2 Experiment Results

We will use three different images (Fig. 1) for the comparison of different edge detectors. Fig. 3 shows the edge detection results of various methods with the fixed scale $s=0.5$ . From top to down of Fig. 2, there are six rows. Each row shows one comparison method. They are Canny, Sobel, DPC, LA, MDPC and LA+MDPC methods, respectively. From the experiment results, we can draw the following conclusions.

1.

First, the mixed method yields decent edge detection results with fewer mistakes, outperforming other algorithms in some cases.
2.

Second, the comparison between the results of LA, DPC, MDPC and the mixed methods also suggest that both local attenuation and local phase are important in edge detection.
3.

Our proposed method MDPC can achieve very good performances in dealing with the details. For the pepper in Fig. 2, we found that our method and canny’s results are similar, we can find the edge of pepper in the results. However, for the shadows of the house and liver in Fig. 2, where the human eye is relatively subtle. Fortulately, the DPC and MDPC methods have found the details in the shadow. However, Canny, Sobel and LA methods cannot give the information about the shadows. Canny uses the Gaussian filter which will make part of these shadows as noise and removed them. While Sobel directly generate the horizontal and vertical differences, because the shadows and the surrounding area is not much difference, which may not find the shadow of the image. By applying the phase based method, these details can be clearly found in our experiment results. This shows that our method can detect the whole smooth region and local small change region. Applications can be useful in places where it is difficult for the human eye to find the details.

Remark 5.1.

For intrinsically 1D signals, we know that edge detection by means of local amplitude maxima is equivalent to edge detection by phase congruency. While, in intrinsically 2D signals, we know that it dose not hold. In [11], the authors said: $``$ We cannot give an exhaustive answer to this question. In this paper, since the behavior of phase and attenuation in intrinsically 2D neighborhoods is still work in progress.” From Eq. (38), we know that difference between the DPC and MDPC methods is ${\rm Vec}\left[(\underline{D}\frac{\underline{v}}{|\underline{v}|})\frac{\underline{v}}{|\underline{v}|}\right]\sin^{2}\theta$ . By experiment, we know that the effect is not obvious.

5.3 Effect of Scale

Monogenic signals at any scale $s>0$ form the monogenic scale-space $M^{2}({\bf R}^{m,+}_{1})$ . The representation of monogenic scale-space is just a monogenic function in the upper half space ${\bf R}^{m,+}_{1}$ . Therefore, considering the monogenic scale space instead of monogenic signal, it has a scale parameter $s>0$ to choose which provides us more analysis tools for different purposes.

We found that when $s$ tends to $0$ , the Poisson integral tends to Hilbert transform. Moreover, in the paper [18], Hilbert transform has been proved to be useful for image edge extraction. In the choice of scale, we compare $s$ from $0.1$ to $5$ , as can be seen in Fig. 2, when $s$ is smaller, the edge is more fine. But too much detail is not good at all, so in the comparative experiment in Fig. 3, we chose the case of 0.5 for $s$ to compare with the various methods.

6 Acknowledgements

The authors acknowledge financial support from the National Natural Science Funds (No. 11401606) and University of Macau No. MYRG2015-00058-L2-FST and the Macao Science and Technology Development Fund FDCT/099/2012/A3.

References

[1] T. Alieva and M.J. Bastiaans. Powers of transfer matrices determined by means of eigenfunctions, J. Opt. Soc. Amer. A, 16(10): 2413-2418 (1999).
[2] F. Brackx, R. Delanghe, and F. Sommen. Clifford Analysis, Pitman 76, Boston-London-Melbourne, 1982.
[3] J. Babaud, A. P. Witkin, M. Baudin, and R. O. Duda. Uniqueness of the Gaussian kernel for scale-space filtering, IEEE Trans. on Pattern Analysis and Machine Intelligence, 8(1): 26-33 (1986).
[4] C. Li, A. Mcintosh, and T. Qian. Clifford algebras, Fourier transforms, and singular convolution operators on Lipschitz surfaces, Rev. Mat. Iberoamericana, 10: 665-721 (1994).
[5] L. Cohen. Time-Frequency Analysis: Theory and Applications, Prentice Hall, Inc. Upper Saddle River, NJ, USA, 1995.
[6] P. Dang, T. Qian, and Z. You, Hardy-Sobolev spaces decomposition in signal analysis, J. Fourier Anal. Appl. , 17(1): 36-64 (2011).
[7] R. Delanghe, F. Sommen, and V. Soucek. Clifford Algebra and Spinor Valued Functions, Kluwer, Dordrecht-Boston-London, 1992.
[8] M. Felsberg. Low-level Image Prossing with the Structure Multivector, Institut Fur Informatik Und Praktische Mathmatik, Christian-Alberchis-Universitat Kiel, 2002.
[9] M. Felsberg, R. Duits, and L. Florack. The Monogenic Scale Space on a Rectanular Domain and it Features, Internertional Journal of Computer Vision, 64(2/3): 187-201 (2005).
[10] M. Felsberg and G. Sommer. The monogenic signal, IEEE Transactions on Signal Processing, 49(12): 3136-3144 (2001).
[11] M. Felsberg and G. Sommer. The Monogenic Scale-Space: A Unifying Approach to Phase-Based Image Processing in Scale-Space, Journal of Mathematical Imaging and Vision, 21(1): 5-26 (2004).
[12] J. Garnet. Bounded Analytic Functions, Academic Press, New York, 1981.
[13] S. L. Hahn. Hilbert Transforms in Signal Processing, Artech House, Norwood, MD, 1996.
[14] T. Iijima. Basic theory of pattern obsercation, In Papers of Technical Group on Automata and Automatic Control, IECE, Japan, (1959).
[15] T. Lindeberg. Scale-Space Theory in Computer Vision, The Kluwer International Series in Engeneering and Computer Science. Kluwer Academic Publishers, Boston, 1994.
[16] Nussenzveig, H. M. Causality and dispersion relations, London: Academic Press.
[17] K-I. Kou and T. Qian. The Paley-Wiener Theorem in ${{\bf R}^{m}}$ with the Clifford Analysis Setting, Journal of Functional Analysis, 189: 227-241 (2002).
[18] S. Pei, J. Huang, J. Ding, and G. Guo. Short response Hilbert transform for edge detection, IEEE Asia Pacific Conference on Circuits and Systems, 340-343 (2008).
[19] Y. Yang and K-I. Kou. Uncertainty principles for hypercomplex signals in the linear canonical transform domains, Signal Processing, 95: 67-75 (2014).
[20] Y. Yang, P. Dang, and T. Qian. Space-Frequency Analysis in Higher Dimensions and Applications, Annali di Mathematica, 194(4): 953-968 (2015).
[21] Y. Yang, T. Qian, and F. Sommen. Phase derivative of monogenic functions in higher dimensional spaces, Complex Analysis and Operator Theory, 6(5): 987-1010 (2012).
[22] P. Kovesi. Image Features From Phase Congruency, Videre: A Journal of Computer Vision Research. MIT Press, 1(3):1-30 (1999).
[23] P. Kovesi. Phase Congruency Detects Corners and Edges, Proceedings DICTA, (2013).
[24] L. Zhang, W. Zuo and D. Zhang. LSDT: Latent Sparse Domain Transfer Learning for Visual Adaptation, IEEE Trans. Image Processing, 25(3): 247-259 (2016).
[25] L. Zhang and D. Zhang. Visual Understanding via Multi-Feature Shared Learning with Global Consistency, IEEE Trans. Multimedia, 18(2): 247-259 (2016).
[26] L. Zhang and D. Zhang. Robust Visual Knowledge Transfer via Extreme Learning Machine-Based Domain Adaptation, IEEE Trans. Image Processing, 25(10): 4959-4973 (2016).