Applied Science, Faculty of
Biomedical Engineering, School of
DSpace
UBCV
Ai, Min
Doctor of Philosophy - PhD
University of British Columbia
Photoacoustic tomography (PAT) is a hybrid imaging modality that combines optical excitation and acoustic detection and overcomes the optical diffusion limit in tissue. Taking the advantages of deep penetration and high resolution, PAT is widely applied in biosciences. In prostate imaging, several PAT systems have been demonstrated to be capable of showing distinct photoacoustic (PA) contrast from malignant tissue or angiogenesis-related to prostate cancer. In this thesis, I explore the potential of applying PAT with transurethral illumination for prostate imaging. Challenges in translating PAT system to clinical imaging for prostate include the lack of sufficient local fluence for deep tissue penetration, risk of over irradiation near the laser-tissue contact surface and limited image reconstruction quality caused by the limited detection view and noisy PA signal. In this thesis, systematic design, optimization, and application of a PAT system with transurethral illumination are conducted. A fiber coupling scheme with a beam homogenizer is demonstrated for coupling high energy pulses in a single multimode fiber, achieved by using a cross cylindrical lens array. The peak power on the fiber tip surface is reduced and thus enhances the coupling performance. Recorded high pulse energy is achieved with high coupling efficiency as well. The high pulse energy can enable deep imaging depth in tissue. With high-energy pulses delivered by the multimode fiber, a transurethral illumination probe is designed, which can illuminate the prostate from the urethra. A parabolic cylindrical mirror reflects the light emitted from the diffusing fiber to achieve a parallel side illumination with doubled fluence. This design is optimized for both high energy pulse illumination and maintaining the laser fluence below the safety limit in tissue. As the detection view for prostate imaging is limited, a variance-reduced stochastic gradient descent (VR-SGD) algorithm is developed to improve the image reconstruction. The algorithm is verified by simulation and experiment of 2D imaging. VR-SGD algorithm demonstrates its capability to reduce the noise level and artifacts generated in the limited-view detection by linear array transducers. Through this study, the PAT system with transurethral illumination is shown to be capable and feasible for prostate imaging.
PHOTOACOUSTIC TOMOGRAPHY FOR PROSTATE IMAGING BY TRANSURETHRAL ILLUMINATION by Min Ai B.Sc., Harbin Institute of Technology, 2011 M.Sc., Harbin Institute of Technology, 2013 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Biomedical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) March 2020 © Min Ai, 2020 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled: Photoacoustic tomography for prostate imaging by transurethral illumination submitted by Min Ai in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The faculty of graduate and postdoctoral studies (Biomedical Engineering) Examining Committee: Prof. Shuo Tang Supervisor Prof. Septimiu E. Salcudean Supervisory Committee Member Prof. Robert Rohling Supervisory Committee Member Prof. Jane Z. Wang University Examiner Prof. Keng Chou University Examiner Additional Supervisory Committee Members: Prof. Purang Abolmaesumi Supervisory Committee Member Supervisory Committee Member iii Abstract Photoacoustic tomography (PAT) is a hybrid imaging modality that combines optical excitation and acoustic detection and overcomes the optical diffusion limit in tissue. Taking the advantages of deep penetration and high resolution, PAT is widely applied in biosciences. In prostate imaging, several PAT systems have been demonstrated to be capable of showing distinct photoacoustic (PA) contrast from malignant tissue or angiogenesis-related to prostate cancer. In this thesis, I explore the potential of applying PAT with transurethral illumination for prostate imaging. Challenges in translating PAT system to clinical imaging for prostate include the lack of sufficient local fluence for deep tissue penetration, risk of over irradiation near the laser-tissue contact surface and limited image reconstruction quality caused by the limited detection view and noisy PA signal. In this thesis, systematic design, optimization, and application of a PAT system with transurethral illumination are conducted. A fiber coupling scheme with a beam homogenizer is demonstrated for coupling high energy pulses in a single multimode fiber, achieved by using a cross cylindrical lens array. The peak power on the fiber tip surface is reduced and thus enhances the coupling performance. Recorded high pulse energy is achieved with high coupling efficiency as well. The high pulse energy can enable deep imaging depth in tissue. With high-energy pulses delivered by the multimode fiber, a transurethral illumination probe is designed, which can illuminate the prostate from the urethra. A parabolic cylindrical mirror reflects the light emitted from the diffusing fiber to achieve a parallel side illumination with doubled fluence. This design is optimized for both high energy pulse illumination and maintaining the laser fluence below the safety limit in tissue. iv As the detection view for prostate imaging is limited, a variance-reduced stochastic gradient descent (VR-SGD) algorithm is developed to improve the image reconstruction. The algorithm is verified by simulation and experiment of 2D imaging. VR-SGD algorithm demonstrates its capability to reduce the noise level and artifacts generated in the limited-view detection by linear array transducers. Through this study, the PAT system with transurethral illumination is shown to be capable and feasible for prostate imaging. v Lay Summary PAT is a hybrid imaging modality which combines optical excitation and acoustic detection. Taking the advantages of deep penetration and high resolution, PAT is widely studied and applied for visualizing tissues and organs. In prostate imaging, the current PAT systems are limited by the lack of sufficient laser power for deep tissue penetration and the risk of over irradiation near the laser-tissue contact surface. Here a PAT system with transurethral illumination is developed. Optimization of the fiber delivery is conducted to couple high energy pulses. The transurethral illumination probe is built, achieving a parallel side illumination. Application of the PAT system is demonstrated in prostate mimicking phantom. To improve the image reconstruction, an iterative reconstruction algorithm based on variance-reduced stochastic gradient descent is also developed and validated. This study shows that PAT system with transurethral illumination has great potential for clinical application. vi Preface All the work presented henceforth was conducted in the Biophotonics Laboratory and Robotics and Control Laboratory at the University of British Columbia (Vancouver campus). A version of Chapter 2 has been published in the following papers: â¢ M. Ai, W. Shu, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, "High energy laser pulse coupling in a multimode fiber for photoacoustic tomography," Photonics West 2016, SPIE, 9708-4H (2016). â¢ M. Ai, W. Shu, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, âDesign of high energy laser pulse delivery in a multimode fiber for photoacoustic tomography,â Opt. Express 25(15), 17713â17726 (2017). I was the lead investigator, responsible for all major areas of concept formation, design, experiment, data collection and analysis, as well as manuscript composition. Weihang Shu was involved in the experiment a data collection. Dr. Tim Salcudean, Dr. Robert Rohling and Dr. Purang Abolmaesumi were involved in concept formation and manuscript composition. Dr. Shuo Tang was the supervisory author on this project and was involved throughout the project in concept formation and manuscript composition. A version of Chapter 3 has been published in the following papers: â¢ M. Ai, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, "Transurethral illumination probe design for deep photoacoustic imaging of prostate," Photonics West 2018, SPIE, 10494-4C (2018). vii â¢ M. Ai, J. Youn, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, âPhotoacoustic tomography for imaging the prostate: a transurethral illumination probe design and application,â Biomed. Opt. Express 10, 2588-2605 (2019). I was the lead investigator, responsible for all major areas of concept formation, design, experiment, data collection and analysis, as well as manuscript composition. Dr. Jong-in Youn was involved in the experiment and contributed to the development of prostate mimicking phantom. Dr. Tim Salcudean, Dr. Robert Rohling and Dr. Purang Abolmaesumi were involved in concept formation and manuscript composition. Dr. Shuo Tang was the supervisory author on this project and was involved throughout the project in concept formation and manuscript composition. A version of Chapter 4 has been published in the following papers: â¢ W. Shu, M. Ai, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, "Image registration for limited-view photoacoustic imaging using two linear array transducers," Photonics West 2015, SPIE, 9323-48 (2015). â¢ W. Shu, M. Ai, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, âBroadening the detection view of 2D photoacoustic tomography using two linear array transducers,â Opt. Express 24, 12755-12768 (2016). A manuscript based on Chapter 4 is under preparation: â¢ M. Ai, J. Cheng, D. Karimi, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, âVariance-reduced stochastic gradient descent algorithm for photoacoustic image reconstruction,â (in preparation). In the first-author journal manuscript, I was the lead investigator, responsible for all major areas of concept formation, algorithm development, experiment, data collection and analysis, as well as manuscript composition. Jiayi Cheng and Dr. Davood Karimi were involved in the concept viii formation and algorithm development. Dr. Tim Salcudean, Dr. Robert Rohling and Dr. Purang Abolmaesumi were involved in concept formation and manuscript composition. Dr. Shuo Tang was the supervisory author on this project and was involved throughout the project in concept formation and manuscript composition. In the second-author papers, I designed the optics setup for the PAT system and involved in the experiment and data collection. ix Table of Contents Abstract .......................................................................................................................................... iii Lay Summary .................................................................................................................................. v Preface ............................................................................................................................................ vi Table of Contents ............................................................................................................................ ix List of Tables ............................................................................................................................... xiii List of Figures ............................................................................................................................... xiv List of Abbreviations .................................................................................................................... xix Acknowledgements ...................................................................................................................... xxi Chapter 1: Introduction ................................................................................................................ 1 1.1 Introduction of photoacoustic imaging ........................................................................ 1 1.2 Basic principles of PAT ............................................................................................... 3 1.3 Implementation of PAT ............................................................................................... 4 1.4 PAT for prostate imaging ............................................................................................ 6 1.5 Objectives and contributions ....................................................................................... 9 1.6 Outline ....................................................................................................................... 11 Chapter 2: Design of high energy laser pulse delivery in a multimode fiber for photoacoustic tomography .................................................................................................................................. 13 2.1 Introduction ............................................................................................................... 13 2.2 Beam homogenization with a cross cylindrical lens array ........................................ 16 x 2.2.1 Principle ............................................................................................................. 16 2.2.2 Theoretical and experimental validation of the focal intensity distribution ...... 21 2.3 Coupling high energy pulses into a multimode fiber for PAT .................................. 25 2.3.1 Experimental setup ............................................................................................ 25 2.3.2 Performance of fiber coupling ........................................................................... 26 2.3.3 Preliminary PAT Imaging ................................................................................. 29 2.4 Summary .................................................................................................................... 31 Chapter 3: Photoacoustic tomography for imaging prostate: a transurethral illumination probe design and application ...................................................................................................... 33 3.1 Introduction ............................................................................................................... 33 3.2 Design of the transurethral illumination probe .......................................................... 37 3.2.1 Configuration of the probe ................................................................................ 37 3.2.2 Fabrication of the diffusing fiber ....................................................................... 39 3.2.3 Design of side reflection mirror to reshape the beam ........................................ 43 3.2.4 Characterization of the transurethral illumination probe ................................... 46 3.3 PAT imaging experimental results ............................................................................ 48 3.3.1 Experimental setup ............................................................................................ 48 3.3.2 Human hairs inside chicken breast tissue .......................................................... 49 3.3.3 Prostate mimicking phantom ............................................................................. 50 3.4 Discussion .................................................................................................................. 53 xi 3.5 Summary .................................................................................................................... 57 Chapter 4: Variance-reduced stochastic gradient descent algorithm for photoacoustic image reconstruction .............................................................................................................................. 59 4.1 Introduction ............................................................................................................... 59 4.2 Method ....................................................................................................................... 62 4.2.1 PAT imaging model ........................................................................................... 62 4.2.2 Variance-reduced stochastic gradient descent method for solving the cost function â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦â¦64 4.3 Numerical simulation ................................................................................................ 68 4.4 Experiments ............................................................................................................... 79 4.4.1 PAT imaging by a single linear array transducer .............................................. 80 4.4.2 PAT imaging by two linear array transducers ................................................... 85 4.5 Summary .................................................................................................................... 88 Chapter 5: Conclusions and future work .................................................................................. 90 5.1 Overview ................................................................................................................... 90 5.2 Summary .................................................................................................................... 90 5.3 Future work ................................................................................................................ 93 5.3.1 Improvement on the transurethral probe ........................................................... 93 5.3.2 Improvement on the reconstruction algorithm .................................................. 94 5.3.3 In vivo prostate imaging .................................................................................... 96 xii 5.4 Conclusion ................................................................................................................. 97 Bibliography .................................................................................................................................. 98 xiii List of Tables Table 1.1 Current PAT systems for in vivo prostate imaging ......................................................... 9 Table 2.1 Summary of PAT systems with fiber delivery .............................................................. 16 Table 2.2 Specifications of the cross cylindrical lens array .......................................................... 22 Table 4.1 Parameter settings for the simulation of VR-SGD algorithm ....................................... 69 xiv List of Figures Figure 2.1 Schematic of the coupling optics. The incident beam from the pulse laser is split into beamlets by a cross cylindrical lens array, and the beamlets are focused by a plano-convex lens at its focal plane. To couple the light into a fiber, the fiber input end is placed at the focal plane of the plano-convex lens. ................................................................................................................... 17 Figure 2.2 (a) Optical model of the coupling scheme with the lens array. The lens array transforms the incident plane wave into multiple spherical wavelets. Constructive interference happens where the propagation angle ð! satisfies ðð ððð! = ðð, which forms the mth diffraction order on the focal plane. (b) 1D plot of the periodic phase modulation function introduced by the lens array. p is the pitch of the lens array. .......................................................................................................... 19 Figure 2.3 The light intensity distribution map (left column) and a corresponding intensity line profile (right column) at the focal plane. (a) and (b) are obtained by theoretical analysis; (c) and (d) are obtained by Zemax simulation; (e) and (f) are obtained by experiment. The location for the line profile is marked by an arrow. ................................................................................................ 23 Figure 2.4 Measured light intensity distribution map (a) and line profile; (b) at the output of the multimode fiber. The location for the line profile is marked by an arrow. ................................... 25 Figure 2.5 Experimental setup of the photoacoustic tomography with fiber delivering high energy pulses. The laser beam from the optical parametric oscillator (OPO) is homogenized and coupled into a multimode fiber to deliver light to phantom. DAQ: data acquisition system. ..................... 26 Figure 2.6 Output power characterization at different wavelengths and fiber core sizes. (a) Averaged pulse energy coupled into a 1000-µm-core-diameter fiber (blue curve). The error bar indicates the energy fluctuation measured over one hundred consecutive pulses. The laser energy measured at the input of the fiber is shown in red curve. (b) Similar power characterization as in xv (a) but with the light coupled into a 1500-µm-core-diameter fiber. (c) Comparison of the coupling efficiency of the two types of fiber. ............................................................................................... 28 Figure 2.7 PAT image of printed dots array. (a) Imaging result without the cover of chicken breast and the effective imaged area in a 10-mm-diameter circle. (b) Imaging result with the cover of 5 mm thick chicken breast and the increased effective imaged area in a 30-mm-diameter circle. .. 31 Figure 2.8 In vivo PAT image of forearm. (a) Overlay of PAT image (colored) on top of ultrasound image (grey scale). The PAT image highlights the skin surface and locations of several blood vessels labeled with numbers. (b) Photograph of the forearm where the blood vessels imaged by PAT are also marked. The dash line shows the place of the cross-sectional image. ..................... 31 Figure 3.1 Design of the transurethral illumination probe. (a) Overview of the transurethral illumination probe. (b) The beam profile from the probe displayed by a 660 nm light. (c) Longitudinal cross-section view of the probe. (d) Axial cross-section view of the probe. ........... 39 Figure 3.2 Schematic illustration of the etching process (left) and photographs of the diffusing fiber (right). (a) Three steps of the etching process. (b) Rough surface of the diffusing fiber. (c) The diffusing fiber with 660 nm visible light coupled in. (d) Illumination pattern on a piece of white paper of the 660 nm light from the diffusing fiber. ............................................................. 41 Figure 3.3 Calculated fluence on a cylindrical illumination surface at various distances. The output energy from the fiber is assumed to be 50 mJ/pulse. The dash lines indicate the position of the fiber surface and the tissue surface, respectively. .................................................................................. 42 Figure 3.4 Comparison of the reflection pattern. (a) and (b) are the ray tracing for the circular cylindrical mirror and the parabolic cylindrical mirror, respectively. (c) and (d) are the simulated light distribution pattern at 5 mm away from the circular cylindrical mirror and the parabolic cylindrical mirror, respectively. ..................................................................................................... 43 xvi Figure 3.5 Experimental results of the beam profile. (a) Illustration of the measurement setup. (b) Changes of the beam profile at various distances for the circular cylindrical mirror and the parabolic cylindrical mirror, respectively. (c) and (d) are the comparison of the FWHM of the beam profile in the x and y directions, respectively. ..................................................................... 46 Figure 3.6 Beam profile from the parabolic cylindrical mirror with and without the 45° end mirror. (a) 2D beam profile without the end mirror (b) 2D beam profile with the end mirror (c) 1D beam profile along x axis. (d) 1D beam profile along y axis. White arrows indicate the positions where the 1D profiles are extracted. ......................................................................................................... 47 Figure 3.7 PAT images of human hair embedded in chicken breast tissue. (a) Without attenuation compensation. (b) With attenuation compensation. ...................................................................... 50 Figure 3.8 PAT imaging of prostate mimicking phantom. (a) Photograph of the phantom where two tube shaped cavities indicate the location of rectum and urethra. (b) Experimental setup. (c) and (d) are the PAT images of the pencil leads embedded in the phantom without, and with attenuation compensation, respectively. ........................................................................................ 52 Figure 3.9 PAT images of hair embedded in the prostate mimicking phantom. (a) Forward illumination where the illumination probe faces the TRUS probe. (b) Backward illumination where the illumination probe faces the opposite side of the TRUS probe. (c) Result combining (a) and (b). .................................................................................................................................................. 53 Figure 4.1 Imaging configuration by using various detection geometry with equivalent number of transducers. (a) Linear array detection. (b) Parallel-shape detection. (c) L-shape detection. (d) Square-shape detection. ................................................................................................................. 69 Figure 4.2 The tree branch phantom and its reconstruction by DAS. (a) The tree branch phantom. The red arrows indicate the location and orientation of the line profile discussed in the following xvii paragraph. (b) Reconstruction result by DAS in the square-shape configuration with 128 transducers. .................................................................................................................................... 71 Figure 4.3 The tree branch phantom reconstruction results by VR-SGD algorithm. The first to forth rows refer to the image from the detection geometry by linear array (a-d), parallel-shape (e-h), L-shape (i-l) and Square-shape (m-p). From the first column to the fourth column, the image is reconstructed by 16, 32, 64 and 128 transducer elements correspondingly. The red lines at each row indicates the orientation of the transducers array in each configuration. ............................... 72 Figure 4.4 The line profiles of the reconstructed image. (a) The horizontal line profile indicated as line 1 in Figure 4.2(a). (b) The vertical line profile indicated as line 2 in Figure 4.2(a). .............. 73 Figure 4.5 The FORBILD head phantom and its reconstruction by DAS. (a) The FORBILD head phantom. The red arrows indicate the location and orientation of the line profile discussed in the following paragraph. Line 1 and 2 are horizontal line profiles and line 3 is a vertical line profile. (b) Reconstruction result by DAS in the square-shape configuration with 128 transducers. ........ 74 Figure 4.6 The FORBILD head phantom reconstruction results by VR-SGD algorithm. The first to forth rows refer to the image from the detection geometry by linear array (a-d), parallel-shape (e-h), L-shape (i-l) and Square-shape (m-p). From the first column to the fourth column, the image is reconstructed by 16, 32, 64 and 128 transducer elements correspondingly. The red lines at each row indicates the orientation of the transducers array in each configuration. ............................... 76 Figure 4.7 The line profiles of the reconstructed image. (a) The first horizontal line profile indicated as line 1 in Figure 4.5(a). (b) The second horizontal line profile indicated as line 2 in Figure 4.5(a). (c) The vertical line profile indicated as line 3 in Figure 4.5(a). ............................ 77 Figure 4.8 RMSE of VR-SGD based PAT reconstruction for the FORBILD head phantom. (a) RMSE for the reconstruction in different detection geometry with 128 transducers. (b) RMSE for xviii the reconstruction by the L-shape and the square-shape detection with different number of transducers. .................................................................................................................................... 78 Figure 4.9 PAT images of printed dots array detected by single linear array transducer. (a) Photograph of the phantom. (b) Initial image by DAS based on ð(ð, ð¡) . (c) Corresponding reconstruction result after VR-SGD for 25 times iteration. ........................................................... 81 Figure 4.10 Comparison of the cost function value versus the number of iterations under different tuning factors and step sizes. (a) Value of the misfit term. (b) Value of the TV term. (c) Value of the total cost function. There are 25 times of outer iterations. For each loop of outer iteration, 128 times of inner iterations is processed. ............................................................................................ 82 Figure 4.11 Comparison of the convergence rate for FGD, SGD and VR-SGD. (a) Normalized cost function versus the number of times that the gradient updates. (b) Normalized cost function versus the time each algorithm costs. For each algorithm, it updates the corresponding gradient 3200 times. .............................................................................................................................................. 83 Figure 4.12 PAT images of pencil leads array detected by a single linear array transducer. (a) Photo of the phantom. (b) Initial image reconstructed by DAS based on ð(ð, ð¡) . Figures (c) The reconstructed image after VR-SGD. .............................................................................................. 84 Figure 4.13 PAT images of a tube phantom detected by two linear array transducers. (a) Photograph of the phantom. (b) Image reconstructed by the first single linear array transducer which is placed on the top of the image. (c) Image reconstructed by the second single linear array transducer which is placed on the right-hand side of the image. (d) Combined image which overlaps the results from the first and second probe. (e) Image reconstructed by VR-SGD using the data from the two probes. ........................................................................................................ 87 xix List of Abbreviations ANSI American National Standards Institute CCD Charged-coupled device CMUT Capacitive micro-machined ultrasonic transducer CT Computed tomography DAQ Data-acquisition system DAS Delay-and-sum DCE-MRI Dynamic contrast-enhanced magnetic resonance imaging DCE-US Dynamic contrast-enhanced ultrasound DRE Digital rectal examination FGD Full gradient descent FWHM Full width of half maximum ICG Indocyanine green LED Light-emitting diode MPE Maximum permissible exposure mpMRI Multiparametric magnetic resonance imaging MRI Magnetic resonance imaging MVD Microvessel density Nd:YAG Neodymium-doped yttrium aluminum garnet OPO Optical parametric oscillator PA Photoacoustic PACT Photoacoustic computed tomography PAM Photoacoustic microscopy PAT Photoacoustic tomography xx PCa Prostate cancer PET Positron emission tomography PLD Pulsed laser diode PSA Prostate-specific antigen PSNR Peak to noise ratio RMSE Root mean square error SGD Stochastic gradient descent SNR Signal-to-noise ratio TR Time reversal TRUS Transrectal ultrasound TV Total variation TVA Total vascular area TVL Total vascular length UBP Universal back-projection US Ultrasound VR-SGD Variance-reduced stochastic gradient descent 1D One-dimensional 2D Two-dimensional 3D Three-dimensional xxi Acknowledgements Six years has passed since I came to Vancouver and joined Biophotonics group as a graduate student in the University of British Columbia. Throughout the study of PhD and writing of this dissertation I have received a great deal of support and assistance. I would first like to express my gratitude to my supervisor, Prof. Shuo Tang, whose expertise is invaluable in the formulating of the research topic and methodology in particular. Thanks for her patience guidance and kind support. Without her instructions and encouragement, I could not finish the project and become independent in research. Special thanks to Prof. Jong-in Youn at the Daegu Catholic University for this valuable help to develop the mimicking prostate phantom. I also thank Dr. Davood Karimi for his advice and discussion about the variance-reduced stochastic gradient descent method. I would like to thank my supervisory committee members Prof. Tim Salcudean, Prof. Robert Rohling and Prof. Purang Abolmaesumi. Thanks for their helpful advices and suggestions on my project and dissertation. I would like to thank the exam committee members and external reviewers. Thanks for their helpful advices and suggestions, I could successfully finalize my dissertation. I would also like to thank anonymous reviewers and all the excellent researchers in the world for having a discussion, question and suggestion. Their comments and replies inspire me a lot during my study and writing of each paper. In addition, I would like to thank all of the lab members and ex-members in the Biophotonics Laboratory at UBC for a pleasant time and great teamwork. The following listing order is based on the time sequence that I met with each of them. They are Jiali Yu, Leo Pan, Tom xxii Lai, Lin Huang, Weihang Shu, Mengzhe Shen, Myeong Jin Ju, Yuan Zhao, Mingyu Sheng, Xin Zhou, Jiayi Cheng, Yonghan Zhou, Ruochen Wen, Qihao Liu, Mukhlasur Rahman Tanvir, Christoph Brandt and Ning Zhang from UBC ECE. Also, I would like to thank Dr. Hamid Moradi and Corey Kelly who are in the photoacoustic imaging project, for sharing ideas and thoughts. I extend my thanks and gratitude to all my friends for genuine help and accompany. Especially, Prof. James B. Pawley and Prof. Christine Pawley recommends me to UBC and help me a lot to get used to living in Vancouver. Thanks Jiayu Luo for the company to overcome the miserable time in my life. Special thanks to my lovely girlfriend Elena Shi, who has always been my side. Finally, I dedicate this dissertation to my parents. Thanks for their endless support and consideration. They encourage me to keep pursuing my path as a researcher. Min Ai The University of British Columbia (Vancouver) October 2019 xxiii Dedication To everyone encounters in my life. To anyone who has ever loved. 1 Chapter 1: Introduction 1.1 Introduction of photoacoustic imaging Photoacoustic (PA) effect or optoacoustic effect is the generation of ultrasound wave in materials due to the absorption of modulated or pulsed light. It was discovered in 1880 by Alexander Graham Bell [1, 2]. In the early work, PA effect was used for spectroscopic analysis [3]. In the late 1990s, it was demonstrated to be potential to realize high resolution and deep penetration imaging by using pulsed laser for excitation and ultrasound transducer for detection [4]. Compared with conventional optical imaging modalities such as confocal microscopy [5], multiphoton microscopy [6] and optical coherence tomography [7], the advent of PA imaging extremely extends the penetration depth from a few millimeters to several centimeters [8, 9], because PA imaging detects ultrasound wave which is much less scattered than optical wave. In addition, PA imaging complements ultrasound imaging by distinguishing materials based on optical absorption which can provide biochemical specificity [10, 11]. Two PA imaging modalities have been developed based on the PA effect, the photoacoustic microscopy (PAM) and photoacoustic tomography (PAT). PAM applies raster scanning of either the pulsed laser beam or the ultrasonic transducer [8, 9]. As the acquired PA signal from each scanning point corresponds to a depth-resolved A-line information, PAM image can be formed directly. The lateral resolution of PAM is determined by the combination of the optical and acoustic foci. The optical diffusion limit is defined as the transport mean free path which is the mean distance after a photon's direction becomes random in tissue. Within the diffusion limit, optical resolution is achieved by tightly focusing the optical beam [12]. Beyond the diffusion limit, loose acoustic focusing determines the lateral resolution due to the weaker acoustic scattering in tissue [13, 14]. Axial resolution is related to the PA response and determined by the central 2 frequency of the detecting ultrasonic transducer. In general, higher central frequency leads to a higher axial resolution [9]. The maximum imaging depth is limited by acoustic attenuation and optical attenuation [13]. With a given optical wavelength, PA signal with lower acoustic frequency can penetrate deeper. Since there is a trade-off between the spatial resolution and the maximum imaging depth, PAM systems apply transducers with different central frequencies based on their specific need. Among the reported PAM systems so far, the lateral resolution varies from to , and the maximum imaging depth ranges from to correspondingly [15]. As a developing imaging modality, PAM involves in many biomedical and preclinical applications including vascular biology [16], neurology [16, 17], ophthalmology [18, 19] and etc.. For PAT or photoacoustic computed tomography (PACT), the beam is expanded and illuminates the whole volume of the sample instead of focusing on a finite region as described in PAM [8, 20]. The PA signals are collected by an array of unfocused transducers or by scanning a single-element transducer to function as an array. As signals are acquired simultaneously at multiple locations around the sample, PAT can achieve high frame rate without the need to perform point-to-point scanning as PAM does [8, 21]. However, an inverse algorithm is needed in order to reconstruct the tomographic image based on the spherical, cylindrical or planar detection geometry. The lateral resolution of PAT is mainly decided by the geometry of the transducer array. Similarly, the axial resolution and maximum imaging depth are both affected by the central frequency of the transducer. Since PAT usually applies a lower frequency than PAM does, deeper image penetration can be achieved which is desired for human clinical imaging [8, 20, 21]. So far, PAT has been studied to detect skin cancer [22], prostate cancer [23-25], breast cancer [26] and etc.. 220 nm560 μm 100 μm 4 cm3 1.2 Basic principles of PAT To generate PA effect efficiently, two conditions, namely thermal confinement and stress confinement, need to be satisfied. This is achieved by using a pulsed laser light in most PA imaging systems. Under the thermal confinement, heat diffusion can be neglected during the excitation. The laser pulse width, , should be shorter than the time scale of the heat dissipation of absorbed energy, , which can be approximately given as, , (1.1) where is the characteristic linear dimension of the absorbing structure, and is the thermal diffusivity of the sample [27]. Under stress confinement, high thermal elastic pressure can build up rapidly. The pulse width, , should be shorter than the time scale for the stress to transit the heated region, , which can be calculated by, , (1.2) where is the speed of sound [28]. For example, to achieve a spatial resolution of , with and as typical values for soft tissues [29], must be less than 30 ns. Hence, most of the current PA imaging systems apply pulses shorter than 15 ns. After meeting the confinement conditions to generate PA signal, an initial pressure can be generated caused by thermal elastic expansion, which can be expressed as [28, 30], , (1.3) where is the isobaric volume expansion coefficient in , is the specific heat in J/ð¾ â ðð, is the absorption coefficient in , is the light fluence in , is the local energy deposition density in , and is the Grüneisen coefficient given as . pttht2 / 4th p TL Dt »pL TDptst/s pL ct =c 50 μmpL =1.5 mm/μsc = 3 21.4 10 cm /TD s-= Â´ pt0p( )20 / p ap c C F Ab µ= =Gb 1K- pCaµ1cm- F 2J/cm aA Fµ=3J/cm G 2 / pc CbG =4 More details related to the propagation of excited PA signal is described in Chapter 4. Since optical absorption can be obtained in biological tissues from endogenous molecules such as hemoglobin, PA imaging is capable of imaging blood vessels and monitoring tumor angiogenesis. 1.3 Implementation of PAT The major components of a PAT system include a pulsed laser (e.g., a nanosecond pulsed laser), an ultrasonic transducer array or a single element transducer with scanning, a data-acquisition system (DAQ), and a computer for synchronization and image reconstruction [8, 20]. A nanosecond pulsed laser is usually applied in PAT imaging to satisfy the thermal confinement and stress confinement such as a neodymium-doped yttrium aluminum garnet (Nd:YAG) laser [31-33], pulsed laser diode (PLD) [34] or light-emitting diode (LED) [35]. Nd:YAG laser with an optical parametric oscillator (OPO) can provide high energy pulses with tunable wavelengths which is the most popular light excitation source applied in PAT experiment. However, the laser is expensive and bulky which is not suitable for a portable system. Many studies have investigated the feasibility of using PLD or LED in PAT which can make the system low-cost and portable. For means of light delivery, multimode fiber and optical fiber bundles have been utilized, which can provide more flexible and guided light delivery than free-space light delivery. More details about fiber delivery of pulsed laser will be introduced in Chapter 2. Compared with a single element transducer which needs mechanical scanning and many times of laser firing to form one frame of image, a transducer array takes the advantage of detecting PA signals in parallel and thus accelerates the imaging speed [20, 21]. Since linear and planar transducer arrays are commonly applied in the clinical ultrasound (US) systems, they are also employed in the US and PAT integrated systems. As the detection angle is limited for either type of those transducer arrays, many customized US transducer arrays have been developed for PAT 5 systems for increasing the data collection angle. In PAT for breast imaging, several types of transducer array have been demonstrated including semi-cylindrical [36], semi-spherical [37] and circular planar arrays [38]. Also, ring-array transducers are used in the PAT system for mouse brain imaging [39]. Except for the US transducers, other sensors are also developed for the detection of PA signal, such as the Fabry-Perot polymer film US sensor which utilizes optical interference to detect acoustic vibration [40]. The image reconstruction of PAT is based on the PA signals acquired by each transducer element. With the known detection geometry of the transducers, all the PA data is back projected to generate the image. Universal back-projection (UBP) [41] and time reversal (TR) [42] are reconstruction methods widely used. UBP has been demonstrated in spherical, cylindrical and planar detection geometries and the general formula can be written as [41], (1.4) where is the initial pressure corresponding to the optical absorption map, denotes the detected signal obtained by each transducer, is the solid angle of the detection surface with respect to the source point at , denotes the speed of sound, and represents the propagation time between the excitation of laser pulse and the receiving of the pressure wave by the US transducer. In UBP, is back-projected onto a spherical surface centered at with a radius of . The projected spherical surface becomes a circle in 2D imaging. In TR, the measured PA signal is propagated backwards into the image space in time reversed order. TR can work for any closed geometry. However, the wave equation needs to be solved by numerical methods such as k-space pseudo spectral method [42], which adds more computational time compared to UBP. Moreover, the detection geometry in clinical applications is usually not ideally ( ) ( ) ( )0 00,1 2 ,sSv tp tp p t t dt= -¶é Ã¹= - WÃª úW ¶Ã« Ã»òr r'rr' r( )0p r' ( ),p tr0Wr' sv t( ),p tr rsv t = -r r'6 enclosed or infinite. In that case, a limited-view detection can cause artifacts in the reconstructed image. In Chapter 4, more details related to PAT image reconstruction under limited detection will be discussed. 1.4 PAT for prostate imaging In the recent years, PAT/US integrated systems have been translated to many clinical studies for breast cancer [26, 37, 43], ovarian cancer [44, 45], prostate cancer (PCa) [23-25]. Prostate cancer is one of the most common cancers and leading cause of death among men [46]. The gold standard for PCa detection is the transrectal ultrasound (TRUS)-guided biopsy. However, the conventional US imaging is unable to provide reliable image of PCa, resulting in over-detection of indolent tumors or under-detection of malignant tumors [47]. To improve the accuracy of the diagnosis, many imaging technologies have been investigated including US elastography [48], magnetic resonance imaging (MRI) [49] and positron emission tomography (PET) [50]. Meanwhile, the development of PAT shows its compatibility to be integrated with the current clinical US system and performs non-invasive real-time imaging. Hemoglobin absorption enables PAT to provide extra optical contrast for blood vasculature, oxygen saturation and associated angiogenesis. By targeting cancer biomarkers such as tumor angiogenesis and hypoxia, PAT is capable to improve the detection of PCa by monitoring the vasculatures around tumors. Other than the potential of imaging PCa, PA spectral analysis has shown its feasibility for characterizing the dimension and concentration of the optical absorbers in tissues [51]. By quantifying the microscopic architecture of tumor tissue, PA spectral analysis has shown a correlation with Gleason score which is another prognostic factor for PCa [52]. Early study by Valeri et al. revealed the capability and potential of PAT to provide contrast between normal tissue and cancerous tissue in a prostatectomy specimen [53]. This result 7 encouraged many further studies of PAT for detecting and staging PCa, guiding biopsy and real-time monitoring of brachytherapy seeds. Several studies on prostate imaging in animal models have been carried out [23-25, 54, 55]. Bell et al. demonstrated in vivo PAT imaging on a canine prostate which was shown to be able to visualize brachytherapy seeds [55]. Moreover, Mohammad et al. showed the effectiveness for PAT to image blood-rich lesion in vivo in a canine prostate which was not detected by US [56]. However, human clinical study is limited until recently. Ishihara et al. firstly applied PAT to identify neurovascular bundle for human prostate in vivo during invasive radical prostatectomy [25]. Later, Ishihara et al. further showed the feasibility of imaging PCa angiogenesis based on the PA signals obtained by a TRUS probe [24]. Recently, Sri-Rajasekhar et al. demonstrated an integrated TRUS/PA system using intrinsic hemoglobin contrast as well as intravenous contrast agent indocyanine green (ICG) during PAT imaging of prostate [54]. PAT demonstrated a distinct PA contrast from the malignant tissue in one patient with advanced PCa, which was confirmed by the result from PET-MRI and biopsy. In a limited number of patients, consistent result was obtained based on the hemoglobin absorption within vasculature of the tumor, which was correlated to pathology results. ICG was shown to enhance the image contrast. A summary of the current PAT systems for prostate imaging is listed in Table 1.1 [23-25, 54, 55, 57]. Several PAT systems have been reported for prostate imaging in animal models but only two PAT systems have been reported for in vivo human prostate imaging [24, 25, 54]. To translate PAT for in vivo human prostate imaging, delivering high-energy laser pulses to the prostate effectively and safely is a challenge. In the reported human prostate imaging [24, 25, 54], transrectal illumination was employed in both systems. Fiber bundles were used to output more energy than using a single multimode fiber. While the rectal tract provides a large space to 8 integrate optical illumination and US detection together, high-level light energy is needed in order to cover the whole prostate when illuminated from the rectal wall. Currently, contrast agent such as ICG has been applied to enhance the contrast. Fluence rate can be further increased to safety allowed fluence to improve the PA signal intensity. On the other hand, transurethral illumination was demonstrated in prostate imaging on animal models [23, 55]. It is considered as an efficient design to achieve deep penetration in PAT prostate imaging. A more uniform illumination with less required energy can be achieved when light is illuminated from the urethral wall than from the rectal wall. As the urethral tract is much narrower than the rectal tract, the narrow tract limits the application of a fiber bundle. A multimode fiber which has a smaller diameter is more suitable for transurethral illumination. In the reported study, a multimode fiber illumination was used to localize brachytherapy seeds with limited imaging area, due to the lower energy that was delivered by the fiber [23]. How to deliver higher energy in a multimode fiber is a challenge. Furthermore, with high laser energy, how to maintain the laser safety limit at the laser-tissue contact surface is another challenge. For imaging prostate, the feasible detection geometry of the transducer array is restricted by the anatomy of prostate. Current approaches for prostate examination use transrectal ultrasound (TRUS). The transrectal view has been used for PAT imaging where TRUS probe is placed near the prostate. However, the detection by TRUS only allows a linear or curvilinear transducer array which cannot enclose the prostate. The limited view caused by the detection geometry poses an ill-conditioned reconstruction with missing features and strong artifacts. Those challenges can hinder the clinical application of PAT imaging of prostate in patients [23-25, 54, 55, 57]. Therefore, improving the design of light illumination and PA signal detection are important issues for in vivo prostate imaging. 9 Table 1.1 Current PAT systems for in vivo prostate imaging PAT Bell, et al. [23, 55] Ishihara, et al. [24, 25] Sri-Rajasekhar, et al. [54] Optical excitation Light source Output energy Fluence Wavelength Delivery Illumination geometry Nd:YAG 6 - 8 mJ/pulse 88 - 99 mJ/cm2 1064 nm multimode fiber Transurethral Alexandrite 60 mJ/pulse 11 mJ/cm2 756 nm Fiber bundle Transrectal Alexandrite - 10 mJ/cm2 800 nm Fiber bundle Transrectal Ultrasonic detection Transducer Detection geometry Center frequency Bandwidth Transrectal probe Linear / curvilinear - 5-9 MHz / 4-8MHz Transrectal probe Convex array 6.5 MHz - Transrectal probe Linear array 5 MHz - PA imaging Imaging depth Imaging object 5 - 30 mm Seeds in canine prostate Up to 20 mm PCa angiogenesis (3 patients) 30 - 40 mm Human prostate (20 patients) 1.5 Objectives and contributions PAT provides rich optical contrast beyond the diffusion limit with high-resolution imaging which is ideal for imaging thick tissues and clinical imaging. In PAT for prostate imaging, the imaging penetration depth is mainly limited by the optical attenuation in prostate tissue. Delivering higher energy pulses leads to deeper penetration depth and higher signal to noise ratio. Meanwhile, the design of a transrectal or transurethral illumination should be compact to fit the dimension of the rectal or urethral tract respectively, which restricts the size of the fiber or fiber bundle. Moreover, clinical application regulates the fluence to meet the safety requirement which conflicts with the demand for high energy. Thus, in vivo imaging of prostate cancer with PAT is currently limited by the lack of sufficient local fluence for deep tissue penetration and the risk of over irradiation near the laser-tissue contact surface. To translate PAT technology to in vivo prostate imaging, a PAT system with transurethral illumination is needed for prostate imaging. 10 The main objective of this dissertation is to investigate the potential of applying PAT to image the prostate. There are several challenges in developing PAT for prostate imaging including the design of light delivery scheme to achieve the coupling of high energy pulse, optimization of the illumination to meet the need of deep penetration and requirement of safe fluence rate, and the development of PAT image reconstruction which tolerates noisy data collected from the tissue. The main contributions I have made include: (1) developed a fiber coupling scheme for high energy laser pulses in PAT imaging; (2) designed a transurethral illumination probe for PAT imaging of prostate; (3) optimized PAT image reconstruction by developing a variance reduced stochastic gradient descent iterative algorithm. First, we proposed a fiber coupling scheme with a beam homogenizer for coupling high energy pulses into a multimode fiber. This scheme can benefit PAT applications that require miniaturized illumination or internal illumination with a small fiber. The beam homogenizer was achieved by using a cross cylindrical lens array, which provided a periodic spatial modulation on the phase of the input light. Thus, the lens array acted as a phase grating which diffracted the beam into a 2D diffraction pattern. Both theoretical analysis and experiments demonstrated that the focused beam was split into a 2D spot array that reduced the peak power on the fiber tip surface and thus enhanced the coupling performance. The theoretical analysis of the intensity distribution of the focused beam was carried out by Fourier optics. Our experiment achieved coupled energy at 48 mJ/pulse and 64 mJ/pulse in a 1000-μm and a 1500-μm-core-diameter fiber, respectively. The corresponding coupling efficiency was 70% and 90%. The high energy pulses delivered by the multimode fiber were further tested for PAT imaging in phantoms. The large improvement in coupling energy can benefit PAT with single fiber delivery to achieve large area imaging and deep penetration detection. 11 Second, we proposed the design of a transurethral illumination probe that addressed the limitations such as insufficient local fluence for deep tissue penetration and the risk of over-irradiation at the laser-tissue contact surface. A diffusing fiber was fabricated, which delivered light in radial illumination. The radial illumination was then reflected and reshaped by a parabolic cylindrical mirror to obtain nearly parallel side illumination with a doubled fluence. The fiber assembly was housed in a 25 Fr cystoscope sheath. A large laser-tissue contact surface area was obtained and the fluence on the tissue surface was kept below the maximum permissible exposure. By imaging a prostate mimicking phantom, a penetration depth of 3.5 cm at 10 mJ/cm2 fluence and 700 nm wavelength was demonstrated. Third, we demonstrated a variance reduced stochastic gradient descent (VR-SGD) algorithm for PAT image reconstruction through simulations and experiments. In the simulation, four detection geometries were studied, including linear array, parallel-shape, L-shape and square-shape arrays. The reconstructed image was improved after applying the iterative method and the corresponding root mean square error (RMSE) was reduced. Moreover, the L-shape detection was found to be comparable to the square-shape detection in terms of the image quality. The experiments using a single linear array transducer and two linear arrays positioned in the L-shape detection geometry were carried out. The VR-SGD algorithm was shown to improve the image contrast and reduce the noise and artifacts when compared to the result obtained by delay and sum (DAS). 1.6 Outline This dissertation is organized as follows: In Chapter 2, the fiber coupling scheme with a beam homogenizer is demonstrated for PAT imaging. A cross cylindrical lens array is applied to achieve beam homogenization. Its focusing 12 pattern is theoretically analyzed by Fourier optics and simulated in Zemax. The output energy is measured, and the corresponding coupling efficiency is investigated. Both theoretical analysis and experiments show the enhancement of the coupling performance. Ex vivo and in vivo imaging is demonstrated to show the PAT imaging with fiber delivery. In Chapter 3, a transurethral illumination probe is designed at the output end of the multimode fiber. Details of the design of the probe are explained, which mainly include fabricating the diffusing fiber, designing the parabolic cylindrical mirror and integrating the fiber with a cystoscope sheath. The illumination pattern of the probe is analyzed and characterized. Through this design, the output beam shows a large laser-tissue contact area and its fluence can be kept below the maximum permissible exposure. PAT imaging is demonstrated on a prostate mimicking phantom with the transurethral illumination probe. In Chapter 4, the VRSGD algorithm for PAT reconstruction is presented. Four types of detection geometries are discussed by using the VR-SGD algorithm including linear array, parallel-shape, L-shape and square-shape detections. Simulations of tree branch phantom and FORBILD head phantom are performed by VR-SGD algorithm and compared with the results obtained by DAS. PAT experiment by single linear array transducer and L-shape detection is demonstrated by applying the VR-SGD algorithm. Finally, Chapter 5 summarizes the work presented. The achievement and the limitation of the current work is discussed. Several suggestions are presented as long-term goal for future work. 13 Chapter 2: Design of high energy laser pulse delivery in a multimode fiber for photoacoustic tomography 2.1 Introduction Photoacoustic tomography (PAT) is a photoacoustic (PA) imaging modality characterized by large area illumination and multi-directional detection [8, 9]. Due to its relative deep penetration capability, PAT has been applied to study skin cancer [22], prostate cancer [23-25], breast cancer [26], abdominal organs [58], and brain activity [59]. PAT uses high energy laser pulses to excite PA signals. These pulses have been delivered traditionally in free space by prisms [22]. However, to translate PAT into clinical applications, delivering the high-energy pulses through an optical fiber is necessary for operational convenience and laser safety. Several prototype PAT systems using optical fiber for light delivery have been developed for imaging cancers [43, 60, 61]. Manohar et al. developed a breast imaging system employing a 588-element, 1-MHz ultrasound array with a 1064 nm pulse laser for excitation [43, 60, 61]. To illuminate the entire breast, a 1-to-9 splitting fiber bundle was applied to deliver nearly 81 mJ of pulse energy to the breast surface. A total 30 cases of breast malignancies from 31 patients were visualized with high contrast, showing PAT as a potential alternative to X-ray mammography for breast cancer imaging [26]. Horiguchi et al. developed a PAT system to monitor periprostatic tissue in prostate cancer patients undergoing radical prostatectomy [25]. Their light delivery was through two illumination apertures which were integrated with a 6.5 MHz transrectal ultrasonic transducers, delivering 60 mJ per pulse. By applying a 756 nm Q-switched pulse laser for illumination, periprostatic microvessels were imaged and located, improving the identification of the neurovascular bundle. Seven patients were imaged in this study. Furthermore, the small fiber diameter enables PAT to perform internal illumination to image internal organs, which can overcome the problem of insufficient light 14 penetration compared with external illumination [25, 45, 59]. For example, Lin et al. developed an illumination for oral cavity to provide sufficient light to image a mouse brain [59]. Using a multimode fiber with 2.8 mm core diameter, nearly 10 mJ energy at 780 nm was delivered, where the internal illumination enabled identification of internal brain structures. A major challenge in delivering high-energy laser pulses is the fiber damage of the tip surface caused by the high peak power density when focusing and coupling light into the fiber [62]. Thus the output energy from optical fiber delivery has been limited, which restricts the illumination area and penetration of PAT. Several measurements at different wavelengths and pulse widths have reported the damage threshold for fused silica fibers [63-65]. Gallais et al. observed a damage threshold at a peak power density of 9.3 GW/cm2 for silica glass, obtained at 1064 nm wavelength and 7 ns pulse duration after 1000 continuous pulses [65]. Robinson and Ilev observed a damage threshold at a peak power density of 3.7 GW/cm2 in a 100-μm-core-diameter fiber and 3.9 GW/cm2 in a 200-μm-core-diameter fiber, where the corresponding pulse energy was 2.9 mJ and 9.7 mJ respectively [66]. However, in the same work, the damage threshold was only 0.86 GW/cm2, corresponding to 13 mJ in pulse energy, in a 700-μm-core-diameter fiber. Since there is not a definite value of damage threshold for different fibers, commercial fibers usually have their own practical and theoretical damage threshold as a guideline for usage. For example, Thorlabs Inc. specifies a practical damage threshold of 1 GW/cm2 and theoretical damage threshold of 5 GW/cm2. Such high peak power occurs when focusing and coupling light into the fiber. A traditional coupling method uses a single lens to focus the collimated beam and couple the focused light into a multimode fiber with a large core diameter. However, for multimode fiber with relatively large core diameter (e.g. 1 mm), it is difficult to find a lens to generate a focal spot in the similar dimension. Therefore, the coupling of pulse laser usually places the fiber behind the focus to obtain 15 an enlarged beam which can match with the core size of the fiber [67, 68]. A fiber bundle can couple more energy than a single fiber due to its larger effective core diameter. Table 2.1 lists several systems that use a multimode fiber or a fiber bundle for pulse laser delivery in PAT. Variations in the coupled pulse energy depend on the different systems and applications. In a multimode fiber, the reported energy is typically limited to be less than 28 mJ. In a fiber bundle, energy up to 81 mJ has been reported [61, 68-70]. Although fiber bundle can couple higher power than multimode fiber, fiber bundle usually has a much larger size. Meanwhile, multimode fiber has relatively smaller size and it can benefit many applications that require miniaturized illumination or internal illumination. In PA system that uses a compact handheld probe to deliver light and collect signal [37], it is possible to reduce the probe dimension by using multimode fiber instead of fiber bundle. In applications which need internal illumination, such as transurethral illumination for prostate imaging [23] or PA endoscopy [58, 71], multimode fiber is necessary because small size is critical. However, the energy that can be coupled into a multimode fiber is still relatively low. In this chapter, the coupling scheme based on beam homogenization is applied to a PAT system with multimode fiber delivery. The goal is to improve upon 28 mJ/pulse in the multimode fiber, in order to develop the PAT system for imaging prostate through transurethral or transrectal illumination. A cross cylindrical lens array and a plano-convex spherical lens are used to homogenize and couple high-energy laser pulses into a multimode fiber. Compared with previous work [72, 73], a detailed theoretical analysis is carried out here on the focusing effect of the lens array based on the theory of Fourier optics. Optical simulations and experiments are also carried out to study the effect of the lens array. The theoretical result, optical simulation, and experimental results match with each other and show that the laser energy is split into multiple focal spots on 16 the fiber tip surface, lowering the risk of fiber damage. In the fiber coupling experiment with the lens array, a 1000-μm-core-diameter fiber can output 48 mJ/pulse energy and a 1500-μm fiber can output 60 mJ/pulse energy, which are much higher than other reported PAT systems with light delivery by a single fiber. The high-energy pulses delivered by the multimode fiber are further applied to PAT imaging on phantom and in vivo imaging. The improvement in fiber output energy enables PAT illumination and imaging of a larger tissue area. Table 2.1 Summary of PAT systems with fiber delivery Ref. Fiber Type Core Diameter Output Energy Coupling Efficiency Laser Wavelength Laser Pulse Width [22] Multimode 400 μm 25 mJ/pulse 50% 760 nm 60 ns [67] Multimode 1000 μm 28 mJ/pulse 85% 740-1100 nm 7 ns [68] Multimode 1500 μm 20 mJ/pulse 80% 700-900 nm 3-5 ns [61] Fiber bundle 9000 μm 81 mJ/pulse Not given 755 /1064 nm 60 ns/5 ns [69, 70] Fiber bundle 5000 μm 68 mJ/pulse 80% 700-900 nm 5 ns 2.2 Beam homogenization with a cross cylindrical lens array 2.2.1 Principle Figure 2.1 shows the schematic of the beam homogenization and fiber coupling system consisting of a cross cylindrical lens array and a plano-convex lens. The cross cylindrical lens array contains two 1D arrays of cylindrical lenses, an x array and a y array, joined together back-to-back, resulting in a 2D lens array effectively. This beam homogenization and fiber coupling system can split the incident laser beam to generate a more distributed intensity profile on the focal plane, and thus reduce the peak energy density on the fiber tip. 17 Figure 2.1 Schematic of the coupling optics. The incident beam from the pulse laser is split into beamlets by a cross cylindrical lens array, and the beamlets are focused by a plano-convex lens at its focal plane. To couple the light into a fiber, the fiber input end is placed at the focal plane of the plano-convex lens. The analysis of the intensity distribution on the focal plane is carried out based on the theory of Fourier optics. The optical model is shown in Figure 2.2. In Figure 2.2, defines the scalar electric field right after the lens array and presents the scalar electric field on the focal plane of the coupling lens. The distance between the lens array and the coupling lens is . The detection plane is the back focal plane of the coupling lens and is the focal length of the coupling lens. A plane wave with constant amplitude is assumed as the input in front of the lens array and can be expressed as, (2.1) where is the amplitude of the incident plane wave and denotes the transmittance function of the lens array. The transmittance function for a single cylindrical lens focusing in the x direction can be expressed as [74], , (2.2) where is the wave number of the incident plane wave and is the focal length of the cylindrical lens. By neglecting the thickness of the cylindrical lens array, its transmittance function can be expressed as, 0 ( , )U x h( , )fU x yd Lf0 ( , )U x h0 ( , ) ( , )clU A tx h x h= ÃA ( , )clt x h22( , ) clxjkff x y e-=k clf18 , (2.3) where is the number of elements of the lens array in both x and y, and denote the central position of the x and y cylindrical lens array, respectively, and is the pitch size of each cylindrical lens, which equals to the width of each lens element. The rectangular function is defined as, . (2.4) The rectangular function in Eq. (2.4) defines the aperture boundary for each cylindrical lens. The first summation term in Eq. (2.3) represents one side of the lens array, which only focuses light in the x direction. The exponential term indicates the phase change caused by each cylindrical lens. The second summation term presents the other side of the lens array, which only focuses light in the y direction. Because the focusing effect along the two directions is independent, the two summation terms can be multiplied directly to obtain the transmittance function of the lens array in 2D. The electric field at the focal plane can be obtained from the electric field using Fourier optics [75], , (2.5) where is the Fourier transform operator. The field at the focal plane is the Fourier transform of the field , multiplied with a phase term related to the distances d and . The 2 2( ) ( )2 2( , ) ( ) ( )n ncl cljk jkf fn ncln nt e rect e rectp px x h hx x h hx h- -- -é Ã¹ é Ã¹- -ê ú ê ú= ÃÃª ú ê úë Ã» Ã« Ã»å ån nx nhp0 if 0.5( ) 0.5 if 0.51 if 0.5xrect x xxì >ï= =íï <î( , )fU x y0 ( , )U x h[ ]2 2( ) 120( , ) ( , )LL Lk djkf j x yf ffLAeU x y e Uj fx hlæ ö+ -ç Ã·è ø= F( )F ( , )fU x y0 ( , )U x h Lf19 detailed derivation of Eq. (2.5) can be found in [75]. The intensity distribution can be obtained as the modulus squared of the field amplitude, . (2.6) Figure 2.2 (a) Optical model of the coupling scheme with the lens array. The lens array transforms the incident plane wave into multiple spherical wavelets. Constructive interference happens where the propagation angle ð½ð satisfies ððððð½ð = ðð , which forms the mth diffraction order on the focal plane. (b) 1D plot of the periodic phase modulation function introduced by the lens array. is the pitch of the lens array. As shown in Figure 2.2, the lens array acts similarly as a phase grating [75, 76]. Each element of the lens array introduces a phase delay according to its transmittance function, . (2.7) 2( , ) ( , )f fI x y U x y=p2 22( , ) ( ) ( )cljkft e rect rectp px h x hx h+-=20 The transmittance function can be viewed as a quadratic approximation to a spherical wave [75]. For an input of plane wave, it is transformed into an array of secondary spherical wavelets after the lens array. The lens array introduces a periodic spatial phase modulation that can be written as, , (2.8) where the comb function is a sequence of delta impulses uniformly distributed in space and defined as [75], . (2.9) A modulation depth is defined as below to describe the strength of the phase modulation, . (2.10) When the input light is coherent, the multiple spherical wavelets lead to multi-beam interference. At certain propagation angles , constructive interference happens, where the path length difference between two adjacent beams is an integer number of wavelengths. The constructive interference condition is similar as the diffraction orders of phase grating. The coupling lens focuses the beams at the same angle to one spot on its focal plane, forming the mth diffraction order. The position of the mth order diffraction spot can be written as: . (2.11) For the effective 2D lens array, the light intensity distribution on the focal plane of the coupling lens forms a 2D diffraction pattern. The period of the diffraction pattern is . (2.12) ( ) ( )2 2, ( ) ( )2 clkrect rect comb combf p p p px h x h x hf x h+ æ ö æ ö é Ã¹= - Ã Ã Ã ÃÃ§ Ã· ç Ã· ê úè ø è ø Ã« Ã»( )( )ncomb x x nd= -åfD( ) ( )2 2 2/ 2 / 22 4cl clk p p kpf ffé Ã¹+Ã« Ã»D = =mqmqmLtan sin Lm L m L mm ff fplq q Ã ÃL = » =Lfpl ÃDL =21 A phase grating can be much more efficient in diffracting light than an amplitude grating. The diffraction efficiency for each order depends on the modulation depth and the shape of the periodic phase modulation function. For small modulation depth, the zero-order beam contains most of the power. As the modulation depth increases, more and more power can be diffracted into the higher order beams. 2.2.2 Theoretical and experimental validation of the focal intensity distribution Both theoretical analysis and experimental measurement have been carried out to study the intensity distribution on the focal plane. The cross cylindrical lens array is obtained from SUSS MicroOptics (Nr.18-00142, Hauterive, Switzerland). Its specifications are listed in Table 2. The focal length of the cylindrical lens array is calculated to be 133 mm. The coupling lens is a plano-convex spherical lens with 75 mm focal length (LA1608, Thorlabs, Newton, United States). The input beam is modeled as a Gaussian beam with a beam waist of 5 mm and a wavelength of 700 nm. The theoretical analysis of the intensity distribution on the focal plane can be calculated based on Eqs. (2.5) and (2.6). The transmittance function of the cross cylindrical lens array can be obtained based on the parameters in Table 2.2, where , and . The calculated intensity distribution map on the focal plane is shown in Figure 2.3(a) and a corresponding intensity line profile from the center position is shown in Figure 2.3(b). A periodic spot pattern can be clearly seen, which covers an area larger than 0.5 mm Ã 0.5 mm. The pattern contains more than a nine by nine spot array, caused by interference. The calculated period based on Eq. (2.12) is 52.5 µm which matches with the period in Figure 2.3(a) and (b). The uneven light intensity distribution among the spots is a property of the phase grating, which depends on the modulation depth and shape of the phase grating [75, 76]. fD9n = 1015 p mµ= 133 clf mm=22 Table 2.2 Specifications of the cross cylindrical lens array Parameter Name Value Number of elements 9Ã9 Pitch 1015 µm Radius of curvature 61 mm ± 5% Focal length 133 mm Thickness 1.2 mm Size 10 mm Ã10 mm ± 0.05 mm Material Fused silica Zemax simulation can also be carried out to simulate the optical system. The Zemax file of the cross cylindrical lens array is provided by SUSS MicroOptics. The light distribution on the focal plane can be simulated by the physical optics analysis module. The 2D intensity distribution map is shown in Figure 2.3(c) and a corresponding line profile is shown in Figure 2.3(d). The Zemax simulation result matches very well with the above analytical result. The number of spots and the spacing between spots are all matched. A small difference is that the spot size obtained by Zemax simulation is larger than the ideal focal spot depicted in Figure 2.3(a). This may have been caused by the finite value of the thickness of the cross cylindrical lens array used in the Zemax simulation, where the lens thickness is neglected in the theoretical calculation in Figure 2.3(a). 23 Figure 2.3 The light intensity distribution map (left column) and a corresponding intensity line profile (right column) at the focal plane. (a) and (b) are obtained by theoretical analysis; (c) and (d) are obtained by Zemax simulation; (e) and (f) are obtained by experiment. The location for the line profile is marked by an arrow. In the experiment, a Fabry-Perot laser operated at 660 nm wavelength (S1FC660, Thorlabs, Newton, United States) is used instead of a nanosecond pulse laser only for the characterization of the beam profile, concerning the potential damage. The collimated beam from the fiber coupled 24 laser is expanded to nearly 9 mm in diameter to match with the output beam size of the nanosecond pulse laser. A charged-coupled device (CCD) (MU300, AmScope, Irvine, United States) monitors the light intensity at the focal plane. The 2D intensity distribution map is shown in Figure 2.3(e) and a line profile is shown in Figure 2.3(f). As we can see, the experimental results match reasonably well with the analytical and Zemax simulation results. With the cross cylindrical lens array, the light distribution on the focal plane is split into a spot array. Compared with the single focal spot generated by a single lens without the lens array, such a spot array distribution can lower the peak power density on the fiber tip and thus enable coupling higher energy into the fiber without damage. For some high energy pulsed lasers used in PAT, âhot spotsâ exist in the beam profile due to beam inhomogeneity, where the power density is higher than the surrounding area. The lens array based beam homogenizer can also reduce the effect of âhot spotâ. If a âhot spotâ exists in the input wave, it will be transformed to spherical wavelet after the lens array. The spherical wavelet contains different propagation angles, which will be focused to different locations on the focal plane of the coupling lens. Therefore, the energy from the âhot spotâ will be redistributed to a relatively large area on the focal plane, reducing the potential of thermal damage. To check whether the beam distribution pattern at the input of the fiber affects the fiber output, the beam profile at the output of the fiber is measured experimentally. The diverged beam output from a 1000-μm-core-diameter fiber (FT1000UMT, Thorlabs, Newton, United States) is imaged by the CCD camera and shown in Figure 2.4. The 2D intensity distribution map is shown in Figure 2.4(a) and an intensity line profile is shown in Figure 2.4(b). A quasi-Gaussian beam is obtained at the fiber output. This shows that the spot pattern at the input end of the fiber can be eliminated after propagating inside the multimode fiber. The elimination of the spot pattern is 25 likely caused by the effect of a large number of spatial modes propagating in the multimode fiber. The speckles observed in the beam profile are also caused by the nature of multimode fiber [32]. The quasi-Gaussian beam profile at the output of the multimode fiber shows that using the cross cylindrical lens array for light coupling does not affect the uniformity of illumination on sample for PA imaging. In fiber bundle, the beam profile is usually discretized because the bundle contains many fibers. Therefore, a multimode fiber can provide a more well-defined beam profile than fiber bundle. Figure 2.4 Measured light intensity distribution map (a) and line profile; (b) at the output of the multimode fiber. The location for the line profile is marked by an arrow. 2.3 Coupling high energy pulses into a multimode fiber for PAT 2.3.1 Experimental setup The cylindrical lens array is used to couple high-energy pulses into a multimode fiber. A PAT system using the multimode fiber to deliver light is developed. A schematic of the PAT system is shown in Figure 2.5. A Q-switched Nd:YAG laser (Surelite OPO Plus SLIII-10, Continuum, San Jose, United States) generates 3-5 ns laser pulses at 10 Hz repetition rate. The laser output is frequency doubled to 532 nm, which then pumps an optical parametric oscillator (OPO). The OPO output wavelength is tunable from 675 nm to 2500 nm and the maximum output energy is 120 mJ. 26 After the OPO, the light is coupled into a multimode fiber (FT1000UMT/FT1500UMT, Thorlabs, Newton, United States) by the cylindrical lens array and the coupling lens. The multimode fiber delivers laser pulses to the sample. The beam diameter after OPO is 9.5 mm. Hence, the lens array can be fully covered. The PA signal is detected by a wide-bandwidth linear transducer array (L14-5/38, Analogic, Richmond, Canada) with 128 channels, a central frequency of 7.2 MHz, and a 70% fractional bandwidth at -6 dB. The received data is then transferred to a data acquisition system (DAQ) and is sent to the ultrasound machine (Ultrasonix MDP, Analogic, Richmond, Canada). Synchronization between the laser firing and data acquisition is achieved by using a signal from the laser to trigger the DAQ. The reconstruction method uses the basic back-projection algorithm [41]. Figure 2.5 Experimental setup of the photoacoustic tomography with fiber delivering high energy pulses. The laser beam from the optical parametric oscillator (OPO) is homogenized and coupled into a multimode fiber to deliver light to phantom. DAQ: data acquisition system. 2.3.2 Performance of fiber coupling To evaluate the coupling performance with the lens array, the pulse energy at the fiber output is measured by an energy sensor (J-50MB-YAG, Coherent, Palo Alto, United States). The output pulse energy is averaged over one hundred consecutive pulses. The laser wavelength is tuned over the spectrum range from 675 nm to 900 nm. Light coupling into multimode fibers with 1000 µm 27 and 1500 µm core diameters are investigated respectively. The fiber length is about 1 to 2 m. By gradually increasing the laser power, the maximum stable output without fiber damage is recorded at each wavelength. The results of the output pulse energy at different wavelengths are shown in Figure 2.6(a) and (b) for the 1000 µm and 1500 µm fiber, respectively. The maximum output for the 1000-μm-core-diameter fiber is about 48 mJ/pulse at 750 nm wavelength. By comparison, the maximum output for the 1500-μm-core-diameter fiber is about 64 mJ/pulse at 700 nm wavelength. In Figure 2.6(a) and (b), the fluctuation of the output energy of the fiber is mainly affected by the energy variation of the laser output of the pulse laser. Figure 2.6(c) shows the coupling efficiency of the two types of fiber. The coupling efficiency is calculated as the ratio of the energy at the output end of the fiber to the energy at the input end of the fiber. The coupling efficiency reaches 70% and 90% for the 1000-μm and 1500-μm-core-diameter fibers, respectively. The lens array and the coupling lens introduce some additional insertion loss, which is measured to be ~10%. The laser output from the OPO is slightly divergent, which affects the intensity distribution at the focal plane and makes the dimension of the focusing pattern slightly larger than the theoretical size estimated for a collimated beam. As a result, the 1500-μm-core-diameter fiber under the same optical setup obtains higher coupling efficiency than the 1000-μm-core-diameter fiber. However, 1500-μm-core-diameter is more rigid than the 1000-μm-core-diameter fiber and fiber with smaller core diameter is more preferable for transurethral illumination in prostate imaging. Hence, the 1000-μm-core-diameter fiber is selected and applied in the following PAT imaging. In comparison, the conventional coupling scheme by using a single lens with 75 mm focal length (LA1608, Thorlabs, Newton, United States) is also tested. To avoid high power density on the optical focal plane, the fiber is placed around 3 mm behind the focal plane to have an equivalent 28 size of spot on its tip surface. The coupling energy and efficiency are found to be much lower than the situation with the lens array. The coupled energy with the single lens scheme is less than 20 mJ/pulse and the coupling efficiency is typically less than 40%. Moreover, the focused light could cause air breakdown at the focus which can indirectly damage the fiber in this situation. Another single lens with longer focal length and a two-lens system for reducing the beam size have also been tried but the improvement is limited. Furthermore, those coupling schemes rely on a good beam quality. However, neither the intensity distribution nor the beam divergence follows the theoretical prediction due to the poor beam quality output from the OPO. Figure 2.6 Output power characterization at different wavelengths and fiber core sizes. (a) Averaged pulse energy coupled into a 1000-µm-core-diameter fiber (blue curve). The error bar indicates the energy fluctuation measured over one hundred consecutive pulses. The laser energy measured at the input of the fiber is shown in red curve. (b) Similar power characterization as in (a) but with the light coupled into a 1500-µm-core-diameter fiber. (c) Comparison of the coupling efficiency of the two types of fiber. 29 In the section 2.2.1, we have shown that the beam homogenization system with the cross cylindrical lens array generates a 2D array of beam spots on the focal plane of the coupling lens. The distributed intensity on the focal plane reduces the peak energy density on the fiber tip. Thus higher energy pulses can be coupled into the multimode fiber without damage. A practical damage threshold is usually found to be 1 GW/cm2, which is equivalent to 40 mJ/pulse with 5 ns pulsewidth for a 1000-μm-core-diameter fiber. Our coupling result exceeds the practical limit but is still far below the theoretical damage threshold, which is 5 GW/cm2 and corresponds to 200 mJ/pulse for a 1000-μm-core-diameter fiber. Further enhancement of the coupling energy relies on improving the beam uniformity on the fiber tip surface and preventing air breakdown. For example, a top-hat beam homogenizer can be used to transform the Gaussian incident laser beam into a near uniform intensity distribution. As the OPO outputs an asymmetric beam shape with significant divergence and high pulse energy, the beam shaping performance of the top-hat beam homogenizer needs more investigation by systematic testing. To avoid air breakdown, the fiber tip can be protected by putting it in a chamber filled with noble gas or non-combustible liquid which will also result in a higher coupling energy. 2.3.3 Preliminary PAT Imaging The high-energy pulses delivered by the multimode fiber are tested for PAT imaging. The 1000-μm-core-diameter fiber is used and the maximum output is about 48 mJ/pulse. A phantom containing a piece of paper with printed dots array placed in gelatin is imaged. The diameter of the dots is 0.1 mm and the spacing between two adjacent dots is 1 mm. The sample is illuminated by the fiber placed 4 cm away, which forms a circular beam around 10 mm in diameter. The PAT image of the dots array is shown in Figure 2.7(a). The excited area is determined by the beam size. Next, a 5 mm thick chicken breast tissue is added on top of the phantom to increase the optical 30 depth and simulate optical diffusion in biological tissue. Figure 2.7(b) shows the PAT image of the dots array with the 5 mm thick chicken breast tissue on top. As the light is diffused inside the chicken breast tissue, the beam size is increased significantly and PA signal from a larger area is excited. The reconstructed image in Figure 2.7(b) shows at least 7 cm2 excitation area. Good signal level is still obtained over this large excitation area, as the signal-to-noise ratio drops from 16.1 dB to 10.5 dB due to the decrease of fluence. In vivo PAT imaging is also performed on the human forearm of a volunteer. The excitation wavelength is set to be 750 nm and the illumination area on the skin is a circle of ~3 cm diameter. The energy density at the skin surface is controlled to be under 10 mJ/cm2, which is lower than the maximum permissible exposure (26 mJ/cm2 at 750 nm) by ANSI standard [77]. A coupling pad is placed between the forearm and the linear transducer array for acoustic coupling. The illumination fiber is placed on the same side with the transducer, forming a diverged circular pattern on the region of interest. Figure 2.8(a) shows the reconstructed PAT image (colored) on top of the ultrasound image (gray scale) and its signal-to-noise ratio (SNR) for PA signal is 15.8 dB. PAT signal is obtained from the skin surface and blood vessels. Three veins can be identified from the PAT image. One shallow capillary with lower contrast appears between the second and the third vein, which is also labeled in Figure 2.8(a). Figure 2.8(b) shows the photography of the imaging location in the forearm and the corresponding blood vessels. These experimental results demonstrate the large-area imaging capability of the PAT system with the high energy pulses delivered by the multimode fiber. 31 Figure 2.7 PAT image of printed dots array. (a) Imaging result without the cover of chicken breast and the effective imaged area in a 10-mm-diameter circle. (b) Imaging result with the cover of 5 mm thick chicken breast and the increased effective imaged area in a 30-mm-diameter circle. Figure 2.8 In vivo PAT image of forearm. (a) Overlay of PAT image (colored) on top of ultrasound image (grey scale). The PAT image highlights the skin surface and locations of several blood vessels labeled with numbers. (b) Photograph of the forearm where the blood vessels imaged by PAT are also marked. The dash line shows the place of the cross-sectional image. 2.4 Summary In this chapter, a fiber coupling scheme with beam homogenization is developed and demonstrated for PAT imaging. It increases the pulse energy that can be coupled into a multimode fiber and reduces the fiber damage. The beam homogenization is achieved by using a cross cylindrical lens 32 array, where the focused pulse energy is distributed over a 2D array of beam spots. Our results demonstrate that by adding the beam homogenization, the pulse energy delivered in a single fiber can be significantly improved. We have obtained the delivery of 48 mJ/pulse by a 1000-μm-core-diameter fiber and 64 mJ/pulse by a 1500-μm-core-diameter fiber. To the best of our knowledge, this is the highest output energy reported to date using a single fiber delivery for PAT imaging. With the increased output energy, the PAT system can illuminate a larger area at the same fluence or provide deeper penetration. 33 Chapter 3: Photoacoustic tomography for imaging prostate: a transurethral illumination probe design and application 3.1 Introduction Prostate cancer (PCa) is the most common cancer and the second leading cause of cancer deaths in North America men [78]. Digital rectal examination (DRE) and serum prostate-specific antigen (PSA) are the current techniques for early detection of PCa. Upon a suspicious DRE and/or a high PSA reading, a transrectal ultrasound (TRUS) guided systematic random biopsy is performed [79]. If there are PCa-positive biopsy cores, a PCa diagnosis is confirmed. However, DRE is limited by low overall sensitivity (37%) [80] and PSA is limited by low specificity (36%) [81]. Because TRUS cannot identify PCa, the systematic sextant biopsy approach leads to repeated biopsies [82, 83]. Meanwhile, various magnetic resonance imaging (MRI) techniques have been investigated for detecting PCa and show a consistent higher accuracy than conventional ultrasound (US) [84]. Combining different modalities to perform multi-parametric MRI (mpMRI) has the advantage of further increasing the accuracy. When mpMRI is available, it can be fused with US in order to guide some of the biopsy sampling to areas that appear suspicious on mpMRI [85]. However, unavoidable movement of the prostate affects the US-MRI registration and therefore the biopsy localization error [86, 87]. Regardless of the techniques, since not all prostate tissue is sampled, there could be underestimation of the Gleason grade which is related to PCa aggressiveness used in PCa patient management [83]. The limitation of the current diagnostic tools reflects a need for a better detection method which can realize targeted biopsy or even PCa diagnosis through in vivo imaging. Meanwhile, the key role of angiogenesis in cancer growth and development in prostate has been revealed and higher microvessel density (MVD) has been reported as a potential prognostic indicator in PCa 34 [88]. This makes tumor angiogenesis a prospective objective of the detection for PCa, which leads to a wide range of new techniques in US [89] and MRI [86] for angiogenesis imaging. For US, both Doppler US and dynamic contrast-enhanced US (DCE-US) have been reported to be capable of monitoring the increase of MVD [88]. However, Doppler US only has high sensitivity to lesions with high Gleason score because it cannot detect microvessels but only large tumor feeding vessels [86]. The procedure of imaging prostate by DCE-US is time-consuming which limits its application for targeted biopsy [89]. Dynamic contrast enhanced MRI (DCE-MRI) has also been reported to image microvessel vascularity in prostate tissue [90]. As prostate is already highly vascularized, DCE-MRI has difficulty in discriminating cancer from prostatitis in the peripheral zone and from the high vascularized benign prostate hyperplasia [91]. Photoacoustic tomography (PAT) is an emerging imaging technique with potential for detecting PCa [9]. The photoacoustic (PA) signal in biological tissue originates from optical absorption in endogenous molecules such as hemoglobin. Upon illumination by a short laser pulse, an acoustic pressure wave is generated due to laser heating and thermal-elastic expansion in tissue. The acoustic wave then propagates towards the tissue surface over the depth determined by the acoustic attenuation inside the tissue. Assuming an attenuation coefficient in soft tissue of 0.5 dB cm-1, the penetration depth for 5 - 10 MHz ultrasonic wave ranges from 4 cm to 2 cm correspondingly [92]. The detection of the PA signal can be achieved by using a commercial US system which allows dual-modality imaging including both optical and ultrasonic contrasts. Previous studies have been carried out to use PAT for PCa detection and localization of brachytherapy seeds [24, 25, 53, 55, 93-96]. In 2003, Andreev et al. reported PAT imaging of human PCa ex vivo, which showed that the contrast difference between normal prostatic tissue and cancerous tissue was sufficient to locate tumors [53]. In another study, Wang et al. demonstrated 35 non-invasive PAT imaging on canine prostate in vivo [93]. A lesion at 10 mm depth was imaged with high sensitivity and contrast-to-noise ratio. In 2016, Ishihara et al. published the work of PAT imaging on human prostate in vivo [25]. The system was designed to detect microvascular structures such as the neurovascular bundle and cancer angiogenesis for assisting the procedure of radical prostatectomy or targeted biopsy. A TRUS probe (6.5 MHz center frequency, 128 elements linear array) was used to receive both PA and US signals. The illumination part delivered 60 mJ/pulse energy at 765 nm wavelength from two convex apertures integrated with the TRUS probe [24]. The estimated fluence on the tissue surface was 11 mJ/cm2. The PA signal was obtained from a depth of 17 mm. However, the strong signal near the rectal wall caused by the high surface fluence hindered the discrimination of the signal from the neurovascular bundle in the same region. With the same system, Ishihara et al. further demonstrated the feasibility of imaging PCa angiogenesis by correlating the PA signal intensity to the total vascular area and total vascular length for prostate cancer [24]. Although the results need to be further confirmed due to the small number of patients, this pilot study shows that PAT has great potential to image PCa. The clinical translation of PAT to prostate intervention still has many challenges. The amplitude of PA signal is proportional to the local fluence which determines the penetration depth and the signal-to-noise ratio (SNR). Therefore, high irradiation energy is always desired in PAT imaging. In the reported systems [55, 94], the integration of light delivery with the TRUS probe has the drawback of light attenuation by the rectal wall. Compared to transrectal illumination, theoretical analysis has shown that transurethral illumination could improve the light delivery [97, 98]. It is more efficient for light to reach the entire volume of the prostate from its center through the urethral tube, than from the rectum. Transurethral illumination was demonstrated in experiment by Bell et al. to image brachytherapy seeds and tested in canine prostate in vivo [23]. Using a 1000-36 µm-core diameter fiber, up to 8 mJ/pulse energy at 1064 nm wavelength was delivered. The fiber tip was polished at 45° angle and together with an end mirror at 45° angle light was directed to side illumination. The fiber and end mirror were housed in a 2 mm diameter quartz tube. The beam size on the tube surface was calculated to be 1.7 mm2 [23]. Since light was emitted from the fiber tip and then reflected to side illumination, the beam size was very small, which resulted in a high fluence at the tube surface. When the pulse energy was 7.5 mJ, the reported fluence was 414 mJ/cm2 at the tube surface. To reduce the fluence, the authors proposed to increase the tube diameter. With a 5.3 mm tube diameter, the authors estimated the fluence to be below 99 mJ/cm2 [23], which was just below the safety limit of 100 mJ/cm2 at 1064 nm wavelength for skin. For translating PAT to clinical applications, the laser safety limit needs to be addressed. In fact, the value of laser safety limit depends on the wavelength and radiation duration. Thermal effect that can cause tissue damage is a main concern when using visible and near infrared light. For a single pulse with pulse duration less than 100 ns, a maximum permissible exposure (MPE) of fluence, , is defined as [77], , (3.1) where is the wavelength. The MPE value decreases exponentially as the wavelength decreases and the value is as low as 20 mJ/cm2 at 700 nm. When the accumulated exposure time by multiple pulses is longer than 10 s, a maximum permissible exposure in terms of average power, , is defined as [77], . (3.2) fluenceMPE( ) ( )( )2 0.7 2fluence 220 10 mJ/cm 0.7 1.05MPE100 mJ/cm 1.05 1.4m mm ml µ l µµ l µ-ì Â´ £ £ï= í£ £ïîlfluenceMPEpowerMPE( ) ( )( )2 0.7 2power 20.2 10 W/cm 0.7 1.05MPE1 W/cm 1.05 1.4m mm ml µ l µµ l µ-ì Â´ £ £ï= í£ £ïî37 Here the average power is the product of the repetition rate and the energy per pulse. Both the requirements on and defined in Eq. (3.1) and (3.2) should be satisfied in order to meet the safety standard. The demand of high energy output for high SNR and the laser safety limit need to be considered together when designing an optimal illumination in PAT. In this chapter, a transurethral illumination probe is reported for PAT imaging of PCa. The probe contains a side-illumination diffusing fiber, and a cylindrical mirror to reflect and shape the beam profile to one side. A method is proposed to fabricate the diffusing fiber to achieve side firing over a few centimetersâ length of the fiber. The 3D printed silver mirror with cylindrical parabolic surface is integrated with the fiber to increase the fluence on one side for deeper penetration. Meanwhile, the fluence on the laser-tissue contact surface is controlled to be under the MPE value. The probe is tested on a tissue phantom which mimics the optical properties of human prostate. A 3 cm optical penetration depth is achieved which meets the need to image the entire prostate. 3.2 Design of the transurethral illumination probe In Chapter 2, a scheme of delivering high energy laser pulses through a multimode fiber is described [32]. A cross cylindrical lens array is placed before the input end of the fiber to homogenize the incident beam which leads to high energy coupling and delivery, e.g. 50 mJ/pulse in a 1000-µm-core diameter multimode fiber. The goal of the illumination probe design is to fully utilize the high output energy while keeping the fluence on tissue to be within the safety limit for both single pulse and long duration exposure. 3.2.1 Configuration of the probe The transurethral illumination probe design is shown in Figure 3.1(a). A 1000-µm-core diameter multimode fiber (Thorlabs, Newton, New Jersey, United States) is fabricated to diffuse light from powerMPE fluenceMPE38 its side surface. The fiber is held by a 3D-printed plastic mount. To increase the fluence and realize side-firing, light exiting from the diffusing fiber is directed to one side of the probe by a cylindrical reflection mirror constructed by 3D-printing. The mirror is glued with the fiber mount where the fiber is aligned to the focus of the cylindrical mirror. A 45° reflection surface at the end of the cylindrical mirror reflects and re-uses the remaining light exiting from the forward direction of the fiber tip to the side direction also. The assembled fiber is then inserted into a 25 Fr cystoscope sheath. The cystoscope sheath can support and protect the optical components when the probe is inserted in the urethra. Figure 3.1(b) shows a photograph of the transurethral illumination probe. Visible light at 660-nm wavelength is coupled into the fiber to display the beam profile. Figure 3.1(c) and 1(d) show the longitudinal and axial views, respectively, of the probe design. The coordinate axes are depicted in Figure 3.1(d), where y is along the long direction of the probe, x is in the perpendicular direction of the probe, and z is the direction of light propagation in the side illumination. The origin of the coordinate system is defined at the center of the fiber. Two different shapes of the cylindrical mirror have been investigated, a circular cylindrical mirror and a parabolic cylindrical mirror. Using the defined coordinates, the surface of the circular cylindrical mirror can be expressed as, , (3.3) where is the radius of the circular cylinder. Under paraxial approximation, the focus of the circular cylindrical mirror is a line parallel to y-axis and its focal distance is . Similarly, the surface of the parabolic cylindrical mirror can be written as, , (3.4) where its focal distance is . 2 22rz x ræ ö- + =ç Ã·è ør/ 2circularf r=22 2x pzp= -/ 2 parabolicf p=39 Figure 3.1 Design of the transurethral illumination probe. (a) Overview of the transurethral illumination probe. (b) The beam profile from the probe displayed by a 660 nm light. (c) Longitudinal cross-section view of the probe. (d) Axial cross-section view of the probe. The rigid surface of the cystoscope sheath is used to support the urethra wall and avoid direct contact between the urethra wall and the optical fiber. The distance between the top surface of the cystoscope and the fiber is considered as the minimal laser-tissue contact distance, e.g. a 25 Fr cystoscope sheath can provide ~5 mm laser-tissue contact distance. In the following sections, the design and analysis of the probe will be described in detail. 3.2.2 Fabrication of the diffusing fiber The normal size of a healthy prostate is approximately 25 ml in volume, or cm in the anterioposterior, transverse, and cephalocaudal prostate diameter, respectively [99]. With advancing age, the prostate size increases and varies in a wide range [100, 101]. For imaging the prostate, the illumination pattern of the transurethral illumination probe should match with the 3 5 3´ Â´40 cephalocaudal dimension along the urethra. Thus, a diffusing fiber that can provide side illumination over a few centimetersâ length is required. Several approaches have been reported to process a diffusing fiber. Mechanical processing methods [57, 102], such as laser micro-machining and abrasion by sandpaper, have a risk to break the denuded fiber during the process, especially for large core diameter fiber, e.g. 1000 µm. Also, manual abrasion by sandpaper requires experience to achieve the uniform roughness of the diffusing fiber. Meanwhile, Kosoglu et al. demonstrated an etching process that obtained a micro-needle of 33 µm diameter and 3 mm length by immersing a 105-µm-core diameter multimode fiber in 50% hydrofluoric acid for 50 minutes [103]. However, the replication of this approach on 1000-µm-core diameter multimode fiber is problematic. It takes much longer to etch the fiber due to the larger core size. Moreover, the etched fiber becomes very fragile since the core size is reduced significantly. To improve the etching process suitable for large core fiber, we have developed a process using glass etching cream (Armour etch, Armour Product, New Jersey, United States) to process large-core-diameter fiber and create a long section of diffusing fiber. Figure 3.2(a) illustrates the procedure in three steps. First, the buffer of the multimode fiber is removed by a stripper. The entire circumference of the cladding is scratched and removed by a sharp knife because the cladding layer does not react with the etching cream. In step 2, the entire bare fiber is etched by etching cream for nearly 4 hours. Due to the slow rate of reaction, the diameter of the bare fiber changes only slightly even after the long etching duration. At last, the cream is washed out by water and an irregular rough surface is formed as shown in Figure 3.2(b). The etched fiber diffuses light circumferentially in 360° angles. Compared with the etching process by hydrofluoric acid, our process is moderate which only creates roughness on the side surface of the fiber without 41 greatly reducing the core size. As a result, the etched part does not break easily. By coupling visible light at 660-nm into the fiber, the diffusing fiber shows a radial light illumination pattern along the etched surface as shown in Figure 3.2(c). Figure 3.2(d) shows the diffusing fiber and its illumination pattern on a piece of paper which has a relatively uniform light distribution. Figure 3.2 Schematic illustration of the etching process (left) and photographs of the diffusing fiber (right). (a) Three steps of the etching process. (b) Rough surface of the diffusing fiber. (c) The diffusing fiber with 660 nm visible light coupled in. (d) Illumination pattern on a piece of white paper of the 660 nm light from the diffusing fiber. The radial illumination pattern along the etched fiber length can be modeled as cylindrical shape. Assuming that the output energy from the fiber is uniformly distributed along the diffusing fiber, the fluence at distance from the fiber position can be estimated as, , (3.5) where is the length of the etched fiber. In the current design, the length of the etched fiber is 2 cm. A longer length of the diffusing fiber can also be fabricated using the same approach. In this model, the beam divergence along the longitudinal direction is neglected. Figure 3 shows how the outputERFluence2outputER lp=Ã Ãl42 fluence drops with increased distance for and in our design. Eq. (3.5) is a rough estimation of light fluence inside tissue where the light distribution is assumed to be homogeneous after multiple scattering events. The region close to the fiber surface has a very high fluence due to the limited contact area. With a 25 Fr cystoscope sheath, the minimum laser-tissue contact distance is estimated as ~5 mm. In Figure 3.3, the fluence is shown to be nearly 8 mJ/cm2 at R=5 mm, which represents the maximum fluence on the laser-tissue contact surface. This value is well below the MPE value of 20 mJ/cm2 (for Î»â¤700 nm). Meanwhile, the fluence drops exponentially as the distance increases. The fluence falls to 2 mJ/cm2 at distance. Due to our design of using a long section of diffusing fiber and the cystoscope sheath, the laser energy is distributed over a relatively large laser-tissue contact area and the fluence expected on the urethral wall is kept well below the laser safety limit MPE. To increase the SNR and penetration depth, approaches of increasing the fluence are also considered. By adding a cylindrical reflection mirror which directs light to just one side of the probe, the fluence can be doubled. Furthermore, nearly one third of the light still exits from the forward direction of the diffusing fiber, which will also be reflected to the side illumination direction. Figure 3.3 Calculated fluence on a cylindrical illumination surface at various distances. The output energy from the fiber is assumed to be 50 mJ/pulse. The dash lines indicate the position of the fiber surface and the tissue surface, respectively. 50 mJ/pulseoutputE = 2 cml =2 cmR =43 3.2.3 Design of side reflection mirror to reshape the beam In 360 degrees radial illumination, the fluence at far distance is relatively low. To increase the fluence, a mirror can be added along the diffusing fiber to reflect the beam to only one side and therefore double the fluence. Considering the geometrical characteristics of the urethra tract, a circular cylindrical mirror and a parabolic cylindrical mirror are tested and compared by simulation and experiment. Based on optical principles, the focus of a circular cylindrical mirror only exists under paraxial condition. Comparatively, a parabolic cylindrical mirror can form an ideal parallel reflected beam. For the diffusing fiber positioned at the focus, the reflected rays by the circular cylindrical mirror and parabolic cylindrical mirror are shown in Figure 3.4(a) and (b), respectively. In Figure 3.4(a), only the paraxial rays are collimated while the rays at large angles are slightly focused. In Figure 3.4(b), all the rays reflected by the parabolic surface become collimated. Figure 3.4 Comparison of the reflection pattern. (a) and (b) are the ray tracing for the circular cylindrical mirror and the parabolic cylindrical mirror, respectively. (c) and (d) are the simulated light distribution pattern at 5 mm away from the circular cylindrical mirror and the parabolic cylindrical mirror, respectively. 44 A Zemax (Zemax LLC, Washington, United States) simulation is carried out to analyze the intensity distribution of the reflected beam by the two mirrors. In the simulation, the parameters are set as for the circular cylindrical mirror and for the parabolic cylindrical mirror. Their focal distances are 1.25 mm and 1 mm, correspondingly. For both mirrors, the width along x-axis is 5 mm and the length along y-axis is 30 mm. With the same width, the height of the two mirrors are slightly different. A 20 mm long radial light source with 360 degrees emission angle is placed at the focus line to simulate the diffusing fiber. Figure 3.4(a) and (b) present the cross-sectional view of the simulation with circular and parabolic cylindrical mirror respectively. Figure 3.4(c) and (d) show the intensity distribution from the two mirrors on an observation plane at 5 mm away from the fiber, which is the expected location of the laser-tissue contact surface. Figure 3.4(c) shows that the reflected beam is not a perfect parallel beam but has a tendency to converge to the center. In comparison, Figure 3.4(d) shows a more uniform distribution pattern. Since the rays that are above the reflection surface are not reflected by the mirror, those rays will still be diverged and super-positioned on the reflected rays from the mirror. For our experiments, a circular cylindrical mirror and a parabolic cylindrical mirror are fabricated by 3D printing (Shapeways, New York, United States) in premium silver which provides high reflectivity over 95%. The fabrication accuracy is around 0.1 mm. The design parameters of the mirrors in the simulation and experiment are the same. The experimental characterizations of the mirrors are shown in Figure 3.5. Figure 3.5(a) illustrates the experimental setup for measuring the light intensity distribution. Visible light at 660-nm is coupled into the probe to display the beam distribution. The beam illumination on a piece of white paper and its intensity profile are recorded by a digital camera [23]. Figure 3.5(b) shows the beam intensity distribution as a function of distance for the circular cylindrical mirror and for the parabolic 2.5 mmr = 2 mmp =45 cylindrical mirror, respectively. The intensity is normalized between 0 to 1. Both the circular and parabolic cylindrical mirrors display a rectangular beam profile. The parabolic cylindrical mirror forms a more uniform distribution than the circular cylindrical mirror. The full width of half maximum (FWHM) of the beam size in the X and Y directions are plotted in Figure 3.5(c) and (d), respectively, as a function of distance away from the probe. By placing the diffusing fiber at the focus of the cylindrical reflection mirror, the divergence of the beam in the X direction is reshaped into a more collimated beam. According to optical principles, the parabolic cylindrical mirror can provide a better collimation than the circular cylindrical mirror. This can be seen in Figure 3.5(b) and (c). With the circular cylindrical mirror, the FWHM of the beam size does not change much for z less than 4 mm and then it diverges much faster than the beam reshaped by the parabolic cylindrical mirror. With the parabolic cylindrical mirror, the beam size diverges slowly. The overall gradual divergence of the beam is because the parabolic cylindrical mirror only reflects and collimates about half of the rays coming from the fiber and the rays that are not reflected by the mirror are still diverged. In the Y direction, the beam size is mainly determined by the length of the diffusing fiber although some divergence is also observed. When the distance between the diffusing fiber and tissue surface is known, the beam size on the laser-tissue contact surface can be estimated based on the beam divergence shown in Figure 3.5. 46 B Figure 3.5 Experimental results of the beam profile. (a) Illustration of the measurement setup. (b) Changes of the beam profile at various distances for the circular cylindrical mirror and the parabolic cylindrical mirror, respectively. (c) and (d) are the comparison of the FWHM of the beam profile in the x and y directions, respectively. 3.2.4 Characterization of the transurethral illumination probe To redirect the light exiting in the forward direction of the fiber to the side direction, a 45 degrees end mirror is added to the distal end of the parabolic cylindrical mirror. The 45 degrees end mirror reflects the residual forward propagating light to the side, which enhances the output energy in the side illumination. Figure 3.6(a) and (b) show the 2D beam intensity distribution with and without the end mirror, respectively, at the screen distance of 5 mm. Figure 3.6(c) and (d) compare the 1D intensity profile in the X and Y direction, respectively. The 1D profiles are extracted from Figure 3.6(a) and (b) at the corresponding positions marked by the white arrows. Since the end mirror strengthens the light intensity at one end by reflecting the light from the fiber tip, Figure 3.6(d) shows a slight non-uniform density distribution where the end mirror is added. In Figure 3.6(c), the central beam profile with the end mirror looks similar with the profile without the end mirror 47 in the x direction. The enhancement of output energy from the illumination probe is also measured by a large area detection sensor (J-50MB-YAG, Coherent, CA, United States). For example, without the end mirror, the pulse energy in the side illumination is measured to be 21 mJ/pulse. With the end mirror, the output energy is improved to 28 mJ/pulse, which represents a 33% increase. From the beam intensity distribution, the illumination area at intensity is found to be around 4 cm2 at 5 mm distance that is assumed to be the closest distance from the fiber to the tissue surface. In the following experiments, this measured beam area will be used to calculate the fluence on the laser-tissue contact surface by using the measured pulsed energy divided by the area. Due to some laser power fluctuation, the fluence value may vary accordingly. Figure 3.6 Beam profile from the parabolic cylindrical mirror with and without the 45° end mirror. (a) 2D beam profile without the end mirror (b) 2D beam profile with the end mirror (c) 1D beam profile along x axis. (d) 1D beam profile along y axis. White arrows indicate the positions where the 1D profiles are extracted. 21 / e48 3.3 PAT imaging experimental results 3.3.1 Experimental setup The PAT system consists a Q-switched Nd:YAG laser (Surelite OPO Plus, SLIII-10, Amplitude, San Jose, United States) and a medical ultrasound system (Ultrasonix MDP, Analogic, Richmond, Canada) [32]. The laser generates 3-5 ns laser pulses at a 10 Hz repetition rate and 532 nm wavelength. The Nd:YAG laser pumps an optical parametric oscillator (OPO) which can output light from 675 to 2500 nm. To couple the laser output into a fiber, a beam homogenizer is applied to reduce the peak power on the fiber tip surface. With this scheme, up to 50 mJ/pulse energy can be coupled into a 1000-µm-core-diameter fiber. For receiving the PA signal, a 5-14 MHz linear transducer array with 38 mm width (L14-5/38, Analogic, Richmond, Canada) is used for general experiment and a 5-9 MHz TRUS transducer with 55 mm width (Analogic, Richmond, Canada) is used for experiment on the prostate mimicking phantom. The received signal is amplified and digitized by a data acquisition system (SonixDAQ, Analogic, Richmond, Canada). PAT images are reconstructed by delay and sum [41]. To characterize the maximum optical power that can be coupled into a 1000-µm-core-diameter fiber, the output from a cleaved fiber without the processing of diffusion tip is measured, and a 50 mJ/pulse energy is obtained. To characterize the optical power that is delivered by the diffusing fiber, the optical power from the side illumination beam and the residual forward emitted beam are measured and summed together, and a total power of 28 mJ/pulse is obtained. Several reasons may cause the power measured from the diffusing fiber to be lower than the unprocessed fiber. The laser coupling into the diffusing fiber may not be optimal because the power is maximized for the residual forward emitted beam, where it may not be the maximum for the side illumination beam simultaneously. Since the light emitted from the diffusing fiber is diffusing over 49 a wide area, the power meter may not capture all the emitted power in all directions. Assuming all the 50 mJ/pulse can be delivered by the diffusing fiber, the fluence on the laser-tissue contact surface would be 13 mJ/cm2, which is still below the safety limit. 3.3.2 Human hairs inside chicken breast tissue The transurethral illumination probe is tested by imaging human hairs embedded in chicken breast tissue. The illumination probe and the transducer are positioned face-to-face with a separation of approximately 6 cm. The fluence on the laser tissue contact surface is around 8.5 mJ/cm2. The reconstructed PAT image of human hair is shown in Figure 3.7(a). The hair can be visualized as deep as 3 cm away from the illumination probe. Laterally, a region with 3 cm width can be detected thanks to the long length of the diffusing fiber. In Figure 3.7(a), some strong PA signals are observed near the probe surface which could be artifacts. The artifacts could be due to reverberation of PA signals or strong laser fluence near the probe. The comparatively strong local fluence near the probe surface could generate stronger PA signal from any absorbers in that region [25]. Higher fluence rate leads to higher SNR, but also generates strong artifacts near the illumination probe. Instead of masking out this artifact as reported [23], an attenuation compensation method is applied to reduce this artifact [104]. In tissue, light intensity attenuates exponentially due to scattering and absorption. Thus, the fluence can be considered to attenuate exponentially, which can be described as [105], , (3.6) where is the effective attenuation coefficient of the turbid medium and is the fluence of the incident beam. The effective attenuation is an approximation for biological tissue when the absorption is dominated by scattering and is given as [106], , (3.7) 0( ) effzI z I e µ-µeffµ 0I( )'3eff a a sµ µ µ µ= +50 where is the absorption coefficient and is the reduced scattering coefficient. To compensate the depth dependent signal, drop due to reduced fluence at deeper penetration depth, an attenuation compensation factor can be multiplied with the received PA signal. In our study, the effective attenuation coefficient of 0.9 cm-1 at 650 nm in chicken breast tissue is used according to literature reports [107]. Figure 3.7(b) shows the reconstructed PAT image after applying the attenuation compensation. The artifact near the illumination probe is effectively removed. Moreover, the contrast for the structures further away from the light source is greatly enhanced. A hair located 4.5 cm away from the illumination probe can be visualized, indicating good penetration depth of the PAT system. Figure 3.7 PAT images of human hair embedded in chicken breast tissue. (a) Without attenuation compensation. (b) With attenuation compensation. 3.3.3 Prostate mimicking phantom To evaluate the potential of the PAT system for prostate imaging, a tissue-mimicking prostate phantom was developed and imaged. The phantom was constructed to contain the specific absorption and scattering coefficients of prostate tissue close to the reported value of human prostate [108-110]. Several studies have characterized the optical properties of human prostate at various wavelengths. Zhu et al. measured the reduced scattering coefficient to be 11.8 ± 8.2 cm-1 aµ'sµeff zeµ51 (range from 1.1 to 44 cm-1) and the absorption coefficient to be 0.4 ± 0.2 cm-1 at 732 nm wavelength from 13 patients [108]. Weersink et al. measured the reduced scattering coefficient to be 3.4 ± 1.6 cm-1 and the absorption coefficient to be 0.4 ± 0.2 cm-1 at 762 nm wavelength from 22 patients [109]. Svensson et al. reported the reduced scattering coefficient to be 8.7 ± 1.9 cm-1 and the absorption coefficient to be 0.5 ± 0.1 cm-1 at 660 nm wavelength from 9 patients [110]. Thus, there seem to have a wide range and large standard deviation due to the difference among patients and the selection of measurement locations. Our phantom recipe includes 0.3 g agar powder, 0.3 g aluminum oxide, and 0.02 ml methylene blue in every 10 ml water. Using the inverse adding-doubling method [111], the reduced scattering coefficient and the absorption coefficient are measured to be 2.8 cm-1 and 0.2 cm-1, respectively. The effective attenuation coefficient and the penetration depth for 1/e decay can be calculated to be 1.34 cm-1 and 0.75 cm [110]. Compared with the values in literature, the coefficients obtained by our recipe are within the range but on the lower end. This phantom will provide a reasonable estimation of the performance of the probe in prostate tissue. The actual penetration for in vivo imaging of prostate will be patient dependent. Two tube-shaped channels are created in the phantom to mimic the rectum and urethra as shown in Figure 3.8(a). Between the two channels, four rows of pencil leads are placed in the middle as shown in Figure 3.8(b). In imaging, the TRUS probe and illumination probe are inserted in the mimicking rectum and urethra channels. The fluence on the phantom is around 10 mJ/cm2. The PAT image is shown in Figure 3.8(c), where two rows of pencil leads closer to the light source can be detected. At larger depth, the fluence is largely attenuated and the SNR is very low. After applying the attenuation compensation, all the four rows of pencil leads can be clearly seen as shown in Figure 3.8(d), indicting an optical penetration of ~3.5 cm. 52 While the reflection mirror design can increase the fluence, the current transurethral illumination probe provides only one-side illumination at each laser shot. To image the structures on both sides of the probe, the transurethral illumination probe can be rotated by 180° to illuminate the other side sequentially. This is tested on a prostate mimicking phantom with inserted human hairs. In the experiment, the position of the TRUS probe is fixed. By rotating the illumination probe at 0° and 180° respectively in the phantom, a forward and a backward illumination is sequentially applied to image the structures on both sides. Figure 3.9(a) and (b) show the human hairs obtained by forward and backward illumination, respectively. The hair can be resolved with better resolution in Figure 3.9(a) than in (b) because the objects are closer to the transducer in the forward illumination. Figure 3.9(c) shows the combined image from both the forward and backward illumination, which is similar to the case as if the fiber probe illuminates at 360°. Figure 3.8 PAT imaging of prostate mimicking phantom. (a) Photograph of the phantom where two tube shaped cavities indicate the location of rectum and urethra. (b) Experimental setup. (c) and (d) are the PAT images of the pencil leads embedded in the phantom without, and with attenuation compensation, respectively. 53 Figure 3.9 PAT images of hair embedded in the prostate mimicking phantom. (a) Forward illumination where the illumination probe faces the TRUS probe. (b) Backward illumination where the illumination probe faces the opposite side of the TRUS probe. (c) Result combining (a) and (b). For imaging the prostate, the fiber probe together with the cystoscope sheath can be inserted through urethra. The direction of the illumination beam can be controlled by rotating the cystoscope sheath to illuminate the different regions around the prostate. For 2D imaging, two rotation angles for a forward and a backward illumination will be sufficient. For 3D imaging, 3-4 rotation angles of the probe may be required in order to cover the full volume of the prostate. For receiving the PA signal, a scheme to image prostate in 3D by using a da Vinci robot to control the scanning of a linear transducer has been reported [112]. 3.4 Discussion Compared to the transurethral illumination probe reported by Bell et al. [23], our design makes several major improvements towards translating the technology for in vivo imaging. A high energy of up to 50 mJ/pulse is coupled into a 1000-µm-core-diameter fiber. The high energy pulse ensures that high SNR and deep penetration can be achieved with this probe in PAT imaging. Furthermore, the laser safety limit is addressed by providing a large illumination area and controlling the fluence on the laser-tissue contact surface to below the MPE requirement. The large illumination area is achieved by processing the fiber tip to create diffusing fiber over a few centimeters length. Thus, most of the light can be diffused to the side illumination direction. A parabolic cylindrical mirror 54 reshapes the beam to one side of the probe for controlling the fluence and illumination pattern. Therefore, both high fluence and laser safety limit are addressed, which can potentially translate PAT for in vivo prostate imaging in the future. There are several improvements that can be considered in the future to make the probe more suitable for prostate imaging in the clinic. Prostate cancer is most commonly found in elderly men, who usually have enlarged prostate caused by benign prostatic hyperplasia [100, 101]. Collins et al. reported that the mean volume of prostate increased from 25 ml for the 40-49 age group, to 40 ml for the 70-79 age group in a study of 472 patients [100]. Eri et al. reported the mean anterioposterior, transverse, and cephalocaudal prostate diameters to be 3.5, 5.3, and 4.9 cm, respectively, in a patient group with average age of 70 [101]. In order to image the entire prostate in elderly patient, a pullback approach of the fiber probe can be applied, or the length of the diffusing fiber can be increased to 5 cm directly. Fabricating a 5 cm diffusing fiber is straightforward by applying the process described in Section 3.1.2 over a longer section of fiber. The method of using etching cream to fabricate the diffusing fiber is not limited on the fiber length. The etched fiber is only slightly reduced in its diameter and the strength of the fiber is maintained. The optical properties of prostate tissue have been measured by several groups. However, the reported values for the absorption coefficient and reduced scattering coefficient vary significantly. The reduced scattering coefficient applied to our mimicking phantom is towards the lower end of the reported values. Therefore, in real prostate tissue, the penetration depth may be lower. Higher optical power could be used to increase the penetration depth. As our current measurement is done at 10 mJ/cm2, there is still space to increase the fluence within the safety limit to achieve deeper penetration for imaging the whole prostate. 55 In the prototype, the probe is inserted in a cystoscope sheath of 25 Fr. For standard diagnostic cystoscopy, 21-22 Fr rigid cystoscope is most frequently used. Intravenous sedation or regional anesthesia is required when rigid cystoscopy is applied. Larger diameter rigid cystoscopy only occurs in some particular cases which requires general anesthesia [113]. To reduce the probe diameter to 21-22 Fr, the 1000-µm-core-diameter fiber can still be used but the size of the parabolic cylindrical mirror needs to be reduced. This will decrease the beam size and result in higher fluence. To keep the fluence below the MPE requirement, the length of the diffusing fiber can be increased to illuminate over a longer beam pattern. The current beam profile is not completely uniform, especially near the tip of the fiber probe where the end mirror reflects the forward emitted light to side illumination. Non-uniform beam pattern could limit the maximum energy that can be delivered to the tissue in order to satisfy MPE requirement everywhere. To reduce the power emitted from the forward direction and to increase the uniformity of the illumination pattern, the length of the diffusing fiber can be increased, where a larger fraction of the light will be diffused to the side illumination and only a smaller fraction will emit from the forward direction of the fiber. A rough reflection surface on the distal end mirror can be processed to create diffuse reflection which can further improve the beam uniformity. The fiber probe is housed in a rigid cystoscope sheath which protects and provides support to the fiber. During insertion of the cystoscope in urethra, urine can flow in the device. Although urine will not likely affect the optical property of the probe, the probe will be in direct contact with urine. Future improvement on the probe can be achieved by adding a quartz tube to protect the diffusing fiber. Similar designs have been reported by Bell et al. in their transurethral illumination catheter [23] and by Li et al. in their internal illumination PAT [114]. Adding a quartz tubing will 56 isolate the fiber from direct contact with the tissue, which can protect the fiber from breaking. Furthermore, it will make sterilization of the probe easier and make it possible to reuse the probe after each operation. With the improved probe, PAT imaging of prostate cancer will be carried out in the future. PAT imaging will be used to detect the neurovasculature bundles and microvessel structures in prostate. Detecting the extension and proximity to critical structures such as the neurovascular bundles during radical prostatectomy is critical for ensuring surgical outcomes and maintaining normal organ functions. Detecting and quantifying the microvessel structures can potentially differentiate prostate cancer from normal tissues. Previously, the functional US and MRI techniques for quantifying MVD have shown limitations due to low sensitivity and difficulty in discriminating cancer from high vascularized benign prostate hyperplasia [100, 101]. Rubin et al. assessed MVD in one hundred radical prostatectomy specimens by comparing pathologic result with intra-observer MVD counting result [115]. The study found that MVD was not associated with Gleanson sum, tumor stage, surgical margin status or seminal vesicle invasion. In a study by Mucci et al. [116], the histopathologic measure of MVD on serial sections of the prostatectomy from 572 patients diagnosed with prostate cancer showed that MVD was not linked to cancer-specific mortality. Instead, the vascular size and irregularity was found to be strong predictors of the fatal disease. With its high sensitivity, PAT could provide more information about the microvessel structures. Horiguchi et al. analyzed the correlation of PA signal intensity with the MVD, total vascular area (TVA) and total vascular length (TVL) measured by histopathology in three prostate cancer patients [24]. The histological examination showed that TVA and TVL were larger and longer in the peritumoral area of the tumor than in the index tumor. The areas with high PA 57 intensity correlated well with regions with greater TVA and TVL. High PA intensity correlated with high MVD region in one patient. Furthermore, other quantitative analysis can also be carried out based on the radio frequency (RF) spectrum of PA signal because the PA signal frequency is sensitive to the size of the absorbers. Sinha et al. and Huang et al. have demonstrated the feasibility of identifying the absorber size by analyzing PA signal in frequency domain [117, 118]. Evaluation of PA signal in frequency domain might show useful information about the size of microvessels. Therefore, PAT has great potential in improving prostate cancer imaging by detecting and analyzing multiple characteristics of the microvessels, including MVD, area and length of microvascular structure, and the size of microvessels. 3.5 Summary A transurethral illumination probe was designed for PAT imaging of prostate. The probe contains a radial illuminating diffusing fiber over 2 cm length. The radial illumination is reflected and reshaped by a parabolic cylindrical mirror to obtain near parallel side illumination with doubled fluence. The fiber assembly is inserted into a 25 Fr cystoscope sheath to provide protection of the fiber and a minimal laser-tissue contact distance of 5 mm. Due to the unique design of the illumination probe, a 4 cm2 laser-tissue contact surface is provided and the fluence is controlled to be below the American National Standards Institute (ANSI) limit. The illumination probe was evaluated using tissue phantoms and prostate mimicking phantoms. A 3.5 cm penetration depth in the prostate mimicking phantom was demonstrated using 10 mJ/cm2 fluence at 700 nm wavelength. An approach of acquiring PAT image from both sides of the probe was also demonstrated by rotating the illumination probe at 0° and 180° for a forward and a backward illumination, respectively. The deep penetration ability makes PAT of the prostate with transurethral 58 illumination feasible. This technique can potentially be used for PCa detection and image-guided targeted biopsy. 59 Chapter 4: Variance-reduced stochastic gradient descent algorithm for photoacoustic image reconstruction 4.1 Introduction As described in previous chapters, PAT achieves deep optical penetration and high ultrasound resolution, which brings plenty of potential applications in biomedical imaging [8, 9, 22-26, 58, 59]. Since PAT is based on the detection of acoustic wave, its imaging reconstruction is realized by collecting the acoustic signal through transducers with various detection geometries. Based on the spatial distribution of the transducers, the detection geometry can be categorized as full-view detection and limited-view detection. In full-view detection, the transducers are placed or scanned completely around the object. In 2D, full-view detection can be obtained by using a full-ring transducer array [119], circularly scanning the transducers [30, 120] or rotating the sample [121]. In 3D, the transducers should be placed or scanned around a complete sphere which is non-practical due to the inhomogeneous complex structure of biological organs. On the other hand, the limited-view detection is very common in PAT applications, e.g. breast imaging with a semi-spherical transducers array for detection [26]. Current customized full-ring or arc-shaped transducers collects data with relatively large detection angles for the reconstruction, but each of them focuses on one specific application because the dimension of the array relates to the size of the sample. Meanwhile commercial linear array transducers are commonly used in clinical ultrasound imaging and widely studied in PAT [31, 112, 121]. Linear array transducers suffer from limited view problem. To overcome the limited-view problem, various methods have been proposed to increase the detection angle [31, 112, 122]. Our group demonstrated a method by using two linear array transducers to improve the detection view in 2D [31]. Images obtained from the two transducers can be merged by the registration and calibration of the transducers' positions. 60 Meanwhile, Moradi et al. presented a 3D imaging approach by using a linear pick-up transducer array scanned on a cylindrical surface by a robotic arm [112, 122], which also improved the image quality by enlarging the detection view. Other than the setup of the transducers, reconstruction algorithm also determines the image quality of PAT. Exact reconstruction algorithm is only possible when the detection surface is enclosed or infinitely unbounded [41, 123-126]. For the linear array transducer, many approximate algorithms are performed in either time domain or frequency domain [127-129] based on an analytical model. As one of the most simple back-projection algorithms, delay-and-sum (DAS) uses each detection position as the center to project the measured PA signal back to the image space as spherical shells and then add them together [130]. Even though the approximate analytical reconstruction algorithms are known for their ease of implementation and low computation cost, artifacts are significant when the detection view is limited, due to insufficient measurements and noise. Without using an explicit priori reconstruction model, deep learning is an emerging technique using artificial neural network to represent nonlinear input-to-output mapping based on the training data [131-133]. Antholzer et al. demonstrated the feasibility of applying deep learning for PAT. Their simulation results showed that the deep learning method achieved equivalent image quality as iterative method but with less computational time [131]. However, in their study the PAT images were required to share similarity with the training data used to adjust the convolutional neural network. Deng et al. applied two deep-learning-based methods, U-net and base model method, to improve the image quality under limited-view detection, and they showed effective suppression of artifacts [132]. Recently, Liu et al. presented a dictionary learning 61 reconstruction for PAT. Their results showed a potential to use fewer transducers to achieve equivalent imaging quality [133]. Iterative reconstruction methods, on the other hand, relies on an accurate forward PA propagation model and the prior knowledge of the reconstructed PA image, such as sparsity [134-142]. Several sparsity representations have been used in reconstruction models such as total variation (TV), Fourier and Wavelet transformation. Compressed sensing is a model-based technique that incorporates sparsity to improve the reconstruction from insufficient measurement. Provost et al. proposed a reconstruction based on compressed sensing for PAT with circular array transducer [134]. Lin et al. implemented the method to linear array transducer for synthetic aperture PAT and reduced the number of measurements [135]. Guo et al. demonstrated compressed sensing-based PAT to image rat brain in vivo and their result showed the effectiveness in reducing artifacts with under sampling [136]. At the expense of reconstruction time, iteration improves the image quality and reduces artifacts. The core of an iterative reconstruction algorithm is to minimize a cost function that usually contains two terms, a misfit term and a regularization term. The misfit term is based on the least-squares error between the projected PA data and the measured PA data. In general, the regularization parameter is a TV term which can prevent over-fitting issue, make the images smooth, and enhance the robustness of the algorithm [138-141]. For minimizing the cost function, full gradient descent (FGD) algorithm computes the full gradient of the cost function by taking the full projection matrix into consideration. However, the computational burden would be huge when the size of the PA data or the image size is large. Paltauf et al. used FGD on relatively small size of simulation and experimental PA data [137]. Since they did not apply a regularization term in the cost function, their method was affected by the insufficient and noisy PA data. Stochastic 62 gradient descent (SGD) algorithm randomly chooses a smaller size of sub-matrix of the projection matrix to compute the gradient and therefore can significantly reduce the computational burden. Zhang et al. used SGD algorithm on PAT where the sub-matrix was selected sequentially from each detection view [138]. Although the SGD algorithm is less costly, the convergence rate of SGD is very poor because the gradient decent direction computed from the sub-matrix in SGD has large variance from that computed from the full matrix in FGD [143]. In order to reduce the effect of the large variance, one approach is to gradually decrease the step size of the gradient decent algorithm [137]. However, this approach requires a careful tuning of the step size and it can only achieve a sub-linear convergence rate. In recent years, variance-reduced stochastic gradient descent (VR-SGD) algorithms have been developed so as to achieve a faster convergence rate by reducing the variance between the stochastic gradient direction and the full gradient direction. VR-SGD algorithms have been applied to various problems, such as machine learning and computed tomography (CT) image reconstruction, and have been shown to have much faster convergence than the SGD schemes [144]. To the best of our knowledge, VR-SGD has not been applied to PAT. In this chapter, VR-SGD is applied to the iterative reconstruction algorithm for PAT. Simulations based on different detection geometries of linear array transducers are conducted. Also, the effect of the number of transducer elements is discussed. Experimental data obtained by single linear array transducer and two combined linear array transducers are reconstructed by the proposed algorithm. 4.2 Method 4.2.1 PAT imaging model In order to apply the iterative PAT image reconstruction method, a discrete model which describes the relationship between the measured PA signal and the optical absorption distribution map needs 63 to be built. The model discussed in this chapter focuses on 2D plane. In response to an optical absorption, the pressure at position r and time t obeys the following wave equation [27], , (4.1) where is the isobaric volume expansion coefficient, is the specific heat and is the speed of sound. is a heating function which is the product of the spatial absorption and a temporal illumination function which can be regarded as an infinitely short pulse . Eq. (4.1) can be solved by using the time-retarded solution as, p(ð, ð¡) = ð½4ðð¶" ððð¡ ; ð´(ðâ²)ððâ²|ð%ð&|'() (4.2) By acting on both sides of Eq (4.2), an integrated data function is obtained as[145], g(ð, ð¡) = ð¡@ ð(ð, ð¡â²)ðð¡â²)* = ð½4ðð¶" ððð¡ ; ð´(ðâ²)ððâ²|ð%ð&|'() (4.3) which represents a scaled acoustic velocity potential [146]. The function is treated as a measurable and known function since it can be derived from and vice versa. Based on the relation between the integrated data function and the optical absorption distribution, a forward propagation operator can be obtained. Next, the optical absorption map will be described as . It can be discretized into a 2D matrix, which refers to m rows and n columns. It is then lexicographically reshaped into a one-dimensional (1D) vector by arranging each column head-to-tail sequentially, leading to a 1D matrix representation which has the size of . A linear array transducer has elements ( , )p tr222 21( , ) ( , ) ( , )pp t p t H tc t C tb¶ ¶Ã - = -¶ ¶r r rb pC c( , )H tr ( )A r( )td0'tt dtò ( , )g tr( , )g tr( , )p tr( )x rm n´X 1mn´ Q64 and the length of the discrete signal for each transducer channel is . For each single transducer element, the discrete version of the model in Eq. (4.3) can be expressed as [138], , (4.4) where represents the measured data by the ith transducer element. The projection matrix has the size of and can be built as [138], , (4.5) where is the distance from the kth pixel on the object to the ith transducer element, and is the time step which equals to the reciprocal of the sampling rate. ð·(ð, ð) is the directivity weighting factor based on the relative position between the kth pixel and the ith transducer element. Hence the full projection map can be expressed as, , (4.6) where the forward projection operator can be written as, , (4.7) The full size of the operator is . The size of the operator depends on the size of the image, the number of transducer elements and the sampling rate. Even for an image with pixels, the size of the operator can be over thousand by thousand. 4.2.2 Variance-reduced stochastic gradient descent method for solving the cost function In the presence of noise, image reconstruction problem with given refers to solve a least-squares cost function [147, 148], ð½(ð) = 12âð´ð â ðºâ++ + ððð(ð) (4.8) T[ ], 1,i iG A X i Q= ÃiG iAT mn´{ } ( ) [ ]( , ) max 1 ( ) ,0 D , , 1,iA j k d c t j i k j T= - D - Ã Ãd X tDX G!G AX=A1 2 ...TQA A A AÃ© Ã¹= Ã« Ã»TQ mn´32 32´G65 where is the entry-wise L2-norm. The existence of the TV term is a means of regularization to avoid over-fitting and is the corresponding tuning parameter. Based on the knowledge of the sparsity of the object, the gradient matrix of the optical absorption map is supposed to be sparse. Thus, the total variation of the object can be used as the regularization term in PAT image reconstruction. The calculation of the total variation can be described as, , (4.9) where refers to the pixel belonging to the ith row and jth column in the 2D image plane. The solution of the estimated object can be denoted as, . (4.10) To find the solution, the FGD algorithm calculates the full gradient of the cost function and gradually updates the value of estimated object along the negative full gradient direction. Although the updating speed can be slightly adjusted by changing the step size, the calculation of the full gradient for the cost function, especially the forward projection operator , takes a long time. In each loop of iteration, the multiplication of and its transpose has to be performed which is time consuming. On the other hand, the SGD algorithm replaces by its sub-matrix in the cost function for each loop of iteration. As the size of is much smaller than the size of , the computational time for each iteration is reduced. However, the variance between the gradient of the sub-matrix and the full gradient could be high [143], which reduces the convergence speed even though with a faster updating rate. In order to balance the tradeoff between convergence rate and updating speed, a VR-SGD algorithm for PAT image reconstruction is proposed as shown in Algorithm 1. The main idea of 2Ãl( ) ( ) ( )2 2, 1, , , 12 2m ni j i j i j i ji jTV X X X X X- -= == - + -åå,i jX*X( )2* 21argmin2XX AX G TV Xlì ü= - +í ýî þAAA iAiA A66 the proposed algorithm is to keep the direction of full gradient during the iteration and a fast updating speed simultaneously. To achieve this goal, the full gradient is updated in the outer loop after Q iterations of the inner loop. By reducing the frequency of computing the full gradient, the updating speed is enhanced. Meanwhile, it also guarantees the updating direction to be consistent with the negative direction of the full gradient without significant deviation. The iteration of the inner loop is to use the knowledge of the full gradient and the difference between the gradient of the currently estimated and the recent updated from the outer loop to build the updating direction which can be expressed as, . (4.11) Unlike the updating described in SGD which only uses the information of a single sub-matrix, sub-matrixes are chosen as a selected set of size to calculate the descent direction. When is larger than 1, the method can reduce the variance of the update direction without calculating the full gradient. The sub-matrixes are chosen based on the probability proportional to their Lipschitz constant, which proves to reach a better convergence [144, 149]. For each sub-matrix , its Lipschitz constant is defined as the largest eigenvalue of . When the direction is provided, the new estimated object in the inner loop can be calculated by using a proximal operator which minimizes the total variation of the updating object. As the result obtained by the iteration gradually converges, the updating step size is also decayed exponentially by after each loop of outer iteration. With a given initial step size and decay constant , the updating process can be expressed as [150, 151], . (4.12) ( )J XÃ !( )J XÃ !kjX X!( ) ( )1kT ki i ji Sd A A X X J XM Ã= - +ÃÃ¥ ! !MkS MMiA Ti iA Aaa a( )1k kj TV j jX X dl a-= -prox67 where denotes the proximal operator used to solve the differentiation of the TV term. is the updated image after jth time of the outer iteration and th time of the inner iteration. is the one updated based on in the inner iteration, respectively. The updating step size in each outer loop of iteration is related to the initial step size following . As the initial state of the iteration affects the whole optimization process, either zero matrix or the reconstructed image by DAS can be used. Since the value from the reconstructed result by DAS is closer to the optical absorption map than the zero matrix, it is chosen as the initial input to the iteration algorithm in the numerical simulations. Algorithm1: VR-SGD algorithm Data: , , , N, for j = 1 to N if j>1 then else end Full gradient direction for k = 1 to Q Select set of size from VR gradient direction end end TVlprox1kjX- ( )1k -kjX1kjX-ja1jj aa a-=0X l M akjX X=!1 /j j aa a -=0X X=!ja a=( )TA AX Gµ = -!!0jX X= !kS M { }1......Q( ) ( )11kT ki i ji Sd A A X X J XM-Ã= - +ÃÃ¥ ! !( )1k kj TV j jX X dl a-= -prox68 4.3 Numerical simulation To verify the proposed algorithm, simulation on several 2D detection geometries are conducted, including linear array, parallel-shape, L-shape and square-shape detections. The transducer configurations are shown in Figure 4.1. In each configuration, the total number of transducers is equivalent. In the linear array detection as shown in Figure 4.1(a), all the number of transducers are placed in one linear array. In the parallel-shape detection as shown in Figure 4.1(b), the depicted transducers are placed in two parallel lines where each side has number of transducers. In the L-shape detection as shown in Figure 4.1(c), two lines of transducers are placed in the perpendicular direction where each side has number of transducers. The L-shape detection can increase the detection view than a linear array detection [31]. In square-shape detection as shown in Figure 4.1(d), the transducers are placed to enclose the region of the sample by forming a square-shape detection where each side has number of transducers. The main parameters used in the simulation are listed in Table 4.1. The physical size of the image is set to be and the corresponding image pixel density is . For each phantom, the four configurations depicted in Figure 4.1 are applied. The sampling rate is 12 MHz and the total number of transducers is set to be 16, 32, 64 and 128 respectively. In the simulation, 30 times of iteration is done for each case. The TV tuning parameter is , the batch size is 5 and the initial step size is 5 with a decay rate at 1.05 for each loop of outer iteration. All the simulations are carried out in Matlab R2018a (Mathworks, Natick, United States) on a laptop with 2.3GHz Intel Core i5 CPU and 16 GB memory. N/ 2N/ 2N/ 4N38 mm 38 mm´ 128 128´Nl 51 10-Â´69 Figure 4.1 Imaging configuration by using various detection geometry with equivalent number of transducers. (a) Linear array detection. (b) Parallel-shape detection. (c) L-shape detection. (d) Square-shape detection. Table 4.1 Parameter settings for the simulation of VR-SGD algorithm Parameter Value Size of the region 38 mm Ã 38 mm Pixels 128 Ã 128 Speed of the sound 1500 m/s Number of elements 16, 32, 64, 128* Sampling rate 12 MHz (TV tuning parameter) M (batch size) 5 N (number of iteration) 30 (initial step size) 5 Decaying constant 1.05 * The number of transducer elements is selected during the simulation l 51 10-Â´aa70 A tree branch phantom as shown in Figure 4.2(a) is used to simulate blood vessels. In the simulation, the measured signal is simulated by using the projection matrix described in Eq. (4.6) in a 2D plane. A noise relative to its maximum intensity is added to . Figure 4.2(b) shows the reconstructed image by DAS for the square-shape detection with 128 elements. Here only the result in the square-shape detection is shown which has better image quality than the results obtained in the other configurations by DAS. Comparatively, Figure 4.3 shows the reconstruction results by the VR-SGD algorithm in different configurations. Figure 4.3(a)-(d) show the reconstructed images by the linear array configuration. Some features oriented perpendicularly to the transducers are missing or blurred. Moreover, the intensity distribution of the tree branch is not homogeneous since features located on the right half of the image are farther away from the transducer. Figure 4.3(e)-(h) show the reconstructed images by the parallel-shape detection. Compared to linear array, more features are distinguishable, and the intensity distribution looks more balanced in the parallel-shape detection. Figure 4.3(i)-(l) show the reconstructed images by the L-shape detection. Compared to the parallel-shape detection, L-shape detection benefits from a larger view angle. In Figure 4.3(l), almost all the artifacts caused by the limited-view issue have been removed after the iteration except for the region located to the bottom right where the sample is farther away from the transducers in both directions. At last, Figure 4.3(m)-(p) show the reconstructed images by the square-shape detection. The square-shape detection provides a full-view detection where the image reconstructed shows the best image quality when the total number of transducer elements are sufficient (e.g. â¥32). Interestingly, when the total number of transducer elements are limited to 16, the linear array actually performs slightly better due to the denser transducer array. In each detection configuration, when the number of transducer elements increases, the contrast of the reconstructed images gets higher and the artifacts G A10 dB G71 caused by the back-projection arcs are reduced. Compared with the image reconstructed by DAS with 128 transducer elements (shown in Figure 4.2(b)), the VR-SGD reconstructed images show much lower background, higher SNR, and significantly reduced artifacts. Figure 4.2 The tree branch phantom and its reconstruction by DAS. (a) The tree branch phantom. The red arrows indicate the location and orientation of the line profile discussed in the following paragraph. (b) Reconstruction result by DAS in the square-shape configuration with 128 transducers. In addition, the line profiles of the reconstructed image for each detection geometry with 128 transducer elements are plotted and compared with the initial pressure distribution in Figure 4.4. The location and orientation where the line profile is taken is indicated in Figure 4.2. Compared to DAS, VR-SGD can make the line profile sharper which means it can reconstruct higher resolution image. Moreover, image reconstructed by the iteration algorithm has lower background which leads to higher SNR. Especially, for the L-shape and square-shape detection geometries, the reconstructed background noise is significantly lower than the linear and parallel-shape detections. 72 Figure 4.3 The tree branch phantom reconstruction results by VR-SGD algorithm. The first to forth rows refer to the image from the detection geometry by linear array (a-d), parallel-shape (e-h), L-shape (i-l) and Square-shape (m-p). From the first column to the fourth column, the image is reconstructed by 16, 32, 64 and 128 transducer elements correspondingly. The red lines at each row indicates the orientation of the transducers array in each configuration. 73 Figure 4.4 The line profiles of the reconstructed image. (a) The horizontal line profile indicated as line 1 in Figure 4.2(a). (b) The vertical line profile indicated as line 2 in Figure 4.2(a). Since the tree branch phantom only provides line structures, a FORBILD head phantom is also simulated as a more complicated phantom [152]. Figure 4.5(a) shows the details of the FORBILD head phantom. The reconstructed image by DAS with 128 transducer elements in a square-shape detection is shown in Figure 4.5(b). The DAS reconstructed image has very low quality with missing details and a strong background. The reconstructed images by VR-SGD with the different detection geometries are displayed in Figure 4.6. Similar trend is also observed in the FORBILD head phantom. The linear array detection has the most artifacts and missing structures as shown in Figure 4.6(a)-(d). The parallel-shape detection shows more features and higher uniformity of signal than the linear array as shown in Figure 4.6(e)-(h). However, only features which have the same orientation as the transducers are well reconstructed, while some features 74 perpendicular to the transducers are still missing. In Figure 4.6(i)-(l), the L-shape detection performs reasonably well and can reconstruct most of the features. The square-shape detection performs the best. Compared with the tree branch phantom, a greater number of transducer elements are required in order to suppress the noise and artifacts as the structure of the phantom gets more complicated. Nevertheless, the images reconstructed by VR-SGD show much higher contrast than DAS in all of the four configurations. Figure 4.5 The FORBILD head phantom and its reconstruction by DAS. (a) The FORBILD head phantom. The red arrows indicate the location and orientation of the line profile discussed in the following paragraph. Line 1 and 2 are horizontal line profiles and line 3 is a vertical line profile. (b) Reconstruction result by DAS in the square-shape configuration with 128 transducers. The corresponding line profiles of the images presented in Figure 4.6 with 128 transducer elements are plotted and compared in Figure 4.7. Unlike the tree branch phantom, there are many continuous regions in the FORBILD head phantom. All the images obtained by VR-SGD show smoother and more accurate result than that by DAS. The edge of the boundary for each region is sharper when the square-shape or the L-shape detection is applied. Figure 4.7(a) shows the horizontal line profile in the top of the reconstructed images. Results obtained by the linear and the parallel-shape detection can resolve the features, but the value deviates from the true value of 75 the initial pressure. Comparatively, results obtained by the L-shape and square-shape detection match well with that of the initial pressure. Among the discussed detection geometries, the square-shape detection leads to the most precise result. The accuracy of the result given by the L-shape detection depends on the sample's location. In Figure 4.7(b), it shows the trend of a decreasing accuracy for the line profile given by the L-shape detection when the location is farther away from the transducer. The horizontal line profile obtained by the L-shape detection matches better with that of the initial pressure on the left side of the image due to the closer distance between the sample and the transducers. Similarly, the vertical line profile plotted in Figure 4.7(c) shows that the top part of the line profile provided by the L-shape detection is more precise than the bottom half as the transducers are placed on the top. Furthermore, the convergence rate and the accuracy of the iteration algorithm is evaluated by comparing the root mean square error (RMSE), which is expressed as, , (4.13) where and denotes the reconstructed and the actual optical absorption distribution, respectively. Figure 4.8(a) compares the RMSE of VR-SGD based PAT reconstruction for the FORBILD head phantom in the different detection geometries. The square-shape detection shows the smallest RMSE. The square-shape detection and the L-shape detection show a faster convergence than the other two detections. Moreover, the RMSE obtained by the L-shape detection is close to the RMSE obtained by the square-shaped detection in terms of the saturated value and convergence rate. Figure 4.8(b) compares the RMSE for the square-shape and the L-shape detection with various number of transducer elements in 30 iterations. The number of transducer elements also affects the convergence rate. As the number of elements increases, the RMSE ( ) ( ) 21 11 , ,m ni jRMSE X i j R i jmn = == -é Ã¹ë Ã»ååX R76 converges faster and saturates at a smaller value. Also, the difference of RMSE between the square-shape and the L-shape detection gets smaller when the number of elements increases. Figure 4.6 The FORBILD head phantom reconstruction results by VR-SGD algorithm. The first to forth rows refer to the image from the detection geometry by linear array (a-d), parallel-shape (e-h), L-shape (i-l) and Square-shape (m-p). From the first column to the fourth column, the image is reconstructed by 16, 32, 64 and 128 transducer elements correspondingly. The red lines at each row indicates the orientation of the transducers array in each configuration. 77 Figure 4.7 The line profiles of the reconstructed image. (a) The first horizontal line profile indicated as line 1 in Figure 4.5(a). (b) The second horizontal line profile indicated as line 2 in Figure 4.5(a). (c) The vertical line profile indicated as line 3 in Figure 4.5(a). 78 Figure 4.8 RMSE of VR-SGD based PAT reconstruction for the FORBILD head phantom. (a) RMSE for the reconstruction in different detection geometry with 128 transducers. (b) RMSE for the reconstruction by the L-shape and the square-shape detection with different number of transducers. Overall compared to DAS, the VR-SGD algorithm shows significant improvement in image quality for the PAT image reconstruction regardless of the detection geometry. The square-shape detection leads to the best reconstruction results among all the described detections due to its full-view geometry. Also, compared with the linear array detection, the parallel-shape and the L-shape detection both show improvement by providing information from another view. The L-shape detection is more advantageous than the parallel-shape detection because it covers more view angles than the latter one does. Although full-view detection is the ideal geometry for the PAT image reconstruction, it requires specialized and customized ring-shape transducers arrays which is not commonly available [153-155]. As linear array transducers are widely used in clinical ultrasound imaging, implementation based on linear array is more approachable and feasible in real applications. A single linear array has limitations in its performance. Two linear arrays in the parallel-shape geometry can also provide good image quality if the arrays are long enough and 79 cover most detection angles of the sample. The L-shape geometry is the optimum configuration with two linear arrays. We have developed an approach to calibrate the relative positions of two transducer arrays [31]. Thus, the positions of the two transducers in the L-shape detection do not need to be predetermined and can be adjusted to fit the geometry of the sample. The L-shape detection is considered as a good balance which increases the view angle and keeps its flexibility to be accessible to different sample geometries. Therefore, there is great potential to apply the L-shape configuration and VR-SGD algorithm for PAT reconstruction. 4.4 Experiments The VR-SGD algorithm is further applied and verified by PAT experiment using a single or two linear array transducers. The details of the setup of the system can be found in our previous study [31, 32]. The system includes a laser source, an optical parametric oscillator (OPO), a data acquisition (DAQ) module, an ultrasound imaging system, and ultrasound linear array transducers. The laser is a Q-switched Nd:YAG laser (Surelite OPO Plus SLIII-10, Continuum, San Jose, United States). The laser wavelength can be tuned from 680 nm to 2500 nm by the OPO. The laser pulse width is 5 ns and the repetition rate is 10 Hz. Two identical linear array transducers are used and each of them has 128 channels with 7.2 MHz center frequency, minimum 70% fractional bandwidth (at â6dB), and 0.3 mm element pitch. The probes are connected to an ultrasound imaging system (SonixMDP, Analogic, Richmond, BC, Canada) through a DAQ module (SonixDAQ, Analogic, Richmond, BC, Canada). The DAQ module acquires and digitizes the pre-beamformed radio-frequency signal from all the channels simultaneously at a sampling rate of 40 MHz and 12-bit resolution. For the experimental data, a normalized DAS reconstruction image or all zero-valued image is used as the initial state for VR-SGD. The acquired PA signal is normalized and processed by Hilbert transformation before performing VR-SGD algorithm. For the VR-SGD 80 algorithm, the batch size is 10, the initial step size is 0.1 and is set to be 1´10-4 for reconstruction of all the experimental data. 4.4.1 PAT imaging by a single linear array transducer VR-SGD is applied to PAT imaging with a single linear array transducer. Figure 4.9 shows the PAT images of a printed dots array detected by a single linear array transducer. A photograph of the phantom is shown in Figure 4.9(a). The diameter of the dots is 0.1 mm which can be treated as ideal point sources for PA imaging. The horizontal and vertical gap between two neighboring rows or columns is around 1 mm. The sample is in 2D plane and the laser illumination shines from the top of the sample plane. The single linear array transducer is aligned in the sample plane and detects PA signals from one side of the phantom. The image reconstructed by DAS from the normalized PA signal ð(ð, ð¡) is displayed in Figure 4.9(b). In the VR-SGD algorithm, the iteration is initialized by the image from DAS. During the iterations, the resulted image is forward projected and compared with g(ð, ð¡) which is the Hilbert transformation of ð(ð, ð¡). The reconstructed image after 25 times of iteration by the VR-SGD algorithm is shown in Figure 4.9(c). Comparing the images before and after applying the iteration algorithm, the image after iteration shows improved reconstruction with reduced artifacts and cleaner background. Most of the dots are reconstructed and the dots show smoother outline without noticeable artifacts. The circular shape is compressed in the direction perpendicular to the linear array transducers due to limited-view detection. Also, the background of the image is clearer which indicates higher SNR. On the other hand, the TV regularization suppresses some information which has relatively low intensity and small size. For applications that rely on this part of information for spectral analysis, the suppressing effect due to the TV regularization needs to be further investigated. l81 Figure 4.9 PAT images of printed dots array detected by single linear array transducer. (a) Photograph of the phantom. (b) Initial image by DAS based on ð(ð, ð). (c) Corresponding reconstruction result after VR-SGD for 25 times iteration. Figure 4.10 shows the cost function versus the number of iterations for the data shown in Figure 4.9. The tuning parameter is selected as λ = 1Ã10-4 or λ = 5Ã10-4 , and the step size is selected as α = 0.025, 0.05, or 0.1, respectively. There are 25 iterations for the outer loop and each outer loop includes 128 iterations of the inner loop. The initial image input obtained by DAS is given as the initial value at the beginning of each curve. To make the results comparable, the cost function is normalized as follows ð½(ð)â = â 12 âð´ð â ðºâ++ + ððð(ð)12 âðºâ++ . (4.14) In Figure 4.10(a), a larger step size leads to a faster convergence of the misfit term. However, the step size needs to be restricted in an appropriate range to avoid diverging situations. In Figure 4.10(b), a larger l results in smaller TV value. Figure 4.10(c) shows the total value of the cost function which combines the misfit term and the multiplication of the TV value with the tuning factor . Since l is comparatively small, the value of the cost function is dominated by the l82 misfit term. In this result, the combination of λ = 1Ã10-4 and α = 0.1 gives the fastest convergence. For other data sets, the selection of the parameters might vary. Figure 4.10 Comparison of the cost function value versus the number of iterations under different tuning factors and step sizes. 0 iteration means the initial input by DAS (a) Value of the misfit term. (b) Value of the TV term. (c) Value of the total cost function. There are 25 times of outer iterations. For each loop of outer iteration, 128 times of inner iterations is processed. Figure 4.11 shows the comparison of the convergence rate of FGD, SGD and VR-SGD algorithms, respectively. Figure 4.11(a) shows the convergence rate versus the number of times that the gradient updates. The number of times that the gradient updates equals to the number of iterations in FGD and SGD. In VR-SGD, there are an outer and an inner iteration loops, respectively. Thus, the number of times that the gradient updates equals to the product of the number of outer loops and the number of iterations processed due to the inner loops. For example, if 25 times of iteration has been run for the outer loop and each loop includes 128 times of update for the gradient of the submatrix, the number of times that the gradient updates is 3200. Since FGD includes the gradient of the full projection matrix, its convergence rate is the fastest at the beginning. In contrast, SGD has a very poor convergence rate. The convergence rate of VR-SGD is slower than that of FGD at first. Thanks to the stable convergence rate, the convergence rate of VR-SGD catches up with that of FGD later. Meanwhile, the convergence rate versus time is shown 83 in Figure 4.11(b) providing the approximate time each algorithm costs. For 3200 times that the gradient updates, SGD only takes around 600 s while FGD takes nearly 6000 s, which means FGD has a faster convergence rate at the expense of computational time and vice versa for SGD. For VR-SGD, it takes a balance between convergence rate and computational time by adjusting the number of submatrices used in computing the gradient. Figure 4.11 Comparison of the convergence rate for FGD, SGD and VR-SGD. (a) Normalized cost function versus the number of times that the gradient updates. (b) Normalized cost function versus the time each algorithm costs. For each algorithm, it updates the corresponding gradient 3200 times. The projection matrix currently used in VR-SGD is limited to 2D. In reality, PA signal is generated and propagates in 3D space. The PAT imaging is confined in a 2D slice with some thickness determined by the focusing of the transducer in the elevation direction. Therefore, the 2D projection matrix as the forward model may have some limitations. For the printed dots array, it is a 2D structure whose forward propagation is described by the projection matrix properly. However, for 3D structures, the projection matrix only takes account of the signal in the 2D plane of the transducers. This difference may cause problems for both DAS and VR-SGD reconstruction. To compare PAT imaging of a 3D structure, Figure 4.12 shows PAT images of pencil lead array 84 which is a 3D structure. Figure 4.12(a) shows the photograph of the pencil lead array. The pencil lead array is inserted into a gelatin phantom. The diameter of each lead is 0.7 mm and the gap between two neighboring leads is 7 mm. The illumination and detection is performed in the cross-sectional plane of the pencil leads. The PA signal is generated and propagates in 3D space. The linear array transducer mainly receives the in-plane signal. However, some out-of-plane signal can also be received by the transducer which can affect the image reconstruction. The initial image is shown in Figure 4.12(b) by the DAS image based on ð(ð, ð¡). The corresponding images after iteration by VR-SGD is shown in Figure 4.12(c). The image after iteration shows improved quality in reduced artifacts and cleaner background. The improvement is less significant compared with the results shown in Figure 4.9, where the pencil leads is a 3D object and the printed dots is a 2D object. The iteration couldnât remove all the artifacts, likely because of the out-of-plane signal. To achieve a better reconstruction for 3D structures, the out-of-plane signal should be considered which requires a projection matrix in 3D. Figure 4.12 PAT images of pencil leads array detected by a single linear array transducer. (a) Photo of the phantom. (b) Initial image reconstructed by DAS based on ð(ð, ð). Figures (c) The reconstructed image after VR-SGD. 85 4.4.2 PAT imaging by two linear array transducers Since the L-shape detection is shown to be a better detection geometry than the linear array, experiments with two linear array transducers in L-shape detection are also performed [31]. Details about the experimental setup and approaches can be found in Ref. 31. The two linear array transducers are assumed to be located in the same imaging plane after visual inspection and manual alignment. The relative positions of the two probes can be measured and calibrated [31]. The space of the reconstructed image is defined in the coordinate system of the first transducer array. The coordinate system of the second transducer array can be calibrated and transformed to the coordinate system of the first array. The transformation can be described by two translational parameters and in the X and Y directions, respectively, and one rotational parameter . To obtain those parameters of the transformation matrix, the relative positions of the two transducers are calibrated by transmitting ultrasound signal from one transducers array while receiving the ultrasound signal with the other [31]. The global rigid transformation matrix can be written as [31, 156], . (4.15) where denote the coordinate values of the second transducer in its own coordinate system, and are the coordinate values of the second transducer transformed to the imaging space (defined in the coordinate system of the first transducer). Based on this process, the projection matrix of the transducer arrays corresponding to the image space can be derived by combining Eqs. (4.5) and (4.15). In the VR-SGD algorithm, the coordinate values of the second linear array are transformed to the space of the optical absorption map. The forward projection operator is then h k q1 ''1 0 cos -sin 00 1 sin cos 01 0 0 1 0 0 1 1ii hj k jq qq q- é Ã¹æ öé Ã¹ é Ã¹ é Ã¹ê úç Ã·ê ú ê ú ê ú= ê úç Ã·ê ú ê ú ê úê úç Ã·ê ú ê ú ê úë Ã» Ã« Ã» Ã« Ã»è ø Ã« Ã»' '( , )i j( ),i j86 applied on both linear arrays. In DAS, two PA images are reconstructed for each linear array and then summed up. Figure 4.13 shows the reconstructed images with the L-shape detection using DAS and VR-SGD. The phantom is a PVC tube filled with methylene blue solution in ethanol as shown in Figure 4.13(a). The two transducers are placed on the top and right-hand side of the object. The reconstructed images by DAS for a single probe are shown in Figure 4.13(b) and (c), respectively. Each image only captures part of the structure which is parallel to the orientation of the transducers due to the limited view. The calibrated transformation parameters for the two probes are = , , and . Figure 4.13(d) shows the combined image which is obtained by transforming the image reconstructed by the second probe to the coordinate defined by the first probe. Although a more complete structure is reconstructed, the magnitude of the background noise is doubled. The image reconstructed by VR-SGD is shown in Figure 4.13(e), which has a higher contrast and clearer background compared with the result obtained by DAS. The side wall of the tube can be distinguished. Comparing DAS with VR-SGD, artifacts caused by the limited-view of the linear array transducer still exist in both reconstructed images, but VR-SGD algorithm reduces the intensity of the artifact. Moreover, the noise background is suppressed significantly within the same dynamic range. Compared to a single probe configuration, the two probe configuration shows more complete structure of the object. q 108.7!50.19 mmh = 7 mmk =87 Figure 4.13 PAT images of a tube phantom detected by two linear array transducers. (a) Photograph of the phantom. (b) Image reconstructed by the first single linear array transducer which is placed on the top of the image. (c) Image reconstructed by the second single linear array transducer which is placed on the right-hand side of the image. (d) Combined image which overlaps the results from the first and second probe. (e) Image reconstructed by VR-SGD using the data from the two probes. Comparing to simulation results, the performance of VR-SGD on experimental data shows improvement but less than the simulation case. In the simulation shown in the previous section, the projection matrix is ideal which accurately describes the relationship between the optical 88 absorption map and the measured data . Thus, the VR-SGD performs very well and can reconstruct the images even in the linear array detection. In experiment, the absolute value of the Hilbert transform of the measured PA signal is taken as . Considering the high frequency fluctuations in the original PA data, a low-pass filter may be applied to improve the convergence of the cost function. Other factors such as the frequency response and the spatial temporal response of the transducers are not included in the forward projection matrix. Thus, the current forward projection matrix may not fully describe the relationship between the optical absorption map and the measured data accurately. Therefore, the experimental result obtained by VR-SGD does not provide as significant improvement as the simulation results. To improve the experimental results, more accurate forward projection matrix will be needed. In the current study, VR-SGD is applied and validated on 2D imaging. When it comes to sources which excites PA signal in 3D space, a corresponding 3D projection matrix will be needed. In the future, 3D PAT imaging of prostate is possible using two probes, a TRUS probe and a pick-up probe. The pick-up probe can be controlled and scanned by a da Vinci robot. The VR-SGD algorithm can be modified and applied for 3D image reconstruction. 4.5 Summary A VR-SGD algorithm for PAT image reconstruction is proposed and verified by 2D simulation and experiment. In the simulation, four different detection configurations are studied by imaging the tree branch and FORBUILD head phantom. The reconstructed results show that the RMSE can be gradually reduced by each time of iteration. Among the proposed detection geometries, images reconstructed by the L-shape detection is comparable to the results obtained by the square-shape detection. L-shape detection is also more achievable in clinical applications. In experiment, VR-SGD algorithm demonstrates its capability to reduce the background noise and artifacts generated ( , )g tr( , )g tr89 by DAS. Overall, the VR-SGD algorithm is shown to perform better than DAS for PAT image reconstruction and using two linear probes can also improve the performance than a single probe. 90 Chapter 5: Conclusions and future work 5.1 Overview PAT is an emerging technology which overcomes the optical diffusion limit based on the PA effect to realize high resolution imaging with deep penetration. PAT has shown great potential to distinguish malignant and benign tissue in human prostate by detecting tumor angiogenesis. In this research, the transurethral illumination for prostate imaging and PAT reconstruction were investigated. The goal was to study the constraints and conditions of performing transurethral illumination. In addition, an iterative reconstruction algorithm was developed when the detection view is limited as in the case for prostate imaging. 5.2 Summary In PAT for prostate imaging, the application is still limited by several challenges, which includes insufficient local fluence for deep penetration, optimizing a compact laser source for transrectal or transurethral illumination, managing the high energy nanosecond laser pulse delivery through multimode fiber, and improving the imaging sensitivity with an advanced reconstruction algorithm. We have been making advancements for translating PAT into clinical prostate imaging in those areas. In Chapter 2, we demonstrated a fiber coupling scheme with a beam homogenizer for coupling high energy pulses into a single multimode fiber. This scheme can benefit prostate imaging which requires miniaturized illumination or internal illumination with a small fiber. A theoretical model was built based on Fourier optics to study the diffraction of the cross cylindrical lens array at the focal plane of a focusing lens. The diffraction pattern was derived in theory and analyzed based on simulation and experiment. The theoretical and experimental analysis reveals the principle where the cross cylindrical lens array distributes the focused pulse energy over a 2D 91 array of beam spots and thus lowers the power density on the fiber tip surface. The scheme increased the pulse energy that was coupled into a multimode fiber and reduced the fiber damage. In experiments, coupled energy as high as 48 mJ/pulse and 64 mJ/pulse has been achieved and the corresponding coupling efficiency was 70% and 90% in a 1000-μm and a 1500-μm-core-diameter fiber, respectively. The high energy pulses delivered by the multimode fiber were further tested for PAT imaging in phantoms. Ex vivo and in vivo PAT imaging showed the capability of the system in imaging biological sample with fiber delivery. Our design not only improves the output energy which provides potential for prostate imaging with large area and deep penetration depth, but also enables the miniaturization of the transurethral illumination due to the small core size of the fiber we applied. In Chapter 3, we introduced the design of a transurethral illumination probe that addressed the limitations such as insufficient local fluence for deep tissue penetration and the risk of over-irradiation at the laser-tissue contact surface. Surface equations of the parabolic and circular cylindrical mirror were derived to estimate the illumination pattern in theory and simulation. A near parallel side illumination was achieved by using the parabolic cylindrical mirror in both simulation and experiment. With a 2 cm diffusing fiber, the fiber assembly was inserted into a 25 Fr cystoscope sheath. Due to the unique design of the illumination probe, a 4 cm2 laser-tissue contact surface was provided and the fluence was controlled to be below the ANSI limit. Moreover, based on the light attenuation model, a compensation was added to the reconstruction to obtain deeper penetration. By imaging a prostate mimicking phantom, a penetration depth of 3.5 cm at 10 mJ/cm2 fluence at 700 nm wavelength was demonstrated. An approach of acquiring PAT image from both sides of the probe was also demonstrated by rotating the illumination probe at 0° and 180° for a forward and a backward illumination, respectively. The deep penetration ability makes 92 PAT of the prostate with transurethral illumination feasible. This technique can potentially be used for PCa detection and image-guided targeted biopsy. In Chapter 4, we demonstrated the variance reduced stochastic gradient descent (VR-SGD) algorithm which was first time applied for PAT image reconstruction. Based on the wave equation of PA signal generation, a 2D forward projection matrix was constructed including the directivity pattern of the transducer. The feasibility of using the integrated measured data for reconstruction was studied through simulation and experiment. The VR-SGD algorithm was shown to enhance the contrast and reduce noise and artifacts when compared to the result obtained by DAS for the linear array detection. Considering the limited view problem caused by linear array detection, the VR-SGD algorithm was also applied to L-shape detection by integrating the transformation matrix information into the forward projection matrix. With L-shape detection, the reconstructed image captured the full structure with improved SNR and reduced artifacts. It demonstrated the feasibility and potential for future clinical application. Overall, the VR-SGD algorithm for PAT image reconstruction is shown to be a powerful reconstruction method and significant improvement can be achieved. Overall, we have carried out thorough investigations on the development and implementation of transurethral illumination PAT for prostate imaging. We have made contributions in the following aspects: (1) developed a fiber coupling scheme for delivering high energy laser pulses in PAT imaging; (2) designed a transurethral illumination probe for PAT imaging of prostate; (3) optimized PAT image reconstruction by developing a variance reduced stochastic gradient descent iterative algorithm. Those improvements are important for translating PAT technique to image prostate cancer in the near future. 93 5.3 Future work The current PAT system still has limitations which can be improved in the future. Several developments and investigations can be done to accelerate the translation of PAT to prostate imaging in the future. 5.3.1 Improvement on the transurethral probe In the hardware development, coupling more energy to the illumination fiber will increase the laser fluence and thus will further increase the image penetration depth and SNR. The image penetration depth is mainly limited by the light attenuation. As the MPE from the laser safety standard allows fluence of 100 mJ/cm2 at 1064 nm wavelength [77], the coupling energy can be further enhanced. Although a larger core diameter fiber can handle higher optical power, the fiber size and rigidity need to be considered for easy insertion into the urethra. When coupling high energy pulses into a small fiber, air break down is one of the main reasons which cause the fiber damage during coupling. In order to minimize this effect, a chamber filled with incombustible solution or dry noble gas such as nitrogen can be integrated with the beam homogenizer as a module. This will also bring much convenience for alignment by pre-aligning the optics inside the module and using standard fiber connector. The size of the current transurethral probe is relatively large. Rigid cystoscope of 21-22 Fr is most frequently used for standard diagnostic cystoscopy, while our current design is packaged in a 25 Fr rigid cystoscope sheath. Therefore, reducing the probe size is needed in the future. We can reduce the size of the parabolic cylindrical mirror. Currently only the parabolic mirror is manufactured by 3D printing technology. In the future, the fiber holder described in Chapter 3 can be fused with the mirror and manufactured by 3D printing as an entity with the same material. The fiber holder will provide support and protection for the fiber and thus the cystoscope sheath could 94 be eliminated. The reduction of the mirror size will decrease the beam size and result in higher fluence. To keep the fluence below the MPE requirement, the length of the diffusing fiber can be increased to illuminate over a longer beam pattern. As PCa is commonly seen in elderly men who have enlarged prostates, extending the diffusing fiber length will also allow the illumination of large prostates. Moreover, our current probe design focuses on the illumination in 2D plane. For volumetric imaging of the prostate, much longer image acquisition time is needed because multiple flashes is required to illuminate the prostate by scanning the illumination probe. As the fluence provided by the transurethral illumination probe is sufficient to penetrate the entire prostate, the reflection mirror can be further engineered to reshape the beam which enhances the beam divergence and increases the detectable angle within a single pulse. This improvement will make the detection more time efficient. 5.3.2 Improvement on the reconstruction algorithm For the VR-SGD algorithm applied in the PAT image reconstruction, the algorithm works very well in the simulation data where the forward project matrix is idea and accurate. However, for the experimental data, the improvement by VR-SGD is reduced, likely because the projection matrix in the real situation differs from the simple model. Thus, the current forward projection matrix may not describe the relationship between the optical absorption map and the measured data accurately. From the aspect of transducer receiving, the directivity pattern and the frequency response of the transducers can affect the generation of the projection matrix. From the aspect of propagating, the inhomogeneity of tissue for either light propagation or sound propagation can lead to a more complicated wave propagation model than the one we currently used. Prior information of transducers should be measured and contributed to the projection matrix including 95 the directivity index and the frequency response. The directivity index reflects the relationship between the received signal intensity of a transducer and the direction of the signal. A weighting factor can be added to Eq. (4.5) based on the angle between each pixel on the image plane and each transducer element's normal line. Since our projection matrix is built in spatial domain, it is complicated to update the matrix with the measured frequency response directly. However, the measured PA data can be processed through deconvolution based on the frequency response. Based on the analysis provided in Chapter4, a low-pass filter may also be needed to filter out the unnecessary high frequency component of the original PA data. Using the processed PA data will give VR-SGD better convergence of the cost function. Meanwhile, the illumination pattern determines how light attenuates as depth increases, resulting in the intensity change of PA signal. Light compensation model can be also added into the projection matrix. With the known beam distribution, the excited PA intensity is proportional to the light intensity. Thus, a light intensity distribution matrix on the image plane can be applied to reflect the change of PA intensity in Eq. (4.5). As the model will become more complicated, corresponding optimization for both hardware and algorithm is expected to accelerate the computational speed. So far, our simulation and experiment only focus on 2D images because our projection matrix is built on the ideal 2D propagation model. A 3D model is required for more accurate reconstruction which can deal with more general reconstruction. Moreover, a TRUS probe and a pick-up probe can be applied to 3D PAT imaging of prostate. The VR-SGD algorithm can be modified for 3D image reconstruction in the future. 96 5.3.3 In vivo prostate imaging The transurethral probe will provide an efficient and uniform light illumination from the center of the prostate. Further investigations will be carried out to study the laser safety for patient use. Once the safe laser power level is identified, PAT application to in vivo imaging of prostate can be carried out. PAT imaging will augment ultrasound imaging to provide complementary contrast from the neurovascular bundles and microvasculature surrounding tumors. For receiving the PA signal, a TRUS probe can receive the signal from the posterior side. To address the limited view problem from a single linear array transducer, a second linear array transducer can receive the PA signal from the anterior side. This can be achieved intra-operatively using a da Vinci robot and a pick-up transducer designed for abdominal ultrasound. The pick-up transducer can be controlled and scanned at different positions by the da Vinci robot to acquire PA signals from a wide range of angles. With both the TRUS and abdominal ultrasound transducer arrays, the view angle of the PAT imaging can be greatly increased. Together with our iterative image reconstruction algorithm, the image quality can be improved significantly. The optimum configuration by using the TRUS probe and abdominal ultrasound has been investigated [112, 122]. Other than 2D imaging, 3D volumetric imaging can potentially be achieved by rotating the TRUS and scanning the pick-up probe. There are several potential applications using the integrated US/PAT system. For real-time image guided biopsy, PAT image will improve PCa visibility and localization and thus enhance the accuracy of biopsies. During surgery, PAT imaging can help detect the location and extension of PCa and neurovascular bundles to assistant the surgeons. The real-time PAT imaging will help surgeons to localize PCa and preserve more nerve tissue during a robotic surgery. With a compact 97 transurethral illumination probe and advanced reconstruction algorithm, PAT imaging will be capable of in vivo prostate imaging to detect PCa during robotic surgery or image guide biopsy. 5.4 Conclusion The thesis studied the feasibility of PAT with transurethral illumination for prostate imaging. Currently, our study manages high-energy pulses delivery in a single multimode fiber which enables the miniaturization of transurethral illumination probe. We have shown that PAT with transurethral illumination is capable to penetrate the depth of 3-4 cm that can cover the entire prostate. However, the knowledge of the optical fluence applied to prostate imaging remains unclear which might be the biggest challenge before its clinical translation. So far, the ANSI safety limit for human skin is used as the guidance to estimate the MPE value in PAT for prostate imaging. Since prostate tissue is different than skin, the optimum fluence for illuminating the prostate needs to be studied. A systematic way to measure the fluence should be defined for better evaluation. We have shown the potential of PAT by using two linear array transducers for detection. In practice, the two probes might not be able to be perpendicular to each other. The feasible angle to position two probes depends on the anatomy around the prostate. The addition of the second transducer array collects PA signal from more directions which will benefit the reconstruction. In the future, a specialized PAT system could be built based on the design of transurethral illumination and two probes detection. After validation on animal models and ex vivo human prostate, it can be applied to in vivo prostate imaging. 98 Bibliography 1. A. G. Bell, "Upon the production and reproduction of sound by light." Journal of the Society of Telegraph Engineers 9(34), 404-426 (1880). 2. A. G. Bell, "LXVIII. Upon the production of sound by radiant energy." The London, Edinburgh, Dublin Philosophical Magazine and Journal of Science 11(71), 510-528 (1881). 3. A. Rosencwaig, Photoacoustics and photoacoustic spectroscopy. Wiley (1980). 4. R. A. Kruger, "Photoacoustic ultrasound." Medical Physics 21(1), 127-131 (1994). 5. T. Wilson, Confocal microscopy. Academic press (1990). 6. W. R. Zipfel, R. M. Williams, and W. W. Webb, "Nonlinear magic: multiphoton microscopy in the biosciences." Nature Biotechnology 21(11), 1369 (2003). 7. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, and C. A. Puliafito, "Optical coherence tomography." Science 254(5035), 1178-1181 (1991). 8. L. V. Wang, "Multiscale photoacoustic microscopy and computed tomography." Nature Photonics 3(9), 503 (2009). 9. L. V. Wang and S. Hu, "Photoacoustic tomography: in vivo imaging from organelles to organs." Science 335(6075), 1458-1462 (2012). 10. R. G. Kolkman, P. J. Brands, W. Steenbergen, and T. G. van Leeuwen, "Real-time in vivo photoacoustic and ultrasound imaging." Journal of Biomedical Optics 13(5), 050510 (2008). 99 11. C. Kim, T. N. Erpelding, L. Jankovic, M. D. Pashley, and L. V. Wang, "Deeply penetrating in vivo photoacoustic imaging using a clinical ultrasound array system." Biomedical Optics Express 1(1), 278-284 (2010). 12. K. Maslov, H. F. Zhang, S. Hu, and L. V. Wang, "Optical-resolution photoacoustic microscopy for in vivo imaging of single capillaries." Optics Letters 33(9), 929-931 (2008). 13. H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, "Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging." Nature Biotechnology 24(7), 848 (2006). 14. K. Maslov, G. Stoica, and L. V. Wang, "In vivo dark-field reflection-mode photoacoustic microscopy." Optics Letters 30(6), 625-627 (2005). 15. J. Yao and L. V. Wang, "Photoacoustic microscopy." Laser Photonics Reviews 7(5), 758-778 (2013). 16. S. Oladipupo, S. Hu, J. Kovalski, J. Yao, A. Santeford, R. E. Sohn, R. Shohet, K. Maslov, L. V. Wang, and J. M. Arbeit, "VEGF is essential for hypoxia-inducible factor-mediated neovascularization but dispensable for endothelial sprouting." Proceedings of the National Academy of Sciences 108(32), 13264-13269 (2011). 17. V. Tsytsarev, S. Hu, J. Yao, K. I. Maslov, D. Barbour, and L. V. Wang, "Photoacoustic microscopy of microvascular responses to cortical electrical stimulation." Journal of Biomedical Optics 16(7), 076002 (2011). 18. W. Song, Q. Wei, T. Liu, D. Kuai, H. F. Zhang, J. M. Burke, and S. Jiao, "Integrating photoacoustic ophthalmoscopy with scanning laser ophthalmoscopy, optical coherence 100 tomography, and fluorescein angiography for a multimodal retinal imaging platform." Journal of Biomedical Optics 17(6), 061206 (2012). 19. S. Hu, B. Rao, K. Maslov, and L. V. Wang, "Label-free photoacoustic ophthalmic angiography." Optics Letters 35(1), 1-3 (2010). 20. L. V. Wang and J. Yao, "A practical guide to photoacoustic tomography in the life sciences." Nature Methods 13(8), 627 (2016). 21. Y. Zhou, J. Yao, and L. V. Wang, "Tutorial on photoacoustic tomography." Journal of Biomedical Optics 21(6), 061007 (2016). 22. J. J. Niederhauser, M. Jaeger, R. Lemor, P. Weber, and M. Frenz, "Combined ultrasound and optoacoustic system for real-time high-contrast vascular imaging in vivo." IEEE Transactions on Medical Imaging 24(4), 436-440 (2005). 23. M. A. L. Bell, X. Guo, D. Y. Song, and E. M. Boctor, "Transurethral light delivery for prostate photoacoustic imaging." Journal of Biomedical Optics 20(3), 036002 (2015). 24. A. Horiguchi, M. Shinchi, A. Nakamura, T. Wada, K. Ito, T. Asano, H. Shinmoto, H. Tsuda, and M. Ishihara, "Pilot study of prostate cancer angiogenesis imaging using a photoacoustic imaging system." Urology 108, 212-219 (2017). 25. A. Horiguchi, K. Tsujita, K. Irisawa, T. Kasamatsu, K. Hirota, M. Kawaguchi, M. Shinchi, K. Ito, T. Asano, and H. Shinmoto, "A pilot study of photoacoustic imaging system for improved realâtime visualization of neurovascular bundle during radical prostatectomy." The Prostate 76(3), 307-315 (2016). 26. T. Kitai, M. Torii, T. Sugie, S. Kanao, Y. Mikami, T. Shiina, and M. Toi, "Photoacoustic mammography: initial clinical results." Breast Cancer 21(2), 146-153 (2014). 101 27. M. Xu and L. V. Wang, "Photoacoustic imaging in biomedicine." Review of Scientific Instruments 77(4), 041101 (2006). 28. V. E. Gusev and A. A. Karabutov, "Laser optoacoustics." NASA STI/Recon Technical Report A 93, (1991). 29. F. A. Duck, Physical properties of tissues: a comprehensive reference book. Academic press (2013). 30. X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, "Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain." Nature Biotechnology 21(7), 803 (2003). 31. W. Shu, M. Ai, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, "Broadening the detection view of 2D photoacoustic tomography using two linear array transducers." Optics Express 24(12), 12755-12768 (2016). 32. M. Ai, W. Shu, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, "Design of high energy laser pulse delivery in a multimode fiber for photoacoustic tomography." Optics Express 25(15), 17713-17726 (2017). 33. M. Ai, J.-i. Youn, S. E. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, "Photoacoustic tomography for imaging the prostate: a transurethral illumination probe design and application." Biomedical Optics Express 10(5), 2588-2605 (2019). 34. P. K. Upputuri and M. Pramanik, "Performance characterization of low-cost, high-speed, portable pulsed laser diode photoacoustic tomography (PLD-PAT) system." Biomedical Optics Express 6(10), 4118-4129 (2015). 102 35. T. J. Allen and P. C. Beard, "High power visible light emitting diodes as pulsed excitation sources for biomedical photoacoustics." Biomedical Optics Express 7(4), 1260-1270 (2016). 36. S. A. Ermilov, T. Khamapirad, A. Conjusteau, M. H. Leonard, R. Lacewell, K. Mehta, T. Miller, and A. A. Oraevsky, "Laser optoacoustic imaging system for detection of breast cancer." Journal of Biomedical Optics 14(2), 024007 (2009). 37. X. L. DeánâBen, T. F. Fehm, M. Gostic, and D. Razansky, "Volumetric handâheld optoacoustic angiography as a tool for realâtime screening of dense breast." Journal of Biophotonics 9(3), 253-259 (2016). 38. R. A. Kruger, R. B. Lam, D. R. Reinecke, S. P. Del Rio, and R. P. Doyle, "Photoacoustic angiography of the breast." Medical Physics 37(11), 6096-6100 (2010). 39. J. Tang, J. E. Coleman, X. Dai, and H. Jiang, "Wearable 3-D photoacoustic tomography for functional brain imaging in behaving rats." Scientific Reports 6, 25470 (2016). 40. E. Zhang, J. Laufer, and P. Beard, "Backward-mode multiwavelength photoacoustic scanner using a planar Fabry-Perot polymer film ultrasound sensor for high-resolution three-dimensional imaging of biological tissues." Applied Optics 47(4), 561-577 (2008). 41. M. Xu and L. V. Wang, "Universal back-projection algorithm for photoacoustic computed tomography." Physical Review E 71(1), 016706 (2005). 42. B. E. Treeby and B. T. Cox, "k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields." Journal of Biomedical Optics 15(2), 021314 (2010). 43. M. Heijblom, D. Piras, F. M. van den Engh, M. van der Schaaf, J. M. Klaase, W. Steenbergen, and S. Manohar, "The state of the art in breast imaging using the Twente 103 Photoacoustic Mammoscope: results from 31 measurements on malignancies." European Radiology 26(11), 3874-3887 (2016). 44. P. D. Kumavor, U. Alqasemi, B. Tavakoli, H. Li, Y. Yang, X. Sun, E. Warych, and Q. Zhu, "Coâregistered pulseâecho/photoacoustic transvaginal probe for real time imaging of ovarian tissue." Journal of Biophotonics 6(6â7), 475-484 (2013). 45. H. S. Salehi, T. Wang, P. D. Kumavor, H. Li, and Q. Zhu, "Design of miniaturized illumination for transvaginal co-registered photoacoustic and ultrasound imaging." Biomedical Optics Express 5(9), 3074-3079 (2014). 46. F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, and A. Jemal, "Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries." CA: a cancer journal for clinicians 68(6), 394-424 (2018). 47. F. H. Schröder, J. Hugosson, M. J. Roobol, T. L. Tammela, S. Ciatto, V. Nelen, M. Kwiatkowski, M. Lujan, H. Lilja, and M. Zappa, "Screening and prostate-cancer mortality in a randomized European study." New England journal of medicine 360(13), 1320-1328 (2009). 48. M. L. Palmeri, T. J. Glass, Z. A. Miller, S. J. Rosenzweig, A. Buck, T. J. Polascik, R. T. Gupta, A. F. Brown, J. Madden, and K. R. J. U. i. m. Nightingale, "Identifying clinically significant prostate cancers using 3-D in vivo acoustic radiation force impulse imaging with whole-mount histology validation." Ultrasound in medicine biology 42(6), 1251-1262 (2016). 49. M. Gayet, A. van der Aa, H. P. Beerlage, B. P. Schrier, P. F. Mulders, and H. Wijkstra, "The value of magnetic resonance imaging and ultrasonography (MRI/US)âfusion biopsy 104 platforms in prostate cancer detection: a systematic review." BJU international 117(3), 392-400 (2016). 50. A. Afshar-Oromieh, C. M. Zechmann, A. Malcher, M. Eder, M. Eisenhut, H. G. Linhart, T. Holland-Letz, B. A. Hadaschik, F. L. Giesel, J. J. E. j. o. n. m. Debus, and m. imaging, "Comparison of PET imaging with a 68 Ga-labelled PSMA ligand and 18 F-choline-based PET/CT for the diagnosis of recurrent prostate cancer." 41(1), 11-20 (2014). 51. G. Xu, I. A. Dar, C. Tao, X. Liu, C. X. Deng, and X. Wang, "Photoacoustic spectrum analysis for microstructure characterization in biological tissue: A feasibility study." Applied physics letters 101(22), 221102 (2012). 52. G. Xu, M. C. Davis, J. Siddiqui, S. A. Tomlins, S. Huang, L. P. Kunju, J. T. Wei, and X. Wang, "Quantifying Gleason scores with photoacoustic spectral analysis: feasibility study with human tissues." Biomedical optics express 6(12), 4781-4789 (2015). 53. V. G. Andreev, A. E. Ponomarev, P. M. Henrichs, M. Motamedi, E. Orihuela, E. Eyzaguirre, and A. A. Oraevsky. "Detection of prostate cancer with optoacoustic tomography: feasibility and modeling." Biomedical Optoacoustics IV. (2003) 54. S.-R. Kothapalli, G. A. Sonn, J. W. Choe, A. Nikoozadeh, A. Bhuyan, K. K. Park, P. Cristman, R. Fan, A. Moini, and B. C. Lee, "Simultaneous transrectal ultrasound and photoacoustic human prostate imaging." Science Translational Medicine 11(507), eaav2169 (2019). 55. M. A. L. Bell, N. P. Kuo, D. Y. Song, J. U. Kang, and E. M. Boctor, "In vivo visualization of prostate brachytherapy seeds with photoacoustic imaging." Journal of Biomedical Optics 19(12), 126011 (2014). 105 56. M. A. Yaseen, S. A. Ermilov, H.-P. F. Brecht, R. Su, A. Conjusteau, M. P. Fronheiser, B. A. Bell, M. Motamedi, and A. A. Oraevsky, "Optoacoustic imaging of the prostate: development toward image-guided biopsy." Journal of Biomedical Optics 15(2), 021310 (2010). 57. B. Bungart, Y. Cao, T. Yang-Tran, S. Gorsky, L. Lan, D. Roblyer, M. O. Koch, L. Cheng, T. Masterson, and J.-X. Cheng, "Cylindrical illumination with angular coupling for whole-prostate photoacoustic tomography." Biomedical Optics Express 10(3), 1405-1419 (2019). 58. J.-M. Yang, K. Maslov, H.-C. Yang, Q. Zhou, K. K. Shung, and L. V. Wang, "Photoacoustic endoscopy." Optics Letters 34(10), 1591-1593 (2009). 59. L. Lin, J. Xia, T. T. Wong, L. Li, and L. V. Wang, "In vivo deep brain imaging of rats using oral-cavity illuminated photoacoustic computed tomography." Journal of Biomedical Optics 20(1), 016019 (2015). 60. S. Manohar, A. Kharine, J. C. van Hespen, W. Steenbergen, and T. G. van Leeuwen, "The Twente Photoacoustic Mammoscope: system overview and performance." Physics in Medicine Biology 50(11), 2543 (2005). 61. W. Xia, D. Piras, M. K. Singh, J. C. van Hespen, T. G. van Leeuwen, W. Steenbergen, and S. Manohar, "Design and evaluation of a laboratory prototype system for 3D photoacoustic full breast tomography." Biomedical Optics Express 4(11), 2555-2569 (2013). 62. S. Seidel and G. Phillipps, "Pulse lengthening by intracavity stimulated Brillouin scattering in a Q-switched, phase-conjugated Nd: YAG laser oscillator." Applied Optics 32(36), 7408-7417 (1993). 106 63. B. C. Stuart, M. D. Perry, R. D. Boyd, J. A. Britten, B. W. Shore, M. D. Feit, and A. M. Rubenchik. "Development of high damage threshold optics for petawatt-class short-pulse lasers." Generation, Amplification, and Measurement of Ultrashort Laser Pulses II. (1995) 64. W. Koechner, Solid-state laser engineering. Springer (2013). 65. L. Gallais, J.-Y. Natoli, and C. Amra, "Statistical study of single and multiple pulse laser-induced damage in glasses." Optics Express 10(25), 1465-1474 (2002). 66. R. A. Robinson and I. K. Ilev, "Design and optimization of a flexible high-peak-power laser-to-fiber coupled illumination system used in digital particle image velocimetry." Review of Scientific Instruments 75(11), 4856-4862 (2004). 67. J. Laufer, C. Elwell, D. Delpy, and P. Beard, "Spatially resolved blood oxygenation measurements using time-resolved photoacoustic spectroscopy." in Oxygen Transport to Tissue XXVII. Springer. 155-160 (2006). 68. L. G. Montilla, R. Olafsson, D. R. Bauer, and R. S. Witte, "Real-time photoacoustic and ultrasound imaging: a simple solution for clinical ultrasound systems with linear arrays." Physics in Medicine Biology 58(1), N1 (2012). 69. D. Razansky, A. Buehler, and V. Ntziachristos, "Volumetric real-time multispectral optoacoustic tomography of biomarkers." Nature Protocols 6(8), 1121 (2011). 70. A. Buehler, E. Herzog, D. Razansky, and V. Ntziachristos, "Video rate optoacoustic tomography of mouse kidney perfusion." Optics Letters 35(14), 2475-2477 (2010). 71. H. He, A. Buehler, and V. Ntziachristos, "Optoacoustic endoscopy with curved scanning." Optics Letters 40(20), 4667-4670 (2015). 107 72. T. Schmidt-Uhlig, P. Karlitschek, M. Yoda, Y. Sano, and G. Marowsky, "Laser shock processing with 20 MW laser pulses delivered by optical fibers." The European Physical Journal-Applied Physics 9(3), 235-238 (2000). 73. T. Schmidt-Uhlig, P. Karlitschek, G. Marowsky, and Y. Sano, "New simplified coupling scheme for the delivery of 20 MW Nd: YAG laser pulses by large core optical fibers." Applied Physics B 72(2), 183-186 (2001). 74. B. E. Saleh and M. C. Teich, Fundamentals of photonics. John Wiley & Sons (2019). 75. J. W. Goodman, Introduction to Fourier optics. Roberts and Company Publishers (2005). 76. R. Magnusson and T. Gaylord, "Diffraction efficiencies of thin phase gratings with arbitrary grating shape." Journal of the Optical Society of America 68(6), 806-809 (1978). 77. A. N. S. Institute, American national standard for safe use of lasers. Laser Institute of America (2007). 78. R. L. Siegel, K. D. Miller, and A. Jemal, "Cancer statistics, 2019." CA: a cancer journal for clinicians 69(1), 7-34 (2019). 79. A. Heidenreich, P. J. Bastian, J. Bellmunt, M. Bolla, S. Joniau, T. van der Kwast, M. Mason, V. Matveev, T. Wiegel, and F. Zattoni, "EAU guidelines on prostate cancer. Part 1: screening, diagnosis, and local treatment with curative intentâupdate 2013." European Urology 65(1), 124-137 (2014). 80. F. H. Schröder, A. B. Kruger, J. Rietbergen, R. Kranse, P. v. d. Maas, P. Beemsterboer, and R. Hoedemaeker, "Evaluation of the digital rectal examination as a screening test for prostate cancer." J Journal of the National Cancer Institute 90(23), 1817-1823 (1998). 81. W. J. Catalona, J. P. Richie, F. R. Ahmann, M. L. A. Hudson, P. T. Scardino, R. C. Flanigan, J. B. DeKernion, T. L. Ratliff, L. R. Kavoussi, and B. L. Dalkin, "Comparison 108 of digital rectal examination and serum prostate specific antigen in the early detection of prostate cancer: results of a multicenter clinical trial of 6,630 men." The Journal of Urology 151(5), 1283-1290 (1994). 82. M. NOGUCHI, T. A. STAMEY, J. E. McNEAL, and C. M. YEMOTO, "Relationship between systematic biopsies and histological features of 222 radical prostatectomy specimens: lack of prediction of tumor significance for men with nonpalpable prostate cancer." The Journal of Urology 166(1), 104-110 (2001). 83. B. Djavan, V. Ravery, A. Zlotta, P. Dobronski, M. Dobrovits, M. Fakhari, C. Seitz, M. Susani, A. Borkowski, and L. Boccon-Gibod, "Prospective evaluation of prostate cancer detected on biopsies 1, 2, 3 and 4: when should we stop?" The Journal of Urology 166(5), 1679-1683 (2001). 84. E. D. Nelson, C. B. Slotoroff, L. G. Gomella, and E. J. Halpern, "Targeted biopsy of the prostate: the impact of color Doppler imaging and elastography on prostate cancer detection and Gleason score." Urology 70(6), 1136-1140 (2007). 85. M. Kongnyuy, A. K. George, A. R. Rastinehad, and P. A. Pinto, "Magnetic resonance imaging-ultrasound fusion-guided prostate biopsy: review of technology, techniques, and outcomes." Current Urology Reports 17(4), 32 (2016). 86. C. M. Hoeks, J. O. Barentsz, T. Hambrock, D. Yakar, D. M. Somford, S. W. Heijmink, T. W. Scheenen, P. C. Vos, H. Huisman, and I. M. van Oort, "Prostate cancer: multiparametric MR imaging for detection, localization, and staging." Radiology 261(1), 46-66 (2011). 87. K. M. Pondman, J. J. Fütterer, B. ten Haken, L. J. S. Kool, J. A. Witjes, T. Hambrock, K. J. Macura, and J. O. Barentsz, "MR-guided biopsy of the prostate: an overview of techniques and a systematic review." European Urology 54(3), 517-527 (2008). 109 88. G. Russo, M. Mischi, W. Scheepens, J. J. De la Rosette, and H. Wijkstra, "Angiogenesis in prostate cancer: onset, progression and imaging." BJU international 110(11c), E794-E808 (2012). 89. A. Postema, M. Mischi, J. De La Rosette, and H. Wijkstra, "Multiparametric ultrasound in the detection of prostate cancer: a systematic review." World Journal of Urology 33(11), 1651-1659 (2015). 90. R. Alonzi, A. R. Padhani, and C. Allen, "Dynamic contrast enhanced MRI in prostate cancer." European Journal of Radiology 63(3), 335-350 (2007). 91. G. J. Jager, E. Ruijter, C. A. van de Kaa, J. J. de la Rosette, G. Oosterhof, J. R. Thornbury, and J. Barentsz, "Local staging of prostate cancer with endorectal MR imaging: correlation with histopathology." American Journal of Roentgenology 166(4), 845-852 (1996). 92. A. Ng and J. Swanevelder, "Resolution in ultrasound imaging." Continuing Education in Anaesthesia Critical Care & Pain 11(5), 186-192 (2011). 93. X. Wang, W. W. Roberts, P. L. Carson, D. P. Wood, and J. B. Fowlkes, "Photoacoustic tomography: a potential new tool for prostate cancer." Biomedical Optics Express 1(4), 1117-1126 (2010). 94. J. L. Su, R. R. Bouchard, A. B. Karpiouk, J. D. Hazle, and S. Y. Emelianov, "Photoacoustic imaging of prostate brachytherapy seeds." Biomedical Optics Express 2(8), 2243-2254 (2011). 95. T. Harrison and R. Zemp, "Coregistered photoacoustic-ultrasound imaging applied to brachytherapy." Journal of biomedical optics 16(8), 080502 (2011). 110 96. K. L. Bell, T. Harrison, N. Usmani, and R. J. Zemp. "Integrated transrectal probe for translational ultrasound-photoacoustic imaging." Photons Plus Ultrasound: Imaging and Sensing 2016. (2016) 97. S. Tang, J. Chen, P. Samant, K. Stratton, and L. Xiang, "Transurethral photoacoustic endoscopy for prostate cancer: a simulation study." IEEE Transactions on Medical Imaging 35(7), 1780-1787 (2016). 98. D.-q. Peng, Y.-y. Peng, J. Guo, and H. Li. "Laser illumination modality of photoacoustic imaging technique for prostate cancer." Journal of Physics: Conference Series. (2016) 99. M. Mitterberger, W. Horninger, F. Aigner, G. M. Pinggera, I. Steppan, P. Rehder, and F. Frauscher, "Ultrasound of the prostate." Cancer Imaging 10(1), 40 (2010). 100. G. Collins, R. Lee, G. McKelvie, A. Rogers, and M. Hehir, "Relationship between prostate specific antigen, prostate volume and age in the benign prostate." British journal of urology 71(4), 445-450 (1993). 101. L. Eri, H. Thomassen, B. Brennhovd, and L. Håheim, "Accuracy and repeatability of prostate volume measurements by transrectal ultrasound." Prostate cancer 5(4), 273 (2002). 102. T. H. Nguyen, Y.-h. Rhee, J.-c. Ahn, and H. W. Kang, "Circumferential irradiation for interstitial coagulation of urethral stricture." Optics Express 23(16), 20829-20840 (2015). 103. M. A. Kosoglu, R. L. Hood, J. H. Rossmeisl Jr, D. C. Grant, Y. Xu, J. L. Robertson, M. N. Rylander, and C. G. Rylander, "Fiberoptic microneedles: Novel optical diffusers for interstitial delivery of therapeutic light." Lasers in Surgery and Medicine 43(9), 914-920 (2011). 111 104. S. Bu, Z. Liu, T. Shiina, K. Kondo, M. Yamakawa, K. Fukutani, Y. Someda, and Y. Asao, "Model-based reconstruction integrated with fluence compensation for photoacoustic tomography." IEEE Transactions on Biomedical Engineering 59(5), 1354-1363 (2012). 105. R. A. Kruger, P. Liu, Y. R. Fang, and C. R. Appledorn, "Photoacoustic ultrasound (PAUS)âreconstruction tomography." Medical Physics 22(10), 1605-1609 (1995). 106. W.-F. Cheong, S. A. Prahl, and A. J. Welch, "A review of the optical properties of biological tissues." IEEE Journal of Quantum Electronics 26(12), 2166-2185 (1990). 107. G. Marquez, L. V. Wang, S.-P. Lin, J. A. Schwartz, and S. L. Thomsen, "Anisotropy in the absorption and scattering spectra of chicken breast tissue." Applied Optics 37(4), 798-804 (1998). 108. T. C. Zhu, A. Dimofte, J. C. Finlay, D. Stripp, T. Busch, J. Miles, R. Whittington, S. B. Malkowicz, Z. Tochner, and E. Glatstein, "Optical Properties of Human Prostate at 732 nm Measured In Vivo During Motexafin Lutetiumâmediated Photodynamic Therapy." Photochemistry Photobiology 81(1), 96-105 (2005). 109. R. A. Weersink, A. Bogaards, M. Gertner, S. R. Davidson, K. Zhang, G. Netchev, J. Trachtenberg, and B. C. Wilson, "Techniques for delivery and monitoring of TOOKAD (WST09)-mediated photodynamic therapy of the prostate: clinical experience and practicalities." Journal of Photochemistry Photobiology B: Biology 79(3), 211-222 (2005). 110. T. Svensson, S. Andersson-Engels, M. Einarsdóttír, and K. Svanberg, "In vivo optical characterization of human prostate tissue using near-infrared time-resolved spectroscopy." Journal of Biomedical Optics 12(1), 014022 (2007). 112 111. S. A. Prahl, M. J. van Gemert, and A. J. Welch, "Determining the optical properties of turbid media by using the addingâdoubling method." Applied Optics 32(4), 559-568 (1993). 112. H. Moradi, S. Tang, and S. E. Salcudean, "Toward intra-operative prostate photoacoustic imaging: configuration evaluation and implementation using the da Vinci research kit." IEEE Transactions on Medical Imaging 38(1), 57-68 (2018). 113. D. Georgescu, E. Alexandrescu, R. Mulţescu, and B. Geavlete, "Cystoscopy and urinary bladder anatomy." in Endoscopic Diagnosis and Treatment in Urinary Bladder Pathology. Elsevier. 1-24 (2016). 114. M. Li, B. Lan, W. Liu, J. Xia, and J. Yao, "Internal-illumination photoacoustic computed tomography." Journal of Biomedical Optics 23(3), 030506 (2018). 115. M. Rubin, M. Buyyounouski, E. Bagiella, S. Sharir, A. Neugut, M. Benson, A. De La Taille, A. Katz, C. Olsson, and R. Ennis, "Microvessel density in prostate cancer: lack of correlation with tumor grade, pathologic stage, and clinical outcome." Urology 53(3), 542-547 (1999). 116. L. A. Mucci, A. Powolny, E. Giovannucci, Z. Liao, S. A. Kenfield, R. Shen, M. J. Stampfer, and S. K. Clinton, "Prospective study of prostate tumor angiogenesis and cancer-specific mortality in the health professionals follow-up study." Journal of Clinical Oncology 27(33), 5627 (2009). 117. S. Sinha, N. A. Rao, B. K. Chinni, and V. S. Dogra, "Evaluation of frequency domain analysis of a multiwavelength photoacoustic signal for differentiating malignant from benign and normal prostates: ex vivo study with human prostates." Journal of Ultrasound in Medicine 35(10), 2165-2177 (2016). 113 118. S. Huang, Y. Qin, Y. Chen, J. Pan, C. Xu, D. Wu, W. Y. Chao, J. T. Wei, S. A. Tomlins, and X. Wang, "Interstitial assessment of aggressive prostate cancer by physioâchemical photoacoustics: An ex vivo study with intact human prostates." Medical Physics 45(9), 4125-4132 (2018). 119. J. Xia, C. Huang, K. Maslov, M. A. Anastasio, and L. V. Wang, "Enhancement of photoacoustic tomography by ultrasonic computed tomography based on optical excitation of elements of a full-ring transducer array." Optics Letters 38(16), 3140-3143 (2013). 120. D. Yang, D. Xing, S. Yang, and L. Xiang, "Fast full-view photoacoustic imaging by combined scanning with a linear transducer array." Optics Express 15(23), 15566-15575 (2007). 121. G. Li, L. Li, L. Zhu, J. Xia, and L. V. Wang, "Multiview Hilbert transformation for full-view photoacoustic computed tomography using a linear array." Journal of Biomedical Optics 20(6), 066010 (2015). 122. H. Moradi, S. Tang, and S. E. Salcudean, "Toward Robot-Assisted Photoacoustic Imaging: Implementation Using the da Vinci Research Kit and Virtual Fixtures." IEEE Robotics Automation Letters 4(2), 1807-1814 (2019). 123. C. Li and L. V. Wang, "Photoacoustic tomography and sensing in biomedicine." Physics in Medicine Biology 54(19), R59 (2009). 124. Y. Xu and L. V. Wang, "Time reversal and its application to tomography with diffracting sources." Physical Review Letters 92(3), 033902 (2004). 125. P. Burgholzer, C. Hofer, G. Paltauf, M. Haltmeier, and O. Scherzer, "Thermoacoustic tomography with integrating area and line detectors." IEEE Transactions on Ultrasonics, Ferroelectrics,Frequency Control 52(9), 1577-1583 (2005). 114 126. P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf, "Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface." Physical Review E 75(4), 046706 (2007). 127. M. Xu, Y. Xu, and L. V. Wang, "Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries." IEEE Transactions on Biomedical Engineering 50(9), 1086-1099 (2003). 128. Y. Xu, D. Feng, and L. V. Wang, "Exact frequency-domain reconstruction for thermoacoustic tomography. I. Planar geometry." IEEE transactions on Medical Imaging 21(7), 823-828 (2002). 129. R. A. Kruger, D. R. Reinecke, and G. A. Kruger, "Thermoacoustic computed tomographyâtechnical considerations." Medical Physics 26(9), 1832-1837 (1999). 130. Z. Wang, J. Li, and R. Wu, "Time-delay-and time-reversal-based robust capon beamformers for ultrasound imaging." IEEE Transactions on Medical Imaging 24(10), 1308-1322 (2005). 131. S. Antholzer, H. Haltmeier, and S. Johannes, "Deep learning for photoacoustic tomography from sparse data." Inverse problems in science and engineering 27(7) 987-1005(2019). 132. H. Deng, X. Wang, C. Cai, J. Luo, and C. Ma, "Machine-learning enhanced photoacoustic computed tomography in a limited view configuration." Advanced Optical Imaging Technologies (2019). 133. F. Liu, X. Gong, L.V. Wang, J. Guan, L. Song, and J. Meng, "Dictionary learning sparse-sampling reconstruction method for in-vivo 3D photoacoustic computed tomography." Biomedical Optics Express 10(4), 1660-1677 (2019). 115 134. J. Provost, and L. Frédéric, "The application of compressed sensing for photo-acoustic tomography." IEEE transactions on medical imaging 28(4), 585-594 (2008). 135. X. Lin, N. Feng, Y. Qu, D. Chen, Y. Shen, and M. Sun, "Compressed sensing in synthetic aperture photoacoustic tomography based on a linear-array ultrasound transducer." Chinese Optics Letters 15(10), 101102 (2017). 136. Z. Guo, C. Li, L. Song, and L. V. Wang, "Compressed sensing in photoacoustic tomography in vivo." Journal of Biomedical Optics 15(2), 021311 (2010). 137. G. Paltauf, J. Viator, S. Prahl, and S. Jacques, "Iterative reconstruction algorithm for optoacoustic imaging." The Journal of the Acoustical Society of America 112(4), 1536-1544 (2002). 138. Y. Zhang, Y. Wang, and C. Zhang, "Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction." Ultrasonics 52(8), 1046-1055 (2012). 139. L. Yao and H. Jiang, "Enhancing finite element-based photoacoustic tomography using total variation minimization." Applied Optics 50(25), 5031-5041 (2011). 140. K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, "Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography." Physics in Medicine Biology 57(17), 5399 (2012). 141. C. Zhang, Y. Zhang, and Y. Wang, "A photoacoustic image reconstruction method using total variation and nonconvex optimization." Biomedical Engineering Online 13(1), 117 (2014). 142. J. Meng, L. V. Wang, D. Liang, and L. Song, "In vivo optical-resolution photoacoustic computed tomography with compressed sensing." Optics Letters 37(22), 4573-4575 (2012). 116 143. L. Bottou, "Stochastic gradient descent tricks." in Neural networks: Tricks of the trade. Springer. 421-436 (2012). 144. D. Karimi and R. K. Ward, "A hybrid stochastic-deterministic gradient descent algorithm for image reconstruction in cone-beam computed tomography." Biomedical Physics Engineering Express 2(1), 015008 (2016). 145. J. Zhang, M. A. Anastasio, P. J. La Rivière, and L. V. Wang, "Effects of different imaging models on least-squares image reconstruction accuracy in photoacoustic tomography." IEEE Transactions on Medical Imaging 28(11), 1781-1790 (2009). 146. P. M. Morse, K. U. Ingard, and R. Beyer, Theoretical acoustics. 1969, American Society of Mechanical Engineers Digital Collection. 147. L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based noise removal algorithms." Physica D: Nonlinear Phenomena 60(1-4), 259-268 (1992). 148. C. R. Vogel, Computational methods for inverse problems. Siam (2002). 149. L. Xiao and T. Zhang, "A proximal stochastic gradient method with progressive variance reduction." SIAM Journal on Optimization 24(4), 2057-2075 (2014). 150. P. L. Combettes and J.-C. Pesquet, "Proximal splitting methods in signal processing." in Fixed-point algorithms for inverse problems in science and engineering. Springer. 185-212 (2011). 151. A. Chambolle, "An algorithm for total variation minimization and applications." Journal of Mathematical Imaging Vision 20(1-2), 89-97 (2004). 152. C. Zhang and Y. Wang, "High total variation-based method for sparse-view photoacoustic reconstruction." Chinese Optics Letters 12(11), 111703 (2014). 117 153. J. Xia, M. R. Chatni, K. I. Maslov, Z. Guo, K. Wang, M. A. Anastasio, and L. V. Wang, "Whole-body ring-shaped confocal photoacoustic computed tomography of small animals in vivo." Journal of Biomedical Optics 17(5), 050506 (2012). 154. J. Gamelin, A. Maurudis, A. Aguirre, F. Huang, P. Guo, L. V. Wang, and Q. Zhu, "A real-time photoacoustic tomography system for small animals." Optics Express 17(13), 10489-10498 (2009). 155. C. Li, A. Aguirre, J. K. Gamelin, A. Maurudis, Q. Zhu, and L. V. Wang, "Real-time photoacoustic tomography of cortical hemodynamics in small animals." Journal of Biomedical Optics 15(1), 010509 (2010). 156. A. A. Goshtasby, Image registration: Principles, tools and methods. Springer Science & Business Media (2012).
Thesis/Dissertation
eng
Biomedical Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Photoacoustic tomography for prostate imaging by transurethral illumination
Text