1.

Introduction

Photoacoustic/optoacoustic imaging (PAI) is an emerging medical imaging technique combining optical contrast and spatial resolution of ultrasound (US).^1,2 PAI is based on the photoacoustic (PA) effect, which provides structural, functional, and potentially the molecular information of tissue.^3,4 In PAI, a laser pulse irradiates the sample/tissue, resulting in local temperature rise, and due to thermoelastic expansion, pressure waves (in the form of ultrasound waves) are generated.⁵ The ultrasound waves (also known as photoacoustic wave) travel within the medium and then are recorded using US transducer. These PA waves are used to obtain the optical absorption map of the inside of the tissue.⁶ In circular geometry photoacoustic tomography (PAT), a single-element/ring array US transducer can be used to acquire the PA signals around the tissue in full 360 deg.^3,7 Recently, low-cost PAT systems are extensively being investigated where pulsed laser diodes are used to make the PA systems more compact, portable, and affordable.^8,9 However, such an imaging system based on a single-element/ring array ultrasound transducer is difficult to translate into clinical applications.¹⁰ Therefore, recently dual-modal clinical ultrasound and photoacoustic imaging have been reported for better clinical translation.^11,12

Clinical ultrasound systems work with linear, convex, or phased array transducers. Here, in this work we focus on one of the clinical applications of sentinel lymph node (SLN) imaging, where a linear array ultrasound probe is usually used.¹³ Sentinel lymph node biopsy (SLNB) is a standard clinical procedure done in breast cancer staging. SLNB replaced axillary lymph node dissection (ALND), where the lymph nodes around the tumur are removed. The removal of nodes is necessary to curb metastasis of cancerous cells. ALND is an unwanted procedure in node-negative patients and it leads to side effects, such as lymphedema, arm weakness, and infections of breast. Hence, in SLNB, biopsy samples of the lymph nodes to which the mammary glands drain first are examined before removing the nodes.^14,15 SLNs are identified by intradermal injection of dyes (such as methylene blue).¹⁶ As the injected methylene blue dye gives a strong PA signal, feasibility of PA imaging for SLN identification has been experimented widely on small animals and in patient trials.^13,17–19

One of the crucial challenges in linear-array-based PAT is the image reconstruction.^20–22 The presence of noise and artifacts in the detected signals and simplifications of the reconstruction algorithms for speed degrade the PA image quality.²³ Indeed, there are some inherent artifacts, caused by the image formation algorithms, observed in the reconstructed images.²⁴ Moreover, due to the limited view of the linear-array transducer (in comparison with the ring-array transducer, which obtains a full 360 deg view) reconstructed image quality is poor.^10,11 In other words, linear-array transducers have a view of approximately only 40 deg, and the image quality is lower compared with circular tomography.^25–27

Beamforming algorithms, which are commonly used in Radar and US imaging, can be used in PAT with some modifications.²⁸ There are some studies conducted to use a single beamforming algorithm for the integrated US/PA imaging systems.^25,29,30 Delay-and-sum (DAS) can be considered as the most common beamforming algorithm in US and PAT.^31–33 To address the in-capabilities in DAS, which is mainly due to its blindness, minimum variance (MV) can be used.^34,35 In MV, all the calculated samples for each point of imaging are weighted adaptively, resulting in a significant resolution, but it should be noticed that the sidelobes will not be well degraded. On the other hand, delay-multiply-and-sum (DMAS) can be used to suppress the noise and sidelobes, and improve the image quality.³⁶ To address the low noise suppression of DMAS at the presence of a high level of noise, double-stage DMAS (DS-DMAS), in which two stages of DMAS is used, is recently introduced for linear-array US and PAT.^37,38 Despite all the improvement gained by DMAS and DS-DMAS, the resolution improvement is not satisfying, compared with MV. The matter of low contrast of MV and low-resolution of DMAS is addressed using the MV-based DMAS (MVB-DMAS), where the combination of these two methods is used for PA image reconstruction.^28,39 Eigenspace-based MV (EIBMV) is introduced for US image formation to degrade the sidelobes and improve the contrast obtained by the MV.^40,41 The concept of MVB-DMAS is applied to the EIBMV and used for linear-array PAT.⁴² Coherence factor (CF) is applied to the MV beamformed signals to improve the resolution and suppress the sidelobes.⁴³ Two modifications of CF are introduced for linear-array PAT to have a lower sidelobes and higher resolution, compared with the conventional CF.^44,45

In this work, an image formation algorithm for linear-array PAT is proposed. The proposed algorithm is a nonlinear (NL) beamformer, and it uses the $p$’th root of the detected PA signals. The main improvement gained by the proposed method is lower sidelobes and noise while the complexity of the algorithm is on the order of DAS [$O(M)$]. The results show (both numerical and experimental) that the proposed algorithm can be an appropriate choice for linear-array PA image formation, especially when a high level of noise affects the image quality.

2.

Methods and Materials

2.1.

Beamforming

Assuming a linear geometry for PA waves detection, the optical absorption distribution of the tissue can be reconstructed using DAS, which can be written as follows:

Eq. (1)

$$y_{\mathrm{DAS}}(k)=\sum _{i=1}^M{x}_i(k-{\mathrm{\Delta }}_i),$$

where $y_{\mathrm{DAS}}(k)$, $k$, $M$ are the output of the beamformer, the time index, and the number of elements of the array, respectively. In addition, $x_i(k)$ and ${\mathrm{\Delta }}_i$ are the detected signals and the corresponding time delay for the detector $i$, respectively. As mentioned before, DAS results in a low-quality image and is commonly used due to its simple implementation. To improve the contrast gained by DAS, DMAS is introduced,³⁶ which is as follows:

Eq. (2)

$$y_{\mathrm{DMAS}}(k)=\sum _{i=1}^{M-1}\sum _{j=i+1}^M{x}_i(k-{\mathrm{\Delta }}_i)x_j(k-{\mathrm{\Delta }}_j).$$

The dimensionally squared problem of Eq. (2) is addressed as follows:³⁶

Eq. (3)

$$x_{ij}^{\prime }(k)=\mathrm{sign}[x_i(k-{\mathrm{\Delta }}_i)x_j(k-{\mathrm{\Delta }}_j)]\sqrt{|x_i(k-{\mathrm{\Delta }}_i)x_j(k-{\mathrm{\Delta }}_j)|},$$

Eq. (4)

$$y_{\mathrm{DMAS}}(k)=\sum _{i=1}^{M-1}\sum _{j=i+1}^M{x}_{ij}^{\prime }(k).$$

The new components appear in the spectrum due to the similar ranges of frequency for $x_i(k-{\mathrm{\Delta }}_i)$ and $x_j(k-{\mathrm{\Delta }}_j)$, and the multiplication procedure in the DMAS algorithm. A bandpass filter is applied on the beamformed output signal to only pass the necessary frequency components while keeping the one centered on $2f_0$ almost unaltered. As a result of the used filter, it is named filtered-DMAS (F-DMAS), extensively evaluated in the reference.³⁶ In the previous publications, the superiority of the DMAS has been proved in the terms of sidelobes, resolution, and contrast. However, all the advantages are achieved at the expense of a higher computational burden. DMAS is a nonlinear beamformer as a result of its correlation procedure. In this paper, it is proposed to use the $p$’th root of the detected signals to improve the contrast of the PA images.⁴⁶ The $p$’th root of the signals is used as the input of the DAS algorithms. The proposed method formula can be written as follows:

Eq. (5)

$$y_{NL\_p}(k)={\{\frac{1}{M}\sum _{i=1}^M\mathrm{sign}[x_i(k-{\mathrm{\Delta }}_i)]\sqrt[p]{|x_i(k-{\mathrm{\Delta }}_i)|}\}}^p.$$

As can be seen in Eq. (5), the $p$’th root of the detected signals is used in the summation procedure. The same as DMAS algorithms, which changes the dimension to the order of ${\mathrm{Volt}}^2$, the dimension of the signals would be changed to ${\mathrm{Volt}}^{1/p}$ in the NL algorithm. This problem is addressed by the ($p$’th) power finally used (bringing back the dimension to the Volt). It should be noticed that for $p=2$, the NL algorithm would be a close approximation of the DMAS beamformer. To illustrate this, consider the expansion of Eq. (5) when $p=2$

Eq. (6)

$$y_{NL\_2}(k)={\{\frac{1}{M}\sum _{i=1}^M\mathrm{sign}[x_i(k-{\mathrm{\Delta }}_i)]\sqrt[2]{|x_i(k-{\mathrm{\Delta }}_i)|}\}}^2={[\frac{1}{M}\sum _{i=1}^M{x}_i^{\prime }(k)]}^2=[\frac{1}{M}\sum _{i=1}^M{x}_i^{\prime }(k)][\frac{1}{M}\sum _{i=1}^M{x}_i^{\prime }(k)]=\frac{1}{M^2}[x_1^{\prime }(k)+x_2^{\prime }(k)+\dots +x_M^{\prime }(k)][x_1^{\prime }(k)+x_2^{\prime }(k)+\dots +x_M^{\prime }(k)].$$

By multiplying the term in the parentheses in Eq. (6), the following equation is generated:

Eq. (7)

$$\frac{1}{M^2}\{[x_1^{\prime 2}(k)+x_1^{\prime }(k)x_2^{\prime }(k)+x_1^{\prime }(k)x_3^{\prime }(k)+\dots +x_1^{\prime }(k)x_M^{\prime }(k)]+[x_2^{\prime }(k)x_1^{\prime }(k)+x_2^{\prime 2}(k)+x_1^{\prime }(k)x_3^{\prime }(k)+\dots +x_1^{\prime }(k)x_M^{\prime }(k)]+[x_M^{\prime }(k)x_1^{\prime }(k)+x_M^{\prime }(k)x_2^{\prime }(k)+\dots +x_M^{\prime }(k)x_{M-1}^{\prime }(k)+x_M^{\prime 2}(k)]\}.$$

As shown in Eq. (7), there are some terms with the power of 2 [i.e., $x_1^{\prime 2}(k)$, $x_2^{\prime 2}(k)$, $x_M^{\prime 2}(k)$]. All these terms can be written as follows:

Eq. (8)

$$\frac{1}{M^2}[x_1^{\prime 2}(k)+x_2^{\prime 2}(k)+\dots +x_M^{\prime 2}(k)]=\frac{1}{M^2}\genfrac{}{}{0ex}{}{\underbrace{\sum _{i=1}^M{x}_i^{\prime 2}(k)}}{\mathrm{a}\mathrm{DAS}\text{on the}p\prime \mathrm{th}\text{power of the signals,}}$$

where, as is clarified in the equation, there is a DAS implemented on the $p$’th power of the signals. The rest of the terms of Eq. (7), those that are not in Eq. (8), can be mathematically written as follows:

Eq. (9)

$$\frac{2}{M^2}\genfrac{}{}{0ex}{}{\underbrace{\sum _{i=1}^{M-1}\sum _{j=i+1}^M{x}_i^{\prime }(k)x_j^{\prime }(k)}}{\mathrm{DMAS}}.$$

Finally, the Eq. (6) is mathematically written as follows:

Eq. (10)

$$\frac{1}{M^2}\genfrac{}{}{0ex}{}{\underbrace{\sum _{i=1}^M{x}_i^{\prime 2}(k)}}{\mathrm{a}\mathrm{DAS}\text{on the}p’\mathrm{th}\text{power of the signals}}+\frac{2}{M^2}\genfrac{}{}{0ex}{}{\underbrace{\sum _{i=1}^{M-1}\sum _{j=i+1}^M{x}_i^{\prime }(k)x_j^{\prime }(k)}}{\mathrm{DMAS}}.$$

As can be seen, Eq. (6) has led to two terms (a DAS and a DMAS). The results in the next section show that Eq. (6) would be a close approximation of DMAS. In fact, the performance of DMAS has been achieved while only the computational complexity of DAS has been expended. A flowchart is shown in Fig. 1 to better illustrate the proposed reconstruction method. It should be noticed that for even $p$-values, the sign of the PA signals will be lost by the final power of $p$, leading to the splitting of the original spectrum into DC components and doubled spectrum band of the original signal. Therefore, at the end of the procedure of image formation, a bandpass filter should be applied on the beamformed PA signals to only pass the necessary frequency components. In the next section, the numerical and experimental results of the proposed method for linear-array PAT are reported.

Fig. 1

The flowchart of the proposed NL beamformer.

2.4.

In-Vivo Imaging of Sentinel Lymph Node

The experimental set-up for the in-vivo imaging of asentinel lymph node is shown in Fig. 2. A laser beam of 1064 nm (the same as the wire-target phantom experiment) is focused into a fiber bundle using a 150-mm converging lens. A beam-splitter is used to reflect a small percentage of the beam to photodiode, which is needed to trigger the clinical ultrasound system. The fiber bundle has 1600 multimode fibers that are separated into two bundles at the output. The two bundles are held across the ultrasound transducer using a custom-designed three-dimensional (3-D) printed holder. The angle of illumination from both the fiber bundles is 15 deg. Adult Sprague Dawley rat of 250 gms is initially anesthetized using a cocktail of ketamine ($85\mathrm{mg}/\mathrm{kg}$) and xylazine ($15\mathrm{mg}/\mathrm{kg}$). The hair in the region of interest is depleted and then inhalation anesthesia of $1\mathrm{L}/\min$ oxygen and 0.75% isoflurane (Euthanex Corp.) is given. The fur on the scalp of the animal is removed using a hair clipper. The hair removal cream (Veet, Reckitt Benckiser) is gently applied to the shaved area for further depletion of the fur. The applied cream is removed after 4 to 5 min using a cotton swab. A chicken breast tissue of 5 mm thickness is placed on the animal to mimic SLN imaging of a human. Black ink is injected in the forepaw of the animal and massaged for the ink to flow into the lymph node [Fig. 2(b)]. The region of interest is illuminated with $20\mathrm{mJ}/{\mathrm{cm}}^2$ of energy, which is within the maximum permissible limit of $100\mathrm{mJ}/{\mathrm{cm}}^2$ for 1064 nm as per American National Standard for Safe Use of Lasers.⁵¹ After the acquisition of PA images, the animal is euthanized with an overdose of pentobarbital. An incision is made to expose the SLN [Fig. 2(c)]. The SLN is excised [Fig. 2(d)] from the animal. All experiments are performed in accordance with the guidelines and regulations approved by the Institutional Animal Care and Use Committee of Nanyang Technological University, Singapore (Animal Protocol Number ARF-SBS/NIE-A0263).

Fig. 2

Experimental setup of in-vivo imaging, (a) experimental setup, (b) photograph of rat used for in-vivo imaging, (c) photograph of the rat exposing the sentinel lymph node after incision, and (d) excised sentinel lymph node stained black.

3.

Results

3.1.

Numerical Study

First, the performance of the proposed method is evaluated for $p=2$ and $p=3$, compared with DAS and DMAS. The reconstructed PA images are shown in Fig. 3, where noise was added to the detected signals considering a signal-to-noise ratio (SNR) of 30 dB. As can be seen in Fig. 3(a), DAS results in high sidelobes, which degrades the image quality, and the effects of the added noise are obvious in the image. DMAS reduces the sidelobes and artifacts, but they still affect the image. NL_2 provides an image the same as DMAS, whereas NL_3 provides a higher noise and sidelobes suppression, which leads to a higher image quality in comparison with other mentioned beamformers. To evaluate the performance of the beamformers in details, the lateral variations of the images shown in Fig. 3 are presented in Fig. 4. The same performance of NL_2 and DMAS is clear considering their lateral variations, which are almost the same, while they outperform DAS in the terms of sidelobes and lateral valley. Considering Fig. 4(a), it can be seen that the level of sidelobes for DAS, DMAS, and NL_2 is about $-25$, $-37$, and $-37\mathrm{dB}$, respectively. On the other hand, NL_3 provides a lower sidelobes and lateral valley. As a result, NL_3 degrades the sidelobes for about 21, 9, and 9 dB, compared with DAS, DMAS, and NL_2, respectively. We have calculated full-width-half-maximum (FWHM) to evaluate the lateral resolution of the beamformers. At the depth of 40 mm, DAS, DMAS, NL_2, and NL_3 lead to a FWHM of 2.05, 1.43, 1.43, and 1.17 mm, respectively. Thus, NL_3 results in lower FWHM, indicating a better resolution, whereas it is not a significant improvement compared with those obtained in previous studies.^28,39 These results indicate that the NL method can be used to provide a high image quality.

Fig. 3

Images of the simulated point-target phantom. (a) DAS, (b) DMAS, (c) NL_2, and (d) NL_3. All the images are shown with a dynamic range of 60 dB. Noise was added to the detected signals considering a SNR of 30 dB.

Fig. 4

Lateral variations of DAS, DMAS, NL_2, and NL_3 at the depths of (a) 35 mm and (b) 45 mm using the images shown in Fig. 3.

3.1.1.

High level of medium noise

To further evaluate the NL beamformer, noise is added to the detected signals considering a SNR of 0 dB, and the reconstructed image is shown in Fig. 5. The presence of high-power noise clearly reduces the quality of the image obtained by DAS, whereas the effects are mitigated using DMAS and NL_2. Considering the image shown in Figs. 5(b) and 5(c), it is demonstrated that the noise still degrades the image quality. However, the effects of the noise are highly suppressed using NL_3, as shown in Fig. 5(d). The lower sidelobes, noise, and lateral valley, gained by the NL_3, in the presence of the powerful noise, are shown in Fig. 6 using the lateral variations.

Fig. 5

Images of the simulated point-target phantom. (a) DAS, (b) DMAS, (c) NL_2, and (d) NL_3. All the images are shown with a dynamic range of 60 dB. Noise was added to the detected signals considering a SNR of 0 dB.

Fig. 6

Lateral variations of DAS, DMAS, NL_2, and NL_3 at the depths of (a) 35 mm and (b) 45 mm using the images shown in Fig. 5.

To quantitatively compare the methods, the signal-to-noise ratio (SNR) has been used, which is as follows:

Eq. (11)

$$\mathrm{SNR}=20{\log }_{10}{P}_{\text{signal}}/P_{\text{noise}},$$

where $P_{\text{signal}}$ and $P_{\mathrm{noise}}$ are difference between maximum and minimum intensities of a rectangular region including a point target [white dashed rectangle in Fig. 5(d)] and standard deviation of the noisy part of the region [green rectangle in Fig. 5(d)], respectively.^37,38 The calculated SNRs are shown in Table 1, at the all depths of imaging. As demonstrated, the proposed method when $p=2$ performs almost the same as DMAS while both of them provide a higher SNR compared with DAS. This is mainly due to higher noise and sidelobes suppression. It should be noticed that the same performance is achieved while the computational burden of the NL_2 algorithm is on the order of magnitude of $M$ [$O(M)$], which is lower than DMAS with [$O(M^2)$]. On the other hand, increasing the root used in the proposed method to 3 (NL_3) would result in a higher noise suppression, which leads to a higher SNR. Consider, for instance, the depth of 45 mm, where DAS, DMAS, NL_2, and NL_3 result in a SNR of about 23.71, 32.07, 32.17, and 38.98 dB, respectively. In other words, NL_3 improves the SNR for about 15.27, 6.91, and 6.91 dB in comparison with DAS, DMAS, and NL_2, respectively.

Table 1

SNR (dB) values at the different depths.

Depth (mm)	DAS	DMAS	NL_2	NL_3
25	32.85	46.50	46.52	59.33
30	30.50	44.74	44.80	57.01
35	29.70	43.30	43.36	54.96
40	25.17	35.5	35.69	44.65
45	23.71	32.07	32.17	38.98
50	21.39	28.55	28.64	34.81

3.1.2.

Effects of the parameter $P$

Here, we evaluate the effects of $p$ on the reconstructed PA images. The reconstructed images using different $P\mathrm{s}$ are shown in Fig. 7, where increasing the amount of $p$ would improve the image quality and degrade the sidelobes. The further evaluation can be conducted using the lateral variations as shown in Fig. 8, where the bigger the amount of $p$, the lower the sidelobes. It can be seen that in each step of increasing $p$, the sidelobes are degraded for about 13 dB. The FWHM improvement is not significant when $p$ is increased while the SNR improvement would be considerable as the sidelobes are degraded. Considering the presented numerical results, one might come to the conclusion that increasing the parameter $p$ would always result in image improvement. However, in the Sec. 4, the effect of parameter $p$ is investigated using experimental data, and it is shown that a large or wrong $p$ would result in image disturbance and signal removal. In other words, when it comes to selecting the best $p$ (providing the highest SNR), we need to consider the main signal removal of the NL beamformer as well as the noise/artifacts removal.

Fig. 7

Images of the simulated point-target phantom. (a) NL_2, (b) NL_3, (c) NL_4, (d) NL_5, and (e) NL_6. All the images are shown with a dynamic range of 60 dB. Noise was added to the detected signals considering an SNR of 30 dB.

Fig. 8

Lateral variations of NL_2, NL_3, NL_4, NL_5, NL_6, NL_7, NL_8, and NL_9 at the depths of 32.5 mm using the images shown in Fig. 7.

3.2.

Wire-Target Phantom Experiment

The reconstructed experimental images are shown in Fig. 9 (zoomed version in Fig. 10), where DAS leads to a low-quality image having a high sidelobes, artifacts, and noise. DMAS and NL_2 result in the same image quality [Figs. 9(b) and 9(c)], but increasing the root of the NL beamformer leads to image quality improvement, as shown in [Figs. 9(e) and 9(f)]. To have a better look at the gained improvement, consider the lateral variations shown in Fig. 11, where the NL_5 method outperforms DAS, DMAS, and NL with lower $p$ (see the circle and arrows). It should be noticed that, as mentioned in the Sec. 2, the performance of the DMAS and NL_2 is almost the same, which is also proved with the experimental data, considering their lateral variations. To quantitatively evaluate the NL beamformer, we have calculated the SNR for the experimental PA images. The SNRs are shown in Table 2 where, at the depth of 11.3 mm NL_5 improves the SNR for about 19.1, 14.5, 14.5, 9.8, and 4.5 dB, compared with DAS, DMAS, NL_2, NL_3, and NL_4, respectively. All the experimental results indicate the superiority of the NL beamformer compared with DAS and DMAS.

Fig. 9

The reconstructed experimental images using (a) DAS, (b) DMAS, (c) NL_2, (d) NL_3, (e) NL_4, and (f) NL_5. A wire-target phantom was utilized. All the images are shown with a dynamic range of 70 dB.

Fig. 10

A zoomed version of the Fig. 9 for better comparison. (a) DAS, (b) DMAS, (c) NL_2, (d) NL_3, (e) NL_4, and (f) NL_5.

Fig. 11

The lateral variations for the PA images shown in Fig. 9.

Table 2

SNR (dB) values for the experimental image shown in Fig. 9.

Depth (mm)	DAS	DMAS	NL_2	NL_3	NL_4	NL_5
7.1	35.8	44.0	44.0	50.7	56.4	59.9
11.3	24.8	29.4	29.4	34.1	39.4	43.9

3.3.

In-Vivo Imaging of Sentinel Lymph Node

The reconstructed in-vivo images are shown in Fig. 12, where the position of the SLN has been indicated with the arrow [see Fig. 12(c)]. As can be seen in Fig. 12(f), the background noise of the image reconstructed by the NL_5 is better suppressed in comparison with other images. In other words, the NL beamformer (for $p>2$) provides higher image quality compared with DAS and DMAS. It should be mentioned that at the wavelength of 1064 nm, the blood vessels have a very poor contrast. Therefore, they do not show up or get buried in the background.

Fig. 12

The reconstructed in-vivo images using (a) DAS, (b) DMAS, (c) NL_2, (d) NL_3, (e) NL_4, and (f) NL_5. All the images are shown with a dynamic range of 60 dB.

4.

Discussion

In this work, we have introduced an image reconstruction algorithm (NL beamformer) for PAT in the case that a linear-array US transducer has been used for data acquisition. We have extensively evaluated the proposed method for photoacoustic imaging in the terms of SNR, FWHM, level of sidelobes, and target separability. In addition, we have evaluated the effects of the parameter p on the reconstructed images, and employed the proposed method for SLN imaging, for the first time. The DAS algorithm considers all the signals the same, which is why it results in a low-quality image having a high sidelobes and artifacts. The DMAS beamformer uses a correlation procedure to suppress the noise and the sidelobes caused by the simple summation in DAS beamformer, and finally results in higher image quality compared with DAS. The main improvements gained by the proposed method are the lower sidelobes and a higher SNR, compared with DMAS, while a lower computational burden is imposed. As can be seen in Figs. 3 and 5, the NL beamformer results in a high-image quality, especially when $p$ is increased. This is because of the suppression of the weak signals while the powerful signals are maintained. To put it more simply, using the $p$’th root of the detected PA signals suppresses the artifacts and noise. Although the powerful samples are maintained, the contribution of the weak samples of the detected PA signals would be lower. It is expected to detect the targets using the powerful samples, and the $p$’th root makes it possible. The lateral variations shown in Figs. 4 and 6, where the proposed method has been evaluated at the presence of high levels of imaging noise and medium heterogeneities, prove the superiority of the NL algorithm. Table 1 shows the higher noise suppression of the NL algorithm that, as was mentioned, results from the usage of the $p$’th root. Despite all the promising numerical results, it was necessary to evaluate the NL beamformer with experimental data. Figures 9 and 11 show the reconstructed PA image when a wire-target phantom is used. The experimental results indicate that the NL beamformer (with $p>2$) outperforms DAS and DMAS, which makes it an appropriate algorithm for PA image reconstruction. The superiority of the NL beamformer is also evaluated with an in-vivo experiment, and the results are shown in Fig. 12. To pass the necessary information and attenuate the DC component generated after the NL_p (with an even p), a bandpass filter was applied by a Tukey window ($\alpha =0.5$) to the beamformed PA signal spectra, covering 4.5 to 11.5 MHz and 11 to 19 MHz for numerical and experimental studies, respectively. The algorithm is independent of the laser wavelength that is used. Thus, in experiments, there is no reason to use only a single wavelength for evaluation. At 1064 nm, the black ink has a strong optical absorption, so it is used as the contrast agent. Any other contrast agents having strong absorption can be used.

The matter of selection of $p$ is of importance in the NL beamformer. As mentioned, using the p’th root degrades the artifacts and noise contribution in the procedure of PA image reconstruction, and the numerical results indicated that the larger the parameter $p$, the higher the image quality. However, if the parameter $p$ is selected so large (larger than the best $p$), it could suppress the signals caused by the target of imaging as well as the noise and artifacts. To illustrate this, consider Fig. 13 where large $p$ is used for reconstruction as well as a low amount of $p$. It can be seen that having a $p$ of 40, or even 30, would suppress the main target of imaging as well as the background noise and sidelobes. In other words, wrong selection of $p$ may result in a disturbed and nonsense image quality, as can be seen in Figs. 13(e) and 13(f). Here, we show the application of the proposed algorithm for SLN imaging, but the algorithm can be used in other applications as well. Due to our investigation, for SLN imaging, selection of $p=12$ would result in the best image quality (can be seen in Fig. 14) without compromising the main target of imaging (SLN). It should be noted that $p=12$ was obtained empirically. An adaptive way of choosing the best $p$, based on the energy or power of the detected PA signals, would be a matter of investigation for our future studies.

Fig. 13

The effects of $p$ on the experimental images shown in Fig. 10. (a) $p=5$, (b) $p=9$, (c) $p=14$, (d) $p=20$, (e) $p=30$, and (f) $p=40$.

Fig. 14

The reconstructed PA image using NL_12, where the yellow arrow shows the target.

Beamformers usually can be utilized in the US and PA image formations. The proposed method changes the dynamic range of the detected signals, and we can use it to improve the image quality. The proposed nonlinear reconstruction will certainly affect the quantitative values of the signal, making it ineligible for a quantitative data analysis. However, PAI is not only about multispectral imaging. There are other applications such as structural imaging as well. Thus, in the case of structural imaging, the proposed algorithm would be useful. As shown in Figs. 11 and 8, even though the pattern of noise is not changed, its level is reduced. We can use this reduction to improve the image quality in PAI. It is worth to mention that the proposed method in this paper, especially when a $p>3$ is selected, may not be an appropriate algorithm for US imaging. This is mainly because of the importance of speckles in US images. To put it more simply, in most of US images, the speckles contain important information for diagnosis and detection. To this end, using $p>3$ for US imaging is not suggested as it would result in speckle removal, which is not desired. In Ref. 36, it has been mentioned that the computational burden of DMAS is on the order of $M^2$ [$O(M^2)$], whereas the order of DAS is $M$ [$O(M)$]. The formula of the NL beamformer is shown in Eq. (5), and in Sec. 2, it is proved that NL_2 could be a close approximation of the DMAS algorithm while it only imposes an order of $M$. This is mainly due to the fact that there is just a summation procedure in the proposed method, and all the improvements by NL are achieved while a lower computational burden is imposed, compared with DMAS. Of note, considering the sign and calculation of p’th root in the computational complexity, the NL_p algorithm imposes a complexity on the order of $M*\log (M)$. The size of the experimental images is $256\times 256\text{pixel}$, where DAS, DMAS, and NL_p take a time of 0.16, 17.9, and 0.26 s to reconstruct the images, respectively. In addition, the size of numerical images is $550\times 200\text{pixel}$, where DAS, DMAS, and NL_p take a time of 0.18, 31.69, and 0.40 s to reconstruct the images, respectively. It deserves to be mentioned that the proposed method can be implemented on a processing chip such as FPGA and be used in the investigation and commercial PA devices. As mentioned above, the selection of parameter $p$ highly affects the performance of the proposed method, especially when it comes to choosing the best $p$. Thus, implementation of the proposed algorithm on an FPGA and designing an adaptive way to select the parameter $p$ can be considered as further works.

5.

Conclusion

In this work, an image formation algorithm (NL_p) was utilized for linear-array PAT. Although DMAS improves the image quality, compared with DAS, it imposes a computational burden of $O(M^2)$. The introduced method uses the $p$’th root of the detected signals, and a summation, resulting in $O(M)$. The NL algorithm was evaluated numerically and experimentally (wire-target and in-vivo SLN imaging). The results indicated that NL outperforms DMAS while it imposes the computational burden of DAS. The reconstructed experimental images of the wire-target phantom showed that, at the depth of 11.3 mm, NL_5 leads to SNR improvement of about 77%, 49%, 49%, 28%, and 11% in comparison with DAS, DMAS, NL_2, NL_3, and NL_4, respectively. In addition, the appropriate $p$ for SLN imaging was obtained to be 12.