INVESTIGATING SIGNAL DENOISING AND ITERATIVE RECOSNTRUCTION ALGORITHMS IN PHOTOACOUSTIC TOMOGRAPHY by Jiayi Cheng Master of Engineering, Tsinghua University, 2015 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2017 © Jiayi Cheng, 2017 ii Abstract Photoacoustic tomography (PAT) is a promising biomedical imaging modality that achieves strong optical contrast and high ultrasound resolution. This technique is based on the photoacoustic (PA) effect which refers to illuminating the tissue by a nanosecond pulsed laser and generating acoustic waves by thermoelastic expansion. By detecting the PA waves, the initial pressure distribution that corresponds to the optical absorption map can be obtained by a reconstruction algorithm. In the linear array transducer based data acquisition system, the PA signals are contaminated with various noises. Also, the reconstruction suffers from artifacts and missing structures due to the limited detection view. We aim to reduce the effect of noise by a denoising preprocessing. The PAT system with a linear array transducer and a parallel data acquisition system (DAQ) has prominent band-shaped noise due to signal interference. The band-shaped noise is treated as a low-rank matrix, and the pure PA signal is treated as a sparse matrix, respectively. Robust principal component analysis (RPCA) algorithm is applied to extract the pure PA signal from the noise contaminated PA measurement. The RPCA approach is conducted on experiment data of different samples. The denoising results are compared with several methods and RPCA is shown to outperform the other methods. It is demonstrated that RPCA is promising in reducing the background noise in PA image reconstruction. We also aim to improve the iterative reconstruction. The variance reduced stochastic gradient descent (VR-SGD) algorithm is implemented in PAT reconstruction. A new forward projection matrix is also developed to more accurately match with the measurement data. Using different evaluation criteria, such as peak signal-to-noise ratio (PSNR), relative root-mean-square of reconstruction error (RRMSE) and line profile comparisons, the reconstructions from various iterative algorithms are compared. The advantages of VR-SGD are demonstrated on both simulation and experimental data. Our results indicate that VR-SGD in combination with the accurate projection matrix can lead to improved reconstruction in a small number of iterations. RPCA denoising and VR-SGD iterative reconstruction have been implemented in PAT. Our results show that RPCA and VR-SGD are promising approaches to improve the image reconstruction quality in PAT. iii Lay Summary This thesis will apply the optimization methods to denoise the PA signals and iteratively reconstruct the PA images. As for denoising the PA signals, we use the RPCA algorithm, which can exact the low-rank matrix from an observation matrix, to remove the low-rank band-shaped noise. With respect to the iterative reconstruction, a new projection matrix, which describes the relation between the PA signal and the image, is developed and the VR-SGD algorithm is implemented. To the best of our knowledge, this is the first time that the RPCA and VR-SGD are used in PAT. iv Preface Part of the work done in Chapter 4 will be submitted. J. Cheng, D. Karimi, M. Ai, T. Salcudean, R. Rohling, P. Abolmaesumi, and S. Tang, “Variance-reduced stochastic gradient descent algorithm for photoacoustic image reconstruction,” in preparation. I conducted all the experiments, analyzed the results, and wrote the manuscript. Min Ai was involved in the helping with conducting experiments. Davood Karimi, Tim Salcudean, Robert Rohling, Purang Abolmaesumi, and Shuo Tang were involved throughout the project in concept formation and manuscript composition. v Table of Contents Abstract .......................................................................................................................................... ii Lay Summary ............................................................................................................................... iii Preface ........................................................................................................................................... iv Table of Contents ...........................................................................................................................v List of Tables ............................................................................................................................... vii List of Figures ............................................................................................................................. viii List of Abbreviations ................................................................................................................... xi Acknowledgements .................................................................................................................... xiii Dedication ................................................................................................................................... xiv Chapter 1: Introduction ............................................................................................................... 1 1.1 Photoacoustic imaging ................................................................................................ 2 1.2 Photoacoustic imaging denoising ............................................................................... 3 1.3 Photoacoustic image reconstruction ........................................................................... 5 1.4 Problem statement and motivation.............................................................................. 6 1.5 Organization of the thesis ........................................................................................... 7 Chapter 2: Photoacoustic imaging background .......................................................................... 8 2.1 Photoacoustic imaging principles ............................................................................... 8 2.2 Experimental photoacoustic imaging system............................................................ 11 Chapter 3: Investigating robust principal component analysis for photoacoustic signal denoising ................................................................................................................................... 13 3.1 Introduction ............................................................................................................... 13 vi 3.2 Related work ............................................................................................................. 14 3.3 Methodology ............................................................................................................. 16 3.3.1 RPCA principle ..................................................................................................... 16 3.3.2 RPCA implementation on PA data ....................................................................... 18 3.4 Experiments and analysis .......................................................................................... 19 Chapter 4: Variance-reduced stochastic gradient descent algorithm for photoacoustic image reconstruction ............................................................................................................................ 28 4.1 Introduction ............................................................................................................... 28 4.2 Theoretical background ............................................................................................ 30 4.2.1 Photoacoustic imaging model ............................................................................... 30 4.2.2 Building the basic discrete forward projection matrix .......................................... 31 4.2.3 Gradient method for solving cost function ........................................................... 32 4.3 Methodology ............................................................................................................. 34 4.3.1 A new improved projection matrix ....................................................................... 34 4.3.2 VR-SGD algorithm ............................................................................................... 35 4.4 Simulation validation ................................................................................................ 38 4.5 Experimental implementation ................................................................................... 43 Chapter 5: Conclusion and future work .................................................................................... 47 5.1 Significance of work ................................................................................................. 47 5.2 Future work and improvement .................................................................................. 48 Bibliography .................................................................................................................................50 vii List of Tables Table 2.1: Summary of the PAT system characterization ............................................................ 12 Table 4.1: Main features of different iterative reconstruction algorithms. .................................. 39 Table 4.2: PSNR comparisons by different algorithms after 60 iterations. .................................. 42 viii List of Figures Figure 2.1: Schematic of the PAT experimental system used in this thesis. ND:YAG is a Surelite II neodymium-doped yttrium aluminum garnet laser; OPO is an optical parametric oscillator [67]. .............................................................................. 11 Figure 3.1: Observation matrix D formation process. .................................................................. 18 Figure 3.2: Data acquisition for different cases. (a) Without transducer and laser. (b) With transducer and without laser. (c) With transducer and the laser is on but blocked. ...................................................................................................................... 20 Figure 3.3: Data comparison and power spectrum analyis. (a) Raw PA signal from one channel and the measured band-shaped noise. (b) Frequency spectrum analysis. ...................................................................................................................... 21 Figure 3.4: Data separation results. (a) Foreground data corresponds to the pure PA data. (b) Background data refers to the band-shaped noise. (c) The observation matrix is the original PA signals. ............................................................................... 21 Figure 3.5: DAS reconstructions and line profile comparisons. (a, b, c, d) DAS reconstruction from the raw PA signals, background subtracted data, SVD denoised data and RPCA denoised data, respectively. (e) Horizontal line profile comparison (along the red arrow). (f) Vertical line profile comparison (along the red arrow). ................................................................................................. 22 Figure 3.6: Data separation results for the dots array data. (a) Foreground data corresponds to the pure PA data. (b) Background data refers to the band-shaped noise. (c) The observation matrix is the original PA data.......................................................... 23 ix Figure 3.7: DAS reconstructions. (a, b, c, d) DAS reconstruction from the raw PA data, background subtracted data, SVD denoised data and RPCA denoised data, respectively. ............................................................................................................... 24 Figure 3.8: Line profile comparisons. (a) Horizontal line profile comparison. (b) Vertical line profile comparison. ............................................................................................. 24 Figure 3.9: Data separation results for the human forearm data. (a) Foreground data corresponds to the pure PA signals. (b) Background data refers to the band-shaped noise. (c) Observation matrix is the original PA signals. .............................. 25 Figure 3.10: Background is the ultrasound image and foreground is PA reconstruction. (a, b, c, d) Foreground is the DAS reconstruction from the raw PA data, background subtracted data, SVD denoised data and RPCA denoised data, respectively. ............................................................................................................... 26 Figure 4.1: PAT simulation setup and signals generation comparison. (a) Simulation setup. (b) PA signals comparison. ........................................................................................ 38 Figure 4.2: PAT reconstructions. (a) Normalized DAS reconstruction with a threshold of 0.16. (b, c, d, e, f) Reconstructions by iterative algorithms: IR, A-SGD, M-SGD and VR-SGD, respectively. ............................................................................... 40 Figure 4.3: Comparisons of PAT reconstructions. (a) RRMSE comparison. (b) Line profile comparison. ................................................................................................................ 40 Figure 4.4: PAT reconstructions for the points array sample and the line profile comparisons. (a) Normalized DAS reconstruction with a threshold of 0.16. (b) M-SGD reconstruction. (c) VR-SGD reconstruction. (d) Horizontal line profile comparison. (e) PSNR1 comparison. ........................................................................ 44 x Figure 4.5: Background is the ultrasound image and foreground PA reconstruction. (a) Foreground is a normalized DAS reconstruction with a threshold of 0.16. (b) Foreground is the M-SGD reconstruction after six iterations. (c) Foreground is the VR-SGD reconstruction after six iterations. ........................................................ 46 xi List of Abbreviations PAT Photoacoustic Tomography PA Photoacoustic RPCA Robust Principal Component Analysis VR-SGD Variance Reduced Stochastic Gradient Descent PSNR Peak Signal-to-Noise Ratio RRMSE Relative Root-Mean-Square of Reconstruction Error CT Computed Tomography MRI Magnetic Resonance Imaging CLSM Confocal Laser Scanning Microscopy OCT Optical Coherence Tomography MPM Multi-Photon Microscopy PAM Photoacoustic Microcopy OA Optoacoustic 3D Three-Dimensional 2D Two-Dimensional 1D One-Dimensional FBP Filtered Back-Projection DAS Delay-and-Sum SNR Signal-to-Noise Ratio OPO Optical Parametric Oscillator DAQ Data Acquisition SVD Singular Value Decomposition xii FGD Full Gradient Descent SGD Stochastic Gradient Descent TV Total Variation EIR Electric Impulse Response A-SGD Projection Matrix A based Stochastic Gradient Descent Algorithm M-SGD Projection Matrix M based Stochastic Gradient Descent Algorithm IR Iterative Algorithm xiii Acknowledgements I would like to express my appreciation to my supervisor Dr. Shuo Tang for her continuous guidance and support. I would also like to thank every single member of the Biophotonics group. Thanks to Min Ai and Weihang Shu for their cooperation and brilliant ideas in the photoacoustic project. Thanks to Lin Huang and Xin Zhou for their selfless assistances and support throughout the years. Thanks to Davood Karimi for his valuable suggestions and help to my paper. I would also like to thank Professor Tim Salcudean, Professor Robert Rohling, Professor Purang Abolmaesumi, and every single member of the Robotic Control laboratories for their help and support. Special thanks are owed to my families, whose have supported me throughout my years of education, both morally and financially. Jiayi Cheng University of British Columbia August 2017 xiv Dedication This thesis is dedicated to my family for the unconditional love every step of the way. 1 Chapter 1: Introduction With the development of science and technology, various imaging tools, such as X-ray, X-ray computed tomography (CT), magnetic resonance imaging (MRI), ultrasonography and various optical imaging modalities, have been widely used for diagnosing diseases noninvasively. As for the X-ray and CT, they can provide deep penetration and high contrast between soft tissues and bones of human body, but the patients have to be exposed to potentially harmful radiations. MRI takes the advantages of magnetic fields, radio waves and field gradient to form images with high resolution and strong contrast, which cannot be achieved by X-ray and CT. However, the cost of MRI is higher and the operation time is longer than those of X-ray and CT. As a highly mobile and relatively low-cost imaging modality, ultrasonography can provide high resolution and good contrast. Since ultrasonography relies on the difference in acoustic impedance to gain contrast, it is hard to differentiate soft tissues which have relatively similar acoustic impedance values. Optical imaging modalities, such as confocal laser scanning microscopy (CLSM), optical coherence tomography (OCT), and multi-photon microscopy (MPM), can offer the highest resolution among all the modalities mentioned above by using light with different wavelengths. The downside of optical imaging is its shallow imaging depth due to the significant scattering of light in tissues. Photoacoustic (PA) imaging, as an emerging technique in biomedical imaging, is based on the detection of acoustic waves that are generated by thermoelastic expansion of light irradiated tissues [1]. As a hybrid modality, PA imaging has the advantages of deep penetration of ultrasound as well as the biochemically specific optical absorption contrast. These advantages have greatly increased the applications of PA imaging in recent years [2-5]. The application areas of PA imaging include skin, prostate, breast, liver, brain and nerve [6-11]. 2 1.1 Photoacoustic imaging PA imaging is mainly based on the PA or optoacoustic (OA) effect which was firstly discovered by Alexander Graham Bell in 1880s [14]. Sound waves were generated from a solid sample irradiated by a beam of modulated sunlight. It was also discovered that irradiated liquid and gas materials could generate acoustic waves [15, 16]. With the development of lasers, the light source for exciting PA signals became more stable and reliable. PA signal detection from gas materials and photochemical studies were conducted in the 1970s and 1980s [17, 18]. In 1980s, Bowen applied the PA effect for imaging the biological information of tissues and proved the corresponding concepts [19]. Later, PA imaging has grown rapidly and many PA imaging studies were published [20-25]. PA effect is based on the generation of ultrasound waves by the absorption of laser energy. Illuminated by a short-pulsed laser, an optically absorbing sample, such as hemoglobin, melanin, water, or lipid, experiences a rapid temperature increase and an abrupt and localized pressure change is produced by thermoelastic expansion. This pressure change results in the PA waves that propagate from the heated region out toward the surface and then can be collected by ultrasound transducers. As for the incident light, the pulse width is usually on the order of nanoseconds. Since the scattering of acoustic waves in tissues is two to three orders of magnitude lower than that of optical waves, PA imaging can achieve deeper penetration. Also, the contrast of PA imaging is based on optical absorption which depends on the specific absorbing molecules and the light wavelength. Generally, PA imaging can be divided into two modes: photoacoustic tomography (PAT) and photoacoustic microcopy (PAM) [12, 13]. PAT images a relatively large tissue area and the image is reconstructed from PA signals collected at different angels by an ultrasound transducer 3 array or a scanning transducer. PAM images a smaller area but with higher resolution than PAT by focusing the laser beam or the ultrasound transducer. This thesis mainly focuses on the usage of PAT. With respect to the signal collection in PAT, either a probe with single element or a specially designed probe with multiple channels can be selected. The single detector often works with a scanning device such as a stepping motor, which makes the system more complicated and time consuming. By using a linear array transducer, the detection process is more flexible and faster. However, linear arrays still suffer from limited detection view, where the linear array cannot collect all the PA signals in a 360º angle in a two-dimensional (2D) imaging configuration. In this thesis, we focus on the usage of a linear array transducer that is widely used in clinics and hospitals. From the PA data, the optical absorption distribution can be recovered by reconstruction algorithms. As the PA signals collected by a linear array transducer are insufficient and always corrupted with noise, the reconstruction may contain strong artifacts and noisy background. Thus a signal denoising method and an improved reconstruction method are required. 1.2 Photoacoustic imaging denoising In PA imaging, the acoustic waves are detected by an ultrasound probe, and then amplified and digitized by a sophisticated electronic system. Normally, the data acquisition process is under a noisy environment and the PA signals are corrupted with additive noises, which are caused by different sources such as electronic thermal noise or electromagnetic interference. Some of those noises, such as the electronic thermal noise, have a random distribution that can be modeled as white Gaussian noise. Other types of noise have also been observed in PA signals, such as the fixed pattern noise caused by electromagnetic interference, which shows specific noise characteristics. In Ref. [43], a band-shaped electronic noise named as laser-induced noise is 4 observed, which is produced by the driving circuitry of the laser and coupled through cables due to electromagnetic interference. This noise is prominent for cases of parallel PA data acquisition and is presented as vertical or horizontal bands. For the reduction of random noise, averaging across multiple frames repetitively acquired from the same location is commonly used. However, this method is based on the collection of multiple frames of data, which is time consuming and not suitable for the target-moving cases. Another commonly used denoising method for the random noise is filtering, but this is not effective when the noise has overlapping bandwidth with the PA signals. Since the spectrum of PA signal depends on the size of the sample and the sample may comprise different sizes, the PA signals are usually broadband and the spectra of the PA signal and the noise may overlap [30]. Thus, a frequency filtering method could remove some useful signal components overlapped with noise. The wavelet method is a promising denoising method and decomposes the signals into a series of basis functions with different coefficients. By thresholding the small values of the coefficients, the corresponding white random noise can be removed [66]. For the reduction of the fixed pattern noise, singular value decomposition (SVD) is used in Ref. [43]. The band-shaped noise can be partly removed by zeroing the singular values. However, this method requires curved signal wavefronts in a 2D data matrix so as to be distinguished from the band-shaped noise. Accordingly, when the targets are too far away from or parallel to the transducer, the SVD denoising method performs poorly. Besides, this method is limited by the cases that the magnitude of noise is lower than that of the PA signals. Therefore, a better denoising method is still needed in order to improve the image reconstruction. 5 1.3 Photoacoustic image reconstruction From the acquired PA signals, the optical absorption map can be obtained by reconstruction algorithms [35-42]. Reconstruction methods can be classified into direct and iterative methods. With respect to the direct methods, since the PA process can be modeled by heat conduction and wave equations, it is intuitive to reconstruct the absorption map by analytically inversing the corresponding equations in either time domain or frequency domain [35, 36]. Besides, the filtered back-projection (FBP) algorithm is widely used because of its efficiency and ability to analytically reconstruct PAT images [37]. Delay-and-sum (DAS) beamforming is a widely used reconstruction algorithm in ultrasound imaging, which is also implemented in PAT [38]. In DAS, the PA signal received by each transducer element is back-projected as an arc based on the wave propagation time. The reconstruction is the superposition of the arcs from all the transducer elements. When the transducer array cannot receive the PA signals from 360º angles completely, significant artifacts can be generated. Even though the direct methods are known for their ease of implementation and low computation cost, their reconstructions are often affected by the insufficient measurements and noise. Iterative methods, on the other hand, start with an initial image estimate, which can be obtained using a direct reconstruction method, and gradually improve the image estimate [39-42]. Iterative methods can improve the image reconstruction by minimizing a cost function which consists of the difference between the image estimate and the measurement data, and a prior information of the PA image. More details about iterative methods will be presented in Chapter 4. 6 1.4 Problem statement and motivation Linear array transducers are widely used in clinics and hospitals. However, it has incomplete detection view in detecting PA signals. The direct reconstructions suffer from artifacts due to the incomplete measurements. PA signals measured by a linear array transducer are also contaminated with band-shaped noise which leads to a noisy background of the reconstruction. This thesis aims to reduce the band-shaped noise and mitigate the artifacts in reconstruction through optimization methods. With respect to the removal of band-shaped noise, a background subtraction method can be easily implemented. However, the background cannot be exactly the same in two measurements and a simple subtraction may also remove some true signal. Another approach to remove the band-shaped noise is SVD. This method requires the magnitude of band-shaped noise is comparable to the largest PA signals from sample, otherwise the quality of the reconstruction could be affected. To address this issue, robust principal component analysis (RPCA) method will be applied to improve the reduction of band-shaped noise. In order to mitigate the artifacts caused by insufficient measurements, iterative reconstruction methods will be taken advantage. Iterative methods are based on a discrete forward projection matrix and a gradient descent updating algorithm. Others have reported a projection matrix but it was not paired with the measurements [40]. Thus a large error was introduced. In this thesis, I will propose a new projection matrix which is well paired with the PA data to achieve a better reconstruction. Full gradient descent updating is very computationally costly. I will implement an improved updating algorithm named variance reduced stochastic gradient descent (VR-SGD) to PA imaging to reduce the computational cost and achieve a high quality reconstruction. 7 1.5 Organization of the thesis Chapter 1: An introduction about PA imaging basics, PA imaging denoising methods and PA imaging reconstruction methods. Also, our motivation which aims to remove the band-shaped noise and artifacts is presented. Chapter 2: An overview of PA imaging background. The PA imaging principles are described in details and the experimental PA imaging system is introduced. Chapter 3: We validate the existence of band-shaped noise for the experimental PA data. We use the RPCA to remove the band-shaped noise and compare the results with other denoising methods. Chapter 4: We present an iterative reconstruction method based on a new forward propagation matrix and the gradient descent based algorithm named VR-SGD. Using different evaluation criteria, the reconstructions from various algorithms are compared for both simulation and experimental data and the advantages of VR-SGD are demonstrated. Chapter 5: A summary of all the work done in this thesis and discussions about the outlook of our project. 8 Chapter 2: Photoacoustic imaging background In this chapter, we will introduce the photoacoustic imaging principles and our experimental PAT system. 2.1 Photoacoustic imaging principles When a sample is illuminated by a short-pulsed laser, its temperature will suddenly increase and an abrupt change of pressure will be caused by thermoelastic expansion. This pressure change results in the PA waves that propagate from the heated region to the medium and then can be collected by an ultrasound probe. As for the light source, the pulse width is usually on the order of nanoseconds. In order to successfully generate the PA signals, two physical conditions, namely the thermal and stress confinements, are required to be met. The thermal confinement limits the temporal duration p of the laser pulse. Since heat dissipation can limit the local temperature increase and inhibit the PA effect, the duration of the laser pulse should be shorter than the heat dissipation th . Thus the heat diffusion can be neglected. The duration th can be calculated as [26]: 2~4pthTLD (2.1) where Lp is the characteristic linear dimension of the excited structure, and TD is the thermal diffusivity of the sample. Normally, the dissipation time is several hundred microseconds. The stress confinement considers the stress relaxation time s and can be given as: 9 pssLv (2.2) where sv is the velocity of sound and the typical value in soft tissue is 1540 m/s. The stress relaxation time is of the order of tens of nanoseconds. Therefore, if a laser pulse width is on the order of nanoseconds, both the thermal and stress confinement conditions can be met. Suppose an optical absorber is embedded in a homogeneous medium and a nanosecond short pulsed laser uniformly illuminates the absorber. Part of the light energy is converted to heat. The initial pressure 0 0 0( , )p r t at position 0r and time 0t generated by the laser illumination can be expressed as [26]: 0 0 02( , ) sa apvFCr t F (2.3) where a is the tissue absorption coefficient, is the coefficient of volume thermal expansion , C is the specific heat capacity and F is the local optical fluence. The terms in the parenthesis make up a dimensionless parameter named the Grüneisen parameter that relates the pressure change 0 0 0( , )p r t to the local energy deposition a F . For example, is ~0.80 in fatty tissue and ~0.20 in blood [28]. The initial pressure from the surface of the heated absorber then propagates in the medium, assumed to be homogeneous. The pressure waves ( , )p r t at position r and time t obeys a wave equation [20, 29-33]: 222 22( , ) ( , )( , )ssvp r t S r tv p r tt C t (2.4) 10 where r is the three-dimensional (3D) position vector, t is the time and ( , )S r t is the heat source. The heat source term ( , )S r t , which describes heat-producing and optical absorption within the medium, can be separated into a purely spatial part and a purely temporal part [20]: 0( , ) ( ) ( )S r t I x r G t (2.5) where 0I is a scaling factor proportional to the incident radiation intensity, ( )G t describes the shape of the irradiating pulse and ( )x r describes the optical absorption properties of the medium. By using a Green’s function, the pressure waves can be solved from Eq. 2.4 and written as: 0( , ) ( )4sr r v tI d drp r t x rC dt t (2.6) where r is the 3D position of the source which generates the pressure waves, and is the optical pulse width. According to Eq. 2.6 the pressure waves can be calculated from the heat source in a 3D domain. If the irradiated source is restricted in a thin slice, then Eq. 2.6 can be expressed in a 2D domain. We can get an intermediate variable ( , )g r t which represents a scaled version of the acoustic velocity potential [34]: 004( , ) ( , ) ( )str r v tCtg r t p r t dt x r drI z (2.7) where z is the thickness of the ‘slice’ which can be treated as a constant. Here r and r are now 2D position vectors. The optical absorption ( )x r at position r can be reconstructed from the acoustic pressure ( , )p r t or the acoustic velocity potential ( , )g r t correspondingly. In PAT imaging, ( , )p r t is the PA signal measured by ultrasound transducer and ( )x r is the optical absorption distribution that we want to reconstruct. 11 2.2 Experimental photoacoustic imaging system Figure 2.1: Schematic of the PAT experimental system used in this thesis. ND:YAG is a Surelite II neodymium-doped yttrium aluminum garnet laser; OPO is an optical parametric oscillator [67]. A previously developed PAT system is used to collect all the data used in this thesis [67]. Fig. 2.1 shows the schematic of the PAT experimental system. The system includes a laser source, an optical parametric oscillator (OPO), a data acquisition (DAQ) module, an ultrasound imaging system, and an ultrasound linear array transducer. The laser is a Q-switched Nd:YAG laser (Surelite II, Continuum, Inc., Santa Clara, CA, USA). 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. The ultrasound linear array transducer has 128 channels with 7.2 MHz center frequency, minimum 70% fractional bandwidth (at −6dB), and 0.3 mm element pitch. The probe is connected to a SonixMDP ultrasound imaging system (Ultrasonix Medical Corporation, Richmond, BC, Canada) through a DAQ module SonixDAQ (Ultrasonix Medical Corporation, Richmond, BC, Canada). The DAQ module is responsible for acquiring and digitizing the pre-beamformed radio-frequency signal from all the channels individually at a sampling rate of 40 MHz and 12-bit resolution. The 12 system can acquire both PAT and traditional ultrasound imaging. Table 2.1 describes the system properties. Table 2.1: Summary of the PAT system characterization System Characterization Frame rate 10 Hz Field of view ~3 cm2 Peak laser energy output 120 mJ Tunable laser wavelengths 680 ~ 2500 nm Axial resolution (3 cm away from transducer) 0.44 mm Lateral resolution (3 cm away from transducer) 0.52 mm To acquire a PAT image, a sample is placed in a water tank or embedded in gelatin phantom for acoustic coupling. Laser light shines on the sample and PA signal is collected by the linear array transducer. For each laser pulse, one frame of the PA signals from the 128 channels of the linear array transducer is collected. Multiple frames can be collected by continuously firing the laser pulses at 10 Hz repetition rate and acquiring the data. From the collected PA signals, PAT images will be reconstructed by DAS or iterative algorithms after denoising. 13 Chapter 3: Investigating robust principal component analysis for photoacoustic signal denoising In this chapter, we will apply the RPCA to remove the band-shaped electronic noise. This noise normally exists in the PA signals collected by a linear array transducer and is presented as vertical and horizontal bands. We first validate the existence of this noise by experiments. Then we implement the RPCA to remove it for different experimental samples and compare the results with other denoising methods. 3.1 Introduction PA signals are corrupted by different types of noises, including random noise and fixed pattern noise. The removal of random noise can be achieved by averaging the repetitive data acquisition, filtering in frequency domain or utilizing wavelet methods. In this chapter, we focus on the reduction of the fixed pattern noise which exists in the PA data collection by a linear array transducer. This noise is caused by radiofrequency (RF) coupling of high frequency electronic signal from the surrounding area to the data acquisition electronics of the PA system. It is often presented as vertical and horizontal bands in a 2D data matrix. Generally, PA signals have relatively curved wavefronts in a 2D data matrix due to ultrasound wave propagation, as long as the sample is not too far away. Thus, by decomposing the raw data matrix through SVD, the different weights of the corresponding singular values can be obtained. Specifically, the first a few singular values can well model the band-shaped noise, but not the PA signals. Therefore, the band-shaped noise can be removed by zeroing the first a few singular values of the measurement matrix when the magnitude of the noise is comparable to the largest PA signals [43]. 14 As a matter of fact, the reason why the band-shaped noise can be removed by zeroing the first a few singular values of the raw data matrix is because the first a few singular values contain most information of the band-shaped noise when the magnitude of band-shaped noise is comparable to the largest PA signals. Besides, since there are many repeated columns and rows in the band-shaped noise matrix, the noise matrix can be regarded as a low-rank matrix. On the other hand, the pure PA signals usually occupy a small region of the data matrix and the data matrix can be considered as sparse. Therefore, the raw PA matrix can be considered as the summation of a low-rank matrix and a sparse matrix that refer to the band-shaped noise and pure PA signal, respectively. Accordingly, in order to remove the band-shaped noise and acquire the pure PA signals, we can take the advantages of the RPCA approach [45]. The basic idea of RPCA is to extract either a low-rank matrix or a sparse matrix from the observation data matrix. There are various applications based on RPCA. This method is popularly used for background subtraction in video surveillance systems, which aims to find the moving objects in a scene. Another application of RPCA is to perform denoising of seismic signals by applying this algorithm in the frequency domain, thereby both the Gaussian and erratic noise can be removed [46]. Due to RPCA's outstanding performance on data extraction, we propose to apply RPCA to denoise the PA signals. Specifically, we apply this method to experimental PA data and compare with other methods. 3.2 Related work As for the removal of the band-shaped noise, the averaging and filtering methods are not appropriate, because the spectrum of this band-shaped noise is broad. In order to remove it, a straightforward way is to subtract the band-shaped noise which is measured in advance. However, 15 the band-shaped noise cannot be exactly the same in two measurements and the subtraction may induce additional noise and degrade the signals. We can also realize the noise reduction by SVD [43]. We assume that our observation matrix D is made up of the band-shaped noise matrix B and the pure PA signal matrix E , denoted as D B E . The dimension of these three matrices is i-by-j, where i is the number of samplings of each transducer element and j is the number of transducer elements. By applying SVD, D is decomposed into three matrices: TDD US V (3.1) where DS is a diagonal matrix made up of singular values and its size is i-by-j, U and V are feature square matrices with a size of i-by-i and j-by-j, respectively. Suppose the matrix B and matrix E can also be decomposed with respect to U and V , and we have: ( ) T T T TD B E B EUS V US V US V U S S V (3.2) where BS and ES are the diagonal matrices corresponding to B and E , respectively, with size i-by-j. Hence, the operation of the matrices can be achieved by manipulating the singular values and then performing inverse SVD. Since the number of singular values is equal to the rank of the matrix and the rank of B is often lower than that of E, zeroing the first a few singular values of D can remove the band-shaped noise and obtain a new diagonal matrix DS . Then the denoised PA data matrix D can be obtained as: TDD U S V (3.3) 16 There are limitations of this denoising method. Zeroing the first a few singular values of D could also remove some signals especially when the magnitude of the signals is larger than the noise. Therefore, it lacks robustness. RPCA can also be used to perform denoising. Assuming that the matrix B and E can be treated as a low-rank and a sparse matrix, respectively, in order to extract the signal E (or remove the noise B) from the observation matrix D, we can formulate it into a constrained optimization problem [50]: * 1, arg m in , . . B B EB E s t D E (3.4) where * denotes the nuclear norm (the sum of singular values of the matrix). Minimizing the nuclear norm of B encourages it to be low rank. 1 is the entry-wise L1-norm and this encourages it to be sparse. The entry-wise norm applies after vectorization of the matrix. is a positive parameter. In this optimization problem, the low rank of B and the sparsity of E are optimized together and it is more robust than the SVD denoising method. By selecting a suitable parameter , which is related to the size of the data matrix, the foreground data E can be correctly extracted from D. RPCA has not been applied to remove the band-shaped noise of PA signals before. By exploiting the advantages of RPCA, we propose a new application of RPCA to denoise PA signals. 3.3 Methodology 3.3.1 RPCA principle In order to apply RPCA, the constrained optimization Eq. 3.4 needs to be solved. Currently, there are two main approaches to tackle this optimization problem. We can either use a first-order method to solve the primal problem or solve it through a dual method [50]. Here we choose the 17 first approach and approximate Eq. 3.4 into an unconstrained optimization by imposing a penalty : 2* 1 2, 1arg min 2 B EB E D B E (3.5) where 2 is the entry-wise L2-norm which measures the difference between D and B+E. Here and are positive parameters corresponding to the weights of the L1-norm of E and the difference, respectively. This problem can be treated as a regularized block multi-convex optimization and block coordinate optimization can be applied [51]. Therefore, Eq. 3.5 can be solved in two steps: 21* 21arg min2 t tBB B D B E (3.6) 21 11 21arg min2 t tEE E D B E (3.7) Eq. 3.6 can be solved by Proximal Gradient and SVD method [52]. Eq. 3.7 can be solved by the Proximal Gradient named soft-threshold [53]. 18 3.3.2 RPCA implementation on PA data Figure 3.1: Observation matrix D formation process. In order to implement RPCA on PA signals, the specific features of PA signals need to be considered compared to conventional RPCA applications. With regard to detecting the moving objects in surveillance systems, in order to fit the nature of RPCA, multiple frames of the video are rearranged into a new observation matrix, where each frame is vectorized by unfolding the data into one column of the observation matrix. Since the background of each frame barely changes and the moving object only takes a relatively small part of the whole image, the observation matrix can be assumed to be made up of a low-rank matrix and a sparse matrix. However, the vectorization destroys certain structures of the input data set and increases the computation burden heavily. For the PA signals, we can use a single B-scan frame of the PA signals as the observation matrix D 19 and it is already made up of a low-rank matrix (band-shaped noise) and a sparse matrix (PA signals). Our implementation avoids the vectorization steps and reduces the computational cost. It also shortens the PA data acquisition process without the need to acquire multiple frames. 3.4 Experiments and analysis The experimental system is described in Chapter 2. An external function generator is used to trigger the laser and the DAQ simultaneously. To investigate where the band-shaped noise is coming from, the data is collected under different situations. In the first scenario, the laser is turned off, the probe is disconnected, and only the DAQ is triggered. The data collected by the DAQ is only the electronic noise from the DAQ itself, as shown in Fig. 3.2(a). This type of noise is random and is likely caused by the thermal noise of the electronics of the DAQ. In the second case, the laser is still off, the linear array transducer is connected with the DAQ, and the DAQ is triggered. The collected data is still noise and the band-shaped noise is clearly observed as shown in Fig. 3.2(b). The band-shaped noise is likely caused by the RF coupling of the high frequency electronic signals from the surrounding area through the linear array transducer. In the third case, the laser is on but not illuminating the sample (light is blocked), the transducer is connected, and the external function generator triggers both the DAQ and the laser. The data is shown in Fig. 3.2(c) and the band-shaped noise is more obvious. The driving circuity of the laser adds another source of RF signal interference to the DAQ of the PA system to increase the band-shaped noise. 20 Figure 3.2: Data acquisition for different cases. (a) Without transducer and laser. (b) With transducer and without laser. (c) With transducer and the laser is on but blocked. Next, we will apply RPCA to remove the band-shaped noise and compare the results with other denoising methods such as simple background subtraction and SVD. For the implementation of background subtraction, background PA signal is collected in advance when the laser beam is blocked and then subtracted from the subsequently collected PA signals to remove the background noise. To evaluate the reconstruction quality, peak signal-to-noise ratio (PSNR) is utilized which is expressed as: 1010 log ( )PSNR MAX MEAN (3.8) where MAX refers to the maximum pixel value of the image and MEAN is the mean value of a region outside the laser illumination area which is considered as the background. The first sample is human hair embedded in an optically and acoustically transparent gelatin phantom. Two pieces of hair are placed across each other and near parallel to the linear 21 array transducer. The frequency spectrum of PA signal from one channel and the band-shaped noise measured in the above-mentioned case 3 are analyzed and the results are shown below: Figure 3.3: Data comparison and power spectrum analysis. (a) Raw PA signal from one channel and the measured band-shaped noise. (b) Frequency spectrum analysis. Figure 3.3 describes the PA signal and band-shaped noise in time domain and frequency domain. It can be clearly seen that the band-shaped noise has a broad bandwidth which overlapped with the PA signal. Therefore, the filtering method is not suitable for removing the band-shaped noise. The data separation results are shown in Fig. 3.4. The obtained background contains mostly the band-shaped noise, successfully separated from the foreground. The PA signal from the hair is relatively high and the magnitude of the band-shaped noise is much lower than the PA signal. Figure 3.4: Data separation results. (a) Foreground data corresponds to the pure PA data. (b) Background data refers to the band-shaped noise. (c) The observation matrix is the original PA signals. 22 Figure 3.5: DAS reconstructions and line profile comparisons. (a, b, c, d) DAS reconstruction from the raw PA signals, background subtracted data, SVD denoised data and RPCA denoised data, respectively. (e) Horizontal line profile comparison (along the red arrow). (f) Vertical line profile comparison (along the red arrow). Different denoising methods are applied on the PA signal and the results are compared, including raw PA signal (no denoising), background subtraction, SVD, and RPCA. The pre-processed PA signals are then reconstructed by DAS and the obtained images are compared in Fig. 3.5. Among the four reconstructions, the RPCA denoised has the best visual quality. The 23 background subtraction in Fig. 3.5(b) shows more artifacts marked by a green arrow. Similarly, strong artifacts are observed and part of the signal is removed in the SVD (marked by green and white arrows, respectively). The line profiles marked by red arrows along the vertical and horizontal directions are compared in Figs. 3.5(e) and 3.5(f), respectively. In the horizontal line-profile, the SVD denoised curve (in green) shows a lot of noise. The curves of the background subtraction and RPCA show much less background noise than SVD. The PSNR of the reconstructions from the raw PA data, background subtracted data, SVD denoised data and RPCA denoised data are 35.94 dB, 40.58 dB, 36.94 dB and 47.22 dB, respectively. RPCA shows a ~11 dB improvement compared to the case from the raw PA data. The second sample is a 10-by-10 printed point array embedded in gelatin phantom. The band-shaped noise is extracted by RPCA and the results are shown in Fig. 3.6. Figure 3.6: Data separation results for the dots array data. (a) Foreground data corresponds to the pure PA data. (b) Background data refers to the band-shaped noise. (c) The observation matrix is the original PA data. 24 Figure 3.7: DAS reconstructions. (a, b, c, d) DAS reconstruction from the raw PA data, background subtracted data, SVD denoised data and RPCA denoised data, respectively. Figure 3.8: Line profile comparisons. (a) Horizontal line profile comparison. (b) Vertical line profile comparison. The results of different denoising methods are compared in Fig. 3.7. In Fig. 3.7(c), the first row is severely distorted and some points in the circled area are missing in SVD. The 25 reconstructions from background subtraction and RPCA are similar. Fig. 3.8 compares the line profiles along the red arrow directions. In the horizontal line profile, the SVD denoised curve (in green) deviates from the original curve (in black) and the fifth and sixth points are distorted. As for the PSNR, the original reconstruction is 38.46 dB, the background subtracted is 40.98 dB, the SVD denoised is 39.18 dB and the RPCA denoised is 41.44 dB. RPCA shows a ~3 dB improvement compared to the original reconstruction without denoising. Figure 3.9: Data separation results for the human forearm data. (a) Foreground data corresponds to the pure PA signals. (b) Background data refers to the band-shaped noise. (c) Observation matrix is the original PA signals. Furthermore, we implement the RPCA to the in vivo imaging of the human forearm which is submerged in a water tank. The incident laser with a wavelength of 750 nm illuminates the skin surface from the top and the transducer is positioned on the same side as the incoming light. The illumination area on the skin is a circle with a diameter of 3 cm. In order to generate the PA signals without damaging the skin, the local light fluence rate on the skin surface is kept below 10 mJ/cm2 26 under the ANSI safety limit [54]. The band-shaped noise is removed by the RPCA and the results are shown in Fig. 3.9. Figure 3.10: Background is the ultrasound image and foreground is PA reconstruction. (a, b, c, d) Foreground is the DAS reconstruction from the raw PA data, background subtracted data, SVD denoised data and RPCA denoised data, respectively. Figure 3.10 shows the reconstructed PA images (colored) after different denoising preprocess, overlaid on top of the corresponding ultrasound image (gray scale). The PA images I only used the data of in vivo imaging and did not conduct the experiments. 27 are shown in log scale. The skin surface and three veins can be seen in the PA images. In Fig. 3.10(a), the PA image from the raw data without denoising shows a very noisy background. In Fig. 3.10(b), the background subtraction shows additional background noise in the circled area. In Fig. 3.10(c), the SVD method still has strong background noise or artifacts. Fig. 3.10(d) shows that RPCA has obtained a relatively clean background. The circled region in Fig. 3.10(d) shows some residual noise but is lower than the other cases. This region is deep in tissue, which shows the effectiveness of RPCA. The PSNR of the original reconstruction, background subtraction one, SVD denoised one and RPCA denoised one are 33.31 dB, 40.03 dB, 36.9 dB and 43.71 dB, respectively. The enhancement of PSNR in RPCA is remarkable, which is ~10 dB. For the in vivo imaging, we can get a better quality of reconstruction by applying RPCA. In summary, RPCA is successfully implemented to remove the band-shaped noise in PA signal and the denoised results are compared with other methods. Since the background of different measurements cannot be exactly the same, a simple subtraction may induce additional noise which is more noticeable when the SNR of the PA data is low. The SVD denoising method requires that the magnitude of noise is comparable to that of the PA signal. In comparison, the RPCA method is applicable whether the SNR is low or high. Therefore, it is demonstrated that RPCA is an effective denoising method for the band-shaped noise and it can be further used for various applications of PAT. 28 Chapter 4: Variance-reduced stochastic gradient descent algorithm for photoacoustic image reconstruction As we mentioned before, PA image reconstruction from insufficient measurements suffers from artifacts and missing structures. In this chapter, we will apply iterative reconstruction method to improve the image quality. Our iterative method is based on a new discrete forward projection matrix and the VR-SGD algorithm. The proposed method is tested on both simulation and experimental data. 4.1 Introduction As an emerging technique in biomedical imaging, PAT achieves strong optical contrast and high ultrasound resolution in a single modality. In PAT, the image of the optical absorption map is reconstructed from the PA signals acquired by ultrasound probe. In this chapter, we focus on the iterative reconstruction algorithms which have advantages over direct reconstructions on the measurements with limited detection view. A necessary step of iterative methods is to model the forward propagation by a projection matrix which projects the PA signal from the optical absorption image. Then the iterative reconstruction process is formulated into a cost function that usually contains a data misfit term and a regularization term. The data misfit term usually describes the agreement of the image estimate with the measurements. A common choice for the data misfit term is the least-squares error between the projected PA data and the measured PA data. Also an absolute value of the error can be adopted so as to acquire a more robust solution. A regularization term reflects a prior knowledge about the image such as image smoothness or sparsity in some 29 transform domain. The regularization term can prevent overfitting issue, make the images smooth, and enhance the robustness of the algorithm. As for minimizing the cost function, the 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 is large. Paltauf et. al. used FGD on relatively small size of simulation and experimental PA data [39]. Since they did not apply a regularization term in the cost function, their method is affected by the insufficient and noisy PA data. Stochastic gradient descent (SGD) algorithm randomly chooses a smaller size of sub-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 [40]. Their cost functions contain a total variation (TV) regularization term. By utilizing analytical models and implementing the SGD algorithm to solve the TV based cost function, PAT image reconstruction was realized in both 2D and 3D space by Huang et. al. and Wang et. al., respectively [41, 42]. Even though 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 [55]. In order to reduce the effect of the large variance, one popular approach is to gradually decrease the step size of the gradient decent algorithm [40]. However, this approach requires a careful tuning of the step size and it can only achieve a sublinear convergence rate. In recent years, 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 CT image reconstruction, and have been shown to have much faster convergence than the SGD schemes [56]. 30 To the best of our knowledge, VR-SGD has not been applied to PAT. In our work, we apply VR-SGD to PAT reconstruction and compare it with other algorithms. As for the modeling the forward propagation, Huang et al. reported a full-wave based approach [41]. Their forward and backward projection operators were based on the exact PA wave equation and built by the k-space pseudospectral method. Inhomogeneous speed-of-sound, mass density, and acoustic attenuation were taken into consideration for the medium. Wang et al. proposed a discrete model which incorporates the spatial and acousto-electric impulse response (EIRs) of an ultrasound transducer [42]. Since their models are based on the wave propagation and consider both the homogeneous and heterogeneous environment, they are very complex. A relatively simple model which was based on the arc projection and the homogeneous medium was used in Ref. [39, 40, 57]. A discrete projection matrix was constructed to project the PA signal from the optical absorption image in the time domain, but there was a mismatch between the projection matrixes computed PA signal and the measured PA signal, which could result in inaccuracy when computing the misfit term in the cost function. We will improve the accuracy of the projection matrix and match it with the measured data. 4.2 Theoretical background 4.2.1 Photoacoustic imaging model In order to apply iterative PAT image reconstruction method, a discrete model which describes the relationship between the measured PA signal and the optical absorption distribution needs to be built. The background of PA imaging model has been presented in Chapter 2. Based on this model, the projection matrix can be developed. 31 4.2.2 Building the basic discrete forward projection matrix The optical absorption map can be treated as a 2D matrix. The size of the absorption image ( )x r is m-by-n, 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. After vectorization, we get a 1D absorption vector x which has the size of mn-by-1. Suppose an ultrasound linear array probe has Q channels and the length of the discrete signal for each channel is T. A forward projection matrix iA relates ( , )ig r t at channel i to x [40]: , 1,i ig A x i Q (4.1) The projection matrix iA has the size of T-by-mn and can be obtained as [40]: ( , ) m ax{1 / ( ) , 0}, 1,i sA j k d v t j j T (4.2) where d is the distance from the kth element on the absorption image to the ith transducer channel, sv is the speed of sound, t is the time step which is equal to the reciprocal of the sampling rate. Here ( , )iA j k represents the weighting factor from the kth element of the absorption image to the jth element of the sampling points of the ith transducer channel. All the measurements at different transducer channels can also be lexicographically reshaped into a vector g . Hence, the forward propagation model can be represented as: g Ax (4.3) where the full projection matrix A can be written as: 1 2 3TQA A A A A (4.4) The size of the full projection matrix A is TQ-by-mn. When the size of the absorption map matrix or the measurements is large, the size of A would be huge. 32 4.2.3 Gradient method for solving cost function According to Eq. 4.3 the absorption image ( )x r is the solution of a system of linear equations. However, the system is underdetermined and the measurements can be very noisy. To solve this problem, we consider iterative methods for reconstructing the absorption image ( )x r from the measurements by minimizing a cost function which is based on the least-squares error between the projected PA data and the measured PA data [58]: 221( )2J x Ax g (4.5) where 2 is the entry-wise L2-norm. Then the estimated image *x is obtained by minimizing the cost function above: 2*21arg min2xx Ax g (4.6) Since the linear array probe is only able to detect a limited-view of the acoustic waves, the solution of minimizing Eq. 4.5 without a regularization term would generate many artifacts. Many different forms of regularization are available for solving inverse problems. The TV regularization promotes piecewise-constant solutions, which has been proved to be a very effective assumption in various imaging processing applications and can be estimated as [59]: 2 2, 1, , , 12 2( ) ( )m ni j i j i j i ji jTV x x x x x (4.7) where ,i jx refers to the ith row and jth column of image ( )x r . After adding a TV regularization term, the cost function can be described as: 33 221( )2J x Ax g TV x (4.8) where is a tuneable parameter which balances the data misfit and the regularization effects. Then a new solution of estimated image *x is: 2*21arg min2xx Ax g TV x (4.9) In FDG algorithm, the cost function is gradually reduced by taking a small step along the negative full gradient direction ( )J x : ( ) ( ) ( )TJ x A Ax g TV x x (4.10) The partial derivative of TV can be calculated by an approximation [40]. Hence, the gradient descent step to update the image estimate takes the following form: 1( ) [ ( ) ( ) ]k k k k T k k kk kx x J x x A Ax g TV x x (4.11) where k and k+1 are iteration numbers and k is the step size. For typical image and measurement sizes, the size of the projection matrix, A, can be very large. Therefore, burden of computing the full gradient, which requires multiplication with A and its transpose, can be very heavy. In SGD algorithm, the cost function is rewritten by using the row sub-matrices of A: 221 11 1 1( ) ( ) { ( )}2Q Qi i ii iJ x j x A x g TV xQ Q (4.12) The corresponding SGD gradient direction ( )kij x and the updating process can be written as: 1( ) [ ( ) ( ) ]k k k k T k k kk i k i i ix x j x x A A x g TV x x (4.13) The sub-matrix iA , corresponding to the ith transducer channel position, can be randomly or sequentially selected from the transducer channels in each iteration. In Ref. [40], the sub-matrices 34 are sequentially selected during the iterations. However, the convergence behavior of SGD is poor. This is because even though ( )ij x is equal to ( )J x in its expected value, its variance is very high [55]. One has to gradually reduce the step size and still the convergence will be very slow [60]. 4.3 Methodology 4.3.1 A new improved projection matrix In Section 4.2.2, we can see that the matrix A projects ( , )g r t from ( )x r . Therefore, the data misfit term in the cost function should be 221 2 Ax g . However, in Ref. [40, 57], the PA data ( , )p r t was directly used to replace the ( , )g r t in the cost function for simplicity. This is not accurate and can cause errors. Next, we will obtain a more accurate matrix to project ( , )p r t from ( )x r . The relationship between ( , )g r t and ( , )p r t can be obtained from Eq. 2.7 of Chapter 2 as: 0 01 ( , ) ( , ) 1 ( , )( , )4 4 I z I zg r t g r t g r tp r tC t t t C t t (4.14) As for the experimental cases, the central frequency of transducer is 7.2 MHz and the pulse width of PA signal is 0.14 µs. Suppose the distance from the sample to transducer is 2 cm and the speed of sound is 1540 m/s, the travelling time t of acoustic waves is 13 us. Since the pulse width of signal is much smaller than the travelling time, ( , )g r tt can be neglected. Therefore, according to Eq. 4.14, an improved forward projection matrix which considers a derivative and a 1 t attenuation matrix is described as: 35 i ip EDg (4.15) where the matrix E and matrix D correspond to the 1 t and ( , )g r t t in Eq. 4.14, respectively. The constant parts in Eq. 4.14 are neglected and the sizes of E and D matrices are both T-by-T. More specifically, the elements ,u we of E are: ,1 ,0,u wu u weelse (4.16) where w refers to sequence number of the measurements and matrix E is diagonal. For the matrix D, a forward difference approximation and a backward difference approximation are adopted in the first and the last row. Other rows use central difference approximation and an example of D with a size of 4-by-4 is shown below: 1 1 0 01 / 2 0 1 / 2 00 1 / 2 0 1 / 20 0 1 1D (4.17) A new modified projection matrix iM relating the acoustic waves of one transducer channel ip to the absorption image ( )x r is given by: i i ip EDA x M x (4.18) 4.3.2 VR-SGD algorithm Algorithm 4.1: The VR-SGD algorithm. is the regularization parameter. 1 iL is the step size. Input: a zero matrix or a DAS reconstruction 0x with a threshold Output: final image Nx after N iterations 36 01 0x x For k = 0 to Q - 1 A round of basic Select 1, ,i Q SGD for image 1( )T ki i id M M x p Stochastic gradient direction estimate x 11 1( )k kTV ix x d L End For j = 2 to N 11kjx x ( )TM Mx p Full gradient direction VR updates 0jx x for image 1kjx For k = 0 to 2Q – 1 Select 1, ,i Q with probability with i iiL L ( )T ki i jd M M x x VR gradient direction 1( )k kj TV j ix x d L End End 1kN Nx x In Section 4.2.3, we have mentioned that SGD has poor convergence. Next, we will introduce the VR-SGD algorithm (shown in Algorithm 4.1) which has good convergence behavior. The first part of this algorithm starts by performing a round of basic proximal SGD updates and acquire an image estimate x , which is known to lead to a faster initial convergence [61, 62]. Then it will switch to the second part: proximal VR-SGD updates. The core idea of the second part is to use the knowledge of the full gradient ( ( ) J x ) and the difference between the gradient of the current image estimate kx and the recent image estimate x to build the update direction: 37 ( ) ( ) ( ) ( ) ( )k T ki i i id j x j x J x M M x x J x (4.19) The corresponding updating process is: 1( ( ) ( ))k k T kk i ix x M M x x J x (4.20) In order to reduce the computational burden, the full gradient ( ) J x ( in Algorithm 4.1) can be updated after 2Q inner iterations [63]. As the estimated image approaches a local optimum, the variance of the difference between ( )kij x and ( )ij x decreases to zero and the direction approaches the full gradient. Therefore, the variance is reduced and the convergence rate is substantially improved compared with the basic SGD [56]. Since the initial state of the iteration affects the whole optimization process, the initial image should be carefully selected. As for the implementation of VR-SGD to PAT reconstruction, we can either use a zero matrix or a DAS reconstructed image as the initial image. The closer of the initial image to the real optical absorption map, the fewer number of iterations the algorithm will need to achieve a certain reconstruction quality. Furthermore, we can use a proximal operator to solve the differentiation of the TV part, instead of using the approximation in [40]. The proximal operator of TV can be expressed as [63]: 221( ) arg min{ ( )}2TVzx z x TV z (4.21) In this thesis, the proximal operator is solved without taking gradients [64]. Besides, the detection views of the inner iteration of the second part are randomly selected with a probability proportional to the Lipschitz constant iL of their gradient [56]. Lipschitz constant iL is the largest eigenvalue of Ti iM M and can be used as the step size (1 iL in Algorithm 4.1), which is known to lead to a faster convergence [62]. 38 4.4 Simulation validation Figure 4.1: PAT simulation setup and signals generation comparison. (a) Simulation setup. (b) PA signals comparison. The accuracy of projection matrices M and A is first compared with k-wave generated PA signal [65]. 2D simulation with a linear array transducer is conducted using Matlab (version 8.6.0.267246) k-wave toolbox on a PC with 3.7 GHz CPU and 16.0 GB memory. The medium is assumed to be homogeneous and the speed of sound is 1540 m/s. Fig. 4.1(a) shows the simulation setup, where the sample is located in the center and the linear array transducer at the top. The source is a 1.536 mm square with uniform density and the linear array transducer has 40 elements. The simulated optical absorption matrix is 40×40 pixels. Fig. 4.1(b) shows the PA signals generated by the M and A projection matrices and k-wave toolbox in black, red and blue curves, respectively. As we can see, the PA signal projected by the M matrix is much closer to the k-wave signal than that projected by the A matrix. Therefore, the matrix M provides a more accurate forward projection from the optical absorption map to the PA signal. Next, we compare the reconstruction results achieved by different iterative algorithms. The simulation setup is similar to Fig. 4.1(a). The optical absorption matrix is increased to 128 ×128 39 pixels and the objects consist of three solid discs and four points in a row. The linear array transducer consists of 128 elements with a pitch size of 0.08 mm to match with experiment. Table 4.1 lists the different algorithms that are compared, and whether they have used the projection matrix A or M, and whether they have used gradient estimation or proximal operator to handle the TV regularization term. The SGD algorithm with projection matrix A, denoted as A-SGD in this chapter, is the method used in Ref [40] and the corresponding update step is: 1( ) [ ( ) ( ) ]k k k k T k k kk i k i i ix x j x x A A x p TV x x (4.22) The update steps of the SGD algorithm with projection matrix M, denoted as M-SGD, can be described as: 1( ) [ ( ) ( ) ]k k k k T k k kk i k i i ix x j x x M M x p TV x x (4.23) Table 4.1: Main features of different iterative reconstruction algorithms. IR A-SGD M-SGD VR-SGD Matrix A √ Matrix M √ √ √ Gradient Estimation √ √ Proximal Operator √ Feature Algorithm 40 Figure 4.2: PAT reconstructions. (a) Normalized DAS reconstruction with a threshold of 0.16. (b, c, d, e, f) Reconstructions by iterative algorithms: IR, A-SGD, M-SGD and VR-SGD, respectively. Figure 4.3: Comparisons of PAT reconstructions. (a) Line profile comparison. (b) RRMSE comparison. When the regularization parameter λ is set to zero in Eq. 4.23, it is a regular iterative (IR) algorithm and the updating step is [39]: 1( ) ( )k k k k T kk i k i i ix x j x x M M x p (4.24) In iterative reconstructions, the regularization parameter needs to be selected carefully. The IR algorithm with 0 is not an effective method for noisy and under-sampled measurements. In the other methods, can be tuned by a strategy mentioned in [40], which requires a gradual 41 decrease of . In our simulation, the initial value of is set as 0.00005 and then an adaptive tuning rule is used as: 0.00005 / ,1 100.000005, 11k kk (4.25) where k is the number of iteration. Figures 4.2(a)-4.2(e) show the reconstructed images with the different algorithms, DAS, IR, A-SGD, M-SGD, and VR-SGD, respectively. In Fig. 4.2(a), it can be clearly seen that the normalized DAS reconstruction with a threshold of 0.16 has strong artifacts due to the incomplete measurement data. With the iterative methods, the artifacts are all reduced compared with DAS. In Fig. 4.2(b), the IR algorithm shows relatively strong artifacts than the other iterative algorithms, because the cost function of IR does not adopt a regularization term and the artifacts cannot be effectively removed. In Fig. 4.2(c), A-SGD only reconstructs part of the structures, only the three big solid discs are recovered and the four small points are missed. This is mainly because the project matrix A does not pair with the PA signals. In Figs. 4.2(d) and 4.2(e), both M-SGD and VR-SGD can recover most of the structures. Thus the structure can be recovered with a good visual quality only when the projection matrix and the PA signals are matched. Furthermore, VR-SGD shows fewer artifacts than M-SGD. The line profiles that go through the middle of the reconstructed objects are compared in Fig. 4.3(a). VR-SGD can recover all the seven objects with sharp edges and the corresponding line profile is smooth. The line profile from A-SGD has sharp edges because the negative parts of the signal are simply removed in each iteration according to Ref. [40]. Since our simulations are noise-free, the line profile of IR has a good shape and the four small points can be recovered. However, the artifacts of IR reconstruction are strong. A quantitative comparison by examining the evolution 42 of the relative root-mean-square of reconstruction error (RRMSE) is shown in Fig. 4.3(b). The equation of RRMSE is described as: 2 21 1 1 1[ ( , ) ( , )] ( , )m n m ni j i jRRM SE x i j R i j R i j (4.26) where x is the reconstructed image and R is the original sample absorption map in the simulation. The RRMSE of A-SGD reconstruction saturates after 5 iterations, but those of the other iterative algorithms keep decreasing with more iterations. After 60 iterations, all the iterative algorithms get saturated. After 60 iterations, the RRMSEs are 0.9875, 0.8536, 0.8478 and 0.8327 for the A-SGD, IR, M-SGD and VR-SGD, respectively. Therefore, the reconstruction by the VR-SGD is the closest to the original image. The corresponding computational time of A-SGD, IR, M-SGD and VR-SGD are 699 s, 619 s, 626 s and 810 s, respectively. Since VR-SGD needs to calculate the full gradient, it costs more time. The peak signal-to-noise ratio (PSNR) is also calculated as: 1010 log ( )PSN R M AX M SE (4.27) where MAX refers to the maximum value of the image and MSE is the mean squared error: 21 11[ ( , ) ( , )]m ni jM SE x i j R i jm n (4.28) Table 4.2: PSNR comparisons by different algorithms after 60 iterations. Algorithm DAS IR A-SGD M-SGD VR-SGD PSNR [dB] 11.48 12.95 11.21 13.06 13.26 The PSNR are compared in Table 4.2. The reconstruction by A-SGD has a PSNR of 11.21 dB and the reconstruction by the proposed VR-SGD has the highest PSNR of 13.26 dB, almost a 2.05 dB improvement. Therefore, the VR-SGD is proved to be a more accurate algorithm than the 43 other algorithms and can reconstruct images with a higher PSNR and a lower RRMSE. It is necessary to mention that since a linear array transducer is used, the parts of the absorbers that are parallel to the transducer array is better reconstructed than the parts that are vertical to the transducer array. 4.5 Experimental implementation The projection matrix M based VR-SGD algorithm is also verified through PAT experiments. Details about the PAT experimental system can be found in Chapter 2. For all the experimental data, a normalized DAS reconstruction with a threshold is used as the initial state of VR-SGD. Also, an adaptive tuning strategy for the regularization parameter is adopted. Since the initial value of is highly related to the size of the reconstructed image, a good initial value of 0.005 is found for an image with a size of 740×500 pixels. Therefore, our strategy for tuning the TV term is: 0.005 / ,1 100.0005, 11k kk (4.29) In order to evaluate the quality of reconstruction, we analyzed the PSNR and named it as PSNR1 so as to separate it from Eq. 4.27. The PSNR1 can be expressed as: 101 10 log ( )PSNR MAX MEAN (4.30) where MAX refers to the maximum value of the image and MEAN is the mean value of a selected region outside the laser illumination area. 44 Figure 4.4: PAT reconstructions for the points array sample and the line profile comparisons. (a) Normalized DAS reconstruction with a threshold of 0.16. (b) M-SGD reconstruction. (c) VR-SGD reconstruction. (d) Horizontal line profile comparison. (e) PSNR1 comparison. A 10-by-10 printed dots array embedded in an optically and acoustically transparent gelatin phantom is imaged. The radius of each dot is 0.1 mm and the spacing between them is 1 mm. The sample plane is parallel to the linear array transducer plane. Figure 4.4(a) shows the normalized DAS reconstruction with a threshold of 0.16 and most of the artifacts are removed. Figures 4.4(b) and 4.4(c) show the M-SGD and VR-SGD reconstructed images, respectively, after six iterations. The artifacts can also be removed and the background is cleaner than DAS. Line profiles along the direction indicated by the arrows are plotted as red, green, and blue curves for DAS, M-SGD, and VRSGD, respectively, in Fig. 4.4(d). In the green and blue curves, 10 peaks can be clearly distinguished. However, the red curve from DAS reconstruction fails to distinguish all 10 dots mainly because of artifacts and background noise. Furthermore, the peaks of the green and blue curves have a narrower width than that of the red, which indicates that the iterative reconstruction 45 methods have better reconstruction abilities than that of the direct method. Fig. 4.4(e) shows the PSNR1 versus iteration. M-SGD saturates after 2 iterations and the corresponding PSNR1 is 68.65 dB. VR-SGD acquires a much higher PSNR1, 82.48 dB, after 3 iterations. On the other hand, the PSNR1 of DAS reconstruction is 38.46 dB, which is relatively low. Similar to our simulations results, the experiment results also show that VR-SGD algorithm can lead to better image reconstruction with less artifacts and a cleaner background. In vivo imaging of the human forearm is also carried out and the image is reconstructed by VR-SGD. The incident laser with a wavelength of 750 nm illuminates the skin surface from the top and the transducer is positioned on the same side as the incoming light. The light fluence rate on the skin surface is kept below 10 mJ/cm2 under the ANSI safety limit [54]. The reconstructed images are shown in Figs. 4.5(a)-4.5(c) for DAS, M-SGD, and VR-SDG, respectively. The PA image is in color and overlaid on top of the gray scale ultrasound image. The reconstructed skin and blood vessels are marked with red arrows. The normalized DAS reconstruction uses a threshold of 0.16. For the iterative reconstructions, the skin and three veins can be clearly recovered with less noise than DAS. The PSNR1 of normalized DAS reconstruction is 34.26 dB, which is relatively low. The M-SGD saturates after 2 iterations and has a PSNR1 of 68.14 dB. The VR-SGD has a PSNR1 of 77.27 dB after 6 iterations. Hence, compared with regular SGD based iterative reconstruction algorithm, the VR-SGD can reconstruct images with a higher PSNR1. Therefore, the proposed VR-SGD algorithm is successfully applied to the in vivo imaging and the artifacts and background noise are significantly reduced. 46 Figure 4.5: Background is the ultrasound image and foreground PA reconstruction. (a) Foreground is a normalized DAS reconstruction with a threshold of 0.16. (b) Foreground is the M-SGD reconstruction after six iterations. (c) Foreground is the VR-SGD reconstruction after six iterations. In this chapter, VR-SGD was implemented to the PA image reconstruction. A new projection matrix was developed to pair up with the measured PA data. Our proposed method was demonstrated to mitigate the artifacts generated be a limited-view detection and improve the image quality. The VR-SGD was compared with other iterative reconstruction algorithms for both simulation and experimental data. The results were quantitatively evaluated by different standards and the advantages of VR-SGD were demonstrated. I only used the data of in vivo imaging and did not conduct the experiments. 47 Chapter 5: Conclusion and future work As a promising hybrid imaging technology in biomedical field, PAT is able to achieve centimeter level of penetration depth and high ultrasound resolution. In recent years, PAT has been increasingly applied to clinical applications. This technique is based on the acquisition of acoustic waves that are generated by thermoelastic expansion caused by illumination of a short pulsed laser. The generated acoustic pressure waves can be collected by different types of transducers which provide either a full or a limited view. Linear array transducers are commonly used in clinical ultrasound applications but the reconstructions always suffer from the artifacts which are caused by the insufficient measurements. Besides, the band-shaped noise is prominent in the radiofrequency data acquired by the linear array transducers and it could affect the reconstructions. This thesis addresses those problems by implementing RPCA to reduce the band-shaped noise and VR-SGD with a new projection matrix to improve the iterative image reconstruction. 5.1 Significance of work A novel application of RPCA is implemented to remove the band-shaped noise in the PAT system. It decomposes the measured PA signal matrix into a low-rank background matrix and a sparse foreground matrix without destroying the structure of the input data. Based on the foreground matrix extracted from RPCA, we can reconstruct an image with less background noise. We tested the feasibility and effectiveness of our proposed method by comparing it with background subtraction and SVD method. RPCA shows better results than the other methods and can acquire a clean background. Thus RPCA has great potential in denoising PA signals. It can be further used in various fields of medical detection and bio-medical imaging applications. 48 Linear array transducers have limited view and suffer from artifacts and missing structures because of the incomplete measurements. A novel application of VR-SGD is implemented for PA imaging to improve the reconstruction with limited view data. To perform VR-SGD, an accurate projection matrix is developed to project PA signals from the optical absorption map. VR-SGD is shown to have better performance than several other iterative methods used in PAT. Using different evaluation criteria, such as PSNR, RRMSE and line profile comparisons, the reconstructions from various algorithms are compared for both simulation and experimental data and the advantages of VR-SGD are demonstrated. VR-SGD algorithm is a good choice for the linear array transducer based PAI reconstructions. The RPCA denoising method and the VR-SGD reconstruction algorithm have been implemented on PAT imaging and great improvement on image quality is demonstrated. Those approaches have a high potential for the clinical applications of PAT. 5.2 Future work and improvement The VR-SGD used in this thesis can remove the artifacts caused by insufficient measurements, but it cannot recover the missing structures. In order to reconstruct the missing structures, we can implement VR-SGD on a two-probe measurement by our PAT system using two linear array ultrasound probes [67]. The iterative reconstruction will be applied on the data from both probes. A transformation matrix will first be acquired between the two probes and it can help rotate the image. In the iterative update of the image through VR-SGD, the image will be reconstructed corresponding to one probe, and then transformed into another view for updating by the other probe. With the two probe data collection and the VR-SGD reconstruction, the missing structures caused by a limited-view detection can potentially be recovered. 49 The reconstruction is currently processed in time domain and is time consuming. In the future, in order to reduce the computation burden, we can conduct the optimization in another domain by Discrete Cosine Transform (DCT). The reconstruction process can be performed by manipulating the DCT coefficients. 50 Bibliography 1. R. A. Kruger, D. R. Reinecke, and G. A. Kruger, “Thermoacoustic computed tomography--technical considerations,” Med. Phys. 26(9), 1832–1837 (1999). 2. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). 3. P. Beard, “Biomedical photoacoustic imaging,” Interface Focus 1(4), 602–631 (2011). 4. L. V. Wang, "Multiscale photoacoustic microscopy and computed tomography," Nat. Photonics 3, 503–509 (2009). 5. L. V. Wang, “Prospects of photoacoustic tomography,” Med. Phys. 35, 5758–5767 (2008). 6. A. Oraevsky, S. Ermilov, K. Mehta, T. Miller, B. Bell, E. Orihuela, and M. Motamedi, “In vivo testing of laser optoacoustic system for image-guided biopsy of prostate,” In Biomedical Optics 60860B (2006). 7. D. Piras, C. Grijsen, P. Schütte, W. Steenbergen, and S. Manohar, “Photoacoustic needle: minimally invasive guidance to biopsy,” J. Biomed. Opt. 18(7), 070502 (2013). 8. G. Xu, Z. X. Meng, J. D. Lin, C. X. Deng, P. L. Carson, J. B. Fowlkes, C. Tao, X. Liu, and X. Wang, “High resolution Physio-chemical Tissue Analysis: Towards Non-invasive In Vivo Biopsy,” Sci. Rep. 6, 16937 (2016). 9. M. A. Lediju Bell, N. P. Kuo, D. Y. Song, J. U. Kang, and E. M. Boctor, “In vivo visualization of prostate brachytherapy seeds with photoacoustic imaging,” J. Biomed. Opt. 19(12), 126011 (2014). 10. W. Xia, D. I. Nikitichev, J. M. Mari, S. J. West, R. Pratt, A. L. David, S. Ourselin, P. C. Beard, and A. E. Desjardins, “Performance characteristics of an interventional multispectral 51 photoacoustic imaging system for guiding minimally invasive procedures,” J. Biomed. Opt. 20(8), 086005 (2015). 11. J. M. Mari, W. Xia, S. J. West, and A. E. Desjardins, “Interventional multispectral photoacoustic imaging with a clinical ultrasound probe for discriminating nerves and tendons: an ex vivo pilot study,” J. Biomed. Opt. 20(11), 110503 (2015). 12. C. P. Favazza, O. Jassim, L. a Cornelius, and L. V. Wang, "In vivo photoacoustic microscopy of human cutaneous microvasculature and a nevus.," J. Biomed. Opt. 16, 16015 (2011). 13. R. A. Kruger, D. R. Reinecke, and G. A. Kruger, “Thermoacoustic computed tomography--technical considerations,” Med. Phys. 26(9), 1832–1837 (1999). 14. A. G. Bell, "Upon the production and reproduction of sound by light," Telegr. Eng. J. Soc. Of 9, 404–426 (1880). 15. J. Tyndall, "Action of an intermittent beam of radiant heat upon gaseous matter," Proc. R. Soc. Lond. 31, 307–317 (1880). 16. W. C. Röntgen, "On tones produced by the intermittent irradiation of a gas," (1881). 17. A. C. Tam, "Applications of photoacoustic sensing techniques," Rev. Mod. Phys. 58, 381 (1986). 18. C. K. N. Patel and A. C. Tam, "Pulsed optoacoustic spectroscopy of condensed matter," Rev. Mod. Phys. 53, 517 (1981). 19. T. Bowen, "Radiation-induced thermoacoustic soft tissue imaging," in 1981 Ultrasonics Symposium (IEEE, 1981), pp. 817–822. 20. R. A. Kruger, P. Liu, Y. “‘Richard’” Fang, and C. R. Appledorn, "Photoacoustic ultrasound (PAUS)—Reconstruction tomography," Med. Phys. 22, 1605–1609 (1995). 52 21. R. O. Esenaliev, A. A. Karabutov, F. K. Tittel, B. D. Fornage, S. L. Thomsen, C. Stelling, and A. A. Oraevsky, "Laser optoacoustic imaging for breast cancer diagnostics: limit of detection and comparison with x-ray and ultrasound imaging," in BiOS’97, Part of Photonics West (International Society for Optics and Photonics, 1997), pp. 71–82. 22. A. A. Oraevsky, V. A. Andreev, A. A. Karabutov, R. D. Fleming, Z. Gatalica, H. Singh, and R. O. Esenaliev, "Laser optoacoustic imaging of the breast: detection of cancer angiogenesis," in BiOS’99 International Biomedical Optics Symposium (International Society for Optics and Photonics, 1999), pp. 352–363. 23. R. A. Kruger, K. D. Miller, H. E. Reynolds, W. L. Kiser Jr, D. R. Reinecke, and G. A. Kruger, "Breast Cancer in Vivo: Contrast Enhancement with Thermoacoustic CT at 434 MHz—Feasibility Study 1," Radiology 216, 279–283 (2000). 24. G. Ku and L. V. Wang, "Scanning thermoacoustic tomography in biological tissue," Med. Phys. 27, 1195–1202 (2000). 25. P. C. Beard, "Photoacoustic imaging of blood vessel equivalent phantoms," in International Symposium on Biomedical Optics (International Society for Optics and Photonics, 2002), pp. 54–62. 26. V. E. Gusev and A. A. Karabutov, "Laser optoacoustics," NASA STIRecon Tech. Rep. A 93, (1991). 27. F. A. Duck, "Physical properties of tissue: a comprehensive reference book. 1990," Lond. UK Acad. (n.d.). 28. D.-K. Yao, C. Zhang, K. Maslov, and L. V. Wang, “Photoacoustic measurement of the Grüneisen parameter of tissue,” J. Biomed. Opt. 19(1), 017007–017007 (2014). 53 29. G. J. Diebold, T. Sun, and M. I. Khan, "Photoacoustic monopole radiation in one, two, and three dimensions," Phys. Rev. Lett. 67, 3384–3387 (1991). 30. C. Li and L. V. Wang, "Photoacoustic tomography and sensing in biomedicine," Phys. Med. Biol. 54, R59 (2009). 31. J. D. Jackson, Classical Electrodynamics (Wiley, 1999). 32. G. B. Arfken, H. J. Weber, and F. E. Harris, Mathematical Methods for Physicists: A Comprehensive Guide (Academic press, 2011). 33. H. Feshbach and P. M. Morse, Methods of Theoretical Physics (McGraw-Hill Interamericana, 1953). 34. J. Zhang, M. A. Anastasio, P. J. La Riviere, and L. V. Wang, “Effects of different imaging models on least-squares image reconstruction accuracy in photoacoustic tomography,” IEEE Trans. Med. Imaging 28, 1781–1790 (2009). 35. M. H. Xu, Y. Xu, and L. H. V. Wang, "Time-domain reconstruction-algorithms and numerical simulations for thermoacoustic tomography in various geometries," IEEE Trans. Biomed. Eng. 50, 1086-1099 (2003). 36. Y. Xu, D. Z. Feng, and L.-H. V. Wang, "Exact frequency-domain reconstruction for thermoacoustic tomography-I: planar geometry," IEEE Trans. Med. Imaging 21, 823-828 (2002). 37. M. Xu, and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71(1), 016706 (2005). 38. Z. Wang, J. Li and R. Wu, “Time-delay-and time-reversal-based robust capon beamformers for ultrasound imaging,” IEEE Trans. Med. Imaging 24(10), 1308-1322 (2005). 54 39. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am. 112, 1536-1544 (2002). 40. Y. Zhang, Y. Wang, and C. Zhang, “Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction,” J. Ultrasonics 52(8), 1046-1055 (2012). 41. C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media,” IEEE Trans. Med. Imaging 32(6), 1097–1110 (2013). 42. K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Bio. 57, 5399–5423 (2012). 43. E. R. Hill, W. Xia, M. J. Clarkson and A. E. Desjardins, “Identification and removal of laser-induced noise in photoacoustic imaging using singular value decomposition,” Opt. Express. 8, 68-77 (2017). 44. P. Adam and P. J. L. Rivière, "Comparison of photoacoustic image reconstruction algorithms using the channelized Hotelling observer," Journal of biomedical optics 18(2), 026009-026009 (2013). 45. E. J. Candès, X. D. Li, Y. Ma and J. Wright, "Robust principal component analysis?," Journal of the ACM (JACM) 58(3), 11 (2011). 46. J. K. Cheng, K. Chen and M. D. Sacchi, "Application of Robust Principal Component Analysis (RPCA) to suppress erratic noise in seismic records," SEG Technical Program Expanded Abstracts, 4646-4651 (2015). 47. Z. Guo, L. Li, and L. V. Wang, “On the speckle-free nature of photoacoustic tomography,” Med. Phys. 36(9), 4084–4088 (2009). 55 48. N. Ryo, R. Yamazaki, and Y. Saijo, "Adaptive Spatial Filtering with Principal Component Analysis for Biomedical Photoacoustic Imaging," Physics Procedia 70, 1161-1164 (2015). 49. L. L. Pan, Photoacoustic Imaging for Prostate Brachytherapy (University of British Columbia, 2014). 50. Z. C. Lin, A. Ganesh, J. Wright, L. Q. Wu, M. M. Chen, and Y. Ma, "Fast convex optimization algorithms for exact recovery of a corrupted low-rank matrix," Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP) 61(6), (2009). 51. Y. Y. Xu, and W. T. Yin, "A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion," SIAM Journal on imaging sciences 6(3), 1758-1789 (2013). 52. J. C. Nash, "The singular-value decomposition and its use to solve least-squares problems," Compact Numerical Methods for Computers: Linear Algebra and Function Minimisation, 30-48 (1990). 53. N. Yuri, “Gradient methods for minimizing composite objective function,” CORE discussion paper, Université Catholique de Louvain, Center for Operations Research and Econometrics (2007). 54. American National Standards Institute, Safe Use of Lasers, ANSI Z136.1-2000, (Laser Institute of America, 2000). 55. L. Bottou, “Stochastic gradient tricks,” In: Montavon, G., Orr, G.B., Müller, K.-R. (eds.) Neural Networks, Tricks of the Trade, Reloaded. Lecture Notes in Computer Science, Vol. 7700, (Springer, 2012), pp. 430–445. 56 56. 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). 57. C. Zhang, Y. Wang, and Y. Wang, “A photoacoustic image reconstruction method using total variation and nonconvex optimization,” Biomed. Eng. Online 13(1), 117 (2014). 58. L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena, 60(1-4), 259-268 (1992). 59. C. R. Vogel, Computational methods for inverse problems (Society for Industrial and Applied Mathematics, 2002). 60. F. Bach and E. Moulines, “Non-strongly-convex smooth stochastic approximation with convergence rate O(1/n),” In Advances in Neural Information Processing Systems, (NIPS, 2013), pp. 773-781. 61. R. Johnson and T. Zhang, “Accelerating stochastic gradient descent using predictive variance reduction,” in Advances in Neural Information Processing Systems 26, C.J.C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, Eds., (Curran Associates, Inc., 2013), pp. 315–323. 62. L. Xiao, and T. Zhang, “A proximal stochastic gradient method with progressive variance reduction,” Siam Journal on Optimization, 24, 2057–2075 (2014). 63. P. L. Combettes, and J.-C. Pesquet, “Proximal splitting methods in signal processing,” Fixed-point algorithms for inverse problems in science and engineering, (Springer, 2011), pp. 185-212. 64. C. Antonin, “An algorithm for total variation minimization and applications,” Journal of Mathematical imaging and vision, 20(1), 89-97 (2004). 57 65. 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-021314 (2010). 66. S. Tzoumas, A. Rosenthal, C. Lutzweiler, D. Razansky, and V. Ntziachristos, “Spatiospectral denoising framework for multispectral optoacoustic imaging based on sparse signal representation,” Med. Phys. 41(11), 113301 (2014). 67. 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(12), 12755–12768 (2016).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Investigating signal denoising and iterative reconstruction...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Investigating signal denoising and iterative reconstruction algorithms in photoacoustic tomography Cheng, Jiayi 2017
pdf
Page Metadata
Item Metadata
Title | Investigating signal denoising and iterative reconstruction algorithms in photoacoustic tomography |
Creator |
Cheng, Jiayi |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | Photoacoustic tomography (PAT) is a promising biomedical imaging modality that achieves strong optical contrast and high ultrasound resolution. This technique is based on the photoacoustic (PA) effect which refers to illuminating the tissue by a nanosecond pulsed laser and generating acoustic waves by thermoelastic expansion. By detecting the PA waves, the initial pressure distribution that corresponds to the optical absorption map can be obtained by a reconstruction algorithm. In the linear array transducer based data acquisition system, the PA signals are contaminated with various noises. Also, the reconstruction suffers from artifacts and missing structures due to the limited detection view. We aim to reduce the effect of noise by a denoising preprocessing. The PAT system with a linear array transducer and a parallel data acquisition system (DAQ) has prominent band-shaped noise due to signal interference. The band-shaped noise is treated as a low-rank matrix, and the pure PA signal is treated as a sparse matrix, respectively. Robust principal component analysis (RPCA) algorithm is applied to extract the pure PA signal from the noise contaminated PA measurement. The RPCA approach is conducted on experiment data of different samples. The denoising results are compared with several methods and RPCA is shown to outperform the other methods. It is demonstrated that RPCA is promising in reducing the background noise in PA image reconstruction. We also aim to improve the iterative reconstruction. The variance reduced stochastic gradient descent (VR-SGD) algorithm is implemented in PAT reconstruction. A new forward projection matrix is also developed to more accurately match with the measurement data. Using different evaluation criteria, such as peak signal-to-noise ratio (PSNR), relative root-mean-square of reconstruction error (RRMSE) and line profile comparisons, the reconstructions from various iterative algorithms are compared. The advantages of VR-SGD are demonstrated on both simulation and experimental data. Our results indicate that VR-SGD in combination with the accurate projection matrix can lead to improved reconstruction in a small number of iterations. RPCA denoising and VR-SGD iterative reconstruction have been implemented in PAT. Our results show that RPCA and VR-SGD are promising approaches to improve the image reconstruction quality in PAT. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-08-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0354460 |
URI | http://hdl.handle.net/2429/62689 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2017-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_september_jiayi_cheng.pdf [ 1.8MB ]
- Metadata
- JSON: 24-1.0354460.json
- JSON-LD: 24-1.0354460-ld.json
- RDF/XML (Pretty): 24-1.0354460-rdf.xml
- RDF/JSON: 24-1.0354460-rdf.json
- Turtle: 24-1.0354460-turtle.txt
- N-Triples: 24-1.0354460-rdf-ntriples.txt
- Original Record: 24-1.0354460-source.json
- Full Text
- 24-1.0354460-fulltext.txt
- Citation
- 24-1.0354460.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0354460/manifest