The Beltrami Flow over Triangulated Manifolds Lucero Lopez-Perez1 , Rachid Deriche1 , and Nir Sochen2 1
Odyssee Project, INRIA, 06902 Sophia Antipolis, France {Lucero.Lopez Perez, Rachid.Deriche}@inria.fr http://www-sop.inria.fr/odyssee 2 Department of Applied Mathematics, University of Tel-Aviv, Tel-Aviv 69978, Israel
[email protected] http://www.math.tau.ac.il/∼sochen/index.html
Abstract. In several image processing applications one has to deal with noisy images defined on surfaces, like electric impulsions or diffusion tensors on the cortex. We propose a new regularization technique for data defined on triangulated surfaces: the Beltrami flow over intrinsic manifolds. This technique overcomes the over - smoothing of the L2 and the stair-casing effects of the L1 flow for strongly noised images. To do so, we locally estimate the differential operators and then perform temporal finite differences. We present the implementation for scalar images defined in 2 dimensional manifolds and experimental results.
1
Introduction
To regularize data on manifolds, differential geometry and calculus of variation are frequently used in two principal ways: The Polyakov action [8–10, 3, 4], which uses an intrinsic-parametric description of the manifold and an explicit form of the metric; and The Harmonic maps [1, 13, 14, 2, 17], which uses an implicit representation of the manifolds. The relation between these two approaches is explained in [18–20], in which is introduced a new approach to perform regularization on manifolds: The Beltrami Flow. It has been shown that the Beltrami flow interpolates between the L2 norm flow and the L1 norm flow, for flat gray-value images [9, 11], and for general non-flat surfaces [18–20]. In [19] is shown its implementation for the case of a manifold represented by a level set surface, but this has not be done for the case of intrinsically represented manifolds, often used in image processing. Even though the level set approach is easier to implement, for certain applications it is better to handle triangulated-based techniques rather than implicit ones. For example when the regularization is an intermediate step in the process to achieve a further goal, and the same triangulation is needed for the input and output of the step.
2
To compute the Beltrami Flow, we propose a finite differences scheme for the time discretization and a spatial average of the differential operator on the vertex of the triangulation. This article is organized as follows: Section 2 reviews the intrinsic formulation of the problem for the case of a 2D surface with scalar data, Section 3 shows how to estimate the differential operator. Section 4 presents the implementation details to compute the flow and some examples and results and we conclude in Section 5.
Fig. 1. Triangulated cortex
2
Intrinsic Formulation
Let S be a 2-dimensional manifold described by an arbitrary parameterization with coordinates u, v, that induce a metric (˜ gij ). It is given a scalar function I : S → IR, the data function. As it is explained in [20], we introduce the manifold Σ that describes I(u, v) as a surface embedded in the 3D fiber bundle M = S × IR Assuming an isometric embedding i.e. dsΣ = dsM , the metric gi,j induced on the section Σ is: ¡ ¢ dsΣ = g˜11 du2 + 2˜ g12 dudv + g˜22 dv 2 + β 2 Iu du2 + 2Iu Iv dudv + Iv2 dv 2 , where β 2 takes into account the differences in dimensions between the spatial directions and the intensity. This metric can be used to calculate the Laplace Beltrami Operator ¡√ ¢ 1 ∆g I = √ Div gG−1 ∇I . (1) g
3
where G is the matrix whose elements are gij , and g = det(G). The equation of motion that results from the Polyakov action [5] is: It = ∆g I It is shown in [20] that this flow becomes the L2 flow if we take the limit when β → 0 and it becomes the L1 flow if we take the limit when β → ∞
3
Estimating the Differential Operators
Let S be a surface parametrized by P (u) = {x(u), y(u), z(u) : u = (u1 , u2 ) = (u, v) ∈ D} Let G = (Gi,j ), Gi,j = hPi , Pj i be the Riemannian metric tensor for this parametrization. Then, for a function I : S → R ∇P I(x) = ∇u,v I(x) =
2 X i,j=1
Gi,j
∂I Pi ∂ui
where (Gi,j ) = G−1 We will use the (u, v)-space, with the metric induced by the local parameterization so that (1) can be expressed as: µq ¶ 1 2 Div 1 + β|∇u,v I| ∇u,v I , I t = ∆g I = p 1 + β|∇u,v I|2 Let us call ∇S I(x) = ∇u,v I(x) for x ∈ S, and 1 1 + β|∇u,v I(x)|2
φ(|∇S (I(x))|) = p
(2)
d We look for a local approximation ∆ g I on the nodes of S . p being a node of the triangulated surface, and φ a function φ : IR → IR. Following the ideas in [22] , we take the spatial mean of this quantity in the area A defined by the triangles inmediatly surrounding p (see fig 2) Z 1 ∇S · (φ(|∇S I(x)|)∇S I(x))dx ∇S · (φ(|∇S I(p)|)∇S I(p)) ≈ A A By the Gauss Theorem and because F is linear on each triangle (and so ∇φ(|∇I(x)|) is constant), Z Z 1 1 ∇S · (φ(|∇S I(x)|)∇S I(x))dx = φ(|∇S I(x)|)∇S I(x) · nu,v dl A A A ∂A Z 1 X = φ(|∇S I(x)|)∇S I(x) · nu,v dl A ∂A∩T Ti ∈A
4
=
1 X φ(|∇S I(x)|)∇S I(x) · [Xi − Xi+1 ]⊥ I∂A∩Ti (x) A Ti ∈A
where the triangle Ti has p, pi , pi+1 as vertex; X, Xi , Xi+1 are the correspondences of the vertex in the (u, v) space; and I∂A∩Ti (x) stands for the indicating function over the set ∂A ∩ Ti . pi+1 pi
θi
γi Ti
p
βi
A
Fig. 2.
Let us now take a fixed triangle Ti . Let Bl , l = 1, 2, 3 be the linear basis functions over the triangle Ti . Because of (B1 + B2 + B3 )(u, v) = 1, ∇S I(x) = I(p)∇S B1 + I(pi )∇S B2 + I(pi+1 )∇S B3 = (I(pi ) − I(p))∇S B2 (u, v) + (I(pi+1 ) − I(p))∇S B3 (u, v) 1 = [(I(pi ) − I(p))(X − Xi+1 )⊥ 2ATi + (I(pi+1 ) − I(p))(Xi − X)⊥ )] and we found 1 (X − Xi+1 )⊥ 2ATi 1 ∇B3 (u, v) = (Xi − X)⊥ 2ATi ∇B2 (u, v) =
where ATi is the area of the triangle Ti . Using the fact that ATi is proportional to the sine of any angle of the triangle, ∇S I(x) · [Xi − Xi+1 ]⊥ = [cot(pi+1 )(I(pi ) − I(p)) + cot(pi )(I(pi+1 ) − I(p))]
5
So we get X φ(|∇\ S I(p)|) ∆\ { φ(|∇\ g I(p) = S I(p)|) · A Ti ∈A
[cot θi (I(pi ) − I(p)) + cot γi (I(pi+1 ) − I(p))] } for x in the triangle Ti , where θi is the internal angle of the node pi+1 and γi is the internal angle of the node pi , and |∇\ S I(p)| =
1 X |(I(pi ) − I(p))(p − pi+1 ) + (I(pi+1 ) − I(p))(pi − p)| 2AT Ti ∈A
Note that the approximation is independent from the parameterization chosen, and it can be used to approximate other differential operators of the form ∇S · (φ(|∇S I(p)|)∇S I(p))
(3)
for φ being a function φ : IR → IR. Remark that if we take β → 0, we recover the expression of [21] which approximates the Laplace Beltrami operator in order to perform isotropic smoothing: lim ∆\ g I(p) = lim
β→0
β→0
X 1 1 [cot θi (I(pi ) {p 2 1 + β|∇u,v I(x)| Ti ∈A 1 + β|∇u,v I(x)|2
p
A
− I(p)) + cot γi (I(pi+1 ) − I(p))] } 1 X = [cot θi (I(pi ) − I(p)) + cot γi (I(pi+1 ) − I(p))] A Ti ∈A
which is the expression found in [21] as ∆\ u,v I(p), that can be written as ∆\ u,v I(p) =
X
wi (I(pi ) − I(p))
Ti ∈A
with
wi =
cot θi + cot βi AT i
using some common trigonometric identities.
4 4.1
Examples and Results Implementation Details
The implementation was done using the 3D visualization package developed at (Hidden for anonymous review). We compute the value of Ipn , the value of I on the vertex p of Σ at the nth iteration based on the values of Ipn−1 , i pi , .., pnb of neigh being the vertex neighboring p.
6
First, we compute the cot θi ,cot γ, and ATi for each node p only once. Then, for each iteration n, we actualize the flow in this manner: Ipn = Ipn−1 + dt∆\ g I(p) X φ(|∇\ S I(p)|) ∆\ { φ(|∇\ g I(p) = S I(p)|) cot θi (I(pi ) − I(p)) A Ti ∈A
+ cot γi (I(pi+1 ) − I(p)) } 4.2
Examples
For the first example (fig. 3) we use a triangulation from the Stanford’s surface data base with the Japanese word for peace as the data function, and we added a gaussian noise with σ = 40. The first stage of the algorithm takes about 2 minutes to be computed, and the iterations a few more minutes, (from 1 to 5 depending on the value of β). For the second example example (fig. 5) a slice of cortex is taken from our database. The scalar data is a variable that needs to be regularized as an intermediate step to obtain a map of the visual areas on the cortex. The scalar data is naturally noisy because of the measurement errors. We already have from the beginning a stair-casing effect due to the differences of resolution between the anatomic and functional brain images. Even though we have no access to the original non-noisy image, we can select a region where we know that the original image has null variance. We can do this because we know that there should be no electric impulsions outside a certain region of the cortex. To compare the performance with different values of β, we applied the flow until the variance measurement in the selected region reaches a fixed threshold. Then we observed how different values of β act on discontinuities on the image for the same amount of noise reduction on the selected region.
5
Conclusions
We have shown a method to compute the Beltrami flow for scalar functions defined on triangulated manifolds using a local approximation of the operator. This procedure can be extended to any operator of the form (3). We have also illustrated the utility of this approach showing synthetic and real examples.
References 1. M. Bertalm´ıo, L. T. Cheng, S. Osher and G. Sapiro, “Variational Problems and Partial Differential Equations on Implicit Surfaces”, Journal of Computational Physics 174 (2001) 759-780.
7
Original noisy image
β = 0 (Isotropic diffusion)
β = 0.01
β = 1 (Anisotropic diffusion)
Fig. 3. Note the stair-casing effect on the anisotropic diffusion regularization and the blurry image obtained with the isotropic one
8
Fig. 4. Triangulation used: 18979 vertex, 37954 triangles
2. T. Chan and J. Shen, “Variational restoration of non-flat image features: Models and algorithms”, SIAM J. Appl. Math., 61 (2000) 1338-1361. 3. R. Kimmel and R. Malladi and N. Sochen, “Images as Embedding Maps and Minimal Surfaces: Movies, Color, Texture, and Volumetric Medical Images”, International Journal of Computer Vision 39(2) (2000) 111-129. 4. N. Sochen and R. Kimmel, “Stereographic Orientation Diffusion”, in proceedings of the 4th Int. Conf. on Scale-Space, Vancouver Canada, October 2001. 5. A. M. Polyakov, “Quantum geometry of bosonic strings”, Physics Letters, 103B (1981) 207-210. 6. L. Rudin, S. Osher and E. Fatemi, “ Non Linear Total Variation Based Noise Removal Algorithms”, Physica D 60 (1992) 259-268. 7. N. Sochen and R. Kimmel and A. M. Bruckstein, “Diffusions and confusions in signal and image processing”, accepted to Journal of Mathematical Imaging and Vision. 8. N. Sochen and R. Kimmel and R. Malladi, “From high energy physics to low level vision”, Report, LBNL, UC Berkeley, LBNL 39243, August, Presented in ONR workshop, UCLA, Sept. 5 1996. 9. N. Sochen and R. Kimmel and R. Malladi, “A general framework for low level vision”, IEEE Trans. on Image Processing, 7 (1998) 310-318. 10. N. Sochen and Y. Y. Zeevi, “Representation of colored images by manifolds embedded in higher dimensional non-Euclidean space”, Proc. IEEE ICIP’98, Chicago, 1998. 11. N. Sochen and Y. Y. Zeevi, “Representation of images by surfaces embedded in higher dimensional non-Euclidean space”, 4th International Conference on Mathematical Methods for Curves and Surfaces, Lillehamer, Norway, July 1997. 12. B. Tang and G. Sapiro and V. Caselles, “Diffusion of General Data on Non-Flat Manifold via Harmonic Map Theory: The Direction Diffusion Case”, International Journal on Computer Vision 36(2) 149-161, 2000.
9
Original noisy image
Area of interest
β = 0 (Isotropic diffusion)
β = 0.0001
β = 0.1
β = 1 (Anisotropic diffusion)
Fig. 5. Eccentricity image of the visual areas on the cortex
10 13. F. Memoli and G. Sapiro and S. Osher, “Solving Variational Problems and Partial Differential Equations, Mapping into General Target Manifolds”, January 2002 UCLA CAM Technical Report (02-04). 14. L. Cheng, P. Burchard, B. Merriman and S. Osher, “Motion of Curves Constrained on Surfaces Using a Level Set Approach”, September 2000 UCLA CAM Technical Report (00-32). 15. M. Clerc, O. Faugeras, R. Keriven, J. Kybic, T. Papadopoulo, “A level set method for the inverse EEG/MEG problem”, Poster No.: 10133, NeuroImage Human Brain Mapping 2002 Meeting 16. P. Kornprobst, R. Deriche, and G. Aubert. “Nonlinear operators in image restoration”. In Proceedings of the International Conference on Computer Vision and Pattern Recognition - pages 325-331.IEEE Computer Society. San Juan, Puerto Rico, June 1997. 17. G. Sapiro, “Geometric Partial Differential Equations and Image Analysis”, Cambridge University Press, January 2001. 18. N. Sochen, R. Deriche, and L. Lopez-Perez, “Variational Beltrami Flows Over Manifolds”, IEEE ICIP 2003, Barcelone 2003. 19. N. Sochen, R. Deriche, and L. Lopez-Perez, “The Beltrami Flow over Implicit Manifolds”, ICCV 2003. 20. N. Sochen, R. Deriche, and L. Lopez-Perez, “Variational Beltrami Flows Over Manifolds”, INRIA Resarch Report 4897, June 2003. 21. Moo K. Chung, Keith J. Worsley, Jonathan Taylor, Jim Ramsay, Steve Robbins, Alan C. Evans, “Diffusion Smoothing on the Cortical Surface”, Human Brain Mapping 2001 Conference. 22. M. Meyer, M. Desbrun, P. Schr¨ oder, and A. H. Barr, Discrete differential-geometry operators for triangulated 2-manifolds. In Hans-Christian Hege and Konrad Polthier, editors, Visualisation and Mathematics III, pages 35-57. Springer-Verlag, Heidelberg, 2003.
