EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 3, Article Number 6192
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

A Lattice Boltzmann Method for Image Inpainting
Inspired by Fluid Dynamics

Yassine Douich1,∗, Hassan Silkan1, Youssef Hanyf2

1 Department of Computer Science, Laboratory LAROSERI, Faculty of Sciences, Chouaib
Doukkali University, El Jadida, Morocco
2 Research Laboratory in Management and Decision Support, AI Data SEED Team,
Ibn Zohr University, Dakhla, Morocco

Abstract. This paper presents a novel image inpainting algorithm based on the Lattice Boltzmann
Method (LBM), inspired by the vorticity-stream formulation of the Navier-Stokes equations. The
proposed method, employing D2Q5 and D2Q9 lattice models, effectively reconstructs missing or
damaged image regions by propagating smoothness information along isophote directions. Experi-
mental evaluations demonstrate that LBM achieves superior results in terms of PSNR, SSIM, and
reduced CPU time compared to classical inpainting techniques, such as Total Variation (TV) regu-
larisation and the Bertalmio-Sapiro-Caselles-Ballaster (BSCB) algorithm. The method’s efficiency
in handling complex boundary conditions and preserving image details highlights its robustness
for advanced image restoration tasks.

2020 Mathematics Subject Classifications: 65N20, 76-XX

Key Words and Phrases: LBM Inpainting, Navier-Stokes Inpainting, Image Reconstruction,
Edge-Preserving Diffusion, Vorticity-Stream Function

1. Introduction

Image inpainting plays a crucial role in image restoration, aiming to accurately and
seamlessly fill in missing or damaged regions of an image. It is widely used for tasks
such as removing unwanted objects, repairing degraded images, and correcting issues like
cracks, scratches, stains, and red-eye effects.

The lattice boltzmann method (LBM), grounded in the principles of microscopic
statistical physics, has emerged as an effective approach for simulating fluid flows and
capturing complex fluid behavior. Unlike traditional fluid dynamics solvers, LBM offers a
more straightforward programming model by simulating particle interactions and deriving
partial differential equations (PDEs) through moment operations on particle distributions.
Although LBM’s application in image processing is relatively recent, its precursor, cellular

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6192

Email addresses: y.douich@ucd.ac.ma (Y. Douich)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 2 of 32

automata, has been extensively applied in this field. Notable examples include studies by
Preston et al. [1], Hernandez et al. [2], Popovici et al. [3], Wongthanavasu et al. [4],
and Rosin [5], which primarily target binary image processing due to the challenges in
designing cellular automata rules for 256 gray-level images. LBM, however, does not
face this restriction, and it has been used to solve equations such as the diffusion and
diffusion-convection equations, simplifying rule design for processing images with 256 gray
levels.

The lattice Boltzmann method (LBM) has shown great promise in image processing.
Jawerth et al. (1999) applied LBM to anisotropic diffusion using the Perona–Malik model
[6], followed by our 2008 study [6, 7], and Chang’s work on denoising [8]. In this paper,
we propose a new LBM framework for image inpainting based on the vorticity–stream
formulation of the Navier–Stokes equations, treating the image as a fluid medium.

Our method couples two LBM schemes for vorticity transport and intensity evolution,
incorporates adaptive anisotropic diffusion for edge preservation, and uses non-equilibrium
boundary conditions to handle complex contours.

The main contribution is a fast, robust LBM scheme tailored for inpainting, delivering
coherent and high-quality restorations beyond traditional methods.

This work is structured as follows: Section 2 reviews related work in image inpainting,
categorizing methods into PDE-based and learning-based approaches while highlighting
gaps in fluid dynamics-inspired techniques. Section 3 formalizes the mathematical foun-
dation, presenting the vorticity-stream formulation of Navier-Stokes equations adapted
for inpainting. Section 4 details the proposed Lattice Boltzmann Method (LBM) frame-
work, including D2Q5/D2Q9 implementations and boundary condition treatments. Sec-
tion 5 evaluates performance through quantitative metrics (PSNR/SSIM) and comparative
benchmarks against TV and BSCB methods. Section 6 concludes with contributions and
outlines future research directions, notably LBM-CNN hybridization.

2. Related Work for Image Inpainting

In this section, we examine two main categories of research in this field: non-learning-
based methods and learning-based methods.

2.1. Non-Learning-Based Methods

Non-learning-based approaches do not require training data and are typically grounded
in mathematical models. Examples include partial differential equations (PDEs) for
isophote propagation [9], computational fluid dynamics-based isophote propagation [10],
and total variation (TV) minimisation [11] and more traditional methods in [12] . These
foundational techniques aim to preserve the structural integrity of the image during the
inpainting process.

A notable subset of these methods, known as diffusion-based approaches, incorpo-
rates concepts such as intra-channel and inter-channel local variances [13], the distance
and direction between the damaged pixel and its neighbouring pixels [14], and the applica-



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 3 of 32

tion of fractional-order derivatives combined with Fourier transforms [15]. These methods
excel in handling diverse inpainting challenges while maintaining computational efficiency.

2.2. Learning-Based Methods

In recent years, deep convolutional neural networks (CNNs) have demonstrated re-
markable success in various computer vision tasks, including image inpainting. Unlike
non-learning-based methods, CNN-based approaches leverage large-scale datasets and ad-
vanced architectures to achieve high-quality results. These methods are particularly effec-
tive for texture reconstruction and ensuring global structural consistency.

Several specialised CNN architectures have been developed to address specific inpaint-
ing challenges, including U-Net [16], Shift-Net, PEN-Net [17], and BICNN [18]. These
architectures are designed for tasks such as texture completion, fine-detail restoration, and
filling large gaps. Additional applications include inpainting X-ray medical images [19],
forensic image restoration [20], video inpainting for object removal [21], and scanned data
recovery [22]. Techniques such as coherent semantic attention layers [23] and pyramidal-
context architectures [24] further enhance the capabilities of CNN-based methods.

Furthermore, specialised models have been introduced for unique applications. For in-
stance, Artist-Net [25] focuses on high-quality artistic restoration, VGG-based approaches
are tailored for UAV data recovery [26], and GRNN-based methods are effective for non-
texture inpainting [27]. Other applications include text removal [28], artwork enhancement
[29], and semantic object removal [30], showcasing the versatility of CNN-based methods.
Emerging techniques, such as blind inpainting [18] and advanced encoder-decoder archi-
tectures that integrate super-resolution, denoising, and inpainting capabilities [31, 32],
continue to push the boundaries of image restoration.

The development and evaluation of image inpainting algorithms rely heavily on publicly
available and large-scale datasets. The choice of dataset categories significantly impacts
the effectiveness of the proposed methods. Common categories include natural images,
artificial images, face images, and others.

This work highlights some of the most widely used datasets for image inpainting, such
as Paris StreetView [33], Places [34], depth image datasets [35], Foreground-aware datasets
[36], Berkeley segmentation [37], and ImageNet [38], among others.

3. Mathematical Formulation

3.1. Navier-Stokes Equations

Building upon the ideas presented in [10], we consider the incompressible, Newtonian
fluid flow governed by the Navier-Stokes equations:

∂u⃗

∂t
+ (u⃗ · ∇)u⃗+∇P = ν∆u⃗, ∇ · u⃗ = 0

where u⃗ is the velocity vector, P is the pressure, and ν is the viscosity. For two-dimensional
flows, we introduce a stream function Ψ, such that u⃗ = ∇⊥Ψ. Defining the vorticity



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 4 of 32

ω = ∇ × u⃗, we obtain the vorticity-stream function formulation of the Navier-Stokes
equations:

ωt + u⃗ · ∇ω = ν∆ω (1)

In the limit of negligible viscosity, i.e., when ν ≈ 0, these equations reduce to the Euler
equations for inviscid flow. The steady-state solution of Eq. (1) then approaches:

u⃗ · ∇ω = ∇⊥Ψ · ∇∆Ψ ≈ 0 (2)

Note that, in two dimensions, the vorticity ω becomes a scalar quantity related to the
stream function via the Laplace operator: ∆Ψ = ω. The inpainting problem in [9] focuses
on propagating smoothness information ∆I in the direction of the isophotes ∇⊥I (see 1).
Let I0 denote the original image, and let Ω represent a region within the image where
inpainting is desired based on the surrounding data. For an image intensity function
I : Ω→ R, the grayscale values within Ω are determined by the smoothness criterion within
Ω and boundary values on ∂Ω, rather than by restoration methods. Mathematically, this
inpainting criterion is expressed as:

∇⊥I · ∇∆I = 0 on Ω (3)

with boundary condition

I = I0 on ∂Ω (4)

A generalization of Eq.(3) was proposed in [10], introducing an anisotropic diffusion
term in the image I. This generalization is formulated as:

It = ∇⊥I · ∇∆I + ν∇ · (g(|∇I|)∇I) (5)

where g allows for anisotropic diffusion, preserving edges while diffusing smoothness
∆I.

3.2. Inpainting Analogy

The notable resemblance between Eq.(2) and the inpainting criterion Eq.(3) suggests
an analogy between 2D incompressible fluid flow and image inpainting. This analogy is
summarized in the table 1, where Ψ in fluid flow is analogous to the image matrix I in
inpainting.

In image processing terms, we now present an analogue to the vorticity-stream function
formulation (1):

∂ω

∂t
+ u⃗.∇ω = ν∇.(g(|∇ω|)∇ω) on Ω, (6)

where g provides anisotropic diffusion of the smoothness ω or enables edge-preserving
diffusion.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 5 of 32

Table 1: Analogy between Fluid Dynamics and Image Inpainting
c

Aspect Fluid Dynamics Image Inpainting

Physical Quantity Stream function (Ψ) Image intensity (I)

Directional Flow Fluid velocity (u = ∇⊥Ψ) Isophote direction (∇⊥I)

Smoothness Operator Vorticity (ω = ∆Ψ) Smoothness (ω = ∆I)

Diffusion Coefficient Viscosity (ν) Anisotropic diffusion (ν)

The image intensity I, which determines the velocity field u⃗ = ∇⊥I, can be recovered
by simultaneously solving the Poisson problem:

∆I = ω, on Ω, I|∂Ω = I0. (7)

For fluid problems with low viscosity ν, stabilising this formulation can be time-
consuming, which may restrict its practical applicability. A viable alternative is to adopt
a pseudo-steady approach, replacing the Poisson equation (7) with a dynamic relaxation
method based on the heat equation:

∂I

∂t
− α(∆I + ω) = 0, α > 0, I|∂Ω = I0, (8)

where α, the thermal diffusivity, serves as a relaxation parameter that regulates the rate
of convergence.

3.3. Boundary Conditions

In the Navier-Stokes-based inpainting method, continuity across the boundary is in-
herently preserved. For instance, when the boundary condition is defined by setting the
velocity field as u⃗ = ∇⊥I along ∂Ω, solving the inpainting equation 6 under this con-
straint results not only in continuous isophotes but also in an image intensity function
that remains smooth across the boundary ∂Ω. Moreover, a second boundary condition of
Dirichlet type is applied for the heat equation 8. Both types of boundary conditions are
treated at the microscopic level using distribution functions, as described in subsection
4.4.

In the subsequent section, we rigorously recover the vorticity–stream function formu-
lation (6) and the parabolic equation (8) via the Chapman–Enskog asymptotic expansion,
utilizing the lattice Boltzmann models {D2Qm}m=5,9 (refer to Figure 3).

4. Lattice Boltzmann Method Framework

The Lattice Boltzmann Method (LBM) operates on a discretized phase space defined
over a uniform DnQm lattice, where m denotes the discrete velocity directions in an n-
dimensional spatial domain. For two-dimensional problems, the D2Q5 and D2Q9 lattice



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 6 of 32

Figure 1: Illustration of the inpainting mechanism where the smoothness information ∆I
propagates along the isophote directions ∇⊥I, from the known boundary ∂Ω toward the
interior of the missing region Ω.

configurations are employed, as they are well-suited for modeling advection-diffusion and
heat transfer phenomena. A summary of these discrete velocity models is presented in
Table 3.

4.1. Macroscopic Equations Derived via Chapman–Enskog Expansion

To derive the macroscopic advection-diffusion equation for vorticity from the meso-
scopic lattice Boltzmann formulation, we begin by considering the target form:

∂ω(x, t)

∂t
+ u⃗ · ∇ω(x, t) = ν∇ · (k∇ω(x, t)), (9)

where ω(x, t) represents the vorticity field, u⃗ is the velocity, ν is the kinematic viscosity,
and k(x) is a positive, spatially varying diffusion coefficient.

The underlying evolution of the distribution function fi in the lattice Boltzmann
method is governed by:

fi(x+ ci∆t, t+∆t) = fi(x, t)−
1

τ
(fi(x, t)− f

(eq)
i (x, t)), i = 0, . . . ,m− 1, (10)

which can be decomposed into two separate steps: Collision:

f∗
i (x, t) = fi(x, t)−

1

τ
(fi(x, t)− f

(eq)
i (x, t)), (11)

Propagation:
fi(x+ ci∆t, t+∆t) = f∗

i (x, t). (12)



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 7 of 32

Here, fi and f∗
i denote the pre- and post-collision distribution functions, respectively,

ci is the discrete lattice velocity in direction i, and τ is the relaxation time.
The equilibrium distribution function is defined as:

f
(eq)
i = δiω(x, t)

(
1 +

ci · u⃗
c2s

)
, (13)

where δi denotes the lattice weights and cs is the lattice sound speed.
The weights δi are given as follows:

• D2Q5 model: δi =
1
5 , i = 0, 1, . . . , 4

• D2Q9 model:

δi =


4
9 i = 0,
1
9 i = 1, 2, 3, 4,
1
36 i = 5, 6, 7, 8.

(14)

The discrete velocities ci are defined as:

• D2Q5 model:

ci =

{
(0, 0) i = 0,√
2c
(
cos
[
(2i− 1)π4

]
, sin

[
(2i− 1)π4

])
i = 1, 2, 3, 4.

(15)

• D2Q9 model:

ci =


(0, 0) i = 0,

c
(
cos
[
(i− 1)π2

]
, sin

[
(i− 1)π2

])
i = 1, 2, 3, 4,√

2c
(
cos
[
(2i− 9)π4

]
, sin

[
(2i− 9)π4

])
i = 5, 6, 7, 8.

(16)

We also enforce the following moment relations:∑
i

fi =
∑
i

f
(eq)
i = ω(x, t), (17)

∑
i

cif
(eq)
i = ω(x, t)u⃗, (18)

∑
i

cicif
(eq)
i = c2sω(x, t)I, (19)

where I denotes the identity matrix and c2s =
2
5c

2.
To obtain the macroscopic behavior, we employ the Chapman–Enskog expansion:

fi = f
(0)
i + εf

(1)
i + ε2f

(2)
i , ∂t = ε∂t1 + ε2∂t2 , ∇ = ε∇1, (20)



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 8 of 32

with f
(0)
i = f

(eq)
i .

Applying a Taylor expansion to Eq. (10) yields:

fi(x+ ci∆t, t+∆t)− fi(x, t) =
2∑

n=1

∆tn

n!
(∂t + ci · ∇)nfi(x, t) +O(∆t3). (21)

Substituting Eqs. (20) and (21) into Eq. (10) and grouping terms by orders of ε, we
obtain:

∂t1f
(0)
i + ci · ∇1f

(0)
i = − 1

τ∆t
f
(1)
i , (22)

∂t2f
(0)
i + (∂t1 + ci · ∇1)f

(1)
i +

∆t

2
(∂t1 + ci · ∇1)

2f
(0)
i = − 1

τ∆t
f
(2)
i . (23)

Defining D1i = ∂t1 + ci · ∇1, Eq. (22) becomes:

D1if
(0)
i = − 1

τ∆t
f
(1)
i . (24)

Substituting Eq. (24) into Eq. (23) yields:

∂t2f
(0)
i +

(
1− 1

2τ

)
D1if

(1)
i = − 1

τ∆t
f
(2)
i . (25)

Summing Eq. (25) over i, and using Eq. (18), we derive:

∂t2ω +∇1 ·

[(
1− 1

2τ

)∑
i

cif
(1)
i

]
= 0. (26)

Applying Eq. (24) and using moment properties, we obtain:∑
i

cif
(1)
i = −τ∆t∂t1ωu⃗− τ∆t∇1(c

2
sω). (27)

Combining the results, the equations at successive time scales become:

∂t1ω +∇1 · (ωu⃗) = 0, (28)

∂t2ω = ∇1 ·
[
(τ − 1

2
)∆t∂t1ωu⃗

]
+∇1 ·

[
c2s(τ −

1

2
)∆t∇1ω

]
. (29)

Combining both contributions and setting νk = c2s(τ − 1
2)∆t, we recover:

∂tω + u⃗ · ∇ω = ν∇ · (k∇ω) + (τ − 1

2
)∆t∇ · (∂tωu⃗). (30)

Since the flow is incompressible (∇ · u⃗ = 0), the convective term simplifies as:

∇ · (ωu⃗) = u⃗ · ∇ω. (31)



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 9 of 32

Thus, the recovered advection equation reads:

∂tω + u⃗ · ∇ω = ν∇ · (k∇ω) + (τ − 1

2
)∆t∇ · (∂tωu⃗)︸ ︷︷ ︸

additional term

. (32)

In the advection-diffusion equation provided in (32), only terms of order O(u) are
considered. An additional term, (τ − 1

2)∆t∇ · (∂tωu⃗), appears in the formulation. For
incompressible flows, where the velocity field satisfies the condition ∇ · u⃗ = 0, this term
simplifies further to (τ− 1

2)∆t ∂t(∇·(ωu⃗)). Under steady-state conditions in incompressible
flows, this term naturally becomes zero [39].

In the case of unsteady incompressible flows, the magnitude of this additional term
scales as O(u2

c2s
). Given that simulations of incompressible flows are typically conducted

with ∥u⃗∥
cs

< 0.1, the term represents a higher-order correction that contributes negligibly
to the overall dynamics. As such, it is commonly neglected without significant impact on
the accuracy of the simulation [39].

4.2. Recovered Heat Equation

To derive the macroscopic heat equation from the underlying lattice Boltzmann for-
mulation, we consider the mesoscopic evolution equation for the temperature-related dis-
tribution function gi(x, t):

gi(x+ ci∆t, t+∆t) = gi(x, t)−
1

τg
(gi(x, t)− g

(eq)
i (x, t))+α∆tSi, i = 0, . . . ,m− 1, (33)

where τg is the relaxation time, α a scaling coefficient, and Si = ϖiω(x, t) is a source term
associated with the vorticity ω(x, t). The equilibrium distribution function is taken as:

g
(eq)
i (x, t) = ϖiI(x, t), i = 0, . . . ,m− 1, (34)

with I(x, t) denoting the image intensity or temperature field and ϖi the lattice weights.
The moments of the distribution function satisfy the following relations:

m−1∑
i=0

gi(x, t) =
m−1∑
i=0

g
(eq)
i (x, t) = I(x, t), (35)

m−1∑
i=0

cig
(eq)
i (x, t) = 0, (36)

m−1∑
i=0

cicig
(eq)
i (x, t) = c2sI(x, t)I2, (37)

where I2 is the identity matrix in two dimensions.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 10 of 32

To derive the macroscopic behavior, we use a Chapman–Enskog multiscale expansion:

gi = g
(0)
i + εg

(1)
i + ε2g

(2)
i , ∂t = ε∂t1 + ε2∂t2 , ∇ = ε∇1, Si = ε2S

(2)
i . (38)

Expanding Eq. (33) using Taylor series:

gi(x+ ci∆t, t+∆t)− gi(x, t) =

2∑
n=1

∆tn

n!
(∂t + ci · ∇)ngi(x, t) +O(∆t3). (39)

Substituting Eqs. (38) and (39) into Eq. (33) and collecting terms by ε-order yields:

∂t1g
(0)
i + ci · ∇1g

(0)
i = − 1

τg∆t
g
(1)
i , (40)

∂t2g
(0)
i + (∂t1 + ci · ∇1)g

(1)
i +

∆t

2
(∂t1 + ci · ∇1)

2g
(0)
i = − 1

τg∆t
g
(2)
i + S

(2)
i . (41)

Using the operator D1i = ∂t1 + ci · ∇1, Eq. (40) becomes:

D1ig
(0)
i = − 1

τg∆t
g
(1)
i . (42)

Substituting into Eq. (41), we get:

∂t2g
(0)
i +

(
1− 1

2τg

)
D2

1ig
(0)
i = − 1

τg∆t
g
(2)
i + S

(2)
i . (43)

Summing over all i and using Eq. (35) and properties of the equilibrium distribution,
we obtain:

∂t1I = 0, (44)

∂t2I +

(
1− 1

2τg

)
∇1 ·

(
m−1∑
i=0

cig
(1)
i

)
= S(2), (45)

where S(2) =
∑
i

S
(2)
i .

To evaluate the flux term, we use Eq. (42):∑
i

cig
(1)
i = −τg∆t

∑
i

ciD1ig
(eq)
i = −τg∆tc2s∇1I. (46)

Substituting into Eq. (45) gives:

∂t2I = c2s(τg −
1

2
)∆t∇2

1I + S(2). (47)

Combining Eqs. (44) and (47), and setting α = c2s(τg−1
2)∆t, we recover the macroscopic

nonlinear heat equation:
∂tI = α∇2I(x, t) + αω(x, t). (48)



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 11 of 32

4.3. Computation of ∇ω, τ , τg, and ∇⊥I

In the macroscopic equation (6), the coefficient g(s) represents a positive, non-increasing
diffusion function. Inspired by the work of Perona and Malik [40], typical choices for g(s),
where s = ∥∇ω∥, include:

g(s) = e−ls2 , or g(s) =
1

1 + ls2
, (49)

where l > 0 is a tunable parameter. In this work, we adopt the second form:

g(s) =
1

1 + l∥∇ω∥2
. (50)

Under the anisotropic diffusion framework where k = g(∥∇ω∥), the relaxation times τ
and τg are obtained from the macro-micro relations:

νk = c2s

(
τ − 1

2

)
∆t, α = c2s

(
τg −

1

2

)
∆t. (51)

Hence, the expressions for τ and τg are:

τ =
1

2
+

νk

c2s∆t
=


1

2
+

5∆tν

2(∆x)2(1 + l∥∇ω∥2)
for D2Q5 lattice,

1

2
+

3∆tν

(∆x)2(1 + l∥∇ω∥2)
for D2Q9 lattice.

(52)

τg =
1

2
+

α

c2s∆t
=


1

2
+

5∆t

2(∆x)2
α for D2Q5 lattice,

1

2
+

3∆t

(∆x)2
α for D2Q9 lattice.

(53)

To compute τ from Eq. (52), one needs to first evaluate the vorticity gradient ∇ω.
This is computed using a central finite difference scheme depending on the lattice model:

∇ω(x, t) =


1

2∆x

4∑
i=0

eiω(x+ ci∆t), for D2Q5 lattice,

1

6∆x

8∑
i=0

eiω(x+ ci∆t), for D2Q9 lattice,

(54)

where ei =
ci
c .

To evaluate the image gradient ∇I, we recall from Eq. (46) that:

∇1I = − 1

τg∆tc2s

m−1∑
i=0

cig
(1)
i . (55)



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 12 of 32

Assuming εg
(1)
i ≈ gi − g

(eq)
i , and substituting from Eq. (34), we obtain:

∇I = − 1

τg∆tc2s

m−1∑
i=0

cigi(x), (56)

where we used the fact that
∑
i

cig
(eq)
i = 0 due to symmetry.

The perpendicular gradient ∇⊥I is then given by:

∇⊥I = − 1

τg∆tc2s

m−1∑
i=0

c⊥i gi(x), x ∈ Ω, (57)

where c⊥i denotes the velocity vector perpendicular to ci.

4.4. Curved Boundary Conditions

In partial differential equation (PDE)-based image inpainting, boundary conditions
play a pivotal role in governing how information propagates into the inpainting domain.
In this work, the region to be inpainted, denoted Ω, is manually selected and often features
a nontrivial, curved boundary ∂Ω that conforms closely to the underlying image structures.

To ensure a coherent and stable reconstruction, two boundary conditions are imposed.
First, a velocity field is prescribed on the boundary by enforcing u⃗ = ∇⊥I|∂Ω, which serves
as a guiding advection flow along ∂Ω. This condition regulates the directional transport
of isophote information into the missing region. Second, a Dirichlet condition I|∂Ω = I0
is applied to the nonlinear heat equation, thereby anchoring the intensity values on the
boundary to those of the known image.

Together, these boundary constraints ensure that the inpainting process yields a smooth
and seamless transition between the reconstructed and known regions, preserving the struc-
tural and tonal coherence of the original image.

4.5. Non-Equilibrium Extrapolation Method for Curved Boundaries

Curved boundary conditions are treated using the non-equilibrium extrapolation tech-
nique. In this approach, the distribution functions at boundary nodes xb are decomposed
into equilibrium and non-equilibrium components:

fī(xb, t) = f
(eq)

ī
(xb, t) + f

(neq)

ī
(xb, t), (58)

where the equilibrium part f
(eq)

ī
(xb, t) is approximated by:

f
(eq)

ī
(xb, t) ≈ f∗

ī (xb, t) ≡ δī ω(xb)

(
1 +

c⃗ī · u⃗(xb)
c2s

)
, i = 0, . . . ,m− 1, (59)



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 13 of 32

Figure 2: Schematic diagram of the node structure used in the Lattice Boltzmann frame-
work. The domain includes fluid nodes (xf , xff ), solid nodes (xs), boundary nodes (xb),
and wall boundary nodes (xw). The lattice spacing is denoted by ∆x, and c⃗i represents
the discrete velocity in the i-th direction, with c⃗ī as its opposite counterpart.

where ω(xb) and u⃗(xb) denote the vorticity function and velocity vector at node xb,
respectively. Since these quantities are generally not known a priori, Guo et al. [41]
proposed the following extrapolation scheme:

f
(neq)

ī
(xb, t) =


fī(xf , t)− f

(eq)

ī
(xf , t), q ≥ 3

4
,

q
(
fī(xf , t)− f

(eq)

ī
(xf , t)

)
+ (1− q)

(
fī(xff , t)− f

(eq)

ī
(xff , t)

)
, q <

3

4
,

(60)

where q =
|xf − xw|
|xf − xb|

=
|xf − xw|

∆x
, and ω(xb) = ω(xf ). The velocity u⃗(xb) is estimated

by:

u⃗(xb) =


1

q
u⃗(xw) +

q − 1

q
u⃗(xf ), q ≥ 3

4
,

3− q

q + 1
u⃗(xw) + (q − 1)u⃗(xf ) + (1− q)

q − 1

q + 1
u⃗(xff ), q <

3

4
.

(61)

In the fluid region Ω, the velocity field u⃗, corresponding to the inpainted region, is
computed from the gradient of the auxiliary function gi as follows:

u⃗(xf ) = −
1

τg∆tc2s

m−1∑
i=0

c⃗⊥i gi(xf ), (62)

u⃗(xff ) = −
1

τg∆tc2s

m−1∑
i=0

c⃗⊥i gi(xff ). (63)

At the wall boundary ∂Ω, the velocity u⃗(xw) is evaluated via a discrete transverse
gradient of the initial image I0:



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 14 of 32

u⃗(xw) = ∇⊥I|∂Ω =


1

2∆x

4∑
i=0

e⃗⊥i I0(xw + c⃗i∆t),

1

6∆x

8∑
i=0

e⃗⊥i I0(xw + c⃗i∆t),

(64)

depending on whether a D2Q5 or D2Q9 stencil is used.
Using the above approximations, the post-collision distribution function at the bound-

ary node xb is given by:

f
′

ī (xb, t) = fī(xb, t)−
1

τ
f
(neq)

ī
(xb, t) = f∗

ī (xb, t) +
τ − 1

τ
f
(neq)

ī
(xb, t). (65)

For scalar transport with Dirichlet boundary conditions, the distribution function
gī(xb, t) is treated analogously:

gī(xb, t) = g
(eq)

ī
(xb, t) + g

(neq)

ī
(xb, t), (66)

where the equilibrium part is approximated by:

g
(eq)

ī
(xb, t) ≈ g∗ī (xb, t) ≡ ϖiI(xf ), i = 0, . . . ,m− 1, (67)

and the non-equilibrium part is extrapolated as:

g
(neq)

ī
(xb, t) =


gī(xf , t)− g

(eq)

ī
(xf , t), q ≥ 3

4
,

q
(
gī(xf , t)− g

(eq)

ī
(xf , t)

)
+ (1− q)

(
gī(xff , t)− g

(eq)

ī
(xff , t)

)
, q <

3

4
.

(68)
Thus, the post-collision distribution for the scalar field becomes:

g
′

ī(xb, t) = g∗ī (xb, t) +
τg − 1

τg
g
(neq)

ī
(xb, t). (69)

Figure 3: Lattice configurations: (a) D2Q9, (b) D2Q5.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 15 of 32

4.6. implementation of the Algorithm

To facilitate the understanding of the proposed Lattice Boltzmann inpainting frame-
work, all relevant parameters and their respective roles are summarised in the following
two tables. Table 2 lists structural and numerical variables, while Table 3 gathers physical
quantities and functional coefficients involved in the image reconstruction process.

Table 2: Structural and numerical parameters involved in the LBM-based inpainting al-
gorithms.

Symbol Name Description

I0 Initial image Input grayscale or RGB image to be restored
M Inpainting mask Binary matrix: 1 inside the missing region, 0 elsewhere
I Inpainted image Output image updated through iterative reconstruction
fi Vorticity distribution LBM distribution function for vorticity
gi Intensity distribution LBM distribution function for intensity

f
(eq)
i Equilibrium (vorticity) Local equilibrium for ω

g
(eq)
i Equilibrium (intensity) Local equilibrium for I
ci Lattice velocities Discrete propagation directions (D2Q5/D2Q9)
cs Sound speed LBM constant for viscosity/diffusivity scaling
δi, ϖi Lattice weights Weights in equilibrium formulations
N Number of iterations Fixed stopping criterion for simulation
∆t Time step Time increment (set to 1)
∆x, ∆y Grid spacing Spatial resolution (set to 1)

Table 3: Physical and functional parameters governing the behaviour of the LBM inpaint-
ing model.

Symbol Name Description

ω Vorticity field Encodes local image smoothness flow
v = ∇⊥I Guidance field Isophote direction used in transport
Si Source term Forcing term derived from ω for the heat equation
τ Relaxation time for ω Controls diffusion rate via viscosity ν
τg Relaxation time for I Controls thermal behaviour via diffusivity α
ν Kinematic viscosity Regulates the smoothing of vorticity
α Thermal diffusivity Controls speed of convergence in I evolution
k Diffusion function Anisotropic weight for edge preservation



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 16 of 32

Algorithm 1 Proposed Lattice Boltzmann Method Algorithm for Grayscale Image In-
painting

Data: Initial image I0, binary mask M indicating the inpainting region
Result: Reconstructed image I

Initialization:
Define the parameters of lattice DnQm ; // ∆x, ∆t, c = ∆x

∆t , velocities ci,wi

and cs
I ← I0 ; // Copy initial image

gi(x)
(0) ← geqi (x, 0) ; // Initialize g with equilibrium

ω(0) ← 0, f
(0)
i ← 0 ; // Initialize distributions

v ← ∇⊥I ; // Perpendicular gradient (guidance field)

for n = 1 to N do
foreach pixel x where M(x) = 1 do

Compute ω(n)(x) ←
4∑

i=0

f
(n)
i (x) Compute ∇ω(n)(x) using Eq. (54) Compute re-

laxation time τ (n)(x) using Eq. (52)
end
Apply boundary conditions using Eq. (65) and Eq. (69)
foreach pixel x where M(x) = 1 do

f
(n)
i (x)← f

eq(n)
i (x) ; // Equilibrium initialization

Collision step:

f
(n)
i (x) ← f

(n)
i (x) − 1

τ (n)(x)

(
f
(n)
i (x)− f

eq(n)
i (x)

)
g
(n)
i (x) ← g

(n)
i (x) −

1

τg

(
g
(n)
i (x)− g

eq(n)
i (x)

)
+ α∆tϖiω

(n)(x)

Streaming step:

f
(n+1)
i (x+ ci∆t)← f

(n)
i (x) g

(n+1)
i (x+ ci∆t)← g

(n)
i (x)

Reconstruction:

I(n+1)(x)←
4∑

i=0

g
(n+1)
i (x)

end

end

I ← I(N+1)(x) ; // Final inpainted image



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 17 of 32

Algorithm 2 Proposed Lattice Boltzmann Method Algorithm for Inpainting Colored
Images

Input: Initial color image I0 with 3 channels (RGB), mask M for each channel
Output: Inpainted color image I

Initialization:
Define the parameters lattice DnQm ; // ∆x, ∆t, c = ∆x

∆t , velocities ci,wi

and cs
I ← I0 ; // Copy initial image into working variable

gi(x)
(0) ← geqi (x, 0) ; // Initialize auxiliary distribution function with

equilibrium

ω(0) ← 0, f
(0)
i ← 0 ; // Initialize macroscopic and microscopic

distributions

v ← ∇⊥I ; // Compute perpendicular gradient (guidance field)

for k = 1 to 3 // Loop over R, G, B channels do
Apply Algorithm 1 to channel I(:, :, k) ; // Inpaint current channel using LBM

end
return Inpainted color image I ; // Final reconstructed RGB image

5. Experimental results and discussion

The proposed image inpainting algorithm was evaluated using a total of six test images,
including three grayscale images (Cameraman, Clown, Girlface) and three colour images
(Boat, Bridge, El Jadida City). To simulate degradation, random lines were superimposed
on the original images, mimicking scratches or missing regions. The inpainting was per-
formed using the Lattice Boltzmann Method (LBM) based on the D2Q5 and D2Q9 lattice
models, (see Fig.3).

To benchmark the performance of the proposed LBM-based inpainting method, we
compare it against two well-established classical approaches: the Total Variation (TV)
regularisation method and the Bertalmio–Sapiro–Caselles–Ballester (BSCB) algo-
rithm. The TV model, introduced by Rudin, Osher, and Fatemi, represents a variational
framework that effectively preserves edges while smoothing homogeneous regions. Mean-
while, the BSCB method is one of the earliest inpainting techniques inspired by fluid dy-
namics, using Navier–Stokes-based transport along isophotes. These two methods have
become canonical baselines in the literature, with TV exemplifying geometric regu-
larisation and BSCB reflecting the physical modeling perspective. Their inclusion
provides a meaningful contrast for evaluating the reconstruction quality, computa-
tional efficiency, and structural coherence of our LBM framework.

The restored images were compared with those obtained using BSCB and TV ap-
proaches. All implementations were carried out in Matlab 2018 on a Windows 10 operat-
ing system. The quality of the restored images was quantitatively assessed using the Peak
Signal-to-Noise Ratio (PSNR) and the Structural Similarity Index Measure (SSIM) . The
PSNR is a standard metric that evaluates the fidelity of the restored image relative to the



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 18 of 32

original, defined by:

PSNR(t) = 10 log10


mn

(m,n)∑
(i,j)=(1,1)

[I((i, j), t)− I0(i, j)]
2

 ,

where m and n denote the dimensions of the image, I((i, j), t) is the pixel value of the
inpainted image at time t, and I0(i, j) is the pixel value of the original image.

The SSIM metric evaluates the similarity between two images by incorporating per-
ceptual factors such as luminance, contrast, and structural information. It is given by:

SSIM(X,Y ) =

[
2µXµY + C1

µ2
X + µ2

Y + C1

]
·
[

2σXY + C2

σ2
X + σ2

Y + C2

]
·

 σXY√
σ2
Xσ2

Y

 ,

where µX and µY are the means, σ2
X and σ2

Y are the variances, and σXY is the covari-
ance of local image patches from images X and Y . The constants C1 and C2 are stabilizing
terms defined as C1 = (K1L)

2 and C2 = (K2L)
2, with L = 255 representing the dynamic

range of pixel intensities, and K1 = 0.01, K2 = 0.03.
In our numerical experiments, the proposed algorithm was applied to both grayscale

and colour images and compared against two classical inpainting techniques: Total Varia-
tion (TV) regularisation and the Bertalmio–Sapiro–Caselles–Ballester (BSCB) algorithm.
To ensure consistency across all methods, a uniform spatial discretisation was adopted,
with ∆x = ∆y = 1. For the TV-based inpainting approach, the time step ∆t was selected
from the set {0.5, 1}, and the stopping criterion was based on a convergence tolerance
ϵ ∈ {0.01, 0.001, 0.0008}, depending on the required reconstruction accuracy. In contrast,
the BSCB method used a time step ∆t ∈ {0.1, 0.5}, chosen to ensure numerical stability
and structural coherence. For the Lattice Boltzmann Method (LBM), a standard configu-
ration with ∆x = ∆y = 1 and ∆t = 1 was employed, as commonly adopted in LBM-based
schemes. Unlike energy-minimisation-based iterative methods, the stopping condition for
the LBM algorithm was based on a fixed number of iterations Niter, generally ranging
between 100 and 2000 depending on the image resolution and the extent of the inpainting
region. In the LBM framework, the parameters ν and α play crucial roles in controlling
the diffusion dynamics. The kinematic viscosity ν, which governs the spread of vorticity in
the advection-diffusion equation, depends on both the relaxation time τ and the spatially
varying diffusion function k, which is often designed to preserve edge information through
anisotropic smoothing. On the other hand, the thermal diffusivityα, which regulates the
temporal evolution of the image intensity field, is determined by the relaxation time τg
in the thermal LBM scheme. For numerical stability, both relaxation times must satisfy
the condition τ > 1

2 and τg > 1
2 , which ensures the positivity of the physical diffusion

coefficients and prevents unphysical oscillations during the simulation. Proper calibra-



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 19 of 32

tion of these parameters is essential to achieve a balance between fast convergence, detail
preservation, and stable inpainting dynamics.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 20 of 32

(a) cameraman

(b) clown

(c) girlface

(d) boat

(e) bridge

(f) El Jadida city

Figure 4: The grid of 3 × 2 masked images includes grayscale images in the left column
and color images in the right column.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 21 of 32

(a) LBM-D2Q5(α = 0.88, ν = 2.5) (b) LBM-D2Q9(α = 0.43, ν = 2)

(c) BSCB Method: ∆t = 0.1 (d) TV Regularization :ϵ = 10−3, ∆t =
1

Figure 5: Inpainting comparison on the Cameraman grayscale image corrupted by a syn-
thetic random line simulating a scratch. The missing region was restored using four
different techniques. (a)–(b): reconstructions via the proposed LBM method with D2Q5
and D2Q9 lattices; (c): BSCB algorithm; (d): Total Variation (TV) regularisation. LBM-
based methods show sharper edge continuity and better restoration of the corrupted region
while requiring significantly less computation time.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 22 of 32

(a) LBM-D2Q5 (α = 0.16, ν = 2.7) (b) LBM-D2Q9(α = 0.23, ν = 0.9)

(c) BSCB method: ∆t = 0.1 (d) TV regularization ϵ = 10−3, ∆t = 1

Figure 6: Inpainting results for the Girlface grayscale image artificially damaged by ran-
domly oriented lines. The missing regions are reconstructed using four different algorithms.
(a)–(b) show the results obtained with the proposed LBM approach using D2Q5 and D2Q9
lattices. (c) and (d) correspond to the BSCB and TV regularisation methods, respectively.
The LBM-D2Q9 model offers enhanced restoration quality, particularly in facial contours,
with minimal blurring and improved structural consistency.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 23 of 32

(a) LBM-D2Q5(α = 0.36, ν = 1) (b) LBM-D2Q9 (α = 0.43, ν = 2.7)

(c) BSCB Method: ∆t = 0.1 (d) TV Regularization: ϵ = 8.10−4 and
∆t = 1

Figure 7: Inpainting results for the Clown grayscale image artificially damaged by a
magenta-colored line simulating a severe scratch. (a)–(b): reconstructions using the pro-
posed LBM method with D2Q5 and D2Q9 lattices; (c)–(d): results from BSCB and TV-
based techniques. LBM approaches, especially D2Q9, exhibit superior capability in restor-
ing textures and curved contours with minimal artifacts.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 24 of 32

(a) LBM-D2Q5(α = 0.44 and ν = 2) (b) LBM-D2Q9(α = 0.23 and ν = 3)

(c) BSCB Method: ∆t = 0.5 (d) TV Regularization :∆t = 1 and ϵ =
0.01

Figure 8: Object removal results on the El Jadida City colour image. A building structure
was manually masked to simulate a missing region. (a)–(b): inpainting results using the
proposed LBM method with D2Q5 and D2Q9 lattices. (c): BSCB method; (d): Total
Variation (TV) regularisation. The LBM approaches achieve structurally consistent and
visually plausible restorations, particularly in complex regions such as architectural lines
and sky gradients.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 25 of 32

(a) LBM-D2Q5:α = 1.68 and ν = 0.5 (b) LBM-D2Q9:α = 0.8 and ν = 2.5

(c) BSCB method with ∆t = 0.1 (d) TV Regularization : ϵ = 0.01 and
∆t = 0.5

Figure 9: Inpainting results for the Boat color image corrupted by a red-colored line
simulating a structural defect. (a)–(b): Reconstructions using the proposed LBM method
with D2Q5 and D2Q9 lattice models, respectively. (c): Result obtained using the BSCB
algorithm. (d): Output from the TV regularization method. The LBM-D2Q5 model
achieves the highest PSNR, effectively restoring fine details and textures with minimal
artifacts and significantly lower computation time compared to classical approaches.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 26 of 32

(a) LBM-D2Q5: α = 0.48 and ν = 1 (b) LBM-D2Q9: α = 0.13 and ν = 2

(c) BSCB method: ∆t = 0.5 (d) TV Regularization: ϵ = 0.01 and
∆t = 0.5

Figure 10: Object removal results on the Bridge colour image. A structural component
was manually masked to simulate a missing region. (a)–(b): inpainting results using the
proposed LBM algorithm with D2Q5 and D2Q9 lattice models. (c): BSCB method; (d):
Total Variation (TV) regularisation. The LBM-based approaches effectively restore fine
contours and background textures, yielding more visually coherent reconstructions than
classical methods.

The results in Tables 4 and 5 demonstrate the superior performance of the Lattice
Boltzmann Method (LBM) compared to traditional methods like Total Variation regu-
larisation (TV reg) and the Bertalmio-Sapiro-Caselles-Ballaster (BSCB) algorithm. For
grayscale images, LBM achieved comparable or slightly better PSNR and SSIM values
while significantly reducing computational time. For instance, the Cameraman image pro-
cessed with the LBM-D2Q5 model achieved a PSNR of 25.30 dB in just 0.93 minutes,



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 27 of 32

Table 4: Quantitative performance comparison for grayscale image inpainting across dif-
ferent methods. Metrics include PSNR (dB), SSIM, and CPU Time (in minutes) for the
Cameraman, Girlface, and Clown images. The LBM-based approaches (D2Q5 and D2Q9)
consistently yield comparable or superior reconstruction quality while significantly reduc-
ing computational cost compared to classical methods (TV regularization and BSCB).

Image Method PSNR (dB) SSIM CPU Time (min)

4*Cameraman TV reg 25.28 0.999 5.10
BSCB 25.27 0.999 5.54

LBM-D2Q5 25.30 0.999 0.93
LBM-D2Q9 25.31 0.999 1.21

4*Girlface TV reg 15.16 0.9996 10.73
BSCB 15.12 0.9996 9.31

LBM-D2Q5 15.16 0.9996 0.62
LBM-D2Q9 15.17 0.9996 1.60

4*Clown TV reg 16.47 0.9995 11.27
BSCB 16.36 0.9995 16.76

LBM-D2Q5 16.48 0.9995 0.64
LBM-D2Q9 16.48 0.9995 1.12

compared to over 5 minutes for both TV reg and BSCB. This efficiency was also ev-
ident for other images, such as Girlface and Clown, where LBM consistently delivered
high-quality inpainting results in a fraction of the time required by classical methods.

Similarly, for colour images, the LBM methods showcased their ability to handle com-
plex inpainting tasks efficiently. The D2Q5 model stood out, achieving the highest PSNR
values with the shortest processing times. For example, the Boat image inpainted with
LBM-D2Q5 reached a PSNR of 37.31 dB in 1.50 minutes, compared to 14.36 minutes
for TV reg (PSNR: 35.25 dB) and 19.48 minutes for BSCB (PSNR: 33.91 dB). The
D2Q9 model also performed well, particularly for higher relaxation times, highlighting the
adaptability of LBM to different conditions. These results confirm LBM as a robust and
efficient approach for image inpainting, offering significant improvements in computational
speed and quality over traditional methods.

The figure 4 highlights the influence of the thermal diffusivity parameter α on the
temporal evolution of PSNR for the ”Girlface” image over 100 iterations using the Lattice
Boltzmann Method (LBM). Smaller α values enable faster convergence and higher final
PSNR, indicating improved inpainting quality due to controlled diffusion that preserves
image details. Conversely, larger α values result in slower convergence and reduced in-
painting performance, likely due to excessive smoothing that diminishes image fidelity.
This underscores the importance of optimising α for achieving high-quality inpainting
results.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 28 of 32

Table 5: Performance metrics for color image inpainting on three datasets: El Jadida
City, Boat, and Bridge. The LBM-D2Q5 method demonstrates superior PSNR values and
drastically reduced CPU times. Both LBM variants outperform classical methods (TV
reg and BSCB) in terms of efficiency and image fidelity, confirming the robustness of the
proposed framework for high-resolution color image restoration.

Image Method PSNR (dB) SSIM CPU Time (min)

4*City TV reg 27.23 1.000 85.56
BSCB 28.58 1.000 74.90

LBM-D2Q5 29.94 1.000 17.01
LBM-D2Q9 29.85 1.000 26.07

4*Boat TV reg 35.25 1.000 14.36
BSCB 33.91 1.000 19.48

LBM-D2Q5 37.31 1.000 1.50
LBM-D2Q9 36.89 1.000 2.10

4*Bridge TV reg 30.92 1.000 52.20
BSCB 29.97 1.000 55.23

LBM-D2Q5 31.96 1.000 7.73
LBM-D2Q9 31.98 1.000 11.60

Figure 11: Temporal evolution of the Peak Signal-to-Noise Ratio (PSNR) for the Girlface
image over 100 iterations, illustrating the impact of the thermal diffusivity parameter α on
the inpainting performance. Lower α values lead to faster convergence and higher PSNR,
reflecting enhanced reconstruction quality through controlled diffusion. In contrast, higher
α values result in slower convergence and diminished fidelity due to excessive smoothing.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 29 of 32

(a) D2Q5 performs better at α = 0.48, τg = 1.7
for D2Q5; α = 0.40, τg = 1.7 for D2Q9

(b) D2Q9 performs better at α = 2.48, τg = 6.7
for D2Q5; α = 2.06, τg = 6.7 for D2Q9

Figure 12: Comparison of PSNR evolution for the Girlface image using D2Q5 and D2Q9
lattice models over 100 iterations under two thermal diffusivity settings. (a) For moder-
ate diffusivity, D2Q5 shows faster convergence and higher final PSNR, indicating better
preservation of image structure. (b) For higher diffusivity, D2Q9 achieves superior PSNR,
suggesting improved performance in scenarios requiring stronger smoothing.

6. Conclusions and future work

In this study, we introduced a novel image inpainting algorithm grounded in the Lat-
tice Boltzmann Method (LBM), drawing inspiration from fluid dynamics through a vor-
ticity–stream formulation. The method was validated on six benchmark images—three
grayscale and three color—artificially degraded using random scratches. Experimental re-
sults show that the proposed LBM approach outperforms classical techniques such as Total
Variation (TV) regularization and the BSCB algorithm in both reconstruction quality and
computational efficiency.

The LBM-D2Q5 model achieved faster processing times than TV and BSCB, reducing
computation from several minutes to under one minute while preserving or improving
PSNR and SSIM. This confirms its ability to deliver high-quality inpainting with low
computational cost.

Nevertheless, limitations exist. The model may underperform in highly textured re-
gions, where deep learning approaches tend to be more effective. In addition, the method’s
sensitivity to parameters such as α and ν necessitates empirical tuning to achieve optimal
performance.

Future work will explore hybrid frameworks combining LBM with deep neural networks
to synergize the physical interpretability and speed of LBM with the representational
power of deep learning. Such integration could enable adaptive parameter adjustment,
enhanced texture reconstruction, and better generalization across diverse image types and



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 30 of 32

damage scenarios.

References

[1] K. Preston, M.J.B. Duff, S. Levialdi, P.E. Norgren, and J. Toriwaki. Basics of cellular
logic with some applications in medical image processing. Proceedings of the IEEE,
67(5):826–856, 1979.

[2] G. Hernandez and H. J. Herrnann. Cellular automata for elementary image enhance-
ment. Graphical Models and Image Processing, 58:82–89, 1996.

[3] A. Popovici and D. Popovici. Cellular automata in image processing. In: MTNS,
2002.

[4] S. Wongthanavasu and R. Sadananda. A ca-based edged operator and its performance
evaluation. Journal of Visual Communication and Image Representation, 14:83–96,
2003.

[5] SP. L. Rosin. Training cellular automata for image processing. IEEE Transaction on
image processing, 15:2076–2087, 2006.

[6] P. Lin B. Jawerth and E. Sinzinger. Lattice boltzmann models for anisotropic diffusion
of images. Journal of Mathematical Imaging and Vision, pages 231–237, 1999.

[7] Z.Z.Yan Y.Chen and Y.H. Qian. An anisotropic diffusion model for medical im-
age smoothing by using the lattice boltzmann method. IFMBE Proceeding of
APCMBE2008, pages 255–259, 2008.

[8] Q. S. Chang and Y. Tong. A lattice boltzmann method for image denoising. IEEE
Transactions on Image Processing, 18:2797–2802, 2009.

[9] V. CASELLES M. Bertamlio, G. SAPIRO and C. BALLESTER. Image inpaint-
ing. Proceedings of the ACM SIGGRAPH Conference on Computer Graphics, page
417–424, 2000.

[10] A. L. Bertozzi M. Bertalmio and G. Sapiro. Navier-stokes, fluid dynamics, and image
and video inpainting. Proceedings of the 2001 IEEE Computer Society Conference on
Computer Vision and Pattern Recognition. CVPR 2001, 1:I–I, 2001.

[11] T T. Chan. Local inpainting models and tv inpainting. SIAM J. Appl. Math,
62:1019–1043, 2001.

[12] Riya Shah, Anjali Gautam, and Satish Kumar Singh. Overview of image inpainting
techniques: A survey. pages 1–6, 2022.

[13] Haodong Li, Weiqi Luo, and Jiwu Huang. Localization of diffusion-based inpaint-
ing in digital images. IEEE Transactions on Information Forensics and Security,
12(12):3050–3064, 2017.

[14] S. Liu Y. Wei. Domain-based structure-aware image inpainting. Signal, Image and
Video Processing10, 5:911–919, 2016.

[15] G. Yang D. Zhang, Z. Liang, Q. Li, L. Li, and X. Sun. A robust forgery detection
algorithm for object removal by exemplar-based image inpainting 77. Signal, Image
and Video Processing10, 10:11823–11842, 2018.

[16] M. Li Z. Yan, X. Li, W. Zuo, and S. Shan. Shift-net: Image inpainting via deep feature



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 31 of 32

rearrangement. in:Proceedings of the European Conference on Computer Vision, pages
1–17, 2018.

[17] H. Chao Y. Zeng, J. Fu and B. Guo. Learning pyramid-context encoder network for
highquality image inpainting. in: Proceedings of the IEEE Conference on Computer
Visionand Pattern Recognition, pages 1486–1494, 2019.

[18] Z. Lin N. Cai, Z. Su, H. Wang, Z. Yang, and B. W.-K. Ling. Blind inpainting using
thefully convolutional neural network. The Visual Computer 33, 2:249–261, 2017.

[19] Ji Zhao, Zhiqiang Chen, Li Zhang, and Xin Jin. Unsupervised learnable sinogram
inpainting network (sin) for limited angle ct reconstruction, 2018.

[20] X. Zhao X. Zhu, Y. Qian, B. Sun, and Y. Sun. A deep learning approach to patch-
based image inpainting forensicsk. Signal Processing: Image Communication, 67:90–
99, 2018.

[21] Z. Yu Liu Y.-L. Chang, W. Hsu, and Vornet. Spatio-temporally consistent video
inpainting for object removal. in: Proceedings of the IEEE Conference on Computer
Vision and Pattern Recognition Workshops, 2019.

[22] S. Iizuka K. Sasaki, E. Simo-Serra, and H. Ishikawa. Joint gap detection and inpaint-
ing of line drawings. in: Proceedings of the IEEE Conference on Computer Vision
and Pattern Recognition, pages 5725–5733, 2017.

[23] Hongyu Liu, Bin Jiang, Yi Xiao, and Chao Yang. Coherent semantic attention for
image inpainting, 2019.

[24] H. Chao Y. Zeng, J. Fu and B. Guo. Learning pyramid-context encoder network for
highquality image inpainting. in: Proceedings of the IEEE Conference on Computer
Vision and Pattern Recognition, pages 1486–1494, 2019.

[25] J. Xiao L. Liao, R. Hu and Z. Wang. Artist-net: Decorating the inferred content with
unified style for image inpainting. IEEE Access, 7:36921–36933, 2019.

[26] G. Wang C. Hsu, F. Chen. High-resolution image inpainting through multiple deep
networks. 2017International Conference on Vision, Image and Signal Processing
(ICVISP), IEEE, pages 76–81, 2017.

[27] F. Yaghmaee V. K. Alilou. Application of grnn neural network in non-texture image
inpainting and restoration. Pattern Recognition Letters, 62:24–31, 2015.

[28] K. Yanai T. Nakamura, A. Zhu and S. Uchida. Scene text eraser,. in: 2017 14th IAPR
InternationalConference on Document Analysis and Recognition (ICDAR) IEEE,
1:832–837, 2017.

[29] J. Cheng P. Xiang, L. Wang, B. Zhang, and J. Wu. A deep network architecture
for image inpainting. in:2017 3rd IEEE International Conferenceon Computer and
Communications (ICCC), IEEE, pages 1851–1856, 2017.

[30] Xiuxia Cai, Song, and Bin. Semantic object removal with convolutional neural net-
work feature-based inpainting approach. Multimedia Systems, 24:597–609, 2018.

[31] Oleksii Sidorov and Jon Yngve Hardeberg. Deep hyperspectral prior: Denoising,
inpainting, super-resolution, 2019.

[32] J. Donahue D. Pathak, P. Krahenbuhl, T. Darrell, and A. A. Efros. Context encoders:
Feature learning by inpainting. n: Proceedings of the IEEE conference on computer
vision and pattern recognition, pages 2536–2544, 2016.



Y. Douich et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6192 32 of 32

[33] A. Gupta C. Doersch, S. Singh, J. Sivic, and A. A. Efros. What makes paris look like
paris? Communications of the ACM 58, 12:103–110, 2015.

[34] A. Khosla B. Zhou, A. Lapedriza, A. Oliva, A. Torralba, and Places. A 10 million
image database for scene recognition. IEEE transactions on pattern analysis and
machine intelligence 40, 6:1452–1464, 2017.

[35] S. Zhang H. Xue and D. Cai. Depth image inpainting: Improving low rank matrix
completion with low gradient regularization. IEEE Transactions on Image Processing
26, 9:4311–4320, 2017.

[36] C. Fowlkes D. Martin, D. Tala J. Malik, and et al. A database of human segmented
natural images and its application to evaluating segmentation algorithms and mea-
suring ecological statistics. Iccv Vancouver, pages 4311–4320, 2001.

[37] H. Su O. Russakovsky, J. Deng, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy,
A. Khosla, M. Bernstein, and et al. Imagenet large scale visual recognition challenge.
International journal of computer vision 115, 3:211–252, 2015.

[38] X. Wang Z. Liu, P. Luo and X. Tang. Large-scale celeb faces attributes (celeba)
dataset. Retrieved August15 2018, 2018.

[39] Bastien Chopard, Jean-Luc Falcone, and Jonas Latt. The lattice boltzmann
advection-diffusion model revisited. European Physical Journal: Special Topics,
171:245–249, 2009.

[40] J PERONA, P.and MALIK. Scale space and edge detection using anisotropic diffu-
sion. IEEE Trans. Pattern Anal. Mach.Intell, pages 629–639, 1990.

[41] Leonid Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based
noise removal algorithms. Physica D: Nonlinear Phenomena, 60:259–268, 11 1992.