Real-time B-Mode Ultrasound Image Simulation and Artifacts Modelling of Needles and Brachytherapy Seeds by Mengchen Zhu B.Eng., Zhejiang University, 2007 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2009 © Mengchen Zhu 2009 Abstract The training of ultrasound guided interventional procedures could benefit from ultrasound simulators that synthesize realistic-looking B-mode images of tissue and interventional tools, such as the needle. For a prostate brachytherapy simulator in particular, both the needle and seed images need to be synthesized. In this thesis, we propose an interpolation-based method for simulating needle and seed B-mode images in real time. We parametrize a needle or a seed B-mode image as a function of its position and orientation. We collect needle and seed images under various spatial configurations in a water-tank using a needle guidance robot. Then we use multi-dimensional tensor-product interpolation to simulate images of needles and seeds with arbitrary poses and positions using collected images. After further processing, the interpolated needle and seed images are superimposed on top of phantom or tissue image backgrounds. The similarity between the simulated and the real images is measured using a correlation metric. The results are also compared to in vivo images. Images in both the transverse and the sagittal planes are simulated. Artifacts associated with the needle and seed are studied in detail by experiments and by physical models. It is demonstrated that the bright band pattern associated with the sagittal needle images could be interpreted using a plate reflection model; the bright tails associated with the transverse needle images could be analyzed by simulating the back-scattered stress field from a sphere. ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Abstract List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Motivations and overview . . . . . . . . . . . . . . . . . 1.2 Background and literature review . . . . . . . . . . . . 1.2.1 Prostate brachytherapy . . . . . . . . . . . . . . 1.2.2 B-mode Ultrasound image simulation . . . . . . 1.2.3 Needle image characterization and visibility . . 1.2.4 Needle image simulation . . . . . . . . . . . . . 1.2.5 Image-based rendering . . . . . . . . . . . . . . 1.2.6 Ultrasound artifacts . . . . . . . . . . . . . . . . 1.2.7 Scattering from cylindrical and spherical objects 1.2.8 Seed detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Needle and brachytherapy seeds in B-mode ultrasound: realtime simulation and experimental studies . . . . . . . . . . 2.1 Transverse needle image simulation and tip visibility study . 2.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Simulation methods . . . . . . . . . . . . . . . . . . . 2.1.3 Data collection method . . . . . . . . . . . . . . . . . 2.1.4 Superimposition and shadow simulation . . . . . . . . 2.1.5 Results and validation . . . . . . . . . . . . . . . . . . 1 1 2 2 2 4 5 6 7 8 10 11 11 11 12 18 20 23 iii Table of Contents 2.2 2.3 2.4 Sagittal needle image simulation . . . . . . . . . . . . . . . . 2.2.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Results and Discussions . . . . . . . . . . . . . . . . . Seed image simulation in the transverse and sagittal plane . 2.3.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Results and Discussions . . . . . . . . . . . . . . . . . Experimental characterization of the transverse shadow forming conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 The bright tail . . . . . . . . . . . . . . . . . . . . . . 2.4.2 The dark shadow . . . . . . . . . . . . . . . . . . . . 2.4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 3 Modelling of the bright band artifact in 3.1 A description of the phenomenon . . . 3.2 Modelling of the phenomenon . . . . . 3.3 Validation of the model . . . . . . . . . 3.3.1 Sound speed . . . . . . . . . . . 3.3.2 Band intensity . . . . . . . . . . 3.4 Discussion . . . . . . . . . . . . . . . . 3.4.1 Incident direction . . . . . . . . 3.4.2 Quarter-wave layer . . . . . . . 3.4.3 Surface roughness . . . . . . . . 3.4.4 Relation to the needle image . . the sagittal plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Modelling of the bright tail and dark shadow artifacts in the transverse plane . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Hypotheses for the bright tail and dark shadow forming mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Boundary value problem for scattering from spherical shells . 4.4 Scattering spectrum simulation . . . . . . . . . . . . . . . . . 4.4.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Transient simulation . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Modelling of the ultrasound pulse . . . . . . . . . . . 4.5.2 Validation of US pulse modelling . . . . . . . . . . . . 4.5.3 Pulse response . . . . . . . . . . . . . . . . . . . . . . 4.5.4 Simulation of transient back-scattering from a spherical shell . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Validations and interpretations of the transient response . . 4.6.1 Wave speed . . . . . . . . . . . . . . . . . . . . . . . 26 26 29 30 30 35 37 37 41 42 43 44 44 47 48 50 53 53 53 53 54 55 55 56 57 60 63 64 65 67 67 67 68 68 iv Table of Contents 4.6.2 4.7 Resemblance to response from needle shaft in the transverse plane . . . . . . . . . . . . . . . . . . . . . . . . 4.6.3 Dark shadow from phase cancellation . . . . . . . . . 4.6.4 Solid vs. shell . . . . . . . . . . . . . . . . . . . . . . Numerical issues . . . . . . . . . . . . . . . . . . . . . . . . . 4.7.1 Convergence of the incident wave series . . . . . . . . 4.7.2 Condition number of the linear system . . . . . . . . 4.7.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 69 69 70 72 72 72 75 5 Implementation of the ultrasound simulator for interactive needle insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Notations and overview . . . . . . . . . . . . . . . . . . . . . 5.2 Re-packing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Transformation from the global coordinate to the probe coordinate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Problem description . . . . . . . . . . . . . . . . . . . 5.3.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Line crossing . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4.1 Problem description . . . . . . . . . . . . . . . . . . . 5.4.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 θ-correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5.1 Problem description . . . . . . . . . . . . . . . . . . . 5.5.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Depth calculation . . . . . . . . . . . . . . . . . . . . . . . . 5.6.1 Problem description . . . . . . . . . . . . . . . . . . . 5.6.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Needle insertion simulator . . . . . . . . . . . . . . . . . . . 79 79 79 80 80 80 81 81 83 83 83 83 85 6 Conclusion . . . . . . . 6.1 Summary . . . . . . 6.2 Contributions . . . 6.3 Future work . . . . . . . . 86 86 87 88 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 76 79 Appendices A A derivation of reflection and transmission coefficients at fluid-solid interface . . . . . . . . . . . . . . . . . . . . . . . . . 98 v Table of Contents A.1 A.2 A.3 A.4 Solid . . . . . . . . Fluid . . . . . . . . Boundary conditions Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 . 99 . 99 . 100 B Implementation of automatic image collection using the needle guidance robot . . . . . . . . . . . . . . . . . . . . . . . . . 102 vi List of Tables 1.1 Comparison of previous works on scattering from spherical and cylindrical solid and shell. The first group studies the generic scattering; the second group studies the scattering in the context of ultrasound contrast agents; the third uses scattering to model the comet tail artifact. . . . . . . . . . . . . . . . . . . . . . samples . . . . . . . . . in . . 9 2.1 2.2 2.3 Sample range and number of samples . . . . Sample range and numbers. . . . . . . . . . The range and numbers of the seed image sagittal and the transverse planes. . . . . . . . . . . . the . . . 19 28 4.1 Material properties used in the simulation . . . . . . . . . . . 64 35 vii List of Figures 1.1 1.2 1.3 1.4 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 A sketch of the geometry of ultrasound guided prostate brachytherapy (reprinted from [67]). The current technique employed by the BC cancer agency involves insertion of only one needle at a time, not as depicted in this figure. . . . . . . . . . . . . . . 3 The needle used for prostate brachytherapy. The arrow on the base indicates the bevel tip direction (reprinted from [30]). 3 The brachytherapy seeds (reprinted from [20]). . . . . . . . . 3 Lead pellets in the liver producing bright tail (reprinted from [45]). 7 Needle insertion process in the water and in the phantom . . Coordinate assignment (not to physical scale). ρ is the roll angle which measures how much the needle rotates around its axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Validation of θ invariance. After rotating the first image by θ, the two images have a correlation score of 0.96. . . . . . . Visualization of the interpolation process. Here we only draw 3 dimensions. Every vertex represents a sample image. The red cross represents the queried configuration (r, l, ρ). The algorithm first picks out the nearest rectangular sample subspace, then it bi-linearly interpolates the four sample images corresponding to the subspace vertices. . . . . . . . . . . . . Data collection setup . . . . . . . . . . . . . . . . . . . . . . . When the needle is not perfectly straight, rotating the needle around its axis will cause the needle image to deviate from the original position. . . . . . . . . . . . . . . . . . . . . . . . The interface used for registration. The pink dot indicates the “theoretical” position of the needle calculated from the needle configurations. User could override it by clicking using the cross-hair. . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparison of simulated needle images in phantom using different superimposition methods with the real image. . . . . . 12 13 14 17 19 20 21 22 viii List of Figures 2.9 2.22 2.23 Comparison of the simulated and the real needle images. In (a) through (e), upper-left quadrant: simulated image in water; upper-right: real image in water; lower-left: simulated image in phantom; lower-right: real image in phantom. Configurations shown in the subfigure captions are of the form (r, l, ρ, α, β). . . . . . . . . . . . . . . . . . . . . . . The needle in water exhibited the bright tail. The needle in phantom did not show bright tail pattern in this experiment. Comparison of simulated transverse needle image and in vivo needle image from a brachytherapy session. . . . . . . . . . . Coordinate assignment for the sagittal image simulation (not to physical scale) . . . . . . . . . . . . . . . . . . . . . . . . . User interface showing different interpolation methods used for image rotation . . . . . . . . . . . . . . . . . . . . . . . . Comparison of simulated and real needle images. In (a) through (d), upper-left quadrant: simulated image in water; upperright: real image in water; lower-left: simulated image in phantom; lower-right: real image in phantom. Configurations shown in the subfigure captions are of the form (x, y, z, ρ, α, β). . . . . . . . . . . . . . . . . . . . . . . . . . . Comparison of simulated and real needle images. In (a) through (e), upper-left quadrant: simulated image in water; upperright: real image in water; lower-left: simulated image in phantom; lower-right: real image in phantom. Configurations shown in the subfigure captions are of the form (x, y, z, ρ, α, β). . . . . . . . . . . . . . . . . . . . . . . . . . . Phantom with a seed mounted on a needle . . . . . . . . . . . Comparison of seed images in the transverse plane . . . . . . Comparison of simulated and real seed images. . . . . . . . . Needle image showing bright tail in a scatterer-free phantom Needle image in water with suspended scatterers. . . . . . . . Needle in a PVC phantom. . . . . . . . . . . . . . . . . . . . The load of the phantom may exert torque on the probe and change the pitch angle. . . . . . . . . . . . . . . . . . . . . . . Phantom fracture makes bright tail disappear. . . . . . . . . Needle in a pig liver; the needle location is pointed out. . . . 3.1 Sagittal B-mode images of needle in water . . . . . . . . . . . 2.10 2.11 2.12 2.13 2.14 2.14 2.15 2.16 2.17 2.18 2.19 2.20 2.21 24 25 26 28 31 32 33 34 35 36 37 38 39 39 40 41 44 ix List of Figures 3.2 3.3 3.4 3.5 3.6 3.7 3.8 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 Bright band when imaging in B-mode an aluminum plate and a steel plate at perpendicular angle with the transducer. The double arrow indicates the incident direction; the small arrow indicates the reverberation between the transducer and the plate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ray tracing of the reflections and transmissions; the lengths of the rays are roughly proportional to the amplitudes of the pressure/traction when the solid is aluminum. The time axis is from the left to the right. . . . . . . . . . . . . . . . . . . . Set-up for metal plate reflection study . . . . . . . . . . . . . RF lines corresponding to the B-mode images in Fig. 3.2; the position of the lines are both at around 3/4 of the B-mode image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparison of the theoretical and experimental results when we assume the RF amplitude is proportional to the pressure. Comparison of the theoretical and experimental results when we assume the RF amplitude is proportional to the sound intensity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Band pattern in the quarter-wave layer when scanning air. . . Band pattern visible in water, phantom and liver. . . . . . . . Scattering by a spherical shell. For simplicity,only three layers are drawn here. All layers are elastic solids. . . . . . . . . . Simulated spectrum of the back-scattered stress from a radius a = 1 mm solid sphere measured at a distance r = 60 mm from the sphere, 200 sample points. . . . . . . . . . . . . . . . Form function of a steel spherical scatterer in water (reprinted from [19] chapter 33). To the left of the added dash line is the Rayleigh scattering region. . . . . . . . . . . . . . . . . . Simulated impulse response corresponding to the back-scattered stress from a radius a = 1 mm solid sphere measured at a distance r = 60 cm from the sphere. The y-axis has the unit of traction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The simulated US pulse and power spectrum. . . . . . . . . . The experimentally measured power spectrum of the received RF signal (reprinted from [18]). . . . . . . . . . . . . . . . . . The simulated stress resulting from the simulated US pulse with the otherwise same parameters as in Fig. 4.5 . . . . . . . The transient pulse response from a spherical shell with a = 1 mm, shell thickness 0.02mm, measured at r = 60 mm. . . . 45 47 48 49 51 52 53 56 57 63 64 65 66 67 68 69 x List of Figures 4.10 The RF line of the back-scattering from the needle shaft. The length of the line is around 50 µs. The later half of the response is contaminated by the reverberation between the needle and the transducer. . . . . . . . . . . . . . . . . . . . . 4.11 The forward incident, forward scattered response, and their sum for a solid sphere with a = 5 mm, at r = 5 mm. The maximums are pointed out. . . . . . . . . . . . . . . . . . . . 4.12 The comparison between a cylindrical shell and a solid. . . . 4.13 The convergence of the incident wave series with varying r. The terms where the series start to converge are pointed out by arrows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 The contour of log10 of the condition numbers of the linear systems of displacement potential coefficients for each mode n and each frequency. . . . . . . . . . . . . . . . . . . . . . . 4.15 The Bessel functions of the first and second kind at ka = 5: jn (5) and nn (5). . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 5.2 5.3 5.4 5.5 5.6 Flow chart of the needle insertion simulation . . . . . . . . . The global frame, image frame, and probe frame. Compare to Fig. 2.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Calculation of the pitch (α) and yaw (β) angles from the probe and needle frame . . . . . . . . . . . . . . . . . . . . . Frames involved in the θ correction . . . . . . . . . . . . . . . The geometry of the needle rotating around the robot wrist rotation centre. . . . . . . . . . . . . . . . . . . . . . . . . . . The needle insertion simulator graphics interface (courtesy of Orcun Goksel). The needle is pointed out. . . . . . . . . . . . 70 71 71 72 73 74 77 78 82 82 84 85 A.1 Geometry of the problem . . . . . . . . . . . . . . . . . . . . 100 B.1 The image acquisition automaton. . . . . . . . . . . . . . . . 102 xi Acknowledgements I would like to first thank my thesis advisor, Dr. Tim Salcudean, who has provided invaluable guidance and feedback on my research throughout my two years at UBC. I am deeply grateful for the freedom he has given me to pursue the research topics I am interested in. I also appreciate his patience and understanding when ideas and results took weeks, even months to take shape. The joy of sharing good results will be remembered and cherished. Thanks also to the staff at ECE, UBC, whose technical support and professional help are essential to the completion of this thesis work. The results obtained in this thesis will be impossible without the support from my brilliant colleagues in the Robotics and Control Lab: thanks to Victor Wang and Diego Prananta for help on the robot; to Ali Baghani for sharing his spherical scattering simulator; to Orcun Goksel for various helpful discussions on the simulator; to Stanley Lee for converting the Matlab code; to Hani Eskandari for initial help on recording the RF signal from the needle and metal plates – without his “careless” act of throwing a steel ruler into my water tank, a whole chapter of this thesis might be missing; to Sam Frew for helpful discussions on transducer modelling and for providing New Zealand indie music which has accompanied me throughout the writing of this thesis (the thesis is therefore best read with an accent); to Jeff Abeysekera for lending me his probe holder and for allowing me to mock his design, which is actually great; to Jing Xiang for letting me jab her precious pig livers (for experimental purpose of course). In addition, help from Reza Zahiri-Azar, Xu Wen, Denis Tran, Sara Mahdave, Ehsan Dehghan, Mehdi Moradi, Ramin Sahebjavaher, Eric Yim, Farheen Taquee, Jack Zheng and Honza Vohradsky are all indispensable at various points of my research and I am eternally grateful to them all. Thanks to my friends here in Vancouver, with whom I have spent two of the best years in my life. The hiking trips, night skiing, dragon boating, ocean swimming, and various parties, live concerts, performances, museumvisits will all be remembered and treasured. Thanks also to my old pals around the globe. Special thanks to Luyi for her emotional support and understanding from across the border. xii Acknowledgements Thanks to my parents for their unconditional support from afar. I am eternally indebted to them. Last but not least, thanks to all the books that have stayed on my bookshelf for the last two years. The joy, strength, and wisdom I have gained from them are immeasurable. xiii — To the Pacific Ocean and the singing cetaceans within. xiv Chapter 1 Introduction 1.1 Motivations and overview Ultrasound imaging has become increasingly popular for guiding medical interventional procedures [14, 21, 53], in which minimally invasive instruments (usually needles) are inserted for various purposes, such as brachytherapy and biopsy. In these procedures, the subjects of ultrasound imaging are no longer solely the biological tissue, but also include various foreign metallic objects such as the needle and seeds. For these objects, the assumptions of the medical ultrasound instruments no longer hold and various artifacts may occur accordingly [38]. Partly due to this complication, in an interventional procedure, judging the precise position of the tool in the ultrasound and guiding the tool to the target is not an easy task. The need for training arises. The traditional apprenticeship training model has several drawbacks such as potential patient discomfort and complications [76]. Simulationbased training is therefore desired. In such simulators, both the ultrasound (visual) and haptics (feel) may be simulated. Between the two, the visual part has been shown to be particularly important for performing tasks such as needle insertion [26]. To achieve realistic ultrasound image simulation for interventional procedures such as needle insertion, both the tool (e.g. needle) and the tissue need to be simulated and their interactions modelled correctly. One goal of this thesis is to develop a real-time image simulation method for the needle and the seed. The method should be compatible with the existing image simulation method for the tissue [27], so that together they constitute an interactive needle insertion image simulator. Another goal of this thesis is to investigate physical models of the artifacts formation in needle/seed B-mode images. These models should be able to explain why some peculiar patterns occur in ultrasound B-mode images of the needle and the seed. An outline of the thesis is as follows: 1 1.2. Background and literature review In Chapter 2, a real-time needle and seed image simulator is developed in the context of prostate brachytherapy. The simulation takes into account all spatial degrees of freedom of the needle and seed and also the thickness of ultrasound imaging plane, therefore could correctly reflect the positioning of the needle and seed. To better understand the bright band/tail artifact apparent in some needle images, in Chapter 3 and 4 two basic scattering problems are studied and computational models developed. The results demonstrate that the bright band/tail could be modelled using linear acoustics. In Chapter 5, the integration of the needle image simulation developed in Chapter 2 and previously developed deformable tissue image simulation [27, 28] is discussed. Algorithms for the integration are described. In the following of this chapter, background knowledge and critical reviews on several topics are given. 1.2 1.2.1 Background and literature review Prostate brachytherapy Reviews of the prostate brachytherapy procedure are available in [13, 47]. In Fig. 1.1 a sketch of the set-up is shown. During the procedure, pre-loaded needles are inserted through a perineal template percutaneously into the pre-planned positions in the prostate under ultrasound guidance using the TRUS probe. The position of the needle is mainly determined by localizing the needle tip in the transverse view of the ultrasound; the sagittal view and manual measurement are used to confirm the needle position. After the implant of the seeds, the distribution may be checked by fluoroscopy. The needle used in the procedure is shown in Fig. 1.2. The needle is 18 gauge, 20 cm long, with a diameter of 1.25 mm; the stylet for depositing the loaded seeds is of the same length, with a diameter of 0.92 mm. Several kinds of seed may be used in the procedure and are shown in Fig. 1.3. In this thesis, the standard seeds are used in the experiments and simulation. The standard seed has a length of 4.6 mm and a diameter of 0.8 mm. 1.2.2 B-mode Ultrasound image simulation We classify the previous works on B-mode ultrasound simulation into three categories: the acoustics-based methods use acoustic models explicitly in the simulation; the pseudo-acoustic models use computer graphics techniques 2 1.2. Background and literature review Figure 1.1: A sketch of the geometry of ultrasound guided prostate brachytherapy (reprinted from [67]). The current technique employed by the BC cancer agency involves insertion of only one needle at a time, not as depicted in this figure. Figure 1.2: The needle used for prostate brachytherapy. The arrow on the base indicates the bevel tip direction (reprinted from [30]). Figure 1.3: The brachytherapy seeds (reprinted from [20]). 3 1.2. Background and literature review to simulate the ultrasound image formation; the interpolation-based models generally involve slicing a pre-scanned 3D ultrasound volume. Acoustics based The most well-known work in this category is the FIELD II simulator [40– 42]. It uses linear acoustic models and the spatial impulse response to simulate the transducer pressure field and the backscattered field from a finite number of point scatterers. Pseudo-acoustic model In [75, 76], computer graphics techniques were used to generate artificial ultrasound images from CT for needle insertion simulation. Shadow, noise, and reflections were simulated. The needle was simply rendered as a bright line and was only simulated for the longitudinal plane. [50, 80, 81] are another series of research in needle insertion simulation. In these works, artificial ultrasound images were generated from CT scans. Speckles in the ultrasound were simulated by Gaussian noise, shadows by ray-casting, and radial blur by alpha blending. Needle rendering was only done in the longitudinal plane, with transparency defined by the intersection volume with the image plane of tunable thickness. In [69], attenuation, scattering, reflection, speckles, and brightness of the ultrasound image were simulated by using various computer graphics techniques and statistic models . Interpolation based The classical way of simulating ultrasound images by interpolation is by slicing a pre-scanned 3D volume, using techniques surveyed in [64]. [1, 22] and those surveyed in [54] all used this approach. 1.2.3 Needle image characterization and visibility To motivate the development of our needle image simulation approach, we review previous clinical studies on needle image characteristics. [16] provides a good review of the clinical aspects of ultrasound imaging of needle insertion, including basic imaging physics, artefacts, transducer, visibility, and guidance techniques. One of the early studies on needle visibility is [11]. It investigated the needle tip echo intensity in water under various incident angles, bevel angles 4 1.2. Background and literature review (see Fig. 1.2), tip opening directions, and needle diameters using a linear probe. It was found echo intensity decreases linearly with the decreasing ultrasound incident angle (from vertical to parallel) and needle diameter, but does not change much with the bevel angle. Tip opening directions were also found to have influence on the echo intensity. In [68], needle visibility was studied using a linear probe in a gelatine phantom. Various kinds of needles were tested. Both the transverse and longitudinal images were collected. The tip and shaft images in the transverse plane were compared. The relative intensities of the needle and subjective visibility were measured. In [58], a linear probe was used to scan the needle in beef gravy phantom. The focus of the study was the effect of incident angle on visibility. Common findings of these previous studies are that the needle image depends highly on the spatial orientation of the needle with respect to the transducer (view-dependency) and that the needle tip is distinctively echogenic. 1.2.4 Needle image simulation Based on the discussions of Section 1.2.3, a realistic needle image simulator should be able to model the view-dependency of the needle images. The view dependency is an important feature for the needle insertion simulator also because in the real procedures, the doctors exploit this dynamic effect to judge the position of the needle, especially in the transverse plane [16]. For example, in prostate brachytherapy [47], doctors insert the needle until they see the needle tip, then rotate the needle to insure that the tip is at the desired transverse position. We examine the previous ultrasound simulation methods discussed in Section 1.2.2 to see if they are suitable for simulating the needle images. We should note that the previous works on the ultrasound simulation are designed primarily for the simulation of soft tissue B-mode images. 1. Acoustics-based simulation, specifically FIELD II: the point scatterer assumption of FIELD II does not hold for the continuous, relatively large objects such as the needle and seed that we are considering in this thesis. The geometry of the larger object needs to be explicitly modelled as is done in chapter 3 and 4. Also, the physically-based approach is slow, thus not suitable for simulations requiring real-time performance. 5 1.2. Background and literature review 2. Pseudo-acoustic methods: the previous works on needle insertion simulation typically opt for this approach. As we have noted, these works render the needle as a bright line without much consideration of the view-dependency of the needle and the echogenic tip. 3. The volume slicing approach cannot deal with the view-dependent features of the needle images. Indeed, looking at the needle from another angle is not equivalent to rotating the original image. It is in fact known in general that volume slicing does not retain the look of ultrasound artifacts [8]. As a conclusion, none of the existing ultrasound simulation techniques can meet our requirement of simulating view-dependent objects in real-time. We propose an multi-dimensional interpolation-based approach that is a natural extension to the classical 3D slicing technique: instead of interpolating in 3D, we parameterize the image using the full 6 degrees of freedom spatial parameters and interpolate in 6D. 1.2.5 Image-based rendering The idea of using images of an object under various poses and positions to simulate new object poses has been explored in the field of image-based rendering in computer graphics. There are two paradigms in computer graphics: physical model based versus image (memory) based [62]. The model based approach is more versatile, but may not achieve the desired effect; also the time required for rendering might be unacceptable for real-time applications. The image based method is an attractive alternative. It achieves 3D rendering by synthesizing new views from the stored views by approximation/interpolation. The image-based rendering for graphics typically involves the interpolation on a subspace of the 7-D plenoptic function. A review can be seen in [79]. In [49], the 7D plenoptic function was reduced to 4D. The method of image collection using a camera gantry is similar in spirit to our robotic image collection method. However, the ray-based parameterization method is not applicable in our case, because in [49] the camera is moving and the light source is fixed, whereas in our case the transducer serves as both the camera and the light source, therefore no relative motion is present. 6 1.2. Background and literature review Figure 1.4: from [45]). 1.2.6 Lead pellets in the liver producing bright tail (reprinted Ultrasound artifacts As we will see in the body of the thesis, various ultrasound artifacts may occur in the needle and seed images. Two artifacts that we are particularly interested in are the bright band/tail and the dark shadow (see e.g. Fig. 2.23). The former is also called the comet tail, ring-down, or the dirty shadow; the latter is also called the acoustic shadow or the clean shadow. In the following papers, these two artifacts are described, the formation mechanism postulated, but the models are largely qualitative. [82] documents the “discovery” of the comet tail artifact. The formation of the artifact was ascribed to the reverberation inside the metal objects. The acoustic shadow was observed to co-exist with the comet tail. In [5], the bright tail from air bubbles was studied. It was found that the tail originates from the fluid “bugle” formed in between two layers of air bubbles. The term “ring-down” artifact comes from this study. In [65], the bright and dark tail were studied experimentally and were shown to be related to the surface roughness and curvature of the insonified object. [45, 46] survey ultrasound artifacts and their formation. Both the bright tail and the dark shadow are examined. See Fig. 1.4 for bright tails generated from lead shots. In [38], artifacts from a steel rod were studied in water and used to interpret artifacts from minimally invasive instruments. The tip of the cylindrical rod was observed to generate distinct comet tail artifact, with band interval proportional to the diameter of the rods. This was postulated to be due to the circumferential wave. 7 1.2. Background and literature review In [34], the scattering from metal spheres was physically modelled to interpret the comet tail artifact. However, it did not follow the classical formulation of the fluid-solid acoustic interaction. For example, no shear wave was considered in the solid. Therefore care should be taken to interpret the results from this simulation study. To conclude, we note that the bright tail artifact is usually explained by reverberation, but no quantitative modelling is known except for [34], which did not model the shear wave in a solid. 1.2.7 Scattering from cylindrical and spherical objects The sphere and the cylinder are the two “canonical” geometries studied in acoustic scattering [19]. Similar to the strategy taken in [34], we hope to model the bright tail and dark shadow physically by modelling the scattering from the spherical or cylindrical objects. A partial listing of several classical and more recent studies on the spherical and cylindrical solid/shell scattering can be seen in Table. 1.1. In these various modelling efforts, boundary value problems are usually formulated to derive the partial wave series representations of the scattered pressure field as a function of the observing position and ka, which is a function of frequency and object radius. The steady-state solution of the scattered field follows. It is also possible to approximate the transient solutions from the steady-state ones. It should be noted that higher level approximations such as the creeping wave formulation (see chapter 4 and 33 of [19]) are also useful for interpreting the scattering. In this formulation, the acoustic waves are interpreted to penetrate into the solid, reflect at the surface, and propagate on the surface at the same time. More theoretical analysis of the interactions between the sound wave and structures including beam, thin plate, spherical and cylindrical shells can be found in [43]. 8 Steady-state/transient steady-state transient & steady-state transient & steady-state transient & steady-state steady-state steady-state transient steady-state steady-state steady-state steady-state steady-state transient Geometry spherical & cylindrical solid spherical solid cylindrical solid & shell spherical solid spherical shell finite cylinder spherical shell spherical & cylindrical shell spherical shell spherical shell spherical shell spherical shell spherical solid Method theory & experiment experiment experiment theory theory theory theory theory theory theory theory theory theory Table 1.1: Comparison of previous works on scattering from spherical and cylindrical solid and shell. The first group studies the generic scattering; the second group studies the scattering in the context of ultrasound contrast agents; the third uses scattering to model the comet tail artifact. Author Faran [23] Hampton et al [32] Barnard et al [9] Hickling [36] Ayres et al [6] Stanton [71] Gaunaurd et al [24] Hasegawa et al [35] Allen et al [2] Ye [78] Hu et al [37] Baghani [7] Harpen [34] 1.2. Background and literature review 9 1.2. Background and literature review 1.2.8 Seed detection The radioactive seeds used in prostate brachytherapy exhibit similar (though weaker) B-mode image pattern to that of the needle. The body of the seed may cast a short bright tail when the incidence is close to orthogonal and a dark shadow when there are scatterers in the background. The bright tail, as will be demonstrated in this thesis, is related to the RF frequency spectrum of the scatterer (seed in this case). It is therefore natural to suggest seed detection methods based on spectrum analysis. This has been demonstrated in several previous works. In [51], singular spectrum analysis was applied to detect seeds. This was later extended to deal with angle-dependency of the seed in [52]. In [20], the backscattering signal from different seed types under different angles were compared; corrugated seeds were shown to be less angle dependent. In [55], vibro-acoustography was demonstrated to be superior in imaging seeds than traditional ultrasound. However the vibro-acoustography technique requires two sound sources driven at slightly different frequencies, thus may not be readily applicable using existing transducers. It should be noted that the above methods were only experimentally validated; no physical scattering model of the seed was developed in these works. We also note that previous literature on underwater target identification (such as [25]) might be relevant for the purpose of seed detection. 10 Chapter 2 Needle and brachytherapy seeds in B-mode ultrasound: real-time simulation and experimental studies In this chapter, a general real time image simulation method is developed to enable 6DOF interaction with an object in ultrasound images. The method is first developed in the context of needle image simulation in the transverse plane and later extended to deal with needle image simulation in the sagittal plane and brachytherapy seed image simulation in both planes. 2.1 2.1.1 Transverse needle image simulation and tip visibility study Introduction In procedures such as the prostate brachytherapy, the needle is inserted perpendicular to the transverse plane of the transducer. By discerning the advancing needle tip in the B-mode ultrasound, the doctor can ensure that the desired spatial position is reached by the needle before the radioactive seeds are deposited. The asymmetry of the needle tip results in changing needle images when the needle position and/or orientation is changed. The doctor exploits this to better localize the position of the needle during the insertion: s/he will actively rotate and/or rock the needle back and forth to judge the position of the needle from the changing needle and surrounding tissue images [53]. See Fig. 2.1 for an example of the appearance of the needle image during the insertion process. In the ultrasound simulation for the needle insertion procedures, it is therefore important to take account of the changes in the needle image as the needle’s spatial configurations change. In this section, we concentrate 11 2.1. Transverse needle image simulation and tip visibility study { tail body (a) The needle tip in the water (b) The needle shaft in the water (c) The needle tip in the phantom (d) The needle shaft in the phantom Figure 2.1: Needle insertion process in the water and in the phantom on simulating the transverse needle images with different needle configurations. The changes in the image of tissue reacting to the needle insertion are investigated and simulated in other works [27, 28]. Several issues regarding the integration with these previous works are discussed in Chapter 5. We further note that the image of the needle in the water resembles that in the phantom/tissue (Fig. 2.1). The basic idea of the simulation is therefore to first simulate the needle images in the water by interpolating sample images collected in water and then superimpose the simulated images onto the tissue background. A validation based on correlation is performed to justify this method. In Section 2.4 the differences between the needle image in water and in phantom/tissue are investigated experimentally. To aid the discussion, we call the echogenic spot of the needle image corresponding to the real needle position the body and the artifact distal to it the tail of the needle image. 2.1.2 Simulation methods We make the assumption that the needle image depends only on the local spatial configuration (position and orientation) of the intersection of the needle with the imaging plane. In this manner, we can represent the planar 12 2.1. Transverse needle image simulation and tip visibility study Imaging plane bevel tip l Y X e Xn r α β {P n , C n } Yn Zn Z Zn θ Y ρ X {Op , C p } e Z TRUS probe Figure 2.2: Coordinate assignment (not to physical scale). ρ is the roll angle which measures how much the needle rotates around its axis. needle image as a function In : R6 → R2 of the six degrees of freedom of the needle. In this section, we concentrate on the simulation of the needle images in the transverse plane, which is the major imaging plane used for procedures such as prostate brachytherapy. The parameterization discussed below is suitable for the curvilinear geometry such as the transverse plane of the transrectal ultrasound (TRUS) transducer; for a linear transducer, the Cartesian coordinates can be used in place of the polar coordinates. Parameterization The parameterization for the needle configuration with respect to the TRUS probe transverse plane is as follows: considering a general beam profile where the imaging plane has a thickness along the Z direction, we define the characteristic imaging plane to be the plane with the largest Z coordinate, as shown in Fig. 2.2. This plane is determined experimentally by inserting the needle until the tip is barely seen. We attach a Cartesian frame {Op , C p } e to the intersection of this imaging plane and the probe axis, and another Cartesian frame {P n , C n } to the intersection of the imaging plane and the e needle axis. Zn coincides with the needle axis and Yn points to the bevel direction. We use (r, θ) to denote the polar coordinates of P n in {Op , C p }, e and l to denote the needle depth. We describe the rotationetransformation between C p and C n using roll-pitch-yaw angles (ρ, α, β). The needle image is a function In (r, θ, l, ρ, α, β) of these variables. We could show experimentally that the needle image is invariant with 13 2.1. Transverse needle image simulation and tip visibility study (a) (b) Figure 2.3: Validation of θ invariance. After rotating the first image by θ, the two images have a correlation score of 0.96. respect to θ, that is, if we keep other spatial parameters the same and vary θ, the image will only rotate by θ but will not change in intensity or shape otherwise (see Fig. 2.3). Mathematically we may write In (r, θ, l, ρ, α, β) = Rθ (In (r, 0, l, ρ0 , α0 , β0 )) (2.1) where Rθ is a spatial 2D rotation transform of the image with respect to the ultrasound probe centre; (ρ, α, β) in the local coordinate and (ρ0 , α0 , β0 ) in the global coordinate are related by the same rotation transform (more on these variables in Section 5.5). This symmetry is due to the curvilinear geometry of the TRUS probe we have used. To obtain images at arbitrary θ, we rotate the image at θ = 0. Thus the function to be sampled and interpolated is In (r, 0, l, ρ0 , α0 , β0 ) : R5 → R2 . Simulation using tensor-product interpolation The speed of interpolation is of major concern. Nearest neighbour and linear interpolation are two of the fastest interpolation methods available [73], but cannot be used directly for a multi-dimensional function In . However, as demonstrated below, we can use them for each dimension separately following an interpolation scheme. Since the images are digital, we can treat them one pixel at a time I : R5 → R. The task of interpolation is therefore to estimate the value of an arbitrary I(r, l, ρ, α, β) using sample values of I. The samples may be scattered or may be collected on a grid. For the former case, interpolation using radial basis functions is the usual choice. However, the computation 14 2.1. Transverse needle image simulation and tip visibility study overhead is huge. We therefore choose to sample over a multi-dimensional grid and use tensor-product interpolation. An introduction to this interpolation method could be seen in [3]; in the following we provide a similar development that is tailored for our needs. Tensor-product interpolation could be understood as a generalization of the bi-linear (2D) interpolation to n-D and to any order polynomial basis functions. It provides a scheme to use 1D piecewise polynomial basis functions for the interpolation of n-D functions. In 1D, we interpret an interpolation process as approximating the original function f (x) being sampled with a new set of basis functions {li (x)}. Following the notation in [3], we can express this process as applying an interpolatory projection operator π to the function f (x). If we sample f (x) on a grid {xi } and use the piecewise polynomials li (x) of order n as the nodal basis functions, we can encapsulate the 1D interpolation process as an approximation: X (πf )(x) = fˆ(x) = f (xi )li (x), π : f → fˆ i As special cases, when n = 0, li (x) are the shifted “box” functions and π corresponds to the nearest neighbour interpolation operation; when n = 1, li (x) are the shifted “hat” functions and π corresponds to the linear interpolation operation. We can show that the error bound decreases with decreasing h and increasing n: ||f −πf ||∞ ≤ C(n) ||f (n+1) ||∞ hn+1 , where h is the mesh size, provided that f has (n+1) bounded continuous derivatives. The extension to the higher dimensions is made possible by defining the “partial” interpolatory projections, which resemble in spirit the partial derivatives. For example if we know the sampled values of a function g : R2 → R on an M × N grid, we can approximate the two dimensions of the original g separately using two vector spaces li1 (x) and lj2 (y), given that the other dimension does not change: (π1 g)(x, y) := M X g(xi , y)li1 (x) i=0 (π2 g)(x, y) := N X g(x, yj )lj2 (y) j=0 Now we can define a new projection by combining the two partial projections 15 2.1. Transverse needle image simulation and tip visibility study π1 and π2 as follows: (π1 π2 g)(x, y) := M X N X g(xi , yj )li1 (x)lj2 (y) i=0 j=0 which has the property: (π1 π2 g)(x, y) = (π2 π1 g)(x, y) (2.2) We can show that the function space spanned by {li1 (x)lj2 (y)} is a vector space, therefore the process is again an interpolation, this time in 2D. This vector space is called a 2D tensor-product space. The approximation error can be shown to have the bound n+1 ∂ m+1 g ∂ g n+1 m+1 h , h (2.3) ||g − π1 π2 g||∞ ≤ C(n) max ∂xn+1 ∞ ∂y m+1 ∞ given that the partial derivatives exist and are bounded and continuous. The extension to n-D is immediate. Under this scheme, interpolating an n-D function I sampled on a grid can be decomposed into n separate steps, each of which involves interpolating on 1D when holding other dimensions constant. To ensure no aliasing occurs, the mesh size h should be chosen to be small enough when the corresponding partial derivative in Eq. 2.3 is large, that is, when I has high frequency components along that dimension. It was found experimentally that the variations of I along r and l have high frequency components; the scheme calls for finer grids in these two dimensions. To meet our speed requirement, only the nearest neighbour and the linear interpolation are used for each dimension. This means we will be using at most two neighbouring sample nodes for each dimension. Also, from Eq. 2.2, we know that the sequence of interpolation does not matter, so we compute the cheaper nearest neighbour first. Operationally, we could write our interpolation process as follows (see Fig. 2.4 for a visualization of the process): 1. Find the enclosing hypercube in the 5-dimensional sample space that contains the queried coordinate: (r, l, ρ, α, β) ∈ [ra , rb ]×[lc , ld ]×[ρe , ρf ]× [αg , αh ] × [βi , βj ]. 2. Find the nearest neighbours for r, α, β; say ra , αh , βj respectively; 16 2.1. Transverse needle image simulation and tip visibility study (r, l, ρ) (ra , l, ρ) ρ (rb , lc , ρe ) (ra , ld , ρe ) (ra , lc , ρe ) ρf r ρe rb ra l ld ρf (ra , l, ρ) lc ρe ld lc Figure 2.4: Visualization of the interpolation process. Here we only draw 3 dimensions. Every vertex represents a sample image. The red cross represents the queried configuration (r, l, ρ). The algorithm first picks out the nearest rectangular sample subspace, then it bi-linearly interpolates the four sample images corresponding to the subspace vertices. 3. In the rectangular subspace [lc , ld ] × [ρe , ρf ]|(ra ,αh ,βj ) ∈ R2 , use linear ˜ This is just bi-linearly interpolation to get the interpolant image I. interpolating I(ra , lc , ρe , αh , βj ), I(ra , lc , ρf , αh , βj ), I(ra , ld , ρe , αh , βj ), and I(ra , ld , ρf , αh , βj ); 4. Now I˜ has polar coordinates (ra , 0) in the imaging plane. Translate and rotate I˜ to the point (r, θ) on the imaging plane to finally get the approximated I(r, θ, l, ρ, α, β). The choice between the nearest neighbour and the linear interpolation for each dimension is based on the fact that during the needle insertion, the 17 2.1. Transverse needle image simulation and tip visibility study needle position on the imaging plane and its orientation are constrained by the tissue, therefore not likely to change much. It is the depth and roll angle that are under the doctor’s control and therefore responsible for the continuous changes of the needle images. Using linear interpolation for l, ρ will enable us to simulate these changes, and using nearest neighbour for other dimensions will speed up the interpolation considerably. Furthermore, using nearest neighbour for r will eliminate false artifacts caused by averaging two images with features at different locations (e.g. the bright tail discussed in Section 2.4.1). It can be shown that the first 3 steps of the above algorithm has a computational cost of O(Np ), where Np is the number of pixels of the images being interpolated. In the current implementation Np = 80 × 350 and is the size of needle images cropped from the original 640 × 480 images (the rest is black background). See Fig. 2.7 for the relative sizes of the needle image and the whole image. The speed of the last step largely depends on the interpolation method used for the spatial image rotations. In Matlab, using the nearest-neighbour will result in real time performance but may introduce aliasing artifacts or “jaggies”. See Fig. 2.13 for this artifact in the sagittal simulation, where it is more pronounced. This is eliminated in the OpenGL implementation (see Chapter 5) where linear interpolation can be used in real time. On a PC with an Intel 2.2GHz dual core CPU, the overall algorithm can achieve an average frame rate of 44 FPS for a resolution of 640 × 480 in a Matlab implementation without major optimization efforts. 2.1.3 Data collection method All data collection was done in a water tank of size LxWxH = 50x30x20cm filled with degassed water at 20. Brachytherapy needles (Bard , 18 gauge) were imaged with a TRUS probe connected to an ultrasound machine (Ultrasonix SonixRP). The machine parameters were set to be the same as in a brachytherapy procedure, except for disabled Time Gain Control, single focus, and sound speed c = cwater . The time gain control was disabled during the data collection because in the water there is very little attenuation with increasing depth; single focus was used because multiple focal points may sometimes cause split image in water. We collected the samples in the water tank rather than in the phantom because in water we can get clear, well-defined needle images, and because it is easier to move from one needle configuration to another without leaving needle tracks. ® 18 2.1. Transverse needle image simulation and tip visibility study Robot Wrist Stepper Needle Water tank TRUS probe Figure 2.5: Data collection setup Variable r l ρ α β Range 3 ∼ 6.5cm 1 ∼ 22mm 0 ∼ 180◦ 0 ∼ 9◦ 0 ∼ 15◦ # of samples 8 12 3 4 3 Table 2.1: Sample range and number of samples A needle guidance robot previously developed [67] was used for the data collection. Since the robot does not move in l or ρ, we controlled these two variables by hand. The robot was programmed to remotely connect to the SonixRP and capture images automatically after each desired configuration was reached. See Appendix B for more implementation details. The data collection set-up is shown in Fig. 2.5; the range and sample numbers of each dimension are shown in Table 2.1. The ranges of α, β were kept small, as in a real brachytherapy procedure. Only positive rotation angles were collected due to symmetry. To ensure that the sampling rate is high enough, some test images were first captured on a fine sample grid; then part of these images were used to interpolate the rest. Correlation was then used to compare the results, and the interpolated images were found to be similar to the real images when the sample numbers in Table 2.1 were used. As an example, for the roll angle ρ, when using images collected at every 180◦ to interpolate images of 19 2.1. Transverse needle image simulation and tip visibility study Water Tank Needle TRUS probe Imaging Plane Figure 2.6: When the needle is not perfectly straight, rotating the needle around its axis will cause the needle image to deviate from the original position. 45◦ ρ intervals, the correlation between the interpolated and the real images was 0.93; when using images collected at every 90◦ , the correlation was 0.97. Therefore, sampling at 90◦ interval is nearly the same as using more samples. Since our interpolation method requires image data samples at the exact sampling coordinate and instrumentation errors such as non-straight needle may cause misalignment (see Fig. 2.6), a rigid landmark based manual registration was performed after the sample collection to make sure that the images being interpolated are aligned. The landmark used was the bottom surface of the needle which was visible or could be inferred in all images. The result of the registration was that all sample images with the same r were brought into alignment along the bottom surface. The interactive interface used for the registration is shown in Fig. 2.7. As a side note, the relative needle depth l could be modified by either moving the probe or moving the needle. The choice will make a difference in the simulation, as discussed in Section 5.6. 2.1.4 Superimposition and shadow simulation To simulate the needle images in the tissue, background images of a gelatine tissue mimicking phantom were collected, onto which the interpolated 20 2.1. Transverse needle image simulation and tip visibility study Figure 2.7: The interface used for registration. The pink dot indicates the “theoretical” position of the needle calculated from the needle configurations. User could override it by clicking using the cross-hair. needle images in water were superimposed by addition. Superimposition by taking the maximum of the foreground and the background has also been tested. See Fig. 2.8 for a comparison. We see that addition with the background generates more realistic looking bright tail; a minor issue might be that we could see the background “through” the tail – this is not easily noticeable, however. In the current implementation we thus use addition as the superimposition method. It was observed in the phantom and in the tissue that dark shadows may occur distal to the main body of the needle image under certain conditions (see Fig. 2.22 for example). The conditions for the dark shadow formation are investigated experimentally in Section 2.4. The formation mechanisms are discussed in Chapter 4. To simulate the dark shadow, it was noted that the extent of the dark shadow resembles that of the bright tail formed by the needle shaft. We therefore apply several image processing techniques on the bright tail to simulate the dark shadow; the steps are as follows: 1. Use Ostu thresholding [59] to segment out the bright tail from the background; 2. Scale the intensity of the bright tail down by multiplying by a small constant so that the bright tail becomes dark; 3. Use Gaussian filtering to blur out the edges; 21 2.1. Transverse needle image simulation and tip visibility study (a) Superimposition by addition. (b) Superimposition by max. (c) Real Figure 2.8: Comparison of simulated needle images in phantom using different superimposition methods with the real image. 22 2.1. Transverse needle image simulation and tip visibility study 4. Superimpose the simulated shadow onto the background using a linear combination of intensity values. A sample result is shown in Fig. 2.9f where it is compared to the real shadow. 2.1.5 Results and validation A user interface has been developed to demonstrate the simulation of needle insertion in the transverse plane. Using this interface the user could change the spatial configurations interactively and observe the simulated images in real time. The interface is similar to the one designed for simulating the sagittal images, as shown in Fig. 2.13. Validation using phantom To validate the simulation method, 43 needle images in both the water tank and a homogeneous gelatine phantom with cellulose scatterers were collected using the same machine settings as used when the needle image samples were collected. The needle configurations were chosen to be different from those previously collected. Blank phantom background images were also collected before each insertion for use in the superimposition. Because images of the same modality are being compared, classical normalized cross-correlation is used as the similarity metric. In this validation, we are mainly concerned with the body part of the needle image, not the tail (see Fig. 2.1b for what the body and tail denote). One major reason is that in procedures such as prostate brachytherapy, the needle tip being imaged typically generates a very localized image (see Fig. 2.11). We do recognize the occurrence of the bright tail for the shaft part of the needle. This will be studied in Section 2.4 through experiments. It is sufficient to note here that our simulation method does simulate the tail part of the needle as could be seen in Fig. 2.8 and Fig. 2.9f. For each of the 43 image pairs, an 80×80 region around the body part of the needle was therefore cropped out and correlations were calculated. The mean correlation between the simulated and the real needle images in the water is corw = 0.86 out of 1, which gives an experimental measure of the error bound in Eq. 2.3. The mean correlation between the simulated (after superimposition) and the real needle images in phantom is corp = 0.79. One issue revealed in this particular experiment is that the bright tail in the phantom is not as visible as in the water (see e.g. Fig. 2.9e; see Fig. 2.10 for a fuller view.) As discussed in Section 2.4, this was possibly due to the instrumentation and the phantom. In subsequent experiments 23 2.1. Transverse needle image simulation and tip visibility study the bright tail invariably shows up in the phantom, thus the real correlation score should be higher. Also note that the correlation measure is sensitive to misalignment – this is why it is usually used for registration purposes. In our case, small misalignment always exists due to errors in instrumentation and manual registration. To achieve a better correlation score, an automatic registration could be performed to maximize the correlation. Needle visibility Since the collected sample database covers nearly all the physically attainable needle configurations during brachytherapy, it could be used to study the needle visibility in the transverse plane. By comparing the relative brightness (intensity) of the needle images collected, we can study the changes in needle visibility with respect to the needle configurations. Note Sim Real (a) (5.5, 6, 0, 0, 4) Sim Real (d) (4.5, 6, 135, 0, 0) Sim Real (b) (6.5, 7.5, 0, 2, 0) Sim Real (e) (5.2, 7.5, 0, 0, 0) Sim Real (c) (4.8, 6, 0, 0, 0) Sim Real (f) Shadow Figure 2.9: Comparison of the simulated and the real needle images. In (a) through (e), upper-left quadrant: simulated image in water; upper-right: real image in water; lower-left: simulated image in phantom; lower-right: real image in phantom. Configurations shown in the subfigure captions are of the form (r, l, ρ, α, β). 24 2.1. Transverse needle image simulation and tip visibility study (a) In water (b) In phantom Figure 2.10: The needle in water exhibited the bright tail. The needle in phantom did not show bright tail pattern in this experiment. that in this experiment when l = 6, it is the needle tip being imaged; when l = 7.5, it is the shaft being imaged. There are several observations we can make from Fig. 2.9: 1. By comparing Fig. 2.9c and 2.9e, we see that when the pitch angle α is small, the visibility increases with deeper insertion (larger l); 2. By comparing Fig. 2.9c and 2.9d, we see that the shape of the image changes with roll angle ρ when imaging the tip of the needle; 3. By comparing Fig. 2.9e and 2.9b, we see that visibility decreases with increasing pitch angle α; this is especially true for the shaft of the needle; 4. By comparing Fig. 2.9c and 2.9a, we see that with the increasing yaw angle β, the visibility of the needle image decreases slightly. Validation using in vivo images The simulation results were also compared with several in vivo needle images collected during a brachytherapy session. To obtain a similar background for superimposition, the needles in the in vivo images were carefully “erased” by copying surrounding background using the image editing software GIMP (GNU Image Manipulation Program). Then the simulated needle images were superimposed. Example results are shown in Fig. 2.11. Since there are no measurements of the needle position and orientation in-vivo, the needle configurations are chosen so that similar images can be obtained. 25 2.2. Sagittal needle image simulation Sim Real Figure 2.11: Comparison of simulated transverse needle image and in vivo needle image from a brachytherapy session. We see that the synthesized images do resemble the in vivo images. The difference may be attributed to the different ultrasound machines being used for data collection (an Ultrasonix machine was used for the sample image collection; a BK Medical machine was used in the hospital.) ® 2.2 ® Sagittal needle image simulation The same method used for simulating the transverse plane needle images can be used for simulating sagittal B-mode needle images. The results are not only useful for the simulation of brachytherapy, but are also generally applicable to other needle insertion scenarios where the needle is viewed in the sagittal (longitudinal) configuration, e.g. in biopsy, or anaesthesia. 2.2.1 Methods We follow the work-flow laid out in the previous section. The major differences lie in the parameterization. We assume the needle is straight and the imaging plane to be of finite thickness (not infinitely thin). Coordinate Assignment Referring to Fig. 2.12, we use the same XY Z coordinate system as was used in the transverse case – this is used for the real brachytherapy procedure. In the sagittal case, however, we use the Cartesian parameterization directly. 26 2.2. Sagittal needle image simulation We write the projection from the configuration coordinates to the image as In (P n , C n ) : R6 → R2 , or more explicitly In (x, y, z, ρ, α, β). The l e in the figure is used for simulating the insertion, which is discussed parameter in Section 2.2.2. For illustration purpose, we draw P 1 and P 2 as the points e where the needle axis enters and exits the imagingeplane; there might also be one or zero intersections if the needle is partially penetrating or entirely residing outside the imaging slab. Note that the thickness of the imaging plane may vary along the elevation (Y ) direction [70]; we are assuming the needle is kept near the focal point so that the plane thickness is nearly constant along the Y direction. The sagittal imaging plane is linear, in contrast with the transverse plane which is curvilinear. We can verify by experiments that when the focal point is kept near the needle, shifts of the needle in the Y and Z directions correspond to translations of the needle image in these directions; no intensity change is noticeable. We could therefore conclude that the images are isotropic along both the Y and Z directions (in the transverse plane the images are isotropic only along θ). Similar to Eq. 2.1, we can write the spatial transform that relates the needle images at different positions as In (x, y, z, ρ, α, β) = Ty (Tz (In (x, y0 , z0 , ρ, α, β))) where Ty , Tz stand for the spatial translation transforms (shifts) of the image from (y0 , z0 ) to (y, z) . This reduces our sample space to R4 . Also the image is digital so we consider one pixel at a time. The function we are interpolating is thus I(x, y0 , z0 , ρ, α, β) : R4 → R. The symmetry also indicates the following relation: I(−x, −β, :) = I(x, β, :) (2.4) where “:” indicates that dimension being the same. We can utilize this fact to collect only positive β. Data Collection The data collection procedures follow that in the transverse case. The range and sample numbers of each dimension is shown in Table. 2.2. We only collect positive β since the negative ones can be obtained by using Eq. 2.4; we collect only positive α because this is the common case in the real procedure. We observe that when x ∈ [−0.7, 0.7] cm the needle tip is inside the imaging slab; the image varies with the changing ρ. When x ∈ [0.7, 1.5] cm the needle tip is outside the imaging slab, therefore the image is ρ-invariant. 27 2.2. Sagittal needle image simulation {P n , C n } Yn e Zn Y Y x X α β {Op , C p } e Zn Imaging plane Xn l P0 P2 e X P1 eρ e Z TRUS probe Z Figure 2.12: Coordinate assignment for the sagittal image simulation (not to physical scale) We thus separate the image collection into these two groups. The collected images have resolution of 488 × 356. Template Preparation The images collected have the appearance of bright line segments with different “projected” pitch angles – the angle between the line segment and the horizontal image axis of the image plane: α0 = arctan( sin α ) cos α · cos β These images cannot be used directly for interpolation because they may not aligned. To align each image, we mark the needle tip position interactively, then translate the image so that the tip aligns with the “theoretical” tip Variable x ρ α β Range -0.7 ∼ 1.5cm 0 ∼ 180◦ 0 ∼ 6◦ 0 ∼ 15◦ # of samples 23 3 3 4 Table 2.2: Sample range and numbers. 28 2.2. Sagittal needle image simulation position calculated from the pre-plan and the pixel dimensions. We then use the calculated α0 to rotate the image around the new tip position by −α0 . After these procedures, all needle images are horizontal and aligned at the theoretical tip positions. Simulation by Interpolation We use the same tensor-product interpolation as used in the transverse case. With a similar consideration for speed and accuracy, we choose to use linear interpolation for x and ρ, and nearest neighbour for α and β. The resulting algorithm is as follows: 1. Find the enclosing hypercube in R4 that contains the queried coordinate: (x, ρ, α, β) ∈ [xa , xb ] × [ρc , ρd ] × [αe , αf ] × [βg , βh ]. 2. Find the nearest neighbours for α, β, say (αe , βh ); 3. In the rectangular subspace [xa , xb ] × [ρc , ρd ]|(αe ,βh ) ∈ R2 , linearly interpolate to get the interpolant I˜ ; 4. Translate I˜ so that the tip is at (y, z) and then rotate the image around the tip by α0 . We finally get the approximated I(x, y, z, ρ, α, β). 2.2.2 Results and Discussions Insertion Simulation To simulate the needle insertion insertion without using a physical user interface, we need to define a needle trajectory. In the current interface (as shown in Fig. 2.13), the user can move the needle tip freely around the achievable configuration space. After the user has moved the needle to the desired initial position and orientation {P 0 , C n }, the user could initiate an e insertion by changing the l parameter, which parametrizes the needle path, given the initial position and orientation. The orientation of the needle remains unchanged under the straight needle assumption. The position of the needle tip thus changes as follows: P (l) = P 0 + l · ẑn e e where ẑn is the normalized direction vector: ẑn = ~zn |~zn | 29 2.3. Seed image simulation in the transverse and sagittal plane and ~zn = (xZn , yZn , zZn ) = (cos α sin β, sin α, cos α cos β) With the configuration shown in Fig. 2.12, the needle will pass sequentially through P 1 , P 2 and reach P n with an increasing l. e e e Validations Several typical images in the sagittal plane comparing the simulated and the real needle images in both water and a gelatine phantom are shown in Fig. 2.14. For example the needle in Fig. 2.14d roughly corresponds to the geometry in Fig. 2.12. As we could see, the intensity and span of the simulated images match quite well with the real ones. In these figures, the needle is pointed out by the arrows; the other bright lines are needle tracks left by previous insertions. As can be seen in Fig. 2.14d, the needle track immediately behind the needle is not clearly visible in the real image, because the needle obstructs the ultrasound. Further work is needed to simulate this kind of effect. It also could be seen that the needle is always tilted a little bit upward in the phantom; this is possibly due to the torque on the transducer resulting from the excess load, as explained further in Section 2.4. 2.3 Seed image simulation in the transverse and sagittal plane In this section, the data collection and interpolation methods for simulating the brachytherapy seed images in both the transverse and the sagittal planes are described. Representative results are presented. 2.3.1 Methods Data Collection The goal is to collect B-mode images of seeds with various orientations and positions with respect to the imaging plane. There are several options: first, we can rotate the probe – this is the method used in [20]. Second, we can change the inclination of the water tank – this is the method used in [52]. Third, we can manipulate the seed orientation directly – this is the method we will use. It will be of great convenience if the robot could be used for manipulating the spatial configurations of the seed. However the robot is designed to hold 30 2.3. Seed image simulation in the transverse and sagittal plane (a) Linear interpolation (b) Nearest neighbour interpolation Figure 2.13: User interface showing different interpolation methods used for image rotation 31 2.3. Seed image simulation in the transverse and sagittal plane Sim Real (a) (0.18, 5, 0, 0, 0, 0) Sim Real (b) (0, 5, 4, 0, 0, 0) Figure 2.14: Comparison of simulated and real needle images. In (a) through (d), upper-left quadrant: simulated image in water; upper-right: real image in water; lower-left: simulated image in phantom; lower-right: real image in phantom. Configurations shown in the subfigure captions are of the form (x, y, z, ρ, α, β). 32 2.3. Seed image simulation in the transverse and sagittal plane Sim Real (c) (−0.25, 5, 0, 135, 2.5, 0); images in water contrast-adjusted for better visibility. Sim Real (d) (0, 3.3, 4, 0, 2, 8); images in water contrast-adjusted for better visibility. Figure 2.14: Comparison of simulated and real needle images. In (a) through (e), upper-left quadrant: simulated image in water; upper-right: real image in water; lower-left: simulated image in phantom; lower-right: real image in phantom. Configurations shown in the subfigure captions are of the form (x, y, z, ρ, α, β). 33 2.3. Seed image simulation in the transverse and sagittal plane Figure 2.15: Phantom with a seed mounted on a needle needles, therefore the problem is to attach the seed rigidly to the robot. After some experiments, we decided to do this by passing a needle through a small piece of phantom that encapsulates a seed. By aligning the needle with the seed and measuring the distance in between, we can control the seed’s position and orientation by moving the needle using the robot. Fig. 2.15 shows the mounting. The quality of the phantom is crucial for the task: a too rigid phantom will break under the needle load; one too soft will wobble when being moved in the water and takes a long time to settle down. In the end, a 10% gelatine concentration proves to work fine. The phantom loses water and becomes fracture-prone very quickly in the fridge, so it should be used instantly after being made. Medi-physics non-corrugated seeds were embedded in the phantom by layering: a layer of gelatine was cooled down first, then the seeds were put on top, a second layer was then poured on top to include the seeds. After cooling, small cubes of phantom were cut off and mounted on the needle. The cube should be large enough to maintain clearance between the needle and the seed, so that we can segment the seed image from the needle later on. The parameterization is the same as for the needle, except that the seed is symmetric around its axis, therefore we do not need the ρ dimension. l is changed by moving the probe; with the help of the encoder mounted inside the stepper knob, 1mm accuracy could be achieved. The range and sample numbers of each dimension are shown in Table. 2.3a and Table. 2.3b. The collected images have a resolution of 488 × 356. ® 34 2.3. Seed image simulation in the transverse and sagittal plane Variable x α β Range -10 ∼ 10mm 0 ∼ 9◦ 0 ∼ 15◦ # of samples 22 4 4 (a) Sagittal sample range and numbers Variable r l α β Range 2.5 ∼ 5.5cm 1 ∼ 17mm 0 ∼ 9◦ 0 ∼ 15◦ # of samples 7 17 4 4 (b) Transverse sample range and numbers Table 2.3: The range and numbers of the seed image samples in the sagittal and the transverse planes. (a) One end of the seed (b) Centre of the seed (c) The other end of the seed Figure 2.16: Comparison of seed images in the transverse plane 2.3.2 Results and Discussions Visibility The general trends of the seed’s visibility agree with that of the needle. A special pattern is that when imaging the two ends of the seed, the transverse images are brighter presumably due to different scattering conditions. See Fig. 2.16 for three consecutive transverse images with changing depths. The first and the third images correspond to the two ends of the seed. It can be noted from the samples collected that: 1. The seed is nearly invisible when the pitch angle α is larger than 9◦ . 2. The seed image does not change much with the change in the yaw 35 2.3. Seed image simulation in the transverse and sagittal plane Sim Real (a) Transverse simulation 1 Sim Real (b) Transverse simulation 2 Sim Real (c) Sagittal simulation Figure 2.17: Comparison of simulated and real seed images. angle β. Validation Both the transverse and the sagittal simulations of the seed images were validated by comparison to in vivo images collected in the hospital. Several examples are shown in Fig. 2.17. The seed locations were verified from the pre-plan so that other artifacts were not mistaken for seeds. The seed image was first artificially removed by the “clone tool” in GIMP as was done in the transverse validation for the needle. Then the simulated seed image was superimposed on the blank background. One of the artifacts that could sometimes be seen in the seed images is 36 2.4. Experimental characterization of the transverse shadow forming conditions Figure 2.18: Needle image showing bright tail in a scatterer-free phantom the dark shadow. Potentially this could be simulated by using the shadow simulation method described in Section 2.1.4. 2.4 Experimental characterization of the transverse shadow forming conditions 2.4.1 The bright tail One issue observed in the phantom validation experiment for the needle was that the bright tail was seen in the water but not seen in the phantom (Section 2.1.5). A dark shadow showed up instead. We therefore test several hypotheses that may explain this difference. Hypothesis 1: The elasticity of the phantom changes the needle response. We test this hypothesis by making a transparent phantom that has the same elasticity with the one used in the validation experiments. The bright tail does show up (Fig. 2.18), rejecting this hypothesis. Hypothesis 2: The scatterers cause the bright tail to disappear in the phantom. This is tested by adding celluloid scatterers in the water. Again, the bright tail appears (Fig. 2.19a), rejecting the hypothesis. Since we know from the water tank experiment that when the pitch angle is not close to zero, the bright tail is not seen, we formulate another hypothesis. 37 2.4. Experimental characterization of the transverse shadow forming conditions (a) Bright tail in water with suspended scatterers. (b) Dark shadow in water with suspended scatterers. Figure 2.19: Needle image in water with suspended scatterers. Hypothesis 3: The bright tail was not seen in the validation experiment because the pitch angle is not zero. We therefore scan a phantom made from PVC with a linear probe using freehand so that we can freely adjust the angle. The phantom has similar elasticity and scatterer distribution with the gelatine phantom used in the validation experiments. The result shows that the bright tail does appear when the probe is orthogonal to the phantom (Fig. 2.20a). This hypothesis is further confirmed by the sagittal validation experiments, which used the same set-up as the transverse validation. As can be seen in Fig. 2.14, the horizontally inserted needle appears to be tilted upwards. The reason this happened is postulated to be due to the excess torque exerted on the probe by the phantom and its casing, as sketched in Fig. 2.21. To test the influence of this factor on the outcome of the original experiment, another gelatine phantom was made and was carefully placed to avoid excessive loading. The bright tail appears as expected (see Fig. 2.22a). Phantom fracture It was also noticed during the course of the experiment that another factor might contribute to the disappearance of the bright tail in the previous 38 2.4. Experimental characterization of the transverse shadow forming conditions (a) The bright tail in a PVC phantom (b) The dark shadow in a PVC phantom Figure 2.20: Needle in a PVC phantom. Needle Phantom τ TRUS probe α Figure 2.21: The load of the phantom may exert torque on the probe and change the pitch angle. 39 2.4. Experimental characterization of the transverse shadow forming conditions (a) The needle shaft in the phantom, before fracture (b) The needle shaft in the phantom, after fracture Figure 2.22: Phantom fracture makes bright tail disappear. validation: after multiple insertions into the same track, the phantom tends to fracture, allowing air to seep into the cracks. This in fact happens very frequently in aged and dry phantoms. As a result around the needle track the acoustic impedance increases dramatically; little energy can be transmitted through the needle surface to create the bright tail. The needle shaft image after phantom fracture is shown in Fig. 2.22b; the needle has the same configuration as that in Fig. 2.22a. This mechanism is confirmed by immersing the phantom under water and repeating the insertion – the bright tail shows up even after the fracture because the air layer is not formed. As for the needle image in the tissue, similar factors determine the appearance of the bright tail, although fracture is not likely to occur. In Fig. 2.23a the bright tail in a pig liver in vitro is apparent when the ultrasound beam is orthogonal to the needle. 40 2.4. Experimental characterization of the transverse shadow forming conditions (a) The bright tail in a pig liver (b) The dark shadow in a pig liver Figure 2.23: Needle in a pig liver; the needle location is pointed out. 2.4.2 The dark shadow Provided that there are scatterers present to create contrast, when the incident angle changes from orthogonal to oblique, the bright tail gradually disappears and exposes the dark shadow distal to the needle. This can be observed both in the phantom (Fig. 2.20b) and in the tissue (Fig. 2.23b). In fact, when the scatterers are present, the dark shadow is also visible in water (Fig. 2.19b). As noted in Section 2.4.1, the dark shadow can also be formed when the phantom fractures in air. This is due to an acoustic impedance change, not due to the needle configuration. It is known that the dark shadow may be less apparent when the focal zone is not set near the object [45]. This is because when the object is not in the focal zone, the incident wave is not planar, thus may reach behind the object and insonify the distal region. However, on the machine and with the brachytherapy settings we use, this effect is not apparent. 41 2.4. Experimental characterization of the transverse shadow forming conditions 2.4.3 Conclusion The appearance of the bright tail in the phantom requires orthogonal incidence, and no fracture in the phantom. Scatterers and phantom elasticity are not significant factors for the appearance of the bright tail. The dark shadow requires scatters and an oblique incidence of the US to needle in order to appear. A high acoustic impedance difference may also cause the dark shadow. As demonstrated in the various experiments, the bright tail appearance in the water and in the phantom/tissue obeys the same rule. With the correct set-up, the bright tail pattern in the phantom does resemble that in the water (see Fig. 2.1); direct superimposition as we used in Section 2.1.4 is therefore justified. To simulate the dark shadow, on the other hand, we need to artificially add the artifact, as is done in Section 2.1.4. When looking closer, one could see that the dark shadow is actually also present when the bright tail occurs, although the former is usually outshone by the later (see Fig. 2.19a for example). An explanation for this is provided by the models developed in the next two chapters. 42 Chapter 3 Modelling of the bright band artifact in the sagittal plane From the needle insertion experiments discussed in the previous chapter, we observe that the bright tail pattern in the transverse plane is a universal feature no matter in which medium the needle/seed is embedded. Also, the dark shadow occurs under certain conditions. We also note that in the sagittal plane, the needle exhibits bright band artifacts (see Fig. 3.1). It is our goal in this chapter and next to gain a more general understanding of scattering from metallic objects; more specifically we wish to understand the bright band/tail and dark shadow forming mechanisms. To this end, physical models of acoustic scattering from two basic geometries – a plate and a sphere – are studied in detail and the results are validated by experiments. Their connections to the needle and seed scattering are also pointed out. In the studies of acoustic scattering, three general approaches are available. For scatterers much larger than the wave length, “geometric acoustics” is an easier and quicker way for modelling. In this case the acoustic wave propagation can be approximated by tracing “rays” using reflection and refraction laws, as will be demonstrated in this chapter. On the other hand, scatterers much smaller than the wavelength can also be studied by using approximations; for example, quantities such as back-scattered intensity has (ka)4 dependency [4] in this case. These ultra small scatterers are known as Rayleigh scatterers. The most complex case which lacks apparent approximations is when the scatterers are at the same scale as the wavelength, as is the case for the needle or the seed. In this case, explicit wave equations need to be used to obtain the solutions, as will be demonstrated in Chapter 4. Furthermore, we note that the first two approaches are simply approximations and simplifications to the wave equation approach. Specifically, as demonstrated in Appendix A, geometric acoustics description can be derived from the wave equations. In this chapter, we study the bright band artifact exhibited by a metal plate and also by the needle in the sagittal plane. 43 3.1. A description of the phenomenon (a) Needle filled with water (b) Needle filled with air Figure 3.1: Sagittal B-mode images of needle in water In the following as a general practice, except explicitly noted, the time part e−iωt of a time harmonic quantity C(ω)e−iωt is omitted. See Appendix A for a short discussion. 3.1 A description of the phenomenon When imaging metal plates with nearly perpendicular incident angles, repetitive secondary reflections can be observed behind the major reflection as shown in Fig. 3.2. See Fig. 3.4 for the imaging set-up. It was our previous (false) belief that because of the large impedance difference between the metal and the water, secondary reflections should be extremely weak and should quickly die out, if appearing at all. With the following analysis, it is shown that the secondary reflections are indeed much weaker, however they do not die out easily, agreeing with the experimental observations. 3.2 Modelling of the phenomenon First, we note that the beam pattern from the transducer can be assumed to be planar around the focal point [4]. In the following, we use the terms “reflection” and “scattering”, “refraction” and “transmission” interchangeably. Because the plate thickness cannot be neglected compared with the wavelength, thin plate theory [43] cannot be applied. We instead model the thick 44 3.2. Modelling of the phenomenon (a) Aluminum plate (b) Steel Figure 3.2: Bright band when imaging in B-mode an aluminum plate and a steel plate at perpendicular angle with the transducer. The double arrow indicates the incident direction; the small arrow indicates the reverberation between the transducer and the plate. 45 3.2. Modelling of the phenomenon plate as two parallel planes separating two media and use methods developed for layered media [12] to investigate the phenomenon. The two media to be studied here are fluid and solid, each governed by a set of wave propagation laws. When assuming linear propagation, the pressure in the fluid and the longitudinal and shear displacement potentials in the solid obey Helmholtz equations. With associated material law and the compatibility and equilibrium boundary conditions, we can solve these equations for pressure in the fluid and stress in the solid. Since we are more interested in the intensity of the propagating fields and do not care too much about the phase, we can denote the amplitude of the pressure and stress fields by “rays” and use ray tracing to study the wave propagation, as is done in geometrical optics. The two important physical quantities needed for ray tracing are the reflection and refraction coefficients V and W , denoting the scaling ratios to the incident wave amplitude. These ratios could be derived by solving the wave equations. In Appendix A, we derive V and W for the case of orthogonal incidence from the fluid to the solid, where V is associated with the reflected pressure amplitude and W with the forward transmitted traction in the solid. We restate the results here: p = pinc + ps = eikz + V1 e−ikz T−n = W1 eikα z where k is the wave number in the fluid, kα that in the solid, z the incident direction; p and T−n are the pressure and forward traction respectively and V = ρs α − ρf c 1 ρs α + ρf c (3.1) W =1+V 1 1 are determined by the acoustic impedance ρc of the materials. Here α is the longitudinal wave speed in the solid. It is possible to show [12], though not attempted in Appendix A, that for orthogonal incidence of longitudinal wave from the solid to the fluid, the reflection and transmission coefficients are given by V = ρf c − ρs α 2 ρs α + ρf c (3.2) W =1+V 2 2 Now we could ray-trace the orthogonal wave propagation in the media comprised of three layers – in our case fluid–solid–fluid; see Fig. 3.3 for the first few reflections/transmissions. 46 3.3. Validation of the model fluid solid W1 W2 W1 W 1 V2 W1 V22 fluid 1 V1 W 1 V2 W 2 Figure 3.3: Ray tracing of the reflections and transmissions; the lengths of the rays are roughly proportional to the amplitudes of the pressure/traction when the solid is aluminum. The time axis is from the left to the right. There are several things worth noticing here: 1. Since mode conversion does not occur (see Appendix A), only the longitudinal wave is present in the solid, thus we could use Eq. 3.2. Also the speed of sound in the solid is that of the longitudinal sound speed. 2. Since the acoustic impedance of the metal is larger than that of the water, V2 is negative, which means that the phase of the reflection from the solid-water interface is inverted. This is analogous to the reflection of wave on a string from the support (see [56] chapter 4). 3. W1 > 1, which means the pressure amplitude is “amplified” in the solid. The quantity that is conserved is the energy, not the pressure: the power flux (sound intensity) is related to the pressure amplitude and the acoustic impedance as follows: I= |p|2 ρc (3.3) (See [56] chapter 6 for the assumptions for this relation), which means that for the same intensity, the solid with higher impedance will have higher pressure. 3.3 Validation of the model To validate the reflection/transmission model in the previous section, metal plates of various thickness and materials are scanned in the water using a 47 3.3. Validation of the model Probe Plate being scanned Support Figure 3.4: Set-up for metal plate reflection study linear ultrasound transducer held perpendicular to the plate; the set-up is shown in Fig. 3.4. As examples, we examine the results from a 6.33 mm thick aluminum plate and a 19.9 mm thick steel plate. In all these experiments, Time Gain Compensation (TGC) is disabled since we are scanning in the water; saturation of the RF signal is avoided by checking that the maximum value in the signal is much less than the resolution of the receiver (215 in our case). Example B-mode images are shown in Fig. 3.2; example raw RF signals (pre-envelope-detected) are shown in Fig. 3.5. 3.3.1 Sound speed As indicated in the previous theoretical discussion, the sound speed in the plate should be that of the longitudinal wave. In order to validate this, the thickness of the plate is measured with calliper and the time interval between the band is estimated by taking the average length over several band. The sound speed in the material can be calculated from the elasticity properties (see Appendix A). Here we quote the values from Onda [60]: the longitudinal wave speed in the aluminum is around 6380 ∼ 6420 m/s and the shear wave speed around 3040 m/s; from the RF line, we estimate the sound speed to be 6330m/s – agreeing with the longitudinal wave value. Similar results are obtained for the steel plate. 48 3.3. Validation of the model 4 1.5 RF: aluminum plate x 10 1 Intensity 0.5 0 −0.5 −1 −1.5 0 500 1000 1500 2000 2500 3000 3500 2500 3000 3500 t (a) Aluminum plate 4 2.5 RF: steel plate x 10 2 1.5 1 Intensity 0.5 0 −0.5 −1 −1.5 −2 −2.5 0 500 1000 1500 2000 t (b) Steel plate Figure 3.5: RF lines corresponding to the B-mode images in Fig. 3.2; the position of the lines are both at around 3/4 of the B-mode image. 49 3.3. Validation of the model 3.3.2 Band intensity The piezoelectric material used in the ultrasound transducer has the following pressure-voltage relationship [17]: pi = eii Ei − cii where p is the pressure in the material, eii is the piezoelectric stress coefficient, Ei is the electric field, cii is the elastic stiffness constant, is the strain in the material, and i denotes direction. We see that the relationship is roughly linear. It is therefore our first guess that the amplitude of the RF line is linearly proportional to the reflection pressure with amplitudes ps1 = V1 , ps2 = W1 V1 W2 , and so on, as in Fig. 3.3: Hypothesis 1: Z Z ps δ(t)dt ∝ |ARF |dt (3.4) We test this hypothesis by integrating the absolute values of the RF signal around each peak, so that we could include all the transient response. We normalize the two sides of Eq. 3.4 so that the amplitudes for the first band are both 1. The results are shown in Fig. 3.6. To test the sensitivity of the experimental results to various factors, we vary the incident power to make sure that the piezoelectric material in the transducer works in the linear region. We also vary the frequency to test the frequency-dependency. Both shows consistently similar relative amplitudes in the experiment, which do not agree with that of the theory. We therefore make a second hypothesis, that the amplitude of the RF signal is proportional to the sound intensity. According to Eq. 3.3, we could write this as: Hypothesis 2: Z Z |ps | δ(t)dt ∝ 2 |ARF |dt (3.5) Again we test it by plotting the left and right sides of the equation for each band and compare. The results are shown in Fig. 3.7, which shows very good agreement between the theoretical and experimental results. Hypothesis 2 is thus supported by the experimental evidence. 50 3.3. Validation of the model Comparison of results: reflections from an aluminum plate 1 theory experiment: power 2, 6.6MHz experiment: power 9, 6.6 MHz experiment: power 15, 6.6MHz experiment: power 15, 10MHz expermient: power 15, 5MHz relative integrated intensity in each band 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 5 6 number of reverbs 7 8 9 10 (a) Aluminum plate Comparison of results: reflections from a steel plate 1 theory experiment relative integrated intensity in each band 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 5 6 number of reverbs 7 8 9 10 (b) Steel plate Figure 3.6: Comparison of the theoretical and experimental results when we assume the RF amplitude is proportional to the pressure. 51 3.3. Validation of the model Comparison of results: reflections from an aluminum plate 1 theory experiment: power 2, 6.6MHz experiment: power 9, 6.6 MHz experiment: power 15, 6.6MHz experiment: power 15, 10MHz expermient: power 15, 5MHz relative integrated intensity in each band 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 5 6 number of reverbs 7 8 9 10 (a) Aluminum plate Comparison of results: reflections from a steel plate 1 theory experiment relative integrated intensity in each band 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 5 6 number of reverbs 7 8 9 10 (b) Steel plate Figure 3.7: Comparison of the theoretical and experimental results when we assume the RF amplitude is proportional to the sound intensity. 52 3.4. Discussion Figure 3.8: Band pattern in the quarter-wave layer when scanning air. 3.4 3.4.1 Discussion Incident direction It is observed that the band pattern is strongest and clearest at one single orientation of the transducer; according to the ray tracing model and the transducer geometry, we expect this orientation to indicate the orthogonal direction of incidence. However, in the experiment carried out, the strongest reflection appears when the metal plate is tilted slightly by a support as in Fig. 3.4. This possibly indicates that the incident direction of the linear transducer used is tilted. Also, if it is indeed at the perfect orthogonal angle that the multiple band pattern occurs, this pattern may potentially be useful for calibrating the probe position with respect to a target. 3.4.2 Quarter-wave layer It is interesting to note that the quarter-wave layer in the transducer used for impedance matching with the tissue exhibits the same band pattern when the transducer is exposed in air (Fig. 3.8). This is easily understood by noting that the air has a much smaller impedance than the water/tissue. It is also interesting to note that for the focal point setting in the figure, the band are not horizontally homogeneous. 3.4.3 Surface roughness Note that for rough surface materials, the band pattern is not seen. Intuitively, this corresponds to the “diffuse reflection” scenario in geometric optics, whereas the one considered above corresponds to the “specular re- 53 3.4. Discussion flection” (and refraction). A quantitative study of scattering from surface roughness can be seen in [56]. 3.4.4 Relation to the needle image As shown in Fig. 3.1, a similar band pattern can be observed when imaging the needle in the sagittal view. In this case, there are 5 layers of media. Our analysis for the plate can be used to help understand the artifact. When the middle layer is air instead of water, less energy is transmitted to the second needle surface due to a larger acoustic impedance difference, leaving weaker band behind. It should be noted that for the needle in the sagittal view the band pattern is usually less pronounced than for the plate, especially in the tissue or in the phantom. This is presumably because most of the energy is directed away from the transducer due to the curvature of the needle surface. In the transverse plane, on the other hand, the bright tail is more apparent because the energy directed away by the curvature is also received by the transducer, especially when the needle is perpendicular to the imaging plane, giving a brighter tail. 54 Chapter 4 Modelling of the bright tail and dark shadow artifacts in the transverse plane 4.1 Introduction When imaging the needle in the transverse plane, both the bright tail and the dark shadow may occur distally to the needle body (Section 2.4). In most cases, the bright tail of the needle image is made up of tightly spaced band which are reminiscent of the band pattern discussed in Chapter 3. For example, the band can be seen in water, phantoms and tissue. See Fig. 4.11 . This distinct bright tail was also observed from small metal spheres such as lead pellets ([82], also Fig. 1.4) and collections of air bubbles [5] in medical ultrasound. We thus seek to model the tail forming mechanism by studying scattering from metal spheres or cylinders. As is noted in Section 1.2.7, scattering from spheres and cylinders is a classical problem and can be modelled by forming the appropriate boundary value problem; the formulation is studied further in Section 4.3. The choice to model the scattering from the sphere in this chapter is largely due to the availability of the simulator [7], in which the displacement potentials can be calculated as a function of the observing position and incident frequency. In the following we develop methods to simulate both the frequency and transient responses of the metal scatterer based on the displacement potential coefficients (c’s in Eq. 4.4) calculated by the simulator. It should be stressed that the basic assumption made in [7] and inherited here is that the media and the scatterer are both made of solid materials and 1 The band pattern is not always seen in the B-mode images of the needle in the gelatine phantom, however it is still present in the RF signal and is presumably obscured by the digital processing of the B-mode image. 55 4.2. Hypotheses for the bright tail and dark shadow forming mechanism (a) The needle shaft in water (b) The bright tail in a PVC phantom (c) The bright tail in a pig liver Figure 4.1: Band pattern visible in water, phantom and liver. therefore can support both longitudinal and shear waves. This is useful to study scatterers in the tissue and phantom, but results in a different formulation from the rest of the underwater acoustics literature (as those reviewed in Section 1.2.7), which assume the media to be fluid. See Section 4.3 for more details on the assumptions. 4.2 Hypotheses for the bright tail and dark shadow forming mechanism Bright tail: We postulate that the bright tail is formed from the backscattered pulse response of the sphere. Because of the elasticity of the sphere, acoustic energy enters the sphere and is gradually emitted back to the transducer, forming long tails. Dark shadow: We postulate that the dark shadow appears because of the phase cancellation (destructive interference) between the incident wave and the forward scattered wave when the scatterer is sufficiently large compared with the wave length. This is stated as a fact in [56]. This effect is amplified in the B-mode imaging because the back-scattered wave from the Rayleigh scatterers in the tissue distal to the metal 56 4.3. Boundary value problem for scattering from spherical shells r θ Φ(2) , Ψ(2) Φ(1) , Ψ(1) (0) (0) Φs , Ψs (0) Φinc Figure 4.2: Scattering by a spherical shell. For simplicity,only three layers are drawn here. All layers are elastic solids. scatterer has to transmit past the metal scatterer once more to reach the transducer. We test these hypotheses by simulating the scattered and incident waves and observe that they indeed yield the wave forms that resemble the bright tail and dark shadow. 4.3 Boundary value problem for scattering from spherical shells The formulation and the solving of the boundary value problem is a standard technique to obtain the frequency and time domain characteristics of the scattered wave from spherical solid and shells. The references in Table 1.1 all use this formulation in one form or another. Here, we follow the assumptions in [7] and give a sketch of the procedure. The major difference between this formulation and other previous works is that all layers are considered to be solid and therefore allowing both longitudinal and shear waves. The geometry of the problem is shown in Fig. 4.2. As shown in section A.1, under the time harmonic assumption, with linear elasticity, the governing dynamic equations for the elastic solid are the Helmholtz wave 57 4.3. Boundary value problem for scattering from spherical shells equations for displacement potentials: Eq. A.2. Repeating them here: ∇ · ∇Φ + kα2 Φ = 0 ~ =0 ~ + k2 Ψ ∇ · ∇Ψ (4.1) β where kα = ωα , kβ = ωβ are the longitudinal and shear wave numbers. Φ and ~ are respectively the longitudinal and shear components of the displaceΨ ~ (the gradient of U ~ is the displacement.): ment potential U ~ = ∇Φ + ∇ × Ψ ~ U (4.2) See Appendix A for the full details. To solve the wave equations 4.1 in spherical coordinates, we separate variables, that is, we write each of the equations in Eq. 4.1 as 3 separate differential equations with variables r, θ, φ respectively. We can solve these separated equations by using classical techniques. 2 In spherical coordinates, when we assume the incident wave to be planar and to propagate along the θ = 0 direction, the axial symmetry with respect to the azimuth angle φ makes one degree of freedom vanish. Eq. 4.2 becomes: ~ = ∇Φ + ∇ × ( ∂Ψ êφ ) U (4.3) ∂θ where Ψ is now a scalar field. By separation of variables, we can write the general solutions for layer m as: X Φ(m) = [j cΦ,(m) jn (kα(m) r) + n cnΦ,(m) nn (kα(m) r)]Pn (cos θ) n n Ψ (m) = X (m) (m) [j cΨ,(m) jn (kβ r) + n cnΨ,(m) nn (kβ r)]Pn (cos θ) n (4.4) n where jn (·) is the nth mode of the spherical Bessel function, nn (·) is that of the spherical Neumann function, Pn (·) is the Legendre polynomial of order n, m indicates the layer. We call the coefficients c the displacement potential coefficients; they are the unknowns we wish to solve for in the end. For different layers and different propagating directions, different subspaces of these general solutions are allowed. We analyze this layer by layer: 2 The separated equations can be seen in [56] chapter 7. Before writing out the solution, we can make some general observations: the separated ordinary differential equations have the form of the eigen problem. Therefore the solutions are eigen functions. They are also orthogonal basis functions that span the solution space, analogous to the case of the eigen vectors being orthogonal basis when the linear system is symmetric [72]. 58 4.3. Boundary value problem for scattering from spherical shells 1. For the outer-most layer, we separate the incident wave and the scattered wave. We assume the incident displacement potential to have unit amplitude and to be irrotational (no Ψ). We see these assumptions are satisfied when Φ has the form of Eq. A.6: (0) (0) Φinc = eikα z We can expand it in spherical coordinates [56]: (0) (0) Φinc = eikα r cos θ = X (2n + 1)in jn (kα(0) r)Pn (cos θ) (4.5) n Φ,(0) Comparing to Eq. 4.4, it can be seen j cn = (2n + 1)in , and the (0) Neumann function does not contribute to this solution. Ψinc = 0 by definition. The scattered solution is to the outward direction. Only a certain combination of the Bessel and Neumann functions, i.e. spherical Hankel function of the first kind, can meet this requirement [39]: h[1] n (kr) := jn (kr) + inn (kr) Φ,(0) In this case, we can use a single coefficient h cn (0) for Ψs . (0) Ψ,(0) for Φs , and h cn 2. For the middle layers, both the inward and outward solutions are allowed. We know that the spherical Hankel function of the second kind: h[2] n (kr) := jn (kr) − inn (kr) represents the inward waves. Therefore the solution could be any linear combinations of the Bessel functions of the first and second kinds; we need all two degrees of freedom of the original solution space for both Φ(1) and Ψ(1) . 3. For the inner-most layer, only the spherical Bessel function is allowed as the basis function. This is because nn (kr)|r=0 = ∞ which should not be allowed. 59 4.4. Scattering spectrum simulation (0) (0) We therefore have 2 unknowns in Φs , Ψs for the outermost layer, 4 unknowns for each of the middle layers, and 2 unknowns for the innermost layer. On average, we have 4 unknowns for each boundary. We thus need 4 boundary conditions for each boundary to solve the system. These boundary conditions are: 1. The normal and tangential displacements must be continuous across the boundaries for each mode: (i) (i+1) Ur,n = Ur,n (i) (i+1) Uθ,n = Uθ,n (4.6) 2. The normal stress must be continuous across the boundaries for each mode: (i) (i+1) = σrr,n σrr,n (i) (i+1) σrθ,n = σrθ,n (4.7) Note that although physically we only require the continuity conditions to be met for the “lumped” variables Ur , Uθ , etc, we require each mode (projection to the basis function) to be continuous. This is because modes are orthogonal and the projections need to be unique. Also note that these boundary conditions are different than those obtained from the fluid-loaded surface (e.g. Eq. A.4 and A.5). Using Eq. 4.3, we can relate Ur,n , Uθ,n to the coefficients in Eq. 4.4 and form Eq. 4.6. We can further use the linear stress-strain relations in spherical coordinates to form Eq. 4.7. The detailed forms of these relations can be found in [7]. After solving the linear system, we obtain the displacement potential coefficients in Eq. 4.4. Thereafter we can calculate the displacement, the strain, and the stress using the same relations used to construct the linear systems. 4.4 Scattering spectrum simulation The physical quantity of interest, based on the discussions in Section 3.3.2, is the acoustic intensity, which is proportional to the pressure squared. It is unknown, however, how the sign is assigned in the RF signal. Therefore, in the following, we consider only pressure or traction, with the knowledge that the real RF signal may have larger fluctuation due to squaring. 60 4.4. Scattering spectrum simulation Referring to Fig. 4.2, the directions we are interested in are θ = π and θ = 0, that is, the back scattering and forward scattering directions. These correspond to the cases in the two hypotheses in Section 4.2. It can be shown, using the procedures developed in Section 4.3, that (0) σrr,n,scattered (t) – the scattered normal stress component in the outermost layer (the tissue or phantom) can be written in terms of the displacement potential coefficients defined in Eq. 4.4 [7]. Omitting the subscript “scattered” and the superscript (0): σrr,n (ω0 t) =σrr,n e−iω0 t Φ Φ = h cΦ n λT1 + 2µT2 + h cΨ n Ψ Ψ λT1 + 2µT2 Pn (cos θ) e−iω0 t (4.8) where T1Φ = − kα2 h[1] n (kα r) n(n − 1) 2kα [1] T2Φ =( − kα2 ) h[1] hn+1 (kα r) n (kα r) + 2 r r (n − 1)n(n + 1) [1] hn (kβ r) T1Ψ = − 2 r2 kβ n(n + 1) [1] (n − 1)n(n + 1) [1] T2Ψ = − hn (kβ r) + hn+1 (kβ r) 2 r r h Ψ The coefficients h cΦ n , cn , as shown in Section 4.3, can be solved from a system of linear equations 4.6, 4.7 obtained from the boundary conditions. (0) (0) Furthermore, the shear stress components σrθ , σrφ can be shown to be proportional to µ, which is very small compared to λ in the phantom (Table. 4.1). Therefore we can assume that in the backward and forward directions, the traction is proportional only to the normal stress σrr 3 . We now have the normal traction as a function of observing position (r, θ), and the incident frequency ω0 , to which the wave number kα , kβ are proportional. In other words, we have a steady-state representation of the traction, corresponding to the response of the system under continuous single-frequency interrogation. Now with the input being the unit incident displacement potential (Eq. 4.5) and the output being the normal traction at the observing point (Eq. 4.8), 3 This is verified by evaluating σrθ and comparing with σrr – the former indeed very small. Legendre polynomials and derivatives are calculated using the “lpn” routine converted from FORTRAN [10] 61 4.4. Scattering spectrum simulation we can use the Fourier transform to reconstruct the spectrum of the linear system response4 . First we note that for a function f (x, t) = f (x)e−iω0 t (with separated spatial and time parts), the Fourier transform is: F{f (x)e−iω0 t } = f (x)2πδ(ω + ω0 ) (4.9) that is, the signal contains only one frequency, with the spatial part as the amplitude. Without loss of generality, we measure the input at z = 0: Φinc (ω0 t) = eikz−iω0 t = e−iω0 t Applying Eq. 4.9: F{Φinc (ω0 t)} = 2πδ(ω + ω0 ) Similarly for the output (Eq. 4.8): F{σrr,n (ω0 t)} = 2π σrr,n δ(ω + ω0 ) The frequency response of the system is therefore H(−jω0 ) = σrr,n (jω0 ) That is, we get the negative part of the frequency response from the positive frequency part of the stress. To obtain the positive part of the frequency response, we note that to be physically meaningful, the stress needs to be real, therefore its Fourier transform needs to be conjugate symmetric [61]: ∗ H(jω0 ) = H ∗ (−jω0 ) = σrr,n (jω0 ) where ∗ denotes complex conjugate. One example of the resulting spectrum of the back-scattered stress from a solid steel sphere – that is, not hollow inside – can be seen in Fig. 4.3. We simulate the solid sphere first since it is the simplest case; the scattering from shell will be simulated in Section 4.5.4. The frequency range we will be studying is from 0 to 6 MHz, corresponding roughly to the frequency range of the ultrasound transducer (see Section 4.5.1 for more details). When studying scattering, the characteristic quantity is ka – the argument to the eigen functions in Eq. 4.4. By noting that ka = 2πa/λ, we see that it is the ratio between the wave length and the dimension of the scatterer that characterizes the scattering. When a = 1 mm, ka ranges from 0 to around 25 and is marked on the top of Fig. 4.3. 4 We can prove the linearity by noting that the inputs are on the right side of the linear system that we construct in Section 4.3, thus the solutions can be written in the form of the linear combinations of the inputs. 62 4.4. Scattering spectrum simulation 16 Amplitude 2 5 x 10 Magnitude, θ=180, r =60a 10 15 ka 20 25 2 1.5 1.5 1 1 0.5 0.5 0 5 x 10 0 6 5 6 16 0 1 2 3 4 Frequency (MHz) Phase, θ=180, r =60a Radians 0 −100 −200 scattered stress −300 0 1 2 3 Frequency (MHz) 4 Figure 4.3: Simulated spectrum of the back-scattered stress from a radius a = 1 mm solid sphere measured at a distance r = 60 mm from the sphere, 200 sample points. 4.4.1 Discussion Several features of the spectrum are worth noticing: 1. The higher frequency generally has larger magnitude. This agrees with the fact in the fluid: the larger the ka, the greater the back-scattered intensity [56], when ka is small. The exact form of the intensity growth, however, does not agree with that of the scatterers in fluid. In underwater acoustics, the quantity that measures the back-scattered intensity is the form function, which is proportional to the pressure. Referring to Fig. 4.4, when ka is small (< 1 generally), the form function amplitude roughly follows (ka)2 – this is the region of Rayleigh scattering [4]. When ka is larger, the form function amplitude fluctuates. Comparing Fig. 4.3 and Fig. 4.4, we see that the simulated traction assuming rigid medium has significantly smaller low frequency components5 . 5 Curiously, the back-scattered longitudinal displacement potential Φ exhibits similar spectrum to Fig. 4.4. Φ could not, however, be linearly related to the traction in solid. 63 4.5. Transient simulation ka Figure 4.4: Form function of a steel spherical scatterer in water (reprinted from [19] chapter 33). To the left of the added dash line is the Rayleigh scattering region. Properties ρ (kg/m3 ) λ µ Phantom 1000 2.25e9 3e4 Steel 7700 9.695e10 7.617e10 Table 4.1: Material properties used in the simulation 2. Both Fig. 4.3 and Fig. 4.4 exhibit isolated dips, which according to the previous studies correspond to the resonant frequency of the surface waves. These resonances destructively interfere and cause nulls in the spectrum (see [19] chapter 4 for an example). 3. The phase of the simulated spectrum is linear. The system response therefore can be completely characterized by the magnitude modulation and a phase shift [61]. As a side note, the elasticity moduli of the tissue/phantom used in the simulation can be estimated as follows: the first Lamé constant λ of the phantom can be calculated from the sound speed in the phantom (around 1500 m/s) and the density of the phantom (around 1 kg/m3 , assuming small shear modulus). We also assume that the Young’s modulus E of the phantom to be somewhere between 1 to 150 kPa [31]. From λ, E, the second Lamé constant µ can be calculated. Table 4.1 summarizes the material properties we use. 4.5 Transient simulation With the knowledge of the system response spectrum as shown in Fig. 4.3, we can conveniently reconstruct the impulse response of the system by taking the inverse Fourier transform. 64 4.5. Transient simulation Real part of the impulse respone, θ=180, r =60a 15 4 x 10 Scaled traction (N/m2) 3 2 1 0 −1 −2 −3 30 35 40 45 50 Time (µ s) 55 60 65 70 Figure 4.5: Simulated impulse response corresponding to the back-scattered stress from a radius a = 1 mm solid sphere measured at a distance r = 60 cm from the sphere. The y-axis has the unit of traction. In discrete time, since we expect a real signal, the 2-sided spectrum (with the negative frequency components) should be conjugate symmetric as noted in the previous section. For the spectrum in Fig. 4.3, we construct the 2-sided spectrum by mirroring the conjugate of the positive spectrum. In this case, the Nyquist rate Fs is 12 MHz; the sampling period in the time domain ∆t = 1/Fs = 0.08 µs. The length of the time axis is related to the sampling number by lt = ∆t × Ns = 33 µs, where Ns = 400 in the 2-sided spectrum6 (as noted in Fig. 4.3, the 1-sided spectrum has 200 sample points). With these parameters, the impulse response of the system in the time domain is obtained using the inverse DFT and is plotted in Fig. 4.5. 4.5.1 Modelling of the ultrasound pulse The real ultrasound pulse is not an ideal impulse, but rather an exponentially decaying sinusoid resulted from the damping of the transducer [17]. We model the pulse displacement potential as: z z Φ(z, t) = sin(ω(t − ))e−α(t− c ) c 6 This lt corresponds to roughly 25 mm in the B-mode. 65 4.5. Transient simulation Decaying sinusoid 0.6 0.4 Displacement 0.2 0 −0.2 −0.4 −0.6 −0.8 0 0.5 1 1.5 2 Time (µ s) 2.5 3 3.5 (a) The simulated ultrasound pulse. Spectrum density 7 6 Power 5 4 3 2 1 0 0 1 2 3 Hz 4 5 6 6 x 10 (b) The power spectrum of the simulated pulse Figure 4.6: The simulated US pulse and power spectrum. where c is the speed of sound in the medium. Neglecting a constant scaling factor and phase change, we can model the displacement x(t) = ∂Φ(z,t) at ∂z z = 0 as: x(t) = sin(ωt)e−αt When the frequency of the sinusoid is taken to be 5 MHz – the typical excitation frequency of the transducer, and α = 0.2, the resulting time domain signal and its power spectrum are seen in Fig. 4.6. 66 4.5. Transient simulation Figure 4.7: The experimentally measured power spectrum of the received RF signal (reprinted from [18]). 4.5.2 Validation of US pulse modelling The power spectrum of the received RF signal of the Ultrasonix machine being used was previously experimentally measured in [18], as shown in Fig. 4.7. We assume this received spectrum does not deviate much from the transmitted power distribution, so that we could use it to validate our simulation. By comparing Fig. 4.7 and Fig. 4.6b, we see the spectra resemble each other. 4.5.3 Pulse response By multiplying the previously simulated system response (Fig. 4.3) with the frequency components of the simulated US pulse, we obtain the frequency spectrum of the output stress. Again using inverse DFT, we obtain the time domain pulse response in Fig. 4.8. Comparing with the impulse response in Fig. 4.5, we see the two are not very different except that the pulse response has a little bit more “blurred” initial response. The residue after the first pulse in Fig. 4.8 could serve as an interpretation of the bright tail artifact discussed earlier. Also as noted in Section 4.1, tightly spaced band can usually be observed in the bright tail. In this simulation we could indeed see some band although they are not obvious. 4.5.4 Simulation of transient back-scattering from a spherical shell As demonstrated in Section 4.3, the scattering from a spherical shell could be formulated as a boundary value problem. The solid sphere case we have been working on is simply a special case of the shell, with the radius of the 67 4.6. Validations and interpretations of the transient response Real part of the pulse respone, θ=180, r =60a 15 4 x 10 3 Scaled traction (N/m2) 2 1 0 −1 −2 −3 −4 30 35 40 45 50 Time (µ s) 55 60 65 70 Figure 4.8: The simulated stress resulting from the simulated US pulse with the otherwise same parameters as in Fig. 4.5 inner-most layer being 0. When this value is nonzero, we are simulating a true shell. An example transient response from a shell is shown in Fig. 4.9. Comparing with Fig. 4.5, we could see that the thin shell exhibits more obvious band pattern. 4.6 Validations and interpretations of the transient response The quantitative validation of the values of the displacement potentials, which are used as the bases for our simulation, can be seen in [7]. Since there are no known previous experimental studies on the back-scattered stress in the tissue/phantom from spheres, the goal of the current validation is only to demonstrate that the transient response gives qualitatively correct results. 4.6.1 Wave speed The back scattered pulse should propagate at around 1500 m/s in the tissue. This could easily be checked by noting in Fig. 4.9 for example, the first pulse appears at around 40 µs, which corresponds to 6cm away from the origin. 68 4.6. Validations and interpretations of the transient response Real part of the pulse respone, θ=180, r =60a 15 6 x 10 Scaled traction (N/m2) 4 2 0 −2 −4 −6 30 35 40 45 50 Time (µ s) 55 60 65 70 Figure 4.9: The transient pulse response from a spherical shell with a = 1 mm, shell thickness 0.02mm, measured at r = 60 mm. 4.6.2 Resemblance to response from needle shaft in the transverse plane The needle shaft is cylindrical in shape. It seems that the best model for it is the infinite cylindrical shell, which is the subject of several previous studies (see section 1.2.7). From these studies it is noted that the responses from a sphere and an infinite cylinder share many common features. This is confirmed by comparing Fig. 4.10, the experimentally recorded needle response, and Fig. 4.9. Both exhibit band patterns in the gradually decaying tail. In the following, we further demonstrate that the spherical simulator could display wave forms that resemble the needle B-mode images . 4.6.3 Dark shadow from phase cancellation As is noted in section 2.4.2, the dark shadow may be observed distally to the body of the needle image. In section 4.2, the dark shadow is hypothesized to be a result of phase cancellation between the incident and the forward scattered wave. Here through simulation, we demonstrate that this is indeed the case when the scatterer size is large: in Fig. 4.11, we sum the near-field incident and forward scattered response from a r = 5 mm sphere. The sum appears 69 4.6. Validations and interpretations of the transient response RF: needle shaft 8000 6000 4000 Intensity 2000 0 −2000 −4000 −6000 −8000 1000 1500 2000 t 2500 3000 Figure 4.10: The RF line of the back-scattering from the needle shaft. The length of the line is around 50 µs. The later half of the response is contaminated by the reverberation between the needle and the transducer. to be much smaller than both the incident and the scattered wave. This can be verified to be due to the π radians phase difference in the scattered and the incident wave. This means that the forward scattering direction is insonified by less acoustic energy than the surrounding media, leaving the region appearing to be dark in the B-mode. We only simulate the near field here because the incident series only converges when r is small, as demonstrated in section 4.7.1. 4.6.4 Solid vs. shell As was noted in the simulation (Fig. 4.9 and Fig. 4.8), the shell creates more pronounced band. This can be confirmed experimentally by comparing the shaft image and the stylet image, the former shell and the latter solid. See Fig. 4.12. 70 4.6. Validations and interpretations of the transient response Real part of the pulse respone, θ=0, r =1a 17 5 x 10 scattered incident sum 4 Scaled traction (N/m2) 3 2 1 0 −1 −2 −3 −4 −5 0 5 10 15 20 Time (µ s) 25 30 35 Figure 4.11: The forward incident, forward scattered response, and their sum for a solid sphere with a = 5 mm, at r = 5 mm. The maximums are pointed out. (a) The shaft (shell). (b) The stylet (solid). Figure 4.12: The comparison between a cylindrical shell and a solid. 71 4.7. Numerical issues Convergence of the incident series 3 r=2mm r=4mm r=6mm r=8mm 2.5 Amplitude 2 1.5 1 0.5 0 0 50 100 150 200 250 Number of terms 300 350 400 Figure 4.13: The convergence of the incident wave series with varying r. The terms where the series start to converge are pointed out by arrows. 4.7 4.7.1 Numerical issues Convergence of the incident wave series The incident wave series in Eq. 4.5 requires more terms to converge with the increasing r, the distance from the origin. In Fig. 4.13, we could see for instance that with r = 4 mm, 120 terms are required for the series to converge, while for r = 6 mm 180 terms are needed. On the other hand, as we will see in the next section, the linear system of displacement potential coefficients is ill-conditioned with large n. This is why we only simulate the near field of the dark shadow in section 4.6.3. 4.7.2 Condition number of the linear system There is a numerical issue inherent in the boundary value problem for scattering from sphere and cylinder when ka is small. The symptom is that when ka is small (lower frequency and/or smaller scatterer), and the mode number n is large, the linear system constructed to solve for the displacement potential coefficients is ill-conditioned. Referring to Fig. 4.14, the condition number increases with n for each frequency and reaches 1016 typically when n > 10. Since 1016 is the machine epsilon for the double-precision floating point numbers, we should suspect that the linear 72 4.7. Numerical issues 272 4 30 27 2 6 30 4 0 17 8 20 27 11 48 16 80 2 20 5 24 10 240 208 208 112 144 48 ka 16 80 15 24 0 20 176 144 112 80 48 16 25 2 4 6 14 17 240 208 272 40 304 60 80 100 120 140 n Figure 4.14: The contour of log10 of the condition numbers of the linear systems of displacement potential coefficients for each mode n and each frequency. systems with larger n are in fact singular, thus cannot be solved in the usual sense. The reason behind this symptom lies in the behaviour of the Bessel function of first kind and second kind with small arguments. As an example, consider one of the boundary conditions used to solve for the coefficients (see Eq. 4.7): (0) (0) (1) σrr,n,scattered − σrr,n = −σrr,n,incident (4.10) (0) As shown in Eq. 4.8, σrr,n,scattered can be written as a linear combination of the Hankel functions of the first kind, which contains both the Bessel functions of the first and second kind. Without loss of generality, consider (1) the solid sphere (only 1 layer); σrr,n in this case contains only the Bessel functions of the first kind, as discussed in Section 4.3. Similarly, the right hand side of Eq. 4.10 contains only the Bessel function of the first kind. Now, what Eq. 4.10 asks for are 4 coefficients (in the case of 1 layer) that will match the left hand side (Hankel + Bessel first) with the right (Bessel first). From Fig. 4.15, we see that this is problematic when ka is small and n is large, because then the Neumann function (thus the Hankel function) is much much larger than the Bessel function; to linearly combine the left hand side to match the right hand side thus requires the coefficients 73 4.7. Numerical issues ka =5 50 10 Bessel first kind Bessel second kind (Neumann) 0 10 −50 10 0 10 20 30 40 50 n Figure 4.15: The Bessel functions of the first and second kind at ka = 5: jn (5) and nn (5). of the Hankel functions to be always close to 0; the coefficients solved for the Bessel part of the left hand side, on the other hand, is unlikely to be precise. To eliminate the possibility of numerical errors in the calculation of the series shown in Fig. 4.15, we compare the results with those obtained from the Matlab version of the FORTRAN “cjylv” routine [10] which is specially designed for large n. The results are similar7 . In the matrix terms, the left hand side of Eq. 4.10 has two very large columns and two very small columns in the case of small ka and relatively large n. The solutions are therefore close-to-zero coefficients for the Hankel part (the scattering part) of the equation and indeterminable coefficients for the Bessel part (the transmitting part) . This numerical issue is therefore unlikely to influence the credibility of the scattering solutions – the scattering coefficients converges to 0 very fast with increasing n and results in nicely converging series of displacement potentials; also the summation in Eq. 4.8 converges after a few terms. Thus the solution method in [7] should be valid for scattering, though there might be issues with the transmitted waves, which is not our present concern. We also note that in the previous literature reviewed in Table 1.1, this 7 This can also be verified by using the asymptotic form of the Neumann function for very small argument [77], which indicates the function diverges with n. 74 4.7. Numerical issues ill-conditioning problem may also occur. For example, in [6], the last column of the matrix formed from the boundary conditions (Eq. (6) of the paper) contains only the Bessel function of the first kind, while the second to last column contains the Neumann function. This will cause problems when ka is reasonably small and n is reasonably large as discussed above. However, this issue was not seen discussed in previous studies. In earlier studies on the solid scatterers such as in [23], the numerical problem might not exist: the solution to the boundary value problem was written out in the closed form series summation containing arctan of Bessel and Neumann functions – the possible infinity was dealt with automatically. 4.7.3 Conclusion In this section, we have demonstrated the simulation of steady-state and transient response from spherical scatterers under ultrasound pulse interrogation. The results indicate that the simulation can replicate and provide models for the bright tail and dark shadow artifacts that we see in the needle images. It is necessary to reiterate the difference between this study and the previous studies: the media the scatterer is embedded in is tissue or phantom, thus the physical quantity we have considered is stress, not pressure. The transient simulation here is derived from the first principles without considering previous studies on the similar topic in water; a comparison may be needed to see the differences. 75 Chapter 5 Implementation of the ultrasound simulator for interactive needle insertion Combined with the works in [27, 28], the needle image simulation methods developed in Chapter 2 can be used in an ultrasound simulator for needle insertion where both the probe position and the needle position are freely adjustable, and the tissue deformation visible. To achieve this, the needle spatial configuration needs to be represented in the probe’s frame and the crossing points calculated. In the following, we break the needle-probe-tissue interaction problem into several subproblems. The resultant procedures of the needle insertion simulation is shown in Fig. 5.1. The image database is prepared off-line and stored as a 1D array. The needle configuration is calculated from geometry and is used as the input to the interpolation method described in Chapter 2. After the interpolation, spatial transforms are executed to bring the needle image to the correct location and the image is superimposed on the (deformed) tissue background. The synthesized image is displayed in real time. The original algorithms were first implemented in Matlab and subsequently converted to C and re-implemented in OpenGL for real-time performance8 . 5.1 Notations and overview Some notation conventions are as follows: we measure every frame with respect to a base frame {O0 , C 0 }, where O0 is the origin and C 0 is the e the points and e vectors are understood to be identity. Without superscripts, 8 The conversion to C from Matlab was done by Stanley Lee; the shift-and-rotate spatial transforms and the superimposition were re-implemented in OpenGL by Orcun Goksel using texture mapping and alpha blending. 76 5.1. Notations and overview Sample preparation, registration Needletissue interaction model Matlab 5D database Repacking 1D database Offline Needle points in global frame Coordinate transform Needle points in probe frame Linecrossing local configuration θ-correction Interpolation Display Superimpose Needle image Shift and rotate Interpolated image background depth calculation OpenGL Figure 5.1: Flow chart of the needle insertion simulation measured in the base frame: we write points as P 1 , P 2 , etc., vector that e eframes as {O , C }, points from P 1 to P 2 as P1 P2 , right-handed coordinate 1 1 e e e in C etc., where C 1 is an orthonormal matrix. We write points measured 1 as P 11 , P 12 , etc., vectors in C 1 as P11 P21 , etc. With this convention, we can e vectors e write measured in the base frame using its representation in C 1 : 1 P1 P2 = C 1 (P 2 − P 11 ) = C 1 · (P11 P21 ).9 e of rotation e The sense is: looking into the frame origin from the rotation axis, counter-clockwise is the positive direction. In the simulation, there are four coordinate systems in use (refer to Fig. 5.2 and Fig. 2.2): 1. The tissue (global) coordinate system {Ot , C 0 }; e 2. The imaging plane coordinate system {Oi , C p } with one axis defining the norm and the other two the image eedge directions; Oi is located e at the bottom centre of the image; 9 A more coherent and general notation for coordinate systems representations and transforms can be found in [66] and [57]; the convention used here could be regarded as a shorthand that is sufficient for our purposes. 77 5.1. Notations and overview Pi e x Pn e z {Ot , C 0 } e y Op e Y P i+1 e X {Oi , C p } e Z Figure 5.2: The global frame, image frame, and probe frame. Compare to Fig. 2.2. 3. The probe coordinate system {Op , C p } located at the probe centre and e of the image; with the same orientation as that 4. The needle/robot coordinate system {P n , C n }. e The image simulation algorithm developed in Chapter 2 (the “interpolation” block) requires as input {P n , C n } measured with respect to {Op , C p }, e that is, the relative position andeorientation of the needle at the intersection with the imaging plane, measured in the probe frame. On the other hand, the needle-tissue interaction model [29] describes the needle as a series of points P i , i = 1, . . . , N in the global frame {Ot , C 0 }. We therefore need two e from the geometry given by steps toecalculate the local needle configuration the needle-tissue interaction model: 1. Transform the points description P i , i = 1, . . . , N from the global e frame to the probe frame (section 5.3); 2. Use a line crossing algorithm to obtain the position and orientation of the needle at the intersection with respect to {Op , C p }(section 5.4). e The 1D database required in the C implementation of the interpolation comes from re-packing of the prepared 5D Matlab database (section 5.2). Also, the configuration parameters need to be adjusted according to θ (section 5.5). After the interpolation, the resulting interpolated image needs to be shifted and rotated, and then superimposed on a tissue background as described in Chapter 2. 78 5.2. Re-packing In the following sections, we describe these various algorithms developed for the tasks in Fig. 5.1. 5.2 Re-packing In the Matlab implementation, the image database is stored as a 5D cell array with each element populated by either a sample image or by a NULL. The NULL indicates the needle image with a certain configuration was not collected during the sample collection. Such case may happen when the needle needs to move into the probe to achieve the desired position. In order to efficiently access the image database in C, we need to convert the 5D array to 1D. Because of the presence of empty elements due to unattainable needle configurations, we need to construct an extra array with flags that mark these unavailable samples. The C implementation thus operates on two 1D arrays: one with continuous storage of the sample images, the other with flags indicating the empty images. 5.3 5.3.1 Transformation from the global coordinate to the probe coordinate Problem description The problem is to convert a point P n from the global coordinate (here e as the base frame) to the probe coordinate. What we know is the image frame {Oi , C p } measured in the global frame {Ot , C 0 }, and the vector Oi Op e in the image frame; this vector liese in the XY plane and can measured be written as (λ, γ, 0)T . In the implementation of needle-tissue interaction model, XY Z is also called nsa. 5.3.2 Solution Referring to Fig. 5.2, in C 0 the base frame: Op Pn = P n − Op e e = P n − (Oi + Op Oi ) e e = P n − (Oi + C p (λ~e1 + γ~e2 )) e e Using the measurement in C p : Op Pn = C p P pn e 79 5.4. Line crossing Equating the above two, P pn = C Tp [P n − (Oi + C p (λ~e1 + γ~e2 ))] e e e = C Tp (P n − Oi ) − (λ~e1 + γ~e2 ) e e 5.4 5.4.1 Line crossing Problem description Given a series of points P i , i = 1, . . . , N , and the imaging plane {Oi , C p }, calculate the intersectione P n between the line segment connectinge the N e that the problem itself is frame independent, points and the plane. Note but we need the result to be in the probe frame, hence the pre-transform in section 5.3. Here we omit the superscript p and assume everything is measured with respect to the probe. 5.4.2 Solution P n is constrained by two relations: e ~ is perpendicular to the vector from the probe origin to the intersec1. Z tion: ~ =0 (P n − Op ) · Z (5.1) e e 2. P i , P n , P i+1 lie on a line: e e e P n − P i = c (P i+1 − P n ) (5.2) e e e e From 5.2: P i + cP i+1 Pn = e e 1+c e Substituting into Eq. 5.1, we can solve for c: (5.3) ~ (Op − P i ) · Z (5.4) e e ~ (P i+1 − Op ) · Z e e the sign of which determines if Pi Pi+1 crosses the imaging plane. Note that ~ =0 c = −1 ⇐⇒ P i = P i+1 , which is unrealistic. Also, (P i+1 − Op ) · Z e e e e means P i+1 is in the imaging plane. Theealgorithm to determine the crossing point is therefore: c= 80 5.5. θ-correction 1. Start from P 0 from the point list; e ~ = 0, P i+1 is the crossing point; return; 2. If (P i+1 − Op ) · Z e e e 3. Calculate c using Eq. 5.4; 4. If c ≥ 0, Pi Pi+1 crosses the imaging plane. The crossing point is given in Eq. 5.3. If c < 0, move on to the next P i+1 . e The above algorithm only gives the position of P n ; the orientation can be obtained by referring to Fig. 5.3. In the figure, e ~ n = Pi Pi+1 Z |Pi Pi+1 | ~ n · Z) ~ Z ~ =Z ~ nz (Z ~ n · X) ~ X ~ =Z ~ nx (Z ~ nz + Z ~ nx = Z ~ nxz Z we have, ~n · Y ~ ) atan ||Ẑn − Ẑnxz || α = 2 sign(Z ||Ẑn + Ẑnxz || ~ ~ n · X) ~ atan |Zn x| β = sign(Z ~ nz | |Z where we use the notation v̂ = ~v /||~v ||.10 The needle is assumed to be rigid for the twist around central axis: ρ angle propagates from the base to the tip. Therefore we obtain ρ directly from the user interface. 5.5 5.5.1 θ-correction Problem description As demonstrated in section 2.1.2, a needle image at an angle θ can be obtained by rotating the image at θ = 0: In (r, θ, l, ρ, α, β) = Rθ (In (r, 0, l, ρ0 , α0 , β0 )). The geometry is shown in Fig. 5.4: we have the roll-pitch-yaw angles ρ0 , α0 , β0 (corresponding to C n ) measured in C p ; we want to know ρ, α, β with respect to the rotated Rθ C p , so that we can retrieve the sample image (taken at θ = 0) with the correct configurations. 10 We have used the more accurate atan of half angle instead of the less accurate acos [44, 66]. 81 5.5. θ-correction Y Z~n ρ α ~ nz Z β ~ Znx ~ nxz Z X Z Figure 5.3: Calculation of the pitch (α) and yaw (β) angles from the probe and needle frame C {P n , C n } e Y Z θ X {Op , C p } e Figure 5.4: Frames involved in the θ correction 82 5.6. Depth calculation 5.5.2 Solution From the geometry, the problem is equivalent to finding ρ, α, β (corresponding to C = R−θ C n ) with respect to C p . The procedure is as follows: 1. Compute C n the rotation matrix from ρ0 , α0 , β0 ; this is straight-forward: C n = RY (β)RX (α)RZ (ρ) where RZ (ρ) means rotating around Z by ρ. 2. Rotate C n : C = RZ (−θ)C n 3. Calculate the new ρ, α, β from the element in C [48]: ρ = atan(r32 /r33 ) q 2 + r2 ) α = atan(−r31 / r32 33 β = atan(r21 /r11 ) 5.6 Depth calculation 5.6.1 Problem description One problem that was not previously addressed is the calculation of the needle/seed depth with respect to the imaging plane given the needle configuration. In Fig. 5.5, after rotating the needle around O by pitch and yaw e angles α and β, the new depth is l0 . The distance between the imaging plane and the robot motion plane is di . During the data collection, l is the 0 ), we wish parameter that we manipulate; in the simulation, given l0 (or l⊥ to know the l used in the experiment, so that we can retrieve the images needed for interpolation. In short, we want to write l as a function of l0 . 5.6.2 Solution To derive this relationship, we need to take into account how l was manipulated during the data collection. On one hand, during the collection of the needle images, l was changed by inserting and retracting the needle with di fixed. On the other hand, during the collection of the seed images, l was changed by moving the probe, that is, changing di with d fixed. These two techniques are different. 83 5.6. Depth calculation O e d ~s l h ′ l⊥ s~′ l′ di Imaging Plane Figure 5.5: The geometry of the needle rotating around the robot wrist rotation centre. From Fig. 5.5, we have: ~s = (0, h, d)T cos β 0 sin β 1 0 0 0 1 0 0 cos α − sin α h s~0 = 0 − sin β 0 cos β 0 sin α cos α d Extracting the third dimension: 0 + di = d cos α cos β + h sin α cos β l⊥ Also we have: 0 l⊥ = l0 cos α cos β l = d − di When di is a constant, that is, when the needle depth l is changed by moving the needle, we have: l= di + l0 cos α cos β − h sin α cos β − di cos α cos β When d is a constant, that is, when l is changed by moving the probe, we have: l = (1 − cos α cos β)d + l0 cos α cos β − h sin α cos β From these two relations, we can obtain l from l0 for both ways of changing l. 84 5.7. Needle insertion simulator Figure 5.6: The needle insertion simulator graphics interface (courtesy of Orcun Goksel). The needle is pointed out. 5.7 Needle insertion simulator The above algorithms and the needle simulation method developed in Chapter 2 are integrated into the needle insertion simulator developed by Orcun Goksel. In this simulator, the virtual needle could be inserted interactively using a haptic device; the needle image is alpha-blended into the background tissue image which deforms under the needle load. 85 Chapter 6 Conclusion 6.1 Summary Real-time simulation of ultrasound images of the needle is an essential building block for needle insertion procedure simulator, which could be beneficial for clinical training of interventional procedures. The physical modelling of scattering from metallic objects is not only useful for physically simulating ultrasound artifacts associated with the needle, but is also potentially helpful for seed detection in brachytherapy. In this thesis, a real-time 6DOF needle and seed B-mode image simulation method has been developed using multi-dimensional interpolation. With the intention to understand needle image artifacts forming mechanisms, the scattering from metal plates and spheres are studied and simulated using physical models. The physicallybased simulations are demonstrated to be able to generate wave forms that resemble the bright band/tail and dark shadow artifacts and are individually validated by comparing to experimental results. The physically-based simulation technique developed in Chapter 3 and 4 has the following drawbacks compared with the real-time simulation technique developed in Chapter 2: 1. The physically-based simulation is slow – generating a single RF line takes several minutes; 2. The physically-based simulation discussed could only deal with “canonical” geometries – the plates and the spheres in our case. More complex techniques exist for modelling special geometries but require considerably more effort. Comparatively, the real-time simulation developed in the previous chapter does not depend on the specific geometry of the scattering object. Indeed, it is demonstrated to work for both the needle and seed using essentially the same method. The physically based simulation has its own merits: 1. The phenomenon could be modelled mathematically, thus shedding light on the internal mechanisms; 86 6.2. Contributions 2. The theoretical and numerical studies may complement and, if carried to enough accuracy, substitute experiments to some degree. Given a “perfect” physically based simulator, we could computationally generate the database needed for the interpolation based simulator without the labour of collecting the images using the robot. However, the current physical simulation is far from this capability. 6.2 Contributions The novelties and contributions of this thesis are: Fast pose-dependent image simulation: The novel image synthesis scheme based on tensor-product interpolation provides an easy and fast way of generating realistic new poses from a finite number of sample images. This technique is especially useful for interactive simulations. Robot-aided ultrasound image collection: The use of a robot to collect pose-dependent sample images proves to be accurate, efficient, and repeatable. Reciprocally, the registration process in the simulation method helps measure and test the accuracy of the robot. Experimental studies of needle visibility and artifacts: The needle visibility and bright tail and dark shadow artifacts are studied experimentally in the transverse plane. It is shown that the bright tail only occurs near orthogonal incidence for the shaft. The needle tip behaves differently from the shaft and may still be visible even if the pitch angle is large. These characteristics will be useful for understanding and improving needle insertion procedures in the transverse plane. Physical models of scattering from plates and spheres in ultrasound: The scattering of ultrasound from metal plates are modelled by ray tracing and validated by water tank experiments. The bright band pattern of the needle in the sagittal plane could be interpreted using this model. It could be inferred from the results that the RF signal of the ultrasound machine is proportional to the acoustic intensity. The scattering from metal spheres are modelled by solving the wave equations numerically. The results indicate that the bright tail artifact of the sphere in B-mode could be understood by the RF line simulated by linear acoustic models. The similar bright tail artifact generated by 87 6.3. Future work the cylindrically shaped needle could be qualitatively analyzed using the spherical simulation. Numerical study of the conditioning problem for spherical scatterer: A potential numerical pitfall that may undermine the accuracy of the spherical scattering simulation is discussed. It is pointed out that the numerical difficulty may not influence the credibility of the scattering solutions, but may be problematic for the transmitting solutions. The problem seems to be present in the previous works on similar topics but was not seen discussed. 6.3 Future work More validation of the transverse image simulation: as noted in Section 2.1.5 and 2.4, the previous validation experiment may have small misalignment and phantom fracture issue. A more careful validation should eliminate these factors and the correlation scores should be higher than those reported in Section 2.1.5. Comparison with dark shadow simulation using computer graphics: in section 2.1.4, a dark shadow simulation method is proposed based on experimental observations. Various other ways to simulate the dark shadow using computer graphics techniques were previously studied, such as the use of ray casting in [50]. A comparison is needed to judge which generates more realistic shadows. Needle bending: if the needle is not rigid, the pitch and yaw angles might be continuously changing when the needle is pushed sideways. From the experiments we observe that the changing pitch angle may cause continuous changes in the bright tail pattern. The current nearest neighbor interpolation method for the pitch angle may not be able to reflect this continuous change. An image database sampled with finer pitch angle intervals and the use of linear interpolation for the pitch angle may be needed to better simulate the bending needle. Integration of the sagittal and seed image simulation into the simulator: the needle insertion simulator discussed in Section 5.7 at present only supports simulation of the needle image in the transverse plane. The extension to the sagittal plane and adding seed depositing simulation requires a moderate amount of effort to convert current Matlab implementation to C++ and OpenGL. 88 6.3. Future work Scattering from infinite and finite cylinders: scattering from sphere is studied in this thesis to explain the physical origin of the bright tail. To be able to quantitatively analyze the scattering from the needle and seed, a cylindrical scattering model needs to be constructed. The study of scattering from an infinite cylinder should parallel that of the spherical scatterer nearly exactly as was shown in various previous publications; the modelling of the finite cylinder and oblique incidence needs more efforts (see Table 1.1). Physically-based B-mode simulation of large metallic scatterers: By computing the scattering stress in all directions, it should be possible to construct the whole B-mode image by assembling multiple RF lines. A complication is that when the received traction is at oblique angle with the transducer, the wave is foreshortened, resulting in weaker signal the further away the transducer is from the backscattering direction. Seed detection: The results in this thesis suggest two ways for seed detection: 1. The collected B-mode seed image database may be used directly as templates for detection; simple volume correlation or more sophisticated machine learning algorithms may be applied. 2. The collected seed RF database may be analyzed and the characteristic frequency components recognized for detection purpose. The spherical scattering model studied in Chapter 4 might be useful to recognize some features of the spectrum, such as the repetitions of nulls and peaks in Fig. 4.3. Boundary value problem conditioning: The conditioning issue for the boundary value problem of spherical scatterer might have been noticed in previous studies. A more thorough literature review on this subject will be helpful. Physically-based simulation of scattering from non-canonical geometries: for non-standard geometries such as the needle tip, the simulation of scattering might or might not be possible by using Finite Element Methods [33, 39, 74]. Boundary element method may also be considered. An introduction is [15]. 89 Bibliography [1] D. Aiger and D. Cohen-Or. Real-Time Ultrasound Imaging Simulation. Real-Time Imaging, 4(4):263–274, 1998. [2] JS Allen, DE Kruse, and KW Ferrara. Shell waves and acoustic scattering from ultrasound contrast agents. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 48(2):409–418, 2001. [3] M.B. Allen and E.L. Isaacson. Numerical Analysis for Applied Science. Wiley-Interscience, 1998. [4] M.E. Anderson and G.E. Trahey. A seminar on k-space applied to medical ultrasound. Department of Biomedical Engineering, Duke University, 2002. [5] L. Avruch and PL Cooperberg. The ring-down artifact. Journal of ultrasound in medicine: official journal of the American Institute of Ultrasound in Medicine, 4(1):21, 1985. [6] VM Ayres, GC Gaunaurd, CY Tsui, and MF Werby. The effects of lamb waves on the sonar cross-sections of elastic spherical shells. International journal of solids and structures, 23(7):937–946, 1987. [7] Ali Baghani. A simulation software for the analysis of ultrasound interaction with multi-layer spherical scatterers. Technical report, University of British Columbia, 2008. [8] J.E. BAILEY, R.O. BUDE, and T. TUTHILL. US artifacts: Effects on out-of-plane US images reconstructed from three-dimensional data sets. Radiology, 218(2):592–597, 2001. [9] GR Barnard and CM McKinney. Scattering of Acoustic Energy by Solid and Air-Filled Cylinders in Water. The Journal of the Acoustical Society of America, 33:226, 1961. [10] B. Barrowes et al. Matlab routines for computation of special functions. Online document, http://ceta. mit. edu/ceta/comp spec func, 2003. 90 Chapter 6. Bibliography [11] S. BONDESTAM and J. KREULA. Needle Tip Echogenicity: A Study with Real Time Ultrasound. Investigative Radiology, 24(7):555, 1989. [12] L.M. Brekhovskikh. Waves in layered media. Academic Pr, 1980. [13] WM Butler, S.O. Center, and WV Wheeling. Review of modern prostate brachytherapy. Engineering in Medicine and Biology Society, 2000. Proceedings of the 22nd Annual International Conference of the IEEE, 1, 2000. [14] É. Cardinal, R.K. Chhem, and C.G. Beauregard. Ultrasound-guided interventional procedures in the musculoskeletal system. Radiologic Clinics of North America, 36(3):597–604, 1998. [15] S. Chandler-Wilde and S. Langdon. Boundary element methods for acoustics. Department of Mathematics, University of Reading, 2007. [16] GA Chapman, D. Johnson, and AR Bodenham. Visualisation of needle position using ultrasonography. Anaesthesia, 61(2):148, 2006. [17] D.A. Christensen. Ultrasonic bioinstrumentation. John Wiley & Sons Inc, 1988. [18] Ultrasonix Medical Corporation. Frequency spectrum analysis. Ultrasonix Research Users Forum. http://research.ultrasonix.com/ viewtopic.php?f=2&t=363&start=0&hilit=spectrum, April 2009. [19] M.J. Crocker. Handbook of acoustics. Wiley-Interscience, 1998. [20] B.J. Davis, R.R. Kinnick, M. Fatemi, E.P. Lief, R.A. Robb, and J.F. Greenleaf. Measurement of the ultrasound backscatter signal from three seed types as a function of incidence angle: Application to permanent prostate brachytherapy. International journal of radiation oncology, biology, physics, 57(4):1174–1182, 2003. [21] GD Dodd 3rd, CC Esola, DS Memel, AA Ghiatas, KN Chintapalli, EK Paulson, RC Nelson, JV Ferris, and RL Baron. Sonography: the undiscovered jewel of interventional radiology. Radiographics, 16(6):1271–1288, 1996. [22] D. dAulignac, C. Laugier, J. Troccaz, and S. Vieira. Towards a realistic echographic simulator. Medical Image Analysis, 10(1):71–81, 2006. [23] J.J. Faran Jr. Sound scattering by solid cylinders and spheres. The Journal of the Acoustical Society of America, 23:405, 1951. 91 Chapter 6. Bibliography [24] G. Gaunaurd and W. Wertman. Transient acoustic scattering by fluidloaded elastic shells. Int. J. Solids Struct., 27:699–711, 1991. [25] GC Gaunaurd, HC Strifors, U.S.N.S.W. Center, and S. Spring. Frequency-and time-domain analysis of the transient resonancescattering resulting from the interaction of a sound pulse withsubmerged elastic shells. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 40(4):313–324, 1993. [26] O. Gerovich, P. Marayong, and A.M. Okamura. The effect of visual and haptic feedback on computer-assisted needle insertion. Computer Aided Surgery, 9(6):243–249, 2004. [27] O. Goksel and S. Salcudean. B-Mode Ultrasound Image Simulation in Deformable 3D Medium. IEEE transactions on medical imaging, 2009. [28] O. Goksel, SE Salcudean, and SP Dimaio. 3D simulation of needletissue interaction with application to prostate brachytherapy. Comput Aided Surg, 11(6):279–88, 2006. [29] O. Goksel, S.E. Salcudean, S.P. DiMaio, R. Rohling, and J. Morris. 3D Needle-Tissue Interaction Simulation for Prostate Brachytherapy. Proceedings of MICCAI, 3749:827–834, 2005. [30] Orcun Goksel. Ultrasound image and 3d finite element based tissue deformation simulator for prostate brachytherapy. Master’s thesis, University of British Columbia, 2004. [31] TJ Hall, M. Bilgen, MF Insana, and TA Krouskop. Phantom materials for elastography. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 44(6):1355–1365, 1997. [32] LD Hampton and CM McKinney. Experimental study of the scattering of acoustic energy from solid metal spheres in water. The Journal of the Acoustical Society of America, 33:664, 1961. [33] I. Harari. A survey of finite element methods for time-harmonic acoustics. Computer Methods in Applied Mechanics and Engineering, 195(1316):1594–1607, 2006. [34] M.D. Harpen. Partial wave analysis of the ultrasound comet tail artifact. Medical Physics, 27:2072, 2000. 92 Chapter 6. Bibliography [35] T. Hasegawa, Y. Hino, A. Annou, H. Noda, M. Kato, and N. Inoue. Acoustic radiation pressure acting on spherical and cylindrical shells. The Journal of the Acoustical Society of America, 93:154, 1993. [36] R. Hickling. Analysis of echoes from a solid elastic sphere in water. J. Acoust. Soc. Am, 34(10):1582–1592, 1962. [37] Y. Hu, S. Qin, and Q. Jiang. Characteristics of acoustic scattering from a double-layered micro shell for encapsulated drug delivery. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 51(7):808–820, 2004. [38] J. Huang, J.K. Triedman, N.V. Vasilyev, Y. Suematsu, R.O. Cleveland, and P.E. Dupont. Imaging Artifacts of Medical Instruments in Ultrasound-Guided Interventions. J Ultrasound Med, 26:1303–1322, 2007. [39] F. Ihlenburg. Finite element analysis of acoustic scattering. Springer Verlag, 1998. [40] J.A. Jensen. FIELD: A Program for Simulating Ultrasound Systems. MEDICAL AND BIOLOGICAL ENGINEERING AND COMPUTING, 34:351–352, 1996. [41] JA Jensen. Simulation of advanced ultrasound systems using Field II. In Biomedical Imaging: Macro to Nano, 2004. IEEE International Symposium on, pages 636–639, 2004. [42] JA Jensen and I. Nikolov. Fast simulation of ultrasound images. In Ultrasonics Symposium, 2000 IEEE, volume 2, 2000. [43] M.C. Junger and D. Feit. Sound, structures, and their interaction. MIT press Cambridge, MA, 1972. [44] W. Kahan. Lectures on computational aspects of geometry. Unpublished manuscripts, pages 3–13, 1983. [45] C.F. Keogh and P.L. Cooperberg. Is It Real or Is It an Artifact. Ultrasound Quarterly, 17(4):201, 2001. [46] RM Kirberger. Imaging artifacts in diagnostic ultrasound: a review. Veterinary radiology and ultrasound: the official journal of the American College of Veterinary Radiology and the International Veterinary Radiology Association (USA), 1996. 93 Chapter 6. Bibliography [47] SEM Langley and R. Laing. Prostate brachytherapy has come of age: a review of the technique and results. BJU International, 89(3):241–249, 2002. [48] S.M. LaValle. Planning algorithms. Cambridge University Press, 2006. [49] M. Levoy and P. Hanrahan. Light field rendering. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pages 31–42. ACM New York, NY, USA, 1996. [50] D. Magee, Y. Zhu, R. Ratnalingam, P. Gardner, and D. Kessel. An augmented reality simulator for ultrasound guided needle placement training. Medical and Biological Engineering and Computing, 45(10):957– 967, 2007. [51] J. Mamou and E.J. Feleppa. Singular spectrum analysis applied to ultrasonic detection and imaging of brachytherapy seeds. The Journal of the Acoustical Society of America, 121:1790, 2007. [52] J. Mamou, S. Ramachandran, and E.J. Feleppa. Angle-dependent ultrasonic detection and imaging of brachytherapy seeds using singular spectrum analysis. The Journal of the Acoustical Society of America, 123:2148, 2008. [53] TA Matalon and B. Silver. US guidance of interventional procedures. Radiology, 174(1):43–47, 1990. [54] H. Maul, A. Scharf, P. Baier, M. Wuestemann, HH Guenter, G. Gebauer, and C. Sohn. Ultrasound simulators: experience with the SonoTrainer and comparative review of other training systems. Ultrasound in Obstetrics and Gynecology, 24(5):581–585, 2004. [55] FG Mitri, BJ Davis, JF Greenleaf, and M. Fatemi. In vitro comparative study of vibro-acoustography versus pulse-echo ultrasound in imaging permanent prostate brachytherapy seeds. Ultrasonics, 49(1):31–38, 2009. [56] P.M.C. Morse and K.U. Ingard. Theoretical Acoustics. Princeton University Press, 1986. [57] R.M. Murray, Z. Li, and S.S. Sastry. A mathematical introduction to robotic manipulation. CRC, 1994. 94 Chapter 6. Bibliography [58] K. Nichols, L.B. Wright, T. Spencer, and W.C. Culp. Changes in Ultrasonographic Echogenicity and Visibility of Needles with Changes in Angles of Insonation. Journal of Vascular and Interventional Radiology, 14(12):1553–1557, 2003. [59] O. Nobuyuki. A threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern, 9:62–6, 1979. [60] Onda Corporation. Acoustic Properties of Solids, 2003. [61] A.V. Oppenheim and R.W. Schafer. Discrete-time signal processing. Prentice-Hall, Inc. Upper Saddle River, NJ, USA, 1989. [62] T. Poggio and R. Brunelli. A Novel Approach to Graphics. AI Memo, 1992. [63] J. Prosise. Programming Windows with MFC (with CD-ROM). Microsoft Press Redmond, WA, USA, 1999. [64] R. Rohling, A. Gee, and L. Berman. A comparison of freehand threedimensional ultrasound reconstruction techniques. Med Image Anal, 3(4):339–59, 1999. [65] JM Rubin, RS Adler, RO Bude, JB Fowlkes, and PL Carson. Clean and dirty shadowing at US: a reappraisal. Radiology, 181(1):231–236, 1991. [66] S.E. Salcudean. Introduction to robotics course notes. [67] S.E. Salcudean, T.D. Prananta, W.J. Morris, and I. Spadinger. A robotic needle guide for prostate brachytherapy. In Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on, pages 2975–2981, 2008. [68] I. Schafhalter-Zoppoth, C.E. McCulloch, and A.T. Gray. Ultrasound visibility of needles used for regional nerve block: An in vitro study. Regional Anesthesia and Pain Medicine, 29(5):480–488, 2004. [69] K. Schwenk, G. Reis, and R.H. van Lengen. Realtime artificial ultrasound images. Submitted to 2005/06 Elsevier Journal Simulation Modelling Practice and Theory, 2005. [70] ML Skolnick. Estimation of ultrasound beam width in the elevation (section thickness) plane. Radiology, 180(1):286–8, 1991. 95 Chapter 6. Bibliography [71] TK Stanton. Sound scattering by cylinders of finite length. II. Elastic cylinders. The Journal of the Acoustical Society of America, 83:64, 1988. [72] G. Strang. Introduction to applied mathematics. Wellesley Cambridge Pr, 1986. [73] P. Thevenaz, T. Blu, and M. Unser. Image interpolation and resampling. Handbook of Medical Imaging, Processing and Analysis, pages 393–420, 2000. [74] L.L. Thompson. A review of finite-element methods for time-harmonic acoustics. The Journal of the Acoustical Society of America, 119:1315, 2006. [75] F.P. Vidal, N. Chalmers, D.A. Gould, A.E. Healey, and N.W. John. Developing a needle guidance virtual environment with patient-specific data and force feedback. International Congress Series, 1281:418–423, 2005. [76] FP Vidal, NW John, AE Healey, and DA Gould. Simulation of ultrasound guided needle puncture using patient specific data with 3D textures and volume haptics. Computer Animation and Virtual Worlds, 19(2), 2008. [77] E. W. Weisstein. Bessel function of the second kind. MathWorld – A Wolfram Web Resource. http://mathworld.wolfram.com/ BesselFunctionoftheSecondKind.html, 2009. [78] Z. Ye. On sound scattering and attenuation of Albunex bubbles. The Journal of the Acoustical Society of America, 100:2011, 1996. [79] C. Zhang and T. Chen. A survey on image-based renderingrepresentation, sampling and compression. Signal Processing: Image Communication, 19(1):1–28, 2004. [80] Y. Zhu, D. Magee, R. Ratnalingam, and D. Kessel. A virtual ultrasound imaging system for the simulation of ultrasoundguided needle insertion procedures. Proceedings of Medical Image Understanding and Analysis (MIUA), 1:61–65, 2006. [81] Y. Zhu, D. Magee, R. Ratnalingam, and D. Kessel. A training system for ultrasound-guided needle insertion procedures. In Medical image 96 computing and computer-assisted intervention: MICCAI, volume 10, page 566, 2007. [82] MC Ziskin, DI Thickman, NJ Goldenberg, MS Lapayowker, and JM Becker. The comet tail artifact. Journal of ultrasound in medicine, 1(1):1, 1982. 97 Appendix A A derivation of reflection and transmission coefficients at fluid-solid interface Here we derive from the first principles the reflection and refraction coefficients with orthogonal incidence from the fluid to the solid. We roughly follow the notation in [39]. A.1 Solid ~ using the Helmholtz In the solid: after decomposing the displacement field U decomposition: ~ = ∇Φ + ∇ × Ψ ~ U (A.1) and assuming zero body force, isotropic material, with Cauchy strain and linear stress-strain relationship, the equilibrium equation in the elastic solid ∇ · σ = ρs ~ ∂2U ∂t2 is equivalent to the wave equations: ∂2Φ = α2 ∇ · ∇Φ ∂t2 ~ ∂2Ψ ~ = β 2 ∇ · ∇Ψ ∂t2 q λ+2µ where σ is the stress tensor, α = is the longitudinal wave speed, ρs q µ β = ρs is the shear wave speed, ρs is the solid density, λ, µ are the Lamé 98 A.2. Fluid constants. Further assuming time harmonic fields11 , these are equivalent to the Helmholtz equations: ∇ · ∇Φ + kα2 Φ = 0 ~ + k2 Ψ ~ =0 ∇ · ∇Ψ (A.2) β where kα = A.2 ω α , kβ = ω β are the longitudinal and shear wave numbers. Fluid In the fluid: assuming small oscillation and non-viscous ideal fluid, the first order time-harmonic wave equation can be derived from the continuity condition and the equation of motion; the result is: ∇ · ∇p + k 2 p = 0 (A.3) where k = ωc is the wave number and c is the speed of sound, which can be determined by the compressibility and density of the fluid [56]. A.3 Boundary conditions On the boundary, two boundary conditions can be used: 1. The compatibility condition: the normal displacement on the surface must be continuous. The fluid Euler equation under time harmonic assumption ρf ω 2 ~u · ~n = ∇p · ~n where ~u is the fluid displacement and ~n is the normal vector, can be expressed in terms of the normal displacement of the solid: ~ · ~n = ∇p · ~n ρf ω 2 U (A.4) 2. The equilibrium condition: the outward normal surface traction of the solid balances the inward pressure of the fluid: (σ · ~n) · ~n = −p (A.5) 11 Time harmonic field variables are in the form of C(ω)e−iωt , where C is the spatial part of the fluctuation. To arrive at this form, we essentially write a field variable as a combination of a DC (steady-state) part and a time-dependent fluctuation part: c + C(ω)e−iωt . The former does not enter the differential equation. Thus we can consider only the latter. 99 A.4. Solution z T−n 0 pinc Solid Fluid ps Figure A.1: Geometry of the problem If we multiply both sides by ~n, we can see this means an increase of traction in the outward normal direction corresponds to an increase of pressure in the inward normal direction. A.4 Solution The geometry of the problem indicates that Cartesian coordinates are best suited for a solution. We choose the incident direction as positive z, as shown in Fig. A.1. The total pressure in the fluid is p = pinc + ps The solutions to the Helmholtz wave equations such as Eq. A.2 and A.3 in the Cartesian coordinates can be obtained by separation of variables; the result is in the form p = X(x)Y (y)Z(z) = ei(kx x+ky y+kz z) The incident pressure which satisfies Eq. A.3 has only the z-component; we normalize everything so that it has the unit amplitude: pinc = eikz (A.6) We further assume the reflected and refracted solutions are also constant along the x, y directions (one should be able to argue that this is the only possibility, similar to the proof of the Snell’s law, but allow me to be a little sloppy here). We can write the transmitted (i.e. refracted) and scattered (i.e. reflected) solutions to the wave equations A.2 and A.3 as: Φ = CΦ eikα z 100 A.4. Solution CΨ1 ~ = CΨ eikβ z Ψ 2 CΨ3 ps = V e−ikz where V is called the reflection coefficient. Applying Eq. A.1 again, we get the displacement: −CΨ2 ikβ eikβ z ~ (z) = CΨ ikβ eikβ z U (A.7) 1 CΦ ikα eikα z Also since the outward norm is ~n = (0, 0, −1)T , the normal traction is (σ · ~n) · ~n = σ33 From the material law: σ33 = 2µ ∂Uz ∂Ux ∂Uy ∂Uz + λ( + + ) = −(2µ + λ)CΦ kα2 eikα z ∂z ∂x ∂y ∂z (A.8) This means that the transmitted stress does not contain the shear component – “mode conversion” does not occur at orthogonal incidence. Substituting Eq. A.7 A.8 and the solutions to the wave equations to Eq. A.4 A.5: ρf ω 2 CΦ kα eikα z = −kV e−ikz + keikz at z = 0 (2µ + λ)CΦ kα2 eikα z = eikz + V e−ikz From the solutions we have: p = V e−ikz = 2µ + λ − αcρf e−ikz = ρs α − ρf c e−ikz s 2µ + λ + αcρf ρs α + ρf c T−n = −σ33 = (1 + V )eikα z = W eikα z (A.9) where T−n is the traction in the positive z direction, i.e. the transmitted pressure; W is the transmission coefficient. Note that the result is frequency independent. The results agree with the form when both media are fluid. This is understandable since mode conversion does not occur. A thorough discussion of wave propagation in the layered media can be seen in the monograph [12]. 101 Appendix B Implementation of automatic image collection using the needle guidance robot The goal is to automatically store and download the B-mode and RF frames from the Ultrasonix machine after each desired needle position is reached. This can be achieved by a finite state machine as shown in Fig. B.1. The “NextPlan” sends motor commands to and receives encoder/caliper readings from all joints as was implemented in [67]. To acquire the images after the desired position is attained, we have implemented “AcquireImg” function using the Ulterius SDK developed by Ultrasonix; this SDK supports remotely connecting to and controlling the Ultrasonix US machine through the Internet. This “AcquireImg” function is implemented as a separate thread which is initiated by “NextPlan”. To synchronize the two threads, MFC CEvent class has been used [63]. Before starting, four events representing the reaching of desired positions by the four motors are reset. After each encoder/caliper reports the desired position is attained, the corresponding event is set. In the “AcquireImg” thread, the image collection function initially blocks on the four events because of the reset status. After these events are all set by the reaching of desired positions, a few frames are frozen on the US machine and transferred OnAcquisitionComplete NextPlan AcquireImg all motors reaching desired positions Figure B.1: The image acquisition automaton. 102 Appendix B. Implementation of automatic image collection using the needle guidance robot to a local folder. After the images of current position are collected, “AcquireImg” posts “AcquisitionComplete” message which results in a call to “NextPlan”, closing the loop. The needle is moved to the next pre-plan position and the whole process repeats itself until the last pre-plan. Note that the events are reset automatically after being set so that the “AcquireImg” will block each time “NextPlan” is in action. The program also supports manual saving, in which after acquiring the image, “NextPlan” is not called, so the user could adjust the needle and save the image manually. This is useful during phantom experiments where the needle cannot be moved inside the phantom and has to be inserted by hand. 103
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Real-time B-mode ultrasound image simulation and artifacts...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Real-time B-mode ultrasound image simulation and artifacts modelling of needles and brachytherapy seeds Zhu, Mengchen 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Real-time B-mode ultrasound image simulation and artifacts modelling of needles and brachytherapy seeds |
Creator |
Zhu, Mengchen |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | The training of ultrasound guided interventional procedures could benefit from ultrasound simulators that synthesize realistic-looking B-mode images of tissue and interventional tools, such as the needle. For a prostate brachytherapy simulator in particular, both the needle and seed images need to be synthesized. In this thesis, we propose an interpolation-based method for simulating needle and seed B-mode images in real time. We parametrize a needle or a seed B-mode image as a function of its position and orientation. We collect needle and seed images under various spatial configurations in a water-tank using a needle guidance robot. Then we use multi-dimensional tensor-product interpolation to simulate images of needles and seeds with arbitrary poses and positions using collected images. After further processing, the interpolated needle and seed images are superimposed on top of phantom or tissue image backgrounds. The similarity between the simulated and the real images is measured using a correlation metric. The results are also compared to in vivo images. Images in both the transverse and the sagittal planes are simulated. Artifacts associated with the needle and seed are studied in detail by experiments and by physical models. It is demonstrated that the bright band pattern associated with the sagittal needle images could be interpreted using a plate reflection model; the bright tails associated with the transverse needle images could be analyzed by simulating the back-scattered stress field from a sphere. |
Extent | 3700235 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-08-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067656 |
URI | http://hdl.handle.net/2429/12611 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_fall_zhu_mengchen.pdf [ 3.53MB ]
- Metadata
- JSON: 24-1.0067656.json
- JSON-LD: 24-1.0067656-ld.json
- RDF/XML (Pretty): 24-1.0067656-rdf.xml
- RDF/JSON: 24-1.0067656-rdf.json
- Turtle: 24-1.0067656-turtle.txt
- N-Triples: 24-1.0067656-rdf-ntriples.txt
- Original Record: 24-1.0067656-source.json
- Full Text
- 24-1.0067656-fulltext.txt
- Citation
- 24-1.0067656.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067656/manifest