Abstract: 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.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
OCIS codes: (170.5120) Photoacoustic imaging; (170.3880) Medical and biological imaging. References and links
1. 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).
2. M. Heijblom, D. Piras, F. M. van den Engh, M. van der Schaaf, J. M. Klaase, W. Steenbergen, and S. Manohar, “The state of the art in breast imaging using the Twente Photoacoustic Mammoscope: results from 31 measurements on malignancies,” Eur. Radiol. 26(11), 3874–3887 (2016).
3. M. Toi, Y. Asao, Y. Matsumoto, H. Sekiguchi, A. Yoshikawa, M. Takada, M. Kataoka, T. Endo, N. Kawaguchi-Sakita, M. Kawashima, E. Fakhrejahani, S. Kanao, I. Yamaga, Y. Nakayama, M. Tokiwa, M. Torii, T. Yagi, T. Sakurai, K. Togashi, and T. Shiina, “Visualization of tumor-related blood vessels in human breast by
photoacoustic imaging system with a hemispherical detector array,” Sci. Rep. 7(1), 41970 (2017). 4. P. J. van den Berg, K. Daoudi, H. J. Bernelot Moens, and W. Steenbergen, “Feasibility of
photoacoustic/ultrasound imaging of synovitis in finger joints using a point-of-care system,” Photoacoustics 8, 8–14 (2017).
5. J. Jo, G. Xu, M. Cao, A. Marquardt, S. Francis, G. Gandikota, and X. Wang, “A functional study of human inflammatory arthritis using photoacoustic imaging,” Sci. Rep. 7(1), 15026 (2017).
6. K. Jansen, M. Wu, A. F. van der Steen, and G. van Soest, “Photoacoustic imaging of human coronary atherosclerosis in two spectral bands,” Photoacoustics 2(1), 12–20 (2014).
7. K. Daoudi, P. J. van den Berg, O. Rabot, A. Kohl, S. Tisserand, P. Brands, and W. Steenbergen, “Handheld probe integrating laser diode and ultrasound transducer array for ultrasound/photoacoustic dual modality imaging,” Opt. Express 22(21), 26365–26374 (2014).
8. J. J. Niederhauser, M. Jaeger, R. Lemor, P. Weber, and M. Frenz, “Combined ultrasound and optoacoustic system for real-time high-contrast vascular imaging in vivo,” IEEE Trans. Med. Imaging 24(4), 436–440 (2005). 9. C. Kim, T. N. Erpelding, L. Jankovic, M. D. Pashley, and L. V. Wang, “Deeply penetrating in vivo
photoacoustic imaging using a clinical ultrasound array system,” Biomed. Opt. Express 1(1), 278–284 (2010). 10. C. Haisch, K. Eilert-Zell, M. M. Vogel, P. Menzenbach, and R. Niessner, “Combined optoacoustic/ultrasound system for tomographic absorption measurements: possibilities and limitations,” Anal. Bioanal. Chem. 397(4), 1503–1510 (2010).
11. K. Sivasubramanian, V. Periyasamy, R. A. Dienzo, and M. Pramanik, “Hand-held, clinical dual mode ultrasound - photoacoustic imaging of rat urinary bladder and its applications,” J. Biophotonics 11(5), e201700317 (2018).
#335295 https://doi.org/10.1364/BOE.9.004613
12. 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. Biophotonics 11(1), e201700061 (2018).
13. M. K. A. Singh, V. Parameshwarappa, E. Hendriksen, W. Steenbergen, and S. Manohar, “Photoacoustic-guided focused ultrasound for accurate visualization of brachytherapy seeds with the photoacoustic needle,” J. Biomed. Opt. 21(12), 120501 (2016).
14. M. Jaeger, L. Siegenthaler, M. Kitz, and M. Frenz, “Reduction of background in optoacoustic image sequences obtained under tissue deformation,” J. Biomed. Opt. 14(5), 054011 (2009).
15. M. Jaeger, J. C. Bamber, and M. Frenz, “Clutter elimination for deep clinical optoacoustic imaging using localised vibration tagging (LOVIT),” Photoacoustics 1(2), 19–29 (2013).
16. H.-M. Schwab, M. F. Beckmann, and G. Schmitz, “Photoacoustic clutter reduction by inversion of a linear scatter model using plane wave ultrasound measurements,” Biomed. Opt. Express 7(4), 1468–1478 (2016). 17. M. K. A. Singh and W. Steenbergen, “Photoacoustic-guided focused ultrasound (PAFUSion) for identifying
reflection artifacts in photoacoustic imaging,” Photoacoustics 3(4), 123–131 (2015).
18. D. Allman, A. Reiter, and M. A. L. Bell, “Photoacoustic Source Detection and Reflection Artifact Removal Enabled by Deep Learning,” IEEE Trans. Med. Imaging 37(6), 1464–1477 (2018).
19. T. Petrosyan, M. Theodorou, J. Bamber, M. Frenz, and M. Jaeger, “Fast scanning wide-field clutter elimination in epi-optoacoustic imaging using comb-LOVIT,” in Ultrasonics Symposium (IUS), 2017 IEEE International (IEEE, 2017), p. 1.
20. H.-M. Schwab and G. Schmitz, “An advanced interpolation approach for photoacoustic clutter reduction based on a linear plane wave scatter model,” in Ultrasonics Symposium (IUS), 2016 IEEE International (IEEE, 2016), pp. 1–4.
21. M. K. A. Singh, M. Jaeger, M. Frenz, and W. Steenbergen, “Photoacoustic reflection artifact reduction using photoacoustic-guided focused ultrasound: comparison between plane-wave and element-by-element synthetic backpropagation approach,” Biomed. Opt. Express 8(4), 2245–2260 (2017).
22. M. K. A. Singh, M. Jaeger, M. Frenz, and W. Steenbergen, “In vivo demonstration of reflection artifact reduction in photoacoustic imaging using synthetic aperture photoacoustic-guided focused ultrasound (PAFUSion),” Biomed. Opt. Express 7(8), 2955–2972 (2016).
23. V. E. Gusev, and A. A. Karabutov, “Laser optoacoustics,” NASA STI/Recon Technical Report A, 93 (1991). 24. A. Oraevsky and A. Karabutov, “Optoacoustic tomography,” Biomedical photonics handbook 34, 1– 34 (2003). 25. A. A. Oraevsky, S. L. Jacques, and F. K. Tittel, “Measurement of tissue optical properties by time-resolved
detection of laser-induced transient stress,” Appl. Opt. 36(1), 402–415 (1997).
26. A. Bashkatov, E. Genina, V. Kochubey, and V. Tuchin, “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D Appl. Phys. 38(15), 2543–2555 (2005).
27. V. V. Tuchin, and V. Tuchin, Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis (SPIE Press Bellingham, 2007).
28. M.-L. Li, J.-T. Oh, X. Xie, G. Ku, W. Wang, C. Li, G. Lungu, G. Stoica, and L. V. Wang, “Simultaneous molecular and hypoxia imaging of brain tumors in vivo using spectroscopic photoacoustic tomography,” Proc. IEEE 96(3), 481–489 (2008).
29. I. Sobel, and G. Feldman, “A 3x3 isotropic gradient operator for image processing, presented at a talk at the Stanford Artificial Project,” in Pattern Classification and Scene Analysis, R. Duda and P. Hart, eds. (John Wiley & Sons, 1968), pp. 271–272.
30. J. Benesty, J. Chen, Y. Huang, and I. Cohen, “Pearson correlation coefficient,” in Noise Reduction in Speech
Processing (Springer, 2009), pp. 1–4.
31. M. Jaeger, S. Schüpbach, A. Gertsch, M. Kitz, and M. Frenz, “Fourier reconstruction in optoacoustic imaging using truncated regularized inverse k-space interpolation,” Inverse Probl. 23(6), S51–S63 (2007).
32. R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express 16(8), 5907–5925 (2008).
33. R. Putz and R. Pabst, Sobotta: Atlas of Human Anatomy (Elsevier, Urban & Fischer, 2006).
34. D. G. Bonett and T. A. Wright, “Sample size requirements for estimating Pearson, Kendall and Spearman correlations,” Psychometrika 65(1), 23–28 (2000).
1. 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
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–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.
2. Theory
2.1. 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–25]:
a
p= Γ Φμ (1)
whereμais the absorption coefficient [cm−1],Φis the light fluence [J/cm2], andΓ(Grüneisen parameter) is a dimensionless parameter and is defined as 2/
P
c C
β
Γ = , whereβ is the thermal expansion coefficient [K−1], c is the speed of sound [m/s], and CP is the isobaric specific heat [J/kgK].
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:
( , , , ).
p= f λ x y z (2)
2.2. Reflection artifacts in photoacoustic imaging
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.
Fig. 1 image water
In a clinic tumor can ref tumor. This c hemoglobin [3 3. Method The principle image represe excitation wav images with m interest (ROI) values reveal following two 1. Absorb spect 2. Both d If the abov spectral respo misidentified at the end of t Figure 2 multiple wav selected for segmentation information is image which recover the sh 1. RA in PAI. (a) e of a phantom (a represents this situ
cal scenario w flect US signal can negatively
3, 28]. e of our metho
ents the genera velength. Expl multiple wavele
), they show th s the spectral o assumptions: bers with ident tral responses d direct and reflec ve two assump onse of their as one reflecti this section.
shows the flow velengths. Of
segmentation step are appli s then used in t is a segmented hape of the rem
) A deep reflector black thread plac uation.
where there are ls generated fr
affect the abil
od is based on ated initial pres
loiting this prin engths of light he same structu responses of a tical optical pr due to differen cted PA signal ptions are fulfil
source (real a ion artifact of wchart of the the acquired which detects ied to all the o the RA correct d image is fur maining feature
r leads to reflectin ced above a plastic
e a few blood rom the blood lity to assess tu
Eq. (2) where ssure as a func nciple with mu t is obtained. A ure of the samp absorbers in th roperties locate nt local light flu s convey the o lled, the spectr absorber) and the other. The method. Imag images, the im image featur other images to
tion step to ide rther processed s giving the fin
ng US waves. (b) c petri dish lid) em
vessels locate vessels causin umors based o
e the pixel va ction of the loc ulti-wavelength As all images ar ple. Studying t he images. Ou ed at different uence. optical properti ral response of two identical e first assumpt
ges of the sam mage giving res. The featur
o obtain their entify and remo d through the d nal corrected im ) An acquired PA mbedded in demi-ed above a tum ng RAs surroun on oxygen satu
alue in the acq cal light fluenc h PAI, a sequen re of the same the changes of ur method relie positions give es of the sourc f RAs is identi absorbers wi tion is discusse mple are acqu
the strongest ures extracted spectral respo ove RAs. The de-segmentatio mage. A -mor. The nding the uration of quired PA ce and the nce of PA region of f the pixel es on the different ce. ical to the ill not be ed further ired with signal is from the onse. This corrected on step to
In the seg Sobel edge de Sobel edge th Figure 3 i represents tw arrows). Prop lead to over-th cannot be de where bottom case. Several We avoid consequence, process. In t surrounding t giving separa absorber is se one absorber, gmentation ste etection algorit hreshold is supp illustrates this wo blood vesse perly segmentin hresholding or tected resultin m features are n features are de over-threshold under-thresho his process, w the peaks to z ate features, se egmented into
they share the
Fig. 2. The f
ep, an automat thm [29] is im ported by Matl segmentation els (upward blu
ng these featur r under-thresho ng in feature l not detectable. etected as one s ding by choosi olding might h we find all pe ero. Under-thr en in Fig. 3(d) a number of f e spectral respo flowchart of the m tic segmentatio mplemented to d ab. method. Figur ue arrows) an res is expected olding. In the c loss. Figure 3( In contrast, F single feature i ing half of the happen. We fur eaks in the im resholding is s
). It might lea features. Howe onse of the abso
method. on algorithm detect image f ure 3(a) is a sa nd their reflect d. However, ap case of over-th (b) shows an Fig. 3(c) shows in this case. threshold calc rther process t mage and then significantly im ad to over-segm ever, since the orber. which is base features. Comp ample PA ima tion (downwar pplying a thres hresholding, we over-threshold s an under-thre culated by Mat
the image with n set a part mproved after mentation in w ese features are
ed on the puting the ge which rd yellow shold can eak edges ding case esholding tlab. As a h a peak-of pixels this step which one e parts of
Fig. 3 vessel thresh thresh To obtain applied to all images, givin compared to e where A and B, σA and σB 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(b). 3. An example of ls (upward blue a holding segmented holding and peak-p
features spect l other images ng that feature
each other usin
B are spectral B are the stand entifying and r on coefficients ions exceeding rger depth and gation and atten
n their real abs he each group emoved by set age of Fig. 3 s corrected ima de-segmentati emoved in the the segmentation arrows) and their d image. (c) An u processed segment
tral response, th s. Of each fea e’s spectral res ng the Pearson ( ,A ρ l responses of dard deviation o removing RAs s from all corre g ρth are grou d have a weake nuation. An ex sorber at least are analyzed b tting the pixel v (a). Features age is a segme on step is app peak-process a
process. (a) The reflection (downw under-thresholding
ted image.
he detected fea ature the max sponse. Spectr correlation coe cov( , ) ) A B A B B σ σ = two features, of A and B resp step, a thresh elation coeffici uped together er signal than th xtra condition i min z Δ , which i based on these values of the R 8, 9, and 11 ented image. T plied. All pixe are recovered g
original image sh ward yellow arrow g segmented imag
atures from the ximum pixel v ral responses o efficient [30]: ) cov( , )A B is th pectively. hold, ρth, is th ients. Features as suspicious he correspondi s used to ident is described at e conditions to RA features to detected in F To recover the els surroundin giving the fina
howing two blood ws). (b) An over-ge. (d) The
under-e sunder-egmunder-entation value is taken of all features he covariance hen applied to with spectral r RAs. In addit ing real absorb tify RAs, that R t the end of thi o identify RAs zero. Figure 4 Fig. 3(d) are shape of the r ng the remaini al corrected im d -n step are from all are then (3) of A and o separate responses tion, RAs ber due to RAs must s section. and thus 4(a) is the removed. remaining ing peaks mage, seen
Fig. 4 The fi Segmentat might be no implemented an object to s than features, the same as th A compar experimental In a highly thus two iden correlation co avoid this, a extra conditio absorbers wit In other word features are a would be con measurements discussed furt 4. Setup Experiments w probe is conn Netherlands) maximum sam research mode The US tr bandwidth of central 64 ele different wave 4. An example of c inal corrected imag
tion benefits ot obtained.
for analyzing study the spect are correlated he analysis use rison of the m results section y scattering me ntical absorbers oefficient excee minimum dist on to assess tha
h the same spe ds, Δzmin defi
assessed as its R nsidered as an
s were perform ther in the disc
were carried o nected to a com
for the acquis mpling frequen e so that raw d ransducer array f 100%. It com ments were us elengths (808, correcting RAs in ge. image feature Therefore, an images. In this tral responses. d to each other. d for segmente method with an n, and further di edium, a local s in that region eding ρth. As tance Δzmin in at whether one ectral response ines a region b RAs. In additio RA. This can med and will b ussion section. out using a han mmercial ultras sition of US a ncy of 50 MH data could be ac y in the handhe mprises 128 el sed. Diode lase 915, 940, and PA images. (a) T analysis, how nother approa s approach, eac Particularly, s The analysis b ed features. nd without se iscussed in the region can hav n might have t
a consequence n the depth bet e feature is an e. The value of below a PA i on, if one abso n be avoided w be reported in
.
ndheld PAI pro sound scanner and PA images Hz with 12 bit
cquired in an e eld probe has a ements with a ers integrated in
980 nm) at a r
The corrected imag
wever, a prop ach without ch pixel of the spectral respon based on the c egmentation w e discussion sec ve light fluenc the same spect e, assumption tween the two RA of the oth f Δzmin is relat image feature orber is over-se with Δzmin. To the experimen obe, depicted i MyLabOne (E s. The scanner digitization. T external PC for a center freque a pitch of 0.3
nto the handhe repetition rate o
ge of Fig. 3(a). (b)
perly segmente segmentation e image is cons nses of all pixe correlation coef will be presente ction. ce nearly homo tral responses, 1 is not approp o features is us her one or two ted to the valu where no oth egmented, one o determine th ntal results sec
in Fig. 5. The Esaote Europe r can acquire This device wa r offline proces ency of 7.5 MH mm. For our eld probe emit of up to 10 kHz ) ed image is also sidered as els, rather fficient is ed in the ogeneous, , giving a priate. To sed as an o separate ue of ρth. her image e segment his Δzmin, ction, and handheld BV, The data at a as used in ssing. Hz with a study the light at 4 z.
Table 1 pr In addition, th and 980 nm, r output. These wavelengths. Offline pr running Matla 5. Experime To demonstra as in vivo. reconstruction In each ex sent repeatedl then acquired rate of 1 kHz RAs. 5.1. 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 Fig. 5 resents specific he angles at th respectively, du e differences b Table 1 wavelength [nm 808 915 940 980 rocessing of da ab R2016b. ental results
ate the feasibili The PA imag n algorithm [31 xperiment, 4 la ly for 100 time d by averaging . The US imag ms was made of tw iner Bio-One G Fig. 6(a). A s e lid was posi mounted on a id 20% (Frese c coupling me efficient of the
[32]. Figure
. Photo and schem
cations of the d he output are 5
ue to diode las between wavel 1. Lasers specific
m] pulse energyat the output 0.96 0.98 0.89 0.82
ata was done
ity of the meth ge reconstruct 1]. ser pulses of 4 es. 4 PA image signal over 10 ge was used to wo black threa GmbH, Germa schematic draw tioned underne a mount (CP02 enius Kabi, Th dium as well a solution was 6(c) shows a matic drawing of th diode lasers w 53.1, 55.6, 47.8 ers placed at d lengths add to ations at a repeti y t [mJ] pulse wid 84.2 88 98.9 94.2 on a PC (Intel hod, we perform tion was done
different wave es at 4 differen 00 pulses. The o verify the loc
ads, with the d any) as an acou wing of a cross
eath one black /M, Thorlabs, e Netherlands) as an optically estimated as
μ
combined PA he handheld probe. working at the r 8, and 50.3 de different stackso the light flue ition rate of 1 kHz dth [ns] Fluenceoutput 1.04 1.01 0.95 0.87 l Core i7 3.41 med experimen e using a Fo elengths follow nt wavelengths e diode lasers w cation of absor diameter of 200 ustic reflector, s-section of th k thread as an Germany) to f ) in demi-wate y scattering ba ' s
μ
= 6 cm−1 at A and US ima . repetition rate o egrees at 808, 9 and a prism at ence variation z. e at the [mJ/cm2] 1 GHz, 8 GB o nts in phantom ourier transfor wed by 1 US pu s and 1 US im were run at a r rbers and corre0-250 mμ , an , with thicknes he phantom is
acoustic refle fixate it in a so er. This solutio ackground. The t the waveleng age illustrating of 1 kHz. 915, 940, t the light n of these of RAM) ms as well rm based ulse were mage were repetition esponding nd a petri ss of 750 shown in ector. The olution of on served e reduced gth of 900 g a
cross-section of this which reflect visualized at features. As t thread above which perhap 4 PA ima intensity of t reflectivity of the US image signal resultin acquired at a for segmentat Fig. 6 Comb 980 n Figure 7(a also Fig. 15, A spectral respo s phantom. Th US waves. T expected pos there was no ab
the lid. In the s is a reconstru ages of the RO the reflections f the petri dish e’s case, the U ng in high in 940 nm wave tion. 6. (a) A phantom u bined PA and US i m). a) shows the s Appendix 1). T onse was norma
e gray color pa The hot color p sitions relative
bsorbers under e PA image, t uction artifact. OI correspondi s was not stro h lid in that cou US transducer ntensity of the elength had the
used for experimen image. (d) 4 PA im
egmented ima The spectral re alized with the
art is the US im part is the PA e to the lid. U
rneath the lid, here is a long ing to 4 wave ong. The expl upling medium r array generat e reflector in e strongest sign nts. (b) A schema mages acquired at ge acquired at sponse of seve e maximum val mage showing A image where Underneath the these features g “tail” of the elengths are sh lanation migh m was not high ted higher pre the acquired U nal therefore t atic cross-section o t 4 wavelengths (8 t 940nm with n eral features is lue). g two surfaces e two black th e lid were so were RAs of absorber abov hown in Fig. 6 ht be that the h. On the other essure compare US image. Th this image was
of the phantom (c) 808, 915, 940, and numbered feat shown in Fig. of the lid hreads are me more the black ve the lid 6(d). The acoustic r hand, in ed to PA he image s selected ) d tures (see 7(b) (the
Fig. 7 and a norma Two blac properties. H absorbers is o This matches spectral respo words, optica reflection con Of the o “similarity” o The correlatio to 1) and thes that they are R 26 are 0.937 sufficient to s likely an RA as the intensit response and Feature 24 26 51 52 54 57 58 64 In this exp These feature image, Fig. 8 image, seen in noise level, it comparison o 7. Image analysis o ll pixels in a featu alized spectral resp
ck threads ma owever, from observable (fea with the first onses of feature al properties o nfirming the sec btained spectr of the response on coefficients se features app RAs of feature and 0.848 res separate feature of feature 57, ty of feature 5
thus the correl Table 2. Correlat 24 26 1 0. 0.979 1 0.843 0. 0.974 0. 0.898 0. 0.937 0. 0.901 0. 0.928 0. 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. Table 2 sho between featu pear at larger d e 57. The corre spectively. A t es 24 and 26 fr however has a 1 is close to th ation. tion coefficients o experiment 6 51 .979 0.843 0.932 .932 1 .912 0.702 .791 0.532 .848 0.618 .798 0.557 .846 0.643 = 0.95 was app quently remov ected image w ature 51 was n iderably affect PA image and
riment. (a) Segme d the maximum v ures. ame materials e difference in e thread and fe escribed in the As, features 54 er were conse on. mutual corre ows the correl ure 57 and featu depths than fea elation coefficie threshold ρth rom features 52 a low correlati he background, of obtained spectr t (see also Data F
52 54 0.974 0.8 0.912 0.7 0.702 0.5 1 0.9 0.963 1 0.978 0.9 0.952 0.9 0.951 0.9 plied identifyin ved from the i was then de-seg not removed. H t the interpretat the final corre
ented image with n value of that featu
s have identi n the spectral features 24 and e method secti 4 and 58, are hi erved in the elations are ca lation coefficie ures 52, 54, 58 ature 57 with a ents of feature in between 0 2, 54, 58 and 6 ion coefficient , that would hi ral responses in th File 1). 57 898 0.937 791 0.848 532 0.618 963 0.978 0.994 994 1 995 0.995 977 0.99 ng features 52, image giving gmented to ob However, as thi
tion of the ima ected image. numbered features ure. (b) Maximum ical optical ab response of t d 26 for anothe ion. On the oth
ighly identical US waves in alculated to d ents of these r 8, and 64 are hi a lower signal e 57 and feature .94 and 0.97 w 64. Feature 51, . This can be e ighly affect the he phantom 58 64 0.901 0.92 0.798 0.84 0.557 0.64 0.952 0.95 0.995 0.97 0.995 0.99 1 0.99 0.992 1 , 54, 58 and 64 a corrected se btain the final
is feature is clo age. Figure 8(c) , m bsorption these two er thread). her hand, l. In other spite of determine esponses. igh (close revealing es 24 and would be , which is explained e spectral 28 46 43 51 77 9 92 4 as RAs. egmented corrected ose to the ) shows a
Fig. 8 segme the ac 5.2. 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 adapte photo The finger analyzed the s 8. Processing fea ented image. (b) T cquired PA image a
ssed the metho are close und he probe and de a) shows an in s a photo depic ed in 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(b) ealing the peri e therefore disc
PAI. Acquired P Atlas of human a ging experiment.
with 4 wavele bed in the phan
ired PA image o d image after de-ected image.
o experiments. kin giving clea used as the US ge of a cross-se imental config respectively. T osteum and bo closed as RAs o PA (a) and US (b) natomy” [33] show
engths and the ntoms section. of the phantom. ( -segmentation. (c) . Our experime ar RAs. Finger S coupling med ection of a hea guration. Acqui The US image one. In the PA
of the skin and
) images of a fing ws a cross-section
acquired imag . The processed
(a) The corrected A comparison of ents focused o rs were placed dium. althy volunteer ired PA and U e is compared image, feature d blood vessels ger. (c) An image n of a finger. (d) A
ges were proce d image with n d f on fingers d ~7 mm r’s finger. US images with Fig. s beneath . e A essed and numbered
Fig. 1 featur The spectr they represen representing t of different sO It is worth features 18 an remove featur Figure 11 removed. Fig. 1 finger th ρ = 0.9 coefficient tab 10. Image analysi res. (b) Spectral res
ral responses o t different chro two blood vess O2 in these bloo h noting that f nd 124 which re 50 will also 1 shows acqui 11. Correcting RA r. (b) The corrected 95 was used i ble). is of an in vivo e sponses of the feat
of features 32 omophores (m sels also have od vessels, or o feature 50 whi are two blood remove feature ired and corre
As in an in vivo im d image. in this experim experiment. (a) A tures. and 56 are ob melanin and blo different spect of different loc ich is an RA d vessels. Incr es 18 and 124 l ected images maging experimen
ment (see Dat
A segmented imag bservably differ ood). Interestin tral responses. cal spectra of th of the skin ha reasing the seg leading to losin where all RA
nt. (a) An acquire
ta File 2 for
ge with numbered
rent from each ngly, features 5 . This might b he excitation li as higher inten gmentation thr ng real feature As were identi ed PA image of a a complete co d h other as 56 and 77 e a result ight. nsity than reshold to es. ified and a orrelation
5.3. Minimum As mentioned two identical same spectral assumption 1 miscorrect on To evalua comparing th optical proper cm−1 at the w Netherlands), Milli-Q wate 0.456, 0.534, measured usi suture wires ( two identical fixated and t Thorlabs, Ger At each Δ processed to c measurement minimum valu with the distan
1 z Δ = 0.68 m 2 z Δ in these m Fig. 1 vertic = 0.68 For a thre bar in Fig. 1 spectral respo coefficient, 3 points would m vertical dist d in the metho absorbers with l responses, g formulated in ne absorber as a ate the correla e spectral resp rties. The med
avelength of 9 and 4x10−4 vo r. The absorp and 0.776 cm ng a photo-sp (USP 3/0, diam absorbers wer the other one rmany) to adju
z
Δ , PA image calculate the c was repeated ues are shown nce. Figure 12 mm and a measu measurements w 2. (a) Correlation al distance. (b) Sp 8 mm and Δz2 = eshold of the co 2(a) is large. onses was cal data points can affect the conf
tance Δzmin od section, in a h a vertical dis giving a corre n section Meth an RA of the ot ation as a fun ponses of two ium was a solu 900 nm [32], 10 olume fraction tion coefficien m−1 at the wave pectrometer (U meter of 0.24 m re used. These was attached st the vertical d es were acquir correlation coe d 4 times at in Fig. 12(a) p (b) shows spec urement at Δz2 were 0.9916 an coefficient of spe pectral responses o 2.1 mm. orrelation coef The reason fo lculated with n give a meani fidence interva a region havin stance (along d lation exceedi hod. The meth
ther one. nction of zΔ , identical abso ution of 2% In 0−4 volume fra of black ecolin nt of this solu elength of 808 UV-2600, Shim mm, Vetsuture two wires wer d to a motori
distance zΔ be red using the
fficient of spe each distance presenting the b ctral responses 2 = 2.1 mm. T nd 0.9126 respe ctral responses of of the two suture w
fficient of 0.95 or this might b only 4 wavel ingful coefficie l [34]. ng nearly homo depth) less than ing ρth, whic hod, therefore, , we performe orbers in a me ntralipid 20% w action of India ne (Royal Tale ution (without , 915, 940, an madzu, The Ne e Nylon, The N re embedded i ized translatio etween the two
four waveleng ctral responses e zΔ . The av behavior of the s of the two wir The correlation
ectively.
f two identical abso wires at two differe
5, Δzmin is 1.7 be that the co lengths. For th ent. However, ogeneous light n Δzmin might ch results in f might miside ed several exp edium mimicki with estimated ink (Royal Ta ens, The Nether
Intralipid) wa nd 980 nm resp etherlands). Tw Netherlands) m in the solution, n stage (MTS o wires. gths. Images w s of the two w verage, maxim e correlation co res of a measu coefficients at
orbers versus their ent distances Δz1 mm. The vert orrelation coeff the Pearson co a small numbe t fluence, t have the failure of entify and periments ing tissue d ' s
μ
= 3.5 alens, The rlands) in as 0.599, pectively, wo black mimicking , one was S50A-Z8, were then wires. The mum and oefficient rement at t Δz1 and r tical error fficient of orrelation er of datathe considered Fig. 1 vivo P mm b pixels Figure 14 image, the co of the phantom Fig. 1 PA im segme correc respec d pixel. 3. RA identificati PA image. (b) The elow the considere s) of the considered
shows the resu rrected image m. The bottom
14. Comparison of mage, the correc entation of the p cted image with s ctively.
ion of the method e correlation coeff
ed pixel (values ab d pixel.
ults using the t with segmenta m images are im
f the method with cted image with
hantom, respectiv segmentation, and
without segmenta ficient map of a pix
bove 0.95 are colo
two approache ation, and the c mages in vivo in
h and without segm segmentation, a vely. (d), (e), and d the corrected im
ation in an in vivo xel in the skin wit ored red). (c) Ident
es. The top ima corrected imag n the same ord
mentation. (a), (b) and the corrected d (f) An acquire mage without segm
o image. (a) An in th others at least 2 tified RAs (yellow
ages are an acq ge without segm der. ), (c) An acquired d image without ed PA image, the mentation in vivo n 2 w quired PA mentation d t e ,
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, Δzmin was set as 2 mm for compensating
the size of features. 6. 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 of-plane artifacts. In contrast, they can be treated under the proposed method if the direct of-plane artifacts also appear in the image. However, both methods cannot handle direct out-of-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, Δzminwas 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 Δzmin than in this study. Additionally, Δzmin can be further
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 7. Conclusio The proposed exploiting lo phantoms and 940 and 980 with no separ compact PAI method for m Appendix 1 Nevertheless, a od might not b rber or another response of thi oach such as P number of data ween the num n the future wo od is a post-pr nd thus reduce mm depth), corr ing in Matlab, e method with 12 GB for th significantly. S of images can utilized for furt
on
d method can cal light flue d in vivo were nm. Results s rate ultrasound system suitabl medical use. : Phantom’s Fig. 15. Se an appropriate r be able to iden r RA. In that ca is feature will AFUSion is ne a points might mber of wavel ork. rocessing meth es the frame ra recting with an respectively, s hout segmenta his case) vecto
Since the imag be reduced res ther acceleratio identify and r ence variation
performed for how the poten
images neede le for clinical u segmented gmented image of reconstruction ntify an RA if ase, they appea
not highly co eeded to identif affect the Pea engths and th hod. Therefore, ate. However, nd without segm showing the po ation, if the m orizing data c ging depth cap sulting in less e on of the metho remove in pla between mu r a proof of con ntial of the me d. In addition, use, demonstra image
f the phantom with
algorithm sho f it appears at ar as one singl orrelate with th fy the RAs. arson correlatio he efficacy of , compared to p for images wi mentation take otential of the m memory of the can be accom pability of this expensive calc od, especially w ane reflection ultiple wavelen ncept, using 4 ethod for corre experiments w ating the practi
h all features numb
ould be conside the same posi le feature and t the RA’s real on coefficient f our method pure PAI, it co ith the size of e ~0.2 seconds method runnin e computing d mplished reduc handheld prob culation. Never without segmen artifacts (RAs ngths. Experim wavelengths: ecting RAs in were carried ou ical applicabili bered. ered. ition with therefore, absorber. [34]. The d will be osts more f 512×64 and ~2-3 ng in real-device is cing time be is ~10 rtheless, a ntation. s) in PAI ments in 808, 915, real-time ut using a ity of this
Appendix 2 Funding European Un 731771) Acknowledg The authors w Twente, for th University of Disclosures The authors d : Segmented Fig. 16 nion’s Horizon gments would like to heir discussion Twente, for th s
declare that ther
d in vivo imag
6. Segmented in viv
n 2020 resear
acknowledge ns. The authors heir help with p re are no confl
ge
vo image with all
ch and innov
David Thomp are also gratef phantom prepar licts of interest
features numbered
ation program
pson and Bart ful to Tom Kno
ration. t related to this d. mme (CVENT, Fischer, Univ op and Wilma article. , H2020-versity of Petersen,
