Rigid Registration for Image-Guided Interventions

1. Rigid Registration for Image-Guided Interventions Ziv YanivComputer Aided Interventions and Medical Robotics, Georgetown University Last Modified September 2008 (Image-guided interventions Tutorial MICCAI’08)

2. Why? • Rigid registration is a basic human activity:

3. Tumor location is barely visible on CT, but clearly visible in fused PET/CT. Why? (in our context) • Combine multiple information sources to improve interventional outcome (e.g. PET/CT). • Alignment of virtual/image space and the physical patient space. • An initial step for non-rigid registration.

4. Sources of Data • Images: • Diagnostic (primarily 3D): anatomical (CT,MR…), functional (PET, SPECT…) • Interventional (primarily 2D): projective (X-ray, video…), tomographic (US, CT-Fluorocsopy) • Digitized points: • 3D digitized points: contact-based digitization, laser range scanning, stereo imaging.

5. 3D/3D Registration

6. Definition Given two data sets defined over the domains with, possibly, additional side information (e.g. pairing between corresponding features) find the rigid transformation: where Specific instances: are the identity map, we have point set to point set registration. map locations to intensity values, we have image to image registration.

7. Paired-Point Methods are the identity map, and point correspondences are known. • Analytic least squares solutions exist (e.g. Arun et al. 1987, Horn 1987), and are optimal if it is assumed that the measurement process is corrupted by additive isotropic i.i.d. zero mean Gaussian noise. • Correspondence is assumed to be correct (breakdown point of 1). • Require that points be detected in both coordinate systems – anatomical landmarks or fiducials. • Fiducials: • accurate localization by design • number of points is independent of the anatomical structures.

8. IGI - Assumptions can Kill • What happens if the data acquisition process doesn’t match our assumptions? • An acquisition which is biased by a rigid transformation can have serious implications as the algorithm cannot detect this situation and will compensate for the bias:

9. Head and Hat • Align two surface representations of the brain segmented from complementary modalities (Pelizzari et al. 1989). • Head – higher resolution data, represented as stacked 2D closed contours. • Hat – lower resolution data, represented as a set of points. • Minimize (Powell’s direction set method) the sum of distances between the head surface and the rays emanating from the head centroid and going through the hat points. • Point matching strategy suited for spherical objects.

10. Head and Hat hat head

11. Iterative Closest/Corresponding Point (ICP) A natural extension of the analytic paired-point algorithms (Besl and McKay 1992, Zhang 1994, Chen and Medioni 1992). • Input: point set P, surface S t>0 - improvement threshold n - maximal number of iterations • Initialization: • Set cumulative transformation, and apply to points. • Find corresponding points and compute similarity (e.g. root mean square distance). • Iterate: • Compute incremental transformation using the current correspondences (e.g. analytic least squares solution), update cumulative transformation and apply to points. • Find corresponding points and compute similarity. If improvement in similarity is less than t, or number of iterations has reached n terminate.

12. ICP – It almost works • Convergence is local, requiring an initial transformation near the optimum. • Point pairing is a computationally expensive operation. • The use of an analytic least squares method for computing the incremental transformations implies that ICP assumes additive isotropic i.i.d. zero mean Gaussian noise, and that there are no outliers.

13. ICP – It almost works • Convergence is local, requiring an initial transformation near the optimum. Improve probability of convergence by: • gross manual alignment of the data • initial gross localization of anatomical landmarks (Ma and Ellis 2003) • use multiple starting points in parameter space (Besl and McKay 1992) • simulated annealing (Penney et al. 2001) • Point pairing is a computationally expensive operation. • The use of an analytic least squares method for computing the incremental transformations implies that ICP assumes additive isotropic i.i.d. zero mean Gaussian noise, and that there are no outliers.

14. ICP – It almost works • Convergence is local, requiring an initial transformation near the optimum. • Point pairing is a computationally expensive operation. Accelerate run time by: • kD spatial data structure (Besl and McKay 1992) • Cache closest points (Simon 1996) • Approximate nearest neighbor searching (Greenspan and Yurik 2003) • Parallel nearest neighbor searching (Langis et al. 2001). • Hierarchical data sampling (Jost and Hugli 2003) • The use of an analytic least squares method for computing the incremental transformations implies that ICP assumes additive isotropic i.i.d. zero mean Gaussian noise, and that there are no outliers.

15. ICP – It almost works • Convergence is local, requiring an initial transformation near the optimum. • Point pairing is a computationally expensive operation. • The use of an analytic least squares method for computing the incremental transformations implies that ICP assumes additive isotropic i.i.d. zero mean Gaussian noise, and that there are no outliers. Replace analytic least squares methods with: • M-estimators (Ma and Ellis 2003, Kaneko et al. 2003) • Least median of squares (Trucco et al. 1999, Masuda and Yokoya 1995) • Weighted least squares (Turk and Levoy 1994, Maurer et al. 1998) • Least trimmed squares (Chetverikov et al. 2005) • Use a fraction of the data (Phillips et al. 2007)

16. Standard Optimization • Replace the ICP match-optimize approach with standard optimization (Fitzgibbon 2003). • Implicitly represent point set/surface using a precomputed uniform distance map. • Minimize using Levenberg-Marquardt with M-estimators instead of sum of squares.

17. Intensity based Registration • Cast registration as an optimization task, requiring three components: • Similarity measure (the objective function). • Optimization algorithm. • Interpolation scheme (intensity values are only given at a discrete set of locations).

18. Popular Similarity Measures

19. Which Similarity Measure • Computational complexity. • Robustness. • Capture range. • Accuracy

20. Which Similarity Measure • Computational complexity. • Theoretical analysis. • Optimize implementation as the similarity measure will be evaluated many times (improve the constant factor). • Robustness. • Theoretical analysis (e.g. SAD more robust than SSD) • Capture range. • Accuracy

21. The blind men and the elephant: Which Similarity Measure • Computational complexity. • Robustness. • Capture range. • Assess registration performance as a function of the similarity measure – evaluates the combination of the similarity measure and the optimization algorithm. • Accuracy • Explore the behavior of the similarity measure as a function of the transformation parameters – axis aligned orthogonal slices through the parameter space

22. Which Similarity Measure • Computational complexity. • Robustness. • Capture range. • Accuracy • Exhaustively sample the parameter space along diameters of a 6D hypersphere centered on the gold standard parameters (Škerl et al. 2006). • Analyze the function behavior based on these observations (e.g. distinctiveness of optimum).

23. Optimization • General purpose optimization methods. • No single algorithm is optimal for all similarity measures (no-free-lunch theorem). • Improve convergence range using heuristic search algorithms (e.g. simulated annealing) – double edged sword. Sometimes the desired pose is a local optimum and not a global one.

24. Interpolation • Linear interpolation – a compromise between accuracy and computational complexity. • May have adverse effects on the similarity measure (Pluim et al. 2000). • B-spline interpolation - best tradeoff between accuracy and computational complexity (Lehmann et al. 1999, Meijering et al. 2001).

25. 2D/3D Registration

26. Definition Given two data sets defined over the domains and camera matrices, , where is a rigid transformation and a perspective projection matrix, find the rigid transformation: where • Closely related to the pose estimation task studied in computer vision. • A single image of a known scene is acquired and we seek the camera’s pose relative to the scene. • Use point and line features which arise in man made environments.

27. Camera Model • X-ray fluoroscopy and endoscopic video are modeled using the standard perspective camera model. • Images are distorted for X-ray fluoroscopy utilizing image intensifiers, and for • standard endoscopes.

28. Iterative Algorithms - Initialization • Manual initialization via visual inspection User interactively explores the parameters space • Clinical setup Intersect bounding pyramids from multiple viewing angles. • Approximate paired-point registration • Localize anatomical landmarks • The same anatomical landmarks (>3) in multiple projection images. • Digitize with a tracked pointer if already exposed. • Localize skin adhesive fiducials.

29. 2D/3D Head and Hat • A variation of the head and hat is used for patient alignment in radiation therapy, X-ray/CT, (Murphy 1997). • A digitally reconstructed radiograph (DRR), simulated X-ray is generated and the optimization is performed in 2D. • The X-ray serves as the head and the DRR as the hat. • Edges (skull) are detected along rays emanating from the center of the X-ray image (head centroid) and going towards the image boundary – both in DRR and X-ray. • The distance between edge points on the rays is minimized, along with the mean grey level value along the ray.

30. 2D/3D ICP • Minimize the distance between back-projected rays and a 3D model, formulation is in the context of computer vision (Wunsch and Hirzinger 1996). • In the medical context this approach requires extracting a 3D surface model and detecting 2D edge points that were generated by the anatomical structure. back-projected rays camera

31. 2D/3D Standard Optimization • Minimize the sum of squared distances between a surface model and back-projected rays, CT/X-ray, (Lavallée and Szeliski 1995). • Surface is implicitly represented using an octree based distance map (rediscovered in computer graphics, Frisken et al. 2000, Adaptively Sampled Distance Fields). • Minimization using Levenberg-Marquardt.

32. Intensity based Registration • Most often X-ray/CT, in some cases X-ray/MR via mapping MR values to CT, and video/CT. • Simulate the 2D image generation process. • Maximize the similarity between the simulated and actual 2D images – choice depends on the realism of the simulation. • Primary challenge - fast generation of simulated images, either X-ray or surface/volume rendering in the case of video.

33. ray casting with mask fluoroscopy Generation of Simulated Images • Use established volume rendering techniques: • Ray casting. • Splatting. • Shear warp. • Image based rendering - Use a set of images to render a new image from a novel camera pose. • In our context, simulate a set of images from specific camera poses and use them to generate all other simulated images (LaRose 2001, Knaan and Joskowicz 2003, Russakoff 2005). • GPU based rendering (Khamene 2006, LaRose 2001). • Standard rendering of polygonal models (Mori et al 2002, Turgeon et al. 2005).

34. Other Approaches • Gradient based registration – directly relate the volumetric gradients to the X-ray image gradients (Livyatan et. al 2003, Tomaževič et al. 2003) : • Perform cone-beam reconstruction from a few X-ray images (2-15) and then we are back to 3D/3D registration (Tomaževič et al. 2006) [Quality of reconstruction is poor, required the introduction of a novel similarity measure].

35. Evaluation Ideal registration algorithm: • Fast: The result is obtained in real time (<0.1sec). • Accurate: The distance between corresponding points after applying the estimated transformation is less than 0.1mm. • Robust: More than half of the data elements must be outliers in order to throw the registration out of reasonable bounds (breakdown point of N/2). • Automatic: No user interaction required. • Reliable: Given the expected input the registration always succeeds.

36. Evaluation • Empirically evaluate using standard data sets with known “gold standard” registration (preferably in vivo): • Retrospective Image Registration Evaluation (in-vivo head CT/MR, PET/MR) [Fitzpatrick]: http://www.insight-journal.org/rire/ • 2D/3D Standardized Evaluation Methodology (in-vitro spine X-ray/MR/CT/CBCT) [Van De Kraats]: http://www.isi.uu.nl/Research/Databases/GS/ • Gold standard data for evaluation of 2D-3D registration (in-vitro spine X-ray/MR/CT) [Tomaževič]: http://lit.fe.uni-lj.si/tools.php?lang=eng

37. Accuracy • Most important criterion. Directly addresses the goal of registration, alignment of corresponding points. • Spatially variant. • When it is not possible to digitize the target to report actual target registration error, report the expected target registration error (e.g. first order approximation [Fitzpatrick et al. 1998], second order approximation [Moghari and Abolmaesumi 2006]).

38. Take Home Message • Be explicit about your assumptions. • Evaluate using all criteria from the definition of the ideal algorithm. • Usability is context dependent, registration is not a goal but a means to an end:

39. Bibliography • Arun KS, Huang TS, and Blostein SD. (1987). “Least-squares fitting of two 3D point sets.” IEEE Transact Pattern Anal Machine Intell, 9(5), 698–700. • Horn BKP. (1987). “Closed-form solution of absolute orientation using unit quaternions.” J Opt Soc Am A, 4(4), 629–642. • Pelizzari CA, Chen GT, Spelbring DR, Weichselbaum RR, and Chen CT. (1989). “Accurate three-dimensional registration of CT, PET, and/or MR images of the brain.” J Comput Assist Tomogr, 13(1), 20–26. • Besl PJ and McKay ND. (1992). “A method for registration of 3D shapes.” IEEE Transact Pattern Anal Machine Intell, 14(2), 239-255. • Zhang Z. (1994). “Iterative point matching for registration of free-form curves and surfaces.” Int J Comput Vision, 13(2), 119–152. • Chen Y and Medioni G. (1992). “Object modeling by registration of multiple range images.” Image Vis Comput, 10(3), 145–155. • Ma B, and Ellis RE. (2003). "Robust registration for computer-integrated orthopedic surgery: laboratory validation and clinical experience." Med Image Anal, 7(3), 237-250. • Penney GP, Edwards PJ, JKing AP, Blackall JM, Batchelor PG, and Hawkes DJ. (2001). "A stochastic iterative closest point algorithm (stochastICP)." Medical Image Computing and Computer-Assisted Intervention (MICCAI 2001), Springer-Verlag, 762-769. • Simon D. (1996). "Fast and Accurate Shape-Based Registration," PhD, Carnegie Mellon University. • Greenspan M, and Yurick M. (2003). "Approximate K-D tree search for efficient ICP." Int. Conference on 3-D Digital Imaging and Modeling, 442-448. • Langis C, Greenspan M, and Godin G. (2001). "The parallel iterative closest point algorithm." International Conference on 3-D Digital Imaging and Modeling, 195-204. • Jost T, and Hugli H. (2003). "A multi-resolution ICP with heuristic closest point search for fast and robust 3D registration of range images." International Conference on 3-D Digital Imaging and Modeling, 427-433.

40. Bibliography • Kaneko S, Kondo T, and Miyamoto A. (2003). "Robust matching of 3D contours using iterative closest point algorithm improved by M-estimation." Pattern Recognition, 36(9), 2041-2047. • Trucco E, Fusiello A, and Roberto V. (1999). "Robust motion and correspondence of noisy 3-D point sets with missing data." Pattern Recognition Letters, 20(9), 889-898. • Masuda T, and Yokoya N. (1995). "A robust method for registration and segmentation of multiple range images." Computer Vision and Image Understanding, 61(3), 295-307. • Turk G, and Levoy M. (1994). "Zippered polygon meshes from range images." Computer graphics and interactive techniques (SIGGRAPH 1994), 311-318. • Maurer CR, Jr., Maciunas RJ, and Fitzpatrick JM. (1998). "Registration of head CT images to physical space using a weighted combination of points and surfaces." IEEE Trans Med Imaging, 17(5), 753-761. • Chetverikov D, Stepanov D, and Krsek P. (2005). "Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm." Image and Vision Computing, 23(3), 299-309. • Phillips JM, Liu R, and Tomasi C. (2007). “Outlier Robust ICP for Minimizing Fractional RMSD." International Conference on 3-D Digital Imaging and Modeling, 427-434. • Fitzgibbon AW. (2003). "Robust registration of 2D and 3D point sets." Image and Vision Computing, 21(13-14), 1145-1153. • Škerl D, Likar B, and Pernus F. (2006). "A protocol for evaluation of similarity measures for rigid registration." IEEE Trans Med Imaging, 25(6), 779-791. • Pluim JPW, Maintz JBA and Viergever MA. (2000). "Interpolation artefacts in mutual information-based image registration." Computer Vision and Image Understanding, 77(2), 211-223. • Lehmann TM, Gonner C, and Spitzer K. (1999). "Survey: Interpolation methods in medical image processing." IEEE Trans Med Imaging, 18(11), 1049-1075. • Meijering EH, Niessen WJ, and Viergever MA. (2001). "Quantitative evaluation of convolution-based methods for medical image interpolation." Med Image Anal, 5(2), 111-126.

41. Bibliography • Murphy MJ. (1997). "An automatic six-degree-of-freedom image registration algorithm for image-guided frameless stereotaxic radiosurgery." Med Phys, 24(6), 857-866. • Wunsch P, and Hirzinger G. (1996). "Registration of CAD-models to images by iterative inverse perspective matching." International Conference on Pattern Recognition, 78-83. • Lavallée S, and Szeliski R. (1995). "Recovering the position and orientation of free-form objects from image contours using 3D distance maps." IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(4), 378-390. • Frisken SF, Perry RN, Rockwood AP, and Jones TR. (2000). “Adaptively sampled distance fields: a general representation of shape for computer graphics.” SIGGRAPH, 249-254. • LaRose DA. (2001). "Iterative X-ray/CT Registration Using Accelerated Volume Rendering," PhD, Carnegie Mellon University, Pittsburgh, PA. • Knaan D, and Joskowicz L. (2003). "Effective intensity-based 2D/3D rigid registration between fluoroscopic X-ray and CT." Medical Image Computing and Computer-Assisted Intervention (MICCAI 2003), Springer-Verlag, Montreal, Canada, 351-358. • Russakoff DB, Rohlfing T, Mori K, Rueckert D, Ho A, Adler JR, Jr., and Maurer CR, Jr. (2005). "Fast generation of digitally reconstructed radiographs using attenuation fields with application to 2D-3D image registration." IEEE Trans Med Imaging, 24(11), 1441-1454. • Khamene A, Bloch P, Wein W, Svatos M, and Saur F. (2006). "Automatic registration of portal images and volumetric CT for patient positioning in radiation therapy." Med Image Anal, 10(1), 96-112. • Mori K, Deguchi D, Sugiyama J, Suenaga Y, Toriwaki J, Maurer CR, Jr., Takabatake H, and Natori H. (2002). "Tracking of a bronchoscope using epipolar geometry analysis and intensity-based image registration of real and virtual endoscopic images." Med Image Anal, 6(3), 321-336. • Turgeon GA, Lehmann G, Guiraudon G, Drangova M, Holdsworth D, and Peters T. (2005). “2D-3D registration of coronary angiograms for cardiac procedure planning and guidance.“ Med Phys, 32(12), 3737-3749. • Livyatan H, Yaniv Z, and Joskowicz L. (2003). "Gradient-based 2D/3D rigid registration of fluoroscopic X-ray to CT." IEEE Trans Med Imaging, 22(11), 1395-1406.

42. Bibliography • Tomazežič D, Likar B, Slivnik T, and Pernus F. (2003). "3-D/2-D registration of CT and MR to X-ray images." IEEE Trans Med Imaging, 22(11), 1407-1416. • Fitzpatrick JM, West JB, and Maurer CR, Jr. (1998). "Predicting error in rigid-body point-based registration." IEEE Trans Med Imaging, 17(5), 694-702. • Moghari MH, and Abolmaesumi P. (2006). "A high-order solution for the distribution of target registration error in rigid-body point-based registration." Ninth International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI'06), Springer-Verlag, Copenhagen, Denmark., 9(Pt 2), 603-611.