1. Introduction

1.1. Generic problem and challenges

Precise and robust estimation of the optical properties of biological tissues is a key issue in the development of spectroscopic and imaging methods and in their applicability to non invasive diagnostics in clinics. Promising UltraViolet - Visible - Near InfraRed spectroscopy and imaging techniques have led to significant results in terms of in vivo biological tissue characterization and diagnosis [1]. Among these methods, Spatially-resolved Diffuse Reflectance Spectroscopy (SDRS) has been investigated for about two decades and proved its potential as a valuable optical biopsy tool [2–11]. The principle of SDRS, as shown in Fig. 1(a), is to place a multi-fiber optical probe in gentle contact with a biological tissue. The excitation fiber of the probe conducts the light of a wide wavelength range source I₀(λ) down to the tissue inside which it is absorbed and multiply scattered. Several receiving fibers collect the backscattered photons at known Source-Detector Separations (SDS) at the probe’s tip. Those collected photons are then converted by a spectrometer into measured (experimental) reflectance intensity spectra I_exp(λ,r) as a function of wavelength λ = [λ₁, .., λ_k, .., λ_K] (vector of a K finite number of sampled wavelengths depending on the spectral range and resolution of the system) and distance r = [r₁, .., r_m, .., r_M] (vector of M number of sampled SDS depending on the probe features). Taking into account the complexity of biological tissues, like sketched in Fig. 1(b), the in vivo quantification at various depths of their optical properties, such as absorption coefficients µ_a(λ), scattering coefficients µ_s(λ) and anisotropy factors g(λ), is intended to provide new key features for characterizing pathological changes arising before evolving towards malignant lesions [10, 12, 13]. The estimation of the optical properties of biological tissues consists of solving an inverse problem which is classically ill-posed (different sets of parameters may lead to closely similar results) and ill-conditioned (numerical instability of the solutions due to finite precision and errors in the data). In the inverse problem framework illustrated in Fig. 1(c), firstly, an analytical or a numerical modeling of the physical phenomena of light-tissue interaction has to be chosen to describe the forward problem appropriately. The latter provides the modeled reflectance intensity spectra I_mod(λ,r;p(λ)) simulated with the tissue optical parameter set p(λ) which can be for instance $\mathbf{p}(\lambda )=\left[{\mu }_s^{1\dots l\dots L}(\lambda ),{\mu }_a^{1\dots l\dots L}(\lambda ),g^{1\dots l\dots L}(\lambda )\right]$, where the superscript 1…l…L indicates the optical parameters of layers 1 to L and where l and L stand for the tissue layer index and the total number of layers respectively, as represented in Fig. 1(b). Secondly, an optimization procedure is applied with the aim of maximizing or minimizing an objective function (also called energy- or cost-function, usually quantifying a distance between experimental data and modeled data) by iteratively varying the values of the parameter set. In Fig. 1(c), the objective function is based on a measure quantifying the similarity between measured intensities I_exp(λ,r) and intensities I_mod(λ,r) simulated with the forward model. The choice of the initial conditions (initial set of p(λ) values) of the method to optimize the objective function has a direct impact on the precision and robustness of the overall estimation scheme. Finally, after several iterations, the output of the algorithm is presented as the optimized optical parameter set $\tilde{\mathbf{p}}(\lambda )=\left[{\tilde{\mu }}_s^{1\dots L}(\lambda ),{\tilde{\mu }}_a^{1\dots L}(\lambda ),{\tilde{g}}^{1\dots L}(\lambda )\right]$ of the biological tissue under investigation.

 

Fig. 1 General context of the work. (a) Simplified scheme of the SDRS system including main parts: light source, probe, tissue and multichannel spectrometer; (b) Schematic representation of the correspondence between experimental (tissue), simulation (numerical model) and physical model configurations; (c) Inverse problem general scheme corresponding to optical parameter estimation based on SDRS measurements.

Download Full Size | PDF

1.2. Previous works on estimation of tissue optical parameters based on SDRS

Over the past two decades, different approaches of light-tissue interaction modeling were proposed, namely based on the analytical solution of the Diffusion Approximation (DA), on semi-empirical models, on Monte Carlo (MC) simulations or on photon migration and look-up table-based models [14, 15]. In comparison to other methods, the MC simulations for light propagation in biological tissues give precise results with acceptable statistical noise level depending on the number of photons launched and the complexity of the model [6, 16, 17]. However its application in the frame of the aforementioned inverse problem to find the exact values of tissue optical parameters gives high error rates, notably due to the stochastic nature of MC simulations [18]. Consequently, the difference between experimental and modeled data (error-function to be minimized) remains large and affects strongly the robustness and precision of the optimization process due to inaccurate convergence. Thus, the first main drawback of the existing solutions to inverse experimental modeling which use stochastic methods as forward model is a probable convergence towards local minima instead of a global one.

Several recently developed approaches allow for solving this inverse problem of optical spectroscopy diagnosis with lower error in optical property estimation. A large number of these solutions implement the Levenberg-Marquardt Algorithm (LMA) well-known for solving nonlinear least-squares problems. However, this second-order optimization method needs first (gradient) and second (Hessian) derivatives to be calculated which is a strong drawback when used in an inverse problem where the forward model is based on a stochastic, i.e. noisy and computation time consuming, model [18, 19].

The look-up table-based method recently proposed by Tunnell et al. [20, 21] allows for determining values of optical coefficients on bi-layer phantoms by searching in an extensive database of intensity spectra (built by performing MC simulations for a large range of (µ_s(λ), µ_a(λ)) pairs), the spectrum which gives the lowest cost-function value, i.e. the closest spectrum to the one measured on the phantom. The problem was solved by considering a fixed value for the reduced scattering coefficient ${{\mu }^{\prime }}_s={\mu }_s(1-g)$ and by estimating the values of the top layer absorption coefficient ( ${\mu }_a^1(\lambda )$), i.e. index 1 indicating the first layer. The optimal value for the latter was obtained with an error of 10% compared to the “ground truth” (exact known) value. For both layers, the reduced scattering coefficient was determined with an average error of 15 % and the absorption coefficient estimation error for the bottom layer ranged from 12 to 25 % depending on the top layer thickness (highest error for the maximum thickness of 550 µm). The optimization scheme was based on the implementation of the LMA. The following mean squared percent error between the measured spectra I_exp(λ,r) and the modeled one I_mod(λ,r;p(λ)), was used as a cost-function f(I_exp,I_mod) consisting only of a data-term:

In Eq. (1)M = 2 (corresponding to SDS = 370 and 740 µm) and p(λ) = [µ_a(λ)]. The other parameters of the model were fixed and equal for both layers: ${\mu }_s(\lambda )=15{\frac{\lambda }{630}}^{-2}c{m}^{-1}$, g = 0.85, refractive index n = 1.4 and top layer variable thickness z ∈ [0–550]µm.

The research group of Ramanujam et al. [22, 23] also used a MC technique for the modeling of a turbid media with single and multiple fluorophores over a wide wavelength range using the LMA based optimization. They applied the following quadratic error function (see Eq. (2)) over wavelengths λ to determine the value of the thickness z of the top layer , as well as those of the scattering (µ_s) and absorption (µ_a) coefficients of the bottom layer:

Other recent approaches on that inverse problem solving proposed the use of one or several neural networks to determine the values of the optical parameters for two-layered tissue models [24, 25]. The use of several neural networks leads to errors in the thickness estimation from 12 to 42% depending on the absorption rate: the higher it is, the higher is the estimated error in values of optical properties. The cost-function for this method is given in Eq. (2).

In 2012, another work of Fredrikson et al. [19] used an energy function not only based on the difference between modeled and experimental SDR intensity spectra, but also taking into account relative differences according to SDS by assigning lower weighting factors σ_m(λ) to longer SDS:

In Eq. (3), σ_m is influenced both by absolute and by relative differences of intensities. $\varphi (I)=\frac{{\overline{I}}_{\lambda \in [535,540],[575,580]}}{{\overline{I}}_{\lambda \in [550,570]}}$ is an additional regularization term (rightmost in this equation) adding importance to the shape of the spectrum by taking into account the double-peak of oxygenated hemoglobin through their mean intensity values Ī. A_m are the weighting coefficient values applied as a function of the different distances m. Here, A_m = 2 and 0.5 for m = 1 and m = M = 2 respectively.

The principal difference between the above mentioned cost-functions is that the one given in Eq. (1) takes into account the intensities at different SDS with the same weight, i.e. relative difference, while Eqs. (2) and (3) use either no weight or different weights for the respective different SDS. These energy functions thereby calculate absolute difference between measured and modeled spectra. For forward models based on a stochastic numerical simulation method, the statistical noise level, which increases with the SDS value, affects the cost-function value at larger SDS.

1.3. Aim and objectives of the present contribution

This contribution proposes a new approach for drastically reducing the error between estimated and ground truth values of optical tissue parameters. To do so, a multidimensional approach is used to solve the inverse problem step in the SDR spectroscopy approach sketched in Fig. 1(c). The first input of the corresponding optimization problem corresponds to the data measurement configuration with our experimental set-up. The latter is described in section 2 which indicates the set-up parameters involved in the optimization process described in this paper. Section 3 details all aspects of the new method proposed. Subsection 3.1 first describes and justifies the chosen forward model providing the second input (simulated data with given µ_s and µ_a parameter values) to the optimization process (the first input being the acquired or ground truth data). To clearly explain our optimization approach, the whole algorithm is gradually explained in subsections 3.2, 3.3 and 3.4. The subsection 3.2 briefly justifies the use of a Particle Swarm Optimization (PSO) method and describes its features and first implementation for spatially resolved data, i.e. for all available spatial point (or SDS distances), but for a given wavelength λ. Section 3.3 shows how the algorithm is modified and trained to take into account the Wavelength-Resolved (WR) data, i.e. both the whole discrete wavelength range and all SDS distances. In section 3.4, this WR-PSO approach is then refined using a shape constraint on the wavelength curves of the parameters µ_s(λ) and µ_a(λ) to be estimated. This final optimization method, called FWR-PSO, is again trained to evaluate the algorithm convergence performance. In order to reduce the computing time during the training procedure, while still taking into account the statistical noise features of the MC based simulations, the combination of the Diffusion Approximation (DA) equation together with white Gaussian noise according to the number of photons is proposed in subsection 3.1. Section 4 presents the results of the experimental validation performed on numerical data, obtained with MC simulations performed for a configuration of multi-fiber probe used in our clinical and experimental research works [26–28]. The significance and scope of these results in terms of parameter estimation precision, number of iterations along the optimization process and overall computational cost are discussed.

2. Characteristics of the experimental set-up for numerical model

The parameters of the present numerical study are based on the geometrical configuration and features of our experimental SDR spectroscopy setup as implemented in previous studies (see details in [18, 28, 30]). Briefly, the white light source is a halogen lamp with maximum temperature of 3700 K and the spectrometer is tuned for the wavelength range [400–900] nm with resolution of 1 nm. For our numerical study, a wavelength range of [500–600] nm with 10 nm step (however the results could be applied to any wavelength range in optical region of interest).

As sketched on Fig. 2, our multi-fiber probe has M = 6 different SDS with distance values r_1..6 = [261,344,500,778,1041,1290] µm, with each fiber having a core refractive index n = 1.54 and a diameter d = 200 µm [28]. The choice of the appropriate ranges of the optical parameters to be adjusted is discussed in subsection 3.3.

 

Fig. 2 Geometrical configuration of the experimental set-up used for the numerical simulations. The excitation light source illuminates the tissue through fiber F₀ of the multi-fiber probe [26] with tissue incident light intensity I₀(λ). SDR intensity spectra I_exp(λ,r) are collected at 6 different SDS corresponding to distances r_1…6 = [261,344,500,778,1041,1290] µm, respectively.

Download Full Size | PDF

3. Methods and solutions proposed for solving the inverse problem

3.1. Forward problem modeling

As mentioned in section 1.1, estimating tissue optical properties is a complex task because of the large dimension of the parameter space (µ_a(λ), µ_s(λ)) which requires a time-expensive iterative procedure. Therefore, the amount of data related to this large parameter space has to be taken into account during both the forward and the inverse modeling steps, so that the global approach presented in Fig. 1(c) leads to reasonable processing times. As also mentioned in subsection 1.2, the MC method is classically used as the reference approach for simulating the forward model. However, the computation time and the statistical noise associated to this simulation approach are main drawbacks in our complex optimization stated above. To address that issue, the modeled spatially resolved DR intensity spectra I_mod(λ,r) will be obtained by the Diffusion Approximation (DA) equation defined in Eq. (4).

This analytical expression can be fastly calculated. To make this analytical calculation closer to the MC simulation results in terms of statistical noise and SDS dependence, our idea was to add similar features of noise. To characterize these noise rates and SDS dependence (according to the number of photons launched), we implemented a CUDA-based parallel version of MC (Intel(R) Xeon(R)CPU E5-2609, 2.4 Ghz, RAM 16 Gb and Video Card - NVidia Tesla C2075), whose basic algorithm description could be found elsewhere [6, 31–33]. For each of the following number of photons, N_f = 10⁷, 10⁸, 10⁹, MC simulations were run 10 times for the configuration given in Fig. 2. The obtained probability distributions were found as normal (Gaussian) and the corresponding standard deviations (SD, corresponding to “noise level” in %) are shown in Fig. 3 as a function of SDS(r), and N_f.

 

Fig. 3 MC simulation noise level (standard deviation in %) as a function of SDS (r) and photon number N_f launched : • for N_f = 10⁷, ▲ for N_f = 10⁸ and for N_f = 10⁹. The continuous, dashed and dotted curves drawn are respective exponential fitting (with corresponding equations listed in the legend) for each N_f.

Download Full Size | PDF

In the proposed approaches described in the next sections, anisotropy factor and refractive index values were fixed at g = 0.8 and n = 1.44 respectively and all the simulations were performed for a one layer numerical model with thickness z = 5 cm. Scattering and absorption coefficient ground-truth spectral functions ${\mu }_s^{ref}(\lambda )$ and ${\mu }_a^{ref}(\lambda )$ used in our model are presented in Fig. 4.

 

Fig. 4 Scattering ${\mu }_s^{ref}(\lambda )$ and absorption ${\mu }_a^{ref}(\lambda )$ coefficients defined as the ground-truth spectral data in the estimation algorithms using WR-PSO (section 3.3) and FWR-PSO (section 3.4) approaches.

Download Full Size | PDF

3.2. PSO-based algorithm applied to spatially resolved data

As mentioned in subsections 1.1 and 3.1, the large parameter space and the noise features of the forward model (depending on the SDS) result in a complex inverse problem to solve with one global minimum and many local minima. Metaheuristic, are now known as efficient approach for such hard optimization problem [36]. Metaheuristics are “global and stochastic” optimization algorithms which demonstrate to efficiently face mutidimensional problems. Among the population-based metaheuristics, most bioinspired algorithms belong to the classes of Evolutionary Algorithms and Particle Swarm Optimization (PSO) [37]. PSO is a swarm-based evolutionary algorithm introduced by Eberhart and Kennedy and used for optimization of continuous non-linear functions [38–41] which leads to faster solution in complex multidimensional problems [18]. A first promising advantage of this algorithm is that it doesn’t require derivative calculation, thus reducing computation complexity and cost. Another advantage lies in the directivity of the algorithm, that is its ability to choose the direction of the next iteration with new starting points based on the previous experience [42]. Its limitations are in the final convergence rate (convergence when close to the end) and the difficulty to adjust the values of the optimization method parameters which are usually tuned empirically. The question of the theoretical analysis of the convergence of PSO is still under investigation especially for multidimensional spaces [36, 43]. To empirically explore the convergence performance in our present study, the results of the PSO-based algorithms proposed in subsections 3.3 and 3.4 are compared in terms of success rate (number of converging results divided by the total number of results) and number of iterations for reaching successful convergence.

The PSO algorithm consists of initializing the optimization process with a population of random solutions, called particles, and searching for optima by updating generations while progressing throughout the search space of probable solutions. PSO stores the coordinates of each particle related to the personal best solution obtained as far by each particle in its cluster or group (personal best value or grBest for “group Best”) and the global best value among all the particles (global best value or gBest). At each iteration step, the acceleration of each particle towards grBest and gBest locations is changed by respective weighting random terms.

The subsequent part of this section describes the basic solving scheme for our inverse problem based on a PSO algorithm adapted for taking into account the spatial resolution of the spectroscopic data (spectral dimension being added in the next subsection). In our case, where µ_a and µ_s are the two parameters to be estimated, a particle corresponds to a point in the (µ_a, µ_s) parameter space of the cost-function. For the dimension D = 2 of the parameter space, our approach considered C = (D+1)² = 9 clusters of P = D+1 = 3 particles. In the rest of the paper, any cluster will be denoted c, while cj refers to a specific cluster with j = [1 : C] the cluster index. In a similar way, any particle will be denoted p, while a specific particle in a cluster c_j will be denoted p_i,j with i = [1 : P] the particle index in a cluster. At start of the optimization process (initial conditions), these 9 clusters were spread all over the search space and the 3 particles in each cluster were placed close to each other with less than 5% difference in their coordinates X = (µ_a, µ_s) in the space of variables. Given the data found in literature on epithelial tissue optical properties [8, 14, 35, 44, 45], the [lower limit:upper limit] range chosen for the solution space of µ_s and µ_a values to be estimated were [100:1000] cm⁻¹ and [1:10] cm⁻¹, respectively. The “error for a cluster c_j”, at one single wavelength λ_k, corresponds to the averaged value of the cost-function ${\chi }_j^2\left({\lambda }_k\right)$ calculated for the 3 particles in this cluster (see Eq. (5)):

Moreover, the optical parameters are averaged inside each cluster thus giving, for each cluster, $grBes{t}_j=\left[{\chi }_j^2,{\mu }_s^j,{\mu }_a^j\right]$. In the C × P-dimensional space, gBest is related to the minimal value of ${\chi }_j^2$ among all the clusters, i.e. gBest = min_j(grBest_j). The use of gBest (instead of the individual cluster values grBest_j) contributes to limit the risk of convergence towards local minima. Thus, after calculating the cost-function ${\chi }_j^2$, all the c clusters of particles are moved towards the one which has the global “best” properties ( ${\mu }_s^{gBest}$, ${\mu }_a^{gBest}$), and each cluster towards its averaged “best” properties ( ${\mu }_s^{grBest,j}$ and ${\mu }_a^{grBest,j}$.

Then, going from iteration iter to iter + 1, each particle with coordinates X^iter( ${\mu }_s^{iter}$, ${\mu }_a^{iter}$) is moved with the speed V^iter+1, to the next position X^iter+1 in the space of variables (µ_s, µ_a). Their updated speed V^iter+i and coordinate X^iter+i are defined by Eq. (6):

For each possible combination of these values, the estimated values ${\tilde{\mu }}_a$ and ${\tilde{\mu }}_s$ were obtained from the PSO-algorithm and compared to the ground truth values ${\mu }_a^{ref}$ and ${\mu }_s^{ref}$, as summarized in Algorithm 1. The best parameter value combination w_min, w_max, a₁, a₂, among all possible combinations of the above, was that leading to ${\tilde{\mu }}_a$ and ${\tilde{\mu }}_s$ which are simultaneously the closest to ${\mu }_a^{ref}$ and ${\mu }_s^{ref}$ respectively, and for iteration iter for which the ${\chi }^2=min\left({\chi }_j^2\right)$ is smaller than a given value of PSO stop criteria ∆ (see Algorithm 1). In our case, convergence was supposed to be reached in less than maxIter = 30 iterations, when χ² is smaller than optimization stop criterion ∆ with values between 5% and 10%. The latter were chosen from the MC noise level results obtained in the subsection 3.1 and shown in Fig. 3 where the error reaches 10 % at the largest SDS and for the N_f = 10⁷ and is lower than 5% for N_f = 10⁸. Therefore, Algorithm 1 was also run for the range of ∆ values from 5% to 10%.

Algorithm 1. PSO-based algorithm applied to Spatially Resolved (SR) data

View Table | View all tables in this article

In practice, for the first iteration (iter = 1), the parameter space was divided into rectangles with lengths $L_{{\mu }_s}$ and $L_{{\mu }_a}$ of 90 cm⁻¹ and 9 cm⁻¹, respectively, the total number of squares being 100. However the resolution of the derived final values of scattering and absorption coefficients does not depend on such preliminary division. Starting with (D + 1)² = 9 clusters of particles, the weight w is reduced during the first iteration steps to make each cluster “choose” the square of local search. Thereby after 10 iterations, it was decided to keep weight w constant as the change in weight from maximum to minimum depends on the number of maximum iterations (see Eqs. (6)). Then, with the increasing number of iteration iter, the transformation of global search was slowed down to local, thus thresIter and maxIter were fixed to the values of 10 and 30 respectively (see Eq. (6)).

3.3. Wavelength Resolved (WR)-PSO approach and PSO parameters adjustment

The aim of this section is to estimate the Wavelength Resolved (WR) values of absorption µ_a(λ) and scattering µ_s(λ) coefficients by extending the aforementioned PSO-based inverse problem algorithm on spatially resolved data to spectral resolution, both parameters being now K-dimensional vectors instead of scalars. In practice, the best compromise in the values of the PSO parameters w_min, w_max, a₁, a₂ and ∆ has to be determined together with the optimization of µ_a(λ) and µ_a(λ) with λ a K = 10 wavelength value vector such as λ = [λ₁,⋯ ,λ_k,⋯ ,λ_K] = [500,510,⋯ ,580,590] nm.

The principle of the Wavelength-Resolved(WR)-PSO method estimating the spectral values of µ_a(λ) and µ_s(λ) is described in Algorithm 2. The latter also includes the systematic search implemented to determine the best combination of values of the PSO parameters a₁, a₂, w_min and w_max (which gave the best estimation of µ_a(λ) and µ_s(λ)) among a large set of values given in Eq. (7).

Algorithm 2. Wavelength-Resolved (WR)-PSO based algorithm

View Table | View all tables in this article

The main steps of the algorithm are the following:

The results are presented in Table 1. Columns 1–4: Selection of 9 combinations of the PSO parameters (particles acceleration parameters a₁ and a₂ towards global gBest and cluster grBest bests, and initial w_Max and final w_Min weights). Columns 5–10: Corresponding performance criteria obtained (Simulation Success Rate SSR, maximum I_max and average < I > number of iterations, errors ${\epsilon }_{{\mu }_s}$ and ${\epsilon }_{{\mu }_a}$) in µ_s(λ) and µ_a(λ) estimations.

Table 1. Results obtained with the WR-PSO algorithm 2, page 13.

View Table | View all tables in this article

It can be observed that the lowest errors (in estimating the scattering and absorption coefficients) ${\epsilon }_{{\mu }_s}=10\%$ and ${\epsilon }_{{\mu }_a}=9\%$ were achieved for the following PSO parameter value combination: a₁ = 0.1, a₂ = 0.6, w_min = 0.4, w_max = 0.7 and for χ² = 7%. However, even if the errors in optical property estimation are comparatively small, it can also be observed that the algorithm converges only once out of five runs of the WR-PSO algorithm. Further refinements are therefore proposed through a complementary approach of the WR-PSO algorithm extension in the subsection 3.4.

3.4. Modified WR PSO-based approach with spectral fitting

To further improve the convergence of the aforementioned WR-PSO-based algorithm, a modified approach of Algorithm 2 was proposed through a “spectral Fitting” strategy (FWR-PSO). Our idea was to take advantage of the spectral dimension of the data to bring a new coherence constraint into the choice of scattering µ_s(λ) and absorption µ_a(λ) values during the optimization process. To do so, an a priori spectral shape was considered for each optical parameter. Then, at each iteration step, each parameter value that is calculated independently at each wavelength is slightly modified to fit to the a priori spectral shape. The spectral shape of the optical parameters to fit µ_s(λ) and µ_a(λ) coefficients was chosen to be exponential, as drawn in Fig. 4. The latter is a common shape for biological tissue optical properties in the UV-visible range such as absorption of melanin and the scattering by small particles according to Rayleigh and Mie theories [14, 46–49]. This FWR-PSO strategy was implemented in Algorithm 3 hereafter in which the fitting procedure (in the rectangular box) was added to the previous WR-PSO Algorithm 2. After fitting ${\mu }_s^{gBest}(\lambda )$ and ${\mu }_a^{gBest}(\lambda )$ with the aforementioned exponential functions ${\mu }_s^{fit}(\lambda )=a_{{\mu }_s}\times e^{-b_{{\mu }_s}\lambda }$ and ${\mu }_a^{fit}(\lambda )=a_{{\mu }_a}\times e^{-b_{{\mu }_a}\lambda }$, a “ ${\chi }_{fit}^2$” value is calculated using Eq. (8) with ${\mu }_s^{fit}(\lambda )$ and ${\mu }_a^{fit}(\lambda )$. Compared to the χ² value obtained in the current iteration iter, if the fitting cost-function is smaller ( ${\chi }_{fit}^2<{\chi }^2$) i.e. the convergence is better, then, X^iter, gBest(λ) and grBest_j(λ) are updated with the fitted values ${\mu }_s^{fit}(\lambda )$ and ${\mu }_a^{fit}(\lambda )$.

Algorithm 3. Fitting modified (FWR)-PSO algorithm

View Table | View all tables in this article

As described in subsection 3.3, the best combinations in the values of the PSO parameters (for the ranges given in Eq. (7)) were also determined together with the optimization of µ_s(λ) and µ_a(λ). To be compared to the WR-PSO algorithm, the same set of performance criteria was calculated for FWR-PSO: the average < I > and maximum I_max number of iterations before convergence, the Simulation Success Rate (SSR) and the mean errors (%) ${\epsilon }_{{\mu }_s}$ and ${\epsilon }_{{\mu }_a}$ obtained for the final estimates ${\tilde{\mu }}_s\left(\lambda \right)$ and ${\tilde{\mu }}_a\left(\lambda \right)$ respectively, calculated by Eq. (9). These values are given and discussed in the next section.

4. Results and discussion

Besides the PSO parameters a₁, b₁, w_min, w_max and ∆ whose combinations (within predefined ranges) were systematically tested, other parameters had to be adjusted to find an appropriate balance between accuracy and computation amount for both WR-PSO (see Section 3.3) and FWR-PSO (see Section 3.4) approaches.

The first parameter was the number of the iterations which relates to the time required to find a solution for the inverse problem. In Fig. 5, < I > and I_max are the mean number of iterations and the maximal number of iterations respectively observed when testing the WR-PSO and the FWR-PSO algorithms. For the WR-PSO method, < I > values are in the range of 10–15] iterations whatever χ², whereas for the FWR-PSO method, < I > clearly decreases with increasing χ², with a mean iteration number lower than 5 for χ² greater than 7%. For the FWR-PSO algorithm, accepting a χ² error ranging in [5, 8]% leads to a mean number of iterations belonging to [5–10]. I_max(FWR-PSO) is also decreasing compared to I_max(WR-PSO) which increases but with a lower number of maximum number of iterations for χ² values below 7%.

 

Fig. 5 Impact of the average < I > and maximum (I_max) number of iterations on the χ² error for the WR-PSO and FWR-PSO-based approaches. Solid curves correspond to fitting functions whose meaning is indicated by their different associated symbols.

Download Full Size | PDF

Figure 6 gives, according to the χ² error, the percentage of tries (referred to as SSR) which enabled the optimization algorithm (WR-PSO or FWR-PSO) to reach a given accuracy (i.e. maximal χ² error) without exceeding a given maximum number of iterations (here maxIter = 30). It can be seen on Fig. 6 that when χ² decreases from 10% to 5%, the SSR diminishes strongly for both WR-PSO and FWR-PSO approaches. However, the SSR value of the FWR-PSO method is systematically better (about 30% more) to the one of the WR-PSO method. For instance, for χ² = 5%, SSR(WR-PSO) is below 5% while SSR(FWR-PSO) is still around 40%.

 

Fig. 6 Impact of the χ² value on the WR-PSO and FWR-PSO convergence success rates (RSS in %).

Download Full Size | PDF

Figure 7 gives the errors ( ${\epsilon }_{{\mu }_s}$ and ${\epsilon }_{{\mu }_a}$) defined in Eq. (9) for the scattering and absorption coefficients which allow for quantifying the performances of the WR-PSO and FWR-PSO algorithms according to the measured χ² errors. As expected, estimation errors for both µ_a and µ_s increase with χ². The errors on µ_a and on µ_s obtained with FWR-PSO are significantly lower than the ones obtained with WR-PSO. The lowest ${\epsilon }_{{\mu }_s}$(FWR-PSO) and ${\epsilon }_{{\mu }_a}$(FWR-PSO), corresponding to χ² = 5%, are below 5% while the lowest ${\epsilon }_{{\mu }_s}$(WR-PSO) and ${\epsilon }_{{\mu }_a}$(WR-PSO) are above 10%.

 

Fig. 7 Dependence on χ² value of estimated optical parameters error in scattering ( ${\epsilon }_{{\mu }_s}$, %) and absorption ( ${\epsilon }_{{\mu }_a}$, %) values during WR-PSO (section 3.3) and FWR-PSO (section 3.4) algorithms convergence

Download Full Size | PDF

As it is seen from the Fig. 7, the error between modeled and experimental data exhibits a linear relationship with the error in optical property estimation, i.e. χ² ∼ ε(p_i). Thus it could be concluded that error between modeled and reference data should be [5, 7] % to have an error in optical properties less than 10%. According to the average < I > and maximum number I_max of iterations, a χ² = 6% error was chosen to reduce computational time (for this accuracy, the rate of success is 60 % and the average number of iterations until the FWR-PSO algorithm converges is 3). More precisely, for 3 iterations, the errors in the optical parameter estimation were ${\epsilon }_{{\mu }_s}=6.2\%$ and ${\epsilon }_{{\mu }_a}=6.6\%$. For these results (with χ² = 6%), the FWR-PSO-based algorithm parameters were: a₁ = 1.6, a₂ = 1.2, w_Min = 0.1 and w_Max = 0.5 (see Table 2). The table is organized as the table 1 on page 14.

Table 2. Results obtained with the FWR-PSO algorithm 3 page 15.

View Table | View all tables in this article

The low success rate of 48 % with high average number of simulations (meanIter = 29), which increases computational time in a drastic way when MC simulations are used as forward model, give lower error (less than 5 %) for both scattering (µ_s) and absorption (µ_s), and with PSO parameter values: a₁ = 0.6, a₂ = 1.7, w_Min = 0.8 and w_Max = 1.3 (χ² = 6 %).

When comparing the results of the two proposed algorithms, it can be noted that the FWR-PSO-based approach has a rate of successful convergence in average three times higher than the one of the WR-PSO-based approach. In a similar way, the ${\epsilon }_{{\mu }_a}$ and ${\epsilon }_{{\mu }_s}$ errors computed for the optical parameter µ_s and µ_a respectively, decreased by a factor of 2 and 2.4 when passing from the WR-PSO-based approach to the improved FWR PSO-based approach. The comparison of the best results given in Fig. 7 for both approaches, shows that the minimum values of errors on µ_s and µ_a were respectively ${\epsilon }_{{\mu }_s}=15\%$ and ${\epsilon }_{{\mu }_a}=12\%$ for the WR-PSO approach, and reduced to ${\epsilon }_{{\mu }_s}=6\%$ and ${\epsilon }_{{\mu }_a}=5\%$ for the FWR-PSO approach which includes the fitting method. When passing from the WR-PSO to FWR-PSO algorithm, the average number of iterations (< I >) was also decreased from 14 down to 6 and 3 (see Table 1 and Table 2). It is obvious that the mean number of iterations is in direct relation with the duration of the simulations. In particular, when the MC method is used in the simulations with the forward model, the computation time of the inverse problem till convergence would be decreased by a factor of 3.

5. Conclusion

The FWR-PSO method allows for determining wavelength depending optical properties. For the proposed noisy DA method, the fitting based FWR-PSO algorithm had been trained for each particular shape of the “fitting curve” chosen to increase the accuracy and speed of convergence. The method presented allows for reducing the number of iterations and contributes to an accurate convergence. Moreover, the proposed optimization algorithm robustness (convergence towards the global minimum) is ensured since local minima are avoided due to several reasons: i) at each iteration, the area in the parameter space, covered by the particles of a cluster, diminishes with the decreasing values of weight w according to acceleration a₂, ii) the area, covered by each of the clusters, also reduces with the decreasing weight w, but with the other acceleration coefficient a₁, and iii) the averaging inside each cluster of the optical property parameters stored in the vectors p(λ) associated to the particles. The use of two best solutions from the whole spectrum of optical properties allows for reducing the number of steps of the optimization process. Indeed, the fit of the curve (exponential decay or other curves) to this best two solutions forces the particles to move towards the global minimum. It is noticeable that the proposed algorithm can be extended to wavelength dependent anisotropy factors g(λ) and refractive indexes n(λ), as well as to multi-layer models with different layer thicknesses. The constraint imposed by the fitting function lead to results with an high accuracy: the optical property parameter estimation is affected by an error of less than 5%. To avoid convergence problems, further analytical analysis is needed: for instance, in the PhD thesis of Schmitt [43], an appropriate balance between the number of iterations and the error in optimized optical parameters was achieved by adjusting the number of clusters and number of particles in each cluster according to the dimension of the parameter space.

Finally, the proposed FWR-PSO method can be used with any forward model, under the condition that the direct model (simulation box in Fig. 1(c)) and the inverse problem (optimization box in Fig. 1(c)) deals with the same optical parameters since the particle and cluster speed intervals have to be chosen accordingly to the type of the parameters. As a perspective of this contribution, it is envisaged to train the algorithm for other wavelength depending optical parameters, such as absorption of oxy- and deoxyhaemoglobin with Protoporphyrin IX. A multi-layered model will also be studied in further work.