Diffuse optical tomography (DOT) is an emerging medical imaging modality that uses nonionizing near-infrared (NIR) light in the range of 600 to 1000 nm as the probing medium with a capability of estimating the pathophysiological state of the tissue under investigation.¹ Multiple scattering of the NIR light in the turbid medium inherently limits the resolution of the reconstructed diffuse optical image to around 0.5 cm,¹ which is at least five times lower than in traditional medical imaging modalities, such as computed tomography (CT) and magnetic resonance imaging (MRI). Recent advances in developing hybrid optical imaging modalities, which combine the traditional imaging modalities (MRI/CT) with diffuse optical imaging, have proven that the diagnostic/therapeutic information provided by hybrid optical imaging is of greater value than that of the stand-alone individual imaging modalities.¹

One of the main aspects of hybrid optical imaging is utilization of structural prior information provided by the traditional imaging modality to either guide or improve the solution. Among methods that guide the solution using the structural information, soft priors and hard priors are prominent.^2–5 Other approaches, such as use of region boundary information,⁵ makes the estimation problem into a two-stage approach, first performing the hard priors, and second incorporating the obtained solution as a initial guess into the traditional reconstruction scheme (no-priors). The soft-priors approach uses the structural information in the regularization scheme to effectively guide the solution.^2–4 In this, the number of optical parameters to be estimated remain the same as in a traditional approach, except that the regularization matrix consists of the structural information. In the hard-priors approach, the number of optical parameters that needs to be reconstructed is reduced to the number of distinct regions identifiable with traditional imaging, leading to substantial reduction in the problem size. For the case of breast imaging, the typical regions that are identifiable are adipose (fatty), fibro-glandular, and tumor,^4,6 resulting in the number of optical parameters to be reconstructed equal to three.

The hard-priors approach makes the diffuse optical image reconstruction problem less ill-posed compared to the traditional approach (including soft priors). Specifically, this converts the problem from underdetermined in nature to overdetermined. Even then, the reconstruction procedures that are adapted in the hard-priors case are the same as those used in the traditional approach, the most popular being Levenberg–Marquardt (LM) minimization scheme.^2–5 LM scheme requires a regularization parameter, which not only controls the convergence (number of iterations required),⁷ but also sometimes leads to biased solutions. Moreover, the LM scheme that is typically adapted to diffuse optical imaging requires calculation of the Jacobian (first-order derivative of the model).⁸

This work aims to show that a gradient-free simplex method,^9,10 which does not require regularization or computation of gradient (or its variant Jacobian), is highly efficient compared to the existing methods for solving the image reconstruction problem with the hard-priors approach. The reconstruction results using the proposed method are compared with traditional methods using numerical and gelatin-phantom experimental cases.

Continuous-wave (CW) light propagation through a highly scattering medium such as biological tissues is modeled using an approximation of radiative transport equation called diffusion equation (DE),¹¹ given by

where $\mathrm{\Phi }(r)$ represents the photon density (real values) at position $r$. The isotropic light source is given by $Q_o(r)$ at $r$. The diffusion coefficient $D(r)$ is defined as $D(r)=1/\{3[{\mu }_a(r)+{\mu }_s^{\prime }(r)]\},$ with ${\mu }_s^{\prime }(r)$ and ${\mu }_a(r)$ representing the reduced scattering coefficient (assumed to be known for CW case) and absorption coefficient, respectively. For the above DE, the type III (Robin-type) boundary condition is used, as it models the refractive index mismatch at the boundary of the tissue. This partial differential equation, along with the boundary condition, is solved using finite element method (FEM),¹² and $\mathrm{\Phi }(r)$ is obtained. The modeled data [$G({\mu }_a)$] for the CW case becomes the natural logarithm of the amplitude data.¹²

The estimation of the optical absorption coefficients (here, equal to the number of regions) involves matching the experimentally measured boundary data ($y$) iteratively with $G({\mu }_a)$ in a least-square sense by varying ${\mu }_a$ in different regions.² Explicitly, the objective function that needs to be minimized with respect to ${\mu }_a$ here is given by

where the dimension of $y$ and ${\mu }_a$ is given by $M$ (number of measurements) and $R$ (number of regions), respectively. Even though $R\ll M$, the minimization problem is still ill-posed (also nonlinear), requiring regularization. The popular LM minimization scheme results in an updated equation:² where $J^T$ represents the transposed Jacobian $J$($=\partial G({\mu }_a)/\partial {\mu }_a$, dimension: $M\times R$), $I$ is the identity matrix, $\mathrm{\Delta }{\mu }_a$ is the update, and the regularization parameter is represented by $\lambda$. The $\lambda$ is systematically decreased with each iteration, here by a factor ${10}^{0.25}$. As $J$ is known to be semipositive definite in the case of hard priors,² the $\lambda$ affects the convergence properties of the solution.⁷ After every iteration, the ${\mu }_a$ is updated by adding $\mathrm{\Delta }{\mu }_a$ followed by recomputing of $J$ and $G({\mu }_a)$ before proceeding to Eq. (3). This procedure is repeated until the change in $\mathrm{\Omega }$ [Eq. (2)] becomes $<2\%$ between successive iterations, ensuring the convergence of the solution.

As the convergence is dependent on $\lambda$ for the LM minimization, a direct search method that can minimize a scalar-valued nonlinear function, such as $\mathrm{\Omega }$ [Eq. (2)], is the Nelder and Mead simplex algorithm.^9,10 This algorithm is extremely popular for the unconstrained minimization problem, which does not require explicit or implicit calculation of the derivatives (or its variant gradients),⁹ making it one of the most efficient techniques. The only requirement for this algorithm is that it involves repeated calculation of $G({\mu }_a)$. Because the objective function ($\mathrm{\Omega }$) has only $R$ number of parameters, the direct search method forms the initial simplex with $R+1$ points, with each point having a dimension of $R\times 1$ using the initial guess of ${\mu }_a$. These form the vertices of the simplex. The simplex method of minimization is achieved through a series of steps deployed repeatedly. Let the maximum value of the objective function ($\mathrm{\Omega }$) be at the point ${\mu }_a^R$, which is replaced by a new point ${\mu }_a^r$ via reflection, i.e.,

where ${\overline{\mu }}_a$ represents the centroid of the the simplex and $\alpha$ is the reflection coefficient (positive value here it is 1). If the $\mathrm{\Omega }({\mu }_a^r)$ lies between $\mathrm{\Omega }({\mu }_a^1)$ (lowest value) and $\mathrm{\Omega }({\mu }_a^R)$, then ${\mu }_a^r$ becomes the ${\mu }_a$ and we terminate the iteration. Otherwise, if $\mathrm{\Omega }({\mu }_a^r)<\mathrm{\Omega }({\mu }_a^1)$, then we expand using where $\chi$ is the expansion coefficient taking a value of 2. If $\mathrm{\Omega }({\mu }_a^e)<\mathrm{\Omega }({\mu }_a^r)$, then we make ${\mu }_a={\mu }_a^e$ if $\mathrm{\Omega }({\mu }_a^e)\ge \mathrm{\Omega }({\mu }_a^r)$, then ${\mu }_a={\mu }_a^r$ and we terminate the iteration. In case the expand condition is not satisfied, then the contract step is taken, similar to earlier $\mathrm{\Omega }$, and using the contraction, ${\mu }_a$ is compared and updated.¹⁰ This follows a shrinkage step.¹⁰ The coefficients that are used here are 0.5 (contract) and 0.5 (shrinkage). Note that these coefficients are universally accepted¹⁰ as the standard and fixed for the method proposed here. More details of the algorithm, along with the flow chart and tie-breaking rules, are given in Ref. 10. The procedure is repeated until the change in ${\mu }_a$ between the successive iterations becomes $<{10}^{-12}$ or the objective function reaches a value $<0.1\%$ of the initial objective function value $[\mathrm{\Omega }({\mu }_{a0})]$. The simplex algorithm requires repeated computation of $\mathrm{\Omega }$ (Eq. 2), which in turn evaluates $G({\mu }_a)$ at every step, for converging to solution. This algorithm is known to provide an optimal solution when the parameter space dimensionality is within the range of data space.¹⁰ So the hard-priors method becomes a good case for deployment of this method.

To effectively assess the reconstruction performance of the proposed method, a numerical experiment involving a typical MRI-NIR human breast case is considered. The target distribution of the same is shown in the left top corner of Fig. 1. The optical properties of the tissue types are ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ for the fatty region; ${\mu }_a=0.015{\mathrm{mm}}^{-1}$ for fibro-glandular; and tumor with ${\mu }_a=0.02{\mathrm{mm}}^{-1}$. The ${\mu }_s^{\prime }$ is assumed to be known and kept constant at $1{\mathrm{mm}}^{-1}$ throughout the imaging domain. The boundary measurements are taken using the 16 light source detector fiber-optic bundles arranged in an equispaced manner along the edge when one fiber delivers the light, then the rest act as detectors, which provides $M=240$ measurements. The mesh used for the data collection is a fine mesh with 4876 nodes (9567 triangular elements), and 1%, 5%, and 10% Gaussian noise was added to mimic the experimental case. For the reconstruction, a coarse mesh with 1969 nodes (3753 triangular elements) was used. The reconstruction results using fatty region as initial guess is given in Fig. 1, except on the top row where the guess of $0.001{\mathrm{mm}}^{-1}$, data noise level of 1% used which is far away from the actual solution. The details of the computational cost for the case with 1% noise are noted as follows. The traditional method converges to solution in three iterations (1.9 s) for the initial guess, being close to the actual solution (second row) with $\lambda =0.01$. If $\lambda =100$ is used in this case, the convergence is achieved in 11 iterations, taking 15.5 s. For the proposed simplex method, total computation time is 2.7 s (47 function evaluations). For the top-row results, the computation time for the traditional LM method is 3.5 s with $\lambda =0.01$ and 28.6 s with $\lambda =100$. The computational time is 12.8 s for the proposed simplex method. For the cases of 5% and 10% noise, the trend in the observed computational time is the same as for the 1% case. Figure 1 results were obtained using $\lambda =100$. The results obtained using the proposed simplex method match within 2% of results that were obtained using the traditional LM method.

Next, experimental data obtained from a cylindrical-shaped three-layered gelatin phantom that mimics the breast structure, of height 25 mm and diameter 86 mm, was used to test the effectiveness of the proposed method. The details of the gelatin phantom are given in Ref. 2 (Sec. 3.3). The reconstructed distributions of ${\mu }_a$ using both traditional LM and the proposed simplex methods are given in Fig. 2 along with the one-dimensional cross-sectional profiles. The number of iterations required for the traditional method was 6 for $\lambda =1$, requiring 3.1 s. For $\lambda =100$, the convergence was achieved in 11 iterations, requiring 14.9 s. The proposed method computational time is 9.5 s (180 function evaluations). This case also shows that the proposed method results are in close agreement with the ones obtained using the traditional LM method.

In both numerical and gelatin-phantom experiments, the total computational time that was reported for the traditional LM method was dependent on $\lambda$. The converged solution is the same for cases when $\lambda$ lies between 0.01 and 100; outside this range, it is far away from the actual solution. Note that there are automated methods for choosing the optimal $\lambda$, but they result in an overhead time that is almost equivalent to performing the reconstruction scheme, making it computationally inefficient.¹³ The proposed method does not require regularization, in turn removing the unwarranted bias introduced in the solution in the traditional LM method. The proposed method solves the minimization problem in an iterative manner, where the number of iterations acts as indirect regularization of the solution space (similar to least-squares QR¹⁴). As the number of iterations is dependent on the stopping criterion, the choice of the same becomes important in this case. More importantly, the proposed method does not require any Jacobian, which might be a time-consuming process even for the image-guided DOT where the detection mechanism is based on CCD cameras resulting in $M=1e7$.¹⁵

In conclusion, we have proposed a new approach for image-guided diffuse optical tomographic image reconstruction that is based on the popular gradient-free Nelder–Mead simplex method. The reconstructed optical property images closely match the ones obtained using traditional image reconstruction methods. The developed algorithm for image-guided DOT is provided as an open source¹⁶ for interested readers.