A Graph Framework for Manifold-valued Data

TL;DR

This paper extends graph-based variational and PDE methods to manifold-valued data, enabling advanced processing of complex data types like DTI and LiDAR through novel operators and numerical algorithms.

Contribution

It introduces a generalized graph framework for manifold-valued data, including calculus and operators for variational models and PDEs, expanding beyond Euclidean domains.

Findings

Successfully applied to synthetic and real-world data

Enabled total variation and Tikhonov regularization on manifolds

Demonstrated effectiveness on DTI and LiDAR data

Abstract

Graph-based methods have been proposed as a unified framework for discrete calculus of local and nonlocal image processing methods in the recent years. In order to translate variational models and partial differential equations to a graph, certain operators have been investigated and successfully applied to real-world applications involving graph models. So far the graph framework has been limited to real- and vector-valued functions on Euclidean domains. In this paper we generalize this model to the case of manifold-valued data. We introduce the basic calculus needed to formulate variational models and partial differential equations for manifold-valued functions and discuss the proposed graph framework for two particular families of operators, namely, the isotropic and anisotropic graph~$p$-Laplacian operators, $p\geq1$. Based on the choice of $p$ we are in particular able to solve optimization problems on manifold-valued functions involving total variation ($p=1$) and Tikhonov ($p=2$) regularization. Finally, we present numerical results from processing both synthetic as well as real-world manifold-valued data, e.g., from diffusion tensor imaging (DTI) and light detection and ranging (LiDAR) data.

Equations

$${\mathcal{D}}_x:=\{y\in \mathcal{M}:y={\exp }_x\xi ,\text{ for some }\xi \in T_x\mathcal{M}\text{ with }\Vert \xi {\Vert }_x<r_x\}.$$

$$D_{V(t)}W(t)=0$$

$$f({\gamma }_{\stackrel{⌢}{x,y}}(t)),t\in [0,d_{\mathcal{M}}(x,y)],$$

$$\begin{array}{rl}f:V& \to \mathcal{M},\\ u& \mapsto f(u).\end{array}$$

$${\mathrm{T}}_f\mathcal{M}:=\dot{\bigcup _{u\in V}}{\mathrm{T}}_{f(u)}\mathcal{M}\subset \mathrm{T}\mathcal{M}.$$

$$\mathcal{H}(V;\mathcal{M})\ \coloneqq\ \{f\colon V\rightarrow\mathcal{M}\bigl{|}f(u)\in\mathcal{D}_{f(v)}\text{ for all }(u,v)\in E\}\ ,$$

$$d_{\mathcal{H}(V;\mathcal{M})}(f,g)\ \coloneqq\ \biggl{(}\sum_{u\in V}d_{\mathcal{M}}^{2}(f(u),g(u))\biggr{)}^{\frac{1}{2}}.$$

$$\begin{array}{rl}{H}_f:E& \to {\mathrm{T}}_f\mathcal{M},\\ (u,v)& \mapsto H_f(u,v)\in {\mathrm{T}}_{f(u)}\mathcal{M}.\end{array}$$

$$\displaystyle\lVert H_{f}\rVert_{\ell_{p,q}(E;\mathrm{T}_{f}\mathcal{M})}\ \coloneqq\ \Biggl{(}\frac{2}{p}\sum_{u\in V}\biggl{(}\sum_{v\sim u}\lVert H_{f}(u,v)\rVert_{f(u)}^{q}\biggr{)}^{\frac{p}{q}}\biggr{)}^{\frac{1}{p}}\ .$$

$$⟨H_f,G_f{⟩}_{\mathcal{H}(E;{\mathrm{T}}_f\mathcal{M})}:=\sum _{(u,v)\in E}⟨H_f(u,v),G_f(u,v){⟩}_{f(u)}.$$

$$\mathcal{H}(E):=\bigcup _{f\in \mathcal{H}(V;\mathcal{M})}\mathcal{H}(E;{\mathrm{T}}_f\mathcal{M}).$$

$$\displaystyle\lVert H_{f}\rVert^{q}_{q,f(u)}\$$

$$\displaystyle\langle f,g\rangle\coloneqq\sum_{u\in V}\biggl{(}\sum_{v\sim u}\bigl{\langle}\log_{f(u)}f(v),\operatorname{PT}_{g(u)\to f(u)}\log_{g(u)}g(v)\bigr{\rangle}_{f(u)}\biggr{)}.$$

$$\lVert f\rVert_{p}^{p}\coloneqq\sum_{u\in V}\biggl{(}\sum_{v\sim u}\bigl{\langle}\log_{f(u)}f(v),\log_{f(u)}f(v)\bigr{\rangle}_{f(u)}\biggr{)}^{\frac{p}{2}},\qquad p\geq 1.$$

$$\mathcal{H}(V;T_f\mathcal{M}):=\{H:V\to T_f\mathcal{M},H(u)\in T_{f(u)}\mathcal{M}\},$$

$${\partial }_{v}f(u):=\sqrt{w(u,v)}{\log }_{f(u)}f(v).$$

$${\partial }_{v}f(u)$$

$$=-{\mathrm{PT}}_{f(v)\to f(u)}({\partial }_{u}f(v)).$$

$$\nabla f(u,v):={\partial }_{v}f(u)=\sqrt{w(u,v)}{\log }_{f(u)}f(v),(u,v)\in E.$$

$$\begin{array}{rl}⟨\nabla f,H_f& {⟩}_{\mathcal{H}(E;{\mathrm{T}}_f\mathcal{M})}=\\ & \sum _{u\in V}\sum _{v\sim u}⟨{\log }_{f(u)}f(v),\frac{\sqrt{w(u,v)}}{2}{H}_f(u,v)-\frac{\sqrt{w(v,u)}}{2}{\mathrm{PT}}_{f(v)\to f(u)}(H_f(v,u)){⟩}_{f(u)}.\end{array}$$

$${⟨\nabla f,H_f⟩}_{\mathcal{H}(E;{\mathrm{T}}_f\mathcal{M})}=\sum _{(u,v)\in E}{⟨\nabla f(u,v),H_f(u,v)⟩}_{f(u)}$$

$$\displaystyle=\ \frac{1}{2}\sum_{(u,v)\in E}\bigl{\langle}\sqrt{w(u,v)}\log_{f(u)}f(v),H_{f}(u,v)\bigr{\rangle}_{f(u)}$$

$$\displaystyle\hskip 56.9055pt+\frac{1}{2}\sum_{(u,v)\in E}\bigl{\langle}\sqrt{w(u,v)}\log_{f(u)}f(v),H_{f}(u,v)\bigr{\rangle}_{f(u)}$$

$$\displaystyle=\ \sum_{(u,v)\in E}\Bigl{\langle}\log_{f(u)}f(v),\frac{\sqrt{w(u,v)}}{2}H_{f}(u,v)\Bigr{\rangle}_{f(u)}$$

$$\displaystyle\hskip 56.9055pt+\sum_{(v,u)\in E}\Bigl{\langle}\operatorname{PT}_{f(v)\rightarrow f(u)}(\log_{f(v)}f(u)),\operatorname{PT}_{f(v)\rightarrow f(u)}\left(\frac{\sqrt{w(v,u)}}{2}H_{f}(v,u)\right)\Bigr{\rangle}_{f(u)}$$

$$\displaystyle=\ \sum_{(u,v)\in E}\Bigl{\langle}\log_{f(u)}f(v),\frac{\sqrt{w(u,v)}}{2}H_{f}(u,v)\Bigr{\rangle}_{f(u)}$$

$$\displaystyle\hskip 56.9055pt+\sum_{(v,u)\in E}\Bigl{\langle}-\log_{f(u)}f(v)),\frac{\sqrt{w(v,u)}}{2}\operatorname{PT}_{f(v)\rightarrow f(u)}\left(H_{f}(v,u)\right)\Bigr{\rangle}_{f(u)}$$

$$=\sum _{u\in V}\sum _{v\sim u}⟨{\log }_{f(u)}f(v),\frac{\sqrt{w(u,v)}}{2}{H}_f(u,v)-\frac{\sqrt{w(v,u)}}{2}{\mathrm{PT}}_{f(v)\to f(u)}(H_f(v,u)){⟩}_{f(u)}.$$

$$\mathrm{div}{H}_f(u):=\frac{1}{2}\sum _{v\sim u}\sqrt{w(v,u)}{\mathrm{PT}}_{f(v)\to f(u)}{H}_f(v,u)-\sqrt{w(u,v)}{H}_f(u,v).$$

$$\mathrm{div}{H}_f(u)=-\sum _{v\sim u}\sqrt{w(u,v)}{H}_f(u,v).$$

$${\Delta }_p^{\mathrm{a}}f(u):=$$

$$=$$

$${\Delta }_p^{\mathrm{i}}f(u):=$$

$$=$$

$$\Delta f(u):={\Delta }_2^{\mathrm{i}}f(u)=\mathrm{div}(\nabla f)(u)=-\sum _{v\sim u}w(u,v){\log }_{f(u)}f(v).$$

Full text

1 Introduction

Variational methods and partial differential equations (PDEs) play a key role for modeling and solving image processing tasks. Typically, in the setting of vector-valued functions the respective continuous formulations are discretized, e.g., by using finite differences or finite elements. These discretization schemes are so popular because they are well-investigated and in general preserve important properties of the continuous models, e.g., conservation laws or maximum principles. Recently, the trend is to establish numerical discretization schemes for functions not living on subsets of Euclidean spaces but on surfaces, which might be given as manifolds in a continuous setting or as finite point clouds in a discrete setting. The discretization of differential operators becomes especially challenging for the latter case as there are a-priori no explicit neighborhood relationships for the raw point cloud data. Currently, there is high interest in translating variational methods and PDEs to discrete data, which are modeled by finite weighted graphs and networks. In fact, any discrete data can be represented by a finite graph in which the vertices are associated to the data and its edges correspond to relationships within the data. This modeling approach has led to interesting applications in mathematical image processing and machine learning so far. One major advantage of using graph-based methods is that they unify local and nonlocal methods. This is due to the reason that one may model data relationships not only based on geometric proximity but also based on similarity of features. Furthermore, graphs can model arbitrary spatial neighborhoods, e.g., when the domain itself is a submanifold of the space the measurements are taken on. In order to translate and solve PDEs on graphs, different discrete vector calculi have been proposed in the literature, e.g., see [45] and references therein. One rather simple concept is to use discrete partial differences [32, 43]. This mimetic approach allows to solve PDEs on both regular as well as irregular data domains by replacing continuous partial differential operators, e.g., gradient or divergence, by a reasonable discrete analogue. By employing this approach many important tools and results from the continuous setting can be transferred to finite weighted graphs.

On the other hand, variational methods from image processing have recently been generalized from real and vectorial data in Euclidean spaces to the case of manifold-valued data, most prominently the total variation (TV) regularization [75, 30, 58, 79] including efficient algorithms [17, 15] and second order differences based models [16, 8]. These manifold-valued images occur for example in interferometric synthetic aperture radar (InSAR) imaging [62, 25, 31], where the measured phase-valued data may be noisy and/or incomplete. Sphere-valued data appears, e.g., in directional analysis [51, 61]. Data items with values in the special rotation group $\operatorname{SO}(3)$ (or a quotient manifold thereof) are used to process orientations (including certain symmetries) e.g., in electron backscattered diffraction imaging (EBSD) [52, 11, 2]. Another application is diffusion tensor imaging (DT-MRI) [14, 36, 66, 78], where the diffusion tensors can be represented as $3\times 3$ symmetric positive definite matrices, which also form a manifold. In general when dealing with covariance matrices the given data is also a set of symmetric positive definite matrices of size $n\times n$. Finally, multivariate Gaussian distributions can be characterized as values on a hyperbolic space, namely the Poincaré half-space [4]. Two examples are shown in Figure 1.

The data of Mt. Vesuvius shown in Fig. 1(a)) obtained in InSAR measurements is available from [67]111see https://earth.esa.int/workshops/ers97/program-details/speeches/rocca-et-al/. The data of a phase value in each pixel is given on the regular rectangular image grid. Each pixel represents the “height modulo wavelength” at the measured position, i.e. a phase. Note that the distance measure for phase values is different from the classical images with values on the real line. This is indicated by the hue color map. Values near both the minimal and maximal value are again close with respect to the distance on the sphere $\mathbb{S}^{1}$. In Figure 1(b)) data measured on the implicitly given surface of the human brain is shown. The data is available within the Camino toolbox [29]222see http://camino.cs.ucl.ac.uk. Within the extracted surface data set, each pixel is given on the implicit surface, see Section 5.3 for details on the extraction, while each pixel is a value representing the diffusion tensor at this point, i.e. a symmetric positive definite matrix. These are visualized using their eigenvectors and eigenvalues as main axes and main axes lengths, respectively. Since the data is given on an implicit surface, the vicinity of pixel or the modeling of pixel being neighbors, can be done using a graph. This graph has the pixel as their nodes and edges to neighboring pixel, e.g. whenever two pixel are below a certain threshold in the $\mathbb{R}^{3}$ the data is measured in.

All the mentioned applications including the two examples from Figure 1 have in common that they resemble manifold-valued images or more general manifold-valued data. This yields discrete energy functionals that map from a tensor product of a manifold with itself (the size of the number of pixels) into the (extended) real line. Starting from the early ’90s [77] several groups worked on optimization methods on manifolds, e.g. Riemannian gradient, subgradient and proximal point methods [35, 34, 7, 47], Newton iterations, e.g., see the book [1] for algorithms on matrix manifolds. In all these methods there are three main challenges that emerge for the case of Riemannian manifolds: (i) there is (in general) no additive group on a Riemannian manifold, (ii) straight lines are replaced by geodesics, i.e., the space itself is curved and this curvature has to be taken care of, and most importantly for variational methods, (iii) there is a notion of global convexity only for the case of Hadamard manifolds. These manifolds have nonpositive sectional curvature and many tools from optimization have been transferred to these spaces, see the monograph [9]. In this paper we combine for the first time differential operators on finite weighted graphs and optimization methods on Riemannian manifolds to construct local and nonlocal methods for manifold-valued data based on variational models and partial differential equations.

1.1 Related work

Due to the growing interest in graph-based modeling in the recent years there have been many sophisticated methods based on nonlocal models and finite weighted graphs. While one part of the community has concentrated on the theory of nonlocal methods [6, 24, 43, 57] and the translation of finite weighted graphs to continuous models [38, 39, 60], other groups investigated graph operators and their applications in the discrete setting. As we are proposing a discrete graph framework for manifold-valued data in this work, we will focus our discussion of related work on the latter research domain. Several groups have put their efforts in the translation of differential operators and related PDEs from the continuous setting to finite weighted graphs and the study of their discrete properties, such as the Ginzburg-Landau functional [40, 26], the Allen-Cahn equation [41], variational $p$-Laplacians [32, 50, 33] and its application to spectral clustering [23], and Hamilton-Jacobi equations [64, 76]. In addition to discrete data defined on regular grids, recent works have investigated the feasibility of modeling data defined on discrete surfaces and raw point clouds using finite weighted graphs [12, 38, 59, 72, 76]. Especially the latter scenario of point cloud processing is a difficult task, since there is a-priori no connectivity given for the acquired points. In particular, in a previous work with Elmoataz et al. in [33] we have studied variants of the graph $p$-Laplacian and the graph $\infty$-Laplacian and discussed their relationship to continuous operators such as the variational $p$-Laplacians, the game $p$-Laplacian or the nonlocal $p$-Laplacian, which occur as models in many domains, e.g., in physics, game theory, biology, or economy. With respect to image processing and machine learning applications the graph $p$-Laplacian and the graph $\infty$-Laplacian have been successfully used for denoising, segmentation, and inpainting of images but also for general data processing and clustering [32, 33, 59, 72].

All the above listed graph-based methods share in common that they can only be applied to real- or vector-valued data. To the best of our knowledge, there is no generalization of finite weighted graphs to manifold-valued data so far. However, for manifold-valued images and 3D data sets, the same tasks arise as in usual image processing, e.g., denoising, inpainting or segmentation. Recently several works tackled these tasks such as [18, 30, 75] for inpainting, or [5, 20, 13] for segmentation of such data. For denoising the TV approach or Rudin-Osher-Fatemi (ROF) [70] model was introduced by [79, 58] and generalized to second order methods in [16, 8, 22, 18]. Furthermore, for the ROF model half-quadratic minimization [15] and the Douglas-Rachford algorithm [17] have led to a significant increase in computational efficiency. Another approach in [53] uses second order statistics and employs a nonlocal denoising method. All the discussed methods share in common that they work on discrete regular grids only, e.g., on pixel or voxel grids.

It is possible to embed every $d$-dimensional manifold into an Euclidean space of at most dimension $2d$ due to the theorem of Whitney [80], or even use an isometric embedding following the theorem of Nash [65]. However, a major disadvantage is that this might increase the dimension of the data, by a factor up to two following Whitney or even more when relying on an isometric embedding. While embedding is sometimes easy, e.g., for the sphere $\mathbb{S}^{d}\subset\mathbb{R}^{d+1}$, it is often beneficial to use the intrinsic metric from an application point of view. For example for measurements on the earth the arc length is needed instead of the Euclidean norm in the embedding space $\mathbb{R}^{3}$. There are two further problems when working in the embedding space. First, minimizing within the embedding space can be done by adding an enforcing constraint to a variational model, e.g., the indicator function $\iota_{\mathcal{M}}$ being [math] within the manifold and $+\infty$ otherwise. These can easily be introduced into real-valued state-of-the-art methods like alternating direction method of multipliers (ADMM) [37, 44]. For example in  [68, 69] a matrix-valued TV functional was introduced using the ADMM for the classical TV functional extended by a such projection. However, when considering the example of symmetric positive definite matrices, the set is open. Hence the constraint can not be enforced by introducing a projection. The authors in the aforementioned papers project onto the $\varepsilon$-relaxed set or the closure and observe numerical convergence. However, from the mathematical viewpoint there remains the open question how to model such embeddings. Furthermore, several methods like infimal convolution [28], which was recently generalized in [19, 21] to Riemannian manifolds, split a signal or data set into several parts. When working in the embedding space directly these parts also have to be regarded in the embedding space and lose their manifold-valued interpretation. For these reasons working intrinsically on manifolds is often beneficial.

1.2 Own contribution

In this paper we introduce a novel graph framework for manifold-valued data. On the one hand our approach generalizes manifold-valued image processing models to arbitrary neighborhoods and discretizations which are modeled by the topology of the graph. Furthermore, this work also introduces a unified framework for both local and nonlocal methods for manifold-valued data processing. For local methods our framework introduces the possibility to process data not only on a pixel grid, but also for the case that the measurements are taken on surfaces. This surface might be explicitly given, e.g., measurements on the surface of a sphere or just implicitly by a point cloud. For the nonlocal case this framework unifies all methods for which vicinity is defined via the similarity of features, e.g., adaptive filtering methods or patch-based distances. On the other hand the data items are manifold-valued and hence a huge variety of data measurement modalities are incorporated in this framework. The proposed framework is very flexible as it consists of three independent parts, namely manifold specific operations, graph construction and operators, and numerical solvers. Each of these parts can be exchanged and enhanced without effecting the other modules. By introducing a comprehensive mathematical framework we also derive a notion of an anisotropic as well as an isotropic manifold-valued graph $p$-Laplacian. For the special case of the manifold being an Euclidean space, the operators in our framework simplify to the real- and vector-valued graph $p$-Laplacians. Thus, our approach can be interpreted as a generalization of well-known discrete graph operators to the manifold-valued setting. From the view of variational modeling in convex optimization we further derive optimality conditions for a family of energy functionals corresponding to denoising of manifold-valued data. Finally, by looking at the connection of these optimality conditions to parabolic PDEs, we derive a simple algorithm to compute minimizers of the respective energy functionals. We demonstrate the flexibility and performance of our approach on a variety of synthetic as well as real world applications and solve all these with the same universal numerical algorithm.

1.3 Organization

In Section 2 we first introduce the basic notation of finite weighted graphs and the needed tools from differential geometry to perform data processing on Riemannian manifolds. Then we introduce our graph framework for manifold-valued functions in Section 3 and define discrete differential operators to perform a huge variety of processing tasks on manifold-valued data. Subsequently, we discuss a special class of problems in Section 4, namely variational denoising tasks. We formulate the respective energy functionals, derive necessary optimimality conditions, and present two simple numerical schemes to compute solutions of these problems. In Section 5 we demonstrate the performance of our approach on several synthetic as well as real-world applications. Finally, we conclude our paper with a discussion and a short outlook to future work in Section 6.

2 Mathematical foundation

In this section we introduce the mathematical concepts and notations needed to define graph operators for manifold-valued data. We begin by giving the basic definitions of finite weighted graphs in Section 2.1. Subsequently, in Section 2.2 we introduce the necessary notation from differential geometry to describe (complete) Riemannian manifolds, tangent spaces, and functions defined on a manifold.

2.1 Finite weighted graphs

Finite weighted graphs are present in many different fields of research, e.g., image processing [32, 26, 81], machine learning [12, 39, 81], or network analysis [55, 64, 73] as they allow to model and process discrete data of arbitrary topology. Furthermore, one is able to translate variational methods and partial differential equations to finite weighted graphs and apply these to both local and nonlocal problems in the same unified framework. Although their exact description is application dependent, there exists a widely used consent of basic concepts and definitions for finite weighted graphs in the literature [32, 41, 43]. In the following we recall these basic concepts and the respective mathematical notation, which we will need to introduce a general graph framework for manifold-valued functions in Section 3 below.

Definition 2.1** (Finite weighted graph).**

A finite weighted (directed) graph $G$ is defined as a triple $G=(V,E,w)$ for which

• •

$V=\{1,\ldots,n\},n\in\mathbb{N}$, is a finite set of indices denoting the vertices,

• •

$E\subset V\times V$* is a finite set of (directed) edges connecting a subset of vertices,*

• •

$w\colon E\rightarrow\mathbb{R}^{+}$* is a nonnegative weight function defined on the edges of the graph.*

Typically, for given application data each graph vertex $u\in V$ models an entity in the data structure, e.g., elements of a finite set, pixels in an image, or nodes in a network. Due to the abstract nature of the graph structure it is also possible to identify sets of entities by a vertex, e.g., when summarizing vertices in hierarchical graphs [48]. It is important to distinguish between abstract data entities modeled by graph vertices and their respective characteristics, which will be modeled by vertex functions defined below. A graph edge $(u,v)\in E$ between a start node $u\in V$ and an end node $v\in V$ models certain relationships between two entities, e.g., geometric neighborhoods, interactions, or similarity relationships, depending on the given application. In our case, we consider directed edges, i.e., $(u,v)\in E\nLeftrightarrow(v,u)\in E$ in general.

Definition 2.2** (Neighborhood).**

A node $v\in V$ is called a neighbor of the node $u\in V$ if there exists an edge $(u,v)\in E$. We abbreviate this as $v\sim u$, which reads as “$v$ is a neighbor of $u$”. If on the other hand $v$ is not a neighbor of $u$, we use $v\not\sim u$. We further define the neighborhood $\mathcal{N}(u)$ of a vertex $u\in V$ as $\mathcal{N}(u)\coloneqq\{v\in V\colon v\sim u\}$. The degree of a vertex $u\in V$ is defined as the amount of its neighbors $\operatorname{deg}(u)=\lvert\mathcal{N}(u)\rvert$.

Finally, the weight function is the central element to model the significance of a relationship between two connected vertices with respect to an application dependent criterion. In many cases the weight function is chosen as similarity function based on the characteristics of the modeled entities, i.e., by the evaluation of vertex functions as defined below. Then the weight function $w$ takes high values for important edges, i.e., high similarity of the involved vertices, and low values for less important ones. In many applications one normalizes the values of the weight function by $w\colon E\rightarrow[0,1]$.

A natural extension of the weight function to the full set $V\times V$ is given by setting $w(u,v)=0$, if $v\not\sim u$ or $u=v$ for any $u,v\in V$. In this case the graph edge set can be characterized as $E=\{(u,v)\in V\times V\colon w(u,v)>0\}$.

In many cases it is preferable to use symmetric weight functions, i.e., $w(u,v)=w(v,u)$. This also implicates that $v\sim u\Rightarrow u\sim v$ holds for all $u,v\in V$. Hence, all directed graphs with symmetric weight function discussed in this paper could be interpreted as undirected graphs, though in the following it is important that each edge $(u,v)\in E$ has a start node $u$ which is different from its end node $v$.

2.2 Riemannian Manifolds

We denote by $\mathcal{M}$ a complete, connected, $m$-dimensional Riemannian manifold and its Riemannian metric by $\langle\cdot,\cdot\rangle_{x}\colon\mathrm{T}_{x}\mathcal{M}\times\mathrm{T}_{x}\mathcal{M}\to\mathbb{R}$, where $\mathrm{T}_{x}\mathcal{M}$ is the tangent space at $x\in\mathcal{M}$. We furthermore denote by $\lVert\cdot\rVert_{x}$ the induced norm and by $\mathrm{T}\mathcal{M}$ the disjoint union of all tangent spaces called the tangent bundle. In the following we introduce the necessary notations and theory in order to extend algorithms and methods for real-valued functions on finite weighted graphs to the case of manifold-valued functions. Further details on this subject can be found, e.g., in the textbooks on manifolds [49, 27, 1].

The completeness of $\mathcal{M}$ implies that any two points $x,y\in\mathcal{M}$ can be joined by a (not necessarily unique) shortest curve. We denote such a curve by $\gamma_{\overset{\frown}{x,y}}\colon[0,L]\to\mathcal{M}$, where $L$ is the length of a shortest curve. We denote by $\dot{\gamma}(t)\coloneqq\frac{d}{dt}\gamma(t)\in T_{\gamma(t)}\mathcal{M}$ the derivative vector field of the curve. We further require, that the curve is parametrized with constant speed, i.e., $\lVert\dot{\gamma}_{\overset{\frown}{x,y}}(t)\rVert_{\gamma_{\overset{\frown}{x,y}}(t)}=1$, $t\in[0,L]$. This generalizes the idea of shortest paths from the Euclidean space $\mathcal{M}=\mathbb{R}^{m}$, i.e., straight lines, to a manifold.

Let $D\colon\Gamma(\mathrm{T}\mathcal{M})^{2}\to\Gamma(\mathrm{T}\mathcal{M})$ denote the covariant derivative corresponding to the Levi-Citiva connection [27, Theorem 3.6], where $\Gamma(\mathrm{T}\mathcal{M})$ is the set of all differentiable vector fields on $\mathcal{M}$. A geodesic also fulfills that the covariant derivative $D_{\dot{\gamma}(t)}\dot{\gamma}=\frac{D}{dt}\dot{\gamma}$ of the (tangent) vector field $\dot{\gamma}(t)\in T_{\gamma(t)}\mathcal{M}$ vanishes. However, (shortest) geodesics are not the only curves possessing this property. For example let us fix on the two-dimensional sphere $\mathcal{M}=\mathbb{S}^{2}\coloneqq\{x\in\mathbb{R}^{3}:\lVert x\rVert=1\}$ two points $x,y\in\mathbb{S}^{2}$ not being antipodal. Then there exists a unique great circle containing the points $x,y$. One of the arcs is the shorter one and yields the geodesic $\gamma_{\overset{\frown}{x,y}}$. Still both great arcs (seen as curves on the manifold) have a vanishing covariant derivative. In the literature, the non-shortest curves with vanishing covariant derivative $D\dot{\gamma}=0$ are often also called geodesics. In the following we refer to geodesics as always being shortest ones. The shortest geodesic might also not be unique, e.g., for two antipodal points $x,y$ on $\mathbb{S}^{2}$ both connecting great half circles are geodesics and both are shortest ones. In this case we will still write the geodesic meaning that any of the shortest geodesics is meant.

The length of the shortest geodesic induces the geodesic distance $d_{\mathcal{M}}\colon\mathcal{M}\times\mathcal{M}\to\mathbb{R}$. From the Theorem of Hopf and Rinow, cf. [49, Theorem 1.7.1], we obtain that for some $\varepsilon>0$ and $\xi\in T_{x}\mathcal{M}$ there exists a unique geodesic $\gamma_{x,\xi}\colon(-\varepsilon,\varepsilon)\to\mathcal{M}$ fulfilling $\gamma_{x,\xi}(0)=x$ and $\dot{\gamma}_{x,\xi}(0)=\xi\in\mathrm{T}_{x}\mathcal{M}$. Furthermore, by the uniqueness the Hopf-Rinow Theorem states that the uniqueness is given for any $\varepsilon$ and hence with $\varepsilon=1$ and $y\coloneqq\gamma_{x,\xi}(1)$ we see that the geodesic $\gamma_{\overset{\frown}{x,y}}$ is just a reparametrization of $\gamma_{x,\xi}$, namely $\gamma_{x,\xi}(t)=\gamma_{\overset{\frown}{x,y}}(\lVert\xi\rVert_{x}t)$.

Definition 2.3** (Exponential and logarithmic maps).**

The exponential map $\exp_{x}\colon\mathrm{T}_{x}\mathcal{M}\to\mathcal{M}$ is defined as $\exp_{x}(\xi)=\gamma_{x,\xi}(1)$. Furthermore, let $r_{x}\in\mathbb{R}^{+}$ denote the injectivity radius, i.e., the largest radius such that $\exp_{x}$ is injective for all $\xi$ with $\lVert\xi\rVert_{x}<r_{x}$. Furthermore, let

[TABLE]

Then the inverse map $\log_{x}\colon\mathcal{D}_{x}\to\mathrm{T}_{x}\mathcal{M}$ is called the logarithmic map and maps a point $y=\gamma_{x,\xi}(1)\in\mathcal{D}_{x}$ to $\xi$.

By the properties of the exponential map it holds that $d_{\mathcal{M}}(x,y)=d_{\mathcal{M}}(x,\gamma_{x,\log_{x}y}(1))=\lVert\log_{x}y\rVert_{x}$.

For a differentiable vector field $V\in\Gamma(\mathrm{T}\mathcal{M})$ and a smooth curve $c(t)$ on the manifold we can further define the vector field along $c$ by $V(t)\coloneqq V(c(t))$. The solution $W\in\Gamma(\mathrm{T}\mathcal{M})$ of

[TABLE]

is called parallel transport of $W(0)\in T_{c(0)}\mathcal{M}$ along the curve $c$. We introduce the special notation $\operatorname{PT}_{x\to y}\colon\mathrm{T}_{x}\mathcal{M}\to\mathrm{T}_{y}\mathcal{M}$ for the parallel transport of a vector $\nu\in\mathrm{T}_{x}\mathcal{M}$ to $y$ along a geodesic connecting $x$ and $y$, i.e., $\operatorname{PT}_{x\to y}(\nu)\coloneqq W(d_{\mathcal{M}}(x,y))\in T_{y}\mathcal{M}$ is the parallel transport $W$ with $W(0)=\nu$ and $c(t)=\gamma_{\overset{\frown}{x,y}}$ evaluated at the end point $t=d_{\mathcal{M}}(x,y)$ of the geodesic, i.e., at $y$. The above introduced definitions are collectively illustrated in Figure 2.

Whenever the manifold is not explicitly given, i.e., there is only an approximate description available, the whole following framework can also be employed using numerical approximations of the exponential and logarithmic maps, hence also using only an approximation to the distance $d_{\mathcal{M}}$ and the parallel transport on the manifold.

Finally, we have to generalize the notion of convexity to a Riemannian manifold. A set $\mathcal{C}\subset\mathcal{M}$ is called convex if for any two points $x,y\in\mathcal{M}$ all minimizing geodesics $\gamma_{\overset{\frown}{x,y}}$ lie in $\mathcal{C}$. Such a set $\mathcal{C}$ is called weakly convex if for all $x,y$ there exists a geodesic $\gamma_{\overset{\frown}{x,y}}$ lying completely in $\mathcal{C}$. Note that convexity might be a local phenomenon on certain manifolds: on the sphere $\mathbb{S}^{2}$ there exists no convex set $\mathcal{C}$ larger than an open half-sphere. Finally, a function $f\colon\mathcal{M}\to\mathbb{R}$ is called (weakly) convex on a (weakly) convex set $\mathcal{C}$ if for all $x,y\in\mathcal{C}$ the composition

[TABLE]

is a convex function. This notion of convexity is sometimes called geodesic convexity.

3 A graph framework for manifold-valued functions

In this section we propose a graph-based framework for processing manifold-valued data. This approach allows us to translate successful variational models and partial differential equations to manifolds. Furthermore, the graph framework allows us to unify local as well as nonlocal methods for manifold-valued data. This leads to a comprehensive discussion of related models and elegant numerical solutions. We begin by introducing the concept of manifold-valued vertex functions and edge functions which take values in respective tangent spaces in Section 3.1. In this setting we also discuss the characteristics of the corresponding function spaces. One fundamental definition will be the notion of the discrete directional derivative for manifold values, which we introduce in Section 3.2. Based on this we introduce first order differential operators, i.e., the weighted local discrete gradient and the weighted local divergence operator. We use these definitions to derive higher order differential operators for manifold-valued vertex functions in Section 3.3, namely a family of isotropic and anisotropic graph $p$-Laplacians. Finally, we discuss the special case in which the manifold $\mathcal{M}$ is simply a Euclidean space $\mathbb{R}^{m}$. In this setting we show that our proposed framework is a generalization of well-known graph operators from the literature.

3.1 Discrete calculus for manifold-valued vertex functions

and tangent edge functions

In the following we introduce the basic definitions and observations for manifold-valued vertex functions and tangent edge functions on a finite weighted graph $G=(V,E,w)$.

Definition 3.1** (Manifold-valued vertex function).**

Let $G=(V,E,w)$ be a finite weighted graph and let $\mathcal{M}$ be a complete, connected, $m$-dimensional Riemannian manifold. We define a manifold-valued vertex function $f$ as

[TABLE]

For a given vertex function $f$ we can define the finite, disjoint union of the induced tangent spaces $\mathrm{T}_{f}\mathcal{M}$ as

[TABLE]

In the following we discuss properties of the function spaces of manifold-valued vertex functions. We define a metric for vertex functions based on mappings into respective tangent spaces employing the logarithmic map.

In the following we require that two adjacent data items $f(u),f(v),(u,v)\in E$ posess the property that $f(v)\in\mathcal{D}_{f(u)}$ and vice versa, such that the logarithmic map $\log_{f(u)}f(v)$ is well defined.

Definition 3.2** (Function spaces on vertex functions).**

*Let $G=(V,E,w)$ be a finite weighted graph. We introduce a metric space of manifold-valued vertex functions denoted

as $(\mathcal{H}(V;\mathcal{M}),d_{\mathcal{H}(V;\mathcal{M})}(\cdot,\cdot))$ on a set of admissible functions*

[TABLE]

with the associated metric, which is for two vertex functions $f,g\in\mathcal{H}(V;\mathcal{M})$ given by

[TABLE]

Here, $d_{\mathcal{M}}(\cdot,\cdot)\colon\mathcal{M}\times\mathcal{M}\rightarrow\mathbb{R}_{\geq 0}$ denotes the geodesic distance between two points on the manifold $\mathcal{M}$. In fact this is the standard metric on the product manifold $\mathcal{M}^{n}$.

We will refer to functions $f\in\mathcal{H}(V;\mathcal{M})$ as data with neighboring items fulfilling the locality property, i.e. neighboring data items with respect to the set of edges are within their injectivity radii. If this is not fulfilled for data given in practice, a denser sampling, i.e. larger graph with smaller edge weights is required.

Based on the notion of vertex functions mapping from the set of vertices $V$ into the manifold $\mathcal{M}$ we are able to define edge functions mapping from the set of edges $E$ into the respective tangent spaces. The edge functions in the vector valued setting, see, e.g. [33], represent finite differences. The analogue of a discrete difference as an approximation of the derivative is given by the logarithmic map on a manifold. The corresponding values are therefor given in the tangent bundle. Hence the definition of edge functions is only meaningful with respect to an associated vertex function as this induces the corresponding tangent spaces.

Definition 3.3** (Tangential edge function).**

Let $G=(V,E,w)$ be a finite weighted graph. We define a tangent edge function $H_{f}$ with respect to a manifold-valued vertex function $f\in\mathcal{H}(V;\mathcal{M})$ as

[TABLE]

Note that these tangent edge functions require directed graphs due to the tangent space $\mathrm{T}_{f(u)}\mathcal{M}$ of the function value $f(u)$ at the start point $u\in V$ of the edge involved. For undirected graphs, we consider their directed analogon with the symmetry property $(u,v)\in E\Leftrightarrow(v,u)\in E$ and can still investigate their tangent edge functions thereon.

We can introduce a family of measures on the values of an edge function $H$ at an edge $(u,v)\in E$ by using the respective vector norms which can be associated to the tangent space $\mathrm{T}_{f(u)}\mathcal{M}$.

Let $\mathcal{H}(E;\mathrm{T}_{f}\mathcal{M})$ denote the space of functions for which $H(u,v)\in T_{f(u)}\mathcal{M}$ holds true for any edge $(u,v)\in E$. Then we define a family of norms, namely for $p,q\geq 1$, by

[TABLE]

We further introduce the short hand notation $\lVert\cdot\rVert_{p,q}\coloneqq\lVert\cdot\rVert_{\ell_{p,q}(E;\mathrm{T}_{f}\mathcal{M})}$ if $\mathcal{M}$, $E$, and $f$ are clear from the context. Note that this includes the special cases of the isotropic norm $\lVert\cdot\rVert_{p,2}$ and the anisotropic norm $\lVert\cdot\rVert_{p,p}$. For the special case $p=q=2$ we obtain a Hilbert space, which is equipped with the inner product

[TABLE]

In order to discuss edge functions independently of an associated vertex function $f\in\mathcal{H}(V;\mathcal{M})$ we define a general Hilbert space as union of all Hilbert spaces of admissible edge functions by

[TABLE]

Then one can interpret an edge function $H\in\mathcal{H}(E)$ as a function $H\colon E\rightarrow\mathrm{T}\mathcal{M}$ with the restriction, that all $H(u,v)$, $v\sim u$, are values in the same tangent space.

Definition 3.4** (Local variation of an edge function).**

Let $G=(V,E,w)$ be a finite weighted graph. We define the local variation of a tangent edge function $H_{f}\in\mathcal{H}(E;\mathrm{T}_{f}\mathcal{M})$ at a vertex $u\in V$ as

[TABLE]

Note that this norm represents the summand for a fixed $u\in V$ within the $\ell_{p,q}$-norm  (1).

If we further employ the parallel transport, i.e., for $F,G\in\mathcal{H}(E)$ with $F\in\mathcal{H}_{f}(E)$, $G\in\mathcal{H}_{g}(E)$ we have for $\tilde{F}(u,v)\coloneqq\operatorname{PT}_{f(u)\to g(u)}F(u,v)$ that $\tilde{F}\in\mathcal{H}_{g}(E)$. This way, we may also obtain a Hilbert space.

We define a symmetric mapping on the space of vertex functions $\mathcal{H}(V;\mathcal{M})$ by

[TABLE]

The symmetry can be seen by noting that both sum run over the same set of nodes for $\langle g,f\rangle$ and a single summand reads $\langle\log_{g(u)}g(v),\operatorname{PT}_{f(u)\to g(u)}\log_{f(u)}f(v)\rangle_{g(u)}$. This can be be parallel transported to $\mathrm{T}_{f(u)}\mathcal{M}$ which does not change its value. The symmetry follows from the symmetry of the inner product in $\mathrm{T}_{f(u)}\mathcal{M}$.

We further equip the space $\mathcal{H}(V;\mathcal{M})$ with a norm, which is for $p=2$ induced by the symmetric map

[TABLE]

Note that this distance can also be employed for arbitrary functions $f$, i.e., also for data not fulfilling the locality property by employing the Riemannian distance $d_{\mathcal{M}}(f(u),f(v))^{p}$ instead of the inner products of the logarithmic maps. Hence indeed the sum of all distances $d_{\mathcal{M}}(f(u),f(v))$ between all nodes for $p=2$, and thus induces a proper norm on the set of vertex functions. However, due to the missing linearity, this norm is not induced by an inner product. We still employ the symbol $\langle\cdot,\cdot\rangle$ for this symmetric map for the sake of simplicity.

Finally, we can also look at tangent vertex functions $H\in\mathcal{H}(V;T_{f}\mathcal{M})$ associated with $f\in\mathcal{H}(V;\mathcal{M})$, i.e.,

[TABLE]

which can be equipped with the metric and its induced norm from the tensor product $\mathrm{T}_{f}\mathcal{M}$ of the tangent spaces $T_{f(u)}\mathcal{M}$, $u\in V$. Similarly we use $\mathcal{H}(V;\mathrm{T}\mathcal{M})$.

3.2 First-order difference graph operators

After defining all necessary functions, function spaces, and measures in the last section, we can now proceed to introduce novel first-order graph operators for manifold-valued data.

Definition 3.5** (Weighted directional derivative).**

Let $G=(V,E,w)$ be a finite weighted graph. We define the weighted directional derivative $\partial_{v}f\colon V\rightarrow\mathrm{T}_{f(u)}\mathcal{M}$ of a manifold-valued vertex function $f\in\mathcal{H}(V;\mathcal{M})$ at $u\in V$ in direction of another vertex $v\in V$ as

[TABLE]

Note that we set $\partial_{v}f(u)=0$, whenever $v\not\sim u$. From this definition we can directly deduce the following lemma.

Lemma 3.6** (Properties of the weighted directional derivative).**

*The weighted directional

derivative is*

1. i)

reflexive, i.e., $\partial_{v}f(u)=0$ for $u=v$ 2. ii)

anti-symmetric under parallel transport, i.e., it holds

[TABLE]

We also define another variant of the weighted directional derivative as an edge function.

Definition 3.7** (Weighted local gradient).**

Let $G=(V,E,w)$ be a finite weighted graph. We define the weighted local gradient $\nabla\colon\mathcal{H}(V;\mathcal{M})\to\mathcal{H}(E)$ of a manifold-valued vertex function $f\in\mathcal{H}(V;\mathcal{M})$ as

[TABLE]

Clearly, we see that $\nabla f\in\mathcal{H}(E;\mathrm{T}_{f}\mathcal{M})$. Furthermore, the following Theorem derives a relationship between the weighted local gradient $\nabla f$ and a corresponding edge function $H_{f}$.

Theorem 3.8**.**

Let $G=(V,E,w)$ be a finite weighted graph with a symmetric edge set $E$, i.e. $(u,v)\in E\Leftrightarrow(v,u)\in E$ and an arbitrary, not necessarily symmetric weight function $w$. Let further $f\in\mathcal{H}(V;\mathcal{M})$ be a vertex function and $H_{f}\in\mathcal{H}(E;\mathrm{T}_{f}\mathcal{M})$ a tangent edge function. Then we have the following relationship

[TABLE]

Proof.

Starting with the left hand side we compute

[TABLE]

∎

We use the relationship deduced in Theorem 3.8 to motivate the definition of a local divergence operator for the respective tangent spaces. Due to the fact that the discussed manifold-valued vertex functions are not mapping into a vector space, we are not able to deduce the divergence as an adjoint operator of the gradient operator as it is the case, e.g., for Euclidean spaces.

Definition 3.9** (Weighted local divergence).**

We define a weighted local divergence operator $\operatorname{div}\colon\mathcal{H}(E)\to\mathcal{H}(V;T\mathcal{M})$ for a tangent edge function $H_{f}\in\mathcal{H}(E;\mathrm{T}_{f}\mathcal{M})$ associated to $f\in\mathcal{H}(V;\mathcal{M})$ at a vertex $u\in V$ as

[TABLE]

Remark 3.10**.**

Let $G=(V,E,w)$ be a finite weighted graph with symmetric weight function, i.e., $w(u,v)=w(v,u)$ for all $(u,v)\in E$. Assuming an edge function $H_{f}\in\mathcal{H}(E;\mathrm{T}_{f}\mathcal{M})$ which is anti-symmetric under parallel transport, i.e., $H_{f}(u,v)=-\operatorname{PT}_{f(v)-f(u)}H_{f}(v,u)$ for all $(u,v)\in E$, we obtain a concise representation of the weighted local divergence operator as

[TABLE]

3.3 A family of graph $p$-Laplace operators for manifold-valued functions

Based on the definitions of the weighted local gradient and the weighted local divergence in Section 3.2, we are able to introduce a family of graph $p$-Laplace operators for manifold-valued functions. Note that for $p,q<1$ in (1) and $p<1$ in (2) we have only quasi-norms; however, we still can define the following operators for this case. For the sake of simplicity we will discuss all operators in the case of finite weighted graphs with symmetric weight function in the following. However, these operators can easily be derived for arbitrary weight functions using the tools introduced in Section 3.2 above.

Definition 3.11** (Graph $p$-Laplacian operators for manifold-valued functions).**

Let $G=(V,E,w)$ be a finite weighted graph with symmetric weight function $w$ and let $0<p<\infty$. For a manifold-valued vertex function $f\in\mathcal{H}(V;\mathcal{M})$ we define the anisotropic graph $p$-Laplacian $\Delta^{\mathrm{a}}_{p}\colon\mathcal{H}(V;\mathcal{M})\rightarrow\mathcal{H}(V;T\mathcal{M})$ at a vertex $u\in V$ as

[TABLE]

and the isotropic graph $p$-Laplacian $\Delta^{\mathrm{i}}_{p}\colon\mathcal{H}(V;\mathcal{M})\rightarrow\mathcal{H}(V;T\mathcal{M})$ at a vertex $u\in V$ as

[TABLE]

For the special case of $p=2$ we obtain in both above definitions an operator $\Delta\colon\mathcal{H}(V;\mathcal{M})\rightarrow H(V;T_{f}\mathcal{M})$, which we denote as graph Laplacian for manifold-valued functions, by

[TABLE]

3.4 Special case of $\mathcal{M}=\mathbb{R}^{m}$

In the following we discuss the special case in which the Riemannian manifold $\mathcal{M}$ is simply the Euclidean space $\mathbb{R}^{m},m\in\mathbb{N}$. This setting is ubiquitous in many applications such as image and point cloud processing. As $\mathcal{M}$ has curvature zero everywhere we can demonstrate that the introduced calculus from Subsection 3.3 leads to equivalent definitions of well-known graph operators from the literature, e.g., as used in [32, 33, 43]. Thus, the here proposed framework can be interpreted as a generalization of the common calculus for vector-valued vertex functions $f\colon V\rightarrow\mathbb{R}^{m}$.

For $\mathcal{M}=\mathbb{R}^{m}$ it gets clear that the respective tangent spaces are equal, i.e., $T_{f(u)}\mathcal{M}=\mathbb{R}^{m}$ for all $u\in V$, and hence all mathematical operations are defined globally in this case. In particular, the maps $\log\colon\mathcal{M}\rightarrow\mathrm{T}\mathcal{M}$ and $\exp\colon\mathrm{T}\mathcal{M}\rightarrow\mathcal{M}$ introduced in Section 2.2 become globally defined linear operators which correspond to vector subtraction and addition in $\mathbb{R}^{m}$. Furthermore, the parallel transport reduces to the identity operator because the tangent spaces coincide with the tangent space at the origin, which is itself isomorphic to $\mathbb{R}^{m}$. These observations lead to a simplification of all introduced operators above, which are be summarized in the following Corollary.

Corollary 3.12**.**

*Let $\mathcal{M}=\mathbb{R}^{m},f\in\mathcal{H}(V;\mathcal{M})$ a vertex function, and $\xi\in\mathrm{T}_{f(u)}\mathcal{M}=\mathbb{R}^{m}$ a (tangent) vector. For vertices $u,v\in V$ we obtain

1. i)

Basic mathematical operations

[TABLE] 2. ii)

First-order differential operators

[TABLE] 3. iii)

Graph $p$-Laplacian operators

[TABLE]

Note that for $\mathcal{M}=\mathbb{R}^{m}$ the weighted local divergence operator becomes a linear operator and thus the relationship in Theorem 3.8 can be further elaborated leading to the well-known fact that the negative divergence operator is the adjoint graph operator of the weighted gradient, see e.g., [32].

Corollary 3.13**.**

Let $G=(V,E,w)$ be a finite weighted graph with symmetric weight function $w$ and let $\mathcal{M}=\mathbb{R}^{m}$. Then for any vertex function $f\in\mathcal{H}(V)$ and any edge function $H\in\mathcal{H}(E)$ the following relationship holds:

[TABLE]

for which the inner products above are the standard inner products of $\mathbb{R}^{m\times n\times n}$ and $\mathbb{R}^{m\times n}$, respectively.

4 Formulation of variational problems

In this section we develop different algorithms to solve mathematical problems for manifold-valued vertex functions on graphs. Based on the introduced graph operators in Section 3 we translate important PDEs and variational models from continuous mathematics to graphs in the tradition of [32, 43] and for the first time formulate image processing problems for manifold-valued functions uniformly in a local and nonlocal setting. These image processing problems include segmentation, inpainting, and denoising of manifold-valued data. Since in this section we are interested in demonstrating the feasibility of our approach, we will restrict ourselves to the latter task. Thus, in the following we will formulate a class of denoising tasks as variational minimization problems on graphs in Section 4.1. Subsequently, we will discuss necessary optimality conditions for minimizers of these problems in Section 4.2, which yield the introduced graph $p$-Laplace operators for manifold-valued functions, $p\geq 1$. To compute respective minimizers we derive two different numerical schemes in Section 4.3.

4.1 Problem formulation

Let $G=(V,E,w)$ be finite weighted graph with symmetric weight function $w\colon E\rightarrow[0,1]$, i.e., $w(u,v)=w(v,u)$ for all $u,v\in V$. Furthermore, let $f_{0}\colon V\to\mathcal{M}$ be a manifold-valued vertex function which models the given (perturbed) data. One possibility is to assume that $f_{0}$ is an observation of the following data formation process, see e.g., [53]:

[TABLE]

Here, $f_{0}$ is modeled as a noisy variant of the unknown data $\hat{f}\in\mathcal{H}(V;\mathcal{M})$, which is altered by some perturbation $\varepsilon\in T_{\hat{f}}\mathcal{M}$. Recovering the original noise-free data $\hat{f}$ from $f_{0}$ in (4.1) is a common example of an inverse problem. In order to guarantee the well-posedness of this problem one may incorporate a-priori knowledge about an unknown solution $f$, e.g., smoothness, and pose the problem from a statistical view as a maximum a-posteriori (MAP) estimation method. Typically, one aims to find minimizers of a convex energy functional $\mathcal{E}\colon\mathcal{H}(V;\mathcal{M})\rightarrow\mathbb{R}$ consisting of data fidelity and regularization terms, i.e.,

[TABLE]

For further details on deriving convex variational problems from a statistical modeling perspective we refer to [71].

In the following we discuss variational models for denoising of manifold-valued vertex functions. We begin introducing a family of anisotropic energy functionals for $p\geq 1$ to be optimized:

[TABLE]

Within the data fidelity term $d_{\mathcal{H}(V;\mathcal{M})}$ denotes the metric of the space of vertex functions $\mathcal{H}(V;\mathcal{M})$ measuring the distance of a function $f$ to the given data $f_{0}$, $\lambda>0$ is a fixed regularization parameter controlling the smoothness of the solution, and the regularization term on the right side of (4) denotes the discrete anisotropic Dirichlet energy as a measure of local variance. The denoising task is now to find a solution $f^{\star}\in\mathcal{H}(V;\mathcal{M})$ to the following optimization problem:

[TABLE]

Furthermore, we are interested in a family of energy functionals of the form

[TABLE]

The difference to the previous model in (4) is the exchange of the regularization term by an isotropic norm. Analogously as before, this results in optimization problems of the form

[TABLE]

Note that the optimization problems (5) and (7) cover two interesting special cases. For $p=2$ both formulations are equivalent and can be interpreted as Tikhonov-regularized denoising problem, which aims to reconstruct smooth solutions $f^{\star}$. For $p=1$ the problem (5) corresponds to anisotropic total variation-regularized denoising, which has been investigated recently in the context of manifold-valued functions in [79, 58, 17] using lifting techniques and the cyclic proximal point algorithm. On the other hand (7) yields for $p=1$ an isotropic total variation-regularized denoising formulation, which has up to now only been tackled by half-quadratic minimization in [15].

4.2 Optimality conditions

To perform denoising we need to find minimizers of the discrete energy functionals introduced in Section 4.1. Note that while these energies are convex for $p\geq 1$ on the Euclidean space $\mathcal{M}=\mathbb{R}^{m}$, in general this only holds locally on manifolds. As mentioned in Section 2.2 even sets might not have a global notion of convexity on manifolds. However, when restricting to the case of Hadamard manifolds, i.e., manifolds of nonpositive sectional curvature, all important properties carry over, i.e., convexity and even lower semi-continuity and coerciveness. For details on optimization on Hadamard manifolds we refer to the monograph [10]. In the case of Hadamard manifolds we may thus conclude the existence of respective minimizers $f^{\star}\in\mathcal{H}(V;\mathcal{M})$ of (5) and (7). Note that for general manifolds the minimizers discussed in the following might only be local minimizers, since there is no general notion of global convexity. Locally, we introduced (geodesic) convexity in Section 2, which we employ in the following.

To compute a minimizer of an energy functional $\mathcal{E}\colon\mathcal{H}(V;\mathcal{M})\rightarrow\mathbb{R}$ we briefly introduce the notion of a derivative and a subdifferential of $\mathcal{E}$, i.e., which fall back to subdifferentials on manifolds. Since one is not able to perform basic arithmetic operations on $\mathcal{H}(V;\mathcal{M})$, e.g., addition of functions or scalar multiplication, a gradient of the energy functional $\mathcal{E}$ has to be defined on the set of tangent spaces $\mathrm{T}\mathcal{M}$, see, e.g., the text book [54].

Definition 4.1** (Differential and gradient).**

Let us denote by $\mathcal{C}^{\infty}(\mathcal{M},\mathbb{R})$ all smooth maps from $\mathcal{M}$ to $\mathbb{R}$. Furthermore, let $F\in\mathcal{C}^{\infty}(\mathcal{M},R)$ and let $\gamma\colon(-\varepsilon,\varepsilon)\to\mathcal{M}$ be a curve with $\gamma(0)=x\in\mathcal{M}$ and $\dot{\gamma}(0)=\xi\in\mathrm{T}_{x}\mathcal{M}$. Then the differential $D_{x}F\colon\mathrm{T}_{x}\mathcal{M}\to\mathbb{R}$ of $F$ is given by

[TABLE]

The gradient $\nabla_{\mathcal{M}}F(x)$ can be characterized by

[TABLE]

However, for the interesting case of $p=1$ in (5) and (7) classical differentiability is too restrictive for obtaining suitable minimizers. For this reason we adapt the notion of a subdifferential from [35] to our setting.

Definition 4.2** (Subdifferential).**

Let $F\colon\mathcal{M}\to\mathbb{R}$ be a (locally) convex function and $x\in\mathcal{M}$. Then the subdifferential of $F$ at $x$ is defined by

[TABLE]

An element $\xi\in\partial F(x)$ of the subdifferential of $F$ at $x$ is called subgradient.

If the subdifferential is a singleton the subgradient is unique and equals the gradient $\nabla_{\!\mathcal{M}}$. Before explicitly deriving the necessary optimality conditions for the introduced variational denoising models on manifold-valued functions in Section 4.3, we would like to recall the following useful result from [3].

Lemma 4.3**.**

Let $y\in\mathcal{M}$ denote a fixed point on the manifold and let $F_{p}\colon\mathcal{M}\to\mathbb{R}$ be a function with $F_{p}(x)=d^{p}(x,y)$ for $p\geq 1$. Then we have

[TABLE]

Furthermore, for $p=1$ and $x=y$ we have $\partial F_{1}(x)=B_{x,1}$, where $B_{x,r}\coloneqq\{\xi\in\mathrm{T}_{x}\mathcal{M}\ :\ \lVert\xi\rVert_{x}\leq r\}\subset\mathrm{T}_{x}\mathcal{M}$.

The second case of the subgradient can be derived by looking at the subdifferential and observing that it only contains [math]. The last subdifferential follows from the definition of the subdifferential, cf. [35],

[TABLE]

in a certain ball $\mathcal{M}_{x}\subset\mathcal{D}_{x}\subset\mathcal{M}$ around $x$ where the exponential map is injective and $\xi=\frac{1}{d_{\mathcal{M}}(x,z)}\log_{x}z$ for $x\neq z$. Note that especially $0\in\partial F_{p}(x)$ for $x=y$ and $p\geq 2$.

Based on this result we are able to derive the conditions for a minimizer of the anisotropic energy functional (5). Note that in this discrete setting the energy functional is in fact a function defined on the (product) manifold $\mathcal{M}^{\lvert V\rvert}$. We derive

[TABLE]

Note that only in the case $p=1$ and $f(u)=f(v)$, $u,v\in V$, the subdifferential in the second summand is not a singleton. As discussed above the subdifferential however contains [math]. Furthermore, in the case $p<2$ and $f(u)=f(v)$ we have defined the anisotropic graph $p-$Laplacian as [math]. Bearing these special cases in mind we can further simplify the optimality conditions as

[TABLE]

This leads to the following PDE on a finite weighted graph as necessary condition for a minimizer of (5):

[TABLE]

Following analogously the argumentation above we are able to derive the necessary conditions for a minimizer of the isotropic energy functional (7) as follows:

[TABLE]

Hence we obtain

[TABLE]

This leads to the following PDE on a finite weighted graph as necessary condition for a minimizer of (7):

[TABLE]

4.3 Numerical optimization schemes

In the following we derive numerical minimization schemes to approximate solutions to the necessary optimality conditions of both the anisotropic and the isotropic denoising problems discussed in Section 4.2. In particular, we will discuss two different iterative methods converging to solutions of the PDEs in (8) and (9). For the sake of simplicity we discuss both approaches collectively by introducing an operator $\Delta_{p}\colon\mathcal{H}(V;\mathcal{M})\rightarrow\mathcal{H}(V;T\mathcal{M})$ as place holder for the anisotropic and isotropic $p$-Laplacian operators. We only distinguish between the two formulations when giving explicit computation formulas at the end of each paragraph.

Explicit scheme by time-discretized evolution.

The basic idea of our first numerical approach is to extend the domain of our vertex functions by an additional artificial time dimension and subsequently compute a steady-state solution of the resulting time-dependent PDEs, see e.g., [74]. Thus, the aim is to find a stationary solution of the parabolic PDEs with initial value conditions of the form

[TABLE]

i.e., we try to compute a solution $f\in\mathcal{H}(V\times(0,\infty);\mathcal{M})$ such that

[TABLE]

Note that finding a stationary solution to (10) directly yields a minimizers of the original energy functionals problems (4) and (6) as the left side becomes zero. Setting $\lambda=0$ the problem (10) covers two well-known special cases. For $p=2$ one gets a diffusion equation for manifold-valued vertex functions which is equivalent to the discrete heat equation for a local neighborhood in the finite weighted graph $G$. For $p=1$ one gets the well-known total variation flow.

To solve the initial value problem (10) we discretize the time derivative on the left hand side of the equation using an forward difference scheme and assuming the right hand side of the equation to be known from the last time step, i.e.,

[TABLE]

Here, $\Delta t$ denotes the step size of the evolution process induced by the time discretization. This gives us the following explicit formulas for the iterative computation of solutions to the anisotropic and isotropic denoising problems,

[TABLE]

and

[TABLE]

respectively.

Note that this Euler forward time discretization realizes a gradient descent minimization algorithm. Although this approach is very simple to implement, it leads for $p<2$ to very strict Courant-Friedrich-Lewy (CFL) conditions in order to guarantee the stability of this numerical scheme. This means that updates $f^{n}\rightarrow f^{n+1}$ can only be performed with very small step size $\Delta t$ and hence convergence to possible solutions may need a huge amount of iterations. To alleviate this problem one may consider higher-order schemes, e.g., Runge-Kutta discretization schemes. A general CFL condition for a broad class of graph operators can be found in [33] for the case $\mathcal{M}=\mathbb{R}^{n}$.

Semi-implicit scheme by linearization.

Another approach, which we discuss in the following, is based on the idea of approximating the original problem (10) in every iteration step and solving a simple linear equation system. To perform this linearization we first reformulate (10) to an equivalent problem by inserting two zeros terms $\log_{f(u)}f(u)=0$.

For the anisotropic case we get

[TABLE]

By introducing

[TABLE]

we can write this more concisely as

[TABLE]

Analogously for the isotropic case and writing

[TABLE]

we obtain

[TABLE]

We can investigate both the isotropic and the anisotropic model in the same manner by introducing $b(u,v)$ which refers to exactly one of the above introduced coefficients. The discussed problems can be linearized by assuming that $b(u,v)$ is known (from the last iteration $f_{n}$) and approximating the introduced zero terms as $\log_{f(u)}f(u)\approx\log_{f_{n+1}(u)}f_{n}(u)$. Hence, we need to solve a linear equation system $Bx=c$. Note that this corresponds to a semi-implicit scheme for the original optimality conditions (8) and (9), respectively. For the case of graph constructions with relatively small neighborhoods, e.g., for local grid graphs or $k$-nearest neighbor graphs with small integers $k$, the resulting matrix $B$ is very sparse. Hence, it is advisable to use iterative methods to approximate solutions of the linear equation system.

By applying Jacobi’s method we get the following relationship in $T_{f_{n}(u)}\mathcal{M}$:

[TABLE]

This leads to the iterative update formula

[TABLE]

for which the terms $b(u,v)$ are given as above.

We would like to give some additional notes on the proposed numerical minimization schemes discussed above.

In order to avoid division by zero we use a relaxation of the terms $d_{\mathcal{M}}^{p-2}$ for $0<p<2$ by adding a small term $\varepsilon=10^{-7}$. Thus, we smooth our regularization functionals at zero. 2. 2.

The prefactor $b_{\mathrm{i}}(u,v)$ of the isotropic graph $p$-Laplacian in (13) is equal for all neighbors $v\sim u$ and hence can be precomputed once per vertex. 3. 3.

Although Jacobi’s method is in general not as efficient as related methods, e.g., the Gauss-Seidel method or successive overrelaxation methods, it has the advantage that one is able to compute all $\log$ terms and all geodesic distances $d_{\mathcal{M}}$ needed in (14) in a vectorial manner, which gives a higher overall computational performance when exploiting vectorization techniques. 4. 4.

For small values of $\lambda$ tending to zero we observed that the proposed semi-implicit scheme tends to be less robust than the proposed explicit scheme. This is likely due to the reason that the employed linearization for Jacobi’s method is not valid in cases with almost no data fidelity. Hence, for small values of the regularization parameter $\lambda$ we will use the explicit scheme in our experiments. 5. 5.

While the derivation in this Section is restricted to $p\geq 1$, the numerical schemes can also be applied for the case $p\geq 0$. The investigation of the quasi-norm cases is a point of future work.

5 Applications

In the following we present several examples illustrating the large variety of problems that can be tackled using the proposed manifold-valued graph framework. We begin by discussing our general experiment setup. We first focus on image domains in Section 5.1, i.e., manifold-valued data defined on regular grids and investigate the evolutionary flow for a local graph modeling spatial pixel neighborhoods. Furthermore, we compare our framework for the special case of nonlocal denoising of phase-valued data to a state-of-the-art method. Finally, we demonstrate a real-world application from denoising surface normals in digital elevation maps from LiDAR data. Subsequently, we model manifold-data measured on samples of an explicitly given surface and in particular illustrate denoising of diffusion tensors measured on a sphere in Section 5.2. Finally, we investigate denoising of real DT-MRI data from medical applications in Section 5.3 both on a regular pixel grid as well as on an implicitly given surface.

All algorithm were implemented in Mathworks Matlab by extending the open source software package Manifold-valued Image Restoration Toolbox (MVIRT)333open source, available at www.mathematik.uni-kl.de/imagepro/members/bergmann/mvirt/.

Patch distances on manifolds.

Similar to the real- or vector-valued case we can employ a patch-based distance function. Let $f\in\mathcal{M}^{n,m}$ be a manifold-valued image of dimension $n\times m$. Let us further denote by $f_{i}^{s}\in\mathcal{M}^{2s+1,2s+1}$ the $2s+1\times 2s+1$-sized patch with center $f_{i}^{s}(s+1,s+1)=f(i_{1},i_{2}),i=(i_{1},i_{2})\in\mathcal{G}\coloneqq\{1,\ldots,n\}\times\{1,\ldots,m\}$. For constructing patches at the boundary of the domain we assume periodic boundary conditions if $i_{1}-s<0,i_{2}<0,i_{1}+s>n,$ or $i_{2}+s>m$.

Then the patch-based similarity measure for manifold-valued data is given by

[TABLE]

Using this similarity measure we can introduce the (nonlocal) $\varepsilon$-ball graph $G_{\text{nl}}$ of an image $f$, where $V=\mathcal{G}$, i.e., every pixel is a vertex and $E_{\varepsilon}=\{(i,j)\in V\times V:1-\operatorname{PSM}(i,j)<\varepsilon\}$. Furthermore, we can also introduce the (nonlocal) $k$-neighbor graph, where $i\in V$ is connected via an edge to its $k$ most similar patches.

Evaluation.

For all our experiments we performed a parameter search with respect to the regularization parameter $\lambda$ in the models (4) and (6) in order to minimize the mean squared error (MSE), which is given by

[TABLE]

for which $\hat{f}$ denotes the original data without any perturbations.

5.1 Manifold-valued images

Evolution equations for the graph $p$-Laplacian.

We first take a look at a toy example of $\mathbb{S}^{2}$-valued data by taking the image S2WhirlImage444created by J. Persch shown in the first column of Figure 3. In this picture a few squares consisting of ‘whirl’ are depicted on a smoothly varying background. For both kinds of whirls, one being on the northern hemisphere clockwise (yellow) and the other on the southern hemisphere anti-clockwise (blue), the center point is the antipodal point, i.e., the south and north pole, respectively. To look at the $p$-Laplace flow we set $\lambda=0$ in the variational denoising models (4) and (6). We construct a local grid graph by connecting each pixel with his four direct neighbors. We then employ the explicit schemes (12) and (11), because the linearization introduced for the Jacobi iterations is not robust for small $\lambda$ and is not valid for $\lambda=0$. We employ Neumann boundary conditions in the case of local neighborhoods.

In Figure 3 the evolution of solutions over time in the explicit scheme is shown for the iterations $i=0$ (starting conditions), $i=1\,000$, $i=10\,000$, and $i=100\,000$. We set the time step size for all experiments of the explicit algorithm for $p\in\{1,2\}$ to $\Delta t=10^{-3}$ and for $p=0.1$ to $\Delta t=10^{-4}$. We further changed $\varepsilon=10^{-4}$ for the last case to improve numerical stability. We observe that for all six experiments the algorithm numerically converges to a constant solution. For the Tikhonov case, $p=2$, (first row) both the anisotropic and isotropic Laplace flow coincide and the convergence to a constant image is slower compared to the other rows, it reaches a constant solution after $i=1\,000\,000$ iterations. Due to the quadratic regularization term all edges are smoothed nicely, as can be best observed in iteration $i=1\,000$ in the first row. In fact, at this point of the iterations, the isotropic results are quite similar, comparing the first, second and fourth row for $p=2,1,0.1$, respectively. They are, however, not equal and tend to different constant solutions. The reason for this smooth appearance is due to the fact that we constructed a local graph consisting of four direct neighbors instead of using only two neighbors, as classically done in image processing.

On the other hand for the anisotropic TV case, $p=1$, (third row) edges are visibly preserved; see again iteration $i=1\,000$. In this case the flow field yields a nearly piecewise constant regions in the intermediate results. Furthermore, the anisotropic TV flow first joins regions having only small jumps (in geodesic distance) in between, see, e.g., the top left whirls at the bottom right corner after $i=1\,000$ iterations.

Finally, for the anisotropic $p=0.1$ case (last row) the intermediate results are even more piecewise constant compared to the anisotropic TV case, yielding a quite sparse gradient as one would expect for $p<1$.

Nonlocal denoising of phase-valued data.

As a first data-driven denoising example let $\mathcal{M}=\mathbb{S}^{1}$, i.e., we are considering an image of size $256\times 256$ pixels, where each pixel of $f_{0}\in(\mathbb{S}^{1})^{256\times 256}$ is a phase value living on the unit circle. This data occurs for example in interferometric synthetic aperture radar (InSAR), where height maps are generated from two parallel measurements with a laser of a certain wavelength. Their inference often results in a noisy phase shift, i.e., a height map “modulo wavelength”. For an example of real-valued data see Figure 1(a). An artificial phase valued image was introduced in [16], see Figure 4(a)) and further used in [53]. For our experiments we use the same data, which has been perturbed by wrapped Gaussian noise with $\sigma_{\text{n}}=0.3$ as shown in Figure 4(b)). The best result with respect to the MSE using second order statistics [53] is shown in Figure 4(c)). Indeed, the $k-$nearest neighbor construction in [53] is similar to ours as the authors employ the same patch distance. However, our edge weight function is different. Based on the $k$-nearest-neighbors for each vertex we set the weight function for the most similar neighbor to $1$ and for the least similar neighbor to [math] and interpolate the weight function linearly in between for the other neighbors. We optimized the involved parameters $s,k$, and $\lambda$ with respect to the MSE in the range from patch size $s\in\{2,\ldots,7\}$, $\lambda\in\tfrac{1}{512}\{1,\ldots,8\}$, and $k\in\{5,\ldots,15\}$ neighbors. We employ our explicit scheme with step size $\Delta t=10^{-4}$. As can be seen in Figure 4 our reconstruction result can compete with the best result from the state-of-the-art method in [53]. While the background and also small patches in the top right corner in Figure 4(c)) suffer from non-constant reconstructions, our approach reconstructs constant regions nearly perfectly due to the TV regularization. Only the the most left patch in blue suffers from fuzzy edges since there were not enough patches being similar in the data. Furthermore some small areas along the yellow part of the ellipsoidal region seem also to suffer from too few similar patches, yielding a slightly worse MSE for our reconstruction.

Denoising surface normals.

Next, we discuss a real world application based on topological surface data from the National Elevation Dataset (NED) [42], which was already used for manifold-valued data processing in [56, 58] with different denoising approaches, e.g., lifting for the TV-based regularization approach proposed in [58]. The provided digital elevation maps (DEM) are generated by light detection and ranging (LiDAR) measurements of earth’s surface topology. The given DEMs consist of measured surface heights and surface patch normals corresponding to the heights. The resolution of the image is quite low and hence one height and surface normal pair is a measurement from a relatively large area of the surface. The normals are usually used to color the data correctly and thus enhance the visual 3D perception of it. The original data often gives a perturbed visual impression (see Figure 5(a))) due to measurement noise and the complexity of the underlying real surface, which is only sparsely sampled during data acquisition. The provided data in the following has a size of only $40\times 40$ pixels, while covering a large area of a mountain formation. The processing task is now to reconstruct the measured surface normals, which live on the manifold $\mathbb{S}^{2}$, such that the visual impression of the DEM is enhanced, while preserving important topological features, e.g., kinks, and non-smooth regions with high curvature.

As can be seen in Figures 5(b)) and 5(c)) using the regular graph Laplacian, i. e. $p=2$, and the explicit scheme with $\Delta t=10^{-4}$, effectively reduces the noise of the surface normals and thus improves the impression of the respective visualizations. However, for $p=2$ one is not able to preserve the predominant kink of the mountain ridge. For this reason, we demonstrate in Figures 5(d)) and 5(e)) the effect of reconstructing the DEMs with a total variation-based regularization functional for $p=1$ in the cases of the anisotropic and isotropic graph $p$-Laplacian, respectively. While the top of the mountain ridge is now well-preserved the resulting normals of the anisotropic TV regularization are nearly piecewise constant, which leads to a rather blocky visual impression in Figure 5(d)). The isotropic TV regularization in Figure 5(e)) on the other hand succeeds in preserving geological details of the surface, while smoothing the surface normals for an improved visual impression. As a final comparison we show the results of non-convex regularization for $p=\tfrac{1}{2}$ in the anisotropic graph $p$-Laplacian. It can be observed that the resulting surface normals field has only very few directions and hence is relatively sparse. This leads to the visual impression of very sharp edged in Figure 5(f)).

5.2 Diffusion tensors on the sphere

To process manifold-valued data defined on a manifold itself, we investigate diffusion tensors defined on the unit sphere $\mathbb{S}^{2}$. For sampling we use the quadrature points of Chebyshev-type, which were computed in [46, Example 6.19] corresponding to polynomial degree $30$, i.e., on $M=480$ points. These are available on the homepage of the software package accompanying the thesis555see homepage.univie.ac.at/manuel.graef/quadrature.php [46]. We assign to each sampling point a symmetric positive definite matrix, which yields the original (unperturbed) dataset shown in Figure 6(a)). The illustrated ellipsoids on the sphere are colored with respect to the geodesic anisotropy, cf. [63, p. 290], employing the colormap viridis. This artificial data set is then perturbed by i.i.d. Gaussian noise, see e.g. [53], with standard deviation $\sigma=\frac{1}{4}$, which yields the noisy data in Fig. 6(b)). Each of the sampling points on the sphere $p_{i}\in\mathbb{S}^{2}$, $i=1,\ldots,480$, constitutes one vertex $v_{i}\in V$ of the finite weighted graph $G=(V,E,w)$. We construct an $\varepsilon-$ball graph by connecting to vertices $v_{i},v_{j}\in V$ by an edge $(v_{i},v_{j})\in E$ if $d_{\mathbb{S}^{2}}(p_{i},p_{j})\leq\varepsilon$ with $\varepsilon=\frac{\pi}{12}$ for this experiment. This results in $3648$ edges, i.e., each node has in average $7.6$ neighbors. The edges are weighted with the arc length distance on $\mathbb{S}^{2}$ of their incident nodes, i.e., $w(i,j)=d_{\mathbb{S}^{2}}(p_{i},p_{j})^{-2}$.

We then employ the aforementioned Jacobi iteration methods with $N=10,000$ iterations to compute a denoised reconstruction of the noisy data using the anisotropic graph $p-$Laplace for $p=1$ and $p=2$ with parameters $\lambda=23$ and $\lambda=80$, respectively. These parameters where chosen among all integers $\lambda\in\mathbb{Z}\cap\{1,\ldots,100\}$ such that the MSE is minimized. The results of this approach are shown in Figure 6(c)) and 6(d)). While the TV regularization for $p=1$ preserves edges, the Tikhonov regularization for $p=2$ introduces smooth reconstructions, which leads to higher values of the MSE.

5.3 Diffusion tensor imaging

A special case of data taking values on the Riemannian manifold $\mathcal{M}=\mathcal{P}(3)$ of symmetric positive definite $3\times 3$-matrices are diffusion tensors as obtained in diffusion tensor magnetic resonance imaging (DT-MRI). For this experiment we use the open data set from the Camino project666see http://camino.cs.ucl.ac.uk [29]. The provided data $f_{0}\in(\mathcal{P}(3))^{112\times 112\times 50}$ of a human brain consists of $50$ transversal slices, each contained in an array of size $112\times 112$. Values outside the brain are masked and take zero-matrices as values.

Nonlocal denoising.

In our first experiment on the Camino data we extract the transversal slice $z=35$, see Figure 7(a)). Analogously to the artificial InSAR experiment with phase-valued data in Section 5.1, we perform nonlocal TV-based denoising based on a $k-$nearest neighbor graph $G=(V,E,w)$ constructed for $k=10$ from patches of size $2s+1$, $s\in\{1,2\}$, cf. (15). We use the same linearly interpolated weight function for nonlocal denoising as already described for the artificial InSAR experiment discussed above. This experiment is the first to apply nonlocal TV-based approaches for manifold-valued data on a real medical dataset. Since the image consists also of masked data values (the surrounding of the brain), those pixels are ignored during computation of the patch-distance the $k-$nearest neighbors. For minimization we employ the semi-implicit minimization scheme using Jacobi iterations from (14). In Figure 7(b))–7(e)) the results of the anisotropic denoising approach are shown with a maximal number $N=100$ of iterations unless the average relative change $\frac{1}{\lvert V\rvert}\sum_{u\in V}d_{\mathcal{M}}(f_{n-1},(u),f_{n}(u))<10^{-7}$ is met earlier. For the isotropic case the results are shown in Figure 7(f))–7(i)). For both experiments we chose $\lambda\in\{10,25\}$, yielding eight results in total. The strongest regularization effects can be observed for $s=1,\lambda=10$, cf. Figure 7(b)) and 7(f)). As can be seen the anisotropic model introduces more smoothing within the nonlocal neighborhood than the isotropic one. Increasing $\lambda$ to $25$, cf. Figure 7(c)) and 7(g)), and hence increasing the importance of the data fidelity term reduces the smoothing but also introduces a single outlier in the middle of the slice. Furthermore, the isotropic case already seems to keep a little bit of noise in the left half of the brain. Increasing the patch size to $s=2$, see Figures 7(d)), 7(h)), 7(e)), and 7(i)), leads to results which are more smoothed within anatomical structures, e.g., the longitudinal fissure in the anisotropic case, $\lambda=10$, in Figure 7(d)) is denoised quite nicely. On the other hand larger patches also tend to oversmooth small details and fine structures within the data.

Denoising on an implicitly given surface.

Furthermore, we can use the full 3D brain data set to demonstrate the advantages of using a graph-based modeling of the data topology. Using our approach we can not only denoise manifold-valued data given on a regular grid, but also on implicitly given surfaces. Let $\mathcal{G}\coloneqq\{1,\ldots,112\}\times\{1,\ldots,112\}\times\{1,\ldots 50\}$ denote the index set of the Camino brain data set. We define the surface of the imaged human brain as those voxels that belong to the brain and additionally have at least two surrounding background pixels in their direct neighborhood. By this we avoid single missing measurements already being treated as surface voxels. We then construct a finite weighted $\varepsilon-$ball graph $G=(V,E,w)$ from of all surface voxels that have a spatial distance less than $\varepsilon$ voxels. The resulting graph consists of $\lvert V\rvert=14,384$ nodes and $\lvert E\rvert=141,978$ edges. Only $23$ isolated open space pixels were excluded to obtain a connected graph.

The resulting set of surface voxels is shown in Figure 8(a)), where one can see the front left side of the measured brain data slightly from above. We then apply the semi-implicit minimization scheme based on Jacobi iterations (14) for the anisotropic $p$-graph Laplacian model (4), $\lambda\in\{10,25,50\}$ and neighborhoods of spatial distance $\varepsilon\in\{2,3\}$ in the second and third row of Figure 8, respectively. For termination of our method we use the same stopping criterion as in the last experiment. As observed before, for $\lambda=10$, cf. Figures 8(b)) and 8(e)), the smoothing effect on the surface is strongest, leading to nearly constant areas for the latter case. For $\lambda=25$ both neighborhood sizes, cf. Figure 8(b)) and  8(e)) still smooth the surface. For the smaller neighborhood size some features already appear to look perturbed (most prominently the green ellipsoids). This is also the case for $\varepsilon=3,\lambda=50$, cf. Figure 8(g)). However, all brain surface regions seem to be more piecewise constant than for $\varepsilon=2,\lambda=25$ in Figure 8(c)). Both results for $\lambda=50$ (high data fidelity) already look quite noisy in the lower part of the data.

6 Conclusion and future work

In this paper we proposed a graph framework for processing of manifold-valued data. Using this framework we translated variational models and partial differential equations defined over manifold-valued functions to discrete domains of arbitrary topology. Modeling the relationship of data entities by finite weighted graphs allows to handle both local as well as nonlocal methods in a unified way. The nonlocal graph $p$-Laplacian provides furthermore a new approach for manifold-valued denoising that leads to results comparable to state-of-the-art methods. As the proposed framework consists of three independent parts it can be easily adapted to different manifolds, other variational models and partial differential equations, or more sophisticated numerical solvers. After the introduction of a family of graph $p$-Laplacian operators in this work, an application of our approach to other processing tasks, e.g., segmentation or inpainting of manifold-valued data, is immediate. We will investigate these applications in future work and also translate additional variational models and partial differential equations to our framework. In particular, the nonlocal Eikonal equation and the $\infty$-Laplacian operator on finite weighted graphs will be interesting to investigate for segmentation and interpolation tasks on manifold-valued data. From a theoretical point-of-view there are various challenging questions to be explored in future works. First, investigating the transition from our proposed discrete model on finite weighted graphs to continuous domains will give further insights in underlying properties and characteristics of the involved equations and their solutions. For this, one has to define the right measures to analyze a sequence of finite weighted graphs which get denser and denser for $|V|\rightarrow+\infty$. Furthermore, using elaborated concepts from differential geometry will help to prove consistency of the involved operators and energy functionals with respect to existing continuous operators, e.g., the Dirichlet energy and the graph $p$-Laplacians. An investigation of the rigorous stability conditions for the numerical algorithms proposed in this work is an interesting topic on its own. Instead of heuristically choosing the time step width $\Delta t$ in the explicit Euler time discretization, one could try to deduce sharp stability conditions. It becomes clear that this will be directly depending on the local curvature of the underlying manifold and hence involves estimations based on further tools from differential geometry, especially the curvature tensor and Christoffel symbols.

Acknowledgements.

The authors would like to thank Jan Lellmann for helpful discussions on the visualization of the NED dataset, Johannes Persch for fruitful discussions, and Gabriele Steidl for carefully reading a preliminary version of this manuscript. RB acknowledges support by the German Research Foundation (DFG) within the project BE 5888/2-1. DT acknowledges his funding by the ERC via grant EU FP7-ERC Consolidator Grant $615216$ LifeInverse.

\AtNextBibliography