Reflection artifact identification in photoacoustic imaging using multi-wavelength excitation

: Photoacoustic imaging has been a focus of research for clinical applications owing to its ability for deep visualization with optical absorption contrast. However, there are various technical challenges remaining for this technique to find its place in clinics. One of the challenges is the occurrence of reflection artifacts. The reflection artifacts may lead to image misinterpretation. Here we propose a new method using multiple wavelengths for identifying and removing the reflection artifacts. By imaging the sample with multiple wavelengths, the spectral response of the features in the photoacoustic image is obtained. We assume that the spectral response of the reflection artifact is better correlated with the proper image feature of its corresponding absorber than with other features in the image. Based on this, the reflection artifacts can be identified and removed. Here, we experimentally demonstrated the potential of this method for real-time identification and correction of reflection artifacts in photoacoustic images in phantoms as well as in vivo using a handheld photoacoustic imaging probe.


Introduction
In the last decade, significant progress has been made for translating photoacoustic imaging (PAI) into clinics [1]. This technique uses the photoacoustic (PA) effect, where materials absorb short pulsed light and generate ultrasound (US) waves. The US waves can be detected using US transducers for reconstructing the absorbing structures. Since in tissue the US waves experience order of magnitude less scattering compared to light, much deeper information can be reconstructed compared to purely optical imaging techniques. Therefore, PAI provides optical absorption contrast and has the ability to image deeper than purely optical imaging techniques at ultrasonic resolution. Exploiting these properties, current research is focusing on investigating its clinical applications such as imaging of breast cancer [2,3], rheumatoid arthritis [4,5], and atherosclerosis [6]. Additionally, multispectral PAI strengthens the advantages of this technique for screening and monitoring human diseases, for instance by examining the oxygen saturation of hemoglobin (sO 2 ) in the lesion [3] or characterizing different tissues [6].
Recent research has focused on developing compact and affordable PAI systems. Integrating the laser source, especially a low cost laser source [7], into commercial handheld US probes for clinical use of PAI systems was proposed [7][8][9][10][11][12]. However, one of the major drawbacks of using a linear US transducer array is the occurrence of reflection artifacts (RAs) due to its limited view angle. As photoacoustically generated US waves propagate in all directions, the US waves propagating away from the US transducer array can be reflected towards the US transducer array by acoustic heterogeneities such as bone and tendon causing RAs in the acquired photoacoustic image. The RAs appear at larger depths than the real absorbers leading to misinterpretation of the acquired images. For clinical usage, real-time correction of the RAs in PAI is of fundamental importance.
RAs are also called in-plane artifacts, one of two types of artifacts (clutter) in PAI. The other type is out-of-plane artifact [13]. The term "plane" represents the imaging plane defined by the US transducer array. Since the laser beam excites a large volume, absorbers which are not in the imaging plane absorb the light and generate signals. If the out-of-plane sensitivity of the transducers is high enough, these absorbers appear in the acquired image resulting in out-of-plane artifacts (direct out-of-plane artifacts). If there is an acoustic reflector located underneath these out-of-plane absorbers, out-of-plane RAs (indirect out-of-plane artifacts) can be present in the acquired image [13]. In this work, we aim to tackle RAs (in-plane artifacts).
Several methods for reducing RAs have been presented [14][15][16][17][18]. Deformation Compensated Averaging (DCA) [14] employs tissue deformation for de-correlating the artifact by slightly palpating the tissue. This technique requires a well-trained person, sufficient deformation of tissue and works for easily deformable tissue. Localized vibration tagging (LOVIT) [15], introduced by Jaeger, uses a similar principle as DCA but using the acoustic radiation force (ARF) aimed at the artifact in the focal region of the ultrasonic beam instead of tissue palpation. This is a promising approach to overcome the disadvantages of DCA. However, it can only reduce artifacts based on the deformation of tissue in the US focal region. This limits the real-time capability and has safety challenges. Recently, LOVIT has been further improved by using multiple foci [19]. Another method exploits the acoustic tissue information by inversion of a linear scatter model using plane wave US measurements [16]. This method has to match PA and US measurements and requires numerous plane wave angles limiting itself to real-time performance. Schwab then introduced an advanced interpolation approach to significantly reduce the required number of plane waves in a linear scattering medium [20]. Allman introduced a convolutional neural network to remove RAs of point-like sources with high accuracy [18]. However, since the network is trained with simulated data, the accuracy might be negatively affected in in vivo situations.
Previously Singh introduced a method, photoacoustic-guided focused ultrasound (PAFUSion), using focused ultrasound or synthetic backpropagation to mimic PA sources and thus identify the RAs [17,21,22]. This method can efficiently reduce the RAs, however it has several limitations: mimicking the PA source is limited by the angular aperture of the US probe; numerous additional US images are needed, challenging real-time artifact reduction; the PA sources (skin, blood vessels) must be perpendicular to the imaging plane that requires demanding alignment effort; the PA signal from the source and the mimicked signal by US must match each other in terms of amplitude and frequency content which might negatively affect the accuracy of the method.
In this paper, we propose a new method where we exploit the use of multispectral PAI for identifying and removing RAs. Imaging with multiple wavelengths, PA spectral responses of the features in the acquired image can be obtained. Our method is based on the assumption that RAs are better correlated with the image features of their corresponding original absorbers than with other features, exposing the suspicious artifacts. In addition, RAs appear at larger depths and have weaker signals than the original image feature. Combining these findings can reveal the RAs and remove them.
To test the method, a handheld probe with integrated diode lasers was used for PAI. These diode lasers emit light at 4 wavelengths (808, 915, 940 and 980 nm). We performed experiments in phantoms and in vivo. Results show that this is a promising method for correcting RAs, potentially in real-time.

Photoacoustic imaging
PAI is an imaging technique using pulsed laser irradiation to generate US waves which are subsequences of pressure changes due to thermal expansion and relaxation. The generated initial pressure is described as [23][24][25]: where a μ is the absorption coefficient [cm −1 ], Φ is the light fluence [J/cm 2 ], and Γ (Grüneisen parameter) is a dimensionless parameter and is defined as Light propagating in the tissue is scattered and absorbed. Since scattering and absorption are strongly dependent on the wavelength, light at different wavelengths reaches different depths [26,27]. Therefore, the light fluence inside the tissue depends on both the excitation wavelength and the position.
The absorption coefficient, a μ , is a wavelength-dependent optical property of the absorber. The generated initial pressure, p , can be rewritten as a function of the excitation wavelength and the local position: Figure 1 illustrates the principle of RAs in PAI. A part of generated US waves (blue) is reflected at the acoustic reflector, seen in Fig. 1(a) The reflected US waves (red) propagating back to the detector resemble a virtual acoustic source, so-called RA, located at a larger depth. Figure 1(b) is a reconstructed PA image of a phantom representing this situation. The phantom was made of a black thread placed above a plastic petri dish lid, and demi-water was used as an acoustic coupling medium. An RA at a larger depth is clearly visible in this image.

Method
The principle image represe excitation wav images with m interest ( In the ide high correlatio with correlati appear at a lar longer propag be deeper than Features in th remove them. RAs are re corrected ima However, this features, the which were re in Fig. 4 Fig. 3(a) Offline pr running Matla

Experime
To demonstra as in vivo. reconstruction In each ex sent repeatedl then acquired rate of 1 kHz RAs.

Phantom
A phantom w dish lid (Grei m μ , seen in Fig. 6(b). The phantom was 3.5% Intralipi as an acoustic scattering coe nm based on   Figure 7(a also Fig. 15 periment, th ρ = es were subseq 8(a). This corre n Fig. 8(b). Fe does not consi f the acquired of a phantom exper ure being assigne ponses of the featu ade of the s Fig. 7(b), the ature 57 for on assumption de e 57 and its RA of the absorb cond assumptio ral responses es.

In vivo
We also asses where bones underneath th Figure 9(a Figure 9(d) is are represente 9(c), adapted the periosteum  Fig. 9(a) from [33], rev m and bone are 9. RAs in in vivo ed from "Sobotta: of an in vivo imag r was imaged same as describ tures in an acqui The final correcte and the final corre od with in vivo derneath the sk emi-water was vivo PA imag cting the experi and Fig. 9 The spectr they represen representing t of different sO It is worth features 18 an remove featur Figure 11 removed.    method without all other pixels Figure 13(b) this pixel is co oefficient map.
identified RA d pixel.  In the corrected images without segmentation, the absorber in the bottom left corner of the phantom was observably shrunk. This could be due to the small amount of data points for the Pearson correlation coefficient.
The "tail" of the absorber in the top right corner of the phantom is likely a reconstruction artifact. A part of it was removed in the corrected image without segmentation. The reason is that the removed pixels in the "tail" were identified as RAs of the absorber.
The threshold for the method without segmentation was applied the same as the threshold for the method with segmentation (0.95). However, min z Δ was set as 2 mm for compensating the size of features.

Discussion
In in vivo imaging, the RAs of the blood vessels look like excess vessels and could be erroneously recognized as angiogenesis or hyper-vascularization, hallmarks of various diseases such as cancer or rheumatoid arthritis. With not much experience of RAs, misdiagnosis might happen. Simple solutions such as thresholding or limiting the imaging depth may be able to remove RAs in cases that the RAs are not accompanied by real image features with the same amplitude or depth range. If such image features exist, however, these simple solutions are not appropriate.
Compared to previously reported methods for reducing RAs, the proposed method offers significant advantages. First of all, the method works automatically and performing it does not require experience or training of the users, as is the case with DCA or PAFUSion in which the users have to hold the probe perpendicularly to the acoustic reflectors. Secondly, no US image is needed. Acquiring US images with multiple plane wave angles is more time consuming and comes at a higher processing expense. Thirdly, as the method does not need US images to detect RAs, matching features between PA and US images (as in PAFUSion) is not necessary resulting in detected RAs being completely removed in this proposed method. Fourthly, unlike deep learning, the method does not require training with various generated PA distribution sizes and geometries, and acoustic characteristic of the sample which might be unknown in in vivo imaging. Finally, the proposed method enhances the advantages of multispectral PAI.
Out-of-plane artifacts (direct and indirect out-plane-artifacts) can appear in the acquired PA images, especially artifacts caused by the skin. PAFUSion does not work for indirect outof-plane artifacts. In contrast, they can be treated under the proposed method if the direct outof-plane artifacts also appear in the image. However, both methods cannot handle direct outof-plane artifacts. Another method for these artifacts is needed. Our future work will focus on a complete method combining this proposed method and a new method for out-of-plane artifacts.
The proposed method exploits the variance of light distribution of different wavelengths. In other words, the local illumination spectrum needs to be different at different depths for the method to work. In a scenario that the absorption and scattering spectrum is flat, resulting in a similar illumination spectrum at different depths, this method will not work. However, this is not likely in clinical imaging where the absorption coefficient varies with wavelength and various tissues [26]. min z Δ , which defines a region below an original PA feature where the method will not correct its RAs, is the main limitation of this method. In our experiments demonstrating the method, min z Δ was 1.7 mm corresponding to the threshold of the correlation coefficient of 0.95. However, this value is not a critical limitation in clinically relevant scenarios. In in vivo imaging, with different layers of tissue such as skin and muscle, optical heterogeneity would be stronger resulting in a smaller min z Δ than in this study. Additionally, min z Δ can be further reduced by selecting wavelengths at which light fluence varies strongly along the depth.
The propo wavelengths. excitation wa principle of th performed av pulse variatio method. How In this m generated init itself and thus in our work. N The metho another absor the spectral r Another appro A small n relation betw investigated in Our metho calculation an pixels (~25 m seconds runni time. For the sufficient (~1 consumption mm, the size o GPU can be u

Conclusio
The proposed exploiting lo phantoms and 940 and 980 with no separ compact PAI method for m