# Isogeometric Analysis For Image Registration and Segmentation Using Optimal Transport Problem

Mustapha Bahari, Abderrahmane Habbal, Ahmed Ratnani

► To cite this version:

Mustapha Bahari, Abderrahmane Habbal, Ahmed Ratnani. Isogeometric Analysis For Image Registration and Segmentation Using Optimal Transport Problem. M3A 2024 - Mathematical Modeling with Modern Applications, Jun 2024, Istanbul, Turkey. pp.131-150, 10.1007/978-3-031-89041-3_7 hal-05399303

HAL Id: hal-05399303

https://inria.hal.science/hal-05399303v1

Submitted on 6 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Isogeometric Analysis For Image Registration and Segmentation Using Optimal Transport Problem

Mustapha BAHARI<fnref n="1" />,<fnref n="2" />, Abderrahmane HABBAL<fnref n="2" />, and Ahmed RATNANI<fnref n="1" />

**Abstract.** We present a fast and high order method for the problem of Image Registration, using Optimal Transport and the Isogeometric Analysis paradigm. Our method is based on the resolution of the Monge-Ampère equation and ensures the one-to-one property. In addition, the use of B-Splines allows to create a map that can be evaluated everywhere, and reduces the number of degrees of freedom needed to store the constructed (gradient) map, by using e.g. high order B-Splines functions.

**Keywords:** B-spline discretization, Optimal Transport, Isogeometric Analysis, Image registration, mass-preservation, mesh generation, r-refinement

## 1 Introduction

In this paper, we adapt our *optimal transport* methods developed for r-refinement and Isogeometric Analysis (see e.g. [2]) to the problem of Image Registration.

Optimal Transport [21, 19] has regain more popularity since the work of C. Villani [24]. It has been used in different areas such as economics, logistics, physics and recently in Machine Learning and Computer Vision to compare histograms [10]. Another interesting field of application is mesh generation and r-refinement [9, 13, 12].

Image registration is a crucial process in the fields of computer vision and medical imaging, aimed at aligning two or more images of the same scene or object acquired at different times, from different perspectives, or using different modalities [8, 16, 22, 26]. The goal is to establish a meaningful spatial correspondence between the pixels of these images, enabling comparisons, analysis, and fusion of information. In [17], the authors method for finding the registration mapping is based on a partial differential equation approach to the minimization of a Kantorovich-Wasserstein distance. The adaptive mesh technique for image registration is widely recognized, refer e.g. to [3, 14, 15]. Our approach relies on using optimal transport to generate adaptive mesh, optimally adapted to a given target image since it is obtained by optimally transporting an initial

<footnote n="1">¹ The UM6P Vanguard Center, Mohammed VI Polytechnic University, Benguerir 43150, Lot 660, Hay Moulay Rachid Morocco</footnote>

<footnote n="2">² Universite Cote d'Azur, Inria, CNRS, LJAD UMR 7351, Parc Valrose 06108 NICE CEDEX 2, France, abderrahmane.habbal@unice.fr,</footnote>
uniform -mesh- density to a density directly computed from the target image which accurately renders the features of interest in image registration.

The goal is to derive an adaptive mapping based on the intensity values $I$ of an MRI image. This mapping should generate an adaptive mesh using a density function derived from $I$, such that the resulting mesh closely resembles a reference image. Consequently, this mapping facilitates comparisons across various (patient) images. Specifically, in intensity-based image registration, we seek a mapping that aligns the density-based intensity $\rho_r$ of the reference image $I^r$ with the density-based intensity $\rho$ of the target image $I^t$, optimizing the similarity between the two images. Therefore, we compute a mapping $\Psi$ that aligns the two densities by solving the optimal transport problem described in the next section.

The overview of the paper is as follows, in Section 2, we give a brief recall on the Monge-Kantorovich problem and the mixed variational formulation method for the Monge-Ampère Equation (MAE). Section 3 is devoted to Isogeometric Analysis (IGA) approach and B-Splines Finite Elements. In Section 4 we describe the moving mesh algorithm for the Allen-Cahn equation. Section 5 is devoted to numerical experiments including image smoothing and segmentation and finally, in Section 6, we draw some conclusions.

## 2 The $L^2$ Monge-Kantorovich problem

We first start by recalling the $L^2$ Monge-Kantorovich problem and the relation with the Monge-Ampère equation. The $L^2$ Monge-Kantorovich problem is stated as follows:

**Definition 1 (L² Monge-Kantorovich problem).** Let $\rho_r$ and $\rho$ be two given densities of equal masses, defined in $\mathbb{R}^d$.

Find a mapping $\mathbf{x}' = \Psi(\mathbf{x})$, $\mathbf{x}, \mathbf{x}' \in \mathbb{R}^d$, that transfers the density $\rho_r$ to $\rho$ and minimizes the transport cost

$$ \mathcal{J}[\Psi] = \int_{\mathbb{R}^d} |\Psi(\mathbf{x}) - \mathbf{x}|^2 \rho_0(\mathbf{x}) \, d\mathbf{x}. \qquad (2.1) $$

The density transfer means that, for all Borelian set $\Omega \subset \mathbb{R}^d$,

$$ \int_{\Psi^{-1}(\Omega)} \rho_r(\mathbf{x}) \, d\mathbf{x} = \int_{\Omega} \rho(\mathbf{x}') \, d\mathbf{x}'. \qquad (2.2) $$

If the mapping $\Psi$ is a smooth one-to-one map, then equation (2.2) can be written as

$$ \rho(\Psi(\mathbf{x})) \det \left( \frac{\partial \Psi}{\partial \mathbf{x}} \right) = \rho_r(\mathbf{x}). \qquad (2.3) $$

It has been proved in [7] that the solution of the $L^2$ Monge-Kantorovich problem exists. Moreover, it is stated that the mapping $\Psi$ is the gradient of a convex
potential function $\phi$. Using equation (2.3), for $\rho_r$ and $\rho$ defined in the computational domain $\Omega_c$ ($\Omega_c = [0, 1]^d$) we get the Monge-Ampère equation (MA)

$$ \rho(\nabla\phi(\mathbf{x})) \det H(\phi(\mathbf{x})) = \rho_r(\mathbf{x}), \quad \mathbf{x} \in \Omega_c, \qquad (2.4) $$

where $H(\phi)$ is the Hessian matrix of the potential $\phi$. In order to solve the Monge-Ampère equation, we need additional boundary conditions. In our case, we must ensure that $\partial\Omega_c$ is mapped onto $\partial\Omega_c$, i.e.

$$ \nabla\phi(\partial\Omega_c) = \partial\Omega_c, $$

which writes:

$$ \nabla\phi \cdot \mathbf{n} = \mathbf{x} \cdot \mathbf{n}, \quad \text{on } \Omega_c. \qquad (2.5) $$

The Monge-Ampère equation (Eq. 2.4) augmented with the boundary conditions (Eq. 2.5), writes

$$ \begin{cases} \rho(\nabla\phi) \det H(\phi) = \rho_r(\mathbf{x}) & \text{in } \Omega_c, \\ \nabla\phi \cdot \mathbf{n} = \mathbf{x} \cdot \mathbf{n} & \text{on } \partial\Omega_c. \end{cases} $$

Hence, our aim is to solve the following nonlinear problem

$$ \begin{cases} \det(H(\phi)) = f(\nabla\phi) & \text{in } \Omega_c, \\ \nabla\phi \cdot \mathbf{n} = \mathbf{x} \cdot \mathbf{n} & \text{on } \partial\Omega_c, \end{cases} \qquad (2.6) $$

where $f(\nabla\phi) = \frac{\rho_r}{\rho(\nabla\phi)}$. Later on in our numerical experiments, we shall take $\rho_r = 1$ which corresponds to uniform mesh. In this case, we do not align two images; instead, we match the given target image.

### Monge-Ampère solver

Our fast solver is based on the Benamou-Froese-Oberman (BFO) method [4]. In their work, the authors proved that when $u \in H^2(\Omega)$ is a solution of the Monge-Ampère equation (Eq. 2.6), then it is a fixed point of the following operator $T$ (Eq. 2.7):

$$ \begin{gathered} T: H^2(\Omega_c) \to H^2(\Omega_c) \\ u \to T(u) = (\Delta)^{-1} \mathcal{G}(\nabla u), \end{gathered} \qquad (2.7) $$

where $\mathcal{G}$ is defined as

$$ \mathcal{G}(\nabla u) = \left( (\nabla \cdot \nabla u)^d + d! (f(\nabla u) - \det(\nabla(\nabla u))) \right)^{\frac{1}{d}}. \qquad (2.8) $$

In practice, since we are considering Neumann boundary conditions on the whole boundary, any approximate system will introduce boundary errors [2]. For this
work, we shall use a mixed variational formulation method. The goal is to integrate Neumann boundary conditions in strong form, hence (2.6) and (2.7) lead to the following problem

$$ 
\left\{ 
\begin{array}{ll} 
\nabla \cdot \mathbf{\Psi} = \mathcal{G}(\mathbf{\Psi}) & \text{in } \Omega_c, \\ 
\mathbf{\Psi} = \nabla \phi & \text{in } \Omega_c, \\ 
\mathbf{\Psi} \cdot \mathbf{n} = \mathbf{x} \cdot \mathbf{n} & \text{on } \partial\Omega_c, \\ 
\phi & \text{convex.} 
\end{array} 
\right. 
\qquad (2.9) 
$$

Now, let us introduce the decomposition $\mathbf{\Psi} = \mathbf{\Psi}_0 + \mathbf{\Psi}_N$ such that $\mathbf{\Psi}_0 \cdot \mathbf{n} = 0$ and $\mathbf{\Psi}_N \cdot \mathbf{n} = \mathbf{x} \cdot \mathbf{n}$. Consequently, we solve the problem (2.9) using the Picard iteration algorithm.

**Algorithm 1: Picard algorithm**

1: Given an initial value $\mathbf{\Psi}^0$.

2: Repeat until convergence: Compute $\mathbf{\Psi}^{n+1} = \mathbf{\Psi}_0^{n+1} + \mathbf{\Psi}_N$ such that $\mathbf{\Psi}_N$ given by $\mathbf{\Psi}_N \cdot \mathbf{n} = \mathbf{x} \cdot \mathbf{n}$, and $\mathbf{\Psi}_0^{n+1}$ solution of

$$ 
\left\{ 
\begin{array}{ll} 
\nabla \cdot \mathbf{\Psi}_0^{n+1} = \mathcal{G}(\mathbf{\Psi}^n) - \nabla \cdot \mathbf{\Psi}_N & \text{in } \Omega_c, \\ 
\mathbf{\Psi}_0^{n+1} = \nabla \phi - \mathbf{\Psi}_N & \text{in } \Omega_c, \\ 
\mathbf{\Psi}_0^{n+1} \cdot \mathbf{n} = 0 & \text{on } \partial\Omega_c. 
\end{array} 
\right. 
\qquad (2.10) 
$$

Finally, the density is developed and adjusted to assimilate the data $I$.

## 2.1 Density function

Given the target image (respectively, the reference image), we first extract the intensity data represented by $I$. Then, we define our density, denoted $\rho$, as follows:

$$ \rho = \frac{\alpha_i \Pi_p(I) + \max(I) - (\alpha_i - 1) \min(I)}{\max(I) - \min(I)}, \qquad (2.11) $$

where $\Pi_p$ represents the $L^2$-projection (or quasi-interpolant operator [23]), ensuring that $\rho \in \mathcal{V}_h$. To ensure stability and convergence, the intensity of a density function $\alpha_i$ has been chosen to be less than 100.

Therefore, employing a multilevel method, the formulation for the adaptive mesh in image registration can be expressed as follows:

## 3 B-splines based Isogeometric Analysis method

Isogeometric analysis (IGA) [18, 11] is an approach to solve partial differential equations using smooth Splines as basis functions. In IGA, the exact geometry representation is introduced within the analysis, as in opposition to the
**Algorithm 2:** Adaptive Mesh for Image Registration

**Require:** Target image, Reference image, initial fine mesh $Q_h$, coarse mesh $Q_H$

**Ensure:** Adaptive mapping $\Psi_I$ corresponding to adaptive mesh for the image.

1. Extract the intensity data $I$ of the image.

2. Compute $\rho$ and $\rho_r$ in $Q_h$ from (2.11).

3. Compute initial guess $\Psi_I$ in $Q_H$ solution of (2.9).

4. While $H > h$, do

1. Update $H \leftarrow H/2$.

2. Prolongate $\Psi_I$ into the new grid $Q_H$.

3. Use $\Psi_I$ as an initial guess in the Picard iteration and compute the updated solution in $Q_h$ from Algorithm 1.

isoparametric approach. The novelty of the IGA paradigm is to develop a close connection between Computer Aided Design (CAD) and Finite Elements Analysis (FEA), with the aim of using CAD models directly and avoiding the time-consuming meshing processes. Since its introduction, IGA has been used in many engineering application, computational physics as well as many other fields.

### A brief recall on B-splines

Given two positive integers $p$ and $n$, we consider the (ordered) open knot vectors

$$ \Xi = \{\xi_1, \xi_2, \dots, \xi_{n+p+1}\}, $$

such that

$$ \xi_1 = \dots = \xi_{p+1} < \xi_{p+2} < \dots < \xi_n < \xi_{n+1} = \dots = \xi_{n+p+1}, $$

where $p$ is the degree of the B-spline and $n$ is the dimension of the associated Schoenberg space. The univariate B-spline basis functions $B_{i,p}$, $i = 1, \dots, n$ are defined recursively by the well known Cox-de Boor recursion formula:

$$ N_i^p(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} N_i^{p-1}(\xi) + \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} B_{i+1}^{p-1}(\xi), \quad (3.1) $$

where

$$ N_{i,0}(\xi) = \begin{cases} 1 & \text{if } \xi \in [\xi_i, \xi_{i+1}), \\ 0 & \text{otherwise.} \end{cases} $$

In (3.1), we adopt the convention $0/0 = 0$. In the multi-dimensional case, B-splines can be constructed using tensor products, $N_{j_1, j_2, \dots, j_d}^{p_1, p_2, \dots, p_d} = N_{j_1}^{p_1} \otimes \dots \otimes N_{j_1}^{p_1}$. B-splines present many interesting properties, among them we find:

1. B-splines are piecewise polynomial of degree $p$.

2. B-Spline has minimal compact support: $\text{supp}(N_i^p) = [\xi_i, \xi_{i+p+1}]$.

3. If $\xi \in [\xi_i, \xi_{i+1})$, then only the B-splines $\{N_{i-p}^p, \dots, N_i^p\}$ are non vanishing at $\xi$.

4. Non-negativity: $N_i^p \ge 0$ in $[\xi_i, \xi_{i+p+1})$.
5. Partition of unity : $\sum_i N_i^p(x) = 1$ for all $x \in (0, 1)$.

6. Local linear independence : B-splines are linearly independent on each interval $[\xi_i, \xi_{i+1})$.

In figure (Fig. 1), we show an example of a family of quadratic B-Splines associated to the knot sequence $\Xi = \{0, 0, 0, 1, 2, 2, 3, 4, 5, 5, 5\}$. As we can see, repeating the internal knots leads to reducing the regularity of the associated B-Splines. For a given knot sequence $\Xi$, the Schoenberg space will be denoted

<figure><img src="image_2.png" /><figcaption>Fig. 1: Quadratic B-Splines generated using the knots sequence $\Xi = \{0, 0, 0, 1, 2, 2, 3, 4, 5, 5, 5\}$.</figcaption></figure>

by $S(\Xi, p)$, but in order to simplify the notations, we shall use $S_p$. In the multi-dimensional case, given $(\Xi_s)_{1\le s\le d}$ and $\mathbf{p} := (p_s)_{1\le s\le d}$, the generated Schoenberg space is defined as the generalized cartesian product of the 1D spaces, i.e. $\mathbf{X}_{1\le s\le d} S_p$.

Using B-Splines, one can create parametric curves by associating a coefficient to each basis function. These coefficients are called Control Points. Let $(\mathbf{P}_i)_{1\le i\le n} \in \mathbb{R}^d$ be a sequence of control points. The B-spline curve in $\mathbb{R}^d$ associated to $\Xi = (\xi_i)_{1\le i\le n+p+1}$ and $(\mathbf{P}_i)_{1\le i\le n}$ is defined by :

$$C(\xi) = \sum_{i=0}^{n} N_{i}^{p}(\xi) \mathbf{P}_{i}.$$

In figure (Fig. 2), we plot two examples of quadratic B-Splines, showing the impact of changing the regularity of the Schoenberg space. Notice that the use of open knot vector leads to an interpolating curve at the endpoints. Moreover, one can construct any conic using Rational B-Splines (see Fig. 3 for the case of a circle). For more details, we refer the reader to [20, 11].

*Variational Formulation for Monge-Ampère using IGA*
The variational formulation of the problem (Eq. 2.10) writes
<figure><img src="image_3.png" /><figcaption>Fig. 2: Quadratic B-Spline curve using the knot vectors (left) $\Xi = \{0, 0, 0, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, 1, 1, 1\}$ (right) $\Xi = \{0, 0, 0, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, 1, 1, 1\}$ with a cusp.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3: Circle as four concatenated NURBS arcs using the knot vector $\Xi = \{0, 0, 0, \frac{1}{4}, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, \frac{3}{4}, 1, 1, 1\}$</figcaption></figure>

Find $(\Psi_0^{n+1}, \phi^{n+1}) \in H_0(\text{div}, \Omega_c) \times L_0^2(\Omega_c)$ such that $\forall (v, \varphi) \in H_0(\text{div}, \Omega_c) \times L_0^2(\Omega_c)$

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \Psi_0^{n+1} \cdot v + \int_{\Omega} \phi^{n+1} \nabla \cdot v = - \int_{\Omega} \Psi_N \cdot v, \\ \displaystyle \int_{\Omega} (\nabla \cdot \Psi_0^{n+1}) \varphi = \int_{\Omega} \mathcal{G}(\Psi^n) \varphi - \int_{\Omega} (\nabla \cdot \Psi_N) \varphi, \end{array} \right. \quad (3.2) $$

using the Galerkin method, we shall consider the discrete version of the problem (3.2), on a subspace $V_h(\text{div}, \Omega_c) \times V_h(L^2, \Omega_c)$ of $H_0(\text{div}, \Omega_c) \times L_0^2(\Omega_c)$. In our case, the finite elements space must be $C^1$, which can be easily achieved using B-splines, especially if one uses the maximum regularity, i.e. $p-1$. In the sequel, we use the same definition the Algorithm 3 [2].

## 4 Allen-Cahn equation using IGA

In this section, we introduce a technique that combines the Allen-Cahn solver with mesh adaptation for smoothing shapes (Images) and image segmentation. The proposed approach (smoothing shapes) has been evaluated within the framework of a finite difference method and has demonstrated commendable performance across numerous examples [25]. This equation is widely used in literature, including image segmentation [5]. Here is the introduction to the Allen-Cahn
equation:

$$ \frac{\partial \Phi}{\partial t} = - \frac{\mathcal{E}'(\Phi)}{\varepsilon^2} + \Delta \Phi + \lambda (f - \Phi) \quad \text{in } \Omega_c, \qquad (4.1) $$

where $\mathcal{E}(\Phi) = 0.25 (\Phi^2 - 1)^2$ is the potential energy function that takes values between -1 and 1, $\varepsilon$ represents the thickness of the transition layer, and $\lambda$ is the fidelity strength coefficient. Here, $f$ denotes the bounded value of a given image $I$, defined by

$$ f = \frac{2I - \max(I) - \min(I)}{\max(I) - \min(I)}. \qquad (4.2) $$

The equation (4.1) is formulated with homogeneous Neumann boundary conditions.

*Allen-Cahn numerical scheme* We begin by introducing an additional quasi-Lagrangian term, making the equation (4.1) equivalent to:

$$ \frac{\partial (\Phi \circ \Psi)}{\partial t} \circ \Psi^{-1} = - \frac{\mathcal{E}'(\Phi)}{\varepsilon^2} + \Delta \Phi + \lambda (f - \Phi) + \left( \nabla \Psi \frac{\partial \Psi}{\partial t} \right) \circ \Psi^{-1} \cdot \nabla \Phi. \quad (4.3) $$

Therefore, we introduce the implicit scheme for the Allen-Cahn equation: given the initial solution $\Phi^n$, we aim to find $\Phi^{n+1} \in H^1(\Omega_c)$ such that

$$ \frac{\Phi^{n+1} - \Phi^n \circ \Psi^n \circ (\Psi^{n+1})^{-1}}{\Delta t} = - \frac{\mathcal{E}'(\Phi^{n+1})}{\varepsilon^2} + \Delta\Phi^{n+1} + \lambda(f - \Phi^{n+1}) + \left( \nabla\Psi^{n+1} \frac{\Psi^{n+1} - \Psi^n}{\Delta t} \right) \circ (\Psi^{n+1})^{-1} \cdot \nabla\Phi^{n+1}. \quad (4.4) $$

The variational formulation is obtained by multiplying equation (4.4) by test function and integrating it over $\Omega_c$. Consequently, upon applying Green's formula, we obtain: Find $\Phi^{n+1} \in H_0^1(\Omega_c)$ such that $\forall v \in H_0^1(\Omega_c)$

$$ \int_{\Omega_c} \Phi^{n+1} v \, dx - \int_{\Omega_c} \Phi^n \circ \Psi^n \circ (\Psi^{n+1})^{-1} v \, dx + \Delta t \int_{\Omega_c} \nabla \Phi^{n+1} \cdot \nabla v \, dx \quad (4.5) $$

$$ + \Delta t \int_{\Omega_c} \frac{\mathcal{E}'(\Phi^{n+1})}{\varepsilon^2} v \, dx + \Delta t \int_{\Omega_c} \lambda(f - \Phi) v \, dx $$

$$ + \int_{\Omega_c} \left( \nabla \Psi^{n+1} \frac{\Psi^{n+1} - \Psi^n}{\Delta t} \right) \circ (\Psi^{n+1})^{-1} \cdot \nabla \Phi^{n+1} v \, dx = 0. \quad (4.6) $$

We recall that the Galerkin approximation space on the adapted mesh is given by $S_{p_1} \times S_{p_2} \circ (\Psi^n)^{-1}$ at each time step $t_n = n\Delta t$. For more details on the Galerkin variational formulation and the numerical scheme based on a quasi-Lagrangian term, see [1]. In this context, the Picard iterations are limited to a single iteration, as they function as time iterations for the Parabolic Monge-Ampère equation.
# 5 Numerical experiments

In this section, we illustrate several tests of our Algorithm 2 using various images extracted from the literature. We fix a tolerance in Picard iteration Algorithm 1 of $L^2$-residual error at $10^{-08}$.

For our numerical simulations, we used a laptop with a quad-core 2.6 Ghz Intel Core i7, running on Ubuntu OS. The code was implemented in pure Python, then accelerated with Pyccel [6], a Python static compiler for Scientific High Performance Computing, that translates Python code to C/Fortran code.

We demonstrate the effectiveness of mesh adaptation by evaluating the following:

$$Qual^2 = \int_{\Omega_c} \frac{||\Psi - x||^2}{\det(\nabla \Psi)} dx \quad (5.1)$$

which quantifies the global displacement of the grid points.

## 5.1 Application to Image Registration

In this section, we illustrate several tests of our Algorithm 2 using various images extracted from the literature in the following examples.

**Example 1 (brain MRI image).** In this example, we provide an example of image registration that aligns two MRI images. In Table 1 we show the number of iterations until convergence, as well as the CPU time taken. As we can see, the number of iterations is stable, and moving from cubic to quartic splines doubles the CPU-time. In Figure 4, we plot the result of deforming the reference image, the generated grid (mesh), and the density function values over the entire computation domain.

<table><caption>Table 1: **Example 1:** CPU-time needed to solve the Monge-Ampère equation and global mesh displacement using bicubic B-Spline.</caption><thead><tr><th rowspan="2"># cells</th><th colspan="3">p = 3</th><th colspan="3">p = 4</th></tr><tr><th>CPU-time (s)</th><th># iter</th><th>Qual</th><th>CPU-time (s)</th><th>#iter</th><th>Qual</th></tr></thead><tbody><tr><td>16</td><td>0.138</td><td>30</td><td>2.656e-02</td><td>0.253</td><td>33</td><td>2.721e-02</td></tr><tr><td>32</td><td>0.625</td><td>30</td><td>2.586e-02</td><td>1.065</td><td>30</td><td>2.554e-02</td></tr><tr><td>64</td><td>2.717</td><td>24</td><td>2.892e-02</td><td>3.54</td><td>19</td><td>2.617e-02</td></tr><tr><td>128</td><td>8.903</td><td>19</td><td>2.621e-02</td><td>14.173</td><td>22</td><td>2.622e-02</td></tr></tbody></table>

**Example 2 (Flower).** In this example, we test our solver with image of flower. In Table 2 we show the number of iterations until convergence, as well as the CPU time taken. In Figure 5, we plot the original image, the generated grid (mesh) as well as the density function values over the whole computation domain. We notice the same remarks as for the previous example. From the image, we can observe that the mesh follows the luminosity (intensity) of the image, with coarse grid in the black regions and fine grids in the white regions.
<figure><img src="image_5.png" /><figcaption>Fig. 4: **Example 1:** (top-left) Brain MRI target image (top-middle) Adaptive meshes and (top-right) Density function $\rho$. (bottom-left) Deformed reference image (bottom-middle) Brain MRI reference image and (bottom-right) Density function $\rho_r$.</figcaption></figure>

<table><caption>Table 2: **Example 2:** Number of Picard's iterations (#iter) and CPU-time.</caption><thead><tr><th rowspan="2"># cells</th><th colspan="2">p = 3</th><th colspan="2">p = 4</th></tr><tr><th>CPU-time (s)</th><th>#iter</th><th>CPU-time (s)</th><th>#iter</th></tr></thead><tbody><tr><td>16</td><td>0.565</td><td>30 3.997e-02</td><td>0.737</td><td>28 3.985e-02</td></tr><tr><td>32</td><td>1.224</td><td>19 3.983e-02</td><td>1.612</td><td>19 3.974e-02</td></tr><tr><td>64</td><td>3.403</td><td>17 3.944e-02</td><td>3.864</td><td>16 3.944e-02</td></tr><tr><td>128</td><td>8.741</td><td>13 3.931e-02</td><td>10.394</td><td>13 3.930e-02</td></tr></tbody></table>

*Example 3 (Grayscale).* In Table 3 we show the number of iterations until convergence, as well as the CPU time taken. In Figure 6, we plot the original image, the generated grid (mesh) as well as the density function values over the whole computation domain. We notice the same remarks as for the previous example. The $C^1$ regularity construction in our density ensures that the changes in adaptivity are smooth.

*Example 4 (Cat).* In this example, we test our solver with an image of a cat. In Table 4 we show the number of iterations until convergence, as well as the CPU time taken. In Figure 7, we plot the original image, the generated grid (mesh) as well as the density function values over the whole computation domain. We notice the same remarks as for the previous example.
<figure><img src="image_6.png" /><figcaption>Fig. 5: **Example 2**: (left) flower (middle) Adaptive meshes and (right) Density function.</figcaption></figure>

<table><caption>Table 3: **Example 3**: Number of Picard's iterations (#iter) and CPU-time.</caption><thead><tr><th rowspan="2"># cells</th><th colspan="3">p = 3</th><th colspan="3">p = 4</th></tr><tr><th>CPU-time (s)</th><th>#iter</th><th>Qual</th><th>CPU-time (s)</th><th>#iter</th><th>Qual</th></tr></thead><tbody><tr><td>16</td><td>0.428</td><td>25 3.545e-02</td><td></td><td>0.727</td><td>30 3.527e-02</td><td></td></tr><tr><td>32</td><td>1.171</td><td>21 3.532e-02</td><td></td><td>1.764</td><td>19 3.530e-02</td><td></td></tr><tr><td>64</td><td>2.789</td><td>14 3.522e-02</td><td></td><td>4.857</td><td>18 3.520e-02</td><td></td></tr><tr><td>128</td><td>7.92</td><td>13 3.519e-02</td><td></td><td>11.247</td><td>11 3.519e-02</td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6: **Example 3**: (left) Grayscale (middle) Adaptive meshes and (right) Density function.</figcaption></figure>

<table><caption>Table 4: **Example 4**: Number of Picard's iterations (#iter) and CPU-time.</caption><thead><tr><th rowspan="2"># cells</th><th colspan="3">p = 3</th><th colspan="3">p = 4</th></tr><tr><th>CPU-time (s)</th><th>#iter</th><th>Qual</th><th>CPU-time (s)</th><th>#iter</th><th>Qual</th></tr></thead><tbody><tr><td>16</td><td>0.458</td><td>34 8.663e-02</td><td></td><td>0.81</td><td>35 8.643e-02</td><td></td></tr><tr><td>32</td><td>1.181</td><td>28 8.902e-02</td><td></td><td>1.954</td><td>26 8.858e-02</td><td></td></tr><tr><td>64</td><td>3.384</td><td>24 8.934e-02</td><td></td><td>5.39</td><td>24 8.934e-02</td><td></td></tr><tr><td>128</td><td>8.899</td><td>17 8.930e-02</td><td></td><td>14.057</td><td>17 8.927e-02</td><td></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Fig. 7: **Example 4**: (left) Cat (middle) Adaptive meshes and (right) Density function.</figcaption></figure>

*Example 5 (Print)*. In this example, we test our solver with an image of Print. In Table 5 we show the number of iterations until convergence, as well as the CPU time taken. In Figure 8, we plot the original image, the generated grid (mesh) as well as the density function values over the whole computation domain. The global displacement of grid points, denoted as "Qual", increases notably when the target image is in black and white or an MRI image compared to others. This observation suggests that we can enhance the resulting mesh by providing the intensity data *I* with more details from the image.

<table><caption>Table 5: **Example 5**: Number of Picard's iterations (#iter) and CPU-time.</caption><thead><tr><th rowspan="2"># cells</th><th colspan="3">p = 3</th><th colspan="3">p = 4</th></tr><tr><th>CPU-time (s)</th><th>#iter</th><th>Qual</th><th>CPU-time (s)</th><th>#iter</th><th>Qual</th></tr></thead><tbody><tr><td>16</td><td>0.577</td><td>47</td><td>9.336e-02</td><td>1.214</td><td>45</td><td>9.371e-02</td></tr><tr><td>32</td><td>1.294</td><td>27</td><td>9.353e-02</td><td>2.43</td><td>24</td><td>9.358e-02</td></tr><tr><td>64</td><td>3.895</td><td>30</td><td>9.788e-02</td><td>7.061</td><td>30</td><td>9.781e-02</td></tr><tr><td>128</td><td>13.289</td><td>27</td><td>9.862e-02</td><td>22.869</td><td>28</td><td>9.865e-02</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8: **Example 5**: (left) Print (middle) Adaptive meshes and (right) Density function.</figcaption></figure>
*Example 6 (Hand MRI)*. In this example, we test our solver by using a hand MRI test case. In Table 6 we show the number of iterations until convergence, as well as the CPU time taken. In Figure 9, we plot the original image, the generated grid (mesh) and the density function values over the whole computation domain. We notice the same remarks as in the previous example.

<table><caption>Table 6: **Example 6**: Number of Picard's iterations (#iter) and CPU-time.</caption><thead><tr><th rowspan="2"># cells</th><th colspan="3">p = 3</th><th colspan="3">p = 4</th></tr><tr><th>CPU-time (s)</th><th>#iter</th><th>Qual</th><th>CPU-time (s)</th><th>#iter</th><th>Qual</th></tr></thead><tbody><tr><td>16</td><td>0.46</td><td>38 1.340e-01</td><td></td><td>0.781</td><td>38 1.332e-01</td><td></td></tr><tr><td>32</td><td>1.138</td><td>25 1.319e-01</td><td></td><td>1.994</td><td>26 1.318e-01</td><td></td></tr><tr><td>64</td><td>2.87</td><td>20 1.317e-01</td><td></td><td>4.413</td><td>17 1.316e-01</td><td></td></tr><tr><td>128</td><td>7.842</td><td>15 1.316e-01</td><td></td><td>11.95</td><td>15 1.315e-01</td><td></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Fig. 9: **Example 6**: (left) hand (middle) Adaptive meshes and (right) Density function.</figcaption></figure>

*Example 7*. In this example, we present additional illustrations in Figure 10. The images are chosen carefully so that the data *I* extracted from the image helps our density function to align significantly with the images.
<figure><img src="image_11.png" /><figcaption>Fig. 10: **Example 5.1:** (left) Original image (middle) Adaptive meshes and (right) Density function.</figcaption></figure>

*Example 8.* In this last example, we use our solver to transport the hand MRI image on a puzzle piece. The method replicates the approach employed in Example 6. In Figure 11, we plot the initial grid in the unit square, the generated grid (mesh) as well as the image of the hand onto a Puzzle piece.

<figure><img src="image_12.png" /><figcaption>Fig. 11: **Example 8:** (left) Initial adapted mesh (middle) its image under the initial mapping and (right) Density function.</figcaption></figure>
## 5.2 Numerical results of segmentation technique applied within the AC equation

For the analysis of numerical results, we define the energy functional as:

$$ \text{Energy} = \int_{\Omega_c} \frac{\mathcal{E}(\Phi)}{\varepsilon^2} + \frac{|\nabla \Phi|^2}{2} + \frac{\lambda}{2} (f - \Phi)^2 dx. \quad (5.2) $$

Additionally, we quantify the L²-error using:

$$ (\mathbf{M}_2)^2 = \int_{\Omega_c} (f - \Phi)^2 dx. \quad (5.3) $$

The following results are given using quadratic B-splines, utilizing 64 elements in each direction, $\lambda = 1000$, $\varepsilon = 2h / (\sqrt{2} \tanh^{-1}(0.9))$, and $\Delta t = \frac{h^2}{7}$.

*Example 9.* In this example, we validate Allen-Cahn solver by examining the fidelity term with two distinct initial and data configurations. We consider the following initial data

$$ f(x, y) = \begin{cases} 1, & \text{if } \sqrt{(x - 0.4)^2 + (y - 0.45)^2} < 0.2, \\ -1, & \text{otherwise} \end{cases} \quad (5.4) $$

and

$$ \Phi(x, y, 0) = \begin{cases} 1, & \text{if } |x - 0.75| < 0.2 \text{ and } |y - 0.4| < 0.2, \\ -1, & \text{otherwise.} \end{cases} \quad (5.5) $$

Our numerical findings validate that with $\lambda = 1000$, the fidelity term significantly influences the solution to closely match the given data, see Figure 12. The energy and error evolution profiles remain consistent between the two solutions obtained using uniform and adaptive meshes Figure 13. The use of adaptive mesh results in an increased number of points at the interface, as illustrated in the right panel of Figure 13. The energy and L²-error consistently decrease over time while approaching the steady state rapidly using the adaptive mesh, as illustrated in Figure 13.

*Example 10 (Image segmentation).* In this example, we demonstrate the application of the Allen-Cahn solver for MRI brain tumor segmentation. We solve the Allen-Cahn equation in a uniform mesh, where the adaptive mesh dynamically follows the variation of the Allen-Cahn solution to detect the brain tumor by surrounding the affected area with adapted meshes, as shown in Figures 14-16. Another instance demonstrates the Allen-Cahn solution and the surrounding mesh encompassing the solid tumor located on the left side of the brain, as shown in Figures 14-16. The effectiveness of the adaptive mesh is demonstrated by the minimization of energy and L²-error in both examples, as shown in Figures 15-17.
<figure><img src="image_13.png" /><figcaption>Fig. 12: **Example 9:** Effect of $\lambda = 1000$. The evolution of approximate solution alongside their adaptive meshes at $t = 0$, $t = 2\Delta t$, $t = 10\Delta t$, and $t = 40\Delta t$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13: **Example 9:** Effect of $\lambda = 1000$. (left) The Energy profile (middle) $M_2$ profile, and (right) projection of approximate solution into $x = 0.5$ at $t = 290\Delta t$.</figcaption></figure>

# 6 Conclusions

We introduced an original approach which relies on using optimal transport to create computational meshes well adapted to solving Allen-Cahn equation in view of performing image registration. The adapted meshes are obtained by transporting an initial uniform density onto a density which accurately accounts for the features of interest in image segmentation. Our numerical methodologies include solving a mixed formulation of Monge-Ampère equation and solving Allen-Cahn nonlinear equation on adapted meshes. Several numerical experiments were provided that corroborate the efficiency of our approach, including challenging test-cases where geometries had singularities or features with different scales.
<figure><img src="image_15.png" /><figcaption>Fig. 14: **Example 10**: (top-left) Exact MRI Image (top-middle) adaptive mesh (top-right) final configuration (bottom) The results at $t = 1\Delta t, 500\Delta t, 1200\Delta t$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15: **Example 9**: (left) The Energy profile (middle) $M_2$ profile, and (right) projection of approximate solution into $x = 0.5$ at $t = 1200\Delta t$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 16: **Example 10**: (top-left) Exact MRI Image (top-middle) adaptive mesh (top-right) final configuration (bottom) The results at $t = 1\Delta t, 500\Delta t, 1200\Delta t$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17: **Example 9**: (left) The Energy profile (middle) $M_2$ profile, and (right) projection of approximate solution into $x = 0.5$ at $t = 1200\Delta t$.</figcaption></figure>
Bibliography

[1] Bahari, M., Habbal, A., Ratnani, A.: Isogeometric analysis based mesh adaptation for time dependent problems. Engineering with Computers pp. 1–24 (2024). DOI 10.1007/s00366-024-02009-8

[2] Bahari, M., Habbal, A., Ratnani, A., Sonnendrücker, E.: Adaptive Iso-geometric Analysis using optimal transport and their fast solvers. Computer Methods in Applied Mechanics and Engineering **418**, 116,570 (2024). DOI https://doi.org/10.1016/j.cma.2023.116570. URL https://www.sciencedirect.com/science/article/pii/S0045782523006941

[3] Barnafi, N., Gatica, G.N., Hurtado, D.E., Miranda, W., Ruiz-Baier, R.: Adaptive Mesh Refinement in Deformable Image Registration: A Posteriori Error Estimates for Primal and Mixed Formulations. SIAM Journal on Imaging Sciences **14**(3), 1238–1272 (2021). DOI 10.1137/20M1364333. URL https://doi.org/10.1137/20M1364333

[4] Benamou, J.D., Froese, B.D., Oberman, A.M.: Two numerical methods for the elliptic Monge-Ampere equation. ESAIM: Mathematical Modelling and Numerical Analysis **44**(4), 737–758 (2010)

[5] Beneš, M., Chalupecký, V., Mikula, K.: Geometrical image segmentation by the Allen–Cahn equation. Applied Numerical Mathematics **51**(2), 187–205 (2004). DOI https://doi.org/10.1016/j.apnum.2004.05.001. URL https://www.sciencedirect.com/science/article/pii/S0168927404000819

[6] Bourne, E., Güçlü, Y., Hadjout, S., Ratnani, A.: Pyccel: a Python-to-X transpiler for scientific high-performance computing. Journal of Open Source Software **8**(83), 4991 (2023). DOI 10.21105/joss.04991. URL https://doi.org/10.21105/joss.04991

[7] Brenier, Y.: Polar factorization and monotone rearrangement of vector-valued functions. Communications on pure and applied mathematics **44**(4), 375–417 (1991)

[8] Brown, L.G.: A survey of image registration techniques. ACM Comput. Surv. **24**(4), 325–376 (1992). DOI 10.1145/146370.146374. URL https://doi.org/10.1145/146370.146374

[9] Budd, C.J., Williams, J.F.: Parabolic Monge–Ampère methods for blow-up problems in several spatial dimensions. Journal of Physics A: Mathematical and General **39**(19), 5425 (2006)

[10] Chan, T., Esedoglu, S., Ni, K.: Histogram based segmentation using wasserstein distances. In: Scale Space and Variational Methods in Computer Vision: First International Conference, SSVM 2007, Ischia, Italy, May 30-June 2, 2007. Proceedings 1, pp. 697–708. Springer (2007)

[11] Cottrell, J., Hughes, T., Bazilevs, Y.: Isogeometric Analysis : toward Integration of CAD and FEA, first edn. John Wiley & Sons, Ltd (2009)

[12] Delzanno, G.L., Chacón, L., Finn, J.M., Chung, Y., Lapenta, G.: An optimal robust equidistribution method for two-dimensional grid adaptation based
on Monge-Kantorovich optimization. Journal of Computational Physics **227**(23), 9841-9864 (2008)

[13] Finn, J.M., Delzanno, G.L., Chacón, L.: Grid generation and adaptation by Monge-Kantorovich optimization in two and three dimensions. In: Proceedings of the 17th International Meshing Roundtable, pp. 551-568. Springer (2008)

[14] Fogtmann, M., Larsen, R.: Adaptive mesh generation for image registration and segmentation. In: 2013 IEEE International Conference on Image Processing, pp. 757-760 (2013). DOI 10.1109/ICIP.2013.6738156

[15] Haber, E., Heldmann, S., Modersitzki, J.: Adaptive Mesh Refinement for Nonparametric Image Registration. SIAM Journal on Scientific Computing **30**(6), 3012-3027 (2008). DOI 10.1137/070687724. URL https://doi.org/10.1137/070687724

[16] Haber, E., Modersitzki, J.: Numerical methods for volume preserving image registration. Inverse Problems **20**(5), 1621-1638 (2004). DOI 10.1088/0266-5611/20/5/018. URL https://dx.doi.org/10.1088/0266-5611/20/5/018

[17] Haker, S., Zhu, L., Tannenbaum, A., Angenent, S.: Optimal Mass Transport for Registration and Warping. Int. J. Comput. Vision **60**(3), 225-240 (2004)

[18] Hughes, T.J., Cottrell, J.A., Bazilevs, Y.: Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer methods in applied mechanics and engineering **194**(39-41), 4135-4195 (2005)

[19] Kantorovich, L.: On the translocation of masses, CR Dokl. Acad. Sci. URSS **37**, 191-201 (1942)

[20] L. Piegl, W.T.: The NURBS Book. Springer-Verlag, Berlin, Heidelberg (1995). Second ed.

[21] Monge, G.: Mémoire sur la théorie des déblais et des remblais. Histoire de l'Académie Royale des Sciences de Paris (1781)

[22] Sotiras, A., Davatzikos, C., Paragios, N.: Deformable medical image registration: A survey. IEEE Transactions on Medical Imaging **32**(7), 1153-1190 (2013). DOI 10.1109/TMI.2013.2265603

[23] Speleers, H., Manni, C.: Effortless quasi-interpolation in hierarchical spaces. Numerische Mathematik **132**, 155-184 (2016). DOI 10.1007/s00211-015-0711-z

[24] Villani, C.: Optimal transport: old and new, vol. 338. Springer (2009)

[25] Wang, J., Li, Y., Choi, Y., Lee, C., Kim, J.: Fast and Accurate Smoothing Method Using a Modified Allen-Cahn Equation. Computer-Aided Design **120**, 102,804 (2020). DOI https://doi.org/10.1016/j.cad.2019.102804. URL https://www.sciencedirect.com/science/article/pii/S0010448519305366

[26] Xie, Z., Farin, G.E.: Image registration using hierarchical b-splines. IEEE Transactions on Visualization and Computer Graphics **10**(1), 85-94 (2004). DOI 10.1109/TVCG.2004.1260760