2.1.

Multi-Spectral Wide-Field Imaging MM Polarimeter

Conization specimens were imaged using the full-field multi-spectral Mueller polarimetric imaging system in a backscattering configuration described in Refs. 31, 33, 34. The main components and operational principle of the optical setup are briefly recalled in this section. Figure 1 shows the schematic layout of the multi-spectral wide-field imaging MM polarimeter used for the data acquisition. Polarization modulation of incident light from a broad band white light source is carried out by the polarization state generator (PSG), which includes a sequential assembly of a linear polarizer and two voltage-driven ferroelectric liquid crystals (FLCs). The polarization state analyzer (PSA) inserted in the detection arm of the instrument contains the same optical components as the PSG but assembled in the reverse order. Each FLC operates in transmission as a quarter wave plate with the orientation of the fast optical axis switching between 0 deg and 45 deg in the laboratory reference framework depending on the applied voltage. A light-emitting diode was used as the source of incoherent white light.

Fig. 1

(a) Schematic presentation and (b) the actual photograph of the experimental setup used to acquire MM images.

The incident light beam illuminates a sample at an incidence angle of $\sim 15\mathrm{deg}$, producing a spot size of $\sim 10\mathrm{cm}$ along the main ellipse axis within the imaging plane. Back-scattered light from a sample passes through the PSA and rotating wheel that contains 40 nm band-pass interference filters enabling the measurements from 450 to 700 nm in 50 nm increments. The image is recorded with a $600\times 800\text{pixels}$ CCD camera (Stingray F080B, Allied Vision, Germany) placed normal to the imaging plane. To measure all 16 elements of MM, one needs to perform at least 16 intensity measurements. In our polarimetric imaging system, four different polarization states of incident light beam are produced sequentially by the PSG and projected onto four different polarization configurations of the PSA, thus resulting in $4\times 4=16$ intensity images detected by the camera. The Eigenvalue Calibration Method by E. Compain³⁵ was used to calibrate the instrument, with the reference samples (two crossed polarizers and a wave plate at 30 deg) placed in a rotating wheel in front of the PSA.

The choice of the optimal polarization states for both PSG and PSA was governed by the minimization of measurement error propagation.³⁶ It was shown that it requires the points representing the polarization states of PSG and PSA on the Poincaré sphere to be the vertices of the platonic solids (a regular tetrahedral in our case).³⁷

2.2.

Histopathological Labels and Mapping

The specimens of uterine cervix from 23 patients with histologically confirmed high grade cervical intraepithelial neoplasia (CIN3 or precancer) were collected in the Kremlin-Bicetre (KB) University Hospital, France. An original method was developed to produce detailed histological mapping superimposed on Mueller polarimetric images.³¹ This method consists of the following steps: (i) all samples are formalin-fixed (according to standard protocol); (ii) shallow parallel incisions at $\sim 3\mathrm{mm}$ separation distance from each other are made by a pathologist for each specimen before polarimetric imaging; (iii) the MM images of cervical specimens are recorded with the wide-field imaging Mueller polarimeter described in previous section; (iv) all cervical specimens are cut by a pathologist into the thick slices along the pre-cut lines as shown in Fig. 2(a); (v) thin histological sections ($\sim 10\mu \mathrm{m}$) are prepared from each slide using a microtome; (vi) the histological sections are annotated by a pathologist as shown in Fig. 2(b); and (vii) the diagnosis obtained on each histological section is repositioned along the pre-cuts on the polarimetric images.

Fig. 2

Illustration of the histological labeling of a cervical specimen; (a) color photo of a sample with the positions of cuts shown by the white dashed lines (see text). The sample was placed on a glass support graded with a step of 5 mm; (b) thin hematoxylin-eosin stained histological section with the color-coded annotated segments of cervical tissue epithelium (healthy, green; CIN3, red; glandular epithelium, cyan; not annotated, black).

An example of the recorded MM images of a cervical specimen is shown in Fig. 3, (i). The experimental MM images were post-processed using the Lu-Chipman decomposition algorithm.²⁹ The corresponding map of scalar retardance with the histology analysis placed along the cut lines is shown in Fig. 3, (ii). It was found that the projection of histologically labeled lines onto the polarimetric images tended to result in small vertical and horizontal displacement errors that may affect the accuracy of polarimetric diagnosis.³¹ To mitigate this effect and to increase the number of labeled polarimetric data, we manually selected the zones that are located between the segments of two adjacent cut lines tagged with the same histology diagnosis [e.g., white-dashed zones between the adjacent gray-colored segments of cut lines in Fig. 3, (ii)] and labeled all pixels of these zones with the same histology diagnosis. This approach was used to create binary masks for both healthy and CIN3 zones for all cervical specimens (shown in Fig. 3, bottom row). Due to the lack of histologically confirmed information for low grade lesions CIN1 and CIN2, we focused our efforts on a binary classification problem (CIN3 versus healthy tissue) using a set of labeled data.

Fig. 3

Example of experimental data with histological labeling with the actual size of each image being $3\times 3\mathrm{cm}$; (I) 16 MM images of a sample. All MM images except $M_{00}$ are normalized by $M_{00}$ pixel-wise. With respect to the color bar, images on the diagonal range between 0 and 1, and the off-diagonal images range between $-0.1$ to 0.1; (II) the cut lines labeled with color-coded histological diagnosis (healthy, green; low grade CIN, yellow; high grade CIN3, red; glandular epithelium, cyan; metaplastic squamous epithelium, blue; not annotated, black) are superimposed on the scalar linear retardance image (in degrees). The semi-transparent white shapes overlaying the image illustrate how masked regions are selected between the adjacent cut lines labeled with the same histological diagnosis (see text). With respect to the color bar, this image ranges from 0 deg to 50 deg. The bottom row depicts examples of the intensity images for three specimens with the superimposed masks of histological diagnosis: (a) specimen with the healthy zones only, (b) specimen with the CIN3 zones only, and (c) specimen with both healthy and CIN3 zones. Pre-cancerous (CIN3) zones are rendered in red, and healthy zones are in green (scheme adapted from Ref. 33).

Previous research has investigated the utility of Mueller polarimetry in the context of detecting pre-cancer in cervical tissue specimens. Rehbinder et al.³¹ obtained an average value of 83% for both sensitivity and specificity. Their work involved using the scalar retardance values calculated from the MM images of 17 formalin-fixed cervical specimens as a decision variable. The work undertaken by Kupinski et al.³³ found that using the scalar retardance values or three smallest eigenvalues of the coherency matrix ($4\times 4$ complex-valued Hermitian semi-definite matrix²³) permitted high performance, achieving AUC scores of 0.93 and 0.94, respectively. The high performance of depolarization metrics based on three smallest eigenvalues of the coherency matrix in discriminating between healthy and pre-cancerous regions was a novel finding because the detection performance of total depolarization from Lu-Chipman was close to that of a random classifier (AUC = 0.5). Using both the scalar retardance and the three smallest eigenvalues enabled an average AUC of 0.95. Interestingly, it was also established that only 6 of the 16 MM elements were needed to achieve such a performance. This indicates that not all MM elements are equally relevant to the pre-cancerous classification task. For instance, Heinrich et al.³⁴ developed a general framework for the use of machine learning approaches in the classification of Mueller polarimetric data. This approach introduces a general metric, called empirical risk, allowing for training the different classifiers, choosing the polarimetric variables of interest, defining the hyperparameters of the classifiers, and comparing the performance of the classifiers. Nevertheless, it is now possible to extend beyond assessing the diagnostic performance of the scalar retardance and depolarization parameters. In particular, machine learning models can be trained directly on the labeled pixels of all images of the MM. The quality of these models’ predictions can be assessed by examining the overlap of the histologically labeled lines and the predicted classification zones.

In this paper, we utilize machine learning methods that require less overhead in their implementation. The first key difference is the manner in which model parameters, such as tree depth or number of neural network layers, are selected. The cross-validation scheme for model parameters implemented in Ref. 34 was not used due to the ability of our classifiers to achieve high performance regardless of changes to default parameter settings in this pre-screening task. Such an exhaustive search through the parameter space for optimal parameters incurs significantly greater computational cost and becomes a case of diminishing returns. In addition, we focus on the raw Mueller data instead of decomposed optical properties. This challenge is motivated by the desire to avoid pre-processing to reduce the computational effort. The approach by Heinrich et al. also used feature selection, reducing the raw data to “superpixels” to mitigate redundancy in the training process. The models that we chose to include in this work are less susceptible to data redundancy; thus this step could be omitted while retaining high performance, thereby optimizing the overall process based on the previously obtained insights.

2.3.

Artificial Intelligence and Machine Learning Techniques

The underlying idea of ML is one of automating the process of building analytical models. Without explicit instruction, systems are required to learn from patterns in data to classify them or predict other information about them. In a biomedical context, this could be learning to predict whether a region of tissue is malignant or not. Supervised ML is a sub-category of ML that uses labeled data to train models.³⁸ This means that a subset of data, or instances, is available that was labeled for the system to learn from, e.g., tissue regions labeled as malignant or not by a pathologist. After training, the resultant model is able to predict these labels when provided with unlabeled data. ML techniques can be used to solve a variety of tasks, one of which is classification wherein the system learns to map a vector into one of multiple classes.³⁹ These classes are discrete labels, such as whether a tissue region is malignant or not in a binary task. There can also be more than two classes, e.g., grades of cancer. Another common task is regression, in which the model typically predicts continuous, numerical values instead of discrete labels.⁴⁰ An example of this is quantifying biomarkers for disease states in tissues.^41–43

In the aforementioned diagnostics context, each of the masked pixels taken from the $4\times 4$ MM is treated as an instance in our dataset. Each instance has 16 features, with each feature corresponding to one of the 16 channels of the MM. Each channel contains different optical information about the pixel. Labels show the class to which a given pixel belongs: to the CIN3 class or the healthy class. In this work, a total of three ML classification techniques are used to carry out the task of predicting healthy and CIN3 masks in the cervical tissue section samples. Below, we provide the necessary background for these techniques including a decision tree (DT), multi-layer perceptron (MLP), and one-dimensional (1D) convolutional neural network (CNN).

DTs date back to the 1960s, when Morgan and Sonquist published their work on the first regression tree algorithm, AID.⁴⁴ Following this, in 1972 Messenger and Mandell published the first classification tree algorithm, THAID.⁴⁵ Breiman et al. combined the strengths of AID and THAID along with several extensions in their 1984 publication on classification and regression trees (CART), one of the most widely recognized DT algorithms to date.⁴⁶ They are typically visualized as trees, with a simple example shown in Fig. 4(a). DTs split the dataset containing the pixels, through a series of choices, into the CIN3 and healthy classes.

Fig. 4

Schematic presentation of the three classification models used in this study: (a) DT, (b) MLP, and (c) 1D CNN.

The trees consist of nodes and branches, such that the root node and internal nodes represent the decisions, the branches (denoted by black arrows) are the outcomes of these decisions, and the leaf nodes contain the class labels.⁴⁷ A decision could be whether a given pixel in a particular channel exceeds a threshold value, e.g., 0.5. Instances that have a value below 0.5 are sent down the left branch to one sub-node, and those above 0.5 will be sent down the right branch to the other sub-node. It is also possible to have non-binary decisions, such that a node maps into more than two sub-nodes. The algorithm chooses to include decisions in the tree that split the feature subspace such that instances with the same class labels are grouped together.⁴⁸ In the context of this paper, this means that the pixels are recursively partitioned into increasingly homogeneous sub-nodes, by grouping CIN3 instances together and healthy instances together as much as possible.

To measure the homogeneity of the class labels, this typically involves minimizing entropy or the Gini index.⁴⁹ These measures essentially state how impure or heterogeneous a collection of instances is with respect to each class label $j$ in a given node $i$. For entropy this involves computing the Shannon entropy of the $n$ classes:

In this work, $n=2$, which are the two class labels, CIN and healthy. Shannon entropy treats the class frequencies of the instances in each partition as a probability $p_j$. Computing the Gini index involves subtracting the sum of the squared probabilities (or proportions) of each class from 1:

MLP is another technique explored in this work. Perceptrons were initially developed as hardware by Frank Rosenblatt⁵⁰ before their algorithmic implementation. They are linear classifiers, meaning they classify inputs by forming a weighted, linear combination yielding an output value or label. Perceptrons are the very basis of MLPs, which are comprised of many perceptrons. Further still, MLPs can be considered the precursor to larger neural networks and fall within an area of AI that aims to solve complex computational tasks by imitating human brain function. The perceptron (otherwise often referred to as a neuron) is the fundamental unit of an MLP or neural network and is shown as a circle in Fig. 4(b). In essence, neurons are computational units that manipulate a weighted input signal to produce an output signal using an activation function. The weights are analogous to the coefficients in a regression model. Each neuron also has a bias, which is similar to the intercept in regression. Activation functions are a function of the weighted inputs and the bias. In a binary classification task, a very simple implementation of this can be a step function in which, if the weighted input exceeds a threshold such as 0.5, the neuron will output a value of 1.0 and otherwise 0.0, with each corresponding to one of the two possible classes. Various other non-linear activation functions can represent more complex relationships between the inputs and output,⁵¹ e.g., “ReLU” or “LeakyRelu.” In terms of their overall architecture, MLPs are made up of an input layer, one or more hidden layers, and an output layer, each containing neurons. The key algorithm underpinning the success of MLPs is back-propagation, which comprises a forward pass and a backward pass.⁵² In the forward pass, the inputs are fed into the MLP with weights initialized to some set of starting values. Each hidden and output neuron operates on its inputs by multiplying them by the initial weights, summing the result, and passing this through a non-linear activation function to produce the predicted outputs. The difference between the true output and the predicted output of the output layer is used when calculating a loss function. During the backward pass, the resultant loss is used to adjust the weights such that the predicted output of the MLP is closer to the true output. When the activation function is differentiable, this update is typically carried out using an optimization algorithm such as gradient descent or “Adam,” which changes the weights in proportion to the negative of the derivative of the error term.⁵³ In simple terms, the purpose of this optimization algorithm is to find the values for the weights that make the loss function as small as possible.

CNNs are a type of artificial neural network that use convolutional layers to extract local features from the data. Similarly to MLPs, CNNs draw their inspiration from the human brain. The structure or connectivity of a CNN is akin to that of the visual cortex.⁵⁴ Each neuron responds to stimuli in one portion of the visual field at a time, known as the receptive field.⁵⁵ A filter, equal in size to this receptive field, is element-wise multiplied with the pixels that it overlaps. This operation is known as a convolution.

Convolutional layers containing this type of processing are able to capture spatial and temporal dependencies between features through the use of filters. As shown in Fig. 4(c), a convolution layer typically uses a plurality of filters and produces multiple feature maps, each with a reduced number of features. With increasing depth, the CNN refines the number of features to those that are highly relevant to the classification task. Typically, early layers in a CNN are dedicated to extracting low-level features such as colors and edges, whereas later layers are responsible for learning higher level, holistic features of an image. Reduction of the size of the feature maps (down-sampling) is obtained either by a suitable stride of the convolution or by means of a pooling layer. Unlike convolutional layers, which contain learnable weights, pooling layers simply apply a given function. At the output of the CNN, a flattening layer can be used to transform a feature map into a 1D vector. This vector can then be classified using a linear weighting and an activation function, similar to that described in MLPs.

3.1.

Data Collection, Pre-Processing, and Quality Tests

MM images of the cervical specimens were acquired at several wavelengths (450, 550, and 600 nm, image resolution $600\times 800\text{pixels}$). Each of the 16 channels contains different polarization information. We did not apply any algorithms of the non-linear data compression (i. e., the decomposition algorithms) and used the recorded MM data to construct two broad datasets. The first one uses the four diagonal elements of MM, and another one uses all 16 values. It is known that, in general, bulk biological tissues demonstrate strong depolarization, moderate linear birefringence, and negligible linear diattenuation.²⁰ In the absence of the anisotropy of the optical refractive index, the depolarization properties of tissue can be estimated using the values of the last three diagonal elements of its MM.⁶⁰ Otherwise, these values can be considered to be a fused metric for both effects, namely, depolarization and linear retardance with the element $M_{00}$ representing the total reflectivity of tissue. The decision to use two different formats of the dataset is dictated by the necessity to ascertain whether the contribution to performance of the off-diagonal elements is sufficient to justify their inclusion. It is speculated that any additional accuracy that they can provide does not offset the computational cost. This is a method of feature selection, wherein features that are not relevant to the classification task or are redundant in that they do not add additional predictive performance from other features are discarded. Therefore, it improves the computation time by minimizing the model complexity.⁶¹

Using the histopathology masks described in Sec. 2.2, we extracted the pixels from the MM images that are known to belong to either CIN3 or healthy zones and labeled these pixels accordingly. We did not consider the low grade CIN, glandular epithelium, and metaplastic squamous epithelium because there was an insufficient number of pixels labeled with each of these histology conditions across the images. As a first step, we restrict ourselves to binary classification, focusing on the correct detection of CIN3 zones, because such lesions require surgical intervention, whereas for other pathological conditions, watchful waiting is recommended.⁷

The unlabeled pixels were not included in the training process. This extraction yielded 57,925 samples of CIN3 pixels and 133,910 healthy pixels, across all cervical tissue samples. Each instance has either four features or 16 features depending on which of the broader aforementioned datasets is used. The ratio between CIN3 and healthy instances translates to roughly 30% of the data belonging to CIN3 (class 1) and 70% belonging to healthy (class 2). Because there are high numbers of instances in both classes, the class imbalance problem is avoided, and there is minimal risk of model bias.

3.2.

Data Quality Assessment

The extracted pixels belonging to each of the two classes then underwent a series of tests. The first test explores whether the samples belonging to each of the two classes can reliably be considered to come from the same distributions. This can provide assurance that there is minimal variability in the samples that is unrelated to the classification task itself, i.e., confounding factors. To determine whether the samples came from the same distribution, a Kruskal–Willis test is applied to the pixels labeled “CIN3” and the pixels labeled “healthy” independently. The Kruskal–Willis test reveals that the instances belonging to each class do indeed follow the same distributions as each other, respectively. Their distributions in terms of mean and standard deviation are plotted in Fig. 6. It is worth noting the scale of the individual elements differences as the diagonal elements vary across a larger range than the non-diagonal elements.

Fig. 6

Mean value and standard deviation of the distributions of labeled pixels in the images of 16 elements of MM for all 23 specimens measured at 550 nm. Standard deviations are denoted by black bars, and the horizontal axis corresponds to the sample ID number.

We do not expect the MM elements to have normal distributions; to validate this and build highly tailored classifiers, normality tests were undertaken. Due to the large size of the dataset used in this work, the Kolmogorov–Smirnov test is used in favor of the more conventional, Shapiro–Wilk test for normality, as the latter is not well-suited for sample sizes containing more than 2000 instances.⁶² It was established that all samples are not normally distributed according to the Kolmogorov–Smirnov test, which concurs with our expectations.

3.3.

Testing and Training Sets

Two approaches for dividing the instances between training and test sets were considered in this work. The first approach involves a uniform, stratified, shuffled split in which 10% of instances are assigned to the test set and the remaining 90% form the training set. This approach is conventional; however, it is limited when entire samples are predicted. For instance, when the model learns from pixels belonging to all samples by nature of a shuffled split, the measured performance becomes an over-estimate of true performance. To address this, a second, less commonly used splitting approach is considered. This is known as leave-one-out cross-validation, wherein the model is trained a total of 24 times, with each sample taking a turn as the test set and the remainder constituting the training set. This will lead to a poorer but more accurate performance than the previous approach as it does not predict on samples that were used to train the model.

3.4.

Implementation and Evaluation of AIML Methods

The purpose of this work is to obtain an ML classifier with high sensitivity and specificity in distinguishing pre-cancerous and healthy regions of cervical tissue sections. To accomplish this task, three different ML models were trained and their performance evaluated. These ML approaches include a DT, multilayer perceptron (MLP), and 1D CNN. The implementation of these approaches involved the use of Python libraries, such as Scikit-learn, Keras, and Tensorflow. The optimal parameter settings for each ML method were obtained through empirical search and are described below.

Our DT implementation uses CART and splitting occurs according to the Gini index [see Eq. (1)], as opposed to Shannon entropy. Although both have comparable performance,⁶³ the former method is less computationally expensive than the latter due to the logarithm calculation [Eq. (2)]. A maximum tree depth of 10 is used to limit the tree depth. There are various motivations for this. One includes that it reduces model complexity, such that the computational time is reasonable. Allowing a model to contain many layers, such that it is highly complex, can result in diminishing returns whereby the performance does not increase significantly but the training time is very high. Furthermore, a phenomenon known as over-fitting can occur. This arises when the model becomes too specific to the training data, such that it achieves very high performance on training data, but cannot generalize well to test data and consequently cannot achieve high test performance. Finally, models that have limited depth are more comprehensible than those without restrictions, by nature of there being many fewer nodes or decisions.

The architecture of our MLP consists of three hidden layers, each with 32 hidden nodes. Following a similar tangent that motivated the design of the DT used in this work, the choice of this architecture for the MLP primarily stems from a desire to balance performance and complexity. If the number of hidden nodes and layers is insufficient, the MLP will risk under-fitting the data due to its lack of complexity, such that it is unable to map the relationships between the features, or channels, and the class labels. If this number is instead too high, the increase in complexity will not necessarily correspond to a proportionate increase in performance and may instead be deleterious, resulting in a worse performance due to over-fitting behavior. Thus, it is a matter of determining the best balance of hidden layers and nodes. We used the rectified linear unit (ReLU) as the activation function,⁶⁴ for it is widely recognized that non-linear functions are able to capture complex relationships that linear activation functions can fail to describe. The Adam optimization algorithm is used for changing the weights according to the cross entropy loss function. It has been demonstrated that Adam optimization is able to outperform conventional optimization methods, e.g., gradient descent,⁶⁵ and cross-entropy is used as the loss function due to its known utility in classification tasks. The training process is allowed to run for up to 100 iterations or weight updates. Early stopping is implemented, however, such that if the performance does not improve by 0.1% between iterations, the training process is terminated to help minimize the risk of over-fitting.

Our architecture of the 1D CNN encompasses one hidden layer containing 16 filters (corresponding to the number of channels), a kernel size of 4, and similar padding. In addition to this, an ReLU activation function, an Adam Optimizer, and a cross-entropy loss function were employed. The training process was allowed to run for 100 epochs, with no early stopping criteria. The application of a CNN in this manner is recognized as somewhat unconventional as typically a CNN operates on a region of pixels instead of one pixel at a time. This means of using a CNN was motivated by the need to make fair comparisons between the different machine learning models and examine whether the 16 channels for a given pixel contain enough information independent of other pixels to be able to classify it.

To evaluate each of the three ML models in their ability to predict pre-cancerous and healthy masks, several evaluation metrics were recorded and analyzed. These metrics include area under the curve (AUC) or accuracy, sensitivity, specificity, and training time. AUC is used for 90:10 train–test split, and accuracy is used for leave-one-out cross-validation. AUC is not applicable for the latter due to the nature of the splitting approach, wherein only three samples have both positive and negative labels, upon which the metric relies. Such an evaluation process pertains to the difference between the true mask values for each pixel and those mask values predicted by the ML models. Two averaging approaches are used. In one, 90% of the data is used for training with the remaining 10% withheld for testing or evaluation purposes. The second approach involves using 23 samples for training, with one being withheld for evaluation, in a leave-one-out cross-validation process. In both approaches, the results were averaged over many runs, with the former using the average of 30 different random seeds and the latter corresponding to the average of a different sample being withheld as the test set each time. An additional evaluation procedure was also explored, whereupon entire samples are predicted. All models were trained using a 2.8 GHz quad-core Intel Core i7 processor.