# Processing Manifold-Valued Images

By Ronny Bergmann, Friederike Laus, Johannes Persch, and Gabriele Steidl

The mathematical notion of a manifold dates back to 1828, when Carl Friedrich Gauss established an important invariance property of surfaces while proving his Theorema Egregium. In his habilitation lecture in 1854, Bernhard Riemann intrinsically extended Gauss’s theory to manifolds of arbitrary dimension, such that they are not dependent upon the embedding in higher dimensional spaces. This is now called a Riemannian manifold. Modern image acquisition methods are able to capture information that is no longer restricted to Euclidean spaces but can be manifold-valued. Such imaging methods include the following:

- Interferometric Synthetic Aperture Radar (InSAR), used in geodesy and remote sensing where each image pixel is on the circle \(\mathbb S^1\)
- Diffusion-tensor magnetic resonance imaging (DT-MRI), which produces images with values in the manifold of symmetric, positive, definite \(3\times 3\) matrices SPD(3), thus mapping directional diffusion processes of molecules—mainly water—in biological tissues
- Electron Backscatter Diffraction (EBSD), which analyzes crystalline materials where the images have pixels in the group of rotations SO(3) modulo the crystal’s symmetry group.

Images produced by these techniques are depicted in Figure 1. Furthermore, manifold-valued images arise when working in color spaces that are more adapted to human color perception than the RGB space, such as HSI (hue, saturation, intensity) or CB (chromaticity, brightness). Here, we deal with the (product)-manifolds \(\mathbb{S}^1\times \mathbb{R}\times \mathbb{R}\) and \(\mathbb{S}^2\times \mathbb{R}\), respectively.

**Figure 1.** Manifold-valued images acquired with different devices, yielding data values given on the following: **1a.** The circle 𝕊^{1}, colored using the hue [11]. **1b.** The manifold of symmetric positive definite 3x3 matrices, illustrated using their eigenvalues and eigenvectors to draw an ellipsoid [3]. **1c.** Orientations SO(3), where the orientation modulo the phase is mapped onto a colorized sphere. 1c courtesy of the Institute of Materials Science and Engineering at the University of Kaiserslautern.

Processing manifold-valued signals and images proposes new challenges in image processing that affect classical tasks like denoising, inpainting, and segmentation. In the real-valued case, variational methods—convex optimization in particular—are well adapted to such large-scale problems. Other successful approaches include nonlocal patch-based methods, which rely on an image self-similarity assumption, and recently-developed learning methods. What follows offers a brief overview of recent results in manifold-valued image processing, obtained by attempting to generalize results from the real to manifold-valued setting.

Tiny inaccuracies that result in noisy data are inevitable when measurements are taken, regardless of whether the measurements are manifold-valued. Modeling an image as a realization of some random variable, whose distribution reflects the circumstances under which the image is taken, is a common method to account for this randomness. In the real-valued case, the standard approach assumes that the noise is additive, white, and Gaussian, which is asymptotically justified by the Central Limit Theorem. Much effort has been spent denoising real-valued images corrupted with white Gaussian noise, but the situation completely changes with manifold-valued images. The definition of a Gaussian distribution is not canonical on a manifold, and different attempts—generalizing different characterizing properties of the real-valued normal distribution—appear in the literature.

Current state-of-the-art denoising methods for real-valued images include nonlocal patch-based approaches, in particular the nonlocal Bayes algorithm of Marc Lebrun, Antoni Buades, and Jean-Michel Morel [9], whose idea was reinterpreted and generalized to manifolds in [8]. The proposed estimation procedure relies heavily on the computation of Riemannian centers of mass — counterparts of the classical mean or expectation. Numerical examples demonstrate the excellent denoising performance of the proposed estimator for different manifolds, including \(\mathbb{S}^1\), the sphere \(\mathbb{S}^2\), and SPD(3). Figure 2 illustrates a result obtained for an orientation field, i.e., an \(\mathbb{S}^2\)-valued image. However, these are academic examples, and part of ongoing research is to examine the acquisition-dependent noise models appearing in applications.

**Figure 2.** Denoising of an artificial 𝕊^{2}-valued image by a patch-based method. Image courtesy of [8].

Variational approaches that are not restricted to denoising generate a restored image as a minimizer of some functional of the form

\[\mathcal{J}(u) = \mathcal{D}(u,f) + \lambda \mathcal{R}(u),\quad \lambda>0,\]

where \(\mathcal{D}\) is a data-fitting term measuring the distance to the given data \(f\), \(\mathcal{R}\) is a regularization term—also called a prior—reflecting the assumed properties of an ideal (clean) image, and \(\lambda>0\) is a factor that balances the influence of the regularizer. A usual data-fitting term for real-valued images is the squared Euclidean distance between \(f\) and \(u\), which is naturally replaced by the squared geodesic distance for manifold-valued images. Choosing appropriate regularization terms is usually more involved; selections should ensure a reduction of noise in the minimizer and the preservation of the image structure. In the Euclidean case, the total variation (TV)—proposed by Leonid Rudin, Stanley Osher, and Emad Fatemi [12]—is demonstrably a powerful, edge-preserving, and convex nonsmooth regularizer. Mariano Giaquinta and Domenico Mucci [7] used Cartesian currents to investigate the notation of TV of functions with values on a manifold, while Evgeny Strekalovskiy and Daniel Cremers [13] first applied the technique to phase-valued images.

In addition to first-order derivatives that occur in the classical TV approach, the incorporation of higher-order derivatives into the model to reduce the staircasing effect caused by TV regularization and adapt to specific applications is also desirable. The definition of second-order spatial differences is not straightforward for manifold-valued images. Following the idea that the second-order difference term in the Euclidean setting reads \(\|x - 2y + z\| = 2 \| \frac12 (x+z) - y\|\), [2] provides a definition on manifolds that uses the geodesic distance from the midpoints of the geodesics connecting \(x\) with \(z\) to \(y\).

One can apply Riemannian optimization methods to compute a minimizer of the resulting functionals \(\mathcal{J}\). These intrinsic methods are often very efficient, since they exploit the underlying geometric structure of the manifold [10]. Various methods have been proposed for smooth functions \(\mathcal{J}\), ranging from simple gradient descents on manifolds to more sophisticated trust region methods. Recently, so-called half-quadratic minimization methods, belonging to the group of quasi-Newton methods and covering (for example) iteratively re-weighted least squares methods, were generalized to manifold-valued images [3].

Unfortunately, the TV regularization term is not differentiable. However, in the Euclidean setting it is convex, meaning that convex analysis tools—including powerful algorithms based on duality theory—are applicable. A prominent example is the alternating direction method of multipliers (ADMM), which is equivalent to the Douglas-Rachford algorithm. Proximal mappings, which one can efficiently compute for special priors appearing in Euclidean image processing tasks, are a central ingredient of these algorithms.

**Figure 3.** Denoising and inpainting of manifold-valued images by a variational method. 3a and 3b courtesy of [2], 3c and 3d courtesy of [3].

Several efforts have recently attempted to translate concepts from convex analysis to manifolds. One can establish a certain theory of convex functions on Hadamard manifolds, i.e., complete, simply-connected Riemannian manifolds of nonpositive sectional curvature, as (for example) or hyperbolic spaces. In particular, one can introduce the (inexact) cyclic proximal point algorithm on these manifolds [1], a method that was also used to minimize the functional with first- and second-order differences [2]. Efficient computation of the proximal mapping of second-order differences utilizes the machinery of Jacobi fields and is restricted to symmetric spaces. Figure 3 shows two denoising results employing first- and second-order differences. Symmetric spaces are characterized by the property that geodesic reflections at points are isometries. Since the classical Douglas-Rachford algorithm relies on point reflections, it was natural to extend this algorithm to symmetric Hadamard manifolds [4]. However, Hadamard spaces do not embody *all* the nice properties of convex analysis. For example, reflections at convex sets are not nonexpansive in general, and although the parallel Douglas-Rachford algorithm shows a very good numerical performance, existing theoretical convergence results remain limited to manifolds with constant nonpositive curvature.

Finally, an important aspect when working with real data is the practical implementability of the developed methods and algorithms. In the spirit of reproducible research, several groups provide their software and toolboxes, e.g., the manifold-valued image restoration toolbox MVIRT, the manifold optimization Manopt package [5], and toolboxes focusing on specific manifolds, such as the MTEX toolbox for EBSD images. All of these packages are available in MATLAB, making it possible to test one’s own ideas and enter this active field of research.

**References**

[1] Bačák, M. (2014). *Convex Analysis and Optimization in Hadamard Spaces*. Vol. 22. Series in Nonlinear Analysis and Applications. Berlin, Germany: De Gruyter.

[2] Bačák, M., Bergmann, R., Steidl, G., & Weinmann, A. (2016). A second order non-smooth variational model for restoring manifold-valued images. *SIAM Journal on Scientific Computing, 38*(1), 567-597.

[3] Bergmann, R., Chan, R.H., Hielscher, R., Persch, J., & Steidl, G. (2016). Restoration of manifold-valued images by half-quadratic minimization. *Inverse Problems and Imaging, 10*(2), 281-304.

[4] Bergmann, R., Persch, J., & Steidl, G. (2016). A parallel Douglas-Rachford algorithm for restoring images with values in symmetric Hadamard manifolds. *SIAM Journal on Imaging Sciences, 9*(3), 901-937.

[5] Boumal, N., Mishra, B., Absil, P.A., & Sepulchre, R. (2014). Manopt, a Matlab toolbox for optimization on manifolds. *Journal of Machine Learning Research, 15*, 1455-1459.

[6] Cook, P.A., Bai, Y., Nedjati-Gilani, S., Seunarine, K.K., Hall, M.G., Parker, G.J., & Alexander, D.C. (2006). Camino: Open-source diffusion-MRI reconstruction and processing. In *Proceedings of the International Society for Magnetic Resonance in Medicine Scientific Meeting and Exhibition*. (pp. 2759). Seattle, WA: International Society for Magnetic Resonance in Medicine.

[7] Giaquinta, M., & Mucci, D. (2007). Maps of bounded variation with values into a manifold: total variation and relaxed energy. *Pure Applied Mathematics Quarterly, 3*(2), 513-538.

[8] Laus, F., Nikolova, M., Persch, J., & Steidl, G. (2007). A nonlocal denoising algorithm for manifold-valued images using second order statistics. *SIAM Journal on Imaging Sciences, 10*(1), 416-448.

[9] Lebrun, M., Buades, A., & Morel, J.-M. (2013). A nonlocal Bayesian image denoising algorithm. *SIAM Journal on Imaging Sciences, 6*(3), 1665-1688.

[10] Ring W., & Wirth, B. (2012). Optimization methods on Riemannian manifolds and their application to shape space. *SIAM Journal on Optimization, 22*(2), 596-627.

[11] Rocca, F., Prati, C., & Guarnieri, A.M. (1997). Possibilities and limits of SAR interferometry. *Image Processing Techniques*, 15-26.

[12] Rudin, L.I., Osher, S., & Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms. *Physica D: Nonlinear Phenomena, 60*(1), 259-268.

[13] Strekalovskiy, E. & Cremers, D. (2013). Total cyclic variation and generalizations. *Journal of Mathematical Imaging and Vision, 47*(3), 258-277.

Ronny Bergmann received his Ph.D. in mathematics from the University of Lübeck in 2013. He is currently a postdoctoral research fellow at the University of Kaiserslautern. Friederike Laus and Johannes Persch received their M.Sc. in mathematics from the University of Kaiserslautern in 2015 and are currently Ph.D. students under the supervision of Gabriele Steidl. Gabriele Steidl earned her Ph.D. and habilitation in mathematics from the University of Rostock. She is currently head of the Mathematical Image Processing Group at the University of Kaiserslautern and a consultant at the Fraunhofer Institute for Industrial Mathematics.