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 bene- fit from ultrasound simulators that synthesize realistic-looking B-mode im- ages 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 simulat- ing needle and seed B-mode images in real time. We parametrize a nee- dle 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 ar- bitrary poses and positions using collected images. After further processing, the interpolated needle and seed images are superimposed on top of phan- tom 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 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivations and overview . . . . . . . . . . . . . . . . . . . . 1 1.2 Background and literature review . . . . . . . . . . . . . . . 2 1.2.1 Prostate brachytherapy . . . . . . . . . . . . . . . . . 2 1.2.2 B-mode Ultrasound image simulation . . . . . . . . . 2 1.2.3 Needle image characterization and visibility . . . . . 4 1.2.4 Needle image simulation . . . . . . . . . . . . . . . . 5 1.2.5 Image-based rendering . . . . . . . . . . . . . . . . . 6 1.2.6 Ultrasound artifacts . . . . . . . . . . . . . . . . . . . 7 1.2.7 Scattering from cylindrical and spherical objects . . . 8 1.2.8 Seed detection . . . . . . . . . . . . . . . . . . . . . . 10 2 Needle and brachytherapy seeds in B-mode ultrasound: real- time simulation and experimental studies . . . . . . . . . . 11 2.1 Transverse needle image simulation and tip visibility study . 11 2.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.2 Simulation methods . . . . . . . . . . . . . . . . . . . 12 2.1.3 Data collection method . . . . . . . . . . . . . . . . . 18 2.1.4 Superimposition and shadow simulation . . . . . . . . 20 2.1.5 Results and validation . . . . . . . . . . . . . . . . . . 23 iii Table of Contents 2.2 Sagittal needle image simulation . . . . . . . . . . . . . . . . 26 2.2.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.2 Results and Discussions . . . . . . . . . . . . . . . . . 29 2.3 Seed image simulation in the transverse and sagittal plane . 30 2.3.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Results and Discussions . . . . . . . . . . . . . . . . . 35 2.4 Experimental characterization of the transverse shadow form- ing conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4.1 The bright tail . . . . . . . . . . . . . . . . . . . . . . 37 2.4.2 The dark shadow . . . . . . . . . . . . . . . . . . . . 41 2.4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 42 3 Modelling of the bright band artifact in the sagittal plane 43 3.1 A description of the phenomenon . . . . . . . . . . . . . . . 44 3.2 Modelling of the phenomenon . . . . . . . . . . . . . . . . . 44 3.3 Validation of the model . . . . . . . . . . . . . . . . . . . . . 47 3.3.1 Sound speed . . . . . . . . . . . . . . . . . . . . . . . 48 3.3.2 Band intensity . . . . . . . . . . . . . . . . . . . . . . 50 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4.1 Incident direction . . . . . . . . . . . . . . . . . . . . 53 3.4.2 Quarter-wave layer . . . . . . . . . . . . . . . . . . . 53 3.4.3 Surface roughness . . . . . . . . . . . . . . . . . . . . 53 3.4.4 Relation to the needle image . . . . . . . . . . . . . . 54 4 Modelling of the bright tail and dark shadow artifacts in the transverse plane . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2 Hypotheses for the bright tail and dark shadow forming mech- anism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3 Boundary value problem for scattering from spherical shells . 57 4.4 Scattering spectrum simulation . . . . . . . . . . . . . . . . . 60 4.4.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 63 4.5 Transient simulation . . . . . . . . . . . . . . . . . . . . . . . 64 4.5.1 Modelling of the ultrasound pulse . . . . . . . . . . . 65 4.5.2 Validation of US pulse modelling . . . . . . . . . . . . 67 4.5.3 Pulse response . . . . . . . . . . . . . . . . . . . . . . 67 4.5.4 Simulation of transient back-scattering from a spher- ical shell . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.6 Validations and interpretations of the transient response . . 68 4.6.1 Wave speed . . . . . . . . . . . . . . . . . . . . . . . 68 iv Table of Contents 4.6.2 Resemblance to response from needle shaft in the trans- verse plane . . . . . . . . . . . . . . . . . . . . . . . . 69 4.6.3 Dark shadow from phase cancellation . . . . . . . . . 69 4.6.4 Solid vs. shell . . . . . . . . . . . . . . . . . . . . . . 70 4.7 Numerical issues . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.7.1 Convergence of the incident wave series . . . . . . . . 72 4.7.2 Condition number of the linear system . . . . . . . . 72 4.7.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 75 5 Implementation of the ultrasound simulator for interactive needle insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.1 Notations and overview . . . . . . . . . . . . . . . . . . . . . 76 5.2 Re-packing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.3 Transformation from the global coordinate to the probe co- ordinate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.3.1 Problem description . . . . . . . . . . . . . . . . . . . 79 5.3.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.4 Line crossing . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.4.1 Problem description . . . . . . . . . . . . . . . . . . . 80 5.4.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.5 θ-correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.5.1 Problem description . . . . . . . . . . . . . . . . . . . 81 5.5.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.6 Depth calculation . . . . . . . . . . . . . . . . . . . . . . . . 83 5.6.1 Problem description . . . . . . . . . . . . . . . . . . . 83 5.6.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.7 Needle insertion simulator . . . . . . . . . . . . . . . . . . . 85 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Appendices A A derivation of reflection and transmission coefficients at fluid-solid interface . . . . . . . . . . . . . . . . . . . . . . . . . 98 v Table of Contents A.1 Solid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 A.2 Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 A.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 99 A.4 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 B Implementation of automatic image collection using the nee- dle 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. . . . . . . . . . . . 9 2.1 Sample range and number of samples . . . . . . . . . . . . . . 19 2.2 Sample range and numbers. . . . . . . . . . . . . . . . . . . . 28 2.3 The range and numbers of the seed image samples in the sagittal and the transverse planes. . . . . . . . . . . . . . . . 35 4.1 Material properties used in the simulation . . . . . . . . . . . 64 vii List of Figures 1.1 A sketch of the geometry of ultrasound guided prostate brachyther- apy (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 1.2 The needle used for prostate brachytherapy. The arrow on the base indicates the bevel tip direction (reprinted from [30]). 3 1.3 The brachytherapy seeds (reprinted from [20]). . . . . . . . . 3 1.4 Lead pellets in the liver producing bright tail (reprinted from [45]). 7 2.1 Needle insertion process in the water and in the phantom . . 12 2.2 Coordinate assignment (not to physical scale). ρ is the roll angle which measures how much the needle rotates around its axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Validation of θ invariance. After rotating the first image by θ, the two images have a correlation score of 0.96. . . . . . . 14 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 sub- space, then it bi-linearly interpolates the four sample images corresponding to the subspace vertices. . . . . . . . . . . . . 17 2.5 Data collection setup . . . . . . . . . . . . . . . . . . . . . . . 19 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. . . . . . . . . . . . . . . . . . . . . . . . 20 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. . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.8 Comparison of simulated needle images in phantom using dif- ferent superimposition methods with the real image. . . . . . 22 viii List of Figures 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: sim- ulated image in phantom; lower-right: real image in phan- tom. Configurations shown in the subfigure captions are of the form (r, l, ρ, α, β). . . . . . . . . . . . . . . . . . . . . . . 24 2.10 The needle in water exhibited the bright tail. The needle in phantom did not show bright tail pattern in this experiment. 25 2.11 Comparison of simulated transverse needle image and in vivo needle image from a brachytherapy session. . . . . . . . . . . 26 2.12 Coordinate assignment for the sagittal image simulation (not to physical scale) . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.13 User interface showing different interpolation methods used for image rotation . . . . . . . . . . . . . . . . . . . . . . . . 31 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. Con- figurations shown in the subfigure captions are of the form (x, y, z, ρ, α, β). . . . . . . . . . . . . . . . . . . . . . . . . . . 32 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. Con- figurations shown in the subfigure captions are of the form (x, y, z, ρ, α, β). . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.15 Phantom with a seed mounted on a needle . . . . . . . . . . . 34 2.16 Comparison of seed images in the transverse plane . . . . . . 35 2.17 Comparison of simulated and real seed images. . . . . . . . . 36 2.18 Needle image showing bright tail in a scatterer-free phantom 37 2.19 Needle image in water with suspended scatterers. . . . . . . . 38 2.20 Needle in a PVC phantom. . . . . . . . . . . . . . . . . . . . 39 2.21 The load of the phantom may exert torque on the probe and change the pitch angle. . . . . . . . . . . . . . . . . . . . . . . 39 2.22 Phantom fracture makes bright tail disappear. . . . . . . . . 40 2.23 Needle in a pig liver; the needle location is pointed out. . . . 41 3.1 Sagittal B-mode images of needle in water . . . . . . . . . . . 44 ix List of Figures 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.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. . . . . . . . . . . . . . . . . . . . 47 3.4 Set-up for metal plate reflection study . . . . . . . . . . . . . 48 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.6 Comparison of the theoretical and experimental results when we assume the RF amplitude is proportional to the pressure. 51 3.7 Comparison of the theoretical and experimental results when we assume the RF amplitude is proportional to the sound intensity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.8 Band pattern in the quarter-wave layer when scanning air. . . 53 4.1 Band pattern visible in water, phantom and liver. . . . . . . . 56 4.2 Scattering by a spherical shell. For simplicity,only three layers are drawn here. All layers are elastic solids. . . . . . . . . . 57 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. . . . . . . . . . . . . . . . 63 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. . . . . . . . . . . . . . . . . . 64 4.5 Simulated impulse response corresponding to the back-scattered stress from a radius a = 1 mm solid sphere measured at a dis- tance r = 60 cm from the sphere. The y-axis has the unit of traction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.6 The simulated US pulse and power spectrum. . . . . . . . . . 66 4.7 The experimentally measured power spectrum of the received RF signal (reprinted from [18]). . . . . . . . . . . . . . . . . . 67 4.8 The simulated stress resulting from the simulated US pulse with the otherwise same parameters as in Fig. 4.5 . . . . . . . 68 4.9 The transient pulse response from a spherical shell with a = 1 mm, shell thickness 0.02mm, measured at r = 60 mm. . . . 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. . . . . . . . . . . . . . . . . . . . . 70 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. . . . . . . . . . . . . . . . . . . . 71 4.12 The comparison between a cylindrical shell and a solid. . . . 71 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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 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. . . . . . . . . . . . . . . . . . . . . . . 73 4.15 The Bessel functions of the first and second kind at ka = 5: jn(5) and nn(5). . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.1 Flow chart of the needle insertion simulation . . . . . . . . . 77 5.2 The global frame, image frame, and probe frame. Compare to Fig. 2.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3 Calculation of the pitch (α) and yaw (β) angles from the probe and needle frame . . . . . . . . . . . . . . . . . . . . . 82 5.4 Frames involved in the θ correction . . . . . . . . . . . . . . . 82 5.5 The geometry of the needle rotating around the robot wrist rotation centre. . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.6 The needle insertion simulator graphics interface (courtesy of Orcun Goksel). The needle is pointed out. . . . . . . . . . . . 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 pa- tience and understanding when ideas and results took weeks, even months to take shape. The joy of sharing good results will be remembered and cher- ished. 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 sup- port 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 Mat- lab 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 writ- ing 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, museum- visits 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 in- terventional 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]. Simulation- based 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 pro- cedures 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 re- views on several topics are given. 1.2 Background and literature review 1.2.1 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 sim- ulate 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 lon- gitudinal 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 sim- ulations requiring real-time performance. 5 1.2. Background and literature review 2. Pseudo-acoustic methods: the previous works on needle insertion sim- ulation 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 interpo- lating 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 interpo- lation 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: Lead pellets in the liver producing bright tail (reprinted from [45]). 1.2.6 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 particu- larly 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 inter- pret 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 spher- ical 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 solu- tions 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 1.2. Background and literature review A ut ho r St ea dy -s ta te /t ra ns ie nt G eo m et ry M et ho d Fa ra n [2 3] st ea dy -s ta te sp he ri ca l & cy lin dr ic al so lid th eo ry & ex pe ri m en t H am pt on et al [3 2] tr an si en t & st ea dy -s ta te sp he ri ca l so lid ex pe ri m en t B ar na rd et al [9 ] tr an si en t & st ea dy -s ta te cy lin dr ic al so lid & sh el l ex pe ri m en t H ic kl in g [3 6] tr an si en t & st ea dy -s ta te sp he ri ca l so lid th eo ry A yr es et al [6 ] st ea dy -s ta te sp he ri ca l sh el l th eo ry St an to n [7 1] st ea dy -s ta te fin it e cy lin de r th eo ry G au na ur d et al [2 4] tr an si en t sp he ri ca l sh el l th eo ry H as eg aw a et al [3 5] st ea dy -s ta te sp he ri ca l & cy lin dr ic al sh el l th eo ry A lle n et al [2 ] st ea dy -s ta te sp he ri ca l sh el l th eo ry Y e [7 8] st ea dy -s ta te sp he ri ca l sh el l th eo ry H u et al [3 7] st ea dy -s ta te sp he ri ca l sh el l th eo ry B ag ha ni [7 ] st ea dy -s ta te sp he ri ca l sh el l th eo ry H ar pe n [3 4] tr an si en t sp he ri ca l so lid th eo ry T ab le 1. 1: C om pa ri so n of pr ev io us w or ks on sc at te ri ng fr om sp he ri ca l an d cy lin dr ic al so lid an d sh el l. T he fir st gr ou p st ud ie s th e ge ne ri c sc at te ri ng ;t he se co nd gr ou p st ud ie s th e sc at te ri ng in th e co nt ex t of ul tr as ou nd co nt ra st ag en ts ; th e th ir d us es sc at te ri ng to m od el th e co m et ta il ar ti fa ct . 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 differ- ent 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 tech- nique 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 Transverse needle image simulation and tip visibility study 2.1.1 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 (a) The nee- dle tip in the water body {tail (b) The nee- dle shaft in the water (c) The nee- dle tip in the phantom (d) The nee- dle 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 configura- tions. 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 {P˜n, Cn} l Yn Xn Zn Y X Z r θ {O˜p, Cp} ρ Imaging plane TRUS probe β α bevel tip Y X Z Zn 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 {O˜ p, Cp}to the intersection of this imaging plane and the probe axis, and another Cartesian frame {P˜n, Cn} to the intersection of the imaging plane and theneedle 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 {O˜ p, Cp},and l to denote the needle depth. We describe the rotation transformation between Cp and Cn 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 lin- ear 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 fol- lowing 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 interpo- lation 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 ba- sis 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 orig- inal 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 pi 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: (pif)(x) = f̂(x) = ∑ i f(xi)li(x), pi : f → f̂ As special cases, when n = 0, li(x) are the shifted “box” functions and pi corresponds to the nearest neighbour interpolation operation; when n = 1, li(x) are the shifted “hat” functions and pi corresponds to the linear interpolation operation. We can show that the error bound decreases with decreasing h and increasing n: ||f−pif ||∞ ≤ 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 l1i (x) and l 2 j (y), given that the other dimension does not change: (pi1g)(x, y) := M∑ i=0 g(xi, y)l1i (x) (pi2g)(x, y) := N∑ j=0 g(x, yj)l2j (y) Now we can define a new projection by combining the two partial projections 15 2.1. Transverse needle image simulation and tip visibility study pi1 and pi2 as follows: (pi1pi2g)(x, y) := M∑ i=0 N∑ j=0 g(xi, yj)l1i (x)l 2 j (y) which has the property: (pi1pi2g)(x, y) = (pi2pi1g)(x, y) (2.2) We can show that the function space spanned by {l1i (x)l2j (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 ||g − pi1pi2g||∞ ≤ C(n) max (∣∣∣∣∣∣∣∣∂n+1g∂xn+1 ∣∣∣∣∣∣∣∣ ∞ hn+1, ∣∣∣∣∣∣∣∣∂m+1g∂ym+1 ∣∣∣∣∣∣∣∣ ∞ hm+1 ) (2.3) 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 rb lcld ρe ρf (r, l, ρ) (ra, lc, ρe)(ra, ld, ρe) (rb, lc, ρe) (ra, l, ρ) (ra, l, ρ) ld lc ρe ρf Figure 2.4: Visualization of the interpolation process. Here we only draw 3 dimensions. Every vertex represents a sample image. The red cross rep- resents 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 interpolation to get the interpolant image Ĩ. This is just bi-linearly 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 Ĩ has polar coordinates (ra, 0) in the imaging plane. Translate and rotate Ĩ 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 contin- uous 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 (Ul- trasonix 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 dur- ing 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 TRUS probe Stepper Needle Robot Wrist Water tank Figure 2.5: Data collection setup Variable Range # of samples r 3 ∼ 6.5cm 8 l 1 ∼ 22mm 12 ρ 0 ∼ 180◦ 3 α 0 ∼ 9◦ 4 β 0 ∼ 15◦ 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 TRUS probeNeedle Imaging Plane Water Tank 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 reg- istration 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 gela- tine 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 configura- tions. 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 back- ground 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) Superimposi- tion by addition. (b) Superimposi- tion by max. (c) Real Figure 2.8: Comparison of simulated needle images in phantom using dif- ferent 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 at- tainable needle configurations during brachytherapy, it could be used to study the needle visibility in the transverse plane. By comparing the rela- tive 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 (b) (6.5, 7.5, 0, 2, 0) Sim Real (c) (4.8, 6, 0, 0, 0) Sim Real (d) (4.5, 6, 135, 0, 0) Sim Real (e) (5.2, 7.5, 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 wa- ter (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 differ- ences 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, Cn) : R6 → R2, or more explicitly In(x, y, z, ρ, α, β). The lparameter in the figure is used for simulating the insertion, which is discussed in Section 2.2.2. For illustration purpose, we draw P˜1 and P˜2 as the pointswhere the needle axis enters and exits the imaging plane; 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 in- tensity 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 im- age 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, Cn} YnXn Zn Y X Z{O˜p, Cp} ρ Imaging plane TRUS probe P˜0 P˜1 P˜2 β α Y X Z Zn x l 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 dif- ferent “projected” pitch angles – the angle between the line segment and the horizontal image axis of the image plane: α′ = 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 Range # of samples x -0.7 ∼ 1.5cm 23 ρ 0 ∼ 180◦ 3 α 0 ∼ 6◦ 3 β 0 ∼ 15◦ 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 α′ to rotate the image around the new tip position by −α′. 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 coordi- nate: (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 in- terpolate to get the interpolant Ĩ ; 4. Translate Ĩ so that the tip is at (y, z) and then rotate the image around the tip by α′. 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, Cn}, the user could initiate aninsertion 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 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. 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 visi- bility. 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 Range # of samples x -10 ∼ 10mm 22 α 0 ∼ 9◦ 4 β 0 ∼ 15◦ 4 (a) Sagittal sample range and numbers Variable Range # of samples r 2.5 ∼ 5.5cm 7 l 1 ∼ 17mm 17 α 0 ∼ 9◦ 4 β 0 ∼ 15◦ 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 phan- tom 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 sus- pended scatterers. (b) Dark shadow in water with sus- pended 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 us- ing freehand so that we can freely adjust the angle. The phantom has similar elasticity and scatterer distribution with the gelatine phan- tom 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 exper- iments, 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. α Phantom TRUS probe Needle τ 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 phan- tom, before fracture (b) The nee- dle 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 ap- pearance 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 ultra- sound 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 in- cident 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 inci- dence, 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 under- standing 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 avail- able. 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 re- fraction 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 ap- proximations 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 equa- tions 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. Specif- ically, 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 Ap- pendix A for a short discussion. 3.1 A description of the phenomenon When imaging metal plates with nearly perpendicular incident angles, repet- itive 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” interchange- ably. Because the plate thickness cannot be neglected compared with the wave- length, 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 orthog- onal 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 + V1e−ikz T−n = W1eikα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 V1 = ρsα− ρfc ρsα+ ρfc W1 = 1 + V1 (3.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 V2 = ρfc− ρsα ρsα+ ρfc W2 = 1 + V2 (3.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 1 W1 V1 W1W2 W1V2 W1V2W2 W1V 2 2 solid fluid fluid 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 Support Plate being scanned 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; satu- ration 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 0 500 1000 1500 2000 2500 3000 3500 −1.5 −1 −0.5 0 0.5 1 1.5 x 104 RF: aluminum plate t In te ns ity (a) Aluminum plate 0 500 1000 1500 2000 2500 3000 3500 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 x 104 RF: steel plate t In te ns ity (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 follow- ing pressure-voltage relationship [17]: pi = eiiEi − cii where p is the pressure in the material, eii is the piezoelectric stress coef- ficient, 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 re- flection pressure with amplitudes ps1 = V1, ps2 = W1V1W2, and so on, as in Fig. 3.3: Hypothesis 1: ∫ 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: ∫ |ps|2δ(t)dt ∝ ∫ |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. Hy- pothesis 2 is thus supported by the experimental evidence. 50 3.3. Validation of the model 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 number of reverbs re la tiv e in te gr at ed in te ns ity in e ac h ba nd Comparison of results: reflections from an aluminum plate 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 (a) Aluminum plate 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 number of reverbs re la tiv e in te gr at ed in te ns ity in e ac h ba nd Comparison of results: reflections from a steel plate theory experiment (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 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 number of reverbs re la tiv e in te gr at ed in te ns ity in e ac h ba nd Comparison of results: reflections from an aluminum plate 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 (a) Aluminum plate 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 number of reverbs re la tiv e in te gr at ed in te ns ity in e ac h ba nd Comparison of results: reflections from a steel plate theory experiment (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 Discussion 3.4.1 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. In- tuitively, 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 displace- ment 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 1The 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 nee- dle 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 formu- lation 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 back- scattered 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 Φ(0)inc Φ(0)s ,Ψ (0) s Φ(1),Ψ(1) Φ(2),Ψ(2) θ r 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 sec- tion 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: ∇ · ∇Φ + k2αΦ = 0 ∇ · ∇~Ψ + k2β~Ψ = 0 (4.1) where kα = ωα , kβ = ω β are the longitudinal and shear wave numbers. Φ and ~Ψ are respectively the longitudinal and shear components of the displace- ment potential ~U (the gradient of ~U is the displacement.): ~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: Φ(m) = ∑ n [jcΦ,(m)n jn(k (m) α r) + ncΦ,(m)n nn(k (m) α r)]Pn(cos θ) Ψ(m) = ∑ n [jcΨ,(m)n jn(k (m) β r) + ncΨ,(m)n nn(k (m) β r)]Pn(cos θ) (4.4) 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 sub- spaces of these general solutions are allowed. We analyze this layer by layer: 2The 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 scat- tered wave. We assume the incident displacement potential to have unit amplitude and to be irrotational (no Ψ). We see these assump- tions are satisfied when Φ has the form of Eq. A.6: Φ(0)inc = e ik (0) α z We can expand it in spherical coordinates [56]: Φ(0)inc = e ik (0) α r cos θ = ∑ n (2n+ 1)injn(k(0)α r)Pn(cos θ) (4.5) Comparing to Eq. 4.4, it can be seen jcΦ,(0)n = (2n + 1)in, and the Neumann function does not contribute to this solution. Ψ(0)inc = 0 by definition. The scattered solution is to the outward direction. Only a certain com- bination 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) In this case, we can use a single coefficient hcΦ,(0)n for Φ (0) s , and hc Ψ,(0) n for Ψ(0)s . 2. For the middle layers, both the inward and outward solutions are al- lowed. 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 We therefore have 2 unknowns in Φ(0)s ,Ψ (0) 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: U (i)r,n = U (i+1) r,n U (i) θ,n = U (i+1) θ,n (4.6) 2. The normal stress must be continuous across the boundaries for each mode: σ(i)rr,n = σ (i+1) rr,n σ (i) rθ,n = σ (i+1) 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 ob- tained 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 θ = pi 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 po- tential coefficients defined in Eq. 4.4 [7]. Omitting the subscript “scattered” and the superscript (0): σrr,n(ω0t) =σrr,n e−iω0t = [ hcΦn ( λTΦ1 + 2µT Φ 2 ) + hcΨn ( λTΨ1 + 2µT Ψ 2 )] Pn(cos θ) e−iω0t (4.8) where TΦ1 =− k2α h[1]n (kαr) TΦ2 =( n(n− 1) r2 − k2α) h[1]n (kαr) + 2kα r h [1] n+1(kαr) TΨ1 =− 2 (n− 1)n(n+ 1) r2 h[1]n (kβr) TΨ2 =− (n− 1)n(n+ 1) r2 h[1]n (kβr) + kβn(n+ 1) r h [1] n+1(kβr) The coefficients hcΦn , hcΨn , as shown in Section 4.3, can be solved from a system of linear equations 4.6, 4.7 obtained from the boundary conditions. Furthermore, the shear stress components σ(0)rθ , σ (0) 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), 3This is verified by evaluating σrθ and comparing with σrr – the former indeed very small. Legendre polynomials and derivatives are calculated using the “lpn” routine con- verted 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ω0t (with separated spatial and time parts), the Fourier transform is: F{f(x)e−iω0t} = f(x)2piδ(ω + ω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(ω0t) = eikz−iω0t = e−iω0t Applying Eq. 4.9: F{Φinc(ω0t)} = 2piδ(ω + ω0) Similarly for the output (Eq. 4.8): F{σrr,n(ω0t)} = 2pi σ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, corre- sponding 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 argu- ment to the eigen functions in Eq. 4.4. By noting that ka = 2pia/λ, 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. 4We 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 0 1 2 3 4 5 6 0 0.5 1 1.5 2 x 1016 Frequency (MHz) Am pl itu de 5 10 15 20 25 0 0.5 1 1.5 2 x 1016 ka Magnitude, θ=180, r =60a 0 1 2 3 4 5 6 −300 −200 −100 0 Frequency (MHz) R ad ia ns Phase, θ=180, r =60a scattered stress 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 com- ponents5. 5Curiously, 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 Phantom Steel ρ (kg/m3) 1000 7700 λ 2.25e9 9.695e10 µ 3e4 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 modula- tion 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 phan- tom 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 30 35 40 45 50 55 60 65 70 −3 −2 −1 0 1 2 3 4 x 1015 Real part of the impulse respone, θ=180, r =60a Time (µ s) Sc al ed tr ac tio n (N /m 2 ) 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, t) = sin(ω(t− z c ))e−α(t− z c ) 6This lt corresponds to roughly 25 mm in the B-mode. 65 4.5. Transient simulation 0 0.5 1 1.5 2 2.5 3 3.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Time (µ s) D is pl ac em en t Decaying sinusoid (a) The simulated ultrasound pulse. 0 1 2 3 4 5 6 x 106 0 1 2 3 4 5 6 7 Spectrum density Hz Po w er (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)∂z at 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 interpreta- tion 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 sim- ulation 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 30 35 40 45 50 55 60 65 70 −4 −3 −2 −1 0 1 2 3 4 x 1015 Real part of the pulse respone, θ=180, r =60a Time (µ s) Sc al ed tr ac tio n (N /m 2 ) 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 30 35 40 45 50 55 60 65 70 −6 −4 −2 0 2 4 6 x 1015 Real part of the pulse respone, θ=180, r =60a Time (µ s) Sc al ed tr ac tio n (N /m 2 ) 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 1000 1500 2000 2500 3000 −8000 −6000 −4000 −2000 0 2000 4000 6000 8000 RF: needle shaft t In te ns ity 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 contam- inated 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 pi 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 0 5 10 15 20 25 30 35 −5 −4 −3 −2 −1 0 1 2 3 4 5 x 1017 Time (µ s) Sc al ed tr ac tio n (N /m 2 ) Real part of the pulse respone, θ=0, r =1a scattered incident sum 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 0 50 100 150 200 250 300 350 400 0 0.5 1 1.5 2 2.5 3 Number of terms Am pl itu de Convergence of the incident series r=2mm r=4mm r=6mm r=8mm 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 Numerical issues 4.7.1 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 scat- tering 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 16 16 16 48 48 48 80 80 80 11 2 11 2 11 2 14 4 14 4 14 4 17 6 17 6 17 6 208 20 8 20 8 20 8 240 24 0 24 0 24 0 272 27 2 27 2 27 2 304 30 4 30 4 n ka 20 40 60 80 100 120 140 5 10 15 20 25 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) rr,n,scattered − σ(1)rr,n = −σ(0)rr,n,incident (4.10) As shown in Eq. 4.8, σ(0)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 the solid sphere (only 1 layer); σ(1)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 0 10 20 30 40 50 10−50 100 1050 n ka =5 Bessel first kind Bessel second kind (Neumann) 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 7This 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 interro- gation. 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 pre- vious 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 consid- ering 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 Chap- ter 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 subse- quently converted to C and re-implemented in OpenGL for real-time per- formance8. 5.1 Notations and overview Some notation conventions are as follows: we measure every frame with respect to a base frame {O˜ 0, C0}, where O˜ 0 is the origin and C0 is theidentity. Without superscripts, the points and vectors are understood to be 8The 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 Needle- tissue interaction model Coordinate transform Line- crossing Interpolation Matlab 5D database Repacking Sample prepa- ration, registration Needle points in global frame Needle points in probe frame depth calcula- tion local configu- ration 1D database Interpolated image Shift and rotate Display Superimpose background Needle image θ-correction Offline 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 thatpoints from P˜1 to P˜2 as P1P2, right-handed coordinate frames as {O˜ 1, C1},etc., where C1 is an orthonormal matrix. We write points measured in C1 as P˜11, P˜12, etc., vectors in C1 as P 11P 12 , etc. With this convention, we canwrite vectors measured in the base frame using its representation in C1: P1P2 = C1(P˜12 − P˜11) = C1 · (P 11P 12 ).9The sense of rotation 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 {O˜ t, C0}; 2. The imaging plane coordinate system {O˜ i, Cp} with one axis definingthe norm and the other two the image edge directions; O˜ i is locatedat the bottom centre of the image; 9A 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 {O˜ t, C0} x y z X Y Z {O˜ i, Cp} O˜p P˜nP˜ i P˜ i+1 Figure 5.2: The global frame, image frame, and probe frame. Compare to Fig. 2.2. 3. The probe coordinate system {O˜ p, Cp} located at the probe centre andwith the same orientation as that of the image; 4. The needle/robot coordinate system {P˜n, Cn}. The image simulation algorithm developed in Chapter 2 (the “interpola- tion” block) requires as input {P˜n, Cn} measured with respect to {O˜ p, Cp},that is, the relative position and orientation 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 {O˜ t, C0}. We therefore need twosteps to calculate the local needle configuration from the geometry given by the needle-tissue interaction model: 1. Transform the points description P˜ i, i = 1, . . . , N from the globalframe 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 {O˜ p, Cp}(section 5.4). 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 θ (sec- tion 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 con- vert 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 op- erates on two 1D arrays: one with continuous storage of the sample images, the other with flags indicating the empty images. 5.3 Transformation from the global coordinate to the probe coordinate 5.3.1 Problem description The problem is to convert a point P˜n from the global coordinate (hereas the base frame) to the probe coordinate. What we know is the image frame {O˜ i, Cp} measured in the global frame {O˜ t, C0}, and the vector OiOpmeasured in the image frame; this vector lies in the XY plane and can 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 C0 the base frame: OpPn = P˜n −O˜ p = P˜n − (O˜ i +OpOi) = P˜n − (O˜ i + Cp(λ~e1 + γ~e2)) Using the measurement in Cp: OpPn = CpP˜pn 79 5.4. Line crossing Equating the above two, P˜pn = CTp [P˜n − (O˜ i + Cp(λ~e1 + γ~e2))] = CTp (P˜n −O˜ i)− (λ~e1 + γ~e2) 5.4 Line crossing 5.4.1 Problem description Given a series of points P˜ i, i = 1, . . . , N , and the imaging plane {O˜ i, Cp},calculate the intersection P˜n between the line segment connecting the Npoints and the plane. Note that the problem itself is frame independent, 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: 1. ~Z is perpendicular to the vector from the probe origin to the intersec- tion: (P˜n −O˜ p) · ~Z = 0 (5.1) 2. P˜ i, P˜n, P˜ i+1 lie on a line: P˜n − P˜ i = c (P˜ i+1 − P˜n) (5.2) From 5.2: P˜n = P˜ i + cP˜ i+11 + c (5.3) Substituting into Eq. 5.1, we can solve for c: c = (O˜ p − P˜ i) · ~Z (P˜ i+1 −O˜ p) · ~Z (5.4) the sign of which determines if PiPi+1 crosses the imaging plane. Note that c = −1 ⇐⇒ P˜ i = P˜ i+1, which is unrealistic. Also, (P˜ i+1 − O˜ p) · ~Z = 0means P˜ i+1 is in the imaging plane.The algorithm to determine the crossing point is therefore: 80 5.5. θ-correction 1. Start from P˜0 from the point list; 2. If (P˜ i+1 −O˜ p) · ~Z = 0, P˜ i+1 is the crossing point; return; 3. Calculate c using Eq. 5.4; 4. If c ≥ 0, PiPi+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. The above algorithm only gives the position of P˜n; the orientation canbe obtained by referring to Fig. 5.3. In the figure, ~Zn = PiPi+1 |PiPi+1| (~Zn · ~Z)~Z = ~Znz (~Zn · ~X) ~X = ~Znx ~Znz + ~Znx = ~Znxz we have, α = 2 sign(~Zn · ~Y ) atan ||Ẑn − Ẑnxz||||Ẑn + Ẑnxz|| β = sign(~Zn · ~X) atan | ~Znx| |~Znz| 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 θ-correction 5.5.1 Problem description As demonstrated in section 2.1.2, a needle image at an angle θ can be ob- tained 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 Cn) measured in Cp; we want to know ρ, α, β with respect to the rotated RθCp, so that we can retrieve the sample image (taken at θ = 0) with the correct configurations. 10We have used the more accurate atan of half angle instead of the less accurate acos [44, 66]. 81 5.5. θ-correction Y X Z β α ~Zn ρ ~Znz ~Znx ~Znxz Figure 5.3: Calculation of the pitch (α) and yaw (β) angles from the probe and needle frame θ {O˜p, Cp} {P˜n, Cn} C Y X Z 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 ρ, α, β (correspond- ing to C = R−θCn) with respect to Cp. The procedure is as follows: 1. Compute Cn the rotation matrix from ρ0, α0, β0; this is straight-forward: Cn = RY (β)RX(α)RZ(ρ) where RZ(ρ) means rotating around Z by ρ. 2. Rotate Cn: C = RZ(−θ)Cn 3. Calculate the new ρ, α, β from the element in C [48]: ρ = atan(r32/r33) α = atan(−r31/ √ r232 + r 2 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 con- figuration. In Fig. 5.5, after rotating the needle around O˜ by pitch and yawangles α and β, the new depth is l′. The distance between the imaging plane and the robot motion plane is di. During the data collection, l is the parameter that we manipulate; in the simulation, given l′ (or l′⊥), we wish 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 l′. 5.6.2 Solution To derive this relationship, we need to take into account how l was manipu- lated 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 di h Imaging Plane d l l′⊥ l′ O˜ ~s ~s′ 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 ~s′ = cosβ 0 sinβ0 1 0 − sinβ 0 cosβ 1 0 00 cosα − sinα 0 sinα cosα 0h d Extracting the third dimension: l′⊥ + di = d cosα cosβ + h sinα cosβ Also we have: l′⊥ = l ′ 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 + l′ cosα cosβ − h sinα cosβ cosα cosβ − di When d is a constant, that is, when l is changed by moving the probe, we have: l = (1− cosα cosβ)d+ l′ cosα cosβ − h sinα cosβ From these two relations, we can obtain l from l′ 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 Chap- ter 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 build- ing 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 help- ful 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 physically- based 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 tech- nique 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 “canon- ical” geometries – the plates and the spheres in our case. More complex techniques exist for modelling special geometries but require consider- ably 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 gener- ate 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 sim- ulations. 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 vis- ibility and bright tail and dark shadow artifacts are studied experi- mentally 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 Sec- tion 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 near- est 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 simula- tor: 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 simu- lation 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 di- rection. 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 charac- teristic 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 geome- tries: for non-standard geometries such as the needle tip, the simu- lation 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 scatter- ing 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 Uni- versity, 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. Inter- national journal of solids and structures, 23(7):937–946, 1987. [7] Ali Baghani. A simulation software for the analysis of ultrasound in- teraction with multi-layer spherical scatterers. Technical report, Uni- versity 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. Ul- trasonix 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 fluid- loaded 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 resonancescat- tering 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 needle- tissue 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, Uni- versity 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 acous- tics. Computer Methods in Applied Mechanics and Engineering, 195(13- 16):1594–1607, 2006. [34] M.D. Harpen. Partial wave analysis of the ultrasound comet tail arti- fact. 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 scatter- ing from a double-layered micro shell for encapsulated drug delivery. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Con- trol, 51(7):808–820, 2004. [38] J. Huang, J.K. Triedman, N.V. Vasilyev, Y. Suematsu, R.O. Cleve- land, 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 COMPUT- ING, 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. Ultra- sound Quarterly, 17(4):201, 2001. [46] RM Kirberger. Imaging artifacts in diagnostic ultrasound: a review. Veterinary radiology and ultrasound: the official journal of the Amer- ican 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 tech- niques, pages 31–42. ACM New York, NY, USA, 1996. [50] D. Magee, Y. Zhu, R. Ratnalingam, P. Gardner, and D. Kessel. An aug- mented reality simulator for ultrasound guided needle placement train- ing. 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 ul- trasonic 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. Ul- trasound in Obstetrics and Gynecology, 24(5):581–585, 2004. [55] FG Mitri, BJ Davis, JF Greenleaf, and M. Fatemi. In vitro compara- tive study of vibro-acoustography versus pulse-echo ultrasound in imag- ing permanent prostate brachytherapy seeds. Ultrasonics, 49(1):31–38, 2009. [56] P.M.C. Morse and K.U. Ingard. Theoretical Acoustics. Princeton Uni- versity 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 Ul- trasonographic 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). Mi- crosoft Press Redmond, WA, USA, 1999. [64] R. Rohling, A. Gee, and L. Berman. A comparison of freehand three- dimensional 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 Au- tomation, 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 ul- trasound 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 resam- pling. 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 ultra- sound guided needle puncture using patient specific data with 3D tex- tures and volume haptics. Computer Animation and Virtual Worlds, 19(2), 2008. [77] E. W. Weisstein. Bessel function of the second kind. Math- World – 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 renderingrepresenta- tion, sampling and compression. Signal Processing: Image Communi- cation, 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 coef- ficients with orthogonal incidence from the fluid to the solid. We roughly follow the notation in [39]. A.1 Solid In the solid: after decomposing the displacement field ~U using the Helmholtz 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∂ 2~U ∂t2 is equivalent to the wave equations: ∂2Φ ∂t2 = α2∇ · ∇Φ ∂2~Ψ ∂t2 = β2∇ · ∇~Ψ where σ is the stress tensor, α = √ λ+2µ ρs is the longitudinal wave speed, β = √ µ ρ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: ∇ · ∇Φ + k2αΦ = 0 ∇ · ∇~Ψ + k2β~Ψ = 0 (A.2) where kα = ωα , kβ = ω β are the longitudinal and shear wave numbers. A.2 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 con- dition and the equation of motion; the result is: ∇ · ∇p+ k2p = 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: ρfω 2~U · ~n = ∇p · ~n (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) 11Time 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 pinc ps T−n Solid Fluid 0 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(kxx+kyy+kzz) 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Ψ1CΨ2 CΨ3 eikβz ps = V e−ikz where V is called the reflection coefficient. Applying Eq. A.1 again, we get the displacement: ~U(z) = −CΨ2ikβeikβzCΨ1ikβeikβz CΦikαe ikαz (A.7) 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 ∂z + λ( ∂Ux ∂x + ∂Uy ∂y + ∂Uz ∂z ) = −(2µ+ λ)CΦk2α eikαz (A.8) This means that the transmitted stress does not contain the shear compo- nent – “mode conversion” does not occur at orthogonal incidence. Sub- stituting Eq. A.7 A.8 and the solutions to the wave equations to Eq. A.4 A.5: { ρfω 2CΦkαe ikαz = −kV e−ikz + keikz (2µ+ λ)CΦk2α e ikαz = eikz + V e−ikz at z = 0 From the solutions we have: ps = V e−ikz = 2µ+ λ− αcρf 2µ+ λ+ αcρf e−ikz = ρsα− ρfc ρsα+ ρfc e−ikz T−n = −σ33 = (1 + V )eikαz = Weikα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 be- cause 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 AcquireImg NextPlan OnAcquisitionComplete all motors reaching de- sired 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, “Ac- quireImg” posts “AcquisitionComplete” message which results in a call to “NextPlan”, closing the loop. The needle is moved to the next pre-plan po- sition 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