Analytical model of optical fluence inside multiple cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium for quantitative photoacoustic imaging

Shengfu Li; Bruno Montcel; Wanyu Liu; Didier Vray

doi:10.1364/OE.22.020500

1. Introduction

The forward model of optical fluence in biological tissues is an important topic in diffuse optical tomography (DOT) [1–3] and in quantitative photoacoustic (PA) imaging [4–6]. There are three types of forward models: (1) Monte Carlo (MC) methods [7–9] are based on a stochastic model and therefore need significant computation time to achieve sufficient precision; (2) numerical methods, e.g., finite element (FEM) [10–16], involve matrix inversion and therefore are computationally costly; (3) analytical methods involve less computation time, and thus should be more suitable for real-time PA imaging. However, analytical modes only exist for simple cases.

Several studies have presented analytical models of optical fluence for homogeneous and inhomogeneous media and these models have been validated with experiments and simulations. Patterson et al. developed analytical models for semi-infinite and finite homogeneous tissue slab [17,18]. Boas et al. developed analytical models for an inhomogeneity embedded in an otherwise homogeneous turbid medium. They introduced boundary conditions that take into account the index of refraction mismatch between the embedded inhomogeneity and the background medium [19,20] and their model has been validated experimentally [21]. Some models were also based on the Born approximation [9], which is limited to small fluctuations of the optical properties. Ripoll et al. developed analytical models based on the Kirchhoff approximation [22,23], assuming that the total intensity at a certain point in the medium is equal to the sum of the incident field and the wave reflected from the plane tangent to the interface. These studies focused on DOT, and therefore the optical fluence is mainly investigated at the detector position. Relatively little work has been done on the analytical model of multiple inhomogeneities embedded in an otherwise homogeneous turbid medium, probably because this model can hardly be integrated in a DOT reconstruction process for tissues. Indeed, optical properties of tissues are very inhomogeneous. Furthermore, for DOT, the optical detectors are positioned onto the tissue surface, and therefore surface boundary conditions are mandatory. However, we believe that this model could be relevant for forward problems and inverse problems of PA imaging, since PA waves are only generated by absorbers inside tissues. The absorbers in tissues are mainly localized in anatomical structures. In particular, the main endogenous contrast for PA imaging is the blood vessels due to the strong absorption of hemoglobin. Hence, considering a model of inhomogeneities of absorption could be more relevant in PA imaging than in DOT.

Blood vessels are roughly cylindrical and have a much stronger absorption coefficient than the background tissues. Furthermore, hemoglobin concentration can be assumed to be uniform inside the vessel. Therefore, the blood vessels can be considered as “cylindrical inhomogeneities.” In this paper we present an analytical model for multiple parallel cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium to analyze the factors that influence optical fluence distribution. The practical applications of this model are limited to the biological tissues where blood vessels can be assumed as “parallel vessels”.

Analytical model for such configuration is a well-known problem in acoustics, electromagnetics and the light scattering of particles [24–31]. However, all these related works focused only on the analytical model corresponding to the physical quantity outside the “inhomogeneities”. It was assumed that a plane wave was incident on the “inhomogeneities” in acoustics, electromagnetics [24,25,27–31]. In the research of the light scattering of cylindrical particles, only single scattering has been considered based on the assumption that the interaction of particles can be neglected [26]. However, in biological tissues, PA waves are generated by the local absorption contrast introduced by “inhomogeneities” in biological tissues, thus the optical fluence inside inhomogeneities is a key point to analyze the PA waves. The objective of this study is to present a model that focuses on the optical fluence both outside and inside the inhomogeneities embedded in turbid medium.

Our model arises out of the solution of the diffusion equation that is the approximation of the Boltzmann transport equation. The model represents the optical fluence distribution in tissues as a function of the reduced scattering and absorption coefficients. The migration of light through biological tissues in this model can be treated as a wave that is called a diffuse photon density wave (DPDW) in literature [19]. In tissues with homogenous optical properties, the incident DPDW of light source propagates unobstructed. However, the presence of “inhomogeneities” in the optical properties leads to a distortion, modeled as a scattered DPDW [20]. Light that diffuses through a medium including “inhomogeneities” can be considered as a superposition of the incident wave and scattered waves [32,33]. In view of high scattering properties of biological tissues, this incident wave cannot be modeled as a plane wave and the multiple scattering of “inhomogeneities” needs to be considered. The interplay between cylindrical inhomogeneities is incorporated by using new boundary conditions that take into account the multiple scattering of inhomogeneities. The model has been compared with the numerical solution obtained by use of NIRFAST software [10,11]. Close agreement of fluence inside inhomogeneities, which is the source of PA signal, permits the use of our model as a forward model needed in PA imaging.

2. Theory

The geometry of the problem is illustrated in Fig. 1.A point light source is incident on cylindrical inhomogeneities embedded in an otherwise homogeneous infinite turbid medium. The inhomogeneities have different optical parameter from background medium.

Fig. 1 Horizontal cross section of the geometry. Cylindrical inhomogeneities of absorption coefficient embedded in an otherwise homogeneous infinite turbid medium. i is index of i^th cylinder and i’ represents an arbitrary remaining cylinders. Cylindrical coordinate systems have their origins at the center of each inhomogeneity. The source is noted by S. P indicates an arbitrary point within the medium, and its coordinates values with respect to each cylinder coordinate system are indicated. The z-axis is along the axis of each inhomogeneity and comes out of the cross section.

Download Full Size | PDF

2.1 Incident waves

In PA imaging the acoustic propagation occurs on a timescale several orders of magnitude longer than the heat deposition. Therefore the time-integrated absorbed power density (i.e., the absorbed energy density) is the quantity of interest [5]. In most biological tissues the optical fluence obeys the diffusion equation [20]. In an infinite homogeneous medium, the time-independent diffusion equation has the form [34],

\nabla^{2} Φ (r) + k^{2} Φ (r) = - \frac{S (r)}{D}

and

k = j \cdot \sqrt{3 \cdot (μ_{s}^{'} + μ_{a}) \cdot μ_{a}}

where

j^{2} = - 1

,

μ_{s}^{'}

is the reduced scattering coefficient,

μ_{a}

is the absorption coefficient,

D = {[3 (μ_{s}^{'} + u_{a})]}^{- 1}

is the diffusion coefficient, and

S (r)

is the source term.

The optical fluence at a detector position $r_{d}$ , $Φ_{i n c} (r_{s}, r_{d})$ , in an infinite homogeneous medium for a point source at $r_{s}$ has been given by other studies [17,19]. In this study we use the incident optical fluence at an arbitrary position r inside the medium,

Φ_{i n c} (r_{s}, r) = \frac{S}{4 π D | r - r_{s} |} \exp (j k | r - r_{s} |)

In the presence of cylindrical inhomogeneities, it is natural to analyze the problem in cylindrical coordinates that are defined with respect to each center of the cylinders, as in Fig. 1. The z-axis of each cylindrical coordinate is the axis of each cylinder. The polar axis of all coordinates is parallel to the x-axis. Important notations are explained in Fig. 1, in particular i is index of i^th cylinder and $i^{'}$ represents an arbitrary remaining cylinders, with respect to the i^th cylinder.

Based on the works of Jackson [35], the cylindrical expansion of $Φ_{i n c} (r_{s}, r)$ with respect to the i^th inhomogeneity center can be written as:

{}^{i}Φ_{i n c} = \frac{S}{2 π^{2} D_{o u t}} \sum_{n = - \infty}^{+ \infty} \int_{0}^{\infty} d p \exp [j \cdot n ({}^{i}θ - {}^{i}θ_{s})] \cos (p \cdot z) I_{n} ({}^{i}r_{<}) K_{n} ({}^{i}r_{>})

D_{o u t}

is the diffusion coefficient of the background.

I_{n}

and

K_{n}

are modified Bessel functions of the first and the second kind, respectively.

^{i} r_{<} = {}^{i}ρ_{<} \cdot {(p^{2} - k_{o u t}^{2})}^{1 / 2}

,

^{i} r_{>} = {}^{i}ρ_{>} \cdot {(p^{2} - k_{o u t}^{2})}^{1 / 2}

,

{}^{i}ρ_{<}

and

{}^{i}ρ_{>}

are, respectively, the smaller and the larger values of

{}^{i}ρ_{s}

and

{}^{i}ρ

, where

{}^{i}ρ_{s}

is the radial distance of the light source and

{}^{i}ρ

is an arbitrary position (indicated by p) in the medium with respect to the i^th inhomogeneity center.

n

and

p

arise from the separation of variables.

k_{o u t}

is the wave number of the background.

{}^{i}θ

and

{}^{i}θ_{s}

are the azimuth angles of a given position in the medium and of the light source, respectively, with respect to the i^th inhomogeneity center.

2.2 Scattered waves

If the individual cylindrical inhomogeneities are sufficiently far from each other, the interplay between cylindrical inhomogeneities can be neglected and the boundary conditions given by Li et al. [36] are valid.

If the cylindrical inhomogeneities are closer, the physical interpretation of their interplay is shown in Fig. 2. $Φ_{i n c}$ , the optical fluence incident on the i^th inhomogeneity, produces the first-order scattered wave, $^{i} Φ_{s c a t}^{(1)}$ (see Fig. 2(a)). In the same manner all the remaining cylinders produce first-order scattered waves, $\sum_{i^{'}}^{i^{'}} Φ_{s c a t}^{(1)}$ , which are incident on the i^th cylinder and produce the second-order scattered wave, $^{i} Φ_{s c a t}^{(2)}$ . Proceeding in this way, one can obtain the l-order scattered waves, $^{i} Φ_{s c a t}^{(l)}$ , with l from 1 to infinity.

Fig. 2 Physical interpretation of the different orders of scattering waves by the i^th inhomogeneity: (a) the incident wave produces the first-order scattered wave by the i^th inhomogeneity; (b) the sum of first-order scattered waves from the remaining cylinders produces the second-order scattered wave by the i^th inhomogeneity.

Download Full Size | PDF

According to this physical interpretation, the second order of scattering, $^{i} Φ_{s c a t}^{(2)}$ , can be interpreted as arising out of the consecutive scattering of the incident wave by two different inhomogeneities. $^{i} Φ_{s c a t}^{(l)}$ arises out of the scattering of the incident wave by l consecutive scattering process. Therefore, $^{i} Φ_{s c a t}^{(l)}$ approaches 0 when l approaches infinity. This enables the truncation by considering only the first t orders of scattering to generate a solution, where t is a positive integer.

On the basis of the general solution given by Walker et al. [20] and of Eq. (3), the l-order scattered wave from the i^th cylinder can be expressed as follows:

{}^{i}Φ_{s c a t}^{(l)} = \sum_{n = - \infty}^{\infty} \int_{0}^{\infty} d p \exp (j \cdot n \cdot {}^{i}θ) \cos (p \cdot z) \cdot {}^{i}B_{n}^{(l)} (p) \cdot K_{n} ({}^{i}x)

And the sum of the first t orders of scattering waves from the i^th cylinder ⁱSC^t, can be expressed as follows:

^{i} S C^{t} = \sum_{n = - \infty}^{\infty} \int_{0}^{\infty} d p \exp (j \cdot n \cdot {}^{i}θ) \cos (p \cdot z) \cdot {}^{i}B_{n}^{t} (p) \cdot K_{n} ({}^{i}x)

With

{}^{i}B_{n}^{t} (p) = \sum_{l = 1}^{t} {}^{i}B_{n}^{(l)} (p) .

Correspondingly, the fluence inside the i^th cylinder has the form,

{}^{i}Φ_{i n}^{t} = \sum_{n = - \infty}^{\infty} \int_{0}^{\infty} d p \exp (j \cdot n \cdot {}^{i}θ) \cos (p \cdot z) \cdot {}^{i}C_{n}^{t} (p) \cdot I_{n} ({}^{i}y)

Here,

{}^{i}x = {}^{i}ρ \cdot {(p^{2} - k_{o u t}^{2})}^{1 / 2}

,

{}^{i}y = {}^{i}ρ \cdot {(p^{2} - {}^{i}k_{i n}^{2})}^{1 / 2}

,

{}^{i}k_{i n}

is the wave number inside the i^th cylinder,

{}^{i}B_{n}^{t} (p)

,

{}^{i}B_{n}^{(l)} (p)

and

{}^{i}C_{n}^{t} (p)

are unknown coefficients that are determined by use of boundary conditions.

2.3 Boundary conditions

As explained by Fig. 2, we can see that ${}^{i}Φ_{i n c}$ produces $^{i} Φ_{s c a t}^{(1)}$ and that $\sum_{i^{'}}^{i^{'}} Φ_{s c a t}^{(1)}$ produces $^{i} Φ_{s c a t}^{(2)}$ . Thus, ${}^{i}Φ_{i n c} + \sum_{i^{'}}^{i^{'}} Φ_{s c a t}^{(1)}$ produces $^{i} S C^{2} = \sum_{l = 1}^{2}^{i} Φ_{s c a t}^{(l)}$ . Proceeding in this way, and considering the first t orders of scattered waves, we can see that ${}^{i}Φ_{i n c} + \sum_{i^{'}} {}^{i^{'}}S C^{t - 1}$ produces $^{i} S C^{t} = \sum_{l = 1}^{t}^{i} Φ_{s c a t}^{(l)}$ . Thus the fluence outside the i^th cylinder can be written as:

{}^{i}Φ_{o u t}^{t} = {}^{i}Φ_{i n c} + {}^{i}S C^{t} + \sum_{i^{'}} {}^{i^{'}}S C^{t - 1}

It is noteworthy that ${}^{i}Φ_{o u t}^{t}$ include three terms that are incident waves, the sum of the first t orders of scattering waves from the cylinder in question and the sum of the first t-1 orders of scattering waves from the remaining cylinders, respectively. Boundary conditions require that: (1) the flux normal to the boundary of the i^th cylinder must be continuous; (2) the optical fluence must be continuous across the boundary of the i^th cylinder. They can be written as:

\begin{array}{l} D_{o u t} {[\frac{\partial}{\partial {}^{i}ρ} {}^{i}Φ_{o u t}^{t}]}_{{}^{i}ρ = {}^{i}a} = {}^{i}D_{i n} {[\frac{\partial}{\partial {}^{i}ρ} {}^{i}Φ_{i n}^{t}]}_{{}^{i}ρ = {}^{i}a} \\ {[{}^{i}Φ_{o u t}^{t}]}_{{}^{i}ρ = {}^{i}a} = {[{}^{i}Φ_{i n}^{t}]}_{{}^{i}ρ = {}^{i}a} \end{array}

Solving the system of linear equations we have,

{}^{i}B_{n}^{t} (p) = \frac{S}{2 π^{2} D_{o u t}} \cdot \frac{D_{o u t} \cdot (S_{n} (a) + S C_{n}) \cdot I_{n}^{'} (y_{b}) \cdot y_{b} - {}^{i}D_{i n} \cdot (S_{n}^{'} (a) + S C_{n}^{'}) \cdot I_{n} (y_{b}) \cdot x_{b}}{D_{o u t} \cdot K_{n}^{'} (x_{b}) \cdot I_{n} (y_{b}) \cdot x_{b} - {}^{i}D_{i n} \cdot K_{n} (x_{b}) \cdot I_{n}^{'} (y_{b}) \cdot y_{b}}

and

{}^{i}C_{n}^{t} (p) = \frac{S}{2 π^{2} D_{o u t}} \cdot \frac{D_{o u t} \cdot (S_{n} (a) + S C_{n}) \cdot K_{n}^{'} (x_{b}) \cdot x_{b} - D_{o u t} \cdot (S_{n}^{'} (a) + S C_{n}^{'}) \cdot K_{n} (x_{b}) \cdot x_{b}}{D_{o u t} \cdot K_{n}^{'} (x_{b}) \cdot I_{n} (y_{b}) \cdot x_{b} - {}^{i}D_{i n} \cdot K_{n} (x_{b}) \cdot I_{n}^{'} (y_{b}) \cdot y_{b}}

Where

x_{b} = {}^{i}a \cdot \sqrt{p^{2} - k_{o u t}^{2}}

,

y_{b} = {}^{i}a \cdot \sqrt{p^{2} - k_{i n}^{2}}

and

z_{b} = {}^{i}ρ_{s} \cdot \sqrt{p^{2} - k_{o u t}^{2}}

.

\begin{array}{l} S_{n} (a) = I_{n} (x_{b}) \cdot \sum_{m = - \infty}^{\infty} \exp (- j \cdot (n - m) \cdot {}^{i}θ_{s}) \cdot I_{n - m} (z_{s}) \cdot K_{m} (z_{b}) \\ S C_{n} = {(- 1)}^{n} \cdot I_{n} (x_{b}) \cdot \sum_{m = - \infty}^{\infty} \exp (- j \cdot (n - m) \cdot {}^{i i^{'}}θ) \cdot K_{n - m} (z_{i^{'}}) \cdot {}^{i^{'}}B_{m}^{t - 1} (p) \end{array}

\begin{array}{l} S_{n}^{'} (a) = I_{n}^{'} (x_{b}) \cdot \sum_{m = - \infty}^{\infty} \exp (- j \cdot (n - m) \cdot {}^{i}θ_{s}) \cdot I_{n - m} (z_{s}) \cdot K_{m} (z_{b}) \\ S C_{n}^{'} = {(- 1)}^{n} \cdot I_{n}^{'} (x_{b}) \cdot \sum_{m = - \infty}^{\infty} \exp (- j \cdot (n - m) \cdot {}^{i i^{'}}θ) \cdot K_{n - m} (z_{i^{'}}) \cdot {}^{i^{'}}B_{m}^{t - 1} (p) \end{array}

Where,

z_{s} = {}^{i}ρ_{s} \cdot \sqrt{p^{2} - k_{o u t}^{2}}

,

z_{i^{'}} = {}^{i i^{'}}r \cdot \sqrt{p^{2} - k_{o u t}^{2}}

,

{}^{i i^{'}}r

and

{}^{i i^{'}}θ

are, respectively, the radial distance and the azimuth angle of the center of the i^th inhomogeneity with respect to the center of the

i^{' t h}

inhomogeneity.

Therefore, the unknown coefficients ${}^{i}B_{n}^{t} (p)$ and ${}^{i}C_{n}^{t} (p)$ can be determined by a recurrence process. The initial terms, ${}^{i}B_{n}^{1} (p)$ , are given in the study of Walker et al. [20], and therefore not written here. Since $^{i} Φ_{s c a t}^{(l)}$ approaches 0 when l approaches infinity, and since the number of cylindrical inhomogeneities, N, is finite, then:

\lim_{t \to \infty} | \frac{{}^{i}B_{n}^{t} (p)}{{}^{i}B_{n}^{t + 1} (p)} | = 1, i = 1, 2, ..., N

Considering Eq. (12), the stop condition of the recurrence process can be written as:

| \frac{{}^{i}B_{n}^{t} (p)}{{}^{i}B_{n}^{t + 1} (p)} - 1 | \leq ε i = 1, 2, ..., N

ε

is an arbitrarily small quantity that depends on the required precision.

Obviously no boundary conditions are defined for the tissue surface because the medium is infinite in the model. This could be problematic for DOT but not for PA imaging because the PA signal is produced inside absorption inhomogeneities.

3. Numerical phantoms

Our aim was to investigate the validity of this new model when used as a forward model in a quantitative PA imaging reconstruction process. The quantity of interest is the optical fluence inside the cylindrical absorption inhomogeneities, which represent the blood vessels. Indeed the PA signal is directly linked to the absorbed energy into the blood vessels, which is proportional to the product of the absorption coefficient and the optical fluence. We used NIRFAST software [10,11] as the gold standard for the characterization of the optical fluence into cylindrical inhomogeneities.

Several numerical phantoms were defined so as to validate the model. They were all based on a cube ( $4 \times 4 \times 4 c m^{3}$ ), and a point light source was positioned 0.1 cm ( $1 / μ_{s}^{'}$ ) from the border. The distance between nodes was taken as 0.05 cm ( $0.5 / μ_{s}^{'}$ ) for NIRFAST simulations. Optical properties of phantoms were close to the typical parameters of biological tissues [37] in the near-infrared region. The reduced scattering coefficient was constant in the whole medium (background and inhomogeneities), $μ_{s}^{'} = 10 c m^{- 1}$ .

The first three sets of phantoms include two inhomogeneities to investigate the influence of different parameters of inhomogeneities (depth, size, and absorption) on the interaction between inhomogeneities. The last two phantoms include multiple inhomogeneities to validate our model for more complex cases.

For each phantom we simulated the optical fluence distribution with NIRFAST. This solution is denoted “NIRFAST” and is used as the gold standard. We then simulated the optical fluence distribution with the analytical model in each phantom, considered in this case as an infinite cube. $ε$ was set to $10^{- 3}$ , which fixed a value of t_max by use of Eq. (13). We evaluated the influence of the orders of scattering wave by simulating different analytical solutions with increasing values of t. The optical fluence distribution corresponding to, t = 1, was denoted “1-order.” We calculated the optical fluence distribution as the value of t increased from 2 to t_max, these solutions were denoted “t-order”. n in our analytical model arises from the separation of variables and should be infinity in theory. In the following simulation results, the values of n ranges from −20 to 20. The upper limit of n can be fixed by a method similar to Eq. (13). In this paper, it was calculated by following equation,

\frac{| Φ_{i n c} (n) - Φ_{i n c} (n + 5) |}{Φ_{i n c} (n + 5)} < 10^{- 5} .

In order to evaluate quantitatively the error of the proposed model compared with the chosen standard, we calculated the mean relative error (MRE) with the following equation:

M R E = \frac{1}{N} \sum_{i = 1}^{n} | \frac{x_{i} - y_{i}}{y_{i}} |

In this paper, MRE measured the optical fluence distribution difference between the result calculated by the analytical model, x, and the reference result calculated by NIRFAST, y. The MRE is calculated only into the inhomogeneities. The MRE is normalized to the local reference result y. Therefore, this parameter evaluates local errors. It is the more relevant way of evaluating the ability of the model for use in PA imaging. Global parameters would neglect deep parts of the inhomogeneities that receive lower optical fluence, and are then the more challenging part of PA imaging.

4. Results and discussion

4.1 Two inhomogeneities with different depth

This set of simulation examples shows the interaction of two inhomogeneities with different depth. Two columns in Fig. 3 show two phantoms and the corresponding results. The position of inhomogeneities (denoted by c1 and c2) and light source are shown in the top rows, Fig. 3(a) and Fig. 3(d). The corresponding depths of the c1 center are 2.1 cm and 2.3 cm, the depth of the c2 center in Fig. 3(a) and Fig. 3(e) is 1.5 cm. The absorption coefficient of all inhomogeneities is 0.8 cm⁻¹. The radius of all inhomogeneities is 0.3 cm. The corresponding fluence distribution along the line y = 0 cm is shown in the second rows of Fig. 3(b) and Fig. 3(e). The curves of MRE versus the order of scattering are given in Fig. 3(c) and Fig. 3(f). The resolved t_max values are 4 and 3, the difference of t_max indicates that a higher order of scattering needs to be considered if the two inhomogeneities are closer. The t_max values were determined as c1 was moved toward the -x direction along the dotted line in Fig. 3(a). Figure 4 shows the curve of t_max versus the depth of the c1 center. This curve shows that a higher order of scattering needs to be considered if the depth of the c1 center decreases. This is due to the fact that the scattering wave is an outgoing wave, which is reduced as the distance of propagation increases.

Fig. 3 Comparison between the analytical model and NIRFAST for two inhomogeneities with different depth. Top row ((a) and (d)): Horizontal cross section of the first numerical phantom of two cylindrical inhomogeneities with different depth (depths of the c1 center in (a) and (e) are 2.1 cm and 2. 3cm, depths of the c2 center in (a) and (e) are 1.5cm) in a homogeneous background (µ_a = 0.2cm⁻¹ and µ’_s = 10cm⁻¹). Second row ((b) and (e)): Horizontal cross section of the optical fluence distribution through y = 0cm; the dashed lines indicate the positions of the inhomogeneities c1 and c2. Third row ((c) and (g)): MRE corresponding to the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm ].

Download Full Size | PDF

Fig. 4 Curve of t_max versus the depth of the c1 center.

Download Full Size | PDF

4.2 Two inhomogeneities with different size

This set of simulation examples shows the interaction of two inhomogeneities with different radius. The position of the inhomogeneities (denoted by c1 and c2) and light source are shown in Fig. 5(a) and Fig. 5(e). The corresponding radius of c1 is 0.5 cm and 0.7 cm, the radius of c2 in Fig. 5(a) and Fig. 5(e) is 0.3 cm. The absorption coefficient of all inhomogeneities is 0.8 cm⁻¹. The corresponding fluence distribution along the line y = −0.5cm is shown in Fig. 5(b) and Fig. 5(f). The fluence distribution along the line y = 0.5cm is shown in Fig. 5(c) and Fig. 5(g). The resolved t_max values are 3 and 5. The curve of MRE versus the order of scattering is given in Fig. 5(d) and Fig. 5(h). The t_max values were determined as the radius of c1 was changed from 0.1 cm to 0.7 cm. Figure 6 shows the curve of t_max versus the radius of c1. This curve shows that the order of scattering to be considered increases with the increase in the radius of c1. This is due to the fact that the interaction of c1 and c2 is improved as the radius of c1 increases. In this set of simulations the center of two inhomogeneities is fixed because the source of scattering waves can be considered as the center of inhomogeneities according to Eq. (5).

Fig. 5 Comparison between the analytical model and NIRFAST for two inhomogeneities with different radius. Top row ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different radius (radius of c1 in (a) and (e) is 0.5cm and 0.7cm, respectively, radius of c2 in (a) and (e) is 0.3cm) in a homogeneous background (µ_a = 0.2cm⁻¹ and µ’_s = 10cm⁻¹). Second row ((b) and (f)): Horizontal cross section of the optical fluence distribution through y = −0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) and (g)): Horizontal cross section of the optical fluence distribution through y = 0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row ((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm ].

Download Full Size | PDF

Fig. 6 Curve of t_max versus the radius of c1.

Download Full Size | PDF

4.3 Two inhomogeneities with different absorption coefficient

This set of simulation examples shows the interaction of two inhomogeneities with different absorption coefficient. The position of the inhomogeneities (denoted by c1 and c2) and light source are shown in Fig. 7(a) and Fig. 7(e). The corresponding absorption coefficients of c1 are 0.4 cm⁻¹ and 0.8 cm⁻¹, and the absorption coefficient of c2 in Fig. 7(a) and Fig. 7(e) is 0.8 cm⁻¹. The radius of all inhomogeneities is 0.3 cm. The corresponding optical fluence distribution along the line y = −0.5 cm is shown in Fig. 7(b) and Fig. 7(f). The optical fluence distribution along the line y = 0.5 cm is shown in Fig. 7(c) and Fig. 7(g). These curves of optical fluence indicate the higher the absorption coefficient of c1 is, the faster the optical fluence decays. This is due to the fact that the rise of absorption results in an increase of wavenumber. The resolved t_max value is 4. The curve of MRE versus the order of scattering is given in Fig. 7(d) and Fig. 7(h). The t_max values were determined as the absorption of c1 was changed from 0.1 cm⁻¹ to 0.9 cm⁻¹. Figure 8 shows the curve of t_max versus the absorption coefficient of c1. This curve shows that t_max is slightly influenced by the absorption value. We did not note any significant influence of the absorption coefficient properties variation.

Fig. 7 Comparison between the analytical model and NIRFAST for two inhomogeneities with different absorption. Top row ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different absorption coefficient (µ_a of c1 in (a) and (e) is 0.4cm⁻¹ and 0.8 cm⁻¹, respectively, µ_a of c2 in (a) and (e) is 0.8 cm⁻¹) in a homogeneous background (µ_a = 0.2cm⁻¹ and µ’_s = 10cm⁻¹). Second row ((b) (f)): Horizontal cross section of the optical fluence distribution through y = −0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) (g)): Horizontal cross section of the optical fluence distribution through y = 0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm ].

Download Full Size | PDF

Fig. 8 Curve of t_max versus absorption coefficient of c1.

Download Full Size | PDF

4.4 Multiple inhomogeneities

In this set of simulations, the first phantom had five identical cylindrical inhomogeneities (denoted c1 to c5 in Fig. 9(a)) embedded within the cube. They had a radius of 0.3 cm and an absorption coefficient µ_a of 0.8 cm⁻¹. The second phantom had three different cylindrical inhomogeneities embedded within the cube. Their radius were set in the range 0.3 cm to 0.5 cm and their absorption coefficient in the range 0.4 cm⁻¹ to 0.8 cm⁻¹ (see Table 1). Figure 9(a) and Fig. 10(a) show, respectively, the horizontal cross sections of the first and second phantoms.

Fig. 9 Comparison between the analytical model and NIRFAST. (a)Horizontal cross section of the first numerical phantom with five cylindrical inhomogeneities (µ_a = 0.8cm⁻¹ and µ’_s = 10cm⁻¹) in a homogeneous background (µ_a = 0.2cm⁻¹ and µ’_s = 10cm⁻¹). (b) Horizontal cross section of the optical fluence distribution through y = −0.7cm; the dashed lines indicate the positions of the inhomogeneities c1 and c3. (c) Horizontal cross section of the optical fluence distribution through y = 0cm; the dashed lines indicate the position of the inhomogeneity c5. (d) MRE corresponding to the area [-1.2cm, 1.2cm] × [-1.2cm, 1.2cm].

Download Full Size | PDF

Table 1. Parameters of the second numerical phantom

View Table

Fig. 10 Comparison between the analytical model and NIRFAST. (a) Horizontal cross section of the second numerical phantom with three cylindrical inhomogeneities (see Table 1 for optical properties) in a homogeneous background (µ_a = 0.2cm⁻¹ and µ’_s = 10cm⁻¹). (b) Horizontal cross section of the optical fluence distribution through y = 0.7cm; the dashed lines indicate the positions of the inhomogeneities c3. (c) Horizontal cross section of the optical fluence distribution through y = 0cm; the dashed lines indicate the positions of the inhomogeneities c1. (d) Horizontal cross section of the optical fluence distribution through y = −0.7cm; the dashed lines indicate the positions of the inhomogeneities c3. (e) MRE corresponding to the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm ].

Download Full Size | PDF

Figure 9(a) shows the horizontal cross sections of the first phantom. The t_max value was found to be 5. The configuration of the inhomogeneities was symmetrical about the line, y = 0. Figure 9(b) and Fig. 9(c) show the optical fluence distribution through the horizontal cross sections corresponding to lines y = −0.7cm and y = 0cm, respectively. We can clearly see that the analytical solutions obtained by only considering first order and first two orders of scattering waves had large errors. This means that the higher-order (larger than second order) scattering waves need to be considered in this case.

The MRE corresponding to the fluence distribution in the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm ] is given in Fig. 9(d). The two remaining inhomogeneities had the same behavior because of the symmetry of the phantom. The curve of MRE versus the order of scattering used in the simulation confirms that a higher order of scattering has to be considered.

In the second phantom we considered a geometry including inhomogeneities with different parameters. The t_max value was found to be 5. Figure 10(b), Fig. 10(c), and Fig. 10(d) show the optical fluence distribution through the lines y = 0.7 cm, y = 0 cm, and y = −0.7 cm, respectively. The analytical solutions obtained by only considering first order and first two orders of scattering waves still had larger errors. Simulation results showed that the analytical solution inside and outside the inhomogeneities was in good agreement with the NIRFAST solution when t = t_max. The curves of the MRE corresponding to the area [-1.2cm, 1.2cm] × [-1.2cm, 1.2cm] versus the number of order, Fig. 10(e), confirm the observation.

These results also demonstrate the relevance of the new analytical model for the investigation of real vasculature, which presents numerous blood vessels with different sizes, depths, and absorption coefficients. Indeed oxy-hemoglobin and deoxy-hemoglobin have very different absorption coefficients, and thus the oxygen saturation of blood vessels greatly influences its absorption coefficient. From the MRE curves we can see that the increasing number of the order of scattering considered, t, optimizes the analytical solution. MRE falls sharply when t goes from 1 to 2, and it falls more slightly when t goes from 2 to t_max. This means that the scattered waves are weakened with increasing orders. As we mentioned previously, this property enables the convergence of the recurrence process.

To compare our analytical model in terms of computation speed, we implemented the algorithms with Matlab (R2011b) and they were run on an ordinary computer (Intel(R) core(TM) i7-2760QM CPU @ 2.40GHz). It took 642 s to run the first example by using NIRFAST software compared with 54 s with the analytical model. This is a 12-fold improvement in computation speed, which means that the analytical model is more likely to be used for real-time imaging when implemented in a faster computing environment.

In all the results, the fluence was given within the interval [-1.2 cm, 1.2 cm] not within the interval [-2 cm, 2 cm], and the fluence distribution was normalized by its maximum value. In the area close to the border, there were larger errors due to different boundary conditions.

The index of refraction is considered as a constant in the whole media. If one wants to consider the variation of the index of refraction, the boundary conditions at the surface of a certain cylinder need to be modified to incorporate Fresnel reflection [20]. If the point source is replaced with a light beam, the source term can be considered as being composed of many point sources, which means that the analytical solution of a beam source can be obtained by an integral of the solution of the point source. It is noteworthy that our model can also be used to research the influence of scattering coefficient on optical fluence distribution, though the scattering coefficient in all numerical phantoms used in simulations has been assumed as a constant in whole medium.

4. Conclusion

In this paper, we proposed an analytical model of the optical fluence distribution for multiple cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium based on new boundary conditions taking into account the scattering waves of the inhomogeneities. The scattering waves were determined by combining a recurrence process with a proposed stop condition. The model was compared with the solution of NIRFAST software for numerical phantoms. The close agreement between the two methods used to simulate the optical fluence distribution into absorption inhomogeneities permits the use of our model as a forward model needed in quantitative PA imaging.

Acknowledgments

This work was performed within the framework of the LABEX CELYA (ANR-10-LABX-0060) and the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program “Investissements d'Avenir” (ANR-11-IDEX-0007) run by the French National Research Agency (ANR).

References and links

1. A. Corlu, R. Choe, T. Durduran, M. A. Rosen, M. Schweiger, S. R. Arridge, M. D. Schnall, and A. G. Yodh, “Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans,” Opt. Express 15(11), 6696–6716 (2007). [CrossRef] [PubMed]

2. V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U.S.A. 97(6), 2767–2772 (2000). [CrossRef] [PubMed]

3. D. A. Boas, T. Gaudette, G. Strangman, X. Cheng, J. J. Marota, and J. B. Mandeville, “The accuracy of near infrared spectroscopy and imaging during focal changes in cerebral hemodynamics,” Neuroimage 13(1), 76–90 (2001). [CrossRef] [PubMed]

4. B. T. Cox, J. G. Laufer, and P. C. Beard, “Quantitative Photoacoustic Image Reconstruction using Fluence Dependent Chromophores,” Biomed. Opt. Express 1(1), 201–208 (2010). [CrossRef] [PubMed]

5. B. Cox, J. G. Laufer, S. R. Arridge, and P. C. Beard, “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt. 17(6), 061202 (2012). [CrossRef] [PubMed]

6. Z. Yuan and H. Jiang, “Quantitative photoacoustic tomography: Recovery of optical absorption coefficient maps of heterogeneous media,” Appl. Phys. Lett. 88(23), 231101 (2006). [CrossRef]

7. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [CrossRef] [PubMed]

8. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef] [PubMed]

9. D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002). [CrossRef] [PubMed]

10. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef] [PubMed]

11. M. Jermyn, H. Ghadyani, M. A. Mastanduno, W. Turner, S. C. Davis, H. Dehghani, and B. W. Pogue, “Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography,” J. Biomed. Opt. 18(8), 086007 (2013). [CrossRef] [PubMed]

12. T. Tarvainen, M. Vauhkonen, and S. R. Arridge, “Gauss–Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transf. 109(17-18), 2767–2778 (2008). [CrossRef]

13. Z. G. Wang, L. Z. Sun, L. L. Fajardo, and G. Wang, “Modeling and reconstruction of diffuse optical tomography using adjoint method,” Commun. Numer. Methods Eng. 25(6), 657–665 (2009). [CrossRef]

14. L. Yao and H. Jiang, “Finite-element-based photoacoustic tomography in time domain,” J. Opt. A, Pure Appl. Opt. 11(8), 085301 (2009). [CrossRef]

15. Z. Yuan and H. Jiang, “Quantitative photoacoustic tomography,” Philos. Trans. R. Soc. Lond. A 367, 3043–3054 (2009).

16. Y. Zhai and S. A. Cummer, “Fast tomographic reconstruction strategy for diffuse optical tomography,” Opt. Express 17(7), 5285–5297 (2009). [CrossRef] [PubMed]

17. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 28(12), 2331–2336 (1989). [CrossRef] [PubMed]

18. A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14(1), 246–254 (1997). [CrossRef] [PubMed]

19. D. A. Boas, M. A. O’Leary, B. Chance, and A. G. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications,” Proc. Natl. Acad. Sci. U.S.A. 91(11), 4887–4891 (1994). [CrossRef] [PubMed]

20. S. A. Walker, D. A. Boas, and E. Gratton, “Photon density waves scattered from cylindrical inhomogeneities: theory and experiments,” Appl. Opt. 37(10), 1935–1944 (1998). [CrossRef] [PubMed]

21. H. O. Di Rocco, D. I. Iriarte, M. Lester, J. Pomarico, and H. F. Ranea-Sandoval, “CW laser transilluminance in turbid media with cylindrical inclusions,” Opt. Int. J. Light Electron Opt. 122(7), 577–581 (2011). [CrossRef]

22. J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchhoff approximation for diffusive waves,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 64(5), 051917 (2001). [CrossRef] [PubMed]

23. J. Ripoll, V. Ntziachristos, and E. N. Economou, “Experimental demonstration of a fast analytical method for modeling photon propagation in diffusive media with arbitrary geometry,” Proc. SPIE 4431, 233–239 (2001). [CrossRef]

24. X.-P. Wu, J.-M. Shi, and J.-C. Wang, “Multiple Scattering by Parallel Plasma Cylinders,” IEEE Trans. Plasma Sci. 42(1), 13–19 (2014). [CrossRef]

25. S.-C. Lee, “Scattering of evanescent wave by multiple parallel infinite cylinders near a surface,” IEEE Trans. Plasma Sci. 147, 252–260 (2014).

26. M. I. Mishchenko, L. D. Travis, and A. Macke, “Scattering of light by polydisperse, randomly oriented, finite circular cylinders,” Appl. Opt. 35(24), 4927–4940 (1996). [CrossRef] [PubMed]

27. D. Felbacq, G. Tayeb, and D. Maystre, “Scattering by a random set of parallel cylinders,” J. Opt. Soc. Am. A 11(9), 2526–2538 (1994). [CrossRef]

28. V. M. Twersky, “Multiple Scattering of Radiation by an Arbitrary Configuration of Parallel Cylinders,” J. Acoust. Soc. Am. 24(1), 42–46 (1952). [CrossRef]

29. C. M. Linton and P. A. Martin, “Multiple scattering by random configurations of circular cylinders: Second-order corrections for the effective wavenumber,” J. Acoust. Soc. Am. 117(6), 3413–3423 (2005). [CrossRef] [PubMed]

30. J. W. Young and J. C. Bertrand, “Multiple scattering by two cylinders,” J. Acoust. Soc. Am. 77, 1190–1195 (1976).

31. W. H. Lin and A. C. Raptis, “Sound scattering by a group of oscillatory cylinders,” J. Acoust. Soc. Am. 77(1), 15–28 (1985). [CrossRef]

32. X. Cheng and D. Boas, “Diffuse optical reflection tomography using continuous wave illumination,” Opt. Express 3(3), 118–123 (1998). [CrossRef] [PubMed]

33. A. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48(3), 34–40 (1995). [CrossRef]

34. A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. 43(5), 1285–1302 (1998). [CrossRef] [PubMed]

35. J. D. Jackson, Classical Electrodynamics, (Wiley, 1975), Chap. 3.

36. X. D. Li, M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Fluorescent diffuse photon density waves in homogeneous and heterogeneous turbid media: analytic solutions and applications,” Appl. Opt. 35(19), 3746–3758 (1996). [CrossRef] [PubMed]

37. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77(4), 041101 (2006). [CrossRef]

Analytical model of optical fluence inside multiple cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium for quantitative photoacoustic imaging

Abstract

1. Introduction

2. Theory

2.1 Incident waves

2.2 Scattered waves

2.3 Boundary conditions

3. Numerical phantoms

4. Results and discussion

4.1 Two inhomogeneities with different depth

4.2 Two inhomogeneities with different size

4.3 Two inhomogeneities with different absorption coefficient

4.4 Multiple inhomogeneities

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (10)

Tables (1)

Equations (18)

Optics Express

Inhomogeneity	Position (x, y)	Radius (cm)	Absorption (cm^-¹)
c1	(0.5, 0)	0.5	0.8
c2	(−0.5, 0.7)	0.4	0.6
c3	(0, −0.7)	0.3	0.4