Open Access
Add to BibSonomy
Abstract
In a photoacoustic (PA) imaging system, the detectors are bandwidth-limited. Therefore, they capture PA signals with some unwanted ripples. This limitation degrades the resolution/contrast and induces sidelobes and artifacts in the reconstructed images along the axial direction. To compensate for the limited bandwidth effect, we present a PA signal restoration algorithm, where a mask is designed to extract the signals at the absorber positions and remove the unwanted ripples. This restoration improves the axial resolution and contrast in the reconstructed image. The restored PA signals can be considered as the input of the conventional reconstruction algorithms (e.g., Delay-and-sum (DAS) and Delay-multiply-and-sum (DMAS)). To compare the performance of the proposed method, DAS and DMAS reconstruction algorithms were performed with both the initial and restored PA signals on numerical and experimental studies (numerical targets, tungsten wires, and human forearm). The results show that, compared with the initial PA signals, the restored PA signals can improve the axial resolution and contrast by 45% and 16.1 dB, respectively, and suppress background artifacts by 80%.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
In the last two decades, photoacoustic imaging (PAI) has been one of the nascent imaging modalities in biomedical imaging [1,2]. The principle of PAI is: 1) A light source irradiates tissue. 2) Based on the tissue absorption coefficient $({\mu _a}),$ which varies at different wavelengths, the tissue absorbs part of the light energy. 3) The absorbed energy warms up the tissue rapidly. 4) This process follows to produce a photoacoustic (PA) wave by the thermoelastic expansion effect. 5) The PA wave can be detected by an ultrasound transducer (UST). By applying proper reconstruction algorithms, a map of the absorption coefficients is obtained from the captured PA signals [3]. Various in-vivo, in-vitro, and ex-vivo applications of PAI demonstrate great potential in clinical applications. As one example, PAI enables high-quality vascular imaging due to the high optical absorption of blood cells [4].
PAI can be implemented in planar, cylindrical, and spherical geometries in different limited-view/full-view arrangements and each has its own clinical/preclinical applications. Cylindrical and spherical geometries can be used for small animal whole-body [5], brain [6], thyroid, carotid [7,8], peripheral arteritis [9], and breast imaging [10]. Planar geometry is appropriate for imaging the skin [11], palm, forearm [12], finger joint, and thigh [13].
In a PAI system, the extension of the point spread function (PSF) along the axial direction depends mainly on the detector bandwidth [14]. If the detector bandwidth is infinite, the PSF extension along the axial direction would be approximated as a Dirac delta function. In practice, for a bandwidth-limited transducer, the PSF extends along the axial direction. However, the extension of the PSF along the lateral direction depends on the solid-angles span of the captured PA signal [14–16]. If the arrays of detectors capture the photoacoustic waves in the whole angular direction (all solid angles are captured), the lateral resolution would be the finest. This condition is met in a full-ring array scenario with a fine pitch [17,18]. In the planar and cylindrical limited-view geometries, due to the limited field-of-view (FOV), all solid angles may not be captured by the detectors. Hence the PSF extends along the lateral direction resulting in a loss of lateral resolution [19]. In a planar PAI system, a real-time reconstruction can be achieved via performing a delay-and-sum (DAS) algorithm [20]. In the reconstruction procedure, the delayed PA signals of all detectors are summed up. As a result, side-lobes and artifacts may appear in the reconstructed image. Furthermore, the above-mentioned blind solid-angles and bandwidth-limited detectors degrade the resolution along both lateral and axial directions.
To overcome the limitations of the DAS algorithm, improved versions called delay-multiply and sum (DMAS) and filtered DMAS (F-DMAS) were proposed [21]. DMAS algorithm multiplies PA signals before summation to improve the image quality in terms of suppressing side-lobes and artifacts. Moreover, an improved version of DMAS, called double stage DMAS (DS-DMAS) is presented in [22]. Adaptive algorithms (e.g., minimum variance (MV)) have also been presented to improve resolution [23]. Although these algorithms improve the resolution along the lateral direction, however, they have no impact on improving the axial resolution. Moreover, Mukaddim and Varghese [24] showed that the MV algorithm does not improve reconstructed image quality especially in in vivo cardiac imaging. The computational cost of the MV algorithm is very high, which prevents its applicability in video-rate imaging applications. Another variation of DAS namely multiple-DAS with enveloping is presented to improve the resolution/contrast [25]. In the presented method in [25], a weighting factor and power three are applied to the output of the DAS algorithm. This power three causes missed data in the reconstructed image, especially in highly detailed areas. The missing data effect occurs when there are small or low-intensity PA absorbers in the ROI. Therefore, they are missed or suppressed in the reconstructed image. Moreover, the limited bandwidth effect of the detectors has not been considered. In recent years, deep learning-based algorithms are also proposed to improve image quality [26,27]. Still, the need for a large number of in-vivo data-set limits their utilization in real applications.
In addition to improving reconstruction algorithms, some weighing factors have been proposed based on DAS, DMAS, and MV algorithms, which can be implemented on the reconstructed image to improve the image quality of ultrasound imaging (USI) and PAI systems. The most prominent one is the coherence factor (CF), which is the ratio of the coherence energy to the total energy of the signals [23,28]. Moreover, phase coherence factor (PCF), sign coherence factor (SCF) [29,30], Hilbert-based coherence factor [31] and signal mean-to-standard-deviation factor (SMSF) [32] were also used to reduce artifacts and improve resolution and image quality. Recently, a delay-and-sum-to-delay-standard-deviation factor is presented to suppress artifacts and improve the lateral resolution of the reconstructed image by the MV algorithm. This weighting factor can be a suitable factor instead of CF [33]. Another fast-weighting factor based on the DAS and DMAS algorithms, named simplified-delay-multiply-and-standard-deviation (SDMASD) factor is presented and accelerated using graphics processing unit (GPU) programming. Appling the SDMASD factor to the output of the DAS and DMAS algorithms will further improve the resolution and contrast of the reconstructed images [34]. One limitation of these weighting-factor-based methods is that they may cause a missing data effect and also are not effective in improving axial resolution. Moreover, some artifact removal factors based on considering the N-shaped PA signal have been proposed for circular-view systems [35,36].
Deconvolution techniques are also presented to compensate for the effect of bandwidth-limited ultrasound detectors by deconvolving the captured PA signal from detectors with a deconvolution kernel. This technique is generally applied iteratively to the captured PA signals from the transducer to compensate for the effect of bandwidth limitation [37]. Although the use of this method in simulations has resulted in excellent results, due to the difficult calculation of the deconvolution kernel in the experiment, they are not satisfactory. Moreover, deconvolution by the Gaussian fitted PSF of the system helps to reach finer lateral resolution. In addition, injection of nano-particles, as contrast agent, enhances the generated PA intensity that helps the small absorbers to be recovered in the reconstructed image [38]. Some frequency compounding methods have been developed to enhance contrast and reduce artifacts in an ultrasound image by averaging at least two uncorrelated sub-band images [39,40]. However, since the attenuation is frequency-dependent, the contrast/resolution cannot be enhanced well, especially in deep tissue.
As stated before, the bandwidth-limited effect of the detectors causes resolution degradation, sidelobe, and artifact appearance in reconstructed PA/US images, especially along the axial direction. In this paper, to compensate for the bandwidth-limited effect of the detectors, a photoacoustic signal restoration algorithm is proposed. The proposed algorithm consists of the following stages; in any captured PA signal, by keeping the locations referring to the position of excited PA absorbers, other weaker ripples of the PA signal are removed. In fact, in an ideal infinite bandwidth structure, the PA signal would be in the form of an N-shaped with no ripples, and the mentioned appeared ripples are due to the limited bandwidth of the detector [41]. This PA signal restoration process results in improving resolution and suppressing artifacts in the reconstructed image. Applying the proposed method to adapted existing reconstruction algorithms and weighting factors, further improves the resolution/contrast of the reconstructed image. The results of considering the restored PA signal instead of the initial PA signal as the input of the DAS and DMAS reconstruction algorithms are represented in the Section 4.
2. Materials and methods
2.1 Effect of a bandwidth-limited transducer
Based on the detector’s properties (low/high cut-off frequency and central frequency), the captured PA signal of a point absorber (delta Dirac function) contains some ripples. Wang et al. [14] expressed the PSF $({\delta (R )} )$ along the axial direction, which is given in (1).
In practice, there is no point absorber, and the PA absorbers have a limited size. In a pulse laser-based photoacoustic tomography (PAT) system, in both thermal and stress confinement conditions, each excited PA absorber will generate a bipolar N-shaped (positive followed by a negative slope edge) PA signal (Fig. 1(a)) [14,42]. In an ideal case (infinite bandwidth detector), the locations between the edges of the N-shaped PA signal (the positive maxima and the negative minimum that comes after) refer to the exact location of the PA absorber. In this case, ripples do not appear in the N-shaped signal. Therefore, a high-resolution image would be reconstructed by performing a reconstruction algorithm. The mentioned limited-size PA absorber can be considered as a set of ideal point absorbers (Fig. 1(b)). Therefore, in a practical case (bandwidth-limited scenario), the captured bandwidth-limited N-shaped signal by the detector can be interpreted as a summation of some PSFs of these sets of point absorbers (Fig. 1(c)), here is called the N-shaped spread function (NSF). For an absorber with a diameter of D, the NSF based on the detector’s high cut-off frequency, would be:
Fig. 1. (a) an N-shaped PA signal originated by a limited size PA absorber with a diameter of D. c is the speed of sound. (b) The N-shaped signal is considered as a set of some delta Dirac functions. (c) The captured bandwidth-limited N-shaped signal by the detector (black solid-line) is a summation of some PSFs of delta Dirac functions in (b). (d) bandwidth-limited N-shaped signal with three different high cut-off frequencies; 6, 8, 11 MHz, and infinite bandwidth. Arrows show the location of the global positive maxima of each high cut-off frequency.
Since the PSF in a bandwidth-limited scenario contains some ripples, these ripples will also appear in the NSF signal, which degrades the image quality and axial resolution. To be more specific, by assuming the NSF as a summation of PSFs, the ripples occurring in each PSF will overlay and affect the main-lobe of other PSFs, and thus the overall shape of the NSF will be changed (see Fig. 1 (c)). In particular, the width of the main-lobe and the location of its global maxima determine the position of the general positive edge. Nevertheless, given that the Dirac functions in the model of the N-shaped signal have different amplitudes, their ripple levels are different (Fig. 1(c)). Specifically, the ripple level of impulses that originates from the center of the N-shaped signal is lower than those that refer to the edges of the N-shaped signal. Therefore, the effect of other ripples, in shifting the maxima PSF which determines the global maxima, is negligible. Based on this, it is expected that the location of the global maxima of the NSF of different band-limited signals (of different sub-band signals of a single band-limited detector) in different frequency bands would be different. This relative shift can be quantitatively calculated by finding the location of the global maxima of each band-limited signal based on (3) and comparing it with the full-bandwidth signal. However, the position of the global maxima of different band-limited signals would be very close (see arrows in Fig. 1(d)). But the spread of the main-lobe is expected to be different. An example of NSF of an excited spherical PA absorber with a diameter of 0.2 mm with different high cut-off frequencies is shown in Fig. 1(d). In this example, high cut-off frequencies are 6 MHz, 8 MHz, and 11 MHz (the bandwidth range of the commercial Verasonics’ transducers; The L11-4 v, 128-element linear array ultrasound probe).
On the other side, the initial captured PA signal has a higher bandwidth than all sub-bands. This higher bandwidth helps to suppress ripple levels and thus low degradation of the original N-shaped signal. Therefore, the location of its global maxima matches the position of the positive edge of the ideal N-shaped signal more accurately. Still, extracting the global positive maxima that refers to the location of the first edge of the absorber from the captured signal (see arrow in Fig. 1(a)) is impossible as there are a lot of local maxima in the captured signal (see arrows in Fig. 1(c)).
In sub-section 2.2, the mentioned features of the sub-band signals are used to extract the approximate location of the global maxima of the NSF from different sub-band signals. This approximate location is then made more accurate by using the location of local maxima in the initial captured PA signal (which has higher bandwidth than all sub-bands). To be more specific, among all local maxima of the initial captured PA signal, a few numbers of them will occur in the extracted approximate region of the PA absorber from the sub-band signals, which results in reducing the number of global maxima candidates. This technique eventually leads to a more accurate estimation of the PA absorber location. Therefore, the proposed method reconstructs the PA absorber with its actual size without any edge extension and unwanted ripples. As a result, improvement in the resolution/contrast compared to the conventional algorithms is achievable.
2.2 PA signal restoration algorithm
As described in sub-section 2.1, our primary goal is to improve the resolution/contrast by restoring the captured PA signals. At first, the initial captured PA signal is divided into several sub-bands. The common region of the positive parts of these sub-bands yields the approximate location of the global positive edge (see arrow in Fig. 1(c)) of the N-shaped signal and thus the PA absorber. The next step is to extract local maxima of the initial PA signal, which occur in this region. Then, calculating the subsequent local minimum and using the edge spread functions (ESF) results in estimating the final location of both edges of each PA absorber. Finally, by maintaining the parts of the signal between noted edges and zeroing others, the restored PA signal is obtained.
The detailed procedure workflow is shown in Fig. 2 and summarized in the following steps:
- 1. Applying sub-band filters
- 2. Extracting common positive region between sub-band signals
- 3. Calculation of the location of the local maxima of the initial PA signal that occurs in the approximate extracted region of step 2.
- 4. Extraction of the subsequent local minimum of each local maxima.
- 5. Consideration of the edge-spread function to extract the location between both edges of each absorber.
- 6. Using step 5 to make a binary mask to be applied to the original signal to remove unwanted ripples and artifacts.
In a typical region of interest, there are different absorbers with different sizes and amplitudes. Therefore, the superposition of generated PA waves from these absorbers is captured by limited-bandwidth detectors. To apply the restoration algorithm, steps 1 to 6 (see Fig. 2) are applied to the captured PA signal by each element in the transducer, separately. Although the related location of each absorber in the time-domain captured PA signal is unknown, by applying steps 2 and 3 to the whole PA signal, the related location of the whole positive edges of all PA absorbers would be extracted simultaneously. Therefore, all locations referring to the positive edges are extracted. Then, steps 4, and 5 should be applied to each extracted positive edge location individually. After that, the superposition of extracted locations in the previous steps will become a binary mask for the whole time-domain PA signal (step 6). Since in the binary mask the value of the locations that refer to all PA absorbers is one, and the other locations’ values are zero, by applying the mask to the initial captured PA signal, the restored signals of all PA absorbers would be obtained simultaneously for each detector. It should be noted that each detector captures the generated PA wave of each absorber in a specific direction. Therefore, the binary mask related to each detector would be different from others and the restoration procedure should be applied to all detectors’ captured signals.
After removing the unwanted parts of the PA signal, the restored PA signal replaces the initial captured PA signal as the input of the reconstruction algorithm. The well-known DAS and DMAS reconstruction algorithms would be good choices to compare the reconstructed image with restored and initial PA signals. A brief overview of the DAS and DMAS is explained in Appendix. The outputs of the DAS reconstruction algorithm based on the initial PA signal (${S_i}$)) and restored PA signal (${S_r}$)) are given in (4a) and (4b).
It is expected that the resolution/contrast of $y_{\textrm{DAS}_\textrm{r}}^{(b )}(. )/y_{\textrm{DMAS}_\textrm{r}}^{(b )}(. )$ would be better than $y_{\textrm{DMAS}_\textrm{i}}^{(b )}(. )/y_{\textrm{DAS}_\textrm{i}}^{(b )}(. ).$ This would be analyzed and discussed via numerical and experimental results in the subsequent sections.
In general, the proposed method is based on the frequency dependent property of the N-shaped PA signal (the dependency of its main lobe spread and ripples on frequency) and also the property of the wide band (the whole spectrum) signal which has a more ideal maxima (with narrower width) but together with a lot of local maxima (ripples). The ripples are excluded by integrating the signals of different sub-bands while the main maxima is made narrower (better axial resolution) based on the wideband signal. This feature is based on the properties of each point absorber and thus a general image with absorbers of different sizes can be considered as containing different groups of absorbers each containing a different number of absorbers. Furthermore, if the intensity level (contrast) of different absorbers is different still this does not affect the output as the sub-bands and wide-band signals are integrated for each specific point-absorber separately. Therefore, the proposed method is object independent, and different excited PA absorbers will not affect the performance of the proposed method.
3. Numerical and experimental study
3.1 Numerical phantom study
The k-Wave MATLAB toolbox [44] was used to generate raw PA data for the numerical study. Two numerical phantoms were synthesized to assay resolution/contrast; one contains nine point absorbers with a diameter of 0.1 mm and an in-between distance of 1 mm along the axial direction called “vertical targets”. The schematic of vertical targets and imaging region is shown in Fig. 4(a). The other synthesized phantom contains two PA absorbers with a diameter of 0.5 mm and five lines with thickness of ∼0.1 mm, which is here called the “hybrid phantom”. The schematic of the hybrid phantom and imaging region is shown in Fig. 5(a). The speed of sound of the medium is 1500 m/s. The ROI was considered 2000 × 1500 grids in the axial and lateral directions with a grid size of 25 µm in both directions. 128 elements with the pitch and kerf size of 12 and 1 grid in the axial and lateral directions respectively, were considered for the arrays of detectors. The properties of Verasonics’ transducers; The L11-4 v were considered as an array of acoustic detectors (central frequency of 6.25 MHz and bandwidth of 70%). A desktop computer with a 64-bit CPU (i9, Intel, USA), 64 GB RAM, and a GPU board (GeForce GTX 2060, NVIDIA, USA) has been used for all computations performed in this work. In our method, like the frequency compounding methods [39,40], the number of sub-bands should be considered in a way to achieve a suitable result. The more sub-bands considered, the better results from restoration algorithm would be obtained. However, as the number of sub-bands increases, the computational cost would also be increased. Moreover, if lots of sub-bands are considered, the improvement of the restoration algorithm would be insignificant. Since at least two sub-bands are required, considering three to four sub-bands (like frequency compounding methods) would be satisfying. In all numerical and experimental studies, three sub-bands are considered. High cut-off frequencies for restoring PA signals are 6 MHz, 8 MHz, and 11 MHz (the bandwidth range of the commercial Verasonics’ transducers; The L11-4 v, 128-element linear array ultrasound probe).
Fig. 3. Schematic of the PAI experimental setup; FB, fiber bundle; L, laser; DAQ, data acquisition; USP, ultrasound probe; WT, water tank; PC, personal computer.
Fig. 4. The numerical results of the vertical targets. (a1) the schematic of vertical targets and imaging region. The light gray box indicates the position of the ultrasound probe. (a2) The maginified region of white box in (a1). (a3) indicates the regions outlined whitin the green box in (a2). The white dashed line in (a3) indicate the locations of CR calculation. (b1, c1) are the results of DAS and DMAS based on the initial PA signal, respectively. (d1, e1) are the results of DAS and DMAS based on restored PA signal, respectively. The pixels outlined in green boxes in (b1-e1) are magnified in (b2-e2). (f) One-dimensional profile along white dashed lines in (b2-e2).
Fig. 5. The numerical results of the hybrid phantom. (a) the schematic of hybrid phantom and imaging region. The light gray box indicates the position of the probe. (b, d) are the results of DAS and DMAS based on the initial PA signal, respectively. (c, e) are the results of DAS and DMAS based on restored PA signal, respectively.
3.2 Experimental phantom and in-vivo studies
Two experimental studies were performed; 1) a lab-made tungsten wires phantom and 2) a forearm of a 30-years old man (see Fig. 7(a)). In lab-made tungsten wires phantom, six tungsten wires were placed as PA absorbers in the optical path of the imaging system (see ROI Fig. 3). The schematic of the experimental setup is shown in Fig. 3. An Nd: YAG pulsed laser with a pulse repetition frequency (PRF) of 20 Hz and pulse width of 6 ns at 650 nm wavelength was used as the excitation source [45]. The laser was guided by a one-to-two fiber bundle. These fiber bundles can excite an ROI with a diameter of 20 mm, approximately.
3.3 Figures of merit
For quantitative performance evaluation of numerical/ experimental studies, we calculated some image metrics, such as full-width at half maxima (FWHM), contrast ratio (CR), contrast-to-noise ratio (CNR), peak signal to noise ratio (PSNR), mean of the standard deviation of the back-ground levels (STD) and structural similarity index measure (SSIM). FWHM, CR, PSNR and SSIM were measured in numerical and phantom experiments. The CNR was measured in the in-vivo study, where the ground truth is not available.
CR and CNR (in dB) are defined as follows:
PSNR (in dB) is expressed as follows:
where, $MSE = \mathop \sum \nolimits_{M,N} \frac{{{{[{{I_G}({M,N} )- {I_R}({M,N} )} ]}^2}}}{{MN}},$ R is the maxima pixel value of the image, M and N are the numbers of pixels in the reconstructed image along the axial and lateral directions, respectively. ${I_G}({M,N} )$ and ${I_R}({M,N} )$ are the intensity of the ground truth and the reconstructed image, respectively. Higher the PSNR value represents a less amount of artifacts in the reconstructed image.4. Numerical and experimental results
4.1 Numerical results of the vertical targets
The synthesized vertical targets are reconstructed with both DAS and DMAS algorithms. The DAS results of the reconstruction with the initial PA signal (DASi) and the restored PA signal (DASr) are shown in Fig. 4(b1, b2) and (c1, c2), respectively. Also, Fig. 4(d1, d2) and (e1, e2) are the DMAS results of the initial PA signal (DMASi) and the restored PA signal (DMASr), respectively. To compare the impact of the restoration process on improving resolution/contrast in different reconstruction algorithms, the region outlined in the white box in Fig. 4(a1) is magnified in Fig. 4(a2) and the reconstructed image whitin white box are shown in Fig. 4(b1-e1). Also, the regions outlined in the green boxes in Fig. 4(b1-e1) are magnified in Fig. 4(b2-e2). From Fig. 4(b-e), it is obvious that the reconstructed image with the restored PA signals has a finer resolution compared to the initial signal in both DAS and DMAS algorithms (see arrows in Fig. 4(b-e)). Quantitatively, the PSNR values calculated for pixels are outlined in the green boxes of Fig. 4(b1-e1). The PSNR values for DASi, DASr, DMASi, and DMASr, are 25.3 dB, 28.2 dB, 29.1 dB, and 32.9 dB, respectively. For the CR and FWHM calculations, one-dimensional profiles along the axial direction pass through the middle of the PA absorbers of Fig. 4(b2-e2) are plotted in Fig. 4(f). By comparing the sidelobe levels of Fig. 4(f) in different scenarios, it can be concluded that the DMAS reconstruction algorithm with the initial PA signals, does not suppress sidelobes along the axial direction (see red circle in Fig. 4(f)). While by using restored PA signals in both DAS and DMAS algorithms, sidelobes suppression is reached (compare arrows and red oval in Fig. 4(f)). The CR values of all different scenarios are calculated based on (5). The average FWHM of all targets in Fig. 4(b2-e2), CR values of one-dimentional profiles (see white dashed lines in Fig. 4(a3, b2-e2), and the PSNR values of Fig. 4(b2-e2) are listed in Table 1.
Table 1. The measured quantitative values of different phantoms
4.2 Numerical results of the hybrid phantom
The results of the hybrid phantom that contain two adjacent point targets and five solid lines are shown in Fig. 5 in dB scale. The DAS results show that the reconstruction with the restored PA signals has fewer artifacts than the initial PA signals (compare white arrows and white dashed circles in Fig. 5(b, c)) and the same is true for the DMAS algorithm (compare white arrows and dashed circles in Fig. 5(d, e)). The SSIM values in all scenarios (Fig. 5(b-e)) are listed in Table 1. For better comparison in terms of contrast of all scenarios, the CR values of a one-dimensional profile along the axial direction passing through the middle of the left-sided PA absorbers (see blue lines in Fig. 5(b-e)) are calculated. The SSIM values of the DAS and DMAS algorithms with restored PA signals reach 0.77 and 0.96, respectively which results in 150% and 81% improvement compared to DAS and DMAS algorithm with initial PA signals, respectively. Moreover, the restored PA signals have improved CR values by 8.7 dB and 16.1 dB compared to the initial PA signals reconstructed with DAS and DMAS algorithms (see Table 1 for more detail).
4.3 Experimental results of the tungsten wires phantom
The reconstructed images of the tungsten wires phantom of different scenarios DASi, DASr, DMASi, and DMASr are shown in Fig. 6(a1-d1) in the dB scale, respectively. For a closer look and comparing the results of using the restored PA signals and initial PA signals, the areas within the blue and green boxes in Fig. 6(a1-d1) are magnified in Fig. 6(a2-d2) and Fig. 6(a3-d3), respectively. The restored PA signals with DAS result in a reconstructed image have lower sidelobes and artifacts than the initial PA signals (see and compare arrows and dashed circles in Fig. 6(a1, b1), Fig. 6(a2, b2) and Fig. 6(a3, b3)). The initial PA signals with DMAS suppress more artifacts than the initial signals with DAS. By using the restored PA signals instead of the initial PA signals in DMAS, more artifacts and sidelobe levels are suppressed in the reconstructed image (compare arrows and dashed circle in Fig. 6(c1, d1), Fig. 6(c2, d2) and Fig. 6(c3, d3)). The PSNR values of Fig. 6(a3-d3) are 26.3 dB, 27.4 dB, 30 dB, and 31.8 dB, respectively. To evaluate the performance of the restoration procedure, the intensities of pixels along the vertical line passing through the middle of the PA absorbers of green boxes in Fig. 6(a1-d1) are plotted in Fig. 6(e). The blue area in Fig. 6(e) is magnified in Fig. 6(f). There is no major difference in terms of resolution (width of -6 dB intensity) in the results of the initial signal with both DAS and DMAS. But using the restored PA signals, instead of the initial PA signals in both DAS and DMAS algorithms, improves the resolution (Compare the green dashed circles in Fig. 6(i)). The mean values of FWHM and CR values for all four points are listed in Table 1.
Fig. 6. The experimental results of the tungsten wires phantom. (a1, c1) are the results of DAS and DMAS based on the initial PA signal, respectively. (b1, d1) are the results of DAS and DMAS based on restored PA signal, respectively. The pixels outlined in blue and green boxes in (a1-d1) are magnified in (a2-d2) and (a3-d3), respectively. (e) One-dimensional profile along the axial direction passing through the middle of the PA absorbers in green boxes. The values outlined in a light blue box in (e) are magnified in (f).
4.4 In-vivo results of the human forearm
To evaluate the performance of restored PA signals in suppressing artifacts and improving the image quality, the forearm of a 30-years old man was imaged (Fig. 7(a)). Figure 7(b) shows the dominant map of a human hand's arteries and nerves. In the forearm, the radial artery (RA) and ulnar artery (UA) are the main arteries. Therefore, those are the main PA absorbers in the forearm. The blue rectangle in Fig. 7(b) and the dashed line in Fig. 7(a) indicate the position of the ultrasound probe to capture PA waves in the in-vivo study. Based on the position of the transducer, it is expected that the UA and maybe some melanocytes are restored well in the reconstructed image. The experimental results of different scenarios are shown in Fig. 7(c-f). In the results, two primary PA absorbers; the ulnar artery (green dashed rectangle in Fig. 7(c-f)) and a melanocyte (see white arrow in Fig. 7(b)), are detectable via DASi, DASr, DMASi, and DMASr methods. In the DAS reconstruction, the reconstructed image with the restored PA signals has much fewer artifacts in the upper part of the skin (see the light blue arrows in Fig. 7(c, d)) and areas near the ulnar artery (compare dashed ovals and rectangles in Fig. 7(c, d)). Also, artifacts within these regions (ovals and rectangles) are suppressed more in DMASr compared to DMASi (compare areas within the dashed ovals and rectangles in Fig. 7(e, f)) because the primary PA absorbers (the ulnar artery and melanocyte) are completely recovered in the DASr and DMASr. Therefore, the results of the human forearm demonstrate that the restoration of the PA signals suppresses artifacts and improves the image quality without missing any data. The average of the STD values measured for the two regions outlined with dashed lines (white dashed rectangle and oval regions) are 0.0080, 0.0034, 0.0051, and 0.0011 for DASi, DASr, DMASi, and DMASr, respectively. The results show that the reconstructed image by DMASr reduces the mean of the STD values by ∼78% and ∼86% compared to DMASi and DASi, respectively. Moreover, the CNR values of the DAS and DMAS algorithms with restored PA signals reach 55.9 dB and 49.4 dB, respectively which results in ∼ 11 dB and 7.5 dB improvement compared to DAS and DMAS algorithm with initial PA signals, respectively (see Table 1 for more detail).
Fig. 7. The experimental results of the right forearm of a 30-years old human hand. (a) a photograph of the right forearm. The purple dashed line indicates the location of the imaging region cross-section. (b) the dominant map of a human hand's arteries and nerves; UN, ulnar nerve, UA, ulnar artery; RA, radial artery. The blue rectangle indicate the position of the linear array probe. (c, d) are the results of DAS and DMAS based on the initial PA signal, respectively. (e, f) are the results of DAS and DMAS based on restored PA signal, respectively.
5. Discussion
In a planar PAI system, because of the blind solid-angles and bandwidth-limited detectors, reconstruction of an image with high contrast and resolutions is challenging. In recent years, to overcome these limitations, some methods have been presented. However, most of these methods only improve the lateral resolution and suppress artifacts and those did not consider the bandwidth-limited effect. Also, some frequency compounding methods have been developed for ultrasound image quality improvement. However, since the attenuation is frequency-dependent, the contrast/resolution cannot be enhanced well. Here, we implement a PA signal restoration algorithm to compensate for the bandwidth-limited effect. After performing the restoration procedure, the restored PA signal is considered as the input of the reconstruction algorithm instead of the initial PA signals. On the other hand, the proposed algorithm is a pre-processing technique to remove unwanted parts from PA signals to achieve a reconstructed image with high axial resolution/contrast. Since the proposed method is developed based on the properties of point absorbers, each absorber in an imaging region can be considered as a region containing several numbers of point absorbers with varying intensities. Therefore, the proposed method is robust, and a real sample containing different absorbers in terms of size and amplitude, will not affect the performance of the proposed method. In the proposed PA signal restoration algorithm, we find the position of absorbers in each PA signal with a binary mask and zero other parts of the PA signal that refer to the ripples and unwanted parts which have of no information. Therefore by zeroing these unwanted parts, the restored PA signal will be filtered which results in higher SNR/CNR values compared to the initial PA signal.
In the conventional method, the initial PA signal was considered as the IQ demodulation of the RF data. But in the restoration algorithm, the location of the binary mask is calculated based on the real parts of the initial signal. Finally, to display a smooth and unipolar reconstructed image, absolute values are taken after reconstructing the image.
Based on the numerical/experimental results for a PA absorber with a diameter of 100 µm, the axial resolution reaches ∼120 µm via the proposed approach. The PA restoration algorithm improves the resolution by ∼ 45% and CR up to 16.1 dB. Moreover, restored PA signals improves the SSIM value ∼150% and 81% in the DAS and DMAS algorithms, respectively. Besides the advantage of the proposed method in improving the axial resolution/contrast, the mean of the STD/CNR values of the human forearm indicates that the DMAS algorithm with restored PA signals suppresses background artifacts by ∼78%/11 dB and ∼80%/7.5 dB compared to DMAS and DAS algorithm with initial PA signals, respectively.
Since the proposed method provides high axial resolution/contrast, in the future, we aim to detect embolism and blood clots in clinical applications where accidental injuries occur. Moreover, by applying the proposed method to multi-spectral photoacoustic imaging, a more accurate and high-contrast map of chromophores can be obtained. Given that the main goal of this paper is PA signal restoration, applying the proposed method to other reconstruction algorithms (e.g., MV, Model-based methods [46,47]) and weighting factors (e.g., CF and SMSF) may further improve the quality of the reconstructed image.
6. Conclusion
In this work, we implemented a photoacoustic signal restoration algorithm to improve resolution and contrast. The numerical and experimental results have shown that considering the restored PA signals instead of the initial PA signals improves resolution/contrast and suppresses background artifacts in the reconstructed image. Due to the simplicity of the restoration algorithm, the restored PA signal would be replaced by the initial PA signal in any reconstruction algorithms.
Appendix
Delay-and-sum (DAS) algorithm
In a DAS algorithm, the captured signals by the detectors are delayed and summed to generate the output refering to the considered pixel. The output of the DAS algorithm for pixel $({x,y} )$ is define as follows:
where, N is the receiving aperture size, ${S_k}$ is the PA signal captured by k-th element with the corresponding time delay, $\Delta {t_{k,x,y}}$.Delay-multiply-and-sum (DMAS) algorithm
In DMAS algorithm, to generate the output refering to the considered pixel, the delayed signals are combinatorially multiplied with each other and summed. The reconstructed image with DMAS algorithm, has higher contrast and finer lateral resolution compared to DAS and for each pixel $({x,y} )\,$ is defined as follows:
Funding
Iran National Science Foundation (97018266).
Acknowledgments
This work was supported by the Iran Natural Sciences Fundation (Grant No. 97018266).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. D. Das, A. Sharma, M. Rajendran, and Pramanik, “Another decade of photoacoustic imaging,” Phys. Med. Biol. 66(5), 05TR01 (2021). [CrossRef]
2. A. B. E. Attia, G. Balasundaram, M. Moothanchery, U. S. Dinish, R. Bi, V. Ntziachristos, and M. Olivia, “A review of clinical photoacoustic imaging: Current and future trends,” Photoacoustics 16, 100144 (2019). [CrossRef]
3. S. Manohar and D. Razansky, “Photoacoustics: a historical review,” Adv. Opt. Photonics 8(4), 586 (2016). [CrossRef]
4. M. Heijblom, W. Steenbergen, and S. Manohar, “Clinical Photoacoustic Breast Imaging: The Twente experience,” IEEE Pulse 6(3), 42–46 (2015). [CrossRef]
5. H.-P. Brecht, R. Su, M. Fronheiser, S. Ermilov, A. Conjusteau, and A. Oraevsky, “Whole-body three-dimensional optoacoustic tomography system for small animals,” J. Biomed. Opt. 14(06), 1 (2009). [CrossRef]
6. L. Nie, Z. Guo, and L. V. Wang, “Photoacoustic tomography of monkey brain using virtual point ultrasonic transducers,” J. Biomed. Opt. 16(7), 076005 (2011). [CrossRef]
7. A. Dima and V. Ntziachristos, “In-vivo handheld optoacoustic tomography of the human thyroid,” Photoacoustics 4(2), 65–69 (2016). [CrossRef]
8. M. Yang, L. Zhao, X. He, N. Su, C. Zhao, H. Tang, W. Li, F. Yang, L. Lin, B. Zhang, R. Zhang, Y. Jiang, and C. Li, “Photoacoustic/ultrasound dual imaging of human thyroid cancers: an initial clinical study,” Biomed. Opt. Express 8(7), 3449–3457 (2017). [CrossRef]
9. S. Ermilov, R. Su, M. Zamora, T. Hernandez, V. Nadvoretsky, and A. Oraevsky, Optoacoustic angiography of peripheral vasculature (SPIE BiOS). (SPIE, 2012).
10. A. P. Rao, N. Bokde, and S. Sinha, “Photoacoustic Imaging for Management of Breast Cancer: A Literature Review and Future Perspectives,” Appl. Sci. 10(3), 767 (2020). [CrossRef]
11. X. L. Dean-Ben and D. Razansky, “Optoacoustic imaging of the skin,” Exp. Dermatol. 30, 1598–1609 (2021). [CrossRef]
12. J. Yang, G. Zhang, W. Chang, Z. Chi, Q. Shang, M. Wu, T Pan, L. Huang, and H. Jiang, “Photoacoustic imaging of hemodynamic changes in forearm skeletal muscle during cuff occlusion,” Biomed. Opt. Express 11(8), 4560–4570 (2020). [CrossRef]
13. T. Vu, D. Razansky, and J. Yao, “Listening to tissues with new light: recent technological advances in photoacoustic imaging,” J. Opt. 21(10), 103001 (2019). [CrossRef]
14. M. Xu and L. V. Wang, “Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction,” Phys. Rev. E 67(5), 056605 (2003). [CrossRef]
15. S. Hakakzadeh, M. Mozaffarzadeh, S. M. Mostafavi, M. Amjadian, Z. Kavehvash, M. Verweij, and N. D. Jong, “Finite Transducer Size Compensation in Two-Dimensional Photoacoustic Computed Tomography,” in 2021 IEEE International Ultrasonics Symposium (IUS), 11-16 Sept. 2021 2021, 1–3.
16. S. Hakakzadeh and Z. Kavehvash, “Image Quality Equations for Focused Transducer in Circular Photoacoustic Computed Tomography,” in 2021 29th Iranian Conference on Electrical Engineering (ICEE), 18-20 May 2021 2021, 917–921.
17. S. Hakakzadeh, M. Mozaffarzadeh, S. M. Mostafavi, Z. Kavehvash, M. Rajendran, N. D. Verweij, M. Jong, and Pramanik, “Multi-angle data acquisition to compensate transducer finite size in photoacoustic tomography,” Photoacoustics 27, 100373 (2022). [CrossRef]
18. Y. Zhang, S. Yang, Z. Xia, R. Hou, B. Xu, L. Hou, J. H. Marsh, J. J. Hou, S. M. R. Sani, X. Liu, and J. Xiong, “Co-optimization method to improve lateral resolution in photoacoustic computed tomography,” Biomed. Opt. Express 13(9), 4621 (2022). [CrossRef]
19. S. Hakakzadeh and Z. Kavehvash, “Blind Angle and Angular Range Detection in Planar and Limited-View Geometries for Photoacoustic Tomography,” in 2021 29th Iranian Conference on Electrical Engineering (ICEE), 2021: IEEE, 922–926.
20. S. Hakakzadeh, S. M. Mostafavi, M. Amjadian, and Z. Kavehvash, “Fast Adapted Delay and Sum Reconstruction Algorithm in Circular Photoacoustic Tomography,” in 2021 28th National and 6th International Iranian Conference on Biomedical Engineering (ICBME), (2021).
21. G. Matrone, A. S. Savoia, G. Caliano, and G. Magenes, “The Delay Multiply and Sum Beamforming Algorithm in Ultrasound B-Mode Medical Imaging,” IEEE Trans. Med. Imaging 34(4), 940–949 (2015). [CrossRef]
22. M. Mozaffarzadeh, A. Mahloojifar, M. Orooji, S. Adabi, and M. Nasiriavanaki, “Double-Stage Delay Multiply and Sum Beamforming Algorithm: Application to Linear-Array Photoacoustic Imaging,” IEEE Trans. Biomed. Eng. 65(1), 31–42 (2018). [CrossRef]
23. S. Park, A. B. Karpiouk, S. R. Aglyamov, and S. Y. Emelianov, “Adaptive beamforming for photoacoustic imaging,” Opt. Lett. 33(12), 1291 (2008). [CrossRef]
24. R. A. Mukaddim and T. Varghese, “Spatiotemporal Coherence Weighting for In Vivo Cardiac Photoacoustic Image Beamformation,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 68(3), 586–598 (2021). [CrossRef]
25. X. Ma, C. Peng, J. Yuan, Q. Cheng, G. Xu, X. Wang, and P.L. Carson, “Multiple Delay and Sum With Enveloping Beamforming Algorithm for Photoacoustic Imaging,” IEEE Trans. Med. Imaging 39(6), 1812–1821 (2020). [CrossRef]
26. J. Kyong Hwan, M. T. McCann, E. Froustey, and M. Unser, “Deep Convolutional Neural Network for Inverse Problems in Imaging,” IEEE Trans. on Image Process. 26(9), 4509–4522 (2017). [CrossRef]
27. S. Jeon, W. Choi, B. Park, and C. Kim, “A Deep Learning-Based Model That Reduces Speed of Sound Aberrations for Improved In Vivo Photoacoustic Imaging,” IEEE Trans. on Image Process. 30, 8773–8784 (2021). [CrossRef]
28. S. Hakakzadeh, S. M. Mostafavi, M. Amjadian, and Z. Kavehvash, “Adapted Coherent Weighting in Photoacoustic Tomography,” in 2021 28th National and 6th International Iranian Conference on Biomedical Engineering (ICBME), 25-26 Nov. 2021 2021, 28–32, doi: 10.1109/ICBME54433.2021.9750310.
29. Y. Zhang, Y. Guo, and W. Lee, “Ultrafast Ultrasound Imaging Using Combined Transmissions With Cross-Coherence-Based Reconstruction,” IEEE Trans. Med. Imaging 37(2), 337–348 (2018). [CrossRef]
30. J. Camacho, M. Parrilla, and C. Fritsch, “Phase Coherence Imaging,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 56(5), 958–974 (2009). [CrossRef]
31. S. Hakakzadeh and Z. Kavehvash, “A Hilbert-based Coherence Factor for Photoacoustic Imaging,” IEEE (2022).
32. Y. Wang, C. Zheng, H. Peng, and Q. Chen, “An adaptive beamforming method for ultrasound imaging based on the mean-to-standard-deviation factor,” Ultrasonics 90, 32–41 (2018). [CrossRef]
33. S. Paul, A. Thomas, and M. S. Singh, “Delay-and-sum-to-delay-standard-deviation factor: a promising adaptive beamformer,” Opt. Lett. 46(18), 4662–4665 (2021). [CrossRef]
34. S. Paul, S. Mulani, N. Daimary, and M. S. Singh, “Simplified-Delay-Multiply-and-Sum-Based Promising Beamformer for Real-Time Photoacoustic Imaging,” IEEE Trans. Instrum. Meas. 71, 1–9 (2022). [CrossRef]
35. S. Hakakzadeh, Z. Kavehvash, and M. Pramanik, “Artifact Removal Factor for Circular-view Photoacoustic Tomography,” in 2022 IEEE International Ultrasonics Symposium (IUS), 10-13 Oct. 2022 2022, 1–4.
36. S. Hakakzadeh, Z. Rajendran, M. Kavehvash, and Pramanik, “A Spatial-domain Factor for Sparse-sampling Circular-view Photoacoustic Tomography,” IEEE JSTQE 1, 1–4 (2023). [CrossRef]
37. J. Prakash, A. S. Raju, C. B. Shaw, M. Pramanik, and P. K. Yalavarthy, “Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography,” Biomed. Opt. Express 5(5), 1363–1377 (2014). [CrossRef]
38. R. S R, K. Chitre, A. Kurup, U. Nongthomba, S. Murty Srinivasula, and M. Suheshkumar Singh, “Simultaneous multiple-level magnification selective plane illumination microscopy (sMx-SPIM) imaging system,” J. Opt. 24(2), 024010 (2022). [CrossRef]
39. J. H. Chang, H. H. Kim, J. Lee, and K. K. Shung, “Frequency compounded imaging with a high-frequency dual element transducer,” Ultrasonics 50(4-5), 453 (2010). [CrossRef]
40. C. Yoon, G.-D. Kim, Y. Yoo, T.-K. Song, and J. H. Chang, “Frequency equalized compounding for effective speckle reduction in medical ultrasound imaging,” Biomedical Signal Processing and Control 8(6), 876–887 (2013). [CrossRef]
41. S. Hakakzadeh, S. M. Mostafavi, and Z. Kavehvash, “Unipolar Back-Projection Algorithm for Photoacoustic Tomography,” in 2022 IEEE International Ultrasonics Symposium (IUS), 10-13 Oct. 2022 2022, 1–4.
42. C. Tian, M. Pei, K. Shen, S. Liu, Z. Hu, and T. Feng, “Impact of System Factors on the Performance of Photoacoustic Tomography Scanners,” Phys. Rev. Appl. 13(1), 014001 (2020). [CrossRef]
43. L. V. Wang and H.-I. Wu, Biomedical optics: principles and imaging, (John Wiley & Sons, 2012).
44. N. A. Rejesh, H. Pullagurla, and M. Pramanik, “Deconvolution-based deblurring of reconstructed images in photoacoustic/thermoacoustic tomography,” J. Opt. Soc. Am. A 30(10), 1994–2001 (2013). [CrossRef]
45. Y. Zhang, Y. Wang, L. Lai, and Wang, “Video-rate Dual-modal Wide-beam Harmonic Ultrasound and Photoacoustic Computed Tomography,” IEEE Trans. Med. Imaging. 41(3), 727–736 (2022). [CrossRef]
46. M. A. Araque Caballero, J. Gateau, X.-L. Dean-Ben, and V. Ntziachristos, “Model-Based Optoacoustic Image Reconstruction of Large Three-Dimensional Tomographic Datasets Acquired With an Array of Directional Detectors,” IEEE Trans. Med. Imaging 33(2), 433–443 (2014). [CrossRef]
47. L. Ding, D. Razansky, and X. L. Dean-Ben, “Model-Based Reconstruction of Large Three-Dimensional Optoacoustic Datasets,” IEEE Trans. Med. Imaging 39(9), 2931–2940 (2020). [CrossRef]
