high - quality images are crucial to accurately diagnose a patient or determine treatment . in addition to requiring the best images possible
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 .
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 .
we consider two standard discretization schemes of nonlinear time - dependent pdes : semi - implicit scheme and fully implicit scheme .
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 .   in recent years
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 .
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 section 2 , an overview of the pde model equation used in describing some of inverse problems in imaging applications
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 .
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 .
most medical imaging problems lead to ill - posed ( inverse ) problems in the sense of hadamard [ 911 ] .
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 , 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 )   to denoise images .
moreover , the authors derive the euler - lagrange equations as a time - dependent pde . in the same years
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 image u(x ) , x   should be embedded into a one - parameter family of derived images , in which fine - scale structures are successively suppressed : 
						 ( 1)ss: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 , defined as follows . 

nonlinear diffusion modelslet x = ( x , y )   and u(x ,  ) , defined in [ 0 , t ]   , be the scale - space representation of the brightness function image u(x ) defined in    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 ) . 
 
let x = ( x , y )   and u(x ,  ) , defined in [ 0 , t ]   , be the scale - space representation of the brightness function image u(x ) defined in    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 .
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 ) [ 1820 ] .
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   which proved existence and uniqueness of viscosity solutions . in this paper
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 - view u0 , centered in the object we are interested in segmenting .
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 i. moreover , as proposed in , 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|)u2+|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 term g(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 : g(v ) = 1/(1 + v/ ) ,  > 0 . in many models , the function g(|u| )
is replaced by its smoothed version g(|gu| ) , 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(|gu|).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 .
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|)u2+|u|2 ) . 

scale discretizationthat 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 . 
 
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,i1)=f[u(x , y, ) ] 
					 or , equivalently 
						 ( 12)u(x , y,i)=u(x , y,i1)+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 = ui1+f(u(x , y,))g[u(x , y, ) ] . 
				
depending on the collocation value , used to evaluate u(x , y ,  ) with respect to the parameter  , inside the f function on the right hand side of ( 12 ) three iterative schemes derive :  
 explicit scheme : ui = g[ui1 ] , that is , the function f is evaluated at ui1 = u(x , y , i1 ) ; 
 semi - implicit scheme : ui = g[ui1 , ui ] , that is , we use ui to discretize the numerator |u| of the fraction u/2+|u|2 .
other quantities are evaluated at ui1 ; 
 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 .
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 ) = ( lx , my),l = 1 ,  , nx , m = 1 , 
, ny , and we use finite volumes to discretize the partial derivatives of u , as in [ 24 , 25 ] . 
  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 ) = ( lx , my),l = 1 ,  , nx , m = 1 , 
, ny , and we use finite volumes to discretize the partial derivatives of u , as in [ 24 , 25 ]
. 
 scale - space discretizationlet 
							 ( 14)uil , m = u(xl , ym,i)nxnynscales 
						 be the vector obtained from the scale - space discretization of the function u , we have the following iteration formulas.explicit scheme : 
										 ( 15)uil , mui1l , m=2+|ui1l , m|2  (g(|i0|)ui1l , m2+|ui1l , m|2)uil , m=(i+[a]i1l , m)ui1l , m  i=1,2,  ,ne , 
									 where , for each i , the matrix [ a]i1  ny and i  ny is the unit matrix , while ne is the scale steps number . let 
							 ( 14)uil , m = u(xl , ym,i)nxnynscales 
						 be the vector obtained from the scale - space discretization of the function u , we have the following iteration formulas.explicit scheme : 
										 ( 15)uil , mui1l , m=2+|ui1l , m|2  (g(|i0|)ui1l , m2+|ui1l , m|2)uil , m=(i+[a]i1l , m)ui1l , m  i=1,2,  ,ne , 
									 where , for each i , the matrix [ a]i1  ny and i  ny is the unit matrix , while ne is the scale steps number .
explicit scheme : 
										 ( 15)uil , mui1l , m=2+|ui1l , m|2  (g(|i0|)ui1l , m2+|ui1l , m|2)uil , m=(i+[a]i1l , m)ui1l , m  i=1,2,  ,ne , 
									 where , for each i , the matrix [ a]i1  ny semi - implicit scheme : 
								 ( 16)uil , mui1l , m=2+|ui1l , m|2  (g(|i0|)uil , m2+|ui1l , m|2)uil , m = ui1l , m+[a]i1l , muil , m(i+[a]i1l , m)uil , m = ui1l , m  i , i=1,  ,nsi 
							 where , for each i , the matrix [ a]i1  ny fully - implicit scheme : 
								 ( 17)uil , mui1l , m=2+|uil , m|2  (g(|i0|)uil , m2+|uil , m|2)uil , m = ui1l , m+[a]il , m(uil , m )  i=1,  ,ni , 
							 where ni is the scale steps number and [ a]i , for each i , is a nonlinear vector operator on ny           in particular , we apply the crank - nicholson scheme   which uses the average of the forward euler method at step i  1 and the backward euler method at step i : 
						 ( 18)uil , mui1l , m=12[2+|uil , m|2  (g(|i0|)uil , m2+|uil , m|2)]+12  [2+|ui1l , m|2  (g(|i0|)ui1l , m2+|ui1l , m|2 ) ] ,  i=1,  ,niuil , m = ui1l , m+[bi1l , mui1l , m+[a]il , m(uil , m ) ] ,  i=1,  ,ni , 
					 where bi1 is a matrix on ny
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 ) .         
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 .
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 .
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 . in the following , we collect these results : 
						 ( 20)explicit  scheme :  p = q=1,semi - implicit  scheme :  p=1 ,  q=2,implicit  scheme :  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 .      
concerning the semi - implicit scheme , we employ krylov subspace methods , which are the most effective approaches for solving large linear systems .
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 ) .
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 all i = 1 ,  , nsi solution of 
							 ( 21)(i+[a]i1l , m)uil , m = ui1l , mhsi1l , muil , m = ui1l , m , 
						 with respect to ui .
hsi1 , for each i , is a matrix ny . 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]i1 is a block pentadiagonal matrix with tridiagonal blocks along the main diagonal and diagonal blocks along the upper and lower diagonals . for all i = 1 ,  , nsi solution of 
							 ( 21)(i+[a]i1l , m)uil , m = ui1l , mhsi1l , muil , m = ui1l , m , 
						 with respect to ui .
hsi1 , for each i , is a matrix ny 
 algorithm i ( implicit scheme)for all i = 1 ,  , ni solution of 
							 ( 22)uil , m = ui1l , m+[bi1l , mui1l , m+[a]il , m(uil , m)]hiil , m(uil , m , ui1l , m)=0 , 
						 with respect to ui .
hii , for each i , is a nonlinear vector operator on ny . for all i = 1 ,  , ni solution of 
							 ( 22)uil , m = ui1l , m+[bi1l , mui1l , m+[a]il , m(uil , m)]hiil , m(uil , m , ui1l , m)=0 , 
						 with respect to ui .
hii , for each i , is a nonlinear vector operator on ny 
 algorithm sifor each scale step , to solve the linear system ( 21 ) , we employ gmres iterative method ( see algorithm 1 ) .
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(kgmressi5n2 ) , 
						 where kgmres is the maximum iterations of gmres ( over the scale steps ) .
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(kgmressi5n2 ) , 
						 where kgmres is the maximum iterations of gmres ( over the scale steps ) .
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 .
the main idea of mg is to remove the smooth error , that can not 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 .
the interpolation operator plv and the restriction operator rlv 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 solution phase consists of a multilevel iteration . the number of recursive calls , which is the number of levels lv , 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 .
looking at the algorithm si in table 1 , the preconditioner p is just alv+1 .   
computational 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)lvn2 + 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 ifor 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)=0hiil , m(unl , m)+ji(unl , m)(un+1l , munl , 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 .
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)=0hiil , m(unl , m)+ji(unl , m)(un+1l , munl , 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 .
forming each element of j which is a matrix of dimension n  n 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 hii(un + 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 .
algorithm 3 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 hii 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 , kgmres is the maximum number of gmres iterations ( computed over newton 's steps and scale steps ) ) , f is the number of evaluations of hii . 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 .
the software has been developed using the high - level software library petsc ( portable extensible toolkit for scientific computations ) ( release 3.1 , march 2010 ) .
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 .
our parallelization strategy is based on domain decomposition : in particular , we adopt the row - block data distribution , which is the standard petsc data distribution .
n / p  n of contiguous rows of matrices of dimension n  n 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 to p = 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 .
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 ,
( 1 ) test imageswe 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 figure 7 ) .
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 ,
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 ) . 
 
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 figure 7 ) .
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 ,
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 ) . 

( 2 ) comparison criteria     we compare the algorithms using the following criteria : 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 ) , here we restrict ourselves to finite point sets because that is all that is necessary for segmentation algorithms . 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)=maxac1  minbc2||ab|| , 
									 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 , efficiency : execution time of ( serial ) algorithms , convergence history : behavior of residuals and iteration numbers of inner solvers , parallel performance : execution time , speedup , and efficiency versus cores number . 
     
we compare the algorithms using the following criteria : 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 ) , here we restrict ourselves to finite point sets because that is all that is necessary for segmentation algorithms . 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)=maxac1  minbc2||ab|| , 
									 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 , efficiency : execution time of ( serial ) algorithms , convergence history : behavior of residuals and iteration numbers of inner solvers , parallel performance : execution time , speedup , and efficiency versus cores number . 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 ) , here we restrict ourselves to finite point sets because that is all that is necessary for segmentation algorithms .
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)=maxac1  minbc2||ab|| , 
									 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 ,   efficiency : execution time of ( serial ) algorithms , convergence history : behavior of residuals and iteration numbers of inner solvers , parallel performance : execution time , speedup , and efficiency versus cores number . 

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 . 
 
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=1010 , 
					 while the preconditioner amg uses the tolerance tol = 10 and 25 as maximum number of amg - levels .   in the algorithm 2
, the stopping criterion of the nonlinear solver uses the tolerance 
						 ( 36)tol=1010 . 
				
regarding the number of scale steps ( nsi and ni ) , taking into account that 
						 ( 37)nsi , i = tsi , 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  
				 we found that this corresponds to t = 0.4 ( see figure 2 ) for test 1 and to t = 2 for test 2 .  

test 1 : synthetic imagein tables 1 , 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 ) , ed = o(10 ) , and ed = o(10 ) , respectively .
hence , we get the following values of the scale step size : 
							 ( 39)i=0.4si=0.16,i=0.04si=0.016,i=0.004si=0.0016 . 
  in tables 1 , 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 ) , ed = o(10 ) , and ed = o(10 ) , respectively .
hence , we get the following values of the scale step size : 
							 ( 39)i=0.4si=0.16,i=0.04si=0.016,i=0.004si=0.0016 . 
				
moreover , concerning the number of scale steps , it follows that  
 ( 40)ed = o(101)ni=0.40.4=1 ,  nsi = int[0.40.16]=3,ed = o(102)nsi=0.40.04=10 ,  nsi=0.40.016=25,ed = o(103)nsi=0.40.04=100 ,  nsi=0.40.016=250 . 
 
				
note that while in the first case , that is , if we require ed = o(10 ) , 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 .  
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 . 

convergence historyconvergence history is illustrated by showing the behavior of relative residuals versus the scale steps ( see figures 4 , 5 , and 6 ) , and by reporting iteration number of the gmres and of newton 's method , respectively , ( see tables 4 , 5 , 6 , and 7 ) .
convergence history is illustrated by showing the behavior of relative residuals versus the scale steps ( see figures 4 , 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 .   in the following ,
we show results corresponding to the segmentation of a melanoma ( see figure 7 ) . 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.  
 parallel performancewe 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 to p = 8 cores of the parallel machine . 
 
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 to p = 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 ) ( i.e. , figures 9 , 10 and 11 ) ,
then same results are shown at scale step size  = 0.016 , corresponding to ed = o(10 ) ( i.e. , figures   12 , 13 and 14 ) . 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% .
execution time of both algorithms reduces to about 2 seconds on eight cores in the first case ( i.e. , ed = ( 10 ) , and to about 10 seconds in the second case ( ed = o(10 ) ) .
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 .   figures 21 , 22 , 23 , 24 , 25 , and 26 show time , speedup , and parallel efficiency of two algorithms in case of test 2 .
finally , we show results on scalability of parallel algorithms . let tp(n ) be the execution time of the parallel algorithm running on p cores for segmenting an image of size n  n. 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 figures 27 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 .
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 error ed is ed = o(10 ) , 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 , of 60% and , on average , of about 80% .
execution time of both algorithms reduces to about 2 seconds on eight cores if ed = ( 10 ) and to about 10 seconds if ed = o(10 ) .
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 .
 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 . 
1. Introduction
2. Diffusion Models Arising in Medical Imaging
3. Numerical Schemes
4. Algorithms and Their Computational Complexity
5. The Multicore-Based Implementation
6. Experiments
7. Conclusions