Non-Rigid Registration Based on Huber Prior : Application to Multiphase Liver CT Volumes Alignment

This study investigates a new non-rigid registration method based on Huber prior. The goal of image registration is to determine a reasonable transformation, such that a transformed version of the first image is similar to the second one. In order to improve the computing efficiency and registration precision, we model the motion of the object by both affine transformation and B-spline based Free Form deformation (FFD). Our innovation is using Huber prior as penalty term of the energy function for better results. The FFD model with Huber prior which has high accuracy and strong robustness can settle the over-smoothing and edge information deficiency problem. We applied the proposed algorithm to contrast-enhanced liver CT volumes registration. Distinctly, the experiment results show that our method gets better results compared to affine transformation and original FFD.


Introduction
Multiphase liver CT exploits the temporal difference in the contrast material delivery to the liver between arterial and portal phases.By obtaining CT volumes at different phases after intravenous bolus injection of contrast material, the hemodynamic characteristics of different hepatic focal lesions, such as hepatocellular carcinoma (HCC) and metastases are displayed (1) .The work of radiologists is searching a set of 2D slice images for the correspondences of different phases and visually comparing them.If there is an effective method for multiphase liver CT registration, identifying a lesion that shows a dynamic enhancement pattern different from that of the background liver (lesion detection), and classifying this enhancement pattern (lesion characterization) will be much easier and more efficient (2) .
However, registration of multiphase liver is a particularly challenging task, because the liver is a highly deformable organ and deformed passively by complex movements of other organs.Traditional registration algorithms, such as affine transformation, Optical Flow and Thin Plate Spline are not accurate enough.In this paper, we introduce Huber function (3,4) as prior term to multiphase liver CT alignment.Although intensity of hepatic regions between phases is different due to contrast agent, we still use SSD (5) as similarity measure due to its less computing time without loss of accuracy.

Free form deformation (FFD)
The goal of image registration in contrast-enhanced liver volumes is to relate any point in the post-contrast enhanced image to the pre-contrast enhanced reference image.Affine transformation only describes the global motion of the liver, which is insufficient.Therefore, we develop a combined transformation which consists of a global transformation and a local transformation (6) .
( , , ) ( ( , , )) T(x, y, z) denotes the deformation of any location (x, y, z).Firstly, we use affine transformation to describe the overall motion of the liver, and then we choose FFD, which is a powerful tool for modeling 3D deformable objects, to model the local deformation.
To define a spline-based FFD, we denote the domain of the image volume as Let Φ denote a , the deformation at these coordinates is computed from the positions of the surrounding 4×4×4 neighborhood of control points (6,7) .
, and l B is the lth basis function of the uniform B-spline (8,9) .

Non-rigid registration based on Huber function
The total objective function for estimating the optimal deformation field is formulated as: ( , ( )) ( ) I and J denote the reference and template image, respectively.λ is the weighting parameter which signifies the tradeoff between the alignment of the two images and the smoothness of the transformation T .We selected SSD as similarity measure for its faster computing speed without loss of accuracy and Huber function as the penalty term.SSD is defined as: where N represents the total voxel number of the volumes.When the two volumes are geometrically aligned, the SSD value reaches its minimum (5) .The traditional penalty term, least-squares (l 2 ) model, is widely known for its drawback to large errors: results in discontinuity points over-smoothing.The l 1 method is robust against large errors (3) .However, when the data contains many small errors, the l 1 solution can be undesirably biased toward a subset of the discontinuity points.This indicates that l 1 is not suitable in general as penalty term.Neither the l 2 nor l 1 model has flexible discriminatory power to recognize and treat differently large errors and small errors.
As Huber function has the advantage of being convex and more or less edge-preserving (10) , we utilize it as prior to solve this problem.Huber function can provide a smoother estimate than l 1 and is more robust to large errors compared with l 1 .
The formula of Huber function and derivative are: Consequently, the penalty term can be rewritten as: The limited memory BFGS method (l-BFGS) is described by Nocedal et al., 1980 (11) .Because the Huber function is not twice continuously differentiable, the Hessian is not computed directly but approximated using a limited Memory BFGS (l-BFGS) to optimize the energy function and get displacement fields (7) .T T (7)   This framework can not only reduce the computation time, but also effectively avoid the local minima (6) .The number of hierarchy levels is mostly dependent on the quality and the original size of the input images.For the experiments presented in this paper, it was sufficient to use a maximum of two levels of hierarchy.

Experiments
We evaluate the effectiveness of the proposed method on two 3D liver data from the same patient scanned at the arterial and the portal phase.The set of CT data scanned at different times were acquired using a Philips Brilliance 64-Slice CT scanner and provided by Tianjin Medical University General Hospital.The size of tested CT data is 512 × 512 × 397, and the voxel size is 0.68 mm × 0.68 mm × 0.50 mm.Considering the computation speed and memory, we reduce the voxel size from 512 × 512 × 397 to 256 × 256 × 132 in the experiment.
Figure 2 shows the 100 th slice in coronal plane, 120 th slice in vertical plane and 80 th slice in horizontal plane of arterial and the portal phase, respectively.
The deformation algorithm started with global affine transformation and then local transformation with control points spacing from 32 (voxel unit) to 16. Threshold γ is set to e -7 and weighting parameter λ is set to 4. Brighter or darker edges can be seen as indicated.
Figure 3 demonstrates the registration results by residual plot.Since the data were scanned at different times during the dynamic contrast enhancement, the residual plot would not be zero even if it's 100% perfectly registered.But we still can tell which result is better by comparing the edges in the plots.If two volumes are perfectly registered, the edges would be exactly aligned, which means there is no too bright or too dark.As indicated by the arrows, we can see obvious brighter or darker edges in group (a), (b) and (c).So it proves that edge information get better retained in our proposed method.
In addition, the mean and variance of the sum of absolute differences (SAD) as well as the correlation coefficient (CC) are calculated to assess the quality of our method.The results are summarized in Table 1.Comparing to conventional FFD, the average and variance of the volumn error was reduced by 6.6% and 22.8%.
To prove our method is robust to noise, we have added Gaussian noise with (0, 0.0001) N to the reference and template volumes.The results shown in Figure .4and Table 2 prove that edge information can still be obtained when adding noise with Gaussian.We select the 90 th slice in horizontal plane and 120 th slice in coronal plane at random.In the proposed algorithm, the average and variance of the volumn error with Gaussian noise was reduced by 4.7% and 20% compared to the conventional FFD.

Conclusions
Our innovation is the introduction of Huber prior to medical image registration.By using it as penalty term, the edge information of the images is better retained.In addition, our method is more robust to noise.In the experimental section, we show our advantages in both quantitative and qualitative analysis.Therefore, we can draw a conclusion that using a hybrid function as penalty term is feasible.

Figure 1
Figure 1 plots Huber function ) (x ρ , least-squares (l 2 ) model 2 ( ) L x and l 1 function 1 ( ) L x in the first row, the first derivative ( ) ' ρ x , 2 ( ) ' L x and 1 ( ) ' L x in the second row.As illustrated, Huber function ensures a smooth transition between l 1 and l 2 at γ .When γ is set to a sufficiently large value, the Huber function approximates the least squares model.Similarly, as γ approaches zero, the Huber function approaches the l 1 function as well.

FFD 1 (
model is parameterized by the coefficients of a set of sparse, uniformly-spaced control points.As the increase of the control points, the more accurate results are obtained, along with more computation time.In order to compromise the accuracy and efficiency, we use the multi-resolution framework in this paper.Let 1,…, L denote a hierarchy of control point meshes at different resolutions.Each level of resolution defines a transformation l T and their sum defines the local transformation T .

Fig. 4 .
Fig. 4. Intensity residual image of noise with Gaussian noise.(a) Initial difference; (b) Affine; (c) Conventional FFD; (d) FFD based on Huber.Brighter or darker edges can be seen as indicated.

Table . 1
: The registration error of sum of absolute differences (SAD) and correlation coefficient (CC) for different types of methods.