Feasibility of identifying reflection artifacts in photoacoustic imaging using two-wavelength excitation.

One of the remaining challenges of bringing photoacoustic imaging to clinics is the occurrence of reflection artifacts. Previously, we proposed a method using multi-wavelength excitation to identify and remove the RAs. However, this method requires at least 3 wavelengths. Here we improve the method further by reducing the required number of wavelengths to 2. We experimentally demonstrate this new method and compare it with the previous one. Results show that this new method holds great feasibility for identifying reflection artifacts in addition to preserving all advantages of the previous method.


Method
In our previous work [4], we imaged the sample with 4 different wavelengths. Extracting information from images corresponding to those wavelengths giving spectral responses of the image features. We assumed that RAs have a similar spectral response as their original real image feature. The similarity was quantified by using the Pearson correlation coefficient between the spectral responses. The Pearson correlation coefficient requires at least 3 data points (3 wavelengths) to give a meaningful coefficient. Therefore, 2 wavelengths are not sufficient with this approach.
In this work, the acquired PA images are first divided by the output energy of the corresponding wavelength. We then obtain the 2-wavelength responses of the image features in the same way as obtaining spectral responses presented in [4]. The 2-wavelength responses are normalized to their maximum. We define slope s of a normalized 2-wavelength response as: s ≡ p(λ 2 ) − p(λ 1 ), where p(λ 1 ) and p(λ 2 ) are the normalized pixel values at 2 wavelengths λ 1 and λ 2 , respectively.
Since the values are normalized to the maximum, either p(λ 1 ) or p(λ 2 ) is 1 resulting in s varying between −1 and 1. Figure 1 shows a normalized 2-wavelength response example of an image feature, the red line. As RAs have a similar spectral response as their original real image feature [4], we assume that the slope of RAs is close to the slope of their real image feature. A threshold, ∆s th , is used to determine the closeness. If an image feature has a slope of s, the slope of its RAs must be in between s − ∆s th and s + ∆s th . The light blue area in Fig. 1 depicts the possible slope of RAs of that image feature. In this work, we use the term "slope difference", ∆s i,j = |s i − s j |, to refer to the difference of the slope of image features i and j. If feature i or feature j is an RA of the other, ∆s i,j <∆s th . On the other hand, since s varies between −1 and 1, ∆s is then between 0 and 2. It is worth reminding that RAs appear at a larger depth and have weaker signal than the original image feature [4]. Combining this with the slope difference, RAs can therefore be identified. Figure 2 shows the flowchart of the method. It is similar to the one presented in [4] except the "Calculate correlation coefficients" step replaced by "Calculate ∆s". Of the 2 acquired images corresponding to the 2 wavelengths, the image having higher signal is selected for segmentation. The segmentation algorithm, based on the Sobel edge detection algorithm, is to detect image features and thus obtain their 2-wavelength responses. These 2-wavelength responses are processed in order to identify RAs.
∆z min in the flowchart is a minimum distance between an RA and its real image feature. This is to avoid the situation when 2 identical absorbers are close to each other, they might have similar 2-wavelength responses, consequently, one might be identified as an RA of the other. In other words, the method only looks for RAs of an image feature in the region below it at least ∆z min . We performed a ∆z min study in our previous work [4], we will apply the same value ∆z min = 1.7 mm in this work. ∆s th is empirically set 0.05 and will be discussed further in the Discussion section.
As reported in [4], the RA identification can also be achieved without using segmentation. In this approach, the slope of each image pixel is compared to all other pixels. Figure 3 shows an example of this approach. Figures 3(a) and 3(b) show 2 PA images of a finger at 808 and 940 nm, respectively. A slope map of all pixels is calculated based on these 2 images, seen in Fig. 3(c). All slopes are then compared to each other pixel by pixel. Further detail of this approach will be described in Section 4.4.

Simulation
In principle, the method is based on the assumption that RAs have a similar 2-wavelength response as their real image feature. In other words, the direct and reflected US waves still carry the optical properties of the absorber along the propagation despite of acoustic attenuation. This assumption may not be valid in cases where the two optical wavelengths generate PA waves in completely different acoustic spectral ranges, and where the broader bandwidth signal suffers significantly more from acoustic attenuation than the narrower bandwidth. To verify this assumption, we performed a simulation using ValoMC [13] which is based on the Monte Carlo and k-Wave toolbox.
We simulated photoacoustic imaging at 2 different wavelengths of a blood vessel above a piece of Perspex in a soft tissue medium. The yellow round and yellow rectangle objects in Fig. 4(a) are the blood vessel with a diameter of 1 mm and the Perspex piece with a thickness of 10 mm. The illumination is vertically from the top center with a beam width of 2 mm and the detection is a line of 64 point detectors at z = 0 mm with a total length of 7.4 mm. As the acoustic attenuation might affect the US waves carrying the optical information of the absorber, especially at high acoustic frequencies, we chose 2 excitation wavelengths to generate signals with different acoustic spectra. The first wavelength is 532 nm. As hemoglobin is strongly absorbing at this wavelength, light can only penetrate into small part of the lumen generating high frequency US waves, see in Fig. 4(b) the initial pressure at 532 nm. The second wavelength is 800 nm to generate lower frequency US waves, see in Fig. 4(c) the initial pressure at 800 nm. We simulated the blood vessel with an oxygen saturation level of 90% which has an absorption coefficient of 233 and 4.3 cm −1 at 532 and 800 nm respectively [14]. The absorption coefficient elsewhere is 0.01 cm −1 . The reduced scattering coefficient of the medium is 12.3 and 8.4 cm −1 at 532 and 800 nm respectively [14] mimicking soft tissue. The blood vessel and the medium have the same speed of sound of 1575 m/s and density of 1065 kg/m 3 [15] while these values in the Perspex piece are 2730 m/s and 1180 kg/m 3 . The acoustic attenuation coefficient in the medium is frequency dependent and modelled as α = af b , where a =0.56 dB/cm/MHz and b = 1 [16].  Figure 4(e) shows the recorded signal at the center detector at the wavelengths of 532 and 800 nm and Fig. 4(f) shows the spectrum of the US waves of the blood vessel and its RA at 532 and 800 nm. The oscillation in the spectrum is similar to the observation in [17]. The spectrum of the signal at the wavelength of 532 nm is broader and thus experiences attenuation differently compared to the signal at the wavelength of 800 nm.
The slopes of the blood vessel's and RA's response are −0.937 and −0.94, respectively. It can be seen that though the reflected US waves experience stronger attenuation due to a longer propagation path (10 mm in this simulation) affecting their acoustic spectrum, they still carry the optical information as in the direct US waves. We repeated the simulation with the Perspex piece repositioned 5 mm deeper resulting in a propagation path of the reflected US waves 20 mm longer than the direct US waves. The slope of the RA's response in this case is −0.934 which is still highly close to the slope of the blood vessel's response. These simulation results strongly enhance our assumption.
We also performed another simulation similar to the one above but without the Perspex layer. A new blood vessel identical to the one above was placed at the position of the RA in the previous simulation. The slopes of the response of the top and the bottom blood vessels are −0.932 and −0.879, respectively. The slope difference is 0.053 and with ∆s th = 0.05, the proposed method still can distinguish these 2 identical absorbers. The slope difference would be even larger for optically different absorbers showing that the method will not mis-identify a real absorber as an RA of another one.
It is worth noting that this simulation did not take into account noise and the spectral response of the detectors and the medium was optically and acoustically homogeneous. In practice, these aspects might affect the 2-wavelength response of the image features.

Experimental results
To experimentally validate the method, we use the same datasets presented in [4]. The data was acquired using a hand-held PAI probe [4]. It consists of laser diodes emitting laser light at 4 different wavelengths 808, 915, 940 and 980 nm with output pulse energy of 0.96, 0.98, 0.89 and 0.82 mJ, respectively. In the experiments. The laser diodes were triggered at a repetition rate of 1 kHz. The transducer array has a center frequency of 7.5 MHz with a bandwidth of 100%. The array contains 128 elements with a pitch of 0.245 mm. The center 64 elements were used in the experiments.
Each dataset contains 4 PA images at the 4 wavelengths and 1 US image averaged over 100 pulses. Of these 4 wavelengths, we select two. In our previous work [4], of all images, the image having the strongest signal was used for segmentation. In this work, the image for segmentation is selected from the images of the 2 wavelengths. Figure 5 shows a phantom experiment. A schematic cross-section and a photo of the phantom are given in Figs. 5(a) and 5(b), respectively. Four black suture wires (USP 4/0, diameter of 0.2 mm, Vetsuture Nylon, The Netherlands) are to mimic four absorbers. We used the same material for the absorbers in order to verify if the method would mis-identify a real absorber as an RA of another absorber. Underneath absorbers 3 and 4 is a Perspex plate with a thickness of 1 mm mimicking an acoustic reflector. The coupling medium was a solution of 3.5% Intralipid 20% in demi-water with an estimated µ s =6 cm −1 at the wavelength of 900 nm [18]. Figure 5(c) shows a PA and US combined image. The gray color part is the US image showing two surfaces of the plate. The hot color part is the PA image at 808 nm. In the PA image, underneath the plate are RAs of the absorbers above the plate.   In the approach using 915 and 940 nm, ∆s 3,5 = 0.0365 and ∆s 4,6 = 0.0492, features 5 and 6 are then identified as an RA of features 3 and 4, respectively. This is in great agreement with the result using 4 wavelengths showing the feasibility of identifying RAs with these two wavelengths. In contrast, in the approach using 940 and 980 nm, most of the slopes are in the light orange area, seen in Fig. 6(c). As a consequence, features 1, 2, 5, 6 are also identified as RAs of feature 3, seen in Fig. 6(d). This result indicates that not all combinations of two wavelengths are appropriate for RA identification.

Phantom
In this phantom study, the image at 940 nm has the highest signal among the selected wavelengths and is used for segmentation in the above approach using 2 wavelengths.

4.2.
In vivo 1 Figure 7 shows the analysis presented in [4]. Figure 7(a) is an acquired PA image showing a cross-sectional view of a finger. The image is then processed for RA removal using 4 wavelengths, seen in Fig. 7(b). Figure 7(c) shows the segmented image with a few features numbered. Features 4, 5 and 6 are RAs of features 2, 3, and 1 respectively [4]. The spectral responses of these features are depicted in Fig. 7(d).   The result of using 2 wavelengths (808 and 940 nm) is similar to the one using 4 wavelengths ( Fig. 7(b)). This enhances the feasibility of identifying RAs using 2 wavelengths shown in the phantom study. However, as mentioned in the phantom study, not all 2-wavelength combinations could work properly. Figure 8(d) shows an example that using 2 wavelengths might miss some RAs.
In this in vivo study, the image at 915 nm is selected for segmentation. However, this image is not included in the approach using 808 and 940 nm. The image at 940 is then used for segmentation in this approach resulting in a slightly different segmented image and feature numbering. This difference does not affect the method. However, for convenience and simplicity of presenting the spectral response, we refer to the same image features, shown in Fig. 7(c), in all approaches.

In vivo 2
In both phantom and in vivo studies presented above, the approach using 808 and 940 nm provides the most reasonable RA identification while any other 2-wavelength combinations either mis-or over-identify RAs. To verify the performance of this combination, we carried out another in vivo experiment with the same protocol as in vivo 1 described in [4]. Figure 9 presents the result of this experiment. Similar to the in vivo 1 study, Fig. 9(a) shows a PA image of a finger. This image is then RA processed using 4 wavelengths ( Fig. 9(b)) or 2 wavelengths (Figs. 9(c)-9(d)). Figures 9(c) and 9(d) are obtained by using 808 and 940 nm and 808 and 980 nm, respectively. In this experiment, the image acquired at 808 nm is used for segmentation, therefore, the result of all approach are presented based on the same image. In the approach using 4 wavelengths, while most of the RAs are removed, one RA, marked with a white arrow, is not. This is probably due to the noise as discussed in [4]. While using 808 and 940 nm works reasonably well in the studies in Sections 4.1 and 4.2, this approach is under performing in this case leaving a few RAs unremoved, seen in Fig. 9(c). However, the approach using 808 and 980 nm provides a corrected image with all RAs identified. It even slightly outperforms the approach using 4 wavelength as the marked RA is also identified. These observations indicate that the approach using 2 wavelengths is strongly affected by the noise and the choice of 2 particular wavelengths is of importance to achieve a reasonable RA identification.

RA identification without segmentation
As described in Section 2, the RA identification can also be achieved without using segmentation. The image processing is the same as presented in [4].   Figure 10(b) shows the spectral slope value of all pixels (see also Fig. 15 in Appendix). From the similarity in spectral slope of the skin and of deeper structures, we can already identify RAs belonging to the skin signal. Figure 10(c) illustrates a pixel (top yellow pixel) in the skin. The slope of the two wavelength spectral response (808 and 940 nm) of this pixel is indicated by the blue line in the white dashed window in Fig. 10(a).
The slope of the reference pixel is compared to other pixels at least 2 mm below it giving a ∆s map, seen in Fig. 10(c) and also in Visualization 1. The red color in Fig. 10(c) indicates ∆s below 0.05 revealing RAs (red pixels) of the considered pixel. Figure 11 shows the results of identifying RAs using 4 (Figs. 11(b) and 11(e)) or 2 (Figs. 11(c) and 11(f)) wavelengths without segmentation. Figures 11(a)   In the phantom experiment, absorbers on the left are partly removed in the method using 4 wavelengths. This is similar to the observation in [4]. The result of 2 wavelengths is slightly better in this aspect. In the in vivo 1 experiment, the results of these 2 methods are in good agreement.
In the approach without using segmentation, most of the noisy pixels at large depths are identified as RAs. This can be also seen in Visualization 1 that while scanning through noisy pixels at the top of the image, noisy pixels at larger depths are identified as RAs of those pixels. At the end of the movie, most of noise pixels are removed. Our hypothesis is that for random noise, there would be always some other noisy pixels having a similar slope. To verify this, we create 2 images (same size as in vivo 1 images) with only random noise (command 'rand" in MATLAB which generates uniformly distributed random numbers), Figs. 12(a) and 12(b), and treat them as 2 PA images. We then calculate a slope map based on these 2 images, seen in Fig. 12(c). Figure 12(d) shows the result with the method without using segmentation. It can be seen that though the 2 initial images contain only noise, most of noisy pixels at large depths are still identified as RAs of their above noise. This result strengthens our hypothesis. It is worth noting that, RAs must have a pixel value smaller than their real absorber. Therefore, real images features at large depths having pixel value higher than the noise level will not be affected by the noisy pixels above them.
We also assess the performance of this method in low signal-to-noise ratio (SNR) images. Figure 13(a) shows again the phantom image with a green and a dashed white rectangle depicting signal and noise windows, respectively. The SNR is defined as: where E(S) and E(N) are the mean square of the pixel values in signal and noise windows, respectively. We manually add random noise to the phantom image using 'rand' in MATLAB, seen in the upper row in Fig. 13. The lower row are RA corrected images. All the images have the same color scale shown at the bottom. We analyze the absorber with the second strongest signal rather than other absorbers with low signals so that we could have a large range of the SNR. At the SNR is 9.4 and lower, RAs of the top right absorber have pixel values close to the noise level. The method can still remove them as RAs of either the top right absorber or above noisy pixels (as discussed above). When reducing the SNR to 5.2 dB and lower, the considered absorber is barely distinguishable from noise, it is then significantly removed. These results show that the method might not work properly at a certain SNR and lower.
To analyze the effect of the SNR on real image features at large depth, we count the number of remaining pixels in the signal window that have a greater pixel value than the noise level determined in the noise window. We repeat this analysis 3 times at each SNR. Figure 14 shows the number of remaining pixels versus the SNR in an error bar plot presenting the average, minimum and maximum values. It can be seen that from the SNR equal and lower than around 6 dB, the real image feature is significantly removed.

Discussion
PAI has great potential in clinical applications. Tackling RAs in PAI is therefore an important step to bring this technique further forwards. The proposed method offers significant advantages over existing methods for identifying RAs. First, the proposed method does not require any training or experience of the user as in [7,9]. Second, it takes advantage of the use of 2-wavelength PAI which is a current focus of PAI [2,11,12]. Third, this method does not need US images to work as in [6][7][8]. Multi-plane wave US imaging comes at a greater time consuming and higher processing expense. Finally, compared to the method using 4 wavelengths [4], the 2-wavelength method is feasible with an even more compact and cost-efficient system.
In this work, we use the ∆z min determined in the previous report [4]. ∆z min defines a region below an original image feature in which the method will not identify any RAs of that feature. This still remains as one of the main limitations of the proposed method.
In the proposed method, ∆s th used to quantify the closeness between the 2-wavelength slopes is a crucial threshold. The larger ∆s th is, the higher the chance to identify RAs the method has. However, it is also a higher chance RA over-misidentification is such as real image features identified as RAs. In this work, ∆s th is set based on our experience to get reasonable results.
Choosing the 2 wavelengths for a reasonable result is of importance as the method must distinguish: (1) different absorbers with the same optical absorption properties located at different locations and (2) different chromophores. For the former, the variance of the local light fluence along the depth (rather than other directions as we use ∆z min ) must be significant enough so that the method will not mis-identify one absorber as an RA of another identical absorber. That means the optical efficient attenuation coefficient, µ eff = 3µ a (µ a + µ s ), of the bulk medium must not be flat at the chosen wavelengths. For the latter, the optical absorption of different chromophores must have different slopes at the chosen wavelengths. These together then require prior knowledge of the imaging sample. However, since the main chromophores in in vivo PAI are melanin and hemoglobin and based on their well-known optical absorption spectra [19], various wavelength combinations can be used. For example, the 2 wavelengths around 800 and 900 nm would be an option. The absorption of melanin in this range has a negative slope while oxygenated hemoglobin and deoxygenated hemoglobin have a positive slope and a relatively flat slope, respectively. Therefore the method could differentiate these chromophores in this range. The µ eff of the bulk tissue mainly depend on the absorption of the blood and scattering the soft tissue. The blood volume fraction in soft tissue is 1-8% [20]. The optical scattering of the soft tissue has a slightly negative slope [14]. If the oxygen saturation level of the blood is 60% or higher, the slope of the µ eff of the bulk medium is always positive with the 2 wavelengths of 800 nm and any wavelength in between 810-900 nm. Therefore the method could distinguish identical absorbers.
It is shown in the experimental result section that in one case, a 2-wavelength combination works properly but another combination does not, and it is opposite in another case. This negatively affects the reliability of the proposed method though it still shows the feasibility of the method for identifying RAs. There would be 2 reasons for this. First, as presented in the experimental results, the slope formed by using 2 wavelengths is highly sensitive to noise. RA mis-or over-identification might happen with the proposed method. This can be seen in Section 4.4 when the SNR of the image is too low. Second, the chosen 2 wavelengths might not differentiate well the optical absorption properties of the different absorbers as discussed in the strategy above. Further improvements in these aspects will be investigated in our future work.
As discussed above, choosing the 2 wavelengths might require prior knowledge about the imaging sample. If the optical properties of the sample change, choosing another 2-wavelength combination might be needed for a reliable RA identification. However, in the wavelength range of 800-900 nm, the change of the absorption of blood at different oxygen saturation levels is little. Using 2 wavelengths in this range would work well for RA identification in in vivo imaging.
As shown in [5], out-of-plane artifacts (with the acoustic source located outside the imaging plane) can also appear in the acquired PA images, these artifacts cannot be tackled using 4 wavelengths. This still holds for the proposed method. An additional method, such as transducer displacement [5,9], is needed to identify these artifacts.
In the approach using segmentation, both methods using 4 or 2 wavelengths take ∼ 0.2-0.3 s time for correcting a PA image with a dimension of 512 × 64 pixels. The processing time in this case is approximately the same since most of the processing time is on the segmentation. However, in the approach without using segmentation, while the method using 4 wavelengths takes ∼ 2 s, the method using 2 wavelengths only takes 0.2-0.3 s. The reason is that calculating the Pearson correlation coefficient is significantly more computationally expensive than comparing the slopes. The processing speed is still slow for real-time imaging however, this is a considerable advantage of the proposed method over the previous one [4].

Conclusion
In this work, we have shown the feasibility of identifying reflection artifacts (RAs) in PAI using 2-wavelength excitation. Results show that a reasonable RA identification can be achieved by using a proper combination of 2 wavelengths. The proposed method offers significant advantages over the previous method using 4 wavelengths [4] in terms of processing time and imaging system compactness and cost.