Stable Backward Diffusion Models that Minimise Convex Energies

Stable Bac kw ard Diffusion Mo dels that Minimise Con v ex Energies Leif Bergerhoff 1 , Marcelo Cárdenas 1 , Joac him W eic kert 1 , and Martin W elk 2 1 Mathematical Image Analysis Group, F aculty of Mathematics and Computer Science, Saarland Univ ersit y , Campus E1.7, 66041 Saarbrück en, Germany {bergerhoff,cardenas,weickert}@mia.uni-saarland.de 2 Institute of Biomedical Image Analysis, Priv ate Universit y for Health Sciences, Medical Informatics and T echnology , Eduard-W allnöfer-Zentrum 1, 6060 Hall/T yrol, Austria martin.welk@umit.at Abstract. The in v erse problem of bac kward diffusion is kno wn to be i ll-posed and highly unstable. Bac kward diffusion processes app ear naturally in image enhancemen t and deblur- ring applications. It is therefore greatly desirable to establish a backw ard diffusion mo del whic h implemen ts a smart stabilisation approac h that can b e used in combination with an easy to handle n umerical sc heme. So far, existing stabilisation strategies in literature require sophisticated n umerics to solve the underlying initial v alue problem. W e deriv e a class of space-discrete one-dimensional backw ard diffusion as gradient descent of energies where we gain stability by imp osing range constraints. Interestingly , these energies are ev en conv ex. F urthermore, we establish a comprehensiv e theory for the time-contin uous ev olution and w e sho w that stabilit y carries o v er to a simple explicit time discretisation of our model. Finally , w e confirm the stability and usefulness of our tec hnique in experiments in which w e enhance the con trast of digital greyscale and colour images. Keyw ords: In v erse Problem · Backw ard Diffusion · Mo delling · Con vex Energy · Gradient Descen t · Con trast Enhancemen t · Image Pro cessing Mathematics Sub ject Classification (2010): 76R50 · 60J60 · 58J65 · 37N30 · 37N40 · 65D18 · 68U10 1 In tro duction F orw ard diffusion pro cesses are well-suited to describ e the smo othing of a giv en signal or image. This pro cess of blurring implies a loss of high frequencies or details in the original data. As a result, the inv erse problem, backw ard diffusion, suffers from deficien t information which are needed to uniquely reconstruct the original data. The in tro duction of noise due to measured data increases this difficult y even further. Consequen tly , a solution to the inv erse problem – if it exists at all – is highly sensitiv e and hea vily dep ends on the input data: Ev en the smallest perturbation in the initial data can ha ve a large impact on the ev olution and ma y cause large deviations. Therefore, it b ecomes clear that bac kward diffusion pro cesses necessitate further stabilisation. Pr evious W ork on Backwar d Diffusion. Already more than 60 years ago, John [26] discussed the quality of a numerical solution to the in verse diffusion problem given that a solution exists, and that it is b ounded and non-negativ e. Since then, a large 2 n umber of differen t regularisation metho ds ha ve ev olved which achiev e stabilit y by b ounding the noise of the measured and the unperturb ed data [49], b y op erator splitting [27], by F ourier regularisation [13], or by a mo dified Tikhonov regulariser [55]. Hào and Duc [23] suggest a mollification metho d where stabilit y for the inv erse diffusion problem follows from a con volution with the Dirichlet kernel. In [24] the same authors pro vide a regularisation metho d for backw ard parab olic equations with time-dep enden t co efficien ts. T ernat et al. [51] suggest lo w-pass filters and fourth- order regularisation terms for stabilisation. Bac kward parab olic differential equations also enjo y high p opularit y in the im- age analysis comm unity where they ha ve e.g. b een used for image restoration and image deblurring resp ectiv ely . The first contribution to bac kw ard diffusion dates bac k to 1955 when Ko v ásznay and Joseph [28] prop osed to use the scaled negative Laplacian for contour enhancement. Gab or [14] observed that the isotropy of the Laplacian op erator leads to amplification of accidental noise at contour lines at the same time it enhances the con tour lines. As a remedy , he prop osed to restrict the con trast enhancemen t to the orthogonal contour direction and – in a second step – suggested additional smo othing in tangen t direction. Linden baum et al. [29] mak e use of a veraged deriv atives in order to improv e the directional sensitiv e filter b y Gab or. How ever, the authors p oin t out that smo othing in only one direction fav ours the emergence of artefacts in nearly isotropic image regions. They recommend to use the Perona-Malik filter [40] instead. F orces of Perona-Malik type are also used b y P ollak et al. [43] who sp ecify a family of ev olution equations to sharp en edges and suppress noise in the context of image segmentation. In [50], ter Haar Romeny et al. stress the influence of higher order time deriv atives on the Gaussian deblurred image. Referring to the heat equation, the authors express the time deriv atives in the spatial domain and approximate them using Gaussian deriv atives. Steiner et al. [47] highlight ho w backw ard diffusion can b e used for feature enhancement in planar curv es. In the field of image pro cessing, a frequently used stabilisation technique con- strains the extrema in order to enforce a maximum-minim um principle. This is e.g. implemen ted in the inv erse diffusion filter of Osher and Rudin [38]. It imp oses zero diffusivities at extrema and applies bac kward diffusion ev erywhere else. The so-called forw ard-and-backw ard (F AB) diffusion of Gilb oa et al. [19] follo ws a slightly differen t approac h. Closely related to the Perona-Malik filter [40] it uses negative diffusivities for a sp ecific range of gradien t magnitudes. On the other hand, it imp oses forw ard diffusion for v alues of low and zero gradien t magnitude. By doing so, the filter pre- v ents the output v alues from explo ding at extrema. Ho wev er, it is worth mentioning that – so far – all adequate implemen tations of inv erse diffusion pro cesses with for- w ard or zero diffusivities at extrema require sophisticated n umerical sc hemes. They use e.g. minmo d discretisations of the Laplacian [38], nonstandard finite difference appro ximations of the squared gradient magnitude [53], and splittings into t wo-pixel in teractions [54]. Another, less p opular stabilisation approac h implies the application of a fidelity term and has b een used to p enalise deviations from the input image [6,46] or from the av erage grey v alue of the desired range [45]. Consequently , b oth the w eights of the fidelit y and the diffusion term control the range of the filtered image. 3 F urther metho ds ac hieve stabilisation using a regularisation strategy built on FFT-based op erators [7,8,9] and by the restriction to p olynomials of fixed finite degree [25]. Mair et al. [31] discuss the w ell-p osedness of deblurring Gaussian blur in the discrete image domain based on analytic n umber theory . In summary , the presen ted methods offer an insigh t in to the c hallenge of handling bac kward diffusion in practice and underline the demand for careful stabilisation strategies and sophisticated n umerical metho ds. In our pap er w e are going to present an alternative approach to deal with bac k- w ard diffusion problems. It prefers smarter mo delling o ver smarter n umerics. T o understand it b etter, it is useful to recapitulate some relations b et w een diffusion and energy minimisation. Diffusion and Ener gy Minimisation. F or the sake of conv enience we assume a one- dimensional ev olution that smo othes an initial signal f : [ a, b ] → R . In this context, the original signal f serv es as initial state of the diffusion equation ∂ t u = ∂ x  g ( u 2 x ) u x  (1) where u = u ( x, t ) represen ts the filtered outcome with u ( x, 0) = f ( x ) . A dditionally , let u x = ∂ x u and assume reflecting boundary conditions at x = a and x = b . Giv en a nonnegative diffusivit y function g , growing diffusion times t lead to simpler represen tations of the input signal. F rom P erona and Malik’s work [40] we know that the smo othing effect at signal edges can b e reduced if g is a decreasing function of the contrast u 2 x . As long as the flux function Φ ( u x ) := g ( u 2 x ) u x is strictly increasing in u x the corresp onding forward diffusion pro cess ∂ t u = Φ 0 ( u x ) u xx in volv es no edge sharp ening. This diffusion can b e regarded as the gradient descent evolution whic h minimises the energy E [ u ] = Z b a Ψ ( u 2 x ) d x. (2) The p oten tial function ˜ Ψ ( u x ) = Ψ ( u 2 x ) is strictly con vex in u x , increasing in u 2 x , and fulfils Ψ 0 ( u 2 x ) = g ( u 2 x ) . F urthermore, the energy functional has a flat minimiser whic h is – due to the strict conv exity of the energy functional – unique. The gradient descen t / diffusion evolution is well-posed and con verges to w ards this minimiser for t → ∞ . Due to this classical emergence of w ell-p osed forw ard diffusion from strictly conv ex energies it seems natural to b eliev e that bac kward diffusion pro cesses are necessarily asso ciated with non-conv ex energies. How ev er, as we will see, this conjecture is wrong. Our Contribution. In our article, we show that a sp ecific class of bac kward diffu- sion pro cesses are gradient descent ev olutions of energies that hav e one unexp ected prop ert y: They are conv ex! Our second innov ation is the incorp oration of a sp ecific constrain t: W e imp ose reflecting boundary conditions in the diffusion c o-domain . This means that in case of greyscale images with an allo w ed grey v alue range of (0 , 255) the o ccurring v alues are mirrored at the b oundary p ositions 0 and 255 . While suc h range constrain ts hav e shown their usefulness in some other context (see e.g. [34]), to our knowledge they hav e nev er b een used for stabilising backw ard diffu- sions. F or our nov el backw ard diffusion mo dels, w e show also a surprising numerical 4 fact: A simple explicit sc heme turns out to b e stable and conv ergent. Last but not least, w e apply our mo dels to the con trast enhancement of greyscale and colour images. This article is a revised version of our conference con tribution [4] whic h w e extend in sev eral asp ects. First, w e enhance our mo del for conv ex bac kward diffusion to sup- p ort not only a global and w eigh ted setting but also a lo calised v ariant. W e analyse this extended mo del in terms of stability and con vergence tow ards a unique min- imiser. F urthermore, w e form ulate a simple explicit sc heme for our newly prop osed approac h which shares all imp ortan t prop erties with the time-contin uous evolution. In this context, we provide a detailed discussion on the selection of suitable time step sizes. Additionally , w e suggest tw o new applications: global con trast enhancemen t of digital colour images and lo cal con trast enhancement of digital grey and colour images. Structur e of the Pap er. In Section 2, we present our mo del for con vex bac kw ard diffusion with range constrain ts. W e describ e a general approac h whic h allo ws to form ulate w eighted lo cal and global ev olutions. Section 3 includes pro ofs for mo del prop erties such as range and rank-order preserv ation as well as conv ergence analysis and the deriv ation of explicit steady-state solutions. Section 4 provides a simple explicit sc heme whic h can b e used to solv e the o ccurring initial-v alue problem. In Section 5, we explain ho w to enhance the global and lo cal con trast of digital images using the prop osed model. F urthermore, w e discuss the relation to existing literature on con trast enhancement. In Section 6, w e draw conclusions from our findings and giv e an outlo ok on future researc h. 2 Mo del Let us no w explore the ro ots of our mo del and derive – in a second step – the particle ev olution which forms the heart of our metho d and whic h is given by the gradien t descen t of a conv ex energy . 2.1 Motiv ation from Sw arm Dynamics The idea b ehind our model go es bac k to the scenario of describing a one-dimensional ev olution of particles within a closed system. Recent literature on mathematical sw arm mo dels employs a pairwise p oten tial U : R n → R to characterise the b e- ha viour of individual particles (see e.g. [10,11,12,15,17] and the references therein). The p oten tial function allows to steer attractive and repulsive forces among sw arm mates. Physically simplified mo dels like [16] neglect inertia and describ e the indi- vidual particle v elo cit y ∂ t v i within a sw arm of size N directly as ∂ t v i = − N X j =1 j 6 = i ∇ U ( | v i − v j | ) , i = 1 , . . . , N , (3) where v i and v j denote particle p ositions in R n . These mo dels are also referred to as first order mo dels. Often they are inspired by biology and describ e long-range 5 attractiv e and short-range repulsiv e b eha viour b et ween sw arm mem b ers. The inter- pla y of attractive and repulsive forces leads to flo c king and allows to gain stabilit y for the whole sw arm. Inv erting this behaviour – resulting in short-range attractiv e and long-range repulsiv e forces – leads to a highly unstable scenario in whic h the sw arm splits up into small separating groups which might nev er reach a p oin t where they stand still. One w ould exp ect that a restriction to repulsive forces only will increase this instability even further. Ho wev er, we will present a mo del which cop es w ell with exactly this situation. In our set-up ev ery particle mov es within the op en 1 l 1 1 r 2 l 2 2 r 3 l 3 3 r 4 l 4 4 r -1 0 1 2 Fig. 1. F our particles with p ositions in (0 , 1) and their reflections at the left and right domain b oundary (lab elled l and r accordingly). P article 2, for example, gets rep elled by the particles 1 , 3 , 4 , 1 l , 2 r , 3 r , 4 r . in terv al (0 , 1) and has an in teraction radius of size 1 . Keeping this in mind, let us briefly examine the t w o main assumptions of the evolution. First, there exist reflec- tions for all particles at the left and righ t domain b oundary . Secondly , the particles rep el eac h other and – furthermore – get rep elled by the reflections. Ho wev er, due to the limited viewing range, only one of the tw o reflections of a certain particle is considered at any given time, namely the one which is closer to the reference particle (see Figure 1). A sp ecial case o ccurs if the reference particle is lo cated at p osition 0.5: the repulsiv e forces of b oth of its o wn reflections equal out. 2.2 Discrete V ariational Mo del W e prop ose a dynamical system which has its ro ots in a spatial discretisation of the energy functional (2). F urthermore, we make use of a decreasing energy function Ψ : R + 0 → R and a global range constrain t on u . The corresp onding flux function Φ is defined as Φ ( s ) := Ψ 0 ( s 2 ) s . Our goal is to describ e the ev olution of one-dimensional – not necessarily distinct – particle p ositions v i ∈ (0 , 1) , where i = 1 , . . . , N . Therefore, we extend the p osition v ector v = ( v 1 , . . . , v N ) T with the additional co ordinates v N +1 , . . . , v 2 N defined as v 2 N +1 − i := 2 − v i ∈ (1 , 2) . This extended p osition v ector v ∈ (0 , 2) 2 N allo ws to ev aluate the energy function E ( v , W ) = 1 4 · 2 N X i =1 2 N X j =1 w i,j · Ψ (( v j − v i ) 2 ) , (4) whic h mo dels the repulsion p oten tial betw een all positions v i and v j . The co ef- ficien t w i,j denotes entry j in ro w i of a constant non-negativ e weigh t matrix W = ( w i,j ) ∈ ( R + 0 ) 2 N × 2 N . It mo dels the imp ortance of the interaction b et ween 6 Ψ a,n ( s 2 ) Ψ 0 a,n ( s 2 ) Φ a,n ( s ) a ·  ( s − 1) 2 n − 1  a · n s · ( s − 1) 2 n − 1 a · n · ( s − 1) 2 n − 1 T able 1. One exemplary class of p enaliser functions Ψ ( s 2 ) for s ∈ [0 , 1] with n ∈ N , a > 0 and corresponding diffusivit y Ψ 0 ( s 2 ) and flux Φ ( s ) functions. p osition v i and v j . All diagonal e lemen ts of the w eight matrix are p ositiv e, i.e. w i,i > 0 , ∀ i ∈ { 1 , 2 , . . . , 2 N } . In addition, w e assume that the weigh ts for all ex- tended p ositions are the same as those for the original ones. Namely , w e hav e w i,j = w i, 2 N +1 − j = w 2 N +1 − i, j = w 2 N +1 − i, 2 N +1 − j (5) for 1 ≤ i, j ≤ N . F or the p enaliser function Ψ w e imp ose several restrictions whic h we discuss subse- quen tly . T able 1 shows one reasonable class of functions Ψ a,n as well as the corre- sp onding diffusivities Ψ 0 a,n and flux functions Φ a,n . In Figure 2 we provide an illus- tration of three functions using a = 1 and n = 1 , 2 , 3 . The p enaliser is constructed follo wing a three step pro cedure. First, the function Ψ ( s 2 ) is defined as a contin- uously differentiable, decreasing, and strictly con vex function for s ∈ [0 , 1] with Ψ (0) = 0 and Φ − (1) = 0 (left-sided deriv ativ e). Next, Ψ is extended to [ − 1 , 1] b y symmetry and to R by p eriodicity Ψ  (2 + s ) 2  = Ψ ( s 2 ) . This results in a p enaliser Ψ ( s 2 ) whic h is contin uously differen tiable ev erywhere except at even integers, where it is still contin uous. Note that Ψ ( s 2 ) is increasing on [ − 1 , 0] and [1 , 2] . The flux Φ is con tinuous and increasing in (0 , 2) with jump discontin uities at 0 and 2 (see Figure 2). F urthermore, we hav e that Φ ( s ) = − Φ ( − s ) and Φ (2 + s ) = Φ ( s ) . Exploiting the prop erties of Ψ allows us to rewrite (4) without the redundan t entries v N +1 , . . . , v 2 N as E ( v , W ) = 1 2 · N X i =1 N X j =1 w i,j ·  Ψ (( v j − v i ) 2 ) + Ψ (( v j + v i ) 2 )  . (6) A gradien t descent for (4) is given by ∂ t v i = − ∂ v i E ( v , W ) = X j ∈ J i 1 w i,j · Φ ( v j − v i ) , i = 1 , . . . , 2 N , (7) where v i no w are functions of the time t and J i 1 := { j ∈ { 1 , 2 , . . . , 2 N } | v j 6 = v i } . (8) Note that for 1 ≤ i, j ≤ N , thus | v j − v i | < 1 , the flux Φ ( v j − v i ) is negativ e for v j > v i and p ositiv e otherwise, th us driving v i alw ays a w ay from v j . This implies that we hav e negative diffusivities Ψ 0 for all | v j − v i | < 1 . Due to the conv exit y of Ψ ( s 2 ) , the absolute v alues of the repulsive forces Φ are decreasing with the distance b et w een v i and v j . W e remark that the jumps of Φ at 0 and 2 are not problematic here, as all p ositions v i and v j in the argument of Φ are distinct by the definition of J i 1 . Let us discuss shortly how the interv al constraint for the v i , i = 1 , . . . , N , is enforced in (4) and (7). First, notice that v 2 N +1 − i for i = 1 , . . . , N is the reflection 7 − 1 − 0 . 5 0 0 . 5 1 1 . 5 2 2 . 5 3 − 1 − 0 . 8 − 0 . 6 − 0 . 4 − 0 . 2 0 s ˜ Ψ ( s ) ˜ Ψ 1 , 1 ( s ) ˜ Ψ 1 , 2 ( s ) ˜ Ψ 1 , 3 ( s ) − 1 − 0 . 5 0 0 . 5 1 1 . 5 2 2 . 5 3 − 100 − 50 0 50 100 s ˜ Ψ 0 ( s ) ˜ Ψ 0 1 , 1 ( s ) ˜ Ψ 0 1 , 2 ( s ) ˜ Ψ 0 1 , 3 ( s ) − 1 − 0 . 5 0 0 . 5 1 1 . 5 2 2 . 5 3 − 2 0 2 s Φ ( s ) Φ 1 , 1 ( s ) Φ 1 , 2 ( s ) Φ 1 , 3 ( s ) Fig. 2. T op: Exemplary p enaliser functions ˜ Ψ 1 , 1 , ˜ Ψ 1 , 2 , and ˜ Ψ 1 , 3 extended to the in terv al [ − 1 , 3] by imposing symmetry and p erio dicity with ˜ Ψ a,n ( s ) := Ψ a,n ( s 2 ) . Middle: Corresp onding diffusivities ˜ Ψ 0 1 , 1 , ˜ Ψ 0 1 , 2 , and ˜ Ψ 0 1 , 3 with ˜ Ψ 0 a,n ( s ) := Ψ 0 a,n ( s 2 ) . Bottom: Correspond ing flux functions Φ 1 , 1 , Φ 1 , 2 , and Φ 1 , 3 with Φ a,n ( s ) = Ψ 0 a,n ( s 2 ) s . 8 of v i on the right interv al b oundary 1 . F or v i and v 2 N +1 − j with 1 ≤ i, j ≤ N and v 2 N +1 − j − v i < 1 there is a repulsiv e force due to Φ ( v 2 N +1 − j − v i ) < 0 that drives v i and v 2 N +1 − j a wa y from the right in terv al b oundary . The closer v i and v 2 N +1 − j come to this b oundary , the stronger is the repulsion. F or v 2 N +1 − j − v i > 1 , we hav e Φ ( v 2 N +1 − j − v i ) > 0 . By Φ ( v 2 N +1 − j − v i ) = Φ  (2 − v j ) − v i  = Φ  ( − v j ) − v i  , this can equally b e interpreted as a repulsion b et w een v i and − v j where − v j is the reflection of v j at the left interv al b oundary 0 . In this case the interaction b et w een v i and v 2 N +1 − j driv es v i and − v j a wa y from the left in terv al b oundary . Recapitulating b oth p ossible cases, it b ecomes clear that every v i is either rep elled from the reflection of v j at the left or at the righ t in terv al b oundary , but never from b oth at the same time. As ∂ t v 2 N +1 − i = − ∂ t v i holds in (7), the symmetry of v is preserved. Dropping the redundan t entries v N +1 , . . . , v 2 N , Equation (7) can b e rewritten as ∂ t v i = X j ∈ J i 2 w i,j · Φ ( v j − v i ) − N X j =1 w i,j · Φ ( v j + v i ) , (9) for i = 1 , . . . , N , where the second sum expresses the repulsions b et ween original and reflected co ordinates in a symmetric wa y . The set J i 2 is defined as J i 2 := { j ∈ { 1 , 2 , . . . , N } | v j 6 = v i } . (10) Equation (9) denotes a form ulation for pure repulsion amongst N different p ositions v i with stabilisation b eing achiev ed through the consideration of their reflections at the domain b oundary . It is worth mentioning that within (6) and (9) we only make use of the first N × N entries of W . In the following, we denote this submatrix by ˜ W and refer to its elements as ˜ w i,j . Given an initial vector f ∈ (0 , 1) N and initialising v i (0) = f i , v 2 N +1 − i (0) = 2 − f i for i = 1 , . . . , N , the gradien t descen t (7) and (9) ev olves v to wards a minimiser of E . 3 Theory Belo w w e pro vide a detailed analysis of the evolution and discuss its main prop erties. F or this purp ose w e consider the Hessian of (6) whose entries for 1 ≤ i ≤ N read ∂ v i v i E ( v , ˜ W ) = X j ∈ J i 2 ˜ w i,j · Φ 0 ( v j − v i ) + N X j =1 ˜ w i,j · Φ 0 ( v j + v i ) , (11) ∂ v i v j E ( v , ˜ W ) =    ˜ w i,j ·  Φ 0 ( v j + v i ) − Φ 0 ( v j − v i )  , ∀ j ∈ J i 2 , ˜ w i,j · Φ 0 ( v j + v i ) , ∀ j ∈ J i 3 , (12) where J i 3 := { j ∈ { 1 , 2 , . . . , N } | v i = v j } . (13) 9 3.1 General Results In a first step, let us inv estigate the well-posedness of the underlying initial v alue problem in the sense of Hadamard [22]. Theorem 1 (W ell-Posedness). L et Ψ = Ψ a,n as define d in T able 1. Then the initial value pr oblem (9) is wel l-p ose d sinc e (a) it has a solution, (b) the solution is unique, and (c) it dep ends c ontinuously on the initial c onditions. Pr o of. The initial v alue problem (9) can b e written as ˙ v ( t ) = f ( v ( t )) := − ∇ v E ( v ( t ) , W ) (14) v (0) = v 0 (15) with v ( t ) ∈ R 2 N and t ∈ R + 0 where we mak e use of the fact that W is a constant w eight matrix. In case f ( v ( t )) is con tinuously differentiable and Lipschitz contin uous all three conditions (a)–(c) hold. Existence and uniqueness directly follow from [39, c hapter 3.1, Theorem 3]. Contin uous dep endence on the initial conditions is guaran teed due to [39, chapter 2.3, Theorem 1] which is based on Gron wall’s Lemma [21]. Thus, let us no w prov e differen tiabilit y and Lipschitz contin uit y of f ( v ( t )) . Differ entiability: Differen tiability follows from the fact that all functions Φ a,n are con tinuously differen tiable. As a consequence the partial deriv atives of (8) w.r.t. v i exist for i = 1 , . . . , 2 N . Lipschitz Continuity: The Gershgorin circle theorem [18] allows to estimate a v alid Lipsc hitz constant L as an upp er b ound of the sp ectral radius of the Jacobian of (9). F or 1 ≤ i ≤ N the entries read ∂ v i ( ∂ t v i ) = − X j ∈ J i 2 ˜ w i,j ·  Φ 0 ( v j − v i ) + Φ 0 ( v j + v i )  − 2 · X j ∈ J i 3 ˜ w i,j · Φ 0 ( v j + v i ) , (16) ∂ v j ( ∂ t v i ) =    ˜ w i,j ·  Φ 0 ( v j − v i ) − Φ 0 ( v j + v i )  , ∀ j ∈ J i 2 , − ˜ w i,j · Φ 0 ( v j + v i ) , ∀ j ∈ J i 3 . (17) The radii of the Gershgorin discs fulfil r i = N X j =1 j 6 = i   ∂ v j ( ∂ t v i )   = X j ∈ J i 2 ˜ w i,j · | Φ 0 ( v j − v i ) − Φ 0 ( v j + v i ) | + X j ∈ J i 3 j 6 = i ˜ w i,j · | Φ 0 ( v j + v i ) | < X j ∈ J i 2 ˜ w i,j · | Φ 0 ( v j − v i ) − Φ 0 ( v j + v i ) | + X j ∈ J i 3 ˜ w i,j · | Φ 0 ( v j + v i ) | =: ˜ r i , i = 1 , . . . , N . (18) 10 Then we ha ve | λ i − ∂ v i ( ∂ t v i ) | < ˜ r i for 1 ≤ i ≤ N where λ i denotes the i -th eigenv alue of the Jacobian of (9). This leads to the b ounds λ i < X j ∈ J i 2 ˜ w i,j ·  | Φ 0 ( v j − v i ) − Φ 0 ( v j + v i ) | − ( Φ 0 ( v j − v i ) + Φ 0 ( v j + v i ))  + X j ∈ J i 3 ˜ w i,j ·  | Φ 0 ( v j + v i ) | − 2 · Φ 0 ( v j + v i )  ≤ X j ∈ J i 2 ˜ w i,j ·  | Φ 0 ( v j − v i ) | + | Φ 0 ( v j + v i ) | + | Φ 0 ( v j − v i ) | + | Φ 0 ( v j + v i ) |  + X j ∈ J i 3 ˜ w i,j ·  | Φ 0 ( v j + v i ) | + 2 · | Φ 0 ( v j + v i ) |  ≤ 4 · L Φ · X j ∈ J i 2 ˜ w i,j + 3 · L Φ · X j ∈ J i 3 ˜ w i,j < 4 · L Φ · N X j =1 ˜ w i,j , i = 1 , . . . , N , (19) where L Φ represen ts the Lipschitz constant of the flux function Φ . Using the same reasoning one can sho w that λ i > − 4 · L Φ · N X j =1 ˜ w i,j , i = 1 , . . . , N . (20) Consequen tly , an upp er b ound for the sp ectral radius – and thus for the Lipsc hitz constan t L of the gradient of (9) – reads L ≤ max 1 ≤ i ≤ N | λ i | < 4 · L Φ · max 1 ≤ i ≤ N N X j =1 ˜ w i,j =: L max . (21) F or our sp ecific class of flux functions Φ a,n a v alid Lipschitz constant L Φ is giv en b y L Φ = a · n · (2 n − 1) · 2 2 n − 2 (22) suc h that we hav e L < 4 · a · n · (2 n − 1) · 2 2 n − 2 · max 1 ≤ i ≤ N N X j =1 ˜ w i,j . (23) This concludes the pro of. u t Next, let us show that no p osition v i can ever reach or cross the interv al b oundaries 0 and 1 . Theorem 2 (A voidance of Range Interv al Boundaries). F or any weighting matrix ˜ W ∈ ( R + 0 ) N × N al l N p ositions v i which evolve ac c or ding to (9) and have an arbitr ary initial value in (0 , 1) do not r e ach the domain b oundaries 0 and 1 for any time t ≥ 0 . 11 Pr o of. Equation (9) can b e written as ∂ t v i = X j ∈ J i 2 ˜ w i,j ·  Φ ( v j − v i ) − Φ ( v j + v i )  − X j ∈ J i 3 ˜ w i,j · Φ (2 v i ) , (24) where 1 ≤ i ≤ N . Notice that for j ∈ J i 2 w e hav e lim v i → 0 + Φ ( v j − v i ) − Φ ( v j + v i ) = 0 , (25) lim v i → 1 − Φ ( v j − v i ) − Φ ( v j + v i ) = 0 , (26) where the latter follows from the p erio dicit y of Φ . Consequently , any p osition v i whic h gets arbitrarily close to one of the domain b oundaries 0 or 1 exp eriences no impact b y p ositions v j with j ∈ J i 2 , and the first sum in (24) gets zero. The definition of Ψ ( s 2 ) implies that Ψ 0 ( s 2 ) < 0 , ∀ s ∈ (0 , 1) , (27) Ψ 0 ( s 2 ) > 0 , ∀ s ∈ (1 , 2) , (28) from whic h it follows for 1 ≤ i ≤ N that − Φ (2 v i ) > 0 , ∀ v i ∈  0 , 1 2  , (29) − Φ (2 v i ) < 0 , ∀ v i ∈  1 2 , 1  . (30) No w remember that ˜ W ∈ ( R + 0 ) N × N and ˜ w i,i > 0 . In combination with (29) and (30) w e get lim v i → 0 + ∂ t v i > 0 and lim v i → 1 − ∂ t v i < 0 , (31) whic h concludes the pro of. u t Let us for a moment assume that the p enaliser function is given by Ψ = Ψ a,n from T able 1. Below, w e pro ve that this implies con vergence to the global minimum of the energy E ( v , ˜ W ) . Theorem 3 (Conv ergence for Ψ = Ψ a,n =1 ). F or t → ∞ , given a p enaliser Ψ a, 1 with arbitr ary a > 0 , any initial c onfigur ation v ∈ (0 , 1) N c onver ges to a unique ste ady state v ∗ which is the glob al minimiser of the ener gy given in (6) . Pr o of. As a sum of conv ex functions, (6) is con vex. Therefore, the function V ( v , ˜ W ) := E ( v , ˜ W ) − E ( v ∗ , ˜ W ) (where v ∗ is the equilibrium p oin t) is a Ly apunov function with V ( v ∗ , ˜ W ) = 0 and V ( v , ˜ W ) > 0 for all v 6 = v ∗ . F urthermore, we hav e ∂ t V ( v , ˜ W ) = − N X i =1  ∂ v i E ( v , ˜ W )  2 ≤ 0 . (32) A ccording to Gershgorin’s theorem [18], one can show that the Hessian matrix of (6) is p ositiv e definite for Ψ = Ψ a, 1 from whic h it follows that E ( v , ˜ W ) has a strict (global) minimum. This implies that the inequalit y in (32) b ecomes strict except in case of v = v ∗ , and guarantees asymptotic Ly apunov stability [30] of v ∗ . Thus, we ha ve conv ergence to v ∗ for t → ∞ . u t 12 R emark 1. Theorem 3 can b e extended to the case of n = 2 and – in a weak er form ulation – to arbitrary n ∈ N . The pro ofs for b oth cases are based on a straight- forw ard application of the Gershgorin circle theorem. F or details w e refer to the supplemen tary material. (a) Giv en that Ψ = Ψ a,n =2 , let us assume that one of the follo wing tw o conditions – v i 6 = 1 2 , or – there exists j ∈ J i 2 for whic h v j 6 = 1 − v i and ˜ w i,j > 0 , is fulfilled for every i ∈ [1 , N ] and t ≥ 0 . Then the Hessian matrix of (6) is p ositiv e definite and con vergence to the strict global minim um of E ( v , ˜ W ) follo ws. (b) F or all p enaliser functions Ψ = Ψ a,n , one can sho w that the Hessian matrix of (6) is p ositiv e semi-definite. This means that our metho d con verges to a global minim um of E ( v , ˜ W ) . How ev er, this minimum do es not ha ve to b e unique. In general, the steady-state solution of (9) dep ends on the definition of the p enaliser function Ψ . Based on (24), and assuming that Ψ = Ψ a,n , a minimiser of E ( v , ˜ W ) necessarily fulfils the equation 0 = X j ∈ J i 2 ˜ w i,j ·  ( v ∗ j − v ∗ i − 1) 2 n − 1 − ( v ∗ j + v ∗ i − 1) 2 n − 1  − X j ∈ J i 3 ˜ w i,j · (2 v ∗ i − 1) 2 n − 1 , (33) where i = 1 , . . . , N . 3.2 Global Mo del If all p ositions v i in teract with each other during the evolution, i.e. ˜ w i,j > 0 for 1 ≤ i, j ≤ N , we sp eak of our mo del as acting glob al ly . Below, w e prov e the existence of weigh t matrices ˜ W for whic h distinct p ositions v i and v j (with i 6 = j ) can never b ecome equal (assuming that the p ositions v i , i = 1 , . . . , N , are distinct for t = 0 ). This implies that the initial rank-order of v i is preserv ed throughout the ev olution. Theorem 4 (Distinctness of v i and v j ). A mong N initial ly distinct p ositions v i ∈ (0 , 1) evolving ac c or ding to (9) , no two ever b e c ome e qual if ˜ w j,k = ˜ w i,k > 0 for 1 ≤ i, j, k ≤ N , i 6 = j . Pr o of. Giv en N distinct p ositions v i ∈ (0 , 1) , equation (9) can b e written as ∂ t v i = N X k =1 k 6 = i ˜ w i,k · Φ ( v k − v i ) − N X k =1 ˜ w i,k · Φ ( v k + v i ) , (34) for i = 1 , . . . , N . W e use this equation to derive the difference ∂ t ( v j − v i ) = ( ˜ w j,i + ˜ w i,j ) · Φ ( v i − v j ) + N X k =1 k 6 = i,j  ˜ w j,k · Φ ( v k − v j ) − ˜ w i,k · Φ ( v k − v i )  − N X k =1  ˜ w j,k · Φ ( v k + v j ) − ˜ w i,k · Φ ( v k + v i )  , (35) 13 where 1 ≤ i, j ≤ N . Assume w.l.o.g. that v j > v i and consider (35) in the limit v j − v i → 0 . Then w e hav e lim v j − v i → 0 ( ˜ w j,i + ˜ w i,j ) · Φ ( v i − v j ) > 0 , (36) if ˜ w j,i + ˜ w i,j > 0 , which ev ery global mo del fulfils by the assumption that ˜ w i,j > 0 for 1 ≤ i, j ≤ N . This follows from the fact that Φ ( s ) > 0 for s ∈ ( − 1 , 0) . F urthermore, w e hav e lim v j − v i → 0 N X k =1 k 6 = i,j  ˜ w j,k · Φ ( v k − v j ) − ˜ w i,k · Φ ( v k − v i )  − N X k =1  ˜ w j,k · Φ ( v k + v j ) − ˜ w i,k · Φ ( v k + v i )  = N X k =1 k 6 = i,j ( ˜ w j,k − ˜ w i,k ) · Φ ( v k − v i ) − N X k =1 ( ˜ w j,k − ˜ w i,k ) · Φ ( v k + v i ) = 0 if ˜ w j,k = ˜ w i,k for 1 ≤ k ≤ N . (37) In conclusion, w e can guarantee for global mo dels with distinct particle p ositions that lim v j − v i → 0 ∂ t ( v j − v i ) > 0 , (38) if ˜ w j,k = ˜ w i,k where 1 ≤ i, j, k ≤ N and i 6 = j . A ccording to (38), v j will alw ays start moving aw a y from v i (and vice versa) when the difference b et ween b oth gets sufficien tly small. Since the initial p ositions are distinct, it follows that v i 6 = v j for i 6 = j for all times t . u t A sp ecial case o ccurs if all en tries of the weigh t matrix ˜ W are set to 1 – i.e. ˜ W = 11 T with 1 := (1 , . . . , 1) T ∈ R N . F or this scenario, w e obtain an analytic steady-state solution whic h is indep enden t of the p enaliser Ψ : Theorem 5 (Analytic Steady-State Solution for ˜ W = 11 T ). Under the as- sumption that ( v i ) is in incr e asing or der, ˜ W = 11 T , and that Ψ ( s 2 ) is twic e c on- tinuously differ entiable in (0 , 2) the unique minimiser of (4) is given by v ∗ = ( v ∗ 1 , . . . , v ∗ 2 N ) T , v ∗ i = ( i − 1 / 2 ) / N , i = 1 , . . . , 2 N . Pr o of. With ˜ W = 11 T , Equation (4) can b e rewritten without the redundan t en tries of v as E ( v ) = N − 1 X i =1 N X j = i +1 Ψ (( v j − v i ) 2 ) + 1 2 · N X i =1 Ψ (4 v 2 i ) + N − 1 X i =1 N X j = i +1 Ψ (( v j + v i ) 2 ) . (39) F rom this, one can v erify b y straigh tforward, albeit length y calculations that ∇ E ( v ∗ ) = 0 . Moreo ver, one finds that the Hessian of E at v ∗ is D 2 E ( v ∗ ) = N X k =1 A k Φ 0  k N  . (40) 14 1 2 3 4 5 6 7 0 1 initial state 1 2 3 4 5 6 7 0 1 steady state Fig. 3. Application of the global mo del to a system of 7 particles with weigh t matrix ˜ W = 11 T . 1 2 3 4 5 6 7 0 1 initial state 1 2 3 4 5 6 7 0 1 steady state Fig. 4. Application of the global mo del to a system of 7 particles with ˜ w i,k = 1 / k for 1 ≤ i, k ≤ N . Here, A k are sparse symmetric N × N -matrices given by A k = 2 I − T k − T − k + H k +1 + H 2 N − k +1 , (41) A N = I + H N +1 , (42) for k = 1 , . . . , N − 1 , where the unit matrix I , the single-diagonal T o eplitz matrices T k , and the single-antidiagonal Hank el matrices H k are defined as I =  δ i,j  N i,j =1 , (43) T k =  δ j − i,k  N i,j =1 , (44) H k =  δ i + j,k  N i,j =1 . (45) Here, δ i,j denotes the Kronec ker sym b ol, δ i,j = 1 if i = j , and δ i,j = 0 otherwise. All A k , k = 1 , . . . , N are w eakly diagonally dominan t with p ositiv e diagonal, thus p ositiv e semidefinite b y Gershgorin’s Theorem. Moreo ver, the tridiagonal matrix A 1 is of full rank, th us ev en p ositiv e definite. By strict conv exit y of Ψ ( s 2 ) , all Φ 0 ( k / N ) are p ositiv e, th us D 2 E ( v ∗ ) is p ositiv e definite. As a consequence, the steady state of the gradien t descen t (9) for any initial data f (with arbitrary rank-order) can – under the condition that ˜ W = 11 T – b e computed directly by sorting the f i : Let σ b e the p erm utation of { 1 , . . . , N } for whic h ( f σ − 1 ( i ) ) i =1 ,...,N is increasing (this is what a sorting algorithm computes), the steady state is giv en by v ∗ i = ( σ ( i ) − 1 / 2 ) / N for i = 1 , . . . , N (cf. Figure 3). u t A dditionally , we presen t an analytic expression for the steady state of the global mo del given a p enaliser function Ψ = Ψ a,n (cf. T able 1) with n = 1 . 15 Theorem 6 (Analytic Steady-State Solution for Ψ = Ψ a,n =1 ). Given N dis- tinct p ositions v i in incr e asing or der and a p enaliser function Ψ = Ψ a,n =1 , the unique minimiser of (4) is given by v ∗ i = i P j =1 ˜ w i,j − 1 2 ˜ w i,i N P j =1 ˜ w i,j , i = 1 , . . . , N . (46) Pr o of. The presented minimiser follows directly from (33). Figure 4 pro vides an illustration of the steady state. u t Finally , and in case all entries of the w eigh t matrix ˜ W are set to 1 , we show that the global mo del conv erges – indep enden tly of Ψ – to a unique steady state: Theorem 7 (Conv ergence for ˜ W = 11 T ). Given that ˜ W = 11 T , any initial c onfigur ation v ∈ (0 , 1) N with distinct entries c onver ges to a unique ste ady state v ∗ for t → ∞ . This is the glob al minimiser of the ener gy given in (6) . Pr o of. Using the same reasoning as in the pro of for Theorem 3 w e kno w that in- equalit y (32) holds. Due to the p ositiv e definiteness of (40) it follo ws that E ( v , ˜ W ) has a strict (global) minimum which implies that the inequality in (32) b ecomes strict except in case of v = v ∗ . This guarantees asymptotic Ly apunov stability of v ∗ and th us conv ergence to v ∗ for t → ∞ . 3.3 Relation to V ariational Signal and Image Filtering Let us now in terpret v 1 , . . . , v N as samples of a smo oth 1D signal u : Ω → [0 , 1] o ver an in terv al Ω of the real axis, taken at sampling p ositions x i = x 0 + i h with grid mesh size h > 0 . W e consider the mo del (4) with w i,j := γ ( | x j − x i | ) , where γ : R + 0 → [0 , 1] is a non-increasing weigh ting function with compact supp ort [0 , % ) . Theorem 8 (Space-Contin uous Energy). Equation (4) c an b e c onsider e d as a discr etisation of E [ u ] = 1 2 Z Ω  W ( u 2 x ) + B ( u )  d x (47) with p enaliser W ( u 2 x ) ≈ C Ψ ( u 2 x ) and b arrier function B ( u ) ≈ D Ψ (4 u 2 ) , wher e C and D ar e p ositive c onstants. R emark 2. (a) The p enaliser W is decreasing and con vex in u x . The barrier function B is conv ex and it enforces the in terv al constrain t on u b y fav ouring v alues u aw a y from the in terv al b oundaries. The discrete p enaliser Ψ generates b oth the p enaliser W for deriv atives and the barrier function B . (b) Note that b y construction of W the diffusivity g ( u 2 x ) := W 0 ( u 2 x ) ∼ Ψ 0 ( u 2 x ) has a singularit y at 0 with −∞ as limit. (c) The cut-off of γ at radius % implies the locality of the functional (47) that can thereby b e linked to a diffusion equation of t yp e (1). With o ut a cut-off, a nonlo cal diffusion equation would arise instead. Pr o of (of The or em 8). W e notice first that v j − v i and v i + v j for 1 ≤ i, j ≤ N are first-order appro ximations of ( j − i ) h u x ( x i ) and 2 u ( x i ) , resp ectiv ely . 16 Derivation of the Penaliser W . Assume first for simplicity that Ψ ( s 2 ) = − κs , κ > 0 is linear in s on [0 , 1] (thus not strictly conv ex). Then we hav e for a part of the inner sums of (4) corresp onding to a fixed i : 1 2  N X j =1 γ ( | x j − x i | ) · Ψ  ( v j − v i ) 2  + 2 N X j = N +1 γ ( | x j − x 2 N +1 − i | ) · Ψ  ( v j − v 2 N +1 − i ) 2   = N X j =1 γ ( | x j − x i | ) · Ψ  ( | v j − v i | ) 2  ≈ − κ h u x ( x i ) N X j =1 γ ( | j − i | h ) · | j − i | = h Ψ  u x ( x i ) 2  N − i X k =1 − i | k | γ ( | k | h ) ≈ hΨ ( u x ( x i ) 2 ) · 2 h 2 Z % 0 z γ ( z )d z =: hC Ψ ( u x ( x i ) 2 ) , (48) where in the last step the sum ov er k = 1 − i, . . . , N − i has b een replaced with a sum ov er k = −b %/h c , . . . , b %/h c , th us introducing a cutoff error for those lo cations x i that are within the distance % from the interv al ends. Summation o ver i = 1 , . . . , N appro ximates R Ω C Ψ ( u 2 x )d x from whic h we can read off W ( u 2 x ) ≈ C Ψ ( u 2 x ) . F or Ψ ( s 2 ) that are non-linear in s , Ψ ( u x ( x i ) 2 ) in (48) is changed in to a weigh ted sum of Ψ  ( k u x ( x i )) 2  for k = 1 , . . . , N − 1 , which still amounts to a decreasing function W ( u 2 x ) that is conv ex in u x . Qualitatively , W 0 then b eha ves the same wa y as b efore. Derivation of the Barrier F unction B . Collecting the summands of (4) that were not used in (48), w e hav e, again for fixed i , 1 2  2 N X j = N +1 γ  | x j − x i |  · Ψ  ( v j − v i ) 2  + N X j =1 γ  | x j − x 2 N +1 − i |  · Ψ  ( v j − v 2 N +1 − i ) 2   = N X j =1 γ  | x j − x i |  · Ψ  ( v i + v j ) 2  ≈  2 h Z % 0 γ ( z )d z + 1  · Ψ (4 u ( x i ) 2 ) =: hD · Ψ (4 u ( x i ) 2 ) , (49) and thus after summation o ver i analogous to the previous step R Ω B ( u ) d x with B ( u ) ≈ D Ψ (4 u 2 ) . u t 17 Similar deriv ations can b e made for patches of 2D images. A p oin t w orth noticing is that the barrier function B is b ounded. This differs from usual contin uous mo dels where such barrier functions tend to infinit y at the interv al b oundaries. Ho wev er, for eac h giv en sampling grid and patc h size the barrier function is just strong enough to prev ent W from pushing the v alues out of the in terv al. 4 Explicit Time Discretisation Up to this p oin t we ha ve established a theory for the time-contin uous ev olution of particle p ositions. In order to b e able to emplo y our mo del in sim ulations and appli- cations we need to discretise (9) in time. Subsequently , we pro vide a simple y et p o w- erful discretisation whic h preserves all imp ortan t prop erties of the time-contin uous mo del. An approximation of the time deriv ative in (9) b y forw ard differences yields the explicit sc heme v k +1 i = v k i + τ · X ` ∈ J i 2 ˜ w i,` · Φ ( v k ` − v k i ) − τ · N X ` =1 ˜ w i,` · Φ ( v k ` + v k i ) , (50) for i = 1 , . . . , N , where τ denotes the time step size and an upp er index k refers to the time k τ . In the follo wing, we deriv e necessary conditions for which the explicit sc heme preserv es the p osition range (0 , 1) and the p osition ordering. F urthermore, w e show conv ergence of (50) in dep endence of τ . Theorem 9 (A voidance of Range In terv al Boundaries of the Explicit Sc heme). L et L Φ b e the Lipschitz c onstant of Φ r estricte d to the interval (0 , 2) . Mor e over, let 0 < v k i < 1 , for every 1 ≤ i ≤ N , and assume that the time step size τ of the explicit scheme (50) satisfies 0 < τ < 1 2 · L Φ · max 1 ≤ i ≤ N N P ` =1 ˜ w i,` . (51) Then it fol lows that 0 < v k +1 i < 1 for every 1 ≤ i ≤ N . Pr o of. In accordance with (24) the explicit sc heme (50) can b e written as v k +1 i = v k i + τ · X ` ∈ J i 2 ˜ w i,` ·  Φ ( v k ` − v k i ) − Φ ( v k ` + v k i )  − τ · X ` ∈ J i 3 ˜ w i,` · Φ (2 v k i ) , (52) where i = 1 , . . . , N . No w assume that 0 < v k i , v k j < 1 and let us examine the con- tribution of the t wo summation terms in (52). W e need to distinguish the follo wing fiv e cases: 1. If v k i = v k j ≤ 1 2 then 2 v k i ∈ (0 , 1] . Thus, 0 ≤ − Φ (2 v k i ) . (53) 2. If 1 2 < v k i = v k j then 2 v k i ∈ (1 , 2) . Thus, using Φ (1) = 0 , | Φ (2 v k i ) | = | Φ (2 v k i ) − Φ (1) | ≤ | 2 v k i − 1 | · L Φ < 2 v k i · L Φ . (54) 18 3. If v k i < v k j then v k j − v k i , v k j + v k i ∈ (0 , 2) . Thus, | Φ ( v k j + v k i ) − Φ ( v k j − v k i ) | ≤ L Φ · 2 v k i . (55) 4. If v k j < v k i ≤ 1 2 then v k j − v k i ∈ ( − 1 , 0) and v k j + v k i ∈ (0 , 1) . Thus, 0 ≤ Φ ( v k j − v k i ) − Φ ( v k j + v k i ) , (56) 0 ≤ − Φ (2 v k i ) . (57) 5. Finally , if v k j < v k i and 1 2 < v k i , using the p eriodicity of Φ we get | Φ ( v k j − v k i ) − Φ ( v k j + v k i ) | = | Φ ( v k j + v k i ) − Φ (2 + v k j − v k i ) | ≤ 2 v k i · L Φ . (58) Com bining (50) with (51) and (53)–(58) we obtain that v k +1 i − v k i = − τ · X ` ∈ J i 2 ˜ w i,` ·  Φ ( v k ` + v k i ) − Φ ( v k ` − v k i )  − τ · X ` ∈ J i 3 ˜ w i,` · Φ (2 v k i ) ≥ − τ · L Φ · 2 v k i · N X ` =1 ˜ w i,` > − v k i , (59) from whic h it directly follows that v k +1 i > 0 , as claimed. The pro of for v k +1 i < 1 is straightforw ard. Assume w.l.o.g. that ˜ v k i := 1 − v k i . F or the reasons given ab o v e, we obtain ˜ v k +1 i > 0 . Consequently , 1 − v k +1 i > 0 and v k +1 i < 1 follo ws. u t Theorem 10 (Rank-Order Preserv ation of the Explicit Sc heme). L et L Φ b e the Lipschitz c onstant of Φ r estricte d to the interval (0 , 2) . F urthermor e, let v 0 i , for i = 1 , . . . , N , denote the initial ly distinct p ositions in (0 , 1) and – in ac c or danc e with The or em 4 – let the weight matrix ˜ W have c onstant c olumns, i.e. ˜ w j,` = ˜ w i,` for 1 ≤ i, j, ` ≤ N . Mor e over, let 0 < v k i < v k j < 1 and assume that the time step size τ use d in the explicit scheme (50) satisfies 0 < τ < 1 2 · L Φ · max 1 ≤ i ≤ N N P ` =1 ˜ w i,` . (60) Then we have v k +1 i < v k +1 j . Pr o of. F or distinct p ositions, (50) reads v k +1 i = v k i + τ · N X ` =1 ` 6 = i ˜ w i,` · Φ ( v k ` − v k i ) − τ · N X ` =1 ˜ w i,` · Φ ( v k ` + v k i ) (61) 19 for i = 1 , . . . , N . Considering this explicit discretisation for ∂ t v i and ∂ t v j w e obtain for i, j ∈ { 1 , 2 , . . . , N } : v k +1 j − v k +1 i = v k j − v k i + τ · ( ˜ w j,i + ˜ w i,j ) · Φ ( v k i − v k j ) + τ · N X ` =1 ` 6 = i,j  ˜ w j,` · Φ ( v k ` − v k j ) − ˜ w i,` · Φ ( v k ` − v k i )  − τ · N X ` =1  ˜ w j,` · Φ ( v k ` + v k j ) − ˜ w i,` · Φ ( v k ` + v k i )  . (62) No w remember that v k i < v k j b y assumption and that – as a consequence – τ · ( ˜ w j,i + ˜ w i,j ) · Φ ( v k i − v k j ) > 0 . (63) Using the fact that ˜ w j,k = ˜ w i,k for 1 ≤ i, j, k ≤ N and that Φ is Lipschitz in (0 , 2) , w e also know that T 1 := τ · N X ` =1 ` 6 = i,j    ˜ w j,` · Φ ( v k ` − v k j ) − ˜ w i,` · Φ ( v k ` − v k i )    = τ · N X ` =1 ` 6 = i,j ˜ w j,` ·    Φ ( v k ` − v k j ) − Φ ( v k ` − v k i )    ≤ τ · L Φ · | v k i − v k j | · N X ` =1 ` 6 = i,j ˜ w j,` , (64) T 2 := τ · N X ` =1    ˜ w j,` · Φ ( v k ` + v k j ) − ˜ w i,` · Φ ( v k ` + v k i )    = τ · N X ` =1 ˜ w j,` ·    Φ ( v k ` + v k j ) − Φ ( v k ` + v k i )    ≤ τ · L Φ · | v k j − v k i | · N X ` =1 ˜ w j,` . (65) Let the time step size τ fulfil (60). Then we can write T 1 + T 2 < 2 · L Φ · 2 · | v k j − v k i | · N X ` =1 ˜ w j,` < v k j − v k i . (66) In com bination with T 1 , T 2 ≥ 0 , it follo ws that T 2 − T 1 ≥ − T 2 − T 1 > − ( v k j − v k i ) , (67) and w e immediately kno w that v k j − v k i − T 1 + T 2 > 0 . T ogether with (62) and (63) w e get 0 < v k +1 j − v k +1 i , as claimed. u t 20 Theorem 11 (Conv ergence of the Explicit Sc heme). L et (6) b e a twic e c on- tinuously differ entiable c onvex function. Then the explicit scheme (50) c onver ges for time step sizes 0 < τ ≤ 1 2 · L Φ · max 1 ≤ i ≤ N N P j =1 ˜ w i,j < 2 L , (68) wher e L Φ denotes the Lipschitz c onstant of Φ r estricte d to the interval (0 , 2) and L r efers to the Lipschitz c onstant of the gr adient of (6) . Pr o of. Con vergence of the gradient metho d to the global minimum of E ( v , ˜ W ) is w ell-known for contin uously differen tiable conv ex functions with Lipschitz contin- uous gradient and time step sizes 0 < τ < 2 /L (see e.g. [33, Theorem 2.1.14]). A v alid Lipsc hitz constant is giv en by L max as defined in (21). Consequen tly , the time step sizes τ need to fulfil (68) in order to ensure conv ergence of (50). The smaller or equal relation results from (21). The latter defines L max > L such that τ = 2 /L max represen ts a v alid time step size. u t R emark 3 (Optimal Time Step Size). The optimal time step size, i.e. the v alue of τ whic h leads to most rapid descen t, is giv en b y τ = 1 /L (see e.g. [33, Section 2.1.5]). Th us, we suggest to use τ = 1 /L max . 5 Application to Image Enhancemen t No w that we ha ve presen ted a stable and con vergen t numerical sc heme, we apply (50) to enhance the contrast of digital greyscale and colour images. Throughout all exp erimen ts w e use Ψ = Ψ 1 , 1 (cf. T able 1 and Figure 2). 5.1 Greyscale Images The application of the prop osed mo del to greyscale images follows the ideas pre- sen ted in [5]. W e define a digital greyscale image as a mapping f : { 1 , . . . , n x } × { 1 , . . . , n y } → [0 , 1] . Note that all grey v alues are mapp ed to the in terv al (0 , 1) to ensure the v alidit y of our mo del b efore pro cessing. The grid p osition of the i -th im- age pixel is giv en by the v ector x i whereas v i denotes the corresp onding grey v alue. Subsequen tly , w e will see that a well-considered choice of the weigh ting matrix ˜ W allo ws either to enhance the glob al or the lo c al contrast of a giv en image. Global Con trast Enhancement F or global con trast enhancement we make use of the global mo del as discussed in Section 3.2. Only the N different o ccurring grey v alues v i – and not their p ositions in the image – are considered. W e let every en try ˜ w i,j of the weigh ting matrix denote the frequency of grey v alue v j in the image. Assuming an 8-bit greyscale image this leads to a w eighting matrix of size 256 × 256 whic h is indep enden t of the image size. As illustrated in Figure 5, global con trast enhancemen t can b e achiev ed in t w o wa ys: As a first option one can use the explicit sc heme (50) to describ e the evolution of all grey v alues v i up to some time t (see column tw o of Figure 5). The amount of con trast enhancemen t gro ws with increasing 21 Original image t = 2 · 10 − 6 Steady state (46) Fig. 5. Global contrast enhancement using Φ = Φ 1 , 1 and greyscale versions of images from the BSDS500 [2]. v alues of t . In our exp erimen ts an image size of 481 × 321 pixels and the application of the flux function Φ 1 , 1 with L Φ = 1 imply an upp er b ound of 1 / (2 · 481 · 321) for τ . Th us, w e can achiev e the time t = 2 · 10 − 6 in Figure 5 in a single iteration. If one is only in terested in an enhanced version of the original image with maximum global con trast there is an alternative, namely the derived steady state solution for linear flux functions (46). The results are shown in the last column of Figure 5. This figure also confirms that the solution of the explicit sc heme (50) con verges to the steady- state solution (46) for t → ∞ . F rom (46) it is clear that this steady state is equiv alent to histogram equalisation. In summary , this means that the application of our global mo del to greyscale images offers an evolution equation histogram equalisation which allo ws to con trol the amount of contrast enhancemen t in an intuitiv e wa y through the time parameter t . Lo cal Con trast Enhancemen t In order to achiev e lo cal contrast enhancemen t w e use our mo del to describ e the evolution of grey v alues v i at all n x · n y image grid p ositions. The change of ev ery grey v alue v i dep ends on all grey v alues within a disk-shap ed neighbourho o d of radius % around its grid p osition x i . W e assume that ˜ w i,j := γ ( | x j − x i | ) , ∀ i, j ∈ { 1 , 2 , . . . , N } , (69) 22 where w e weigh t the spatial distance | x j − x i | b y a function γ : R + 0 → [0 , 1] with compact supp ort [0 , % ) which fulfils γ ( x ) ∈ (0 , 1] , if x < %, γ ( x ) = 0 , if x ≥ %. (70) The choice of γ is application dep enden t. Ho wev er, it usually makes sense to define γ ( x ) as a non-increasing function in x . Possible choices are e.g. γ 1 ( x ) = ( 1 , if x < %, 0 , else , (71) γ 2 ( x ) =      1 − 6 x 2 % 2 + 6 x 3 % 3 , if 0 ≤ x < % 2 , 2 · (1 − x % ) 3 , if % 2 ≤ x < %, 0 , else , (72) whic h are b oth sketc hed in Figure 6. When applying this lo cal mo del to images x % γ 1 ( x ) 1 x % γ 2 ( x ) 1 Fig. 6. Bo x function γ 1 and scaled cubic B-spline γ 2 . w e make use of mirroring b oundary conditions in order to av oid artefacts at the image b oundaries. Figure 7 pro vides an example for lo cal contrast enhancement of digital greyscale images. Again, we describ e the grey v alue ev olution with the explicit scheme (50). F urthermore, we use γ 1 to mo del the influence of neigh b ouring grey v alues. As is eviden t from Figure 7, increasing the v alues for t go es along with enhanced lo cal contrast. 5.2 Colour Images Based on the assumption that our input data is given in sRGB colour space [48] (in the following denoted b y R GB) w e represen t a digital colour image by the mapping f : { 1 , . . . , n x } × { 1 , . . . , n y } → [0 , 1] 3 . Subsequently , our aim is the contrast en- hancemen t of digital colour images without distorting the colour information. This means that w e only w ant to adapt the luminanc e but not the chr omaticity of a giv en image. F or this purp ose, w e con vert the giv en image data to YCbCr colour space [44, Section 3.5] since this representation provides a separate luminance c han- nel. Next, we p erform contrast enhancement on the luminance channel only . Just as for greyscale images we map all Y-v alues to the in terv al (0 , 1) to fulfil our mo del 23 Original image t = 2 · 10 − 5 t = 4 · 10 − 5 Fig. 7. Lo cal con trast enhancemen t using Φ = Φ 1 , 1 , γ = γ 1 , % = 60 , and greyscale versions of images from the BSDS500 [2]. requiremen ts. After enhancing the con trast, we transform the colour information of the image bac k to RGB colour space. A t this p oin t it is important to men tion that the colour gam ut of the R GB colour space is a subset of the YCbCr colour gam ut and during the con version pro cess of colour co ordinates from YCbCr to RGB colour space the so-called colour gam ut problem may o ccur: Colours from the YCbCr colour gamut ma y lie outside the RGB colour gamut and th us cannot b e represented in R GB colour co ordinates. Naik and Murth y [32] state that a simple clipping of the v alues to the b ounds creates undesired shift of hue and ma y lead to colour artefacts. In order to av oid the colour gam ut problem we adapt the ideas presen ted by Nikolo v a and Steidl [34] which are based on the intensit y representation of the HSI colour space [20, Section 6.2.3]. Using the original and enhanced intensities, they define an affine colour mapping and transform the original R GB v alues. This preserves the hue and results in an enhanced RGB image. It is straightforw ard to show that their algorithms are v alid for an y intensit y ˆ f of type ˆ f = c r · r + c g · g + c b · b, (73) with c r + c g + c b = 1 and c r , c g , c b ∈ [0 , 1] , where r , g , and b denote RGB colour co ordinates. Th us, they are applicable to the luminance represen tation of the YCbCr colour space, to o, i.e. c r = 0 . 299 , c g = 0 . 587 , c b = 0 . 114 . Tian and Cohen make use of the same idea in [52]. As in [34], our result image is a conv ex com bination of the 24 outcomes of a multiplicativ e and an additive algorithm (see [34, Algorithm 4 and 5]) with co efficien ts λ and 1 − λ for λ ∈ [0 , 1] . During our exp erimen ts w e use a fixed v alue of λ = 0 . 5 (for details on ho w to choose λ w e refer to [34]). An ov erview of our strategy for contrast enhancement of digital colour v alue images is given in Figure 8. Input Image sR GB Con v erted Image YCbCr Enhanced Image YCbCr Multiplicativ e Enhanced Image sR GB A dditiv e Enhanced Image sR GB Output Image sR GB sR GB to YCbCr Bac kw ard Diffusion Affine Colour T ransform Con v ex Com bination Fig. 8. Pro cedure of contrast enhancemen t for digital colour images following [34]. Global Con trast Enhancemen t Again, we apply the global mo del from Section 3.2 in order to ac hieve global con trast enhancement. As mentioned b efore, we con- sider the N different o ccurring Y-v alues of the YCbCr represen tation of the input image and denote them b y v i (similar to Section 5.1 w e neglect their p ositions in the image). Ev ery en try of the w eighting matrix ˜ w i,j con tains the n umber of o ccurences of the v alue v j in the Y-channel of the image. It becomes clear that the application of our mo del – in this setting – basically comes down to histogram equalisation of the Y-c hannel. Figure 9 shows the resulting RGB images after global contrast enhance- men t. Similar to the greyscale scenario, we can either apply the explicit scheme (50) or – for Φ = Φ a, 1 – estimate the steady state solution following (46). F or the first case the amoun t of con trast enhancement grows with the p ositiv e time parameter t . The second column of Figure 9 sho ws the results for Φ = Φ 1 , 1 giv en time t . The corresp onding steady state solutions are illustrated in the last column of Figure 9. Lo cal Con trast Enhancement In a similar manner – and adapting the ideas from Subsection 5.1 – w e ac hieve lo cal con trast enhancemen t in colour images. F or this purp ose w e describ e the evolution of Y-v alues v i at all n x · n y image grid p osi- tions using a disk-shap ed neighbourho o d of radius % around the corresp onding grid p ositions x i . The entries of the weigh ting matrix ˜ W follow (69). In com bination with mirrored b oundary conditions the explicit scheme (50) allo ws to increase the 25 Original image t = 5 · 10 − 7 Steady state (46) Fig. 9. Global contrast enhancement using Φ = Φ 1 , 1 , λ = 0 . 5 , and images from [1]. Original image t = 1 · 10 − 5 t = 2 · 10 − 5 Fig. 10. Lo cal contrast enhancemen t using Φ = Φ 1 , 1 , γ = γ 1 , % = 60 , λ = 0 . 5 , and images from [1]. 26 lo cal contrast of an image with gro wing t . Figure 10 sho ws exemplary results for Φ = Φ 1 , 1 and γ = γ 1 (cf. (71)). Note, how well – in comparison to the global set-up in Figure 9 – the structure of the do or gets enhanced while the details of the do or knob are preserv ed. The differences are even larger in the second image: F or b oth the coup le in the foreground and the background scenery , contrast increases whic h implies visibilit y also for larger times t . 5.3 P arameters In total, our mo del has up to six parameters: Φ , τ , t , λ , % , and γ . During our exp erimen ts w e hav e fixed Φ ( s ) to the linear flux function Φ 1 , 1 ( s ) and λ to 0 . 5 . V alid b ounds for the time step size τ are giv en in Theorems 9 – 11. F rom the theory in Section 3 and the subsequent experiments on greyscale and colour images it b ecomes clear that the amount of con trast enhancemen t gro ws with the diffusion time. Thus, it remains to discuss the influence of the parameters % and γ . W e found out that the neigh b ourho od radius % affects the diffusion time and con trols the amount of p erceiv ed lo cal con trast enhancement, i.e. it steers the lo calisation of the con trast enhancemen t pro cess. Whereas small radii lead to high contrast in already small image areas, the size of image sections with high con trast increases with % . F or sufficien tly large v alues of % global histogram equalisation is approximated. Another in teresting p oin t is the c hoice of the weigh ting function γ . Ov erall, choosing γ = γ 1 leads to more homogeneous con trast enhancemen t resulting in smo other p erception. F or γ = γ 2 the fo cus alwa ys lies on the neighbourho o d cen tre whic h implies ev en more enhancement of lo cal structures than in the preceding case. In summary , γ 2 leads to more enhancemen t which, how ev er, also creates undesired effects in smo oth or noisy regions. Thus, w e prefer γ 1 o ver γ 2 . F urther exp erimen ts which visualise the effect of the parameters can b e found in the supplementary material. 5.4 Related W ork from an Application Perspective No w that w e ha ve demonstrated the applicabilit y of our mo del to digital images we w ant to discuss briefly its relation to other existing theories in the context of image pro cessing. As mentioned in Section 5.1, applying th e global mo del – with the entries of ˜ W represen ting the grey v alue frequencies – is identical to histogram equalisation (a common formulation can e.g. b e found in [20]). F urthermore, there exist other closely related histogram sp ecification tec hniques – suc h as [45,35,36] – which can hav e the same steady state. If we compare our ev olution with the histogram mo dification flo w in tro duced by Sapiro and Caselles [45], w e see that their flow can also b e translated in to a combination of repulsion among grey-v alues and a barrier function. How ever, in [45] the repulsiv e force is constant, and the barrier function quadratic. Thus, they cannot b e derived from the same kind of in teraction b et ween the v i and their reflected coun terparts as in our pap er. Referring to Section 5.1, there also exist well-kno wn approaches whic h aim to en- hance the lo cal image contrast such as adaptiv e histogram equalisation – see [42] and the references therein – or contrast limited adaptive histogram equalisation [56]. The 27 latter technique tries to ov ercome the ov er-amplification of noise in mostly homoge- neous image regions when using adaptive histogram equalisation. Both approaches share the basic idea with our approach in Section 5.1 and p erform histogram equal- isation for eac h pixel, i.e. the mapping function for ev ery pixel is determined using a neigh b ourho od of predefined size and its corresp onding histogram. Another related researc h topic is the ric h field of colour image enhancement whic h w e broac h in Section 5.2. A short review of existing metho ds – as well as t wo new ideas – is presen ted in [3]. Therein, Bassiou and K otrop oulus also men tion the colour gam ut problem for metho ds whic h p erform con trast enhancement in a differen t colour space and transform colour co ordinates to R GB afterwards. Of particular in terest are the publications by Naik and Murthy [32] and Nikolo v a and Steidl [34] whose ideas are used in Section 5.2. Both of them suggest – based on an affine colour transform – strategies to ov ercome the colour gamut problem while av oiding colour artefacts in the resulting image. A recen t approach which also makes use of these ideas is presented b y Tian and Cohen [52]. Ojo et al. [37] mak e use of the HSV colour space to av oid the colour gamut problem when enhancing the contrast of colour images. A v ariational approach for contrast enhancement which tries to appro ximate the hue of the input image was recently published b y Pierre et al. [41]. 6 Conclusions and Outlo ok In our pap er we ha ve presen ted a mathematical mo del which describ es pure back- w ard diffusion as gradient descent of strictly con vex energies. The underlying evolu- tion makes use of ideas from the area of collective b eha viour and – in terms of the latter – our mo del can b e understo o d as a fully repulsiv e discrete first order sw arm mo del. Not only it is surprising that our mo del allo ws bac kward diffusion to b e form ulated as a conv ex optimisation problem but also that it is sufficient to imp ose reflecting b oundary conditions in the diffusion co-domain in order to guarantee sta- bilit y . This strategy is contrary to existing approaches which either assume forward or zero diffusion at extrema or add classical fidelit y terms to av oid instabilities. F ur- thermore, discretisation of our mo del do es not require sophisticated numer ics. W e ha ve prov en that a straigh tforward explicit sc heme is sufficien t to preserv e the stabil- it y of the time-con tinuous ev olution. In our exp erimen ts, w e sho w that our mo del can directly b e applied to contrast enhancement of digital greyscale and colour images. W e see our con tribution mainly as an example of stable mo delling of backw ard parab olic ev olutions that create neither theoretical nor n umerical problems. W e are convinced that this concept has far more widespread applications in in v erse problems, image pro cessing, and computer vision. Exploring them will b e part of our future researc h. A ckno wledgements Our research on intrinsic stability prop erties of sp ecific backw ard diffus ion pro cesses go es back to 2005, when Joachim W eic kert was visiting Mila Nik olov a in Paris. Un- fortunately the diffusion equation considered at that time did not ha ve the exp ected prop erties, such that it to ok one more decade to come up with more appropriate 28 mo dels. Leif Bergerhoff wishes to thank Antoine Gautier and Peter Oc hs for fruitful discussions around the topic of the pap er. This pro ject has receiv ed funding from the DF G Cluster of Excellence Multimo dal Computing and Inter action and from the Eu- rop ean Union’s Horizon 2020 researc h and innov ation programme (gran t agreemen t No 741215, ER C Adv anced Gran t INCOVID ). This is gratefully ackno wledged. References 1. Ko dak Lossless True Color Image Suite. http://www.r0k.us/gr aphics/ko dak/ , last visited August 31, 2018 2. Arb elaez, P ., Maire, M., F owlk es, C., Malik, J.: Contour detection and hierarchical image segmentation. IEEE T ransactions on Pattern Analysis and Machine In telligence 33 (5), 898–916 (2011) 3. Bassiou, N., Kotropoulos, C.: Color image histogram equalization by absolute discounting back-off. Computer Vision and Image Understanding 107 (1), 108–122 (Jul 2007) 4. Bergerhoff, L., Cárdenas, M., W eick ert, J., W elk, M.: Mo delling stable backw ard diffusion and repul- siv e sw arms with conv ex energies and range constrain ts. In: Pelillo, M., Hanco c k, E. (eds.) Energy Minimization Metho ds in Computer Vision and Pattern Recognition, 11th International Conference, EMMCVPR 2017, V enice, Italy . Lecture Notes in Computer Science, v ol. 10746, pp. 409–423. Springer, Cham (Mar 2018) 5. Bergerhoff, L., W eic kert, J.: Modelling image pro cessing with discrete first-order swarms. In: Pilla y , N., Engelbrec h t, P .A., Abraham, A., du Plessis, C.M., Snášel, V., Muda, K.A. (eds.) A dv ances in Nature and Biologically Inspired Computing, vol. 419, pp. 261–270. Springer, Cham (2016) 6. Carasso, A., Sanderson, J., Hyman, J.: Digital remov al of random media image degradations by solving the diffusion equation backw ards in time. SIAM Journal on Numerical Analysis 15 (2), 344–367 (Apr 1978) 7. Carasso, A.S.: Comp ensating op erators and stable backw ard in time marching in nonlinear parabolic equations. GEM - International Journal on Geomathematics 5 (1), 1–16 (Apr 2014) 8. Carasso, A.S.: Stable explicit time marching in well-posed or ill-p osed nonlinear parab olic equations. In v erse Problems in Science and Engineering 24 (8), 1364–1384 (2016) 9. Carasso, A.S.: Stabilized Richardson leapfrog scheme in explicit stepwise computation of forw ard or bac kw ard nonlinear parab olic equations. Inv erse Problems in Science and Engineering 25 (12), 1719– 1742 (2017) 10. Carrillo, J.A., F ornasier, M., T oscani, G., V ecil, F.: Particle, kinetic, and hydrodynamic mo dels of sw arming. In: Naldi, G., P areschi, L., T oscani, G. (eds.) Mathematical Modeling of Collective Beha vior in So cio-Economic and Life Sciences, pp. 297–336. Birkhäuser, Boston (2010) 11. Ch uang, Y.L., Huang, Y.R., D’Orsogna, M.R., Bertozzi, A.L.: Multi-vehicle flo c king: Scalability of co operative control algorithms using pairwise p oten tials. In: 2007 IEEE International Conference on Rob otics and Automation, ICRA 2007, 10-14 April 2007, Roma, Italy . pp. 2292–2299 (2007) 12. D’Orsogna, M.R., Chuang, Y.L., Bertozzi, A.L., Chay es, L.S.: Self-prop elled particles with soft-core in teractions: P atterns, stabilit y , and collapse. Physical Review Letters 96 , 104302 (Mar 2006) 13. F u, C.L., Xiong, X.T., Qian, Z.: F ourier regularization for a bac kward heat equation. Journal of Mathematical Analysis and Applications 331 (1), 472–480 (2007) 14. Gab or, D.: Information theory in electron microscopy . Laboratory In vestigation 14 , 801–807 (Jun 1965) 15. Gazi, V., Fidan, B.: Co ordination and control of m ulti-agent dynamic systems: Mo dels and approac hes. In: Sahin, E., Sp ears, W.M., Winfield, A.F.T. (eds.) Swarm Rob otics, Lecture Notes in Computer Science, v ol. 4433, pp. 71–102. Springer, Berlin (2007) 16. Gazi, V., Passino, K.M.: Stabilit y analysis of sw arms. IEEE T ransactions on Automatic Con trol 48 (4), 692–697 (Apr 2003) 17. Gazi, V., Passino, K.M.: Stability analysis of so cial foraging swarms. IEEE T ransactions on Systems, Man, and Cyb ernetics, Part B: Cybernetics 34 (1), 539–557 (F eb 2004) 18. Gersc hgorin, S.: Üb er die Abgrenzung der Eigenw erte einer Matrix. Bulletin de l’A cadémie des Sciences de l’URSS. Classe des sciences mathématiques et naturelles pp. 749–754 (1931) 19. Gilb oa, G., So c hen, N.A., Zeevi, Y.Y.: F orw ard-and-bac kward diffusion pro cesses for adaptive image enhancemen t and denoising. IEEE T ransactions on Image Pro cessing 11 (7), 689–703 (2002) 20. Gonzalez, R., W o ods, R.: Digital Image Pro cessing. Pren tice-Hall, Upper Saddle River, NJ, USA, third edn. (2008) 21. Gron w all, T.H.: Note on the deriv atives with resp ect to a parameter of the solutions of a system of differen tial equations. Annals of Mathematics 20 (4), 292–296 (Jul 1919) 29 22. Hadamard, J.: Sur les problèmes aux dériv ées partielles et leur signification physique. Princeton Uni- v ersit y Bulletin 13 , 49–52 (1902) 23. Hào, D.N., Duc, N.V.: Stabilit y results for the heat equation backw ard in time. Journal of Mathematical Analysis and Applications 353 (2), 627–641 (May 2009) 24. Hào, D.N., Duc, N.V.: Stability results for backw ard parab olic equations with time-dep enden t co effi- cien ts. In v erse Problems 27 (2), 025003 (Jan 2011) 25. Hummel, R.A., Kimia, B.B., Zuc k er, S.W.: Deblurring Gaussian blur. Computer Vision, Graphics, and Image Pro cessing 38 (1), 66–80 (Apr 1987) 26. John, F.: Numerical solution of the equation of heat conduction for preceding times. Annali di Matem- atica Pura ed Applicata 40 (1), 129–142 (Dec 1955) 27. Kirkup, S., W adsworth, M.: Solution of in verse diffusion problems b y op erator-splitting metho ds. Applied Mathematical Mo delling 26 (10), 1003–1018 (2002) 28. K o v ászna y , L.S.G., Joseph, H.M.: Image processing. Pro ceedings of the IRE 43 (5), 560–570 (Ma y 1955) 29. Linden baum, M., Fisc her, M., Bruckstein, A.M.: On Gabor’s con tribution to image enhancemen t. P attern Recognition 27 (1), 1–8 (Jan 1994) 30. Ly apuno v, A.M.: The general problem of the stability of motion. International Journal of Control 55 (3), 531–534 (1992) 31. Mair, B.A., Wilson, D.C., Reti, Z.: Deblurring the discrete Gaussian blur. In: Proceedings of the W orkshop on Mathematical Methods in Biomedical Image Analysis. pp. 273–277. IEEE, Los Alamitos (Jun 1996) 32. Naik, S.K., Murthy , C.A.: Hue-preserving color image enhancement without gamut problem. IEEE T ransactions on Image Pro cessing 12 (12), 1591–1598 (Dec 2003) 33. Nestero v, Y.: Introductory Lectures On Conv ex Optimization. Springer, New Y ork (2004) 34. Nik olo v a, M., Steidl, G.: F ast hue and range preserving histogram specification: Theory and new algorithms for color image enhancement. IEEE T ransactions on Image Pro cessing 23 (9), 4087–4100 (Sep 2014) 35. Nik olo v a, M., Steidl, G.: F ast ordering algorithm for exact histogram sp ecification. IEEE T ransactions on Image Pro cessing 23 (12), 5274–5283 (Dec 2014) 36. Nik olo v a, M., W en, Y.W., Chan, R.: Exact histogram sp ecification for digital images using a v ariational approac h. Journal of Mathematical Imaging and Vision 46 (3), 309–325 (Jul 2013) 37. Ojo, J.A., Solomon, I.D., Adeniran, S.A.: Colour-preserving contrast enhancement algorithm for im- ages. In: Chen, L., Kap o or, S., Bhatia, R. (eds.) Emerging T rends and Adv anced T echnologies for Computational Intelligence: Extended and Selected Results from the Science and Information Confer- ence 2015. pp. 207–222. Springer, Cham (2016) 38. Osher, S., Rudin, L.: Shocks and other nonlinear filtering applied to image processing. In: T escher, A.G. (ed.) Applications of Digital Image Pro cessing XIV, Pro ceedings of SPIE, v ol. 1567, pp. 414–431. SPIE Press, Bellingham (1991) 39. P erk o, L.: Differen tial Equations and Dynamical Systems. No. 7 in T exts in Applied Mathematics, Springer, third edn. (2001) 40. P erona, P ., Malik, J.: S ca le space and edge detection using anisotropic diffusion. IEEE T ransactions on P attern Analysis and Machine Intelligence 12 , 629–639 (1990) 41. Pierre, F., Aujol, J.F., Bugeau, A., Steidl, G., T a, V.T.: V ariational contrast enhancement of gray-scale and R GB images. Journal of Mathematical Imaging and Vision 57 (1), 99–116 (Jan 2017) 42. Pizer, S.M., Am burn, E.P ., Austin, J.D., Cromartie, R., Geselowitz, A., Greer, T., ter Haar Romeny, B.M., Zimmerman, J.B., Zuiderv eld, K.: A daptiv e histogram equalization and its v ariations. Computer Vision, Graphics, and Image Pro cessing 39 (3), 355–368 (1987) 43. P ollak, I., Willsky , A.S., Krim, H.: Image segmentation and edge enhancement with stabilized inv erse diffusion equations. IEEE T ransactions on Image Pro cessing 9 (2), 256–266 (F eb 2000) 44. Pratt, W.K.: Digital Image Pro cessing. Wiley & Sons, New Y ork, third edn. (2001) 45. Sapiro, G., Caselles, V.: Histogram modification via differential equations. Journal of Differential Equations 135 , 238–268 (1997) 46. So c hen, N.A., Zeevi, Y.Y.: Resolution enhancement of colored images by in verse diffusion pro cesses. In: Pro c. 1998 IEEE International Conference on Acoustics, Speech and Signal Pro cessing. pp. 2853–2856. Seattle, W A (May 1998) 47. Steiner, A., Kimmel, R., Bruckstein, A.M.: Planar shap e enhancement and exaggeration. Graphical Mo dels and Image Pro cessing 60 (2), 112–124 (Mar 1998) 48. Stok es, M., Anderson, M., Chandrasek ar, S., Motta, R.: A standard default color space for the internet - sR GB (v ersion 1.10). https://www.w3.or g/Gr aphics/Color/sRGB (Nov 1996), last visited August 31, 2018 49. T autenhahn, U., Schröter, T.: On optimal regularization metho ds for the backw ard heat equation. Zeitsc hrift für Analysis und ihre Anw endungen 15 (2), 475–493 (1996) 30 50. ter Haar Romeny, B.M., Florack, L.M.J., de Swart, M., Wilting, J., Viergev er, M.A.: Deblurring Gaussian blur. In: Bo okstein, F.L., Duncan, J.S., Lange, N., Wilson, D.C. (eds.) Mathematical Metho ds in Medical Imaging II I, Proceedings of SPIE, vol. 2299, pp. 139–148. SPIE Press, Bellingham (Jul 1994) 51. T ernat, F., Orellana, O., Daripa, P .: T wo stable metho ds with numerical experiments for solving the bac kw ard heat equation. Applied Numerical Mathematics 61 (2), 266–284 (F eb 2011) 52. Tian, Q.C., Cohen, L.D.: Color consistency for photo collections without gam ut problems. In: Am- saleg, L., Guðmundsson, G., Gurrin, C., Jónsson, B., Satoh, S. (eds.) MultiMedia Modeling: 23rd In ternational Conference, MMM 2017, Reykjavik, Iceland, January 4-6, 2017, Proceedings, Part I. pp. 90–101. Springer, Cham (Jan 2017) 53. W elk, M., Gilboa, G., W eic k ert, J.: Theoretical foundations for discrete forward-and-bac kward diffusion filtering. In: T ai, X.C., Mørk en, K., Lysaker, M., Lie, K.A. (eds.) Scale Space and V ariational Methods in Computer Vision, Lecture Notes in Computer Science, vol. 5567, pp. 527–538. Springer, Berlin (2009) 54. W elk, M., W eick ert, J., Gilb oa, G.: A discrete theory and efficien t algorithms for forward-and-bac kw ard diffusion filtering. Journal of Mathematical Imaging and Vision 60 (9), 1399–1426 (Nov 2018) 55. Zhao, Z., Meng, Z.: A mo dified Tikhono v regularization metho d for a backw ard heat equation. In v erse Problems in Science and Engineering 19 (8), 1175–1182 (2011) 56. Zuiderv eld, K.: Contrast Limited Adaptiv e Histogram Equalization. In: Hec kb ert, P .S. (ed.) Graphics Gems IV, pp. 474–485. Academic Press Professional, Inc., San Diego, CA, USA (1994)