Application of Reconstruction and Optimization Algorithms in Optical Tomography

Optical tomography is a non-invasive technique that uses visible or near infrared radiation to analyze biological tissues. Researchers take immense attention towards advancement in optical tomography because of its low cost and an advantage of providing anatomical information. Based on the information of optical characteristics, forward and inverse problem of tomography are solved. In this research, finite element method is employed for forward problem and gradient-based optimization algorithm is developed for inverse problem of optical tomography. It is found from simulations that information about imaging is processed more distinctly and in less computational time. Normal and abnormal conditions in imaging are readily distinguished. Simulations are carried out in Matlab. Different scenarios are developed and are simulated to validate the performance of reconstruction and optimization algorithms in optical tomography.


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

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

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.

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

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

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

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 ( ) , a b , which is equal to cosine of scattering angle. Anisotropy is then defined as mean cosine of scattering angle

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 ( ) , , I r s t at a time t will be given: Explanation of each term of radiative equation is given in Table 1.

P1 and PN Approximation
N P Approximation is derived if first N spherical harmonics are taken, which gives ( ) 2 1 N + coupled PDEs. Equation (3.14) and Equation (3.15) has also been expanded into spherical harmonics. Substitution of these terms into equations for 1 P approximation leads to following expressions:

Diffusion Approximation
By making further assumption to neglect effect of scattering, it has been assumed    that source is isotropic (i.e. its properties are independent of medium geometry) and time rate of change of photon flux is too slow.
Therefore, by putting these values in Equation (3) following equation is obtained: By substituting expression (5) in expression (2) new expression is given as: Equation (6) is modified radiative transfer equation and is also an energy conservation equation, called diffusion equation.
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 s

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

Image Reconstruction
Objective of image reconstruction in optical tomography is to determine optical properties (i.e. a µ and s µ ) 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.

Forward Problem
Forward problem can be modelled by mathematical expression given below: 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.

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: 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.

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.

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.

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

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.

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 itera-

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 0 iter = and repeat it to find solution for every source and measurement position. Then measure k x and x δ . Compute Jacobian Matrix ( )Ĵ x .
Set α to initial value for regularization. Find a pre-estimate of is approaches to desired value.

Experiments
Practical implementation of above discussed algorithms is simulated in MATLAB and result is shown in upcoming different type of scenarios.

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.

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.

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        chart. Reconstruction algorithm constructs image at iteration number 61 according to stopping criteria given in Table 6. World Journal of Engineering and Technology   Image quality is again tested after image reconstruction as shown in Table 7. Test showed that with three inclusions this algorithm also reconstructs well World Journal of Engineering and Technology 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.

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 aniso-tropic, 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.