The purpose of this section is to describe the registration module of our developed system to implement and illustrate the pertinence of our proposed anatomical structure constraints based non-rigid registration method [15].
5.4.1 Non-rigid Registration Method with Anatomical Struc-ture Constraints
As already explained above, to evaluate LT of HCC, we need to register the treated region onto the tumor using our proposed non-rigid registration method. The objective of image registration is to find the optimal transformation parametersµ, which maps every point in the postoperative liver (after LT) which is called the
moving imageM onto its corresponding point in the preoperative liver (before LT) which is called the fixed imageR, so that the similarity function, which measures the difference between M and R can be minimized. Then, it can automatically transform the treated region to the optimal position using the same transformation parameters.
The input images are given as two 3D discrete signals: the moving imageM and the fixed imageRwith intensitiesfm(x) andfr(x), respectively, wherex∈I ⊂ZN, and I is a 3D discrete interval representing the set of all voxel coordinates in the image. The transformed moving image W can be obtained by using the optimal transformation parameters µ, that can be defined as fw(x) = fm(x;µ).
Nevertheless the majority of the image registration methods can be viewed as a combination of the four distinctive components: transformation, interpolation, the similarity metric based on cost function, optimization. Among the four compo-nents, the similarity metric based on cost function is the most important significant critical for the registration algorithm, as it defines the goal of optimization and measures how well the transformed moving image W matches the fixed image R.
In the following, the non-rigid registration algorithm is described in detail.
A. Similarity metric
As explained above, we estimate the parameter of transformation by optimizing a cost function (similarity metric) in the processing of registration. In general, the solution to the registration problem can be computed by optimizing the following cost function [16, 17, 18]:
E(µ) =Esim(R, W) = Esim(fr, fm·T(µ)) =Esim(fr(x), fm(x;µ)) (5.1) According to the intensity-based registration process, the aim is to minimize (or maxmize) the similarityEsim between the two images. There are several image metrics that are used in voxel similarity-based image registration [5, 6, 7, 8, 9]: cor-relation coefficient, sum of squared differences (SSD) or mutual information (MI).
SSD is one of the simplest voxel similarity based cost functions which should be minimized during registration. SSD is based on the strict assumption that the two images only differ from the Gaussian noise, but this is not true for medical images, especially for multi-modal images. Therefore, this method can be only used for mono-modal registration. However, MI [10, 11] is the most widely used for medical image similarity measures that no limit to the image content about the modalities (mono-modal and multi-modal). MI can measure the statistical dependence or information redundancy between the image intensities of corresponding voxels in both images, which is assumed to be maximal if the images are registered.
However, the biggest serious drawback of MI is the absence of spatial infor-mation [12, 16, 19, 20], due to it’s definition being based on Shannon’s Entropy,
which assumes each pixels is independent of its neighbors. Hence, intensity rela-tionships in one region can occasionally cause errors to occur in another region where the intensity relationships are vary considerable. In [13], Loeckx et al. de-scribed that MI-based non-rigid registration might reduce smaller image details or have problems with spatially-varying intensity inhomogeneity due to a bias field.
The classical non-rigid registration based on intensity exhibits an inherent weakness in that it destroys internal structures of the volume. Hence, in this study, we proposed a novel non-rigid registration algorithm using an anatomical structure term to constrain the deformation, which promised to outperform the classical intensity-based non-rigid registration method in terms of maintaining the anatomical structural stability.
Based on intensity similarity metric, a different approach is used to constrain the anatomical structure to be stable by adding a landmark penalty term. On the basis of minimizing the cost function, the proposed similarity metric MI-LC is defined as:
E(µ) = Esim(R, W) +λLC(Lr, Lw)
=Esim(fr, fm·T(µ)) +λLC(Lr, Lm·T(µ))
=Esim(fr(x), fm(x;µ)) +λLC(Lr, Lm;µ)
(5.2) whereLr and Lm are the fixed landmarks and moving landmarks, respectively.
Lw =Lm;µ is the corresponding transformed moving landmarks. The first term is the driving force behind the intensity-based registration process, we can evaluate the similarity by computing SSD or MI metric. Due to SSD having limitations, it can be only used for the mono-modal registration. Hence, for future applications, we apply MI into our method for the multi-modal images. For the second term, LCis a landmark penalty term with the weightλused to constrain the anatomical structure. The definition of mutual information would be explained in [12].
Considering that a landmark penalty term is introduced in the cost function for guaranteeing the anatomical structure stability. To be reliable, the landmarks are taken as being clearly discernible geometrical features of the medical image. In this study, we have have had the following four sets of landmarks between before LT and after LT manually defined by the medical experts. These landmarks are shown in Fig 5.3.
Combining the anatomical structure constraint term into the mutual informa-tion (MI), the cost funcinforma-tion is as follows:
E(µ) =−IM I(R, W) +λLC(Lr, Lw)
=−IM I(fr(x), fm(x;µ)) +λLC(Lr, Lm;µ) (5.3) Where IM I is the mutual information value between two livers before LT and after LT.
We introduce an additional information based on the anatomical structure constraint term that can be written into:
Figure 5.3: Illustration of landmark positions before LT and after LT. (a) the branch of the right hepatic (the first one is before LT and last one is after LT);
(b) the umbilical portion; (c) the first branch of the portal vein; (d) the second branch of the portal vein.
LC(Lr, Lw) = n1
L
nL
P
i=1
[(Lr)i−(Lw)i]2 = n1
L
nL
P
i=1
(Lr)i−(Lm;µ)i2
(5.4) Our proposed similarity metric (MI-LC) is defined using the above concepts as:
E(µ) =−IM I(fr(x), fm(x;µ)) + n1
L
nL
P
i=1
(Lr)i−(Lm;µ)i2
(5.5) Notice that with our proposed method we are not neglecting the spatial in-formation of the image, as a landmark penalty term is introduced into the 3D intensity distributions and optimized together. The similarity calculation is based on which kind of interpolation method is used, so the B-spline interpolation [20]
is adopt in this work.
B. Transformation
In general, depending on the global deformation alone is not sufficient for image registration. Therefore, a method that can capture the local deformation was needed. To achieve better accuracy, the transformation in our algorithm comprised
of a global and a local transformation. Our work aimed to attain a local warp that is not affected the stability of the structure. The global transformation describe the global or large motions on the two volumes, such as rotation, translation and scaling; while the local transformation describe the local or detailed motions, such as deformation of the tissues. The parameters of global transformation are first determined by registration and the results are the initial values to estimate the parameters of local transformation. By such a combination, the deformation is estimated after the large differences on the volumes are eliminated. This not only ensures the estimated deformation is more reasonable, but also makes the registration to end earlier.
T(x) = TGobal(x) +TLocal(x) (5.6) where x= [x, y, z]T is the coordinate of a 3D point.
The simplest transformation of the global deformation method is the affine transformation [4, 6, 7, 8, 9]. For the 3D case, it can be definite by:
TGobal(x) =Ax+t=
a1 a2 a3 a4 a5 a6 a7 a8 a9
x y z
+
tx ty tz
(5.7)
whereAis the linear transformation matrix. tis the translation vector [tx, ty, tz]T along each axis. Hence, the global transformation is parameterized by 12 degrees of freedom.
B-spline functions [3, 4, 21] have been widely used to represent deformations, because of their compact support, computational simplicity, good approximation properties, and implicit smoothness. Hence, in our registration process, the cubic B-spline function based free form deformation (FFD) model [3, 5, 19] is chose to describe the localized characteristics or detailed motions.
The FFD model is parameterized by the coefficients of a set of sparse, uniformly-spaced control points. We denote Φ as the lattice which is an integer grid of control points with uniform spacing ρ = (ρx, ρy, ρz). Given the coefficients of the control point are denoted as κi,j,k, the deformation of each point can be calculated from these coefficients by cubic B-spline interpolation, according to Equation 5.8.
Let TLoacl(x;κ) be a B-spline transformation function. It is represented as:
Tloacl(x;κ) = x+
3
P
l=0 3
P
m=0 3
P
n=0
Bn(u)Bm(v)Bl(w)κi+n,j+m,k+l (5.8) where i=bx/ρxc −1,j =by/ρyc −1, k=bz/ρzc −1, andu=x/ρx− bx/ρxc, v =y/ρy− by/ρyc, w=z/ρz− bz/ρzc.
Bl,Bm,Bnare the basis function of the cubic B-spline kernel and the definition are given by:
B0(u) = (−u3+3u62−3u+1) B1(u) = (3u3−6u6 2+4) B2(u) = (−3u3+3u62+3u+1) B3(u) = u63
(5.9)
κ can be regarded as the parameters of the local transformation and the de-grees of non-rigid deformation which can be modeled depends essentially on the resolution of the mesh of control points. A small spacing of control points allowed modeling of highly local non-rigid deformations. In our application, the spacing on each axis is equal to 16-voxel-width, so the number of the control points ap-proximate 4096 and the number of the parameters is about 12288. According to Equation 5.8 and Equation 5.9, the set of parameters for the transformation is µ ={A,t, κ}. Using this parameters of transformation, the transformed moving image fw(x) = fm(T(x;µ)) can be computed.
C. Optimization
In our application, in order to minimize the cost function E(µ), the L-BFGS-B optimization method [26] are applied to find the optimal transformation parame-ters. The parameters of global and local transformation are optimized separately.
Initially, rigid registration is used, the parameters of affine transformation are adjusted by L-BFGS-B optimization method. When this process is finished, the parameters of global transformation are set to be fixed. Then, the L-BFGS-B opti-mization method is applied to find the optimal parameters of local transformation in nonrigid registration process.
5.4.2 Implementation of Registration Module
As descried in Fig 5.4, the main steps to implement this registration module are:
(i) Import 6 data: the segmented liver, segmented tumor and selected landmarks before LT; the segmented liver, segmented treated region and selected landmarks after LT, (ii) Set the constrained weightλ, (iii) Register the segmented liver before LT on the segmented liver data after LT using our proposed registration method and obtain the optimal transformation, (iv) Automatically transform the treated region to the optimal position using the same transformation parameter.