1. Introduction

Polarimetry comprises a set of non-invasive and non-contact testing techniques that provide information about the optical properties of a sample [1]. Those techniques appear as very useful tools to be applied in different research fields, such as material classification [2], gas sensing [3], astronomy [4], biomedicine [5–8], and remote sensing [9] but among all, for the purpose of this work, they are an especially useful for analyzing and characterizing biological samples [10–12]. Biological tissues are turbid media that strongly scatter light and present certain characteristic degree of anisotropy due to their ultrastructure components (e.g., the dense collagen and elastin matrix with oriented fibers, etc.) [13] so these properties lead to certain polarimetric response (depolarization, retardance, etc.). In several pathologies, morphological changes occur that alter the polarization properties of tissues [14]. Therefore, polarimetry constitutes a powerful and promising optical tool for the study and classification of biological samples [11,15]. Moreover, polarimetry can be combined with other optical techniques such as multispectral imaging [15–21] or microscopy [22–24].

The Mueller-Stokes (M-S) is a widespread mathematical formalism which allows to characterize the state of polarization of light beams (Stokes vector) and the polarimetric properties of samples (Mueller matrices, M). The Mueller matrix, MM, describes the effect of a material media, such as an optical element, or a sample, upon the polarization state of light by generating an output Stokes vector, S_out, from an input light state of polarization, S_in, for linear processes [25,26]. One interesting property of the (M-S) formalism is that it is based on radiometric measurements and therefore is easy to be implemented to interpret experimental data. In addition, unlike other formalism used to describe polarization, as for instance, the Jones formalism, the M-S is valid to deal with fully polarized light, but also with depolarization content. In addition, the Mueller matrix provides the intensity and the complete set of polarization properties, diattenuation and polarizance (3-dimensional vectors D and P, respectively, in Eq. (1)), retardance and depolarization (both mixed in the 3 × 3 m submatrix) [27–30], which are sensitive to the microarchitecture of the tissue [16].

The polarimetric characteristics are not related to a particular element of the matrix M but they appear spread (or entangled) in multiple M elements in a complex way. When a physical model is not applicable, the polarimetric properties can be evaluated after separating the original M into a set of simpler matrices. The process of separation is known as matrix decomposition. In practice, there exists different strategies to perform M decomposition which can divided in three classes; product, sum and differential decompositions. The choice of the decomposition method is not trivial, and in general, it is specific to the sample being probed, and sometimes different decompositions may lead to comparable results. The best practice is to apply different decompositions and then to chose the one which performs the best in terms of computation time and reliability of the obtained results. For instance for the present work, after testing the symmetric decomposition [31] and the Lu-Chipman decomposition [32] we chose the later because it was more numerically efficient. Depolarization has been considered as a deletary effect which tends to screen or to erase the polarimetric properties (traditionally considered as the desired observables). Recently the introduction of different depolarization metrics [33–35], and in particular the depolarization spaces [15], show that there may be an interest of considering depolarization at the same level than the polarization properties. The different depolarization metrics are related to the way in which the purity of polarization of the illuminating beam is lost after interaction with the sample. In biological tissues, light scattering and the non uniform distribution of polarization properties along the probed samples, are the main causes of depolarization. Therefore depolarization metrics may unveil differences in the way light scattered or the polarization modified. Since scattering and polarization properties are linked to the ultrastructures present in the tissues, the depolarization metrics can serve to unveil differences or characteristics in the tissues, overlooked when probed with non polarized light, and which can present a true interest for biological / medical purposes or characterization. Although we focused the attention to the physical properties of the medium, other metrics such as the higher order moments of polarimetric observables recently discussed in the literature [36–39] are of potential interest and could be used in analogous procedures. Among the known depolarization metrics we chose to work with the indices of polarimetric purity (IPPs) [11,12,20]. The IPPs can be calculated from the combination of the eigenvalues (λ_i) of the covariance matrix H(M) associated to M [40]. These IPPs provide well-defined polarimetric magnitudes with easier physical interpretation and further synthetization of sample’s depolarizing mechanisms.

In the present study, we develop a predictive optical model for tissue classification based on the statistical analysis of different polarimetric metrics, directly derived from the experimental Mueller matrix of a collection of several ex-vivo chicken samples. The efficiency of the developed models, when classifying different biological tissues, is given in terms of sensitivity and specificity. We also provide the potential of the developed models for tissue characterization imaging techniques. This study could be of interest for biomedical applications such as for guided systems in surgical procedures and for early detection of some diseases.

The outline of the manuscript is as follows. In section 2, we briefly review the mathematical background of the polarimetric metrics selected for the study. In section 3, we provide the experimental methodology and the pathological description of the four different studied tissues: tendon, muscle, myotendinous junction and bones (subsection 3.1). We also detail the description of the statistical analysis conducted on polarimetric data (subsection 3.2). In section 4, we describe the construction of an optical probabilistic model suitable for the classification of biological samples, which is based on a logistical regression applied to polarimetric data, and we discuss the validity of developed models in terms of sensitivity and specificity, as well as, their suitability for imaging classification of tissues. Finally, the main conclusions of the study are provided in section 5.

2. Methods

In this section, we present the mathematical background necessary to build the proposed tissue recognition model. The model is based on a set of metrics derived from experimental Mueller matrices of biological samples. Those observables allow us to obtain a clear interpretation of some physical properties of the samples. We also provide the experimental methodology followed for the completion of the study.

2.1 Polarimetric metrics

Recalling the M expression in Eq. (1) we define the diattenuation vector D as the metric that describes the property of an optical element by which the intensity of the exiting beam depends on the polarization state of the incident beam. Analogously, the polarizance vector P describes the polarization capability of a sample, when illuminated by an unpolarized input state [41]. Such magnitudes can be computed and described directly from the combination of M coefficients by using:

Other polarimetric characteristics of a sample are also encoded in its Mueller matrix in a complex way. This information can be synthesized in a product of three pure Mueller matrices, which are functions of well-defined polarimetric magnitudes for an easier physical interpretation, by employing the so-called Lu-Chipman decomposition [32]:

Therefore, the information related to the depolarization, retardance and diattenuation properties of samples is encoded in different scalar polarimetric observables derived from those pure matrices. In particular, the total retardance R, the linear retardance δ and the optical rotation Ψ are obtained from the pure retarder matrix M_R in this way [28]:

The total retardance (R) allows the description of the global behavior of a general retarder (i.e., an elliptical retarder), as it combines the linear retardance and the optical rotation [19]. However, in the context of tissue discrimination, to include individual characteristics of the retarder constituents, as the linear retardance (δ) or the rotation applied to the linear retarder (Ψ), may be of interest. Thereby, in the present work, all them will be analyzed individually in order to consider all the contributions.

Further description of a sample can be achieved by studying its depolarizing properties. To this aim, it is worth using the indices of polarimetric purity (IPP) [11,12,20]: three invariant depolarizing indicators that contain complete and detailed information of the capability of samples to induce randomness to input polarization states [29]. The nature of the IPPs lies in the statement that any depolarizer response can be synthesized as an incoherent sum of four non-depolarizing components whose relative statistical weights are performed by combinations of IPPs. In this way, we can discriminate different types of depolarizers by only looking at the statistical weight of each pure component.

The IPP observables are obtained from the four eigenvalues λ_i (which fulfill the condition ${\lambda _0} \ge {\lambda _1} \ge {\lambda _2} \ge {\lambda _3}$) of the covariance matrix H(M), a positive semidefinite Hermitian matrix whose elements are defined through linear combinations of the Mueller matrix [20]. Particularly, the IPPs (real magnitudes P₁, P₂ and P₃) are defined as:

In addition, the depolarization index ${P_\Delta }$, which estimates the overall depolarization of a M, is built from IPPs as follows:

These indicators characterize a pure (non-depolarizing) system when ${P_1} = {P_2} = {P_3} = {P_\Delta } = 1$ is found [42]. Due to eigenvalue restrictions, IPPs are limited by the $0 \le {P_1} \le {P_2} \le {P_3} \le 1$ inequalities. The ideal depolarizer is described for ${P_1} = {P_2} = {P_3} = {P_\Delta } = 0$. Moreover, IPPs allow synthetizing the global information provided by ${P_\Delta }$, allowing to transform a 1-dimensional space into a 3-dimensional information space.

3. Experimental methodology and statistical analysis

In this section, we provide the experimental methodology and the description of the inspected tissue samples (subsection 3.1) as well as the description of the data analysis conducted on Mueller matrix-based polarimetric data (subsection 3.2). The study provided in this section lays the foundation for the subsequent construction of an optical model suitable for the classification of biological samples, that will be discussed in section 4.

3.1 Experimental methodology and sample description

In the current study, the experimental Mueller matrices of 157 biological samples are measured at three different wavelengths covering the visible range (625 nm, 530 nm and 470 nm), since different wavelengths are associated with different light penetration capability in tissues [43]. Such measurements are performed in a scattering configuration (capturing diffuse light) by using a complete Mueller imaging polarimeter (described in section 1 of Supplement 1) based on parallel-aligned liquid crystal retarders [11,12,15,25].

The inspected tissues were obtained from 25 different ex-vivo specimens of chicken thighs provided by the same commercial brand. Tissues were taken from unambiguous anatomic locations under the supervision of a pathologist: we disinserted the muscles in order to dissect skeletal muscle, tendon and myotendinous junction as well as bone and worked with the same decomposition conditions. To do so, two samples per each tissue type were obtained (a total of 23 bones, 50 muscles, 50 myotendinous junctions and 34 tendons) and subsequently frozen at $-16^{\circ}C$. Three hours before measuring, the corresponding samples were systematically defrosted in order to equalize the decomposition stage. Also, histological sections from the same regions that were studied were collected and stained using a Masson Trichromic technique, which allows a clear-cut differentiation of all connective tissue (and otherwise mesodermal as bone and muscle). The detailed process and the histological images (Figs. S2 and S3) are provided in the section 2 of Supplement 1.

An example of an image of each one of the four studied tissues is shown in Fig. 1. Soft tissues (muscle, tendon and myotendinous junction) and bone present different properties owing to their physiologic functionality and structure. In the following, we describe a physiological analysis of each soft tissue in order to well-differentiate its structure and biological components. The skeletal muscle (Fig. 1(a)) is composed by contractile myofibril chains bundled into fascicles sheathed by perimysium (connective tissue which contains type-I collagen fibers [44]) and packed within a collagen-rich epimysium layer, which protects muscle fibers from friction. By comparison, tendons (Fig. 1(b)) are composed by dense collagen (60-80% type-I collagen) fibers arranged into parallel fascicles by following the same orientation as muscle fibers [45,46] which are completely enclosed by paratenon (a thick fibroadipose layer) and fascia (collagen-rich tissue that covers both muscle and tendon) [47,48]. In turn, myotendinous junction (Fig. 1(c)) is a variable combination of both previously mentioned soft tissues, which are progressively mixed along the transition muscle-tendon: contractile myofibrils and collagen fibers are bundled into separated fascicles [49] covered by fasciae. Each fascicle contains fibers of single type (either contractile or collagenous) but fascicles of each type are intermingled in varying proportions. Finally, bones (Fig. 1(d)) are composed by a dense matrix of collagen fibers with a varying arrangement that depends on the type of bony tissue examined. The cortical bony tissue of a long bone like the one here examined is formed by parallel collagen fibers arranged following the main axis of the bone and interrupted by lacunae (spaces containing cells) at more or less regular distances. Around lacunae collagen fibers are distributed concentrically. Mineral deposits are located in the collagen fibers and, particularly, in the spaces between adjacent collagen fibers.

 

Fig. 1. Intensity image measured at 625 nm of a) muscle, b) tendon, c) myotendinous junction and d) bone tissues of a given chicken thigh. Images correspond to an area of 1.1 × 1.1 cm².

Download Full Size | PDF

Obtained experimental Mueller images correspond to a region of interest (ROI) of 512 × 512 pixels which correlates to an area of 1.1 × 1.1 cm² (the images of the experimental Mueller matrices, MMs, as well as MMs of the decomposed parameters according to section 2 are presented on Figs. S4 to S11 on Section 3 of Supplement 1). The polarimetric analysis is performed for a pure tissue region, which exclusively contains the pixels related to the desired type of tissue and from which the background has been mostly removed. By using the measured experimental Mueller images, polarimetric images corresponding to different polarimetric metrics were calculated, according to the equations provided in section 2. To suppress the possible specular reflection effects on camera and to perfectly isolate the desired type of tissue, we selected, from the original 512 × 512 pixels polarimetric images, a specific sub-image with a sub-ROI of 150 × 150 pixels (which does not include the saturated -specular reflection originated- regions) from which we compute the mean value and standard error for the different polarimetric metric sub-images. The obtained standard errors of the mean for all the metrics calculated are within a value ranging between 0.006% and 0.02%, except for the optical rotation Ψ, whose highest standard errors are ranging from 0.3% to 3.4%: those optical rotation errors could be related with the heterogeneity of the tissue fibers orientations through the whole sample. The general low values of the standard errors suggest that the computed mean values are reliable and consistent, probably due to the large number of pixels (150 × 150) measured per sample.

3.2 Statistical treatment

The object of the study is a data matrix composed of 157 tissues cases, analysed through 27 predictor variables (9 polarimetric metrics measured at three wavelength channels; 625 nm, 530 nm, and 470 nm) and one classifier variable with four categories (the type of biological tissues: muscle, tendon, myotendinous junction and bone). The 9 selected polarimetric metrics, calculated from the experimental Mueller matrix images according to section 2, are: P_Δ, P₁, P₂, P₃, P, R, D, δ and Ψ. Under this scenario, the whole data matrix corresponds to 157 x (27 + 1) items. It is important to note that there are not conditional relations when experimentally measuring the 157 cases, so they can be considered as independent observations.

In the following, two statistical studies we applied on measured polarimetric data are described: (1) an exploratory univariant statistical analysis of the polarimetric data distributions profile; and (2) a factor analysis on data matrix with principal components extraction, based on the results obtained with the exploratory analysis, to obtain a more adequate basis to construct a logistic model capable to discriminate between the four studied tissue categories. These two studies are provided in subsections 3.2.1 and 3.2.2, respectively. Note that the different statistical tests, following discussed, applied on experimental data were conducted by using the SPSS software.

3.2.1. Exploratory univariant statistical analysis on polarimetric data distributions

The exploratory univariant statistical analysis applied on the data matrix above-mentioned was performed at a global level and according to classificatory variable groups (the type of tissues). Particularly, data distribution is analyzed by studying the central tendency, the dispersion and asymmetry behavior, outlier’s identification and graphical differences between samples sorted by type of tissue characterization from Boxplot [50,51], as well as by checking normal distribution fit assumptions (assess of normality charts and contrast tests Kolmogorov-Smirnov and Shapiro-Wilk [52–55]). The conducted analysis proved that, in general, neither these data distributions nor their transformations (according to Tuckey’s ladder of powers [56]) fit the normal distribution behavior. Particularly, primary statistics show that an asymmetric data cloud cannot be fitted into a normal distribution. In this context, we observe that the 27 metrics present asymmetry: skewness values different from 0 (associated with normal distribution) indicate the mean position with respect to the median. The closest skewness value to zero is -0.111 (linear retardance, δ, distribution at 625 nm) meanwhile the highest one, 2.257, is achieved by polarizance, P, measured at 530 nm. Only 7 from the 27 measured polarimetric distributions can be associated to normal distributions, with skewness values close to zero (they present values lower than |0.840|), and they are mostly found at 625 nm. The remaining 20 polarimetric distributions present higher skewness values, being the polarimetric distributions based on the 470 nm channel those with the clearest asymmetric tendency.

Additionally, regarding the Boxplot analysis [50,51], we evaluate the total percentage of outliers (mild and extreme) to show that there is no relevant measurement bias: among the different polarimetric distributions analyzed, the outlier values range from 0% to 5.09% of the total distribution data. Because we have not done any normality assumption on measured data and because the subsequent statistical analysis is based on robust statistical techniques, we conduct a conservative approach and those detected outlier values are not removed from the study.

In Fig. 2, we present, as a representative example to illustrate the small number of outliers associated to the observed data distributions, the Boxplot for the Index of Purity P₂ distributions measured at 625 nm illumination channel and corresponding to the four studied tissues. The few observed outliers are represented in Fig. 2 by circles (mild) and stars (extremes). What is more, the asymmetric tendency of the distributions is suggested by the relative position of the median, which is not equidistant from the first and third quartiles (they are not placed at the middle of boxes). In addition, the tissue discrimination capability of the metrics is also suggested by the boxplot analysis as medians do not fit within the boxes of other tissues. For instance, Fig. 2 illustrates the potential of P₂ to discriminate muscle among the remaining tissues (see dashed red line). Similar tendency is obtained when analyzing 530 nm and 470 nm illumination wavelength (represented in Supplement 1, Figs. S12 and S13, respectively). These results can be extrapolated to the remaining M-metrics studied in this work.

 

Fig. 2. Boxplot of P₂ index for all tissues measured at 625 nm illumination channel. Red-dashed line visually represents the potential of the metric to discriminate muscle among remaining tissue types: the median of the muscle box does not fit within the other tissues’ boxes. Boxplot points out the low quantity of outliers on data distributions (which can be extrapolated to the remaining metrics): mild and extreme values are represented by circles and stars, respectively.

Download Full Size | PDF

The normality analysis above conducted on polarimetric distributions proved that subsequent data analysis must be done by following non-parametric statistics. Under this scenario, the different polarimetric distributions, sorted by type of tissue, were studied by performing a non-parametric homogeneity multiple contrast analysis: the Kruskal-Wallis test [57,58]. This test allows us to determine if there are significant statistical differences between two or more data distributions. Therefore, in our case, the outputs are connected with the tissue discriminatory potential of the studied metrics: p-values lower than 0.05 confirm, with a 95% of confidence, the capability of the polarimetric indicator to discriminate the origin of the data between different pair of tissues.

Concerning the Kruskal-Wallis test results (whole data is presented in Table S1 of Supplement 1) we confirm that P_Δ, IPPs (P₁, P₂ and P₃) and polarizance P, are the most sensitive and interesting polarimetric indicators for the studied tissue characterization, while the optical rotation Ψ is not capable to differentiate between any pair of tissues. Particularly, polarizance P seems to be the only metric capable to discriminate between myotendinous junction and bone pair. With regards to retardance indicators, they also provide some predictive potential. From this, we conclude that the main polarimetric characteristics of samples (selective absorption, retardance and depolarization) present certain signatures in the studied tissues, and therefore, the ensemble of polarimetric indicators selected seems to be adequate to be used in optical models to discriminate between the studied tissues.

The univariant analysis results above-discussed, particularly the non-normal tendency of studied polarimetric distribution, manifest the impossibility to implement our solution based on techniques such as discriminant analysis: this requires, in addition to the assumption of equality of original metrics, sorted by type of tissue, variance-covariance matrices, that these metrics follow a multivariant normal distribution behavior. As a consequence, our proposal considers an elaborated model based on a more robust technique, the logistic regression, which does not require these assumptions. At the same time, a dimension reduction of the original matrix in pursuit of non-correlated predictors (that minimize or delete the multicollinearity that generates instability on models) is aimed. Particularly, this solution considers:

The factor analysis is described in subsection 3.2.2, whereas the implementation of the statistical models and their evaluation are described in section 4.

3.2.2. Factor analysis: principal components extraction

The goal of this work is to implement multivariant probabilistic models, particularly a binary logistic regression, with the capability of determining whether a certain sample corresponds to any of the chicken tissue categories studied in this work: muscle, bone, tendon and connective tissue. Because the application of logistic regression on data does not necessarily converge when using a non-orthogonal basis, it is convenient to obtain a set of uncorrelated (orthogonal) predictors for the regression formulation which, in turn, define a reduced-dimension space. To do so, when performing factor analysis, we choose the principal components extraction as it always provides a solution in which the output factors compound a set of independent (uncorrelated) variables obtained by linear combinations of the original indicators. The extracted principal components maintain the same information of the original variables while keeping the maximum possible variance: they constitute the ideal option to implement on predictive models. It is important to remark that, because of the uncorrelation of these new components, the multicollinearity problem in a regression model is avoided.

Before conducting principal components analysis (PCA) on our particular polarimetric data matrix, we have previously computed the Bartlett’s test of Sphericity [79] and the Kaiser-Meyer-Olkin Measure of Sampling Adequacy (KMO) [80] to ensure the suitability of factor analysis for the particular studied polarimetric data. Particularly, Bartlett’s test of Sphericity is based on the null hypothesis according to which the population correlation matrix is equal to the identity matrix: in such case, the factor model is not appropriate to treat the input data [79]. On the other hand, KMO allows evaluating the degree in which each variable is predictable from the remaining ones: obtained values larger than 0.5 indicate that at least one common factor underlying the observed variables exists [80]. The output of both tests shows the suitability of applying factor analysis (PCA) to our studied data matrix. Bartlett’s test clearly rejects the null hypothesis (p-value 0.000). Likewise, Kaiser-Meyer-Olkin outputs a value of 0.650, larger than the typical cut-off, meaning that there’s a common subjacent factor for the studied variables. Because a correlation between the variables exists, the factor model is an appropriate methodology for treating the data. Hence, considering the output of both Bartlett’s and KMO tests, they all suggest that PCA may suit our particular data set when being applied.

Once the Principal Components have been calculated by using the SPSS software on the data matrix, we have computed the scree plot [71] (Fig. 3) which allows us to select a reasonable number of components to be used for the predictive models. In particular, we decided to hold all components preceding the sedimentation zone: 10 components explain more than 90% of the original metrics variance.

 

Fig. 3. Scree plot of the principal component analysis.

Download Full Size | PDF

The eigenvalues of those components, together with the data variance (in %) and the cumulative variance (in %) are shown in Table S2 of the Supplement 1. In this way, PCA reduces the 27-dimension space of M-metrics into a new 10-dimension space by identifying the original variables set underlying dimensions. The extracted factors cannot be directly measured but they enable the data structuration around a reduced number of variability axis. As the new factors are mathematical constructions, they do not necessarily have to be connected to physical properties, but they can be further interpreted by analyzing their dependence with polarimetric metrics. From a geometrical point of view, the principal components correspond to axes perpendicular to each other that better fit with the point cloud that sets the data matrix.

4. Predictive model construction: results and discussion

The selected principal components encode the polarimetric information so we can write each component as the linear combination of the 27 studied polarimetric metrics weighed by some constants provided by the so-called component score coefficient matrix, shown in Table S3 of Supplement 1. The fact that the component score coefficient matrix allows expressing the different polarimetric metrics as a linear combination of the principal components, or vice versa, is determinant because this situation defines the basis for building a probabilistic model. In a general way, the principal components (written as $\widehat C$) are obtained by multiplying the M-metrics 27×1-dimension vector, $\widehat P$, by the transposed component score matrix, $\widehat {CS}$ (a 27×10-dimension data matrix in Table S3):

To illustrate the correspondence between the principal components and the physical variables, we represent a two-dimensional space plot where the 27 indicators are scattered as a function of the two main principal components, C₁ and C₂. The obtained results are shown in Fig. 4.

 

Fig. 4. The plot showing principal components C₁ against C₂ represents the correlation coefficients between the 27 polarimetric indicators and the two first extracted principal components. The notation of the polarimetric indicators names is composed of the word that represents the measured parameter: P1, P2, P3 and PA (IPPs and P_Δ, respectively), D and P (diattenuation and polarizance), R and Delta (global and linear δ retardance) and Phi (optical rotation Ψ), followed by R, G or B (corresponding to red, green and blue measured wavelength, respectively).

Download Full Size | PDF

Note that the metrics are clustered in four main areas (highlighted by black circles) and more interestingly, those clusters mainly coincide with their physical characteristics regardless its measured wavelength: retardance related indicators (total and linear retardance, R and δ respectively), depolarization (IPPs and P_Δ), dichroism (diattenuation and polarizance, D and P, respectively) and optical rotation Ψ. By comparing the weight or influence of the two first principal components when analyzing polarimetric indicators we find that IPPs and P_Δ (scattering, depolarization) are mostly described by C₁ as they are clustered with weights around 0.5 and 1.0 meanwhile the weights of C₂ are between 0 and 0.5. When looking at how retardance indicators are grouped, we find that the presence of the C₂ principal component (between 0.5 and 1.0) is higher than the influence on the C₁ component (for the last one, its weight value is close to 0). For diattenuation and polarizance, they are both influenced by the two main principal components but they are especially represented by the C₁ component as most dispersion is observed in such a component. Finally, optical rotation principal components C₁ and C₂ values are both near 0 so they do not provide enough information to describe Ψ (the two first principal components are not able to explain optical rotation so we expect to find weights different from zero when analyzing the remaining components).

Pointing out some interesting facts: plotting principal components against each other provides us a very intuitive way to understand in which way those components are connected with actual physical information of samples (the 27 polarimetric indicators), and which one of those polarimetric indicators (or collection of indicators) carry more data variability in the original raw data matrix. Interestingly, the fact that physical variables of the same type (for instance retardance) are grouped in well-defined areas of the C₁-C₂ space (or of other combinations of 2-dimensional C_i spaces) indicates that different origin in data variability (orthogonal dimensions) is well connected with actual different physical structures present in the studied tissue samples.

The discussion of the principal component’s space characteristics (such as its discriminatory potential, sensitivity and specificity) and the predictive model construction of each type of tissue and its predictive features are described in subsections 4.1 and 4.2, respectively.

4.1 Discriminatory potential of the principal components space

In this section, we are interested in studying the the performance of 10-principal components space as tissue classifiers. With this aim, the Receiver Operating Characteristic (ROC) curve analysis [75–78] is carried out: it describes the performance of a classifier (an algorithm or a particular variable) when classifying measures into two categories by plotting the true positive rate (TPR), or sensitivity, against the false positive rate (FPR), or 1-specificity, for multiple threshold values of the classifier. For each classifier and measure, the combination of what the sample actually is (real value, positive or negative) and how it is classified by the model (prediction; positive or negative) results in four possible outcomes [76,77]. The mathematical relations connecting these four criteria define the sensitivity and specificity, which are the base of ROC curves, and are given by:

By using the SPSS software, we computed the ROC curve for all the 10-principal components and for the four different studied tissues (a specific dichotomic predictive model is constructed for each studied tissue, so four models are constructed). As an example, the ROC curve of the principal component C₁ for each biological tissue is shown in Fig. 5. Note that the different threshold are all the possible values in C₁. However, all principal components present statistical significance in ROC curves for one or more tissues, which bolsters the interest in keeping them in the process.

 

Fig. 5. ROC curve of the principal component C₁ for (a) muscle, (b) tendon, (c) myotendinous junction and (d) bone.

Download Full Size | PDF

Usually, to compare the performance of classifiers (in our case, to estimate which component provides better sensitivity-specificity values for a particular tissue), the area under the curve (AUC) is calculated, with values ranging from 0, when the variable has no predictive capability, to 1, with 100% both sensitivity and specificity. The AUC values of all principal components when discriminating each biological tissue against the remaining three are summarized in Table 1.

Table 1. AUC of 10 first principal components for all tissues.

View Table | View all tables in this article

With regards to the exposed information, we see how most of the principal components provide certain discriminative potential, with values larger than 0.5. The components with the largest AUC for a particular tissue are highlighted in Table 1. Notice that C₁ present significant discrimination values in all tissues, C₂ and C₃ mostly provide discriminatory information for tendon, C₄ and C₅ focus on muscle and tendon discrimination, C₆ on tendon and bone, C₇ on bone, C₈ keeps the information about muscle, tendon and bone and, finally, C₉ and C₁₀ provide discrimination for the myotendinous junction.

Interestingly, it has been noticed that the principal components predictive potential (given by the large area under the curve, AUC, values exposed in Table 1) is different for distinct tissues. Such behavior was expected as every component is written as a linear combination of the polarimetric indicators (rows in Table S3 in Supplement 1), and the metric values strongly depend on the composition characteristics of each biological tissue. Therefore, considering the structural differences between them, some polarimetric indicators and, consequently, some principal components, are expected to have a higher capacity to discriminate a particular tissue between the rest.

4.2 Predictive model construction

In this section, the determination of the predictive model for tissue recognition is described. The constructed model is based on a probabilistic function achieved by fitting the experimental data (array of principal components values for all chicken samples) to a sigmoid function, more specifically, to a logistic function. In other words, we perform a non-linear data model transformation into a linear one by means of a logistic regression fit of data. The logistic function with the principal components as variables and with a curve maximum value equal to 1 is written as follows,

We have chosen a logistic function because it is valid for non-parametric data and it does not require the relation between the predictors and the probability of a target outcome to be linear. Moreover, despite other fitting methods, as it is the case of those based on the ordinary least square regression (OLS), the logistic function neither requires the residuals to be normally distributed and exhibit constant variance. Particularly, we have previously verified the non-dependency of the model predictors (principal components) and their connection with the dichotomic dependent variable. Moreover, by performing the PCA, the multicollinearity between predictors have been removed so stable predictive models and convergent iteration processes are achieved. Under this scenario, four logistic regression functions are designed: one model for each studied biological type of tissue. As seen in Eq. (13) the obtained probability function depends on some of the principal components factors which, in turn, depend on the polarimetric indicator values measured from the image sample. The logistic regression fit on the experimental principal components (obtained from the experimental polarimetric indicators according to the relations provided in the component score coefficient matrix in Table S3 in Supplement 1) is conducted by using the SPSS software. By applying a stepwise regression approach (backward elimination) and using the Wald estimator throughout multiple steps, the principal components are removed until only the most significant ones remain [74].

This routine gave us the values for the weights in Eq. (13) and, accordingly, the obtained probabilistic functions for the four studied tissues are exhibit below:

To illustrate the goodness-of-fit for the four regressions, the Hosmer-Lemeshow significance [79,80] and the R² of Nagelkerke [81,82] indicators are provided in Table S4 of Supplement 1.

Afterwards, to study the efficiency of the probabilistic models in terms of sensitivity and specificity, the associated ROC curves were also computed for each one of the four models. They are provided in Fig. 6, in which their corresponding AUC are also indicated. We see how the AUC are significantly larger than those we obtained when representing the principal components (see Fig. 5 as an example), so we see how the discriminatory potential of the probabilistic functions clearly overcome the obtained with the principal components by themselves. We want to highlight that the AUC values obtained for muscle, tendon and bone tissues models, are significantly high (0.92, 0.95 and 0.89 respectively), providing the efficiency of these models for tissue discrimination. In turn, the model providing lower discriminatory potential is that associated to the myotendinous junction (AUC=0.79). This situation agrees with the fact that only C₁, C₉ and C₁₀ components present discriminatory potential for the myotendinous junction tissue (see Table 1), and the two latter, represent a very small data variation of the whole data matrix (2.350% and 1.885% of the variance (Table S2), respectively).

 

Fig. 6. ROC curve of the probabilistic model for (a) tendon, (b) muscle, (c) myotendinous junction and (d) bone.

Download Full Size | PDF

Note that models in Eqs. (14–17) provide as output a real number between 0 and 1, associated to the probability of a given tissue to belong to a particular tissue category. To build a dichotomic model determining if a tissue belongs or not to a particular category, a threshold (specific probability value) must be set in a way that values above/below such threshold are associated with Yes/No answers (i.e., is the unknown sample a bone? Yes-No; is it a tendon? Yes-No, etc.). However, each possible threshold (probability), as can be seen from ROC curves in Fig. 6, will lead to a particular pair of sensitivity and specificity values for the constructed models. In our work, the criterion used to set such threshold is by using the Youden’s Index [83] of the ROC curves in Fig. 6. According to this criterion, the optimal cut-off for our models corresponds to the farthest point of the probabilistic function in the ROC curves (blue curve in Fig. 6), from the ROC diagonal (red line in Fig. 6) along Y axis, which is calculated as the maximum value, d_max, of d = sensitivity + specificity – 1.

The computed Youden’s index [83] for each model, and its corresponding sensitivity and specificity values, are shown in Table 2. Note that the models for the muscle, tendon and bone present high values of specificity and sensitivity (values higher than 80% in all the cases). In turn, the worst obtained result is that for the myotendinous junction tissue specificity, which descends to a value of 71% (in agreement with the limited potential of the myotendinous model discussed in the previous analysis).

Table 2. AUC, Youden’s index, sensitivity and specificity of each predictive model.

View Table | View all tables in this article

Therefore, the sensitivity and specificity obtained results for the constructed predictive models (data in Table 2) quantify the efficiency of those models to categorize the studied organic tissues. What is more, the current study highlights the potential of predictive models based on polarimetric data to discriminate between animal tissues, and the methods provided could be applied in the future, for instance, for medical applications.

Following this same idea, to show the potential of the method to be applied to imaging techniques to discriminate between tissues, we encoded an algorithm that computes the probabilistic function for each pixel of the measured sample. In particular, for an arbitrary chicken tissue, the Mueller matrix image is measured at the three studied wavelengths (625 nm, 530 nm, 470 nm), from which the 27 polarimetric images are calculated. From this polarimetric image database, the probability image of the arbitrary tissue to be categorized as a particular category (muscle, bone, tendon and myotendinous junction) is calculated according to the logistic probability functions in Eqs. (14–17). The algorithm outputs which is the probability of the analyzed pixel to be recognized as a particular tissue, so the probability image is constructed, with values ranging between 0 and 1.

As an example, in Fig. 7 we provide the probability function images, corresponding to the four predictive models when analyzing an arbitrary sample of a chicken tendon. Figure 7 presents the intensity image of the tendon (M00, Fig. 7(a)), and the output images of the probabilistic models for the recognition of muscles (Fig. 7(b)), tendons (Fig. 7(c)), myotendinous junctions (Fig. 7(d)) and bones (Fig. 7(e)). Probability images are given in grayscale, in which white is associated to the maximum probability value (1) and black to the minimum (0). By comparison, the highest probabilistic values (closer to 1) are clearly obtained by the tendon predictive model (Fig. 7(c)), recognizing almost all the pixels along the tendon tissue image, and thus, perfectly classifying the sample as a tendon. In turn, the probability images for muscle (Fig. 7(b)), myotendinous junction (Fig. 7(d)), and bone (Fig. 7(e)) are based, despite some clustered high output pixels, on significantly smaller probabilistic values than those present in tendon recognition image (Fig. 7(c)). Therefore, those models do not assign the tendon sample as muscle, myotendinous junction or bone, respectively. Among these three probability images, the larger values are obtained when using the myotendinous junction model, providing that this model prediction is the less sensitive and specific one, in agreement with previous discussions on model efficiency.

 

Fig. 7. Intensity image M00 (a) and probability image outcome when applying the Muscle-model (b), Tendon-model (c), Myotendinous junction-model (d) and Bone-model (e) on an arbitrary chicken tendon sample. The gray level bars, placed to the right of the corresponding probability function images, defines the probability of the pixel to be recognized as a particular tissue, in a range between one (white) or zero (black).

Download Full Size | PDF

Analogously, we provide a second example, for the study of an arbitrary muscle sample. Results are provided in Fig. 8, for the intensity image (M00, Fig. 8(a)), and the four predictive model output images (muscle, tendon, myotendinous junction and bone recognition, Figs. 8(b-e), respectively). The obtained predictive results demonstrate the high capability of the muscle-recognition model: the highest overall values for the output probability image correspond to the muscle recognition model (Fig. 8(b)), perfectly classifying the muscle tissue. When analyzing the probability distribution for the three remaining models, a good discriminating potential is also demonstrated: the sample’s output shows low probability values for classifying the muscle as a tendon (Fig. 8(c)), a myotendinous junction (Fig. 8(d)) or a bone (Fig. 8(e)). Once again, and as expected by the previous analysis (sensitivity and specificity model values in Table 2), the worst predictive efficiency is related to the myotendinous junction model, leading to larger probability values than the tendon and bone models.

 

Fig. 8. Intensity image M00 (a) and probability function of muscle (b), tendon (c), myotendinous junction (d) and bone (e) for chicken muscle measurements. The gray level bars, placed to the right of the corresponding probability function images, defines the probability of the pixel to be recognized as a particular tissue, in a range between one (white) or zero (black).

Download Full Size | PDF

Note that the tendon and muscle study cases (results in Figs. 7 and 8) highlight the strong qualitative and quantitative discriminative potential of the models to be applied for tissue recognition imaging techniques. For the sake of completeness, the remaining tissue study cases, bone and myotendinous junction recognition test, can be found in section 6 of Supplement 1. These examples (Fig. 7, Fig. 8 and Fig. S14 and Fig. S15 of Supplement 1) point in the direction that the proposed models are suitable for tissue classification and recognition, and the corresponding sensitivity and specificity outcomes will be limited by data provided in Table 2.

5. Conclusions

Four predictive models based on the measure of diverse polarimetric metrics derived from the experimental MM of four ex-vivo chicken tissues (bone, tendon, muscle and myotendinous junction tissue) have been designed for tissue recognition and provided in this manuscript. The predictive models start by studying the distribution of the experimental data and their different statistical origins. Obtained results proved that we deal with non-parametric data with significant potential to discriminate between studied tissues. Afterwards, with the idea of constructing a robust non-parametric predictive probabilistic model, we applied a Principal Component Analysis on data, which allowed us to reduce the dimension of the data space (from 27-dimension to 10-dimensions) but dealing with maxima information as well as to obtain an orthogonal basis of parameters. By analyzing the connection between the principal components basis and the original polarimetric metrics, we realized that different variability directions of data in the principal components space can be roughly associated by different physical origins of the tissue structures (dichroism, retardance and depolarization).

Subsequently, based on the computed principal components basis, we constructed a logistic regression on data, leading to four probabilistic models (a model for each one of the four studied tissues) able to categorize if an arbitrary tissue belongs to a certain tissue category. The sensitivity and specificity values reached for each one of the four models are respectively: 85.3% and 93.5% (Tendon-model), 86.0% and 88.8% (Muscle-model), 82.6% and 80.6% (Bone-model) and 82.0% and 71.0% (Myotendinous junction-model), proving the potential of the method.

Finally, to highlight the suitability of the provided methods to be applied in biological samples imaging recognition techniques, we computed an algorithm, based on the constructed probabilistic models, that provides the probability, between 0 to 1, of a given pixel to be categorized as a particular tissue, and this is done for all the pixels of a given tissue image. In short, the procedure outputs four different probability images, corresponding to the probability of an arbitrary tissue to be muscle, bone, tendon or myotendinous junction tissue. The method provided complementary visual interpretation for tissue recognition, and a satisfactory categorization of the analyzed tissue image, according to the sensitivities and specificities stated above.

The proposed non-invasive methods discussed in this work could be applied, by conducting the required statistical data feeding, in multiple biomedical scenarios, for example, for the early diagnosis of some pathologies.