Restoration of Manifold-Valued Images by Half-Quadratic Minimization

The paper addresses the generalization of the half-quadratic minimization method for the restoration of images having values in a complete Riemannian manifold. We recall the half-quadratic minimization method using the notation of the c-transform and adapt the algorithm to our special variational setting. We prove the convergence of the method for Hadamard spaces. Extensive numerical examples for images with values on spheres, in the rotation group SO(3) and in the manifold of positive definite matrices demonstrate the excellent performance of the algorithm. In particular, the method with SO(3)-valued data shows promising results for the restoration of images obtained from Electron Backscattered Diffraction which are of interest in material science.


1.
Introduction. Many edge-preserving variational methods for the denoising or inpainting of real-valued images utilize the following model: let G := {1, . . . , n} × {1, . . . , m} be the image grid, ∅ = V ⊆ G and N (i) + := (i 1 + 1, i 2 ), (i 1 , i 2 + 1) the set of right and upper neighbors of pixel i ∈ G, where we impose mirror boundary conditions. From corrupted image values f : V → R we want to restore the original image u 0 : G → R as a minimizer of one of the following energy functionals where λ > 0 is a regularization parameter and ϕ : R ≥0 → R ≥0 . For V = G this is a typical denoising model in the presence of additive Gaussian noise. Otherwise, the model can be used for inpainting the missing image values in G\V. Throughout this paper, we consider even functions ϕ : R → R ≥0 . For ϕ(t) := |t|, the models (1) and (2) are discrete variants of the anisotropic and isotropic Rudin-Osher-Fatemi model [42], respectively. Then the regularization term is often referred as anisotropic/isotropic discrete total variation (TV) regularization in resemblance to its functional analytic counterpart. Remark 1. More generally one may consider G as vertices of a graph with edge set E := (i, j) : i ∈ G, j ∈ N (i) for some appropriate neighborhoods N (i). This is for example useful in nonlocal means approaches. Then the regularizing term sums over the edge set E. For simplicity we restrict our attention to the special neighborhoods N (i) + here.
Starting with the inaugural work [21,22], a huge number of papers have examined the so-called half-quadratic minimization methods for solving the above restoration problems with various functions ϕ as well as other optimization problems. We only mention the ARTUR algorithm in [16]. Basically, the original problem is reformulated into an augmented one which is quadratic with respect to the image and separable with respect to additional auxiliary variables. Then an alternating minimization process is applied whose steps allow an efficient computation. Half-quadratic minimization is connected with other well-known minimization approaches. We only mention the relation to EM algorithms [14], quasi-Newton minimization [4,33] and gradient linearization algorithms [33]. A gradient linearization method was used in particular in [53,54] to minimize an approximate total variation regularization (2) with ϕ(t) := √ t 2 + ε 2 for ε 1. It was called the "lagged diffusivity fixed point iteration" and the authors mention that it amounts to apply the multiplicative form of half-quadratic minimization to this ϕ. Finally, there is a relation to iteratively reweighted least squares methods [18,30]. For a convergence analysis of half-quadratic minimization methods for convex functions ϕ we refer to [16,34] and for weaker convergence results for nonconvex ϕ to [19].
In many applications signals or images having values in a manifold are of interest. Circle-valued images appear in interferometric synthetic aperture radar [13,20] and various applications involving the phase of Fourier transformed data. Images with values in S 2 play a role when dealing with 3D directional information [27,29,51] or in the processing of color images in the chromaticity-brightness (CB) setting [15]. The motion group and the rotation group SO(3) were considered in tracking, (scene) motion analysis [38,40,50] and in the analysis of back scatter diffraction data [8]. Finally, images with values in the manifold of positive definite matrices appear in DT-MRI [36,45,55,57] and whenever covariance matrices are adjusted to image pixels, see, e.g., [50].
Recently a TV-like model for circle-valued images was introduced in [46,47]. For manifold-valued image restoration such an approach was proposed in [31], where the problem was reformulated as a multilabel optimization problem which was handled using convex relaxation techniques. Another method suggested in [56] employs cyclic and parallel proximal point algorithms and does not require labeling and relaxation techniques. This approach was generalized by including second order differences for circle-valued images in [9,10], for coupled circle and real-valued images in [11] and for Riemannian manifolds in [5]. A restoration method which circumvents the direct work with manifold-valued data by embedding the matrix manifold in the appropriate Euclidean space and applying a back projection to the manifold was suggested in [41]. Recently, an iteratively reweighted least squares method for the restoration of manifold-valued images was suggested in [24]. This method can be seen as multiplicative half-quadratic minimization method for the special function ϕ(t) := √ t 2 + ε 2 . In this paper we adopt the idea of half-quadratic minimization for general functions ϕ for the restoration of manifold-valued images. We prefer the notation of the c-transform known from optimal transport to recall the basic half-quadratic minimization approach. Then we describe the algorithm for our problems of denoising or inpainting of manifold-valued images both in the anisotropic and isotropic case. Here we focus on the multiplicative half-quadratic minimization method. Convergence of the algorithm can be shown for images having entries in an Hadamard space. The manifold of positive definite matrices is such an Hadamard manifold. We provide several applications of the algorithm as denoising of phase-valued images, restoration of color images with disturbed chromaticity or of 3D directions, and improvement of images obtained from electron backscatter diffraction of a Magnesium sample.
The outline of the paper is as follows: In Section 2 we propose our variational model and show how to handle it by the half-quadratic minimization approach. A convergence proof for the algorithm applied on Hadamard manifolds is given in Section 3. Section 4 shows various numerical examples. Finally, Appendix A contains the proofs and Appendix B lists special quantities necessary for the numerical computations.
2. Half-quadratic minimization. Let M be a complete, connected n-dimensional Riemannian manifold with geodesic distance d : M × M → R ≥0 . Now we consider manifold-valued images. More precisely, from corrupted image values f : V → M we want to restore the original manifold-valued image u 0 : G → M as a minimizer of one of the following energy functionals with λ > 0 being a regularization parameter and ϕ : R ≥0 → R ≥0 . As before we set ϕ(t) := ϕ(−t) for t < 0 and consider ϕ as a function defined on the whole real axis. For ϕ(t) := |t|, the second sum in the regularizer of J ν is just d(u i , u j ) j∈N (i) + ν , ν ∈ {1, 2}. Then J 1 resembles the setting in [56] and is related to the anisotropic ROF functional (1) and J 2 yields the approach [31] related to the isotropic case (2).
In this paper we will consider smooth regularization terms, i.e., even, differentiable functions ϕ. We will compute a minimizer of J ν , ν ∈ {1, 2}, by half-quadratic minimization methods.
In the following, we briefly recall the reformulation idea of half-quadratic minimization using the concept of the c-transform and apply it to (3) and (4). The c-transform of functions defined on metric spaces is used in connection with optimal transport problems see, e.g., [52, p. 86f] and seems also to be an appropriate approach here. Given a function c : R × R → R, the c-transform of a function ϕ : R → R is defined by By this definition we see that Recall that h * * = h if and only if h is lower semi-continuous (lsc) and convex. In the multiplicative and additive half-quadratic methods we use the functions respectively. The multiplicative setting was introduced in [21] and the additive one in [22]. The quadratic cost function c(t, s) in the additive setting is also handled in optimal transport topics. Note that the cost function c(t, s) in the multiplicative case is not bounded from below in s. The following proposition is crucial for the half-quadratic reformulation. Proposition 1. Let ϕ : R → R ≥0 be an even, differentiable function and let c : R × R → R be given by (5) and (6), respectively.
Note that in the multiplicative case ψ(s) = −∞ for s < 0, so that we can restrict our attention in the infimum in (8) to s ≥ 0. Further, the assumption ϕ (t) ≥ 0, t ≥ 0, is in particular fulfilled if ϕ(t) is convex and ϕ (0+) ≥ 0. The proof can be given following for example the lines in [16,34]. However, since the assumptions in these papers are slightly different and the c-transform notation is not used, we add the proof in the appendix A to make the paper self-contained. The functions ϕ which fulfill the conditions in Proposition 1 and which were used in our numerical test are listed in Table 1. Further examples are collected in [16,33,34].
In the following we assume that ϕ fulfills the assumptions of Proposition 1. Now the idea is to replace ϕ in (3), resp., (4) by the expression in (8) and to consider where we have used the notation v := (v i,j ) i,j∈G in the anisotropic case and v := (v i ) i∈G in the isotropic case. By Proposition 1 i), minimizing J ν over u and v gives the same solutions for u as just minimizing J ν over u. More precisely we give the following remark.
Ifû is a minimizer of J ν , ν = 1, 2, then û, s(dû) is a minimizer of J ν , ν = 1, 2, and conversely. In particular, ifû = arg min u J ν u, s(û) holds true, thenû is a minimizer of J ν . Now we can apply an alternating minimization over v i,j ∈ R, resp., v i ∈ R and u ∈ M n×m =: M: Clearly, under the assumptions of Proposition 1 we have We want to work with the differentiable function d 2 in the second iteration (15).
Since the additive reformulation leads to a non-differentiable function d 2 + sd, we restrict our attention in the following to the multiplicative case. Minimization with respect to v. The minimization over v in (14) can be done separately for the v i,j or v i . By Proposition 1 a minimizer is given by where s is defined as in (11). Note that only for d u = 0 in the isotropic case a larger value also could be taken as a minimizer.
If ϕ fulfills the assumptions of Proposition 1 iii), then v (k+1) ∈ 0, ϕ (0+)/2 . Minimization with respect to u. The minimization over u in (15) is equivalent to finding the minimizer of respectively. We can apply, e.g., a gradient descent or a Riemann-Newton method, see [1]. Both methods are described in the following for our setting and were implemented. We need the following notation. Let T x M denote the tangential space of M at x ∈ M and ·, · x : T x M × T x M → R the Riemannian metric with induced norm · x . Let γ x,ξ (t), x ∈ M, ξ ∈ T x M be the minimal geodesic starting from γ x,ξ (0) = x withγ x,ξ (0) = ξ. Then the exponential map exp x : T x M → M is given by exp x ξ = γ x,ξ (1). The inverse exponential map denoted by log x = exp −1 x : M → T x M is locally well-defined. For the manifolds used in our numerical examples, namely the sphere S n , n ∈ N, the SO(3) and the manifold P(r), r ∈ N, of symmetric positive definite r × r matrices, the specific maps are given in the Appendix B.
be the Riemannian gradient and the Hessian of F : Considering the image u ∈ M, we abbreviate the gradient and the Hessian of a function F : M → R defined on the product manifold M also by grad F and Hess F , resp., since its use becomes clear from the context. Gradient descent method. The gradient descent method computes, starting with u (0) := u (k) , iteratively (18), resp. (19) or by Newton's approach (20); until stopping criterion is reached; with appropriate step sizes t r > 0. The gradient grad J ν,v (k+1) is given by

Algorithm 1 Image Restoration by Half-Quadratic Minimization (multiplicative)
Riemann-Newton method. Alternatively we can use a Riemann-Newton method to compute a minimizer. Finding a descent direction η r ∈ Tũ(r) M with Newton's method is done for ν ∈ {1, 2} by solving the system of equations Then we update, starting withũ (0) := u (k) , iterativelỹ The whole half-quadratic minimization method for our problem is given in Algorithm 1.
In [33] it was shown that the multiplicative half-quadratic minimization is equivalent to the quasi-Newton descent method, and therefore we expect that its performance is better than the simple gradient descent method. Let us comment on this for the manifold-valued setting.
Remark 3. (Relation of half-quadratic minimization to gradient descent and (quasi) Newton methods) We restrict our attention to the case ν = 1. Similar considerations can be done for ν = 2. The gradient of the initial functional J 1 in (3) is given (11) can be rewritten as Hence a gradient descent algorithm applied to the initial functional J 1 coincides with the half-quadratic method if we perform only one step in the gradient descent method (17) to obtain an update of u. More than one gradient descent step in (17) leads to a linearized gradient descent of J 1 . If we perform a Riemann-Newton step (20) to update u we have a quasi-Newton method for is fixed in the half-quadratic update step of u which is not the case in (21). This simplifies the computation of the Hessian in the half-quadratic approach.
3. Convergence in hadamard manifolds. We start with a general remark.
if ν = 2 (isotropic case). By construction we have for the iterates produced by Algorithm 1 that so that the sequence J ν u (k) , v (k) k∈N is monotonically decreasing. By (8) and since ϕ is nonnegative, the function J ν is bounded from below by zero and the sequence J ν u (k) , v (k) k converges to some b ν . This holds also true for J ν u (k) by (16). If M is compact as in the case of spheres or SO(3), then u (k) k is clearly bounded. If M is an Hadamard space as defined in the next subsection and ϕ is coercive, then u (k) k is bounded since J ν is by Proposition 2 coercive. In these cases u (k) , v (k) k is also bounded and therefore there exists a subsequence For Hadamard spaces many results on the convergence of algorithms carry directly over from the Hilbert space setting. This is in particular true for the halfquadratic minimization algorithm. In this section we summarize these results for convex functions ϕ and data in Hadamard spaces.
We start by recalling some basic facts. A curve γ : and strictly convex if we have a strict inequality for all 0 < t < 1. An Hadamard space is a complete metric space (H, d) with the property that any two points x, y are connected by a geodesic and the following condition holds true (23) d for any x, y, v, w ∈ X. Inequality (23) implies that Hadamard spaces have nonpositive curvature [3,39] and Hadamard spaces are thus a natural generalization of complete simply connected Riemannian manifolds of nonpositive sectional curvature, the so-called Hadamard manifolds. For more details, the reader is referred to [6,26]. Unfortunately, the spheres and the rotation group are not Hadamard manifolds, while the symmetric positive definite matrices have this nice property. The following facts can be shown similarly as in R d , see [ Then we obtain the following proposition whose simple proof is added for convenience in Appendix A. Under the assumptions of Proposition 1 iii), we have in our algorithm that v (k) > 0. Then, we see similarly as in the proof of Proposition 2 that the functionals J ν,v (k) , ν = 1, 2, are coercive and strictly convex. Thus the minimizer u (k) exists and is unique.
In the case V = G we further assume that ϕ is strictly convex. Then the sequence u (k) k∈N generated by Algorithm 1 converges to the minimizer of J ν , ν = 1, 2. The proof which follows standard arguments is given in the Appendix A. Note that the assumptions of the theorem are fulfilled for the first two functions in Table 1.
Remark 5. The assumptions in the theorem include those in Proposition 1. Note that the convexity of ϕ and ϕ (0+) > 0 imply that ϕ (t) > 0 for all t > 0. Additionally the continuity of ϕ is required to make the function s in (11) continuous and the (strict) convexity of ϕ to make the objective function strictly convex. Since ϕ is convex, its derivative is increasing. Together with ϕ (0+) > 0 this implies that ϕ (t) > 0 for all t > 0 so that ϕ is increasing for t > 0 and coercive.

Numerical examples.
In this section we demonstrate the performance of Algorithm 1 for the functions ϕ from Table 1. These functions are known for their edge-preserving properties. Note that ϕ 1 was used in the "lagged diffusivity fixed point iteration" [53] for real-valued images and in the iteratively re-weighted least squares method [24] for S 2 -valued and P(3)-valued images. The function ϕ 2 is a Moreau envelope of the absolute value function, also known as the Huber function.
Both functions are convex. The non-convex function ϕ 3 was used for edge-preserving restoration of real-valued images in [16,33]. Unless stated otherwise we use the anisotropic approach and Newton's method in our implementations. Although neither the spheres nor the rotation group are 4.1. S 1 -valued data. We start with the one-dimensional signal in Fig. 1 to show how the different functions ϕ from Table 1 perform and how the parameter ε influences the results. The original signal in Fig. 1 (a) was obtained from f (x) = 8πx 2 by sampling with size 0.01 and unwrapping modulo 2π such that the data are represented in [−π, π). Then wrapped Gaussian noise with standard deviation σ = 0.3 was added. Using ϕ 1 to restore the signal gives relatively largest error in Fig. 1 (b). The Huber function ϕ 2 and the exponential function ϕ 3 show better results in Fig. 1 (c) and 1 (d). The regularization parameter λ was adapted to get the best error where N = 101. Making ε in the Huber function larger leads to the smoother result in Fig.1 (e) which approximates the original signal only well at the beginning of the signal. In Fig. 1 (f) we choose a larger ε in the function ϕ 3 with the effect that edges of smaller height are smoothed and we have a staircasing effect for nearly equally high ascents.
Next we want to demonstrate the difference between the anisotropic (12) and isotropic (13) half-quadratic minimization methods. To this end, the function atan2(x, y) was sampled over a regular grid − 1 2 , 1 2 2 with grid size 1 128 , resulting in Fig. 2 (a). Then we corrupt the image by removing a circular region from the  center as shown in Fig. 2 (d). Using the anisotropic functional leads to Fig. 2 (b), where we observe artifacts in vertical and horizontal directions. The image produced by applying the isotropic functional in Fig. 2 (c) does not have this problem. This effect is also illustrated by the error plots in Figs. 2 (e) and 2 (f).

S 2 -valued data.
In our first example we denoise color images in the chromaticity and brightness space. For an RGB image the brightness is given by the real positive numbers b := R 2 + G 2 + B 2 1 2 and the chromaticity by the S 2 -values c := (R, G, B)/b. We want to mention that the first two examples consider problems which have values on the positive octant; only in the last example we look at data covering the whole sphere. We compare half-quadratic minimization with the different functions ϕ and the TV approach from [56] in Fig. 3. We took the image "Peppers" 1 in Fig. 3 (a) and added Gaussian noise with standard deviation σ = 0.1 to all three color channels in the RGB model. For denoising in the chromaticitybrightness model we optimized λ with respect to best PSNR for both channels separately using a grid search on 1 100 N for ϕ 1 and ϕ 3 , and for ϕ 2 on N. Furthermore, we optimized ε first in order of magnitude k = 10 j and refined the search on    more to staircasing. Using the exponential function ϕ 3 yields a worse PSNR and hence does not compete with the previous two functions; the edges are preserved but within the more constant regions some noise is left. The bright spot (see upper magnification) appears also too smooth. This originates from the too smooth transitions in the brightness which are not detected as edges. The TV regularization introduces staircasing and it is not able to reduce the noise in the dark area (lower magnification).
In the second example we use half-quadratic minimization for colorization in the chromaticity-brightness space. We assume that the brightness of the image is known, but 99 percent of the chromaticity information is lost. The original image is shown in Fig. 4 (a) and its corrupted version in Fig. 4 (b). For inpainting the chromaticity we have used a nearest neighbor initialization. With the regularizing function ϕ 1 we obtain the result depicted in Fig. 4 (c). We compare this with Fig. 4 (d) which is obtained by using the chromaticity colorization method in [37] which we have implemented for comparison.
Our final experiment shows the smoothing of 3D directions in the synthetic image in Fig. 5 (b). We use half-quadratic minimization with ϕ 2 to obtain Fig. 5 (c). The original pattern is again visible.

P(3)-valued data.
Our first example illustrates the inpainting capabilities of the half-quadratic minimization method by an artificial example. The P(3)-valued image of size 16×16 in Fig. 6 (a) has a jump at 2 3 in x direction. We destroy a center square of size 12 × 12, see Fig. 6 (b). We are able to reconstruct the inpainting area nearly perfectly by using ϕ 1 and ε = 10 −3 , see Fig. 6 (c). Nevertheless decreasing   ε introduces more and more staircasing, cf. Fig. 6 (d). This resembles the TV case shown in Fig. 6 (e), using the model from [56]. An important application of P(3)-valued image denoising is Diffusion Tensor Magnetic Resonance Imaging (DT-MRI). The Camino project 2 [17] provides a DT-MRI dataset of the human head which is freely available. 3 The complete data is given as a 3D imagef = f i,j,k ∈ P(3) 112×112×50 , where we apply the halfquadratic minimization to each of the traversal planes k ∈ {1, . . . , 50}. The original   dataset, cf. Fig. 7 (a), is plotted using the anisotropy index relative to the Riemannian distance [32] normalized onto [0, 1) and colored in hue. The half-quadratic minimization is used with ϕ 1 and the parameters λ = 0.1, ε = 10 −3 and a maximum change between two successive iterations being 10 −12 as a stopping criterion. We obtain the result shown in Fig. 7 (b). For the complete dataset of 168,169 nonzero matrices, the algorithm needed 2,492 seconds to compute the result.

SO(3)-valued data.
Processing images with SO(3)-valued entries is fundamental in the analysis of polycrystalline materials by means of Electron Backscattered Diffraction (EBSD), cf. [2,28]. Since the microscopic grain structure affects  macroscopic attributes of materials such as ductility, electrical and lifetime properties, there is a growing interest in the grain structure of crystalline materials such as metals and minerals. EBSD provides us for each position on the surface of a specimen with a so called Kikuchi pattern, which allows the identification of the structure (material index) and the orientation of the crystal at this position relative to a fixed coordinate system (SO(3) value). Since the atomic structure of a crystal is invariant under its specific symmetry group S ⊂ SO(3) the orientation is only given as an equivalence class [m 0 ] = {m 0 s | s ∈ S} ∈ SO(3)/S, m 0 ∈ SO(3). Fig. 8 (a) displays a typical EBSD image consisting of lattice orientations of deformed Magnesium collected by [43]. Each pixel of the image corresponds to a position on the surface of a Magnesium specimen. The color of the pixels is chosen corresponding to the orientation measured at this position according to the following color mapping: for a fixed vector r ∈ S 2 we consider the mapping Φ : SO(3)/S → S 2 /S, [m] → [m −1 r]. Next we colorize the quotient S 2 /S as it is depicted in Fig. 8 (b). From the colorization scheme the symmetry group S ⊂ SO(3) of Magnesium becomes visible which has six rotations with respect to a 6-folded axis (kπ/3, k = 1, . . . , 6 rotations around c direction) and six rotations with respect to 2-folded axis a 1 , a 2 perpendicular to that.
EBSD images usually consist of regions with similar orientations called grains. For certain macroscopic properties the pattern of orientations within single grains is of importance [7], e.g., for the computation of geometrically necessary dislocations [35,49] the gradient of the rotations within single grains has to be determined. As the rotation determination by Kikuchi patterns is sometimes fragile, the rotation valued images determined by EBSD are often corrupted by noise and suffer from missing data so that denoising and inpainting techniques have to be applied [25]. For detecting grains in the raw EBSD data we applied a thresholding algorithm [8]. Fig. 9 (a) displays a single grain with its rotations. Since the rotations vary very little within a single grain we applied a sharper colorization, cf. Fig. 8 (d) to make the noise and the rotation gradient visible.
We want to apply half-quadratic minimization to denoise EBSD images. Since crystallographic symmetry groups are finite the quotient SO(3)/S is locally isomorphic to SO(3). In particular, the formulas given in Appendix B.3 can be applied. In a first experiment we apply half-quadratic minimization using ϕ 1 to the rotationvalued image depicted in Fig. 9 (a) which leads to the smooth image in Fig. 9 (b).  In a second experiment we randomly removed 30% of the data shown in Fig. 9 (c). Using half-quadratic minimization for jointly inpainting and denoising the image we obtain the result shown in Fig. 9 (d) which looks very similar to those in Fig. 9 (b). In a third experiment we applied half-quadratic minimization simultaneously to several grains. The challenge from the mathematical point of view is that SO (3) is not an Hadamard manifold. Convergence can be guaranteed only locally, which is the case for single grains but may be not true when considering several grains simultaneously. From the practical point of view this case is especially interesting as missing data usually occur at grain boundaries, i.e., between grains. However, for our data we have got promising results. Fig. 10 (a) shows the grain from the previous example (pink color) and two other grains in its neighborhood. Note that the top middle area (light green) and top right area (brown-green) belong to the same grain. Pixels with missing data are plotted white. Half-quadratic minimization restoration with ϕ 1 improves the image as can be seen in Fig. 10 (b). For our last experiment we again randomly remove 30% of the data, cf. Fig 10 (c). Restoration with ϕ 1 leads to the result in Fig. 10 (d), which is again hardly to distinguish from Fig. 10 (b). Using ϕ 3 leads to even better results as depicted in Fig. 10 (e). We can adjust the smoothing in such a way that the edge distinguishing two grains is not smoothed, while smaller rotation changes are smoothed. This is advantageous in the large top grain. Finally Fig. 10 (f) shows the zoom to the (pink) grain from Fig. 9 with adapted color map.

5.
Conclusions. We adapted the principles of half-quadratic minimization to the setting of complete, connected Riemannian manifolds. In particular, the notation of the c-transform provides an interesting point of view. For Hadamard manifolds we proved the existence and uniqueness of the minimizer of the corresponding functionals as well as the convergence for the alternating minimization algorithm under moderate assumptions. The multiplicative half-quadratic minimization method resembles a quasi-Newton method [33] and appears to be very efficient in our numerical examples. There are numerous applications of the approach. In this paper,  images having values in a manifold such as the 2-sphere or the symmetric positive definite matrices were denoised. The method was also used for inpainting missing information into images consisting of either rotation matrices or symmetric positive definite matrices. In the chromaticity-brightness color model, the inpainting technique was applied to the task of colorization. The method has further potential in EBSD.
Topics of future research are the derivation of convergence proofs for more general manifolds under special assumptions on the local behavior of the data. Furthermore, different data terms must be included for other applications, and the inclusion of higher order differences into the regularization term of the model is of interest. (c) Grains with 30% lost data.
(f) Grain from Fig. 9. By assumption on Φ we know that Φ = Φ * * which implies is continuous, convex and by (10) coercive so that the global minimizer in (7) is attained for 0 = h (t) = at − ϕ (t) − s, i.e., for t, at − ϕ (t) which proves ii). 2. In the multiplicative case we obtain, since ϕ is even, We have so that we can restrict our attention to t ≥ 0. By assumption, Φ is convex and lsc. Thus, Φ = Φ * * and we obtain for t ≥ 0 that This yields i). To see ii) we first note that condition (9) implies for s ≥ 0 that the objective function in (7) is coercive such that the infimum is attained. For s < 0, we have ϕ c (s) = −∞. For s ≥ 0, we obtain arg min t≥0 t 2 s − ϕ(t) = arg min r≥0 rs − ϕ( √ r) 1 2 By assumption on Φ, the function h(r) solution is positive or for r = 0 if lim r→0+ h (r) ≥ 0, i.e., s ≥ 1 2 ϕ (0+). Therefore (t, s) = 0, 1 2 ϕ (0+) is a solution. Finally, the concavity of ϕ( √ t) for t ≥ 0 implies that ϕ ( √ t)/(2 √ t) and thus ϕ(t)/(2t) is decreasing. Under the additional assumption in iii) we get s ∈ 0, Then F (u) := 1 2 i∈V d 2 (u i , f i ) goes to infinity and the functionals J ν , ν = 1, 2, are coercive. In the case V = G choose i 0 ∈ V and let u 0 be the constant image with entries f i0 . Let d(u, u 0 ) 2 → ∞. Assume that J ν (u) remains finite, so that in particular d(f i0 , u i0 ) and d(u i , u j ), j ∈ N (i) + , i ∈ G are finite. By the construction of the neighborhoods N (i) + there exists for every j ∈ G a path i 0 , i 1 , . . . , i kj = j with i l+1 ∈ N (i l ) + and Since the right-hand side remains finite this contradicts d(u, u 0 ) 2 → ∞. Hence J ν ν = 1, 2 are coercive. ii) By (D2) we have that F (u) is convex and strictly convex if V = G. If a function h : H κ → R is convex, then, for any geodesic γ : [0, 1] → H κ joining x, y ∈ H κ , we obtain since ϕ is increasing and convex that (24) ϕ so that ϕ•h is convex. If ϕ is strictly convex the last inequality is strong so that ϕ•h is strictly convex. With h := h ij = d(u i , u j ) : H 2 → R this implies by (D1) that J 1 is convex, resp. strictly convex. Concerning J 2 notice that the convexity of h i (x 0 , x i ), i = 1, . . . , κ − 1, on H 2 implies convexity of h(x 0 , . . . , h κ−1 ) := and by the Schwarz inequality In the case V = G, strict convexity follows by the strict convexity of the data term i∈G d 2 (f i , u i ) and for strictly convex ϕ by the strict convexity in (24).
Proof of Theorem 3.2. By Remark 4 we know that lim k→∞ J ν u (k) , v (k) =:b and that there exists a subsequence u (kj ) , v (kj ) j which converges to some (ū,v).
Assume that lim k→∞ u (k) →ū. Then there exists ε > 0 such that infinitely many u (ki) not contained in the open ball B ε (ū) of radius ε centered atū. Since u (ki) is bounded, there exists a convergent subsequence u (ki j ) which converges to some point u * =ū in the closed set M\B ε . Then lim j→∞ J ν u (ki j ) = J ν (u * ) = J(ū) which contradicts the fact that J ν has a unique minimizer.
Let P(r) be the manifold of symmetric positive definite r × r matrices. It has the dimension dim P(r) = n = r(r+1)

2
. The tangent space of P(r) at x ∈ P(r) is given by T x P(r) = {x} × Sym(r) := x 1 2 ηx 1 2 : η ∈ Sym(r) , in particular T I P(r) = Sym(r), where I denotes the r × r identity matrix. The Riemannian metric on T x P reads where ·, · denotes the matrix inner product (25). The geodesic distance is given by d P (x 1 , x 2 ) := Log x For η ∈ T x SO(3) the exponential map at x ∈ SO(3) and (locally) its inverse are defined as exp x (η) := x Exp x T η , resp., log x1 x 2 = x 1 Log x T 1 x 2 . The SO(3) can be parametrized in various ways. Due to the form of our data we prefer to use quaternions for the representation and the similarity of SO(3) to the group S 3 , see, e.g., [23]: We decompose the unit quaternions q = s, v T T ∈ S 3 into a real part s ∈ R and a vector part v ∈ R 3 . The multiplication of two quaternions q 1 , q 2 ∈ S 3 is given by * (q 1 , q 2 ) = 2 arccos | q 1 , q 2 |. The exponential map of η ∈ T p S 3 * is given by exp q η := sgn s s v , (s, v) T := q cos η 2 + η η 2 sin η 2 and the logarithmic map at q 1 of q 2 by log q1 q 2 := q 2 − q 2 , q 1 q 1 q 2 − q 2 , q 1 q 1 2 arccos | q 1 , q 2 | sgn q 1 , q 2 .