Translator Disclaimer
27 July 2018 Efficient nonlinear beamformer based on P’th root of detected signals for linear-array photoacoustic tomography: application to sentinel lymph node imaging
Author Affiliations +
In linear-array transducer-based photoacoustic (PA) imaging, B-scan PA images are formed using the raw channel PA signals. Delay-and-sum (DAS) is the most prevalent algorithm due to its simple implementation, but it leads to low-quality images. Delay-multiply-and-sum (DMAS) provides a higher image quality in comparison with DAS while it imposes a computational burden of O  (  M2  )  . We introduce a nonlinear (NL) beamformer for linear-array PA imaging, which uses the p’th root of the detected signals and imposes the complexity of DAS [O  (  M  )  ]. The proposed algorithm is evaluated numerically and experimentally [wire-target and in-vivo sentinel lymph node (SLN) imaging], and the effects of the parameter p are investigated. The results show that the NL algorithm, using a root of p (NL_p), leads to lower sidelobes and higher signal-to-noise ratio compared with DAS and DMAS, for (p  >  2). The sidelobes level (for the wire-target phantom), at the depth of 11.4 mm, are about −31, −52, −52, −67, −88, and −109  dB, for DAS, DMAS, NL_2, NL_3, NL_4, and NL_5, respectively, indicating the superiority of the NL_p algorithm. In addition, the best value of p for SLN imaging is reported to be 12.



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.5 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.6 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.10 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.13 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).16 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,1719

One of the crucial challenges in linear-array-based PAT is the image reconstruction.2022 The presence of noise and artifacts in the detected signals and simplifications of the reconstruction algorithms for speed degrade the PA image quality.23 Indeed, there are some inherent artifacts, caused by the image formation algorithms, observed in the reconstructed images.24 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.2527

Beamforming algorithms, which are commonly used in Radar and US imaging, can be used in PAT with some modifications.28 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.3133 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.36 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.42 Coherence factor (CF) is applied to the MV beamformed signals to improve the resolution and suppress the sidelobes.43 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.


Methods and Materials



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)

where yDAS(k), k, M are the output of the beamformer, the time index, and the number of elements of the array, respectively. In addition, xi(k) and Δ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,36 which is as follows:

Eq. (2)


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

Eq. (3)


Eq. (4)


The new components appear in the spectrum due to the similar ranges of frequency for xi(kΔi) and xj(kΔ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 2f0 almost unaltered. As a result of the used filter, it is named filtered-DMAS (F-DMAS), extensively evaluated in the reference.36 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.46 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)


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 Volt2, the dimension of the signals would be changed to Volt1/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)


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

Eq. (7)


As shown in Eq. (7), there are some terms with the power of 2 [i.e., x12(k), x22(k), xM2(k)]. All these terms can be written as follows:

Eq. (8)

1M2[x12(k)+x22(k)++xM2(k)]=1M2i=1Mxi2(k)aDASon thepthpower 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)


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

Eq. (10)

1M2i=1Mxi2(k)aDASon thepthpower of the signals+2M2i=1M1j=i+1Mxi(k)xj(k)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.



Numerical Study

The K-wave MATLAB toolbox is used to perform the numerical study.47 Six pairs of point targets, having a radius of 0.1 mm, are positioned at the depth of 25 mm until 50 mm, separated 5 mm in the axial axis and 4 mm in the lateral axis, and two single-point targets are positioned at the depth of 32.5 and 42.5 mm, respectively. A linear-array having M=128 elements operating at 4-MHz central frequency and 77% fractional bandwidth is used to detect the PA signals. The speed of sound is assumed to be 1540  m/s during simulations. The sampling frequency is 50 MHz.


Wire-Target Phantom Experiment

Experimental PA images were acquired using the clinical ultrasound system (E-CUBE 12R, Alpinion, South Korea). The ultrasound transducer has 128 elements over a length of 3.85  cm×1  cm, center frequency of 8.5 MHz, and the fractional bandwidth is 95%. A laser beam of 532 nm from Nd:YAG pump laser (Continuum, Surelite Ex, San Jose, California) is used for excitation. The laser has a pulse repetition rate of 10 Hz.13,48,49 A lens of focal length 50  mm is used to diverge the laser beam to illuminate the target, which is four pencil leads of diameter 0.5 mm positioned over an area of 25  mm×25  mm. About 1% of the laser beam is reflected to a photodiode, which is used as a trigger signal to synchronize the ultrasound system and the laser excitation.50 For every trigger, the ultrasound system acquires channel data from 64 elements of the array transducer. Hence, the PA images are acquired at a frame rate of 5  frames/s. Acquired radio frequency data were saved in the local machine and used later for testing the reconstruction algorithm. The pencil leads (Mountain peak, China) were pinned in clay and immersed in water. The laser passed through the wall of the water tank before reaching the target. Leads were at a depth of 5 to 15 mm from the transducer.


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  mg/kg) and xylazine (15  mg/kg). The hair in the region of interest is depleted and then inhalation anesthesia of 1  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  mJ/cm2 of energy, which is within the maximum permissible limit of 100  mJ/cm2 for 1064 nm as per American National Standard for Safe Use of Lasers.51 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.





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  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.



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)

where Psignal and Pnoise 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(M2)]. 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)DASDMASNL_2NL_3


Effects of the parameter P

Here, we evaluate the effects of p on the reconstructed PA images. The reconstructed images using different Ps 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.



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)DASDMASNL_2NL_3NL_4NL_5


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.




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 (α=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 M2 [O(M2)], 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×256  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×200  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.



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(M2). 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.


Authors have no relevant financial interests in the manuscript and no other potential conflicts of interest to disclose.


The authors VP and MP acknowledge Dr. Dienzo Rhonnie Austria for his help in animal experiments. The research is supported by Tier 2 grant funded by the Ministry of Education in Singapore (ARC2/15: M4020238 to M.P.).



L. V. Wang and J. Yao, “A practical guide to photoacoustic tomography in the life sciences,” Nat. Methods, 13 (8), 627 –638 (2016). JBOPFO 1083-36681548-7091 Google Scholar


M. Pramanik et al., “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (ta) and photoacoustic (pa) tomography,” Med. Phys., 35 (6), 2218 –2223 (2008). MPHYA6 0094-2405 Google Scholar


L. Li et al., “Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution,” Nat. Biomed. Eng., 1 (5), 0071 (2017). Google Scholar


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


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). Google Scholar


Q. Sheng et al., “A constrained variable projection reconstruction method for photoacoustic computed tomography without accurate knowledge of transducer responses,” IEEE Trans. Med. Imaging, 34 (12), 2443 –2458 (2015). ITMID4 0278-0062 Google Scholar


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


A. Hariri et al., “Development of low-cost photoacoustic imaging systems using very low-energy pulsed laser diodes,” J. Biomed. Opt., 22 (7), 075001 (2017). JBOPFO 1083-3668 Google Scholar


P. K. Upputuri and M. Pramanik, “Dynamic in vivo imaging of small animal brain using pulsed laser diode-based photoacoustic tomography system,” J. Biomed. Opt., 22 (9), 090501 (2017). JBOPFO 1083-3668 Google Scholar


P. K. Upputuri and M. Pramanik, “Recent advances toward preclinical and clinical translation of photoacoustic tomography: a review,” J. Biomed. Opt., 22 (4), 041006 (2016). Google Scholar


K. Sivasubramanian et al., “Hand-held, clinical dual mode ultrasound-photoacoustic imaging of rat urinary bladder and its applications,” J. Biophotonics, 11 (5), e201700317 (2018). JBOPFO 1083-3668 Google Scholar


J. Kim et al., “Programmable real-time clinical photoacoustic and ultrasound imaging system,” Sci. Rep., 6 35137 (2016). SRCEC3 2045-2322 Google Scholar


K. Sivasubramanian, V. Periyasamy and M. Pramanik, “Non-invasive sentinel lymph node mapping and needle guidance using clinical handheld photoacoustic imaging system in small animal,” J. Biophoton., 11 (1), e201700061 (2018). Google Scholar


K. K. Swenson et al., “Comparison of side effects between sentinel lymph node and axillary lymph node dissection for breast cancer,” Ann. Surg. Oncol., 9 (8), 745 –753 (2002). Google Scholar


K. V. Ballman, L. M. McCall and A. E. Giuliano, “Axillary vs sentinel lymph node dissection in women with invasive breast cancer-reply,” J. Am. Med. Assoc., 319 (3), 306 –307 (2018). Google Scholar


H. Deka et al., “A study of sentinel lymph node biopsy with methylene blue dye in early carcinoma of breast,” J. Evol. Med. Dent. Sci., 6 (21), 1701 –1705 (2017). Google Scholar


M. Pramanik et al., “In vivo carbon nanotube-enhanced non-invasive photoacoustic mapping of the sentinel lymph node,” Phys. Med. Biol., 54 (11), 3291 –3301 (2009). PHMBA7 0031-9155 Google Scholar


T. N. Erpelding et al., “Sentinel lymph nodes in the rat: noninvasive photoacoustic and us imaging with a clinical us system,” Radiology, 256 (1), 102 –110 (2010). RADLAX 0033-8419 Google Scholar


A. Garcia-Uribe et al., “Dual-modality photoacoustic and ultrasound imaging system for noninvasive sentinel lymph node detection in patients with breast cancer,” Sci. Rep., 5 15748 (2015). SRCEC3 2045-2322 Google Scholar


K. S. Uddin et al., “Two step imaging reconstruction using truncated pseudoinverse as a preliminary estimate in ultrasound guided diffuse optical tomography,” Biomed. Opt. Express, 8 (12), 5437 –5449 (2017). BOEICL 2156-7085 Google Scholar


Y. Lou et al., “Generation of anatomically realistic numerical phantoms for photoacoustic and ultrasonic breast imaging,” J. Biomed. Opt., 22 (4), 041015 (2017). JBOPFO 1083-3668 Google Scholar


M. A. Anastasio, P. C. Beard and Q. Zhu, “Special section guest editorial: photoacoustic imaging and sensing,” J. Biomed. Opt., 22 (4), 041001 (2017). JBOPFO 1083-3668 Google Scholar


Y. Han, V. Ntziachristos and A. Rosenthal, “Optoacoustic image reconstruction and system analysis for finite-aperture detectors under the wavelet-packet framework,” J. Biomed. Opt., 21 (1), 016002 (2016). JBOPFO 1083-3668 Google Scholar


A. Rosenthal et al., “Efficient framework for model-based tomographic image reconstruction using wavelet packets,” IEEE Trans. Med. Imaging, 31 (7), 1346 –1357 (2012). ITMID4 0278-0062 Google Scholar


E. Mercep et al., “Hybrid optoacoustic tomography and pulse-echo ultrasonography using concave arrays,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 62 (9), 1651 –1661 (2015). ITUCER 0885-3010 Google Scholar


K. Daoudi et al., “Handheld probe integrating laser diode and ultrasound transducer array for ultrasound/photoacoustic dual modality imaging,” Opt. Express, 22 (21), 26365 –26374 (2014). OPEXFF 1094-4087 Google Scholar


H. Ke et al., “Performance characterization of an integrated ultrasound, photoacoustic, and thermoacoustic imaging system,” J. Biomed. Opt., 17 (5), 056010 (2012). Google Scholar


M. Mozaffarzadeh, A. Mahloojifar and M. Orooji, “Medical photoacoustic beamforming using minimum variance-based delay multiply and sum,” Proc. SPIE, 10335 1033522 (2017). PSISDG 0277-786X Google Scholar


H. K. Zhang et al., “Synthetic-aperture based photoacoustic re-beamforming (spare) approach using beamformed ultrasound data,” Biomed. Opt. Express, 7 (8), 3056 –3068 (2016). BOEICL 2156-7085 Google Scholar


L. Zuoran et al., “The photoacoustic tomography system based on medical ultrasound array,” in Int. Conf. on Photonics and Imaging in Biology and Medicine, W3A–56 (2017). Google Scholar


D. Yang et al., “Fast full-view photoacoustic imaging by combined scanning with a linear transducer array,” Opt. Express, 15 (23), 15566 –15575 (2007). OPEXFF 1094-4087 Google Scholar


T. Chaigne et al., “Light focusing and two-dimensional imaging through scattering media using the photoacoustic transmission matrix with an ultrasound array,” Opt. Lett., 39 (9), 2664 –2667 (2014). OPLEDP 0146-9592 Google Scholar


Y. Wang et al., “Slit-enabled linear-array photoacoustic tomography with near isotropic spatial resolution in three dimensions,” Opt. Lett., 41 (1), 127 –130 (2016). OPLEDP 0146-9592 Google Scholar


J.-F. Synnevag, A. Austeng and S. Holm, “Benefits of minimum-variance beamforming in medical ultrasound imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 56 (9), 1868 –1879 (2009). ITUCER 0885-3010 Google Scholar


R. Paridar et al., “Photoacoustic image formation based on sparse regularization of minimum variance beamformer,” Biomed. Opt. Express, 9 (6), 2544 –2561 (2018). Google Scholar


G. Matrone et al., “The delay multiply and sum beamforming algorithm in ultrasound b-mode medical imaging,” IEEE Trans. Med. Imaging, 34 (4), 940 –949 (2015). ITMID4 0278-0062 Google Scholar


M. Mozaffarzadeh et al., “Double-stage delay multiply and sum beamforming algorithm: application to linear-array photoacoustic imaging,” IEEE Trans. Biomed. Eng., 65 (1), 31 –42 (2018). IEBEAX 0018-9294 Google Scholar


M. Mozaffarzadeh et al., “Double-stage delay multiply and sum beamforming algorithm applied to ultrasound medical imaging,” Ultrasound Med. Biol., 44 (3), 677 –686 (2018). USMBA3 0301-5629 Google Scholar


M. Mozaffarzadeh et al., “Linear-array photoacoustic imaging using minimum variance-based delay multiply and sum adaptive beamforming algorithm,” J. Biomed. Opt., 23 (2), 026002 (2018). JBOPFO 1083-3668 Google Scholar


S. Mehdizadeh et al., “Eigenspace based minimum variance beamforming applied to ultrasound imaging of acoustically hard tissues,” IEEE Trans. Med. Imaging, 31 (10), 1912 –1921 (2012). ITMID4 0278-0062 Google Scholar


X. Zeng et al., “Correspondence-beam-domain Eigenspace-based minimum variance beamformer for medical ultrasound imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 60 (12), 2670 –2676 (2013). ITUCER 0885-3010 Google Scholar


M. Mozaffarzadeh et al., “Eigenspace-based minimum variance adaptive beamformer combined with delay multiply and sum: experimental study,” Proc. SPIE, 10467 1046717 (2018). PSISDG 0277-786X Google Scholar


S.-L. Wang and P.-C. Li, “MVDR-based coherence weighting for high-frame-rate adaptive imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 56 (10), 2097 –2110 (2009). ITUCER 0885-3010 Google Scholar


M. Mozaffarzadeh, M. Mehrmohammadi and B. Makkiabadi, “Image improvement in linear-array photoacoustic imaging using high resolution coherence factor weighting technique,” (2017). Google Scholar


M. Mozaffarzadeh et al., “Enhanced linear-array photoacoustic beamforming using modified coherence factor,” J. Biomed. Opt., 23 (2), 026005 (2018). JBOPFO 1083-3668 Google Scholar


M. Polichetti et al., “A computationally efficient nonlinear beamformer based on p-th root signal compression for enhanced ultrasound b-mode imaging,” in IEEE Int. Ultrasonics Symp. (IUS), 1 –4 (2017). Google Scholar


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). JBOPFO 1083-3668 Google Scholar


K. Sivasubramanian, V. Periyasamy and M. Pramanik, “Hand-held clinical photoacoustic imaging system for real-time non-invasive small animal imaging,” J. Visual. Exp., 128 e56649 (2017). Google Scholar


K. Sivasubramanian and M. Pramanik, “High frame rate photoacoustic imaging at 7000 frames per second using clinical ultrasound system,” Biomed. Opt. Express, 7 (2), 312 –323 (2016). BOEICL 2156-7085 Google Scholar


K. Sivasubramanian et al., “Optimizing light delivery through fiber bundle in photoacoustic imaging with clinical ultrasound system: Monte Carlo simulation and experimental validation,” J. Biomed. Opt., 22 (4), 041008 (2017). JBOPFO 1083-3668 Google Scholar


Laser Institute of America, Orlando, Florida, (2007). Google Scholar


Moein Mozaffarzadeh received his BSc degree in electrical engineering from Babol Noshirvani University of Technology (Mazandaran, Iran), in 2015, and his MSc degree in biomedical engineering from Tarbiat Modares University (Tehran, Iran), in 2017. Currently, he is a research assistant at a research center for biomedical technologies and robotics, institute for advanced medical technologies (Tehran, Iran). His current research interests include photoacoustic image reconstruction, ultrasound beamforming, and biomedical imaging.

Vijitha Periyasamy joined Wayne State University since 2015-May as a PhD student in functional and molecular ultrasound research laboratory ( His interesting areas are medical imaging, object detection, and data mining. He has a strong background in computer science pattern recognition and computer graphics. He holds two bachelors in computer science and post- and telecommunication. He also received a master in computer science from Wayne State University 2017. His current research interests include endocavity ultrasound and photoacoustic for fetal and maternal care.

Manojit Pramanik received his PhD in biomedical engineering from Washington University in St. Louis, Missouri, USA. Currently, he is an assistant professor of the School of Chemical and Biomedical Engineering, Nanyang Technological University, Singapore. His research interests include the development of photoacoustic/thermoacoustic imaging systems, image reconstruction methods, clinical application areas, such as breast cancer imaging, molecular imaging, contrast agent development, and Monte Carlo simulation of light propagation in biological tissue.

Bahador Makkiabadi received his BSc degree in electronics engineering from Shiraz University (Shiraz, Iran), his MSc degree in biomedical engineering from Amirkabir University of Technology (Tehran, Iran), and his PhD in biomedical engineering from the University of Surrey (Guildford, Surrey, UK). Currently, he is working at Research Center for Biomedical Technologies and Robotics (RCBTR), Institute for Advanced Medical Technologies (IAMT), Tehran University of Medical Sciences, Tehran, Iran. His research interests include blind source separation, advanced array signal processing for medical applications, and biomedical imaging.

© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2018/$25.00 © 2018 SPIE
Moein Mozaffarzadeh, Vijitha Periyasamy, Manojit Pramanik, and Bahador Makkiabadi "Efficient nonlinear beamformer based on P’th root of detected signals for linear-array photoacoustic tomography: application to sentinel lymph node imaging," Journal of Biomedical Optics 23(12), 121604 (27 July 2018).
Received: 19 April 2018; Accepted: 13 June 2018; Published: 27 July 2018

Back to Top