Multigrid-based reconstruction algorithm for quantitative photoacoustic tomography

: This paper proposes a multigrid inversion framework for quantitative photoacoustic tomography reconstruction. The forward model of optical fluence distribution and the inverse problem are solved at multiple resolutions. A fixed-point iteration scheme is formulated for each resolution and used as a cost function. The simulated and experimental results for quantitative photoacoustic tomography reconstruction show that the proposed multigrid inversion can dramatically reduce the required number of iterations for the optimization process without loss of reliability in the results.


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 ( ) m H r that can be obtained using conventional PA tomography [1][2][3]. The process can be written as: is obtained by solving the diffusion equation. Analytical solutions exist for simple cases [4][5][6], and numerical solutions are used for more realistic cases, using for example the finite difference method [7][8][9] or the finite element method [10][11][12]. The main optimization schemes used in the bibliography embrace fixedpoint iteration and gradient-based optimization [13][14][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][17][18]. The obvious advantages of solving the diffusion equation using multigrid methods can be found and understood in associated references [16][17][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:

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.

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.
The multigrid principle can be found in associated references [16][17][18]. In this scheme, the problem is solved in multiple grids with different spacing. Let x (0) 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 ( ) The multigrid solution of Eq. (1) can be understood as multiple recursive calls to the twogrid 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 ( ) q x has been found for Eq. (2).
(2) Then the residual corresponding to ( ) q x can be expressed as ( ) x is corrected by the defect correction, ( ) q e , which can be obtained by solving Eq. (3) at the (q + 1) th grid and using the linear interpolation operation, ( ) 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).  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.

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 timeindependent diffusion equation is written as [11,27,28]: 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 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 ( ) 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 ( ) q a μ can be expressed as, According to multigrid idea, ( ) q a μ is improved or corrected using the solution at the (q + 1) th grid in order to achieve a faster convergence.
Let ( 1) ( 1) q q a e μ + + + be the exact solution at the (q + 1) th grid; then one can define the residual for (q + 1) th grid as, C H e H r μ μ If ( 1) q r + is chosen to be the decimation of the residual of the (q) th grid, Here, n is the number of iterations. It is noteworthy that (q 1) + Φ , (q 1) H + and (q) H are to be solved by multigrid methods as mentioned in section A and B. The solution for ( 1) q a μ + , noted by ( 1) q a μ + , can be found using Eq. (10) after several iterations. Then, according to Eq. (9) and multigrid idea, one can obtain ( 1) ( 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.

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 2 ) 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 −1 and 0.1cm −1 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).  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 -1 . 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, n E is the residual of the n th iteration. 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   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). 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].

Experimental results and discussion
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 2 ). 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.