Non-rigid MR-TRUS image registration for image-guided prostate biopsy using correlation ratio-based mutual information
© The Author(s) 2017
Received: 10 September 2016
Accepted: 27 December 2016
Published: 10 January 2017
To improve the accuracy of ultrasound-guided biopsy of the prostate, the non-rigid 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 MR-TRUS image registration. However, the use of MI has been challenged because of intensity distortion, noise and down-sampling. Hence, we need to improve the MI measure to get better registration effect.
We present a novel two-dimensional non-rigid MR-TRUS registration algorithm that uses correlation ratio-based 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.
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.
With more accurate matching performance and lower sensitivity to noise and down-sampling, 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.
KeywordsProstate Needle biopsy Non-rigid registration CRMI ISGD
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 . 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 . Transrectal ultrasound-guided (TRUS) prostate biopsy is the most common means of diagnosing prostate cancer when an individual exhibits high blood levels of prostate-specific antigen (PSA). Although TRUS has many advantages, including real-time detection, low cost, and easy operation, its poor image quality and lack of clear contrast between malignant and normal tissue lead to false-negative rates of up to 30% for systematic sextant biopsies . In contrast, magnetic resonance (MR) is the most sensitive imaging modality for observing anatomical structures and locating prostate tumor. Therefore, the registration of pre-operative MR images onto inter-operative TRUS images is of important clinical significance for improving biopsy accuracy.
Pre-operative and inter-operative 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 . To compensate for these movements, many non-rigid registration algorithms have been proposed over the past 20 years [5–8]. Mutual information (MI) is a widely used similarity criterion in multi-modal image registration and has been independently proposed by Collignon  and Viola . Moradi  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 B-splines 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  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 L-BFGS optimization method. However, several recent studies have shown that MI-based 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  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. CMI-based registration resulted in a significant improvement in theoretical, phantom and clinical data compared with MI-based 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 . Another kind of method assigns different weights to pixels using feature detection operators, e.g., the saliency measure  and the Harris corner detector . 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 information-based method termed correlation ratio-based 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  to an stochastic gradient descent (SGD) optimizer, which increases both the random disturbance and the probability of escaping the local extremum.
On the other hand, MI requires enough samples to compute the probability density functions. This limits the application of MI in multi-resolution 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 has been applied to register multi-modal 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 non-linear 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.
Although mapping all intensities between a float image and a reference image using one function may be over-constrained, 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 MI-based registration as prior information to CR-based registration. However, this approach is not advisable because it might degrade the performance into the original CR-based 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 scatter-plot 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 scatter-plot dispersion of intensity values.
For an arbitrary control point u i,j , we derive the derivative of the three terms independently.
Klein  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  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.
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 mm2; the size of the TRUS images is 383 × 691 pixels, and the spatial resolution is 0.1667 × 0.1667 mm2. To improve registration accuracy, the TRUS images are resampled to have the same pixel size as the MR images, and a zero-padding 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 Chan-Vese model. For the prostate images of all 12 patients, the corresponding anatomical landmarks are selected by clinical expert.
We employ the common metrics  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  which is widely used in 3D Slicer and CMI  which exhibits higher registration accuracy over MI in many cases. All three algorithms use the same multi-resolution 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 scatter-plot 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 scatter-plot 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 down-sampling and is more suitable for multi-resolution 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 closed-loop 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 non-rigid 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 real-time registration for TRUS-guided prostate biopsy.
prostate specific antigen
label map registration frame
iterative closest point
target registration error
normalized mutual information
dice similarity coefficient
conditional mutual information
signal to noise ratio
correlation ratio-based mutual information
stochastic gradient descent
adaptive stochastic gradient descent
improved stochastic gradient descent
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.
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.
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.
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.
- Torre LA, Bray F, Siegel RL, Ferlay J, Lortet-Tieulent 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 ultrasound-guided 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 spline-based non-linear 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. Real-time MRI-TRUS 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 image-guided 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 microbubble-mediated ultrasound-MRI 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 ultrasound-guided 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 image-guided 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 B-splines 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 free-form 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. Parzen-window 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 CT-ultrasound registration for diagnostic imaging and image-guided 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 MR-ultrasound registration for image-guided 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 B-splines. 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