Application of Reconstruction and Optimization Algorithms in Optical Tomography ()
1. Introduction
Medial imagining is recently the most growing field of research. It is a technique and a process of creating image of inner body for analysis and takes action to improve medical disorder through visual representation of some organs and tissues. In past, most reliable way of diagnosing disease was biopsy in which a sample of tissue is taken from body to examine it more closely. Invasive surgery was needed to take tissue sample. This process was painful as well as harmful. Therefore, there is an increasing need for new diagnostic techniques that are capable of absorbing tissue data without any surgery. First technique came to know is X-ray. Later on, there are several medical techniques are discovered such as MRI, CT scan, PET and Nuclear Medicine etc. They all are very powerful techniques but have also harmful effects on body due to high radiations and contrast agents’ absorption on body and cannot be used on regular basis. On the other hand, in some cases like breast cancer diagnosis or brain tumors diagnosis CT scan or ultrasounds are very costly and need more expertise as well as provides relatively poor precision in certain cases. Therefore, scientists aimed to search for an imaging technique that can non-invasively differentiate between malignant and healthy tissues [1] - [5].
Optical tomography is the best alternative to all of these techniques. It is a non-invasive technique that uses visible or near infrared radiation to analyze biological media. Apart from different medical imaging technology, the scientists took immense attention because of its low cost, with an advantage of providing anatomical information. Practically in optical tomography light is guided by fiber optics (which are placed very close to subject) to surface of object under test and detecting fibers (which are also very close to the object) are used to measure trans-illuminated or backscattered light [5] - [10]. This trans-illuminated light is scattered through multiple tissues layers and structures present inside tissues. This trans-illuminated light then converted into a series of voltages to amplify, filtered and digitized to construct an image [11] [12] [13] [14]. The data are developed either through time varying intensities or steady state complex intensities of light that have measurable amplitude and phase. The image reconstruction with optical tomography was a great challenge for scientists as it has much crosstalk noises. A lot of work is needed to implement this technique in more fields on its true basis. In this paper, an attempt has been made to give an idea to optimize the constructed image to make this technique more useful [14] - [19].
2. Optical Tomography
Optical tomography (OT) is also a non-invasive technique in which we use visible or near infrared radiation to analyze biological media. Apart from different medical imaging technology, it has captured attention of scientists and engineers due to its low cost alternative to other imaging techniques, with an advantage of providing anatomical information [1].
Practically, in optical tomography, light is guided by fiber optics (which are placed very close to object) to surface of object and detecting fibers (which are also very close to object) are used to measure trans-illuminated or backscattered light. This trans-illuminated light is scattered through multiple tissue layers and structures present inside tissues. Trans-illuminated light then converted into a series of voltages to amplify, filtered and digitized to construct an image. Data is developed either through time varying intensities or steady state complex intensities of light that has measurable amplitude and phase [1] [2] [3] [4].
Optical properties of tissue are mainly confined in scattering and absorption coefficients of light and are exploited to provide qualitative and quantitative images. Differences in optical properties between healthy and irrational tissues provide different contrasts when imaging disease. Healthy tissues properties are already measured and stored to compare with unhealthy tissues. Medical applications of optical tomography are brain imaging, skin and breast cancer imaging. Osteoarthritis, rheumatoid arthritis and diabetes treatment are also done with optical tomography [1] - [5].
3. Photon Theory
3.1. Properties of Light
Optical properties of tissues are generally describing in terms of absorption, scatter, anisotropy and refractive index. All of these parameters are dependent on wavelength of light penetrating in tissues.
3.1.1. Absorption of Light
When infrared radiation is incident on medium, resonance will occur about natural frequencies by which energy is transferred from incident field to medium and its amplitude of vibration is significantly increased. Though lifetime of excited state is around seven to ten seconds, atoms or molecules will usually lose their energy by colliding with one another within ten to twelve seconds, so raising kinetic energy of other particles involved in collisions. Hence, energy-related with incident field is most often dissipated as heat within medium. This process is known as absorption [1].
3.1.2. Refractive Index
In optics, refractive index of material defines how light propagates through that medium. Refractive index is defined as ratio between speed of light in vacuum to speed of light in medium Refractive index of medium is depending on number of particles/molecules per unit volume and polarize-ability of medium [1].
3.1.3. Scattering of Light
When an incident light propagates through medium, charged particles present in a medium are set into oscillatory motion due to electric field of incident wave and some light particles re-emit from medium with same frequency as incident wave particles, this is called scattering of light [1].
3.1.4. Anisotropy
In scattering photons do not travel in a straight line in medium and can be scattered anywhere in any angle in three dimensions. Phase function can be stated as a function of scalar product of two unit vectors which are formed by initial and final directional vectors
, which is equal to cosine of scattering angle. Anisotropy is then defined as mean cosine of scattering angle [1] [2] [3].
3.2. Photon Transport Model
Photon propagates through medium in a deterministic and foreseeable way which can be described by Boltzmann transport equation. It is also called a Radiative Transfer Equation (RTE). Basically, it is derived by considering effect of gains and losses of intensity of photon travelling in mediums such as tissues [1] - [5]. Radiative transfer equation (RTE) for photons which are travelling from point r in direction
with a specific intensity
at a time t will be given:
(1)
Explanation of each term of radiative equation is given in Table 1.
P1 and PN Approximation
Approximation is derived if first N spherical harmonics are taken, which gives
coupled PDEs. Equation (3.14) and Equation (3.15) has also been expanded into spherical harmonics. Substitution of these terms into equations for
approximation leads to following expressions:
(2)
(3)
Here
is isotropic component of source,
is anisotropic dipole-like component of source.
and
[1] [2] [3] [4].
3.3. Diffusion Approximation
By making further assumption to neglect effect of scattering, it has been assumed
Table 1. Terms of radiative transfer equation.
that source is isotropic (i.e. its properties are independent of medium geometry) and time rate of change of photon flux is too slow.
(4)
Therefore, by putting these values in Equation (3) following equation is obtained:
(5)
By substituting expression (5) in expression (2) new expression is given as:
(6)
Equation (6) is modified radiative transfer equation and is also an energy conservation equation, called diffusion equation.
(7)
Diffusion approximation assumes that source’s angular dependence, specific intensity of source and phase function of scattering parameter are uniform in order to modelled as first order spherical harmonics. To achieve these conditions, this will be defined by assuming
[1] - [10].
3.4. Boundary Conditions
To make diffusion equation specific to boundary value problem boundary condition has been applied on it. First condition is that no photons travel in an inward direction at boundary except source photons. Diffusion equation also cannot satisfy this condition precisely, rather it is assumed that total inward directed current is zero [11] - [15].
4. Image Reconstruction
Objective of image reconstruction in optical tomography is to determine optical properties (i.e.
and
) of a theoretical model which produces the same results that are obtained experimentally. According to problem statement of finding parameters we divide it into two different problems one is called forward problem and other is called inverse problem. In forward problem, optical properties are known and modelling of physical system is done for analysis. Whereas in inverse problem optical properties and geometry of medium are calculated.
4.1. Forward Problem
Forward problem can be modelled by mathematical expression given below:
(8)
where y denotes measurements which are constructed, x denotes optical properties of medium and J is sensitivity matrix which is also known as Jacobian of forward operator or weight matrix. y is calculated by knowing geometry and optical properties of medium, location of sources and detectors which are placed around medium and transport of light according to diffusion equation.
Objective of forward problem is to extract model data for comparison with experimental data during image reconstruction, or it may use to generate simulated data to test reconstruction techniques or expected outcomes of a photon study. Sensitivity matrix can be calculated numerically by using Finite Element Method (FEM) or analytically by Monte Carlo Method (MCM), or experimentally by placing a small absorber (and/or a scatterer) at discrete positions within domain.
4.2. Inverse Problem
In inverse problem, it is recommended to find unknown functions occurring in formulation of forward problem from data. This can be mathematically expressed by equation:
(9)
In this equation, y represents data types that are extracted in forward model, x represents optical properties of medium which are unknown. To solve this equation, we first invert sensitivity or Jacobian matrix which is calculated earlier in forward model. There are two characteristics of sensitivity matrix, which make difficulties in inversion process.
4.2.1. Ill-Posedness
An ill-posed problem has no solutions in desired class or has many solution, or solution procedure is unstable [15] [16] [17] [18]. In this research case, this means that there are a smaller number of independent measurements than unknown pixel values.
4.2.2. Ill-Conditionedness
A system or matrix is called ill conditioned if some small perturbation/changes in system causes a relatively large or very large change in exact solution [1]. This condition is typically resulting in amplification of both measurement errors and numerical errors in inversion process.
According to ill-posed and ill-conditioned system type, image reconstruction is also dived into two major types one is called linear image reconstruction and other is called nonlinear image reconstruction.
5. Numerical and Simulation Modelling
Simulation has been done using MATLAB algorithms. In formulation of forward and inverse modeling in most medical imaging technique, effects of scattering are neglected and path of detected photon is considered to be a straight line between a chosen source-detector pair. For low scattering environments, reconstruction problem can therefore be formulated as a system of linear equations and is directly derived from Boltzmann transport equation. For highly scattering environments, gradual impact of scattering with distance makes scattering problem notably more complex and non-linear. Solution therefore does not only involve calculation of optical properties but must also consider scattering effects inside medium. Thus, instead of apply commonly applicable direct method for solving such problems, inversion techniques are considered to evaluate which material properties would have created original signal [14] - [19].
5.1. Image Reconstruction Using Gradient Based Optimization Methods
Basic algorithms like least square (LS) methods assume a Gaussian approximation for distribution of a random variable. To obtain best fit solution for an inverse model estimation and forward model data, minimum of error function is needed which can be achieved when gradient of error function approaches to zero. Best fit measurements are then approximated by an iterative search to find a best solution to a given problem. Popular optimization techniques are the steepest descent method, conjugate gradient method and second order Levenberg Marquardt method. Method used in this paper to optimize optical tomography image reconstruction is conjugate gradient (CG) method.
5.1.1. Conjugate Gradient (CG) Methods
In this method, usually particular system of equations is tried to solve which have symmetrical and positive definite matrix. It is frequently used as an iterative search algorithm, to solve such sparse systems, which are too large to handle by a direct method. Conjugate gradient algorithms are gradient methods that perform searches along conjugate directions. Although error function decreases very quickly along negative gradient direction, but it is not necessarily that it produces fastest convergence.CG methods are divided into two categories according to nature of system of equations, i.e., Linear Conjugate Gradient (LCG) Method and Non-Linear Conjugate Gradient (NLCG) Method [10] [11] [12] [13].
5.1.2. Optimized Image Reconstruction
In light of whole discussion above, detailed algorithms for solution to optical tomography problem is summed up below:
Divide boundary domain say
area into a set of specific finite elements. Set initial nodal parameter and perturbation value. Set iterative index
and repeat it to find solution for every source and measurement position. Then measure
and
. Compute Jacobian Matrix
.
Set
to initial value for regularization.
Find a pre-estimate of
by solving Equation (4.22) in iterative SCG algorithm.
Compute residual
;
If
then increase
;
If
then k = k + 1;
Loop continue until
is approaches to desired value.
5.2. Experiments
Practical implementation of above discussed algorithms is simulated in MATLAB and result is shown in upcoming different type of scenarios.
5.2.1. Experiments/Simulations #01
In first simulation, conjugate gradient method is used to construct image. To generate forward model, MATLAB algorithm is used to define target image. Parameters that have been kept are shown in Table 2.
Target image which generated is shown in Figure 1 has only one inclusion with absorption coefficient differ than entire image located at center of medium. Image generated is based on FEM. Mesh generated for forward model has radius 25 mm. There are 16 source and measurement pairs are arranged around medium. Effect of each source on inclusion is shown in Figure 2 which is constructed using conjugate gradient squared algorithm. The complete parameters are shown in Table 2.
Objective function converges at iteration number 90 according to stopping criteria given in Table 2. Figure 3 and Figure 4 showed that actual data and reconstructed data having one inclusion in middle of image. Image quality test results are shown in Table 3.
Table 3. Image quality test results for experiment 01.
Figure 1. Experiment 01 forward model-target data generated by mathematical modeling of forward problem.
Figure 2. Effect of each source on the inclusion—source points and their effects-photon density chart w.r.t every source.
Figure 3. Forward and reconstructed data in grey image of experiment 01—which are similar to other medical imaging techniques like X-ray.
5.2.2. Experiments/Simulations #02
In second experiment, image is constructed again using NCG method but adding two inclusions having different absorption coefficients parameters and have more source and detectors shown in Table 4. Figure 5-8 showed actual data and reconstructed data with two inclusion having different absorption coefficient, sources and measurements data chart. Reconstruction algorithm contracts image at iteration number 88 according to stopping criteria given in Table 4. Image quality is also tested after image reconstruction as shown in Table 5. Test showed that with two inclusions this algorithm also reconstructs well enough information from forward model.
5.2.3. Experiments/Simulations #03
In third experiment, image is constructed again using NCG method but adding three inclusions having different absorption coefficients parameters and have more source and detectors shown in Table 6.
Figures 9-12 showed actual data and reconstructed data with two inclusions having different absorption coefficient, sources and measurements data
Table 4. Experiment #02 parameters.
Table 5. Image quality test results for experiment 02.
Figure 4. Experiment-01-image of forward problem (target data) and inverse problem (reconstructed data)—the left side image showed the distribution of internal optical properties having one inclusion—the right-side image is reconstructed using CGS method with continuous source by putting freq. ω = 0 Hz.
Figure 5. Forward model data for experiment 02.
Figure 6. Photon density charts w.r.t sources and detectors arrangement for experiment 02.
Figure 7. Experiment-02-image of target data and reconstructed data—the left side image showed the distribution of internal optical properties having different optical inclusions—the right-side image is reconstructed and showed the extracted optical inclusions.
Figure 8. Grey image data of experiment 02 using optical tomography.
chart. Reconstruction algorithm constructs image at iteration number 61 according to stopping criteria given in Table 6.
Figure 9. Forward model for experiment 03.
Figure 10. Photon density charts w.r.t sources and detectors arrangement for experiment 03.
Figure 11. Grey image data of example 03 showing black and white image cintrast for forward and inverse model data.
Image quality is again tested after image reconstruction as shown in Table 7. Test showed that with three inclusions this algorithm also reconstructs well
Figure 12. Experiment-03-image of target data and reconstructed data—the left side image showed the distribution of internal optical properties having different optical inclusions—the right-side image is reconstructed and showed the extracted optical inclusions.
Table 7. Image quality test of reconstructed data of experiment 03.
enough information from forward model. Moreover, result showed that this algorithm successfully reconstructs forward problem with very less time and gave best approximated solution of diffusion equation.
6. Conclusion
Experiments which are done in pervious examples have very encouraging results demanding a very deep study in this field. It also showed that image reconstruction of optical tomographic problems is much closed to other medical imaging techniques but has less simulation time. Moreover, image reconstruction of optical media has more anatomical data instead of other imaging techniques due to different absorptions and scattering properties of medium or layers of tissues as seen by colored images in simulations. It provides us full details of layers due to change of absorption and scattering coefficient of different layers.
Recommendation and Future Work
An efficient study is needed to fully understand limitations of reconstruction algorithms. Factors that should be considered include geometry, location of sources and detectors and proper use of a priori information.
Geometry used in given examples is very simple in which medium’s material properties are presumed to be isotropic whereas biological tissue has an anisotropic, inhomogeneous substantial composition, which can confine photon transport in certain directions. Moreover, examples/simulations are presented in this paper have been modelled as simple two dimensional geometric approximations of biological media. 3D geometries should be considered to increase realization. In order to get more accurate results, optical tomography should ultimately be used in conjunction with another imaging technique to monitor evolution of disease. Therefore, it is very advantageous to get some experimental data from MRI or CT scans so that more approximate geometry of a given tissue region may be considered.
As described earlier, in real optical imaging material properties are anisotropic, uncertainties related to scattering coefficient should be considered. Regularization and optimization techniques should therefore be replaced with a more appropriate regularization and optimization schemes such as particle swarm optimization (PSO) and genetic algorithms (GA), so that optical tomographic inversion calculations could provide more accurate and unique solutions.
In real medical imaging problem, we know that optical properties of biological tissues are changed with every new tissue layer thus, the refractive index mismatches between boundary regions/layers are largely unknown and therefore require further experimental investigation to estimate and regularize a prior information. Further, in Pakistan there is no such advanced laboratory to research on optical imaging till now, so an implementation of a prototype for optical tomography could be done.
Acknowledgements
The authors are grateful to the Electrical Engineering Department of Air University, Islamabad, Pakistan and Electrical Engineering Department of University of Lahore, Lahore, Pakistan for providing us the facility for the conduction of research in medical imaging.