Improved pertubation method for absorption map reconstruction in diffusion optical tomography

The perturbation method is an effective approach for optical tomography reconstruction, however its success depends to a great extent on how close initial estimates are to the actual solutions. In addition, the linear perturbation method can only be applied to the reconstruction of differences between two similar states. To overcome these limitations, we present a pre-scaling technique applied to the qualitative reconstruction of the absorption map. In this method, the initial estimate of the absorption coefﬁcient is scaled to an appropriate value before it is employed in iteration. The scaled estimate leads to the forward solution on the boundary close to the measured data. In the simulated experiments, reconstructions were performed with and without the pre-scaling technique. The results demonstrate that the proposed technique extends the selection of the initial estimate of optical properties, and makes it feasible that the absolute value of optical properties can also be used in the linear perturbation approach. [DOI: 10.2971/jeos.2007.07015]


INTRODUCTION
Diffusion optical tomography (DOT) is a relatively new but rapidly growing technique, which has shown its powerful potential in medical imaging [1]- [3].Currently the main applications of DOT are monitoring cerebra oxygenation for neonates and detecting breast cancer.
Tissue is a highly scattering medium and its interaction with photons is highly nonlinear, which makes the reconstruction of DOT images a very challenging task.Many investigators have made significant contributions to this field, and various reconstruction algorithms have been developed [4]- [11].The perturbation method is one of the most effective approaches, however its success depends considerably on how close the initial estimate is to the actual solution [4].When they deviate from each other significantly, the perturbation method fails, so selecting good initial estimates of optical properties generally requires prior information.It is not easy in any case.The linear perturbation method is typically used in the reconstruction of difference imaging where measurements are taken before and after a small change in the optical properties [1].In this case, two sets of data are generally required.If the reconstruction is performed using the absolute value, i.e. one set of acquired data, the full non-linear problem must be considered.
A novel nonlinear perturbation approach was proposed by A. Torricelli et al [12,13].It is based on the application of "Padé approximants" and can empirically extend the linear perturbation method.In this paper, we propose another approach, which also extends the linear perturbation method by empirically selecting the initial estimates of optical parameters.In addition, it makes it feasible that the linear perturbation method can be done using absolute values.This approach is simple but useful to primarily locate a tumor in tissue.Furthermore the requirement of prior information about the optical properties is not necessary.
The difference between the absorption coefficients of diseased and healthy tissue is meaningful in clinical diagnosis.In many practical cases, the calculation of the exact absorption coefficient distribution in tissue is not necessary since the outline of the absorption contrast map between the diseased and the healthy tissue is sufficient for clinical diagnosis.Thus in many situations a qualitative, rather than a quantitative, absorption map is clinically sufficient.
To achieve this objective, we present a pre-scaling technique for the qualitative reconstruction of tissue absorption maps.In this method, we scale the initial estimate of the absorption coefficient to an appropriate value, which consequently leads to forward solutions, namely the derived values on the boundary, that are close to the measured data.The scaled estimate is then iteratively updated as in the conventional linear perturbation method.
A number of different sets of initial estimates were employed in simulated experiments to assess the effectiveness of the proposed method.Further simulations were also run in which noise was added to test the effect on the reconstruction results.

OPTICAL MODEL IN TISSUE
As photons migrate in tissue they undergo multiple scattering as well as absorption.In the near-infrared (NIR) region scatter-ing is the dominant interaction.The process of photon transport in tissue can be represented using the radiation transport equation (RTE) [1], however, solving the RTE is extremely computationally expensive.Under certain assumptions, the P 1 approximation to the RTE derives the Diffusion Equation (DE), which is more widely used to model light transport in tissue [1,4,14].
Optical imaging systems can be broadly divided into three categories [1]: continuous-wave (CW), time-domain (TD) and frequency-domain (FD).Accordingly, the diffusion equation has three versions.The time and frequency-domain versions have been derived in detail in earlier literature [14].In this paper, the CW model is used, therefore the diffusion equation can be written as: where r is the location in the tissue domain Ω, Φ(r) is the photon density distribution, q(r) is the source term, µ a (r) is the absorption coefficient and D is the diffusion coefficient given by: where µ s = (1 − g)µ s is the reduced scattering coefficient, µ s is the scattering coefficient and g is the anisotropic factor.The spatially dependent coefficients D(r) and µ a (r) are the two main optical properties that reflect the function of the diseased and healthy tissue.It is the objective of DOT reconstruction to reconstruct these.If the source q(r) is a collimated incident beam, it is treated as a spatial Dirac function, namely a "point source" under the surface ∂Ω at a depth of one mean free length [6], where q(r) = q 0 δ(r − r s ), r s is the location of the equivalent point source and q 0 is the strength of the source.
The Robin boundary condition in the steady-state case is usually employed [15] and hence the measured photon flux Γ(ξ) on the boundary can be expressed as: where n is the outward normal at the site ξ on the boundary ∂Ω.In the CW case Eq. ( 1), along with boundary condition (2), is the most commonly used forward model for DOT and is the model that will be used in this paper.

PERTUBATION METHOD FOR DOT INVERSE PROBLEM
The definitions of the forward and inverse problems for DOT were presented in the literature [14].The forward problem is often expressed by a general nonlinear forward operator as y = F[p(r)] and the inverse by a general inverse operator as where p is the parameter in domain Ω and y is the measured data set on the boundary ∂Ω.
The forward problem can be solved by stochastic or deterministic approaches [16].The finite-element method (FEM) is a very powerful deterministic approach as it is suitable for arbitrary geometries and highly inhomogeneous parameter distributions [6].Using FEM, the forward problem can be reduced to a linear matrix equation of the form [A]{Φ} = {b}.Standard methods, such as the Cholesky decomposition or conjugated gradient method can then be employed.
The inverse problem can be solved by analytical, backprojection, linear or nonlinear methods [14].The perturbation method is a linearisation method based on a Taylor expansion.
If the estimate p is close to the true solution p, the forward operator y = F[p(r)] can be expanded as a Taylor series around p [4]: where F and F are the first-and second-order Fréchet derivatives of the forward operator F. In the discrete case F is represented by J, the Jacobian matrix, and F is represented by H, the Hessian matrix.Neglecting terms above first-order and letting ∆y = (y − ŷ) and ∆p = (p − p), Eq. ( 3) becomes: The linear perturbation method is based on Eq. ( 4), in which the parameter p is iteratively updated to approach the true solution p. Unfortunately its success is significantly dependent on how close the initial estimate is to the true solution [4,16].
In addition, this method is suitable for relative rather than absolute image reconstruction [1].
To overcome these limitations, we present a pre-scaling technique that is used to qualitatively reconstruct the absorption map.Selection of initial estimates is more flexible when using this method.
The reverse process is usually interpreted as the optimization of an objective function: where (Γ j,i ) mea is the measurement on the boundary ∂Ω at the i th detector due to the j th source, (Γ j,i ) cal is the predicted value at the same site, S and M are the number of sources and detectors respectively.The inverse process is to iteratively adjust the optical parameters of the tissue so that the predicted values (Γ j,i ) cal derived by the forward problem gradually approach the measured data (Γ j,i ) mea .If the initial estimate p deviates from its true solution p considerably, the derived value (Γ j,i ) cal will accordingly deviate considerably from the measurement (Γ j,i ) mea .In this case, the conventional perturbation method can not be applied.The idea of the pre-scaling technique is to adjust the initial estimate p according to the difference between (Γ j,i ) cal and (Γ j,i ) mea so as to make it closer to p before the conventional perturbation method is used.

Theoretical foundation for the pre-scaling technique
When the diffusion equation is employed to model photon transport in tissue and CW measured data are employed in optical tomography, the diffusion and absorption coefficients cannot be recovered simultaneously [17,18].The absorption coefficient µ a is, however, more meaningful than the scattering coefficient µ s in clinic diagnosis, so in many cases only the distribution of the absorption coefficient in tissue is recovered.
Before deriving the pre-scaling technique, let's reveal the relationship between the absorption coefficient and the boundary measurement.In other words we want to make clear how the boundary measurements Γ(ξ) vary when the absorption coefficient µ a varies.Using this knowledge we can then investigate how to scale the initial estimate of μa to an appropriate value.
For simplicity, it is assumed that the tissue domain Ω is filled with an homogeneous object (g = 0) with diffusion coefficient D and absorption coefficient µ a .Furthermore since D = 1/[3(µ a + µ s )]and the diffusion Eq. ( 1) is based on the assumption µ a µ s , it is reasonable to say that D does not depend on µ a .
Fixing the other quantities, such as the source intensity q and the diffusion coefficient D we observe by means of a numerical experiment how the boundary measurement Γ(ξ) varies when the absorption coefficient µ a varies from 0.001 to 0.4.The curve of Γ(ξ) against µ a at any location on the boundary is illustrated in Figure 1a.The result is similar to an exponentially decaying curve.This can be partially interpreted as the fact that energy attenuation during light transport in tissue is governed by exponential decay [19]: Φ(r )/Φ 0 (r) = exp(−µ tr z), where µ tr = µ a + µ s and z is the distance that photons travel from location r to r .Curves that have a similar shape shown in Figure 1a can always be derived for any given diffusion coefficient D and source intensity q when µ a is varied.
As shown in Figure 1b, the measurement Γ(ξ) can be normalised by multiplying by a factor α. Similarly, the absorption coefficient µ a can be normalized by multiplying by β.The areas s 1 = βµ a1 • αΓ 1 (ξ) and s 2 = βµ a2 • αΓ 2 (ξ) can be considered to be roughly equal i.e. s 1 ≈ s 2 thus µ a1 • Γ 1 (ξ) ≈ µ a2 • Γ 2 (ξ).Letting Γ 1 (ξ) = Γ j,i mea , Γ 2 (ξ) = Γ j,i cal , µ a2 = μa and µ a1 = µ a , the following equation can be derived: where [λ j,i ] is defined as the gain factor due to the j th source and the i th dector.Eq. ( 6) is the mathematical basis of the pre-scaling technique.Considering all M detectors and all S sources, the mean gain factor λ mean , which more accurately reflects the difference between the estimate μa and the true so-lution µ a , can be derived: If λ mean > 1, it means that the estimate μa is larger than the true solution µ a and hence μa should be decreased, otherwise μa should be increased.The estimate should thus be scaled using the formula: μa := μa /λ mean .The estimate μa is repeatedly adjusted until λ mean ≈ 1, which means that the predicted value Γ j,i cal on the boundary is roughly equal to the measurement Γ j,i mea , and the estimate μa is close to the true solution µ a .The scaled μa can then be used in the conventional perturbation method to be iteratively updated.

Improved perturbation method for absorption map reconstruction
We present the improved perturbation method for the reconstruction of absorption maps in algorithm 1.To reduce the intrinsic ill-posedness of the inverse problem, the Tikhonov regularization is employed [20,21].

Algorithm 1. Improved perturbation method for the reconstruction of absorption maps
Begin: 1. Set the tolerances ε needed to terminate the iterations and ε λ > 0 needed to pre-scale the initial absorption coefficient, the initial λ mean > 1 + ε λ , and the initial estimates of parameters p(including μa and μs or D).

Compute the mean gain factor λ
Update the absorption coefficient using μa : = μa /λ mean .

End begin
In algorithm 1, the loop consisting of steps 2 and 3 is the prescaling process.Its function is to scale the initial estimate of absorption coefficient μa to the same magnitude of the true solution µ a .Neglecting these two steps, algorithm 1 reduces to the conventional perturbation method.
In the reverse process, the Jacobian construction is the most computationally costly.Various methods to calculate the matrix J, such as the standard, joint source, direct and perturbation methods have been developed [14,15,22,23].In this paper, the joint source method is employed.
In actual applications, the noise model should be considered.Inverse problem are intrinsically ill-posed, which means slight fluctuations in the measured data contaminated by noise might cause the reconstructed result to drastically deviate from its true solution.There are three main kinds of noise in a photodetector system [24]: thermal noise, shot noise, and relative intensity noise.For ideal photodetectors, thermal noise is negligible and shot noise deriving from photodetector current is dominant.At low optical intensities shot noise is governed by a Poisson distribution, however when the intensity is significantly large, it is approximately governed by a Gaussian distribution [11].The noise model for measured data y j,i was developed in detail in literature [25].If each measured data y j,i is independently corrupted by Gaussian noise with standard deviation σ j,i , the objective function ( 5) should be rewritten as follows [4]: The inverse problem then reduces to a least-squares (LS) minimization procedure.Due to the fundamental limitation of linearised algorithms [26] in this situation, full-nonlinear schemes or statistical methods [27], such as Bayesian framework [11,28] and the Monte Carlo method [29], are required.Analysis of the noise model is not the task of this paper.In this article we only simply add Gaussian noise to the simulated measurements to test the effect of noise on the reconstruction results.

SIMULATED EXPERIMENTS AND RESULTS
The geometry of the model used in simulated experiments is illustrated in Figure 2. The circular domain Ω of diameter 4cm is full of homogeneous tissue-like medium, which is isotropic, namely, the anisotropic factor g = 0. Inside the homogeneous medium a circular absorber of diameter 0.2cm is embedded, at the position with coordinate (1.2, 0) (unit: cm) where the coordinate origin is taken at the center of the circular domain Ω.
The source fibers and detector fibers are arranged uniformly around the boundary ∂Ω, as illustrated in Figure 2a.The practical model is illustrated in Figure 2b, in which the domain Ω is divided into 1644 elements jointed at 887 vertex nodes.Some elements are selected as heterogeneities and are painted black in the FEM mesh graphics.
Since the reconstruction of µ a is done with arbitrary selection of D, the reconstructed absorption coefficient distribution is not the exact absorption coefficient map, but a contrast absorption map.Its function is to outline the difference between the homogeneity, such as healthy tissue, and the heterogeneity, such as diseased tissue.
It should be stated that in our work µ a and D no longer represent physical optical properties of tissue, but are merely mathematical quantities used for reconstruction purposes.In conventional perturbation methods, the selection of the initial estimates of µ a and D should be limited to regions near their actual physical value by use of prior information.In this paper, the chosen initial parameters can be completely out of these bounds.
In the simulated experiments we assumed the true parameters to be as follows: µ s = 2 and µ a = 0.025 for the background medium, namely the homogeneous tissue-like medium (consequentlyD = 0.165); µ s = 4 and µ a = 0.5 for the heterogeneity, namely the embedded absorber (D = 0.741) and g=0 for both background medium and embedded absorber.The simulated data were derived from the forward problem of the model shown in Figure 2 and based on the above mentioned parameters.They were noise-free data.Only one set of data that was acquired after the absorber had been embedded into the background medium was used.
The initial estimates were selected as: μa = 0.2 and D = 4.We can see that they are both far from the true background solutions: µ a = 0.025 and D = 0.165.Figure 3 illustrates the reconstruction results using the perturbation method without the pre-scaling technique, i.e., the conventional perturbation method.The figure demonstrates that the conventional perturbation method is unable to obtain the desired results through the above initial estimates.
When the pre-scaling technique was employed, by steps 2 and 3 of algorithm 1, the estimate of μa was adjusted to the same magnitude of the true solution µ a , which is shown in column one (the column with µ a = 0.025 and μa = 0.2) of Table 1.
The initial mean gain factor λ mean , written as λ for short, was selected as 2, and the tolerance ε λ for pre-scaling process was selected as 0.001.As shown in Table 1, the absorption coefficient μa was quickly (only three iterations) scaled from its initial estimate 0.2 to 0.0312, which is of the same magnitude as the true solution 0.025, and then the scaled estimate can be used in the conventional perturbation method to be iteratively updated.The reconstruction results are illustrated in Figure 4.As the iterations increase, the reconstruction image gradually approaches the correct solution.The curve of the objective function is illustrated in Figure 4d.The speed of the convergence depends on the factor α when the Tikhonov Regularization (step 5 of Algorithm 1) is performed.
In the above instance, μa µ a .Now let's investigate the other extreme situation: μa µ a , in which μa = 0.001.The prescaling process of μa is shown in column two of Table 1.It was also scaled to 0.0312.
We changed the absorption coefficient µ a of the background medium to 0.075 and maintained the other quantities.The initial estimate of μa was selected as 0.2 and 0.001 respectively.We can observe that by the pre-scaling technique, they were all scaled to the value 0.0748, which is of the same magnitude as the true solution 0.075.The pre-scaling processes are shown in Table 2.Moreover, we have tried various other models to test the proposed technique.For example we changed the location and size of the heterogeneity, and then various sets of estimates were employed.The above mentioned conclusion can also be reached.
This suggests that the pre-scaling technique is an effective improvement of the conventional linear perturbation method for qualitative reconstruction of the absorption map.
To test the effect of noise on the reconstructed results, we added white Gaussian noise to the simulated data.The data with 1%, 5%, 10% and 15% noise were used for reconstruction respectively.In Figure 5 the reconstructed results after 60 iterations are shown.In this simulation, the parameters were the same as in Figure 4, i.e., the true parameters were: background medium with µ s = 2 and µ a = 0.025 (D = 0.165) and heterogeneity with µ s = 4 and µ a = 0.5(D = 0.741).The initial estimates were also selected as: µ a = 0.2 and D=4.The pre-scaling technique was used.
Observing the results shown in Figure 5, we can see that, the larger the noise, the more gravely the reconstructed embedded absorber deviates from its true solution.

DISCUSSION
Reconstruction of diffuse optical images is a non-unique, illposed and underdetermined problem since the number of the unknown image pixels is much bigger than the number of the known measured data.This significantly influences the image quality, particularly the spatial resolution.In addition, the perturbation method described by Eq. ( 4) is a linear approximation to the origin problem (3), in which the terms above firstorder are neglected.This makes the image quality intrinsically no better than the image reconstructed by full-nonlinear methods, especially when the terms above the first-order in Eq. ( 3) play a relatively important role.Furthermore, the measurement noise is another important factor that corrupts the image quality.One way to improve the resolution of the reconstructed image is to make use of prior information as much as possible, for example, take advantage of the anatomical imaging or the physiology information.Another way to improve the image quality is to reduce the noise contaminating the measured data.Knowing the derivation of the noise and its model is very useful to make the measured signal purer.

FIG. 3
FIG. 3 Reconstruction results without pre-scaling technique, in the case of the initial estimates µ a = 0.2, D = 4.

FIG. 4
FIG. 4 Reconstruction results with pre-scaling technique, in the case of the initial estimates µ a = 0.2,D = 4.

FIG. 5
FIG. 5 Reconstruction results with noise-polluted data, in the case of the initial estimates µ a = 0.2, D = 4, iteration 60