Multigrid-based reconstruction algorithm for quantitative photoacoustic tomography

Shengfu Li; Bruno Montcel; Zhen Yuan; Wanyu Liu; Didier Vray

doi:10.1364/BOE.6.002424

1. Introduction

Quantitative photoacoustic (PA) imaging consists in producing the spatial distribution of the optical properties of tissue, including absorption coefficient $μ_{a}$ and reduced scattering coefficient $μ_{s}^{'}$ , from a measured absorbed energy distribution $H_{m} (r)$ that can be obtained using conventional PA tomography [1–3]. The process can be written as:

({\hat{μ}}_{a}, {\hat{μ}}_{s}^{'}) = \underset{(μ_{a}, μ_{s})}{\arg \min} ({‖ H_{m} (r) - C \cdot H (μ_{a}, μ_{s}^{'}, r) ‖}^{2})

Here,

H (μ_{a}, μ_{s}^{'}, r)

is the forward model of absorbed energy, C is a calibration factor related to the acquisition system, r refers to an arbitrary point in the medium,

{\hat{μ}}_{a}

and

{\hat{μ}}_{s}^{'}

are the recovered absorption and reduced scattering coefficient, respectively.

There are two key points for quantitative PA reconstruction: (1) the forward model of absorbed energy to obtain $H (μ_{a}, μ_{s}^{'}, r)$ and (2) the optimization scheme to minimize the errors between the absorbed energy measured and the output of the forward model of absorbed energy. Usually, $H (μ_{a}, μ_{s}^{'}, r)$ is obtained by solving the diffusion equation. Analytical solutions exist for simple cases [4–6], and numerical solutions are used for more realistic cases, using for example the finite difference method [7–9] or the finite element method [10–12]. The main optimization schemes used in the bibliography embrace fixed-point iteration and gradient-based optimization [13–15]. All these inversion works were implemented on a fixed fine grid. Choosing a fixed fine grid can reduce the errors of the forward model numerical solution and enhance the resolution. However, resolving the forward model and the inverse problem at a fine resolution leads to high computational complexity and slow convergence.

In this paper, we propose a framework of multigrid inversion for quantitative PA tomography reconstruction in which both the forward model of absorbed energy and the optimization scheme are expressed at multiple resolutions, hereafter called the “pure multigrid inversion algorithm.” Multigrid methods have been applied successfully to solve partial differential equations (PDEs) [16–18]. The obvious advantages of solving the diffusion equation using multigrid methods can be found and understood in associated references [16–18]. The essential aspect of the multigrid methods is to solve a given problem using many discretization scales, which provides rapid convergence rates. Relatively little work has been done on multigrid inverse schemes. Dreyer et al. applied multigrid methods to optimization [19] to solve linear-like problems. However, the inverse problem in quantitative PA reconstruction is a nonlinear problem and therefore this work cannot be extended to solving the present problem. Ye et al. have proposed a multigrid algorithm for optical diffusion tomography [20]. Their algorithm is effective in terms of convergence speed and reliability. However, they did not use multigrid ideas directly and instead made a linear approximation of a nonlinear function using Newton’s method and then using multigrid ideas to solve the linear system.

In the present pure multigrid inversion algorithm, both the forward and inverse problems are solved at different grid resolutions. The computation time is substantially lowered by reducing the required number of iterations and evaluating the forward model and optimization scheme at multiple resolutions. In this paper, it is assumed that the reduced scattering coefficient in the cost function has been known, as proposed in other studies [11,21,22]; only the absorption coefficient needs to be reconstructed and then the nonuniqueness problem existing in a single measurement can be overcome [22,23]. The reconstruction process of this paper can be written as:

{\hat{μ}}_{a} = \underset{μ_{a}}{\arg \min} ({‖ H_{m} (r) - C \cdot H (μ_{a}, r) ‖}^{2})

2. Methods and materials

As mentioned in the introduction, the forward model of absorbed energy, coupled with an optimization scheme, can be used to reconstruct optical parameters. In this section, the two key points will be stated. The multigrid solver of the diffusion equation is used as a forward model, and the multigrid idea is used for optimization.

2.1 Forward model and optimization method

A. Outline of multigrid methods

Multigrid methods are usually used as solvers for linear equations, like Eq. (1), usually representing a discretization form of a differential equation.

A x = b

The multigrid principle can be found in associated references [16–18]. In this scheme, the problem is solved in multiple grids with different spacing. Let x⁽⁰⁾ denote the solution of Eq. (1) at the finest grid and let x^(q) be the solution at an arbitrary coarser grid, here q is a positive integer and (q) represents the number of grid levels . Grid spacing corresponding to the (q)^th grid is 2^q times the spacing of the finest grid. The map of the related quantity from (q)^th to (q + 1)^th (from (q + 1)^th to (q)^th) is based on the linear decimation matrix, noted

I_{(q)}^{(q + 1)}

(linear interpolation matrix, noted

I_{(q + 1)}^{(q)}

).

The multigrid solution of Eq. (1) can be understood as multiple recursive calls to the two-grid algorithm. Therefore, the basic principle of the two-grid algorithm is introduced first for (q)^th and (q + 1)^th grids and then the recursive procedure to be used in this paper can be explained.

Suppose an approximation of a solution for (q)^th grid ${\hat{x}}^{(q)}$ has been found for Eq. (2).

A^{(q)} x^{(q)} = b^{(q)}

Then the residual corresponding to ${\hat{x}}^{(q)}$ can be expressed as $r^{(q)} = b^{(q)} - A^{(q)} {\hat{x}}^{(q)}$ . ${\hat{x}}^{(q)}$ is corrected by the defect correction, $e^{(q)}$ , which can be obtained by solving Eq. (3) at the (q + 1)^th grid and using the linear interpolation operation, $e^{(q)} = I_{(q + 1)}^{(q)} e^{(q + 1)}$ .

A^{(q + 1)} (e^{(q + 1)}) = I_{(q)}^{(q + 1)} r^{(q)}

Thus ${\hat{x}}^{(q)}$ is improved by Eq. (4).

{\hat{x}}^{(q)} \leftarrow {\hat{x}}^{(q)} + e^{(q)}

In this paper, we use the recursion known as the V-cycle multigrid [17], as shown in Fig. 1, in which a four-level grid is used as an example to illustrate an iteration procedure. To implement an iteration, the V-cycle multigrid algorithm moves from fine to coarse grids (from grid (0) to grid (3)) by recursive calls to Eq. (2) and Eq. (3), and then backtracks to a coarse grid (from grid (3) to grid (0)) by recursive calls to Eq. (4).

Fig. 1 An iteration procedure of the V-cycle. The grid levels are noted (0) to (3). The red line represents the path that the multigrid algorithm traverses. The dotted lines represent the grid with different spacing.

Download Full Size | PDF

Figure 1 shows an iteration procedure. In fact, multiple iterations that move back and forth from fine to coarse grids are required to obtain a solution with a certain precision, as shown in Fig. 2, in which the number of updates is noted by 1, 2,…. The maximum number of updates is usually determined by the precision required. It is noteworthy that a normal scheme obtains solutions only using (0)^th grid in Fig. 2.

Fig. 2 Multiple iterations to obtain a solution. The grid levels are noted (0) to (3). The red line represents the multigrid algorithm recursively moving back and forth. The dotted lines represent the grid with different spacing. The position of each update of x⁽⁰⁾ is noted by black dots along the green line with an arrow. The number of updates is noted by 1,2,….

Download Full Size | PDF

B. Forward model

The most general PDE-based model for light transport in a turbid medium is the radiative transfer equation (RTE) [24,25]. However, to solve the RTE, one needs to consider the discretization in both spatial and angular spaces. Also, by assuming that reduced scattering coefficients are much larger than absorption coefficients, the RTE can be reduced to the diffusion equation [14,26], in which one needs to consider only the discretization in the spatial space. Fortunately, in most human tissues, this assumption is valid [5].

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 [14]. The time-independent diffusion equation is written as [11,27,28]:

\nabla \cdot D (r) \nabla Φ (r) - μ_{a} (r) Φ (r) = - S (r)

where

Φ (r)

is the optical fluence at a point in medium r,

μ_{s}^{'} (r)

is the reduced scattering coefficient,

μ_{a} (r)

is the absorption coefficient,

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

is the diffusion coefficient and

S (r)

is the source term.

To solve the numerical solution of Eq. (5), an appropriate boundary condition is needed. Here, the Dirichlet boundary condition is adopted by letting the fluence be equal to zero on an extrapolated boundary (at a distance 2 × D from the medium boundary, here D is the diffusion coefficient). Once the solution of Eq. (5) $Φ (r)$ , is obtained, the absorbed energy can be determined by $H (μ_{a}, r) = μ_{a} (r) Φ (r)$ . In this paper, Eq. (5) is firstly discretized using the finite difference method as references [7], and then $Φ (r)$ is solved by combining the multigrid method in section A with Gauss-Seidel method [18].

C. Multigrid optimization

For the sake of convenience and considering that r refers to a point of the medium, the model of absorbed energy is noted $H (μ_{a})$ and the measured absorbed energy is noted $H_{m}$ . Then the essence of reconstruction solves the following equation: $H_{m} = C \cdot H (μ_{a}) = C \cdot μ_{a} Φ (μ_{a})$ .

Let $μ_{a}^{(0)}$ denote the solution of the inverse problem at the finest grid and let $μ_{a}^{(q)}$ be the solution at an arbitrary coarser grid, here q is a positive integer. Grid spacing corresponding to the (q)^th grid is 2^q times the spacing of the finest grid.

An approximation of the solution for the (q)^th grid ${\hat{μ}}_{a}^{(q)}$ can be found by use of Eq. (6) after several iterations.

{}^{(n + 1)}μ_{a}^{(q)} = \frac{H_{(m)}^{(q)}}{C \cdot Φ^{(q)} ({}^{(n)}μ_{a}^{(q)})}

Here, n is the number of iterations. Note that $Φ^{(q)}$ is to be solved by multigrid methods as mentioned in section A and B. The residual corresponding to ${\hat{μ}}_{a}^{(q)}$ can be expressed as, $r^{(q)} = H_{m}^{(q)} - C \cdot H^{(q)} ({\hat{μ}}_{a}^{(q)})$ . According to multigrid idea, ${\hat{μ}}_{a}^{(q)}$ is improved or corrected using the solution at the (q + 1)^th grid in order to achieve a faster convergence.

Let $μ_{a}^{(q + 1)} + e^{(q + 1)}$ be the exact solution at the (q + 1)^th grid; then one can define the residual for (q + 1)^th grid as,

C \cdot (H^{(q + 1)} (μ_{a}^{(q + 1)} + e^{(q + 1)}) - H^{(q + 1)} (μ_{a}^{(q + 1)})) = r^{(q + 1)}

If $r^{(q + 1)}$ is chosen to be the decimation of the residual of the (q)^th grid, $r^{(q + 1)} = I_{(q)}^{(q + 1)} r^{(q)}$ , then $r^{(q + 1)} = I_{(q)}^{(q + 1)} (H_{m}^{(q)} - C \cdot H^{(q)} ({\hat{μ}}_{a}^{(q)}))$ , and one can obtain Eq. (8) from Eq. (7).

C \cdot H^{(q + 1)} (I_{(q)}^{(q + 1)} {\hat{μ}}_{a}^{(q)} + e^{(q + 1)}) = C \cdot H^{(q + 1)} (I_{(q)}^{(q + 1)} {\hat{μ}}_{a}^{(q)}) + I_{(q)}^{(q + 1)} (H_{(m)}^{(q)} - C \cdot H^{(q)} ({\hat{μ}}_{a}^{(q)}))

Letting

μ_{a}^{(q + 1)} = I_{(q)}^{(q + 1)} {\hat{μ}}_{a}^{(q)} + e^{(q + 1)}

one can obtain Eq. (10) from Eq. (8) and Eq. (9),

{}^{(n + 1)}μ_{a}^{(q + 1)} = \frac{C \cdot H^{(q + 1)} (I_{(q)}^{(q + 1)} {\hat{μ}}_{a}^{(q)}) + I_{(q)}^{(q + 1)} (H_{(m)}^{(q)} - C \cdot H^{(q)} ({\hat{μ}}_{a}^{(q)}))}{C \cdot Φ^{(q + 1)} ({}^{(n)}μ_{a}^{(q + 1)})}

Here, n is the number of iterations. It is noteworthy that

Φ^{(q + 1)}

,

H^{(q + 1)}

and

H^{(q)}

are to be solved by multigrid methods as mentioned in section A and B. The solution for

μ_{a}^{(q + 1)}

, noted by

{\hat{μ}}_{a}^{(q + 1)}

, can be found using Eq. (10) after several iterations. Then, according to Eq. (9) and multigrid idea, one can obtain

e^{(q + 1)} = {\hat{μ}}_{a}^{(q + 1)} - I_{(q)}^{(q + 1)} {\hat{μ}}_{a}^{(q)}

and

{\hat{μ}}_{a}^{(q)} \leftarrow {\hat{μ}}_{a}^{(q)} + I_{(q + 1)}^{(q)} ({\hat{μ}}_{a}^{(q + 1)} - I_{(q)}^{(q + 1)} {\hat{μ}}_{a}^{(q)})

Multigrid optimization can be implemented by recursively applying Eq. (6), Eq. (10) and Eq. (11). The recursive process was explained in section A.

In the following simulations and experiments, Eq. (12) is used to measure the residual of the recovered solution.

E = \frac{‖ H_{m} - C \cdot H ({\hat{μ}}_{a}) ‖}{‖ H_{m} ‖}

2.2 Simulation and experimental tests

The proposed algorithm described above was first validated with a 2D numerical phantom. The geometry is shown in Fig. 3(a) where three circular inhomogeneities (dark red discs) were embedded in an otherwise homogeneous medium (5 × 5 cm²) and a wide light beam (red arrows) was incident on the medium from the top forming an illumination pattern 5 cm in diameter (pink disc). This geometry was used to simulate a horizontal slice with a depth greater than $1 / μ_{s}^{'}$ from the top of phantom in which three parallel cylindrical inhomogeneities are embedded in an otherwise homogeneous medium. Figure 3(b) shows the relative positions of inhomogeneities that are noted by c1, c2 and c3. The parameters of the inhomogeneities are shown in Table 1. The reduced scattering and absorption coefficient were 10 cm⁻¹ and 0.1cm⁻¹ in the background, respectively. Figure 3(c) shows the absorbed energy distribution that was produced by NIRFAST [29,30] for steady-state case, the distance between nodes was taken as 0.05 cm. In practice, one can obtain the absorbed energy from a conventional PA reconstruction. Considering the position of absorbers, the absorbed energy distribution corresponding to the area [−1 cm, 1 cm ] × [−1 cm, 1 cm ] was used as testing data, as shown in Fig. 3(c). c3 has the same absorption coefficient as the background. Thus, c3 is not clear in Fig. 3(c).

Fig. 3 The geometry and simulated absorbed energy. (a) The geometry. (b) Relative positions of inhomogeneities, c1, c2 and c3. (c) The absorbed energy distribution produced by NIRFAST.

Download Full Size | PDF

Table 1. Parameters of the numerical phantom

View Table

The number of grid levels was taken as 2 (the number of grid levels will be discussed below), and the absorption coefficient was recovered with the multigrid scheme and the fixed grid method using a constant absorption coefficient, 0.05 cm^-¹. The recovered results were compared, as shown in Fig. 4, where Fig. 4(A) shows the recovered reconstruction map. Close agreement between the true and recovered absorption coefficient in Fig. 4(B) and in Fig. 4(C) indicates that the multigrid algorithm can be used to reconstruct the absorption coefficient as a normal scheme (based on the fixed grid). Figure 4(E) and Fig. 4(F) show the convergence of the residual, which has been defined in Eq. (12) for the multigrid inversion algorithm and the fixed grid algorithm, as a function of the number of iterations and the CPU time, respectively. The multigrid algorithm reduced the number of required iterations substantially. In this simulation example, the number of iterations for the multigrid approach was approximately half the number for the fixed grid algorithm. The multigrid algorithm was ~1.7 times faster than the fixed algorithm. In this simulation, Eq. (13) was used as a stop condition of the iteration process in the multigrid and fixed grid algorithm. Here, $E^{n}$ is the residual of the n^th iteration.

Fig. 4 Simulation tests of quantitative reconstruction of the absorption coefficient using the fixed grid and multigrid algorithm. (A) Recovered absorption coefficient. True (black line) and recovered absorption coefficient profile (red and green lines) plotted along the lines x = −0.5 cm (B), x = 0.25 cm (C) and x = 0.5 cm (D). The curves of relative errors versus the number of iterations (E) and CPU time (F) for the multigrid and fixed grid algorithms.

Download Full Size | PDF

\frac{‖ E^{n} - E^{n - 1} ‖}{‖ E^{n - 1} ‖} < 0.1

The reconstruction algorithm described above was also tested using experimental data. The corresponding experimental system has been given in previous papers [13]. Briefly, a pulsed light beam emitted by a Nd:YAG laser (wavelength: 532 nm, pulse duration: 6 ns) was sent to the phantom. The light beam was incident on the phantom from the top, as shown in Fig. 5. The incident laser beam diameter was 5.0 cm at the phantom surface. An ultrasound transducer (central frequency, 1 MHz) and the phantom were immersed in a water tank. The transducer was rotated relative to the center of the tank. The right part of Fig. 5 shows the top-view photograph of the phantom where three cylindrical inclusions were embedded in an otherwise homogeneous medium. The background phantom had $μ_{a} = 0.1 c m^{- 1}$ and $μ_{s}^{'} = 10 c m^{- 1}$ while the top left target had $μ_{a} = 0.2 c m^{- 1}$ and $μ_{s}^{'} = 10 c m^{- 1}$ , the top right target had $μ_{a} = 0.2 c m^{- 1}$ and $μ_{s}^{'} = 20 c m^{- 1}$ , and the bottom right target had $μ_{a} = 0.1 c m^{- 1}$ and $μ_{s}^{'} = 20 c m^{- 1}$ .

Fig. 5 Schematic of experimental setup (left) and top-view photograph of the phantom (right).

Download Full Size | PDF

3. Experimental results and discussion

Figure 6 shows the absorbed energy recovered from the acquired signals using the time reversal algorithm [31,32]. The calibration factor (C) was determined by dividing the mean value of H_m by the mean value of H in the area [−1 cm, 1 cm ] × [−1 cm, 1 cm ], as shown in Fig. 6. H was obtained by calculating the absorbed energy with the forward models used in the optimization scheme for the experimental phantom. Figure 7(B) and Fig. 7(C) state that the multigrid scheme can be used to reconstruct the absorption coefficient as the normal scheme (based on a fixed grid). The number of grid levels was also taken as 2 (the number of grid levels will be discussed below). Figure 7(E) and Fig. 7(F) show the convergence of the residual for the multigrid inversion algorithm and the fixed grid algorithm, as a function of the number of iterations and CPU time, respectively. The multigrid algorithm also substantially reduced the number of iterations required, the multigrid algorithm took approximately half the iterations as the fixed grid algorithm under adopting two grid levels. The per-iteration times of the multigrid method were longer than the fixed grid algorithm: the faster convergence of the multigrid algorithm results from the many fewer iterations required. Finally, the multigrid algorithm was approximately 1.6 times faster than the fixed algorithm to meet the condition expressed in Eq. (13).

Fig. 6 Recovered absorbed energy distribution corresponding to the area [−1 cm, 1 cm ] × [−1 cm, 1 cm ].

Download Full Size | PDF

Fig. 7 Experimental validation of quantitative reconstruction of the absorption coefficient using the fixed grid and multigrid algorithm. (A) Recovered absorption coefficient. True (black line) and recovered absorption coefficient profile (red and green lines) plotted along the lines x = −0.5 cm (B), x = 0.25 cm (C) and x = 0.5 cm (D). The curves of relative errors versus the number of iterations (E) and CPU time (F) for the multigrid and fixed grid algorithm.

Download Full Size | PDF

Theoretically, the greater the number of grid levels is, the faster the convergence is, as shown in Fig. 8(A). However, the limit of grid levels is determined by the finest grid. These curves shown in Fig. 8 were obtained with the data used in the simulation test. For experimental data, one can obtain a set of similar curves. Figure 8(B) shows six grid levels and took three times less time than the fixed algorithm to meet the condition expressed in Eq. (13). Besides the number of grid levels, the recursive form is another key point for multigrid methods. In this paper, we used the V-cycle form, as shown in Fig. 1, one can use other recursive forms, such as W-cycle [18].

Fig. 8 Comparison of the convergence speed of multigrid-based algorithms with different numbers of grid levels.

Download Full Size | PDF

Both the simulation and experimental results show that the proposed algorithm can be used to recover the absorption coefficient distribution. Although the recovered results shown in the experiments and simulations were based on the data produced by a single optical wavelength, the extension of the method to multiple wavelengths may allow access to imaging physiological parameters such as the total hemoglobin concentration and oxygen saturation (SO₂). The authors believe that the advantages of the multigrid algorithm would also be relevant for working with multiple wavelengths.

It is also noteworthy that the proposed method can be easily extended into 3D applications, though the simulation and experimental results are obtained for 2D cases. This is due to the fact that all the equations in this paper can be easily extended into 3D cases.

The last comment is about the absorbed energy, shown in Fig. 6, which was obtained using the time reversal algorithm with a constant sound velocity. For targets embedded in acoustically heterogeneous medium, one should select other algorithms that can take spatially varying sound velocity into account.

In summary, we developed a quantitative photoacoustic tomography reconstruction algorithm in which both the forward model and the inverse scheme are solved by a method based on a multigrid idea. This scheme can improve the efficiency of the algorithm since fewer iterations are required compared with the algorithm based on a single fixed grid. The simulation and experimental results show that the multigrid inversion scheme is more efficient than the normal scheme and can obtain a reliable result that has the same precision as with the fixed grid algorithm.

Acknowledgments

This study was conducted within the framework of the LABEX CELYA (ANR-10-LABX-0060) and the LABEX PRIMES (ANR-11-LABX-0063) of the University of Lyon, within the “Investissements d'Avenir” (ANR-11-IDEX-0007) program of the French National Research Agency (ANR). Shengfu Li’s work was also supported by the China Scholarship Council (CSC). The authors would like to thank the ITN OILTEBIA project funded by the European Community's FP7 (Grant Agreement Number 317526) for partial financial support.

References and links

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

2. G. Bal and K. Ren, “Multi-source quantitative PAT in diffusive regime,” Inverse Probl. 27, 075003 (2011). [CrossRef]

3. M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71(1), 016706 (2005). [CrossRef] [PubMed]

4. S. Li, B. Montcel, W. Liu, and D. Vray, “Analytical model of optical fluence inside multiple cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium for quantitative photoacoustic imaging,” Opt. Express 22(17), 20500–20514 (2014). [CrossRef] [PubMed]

5. 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]

6. 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. USA 91(11), 4887–4891 (1994). [CrossRef] [PubMed]

7. S. L. Jacques and B. W. Pogue, “Tutorial on diffuse light transport,” J. Biomed. Opt. 13(4), 041302 (2008). [CrossRef] [PubMed]

8. T. Tanifuji and M. Hijikata, “Finite difference time domain (FDTD) analysis of optical pulse responses in biological tissues for spectroscopic diffused optical tomography,” IEEE Trans. Med. Imaging 21(2), 181–184 (2002). [CrossRef] [PubMed]

9. 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]

10. D. Razansky and V. Ntziachristos, “Hybrid photoacoustic fluorescence molecular tomography using finite-element-based inversion,” Med. Phys. 34(11), 4293–4301 (2007). [CrossRef] [PubMed]

11. Z. Yuan, Q. Wang, and H. Jiang, “Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method,” Opt. Express 15(26), 18076–18081 (2007). [CrossRef] [PubMed]

12. J. Laufer, D. Delpy, C. Elwell, and P. Beard, “Quantitative spatially resolved measurement of tissue chromophore concentrations using photoacoustic spectroscopy: application to the measurement of blood oxygenation and haemoglobin concentration,” Phys. Med. Biol. 52(1), 141–168 (2007). [CrossRef] [PubMed]

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

14. 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]

15. B. T. Cox, S. R. Arridge, and P. C. Beard, “Gradient-based quantitative photoacoustic image reconstruction for molecular imaging,” Proc. SPIE 6437, 64371T (2007).

16. A. K. Garg, I. Kumar, S. Bansal, and I. Goyal, “Multigrid approach for solving elliptic type partial differential equations,” Int. J. Sci. Res. 3, 473–475 (2014).

17. H. Gao and H. Zhao, “A fast forward solver of radiative transfer equation,” (ftp://ftp.math.ucla.edu/pub/camreport/cam09-94.pdf). [CrossRef]

18. W. L. Briggs and V. E. Henson, “A Multigrid Tutorial,” (http://computation.llnl.gov/casc/people /henson/mgtut/welcome.html).

19. T. Dreyer, B. Maar, and V. Schulz, “Multigrid optimization in applications,” Int. J. Comp. Appl. Math. 120, 67–84 (2000).

20. J. C. Ye, C. A. Bouman, K. J. Webb, S. Member, R. P. Millane, and S. Member, “Nonlinear Multigrid Algorithms for Bayesian Optical Diffusion Tomography,” IEEE Trans. Image Process. 10(6), 909–922 (2001). [CrossRef]

21. 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]

22. B. T. Cox, S. R. Arridge, K. P. Köstli, and P. C. Beard, “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt. 45(8), 1866–1875 (2006). [CrossRef] [PubMed]

23. B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A 26(2), 443–455 (2009). [CrossRef] [PubMed]

24. 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]

25. 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]

26. 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]

27. 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]

28. 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]

29. 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]

30. 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]

31. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15(2), 021314 (2010). [CrossRef] [PubMed]

32. P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf, “Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 75(4), 046706 (2007). [CrossRef] [PubMed]

Inhomogeneity	Position (x,y)	Radius(cm)	Absorption(cm^-¹)	Scattering(cm^-¹)
C1	(−0.5, −0.1)	0.2	0.2	10
C2	(0.25, 0.25)	0.2	0.2	20
C3	(0.5, −0.5)	0.2	0.1	20

Multigrid-based reconstruction algorithm for quantitative photoacoustic tomography

Abstract

1. Introduction

2. Methods and materials

2.1 Forward model and optimization method

A. Outline of multigrid methods

B. Forward model

C. Multigrid optimization

2.2 Simulation and experimental tests

3. Experimental results and discussion

Acknowledgments

References and links

Cited By

Figures (8)

Tables (1)

Equations (15)

Biomedical Optics Express