Yves Bourgault

Charles Pierre

Comparing the bidomain

Keywords: February 9, 2011 Electro-cardiology, bidomain and monodomain models, excitation process and activation time, reaction diffusion equations, numerical simulations, anisotropic diffusion. 1

come

Introduction

The bidomain model [START_REF] Tung | A bidomain model for describing ischemic myocardial D-D properties[END_REF][START_REF] Krassowska | Homogenization of syncytial tissues[END_REF][START_REF] Ambrosio | On the asymptotic behaviour of anisotropic energies arising in the cardiac bidomain model[END_REF][START_REF] Colli-Franzone | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF][START_REF] Veneroni | Reaction-diffusion systems for the microscopic cellular model of the cardiac electric field[END_REF] is currently considered as the most accurate and physiologically founded description for the electrical cardiac behaviour and is widely used to simulate action potential spreading in the myocardium as well as electrocardiograms. Its mathematical formulation reads a system of two parabolic reaction diffusion equations, or equivalently one parabolic reaction diffusion equation coupled with one elliptic equation. This system is coupled with an ODE system describing cell membrane kinetics. The mathematical properties of the bidomain model are quite delicate to study. In [START_REF] Colli-Franzone | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF] the bidomain model is studied under the form of a degenerate system of parabolic reaction diffusion equations. It has been shown in [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF] that it can also be reformulated into one parabolic semi-linear PDE but including non locality in space. These structural properties (degeneracy in one hand, non locality on the other) bring numerical difficulties. One ill-conditioned linear system inversion is required per time step. Moreover cardiac action potential involving fast space and time potential variations, fine space and time grids must be considered. For these two reasons simulating the cardiac electrical activity with the bidomain model has a very high cost, and many efforts towards the reduction of this cost have been made, e.g. [START_REF] Colli Franzone | A parallel solver for reaction-diffusion systems in computational electrocardiology[END_REF][START_REF] Ethier | Semi-implicit time discretization schemes for the bidomain model[END_REF][START_REF] Colli Franzone | Adaptivity in space and time for reaction-diffusion systems in electrocardiology[END_REF]. The monodomain model is a simplification of the bidomain model reading a single parabolic reaction diffusion equation (still coupled with the same ODE system modelling cell membrane). Although this simplification has no mathematical general justification, and although the monodomain model lacks physiological foundation, it is commonly used in electro-cardiology: firstly because it obviously lead to much lower computational efforts than the bidomain model (as analysed in [START_REF] Sundnes | On the computational complexity of the bidomain and the monodomain models of electrophysiology[END_REF]). The second reason motivating the interest for the monodomain model is that, as an approximation of the bidomain model, it may serve to improve numerical scheme efficiency for the bidomain model [START_REF] Gerardo-Giorda | Optimized schwarz coupling of bidomain and monodomain models in electrocardiology[END_REF] or to build powerful preconditioners [START_REF] Gerardo-Giorda | A model-based block-triangular preconditioner for the bidomain system in electrocardiology[END_REF][START_REF] Pierre | Preconditioning the coupled heart and torso bidomain model with an almost linear complexity[END_REF]. The quality of the monodomain model approximation cam moreover be numerically optimised as developed in [START_REF] Nielsen | Optimal monodomain approximations of the bidomain equations[END_REF].

In comparison with the amount of papers dealing with the bidomain and the monodomain models, quite few studies measuring the discrepancy between these two models are available. In [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models[END_REF][START_REF] Mc | Simulation of qrst integral maps with a membrane-based computer heart model employing parallel processing[END_REF], a precise approximation of the bidomain model, namely the adapted monodomain model (or modified monodomain), is stated and compared with the bidomain model. In [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models[END_REF], Colli-Franzone et al. compared the two models on an academic 3D test case, a small slab of cardiac tissue, including anisotropy (orthotropic or axysymmetric). Their comparison is based on activation time, recovery time and action potential duration measurements. They observed a strong qualitative agreement between the two models and a noticeable quantitative difference (5 to 10 % of relative error on the activation time e.g.). In [START_REF] Mc | Simulation of qrst integral maps with a membrane-based computer heart model employing parallel processing[END_REF], Potse et al. compared the two models on a complete 3D human hear. The discrepancy between the models is here reported to be small but not precisely measured (in terms of activation time, recovery time and epicardial potential measurements). Although a strong quantitative agreement for the epicardial potential depicted in the paper is clear, quantitative differences are again visible.

Meanwhile, Potse et al. in [START_REF] Mc | Simulation of qrst integral maps with a membrane-based computer heart model employing parallel processing[END_REF] brought to the fore a fundamental question: "when comparing the bidomain and monodomain models, how to distinguish in the numerical results the amount of error due to the model choice from the amount of error caused by discretisation ?" Precisely in [START_REF] Mc | Simulation of qrst integral maps with a membrane-based computer heart model employing parallel processing[END_REF], potential wave velocities have been measured for both models on two meshes of 0.1 and 0.2 mm resolution. On the finest grid the relative error on the longitudinal velocity predicted by the two models is of 2.5 %. The same longitudinal velocity predicted by the bidomain model on the two grids differ of 5.0 % (also of relative error). The error caused by discretisation then is twice larger than the measured discrepancy between the models. From the methodological point of view, no conclusion can be drawn in such a situation concerning a quantification of the discrepancy between the bidomain and monodomain models.

This example also illustrate that even on very fine meshes the discretisation error might be as high as the model discrepancy. Appreciating this error is not easy. From the theoretical point of view, a posteriori error estimators for the quite complex bidomain model (involving enhanced reaction terms in practice) and especially designed for physiologically oriented criterion such as activation time or APD (action potential duration) have not be developed yet. From the numerical point of view, to the author knowledge very few studies intending to evaluate the discretisation errors through convergence analysis are available. Such a convergence study has been presented in [START_REF] Andreianov | Convergence of DDFV schemes for the bidomain cardiac model[END_REF] concerning the approximation of the bidomain model by Discrete Duality Finite Volume . Though the model setting was rather academical (involving in particular a very simplified cell membrane description), it illustrates that the convergence in terms of activation time is quite slow.

In this paper we propose a numerical methodology to evaluate the discrepancy between the bidomain and adapted monodomain models at the continuous level. Because of its particular physiological importance the comparison criterion is here set to the activation time, though the method used could also be applied to other criteria. The considered method is a convergence analysis for the two models allowing to evaluate the discretisation error and thus to get rid of it when evaluating the model discrepancy. The two following purposes of same practical importance in biomedical engineering are then addressed in this paper:

(1) measuring the discrepancy between the bidomain and monodomain models, (2) evaluating the discretisation error and providing necessary grid resolution in order to obtain accurate simulations of cardiac spreading of excitation.

Two different test cases are considered in dimension 2, both involving complex ionic dynamics on the cell membrane and tissue anisotropy. A first test for an academic geometry on the unit square. A second test case on a 2D cut of a segmented human heart. For both test cases the Luo and Rudy class II model of cell membrane [START_REF] Luo | A Dynamic Model of the Cardiac Ventricular Action Potential I. Simulations of Ionic Currents and Concentration Changes[END_REF] (for mammalian ventricular cells) will be considered. The paper is organised as follows. The bidomain and the adapted monodomain models are formulated in Sec. 2. In Sec. 3 is presented the methodology including the numerical implementation of the models and the two tests case presentation. The simulation results are given in Sec. 4. Results are discussed and conclusions are drawn in Sec. 5.

Models

The heart of a living organism is assumed to occupy a fixed domain Ω that is a bounded open subset of R d , d = 2, 3. At the macroscopic scale the cardiac tissue is considered as the superimposed of the intra-cellular (i) and extra-cellular (e) media. The bidomain and adapted monodomain models presented here describe the heart electrical activity at this scale.

For a derivation of the bidomain model from physiological considerations at the microscopic scale we refer to [START_REF] Tung | A bidomain model for describing ischemic myocardial D-D properties[END_REF][START_REF] Krassowska | Homogenization of syncytial tissues[END_REF][START_REF] Ambrosio | On the asymptotic behaviour of anisotropic energies arising in the cardiac bidomain model[END_REF][START_REF] Colli-Franzone | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF][START_REF] Veneroni | Reaction-diffusion systems for the microscopic cellular model of the cardiac electric field[END_REF].

These two models involve two electrical potentials: the intra-and extra-cellular potentials u i , u e : Q → R, where Q denotes the time-space cylinder (0, T ) × Ω. Their difference is referred to as the transmembrane potential : v := u i -u e : Q → R. The heart tissue has a fibrous organization into muscular fibres. This causes anisotropy for the electrical conductivity described by two conductivity tensors σ i (x) and σ e (x) at point x ∈ Ω. Introducing the conductivities g l i,e , g t i,e longitudinal and transverse to the fibres, these tensors read:

σ i (x) = Diag(g l i , g t i )
, σ e (x) = Diag(g l e , g t e ), (1) in a moving system of coordinates whose principal direction is given by the fibre orientation at point x. The fibres are moreover assumed to be tangent to the heart boundary ∂Ω.

In this paper, we consider the case of an isolated heart (the interaction with the surrounding tissues is neglected). This insulation assumption reads a zero flux boundary condition on u i and u e physically meaning that no current flows out of the heart. The fibres being tangent to ∂Ω, this condition is equivalent with:

∇u i • n = ∇u e • n = 0 on ∂Ω, (2)
and so, ∇v • n = 0 on ∂Ω.

2.1. Bidomain model. The bidomain model reads the three following equations, for (t, x) ∈ Q:

(4)      div((σ i (x) + σ e (x))∇u e ) = -div(σ i (x)∇v), χ (c∂ t v + I ion (v, w) -I st (t, x)) = div(σ i (x)∇(u e + v)), ∂ t w = g(v, w).
In equation 2, c denotes the membrane surface capacitance, χ is the rate of cell membrane surface per unit volume (homogenization parameter), I st : Q → R is the stimulation current (source term). I ion (v, w) (reaction term) denotes the surface ionic current distribution on the membrane: its evolution is controlled by the gating variable w : Q → R N via the ODE system in line 3. The definitions of I ion and of g are fixed by the chosen ionic model in 2.4. Equations in (4) are closed by the boundary condition [START_REF] Andreianov | Convergence of DDFV schemes for the bidomain cardiac model[END_REF] and by an initial datum imposed on v and w:

(5) v(0, x) = v 0 (x), w(0, x) = w 0 (x), x ∈ Ω.
Clearly, the bidomain model equations ( 4), ( 2) and ( 5) are invariant under the simultaneous change of u e ,, u i into ue + k, u i + k for k ∈ R: we therefore impose the normalization condition [START_REF] Clayton | Models of cardiac tissue electrophysiology: Progress, challenges and open questions[END_REF] Ω u e (t, •)dx = 0. 2.2. Adapted monodomain model. The adapted monodomain model is a simplification of the bidomain model where the transmembrane potential is simply defined by a parabolic reaction-diffusion equation, for (t, x) ∈ Q: [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF] χ (c∂ t v + I ion (v, w) -I st (t, x)) = div(σ m (x)∇v),

∂ t w = g(v, w),
with boundary condition (3) and initial condition [START_REF] Cai | The finite volume element method for diffusion equations on general triangulations[END_REF]. The conductivity tensor σ m is defined as the harmonic mean between σ i and σ e , ( 8)

σ -1 m (x) := σ -1 i (x) + σ -1 e (x)
. Again, since the fibres are tangent to the domain boundary, then σ m (x)n and n have the same direction on ∂Ω, so that (3) precisely is the classical homogeneous Neumann boundary condition on ∂Ω associated with the operator div(σ m ∇•).

The potential v : Q → R being defined by equations ( 7) (3) ( 5), u e can be recovered using the second equation in (4), this will not be discussed here.

Comments.

In the adapted monodomain model framework, v is independent from u e ; whereas for the bidomain model, v and u e are strongly coupled and in general none of these two quantities can be computed independently from the other one. This illustrates a complete difference of nature between these two models which we briefly discuss here.

In [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF] it has been showed that the bidomain model can be reformulated in terms of v only as follows:

χ (c∂ t v + I ion (v, w) -I st (t, x)) = Av, ∂ t w = g(v, w),
where A denotes the harmonic mean between the two elliptic operators div(σ i ∇•) and div(σ e ∇•). This operator is non-local in general, in the following sense: Au(x) is not determined by the values of u in a neighbourhood of the point x. Precisely, the definition of Au requires an elliptic problem inversion. The strong coupling between v and u e precisely relies on this fact. However in some particular cases, namely in dimension 1 or in case of equal anisotropy ratio (i.e. if σ e (x) = k σ i (x) for some constant k ∈ R), one has the equality A = div(σ m ∇•) for the tensor σ m in [START_REF] Colli Franzone | Adaptivity in space and time for reaction-diffusion systems in electrocardiology[END_REF]. In these cases, the two models coincide: they predict exactly the same transmembrane potential v. These cases are exceptional. In general A is non local and the two models do not match. The adapted monodomain model thus is based on the heuristic approximation A div(σ m ∇•). Its motivation is to benefit from the numerical facilities provided by this approximation. On the contrary of the bidomain model, no physiological interpretation for [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF] has to be sought. We point out that the equal anisotropy ratio condition do not fit with the experimental data. The conductivities considered in this paper (introduced in the following subsection and given in Tab. 1) induce anisotropy ratios between the longitudinal and transverse directions of 9.0 and 2.0 for the intra and extra-cellular media respectively. 2.4. Model settings. In ( 4) and ( 7), the reaction terms I ion and g are fixed by choosing an ionic model: such a model describes the ion transfer across the cell membrane due to the cell metabolism. The gating variable w being here aimed to characterise the state of the cell membrane. We consider in this paper the Luo and Rudy II model [START_REF] Luo | A Dynamic Model of the Cardiac Ventricular Action Potential I. Simulations of Ionic Currents and Concentration Changes[END_REF] that has been developed for mammalian ventricular cells.

Model parameters

Value Unit Cell membrane surface-to-volume ratio

χ ∈ [1000, 2000] [cm -1 ] Membrane surface capacitance c = 1.0 [µ F/cm 2 ] Longitudinal intra-cellular conductivity g l i = 1.741 [mS/cm] Transverse intra-cellular conductivity g t i = 0.1934 [mS/cm] Longitudinal extra-cellular conductivity g l e = 3.

[mS/cm] Transverse extra-cellular conductivity

g t e = 1.970 [mS/cm] Table 1. Bidomain and adapted monodomain model parameters All along this paper, the reference length and time scales are respectively set to the cm and ms. The macroscopic cardiac tissue parameters are displayed in Tab. 1. In this table F and S stand for the Faraday and the Siemens units. Cardiac tissue conductivities are subject to strong individual variabilities [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF][START_REF] Ms | Extracellular Discontinuities in Cardiac Muscle Evidence for Capillary Effects on the Action Potential Foot[END_REF]: ranges for these values are available in the review paper [START_REF] Clayton | Models of cardiac tissue electrophysiology: Progress, challenges and open questions[END_REF]. In the present paper, conductivities are taken from [START_REF] Guyader | Extracellular measurement of anisotropic bidomain myocardial conductivities. I. theoretical analysis[END_REF]. The longitudinal to transverse conductivity ratio is of 9.0 for the intra-cellular media and of 2.0 for the extra-cellular media: thus clearly avoiding here the equal anisotropy ratio where bidomain and monodomain models do coincide, see Sec. 2.3. These conductivities induce an axial to transverse velocity ratio of 2.6. Setting χ to 1800 , axial and transverse velocities are of 0.5 and 0.19 m/s respectively, which data are in agreement with [START_REF] Ms | Extracellular Discontinuities in Cardiac Muscle Evidence for Capillary Effects on the Action Potential Foot[END_REF]. A wide range of values for the cell membrane surface to volume ratio χ also is available in the literature. In [START_REF] Spach | Electrophysiological Effects of Remodeling Cardiac Gap Junctions and Cell Size[END_REF] the values of 3 300 and 6 600 cm -1 is reported for a human adult and neonatal heart respectively. A range of 2400 to 8400 cm -1 is provided in [START_REF] Clayton | Models of cardiac tissue electrophysiology: Progress, challenges and open questions[END_REF]. Beyond physiological questioning about valuing χ we only focus In this paper on the impact of its value for numerical computations. To measure this impact, a range of χ ∈ [1000, 2000] has been adopted here.

Methodology

3.1. Numerical methods. The following weak formulation of the bidomain model is considered, ∀ψ ∈ H 1 (Ω):

Ω (σ i + σ e )∇u e • ∇ψdx + Ω σ i ∇v • ∇ψdx = 0, ( 9
)
χc∂ t Ω vψdx + χ Ω (I ion (v, w) -I st (x, t))ψdx = Ω σ i ∇(u e + v) • ∇ψdx. (10)
To this system is added the normalisation condition [START_REF] Clayton | Models of cardiac tissue electrophysiology: Progress, challenges and open questions[END_REF] to ensure uniqueness on u e . This formulation is discretised in time using an Euler semi-implicit scheme: implicit for the diffusion and explicit for the reaction. Spatial discretisation has been led using Control Volumes Finite Elements (CVFE, see eg [START_REF] Cai | The finite volume element method for diffusion equations on general triangulations[END_REF]). One large linear system (symmetric and non-negative) has to be inverted at each time step. These inversions are the main numerical difficulty for the bidomain model: precise details on the implementation, and on the preconditioning of the bidomain model are available in [START_REF] Pierre | Preconditioning the coupled heart and torso bidomain model with an almost linear complexity[END_REF].

The adapted monodomain model implementation is much simpler. We consider the classical weak formulation of the parabolic equation [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF]. It is discretised as for the bidomain model: using an Euler semi-implicit scheme in time and a CVFE discretisation in space.

Remark 1. We eventually discuss the choice of the spatial discretisation. We numerically experimented that the result quality is highly sensitive to this choice. We compared 3 methods of order 2 (with respect to the mesh size) for the elliptic problem discretisation: P1 finite element, Discrete Duality Finite Volumes (DDFV, see eg [START_REF] Coudiere | 2D/3D DDFV scheme for anisotropic-heterogeneous elliptic equations, application to electrograms simulation from medical data[END_REF][START_REF] Andreianov | Convergence of DDFV schemes for the bidomain cardiac model[END_REF] as applied to the bidomain model) and Control Volumes Finite Elements (CVFE, see eg [START_REF] Cai | The finite volume element method for diffusion equations on general triangulations[END_REF]). The CVFE method gave the best results and will be used here. For the CVFE method the nodes are located at the mesh vertices, allowing to associate to the numerical solution a P1 function on the mesh elements, which property will be useful in the next subsection 3.2. Estimation of the errors. One searches to measure the discrepancy between the bidomain and the adapted monodomain model at the continuous level. For this one needs to get rid of the errors caused by the discretisation. We proceed as follows:

-Firstly, we numerically study the convergence of the numerical schemes for both the bidomain and the adapted monodomain models, using series of successively refined meshes. An a posteriori error estimate then allows to evaluate the error induced by discretisation. -Secondly, the error between the (discrete) bidomain and adapted monodomain models is computed on each mesh. Once this error becomes wider than the discretisation errors measured in step one, one can deduce the discrepancy between the (continuous) bidomain and adapted monodomain models. Because of its particular physiological importance, we adopt an error criterion based on the activation time. Activation time φ : Ω → R is defined at each point x as the time φ(x) such that v(φ(x), x) = v s for a threshold value v s set here to -20 mV . The value φ(x) tells us at what time the excitation wave reaches the point x. Activation time for test cases 1 and 2 are depicted on Figs. 2 and5 respectively.

Let us denote by T n the series of meshes. We consider the numerical activation time φ n computed relatively to T n . We compare φ n-1 and φ n by introducing the projection φn-1 of φ n-1 on the finer level T n . For this φ n-1 is considered as a P1 function (i.e. continuous piecewise affine) on the elements of T n-1 . The values of this function are computed on the vertices of T n , which data allows to define the considered projection φn-1 as a P1 function on the elements of T n . The following relative error in L 2 norm is computed, as well as the (logarithmic) convergence rate r n 2 with respect to the mesh number of vertices denoted D n :

(11) e n 2 := φn-1 -φ n L 2 (Ω) , r n 2 := log(e n-1 2 /e n 2 ) log(D n-1 /D n ) . An error estimator ε n 2 , φ n -φ ∞ L 2 φ ∞ L 2 ≤ ε n 2 ,
on the relative error between the discrete activation time mapping φ n on T n and the activation time mappings φ ∞ for the continuous problem is defined as follows. The asymptotic behaviour of e n 2 is numerically analysed using r n 2 . Assuming that e n 2 has a geometrical behaviour with geometrical growth rate q, ie e n 2 qe n-1 2 , then we can derive the upper bound:

φ n -φ ∞ L 2 (Ω) ≤ φ n -φ n+1 L 2 (Ω) + φ n+1 -φ n+2 L 2 (Ω) + . . . ≤ q 1 -q φ n -φ n-1 L 2 (Ω) .
We thus will define

(12) ε n 2 := q 1 -q e n 2 φ M L 2 (Ω) ,
with M the index of the finest mesh T M . We precise that saying that e n 2 has an asymptotic geometrical behaviour is equivalent with saying that r n 2 goes to a constant r, in which case q can be deduced from r.

We define the (relative) discrepancy between the bidomain and adapted monodomain models as:

(13) δ n 2 := φ b n -φ m n L 2 (Ω) φ b n L 2 (Ω)
, where b and m stands for bidomain and adapted monodomain respectively.

3.3. Test cases.

3.3.1.

Unit square. The following academic test case is considered. The domain Ω is set to [0, 1] 2 , recalling that the cm is the reference length scale. A series of 6 meshes (T n ) n=1,6 of this domain, from 177 to 160 385 vertices, has been considered. These 6 meshes are triangle regular meshes that have been generated using the Matlab mesher: the (n + 1) th mesh is generated by cutting each triangle into 4 triangles. Thus, recalling that D n denotes the number of vertices of T n , we have D n+1 = 4D n . On this geometry the spreading of excitation potential waves is simulated as depicted on Fig. 2. Excitation is initiated by applying a centered stimulation during a short period of time, precisely: I st (x, t) = 52 µA for 5 < t < 5.1 and |x -x 0 | < 0.1 (x 0 denoting the centre of Ω) and I st (x, t) = 0 otherwise. Stimulation is applied at time t = 5 rather than t = 0 to ensure the ionic model to have reached its equilibrium rest state, which state might not be reached initially because of the large number of parameters in the Luo and Rudy II model and to possible rounding errors on these parameters.

The model parameters are set following the data in Tab. 1: three different values for χ will be considered: χ = 1000, χ = 1500 or χ = 2000 cm -1 . Initial condition ( 5) is taken as the rest state for the ionic model: I ion (v 0 , w 0 ) = 0 = g(v 0 , w 0 ). The domain Ω is assumed to be composed of a bundle of parallel horizontal muscular fibres, the following homogeneous tensors are defined: [START_REF] Ethier | Semi-implicit time discretization schemes for the bidomain model[END_REF] σ i (x) = σ i := Diag(g l i , g t i ) , σ e (x) = σ e := Diag(g l e , g t e ) . 3.3.2. 2D cut of a segmented heart. We now consider a realistic 2D setting depicted on Fig. 1. The geometry has been obtained from the segmentation of a medical image of a human heart [START_REF] Rousseau | Geometrical modeling of the heart[END_REF]: it represents a horizontal slice of the two ventricles with a resolution of 0.3 mm. The geometry has a surface of 27 cm 2 . Four meshes (T n ) n=1,4 of this geometry have been built using the algorithm Distmesh 1 in [START_REF] Persson | A simple mesh generator in matlab[END_REF]. These meshes precisely counting D 1 =117 285, D 2 =279 447, D 3 =551 484 and D 4 =1 110 270 vertices. Anisotropy is given from the fibrous organisation of the tissue following [START_REF] Ambrosio | On the asymptotic behaviour of anisotropic energies arising in the cardiac bidomain model[END_REF]. A set of fibres rotating around the two ventricles and remaining tangent to the boundary has been built. Stimulation is initiated at 4 sites (two on each ventricle) on the endocardium. A stimulation current of 52 µA is applied during 0.1 ms at these locations at time 20 ms and 25 ms on the left ventricle and right ventricle respectively. The stimulation site locations as well as the delay of 5 ms between the stimulation of the left and of the right ventricles is aimed to mimic experimental measurements on a perfused heart in [START_REF] Durrer | Total excitation of the isolated human heart[END_REF]. Stimulation is initiated at time t=20 ms rather than t = 0 for the same reason as in the Unit Square test case.

The model parameters are set following the data in Tab. 1, χ is here set to χ = 1500 cm -1 . the adapted monodomain model and for χ = 1500 cm -1 on the finest mesh. On Fig. 2 are also depicted the resulting activation time mappings for the bidomain and adapted monodomain models. As a first observation, activation time predicted by both models are in qualitative agreement, isolines displaying the same expected elliptical shape. A first and simple quantitative comparison criterion on the propagation velocities also shows the concordance between the two model predictions: the total time for the depolarisation wave to recover the whole domain is of 21.6 and 21.9 ms for the bidomain and monodomain models. Although the bidomain model appears to be slightly faster, the relative error between these figures is of 1.4 %. This estimation can be reinforced by comparing the relative error between the L 2 norm of the activation time mappings:

Results

(15) e 0 = φ mono L 2 -φ bid L 2 φ bid L 2 .
Computed on the finest mesh, we get e 0 = 0.85 10 -2 , 0.97 10 -2 and 1.07 10 -2 for χ = 1000, 1500 and 2000 respectively. Numerical results on the convergence e n 2 and convergence rate r n 2 of the activation time mappings in L 2 norm are displayed in Tab. 2 and Tab. 3 for the bidomain and adapted monodomain models respectively. For both models the convergence of the activation time mapping is clearly observed; it is moreover reasonable to assume an order -1 convergence with respect to the mesh number of vertices D n . This order -1 convergence justifies the definition of the error estimator ε n 2 in equation ( 12) taking

χ = 1000 cm -1 χ = 2000 cm -1
Figure 3. Comparison of the discretisation error ε n 2 (defined in [START_REF] Coudiere | 2D/3D DDFV scheme for anisotropic-heterogeneous elliptic equations, application to electrograms simulation from medical data[END_REF]) between the continuous and discrete solution and of the relative discrepancy δ n 2 (defined in ( 13)) between the bidomain and monodomain models activation times in L 2 norm. The comparison is depicted simultaneously for the bidomain and monodomain models and for two values of χ. q = 1/4, since D n is here multiplied by 4 each time n is incremented of 1.

Results on the error estimator ε n 2 are given in Tab. 4 for both models. They also are displayed on Fig. 3 for the two values of χ = 1000 and χ = 2000. Three immediate conclusions can be drawn from these figures. Firstly the discretisation errors are quite similar for the bidomain and the monodomain models. Secondly the discretisation error increases with χ. Thirdly, obtaining a precise approximation using physiological values of the parameter χ necessitates to consider very fine meshes. These three points will be discussed deeper on in the conclusion section 5.

The relative discrepancy δ n 2 between the (discrete) bidomain and adapted monodomain models is quantified in Tab. 5. Numerical results in Tabs. 4 and 5 are summarised on Fig. 3 for χ = 1000 and χ = 2000. These results clearly indicate that the discretisation error becomes smaller than the discrepancy between the two models on the two finest meshes. More precisely the discrepancy between the two models is larger than the discretisation error with roughly a factor 10 on the finest mesh. Thus the discrepancy between the bidomain and the adapted monodomain models on this test case appears to be smaller than 1%. This estimation can be refined by examining more carefully the data in Tab. 5. The sequence δ n 2 decreases towards its limit for the three considered values of χ, which limit we can try to extrapolate from the data. On Fig. 4 is depicted δ n 2 as a function of 1/D n (D n still denoting the number of vertices of the mesh T n ) in decimal logarithmic scale. For all the three considered values of χ the plots display an almost linear behaviour. The slopes have been evaluated using a least square approximation taking into account the 5 last data points only (i.e. skipping δ 1 2 ). The results respectively are of 2.98, 3.10 and 3.20 for χ = 1000, 1500 and 2000. This means a convergence of order -0.33, -0.32 and -0.31 respectively of δ n with respect to D n . From this we can conclude to a convergence of order -1/3 of δ n 2 with respect to the mesh number of vertices. The limit of δ n 2 is then extrapolated by considering the intersection of this least square best linear approximation of the 5 last data points with the y axis. The extrapolated limits for δ n 2 are of 0.26 %, 0.29 % and 0.30 % for χ = 1000, 1500 and 2000 respectively.

We eventually point out that although these extrapolated limit values are 2 to 3 times larger than the estimated discretisation errors in Tab 4 on the finest mesh, our numerical results remain too rough to provide a precise approximation such as 0.3 %. We can however conclude that on this test case, for the three considered values of χ, the relative discrepancy in L 2 norm between the activation time mapping predicted by the bidomain and the adapted monodomain models is 2 to 3 times smaller than 1 %. 4.2. 2D cut of a segmented heart. The spreading of transmembrane potential wave v across the 2D segmented heart geometry is depicted on Fig. 5 using the bidomain model (and for χ = 1500 cm -1 , which value is fixed in this section). The activation time mappings simulated using the bidomain and adapted monodomain models are displayed on the same Fig. 5. As for the first test case, these activation times are in good qualitative agreement. The total amount of time for the depolarisation wave to recover the whole domain is of 77.48 and 79.34 ms for the bidomain and monodomain models. The bidomain model again appears as being slightly faster, the relative error between these two times is of 2.4 %. Computing the relative error between the L 2 norms of the activation time mappings e 0 defined in [START_REF] Gerardo-Giorda | A model-based block-triangular preconditioner for the bidomain system in electrocardiology[END_REF] here gives e 0 = 1.18 10 -2 .

The numerical results on the convergence are given in Tab. 6. We here observe a convergence block between the third and the fourth mesh. Although the L 2 distance φn -φ n-1 L 2 for n = 2, 3 normally decreases (with a convergence rate r 3 2 -0.6 that could be in good agreement with a -1 order of convergence with respect to D n as observed for the previous test case), this distance stops decreasing between meshes 3 and 4. This convergence block moreover very similarly occurs for the bidomain and Bidomain Monodomain adapted monodomain models. Several assumptions can be advanced to explain this phenomenon, that all rely on the quite complex setting of this test case. The numerical geometry (the union of all triangles of the mesh T n ) is not fixed. Practically, the inputs for the mesh construction using the software Distmesh are the mesh size and a level set function. This level set function itself is numerically defined on the underlying segmentation grid. Because of the relative complexity of the segmented domain, T n and T n+1 may not exactly recover the same region, so inducing boundary errors for the projection φn of φ n on T n+1 . The anisotropy tensors are not analytically defined and so do not converge to some limit tensor: the anisotropy has been defined on the image segmentation grid and simply projected on each meshes. This also is a plausible sourde of error inducing the convergence block. Despite this convergence block, the discretisation error estimator ε n 2 however has been computed as in the previous test case (that is assuming a -1 order convergence with respect to the mesh number of vertices), the results also are given in Tab 6. From these figures, it seems acceptable to bound the relative discretisation error on the finest meshes to a few percent: 3 to 4 %. This range of discretisation error has to be put into perspective with the first test case. The domain surface here is of 27 cm 2 and the finest mesh approximately counts 40 000 vertices per cm 2 , and thus has the same resolution as mesh 5 of the unit square. The discretisation error was for this mesh however 6 to 8 times smaller than presently. The complex setting of this test case then induces a supplementary and significant amount of discretisation error. The relative discrepancy δ n 2 between the bidomain and adapted monodomain models itself do present a much finer behaviour that is depicted on Fig. 4, data are precisely reported in Tab. 7. In decimal logarithmic scale, δ n 2 displays a linear dependence with 1/D n , with slope 4.9 (using a least square best approximation of the data with a linear mapping). Thus δ n 2 seems to converge towards its limit with an order -1/5 convergence with respect to the mesh number of vertices. Here again we observe a slow down in the convergence as compared to the first test case where the convergence (already quite slow) had been evaluated to be of order -1/3 with D n . Extrapolating the limit from this data provides us a discrepancy between the continuous bidomain and adapted monodomain models of 0.2 %. Because of the convergence block previously discussed, and to the amount of discretisation error evaluated to 3 to 4 %, it is not reasonable to raise such a conclusion. From the quality of the result on δ n 2 , we however conclude to a discrepancy of order 1 % or even below 1 % between the activation time mappings of the continuous bidomain and adapted monodomain models.

Discussion

This paper was aimed to present a precise and numerically documented comparison of the bidomain and adapted monodomain models for the simulation of cardiac electro-physiology. The models definition, their numerical implementation and the general methodology for their comparison has been clearly stated. The comparison methodology relies on two points: the definition of a physiologically relevant comparison criterion (based on activation time comparison) and a numerical convergence analysis of the two models to get rid of discretisation errors. Comparison has been led on two test cases both involving complex ionic dynamics and realistic macroscopic tissue parameter settings: the first one for an academic unit square geometry and the second one for a complex 2D cut of segmented human ventricles.

Two main conclusions are drawn by analysing these two test case results,

(1) the amount of error caused by discretisation for commonly used grid resolutions is quite high, (2) the discrepancy between the two model predictions is very small, of order 1 % or even below 1 % of relative error in L 2 norm.

5.1. Discretisation error. The (relative) discretisation error on the activation time has been evaluated on the two test cases. For both test cases, this amount of error has been quantified to be quite large even for fine grids (i.e. fine when compared to the mesh size commonly used in biomedical engineering). This error moreover increases when going towards complex settings. On the contrary, skipping from the bidomain to the monodomain model has no real influence on this error. To the author appreciation, once the ionic model set, the main factor of discretisation error is the value of the parameter χ. The influence of χ on the discretisation error has been evaluated for the unit square test case. As expected, the highest χ, the largest the discretisation error. Indeed, for both the bidomain and the monodomain models, χ could be rescaled to the value of 1 by changing the space variable x into x √ χ. Therefore the mesh number of vertices has to be proportional to χ d/2 (with d the space dimension). More phenomenologically, the greatest χ the sharpest the wave fronts: thus more points are needed to obtain a good approximation inside the wave front location. The number of points needed transversally to the wave front then is proportional to √ χ.

The physiological range of values for χ (from 2400 to 8400 cm -1 in [START_REF] Clayton | Models of cardiac tissue electrophysiology: Progress, challenges and open questions[END_REF]) then imply very important computational efforts. Numerical simulations based on the bidomain or the monodomain models commonly depreciate this value (from 1000 to 2000 in this paper, χ =1000 in [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models[END_REF][START_REF] Mc | Simulation of qrst integral maps with a membrane-based computer heart model employing parallel processing[END_REF], χ = 500 in [START_REF] Smith | Altered t wave dynamics in contracting cardiac model[END_REF][START_REF] Buist | The effect of torso impedance on epicardial and body surface potentials: a modeling study[END_REF], e.g.). As a simple illustration: to predict activation times with a tolerance relative error of 5 % using the (physiologically still under estimated) value χ = 2000 cm -1 necessitates to require a mesh resolution of 0.1 mm on the unit square test case or even less for more complex simulation. This would mean more than 1 million of points per cm 3 in dimension 3. Asking for a 1 % relative error would increase this figure to 8 million points per cm 3 (0.05 mm resolution grid). Although these estimations are quite rough, context depending and could be decreased using more accurate numerical methods (increasing the order e.g.), we conclude that running accurate simulations of cardiac electrical activity remains challenging in this context for physiological values of the parameter χ.

Model comparison.

For the two test cases the discrepancy has been evaluated to be of order 1 % in terms of L 2 activation time relative error.

The first test case involved a simple anisotropy type for the tissue and a simple central excitation that resulted into basic patterns of excitation: elliptical isolines for the activation time. The discrepancy in this case has been evaluated to be even smaller than 1 % and was bounded above by 1/2 or 1/3 % (uncertainties remaining due to numerical noise). The second test case involved a much more sophisticated setting in terms of geometry, anisotropy and stimulation settings inducing enhanced patterns of excitation. Although the discrepancy analysis was made harder in this case, we concluded to an order 1 % or even below of relative error between the two models.

The following two points are interesting to notice. Firstly the discrepancy was measured for three values of χ in the first test case: 1000, 1500 and 2000 cm -1 . In first approximation the influence of χ is not visible. Secondly the discretisation imply an over-estimation of the model discrepancy: in Tabs. 5 and 7, δ n 2 decreases with the mesh size. Although we obtained a comparable order of discrepancy between the models for two quite different test cases, our reported order of 1 % of error cannot be generalised and its dependence with the model parameters still need to be evaluated. As already mentioned χ has no important influence on the model discrepancy. However the conductivities might have a great impact since in case of equal anisotropy the discrepancy is of zero. Considering the following ratio ρ: ρ = g l i /g t i g l e /g t e , equal to 1 in case of equal anisotropy ratio and satisfying ρ > 1 in physiological settings. The model discrepancy is expected to grow up with ρ.

The only reported bidomain and monodomain discrepancy in terms of activation time comparison available in the literature (and also computed using L 2 norms) is available in [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models[END_REF]. In this paper Colli-Franzone et al. measured this discrepancy on a 3D slab of tissue either with axisymmetric or orthotropic anisotropy. They measured a significantly higher discrepancy than we did: 6.3 % and 11.3 % for the two different anisotropic settings respectively. Although they considered a different ionic model (Luo and Rudy type I model) and a different value of χ (1000 cm -1 ) we think that the observed difference rather lies in the two following points. Firstly the conductivities in [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models[END_REF] imply a value of ρ = 7.1 whereas ρ = 4.5 here: we then are further from the equal anisotropy ratio in this paper. Secondly the mesh size in [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models[END_REF] is of 0.1 mm. The discretisation error on such a fine grid might remain high in a complex 3D anisotropic setting leading to an over estimation of the model discrepancy.

Conclusion.

The discrepancy between the (continuous) bidomain and adapted monodomain models has been evaluated to be of order 1 % in terms of relative error on the activation time. This is quite small: such a range of error is negligible for a field (biological sciences) where the tolerance on the data errors is much larger. Meanwhile the amount of error induced by the discretisation shall be even higher than the model discrepancy for commonly used mesh resolution in biomedical engineering. Therefore, to the authors feeling, the monodomain model should be preferred to the bidomain model to simulate patterns of excitation in the cardiac tissue taking advantage of its lightest implementation and computational cost.

Figure 1 .

 1 Figure 1. Test case 2 settings: stimulation sites location (above left) and fibrous structure of the tissue (above right).

4. 1 .Figure 2 .

 12 Figure 2. Above: spreading of transmembrane potential wave after stimulation (for the monodomain model, χ = 1500). Potential maps are represented 3, 8 and 12 ms after the stimulation. Below: comparison of the activation time mappings predicted by the bidomain (left) and adapted monodomain (right) models, on the finest mesh and also for χ = 1500 cm -1 .

Figure 4 .

 4 Figure 4. Extrapolation of the limit for δ n 2 as n → +∞. The discrepancy between the two models is depicted as a function of 1/D n in decimal logarithmic scale. Left: unit square test case, for the three considered values of χ. Right: 2D cut of a segmented heart test case.

Figure 5 .

 5 Figure 5. Above: spreading of transmembrane potential wave across the segmented heart geometry, simulated using the bidomain model on the finest mesh, χ = 1500 cm -1 . Potential maps are represented 16, 32 and 48 ms after the stimulation. Below: comparison of the activation time mappings predicted by the bidomain (left) and adapted monodomain (right) models, also on the finest mesh. Isolines (in black) are separated by 10 ms. Stimulation is initiated at time t = 20 ms.

DistMesh: "a simple mesh generator in Matlab", http://wwwmath.mit.edu/ persson/mesh/

 6. 2D cut of a segmented heart. Convergence of the activation time φ n for the bidomain and the monodomain models. The error indicator e n 2 and convergence rate r n 2 being defined in [START_REF] Colli-Franzone | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF]. The discretisation error estimator ε n 2 is defined in [START_REF] Coudiere | 2D/3D DDFV scheme for anisotropic-heterogeneous elliptic equations, application to electrograms simulation from medical data[END_REF]. 7. 2D cut of a segmented heart. Discrepancy between the (discrete) bidomain and adapted monodomain models on the heart slice geometry (χ = 1500 cm -1 ).

(Yves Bourgault) Department of mathematics and statistics, University of Ottawa.

E-mail address: ybourg@uottawa.ca (Charles Pierre) Laboratoire de Mathématiques et Applications Université de Pau et du Pays de l'Adour av. de l'Université BP 1155 64013 PAU Cedex -FRANCE E-mail address: charles.pierre@univ-pau.fr


Junle Wang

Marcus Barkowsky

Vincent Ricordel

Patrick Le Callet

Quantifying how the combination of blur and disparity affects the perceived depth

Keywords: 3DTV, Depth perception, Stereoscopic display, Gaussian blur, Subjective experiment

The influence of a monocular depth cue, blur, on the apparent depth of stereoscopic scenes will be studied in this paper. When 3D images are shown on a planar stereoscopic display, binocular disparity becomes a pre-eminent depth cue. But it induces simultaneously the conflict between accommodation and vergence, which is often considered as a main reason for visual discomfort. If we limit this visual discomfort by decreasing the disparity, the apparent depth also decreases. We propose to decrease the (binocular) disparity of 3D presentations, and to reinforce (monocular) cues to compensate the loss of perceived depth and keep an unaltered apparent depth. We conducted a subjective experiment using a twoalternative forced choice task. Observers were required to identify the larger perceived depth in a pair of 3D images with/without blur. By fitting the result to a psychometric function, we obtained points of subjective equality in terms of disparity. We found that when blur is added to the background of the image, the viewer can perceive larger depth comparing to the images without any blur in the background. The increase of perceived depth can be considered as a function of the relative distance between the foreground and background, while it is insensitive to the distance between the viewer and the depth plane at which the blur is added.

INTRODUCTION

Recently, stereoscopic image and video production is gaining an increasing amount of attention. The displays nowadays used to show these 3D productions are usually planar so that binocular disparity on the display plane becomes a preeminent depth cue enabling viewers to perceive depth. As a binocular cue, disparity is stable, but on the current planar displays, it induces a conflict between accommodation and vergence of the eyes. This conflict is usually considered as one main reason for visual discomfort, especially when the disparity is large. If we decrease this binocular cue, disparity, to limit the visual discomfort, the apparent depth also decreases.

Several other depth cues besides binocular disparity affect also the apparent depth. We propose to decrease the (binocular) disparity of 3D presentations, and to reinforce (monocular) cues to compensate the loss of perceived depth and keep an unaltered apparent depth. The limitation of depth-of-field of human eyes causes blur in the retinal image which is known as an important monocular depth cue. Some previous investigations have shown clear contributions of blur to depth perception [START_REF] Pentland | A new sense for depth of field[END_REF][START_REF] Watt | Focus cues affect perceived depth[END_REF][START_REF] Hoffman | Focus information is used to interpret binocular images[END_REF][START_REF] Held | Using blur to affect perceived distance and size[END_REF] , while others showed that blur has either no effect [START_REF] Mather | Depth cue integration: stereopsis and image blur[END_REF] or only some qualitative effects on perceived depth ordering [START_REF] Marshall | Occlusion edge blur: a cue to relative visual depth[END_REF][START_REF] Mather | Image blur as a pictorial depth cue[END_REF][START_REF] Palmer | Edge-Region Grouping in Figure-Ground Organization and Depth Perception[END_REF] . However, most of the previous literature measures the influence of the blur cue by adding blurriness only at the screen plane, while the depth cues' magnitude varies as a function of distance. Another limitation in most of the previous studies was the experiment apparatus. CRT monitors were used and this type of displays has some limitations: (1) the surface containing stimuli was slightly curved, (2) the stimuli's virtual distance was affected by refraction due to the front glass plate, and (3) the screen was usually not large enough to cover a favorable field-of-view. Depending on the monitor used, the speed of the phosphor changes might also have caused crosstalk.

In our study, we conducted the subjective experiment using a state-of-art stereoscopic display system, the Samsung SyncMaster 2233RZ which is a 22.5-inch 1680*1050@120Hz wide-screen LCD monitor working with active shutter glasses from NVidia. The stimuli used in the experiment contained a background plane and a single object in the foreground, both of which were chosen closer to natural content compared to the stimuli used in previous publications. The image of a butterfly was used as the foreground object, since it is spatially complex enough, containing regions with both low and high frequency.

In our experimental observations, observers viewed stimuli in a two alternative forced choice (2AFC) task, being required to select the stimulus with largest depth interval between foreground and background. Two sources of perceived depth are used: disparity and blur. The perceived depth from disparity stems from the difference of disparity between the foreground object and the background. The perceived depth from blur stems from the amount of blur introduced to the background by convolution with a Gaussian kernel. Both the absolute position and relative distance between the foreground and background stay as a free parameter. This setup is able to evaluate how the combination of disparity and blur affects the perceived depth of objects located at different distance.

This paper is organized as follow: Section 2 introduces the geometry about how the binocular disparity and the defocus blur are created. Section 3 describes the experimental methods. Section 4 presents the results and the analyses of the experiments. Some general conclusions are presented in Section 5.

DISPARITY AND DEFOCUS BLUR

When people fixate an object in a three-dimensional scene, they can perceive defocus blur and binocular disparity simultaneously. Actually, the creations of defocus blur and binocular disparity in the retinal image have the same fundamental geometry. We first consider blur.

An ideal thin lens focuses parallel rays to a point on the opposite side of the lens. The distance between this point and the lens is the focal length, . Light rays emanating from a distance in front of the lens will be focused to another point at distance at the opposite side of the lens. The relationship between these distances derives from the thin-lens equation:

1 1 1 ,
When the lens focuses at an object at distance , the image of this object will be formed at the image plane at distance . This object is now in focus, while objects at other distance become out of focus and hence generate blurred images on the image plane (Figure 1). We express the amount of blur in the image plane by the blur circle diameter , which can be computed by the equation:

1 ,
where A is the diameter of the lens (pupil) and is the posterior nodal distance (approximately the distance from the pupil to the retina). Note that the human eyes have imperfect optics, and this equation does not incorporate any of the eye's aberrations. Nevertheless, previous researches showed that this equation can provide an accurate approximation of blur when the eye is defocused [START_REF] Cheng | A population study on changes in wave aberrations with accomodation[END_REF] . Based on this calculation, blur of the retinal image can be modeled as a twodimensional Gaussian function. We consider the blur circle radius (c 1 ) as the radius of decay to exp(-0.5), i.e. the standard deviation sigma of the Gaussian function. Thus, given A and S 0 , the magnitude of blur varies according to the focus distance and the distance between the object focused and the object defocused . The magnitude of blur increases as the distance from the fixated object increases. When an object is located outside of the depth of field, this object will be projected in the eye with a perceivable magnitude of blur. This relation enables blur as a cue of depth.

Based on the same geometry (Figure 2), disparity provides also depth information to the human brain. When the two eyes separated by a distance p converge on an object at distance d o , another object at distance d 1 creates images with an angular disparity δ which can be obtained by the following equation using small-angle approximation: We can consider that disparity is caused by differing two vantage points of two cameras (eyes), while blur is caused by differing the vantage point at two positions of one camera (eye). By converting the equation of blur circle radius into angular units, and rearranging both equations, we can find the relationship between the magnitudes of disparity and blur:

|δ| .

Therefore, the magnitudes of blur and disparity created by objects in a three-dimensional scene are proportional to each other. This relationship generally holds because accommodation and vergence are coupled in the real world, namely the eyes accommodate and converge to the same distance. It suggests that the human visual system might use information from both the two factors in correlating the two eyes' images. However, when people watch 3D images or videos displayed on a planar stereoscopic display, there exists a conflict between accommodation and vergence. No matter how the vergence changes, the image is always perceived sharpest when the eyes accommodate on the screen plane. This conflict suggests that a manual creation of blur on the focal plane (the screen) might inhibit the conflict and affect the perceived depth.

EXPERIMENT

Method

Subjects

Thirty-five subjects participated in the experiment. Twelve subjects are male, twenty-three are female. The subjects did not know about the purpose of the experiment. Subjects ranged in age from 17 to 40 years. All had either normal or corrected-to-normal visual acuity. They were tested for visual acuity using a Snellen Chart, for depth acuity using a Randot Stereo Test and for color vision using Ishihara plates.

Apparatus and stimuli

Stimuli were displayed on a Samsung 22.5-inch LCD screen (1680 * 1050 at 120 Hz). Each screen pixel subtended 65.32 arcsec. The display yielded a maximum luminance of about 50 cd/m2 when watched through the activated shutter glasses. Stimuli were viewed binocularly through the Nvidia active shutter glasses (Nvidia 3D Vision kit) at a distance of approximately 90 cm. The peripheral environment luminance was adjusted to about 44 cd/m2. When seen through the eye-glasses, this value corresponded to about 7.5 cd/m2 and thus to 15% of the screen's maximum brightness as specified by ITU-R BT.500.

Each stimulus consisted of one single object in the foreground and a background. A 400 * 400 pixels butterfly image was used as the foreground object. This stimulus was easy to accommodate because it was spatially complex and therefore contained a wide range of spatial frequencies from low to high. The butterfly was shown with a horizontal offset between the left and the right view in order to create a disparity cue. The magnitude and direction of this offset varied to supply a variety of near or far disparity. At the 90 cm viewing distance, the background plane subtended 29.8 * 18.9 arcdeg, and the foreground object subtended 7.2 arcdeg. The background plane was a photo of a flowerbed. The background contained a great amount of textures, contained no distinct region of interest distracting the observers, and looked natural. When required, the background was spatially blurred by applying a Gaussian blur kernel with a 5-pixel radius. This amount of blur equals to the blur created by supposing that the focused object (foreground) and the defocused object (background) are at a distance of 13.7 cm in front of the screen and 19.8 cm behind the screen respectively, both of which are the limitations of the comfortable viewing zone derived from the 90 cm viewing distance. Note that all the blurred background were with the same amount of blur.

Design and procedure

In each trial, a pair of stimuli were shown to the subjects. One stimulus contained a blurred background (BB) and a sharp foreground object, while the other stimulus contained a sharp background (SB) and also a sharp foreground object (Figure 3). The backgrounds of both stimuli were positioned at the same depth. Two parameters of the stimuli were varied: the absolute position (D a ) of the background plane and the relative distance (D r ) of the foreground object to the background. The selection of absolute distance ranged from -19.7 cm to 6.6 cm in steps of 6.6 cm (negative values denote the positions behind the screen plane, and positive values denote the positions in front of the screen plane). The relative distance ranged from 0 cm to approximately 33 cm. All the positions of both background and foreground objects are selected considering the limitations of the comfortable viewing zone [START_REF] Chen | New requirements of subjective video quality assessment methodologies for 3DTV[END_REF] . There are thus twenty combinations of absolute position and relative distance in total for the BB-stimuli. Once the absolute position and the relative distance were selected for each BB-stimulus, then we paired the BB-stimulus with a set of 7 or 8 SB stimuli which are with relative distances ranging from 4 cm less than the BB stimulus to 8 cm larger than the reference stimulus. These steps were chosen considering both the depth rendering ability of the screen and the depth perception ability of the observers. One trial contains one BB-stimulus and one SB-stimulus. We had 155 trials in the experiment. This setup is shown in Figure 4, and the parameters are shown in Table 1. 1.

The subjective experiment was conducted using a two-alternative forced choice task. For each trial, one of the 155 conditions was chosen randomly and a pair of stimuli were displayed. Observers were asked to look at the butterfly and then determine in the trial whether the BB stimulus contained a larger depth interval between the butterfly and the background than the SB stimulus. As a two monitor screen setup was technically not possible, observers were able to control the displaying of stimuli by means of a key press to switch from one stimulus to the other. When the observer switched between the stimuli, a 700 ms grey interval was shown in order to avoid memorization by the observers of the "exact" positions of the foregrounds. For each trial, the total observation time and the number of switches were not limited.

D The selection of parameters D a and D r for both BB-stimuli and SB-stimuli. Note that there are only 7 possible selections with the D r_BB = 0 cm cases, because putting the foreground behind the background will cause trouble of fusion.

RESULT AND ANALYSIS

For each condition, some observers considered the BB-stimulus as having a larger depth interval (between the foreground and the background), while the other observers chose the SB-stimulus. We measure the proportion of 'BBstimulus contains a larger depth interval' responses, and plot the data as a function of the disparity difference between the D r in the BB-stimulus (D _BB ) and the D r in the SB-stimulus (D _SB ). The cumulative Weibull function was used as the psychometric function. The disparity difference corresponding to the 50% point can be considered as the Point of Subjective Equality (PSE). When measuring the disparity difference at that point, the increase of perceived depth is obtained. In total, by filtering out the data of 7 observers who made decisions in the test quite differently from other observers, 28 observations of each conditions were included in the computation. An example pattern of response and the fitted psychometric function is shown in Figure 5. In this trial, we consider D r_BB = 6.6 cm and D a = -19.7 cm, -13.2 cm, -6.6 cm, 0 cm, 6.6 cm. An equal apparent depth is reached at -220 arcsec, which corresponds to the second point in Figure 7 (right).

According to the setup of the experiment, blur was added to the backgrounds located at different absolute depth positions in different trials, while the depth interval between the background and the foreground object stayed also as a free parameter. We thus plot the curve of the PSE as a function of the background disparity representing the depth of the blurred background in Figure 6. Each of the five points on the curve is obtained by considering all the possible depth intervals (D _BB ) which are with the same background depth (D a ). The steps are slightly different because the same distances equal to various disparity in angular units depending on the viewing distance.

As we can see in the Figure 6, the PSE curve shows a clear offset from the unit slope, indicating an increase of perceived depth caused by blur in binocular vision. Note that the increase is almost constant (approximately 180 arcsec) regardless of the change of absolute depth of the background. This means that the increase of perceived depth caused by blur is insensitive to the disparity of the blurred background, e.g. the absolute position. This phenomenon makes sense because all the blur actually exists only on the screen plane where human eyes accommodate. We plot also the curve of the PSE as a function of the depth interval (in length unit) between the foreground and background in Figure 7 (left). Each of the six points on the curve is obtained by considering all the possible background depth (D a ) which are with the same depth interval (D _BB ). 0 cm means that the foreground and the background are at the same depth without any difference of disparity between them. According to the viewing geometry, the same depth intervals stand for slightly different disparity, therefore we plot the increase of perceived depth in disparity as a function of depth interval in length unit.

Figure 7 (right) shows also a clear increase of perceived depth created by the influence of blur. In the figure, we can find that the enhancement of perceived depth increases with increasing depth interval between the fixated sharp foreground object and the defocused blurred background. When the foreground object is close to the background (the relative distance is 0 cm or 6.6 cm), the influence of blur is relatively small. The 0 cm depth interval comes with a slightly larger increase of perceived depth than 6cm depth interval. The reason might be because when the depth interval is 0cm, it seems to the observers that they are looking at a planar image with a sharp foreground and blurred background. For this kind of image, intentional blur has been introduced by photographers for a long time to induce some illusions of the existence of depth between the foreground and the background. When disparity differences appear in the image, observers start to notice that they are facing a stereoscopic image. This might be the reason why the 0cm depth interval shows a slightly larger increase of perceived depth than the 6cm depth interval. Starting from the second point (6.6 cm depth interval), the enhancement is monotonic. As shown by the curves in Figure 7, the overall tendency of how the depth interval (D _BB ) and the absolute position of the background (D a ) affect the perceived depth seems easy to understand. However, we split the data from Table 1 and computed the PSEs of every combination of D _BB and D a , we then found that the variations of PSEs are with some uncertainties. In Figure 8, we plot a surface (consisting of blue points) which is obtained from the curve in Figure 7 (right), while the red points come from the PSEs which are computed individually by each combination of D _BB and D a .

The figure shows that the changes of perceived depth do not vary strictly according to the tendency we described previously. There exists a lot of variability at certain points. Namely, despite the overall tendency of how the blur affects perceived depth at different relative distance and absolute position is known, the predictions of perceived depth of objects at certain positions or the perceived depth of certain observers remain uncertain. Several types of analysis were performed in order to learn about the source of these aberrations. However, there was no indication found concerning a bi-or multimodal distribution of the observers. So the variance that is found in our experiment concerning the individual data points has to be considered as measurement noise for the moment.

CONCLUSIONS

The influence of defocus blur on the perceived depth in a stereoscopic scene and its relationship with binocular disparity were studied in this paper. The experimental result indicates that image blur makes contributions to the impression of depth perceived in stereoscopic images, when measured against depth perceived in stereoscopic images without blur. The increase of depth can be considered as a function of the relative distance between the fixated foreground ground object and the blurred background, while this increase is insensitive to the distance between the viewer and the depth plane at which the blur is added.

The feasibility of enhancing the perceived depth by reinforcing a monocular cue, namely defocus blur, provides an interesting way to deal with the conflict between accommodation and vergence when 3D images are shown on a planar stereoscopic display. Generally, the foreground object popping out of the screen is the most important object in the scene, while the large disparity of the object may lead to visual discomfort when it is actually fixated by the observer. Our results show that it is feasible to decrease the disparity of this object without losing its pop-out effect by adding some blur on its background.

Figure 1 :

 1 Figure 1: Schematic diagram of the generation of defocus blur.

Figure 2 :

 2 Figure 2: Schematic diagram of the generation of binocular disparity.

Figure 3 :

 3 Figure 3: Example stimuli: the left views of an SB-stimulus and a BB-stimulus.

Figure 4 :

 4 Figure 4: Schematic diagram of the experiment setup.The blue planes represent the (foreground) depth planes in which the butterfly is located. The red plane represents the depth plane in which the background is located. We named the distance between background plane and the screen as Absolute Position (D a ), while the interval between the foreground plane and the background plane as the Relative Distance (D r ). Both these two distances were free parameters of the design of stimuli. The variation range of these two planes always stayed in the comfortable viewing zone. The first figure shows the BB-stimulus case, while the second shows the SB-stimulus case. Each BB-stimulus is paired with one of 7 or 8 SB-stimuli to create one trial. The selection of the parameters D a and D r for both BB-stimuli and SB-stimuli are presented in Table1.

Figure 5 :

 5 Figure 5: An example pattern of the proportion of observers' responses and the fitted psychometric function. In this trial,we consider D r_BB = 6.6 cm and D a = -19.7 cm, -13.2 cm, -6.6 cm, 0 cm, 6.6 cm. An equal apparent depth is reached at -220 arcsec, which corresponds to the second point in Figure7(right).

Figure 6 :

 6 Figure 6: The PSE curve shows subjective matches between SB-stimuli and BB-stimuli. Negative values denote far disparities. The dashed line at unit slope indicates physically correct disparity matches.

Figure 7 :

 7 Figure 7: Results from the experiment. The left figure shows the increase of perceived depth as a function of absolute position, while the right figure shows the increase of perceived depth as a function of relative distance. Note that each point on the curves might be obtained from different number of observations.

  Depth (arc sec) Relative Distance (cm) Increase of Perceived Depth (arc sec)

Figure 8 :

 8 Figure 8: The surface consisting of blue points is generated from the curve in Figure 7 (right). We ignore the variance in the curve in Figure 7 (left), because the variance is relatively small. The red points are obtained from the PSEs which are computed individually by each combination of relative distance (D _BB and absolute position (D a ).

Table 1 :

 1

	a (cm)	D r_BB (cm)				D r_SB (cm)
	-19.7	0		0	1.2	2.3	3.4	4.5	5.6	6.6
	-19.7	6.6	3.4	6.6	7.6	8.6	9.6 10.6 11.5 13.4
	-19.7	13.2	11.1 13.2 14.3 15.6 16.9 18.1 19.7 21.7
	-19.7	19.7	17.3 19.7 20.9 22 23.5 24.9 26.6 28.3
	-19.7	26.3	24.9 26.3 27	28 29.3 30.5 32 33.4
	-19.7	32.8	30.5 32.8 34 35.1 36.2 37.5 38.7 39.9
	-13.2	0		0	1	2	3	4	4.9	5.9
	-13.2	6.6	4.5	6.6	7.7	8.6	9.8 11.1 12.7 14.3
	-13.2	13.2	11.5 13.2 14.3 15.4 16.5 17.6 18.7 19.7
	-13.2	19.7	18.3 19.7 20.7 21.4 22.3 23.3 24.5 26
	-13.2	26.3	24.5 26.3 27.1 27.9 29 30.1 31.4 32.6
	-6.6	0		0	0.9	1.8	2.7	3.5	4.4	5.2
	-6.6	6.6	3.5	6.6	8	9.1 10.2 11.6 13.4 15.4
	-6.6	13.2	12 13.2 14.4 15.4 16.3 17.3 18.5 19.6
	-6.6	19.7	16.9 19.7 20.8 21.6 22.7 23.8 25.3 26.8
	0	0		0	0.8	1.6	2.3	3.4	4.5	5.6
	0	6.6	4.5	6.6	7.6	8.6	9.5 10.8 12.3 14.3
	0	13.2	11.4 13.2 13.7 14.3 15.1 15.9 17.2 18.7
	6.6	0		0	1	2	2.9	3.9	4.8	6.3
	6.6	6.6	4.2	6.6	7.1	8	8.8	9.9 11.2 12.9