# Fast and stable schemes for non-linear osmosis filtering

L.Calatroni111Universite Côte d’Azur, CNRS, Inria, Laboratoire I3S, Sophia-
Antipolis, France S.Morigi222Department of Mathematics, University of Bologna,
Bologna, Italy. S.Parisotto333Department of Applied Mathematics and
Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge, CB3
0WA, UK G.A.Recupero444Department of Mathematics, University of Bologna,
Bologna, Italy.

###### Abstract

We consider a non-linear variant of the transport-diffusion osmosis model for
solving a variety of imaging problems such as shadow/soft-light removal and
compact data representation. The non-linear behaviour is encoded in terms of a
general scalar function $g$ with suitable properties, which allows to balance
the diffusion intensity on the different regions of the image while preventing
smoothing artefacts. For the proposed model, conservation properties
(intensity and non-negativity) are proved and a variational interpretation is
showed for specific choices of $g$. Upon suitable spatial discretisation, both
an explicit and a semi-implicit iterative scheme are considered, for which
convergence restrictions and unconditional stability are proved, respectively.
To validate the proposed modelling and the computational speed of the
numerical schemes considered, we report several results and comparisons for
the problem of shadow/light-spot removal and compact data representation,
showing that artefact-free and computationally efficient results are obtained
in comparison to standard linear and anisotropic models, and state-of-the art
approaches.

###### keywords:

Osmosis filtering , Non-linear PDEs , Semi-implicit finite-difference schemes
, Shadow removal , Light-spot removal , Compact data representation

††journal: Computers & Mathematics with Applications

##  1 Introduction

The use of Partial Differential Equations (PDEs) in the context of image
processing and computer vision is nowadays pretty well-established. Linear and
non-linear diffusion PDEs have been indeed considered over the last decades
for several several image processing problems, such as denoising,
deconvolution, inpainting and segmentation, see, e.g., [1, 2, 3, 4]. For these
tasks, standard approaches resorting to the use of parabolic (anisotropic)
non-linear diffusion PDEs of second or higher-order have been used, and their
interpretation as gradient flows of suitable non-smooth regularisation
functionals such as the edge-preserving Total Variation (TV) semi-norm have
been drawn. In this work, we concentrate on a non-linear PDE where a non-
linear diffusion term is combined with a further transport term, which finds
its motivations into the modelling of the physical process of _osmosis_.

In biology, such term indicates a non-symmetric transport process where
molecules pass through a semipermeable membrane in such a way that at its
steady state, the liquid concentrations on both sides of the membrane differ.
Osmosis is, for instance, the primary mechanism for transporting water inside
and outside living cells, and it stands at the very basis of many medicine and
engineering processes. Due to its non-trivial (i.e. non-flat) steady states,
such process can be seen indeed as the non-symmetric counterpart of standard
diffusion processes in the sense that, during evolution, the probability of
moving from inside to outside the cell through the membrane is not equal to
the probability of performing the reverse process [5].

The corresponding interpretation of the osmosis process in the context of
image analysis involves the study of how intensity values at each pixel are
propagated within the image domain. Previous works have considered osmosis
filtering in the context of imaging, either in linear isotropic [6] or
anisotropic [7] form, proposing also some efficient numerical schemes based on
standard finite-difference discretisation [8] and/or operator splitting [9].
For all these models, the combination of diffusion with transport is made
thanks to the use of a drift term defined in terms of a vector field which,
intuitively, models the non-symmetric nature of the process. Interestingly,
such PDEs are all in divergence form, thus allowing for provable conservation
of average intensity of the initial image as well as non-negativity.
Furthermore, upon a particular choice of the drift term, an explicit
expression of the non-flat steady state can be found. When applied to imaging
tasks aiming at balancing intensity information (i.e. for shadow/light-spot
removal problems and/or for image fusion), the advantage of using an osmosis
approach is that it can integrate mismatching image data in a seamless way
(see also [10] for a variational osmosis-based model for image fusion showing
the advantage with-respect-to state of the art methods). Recently, a more
sophisticated approach has been proposed in [11] where a non-local version of
the osmosis operator along with an anisotropic regularisation term have been
integrated in a variational formulation. This model achieves good performance
in shadow removal, but at the price of a high computational cost.

Due to its differential structure, osmosis filtering shows similarities to
gradient domain methods introduced in computer graphics for a wide number of
surface editing and reconstruction applications, such as the famous Poisson
editing, see e.g., [12, 13]. Compared to these approaches, image osmosis
offers the additional advantage of being invariant under multiplicative
illumination changes and provide very good results in shadow removal and image
cloning. However, due to the underlying linear smoothing enforced by the
diffusion terms, standard osmosis approaches typically suffer from some
reconstruction artefacts (typically, oversmoothing) which prevent their
practical use. For reference, we further refer the reader to [14, 15, 16, 17,
18] where alternative approaches to shadow removal are used, based, for
instance, of suitable mapping into appropriate colour spaces and/or learning
approaches.

#### Contribution

In this paper, we present a non-linear version of the image osmosis model,
which avoids standard smoothing artefacts of linear models thanks to the use
of a non-linear diffusivity function promoting isotropic edge-stopping
diffusion. Analytically, we prove that the proposed evolution model maintains
the same conservation properties of the original linear model, and that it can
be interpreted as the gradient flow of a suitable non-smooth energy for
specific choices of the drift term. We then show that analogous properties
hold also at a discrete level, whenever appropriate finite-difference
discretisation stencils and explicit and semi-implicit schemes are used. For
those schemes, we prove conditional and unconditional stability, respectively.
Such property is particularly interesting as it allows for the fast
computation of the desired steady-state solution. Thanks to such stability
property, the proposed model is fully automatic as no hyperparameters need to
be tuned in it. The proposed model and algorithms are validated on exemplar
image processing tasks, such as shadow removal, light removal and compact data
representation, and several comparisons with other state-of-the-art as well as
alternative osmosis filtering models are given.

#### Structure of the paper

In Section 2, a short review of existing osmosis filtering models and of their
use in applications is given. In Section 3, the proposed continuous non-linear
osmosis model is proposed and its conservation and variational properties are
proved. In Section 4, spatial and temporal discretisation schemes are studied,
in particular from the point of view of consistency with the aforementioned
properties in the continuous setting and in terms of stability and
convergence. In Section 5, several results for three different imaging
applications (shadow removal, light removal and compact data representation)
are showed, confirming that the use of the proposed non-linear model improves
upon linear and state-of-the art methods. Some conclusions are reported in
Section 6.

#### Notation

We will denote the set of non-negative and positive real numbers by
${\mathbb{R}}_{+}$ and ${\mathbb{R}}_{+}^{*}$, respectively. In the following
Sections 2 and 3 we will use the bold notation to denote vectorial fields in
${\mathbb{R}}^{N}$. To avoid unnecessary heavy notation in the discretised
setting introduced starting from Section 4, we will use standard (unbold)
notation to denote both vectors in ${\mathbb{R}}^{N}$ and matrices in
${\mathbb{R}}^{N\times N}$.

##  2 Linear osmosis filtering

The original image osmosis model firstly introduced in [6] in the context of
imaging is a parabolic, linear and isotropic transport-diffusion PDE defined
on a rectangular domain $\Omega\subset{\mathbb{R}}^{2}$ for an initial
positive image $f:\Omega\rightarrow{\mathbb{R}}_{+}^{*}$ (the extension to the
vectorial case is immediate). Its transport component is defined in terms of a
drift vector field ${\textbf{d}}:\Omega\rightarrow{\mathbb{R}}^{2}$ whose
expression will be made more precise in the following. Assuming homogeneous
Neumann boundary conditions, the model reads:

$\displaystyle\begin{cases}\partial_{t}u={\textnormal{div}}\left(\,\nabla
u-{\textbf{d}}u\right)\quad&\text{on }\Omega\times(0,T],\\\ \left<\nabla
u-{\textbf{d}}u,{\textbf{n}}\right>=0&\text{on }\partial\Omega\times(0,T];\\\
u(x,0)=f(x)&\text{on }\Omega,\end{cases}$ (1)

where $T\in(0,+\infty)$ is a final time, and n is the outer normal vector to
the image boundary $\partial\Omega$. In [6], it has been proved that the
evolution process (1) preserves the average mass (grey value) of the image as
well as its non-negativity. Moreover, if the drift d is defined in a
_canonical_ form, i.e. in terms of a given reference image
$v:\Omega\rightarrow{\mathbb{R}}_{+}^{*}$ as:

${\textbf{d}}:=\nabla\log v,$ (2)

then it can be showed that the steady-state solution
$w:\Omega\to{\mathbb{R}}_{+}$ is nothing but a multiplicative rescaling of
$v$, that is $w=cv$, with $c>0$ [6, Proposition Ic)], see Figure 1 for a
visual illustration. It can be shown that analogous properties hold upon
suitable finite-difference discretisation [8].

(a) $f$

(b) $v$

(c) d

(d) $w$

Figure 1: Evolution of the osmosis model (1). From left to right: (constant)
initial image $f$, reference image $v$, vector representation of drift d (with
zoom) and steady state $w$.

An anisotropic linear version of (1), with analogous conservation and
convergence properties, has been developed in [7] in order to favour
directional linear diffusion. Similarly as in [1], this has been done by
encoding local directional information of the image in a positive semi-
definite symmetric tensor W. The corresponding anisotropic model takes the
form:

$\displaystyle\begin{cases}\partial_{t}u={\textnormal{div}}\left({\textbf{W}}\,(\nabla
u-{\textbf{d}}u)\right)\quad&\text{on }\Omega\times(0,T],\\\
\left<{\textbf{W}}(\nabla u-{\textbf{d}}u),n\right>=0&\text{on
}\partial\Omega\times(0,T];\\\ u(x,0)=f(x)&\text{on }\Omega.\end{cases}$ (3)

For both models (1) and (3), there is of course no clear interest in observing
convergence to a rescaled version $w$ of an image which is need in advance to
define the model. However, upon slight modifications of the drift term d,
osmosis models can be applied to many different problems such as shadow
removal, compact data representation and image cloning, as shown in [6]. In
the following, we review how these tasks can be tackled by osmosis.

###  2.1 Applications of the osmosis filtering

Osmosis models in the form (1) and (3) can be easily modified for several
imaging applications. To do so, a suitable decomposition of the image domain
$\Omega$ in terms of a disjoint union of subsets is required. Let thus
$\Omega=\Omega_{in}\sqcup\Omega_{b}\sqcup\Omega_{out}$ be a given partition of
$\Omega$, where $\Omega_{b}$ denotes a possibly expanded boundary region
separating $\Omega_{in}$ from $\Omega_{out}$. Given the reference image $v$, a
simple modification of the drift term d consists in forcing it to vanish in
correspondence of all points of $\Omega_{b}$ so that an artificial set of
discontinuities is created. In mathematical terms, this means choosing for all
$x\in\Omega$:

${\textbf{d}}(x):=\left(\nabla(\log
v)~{}\chi_{\Omega_{b}^{C}}\right)(x)=\begin{cases}\frac{\nabla
v(x)}{v(x)}\quad&\text{if }x\in\Omega_{in}\cup\Omega_{out},\\\ 0&\text{if
}x\in\Omega_{b},\end{cases},$ (4)

where for $S\subset\Omega$, $\chi_{S}:\Omega\to\left\\{0,1\right\\}$ denotes
the characteristic function of the subset $S$. Practically, the role of the
characteristic function is played upon discretisation by a mask matrix
$M\in\left\\{0,1\right\\}^{|\Omega|}$ whose null pixels identify the region
$\Omega_{b}$. As d is set to zero on $\Omega_{b}$, the transport term vanishes
therein, hence equations (1) and (3) induce linear isotropic/anisotropic
diffusion, respectively. Note that the preliminary detection of $\Omega_{b}$
can be obtained, for instance, by standard image segmentation techniques, such
as, e.g., edge detection [3].

#### Shadow removal

For shadow removal problems the given image $f$ is characterised by the
presence of a constant (i.e. with no space-variant intensity) shadowed region
$S\subset\Omega$, see Fig. 2(a). The mask $M$ is thus associated here to a
preliminary detection of the shadow boundary $\Omega_{b}$ (for instance by
using [16, 19]), see Fig.2(b). By choosing as a reference image $v=f$ in (4),
we can thus evolve model (1) as done, e.g., in [6, 8, 9], to obtain the result
in Fig.2(c).

Due to its intrinsic transport-diffusion properties, we observe that the
linear osmosis filtering does correctly perform shadow removal by balancing
the intensity between the shadowed/un-shadowed regions of the image. However,
due to the choice (4), the drift term d vanishes on the $\Omega_{b}$, hence,
pure Laplace diffusion is enforced therein. By considering a non-linear
version of (1) as described in the following sections, we will improve upon
such drawback and obtain the more accurate result reported in Fig.2(d).
Furthermore, in order to address the more realistic situations of shadows with
blurred (i.e. not sharp) boundaries as it is observed in natural images, we
will further show how the proposed model can be applied also in case of soft-
edged shadows, see Section 5.1.

(a) $f=v$

(b) mask $M$

(c) $w$ (linear)

(d) $w$ (non-linear)

Figure 2: Shadow removal using linear (1) and non-linear (proposed) osmosis
filtering.

#### Spot-light removal

A problem analogous to shadow removal is the spot-light removal problem where,
similarly, the task is to remove a light spot from a given image $f$ (Fig.
3(a)). Given the segmentation of the constant spot light boundary (Fig. 3(b)),
the use of osmosis removes the spot by balancing the intensity between the two
regions in the image. The results obtained by applying standard linear and the
proposed non-linear osmosis filtering are shown in Fig. 3(c)-3(d).

(a) $f=v$

(b) mask $M$

(c) $w$ (linear)

(d) $w$ (non-linear)

Figure 3: Spot-light removal using linear and non-linear (proposed) osmosis
filtering (1), with zoom.

#### Compact image representation

The task of compact image (or, generally, _data_) _representation_ consists in
representing a given image with as little information as possible using some
limited, but significant, image content (see, e.g. [20], [21]). Such problem
has been previously considered in [6] for the linear osmosis model: there, a
greyscale image is reconstructed by osmosis evolution by starting from the
intensity values on pre-detected edges which is encoded in the definition of
the drift-term d so that:

$\displaystyle{\textbf{d}}:=\left(\frac{\nabla
v}{v}\chi_{\Omega_{b}}\right)(x)=\begin{cases}\frac{\nabla
v(x)}{v(x)}\quad&\text{if }x\in\Omega_{b},\\\ 0&\text{if
}x\in\Omega\setminus\Omega_{b},\end{cases}$ (5)

where $\Omega_{b}$ denotes here the set of the edges. The task is thus a sort
of interpolation process where starting from $f$ (Fig. 4(a)) and using only
the data assigned on the mask $M$ (Fig. 4(b)) to define the drift term as in
(5), the osmosis evolution produces the result in Fig. 4(c). Fig. 4(d) is the
ground-truth image, reported here for comparison.

(a) $f=$ const

(b) mask $M$

(c) $w$ (non-linear)

(d) Reference $v$

Figure 4: Compact data representation using the proposed non-linear osmosis
model.

##  3 Analysis of the continuous model

Given a rectangular image domain $\Omega\subset{\mathbb{R}}^{2}$ with boundary
$\partial\Omega$ and a finite time $T>0$, let
${\textbf{d}}:\Omega\to{\mathbb{R}}^{2}$ be a given drift vector field, $f\in
L^{\infty}(\Omega;{\mathbb{R}}_{+}^{*})$ a positive greyscale image (the
extension to the RGB case is straightforward) and
$g:\Omega\times(0,T]\to{\mathbb{R}}_{+}^{*}$ be a positive (non-linear)
diffusivity function driving the evolution process. The proposed non-linear
osmosis model is given by the following drift-diffusion PDEs:

$\partial_{t}u(x,t)={\textnormal{div}}\left(g(x,t)\,\left(\nabla
u(x,t)-{\textbf{d}}(x)u(x,t)\right)\right)\quad\text{for
}(x,t)\in\Omega\times(0,T],$ (6)

which we endow with homogenous Neumann boundary conditions and initial
condition

$\displaystyle\begin{cases}\left<\nabla
u-{\textbf{d}}u,{\textbf{n}}\right>=0&\text{on }\partial\Omega\times(0,T];\\\
u(x,0)=f(x)&\text{on }\Omega.\end{cases}$ (7)

The model (6) favors diffusion and transport according to the function
$g(\cdot)$. The following result shows that thanks to its divergence form,
model (6) enjoys standard conservation properties for any $g$.

######  Proposition 3.1 (Conservation properties).

Any solution $u:\Omega\times(0,T]\rightarrow{\mathbb{R}}$ of the non-linear
osmosis process (6) satisfies the following properties:

  1. 1.

The average grey value is preserved:

$\frac{1}{|\Omega|}\int_{\Omega}u(x,t)dx=\frac{1}{|\Omega|}\int_{\Omega}f(x)dx=:\mu(f)\qquad\forall
t>0;$ (8)

  2. 2.

The evolution preserves non-negativity:

$u(x,t)\geq 0\quad\forall x\in\Omega,\forall t>0.$ (9)

###### Proof.

1) Define, for all $t\in(0,T]$, $m(t):=(\int_{\Omega}u(x,t)dx)/|\Omega|$.
Then, by (6),

$\displaystyle\frac{dm}{dt}$
$\displaystyle=\frac{1}{|\Omega|}\int_{\Omega}\partial_{t}u(x,t)dx=\frac{1}{|\Omega|}\int_{\Omega}{\textnormal{div}}\left(g(x,t)(\nabla
u(x,t)-{\textbf{d}}(x)u(x,t))\right)dx=$
$\displaystyle=\frac{1}{|\Omega|}\int_{\partial\Omega}\left<g(s,t)(\nabla
u(s,t)-{\textbf{d}}(s)u(s,t)),{\textbf{n}}(s)\right>ds=0,\qquad\forall t>0$

thanks to the divergence theorem and by imposing Neumann boundary conditions.

2) Let now $\tau>0$ be the smallest time such that $\min_{x}u(x,t)=0$. Suppose
that this minimum is obtained at a $\xi\in\text{int}(\Omega)$. Then
$u(\xi,\tau)=0$ and $\nabla u(\xi,\tau)=0$. Computing the time derivative of
$u$ at point $(\xi,\tau)$ gives:

$\displaystyle\partial_{t}u(\xi,\tau)$
$\displaystyle=\left({\textnormal{div}}(g(\cdot,\cdot)(\nabla
u(\cdot,\cdot)-{\textbf{d}}(\cdot)u(\cdot,\cdot))\right)(\xi,\tau)=$
$\displaystyle=g(\xi,\tau)\Delta u(\xi,\tau)+\left<\nabla u(\xi,\tau),\nabla
g(\xi,\tau)\right>-{\textnormal{div}}({\textbf{d}}(\cdot)g(\cdot,\cdot)u(\cdot,\cdot))(\xi,\tau)=$
$\displaystyle=g(\xi,\tau)\Delta
u(\xi,\tau)+0-({\textnormal{div}}~{}{\textbf{d}}(\cdot))(\xi,\tau)\cdot(u(\xi,\tau)g(\xi,\tau))-\left<{\textbf{d}}(\cdot),\nabla(g(\cdot,\cdot)u(\cdot,\cdot))\right>(\xi,\tau)=$
$\displaystyle=g(\xi,\tau)\Delta
u(\xi,\tau)-0-\left<{\textbf{d}}(\cdot),\nabla
u(\cdot,\cdot)~{}g(\cdot,\cdot)+u(\cdot,\cdot)\nabla
g(\cdot,\cdot)\right>(\xi,\tau)=$ $\displaystyle=g(\xi,\tau)\Delta
u(\xi,\tau)-0=g(\xi,\tau)\Delta u(\xi,\tau),$

by standard properties of the divergence operator. Then, at $(\xi,\tau)$ the
evolution behaves as the diffusion equation $u_{t}=g\Delta u$, where $g$ is
positive. Hence, the operator $\mathcal{L}:=g\Delta$ is elliptic so the
standard minimum/maximum principle can be applied. This tells us that for any
$t\geq\tau$ the solution of the non-linear model remains non-negative, as
desired. ∎

We now propose a special choice of the diffusivity function
$g:\Omega\times(0,T]\to{\mathbb{R}}_{+}^{*}$ and set for all $x\in\Omega$ and
$t\in(0,T]$:

$g(x,t):=g(|{\textbf{s}}(x,t)|):=\frac{1}{|{\textbf{s}}(x,t)|},\quad\text{with
}\;{\textbf{s}}(x,t):=\nabla u(x,t)-{\textbf{d}}(x)u(x,t),$ (10)

which will be used to characterise the steady states of (6) as minimisers of a
suitable energy functional. To show that, we first need the following lemma.

######  Lemma 3.2.

Let $\Omega$ be an open rectangular domain in ${\mathbb{R}}^{2}$ and
$\mathbf{w}\in C^{1}(\Omega,{\mathbb{R}}^{2})\cap
C^{0}(\overline{\Omega},{\mathbb{R}}^{2})$. Then:

$\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,z\;ds-\int_{\Omega}{\textnormal{div}}(\mathbf{w})\,z\;dx=0\quad\forall\;z\in
C^{0}(\overline{\Omega})$ (11)

if and only if

$\begin{cases}\left<\mathbf{w},{\textbf{n}}\right>=0\quad
on\;\partial\Omega\\\ {\textnormal{div}}(\mathbf{w})=0\quad
on\;\Omega.\end{cases}$ (12)

###### Proof.

Eq. (12) implies eq. (11) obviously. For the inverse implication, consider the
two alternative cases: either ${\textnormal{div}}(\mathbf{w})(x)=0$ for all
$x\in\Omega$, or there exists a point $x_{0}\in\Omega$ s.t.
${\textnormal{div}}(\mathbf{w})(x_{0})\neq 0$. We prove by contradiction that
the second case is not possible. Suppose that
${\textnormal{div}}(\mathbf{w})(x_{0})>0$. Since $\mathbf{w}\in
C^{1}(\Omega;{\mathbb{R}}^{2})$, then there exists a neighbourhood
$B(x_{0},r)\subset\Omega$ s.t. ${\textnormal{div}}(\mathbf{w})(x)>0$ for all
$x\in B(x_{0},r)$. Now, let $z\in C^{0}(\overline{\Omega})$ be a function that
satisfies (11). Consider

$z^{*}:=z+\Phi(\chi_{B(x_{0},r)})\in C^{0}(\overline{\Omega})$

where $\Phi(\chi_{B(x_{0},r)})$ is a non-negative smooth regularization of
$\chi_{B(x_{0},r)}$, with support in $B(x_{0},r)$. Then, we have

$\displaystyle\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,z^{*}\;$
$\displaystyle ds-\int_{\Omega}{\textnormal{div}}(\mathbf{w})\,z^{*}\;dx=$
$\displaystyle=\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,z\;ds-\int_{\Omega}{\textnormal{div}}(\mathbf{w})\,z\;dx-\int_{\Omega}{\textnormal{div}}(\mathbf{w})\Phi(\chi_{B(x_{0},r)})\,dx=$
$\displaystyle=0-\int_{\Omega}{\textnormal{div}}(\mathbf{w})\Phi(\chi_{B(x_{0},r)})\,dx=$
$\displaystyle=-\int_{B(x_{0},r)}{\textnormal{div}}(\mathbf{w})\Phi(\chi_{B(x_{0},r)})\,dx<0$

because $supp(\Phi(\chi_{B(x_{0},r)}))\subset B(x_{0},r)$, $\Phi(x)\geq 0$ and
${\textnormal{div}}(\mathbf{w})>0$ in $B(x_{0},r)$. We have thus reached a
contradiction, as for $z^{*}\in C^{0}(\overline{\Omega})$ equation (11) is not
verified. Therefore, it must be ${\textnormal{div}}(\mathbf{w})(x)=0$ for all
$x\in\Omega$. We claim that this implies also that

$\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,z\;ds=0\quad\forall\;z\in
C^{0}(\overline{\Omega},{\mathbb{R}}).$ (13)

By contradiction, consider a point $x_{1}\in\partial\Omega$ s.t.
$\left<\mathbf{w},{\textbf{n}}\right>(x_{1})\neq 0$ (suppose positive). Since
$\mathbf{w}\in C^{0}(\overline{\Omega};{\mathbb{R}}^{2})$, there exists
$B(x_{1},r)$ s.t. $\left<\mathbf{w},{\textbf{n}}\right>(x)>0$ for all $x\in
B(x_{1},r)\cap\partial\Omega$. Let now $z\in C^{0}(\overline{\Omega})$ be a
function satisfying (13) and define
$z^{*}:=z+\tilde{\Phi}(\chi_{B(x_{1},r)})$, where $\tilde{\Phi}(\cdot)$ is as
a smooth regularisation of $\chi_{B(x_{1},r)}$ with support in $B(x_{1},r)$.
Then,

$\displaystyle\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,z^{*}\;ds$
$\displaystyle=\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,z\;ds+\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,\Phi(\chi_{B(x_{1},r)})\;ds=$
$\displaystyle=0+\int_{\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,\Phi(\chi_{B(x_{1},r)})\;ds=$
$\displaystyle=\int_{B(x_{1},r)\cap\partial\Omega}\left<\mathbf{w},{\textbf{n}}\right>\,\Phi(\chi_{B(x_{1},r)})\;ds>0$

because $supp(\Phi(\chi_{B(x_{1},r)}))\cap\partial\Omega\subset
B(x_{1},r)\cap\partial\Omega$, $\Phi(x)\geq 0$ and
$\left<\mathbf{w},{\textbf{n}}\right>(x)>0$ in $B(x_{1},r)\cap\partial\Omega$.

We thus built a function $z^{*}$ that does not satisfy (13), which is of
course a contradiction. Therefore, $\left<\mathbf{w},{\textbf{n}}\right>(x)=0$
for all $x\in\partial\Omega$. ∎

Thanks to this lemma, we are now able to interpret the non-linear PDE model
(6) as the gradient flow of a suitably defined energy functional $E(u)$ in the
special case when the diffusivity function $g$ is defined as in (10).

######  Proposition 3.3.

Let $v\in H^{1}(\Omega;{\mathbb{R}}_{+}^{*})$ be a given reference image,
${\textbf{d}}=\nabla(\log v)$ the canonical drift vector field associated to
$v$ and $g$ be defined as in (10). Then, a function
$u^{*}:\Omega\times(0,T]\to{\mathbb{R}}_{+}$ satisfies the following steady
state equation

${\textnormal{div}}\left(g(|\nabla u^{*}-{\textbf{d}}u^{*}|)\left(\nabla
u^{*}-{\textbf{d}}u^{*}\right)\right)=0$ (14)

with boundary condition $\left<\nabla
u^{*}-{\textbf{d}}u^{*},{\textbf{n}}\right>=0$, if and only if $u^{*}$ is a
stationary point of the energy functional

$E(u):=\int_{\Omega}\left|\nabla\left(\frac{u}{v}\right)(x)\right|dx.$ (15)

Moreover, $u^{*}$ is a multiplicative rescaling of $v$ where the rescaling
constant is the ratio between the average gray value of $f=u(\cdot,0)$ and
$v$, i.e.:

$u^{*}(x)=\frac{\mu(f)}{\mu(v)}v(x),\quad\forall x\in\Omega.$

###### Proof.

Conventionally, for $x\in\Omega$ and $t\in(0,T]$ such that $\nabla
u(x,t)-{\textbf{d}}(x)u(x,t)=0$, we set

$\frac{\nabla u(x,t)-{\textbf{d}}(x)u(x,t)}{|\nabla
u(x,t)-{\textbf{d}}(x)u(x,t)|}=0.$

Let us now consider any test function $z\in C^{\infty}_{c}(\Omega;R_{+}^{*})$.
We have

$\displaystyle\left(\frac{d}{d\eta}E(u+\eta z)\right)$
${}_{\big{|}{\eta=0}}=\left(\int_{\Omega}\frac{d}{d\eta}\left|\nabla\left(\frac{u+\eta
z}{v}\right)\right|dx\right)_{\big{|}{\eta=0}}=$
$\displaystyle=\left(\int_{\Omega}\left<\frac{\nabla(\frac{u+\eta
z}{v})}{|\nabla(\frac{u+\eta
z}{v})|},\nabla\left(\frac{z}{v}\right)\right>dx\right)_{\big{|}{\eta=0}}=$
$\displaystyle=\int_{\Omega}\left<\frac{\nabla(\frac{u}{v})}{|\nabla(\frac{u}{v})|},\nabla\left(\frac{z}{v}\right)\right>dx=$
$\displaystyle=\int_{\Omega}\left<\frac{\nabla u-{\textbf{d}}u}{|\nabla
u-{\textbf{d}}u|},\nabla\left(\frac{z}{v}\right)\right>dx=$
$\displaystyle=\int_{\Omega}\left[{\textnormal{div}}\left(\frac{\nabla
u-{\textbf{d}}u}{|\nabla
u-{\textbf{d}}u|}\cdot\frac{z}{v}\right)-\frac{z}{v}\,{\textnormal{div}}\left(\frac{\nabla
u-{\textbf{d}}u}{|\nabla u-{\textbf{d}}u|}\right)\right]dx=$
$\displaystyle=\int_{\partial\Omega}\left<\frac{\nabla
u-{\textbf{d}}u}{|\nabla
u-{\textbf{d}}u|}\cdot\frac{z}{v},{\textbf{n}}\right>ds-\int_{\Omega}\frac{z}{v}\,{\textnormal{div}}\left(\frac{\nabla
u-{\textbf{d}}u}{|\nabla u-{\textbf{d}}u|}\right)dx$ (16)

The function $u^{*}$ satisfying (14) and the homogeneous Neumann boundary
condition is thus a stationary point of the energy $E(u)$, i.e. (3) vanishes
at $u^{*}$ and the reverse holds thanks to Lemma 3.2.

To conclude the proof, we have that it is easy to observe that every image
$u^{*}=cv$, with $c\in{\mathbb{R}}_{+}^{*}$, is a stationary point of the
energy $E(\cdot)$. Hence, such functions are solutions of the steady-state
equation (14). Since the osmosis evolution preserves the average grey value by
Proposition 3.1, one can further show that:

$\mu(f)=\mu(u^{*})=\frac{1}{|\Omega|}\int_{\Omega}cv(x)dx=c\mu(v),$

whence $c=\mu(f)/\mu(v)>0$. ∎

An analogous result can be proved as a corollary of the previous proposition
by allowing the further dependence on a smoothing parameter $p\in[1,2)$. In
Section 5, we will validate our model for different values of $p$, evaluating
the practical effects of its choice on exemplar tests.

######  Corollary 3.4.

Under the same assumptions of Proposition 3.3 and defining for $p\in[1,2)$ the
diffusivity term $g_{p}$ as

$g_{p}(x,t):=g_{p}(|\mathbf{s}(x,t)|):=\frac{1}{|\mathbf{s}(x,t)|^{p}},\qquad\forall(x,t)\in\Omega\times(0,T],$
(17)

then the steady-state of the energy functional

$E_{p}(u):=\int_{\Omega}\frac{v^{1-p}}{2-p}\left|\nabla\left(\frac{u}{v}\right)(x)\right|^{2-p}\;dx$
(18)

satisfies:

${\textnormal{div}}\left(g_{p}(|\nabla u^{*}-{\textbf{d}}u^{*}|)\left(\nabla
u^{*}-{\textbf{d}}u^{*}\right)\right)=0.$ (19)

A possible direction of future work consists in investigating the $p$-osmosis
model (19) so as to extend the osmosis filtering to a spatially adaptive
filtering, as proposed for image denoising, e.g., in [22].

##  4 Model discretisation using finite difference schemes

In this section we introduce fast and stable computational methods for the
nonlinear osmosis PDE model (6) by rewriting it using the linearity of the
divergence operator as follows

$\displaystyle\partial_{t}u={\textnormal{div}}\left(g(u)\nabla
u\right)-{\textnormal{div}}\left(g(u){\textbf{d}}u\right),$ (20)

along with the boundary conditions and initial condition $f\in
L^{\infty}(\Omega;{\mathbb{R}}_{+}^{*})$ as in (7).

We use Finite Difference Method (FDM) to approximate the solution of (20) on a
regular grid $\Omega_{h}\subset\Omega$ approximating the space domain $\Omega$
of the PDE, where $h>0$ denotes the spatial grid size along both horizontal
and vertical directions. We first apply a semi-discretisation in space on
$\Omega_{h}$, and then a full discretisation in time on $(0,T]$ using both
explicit and semi-implicit strategies, with the final objective of providing a
numerical scheme which shows stability to any fixed time-step discretisation
parameter.

In the following, we thus assume $u\in{\mathbb{R}}^{m_{x}\times m_{y}}$ to be
a discretised image defined on $\Omega_{h}$ of size $N:=m_{x}m_{y}$ and we
denote by $u_{i,j}$ the value of $u$ at pixel $(i,j)$ for $i=1,\ldots,m_{x}$,
$j=1,\ldots,m_{y}$. To avoid singularities appearing due to the special choice
of the diffusivity function $g$ in (10), we will use an $\epsilon$-regularised
version defined in terms of $0<\epsilon\ll 1$ as
$g(s)=1/|\mathbf{s}|_{\epsilon}$ with $\mathbf{s}=\nabla u-{\textbf{d}}u$,
with

$|\mathbf{s}|_{\epsilon}:=\sqrt{(s^{x})^{2}+(s^{y})^{2}+\epsilon}$

where $s^{x}$ and $s^{y}$ denote the horizontal and vertical components of
$\mathbf{s}$, respectively, and the smoothing parameter $\epsilon>0$ ensures
that $g(s)$ is well-defined at each point. Considering the canonical drift
term ${\textbf{d}}=\nabla\log v=\nabla v/v$, the value of $g$ at $(i,j)$ is
approximated using central finite difference schemes for the discretization of
the partial first-order derivatives by

$\displaystyle g_{i,j}$ $\displaystyle=g(s_{ij})$ (21) $\displaystyle s_{ij}$
$\displaystyle=\left(\frac{u_{i+1,j}-u_{i-1,j}}{2h}-\frac{v_{i+1,j}-v_{i-1,j}}{2h}\cdot\frac{u_{i,j}}{v_{i,j}}\;,\frac{u_{i,j+1}-u_{i,j-1}}{2h}-\frac{v_{i,j+1}-v_{i,j-1}}{2h}\cdot\frac{u_{i,j}}{v_{i,j}}\right).$

Figure 5: Finite difference stencil: the pixel neighbourhood of $(i,j)$ are
denoted by south (S), north (N), east (E), west (W).

###  4.1 Semi-discretisation in space

The spatial differential operators in (20), will be discretised by locating
the values of $\nabla u$, ${\textbf{d}}=(d^{x},d^{y})$, $u$ and the
$\epsilon$-regularised diffusivity function $g$ on the edges between two
adjacent pixels, following the stencil illustrated in Fig. 5. Hence, the semi-
discretisation in space of (20) at $(i,j)$ applying second-order central
finite differences reads:

$\displaystyle(\partial_{t}$ $\displaystyle
u)_{ij}=\frac{1}{h}\Big{[}\left[\nabla_{N}u\cdot g_{N_{i,j}}+\nabla_{S}u\cdot
g_{S_{i,j}}+\nabla_{E}u\cdot g_{E_{i,j}}+\nabla_{W}u\cdot g_{W_{i,j}}\right]+$
$\displaystyle-\left[d_{N}^{y}\cdot g_{N_{i,j}}\cdot
u_{N_{i,j}}+d^{y}_{S}\cdot g_{S_{i,j}}\cdot u_{S_{i,j}}+d_{E}^{x}\cdot
g_{E_{i,j}}\cdot u_{E_{i,j}}+d_{W}^{x}\cdot g_{W_{i,j}}\cdot
u_{E_{i,j}}\right]\Big{]}$ (22)

where each term is defined as in Table 1.

$\nabla_{N}u=\frac{u_{i+1,j}-u_{i,j}}{h}$, |  $\nabla_{S}u=\frac{u_{i-1,j}-u_{i,j}}{h}$, |  $\nabla_{E}u=\frac{u_{i,j+1}-u_{i,j}}{h}$, |  $\nabla_{W}u=\frac{u_{i,j-1}-u_{i,j}}{h}$,
---|---|---|---
$g_{N}=\frac{g_{i+1,j}+g_{i,j}}{2}$, |  $g_{S}=\frac{g_{i-1,j}+g_{i,j}}{2}$, |  $g_{E}=\frac{g_{i,j+1}+g_{i,j}}{2}$, |  $g_{W}=\frac{g_{i,j-1}+g_{i,j}}{2}$,
$u_{N}=\frac{u_{i+1,j}+u_{i,j}}{2}$, |  $u_{S}=\frac{u_{i-1,j}+u_{i,j}}{2}$, |  $u_{E}=\frac{u_{i,j+1}+u_{i,j}}{2}$, |  $u_{W}=\frac{u_{i,j-1}+u_{i,j}}{2}$,
$d_{N}^{y}=\frac{2(v_{i+1,j}-v_{i,j})}{h(v_{i+1,j}+v_{i,j})}$, |  $d_{S}^{y}=\frac{2(v_{i-1,j}-v_{i,j})}{h(v_{i-1,j}+v_{i,j})}$, |  $d_{E}^{x}=\frac{2(v_{i,j+1}-v_{i,j})}{h(v_{i,j+1}+v_{i,j})}$, |  $d_{W}^{x}=\frac{2(v_{i,j-1}-v_{i,j})}{h(v_{i,j-1}+v_{i,j})}$.
Table 1: FD discretization schemes.

Note that the discretisation (4.1) also holds on boundary pixels, having
preliminarily set ${\textbf{d}}=0$ on the boundary edges and mirroring the
image therein in order to satisfy the Neumann homogeneous boundary condition
(6).

We thus rewrite (4.1) by incorporating the discretisations reported in Table 1
and collecting the coefficients of the pixel $(i,j)$ and its neighbours:

$\displaystyle(\partial_{t}u)_{i,j}$
$\displaystyle=\,u_{i-1,j}\,\frac{g_{S}}{h}\left(\frac{1}{h}-\frac{d^{y}_{S}}{2}\right)+u_{i+1,j}\,\frac{g_{N}}{h}\left(\frac{1}{h}-\frac{d^{y}_{N}}{2}\right)+$
$\displaystyle+$ $\displaystyle
u_{i,j-1}\,\frac{g_{W}}{h}\left(\frac{1}{h}-\frac{d^{x}_{W}}{2}\right)+u_{i,j+1}\,\frac{g_{E}}{h}\left(\frac{1}{h}-\frac{d^{x}_{E}}{2}\right)+$
(23) $\displaystyle+$ $\displaystyle
u_{i,j}\Bigg{[}\frac{g_{S}}{h}\left(-\frac{1}{h}-\frac{d^{y}_{S}}{2}\right)+\frac{g_{N}}{h}\left(-\frac{1}{h}-\frac{d^{y}_{N}}{2}\right)+\frac{g_{W}}{h}\left(-\frac{1}{h}-\frac{d^{x}_{W}}{2}\right)+\frac{g_{E}}{h}\left(-\frac{1}{h}-\frac{d^{x}_{E}}{2}\right)\Bigg{]}.$

In order to represent (4.1) in matrix-vector form, we first rearrange into a
single vector $u(t)=(u_{1}(t),\ldots,u_{N}(t))^{T}\in{\mathbb{R}}^{N}$ the
values of the evolving image $u(T)$ at each pixel $i=1,\ldots,N$. Then, we
write semi-discretised-in-space model which leads to the following Cauchy
problem:

$\begin{cases}\frac{du}{dt}=A(u)\,u(t),\\\ u(0)=f,\end{cases}$ (24)

where $A(u)\in{\mathbb{R}}^{N\times N}$ is a non-symmetric sparse matrix, with
only five non-zero diagonals and with entries given by

$\displaystyle
a_{ij}:=\begin{cases}{\displaystyle\sum_{k\in\mathcal{N}(i)}\frac{g_{k}+g_{i}}{2h^{2}}\left(-1-\frac{v_{k}-v_{i}}{v_{k}+v_{i}}\right)}&\text{if
}i=j,\\\
{\displaystyle\frac{g_{j}+g_{i}}{2h^{2}}\left(1-\frac{v_{j}-v_{i}}{v_{j}+v_{i}}\right)}\quad&\text{if
$i\neq j$, $j\in\mathcal{N}(i)$},\\\ 0&\text{if $i\neq j$,
$j\not\in\mathcal{N}(i)$.}\end{cases}$ (25)

where $\mathcal{N}(i)$ represents the set of the neighbours in directions
$\\{N,E,O,W\\}$ of the pixel $i$ and $g_{i}$ is evaluated at pixel $i$
according to formula (21).

The next proposition shows some useful properties of the matrix $A(u)$.

######  Proposition 4.1.

For every vector $u\in{\mathbb{R}}^{N}$, let $A(u)\in{\mathbb{R}}^{N\times N}$
be the matrix with entries given in (25). Then, the following properties hold:

  1. 1.

$A(u)$ has non-negative off-diagonals,

  2. 2.

all column sums of $A(u)$ are 0,

  3. 3.

$A(u)$ is irreducible,

  4. 4.

$A(u)v=0$.

###### Proof.

1) According to (25), if pixel $i$ and pixel $j$ are not adjacent, then the
entry $a_{ij}$ is zero. On the other hand, as $g_{i},g_{j},v_{i},v_{j}>0$,
then $a_{ij}>0,\,\forall i,j$.

2) Considering an arbitrary column $l$ of $A$, we have:

$\displaystyle\sum_{i=1}^{N}a_{i,l}$
$\displaystyle=a_{l,l}+\sum_{j\in\mathcal{N}(l)}a_{j,l}=$
$\displaystyle=\sum_{j\in\mathcal{N}(l)}\frac{g_{j}+g_{l}}{2h^{2}}\left(-1-\frac{v_{j}-v_{l}}{v_{j}+v_{l}}\right)+\sum_{j\in\mathcal{N}(l)}\frac{g_{l}+g_{j}}{2h^{2}}\left(1-\frac{v_{l}-v_{j}}{v_{l}+v_{j}}\right)=0.$

3) Since $a_{i,j}>0$ for all adjacent $i,j$, so the directed graph $G_{A}$
associated to $A$ is strongly connected and $A$ is irreducible.

4) The value of the $ij$ coordinate of the vector $A(u)v$ is obtained from the
right hand side of the expression (4.1), substituting $u$ with $v$. Writing
explicitly every term as described in Table (1), it is easy to see that they
cancels each other. Therefore $A(u)v=0$. ∎

The properties showed in the lemma above are crucial for the following.
Typically, they are proved in the context of symmetric space-discretisation
matrices (see, e.g., [1]), but have been shown to hold also in the case of
non-symmetric diffusion-transport operators such as the osmosis one, see,
e.g., [8, 7].

###  4.2 Fully-discretisation in time

We present two possible fully discretised approximations of (24), obtained by
applying Euler integration methods in time. Let $\tau>0$ denote a uniform
time-step and for $k=0,1,\ldots$ let $u^{k}$ be the approximation of $u$ at
discrete time $k\tau\in(0,T]$.

  * 1.

Explicit discretisation: We first consider the simplest space discretisation
at the $k$th discrete time step that leads to the following full explicit
scheme

$\frac{u^{k+1}-u^{k}}{\tau}=A(u^{k})u^{k}\qquad\Longrightarrow\qquad
u^{k+1}=(I+\tau A(u^{k}))u^{k}.$ (E)

  * 2.

Semi-implicit discretisation: Alternatively, we can consider a semi-implicit
scheme where the nonlinear terms of the equation are treated from the previous
time step, while the linear ones are considered at the current time step; this
leads to the semi-implicit scheme

$\frac{u^{k+1}-u^{k}}{\tau}=A(u^{k})u^{k+1}\qquad\Longrightarrow\qquad
u^{k+1}=(I-\tau A(u^{k}))^{-1}u^{k}.$ (S.I.)

In both cases, the time-stepping process approximating (24) can be represented
in terms of a matrix $P\in{\mathbb{R}}^{N\times N}$ as

$\begin{cases}u^{k+1}=P(u^{k})\,u^{k},\\\ u^{0}=f.\end{cases}$ (26)

where

$P=I+\tau A,\quad\text{in}\quad\text{\eqref{eq:explicit}},\qquad P=(I-\tau
A)^{-1}\quad\text{in}\quad\text{\eqref{eq:semi-implicit}}.$ (27)

The following proposition provides sufficient conditions on matrix $P$ such
that the evolution process (26) is stable.

######  Lemma 4.2.

Let $f\in({\mathbb{R}}_{+}^{*})^{N}$ and let assume that
$P(u^{k})\in{\mathbb{R}}^{N\times N}$, defined in (27), satisfies the
following assumptions for all $k\in\mathbb{N}$:

  1. (A1)

all column sums of $P$ are 1,

  2. (A2)

$P$ is non-negative.

Then the iterative process (26) is stable with respect to the $1$-norm, i.e.
$\|u^{k}\|_{1}\leq\|f\|_{1}$ $\forall k\in\mathbb{N}$.

###### Proof.

For all $k$, we recursively have

$\displaystyle u^{k}$ $\displaystyle=P(u^{k-1})u^{k-1}=P(u^{k-1})\cdot
P(u^{k-2})\cdot\ldots\cdot P(u^{1})\cdot P(f)\cdot f$
$\displaystyle\|u^{k}\|_{1}$
$\displaystyle\leq\|P(u^{k-1})\|_{1}\cdot\|P(u^{k-2})\|_{1}\cdot\ldots\cdot\|P(u^{1})\|_{1}\cdot\|P(f)\|_{1}\cdot\|f\|_{1}.$

Since $P(u^{k})$ has column sums equal to one, then $\|P(u^{j})\|_{1}=1$ for
all $j=1,\ldots,k-1$ and $\|u^{k}\|_{1}\leq 1\cdot\|f\|_{1}$ for all
$k\in\mathbb{N}$. ∎

We are interested in verifying if the conservation properties of the non-
linear osmosis model proved in Proposition 3.1 for the continuum setting hold
in the discrete case. The following lemma adds two other conditions on the
matrix $P$ needed for the process to preserve the average mass and the non-
negativity. Furthermore, it provides a useful information on the localisation
of the eigenvalues of $P$.

######  Lemma 4.3.

Let $f\in({\mathbb{R}}_{+}^{*})^{N}$ and consider the iterative process (26),
where, for all $k$, the matrix $P(u^{k})\in{\mathbb{R}}^{N\times N}$ satisfies
both the assumptions (A1)-(A2) in Lemma 4.2 and the following

  * (A3)

$P$ is irreducible,

  * (A4)

$P$ has only positive diagonal entries.

Then we have

  1. (1)

the average grey value is preserved:

$\sum_{i=1}^{N}u_{i}^{k}=\sum_{i=1}^{N}f_{i},$

  2. (2)

the evolution preserves non-negativity: $u_{i}^{k}\geq 0$ for all $i=1,...,N$,
$k>0$,

  3. (3)

$P$ has a simple eigenvalue $\lambda=1$, with a corresponding positive
eigenvector, and the other eigenvalues have absolute value strictly less than
one.

###### Proof.

(1) Thanks to the hypotheses (A1) of Lemma 4.2, we can write

$\sum_{i=1}^{N}u_{i}^{k}=\sum_{i=1}^{N}\sum_{j=1}^{N}P_{ij}u_{j}^{k-1}=\sum_{j=1}^{N}\left(\sum_{i=1}^{N}P_{ij}\right)u_{j}^{k-1}=\sum_{j=1}^{N}u_{j}^{k-1}=\ldots=\sum_{j=1}^{N}f_{j}.$

(2) Applying the first step to the non-negative initial image $f$, we have

$u_{i}^{1}=p_{ii}f_{i}+\sum_{i\neq j}p_{ij}f_{j}\geq 0\qquad i=1,\ldots,N$

since $f_{i},f_{j},p_{ii},p_{ij}\geq 0$. This observation can be applied
recursively for all $k$.

(3) In Lemma 4.2 we proved that $\|P\|_{1}=1$ and, therefore, the eigenvalues
$\lambda$ of $P$ satisfies $|\lambda|\leq 1$. Moreover, $\lambda=1$ is an
eigenvalue of $P$. In fact, if we define the vector $e=(1,\ldots,1)$, then
$P^{T}e^{T}=e^{T}$, because column-wise $P$ sums up to one. By Gerschgorin’s
theorem, we know that the eigenvalues lie within the set $\Lambda$ of all
Gerschgorin’s disk:

$\Lambda:=\bigcup_{j=1}^{N}B_{j}\qquad\text{with}\qquad
B_{j}:=\left\\{z\in\mathbb{C}\,\Big{|}\,|z-p_{jj}|\leq\sum_{i\neq
j}|p_{ij}|\right\\}.$

Now, since $P$ is non-negative with unit column sums and positive diagonal
entries, the right hand side of the last inequality becomes

$\sum_{i\neq j}|p_{ij}|=\sum_{i\neq j}p_{ij}=1-p_{jj}\in[0,1).$

Therefore $\lambda=1$ is the only eigenvalue of $P$ s.t. $|\lambda|=1$, while
for the others $|\lambda|<1$. Moreover, $P$ has spectral radius $\rho(P)=1$.
By Perron-Frobenius theorem (see, e.g., Th.8.4.4 in [23]) we can thus conclude
as the spectral radius of a $P$ is a simple eigenvalue associated to a
positive eigenvector. ∎

Up to this point, we have defined sufficient conditions on the matrix $P$ to
make the evolution process (26) stable for any time step $k$, and guarantee
the desired conservation properties. Specifically, we stated that $P$ must
satisfy the assumptions (A1)-(A4) for all $k$.

With the following theorem, we show that both the explicit (E) and the semi-
implicit (S.I.) time-discretisation schemes verify such conditions.

######  Theorem 4.4 (Numerical Stability).

Let $f\in({\mathbb{R}}_{+}^{*})^{N}$ and for $u\in{\mathbb{R}}^{N}$ let $A$ be
the matrix with entries $a_{ij}$, with $i,j=1,\ldots,N$, defined in (25). Then

  1. 1.

the explicit scheme (E) is conditionally stable for

$\tau\leq 1/\max\\{|a_{11}|,...,|a_{NN}|\\}$ (28)

  2. 2.

the semi-implicit scheme (S.I.) is unconditionally stable for any $\tau>0$.

###### Proof.

Thanks to Lemma 4.2, we simply need to show that, in both cases, $P$ is non-
negative and has column sums equal to 1.

1) In the explicit scheme, $P=I+\tau A$, where $A$ has non-negative off-
diagonals and column sums zero. Thus $P$ has column sums 1 for any $\tau$ and
is non-negative only if condition (28) is satisfied.

2) In the semi-implicit scheme, $P=(I-\tau A)^{-1}$. $A$ has column sums zero
and for all $j=1,...,N$ we have

$\displaystyle-a_{jj}=\sum_{i\neq j}a_{ij}\quad\Longrightarrow\quad 1-\tau
a_{jj}>\tau\sum_{i\neq j}a_{ij}.$

Hence, the matrix $I-\tau A$ is strictly column diagonally dominant and non-
singular. Moreover, it has non-positive off-diagonal entries. Therefore,
$I-\tau A$ is a non-singular M-matrix (see [24], Th. 6.2.3 ($C_{10}$)) and its
inverse has only strictly positive entries (see [24], Th. 6.2.7). This
satisfies (A1) of Lemma 4.2. Lastly, all column sums of $P=(I-\tau A)^{-1}$
are 1. In fact, if we define $e=(1,\ldots,1)$, we have $e=e(I-\tau A)$ and,
equivalently, $e=e(I-\tau A)^{-1}$. ∎

The following result investigates under which sufficient conditions the
discretised solution obtained by solving (26) satisfies the conservation
properties introduced in Proposition 3.1 along the time iterations $k$ and how
a steady-state can be computed.

######  Theorem 4.5.

Let $f\in({\mathbb{R}}_{+}^{*})^{N}$. Then, the solution $u^{k}$ computed at
discrete time $k\tau$ with $g$ defined as in (21) preserves its mass and non-
negativity for any $k\geq 1$. Moreover, it converges to the steady-state
solution $u^{*}$ defined by

$u^{*}\,\,=\,\,\frac{\mu_{f}}{\mu_{v}}v,$ (29)

if

  1. 1.

$\tau\leq 1/\max\\{|a_{1,1}|,...,|a_{N,N}|\\}$, when the explicit scheme (E)
is applied

  2. 2.

$\tau>0$, when the semi-implicit scheme (S.I.) is applied.

###### Proof.

We start by proving that, in both cases, the assumptions of Lemma 4.3 are
verified. We have already verified in the proof of Theorem 4.5 when $P$
satisfies (A1) and (A2). Two other conditions, namely (A3) and (A4), still
need to be checked.

1) We know that $a_{ii}<0$ for all $i$ because $A$ has column sums zero and
non-negative off-diagonal entries. Therefore, the condition on $\tau$ ensures
that the diagonal entries of $P$ are strictly positive: $1+\tau
a_{ii}>0\;\forall i=1,...,N$. $P=I+\tau A$ is irreducible because it has the
same directed graph of $A$, which is irreducible.

2) In Theorem 4.5 we showed that $P=(I-\tau A)^{-1}$ has only strictly
positive entries. Therefore, its diagonal is positive and it is irreducible
because its directed graph is strongly connected. Hence, independently on
which time discretization we chose, explicit or semi-implicit, we can apply
Lemma 4.3, eventually with restrictions on $\tau$, which guarantees that the
evolution verifies the conservation properties. Moreover, the matrices
$P(u^{k})$ have eigenvalues with magnitude strictly less than one, with the
exception of the simple eigenvalue $\lambda=1$, that has a corresponding
positive eigenvector. Now, we observe that this eigenvector is the same for
all matrices $P(u^{k})$. In fact, since $A(u^{k})v=0$, we have:

  1. 1.

in the explicit case $P(u^{k})v=(I+\tau A(u^{k}))v=v$;

  2. 2.

in the semi-implicit case $P(u^{k})^{-1}v=(I-\tau A(u^{k}))v=v$ and, hence,
$P(u^{k})v=v.$

Therefore, the evolution $u^{k+1}=P(u^{k})u^{k}$ attenuates all components
outside the eigenspace of $\lambda=1$ to zero and the process converges to a
positive eigenvector $w$ with the same average grey value as the initial
vector $f$. In other terms, as $k\rightarrow\infty$, we have

$u^{k}\longrightarrow w=\frac{\mu_{f}}{\mu_{v}}v$

∎

Numerically, the proposed semi-implicit scheme (S.I.) can be solved at each
iteration by means of any efficient preconditioned linear iterative solver
suitable for sparse, diagonally dominant M-matrices, such as e.g., SOR
(successive over-relaxation) method.

##  5 Numerical results

As described in Section 2.1, upon specific choices of the drift term d in (6),
the nonlinear osmosis model can be applied to the shadow removal (see Section
5.1), light-spot removal (see Section 5.2) and the compact image
representation (see Section 5.3) task. We will compare qualitatively and
quantitatively the results by our approach with the ones obtained by other
osmosis-based methods as well as with state-of-the-art approaches.
Furthermore, some considerations on the computational efficiency of the
proposed models will be made. Reconstruction quality will be measured in terms
of the well-known Structural Similarity Index (SSIM $\in[0,1]$) error between
the original non-shadowed image and the reconstructed image, due to its
intrinsic dependence of luminance, contrast and structure features. Note that
the explicit iterative scheme (E) is practically unusable for most of the
examples considered, due to the constraint (28) on the time-step $\tau$ which
for the following examples forced constrains of the type $\tau<10^{-8}$ which
make the use of such scheme very limited in practice. On the contrary, the
unconditional stability of the semi-implicit scheme (S.I.) allows for the use
possibly large values of the time-step. In the examples, we used
$\tau=10^{3}$, which provides convergence and good level of accuracy in few
iterations. For the processing of colour images, the proposed model is run on
each R, G, B channel independently. All algorithms are tested on a AMD Ryzen 5
3450U processor (2.10GHz) with 16GB RAM using MATLAB R2021a.

###  5.1 Shadow removal

Recalling the shadow removal problem described in Section 5.1, we test the
non-linear model by setting the starting point of the osmosis evolution and
the reference image as $u^{0}=v=f$, the given image. We remind that for such
problem, the drift term d is defined as in (4). This choice corresponds to
consider two different evolutions in the regions composing the image domain
which explicitly read

$\displaystyle\begin{cases}\partial_{t}u=\displaystyle{{\textnormal{div}}\left(\frac{\nabla
u}{|\nabla u|}\right)}\quad&\text{on }\Omega_{b}\times(0,T]\\\ \\\
\partial_{t}u=\displaystyle{{\textnormal{div}}\left(\frac{\nabla
u-\frac{\nabla f}{f}u}{\left|\nabla u-\frac{\nabla
f}{f}u\right|}\right)}&\text{on }(\Omega_{in}\cup\Omega_{out})\times(0,T]\\\
\end{cases}$ (30)

with boundary conditions defined as in (7). Note that choosing d in a
canonical form according (4) allows to preserve the image content outside
$\Omega_{b}$, while diffusing non-linearly, in the form of a Total Variation
flow, on $\Omega_{b}$ itself. The intensity balancing effect in the regions
$\Omega_{in}$ and $\Omega_{out}$ depends solely by the evolution equation on
$\Omega_{b}$, the region where the drift term is set to zero. Moreover, the
diffusivity function $g$ allows to tune the diffusion intensity on the pixels
in $\Omega_{b}$ depending on the value of $\nabla u$. In particular, $g$ slows
down the diffusion on the pixels with high gradient values, which correspond
to edges that need to be preserved.

We test the performance of the proposed non-linear shadow removal model both
in case of hard shadows and for more realistic soft-edged shadows.

In the case of hard shadows, synthetic examples $f$ can be easily generated by
a constant multiplicative rescaling of a ground-truth image $f^{*}$ only in
correspondence of a bounded region $S\subset\Omega$. Denoting by $c\in(0,1)$
the (unknown) constant loss of luminosity on $S$, we can construct first the
shadow image by setting it as:

$s=\begin{cases}c\quad&\text{on}\;S\\\ 1&\text{on}\;\Omega\setminus
S,\end{cases}$ (31)

and obtain a shadowed version $f$ of $f^{*}$ by Hadamard (point-wise)
multiplication

$f=\begin{cases}f^{*}\odot s&\text{(hard shadow)}\\\
f^{*}\odot\left(G_{\sigma}*s\right)&\text{(soft shadow)}\end{cases}$ (32)

where ‘$*$’ denotes the standard convolution product and $G_{\sigma}$ is a
Gaussian convolution kernel with standard deviation $\sigma>0$. Such
convolution allows to generate an area of penumbra around shadow boundaries,
which is often the case in real-world images where the shadow edges can not be
precisely distinguished.

Fig. 6 shows the results of the proposed nonlinear osmosis model solved by
means of the semi-implicit scheme on two hard-shadowed images. For such
images, it is easy to distinguish whether a pixel belongs to the shadowed or
to the un-shadowed region. Hence, we can use a thin (2 pixel wide) mask
covering one pixel from each side of the boundary between the two regions.
This choice is sufficient to accurately remove the shadow while preserving all
features on the shadow boundary.

|  |
---|---|---
|  |
(a) Input $f=v$ | (b) mask 2px | (c) output
Figure 6: Shadow removal on images with hard shadows,.

Fig. 7 shows the results obtained for the soft-edged shadow removal problem in
correspondence of two different masks (wide and thin) and three input images
where shadows are characterised by penumbra boundaries of increasing dimension
(from left to right) corresponding to different increasing values of
$\sigma\in\\{0,1,2\\}$ in (32).

| $\sigma=0$ | $\sigma=1$ | $\sigma=2$
---|---|---|---
|  |  |
6px |  |  |
|  |  |
2px |  |  |
|  |  |
Figure 7: Shadow removal on input images $f$ corrupted by a constant shadow
convolved with convolution kernel $G_{\sigma}$ with increasing standard
deviation $\sigma$. Two different mask thicknesses, 6px and 2px, are
considered.

As seen before, the thin mask is a good choice for the case of hard shadows
(i.e. $\sigma=0$), but it becomes less effective as $\sigma$ increases, i.e.
as the penumbra region becomes wider, as it cannot cover appropriately the
boundary between the two regions $\Omega_{in}$ and $\Omega_{out}$. In these
cases, it is necessary to enlarge the shadow mask. Using a wider one, we
observe that the shadows are properly removed and features are well-preserved,
in particular near the pixels with limited intensity variation.

Note, however, that using a large mask may result in some drawbacks, as shown
in Fig. 8, where a wide (6 pixel) mask is used. As this choice enlarges the
region where the non-linear TV-type diffusion is induced, typical TV drawbacks
are observed. For instance, we observe in Fig. 8(b) that some artefacts appear
in the attempt of connecting circular structures (see the boundary between the
yellow and the green region). Furthermore, as it can be observed in Fig. 8(d),
when the shadow mask is too large some undesired connectivity effect is
observed. As it is well-known, the use of TV as a regulariser in inpainting
problems forces the minimisation of the length of level lines causing
artefacs, a property that can be seen through the use of the standard coarea
formula [4].

(a) Input $f=v$

(b) Output

(c) Input $f=v$

(d) Output

Figure 8: Drawbacks of shadow removal on images with soft shadows, using a 6
pixels wide mask.

A possible remedy for these two drawbacks consists in employing a different
diffusivity function $g$. As we remarked in Corollary 3.4, a slight variation
of the non-linear model depending on a diffusivity function $g_{p}$ depending
on a constant parameter $p\in[1,2)$ as in (17) can be considered. The use of
higher values of $p$ helps to better identify the regions of the shadow
boundary that are crossed by the underlying edges of the image. In these
regions, since $g_{p}(u)=1/|\nabla u|^{p}$, we have the effect of slowing down
the diffusion, thus favouring edge preservation.

|  |  |
---|---|---|---
|  |  |
|  |  |
Figure 9: Results obtained using the diffusivity function $g_{p}$ with $p=1$
(second row) and $p=1.9$ (third row).

In Figure 9 we show some results with $p=1$ (second row) and $p=1.9$ (third
row) on the images of Fig. 8 and on two other examples. The use of a possible
space-variant strategy adapting the value of $p$ to the local image content
would mitigate the previously observed drawbacks and, as such, is an
interesting direction of future research.

In the following, we compare qualitatively and quantitatively our results with
the ones obtained by means of the standard linear osmosis model (1) and its
anisotropic variant (3) proposed in [7] to limit the smoothing effect due to
Laplace diffusion on the shadow boundary.

In Fig. 10, some examples of hard shadow removal are shown. We notice that the
linear model is effective in removing the shadow, but produces visible
blurring on the shadow boundary. The anisotropic model is able to connect well
the underlying features of the image, but does not well balance the image
intensity between the inside and the outside of the shadow. Our model works
well for both tasks, achieving very good SSIM values in comparison with the
original ground truth images.

Input $f=v$ | Linear [6] | Anisotropic [7] | Non-linear
---|---|---|---
|  |  |
(a) | 0.9943 | 0.9966 | 0.9979
|  |  |
(b) | 0.9624 | 0.9544 | 0.9951
Figure 10: Comparison between linear, anisotropic and nonlinear osmosis models
on hard-shadowed images, with associated SSIM values.

In Fig. 11 we compare the results obtained by the three different models on
soft-shadowed images and real images. We observe that around the mask region
our model performs as good as the anisotropic one. However, while the
anisotropic model tends to make the whole image slightly blurred (due to the
smoothing kernels required to defined the diffusion tensor W), the non-linear
model preserves the sharpness of the image everywhere.

Input $f=v$ | Linear [6] | Anisotropic [7] | Non-linear
---|---|---|---
|  |  |
(c) | 0.9712 | 0.9937 | 0.9931
|  |  |
(d) | 0.9412 | 0.9464 | 0.9538
|  |  |
(e) |  |  |
|  |  |
(f) |  |  |
|  |  |
(g) |  |  |
Figure 11: Results of osmosis models on artificially soft-shadowed images and
natural images, with associated SSIM when computable.

We now compare the computational efficiency required to compute the numerical
solution of the non-linear model by using the semi-implicit scheme (S.I.) with
the linear [6] (solved with a fully implicit scheme with $\tau=10^{3}$), and
the anisotropic model [7] (solved by exponential integration).

| Linear [6] | Anisotropic [7] | Non-linear
---|---|---|---
(a) | 3.98 | (7) | 27.39 | (4) | 1.25 | (3)
(b) | 6.44 | (19) | 29.71 | (5) | 1.45 | (4)
(c) | 11.54 | (20) | 20.77 | (2) | 1.80 | (3)
(d) | 8.66 | (26) | 18.70 | (3) | 3.53 | (8)
(e) | 31.63 | (9) | 1307.06 | (10) | 10.37 | (3)
(f) | 105.55 | (30) | 2592.00 | (18) | 27.08 | (7)
(g) | 19.14 | (9) | 723.48 | (10) | 8.62 | (4)
Table 2: Execution time (in seconds) and number of iterations till convergence
(in brackets) for results in Fig. 10 (a,b) and Fig. 11 (c,d,e,f,g).

In this respect, we report in Table 2 the execution time (in seconds) and the
number of iterations (in brackets) for the three models when applied to the
images of Figure 10 and 11, using as same stopping criteria a mean squared
error w.r.t. to the ground truth smaller then a given tolerance for
artificially corrupted images (a,b,c,d), and a sufficiently small ($10^{-3}$)
relative change for real-world images (e,f,g). From a visual inspection of
Fig. 11, we can observe that the results obtained by the anisotropic model are
qualitatively comparable with the non-linear ones, although they require a
much more significant computational effort, see Table 2. This is due to the
more sophisticated nature of the anisotropic model (3) which requires the
computation of the matrix field W and its possible update at every iteration.
Interestingly, the proposed non-linear model appears faster than the linear
model in reaching convergence. This is due to the action of the diffusivity
term $g$, which accelerates the osmosis process on
$\Omega_{in}\cup\Omega_{out}$.

Finally, we provide an extensive list of comparisons with alternative
approaches proposed for shadow removal Illuminant invariant approaches are
used, e.g., in [14, 15], where an ”intrinsic” image is extracted from the
input image and then shadow removal is treated as reintegration problem, using
this image as a guide to derive gradient information or similarity measures
between shadowed and lighted regions. A variational approach to shadow and
spot-light removal has been proposed in [25], in the context of image
decomposition. Colour transfer techniques, such as [26, 27], use gradient
information to recover the shadow-free image by solving a Poisson equation.
Other methods consists in scaling the shadow region by a suitable factor,
modelling the shadow as uniform [28] or even non-uniform [29]. Another class
of solvers treat the task as a matting problem, considering the shadow pixels
as foreground and the lit pixels as background [30, 31].

|  |  |  |  |
---|---|---|---|---|---
|  | 0.9172 | 0.9603 | 0.9667 | 0.9676
|  |  |  |  |
|  | 0.7816 | 0.9482 | 0.9557 | 0.9726
|  |  |  |  |
|  | 0.8288 | 0.9092 | 0.9350 | 0.9488
|  |  |  |  |
|  | 0.8335 | 0.9312 | 0.9265 | 0.9650
|  |  |  |  |
|  | 0.8106 | 0.9500 | 0.9531 | 0.9733
|  |  |  |  |
|  | 0.6657 | 0.8174 | 0.8134 | 0.9008
Original | Input $f=v$ | [17] | [18] | [11] | ours
Figure 12: Comparison between shadow removal approaches, with associated SSIM
values.

In Fig. 12 we report a comparison between the results obtained by the proposed
non-linear osmosis model and three of the most recent and successful methods
for shadow removal555The test images have been downloaded from the ISTD
dataset, available at https://github.com/DeepInsight-PCALab/ST-CGAN.: a deep
learning approach, [18], a region-based shadow detection and removal [17] and,
finally, a variational approach based on a non-local version of the osmosis
model, [11] followed by a further contrast correction step as described by the
authors. We observe that our method produces visually pleasant results
compared to state-of-the-art methods and more sophisticated variants of the
standard osmosis model [11] further requiring a contrast correction step.
Quantitatively, it outperforms all other approaches in terms of SSIM values.

###  5.2 Spot-light removal

Recalling Section 2.1, an analogous use of the osmosis model (30) can be done
for the spot-light removal problem. Here, a synthetic spot-lighted image $f$
is obtained by Hadamard product $f\odot(G_{\sigma}*\ell)$ between a light-
spot-free image $f^{*}$ and a light image (with bounded lighted region
$L\subset\Omega$ ) $\ell$ defined as

$\ell=\begin{cases}c\quad&on\quad L\\\ 1&on\quad\Omega\setminus L,\end{cases}$
(33)

with $c>1$. Positive values of the standard deviation $\sigma$ produce soft-
lighted images, which require a larger mask. Fig. 13 shows the results for
four increasing levels of $\sigma$. Analogous comments to the ones made in the
case of shadow removal applications can be done in this case, too.

$\sigma=0$ |  |  |  |
---|---|---|---|---
$\sigma=2$ |  |  |  |
$\sigma=5$ |  |  |  |
$\sigma=8$ |  |  |  |
| (a) Input | (b) Results | (c) Input | (d) Results
Figure 13: Spot-light removal results. Input images were obtained by
convolving the light-free image with a Gaussian kernel with variance $\sigma$,
whose values are reported on the left column.

###  5.3 Compact data representation

As recalled in Section 5.3, image osmosis can also be applied to the problem
of compact data representation, by means of a proper definition of the drift
term. Given an image $v$ with average grey value $\mu_{v}$, we thus consider
$M=\Omega_{b}$ to be a pre-computed edge-mask via standard segmentation
algorithms (Canny, Sobel, …) and define the drift term as in (5). Explicitly,
the two equations guiding the osmosis evolution in this case then read:

$\displaystyle\begin{cases}\partial_{t}u={\textnormal{div}}\left(\frac{\nabla
u}{|\nabla u|}\right)\quad&\text{on }(\Omega\setminus\Omega_{b})\times[0,T]\\\
\\\ \partial_{t}u={\textnormal{div}}\left(\frac{\nabla u-\frac{\nabla
v}{v}u}{\left|\nabla u-\frac{\nabla v}{v}u\right|}\right)&\text{on
}\Omega_{b}\times[0,T]\end{cases}$ (34)

with initial and boundary conditions defined as in (7). By setting
$f\equiv\mu_{v}$ on $\Omega$ and making use of the actual values of the
reference image $v$ only on edges to define the drift term d on $\Omega_{b}$,
by evolving (34) the information is diffused on $\Omega\setminus\Omega_{b}$,
which can be thought of as the piecewise-constant region of the image. The
process thus converges to an approximation $w$ of the reference image $v$.

|  |  |
---|---|---|---
|  | 0.8581 | 0.8632
|  |  |
|  | 0.9432 | 0.9626
|  |  |
|  | 0.9058 | 0.9453
|  |  |
|  | 0.8827 | 0.9160
|  |  |
|  | 0.9671 | 0.9814
Input $v$ | Edge mask | Linear [8] | Non-linear
Figure 14: Compact data representation results with relative SSIM values. Edge
mask has been computed by standard Canny segmentation.

Fig. 14 compares the results obtained via linear and non-linear osmosis. In
general, the reconstructed images seem to have a lower contrast than the
original ones, but the non linear model favours a better approximation, as
confirmed by the SSIM values, in particular on strongly piecewise-constant
images, such as cartoons.

##  6 Conclusions

We presented a non-linear extension of the linear osmosis model originally
proposed in [6]. The non-linearity is encoded in the model through a
diffusivity term $g$ whose definition depends on the magnitude of the
diffusion-transport term, which allows the model to balance the intensity
values inside and outside the region of interest and, at the same time, to
preserve the underlying features on the boundary (i.e., edges), thus
preventing smoothing artefacts.

Forward-Euler discretisation is used to approximate the time evolution. As far
as spatial discretisation is concerned, suitable finite difference schemes are
defined. A semi-implicit scheme is used to deal with nonlinear term in a
computationally efficient way.

The proposed model enjoys conservation properties such average grey value and
non-negativity preservation, both in its continuous formulation and in its
discrete form. Unconditional stability of the semi-implicit iterative scheme
can be proved using standard scale-space properties. A fast convergence to the
steady-state solution can thus be obtained.

The efficiency and the numerical accuracy of the presented schemes are
validated through three imaging applications: shadow and spot-light removal,
and compact data representation. The results are accurate, can be computed
efficiently and outperform the ones obtained by alternative osmosis models and
state-of-the-art approaches both visually and in terms of SSIM values. The
simplicity of the model guarantees a low computational cost and, furthermore,
a fast computation. Moreover, and most importantly, it only requires a
preliminary segmentation of the region of interest (shadow, light-spot, edge
mask) without any other tuning of model hyperparameter, which makes it,
essentially, fully automatic.

An interesting future research direction is the study of a combined non-linear
and anisotropic osmosis model as well as the possible extension to higher-
order differential order so as to reduce the observed reconstruction
drawbacks.

## Acknowledgements

LC acknowledges the support received by the I3S laboratory of Sophia-
Antipolis, France, the CNRS PRIME project Imag’In and and by the EU H2020 RISE
projects NoMADS, GA 777826. Research by SM was supported in part by the
National Group for Scientific Computation (GNCS-INDAM), Research Projects
2021. SP acknowledges the support from the Leverhulme Trust Research Project
Grant (RPG-2018-121) “Unveiling the Invisible - Mathematics for Conservation
in Arts and Humanities”. GR acknowledges the support received by University of
Bologna and the I3S laboratory of Sophia-Antipolis, France.

## References

  * [1] J. Weickert, Anisotropic Diffusion in Image Processing, B.G. Teubner, Stuttgart, 1998.
  * [2] G. Aubert, P. Kornprobst, Mathematical problems in image processing: Partial Differential Equations and the Calculus of Variations, Vol. 147, Springer-Verlag, New York, 2006.  doi:10.1007/978-0-387-44588-5.
  * [3] T. Chan, J. Shen, Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods, Society for Industrial and Applied Mathematics, Philadelphia, 2005.  doi:10.1137/1.9780898717877.
  * [4] C.-B. Schönlieb, Partial Differential Equation Methods for Image Inpainting, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 2015.  doi:10.1017/CBO9780511734304.
  * [5] K. Hagenburg, M. Breuß, J. Weickert, O. Vogel, Novel schemes for hyperbolic pdes using osmosis filters from visual computing, 2011, pp. 532–543.  doi:10.1007/978-3-642-24785-9_45.
  * [6] J. Weickert, K. Hagenburg, M. Breuß, O. Vogel, Linear osmosis models for visual computing, in: A. Heyden, F. Kahl, C. Olsson, M. Oskarsson, X.-C. Tai (Eds.), Energy Minimization Methods in Computer Vision and Pattern Recognition, Springer Berlin Heidelberg, Berlin, Heidelberg, 2013, pp. 26–39.
  * [7] S. Parisotto, L. Calatroni, M. Caliari, C.-B. Schönlieb, J. Weickert, Anisotropic osmosis filtering for shadow removal in images, Inverse Problems 35 (5) (2019) 054001.  doi:10.1088/1361-6420/ab08d2.
  * [8] O. Vogel, K. Hagenburg, J. Weickert, S. Setzer, A fully discrete theory for linear osmosis filtering, in: A. Kuijper, K. Bredies, T. Pock, H. Bischof (Eds.), Scale Space and Variational Methods in Computer Vision, Springer Berlin Heidelberg, Berlin, Heidelberg, 2013, pp. 368–379.  doi:10.1007/978-3-642-38267-3_31.
  * [9] L. Calatroni, C. Estatico, N. Garibaldi, S. Parisotto, Alternating direction implicit (adi) schemes for a pde-based image osmosis model, Journal of Physics: Conference Series 904 (2017) 012014.  doi:10.1088/1742-6596/904/1/012014.
  * [10] S. Parisotto, L. Calatroni, A. Bugeau, N. Papadakis, C.-B. Schönlieb, Variational osmosis for non-linear image fusion, IEEE Transactions on Image Processing 29 (2020) 5507–5516.  doi:10.1109/TIP.2020.2983537.
  * [11] S. Benalia, M. Hachama, A nonlocal method for image shadow removal, Computers & Mathematics with Applications 107 (2022) 95–103.  doi:10.1016/j.camwa.2021.12.023.
  * [12] P. Pérez, M. Gangnet, A. Blake, Poisson image editing, ACM Trans. Graph. 22 (3) (2003) 313–318.  doi:10.1145/882262.882269.
  * [13] R. Fattal, D. Lischinski, M. Werman, Gradient domain high dynamic range compression, ACM Trans. Graph. 21 (3) (2002) 249–256.  doi:10.1145/566654.566573.
  * [14] G. Finlayson, S. Hordley, C. Lu, M. Drew, On the removal of shadows from images, IEEE Transactions on Pattern Analysis and Machine Intelligence 28 (1) (2006) 59–68.  doi:10.1109/TPAMI.2006.18.
  * [15] G. Finlayson, M. Drew, C. Lu, Entropy minimization for shadow removal, International Journal of Computer Vision 85 (1) (2009) 35–57.  doi:10.1007/s11263-009-0243-z.
  * [16] X. Chunxia, S. Ruiyun, X. Donglin, M. Kwanâ-Liu, Fast shadow removal using adaptive multi-scale illumination transfer, Computer Graphics Forum 32 (8) (2013) 207–218.  doi:10.1111/cgf.12198.
  * [17] R. Guo, Q. Dai, D. Hoiem, Paired regions for shadow detection and removal, IEEE Transactions on Pattern Analysis and Machine Intelligence 35 (12) (2013) 2956–2967.  doi:10.1109/TPAMI.2012.214.
  * [18] H. M. Le, D. Samaras, Shadow removal via shadow image decomposition, CoRR abs/1908.08628 (2019).  arXiv:1908.08628.
  * [19] M. Baba, N. Asada, Shadow removal from a real picture, in: ACM SIGGRAPH 2003 Sketches & Applications, SIGGRAPH ’03, ACM, New York, NY, USA, 2003, pp. 1–1.  doi:10.1145/965400.965488.
  * [20] S. Carlsson, Sketch based coding of grey level images, Signal Process. 15 (1) (1988) 57–83.  doi:10.1016/0165-1684(88)90028-X.
  * [21] M. Mainberger, A. Bruhn, J. Weickert, S. Forchhammer, Edge-based compression of cartoon-like images with homogeneous diffusion, Pattern Recognition 44 (2011) 1859–1873.  doi:10.1016/j.patcog.2010.08.004.
  * [22] G. Baravdish, Y. Cheng, O. Svensson, F. Åström, Extension of p-laplace operator for image denoising, in: L. Bociu, J.-A. Désidéri, A. Habbal (Eds.), System Modeling and Optimization, Springer International Publishing, Cham, 2016, pp. 107–116.  doi:10.1007/978-3-319-55795-3_9.
  * [23] R. A. Horn, C. R. Johnson, Matrix Analysis, Cambridge University Press, 1985.  doi:10.1017/CBO9780511810817.
  * [24] A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Society for Industrial and Applied Mathematics, 1994.  doi:10.1137/1.9781611971262.
  * [25] M. Huska, S. H. Kang, A. Lanza, S. Morigi, A variational approach to additive image decomposition into structure, harmonic, and oscillatory components, SIAM Journal on Imaging Sciences 14 (4) (2021) 1749–1789.  doi:10.1137/20M1355987.
  * [26] F. Liu, M. Gleicher, Texture-consistent shadow removal, in: D. Forsyth, P. Torr, A. Zisserman (Eds.), Computer Vision – ECCV 2008, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, pp. 437–450.  doi:10.1007/978-3-540-88693-8_32.
  * [27] Y. Shor, D. Lischinski, The shadow meets the mask: Pyramid‐based shadow removal, Comput. Graph. Forum 27 (2008) 577–586.  doi:10.1111/j.1467-8659.2008.01155.x.
  * [28] C. Fredembach, G. Finlayson, Simple shadow removal, in: 18th International Conference on Pattern Recognition (ICPR’06), Vol. 1, 2006, pp. 832–835.  doi:10.1109/ICPR.2006.1054.
  * [29] E. Arbel, H. Hel-Or, Texture-preserving shadow removal in color images containing curved surfaces, in: 2007 IEEE Conference on Computer Vision and Pattern Recognition, 2007, pp. 1–8.  doi:10.1109/CVPR.2007.383081.
  * [30] T.-P. Wu, C.-K. Tang, A bayesian approach for shadow extraction from a single image, in: Tenth IEEE International Conference on Computer Vision (ICCV’05) Volume 1, Vol. 1, 2005, pp. 480–487 Vol. 1.  doi:10.1109/ICCV.2005.4.
  * [31] T.-P. Wu, C.-K. Tang, M. S. Brown, H.-Y. Shum, Natural shadow matting, ACM Trans. Graph. 26 (2) (2007) 8–es.  doi:10.1145/1243980.1243982.