Open Access
Add to BibSonomy
Abstract
The estimation of skin optical properties by means of inverse problem solving from spatially resolved diffuse reflectance (SR-DR) spectra is one way to exploit the acquired clinical signals. This method requires the comparison between the experimental spectra collected with a medical device, and spectra generated by the photons transport numerical simulations. This comparison is usually limited to spectral shape due to the absence of intensity standardization of the experimental DR spectra. This study proposes to theoretically (using photometric calculation) and experimentally (from experimental spectra acquired on optical phantom) establish a corrective factor to obtain common intensity unit for experimental and simulated signals.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
In the context of non-destructive optical characterization of biological tissues, Diffuse Reflectance Spectroscopy (DRS) has been widely used for several years [1–4] for cancer detection. It consists in back-reflected intensity spectra measurement. Such spectra carry information from light-tissue interactions acquired after a broadband light illumination, from which the optical properties (related to absorption and scattering) as well as structural information about the probed medium can be extracted and analyzed. Wavelength-dependent absorption and scattering events stochastically determine the trajectories of photons in the medium, so that some of them may be re-emitted toward the surface and potentially captured by a spectrometer. The DRS non-invasive nature, its sensitivity to metabolic and structural changes in tissues and its clinical applicability are very well suited to in vivo diagnostics. To more finely analyze the medium under the probe, several Detection Fibers (DF) located at different distances from the Source Fiber (SF) can be jointly used. Such fibered systems are characterized as Spatially Resolved DR (SR-DR) spectroscopy devices and allows for depth-resolved acquisitions since photons detected at a short/large Source-Detector (SD) distance result from shallow/deep trajectories into the biological tissues under evaluation [5]. The interest of the spatial resolution is all the greater for layered organ such as skin (made of epidermis, dermis and their sub-divisions). That is why the literature on this subject points numerous applications of SR-DR spectroscopy to in vivo characterization of cutaneous pathological states including melanoma [6–8] or carcinoma [9,10], and more generally to the characterization of skin layers optical properties (OP) [11,12].
A way to get these OP is to solve the inverse problem that consists in establishing a simulation that considers both the device features (geometrical and optical) and a skin multi-layer model, in which optical (e.g. absorption and scattering coefficients, anisotropy factor) and geometrical (e.g. layer thicknesses) features can be estimated through an optimization process to fit clinically acquired “target” SR-DR spectra [13–16]. This fitting procedure then requires a comparison between clinical (experimental) and numerical (simulated) spectra, that is often quantified by a least squares error cost function.
This comparison is not trivial because obtaining DR spectra intensities in absolute values, defined as the fraction of collected over emitted light and expressed in photons ratio, still presents difficulties for experimental signals. Indeed, the raw spectra in photon counts provided by the spectrometer, is dependent on the broadband white light source spectral features, on the spectral responses of the probed medium and on the components of the optoelectronic acquisition chain i.e. fibers, grating, and photo-detector. To compensate for the light source spectral shape and for the acquisition chain spectral response, which can also vary over time, it is an established practice to carry out a calibration measurement on a spectrally flat reflectance standard, made of spectralon [17,18]. For the same illumination condition, by dividing raw spectra acquired on skin by raw spectra acquired on spectralon, the resulting spectral shape depends only on the probed medium spectral response (i.e. spectral shape correction). However, the intensity unit of this spectrum is difficult to interpret and not only depends on the optical probe geometrical features (SD fiber distances, fibers core diameter and numerical apertures), but also on the calibration measurement set up: distance between the optical fiber probe tip and the spectralon surface.
This lack of standardization (i.e. absolute intensity correction) in experimental DR spectra limits any comparison with simulated ones and it only become possible after normalizing both spectra to bring them back to a common scale (imposing that both maximums correspond to 1). More generally, any comparison of experimental data requires normalization, especially when confronting data to those provided by other research teams raising the issue of inter-devices comparison. Therefore, the definition of a corrective factor allowing to standardize any experimental SR-DR spectrum to absolute magnitude is at stake.
To the best of our knowledge, contributions that explain how to absolute compare simulated and experimental spectra are rare and only deal with experimental values of the calibration factor without dealing with calculation methods. Thueler et al. [17] developed a two-step calibration procedure using a (i) spectrally flat reflectance standard and (ii) a solid turbid siloxane phantom of known optical properties to get absolute reflectance spectra of stomach tissues. The optical properties were estimated at 675 $nm$ by combining frequency-domain and spatially resolved measurements. Those estimates were used as parameters in Monte Carlo (MC) simulations which allow to get the absolute diffuse reflectance, and thus the calibration factor to finally obtain spectra in the same unit as the simulated ones. A similar approach was also considered by Naglic et al. [12]. The use of turbid phantoms [19] to standardize the intensity unit of experimental spectra has the advantage of being more representative of biological media but requires the adaptations of MC photon transport algorithm to geometrically match the acquisitions of an existing device. Other experimental methods, presented by Zhang et al. [20] have been investigated to get rid of the use of MC simulations, but are intended for single fiber reflectance spectroscopy calibration and are therefore incompatible with spatially-resolved configurations i.e. several SD separations. The current study will present a unified approach of experimental calibration factor calculation for SR-DR spectroscopy using a liquid optical phantom of known optical properties and an adapted photon transport simulation. This experimental factor will be compared with theoretical calculations of the same factor using photometry, that only requires the knowledge of some optical parameters of the probe fibers and the geometry of the reflectance measurement configuration on the spectrally flat reflectance standard.
2. Materials and methods
2.1 Geometrical and spectral description of the SR-DR spectroscopy device and uncertainties
This section introduces the technical characteristics of the SpectroLive device [21] used for SR-DR spectra acquisitions on optical phantoms and on human skin in vivo. The multiple fiber optic (MFO) probe is made of a central source fiber (SF) featuring a radius $r_{SF}=300~\mu m$, numerical aperture $NA_{SF}=0.22$ and refractive index $n_{DF}=1.47$ delivering broadband $[365-765]~nm$ white light (WL) illumination to the tissue in contact. 24 detection fibers (DF) each featuring a radius $r_{DF}=100~\mu m$ and $NA_{DF}=0.22$, are equally distributed on 4 concentric circles around the central SF corresponding to the 4 center-to-center distances $D_n$ ($n\in \{1,\ldots,4\}$) ranging approximately from 400 to 1000 $\mu m$ by 200 $\mu m$ step. Real center-to-center $D_n$ distances measured directly on the probe tip, resulting from the self-organization of the fibers during the manufacturing process of the bundle, are provided in Table 1. The photons reaching one of the 6 collection fibers of the $D_n$ distance are gathered to create the diffuse reflectance signal $DR_{D_n}(\lambda )$. A schematic representation of the acquisition geometry is shown in Fig. 1, while the numerical values of the different parameters used are summarized in the Table 1.
Fig. 1. Geometrical disposition of the Source Fiber (SF) and Detection Fibers (DF) at the tip of the spectroscopic probe used for spatially resolved diffuse reflectance measurements.
Table 1. Values of the geometrical and spectral source-detector features of the SR-DRS Spectrolive device (also used for numerical simulations).a
In practice, the DF of the SR-DR probe are self-arranged around the SF and held together with an epoxy resin subsequently polished to ensure a clean glass-skin interface between the end of the fibers and the probed medium. This process leads to real geometrical features with inherent variabilities in (i) the center-to-center SD distance $D_n$ (characterized by $\sigma _{Dn}$) and (ii) their vertical orientations (characterized by $\alpha$ angle) relative to the normal of the interface between the SF and the medium under the probe tip. Both were estimated as $\sigma _{Dn}=\pm 30~\mu m$ and $\alpha =\pm {2}^{\circ}$, respectively. In order to consider the impact of these variations in the calculation of the calibration factor, three geometric configurations have been considered for the MFO probe, as represented in Fig. 2. The first configuration, or “overestimated” case, considers for both SD distance and optical axis orientation the maximal values of variability in favor of an increase of the detected signal, i.e. a shorter effective SD distance $D_{n}^{over}=D_n-\sigma _{Dn}$ and converging optical axes defined by $\alpha ^{over}=-\alpha$. The “perfect” case matches the ideal probe feature $D_{n}^{perfect}=D_n$, without any angle between SF and DF axes i.e. $\alpha ^{perfect}= {0}^{\circ}$. Finally, the “underestimated” case considers the maximal uncertainties in favor of a decrease of the detected signal, i.e. a longer effective SD distance $D_{n}^{under}=D_n+\sigma _{Dn}$ and diverging optical axes, $\alpha ^{under}=\alpha$.
Fig. 2. Representation of the geometric uncertainties of the probe tip for $D_n$ SD distance. Left, "overestimated" case: DF and SF are closer ($D_n-\sigma _{Dn}$) and their optical axes converge ($\alpha = -{2}^{\circ}$). Middle, "perfect" case: ideal probe with perfect $D_n$ SD distance and parallel optical axes. Right, "underestimated" case: DF and SF are furthest away ($D_n+\sigma _{Dn}$) and their optical axes diverge ($\alpha = +{2}^{\circ}$).
2.2 Current calibration and measurement procedure
For the reasons mentioned above, it is common practice to firstly calibrate the device using reflectance standard. The process to obtain a $DR_{Dn}(\lambda )$ spectra consists of the three following steps:
- 1. Measurement on a reflectance standard: The MFO probe is vertically placed in a light-isolated box at a distance $z$ ($z=10000~\mu m$ in our case) from a spectralon reflectance standard (SRS-99, Labsphere). The latter provides a UV-Vis-NIR spectrally flat reflectance of a fraction $R$ of incident WL in the entire half-space on the source side (Lambertian-like behavior). A part of the reflected photons is then captured by the detection fibers located at all $D_n$ distances and collected by spectrometers to generate the reflectance “raw” spectra $S_{Dn}^{stand}(\lambda,z,R)$ whose unit are in photon counts.
- 2. Measurement on skin: The MFO probe is applied on skin surface with gentle contact, and same incident WL illumination is provided. Raw spectra $S_{Dn}^{skin}(\lambda )$ are then measured (also expressed in photon counts).
- 3. DR spectra calculation: From the two previous measurements, the diffuse reflectance spectra $DR_{D_n}(\lambda,z,R)$ are calculated as:$$DR_{D_n}(\lambda,z,R)=R\frac{S_{Dn}^{skin}(\lambda)-S_{Dn}^{dark}(\lambda)}{S_{Dn}^{stand}(\lambda,z,R)-S_{Dn}^{dark}(\lambda)},$$where $S_{Dn}^{dark}(\lambda )$ is the dark photon counts spectra for both acquisitions on spectralon and skin.
In Eq. (1), the intensity unit of $DR(\lambda )$ at SD distance $D_n$ is the photon count ratio between the light collected by the probe in contact with skin and the light collected from the Lambertian standard of reflectivity $R$ by the probe located at a given distance $z$ between probe tip and standard surface. The following text will refer to this signal ratio as “experimental unit”. This dimensionless non-elementary unit is difficult to interpret. This problem may be solved by determining a way to transform these experimental spectra into spectra expressed in absolute unit, i.e. the ratio between the light received by the DF and sent by the SF. In practice, it consists in introducing a standardization factor $c_{D_n}(z,R)$ such that:
The latter factor $c_{D_n}(z,R)$ thus corresponds to the ratio between collected and emitted light during calibration measurement i.e. the normal illumination of a Lambertian standard of reflectance $R$ at a distance $z$. The next sections present the two approaches proposed to theoretically calculate this value: the first one based on radiometric calculations and the second one using MC simulation based (Ray-tracing) calculations.
2.3 Geometry of the calibration measurement
This section describes geometrical configuration of illumination and reflected light distributions during the calibration measurement with a MFO probe. A schematic representation of the successive steps is given in Fig. 3. The incident flux $F_{inc}$ emitted by the SF propagates in a light cone according to the $NA_{SF}$ divergence and illuminates a given $A_{SRS-99}(z)$ area on the standard surface. This area, increasing with $z$, is then considered as a secondary Lambertian light source emitting the reflected flux $F_{ref}=RF_{inc}$ in all the upper half-space. Only a part of $F_{ref}$ is collected by detection fibers located at $D_n$ from SF satisfying the angular detection conditions of $NA_{SF}$. This flux is called $F_{det}(z,R)$ and the light ratio, i.e. $c_{D_n}(z,R)$ is then defined as:
Fig. 3. Schematic (not scaled) representation of the calibration geometry. Left, from top to bottom: illumination geometry with geometry features of the probe tip, side view of the incident light cone flux $F_{inc}$ illuminating area $A_{SRS-99}(z)$ on the reflectance standard, top view. Right: Reflectance geometry with side view of the reflected light distribution showing the half-space reflected flux $F_{ref}$ and the detected part of it $F_{det}(z,R)$ (here by DF4).
2.4 Photometric calculation of the correction factor
In this section, the theoretical calculation of the $c_{D_n}(z,R)$ factor is presented based on radiometry, developed hereafter in 3 steps. For this purpose, a schematic representation of the reflected light geometry is provided in Fig. 5.
- 1. Given that only the light ratio is of interest (see Eq. (3)), the incident flux emitted by the SF is arbitrarily considered to be unitary, i.e. $F_{inc}=1~W$. $F_{inc}$ is in-homogeneously spread according to the angular distribution at the SF output over a disk area $A_{SRS-99}(z)$ of radius $r(z)=r_{SF}+z\cdot \tan (\arcsin (NA_{SF}))$ on the surface of the reflectance standard. An image of this real flux distribution pattern exited from the SF and reaching the reflectance standard surface is shown in the left part of Fig. 4. The latter was acquired by imaging with a camera the light at the SF output through a thin and neutral (i.e. spectrally flat) scattering sheet. The flux quantified by surface unit, i.e. the irradiance, around point $P(r,\theta,0)_{\mathcal {R}}$ in the cylindrical coordinates system associated to the mark $\mathcal {R}(O,\vec{u_x},\vec{u_y},\vec{u_z})$ with origin $O$ in the center of $A_{SRS-99}(z)$ is defined as the ratio between elementary flux and the disk area element $dA=rdrd\theta$. Thus, the irradiance in each point $P$ on the reflectance standard can be written:
- 2. Luminance $L(P,\vec{u},z)$ defines the quantity of the light emitted from a point $P$ in the direction $\vec{u}$. In the case of a Lambertian light source like the one re-emitted by the SRS-99 spectralon surface (see schematic representation in the right part of Fig. 4), the luminance is (i) equal in all the directions i.e. $L(P,\vec{u},z)=L(P,z)$ and (ii) distributed in the entire half-space. Combining both conditions, the luminance is related to the aforementioned irradiance $I(P,z)$ and standard reflectance $R$ by the formula:
- 3. Then, an infinitesimal integral calculation is performed to determine the fraction of the flux re-emitted by the standard which can be captured by detection fiber. Analogously to the elementary surface $dA$, we introduce the surface element $dA'=r'dr'd\theta '$ (see Fig. 5) located on the collecting surface $S_{DF}$ of the DF, i.e. inside a disk centered at $O'(0,D_n,z)_{\mathcal {R}}$ (in Cartesian coordinates) around the point $P'(r',\theta ',0)_{\mathcal {R}'}$ defined by its cylindrical coordinates in the secondary system $\mathcal {R}'(O',\vec{u_x},\vec{u_y},\vec{u_z}))$. From the radiometry laws, the infinitesimal flux $dF_{det}^{rad}$ coming from the surface $dA$ and captured by the surface $dA'$ (see Fig. 5) is given by the formula: where $\phi$ is the angle between the normal to the surface and the direction $\vec{u}_{PP'}$ provided by the line connecting the source point $P$ and the detection point $P'$, and $d\Omega$ is the solid angle defined as the ratio between the projected surface of $dA'$ on $\vec{u}_{PP'}$ direction and the square of distance $d$ between $P$ and $P'$ i.e. $d\Omega =dA'\cos (\phi )/d^2$. The angular detection condition related to numerical aperture and the Fresnel transmission are not considered in Eq. (6). To address this, the Fresnel transmission coefficient in intensity $R_F$ has to be introduced as:$$R_F(\phi)=\frac{n_{DF}\cos(\phi')}{n_{air}\cos(\phi)}t^2 \quad \textrm{where} \quad t=\frac{2n_{air}\cos(\phi)}{n_{air}\cos(\phi)+n_{DF}\cos(\phi')},$$with $\phi '$ linked to $\phi$ according to the refraction Snell–Descartes law ($n_{air}\sin (\phi )=n_{DF}\sin (\phi ')$). The elementary detected flux $dF_{det}$ that really enters the SF is finally linked to $dF_{det}^{rad}$ by the formula:
This case disjuncture enables the collection of photons that verify the angular acceptability of the fiber. Considering all the surface elements of the source and the detection surfaces, the detected flux $F_{det}$ is obtained after double integrals on $A_{SRS-99}(z)$ and $S_{DF}$ surfaces:
Fig. 4. Distribution of the flux on the reflectance standard at the fiber output (left) and schematic representation of the angular emission of the photon flux for a Lambertian source (right).
Fig. 5. Geometrical representation of the reflected light configuration with the parameters involved in radiometry calculation of $c_{D_n}(z,R)$.
In practice, the numerical calculation of the double integrals in Eq. (4) was performed on Matlab with a discretization of the surfaces $A_{SRS-99}(z)$ and $S_{DF}$ into 900 elements corresponding to 30 angular and 30 radial divisions.
2.5 Numerical simulation-based calculation of the correction factor
In order to cross-validate the radiometric calculations undertaken so far with numerical simulations of photon propagation, a Monte Carlo (MC)-based algorithm of Ray-Tracing (described in Algorithm 1 and available in Appendix Section 5) was designed to compute $c_{D_n}(z,R)$. Its comparison with the photometric calculation is done in section 3.1.
2.6 Experimental validation
2.6.1 General presentation of the experimental procedure
This section describes the experimental approach developed to calculate the correction factor $c_{D_n}^{exp}(z,R)$ from SR-DR measurements on optical phantoms combined to MC-based photon transport simulation. The experiment consisted in acquiring SR-DR spectra of a mono-layer liquid optical phantom with known OP (see section 2.6.2). In Eq. (2) the latter $c_{D_n}^{exp}(z,R)$ was defined as the ratio between the spectrum of DR (expressed in photon ratio) and those calculated with the current method described in section 2.2 for the same probed medium, i.e. $DR_{D_n}(\lambda,z,R)$. The MFO probe was fixed on a micrometer translation stage allowing for the varying of the distance $z$ between the probe tip and the SRS-99 reflectance standard surface. In order to obtain SR-DR spectra in photon ratio unit, a CudaMCML light transport simulation [22] was adapted to match the acquisition geometry of the experimental device. From the MC simulations, the output unit of simulated spectra $DR_{Dn}^{sim}(\lambda )$ is intrinsically the photon ratio between detected and emitted light (or photon weights). The geometrical and optical parameters of the optical phantom required for this simulation were its thickness, the absorption $\mu _a(\lambda )$ and scattering $\mu _s(\lambda )$ coefficients, as well as the anisotropy factor $g(\lambda )$. For MC simulations, a large number ($10^9$) of photons were launched, thus ensuring a statistical variance of less than 0.5% on the SR-DR simulated spectra, whatever the SD distance and whatever the wavelength. The simulation finally provides the $DR_{Dn}^{sim}(\lambda )$ spectra, allowing to obtain the following experimental factor:
2.6.2 Optical phantom description
A 500 $mL$ liquid optical phantom made of Methylen Blue (MB, M9140, SigmaAldrich, Missouri, USA) at a concentration of 8 $\mu M$ and 4% (vol/vol) Intralipid20% (68890-65-3, SigmaAldrich, Missouri, USA) was poured into a 50 $mm$ side black square silicone mold. The associated scattering OP ($\mu _s(\lambda )$ and $g(\lambda )$ linked to the intralipids concentration) were calculated according to Aernouts et al. [23] paper, while the absorption coefficient ($\mu _a(\lambda )$ imposed by the MB concentration) was calculated from OMLC database [24]. These optical properties are shown in Fig. 12 in appendix 6.
2.6.3 SpectroLive acquisition
For sake of repeatability, ten $S_{Dn}^{skin}(\lambda )$ spectra were acquired on the optical phantom and averaged by repetitively immersing the SpectroLive device optical probe into the liquid phantom (see right picture in Fig. 6). Concerning the signals acquired on the reflectance standard $S_{Dn}^{stand}(\lambda,z,R)$, 76 raw spectra were collected for $z$ distances in a range $[0-15000]~\mu m$ with a 200 $\mu m\;\Delta z$ step (see left picture in Fig. 6). Those acquisitions allow then to finally get the $DR_{D_n}(\lambda,z,R)$ spectra, presented and discussed in the following results and discussion section.
Fig. 6. The two measurements $S_{Dn}^{stand}(\lambda,z,R)$ (left) and $S_{Dn}^{skin}(\lambda )$ (right) allowing the $DR_{D_n}(\lambda,z,R)$ experimental calculation. The distance $z$ between the standard surface and the probe tip was varied using micrometer translation stage.
3. Results and discussion
3.1 Comparison between results from photometric theory and Monte Carlo-based simulations
The curves of the correction factor $c_{D_n}(z,R=0.99)$ obtained using radiometric calculation (see section 2.4) and MC ray-tracing approaches (see section 2.5 and appendix 5) are plotted in left Fig. 7 (continuous line for radiometry and dashed-line for MC) for the 4 $D_n$ distances of the MFO probe and for the distance $z$ in the range $[0-2.10^4]~\mu m$ between the probe tip and the reflectance standard surface. $c_{D_n}(z,R=0.99)$ curves feature with a zero values for $z<z_{Dn}^{min}$ corresponding to minimal distances allowing for the optical flux to enter the DF, before reaching max value $c_{D_n}^{max}$ at “optimal” distances $z_{Dn}^{max}$. This position shifts toward greater $z$ when $D_n$ increases. The calculation of the the standardisation factor using the ray-tracing approach is based on the geometrical approximation of considering annular detected surfaces rather than real point DF (see Algorithm 1). The latter approximation allowed to have very low MC statistical noise while saving computation time, by considering that the light distribution was uniform on the whole annular surface (see Fig. 1). Error $\epsilon ^c_{Dn}(z)$, shown in right Fig. 7, corresponding to the difference between $c_{D_n}(z,R=0.99)$ obtained from radiometric and MC-based calculations, was introduced to quantify the impact of this assumption on the factor calculating. As expected, its impact is limited (at maximum 2.17% in relative error) and relatively higher for small SD distances, with $\epsilon _{D1}\approx 1\times 10^{-4}$ at $z_{Dn}^{max}$ for $D_1$ whereas $\epsilon _{D4}\approx 1\times 10^{-5}$ (or 1%) for $D_4$. Despite these errors, the two approaches validate each other. Consequently, $c_{Dn}(z,R)$ results obtained by photometry calculation will be considered as reference values.
Fig. 7. Left: Photometric (line, see section 2.4) and Monte Carlo (dashed-line, see section 2.5) calculations of the $c_{D_n}(z,R)$ factor. Right: Error $\epsilon _{Dn}$ between radiometric and MC-based calculations. Signal to noise ratio in $\epsilon ^c_{Dn}(z)$ decreases with $D_n$, this is explained by the stochastic noise in MC calculation.
3.2 Impact of the calibration measurement on DR spectra
Following the experimental protocol described in section 2.6.3, the DR spectra obtained from the SRS-99 standard and expressed in the experimental unit by applying the formula given in Eq. (1) are shown in Fig. 8, as a function of $D_n$ SD distances and for different distances $z$ between the reflectance standard surface and the MFO probe tip. Spectra acquired for $z<z_{D4}^{min}=3800~\mu m$ are not represented, because of infinite values caused by the very weak, i.e. null $S_{Dn}^{stand}(\lambda,z,R)$ signals. Indeed, when $z$ is low, the angular conditions for light entrance into the DF are not met (see Eq. (8)). Weak levels of signals appear also at slightly larger vertical distances $z$ (i.e. $z<6000~\mu m$) especially for $D_4$, which causes noisy DR experimental spectra (see $DR_{D4}(\lambda )$ for $z=3800~\mu m$ in Fig. 8). For large distances (i.e. $z>10000~\mu m$), the low intensity level of spectra is also responsible for noisy $DR_{Dn}(\lambda )$ curves. Thus, to obtain spectra with acceptable noise for all the $D_n$ distances simultaneously, an optimal $z$ calibration distance can be deduced from these obseravtion, within the range $[6000-10000]~\mu m$.
Fig. 8. Diffuse reflectance spectra $DR_{D_n}(\lambda,z,R=0.99)$ in their current experimental unit obtained from SRS-99 standard for 4 SD distances $D_{1,2,3,4} = {400,600,800,1000}~\mu m$ and for 12 standard-to-probe calibration distances $z = {3800,\ldots,14800}~\mu m$.
3.3 Experimental calculation of the standardization factor
Following the protocol described in section 2.6, the final step to perform is the term-to-term quotient of the experimental spectra $DR_{D_n}(\lambda,z,R)$ (the non-normalized ones shown in Fig. 8 as a function of $z$) to the simulated ones $DR_{Dn}^{sim}(\lambda )$ (also non-normalized) using Eq. (10) to get $c_{Dn}^{exp}(z,R)$. A box-plot approach (mean $\pm$ standard deviation) was chosen to consider the wavelength dependency during the calculation of this experimental factor. The resulting curves are shown in Fig. 9 for the 4 $D_n$ SD distances, with the theoretical photometric calculation (cf. section 2.4) in continuous line. In order to quantify the difference between the calculation resulting from the theoretical model (for the “perfect” case, cf. Fig. 2) and the experimental reality, the error $\epsilon ^c_{Dn}(z)=c_{Dn}^{exp}(z,R)-c_{Dn}(z,R)$ was introduced, and plotted in a top insert for each SD distances.
Fig. 9. Comparison between experimental and theoretical values of $c_{D_n}(z,R=0.99)$ for the 4 SD distances. The box-plot experimental representation synthesizes wavelengths dependency for the quotient calculation. Continuous line is for the “perfect” case, whereas dashed-lines are for “underestimated” and “overestimated” cases (cf. Fig. 2) which gives the spread of this factor by considering the possible geometric variabilities of the real optical probe. Top insert provides the error $\epsilon ^c_{Dn}(z)$ between the mean of the experimental value and the theoretical “perfect” case for the factor calculation.
A good correspondence was obtained between experimental and theory “perfect” case. We can notice that maximal error are caused by a delay (i.e. z-shift) of the experimental curves around $z=z_{Dn}^{max}$, of value $5\times 10^{-4}$ ($\approx$ 10%) for $D_1$, $3.1\times 10^{-4}$ ($\approx$ 12.5%) for $D_2$, $3\times 10^{-4}$ ($\approx$ 20%) and $2.6\times 10^{-4}$ ($\approx$ 25%) for $D_4$. The error decreases with increasing of the distance $z$ between the probe tip and the reflectance standard. By summing the absolute relative error of the 4 $D_n$ distances, we identified the distance $z=9200~\mu m$ which minimizes the differences between the radiometric model (“perfect” case) and the mean experimental value, all distances considered. A first source of this error $\epsilon ^c_{Dn}(z)$ is related to the wavelength independence assumption. Indeed, in theoretical calculations, none of the parameters $NA_{DF/SF}$ (numerical aperture of DF and SF), $n_{DF}$ (DF refractive index) and $R$ (reflectance of the standard) considered any chromatic behavior. A second explanation is the geometrical variabilities featuring the real probe tip, that are the effective SD distances and the optical axe orientation (as already described in section 2.1). To quantify the impact of the latter varaitions on the theoretical value of $c_{D_n}(z,R=0.99)$, two additional calculations corresponding to “overestimated” and “underestimated” probe configurations (see Fig. 2) were performed and plotted in black dashed-line in Fig. 9. It is noticeable that those two possible geometrical errors have a strong impact on the shape (start of growth $z_{Dn}^{min}$, upward slope, peak position $z_{Dn}^{max}$ and downward slope) and on the amplitude $c_{Dn}^{max}$ of the calibration factor curve. In particular, one can notice that the experimental calculation of the standardization factor better corresponds to “underestimated” case (larger effective SD distance and divergent optical axes) for $D_3$ SD distance. The latter are to be related to the corresponding geometrical variability featuring the DF located farther from the SF. The proposed method thus potentially allows to characterize and take into account geometrical features of any MFO probe.
3.4 Correction of the experimental DR spectra using the photometric standardization factor
Knowing the calibration height $z$ and the photometric standardization factor $c_{D_n}(z,R=0.99)$ and applying the formula described in Eq. (2) allows to finally get experimental DR spectra with absolute intensity, i.e. photons ratio unit. The Fig. 10 shows the corrected experimental spectra at the 4 $D_n$ SD distances for $z=8000~\mu m$ acquired in the liquid phantom described in section 2.6.2 as well as the numerical ones resulting from the CudaMCML photon transport simulation for a numerical medium with optical properties provided in appendix 6. The “perfect” case (see Fig. 9) photometric standardization factor $c_{Dn}(z=8000,R=0.99)$ was used to scale the $DR_{D_n}(\lambda,z=8000,R=0.99)$ experimental spectra. The global fitting between experimental corrected and numerical DR spectra can be noticed, in particular in the near-IR spectrum, while differences reaching 15% happen in the near-UV range. Although turbid phantom-based approaches coupled with MC photon transport simulation adapted to specific spectroscopic device may allow to obtain better fitting between experimental and simulated spectra [12,17], those involve a consequent experimental step composed of (i) the realization of a turbid phantom, (ii) the accurate optical characterization of the latter, and especially (iii) the implementation of MC simulation including accurate features of the real experimental device. The theoretical alternative proposed here only requires the knowledge of the probe geometrical properties (SD distances, radius and numerical aperture of the emitting and detecting fibers) to led the calculation described in section 2.4 and finally obtaining a scalar standardization factor to get DR spectra at several SD separations with absolute intensity unit. Although it is not the ambition here, the consideration of geometric uncertainties (see Fig. 2 and Fig. 9) in the theoretical calculation of the standardization factor and in the geometric model used in the MC simulation would further reduce the gaps between simulated and experimental spectra.
Fig. 10. Comparison between simulated and experimental corrected DR spectra for acquisition in the liquid phantom made of intralipids ans methylen blue described in section 2.6.2. The experimental spectra was calculated using Eq. (1) with measurement on a reflectance standard such that $z=8000~\mu m$, before being corrected according to Eq. (2), i.e. multiplied by $c_{Dn}(z=8000,R=0.99)$ “perfect” case (see Fig. 9).
4. Conclusion
The ambition of this contribution was to define a correction factor able to link the unit of experimentally acquired DR spectra to spectra expressed in photons ratio, i.e. the fraction of detected over emitted light. For the first, the current study provides a method to not only compare the spectral shapes of the DR signal (between simulated and experimental spectra, or experimental ones from several devices), but also their amplitude. Obtaining this theoretical factor leads us to quantify the fraction of light collected by a DF during calibration measurement, i.e. the illumination of spectralon standard (flat spectral behavior) at a given $z$ vertical distance. An experimental protocol involving (i) a liquid optical phantom of known optical properties, (ii) a simulation of DR spectra acquisition faithful to the real device, and (iii) finally real acquisitions in this same optical phantom was designed and followed to get an experimental value of such a factor. Taking into account geometric approximations in the theoretical calculation, the similarity of the experimental and theoretical curves make it possible to validate the calculation model of the factor.
5. Appendix: Monte Carlo-based algorithm of Ray-Tracing to compute the correction factor
This appendix refers to the calculation of $c_{D_n}(z,R)$ using a MC-based Ray-Tracing, as described in section 2.5. The pseudo-code algorithm and its explanation are developed hereafter:
Algorithm 1. MC simulation to determine the fraction of light reaching DF.
Some details are provided here to explain the light re-emission direction $\vec{u}$ random drawing to simulate a Lambertian bidirectional reflectance distribution [25]. The random draw which allows to reproduce such a reflection was performed using two pseudo random numbers $\theta$ (azimuthal angle) and $\phi$ (polar angle) obtained from two $[0-1]$ uniform variables $u_{\theta }$ and $u_{\phi }$ such as:
The cosine behavior in $\phi$ distribution corresponding to Lambertian source was validated in right sub-figure of Fig. 11. The consideration of annular surface of collection (see Fig. 1) rather than point circular ones, identical to the real device, was motivated by the need to optimize the number of photons launched and the computing duration time. Indeed, increasing the collection surface allows to collect more photons relative to the surface of a single DF. By taking care to multiply the number of photons detected by a correction factor defined as the ratio between the DF area and the associated annular area, this finally provides a geometry almost equivalent to the real device while reducing the number of photons to launch and thus the computation time while keeping low level of statistical noise.
Fig. 11. Light flux distribution on the reflectance standard used for MC calculation (left), representation of the random drawing of the unit directing vector $\vec{u}$ appearing in Algorithm 1 and Eq. (11) (middle), and verification of the Lambertian cosine behavior in $\phi$ distribution (right).
6. Appendix: Optical properties of the liquid phantom used for the experimental validation
This appendix presents in Fig. 2.6.2 the absorption and scattering optical properties of the liquid optical phantom made of intralipids and methylen blue described in section 2.6.2.
Fig. 12. Optical properties of the liquid phantom described in section 2.6.2: absorption coefficient $\mu _a(\lambda )$ (left) ; scattering coefficient $\mu _s(\lambda )$ (middle) ; anisotropy factor $g(\lambda )$ (right).
Funding
Université de Lorraine (PhD grant R01PJYRX-PHD19-COLA-CRAN); Ligue Contre le Cancer (CPER IT2MP, plateforme IMTI); European Regional Development Fund (CPER IT2MP, plateforme IMTI); Contrat de Plan Etat-Région Grand Est 2015-2020 (CPER IT2MP, plateforme IMTI)); Agence Nationale de la Recherche Spec-LCOCT project grant (ANR-21-CE19-0056).
Acknowledgments
The authors want to thank the high Performance Computing resources were partially provided by the EXPLOR centre hosted by the Université de Lorraine (project 2018AM2IX0877).
Disclosures
The authors declare no conflicts of interest.
Data availability
No data were generated or analyzed in the presented research.
References
1. T. M. Bydlon, R. Nachabé, N. Ramanujam, H. J. C. M. Sterenborg, and B. H. W. Hendriks, “Chromophore based analyses of steady-state diffuse reflectance spectroscopy: current status and perspectives for clinical adoption,” J. Biophotonics 8(1-2), 9–24 (2015). [CrossRef]
2. A. N. Bashkatov, E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: A review,” J. Innovative Opt. Health Sci. 04(01), 9–38 (2011). [CrossRef]
3. N. Kollias, I. Seo, and P. R. Bargo, “Interpreting diffuse reflectance for in vivo skin reactions in terms of chromophores,” J. Biophotonics 3(1-2), 15–24 (2009). [CrossRef]
4. R. H. Wilson, K. P. Nadeau, F. B. Jaworski, B. J. Tromberg, and A. J. Durkin, “Review of short-wave infrared spectroscopy and imaging methods for biological tissue characterization,” J. Biomed. Opt. 20(3), 030901 (2015). [CrossRef]
5. G. Greening, A. Mundo, N. Rajaram, and T. J. Muldoon, “Sampling depth of a diffuse reflectance spectroscopy probe for in-vivo physiological quantification of murine subcutaneous tumor allografts,” J. Biomed. Opt. 23(08), 1 (2018). [CrossRef]
6. A. Garcia-Uribe, L. V. Wang, J. Zou, E. B. Smith, M. Duvic, and V. Prieto, “In-vivo characterization of optical properties of pigmented skin lesions including melanoma using oblique incidence diffuse reflectance spectrometry,” J. Biomed. Opt. 16(2), 020501 (2011). [CrossRef]
7. L. Pires, V. Demidov, I. A. Vitkin, V. S. Bagnato, C. Kurachi, and B. C. Wilson, “Optical clearing of melanoma in vivo : characterization by diffuse reflectance spectroscopy and optical coherence tomography,” J. Biomed. Opt. 21(8), 081210 (2016). [CrossRef]
8. M. Amouroux, W. Blondel, F. Granel-Brocard, F. Marchal, and F. Guillemin, “A preliminary study on skin phantoms to test spatially resolved-Diffuse Reflectance Spectroscopy as a tool to help diagnose cutaneous melanoma: A non-invasive measurement of Breslow index,” Bio-Med. Mater. Eng. 18(4-5), 339–343 (2008). [CrossRef]
9. G. Khairallah, W. Blondel, P. Rakotomanga, C. Soussen, A. Delconte, F. Plenat, F. Marchal, and M. Amouroux, “Impact of standardization of autofluorescence and diffuse reflectance spectra on diagnosis accuracy of optical spectroscopy used for skin carcinomas diagnosis,” in In Saratov Fall Meeting SFM'18 International Symposium”Optics and Biophotonics-V, (2018).
10. M. Amouroux, G. Díaz-Ayil, and W. C. P. M. Blondel, “Classification of ultraviolet irradiated mouse skin histological stages by bimodal spectroscopy: multiple excitation autofluorescence and diffuse reflectance,” J. Biomed. Opt. 14(1), 014011 (2009). [CrossRef]
11. A. Koenig, B. Roig, J. L. Digabel, G. Josse, and J.-M. Dinten, “Accessing deep optical properties of skin using diffuse reflectance spectroscopy,” in Clinical and Biomedical Spectroscopy and Imaging IV, vol. 9537, eds., International Society for Optics and Photonics (SPIE, 2015), pp. 42–49.
12. P. Naglič, F. Pernuš, B. Likar, and M. Bürmen, “Estimation of optical properties by spatially resolved reflectance spectroscopy in the subdiffusive regime,” J. Biomed. Opt. 21(9), 095003 (2016). [CrossRef]
13. I. Fredriksson, M. Larsson, and T. Strömberg, “Inverse Monte Carlo method in a multilayered tissue model for diffuse reflectance spectroscopy,” J. Biomed. Opt. 17(4), 047004 (2012). [CrossRef]
14. 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 (2014). [CrossRef]
15. G. M. Palmer and N. Ramanujam, “Monte carlo-based inverse model for calculating tissue optical properties. part i: Theory and validation on synthetic phantoms,” Appl. Opt. 45(5), 1062–1071 (2006). [CrossRef]
16. H. Cen, R. Lu, and K. Dolan, “Optimization of inverse algorithm for estimating the optical properties of biological materials using spatially-resolved diffuse reflectance,” Inverse Probl. Sci. Eng. 18(6), 853–872 (2010). [CrossRef]
17. P. Thueler, I. Charvet, F. P. Bevilacqua, M. S. Ghislain, G. Ory, P. Marquet, P. Meda, B. Vermeulen, and C. D. Depeursinge, “In vivo endoscopic tissue diagnostics based on spectroscopic absorption, scattering, and phase function properties,” J. Biomed. Opt. 8(3), 495–503 (2003). [CrossRef]
18. B. Yu, H. L. Fu, and N. Ramanujam, “Instrument independent diffuse reflectance spectroscopy,” J. Biomed. Opt. 16(11), 1–12 (2011). [CrossRef]
19. P. Lemaillet, C. C. Cooksey, J. Hwang, H. Wabnitz, D. Grosenick, L. Yang, and D. W. Allen, “Correction of an adding-doubling inversion algorithm for the measurement of the optical parameters of turbid media,” Biomed. Opt. Express 9(1), 55–71 (2018). [CrossRef]
20. X. U. Zhang, A. L. Post, D. J. Faber, T. G. van Leeuwen, and H. J. C. M. Sterenborg, “Single fiber reflectance spectroscopy calibration,” J. Biomed. Opt. 22(10), 1 (2017). [CrossRef]
21. W. Blondel, A. Delconte, G. Khairallah, F. Marchal, A. Gavoille, and M. Amouroux, “Spatially-resolved multiply-excited autofluorescence and diffuse reflectance spectroscopy: Spectrolive medical device for skin in vivo optical biopsy,” Electronics 10(3), 243 (2021). [CrossRef]
22. V. Colas, W. Blondel, G. Khairallah, C. Daul, and M. Amouroux, “Proposal for a skin layer-wise decomposition model of spatially-resolved diffuse reflectance spectra based on maximum depth photon distributions: A numerical study,” Photonics 8(10), 444 (2021). [CrossRef]
23. B. Aernouts, E. Zamora-Rojas, R. V. Beers, R. Watté, L. Wang, M. Tsuta, J. Lammertyn, and W. Saeys, “Supercontinuum laser based optical characterization of intralipid® phantoms in the 500-2250 nm range,” Opt. Express 21(26), 32450–32467 (2013). [CrossRef]
24. S. Prahl, “A compendium of tissue optical properties,” Http://omlc.ogi.edu/spectra/.
25. A. V. Prokhorov, S. N. Mekhontsev, and L. M. Hanssen, “Monte carlo modeling of an integrating sphere reflectometer,” Appl. Opt. 42(19), 3832–3842 (2003). [CrossRef]
