Reducing artifacts in photoacoustic imaging by using multi-wavelength excitation and transducer displacement

Reducing artifacts in photoacoustic imaging

by using multi-wavelength excitation and

transducer displacement

H

N

Y. N

*

W

S

Biomedical Photonic Imaging Group, Faculty of Science and Technology, University of Twente, P.O. Box 217, 7500 AE, Enschede, The Netherlands

*h.n.nguyen@utwente.nl

Abstract: The occurrence of artifacts is a major challenge in photoacoustic imaging. The

artifacts negatively affect the quality and reliability of the images. An approach using multi-wavelength excitation has previously been reported for in-plane artifact identification. Yet, out-of-plane artifacts cannot be tackled with this method. Here we propose a new method using ultrasound transducer array displacement. By displacing the ultrasound transducer array axially, we can de-correlate out-of-plane artifacts with in-plane image features and thus remove them. Combining this new method with the previous one allows us to remove potentially completely both in-plane and out-of-plane artifacts in photoacoustic imaging. We experimentally demonstrate this with experiments in phantoms as well as in vivo.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement 1. Introduction

Recent research has shown numerous potential clinical applications of photoacoustic imaging (PAI) [1–3]. This imaging technique is based on the photoacoustic (PA) effect. Samples are illuminated using short pulsed laser light. The local absorption of light generates ultrasound (US) waves which are then detected by a US transducer. PA images are reconstructed from the detected signals providing localized information about optical absorption properties of the samples. In clinical applications, the obtained information of endogenous chromophores such as hemoglobin helps diagnosing early stages of various diseases [2,4–6].

A typical PAI system basically consists of a light source and US transducer array. The transducer array can be classified as one-dimensional or dimensional [7]. While the two-dimensional transducer array provides 3D images, it requires significant users’ effort and experience for acquiring and interpreting these 3D images [8]. Additionally, two-dimensional transducer arrays and the associated scanners are unaffordable for many clinical applications [8]. One-dimensional transducer arrays, in contrast, are widely used for clinical studies [9,10], so from the point of view of clinical translation the incorporation of PAI in a one-dimensional array is preferred.

Several compact and low-cost PAI systems for clinical use have been developed. Integrating a laser source into a handheld US probe stands out among the approaches [9–12]. However, the occurrence of artifacts related to acoustic inhomogeneity of the tissue (clutter) is a major drawback of using a linear US transducer array. The artifacts aimed in this work include in-plane artifacts (IPAs), also called reflection artifacts, and out-of-plane artifacts (OPAs). While IPAs are caused by signals being reflected inside the imaging plane, OPAs are caused by absorbers located outside the imaging plane of the transducer array [13,14]. These artifacts appear as real image features, such as blood vessels, in the acquired image leading to misinterpretation. Therefore, correcting artifacts in PAI is of importance.

Previously, we proposed a method, photoacoustic-guided focused ultrasound (PAFUSion), to reduce IPAs [13,15]. This method has several limitations: a large number of US images are needed; the imaging plane has to be perpendicular to the PA sources; it is limited to reflectors within the angular aperture of the US probe. We then proposed another #364393 https://doi.org/10.1364/BOE.10.003124 Journal © 2019 Received 5 Apr 2019; revised 22 May 2019; accepted 27 May 2019; published 3 Jun 2019

method using method, the s response of th correlated wit due to longer signals than t However, bot Several ap Averaging (D employ tissue almost compl work for eas differently fro frames is requ In this wo US transducer at their origin sequence of i combine this and remove b We demo promising app 2. Theory 2.1. Artifacts Figure 1 illus which severa represents the absorbers in t imaging plan imaging plane array at that image as an O these absorbe appear in the absorber, upw absorber, dow Fi g multi-wavelen ample is image he image featu th the true im r propagation o their correspon th this method pproaches redu DCA) and Loc

e deformation etely OPAs. H ily deformable om in-plane fe uired in these m ork, we propos r array. During nal position rel images acquire method with oth IPAs and O onstrate the m

proach for corr s in photoacou

strates the prin al artifacts are e imaging plane this configurati e. Since the li e may absorb angle is high OPA, also calle ers are one or acquired imag ward dashed wnward dashed ig. 1. Artifacts in P ngth excitation ed with laser li ures. It is assum mage features o of reflected US nding real fea and PAFUSion ucing OPAs hav calized Vibrati

and motion tra However, these e tissue. Addi features might methods. e a method for g the displacem lative to the in ed during the d our previous w OPAs. method in pha recting artifacts ustic imaging nciple of artifa present in an e (orange plane ion: one is in th ight source ex

the light and g enough [19], ed direct OPA more acoustic ge: IPA, also blue arrow) a d yellow arrow)

PAI. (a) Configura

n [14] which ca ight at multiple med that the sp of their corresp S waves, IPAs atures. These f n do not work f ve been reporte ion Tagging ( acking to de-c methods stron itionally, the a be incorrect. r suppression o ment, OPAs mo nitial position o displacement c work using mu antoms and in s in PAI. facts in PAI. F n acquired im e in Fig. 1(a)) he imaging pla xcites a large v generate signal the out-of-pla A (downward y c reflectors. A called reflectio and indirect O ). ation resulting in ar an overcome th e wavelengths pectral respons ponding real a s appear at larg findings, there for OPAs. ed such as Def (LOVIT) [16– correlate the O ngly rely on m assumption tha Furthermore, of the OPAs by ove up while in of the transduc can therefore r ulti-wavelength n vivo. Result

Figure 1(a) sho mage, Fig. 1(b)

of the transduc ane and the oth volume, absorb ls. If the sensi ane absorber a ellow arrow). As a consequen on artifact (ref OPA (reflectio artifacts. (b) Acqui hose limitation to give the PA se of the IPAs absorbers. Add ger depths wit efore, can reve formation Com –18]. These ap OPAs. They can motion tracking

at OPAs are d a large numb y axially displ n-plane feature cer array. Com remove OPAs. th excitation to ts show that

ows a configu ). This acquire cer array. Ther her one is outsi ber located ou itivity of the tr appears in the In addition, un nce, two more

flection of the on of the out ired PA image. ns. In this A spectral is highly ditionally, th weaker eal IPAs. mpensated pproaches n remove and only deformed er of US lacing the es remain mparing a We then o identify this is a uration in ed image re are two ide of the utside the ransducer acquired nderneath e artifacts e in-plane t-of-plane

Only one configuration misinterpretat human spots cause indirect 3. Method 3.1. Transdu The principle image. OPAs transducer arr exact the sam the coordinate plane feature different beha Figure 2 s represent elev configuration seen in Fig. 2 Fig. 2(d). An Fig. 2 Acqui At the init plane absorbe acquired imag between the o appears at a s the same posi and potentiall Denote s and _' 2 o r = x feature (upwa . However, th tion. In clinica and birthmark t OPAs. ucer array disp e of this metho

appear at dep ray. If the tran me distance as t e system of th s stay at the s aviors, OPAs c shows a situatio vation, lateral, depicted in F 2(b). Another P OPA is presen 2. Displacing the ired PA image of c tial position of er and the tran ge, seen in Fig out-of-plane ab smaller depth t ition. Exploitin ly removed. as the shift of 2 (z_o z) + + Δ . s s

ard blue arrow) hree more vis al imaging, dire ks. Skin and bo placement od for identify pths equal to th sducer array is the displaceme he initial positi same position an be different on of one in-pl and axial axe ig. 2(a). The tr PA image is ac nt in these two e transducer array configuration (a). ( f the transduce nsducer is r. T g. 2(c). After bsorber and the than in the prev ng these differ f the OPA. It is s , therefore, is ( )Δ = Δ +z z x ) should be pre sible features ect OPAs migh one layers whi

ying OPAs is b he distance of s axially displa ent, while the O

ion (before the whereas the O tiated from the lane and one o s respectively. ransducer arra cquired at the acquired image

y. (a) Initial, and (d) Acquired PA im er array, Fig. 2 The OPA, ther

the transducer e transducer arr vious image. H rent behaviors, s determined as s a function of 2 2 2 o o o x +z − x + esent in the co which are a ht come from t ich are strong

based on their f their correspo aced up, in-pla OPAs experien e displacemen OPAs move u e in-plane featu out-of-plane ab . Figure 2(c) s ay is then axial new position o es. d (b) displaced mage of configura 2(a), the distan refore, appears r array being s ray is r'. Sinc However, the i , the OPA can s s r= + Δ −z r z Δ and can be 2 (z_o z) . + + Δ orrect PA imag artifacts might

the skin surfac acoustic refle r depth in the onding absorbe ane features mo nce in a lesser e nt) of the prob up. By exploit ures. bsorber. x, y, an shows a PA im lly shifted up of the transduc configuration. (c) ation (b). nce between th s at a depth of shifted up, the ce r'< + Δr z, in-plane featur n therefore be i ' r , where r= e rewritten as: ge of this t lead to ce such as ectors can acquired ers to the ove down extent. In e, the in-ting these nd z axes mage of a with zΔ , cer array, ) he out-of-f r in the e distance the OPA re stays at identified 2 2 o o x +z (1)

Images a displacement features and b and out-of-pla Fig. 3(c) resp where blue an in-plane featu In the se respectively. seen in Fig. 3 seen in Fig. completely th initial position image, therefo Fig. 3 respec the di To estima constant duri ( T) s zΔ =a. S Solving th Equation dependent on o x , and the ax Figure 4 i larger z is, to larger x is, t_o cquired along to have the s background. Th ane features on ectively. Figur nd yellow repre ure remains at t egmented ima Since the OPA (c), if we multi 3(d). The in-p he OPA in thi

n. At Δ = Δz z_T

ore it is comple

3. In-plane and out ctively in (a), as an

splacement. ate ΔzT, we as

ing the displa Substituting Eq Δ his equation, Δ (3) holds whe the OPA’s ax xial distance be llustrates the d the larger ΔzT the smaller Δz_T g the displace ame size. Ass he behavior al n the dashed an re 3(b) and Fig esent backgrou the same positi ages, backgrou A moves along iply segmented plane feature, s way, it is re T and beyond, etely removed. t-of-plane features n effect of the tran ssume that the acement. The q. (1) into this g 2 2 T o o z x z Δ + + T z Δ is determin

(

ial size, a, the etween the OP dependence of is; the larger T is, seen in Fi

ement are cro ume that they ong transducer nd dotted lines, g. 3(c) are a seq und and foregro ion while out-o und and foreg g zΔ into the d images all tog in contrast, is equired that th

, the entire OP .

s behavior, (b) and nsducer array disp e axial size of shift of the O gives: 2 ₍ o o x z − + + ned as:

(

T

z

Δ on the 3 p r a is, the larg

ig. 4(b).

opped out the y are perfectly r array displac , Fig. 3(a), is d quence of segm ound (features) of-plane feature ground pixels background o gether the OPA

not affected. he OPA compl PA is in the ba d (c) on the dotted placement. (d) OPA the OPA, a ( OPA at zΔ = 2 ) . T z a Δ =

)

ve. This equat ween the OPA nsducer array, z parameters. It ger ΔzT is, see

e extra part segmented to cement, zΔ , of depicted in Fig. mented acquire ). It can be see e moves up. s are set as of the previous A is gradually To be able to letely moves o ackground of t

d and dashed lines A correcting along (seen in Fig. 3

T

z

= Δ must th

tion shows tha and the imagi

o z . is worth notin en in Fig. 4(a) from the o separate f in-plane . 3(b) and ed images en that the 0 and 1 s images, removed, o remove out of its the initial s g 3(c)), is a en fulfill (2) (3) at ΔzT is ing plane, g that the ); and the

Above is a initial position are various O insufficient. O a certain zΔ w displacement corrected ima The meth algorithm, w segmentation To estima axial size of t in reality. T displacement section. 3.2. Combin We have prev wavelength e features can b Pearson corre revealed. How Fig. 4. D a study for a s n and another o PAs, Δz_T is th On the other ha when all OPAs since there is n age.

hod is summa which is based

step are the sa

Fig. 5. F ate ΔzT, we ha

the OPA does n This will be f

distance, zΔ ,

ing two metho

viously reporte excitation [14] be obtained. T elation coeffic wever, this m Dependence of Δz ingle OPA. Tw one at ΔzT or b hen different fo and, the larger

s are removed, no more OPA. arized in the d on the Sob ame as presente

Flow chart of the tr ave assumed th

not change dur further discus , might be lim

ods

ed a method fo . By imaging The spectral res cient [21]. Ba ethod does no T z on the 3 paramet wo images are beyond) to rem or different OP z

Δ is, the mor the corrected i At this point, following flo bel edge dete ed in [14].

ransducer array di hat acquired im

ring the displa sed in the D mited, this will

r identifying IP with multiple sponses are th ased on the co ot work for O ters zo, a , and x sufficient for move the OPA. PAs. In that cas re OPAs can be

image does no the process can ow chart, Fig.

ection algorith

isplacement metho mages are perfe acement. Howe Discussion sec also be addre PAs (reflection e wavelengths hen correlated t orrelation coe OPAs. In this w o x. this method (o . In a scenario se, two images

e removed. Ho ot change despi an be finished g . 5. The segm hm [20], and od. ectly segmente ever, this is not ction. In addi essed in the D n artifacts) usi , spectral resp to each other efficients, IPA work, we com one at the that there might be owever, at ite further giving the mentation the de-d ande-d the t the case ition, the Discussion ing multi-ponses of using the s can be mbine the

previous met complete meth The previ without segm significantly m segmentation segmentation It is neces wavelengths, since we com for segmentat 1. Image 2. Image 3. Segme 4. Correc 5. Correc 4. Setup Figure 6 sho commercial U in [14]. The U of 0.24 mm. I central 64 elem The handh Germany) wh (Laser2000, T ranges (680-9 coupled into GmbH, Germ generator (AF externally trig Fig. 6 optica In our exp used the mult energy at the

thod with the hod to remove ious method fo mentation does

more data proc is needed in in this work. sary to remind

the image that mbine the two m

tion. The two m with multiple along the trans ent the image a ct IPAs. ct OPAs giving ows the exper US scanner My US transducer It has a center ments were use held probe is m hich can trans The Netherlan 960 nm or 12 a custom-mad many) with a c

FG 3102, Tek gger the laser. T

6. Schematic drawi al fiber in a fixed p periments, lase ti-wavelength output of the fi transducer a e all artifacts in for IPAs can w s not rely on

cessing time an n the method d that in the pre

t had the highe methods, the im methods are com wavelengths ( sducer array di acquired at λn, g the final corre rimental setup yLabOne (Esao array in the ha frequency of 7 ed. mounted on a m slate along the nds) which can 00-2400 nm) de multi-mode ore diameter o ktronix, Germ The hardware c

ing of the setup w position. er light at wave method with fiber was 4.4 ± array displacem n PAI. work with or w the segmentat nd expensive c for OPAs d evious method, est signal was mage acquired mbined as the 1... n λ λ ). isplacement, Δ z Δ = 0 mm. ected image. . A handheld ote Europe BV, andheld probe 7.5 MHz with a motorized trans e axial axis. T n emit laser li at a repetition e optical fiber of 6.5 mm. Th many) is used components ar

with the probe mov elengths in the 8 wavelengths 0.4 mJ. ment method without segme ation algorithm calculation. On described abov , of the acquire selected for s d at the mutual following step z Δ , at a fixed w US array pro , The Netherla comprises 128 a bandwidth of slation stage (M The light sour

ight at a tunab n rate of 20 H bundle (Light he fiber bundle to synchroni re controlled by

vable in the vertica e range of 680 s rather than 4

described abo entation. The m, however, it n the other han ve. Therefore, ed images with segmentation. H l wavelength is ps: wavelength λ_n. obe is connec ands) which is p 8 elements wit f 66%. In our s MTS50A-Z8, T rce is an Opo able wavelengt Hz. Laser ligh tGuideOptics e is fixated. A ize the system y a LabVIEW

al position and the 0-960 nm was u 4 as in [14]. T ove as a approach t requires nd, image , we use h different However, s selected . cted to a presented th a pitch study, the Thorlabs, olette 532 th in two ht is then Germany A function m and to program. e used. We The pulse

5. Experime We performed experiment, a different wave at every 1 m reconstructed 100 pulses an Since PA of 790 nm. wavelengths w [14]. Compared to higher lase to represent e 5.1. Phantom A phantom w black thread Germany), as acoustic refle One absorber the probe’s ed a suspension wavelength o and US imag surfaces and r in-plane absor Underneath th IPAs of the in expected posi image are som of the out-of-p Fig. 7 elevat Figure 8 (reflections of retained (dash response of th ental results d experiments a total number elengths (720: mm of total 5 from 61 angle nd reconstructed images during The image was at 790 nm d to images in r pulse energy levation, latera m 1 was made of tw embedded in s shown in Fig ctor. Figure 7( was placed un dge to make su of 2% Intral f 750 nm [23] e is shown in reverberations rber (upward b he lid are som n-plane absorb ition relative to me features (do plane absorber

7. phantom 1. (a) tion view of the ex shows the res f the in-plane a hed blue arro hese features is in both phanto r of 13 image 10:790 nm), an mm displacem es (−30:1:30°) d using a Fouri g the transduce selected for s m instead of th our previous s . In order to pr al and axial axe wo black abso

agarose (1.5% g. 7(a). The pe (b) is a schema nderneath the p ure it was outsi lipid 20% in mimicking sc Fig. 7(c). The of the petri dis blue arrow) is v me more image ber. A direct O o the lid and th ownward dashe r. Two absorbers em xperiment configur sult of correcti absorber) are r ws, seen in F s strongly affec

oms and in viv es were acquir nd 5 PA image ment of the tra

of plane waves ier transform b r array displac segmentation he image givin study [14], the resent the resul es respectively orbers with a th

%) in a petri etri dish lid of atic elevation v probe and the o ide of the imag demi-water w cattering in hum e gray color pa

sh lid. The hot visualized at an e features (upw OPA (downwar he in-plane abso ed yellow arro mbedded in agaros ration. (c) Combin ing IPAs usin removed. Howe Fig. 8(b)). The cted by noise w vo to demonstr red: 1 US ima es at a wavelen ansducer array s. The PA ima based reconstru cement were ac in correcting ng the highest signal-to-nois lts, x, y, and z y. thickness of ~0 i dish lid (Gr f 750 mμ thi view of the exp other one was p ging plane. The with estimated

man soft tissue art represents t color part is t n expected pos ward dashed b rd yellow arro orber. At the b w) which are p se in a petri dish ned PA and US im ng 8 waveleng ever, some par e reason woul which was also

rate the method age, 8 PA ima ngth of 790 nm y). The US im ages were avera uction algorithm cquired at a wa

IPAs using t signal as pre se ratio is impr axes are used i

0.7 mm cut ou reiner Bio-One ickness was us periment confi positioned ~3-e coupling m~3-ed d ' s

μ

probably indire lid. (b) Schematic mage. gths. Most of rts of the reflec ld be that the pointed out in d. In each ages at 8 (1 image mage was aged over m [22]. avelength multiple esented in roved due in figures ut from a e GmbH, sed as an figuration. 4 mm off dium was cm−1 at a bined PA showing where the to the lid. which are ible at an ner of the ect OPAs c the IPAs ctions are e spectral n [14].

The acqui 9(b) and Fig. marked with Visualization shift of indir transducer ar resolution of also slightly imaging plan transducer arr Fig. 9 respec Figure 10 displacement, is not comple Fig. 8. Correc ired image was 9(c) respectiv a dotted line 1. The out-of rect OPAs is l rray. The size

the transducer changes. The e or the displa ray produced so

9. In-plane and out ctively in (a), as an 0 presents th , see also in Vi etely removed s

cting IPAs. (a) Ac s then OPA co vely show the

and a dashed f-plane feature less pronounc e of in-plane r array along d reason might acement was n ome positionin t-of-plane features n effect of the tran he result of isualization 2. showing that Δ cquired PA image. orrected with d behavior of th line in Fig. 9 es clearly shift ed since they features can depth. In addit be that the in not perfectly a ng error. s behavior, (b) and nsducer array displ

correcting OP An OPA (dash T z Δ should be (b) IPA corrected displacing the t he out-of-plan 9(a) during the

t up, especiall have a large be seen chan tion, the positi n-plane absorb along the axial

d (c) on the dotted lacement. PAs with 5

hed yellow arr e larger than 5 d image. transducer arra ne and in-plane e displacement ly the direct O er axial distanc nging due to ion of in-plane ber was slightl

l axis or displ

d and dashed lines mm transduc row, seen in Fi mm for this O ay. Figure e features t, seen in OPA. The ce to the the axial e features ly off the acing the s cer array ig. 10(b)) OPA to be

completely re of them is con Figure 11 the acquired P final corrected The displa more IPAs a displacement. 5.2. Phantom Another phan with an in-pla agarose (1.5% used as the co where one ab emoved. Additi nsequently ove Fig. 10. Correc shows the resu PA image, Fig d image with o acing transduce are removed . m 2

ntom was made ane feature. Tw %) in a petri dis oupling medium bsorber was in ionally, as ther ercorrected.

cting OPAs. (a) Ac ult of combinin . 11(a), most o only one true in er array metho using this m Fig. 11. F e to mimic a s wo black absor sh, Fig. 12(a). m. The two abs

the imaging p

re is some sligh

cquired PA image ng the two met of all artifacts a n-plane absorbe od does not wo

ethod due to

Final corrected ima

situation that a rbers (the same A solution of 2 sorbers were p plane and the o

ht movement o

. (b) OPA correcte thods for IPA are identified a er, seen in Fig. rk for IPAs. In their slight age. an OPA appea e as in phantom 2% Intralipid 2 placed under th

other one was

of in-plane feat

ed image. and OPA. Com and removed g . 11(b). n this experime

movement du

ars at the same m 1) were emb 20% in demi-w he probe as in F outside of the tures, part mpared to giving the ent, a few uring the e position bedded in water was Fig. 12(b) e imaging

plane and the acquired PA i Figure 12( position, zΔ feature is a re starts to appe at the same po Fig. 1 Schem displa Figure 13( The OPA, as position. Thes in the acquir artifact marke of reconstruc curved, its par removed, seen

Fig. 1 of fea

distance of the image, the two

(c) shows the P = 0 mm, only econstruction a ar and move u osition, seen in 12. Phantom 2. ( matic elevation vi acement.

(a) presents the expected, mov se behaviors ca red image. Sim ed in the dotted tion artifact is rt in the dotted n in Fig. 13(a). 3. OPA corrected atures in dashed an ese absorbers t absorbers app PA images dur y one feature artifact. When up along with t n Visualization

(a) Two black ab iew of the experi e result of this ves up away fro an be seen clea milarly, Fig. 1 d blue line in t s dependent on d line moves do .

d image. (a) Acqui nd dotted lines in (a

to the probe wa eared at the sa ring the transdu

is visible. Th the probe is li the displaceme 3. bsorbers embedde iment configuratio experiment wh om the real in-arly in Fig. 13 3(c) illustrates the acquired im n the distance own. As a resu

ired and OPA corr a) respectively.

as the same. A ame position.

ucer array disp he curve at the ifted up, zΔ > ent while the in

ed in agarose in on. (c) Acquired here the OPA -plane feature w

(b) which show s the behavior mage. It is wor to the probe. ult, this reconst

rected images. (b) As a consequen placement. At t e same positio > 0 mm, anothe n-plane feature a petri dish. (b) images along the is completely which stays at ws the dashed r of the recon rth noting that Since it beco truction artifact ) and (c) Behavior nce, in the the initial on of this er feature e remains ) e removed. the same blue line nstruction this type omes less t is partly r

It is notab a single featu of the two fea feature canno

5.3. In vivo

We also asses at fingers whi the imaged fin Figure 14 an ink mark w 14(a). The w experiments. Acquired principle, they and bone lay However, in t the position in OPAs. Fig. 1 spot. ( Figure 15 most of the I reported in [

ble that when an re, the recorde atures. Though t be recovered ss our method

ich would give nger to mimic

shows the con was placed in white line in

The ink mark w PA and US im y show similar yers while the

this PA image n agreement w 14. In vivo experim (c) Acquired PA im shows IPA an IPAs caused b 14]. Two featu n OPA appear ed amplitude o h the OPA can

. This will be d with in vivo ex e clear IPAs. I a human spot w nfiguration of a water tank Fig. 14(b) d was a few mill mages are pres r structures as o PA image sh e, there are a f ith the ink mar

ment. (a) Experim mage. (d) Acquire nd OPA correc by the bone la ures, f1 and f s at the same p f that overlapp be removed, t discussed furth xperiments. In In addition, we which would g the in vivo ex filled with dem depicts approx

imeters outside sented in Fig. observed in [14 hows the skin,

few more featu rk’s position on

ment configuration d US image. cted images. In

ayers are rem f2 in Fig. 15(a

position with a ped feature is a the true amplitu her in the Discu n these in vivo e put a black i give strong OPA

xperiments. A mi-water as a ximately the e of the imagin 14(c) and Fig 4] where the U superficial bl ures (dashed b n the skin surfa

n. (b) Ink mark mi n the IPA corre moved. This m a), are of inte

a real in-plane f a sum of the am tude of the real ussion section. experiments, w ink mark on th As [25]. volunteer’s fin coupling med imaging plan ng plane. g. 14(d) respec US image show lood vessels a lue circle) app face. These are

imicking a human ected image, F atches with th erest. Observat feature as mplitudes l in-plane we aimed he skin of nger with dium, Fig. ne in the ctively. In s the skin and IPAs. pearing at probably n ig. 15(b), he results tion from

correcting OP up and f2 rem by the skin, w is not, in the I close to the n can be seen i transducer arr as mentioned reconstruction addition, som reason is tha displacement. Combing 6. Discussio Correcting art various advan needed as in require a larg Third of all, tissue as in D PA shows that mains at the sam

while f2 could b IPA corrected i oise. On the ot in Visualizatio ray displaceme d in the phanto n artifact curve me real in-plan t there was so . correcting IPA on tifacts is of im ntages over pr the case of L ge number of U deforming tiss CA. As a cons f1 is an OPA me position). f1 be a reflection image. The rea ther hand, ther on 4 in which ent. Therefore, om 2 experim e of the skin is ne features are ome slight mo Fig. 15. IPA a A and OPA give

Fig. 16. F mportance for r reviously repo LOVIT and DC US images whi sue by applyin sequence, the d and f2 is an IP 1 could be an in

of the skin sig ason can be tha re are some OP some features the true intens ment, and thus

s also removed party remove ovement of th

and OPA corrected es the final cor

Final corrected ima

reliable imagin orted methods. CA [16–18]. S ile no US imag ng force migh detected signals PA (seen in Vi ndirect OPA o gnal. However, at the intensity PAs at the sam s move up aw sity of f2 is aff it is not iden d with the OPA ed in the OPA he finger durin d images. rrected image, age. ng. Our new m . First of all, Second of all, ge is needed in ht affect the op s might not rep

isualization 4, of the ink mark , f1 is removed of these two f me position with way from f2 d fected by these ntified as an I A correcting m A corrected im ng the transdu Fig. 16.

method for OPA no motion tra these existing n the proposed ptical properti present truly th f1 moves k reflected d while f2 features is h f2. This during the e features, IPA. The method. In mage. The ucer array As offers acking is methods d method. ies of the he source.

In the proposed method, correcting OPAs can be done without tissue deformation. Lastly, deforming tissue by focusing strong pressure, as is done in LOVIT, might violate US safety.

The proposed method for correcting OPA relies on segmentation. In this work, a simple segmentation approach was used for a low calculation cost and time consumption. However, it might not segment images properly giving inaccurate axial dimension of OPAs. As a result, OPAs might not be correctly removed. Over-segmentation also might happen as pointed out in [14]. A more effective segmentation algorithm should be considered for a better performance.

In our experiments, while the probe was displaced, the fiber bundle remained fixated. The purpose of this was to maintain the signal strength of image features. However, if the light source is displaced with the probe, the laser beam is also repositioned and thus excites different tissue volumes. Acquired images along the displacement might show different structures resulting in miscorrecting.

In clinical applications, the displacement distance, Δz, might be limited. Depending on the location of out-of-plane absorbers, ΔzT might not be achieved, as discussed in section 3.1. As a result, OPAs are not completely removed, in which case another approach is needed. However, our results show that within 5 mm displacement, OPAs can be completely removed for a large range of locations and axial dimensions.

In this work, out-of-plane absorbers were positioned outside of the imaging plane elevationally. Lateral out-of-plane absorbers can also cause OPAs. The principle of the proposed method still holds for these OPAs. The quantity x_o, used to describe out-of-plane absorbers in section 3.1, in this situation will be the lateral distance between the out-of-plane absorber and the imaging plane. Therefore, lateral OPAs can be identified and removed.

Axially displacing the probe in essence is to adjust the distance between the probe and in-plane and out-of-in-plane absorbers. Displacing the probe in other directions might also be able to de-correlate in-plane and out-of-plane image features. In a configuration as shown in phantom 1, if the probe is elevationally displaced in the direction to the out-of-plane absorber, real in-plane features will move down and OPAs will move up. However, in a scenario that there is another out-of-plane absorber in the other side of the imaging plane. OPAs of this out-of-plane absorber will also move down. Therefore, elevationally displacing the probe in both directions is required to identify all OPAs. The amount of work is double compared to using axial displacement. Nevertheless, displacing the probe in other manners and comparing with the proposed method will be investigated in our future work.

In a situation that an OPA appears at the same position with a real in-plane feature as a single feature, the image value is a sum of the OPA and the in-plane features. Displacing the transducer array can separate these two and remove the OPA. However, true image value of the real in-plane feature cannot be recovered. Interpolating or extrapolating image values of the OPA along the displacement might be able to estimate its value at the superposition. True image value of the real in-plane feature can, therefore, be recovered. Additionally, in the proposed method, OPAs are removed by setting their pixel values to 0. This might also remove the background information behind the OPAs. If the image value of the OPAs can be estimated, the background information can be retained while removing OPAs by subtracting the recorded value by the estimated one. This will be investigated in our future work.

In this work, the volunteer had to keep the finger still for ~5 minutes. Slight movements were inevitable resulting in some miscorrection. This is not ideal for clinical applications. However, the long experiment time was due to technical limitations. In particular, the translating stage was slow. It took ~2 minutes to acquire in total 5 PA images along 5 mm of the probe displacement. Using a higher speed translating stage will significantly reduce the experiment time. The acquiring data process with 8 wavelengths took ~2 minutes. This was due to the laser pulse repetition rate of 20 Hz. Using a high repetition rate laser would potentially achieve real-time artifact correction as shown in [14].

7. Conclusion

We have proposed a new method to remove out-of-plane artifacts exploiting different behaviors of out-of-plane artifacts and in-plane image features by axially displacing the transducer array. Combining this new method with our previous method for in-plane artifacts using multiple wavelengths [14], in-plane and out-of-plane artifacts in photoacoustic imaging can be identified and thus removed. Experiments in phantoms and in vivo were carried out to evaluate the combination of the two methods as a proof of concept. Results show the potential of this combined method for providing true photoacoustic images with no ultrasound images needed. In addition, a handheld probe suitable for clinical applications was used in the experiments bringing this method a step forwards to clinical translation.

Funding

European Union’s Horizon 2020 research and innovation programme (CVENT, H2020-731771)

Acknowledgements

The authors would like to thank Altaf Hussain for his discussions. The authors are also grateful to Johan van Hespen and Tom Knop, University of Twente, for their help with experiment setup.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. K. S. Valluru, K. E. Wilson, and J. K. Willmann, “Photoacoustic Imaging in oncology: translational preclinical and early clinical experience,” Radiology 280(2), 332–349 (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. Heijblom, W. Steenbergen, and S. Manohar, “Clinical photoacoustic breast imaging: the Twente experience,” IEEE Pulse 6(3), 42–46 (2015).

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

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

7. B. W. Drinkwater and P. D. Wilcox, “Ultrasonic arrays for non-destructive evaluation: A review,” NDT Int.

39(7), 525–541 (2006).

8. C. D. Herickhoff, M. R. Morgan, J. S. Broder, and J. J. Dahl, “Low-cost volumetric ultrasound by augmentation of 2D systems: Design and prototype,” Ultrason. Imaging 40(1), 35–48 (2018).

9. M. K. A. Singh, W. Steenbergen, and S. Manohar, “Handheld probe-based dual mode ultrasound/photoacoustics for biomedical imaging,” in Frontiers in Biophotonics for Translational Medicine(Springer, 2016), pp. 209–247. 10. 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).

11. 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). 12. 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).

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

14. H. N. Y. Nguyen, A. Hussain, and W. Steenbergen, “Reflection artifact identification in photoacoustic imaging using multi-wavelength excitation,” Biomed. Opt. Express 9(10), 4613–4630 (2018).

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

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

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

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

19. M. F. Beckmann, H.-M. Schwab, and G. Schmitz, “Optimizing a single-sided reflection mode photoacoustic setup for clinical imaging,” in Ultrasonics Symposium (IUS), 2015 IEEE International (IEEE, 2015), pp. 1–4. 20. I. Sobel, and G. Feldman, “A 3x3 isotropic gradient operator for image processing,” a talk at the Stanford

Artificial Project in, 271–272 (1968).

21. J. Benesty, J. Chen, Y. Huang, and I. Cohen, “Pearson correlation coefficient,” in Noise Reduction in Speech Processing(Springer, 2009), pp. 1–4.

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

23. R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express 16(8), 5907–5925 (2008).

24. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). 25. S. Preisser, G. Held, H. G. Akarçay, M. Jaeger, and M. Frenz, “Study of clutter origin in in-vivo epi-optoacoustic