Recent Advances in Denoising of Manifold-Valued Images

TL;DR

This paper reviews recent variational model-based methods and algorithms for denoising manifold-valued images, emphasizing practical implementations and software tools for real-world applications.

Contribution

It provides a comprehensive overview of recent advances in variational denoising techniques for manifold-valued data, including algorithm generalizations and software resources.

Findings

Algorithms like subgradient, half-quadratic, cyclic proximal point, and Douglas-Rachford are adapted for manifolds.

Software packages such as Manopt and MVIRT facilitate practical implementation.

Recent methods improve denoising quality for manifold-valued signals and images.

Abstract

Modern signal and image acquisition systems are able to capture data that is no longer real-valued, but may take values on a manifold. However, whenever measurements are taken, no matter whether manifold-valued or not, there occur tiny inaccuracies, which result in noisy data. In this chapter, we review recent advances in denoising of manifold-valued signals and images, where we restrict our attention to variational models and appropriate minimization algorithms. The algorithms are either classical as the subgradient algorithm or generalizations of the half-quadratic minimization method, the cyclic proximal point algorithm, and the Douglas-Rachford algorithm to manifolds. An important aspect when dealing with real-world data is the practical implementation. Here several groups provide software and toolboxes as the Manifold Optimization (Manopt) package and the manifold-valued image restoration toolbox (MVIRT).

Figures

Tables

Table 1. Table 1: Functions φ 𝜑 \varphi fulfilling the assumptions of Proposition 4.2 .

	$\phi \left(t\right)$	$s\left(t\right)$
${\phi }_1\left(t\right)$	$\sqrt{t^2+{\epsilon }^2}$	$\frac{1}{2\sqrt{t^2+{\epsilon }^2}}$
${\phi }_2\left(t\right)$
${\phi }_3\left(t\right)$	$1-\exp \left(-{\epsilon }^2{t}^2\right)$	${\epsilon }^2\exp \left(-{\epsilon }^2{t}^2\right)$

Equations

$$\operatorname{dist}^{2}(x,y)\coloneqq\Big{(}\sum_{j=1}^{n}\operatorname{dist}^{2}(x_{j},y_{j})\Big{)}^{\frac{1}{2}}.$$

$${\exp }_x(\xi ):={\gamma }_{x;\xi }(1).$$

$${\log }_x:={\exp }_x^{-1}:\mathcal{M}\backslash \mathcal{C}(x)\to T_x\mathcal{M}.$$

$$\mathrm{dist}(x,y)=⟨{\log }_x(y),{\log }_x(y){⟩}_x^{\frac{1}{2}}=\Vert {\log }_x(y){\Vert }_x.$$

$$DF(x):T_x\mathcal{M}\to T_{F(x)}\mathcal{N},\xi \mapsto DF(x)[\xi ],$$

$$D(G\circ F)(x)[\xi]=DG\bigl{(}F(x)\bigr{)}\bigl{[}DF(x)[\xi]\bigr{]}.$$

$$⟨{\mathrm{grad}}_{\mathcal{M}}f(x),\xi {⟩}_x:=Df(x)[\xi ]$$

$$R_{x}(x)=x\quad\text{and}\quad D\bigl{(}R_{x}(x)\bigr{)}=-I.$$

$$P_{x\to y}\xi :=X(1),$$

$$P_{x\to y}^{\mathrm{P}}(\xi)\coloneqq-\log_{y}\Bigl{(}\gamma\Bigl{(}\exp_{x}(\xi),\gamma\bigl{(}x,y;\tfrac{1}{2}\bigr{)};2\Bigr{)}\Bigr{)}\in T_{y}\mathcal{M},\$$

$$P_{x\to y}^{\mathrm{S}}(\xi)\coloneqq\log_{y}\Bigl{(}\gamma\Bigl{(}x,\gamma\bigl{(}y,\exp_{x}(\xi);\tfrac{1}{2}\bigr{)};2\Bigr{)}\Bigr{)}\in T_{y}\mathcal{M}.$$

$$DF(x)[\xi ]=\sum _{k=1}^d⟨\xi ,{\Xi }_k(0){⟩}_x\alpha ({\kappa }_k){\Xi }_k(T),$$

$$\alpha (\kappa ):=\{\begin{array}{ll}\cosh (\sqrt{-\kappa })& \kappa <0,\\ 1& \kappa =0,\\ \cos (\sqrt{\kappa })& \kappa >0.\end{array}$$

$$\alpha (\kappa ):=\{\begin{array}{ll}-\sqrt{-\kappa }\frac{\cosh (\sqrt{-\kappa })}{\sinh (\sqrt{-\kappa })}& \kappa <0,\\ -1& \kappa =0,\\ -\sqrt{\kappa }\frac{\cos (\sqrt{\kappa })}{\sin (\sqrt{\kappa })}& \kappa >0.\end{array}$$

$$\alpha (\kappa ):=\{\begin{array}{ll}\frac{\sqrt{-\kappa }}{\sinh (\sqrt{-\kappa })}& \kappa <0,\\ 1& \kappa =0,\\ \frac{\sqrt{\kappa }}{\sin (\sqrt{\kappa })}& \kappa >0.\end{array}$$

$$\alpha(\kappa)\coloneqq\begin{cases}\frac{\sinh\bigl{(}\sqrt{-\kappa}(1-\tau)\bigr{)}}{\sinh(\sqrt{-\kappa})}\quad&\kappa<0,\\ 1-\tau\quad&\kappa=0,\\ \frac{\sin\bigl{(}\sqrt{\kappa}(1-\tau)\bigr{)}}{\sin(\sqrt{\kappa})}\quad&\kappa>0.\end{cases}$$

$$\alpha (\kappa ):=\{\begin{array}{ll}\frac{\sinh (\sqrt{-\kappa }\tau )}{\sinh (\sqrt{-\kappa })}& \kappa <0,\\ \tau & \kappa =0,\\ \frac{\sin (\sqrt{\kappa }\tau )}{\sin (\sqrt{\kappa })}& \kappa >0.\end{array}$$

$$\alpha (\kappa )=\{\begin{array}{ll}\frac{\sinh (\sqrt{-\kappa })}{\sqrt{-\kappa }}& \kappa <0,\\ 1& \kappa =0,\\ \frac{\sin (\sqrt{\kappa })}{\sqrt{\kappa }}& \kappa >0,\end{array}$$

$$(DF{)}^{\ast }(x)[w]=\sum _{k=1}^d⟨w,{\Xi }_k{⟩}_{F(x)}{\alpha }_k{\xi }_k,w\in T_{F(x)}\mathcal{M}.$$

$$\varphi\bigl{(}\gamma_{\overset{\frown}{x,y}}(t)\bigr{)}\leq(1-t)\varphi(x)+t\varphi(y),$$

$$\varphi\bigl{(}\gamma_{\overset{\frown}{x,y}}(t)\bigr{)}\leq(1-t)\varphi(x)+t\varphi(y)-\kappa t(1-t)\operatorname{dist}^{2}(x,y).$$

$$\mathrm{dom}(\phi ):=\{x\in \mathcal{H}:\phi (x)<\infty \},$$

$$\mathcal{J}(u):=\mathcal{D}(u;f)+\alpha \mathcal{R}(u),$$

$$\mathcal{D}(u;f):=\frac{1}{2}\sum _{i\in G}\mathrm{dist}(u_i,f_i{)}^2.$$

$${\nabla }_x{u}_i:=\{\begin{array}{ll}\text{if}{\log }_{u_i}{u}_{i+(1,0)}& i+(1,0)\in G,\\ 0& \text{otherwise},\end{array}$$

$$\mathrm{TV}(u)$$

$$\displaystyle=\sum_{i\in\operatorname{{\mathcal{G}}}}\biggl{(}\sum_{j\in\operatorname{\mathcal{N}}(i)}\operatorname{dist}(u_{i},u_{j})^{p}\biggr{)}^{\frac{1}{p}},\quad p\in\{1,2\},$$

$${\mathrm{TV}}_{\phi }(u):=\sum _{i\in G}\sum _{j\in N(i)}\phi \left(\mathrm{dist}(u_i,u_j)\right),$$

$$\operatorname{TV}_{\varphi}(u)\coloneqq\sum_{i\in\operatorname{{\mathcal{G}}}}\varphi\biggl{(}\Bigl{(}\sum_{j\in\operatorname{\mathcal{N}}(i)}\operatorname{dist}(u_{i},u_{j})^{2}\Bigr{)}^{\tfrac{1}{2}}\biggr{)}.$$

$${\mathrm{d}}_2(x,y,z):=\underset{c\in {\mathcal{C}}_{x,z}}{\min }\mathrm{dist}(c,y),$$

$${\mathrm{d}}_{1,1}(x,y,z,w):=\underset{c\in {\mathcal{C}}_{x,z},\tilde{c}\in {\mathcal{C}}_{y,w}}{\min }\mathrm{dist}(c,\tilde{c}).$$

Full text

1 INTRODUCTION

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 Eregium. In his habilitation lecture in 1854, Bernhard Riemann intrinsically extended Gauss’s theory making manifolds independent of their embedding in higher dimensional spaces. This is now called a Riemannian manifold. Nowadays modern signal and image acquisition methods are able to capture information that is no longer restricted to Euclidean spaces but can be manifold-valued. Sophisticated models for human color perception involve non-Euclidean settings. Moreover, it is often advantageous to model information from large data as points on a certain manifold. Here are some typical examples.

Interferometric Synthetic Aperture Radar (InSAR). An InSAR image is obtained by calculating the phase difference of two Synthetic Aperture Radar (SAR) images of an area taken at different positions or times [23]. It can be used to measure elevation when the measurements are taken at the same time or to detect millimeter-scale deformations over days or years. It has applications in the geophysical monitoring of natural hazards, for example earthquakes, volcanoes and landslides, and in structural engineering, in particular monitoring of subsidence and structural stability. InSAR produces phase-valued images, i.e., in each pixel the measurement lies on the circle $\mathbb{S}^{1}$, see Figure 1 (left).

Image Color Spaces. The usual RGB color space has a vector space structure, but there exist models which are physically better suited for the human color perception as they encode luminance independent of color. Examples are the hue-saturation-value (HSV) color space or the related Ich space, where the hue has its value on the circle $\mathbb{S}^{1}$, as well as the chromaticity-brightness (CB) color space, where the chromaticity has values on the positive octant of the sphere $\mathbb{S}^{2}$. For a recent geometric model of brightness perception we refer to [8].

Electron Backscatter Diffraction (EBSD). EBSD is a microstructural crystallography characterization technique used to analyze the microscopic structure of polycrystalline materials such as metals and minerals [7]. Each point of a specimen is radiated by an electron beam and the diffraction pattern is measured, which gives information on the phase and the crystal orientation, a value on the rotation group $\operatorname{SO}(3)$. Regions of similar orientation are called grains. Material scientists are interested in the grain structure of the specimen, as it affects macroscopic properties such as ductility, electrical and lifetime properties. Since the atomic structure of a crystal is invariant under the symmetry of its atomic lattice, i.e., a symmetry group $S\subset\operatorname{SO}(3)$, the images obtained by EBSD have pixel values in $\operatorname{SO}(3)/S$, see Figure 1 (right). The software MTEX [6] is designed for processing EBSD data.

Diffusion Tensor Magnetic Resonance Imaging (DT-MRI). In DT-MRI the diffusion of water molecules perpendicular to a magnetic field is measured in biological tissue. Taking at least six different data sets measured with different magnetic fields, a DT-MRI image with values in the manifold $\mathcal{P}(3)$ of symmetric positive definite $3\times 3$ matrices is computed at each pixel, see Figure 1 (middle). The diffusion of water is influenced by the structure of the tissue. Hence, the knowledge about the diffusion can be used to distinguish between diseased and healthy parts. DT-MRI is a non-invasive and in-vivo technique, which is extensively used in neurology, but can also be applied to detect defects in other tissue, like muscles.

Covariance Matrices in Texture Analysis and Brain Computer Interfaces. Textures form a special class of images and appear at the same time as an important image feature. There are different ways to model textures, one possibility is to encode the (local) dependence structure as a covariance matrix. This has been used, e.g., in [57], where textures are characterized based on the empirical covariance of certain features, for example color intensities and derivatives of different orders. Under the assumption that the empirical covariance matrices are non-degenerated we are again faced with the manifold $\mathcal{P}(d)$ of symmetric positive definite $d\times d$ matrices, where $d$ equals the number of features. In Brain Computer Interfaces (BCI) approaches, EEG curves related to different human activities are measured at different brain positions at the same time. Their covariance matrices can be used to analyze the corresponding activity, see, e.g. [64].

Beyond these applications, images with values in $\mathbb{S}^{2}$ appear when dealing with 3D directional information [41] as well as in the analysis of liquid crystals [2]. The rotation group $\operatorname{SO}(3)$ and the special Euclidean group $\operatorname{SE}(3)$ are considered in tracking and (scene) motion analysis [54, 58].

Processing manifold-valued signals and images is a new challenge that affects classical tasks like image restoration (denoising, inpainting), segmentation and clustering [12], registration and large deformation diffeomorphic mapping (LDDMM) or metamorphosis between different images, see, e.g., [47, 48, 65].

This chapter focuses on variational denoising methods for manifold-valued images. The simplest idea is to embed the manifold into the Euclidean space and to apply the Euclidean models with the constraint that the image values have to lie in a manifold. Recall that by Whitney’s theorem every smooth $d$-dimensional manifold can be smoothly embedded into an Euclidean space of dimension $2d$. Such an approach was given e.g., in [55]. The advantage is that optimization algorithms in Euclidean spaces can be applied, where the models are in general non-convex due to the constraints. So-called lifting schemes were proposed for circular-valued data in [27] and for more general manifold-valued images in [44]. There, the problem is reformulated as a multilabel optimization problem which is approached using convex relaxation techniques. We also like to mention that current state-of-the-art methods for denoising of real-valued images include stochastic nonlocal patch-based approaches, as, e.g., the nonlocal Bayes’ algorithm [43]. A generalization of this minimum mean square estimator (MMSE) based method to manifold-valued images was proposed in [42].

This chapter deals with spatially discrete, intrinsic models and algorithms. Note that in [35, 36], the notion of total variation of spatially continuous functions having their values on a manifold was investigated where the authors apply the theory of Cartesian currents. For a spatial continuous setting, the reader may also consult the recent paper [25].

An important aspect when working with real data is the practical implementation of the developed methods and algorithms. In the spirit of reproducible research, several groups provide their software and toolboxes, e.g., the Manifold Optimization (Manopt) package [18] and the manifold-valued image restoration toolbox (MVIRT [9].

The outline of this chapter is as follows: Starting with the necessary preliminaries in Section 2, we review several denoising models for manifold-valued images in Section 3. Appropriate minimization algorithms are discussed in Section 4. Numerical examples in Section 5 illustrate the proposed methods.

2 Preliminaries on Riemannian Manifolds

2.1 General Notation

Throughout this chapter, let $\mathcal{M}$ be a connected, complete $d$-dimensional Riemannian manifold. By $T_{x}\mathcal{M}$ we denote the tangent space of $\mathcal{M}$ at $x\in\mathcal{M}$ with the Riemannian metric $\langle\cdot,\cdot\rangle_{x}$ and corresponding norm $\|\cdot\|_{x}$. Further, let $T\mathcal{M}$ be the tangent bundle of $\mathcal{M}$. By $\operatorname{dist}\colon\mathcal{M}\times\mathcal{M}\rightarrow\mathbb{R}_{\geq 0}$ we denote the geodesic distance on $\mathcal{M}$. Let $\mathcal{M}^{n}$ be the product or $n$-fold power manifold with product distance

[TABLE]

Let $\gamma_{\overset{\frown}{x,y}}:[0,1]\rightarrow{\mathcal{M}}$ be a (not necessarily shortest) geodesic connecting $x,y\in\mathcal{M}$. Further, we apply the notation $\gamma_{{x;\xi}}$ to characterize the geodesics by its starting point $\gamma_{{x;\xi}}(0)=x$ and direction $\dot{\gamma}_{{x;\xi}}(0)=\xi\in T_{x}{\mathcal{M}}$. Note that the geodesic $\gamma_{\overset{\frown}{x,y}}$ is unique on manifolds with nonpositive curvature. The exponential map $\exp_{x}\colon T_{x}\mathcal{M}\to\mathcal{M}$ is defined by

[TABLE]

Since $\mathcal{M}$ is connected and complete, we know by the Hopf-Rinow theorem that the exponential map is indeed defined on the whole tangent space. The exponential map realizes a local diffeomorphism from a neighborhood $\mathcal{D}_{T}(0_{x})$ of the origin $0_{x}$ of $T_{x}\mathcal{M}$ into a neighborhood of $x\in\mathcal{M}$. More precisely, extending the geodesic $\gamma_{{x;\xi}}$ from $t=0$ to infinity is either minimizing $\operatorname{dist}(x,\gamma_{{x;\xi}}(t))$ all along or up to a finite time $t_{0}$ and not any longer afterwards. In the latter case, $\gamma_{{x;\xi}}(t_{0})$ is called cut point and the set of all cut points of all geodesics starting from $x$ is the cut locus $\mathcal{C}(x)$. This allows to define the inverse exponential map, also known as logarithmic map as

[TABLE]

Then, the Riemannian distance between $x,y\in\mathcal{M}$, for $y\notin\mathcal{C}(x)$, can be written as

[TABLE]

Let $F\colon\mathcal{M}\to\mathcal{N}$ be a smooth mapping between manifolds and $\xi\in T_{x}\mathcal{M}$. The linear mapping

[TABLE]

is called differential of $F$ at $x\in\mathcal{M}$. Let $F\colon\mathcal{M}_{1}\to\mathcal{M}_{2}$ and $G\colon\mathcal{M}_{2}\to\mathcal{M}_{3}$ be two smooth mappings. Then the differential of their concatenation $G\circ F$ applied to $\xi\in T_{x}\mathcal{M}_{1}$ is given by the chain rule

[TABLE]

For a function $f\colon\mathcal{M}\to\mathbb{R}$, the Riemannian gradient $\operatorname{grad}_{\mathcal{M}}$ is defined by

[TABLE]

for all $\xi\in T_{x}\mathcal{M}$. A mapping $R_{p}\colon\mathcal{M}\rightarrow\mathcal{M}$ is called geodesic reflection at $x\in\mathcal{M}$, if

[TABLE]

For $\mathcal{M}=\mathbb{R}^{n}$ we simply have $R_{p}(x)=2p-x$. A connected Riemannian manifold $\mathcal{M}$ is (globally) symmetric if the geodesic reflection at any point $x\in\mathcal{M}$ is an isometry of $\mathcal{M}$, i.e. $\operatorname{dist}\left(R_{p}(x),R_{p}(y)\right)=\operatorname{dist}(x,y)$ for all $x,y\in\mathcal{M}$. All manifolds considered in this chapter are symmetric ones.

Let ${\mathcal{X}}(\mathcal{M})$ be the set of smooth vector fields on $\mathcal{M}$. Given a curve $\gamma\colon[0,1]\to\mathcal{M}$, we denote by ${\mathcal{X}}(\gamma)$ the set of smooth vector fields along $\gamma$, i.e., $X\in{\mathcal{X}}(\gamma)$ is a smooth mapping $X\colon[0,1]\to T\mathcal{M}$ with $X(t)\in T_{\gamma(t)}\mathcal{M}$. A vector field $X\in{\mathcal{X}}(\gamma)$ is called parallel to $\gamma\colon[0,1]\to\mathcal{M}$, if the covariant derivative along $\gamma$ fulfills $\frac{D}{\,\mathrm{d}t}X=0$ for all $t\in[0,1]$. We define the parallel transport of a tangent vector $\xi\in T_{x}\mathcal{M}$ to $T_{y}\mathcal{M}$ by

[TABLE]

where $X\in\mathcal{X}(\gamma_{\overset{\frown}{x,y}})$ is the vector field parallel to a minimizing geodesic $\gamma_{\overset{\frown}{x,y}}$ with $X(0)=\xi$. There exist analytical expressions of the parallel transport for few manifolds as spheres or positive definite matrices. However, the parallel transport can be locally approximated by Schild’s ladder [29, 40] or by the pole ladder [46]. Recently, it was shown that for connected, complete, symmetric Riemannian manifold, the pole ladder coincides with the parallel transport along geodesics [50]. Therefore, we prefer the pole ladder approach. Given $x,y\in\mathcal{M}$, the pole ladder transports $\xi\in T_{x}\mathcal{M}$ to $T_{y}\mathcal{M}$ by

[TABLE]

where we use the notation $\gamma(x,y;t)\coloneqq\gamma_{\overset{\frown}{x,y}}(t)$. For comparison, Schild’s ladder transports as follows:

[TABLE]

Both transport schemes are illustrated in Figure 2.

In our minimization algorithms, we will need the Riemannian gradient of special functions, in particular of those appearing in the pole ladder (10). These gradients can be derived from differentials, see (7), and can be computed for symmetric Riemannian manifolds using the theory of Jacobi fields. The following lemma collects the final results which can be partially found in [5, 20, 28]. For the complete proof we refer to [51].

Lemma 2.1

Let $\mathcal{M}$ be a symmetric Riemannian manifold and $F$ one of the functions i) - v) below with parameter $y$, resp. $u$, together with the coefficient map $\alpha\colon\mathbb{R}\rightarrow\mathbb{R}$ and parameter $T$. Then the differential $DF(x)$ at $x\in\mathcal{M}$ is given for all $\xi\in T_{x}\mathcal{M}$ by

[TABLE]

where $\{\Xi_{k}\}_{k=1}^{d}$ denotes a parallel transported orthogonal frame along the geodesic $\gamma$ with $\gamma(0)=x$ and $\gamma(1)=y$, where $y\coloneqq\exp_{x}(u)$ if $F$ depends on $u$. Further, the frame diagonalizes the Riemannian curvature tensor $R(\cdot,\dot{\gamma})\dot{\gamma}$ with respective eigenvalues $\kappa_{k}$, $k=1,\dots,d$. The functions $F$ and $\alpha$ are given as follows:

1. i)

For $F\coloneqq\exp_{\cdot}(u)$, we have $T=1$, $y\coloneqq\exp_{x}(u)$ and

[TABLE] 2. ii)

For $F\coloneqq\log_{\cdot}(y)$, we have $T=0$ and

[TABLE] 3. iii)

For $F\coloneqq\log_{y}(\cdot)$, we have $T=1$ and

[TABLE] 4. iv)

For $F\coloneqq\gamma_{\overset{\frown}{\cdot,y}}(\tau)$, we have $T=\tau$ and

[TABLE] 5. v)

For $F\coloneqq\gamma_{\overset{\frown}{y,\cdot}}(\tau)$, we have $T=1-\tau$ and

[TABLE] 6. vi)

Finally, we obtain for $F\coloneqq\exp_{x}(\cdot)$ with

[TABLE]

and $T=1$ that the differential $DF(u)$ of $F$ at $u\in T_{x}\mathcal{M}$ is given by (12), where we have to replace $x\in\mathcal{M}$ by $u\in T_{x}\mathcal{M}$ and to set $y\coloneqq\exp_{x}(u)$.

The adjoint operator $(DF)^{*}(x)\colon T_{F(x)}\mathcal{M}\to T_{x}\mathcal{M}$ of (12), which is also needed for the computation of the gradients, is given by

[TABLE]

2.2 Convexity and Hadamard Manifolds

A subset $C\subseteq\mathcal{M}$ is called weakly (strongly) convex if for all points $x,y\in C$ there exists a (unique) geodesic $\gamma_{\overset{\frown}{x,y}}$ of minimal length which is contained entirely in $C$. Let $\varphi\colon\mathcal{M}\supset C\to\mathbb{R}\cup\{\infty\}$ be a real-valued function on a weakly convex set $C$, then $\varphi$ is called convex if

[TABLE]

for all $x,y\in C$. The function $\varphi$ is called strictly convex, if the above equation holds strictly for all $t\in(0,1)$. A function $\varphi\colon\mathcal{M}\supseteq C\to\mathbb{R}\cup\{\infty\}$ is $\kappa$-strongly convex if

[TABLE]

Simply connected, complete Riemannian manifolds of nonpositive sectional curvature are called Hadamard manifolds. We denote them by $\mathcal{H}$. Examples are the manifold of positive definite matrices of fixed size with the affine invariant metric or hyperbolic spaces. An important property of Hadamard manifolds is, that the distance function $\operatorname{dist}(\cdot,\cdot)$ is jointly convex, which makes $\operatorname{dist}^{2}(\cdot,\cdot)$ strictly convex. Further, $\operatorname{dist}^{2}(\cdot,y)$ is a 1-strongly convex function on a Riemannian manifold if and only if $\mathcal{H}$ is a Hadamard manifold. The domain of a function $\varphi\colon\mathcal{H}\to\mathbb{R}\cup\{\infty\}$ is defined by

[TABLE]

in general we work with proper functions, i.e., $\operatorname{dom}(\varphi)\neq\emptyset$. A function $\varphi$ is called lower semi continuous (lsc) if the set $\{x\in\mathcal{H}:\varphi(x)\leq c\}$ is closed for all $c\in\mathbb{R}$. Whenever $\operatorname{dist}(x,x_{0})\to\infty$ for some $x_{0}\in\mathcal{H}$, the function $\varphi$ is called coercive if $\varphi(x)\to\infty$. Concerning minimizers of convex functions, the next theorem summarizes some basic facts. A proper, convex, lsc functions $\varphi\colon\mathcal{H}\rightarrow\mathbb{R}\cup\{\infty\}$ has a minimizer if it is coercive and unique minimizer if it is in addition strongly convex. For further information on more general Hadamard spaces and the basics of convex analysis therein, we refer to [4].

3 Intrinsic Variational Restoration Models

We consider images as mappings from the image grid $\operatorname{{\mathcal{G}}}=\{1,\dots,n_{1}\}\times\{1,\dots,n_{2}\}$ to a Riemannian manifold $\mathcal{M}$. Let $n\coloneqq n_{1}n_{2}$. Given a corrupted image $f\colon\operatorname{{\mathcal{G}}}\to\mathcal{M}$, variational models generate a restored image $u\colon\operatorname{{\mathcal{G}}}\to\mathcal{M}$ as a minimizer of a functional of the form

[TABLE]

where $\mathcal{D}(\cdot;f)$ denotes the data-fitting term, $\alpha>0$ the regularization parameter, and $\mathcal{R}$ the regularization term or prior.

For real-valued images typically a squared Euclidean distance is chosen as data-fitting term. For manifold-valued data, we can just use the squared distance on the manifold

[TABLE]

Depending on the prior knowledge, several regularization terms were proposed in the Euclidean setting. The (discretized) total variation ($\operatorname{TV}$) introduced by Rudin, Osher, and Fatemi [56] is demonstrably a powerful, edge-preserving, non-smooth, and convex regularizer. It sums up the norms of the gradients at the image points. A natural way to define the discrete gradient $\nabla\coloneqq(\nabla_{x},\nabla_{y})^{\mbox{\tiny{T}}}$ on manifold-valued images is as a vector field in the corresponding tangent spaces $(T_{u_{i}}\mathcal{M})^{2}$ with the difference operators

[TABLE]

and similarly for $\nabla_{y}$. Now, the $\operatorname{TV}$ regularizer on manifold-valued images becomes

[TABLE]

with $\operatorname{\mathcal{N}}(i)\coloneqq\{i+(1,0),i+(0,1)\}\cap\operatorname{{\mathcal{G}}}$. Here, $p=1$ is used for the anisotropic and $p=2$ for the isotropic model. This setting was considered in [44, 63]. If $\mathcal{M}=\mathcal{H}$ is an Hadamard manifold, then functional consisting of the data term (24) and the prior (26) has the same properties as its real-valued version, i.e., it is strongly convex and coercive and hence there exists a unique minimizer. The TV functional (26) is not differentiable. To apply minimization algorithms for differentiable functions we can recast it, using an even, differentiable function $\varphi:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$, in the anisotropic case as

[TABLE]

and in the isotropic one as

[TABLE]

Typical functions $\varphi$ are the Huber function and $\varphi(x)\coloneqq\sqrt{x^{2}+\varepsilon^{2}}$ with a small $\varepsilon$.

The minimizers of the TV regularized functionals prefer piecewise constant functions, a behavior called staircasing. To avoid such artifacts second order differences were incorporated into the regularizer. In the manifold-valued setting, we have to find a counterpart of such differences. A generalization of the anisotropic so-called second order $\operatorname{TV}$ term was given in [5]. Observing that in the Euclidean case the absolute second order difference of $x,y,z\in\mathbb{R}^{d}$ can be rewritten as $\lvert x-2y+z\rvert=2\lvert\frac{1}{2}(x+z)-y\rvert$, a counterpart for $x,y,z\in\mathcal{M}$ is defined as

[TABLE]

where $\mathcal{C}_{x,z}$ denotes the set of midpoints $\gamma_{\overset{\frown}{x,z}}(\frac{1}{2})$ of all geodesics connecting $x$ and $z$. Note that the geodesic $\gamma_{\overset{\frown}{x,z}}$ is unique on Hadamard manifolds. Similarly, second order mixed differences were defined for $x,y,z,w\in\mathcal{M}$ as

[TABLE]

We emphasize that the absolute second order difference $\mathrm{d}_{2}$ is not convex in $x$ and $z$ on Hadamard manifolds. Now, we can introduce the absolute value of the second order difference in $x$-direction as

[TABLE]

and similarly in $y$-direction. As absolute value of the mixed differences we use

[TABLE]

and similarly for $\mathrm{d}_{yx}$. Then we define

[TABLE]

where $p=1$ is used for the anisotropic model and $p=2$ for the isotropic model. In the regularizer, the $\operatorname{TV}$ and $\operatorname{TV}_{2}$ terms can appear separately or in a coupled way. Actually, their addition

[TABLE]

was considered in [5, 17]. Alternatively, couplings which generalize the infimal convolution approach [24] to the manifold-valued setting were proposed in [11, 13]. In the Euclidean setting, the infimal convolution is related to the total generalized variation (TGV) approach of Bredies et al. [21]. Recently, Bredies et al. [20] came also up with a TGV model for manifold-valued images, see also [60] for DT-MRI. In the following we present a different TGV approach from [13]. For the relation between both model see [13, Remark 5.1]. In the Euclidean setting, the (discrete) TGV regularizer reads as

[TABLE]

where $\widetilde{\nabla}$ denotes a certain symmetric backward difference operator. For a vector field $\xi=(\xi_{i})_{i\in\operatorname{{\mathcal{G}}}},\ \xi_{i}\in(T_{u_{i}}\mathcal{M})^{2}$, the distance between $\xi$ and $\nabla u$ in the first summand is given by

[TABLE]

To compute the backward differences $\widetilde{\nabla}\xi$ in the second summand, we need to compare tangent vectors from different tangent spaces. To this end, we apply the parallel transport between the tangent spaces via the pole ladder. Since the corresponding expression in (10) contains only exponential and logarithmic maps, the differentials required in the minimization procedure can be calculated using the chain rule and Lemma 2.1. More precisely, we define backward differences of a vector field $\zeta\in T_{u}\mathcal{M}^{n}$ in $x$-direction by

[TABLE]

and similarly in $y$-direction. Then we define

[TABLE]

Note that due to simplified computations, this definition differs slightly from the symmetric arrangement of the backward differences in the Euclidean TGV setting. Now, we can define a TGV regularizer for manifold-valued images as

[TABLE]

4 Minimization Algorithms

To compute a minimizer of our functionals, Riemannian optimization methods can be applied. These intrinsic methods are often very efficient since they exploit the underlying geometric structure of the manifold, see e.g. [1, 52]. For smooth functionals, various methods have been proposed, reaching from simple gradient descents on manifolds to more sophisticated trust region or (quasi) Newton methods. We start by recalling the gradient decent algorithm or more precisely the subgradient algorithm which can also be applied for the minimization of non-differentiable functions.

4.1 Subgradient Descent

The subdifferential of a convex function $\varphi\colon\mathcal{M}\to(-\infty,+\infty]$ at $x\in\operatorname{dom}f$ is defined by

[TABLE]

see, e.g., [45] or [59] for finite functions $\varphi$. For any $x\in\operatorname{int}(\operatorname{dom}\varphi)$, the subdifferential is a nonempty convex and compact set in $T_{x}\mathcal{M}$. If the Riemannian gradient $\operatorname{grad}_{\mathcal{M}}\varphi(x)$ of $\varphi$ in $x\in\mathcal{M}$ exists, then \partial\varphi(x)=\bigl{\{}\operatorname{grad}_{\mathcal{M}}\varphi(x)\bigr{\}}. Further, we see from the definition that $x\in\mathcal{M}$ is a global minimizer of $\varphi$ if and only if $0\in\partial\varphi(x)$.

Let $\varphi\colon\mathcal{M}\to\mathbb{R}$ be a convex function and $x^{(0)}\in\mathcal{M}$ the starting point. Given a sequence $(\tau_{r})_{r\in\mathbb{N}}$ of nonegative numbers, the subgradient algorithm iterates

[TABLE]

We have the following convergence result, see [31].

Theorem 4.1** **(Convergence of subgradient algorithm)

Let $\mathcal{M}$ be a Riemannian manifold with non-negative sectional curvature, $\varphi\colon\mathcal{M}\to\mathbb{R}\cup\{\infty\}$ a convex function which has a minimizer, and $(\tau_{r})_{r\in\mathbb{N}}$ a sequence of positive numbers in $\ell^{2}\backslash\ell^{1}$. Then the sequence $\{x^{(r)}\}_{r\in\mathbb{N}}$ generated by subgradient algorithm converges to a minimizer of $\varphi$.

For manifolds with curvature bounded from below the subgradient algorithm converges if the iterates stay in bounded sets, see [61] or [62].

4.2 Half-Quadratic Minimization

Half-quadratic minimization methods belonging to the group of quasi-Newton methods are efficient minimization algorithms for the functionals with differentiable anisotropic or isotropic regularizers $\operatorname{TV}_{\varphi}$. These methods, which cover iteratively re-weighted least squares methods, were recently generalized to manifold-valued images [10, 37]. There exist additive and multiplicative versions of the method, see [33, 34, 49]. Here, we focus on the multiplicative one for the isotropic $\operatorname{TV}_{\varphi}$ regularizer, i.e., we want to minimize

[TABLE]

We consider the method based on the so-called $c$-transform. Given a function $c\colon\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$, the $c$-transform of a function $\varphi\colon\mathbb{R}\rightarrow\mathbb{R}$ is defined by

[TABLE]

We see immediately that $\varphi(t)+\varphi^{c}(s)\leq c(t,s)$. For $c(t,s)\coloneqq-st$, the function $\varphi^{c}=-(-\varphi)^{*}$ is just the Fenchel transform of $\varphi$. We need the following proposition, see [10].

Proposition 4.2

Let $\varphi\colon\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ be an even, differentiable function and $c(t,s)\coloneqq t^{2}s$.

1. i)

If the function $\Phi(t)\coloneqq-\varphi(\sqrt{t})$ for $t\geq 0$ and $\Phi(t)\coloneqq+\infty$ for $t<0$ is convex, then $\varphi=\varphi^{cc},$ i.e., for $\psi(s)\coloneqq\varphi^{c}(s)$ it holds

[TABLE] 2. ii)

If in addition $\lim_{t\rightarrow\infty}\frac{\varphi(t)}{t^{2}}\rightarrow 0,$ and $\varphi^{\prime}(t)\geq 0$ for $t\geq 0$ and $\varphi^{\prime\prime}(0+)\coloneqq\lim_{t\rightarrow 0+}\frac{\varphi^{\prime}(t)}{t}$ exists, then the infimum in (39) is attained for the tuple (t,s)=\bigl{(}t,s(t)\bigr{)} with

[TABLE]

These pairs fulfill $\varphi(t)+\psi(s)=c(t,s)$. The choice is unique except for $t=0$, where any $s$ larger than $\frac{\varphi^{\prime\prime}(0+)}{2}$ is also a solution. 3. iii)

If $\varphi^{\prime}(t)>0$ for $t>0$ and $\varphi^{\prime\prime}(0+)>0$, then $s(t)\in(0,\frac{\varphi^{\prime\prime}(0+)}{2}]$ for all $t>0$.

Functions $\varphi$ fulfilling the assumptions of Proposition 4.2 are listed in Table 1.

Then, replacing $\varphi$ in $\operatorname{TV}_{\varphi}$ in the isotropic setting in (29) by the expression in (39), we can minimize instead of ${\mathcal{J}}_{\varphi}(u)$ the functional

[TABLE]

where \mathrm{d}_{i}\coloneqq\bigl{(}\sum_{j\in\operatorname{\mathcal{N}}(i)^{+}}\operatorname{dist}^{2}(u_{i},u_{j})\bigr{)}^{\frac{1}{2}}. We apply alternating minimization over $v=(v_{i})_{i\in\operatorname{{\mathcal{G}}}}\in\mathbb{R}^{n}$ and $u\in\mathcal{M}^{n}$ and obtain together with (40) the following iterations:

[TABLE]

The minimization over $u$ means to find a minimizer of

[TABLE]

Here, we can apply, e.g., a gradient descent or a Riemann-Newton method, see [1]. Concerning the convergence of the algorithm we have the following theorem, see [10].

Theorem 4.3** **(Convergence of half-quadratic minimization)

Let $\mathcal{H}$ be an Hadamard manifold and let $\varphi\colon\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ fulfill the assumptions of Proposition 4.2. Then, the sequence $\{u^{(r)}\}_{r\in\mathbb{N}}$ generated by (44) converges to the uniquely determined minimizer of ${\mathcal{J}}_{\varphi}$.

4.3 Proximal Point and Douglas-Rachford Algorithm

In the Euclidean setting, tools from convex analysis, in particular powerful algorithms based on duality theory, were successfully applied to minimize the proposed functionals. A prominent example is the alternating directions method of multipliers (ADMM), which is equivalent to the Douglas-Rachford algorithm. A central ingredient of these algorithms are proximal mappings which can be efficiently computed for special regularization terms appearing in Euclidean image processing tasks, see [19, 22]. Recently, several attempts have been made to translate these concepts to manifolds and it turns out that on Hadamard manifolds a certain theory of convex functions can be established. For example, the (inexact) cyclic proximal point algorithm can be introduced on these manifolds [4], and this method was also used to minimize the functional with first and second order $\operatorname{TV}$ regularizers in [5, 15, 20, 63]. Since the classical Douglas-Rachford algorithm relies on point reflections, it was natural to extend this algorithm to symmetric Hadamard manifolds [16].

4.3.1 Proximal Mapping

For $\lambda>0$ and a proper, convex, lsc function $\varphi\colon\mathbb{R}^{m}\rightarrow(-\infty,+\infty]$, the proximal mapping $\operatorname{prox}_{\lambda\varphi}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ defined by

[TABLE]

is uniquely determined. The counterpart on manifolds reads for $\varphi:\mathcal{M}^{m}\rightarrow[-\infty,+\infty]$ as

[TABLE]

Indeed, for proper, convex, lower semi-continuous functions on Hadamard manifolds $\mathcal{M}=\mathcal{H}$, the above minimizer is uniquely determined [38]. Moreover, the proximal operator is nonexpansive. For the distance functions appearing in the sums of our models $\varphi_{0}^{p}(x)\coloneqq\operatorname{dist}(x,y)^{p}$, $p\in\{1,2\}$ and $\varphi_{1}^{p}(x,y)\coloneqq\operatorname{dist}(x,y)^{p}$, $p\in\{1,2\}$, the proximal mapping can be given analytically, see [32, 63]. For our absolute second order differences \varphi_{2}(x,y,z)\coloneqq{\rm d}_{2}\big{(}x,y,z\big{)} the proximal mapping can be computed numerically on certain manifolds by the (sub)gradient descent algorithm and Lemma 2.1 as outlined in [5]. An analytical expression for the proximal mapping of $\varphi_{2}$ on the sphere $\mathbb{S}^{1}$ was given in [15].

The results are summarized in the following lemmas.

Proposition 4.4** **(Proximal mapping of distance functions)

Let $\mathcal{H}$ be an Hadamard manifold, $\lambda>0$ and $y\in\mathcal{H}$.

• i)

The proximal mappings of $\varphi_{0}^{p}(x)\coloneqq\frac{1}{p}\operatorname{dist}(\cdot,y)^{p}$, $p\in\{1,2\}$, are given by

[TABLE]

• ii)

The proximal mappings of $\varphi_{1}^{p}(x,y)\coloneqq\operatorname{dist}(x,y)^{p}$, $p\in\{1,2\}$, are given by

[TABLE]

To give the analytical expressions for the proximal mappings of $\varphi_{\nu}$, $\nu\in\{0,1,2\}$, on $\mathbb{S}^{1}$, we represent its elements by the angles in $[-\pi,\pi)$. For $a\in\mathbb{R}$, we denote by $(a)_{2\pi}\in[-\pi,\pi)$ those number for which there exists $k\in\mathbb{Z}$ such that $a+2\pi k=(a)_{2\pi}$.

Proposition 4.5** **(Proximal mapping of distance functions on $\mathbb{S}^{1}$)

Let $w_{1}\coloneqq(-1,1)^{\mbox{\tiny{T}}}$, $w_{2}\coloneqq(1,-2,1)^{\mbox{\tiny{T}}}$ and $s_{\nu}\coloneqq\operatorname{sgn}(\langle x,w_{\nu}\rangle)_{2\pi}$, $\nu\in\{1,2\}$. Then, for $\nu\in\{1,2\}$, the following holds true:

• i)

If $\lvert(\langle x,w_{\nu}\rangle)_{2\pi}\rvert<\pi$, then

[TABLE]

• ii)

If $\lvert(\langle x,w_{\nu}\rangle)_{2\pi}\rvert=\pi$, then the proximal mapping is two-fold

[TABLE]

• iii)

If $\lvert(\langle x,w_{\nu}\rangle)_{2\pi}\rvert<\pi$, then

[TABLE]

• iv)

If $\lvert(\langle x,w_{\nu}\rangle)_{2\pi}\rvert=\pi$, then the proximal mapping is two-fold

[TABLE]

• v)

Finally,

[TABLE]

where

[TABLE]

4.3.2 Cyclic Proximal Point Algorithm

Our functionals have the general form

[TABLE]

where appropriate splittings into the $K$ summands must be determined. Given a starting point $x^{(0)}\in\mathcal{M}$ the cyclic proximal point algorithm (CPPA) iterates

[TABLE]

We have the following convergence result from [5].

Theorem 4.6** **(Convergence of cyclic PPA)

Let ${\mathcal{H}}$ be a Hadamard manifold and $\varphi_{k}\colon{\mathcal{H}}\rightarrow\mathbb{R}$, $k=1,\ldots,K$, convex continuous functions such that ${\mathcal{J}}$ attains a (global) minimum. Assume that there exist $p\in{\mathcal{H}}$ and $C>0$ such that for each $k=1,\dots,K$ and all $x,y\in{\mathcal{H}}$ we have

[TABLE]

Then the sequence $\{x^{(r)}\}_{r\in\mathbb{N}}$ with non-negative $\{\tau_{r}\}_{r\in\mathbb{N}}\in\ell^{2}\backslash\ell_{1}$ converges for every starting point $x^{(0)}$ to a minimizer of ${\mathcal{J}}$.

The result can be generalized for the inexact cyclic PPA which iteratively generates the points $x^{(r+\tfrac{k}{K})}$, $k=1,\ldots,K$, $r\in\mathbb{N}_{0}$, fulfilling

[TABLE]

where $\{\varepsilon_{r}\}_{r\in\mathbb{N}_{0}}$ is a given sequence of positive reals with $\sum_{r=1}^{\infty}\varepsilon_{r}<\infty$, see [5].

4.3.3 Douglas-Rachford Algorithm for Symmetric Hadamard Spaces

The Douglas-Rachford (DR) algorithm relies on reflections. In the Euclidean setting, the reflection of a proper, convex, lsc function $\varphi\colon\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ is defined as

[TABLE]

It is a nonexpansive operator on $\mathbb{R}^{d}$ with respect to the Euclidean norm. Given two proper, convex, lsc functions $\varphi,\psi\colon\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$, the DR algorithm aims to solve

[TABLE]

by iterating a starting point $t^{(0)}$ as follows:

[TABLE]

Note that $x^{(r)}$ must be only computed in the final step of the algorithm. It is known that the DR algorithm converges if $\mathrm{ri}(\operatorname{dom}\varphi)\cap\mathrm{ri}(\operatorname{dom}\psi)\not=\emptyset$, a minimizer exists, $\eta>0$ and $\sum_{r\in\mathbb{N}}\tau_{r}(1-\tau_{r})=+\infty$. The DR algorithm can be considered a special case of the Krasnoselski–Mann iteration

[TABLE]

with $T\coloneqq R_{\eta\varphi}R_{\eta\psi}$. It is well-known that the sequence of iterates to a fixed point of $T$ if $T$ is nonexpansive. This is clearly the case for our setting since reflections $R_{\varphi}$ are nonexpansive and the concatenation of nonexpansive functions is nonexpansive again.

For two points $x,a\in\mathcal{H}$ on a symmetric Hadamard manifold, the geodesic reflection (8) can be written as $R_{p}(x)=\exp_{p}(-\log_{p}x)$. Then, the geodesic reflection of a proper, convex, lsc function $\varphi\colon\mathcal{H}^{n}\rightarrow(-\infty,+\infty]$ is the mapping

[TABLE]

Now in order to minimize

[TABLE]

the DR algorithm can be generalized as follows:

[TABLE]

Again, the algorithm can be seen as special case of the Krasnoselski–Mann iteration

[TABLE]

It was proved in [39], see also [3, Theorem 6.2.1] that such iteration converges to a fixed point of $T$, if $T$ is nonexpansive, has a nonempty fixed point set and $\sum_{r\in\mathbb{N}}\tau_{r}(1-\tau_{r})=+\infty$. Unfortunately, reflections at proper, convex, lsc functions on Hadamard manifolds are in general not nonexpansive. However, for the distance functions involved in our functionals nonexpansivness is guaranteed by the following theorem.

Theorem 4.7** **(Reflections at distance functions)

For an arbitrary fixed $a\in{\mathcal{H}}$ and $\varphi(x)\coloneqq\operatorname{dist}^{p}(a,x)$, $p\in\{1,2\}$, the geodesic reflection ${R}_{\eta\varphi}$, $\eta>0$ is nonexpansive. For $\varphi(x,y)\coloneqq\operatorname{dist}^{p}(x,y)$, $p\in\{1,2\}$, the geodesic reflection $R_{\eta\varphi}$, $\eta>0$, is nonexpansive.

In general we have more than two summands, i.e., we are interested in the minimization of (53) where $\varphi_{k}\colon\mathcal{H}^{n}\rightarrow(-\infty,+\infty]$, $k=1,\ldots,K,$ are proper, convex, lsc functions. Here the trick is to rewrite the functional as the sum of two special components

[TABLE]

where $\Phi(x)\coloneqq\sum_{k=1}^{K}\varphi_{k}(x_{k})$, $x\coloneqq(x_{k})_{k=1}^{K},$ and

[TABLE]

Obviously, D is a nonempty, closed convex set so that its indicator function is proper, convex and lsc, see [3, p. 37]. Now, the DR algorithm can be formulated as

[TABLE]

Note that the second step is indeed only necessary in the final iteration. Concerning the last step note that

[TABLE]

The minimizer of the sum is the so-called Karcher mean, which can be efficiently computed on Hadamard manifolds using the gradient descent algorithm or the cyclic proximal point algorithm.

Concerning the convergence of the parallel DR algorithm, for our setting, ${R}_{\eta\Phi}$ is nonexpansive since it contains only geodesic reflections of the distance functions in Theorem 4.7. Unfortunately, in symmetric Hadamard manifolds, geodesic reflections corresponding to orthogonal projections onto convex sets are in general not nonexpansive. This is also true for our special set D. The situation changes if we consider manifolds with constant curvature $\kappa$. Here those reflections are nonexpansive, see [30, 51]. So, in summary, although the parallel DR algorithm showed very good numerical performance on general Hadamard manifolds in [16], theoretical convergence results remain up to now limited to manifolds with constant non-positive curvature.

5 Numerical Examples

In this section, we give some illustrative numerical examples. The experiments are carried out using Matlab 2017a and the MVIRT toolbox [9]222Open source, available at ronnybergmann.net/mvirt/. As a quality measure we use the mean squared error (MSE) defined by

[TABLE]

where $u_{0}$ denotes the original image. The parameters in the models were obtained via a grid search with respect to the optimal $\epsilon$ and can be found in detail in the respective papers [10, 13]. In Figures 3 and 4 we compare the performance of our variational models with different regularizers, where TV-TV2 denotes the additive coupling of $\operatorname{TV}$ and $\operatorname{TV}_{2}$. For comparison we added the results obtained with the patched-based methods from [14], namely with nonlocal means (NL-means) and nonlocal MMSE (NL-MMSE). In brackets we give the corresponding error value $\epsilon$.

Figures 5 and 6 show denoising results obtain by the half quadratic minimization. Figures 5 contains results for the different functions $\varphi$ in Table 1 and for the nonsmooth $\operatorname{TV}$ regularizer. Here the parameters are optimized with respect to the PSNR of the RGB images and the numbers in brackets refer to the PSNR. In Figure 6 we present a denoising result for the 3D DT-MRI image from Figure 1.

8 RECOMMENDED READING

• [1]

P.-A. Absil, R. Mahony, and R. Sepulchre.

Optimization Algorithms on Matrix Manifolds.

Princeton University Press, Princeton, 2008.

• [2]

F. Alouges.

A new algorithm for computing liquid crystal stable configurations: the harmonic mapping case.

SIAM Journal on Numerical Analysis, 34(5):1708–1726, 1997.

• [3]

M. Bačák.

Computing medians and means in Hadamard spaces.

SIAM Journal on Optimization, 24(3):1542–1566, 2014.

• [4]

M. Bačák.

Convex Analysis and Optimization in Hadamard Spaces, volume 22 of De Gruyter Series in Nonlinear Analysis and Applications.

De Gruyter, Berlin, 2014.

• [5]

M. Bačák, R. Bergmann, G. Steidl, and A. Weinmann.

A second order non-smooth variational model for restoring manifold-valued images.

SIAM Journal on Scientific Computing, 38(1):A567–A597, 2016.

• [6]

F. Bachmann and R. Hielscher.

MTEX – MATLAB toolbox for quantitative texture analysis.

http://mtex-toolbox.github.io/, 2005–2016.

• [7]

F. Bachmann, R. Hielscher, and H. Schaeben.

Grain detection from 2d and 3d EBSD data—specification of the MTEX algorithm.

Ultramicroscopy, 111(12):1720–1733, 2011.

• [8]

T. Batard and M. Bertalmío.

A geometric model of brightness perception anf its application to color image correction.

Journal of Mathematical Imaging and Vision, 60(6):849–881, 2018.

• [9]

R. Bergmann.

MVIRT, a toolbox for manifold-valued image restoration.

In IEEE International Conference on Image Processing, IEEE ICIP 2017, Beijing, China, September 17–20, 2017, 2017.

• [10]

R. Bergmann, R. H. Chan, R. Hielscher, J. Persch, and G. Steidl.

Restoration of manifold-valued images by half-quadratic minimization.

Inverse Problems and Imaging, 10(2):281–304, 2016.

• [11]

R. Bergmann, J. H. Fitschen, J. Persch, and G. Steidl.

Infimal convolution coupling of first and second order differences on manifold- valued images.

In F. Lauze, Y. Dong, and A. B. Dahl, editors, Scale Space and Variational Methods in Computer Vision 2017, page 447–459. Springer, Cham, 2017.

• [12]

R. Bergmann, J. H. Fitschen, J. Persch, and G. Steidl.

Iterative multiplicative filters for data labeling.

International Journal of Computer Vision, 123(3):435–453, 2017.

• [13]

R. Bergmann, J. H. Fitschen, J. Persch, and G. Steidl.

Priors with coupled first and second order differences for manifold-valued image processing.

Journal of Mathematical Imaging and Vision, 60(9):1459–1481, 2018.

• [14]

R. Bergmann, F. Laus, J. Persch, and G. Steidl.

Manifold-valued image processing.

SIAM News, 50(8):1,3, 2017.

• [15]

R. Bergmann, F. Laus, G. Steidl, and A. Weinmann.

Second order differences of cyclic data and applications in variational denoising.

SIAM Journal on Imaging Sciences, 7(4):2916–2953, 2014.

• [16]

R. Bergmann, J. Persch, and G. Steidl.

A parallel Douglas–Rachford algorithm for restoring images with values in symmetric Hadamard manifolds.

SIAM Journal on Imaging Sciences, 9(3):901–937, 2016.

• [17]

R. Bergmann and A. Weinmann.

Inpainting of cyclic data using first and second order differences.

In Energy Minimization Methods in Computer Vision and Pattern Recognition, pages 155–168. Springer, 2015.

• [18]

N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre.

Manopt, a matlab toolbox for optimization.

Journal of Machine Learning Research, 15:1455–1459, 2014.

• [19]

S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.

Distributed optimization and statistical learning via the alternating direction method of multipliers.

Foundations and Trends in Machine Learning, 3(1):101–122, 2011.

• [20]

K. Bredies, M. Holler, M. Storath, and A. Weinmann.

Total generalized variation for manifold-valued data.

SIAM Journal on Imaging Sciences, 11(3):1785–1848, 2018.

• [21]

K. Bredies, K. Kunisch, and T. Pock.

Total generalized variation.

SIAM Journal on Imaging Sciences, 3(3):492–526, 2010.

• [22]

M. Burger, A. Sawatzky, and G. Steidl.

First order algorithms in variational image processing.

In R. Glowinski, S. Osher, and W. Yin, editors, Operator Splittings and Alternating Direction Methods. Springer, 2017.

• [23]

R. Bürgmann, P. A. Rosen, and E. J. Fielding.

Synthetic aperture radar interferometry to measure earth’s surface topography and its deformation.

Annual Review of Earth and Planetary Sciences, 28(1):169–209, 2000.

• [24]

A. Chambolle and P.-L. Lions.

Image recovery via total variation minimization and related problems.

Numerische Mathematik, 76(2):167–188, 1997.

• [25]

R. Ciak, M. Hirzmann, and O. Scherzer.

Regularization with metric double integrals of functions with values in a set of high-dimensional vectors.

Journal of Mathematical Imaging and Vision, accepted, 2019.

• [26]

P. A. Cook, Y. Bai, S. Nedjati-Gilani, K. K. Seunarine, M. G. Hall, G. J. Parker, and D. C. Alexander.

Camino: Open-source diffusion-MRI reconstruction and processing.

In 14th Scientific Meeting of the International Society for Magnetic Resonance in Medicine, page 2759, Seattle, WA, USA, 2006.

• [27]

D. Cremers and E. Strekalovskiy.

Total cyclic variation and generalizations.

Journal of Mathematical Imaging and Vision, 47(3):258–277, 2013.

• [28]

M. P. do Carmo.

Riemannian Geometry, volume 115.

Birkhäuser, Basel, 1992.

Tranlated by F. Flatherty.

• [29]

J. Ehlers, F. A. E. Pirani, and A. Schild.

The geometry of free fall and light propagation.

In L. O’Reifeartaigh, editor, General Relativitiy, pages 63–84. Oxford University Press, 1972.

• [30]

A. Fernández-León and A. Nicolae.

Averaged alternating reflections in geodesic spaces.

Journal of Mathematical Analysis and Applications, 402(2):558 – 566, 2013.

• [31]

O. P. Ferreira and P. R. Oliveira.

Subgradient algorithm on Riemannian manifolds.

Journal of Optimization Theory and Applications, 97(1):93–104, Apr 1998.

• [32]

O. P. Ferreira and P. R. Oliveira.

Proximal point algorithm on Riemannian manifolds.

Optimization, 51(2):257–270, 2002.

• [33]

D. Geman and G. Reynolds.

Constrained restoration and the recovery of discontinuities.

IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(3):367–383, 1992.

• [34]

D. Geman and C. Yang.

Nonlinear image recovery with half-quadratic regularization.

IEEE Transactions on Image Processing, 4(7):932–946, 1995.

• [35]

M. Giaquinta, G. Modica, and J. Souček.

Variational problems for maps of bounded variation with values in $S^{1}$.

Calculus of Variation, 1(1):87–121, 1993.

• [36]

M. Giaquinta and D. Mucci.

Maps of bounded variation with values into a manifold: total variation and relaxed energy.

Pure and Applied Mathematics Quarterly, 3(2):513–538, 2007.

• [37]

P. Grohs and M. Sprecher.

Total variation regularization on Riemannian manifolds by iteratively reweighted minimization.

Information and Inference: A Journal of the IMA, 5(4):353–378, 2016.

• [38]

J. Jost.

Nonpositive Curvature: Geometric and Analytic Aspects.

Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, 1997.

• [39]

B. Kakavandi.

Weak topologies in complete CAT(0) metric spaces.

Proceedings of the American Mathematical Society, 141(3):1029–1039, 2013.

• [40]

A. Kheyfets, W. A. Miller, and G. A. Newton.

Schild’s ladder parallel transport procedure for an arbitrary connection.

International Journal of Theoretical Physics, 39(12):2891–2898, Dec 2000.

• [41]

R. Kimmel and N. Sochen.

Orientation diffusion or how to comb a porcupine.

Journal of Visual Communication and Image Representation, 13(1):238–248, 2002.

• [42]

F. Laus, M. Nikolova, J. Persch, and G. Steidl.

A nonlocal denoising algorithm for manifold-valued images using second order statistics.

SIAM Journal on Imaging Sciences, 10(1):416–448, 2017.

• [43]

M. Lebrun, A. Buades, and J.-M. Morel.

A nonlocal Bayesian image denoising algorithm.

SIAM Journal on Imaging Sciences, 6(3):1665–1688, 2013.

• [44]

J. Lellmann, E. Strekalovskiy, S. Koetter, and D. Cremers.

Total variation regularization for functions with values in a manifold.

In IEEE International Conference on Computer Vision, pages 2944–2951, 2013.

• [45]

C. Li, B. S. Mordukhovich, J. Wang, and J.-C. Yao.

Weak sharp minima on Riemannian manifolds.

SIAM Journal on Optimization, 21(4):1523–1560, 2011.

• [46]

M. Lorenzi and X. Pennec.

Efficient parallel transport of deformations in time series of images: From Schild’s to pole ladder.

Journal of Mathematical Imaging and Vision, 50(1):5–17, Sep 2014.

• [47]

E. Muñoz-Moreno, R. Cárdenes-Almeida, and M. Martín-Fernández.

Review of techniques for registration of diffusion tensor imaging.

• [48]

S. Neumayer, J. Persch, and G. Steidl.

Morphing of manifold-valued images inspired by discrete geodesics in image spaces.

SIAM Journal on Imaging Sciences, 11(3):1898–1930, 2018.

• [49]

M. Nikolova and M. K. Ng.

Analysis of half-quadratic minimization methods for signal and image recovery.

SIAM Journal on Scientific Computing, 27(3):937–966, 2005.

• [50]

X. Pennec.

Parallel transport with pole ladder: a third order scheme in affine connection spaces which is exact in affine symmetric spaces.

ArXiv Preprint, 1805.1143, 2018.

• [51]

J. Persch.

Recent Advances in Denoising of Manifold-Valued Images.

TU Kaiserslautern, 2018.

• [52]

W. Ring and B. Wirth.

Optimization methods on Riemannian manifolds and their application to shape space.

SIAM Journal on Optimization, 22(2):596–627, 2012.

• [53]

F. Rocca, C. Prati, and A. M. Guarnieri.

Possibilities and limits of SAR interferometry.

ESA SP, pages 15–26, 1997.

• [54]

G. Rosman, M. Bronstein, A. Bronstein, A. Wolf, and R. Kimmel.

Group-valued regularization framework for motion segmentation of dynamic non-rigid shapes.

In Scale Space and Variational Methods in Computer Vision, pages 725–736. Springer, 2012.

• [55]

G. Rosman, X.-C. Tai, R. Kimmel, and A. M. Bruckstein.

Augmented-Lagrangian regularization of matrix-valued maps.

Methods and Applications of Analysis, 21(1):121–138, 2014.

• [56]

L. I. Rudin, S. Osher, and E. Fatemi.

Nonlinear total variation based noise removal algorithms.

Physica D, 60(1):259–268, 1992.

• [57]

O. Tuzel, F. Porikli, and P. Meer.

Region covariance: A fast descriptor for detection and classification.

In European Conference on Computer Vision, pages 589–600. Springer, 2006.

• [58]

O. Tuzel, F. Porikli, and P. Meer.

Learning on Lie groups for invariant detection and tracking.

In CVPR 2008, pages 1–8. IEEE, 2008.

• [59]

C. Udrişte.

Convex Functions and Optimization Methods on Riemannian Manifolds, volume 297 of Mathematics and its Applications.

Kluwer Academic Publishers Group, Dordrecht, 1994.

• [60]

T. Valkonen, K. Bredies, and F. Knoll.

Total generalized variation in diffusion tensor imaging.

SIAM Journal on Imaging Sciences, 6(1):487–525, 2013.

• [61]

X. M. Wang.

Subgradient algorithms on Riemannian manifolds of lower bounded curvatures.

Optimization, 67(1):179–194, 2018.

• [62]

X. M. Wang, C. Li, and J. C. Yao.

Subgradient projection algorithms for convex feasibility on Riemannian manifolds with lower bounded curvatures.

Journal of Optimization Theory and Applications, 164(1):202–217, Jan 2015.

• [63]

A. Weinmann, L. Demaret, and M. Storath.

Total variation regularization for manifold-valued data.

SIAM Journal on Imaging Sciences, 7(4):2226–2257, 2014.

• [64]

O. Yair, M. Ben-Chen, and R. Talmon.

Parallel transport on the cone manifold of spd matrices for domain adaptation.

ArXiv Preprint, 1807.10479, 2018.

• [65]

P. Zhang, M. Niethammer, D. Shen, and P.-T. Yapa.

Large deformation diffeomorphic registration of diffusion-weighted imaging data.

Medical Image Analysis, 18(8):1290–1298, 2018.