- Email: [email protected]
Reconstruction of optical absorption coefficient distribution in intravascular photoacoustic imaging
Computers in Biology and Medicine 97 (2018) 37–49
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/compbiomed
Reconstruction of optical absorption coefficient distribution in intravascular photoacoustic imaging Sun Zheng *, Zheng Lan Department of Electronic and Communication Engineering, North China Electric Power University, Baoding 071003, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Intravascular photoacoustics (IVPA) Inverse problem Quantitative intravascular photoacoustics (qIVPA) Optical absorption coefficient Optical inversion
Intravascular photoacoustic (IVPA) imaging can discriminate between normal arterial tissues and atheromatous plaques in images, particularly images of lipid-rich plaques with high spatial resolution and optical contrast. However, conventional IVPA only recovers the deposited optical energy that is the product of the tissue optical absorption coefficient and local optical fluence. Herein, a one-step model-based method for single-wavelength IVPA imaging system is proposed. The proposed method directly reconstructs the optical absorption coefficient from the boundary measurement of acoustic pressure. To obtain the theoretical optical deposition, light illumination and transport in multilayered vessel-wall tissues are numerically modelled with the diffusion approximation to radiative transfer equation. Then, the generation and propagation of photoacoustic (PA) waves in acoustically homogeneous vessel-wall tissues are modelled by using the PA wave equation. The theoretical acoustic pressure series is thus obtained. Finally, the optical absorption coefficient is iteratively updated from an initial guess by minimising a nonlinear least-square error function to quantify the difference between the measured and theoretical acoustic pressure with split Bregman optimisation based on total-variation regularisation. The numerical simulation experiments demonstrated that optical absorption coefficient maps can be directly recovered on the basis of the acoustic boundary measurement of vascular cross-sections. Compared with the optical deposition map, the optical absorption coefficient map can provide more reliable qualitative and quantitative information on the morphology and optical properties of the imaged arteries. The proposed method enables the accurate and reliable evaluation of atheromatous lesions and the early identification of vulnerable plaques.
1. Introduction Intravascular photoacoustic (IVPA) imaging is a typical application of endoscopic photoacoustic tomography (PAT). It has recently emerged as a complementary imaging modality to intravascular ultrasound (IVUS) for the accurate evaluation of arterial atherosclerotic lesions and the effective identification of vulnerable plaques [1]. This approach can visualise the structure and functional composition of arterial vessel walls and atheromatous plaques on the basis of the characteristics of tissue optical absorption. During IVPA intervention, a catheter integrated with a multimode optical fibre, an optical reflection element and an ultrasonic detector is inserted and pushed to the distal end of the target vascular lumen. During the drawback of the catheter, a probe mounted on the catheter tip emits modulated light that is absorbed by the surrounding tissues and then converted to heat. Periodic heat flow causes the thermal expansion of
tissues, thereby exciting a wide band of ultrasonic waves, i.e., photoacoustic (PA) waves. The PA waves are transmitted rapidly to the surface and are finally collected by the ultrasonic detector. The distribution of the absorbed optical energy density can be reconstructed from the timedependent boundary measurement of acoustic pressure through acoustic inversion. To reveal the functional components of normal and pathological tissues, the optical parameters of the tissues are recovered by solving the optical inverse problem [2]. IVPA exhibits the high spatial resolution of IVUS for imaging deep tissues [3] and the main advantages of other optical imaging techniques, such as limited tissue damage, excellent optical contrast and functional and molecular imaging. Important progress has been made in improving thermal IVPA [4], spectroscopic IVPA [5], plasmonic IVPA [6] and high-speed IVPA [7] systems. Moreover, integrated IVPA/IVUS imaging catheter [8] and novel ultrasonic detectors [9] have been developed and utilised in the characterisation of lipid-rich plaques [10] and the detection of
* Corresponding author. P.O. Box 21, North China Electric Power University, Baoding 071003, China. E-mail address: [email protected] (S. Zheng). https://doi.org/10.1016/j.compbiomed.2018.04.012 Received 27 December 2017; Received in revised form 16 April 2018; Accepted 16 April 2018 0010-4825/© 2018 Elsevier Ltd. All rights reserved.
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
by Ref. [30]. In this setting, the pulsed laser source and an ultrasonic transducer array are mounted on a same rotating frame, and the distribution map of the optical absorption is directly recovered from the measured acoustic pressure series. A forward model was constructed by Ref. [31] to combine an optical transport model with an acoustic propagation model. This model simultaneously and directly estimates absorption and scattering coefficients from the PA pressure time series through a Bayesian approach. The scanning aperture of IVPA is enclosed in a geometrically complex lumen. Well-developed qPAT methods cannot be directly applied in qIVPA. This study aims to design a one-step qIVPA method for an IVPA system with single optical excitation at a single wavelength. The transversal map of the absorption coefficient distribution on a vascular crosssection was directly recovered from the acoustic pressure time series generated by multilayered arterial wall tissues and atherosclerotic plaques and then measured using the ultrasonic detector. A forward model based on diffusion approximation (DA) to the radiative transfer equation (RTE) was constructed to simulate theoretical optical deposition in acoustically homogenous tissues with known scattering and Gruneisen coefficients. The PA pressure was modelled by explicitly solving the PA wave equation with a finite difference time domain (FDTD) algorithm. The absorption coefficient of a vascular cross-section was iteratively updated from an initial guess by minimising a squared error function that quantifies the difference between the measured and modelled acoustic pressure series with split Bregman optimisation based on total variation (TV) regularisation. The performance of the proposed method was demonstrated with computer-simulated arterial vessel-mimicking phantoms. The ability of the Monte Carlo (MC) simulation to model light transport was compared with that of DA in accordance with the validation results. The accuracy of L–M and split Bregman optimisation based on TV- and L2-regularisation was also evaluated and compared.
macrophages in atherosclerotic inflammation [11]. However, the integration, mobility, operability, imaging speed and reliability of the IVPA imaging system still cannot meet the standards of clinical applications. The acoustic inverse problem in IVPA involves the reconstruction of the initial acoustic pressure or the deposited optical energy from the time-dependent boundary measurement of acoustic pressure by using filtered back projection, time reversal and iterative algorithms [12–15]. Optical inversion, the so-called quantitative IVPA (qIVPA), is essentially the inversion of a light transport operator and involves the recovery of optical property parameters, particularly the absorption and scattering coefficient [16–18] and thermal expansion parameter (Gruneisen coefficient) [19], from the time-dependent measurement of acoustic pressure or optical deposition recovered through acoustic inversion. The acoustic and optical inverse problems can be solved separately because the time of acoustic propagation in biological soft tissues is approximately three orders of magnitude longer than that of laser illumination and optical energy absorption by tissues [19]. The optical absorption is determined by the local absorption coefficient and light fluence and cannot directly reflect the optical properties of the tissues, which are closely related to chemical composition. The optical properties of normal tissues drastically differ from those of diseased tissues. Therefore, qIVPA can provide more accurate and reliable evaluation for the early identification of vulnerable atherosclerotic plaques than ordinary IVPA. However, to the best of our knowledge, no special methods have been proposed for qIVPA thus far. Current quantitative PAT (qPAT) methods often involve two steps [20–23]. Firstly, the distribution map of initial acoustic pressure is reconstructed from acoustic boundary measurements through acoustic inversion. Optical deposition measurements are determined on the basis of PA images. Secondly, a forward model of light transport in the tissues is established to obtain the theoretical optical deposition. Finally, the map of the optical absorption and scattering coefficient is recovered by minimising the squared error between the measured and theoretical optical deposition. Multiple-illumination PAT can alleviate the ill-posedness of the optical inversion caused by absorption-scattering nonuniqueness. Furthermore, the distribution map of the speed of sound in an acoustically inhomogeneous medium can be simultaneously recovered [20]. The main limitation of two-step qPAT methods is that the accuracy of the optical map highly depends on the optical deposition reconstructed from the boundary PA measurements. Errors are inevitably generated during acoustic inversion given the limited bandwidth of the ultrasonic detector in contrast to the ultrawide PA frequency response; the orientation, size and shape of the imaging target with respect to the detecting aperture; the effect of optical scattering and the acoustic inhomogeneity and complexity of the medium. In addition, the qPAT two-step methods are limited by the cost of imaging systems based on an optical method, e.g., diffuse optical tomography [24], optical contrast agent [25] or acousto–optic modulation [26] to measure the optical fluence [27]. Furthermore, several artefacts in the optical deposition map may reduce the accuracy of optical parameters. The ability of direct one-step optical map recovery from PA boundary measurements to resolve the above issues has been studied recently [27–31]. For example, a calibration-free qPAT method proposed by Ref. [27] needs only the normalized acoustic boundary measurements and the laser source intensity as the input data. These requirements eliminate the influence of detector sensitivity. The method reported by Ref. [28] allows the simultaneous reconstruction of the absorption coefficient and speed of sound from multiple acoustic pressure datasets generated from multiple optical illuminations. In regions where the changes in the absorption coefficient and speed of sound are relatively small, the perturbation to the absorption coefficient around the known backgrounds is calculated with Born approximation and Landweber iteration. The fully nonlinear inverse problem is solved with an optimal control method, where the objective function is minimised with the Levenberg–Marquardt (L–M) algorithm to obtain the optimal solution of the absorption coefficient. A rotating measurement setting was designed
2. Method 2.1. Construction of the imaging model Geometrical models of arterial cross-sections containing atherosclerotic plaques are constructed to illustrate the proposed method. The models are shown in Fig. 1. The cross-section is divided into M uniform sectors, and the mth measuring angle of the detector is θm ¼ 360(i1)/M (see Fig. 2). The angles of the corresponding imaging region are in the range of [θma, θmb], where m ¼ 1,2, …,M, θma ¼ θm 180/M and θmb ¼ θm þ 180/M. A θ‒l polar coordinate system centred at the lumen is built. In this system, θ and l represent the polar angle and polar radius, respectively. As shown in Fig. 3, the ultrasonic detector mounted on the catheter tip is located at the lumen centre and collects PA waves generated by the surrounding tissues along a circle parallel to the imaging plane in the closed lumen [32]. The lumen, vessel-wall intima/media, plaques and adventitia surround the catheter centre. The detector is simplified to be an ideal point-like detector, and the aperture effect related to its size is ignored. 2.2. Simulation of modelled PA pressure 2.2.1. Formation of the forward problem A single short-pulsed laser at a single wavelength illuminates the vessel wall from the lumen centre. The optical energy absorbed by the surrounding tissues is expressed as follows: HðrÞ ¼ μa ðrÞΦðr; μa ðrÞ; μs ðrÞÞ;
(1)
where r 2 Ω represents a point in a bounded imaging domain Ω with a smooth boundary ∂Ω; H(r) is the absorbed energy density at r; μa(r) and μs(r) are the wavelength-dependent optical absorption coefficient and scattering coefficient at r satisfying the Lipschitz continuity and Φ is the 38
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Fig. 1. Cross-sectional geometry of arterial vessels including (a) thin-cap fibroatheroma (TCFA) and lipid plaques, (b) TCFA and calcified plaques, and (c) TCFA, fibrofatty and calcified plaques.
μ0s ðrÞ≫μa ðrÞ;
(3)
the RTE is simplified as a diffusion equation (DE) represented with spherical harmonics where
μs 'ðrÞ ¼ μs ðrÞð1 gÞ;
(4)
is the reduced scattering coefficient at r, and g is the anisotropy factor. The time-independent DE satisfying the Robin boundary condition is represented as follows [17,34]: r⋅kðrÞrΦðr; tÞ þ μa ðrÞΦðr; tÞ ¼ Qðr; tÞ;
(5)
where Φ(r, t) is the light fluence at r and time t, Q(r, t) is the light source and k(r) is the optical diffusion coefficient. k(r) is expressed as follows:
Fig. 2. Multi-layered vessel wall tissue in a θ‒l polar coordinate system.
kðrÞ ¼
1 1 : 3½μa ðrÞ þ μs 'ðrÞ 3μ0s ðrÞ
(6)
The time-domain Robin boundary condition is written as follows [17]: cΦðr; tÞ þ 2kðrÞ
1 þ Rf _ n ⋅rΦðr; tÞ ¼ 0; 1 Rf
(7)
_
where r 2 ∂Ω, n is the outer normal vector of the tissue surface at ∂Ω, c is the speed of light in tissues and Rf is the internal reflection coefficient of diffuse transmission. Rf is approximated with a polynomial in accordance with Fresnel's law of reflection, which is expressed as follows: Rf 1:4399n2 þ 0:7099n1 þ 0:6681 þ 0:0636n;
where n is the refractive index of the target tissue relative to the background. A narrow beam of collimated light is normally incident upon the tissue surface and is approximated as an internal photon flux that propagates at the tissue boundary. The directivity of the light basically disappears when the penetration depth is approximately one photon free path in length because of the strong scattering effect of tissues. Light propagation becomes isotropic, and the source is approximately an undirected point at the penetration depth. The point-like source is represented with a δ-function with respect to time and space:
Fig. 3. Schematics of IVPA imaging.
light fluence (i.e., optical radiance energy) related to r, μa(r) and μs(r). The initial acoustic pressure at r is determined as follows: b pðrÞ ¼ ΓðrÞHðrÞ;
(8)
(2)
b is the conversion efficiency of heat energy to pressure and where Γ represents the thermodynamic properties of the medium. The acoustic pressure time series p(r, t) is finally measured by the detector during rotation inside the lumen.
Qðr; tÞ ¼ δðr r s ; tÞ;
(9)
where rs is the location of the light source in the θ‒l polar coordinate system. The collimated light source model is included in the Robin boundary condition in Eq. (7):
2.2.2. Forward model of light transport Light transport in turbid media, such as biological soft tissues, is usually modelled by the RTE, an integro–differential equation describing photon distribution in the phase space [33]. RTE can accurately describe photon migration in tissues, especially in nondiffusion regions. The DA is usually used to simplify the solution to the RTE because the direct solution may be highly computationally complex. Assuming
cΦðr; tÞ þ 2kðrÞ
1 þ Rf _ n ⋅rΦðr; tÞ ¼ 4Qs wðrÞδðtÞ; 1 Rf
8r 2 ∂Ω;
(10)
where Qs is the intensity of the light source, and w(r) is the weighting 39
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Fig. 4. Meshing of the imaging region. (a) Interior point; (b) Boundary point; (c) Boundary point facing to the light source.
location function, wðrÞ ¼ δðr r s Þ:
where q is the complex frequency factor of LT, Φ(r,rs,q) is the LT of the time-domain optical radiance energy at r and time t and ∂Ω0 is the tissue surface facing the source. If q ¼ 0, DE is a steady-state diffusion equation.
(11)
2.2.3. FDTD expression of DE The generalised pulse spectrum technique [35] is applied to simplify the solution of DE. The Laplace transform (LT) of a time-domain curve at a complex frequency is used instead of the whole time-domain curve. The forward problem is then transformed into a two-dimensional (2-D) problem that is time-related but is not time-dependent. As mentioned above, the collimated light source model is integrated into the Robin boundary condition. Thus, the source term of the original DE is 0. The DE, Robin boundary condition and Robin boundary condition including the source item are expressed in complex frequency domain as: 8 r⋅kðrÞrΦðr; r s ; qÞ ½c⋅μa ðrÞ þ qΦðr; r s ; qÞ ¼ 0 > > > > > > < cΦðr; r ; qÞ ¼ 2kðrÞ 1 þ Rf _ n ⋅rΦðr; r s ; qÞ s 1 Rf > > > > 1 þ Rf _ > > n ⋅rΦðr; r s ; qÞ ¼ 4Qs δðr r s Þ : cΦðr; r s ; qÞ þ 2kðrÞ 1 Rf
2.2.4. FDTD solution of optical deposition As shown in Fig. 4, a smooth region D 2 Ω is meshed with unstructured rectangular grids represented by nodes in the θ‒l coordinate system, where hθ and hl are the steps along the positive θ- and l-axis, respectively, and (ihθ, jhl) is the coordinate of a point P 2 D. A relatively narrow beam of collimated light propagates along the positive l-axis. Both ends of the first equation in Eq. (12) are integrated in D to yield the following equation: ∫k
∂D
r2Ω r 2 ∂Ω
∂Φ _ ds ∬ ðμa c þ qÞΦdθdl ¼ 0; ∂n D
8D 2 Ω:
(13)
where s2∂Ω. Then, the optical radiance energy at each node is obtained by solving the FDTD expression of DE. In the case of an interior node P 2 D in Fig. 4(a), the differential expression of Eq. (13) is expressed as follows:
;
r 2 ∂ Ω0 (12)
Table 1 Parameters of vessel phantom a. Tissue name
Tissue composition
Average refractive index
Absorption coefficient (cm1)
Scattering coefficient (cm1)
Anisotropy factor
Thickness (cm)
Adventitia Intima/Media Lipid pool TCFA Lumen
Connective tissues Muscular tissues Lipid Fibrous tissues Blood
1.41 1.41 1.46 1.46 1.34
0.20 0.20 0.60 0.01 0.99
20 20 500 80 800
0.80 0.80 0.80 0.75 0.99
0.1330 0.2520 0.2267 (middle) 0.0100 0.1150
Table 2 Parameters of vessel phantom b. Tissue name
Tissue composition
Average refractive index
Absorption coefficient (cm1)
Scattering coefficient (cm1)
Anisotropy factor
Thickness (cm)
Adventitia Intima/Media Calcified plaque TCFA Lumen
Connective tissue Muscular tissue Calcium Fibrous tissue Blood
1.41 1.41 1.46 1.46 1.34
0.20 0.20 0.05 0.01 0.99
20 20 80 80 800
0.80 0.80 0.80 0.75 0.99
0.1250 0.2600 0.2175 (middle) 0.0100 0.1150
Table 3 Parameters of vessel phantom c. Tissue name
Tissue composition
Average refractive index
Absorption coefficient (cm1)
Scattering coefficient (cm1)
Anisotropy factor
Thickness (cm)
Adventitia Intima/Media Fibro-fatty plaque Calcified plaque TCFA Lumen
Connective tissue Muscular tissue Fibro-fatty tissue Calcium Fibrous tissue Blood
1.41 1.41 1.46 1.46 1.46 1.34
0.20 0.20 0.45 0.05 0.01 0.99
20 20 80 80 80 800
0.80 0.80 0.80 0.80 0.75 0.99
0.1421 0.2150 0.0619 (middle) 0.0715 (middle) 0.0100 0.1429
40
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
2.2.5. FDTD solution of PA pressure The lossless propagation of PA waves in an acoustically homogeneous medium is physically described with the PA wave equation as [13,36]:
hl hθ hl kði þ 1; jÞΦði þ 1; jÞ þ kði; j þ 1ÞΦði; j þ 1Þ þ kði 1; jÞΦði 1; jÞ hθ hl hθ hθ hl hθ hl þ kði; j 1ÞΦði; j 1Þ kði þ 1; jÞ þ kði; j þ 1Þ þ kði 1; jÞ hl hθ hl hθ hθ þ kði; j 1Þ þ 2½μa ði; jÞc þ qhθ hl Φði; jÞ hl
r2 pðr; tÞ
¼ 0;
where Φ(i, j) is the optical radiance energy at the node (ihθ, jhl); μa(i, j) is the absorption coefficient at (ihθ, jhl) and k(iþ1, j), k(i, jþ1), k(i‒1, j) and k(i, j‒1) are the optical diffusion coefficient at ((iþ1)hθ, jhl), (ihθ, (jþ1)hl), ((i‒1)hθ, jhl) and (ihθ, (j‒1)hl), respectively. In the case of a boundary node P 2 D in Fig. 4(b), the second equation in Eq. (12) is rewritten as follows:
∂Φ c 1 Rf Φ: _ ¼ 2 1 þ Rf ∂n
(20)
where p(r, t) represents the initial acoustic pressure at r2Ω and time t, cs is the speed of sound, β is the thermal expansion coefficient, CP is the specific heat and δ(t) is the unit impulse function representing the instantaneous laser pulse. In accordance with the relationship among acoustic pressure, particle vibration velocity and medium density, Eq. (20) is rewritten as follows:
(14)
k
1 ∂2 pðr; tÞ β ¼ HðrÞ δ'ðtÞ; c2s ∂t2 Cp
8 ∂p ∂vθ ∂vl βc2 > > þ þ s HðrÞδðtÞ > ¼ ρc2s > > ∂t ∂θ ∂l CP > > > < ∂vθ 1 ∂p ; ¼ > ∂t ρ ∂θ > > > > > ∂vl > 1 ∂p > : ¼ ∂t ρ ∂l
(15)
Thus, the differential expression of Eq. (13) is ½Φði 1; jÞ Φði; jÞ ½Φði; j 1Þ Φði; jÞ jBEj þ kði; j 1Þ jABj hθ hl c 1 Rf Φði; jÞjEAj ½μa ði; jÞc þ qΦði; jÞjDj 2 1 þ Rf
(21)
where ρ is the medium density and vθ and vl are the particle vibration velocity along the θ-axis and l-axis, respectively. The acoustic pressure time series at θm is obtained by solving the FDTD expression of Eq. (21) [32]:
kði 1; jÞ
¼ 0; (16)
8 > > > > > > > 1 ðnþ12Þ 1 1 1 ðnþ12Þ 1 1 βc2 ði; jÞΔt ðnþ12Þ ðnþ12Þ > > pnþ1 ði; jÞ ¼ pn ði; jÞ ρc2s ði; jÞΔt v ; j v ; j þ v v þ s Hði; jÞδn i þ i i; j þ i; j > θ θ l l > Δθ 2 2 Δl 2 2 CP > > > > < 1 1 Δt ðn1Þ ðnþ1Þ ; ½pn ði; jÞ pn ði 1; jÞ vθ 2 i ; j ¼ vθ 2 i ; j > 2 2 ρ Δθ > > > > > > 1 1 Δt ðn1Þ ðnþ1Þ > > ¼ vl 2 i; j ½pn ði; jÞ pn ði; j 1Þ vl 2 i; j > > 2 2 ρΔl > > > > :
where (i, j) is the polar coordinate of the measuring location r; Δθ and Δl are the unit steps along the θ- and l-axis, respectively (here, Δθ ¼ Δl); Δt is the temporal step; n represents the discrete time point; cs(i, j) is the
where jABj, jBEj and jEAj are the lengths of AB, BE and EA, respectively, and jDj is the area of D. For a boundary node P 2 D facing the light source in Fig. 4(c), the third equation in Eq. (12) is rewritten as: k
∂Φ 1 1 Rf c 1 Rf ¼ 4Qs δðr r s Þ Φ: b 2 1 þ Rf 2 1 þ Rf ∂n
(22)
(17)
Thus, the differential expression of Eq. (13) is ½Φði þ 1; jÞ Φði; jÞ ½Φði; j þ 1Þ Φði; jÞ jBEj þ kði; j þ 1Þ jABj hθ hl c 1 Rf Φði; jÞjEAj ½μa ði; jÞc þ qΦði; jÞjDj 2 1 þ Rf 1 1 Rf : ¼ jEAj4Qs 2 1 þ Rf (18)
kði þ 1; jÞ
Finally, the theoretical absorbed energy distribution H(r) is obtained in accordance with HðrÞ ¼ μa ðrÞΦðrÞ:
(19)
Fig. 5. Change curve of the optical radiation energy (Fzr) and absorbed energy (Azr) with respect to the tissue thickness l along θ1 ¼ 0 and θ2 ¼ 45 for vessel phantom a. 41
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
speed of sound at (i, j); pn ði; jÞ represents the acoustic pressure at (i, j) and ðnÞ
2.3. Reconstruction of the optical absorption coefficient
ðnÞ
time n; vθ ði; jÞ and vl ði; jÞ are the vibration velocity along the θ-axis and l-axis at (i, j) and time n and H(i, j) is the absorbed energy density at (i, j).
2.3.1. Formation of the optical inverse problem The problem of recovering the optical absorption coefficient of wall tissues with known scattering and Gruneisen coefficient from PA pressure
Fig. 6. Change curve of the optical radiation energy (Fzr) and absorbed energy (Azr) with respect to the tissue thickness L along θ ¼ 0 –50 for (a) phantom a, (b) phantom b and (c) phantom c.
42
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
kXkTV ¼ jMXj;
measurements is expressed as the minimisation of a nonlinear leastsquare error:
2
b μ a ðrÞ ¼ argmin pðmÞ ðr; tÞ C⋅pðr; t; μa ðrÞ ; μa
where M is a sparse Ne N matrix with grid length as entries, Ne is the iteration time, N is the number of grids and j⋅j is the conventional L1 norm. The simplified TV-regularisation in Eq. (26) can be solved with the L–M algorithm by referring to the subgradient solution of L1regularisation [18,41]:
(23)
where b μ a ðrÞ is the estimated absorption coefficient at r, p(m)(r, t) and p(r, t, μa(r)) represent the measured and modelled acoustic pressure at r and time t and C is a calibration factor related to the acquisition system. The optimisation of Eq. (23) is a nonlinear ill-posed problem. Illposedness can be alleviated through regularisation. Current regularisation methods include L1, L2, Tikhonov and TV [37]. In most research on qPAT, optical parameters are usually assumed to be piecewise constants [38–40]. This assumption is also true for arterial vessels. Hence, in theory, TV-regularisation overperforms the L2-and Tikhonov-regularisation methods [18,41] for qIVPA. For the convenience of expression, we define X ¼ μa(r) and
(26)
∂ kXk1 ¼ signðXÞ; ∂X
(27)
where sign(⋅) is a sign function. 2.3.2. Bregman optimisation The Bregman [18] iteration expression of Eq. (25) is expressed as follows: Xnþ1 ¼ argmin½gðXÞ;
(28)
ðX;dÞ
f ðXÞ ¼ pm ðrÞ C⋅pðr; XÞ : SðXÞ ¼ f ðXÞT f ðXÞ
(24) where 8 α > gðXÞ ¼ SðXn Þ þ kMXn dn νn k22 hpn ; Xn i > > > 2 > <
MXnþ1 vn η ; max jMXnþ1 vn j; dnþ1 ¼ > > α jMXnþ1 vn j > > > : vnþ1 ¼ vn þ dnþ1 MXnþ1
The objective function to be optimised is rewritten as follows: b ¼ argmin SðXÞ þ ηkXkTV ; X
(25)
X
b is the estimate of X and η is the regularisation parameter. At the where X unstructured grids, the TV-regularisation term is simplified as follows [18]:
(29)
where α is the penalty parameter, and vnþ1 is the residual produced in the nth iteration and v0 ¼ 0.
Fig. 7. Distribution map of the optical radiation energy for (a) phantom a, (b) phantom b and (c) phantom c.
Fig. 8. Distribution map of the absorbed energy for (a) phantom a, (b) phantom b and (c) phantom c. 43
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
8 T yn sTn yn sTn sn sTn > 1 1 > > B ¼ I B I þ > nþ1 n T T < yn sn y sn yT sn
The detailed steps to solve Eq. (28) are as follows. Step 1 Set the initial guess of X, i.e., X0, n ¼ 0 and B1 0 ¼ I. Step 2 Update B1 n in accordance with
n
> sn ¼ Xnþ1 Xn > > > : yn ¼ rgnþ1 rgn
n
:
(30)
Fig. 9. Distribution map of the initial acoustic pressure for (a) phantom a, (b) phantom b and (c) phantom c.
Fig. 10. Distribution map of the absorption coefficient reconstructed through Bregman optimisation based on TV-regularisation for (a) phantom a, (b) phantom b and (c) phantom c.
Fig. 11. Change curve of the absorption coefficient and absorbed energy along a radius across the plaque centre (depicted with an arrow in Fig. 10) with respect to the tissue thickness l for (a) phantom a, (b) phantom b and (c) phantom c. ‘rAzr’, ‘Ua’ and ‘rUa’ represent the normalized absorbed energy, predefined absorption coefficient (initial guess), and estimated absorption coefficient through Bregman optimisation. 44
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Then, determine the search direction wn ¼ B1 n rgn where rgn is the gradient of g(X) at Xn
α
rgn ¼ rSn þ M ðMX dn vn Þ pn : 2 T
Step 3 Search for Xnþ1 along wn Xnþ1 ¼ Xn an B1 n rgn ;
(31)
(32)
where an is the step length. Step 4 If krgn k2 < 0:01, then the iteration stops; otherwise, return to Step 2.
Table 4 Measured and actual central thicknesses of the plaques.
3. Results
Vessel phantom
Plaque
Actual thickness (cm)
Measured thickness (cm)
a b c
Lipid pool Calcified plaque Fibro-fatty plaque
0.2267 0.2175 0.0619
0.23 0.21 0.07
The performance of the proposed method was validated with computer-simulated arterial vessel-mimicking phantoms. Several results are presented in this section. Tables 1–3 list the optical parameters and size of the vessel phantoms shown in Fig. 1 under the temperature of 25 C and wavelength of 532 nm in reference to [42]. These parameters
Fig. 12. Polar map of the absorbed energy distribution. (a) Phantom a (DE); (b) Phantom b (DE); (c) Phantom c (DE); (d) Phantom a (MC); (e) Phantom b (MC); (f) Phantom c (MC).
Fig. 13. Change curve of the normalized absorbed energy simulated with MC and DE and the predefined absorption coefficient for (a) phantom a (θ ¼ 38 ), (b) phantom b (θ ¼ 308 ) and (c) phantom c (θ ¼ 18 ). 45
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
of the positive or negative amplitude of PA signals at tissue boundaries during the discrete solution of the wave equation with FDTD. Fig. 10 shows the distribution map of the absorption coefficient recovered through Bregman optimisation based on TV-regularisation. The boundaries of the plaques are illustrated in Fig. 10. In the figure, the boundaries of the mixed plaques are emphasised. However, the plaque boundaries are blurred in Fig. 8, which shows that the shapes and sizes of all the plaques are consistent with the actual geometry of the plaques. Consequently, optical deposition distribution can only present a rough reference for the identification of atherosclerotic lesions. The position, size, shape and optical properties of the plaques provided by the absorption coefficient distribution are more accurate than those provided by the absorbed energy distribution. This conclusion is supported by Fig. 11, where the distance between the adjacent troughs of the curves determines the size of a plaque. To quantitatively evaluate the accuracy of absorption coefficient recovery, the central thicknesses of the plaques shown in Fig. 10 were measured and compared with the real values listed in Tables 1–3. The results shown in Table 4 imply that the measured thicknesses of phantom a and b are consistent with the actual thicknesses of phantoms a and b. For phantom c, the measured thickness of the fibro–fatty plaque near the wall is consistent with the true value. However, the thickness of the calcified plaque distant from the lumen is larger than the actual value.
Fig. 14. Change curve of the absorption coefficient with respect to the tissue thickness l along θ1 ¼ 0 (normal vessel wall tissue) and θ2 ¼ 45 (including plaque tissue) for phantom a. ‘Ua’ and ‘rUa’ represent the predefined (initial guess) and reconstructed absorption coefficient. (a) θ1 ¼ 0 , Bregman based on L2; (b) θ1 ¼ 0 , Bregman based on TV; (c) θ2 ¼ 45 , Bregman based on L2; (d) θ2 ¼ 45 , Bregman based on TV.
were used in the forward simulation experiments and were also the initial guesses in the optimisation of Eq. (23). Figs. 5–8 show the distribution of the optical radiation energy and absorbed energy obtained through the forward simulation. The results shown in Fig. 5 are consistent with the geometry and tissue components shown in Fig. 1(a), where the lumen, intima/media and adventitia are found along θ1 ¼ 0 and the lumen, TCFA, lipid pool, intima/media and adventitia are found along θ1 ¼ 45 . The results shown in Figs. 5–8 indicate that the optical radiation energy is related to the absorption and scattering coefficients. However, the absorbed energy is dependent on the absorption coefficient. The radiation and absorbed energy attenuate to 0 with increasing detection depth. The optical deposition in the lumen and lipid and fibro–fatty plaques is higher than that in other tissues given the larger absorption coefficients of the lumen and plaques. Moreover, the photons absorbed by the lumen may be scattered to other tissues because the lumen is filled with blood, which has a large scattering coefficient. The intima/media absorbs partial light energy, and partial photons scatter to other tissues. The optical deposition in the adventitia is nearly 0. Fig. 9 shows the distribution map of the initial acoustic pressure after pseudocolour coding. The ring-down artefacts surrounding the catheter and located between the lumen and wall are caused by the normalisation
4. Discussion 4.1. Comparison of DE and MC The absorbed energy distributions simulated with DE and MC, a random statistical method for solving the RTE [43,44] were compared. The results are shown in Figs. 12 and 13. The figures indicate that the optical deposition obtained with DE distinguishes the range of plaque tissues more accurately than that obtained with MC. 4.2. Comparison of L2-and TV-regularisation The absorption coefficients recovered via Bregman optimisation based on L2- [18] and TV-regularisation were compared, and the results are shown in Figs. 10, 14 and 15. The curves in Fig. 14 exhibit noticeable peaks at the boundary of the lumen–intima/media, intima/media–plaque and plaque–adventitia. The calcified plaque can be identified and localised through TV- and L2-regularisation reconstruction. The quantification results listed in Table 5 suggest that the iteration times of TV-regularisation are lower than those of L2-regularisation although the estimated absorption coefficients of the two methods are similar.
Fig. 15. Distribution map of the absorption coefficient recovered via Bregman optimisation based on L2-regularisation for (a) phantom a, (b) phantom b and (c) phantom c. 46
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
or L2-regularisation are provided in this section. The L–M iteration expression of Eq. (25) is
4.3. Comparison of L–M and Bregman optimisation The results for the comparison of the absorption coefficient distribution recovered via L–M [28] and Bregman optimisation based on TV-
Table 5 Absorption coefficients recovered via L–M or Bregman optimisation based on L2-or TV-regularisation for phantom a. No.
Measuring location (θ,l)
Predefined
μa
L-M
2 3 4 5 6 7 8 9 10
(45 ,0.05) (0 , 0.05) (45 ,0.10) (0 , 0.10) (45 ,0.115) (0 ,0.115) (45 ,0.125) (0 ,0.125) (45 ,0.15) (0 ,0.15) (45 ,0.20) (0 ,0.20) (45 ,0.25) (0 ,0.25) (45 ,0.367) (0 ,0.367) (45 ,0.45) (0 ,0.45) (45 ,0.50) (0 ,0.50)
0.99 0.99 0.99 0.99 0.99 0.99 0.10 0.20 0.60 0.20 0.60 0.20 0.60 0.20 0.60 0.20 0.20 0.20 0.20 0.20
TV
L2
Estimated
μa
Iteration times
0.72 0.73 0.98 0.98 0.92 0.92 0.12 0.43 0.45 0.24 0.58 0.21 0.59 0.17 0.48 0.18 0.23 0.18 0.21 0.19
28 25 26 25 25 24 26 25 25 23 26 24 23 24 28 26 25 23 26 25
Estimated 1
Bregman
L2
μa
Iteration times
0.74 0.75 0.99 0.99 0.92 0.93 0.11 0.40 0.45 0.21 0.55 0.20 0.60 0.19 0.49 0.19 0.21 0.19 0.20 0.20
13 12 15 10 14 11 12 13 12 12 13 10 12 11 11 12 10 10 11 10
Estimated
TV
μa
Iteration times
Estimated
μa
Iteration times
0.74 0.75 0.98 0.98 0.93 0.92 0.11 0.42 0.47 0.22 0.56 0.21 0.60 0.18 0.49 0.19 0.21 0.19 0.20 0.19
27 24 28 26 24 25 25 24 25 24 27 25 25 24 28 26 26 24 25 24
0.75 0.76 0.99 0.99 0.93 0.93 0.10 0.40 0.48 0.20 0.56 0.20 0.60 0.19 0.50 0.20 0.20 0.20 0.20 0.20
14 12 15 10 13 13 11 12 12 11 12 10 13 10 10 11 11 10 10 10
Fig. 16. Distribution map of the absorption coefficient recovered via L–M optimisation. (a) Phantom a (L–M based on L2); (b) Phantom b (L–M based on L2); (c) Phantom c (L–M based on L2); (d) Phantom a (L–M based on TV); (e) Phantom b (L–M based on TV); (f) Phantom c (L–M based on TV). 47
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Fig. 17. Change curve of the absorption coefficient along a radius across the plaque centre (depicted with an arrow in Fig. 10) with respect to the tissue thickness l for (a) phantom a, (b) phantom b and (c) phantom c. ‘Ua’, ‘LM’ and ‘Bregman’ represent the predefined absorption coefficient (initial guess), recovered absorption coefficient via L–M and Bregman optimisation, respectively.
8 2 > < ΔXnþ1 ¼ argmin krSðXn Þ⋅ΔX ½SðXÞ SðXn Þ þ ωn k2 þ ηkΔXkTV ΔX
ω ¼ ωn þ ½SðXÞ SðXn ÞSðXÞ SðXn Þ rSðXn Þ⋅ΔXnþ1 > : nþ1 Xnþ1 ¼ Xn þ ΔXnþ1
and to discriminate between normal arterial tissues and lipid-rich plaques on the basis of the strong absorption of lipids has attracted increased attention [45,46]. At 1.2 μm, the IVPA system requires the laser energy of at least 1.2 mJ to distinguish between lipids and other tissues. However, it only requires the laser energy of 0.4 mJ at 1.7 μm. Therefore, IVPA images with sufficient contrast can be produced with less energy at 1.7 μm than at other NIR wavelengths [47]. Future works will demonstrate the validity of the proposed method at 1.7 μm by adjusting the optical parameter setting in the forward simulation. Moreover, a simple photon propagation pattern was adopted in the forward model. The wall tissues were assumed to be uniformly illuminated in a 2-D space parallel to the imaging plane by ignoring the axial influence of the narrow beam laser pulse and nonuniform distribution of light in the medium. The possible changes in the Gruneisen coefficient and acoustic inhomogeneity were ignored. The ultrasonic detector was simplified as an ideal point-like detector by ignoring the aperture effect related to its size. The simplification introduced differences between the theoretical and measured PA pressure and imposed possible limitations on the practical applications of the proposed method.
; (33)
where ΔX is the increment of X. The detailed steps to solve Eq. (33) are as follows. Step 1 Set the initial guess of X, i.e., X0, n ¼ 0, ω02(0,1), τ > 1, v > 1 and ε > 0. Step 2 Determine f(Xn), S(Xn), rf ðXn Þ and rSðXn Þ ¼ ðrf ðXn ÞÞT ⋅f ðXn Þ. Step 3 Solve ½Q þ τIΔX ¼ rSðXn Þ to obtain ΔX, where Q ¼ ðrf ðXn ÞÞT ⋅rf ðXn Þ and I is the unit matrix. Step 4 Calculate Xnþ1 ¼ Xn þ ΔX. If the absolute error is less than ε, the iteration ends; otherwise, return to Step 5. Step 5 If SðXnþ1 Þ < SðXn Þ þ ωn ½rSðXn ÞT ΔX, then τ←τ/v and proceed to Step 6; otherwise, τ←τ⋅v and return to Step 3. Step 6 k←kþ1 and return to Step 2. The results are shown in Figs. 10, 15 and 16. The visualisation of the images reconstructed with TV- and L2-regularisation combined with the same optimisation method is similar. However, by referring to the actual geometry of the phantoms, the Bregman reconstruction provides visually better images than the L–M reconstruction based on the same regularisation. The absorption coefficients reconstructed with L–M and Bregman optimisation were compared with the predefined values listed in Tables 1–3. Fig. 17 indicates the similarity shared by L–M and Bregman optimisation. However, the iteration times of Bregman optimisation are higher than those of L–M optimisation. Considering the precision and computational cost, L–M may be more suitable than Bregman optimisation for the individual recovery of absorption coefficients, whereas Bregman optimisation may be a better choice than L–M for the simultaneous estimation of absorption and scattering coefficients.
5. Conclusion In this preliminary study, a one-step qIVPA method for the singlewavelength IVPA system was designed. The method successfully recovered the transversal maps of optical absorption coefficient distribution on vascular cross-sections from PA boundary measurements. To implement fully quantitative optical inversion, the model-based inversion utilised the DA to the RTE to model light transport in a medium with a known scattering coefficient and accounted for acoustic propagation in acoustically homogeneous lossless medium. The L–M or Bregman optimisation based on TV-regularisation was applied to minimise a nonlinear leastsquare error function to quantify the difference between the measured and modelled acoustic pressure series. The experimental results obtained for computer-simulated arterial phantoms indicate that the lumen, TCFA, lipid plaque, calcified plaque, fibro–fatty plaque and vessel-wall intima/ media/adventitia can be clearly identified in an absorption coefficient distribution map. Moreover, the location, size and contour of the TCFA, lipid plaque, calcified plaque and fibro–fatty plaque can be clearly distinguished.
4.4. Other discussions The theoretical acoustic pressure generated by vessel-wall tissues were obtained through the forward simulation of excitation laser transport and tissue absorption. A narrow beam of collimated excitation light was assumed to be normally incident upon the wall tissues through the blood. The influence of luminal blood on light transport was ignored. In practical applications, haemoglobin absorption is high in the visible range (410–680 nm). Nevertheless, IVPA can discriminate between normal and atheromatous tissues at these wavelengths. In recent years, the use of experimental IVPA systems at near-infrared excitation wavelengths, such as 1.2 or 1.7 μm, to minimise the effect of blood absorption
Conflicts of interest None Declared. Acknowledgments This work was supported by National Nature Science Foundations of 48
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
China (no. 61372042).
[23] O. Nyk€anen, A. Pulkkinen, T. Tarvainen, Quantitative photoacoustic tomography augmented with surface light measurements, Biomed. Optic Express 8 (10) (2017) 4380–4395. [24] A.Q. Bauer, R.E. Nothdurft, T.N. Erpelding, et al., Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography, J. Biomed. Optic. 16 (9) (2011), 096016. [25] J.R. Rajian, P.L. Carson, X. Wang, Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent, Optic Express 17 (6) (2009) 4879–4889. [26] K. Daoudi, A. Hussain, E. Hondebrink, et al., Correcting photoacoustic signals for fluence variations using acousto-optic modulation, Optic Express 20 (13) (2012) 14117–14129. [27] Y. Zhen, J. Huabei, A calibration-free, one-step method for quantitative photoacoustic tomography, Med. Phys. 39 (11) (2012) 6895–6899. [28] T. Ding, K. Ren, S. Vallelian, A one-step reconstruction algorithm for quantitative photoacoustic imaging, Inverse Probl. 31 (9) (2015) 65003–65023. [29] M. Haltmeier, L. Neumann, S. Rabanser, Single-stage reconstruction algorithm for quantitative photoacoustic tomography, Inverse Probl. 31 (6) (2015) 65005–65028. [30] G. Bal, A. Moradifam, Photo-acoustic tomography in a rotating measurement setting, Inverse Probl. 32 (10) (2016), 105012. [31] A. Pulkkinen, B.T. Cox, S.R. Arridge, et al., Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography, IEEE Trans. Med. Imag. 35 (11) (2016) 2497–2508. [32] Z. Sun, Y. Yuan, D. Han, A computer-based simulator for intravascular photoacoustic images, Comput. Biol. Med. 81 (2017) 176–187. [33] T. Tarvainen, Computational Methods for Light Transport in Optical Tomography, University of Kuopio, Kuopio, 2006. [34] S. Li, B. Montcel, Z. Yuan, et al., Multigrid-based reconstruction algorithm for quantitative photoacoustic tomography, Biomed. Optic Express 6 (7) (2015) 2424–2434. [35] F. Gao, Y. Tanikawa, H. Zhao, et al., Semi-three-dimensional algorithm for timeresolved diffuse optical tomography by use of the generalized pulse spectrum technique, Appl. Optic. 41 (34) (2002) 7346–7358. [36] L.V. Wang, Biomedical Optics: Principles and Imaging. Wang Lihong. Photoacoustic Tomography, John Wiley and Sons, New Jersey, 2007, pp. 283–321. [37] H. Gao, H. Zhao, Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization, Optic Express 18 (3) (2010) 1854–1871. [38] W. Neater, O. Scherzer, Quantitative photoacoustic tomography with piecewise constant material parameters, SIAM J. Imag. Sci. 7 (3) (2014) 1755–1774. [39] E. Beretta, M. Muszkieta, W. Naetar, et al., A Variational Method for Quantitative Photoacoustic Tomography with Piecewise Constant Coefficients, arXiv: 1509.04982, 2015. [40] A. De Cezaro, A. Leit~ao, X. Tai, On piecewise constant level-set (PCLS) methods for the identification of discontinuous parameters in ill-posed problems, Inverse Probl. 29 (29) (2012), 015003. [41] H. Gao, H. Zhao, Multilevel bioluminescence tomography based on radiative transfer equation part 2: total variation and l1 data fidelity, Optic Express 18 (3) (2010) 2894–2912. [42] S.L. Jacques, Optical properties of biological tissues: a review, Phys. Med. Biol. 58 (11) (2013) R37–R61. [43] W. Li, L. Lin, L. Bao, et al., Monte Carlo simulation of photon migration in multicomponent media, Opt. Quant. Electron. 47 (7) (2015) 1919–1931. [44] L. Wang, S.L. Jacques, L. Zheng, MCML-Monte Carlo modeling of light transport in multi-layered tissues, Comput. Meth. Progr. Biomed. 47 (2) (1995) 131–146. [45] T.J. Allen, P.C. Beard, Photoacoustic characterisation of vascular tissue at NIR wavelengths, in: A.A. Oraevsky, L.V. Wang (Eds.), Proceedings of SPIE International Conference on Photons Plus Ultrasound: Imaging and Sensing 2009, vol. 7177, 25 Feb. 2009 id: 71770A. [46] N. Dana, T. Sowers, A. Karpiouk, et al., Optimization of dual-wavelength intravascular photoacoustic imaging of atherosclerotic plaques using Monte Carlo optical modeling, J. Biomed. Optic. 22 (10) (2017) 1–12. [47] K. Jansen, M. Wu, A.F.W. van der Steen, et al., Photoacoustic imaging of human coronary atherosclerosis in two spectral bands, Photoacoustics 2 (1) (2014) 12–20.
References [1] W. Min, F.W.S. Antonius, R. Evelyn, et al., Emerging technology update intravascular photoacoustic imaging of vulnerable atherosclerotic plaque, Intervent Cardiol. 11 (2) (2016) 120–123. [2] J. Yao, L.V. Wang, Breakthroughs in photonics photoacoustic tomography, IEEE Photon. J. 6 (2) (2014), 0701006. [3] S.Y. Nam, E. Chung, L.J. Suggs, et al., Combined ultrasound and photoacoustic imaging to noninvasively assess burn injury and selectively monitor a regenerative tissue-engineered construct, Tissue Eng. C Meth. 21 (6) (2015) 557–566. [4] B. Wang, S. Emelianov, Thermal intravascular photoacoustic imaging, Opt. Soc. Am. 2 (11) (2011) 3072–3078. [5] D. Verya, W. Min, J. Nico de, et al., Frequency analysis of the photoacoustic signal generated by coronary atherosclerotic plaque, Ultrasound Med. Biol. 42 (8) (2016) 2017–2025. [6] B. Wang, E. Yantsen, T. Larson, et al., Plasmonic intravascular photoacoustic imaging for detection of macrophages in atherosclerotic plaques, Nano Lett. 9 (6) (2009) 2212–2217. [7] W. Min, S. Geert, L. Matija, et al., Real-time volumetric lipid imaging in vivo by intravascular photoacoustics at 20 frames per second, Biomed. Optic Express 8 (2) (2017) 943–953. [8] X. Bai, X. Gong, W. Hau, et al., Intravascular optical-resolution photoacoustic tomography with a 1.1 mm diameter catheter, PLoS One 9 (3) (2014) e92463. [9] Y. Liang, L. Jin, L. Wang, et al., Fiber-laser-based ultrasound sensor for photoacoustic imaging, Sci. Rep. 7 (2017) 40849. [10] J. Zhang, S. Yang, X. Ji, et al., Characterization of lipid-rich aortic plaques by intravascular photoacoustic tomography: ex vivo and in vivo validation in a rabbit atherosclerosis model with histologic correlation, J. Am. Coll. Cardiol. 64 (4) (2014) 385–390. [11] H. Qin, T. Zhou, S. Yang, et al., Gadolinium (III)-gold nanorods for MRI and photoacoustic imaging dual-modality detection of macrophages in atherosclerotic inflammation, Nanomedicine 8 (10) (2013) 611–1624. [12] Y.L. Sheu, C.Y. Chou, B.Y. Hsieh, et al., Application of limited-view image reconstruction method to intravascular photoacoustic tomography, in: A.A. Oraevsky, L.V. Wang (Eds.), Proceedings of SPIE International Conference on Photons Plus Ultrasound: Imaging and Sensing 2010, vol. 7564, 24 Feb. 2010 id: 75640B. [13] Y.L. Sheu, C.Y. Chou, B.Y. Hsieh, Image reconstruction in intravascular photoacoustic imaging, IEEE Trans. Ultrason. Ferroelectrics Freq. Contr. 58 (2011) 2067–2077. [14] Z. Sun, D. Han, J. Wang, Review of image reconstruction for intravascular photoacoustic imaging, Opto-Electronic Eng. 42 (3) (2015) 20–27. [15] Z. Sun, D. Han, Y. Yuan, 2-D image reconstruction of photoacoustic endoscopic imaging based on time-reversal, Comput. Biol. Med. 76 (2016) 60–68. [16] N. Song, C. Deumic, S.A. Da, Considering sources and detectors distributions for quantitative photoacoustic tomography, Biomed. Optic Express 5 (11) (2014) 3960–3974. [17] H. Gao, S. Osher, H. Zhao, Mathematical Modeling in Biomedical Imaging II, Springer, Berlin/Heidelberg, 2012, pp. 131–158. [18] H. Gao, H. Zhao, S. Osher, Bregman methods in quantitative photoacoustic tomography, CAM Rep. 30 (6) (2010) 3043–3054. [19] B. Cox, J.G. Laufer, S.R. Arridge, et al., Quantitative spectroscopic photoacoustic imaging: a review, J. Biomed. Optic. 17 (6) (2012), 061202. [20] X. Zhang, W. Zhou, X. Zhang, et al., Forward-backward splitting method for quantitative photoacoustic tomography, Inverse Probl. 30 (2014), 125012. [21] R. Hochuli, S. Powell, S. Arridge, et al., Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance, J. Biomed. Optic. 21 (2016), 126004. [22] M. Venugopal, E.P. Van, S. Manohar, et al., Quantitative photoacoustic tomography by stochastic search: direct recovery of the optical absorption field, Optic Lett. 41 (18) (2016) 4202–4205.
49
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/compbiomed
Reconstruction of optical absorption coefficient distribution in intravascular photoacoustic imaging Sun Zheng *, Zheng Lan Department of Electronic and Communication Engineering, North China Electric Power University, Baoding 071003, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Intravascular photoacoustics (IVPA) Inverse problem Quantitative intravascular photoacoustics (qIVPA) Optical absorption coefficient Optical inversion
Intravascular photoacoustic (IVPA) imaging can discriminate between normal arterial tissues and atheromatous plaques in images, particularly images of lipid-rich plaques with high spatial resolution and optical contrast. However, conventional IVPA only recovers the deposited optical energy that is the product of the tissue optical absorption coefficient and local optical fluence. Herein, a one-step model-based method for single-wavelength IVPA imaging system is proposed. The proposed method directly reconstructs the optical absorption coefficient from the boundary measurement of acoustic pressure. To obtain the theoretical optical deposition, light illumination and transport in multilayered vessel-wall tissues are numerically modelled with the diffusion approximation to radiative transfer equation. Then, the generation and propagation of photoacoustic (PA) waves in acoustically homogeneous vessel-wall tissues are modelled by using the PA wave equation. The theoretical acoustic pressure series is thus obtained. Finally, the optical absorption coefficient is iteratively updated from an initial guess by minimising a nonlinear least-square error function to quantify the difference between the measured and theoretical acoustic pressure with split Bregman optimisation based on total-variation regularisation. The numerical simulation experiments demonstrated that optical absorption coefficient maps can be directly recovered on the basis of the acoustic boundary measurement of vascular cross-sections. Compared with the optical deposition map, the optical absorption coefficient map can provide more reliable qualitative and quantitative information on the morphology and optical properties of the imaged arteries. The proposed method enables the accurate and reliable evaluation of atheromatous lesions and the early identification of vulnerable plaques.
1. Introduction Intravascular photoacoustic (IVPA) imaging is a typical application of endoscopic photoacoustic tomography (PAT). It has recently emerged as a complementary imaging modality to intravascular ultrasound (IVUS) for the accurate evaluation of arterial atherosclerotic lesions and the effective identification of vulnerable plaques [1]. This approach can visualise the structure and functional composition of arterial vessel walls and atheromatous plaques on the basis of the characteristics of tissue optical absorption. During IVPA intervention, a catheter integrated with a multimode optical fibre, an optical reflection element and an ultrasonic detector is inserted and pushed to the distal end of the target vascular lumen. During the drawback of the catheter, a probe mounted on the catheter tip emits modulated light that is absorbed by the surrounding tissues and then converted to heat. Periodic heat flow causes the thermal expansion of
tissues, thereby exciting a wide band of ultrasonic waves, i.e., photoacoustic (PA) waves. The PA waves are transmitted rapidly to the surface and are finally collected by the ultrasonic detector. The distribution of the absorbed optical energy density can be reconstructed from the timedependent boundary measurement of acoustic pressure through acoustic inversion. To reveal the functional components of normal and pathological tissues, the optical parameters of the tissues are recovered by solving the optical inverse problem [2]. IVPA exhibits the high spatial resolution of IVUS for imaging deep tissues [3] and the main advantages of other optical imaging techniques, such as limited tissue damage, excellent optical contrast and functional and molecular imaging. Important progress has been made in improving thermal IVPA [4], spectroscopic IVPA [5], plasmonic IVPA [6] and high-speed IVPA [7] systems. Moreover, integrated IVPA/IVUS imaging catheter [8] and novel ultrasonic detectors [9] have been developed and utilised in the characterisation of lipid-rich plaques [10] and the detection of
* Corresponding author. P.O. Box 21, North China Electric Power University, Baoding 071003, China. E-mail address: [email protected] (S. Zheng). https://doi.org/10.1016/j.compbiomed.2018.04.012 Received 27 December 2017; Received in revised form 16 April 2018; Accepted 16 April 2018 0010-4825/© 2018 Elsevier Ltd. All rights reserved.
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
by Ref. [30]. In this setting, the pulsed laser source and an ultrasonic transducer array are mounted on a same rotating frame, and the distribution map of the optical absorption is directly recovered from the measured acoustic pressure series. A forward model was constructed by Ref. [31] to combine an optical transport model with an acoustic propagation model. This model simultaneously and directly estimates absorption and scattering coefficients from the PA pressure time series through a Bayesian approach. The scanning aperture of IVPA is enclosed in a geometrically complex lumen. Well-developed qPAT methods cannot be directly applied in qIVPA. This study aims to design a one-step qIVPA method for an IVPA system with single optical excitation at a single wavelength. The transversal map of the absorption coefficient distribution on a vascular crosssection was directly recovered from the acoustic pressure time series generated by multilayered arterial wall tissues and atherosclerotic plaques and then measured using the ultrasonic detector. A forward model based on diffusion approximation (DA) to the radiative transfer equation (RTE) was constructed to simulate theoretical optical deposition in acoustically homogenous tissues with known scattering and Gruneisen coefficients. The PA pressure was modelled by explicitly solving the PA wave equation with a finite difference time domain (FDTD) algorithm. The absorption coefficient of a vascular cross-section was iteratively updated from an initial guess by minimising a squared error function that quantifies the difference between the measured and modelled acoustic pressure series with split Bregman optimisation based on total variation (TV) regularisation. The performance of the proposed method was demonstrated with computer-simulated arterial vessel-mimicking phantoms. The ability of the Monte Carlo (MC) simulation to model light transport was compared with that of DA in accordance with the validation results. The accuracy of L–M and split Bregman optimisation based on TV- and L2-regularisation was also evaluated and compared.
macrophages in atherosclerotic inflammation [11]. However, the integration, mobility, operability, imaging speed and reliability of the IVPA imaging system still cannot meet the standards of clinical applications. The acoustic inverse problem in IVPA involves the reconstruction of the initial acoustic pressure or the deposited optical energy from the time-dependent boundary measurement of acoustic pressure by using filtered back projection, time reversal and iterative algorithms [12–15]. Optical inversion, the so-called quantitative IVPA (qIVPA), is essentially the inversion of a light transport operator and involves the recovery of optical property parameters, particularly the absorption and scattering coefficient [16–18] and thermal expansion parameter (Gruneisen coefficient) [19], from the time-dependent measurement of acoustic pressure or optical deposition recovered through acoustic inversion. The acoustic and optical inverse problems can be solved separately because the time of acoustic propagation in biological soft tissues is approximately three orders of magnitude longer than that of laser illumination and optical energy absorption by tissues [19]. The optical absorption is determined by the local absorption coefficient and light fluence and cannot directly reflect the optical properties of the tissues, which are closely related to chemical composition. The optical properties of normal tissues drastically differ from those of diseased tissues. Therefore, qIVPA can provide more accurate and reliable evaluation for the early identification of vulnerable atherosclerotic plaques than ordinary IVPA. However, to the best of our knowledge, no special methods have been proposed for qIVPA thus far. Current quantitative PAT (qPAT) methods often involve two steps [20–23]. Firstly, the distribution map of initial acoustic pressure is reconstructed from acoustic boundary measurements through acoustic inversion. Optical deposition measurements are determined on the basis of PA images. Secondly, a forward model of light transport in the tissues is established to obtain the theoretical optical deposition. Finally, the map of the optical absorption and scattering coefficient is recovered by minimising the squared error between the measured and theoretical optical deposition. Multiple-illumination PAT can alleviate the ill-posedness of the optical inversion caused by absorption-scattering nonuniqueness. Furthermore, the distribution map of the speed of sound in an acoustically inhomogeneous medium can be simultaneously recovered [20]. The main limitation of two-step qPAT methods is that the accuracy of the optical map highly depends on the optical deposition reconstructed from the boundary PA measurements. Errors are inevitably generated during acoustic inversion given the limited bandwidth of the ultrasonic detector in contrast to the ultrawide PA frequency response; the orientation, size and shape of the imaging target with respect to the detecting aperture; the effect of optical scattering and the acoustic inhomogeneity and complexity of the medium. In addition, the qPAT two-step methods are limited by the cost of imaging systems based on an optical method, e.g., diffuse optical tomography [24], optical contrast agent [25] or acousto–optic modulation [26] to measure the optical fluence [27]. Furthermore, several artefacts in the optical deposition map may reduce the accuracy of optical parameters. The ability of direct one-step optical map recovery from PA boundary measurements to resolve the above issues has been studied recently [27–31]. For example, a calibration-free qPAT method proposed by Ref. [27] needs only the normalized acoustic boundary measurements and the laser source intensity as the input data. These requirements eliminate the influence of detector sensitivity. The method reported by Ref. [28] allows the simultaneous reconstruction of the absorption coefficient and speed of sound from multiple acoustic pressure datasets generated from multiple optical illuminations. In regions where the changes in the absorption coefficient and speed of sound are relatively small, the perturbation to the absorption coefficient around the known backgrounds is calculated with Born approximation and Landweber iteration. The fully nonlinear inverse problem is solved with an optimal control method, where the objective function is minimised with the Levenberg–Marquardt (L–M) algorithm to obtain the optimal solution of the absorption coefficient. A rotating measurement setting was designed
2. Method 2.1. Construction of the imaging model Geometrical models of arterial cross-sections containing atherosclerotic plaques are constructed to illustrate the proposed method. The models are shown in Fig. 1. The cross-section is divided into M uniform sectors, and the mth measuring angle of the detector is θm ¼ 360(i1)/M (see Fig. 2). The angles of the corresponding imaging region are in the range of [θma, θmb], where m ¼ 1,2, …,M, θma ¼ θm 180/M and θmb ¼ θm þ 180/M. A θ‒l polar coordinate system centred at the lumen is built. In this system, θ and l represent the polar angle and polar radius, respectively. As shown in Fig. 3, the ultrasonic detector mounted on the catheter tip is located at the lumen centre and collects PA waves generated by the surrounding tissues along a circle parallel to the imaging plane in the closed lumen [32]. The lumen, vessel-wall intima/media, plaques and adventitia surround the catheter centre. The detector is simplified to be an ideal point-like detector, and the aperture effect related to its size is ignored. 2.2. Simulation of modelled PA pressure 2.2.1. Formation of the forward problem A single short-pulsed laser at a single wavelength illuminates the vessel wall from the lumen centre. The optical energy absorbed by the surrounding tissues is expressed as follows: HðrÞ ¼ μa ðrÞΦðr; μa ðrÞ; μs ðrÞÞ;
(1)
where r 2 Ω represents a point in a bounded imaging domain Ω with a smooth boundary ∂Ω; H(r) is the absorbed energy density at r; μa(r) and μs(r) are the wavelength-dependent optical absorption coefficient and scattering coefficient at r satisfying the Lipschitz continuity and Φ is the 38
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Fig. 1. Cross-sectional geometry of arterial vessels including (a) thin-cap fibroatheroma (TCFA) and lipid plaques, (b) TCFA and calcified plaques, and (c) TCFA, fibrofatty and calcified plaques.
μ0s ðrÞ≫μa ðrÞ;
(3)
the RTE is simplified as a diffusion equation (DE) represented with spherical harmonics where
μs 'ðrÞ ¼ μs ðrÞð1 gÞ;
(4)
is the reduced scattering coefficient at r, and g is the anisotropy factor. The time-independent DE satisfying the Robin boundary condition is represented as follows [17,34]: r⋅kðrÞrΦðr; tÞ þ μa ðrÞΦðr; tÞ ¼ Qðr; tÞ;
(5)
where Φ(r, t) is the light fluence at r and time t, Q(r, t) is the light source and k(r) is the optical diffusion coefficient. k(r) is expressed as follows:
Fig. 2. Multi-layered vessel wall tissue in a θ‒l polar coordinate system.
kðrÞ ¼
1 1 : 3½μa ðrÞ þ μs 'ðrÞ 3μ0s ðrÞ
(6)
The time-domain Robin boundary condition is written as follows [17]: cΦðr; tÞ þ 2kðrÞ
1 þ Rf _ n ⋅rΦðr; tÞ ¼ 0; 1 Rf
(7)
_
where r 2 ∂Ω, n is the outer normal vector of the tissue surface at ∂Ω, c is the speed of light in tissues and Rf is the internal reflection coefficient of diffuse transmission. Rf is approximated with a polynomial in accordance with Fresnel's law of reflection, which is expressed as follows: Rf 1:4399n2 þ 0:7099n1 þ 0:6681 þ 0:0636n;
where n is the refractive index of the target tissue relative to the background. A narrow beam of collimated light is normally incident upon the tissue surface and is approximated as an internal photon flux that propagates at the tissue boundary. The directivity of the light basically disappears when the penetration depth is approximately one photon free path in length because of the strong scattering effect of tissues. Light propagation becomes isotropic, and the source is approximately an undirected point at the penetration depth. The point-like source is represented with a δ-function with respect to time and space:
Fig. 3. Schematics of IVPA imaging.
light fluence (i.e., optical radiance energy) related to r, μa(r) and μs(r). The initial acoustic pressure at r is determined as follows: b pðrÞ ¼ ΓðrÞHðrÞ;
(8)
(2)
b is the conversion efficiency of heat energy to pressure and where Γ represents the thermodynamic properties of the medium. The acoustic pressure time series p(r, t) is finally measured by the detector during rotation inside the lumen.
Qðr; tÞ ¼ δðr r s ; tÞ;
(9)
where rs is the location of the light source in the θ‒l polar coordinate system. The collimated light source model is included in the Robin boundary condition in Eq. (7):
2.2.2. Forward model of light transport Light transport in turbid media, such as biological soft tissues, is usually modelled by the RTE, an integro–differential equation describing photon distribution in the phase space [33]. RTE can accurately describe photon migration in tissues, especially in nondiffusion regions. The DA is usually used to simplify the solution to the RTE because the direct solution may be highly computationally complex. Assuming
cΦðr; tÞ þ 2kðrÞ
1 þ Rf _ n ⋅rΦðr; tÞ ¼ 4Qs wðrÞδðtÞ; 1 Rf
8r 2 ∂Ω;
(10)
where Qs is the intensity of the light source, and w(r) is the weighting 39
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Fig. 4. Meshing of the imaging region. (a) Interior point; (b) Boundary point; (c) Boundary point facing to the light source.
location function, wðrÞ ¼ δðr r s Þ:
where q is the complex frequency factor of LT, Φ(r,rs,q) is the LT of the time-domain optical radiance energy at r and time t and ∂Ω0 is the tissue surface facing the source. If q ¼ 0, DE is a steady-state diffusion equation.
(11)
2.2.3. FDTD expression of DE The generalised pulse spectrum technique [35] is applied to simplify the solution of DE. The Laplace transform (LT) of a time-domain curve at a complex frequency is used instead of the whole time-domain curve. The forward problem is then transformed into a two-dimensional (2-D) problem that is time-related but is not time-dependent. As mentioned above, the collimated light source model is integrated into the Robin boundary condition. Thus, the source term of the original DE is 0. The DE, Robin boundary condition and Robin boundary condition including the source item are expressed in complex frequency domain as: 8 r⋅kðrÞrΦðr; r s ; qÞ ½c⋅μa ðrÞ þ qΦðr; r s ; qÞ ¼ 0 > > > > > > < cΦðr; r ; qÞ ¼ 2kðrÞ 1 þ Rf _ n ⋅rΦðr; r s ; qÞ s 1 Rf > > > > 1 þ Rf _ > > n ⋅rΦðr; r s ; qÞ ¼ 4Qs δðr r s Þ : cΦðr; r s ; qÞ þ 2kðrÞ 1 Rf
2.2.4. FDTD solution of optical deposition As shown in Fig. 4, a smooth region D 2 Ω is meshed with unstructured rectangular grids represented by nodes in the θ‒l coordinate system, where hθ and hl are the steps along the positive θ- and l-axis, respectively, and (ihθ, jhl) is the coordinate of a point P 2 D. A relatively narrow beam of collimated light propagates along the positive l-axis. Both ends of the first equation in Eq. (12) are integrated in D to yield the following equation: ∫k
∂D
r2Ω r 2 ∂Ω
∂Φ _ ds ∬ ðμa c þ qÞΦdθdl ¼ 0; ∂n D
8D 2 Ω:
(13)
where s2∂Ω. Then, the optical radiance energy at each node is obtained by solving the FDTD expression of DE. In the case of an interior node P 2 D in Fig. 4(a), the differential expression of Eq. (13) is expressed as follows:
;
r 2 ∂ Ω0 (12)
Table 1 Parameters of vessel phantom a. Tissue name
Tissue composition
Average refractive index
Absorption coefficient (cm1)
Scattering coefficient (cm1)
Anisotropy factor
Thickness (cm)
Adventitia Intima/Media Lipid pool TCFA Lumen
Connective tissues Muscular tissues Lipid Fibrous tissues Blood
1.41 1.41 1.46 1.46 1.34
0.20 0.20 0.60 0.01 0.99
20 20 500 80 800
0.80 0.80 0.80 0.75 0.99
0.1330 0.2520 0.2267 (middle) 0.0100 0.1150
Table 2 Parameters of vessel phantom b. Tissue name
Tissue composition
Average refractive index
Absorption coefficient (cm1)
Scattering coefficient (cm1)
Anisotropy factor
Thickness (cm)
Adventitia Intima/Media Calcified plaque TCFA Lumen
Connective tissue Muscular tissue Calcium Fibrous tissue Blood
1.41 1.41 1.46 1.46 1.34
0.20 0.20 0.05 0.01 0.99
20 20 80 80 800
0.80 0.80 0.80 0.75 0.99
0.1250 0.2600 0.2175 (middle) 0.0100 0.1150
Table 3 Parameters of vessel phantom c. Tissue name
Tissue composition
Average refractive index
Absorption coefficient (cm1)
Scattering coefficient (cm1)
Anisotropy factor
Thickness (cm)
Adventitia Intima/Media Fibro-fatty plaque Calcified plaque TCFA Lumen
Connective tissue Muscular tissue Fibro-fatty tissue Calcium Fibrous tissue Blood
1.41 1.41 1.46 1.46 1.46 1.34
0.20 0.20 0.45 0.05 0.01 0.99
20 20 80 80 80 800
0.80 0.80 0.80 0.80 0.75 0.99
0.1421 0.2150 0.0619 (middle) 0.0715 (middle) 0.0100 0.1429
40
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
2.2.5. FDTD solution of PA pressure The lossless propagation of PA waves in an acoustically homogeneous medium is physically described with the PA wave equation as [13,36]:
hl hθ hl kði þ 1; jÞΦði þ 1; jÞ þ kði; j þ 1ÞΦði; j þ 1Þ þ kði 1; jÞΦði 1; jÞ hθ hl hθ hθ hl hθ hl þ kði; j 1ÞΦði; j 1Þ kði þ 1; jÞ þ kði; j þ 1Þ þ kði 1; jÞ hl hθ hl hθ hθ þ kði; j 1Þ þ 2½μa ði; jÞc þ qhθ hl Φði; jÞ hl
r2 pðr; tÞ
¼ 0;
where Φ(i, j) is the optical radiance energy at the node (ihθ, jhl); μa(i, j) is the absorption coefficient at (ihθ, jhl) and k(iþ1, j), k(i, jþ1), k(i‒1, j) and k(i, j‒1) are the optical diffusion coefficient at ((iþ1)hθ, jhl), (ihθ, (jþ1)hl), ((i‒1)hθ, jhl) and (ihθ, (j‒1)hl), respectively. In the case of a boundary node P 2 D in Fig. 4(b), the second equation in Eq. (12) is rewritten as follows:
∂Φ c 1 Rf Φ: _ ¼ 2 1 þ Rf ∂n
(20)
where p(r, t) represents the initial acoustic pressure at r2Ω and time t, cs is the speed of sound, β is the thermal expansion coefficient, CP is the specific heat and δ(t) is the unit impulse function representing the instantaneous laser pulse. In accordance with the relationship among acoustic pressure, particle vibration velocity and medium density, Eq. (20) is rewritten as follows:
(14)
k
1 ∂2 pðr; tÞ β ¼ HðrÞ δ'ðtÞ; c2s ∂t2 Cp
8 ∂p ∂vθ ∂vl βc2 > > þ þ s HðrÞδðtÞ > ¼ ρc2s > > ∂t ∂θ ∂l CP > > > < ∂vθ 1 ∂p ; ¼ > ∂t ρ ∂θ > > > > > ∂vl > 1 ∂p > : ¼ ∂t ρ ∂l
(15)
Thus, the differential expression of Eq. (13) is ½Φði 1; jÞ Φði; jÞ ½Φði; j 1Þ Φði; jÞ jBEj þ kði; j 1Þ jABj hθ hl c 1 Rf Φði; jÞjEAj ½μa ði; jÞc þ qΦði; jÞjDj 2 1 þ Rf
(21)
where ρ is the medium density and vθ and vl are the particle vibration velocity along the θ-axis and l-axis, respectively. The acoustic pressure time series at θm is obtained by solving the FDTD expression of Eq. (21) [32]:
kði 1; jÞ
¼ 0; (16)
8 > > > > > > > 1 ðnþ12Þ 1 1 1 ðnþ12Þ 1 1 βc2 ði; jÞΔt ðnþ12Þ ðnþ12Þ > > pnþ1 ði; jÞ ¼ pn ði; jÞ ρc2s ði; jÞΔt v ; j v ; j þ v v þ s Hði; jÞδn i þ i i; j þ i; j > θ θ l l > Δθ 2 2 Δl 2 2 CP > > > > < 1 1 Δt ðn1Þ ðnþ1Þ ; ½pn ði; jÞ pn ði 1; jÞ vθ 2 i ; j ¼ vθ 2 i ; j > 2 2 ρ Δθ > > > > > > 1 1 Δt ðn1Þ ðnþ1Þ > > ¼ vl 2 i; j ½pn ði; jÞ pn ði; j 1Þ vl 2 i; j > > 2 2 ρΔl > > > > :
where (i, j) is the polar coordinate of the measuring location r; Δθ and Δl are the unit steps along the θ- and l-axis, respectively (here, Δθ ¼ Δl); Δt is the temporal step; n represents the discrete time point; cs(i, j) is the
where jABj, jBEj and jEAj are the lengths of AB, BE and EA, respectively, and jDj is the area of D. For a boundary node P 2 D facing the light source in Fig. 4(c), the third equation in Eq. (12) is rewritten as: k
∂Φ 1 1 Rf c 1 Rf ¼ 4Qs δðr r s Þ Φ: b 2 1 þ Rf 2 1 þ Rf ∂n
(22)
(17)
Thus, the differential expression of Eq. (13) is ½Φði þ 1; jÞ Φði; jÞ ½Φði; j þ 1Þ Φði; jÞ jBEj þ kði; j þ 1Þ jABj hθ hl c 1 Rf Φði; jÞjEAj ½μa ði; jÞc þ qΦði; jÞjDj 2 1 þ Rf 1 1 Rf : ¼ jEAj4Qs 2 1 þ Rf (18)
kði þ 1; jÞ
Finally, the theoretical absorbed energy distribution H(r) is obtained in accordance with HðrÞ ¼ μa ðrÞΦðrÞ:
(19)
Fig. 5. Change curve of the optical radiation energy (Fzr) and absorbed energy (Azr) with respect to the tissue thickness l along θ1 ¼ 0 and θ2 ¼ 45 for vessel phantom a. 41
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
speed of sound at (i, j); pn ði; jÞ represents the acoustic pressure at (i, j) and ðnÞ
2.3. Reconstruction of the optical absorption coefficient
ðnÞ
time n; vθ ði; jÞ and vl ði; jÞ are the vibration velocity along the θ-axis and l-axis at (i, j) and time n and H(i, j) is the absorbed energy density at (i, j).
2.3.1. Formation of the optical inverse problem The problem of recovering the optical absorption coefficient of wall tissues with known scattering and Gruneisen coefficient from PA pressure
Fig. 6. Change curve of the optical radiation energy (Fzr) and absorbed energy (Azr) with respect to the tissue thickness L along θ ¼ 0 –50 for (a) phantom a, (b) phantom b and (c) phantom c.
42
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
kXkTV ¼ jMXj;
measurements is expressed as the minimisation of a nonlinear leastsquare error:
2
b μ a ðrÞ ¼ argmin pðmÞ ðr; tÞ C⋅pðr; t; μa ðrÞ ; μa
where M is a sparse Ne N matrix with grid length as entries, Ne is the iteration time, N is the number of grids and j⋅j is the conventional L1 norm. The simplified TV-regularisation in Eq. (26) can be solved with the L–M algorithm by referring to the subgradient solution of L1regularisation [18,41]:
(23)
where b μ a ðrÞ is the estimated absorption coefficient at r, p(m)(r, t) and p(r, t, μa(r)) represent the measured and modelled acoustic pressure at r and time t and C is a calibration factor related to the acquisition system. The optimisation of Eq. (23) is a nonlinear ill-posed problem. Illposedness can be alleviated through regularisation. Current regularisation methods include L1, L2, Tikhonov and TV [37]. In most research on qPAT, optical parameters are usually assumed to be piecewise constants [38–40]. This assumption is also true for arterial vessels. Hence, in theory, TV-regularisation overperforms the L2-and Tikhonov-regularisation methods [18,41] for qIVPA. For the convenience of expression, we define X ¼ μa(r) and
(26)
∂ kXk1 ¼ signðXÞ; ∂X
(27)
where sign(⋅) is a sign function. 2.3.2. Bregman optimisation The Bregman [18] iteration expression of Eq. (25) is expressed as follows: Xnþ1 ¼ argmin½gðXÞ;
(28)
ðX;dÞ
f ðXÞ ¼ pm ðrÞ C⋅pðr; XÞ : SðXÞ ¼ f ðXÞT f ðXÞ
(24) where 8 α > gðXÞ ¼ SðXn Þ þ kMXn dn νn k22 hpn ; Xn i > > > 2 > <
MXnþ1 vn η ; max jMXnþ1 vn j; dnþ1 ¼ > > α jMXnþ1 vn j > > > : vnþ1 ¼ vn þ dnþ1 MXnþ1
The objective function to be optimised is rewritten as follows: b ¼ argmin SðXÞ þ ηkXkTV ; X
(25)
X
b is the estimate of X and η is the regularisation parameter. At the where X unstructured grids, the TV-regularisation term is simplified as follows [18]:
(29)
where α is the penalty parameter, and vnþ1 is the residual produced in the nth iteration and v0 ¼ 0.
Fig. 7. Distribution map of the optical radiation energy for (a) phantom a, (b) phantom b and (c) phantom c.
Fig. 8. Distribution map of the absorbed energy for (a) phantom a, (b) phantom b and (c) phantom c. 43
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
8 T yn sTn yn sTn sn sTn > 1 1 > > B ¼ I B I þ > nþ1 n T T < yn sn y sn yT sn
The detailed steps to solve Eq. (28) are as follows. Step 1 Set the initial guess of X, i.e., X0, n ¼ 0 and B1 0 ¼ I. Step 2 Update B1 n in accordance with
n
> sn ¼ Xnþ1 Xn > > > : yn ¼ rgnþ1 rgn
n
:
(30)
Fig. 9. Distribution map of the initial acoustic pressure for (a) phantom a, (b) phantom b and (c) phantom c.
Fig. 10. Distribution map of the absorption coefficient reconstructed through Bregman optimisation based on TV-regularisation for (a) phantom a, (b) phantom b and (c) phantom c.
Fig. 11. Change curve of the absorption coefficient and absorbed energy along a radius across the plaque centre (depicted with an arrow in Fig. 10) with respect to the tissue thickness l for (a) phantom a, (b) phantom b and (c) phantom c. ‘rAzr’, ‘Ua’ and ‘rUa’ represent the normalized absorbed energy, predefined absorption coefficient (initial guess), and estimated absorption coefficient through Bregman optimisation. 44
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Then, determine the search direction wn ¼ B1 n rgn where rgn is the gradient of g(X) at Xn
α
rgn ¼ rSn þ M ðMX dn vn Þ pn : 2 T
Step 3 Search for Xnþ1 along wn Xnþ1 ¼ Xn an B1 n rgn ;
(31)
(32)
where an is the step length. Step 4 If krgn k2 < 0:01, then the iteration stops; otherwise, return to Step 2.
Table 4 Measured and actual central thicknesses of the plaques.
3. Results
Vessel phantom
Plaque
Actual thickness (cm)
Measured thickness (cm)
a b c
Lipid pool Calcified plaque Fibro-fatty plaque
0.2267 0.2175 0.0619
0.23 0.21 0.07
The performance of the proposed method was validated with computer-simulated arterial vessel-mimicking phantoms. Several results are presented in this section. Tables 1–3 list the optical parameters and size of the vessel phantoms shown in Fig. 1 under the temperature of 25 C and wavelength of 532 nm in reference to [42]. These parameters
Fig. 12. Polar map of the absorbed energy distribution. (a) Phantom a (DE); (b) Phantom b (DE); (c) Phantom c (DE); (d) Phantom a (MC); (e) Phantom b (MC); (f) Phantom c (MC).
Fig. 13. Change curve of the normalized absorbed energy simulated with MC and DE and the predefined absorption coefficient for (a) phantom a (θ ¼ 38 ), (b) phantom b (θ ¼ 308 ) and (c) phantom c (θ ¼ 18 ). 45
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
of the positive or negative amplitude of PA signals at tissue boundaries during the discrete solution of the wave equation with FDTD. Fig. 10 shows the distribution map of the absorption coefficient recovered through Bregman optimisation based on TV-regularisation. The boundaries of the plaques are illustrated in Fig. 10. In the figure, the boundaries of the mixed plaques are emphasised. However, the plaque boundaries are blurred in Fig. 8, which shows that the shapes and sizes of all the plaques are consistent with the actual geometry of the plaques. Consequently, optical deposition distribution can only present a rough reference for the identification of atherosclerotic lesions. The position, size, shape and optical properties of the plaques provided by the absorption coefficient distribution are more accurate than those provided by the absorbed energy distribution. This conclusion is supported by Fig. 11, where the distance between the adjacent troughs of the curves determines the size of a plaque. To quantitatively evaluate the accuracy of absorption coefficient recovery, the central thicknesses of the plaques shown in Fig. 10 were measured and compared with the real values listed in Tables 1–3. The results shown in Table 4 imply that the measured thicknesses of phantom a and b are consistent with the actual thicknesses of phantoms a and b. For phantom c, the measured thickness of the fibro–fatty plaque near the wall is consistent with the true value. However, the thickness of the calcified plaque distant from the lumen is larger than the actual value.
Fig. 14. Change curve of the absorption coefficient with respect to the tissue thickness l along θ1 ¼ 0 (normal vessel wall tissue) and θ2 ¼ 45 (including plaque tissue) for phantom a. ‘Ua’ and ‘rUa’ represent the predefined (initial guess) and reconstructed absorption coefficient. (a) θ1 ¼ 0 , Bregman based on L2; (b) θ1 ¼ 0 , Bregman based on TV; (c) θ2 ¼ 45 , Bregman based on L2; (d) θ2 ¼ 45 , Bregman based on TV.
were used in the forward simulation experiments and were also the initial guesses in the optimisation of Eq. (23). Figs. 5–8 show the distribution of the optical radiation energy and absorbed energy obtained through the forward simulation. The results shown in Fig. 5 are consistent with the geometry and tissue components shown in Fig. 1(a), where the lumen, intima/media and adventitia are found along θ1 ¼ 0 and the lumen, TCFA, lipid pool, intima/media and adventitia are found along θ1 ¼ 45 . The results shown in Figs. 5–8 indicate that the optical radiation energy is related to the absorption and scattering coefficients. However, the absorbed energy is dependent on the absorption coefficient. The radiation and absorbed energy attenuate to 0 with increasing detection depth. The optical deposition in the lumen and lipid and fibro–fatty plaques is higher than that in other tissues given the larger absorption coefficients of the lumen and plaques. Moreover, the photons absorbed by the lumen may be scattered to other tissues because the lumen is filled with blood, which has a large scattering coefficient. The intima/media absorbs partial light energy, and partial photons scatter to other tissues. The optical deposition in the adventitia is nearly 0. Fig. 9 shows the distribution map of the initial acoustic pressure after pseudocolour coding. The ring-down artefacts surrounding the catheter and located between the lumen and wall are caused by the normalisation
4. Discussion 4.1. Comparison of DE and MC The absorbed energy distributions simulated with DE and MC, a random statistical method for solving the RTE [43,44] were compared. The results are shown in Figs. 12 and 13. The figures indicate that the optical deposition obtained with DE distinguishes the range of plaque tissues more accurately than that obtained with MC. 4.2. Comparison of L2-and TV-regularisation The absorption coefficients recovered via Bregman optimisation based on L2- [18] and TV-regularisation were compared, and the results are shown in Figs. 10, 14 and 15. The curves in Fig. 14 exhibit noticeable peaks at the boundary of the lumen–intima/media, intima/media–plaque and plaque–adventitia. The calcified plaque can be identified and localised through TV- and L2-regularisation reconstruction. The quantification results listed in Table 5 suggest that the iteration times of TV-regularisation are lower than those of L2-regularisation although the estimated absorption coefficients of the two methods are similar.
Fig. 15. Distribution map of the absorption coefficient recovered via Bregman optimisation based on L2-regularisation for (a) phantom a, (b) phantom b and (c) phantom c. 46
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
or L2-regularisation are provided in this section. The L–M iteration expression of Eq. (25) is
4.3. Comparison of L–M and Bregman optimisation The results for the comparison of the absorption coefficient distribution recovered via L–M [28] and Bregman optimisation based on TV-
Table 5 Absorption coefficients recovered via L–M or Bregman optimisation based on L2-or TV-regularisation for phantom a. No.
Measuring location (θ,l)
Predefined
μa
L-M
2 3 4 5 6 7 8 9 10
(45 ,0.05) (0 , 0.05) (45 ,0.10) (0 , 0.10) (45 ,0.115) (0 ,0.115) (45 ,0.125) (0 ,0.125) (45 ,0.15) (0 ,0.15) (45 ,0.20) (0 ,0.20) (45 ,0.25) (0 ,0.25) (45 ,0.367) (0 ,0.367) (45 ,0.45) (0 ,0.45) (45 ,0.50) (0 ,0.50)
0.99 0.99 0.99 0.99 0.99 0.99 0.10 0.20 0.60 0.20 0.60 0.20 0.60 0.20 0.60 0.20 0.20 0.20 0.20 0.20
TV
L2
Estimated
μa
Iteration times
0.72 0.73 0.98 0.98 0.92 0.92 0.12 0.43 0.45 0.24 0.58 0.21 0.59 0.17 0.48 0.18 0.23 0.18 0.21 0.19
28 25 26 25 25 24 26 25 25 23 26 24 23 24 28 26 25 23 26 25
Estimated 1
Bregman
L2
μa
Iteration times
0.74 0.75 0.99 0.99 0.92 0.93 0.11 0.40 0.45 0.21 0.55 0.20 0.60 0.19 0.49 0.19 0.21 0.19 0.20 0.20
13 12 15 10 14 11 12 13 12 12 13 10 12 11 11 12 10 10 11 10
Estimated
TV
μa
Iteration times
Estimated
μa
Iteration times
0.74 0.75 0.98 0.98 0.93 0.92 0.11 0.42 0.47 0.22 0.56 0.21 0.60 0.18 0.49 0.19 0.21 0.19 0.20 0.19
27 24 28 26 24 25 25 24 25 24 27 25 25 24 28 26 26 24 25 24
0.75 0.76 0.99 0.99 0.93 0.93 0.10 0.40 0.48 0.20 0.56 0.20 0.60 0.19 0.50 0.20 0.20 0.20 0.20 0.20
14 12 15 10 13 13 11 12 12 11 12 10 13 10 10 11 11 10 10 10
Fig. 16. Distribution map of the absorption coefficient recovered via L–M optimisation. (a) Phantom a (L–M based on L2); (b) Phantom b (L–M based on L2); (c) Phantom c (L–M based on L2); (d) Phantom a (L–M based on TV); (e) Phantom b (L–M based on TV); (f) Phantom c (L–M based on TV). 47
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
Fig. 17. Change curve of the absorption coefficient along a radius across the plaque centre (depicted with an arrow in Fig. 10) with respect to the tissue thickness l for (a) phantom a, (b) phantom b and (c) phantom c. ‘Ua’, ‘LM’ and ‘Bregman’ represent the predefined absorption coefficient (initial guess), recovered absorption coefficient via L–M and Bregman optimisation, respectively.
8 2 > < ΔXnþ1 ¼ argmin krSðXn Þ⋅ΔX ½SðXÞ SðXn Þ þ ωn k2 þ ηkΔXkTV ΔX
ω ¼ ωn þ ½SðXÞ SðXn ÞSðXÞ SðXn Þ rSðXn Þ⋅ΔXnþ1 > : nþ1 Xnþ1 ¼ Xn þ ΔXnþ1
and to discriminate between normal arterial tissues and lipid-rich plaques on the basis of the strong absorption of lipids has attracted increased attention [45,46]. At 1.2 μm, the IVPA system requires the laser energy of at least 1.2 mJ to distinguish between lipids and other tissues. However, it only requires the laser energy of 0.4 mJ at 1.7 μm. Therefore, IVPA images with sufficient contrast can be produced with less energy at 1.7 μm than at other NIR wavelengths [47]. Future works will demonstrate the validity of the proposed method at 1.7 μm by adjusting the optical parameter setting in the forward simulation. Moreover, a simple photon propagation pattern was adopted in the forward model. The wall tissues were assumed to be uniformly illuminated in a 2-D space parallel to the imaging plane by ignoring the axial influence of the narrow beam laser pulse and nonuniform distribution of light in the medium. The possible changes in the Gruneisen coefficient and acoustic inhomogeneity were ignored. The ultrasonic detector was simplified as an ideal point-like detector by ignoring the aperture effect related to its size. The simplification introduced differences between the theoretical and measured PA pressure and imposed possible limitations on the practical applications of the proposed method.
; (33)
where ΔX is the increment of X. The detailed steps to solve Eq. (33) are as follows. Step 1 Set the initial guess of X, i.e., X0, n ¼ 0, ω02(0,1), τ > 1, v > 1 and ε > 0. Step 2 Determine f(Xn), S(Xn), rf ðXn Þ and rSðXn Þ ¼ ðrf ðXn ÞÞT ⋅f ðXn Þ. Step 3 Solve ½Q þ τIΔX ¼ rSðXn Þ to obtain ΔX, where Q ¼ ðrf ðXn ÞÞT ⋅rf ðXn Þ and I is the unit matrix. Step 4 Calculate Xnþ1 ¼ Xn þ ΔX. If the absolute error is less than ε, the iteration ends; otherwise, return to Step 5. Step 5 If SðXnþ1 Þ < SðXn Þ þ ωn ½rSðXn ÞT ΔX, then τ←τ/v and proceed to Step 6; otherwise, τ←τ⋅v and return to Step 3. Step 6 k←kþ1 and return to Step 2. The results are shown in Figs. 10, 15 and 16. The visualisation of the images reconstructed with TV- and L2-regularisation combined with the same optimisation method is similar. However, by referring to the actual geometry of the phantoms, the Bregman reconstruction provides visually better images than the L–M reconstruction based on the same regularisation. The absorption coefficients reconstructed with L–M and Bregman optimisation were compared with the predefined values listed in Tables 1–3. Fig. 17 indicates the similarity shared by L–M and Bregman optimisation. However, the iteration times of Bregman optimisation are higher than those of L–M optimisation. Considering the precision and computational cost, L–M may be more suitable than Bregman optimisation for the individual recovery of absorption coefficients, whereas Bregman optimisation may be a better choice than L–M for the simultaneous estimation of absorption and scattering coefficients.
5. Conclusion In this preliminary study, a one-step qIVPA method for the singlewavelength IVPA system was designed. The method successfully recovered the transversal maps of optical absorption coefficient distribution on vascular cross-sections from PA boundary measurements. To implement fully quantitative optical inversion, the model-based inversion utilised the DA to the RTE to model light transport in a medium with a known scattering coefficient and accounted for acoustic propagation in acoustically homogeneous lossless medium. The L–M or Bregman optimisation based on TV-regularisation was applied to minimise a nonlinear leastsquare error function to quantify the difference between the measured and modelled acoustic pressure series. The experimental results obtained for computer-simulated arterial phantoms indicate that the lumen, TCFA, lipid plaque, calcified plaque, fibro–fatty plaque and vessel-wall intima/ media/adventitia can be clearly identified in an absorption coefficient distribution map. Moreover, the location, size and contour of the TCFA, lipid plaque, calcified plaque and fibro–fatty plaque can be clearly distinguished.
4.4. Other discussions The theoretical acoustic pressure generated by vessel-wall tissues were obtained through the forward simulation of excitation laser transport and tissue absorption. A narrow beam of collimated excitation light was assumed to be normally incident upon the wall tissues through the blood. The influence of luminal blood on light transport was ignored. In practical applications, haemoglobin absorption is high in the visible range (410–680 nm). Nevertheless, IVPA can discriminate between normal and atheromatous tissues at these wavelengths. In recent years, the use of experimental IVPA systems at near-infrared excitation wavelengths, such as 1.2 or 1.7 μm, to minimise the effect of blood absorption
Conflicts of interest None Declared. Acknowledgments This work was supported by National Nature Science Foundations of 48
S. Zheng, Z. Lan
Computers in Biology and Medicine 97 (2018) 37–49
China (no. 61372042).
[23] O. Nyk€anen, A. Pulkkinen, T. Tarvainen, Quantitative photoacoustic tomography augmented with surface light measurements, Biomed. Optic Express 8 (10) (2017) 4380–4395. [24] A.Q. Bauer, R.E. Nothdurft, T.N. Erpelding, et al., Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography, J. Biomed. Optic. 16 (9) (2011), 096016. [25] J.R. Rajian, P.L. Carson, X. Wang, Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent, Optic Express 17 (6) (2009) 4879–4889. [26] K. Daoudi, A. Hussain, E. Hondebrink, et al., Correcting photoacoustic signals for fluence variations using acousto-optic modulation, Optic Express 20 (13) (2012) 14117–14129. [27] Y. Zhen, J. Huabei, A calibration-free, one-step method for quantitative photoacoustic tomography, Med. Phys. 39 (11) (2012) 6895–6899. [28] T. Ding, K. Ren, S. Vallelian, A one-step reconstruction algorithm for quantitative photoacoustic imaging, Inverse Probl. 31 (9) (2015) 65003–65023. [29] M. Haltmeier, L. Neumann, S. Rabanser, Single-stage reconstruction algorithm for quantitative photoacoustic tomography, Inverse Probl. 31 (6) (2015) 65005–65028. [30] G. Bal, A. Moradifam, Photo-acoustic tomography in a rotating measurement setting, Inverse Probl. 32 (10) (2016), 105012. [31] A. Pulkkinen, B.T. Cox, S.R. Arridge, et al., Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography, IEEE Trans. Med. Imag. 35 (11) (2016) 2497–2508. [32] Z. Sun, Y. Yuan, D. Han, A computer-based simulator for intravascular photoacoustic images, Comput. Biol. Med. 81 (2017) 176–187. [33] T. Tarvainen, Computational Methods for Light Transport in Optical Tomography, University of Kuopio, Kuopio, 2006. [34] S. Li, B. Montcel, Z. Yuan, et al., Multigrid-based reconstruction algorithm for quantitative photoacoustic tomography, Biomed. Optic Express 6 (7) (2015) 2424–2434. [35] F. Gao, Y. Tanikawa, H. Zhao, et al., Semi-three-dimensional algorithm for timeresolved diffuse optical tomography by use of the generalized pulse spectrum technique, Appl. Optic. 41 (34) (2002) 7346–7358. [36] L.V. Wang, Biomedical Optics: Principles and Imaging. Wang Lihong. Photoacoustic Tomography, John Wiley and Sons, New Jersey, 2007, pp. 283–321. [37] H. Gao, H. Zhao, Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization, Optic Express 18 (3) (2010) 1854–1871. [38] W. Neater, O. Scherzer, Quantitative photoacoustic tomography with piecewise constant material parameters, SIAM J. Imag. Sci. 7 (3) (2014) 1755–1774. [39] E. Beretta, M. Muszkieta, W. Naetar, et al., A Variational Method for Quantitative Photoacoustic Tomography with Piecewise Constant Coefficients, arXiv: 1509.04982, 2015. [40] A. De Cezaro, A. Leit~ao, X. Tai, On piecewise constant level-set (PCLS) methods for the identification of discontinuous parameters in ill-posed problems, Inverse Probl. 29 (29) (2012), 015003. [41] H. Gao, H. Zhao, Multilevel bioluminescence tomography based on radiative transfer equation part 2: total variation and l1 data fidelity, Optic Express 18 (3) (2010) 2894–2912. [42] S.L. Jacques, Optical properties of biological tissues: a review, Phys. Med. Biol. 58 (11) (2013) R37–R61. [43] W. Li, L. Lin, L. Bao, et al., Monte Carlo simulation of photon migration in multicomponent media, Opt. Quant. Electron. 47 (7) (2015) 1919–1931. [44] L. Wang, S.L. Jacques, L. Zheng, MCML-Monte Carlo modeling of light transport in multi-layered tissues, Comput. Meth. Progr. Biomed. 47 (2) (1995) 131–146. [45] T.J. Allen, P.C. Beard, Photoacoustic characterisation of vascular tissue at NIR wavelengths, in: A.A. Oraevsky, L.V. Wang (Eds.), Proceedings of SPIE International Conference on Photons Plus Ultrasound: Imaging and Sensing 2009, vol. 7177, 25 Feb. 2009 id: 71770A. [46] N. Dana, T. Sowers, A. Karpiouk, et al., Optimization of dual-wavelength intravascular photoacoustic imaging of atherosclerotic plaques using Monte Carlo optical modeling, J. Biomed. Optic. 22 (10) (2017) 1–12. [47] K. Jansen, M. Wu, A.F.W. van der Steen, et al., Photoacoustic imaging of human coronary atherosclerosis in two spectral bands, Photoacoustics 2 (1) (2014) 12–20.
References [1] W. Min, F.W.S. Antonius, R. Evelyn, et al., Emerging technology update intravascular photoacoustic imaging of vulnerable atherosclerotic plaque, Intervent Cardiol. 11 (2) (2016) 120–123. [2] J. Yao, L.V. Wang, Breakthroughs in photonics photoacoustic tomography, IEEE Photon. J. 6 (2) (2014), 0701006. [3] S.Y. Nam, E. Chung, L.J. Suggs, et al., Combined ultrasound and photoacoustic imaging to noninvasively assess burn injury and selectively monitor a regenerative tissue-engineered construct, Tissue Eng. C Meth. 21 (6) (2015) 557–566. [4] B. Wang, S. Emelianov, Thermal intravascular photoacoustic imaging, Opt. Soc. Am. 2 (11) (2011) 3072–3078. [5] D. Verya, W. Min, J. Nico de, et al., Frequency analysis of the photoacoustic signal generated by coronary atherosclerotic plaque, Ultrasound Med. Biol. 42 (8) (2016) 2017–2025. [6] B. Wang, E. Yantsen, T. Larson, et al., Plasmonic intravascular photoacoustic imaging for detection of macrophages in atherosclerotic plaques, Nano Lett. 9 (6) (2009) 2212–2217. [7] W. Min, S. Geert, L. Matija, et al., Real-time volumetric lipid imaging in vivo by intravascular photoacoustics at 20 frames per second, Biomed. Optic Express 8 (2) (2017) 943–953. [8] X. Bai, X. Gong, W. Hau, et al., Intravascular optical-resolution photoacoustic tomography with a 1.1 mm diameter catheter, PLoS One 9 (3) (2014) e92463. [9] Y. Liang, L. Jin, L. Wang, et al., Fiber-laser-based ultrasound sensor for photoacoustic imaging, Sci. Rep. 7 (2017) 40849. [10] J. Zhang, S. Yang, X. Ji, et al., Characterization of lipid-rich aortic plaques by intravascular photoacoustic tomography: ex vivo and in vivo validation in a rabbit atherosclerosis model with histologic correlation, J. Am. Coll. Cardiol. 64 (4) (2014) 385–390. [11] H. Qin, T. Zhou, S. Yang, et al., Gadolinium (III)-gold nanorods for MRI and photoacoustic imaging dual-modality detection of macrophages in atherosclerotic inflammation, Nanomedicine 8 (10) (2013) 611–1624. [12] Y.L. Sheu, C.Y. Chou, B.Y. Hsieh, et al., Application of limited-view image reconstruction method to intravascular photoacoustic tomography, in: A.A. Oraevsky, L.V. Wang (Eds.), Proceedings of SPIE International Conference on Photons Plus Ultrasound: Imaging and Sensing 2010, vol. 7564, 24 Feb. 2010 id: 75640B. [13] Y.L. Sheu, C.Y. Chou, B.Y. Hsieh, Image reconstruction in intravascular photoacoustic imaging, IEEE Trans. Ultrason. Ferroelectrics Freq. Contr. 58 (2011) 2067–2077. [14] Z. Sun, D. Han, J. Wang, Review of image reconstruction for intravascular photoacoustic imaging, Opto-Electronic Eng. 42 (3) (2015) 20–27. [15] Z. Sun, D. Han, Y. Yuan, 2-D image reconstruction of photoacoustic endoscopic imaging based on time-reversal, Comput. Biol. Med. 76 (2016) 60–68. [16] N. Song, C. Deumic, S.A. Da, Considering sources and detectors distributions for quantitative photoacoustic tomography, Biomed. Optic Express 5 (11) (2014) 3960–3974. [17] H. Gao, S. Osher, H. Zhao, Mathematical Modeling in Biomedical Imaging II, Springer, Berlin/Heidelberg, 2012, pp. 131–158. [18] H. Gao, H. Zhao, S. Osher, Bregman methods in quantitative photoacoustic tomography, CAM Rep. 30 (6) (2010) 3043–3054. [19] B. Cox, J.G. Laufer, S.R. Arridge, et al., Quantitative spectroscopic photoacoustic imaging: a review, J. Biomed. Optic. 17 (6) (2012), 061202. [20] X. Zhang, W. Zhou, X. Zhang, et al., Forward-backward splitting method for quantitative photoacoustic tomography, Inverse Probl. 30 (2014), 125012. [21] R. Hochuli, S. Powell, S. Arridge, et al., Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance, J. Biomed. Optic. 21 (2016), 126004. [22] M. Venugopal, E.P. Van, S. Manohar, et al., Quantitative photoacoustic tomography by stochastic search: direct recovery of the optical absorption field, Optic Lett. 41 (18) (2016) 4202–4205.
49