Diffusion generated methods for denoising target-valued images

We consider the inverse problem of denoising an image where each point (pixel) is an element of a target set, which we refer to as a target-valued image. The target sets considered are either (i) a closed convex set of Euclidean space or (ii) a closed subset of the sphere such that the closest point mapping is defined almost everywhere. The energy for the denoising problem consists of an $L^2$-fidelity term which is regularized by the Dirichlet energy. A relaxation of this energy, based on the heat kernel, is introduced and the associated minimization problem is proven to be well-posed. We introduce a diffusion generated method which can be used to efficiently find minimizers of this energy. We prove results for the stability and convergence of the method for both types of target sets. The method is demonstrated on a variety of synthetic and test problems, with associated target sets given by the semi-positive definite matrices, the cube, spheres, the orthogonal matrices, and the real projective line.


Introduction
Let Ω ⊂ R d be a Euclidean set with smooth boundary. We consider noisy target-valued data, f : Ω → R k , that takes values in (or near) a certain target set, T ⊂ R k . We assume that either (i) T ⊂ R k is a closed convex set or (ii) T ⊂ R k is a closed subset of the unit sphere, S k−1 , such that the closest point mapping, Π T , is defined almost everywhere; see Section 3.2.
Our goal will be to find a smooth map, u : Ω → T , that approximates the data, f . For α > 0, we consider the general inverse problem, (1) inf For a vector valued function u, the gradient in E α should be interpreted component-wise. The parameter α controls the tradeoff between the 'smoothness of u' and the 'fidelity to the data'; in the limit α 0, we obtain u ≡ f on Ω. For some target sets, T , it is difficult, either analytically or computationally, to handle the constraint that u take values in T . (This is why we didn't specify the class of functions u : Ω → T to take the infimum in (1).) k T L(x) comment section k R k 0 harmonic function k T ⊂ R k convex 1 2 dist 2 (x, T ) convex set-valued field §4.3 n 2 SPD(n) 1 2 dist 2 (x, SPD(n)) SPD matrix-valued field §4.4, §4.5 1 {±1} 1 4 (x 2 − 1) 2 Allen-Cahn 2 S 1 1 4 (|x| 2 − 1) 2 Ginzburg-Landau §4.3 k S k−1 1 4 (|x| 2 − 1) 2 sphere-valued field §4.1, §4.2 n 2 O(n) One penalization approach. Suppose, for the target set T , there exists a smooth auxiliary function, L : R k → R + , such that T = L −1 (0), i.e., T is the zero-level set and set of global minimizers of the non-negative function L. In this case, since T = arg min x L(x) ⊂ R k , we can use the function L to penalize when u does not take values in T . For ε > 0, one may consider the relaxation of (1), inf Together, the first and third terms of F α,ε define a prior; the image is assumed to take values in T and be smooth. Clearly, for ε small, minimizing sequences must take values very near T . When there is no data present, i.e., α → ∞, the energy in (2) simplifies to the geometric problem where F ∞,ε (u) = 1 2 Ω |∇u(x)| 2 dx + 1 ε 2 Ω L (u(x)) dx.
Energies of this general form, as well as the language "target set", appear, e.g., in [RSK89]. It is difficult to prove general theorems about when solutions of (2) or (3) converge to solutions of (1) (if they exist!) as ε 0 for general target sets, T ; generally each target set is treated on a case-by-case basis.
In what follows, we describe in more detail a few choices of target set, T , in (1), their associated auxiliary functions in (2) and (3), along with numerous applications; a summary of various choices of T is given in Table 1. Our intent is to motivate an alternative relaxation of (1), discussed below, which can be applied to all of these choices of T .
Target-valued maps and applications in imaging, inverse problems, and geometry.
Subsets of the Euclidean sphere. For T = {±1} = S 0 ⊂ R 1 , the auxiliary function L(x) = 1 4 (x 2 −1) 2 can be used. The energy, F ∞,ε in (3), corresponds to the Allen-Cahn equation [AC79]. Modica and Mortola showed that a minimizing sequence (u ε ) of εF ∞,ε converges (along a subsequence) to χ D − χ Ω\D in L 1 for some D ⊂ Ω. Furthermore, εF ∞,ε (u ε ) → 2 . For small ε > 0, the gradient flow of F ∞,ε approximates mean curvature flow. This energy serves as a building block for a variety of pattern formation models. When a data fidelity term is added, as in (2), we obtain the Rudin-Osher-Fatemi (ROF) functional, When used for image denoising, the total variation term serves as an 'edge preserving regularizer' [ROF92]. For T = S 1 ⊂ R 2 , the auxiliary function L(x) = 1 4 (|x| 2 − 1) 2 can be used. The energy in (3) corresponds to the Ginzburg-Landau equation [BBH94] and has applications in superconductors and superfluids. In imaging, the target set T = S 1 naturally arises when one tries to recover spatially dependent phase information, such as in Interferometric Synthetic Aperture Radar (In-SAR) [RPF97;Kam06]. Here, the phase difference between an interferometric pair of SAR images, obtained from slightly different camera angles, can be used to construct very accurate elevation maps. The solution of (3) with additional boundary conditions imposed can also be used to design d = 2-dimensional cross fields, which have applications in, e.g., computer graphics and quad mesh generation [VO17]. Finally, this problem is related to the simplification of vector fields for visualization [Skr+15].
For T = S 2 ⊂ R 3 , the gradient flows of these energies are related to the heat flow of harmonic maps to S 2 [EW01]. These equations can be obtained as simplifications of the Landau-Lifschitz equation describing non-equilibrium magnetism. This is also a simplification of the energy (the one-constant approximation) appearing in the Oseen-Frank theory for liquid crystals, where the field represents the preferred direction of molecular alignment [MZ09;Bal17].
The problem in (3) with T = O(n) ⊂ {x ∈ R n×n : x F = 1} was recently studied by the authors in [OW17]. Here, it is natural to associate the auxiliary function L(x) = 1 4 x t x − I n 2 F , where · F denotes the Frobenius norm. Since O(1) ∼ = S 0 , this energy reduces to the Allen-Cahn energy when n = 1. Recalling that O(n) = SO(n) SO − (n) and SO(2) ∼ = S 1 , for n = 2, the gradient flow of this energy with initial conditions taken in SO(2) reduces to the Ginzburg-Landau gradient flow [OW17]. This energy can be considered as a model problem for crystallography, where one considers a field that takes values in SO(3) modulo the symmetry group of the crystal. This is also a model problem for the three-dimensional cross field design problem [VO17]. Finally, this problem is related to problems in rigid motion planning, where one tries to find a time-dependent trajectory, u : [t 1 , t 2 ] → T , where the function takes values in a set that describes admissible rigid motions, such as, e.g., T = SO(3) [SS00].
Convex sets. When the target set T is a convex set of R k , we can generally take the auxiliary function to be L(x) = 1 2 dist 2 (x, T ) = 1 2 min y∈T dist 2 (x, y). For RGB images, the image takes values in the cube, T = [0, 1] 3 ⊂ R 3 , as further discussed in Section 4.3.
For T = SPD(n), the set of semi-positive definite (SPD) matrices, the inverse problem in (2) appears in diffusion tensor magnetic resonance imaging (MRI) [Wan16;Len+09]. This application will be further discussed in Section 4.5.
Other target sets. For the coordinate axis, T = Σ k := {x ∈ R k : i =j x 2 i x 2 j = 0}, the minimizer of (3) with an additional norm constraint gives Dirichlet partitions of Ω in the limit as ε → 0; see [CL07;WO18]. Recently, Dirichlet partitions have been used for image segmentation and data clustering [OWO14; ZO16; OR17].
The target set, T = RP n is related to the Landau-de Gennes model, where the field describing the local orientation of a crystal is described by a Q-tensor [MZ09;Bal17]. While this theory was originally used to describe nematic liquid crystals, it has also been used to describe the orientations of RNA and carbon nanotubes. Thinking of real projective space as the quotient space obtained from the n-sphere after identifying antipodal points, fields with values in T = RP 1 are referred to as line fields, where a pair of antipodal directions is assigned to each point. This application will be further discussed in Section 4.7.
Algorithm 1: A diffusion generated method for approximating minimizers of the energy in (4).
Input: Let τ, λ > 0. Set Ω ∈ R d , the target space as T ∈ R k , the data as f ∈ L ∞ (Ω; R k ) and the initial guess as u 0 = f . Output: A function u ∈ L ∞ (Ω; T ) that approximately minimizes (4). Set s = 1. while not converged do 1. Diffusion Step. Extend u s−1 and f to R d \ Ω by zero. Solve the initial value problem for the free space diffusion equation until time τ : Step. Define u s ∈ L ∞ (Ω; T ) by point-wise applying Π T toũ, Set s = s + 1.
Results. In this paper, we derive and study an alternative relaxation of (1) than (2) based on the heat kernel. Namely, for τ > 0, λ ∈ (0, 1), and f ∈ L ∞ (Ω; R k ), we consider the relaxation of (1) given by Here, ·, · , is the L 2 (Ω; R k )-inner product and e ∆τ denotes the solution operator for the free space diffusion equation at time τ . If u is a vector-valued field, the diffusion operator is understood to be applied component-wise. Here, the parameter τ measures the relaxation of the problem and the parameter λ = τ α controls the data fidelity. The two terms in (4b) come from relaxing the two terms in the energy in (1). This is explained in Section 2, together with conditions on the target set T such that (4) is well-defined, and interpretations of (4). The second term of E λ,τ in (4b) is similar to the region-based active contour model in [Li+08], where the idea is to use a nonlocal fidelity term to characterize the image better.
The main contribution of this paper is to derive and analyze a diffusion generated method to solve (4) for a wide class of target sets, T , including those discussed in Table 1. The proposed algorithm is given in Algorithm 1. The Algorithm consists of taking a convex combination of the previous time step and the data, diffusing until time τ , and applying a map, Π T , point-wise to the resulting function. Here Π T is the convex projection onto the target set, T , in the case that T is convex and the closest point mapping otherwise. A derivation of Algorithm 1 for the target sets considered, as well as convergence properties of the algorithm are given in Section 3.
Algorithm 1 is conceptually simple, computationally efficient, easy to implement, and applicable to a broad class of problems. Algorithm 1 can be interpreted as a splitting method for (2); see Section 3.3. However, we prefer to interpret Algorithm 1 in terms of (4) since neither rely on the auxiliary function L as in (2). There are two extremes for Algorithm 1. If λ = 0, we ignore the data and find an approximate harmonic function with values in T . This is similar to the geometric problem (3) discussed above. If λ = 1, we simply dampen the highly oscillatory terms in f and apply the mapping Π T once.
Remark 1.1. The argument in [LY18,Section 4.1] shows that the output of Algorithm 1 satisfies homogeneous Neumann boundary conditions on ∂Ω. Dirichlet boundary conditions are also discussed by Laux and Yip, but we don't impose these conditions in the present work.
Previous work on diffusion generated methods. The proposed method (Algorithm 1) falls into the class of diffusion generated methods (DGMs). DGMs were first introduced for T = {±1} and showed to be associated with mean curvature flow in [MBO93;MBO92]. DGMs have also been generalized to generate high order geometric motions, such as Wilmore and surface diffusion flows, in [ERT08]. In [ERT10], the authors used the diffusion of the distance function to generate mean curvature flow where the thresholding step was replaced by redistancing. DGMs were recently shown to be stable and generalized to multiphase mean curvature flow in [EO15] and applied to wetting problems in [XWW17]. DGMs have been used for inverse problems for T = {±1} in [Wan+17 ;ET06]. The convergence rate of a DGM to a stationary point was proven in [OW17]. DGMs for T = S 1 were introduced in [Ruu+01], used for quad mesh generation in [VO17], and proven to be convergent in [LY18]. DGMs for T = S 2 were also studied in [EW01]. Finally, DGMs for T = O(n) was introduced and studied in [OW17].
Outline. In Section 2, we describe properties of (4). In Section 3, we derive and study Algorithm 1 for the two types of target sets considered. In Section 4, we use the proposed methods to study a variety of numerical experiments. We conclude in Section 5 with a discussion.

Properties and interpretation of the relaxed problem, (4)
In this section, we motivate and derive properties of the energy, E λ,τ , in (4), show the existence of solutions to (4), and give two interpretations of E λ,τ .
2.1. Motivation and properties of E λ,τ . For u 0 ∈ L 1 (Ω; R k ), we write e ∆τ u 0 to denote the solution to the free space heat equation with initial condition u 0 at time τ , Here · is interpreted as the dot product in R k or the Frobenius inner product if u and v are matrix valued fields.
Then there exists a solution u ∈ L 2 (Ω; T ) that attains the infimum in (4).
Proof. We use the direct method in the calculus of variations to establish existence. By Lemma 2.1(i), the functional is non-negative on L 2 (Ω; T ). Let (u j ) j∈N ⊂ L ∞ (Ω; T ) be a minimizing sequence, i.e., E λ,τ (u j ) → E λ,τ = inf u∈L 2 (Ω;T ) E λ,τ (u) as j → ∞. By Lemma 2.1(vi), the minimizing sequence is bounded in L 2 (Ω; T ). By, e.g., [Eva90, Thm 1.1.2], there exists a subsequence, which we continue to denote by (u j ) j∈N and u ∈ L 2 (Ω; T ), such that u j u . By the continuity of E λ,τ with respect to the weak L 2 (Ω; T ) topology (see Lemma 2.1(ii)), which shows that u attains the infimum.
2.3. Fourier interpretation of the relaxed energy. In this section, we use the Fourier transform to describe the sense in which the minimizer in (4) achieves a balance between smoothness and fidelity to the data, f . Recall that the solution to the diffusion equation with initial data, u(x, t) = u 0 (x) can be expressed using the Fourier transform, u( It follows that the energy in (4b) can be written This transformation shows that any test function u with small energy has the following two properties: • The Fourier transform of u, given byû(k), should be small for |k| τ −1/2 . • The Fourier transform of the residual, u − f , should be small for |k| τ −1/2 .

2.4.
Perimeter interpretation of the relaxed energy. In this section, for convenience, we discuss the target set T = {0, 1}, which can be obtained by shifting and rescaling T = {±1}. In [EO15], Esedoglu and Otto use k indicator functions to implicitly represent each domain and the interfaces, γ ij , between D i and D j . When τ 1, the area of γ ij can be approximated by is the Gaussian kernel; see also [AB98;Mir+07]. Up to a constant, this is equivalent to the regularity term in (5). The expression in (7) was shown to Γ-converge to the area of γ ij when τ 0 in [AB98; Mir+07; EO15]. In [EO15], based on this approximation, Esedoglu and Otto successfully generalized the original MBO method to a general threshold dynamics method to model the multiphase mean curvature flow by using a relaxation and linearization procedure. This procedure provides a proof of unconditional stability and consistency of the algorithm. In [LO16], the algorithm was rigorously proved to converge to multiphase mean curvature flow with an angle constraint at the multiple junction when τ 0. A convergence proof for T = S 1 is given in [LS17].

Derivation and properties of the diffusion generated algorithm
In the following two subsections, we separately derive a diffusion generated method when the target set is (i) a convex set or (ii) a closed subset of the unit sphere. Both derivations lead to Algorithm 1. In Section 3.3, we give an energy splitting interpretation of Algorithm 1.
3.1. The target set is a closed convex set. When the target set, T , is convex, we directly use the gradient projection algorithm (see, e.g., [Ber15]) for a fixed time step size, 1, to find the solution of (4). That is, if u n is the solution at the n-th iteration, we define the n + 1-th iteration, u n+1 , by Here, is the variation of E λ,τ (u) with respect to u at u = u n ; see Lemma 2.1(iii). Direct calculation gives Here, we take a convex combination of the data and current iterate, solve the diffusion equation until time τ , and project into the target set, T . Since e ∆τ ((1 − λ)u n + λf is a C ∞ (Ω; R k ) function, the projection step is equivalent to where Π T (v) is the point-wise convex projection of v ∈ R k into the target set T . The algorithm is summarized in Algorithm 1.
The following theorem gives a convergence result for Algorithm 1 for T convex.
Using the definition of the method in (8), we obtain By the non-expansiveness of the convex projection, it follows that Now adding the two inequalities Combining this with the inequality where L is the Lipschitz constant computed in Lemma 2.1(iv), we obtain the inequality Using this inequality in (9), we obtain On one hand, by summing (10) from n = 0 to N and letting N → ∞, we obtain We conclude that On the other hand, from (10), we also have Hence, we are led to the conclusion that the sequence { u n − u } is convergent hence bounded, implying that { u n } is also bounded. Thus, {u n } n∈N weakly converges in L 2 (Ω; T ) to aũ ∈ L 2 (Ω; T ), i.e., lim n→∞ u n −ũ, v = 0, ∀v ∈ L 2 (Ω; T ). To see thatũ = u , prove strong convergence, and obtain the convergence rate, we use the non-expansiveness of the projection to obtain The desired statement then follows.
3.2. The target set is a closed subset of the sphere. We consider a target set, T , satisfying the following properties: (1) T is a closed subset of the sphere, S k−1 , i.e., T ⊆ {x : |x| = 1} ⊂ R k .
(2) There exists a measure zero set, N ⊂ R k , such that for every point in R k \ N , we can define the closest point map, Π T : R k \ N → T , which takes points to their closest point in T , (3) We define the closed convex set B = conv(T ) to be the convex hull of T .
Example. For T = O(n) ⊂ R n×n , N is the set of singular n × n matrices, and where A = U ΣV t is the singular value decomposition of A. We have B = {A ∈ R n×n : A s ≤ 1} where · s is the spectral norm. See [OW17] for further details.
Since Ω ⊂ R d is compact, L ∞ (Ω; T ) ⊂ L 2 (Ω; T ). If T ⊂ R k is compact, then the following Lemma shows that the converse holds.
Proof. If T is compact, then there exists M > 0 such that |z| ≤ M for every z ∈ T . Let f ∈ L p (Ω; T ) for p ≥ 1. Then ess.im(f ) ⊂ T implies that for all z ∈ R k \ T , there exists ε > 0 such that In particular, we have that µ ({x ∈ Ω : |f (x)| > M }) = 0 which shows that which implies that f ∈ L ∞ (Ω; T ).
The existence to a solution in (11a) follows from Theorem 2.2. The following Lemma follows from calculations similar to as in the proof of Lemma 2.1.
The iterates in Lemma 3.5 are equivalent to those in Algorithm 1. The following theorem gives the stability of Algorithm 1. The proof can be adapted directly from [OW17,Prop. 4.3].
Remark 3.7. Regarding the condition assumed in Theorem 3.6, in practice, we do not observe that e ∆τ ((1 − λ)u n + λf )(x) ∈ N for any x ∈ Ω or n ∈ N. If this condition did occur, a random closest point could be assigned by Π T . 3.3. Energy splitting interpretation of Algorithm 1. In this section, we interpret Algorithm 1 as a splitting method for (2).
We then evolve u n by the gradient flow of E 1 until time t = τ , (16) v t = ∆v, v| t=0 = u n to obtain u n = v(τ ). Finally, we consider E 3 and set We can directly solve (15), by setting the first variation to zero, and solving for u to obtain This together with (16) gives the diffusion step in Algorithm 1 for a modified convex combination parameter.
As ε 0, the solution to (17) is obtained by taking a function u ∈ L 2 (Ω; T ) such that u = arg min which is precisely the point-wise projection, This is the projection step in Algorithm 1.

Computational examples
In this Section, we demonstrate the diffusion generated method in Algorithm 1, developed in Section 3, on several synthetic and test numerical experiments. Several of the numerical experiments considered are from [Bač+16], which can provide a comparison.
For p ∈ M and ξ ∈ T p M, let Γ p,ξ be the unique geodesic starting from Γ p,ξ (0) = p withΓ p,ξ (0) = ξ. Then we define exp p : T p M → M by exp p (ξ) = Γ p,ξ (1). Then, one spherical lemniscate curve can be obtained by γ S (t) = exp p (γ(t)) with p = (0, 0, 1). We discretize the curve in the parameter space of t at t i := 2πi 511 , i = 0, . . . , 511 to get an S 2 -valued signal, which we denote by {f o,i } 511 i=0 . We add noise to the data by taking f i = exp f o,i (η i ) with η i = (η i1 , η i2 ), where η i1 and η i2 are independent Gaussian noises with standard deviation of π 30 . In Figure 1(a), the blue markers indicate the original data and the red markers indicate the noisy data. In this example, we take the target set, T , to be S 2 = {x ∈ R 3 : |x| = 1}. Then, we have N = {x ∈ R 2 : x = 0} and Π T (x) = x |x| if x = 0. Applying Algorithm 1, we get the denoised results shown in Figures 1(b)-1(d) with a fixed τ = 10 −3 and λ = 0.05, 0.1, and 0.15, respectively. Since the original image is periodic, we solve the diffusion equation in Algorithm 1 with the periodic boundary condition.
We observe that the denoised results very closely match the original data and that the results are relatively insensitive to the parameter λ. All of these simulations were completed within 0.01 seconds. We sample both dimensions, t and s, with n = 64 points to obtain a discrete vector field f o ⊂ S 2 which is given in Figure 2(a). Similar to the noise in 4.1, we add the Gaussian noise in the tangential plane of each point with standard deviation 4π 45 ; the noisy data is displayed in Figure 2(b). In this example, we take the target set, T , to be S 2 = {x ∈ R 3 : |x| = 1}. Then, we have N = {x ∈ R 3 : x = 0} and Π T (x) = x |x| if x = 0. Since the original image is periodic, we solve the diffusion equation in Algorithm 1 with the periodic boundary condition.
The results of Algorithm 1 with τ = 10 −3 and λ = 0.05, 0.1, 0.15, and 0.2 are displayed in Figure 2(c)-(f). The numerical results show that Algorithm 1 is robust with respect to the value of the parameter λ and it is applicable to recover the original data image in Figure 2(a). It is also observed that the denoised image is slightly smoothed when decreasing the value of λ. Again, the algorithm performs very efficiently; all simulations can be done within 0.2 seconds. 4.3. Example: the 'peppers' image. Following [Bač+16, §5.1], we consider the denoising of the 'peppers' image, shown in Figure 3(a). It is distorted with Gaussian noise in each of the red, green, and blue (RGB) channels with standard deviation as 0.1, as show in Figure 3 We first consider the image represented in RGB channels. Here, the target set is the unit cube, T = [0, 1] 3 ⊂ R 3 . Since this is a convex set, we let Π T be the convex projection and apply Algorithm 1 to denoise Figure 3 We use a relatively large value of λ compared to other numerical experiments since the original image is non-smooth. We observe that the denoised images are slightly blurred, but the results are robust to changes in the parameter λ. We use the peak signal-to-noise ratio (PSNR) to evaluate the quality of the denoised image. The PSNR obtained in [Bač+16] for both RGB and HSV channels ranges between 28.16 to is 31.24, which are slightly better than the values obtained using this much simpler method. Here, we define where R x i ,x j (t) is the rotation matrix in the x i , x j -plane with angle t and δ a,b = 1 if a > b 0 else .  We discretize the parameter space (s, t) with 25×25 grid points to obtain a 25×25 matrix-valued image f = {f i,j } 25 i,j=1 ⊂ SPD(3). The SPD matrix is visualized in Figure 6(a) by the corresponding   ellipsoid at each pixel location. The noisy data in Figure 6(b) is generated by adding Rician noise with standard deviation 0.03,f = AA T . Here, A = R + √ 0.03B where R T R = f is the Cholesky factorization of f and B is a 3 × 3 upper triangular matrix with each element being a random number from the standard normal distribution. Since f is symmetric positive definite, R is well defined in the Cholesky factorization. We note that adding noise in such way implies thatf is symmetric positive definite.
In this example, we take the target set to be the group of 3 × 3 symmetric positive definite matrices, SPD(3), and solve the free space heat diffusion equation in Algorithm 1. Then, we use the mapping where U ΣV T is the singular value decomposition of A and Σ + (i, j) = max(Σ(i, j), 0).        As in the example in Section 4.4, we take the target set to be the group of 3 × 3 symmetric positive definite matrices, SPD(3), and solve the free space heat diffusion equation in Algorithm 1. Figure 9 displays the reconstructions obtained using Algorithm 1 of the data in Figure 8 with τ = 10 −4 and λ = 0.1, 0.2, and 0.3. The first row displays the reconstructed data in Ω 0 and the second row displays the subset Ω 1 of the corresponding reconstructed data. The results are very similar to those in [Bač+16]. 4.6. Example: an RP 1 -valued image. In this section, we denoise a real projective line (RP 1 )valued image. It is not clear that real projective spaces, RP n , for n > 1 can be cast within the framework of the proposed methods. However, due to the topological equivalence of RP 1 with the circle, S 1 , we can study real projection line-valued images as follows.  We recall that RP n−1 can be viewed as identifying antipodal points of the unit n-sphere, S n−1 ⊂ R k . Thus, identifying R 2 with C, we can uniquely represent each element {x, −x} ∈ RP 1 , by its squared value, x 2 ∈ S 1 ⊂ C. We then denoise this representation map using Algorithm 1 and the mapping Π T = x |x| . Finally, by taking the two square roots of the denoised representation map, we obtain a denoised RP 1 -valued image. A similar strategy is employed in [VO17] for cross-valued fields.
We define a synthetic RP 1 -valued image as follows. First define z(x, y) = ( √ ix − y), ( √ ix − y) where (·) and (·) denote the real and imaginary parts. Now we define the line field, Here, [z] + = (|z 1 |, z 2 ) ∈ S 1 ∩ {z : z 1 > 0} and [z] + = −[z] + . We discretize the parameter space (x, y) by 20 × 20 grid points to give a discretized line field, {f i,j } 20 i,j=1 ⊂ RP 1 , which is displayed in Figure 10(a). We then add Gaussian noise with standard deviation 0.3 to (x i , y i ) point-and component-wise to have (x i ,ỹ i ). Noisy data is then generated byf i,j = f (x i ,ỹ j ) and plotted in Figure 10(b).
Proceeding as described above, we view the noisy line field as taking values in C and point-wise square the values to obtain the representation field,f 2 ⊂ S 1 . We apply Algorithm 1 tof 2 with target set, T = S 1 . We then take the point-wise square root of the resulting field to obtain the denoised line field. The results are displayed in Figure 11 with τ = 10 −2 and λ = 0.05, 0.1, and 0.15. We observe that the reconstructed images are very close to the original data and the index +1/2 singularity is well-preserved.     ) and (e). We apply Algorithm 1 on the 'squared field' {f 2 ij } and display the denoised results in Figures 12(c) and 12(f). We observe that the denoised line field is a good model of the original fingerprint. In Figure 12(c), the blue dot indicates a singularity with index 1/2. In Figure 12(f), there are two singularities: the one indicated by the blue dot has index 1 and the one indicated by the green dot has index −1/2.
In both experiments, we choose τ = 10 −2 and λ = 0.15. We extracted the line field at 38 × 28 points for 12(a) and 35 × 25 points for 12(d). Both simulations were done in 5 × 10 −3 seconds, which demonstrates the efficiency of the proposed algorithm.

Discussion
In this paper, we introduced and analyzed a nonlocal energy for denoising target-valued images. We derived a diffusion generated method to minimize the energy and performed a variety of numerical experiments to show that the method is efficient, stable, and applicable to a wide variety of target-sets. There are a variety of interesting future directions for this work.  The closest comparison for our numerical results can be found in [WDS14;Bač+16]. The models developed in these papers are nonsmooth variational models which include total variation or second order differences in the regularization term. While we expect that these methods preserve edges better than the proposed method, the results in the numerical experiments are visually very similar. However, due to the simplicity of viewing the target set in an ambient Euclidean space, our methods should be faster. As with any inverse problem, the 'best' method depends on the structure of the image and the noise as well as the size of the data. More work should be done to understand the statistical framework for which these methods are consistent and robust estimators.
In this method, we have taken the domain, Ω, to be a Euclidean set. It would be very interesting to consider the case when Ω is a graph and the energy (4) is formulated using the analogous graph operators [Gen+14;BT18].
In this work, we have only looked at the denoising problem for target-valued images. Other image analysis tasks for target-valued images, including inpainting, segmentation, and registration, could be handled using similar techniques.