# Numerical Solution of Diffusion Models  in Biomedical Imaging on Multicore Processors

Authors: Luisa D'Amore; Daniela Casaburi; Livia Marcellino; Almerico Murli
Journal: International Journal of Biomedical Imaging
 (2011)

---

## Abstract

In this paper, we consider nonlinear partial differential equations (PDEs) of diffusion/advection type underlying most problems in image analysis. As case study, we address the segmentation of medical structures. We perform a comparative study of numerical algorithms arising from using the semi-implicit and the fully implicit discretization schemes. Comparison criteria take into account both the accuracy and the efficiency of the algorithms. As measure of accuracy, we consider the Hausdorff distance and the residuals of numerical solvers, while as measure of efficiency we consider convergence history, execution time, speedup, and parallel efficiency. This analysis is carried out in a multicore-based parallel computing environment.


---

## Body

## 1. Introduction

High-quality images are crucial to accurately diagnose a patient or determine treatment. In addition to requiring the best images possible, safety is a crucial consideration. Many imaging systems use X-rays to provide a view of what is beneath a patient’s skin. X-ray radiation levels must be kept at a minimum to protect both patients and staff. As a result, raw image data can be extremely noisy. In order to provide clear images, algorithms designed to reduce noise are used to process the raw data and extract the image data while eliminating the noise. In video imaging applications, data often have to be processed at rates of 30 images per second or more. Filtering noisy input data and delivering clear, high-resolution images at these rates require tremendous computing power. This gave rise to the need of developing high-end computing algorithms for image processing and analysis which are able to exploit the high performance of advanced computing machines.In this paper, we focus on the computational kernels which arise as basic building blocks of the numerical solution of medical imaging applications described in terms of partial differential equations (PDEs) of parabolic/hyperbolic type. Such PDEs arise from the scale-space approach for description of most inverse problems in imaging [1]. One of the main reasons for using PDEs to describe image processing applications is that PDE models preserve the intrinsic locality of many image processing operations. Moreover, we can rely on standard and up-to-date literature and software about basic computational issues arising in such case (such as the construction of suitable discretization schemes, the availability of a range of algorithmic options, and the reuse of software libraries that allow the effective exploitation of high-performance computing resources). Finally, PDEs appear to be effectively implemented on advanced computing environments [2].We consider two standard discretization schemes of nonlinear time-dependent PDEs: semi-implicit scheme and fully implicit scheme [3]. The former leads to the solution of a linear system at each time (scale) step, while the computational kernel of the fully implicit scheme is the solution of a nonlinear system, to be performed at each time (scale) step. Taking into account that we aim to solve such problems on parallel computer in a scalable way, in the first case, we use, as linear solver, Krylov iterative methods (GMRES) with algebraic multigrid preconditioners (AMG) [4, 5]. Regarding the fully implicit scheme, we use the Jacobian-Free Newton-Krylov (JFNK) method as nonlinear solver [6].In recent years, multicore processors are becoming dominant systems in high-performance computing [7]. We provide a multicore implementation of numerical algorithms arising from using the semi-implicit and the implicit discretization schemes of nonlinear diffusion models underlying most problems in image analysis. Our implementation is based on parallel PETSc (Portable Extensible Toolkit for Scientific Computation) computing environment [8]. Parallel software uses a distributed memory model where the details of intercore communications and data managements are hidden within the PETSc parallel objects.The paper is organized as follows. In Section2, an overview of the PDE model equation used in describing some of inverse problems in imaging applications will be given. Then, the segmentation problem of medical structures is discussed. Numerical approach will be introduced in Section 3. Section 4 is devoted to the discussion of numerical algorithms based on semi-implicit and implicit numerical schemes. In Section 6, we describe the experiments that we carried out to show both the accuracy and the performance of these algorithms, while Section 7 concludes the work.


## 2. Diffusion Models Arising in Medical Imaging

The task in medical imaging is to provide in a noninvasive way information about the internal structure of the human body. The basic principle is that the patient is scanned by applying some sort of radiation and its interaction with the body is measured. This result is the data, whose origin has to be identified. Hence, we face an inverse problem. Most medical imaging problems lead to ill-posed (inverse) problems in the sense of Hadamard [9–11]. A standard approach for dealing with such intrinsic instability is to use additional information to construct families of approximate solution. This principle characterizes regularization methods that, starting from the milestone Tikhonov regularization [12], are now one of the most powerful tools for solution of inverse ill-posed problems. In 1992, Rudin et al. introduced the first nonquadratic regularization functional (i.e., the total variation regularization) [13] to denoise images. Moreover, the authors derive the Euler-Lagrange equations as a time-dependent PDE. In the same years Perona and Malik introduced the first nonlinear multiscale analysis [14].Scale-space theory has been developed by the computer vision community to handle the multiscale nature of image data. A main argument behind its construction is that if no prior information is available about what are the appropriate scales for a given data set, then the only reasonable approach for a vision system is to represent the input data at multiple scales. This means that the original imageu(x),x∈ℜ2 should be embedded into a one-parameter family of derived images, in which fine-scale structures are successively suppressed:
(1)SSτ:τ∈R⟶u(x,τ).
A crucial requirement is that structures at coarse scales in the multiscale representation should constitute simplifications of corresponding structures at finer scales—they should not be accidental phenomena created by the method for suppressing fine-scale structures. A main result is that if rather general conditions are imposed on the types of computations that are to be performed, then convolution by the Gaussian kernel and its derivatives is singled out as a canonical class of smoothing transformations [15, 16].A strong relation between regularization approaches and the scale-space approach exists via the Euler-Lagrange equation of regularization functionals: it consists of a PDE of parabolic/hyperbolic (diffusion/advection) type [17], defined as follows.Nonlinear Diffusion Models
Letx=(x,y)∈ℜ2 and u(x,τ), defined in [0,T]×Ω, be the scale-space representation of the brightness function image u(x) defined in Ω⊂ℜ2 describing the real (and unknown) object and u0(x) the observed image (the input data). Let us consider the following PDE problem:
(2)∂u(x,τ)∂τ=|∇u|∇⋅(g(u)∇u|∇u|)τ∈[0,T],(x,y)∈Ωu(0,x)=u0(x)τ=0,(x,y)∈Ω.[0,T] is the scale (time) interval; g(v) is a nonincreasing real valued function (for v>0) which tends to zero as v→∞. Initial and boundary conditions will be provided according to the problem to be solved (denoising, segmentation, deblurring, registration, and so on).Equations in (2) describe the motion of a curve (a moving front) with a speed depending on a local curvature. Such equations, known as level set equations, were first introduced in [18]. The original idea behind the level set method was a simple one. Given an interface Γ in ℜn of codimension one (i.e., its dimension is n-1), bounding an (perhaps multiply connected) open region Ω, we wish to analyze and compute its subsequent motion under a velocity field v⃗. This velocity can depend on position, time, the geometry of the interface (e.g., its normal or its mean curvature), and the external physics. The idea, as devised in 1988 by Osher and Sethian, is merely to define a smooth (at least Lipschitz continuous) function ϕ(x,t), that represents the interface Γ as the set where ϕ(x,t)=0. Thus, the interface is to be captured for all later time, by merely locating the set Γ(t) for which ϕ vanishes. The motion is analyzed by convecting the ϕ values (levels) with the velocity field v. This elementary equation is
(3)∂ϕ∂t+v⋅∇ϕ=0.
Actually, only the normal component of v is needed:
(4)vN=v⋅∇ϕ|∇ϕ|,
and the motion equation becomes
(5)∂ϕ∂t+vN⋅|∇ϕ|=0.
Taking into account that the mean curvature of Γ(t) is
(6)cur=-∇⋅(∇ϕ|∇ϕ|),
equation (5) describes the motion of Γ(t) under a speed vN proportional to its curvature cur (Mean Curvature Motion, MCM equation) [18–20]. This basic model has received a lot of attention because of its geometrical interpretation. Indeed, the level sets of the image solution or level surfaces in 3D images move in the normal direction with a speed proportional to their mean curvature. In image processing, equations like (5) arise in nonlinear filtration, edge detection, image enhancement, and so forth, when we are dealing with geometrical features of the image-like silhouette of object corresponding to level line of image intensity function. Finally, the level set approach instead of explicitly following the moving interface itself takes the original interface Γ and embeds it in higher dimensional scalar function u, defined over the entire image domain. The interface Γ is now represented implicitly as the zeroth level set (or contour) of this function, which varies with space and time (scale) using the partial differential equation in (2), containing terms that are either hyperbolic or parabolic. The theoretical study of the PDE was done by [21] which proved existence and uniqueness of viscosity solutions.


### 2.1. A Case Study: Image Segmentation

In this paper, we use equations (2) for image segmentation. The task of image segmentation is to find a collection of nonoverlapping subregions of a given image. In medical imaging, for example, one might want to segment the tumor or the white matter of a brain from a given MRI image.The idea behind level set (also known implicit active contours, or implicit deformable models) for image segmentation is quite simple. The user specifies an initial guess for the contour, which is then moved by image-driven forces to the boundaries of the desired objects. More precisely, the input to the model is a user-defined point-of-viewu0, centered in the object we are interested in segmenting. The output is the function u(x,τ). Function u(x,τ) in (2) is the segmentation function,u0 represents the initial contour (initial state of the segmentation function), and the image to segment is I0. Moreover, as proposed in [22], instead of following evolution of a particular level set of u, the PDE model follows the evolution of the entire surface of u under speed law dependent on the image gradient, without regard to any particular level set. Suitably chosen, this flow sharpens the surface around the edges and connects segmented boundaries across the missing information. In [22, 23], the authors formalized such model as the Riemannian mean curvature flow where the variability in the parameter ϵ also improves the segmentation process and provides a sort of regularization. Thus, (2) becomes
(7)∂u(x,τ)∂τ=ϵ2+|∇u|2∇⋅(g(|∇I0|)∇uϵ2+|∇u|2)τ∈[0,T],x=(x,y)∈Ωu(x,0)=u0(x)τ=0,x=(x,y)∈Ωu(τ,x)=0τ∈[0,T],x=(x,y)∈∂Ω
accompanied with initial condition u0 and zero Dirichlet boundary conditions. Regarding u0, it is usually defined as a circle completely enclosed inside the region that one wish to segment.The termg(v), called edge detector, is a nonincreasing real function such that g(v)→0 while v→∞, and it is used for the enhancement of the edges. Indeed, it controls the speed of the diffusion/regularization: if ∇u has a small mean in a neighborhood of a point x, this point x is considered the interior point of a smooth region of the image and the diffusion is therefore strong. If ∇u has a large mean value on the neighborhood of x,  x is considered an edge point and the diffusion spread is lowered, since g(v) is small for large v. A popular choice in nonlinear diffusion models is the Perona and Malik function [14]: g(v)=1/(1+v2/β), β>0. In many models, the function g(|∇u|) is replaced by its smoothed version g(|∇Gσ*u|), where Gσ is a smoothing kernel, for example, the Gauss function, which is used in presmoothing of image gradients by the convolution. For shortening notations, we will use abbreviation
(8)g=g(|∇Gσ*u|).In conclusion, we use the Riemannian mean curvature flow, as model equation of the segmentation of medical structures: given I0, the initial image and u0 equals to a circle contained inside an object of the image I0, we are interested in segmenting, we compute u(x,τ) by solving (7). The level sets of u(x,τ), at steady state, provide approximations of the contour to detect.


## 2.1. A Case Study: Image Segmentation

In this paper, we use equations (2) for image segmentation. The task of image segmentation is to find a collection of nonoverlapping subregions of a given image. In medical imaging, for example, one might want to segment the tumor or the white matter of a brain from a given MRI image.The idea behind level set (also known implicit active contours, or implicit deformable models) for image segmentation is quite simple. The user specifies an initial guess for the contour, which is then moved by image-driven forces to the boundaries of the desired objects. More precisely, the input to the model is a user-defined point-of-viewu0, centered in the object we are interested in segmenting. The output is the function u(x,τ). Function u(x,τ) in (2) is the segmentation function,u0 represents the initial contour (initial state of the segmentation function), and the image to segment is I0. Moreover, as proposed in [22], instead of following evolution of a particular level set of u, the PDE model follows the evolution of the entire surface of u under speed law dependent on the image gradient, without regard to any particular level set. Suitably chosen, this flow sharpens the surface around the edges and connects segmented boundaries across the missing information. In [22, 23], the authors formalized such model as the Riemannian mean curvature flow where the variability in the parameter ϵ also improves the segmentation process and provides a sort of regularization. Thus, (2) becomes
(7)∂u(x,τ)∂τ=ϵ2+|∇u|2∇⋅(g(|∇I0|)∇uϵ2+|∇u|2)τ∈[0,T],x=(x,y)∈Ωu(x,0)=u0(x)τ=0,x=(x,y)∈Ωu(τ,x)=0τ∈[0,T],x=(x,y)∈∂Ω
accompanied with initial condition u0 and zero Dirichlet boundary conditions. Regarding u0, it is usually defined as a circle completely enclosed inside the region that one wish to segment.The termg(v), called edge detector, is a nonincreasing real function such that g(v)→0 while v→∞, and it is used for the enhancement of the edges. Indeed, it controls the speed of the diffusion/regularization: if ∇u has a small mean in a neighborhood of a point x, this point x is considered the interior point of a smooth region of the image and the diffusion is therefore strong. If ∇u has a large mean value on the neighborhood of x,  x is considered an edge point and the diffusion spread is lowered, since g(v) is small for large v. A popular choice in nonlinear diffusion models is the Perona and Malik function [14]: g(v)=1/(1+v2/β), β>0. In many models, the function g(|∇u|) is replaced by its smoothed version g(|∇Gσ*u|), where Gσ is a smoothing kernel, for example, the Gauss function, which is used in presmoothing of image gradients by the convolution. For shortening notations, we will use abbreviation
(8)g=g(|∇Gσ*u|).In conclusion, we use the Riemannian mean curvature flow, as model equation of the segmentation of medical structures: given I0, the initial image and u0 equals to a circle contained inside an object of the image I0, we are interested in segmenting, we compute u(x,τ) by solving (7). The level sets of u(x,τ), at steady state, provide approximations of the contour to detect.


## 3. Numerical Schemes

Nonlinear PDE in (7) can be expressed in a compact way as
(9)∂u(x,y,τ)∂τ=F[u(x,y,τ,∇u(x,y,τ),I0)],
where
(10)F=ϵ2+|∇u|2∇⋅(g(|∇I0|)∇uϵ2+|∇u|2).Scale Discretization
That is discretization with respect toτ. If [0,T] is the scale interval and nscales is the number of scale steps, we denote by τi the ith scale-step for all i=1,…,nscales, so that τi+1=τi+Δτ, where Δτ=T/nscales is the step-size.Using the Euler forward finite difference scheme to discretize the scale derivative on the left hand side of (9), we get(11)u(x,y,τi)-u(x,y,τi-1)Δτ=F[u(x,y,τ)]
or, equivalently
(12)u(x,y,τi)=u(x,y,τi-1)+Δτ⋅F[u(x,y,τ)].
Let us denote as ui=u(x,y,τi), i=1,…,nscales, the function u evaluated at τi. Equation (12) is rewritten as
(13)ui=ui-1+Δτ⋅F(u(x,y,τ))≡G[u(x,y,τ)].Depending on the collocation value, used to evaluateu(x,y,τ) with respect to the parameter τ, inside the F function on the right hand side of (12) three iterative schemes derive: (i)
explicit scheme: ui=G[ui-1], that is, the function F is evaluated at ui-1=u(x,y,τi-1);(ii)
semi-implicit scheme: ui=G[ui-1,ui], that is, we use ui to discretize the numerator |∇u| of the fraction ∇u/ϵ2+|∇u|2. Other quantities are evaluated at ui-1;(iii)
implicit scheme: ui-G[u(i)]=0, that is, the function F is evaluated at ui.In summary, the difference between the semi-implicit and the implicit scheme relies on the scale discretization of the term|∇u| at the numerator of ∇u/ϵ2+|∇u|2 inside the function F. This term controls the diffusion process, and it plays the role of edge-enhancement. If we consider the three-dimensional (3D) domain ΩT=Ω×[0,T], the semi-implicit scheme employs a sort of 2D + 1 discretization of ΩT proceeding along nscales two dimensional (2D) slices each one obtained at τ≡τi, while the fully implicit scheme uses a fully 3D discretization of ΩT. This difference suggests that the fully implicit scheme may provide a more accurate edge detection than the semi-implicit scheme. This difference is highlighted by considering their discretization errors.Space Discretization
That is discretization with respect to(x,y). If Ω is the space domain, we introduce a rectangular uniform grid on Ω consisting of Nx×Ny (for simplicity we assume that Ω is a rectangular of dimension 1×1; this means that hx=1/Nx and hy=1/Ny), nodes (xi,yj)=(lΔx,mΔy),l=1,…,Nx, m=1,…,Ny, and we use finite volumes to discretize the partial derivatives of u, as in [24, 25].Scale-Space Discretization
Let(14)uil,m=u(xl,ym,τi)∈RNx×Ny×nscales
be the vector obtained from the scale-space discretization of the function u, we have the following iteration formulas.(i)
Explicit scheme:(15)uil,m-ui-1l,mΔτ=ϵ2+|∇ui-1l,m|2∇⋅(g(|∇I0|)∇ui-1l,mϵ2+|∇ui-1l,m|2)⟺uil,m=(I+Δτ[A]i-1l,m)ui-1l,m∀i=1,2,…,NE,
where, for each i, the matrix [A]i-1l,m∈ℜNx2×Ny2 and I∈ℜNx2×Ny2 is the unit matrix, while NE is the scale steps number.(i)
Semi-implicit scheme:(16)uil,m-ui-1l,mΔτ=ϵ2+|∇ui-1l,m|2∇⋅(g(|∇I0|)∇uil,mϵ2+|∇ui-1l,m|2)⟺uil,m=ui-1l,m+Δτ[A]i-1l,muil,m⟺(I+Δτ[A]i-1l,m)uil,m=ui-1l,m∀i,i=1,…,NSI
where, for each i, the matrix [A]i-1l,m∈ℜNx2×Ny2 and I∈ℜNx2×Ny2 is the unit matrix and NSI is the scale steps number.(ii)
Fully-implicit scheme:(17)uil,m-ui-1l,mΔτ=ϵ2+|∇uil,m|2∇⋅(g(|∇I0|)∇uil,mϵ2+|∇uil,m|2)⟺uil,m=ui-1l,m+Δτ[A]il,m(uil,m)∀i=1,…,NI,
where NI is the scale steps number and [A]il,m, for each i, is a nonlinear vector operator on ℜNx2×Ny2, depending on uil,m.In particular, we apply the Crank-Nicholson scheme [3] which uses the average of the forward Euler method at step i-1 and the backward Euler method at step i:
(18)uil,m-ui-1l,mΔτ=12[∇uil,mϵ2+|∇uil,m|2ϵ2+|∇uil,m|2∇⋅(g(|∇I0|)∇uil,mϵ2+|∇uil,m|2)]+12⋯⋯⋅[∇ui-1l,mϵ2+|∇ui-1l,m|2ϵ2+|∇ui-1l,m|2∇⋅(g(|∇I0|)∇ui-1l,mϵ2+|∇ui-1l,m|2)],∀i=1,…,NI⟺uil,m=ui-1l,m+Δτ[Bi-1l,mui-1l,m+[A]il,m(uil,m)],∀i=1,…,NI,
where Bi-1l,m is a matrix on ℜNx2×Ny2, depending on ui-1l,m, and [A]il,m is a nonlinear vector operator on ℜNx2×Ny2, depending on uil,m.


## 4. Algorithms and Their Computational Complexity

The effectiveness of these schemes depends on a suitable balance between accuracy (scale-space discretization error), number of flop/s per iteration (algorithm complexity), and the total execution time needed to reach a prescribed accuracy (software performance).Let us denote the discretization errorEd. It is
(19)Ed=O(hxp)+O(hyp)+O(Δτq).Explicit scheme is accurate at the first order both with respect to scale and space, that is,p=q=1; anyway, it is the one straightforwardly computable. The computational kernel is a matrix-vector product, at every scale step. This scheme requires very small time steps in order to be stable (CFL (Courant-Friedrich-Levy) condition that guarantees the stability of the evolution), and its use is limited rather by its stability than accuracy. This constraint is practically very restrictive, since it typically leads to the need for a huge amount of iterations [3]. Semi-implicit scheme is absolutely stable for all scale steps. The accuracy, in terms of discretization error with respect to both scale and space, is of the first order, because p=1, q=2 [24, 25]. Crank-Nicholson provides a discretization error of second order, that is, p=q=2, but it requires extra computations, leading to a nonlinear system of equations, at every time step, while stability is ensured for all scale steps [3]. In the following, we collect these results:
(20)explicitscheme:p=q=1,semi-implicitscheme:p=1,q=2,implicitscheme:p=q=2.Then, the fully implicit scheme provides an order of accuracy greater than that provided by the others. This difference may be important in those applications of image analysis where the edges are fundamental to recognize some pathologies.Algorithm complexity of these schemes depends on the choice of the numerical solver. Concerning the semi-implicit scheme, we employ Krylov subspace methods, which are the most effective approaches for solving large linear systems [5]. In particular, we use Generalized Minimal RESidual method (GMRES) equipped with Algebraic multigrid (AMG) preconditioner. Such techniques are convenient because they require as input only the system matrix corresponding to the finest grid. In addition, they are suitable to implement in a parallel computing environment. For the fully implicit scheme we use the Jacobian Free Newton Krylov Method (JFNK) [6]. JFNK methods are synergistic combinations of Newton-type methods for superlinearly convergent solution of nonlinear equations and Krylov subspace methods for solving the Newton correction equations. The link between the two methods is the Jacobian-vector product, which may be probed approximately without forming and storing the elements of the true Jacobian.Let us briefly describe the numerical algorithms that we are going to implement, which are based on the semi-implicit and the implicit discretization schemes, together with their complexity.Algorithm SI (Semi-Implicit Scheme)
For alli=1,…,NSI solution of
(21)(I+Δτ[A]i-1l,m)uil,m=ui-1l,m⟺HSi-1l,muil,m=ui-1l,m,
with respect to uil,m. HSi-1l,m, for each i, is a matrix ∈ℜNx2×Ny2. As space derivative we use the 2nd order finite covolume discretization scheme (see [24, 25] for convergence, consistence and stability). By this way, matrix [A]i-1l,m is a block pentadiagonal matrix with tridiagonal blocks along the main diagonal and diagonal blocks along the upper and lower diagonals.Algorithm I (Implicit Scheme)
For alli=1,…,NI solution of
(22)uil,m=ui-1l,m+Δτ[Bi-1l,mui-1l,m+[A]il,m(uil,m)]⟺HIil,m(uil,m,ui-1l,m)=0,
with respect to uil,m. HIil,m, for each i, is a nonlinear vector operator on ℜNx2×Ny2.Algorithm SI
For each scale step, to solve the linear system (21), we employ GMRES iterative method (see Algorithm 1). Computational kernel of GMRES is a matrix-vector product. Taking into account the structure of the coefficient matrix (we assume that Nx=Ny=N, then h=1/N=hx=hy), the computational cost of GMRES is
(23)TGMRES(N2)=O(kGMRESSI⋅5N2),
where kGMRESSI is the maximum iterations of GMRES (over the scale steps). Computational complexity of Algorithm SI is
(24)TSI-GMRES(N2)=O(NSI⋅kGMRESSI⋅5N2).Algorithm 1:Preconditioned GMRES for solving a linear system Ax=b. Input: A (matrix coefficient), b (right hand side), P (preconditioner). Output: xm, approximate solution at the mth step. For all iteration, a matrix-vector product is required.
Computer0=P(b-Ax0), β:=∥r0∥2 and v1:=r0/βDefine the(m+1)×mHk={hij}1≤i≤m+1,1≤j≤m. Set Hm=0.fork=1 to m doComputewj:=PAvjFori=1 to k do:hij:=(wj,vi)wj:=wj-hijviend forhj+1,j=∥wj∥2. If hj+1,j=0 set m:=j and go to 12vj+1=wj/hj+1,jend forComputeym the minimiser of ∥βe1-Hmy∥2Setxm:=x0+Vmym.Algebraic multigrid (AMG) method follows the main idea of (geometric) multigrid (MG), where a sequence of grids is constructed from the underlying geometry with corresponding transfer operators between the grids [26]. The main idea of MG is to remove the smooth error, that cannot be eliminated by relaxation on the fine grid, by coarse-grid correction. The solution process then as usual consists of presmoothing, transfer of residuals from fine to coarse grids, interpolation of corrections from coarse to fine levels, and optional postsmoothing. In contrast to geometric multigrid, the idea of AMG is to define an artificial sequence of systems of equations decreasing in size. We call these equations coarse-grid equations. The interpolation operator Plvl,m and the restriction operator Rlvl,m define the transfer from finer to coarser grids and vice versa. Finally, the operator on the coarser grid at level lv+1 is defined by
(25)Alv+1l,m=Rlvl,mAlvl,mPlvl,m.The AMG method consists of two main parts, the setup phase and the solution phase. During the setup phase, the coarse-grids and the corresponding operators are defined. The solution phase consists of a multilevel iteration. The number of recursive calls, which is the number of levelslv, depends on the size and structure of the matrix. For our case, we use the V-cycle pattern with the FALGOUT-CLJP coarse grid selection [27]. Looking at the Algorithm SI in Table 1, the preconditioner P is just Alv+1.Table 1
Comparisons between two algorithms:Ed=O(10-1).
AlgorithmΔτNumber of scale stepsdHExecution time (secs)SI0.1630.94929.262I0.410.94987.675Computational cost of each iteration of GMRES is that of the AMG preconditioner plus the matrix-vector products:(26)TAMG+GMRES(N2)=O(kGMRESSI(TAMG(lv)︸lv⋅N2+5N2)),
then, we get:
(27)TSI+AMG+GMRES(N2)=O(NSIkGMRESSIlvN2).Following picture shows a schematic description of Algorithm SI that emphasizes its main steps and the most time consuming operation, that is, the matrix vector products needed at each step of GMRES.Algorithm I
For each scale step, to solve the nonlinear equations (22), we employ the Jacobian-Free Newton-Krylov (JFNK) method. JFNK is a nested iteration method consisting of at least two and usually four levels. The primary levels, which give the method its name, are the loop over the Newton method:
(28)HIil,m(un+1l,m)=0⟺HIil,m(unl,m)+Ji(unl,m)(un+1l,m-unl,m)=0,
and the loop building up the Krylov subspace out of which each Newton step is drawn:
(29)Ji(unl,m)δunl,m=-HIil,m(unl,m),un+1l,m=unl,m+δunl,m,
Outside of the Newton loop, a globalization method is often required. We use line search.Forming each element ofJ which is a matrix of dimension N2×N2 requires taking derivatives of the system of equations with respect to u. This can be time consuming. Using the first order Taylor series expansion of HIil,m(unl,m+ρv), it follows that(30)Ji(unl,m)δunl,m=[HIil,m(unl,m+ρδunl,m)-HIil,m(unl,m)]ρ+O(ρ2),
where ρ is a perturbation. JFNK does not require the formation of the Jacobian matrix; it instead forms a result vector that approximates this matrix multiplied by a vector. This Jacobian-free approach has the advantage to provide the quadratic convergence of Newton method without the costs of computing and storing the Jacobian. Conditions are provided on the size of ρ that guarantee local convergence.Algorithm3 shows a schematic description of Algoritm I that emphasizes its main steps and the most time consuming operation, that is, evaluations of the nonlinear operator HIil,m at each innermost step.Algorithm complexity of JFNK is(31)TJFNK(N2)=O(NNewkGMRESI[f+O(N2)]),where NNew is the maximum number of Newton’s steps, over the scale steps, kGMRESI is the maximum number of GMRES iterations (computed over Newton’s steps and scale steps)), f is the number of evaluations of HIil,m. Finally, we get
(32)TI+JFNK(N2)=O(NINNewkGMRESI[f+O(N2)]).A straightforward comparison between the algorithm complexity of these algorithms shows that Algorithm SI asymptotically seems to be comparable with respect to Algorithm SI. Of course, the performance analysis must also take into account the efficiency of these two schemes in a given computing environment. Next section describes the PETSc-based implementation of these algorithms that we have developed in a multicore computing environment.


## 5. The Multicore-Based Implementation

The software has been developed using the high-level software library PETSc (Portable Extensible Toolkit for Scientific Computations) (release 3.1, March 2010) [8]. PETSc provides a suite of data structures and routines as building blocks for the implementation of large-scale codes to be used in scientific applications modeled by partial differential equation. PETSc is flexible: its modules, that can be used in application codes written in Fortran, C, and C++, are developed by means of object-oriented programming techniques.The library has a hierarchical structure: it relies on standard basic computational (BLAS, LAPACK) and communication (MPI) kernels and provides mechanism needed to write parallel application codes. PETSc transparently handles the moving of data between processes without requiring the user to call any data transfer function. This includes handling parallel data layouts, communicating ghost points, gather, scatter and broadcast operations. Such operations are optimized to minimize synchronization overheads.Our parallelization strategy is based on domain decomposition: in particular, we adopt therow-block data distribution, which is the standard PETSc data distribution. Row-block data distribution means that blocks of dimensions N2/p×N of contiguous rows of matrices of dimension N2×N2 are distributed among contiguous processes. By the same way, vectors of size N are distributed among p processors as blocks of size N/P. Such partitioning has been chosen because overheads, due to redistribution before the solution of the linear systems, are avoided. Further, row-block data distribution introduces a coarse grain parallelism which is best oriented to exploit concurrency of multicore multiprocessors because it does not require a strong cooperation among computing elements: each computing element has to locally manage the blocks that are assigned to it.The computing platform that we consider is made of 16 blades (1 blade consisting of 2 quad core Intel Xeon E5410@2.33 GHz) Dell PowerEdge M600, equipped with IEEE double precision arithmetic. Because high performance technologies can be employed in medical applications only to the extent that the overall cost of the infrastructure is affordable, and because we consider single images of medium size, we show results obtained by using 1 blade, that is, we run the parallel algorithms on up top=8 cores of a single blade. Of course, in case of multiple images or sequences of images, the use of a greater number of cores may be interesting.


## 6. Experiments

In this section we present and discuss computational results obtained by implementing Algorithm I and Algorithm SI in a multicore parallel computing machine. Before illustrating experimental results, let us briefly describe the choice of(1)
test images,(2)
comparison criteria,(3)
parameters selection.(1) Test Images
we have carried out many experiments in order to analyze the performance of these algorithms. Here we show results concerning the segmentation of a (malignant) melanoma (see Figure7) [28]. Epiluminescence microscopy (ELM) has proven to be an important tool in the early recognition of malignant melanoma [29, 30]. In ELM, halogen light is projected onto the object, thus rendering the surface translucent and making subsurface structures visible. As an initial step, the mask of the skin lesion is determined by a segmentation algorithm. Then, a set of features containing shape and radiometric features as well as local and global parameters is calculated to describe the malignancy of the lesion.  In order to better validate computed results and to analyze the software performance, we first consider a synthetic test image simulating the object we are interested in segmenting (see Figure 1).Figure 1
Test 1. Image size:840×840, simulated image and its contour to compute.(2) Comparison Criteria
We compare the algorithms using the following criteria:(a)
distance from original solution. As measure of the difference between two curves, we use the Hausdorff distance measured between the computed curve and the original one. It is well known that the Hausdorff distance is a metric over the set of all closed bounded sets (see [31]), here we restrict ourselves to finite point sets because that is all that is necessary for segmentation algorithms [32]. Given two finite point sets C1 and C2, the Hausdorff distance dH between the sets C1 and C2 is defined as follows:
(33)dH(C1,C2)=max⁡{h(C1,C2),h(C2,C1)},h(C1,C2)=max⁡a∈C1min⁡b∈C2‖a-b‖,
where ∥·∥ is the euclidean norm. It identifies the point a∈C1 that is fastest from any point of C2 and vice versa, then it keeps the maximum,(b)
efficiency: execution time of (serial) algorithms,(c)
convergence history: behavior of residuals and iteration numbers of inner solvers,(d)
Parallel performance: execution time, speedup, and efficiency versus cores number.(3) Parameters Selection
We setK=1.0 and ϵ=1.0. Let us explain how we select the values of the scale step size and the number of scale steps. Regarding Δτ, its value is chosen according to that required to Ed. Taking into account that, in Algorithm SI, Ed is accurate at the first order with respect to ΔτSI, while it is accurate at the second order with respect to ΔτI, in Algorithm I, by requiring that the discretization error is about the same, we get
(34)Ed=O(ΔτSI)=O(ΔτI2)⟹ΔτSI=ΔτI.Finally, in Algorithm SI, the stopping criterion of the linear solver (GMRES) uses the tolerance(35)TOL=10-10,
while the preconditioner AMG uses the tolerance TOL=10-7 and 25 as maximum number of AMG-levels.  In the Algorithm 2, the stopping criterion of the nonlinear solver uses the tolerance
(36)TOL=10-10.Algorithm 2:Sketch of Algorithm SI.
fori=1 to NSI dofork=1 to kGMRES doforlv=1 to levelsAMG domatrix-vector products involvingHSil,m, Alv+1l,m and uil,m(k,lv)end forAlgorithm 3:Sketch of Algorithm I.
fori=1 to NI doforn=1 to NNew dofork=1 to kGMRES doevaluateHIil,m, uil,m(n,k)end forRegarding the number of scale steps (NSI and NI), taking into account that
(37)NSI,I=TΔτSI,I,
its choice depends on ΔτSI,I and on the value of τ≡T, that is, the value of the scale parameter corresponding to steady state of the segmentation function u(x,τ), solution of the PDE model. To check the steady state, we require that the residuals, corresponding to different scale steps, reach the tolerance (38)TOL=10-9.We found that this corresponds toT=0.4 (see Figure 2) for Test 1 and to T=2 for Test 2.Figure 2
Test 1: The 3D visualization of the segmentation functionu(τ,x) at steady state T=0.4.Test 1: Synthetic Image
In Tables1, 2, and 3, we show the Hausdorff distance and execution time by requiring that discretization error is of the first, second, and third order, that is, Ed=O(10-1), Ed=O(10-2), and Ed=O(10-3), respectively. Hence, we get the following values of the scale step size:
(39)ΔτI=0.4⟹ΔτSI=0.16,ΔτI=0.04⟹ΔτSI=0.016,ΔτI=0.004⟹ΔτSI=0.0016.Table 2
Comparisons between two algorithms:Ed=O(10-2).
AlgorithmΔτNumber of scale stepsdHExecution time (secs)SI0.016250.641261.487I0.04100.949871.67Table 3
Comparisons between two algorithms:Ed=O(10-3).
AlgorithmΔτNumber of scale stepsdHExecution time (secs)SI0.00162500.5993356.926I0.0041000.9492356.335Moreover, concerning the number of scale steps, it follows that(40)Ed=O(10-1)⟹NI=0.40.4=1,NSI=int⁡[0.40.16]=3,Ed=O(10-2)⟹NSI=0.40.04=10,NSI=0.40.016=25,Ed=O(10-3)⟹NSI=0.40.04=100,NSI=0.40.016=250.Note that while in the first case, that is, if we requireEd=O(10-1), these two algorithms are quite numerically equivalent, both in terms of execution time and of the computed result; as discretization error decreases, Algorithm I appears to be more robust than Algorithm SI, in the sense that Algorithm I reaches the steady state with high accuracy (the Hausdorff distance is of 95%), while the computed results of Algorithm SI are less accurate, even though the execution time of Algorithm I sometimes slightly increases. These results suggest that if it needs to get an accurate and reliable result, Algorithm I should be preferable. Figure 3 show segmentation results. Finally, note that the execution time of these two algorithms asymptotically is the same, as stated by the analysis of computational cost carried on in Section 4.Figure 3
Test 1: Comparisons between segmentation results. On the left Algorithm I. On the right Algorithm SI. First row:ΔτI=0.4, ΔτSI=0.16. Second row: ΔτI=0.04, ΔτSI=0.16. Third row: ΔτI=0.004,  ΔτI=0.0016.Convergence History
Convergence history is illustrated by showing the behavior of relative residuals versus the scale steps (see Figures4, 5, and 6), and by reporting iteration number of the GMRES and of Newton’s method, respectively, (see Tables 4, 5, 6, and 7). We consider ΔτSI=0.16, 0.016, and ΔτI=0.04.Table 4
Convergence history of Algorithm SI.ΔτSI=0.16. First column reports the scale step number, second column reports the number of GMRES iterations at each scale step, and last column reports the maximum number of AMG levels at each GMRES iteration.
ikgmresilvi186265353Table 5
Convergence history of Algorithm SI.ΔτSI=0.016. On the first column, we denote the scale step number, and the symbol p-q is used to denote the steps ranging from the p-th until to the q-th. Second column reports the number of GMRES iterations at each scale step, third column reports the maximum number of AMG levels at each GMRES iteration.
ikgmresilvi1562-3554545–7448-84310–133314–203221-223123–2521Table 6
Convergence history of Algorithm SI.ΔτI=0.4. First column reports the scale step number, second column reports the number of Newton's steps at each scale step, and last column reports the number of GMRES iterations at each Newton's step.
iNNewikGMRESi1109, 9, 8, 8, 8, 7, 7, 7, 6, 6Table 7
Convergence history of Algorithm SI.ΔτI=0.04. First column reports the scale step number, second column reports the number of Newton’s steps at each scale step, and last column reports the number of GMRES iterations at each Newton's step.
iNNewikGMRESi1109, 9, 8, 8, 7, 7, 7, 7, 6, 62109, 8, 8, 8, 7, 7, 7, 6, 6, 6399, 8, 8, 7, 7, 6, 6, 6, 6498, 8, 8, 7, 7, 6, 6, 6, 5598, 8, 7, 7, 7, 6, 6, 5, 5678, 7, 7, 7, 6, 5, 5777, 7, 7, 7, 6, 5, 5867, 7, 6, 6, 6, 5967, 7, 6, 6, 5, 51057, 6, 6, 5, 4Figure 4
Algorithm SI. Behavior of the relative residual∥ri∥2/∥r0∥2 versus 3 scale steps. ΔτSI=0.16.Figure 5
Algorithm SI. Behavior of the relative residual∥ri∥2/∥r0∥2 versus 25 scale steps. ΔτSI=0.016.Figure 6
Algorithm I. Behavior of the relative residual∥ri∥2/∥r0∥2 versus 10 scale steps. ΔτI=0.04.Figure 7
Test 2: ELM of a melanoma. Image size is840×840.In the following, we show results corresponding to the segmentation of a melanoma (see Figure7). As expected, because this is a real image, the steady state is reached at a scale greater than that of the synthetic test image, that is, T=2, thus both algorithms require a greater number of scale steps to the reach the steady state. Tables 8, 9, and 10, and Figure 8 compare results of Algorithm SI and Algorithm I.Table 8
Test 2: Comparisons between two algorithms:Ed=O(10-1).
AlgorithmΔτNumber of scale stepsExecution time (secs)SI0.161341.125I0.4547.893Table 9
Test 2: Comparisons between two algorithms:Ed=O(10-2).
AlgorithmΔτNumber of scale stepsExecution time (secs)SI0.016125417.56I0.0450455.33Table 10
Test 2: Comparisons between two algorithms:Ed=O(10-3).
AlgorithmΔτNumber of scale stepsExecution time (secs)SI0.001612504828.5I0.0045004663.7Figure 8
Test 2: Comparisons between segmentation results. On the left Algorithm I. On the right Algorithm SI. First row:,ΔτI=0.4, ΔτSI=0.16. Second row: ΔτI=0.04, ΔτSI=0.16. Third row: ΔτI=0.004,ΔτI=0.0016.Parallel Performance
We show the performance of the multicore-based parallel algorithms and their scalability as the number of cores increases. We run the parallel algorithms using up top=8 cores of the parallel machine.Following Figures report execution time, speedup, and efficiency of Algorithm SI, at scale step sizeΔτ=0.16 (i.e., Ed=(10-1) (i.e., Figures 9, 10 and 11), then same results are shown at scale step size Δτ=0.016, corresponding to Ed=O(10-2) (i.e., Figures  12, 13 and 14).Figure 9
Test 1: Algorithm SI: Total execution time versus the number of cores.Δτ=0.16. 3 scale steps.Figure 10
Test 1: Algorithm SI: Speedup versus the number of cores.Δτ=0.16. 3 scale steps.Figure 11
Test 1: Algorithm SI: Parallel efficiency versus the number of cores.Δτ=0.16. 3 scale steps.Figure 12
Test 1: Algorithm SI: Total execution time versus the number of cores.Δτ=0.016.Figure 13
Test 1: Algorithm SI: Speedup versus the number of cores.Δτ=0.016.Figure 14
Test 1: Algorithm SI: Parallel efficiency versus the number of cores.Δτ=0.016.Finally, we report execution time, speedup, and efficiency of Algorithm I, at scale stepΔτ=0.4 (i.e., Figures 15, 16 and 17) and Δτ=0.04 (i.e., Figures 18, 19, and 20), respectively. Note that parallel efficiency of both algorithms always is, at least, of 60% and, on average, of about 80%. In particular, parallel efficiency of Algorithm I is about 90%. Execution time of both algorithms reduces to about 2 seconds on eight cores in the first case (i.e., Ed=(10-1), and to about 10 seconds in the second case (Ed=O(10-2)). This means that, in a multicore computing environment, both algorithms provide the requested solution within a response time that can be considered quite acceptable in medical imaging applications and, in particular, that Algorithm I is competitive with Algorithm SI.Figure 15
Test 1: Algorithm I: Total execution time versus the number of cores.ΔτI=0.4.Figure 16
Test 1: Algorithm I: Speedup.ΔτI=0.4.Figure 17
Test 1: Algorithm I: Efficiency.ΔτI=0.4.Figure 18
Test 1: Algorithm I: Total execution time versus the number of cores.ΔτI=0.04.Figure 19
Test 1: Algorithm I: Speedup.ΔτI=0.04.Figure 20
Test 1: Algorithm I: Efficiency.ΔτI=0.04.Figures21, 22, 23, 24, 25, and 26 show time, speedup, and parallel efficiency of two algorithms in case of Test 2. We only consider the case of ΔτSI=0.16 and ΔτI=0.4.Figure 21
Test 2: Algorithm SI: Total execution time versus the number of cores.ΔτI=0.16.Figure 22
Test 2: Algorithm SI: Speedup.ΔτI=0.16.Figure 23
Test 2: Algorithm SI: Efficiency.ΔτI=0.16.Figure 24
Test 2: Algorithm I: Total execution time versus the number of cores.ΔτI=0.4.Figure 25
Test 2: Algorithm I. Speed up.ΔτI=0.4.Figure 26
Test 2: Algorithm I. Efficiency.ΔτI=0.4.Finally, we show results on scalability of parallel algorithms. LetTp(N) be the execution time of the parallel algorithm running on p cores for segmenting an image of size N2×N2. We measure the scalability of these algorithms by measuring Tp(N) as N varies, once p is fixed, and by measuring Tp(N) as N and p grow. We note that, in case of Algorithm SI, the scaling factor is
(41)Tp(2N)Tp(N)≃4.2,
while, for Algorithm I, we get as scaling factor
(42)Tp(2N)Tp(N)≃2.3.This means that Algorithm I scales better than Algorithm SI. In Figures27 and 28 (for Algorithm SI) and Figure 29 and 30 (for Algorithm SI), we report Tp(N) as N and p varies. In particular, each point of the graph refers to the execution time of the parallel algorithm at (p=k·p1,N=k·N1), where p1=1, N1=210, and k=1,2,3,4.Figure 27
Scalability of parallel Algorithm SI, asN is fixed and p varies, and N=210,420,840,1680. Each line refers to the execution time of the algorithm at a fixed value of N.Figure 28
Scalability of parallel Algorithm SI, asN and p vary. Each point of the graph refers to the execution time of the parallel semi-implicit algorithm at (p=k·p1, N=k·N1), where p1=1, N1=210, and k=1,2,3,4.Figure 29
Scalability of parallel Algorithm I, asN is fixed and p varies, and N=210,420,840,1680. Each line refers to the execution time of the algorithm at a fixed value of N. ΔτI=0.4.Figure 30
Scalability of parallel Algorithm I, asN and p vary. Each point of the graph refers to the execution time of the parallel implicit algorithm at (p=k·p1, N=k·N1), where p1=1, N1=210 and k=1,2,3,4. ΔτI=0.4.


## 7. Conclusions

A straightforward comparison between the semi-implicit and the fully implicit discretization schemes of nonlinear PDE of parabolic/hyperbolic type states that fully implicit discretization usually leads to too expensive algorithms. In this paper, we provide a multicore implementation of two numerical algorithms arising from using these two discretization schemes: semi-implicit (Algorithm SI) and fully implicit (Algorithm I). Taking into account that we aim to solve such problems on parallel computer in a scalable way, in the first case, we use, as linear solver, Krylov iterative methods (GMRES) with algebraic multigrid preconditioners (AMG). Regarding the fully implicit scheme, we use the Jacobian-Free-Newton Krylov (JFNK) method as nonlinear solver.We compare these two algorithms using different metrics measuring both the accuracy and the efficiency.  We note that if we require that the discretization errorEd is Ed=O(10-1), these two algorithms are quite numerically equivalent, both in terms of execution time and of the computed result; while, as discretization error decreases, Algorithm I appears to be more robust than Algorithm SI, in the sense that Algorithm I reaches the steady state with high accuracy (the Hausdorff distance is of 95%), while the computed results of Algorithm SI are less accurate, even though the execution time of Algorithm I sometimes slightly increases. These results suggest that if it needs to get accurate and reliable results, Algorithm I should be preferable.The parallel efficiency of both algorithms always is, at least, of60% and, on average, of about 80%. In particular, parallel efficiency of Algorithm I is of about 90%. Execution time of both algorithms reduces to about 2 seconds on eight cores if Ed=(10-1) and to about 10 seconds if Ed=O(10-2). This means that, in a multicore computing environment, Algorithm I is competitive with Algorithm SI.In conclusion, our results suggest that if it is required high accuracy of the computed solution in a suitable turnaround time, using a multicore computing environment fully implicit scheme provides an accurate and reliable solution within a response time of few seconds, quite acceptable in medical imaging applications, such as computer-aided-diagnosis.



---

*Source: 680765-2011-10-10.xml*