key: cord-0327893-50ldj8j7
authors: Parisotto, Simone; Calatroni, Luca; Caliari, Marco; Schonlieb, Carola-Bibiane; Weickert, Joachim
title: Anisotropic osmosis filtering for shadow removal in images
date: 2018-09-17
journal: nan
DOI: 10.1088/1361-6420/ab08d2
sha: 9113e69103bd93275f8803a8ea41375c713aa9a3
doc_id: 327893
cord_uid: 50ldj8j7

We present an anisotropic extension of the isotropic osmosis model that has been introduced by Weickert et al.~(Weickert, 2013) for visual computing applications, and we adapt it specifically to shadow removal applications. We show that in the integrable setting, linear anisotropic osmosis minimises an energy that involves a suitable quadratic form which models local directional structures. In our shadow removal applications we estimate the local structure via a modified tensor voting approach (Moreno, 2012) and use this information within an anisotropic diffusion inpainting that resembles edge-enhancing anisotropic diffusion inpainting (Weickert, 2006, Gali'c, 2008). Our numerical scheme combines the nonnegativity preserving stencil of Fehrenbach and Mirebeau (Fehrenbach, 2014) with an exact time stepping based on highly accurate polynomial approximations of the matrix exponential. The resulting anisotropic model is tested on several synthetic and natural images corrupted by constant shadows. We show that it outperforms isotropic osmosis, since it does not suffer from blurring artefacts at the shadow boundaries.

The use of partial differential equations (PDEs) has a long tradition in mathematical image processing. In particular, PDEs based on transport and diffusion mechanisms have been considered to model several image reconstruction models suitable for image enhancement, denoising, deblurring, inpainting and segmentation. We refer the reader to the review [15] and the monographs [3, 7, 32, 34, 37] for further references.

Among these models, a very special place is occupied by diffusive PDEs encoding anisotropy, i.e. favouring diffusion along some specific directions only. In [37] nonlinear arXiv:1809.06298v1 [math.AP] 17 Sep 2018 diffusion PDEs with space-variant diffusion tensors are studied. For a regular image domain Ω ⊂ R 2 and a stopping time T > 0, given a degraded image f ∈ L ∞ (Ω; R) and two smoothing parameters ρ, σ > 0, the anisotropic diffusion model in [37] looks for a solution u in a suitable function space satisfying the following initial value problem:

on Ω × (0, T ],

on Ω,

where n is the outward normal unitary vector on ∂Ω and D is a non-constant diffusion tensor satisfying suitable regularity conditions, with eigenvectors inherited from the socalled structure tensor J ρ (∇u σ ). It encodes local directional information of u σ (that is, the image u convolved with a Gaussian kernel of standard deviation σ). More precisely, its eigenvectors point in the directions of largest and smallest contrast averaged over a Gaussian smoothing scale ρ, and the corresponding eigenvalues measure this contrast. The diffusion tensor uses the same eigenvectors, and its eigenvalues are functions of the eigenvalues of the structure tensor. Depending on the application, different models such as edge-enhancing anisotropic diffusion or coherence-enhancing anisotropic diffusion have been proposed [37] Edge-enhancing anisotropic diffusion has been adapted to inpainting problems in [39] . It is particularly useful for sparse inpainting problems encoutered e.g. in inpainting-based compression applications where it outperforms other PDE approaches [13, 33] . Beyond anisotropic PDEs, non-smooth anisotropic regularisers for variational imaging models are also considered. In [14, 21, 22, 29, 28] , for instance, directionality is used to define the anisotropic Total-Variation and Total-Generalised-Variation functionals. This is classically done by considering a re-parametrised version of the gradient operator depending on the local orientation of the image, which allows to enforce diffusion along certain directions only. In terms of variational models, one replaces the squared norm ∇ u∇u by a quadratic form of type ∇ uD∇u. This has a very long tradition in image analysis [26] .

The explicit dependence of these anisotropic models on local terms such as position and local directions of the image makes the analysis of these models more challenging [37, 14, 22, 29] . Also from a numerical point of view, the design of suitable schemes enforcing anisotropy is a non-trivial task since it requires the use of appropriate stencils that perfectly adapt to the local image structure. Many methods have been proposed; see e.g. the unifying framework in [40] and the references therein. While most stencils lead to L 2 -stable schemes, only a few of them allow to preserve nonnegativity and L ∞ stabilty [37, 25, 10] . A rather sophisticated representative among them is the stencil of Fehrenbach and Mirebeau [10] which relies on lattice basis reduction ideas.

In this work we consider a transport-diffusion PDE describing the physical phenomenon of osmosis for imaging applications. Compared to standard plain diffusion models, the model considered therein considers an additional drift term, making the process not symmetric (see [17] for the physical interpretation). For a regular domain Ω ⊂ R 2 , a given vector field d : Ω → R 2 and a given image f ∈ L ∞ (Ω; R), the isotropic osmosis model reads [38]      u t = ∆u − div(du) on Ω × (0, T ],

Differently from plain diffusion models, osmosis steady states are non-constant. In particular, if d is defined in terms of a given image v > 0 as d := ∇ log v, convergence to a rescaled version of v can be proven [38] . This is called the integrable or compatible case. In [38, 36] several imaging applications based on (2) or a slight modification thereof are studied. One of them is the shadow removal problem, see [38, Section 4.2] , as considered in this paper.

Shadow removal. The problem of shadow removal from a given image f : Ω → R + consists in removing the shadow appearing in f while preserving the image geometry and texture underneath. We will assume in the following constant shadows, i.e. where image intensity values inside and outside the shadow region are in relation with each other up to an (unknown) multiplicative constant. This problem is of great interest in computer vision as it often represents a preprocessing step in several segmentation, tracking and face recognition tasks where shadows are removed to avoid false detections/artefacts in the subsequent image processing. We refer the reader to [8] for a review on the existing models for shadow removal in images.

For our purposes, a mathematical formulation of the shadow removal problem can be obtained decomposing the image domain Ω as

where Ω out , Ω sb and Ω in are the unshadowed region, the shadow boundaries and the shadowed region of the image, respectively, see Figure 1 for an example.

Provided that a decomposition as in (3) is given, which for real images may a challenging problem on its own regard [8] , the osmosis model (2) can be easily adapted to solve the shadow removal problem by simply defining the vector field d in (2) in terms of a shadowed image f as d := ∇ log f on Ω in ∪ Ω out and d = 0 on Ω sb . The continuous Figure 1 : Decomposition of Ω as in (3) into (b) unshadowed region, (c) shadow boundaries and (d) shadowed region for a discrete shadowed image f . osmosis model adapted to shadow removal then reads

The evolution on the shadow boundary Ω sb can be interpreted as an inpainting step where information is propagated from Ω out to Ω in over Ω sb . Due to the action of the Laplace operator on Ω sb , image structures in Ω in ∪ Ω out are isotropically diffused on Ω sb , resulting in a shadowless, but blurred inpainting result on Ω sb . To overcome this a post-processing inpainting step is commonly applied, as for instance in Figure 2 .

(a) Isotropic osmosis [38] (b) Post-processing inpainting step [2] (c) Zoom of Figure 2a Figure 2b to remove the blurring artefacts due to Laplace inpainting on Ω sb in (4).

Note that in natural images several acquisition and/or compression artefacts may render the automatic segmentation of the shadow boundary very challenging. On the other hand, its accurate manual selection may be very tedious. In many practical examples, a rough selection of Ω sb is therefore performed manually by using a brush whose possibly large thickness may badly affect the result of the model (4) (see Figure  4 ) due to the Laplace blurring artefacts discussed above.

Vogel et al. [36] have presented a discrete osmosis theory and have proven that explicit and implicit finite difference discretisations satisfy its requirements. Different splitting schemes have been considered in [5, 27] and have been applied to imagery for cultural heritage conservation.

In this paper we extend the isotropic osmosis model (2) and its shadow removal application (4) to a model that features anisotropic diffusion in the flavour of (1) and adapts concepts from anisotropic diffusion inpainting [39, 13] . This will be implemented by incorporating local directionality depending on image orientations. We show that with this modification we can improve the solution of the shadow removal problem, in particular overcoming blurring artefacts in the region around the shadow boundary, without additional post-processing steps. In the integrable case, the resulting driftdiffusion PDE can be derived as the gradient flow of a suitable energy depending on local gradient information. To estimate such local directionality, we adapt the tensor voting framework proposed in [16] . For the numerical solution we combine a numerical time-stepping method based on exponential integrator techniques with the nonnegative space discretisation of Fehrenbach and Mirebeau [10] . Our model is validated on several synthetic and natural images affected by constant shadows. Results show good lightbalance properties and, compared to plain isotropic osmosis models, avoid the smoothing artefacts on the shadow boundary. An illustrative example of the performance of our model is reported in Figure 3 . [38, 36] and our anisotropic one to solve the shadow removal problem.

Organisation of the paper. In Section 2 we introduce the anisotropic osmosis model and study analytically its properties. Then, in Section 3 we study space and time discretisation schemes for the anisotropic model. Finally, in Section 4 we show the application of the anisotropic model to solve the shadow removal problem.

We present in this section a variation of the classical osmosis model (2) encoding local directional information of the image in the diffusion term, propagating geometric structures dominantly along locally preferred directions. For this reason we call our model anisotropic osmosis model in contrast to the model (2) which we refer to as isotropic osmosis model. In what follows we introduce the general form of our anisotropic osmosis model and state some properties of solutions that are inherited from the isotropic model. For specific choices of anisotropy we also show connections of the anisotropic osmosis model to anisotropic diffusion-based inpainting methods such as edge-enhancing anisotropic diffusion [39] . Out of these specific instances we derive our proposed anisotropic osmosisinpainting model for shadow removal.

Let us define an anisotropic osmosis energy as follows.

Definition 2.1 (Anisotropic osmosis energy). Let Ω ⊂ R 2 , u, v ∈ H 1 (Ω; R + ) be two positive images and let W : Ω → R 2 be a positive semi-definite symmetric matrix field. We define the anisotropic osmosis energy of u with respect to W and the reference image v as

We will also use the following alternative notation for E:

where e W := e, We .

Remark 2.2 (Isotropic case). If W is the identity matrix, then (5) corresponds to the isotropic osmosis energy considered in [38] .

Next we define an anisotropic osmosis evolution whose steady state minimises our anisotropic osmosis energy. 

Proof. For any test function ϕ ∈ C ∞ c (Ω), we compute the optimality condition holding for any critical point u of E. We get:

where we have applied the divergence theorem and the Neumann boundary conditions in (7) . Due to the positivity of v and since ϕ is compactly supported in Ω, we have that for any x ∈ Ω:

By definition of d = ∇v v , we note that the above corresponds to the following PDE: div (W (∇u − du)) = 0, which is the steady state of (7).

Similar to the isotropic osmosis PDE (2), the anisotropic model (7) enjoys some properties which makes it amenable for imaging applications.

Theorem 2.4. The solution u : Ω → R of the anisotropic osmosis model (7) satisfies the following properties:

(i) conservation of the average grey value:

(ii) preservation of non-negativity:

u(x, t) ≥ 0, for all x ∈ Ω and t > 0;

(iii) non-constant steady states: The steady state of (7) is given by

Proof. We follow [38] and prove statements (i)-(iii) in turn.

(i) Let µ u (t) := 1 |Ω| Ω u(x, t) dx be the average grey value of the image u at time t ≥ 0. Applying the divergence theorem and the homogeneous Neumann boundary conditions in (7) we obtain:

and the statement follows.

(ii) Assume that T > 0 is the smallest time such that min x,t u(x, t) = 0, and that this minimum is attained in some inner point µ ∈ Ω. Since we have ∇u(µ, T ) = u(µ, T ) = 0, we deduce:

that is in correspondence to the point (µ, T ), the anisotropic osmosis equation behaves like the diffusion equation

with non-constant positive semidefinite diffusivity matrix W. For such equations a generalisation of the weak minimum/maximum principle (which holds in its classical form only for positive definite, i.e. elliptic, operators) holds true, see, e.g. [30] . Therefore, for any t ≥ T the solution of the anisotropic model remains non-negative and since the solution has been further assumed to be strictly positive for all t < T , including t = 0 due to the positivity of f , we have that it stays actually non-negative for any t ≥ 0.

(iii) For any c ∈ R, the function w := cv endowed with the Neumann-type boundary conditions in (7) solves the steady state equation of the system anisotropic PDE, since there trivially holds:

Due to mass and non-negativity conservation, the process is then forced to a nonnegative steady state solution w of such a form. Furthermore, the constant c ∈ R can be easily found by noticing

Anisotropic diffusion inpainting with a diffusion tensor has been introduced in [39] and applied successfully for inpainting-based compression [13, 33] . It exploits the edgeenhancing anisotropic diffusion filter that has been proposed for image denoising [37] . In order to propagate structures from specified image regions into inpainting regions, one uses the differential operator div(D(∇u σ )∇u), where u σ denotes the convolution of u with a Gaussian of standard deviation σ. The diffusion tensor D has eigenvectors perpendicular and parallel to ∇u σ . Its corresponding eigenvalues are given by

with some contrast parameter λ > 0. Thus, the goal is to inpaint fully in the direction of an oriented structure, and to reduce the inpainting perpendicular to a structures, if its contrast is large. Processes of this type can inpaint edge-like structures even when the specified data are sparse and the gaps to be bridged are large [33] . However, they are not well-suited to shadow removal problems, since the shadow boundaries create unphysical edges. Therefore, we will have to modify these ideas such that the local structure directions become more robust w.r.t. shadow boundaries. To this end, we will consider and modify more refined structure descriptors such as tensor voting. This will be done next.

In this section we present a work-flow which is locally insensitive to the light jump produced by the shadow. This will serve us to force anisotropy on Ω sb along suitable directions.

A standard way to provide an estimate of the local structure orientation in an image u consists in computing the eigenvector e 1 associated to the leading eigenvalue λ 1 of the structure tensor J ρ (u) [37] associated to u. By fixing σ, ρ > 0 to be the pre-and post-smoothing parameters, we recall that the structure tensor J ρ (u) of an image u is defined as:

where u σ := K σ * u is a smoothed version of the image u and K σ , K ρ are Gaussian convolution kernels. Usually, σ and ρ are chosen so that σ ρ, where σ is associated to the noise scale and ρ integrates the orientation information. Denoting by λ 1 ≥ λ 2 the eigenvalues of J ρ (u) and by e 1 , e 2 the associated eigenvectors, a classification of the different structural image information in terms of the size of λ 1 and λ 2 can be made following, e.g., [18, 11, 20] . For any x ∈ Ω we thus have:

, then x is likely to belong to a homogeneous region;

, then x is likely to lie on an edge;

, then x is likely to be a corner point.

Tensor Voting has been originally introduced in [16] for extracting curves in noisy images by means of the grouping of local features consistent in a neighbourhood of the measurements. Such framework improves the robustness of structure tensor estimation in presence of noise and image artefacts [24] . Assuming that a generic 2-tensor B in R 2 has the following matrix representation:

where λ 1 , λ 2 are the eigenvalues associated to the eigenvectors e 1 and e 2 , respectively, we have that (12) can be equivalently rewritten as

We can now distinguish the two following quantities:

is called saliency or stickness. It provides an estimation of the confidence on the direction e 1 and it is also called orientation certainty or anisotropy measure;

• λ 2 is called ballness and it measures the size of the minor-axis of the anisotropy ellipse. Since it is in some sense a measure of how the estimation of the main estimated direction is contradicted, it is often called orientation uncertainty or junctionness.

The tensor voting operation consists in adding, at each iteration, the contribution of neighbourhood tensors for each point in the domain, resulting in an enhanced tensor field due to the presence of saliency parameter. In [12] the authors show that the complexity of the original approach may be time-consuming, even for small images. Thus, they propose an efficient computation of the tensor voting framework based on steerable filters theory, i.e. complex-valued convolutions.

For the shadow removal application we are considering, the estimation of the structure direction e 1 may be affected by the shadow edges, which do not correspond to the actual underlying image structures. The presence of such edges make the use of either the structure tensor or the tensor voting very challenging for estimating the local directions.

To circumvent this problem, we propose a modification of the tensor voting framework by means of interpreting shadow edges as bias in the estimation of the structure. From the given initial blurred image u σ , we firstly compute the local orientation θ of the gradient from the eigenvector e 2 . Secondly, we compute the saliency of the given shadowed image. Finally we apply the tensor voting framework with a modified saliency and orientation on Ω sb so as to mark the shadow boundaries as bias: there, the saliency is zeroed and the orientation is initialised at random. Details on the algorithm and results are presented in Section 4.2.1.

Let f : Ω → R + be a positive greyscale image with a constant shadow and let Ω be decomposed as in (3) . We propose the following structure-preserving osmosis model for solving the shadow removal problem:

on Ω,

Here we define the discontinuous vector field d and the discontinuous matrix field W as

where I denotes the 2 × 2 identity matrix, and e 1 and e 2 are the directions from our modified tensor voting applied to the initial image f . This means, an isotropic osmosis evolution is performed in Ω c sb , while an anisotropic inpainting is performed on Ω sb . In other words, classical osmosis balances image intensity in the shadowed region with respect to the unshadowed regions, while the nonlinear interpolation preserves structures and avoids blurring on the shadow boundary. The proposed model performs the osmosis and the inpainting step jointly and avoids any post-processing step.

We now discuss appropriate discretisations of the anisotropic osmosis model that are consistent with the continuous model (7) and preserve some of its properties as stated in Theorem 2.4. To do so, we consider a discrete rectangular image domain with M × N pixels. Let S := M N . The given positive initial image f is then defined as a vector in R S + . For a given grid step size h > 0, we denote by u = (u i,j ) i,j the approximation of the function u and with u i,j its approximated value in suitable discretisation nodes ((i − 1 2 )h, (j − 1 2 )h) with i = 1, . . . , M and j = 1, . . . , N . Similarly, for k ≥ 0 we denote by u k i,j the value of u i,j at the time node t k = kτ , where τ is the time step size. Also, for x ∈ Ω, we denote by λ 1 , λ 2 ∈ R S + the discretised eigenvalues λ 1 (x) and λ 2 (x), while by θ ∈ [0, 2π) S the discretised orientation θ(x).

A discrete theory for osmosis models has been established in [36] . Since it is also applicable to the anisotropic setting, we report here the general result [36, Proposition 1] and list in the following some discrete solvers fulfilling its assumptions. . For a given f ∈ R S + , consider the fully-discretised problem:

where the (unsymmetric) matrix P ∈ R S×S is an irreducible, non-negative matrix with strictly positive diagonal entries and unitary column sum . Then the following properties hold true:

(i) The evolution preserves positivity and the average grey value of f ;

(ii) The eigenvector of P associated to eigenvalue 1 is the unique steady state for k → ∞.

As shown in [36] for the isotropic osmosis model, standard explicit and implicit time finite difference schemes fit this framework, the former being subject to time step size restrictions, the latter being unconditionally stable. Furthermore, both schemes converge to the space discretisation of the elliptic steady state for every stable time step size. For the implicit scheme with a BiCGStab solver, Vogel et al. [36] report speed-ups of two orders of magnitude compared to the explicit method. In [5] a Peaceman-Rachford splitting is shown to satisfy Theorem 3.1 under a time step size restriction, while the additive operator splitting (AOS) considered in [27] fulfils Theorem 3.1 for all time step sizes. Splitting schemes such as the AOS method, however, do not converge to the space-discrete elliptic steady state solution unless the time step size goes to zero. In practice, keeping this time splitting error under control imposes bounds on the time step size.

Note that in contrast to the theory for fully discrete diffusion filters [37] , the discrete osmosis theory in Theorem 3.1 does not require a symmetric matrix. Since it does also not involve any isotropy assumption, it is basically also applicable to anisotropic osmosis processes, if suitable space discretisations are employed. This, however, is not straightforward and will be discussed in the sequel.

In what follows, we describe the discretisation of the weighting matrix W and the differential operators div and ∇ for a grey-scale image of height M and width N unrolled as u ∈ R S , with S = M N , that defines the spatial discretisation matrix A for the semi-discrete osmosis problem u (t) = Au(t), for t ∈ (0, T ],

where T > 0 is a positive final time.

While diffusion processes are stable in many aspects, e.g. in terms of decreasing L 2 norms and decreasing L ∞ norms, the stability of osmosis processes is restricted to essentially one key property: the preservation of nonnegativity. Thus, any suitable space discretisation for osmosis filters should guarantee that it is nonnegativity preserving. For the matrix A this implies that all off-diagonal elements must be nonnegative. While this is easily satisfied for standard discretisations of isotropic processes [36] , it becomes much more challenging for anisotropic approaches: The drift term is fairly unproblematic and can be handled e.g. with classical upwind discretisations. However, most discretisations of the diffusion term div(W∇u) are only stable in the L 2 norm [40] . Thus, they cannot guarantee preservation of nonnegativity. Nonnegativity-preserving discretisations can be found in [37, 25, 10] . In our paper, we use the stencil of Fehrenbach and Mirebeau [10] , since it is a fairly sophisticated nonnegativity-preserving method that has been reported to give good results. This anisotropic diffusion discretisation relies on lattice basis reduction ideas and is called AD-LBR. Let us sketch its underlying ideas.

In [10] , the authors tackle the minimisation of an anisotropic energy which in our notation reads:

where e W := e, We , for any e ∈ R d and W is symmetric and positive definite. The idea is to introduce a discretisation E h of the energy above on the discretised domain Ω h , h > 0 via a sum of weighted squared differences of u ∈ L 2 (Ω h ; R), i.e.:

where V (x) ⊂ Z d is the stencil and γ x (e) ≥ 0 are the associated weights. The key step is the linearisation of u(x + he) as u(x) + ∇u(x), he , which shows that for each x ∈ Ω h and smooth u one can write:

which turns out to be equivalent to require the condition W = e∈V (x) γ x (e)ee T . via the following Lemma [10, Lemma 1].

Lemma 3.2. Let e 0 , e 1 , e 2 ∈ R 2 be such that e 0 + e 1 + e 2 = 0 and |det(e 1 , e 2 )| = 1.

Then for any symmetric positive definite matrix W, there holds:

under the convention e 3+i := e i Actually, for any dimension d ≤ 3 and any symmetric positive definite d × d matrix M, there always exists a family (e i ) i∈I of vectors in R d such that e i , Me j ≤ 0 for any i = j. Such a family is called M-obtuse [9] . Thus, by taking (e 0 , e 1 , e 2 ) as a W-obtuse superbase of Z 2 (i.e. a basis of Z 2 with |det(e 0 , e 1 , e 2 )| = 1 such that e 0 + e 1 + e 2 = 0), a stencil V (x) := {e 0 , e 1 , e 2 , −e 0 , −e 1 , e 2 } can be used to write explicitly the coefficients γ x in (18) for 0 ≤ i ≤ 2 as done in [10, Equation (11)]:

The resulting stencil is shown to be independent of the choice of the superbase [10, Lemma 11] and it is orientated along the preferred diffusion direction given by W.

Also, the AD-LBR scheme is sparse, i.e. it has a limited support of 6 points for twodimensional images.

In the AD-LBR discretisation, an important role is played by the anisotropic ratio κ ∈ [1, ∞), which measures the geometrical shape of the ellipse associated to the eigendecomposition of the anisotropic matrix W: the closer κ is to 1, the more W is similar to the identity matrix I. The computation of the stencil has a logarithmic cost in the anisotropy ratio κ of the diffusion tensor, making AD-LBR appealing for applications. Also, by fixing a direction θ all over the domain Ω for W, the anisotropy ratio κ can be related to the eigenvalues of W; see [10, Equation 60 ]. For instance, if W has eigenvalues ε and 1 with 0 < ε 1, then ε = 1/κ 2 . For completeness, we report in Table 1 the AD-LBR stencil for different choices of ε and a fixed angle θ = 2π/3 that characterises the second eigenvector e 2 = (cos θ, sin θ) of W; cf. also [10, Table 1 , Table 2 ]. Table 1 : AD-LBR stencil for the discretisation of div(W∇( · )): Different choices of ε are presented with fixed angle θ = 2π/3 for the first eigenvector. In bold we denote the (i, j) entry. As expected, we observe that in the case of strong anisotropy, the stencil becomes aligned in θ direction. Remark 3.3. By construction, the AD-LBR space discretisation matrix A has 0 column sum and non-negative off-diagonal entries. Moreover, our numerical experiments give strong evidence that A is also irreducible †. However, a formal proof of the irreducibility of A is left for future research.

In this section we see how it is possible to solve the dynamical system (17), through a highly accurate approximation of the exact solution

First of all, we notice that it is not necessary to compute the large and dense matrix exp(T A) explicitly. It is sufficient to compute only its action on the initial solution f . Polynomial methods (see, for instance, [31, 6, 1] , approximate the action of the exponential by a polynomial of a certain degree applied to the initial vector. They do not require to solve a linear system of equations: usually they scale the matrix and approximate the solution by an iterative procedure like u k+1 = p m k (α k τ A)u k , k = 0, 1, . . . , K − 1, u 0 = f †We applied the Tarjan's algorithm finding the strongly connected components of a directed graph [35] . Code is freely available at MATLAB central: mathworks.com/matlabcentral/fileexchange/50707. where p m k is a polynomial of degree m k which approximates the exponential function and

After the last iteration, u K ≈ u(T ). Such an iterative scheme is usually needed in order to achieve a sufficiently accurate result. For Krylov methods (like [31] ), it is also necessary in order to keep the computational cost as low as possible, since the cost to produce p m k requires O(m 2 k ) scalar products with vectors of the size of f . Among the polynomial methods, the truncated Taylor series [1] and the interpolation in the Newton form at Leja points [6] are able to bound the relative backward error. This means that they construct an approximation in K iterations with time step size τ = T /K:

The tolerance tol can be chosen as small as desidered. The typical value for double precision arithmetic is 2 −53 and in this sense the time integration is said "exact". The polynomial p m (z) is either the truncated Taylor series of e z about a point which depends on the spectrum of the matrix or the interpolation polynomial of e z at real Leja points on an interval related to the spectrum of the matrix.

The choice of the parameters K and m can be done by simply considering the 1norm of T A, or estimates of (T A) q 1/q 1 for small values of q, or, only for interpolation at Leja points, estimates on the -pseudospectra of T A. From the implementation point of view, both algorithms simply require one matrix-vector product and one vector update at each degree elevation, and therefore the cost is O(m). When the spectrum of T A has a skinny shape, either horizontal or vertical, usually interpolation at Leja points performs better (see [6] ).

Since this method can be configured to be exact in time up to machine precision, it is basically possible to reach the final time T in a single time step of size τ = T . However, we prefer to use multiple time steps, since the steady state T is in general not known a priori. On the other hand, there is no restriction on the time step τ and it is therefore possible to implement any desired strategy for steady state detection (based, for instance, on the comparison of the solutions at two successive time steps). In particular, also variable step size τ k implementations are possible without any restriction given by the stability or the computational cost.

The approximation of the action of the matrix exponential to a vector can be used also in the so called exponential integrators (we refer to the survey paper [19] ) when solving general (non-linear) stiff ordinary differential equations.

A natural question in this context is whether the matrix P := exp(τ A) satisfies the properties of Theorem 3.1 for suitable matrices A. We answer this question with the following proposition, for which a simple lemma is required. Proof. For every i, j = 1, . . . , N , we write each element b i,j in terms of the elements of C and D. For every column j we have:

Proposition 3.5. Let A be an irreducible matrix, with column sum 0 and non-negative off-diagonal entries. Then the (non-symmetric) matrix P := exp(τ A) is an irreducible positive matrix with column sum 1.

Proof. We firstly show that P has column sum 1. To this end, we use the expansion

By hypothesis, A has column sum 0. Then the column sum of A k is 0 for every k ≥ 1 by Lemma 3.4. Thus, the matrix P = exp(τ A) has column sum 1.

In order to show that all elements of P are positive, we rewrite A as

where D = diag(A), I is the identity matrix, N = A − diag(A) is a non-negative matrix and α := min i a i,i − γ, for any γ > 0. We have that P can be expressed as the product

where τ (D − αI + N) is a non-negative matrix with positive diagonal and which is also irreducible. In fact, being that false, there would exist an extra-diagonal element in position (i, j) such that (D − αI + N) k i,j = 0 for all k > 0. But then, (A) k ij = (αI + (D − αI + N)) k i,j = 0 for all k > 0, meaning that A is reducible, which is false by hypothesis. So, for any pair (i, j) there exists a k such that (τ (D − αI + N)) k i,j > 0. Since all the powers of τ (D − αI + N) appear in the series which defines exp(τ (D − αI + N)), we conclude that all its entries are positive. Finally, exp(τ A) is obtained by scaling with positive scalars the rows of exp(τ (D − αI + N)). Therefore, exp(τ A) is a positive (and thus irreducible) matrix.

Recalling Remark 3.3, Proposition 3.5 shows that the AD-LBR space discretisation in combination with our "exact" polynomial time discretisation leads to a fully discrete anisotropic osmosis scheme that satisfies all assumptions of the discrete theory (subject to the missing formal irreducibility proof).

In this section we present several numerical examples showing the application of the isotropic and anisotropic osmosis model to solve the shadow removal problem in synthetic and real-world images. We apply the anisotropic model (13) to images affected by almost constant shadows, in order to perform jointly the shadow removal and the inpainting procedure on the shadow edges. Since in general the ground truth is not available for this problem, the quality of the reconstruction of the anisotropic model in comparison the isotropic approach is assessed by visual inspection.

On the thickness of the shadow boundary. In the following numerical experiments we will assume for simplicity that a rough segmentation of the shadow boundary is provided beforehand. For natural real images, this may be a quite challenging task since, due to the possible presence of noise, blur and/or compression artefacts, such region may be not sharp and presents a blurred transition zone from the outside to the inside area of the shadow. As a consequence, standard segmentation methods based, for instance, on edge detection methods may fail. In fact, the task of shadow segmentation has been addressed on its own regard in previous works where brightness-based [4] or clustering [8] methods have been applied. However, in many practical situations, the shadow is often roughly detected manually by the user using a brush including pixels both from the inside and the outside of the shadowed area. Such region corresponds to Ω sb in (14) .

In our experiments we have observed that if this region is chosen to be too small (i.e. smaller than the whole transition area between non-shadowed and shadowed area), then the shadow removal result is fairly poor. On the contrary, in general, selecting a thicker shadow boundary produces better results. Thick boundaries favour the use of the anisotropic model (13)-(14) over the isotropic one (4): as we seen above, the action of the homogeneous diffusion smoothing on a large region Ω sb is not able to preserve the underlying image structures. Figure 4 illustrates these findings. We compare the result obtained applying the isotropic model (4) on a real image where the thickness of the shadow boundary is chosen differently. We clearly observe that a thicker Ω sb corresponds to a better removal of the shadow for the same large final time T .

We will now show the numerical results obtained when the anisotropic model (13)- (14) was applied to a variety of both synthetic and natural images.

Pseudocode. Algorithm 1 describes the main steps required to solve the joint anisotropic osmosis problem after the estimation of local structures in the given image. For the time integration we used the expleja solver ‡, which is the companion software of [6] for computing the action of the exponential matrix on a vector. ‡Freely available at: bitbucket.org/expleja/expleja 

Firts we apply the anisotropic osmosis model to noise-free synthetic images, where the direction of the gradient z is known a priori. The purpose of this synthetic experiment is to check if our approach is able to effectively remove constant shadows. In Figure 5 we show the results obtained for an image with parallel greyscale stripes with θ = 65 • orientation and for colour concentric circles, whose θ is chosen to be as the angle drawn with tangent to the circumferences. Both images are corrupted with an almost constant shadow but a transition zone on the shadow borders. We compare the anisotropic model described in Section 3.2 with the isotropic osmosis model (4) , which results in an homogenous diffusion inpainting at the shadow boundary. The time discretisation is performed as described in Section 3.3. In our visual comparison of both methods, we choose the final time T = 10000 and the time-step τ = 100, and we proceeded using Algorithm 1. In this experiment we also provide a visual representation of the local orientation angle θ inside the shadow boundary Ω sb . It becomes obvious that the anisotropic shadow removal method shows clear advantages at the shadow boundaries, since it does not suffer from blurring artefacts. 

In order to apply the anisotropic model to natural real-world images, the estimation of the discrete local orientation θ(x) becomes crucial. Therefore, we show in Section 4. For the tensor voting, we used the MATLAB implementation of [12] from the companion software § of [23] , a literature review on tensor voting. Also, since tensor voting depends on local neighbourhoods, we use a multi-resolution strategy as is described in Algorithm 2. // zeroing and randomizing the data on the shadow edges Ω sb saliency loc = saliency loc.*mask; orientation loc = orientation loc.*mask + 2π.* rand(size(f (:, :, c))).*(1-mask);

foreach k in scales do

[ saliency, ballness, orientation ] = vote ( saliency loc, orientation loc, s k ); λ 1 = saliency +ballness; λ 2 = ballness; e 1 = ( cos(orientation), sin(orientation)); e 2 = (-sin(orientation), cos(orientation)); TVF = TVF + eigen to tensor ( e 1 , e 2 , λ 1 , λ 2 ); end end [ e 1 , e 2 , λ 1 , λ 2 ] = tensor to eigen ( TVF ); θ = xy2ij (e 2 ); // return the local orientation in (i, j) coordinates return θ

In Figure 6 we compare the local direction estimation by means of the structure §Freely available at MATLAB central: mathworks.com/matlabcentral/fileexchange/47398 tensor and the tensor voting approach applied on the shadowed image in Figure 6a . We plot the magnitude of the leading eigenvalue of the structure tensor in Figure 6b and of the one computed for the tensor voting Algorithm 2 in Figure 6c . These images clearly show that the estimation via tensor voting is less sensitive to the false edges introduced by the shadow boundaries. This is reflected in the plot of the leading directions, too:

The directions computed via the structure tensor in Figure 6d are visibly affected by the light jumps while the ones computed with tensor voting in Figure 6e can still connect the structures from outside to inside the shadow with no additional edges.

(a) Input f (b) STF, λ 1 (c) TVF, λ 1 (d) STF, e 1 (e) TVF, e 1 Figure 6 : Comparison: structure tensor framework (STF) with (σ, ρ) = (0.5, 4) versus tensor voting framework (TVF) with multi-scales (5, 10, 15) and σ = 0.5. We plot the main direction e 1 and its associated eigenvalue λ 1 .

Results on real images. We now combine the proposed tensor voting framework with the anisotropic osmosis model to remove constant shadows by means of an anisotropic drift-diffusion model. In the following experiments we use the final time T = 100000 and the time step size τ = 1000. In practice the shadow removal is accomplished almost completely already for t T . However, for a better approximation of the steady state, we use the final time T for comparison.

In Figure 7 we show a zoom of the results from Figure 3 ¶, presented as motivation for this work. Note that in this case the shadow is artificially added as a multiplicative rescaling factor c ∈ (0, 1). The estimation of the direction on the shadow boundary is computed via the proposed Algorithm 2.

(a) Shadowed f Orientation angle θ Isotropic Anisotropic Figure 7 : Zoom of the results for the shadowed image in Figure 3 . ¶Figure 3a (zoomed in Figure 7a ) courtesy of R. D. Kongskov. In Figure 8 we apply the isotropic and the anisotropic model to a several shadowed images affected by natural constant shadows . Zoomed details can be found in Figure 9 .

Here the complexity of local image structures makes the application of the anisotropic osmosis model harder. Thus, we use the tensor voting algorithm 2 to estimate the local directions connecting the structures from inside to the outside of the shadow region. Also in these real-world scenarios, we observe the superiority of the anisotropic approach: structures are interpolated reliably across the light jump introduced by the shadow. 

We note that our anisotropic approach removes shadows effectively both for synthetic and real-world images. Structures are propagated correctly over the shadow boundary.

Since the AD-LBR stencil requires a positive definite matrix W as input, we need to choose ε > 0, which may lead to some over-smoothing in the orthogonal direction. Also, we noticed some over-smoothing effect in Ω out and Ω in for the AD-LBR scheme, e.g. in Figure 9f , whose understanding is a matter of future research.

In terms of efficiency, the computational time needed for computing the action of the exponential matrix onto a vector, that is the product exp(τ A)u 0 , largely depends on the choice of the time-step τ , the final time T and the size of the images. Although it is possible to directly choose τ = T and proceed in a single step, we prefer multiple steps. In this way, by comparing two successive solutions u k+1 and u k , it is possible to detect whether the evolution is sufficiently close to its steady state. Our solvers that are exact in time offer additional advantages over classical inexact methods such as explicit and implicit schemes when one is also interested in good approximations of intermediate results.

In this work, we have generalised isotropic osmosis filtering introduced in [36, 38] to the anisotropic setting. This was achieved by introducing a weight matrix whose directional information was extracted from a modified tensor voting approach [24] . When applied to the shadow removal problem, the anisotropic model acts as an inpainting interpolator on the shadow boundary. It is close in spirit to inpainting with edge-enhancing anisotropic diffusion [39, 13, 33] .

From a numerical point of view, we have combined the nonnegativity preserving anisotropic diffusion stencil of Fehrenbach and Mirebeau [10] with techniques based on exponential integration [6] : We argued that the fully discrete model satisfies the discrete properties studied in [36] in a general setting for osmosis. We tested the proposed model for synthetic and real examples, showing that the generalised model acts as a combined osmosis-inpainting model for shadow removal problems, thus avoiding any undesirable post-processing inpainting step.

Future work will address the investigation on the applicability of the proposed anisotropic model to more general imaging applications.

Computing the action of the matrix exponential with an application to exponential integrators

A variational framework for exemplar-based image inpainting

Mathematical problems in image processing: Partial Differential Equations and the Calculus of Variations

Shadow removal from a real picture

Sketches & Applications, SIGGRAPH '03

Alternating direction implicit (ADI) schemes for a PDE-based image osmosis model

The Leja method revisited: Backward error analysis for the matrix exponential

Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods

Fast shadow removal using adaptive multi-scale illumination transfer

Low-dimensional lattices. vi. Voronoi reduction of three-dimensional lattices

Sparse non-negative stencils for anisotropic diffusion

A Feature Based Correspondence Algorithm for Image Matching

An efficient method for tensor voting using steerable filters

Image compression with anisotropic diffusion

Anisotropic total variation filtering

A review of PDE models in image processing and image analysis

Inferring global perceptual contours from local features

Novel schemes for hyperbolic pdes using osmosis filters from visual computing

A combined corner and edge detector

Exponential integrators

Analyzing oriented patterns

Directional total generalized variation regularization for impulse noise removal

Directional total generalized variation regularization

Perceptual grouping by tensor voting: a comparative survey of recent approaches

Adaptation of tensor voting to image structure estimation

Consistent positive directional splitting of anisotropic diffusion

An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences

Digital cultural heritage imaging via osmosis filtering

Higher order total directional variation. Part I: Imaging applications (forthcoming)

Higher order total directional variation

Maximum Principles in Differential Equations

Analysis of some Krylov subspace approximations to the matrix exponential operator

Geometric Partial Differential Equations and Image Analysis

Understanding, optimising, and extending data compression with anisotropic diffusion

Partial Differential Equation Methods for Image Inpainting. Cambridge Monographs on Applied and Computational Mathematics

Depth-first search and linear graph algorithms

Scale Space and Variational Methods in Computer Vision

Anisotropic Diffusion in Image Processing

Linear osmosis models for visual computing

Tensor field interpolation with PDEs

L2-stable nonstandard finite differences for anisotropic diffusion

The authors thank J.M. Mirebeau for the insightful discussions and comments on the adaptation of the code from [10] to our problem and the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program "Variational Methods and Effective Algorithms for Imaging and Vision" when work on this paper was undertaken. This work was supported by EPSRC grant number EP/K032208/1.