Soft Tissue Image Reconstruction Using Diffuse Optical Tomography

Abstract

Diffuse optical tomography (DOT) is favorable to analyze physical records in organic tissue with a specific purpose by means of a method related to the forward problem and the inverse solution. This study develops morphological soft tissue realization using an image reconstruction algorithm constructed on multifrequency DOT in Near-Infra-Red (NIR) wavelength. Forward problem solves the Diffusion Equation to compute the optical flux distributed in the phantom geometrical model. Inverse solution, the image is reconstructed using the absorption and reduced scattered coefficients under different boundary conditions. The inverse image reconstruction algorithm is tested for several simulation, with variation in background contrast ratios for different frequencies are simulated. The image reconstruction in DOT eliminates spatial resolution by optimizing source-detector separation and modulation intensities of the source.

1. Introduction

Forward problems are used to explain the propagation of photons within a tissue and to calculate the optical flux at the tissue boundary. The reconstruction of the tissue image is the inverse problem, from light measurements at the boundary phantom surface, tissue absorption, scattering factors, and optical flow [1]. The inverse problem is difficult to solve due to the issue of fairness. This indicates that the problem is not properly configured. Appropriate problem characteristics include the existence of a solution, the uniqueness of the solution, and the constant reliance on the data [2]. The third property determines the stability of the solution and is important for determining the inverse problem. The ill-posedness problem occurs when the problem solution does not depend on the data indefinitely. Small changes to the data can make a big difference in the solution in this case. Regularization method is used to solve this problem, which is a regularization method that introduces additional information in order to create well-posed data [3, 4]. Diffuse optical imaging [5, 6, 7, 8] is a technique that uses an MRI scan and X-rays generate spatially decomposed images and uses high-resolution complementary structural information to improve low-resolution functional images. A set of fiber optics has been connected the object’s boundary in experimental systems. The light source was a near-infrared (NIR) laser source that was diffused on the phantom, the scattered rays were measured with a photodetector [9]. Regularization method [10] is used to remove the ill-posedness, with the Levenberg–Marquardt method (LM) being one of the most commonly used methods. Following the regularization process, the Split Bregman reconstruction method [10, 11, 12] is used to reconstruct a soft tissue image. The sparsity regularization technique for image reconstruction in DOT is described by Bo Bi et al. [9]. Gehre et al. [13] investigates the possibility of sparsity constraints in the inverse problem of deriving distributed conductivity from critical potential measurements in electrical impedance tomography (EIT). Chamorro et al. [14] proposed an Algebraic reconstruction technique—Split Bregman (ART-SB) algorithm solved the L1-regularized problem. Wang et al. evaluated the Split Bregman iteration algorithm for the L1 norm regularization inverse problem in electrical impedance tomography. Figueiredo et al. [15] investigated the use of Split Bregman iterative algorithms for the L1-norm regularized inverse problem of electrical impedance constrained quadratic programming ill-tomography formulation.

2. Inverse model solution

In most cases, the measurement data in DOT reconstruction is derived from the numerical solution of the forward problem. Regularization techniques are used to eliminate the obtrusive inverse problem variables. The measurement technology of optical devices is so limited that their existence cannot be accurately determined from all angles. Instead, it gets the average contact angle data for the phantom. The purpose is to reconstruct the image from known scattering and absorption coefficients, which are assumed to be known. The reconstructed result is obtained by comparing the true value to the measured value. A variety of practical reconstruction algorithms have been developed in tomography to implement the process of reconstructing a 3D object from a projection [16]. These algorithms are mainly based on the mathematics of statistical knowledge of the data acquisition process and the geometry of the data imaging system [17].

2.1 Levenburg-Marquardt method

The inverse problem is used to reconstruct an image in the following ways by estimating the scattering coefficient, absorption coefficient, and optical flux [18, 19]. The noise level is present in the actual measurement; both actual measurement data and actual data are shown here [9]. The following nonlinear equation is used to solve the inverse problem of DOT for i = 1,…..,s. Assume you know the total attenuation coefficient.

$$F_i\left({\mu }_t,{\mu }_s\right)=M_i^{\delta },\left({\mu }_t,{\mu }_s\right)\in D$$E1

The inverse problem of DOT is inappropriate and uses regularization techniques to reconstruct the image [8], i.e. the Tikhonov functional feature that is minimized for the coefficient. R (μ_S) is a penalty function for regularization [4]. By analyzing the minimization problem,

$$J\left({\mu }_s\right)=\frac{1}{2}\sum _{i=1}^s{‖F_i\left({\mu }_s\right)-M_i‖}_{L^2\left(\partial X\right)}^2+\mathit{\alpha R}\left({\mu }_s\right)$$E2

Over the set,

$$Q_{\mathit{ad}}=\left\{{\mu }_s\in L^{\infty }\left(X\right)\right\}$$E3

$${}_{{\mu }_s\in Q_{\mathit{ad}}}{}^{\inf }_J\left({\mu }_s\right)$$E4

The standard reconstruction method is considered using Eqs. (3) and (4).

2.1.1 Standard reconstruction

Traditional norm-squared penalties are believed to reduce the following functions,

$$R\left({\mu }_s\right)=\frac{1}{2}{‖\nabla {\mu }_s‖}_{L^2\left(X\right)}^2$$E5

$$J\left({\mu }_s\right)=\frac{1}{2}\sum _{i=1}^s{‖F_i\left({\mu }_s\right)-M_i^{\delta }‖}_{L^2\left(\partial x\right)}^2+\frac{\alpha }{2}{‖{\mu }_s-{\mu }_s^{\ast }‖}_{L^2\left(X\right)}^2$$E6

In the inverse problem of DOT, the Levenberg–Marquardt regularization method [20] is used. The forward operator is linearized around the initial estimation for each;

$$F_i\left({\mu }_s\right)=F_i\left({\mu }_s^0\right)\left({\mu }_s-{\mu }_s^0\right)+R\left({\mu }_s^0;i\right)$$E7

where Eq. (7) denotes the Taylor remainder for the linearization around and the Frechet derivative is obtained by substituting the above equation and ignoring the higher-order remainder [13].

$$\underset{{\mu }_s\in D}{\mathit{inf}}\frac{1}{2}\sum _{i=1}^s{‖F_i\left({\mu }_s^0\right)+F_i^{\prime }\left({\mu }_s^0\right)\left({\mu }_s-{\mu }_s^0\right)-M_i^{\delta }‖}_{L^2\left(\partial x\right)}^2$$E8

The Euler equation for discrete problems is

$$\sum _{i=1}^s{F}_i^{\prime }{\left({\mu }_s^0\right)}^{\ast }\left(F_i\left({\mu }_s^0\right)+F_i^{\prime }\left({\mu }_s-{\mu }_s^0\right)-M_i^{\delta }\right)+\alpha \left({\mu }_s-{\mu }_s^0\right)=0$$E9

Solving this (9) yields the final Equation [16].

That is,

$$\left(\sum _{i=1}^s{F}_i^{\prime }{\left({\mu }_s^0\right)}^{\ast }{F}_i^{\prime }\left({\mu }_s^0\right)+\mathit{\alpha I}\right)\left({\mu }_s-{\mu }_s^0\right)=-\sum _{i=1}^s{F}_i^{\prime }{\left({\mu }_s^0\right)}^{\ast }\left(F_i\left({\mu }_s^0\right)-M_i^{\delta }\right)$$E10

where I is the identity matrix to solve the new estimate of ${\mu }_s$ based on the initial guess ${\mu }_s^0$.

2.2 Sparsity reconstruction

Sparsity reconstruction function can be minimized as

$$J\left({\mu }_s\right)=\frac{1}{2}\sum _{i=1}^s{‖F_i\left({\mu }_s\right)-M_i^{\delta }‖}_{L^2\left(\partial x\right)}^2+\frac{\alpha }{2}{‖{\mu }_s‖}_{l_1}$$E11

Such that d= ${\mathit{\mu }}_{\mathit{s}}$. Decouple the L1 and L2 portion in (11). The constrained problem [10].

where $\beta \succ 0$ is a split parameter, and iteratively, the next subproblem can be solved as [3, 4].

$$\begin{array}{l}\left({\mu }_s^k,d^k\right)=\mathit{\text{argmin}}\frac{1}{2}\sum _{i=1}^s{‖F_i\left({\mu }_s\right)-M_i^{\delta }‖}_{L^2\left(\partial x\right)}^2+\alpha {‖d‖}_{l_1}+\frac{\beta }{2}{‖d-{\mu }_s-b_d^{k-1}‖}_2^2,\\ b_d^k=b_d^{k-1}+{\mu }_s^k-d^k\end{array}$$E12

By dividing the minimization of (12) and d separately, the sub-problem can be minimized.

1. Consider

   $${\mu }_s^k=\mathit{\text{argmin}}\frac{1}{2}\sum _{i=1}^s{‖F_i\left({\mu }_s\right)-M_i^{\delta }‖}_{L^2\left(\partial x\right)}^2+\frac{\beta }{2}{‖d^{k-1}-{\mu }_s-b_d^{k-1}‖}_2^2$$E13

2. Consider

   $$d^k=\underset{d}{\mathit{\text{argmin}}}{‖d‖}_{l_1}+\frac{\beta }{2}{‖d-{\mu }_s^k-b_d^{k-1}‖}_2^2$$E14

3. Consider

   $$b_d^k=b_d^{k-1}+{\mu }_s^k-d^k$$E15

Minimize ${\mathit{\mu }}_{\mathit{s}}^{\mathit{k}}\mathit{,}{\mathit{d}}^{\mathit{k}}$ as

$${\mu }_s^k=\underset{{\mu }_s}{\mathit{\text{argmin}}}\frac{1}{2}\sum _{i=1}^s{‖F_i\left({\mu }_s^{k-1}\right)+F_i^{\prime }\left({\mu }_s^{k-1}\right)\left({\mu }_s-{\mu }_s^{k-1}\right)-M_i^{\delta }‖}_2^2+\frac{\beta }{2}{‖d^{k-1}-{\mu }_s-b_d^{k-1}‖}_2^2$$E16

The variational equation is given as

$$\begin{array}{l}\left(\sum _{i=1}^s{F}_i^{\prime }{\left({\mu }_s^{k-1}\right)}^{\ast }{F}_i^{\prime }\left({\mu }_s^{k-1}\right)+\mathit{\beta I}\right)\left({\mu }_s-{\mu }_s^{k-1}\right)\\ =\beta \left(d^{k-1}-{\mu }_s^{k-1}-b_d^{k-1}\right)+\sum _{i=1}^s{F}_i^{\prime }{\left({\mu }_s^{k-1}\right)}^{\ast }\left(F_i\left({\mu }_s^{k-1}\right)-M_i^{\delta }\right)\end{array}$$E17

L1 is solved efficiently by the contraction operator; that is

$${\mathit{d}}^{\mathit{k}}\mathbf{=}\mathit{\text{shrink}}\left({\mathit{\mu }}_{\mathit{s}}^{\mathit{k}}\mathbf{+}{\mathit{b}}_{\mathit{d}}^{\mathit{k}\mathbf{-}\mathbf{1}},\frac{\mathit{\alpha }}{\mathit{\beta }}\right)$$E18

Where the shrinkage operator is defined as

$$\mathit{\text{Shrink}}\left(x,t\right)=sign\left(\mathit{x}\right)max\left(|\mathit{x}|-\mathrm{t},0\right)=\left\{\begin{array}{c}\mathit{x}-\mathit{t},\mathit{x}\ge \mathit{t},\\ \mathrm{0},\mid \mathit{x}\mid <\mathit{t},\\ \mathit{x}+\mathit{t},\mathit{x}\le -\mathit{t}\end{array}\right.$$E19

The rate of spilt Bregman is highly dependent on the rate of dissolution ${\mathit{F}}_{\mathit{i}}^{\prime }\left({\mathit{\mu }}_{\mathit{s}}\right)$.

3. Levenberg: Marquardt algorithm

Because the DOT image has poor spatial resolution due to severe ill-posedness, the regularization technique is used in conjunction with reconstruction algorithms to reconstruct the images.

Algorithm

Input: Set the initial estimation ${\mu }_s^0$; The regularization parameter α,β.

Output: Approximate minimizer ${\mu }_s$.

for k = 1,……….k do.

for i = 1,……….s do.

 i. Compute the Frechet derivative $F_i^{\prime }$(${\mu }_s^k$), and $F_i^{\prime }$(${\mu }_s^k$)^∗.

end for.

 i. Compute ${\sum }_{i=1}^s{F}_i^{\prime }$(${\mu }_s^k$)^∗$F_i^{\prime }$(${\mu }_s^k$) + αI, and −${\sum }_{i=1}^s{F}_i^{\prime }$(${\mu }_s^k$)^∗($F_i^{\prime }$(${\mu }_s^k)$-$M_i^{\delta }$).

 ii. Update ${\mu }_s^{k+1}$ by solving the linearization problem

$$\left({\sum }_{i=1}^s{F}_i^{\prime }\left({\mu }_s^k\right){}^{\ast }{F}_i^{\prime }\left({\mu }_s^k\right)+\mathrm{\alpha }\mathrm{I}\right)\left({\mu }_s-{\mu }_s^k\right)=-{\sum }_{i=1}^s{F}_i^{\prime }\left({\mu }_s^k\right){}^{\ast }\left(F_i^{\prime }\left({\mu }_s^k\right)-M_i^{\delta }\right)$$

 iii. Check the stopping criterion.

end for.

4. Spilt Bregman algorithm

Input: set the initial guess${\mathit{\mu }}_{\mathbf{0}}^{\mathit{s}}$;

Regularization parameters $\mathit{\alpha }\mathit{>}\mathbf{0}\mathit{,}\mathit{\beta }\mathit{>}\mathbf{0}\mathit{;}{\mathit{d}}^{\mathbf{0}}\mathit{=}{\mathit{b}}_{\mathit{d}}^{\mathbf{0}}\mathit{=}\mathbf{0}$ and margins of error $\mathit{\epsilon }$.

Output: Outputs an approximate value${\mathit{\mu }}_{\mathit{s}}\mathit{=}{\mathit{\mu }}_{\mathit{s}}^{\mathit{k}}$.

While $‖{\mathit{\mu }}_{\mathit{s}}^{\mathit{k}}\mathit{-}{\mathit{\mu }}_{\mathit{s}}^{\mathit{k}\mathit{-}\mathbf{1}}‖\mathit{>}\mathit{\epsilon }$ do

1. For each $\mathbf{1}\mathit{\le }\mathit{i}\mathit{\le }\mathit{s}$, calculate (17) to be acquired ${\mathit{\mu }}_{\mathit{s}}^{\mathit{k}}$;

2. ${\mathit{d}}^{\mathit{k}}\mathbf{=}\mathit{\text{shrink}}\left({\mathit{\mu }}_{\mathit{s}}^{\mathit{k}}\mathbf{+}{\mathit{b}}_{\mathit{d}}^{\mathit{k}\mathbf{-}\mathrm{1}},\mathit{\alpha }/\mathit{\beta }\right)\mathbf{;}$

3. $\text{Calculate}{b}_d^k=b_d^{k-1}+{\mu }_s^k-d^k$

4. ${\mathit{\mu }}_{\mathit{s}}={\mathit{\mu }}_{\mathit{s}}^{\mathit{k}}$

end while.

The split Bregman method entails locating the Fréchet derivative, which is nothing more than the first order derivative function. The algorithm is built with regularization parameters in mind.

5. Simulation result of image reconstruction

The reconstructed image can be obtained using the Levenberg–Marquardt algorithm by providing optical flux, scattering coefficient, and absorption coefficient values, which are then compared to distinguish between normal soft tissue and cancer-affected tissue. The forward mesh has more nodes to extract all of the optical parameters of the tissue, whereas the inverse mesh has fewer nodes for reconstruction. As shown in Figure 1, the forward mesh has 1097 nodes and 2095 elements, and the reverse mesh has 286 nodes and 522 elements.

Figure 1.

Mesh diagram of inverse problem.

Figure 2 depicts the reconstruction based on absorption and scattering coefficient measurements. The image is reconstructed using optical properties of human tissue such as absorption coefficient and scattering coefficient (and (r) = 2). The reconstructed image is based on the absorption and scattering coefficient values. It is possible to predict normal tissue and cancer-affected tissue by examining the reconstructed image with absorption and scattering coefficients. When the absorption and scattering coefficients are higher, the tumor is classified as malignant or benign soft tissue tumor.

Figure 2.

Levenberg–Marquardt regularization and standard reconstruction.

The split Bregman algorithm with sparsity regularization efficiently solves the DOT image reconstruction problem. Figure 3 depicts a Spilt Bregman reconstruction from scattering coefficients. The scattering coefficient value distinguishes the variation of abnormal tissue to normal tissue. According to Figure 3, the abnormal tissue scattering coefficient ranges from 150 to 210, whereas the normal tissue scattering coefficient ranges from 0 to 20. When compared to other reconstruction algorithms, the Bregman algorithm produces more accurate results. The reconstructed image’s resolution is determined by calculating the signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), relative solution error norm (RE), and CPU time. SNR is calculated as follows:

Figure 3.

Spilt Bregman regularization with sparsity reconstruction.

$$\mathit{SNR}=10{\mathit{log}}_{10}\left(\frac{P_{\mathit{\text{signal}}}}{P_{\mathit{\text{noise}}}}\right)$$E20

The CNR is a metric used to assess image quality. The mean and standard deviation values are used to calculate it. CNR is calculated as follows:

$$\mathit{CNR}=\frac{S_A-S_B}{{\sigma }_0}$$E21

where are the image signal intensities and is the standard deviation of pure image noise. The Relative solution error norm is computed as follows:

$$\mathit{E}=\frac{{‖{\mathit{\mu }}_{\mathit{S}}-{\mathit{\mu }}_{\mathit{S}}^{\mathit{\text{true}}}‖}_{\mathbf{2}}}{{‖{\mathit{\mu }}_{\mathit{s}}^{\mathit{\text{true}}}‖}_{\mathbf{2}}}$$E22

Table 1 compares parameters used to evaluate the performance of reconstruction algorithms. The Split Bregman method has a higher SNR than the Gauss Newton algorithm and improves CNR more than the Gauss Newton method. To achieve better performance, the RE of a reconstructed image should be low. Because the Gauss Newton method has a high RE value, it is not an optimal solution for image reconstruction. Finally, when compared to the Gauss Newton method, the Split Bregman method requires less CPU time to execute. The graph of the performance analysis of the Split Bregman and Gauss Newton algorithms is shown in Figure 4.

Parameters	Split Bregman method	Gauss Newton
SNR	9.2327	4.3402
CNR	66.947	39.743
RE	0.0508	0.2141
CPU time (s)	72.231	75.197

Table 1.

Performance analysis of reconstruction algorithms.

Figure 4.

Performance analysis of reconstruction algorithms.

6. Conclusion

The solution to diffuse light transport through tissues is provided by iterative non-linear reconstruction of diffuse optical tomography using a finite element forward model. The efficiency of the forward solver has a significant impact on reconstruction performance and reconstruction time, which is critical in making optical tomography a viable imaging modality in clinical diagnosis. Standard regularization (Levenburg-Marquadt) with a small anisotropy factor identifies the scattering coefficient better than sparsity regularization in the inverse model. Sparsity regularization (Split Bregman) localizes the inclusion position and has high anisotropy factor g while forward-peaking region. The absorption and scattering coefficient values of the reconstructed it is analyzed to determine the difference between normal soft tissue and cancerous tissue. Increasing the number of measurements by adding more photo detectors is one way to improve the quality of a reconstructed image. Finally, a regularization technique is used to remove the ill-posedness problem, and a Split Bregman reconstruction algorithm is used to achieve a high-resolution image.