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.

Share and Cite:

Hashmi, U. , Muzzammel, R. , Arshad, R. and Mehmood, S. (2020) Application of Reconstruction and Optimization Algorithms in Optical Tomography. World Journal of Engineering and Technology, 8, 188-203. doi: 10.4236/wjet.2020.82016.

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 ( a , b ) , 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 s ^ with a specific intensity I ( r , s ^ , t ) at a time t will be given:

( 1 c t + s ^ + μ t ( r ) ) I ( r , s ^ , t ) = μ s ( r ) 4 π 0 p ( s ^ s ^ ) I ( r , s ^ , t ) d s ^ + Q ( r , s ^ , t ) (1)

Explanation of each term of radiative equation is given in Table 1.

P1 and PN Approximation

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

( 1 c t + μ a ( r ) ) Φ ( r , t ) + J ( r , t ) = Q o ( r , t ) (2)

( 1 c t + 1 3 κ ( r ) ) J ( r , t ) + 1 3 Φ ( r , t ) = Q 1 ( r , t ) (3)

Here Q o is isotropic component of source, Q 1 is anisotropic dipole-like component of source. κ = Diffusion coefficient = 1 3 ( μ a + μ s ) and μ s = ( 1 p 1 ) μ s [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.

q 1 = 0 & J t = 0 (4)

Therefore, by putting these values in Equation (3) following equation is obtained:

J ( r , t ) = κ ( r ) Φ ( r , t ) (5)

By substituting expression (5) in expression (2) new expression is given as:

( 1 c t + μ a ( r ) ) Φ ( r , t ) κ ( r ) Φ ( r , t ) = Q o ( r , t ) (6)

Equation (6) is modified radiative transfer equation and is also an energy conservation equation, called diffusion equation.

1 c × Φ ( r , t ) t Flow Energy + μ a ( r ) × Φ ( r , t ) Absorption Energy + ( κ ( r ) × Φ ( r , t ) ) Diffusion Energy = q o ( r , t ) Source Energy (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 μ s μ a [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. μ 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.

4.1. Forward Problem

Forward problem can be modelled by mathematical expression given below:

y = J x (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:

x = J 1 y (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 i t e r = 0 and repeat it to find solution for every source and measurement position. Then measure x k and δ x ^ . Compute Jacobian Matrix J ^ ( x ^ ) .

Set α to initial value for regularization.

Find a pre-estimate of x k + 1 by solving Equation (4.22) in iterative SCG algorithm.

Compute residual k + 1 ;

If k + 1 k then increase α ;

If x k + 1 = x ( r ) k + 1 then k = k + 1;

Loop continue until k + 1 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 2. Example #01 parameters.

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.

Table 6. Experiment 03 parameters.

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.


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.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.


[1] Muzzammel, R. and Ehsan, M. (2015) Forward and Inverse Problem Formulation of Optical Tomography Based on Equation of Radiative Transfer. CJAST, 8, 194-203.
[2] Hielscher, A.H., Bluestone, A.Y., Abdoulaev, G.S., Klose, A.D., Lasker, J., Stewart, M., Netz, U. and Beuthan, J. (2002) Near-Infrared Diffuse Optical Tomography. Disease Markers, 18, 313-337.
[3] Pifferi, A., Contini, D., Mora, A.D., Farina, A., Spinelli, L. and Torricelli, A. (2016) New Frontiers in Time-Domain Diffuse Optics: A Review. The Journal of Biomedical Optics, 21, Article ID: 091310.
[4] Puszka, L.H., Planat-Chrétien, A., Koenig, A., Derouard, J. and Dinten, J.-M. (2013) Time-Domain Reflectance Diffuse Optical Tomography with Mellin-Laplace Transform for Experimental Detection and Depth Localization of a Single Absorbing Inclusion. Biomedical Optics Express, 4, 569-583.
[5] Durduran, T., Choe, R., Baker, W.B. and Yodh, A.G. (2010) Diffuse Optics for Tissue Monitoring and Tomography. Reports on Progress in Physics, 73, Article ID: 076701.
[6] Schweiger, M., Arridge, S.R., Hiraoka, M. and Delpy, D.T. (1995) The Finite Element Method for the Propagation of Light in Scattering Media: Boundary and Source Conditions. Medical Physics, 22, 1779-1792.
[7] Pogue, B.W., Patterson, M.S., Jiang, H. and Paulsen, K.D. (1995) Initial Assessment of a Simple System for Frequency Domain Diffuse Optical Tomography. Physics in Medicine & Biology, 40, 1709-1729.
[8] Arridge, S.R., Schweiger, M., Hiraoka, M. and Delpy, D.T. (1993) A Finite Element Approach for Modeling Photon Transport in Tissue. Medical Physics, 20, 299-309.
[9] Arridge, S.R. (1999) Optical Tomography in Medical Imaging. Inverse Problems, 15, R41-R93.
[10] Klibanov, M.V. and Lucas, T.R. (1999) Numerical Solution of a Parabolic Inverse Problem in Optical Tomography Using Experimental Data. SIAM Journal on Applied Mathematics, 59, 1763-1789.
[11] Klibanov, M.V., Lucas, T.R. and Frank, R.M. (1997) A Fast and Accurate Imaging Algorithm in Optical/Diffusion Tomography. Inverse Problems, 13, 1341-1361.
[12] Tang, X.L., Han, C.P. and Liu, X. (2019) Super-Resolution X-Ray Luminescence Optical Tomography Imaging. Proceedings of SPIE 11190, Optics in Health Care and Biomedical Optics IX, 20 November 2019, 111902M.
[13] Byrd, B.K., Filan, C.E., Folaron, M.R., Strawbridge, R.R., Wirth, D.J. and Davis, S.C. (2019) Multi-Angle Projection Imaging of Short-Wave Infrared (SWIR) Fluorescence for Small Animal Optical Tomography. Proceedings of SPIE 10862, Molecular-Guided Surgery: Molecules, Devices, and Applications V, 7 March 2019, 108620E.
[14] Lyons, A., Boccolini, A., Tonolini, F., Repetti, A., Henderson, R., Wiaux, Y. and Faccio, D.F.A. (2019) Computational Time-of-Flight Diffuse Optical Tomography. Nature Photonics, 13, 575-579.
[15] Yoo, J., et al. (2020) Deep Learning Diffuse Optical Tomography. IEEE Transactions on Medical Imaging, 39, 877-887.
[16] Abdoulaev, G.S., Ren, K. and Hielscher, A.H. (2005) Optical Tomography as a PDE-Constrained Optimization Problem. Inverse Problems, 21, 1507-1530.
[17] Chalia, M., Dempsey, L.A., Cooper, R.J., et al. (2019) Diffuse Optical Tomography for the Detection of Perinatal Stroke At The Cot Side: A Pilot Study. Pediatric Research, 85, 1001-1007.
[18] Lo, P.A. and Chiang, H.K. (2019) Three-Dimensional Fluorescence Diffuse Optical Tomography Using the Adaptive Spatial Prior Approach. Journal of Medical and Biological Engineering, 39, 827-834.
[19] Arridge, S.R. and Hebden, J.C. (1997) Optical Imaging in Medicine II: Modelling and Reconstruction. Physics in Medicine and Biology, 42, 841-853.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.