extfieldsofstudy:
Mathematics
Computer Science
Medicine

s2fieldsofstudy:
Computer Science
Mathematics

text:
Inverse Diffusion Curves using Shape Optimization

The inverse diffusion curve problem focuses on automatic creation of diffusion curve images that resemble user provided color fields. This problem is challenging since the 1D curves have a nonlinear and global impact on resulting color fields via a partial differential equation (PDE). We introduce a new approach complementary to previous methods by optimizing curve geometry. In particular, we propose a novel iterative algorithm based on the theory of shape derivatives. The resulting diffusion curves are clean and well-shaped, and the final image closely approximates the input. Our method provides a user-controlled parameter to regularize curve complexity, and generalizes to handle input color fields represented in a variety of formats.

Introduction
Vector graphic images remain invaluable for a broad range of 2D applications because of their resolution independence, compactness of representation, and powerful editability. Recently, diffusion curve images  further improve the expressiveness of vector graphics by providing flexible and easyto-manipulate smooth gradients, and since then inspired a variety of novel applications [Jeschke et al. 2009;Takayama et al. 2010;Sun et al. 2012;Ilbery et al. 2013;Sun et al. 2014]. Defined along the curves, colors are diffused across the image by a Poisson or Laplace reconstruction, and their smoothness can be further controlled by the curve's blurriness through post-processing. Although efficient rendering of diffusion curve images has been well explored, its inverse problem of creating diffusion curves automatically given desired target images remains challenging.
The inverse diffusion curve problem is difficult because even though the curves themselves are 1D, their impact on the final image is nonlinear and global over a 2D domain through a partial differential equation (PDE), Laplace's equation. The geometry of curves largely determines the reconstruction quality. Previous methods have used local heuristics to obtain curve geometry. They place curves at locations indicated by edge detectors applied to the target image ] and its Laplacian or bi-Laplacian [Xie et al. 2014]. While these heuristics work well around sharp edges such as object boundaries, they have difficulty handling color variations in regions that are smooth yet visually rich.
We introduce a new approach to the inverse diffusion curve problem. Complementary to existing methods, our approach solves for curve geometry through a global optimization that takes into account the curves' full impact on the PDE-based color field. To achieve this, we characterize how modifications to a diffusion curve can reduce a global cost function determined by the solution of a PDE with the curves acting as boundary conditions.
Our method is grounded on the theory of shape optimization [Sokolowski and Zolésio 1992]. Given a color field, it computes curve geometry by minimizing a measure of the color reconstruction residual. Starting from an initial set of curves, it iteratively evolves their shapes toward an optimal configuration. Mathematically, the curves are treated as continuous functionals. This allows them to deform arbitrarily, enabling full exploration of possible curve configurations. This iterative process is similar to a surface normal flow: our curve evolution at every iteration is guided by the Fréchet derivative of the residual function with respect to the curve's boundary velocity, leading to an efficient "gradient descent" of the residual. This method is mathematically clean and easy to implement: all computations at each iteration boil down to solving a Laplace equation and a Poisson equation.
Building on our curve placement algorithm, we introduce a complete pipeline for solving the inverse diffusion curve problem. Our method generates curves in a clean and concise way, and the resulting images can accurately capture complex color variations of input color fields (see Figures 1, 14, 15, and 16 as well as the supplementary images).
We demonstrate that our method promises practical applications beyond pixel image vectorization. For instance, it enables automatic rendering of vector graphic images from 3D geometries, analogous to the traditional pixel image rendering. Further, using our algorithm, one can directly transform other formats of vector graphics, such as gradient meshes, into diffusion curve images without rasterizing the input.

Figure 1:
We present a new approach to automatically build diffusion curve images approximating provided color fields (a). Starting from a set of boundary curves (black strokes) indicating color jump discontinuities (b), our method iteratively adds curves (dark yellow strokes) and refines their shapes in an optimized manner (c, d). The resulting image (e) accurately matches the input.

Introduction
Vector graphic images remain invaluable for a broad range of 2D applications because of their resolution independence, compactness of representation, and powerful editability. Recently, diffusion curve images  further improve the expressiveness of vector graphics by providing flexible and easyto-manipulate smooth gradients, and since then inspired a variety of novel applications [Jeschke et al. 2009;Takayama et al. 2010;Sun et al. 2012;Ilbery et al. 2013;Sun et al. 2014]. Defined along the curves, colors are diffused across the image by a Poisson or Laplace reconstruction, and their smoothness can be further controlled by the curve's blurriness through post-processing. Although efficient rendering of diffusion curve images has been well explored, its inverse problem of creating diffusion curves automatically given desired target images remains challenging.
The inverse diffusion curve problem is difficult because even though the curves themselves are 1D, their impact on the final image is nonlinear and global over a 2D domain through a partial differential equation (PDE), Laplace's equation. The geometry of curves largely determines the reconstruction quality. Previous methods have used local heuristics to obtain curve geometry. They place curves at locations indicated by edge detectors applied to the target image ] and its Laplacian or bi-Laplacian [Xie et al. 2014]. While these heuristics work well around sharp edges such as object boundaries, they have difficulty handling color variations in regions that are smooth yet visually rich.
We introduce a new approach to the inverse diffusion curve problem. Complementary to existing methods, our approach solves for curve geometry through a global optimization that takes into account the curves' full impact on the PDE-based color field. To achieve this, we characterize how modifications to a diffusion curve can reduce a global cost function determined by the solution of a PDE with the curves acting as boundary conditions.
Our method is grounded on the theory of shape optimization [Sokolowski and Zolésio 1992]. Given a color field, it computes curve geometry by minimizing a measure of the color reconstruction residual. Starting from an initial set of curves, it iteratively evolves their shapes toward an optimal configuration. Mathematically, the curves are treated as continuous functionals. This allows them to deform arbitrarily, enabling full exploration of possible curve configurations. This iterative process is similar to a surface normal flow: our curve evolution at every iteration is guided by the Fréchet derivative of the residual function with respect to the curve's boundary velocity, leading to an efficient "gradient descent" of the residual. This method is mathematically clean and easy to implement: all computations at each iteration boil down to solving a Laplace equation and a Poisson equation.
Building on our curve placement algorithm, we introduce a complete pipeline for solving the inverse diffusion curve problem. Our method generates curves in a clean and concise way, and the resulting images can accurately capture complex color variations of input color fields (see Figures 1, 14, 15, and 16 as well as the supplementary images).
We demonstrate that our method promises practical applications beyond pixel image vectorization. For instance, it enables automatic rendering of vector graphic images from 3D geometries, analogous to the traditional pixel image rendering. Further, using our algorithm, one can directly transform other formats of vector graphics, such as gradient meshes, into diffusion curve images without rasterizing the input.

Related Work
Diffusion curves ] represent a color field by diffusing the colors defined along control curves over the entire image plane. The diffusion process is described by Laplace's equa-tion solved using a finite volume method. Later, solving Laplace's equation was improved using a multigrid method [Jeschke et al. 2009], triangle mesh interpolation [Pang et al. 2012], Boundary Element method [Sun et al. 2012], 2D ray tracing [Prévost et al. 2014], and Fast Multipole method ].
Inverse diffusion curve problem. Our work focuses on the inverse problem of diffusion curves. Previously, Orzan et al. [2008] proposed to place diffusion curves along edges extracted from input images using the Canny detector [Canny 1986]. Jeschke et al. [2011] introduced a technique to improve curve colorings. Xie et al. [2014] further improved this method by detecting edges in a Laplacian (and/or bi-Laplacian) domain and constructing curves hierarchically. They solve the Laplacian and bi-Laplacian weights using least-squares fitting. In all methods, diffusion curves are placed along the detected edges, and never moved or added in continuous color regions. These methods then rely on optimizing curve coloring for better accuracy.
We introduce a fundamentally complementary solution to the inverse diffusion curve problem. Instead of predetermining curve geometry and optimizing their coloring, we propose doing the opposite by first optimizing the geometry and then determining the coloring accordingly. We demonstrate that with a very simple coloring scheme, our method outperforms prior methods under many situations ( §6.2). Furthermore, our approach accepts input color fields beyond pixel images.
Extensions of diffusion curves. Several methods have been proposed to extend the expressiveness of diffusion curves. Sun et al. [2014] enabled fast diffusion curve cloning and multi-layer composition. Finch et al. [2011] introduced a higher-order notion of smoothness: the colors are defined using a 4th-order linear elliptic PDE rather than a Laplace equation. To accelerate the color evaluation, Boyé et al. [2012] developed a vectorial solver using the Finite Element Method, and Sun et al. [2012] proposed a boundary element based formula, which was later improved in [Ilbery et al. 2013] to handle both Laplacian and bi-Laplacian curves in a unified framework. Higher-order curves offer greater flexibility than the standard diffusion curves, but their inverse problems are more difficult and remain unsolved. In this paper, we focus on the inverse problem for original diffusion curves and discuss potential extension to higher-order domains in §6.3.
Theory and applications of shape optimization. We build our curve optimization on the theoretical foundation of shape optimization [Sokolowski and Zolésio 1992;Haslinger et al. 2003], a subfield of optimal control theory. Mathematically, it solves the problem of finding a bounded set Ω to minimize a continuous functional on Ω. The core idea of shape optimization has been used for image segmentation since the seminal work of [Kass et al. 1988;Mumford and Shah 1989]. It is also related to surface gradient flow widely studied in geometry processing [Schneider and Kobbelt 2001;Crane et al. 2013]. In areas outside of computer graphics, shape optimization has been used to enhance mechanical structures such as airfoils [Mohammadi et al. 2001] and photonic crystals [Burger et al. 2004]. It has also been used in computer vision for image segmentation (e.g., [Herbulot et al. 2006;Jung et al. 2012]). To our knowledge, shape optimization has not yet been applied in vector graphics. In this paper, we solve a shape optimization problem with a PDE constraint ( §5.1), which is significantly more challenging than a conventional shape optimization problem.

Background and Overview
We start by briefly revisiting the mathematical formulation of diffusion curve images. We then present the main focus of this work, the inverse diffusion curve problem, and overview our proposed solution.
Diffusion curve images. In a diffusion curve image, as originally formulated in Jeschke et al. 2009], the color field u is a harmonic function, satisfying a Laplace equation with a Dirichlet boundary condition: where the boundary consists of the entire set of diffusion curves; C and C r specify the colors on the left and right side of each curve, respectively. Typically, both the shapes of the curves and their left-and right-side colors are specified by the user, and the entire color field is uniquely determined by solving the Laplace equation (1).
Since its invention, diffusion curves have been augmented. Orzan et al. [2008] proposed to apply per-pixel blurring to the rasterized image of u, the solution of (1). Finch et al. [2011] further extended to diffuse colors using higher-order elliptic PDEs such as the biharmonic equations.
Inverse diffusion curve problem. While plenty of extensions of the forward diffusion curve problem have been proposed, largely underexplored is the inverse problem, one that computes a set of diffusion curves such that the resulting vector image closely resembles a user-provided 2D color field. In this paper, we address this inverse problem. In particular, we note that the inverse problem involves two subproblems: • Curve geometry. To build a diffusion curve image, one needs to decide where to place the curves (namely, to determine ).
• Curve coloring. Given the curve geometry, the colors on both sides of each curve (namely C and C r ) need to be specified.
As discussed in §2, recent work [Jeschke et al. 2011;Xie et al. 2014] has largely focused on optimizing curve coloring with their geometry predetermined (using edge detection). In contrast, we focus on a complementary problem, the problem of directly optimizing curve geometry. We demonstrate (in §6) that curves with optimized geometries generally yield higher-quality of reconstructions, regardless of the curve coloring schemes.

Overview
Pipeline. We develop a complete pipeline for automatic creation of diffusion curve images ( §4). Figure 2 shows an overview of our pipeline. We take as input a color field I allowing to query for color values at for all x ∈ Ω (where Ω denotes the image domain). Starting with extracting a set of boundary curves ( §4.1) indicating jump discontinuities in I, our method generates a set of curves as "initial guesses" ( §4.2) which are then deformed by our core curve optimization algorithm ( §5) to minimize reconstruction error ( §4.3). Lastly, we post-process the deformed curves ( §4.4) to generate final diffusion curve images. (Section 4.1) Figure 2: Our pipeline. The input to our method is a 2D color field. After obtaining a set of boundary curves indicating color jump discontinuities ( §4.1), our method constructs a set of initial curves ( §4.2) and optimizes their shapes and trajectories ( §4.3 and §5). Finally, we post-process the optimized curves and obtain the resulting diffusion curve image ( §4.4).
Curve optimization. As one of our main contributions, the key component of our pipeline is a curve optimization algorithm that deforms diffusion curves to minimize the reconstruction error. In our algorithm, curves are 1D continuous geometries discretized as polylines, while the cost functional includes an integral over a 2D image region. The optimization problem is constrained by a Laplace equation for color diffusion. Our general idea of solving this optimization problem is using iteration: given a set of curves 0 , we construct a velocity field v so that after deforming 0 according to v for a small time step t, the updated set of curves t leads to improved reconstruction accuracy. This process is repeated until convergence. Conceptually, this routine is similar to gradient descent-in each iteration, a local adjustment is made to obtain a better solution. Our main challenge is to derive a proper form for the derivative of the cost function with respect to the shapes of the curves and use it to construct the velocity field v.
In what follows, we will first provide a complete exposition of our end-to-end pipeline (in §4 and outlined in Figure 2), and then describe in details the algorithms and mathematical derivations for the curve optimization (in §5).

Boundary Curves
Provided an input color field I, we start the pipeline by obtaining a set of boundary curves ∂ Ω indicating the outer boundary and jump discontinuities of I. 1 An example color field and corresponding boundary curves are shown in Figure 3-ab. In practice, we obtain the boundary curves ∂ Ω depending on specific representation of the input color field I: • Pixel images. A common way to represent color fields is using standard pixel images. The boundary curves, however, are not uniquely defined in this case. To obtain these curves in practice, we use Canny edge detection similar to ].
• 3D renderings. If the color field is defined by the rendering of a 3D scene, the boundary curves can be obtained by extracting object contours.
• Other vector formats. For input color fields represented in other vector formats (e.g., gradient mesh), ∂ Ω can be determined directly based on the underlying vector representation (e.g., triangle edges).
Please refer to §6.2 for more details and experimental results on boundary curve computation. 1 The necessity of boundary curves is explained in §5.

Curve Initialization
Desired properties. Similar to gradient descent methods, our curve optimization algorithm takes an initial guess to start with. For ensuring high-quality optimization results, there are a few properties required of the initial curves: 1. Easy to compute. The curve initialization step should not require intense computation: we rely on the optimization step to refine the shapes of these curves.
2. Good coverage. The initial curves should provide a good coverage to the full image domain Ω, so that the optimization is less prone to local optima.
3. Being well-shaped. The initial curves need to be well-shaped. For example, they should have low complexities and not selfintersect or collide with the boundary curves.
To achieve these properties, we use iso-contours of the residual field as the initial curves: In (2), u 0 is given by the diffusion curve image using only the input boundary curves ∂ Ω. These curves can be computed easily from a set of iso-values (Property 1). In addition, as long as the iso-values are distributed properly, the resulting iso-contours will provide a good coverage to the image domain while being well shaped (Properties 2 and 3). For example, to make the initial curves never intersecting with ∂ Ω, we can simply pick strictly positive iso-values as R 0 (∂ Ω; x) = 0 for all x ∈ ∂ Ω.
Our approach. To choose a set of properly distributed isovalues, we start with sampling a set of points in Ω and stacking the residual values (2) at these points into a vector r 0 in ascending order (lines 2 and 3 in Algorithm 1). The resulting vector r 0 provides a picture on the distribution of residuals. We adopt two complementary schemes, global and local, to set iso-values using r 0 , and thereby obtain the iso-contours:

Algorithm 1 Diffusion curve initialization
Require: Color field I (defined on Ω) and boundary curves ∂ Ω 1: procedure CURVEINIT(scheme, ∂ Ω, I, Ω) return iso-contours with iso-values specified in A 11: end procedure Figure 4: Fitting r 0 with 9 components (indicated with purple dots) using a piecewise linear functionr with 2 pieces (m = 1). The value ofr at its internal boundary is selected as iso-value.
• Global. The global scheme constructs a relatively large set of initial curves over the entire domain Ω. Assume that the number of iso-values m is given. Ideally, we would like to find m values such that the consequent iso-contours optimally capture the structure of 2D residual field R 0 . In practice, we solve this problem approximately and rely on our curve optimization algorithm to refine the curves. Particularly, we solve a well-studied 1D problem [Ramer 1972]: to fit a piecewise linear functionr with m+1 pieces that closely describes r 0 (interpreted as a polyline). Then, the values ofr at its m internal piece boundaries are used as iso-values (see Figure 4).
• Local. The local scheme, in contrast to the global one, adds curves locally in regions with high approximation error. In this case, we use only one iso-value determined based on the maximal sampled residual (line 8 of Algorithm 1).
In our curve placement algorithm (detailed in §4.3), we use the global scheme at the beginning to ensure that the initial curves provide a good coverage to the domain Ω (Property 2). Then, the local scheme is applied iteratively to add small sets of curves in high-residual areas. The combination of both schemes offers sufficient approximation accuracy without introducing unnecessarily complex curves (Property 3). We find that m = 2 works well in our experiments.

Curve Placement
Given the initial curves generated by Algorithm 1, our curve optimization algorithm iteratively refines their trajectories to reduce reconstruction errors and finalize curve geometry. We postpone the details of this algorithm (Algorithm 3) and its derivations until §5, but present here the complete curve placement steps (Algorithm 2).
The curve placement process is built on the algorithms of curve initialization and optimization schemes. It takes as input the target color field I defined on domain Ω, the previously obtained boundary curves ∂ Ω, and a tolerance ε 0 on reconstruction error. Based on ∂ Ω, we partition the domain Ω into a number of connected components and process them individually in parallel (line 2 of Algorithm 2). Figure 5 illustrates example boundaries and resulting partitionings. Notice that our approach allows boundary curves to exist inside individual components (e.g., Fig-(  ure 5-b)-these curves will remain fixed throughout the entire pipeline.
For each connected component C, our approach generates diffusion curves via several passes, each of which involves initializing a set of curves (Algorithm 1) and optimizing their shapes (Algorithm 3 in §5). In the first pass, we start with initial curves constructed using the global scheme (lines 5 and 6). After this pass, if the approximation error remains beyond a tolerance ε 0 , additional passes are used in which new curves are initialized using the local scheme (lines 8 to 10). After the error drops below the threshold, we perform a final pass (line 12) in which all curves created in previous passes are optimized together. Finally, we post-process the resulting curves to remove redundant curve segments (line 13 and §4.4). An example of this curve placement process is illustrated in Figure 6.

Curve Post-Processing
Lastly, we post-process the curves returned by Algorithm 2 and generate the final diffusion curve image.
Curve Coloring. Notice that Algorithm 2 returns optimized curve geometry instead of actual diffusion curves. Thus, to turn into a set of diffusion curves, their coloring, namely colors on both sides of each curve, needs to be provided. This corresponds to specifying the values of C and C r in (1).

Algorithm 2 Diffusion curve placement
Require: Color field I (defined on Ω) and boundary curves ∂ Ω 1: procedure CURVEPLACEMENT(∂ Ω, I, Ω, ε 0 ) 2: partition Ω into connected components 3: for each component C do As aforementioned, this curve coloring step is completely orthogonal and complementary to our core technique (Algorithm 2). Thus, in the rest of this paper, we use a simple scheme which directly sample color values on both sides of each curve from the input color field I. That is, for any x ∈ , we set where n and n r respectively denote normal directions pointing left and right side of a point x on a curve (thus, n l = −n r ) and δ is a small positive number that can be set to the size of one pixel when I is represented as a pixel image. Our experiments demonstrate that this simple scheme can yield high-quality results thanks to our optimized curve geometry ( §6.2). In §6.3, we show that more advanced coloring techniques can further improve reconstruction accuracy.
Removing redundant curve segments. As mentioned in §5.4, we represent diffusion curves as polylines consisting of a number of line segments. Some of these segments, however, may be unnecessary. Note that the colors across a line segment are continuous because of the boundary condition (5) on . If the color gradient normal to a segment is also continuous across, then the segment as a boundary has no influence on the solution color field u. A mathematical explanation is in §3 of the supplementary document. Precisely, a normal gradient is continuous when is zero. In practice, we solve u using the Finite Element Method ( §5.4) and check if |d n (x)| at the center point x of each segment is below a threshold. If so, we mark the segment as unnecessary.
Lastly, for each curve output by Algorithm 2, we remove a largest set of connected redundant segments to avoid breaking the curve into many small disconnected components.
To transform the final polyline into a standard diffusion curve made from end-to-end connected Bézier curves, we adopt the Potrace algorithm [Selinger 2003], which was also used in ].
Per-pixel blurring (optional). The curves placed in a smooth color region have continuous color values across the curves. However, since these curves serve as boundaries in the Laplace solve, color gradients may not necessarily remain continuous across curves generated by Algorithm 2. Such gradient discontinuities can sometimes lead to noticeable artifacts [Finch et al. 2011]. Thus, our pipeline includes an optional step following the original framework of diffusion curves ] to perform per-pixel blurring on the rasterized image. The size of blur kernel at each pixel is determined by another Laplace equation: where K 0 (x) gives the desired kernel size along the curves. In particular, we set K 0 (x) = 0 for all x ∈ ∂ Ω since the boundaries and discontinuities in the input color field should never be blurred. For x ∈ , the value d n (x) indicates the magnitude of the gradient domain discontinuity. Thus, we set K 0 (x) = b |d n (x)| a for all x ∈ , where a and b are two global parameters. In our implementation, we set a = 0.2 and b to 5% of the longest axis of Ω's bounding box.
Notice that more advanced curve coloring techniques, such as [Xie et al. 2014], may optimize color gradients across the curves, largely removing gradient discontinuity artifacts. In this case, per-pixel blurring is unnecessary (see §6.3).

Color value (input)
Domain outer boundary (input) Diffusion curves (input/output) Figure 7: Input and output of our curve optimization algorithm. Input: color field I, domain Ω and its outer boundary ∂ Ω, initial diffusion curves ; output: refined curves .

Diffusion Curve Optimization
We now detail the core of our pipeline, the optimization of diffusion curve geometries to approximate a given 2D color field. We first describe an algorithm minimizing the approximation error of diffusion curve images ( §5.1-5.2), and then extend it to balance accuracy against curve length ( §5.3). Lastly, we provide implementation details ( §5.4), followed by the discussions of further extensions ( §5.5).
We introduce Shape Optimization [Sokolowski and Zolésio 1992] to formulate the inverse diffusion curve problem. While building our approach on existing shape optimization concepts and theories ( §5.2), we also develop a new formula for regularizing curve length ( §5.3). Please refer to the supplementary document for complete derivations and a review of related background.
Curve Optimization in a Nutshell. The major steps of our approach are outlined in Algorithm 3. Its input includes the color field I, a 2D closed domain Ω over which I is defined, and a set of initial curves in Ω (Figure 7). In this section, the color field I is treated as a black box, allowing I(x) and ∇I(x) to be evaluated for any x ∈ Ω. Our curve optimization algorithm then iteratively refines the curves by changing their shapes (i.e., the trajectories) and topologies to obtain better approximation. The resulting diffusion curve image consists of the optimized curves and the domain boundary ∂ Ω. During the optimization process, the colors along both sides of these curves (i.e., C and C r of the Laplace equation (1)) are sampled from the given color filed I, and the approximated color value u(x) for all x ∈ Ω is determined according to the equation (1) with the Dirichlet boundary condition u(x) = I(x), ∀x ∈ ∪ ∂ Ω.
We note that rather than sampling color values along the curves, prior methods [Jeschke et al. 2011;Xie et al. 2014] post-optimize color values after the curves are determined. We will discuss the extension of our method to incorporate their post-optimization later ( §5.5) and examine it in our experiments ( §6.3).

PDE-Constrained Optimization Problem
Formally, our iterative curve optimization process minimizes a cost functional defined as the L 2 residual of the color approximation, Here u is the color field determined by diffusion curves. We write as a parameter of R to emphasize the dependence of the residual on through the Dirichlet boundary condition (5). Since u is the solution of the Laplace equation (1), we are concerned with an optimization problem with a PDE constraint, min R(Ω; ) s.t. u satisfies the Laplace eqn. .
PDE-constrained optimization problems are known to be challenging in general [Pinnau and Ulbrich 2008]. In our problem (7),

After timestep
Higher cost Lower cost Figure 8: Curve optimization. Given a set of curves 0 , we construct a velocity field v so that if one deforms 0 according to v, the resulting curves t provide a lower cost.
the optimization variables are the shapes of diffusion curves, that is, the spatial trajectories and topologies of the curves. Ideally, a curve can have an arbitrarily continuous trajectory, and therefore needs to be represented using a continuous functional rather than using individual and discrete parameters. More importantly, the error residual R(Ω; ) depends on the optimization variables (the curves) through the Laplace equation (1) in a complex manner: any local change to the curves has a global impact, one that changes u over the entire domain Ω, which further affects the residual via (6).

Gradient-Descent Solver
We propose a new approach for solving the curve optimization problem (7), following the general spirit of gradient descent.
Starting from a set of initial curves, our approach iteratively decreases the residual (6) by adjusting their shapes. Throughout, a fundamental difficulty we need to address is the computation of the residual's "gradient" with respect to the shapes of the curves, as the conventional gradient in terms of continuous curves is undefined.
We develop our method from the perspective of functional analysis: in each gradient-descent step, we first construct a velocity field v on the curves, specifying v(x) for all x ∈ (Figure 8). We then use v(x) to deform the curves, analogous to a (2D) surface flow in geometry processing [Sethian et al. 2003;Brakke 1992]. In other words, we evolve the curves via a single step of the forward Euler method of integratingẋ = v(x), ∀x ∈ .
In this subsection, we present the details of computing such a v that after deforming the curves accordingly, the residual is guaranteed to decrease (Lines 5-7 of Algorithm 3). Briefly speaking, we will first assume that v(x) is known and analytically express how much the residual would change if the curve is deformed according to v. This analytical expression allows us to formulate the condition of v resulting in a decrease of the residual, and thereby provides us a recipe for computing v.

Algorithm 3 Gradient-descent diffusion curve optimization
Require: initial curves 0 , color field I on Ω with boundary ∂ Ω 1: procedure CURVEOPT( 0 , ∂ Ω, I, Ω) compute v n (x) = −B R (x) using (12) 8: forward-Euler curve advancement, ← + v n t 9: evaluate R using (6), and update its change ∆R 10: end while 11: return current curves 12: end procedure Fréchet derivative as a linear form. Given a domain Ω and a set of initial curves 0 , we consider a general cost functional, where y is continuous on Ω and may depend on the choice of 0 . Our residual (6) takes the form y(x; 0 ) = 1 2 (u(x) − I(x)) 2 and depends on 0 via the Laplace solution u. Assuming a known v, we introduce the Fréchet derivative [Coleman 2012] of C with respect to v. Let t denote the curves evolved according to v after an infinitesimal time period of t, that is, x → x + v(x) t for all x ∈ 0 (Figure 8). The Fréchet derivative of C is a linear form of v satisfying that dC(Ω; 0 ) = lim Conceptually, this derivative measures how quickly the cost functional C changes as we deform the curves using v infinitesimally. According to Hadamard-Zeloésio Structure Theorem [Delfour and Zolésio 2011], such a linear form always exists when Ω, 0 and v are sufficiently regular, which is usually the case in practice. For our cost functional (6), we further reduce the Fréchet derivative into a linear form expressed as a boundary integral where v n (x) := v(x) · n(x) denotes the normal velocity on the curves (Figure 8), Γ 0 = ∂ Ω∪ 0 includes both the domain boundary ∂ Ω and all the inner curves 0 (see Figure 7), and B is another function independent from v but related to the specific integrand y. In the rest of this subsection, we aim to derive a formula to evaluate B(x) for any x ∈ 0 .
Once B is known, setting guarantees a negative derivative value in (9) (assuming that B does not vanish everywhere on 0 ). This provides a formula of constructing v n , which we then apply to deform the curves 0 .
With a sufficiently small timestep size t, the deformed curves t , computed by x + v n (x) t for all x ∈ 0 , is guaranteed by construction to yield a smaller residual value and thus a better approximation of I. Computational Recipe. Shape Optimization Theory has provided a simple recipe of computing B for our particular cost functional (6). Here we simply present the formulas. Please see Appendix A for an outline of the derivation and §3 of the supplementary document for more details.
We first solve the Laplace equation (1) to compute u(x), which in turn allows us to construct a Poisson equation with a Dirichlet boundary condition, Next, the solution p of this equation, together with u, allows the computation of B(x) in a simple form Combining (12) and (10) computes the normal velocity v n , the velocity that can deform the curves 0 and lead to a decrease of the approximation residual (6). This computation is performed at each gradient-decent step, and the optimization process stops when the residual change drops below a threshold ε. Figure 10 illustrates the optimization process with synthetic examples.

Regularizing Curve Complexity
So far, our optimization problem (7) focuses solely on minimizing the L 2 residual (6). However, because the L 2 error along a curve is always zero due to the boundary condition (5), one simple way to yield a very low residual is to use space-filling curves. Indeed, if we start with one curve in a complex color region, it becomes zigzag after running the optimization for many iterations (Figure 9-a). While the numerical residual is low for such curves, their largely increased geometric complexity may be undesirable for certain applications (such as vector graphics editing). Thus, we propose an extension to the cost functional (6), providing users the flexibility to trade approximation accuracy for simpler curves. To this end, we add a regularization term to (6) to penalize the total length of the curves: where α is a user-specified scalar controlling the strength of regularization. It can be shown that similar to (9), the Fréchet derivative of the second term is also a linear form of v n . Let where κ(x) measures the curvature of a point x on the curves. This formula has been used to derive the mean curvature flow [Mantegazza 2011] in geometry processing. It is also a special case of the Fréchet derivative of a general boundary integral (see §1.3 of the supplementary document). Following the derivation of B in §5.2, we obtain the normal velocity for decreas-ingR(Ω; ), that is, v n (x) = −B R (x) − ακ(x). With this slightly different velocity formula, the entire optimization algorithm remains the same as before. In addition, the user is able to control the complexity of resulting curves by adjusting the strength of regularization (Figure 9).

Implementation Details
We now present implementation details of Algorithm 3, wherein two major steps are solving the Laplace equation (1) and the Poisson's equation (11). Both PDEs have Dirichlet boundary conditions defined on the boundary of Ω and the optimized curves (recall (5)). Since we also need to evaluate the domain integral over Ω during the iterations (Line 9 of Algorithm 3), we triangulate the entire domain of Ω and use the Finite Element Method [Zienkiewicz and Morice 1971] for both solves, while other numerical solvers (e.g., the Boundary Element Method) could also be applied.
Finite element discretization. We discretize the boundary and optimized curves into piecewise linear segments, and represent them using polylines. The velocity v n is discretized and stored at every vertex along the polylines. We use the package Triangle [Shewchuk 1996] to triangulate the domain Ω (Line 4 of Algorithm 3). The resulting triangle mesh is then used in the finite element solves. The computation of curves' normal velocity in (12) involves boundary normal derivatives of the finite element solutions (i.e., ∂ p/∂ n and ∂ u/∂ n). We choose the secondorder finite element basis, as it offers higher accuracy especially near the boundary (see §4 of the supplementary document).
Curve tracking. Advancing the curves using the computed normal velocity (i.e., computing t given 0 and v n in Figure 8) is a typical yet nontrivial surface tracking problem. We use a recently developed explicit tracking approach [Brochu and Bridson 2009], which advances the vertices on curve polylines using explicit forward Euler method, and then carefully remeshes the polylines to ensure correct topology changes and a collision-free state.
Timestep size t. To ensure robust curve tracking, we dynamically set the timestep size t for the forward-Euler curve advancement (Line 8 in Algoirthm 3). We start with choosing a t value such that the vertex displacement v n (x) t, ∀x ∈ would not collapse any polyline segment on . This ensures that possible topology changes can be robustly processed. From this starting value, we iteratively halve t until the residual value (after a step of curve deformation) decreases.

Discussions
Measuring geometric complexities. In Equation (13) and the rest of this paper, we use the total length of all diffusion curves to measure their geometric complexity. Depending on specific applications, there may exist other metrics more suitable to user needs. As an example, in §2 of the supplementary document, we discuss another possible measure which can also be incorporated in our curve optimization framework. Coloring schemes. As described at the beginning of this section, given the curve geometry , we specify colors on both sides of each curve by directly sampling color values from the input color field I. Alternatively, prior work [Jeschke et al. 2011;Xie et al. 2014] propose to post-optimize curve colors for better reconstruction accuracy. Our method can easily adopt this approach, post-optimizing the colors after the curves are optimized. We implemented this approach and present the results in §6.3.
Higher-order domains and curves. While our approach focuses on solving the inverse problem of the standard (first-order) diffusion curves, it can be also applied to higher-order domains. For instance, as demonstrated in §6.3, we can feed ∇I instead of I to Algorithm 3 to compute curves offering a higher order of smoothness.
In principle, it is possible to generalize (6) with u(x) directly given by higher-order (e.g., biharmonic diffusion) curves. However, this would dramatically complicate the form of the Fréchet

Figure 10: Synthetic validation of our diffusion curve optimization algorithm. The left-most column contains input color fields given by one closed curve (top) and one open curve (bottom). a varying number of iterations. Our resulting curves precisely match the original.
derivate (9) and thus that of the velocity field (10), causing them significantly more difficult to evaluate numerically. Therefore, we leave the extension of (6) to higher-order curves as future work.

Results
We first ( §6.1) show experimental results demonstrating the validity of our curve optimization algorithm as well as how the regularization behaves in practice. Then, in §6.2, we show reconstructed diffusion curve images using input color fields represented in three forms: pixel image, 3D renderings, and gradient meshes. In addition, we show preliminary results motivating possible future applications.

Experimental Results
Synthetic validations. We design two synthetic tests (Figure 10) to validate our diffusion curve optimization algorithm (Algorithm 3). In these tests, the input color fields I are themselves diffusion curve images with continuous colors. In this case, the optimal is simply the set of diffusion curves used to generate I. Although the shape optimization problem is in general non-convex, our method successfully finds the optimal solutions for both closed and open curves. Please see the accompanying video for curve deformation animations.
Regularizing curve complexity. As discussed in §5.3, our method is able to balance resulting accuracy and curve complexity by varying the strength α of regularization in (13). Figure 11 shows how α influences resulting curves generated with our full pipeline (Algorithm 2). Figure 11-a has simpler curves (due to greater α), but some highlights at the bottom-left of the image are absent. Figure 11-b, on the other hand, provides lower approximation error but at the cost of greater curve complexity vious edge detection based methods Xie et al. 2014]. 2 The parameters for each method are selected such that the resulting curves have approximately identical complexities. Our method outperforms previous ones when handling smoothly varying color fields. Notice that, in the bottom two rows (i.e., Flower and Jade), [Xie et al. 2014] has slightly higher approximation errors (measured in RMSE), because it requires higher curve complexities to work properly in these cases. Please see the supplemental material for more comparisons.
3D renderings. Another kind of color field common to computer graphics applications is renderings of 3D scenes. Our approach can be used to approximate these color fields with diffusion curve images. In this case, we represent I as high-resolution pixel images and obtain the boundary curves ∂ Ω directly using object contours. 3 Figure 15 demonstrates diffusion curve images generated using our method from 3D renderings, where slightly higher curve complexities (compared to Figure 14) are permitted to ensure low reconstruction errors. Our method successfully captures detailed appearances: from glossy surfaces to smooth shadow boundaries.

Gradient meshes.
Our approach can also generate diffusion curve images directly from input color fields I represented using other vector formats. We demonstrate this using input color fields represented as gradient meshes in SVG format [Bah 2011] where each mesh grid is a Coons Patch [Coons 1967]. In this case, the color I(x) and gradient ∇I(x) can be evaluated analytically for any x ∈ Ω, and the boundary curves ∂ Ω are simply the  (resulting from a lower α). In all our results, the complexity is numerically defined as total curve length normalized, so that the longest axis of each image's bounding box has unit length.

Main Results
We now show diffusion curve images generated using our method (Algorithm 2). All our results utilize the per-pixel blurring described in §4.4. Please refer to the supplemental material for unblurred versions.
Theoretically, our approach does not require the input color field I to have any particular representation or discretization. In practice, we demonstrate such flexibility using three types of input: pixel images, 3D renderings, and gradient meshes. The execution time for generating each of these results is summarized in Table 1. The supplemental video contains animations demonstrating the creation process for these results.
Pixel images. One common way to represent a color field I is to use standard pixel images. In this case, I(x) is evaluated using bilinear interpolation, and ∇I(x) using finite difference. As stated in §4.1, we perform edge detection to obtain the boundary curves Conclusion. This paper introduces a novel solution to the inverse diffusion curve problem. The key component of our approach is a curve optimization algorithm that iteratively deforms a set of diffusion curves in a way that the reduction of approximation error is guaranteed. Based upon the core algorithm, we develop a full pipeline that takes an input color field plus a set of boundary curves and produces an image with well-shaped and clean curves that closely matches the input. Our approach offers the generality to take input presented in different formats, which we demonstrate using three types: pixel images, 3D renderings, and gradient meshes.

A Brief Derivation of Equation 12
We now briefly outline the derivation of (12). This derivation has been developed in Shape Optimization Theory. We therefore refer to the supplementary document and Chapter 10.6 of the book [Sokolowski and Zolésio 1992] for a detailed exposition. First, the Fréchet derivative of C with a general integrand function y can be expressed as where y is the so-called shape derivative of y under a given velocity field v, defined as The concept of shape derivative is very much analogous to those in continuum mechanics [Bonet and Wood 1997], which has been widely adopted for creating computer animations. In particular, y is equivalent to the material derivative, measuring the change of y in the undeformed (material) space, while y indicates the derivative value in the deformed space, that is, the change rate of y due to the boundary changes only. See §1.2 of the supplementary document for a rigorous mathematical definition ofẏ.
We notice that because of the Dirichlet boundary condition (5), the approximated color u agrees with the input I for all x on the boundary (i.e. x ∈ Γ 0 ). Therefore, when y(x; 0 ) = (u(x) − I(x)) 2 , the second integral term in (15) vanishes, leaving only the first domain integral term, Here the shape derivative u follows the same definition as in (16); we use the fact that the shape derivatives, just like the conventional ones, satisfy the chain rule.
One can prove that the shape derivative u satisfies another Laplace equation (see §1.5 of the supplementary document), where the Dirichlet boundary condition is determined by the normal derivative of both the provided (i.e., I) and approximated (i.e., u) color fields on the same boundaries. Since the Laplacian operator is self-adjoint and u is used in an integral (17), instead of solving the Laplace equation (18), we solve its adjoint problem (11), whose solution enables us to transform the domain integral (17) into a desired boundary integral (9), because ∂ Ω required by Algorithm 2. Although color discontinuities are not well defined for standard pixel images, our method in practice is robust on the choice of boundary curves. Figure 13 shows three examples with boundary curves detected using Canny detector with three thresholds. Notice how missing boundaries (when increasing the threshold) are handled by additional curves generated by Algorithm 2. All our results for standard pixel input used thresholds between 0.1 and 0.2. Figure 14 contains diffusion curve images reconstructed from pixel input (Figures 3 and 12) using Algorithm 2 as well as previous edge detection based methods Xie et al. 2014]. 2 The parameters for each method are selected such that the resulting curves have approximately identical complexities. Our method outperforms previous ones when handling smoothly varying color fields. Notice that, in the bottom two rows (i.e., Flower and Jade), [Xie et al. 2014] has slightly higher approximation errors (measured in RMSE), because it requires higher curve complexities to work properly in these cases. Please see the supplemental material for more comparisons.
3D renderings. Another kind of color field common to computer graphics applications is renderings of 3D scenes. Our approach can be used to approximate these color fields with diffusion curve images. In this case, we represent I as high-resolution pixel images and obtain the boundary curves ∂ Ω directly using object contours. 3 Figure 15 demonstrates diffusion curve images generated using our method from 3D renderings, where slightly higher curve complexities (compared to Figure 14) are permitted to ensure low reconstruction errors. Our method successfully captures detailed appearances: from glossy surfaces to smooth 2 We thank the authors of [Xie et al. 2014] for confirming the correctness of results in Figure 14 generated with their approach. 3 In this example, we assume all objects to be homogeneous. To handle heterogeneity, ∂ Ω needs to include color jumps across object surfaces.

Input Format Scene Time Scene Time
Pixel images (Figure 14) shadow boundaries.
Gradient meshes. Our approach can also generate diffusion curve images directly from input color fields I represented using other vector formats. We demonstrate this using input color fields represented as gradient meshes in SVG format [Bah 2011] where each mesh grid is a Coons Patch [Coons 1967]. In this case, the color I(x) and gradient ∇I(x) can be evaluated analytically for any x ∈ Ω, and the boundary curves ∂ Ω are simply the mesh boundaries. As shown in Figure 16, our method directly creates diffusion curve images closely approximating the gradient meshes, without having to rasterize the input into pixel formats.

Additional Results
Coloring optimization. As discussed in §5.5, our technique is orthogonal and complementary to coloring optimization techniques which take diffusion curve geometry as input and optimizes colors (and color gradients) at both sides of the curves. These techniques can be used to replace our simple coloring scheme which samples color values directly ( §5.1) and performs per-pixel blurring ( §4.4). Figure 18 demonstrates that our optimized curve geometry can be coupled with the coloring optimization scheme introduced by Xie et al.
[2014] to further reduce reconstruction errors. Notice that, since this technique explicitly optimizes color gradients across each curve, we did not perform per-pixel blurring (as stated in §4.4).
Higher-order domain. As discussed in §5.5, our approach can be applied to higher-order domains for generating curves with higher-order smoothness. For instance, as shown in Figure 19, our pipeline can be applied to color gradients rather than original color values. In other words, given a RGB image I, we can use ∇I, a six-channel image, as the input to Algorithm 2. Given our reconstructed gradient image, we solve an additional least square problem to recover the final image.
However, as observed by Xie et al. [2014], we found that for natural images, solving the optimization at higher-order domains normally does not lead to better approximation accuracy under similar curve complexities. This is because higher-order domains are generally filled with significantly more high-frequency contents that require complex (almost space-filling) curve geometry to accurately reconstruct.
Animated result. Lastly, we show preliminary results to motivate future applications of our approach. Since our method optimizes the shape of diffusion curves iteratively, it is suitable for generating animated results from a sequence of gradually changing input color fields. The basic idea is curve reusing: taking optimized curve geometry from one frame as the initial configuration to "warm start" the next one. Figure 20 and the accompanying video show a proof-of-concept example. The input is the relighting (i.e., the object stays static while the light source moves) of a shiny torus knot. In this case, the boundary curves keep unchanged throughout all frames, and optimized curve geometry from one frame remains valid for all other frames. Previous methods Xie et al. 2014] cannot easily enforce curve coherence across different frames, leading to temporally noisy animations. By modifying the curve initialization step in Algorithm 2 to reuse optimized curve geometry, we are able to accelerate the optimization process by 2.3×, and the resulting animation has lower approximation error and little noise. Please see the supplementary video for full animations.

Limitation and Conclusion
Limitation. Our approach has a few limitations that can inspire future work. First, it requires the color field to be C 0 continuous everywhere except at the given boundaries. Robustly finding clean boundary curves, however, can be challenging. Second, if the color field contains spatially high-frequency features, very fine triangulation may be needed to fully resolve them, slowing down our optimization process.
Conclusion. This paper introduces a novel solution to the inverse diffusion curve problem. The key component of our approach is a curve optimization algorithm that iteratively deforms a set of diffusion curves in a way that the reduction of approximation error is guaranteed. Based upon the core algorithm, we develop a full pipeline that takes an input color field plus a set of boundary curves and produces an image with well-shaped and clean curves that closely matches the input. Our approach offers the generality to take input presented in different formats, which we demonstrate using three types: pixel images, 3D renderings, and gradient meshes.

A Brief Derivation of Equation 12
We now briefly outline the derivation of (12 where y is the so-called shape derivative of y under a given velocity field v, defined as The concept of shape derivative is very much analogous to those in continuum mechanics [Bonet and Wood 1997], which has been widely adopted for creating computer animations. In particular,ẏ is equivalent to the material derivative, measuring the change of y in the undeformed (material) space, while y indicates the derivative value in the deformed space, that is, the change rate of y due to the boundary changes only. See §1.2 of the supplementary document for a rigorous mathematical definition ofẏ.
We notice that because of the Dirichlet boundary condition (5), the approximated color u agrees with the input I for all x on the boundary (i.e. x ∈ Γ 0 ). Therefore, when y(x; 0 ) = (u(x) − I(x)) 2 , the second integral term in (15) Here the shape derivative u follows the same definition as in (16); we use the fact that the shape derivatives, just like the conventional ones, satisfy the chain rule.
One can prove that the shape derivative u satisfies another Laplace equation (see §1.5 of the supplementary document), where the Dirichlet boundary condition is determined by the normal derivative of both the provided (i.e., I) and approximated (i.e., u) color fields on the same boundaries. Since the Laplacian operator is self-adjoint and u is used in an integral (17), instead of solving the Laplace equation (18), we solve its adjoint problem (11), whose solution enables us to transform the domain integral (17) into a desired boundary integral (9), because where the last equality follows the integration by parts and Green's formula. Further, the last two integral terms vanish, following the fact that both p (according to the boundary condition of (11)) and ∇u (according to (18)) are zero on the boundary Γ 0 . Eventually, only the first integral term in the last expression remains, yielding the Fréchet derivative of R as a boundary linear form of v n , that is, dR(Ω; 0 ) = Γ 0 B(x) v n (x) dΓ with B(x) expressed as in (12).

Figure 14:
Comparisons between diffusion curve images generated by our approach and previous methods using pixel images (top row) as input. The parameters are adjusted so that the resulting curves generated by each method have roughly identical complexities. Our approach not only yields lower approximation error (measured in RMSE), but also generates better shaped and relatively simple curves. Figure 14: Comparisons between diffusion curve images generated by our approach and previous methods using pixel images (top row) as input. The parameters are adjusted so that the resulting curves generated by each method have roughly identical complexities. Our approach not only yields lower approximation error (measured in RMSE), but also generates better shaped and relatively simple curves.

Reference
Ours Curves Reference Ours Curves

Twill Knots
Reference Ours Curves Reference Ours Curves

Wobble Chess
Reference Ours Curves Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Apple Tomato
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves

Mango Candle
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries. Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Cornell Box Carved
Reference Ours Curves Reference Ours Curves

Twill Knots
Reference Ours Curves Reference Ours Curves

Wobble Chess
Reference Ours Curves Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Apple Tomato
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves

Mango Candle
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries. Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Apple Tomato
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves

Mango Candle
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries. Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries.


Using local redundancy to improve GNSS absolute positioning in harsh scenario

In GNSS context, absolute positioning is a widespread operational mode, largely used in many field activities as automotive, aerospace and ships navigation. The functional model of absolute positioning, relating pseudorange measurements to unknowns, is well defined, while its stochastic model, defining the behavior of measurement error, is currently under investigation and is the core of this work.A realistic model of measurement error can be used to define a weighting scheme, able to reduce the positioning errors; a weighted approach is mainly necessary in signal-degraded scenario, where pseudorange accuracies are significantly different and equally weighting all the measurements would yield very large errors.The quality indicators, commonly adopted to define a pseudorange error model, are the signal-to-noise ratio and the satellite elevation angle. In this work, an additional indicator, the redundancy number, is introduced.The redundancy numbers are the diagonal elements of the redundancy matrix and they represent the degree of controllability of the measurements; a large value of the redundancy number corresponds to a well-controlled measurement, while a small one corresponds to a leverage observation, with a high potential to influence the solution.The benefits of using the proposed weighting schemes is demonstrated in urban canyon and with single (GPS only) or multiple (GPS/Glonass) constellations

INTRODUCTION
In GNSS context, scenarios characterized by bad signal quality and/or bad satellite geometry, usually due to natural or artificial obstacles, are considered harsh environments. Mountainous or highly urbanized areas are examples of harsh environments. In urban scenarios, multipath interference and non-line-of-sight (NLOS) phenomena typically affect the GNSS signals. The multipath interference happens when a signal is received through multiple paths; the combination of multiple signals causes a distortion of the correlation function (between received and locally generated signals) and yields a range error up to several tens of meters [1]. The NLOS reception happens when the direct signal is blocked and only reflected signals are received; this phenomenon is very common in urban scenario, yielding range errors of even km order [2]. The NLOS and multipath phenomena are even more widespread in high sensitivity receivers, which are able to receive very weak signals.
Absolute positioning, also called Single Point Positioning (SPP), is the most common GNSS operational mode. In fact, SPP mode is usually adopted by mass-market receivers, which are mounted on tablet, smartphones, car navigation devices, UAVs [3]- [5]; on these devices, SPP is often combined with other information sources, such as maps, inertial sensors or vision systems [6], [7], [8].
The aforementioned phenomena are responsible for the presence of blunders among the measurements, which can cause large position error. Two approaches are usually adopted to cope with this problem: • Receiver Autonomous Integrity Monitoring (RAIM) techniques, able to detect and reject blunders [9]- [11]; • Robust estimation techniques, able to absorb the effects of the blunders [12]. In this work, it has been followed a different approach, consisting in defining a suitable weighting scheme, able to deweighting potentially dangerous measurements. A capillary survey has been carried out on the existing measurement weighting schemes, which are tested in urban environment; afterwards, an additional quality criterion is proposed in order to limit the effect of potential blunders.
In literature, two quality indicators are typically associated to the GNSS measurements, in order to differently weighting them: satellite elevation angle and carrier-to-noise ratio.
Satellite elevation (El) is considered a quality indicator of a measurement, because satellites at low elevation are usually noisier, due to typical behaviour of multipath, tropospheric and ionospheric errors [1].
The Signal-to-Noise ratio (S/N), or better the Carrier-to-Noise ratio (C/N 0 ), is the ratio between the carrier power and the noise power per unit bandwidth and is usually expressed in Decibel-Hertz; C/N 0 is a measure of signal strength, so it is a suitable indicator of measurement quality [10].
The El and C/N 0 indicators are adoptable individually or synergistically to define a model of measurement error variance; a weighting scheme can be obtained, inverting the variances of the measurement errors.
The positioning solution is differently influenced by each measurement of a dataset; indeed, some measurements, called leverage observations, have high potential to influence the solution, while others are less influent. If a leverage observation is affected by a gross error, a very large solution error could occur; for this reason, the leverage observations are critical. It can be demonstrated that the concept of leverage observation is related to the local redundancy [13], which is well represented by the diagonal elements of redundancy matrix, called redundancy numbers.
In urban scenario, the satellite geometry is often weak and the presence of blunders is very common; in this conditions, a blunder on a leverage observation could strongly degrade the solution.
In this work, an innovative weighting scheme, which takes into account the local redundancy information, is proposed and its effectiveness is shown in a very harsh environment like an urban canyon. The proposed scheme can be superimposed to the existing weighting strategies, which involve the indicators El and C/N 0 , and the benefits of this combined approach are demonstrated.
The proposed weighting scheme is applied to GPS only and to GPS/Glonass cases, demonstrating its usefulness in both configurations.
Below, the SPP operational mode is briefly described; afterwards, the existing weighting schemes, based on El and/or C/N 0 , are described. Then the redundancy matrix, whose diagonal elements are the redundancy numbers, is defined and is linked up to the concept of local redundancy. The possible contribution of the redundancy matrix to a weighting strategy is shown. At the end, the proposed weightings are applied to real data and the results are discussed.

SINGLE POINT POSITIONING
In SPP, the functional model, defining the relationship between measurements and unknowns, is where is the vector of measurements, defined as the difference between measured and computed (with a priori information) pseudorange, is the design matrix, ∆ is the state vector, containing the corrections to update the receiver coordinates and clock offset, ε is the measurement error vector.
Equation (1), also referred to as measurement model, is a linearized equation; the receiver coordinates and clock offset are obtained correcting the a priori information with the estimated state vector.
The stochastic model, expressing the uncertainty of the measurement model, is: where is the variance-covariance matrix of the measurement errors; =1,…, 2 are the measurement error variances; is the number of measurements.
The matrix is diagonal, under the usual assumption of independent measurement errors.
The number of measurements is usually larger than the number of unknowns and the equation (1) is solved using weighted least squares (WLS) method.
WLS solution is: where is the weighting matrix, usually set equal to the inverse of .

WEIGHTING SCHEMES REVIEW
Most of weighting schemes in literature are based on C/N0, satellite elevation El or both; hence, the weighting strategies review is organized according to the considered quality indicators.

Elevation-based weighting schemes
The most widespread variance models for GPS carrier phase, based on satellite elevation, depend on sin ( ) or sin 2 ( ) functions [14], [15]; despite being designed for carrier phase, they were simply adapted for pseudorange observation [16], [17] as shown below: where σ 0 2 is the pseudorange error variance when = 90°. Slight modifications to the previously mentioned models are proposed in literature [18], maintaining the concept that low elevation satellites correspond to noisy measurements: where , , are parameters to be defined according to the receiver and the operational scenario.
A more complex variance model, which takes into account the contributions from several measurement error sources, is defined in [19], [20]; it is designed for Wide Area Augmentation System (WAAS) and it is shown below: where = 0.5 is the variance of the supplied tropo-free and iono-free pseudorange correction, = 0.5 is the variance of the vertical ionosphere correction, is the obliquity factor, projecting vertical ionospheric variance into receiversatellite direction, = 0.22 is the receiver noise variance, 45 = 0.22 is the variance of the multipath contribution at 45 degrees, = 0.15 is the variance of the vertical tropospheric delay estimate.
The model (8) has been adapted to the case of GPS SPP in [21] as shown below: where 2 is the User Range Accuracy (URA) related to the satellite ephemeris and clock, broadcast in the GPS navigation message, 2 is the accuracy related to ionosphere delay after Klobouchar model application, 2 is the accuracy related to troposphere error after correction model application, 2 is the accuracy related to multipath error.

C/N 0 -based weighting schemes
In [14], a variance model for carrier phase observation is proposed; it is re-used and adapted in [2] for pseudorange measurement. The model is exclusively based on C/N0 and has the following form: where c is a constant, empirically defined.
The model of Hartinger and Brunner [15] has been modified by Wieser and Brunner [22], introducing an additive term, and making the model parameters estimation in view of their application to carrier phase observable. Kuusniemi et al. [23], [24] adapt the Wieser and Brunner [22] model to pseudorange, tuning the parameters according to actual measurement errors and used equipment; the obtained variance model is: where a and b are the model parameters.
Kuusniemi [10] has estimated the parameters of model (11) for different scenarios (lightly and heavily degraded signal environment).
The Hartinger and Brunner [14] model has been the basis for Brunner et al. (1999) to build a carrier phase variance model, called sigma-Δ, whose expression is very similar to (10), but with at the exponent a correction term proportional to the difference between the measured and the predicted (according to satellite elevation) C/N0. Li and Wu [26] have modified the model (11), adopting Euler's number as base of the exponentiation and the difference between the measured C/N0 and its minimum threshold as exponent: The parameters of the model are a, b and k and are empirically determined.
Aminian [27] proposed an alternative expression for the pseudorange error variance without exponential form, as shown below: where 0 2 and / 0( ℎ) are respectively the theoretic variance and carrier-to-noise ratio in case of satellite at zenith and open-sky scenario.

Weighting schemes based on both satellite elevation and C/N 0
In literature, few weighting schemes, based on both satellite elevation and C/N0, are available.
In [28] the following model is developed: where The term s 1 is a threshold, used as comparison for the measured C/N0; if / 0 ≥ 1 , the measurement is considered accurate and the corresponding weight is set to 1 (the maximum possible value for a weight), otherwise a complex formula is adopted for the weight. The parameters in (14) and (15) are determined empirically as 0 = 10, = 30, = 30.
In [29] a model, based on both satellite elevation and C/N0, designed for land navigation, is defined: The model (16) is clearly a fusion between the elevationbased model (5) and the C/N0-based model (10). The constant k is 1 if the signal is received in line-of-sight (LOS) and 2 or +∞ if the signal is received after reflection by obstacles surrounding the antenna (NLOS); in [29], a fish-eye camera is used to distinguish LOS and NLOS signals. In this work, a modified model (16) is used, where k is always equal to 1.

REDUNDANCY MATRIX
The residuals , defined as the difference between the actual measurements and their estimated values (�), are an indicator of the measurement mutual agreement [13]: It can be simply demonstrated that the relationship between the measurement errors and the LS residuals, considering equally weighted measurements, is: The matrix is called "Redundancy Matrix" and has the following expression: The matrix has rows and columns, its trace is equal to the total redundancy (or the degree of freedom) of the equation, that is ( − ). The i-th diagonal element of , with = 1, … , is called "redundancy number" and is the contribution of the i-th measurement to the total redundancy [31]. The redundancy number assumes values between 0 and 1. Small values of (near 0) correspond to measurements providing little contribution to total redundancy and so hardly controlled; on the other hand, approximately equal values of are desirable, meaning that every measurement is controllable.
Measurements with small values are leverage observations and have high potential to influence the solution; if a blunder or a large bias is present on a leverage observation, harmful effects can be evident on the positioning. In this context, de-weighting leverage observations could be a successful strategy, above all in environments, as urban canyon, characterized by frequent gross errors.
In Table 1, an example of GPS geometry, in terms of satellite elevations ( ) and azimuths ( ), is shown; from and , the design matrix can be computed and, consequently, the redundancy matrix is obtained. The diagonal elements of are the redundancy numbers.
In Figure 1, the sky-plot of the satellite geometry, detailed in Table 1, is shown; all the listed satellites are considered, so the number of visible GPS satellites is 10 and the PDOP (Position Dilution of Position) value is 1.59. The total redundancy of the measurement model is 6 (= − ), therefore the average redundancy number is 0.6. The sky-plot is accompanied by the values of the redundancy numbers and it is evident that all the values are near the average.
In Figure 2, the sky-plot of the satellite geometry, detailed in Table 1, is shown, but satellites with PRNs 17, 19 and 25 are excluded; consequently, the number of visible GPS satellites is 7 and the PDOP value raises to 2.15. The total redundancy is 3, so the average redundancy number is about 0.43. The redundancy numbers are all around or above the average value, except the one corresponding to satellite 22, which is very low and so can be considered a leverage observation. From Figure  2, it is evident a geometric interpretation of leverage observation, for GPS SPP, as a measurement corresponding to a satellite "isolated" with respect to the others.
In case of multi-constellation approach, a redundancy number is null when a single satellite from a system is present. For instance, in a GPS/Glonass combination, if there is only one GPS measurement and several Glonass, the redundancy number corresponding to the unique GPS observation is 0 and the observation is uncontrollable (because it is the only one able to estimate the clock offset to GPS time). The same concept is valid for any GNSS combination, including the European GNSS, Galileo [30], or the Chinese one BeiDou.

CONTRIBUTION OF REDUNDANCY MATRIX TO WEIGHTING SCHEMES
The diagonal elements of the redundancy matrix, the redundancy numbers , represent the local reliability of the measurement model and so the controllability of the measurements [31].
Measurements with low are critical for the solution, because they can strongly influence it. For this reason, an effective weighting strategy should take into account the redundancy number, above all in scenarios with high probability of blunders.
In [32] the redundancy number is used, together with several other parameters, as an indicator of the measurement quality for robust estimation. In [10] the extra-diagonal elements of the redundancy matrix are used to evaluate the correlation among the measurements, in order to avoid erroneous blunder rejections. Recently, Falco et al. [33] adopted a purely geometric strategy to identify measurements which could strongly influence the solution, in order to de-weight them.
In this paper, the information contained in the redundancy matrix is used to correct the original weighting matrix, which could be defined according to one of the previously described methods. In general, the weight , associated to the i-th measurement, is obtained inverting the pseudorange error variance 2 .
The contribution of a redundancy number to the weighting strategy could be introduced making proportional to as shown below: The redundancy numbers are 0 in case of lack of redundancy and the scheme (20) takes into account this possibility.  Table 1. Example of GPS geometry, in terms of satellite elevation and azimuth.
The strategy is tested on a multi-constellation static data collection, carried out in urban scenario, and demonstrated its efficiency in both GPS only and GPS/Glonass cases, exhibiting significant improvements on positioning accuracy.
The same strategy could be simply extended to GNSS velocity estimation [34].

TEST
To demonstrate the benefits of including the redundancy number into a weighting scheme in difficult environment, a static data collection is carried out in urban scenario; specifically, the test location is the "Centro Direzionale di Napoli" (CDN), a district of Naples (Italy) characterized by several skyscrapers. A specific point in CDN is surveyed with topographical methods, in order to determine its coordinates with mm order. The used equipment consists in a NVS receiver, a single frequency, double constellation (GPS/GLONASS) and high-sensitivity device, connected to a patch antenna; the device is placed on the defined point to collect pseudorange measurements for some hours; the collected data are processed in single point to test several weighting schemes.
In Figure 3, the considered point (indicated as P) is shown and the urban canyon scenario is evident.
The data collection was carried out on 1 st July 2016; the NVS receiver was placed on P point for about 2 hours. In Figure 4, the number of GPS/Glonass visible satellites and the corresponding PDOP values during the session are shown.
The number of GPS visible satellites ranges between 0 and 8, with an average of 6.4. The maximum number of Glonass visible satellites is 5, while the maximum number of combined GPS/Glonass satellites is 12.
The GPS solution availability, defined as the time percentage where the solution can be computed (at least 4 satellites should be visible), is about 99.3%, while the GPS/Glonass solution availability is about 99.5%.
Despite the urban canyon scenario, the solution availability is high, owing to the use of a high-sensitivity receiver, able to acquire very weak signals. The rapid oscillations of the number of visible satellites is due to frequent signal loss, typical of hostile environment.
In GPS only case, the PDOP average value is about 2.7, with maximum and minimum respectively about 15 and 1.6; in GPS/Glonass case, minimum, maximum and mean values are respectively 1.5, 17 and 2.4. The maximum PDOP value in GPS/Glonass case happens when GPS satellites are less than four, so in that time window GPS DOP cannot be computed.

RESULTS
Several weighting strategies are described in section 2; among them, some techniques are chosen to process data collected in urban scenario, in order to compare their performance.
The considered configurations, implementing GNSS SPP with a distinct weighting scheme, are: • The one using weights defined by formula (5), shortly indicated as ELV1; • The one using weights defined by formula (9), shortly indicated as ELV2; • The one using weights defined by formula (11), shortly indicated as CN01; • The one using weights defined by formula (13), shortly indicated as CN02; • The one using weights defined by formula (14), shortly indicated as ELVCN01; • The one using weights defined by formula (16), shortly indicated as ELVCN02. The previously mentioned configurations are compared with the baseline one, where all the measurements are equally weighted, shortly indicated as EQW.
The figure of merits adopted for the comparison are the mean, RMS and maximum errors of both horizontal and vertical components of position.
In Table 2 the configuration performance, in case of GPS only, are summarized.   The results obtained with the EQW configuration is clearly inaccurate; in fact, the errors, on both components and for each figure of merit, are very large owing to the frequent blunders in the dataset. The considered weighting schemes reduce drastically the errors; the models based uniquely on satellite elevation (ELV1 and ELV2) show better performance with respect to models based uniquely on carrier-to-noise ratio (CN01 and CN02). The methods which consider both quality indicators (ELVCN01 and ELVCN02) demonstrate significant improvements relatively to the other configurations. In particular, the model introduced by Tay and Marais [29], i.e. ELVCN02 configuration, shows the best performance for the considered data session; indeed RMS and mean errors are reduced, with respect to EQW configuration, between 58% and 68%, and the maximum horizontal and vertical errors are reduced respectively by 63% and 45%.
The considered weighting schemes are augmented, using information embedded in the redundancy matrix as described in previous sections; the configurations using this information are shortly indicated as RDM. For instance, a configuration using the elevation-based model (5), augmented with redundancy matrix-based approach, is shortly indicated as ELV1+RDM. In Table 3, the performance of the configurations, obtained combining the considered classical weighting schemes with redundancy matrix information, are resumed. Moreover, in the table, the percentage improvements, with respect to corresponding configurations without redundancy matrix use, are reported alongside the errors. It is evident that the introduction into the weighting schemes of information from redundancy matrix brings to significant improvements. The best configuration is ELVCN02+RDM, characterized by horizontal mean, RMS and maximum errors respectively about 4.5, 8.8 and 87.2 meters, a very good performance considering the strongly degraded scenario.
The same analysis, carried out for GPS only, is provided for GPS/Glonass multi-constellation case. In Table 4 the GPS/Glonass configuration performance are summarized. The addition of Glonass measurements to GPS ones does not improve all the considered figure of merits, as Glonass measurements are often very noisy; in detail, GPS/Glonass configurations have better performance in horizontal plane, but worse performance in vertical direction. The errors of EQW configuration are very large, because there are several blunders among the measurements, from both used systems; the Glonass use allows a reduction in maximum horizontal of about 130 meters. As for GPS only case, the weighting schemes based uniquely on satellite elevation (ELV1 and ELV2) demonstrate better performance than models based uniquely on carrier-tonoise ratio (CN01 and CN02); models exploiting both quality indicators show the best results.
In Table 5, the figure of merits of the GPS/Glonass configurations, using classical weighting schemes augmented with redundancy matrix information, are resumed; in the table, the percentage improvements, with respect to corresponding configurations without redundancy matrix use, are reported alongside the errors. The benefits of using redundancy matrix in weighting schemes are evident even in GPS/Glonass case; indeed, all the considered errors are strongly reduced. The best configuration is ELVCN02+RDM, characterized by horizontal mean, RMS and maximum errors respectively about 4.5, 8.1 and 73.6 meters, a very good performance considering the strongly degraded scenario and slightly lower with respect to GPS best configuration.
The use of information derived from redundancy matrix demonstrates very good performance for both GPS only and GPS/Glonass cases. The error reductions are impressive for each considered configuration and figure of merit. In particular, the best configuration is ELVCN02+RDM, using / 0 , satellite elevation and local redundancy as parameters for weighting; the obtained performances are characterized by mean and RMS horizontal errors respectively below 5 and 10 meters, and maximum horizontal error about 80 meters.
The considered approach, taking into account the information contained into the redundancy matrix, demonstrates its effectiveness, owing to the ability to deweighting the geometrically dangerous measurements.

CONCLUSIONS
In GNSS context, some environments, such as urban canyons, are critical, owing to the likely presence of blunders among the measurements, which yields very large position errors. In these conditions, it is essential to adopt a suitable method to differently weight the measurements, in order to limit the blunder effects. In literature, several schemes exist, based on two quality indicators: satellite elevation and carrierto-noise ratio.
In this paper, the most representative weighting schemes for GNSS absolute positioning are tested in urban scenario and their performance are compared. The best results are obtained with the weighting strategies based on both the abovementioned quality indicators.
The concept of local redundancy refers to the degree of controllability of the observations; for instance, a measurement with low local redundancy is difficultly controlled and, if affected by gross errors, it strongly influences the solution. The redundancy matrix relates the measurement errors to the residuals and its diagonal elements are the redundancy numbers, representing the local redundancy. In this work, a strategy to include the local redundancy into the existing weighting schemes is suggested. In the proposed approach, each weight value is set proportionally to the corresponding local redundancy number. In this way, the most influencing measurements (i.e. the leverage measurements) are de-weighted and their effects on the solution are limited.
The considered weighting schemes are applied to process in SPP a 2-hour dataset, collected in urban scenario; both GPS only and multi-constellation GPS/Glonass configurations are analysed. The obtained results demonstrate the effectiveness of the proposed weighting strategy; indeed, the position errors, using the weighting scheme based on redundancy numbers, are significantly reduced for all the considered configurations.
The best performances are attained with a weighting scheme, comprising the two classical quality indicators, satellite elevation and carrier-to-noise ratio, and the redundancy number; in single point positioning, a horizontal accuracy below 10 meters is obtained in a severe urban scenario.