# A Multi-Model Stereo Similarity Function Based on Monogenic Signal Analysis in Poisson Scale Space

Authors: Jinjun Li; Hong Zhao; Chengying Shi; Xiang Zhou
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

A stereo similarity
function based on local multi-model monogenic
image feature descriptors (LMFD) is proposed to
match interest points and estimate disparity map
for stereo images. Local multi-model monogenic
image features include local orientation and
instantaneous phase of the gray monogenic
signal, local color phase of the color monogenic
signal, and local mean colors in the multiscale
color monogenic signal framework. The gray
monogenic signal, which is the extension of
analytic signal to gray level image using Dirac
operator and Laplace equation, consists of local
amplitude, local orientation, and instantaneous
phase of 2D image signal. The color monogenic
signal is the extension of monogenic signal to
color image based on Clifford algebras. The
local color phase can be estimated by computing
geometric product between the color monogenic
signal and a unit reference vector in RGB color
space. Experiment results on the synthetic and
natural stereo images show the performance of
the proposed approach.


---

## Body

## 1. Introduction

Shape and motion estimation from stereo images has been one of the core challenges in computer vision for decades. The robust and accurate computation of stereo depth is an important problem for many visual tasks such as machine vision, virtual reality, robot navigation, simultaneous localization and mapping, depth measurements, and 3D environment reconstruction. Most of conventional approaches, such as intensity-based or correlation-based matching, feature-based matching, and matching function optimization techniques, estimate the disparity only based on local intensity and feature between stereo images so that the results may be susceptible to level shift, scaling, rotation, and noise [1, 2].To overcome these drawbacks, we propose a new method for establishing spatial correspondences between a pair of color images. Unlike classical stereo-matching method based on brightness constancy assumption and phase congruency constraint, we match feature points and estimate disparity map between stereo images based on a new local multi-model monogenic image feature descriptors in the Color Monogenic Signal framework [3, 4]. We firstly introduce the monogenic signal [5] of 2D gray level image using Dirac operator and Laplace equation and extract local amplitude, local orientation, and instantaneous phase information in multiscale space [6]. At the same time, the 2D monogenic signal is extended to color image, and the color monogenic signal is introduced based on Clifford algebras. The local color phase is estimated by computing geometric product between the color monogenic signal and a unit reference vector in RGB color space with Clifford algebras [3, 4]. Then we focus on defining new local multi-model monogenic image feature descriptors which contain local geometric (local orientation), structure (instantaneous phase), and color (local color phase and color values) information in the Color Monogenic Signal framework. Based on the proposed image feature descriptors, a stereo similarity function between two primitives is also defined to solve stereo correspondence problem. Finally, we test the performance of the proposed approach on the synthetic and natural stereo images, and experiment results are given in detail.


## 2. Monogenic Signal Analysis in Poisson Scale Space

### 2.1. Modeling 2D Image Signal

Based on the results of Fourier theory and functional analysis, we assume that each 2D signalf∈L2(ℝ)∩L1(ℝ) can be locally modeled by a superposition of arbitrarily orientated one-dimensional cosine waves [6]:(2.1)(ps*f)(x,y)=∑v=1nav(x,y,s)cos⁡⁡(〈x,y,o̅v(x,y,s)〉+ϕv(x,y,s))
with * as the convolution operator and the orientation o̅ν(x,y,s)=[cosθv(x,y,s),sinθv(x,y,s)]T. Note that each cosine wave is determined with the same amplitude and phase information. The Poisson convolution kernel reads(2.2)ps(x,y)=s2π(s2+x2+y2)3/2.
For a certain scale space parameters s∈ℝ+, the Poisson kernel acts as a low pass filter on the original signal f. The Poisson scale space is naturally related to the generalized Hilbert transform by the Cauchy kernel. To filter a frequency interval of interest, the difference of Poisson (DoP) kernel will be used in practice:
(2.3)psf,sc(x,y)=psf(x,y)-psc(x,y)
with sc>sf>0 and sc as the coarse scale parameter and sf as the fine scale parameter. The filtered signal is defined by convolution with the difference of Poisson kernel which will be used to analyze the original with the DoP operator to consider only a small passband of the original signal spectrum. Without loss of generality the signal model in (2.1) degrades locally at the origin (x,y)=(0,0) of a local coordinate system to(2.4)fp(x,y,s)=∑v=1na(x,y,s)cos⁡ϕ(x,y,s).
In case of image analysis lines, edges, junctions, and corners can be models in this way. The signal processing task is now to determine the local amplitude a(x,y,s), the local orientation θv(x,y,s), and the local phase ϕ(x,y,s) for a certain scale space parameter s and a certain location (x,y). This problem has been already solved for one-dimensional signals by the classical analysis [7] by means of the Hilbert transform [8] and for intrinsically one-dimensional signals [9] by the two-dimensional monogenic signal by means of the generalized first-order Hilbert transforms [10].


### 2.2. The Analytic Signal

Lets: ℝ→ℝ be a real-valued signal, and let f: ℝ→ℝ2,0 [3] be a vector-valued signal such that f(x)=s(x)e2 [4]. The purpose is to construct a function fulfilling the Dirac equations whose real part is real-valued signal. It is equivalent to find the solution of a boundary value problem of the second kind (a Neumann problem):
(2.5)Δu=∂2u∂x2+∂2u∂y2=0,ify>0,e2∂u∂y=f(x),ify=0,
withD=e1(∂/∂x)+e2(∂/∂y) and  Δ=D2.The first equation in (2.5) is the 2D Laplace equation restricted to the open domain y>0. The second equation is called the boundary condition and the basis vector. e2 is coherent with the embedding of complex functions as fields (the real part is embedded as the e2-component). Using the fundamental solution of the 2D Laplace equation, the solution of the problem leads to(2.6)fA(x,y)=p1*f(x,y)+p1*h1*f(x,y),
where p1=y/π(x2+y2) is the 1D-Poisson kernel and h1=(1/πx)e12 is the Hilbert kernel. The variable y is a scale parameter, and, taking it equal to zero, the classical analytic signal can be obtained.


### 2.3. Instantaneous Phase of the Monogenic Signal

Following the previous construction of the analytic signal, Michael Felsberg and Gerald Sommer has proposed an extension to 2D signals and defined a monogenic signal which is the combination of a gray image with its Riesz transform [5]. Let s:ℝ2→ℝ be a real-valued signal and f:ℝ2→ℝ3,0 a vector-valued signal, and {ei}(i=1,2,3) is the orthonormal basis of ℝ3 such that f(x,y)=f3(x,y)e3. According to the 3D Laplace equation restricted to the open half-space z>0 and the boundary condition of the second kind, we can obtain the monogenic signal as follows:(2.7)fM(x,y,z)=p2*f(x,y)+p2*hx*f(x,y)+p2*hy*f(x,y),
where p2=z/2π(x2+y2+z2)3/2 is a 2D Poisson kernel and hR=(hx,hy)=(xe1/2π(x2+y2)3/2,ye2/2π(x2+y2)3/2) is the Riesz kernel, extension on 2D of the Hilbert kernel, z is equal to the  scale space parameter s. Now let fM(x,y,s)=fp(x,y,s)+fx(x,y,s)+fy(x,y,s) and fq(x,y,s)=[fxfy](x,y,s). In case of intrinsic dimension one signals (The intrinsic dimension expresses the number of degrees of freedom necessary to describe local structure. Constant signals without any structure are of intrinsic dimension zero (i0D), arbitrary orientated straight lines and edges are of intrinsic dimension one (i1D), and all other possible patterns such as corners and junctions are of intrinsic dimension two (i2D). In general i2D signals can only be modeled by an infinite number of superimposed i1D signals.) (i.e., n=1 in (2.1)), we can obtain:
(2.8)(fp,fx,fy)=a(cos⁡ϕ,sinϕcos⁡θ,sinϕsinθ),
where a(x,y,s), θ(x,y,s), and ϕ(x,y,s) represent the local amplitude, local orientation, and instantaneous phase, respectively:(2.9)a(x,y,s)=fp2(x,y,s)+fx2(x,y,s)+fy2(x,y,s),θ(x,y,s)=arctan(fy(x,y,s)fx(x,y,s)),ϕ(x,y,s)=arctan(fx2(x,y,s)+fy2(x,y,s)fp(x,y,s)).
However, we do not know the correct signal of the phase since it depends on the directional sense of θ(x,y,s). The best possible solution is to project it onto (cosθ,sinθ) as (2.10)ϕ̂(x,y,s)=fq(x,y,s)|fq(x,y,s)|arg⁡(fp(x,y,s)+i|fq(x,y,s)|).


### 2.4. Local Color Phase of the Color Monogenic Signal

In 2009, Demarcq et al. constructed a scale-space signal for color images seen as vectors inℝ5,0 [4]. Let s:ℝ2→ℝ3 be a real-valued signal and f:ℝ2→ℝ5,0 a vector-valued signal, and {ei}(i=1,2,3,4,5) is the orthonormal basis of ℝ5 such that f(x1,x2)=f3(x1,x2)e3+f4(x1,x2)e4+f5(x1,x2)e5. Then a color image is decomposed in the RGB space represented as the subspace spanned by {e3,e4,e5}. Now we need to find a function which is monogenic and the e3-, e4- and e5-component, of which are the components r, g, and b, respectively.(2.11)Δu=∂2u∂x12+∂2u∂x22+∂2u∂x32+∂2u∂x42+∂2u∂x52=0,e3∂u∂x3+e4∂u∂x4+e5∂u∂x5=f(x)
with D=∑i=15ei(∂/∂xi) and Δ=D2. A scale-space signal which has independent scales in each component (f3,f4,f5)=(r,g,b) can be defined by splitting the problem into three boundary value problems in ℝ5,0 as follows:
(2.12)∂2u∂x12+∂2u∂x22+∂2u∂xi2=0,ifxi>0,i=3,4,5,ei∂u∂xi=fi(x1,x2)ei,ifxi=0,i=3,4,5.Each solution of the system in (2.12) leads to monogenic functions S1, S2, S3: they satisfy the Dirac equation in each subspace Ei=span{e1,e2,ei}(i=3,4,5) and consequently the Dirac equation in ℝ5,0(DSi=0). Let fC=S1+S2+S3, then fC is still monogenic in ℝ5,0 (i.e., DfC=0) and satisfies the boundary conditions. So  the scale-space color monogenic signal can be obtained by using the Dirac operator and the Laplace equation as follows:(2.13)fc=∑i=15Aiei=∑i=35p2i*hx*fi+∑i=35p2i*hy*fi+∑i=35p2i*fiei,
where p2i=xi/(2π(x12+x22+xi2)3/2)(i=3,4,5) is a 2D Poisson kernel and hR=(hx,hy) is the Riesz kernel.  As for the analytic or monogenic signal, a color image f can be represented in terms of local amplitude and local phase. Now our proposal is to use the geometric product in order to compare two vectors in ℝ5,0.In the Clifford algebra of the Euclidean vector spaceℝn, the product of two vectors a and b,  embedded in ℝn,0, is given by [3, 4](2.14)ab=a⋅b+a∧b,
where a·b is the inner product and a∧b, the wedge product of  a and b, is a bivector. This product is usually called the geometric product of  a and b. The inner product is symmetric, and the wedge product is skew symmetric. If V=ue1+ve2+ae3+be4+ce5∈R5,0 is a chosen vector containing structure information (u,v) and color information (r,g,b), then the geometric product of fC and Vcan be given by(2.15)fCV=fC⋅V+fC∧V=〈fCV〉0+〈fCV〉2=〈fCV〉0+〈fCV〉2|〈fCV〉2||〈fCV〉2|,
where 〈·〉0 denotes the scalar part, 〈·〉2 the bivector part, and |·| the magnitude of the bivector part [6]. According to Clifford algebras, the geometric product reveals the relationship between bivectors and complex numbers [3]. This means that we can form the equivalent of a complex number, fCV=〈fCV〉0+i|〈fCV〉2|, by combining a scalar and a unit bivector. The local color phase can be computed as follows:(2.16)φc=arg(fCV)=arctan(|〈fCV〉2|〈fCV〉0).This phase describes the angular distance betweenfC and a given vector V in R5,0, that is, it gives a correlation measure between a pixel fitted with color and structure information and a vector containing chosen color and structure.


## 2.1. Modeling 2D Image Signal

Based on the results of Fourier theory and functional analysis, we assume that each 2D signalf∈L2(ℝ)∩L1(ℝ) can be locally modeled by a superposition of arbitrarily orientated one-dimensional cosine waves [6]:(2.1)(ps*f)(x,y)=∑v=1nav(x,y,s)cos⁡⁡(〈x,y,o̅v(x,y,s)〉+ϕv(x,y,s))
with * as the convolution operator and the orientation o̅ν(x,y,s)=[cosθv(x,y,s),sinθv(x,y,s)]T. Note that each cosine wave is determined with the same amplitude and phase information. The Poisson convolution kernel reads(2.2)ps(x,y)=s2π(s2+x2+y2)3/2.
For a certain scale space parameters s∈ℝ+, the Poisson kernel acts as a low pass filter on the original signal f. The Poisson scale space is naturally related to the generalized Hilbert transform by the Cauchy kernel. To filter a frequency interval of interest, the difference of Poisson (DoP) kernel will be used in practice:
(2.3)psf,sc(x,y)=psf(x,y)-psc(x,y)
with sc>sf>0 and sc as the coarse scale parameter and sf as the fine scale parameter. The filtered signal is defined by convolution with the difference of Poisson kernel which will be used to analyze the original with the DoP operator to consider only a small passband of the original signal spectrum. Without loss of generality the signal model in (2.1) degrades locally at the origin (x,y)=(0,0) of a local coordinate system to(2.4)fp(x,y,s)=∑v=1na(x,y,s)cos⁡ϕ(x,y,s).
In case of image analysis lines, edges, junctions, and corners can be models in this way. The signal processing task is now to determine the local amplitude a(x,y,s), the local orientation θv(x,y,s), and the local phase ϕ(x,y,s) for a certain scale space parameter s and a certain location (x,y). This problem has been already solved for one-dimensional signals by the classical analysis [7] by means of the Hilbert transform [8] and for intrinsically one-dimensional signals [9] by the two-dimensional monogenic signal by means of the generalized first-order Hilbert transforms [10].


## 2.2. The Analytic Signal

Lets: ℝ→ℝ be a real-valued signal, and let f: ℝ→ℝ2,0 [3] be a vector-valued signal such that f(x)=s(x)e2 [4]. The purpose is to construct a function fulfilling the Dirac equations whose real part is real-valued signal. It is equivalent to find the solution of a boundary value problem of the second kind (a Neumann problem):
(2.5)Δu=∂2u∂x2+∂2u∂y2=0,ify>0,e2∂u∂y=f(x),ify=0,
withD=e1(∂/∂x)+e2(∂/∂y) and  Δ=D2.The first equation in (2.5) is the 2D Laplace equation restricted to the open domain y>0. The second equation is called the boundary condition and the basis vector. e2 is coherent with the embedding of complex functions as fields (the real part is embedded as the e2-component). Using the fundamental solution of the 2D Laplace equation, the solution of the problem leads to(2.6)fA(x,y)=p1*f(x,y)+p1*h1*f(x,y),
where p1=y/π(x2+y2) is the 1D-Poisson kernel and h1=(1/πx)e12 is the Hilbert kernel. The variable y is a scale parameter, and, taking it equal to zero, the classical analytic signal can be obtained.


## 2.3. Instantaneous Phase of the Monogenic Signal

Following the previous construction of the analytic signal, Michael Felsberg and Gerald Sommer has proposed an extension to 2D signals and defined a monogenic signal which is the combination of a gray image with its Riesz transform [5]. Let s:ℝ2→ℝ be a real-valued signal and f:ℝ2→ℝ3,0 a vector-valued signal, and {ei}(i=1,2,3) is the orthonormal basis of ℝ3 such that f(x,y)=f3(x,y)e3. According to the 3D Laplace equation restricted to the open half-space z>0 and the boundary condition of the second kind, we can obtain the monogenic signal as follows:(2.7)fM(x,y,z)=p2*f(x,y)+p2*hx*f(x,y)+p2*hy*f(x,y),
where p2=z/2π(x2+y2+z2)3/2 is a 2D Poisson kernel and hR=(hx,hy)=(xe1/2π(x2+y2)3/2,ye2/2π(x2+y2)3/2) is the Riesz kernel, extension on 2D of the Hilbert kernel, z is equal to the  scale space parameter s. Now let fM(x,y,s)=fp(x,y,s)+fx(x,y,s)+fy(x,y,s) and fq(x,y,s)=[fxfy](x,y,s). In case of intrinsic dimension one signals (The intrinsic dimension expresses the number of degrees of freedom necessary to describe local structure. Constant signals without any structure are of intrinsic dimension zero (i0D), arbitrary orientated straight lines and edges are of intrinsic dimension one (i1D), and all other possible patterns such as corners and junctions are of intrinsic dimension two (i2D). In general i2D signals can only be modeled by an infinite number of superimposed i1D signals.) (i.e., n=1 in (2.1)), we can obtain:
(2.8)(fp,fx,fy)=a(cos⁡ϕ,sinϕcos⁡θ,sinϕsinθ),
where a(x,y,s), θ(x,y,s), and ϕ(x,y,s) represent the local amplitude, local orientation, and instantaneous phase, respectively:(2.9)a(x,y,s)=fp2(x,y,s)+fx2(x,y,s)+fy2(x,y,s),θ(x,y,s)=arctan(fy(x,y,s)fx(x,y,s)),ϕ(x,y,s)=arctan(fx2(x,y,s)+fy2(x,y,s)fp(x,y,s)).
However, we do not know the correct signal of the phase since it depends on the directional sense of θ(x,y,s). The best possible solution is to project it onto (cosθ,sinθ) as (2.10)ϕ̂(x,y,s)=fq(x,y,s)|fq(x,y,s)|arg⁡(fp(x,y,s)+i|fq(x,y,s)|).


## 2.4. Local Color Phase of the Color Monogenic Signal

In 2009, Demarcq et al. constructed a scale-space signal for color images seen as vectors inℝ5,0 [4]. Let s:ℝ2→ℝ3 be a real-valued signal and f:ℝ2→ℝ5,0 a vector-valued signal, and {ei}(i=1,2,3,4,5) is the orthonormal basis of ℝ5 such that f(x1,x2)=f3(x1,x2)e3+f4(x1,x2)e4+f5(x1,x2)e5. Then a color image is decomposed in the RGB space represented as the subspace spanned by {e3,e4,e5}. Now we need to find a function which is monogenic and the e3-, e4- and e5-component, of which are the components r, g, and b, respectively.(2.11)Δu=∂2u∂x12+∂2u∂x22+∂2u∂x32+∂2u∂x42+∂2u∂x52=0,e3∂u∂x3+e4∂u∂x4+e5∂u∂x5=f(x)
with D=∑i=15ei(∂/∂xi) and Δ=D2. A scale-space signal which has independent scales in each component (f3,f4,f5)=(r,g,b) can be defined by splitting the problem into three boundary value problems in ℝ5,0 as follows:
(2.12)∂2u∂x12+∂2u∂x22+∂2u∂xi2=0,ifxi>0,i=3,4,5,ei∂u∂xi=fi(x1,x2)ei,ifxi=0,i=3,4,5.Each solution of the system in (2.12) leads to monogenic functions S1, S2, S3: they satisfy the Dirac equation in each subspace Ei=span{e1,e2,ei}(i=3,4,5) and consequently the Dirac equation in ℝ5,0(DSi=0). Let fC=S1+S2+S3, then fC is still monogenic in ℝ5,0 (i.e., DfC=0) and satisfies the boundary conditions. So  the scale-space color monogenic signal can be obtained by using the Dirac operator and the Laplace equation as follows:(2.13)fc=∑i=15Aiei=∑i=35p2i*hx*fi+∑i=35p2i*hy*fi+∑i=35p2i*fiei,
where p2i=xi/(2π(x12+x22+xi2)3/2)(i=3,4,5) is a 2D Poisson kernel and hR=(hx,hy) is the Riesz kernel.  As for the analytic or monogenic signal, a color image f can be represented in terms of local amplitude and local phase. Now our proposal is to use the geometric product in order to compare two vectors in ℝ5,0.In the Clifford algebra of the Euclidean vector spaceℝn, the product of two vectors a and b,  embedded in ℝn,0, is given by [3, 4](2.14)ab=a⋅b+a∧b,
where a·b is the inner product and a∧b, the wedge product of  a and b, is a bivector. This product is usually called the geometric product of  a and b. The inner product is symmetric, and the wedge product is skew symmetric. If V=ue1+ve2+ae3+be4+ce5∈R5,0 is a chosen vector containing structure information (u,v) and color information (r,g,b), then the geometric product of fC and Vcan be given by(2.15)fCV=fC⋅V+fC∧V=〈fCV〉0+〈fCV〉2=〈fCV〉0+〈fCV〉2|〈fCV〉2||〈fCV〉2|,
where 〈·〉0 denotes the scalar part, 〈·〉2 the bivector part, and |·| the magnitude of the bivector part [6]. According to Clifford algebras, the geometric product reveals the relationship between bivectors and complex numbers [3]. This means that we can form the equivalent of a complex number, fCV=〈fCV〉0+i|〈fCV〉2|, by combining a scalar and a unit bivector. The local color phase can be computed as follows:(2.16)φc=arg(fCV)=arctan(|〈fCV〉2|〈fCV〉0).This phase describes the angular distance betweenfC and a given vector V in R5,0, that is, it gives a correlation measure between a pixel fitted with color and structure information and a vector containing chosen color and structure.


## 3. Local Multi-Model Monogenic Image Feature Descriptors

We make use of a visual representation [11] which is motivated by feature processing in the human visual system and define new local multi-model monogenic image descriptors which give an explicit and condensed representation of the local image signal as follows:(3.1)ε=(X,ϕ̂,φc,C).
In fact, this representation performs a considerable condensation of information in a local image patch of n×n(n∈ℝ+,n>1) pixels. The symbol X=(x,y) represents central coordinates of the local image patch, ϕ̂ is the instantaneous phase of the gray monogenic signal, φc is the local color phase, and C=(r,g,b) is the color values in RGB color space.Based on the local multi-model image feature descriptors, a stereo similarity function between two primitives is the weighted sum of squared differences of instantaneous phase, local color phase, and color vector for a pair of stereo images in the local patch:(3.2)e(εl,εr)=αdϕ(ϕ̂l,ϕ̂r)+βdφ(φl,φr)+γdC(Cl,Cr)=α|ϕ̂l-ϕ̂r|+β|φcl-φcr|+γ∑i∈{r,g,b}(cli-cri),
where l,r represent the left and right images, respectively, dϕ∈[0,π) is the distance measurement of instantaneous phase ϕ∈[-π,π), dφ∈[0,π] is the distance measurement of local color phase φc∈[0,π], and dc∈[0,3] is the distance measurement of color vector with C∈[0,1]×[0,1]×[0,1] and cR+cG+cB=1 in RGB color space. The symbols α,β,γ are weighted coefficients with α,β∈[0,1] and γ∈[0,0.5]. In order to achieve a better coherence with the real scene, we use an adaptive support-weight technique in the local patch [12]. The support weight for each pixel in the window is calculated based on the Gestalt Principles, which state that the grouping of pixels should be based on spatial proximity and chromatic similarity. The original formula proposed is given as follows [13]:(3.3)w(x,y,m,n)=exp⁡(-(Δcx,y,m,nγc+Δqx,y,m,nγg)),(3.4)E(xl,yl,xr,yr)=∑m,n∈[-p,p]wl(xl,yl,m,n)⋅wr(xr,yr,m,n)⋅e(εl,εr)∑m,n∈[-p,p]w(x,y,m,n)⋅wr(xr,yr,m,n),
where (xl,yl) and (xr,yr) are the pixel of interest in left and right images, (m,n) is the pixel offset within the local window, w(x,y,m,n) represents the weight of neighboring pixel (x+m,y+n), Δcx,y,m,n is the color difference between pixel (x,y) and   (x+m,y+n), Δqx,y,m,n is the Euclidean distance between pixel (x,y) and (x+m,y+n), γc and γg are user defined parameters and, E(xl,yl,xr,yr) is the aggregated cost between pixel (xl,yl) in the left image and (xr,yr) in the right image.In the proposed feature descriptors, there are several merits. Firstly, the instantaneous phase contains local orientation (or geometric information) and information about contrast transition so that it describes an intrinsically one-dimensional structure in a grey level image, that is, an image structure that is dominated by one orientation. Examples of different contrast transitions are a dark/bright (bright/dark) edge or a bright (dark) line on dark (bright) background. Of course, there is a continuum between these different grey level structures. The instantaneous phase as an additional feature allows us to take this information into account (as one parameter in addition to orientation) in a compact way [11]. Secondly, local color phase describes the angular distance between color monogenic signal fc of color image and a given color vector V in R5,0 and gives a correlation measure between a pixel fitted with color and structure information and a vector containing chosen color and structure [14]. Finally, the color vector indicates the mean color structure of local image because color is also an important cue to improve stereo matching. Because the stereo similarity function with the local multi-model feature descriptors contains local geometric, structure, and color information, it is much more robust against noise and brightness change than others in feature matching and 3D reconstruction.


## 4. Experimentation Results

Once the similarity function is given, minimization process can be performed to find the optimal disparity. In order to reduce noise sensitivity and simultaneously achieve higher efficiency, the multiscale space and a winner-take-all technique are employed to optimize the disparity map for stereo matching. We test the performance of the proposed local multi-model feature descriptors and similarity function on the synthetic and natural stereo images. On the first step, we estimate and compare the disparity maps of a pair of synthetic images using each of three distance measurement and the proposed stereo similarity function. On the second step, we reconstruct 3D shape and appearance of natural object and scene combing the proposed stereo similarity function and multiview stereo technique [15].


### 4.1. Disparity Estimation Experiment

In the first experiment, we chose a pair of color images (cloth1) from the website [16]; the left image, right image, and ground disparity map of the stereo pair are showed in Figures 1(a)–1(c), respectively.Stereo pair of Cloth1. (a) Left image, (b) right image, and (c) ground disparity map.
(a)(b)(c)Firstly, the left and right color images are converted to gray level image. The gray monogenic signals withs=3 scales are computed, and the instantaneous phases are estimated for both gray images in Poisson scale space. Figure 2 show three instantaneous phase maps of the left image. At the same time, the scale-space color monogenic signals with s=3 scales are computed, and a chosen reference vector V=(e3+e4+e5)/3 with local geometric structure (u,v)=(0,0) and a unit vector (a,b,c)=(1,1,1)/3 in RGB color space are used to estimate the local color phase for both images [14]. Figure 3 also shows three color phase maps of the left image.The instantaneous phase maps. (a) First scale, (b) second scale, and (c) third scale.
(a)(b)(c)Local-color-phase maps. (a) First scale, (b) second scale, and (c) third scale.
(a)(b)(c)Secondly, the proposed stereo similarity function based on local multi-model monogenic image descriptors is employed to compute the weighted sum of squared differences of instantaneous phase, local color phase, and color vector. The scale-space color monogenic framework [14] and a local optimization technique [17] are employed to optimize the disparity map according to the adaptive support-weight cost aggregation as in (3.4). In the scale space, the disparity field at the coarser scale is used as starting guesses for estimation at the next finer scale, then a subpixel disparity can be obtained. For example, we estimate the disparity map at the third scale using the cost function (3.2) and (3.4). The disparities at the third scale directly subtending the second scale under scrutiny are all used as candidate starting guesses. The candidate leading to best match is accepted for the regularization step. Estimation proceeds in this fashion, decrementing scale from coarser to finer, with matching following by regularization, until the finest level of detail is reached. Figure 4 shows the dense disparity map based on the proposed algorithm. Figures 4(a)–4(d) show four estimated disparity maps with the different weighted coefficients (α,β,γ).Experiment results for disparity estimation with the different weighted coefficients. (a)α=0,β=0,γ=0.5; (b) α=1,β=0,γ=0; (c) α=0,β=1,γ=0; (d) α=1,β=1,γ=0.5.
(a)(b)(c)(d)Thirdly, a comparison is performed to further validate the claims about the performance of the proposed algorithm. This comparison is performed between the proposed algorithm and a number of selected algorithms from the literature, which are Sum of Absolute Differences (SAD) [18], Sum of Squared Differences (SSD) [18], Graph Cuts (GC) [18], Discrete Wavelet (DWT) [19], Complex Discrete Wavelet (CWT) [19], and Quaternion Wavelet (QWT) [20]. To create a better understanding of the comparison, statistic root mean square error (RMSE) and percent-age of bad disparities (PBD) are calculated as [19](4.1)RMSE=1N∑i,j|dE(i,j)-dG(i,j)|2,PBD=1N∑i,j|dE(i,j)-dG(i,j)|>ζ,
where dE and dG are the estimated and ground truth disparity maps and N is the total number of pixels in an image, whereas ζ represents the disparity error tolerance. The statistics RMSE and PBD related to all of the above algorithms is presented in Table 1. As it can be seen in Table 1, the proposed algorithms with the weighted coefficients α=1, β=1, γ=0.5 obtain better results (RMSE = 1.82 and PDB = 0.09) than others.Table 1
The comparison of the estimated disparity with the different algorithms.
AlgorithmsThe selected algorithmsThe proposed algorithmsSADSSDGCCWTDWTQWTα=0α=1α=0α=1β=0β=0β=1β=1γ=0.5γ=0γ=0γ=0.5RMSE4.233.782.252.693.172.533.242.611.981.82PBD0.370.310.150.190.290.180.300.200.130.09


### 4.2. The Robustness against Noise and Brightness

Now we aim to investigate the robustness of the proposed approach against noise and brightness change. In a first step, Gaussian noises withσ=0.05,0.10,0.15,0.20,0.25,0.30 are added to both color images, respectively, and the disparity maps are estimated using the proposed method with four different weighted coefficients. The statistic root mean square error (RMSE) and percentage of bad disparities (PBD) for each disparity map with noise also are calculated using (4.1). In a second step, the brightness values with I=-30,-20,-10,10,20,30 are added to the right image. We also compute the disparity maps, RMSE and PBD, with four different weighted coefficients.  The experiment results in Table 2 show that the proposed approach with =1,β=1, γ=0.5 are insusceptible to noise and brightness change.Table 2
Experiment results for testing robustness against noise and brightness change.
Noise/brightness changeα=0,β=0,γ=0.5α=1,β=0,γ=0α=0,β=1,γ=0α=1,β=1,γ=0.5RMSEPBDRMSEPBDRMSEPBDRMSEPBDNone3.240.302.610.201.980.131.820.09σ=0.053.460.342.730.221.990.141.890.10σ=0.103.730.372.890.242.130.151.960.11σ=0.153.980.412.990.262.180.162.040.12σ=0.204.410.463.170.282.290.182.120.13σ=0.255.090.513.290.312.410.202.190.15σ=0.306.120.593.500.352.530.222.260.18I=-306.470.622.860.282.300.192.070.14I=-205.350.542.670.252.120.161.960.12I=-104.050.412.430.231.970.141.890.10I=104.020.392.390.221.960.131.860.10I=205.100.492.460.242.140.151.940.11I=306.590.612.720.272.270.182.090.14


### 4.3. 3D Reconstruction Experimentation

In the last experiment, the proposed approach is used to reconstruct a set of oriented points for 3D natural scenes by multiview stereopsis. We firstly captured stereo images with the size of1024×768 around 3D object and scene (a static color paper cup and a deforming color texture paper) under the condition of natural light by two color cameras (AVT Guppy F146C). We also calibrate the intrinsic and external parameters for all of cameras and images and estimated the corresponding epipolar lines among images. Feature points for all images are detected by using Speeded-Up Robust Features (SURF) with blob response threshold of 1000 [21] and are matched by using the corresponding epipolar constraints and the proposed similarity function with α=1, β=1, γ=0.5 in (3.2) and (3.4). According to these init matched feature points, we reconstruct a set of 3D sparse oriented points for object and scene similar to multiview stereopsis in [15]. Then we expand and filter it to a set of robust 3D dense oriented points with each image cell of 2×2 pixels and reconstruct. We reconstruct three-dimensional surface of the object. For the static color paper cup, we capture 16 images and reconstruct the 3D dense oriented points and shape. In Figure 5, there are the captured images, the 3D dense oriented points, and three-dimensional surface at the 1st, 7th, and 13th view, respectively.  For the deforming color texture paper, we capture the stereo sequence using two cameras at the rate of 7 frames each second and reconstruct the time-varying 3D shape and motion. Figure 6 shows the captured image of one camera and the time-varying 3D shape at the 0th, 20th, 40th, 80th, and 100th frames.The 3D dense oriented points reconstruction for the color paper cup. (a)–(c) The captured images, (d)–(f) the 3D dense oriented points, and (g)–(i) the three-dimensional surface at the 1st, 7th, and 13th view.
(a)(b)(c)(d)(e)(f)(g)(h)(i)The time-varying shape and motion estimation for the color texture paper. (a)–(f) The captured images and (g)–(l) the three-dimensional surface at the 0th, 20th, 40th, 60th, 80th, and 100th frames.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)To further validate the performance of the proposed method, a new set of 3D dense patches for the static color paper cup also is reconstructed by the SSD-based photometric discrepancy function in [15]. And, for both algorithms, we calculated the mean feature points (MFP) and percentage of bad disparities (PBD) of each image, the number of 3D oriented points at the init matching (MAT_3D), expanding (EXP_3D) and filtering (FIL_3D) stages, and the percentage of  image cells not to be reconstructed (PNR). Statistical results are shown in Table 3.Table 3
Experiment results for matching features and reconstructed 3D oriented points.
AlgorithmMFPPBDMAT_3DEXP_3DFIL_3DPNRThe proposed method16320.184699102417993367.26%SSD-based16320.3835788796816916.32%


## 4.1. Disparity Estimation Experiment

In the first experiment, we chose a pair of color images (cloth1) from the website [16]; the left image, right image, and ground disparity map of the stereo pair are showed in Figures 1(a)–1(c), respectively.Stereo pair of Cloth1. (a) Left image, (b) right image, and (c) ground disparity map.
(a)(b)(c)Firstly, the left and right color images are converted to gray level image. The gray monogenic signals withs=3 scales are computed, and the instantaneous phases are estimated for both gray images in Poisson scale space. Figure 2 show three instantaneous phase maps of the left image. At the same time, the scale-space color monogenic signals with s=3 scales are computed, and a chosen reference vector V=(e3+e4+e5)/3 with local geometric structure (u,v)=(0,0) and a unit vector (a,b,c)=(1,1,1)/3 in RGB color space are used to estimate the local color phase for both images [14]. Figure 3 also shows three color phase maps of the left image.The instantaneous phase maps. (a) First scale, (b) second scale, and (c) third scale.
(a)(b)(c)Local-color-phase maps. (a) First scale, (b) second scale, and (c) third scale.
(a)(b)(c)Secondly, the proposed stereo similarity function based on local multi-model monogenic image descriptors is employed to compute the weighted sum of squared differences of instantaneous phase, local color phase, and color vector. The scale-space color monogenic framework [14] and a local optimization technique [17] are employed to optimize the disparity map according to the adaptive support-weight cost aggregation as in (3.4). In the scale space, the disparity field at the coarser scale is used as starting guesses for estimation at the next finer scale, then a subpixel disparity can be obtained. For example, we estimate the disparity map at the third scale using the cost function (3.2) and (3.4). The disparities at the third scale directly subtending the second scale under scrutiny are all used as candidate starting guesses. The candidate leading to best match is accepted for the regularization step. Estimation proceeds in this fashion, decrementing scale from coarser to finer, with matching following by regularization, until the finest level of detail is reached. Figure 4 shows the dense disparity map based on the proposed algorithm. Figures 4(a)–4(d) show four estimated disparity maps with the different weighted coefficients (α,β,γ).Experiment results for disparity estimation with the different weighted coefficients. (a)α=0,β=0,γ=0.5; (b) α=1,β=0,γ=0; (c) α=0,β=1,γ=0; (d) α=1,β=1,γ=0.5.
(a)(b)(c)(d)Thirdly, a comparison is performed to further validate the claims about the performance of the proposed algorithm. This comparison is performed between the proposed algorithm and a number of selected algorithms from the literature, which are Sum of Absolute Differences (SAD) [18], Sum of Squared Differences (SSD) [18], Graph Cuts (GC) [18], Discrete Wavelet (DWT) [19], Complex Discrete Wavelet (CWT) [19], and Quaternion Wavelet (QWT) [20]. To create a better understanding of the comparison, statistic root mean square error (RMSE) and percent-age of bad disparities (PBD) are calculated as [19](4.1)RMSE=1N∑i,j|dE(i,j)-dG(i,j)|2,PBD=1N∑i,j|dE(i,j)-dG(i,j)|>ζ,
where dE and dG are the estimated and ground truth disparity maps and N is the total number of pixels in an image, whereas ζ represents the disparity error tolerance. The statistics RMSE and PBD related to all of the above algorithms is presented in Table 1. As it can be seen in Table 1, the proposed algorithms with the weighted coefficients α=1, β=1, γ=0.5 obtain better results (RMSE = 1.82 and PDB = 0.09) than others.Table 1
The comparison of the estimated disparity with the different algorithms.
AlgorithmsThe selected algorithmsThe proposed algorithmsSADSSDGCCWTDWTQWTα=0α=1α=0α=1β=0β=0β=1β=1γ=0.5γ=0γ=0γ=0.5RMSE4.233.782.252.693.172.533.242.611.981.82PBD0.370.310.150.190.290.180.300.200.130.09


## 4.2. The Robustness against Noise and Brightness

Now we aim to investigate the robustness of the proposed approach against noise and brightness change. In a first step, Gaussian noises withσ=0.05,0.10,0.15,0.20,0.25,0.30 are added to both color images, respectively, and the disparity maps are estimated using the proposed method with four different weighted coefficients. The statistic root mean square error (RMSE) and percentage of bad disparities (PBD) for each disparity map with noise also are calculated using (4.1). In a second step, the brightness values with I=-30,-20,-10,10,20,30 are added to the right image. We also compute the disparity maps, RMSE and PBD, with four different weighted coefficients.  The experiment results in Table 2 show that the proposed approach with =1,β=1, γ=0.5 are insusceptible to noise and brightness change.Table 2
Experiment results for testing robustness against noise and brightness change.
Noise/brightness changeα=0,β=0,γ=0.5α=1,β=0,γ=0α=0,β=1,γ=0α=1,β=1,γ=0.5RMSEPBDRMSEPBDRMSEPBDRMSEPBDNone3.240.302.610.201.980.131.820.09σ=0.053.460.342.730.221.990.141.890.10σ=0.103.730.372.890.242.130.151.960.11σ=0.153.980.412.990.262.180.162.040.12σ=0.204.410.463.170.282.290.182.120.13σ=0.255.090.513.290.312.410.202.190.15σ=0.306.120.593.500.352.530.222.260.18I=-306.470.622.860.282.300.192.070.14I=-205.350.542.670.252.120.161.960.12I=-104.050.412.430.231.970.141.890.10I=104.020.392.390.221.960.131.860.10I=205.100.492.460.242.140.151.940.11I=306.590.612.720.272.270.182.090.14


## 4.3. 3D Reconstruction Experimentation

In the last experiment, the proposed approach is used to reconstruct a set of oriented points for 3D natural scenes by multiview stereopsis. We firstly captured stereo images with the size of1024×768 around 3D object and scene (a static color paper cup and a deforming color texture paper) under the condition of natural light by two color cameras (AVT Guppy F146C). We also calibrate the intrinsic and external parameters for all of cameras and images and estimated the corresponding epipolar lines among images. Feature points for all images are detected by using Speeded-Up Robust Features (SURF) with blob response threshold of 1000 [21] and are matched by using the corresponding epipolar constraints and the proposed similarity function with α=1, β=1, γ=0.5 in (3.2) and (3.4). According to these init matched feature points, we reconstruct a set of 3D sparse oriented points for object and scene similar to multiview stereopsis in [15]. Then we expand and filter it to a set of robust 3D dense oriented points with each image cell of 2×2 pixels and reconstruct. We reconstruct three-dimensional surface of the object. For the static color paper cup, we capture 16 images and reconstruct the 3D dense oriented points and shape. In Figure 5, there are the captured images, the 3D dense oriented points, and three-dimensional surface at the 1st, 7th, and 13th view, respectively.  For the deforming color texture paper, we capture the stereo sequence using two cameras at the rate of 7 frames each second and reconstruct the time-varying 3D shape and motion. Figure 6 shows the captured image of one camera and the time-varying 3D shape at the 0th, 20th, 40th, 80th, and 100th frames.The 3D dense oriented points reconstruction for the color paper cup. (a)–(c) The captured images, (d)–(f) the 3D dense oriented points, and (g)–(i) the three-dimensional surface at the 1st, 7th, and 13th view.
(a)(b)(c)(d)(e)(f)(g)(h)(i)The time-varying shape and motion estimation for the color texture paper. (a)–(f) The captured images and (g)–(l) the three-dimensional surface at the 0th, 20th, 40th, 60th, 80th, and 100th frames.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)To further validate the performance of the proposed method, a new set of 3D dense patches for the static color paper cup also is reconstructed by the SSD-based photometric discrepancy function in [15]. And, for both algorithms, we calculated the mean feature points (MFP) and percentage of bad disparities (PBD) of each image, the number of 3D oriented points at the init matching (MAT_3D), expanding (EXP_3D) and filtering (FIL_3D) stages, and the percentage of  image cells not to be reconstructed (PNR). Statistical results are shown in Table 3.Table 3
Experiment results for matching features and reconstructed 3D oriented points.
AlgorithmMFPPBDMAT_3DEXP_3DFIL_3DPNRThe proposed method16320.184699102417993367.26%SSD-based16320.3835788796816916.32%


## 5. Conclusions

In the paper, we propose a stereo similarity function based on local multi-model monogenic image feature descriptors to solve stereo-matching problem. Local multi-model monogenic image features include local orientation and instantaneous phase of the gray monogenic signal, local color phase of the color monogenic signal, and local mean colors in the multiscale color monogenic signal framework. The gray monogenic signal, which is the extension of analytic signal to gray level image using Dirac operator and Laplace equation, consists of local amplitude, local orientation, and instantaneous phase of 2D image signal. The color monogenic signal is the extension of monogenic signal to color image based on Clifford algebras. The local color phase can be estimated by computing geometric product between the color monogenic signal and a unit reference vector in RGB color space. Experiment results on the synthetic and natural stereo images show the performance of the proposed approach. But there is a shortcoming that the proposed method needs much more run time and storage space. So we will be devoted to improving the efficiency and evaluating the proposed method in the future work.



---

*Source: 202653-2011-05-16.xml*