Nonrigid MRTRUS image registration for imageguided prostate biopsy using correlation ratiobased mutual information
 Lun Gong^{1, 2},
 Haifeng Wang^{3},
 Chengtao Peng^{4},
 Yakang Dai^{1},
 Min Ding^{5},
 Yinghao Sun^{3},
 Xiaodong Yang^{1} and
 Jian Zheng^{1}Email author
DOI: 10.1186/s1293801603085
© The Author(s) 2017
Received: 10 September 2016
Accepted: 27 December 2016
Published: 10 January 2017
Abstract
Background
To improve the accuracy of ultrasoundguided biopsy of the prostate, the nonrigid registration of magnetic resonance (MR) images onto transrectal ultrasound (TRUS) images has gained increasing attention. Mutual information (MI) is a widely used similarity criterion in MRTRUS image registration. However, the use of MI has been challenged because of intensity distortion, noise and downsampling. Hence, we need to improve the MI measure to get better registration effect.
Methods
We present a novel twodimensional nonrigid MRTRUS registration algorithm that uses correlation ratiobased mutual information (CRMI) as the similarity criterion. CRMI includes a functional mapping of intensity values on the basis of a generalized version of intensity class correspondence. We also analytically acquire the derivative of CRMI with respect to deformation parameters. Furthermore, we propose an improved stochastic gradient descent (ISGD) optimization method based on the Metropolis acceptance criteria to improve the global optimization ability and decrease the registration time.
Results
The performance of the proposed method is tested on synthetic images and 12 pairs of clinical prostate TRUS and MR images. By comparing label map registration frame (LMRF) and conditional mutual information (CMI), the proposed algorithm has a significant improvement in the average values of Hausdorff distance and target registration error. Although the average Dice Similarity coefficient is not significantly better than CMI, it still has a crucial increase over LMRF. The average computation time consumed by the proposed method is similar to LMRF, which is 16 times less than CMI.
Conclusion
With more accurate matching performance and lower sensitivity to noise and downsampling, the proposed algorithm of minimizing CRMI by ISGD is more robust and has the potential for use in aligning TRUS and MR images for needle biopsy.
Keywords
Prostate Needle biopsy Nonrigid registration CRMI ISGDBackground
Prostate cancer is among the most common diseases in males and exhibits a large and increasing morbidity in many countries. As the fifth leading cause of cancer death worldwide, prostate cancer affected approximately 1.1 million men worldwide in 2012 [1]. In China, for example, according to statistics, approximately 49,000 new cases of prostate cancer were reported, ranking ninth in terms of cancer incidence in men in 2011 [2]. Transrectal ultrasoundguided (TRUS) prostate biopsy is the most common means of diagnosing prostate cancer when an individual exhibits high blood levels of prostatespecific antigen (PSA). Although TRUS has many advantages, including realtime detection, low cost, and easy operation, its poor image quality and lack of clear contrast between malignant and normal tissue lead to falsenegative rates of up to 30% for systematic sextant biopsies [3]. In contrast, magnetic resonance (MR) is the most sensitive imaging modality for observing anatomical structures and locating prostate tumor. Therefore, the registration of preoperative MR images onto interoperative TRUS images is of important clinical significance for improving biopsy accuracy.
Preoperative and interoperative prostate shape may suffer deformation due to extrusion of the ultrasonic probe, inflation of the anorectal coil inside the rectum during the MR scanning, and alterations in patient position [4]. To compensate for these movements, many nonrigid registration algorithms have been proposed over the past 20 years [5–8]. Mutual information (MI) is a widely used similarity criterion in multimodal image registration and has been independently proposed by Collignon [9] and Viola [10]. Moradi [11] created a label map registration frame (LMRF) that aligned TRUS and MR images by using 3D Slicer, which first used the iterative closest point (ICP) method to rigidly align the outlines of the two images and then combined MI and Bsplines to elastically register binary label maps from the two images obtained by manual contouring. Although LMRF aligned contours with high accuracy and could be conveniently implemented in 3D Slicer, it used only label maps to correct the local deformation and ignored pixel intensity, resulting in a high target registration error (TRE) of 3.6 ± 1.7 mm. Mitra [12] utilized directional quadrature filter pairs to convert TRUS and MR images into texture images and then used normalized mutual information (NMI) as a similarity criterion to register the texture images. The average TRE and mean Dice Similarity coefficient (DSC) were 2.64 ± 1.47 mm and 0.943 ± 0.039, but the computation time exceeded 797 s due to the use of 4 quadrature convolutions and the LBFGS optimization method. However, several recent studies have shown that MIbased registration can be improved in certain cases. One improvement is to calculate MI over a set of overlapping image blocks to include spatial information. Loeckx [13] proposed the conditional mutual information (CMI), which considered the spatial location in the reference image of each joint intensity pair as the priori condition and calculated conditional entropies between the intensities given the spatial distribution. CMIbased registration resulted in a significant improvement in theoretical, phantom and clinical data compared with MIbased registration, but it was approximately 15 times slower than MI because it regarded the block label as the third channel of the joint histogram. Furthermore, MI assumes that all pixels in the overlapping area affect the calculations equally, but it is clear that different pixels contribute differently to the computation of MI [14]. Another kind of method assigns different weights to pixels using feature detection operators, e.g., the saliency measure [14] and the Harris corner detector [15]. But this method is often unable to extract effective features from TRUS images due to the low signal to noise ratio (SNR).
In this work, we design a mutual informationbased method termed correlation ratiobased mutual information (CRMI) that includes the functional dependence of intensity values. MI only includes the correspondence of intensity classes to correct for the deformation of location but ignores the possible relationship between intensity values. We also deduce the analytical derivative of CRMI in detail. Furthermore, to improve the global searching ability and reduce the registration time, we introduce Metropolis acceptance criteria [16] to an stochastic gradient descent (SGD) optimizer, which increases both the random disturbance and the probability of escaping the local extremum.
Methods
MI
On the other hand, MI requires enough samples to compute the probability density functions. This limits the application of MI in multiresolution cases because the error between the estimation of the joint density and the real density will increase at low resolution. Furthermore, MI establishes a purely qualitative statistical relationship between the intensity classes, which means that the joint density is sensitive to noise.
CR
CR has been applied to register multimodal images in [20–22]. Assuming a functional relationship between the intensities of float images and reference images is not an insurmountable obstacle for rigid or affine registration but is too restrictive to conform to the assumption for nonlinear registration, especially when the image pairs have completely different presentations for the same anatomical structure. Therefore, CR is only applicable to rigid registration in most cases.
CRMI
Although mapping all intensities between a float image and a reference image using one function may be overconstrained, it is reasonable to assume that the intensities of the aligned pixel pairs can be mapped by the same function. It seems feasible to combine the functional dependence of the intensity values with MI by utilizing the deformation parameters of MIbased registration as prior information to CRbased registration. However, this approach is not advisable because it might degrade the performance into the original CRbased method if the alignment obtained by MI is not reasonably accurate. Furthermore, the registration time is excessive due to the need to separately minimize MI and CR. Therefore, we extend the location alignment to include the functional mapping between the intensity pairs by multiplying MI with CR. The scatterplot dispersion of intensity values is measured by combining the joint density of the intensity classes with the mapping relationship of the intensity values. We not only correct the location deformation, but also decrease the influence of intensity bias on the joint density.
Figure 2d shows the translational results of CRMI with respect to the horizontal shift. It can be observed that the registration function of CRMI has the only one global minimum that agrees with the correct transformation parameters. This finding demonstrates that CRMI is more reliable for measuring the scatterplot dispersion of intensity values.
Gradient
For an arbitrary control point u _{ i,j }, we derive the derivative of the three terms independently.
Optimization
Klein [27] automatically calculated these parameters based on the distribution of the voxel displacements using the adaptive stochastic gradient descent method (ASGD); however, this method required computation of the Jacobin matrix of the deformation model with respect to the transformation parameters, which would be unacceptable if the number of transformation parameters was large. Qiao [28] proposed fast ASGD on the basis of ASGD to estimate the three parameters. This approach simplified the calculation of the bias and variance of the voxel displacements by using frequency statistics, a method that improved the efficiency at the cost of some precision. In addition, SGD requires a sufficient number of iterations to achieve convergence. The maximum number of iterations is also difficult to select. If the number is large, computation time will be increased; otherwise, convergence might not be achieved. Therefore, it is important to achieve the optimal efficiency and precision by setting up an appropriate stop condition.
Experiments and results
In this section, we validate the effective performance of the proposed method using synthetic data and clinical prostate data. The proposed method is implemented in C++ and run on a PC equipped with an Intel Core i7 3.60 GHz CPU and 8 GB RAM.
Synthetic data
Clinical prostate data
The proposed method is validated on prostate TRUS and MR images representing 12 patients, which are obtained from Shanghai Changhai Hospital. The size of the MR images is 384 × 384 pixels, and the spatial resolution is 0.4688 × 0.4688 mm^{2}; the size of the TRUS images is 383 × 691 pixels, and the spatial resolution is 0.1667 × 0.1667 mm^{2}. To improve registration accuracy, the TRUS images are resampled to have the same pixel size as the MR images, and a zeropadding strategy is applied to the resampled TRUS images such that the images are the same size with the MR images. The prostate MR and TRUS images are segmented using the ChanVese model. For the prostate images of all 12 patients, the corresponding anatomical landmarks are selected by clinical expert.
We employ the common metrics [12] such as DSC, HD, TRE and registration time to evaluate the registration performance. DSC and HD measure the global registration accuracy, and TRE describes how well the local details are aligned. A high value of DSC means that the prostate region on a moving image maps well to that on a fixed image. A low value of HD indicates good correspondence between the coordinates of two outlines. The lower the value of TRE is, the better the local registration accuracy has been achieved.
Then we compare the proposed method against two other algorithms: LMRF [11] which is widely used in 3D Slicer and CMI [13] which exhibits higher registration accuracy over MI in many cases. All three algorithms use the same multiresolution framework to improve search efficiency, and the number of bins is set to 64.
As shown in Fig. 6, the average TRE values are 3.78 ± 1.47, 3.19 ± 1.28, and 2.58 ± 1.08 mm for LMRF, CMI and CRMI, respectively, showing that the local homologous structures that are optimized using CRMI have the highest alignment accuracy. The average HD value of CRMI is also significantly reduced. Although the average DSC value of CRMI is not significantly better than that of CMI, it is increased compared with that of LMRF. The improvement in these global measures validates that the new similarity criterion aligns the global deformation more reliably at low resolution. The registration time consumed by CRMI is approximately 16 times less than that consumed by CMI and is similar to LMRF, although CRMI requires the extra calculation of CR.
For most samples, the three metrics of DSC, HD and TRE of CRMI are the best. For sample No. 1, it is easy to observe that the correspondence between the two outlines aligned using CRMI is the optimal (see No. 1 in Fig. 8), and the distances between the landmarks are the least (see No. 1 in Fig. 9). For sample No. 12, it is clear that the registration accuracy of CRMI is higher than CMI, whereas LMRF fails to correct the displacements, thus further confirming that the new similarity criterion is better at evaluating the scatterplot dispersion of the intensities in the presence of noise and intensity distortion.
However, some cases remain that should be considered. For sample No. 3, the global measures obtained using CRMI are the best, but the TRE value is not the least. As the partial enlarged drawing of No. 3 shown in Fig. 9, the distances between the landmarks obtained using CRMI are larger than LMRF. The reason might be that CRMI tends to search more accurate global movement parameters at low resolution because of the sample robustness; however, it falls into a local extremum at high resolution due to the inappropriate step factor sk. For samples No. 11, the checkerboard image obtained using CRMI is more approximate, and the lower left of the fusion image exhibits a smaller overlap rate than CMI, corresponding to the bar graphs shown for DSC and HD. It occurs because the parameter T _{0} decides whether ISGD receives the current solution or not is too large; therefore, the current result will be received even if the cost function changes a lot, eventually leading to premature convergence at low resolution. Unexpectedly, as the checkerboard images of sample No. 1 shown in Fig. 8, the checkerboard image obtained using CMI is smoother than LMRF, while the HD values are opposite. This occurs because the maximum distance between the outlines is not observable when the number of checkerboards is inappropriate.
The sample No. 7 is an extreme case that the values of DSC, HD and TRE obtained using CRMI are far worse than their average values. Considering that the other two algorithms also fail to align this group of images, the reason might be that the images are seriously affected by noise and intensity distortion, so that the intensities could not satisfy the assumption of intensity class correspondence or the functional mapping between aligned pixel pairs.
Discussion and conclusion
In this study, we propose and verify a new similarity metric termed CRMI for 2D MR and TRUS prostate alignment. This method corrects the scatterplot dispersion of the intensities caused by the deformation of location and intensity bias. Obtaining knowledge about the functional mapping between the intensities of float images and reference images is crucial when searching for the optimal alignment. We incorporate the functional dependence of the intensities quantified by CR into the intensity class correspondence held by MI. Experiments based on simulated data demonstrate that CRMI is more smooth and reliable for aligning images with intensity distortion. The new metric is also robust to noise due to its calculation of the functional correlation in L ^{2} space, which indicates that the intensity values of noise will be replaced by those of their neighboring samples because we map the pixel pairs with the smallest distance in intensity space. Furthermore, the computational error between the estimation of the joint density and the real density caused by insufficient sampling will be corrected during the searching of intensity mapping, which ensures that CRMI is not sensitive to downsampling and is more suitable for multiresolution optimization strategy.
Another innovation presented here is the improvement of the global searching ability of SGD with less computation time. We use the exact value of the cost function within a closedloop system to select a suitable step size for each iteration as well as the stop criteria. Therefore, the improved SGD optimizer can adaptively determine the optimal number of iterations for different data, which will avoid hemstitching phenomenon due to the use of redundant iterations. An important parameter to determine the best step is the step factor sk. At low resolution, this value is set to 0.7 to save time. At medium and high resolution, this value is set to 0.17 to capture a more accurate step size. While searching for the minimum of the cost function, we not only accept smaller values but also consider relatively large values at a certain probability level by using Metropolis acceptance criteria. Therefore, the improved SGD optimizer is more likely to jump out of the local extremum because of the additional random disturbance. The constant T _{0} directly affects the global searching ability of ISGD. Considering that the exact cost function value changes within a value of 0.001, this constant is set to 0.0006, 0.00054 and 0.0008 for low, medium and high resolution, as experimentally determined. Comparative experiments using clinical prostate data verify that ISGD yields more accurate optimization performance and requires less registration time.
After comparing two nonrigid registration algorithms of LMRF and CMI, the proposed method of minimizing CRMI using ISGD provides the largest average overlap rate and the least mean HD value. The alignment accuracy of the landmarks also has a significant improvement. The average registration time is similar to that using LMRF, which is approximately 16 times less than that for CMI. These improvements confirm that the proposed method is more suitable for aligning 2D MR and TRUS images. The main limitation is that a few special cases may not fulfil the assumption of the functional mapping between aligned pixel pairs. The proposed algorithm may be extended to 3D MR and TRUS registration which can provide a wide range of information about the prostate, and the parallel computing capability of GPU can be exploited to realize realtime registration for TRUSguided prostate biopsy.
Abbreviations
 TRUS:

transrectal ultrasound
 PSA:

prostate specific antigen
 MR:

magnetic resonance
 MI:

mutual information
 LMRF:

label map registration frame
 ICP:

iterative closest point
 TRE:

target registration error
 NMI:

normalized mutual information
 DSC:

dice similarity coefficient
 LBFGS:

LimitedBroydenFletcherGoldfarbShanno
 CMI:

conditional mutual information
 SNR:

signal to noise ratio
 CRMI:

correlation ratiobased mutual information
 SGD:

stochastic gradient descent
 CR:

correlation ratio
 ASGD:

adaptive stochastic gradient descent
 ISGD:

improved stochastic gradient descent
 HD:

Hausdorff distance
Declarations
Authors’ contributions
Study concept and design LG, CP and MD; drafting of the manuscript LG; critical revision of the manuscript for important intellectual content LG and JZ; obtained funding JZ and HW; administrative, technical, and material support LG, JZ, HW and YS; selecting the anatomical landmarks HW; study supervision YD and XY. All authors read and approved the final manuscript.
Acknowledgements
The authors would like to thank colleagues Dr. Cheng Zhang and Dr. Ming Li, Shanghai Changhai Hospital for providing the clinical data. The authors are also grateful to anonymous reviewers for their valuable comments and suggestions.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
The datasets during the current study are available from the corresponding author on the reasonable request.
Consent for publication
The data for experiments has obtained the consent of the patients and Shanghai Changhai hospital, and has been allowed for publication.
Ethics approval and consent to participate
All data has been approved by the Medical Ethics Committee of Suzhou Institute of Biomedical Engineering and Technology, Chinese Academy of Science, and has been allowed for carrying out experiments.
Funding
This work is supported in part by the National Key Research and Development Program of China (Nos. 2016YFC0104500, 2016YFC0104505), the Youth Innovation Promotion Association CAS (No. 2014281), National Natural Science Foundation of China (No. 61201117), the Science and Technology Program of Suzhou (No. ZXY2013001), the Shanghai Health and Family Planning Commission research projects (No. 201540812), the Shanghai Changhai Hospital New Medical Technology Breeding Program (No. NT201507).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Torre LA, Bray F, Siegel RL, Ferlay J, LortetTieulent J, Jemal A. Global cancer statistics, 2012. CA Cancer J Clin. 2015;65(2):87–108.View ArticleGoogle Scholar
 Hao J, Zhao P, Chen WQ. Chinese cancer registration report 2015. Beijing: Press of Military Medical Sciences; 2015.Google Scholar
 Norberg M, Egevad L, Holmberg L, Sparen P, Norlen B, Busch C. The sextant protocol for ultrasoundguided core biopsies of the prostate underestimates the presence of cancer. Urology. 1997;50(4):562–6.View ArticleGoogle Scholar
 Mitra J, Kato Z, Marti R, Oliver A, Llado X, Sidibe D, Ghose S, Vilanova JC, Comet J, Meriaudeau F. A splinebased nonlinear diffeomorphism for multimodal prostate registration. Med Image Anal. 2012;16(6):1259–79.View ArticleGoogle Scholar
 Xu S, Kruecker J, Turkbey B, Glossop N, Singh A, Choyke P, Pinto P, Wood B. Realtime MRITRUS fusion for guidance of targeted prostate biopsies. Comput Aided Surg. 2008;13(5):255–64.View ArticleGoogle Scholar
 Hu Y, Ahmed HU, Taylor Z, Allen C, Emberton M, Hawkes D, Barratt D. MR to ultrasound registration for imageguided prostate interventions. Med Image Anal. 2012;16(3):687–703.View ArticleGoogle Scholar
 Hou M, Chen C, Tang D, Luo S, Yang F, Gu N. Magnetic microbubblemediated ultrasoundMRI registration based on robust optical flow model. Biomed Eng Online. 2015;14(Suppl 1):S14.View ArticleGoogle Scholar
 Natarajan S, Marks LS, Margolis DJ, Huang J, Macairan ML, Lieu P, Fenster A. Clinical application of a 3D ultrasoundguided prostate biopsy system. Urol Oncol. 2011;29(3):334–42.View ArticleGoogle Scholar
 Collignon A, Maes F, Delaere D, Vandermeulen D, Suetens P, Marchal G. Automated mutimodality image registration using information theory. Inf Process Med Imaging. 1995;3:263–74.Google Scholar
 Viola P, Wells WM. Alignment by maximization of mutual information. 5th International Conference on Computer Vision 1995. p. 16–23.
 Moradi M, Janoos F, Fedorov A, Risholm P, Kapur T. Two solutions for registration of ultrasound to MRI for imageguided prostate interventions. International Conference of the IEEE Engineering in Medicine and Biology Society 2012. p. 1129–32.
 Mitra J, Marti R, Oliver A, Llado X, Ghose S, Vilanova JC, Meriaudeau F. Prostate multimodality image registration based on Bsplines and quadrature local energy. Int J Comput Assist Radiol Surg. 2012;7(3):445–54.View ArticleGoogle Scholar
 Loeckx D, Slagmolen P, Maes F, Vandermeulen D, Suetens P. Nonrigid image registration using conditional mutual information. IEEE Trans Med Imaging. 2010;29(1):19–29.View ArticleGoogle Scholar
 Luan HX, Qi FH, Xue Z, Chen LY, Shen DG. Multimodality image registration by maximization of quantitative–qualitative measure of mutual information. Pattern Recogn. 2008;41:285–98.MATHView ArticleGoogle Scholar
 Woo J, Stone M, Prince JL. Multimodal registration via mutual information incorporating geometric and spatial context. IEEE Trans Image Process. 2015;24(2):757–69.View ArticleGoogle Scholar
 Kirkpatrick S, Gelatt C, Vecchi M. Optimization by simulated annealing. Science. 1983;220(4598):671–80.MathSciNetMATHView ArticleGoogle Scholar
 Rueckert D, Sonoda L, Hayes C, Hill D, Leach M, Hawkes D. Nonrigid registration using freeform deformations: application to breast MR images. IEEE Trans Med Imaging. 1999;18(8):712–21.View ArticleGoogle Scholar
 Thevenaz P, Unser M. Optimization of mutual information for multiresolution image registration. IEEE Trans Med Imaging. 2000;9(12):2083–99.MATHGoogle Scholar
 Xu R, Chen YW, Tang SY, Morikawa S, Kurumi Y. Parzenwindow based normalized mutual information for medical image registration. IEICE Trans Inf Systems. 2008;91(1):132–44.View ArticleGoogle Scholar
 Alexis R, Grégoire M, Xavier P, Nicholas A. Multimodal image registration by maximization of the correlation ratio. Institut National De Recherche En Informatique Et En Automatique; 1998.
 Wein W, Brunke S, Khamene A, Callstrom MR, Navab N. Automatic CTultrasound registration for diagnostic imaging and imageguided intervention. Med Image Anal. 2008;12(5):577–85.View ArticleGoogle Scholar
 Milko S, Melvaer EL, Samset E, Kadir T. Evaluation of bivariate correlation ratio similarity metric for rigid registration of US/MR images of the liver. Int J Comput Assist Radiol Surg. 2009;4(2):147–55.View ArticleGoogle Scholar
 Pluim JPW, Maintz JBA, Viergever MA. Image registration by maximization of combined mutual information and gradient information. IEEE Trans Image Process. 2000;19(8):809–14.View ArticleGoogle Scholar
 Rivaz H, Chen S, Collins D. Automatic deformable MRultrasound registration for imageguided neurosurgery. IEEE Trans Med Imaging. 2015;34(2):365–80.View ArticleGoogle Scholar
 Klein S, Staring M, Pluim JPW. Evaluation of optimization methods for nonrigid medical image registration using mutual information and Bsplines. IEEE Trans Image Process. 2007;16(12):2879–90.MathSciNetView ArticleGoogle Scholar
 Klein S, Staring M, Pluim JPW. Comparison of gradient approximation techniques for optimisation of mutual information in nonrigid registration. Medical Imaging: Image Process; 2005.View ArticleGoogle Scholar
 Klein S, Pluim JPW, Staring M, Viergever MA. Adaptive stochastic gradient descent optimisation for image registration. Int J Comput Vision. 2009;81(3):227–39.View ArticleGoogle Scholar
 Qiao Y, Lew BV, Lelieveldt BPF, Staring M. Fast automatic step size estimation for gradient descent optimization of image registration. IEEE Trans Med Imaging. 2016;35(2):391–403.View ArticleGoogle Scholar