Photoacoustic image improvement based on a combination of sparse coding and filtering

Ebrahim Najafzadeh; Parastoo Farnia; Saeedeh N. Lavasani; Maryam Basij; Yan Yan; Hossein Ghadiri; Alireza Ahmadian; Mohammad Mehrmohammadi

doi:10.1117/1.JBO.25.10.106001

7 October 2020 Photoacoustic image improvement based on a combination of sparse coding and filtering

Ebrahim Najafzadeh, Parastoo Farnia, Saeedeh N. Lavasani, Maryam Basij, Yan Yan, Hossein Ghadiri, Alireza Ahmadian, Mohammad Mehrmohammadi

Author Affiliations +

Journal of Biomedical Optics, Vol. 25, Issue 10, 106001 (October 2020). https://doi.org/10.1117/1.JBO.25.10.106001

Abstract

Significance: Photoacoustic imaging (PAI) has been greatly developed in a broad range of diagnostic applications. The efficiency of light to sound conversion in PAI is limited by the ubiquitous noise arising from the tissue background, leading to a low signal-to-noise ratio (SNR), and thus a poor quality of images. Frame averaging has been widely used to reduce the noise; however, it compromises the temporal resolution of PAI.

Aim: We propose an approach for photoacoustic (PA) signal denoising based on a combination of low-pass filtering and sparse coding (LPFSC).

Approach: LPFSC method is based on the fact that PA signal can be modeled as the sum of low frequency and sparse components, which allows for the reduction of noise levels using a hybrid alternating direction method of multipliers in an optimization process.

Results: LPFSC method was evaluated using in-silico and experimental phantoms. The results show a 26% improvement in the peak SNR of PA signal compared to the averaging method for in-silico data. On average, LPFSC method offers a 63% improvement in the image contrast-to-noise ratio and a 33% improvement in the structural similarity index compared to the averaging method for objects located at three different depths, ranging from 10 to 20 mm, in a porcine tissue phantom.

Conclusions: The proposed method is an effective tool for PA signal denoising, whereas it ultimately improves the quality of reconstructed images, especially at higher depths, without limiting the image acquisition speed.

1. Introduction

In the last two decades, photoacoustic imaging (PAI) as a non-invasive hybrid imaging modality has been used in a wide range of preclinical and clinical applications, such as functional brain mapping,¹^,² molecular imaging,³^,⁴ cancer diagnosis and staging,⁵^–⁸ tissue vasculature imaging,⁹^–¹² guiding interventional procedures,¹³^,¹⁴ and dental health.¹⁵^,¹⁶ PAI detects the optical absorption contrast in tissue through the conversion of light to heat and thermoelastic effect, leading to the generation of acoustic waves.¹⁷^–¹⁹ When the tissue is illuminated by short light pulses, the endogenous chromophores, such as hemoglobin, generate a photoacoustic (PA) signal due to their optical absorption.²⁰^,²¹ In this procedure, the light energy is transformed into acoustic waves, and the efficacy of this conversion is often affected by the presence of noise arising from the surrounding background.²²^,²³ Therefore, the PA signal is often mixed by background noise, including thermal-acoustic noise in the medium as well as the transducer and electronic noises.²⁴ White Gaussian noise is one of the most common models for these types of randomly distributed thermal and electronic noise.²⁵ However, there is another component of noise arising from the light attenuation phenomena, caused by scattering characteristics of the tissue.²⁶ Furthermore, other types of noise disturb the PA signals, such as the fixed-pattern noise caused by electromagnetic interference. The combination of these different types of noise in the PA signal leads to a low signal-to-noise ratio (SNR) and results in a poor quality reconstructed PA image.²⁷^–³⁰

Several previously reported studies have attempted to improve the image reconstruction algorithms to achieve noise- and artifact-free PA images. However, acquiring high-quality PA images from noisy signals requires an effective denoising technique,³¹^–³⁶ prior to utilizing reconstruction algorithms. The most commonly used technique to reduce the noise level and improve the SNR of the PA signal is frame averaging, where the SNR improvement is proportional to the square root of the number of averaged frames.³⁷ However, the signal averaging method needs to acquire multiple frames, which is time-consuming and affects the frame rate of PA imaging.²⁹ Alternatively, the adaptive filtering method without any prior knowledge requirement was proposed for low-energy pulse laser diodes PA signal enhancement.²⁸ They average fewer frames, in comparison with conventional averaging techniques, which leads to shorter acquisition time.²⁸ Notwithstanding, the requirement of frames averaging in the latter method compromises the imaging speed.

Linear time-variant filtering techniques have been widely used to improve the SNR of the PA signal. Generally, these filtering techniques suffer from the inability of low-pass and bandpass filters when signal and noise share similar frequency spectrum.³⁸ A commonly used PA signal denoising technique is the wavelet denoising method.³⁹^–⁴³ Wavelet-based denoising methods face some challenges, such as choosing an appropriate basis function, the optimum number of wavelet decomposition levels, and especially choosing an optimum threshold value.²⁸ While there are solutions for these drawbacks, these solutions are often complicated and computationally expensive.⁴⁴^,⁴⁵ Furthermore, different denoising methods based on empirical mode decomposition (EMD)⁴⁶ were proposed to improve the SNR of the signal. The EMDs decompose the PA signal into several intrinsic mode functions, which should be selected expertly and used for signal denoising. In this regard, the combination of the EMD method and mutual information (MI) was proposed to denoise PA signals.²⁹ Although this method outperforms conventional wavelet and bandpass filtering methods in the term of SNR, it works based on an inaccurate assumption that high-frequency intrinsic mode functions contain considerable noise, and low-frequency intrinsic mode functions contain the majority of useful signals. Additionally, it is not an appropriate technique in the process of real-time PAI imaging, considering EMD and MI consumed time.

Since the size of the light absorber defines the spectrum content of the PA signals, these signals are usually broadband and cannot be considered as signals with specific frequency bands.⁴⁷ On the other hand, the PA signal could be considered as a sparse signal with a sparse derivative. Therefore, one can model the PA signal as a sum of two components of low frequency and sparse.

In this study, for the first time, we proposed an approach for PA signals denoising based on a combination of low-pass filtering and total variation (TV) denoising, allowing for using a hybrid alternating direction method of multipliers (ADMM) in the optimization processes. Since the proposed method can be defined as a combination of low-pass filtering and sparse coding, we call it LPFSC.

The rest of this paper is organized as follows. Section 2 discusses the theory of TV denoising and ADMM methods, which are used in this paper. Section 3 describes the proposed denoising method, validation studies, and evaluation criteria. The experiments conducted to evaluate the performance of the proposed approach are described in Sec. 3, and the results are discussed in Sec. 4, followed by conclusions in Sec. 5.

2. Theoretical Background

2.1.

Principles of Total Variation Denoising

Recently, the TV approach that promotes the sparsity of signals in the gradient domain has attracted significant attention in signal denoising applications.⁴⁸ The goal of TV denoising technique is to efficiently estimate and recover the desired $N$ -point signal $S = {s (n)}_{n = 1}^{N}$ with the sparse or sparse-derivative representation from the measured noisy signal $x$ , which is defined as

Eq. (1)

x (n) = s (n) + w (n),

where

w (n)

is considered as additive Gaussian noise with the variance of

σ^{2}

. TV denoising could be defined as the constrained minimization problem of a non-differentiable cost function in terms of the

l_{1}

norm as below:

Eq. (2)

\arg \min_{s} {‖ D s ‖}_{1}, subject to, {‖ x - s ‖}_{2}^{2} \leq N σ^{2},

where

D

as the first-order difference matrix is of size

N \times (N - 1)

and

D s

is the first-order difference of an

N

-point signal

s (n)

.⁴⁹ With proper regularization parameter selection, Eq. (2) could be converted to the unconstraint problem as

Eq. (3)

T v d (x, λ) = \arg \min_{s} {\frac{1}{2} {‖ x - s ‖}_{2}^{2} + λ {‖ D s ‖}_{1}} .

In this optimization problem, the regularization parameter $λ$ plays a significant role and controls the degree of smoothing, so that when the $λ = 0$ , there is no smoothing and the result is the same as minimizing the sum of squares. On the other hand, increasing $λ$ assigns a higher weight to the second term of $T v d (x, λ)$ , which measures the oscillation of the desired signal $s (n)$ and makes the solution $s (n)$ piecewise. Although there are different algorithms to solve the TV denoising problem, majorization–minimization (MM) is found suitable to solve the optimization problem, which is hard to solve directly.⁵⁰

2.2.

Alternating Direction Method of Multipliers

ADMM is a simple but powerful algorithm to solve a convex optimization problem by breaking it into smaller subproblems.⁵¹^,⁵² The ADMM algorithm is designed to solve the separable convex problems of the form:

Eq. (4)

\min f (x) + g (y), subject to A x + B y = c,

where

x \in R^{n}

,

y \in R^{m}

,

A \in R^{p \times n}

, and

B \in R^{p \times m}

. The augmentation Lagrangian in Eq. (4) can be written as

Eq. (5)

L_{ρ} (x, y, λ) = f (x) + g (y) + λ^{T} (A x + B y - c) + (\frac{ρ}{2}) {‖ A x + B y - c ‖}_{2}^{2},

where

ρ

is the penalty parameter, which is considered positive, and

λ

is the Lagrangian multiplier. Equation (5) is solved using three steps:

x

-minimization and

y

-minimization, which are split into

N

separate problems, and an updating step for multiplier

λ

as follows:

Eq. (6)

x^{k + 1} ≔ \arg \min_{x} L_{ρ} (x, y^{k}, λ^{k}), y^{k + 1} ≔ \arg \min_{y} L_{ρ} (x^{k + 1}, y, λ^{k}), λ^{k + 1} ≔ λ^{k} + ρ (A x^{k + 1} + B y^{k + 1} - c) .

3. Method and Materials

3.1.

Low-Pass Filtering and Sparse Coding

To solve the problem of PA signal denoising, the PA signal is modeled as the measured noisy signal $x (n)$ , which is defined as

Eq. (7)

x (n) = s_{lowfrq} (n) + s_{sparse} (n) + w (n),

where the desired signal

s (n)

includes two main components, a low-frequency component and a sparse or a sparse-derivative component. Here,

s_{lowfrq} (n)

represents the low-frequency component of the desired signal,

s_{sparse} (n)

represents the sparse or the sparse-derivative components of the desired signal, and

w (n)

is considered as additive Gaussian noise with the variance of

σ^{2}

. Since we are looking for efficient estimation of

s_{sparse} (n)

and

s_{lowfrq} (n)

, considering Eq. (7), in the following of

s_{sparse} (n)

estimation, one can estimate

s_{lowfrq} (n)

as follows:

Eq. (8)

s_{lowfrq} (n) \approx lowpass [x (n) - {\hat{s}}_{sparse} (n)],

By replacing $s_{lowfrq} (n)$ based on Eq. (7), Eq. (8) could be written as

Eq. (9)

x (n) - w (n) - s_{sparse} (n) \approx lowpass [x (n) - {\hat{s}}_{sparse} (n)],

where

s_{sparse} (n)

and

{\hat{s}}_{sparse} (n)

are approximately equal. Therefore, Eq. (8) could be modified to

Eq. (10)

[x (n) - {\hat{s}}_{sparse} (n)] - lowpass [x (n) - {\hat{s}}_{sparse} (n)] \approx w (n),

Considering the assumption that the frequency response of the low-pass filter is approximately zero-phase, we can conclude that

Eq. (11)

highpass [x (n) - {\hat{s}}_{sparse} (n)] \approx w (n),

where high pass refers to the high-pass filter. To achieve a computationally efficient approach, a zero-phase non-causal recursive high-pass filter which is proposed in Ref. 49 was deployed in our method. Since the PA signal is sparse and has a sparse derivative, the cost function of the optimization problem contains a linear combination of two regularization parameters, which promote piecewise smooth and sparse solutions. Therefore, the denoising problem, shown in Eq. (3), can be expressed as the unconstraint minimization problem of a non-differentiable cost function in terms of

l_{1}

norm as below:

Eq. (12)

\arg \min_{s, s_{sparse}} {\frac{1}{2} {‖ highpass (x - s_{sparse}) ‖}_{2}^{2} + λ_{0} {‖ s ‖}_{1} + λ_{1} {‖ D s ‖}_{1}} .

There are many solutions, such as MM and ADMM, to solve Eq. (12). Since the denoising process and image reconstruction speed are vital to achieve real-time PA imaging, we proposed to use the hybrid consensus ADMM method⁵³ to achieve a linear convergence and accelerating the conventional ADMM.

3.2.

Validation Studies

3.2.1.

In-silico study

A simulation study generated by the $k$ -wave toolbox in MATLAB^® (Mathworks, Massachusetts) was performed to evaluate the performance of LPFSC in PA signal denoising.⁵⁴ The initial pressure distribution is given by a $512 \times 512 pixel$ image representative of a vascular structure. A 10-mm square grid was created, and the circular array detectors with 9-mm diameter and 60 elements, evenly spaced, were located around the region of interest (ROI) to receive the propagated PA wave from the object. We considered the center of each sensor as a point source. The sound speed was considered to be $1500 m / s$ . The corresponding time array has 1019 data points that are 9.259 ns apart from each other (108-MHz sampling frequency). The input size of the PA signal was assigned $1019 \times 60$ , and the reconstructed images have $128 \times 128 pixels$ .

3.2.2.

Experimental PA data acquisition setup

To further evaluation of the proposed PA signal denoising method and its effects on the quality of reconstructed PA images, phantom experiments were performed. The designed phantoms and the imaging setups are shown in Fig. 1. The first phantom contains two light-absorbing filaments with a diameter of $150 μ m$ that were placed 1 mm apart from each other inside a water tank [Fig. 2(b)]. An Nd: YAG/OPO nanosecond pulsed laser (Phocus core system, OPOTEK Inc., Carlsbad, California) with the pulse repetition rate of 10 Hz at wavelengths of 680 nm was used to illuminate the phantom. An ultrasound scanner (US) (Vantage 128TM, Verasonics Inc., Kirkland, Washington) with a 128 elements linear array transducer (L114v, Verasonics, Inc., Kirkland, Washington) was used to receive the propagated PA RF data [Fig. 1(a)].

Fig. 1

Schematic of (a) the experimental setup used for the PA imaging of (b) two light-absorbing filaments that were placed inside a water tank and (c) cross-section view of blood-filled tubes embedded in a porcine tissue phantom.

Fig. 2

The denoising results of simulated PA signal in four different levels of noise. Noisy signals are shown in red and denoised signals are depicted in blue. (a) The original simulated PA signal detected on a detector #50, (b) SNR: 10 dB (c) SNR: 5 dB, (d) SNR: $- 5 dB$ , and (e) SNR: $- 10 dB$ . The horizontal and vertical axes indicate the time coordinate and amplitude of the recovered signal in different noise levels, respectively.

The second phantom contains three polytetrafluoroethylene tubes (1-mm diameter), filled with human blood and embedded within a 30-mm-thick porcine tissue background. Blood-filled tubes were placed inside the porcine tissue background at different depths from about 10 to 20 mm with 5-mm increments. PA acquisition was performed with the laser energy of $3 mJ / pulse$ at the wavelength of 680 nm. A fiber bundle with a diameter of 20 mm used for guiding the laser light to the tissue. A 64-element phased-array US endoscopic transducer, with an active aperture of 9-mm long, was coupled to the phantom to acquire acoustic signals and provide high-resolution sector images [Fig. 1(c)].

3.3.

PA Signal Evaluation and Analysis

To evaluate the proposed approach, evaluation criteria such as peak signal-to-noise (PSNR) for signals, structural similarity index (SSIM), and contrast-to-noise ratio (CNR) for estimated images were used.

The PSNR as a common criterion to measure the quality of signal denoising based on the maximum possible value in the signal and mean square differences between denoised and reference signals is expressed in term of the logarithmic decibel scale (dB) as below:

PSNR = 20 \log_{10} (\frac{S_{\max}}{\sqrt{MSE}}), MSE = \frac{1}{M} \sum_{m = 0}^{M - 1} {[S_{original} (m, 1) - S_{denoised} (m, 1)]}^{2},

where MSE is defined as a mean-square-error,

S_{original}

and

S_{denoised}

are original and denoised signals in size of

M \times 1

, respectively, and

S_{\max}

is maximum possible value in signals.

We created PA images that represent an optical absorption distribution map of the targets via the conventional delay-and-sum (DAS) approach as the most commonly used reconstruction method in the PAI area.⁵⁵ For reconstructed images, the SSIM (in a scale of 0 to 1) as one of the most common criteria for image quality assessment and for evaluating the similarity of images (i.e., reference image and reconstructed images) is defined⁵⁶ as

SSIM = \frac{(2 μ_{Original} μ_{estimated} + c_{1}) (2 σ_{original, estimated} + c_{2})}{\sqrt{(μ_{original}^{2} + μ_{estimated}^{2} + c_{1}) (σ_{original}^{2} + σ_{estimated}^{2} + c_{2})}},

where

μ_{original}

and

μ_{estimated}

are the mean of the original and estimated images, respectively, and also

σ_{original}^{2}

and

σ_{estimated}^{2}

are the variances of the original and estimated images, respectively. It is worth to mention that

σ_{original, estimated}

is the covariance between the original and estimated image. The values of

c_{1}

and

c_{2}

are considered as constant values to avoid instability when the sum square of means or variances are very close to zero.

Finally, the CNR is a common criterion to determine image quality, especially in denoising processes with below definition:

CNR = 20 \log (\frac{| S_{i} - S_{O} |}{σ_{O}}),

where

S_{i}

and

S_{O}

are the average intensity inside and outside of the objects, respectively. The

σ_{O}

represents the standard deviation of the background. The background was defined as the pixels located inside the green dashed rectangular region selected in each set of PA images. For phantom studies, we considered the average of all frames as a ground truth (reference image) in each experiment.

4. Results and Discussion

To assess the proposed PA signal denoising method, we validated our method on numerical vessel phantom and experimental data of phantoms. The simulated PA signal, which was generated by the $k$ -wave is shown in Fig. 2(a). Four different levels of additive Gaussian white noise with SNR levels 10, 5, $- 5$ , and $- 10 dB$ were added to the original simulated clean PA signal, and results of the proposed denoising method for the simulated noisy signals are shown in Figs. 2(b)–2(e). When the noise level was increased (for example, in SNR level $- 10$ ), the PA signal peak was almost buried in added noise. As we have shown in Fig. 2, the LPFSC could suppress the noise and reconstruct the peak of the original signal, indicating the ability of the TV approach to recover the sparse or sparse-derivative signals. In the worst case (SNR level $- 10 dB$ ), the difference between the original signal peak and the denoised signal peak is about 4%. This variance considering increasing the peak of noisy signal about 90% in comparison with the original signal is negligible. TV-based denoising is the most appropriate method for piecewise constant signals and preserves sharp edges in the underlying signal without requiring any step-size parameter as the amount of peak for the denoised signal.⁵⁷ It is worth mentioning that the proposed method could recover denoised signal until the SNR decreased to $- 15 dB$ .

For further evaluation, the LPFSC signal denoising method was compared with two well-known and widely used approaches of averaging and wavelet-based denoising methods. As shown in Fig. 3, the PSNR of the LPFSC method was compared to the frame averaging (using 20 frames) and wavelet-based denoising methods, at different noise levels. For PA signal denoising by wavelet method, a commonly used Symlet 6 wavelet with six-level of decomposition and Stein’s Unbiased Risk Estimation threshold (SURE threshold) were selected.²⁹

Fig. 3

Comparison of three different denoising methods: wavelet-based signal denoising, the averaging method using 20 frames, and LPFSC in terms of PSNR improvements at five different noise levels. The results clearly demonstrate the superior performance of LPFSC compared to wavelet and averaging.

The quantitative results obtained with the simulations show that the LPFSC compensates the low SNR of PA signal and outperforms the competing wavelet and averaging denoising methods in terms of PSNR by 24% and 26%, respectively, across all simulated noise levels.

The reconstructed images of denoised simulated PA signals with considering SNR of $- 10 dB$ for three methods of averaging, wavelet-based denoising, and LPFSC are shown in Fig. 4. The wavelet method has not been successful in recovering the small size vessels, which are corrupted by the noise. In addition, the averaging method using 20 frames cannot fully clean the signal due to presences of the coherent noise; however, our method using one frame could fully recover the original image and its details.

Fig. 4

(a) Illustration of the phantom used for in-silico studies. Reconstructed images of: (b) simulated PA phantom, (c) noisy image when SNR is considered $- 10 dB$ , (d) PA image using the LPFSC denoising, (e) PA image using the wavelet denoised, and (f) PA image using the averaging of 20 frames for noise reduction.

For the quantitative evaluation of the reconstructed images of denoised simulated PA signals, two SSIM and CNR criteria were used. Comparison between the results of three different methods, including wavelet-based signal denoising, averaging, and LPFSC in terms of SSIM and CNR with different levels of noise are shown in Table 1. On average (for different noise levels), the proposed denoising approach offers better CNR up to about 28% and 30%, and higher SSIM of about 23% and 24% in comparison with averaging and wavelet denoising methods, respectively. With increasing SNR levels from 5 to 20 dB and the reduction of noise level, the performance of all methods are improved. Although, for the lower SNR such as 5 dB, LPFSC outperforms wavelet and averaging denoising method by improvement about 58% and in terms of SSIM.

Table 1

Performance of three different signal denoising methods: wavelet-based, averaging method, and LPFSC across different levels of noise.

Level of SNR (dB)	Wavelet		Averaging		LPFSC
Level of SNR (dB)	SSIM	CNR (dB)	SSIM	CNR (dB)	SSIM	CNR (dB)
5	0.59	65.01	0.58	67.23	0.93	87.61
10	0.78	72.53	0.79	72.98	0.95	94.53
15	0.86	75.15	0.87	76.23	0.98	96.70
20	0.86	76.66	0.91	77.12	0.99	98.21
$Mean \pm std$	$0.77 \pm 0.13$	$72.33 \pm 5.17$	$0.78 \pm 0.14$	$73.39 \pm 4.47$	$0.96 \pm 0.03$	$94.26 \pm 4.68$

During experimental studies, we evaluated our method in a set of phantom studies in which two small-sized absorbers were placed inside a water tank, illuminated from the side and images from the top using a US transducer operating at frequencies between 4 and 11 MHz. The PA signal denoising results for the acquired signals of different detectors from the 20th frames of phantom data are shown in Fig. 5. The experimental results of PA signal denoising show significant improvement for the PSNR of PA signal of wires phantom about 35% by recovering peaks of original signal and reduction of noise.

Fig. 5

(a) The PA signal denoising of experimental phantom via LPFSC and (b) the original reconstructed image of one frame. The reconstructed images of phantom via three different signal denoising methods: (c) the LPFSC method (one frame), (d) the wavelet denoising method (one frame), and (e) the averaging all 20 frames. The background was defined as the pixels located inside the green dashed rectangular region.

Additionally, reconstructed images of denoised PA signals through three different methods are shown in Fig. 5 and are compared in Table 2, in terms of SSIM, CNR, axial and lateral full width half maximum (FWHM) for experimental phantom data. The selection of FWHM is to demonstrate the geometrical accuracy of reconstructed PA images. The experimental results of PA signal denoising prove that the LPFSC beats the performance of averaging denoising method (using all 20 frames) as the reference in the term of SSIM and CNR criteria by the 17% and 38%. The LPFSC method provides a closer to real size reconstructed PA image with a mean of axial FWHM of $0.77 \pm 0.12 mm$ and lateral FWHM of $1.05 \pm 0.04 mm$ , for two objects in experimental data. In comparison with the reconstructed image of one frame without any signal denoising with axial FWHM of $0.85 \pm 0.18$ and lateral FWHM of $1.08 \pm 0.05 mm$ , our method could preserve lateral FWHM and improve axial FWHM about 10%. Whereas the wavelet denoising method did not improve the lateral and axial FWHM, and the averaging method was not effective for axial FWHM. It is worth to mention that both of these methods improve CNR of the reconstructed image in comparison with the reconstructed image of one frame with CNR 19.32 dB. Therefore, one of the main improvements gained by the proposed LPFSC method is having reduced geometrical distortion as well as a high contrast at the same time. More importantly, the LPFSC method offers using one frame in comparison with the averaging method, which requires 20 frames. Using the hybrid ADMM, which stands out as efficient and easily implementable, leads to fast convergence of our method.

Table 2

Comparison of three methods of wavelet, averaging, and LPFSC method in terms of SSIM, CNR, and lateral and axial FWHM for the reconstructed image of the denoised signal.

Criteria	Wavelet	Averaging	LPFSC
SSIM	0.79	0.81	0.95
CNR (dB)	23.72	25.36	35.17
Lateral FWHM (mm) (mean + std)	$1.20 \pm 0.06$	$1.11 \pm 0.09$	$1.05 \pm 0.04$
Axial FWHM (mm) (mean + std)	$1.12 \pm 0.16$	$0.98 \pm 0.06$	$0.77 \pm 0.12$

Finally, we evaluated the performance of the proposed denoising method with data acquired by an endoscopic probe and from the porcine tissue phantom. Since there was no significant difference between the results of frame averaging and wavelet method in in-silico study, the reconstructed PA images of denoised signal with frame averaging and LPFSC methods are compared in Fig. 6. The images form by averaging of 10 frames were used as a ground truth or the reference image. All constructed images from denoised signal using LPFSC and traditional frame averaging with the same number of frames ( $n = 1$ to 7) were compared to the ground-truth reference to calculate the SSIM. The number of averaged frames varied between one and seven. First and second rows related to the results of the averaging method and our proposed method, respectively. The walls of blood-filled tubes are generating strong PA signals that are visually presented as two parallel bright lines in PA images. To calculate CNR, an ROI located within the tube [as shown in Fig. 6(e)] was used to measure the signal from the target. To further quantify the image enhancements, CNR and SSIM of frame averaging and LPFSC using five frames are compared in Fig. 6(h). The use of five frames averaging in conjunction with LPFSC markedly reduced the background noise, which improved the mean of CNR and SSIM by 32% and 31%, respectively.

Fig. 6

Reconstructed PA images of denoised signals of a porcine/blood inclusion phantom. (a)–(c) Frame averaging with 1, 5, and 7 averaged frames, respectively. (d)–(f) LPFSC with 1, 5, and 7 averaged frames, respectively. (g) Averaging 10 frames used as a reference. (h) The comparison of CNR and SSIM of two methods using five frames for inclusions located at different depths. The black arrows indicate objects located at depths 1, 2, and 3. Depths 1, 2, and 3 are about 20, 15, and 10 mm, respectively. The background was defined as the pixels located inside the green dashed rectangular region. The ROI for the target (object) is indicated with black dashed rectangular boxes in (e).

Evaluation of CNR and SSIM for different numbers of averaged frames and at different imaging depths are shown in Fig. 7. As anticipated, CNR and SSIM parameters are improved with increasing the number of averaged frames. However, the proposed denoising method provides higher CNR and SSIM using the same number of averaged frames for objects placed at different depths. The mean of CNR using LPFSC with one frame was improved 63% in comparison with averaging for the objects evaluated at three different depths. A closer look at the results reveal that the mean CNR using the proposed denoising method using three averaged frames is higher than averaging only and using seven averaged frames. This comparison clearly shows the potential of using the proposed method to enhance the quality of PA images with a smaller number of averaged frames, which can lead to a preserving image quality at higher imaging speed. Our results indicated the LPFSC with three frames averaging is capable of detecting the objects located at depth 3 with SSIM of 0.91, which shows 90% improvement compared to the standard averaging method. The average measured SSIM for objects located at three different depths was calculated as $0.88 \pm 0.15$ when LPFSC is utilized. This result indicates 33% improvement compared to averaging alone.

Fig. 7

Porcine tissue results: comparison between the LPFSC with averaging using different number of averaged frames ranging from 1 to 7 in terms of CNR and SSIM: (a) depth 1, (b) depth 2, and (c) depth 3. CNR and SSIM parameters are evaluated and the results indicate the superior performance of LPFSC compared to averaging only.

In this study, a range of cutoff frequencies was experimentally tested on data sets. The accepted range of cutoff frequencies was concluded to be between 0 and $0.3 \frac{cycles}{sample}$ for different data sets, experimentally. For in-silico and first experimental phantom studies, the optimal value of cutoff frequency was selected $0.1 \frac{cycles}{sample}$ and for porcine tissue data, it was selected $0.15 \frac{cycles}{sample}$ . Also, we consider equal weight for the TV and LPF; therefore, the lambda value was selected equal to be 1, experimentally. Changing lambda within the range of 0.8 to 1.2 showed acceptable results for PA signal denoising. For lambda smaller than 0.8, the amplitude of denoised signal is decreased compared to the original signal. The lambda bigger than 1.2 downgraded the performance of LPFSC. The performance of LPFSC could also be affected by the adjustment of low-pass filter cutoff frequency. However, small deviations from the ideal cutoff frequency up to 20% could be compensated by the sparse denoising part of the proposed algorithm.

5. Conclusion

The efficiency of PA imaging is routinely limited by the presence of background noise and suffering from low SNR, which resulted in the poor quality of the reconstructed images. Since the PA signals can be modeled as a sum of two low frequency and sparse components, we proposed a denoising approach that simultaneously estimates a low-pass and a sparse signal from an acquired noisy signal based on TV optimization approach, and using hybrid ADMM. Both in-silico and experimental work on tissue mimics were used to evaluate the performance of the proposed technique. The results demonstrated that the LPFSC method possess a superior performance to compensate the low SNR PA signals and offered a better CNR and SSIM for the reconstructed images compared to the frame averaging method. This comparison clearly shows the potential of using our proposed method to enhance the quality of PA images while maintaining high-speed imaging which is an essential need in many applications of PA imaging. In other words, the frame rate of PA images, which is always a challenge in real-time PA imaging can be significantly improved.

Disclosures

The authors declare that they have no conflicts of interest.

References

1.

X. Wang et al., “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol., 21 (7), 803 (2003). https://doi.org/10.1038/nbt839 NABIF9 1087-0156 Google Scholar

2.

M. Nasiriavanaki et al., “High-resolution photoacoustic tomography of resting-state functional connectivity in the mouse brain,” Proc. Natl. Acad. Sci. U. S. A., 111 (1), 21 –26 (2014). https://doi.org/10.1073/pnas.1311868111 Google Scholar

3.

J.-W. Kim et al., “Golden carbon nanotubes as multimodal photoacoustic and photothermal high-contrast molecular agents,” Nat. Nanotechnol., 4 (10), 688 (2009). https://doi.org/10.1038/nnano.2009.231 NNAABX 1748-3387 Google Scholar

4.

A. De La Zerda et al., “Carbon nanotubes as photoacoustic molecular imaging agents in living mice,” Nat. Nanotechnol., 3 (9), 557 (2008). https://doi.org/10.1038/nnano.2008.231 NNAABX 1748-3387 Google Scholar

5.

M. Basij et al., “Miniaturized phased-array ultrasound and photoacoustic endoscopic imaging system,” Photoacoustics, 15 100139 (2019). https://doi.org/10.1016/j.pacs.2019.100139 Google Scholar

6.

S. S. Alshahrani et al., “All-reflective ring illumination system for photoacoustic tomography,” J. Biomed. Opt., 24 (4), 046004 (2019). https://doi.org/10.1117/1.JBO.24.4.046004 JBOPFO 1083-3668 Google Scholar

7.

M. Mehrmohammadi et al., “Photoacoustic imaging for cancer detection and staging,” Curr. Mol. Imaging, 2 (1), 89 –105 (2013). https://doi.org/10.2174/2211555211302010010 Google Scholar

8.

M. Basij et al., “Development of an ultrasound and photoacoustic endoscopy system for imaging of gynecological disorders,” in IEEE Int. Ultrasonics Symp., 1 –4 (2018). Google Scholar

9.

E. Najafzadeh et al., “Application of multi-wavelength technique for photoacoustic imaging to delineate tumor margins during maximum-safe resection of glioma: a preliminary simulation study,” J. Clin. Neurosci., 70 242 –246 (2019). https://doi.org/10.1016/j.jocn.2019.08.040 Google Scholar

10.

J. G. Laufer et al., “In vivo preclinical photoacoustic imaging of tumor vasculature development and therapy,” J. Biomed. Opt., 17 (5), 056016 (2012). https://doi.org/10.1117/1.JBO.17.5.056016 JBOPFO 1083-3668 Google Scholar

11.

K. M. Kempski et al., “In vivo photoacoustic imaging of major blood vessels in the pancreas and liver during surgery,” J. Biomed. Opt., 24 (12), 121901 (2019). https://doi.org/10.1117/1.JBO.24.12.121905 JBOPFO 1083-3668 Google Scholar

12.

E. Najafzadeh et al., “Evaluation of multi-wavelengths LED-based photoacoustic imaging for maximum safe resection of glioma: a proof of concept study,” Int. J. Comput. Assist. Radiol. Surg., 15 1053 –1062 (2020). https://doi.org/10.1007/s11548-020-02191-2 Google Scholar

13.

Y. Yan et al., “Photoacoustic imaging for image-guided endovenous laser ablation procedures,” Sci. Rep., 9 (1), 1 –10 (2019). https://doi.org/10.1038/s41598-018-37588-2 SRCEC3 2045-2322 Google Scholar

14.

C. Kim et al., “Handheld array-based photoacoustic probe for guiding needle biopsy of sentinel lymph nodes,” J. Biomed. Opt., 15 (4), 046010 (2010). https://doi.org/10.1117/1.3469829 JBOPFO 1083-3668 Google Scholar

15.

S. Arabpou et al., “Detection of early stages dental caries using photoacoustic signals: the simulation study,” Front. Biomed. Technol., 6 35 –40 (2019). https://doi.org/10.18502/fbt.v6i1.1101 Google Scholar

16.

C. Moore et al., “Photoacoustic imaging for monitoring periodontal health: a first human study,” Photoacoustics, 12 67 –74 (2018). https://doi.org/10.1016/j.pacs.2018.10.005 Google Scholar

17.

L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science, 335 (6075), 1458 –1462 (2012). https://doi.org/10.1126/science.1216210 SCIEAS 0036-8075 Google Scholar

18.

P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, 1 (4), 602 –631 (2011). https://doi.org/10.1098/rsfs.2011.0028 Google Scholar

19.

L. V. Wang and J. Yao, “A practical guide to photoacoustic tomography in the life sciences,” Nat. Methods, 13 (8), 627 (2016). https://doi.org/10.1038/nmeth.3925 1548-7091 Google Scholar

20.

L. V. Wang and L. Gao, “Photoacoustic microscopy and computed tomography: from bench to bedside,” Annu. Rev. Biomed. Eng., 16 155 –185 (2014). https://doi.org/10.1146/annurev-bioeng-071813-104553 ARBEF7 1523-9829 Google Scholar

21.

M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 (4), 041101 (2006). https://doi.org/10.1063/1.2195024 RSINAK 0034-6748 Google Scholar

22.

L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron., 14 (1), 171 –179 (2008). https://doi.org/10.1109/JSTQE.2007.913398 IJSQEN 1077-260X Google Scholar

23.

A. M. Winkler, K. I. Maslov and L. V. Wang, “Noise-equivalent sensitivity of photoacoustics,” J. Biomed. Opt., 18 (9), 097003 (2013). https://doi.org/10.1117/1.JBO.18.9.097003 JBOPFO 1083-3668 Google Scholar

24.

J. Yao and L. V. Wang, “Sensitivity of photoacoustic microscopy,” Photoacoustics, 2 (2), 87 –101 (2014). https://doi.org/10.1016/j.pacs.2014.04.002 Google Scholar

25.

B. Stephanian et al., “Additive noise models for photoacoustic spatial coherence theory,” Biomed. Opt. Express, 9 (11), 5566 –5582 (2018). https://doi.org/10.1364/BOE.9.005566 BOEICL 2156-7085 Google Scholar

26.

S. Mahmoodkalayeh et al., “Low temperature-mediated enhancement of photoacoustic imaging depth,” Sci. Rep., 8 (1), 4873 (2018). https://doi.org/10.1038/s41598-018-22898-2 SRCEC3 2045-2322 Google Scholar

27.

S. Telenkov and A. Mandelis, “Signal-to-noise analysis of biomedical photoacoustic measurements in time and frequency domains,” Rev. Sci. Instrum., 81 (12), 124901 (2010). https://doi.org/10.1063/1.3505113 RSINAK 0034-6748 Google Scholar

28.

R. Manwar et al., “Photoacoustic signal enhancement: towards utilization of low energy laser diodes in real-time photoacoustic imaging,” Sensors, 18 (10), 3498 (2018). https://doi.org/10.3390/s18103498 SNSRES 0746-9462 Google Scholar

29.

M. Zhou et al., “A noise reduction method for photoacoustic imaging in vivo based on EMD and conditional mutual information,” IEEE Photonics J., 11 (1), 1 –10 (2019). https://doi.org/10.1109/JPHOT.2019.2897748 Google Scholar

30.

P. Farnia et al., “Dictionary learning technique enhances signal in LED-based photoacoustic imaging,” Biomed. Opt. Express, 11 (5), 2533 –2547 (2020). https://doi.org/10.1364/BOE.387364 BOEICL 2156-7085 Google Scholar

31.

M. Erfanzadeh and Q. Zhu, “Photoacoustic imaging with low-cost sources; a review,” Photoacoustics, 14 1 –11 (2019). https://doi.org/10.1016/j.pacs.2019.01.004 Google Scholar

32.

C. Huang et al., “Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imaging, 32 (6), 1097 –1110 (2013). https://doi.org/10.1109/TMI.2013.2254496 ITMID4 0278-0062 Google Scholar

33.

S. Antholzer, M. Haltmeier and J. Schwab, “Deep learning for photoacoustic tomography from sparse data,” Inverse Prob. Sci. Eng., 27 (7), 987 –1005 (2019). https://doi.org/10.1080/17415977.2018.1518444 Google Scholar

34.

N. Davoudi, X. L. Deán-Ben and D. Razansky, “Deep learning optoacoustic tomography with sparse data,” Nat. Mach. Intell., 1 453 –460 (2019). https://doi.org/10.1038/s42256-019-0095-3 Google Scholar

35.

E. M. A. Anas et al., “Enabling fast and high quality LED photoacoustic imaging: a recurrent neural networks based approach,” Biomed. Opt. Express, 9 (8), 3852 –3866 (2018). https://doi.org/10.1364/BOE.9.003852 BOEICL 2156-7085 Google Scholar

36.

P. Farnia et al., “High-quality photoacoustic image reconstruction based on deep convolutional neural network: towards intra-operative photoacoustic imaging,” Biomed. Phys. Eng. Express, 6 1 –10 (2020). https://doi.org/10.1088/2057-1976/ab9a10 Google Scholar

37.

J. Li et al., “A review of signal enhancement and noise reduction techniques for tunable diode laser absorption spectroscopy,” Appl. Spectrosc. Rev., 49 (8), 666 –691 (2014). https://doi.org/10.1080/05704928.2014.903376 APSRBB 0570-4928 Google Scholar

38.

J. Sun et al., “Photoacoustic wavefront shaping with high signal to noise ratio for light focusing through scattering media,” Sci. Rep., 9 (1), 1 –10 (2019). https://doi.org/10.1038/s41598-019-40919-6 SRCEC3 2045-2322 Google Scholar

39.

S. H. Holan and J. A. Viator, “Automated wavelet denoising of photoacoustic signals for circulating melanoma cell detection and burn image reconstruction,” Phys. Med. Biol., 53 (12), N227 (2008). https://doi.org/10.1088/0031-9155/53/12/N01 PHMBA7 0031-9155 Google Scholar

40.

S. H. Holan and J. A. Viator, “Automated wavelet denoising of photoacoustic signals for burn-depth image reconstruction,” Proc. SPIE, 6437 643719 (2007). https://doi.org/10.1117/12.701228 PSISDG 0277-786X Google Scholar

41.

M. Zhou et al., “Wavelet de-noising method with adaptive threshold selection for photoacoustic tomography,” in 40th Annu. Int. Conf. IEEE Eng. Med. and Biol. Soc., 4796 –4799 (2018). Google Scholar

42.

S. Tzoumas et al., “Spatiospectral denoising framework for multispectral optoacoustic imaging based on sparse signal representation,” Med. Phys., 41 (11), 113301 (2014). https://doi.org/10.1118/1.4893530 MPHYA6 0094-2405 Google Scholar

43.

E. Najafzadeh et al., “Development of an algorithm for reducing scatter distribution in fluorescence reflectance images using Wiener filter in wavelet domain,” in World Congr. Med. Phys. and Biomed. Eng., 1869 –1872 (2009). Google Scholar

44.

C. B. Smith, S. Agaian and D. Akopian, “A wavelet-denoising approach using polynomial threshold operators,” IEEE Signal Process. Lett., 15 906 –909 (2008). https://doi.org/10.1109/LSP.2008.2001815 IESPEJ 1070-9908 Google Scholar

45.

J. Xu et al., “Adaptive wavelet threshold denoising method for machinery sound based on improved fruit fly optimization algorithm,” Appl. Sci., 6 (7), 199 (2016). https://doi.org/10.3390/app6070199 Google Scholar

46.

Y. Lei et al., “A review on empirical mode decomposition in fault diagnosis of rotating machinery,” Mech. Syst. Signal Process., 35 (1–2), 108 –126 (2013). https://doi.org/10.1016/j.ymssp.2012.09.015 Google Scholar

47.

C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol., 54 (19), R59 (2009). https://doi.org/10.1088/0031-9155/54/19/R01 PHMBA7 0031-9155 Google Scholar

48.

S. Vujović et al., “Sparse analyzer tool for biomedical signals,” Sensors, 20 (9), 2602 (2020). https://doi.org/10.3390/s20092602 SNSRES 0746-9462 Google Scholar

49.

I. W. Selesnick et al., “Simultaneous low-pass filtering and total variation denoising,” IEEE Trans. Signal Process., 62 (5), 1109 –1124 (2014). https://doi.org/10.1109/TSP.2014.2298836 ITPRED 1053-587X Google Scholar

50.

M. A. Figueiredo, J. M. Bioucas-Dias and R. D. Nowak, “Majorization–minimization algorithms for wavelet-based image restoration,” IEEE Trans. Image Process., 16 (12), 2980 –2991 (2007). https://doi.org/10.1109/TIP.2007.909318 IIPRE4 1057-7149 Google Scholar

51.

B. Wahlberg et al., “An ADMM algorithm for a class of total variation regularized estimation problems,” IFAC Proc., 45 (16), 83 –88 (2012). https://doi.org/10.3182/20120711-3-BE-2027.00310 IPSEET Google Scholar

52.

J. Eckstein and W. Yao, “Augmented Lagrangian and alternating direction methods for convex optimization: a tutorial and some illustrative computational results,” RUTCOR Res. Rep., 32 (3), 44 (2012). Google Scholar

53.

M. Ma, A. N. Nikolakopoulos and G. B. Giannakis, “Hybrid ADMM: a unifying and fast approach to decentralized optimization,” EURASIP J. Adv. Signal Process., 2018 (1), 73 (2018). https://doi.org/10.1186/s13634-018-0589-x Google Scholar

54.

B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 (2010). https://doi.org/10.1117/1.3360308 JBOPFO 1083-3668 Google Scholar

55.

S. Jeon et al., “Real-time delay-multiply-and-sum beamforming with coherence factor for in vivo clinical photoacoustic imaging of humans,” Photoacoustics, 15 100136 (2019). https://doi.org/10.1016/j.pacs.2019.100136 Google Scholar

56.

Z. Wang et al., “Image quality assessment: from error visibility to structural similarity,” 13 (4), 600 –612 (2004). Google Scholar

57.

I. Selesnick, “Total variation denoising (an MM algorithm),” NYU Polytech. Sch. Eng. Lect. Notes, 32 1 –13 (2012). Google Scholar

Biography

Ebrahim Najafzadeh received his PhD in medical physics from Tehran University of Medical Science (TUMS), Tehran, Iran, in 2020. He has completed his undergraduate and graduate program on radiology and medical physics at Iran University of Medical Science and TUMS, respectively. His research interests are photoacoustic, ultrasound and optical imaging, and signal and image processing. His current research focuses on design of preclinical and clinical applications for LED-based photoacoustic imaging.

Parastoo Farnia received her MSc and PhD degrees in biomedical engineering from Tehran University of Medical Science, Iran, in 2013 and 2020. Her research is focused on multimodal registration like ultrasound-MRI and photoacoustic-MRI. Her research interests include multimodal image registration, sparse representation, ultrasound, and photoacoustic imaging.

Saeedeh N. Lavasani received her BA degree in biomedical engineering from the Islamic Azad University in 2011 and an MSc degree in biomedical engineering from Tehran University of Medical Science in 2014. She completed her PhD in biomedical engineering at Shahid Beheshti University of Medical Science in 2020. Her research interests include image registration, image segmentation, and photoacoustic imaging.

Maryam Basij received her BS and MS degrees in biomedical engineering from Isfahan University, Iran in 2008 and 2012, respectively. Currently, she is pursuing her PhD program in biomedical engineering with Functional and Molecular Ultrasound Research Laboratory, Wayne State University, Michigan, USA, focusing on the development of combined medical imaging technologies for advanced screening, diagnostics, and therapeutics. Besides, she worked on various biomedical images such as IVUS, breast MRI, karyotype images, nanotechnology images, etc.

Yan Yan is a PhD student at Wayne State University. His areas of interest are medical imaging, object detection, and data mining. He had strong background in computer science pattern recognition and computer graphics. He holds a bachelor’s degree from Chongqing University in computer science and telecommunication.

Hossein Ghadiri received his PhD in medical physics from Iran University of Medical Sciences, Tehran, Iran, in 2014. Since then, he had been a faculty member of medical physics and biomedical engineering at Tehran University of Medical Sciences. He is the head of Advanced X-ray Imaging Systems Lab (AXIS) with main research areas including advanced x-ray and micro-CT image analysis.

Alireza Ahmadian is a professor at Tehran University of Medical Sciences, Department of Biomedical Systems & Biophysics. He received his PhD in biomedical image processing from Imperial College, London, followed by a postdoctoral position at Kings College London. Currently, he is the head of school of advanced technologies in medicine and also the Research Center for Biomedical Technologies and Robotics. He is the founder of the knowledge based company, Parsiss, working in image-guided surgery systems.

Mohammad Mehrmohammadi is an associate professor of biomedical/electrical and computer engineering/obstetrics and gynecology at Wayne State University (Michigan, USA). He received his PhD in biomedical engineering from the University of Texas at Austin. Prior to joining WSU, Mohammad worked at Mayo Clinic College of Medicine as a senior research fellow. His research interests include medical ultrasound, photoacoustic imaging, nanobiotechnology, and tissue elastography for various applications including early cancer detection and enhancing maternal, and fetal care.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Ebrahim Najafzadeh, Parastoo Farnia, Saeedeh N. Lavasani, Maryam Basij, Yan Yan, Hossein Ghadiri, Alireza Ahmadian, and Mohammad Mehrmohammadi "Photoacoustic image improvement based on a combination of sparse coding and filtering," Journal of Biomedical Optics 25(10), 106001 (7 October 2020). https://doi.org/10.1117/1.JBO.25.10.106001

Received: 2 June 2020; Accepted: 16 September 2020; Published: 7 October 2020

Access the abstract

JOURNAL ARTICLE
16 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 13 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Denoising

Signal to noise ratio

Interference (communication)

Image quality

Photoacoustic spectroscopy

Wavelets

Image compression

1.

Introduction

2.

Theoretical Background

2.1.

Principles of Total Variation Denoising

Eq. (1)

Eq. (2)

Eq. (3)

2.2.

Alternating Direction Method of Multipliers

Eq. (4)

Eq. (5)

Eq. (6)

3.

Method and Materials

3.1.

Low-Pass Filtering and Sparse Coding

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

3.2.

Validation Studies

3.2.1.

In-silico study

3.2.2.

Experimental PA data acquisition setup

Fig. 1

Fig. 2

3.3.

PA Signal Evaluation and Analysis

4.

Results and Discussion

Fig. 3

Fig. 4

Table 1

Fig. 5

Table 2

Fig. 6

Fig. 7

5.

Conclusion

Disclosures

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years