Translator Disclaimer
19 September 2018 Validation of optical properties quantification with a dual-step technique for biological tissue analysis
Author Affiliations +
Abstract
To approach wide-field optical properties quantification in real heterogeneous biological tissue, we developed a Dual-Step setup that couples a punctual diffuse reflectance spectroscopy (DRS) technique with multispectral imaging (MSI). The setup achieves wide-field optical properties assessment through an initial estimation of scattering with DRS, which is used to estimate absorption with MSI. The absolute quantification of optical properties is based on the ACA-Pro algorithm that has been adapted both for DRS and for MSI. This paper validates the Dual-Step system not only on homogeneous Intralipid phantoms but also on a heterogeneous gelatine phantom with different scattering and absorbing properties.

1.

Introduction

Spectral Imaging technology has arisen in the biomedical field as a powerful analytical tool because of its capability in providing spectral information in a large field of view (FOV), which considerably improves the accuracy and speed of clinical diagnosis. It combines imaging and spectroscopy to obtain simultaneous two-dimensional (2-D) spatial and spectral information that are saved in the so-called “three-dimensional (3-D) data cube.” An interesting review1 details the existing spectral imaging techniques with their advantages and limitations.

The spectral imaging approach addresses two main clinical requirements: the consideration of a wide spatial FOV, which optimally reduces the measurement time of the entire region of interest, and a noncontact setup modality that allows no-touch and sterile measurements of sensitive samples, such as injured tissue.

Yet, the quantification of optical properties to achieve optimal accuracy of diagnosis remains the main challenge to be solved. Indeed, the use of nonpunctual illumination leads to a mixed effect of absorption and scattering properties that cannot be separated from the detected diffuse reflectance. For a given reflectance value, several optical properties couples can explain the measurement. Thus, no quantification of scattering and absorption is possible. Different strategies have been used in the literature to analyze these detected multispectral reflectance data.

A first approach consists of working on raw data directly. This is the case of narrow band imaging, which is a form of multispectral imaging (MSI), first described as an “electronic chromoendoscopy” technique.2 Even though this technique improves visual perception of pathological tissue with respect to white-light endoscopy, it achieves no quantification of optical properties. The accuracy of the diagnosis is, therefore, subjective as it mainly depends on the user’s experience.

Remarkable work has been carried out in the development of a snapshot spectral imager that is able to acquire fast images of reflectance and autofluorescence, appropriate for clinical analysis of the oral cavity.3 The technique works between 471 and 667 nm and covers a 3- to 5-cm2 FOV. The results obtained are limited to spectral intensity ratios and fluorescence spectra that are not corrected from absorption distortions. In other words, the lack of absolute quantification of optical properties limits the accuracy of the diagnosis.

A second approach supposes one of the optical properties (mainly the scattering coefficient) with some a priori knowledge, to access the second one (i.e., absorption coefficient). As a consequence of making suppositions on scattering, only relative differences of absorption between different media can be derived. This is the case of various studies47 that suppose a unique a priori homogeneous scattering property to obtain relative absorption wide field maps. A similar procedure for absorption difference quantification has been proposed8 through the supposition of tissue water content and scattering-dependent variables that define a scattering model. However, these techniques are limited by the accuracy of the scattering model and to homogeneous samples. As soon as some scattering heterogeneities exist, absorption estimation error undoubtedly increases. Moreover, the structural information difference, given by these scattering heterogeneities carrying valuable diagnosis information, is overlooked.

A new approach of wide-field quantification, established in the spatial frequency domain and designated as spatial frequency domain imaging (SFDI), has been developed.9 The Fourier transform of the diffusion model10 is used to derive an expression of diffuse reflectance in terms of illumination spatial frequency.9 This model is compared with the measurements of diffuse reflectance resulting from different illumination spatial frequencies to allow the discrimination and wide-field quantification of both optical properties. Hence, maps of physiological properties in turbid media can be obtained.9,11 The advantages of this technique are many: quantitative scattering and absorption coefficients, inspection of a large FOV, and fast measurement and treatment, facilitating its use under clinical conditions. Therefore, SFDI can be considered as a reference technique of wide field quantification. The spatial resolution of optical properties estimated with SFDI has been analyzed using phantoms with different sizes of surface inclusions.12 It was shown that the quantification of optical properties depends on the inclusion size. Inclusions had to present a scattering contrast of at least 28% to achieve a resolution of 1.25 mm. Likewise, the absorption contrast had to be of at least 29% to detect inclusions of 5-mm size. This gives an idea of the spatial resolution of SFDI quantification, limited by the mixed effect of neighboring pixels on the detected diffuse reflectance. For the last years, the technique has also been investigated toward a depth resolution of quantified optical properties.13,14 Even if the technique seems very promising, absolute quantitative and spatially resolved optical properties estimation is not optimal.

A fourth approach to achieve a wide field quantitative map consists of scanning, across the entire region, quantitative punctual spectroscopy measurements. Indeed, contrary to spectral imaging techniques, the diffuse reflectance spectroscopy (DRS) technique is capable of quantifying optical properties but is limited to a small FOV due to its punctual measurement modality. DRS quantification is based on the spatial confinement of punctual illumination and diffuse reflectance detection, which allows the estimation of the optical path length. Hence, the optical properties can be separated and quantified.15 To obtain spatially quantified 2-D maps several Whiskbroom spectral imaging setups have been developed for in vivo imaging, with a compromise between spatial resolution and acquisition time.16,17 The Whiskbroom techniques spatially scan punctual illumination all throughout the sample. For each point, the diffuse reflectance is measured through concentric fibers placed at a single distance from the source. This is equivalent to a scanning DRS technique. However, to obtain a smooth “3-D data” cube, this modality requires complex control of the spatial scan and a considerable amount of time, which is impractical for clinical application. The acquisition time challenge has been addressed with the development of a multiplexing optical fiber-based tool under contact modality for ex vivo samples.18 This technique corresponds to a parallel Whiskbroom DRS scan that achieves optical properties quantification maps in a single acquisition, largely reducing acquisition time. However, the system is affected by contact modality limitations such as pressure-dependent coupling.

The accurate quantification of both optical properties on an image obtained with wide field MSI of a heterogeneous medium (such as tissue), therefore, remains a scientific challenge to improve the accuracy of clinical diagnosis. To address this issue, we propose an alternative method, namely the Dual-Step technique, which combines the quantification capacity of DRS with the spatial wide FOV offered by MSI. On the one hand, we promote the use of a Non-Contact DRS instrumental setup to perform quantitative and noncontact measurements convenient for clinical application. On the other hand, we use the MSI technique to allow fast wide field measurements.

Similar instrumental combinations have been developed for endoscopy offering wide field imaging and punctual Diffuse, Raman, and fluorescence spectroscopies, to improve the sensitivity and specificity of cancer diagnosis.19,20 The punctual DRS measurement is used to measure spectral reflectance differences of individual sites in the imaged FOV; however, no optical property quantification maps are considered.

The principle of the Dual-Step technique that we propose goes a stage further with the wide-field MSI absorption quantification. The latter is based on estimations of common scattering properties with punctual DRS measurements on specific zones of the imaged sample. These DRS estimated scattering coefficients are used to quantify absorption over the whole image. Contrary to the techniques that suppose a single modeled scattering coefficient and hinder the absolute quantification of the absorption coefficient,47 our system intends to quantitatively estimate both absolute optical properties over a wide FOV. To do so, we have developed the exhaustive adaptive calibration algorithm (ACA)-Pro that has been adapted both for punctual DRS measurements21 and to wide field MSI images.

We describe the instrumental Dual-Step setup in Sec. 2, outline the quantitative ACA-Pro approach for MSI in Sec. 3, and show the results that validate the technique in Sec. 4. Finally, we discuss the results and conclude with the potential of the Dual-Step approach in Sec. 5.

2.

Measurement Set-Up and Data Preprocessing

As a proof of concept, we implemented the Dual-Step approach as the sequential acquisition of Non-Contact DRS and MSI measurements on the same sample via a translation table (see Fig. 1).

Fig. 1

Scheme and picture of the Dual-Step technique instrument.

JBO_23_9_096002_f001.png

The Non-Contact DRS is detailed in Sec. 2.1. The spatial resolution of the DRS probe, given by the concentric fibers placed at several distances from the source, allows the estimation of scattering and absorption separately. For the proposed Dual-Step technique, we focus on the scattering coefficient DRS estimation only, which is then used by the MSI system, described in Sec. 2.2, to obtain a 2-D absorption map.

2.1.

Non-Contact Diffuse Reflectance Spectroscopy

The Non-Contact DRS system that we use is detailed in a previous study21 and its design preserves the same geometrical dimensions of the initial Contact DRS system described. Measurements performed with both systems inspect the same sample volume to allow their direct comparison: Contact DRS is considered as the reference measurement system to validate our Non-Contact DRS results.

The illumination is provided by a Tungsten Halogen (T–H) Lamp HL2000 Ocean Optics and the detection is ensured by a QE65000 Ocean Optics spectrophotometer, cooled to 15°C to reduce dark noise. The spectral working range of the entire instrument lies between 470 and 880 nm. The probe features a central illumination fiber (Ø=500  μm) and six detection rings [each composed of seven concentric fibers (Ø=100  μm)] at six different distances, ranging from 300 to 2488  μm (center to center). The seven fibers of each detection ring are clustered in a specific bundle referred to as F6 to F1, with F6 and F1 measuring the closest and furthest ring to the source, respectively. F7 is a specific fiber that enables a direct measurement of the illumination spectrum. A step-wise measurement of each Fi signal by the spectrometer provides the diffuse reflectance at all source-detector distances individually and sequentially. The noncontact modality of the measurement is ensured by placing an achromatic doublet pair between the probe and the sample as shown in Fig. 1(a) (right). The specific treatment of these Non-Contact DRS measurements has been precisely described21 and will not be detailed here.

2.2.

Multispectral Imaging

The MSI instrumental setup is shown in Fig. 1(a) (left). The system makes use of a 12-bit monochrome VGA PixelFly CCD detector, matrix of 640×480  pixels with a pixel size of 9.9  μm, which allows good spatial resolution. A Xenoplan 2.0/28 objective lens ensures a FOV of 67×50  mm2 (corresponding to a magnification G=0.095) at the focal plane, where the sample is placed. The lamp is a KL2500 LCD Schott source set for a broad spectral distribution (3000 K of temperature). The system sequentially acquires data with a limited number of spectral bands leading to a good tradeoff between spectral sampling and acquisition time. Namely, the source is coupled with a filter wheel that defines four specific illumination wavelengths: 500, 550, 600, and 700 nm having a 10-nm FWHM. The illumination of the sample is achieved through the use of a six-LED illumination ring that is inclined to avoid specular reflection of imaged flat samples.

The multispectral images, noted S, are preprocessed according to a background (Soffset) subtraction and corrections of integration time (t) and flat field (FF). The flat-field image (Sflat) is acquired with a Spectralon® standard SR99%. Background Soffset and signal S images are acquired with the same integration time t. The preprocessed or corrected image is referred to as Sc and calculated according to

Eq. (1)

FF=SflatSoffsetmax(SflatSoffset)Sc=SSoffsetFF·t

Depending on the required signal-to-noise ratio, spatial resolution, and calculation time, the Sc image is averaged on neighboring pixels according to Eq. (2), with Tx=Ty=1 leading to the primary image spatial resolution

Eq. (2)

Sa=x=1Txy=1TySc(x,y)Tx·Ty.

Finally, the resulting averaged image Sa is normalized by an indirect measurement of the illumination intensity, namely the intensity of a fixed plastic corner in the FOV (see Fig. 1, left). This is done to correct the fluctuations of the illumination intensity from one measurement to another and thereby ensure repeatability and accuracy of the quantification analysis. The resulting normalized image is noted as   SN.

2.3.

Phantom set and Dual-Step Acquisition Protocol

2.3.1.

Phantom set

The data used in this study are acquired on phantoms made of deionized water, different concentrations of black ink (Rotring) to mimic absorption, and different concentrations (%) of a fat emulsion Intralipid® (IL) 200 mg to mimic Mie scattering.22 Ten reference phantoms having two different reference μs scattering properties (ILref=1%, ILref=1.5%) and five different reference μa absorption properties (μa,ref=0.44,1.15,3.12,6.61,11.45  cm1 at 550 nm), are considered in this study.

2.3.2.

Acquisition Protocol

The Dual-Step technique performs a sequential measurement acquisition with Non-Contact DRS and MSI, coupled via a translation stage [see Fig. 1(b)], of the sample placed at the common focal plane. The whole system is isolated inside a black-cardboard structure and reflecting components are covered with a black cloth to ensure no uncontrolled reflections. Thus, DRS and MSI data are acquired under the same ambient darkness conditions as the ones used to measure FF, inherent parasite reflections, and offset signals used in the preprocessing of each modality of the Dual-Step system.

3.

Quantification Processing

The Dual-Step technique relies on the estimation of the scattering coefficient μs with Non-Contact DRS measurements, used to further estimate the absorption coefficient μa with MSI measurements.

3.1.

Scattering Coefficient Estimation by Non-Contact Diffuse Reflectance Spectroscopy

Both absorption and scattering coefficients estimation through Non-Contact DRS are performed according to a Monte Carlo (MC) look-up table (LUT)-based approach and the adaptive calibration algorithm ACA-Pro, detailed in a previous study.21

In this work, we define the depth of field (DOF) as the focus range achieving an acceptable spatial resolution that allows accurate estimation of optical properties with Non-Contact DRS. As the scattering property is estimated with the closest fibers to the source, its estimation is more sensitive to the focus range than that of the absorption property, estimated with further fibers.21 Therefore, the Non-Contact DRS scattering estimation error is used to determine the DOF of the Dual-Step system. This is done according to the average error calculated between the reference μs (μs,ref) of a phantom fixed at the focal plane; and μs estimations (μs^) of the same phantom placed at different object planes.

An acceptable error of 4.1% on μs^ is achieved for a DOF of 1.2 mm, as shown in Fig. 2.

Fig. 2

Estimated averaged μs^ (blue line) using a fixed focused reference phantom (black dotted line) and measurements at object plane range of 1.2 mm. The average μs^ error, Er=4.1%, considers the entire working spectrum.

JBO_23_9_096002_f002.png

However, this DOF (1.2 mm) leads to an error on absorption estimation of about 27%, which is unusable in the real context of biological samples. Thus, the Dual-Step technique uses Non-Contact DRS to estimate only the scattering coefficient, noted μs^. All the results in this work are obtained with a scattering coefficient μs^ estimated by Non-Contact DRS.

3.2.

Absorption Coefficient Estimation by Multispectral Imaging

As explained in the introduction, MSI cannot quantify both absolute optical properties. However, quantification of one optical property is possible if the other is known. Consequently, the Dual-Step technique uses the scattering μs^ estimations of Non-Contact DRS to allow absorption estimation with MSI. For the latter, we use an LUT specific to MSI, built according to the diffusion model described in a previous study.10 This model has been used for various MSI techniques4 including the reference SFDI.9 The model expresses diffuse reflectance R, originated from a uniform illumination, as a function of the wavelength and the albedo a, which combines both optical properties μa and μs. We use the equation defined by Cuccia et al.9 under the assumption of a planar illumination

Eq. (3)

R=3Aa(3(1a)+1)(3(1a)+3A),
where A is the proportionality constant A=1Reff2(1+Reff), Reff is the effective reflection coefficient Reff0.0636n+0.668+0.71n+1.44n2, n is the refractive index equal to 1.33 for water-based phantoms and equal to 1.37 for tissue samples, and a is the reduced albedo a=μs/(μs+μa).

Figure 3 (bottom left) shows the modeled diffuse reflectance RLUT in terms of optical properties (μs, μa) saved in an LUT.

Fig. 3

Flowchart describing the absorption quantification method of MSI. On the bottom left, RLUT proper of MSI computed with the diffusion model [Eq. (3)]. A preliminary calibration step (in blue) corrects the treated measure SN, (in pink) from the instrumental response. The comparison of the calibrated signal (RCF,) to the entire LUT (RLUT) via a minimization procedure (in black) leads to the identification of the possible optical properties couples μs-μa (red line in bottom right plot). With the given μs^ previously estimated with DRS, estimation of absorption μa˜ is possible with MSI (green arrow).

JBO_23_9_096002_f003.png

The general method for MSI absorption quantification relies on a calibration process (blue path in Fig. 3) that uses the measurement of a single reference phantom with known optical properties μs,ref and μa,ref. The measured reflectance signal Sref of a given reference phantom is treated according to Eqs. (1) and (2), respectively. The comparison of SN,ref with RrefLUT (modeled reflectance for the same reference phantom) determines a correction ratio CF. This CF calibrates the treated signal SN of a measured unknown sample (having unknown optical properties) from the instrumental response. The resulting calibrated signal RCF, is compared with the entire RLUT to identify the possible couples of μsμa optical properties of the unknown phantom. These couples correspond to the minimal difference between RCF and RLUT (red line of bottom right plot in Fig. 3). From these possible optical properties μsμa couples, and for a given μs^, previously estimated with DRS, the absorption coefficient μa˜ can be determined with MSI (green arrow of Fig. 3).

The improvement in the wide-field absorption estimation accuracy through the use of Non-Contact DRS estimated μs^, compared to theoretically modeled scattering properties,47 has been quantified and validated in a previous study.23

3.3.

Adaptive Calibration Algorithm Procedure for Absorption Quantification

The model expressed by Eq. (3) makes the approximation that scattering properties dominate over absorption properties. One should keep in mind that this is not always the case for skin, especially for dark skin phototypes and at λ<600  nm due to the high absorption of melanin and hemoglobin.24

Randeberg et al.25 approximated the error between the diffusion model and a Monte Carlo simulation for λ=400 to 850 nm and reduced it through the application of a scaling factor, as a function of the absorption coefficients, on the diffusion model.

In the first stage of this work, we did not consider the influence of optical properties on the calibration factor CF and used a single reference phantom defining a single CF (Sec. 3.2). However, this CF highly depends on the scattering coefficient (IL  ref) and the absorption coefficient (μa,ref) of the reference phantom (see Fig. 5). Hence, to further improve absorption estimation results we integrate an optical property-dependent correction through the MSI ACA-Pro approach that makes use of several reference phantoms with different known reference optical properties (μs,ref, μa,ref). These reference phantoms are used to obtain a reference base of various CFs, noted as CF(μs,ref,μa,ref).

Figure 4 shows the flowchart of the ACA-Pro μa quantification process applied to the MSI system in combination with Non-Contact DRS, according to the Dual-Step technique. The ACA-Pro algorithm makes use of an optimized number of reference phantoms having known optical properties (μs,ref, μa,ref) that cover the optical range of study. The measured and treated signals SN,ref of these phantoms are compared with their corresponding modeled diffuse reflectance RrefLUT and their ratio is calculated to create the MSI CF reference base CF(μs,ref,μa,ref).

Fig. 4

Flowchart describing the absorption quantification (μa^) method of MSI with the ACA-Pro approach (in blue), based on the Non-Contact DRS scattering estimation μs^.

JBO_23_9_096002_f004.png

The aim of ACA-Pro is to calculate the optimal CF, noted as CFopt(μs^,μa^), achieving optimal MSI estimation of absorption property μa^. This process starts with the estimation of μs^ with Non-Contact DRS.21

The DRS μs^ estimation defines a reduced set (group) of possible CFs from the entire MSI reference base CF(μs^,μa,ref), i.e., μs^ defines the CF for all the couples μs^μa,ref. If μs^ is included in the reference base CF(μs,ref,μa,ref), i.e., μs^=μs,ref, then the CF(μs^,μa,  ref) set is directly selected. In practice, this is not possible due to estimation error of Non-Contact DRS.21 In other words, even if the theoretical scattering property of the unknown phantom, μs,theo, is equal to μs,ref, real Non-Contact DRS measurements of the unknown phantom will estimate μs^μs,ref. As μs^ is not included in the reference base CF(μs,ref,μa,ref), i.e., μs^μs,ref, the CF(μs^,μa,ref) set is found through interpolation.

The selected CFs from the CF(μs^,μa,ref) set are used to calculate the corresponding MSI absorption estimations μa˜ and related estimation difference |μa˜μa,ref|. ACA-Pro selects the optimal CF, CFopt(μs^,μa^), that estimates optimal μa˜, noted μa^, with minimal difference.

This method supposes an ideal situation in which the theoretical absorption property of the unknown phantom, noted as μa,theo, is equal to one of the μa,ref of the CF(μs^,μa,ref) set, i.e., μa,theo=μa,ref. This ideal case is used to test the ability of ACA-Pro in estimating μa˜ corresponding to a given CF and selecting CFopt(μs^,μa^) that achieves minimal difference |μa˜μa,ref| and hence, minimal estimation error |μa^μa,theo| (Sec. 4.1.1). This resulting minimal estimation error is only due to measurement acquisition noise and corresponds to the ideal ACA-Pro performance.

In a more realistic context, the unknown phantom does not have the same absorption properties as any of the reference phantoms in the CF(μs^,μa,ref) set, i.e., μa,theoμa,ref. Therefore, ACA-Pro includes an iterative interpolation method of the CF. Given a single random CF from the CF(μs^,μa,ref) set estimating an initial μa,0˜, the CF is interpolated (within the set) to CF(μs^,μa,0˜), estimating μa,1˜. Hence, the CF(μs^,μa,i+1˜) is iteratively updated according to the previous estimation μa,i˜ until convergence, i.e., μa,i+1˜=μa,i˜. We note the final interpolated CF as CFinterp, which corresponds to CFopt(μs^,μa^) and estimates μa^ with minimal estimation error |μa^μa,theo|. Experimental results are shown in Sec. 4.1.2.

4.

Results

The results presented in Sec. 4.1 validate the use of MSI ACA-Pro for absorption quantification with the Dual-Step technique on homogeneous phantoms. Section 4.1.1 defines the accuracy limit of ACA-Pro from the results obtained for the ideal situation: μa,theo of the unknown phantom is equal to one of the μa,  ref of the CF(μs^,μa,ref) set, i.e., μa,theo=μa,ref. Section 4.1.2 presents an example of the general case: μa,theo of the unknown phantom is different from all μa,ref of the CF(μs^,μa,ref) set, i.e., μa,theoμa,ref. For this general case, the ACA-Pro estimation follows the flow chart shown in Fig. 4.

The influence of the DOF on the result accuracy is discussed in Sec. 4.2. In Sec. 4.3, we explain how the method is extended to heterogeneous phantoms and focus on the Non-Contact DRS scattering segmentation step preliminary to the wide-field MSI absorption estimation.

4.1.

Validation of ACA-Pro for Absorption Quantification with the Dual-Step Technique on Homogeneous Phantoms

A 10-reference base CF(μs,ref,μa,ref) is built with measurements taken on 10 different reference phantoms (described in Sec. 2.3.1) and four wavelengths λ (500, 550, 600, and 700 nm). Figure 5 shows the resulting 10-reference base CF(μs,ref,μa,ref) according to λ and reference optical properties. Remark the dependence of CF(μs,ref,μa,ref) on μa,ref, originating from the lower accuracy of the diffusion model at high absorption values. This dependence justifies the use of the μa-based ACA-Pro calibration strategy.

Fig. 5

Ten-reference base CF(μs,ref,μa,ref) according to illumination λ built with ten-reference phantoms having two different scattering properties (ILref=1  % and ILref=1.5%) and five different absorption properties μa,ref=0.44,1.15,3.12,6.61, and 11.45  cm1 at 550 nm.

JBO_23_9_096002_f005.png

4.1.1.

Ideal case: choice of CFopt(μs^,μa^), for unknown phantoms included in the CF reference base

In this section, we analyze the MSI μa-based ACA-Pro performance in the absorption estimation μa˜ and corresponding estimation error for any given CF. We also aim at defining the minimal error achievable with ACA-Pro applied to our measurement set-up. For this, we estimate the absorption of five unknown phantoms having theoretical properties equal to those of five reference phantoms used for the 10-reference base CF(μs,ref,μa,ref): ILtheo=1%, μa,theo=0.44,1.15,3.12,6.61, and 11.45  cm1 at 550 nm. The absorption coefficient of all five measured unknown phantoms is estimated with the entire 10-reference base CF(μs,ref,μa,ref), i.e., 10 CFs.

Moreover, we look at the influence of μs^ on the absorption estimation μa˜ by comparing the results obtained with the “wrong” CF(μs^=μs,refIL=1.5%μs,theoIL=1%,μa,ref) and “correct” CF(μs^=μs,refIL=1%=μs,theoIL=1%,μa,ref) reduced sets.

For each measured unknown phantom, the resulting absorption difference |μa˜  μa,ref| of μa˜ estimated with all 10 CFs of the 10-reference base CF(μs,ref,μa,ref) is shown in Fig. 6 for 550 nm. Notice the expected influence of μs^ on μa˜: the absorption difference is higher for the “wrong” CF(μs^=μs,refIL=1.5%μs,theoIL=1%,μa,ref) than for the “correct” CF(μs^=μs,refIL=1%=μs,theoIL=1%,μa,ref). The μa,ref values of the CF also have a coherent effect on the absorption difference: the closer μa,ref is to μa,theo, the lower the difference is. As expected, the minimal absorption difference is obtained with the CF corresponding to the optical properties of the unknown phantom CF(μs^=μs,refIL=1%=μs,theoIL=1%,μa,ref=μa,theo). ACA-Pro converges on this CF and defines it as CFopt(μs^,μa^) that estimates μa^ with minimal error |μa^μa,theo| corresponding to measurement acquisition noise only.

Fig. 6

Spatially averaged |μa˜μa,ref| difference at 550 nm of unknown phantoms having common ILtheo=1% and different μa,theo (shown in the x-axis), calibrated with the 10-reference base CF(μs,ref,μa,ref) (with reference optical properties stated in the legend) and the optimal CFopt(μs^,μa^) chosen by MSI ACA-Pro.

JBO_23_9_096002_f006.png

This optimized estimation procedure, i.e., definition of CFopt(μs^,μa^) is performed for each pixel value in the image to estimate wide field μa^. Figure 7 shows an example of wide-field μa^ for the lowest absorbing unknown phantom considered (μa,theo=0.44  cm1 at 550 nm). Notice that all μa^ pixels are in accordance with the theoretical value μa,theo and present an average error of 4% for all considered λ.

Fig. 7

Wide-field multispectral μa^ quantification for an unknown phantom having μa,theo specified by the central figure and marked by the black arrows on the colorbars. These results are obtained with Tx=Ty=10  pixels [see Eq. (2)]. The image size is 33×33  pixels corresponding to a FOV of 34.4×34.4  mm2.

JBO_23_9_096002_f007.png

Figure 8 shows the minimal spatially averaged μa errors of μa^ estimated with CFopt(μs^,μa^) for all unknown phantoms μa,theo and λ. Notice that all minimal mean errors are <5.5% (see 600 nm). These results comprehensively validate the wide-field   μa calibration procedure of MSI ACA-Pro for the range of optical properties considered.

Fig. 8

Spatially averaged |μa^μa,theo| estimation errors of unknown phantoms having common ILtheo=1% and different μa,theo calibrated with corresponding CFopt(μs^,μa^) for all considered λ. The error is inferior to 5.5% (see 600 nm). The green circle shows the error to be compared with that achieved with the interpolation strategy (Sec. 4.1.2).

JBO_23_9_096002_f008.png

4.1.2.

General case: choice of CFopt(μs^,μa^), for unknown phantoms not included in the CF reference base. Interpolation strategy

We consider a general case to illustrate and analyze the iterative interpolation strategy of MSI μa-based ACA-Pro that calculates CFinterp=CFopt(μs^,μa^).

In this case, ACA-Pro uses a reduced two-reference base CF(μs,ref,μa,ref), built only with two reference phantoms having ILref=1  % and μa,ref=1.15 and 6.61  cm1 at 550 nm. The measured unknown phantom has ILtheo=1% and μa,theo=3.12  cm1 at 550 nm. We use the real Non-Contact DRS μs estimation, μs^μs,  ref (refer to Sec. 3.3). Notice that the μa,theo=3.12  cm1 of the unknown phantom does not correspond to any of the μa,ref in the two-reference base CF(μs,ref,μa,ref), i.e., μa,theoμa,ref.

A first random CF from the CF(μs^,μa,ref) set estimates μa,0˜ that defines the next CF, CF(μs^μa,0˜), through interpolation of the two-reference base/set. CF(μs^μa,0˜) then estimates μa,1˜. After various iterations   μa  i+1˜ stabilizes and converges to μa,  i˜, becoming μa^, estimated with CFopt(μs^,μa^) (refer to Fig. 4).

Figure 9 shows the wide-field μa˜ estimations at 550 nm with the two CFs of the reference base: CFμa,ref=1.15  cm1 and CFμa,ref=6.61  cm1. These μa˜ estimations are compared with the μa^ estimated with the interpolated CFinterp=CFopt(μs^,μa^). The average wide-field |μa˜  μa,theo| estimation errors obtained with all three CFs are shown at the bottom right of Fig. 9. Remark that CFinterp achieves an averaged minimal wide-field |μa^μa,theo| estimation error of 2.6%, which is comparable with that obtained with the ideal CF(μs^,μa,ref=μa,theo): error=1.7% (see Sec. 4.1.1, green circle of Fig. 8 at 550 nm for the same μa,theo). This result validates the iterative interpolation strategy of μa-based ACA-Pro applied to MSI.

Fig. 9

Wide-field μa estimation at 550 nm of an unknown phantom having μa,theo=3.12  cm1 calibrated with the two-reference base CF(μs,ref,μa,ref) built with two-reference phantoms only having μa,ref=1.15 and 6.61  cm1. (a) μ˜a estimation with CF corresponding to μa,ref=1.15  cm1. (b) μ˜a estimation with CF corresponding to μa,ref=6.61  cm1. (c) μa^ estimation with CFinterp=CFopt(μs^,μa^). (d) Mean |μa˜μa,theo| estimation errors with CFμa,ref=1.15  cm1, CFμa,ref=6.61  cm1, and |μa^μa,theo| for CFinterp. The results given in this figure are obtained with Tx=Ty=10  pixels [see Eq. (2)]. The image size is 33×33  pixels corresponding to an FOV of 34.4×34.4  mm2.

JBO_23_9_096002_f009.png

An optimal choice of the minimal amount of reference phantoms used to build the general reference base CF(μs,ref,μa,ref) is necessary to optimally benefit from the interpolation strategy of ACA-Pro. This choice directly depends on the range of interest of optical properties.

Contrarily to what is expected, Figs. 7 and 9 show that the estimated absorption map is not spatially homogeneous for the measured homogeneous unknown phantoms. This artifact is due to the fact that, for the given results, a unique CF, computed in the middle of the image, is used to correct the entire image. Therefore, a spatially local CF correction will be implemented in prospective work to correct this drawback.

4.2.

Influence of the Depth of Field on the Absorption Estimation Accuracy

We mentioned in Sec. 3.1 that an error of 4.1% on μs^ estimation is achieved for a DOF of 1.2 mm with Non-Contact DRS. Furthermore, we quantified the error on μa (27%) and concluded that Non-Contact DRS is only quantitative at the focal plane for both optical parameters (1% error for μs and 2.4% for μa).26

In this work, we quantify the influence of the DOF on the μa^ estimation accuracy by MSI. In this system, the measured reflectance depends linearly on the distance between the sample and the objective lens. We analyze this effect with a test unknown phantom having the same optical properties as the single reference phantom (μa,ref=0.44  cm1 at 550 nm and ILref=1%), by estimating μa^ according to the procedure described in Sec. 3.2. The scattering coefficient of the phantom, estimated by Non-Contact DRS (μs^) and the calculated MSI CFopt(μs^,μa^) are obtained at the focal plane. This CFopt(μs^,μa^) is then used to estimate μa^ with MSI on the same phantom placed at +/3.5  mm around the focal plane. The resulting MSI error μa^μa,theo, due to the reflectance intensity variation related to the plane height difference, is shown in Fig. 10. As expected, the μa^ estimation decreases with height, since the reflectance intensity increases as the sample approaches the source. Notice that for an error of 4.1%, the DOF is around 2.5 mm. Thus, because the DOF is more constraining for μs^, we consider only a DOF of 1.2 mm.

Fig. 10

μa^μa,theo error originated from the reflectance variation at different planes (+/3.5  mm) at 550 nm.

JBO_23_9_096002_f010.png

Therefore, the final DOF of the Dual-Step technique is determined by the most constraining modality, i.e., the Non-Contact DRS. The defined DOF of 1.2 mm obtains an error inferior to 4.1% for Non-Contact DRS μs^ estimation (see Fig. 2). The effect of the 4.1% Non-Contact DRS μs error on the μa^ MSI estimation (5.5% error), has been quantified as a maximum error variation of ±4.03% (for all studied wavelengths: 500, 550, 600, and 700 nm). Further improvement in the μa^ MSI estimation will be achieved in future work by integrating a height-calibration (based on Fig. 10) after a relief analysis of the sample.

4.3.

Optical Properties Quantification Procedure on Heterogeneous Phantoms

After validating with homogeneous phantoms, and to approach clinical measurement of heterogeneous biological samples, we tested our technique with a heterogeneous phantom. For this, we used a homemade gelatine-based phantom with different optical properties shown in Fig. 11.

Fig. 11

Homemade gelatine-based heterogeneous phantom composed of different unknown optical properties μa and μs.

JBO_23_9_096002_f011.png

Absorption of each zone depends on the different concentrations of VAHINÉ red ink or HP-printer blue ink used. Scattering depends on the amount of scattering milk and edible gelatine used. Blue shapes (circle and corners) are made of the same milk-gelatine mixture and two different concentrations of blue ink defined as “dark” for the bottom left corner and the bottom right circle and as “light” for the top right corner (see Fig. 11). The pink background is made of a different milk-gelatine mixture and of a low concentration of red ink. Note that, because the blue shapes have been molded before the pink ones, some overlap of the blue shapes with the pink background is present at border zones.

The 2-D MSI absorption estimation relies on the preliminary DRS estimation of local μs^. In case of heterogeneous phantoms, this scattering coefficient has to be measured in each zone of common scattering properties. Figure 11 shows the different zones that have been measured (crossed black circles on blue corners and pink background) with our reference Contact DRS. Figure 12 shows the Contact DRS μs and μa estimations, specific to the different zones and used as the optimal scattering estimation μs^ and the theoretical absorption μa,theo (see Fig. 12) in our MSI ACA-Pro algorithm. Note that the measured blue corners have the same μs^, which is lower to that of the pink background.

Fig. 12

Contact DRS estimated (a) μs and (b) μa properties of the different zones (background and light blue or dark blue shapes) on the heterogeneous gelatine-based phantom.

JBO_23_9_096002_f012.png

Two filters at 550 and 600 nm were chosen for absorption quantification with MSI as they highlight a clear difference between the two chromophores to be analyzed (blue and red ink, see Fig. 12). Images with these filtered lights are acquired and treated as explained in Sec. 2.2. Furthermore, the small specular reflections from the original image are filtered with an algorithm that accords the minimal value of a mask to the central pixel. Filtered images are shown in Fig. 13. Notice the presence of holes (due to air bubbles during the preparation of the phantoms); they will not be considered for quantification.

Fig. 13

Treated images of the heterogeneous phantom under 550 and 600 nm filtered illumination. The image size is 330×330  pixels corresponding to a FOV of 34.4×34.4  mm2.

JBO_23_9_096002_f013.png

The first step in the quantification of wide-field μa^ with MSI is to give a DRS μs^ estimation value to each pixel of the image. The μs^ wide-field map is built with an intensity-based segmentation procedure, based on the high contrast of the 600 nm image and the Otsu method. DRS μs^ values [Fig. 12(a)] are set for each segmented homogeneous zone (background and shapes) as shown in Fig. 14.

Fig. 14

Wide-field μs^ map for 550 and 600 nm. The image size is 330×330  pixels corresponding to an FOV of 34.4×34.4  mm2.

JBO_23_9_096002_f014.png

This quantitative multispectral μs^ map together with the 10-reference base CF(μs,ref,μa,ref) (refer to Fig. 5) is used to calibrate with CFinterp=CFopt(μs^,μa^) single pixels of images MSI ACA-Pro (refer to Sec. 3.3). This results in the optimal wide-field μa^ quantification shown in Fig. 15. Notice that these quantitative μa^ maps present values that are in accordance with μa,theo estimated with Contact DRS [Fig. 12(b)]. Therefore, these results ultimately validate the wide-field μa^ quantification of unknown heterogeneous samples, based on a priori built μs^ map and the MSI ACA-Pro calibration procedure.

Fig. 15

Wide-field μa^ at 550 and 600 nm shown with two different scales (left to right) to enhance visualization. The results given in this figure are obtained with Tx=Ty=3  pixels [refer to Eq. (2)]. The image size is 110×110  pixels corresponding to an FOV of 34.4×34.4  mm2.

JBO_23_9_096002_f015.png

The analysis of quantitative μa^ maps at 550 and 600 nm is taken a step further with the calculation of blue and red ink concentrations. This is based on linear combinations of spectral unmixing that express the estimated μa^ at each wavelength (μa,550 and μa,600) as a linear combination of the relative extinction coefficient ϵR of red and blue inks multiplied by their corresponding concentrations Cred and Cblue [Eq. (4)]. In this specific case, no theoretical quantification of ϵR is available so it has been replaced by the reference Contact DRS μa signatures of the background and the light blue zones [refer to Fig. 12(b)]

Eq. (4)

μa,550=μa,red,550·Cred,550+μa,550·Cblueμa,600=μa,red,600·Cred,600+μa,600·Cblue,
where Cred and Cblue are so calculated for each pixel to build the concentration maps shown in Fig. 16. Notice the correct distribution of red ink on the pink background mixture and the blue ink on the shapes. Mixtures of red and blue ink are present at the border and the surface of shapes, as expected. Areas corresponding to holes are not quantitative because of the uncorrected shadow effect.

Fig. 16

Relative concentration of (a) red ink (b) blue ink, shown with two different scales for visual enhancement. The results are obtained with Tx=Ty=3  pixels [see Eq. (2)]. The image size is 110×110  pixels corresponding to a FOV of 34.4×34.4  mm2.

JBO_23_9_096002_f016.png

5.

Discussion and Conclusion

The Dual-Step imaging technique combines Non-Contact DRS and MSI to quantify wide-field absolute scattering and absorption properties through noncontact measurements. The scattering properties μs^ are estimated with local Non-Contact DRS measurements performed on the homogenous zones of interest. Based on these μs^ estimations, 2-D absorption maps μa^ are obtained with MSI. The Dual-Step technique relies on a thorough instrumental calibration performed by the developed ACA-Pro algorithm, which has been extended for wide-field absorption quantification in this work. This includes the use of a correction factor CF reference base, CF interpolation, and an additional factor scaling the illumination intensity fluctuations between images.

The wide-field quantification of the Dual-Step technique is validated with a range of well-characterized homogeneous Intralipid phantoms. Using Non-Contact DRS to estimate μs^ with an error inferior to 4.1%, the DOF of the technique lies within 1.2 mm (see Fig. 2). This leads to a maximum average μa^ relative error of 5.5%±4.03% for the studied range in this work with the MSI ACA-Pro calibration. Promising results have been obtained when using the Dual-Step quantitative approach to estimate optical wide-field optical properties on different abdominal ex-vivo human skin samples and in-vivo specific rat models (bicolored and inflammatory zones).26 Nevertheless, the DOF may be a limitation for biological samples where the sample curvature is higher than 1.2 mm leading to less accurate μs^ and μa^ estimation. Prospective work should focus in controlling the curvature of the sample to improve quantification accuracy on biological samples.

Disclosures

The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.

References

1. 

G. Li et al., “Review of spectral imaging technology in biomedical engineering: achievements and challenges,” J. Biomed. Opt., 18 (10), 100901 (2013). https://doi.org/10.1117/1.JBO.18.10.100901 JBOPFO 1083-3668 Google Scholar

2. 

K. Gono et al., “Appearance of enhanced tissue features in narrow-band endoscopic imaging,” J. Biomed. Opt., 9 (3), 568 –577 (2004). https://doi.org/10.1117/1.1695563 JBOPFO 1083-3668 Google Scholar

3. 

N. Bedard et al., “Multimodal snapshot spectral imaging for oral cancer diagnostics: a pilot study,” Biomed. Opt. Express, 4 (6), 938 –949 (2013). https://doi.org/10.1364/BOE.4.000938 BOEICL 2156-7085 Google Scholar

4. 

A. Bjorgan, M. Milanic and L. L. Randeberg, “Estimation of skin optical parameters for real-time hyperspectral imaging applications,” J. Biomed. Opt., 19 (6), 066003 (2014). https://doi.org/10.1117/1.JBO.19.6.066003 JBOPFO 1083-3668 Google Scholar

5. 

A. Basiri et al., “Use of a multi-spectral camera in the characterization of skin wounds,” Opt. Express, 18 (4), 3244 –3257 (2010). https://doi.org/10.1364/OE.18.003244 OPEXFF 1094-4087 Google Scholar

6. 

A. Vogel et al., “Using noninvasive multispectral imaging to quantitatively assess tissue vasculature,” J. Biomed. Opt., 12 (5), 051604 (2007). https://doi.org/10.1117/1.2801718 JBOPFO 1083-3668 Google Scholar

7. 

K. J. Zuzak et al., “Visible reflectance hyperspectral imaging: characterization of a noninvasive, in vivo system for determining tissue perfusion,” Anal. Chem., 74 (9), 2021 –2028 (2002). https://doi.org/10.1021/ac011275f ANCHAM 0003-2700 Google Scholar

8. 

S. L. Jacques, R. Samatham and N. Choudhury, “Rapid spectral analysis for spectral imaging,” Biomed. Opt. Express, 1 (1), 157 –164 (2010). https://doi.org/10.1364/BOE.1.000157 BOEICL 2156-7085 Google Scholar

9. 

D. J. Cuccia et al., “Quantitation and mapping of tissue optical properties using modulated imaging,” J. Biomed. Opt., 14 (2), 024012 (2009). https://doi.org/10.1117/1.3088140 JBOPFO 1083-3668 Google Scholar

10. 

L. O. Svaasand et al., “Tissue parameters determining the visual appearance of normal skin and port-wine stains,” Lasers Med. Sci., 10 (1), 55 –65 (1995). https://doi.org/10.1007/BF02133165 Google Scholar

11. 

S. Gioux et al., “First-in-human pilot study of a spatial frequency domain oxygenation imaging system,” J. Biomed. Opt., 16 (8), 086015 (2011). https://doi.org/10.1117/1.3614566 JBOPFO 1083-3668 Google Scholar

12. 

A. M. Laughney et al., “System analysis of spatial frequency domain imaging for quantitative mapping of surgically resected breast tissues,” J. Biomed. Opt., 18 (3), 036012 (2013). https://doi.org/10.1117/1.JBO.18.3.036012 JBOPFO 1083-3668 Google Scholar

13. 

R. B. Saager et al., “Method for depth-resolved quantitation of optical properties in layered media using spatially modulated quantitative spectroscopy,” J. Biomed. Opt., 16 (7), 077002 (2011). https://doi.org/10.1117/1.3597621 JBOPFO 1083-3668 Google Scholar

14. 

D. M. Burmeister et al., “Utility of spatial frequency domain imaging (SFDI) and laser speckle imaging (LSI) to non-invasively diagnose burn depth in a porcine model,” Burns, 41 (6), 1242 –1252 (2015). https://doi.org/10.1016/j.burns.2015.03.001 BURND8 0305-4179 Google Scholar

15. 

G. Zonios and A. Dimou, “Modeling diffuse reflectance from semi-infinite turbid media: application to the study of skin optical properties,” Opt. Express, 14 (19), 8661 –8674 (2006). https://doi.org/10.1364/OE.14.008661 OPEXFF 1094-4087 Google Scholar

16. 

S. F. Bish et al., “Handheld diffuse reflectance spectral imaging (DRSi) for in-vivo characterization of skin,” Biomed. Opt. Express, 5 (2), 573 –586 (2014). https://doi.org/10.1364/BOE.5.000573 BOEICL 2156-7085 Google Scholar

17. 

C. C. Yu et al., “Quantitative spectroscopic imaging for non-invasive early cancer detection,” Opt. Express, 16 (20), 16227 –16239 (2008). https://doi.org/10.1364/OE.16.016227 OPEXFF 1094-4087 Google Scholar

18. 

B. S. Nichols et al., “A quantitative diffuse reflectance imaging (QDRI) system for comprehensive surveillance of the morphological landscape in breast tumor margins,” PLoS One, 10 (6), e0127525 (2015). https://doi.org/10.1371/journal.pone.0127525 POLNCL 1932-6203 Google Scholar

19. 

H. Zeng et al., “Combining field imaging endoscopy with point analysis spectroscopy for improving early lung cancer detection,” Conf. Proc. IEEE Eng. Med. Biol. Soc., 2008 1849 –1850 (2008). https://doi.org/10.1109/IEMBS.2008.4649540 Google Scholar

20. 

L. Kan, “Multimodal optical spectroscopy and imaging for improving cancer detection in the head and neck at endoscopy,” (2012). Google Scholar

21. 

V. Sorgato et al., “ACA-Pro: calibration protocol for quantitative diffuse reflectance spectroscopy. Validation on contact and noncontact probe-and CCD-based systems,” J. Biomed. Opt., 21 (6), 065003 (2016). https://doi.org/10.1117/1.JBO.21.6.065003 JBOPFO 1083-3668 Google Scholar

22. 

S. L. Jacques, “Optical properties of Intralipid™, an aqueous suspension of lipid droplets,” (1998) http://omlc.ogi.edu/spectra/intralipid/index.html Google Scholar

23. 

V. Sorgato et al., “Wide-field absolute quantification of absorption in turbid media,” in Clinical and Translational Biophotonics, JM3A-32 (2016). https://doi.org/10.1364/CANCER.2016.JM3A.32 Google Scholar

24. 

S. L. Jacques, “Skin optics,” (1998) http://omlc.org/news/jan98/skinoptics.html Google Scholar

25. 

L. L. Randeberg et al., “Performance of diffusion theory vs monte carlo methods,” Proc. SPIE, 5862 58620O (2005). https://doi.org/10.1117/12.633028 Google Scholar

26. 

V. Sorgato et al., “Quantitative 2D maps of optical properties reconstruction: pre-clinical results on rats,” Proc. SPIE, 10412 104120P (2017). https://doi.org/10.1117/12.2286004 PSISDG 0277-786X Google Scholar

Biography

Veronica Sorgato received her PhD in biomedical optics titled “Novel multispectral imaging technique for the wide-field quantification of optical properties” from CEA (Atomic Energy and Alternative Energies Commission, Grenoble, France). She started the PhD in October 2013 after graduating with an MSc degree in biomedical engineering and a BSc degree in electronic engineering. After her PhD, she specialized in the field of medical physics.

Michel Berger works as an engineer at CEA (Atomic Energy and Alternative Energies Commission, France) and has expertise in design, characterization, and implementation of optical systems. For more than 15 years, he has worked on optical systems for clinical or preclinical applications, such as fluorescence reflectance imaging, continued diffuse optical tomography, time-resolved diffuse optical tomography, and spatial resolved diffuse reflectance spectroscopy.

Charlotte Emain worked for two years from 2010 to 2012 at Teem Photonics in the composition of lasers. Since 2012, she has worked as an optical technician in the LISA laboratory of CEA (Atomic Energy and Alternative Energies Commission, France). She has an expertise in optical instrumental design and setup assembly.

Christine Vever-Bizet, a CNRS research scientist, after finishing her PhD in biophysics worked in the Pharmaceutical Louis Pasteur University (Strasbourg, France) and in the Biophysics Laboratory of the MNHN, Paris, France. Her research was devoted to the study of tumor cell or virus photoinactivation. Since 2004, she has joined the Pierre and Marie Curie University, Paris, France, where her studies are focused on fluorescence diagnosis based on the autofluorescence modifications of the extracellular matrix proteins.

Jean-Marc Dinten is a senior scientist at the Biology and Health Division in CEA-LETI. For more than 20 years, he has been developing medical image processing and reconstruction algorithms associated with the development of innovative x-rays and optical imaging systems. Currently, he heads the Imaging Readout Systems Laboratory, which develops new optical imaging systems for health and biology applications.

Geneviève Bourg-Heckly holds a PhD in physics. After 12 years of professional experience in laser and medical imaging industries, she joined the Pierre and Marie Curie University, Paris, France, in 1993 to develop projects in biomedical optics. Her research interests are in the field of noninvasive diagnosis methods based on optical spectroscopy and imaging. Over the last years, her studies have mainly focused on the development of single-photon and multiphoton endomicroscopy.

Anne Planat-Chrétien has been project manager in the field of optical imaging applied to biology and health since 2000 at CEA-LETI, France. She is a PhD engineer in image processing modeling, simulation, and algorithms. During the last decade, she has extended her expertise to system analysis, tests monitoring, and system characterization fields. Currently, she is leading the thematic on deep reconstruction of endogenous properties of biological tissues for both clinical and preclinical applications.

© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2018/$25.00 © 2018 SPIE
Verónica Sorgato, Michel Berger, Charlotte Emain, Christine Vever-Bizet, Jean-Marc Dinten, Geneviève Bourg-Heckly, and Anne Planat-Chrétien "Validation of optical properties quantification with a dual-step technique for biological tissue analysis," Journal of Biomedical Optics 23(9), 096002 (19 September 2018). https://doi.org/10.1117/1.JBO.23.9.096002
Received: 30 August 2017; Accepted: 6 August 2018; Published: 19 September 2018
JOURNAL ARTICLE
14 PAGES


SHARE
Advertisement
Advertisement
Back to Top