## Abstract

This paper presents a new approach to estimate optical properties (absorption and scattering coefficients *µ _{a}* and

*µ*) of biological tissues from spatially-resolved spectroscopy measurements. A Particle Swarm Optimization (PSO)-based algorithm was implemented and firstly modified to deal with spatial and spectral resolutions of the data, and to solve the corresponding inverse problem. Secondly, the optimization was improved by fitting exponential decays to the two best points among all clusters of the “particles” randomly distributed all over the parameter space (

_{s}*µ*) of possible solutions. The consequent acceleration of all the groups of particles to the “best” curve leads to significant error decrease in the optical property estimation. The study analyzes the estimated optical property error as a function of the various PSO parameter combinations, and several performance criteria such as the cost-function error and the number of iterations in the algorithms proposed. The final one led to error values between ground truth and estimated values of

_{s}, µ_{a}*µ*and

_{s}*µ*less than 6%.

_{a}© 2016 Optical Society of America

## 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*_{0}(*λ*) 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

*λ*= [

*λ*

_{1}, ..,

*λ*] (vector of a

_{k}, .., λ_{K}*K*finite number of sampled wavelengths depending on the spectral range and resolution of the system) and distance

*r*= [

*r*

_{1}, ..,

*r*] (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

_{m}, .., r_{M}*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 ),\phantom{\rule{0.2em}{0ex}}{\mu}_{a}^{1\dots l\dots L}(\lambda ),\phantom{\rule{0.2em}{0ex}}{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 ),\phantom{\rule{0.2em}{0ex}}{\tilde{\mu}}_{a}^{1\dots L}(\lambda ),\phantom{\rule{0.2em}{0ex}}{\tilde{g}}^{1\dots L}(\lambda )\right]$ of the biological tissue under investigation.

#### 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*) consisting only of a data-term:

_{exp},I_{mod}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{{\scriptscriptstyle \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 (

*µ*) coefficients of the bottom layer:

_{a}**p**(

*λ*) = [

_{k}*µ*(

_{s}*λ*),

_{k}*µ*(

_{a}*λ*),

_{k}*z*], while

*g*and

*n*are fixed values. The errors of the final values estimated were around 20% and 15% for

*µ*and

_{a}*µ*respectively. For this two-layered model, the top-layer thickness was defined with 20% error.

_{s}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*are the weighting coefficient values applied as a function of the different distances

_{m}*m*. Here,

*A*= 2 and 0.5 for

_{m}*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

*µ*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

_{a}*λ*. 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.

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

^{7}, 10

^{8}, 10

^{9}, 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}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.

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

*µ*are the two parameters to be estimated, a particle corresponds to a point in the (

_{s}*µ*) parameter space of the cost-function. For the dimension

_{a}, µ_{s}*D*= 2 of the parameter space, our approach considered

*C*= (

*D*+1)

^{2}= 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*will be denoted

_{j}*p*with

_{i,j}*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*= (

*µ*) 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

_{a}, µ_{s}*µ*and

_{s}*µ*values to be estimated were [100:1000] cm

_{a}^{−1}and [1:10] cm

^{−1}, respectively. The “error for a cluster

*c*”, at one single wavelength

_{j}*λ*, 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)):

_{k}**p**

*(*

_{i}*λ*) = [

_{k}*µ*(

_{s}*λ*),

_{k}*µ*(

_{a}*λ*)] the optical parameters to be estimated for the

_{k}*i*particle

_{th}*p*of cluster

_{i,j}*c*, and ${\mathbf{p}}^{ref}=\left[{\mu}_{s}^{ref}({\lambda}_{k}),{\mu}_{a}^{ref}({\lambda}_{k})\right]$ the ground-truth value parameter vector.

_{j}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*). The use of

_{j}*gBest*(instead of the individual cluster values

*grBest*) contributes to limit the risk of convergence towards local minima. Thus, after calculating the cost-function ${\chi}_{j}^{2}$, all the

_{j}*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 (

*µ*). Their updated speed

_{s}, µ_{a}*V*

^{iter}^{+}

*and coordinate*

^{i}*X*

^{iter}^{+}

*are defined by Eq. (6):*

^{i}*a*

_{1}and

*a*

_{2}are coefficients weighting the acceleration of the particles to the

*gBest*and

*grBest*, respectively. To accelerate each cluster simultaneously in the same direction, two random numbers are fixed inside each cluster. The random number generators

_{j}*rand*

_{1}and

*rand*

_{2}provide uniformly distributed values between 0 and 1 according to the values of (

*grBest*)and (

_{j}− X^{iter}*grBest − X*), respectively.

^{iter}*w*is a weighting function in which

*w*and

_{Max}*w*are the initial and the final weights, and

_{Min}*iter*the current iteration number running from 1 to

*maxIter*, the maximum number of iterations.

*thresIter*stands for a threshold number of iterations between global and local search strategies. During the first iterations (

*iter < thresIter*), the weighting function results in high

*w*values (high inertia) favoring exploration i.e. assisting to global search. For

*iter > thresIter*,

*w*decreases (low inertia) which favors exploitation i.e. forcing the algorithm to go from global to local search. It is worth to mention that the number of clusters was chosen as a trade-off between the computational cost during one iteration step and the number of iterations. Thus the higher is the number of clusters the longer is each iteration step and the higher is the probability of clusters to find the proper solution. The adjustment of the values of the PSO parameters

*w*,

_{Max}*w*,

_{Min}*a*

_{1},

*a*

_{2},

*thresIter*was empirically performed from the following ranges of values:

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*

_{1},

*a*

_{2}, 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

*χ*

^{2}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*= 10

_{f}^{7}and is lower than 5% for

*N*= 10

_{f}^{8}. Therefore, Algorithm 1 was also run for the range of ∆ values from 5% to 10%.

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 ^{−}*

^{1}and 9

*cm*

^{−}^{1}, 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)

^{2}= 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*

_{1},

*a*

_{2}and ∆ has to be determined together with the optimization of

*µ*(

_{a}*λ*) and

*µ*(

_{a}*λ*) with

*λ*a

*K*= 10 wavelength value vector such as

*λ*= [

*λ*

_{1},⋯ ,

*λ*] = [500,510,⋯ ,580,590] nm.

_{k},⋯ ,λ_{K}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*

_{1},

*a*

_{2},

*w*and

_{min}*w*(which gave the best estimation of

_{max}*µ*(

_{a}*λ*) and

*µ*(

_{s}*λ*)) among a large set of values given in Eq. (7).

The main steps of the algorithm are the following:

- The tissue optical model parameters are initialized as follows:
*g*= 0.8 and*n*= 1.44. The ranges of the tested values of PSO parameters*a*_{1},*a*_{2},*w*are defined in Eq. (7);_{min}w_{max} - The best PSO parameter combination, at each wavelength
*λ*, is chosen according to the minimum error defined in Eq. (8):_{k} *maxIter*= 30 iterations were run for each error calculation in a range*χ*^{2}= 5–10 % with 1% step;- The following performance criteria were calculated:
- The average number of iterations <
*I*> and the maximum number of iterations*I*before convergence is achieved;_{max} - The Simulation Success Rate (SSR) in percent defined as the number of successful simulations divided by
*maxIter*(an optimization was counted as successful when convergence was achieved before reaching*maxIter*); - The mean errors ${\epsilon}_{{\mu}_{s}}$ and ${\epsilon}_{{\mu}_{a}}$ in percent obtained for the final estimates ${\tilde{\mu}}_{s}\left(\lambda \right)$ and ${\tilde{\mu}}_{a}\left(\lambda \right)$ respectively, calculated as:

The results are presented in Table 1. Columns 1–4: Selection of 9 combinations of the PSO parameters (particles acceleration parameters *a*_{1} and *a*_{2} towards global *gBest* and cluster *grBest* bests, and initial *w _{Max}* and final

*w*weights). Columns 5–10: Corresponding performance criteria obtained (Simulation Success Rate SSR, maximum

_{Min}*I*and average <

_{max}*I*> number of iterations, errors ${\epsilon}_{{\mu}_{s}}$ and ${\epsilon}_{{\mu}_{a}}$) in

*µ*(

_{s}*λ*) and

*µ*(

_{a}*λ*) estimations.

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*_{1} = 0.1, *a*_{2} = 0.6, *w _{min}* = 0.4,

*w*= 0.7 and for

_{max}*χ*

^{2}= 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

*χ*

^{2}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 )$.

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

_{max}## 4. Results and discussion

Besides the PSO parameters *a*_{1}, *b*_{1}, *w _{min}*,

*w*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.

_{max}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

*χ*

^{2}, whereas for the FWR-PSO method, <

*I*> clearly decreases with increasing

*χ*

^{2}, with a mean iteration number lower than 5 for

*χ*

^{2}greater than 7%. For the FWR-PSO algorithm, accepting a

*χ*

^{2}error ranging in [5, 8]% leads to a mean number of iterations belonging to [5–10].

*I*(FWR-PSO) is also decreasing compared to

_{max}*I*(WR-PSO) which increases but with a lower number of maximum number of iterations for

_{max}*χ*

^{2}values below 7%.

Figure 6 gives, according to the *χ*^{2} 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 *χ*^{2} error) without exceeding a given maximum number of iterations (here *maxIter* = 30). It can be seen on Fig. 6 that when *χ*^{2} 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 *χ*^{2} = 5%, *SSR*(WR-PSO) is below 5% while *SSR*(FWR-PSO) is still around 40%.

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 *χ*^{2} errors. As expected, estimation errors for both *µ _{a}* and

*µ*increase with

_{s}*χ*

^{2}. The errors on

*µ*and on

_{a}*µ*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

_{s}*χ*

^{2}= 5%, are below 5% while the lowest ${\epsilon}_{{\mu}_{s}}$(WR-PSO) and ${\epsilon}_{{\mu}_{a}}$(WR-PSO) are above 10%.

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. *χ*^{2} ∼ *ε*(*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*of iterations, a

_{max}*χ*

^{2}= 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

*χ*

^{2}= 6%), the FWR-PSO-based algorithm parameters were:

*a*

_{1}= 1.6,

*a*

_{2}= 1.2,

*w*= 0.1 and

_{Min}*w*= 0.5 (see Table 2). The table is organized as the table 1 on page 14.

_{Max}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 (

*µ*), and with PSO parameter values:

_{s}*a*

_{1}= 0.6,

*a*

_{2}= 1.7,

*w*= 0.8 and

_{Min}*w*= 1.3 (

_{Max}*χ*

^{2}= 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

*µ*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

_{a}*µ*and

_{s}*µ*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 (<

_{a}*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*_{2}, ii) the area, covered by each of the clusters, also reduces with the decreasing weight *w*, but with the other acceleration coefficient *a*_{1}, 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.

## Acknowledgments

This work was conducted in the frame of a co-supervised PhD thesis between the Université de Lorraine (CRAN UMR 7039 UL-CNRS,Nancy, France) and the GPI RAS (Moscow, Russia). M. Kholodtsova would like to thank the Conseil Régional de Lorraine and the French Embassy in Moscow for the financial support of her Ph.D. grant.

## References and links

**1. **E. Drakaki, T. Vergou, C. Dessinioti, A. Stratigos, C. Salavastru, and C. Antoniou, “Spectroscopic methods for the photodiagnosis of nonmelanoma skin cancer,” J. Biomed. Opt. **18**(6), 1–10 (2013).

**2. **M. Canpolat and R. M. Judith, “High-angle scattering events strongly affect light collection in clinically relevant measurement geometries for light transport through tissue,” Phys. Med. Biol. **45**(5), 1127 (2000).

**3. **I. J. Bigio and S. G. Bown, “Spectroscopic sensing of cancer and cancer therapy–current status of translational research,” Cancer Biol. Ther. **3**(3), 259–267 (2004).

**4. **K. Badizadegan, V. Backman, C. Boone, C. Crum, R. Dasari, I. Georgakoudi, K. Keefe, K. Munger, S. Shapshay, E. Sheets, and M. Feld, “Spectroscopic diagnosis and imaging of invisible pre-cancer,” Faraday Discuss. **126**, 265–279 (2004).

**5. **B. W. Murphy, R. J. Webster, B. A. Turlach, C. J. Quirk, C. D. Clay, P. J. Heenan, and D. D. Sampson, “Toward the discrimination of early melanoma from common and dysplastic nevus using fiber optic diffuse reflectance spectroscopy,” J. Biomed. Opt. **10**(6), 064020 (2005).

**6. **D. Arifler, C. MacAulay, M. Follen, and R. Richards-Kortum, “Spatially resolved reflectance spectroscopy for diagnosis of cervical precancer: Monte Carlo modeling and comparison to clinical measurements,” J. Biomed. Opt. **11**(6), 1–16 (2006).

**7. **R. R. Allison and C. H. Sibata, “Photodiagnosis for cutaneous malignancy: A brief clinical and technical review,” Photodiag. Photodyn. Ther. **5**(4), 247–250 (2008).

**8. **E. Borisova, P. Troyanova, P. Pavlova, and L. Avramov, “Diagnostics of pigmented skin tumors based on laser-induced autofluorescence and diffuse reflectance spectroscopy,” Quantum Electron. **38**(6), 597–605 (2008).

**9. **D. J. Evers, B. H. W. Hendriks, G. W. Lucassen, and T. J. M. Ruers, “Optical spectroscopy: current advances and future applications in cancer diagnostics and therapy,” Future Oncol. **8**(3), 307–320 (2012).

**10. **M. A. Calin, S. V. Parasca, R. Savastru, M. R. Calin, and S. Dontu, “Optical techniques for the noninvasive diagnosis of skin cancer,” J. Cancer Res. Clin. Oncol. **139**(7), 1083–1104 (2013).

**11. **L. Cunha, I. Horvath, S. Ferreira, J. Lemos, P. Costa, D. Vieira, D. S. Veres, K. Szigeti, T. Summavielle, D. Mathe, and L. F. Metello, “Preclinical imaging: an essential ally in modern biosciences,” Mol. Diagn. Ther. **18**(2), 153–173 (2014).

**12. **Y.-F. Dong, Q.-P. Lu, H.-Q. Ding, and H.-Z. Gao, “Study on the best detector-distance of noninvasive biochemical examination by Monte Carlo simulation,” Spectrosc. Spect. Anal. **34**(4), 942–946 (2014).

**13. **S. Rogalla and C. H. Contag, “Early cancer detection at the epithelial surface,” Cancer J. **21**(3), 179–187 (2015).

**14. **M. J. C. V. Gemert, S. L. Jacques, H. J. C. M. Sterenborg, and W. M. Star, “Skin optics,” IEEE Trans. Biomed. Eng. **36**(12), 1146–1154 (1989).

**15. **R. H. Wilson and M.-A. Mycek, “Models of light propagation in human tissue applied to cancer diagnostics,” Technol. Cancer Res. Treat. **10**(2), 121–134 (2011).

**16. **E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. **13**(6), 060504 (2008).

**17. **E. Alerstam, W. C. Y. Lo, T. D. Han, J. Rose, S. Andersson-Engels, and L. Lilge, “Next-generation acceleration and code optimization for light transport in turbid media using GPUs,” Biomed. Opt. Express **1**(2), 658–675 (2010).

**18. **M. N. Kholodtsova, V. B. Loschenov, C. Daul, and W. Blondel, “Pre-processing method to improve optical parameters estimation in Monte Carlo-based inverse problem solving,” in *Biophotonics: Photonic Solutions for Better Health Care IV*, J. Popp, V. V. Tuchin, D. L. Matthews, F. S. Pavone, and P. Garside, eds., Proc. SPIE9129, 91291Q (2014)

**19. **I. Fredriksson, L. Marcus, and T. Stromberg, “Inverse Monte Carlo method in a multilayered tissue model for diffuse reflectance spectroscopy,” J. Biomed. Opt. **17**(4), 047004 (2012).

**20. **M. Sharma, R. Hennessy, M. K. Markey, and J. W. Tunnell, “Verification of a two-layer inverse Monte Carlo absorption model using multiple source-detector separation diffuse reflectance spectroscopy,” Biomed. Opt. Express **5**(1), 40–53 (2013).

**21. **Y. Zhang, J. Zhu, W. Cui, W. Nie, J. Lie, and Z. Xu, “Threshold thickness for applying diffusion equation in thin tissue optical imaging,” Opt. Commun. **325**, 95–99 (2014).

**22. **Q. Liu and N. Ramanujam, “Sequential estimation of optical properties of a two-layered epithelial tissue model from depth-resolved ultraviolet-visible diffuse reflectance spectra,” Appl. Opt. **45**(19), 4776–4790 (2006).

**23. **C. Liu, N. Rajaram, K. Vishwanath, T. Jiang, G. M. Palmer, and N. Ramanujam, “Experimental validation of an inverse fluorescence Monte Carlo model to extract concentrations of metabolically relevant fluorophores from turbid phantoms and a murine tumor model,” J. Biomed. Opt. **17**(7), 1–15 (2012).

**24. **Q. Wang, D. Le, J. Ramella-Roman, and J. Pfefer, “Broadband ultraviolet-visible optical property measurement in layered turbid media,” Biomed. Opt. Express **3**(6), 1226–1240 (2012).

**25. **M. Juger, F. Florian, and K. Alwin, “Application of multiple artificial neural networks for the determination of the optical properties of turbid media,” J. Biomed. Opt. **18**(5), 1–9 (2013).

**26. **E. Pery, “Biomodal spectroscopy in elastical diffusion and spatially resolved autofluorescence: instrumentation, modeling of light-tissue interactions and application to biological tissue characterization *ex vivo* and *in vivo* for cancer detection (Spectroscopie bimodale en diffusion elastique et autofluorescence resolue spatialement: instrumentation, modelisation des interactions lumiere-tissus et application a la caracterisation de tissus biologiques *ex vivo* et *in vivo* pour la detection de cancer),” Ph.D. thesis, Ecole doctorale IAEM Lorraine, Institut National Polytechnique de Lorraine, France (PhD thesis in French) (2007).

**27. **E. Pery, W. C. P. M. Blondel, C. Thomas, and F. Guillemin, “Monte Carlo modeling of multilayer phantoms with multiple fluorophores: simulation algorithm and experimental validation,” J. Biomed. Opt. **14**(2), 024048 (2009).

**28. **M. N. Kholodtsova, I. S. Samsonova, W. C. P. M. Blondel, and V. B. Loschenov, “Metal nanoparticles of different shapes influence on optical properties of multilayered biological tissues,” in *Medical Laser Applications and Laser-Tissue Interactions VII*, L.D. Lilge and R. Sroka, eds., Proc. SPIE9542, 954205 (2015).

**29. **M. N. Kholodtsova, V. B. Loschenov, C. Daul, and W. C.P.M. Blondel, “Particle swarm optimisation algorithm for Monte Carlo-based inverse problem solving,” in *Proceedings of IEEE conference on Laser Optics*, (Institute of Electrical and Electronics Engineers, 2014), p. 1.

**30. **M. N. Kholodtsova, P. V. Grachev, T. A. Savelieva, N. A. Kalyagina, W. Blondel, and V. B. Loschenov, “Scattered and fluorescent photon track reconstruction in a biological tissue,” Int. J. Photoenergy **2014**1–7 (2014).

**31. **S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch, “A Monte Carlo Model of Light Propagation in Tissue,” in *Dosimetry of Laser Radiation in Medicine and Biology*, Proc. SPIE, 102–111 (1989).

**32. **L. Wang, S. L. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. **47**(2), 131–146 (1995).

**33. **C. Zhu, Q. Liu, and N. Ramanujam, “Effect of fiber optic probe geometry on depth-resolved fluorescence measurements from epithelial tissues: a Monte Carlo simulation,” J. Biomed. Opt. **8**(2), 237–247 (2003).

**34. **S. Feng, F.-A. Zeng, and B. Chance, “Photon migration in the presence of a single defect: a perturbation analysis,” Appl. Opt. **34**(19), 3826–3837 (1995).

**35. **A. N. Bashkatov, E. A. Genina, V. I. Kochubey, and V. V. Tuchin, “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D: Appl. Phys. **38**(15), 2543 (2005).

**36. **I. Boussaid, J. Lepagnot, and P. Siarry, “A survey on optimization metaheuristic,” Inform. Sciences **237**, 81–117 (2013).

**37. **A. El Dor, M. Clerc, and P. Siarry, “A multi-swarm PSO using charged particles in a partitioned search space for continuous optimization,” Comput. Optim. Appl. , **53**, 271–295 (2012).

**38. **J. Kennedy and R. Eberhart, “Particle swarm optimization,” in *Proceedings of IEEE International Conference on Neural Networks* (IEEE, 1995), pp. 1942–1948.

**39. **Y. Shi and R. Eberhart, “A modified particle swarm optimizer,” in *Proceedings of IEEE International Conference on Evolutionary Computation* (IEEE, 1998), pp. 69–73.

**40. **J. Kennedy and R. Mendes, “Population structure and particle swarm performance,” in Proceedings of the IEEE Congress on Evolutionary Computation (IEEE, 2001), 2, 1671–1676.

**41. **D. Tian, “A Review of Convergence Analysis of Particle Swarm Optimization,” Int. J. Grid. Distr. Comput. **6**(6), 117–128 (2013).

**42. **Q. Wang, D. Le, J. Ramella-Roman, and J. Pfefer, “Broadband ultraviolet-visible optical property measurement in layered turbid media,” Biomed. Opt. Express **3**(6), 1226–1240 (2012).

**43. **B. I. Schmitt, “Convergence Analysis for Particle Swarm Optimisation,” PhD thesis, FAU University, Erlangen (2015).

**44. **A. N. Bashkatov, A N. E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues–a review,” J. Innov. Opt. Health Sci. **4**(1), 9–38 (2011).

**45. **S. L. Jacques, “Optical properties of biological tissues: A review,” Phys. Med. Biol. **58**(11), R37–R61 (2013).

**46. **G. Mie, “Contributions to the optics of turbid media, particularly of colloidal metal solutions (Beitrage zur Optik truber Medien, speziell kolloidaler Metallosungen),” Annalen der Physik **330**(3), 377–445 (1908) (In German).

**47. **A. Ishimaru, “Wave propagation and scattering in random media and rough surfaces,” in Proceedings of IEEE Photovoltaic Specialists Conference, (IEEE, 1991), pp. 1359–1366.

**48. **A. N. Bashkatov, E. A. Genina, V. I. Kochubey, M. M. Stolnitz, T. A. Bashkatova, O. V. Novikova, A. Y. Peshkova, and V. V. Tuchin, “Optical properties of melanin in the skin and skin-like phantoms,” in *Controlling Tissue Optical Properties: Applications in Clinical Study*, V.V. Tuchin, ed., Proc. SPIE, 4162 (2000).

**49. **D. K. Sardar, M. L. Mayo, and R. D. Glickman, “Optical characterization of melanin,” J. Biomed. Opt. **6**(4), 404–411 (2001).