A Multi-View Semi-supervised learning method for knee joint cartilage segmentation combining multiple feature descriptors and image modalities

ABSTRACT

Multi-atlas based segmentation techniques constitute an effective approach in the automatic segmentation of medical images. Existing methods usually rely on single spectral descriptors extracted from a specific imaging modality. In this paper, we propose the Multi-View Knee Cartilage Segmentation (MV-KCS) approach, for segmenting the knee joint articular cartilage from MR images. Operating under the Semi-supervised learning framework, MV-KCS leverages spectral content from multiple feature spaces by constructing sparse graphs for each view individually, and aggregating them via optimisation to obtain a common data graph. In We consider two multi-view scenarios: in the former case views correspond to multiple feature descriptors, while on the latter, the views correspond to multiple image modalities. We propose two effective labelling schemes, implementing label propagation from the atlas library to the target image. The proposed methodology is applied to the publicly available Osteoarthritis Initiative repository. We devise a comprehensive experimental design to validate different test cases, comparing single-feature vs multi-features, multi-features vs feature stacking and multi-features vs multi-modalities. Comparative results and statistical analysis reveal that the proposed MV-KCS provides enhanced performance ($\mathit{DSC}=92.56\mathrm{\%}\left(\mathit{FC}\right),89.91\mathrm{\%}\left(\mathit{TC}\right)$), outperforming a series of patch-based approaches, six recent state-of-the-art deep supervised models and three deep semi-supervised ones, in terms of both classification and volumetric measures.

KEYWORDS:

• Osteoarthritis (OA)
• magnetic resonance imaging (MRI) segmentation
• multi-view learning
• multi-atlas patch-based
• semi-supervised learning (SSL)

1. Introduction

Osteoarthritis (OA) is one of the most widespread joint diseases worldwide, affecting as much as 10% of men and 18% of women over the age of 60 (Felipe and McCombie 2002). It is a complex condition whose onset and progression depends on a multitude of socioeconomic, genetic, biomechanical and other systemic factors (Felson 2000; Glyn-Jones et al. 2015). The leading causes of the condition are thought to be 1) bone misalignments caused by either congenital or pathogenic factors, 2) excess body weight, 3) loss of strength and muscle mass and 4) nerve damage of peripheral nerves. It primarily affects the weight-bearing joints of the human body, most commonly the knee joint, causing pain and varying degrees of disability in the patients. The pain and physical discomfort induced by the condition significantly affects the patients’ social functioning and mental health, further diminishing their quality of life. Recent rising trends in life expectancy and ageing populations, especially in the developed countries, are expected to place osteoarthritis among the leading causes of disability in the coming years. In addition to the severe negative effect in the everyday lives of the affected populations, this has the capacity to place a significant burden in national health systems. To that end, research towards a better understanding and treatment of OA is a crucial issue.

Magnetic Resonance Imaging (MRI) provides an indispensable tool in the evaluation of cartilage volume and thickness, thus allowing a robust qualitative and quantitative analysis of the cartilage morphology. However, manual delineations performed by expert radiologists and clinicians, are time-consuming and also susceptible to high inter- and intra-rater variability, highlighting the need for efficient and robust automatic methods of achieving reliable segmentation.

1.1. Related work

In the past decade, considerable research has been conducted towards fully automating the process of cartilage segmentation from MRI images. The thin structure of cartilage, coupled with the inter-slice variability of femoral and tibial cartilage shape, pose significant challenges in achieving the above goal. Most methods suggested to address these issues can be broadly distinguished into five groups, including Region-growing methods, Statistical Shape models, Graph-based methods and finally, classical machine learning & deep learning ones. A thorough review of the existing knee cartilage segmentation literature can be found in (Ebrahimkhani et al. 2020).

1.1.1. Region-based methods

Region growing algorithms operate on a predetermined set of seed points, expanding the initial regions by incorporating neighbouring voxels under a similarity measure. In most applications (Pakin et al. 2002; Öztürk and Albayrak 2016) they are utilised as an initial step to acquire a pre-segmentation mask, progressively refined until the desired result is obtained. Region growing methods are fairly easy to implement, but computationally expensive. Moreover, the required user input in the initial stages of segmentation hinders the widespread application of such methods in a fully automated setting.

1.1.2. Statistical shape models

Statistical Shape Models (SSM) and Active Appearance Models (AAM) enjoy a wide use in knee joint segmentation applications. Shape is a primary feature of rigid anatomical structures and can be successfully utilised as an initial guiding tool for the subsequent complete delineation of the corresponding structure. Finally, SSMs can also feature as shape regularisers at the end of a segmentation pipeline, performing post-processing corrections (Ambellan et al. 2019). While adequately successful in applications involving a small number of subjects, SSMs tolerate over-restrictive shape variation, often requiring a large number of landmarks to deal with larger datasets.

1.1.3. Graph-based methods

Graph-based methods treat segmentation as an energy cost function optimisation problem, where an initial global graph is split into multiple subgraphs under certain constraints, each made to correspond to an object of interest in the image (Yin et al. 2010). Overall, graph-based methods have strong theoretical foundations and achieve desirable results, but similar to the Region-Growing ones require a user’s input.

1.1.4. Classical machine learning methods

This family of methods treats knee joint segmentation as a supervised classification task, estimating the label of each voxel from a set of features extracted from the image. The large number of possible combinations of features and classifiers allows for great variety of applications (Folkesson et al. 2007; Zhang et al. 2013). However, these types of models typically suffer from low generalisability, since the extracted features are usually tailored to a specific dataset.

1.1.5. Deep learning-based methods

In recent years, deep learning methods, primarily involving Convolutional Neural Networks (CNNs) are steadily gaining popularity in medical image applications, in part due to the fact that CNNs are capable of learning an appropriate set of features automatically. In (Prasoon et al. 2013), a deep model is presented, composed of three 2D CNNS, each one corresponding to the three orthogonal views of an MRI (sagittal, planar, coronal). A 2D U-net model is reported in (Norman et al. 2018) for the segmentation of cartilage and meniscus, while a variant of the SegNet (Badrinarayanan et al. 2017) with an application to bone and cartilage segmentation is presented in (Liu et al. 2018). A combination of 2D and 3D CNNs, coupled with an SSM regularisation step is reported in (Ambellan et al. 2019), while applications featuring fully 3D networks constitute only a recent development (Zhou et al. 2018; Dai et al. 2021; Peng et al. 2022), due to the vast computational load incurred in such cases. Although deep learning methods achieve overall attractive segmentation results, the lack of large-scale annotated medical image repositories might lead to overfitting, with the researchers having to rely on fine tuning of pre-trained CNNs, or artificially augmenting the existing available datasets.

1.2. Semi-supervised deep learning-based methods

Semi-supervised deep-learning methods constitute a promising alternative to standard deep-learning models in the face of sparse annotated data, as is the usual case in medical image segmentation applications, while at the same time showcasing substantially faster computation times. The authors in (Zhang et al. 2017) propose a deep adversarial network (DAN) comprising two sub-modules: a segmentation network that provides label maps and an evaluation network to assess the quality of said segmentations. In (Yu et al. 2019), the uncertainty aware mechanism is utilised, employing a student model and teacher model, whereby the student model attempts to minimise the segmentation loss by utilising labels supplied by the teacher, while the teacher continuously estimates an uncertainty map to bias the student towards harnessing information from the more reliable annotated targets. Finally, in (Luo et al. 2021) the authors design a novel uncertainty-aware mechanism, expanding the previous work of (Tarvainen and Valpola 2017), by designing a dual-task deep network that simultaneously learns a pixel-wise segmentation map and a level-set representation of the target.

1.3. Multi-atlas patch-based methods

Multi-atlas patch-based methods segment a target MRI $\mathcal{T}$ by utilizing a collected Atlas Library $\mathcal{A}=\mathcal{L}\mathcal{B}(\mathcal{T})=\{{\mathcal{A}}_i,{\mathcal{L}}_i{\}}_{\mathit{i}=1}^{{\mathit{n}}_{\mathcal{A}}}$, consisting of $n_{\mathcal{A}}$ subject MRIs $({\mathcal{A}}_i)$ and their corresponding labeled masks $({\mathcal{L}}_i)$. Voxels in the target image $\mathcal{T}$ are automatically assigned labels by harnessing the available information from library $\mathcal{A}$. A necessary condition for the successful application of this process is that the atlas images in $\mathcal{A}$ and the target image $\mathcal{T}$, all share a common coordinate space. This is commonly achieved by initially registering all atlases ${\mathcal{A}}_i\in \mathcal{A}$ to the target $\mathcal{T}$, usually via an affine transformation or a deformable one (Hajnal et al. 2001).

These methods typically proceed as follows: For each voxel $\mathbf{x}\in \mathcal{T}$, a search volume $\mathit{N}(\mathbf{x})$ of size $N_s=(\mathit{n}\times \mathit{n}\times \mathit{n})$ with $\mathbf{x}$ at its centre is defined, and each voxel $\mathbf{y}\in \mathit{N}(\mathbf{x})$ yields a patch $\mathit{P}$ of size $P_s=(\mathit{p}\times \mathit{p}\times \mathit{p}),\mathit{p}<\mathit{n}$. Each of these patches gives rise to a corresponding patch vector $\mathbf{p}(\mathbf{y})\in {\mathbb{R}}^{P_s}$, comprising the intensity values of all voxels in $\mathit{P}$. The collection of all patch vectors extracted from the search volumes within each atlas $\{{\mathcal{A}}_i{\}}_{\mathit{i}=1}^{{\mathit{n}}_{\mathcal{A}}}$, forms the patch library corresponding to the target voxel $\mathbf{x}$, ${\mathbf{P}}_{\mathcal{L}}=\{{\mathbf{p}}_{{\mathcal{A}}_i}(\mathbf{y}),\mathrm{\forall }\mathbf{y}\in \mathit{N}(\mathbf{x}),\mathit{i}=(1,\dots ,n_{\mathcal{A}})\}$, ${\mathbf{P}}_{\mathcal{L}}\in {\mathbb{R}}^{P_s\times N_s{n}_{\mathcal{A}}}$.

Under the assumption that spectral and label similarity across patches are encoded by the same underlying function, the goal is to express ${\mathbf{p}}_{\mathcal{T}}(\mathbf{x})$ as a linear decomposition of the patches in ${\mathbf{P}}_{\mathcal{L}}$, via coefficients $\mathbf{w}\in {\mathbb{R}}^{N_s{n}_{\mathcal{A}}}$ reflecting that similarity, then reconstruct the label of ${\mathbf{p}}_{\mathcal{T}}(\mathbf{x})$ as a combination of the corresponding atlas patches, utilizing those same coefficients. Typical examples of this approach are found in (Rousseau et al. 2011; Zhang et al. 2012), where $(\mathbf{w})$ is computed via Non-local Means (NLM) filtering (Buades et al. 2011) and Sparse Coding SC optimisation (Zhang et al. 2015), correspondingly.

Patch-based methods overcome many of the challenges encountered by previous approaches, while simultaneously achieving state-of-the-art segmentation performance, as long as a suitably large and variable atlas library ${\mathbf{P}}_{\mathcal{L}}$ is chosen. However, the voxel-wise manner in which the segmentation process is carried out incurs large computational costs, since the full segmentation of a target image entails the exhaustive extraction of a large number of atlas patches, for each single target voxel.

1.4. Multi-view semi-supervised learning

Multi-view data are widespread among real world applications. In most cases, multiple views on the same dataset result either from a multitude of measuring methods, or from expressing a single view with distinct sets of features. Multi-view methods aim to learn a function for each view separately and then jointly optimise them in order to boost generalisation performance (Sun 2013).

As the availability of data, as well as the ease of their acquisition is steadily increasing across multiple fields, multi-view learning algorithms are becoming more prevalent in the computational intelligence community. Successful applications of multi-view methods have been reported for numerous machine learning tasks, such as transfer learning, dimensionality reduction, clustering, semi-supervised learning and multi-task learning. A more comprehensive list of recent advances and the new challenges arising, can be found in (Zhao et al. 2017).

In this study, we exclusively focus on applying multi-view learning (MVL) under the semi-supervised (SSL) framework. Adopting the manifold assumption, that the available data lie in a lower dimensional manifold embedded in a higher dimensional feature space, SSL aims to leverage the spectral information of few labelled data to predict the labels of the unlabelled ones (Belkin et al. 2006). Since direct access to the underlying manifold structure is impossible, a proxy is often used through the construction of a weighted spectral affinity graph $\mathit{G}(\mathit{V},\mathit{E},\mathit{W})$, where the vertices $\mathit{V}$ correspond to the data points and edges $\mathit{E}$ connect those vertices that display some similarity under a given measure, encoded in an adjacency matrix $\mathit{W}$ (Belkin and Niyogi 2005). SSL is particularly useful in applications where there is an abundance of unlabelled data, but acquisition of the corresponding labels is a costly and tedious process. Medical image segmentation falls exactly in this category, whereby the manual annotation of an MRI by experts may induce considerable workload.

The adoption of multi-view learning is a natural extension of the SSL framework, with the potential to greatly boost the latter’s performance. Given $\mathcal{V}$ distinct views of the same dataset, Multi-view Semi-supervised learning (MV-SSL) assumes the existence of $\mathcal{V}$ distinct manifolds, and constructs $\mathcal{V}$ different weighted graphs $\{\mathit{G}(\mathit{V},\mathit{E},\mathit{W}{)}^{(\mathit{v})}{\}}_{v=1}^{\mathcal{V}}$ corresponding to each one. Various combination schemes have been proposed for the fusion of those multiple graphs into a single entity (Sun 2011; Karasuyama and Mamitsuka 2013). The strength of MV-SSL stems from its ensemble-like rationale. Multiple views individually encode different spectral information of the data. Most often, some aspects of the data variability ignored by a particular view may be captured by another one and vice-versa. Hence, their MV integration boosts significantly the robustness of the overall model.

1.5. Outline of proposed method

Recently, we have proposed a multi-atlas patch-based approach for automatic knee cartilage segmentation utilising the SSL principles (Chadoulos et al. 2022). The MV-KCS method proposed herein integrates effectively the frameworks of MV and SSL in a unifying scheme, thus leveraging the strengths of both to enhance the classification accuracy. In that respect, MV-KCS is an extension of our previous work which can be regarded as a single-view approach. Figure 1 provides a succinct illustration of the algorithmic procedure. A summary of the salient features and contributions of our proposed framework are listed below [noitemsep]

• Multi-view Semi-supervised learning (MV-SSL): Under the multi-atlas setting, we classify target voxels by utilising both labelled data from the atlas library $\{\mathcal{A},\mathcal{L}\}$ and unlabeled ones from the target image $\mathcal{T}$. Moreover, we further embed the SSL framework within the broader scope of MVL, thus enhancing the robustness of our method via leveraging multiple distinct views, each one capturing a different aspect of the data.

• Sparse Coding (SC): The multiple weighted affinity graphs corresponding to each distinct view are constructed via a Sparse Coding optimisation process, whereby each voxel is represented as a linear combination of its sparse neighbours. The incorporation of sparsity aims at reducing noise and enhancing the robustness of the proposed method.

• Voxels’ labeling: Adapting to the MVL framework, we propose two multi-view graph-based labelling schemes, namely, Regional Label Propagation (MV-RegLP) and Hybrid Label Propagation (MV-HyLP) for the classification of the target voxels. The former considers the voxel descriptors at the regional level, while the latter one integrates further a finer level of description by leveraging the notion of Label Map Estimates.

• Out-of-sample labeling: An iterative sampling process successively generates out-of-sample batches of target voxels ${\mathbf{X}}_o$ through progressive densification of a 3D mesh. In addition, local data are generated by sampling the atlas library in spatially correspondent locations. The resulting labeling scheme leverages information by considering voxel similarity along two axes, namely, global – local and spectal – spatial.

• Experimental setup: The proposed MV-KCS considers two cases of multiple views, namely, different feature sets (descriptors), and different MRI acquisition protocols (sequences). Furthermore, we investigate several test cases, such as single-feature vs multi-features, multi-features vs feature stacking and multi-features vs multi-modalities, with the aim to demonstrate the efficacy of the multi-view setting.

Figure 1. General flowchart of the proposed method.

The remainder of this paper is organised as follows: Section 2 focuses on the preliminary steps undertaken with regards to image preprocessing, registration and atlas selection. In Section 8.3, we present the two sets of views examined in this study, namely, the feature sets extracted from popular image descriptors and the different MRI sequences. Section 3 describes the adopted multi-view aggregation framework while Sections 4, 6 and 7 detail the proposed MV-KCS algorithm. Details relevant to the data and model hyperparameters, as well as the experimental layout are presented in Section 8, while Section 7 deals with the presentation and discussion of the experimental results. Finally, Section 8 concludes this study with an overview on the key findings and a number of possible future extensions.

2. Background

We treat the knee cartilage segmentation as a multi-class classification ($\mathit{c}=5$): namely $\{\mathit{Background}:0,\mathit{FemoralBone}:1,:2,\mathit{TibialBone}:3,\mathit{TibialCartilage}:4\}$. Our primary focus lies on the successful segmentation of the cartilage structure (femoral & tibial) while a secondary focus is placed on the bone structure.

2.1. Image preprocessing

Knee cartilage segmentation is a challenging task, owning mainly to the poor separation of the articular cartilage and the surrounding tissues (excluding bone), as well as the fact that the intensity profiles and texture of image regions occupied by those structures exhibit close similarities. This problem is further accentuated by the inter-subject variability inherent in magnetic resonance imaging. To ameliorate these shortcomings, we pre-process each MRI as follows:

1. Curvature flow filtering: A curvature-driven flow filter (Sethian 1999) is applied to each MRI, preserving surface boundaries between adjacent structures while simultaneously smoothing homogeneous image regions.

2. Inhomogeneity correction: N3 intensity nonuniformity bias field correction (Sled et al. 1998) is applied to each image to reduce the intra-subject variability in intensity profiles across similar structures.

3. Intensity standardization: All MRI histograms are standardised to a common template, according to the method described in (Nyul and Udupa 2000), in an effort to reduce inter-subject variability.

4. Non-local-means denoising: A final filtering step is performed to account for any leftover artefacts from the previous steps and to further reduce noise. We opted for non-local-means patch-based denoising (Buades et al. 2011), due to its robust performance and widespread use in medical image applications. As a final step, all image intensities are rescaled to $[0,100]$.

2.2. Registration & atlas selection

An essential component of any multi-atlas approach is the atlas library. Given a target image $\mathcal{T}$ to be segmented, a series of pairwise transformations is performed, registering all atlases $\{{\mathcal{A}}_i{\}}_{\mathit{i}=1}^{{\mathit{n}}_{\mathcal{A}}}$ to $\mathcal{T}$. To avoid the computational cost associated with the use of deformable models, here we employ an affine transformation. This registration aligns all atlases to the target image coordinate space, covering all linear deformations such as translations, rotations, shear, scale, etc. The derived transformations are also applied to the corresponding label masks of each atlas ${\mathcal{L}}_i$, to maintain the spatial correspondence between atlas image and mask pairs. The result of the above process is the construction of the atlas library $\{{\mathcal{A}}_i^{\mathcal{T}},{\mathcal{L}}_i^{\mathcal{T}}{\}}_{\mathit{i}=1}^{{\mathit{N}}_{\mathcal{A}}}$ registered to $\mathcal{T}$.

Next, in an effort to reduce the required number of voxels to be segmented, we define a Region Of Interest (ROI), covering the entire cartilage structure and its surrounding volume. A pre-segmentation binary mask is derived via a Majority Voting filter, operating on the registered atlas masks, by discarding the bone (Femoral & Tibial) labels and regarding the corresponding cartilage ones as a single unit. Finally, a binary morphological dilation filter expands this mask, yielding the resulting cartilage ROI for the image $\mathcal{T}$, as illustrated in Figure 2. The ROI effectively defines the sampling volume for the target image $\mathcal{T}$ and its associated atlas library $\{{\mathcal{A}}_i,{\mathcal{L}}_i{\}}_{\mathit{i}=1}^{{\mathit{N}}_{\mathcal{A}}}$, disregarding all exterior voxels. The final morphological dilation step guarantees that the ROI encompasses all of cartilage voxels across the target and atlas images.

Figure 2. Each atlas ${\mathcal{A}}_i$ is affinely registered to the target image $\mathcal{T}$, and the obtained transformation is subsequently applied to the corresponding masks ${\mathcal{L}}_i$. A majority voting filter (MV) produces an initial pre-segmentation mask $({\tilde{\mathcal{L}}}^{\mathcal{T}})$ of cartilage-only labels and the Region of Interest (ROI) is finally defined as the morphological dilation of ${\tilde{\mathcal{L}}}^{\mathcal{T}}$.

Given that non-linear deformations cannot be effectively handled by affine transformations, we perform a final atlas selection step by retaining in the atlas library only a subset of the registered atlases that exhibit good spatial alignment with the target image. To this end, a measure of spatial misalignment is calculated for every pair $\{\mathcal{T},{\mathcal{A}}_i^{\mathcal{T}}\}$ using the Mean Squared Difference $(\mathit{MS}{D}_i^{(ROI)})$, where the voxels participating in the computation are exclusively those belonging in the ROI. The first $n_{\mathcal{A}}$ atlases with the least scores are included in the library, while the rest are discarded (Algorithm 1).

Algorithm 1:

Affine Registration & Atlas Selection

Data: Target $\mathcal{T}$, Atlases $\{{\mathcal{A}}_i,{\mathcal{L}}_i{\}}_{\mathit{i}=1}^{{\mathit{N}}_{\mathcal{A}}}$, $n_{\mathcal{A}}$ number of atlases to select

Result: Registered selected atlases $\{{\mathcal{A}}_i^{(\mathcal{T})},{\mathcal{L}}_i^{(\mathcal{T})}\}$, ROI

1. for $\mathit{i}\leftarrow 1$ to $n_{\mathcal{A}}$ do

2. ${\mathcal{A}}_i^{(\mathcal{T})},{\mathcal{L}}_i^{(\mathcal{T})}=$ $\mathit{Register}(\mathcal{T},{\mathcal{A}}_i$);

3. ${\tilde{\mathcal{L}}}^{(\mathcal{T})}=\mathit{MV}(\{{\mathcal{L}}_i^{(\mathcal{T})}\});$

4. $\mathit{ROI}=\mathit{BinaryMorphologicalDilation}({\tilde{\mathcal{L}}}^{(\mathcal{T})});$

5. $\mathit{MSE}=\mathit{mse}(\mathcal{T}(\mathbf{x}),{\mathcal{A}}_i^{(\mathcal{T})}(\mathbf{x})),\mathrm{\forall }\mathbf{x}\in \mathit{ROI};$

6. select first $n_{\mathcal{A}}$ atlases corresponding to $n_{\mathcal{A}}$ smallest $\mathit{MSE}$s

3. Multi-view acquisition

3.1. Patch/region description

A robust and distinctive characterisation of the voxels’ information content is essential for cartilage segmentation. In that regard, the spatial scale considered is crucial; a large spatial scale suggests that the resulting features cannot encode meaningful information, while a small one leads to noisy features. Here, we introduce two units of spatial aggregation, namely, the patch and the Regional Search Volume (SV). The patch ${\mathbf{p}}_i=\mathbf{p}({\mathbf{x}}_i)$ is defined as a $5\times 5\times 5$ volume centered around a voxel ${\mathbf{x}}_i$. Each patch ${\mathbf{p}}_i$ is characterized by a feature vector ${\mathbf{x}}_i={\mathbf{h}}_i=f_{\mathit{enc}}({\mathbf{p}}_i)\in {\mathbb{R}}^q$, with $\mathit{q}$ representing the feature dimensionality and $f_{\mathit{enc}}$ the associated encoding function of the descriptor.¹ Correspondingly, each feature vector is associated with a label vector ${\mathbf{y}}_i=[y_{\mathit{i},1},\dots ,y_{\mathit{i},\mathit{c}}{]}^T\in {\mathbb{R}}^c$. We distinguish between hard label vectors, derived from atlas images $({\mathbf{y}}_i^{\mathcal{A}})$ and soft ones, derived from the target image $({\mathbf{y}}_i^{\mathcal{T}})$.

The Regional Search Volume (SV) ${\mathcal{R}}_i$ is defined as a $15\times 15\times 15$ volume, comprising $3\times 3\times 3$ patches, with the central patch corresponding to the one associated with voxel ${\mathbf{x}}_i$. ${\mathcal{R}}_i$ may be viewed as the set of all patches contained within its spatial region, ${\mathcal{R}}_i=\{{\mathbf{p}}_{i_1},\dots ,{\mathbf{p}}_{{\mathit{i}}_{27}}\}$ where ${\mathbf{p}}_{i_1}$ corresponds to the central voxel ${\mathbf{x}}_i$ of ${\mathcal{R}}_i$. The feature descriptor associated with the whole region is formed as a weighted aggregate of all the individual feature vectors of its constituent patches, ${\mathbf{H}}_i={\sum }_j{w}_{h_j}{\mathbf{h}}_j^{(i)},{\mathbf{H}}_i\in {\mathbb{R}}^q$, where the weight associated to each patch is determined according to its city-block distance from the central one (Figure 3). Finally, the label vectors are concatenated in a row-wise manner to produce the corresponding label matrix of ${\mathcal{R}}_i$, ${\mathbf{Y}}_i=[{\mathbf{y}}_{i_1},\dots ,{\mathbf{y}}_{{\mathit{i}}_{27}}{]}^T\in {\mathbb{R}}^{27\times \mathit{c}}$.

Figure 3. Schematic representation of the feature description process for a regional SV ${\mathcal{R}}_i$, centered around a voxel ${\mathbf{x}}_i$. Each SV comprises $3\times 3\times 3$ patches, which in turn consist of $5\times 5\times 5$ voxels. For all patches ${\mathbf{p}}_{\mathbf{j}}^{\mathbf{(}\mathbf{i}\mathbf{)}}\in {\mathcal{R}}_i$, a feature vector ${\mathbf{h}}_j^{(i)}$ is calculated according to an encoding function (HOG, LBP, etc.), here denoted as $f_{\mathit{enc}}(\cdot )$, where ${\mathbf{h}}_j^{(i)}$ serves as the individual descriptor for ${\mathbf{p}}_j^{(i)}$. Finally, the whole ${\mathcal{R}}_i$ is characterized by ${\mathbf{H}}_i$ as the weighted aggregation of all its constituent descriptors.

3.2. Feature descriptors

In this section we introduce our main multi-view case, involving the use of multiple feature descriptors to characterise a voxel and/or SV. Feature descriptors are selected so that they exhibit minimum redundancy, while simultaneously offering complementary information.

3.2.1. Histograms of Oriented Gradients (HOG)

Histograms of Oriented Gradients (HOG) (Navneet and Triggs 2020) are a state-of-the-art descriptor with established performance across multiple computer vision tasks, among the medical imaging ones (Sarwinda and Bustamam 2018). Here, we present a fully 3D variant of the HOG descriptor, based on (Kläser et al. 2008). Our own implementation is modified to fit the feature extraction framework presented in Section 3.1. For a Regional SV ${\mathcal{R}}_i$, we compute the mean gradient vector of each constituent patch ${\mathbf{p}}_j^{(i)}\in {\mathcal{R}}_i$, ${\bar{\mathbf{g}}}_{\mathit{j}}^{(\mathit{i})}\in {\mathbb{R}}^3$. To avoid the unequal bin size of the traditional 2D histogram over the surface of a sphere (Rister et al. 2017), we project ${\bar{\mathbf{g}}}_{\mathit{j}}^{(\mathit{i})}$ onto the vertices of a regular polyhedron, obtaining the orientation histogram ${\mathbf{q}}_j^{(i)}=\mathbf{C}{\bar{\mathbf{g}}}_{\mathit{j}}^{(\mathit{i})}$, with $\mathbf{C}=[{\mathbf{c}}_1,\dots ,{\mathbf{c}}_q{]}^T\in {\mathbb{R}}^{\mathit{q}\times 3}$ the matrix with the polyhedron vertices’ coordinates. ${\mathbf{q}}_j^{(i)}$ is subsequently normalized and thresholded to yield the patch feature descriptor ${\mathbf{h}}_j^{(HOG,i)}\in {\mathbb{R}}^q$. The Regional SV’s descriptor ${\mathbf{H}}_i^{(HOG)}$ is finally obtained as illustrated in Figure 3. We found that $\mathit{q}=42$ provides a good balance between descriptor size and performance.

3.2.2. Local Binary Patterns (LBP)

Local Binary Patterns (LBP) (Ojala et al. 1996) enjoy a widespread use owning to their robustness and ease of implementation. Based on (Banerjee et al. 2013), we propose a 3D extension on the traditional 2D LBP by utilising the concept Spherical Harmonics, thus circumventing the difficult problem of uniform sampling on the surface of a sphere and simultaneously achieving rotational invariance. The information content of a patch ${\mathbf{p}}_i\in {\mathcal{R}}_i$ is characterized through the encoding function $\mathbf{f}=f_{\mathit{lbp}}=[(\mathit{s}(g_0-g_c),\mathit{s}(g_1-g_c),\cdots ,(g_{\mathit{P}-1}-g_c)]\in {\mathbb{R}}^P$, where $g_c$ the intensity of the central voxel, $g_i,\mathit{i}=1,\cdots ,\mathit{P}-1$ the intensities of the $\mathit{P}$ neighbors sampled on the vertices of a regular polyhedron at radius $\mathit{R}$ and $\mathit{s}(\cdot )$ the step function. Sampling the harmonic basis functions’ values at the locations of the polyhedron vertices we obtain ${\mathbf{Y}}_l^m\in {\mathbb{R}}^{P\times l^2}$ (degree $\mathit{l}$, order $\mathit{m}$) and proceed to extract the expansion coefficients of $\mathbf{f}$ via the following inner product ${\mathbf{c}}_l^m={\mathbf{f}\mathbf{Y}}_l^m$. $\mathbf{f}$ can now be reconstructed as a band-limited expansion through with $\mathit{l}$ components, via $\widetilde{{\mathbf{f}}_k}={\sum }_{m=-k}^k{c}_l^m{Y}_l^m$ and the patch descriptor ${\mathbf{h}}_j^{(i)}$ is finally obtained by taking the $\mathit{l}$-2 norm of the expansion’s frequency components, ${\mathbf{h}}_j^{(i)}=\mathit{LB}{P}_{P,R}^{ri3D}=\{{\widetilde{\parallel \mathbf{f}}}_0\parallel ,\parallel \widetilde{{\mathbf{f}}_1}\parallel ,\cdots ,\text{break}\parallel {\tilde{\mathbf{f}}}_{\mathit{l}-1}\parallel \}$. As for HOG, we choose a polyhedron with $\mathit{P}=\mathit{q}=42$ vertices, while radius is set at $\mathit{R}=2$.

3.2.3. Gray-Level Co-occurrence Matrix (GLCM) & Gray-Level Run-Length Matrix (GLRLM)

Features based on the computation of Gray-Level Co-occurrence Matrices (GLCM) (Haralick et al. 1973) and Gray-Level Run-Length Matrices (GLRLM) (Galloway 1975) are among the most commonly used for the description of image texture. Given an image with $\mathit{L}$ levels, both methods rely on the creation of a $\mathit{L}\times \mathit{L}$ matrix: GLRLM yields a matrix whose $(\mathit{i},\mathit{j})$ element denotes the number of occurrences of the intensity value $\mathit{i}$ for runs of length $\mathit{j}$ at a given direction $\mathit{\theta }$, while in GLCM, the element $(\mathit{i},\mathit{j})$ stores the number of times the pair of intensity values $\mathit{i},\mathit{j}$ is present along a given direction $\mathit{\theta }$ and at a distance $\mathit{d}$. The extension to the 3D case is possible through the inclusion of an additional direction $\mathit{\varphi }$. A typical arrangement of directions and the corresponding displacement vectors for the computation of GLCM and GLRLM is found in Table 1.

Table 1. Displacement vectors for GLCM and GLRLM for volumetric data.

Direction	Displacement Vector
$(0,\mathit{\pi }/2)$	$(\mathit{d},0,0)$
$(\mathit{\pi }/4,\mathit{\pi }/2)$	$(\mathit{d},\mathit{d},0)$
$(\mathit{\pi }/2,\mathit{\pi }/2)$	$(0,\mathit{d},0)$
$(3\mathit{\pi }/4,\mathit{\pi }/2)$	$(-\mathit{d},\mathit{d},0)$
$(\cdot ,0)$	$(0,0,\mathit{d})$
$(0,\mathit{\pi }/4)$	$(\mathit{d},0,\mathit{d})$
$(\mathit{\pi }/4,\mathit{\pi }/2-arcsin\sqrt{1/3}$	$(\mathit{d},\mathit{d},\mathit{d})$
$(\mathit{\pi }/2,\mathit{\pi }/4)$	$(0,\mathit{d},\mathit{d})$
$(3\mathit{\pi }/4,\mathit{\pi }/2-arcsin\sqrt{1/3})$	$(-\mathit{d},\mathit{d},\mathit{d})$
$(0,3\mathit{\pi }/4)$	$(\mathit{d},0,-\mathit{d})$
$(\mathit{\pi }/4,\mathit{\pi }/2-arcsin\sqrt{1/3})$	$(\mathit{d},\mathit{d},-\mathit{d})$
$(\mathit{\pi }/2,3\mathit{\pi }/4)$	$(0,\mathit{d},-\mathit{d})$
$(3\mathit{\pi }/4,arcsin\sqrt{1/3})$	$(-\mathit{d},\mathit{d},-\mathit{d})$

In our case, GLCM and GLRLM matrices are calculated for every patch ${\mathbf{p}}_j^{(i)}\in {\mathcal{R}}_i$ to be characterized. A set of features is then calculated from each matrix, yielding the corresponding feature vectors for the patch under consideration, ${\mathbf{h}}_j^{(i),GLCM}\in {\mathbb{R}}^{{\mathit{q}}_{\mathit{GLCM}}}$ and ${\mathbf{h}}_j^{(i),GLRLM}\in {\mathbb{R}}^{{\mathit{q}}_{\mathit{GLRLM}}}$.

3.2.4. Hand-Crafted Geometric Features (HCGF)

Finally, inspired by (Folkesson et al. 2007), we consider a final set of geometric features, heavily relying on the gradients of voxels’ intensities. Denoting $\mathit{I}$ the intensity of a voxel and $I_x,I_y,I_z$ the partial derivatives of $\mathit{I}$ along the $-\mathit{x},-\mathit{y},-\mathit{z}$ axes, we compute the following features for every voxel with a given patch ${\mathbf{p}}_i\in {\mathcal{R}}_i$: spatial coordinates $(\mathit{x},\mathit{y},\mathit{z})$, intensity $\mathit{I}$, first order partial derivatives $I_x,I_y,I_z$, eigenvalues and eigenvectors of the Hessian matrix $\mathit{H}=[I_{\mathit{jk}}]$ and the corresponding ones of the structure tensor $\mathit{T}=[I_j{I}_k],\mathit{for}\mathit{j},\mathit{k}\in \{\mathit{x},\mathit{y},\mathit{x}\}$. This results in a $\mathit{q}=36$ dimensional feature vector for every voxel. To obtain a single feature descriptor for the whole patch, and to avoid the curse of dimensionality, we follow the same approach as in the regional descriptors, taking the weighted sum of all feature vectors, where the weights correspond to the city-block distances from the central voxel.

The above features provide a diversified description of the image content, covering many aspects of the data variability. In particular, HOG and HCGF are well suited for capturing the shape and orientation of anatomical structures, while LBP, GLCM and GLRLM are excellent at describing their texture. Each descriptor gives rise to a different view on the available data ($\mathcal{V}=5$, Figure 4).

Figure 4. Schematic representation of multi-view feature extraction for regional SV ${\mathcal{R}}_i$. Each one of the $\mathit{j}=1,\cdots ,27$ patches is supplied with a feature vector for each descriptor(view) employed. The regional feature vector for each view is then calculated according to the process outlined in Figure 3. The overall description of ${\mathcal{R}}_i$ is the collection of all $\{{\mathbf{H}}_{\mathbf{i}}^{\mathbf{(}\mathbf{v}\mathbf{)}}{\}}_{v=1}^{\mathcal{V}}$.

3.3. MRI sequences

We also consider a secondary multi-view case, whereby multiple views are associated with different MRI acquisition protocols (sequences). The sequences used in our experiments are the Dual Echo Steady State (DESS), T2-weighted (T2) and SPoiled Gradient Recalled echo (SPGR). DESS is the most widely used acquisition protocol for musculoskeletal and joint imaging, due to the adequate delineation it offers between cartilage, bone and some of the surrounding tissue. SPGR, while not yielding a similar quality of delineation as DESS, produces images that are characterised by greater homogeneity within the same anatomical region. Finally, T2 sequences are also widely applicable in musculoskeletal imaging and can offer even better contrast between cartilage and surrounding tissue (muscle, fat, etc.) (Peterfy et al. 2008). Figure 5 showcases the same cross-slice of a single subject across the three different acquisition protocols.

Figure 5. DESS, SPGR and T2 MRI sequences for the same subject. The differences in intensity profiles is apparent, especially between sequence T2 and the rest, suggesting that features learned from these images will yield sufficiently different descriptions of the image information content.

An important preprocessing step in the multi-sequence multi-view case is the intra-subject alignment of the three images. Since each sequence follows a distinct acquisition protocol, there are differences regarding the voxel spacing, image size, direction and origin. However, since the underlying anatomy is the same, non-linear deformations are not an issue (Schneider et al. 2008). Regarding DESS as the reference sequence for each subject, a simple resampling and interpolation operation is sufficient to align all sequences in the same coordinate space. Next, SPGR and T2 sequences undergo the same exact preprocessing steps outlined in Section 2.1.

4. The AMUSE model

4.1. Background on Graph-based SSL

Given a set of $\mathit{n}$ voxels, $n^{\mathcal{A}}$ of which are sampled from an atlas library (labelled) and the remaining $n^{\mathcal{T}}=\mathit{n}-n^{\mathcal{A}}$ are sampled from a target image (unlabelled), $\mathcal{V}$ feature sets are extracted, each corresponding to a different view. The obtained dataset is represented in a column-wise manner as $\{{\mathbf{X}}^{(\mathit{v})}=[{\mathbf{x}}_1^{(v)},\cdots ,{\mathbf{x}}_n^{(v)}],\text{break}{\mathbf{X}}^{(\mathit{v})}\in {\mathbb{R}}^{{\mathit{q}}^{(\mathit{v})}\times \mathit{n}}{\}}_{v=1}^{\mathcal{V}}$, with the corresponding labels being one-hot encoded into the label matrix $\mathbf{Y}=[{\mathbf{Y}}^{\mathcal{A}};{\mathbf{Y}}^{\mathcal{T}}]\in {\mathbb{R}}^{\mathit{n}\times \mathit{c}}$ (components of ${\mathbf{Y}}^{\mathcal{T}}$ are initially set to zero).

Adopting the Gaussian Fields and Harmonic Functions approach(GFHF) (Zhu et al. 2003), an affinity matrix ${\mathbf{W}}^{(\mathit{v})}\in {\mathbb{R}}^{\mathit{n}\times \mathit{n}}$ is calculated for each view, based on a pre-defined similarity measure, corresponding to a graph ${\mathcal{G}}^{(\mathit{v})}=(V^{\mathcal{V}},E^{\mathcal{V}})$, whose vertices correspond to data points and the edge weights encode the spectral similarity of those points within the context of the specific view. The Laplacian of each graph is derived as $L^{(\mathit{v})}={\mathbf{D}}^{(\mathit{v})}-{\mathbf{W}}^{(\mathit{v})}$ where ${\mathbf{D}}^{(\mathit{v})}$ is a diagonal matrix, with $D_{ii}^{(v)}={\sum }_{j=1}^n{W}_{ij}^{(v)}$ (Mohar 1991).

In the single-view case $(\mathcal{V}=1)$, GFHF utilizes the affinity matrix $\mathbf{W}$ and the corresponding Laplacian $\mathbf{L}$ to propagate labels from labeled data to unlabeled ones, under the constraint that already labeled data will not be affected. The underlying principle lies on the manifold assumption,i.e., data in a high-dimensional space share the same labels if they are spectrally similar within a lower-dimensional manifold embedded in that space (for which graph $\mathcal{G}$ serves as a proxy.) This idea can be rigorously defined in the following optimization objective

(1) $\underset{F}{min}\propto \sum _{\mathit{i}=1}^{{\mathit{n}}_{\mathcal{A}}}\parallel {\mathbf{f}}_i-{\mathbf{y}}_i{\parallel }_2^2+\sum _{\mathit{i},\mathit{j}=1}^n{W}_{\mathit{ij}}\parallel {\mathbf{f}}_i-{\mathbf{f}}_j{\parallel }_2^2$(1)

where $\mathbf{F}\in {\mathbb{R}}^{\mathit{n}\times \mathit{c}}$ denotes the learned labels, the first term above corresponds to the loss, and the second one acts as the manifold regulariser. Casting (1) as a trace minimisation problem, we get

(2) $\begin{array}{cc}\underset{F}{min}& \mathit{Tr}({\mathbf{F}}^T\mathbf{LF})\\ \mathrm{s}.\mathrm{t}& {\mathbf{F}}^{\mathcal{A}}={\mathbf{Y}}^{\mathcal{A}}\end{array}$(2)

4.2. Multiple-view integration

In the multi-view setting, given multiple graphs ${\mathcal{G}}^{(\mathit{v})}$ and the corresponding affinity matrices ${\mathbf{W}}^{(\mathit{v})}$, a straightforward way to integrate the multiple views is via a linear combination, yielding:

(3) $\begin{array}{cc}\underset{\mathbf{F},\alpha }{min}& \sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v\mathit{Tr}({\mathbf{F}}^T{\mathbf{L}}^{(\mathit{v})}\mathbf{F})\\ \mathrm{s}.\mathrm{t}& \mathbf{S}=\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathcal{V})}{\mathbf{F}}^{(\mathcal{A})}={\mathbf{Y}}^{\mathcal{A}},{\alpha }^T\mathbf{1}=1,\mathit{\alpha }\ge \mathbf{0}\end{array}$(3)

where ${\mathbf{L}}_{\mathbf{S}}$ corresponds to the Laplacian of the explicit combined affinity matrix $\mathbf{S}={\sum }_{v=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathit{v})}$. The optimal ${\alpha }^{\ast }$ resulting from the above problem is often too sparse to effectively take advantage of the information content across the multiple views. To this end, the equality constraint is relaxed to a quadratic penalty term, as follows:

(4) $\begin{array}{cc}\underset{\mathbf{F},\alpha ,\mathbf{S}}{min}& \mathit{Tr}({\mathbf{F}}^T{\mathbf{L}}^{(\mathit{v})}\mathbf{F})+\mathit{\lambda }\parallel \mathbf{S}-\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathit{v})}{\parallel }_F^2\\ \mathrm{s}.\mathrm{t}.& {\mathbf{F}}^{(\mathcal{A})}={\mathbf{Y}}^{\mathcal{A}},{\alpha }^T\mathbf{1}=1,\mathit{\alpha }\ge \mathbf{0}\\ \mathit{\hspace{0.17em}}& \mathbf{S1}=\mathbf{1},\mathbf{S}\ge \mathbf{0}\end{array}$(4)

4.3. AMUSE optimization

Since the objective function in(4) cannot be efficiently solved for all three variables simultaneously, an alternating optimisation process is proposed, whereby two out of the three variables are fixed, while the third one is updated, in an iterative fashion. Algorithm 2 summarises the steps outlined in the following paragraphs.

4.3.1. Update F,fixing S, α

With $\mathbf{S},\mathit{\alpha }$ fixed, the problem is equivalent to the single view case(2), where $\mathbf{L}={\mathbf{L}}_{\mathbf{S}}$. Compartmentalizing the Laplacian into a $2\times 2$ block matrix format and separating the labelled and unlabelled entries, we get the following optimisation objective:

(5) $\begin{array}{cc}\underset{\mathbf{F},\alpha }{min}& \mathit{Tr}\left({\left[\begin{array}{c}{\mathbf{F}}^{\mathcal{A}}\\ {\mathbf{F}}^{\mathcal{T}}\end{array}\right]}^{\mathit{T}}\left[\begin{array}{cc}{\mathbf{L}}_{\mathcal{A}\mathcal{A}}& {\mathbf{L}}_{\mathcal{A}\mathcal{T}}\\ {\mathbf{L}}_{\mathcal{T}\mathcal{A}}& {\mathbf{L}}_{\mathcal{T}\mathcal{T}}\end{array}\right]\left[\begin{array}{c}{\mathbf{F}}^{\mathcal{A}}\\ {\mathbf{F}}^{\mathcal{T}}\end{array}\right]\right)\\ \mathrm{s}.\mathrm{t}& {\mathbf{F}}^{\mathcal{A}}={\mathbf{Y}}^{\mathcal{A}}\end{array}$(5)

Taking the partial derivative of the objective with respect to ${\mathbf{F}}^{\mathcal{T}}$ and setting to zero, the following solution is obtained:

(6) ${\mathbf{F}}^{\mathcal{T}}=-({\mathbf{L}}_{\mathcal{T}\mathcal{T}}{)}^{-1}{\mathbf{L}}_{\mathcal{T}\mathcal{A}}{\mathbf{Y}}^{\mathcal{A}}$(6)

Finally, $\mathbf{F}$ is updated as

(7) $\mathbf{F}=\left[\begin{array}{c}{\mathbf{Y}}^{\mathcal{A}}\\ {\mathbf{F}}^{\mathcal{T}}\end{array}\right]$(7)

4.3.2. Update α, fixing F,S

The objective of (4) for this round assumes the following formulation:

(8) $\begin{array}{cc}\underset{\mathbf{S}}{min}& \mathit{Tr}({\mathbf{F}}^T{\mathbf{L}}_{\mathbf{S}}\mathbf{F})+\mathit{\lambda }{∥\mathbf{S}-\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathit{v})}∥}_{\mathit{F}}^2\\ \mathrm{s}.\mathrm{t}& \mathbf{S1}=\mathbf{1},\mathbf{S}\ge 0\end{array}$(8)

which can be further algebraically manipulated as follows:

$\mathit{Tr}({\mathbf{F}}^T{\mathbf{L}}_{\mathbf{S}}\mathbf{F})+\mathit{\lambda }{∥\mathbf{S}-\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathit{v})}∥}_{\mathit{F}}^2$

$=\frac{1}{2}\sum _{\mathit{i},\mathit{j}=1}^n{D}_{\mathit{ij}}{s}_{\mathit{ij}}+\mathit{\lambda }\sum _{\mathit{i},\mathit{j}=1}^n{\left({\mathit{s}}_{\mathit{ij}}-\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{W}_{ij}^{(v)}\right)}^2$

$=\mathit{\lambda }\sum _{\mathit{i},\mathit{j}=1}^n{\left({\mathit{s}}_{\mathit{ij}}-\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{W}_{ij}^{(v)}+\frac{{\mathit{D}}_{\mathit{ij}}}{4\mathit{\lambda }}\right)}^2+\mathrm{constant}$

where $D_{\mathit{ij}}=\parallel {\mathbf{f}}_i-{\mathbf{f}}_j{\parallel }_2^2$. The simplex constraint ${\sum }_{j=1}^n{s}_{\mathit{ij}}=1$ essentially means that the rows of $\mathbf{S}$ can be treated independently of one another, leading to the below optimisation problem for $\mathit{i}=1,\cdots ,\mathit{n}$

(9) $\begin{array}{cc}\underset{{\mathbf{S}}_i}{min}& \sum _{\mathit{j}=1}^n{\left({\mathit{s}}_{\mathit{ij}}-\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{W}_{\mathit{ij}}+\frac{{\mathit{D}}_{\mathit{ij}}}{4\mathit{\lambda }}\right)}^2\\ \mathrm{s}.\mathrm{t}& \sum _{\mathit{j}=1}^n{s}_{\mathit{ij}}=1,s_{\mathit{ij}}\ge 0\end{array}$(9)

Denoting $T_{\mathit{ij}}={\sum }_{v=1}^{\mathcal{V}}{\alpha }_v{W}_{\mathit{ij}}-\frac{{\mathit{D}}_{\mathit{ij}}}{4\mathit{\lambda }}$ for notational brevity, the optimisation problem for each row $\mathit{i}$ of $\mathbf{S}$ can finally be formulated as

(10) $\begin{array}{cc}\underset{\mathbf{S}}{min}& \frac{1}{2}\parallel {\mathbf{S}}_i-{\mathbf{T}}_i{\parallel }_2^2\\ \mathrm{s}.\mathrm{t}& {\mathbf{1}}^T{\mathbf{S}}_i=1,s_{\mathit{ij}}\ge 0\end{array}$(10)

This constitutes a projection problem and can be efficiently solved utilising an accelerated projected gradient method. The details of the solver’s implementation can be found in (Huang et al. 2015).

Algorithm 2:

AMUSE Optimization Algorithm

Data: affinity matrices $\{{\mathbf{W}}^{(\mathit{v})}\in {\mathbb{R}}^{\mathit{n}\times \mathit{n}}{\}}_{v=1}^{\mathcal{V}}$, label matrix ${\mathbf{Y}}^{\mathcal{A}}\in {\mathbb{R}}^{\mathit{n}\times \mathit{c}}$, regularisation parameter $\mathit{\lambda }\in \mathbb{R}$

Result: affinity matrix of multi-view graph $\mathbf{S}\in {\mathbb{R}}^{\mathit{n}\times \mathit{n}}$, label matrix of unlabelled data ${\mathbf{F}}^{\mathcal{T}}\in {\mathbb{R}}^{\mathcal{T}\times c}$, linear combination parameters $\mathit{\alpha }\in {\mathbb{R}}^{\mathcal{V}}$

1. Initialize ${\mathbf{S}}_0=\frac{1}{\mathcal{V}}{\sum }_{v=1}^{\mathcal{V}}{\mathbf{W}}^{(\mathit{v})}$

2. Initialize ${\alpha }_v=\frac{1}{\mathcal{V}}\mathrm{\forall }\mathit{v}=1,\dots \mathcal{V}$;

3. while convergence criterion not met do

4. Update $\mathbf{F}$ with Eq. (7);

5. Update $\mathbf{S}$ with Eq. (10);

6. Update $\mathit{\alpha }$ with Eq. (14);

4.3.3. Update S, fixing F,α

With $\mathbf{F},\mathit{\alpha }$ fixed, the original problem (4) becomes

(11) $\begin{array}{cc}\underset{\alpha }{min}& {∥\mathbf{S}-\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathit{v})}∥}_{\mathit{F}}^2\\ \mathrm{s}.\mathrm{t}& {\alpha }^T\mathbf{1}=1,\mathit{\alpha }\ge 0\end{array}$(11)

and can be reformulated in a more compact form as

(12) $\begin{array}{cc}\underset{\alpha }{min}& \parallel \mathit{vec}(\mathbf{S})-\mathrm{\Lambda }\mathit{\alpha }{\parallel }_2^2\\ \mathrm{s}.\mathrm{t}& {\alpha }^T\mathbf{1}=1,\mathit{\alpha }\ge 0\end{array}$(12)

where $\mathit{vec}(\cdot )$ is the vectorisation (flattening) operator, $\mathit{vec}(\mathbf{S})\text{break}\in {\mathbb{R}}^{n^2}$ and $\mathrm{\Lambda }\in {\mathbb{R}}^{n^2\times \mathcal{V}}$ is

$\mathrm{\Lambda }=\left[\begin{array}{ccc}|& \mathit{\hspace{0.17em}}& |\\ \mathit{vec}({\mathbf{W}}^{(1)})& \cdots & \mathit{vec}({\mathbf{W}}^{(\mathcal{V})})\\ |& \mathit{\hspace{0.17em}}& |\end{array}\right]$

Observing that $\mathrm{\Lambda \alpha }=\mathit{vec}({\sum }_{v=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathit{v})})$, the above problem can finally be cast in the form of a Quadratic Program (QP), as showcased below:

(13) $\begin{array}{cc}\underset{\alpha }{min}& {\alpha }^T{\mathrm{\Lambda }}^T\mathrm{\Lambda \alpha }-2\mathit{vec}{(\mathbf{S})}^T\mathrm{\Lambda }\mathit{\alpha }=\underset{\alpha }{min}{\alpha }^T\mathbf{P}\mathit{\alpha }-{\mathbf{q}}^T\mathit{\alpha }\\ \mathrm{s}.\mathrm{t}& {\alpha }^T\mathbf{1}=1,\mathit{\alpha }\ge 0\end{array}$(13)

where $\mathbf{P}={\mathrm{\Lambda }}^T\mathrm{\Lambda }$ and $\mathbf{q}=2{\mathrm{\Lambda }}^T\mathit{vec}(\mathbf{S})$. Introducing an auxiliary $\mathit{\beta }\in {\mathbb{R}}^{\mathcal{V}}$ variable and an additional equality constraint, problem (13) becomes

(14) $\begin{array}{cc}\underset{\mathit{\alpha },\mathit{\beta }}{min}& {\alpha }^T\mathbf{P}\mathit{\alpha }-{\mathbf{q}}^T\mathit{\beta }\\ \mathrm{s}.\mathrm{t}& {\alpha }^T\mathbf{1}=1,\mathit{\alpha }\ge 0,\mathit{\alpha }=\mathit{\beta }\end{array}$(14)

This problem can be efficiently solved using the Alternating Direction Method of Multipliers(ADMM) (Boyd et al. 2010), whereby the original problem is split into optimising $\mathit{\alpha }$ and $\mathit{\beta }$ independently. The exact derivations for the final update rules can be found in (Nie et al. 2020).

5. Multi-view graph constructions

In this section, we describe the global data sampling process and the sparse graph construction of the multiple views.

5.1. Data sampling

A spatially stratified sampling process is used to yield the two global datasets of labelled ${\mathbf{X}}_g^{\mathcal{A}}$ and unlabeled ${\mathbf{X}}_g^{\mathcal{T}}$, respectively. The first step involves the application of KMeans clustering on the spatial coordinates of all voxels ${\mathbf{x}}^{\mathcal{T}}\in \mathit{RO}{I}^{\mathcal{T}}$, yielding $n_g^{\mathcal{T}}$ cluster centers. After interpolating those centers to the closest grid voxels, we obtain a global target (unlabeled) dataset ${\mathbf{X}}_g^{\mathcal{T}}=\{{\mathbf{x}}_g^{\mathcal{T}}(\mathit{i}){\}}_{i=1}^{n_g^{\mathcal{T}}}$. The above step guarantees a sufficient coverage of the entire cartilage ROI.

Next, the same process is repeated, this time clustering the coordinates of voxels in ${\mathbf{X}}_g^{\mathcal{T}}$ into $n_s^{\mathcal{T}}$ clusters $(n_s^{\mathcal{T}}<n_g^{\mathcal{T}})$, thus forming the subset ${\mathbf{X}}_{g_s}^{\mathcal{T}}\subset {\mathbf{X}}_g^{\mathcal{T}}$. Having established a common coordinate space between the target image $\mathcal{T}$ and Atlas Library $\{{\mathcal{A}}_i^{\mathcal{T}},{\mathcal{L}}_i^{\mathcal{T}}{\}}_{\mathit{i}=1}^{{\mathit{n}}_{\mathcal{A}}}$ through registration (Section 2.2), all atlases ${\mathcal{A}}_i$ are sampled at spatially correspondent locations specified by the coordinates of ${\mathbf{X}}_{g_s}^{\mathcal{T}}$. A global atlas (labeled) dataset ${\mathbf{X}}_g^{\mathcal{A}}=\{{\mathbf{x}}_g^{{\mathcal{A}}_i}(\mathit{j}),\mathit{j}:{\mathbf{x}}_{g_s}^{\mathcal{T}}(\mathit{j})\in {\mathbf{X}}_{g_s}^{\mathcal{T}}{\}}_{\mathit{i}=1}^{{\mathit{n}}_{\mathcal{A}}},(n_g^{\mathcal{A}}=n_s^{\mathcal{T}}\cdot n_{\mathcal{A}})$ is thus formed, along with the corresponding label matrix ${\mathbf{Y}}_g^{\mathcal{A}}\in {\mathbb{R}}^{n_g^{\mathcal{A}}\times \mathit{c}}$. In restricting the sampling of atlas voxels to spatial locations corresponding to already sampled target ones, we maintain the patch-based directive: that is, estimating the target voxels’ labels from those of spatially aligned ones in the atlas library.

The cardinalities $n_g^{\mathcal{T}}=|{\mathbf{X}}_g^{\mathcal{T}}|$ and $n_g^{\mathcal{A}}=|{\mathcal{X}}_g^{\mathcal{A}}|$ of the unlabelled and labelled datasets respectively, are properly chosen to balance between improved accuracy (larger values) and moderate computational demands (smaller values). In our experiments, we consider a fixed number of global samples $(n_g=n_g^{\mathcal{T}}+n_g^{\mathcal{A}})$.

5.2. Views sparse graph coding

In constructing the graphs encoding data connectivity for each view, we incorporate the concept of sparse neighbourhood (Zang and Zhang 2012). Each voxel ${\mathbf{x}}_i^{(v)}$ is linearly reconstructed from its respective Local Neighborhood $\mathit{LN}({\mathbf{x}}_i^{(v)})$, which comprises the set of $\mathit{k}$-nearest neighbors from the global dataset ${\mathbf{X}}_g^{(v)}=\text{break}[{\mathbf{X}}_g^{\mathcal{T}(v)},{\mathbf{X}}_g^{\mathcal{A}(v)}]\in {\mathbb{R}}^{q^v\times n_g}$. Let ${\mathbf{A}}_i^{(v)}\in {\mathbb{R}}^{q^v\times \mathit{k}},{\mathbf{A}}_i^{(v)}\subset \{{\mathbf{X}}_g^{(v)}\mathrm{\setminus }{\mathbf{x}}_i^{(v)}\}$ denote the subset of $\mathit{k}$-nearest neighbors of ${\mathbf{x}}_i^{(v)}$. The entire sparse graph ${\mathcal{G}}_g^{(v)}$ associated to view $\mathit{v}=1,\dots ,\mathcal{V}$, is constructed via a series of sparse optimizations (Zhang et al. 2015):

(15) $\begin{array}{cc}\underset{{\mathbf{w}}_i^{(v)}}{min}& \parallel {\mathbf{w}}_i^{(v)}{\parallel }_1,\mathrm{\forall }{\mathbf{x}}_i^{(v)}\in {\mathbf{X}}_g^{(v)}\\ \mathrm{s}.\mathrm{t}& \frac{1}{2}\parallel {\mathbf{x}}_i^{(v)}-{\mathbf{A}}_i^{(v)}{\mathbf{w}}_i^{(v)}{\parallel }_2^2\le \mathit{\epsilon }\end{array}$(15)

where $\mathit{\epsilon }>0$ typically set to a small value. Eq.(15) is a $l_1$-minimisation problem under a quadratic constraint, which can be efficiently solved via coordinate descent. Each resulting vector ${\mathbf{w}}_i^{(v)}\in {\mathbb{R}}^k$ is projected to a unit $l_1$-norm $(\parallel {\mathbf{w}}_i{\parallel }_1=1)$ and resized to ${\mathbb{R}}^{n_g}$ with zero-padding. The edge weights for graph ${\mathcal{G}}_g^{(v)}$ form the affinity matrix for view $\mathit{v}$, ${\mathbf{W}}^{(\mathit{v})}=[{\mathbf{w}}_1^{(v)},\dots ,{\mathbf{w}}_{n_g}^{(\mathit{v})}]\in {\mathbb{R}}^{n_g\times n_g}$. Algorithm 3 summarises the above procedure.

Algorithm 3:

Sparse Graph Construction

Data: labelled atlas voxels ${\mathbf{X}}_g^{\mathcal{A}(v)}$, unlabelled target voxels ${\mathbf{X}}_g^{\mathcal{T}(v)}$, number of neighbours $\mathit{k}$

Result: Sparse Graphs’ edge weights ${\mathbf{W}}^{(\mathit{v})}$

1. for $\mathit{v}\in \{1,\cdots ,\mathcal{V}\}$ do

2. identify $\mathit{k}$-nearest neighbours $\mathrm{\forall }{\mathbf{x}}^{(\mathit{v})}\in {\mathbf{X}}_g^{(v)}$;

3. for each ${\mathbf{x}}_i^{(v)}\in {\mathbf{X}}_g^{(v)}$ do

4. solve for ${\mathbf{w}}_i^{(v)}$ (15);

5. ${\mathbf{W}}^{(\mathit{v})}=[\{{\mathbf{w}}_i^{(v)}{\}}_{\mathit{i}=1}^{n_g}{]}^T;$

In all of the above, each voxel ${\mathbf{x}}_i^{(v)}$ is described by its associated Regional SV descriptor corresponding to a distinct view (HOG, LPB, etc). Each weight $w_{ij}^{(v)}$ encodes the spectral similarities of the Regional SVs ${\mathcal{R}}_i$ and ${\mathcal{R}}_j$ for view $\mathit{v}$.

6. Proposed multi-view labeling schemes

We now propose two label propagation mechanisms, the MV-RegLp and MV-HyLP methods. The first one considers voxel descriptors’ similarity at the regional level, while the second one integrates an additional spatial scale by incorporating the notion of Label Map Estimates.

6.1. Multi-view regional label propagation (MV-RegLP)

Regional Label Propagation (MV-RegLP) operates exclusively on the regional level, whereby each voxel ${\mathbf{x}}_i$ is characterized by its associated ${\mathcal{R}}_i$ multi-view feature descriptors $\{{\mathbf{H}}_i^{(v)}\in {\mathbb{R}}^{q^v}{\}}_{v=1}^{\mathcal{V}}$, and assigned a label ${\mathbf{y}}_i$, (the label of its central voxel ${\mathbf{x}}_i$).

The labelling of MV-RegLP proceeds along the following steps:

1. Compute the affinity matrices $\{{\mathbf{W}}^{(\mathit{v})}{\}}_{v=1}^{\mathcal{V}}$ for each view through Algorithm 3.

2. Solve the multi-view optimisation problem (4) via the AMUSE (Algorithm 2) to obtain the combined affinity matrix

   (16) ${\mathbf{W}}_g=\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{\mathbf{W}}^{(\mathit{v})}\in {\mathbb{R}}^{n_g\times n_g}$(16)

   along with the multi-view coefficients $\mathit{\alpha }\in {\mathbb{R}}^{\mathcal{V}}$.

3. Obtain the Laplacian of ${\mathbf{W}}_g$ as $\mathbf{L}=\mathbf{I}-{\mathbf{W}}_g$ to yield the label propagation weights

   (17) ${\tilde{\mathbf{W}}}_{\mathit{g}}=-{\mathbf{L}}_{\mathcal{T}\mathcal{T}}^{-1}{\mathbf{L}}_{\mathcal{T}\mathcal{A}}\in {\mathbb{R}}^{n_g^{\mathcal{T}}\times n_g^{\mathcal{A}}}$(17)

   which encode the regional transfer mechanism from labelled to unlabelled data.

4. Target voxels ${\mathbf{x}}_g^{\mathcal{T}}\in {\mathbf{X}}_g^{\mathcal{T}}$ assume their labels via

   (18) ${\hat{\mathbf{Y}}}_{\mathit{g}}^{\mathcal{T}}={\tilde{\mathbf{W}}}_{\mathit{g}}{\mathbf{Y}}_g^{\mathcal{A}}$(18)

   where ${\mathbf{Y}}_g^{\mathcal{A}}=[{\mathbf{y}}_1^{\mathcal{A}}\dots {\mathbf{y}}_{n_g^{\mathcal{A}}}^{\mathcal{A}}{]}^T\in R^{(n_g^{\mathcal{A}}\times \mathit{c})}$ and ${\hat{\mathbf{Y}}}_{\mathit{g}}^{\mathcal{T}}=[{\mathbf{y}}_1^{\mathcal{T}}\text{break}\dots {\mathbf{y}}_{n_g^{\mathcal{T}}}^{\mathcal{T}}{]}^T\in R^{(n_g^{\mathcal{T}}\times \mathit{c})}$.

6.2. Multi-view Hybrid Label Propagation (MV-HyLP)

The feature descriptors ${\mathbf{H}}_i^{(v)},\mathit{v}=(1,\cdots ,\mathcal{V})$ associated with a SV ${\mathcal{R}}_i$ provide a voxel’s ${\mathbf{x}}_i$ characterization at a broad spatial scale. To leverage the information encoded within the constituent patches at a finer spatial level, for each view $\mathit{v}=1,\cdots ,\mathcal{V}$ we define the patch-to-patch spectral similarity for ${\mathbf{p}}_i,{\mathbf{p}}_j$ via the chi-squared kernel (Zhang et al. 2006):

(19) $s^{(\mathit{v})}(\mathit{i},\mathit{j})=s^{(\mathit{v})}({\mathbf{p}}_i,{\mathbf{p}}_j)=-\sum _{\mathit{k}=1}^q\frac{{(h_{i_k}^{(\mathit{v})}-h_{j_k}^{(\mathit{v})})}^2}{h_{i_k}^{(\mathit{v})}+h_{j_k}^{(\mathit{v})}}$(19)

Next, a region-to-region label transfer mechanism is developed, whereby the label estimates ${\hat{\mathbf{y}}}_{i_k}^{\mathcal{T}(\mathit{v})}$ of target patches ${\mathbf{p}}_{i_k}^{\mathcal{T}(\mathit{v})}\in {\mathcal{R}}_i^{\mathcal{T}(v)}$, are computed as a weighted combination of all label vectors ${\mathbf{y}}_{j_l}^{\mathcal{T}(\mathit{v})}\in {\mathcal{R}}_j^{\mathcal{A}}$:

(20) ${\hat{\mathbf{y}}}_{i_k}^{\mathcal{T}(\mathit{v})}=\frac{{\sum }_{l=1}^{27}{s}^{(\mathit{v})}(i_1,j_l)({\mathbf{y}}_{j_l}^{\mathcal{A}}{)}^T}{{\sum }_{l=1}^{27}{s}^{(\mathit{v})}(i_1,j_l)}={\tilde{\mathbf{s}}}_{\mathit{ij}}^{(\mathit{v})}(i_1,\cdot ){\mathbf{Y}}_j^{\mathcal{A}}$(20)

where ${\tilde{\mathbf{s}}}_{\mathit{ij}}^{(\mathit{v})}(i_1,\cdot )=[{\tilde{s}}^{(\mathit{v})}(i_1,j_1),\dots ,{\tilde{s}}^{(\mathit{v})}(i_1,j_{27})]\in {\mathbb{R}}^{27}$ is the normalized to unit sum patch similarity vector and ${\mathbf{Y}}_j^{\mathcal{A}}$ the label matrix associated with ${\mathcal{R}}_j^{\mathcal{A}}$. Then, the Label Map Estimate (LME) ${\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}(\mathit{v})}(\mathit{j})\in {\mathbb{R}}^{27\times \mathit{c}}$ of ${\mathcal{R}}_i^{\mathcal{T}}$ is produced as:

(21) ${\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}(\mathit{v})}(\mathit{j})={\tilde{\mathbf{S}}}_{\mathit{ij}}^{(\mathit{v})}{\mathbf{Y}}_j^{\mathcal{A}},\mathit{v}=(1,\dots ,\mathcal{V})$(21)

where ${\tilde{\mathbf{S}}}_{\mathit{ij}}^{(\mathit{v})}=[{\tilde{\mathbf{s}}}_{\mathit{ij}}^{(\mathit{v})}(i_1,\cdot ),\dots ,{\tilde{\mathbf{s}}}_{\mathit{ij}}^{(\mathit{v})}(i_{27},\cdot )]\in {\mathbb{R}}^{27\times 27}$ encodes the region-to-region similarity between ${\mathcal{R}}_i^{\mathcal{T}}$ and ${\mathcal{R}}_j^{\mathcal{A}}$, with the index $\mathit{j}$ identifying the Regional SV ${\mathcal{R}}_j^{\mathcal{A}}$ as the source of label information for ${\mathcal{R}}_i^{\mathcal{T}}$. The above equation provides an effective way of encoding the patch-level label distribution across regional SVs.

The MV-HyLP method integrates label information at two spatial scales, proceeding as follows:

1. The first step is identical to MV-RegLP, obtaining the sparse graphs ${\mathcal{G}}_g^{(v)}$ and solving 4 via AMUSE (Algorithm 2) to obtain the label propagation weights ${\tilde{\mathbf{W}}}_{\mathit{g}}$. This step provides the regional-level description on the data.

2. For each view $\mathit{v}=1,\dots ,\mathcal{V}$, a Regional SV ${\mathcal{R}}_i^{\mathcal{T}}$ is associated with a collection of SVs $\{{\mathcal{R}}_j^{\mathcal{A}},\mathit{j}:{\mathbf{x}}_j^{\mathcal{A}}\in \mathit{LN}({\mathbf{x}}_i^{\mathcal{T}})\}$. Each ${\mathcal{R}}_j^{\mathcal{A}}$ produces a label map estimate ${\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}(\mathit{v})}(\mathit{j})$ which is computed using (21). Complementing the previous step, this one offers a finer level of description within each Regional SV, by leveraging the patch-to-patch similarity matrices ${\tilde{\mathbf{S}}}_{\mathit{ij}}^{\mathit{v}}$

3. Aggregate the view-specific LMEs ${\mathbf{Y}}_i^{\mathcal{T}(v)}$ over the different views using the multi-view coefficients $\mathit{\alpha }$ to obtain:

   (22) ${\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}}(\mathit{j})=\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{\mathbf{Y}}_i^{\mathcal{T}(v)}(\mathit{j})$(22)

4. Translating the manifold assumption from points – labels to Regional SVs – LMEs, aggregate the derived label estimates for ${\mathcal{R}}_i^{\mathcal{T}}$, into a final label map, by utilizing the label propagation coefficients of ${\tilde{\mathbf{W}}}_{\mathit{g}}$:

   (23) ${\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}}=\sum _{j\in LN({\mathbf{x}}_g^{\mathcal{T}}(\mathit{i}))}{\tilde{w}}_{\mathit{ij}}{\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}}(\mathit{j})$(23)

In view of (21) and (22) ${\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}}$ can be expressed in a more compact matrix form as:

(24a) ${\hat{\mathbf{Y}}}_{\mathit{i}}^{\mathcal{T}}=\sum _j{\mathbf{Q}}_{\mathit{ij}}{\mathbf{Y}}_j^{\mathcal{A}}$(24a)

(24b) ${\mathbf{Q}}_{\mathit{ij}}={\tilde{w}}_{\mathit{ij}}\sum _{\mathit{v}=1}^{\mathcal{V}}{\alpha }_v{S}_{ij}^{(v)}$(24b)

The terms ${\mathbf{Q}}_{\mathit{ij}}\in {\mathbb{R}}^{(27\times 27)}$ in (24b) provide the label transfer between ${\mathcal{R}}_i^{\mathcal{T}}$ and the target region ${\mathcal{R}}_j^{\mathcal{A}}$. The first term ${\tilde{w}}_{\mathit{ij}}$ represents their affinity at the regional level, while the second one is a weighting average of the region-to-region similarities across the different views.

(5) Finally, subsuming the label estimates of all target voxels in ${\mathbf{X}}_g^{\mathcal{T}}$, we obtain,

(25) $[{\hat{\mathbf{Y}}}_{\mathit{g}}^{\mathcal{T}}]=\mathbf{Q}[{\mathbf{Y}}_g^{\mathcal{A}}]$(25)

where $\mathbf{Q}=[{\mathbf{Q}}_{\mathit{ij}}]\in {\mathbb{R}}^{(27n_g^{\mathcal{T}}\times 27n_g^{\mathcal{A}})}$, and $[{\hat{\mathbf{Y}}}_{\mathit{g}}^{\mathcal{T}}]=[({\hat{\mathbf{Y}}}_1^{\mathcal{T}}{)}^T,\text{break}\mathit{ldots},\text{(}{\hat{\mathbf{Y}}}_{n_g^{\mathcal{T}}}^{\mathcal{T}}{)}^T{]}^T\in {\mathbb{R}}^{(27\cdot n_g^{\mathcal{T}}\times \mathit{c})}$, $[{\hat{\mathbf{Y}}}_{\mathit{g}}^{\mathcal{A}}]=[({\hat{\mathbf{Y}}}_1^{\mathcal{A}}{)}^T,\dots ,({\hat{\mathbf{Y}}}_{n_g^{\mathcal{A}}}^{\mathcal{T}}{)}^T{]}^T\in \text{break}{\mathbb{R}}^{(27\cdot n_g^{\mathcal{A}}\times \mathit{c})}$ are the label matrices of all Regional SVs for each voxel ${\mathbf{x}}_g^{\mathcal{T}}\in {\mathbf{X}}_g^{\mathcal{T}}$ and ${\mathbf{x}}_g^{\mathcal{A}}\in {\mathbf{X}}_g^{\mathcal{A}}$, respectively.

7. Out-of-sample label propagation

Stage-2 of the proposed method deals with induction on the remaining unlabelled target voxels. As in stage-1, out-of-sample voxels assume their labels via sparse reconstruction from their closest neighbours in the in-sample set. The iterative generation of out-of-sample batches relies on the creation of a 3D tetrahedral mesh, which allows the sampling of voxels at locations dependent on the ones of previously labelled samples. The mesh is progressively densified for each successive iteration until the entire target ROI is classified. The out-of-sample labelling proceeds as follows:

1. Out-of-sample batch generation: At iteration $\mathit{t}=0$, a 3D Delaunay tetrahedral mesh is constructed, with vertices corresponding to voxels in ${\mathbf{X}}_g$. Beginning stage-2 $(\mathit{t}=1,2,\dots )$, the centroids of the mesh tetrahedra are extracted and interpolated to the closest grid point to yield ${\mathbf{X}}_o^{(v)}=\{{\mathbf{x}}_o^{(v)}(\mathit{i})\}\in {\mathbb{R}}^{q^v\times n_o}$, with $n_o$ the batch size, for views $\mathit{v}=1,\dots ,\mathcal{V}$.

2. Local dataset generation: For every ${\mathbf{x}}_o^{(v)}\in {\mathbf{X}}_o^{(v)}$, the four parent vertices in the target image are collected to form the local dataset ${\mathbf{X}}_l^{\mathcal{T}(v)}=\{{\mathbf{x}}_o^{\mathcal{T}(v)}(\mathit{i}),\mathit{i}=1,\dots ,n_l^{\mathcal{T}}\}\text{break}\in {\mathbb{R}}^{q^v\times n_l^{\mathcal{T}}}$. The respective spatially correspondent voxels of those in ${\mathbf{X}}_l^{\mathcal{T}}$ and ${\mathbf{X}}_o$ (centroids & parent vertices) in the atlas library are also collected, to form a second local dataset ${\mathbf{X}}_l^{\mathcal{A}(v)}\in {\mathbb{R}}^{q^v\times n_l^{\mathcal{A}}}$. The label matrices associated with these two local datasets are, accordingly, ${\hat{\mathbf{Y}}}_{\mathit{l}}^{\mathcal{T}}\in {\mathbb{R}}^{(n_l^{\mathcal{T}}\times \mathit{c})}$ and ${\mathbf{Y}}_l^{\mathcal{A}}\in {\mathbb{R}}^{(n_l^{\mathcal{A}}\times \mathit{c})}$.

3. Out-of-sample voxel embedding: The available global ${\mathbf{X}}_g=[{\mathbf{X}}_g^{\mathcal{T}}{\mathbf{X}}_g^{\mathcal{A}}]$, ${\hat{\mathbf{Y}}}_{\mathit{g}}=[({\hat{\mathbf{Y}}}_{\mathit{g}}^{\mathcal{T}}{)}^T,({\mathbf{Y}}_g^{\mathcal{A}}{)}^T{]}^T$ and local ${\mathbf{X}}_l=[{\mathbf{X}}_l^{\mathcal{T}}{\mathbf{X}}_l^{\mathcal{A}}]$, ${\hat{\mathbf{Y}}}_{\mathit{l}}=[({\hat{\mathbf{Y}}}_{\mathit{l}}^{\mathcal{T}}{)}^T,({\mathbf{Y}}_l^{\mathcal{A}}{)}^T{]}^T$ datasets yield complementary information on the data: the former geared towards capturing a larger portion of the class variability across the whole ROI, while the latter contains patterns that are specific to smaller local regions. The existence of multiple views on the data lends yet another dimension on the overall description of the data, rendering it more complete. Leveraging all of the above components, we connect the out-of-sample batch ${\mathbf{X}}_o$ to both the global and local datasets via the learned sparse graphs. Particularly, for each view $(\mathit{v}=1,\dots \mathcal{V})$, the sparse graphs $\{{\mathcal{G}}_{og}^{(v)}({\mathbf{X}}_o^{(v)},{\mathbf{X}}_g^{(v)}),{\mathcal{G}}_{ol}^{(v)}({\mathbf{X}}_o^{(v)},{\mathbf{X}}_l^{(v)})\}$ are constructed according to (15), where each ${\mathbf{x}}_o^{(v)}\in {\mathbf{X}}_o^{(v)}$ is decomposed into a linear combination of its sparse neighbours in ${\mathbf{X}}_g^{(v)},{\mathbf{X}}_l^{(v)}$, respectively.

   The resulting affinity matrices $\{{\tilde{\mathbf{W}}}_{\mathit{og}}^{(\mathit{v})}\in {\mathbb{R}}^{n_g\times n_o},\text{break}{\tilde{\mathbf{W}}}_{\mathit{ol}}(\mathit{v})\in {\mathbb{R}}^{n_l\times n_o}{\}}_{v=1}^{\mathcal{V}}$, with columns restricted to unit $l_1$-norm, are then aggregated via the coefficients $\mathit{\alpha }$ to yield the combined global and local multi-view matrices ${\tilde{\mathbf{W}}}_{\mathit{og}}={\sum }_{v=1}^{\mathcal{V}}{\alpha }_v{\tilde{\mathbf{W}}}_{\mathit{og}}^{(\mathit{v})}$ and ${\tilde{\mathbf{W}}}_{\mathit{ol}}={\sum }_{v=1}^{\mathcal{V}}{\alpha }_v{\tilde{\mathbf{W}}}_{\mathit{ol}}^{(\mathit{v})}$ as illustrated in Figure 6.

   Figure 6. Out-of-sample label propagation: two sparse graphs ${\tilde{\mathbf{W}}}_{\mathit{og}}^{(\mathit{v})},{\tilde{\mathbf{W}}}_{\mathit{ol}}^{(\mathit{v})}$ are constructed for each view $\mathit{v}$ of the out-of-sample batch ${\mathbf{X}}_o^{(v)}$, connecting it to the global ${\mathbf{X}}_g^{(v)}$ and local ${\mathbf{X}}_l^{(v)}$ datasets respectively. The multi-view coefficients $\mathit{\alpha }$ are used to yield the final learned sparse graphs facilitating the labelling of voxels in ${\mathbf{X}}_o$.

   

   Depending on the label transfer mechanism employed, the induction on out-of-sample voxels proceeds by extracting two sets of labels (global + local):

   (26a) ${\hat{\mathbf{Y}}}_{\mathit{og}}={\tilde{\mathbf{W}}}_{\mathit{og}}^{\mathit{T}}{\mathbf{Y}}_g,{\hat{\mathbf{Y}}}_{\mathit{ol}}={\tilde{\mathbf{W}}}_{\mathit{ol}}^{\mathit{T}}{\mathbf{Y}}_l(\mathit{MV}-\mathit{RegLP})$(26a)
   (26b) $[{\hat{\mathbf{Y}}}_{\mathit{og}}]={\mathbf{Q}}_{\mathit{og}}[{\mathbf{Y}}_g],[{\hat{\mathbf{Y}}}_{\mathit{ol}}]={\mathbf{Q}}_{\mathit{ol}}[{\mathbf{Y}}_l](\mathit{MV}-\mathit{HyLP})$(26b)

   with ${\mathbf{Q}}_{\mathit{og}},{\mathbf{Q}}_{\mathit{ol}}$ derived according to (24b). Finally, the complete label estimation is formulated by considering the convex combination of both the globally and locally derived components:

   (27a) ${\hat{\mathbf{Y}}}_{\mathit{o}}=\mathit{\mu }{\tilde{\mathbf{W}}}_{\mathit{og}}^{\mathit{T}}{\mathbf{Y}}_g+(1-\mathit{\mu }){\tilde{\mathbf{W}}}_{\mathit{ol}}^{\mathit{T}}{\mathbf{Y}}_l(\mathit{MV}-\mathit{RegLP})$(27a)
   (27b) $[{\hat{\mathbf{Y}}}_{\mathit{o}}]=\mathit{\mu }{\mathbf{Q}}_{\mathit{og}}[{\mathbf{Y}}_g]+(1-\mathit{\mu }){\mathbf{Q}}_{\mathit{ol}}[{\mathbf{Y}}_l](\mathit{MV}-\mathit{HyLP})$(27b)

   where $\mathit{\mu }\in [0,1]$ adjusts the contribution of global versus local label information, respectively.

4. Mesh densification & termination: At the end of each iteration, the 3D mesh is densified by adding the current centroids to its list of vertices, drawing new edges between child(centroid) and parent nodes. A new round of out-of-sample label propagation begins by locating the centroids of the new refined tetrahedra. This process terminates when a certain percent of those tetrahedra ($\approx 90\mathrm{\%}$) fits entirely in a $5\times 5\times 5$ cell(Algorithm 4).

5. Majority Voting & Post-processing: The remaining unlabelled voxels dispersed through the ROI are guaranteed to be surrounded by multiple labelled ones, owning to the progressive densification of the initial mesh. Thus, a simple Majority Voting rule over the vertices of the encompassing tetrahedra is sufficient to yield accurate classification results. Finally, a binary morphological closing – opening is performed to smooth the surfaces of the segmented structures, followed by a pass through a neighbour voting filter, to rectify potential isolated misclassifications.

Algorithm 4:

Iterative sampling procedure

Data: Global dataset ${\mathbf{X}}_g^{\mathcal{T}}$

1. Construct 3D mesh $\mathcal{D}\mathcal{T}$ via triangulation of ${\mathbf{X}}_g^{\mathcal{T}}$;

2. initialise iteration count $\mathit{t}=1$;

3. while termination condition not met do

4. compute new batch of centroids ${\mathbf{X}}_o$;

5. insert ${\mathbf{X}}_o$ in $\mathcal{D}\mathcal{T}$;

6. $\mathit{t}=\mathit{t}+1;$

8. Experimental setup

8.1. Image data

The MRI sequences used in this study were acquired from the publicly available Osteoarthritis Initiative (OAI) repository (Peterfy et al. 2008). The entirety of the dataset is utilised (507 images), accounting for the whole spectrum of joint degradation with regards to the Kellgren-Lawrence (K – L) grade (Biscaldi et al. 2018). For the multi-view setting involving the use of multiple feature descriptors, we use the sagittal 3D Dual-Echo-Steady-State (3D-DESS) sequence with water excitation, while for the multi-view setting dealing with different MRI sequences, we consider the 3D-DESS, T2-weighted (T2) and SPoiled Gradient Recalled echo (SPGR). In all comparative experiments and test cases, the problem is cast as a 5-class classification: Background Tissue (BT): 0, Femoral Bone (FB): 1, Femoral Cartilage (FC): 2, Tibial Bone (TB): 3 and Tibial Cartilage (TC): 4.

8.2. Evaluation measures

Segmentation performance is evaluated using three well-established volumetric measures, namely, the Dice Similarity Coefficient (DSC), the Volumetric Difference (VD) and the Volume Overlap Error (VOE). Denoting $\mathbf{Y}$ as the ground truth label map and $\hat{\mathbf{Y}}$ the estimated one, they are defined as:

(28a) $\mathit{DSC}=100\frac{|\mathbf{Y}\cap \hat{\mathbf{Y}}|}{|\mathbf{Y}|+|\hat{\mathbf{Y}}|}$(28a)

(28b) $\mathit{VOE}=100\left(1-\frac{\mathit{DSC}}{200-\mathit{DSC}}\right)$(28b)

(28c) $\mathit{VD}=100\frac{|\hat{\mathbf{Y}}|-|\mathbf{Y}|}{|\mathbf{Y}|}$(28c)

Since the majority of voxels in a given MRI belong to either Background or Bone classes, we also include the typical classification measures Recall and Precision, to evaluate the sensitivity of our proposed methodology. All measures refer exclusively to the image content delineated by the respective cartilage ROIs.

8.3. Hyperparameter setting & sensitivity analysis

The proposed MV-RegLP and MV-HyLP, integrate multiple components, the most notable being multi-atlas, sparse representation, semi-supervised learning and finally, multi-view learning. Each component’s contribution hinges upon one or more hyperparameters that control its behaviour: the number of selected atlases $n_{\mathcal{A}}$ (Section 2.2), the percentage $\mathit{ss}{l}_{\mathit{perc}}=100\frac{n_g^{\mathcal{A}}}{n_g}$ of labeled data included in the learning process (Section 5.1), the number of $\mathit{k}$ nearest neighbors in constructing the sparse graphs (Section 5.2), the structural regularization term $\mathit{\lambda }$ in the multi-view semi-supervised objective (Section 4.2) and finally, the convex combination parameter $\mathit{\mu }$ controlling the trade-off between global and local datasets in Equation 27(27a) ${\hat{\mathbf{Y}}}_{\mathit{o}}=\mathit{\mu }{\tilde{\mathbf{W}}}_{\mathit{og}}^{\mathit{T}}{\mathbf{Y}}_g+(1-\mathit{\mu }){\tilde{\mathbf{W}}}_{\mathit{ol}}^{\mathit{T}}{\mathbf{Y}}_l(\mathit{MV}-\mathit{RegLP})$(27a) .

Let $(\mathcal{A},\mathcal{L})=\{{\mathcal{A}}_i,{\mathcal{L}}_i{\}}_{i=1}^n$ be the entire image data set $(\mathit{n}=507)$. The proposed MV-KCS follows the local modeling approach of multi-atlas segmentation, creating a particular model for each target image $(\mathcal{T},{\hat{\mathcal{L}}}_{\mathit{T}})=\mathit{f}(\mathcal{L}\mathcal{B}(\mathcal{T}))$, where $\mathcal{L}\mathcal{B}(\mathcal{T})=\{{\mathcal{A}}_j^{\mathcal{T}},{\mathcal{L}}_j^{\mathcal{T}}{\}}_{\mathit{j}=1}^{{\mathit{n}}_{\mathcal{A}}}$ denotes the corresponding atlas library and $\mathit{f}(\cdot )$ signifies the segmentation mapping function implemented by MV-KCS. To this end, we devise a suitable 5-fold cross validation scheme for the hyperparameter selection and the evaluation of segmentation performance, respectively.

The optimal hyperparameter selection proceeds as follows:

1. Each testing fold $({\mathcal{A}}_t^{(r)},{\mathcal{L}}_t^{(r)})=({\mathcal{A}}_{\mathit{t},\mathit{i}},{\mathcal{L}}_{\mathit{t},\mathit{i}}{)}_{i=1}^m$, $\mathit{r}=1,\dots 5$, $\mathit{m}=\mathit{n}/5$, is used as a source for selecting target images, ${\mathcal{T}}_i^{(r)}={\mathcal{A}}_{t,i}^{(r)},{\mathcal{L}}_{{\mathcal{T}}_i}^{(\mathit{r})}={\mathcal{L}}_{t,i}^{(r)}$. The target labels ${\mathcal{L}}_{{\mathcal{T}}_i}^{(\mathit{r})}$ are considered as unknown, and hence, they are not used in the model’s construction.

2. For each ${\mathcal{A}}_{t,i}^{(r)},\mathit{i}=1,\dots ,\mathit{m}$, we select the nearest neighbor image ${\mathcal{A}}_{\nu ,i}^{(r)}\in \{\mathcal{A}\mathrm{\setminus }{\mathcal{A}}_t^{(r)}\}$ which satisfies the condition ${\mathcal{A}}_{\nu ,i}^{(r)}=min\mathit{MSE}({\mathcal{A}}_{t,i}^{(r)},{\mathcal{A}}_{\nu ,i}^{(r)})$ (Section 2.2). A collection of images $({\mathcal{A}}_{\nu }^{(r)},{\mathcal{L}}_{\nu }^{(t)})=\{{\mathcal{A}}_{\nu ,i}^{(r)},{\mathcal{L}}_{v,i}^{(r)}{\}}_{i=1}^m$ forms the validation data set with known labels, used to validate the models in $({\mathcal{A}}_t^{(r)},{\mathcal{L}}_t^{(t)})$

3. For each ${\mathcal{A}}_{\nu ,i}^{(r)}\in {\mathcal{A}}_{\nu }^{(r)},\mathit{i}=1,\dots ,\mathit{m}$, we select its corresponding library $\mathcal{L}\mathcal{B}({\mathcal{A}}_{\nu ,i}^{(r)})=\{{\mathcal{A}}_{tr,j}^{(r)},{\mathcal{L}}_{tr,j}^{(r)}{\}}_{\mathit{j}=1}^{{\mathit{n}}_{\mathcal{A}}}$, drawn from the training part ${\mathcal{A}}_{tr,j}^{(r)}\in \{\mathcal{A}\mathrm{\setminus }({\mathcal{A}}_{t,i}^{(r)}\cup {\mathcal{A}}_{\nu ,i}^{(r)})\}$.

4. Next, we apply MV-KCS to build the models $({\mathcal{A}}_{\nu ,i}^{(r)},{\hat{\mathcal{L}}}_{\mathit{\nu },\mathit{i}}^{(\mathit{t})})=\mathit{f}(\mathcal{L}\mathcal{B}({\mathcal{A}}_{\nu ,i}^{(r)})),\mathit{i}=1,\dots ,\mathit{m}$ and $\mathit{r}=1,\dots ,5$ for different combinations in the parameter grid space. The optimal hyperparameter assumes the value providing the minimum overall error on the validation parts ${\sum }_{\mathit{i},\mathit{r}}\mathcal{D}\mathcal{S}\mathcal{C}({\mathcal{L}}_{\nu ,i}^{(r)},{\hat{\mathcal{L}}}_{\mathit{\nu },\mathit{i}}^{(\mathit{t})})$

Having determined the optimal parameters, the segmentation performance is evaluated using the traditional $5$-fold cross validation. Concretely, for each ${\mathcal{T}}_i^{(r)}\in {\mathcal{A}}_t^{(r)},\mathit{i}=1,\dots ,\mathit{m}$, we select the corresponding libraries $\mathcal{L}\mathcal{B}({\mathcal{T}}_i^{(r)})=\{{\mathcal{A}}_{tr,j}^{(r)},{\mathcal{L}}_{tr,j}^{(r)}{\}}_{\mathit{j}=1}^{{\mathit{n}}_{\mathcal{A}}}$, where the atlases are now drawn from the training parts ${\mathcal{A}}_{tr,i}^{(r)}\in \{\mathcal{A}\mathrm{\setminus }{\mathcal{A}}_t^{(r)}\}$. Then, models are alternately built via MV-KCS for different testing folds $({\mathcal{T}}_i^{(r)},{\hat{\mathcal{L}}}_{{\mathcal{T}}_i}^{(\mathit{r})})=\mathit{f}(\mathcal{L}\mathcal{B}({\mathcal{T}}_i^{(r)})),\text{break}\mathit{i}=1,\dots ,\mathit{m},\mathit{r}=1,\dots ,5$ and the results are averaged across all images in the data set. The parameter validation and performance evaluation are conducted under the MV-HyLP approach with the default multi-view case (multiple descriptors).

8.4. Single-view vs multi-view test cases

We consider a series of experimental cases, organised in three groups, to validate the effect of fusing multiple views on the performance of our models. [noitemsep]

1. Single-view vs Multi-view: This first group consists of a series of 4 test cases, as outlined in Table 2. Our primary concern in performing this test is two-fold: First, compare the trivial single descriptor case against the cases where multiple features are integrated via AMUSE (Section 3). Secondly, monitor the performance improvements obtained by the various feature subsets, each encoding a different aspect of the image content. The inclusion of features is designed with the aim to successively capture information pertaining to shape (HOG), shape and texture (HOG+LBP) and finally shape, texture and geometry (HOG+LBP+GLCM+GLRLM+HCGF). The descriptors GLCM and GLRLM are included together due to their similarity in the computation.

2. Stacked vs Fused Features: This experiment compares the fusion of HOG,LBP and HCGF as carried out by MV-KCS, against the single-view case where those three features are stacked together. This test case is important to determine whether the additional computational cost pertaining to the optimal MV learning is worth paying, or a simple feature stacking suffices. We chose this particular subset of descriptors to evaluate this idea, since they encode substantially different information content (shape, texture, geometry).

3. Multi-feature vs Multi-modal views: In this case we evaluate our two multi-view settings, elaborated in Sections 3.2,3.3. Specifically, we aim to assess how the setting involving multiple feature descriptors compares against the one associated with multiple imaging sequences. The competing models make use of all the available features and modalities respectively. For the multi-feature model, the DESS sequence is set as the default, while the HOG descriptor is utilised in the multi-sequence scenario.

Table 2. Test cases examined for the first group of experiments.

	Descriptor
	HOG	LBP	GLCM	GLRLM	HCGF
Case 1	✓
Case 2	✓	✓
Case 3	✓	✓	✓	✓
Case 4	✓	✓	✓	✓	✓

All cases utilise the MV-HyLP labelling method. The hyperparameters are set to the values determined by sensitivity analysis.

8.5. Competing methods

The efficacy of our proposed framework is evaluated against two supervised patch-based methods, six supervised and three semi-supervised deep learning models, respectively.

8.5.1. Patch-based methods

The Patch-based Sparse Coding (PB_SC) (Zhang et al. 2012) is a popular approach in medical image segmentation. For consistency reasons, $N_s$ and $P_s$ are taken as $15\times 15\times 15$ and $5\times 5\times 5$, respectively. In addition, we implement a variant termed PB_SC^(stacked), whereby each patch is characterized by a stacked feature vector of HOG,LBP,HCGF descriptors. The rationale for the choice is similar to the one previously outlined in Section 8.4. Aside from the feature description of each patch, the construction of the patch library ${\mathbf{P}}_{\mathcal{L}}$ is identical to the original paper.

A more recent variant on the patch-based family of methods is the Patch-based Joint-Label-Fusion (PB-JLF (Nikolopoulos et al. 2020)). Here, the authors opt to perform a more computational intensive registration process, by incorporating an additional deformable transformation after the initial affine step, achieving better correspondence between the respective image structures. In contrast to the two standard patch-based methods described in the above paragraph, all images are mapped to a common template space, offering a substantial boost in speed over alternative methods that operate on each target image space independently. The segmentation map for each target image is obtained via a joint label fusion mechanism, where each atlas is assigned a weight according to a voting process. The final results are obtained after an inverse transformation to the original target image space. The method can be extended in a straight-forward manner to incorporate additional image modalities. In the following experiments, the multi-modal setting will be the default case considered for comparison.

8.5.2. Deep learning-based methods

Since deep learning is a prevalent approach, increasingly gaining popularity in recent years, we opt to compare our method against eight deep networks, developed to specifically tackle medical image segmentation problems. The following paragraphs contain a brief summary of our own implementations of each model under consideration.

8.5.2.1. Triplanar CNN (Prasoon et al. 2013)

This method consists of a triad of CNNs, each operating on a 2D plane ($-\mathit{xy},-\mathit{xz},-\mathit{yz}$) corresponding to the sagittal, coronal and axial views, respectively. For each voxel to be classified, a 2D patch centered around it is extracted from each of the three planes and is fed to the respective network. Each of the three branches performs an independent 2D segmentation task, with the final layer combining the individual outputs in order to produce the final 3D segmentation map.

8.5.2.2. SegNet (Badrinarayanan et al. 2017)

This is a convolutional Encoder-Decoder whose final layer performs classification on a pixel-wise manner. The topology of the Encoder part of the network is borrowed from VGG16 (Simonyan and Zisserman 2015), a state-of-the-art network developed by Google Deep Mind. To effectively deal with the issue of having a massive network but relatively few data for learning, we resort to transfer learning (Tajbakhsh et al. 2016) and initialise the Encoder weights to the values of the corresponding ones of VGG16. Those parameters are kept fixed; thus, the only trainable part of the network is the Decoder and the final softmax layer. SegNet receives the available imaging modalities (DESS,T2,SPGR) as three input channels, effectively allowing us to leverage its architecture for multi-view learning, as each channel accepts a slice of one of the three sequences.

8.5.2.3. Multi-Modal (MM) CNN-FC

This deep learning model is based on the architecture proposed in (Guo et al. 2019), for the segmentation of multi-modal medical images. The authors suggest a simple arrangement of a series of 5 convolutional layers, followed by 3 fully connected ones, and a final softmax layer to perform multi-class classification. They distinguish between three types of fusion architectures, based on the specific location in the network where the fusion takes place. Here, we implement a Type-II fusion network, where initially the three available modalities are fed as different input channels. Feature fusion takes place at the first layer of the fully connected part of the network, meaning that features are learnt from each modality separately by their respective convolutional parts.

8.5.2.4. DenseVoxNet (Yu et al. 2017)

This is a convolutional network proposed for cardiovascular MRI segmentation. It consists of a downsampling and an upsampling part, with skip connections from each layer to all its subsequent ones, thus improving the information flow across the entire network and achieving more robust supervision overall. In our experiments, DenseVoxNet was trained using the same initialisation scheme for the parameters’ values as the one described in the original paper.

8.5.2.5. VoxResNet (Chen et al. 2016)

This is a deep residual network consisting of a series of stacked residual modules (VoxRes modules), with each one performing a series of Batch Normalization and Convolution operations, and containing a skip connection for the module’s input to its output (skip connection). By concatenating the different imaging modalities and feeding them as input to the network, the complementary information they provide is jointly fused in an implicit way.

8.5.2.6. KCB-Net (Peng et al. 2022)

This is a recently proposed network performing cartilage and bone segmentation. An initial representative slice selection step, in all major orientations (coronal, axial, sagittal), is used to single out the most representative ones on which expert annotations will be performed. Three base learners (2D modules) are then trained on the annotated slices and pseudo labels are assigned to the remaining ones, and a 3D module is subsequently trained. Following that, a variant of the DenseVoxNet is used as a meta-learner in order to obtain the full segmentation. Finally, a post-processing step incorporating a fine-tuning network (Xie et al. 2022) is employed to enhance the surface delineation between adjacent cartilage and bone surfaces.

8.5.2.7. CAN3D

This is an extension of the Context Aggregation Network proposed in (Dai et al. 2021). Utilizing a series of progressively dilated convolution operators, it can efficiently aggregate multi-scale information by extracting dense features and progressively expanding its receptive field, allowing for the final voxelwise classification to occur in full resolution. In order to deal with highly imbalanced class data (as in our case), CAN3D employs a hybrid loss function by combining the Dice Similarity Coefficient (DSC) with the Dice Squared Focal Loss (DSF), a variant of DSC with MSC built inside.

8.5.2.8. 3D CNN + SSM

The work presented in (Ambellan et al. 2019) features a segmentation framework that combines 2D and 3D CNNs in conjunction with a Statistical Shape Model (SSM). The overall workflow consists of two branches, where the first one is tasked with the segmentation of the femoral and tibial bone structures, facilitated through a cascade of CNN and SSM steps. The resulting masks are then passed on to the second branch, where another 3D CNN finalises the segmentation by producing the respective maps for the femoral and tibial cartilage structures. All of the above steps are performed separately for both structures.

In all the above cases, the image dataset was split into three disjoint subsets, used for training, validation and testing respectively $(60\mathrm{\%}-20\mathrm{\%}-20\mathrm{\%})$. All the images were pre-processed as described in Section 2.1. The training parameters (learning rate, optimisation algorithm, etc.) and validation schemes of the above models were chosen based on the respective original papers.

8.5.2.9. DAN

(Zhang et al. 2017) The Deep Adversarial Network (DAN) exploits information from both labelled and unlabelled images, by utilising two sub-modules: (1) a segmentation network (SN) tasked with providing annotated label maps and (2) an evaluation network (EN) that continuously ranks the performance achieved by SN. The training process consists of an iterative procedure whereby SN learns to successively output more accurate segmentation maps, guided by the grading its performance on unlabeled images by EN. Essentially, the EN module attempts to”teach” the student SN what an appropriate segmentation map looks like, thus rendering the overall model capable of utilising unlabelled data to obtain enhanced segmentation quality.

8.5.2.10. UA-MT

The Uncertainty-Aware Self-Ensembling model presented in (Yu et al. 2019) is another variant in the teacher-student overall setting. To deal with the possibility of unreliable inputs from the teacher, the uncertainty-aware mean teacher (UA-MT) framework is proposed, where the student model is gradually guided to learn from reliable targets generated by the teacher, whilst the teacher estimates the uncertainty of each generated prediction. This process is exploited as a means of filtering out unreliable outputs from the student, which in turn encourages the teacher to supply higher-quality targets. Regarding the implementation details, the core for both teacher and student is designed around the topology of th V-net (Milletari et al. 2016), by removing the residual connections in each separate convolutional block, and adapting the rest of the model as a Bayesian network, in order to render it capable of estimating the uncertainty.

8.5.2.11. SS-DTC

The authors in (Luo et al. 2021) present a dual-task-consistency semi-supervised framework that can jointly generate a voxel-wise segmentation map and a geometry-aware level-set representation of the target. A differentiable task transform layer is employed to convert said representation into an approximated segmentation map where subsequently, a combined loss function is utilised in order to facilitate the minimisation of discrepancy between the predicted segmentation map and the probability map converted from the level set function. The aforementioned mechanism allows the model to harness unlabelled data in order to boost the overall segmentation performance.

For the three semi-supervised methods described in the paragraphs above, we follow the same setting as in our own proposed methods, regarding the percentage of labelled vs unlabelled MRIs used. The optimal parameters for each model were chosen according to the validation scheme proposed in paragraph 8.3.

9. Experimental results & discussion

9.1. Parameter sensitivity analysis

In this section, we validate the performance effect of the five critical hyperparameters $(n_{\mathcal{A}},\mathit{ss}{l}_{\mathit{perc}},\mathit{k},\mathit{\lambda },\mathit{\mu })$. While varying one parameter, the rest remain fixed at their default values (Section 8.3).

9.1.1. Number of selected atlases n_A

The number of selected atlases $n_{\mathcal{A}}$ controls the size of the atlas library $\mathcal{A}$ which contributes labeled data at both stages of label propagation. During the transductive stage (stage-1), all atlases are sampled in a stratified manner at locations specified by the global target dataset, ${\mathbf{X}}_g^{\mathcal{T}}$, with the number of sampled voxels per atlas being inversely proportional to the number of available atlases, $n_g^{\mathcal{A}}=\mathit{ss}{l}_{\mathit{perc}}\cdot n_g^{\mathcal{T}}/n_{\mathcal{A}}$ (Section 5.1). At the inductive stage (stage-2) however, this relationship changes as the iterative sampling process via the 3D mesh (Section 6) yields a number of voxels directly proportional to $n_{\mathcal{A}}$.

Figure 7 showcases the effect of $n_{\mathcal{A}}\in \{5,10,20\}$ in terms $\mathit{DSC}$ score for the two cartilage classes (FC, TC). In both cases $n_{\mathcal{A}}=10$ yields the best results, while $n_{\mathcal{A}}=5$ and $n_{\mathcal{A}}=20$ induce a slight drop in performance for the FC class. The same effects are observed for the TC class, albeit more pronounced. The concept of bias-variance trade-off can be useful in interpreting this behavior: a small number of atlases implies that more voxels are sampled from each atlas and that labeling errors emanating from potential misalignments between the target and an atlas image, may not be efficiently countered by the remaining atlases. On the other hand, a large $n_{\mathcal{A}}$ suggests that smaller number of voxels is sampled from each atlas, greatly increasing the variance and noise in ${\mathbf{X}}_g^{\mathcal{A}}$. It seems that the optimal strategy is to strike a balance between the two options; choosing enough atlases to cover as much label information content as possible, but not too many as to render the samples noisy.

Figure 7. Effect of number of atlases on $\mathit{DSC}$ score.

9.1.2. Percentage of labelled voxels ssl_perc

The parameter $\mathit{ss}{l}_{\mathit{perc}}$ effectively controls the amount of supervision available to our method. For a fixed number of atlases, a greater $\mathit{ss}{l}_{\mathit{perc}}$ means that more voxels will be drawn from each atlas to form the initial labeled global dataset. We test the effect of this parameter by evaluating models with $\mathit{ss}{l}_{\mathit{perc}}\in \{10\mathrm{\%},20\mathrm{\%},30\mathrm{\%}\}$ (Figure 8).

Figure 8. Effect of amount of supervision on $\mathit{DSC}$ score.

As can be seen, the segmentation performance for both FC and TC is directly proportional to the amount of supervision allowed. As expected, increasing $\mathit{ss}{l}_{\mathit{perc}}$ yields better results, since more class-specific content is available to the learning algorithm. However, the accompanying computational cost rises considerably, suggesting that a more moderate amount of supervision should be considered.

9.1.3. Number of nearest neighbours $\mathit{k}$

The number $\mathit{k}$ of nearest neighbors is a significant parameter controlling the construction of the sparse graphs in both stages of the label propagation. Figure 9 shows the variation of $\mathit{DSC}$ score for neighborhood sizes $\mathit{k}\in \{10,25,50\}$.

Figure 9. Effect of number of sparse neighbors on $\mathit{DSC}$ score.

It seems that $\mathit{k}$ yields better results when assuming values at the lower range, with the best results obtained for $\mathit{k}=10$, while $\mathit{k}=50$ performs markedly worse. This result can be also explained by considering the bias-variance trade-off. A large neighbourhood (high variance) allows each voxel to be connected with neighbours sharing similar spectral presentation, which might though belong to a different class (i.e. Femoral & Tibial cartilage). Limiting $\mathit{k}$ forces graph construction to draw edges between only the most similar of voxels. Contrary to the selection of $n_{\mathcal{A}}$, here the choice of incurring a high bias (small $\mathit{k}$) seems to work best. A plausible explanation lies on the multi-view learning mechanism; each view of a specific neighbourhood may be characterised by a high bias, but their combination seems to have a regularisation effect nullifying that bias.

9.1.4. Structural regularisation parameter $\mathit{\lambda }$

Parameter $\mathit{\lambda }$ controls the amount of structural regularization in AMUSE (4). In the extreme case where $\mathit{\lambda }=0$, the problem degenerates to the single-view case, while for $\mathit{\lambda }\to \mathrm{\infty }$, the objective is trivialized to computing the combined multi-view graph as a linear combination of the single-view ones. Figure 10, highlights its effect on the $\mathit{DSC}$ score for both cartilage classes, when assuming values in the logarithmic scale $\mathit{\lambda }\in \{0.01,1,100\}$.

Figure 10. Effect of structural regularization parameter on $\mathit{DSC}$ score.

It is evident that increasing the value of $\mathit{\lambda }$ yields better segmentation accuracy for both FC and TC, up to a certain point, while performance drops markedly for extremely small values. Interestingly, the segmentation accuracy on FC seems to diminish for the larger values of $\mathit{\lambda }$. Since the increase of $\mathit{\lambda }$ strengthens the contribution of the combined graph learning, it seems that the multi-view component of MV-KCS is essential to the overall performance. Care must be taken however, since past a certain point the objective (4) is dominated by the second term and the solution is trivialised.

9.1.5. Convex combination parameter $\mathit{\mu }$

Finally, we validate the effect of the convex combination parameter $\mathit{\mu }$ (27), which controls the relative contribution of local and global datasets ${\mathbf{X}}_l,{\mathbf{X}}_g$ respectively (stage-2).

As can be seen from Figure 11, the performance follows a somewhat symmetrical pattern, achieving the best results for both classes in the $[0.4,0.6]$ range. This suggests that a balance between local and global information is the optimal strategy. Setting $\mathit{\mu }$ to either $\mathit{\mu }=1$ or $\mathit{\mu }=0$, entails that the label transfer mechanism relies exclusively on either global or local information, respectively. In the former case, any potential errors in the initial labeling of the voxels in ${\mathbf{X}}_g^{\mathcal{T}}$ will keep propagating to the out-of-sample batches, without the regularizing effect of ${\mathbf{X}}_l^{\mathcal{A}}$. Further, in the latter case, the rich label content of ${\mathbf{X}}_g^{\mathcal{A}}$ is completely discarded.

Figure 11. Effect of convex combination parameter on $\mathit{DSC}$ score.

9.2. Test cases – results

Table 3 provides a thorough evaluation of our proposed labelling method MV-HyLP under all testing cases outlined previously. The reported results refer to both Cartilage classes (FC, TC), while the performance quality is assessed in terms of Recall and Precision, as well as the volumetric indices DSC, VD, and VOE.

Table 3. Performance comparison (mean and std. deviations) of all test cases for the proposed method of MV-HyLP. In case #1, GLM denotes the combination of GLCM, GLLM and All, the inclusion of all feature descriptors. Stacked and fused in case #2 refer to the list of feature descriptors (HOG, LBP, HCGF) being concatenated and fused through AMUSE, respectively. Finally, in case #3, Multi-feature corresponds exactly to the last row in case #1.

		Femoral Cartilage					Tibial Cartilage
Method	Case	Recall	Precision	DSC	VOE	VD	Recall	Precision	DSC	VOE	VD
Group #1	HOG	$88.81\mathrm{\%}$	$88.17\mathrm{\%}$	$88.56\mathrm{\%}$	$23.19\mathrm{\%}$	$6.02\mathrm{\%}$	$85.72\mathrm{\%}$	$83.91\mathrm{\%}$	$84.03\mathrm{\%}$	$27.13\mathrm{\%}$	$9.09\mathrm{\%}$
		$(\pm 0.212)$	$(\pm 0.227)$	$(\pm 0.191)$	$(\pm 0.208)$	$(\pm 0.265)$	$(\pm 0.176)$	$(\pm 0.134)$	$(\pm 0.162)$	$(\pm 0.207)$	$(\pm 0.181)$
	HOG+LBP	$91.97\mathrm{\%}$	$90.40\mathrm{\%}$	$89.67\mathrm{\%}$	$22.35\mathrm{\%}$	$6.02\mathrm{\%}$	$85.88\mathrm{\%}$	$85.36\mathrm{\%}$	$85.79\mathrm{\%}$	$25.29\mathrm{\%}$	$7.79\mathrm{\%}$
		$(\pm 0.055)$	$(\pm 0.049)$	$(\pm 0.029)$	$(\pm 0.035)$	$(\pm 0.038)$	$(\pm 0.046)$	$(\pm 0.27)$	$(\pm 0.021)$	$(\pm 0.033)$	$(\pm 0.016)$
	HOG+LBP+GLM	$92.24\mathrm{\%}$	$89.76\mathrm{\%}$	$90.37\mathrm{\%}$	$21.34\mathrm{\%}$	$5.84\mathrm{\%}$	$86.89\mathrm{\%}$	$87.75\mathrm{\%}$	$86.54\mathrm{\%}$	$24.54\mathrm{\%}$	$6.67\mathrm{\%}$
		$(\pm 0.052)$	$(\pm 0.051)$	$(\pm 0.026)$	$(\pm 0.028)$	$(\pm 0.034)$	$(\pm 0.047)$	$(\pm 0.023)$	$(\pm 0.018)$	$(\pm 0.025)$	$(\pm 0.011)$
	All	$93.24\mathrm{\%}$	$91.26\mathrm{\%}$	$91.76\mathrm{\%}$	$20.65\mathrm{\%}$	$5.17\mathrm{\%}$	$90.89\mathrm{\%}$	$89.79\mathrm{\%}$	$88.65\mathrm{\%}$	$22.95\mathrm{\%}$	$6.05\mathrm{\%}$
		$(\pm 0.052)$	$(\pm 0.051)$	$(\pm 0.026)$	$(\pm 0.028)$	$(\pm 0.034)$	$(\pm 0.047)$	$(\pm 0.023)$	$(\pm 0.018)$	$(\pm 0.025)$	$(\pm 0.011)$
Group #2	Stacked	$88.70\mathrm{\%}$	$90.10\mathrm{\%}$	$89.56\mathrm{\%}$	$21.95\mathrm{\%}$	$55.90\mathrm{\%}$	$76.12\mathrm{\%}$	$80.90\mathrm{\%}$	$85.43\mathrm{\%}$	$5.16\mathrm{\%}$	$29.05\mathrm{\%}$
		$(\pm 0.195)$	$(\pm 0.214)$	$(\pm 0.158)$	$(\pm 0.223)$	$(\pm 0.144)$	$(\pm 0.098)$	$(\pm 0.117)$	$(\pm 0.154)$	$(\pm 0.174)$	$(\pm 0.139)$
	Fused	$77.90\mathrm{\%}$	$87.30\mathrm{\%}$	$90.67\mathrm{\%}$	$12.00\mathrm{\%}$	$29.97\mathrm{\%}$	$91.50\mathrm{\%}$	$93.60\mathrm{\%}$	$87.12\mathrm{\%}$	$86.20\mathrm{\%}$	$11.18\mathrm{\%}$
		$(\pm 0.098)$	$(\pm 0.089)$	$(\pm 0.045)$	$(\pm 0.038)$	$(\pm 0.063)$	$(\pm 0.038)$	$(\pm 0.041)$	$(\pm 0.057)$	$(\pm 0.061)$	$(\pm 0.044)$
Group #3	Multi-Modal	$90.02\mathrm{\%}$	$90.71\mathrm{\%}$	$89.93\mathrm{\%}$	$18.96\mathrm{\%}$	$5.85\mathrm{\%}$	$85.49\mathrm{\%}$	$87.59\mathrm{\%}$	$87.31\mathrm{\%}$	$8.78\mathrm{\%}$	$21.76\mathrm{\%}$
		$(\pm 0.232)$	$(\pm 0.247)$	$(\pm 0.191)$	$(\pm 0.218)$	$(\pm 0.241)$	$(\pm 0.168)$	$(\pm 0.137)$	$(\pm 0.185)$	$(\pm 0.202)$	$(\pm 0.153)$
	All	$94.04\mathrm{\%}$	$93.16\mathrm{\%}$	$92.56\mathrm{\%}$	$15.16\mathrm{\%}$	$5.12\mathrm{\%}$	$91.08\mathrm{\%}$	$89.98\mathrm{\%}$	$89.91\mathrm{\%}$	$19.67\mathrm{\%}$	$5.85\mathrm{\%}$
		$(\pm 0.052)$	$(\pm 0.051)$	$(\pm 0.026)$	$(\pm 0.028)$	$(\pm 0.034)$	$(\pm 0.047)$	$(\pm 0.023)$	$(\pm 0.018)$	$(\pm 0.025)$	$(\pm 0.011)$

9.2.1. Case #1: single-view vs multi-view

This case highlights the effects in performance when incorporating progressively additional feature descriptors in MV-KCS. The first model corresponds to the trivial single-view case, utilising solely the HOG descriptor. HOGs are very powerful in encoding information content pertaining to shape and orientation. Since the anatomy of the cartilage has a very specific geometry, being locally approximated by a thin 3D plate, HOG manages to yield a satisfactory segmentation performance $(\mathit{DSC}=(88.56\mathrm{\%},84.03\mathrm{\%}))$. Although distinctive, HOG is not the best choice in capturing texture information, a salient feature in biomedical images.

This limitation is confronted by incorporating LBP descriptors as an additional view. LBPs are excellent in capturing texture information, making them an ideal synergist to their HOG counterparts. This point is succintly reflected in the improved segmentation score achieved $(\mathit{DSC}=(89.67\mathrm{\%},85.79\mathrm{\%}))$. Under the multi-view setting, each voxel now participates in two single-view sparse graphs, with its spectral neighborhood in each case reflecting a difference notion of similarity. Voxels of different classes exhibiting similar HOG descriptors may assume substantially different characterizations under the LBP view, and vice versa. In such a case, the affinity matrix $\mathbf{S}={\alpha }_{(\mathit{HOG})}{\mathbf{W}}^{(\mathit{HOG})}\text{break}+{\alpha }_{(\mathit{LBP})}{\mathbf{W}}^{(\mathit{LBP})}$ will not be incentivized to maintain each one in the neighborhood of the other, and a potential source of erroneous label propagation can be effectively eliminated. On the other hand, voxels with high affinity across all views, are more likely to share the same labels. This assumption is reflected in the construction of the combined graph, reinforcing propagation of label content from one to the other.

The above discussion remains intact when adding two, and finally three more descriptors (views). Incorporating together the GLCM and GLRLM features, we observe a slight increase in performance $(\mathit{DSC}=(90.37\mathrm{\%},86.54\mathrm{\%}))$ for both classes of interest. Since GLCM and GLRLM mainly capture texture, as do LBP, and are variants of a similar underlying encoding, the additional views are not distinctive enough to provide a remarkable improvement. On the other hand, the inclusion of HCGF yields a more noticeable effect $(\mathit{DSC}=(91.76\mathrm{\%},88.65\mathrm{\%}))$. This effect can be attributed to the structure of the HCGF descriptor, which incorporates multiple sources of low-order statistics (location, intensity, first & second order derivatives etc.) in its encoding. Thus, it is allowed to offer a substantially different voxel characterisation relative to the preceding descriptors. Figure 12 showcases the gradual improvements in segmentation accuracy with the progressive incorporation of additional feature descriptors, for a specific subject.

Figure 12. Segmentation results for femoral (FC) and tibial (TC) cartilage for the four instances of case #1 (left to right: Ground Truth, HOG, HOG & LBP, HOG & LBP & GLCM, all features). The first part of the figure illustrates a case of successful application of MV-KCS, while the second part presents in instance characterized by poor performance. In both cases, however, the positive effect of incorporating additional views is noticeable.

The results seem to indicate that multi-view learning is a powerful tool, significantly boosting the performance. Nevertheless, care must be taken so that the selected views are complementary to each other.

Finally, to establish that the reported results are statistically different, the Friedman statistical test (Friedman 1937) is performed. The test yields a ${\chi }^2$ value of 312.6 with 3 degrees of freedom, which corresponds to a $\mathit{p}$-value $6.23\cdot {10}^{-09}$ rejecting the null hypothesis of statistical equivalence. Following that, a post-hoc Nemenyi test (Nemenyi 1963) is carried out, to evaluate each pair individually. The results in Table 4 $(\mathit{\alpha }=0.05)$, suggest that statistically, the inclusion of additional views further boosts segmentation performance.

Table 4. Pairwise comparisons using the Nemenyi test for the cartilage $\mathit{DSC}$ results of the test cases in Table 3. Reported are the corresponding $\mathit{p}$-values. The binary indices in the subscript of each method’s name MV-HyLP_xxxxx indicate the inclusion or exclusion of the respective feature descriptor HOG, LBP, GLCM, GLRLM, HCGF..

	MV-HyLP10000	MV-HyLP11000	MV-HyLP11110
MV-HyLP11000	$2.81\times {10}^{-03}$	-	-
MV-HyLP11110	$1.05\times {10}^{-05}$	$3.75\times {10}^{-03}$	-
MV-HyLP11111	$4.28\times {10}^{-05}$	$2.25\times {10}^{-04}$	$1.13\times {10}^{-02}$

9.2.2. Case #2: feature stacking vs multi-view graph learning

The purpose of this experiment is to evaluate the two different apporaches to combine feature descriptors: the multi-view learning implemented by in MV-KCS, as contrasted to the simple stacking of multiple features. To keep the dimensionality of the stacked feature vector at a reasonable scale, the experiments consider the least redundant feature subset, (HOG, LBP, HCGF). The results are reported in the second part of Table 3. It is clear that combining the views in a more structured manner yields superior performance $(\mathit{DSC}=(90.67\mathrm{\%},87.12\mathrm{\%}))$, for both FC and TC classes. This is due to the fact that simple stacking of feature vectors most likely introduces considerable noise, while also distorting the similarity between data. This last notice hinges heavily on the size of the feature vectors: in a high dimensional space, there is a great number of ways to define similarity between voxels, and only a very small subset of those correspond to label similarity. Therefore, most likely a learning algorithm identifies the ’wrong’ type of similarity, adversely affecting performance $(\mathit{DSC}=(89.65\mathrm{\%},85.43\mathrm{\%}))$. By combining multiple views via MV-KCS, through structural regularization on a combined graph, we avoid the dimensionality issue and allow for a more natural integration of the available information content supplied by each view. The superiority of the multi-view setting against feature stacking is statistically verified by performing the Wilcoxon signed rank (Wilcoxon 1946) test. The obtained $\mathit{p}$-value $(\mathit{p}=4.8\cdot {10}^{-03})$ confirms this conclusion.

9.2.3. Case #3: multiple features vs multiple modalities

This final test case serves to contrast the two different multi-view settings against each other. Specifically, we examine the merit of integrating multiple feature descriptors compared to utilising multiple imaging sequences (results in third part of Table 3). The Multi-feature case exhibits the best performance $(\mathit{DSC}=(92.56\mathrm{\%},89.91\mathrm{\%}))$, exhibiting a noticeable overhead in comparison to the Multi-modal $(\mathit{DSC}=(89.93\mathrm{\%},87.31\mathrm{\%}))$. Generating multiple data views via distinctive encodings seems to be the optimal strategy. On the contrary, multiple modalities, using a single spectral descriptor yield redundant representations, thus limiting their ability to provide an effective information fusion.

9.3. Comparative results

The comparative results in Table 5 highlight that both our methods MV-RegLP and MV-HyLP report superior results against the competing methodologies, across all examined measures. This effect is observed in both classes, with FC showcasing a slightly augmented performance boost compared to TC.

Table 5. Summary of segmentation performance measures (means and std. deviations) on the two cartilage classes of our proposed method MV-RegLP and its modification MV-HyLP, compared to Patch-Based Sparse-Coding (PB_SC), Patch-Based Sparse-Coding with stacked features (PB_SC^(stacked)), Patch-Based joint label fusion (PB-JLF), Triplanar CNN, SegNet, MM-CNN-FC, VoxResNet, DenseVoxNet, KCB-Net, CAN3D, 3D-CNN + SSM, DAN, UA-MT and SS-DTC. The hyperparameter values for MV-RegLP and MV-HyLP are those obtained after parameter sensitivity analysis in previous section 8.3.

	Femoral Cartilage					Tibial Cartilage
Method	Recall	Precision	$\mathit{DSC}$	$\mathit{VOE}$	$\mathit{VD}$	Recall	Precision	$\mathit{DSC}$	$\mathit{VOE}$	$\mathit{VD}$
PBSC	$83.51\mathrm{\%}$	$82.65\mathrm{\%}$	$82.23\mathrm{\%}$	$29.97\mathrm{\%}$	$11.01\mathrm{\%}$	$79.76\mathrm{\%}$	$81.45\mathrm{\%}$	$78.85\mathrm{\%}$	$34.28\mathrm{\%}$	$11.79\mathrm{\%}$
	$(\pm 0.065)$	$(\pm 0.042)$	$(\pm 0.047)$	$(\pm 0.062)$	$(\pm 0.074)$	$(\pm 0.033)$	$(\pm 0.022)$	$(\pm 0.021)$	$(\pm 0.044)$	$(\pm 0.043)$
PBSC^(stacked)	$84.12\mathrm{\%}$	$84.67\mathrm{\%}$	$84.35\mathrm{\%}$	$26.96\mathrm{\%}$	$7.35\mathrm{\%}$	$82.84\mathrm{\%}$	$83.67\mathrm{\%}$	$81.91\mathrm{\%}$	$30.12\mathrm{\%}$	$10.92\mathrm{\%}$
	$(\pm 0.068)$	$(\pm 0.041)$	$(\pm 0.044)$	$(\pm 0.062)$	$(\pm 0.078)$	$(\pm 0.039)$	$(\pm 0.016)$	$(\pm 0.027)$	$(\pm 0.033)$	$(\pm 0.034)$
PB-JLF	$84.61\mathrm{\%}$	$85.03\mathrm{\%}$	$84.92\mathrm{\%}$	$26.21\mathrm{\%}$	$7.28\mathrm{\%}$	$83.15\mathrm{\%}$	$83.92\mathrm{\%}$	$82.37\mathrm{\%}$	$29.41\mathrm{\%}$	$9.84\mathrm{\%}$
	$(\pm 0.102)$	$(\pm 0.065)$	$(\pm 0.081)$	$(\pm 0.104)$	$(\pm 0.092)$	$(\pm 0.061)$	$(\pm 0.048)$	$(\pm 0.036)$	$(\pm 0.067)$	$(\pm 0.041)$
Triplanar CNN	$85.04\mathrm{\%}$	$85.27\mathrm{\%}$	$84.98\mathrm{\%}$	$25.77\mathrm{\%}$	$7.11\mathrm{\%}$	$83.44\mathrm{\%}$	$84.07\mathrm{\%}$	$83.19\mathrm{\%}$	$28.96\mathrm{\%}$	$8.41\mathrm{\%}$
	$(\pm 0.142)$	$(\pm 0.205)$	$(\pm 0.071)$	$(\pm 0.055)$	$(\pm 0.105)$	$(\pm 0.082)$	$(\pm 0.035)$	$(\pm 0.049)$	$(\pm 0.033)$	$(\pm 0.040)$
MM-CNN-FC	$87.18\mathrm{\%}$	$87.48\mathrm{\%}$	$85.09\mathrm{\%}$	$23.73\mathrm{\%}$	$6.65\mathrm{\%}$	$85.22\mathrm{\%}$	$85.07\mathrm{\%}$	$84.12\mathrm{\%}$	$26.79\mathrm{\%}$	$8.01\mathrm{\%}$
	$(\pm 0.134)$	$(\pm 0.259)$	$(\pm 0.104)$	$(\pm 0.076)$	$(\pm 0.088)$	$(\pm 0.067)$	$(\pm 0.072)$	$(\pm 0.058)$	$(\pm 0.031)$	$(\pm 0.029)$
VoxResNet	$87.88\mathrm{\%}$	$88.67\mathrm{\%}$	$85.81\mathrm{\%}$	$23.05\mathrm{\%}$	$6.56\mathrm{\%}$	$86.76\mathrm{\%}$	$85.91\mathrm{\%}$	$84.52\mathrm{\%}$	$26.13\mathrm{\%}$	$7.89\mathrm{\%}$
	$(\pm 0.186)$	$(\pm 0.221)$	$(\pm 0.078)$	$(\pm 0.065)$	$(\pm 0.104)$	$(\pm 0.078)$	$(\pm 0.043)$	$(\pm 0.039)$	$(\pm 0.028)$	$(\pm 0.036)$
DenseVoxNet	$88.04\mathrm{\%}$	$88.74\mathrm{\%}$	$87.12\mathrm{\%}$	$22.43\mathrm{\%}$	$6.49\mathrm{\%}$	$87.18\mathrm{\%}$	$86.26\mathrm{\%}$	$85.05\mathrm{\%}$	$25.76\mathrm{\%}$	$7.43\mathrm{\%}$
	$(\pm 0.156)$	$(\pm 0.205)$	$(\pm 0.03)$	$(\pm 0.042)$	$(\pm 0.098)$	$(\pm 0.055)$	$(\pm 0.029)$	$(\pm 0.033)$	$(\pm 0.017)$	$(\pm 0.021)$
SegNet	$88.24\mathrm{\%}$	$88.91\mathrm{\%}$	$87.22\mathrm{\%}$	$22.14\mathrm{\%}$	$6.44\mathrm{\%}$	$87.51\mathrm{\%}$	$86.84\mathrm{\%}$	$85.71\mathrm{\%}$	$25.23\mathrm{\%}$	$7.28\mathrm{\%}$
	$(\pm 0.156)$	$(\pm 0.204)$	$(\pm 0.042)$	$(\pm 0.048)$	$(\pm 0.121)$	$(\pm 0.076)$	$(\pm 0.041)$	$(\pm 0.047)$	$(\pm 0.028)$	$(\pm 0.025)$
CAN3D	$88.75\mathrm{\%}$	$89.67\mathrm{\%}$	$88.32\mathrm{\%}$	$21.83\mathrm{\%}$	$6.22\mathrm{\%}$	$87.85\mathrm{\%}$	$87.13\mathrm{\%}$	$86.68\mathrm{\%}$	$24.47\mathrm{\%}$	$7.13\mathrm{\%}$
	$(\pm 0.187)$	$(\pm 0.205)$	$(\pm 0.047)$	$(\pm 0.055)$	$(\pm 0.128)$	$(\pm 0.071)$	$(\pm 0.044)$	$(\pm 0.043)$	$(\pm 0.032)$	$(\pm 0.029)$
KCB-Net	$89.74\mathrm{\%}$	$90.12\mathrm{\%}$	$88.92\mathrm{\%}$	$20.13\mathrm{\%}$	$6.15\mathrm{\%}$	$88.12\mathrm{\%}$	$87.46\mathrm{\%}$	$87.96\mathrm{\%}$	$23.90\mathrm{\%}$	$7.04\mathrm{\%}$
	$(\pm 0.149)$	$(\pm 0.185)$	$(\pm 0.031)$	$(\pm 0.042)$	$(\pm 0.098)$	$(\pm 0.055)$	$(\pm 0.029)$	$(\pm 0.033)$	$(\pm 0.017)$	$(\pm 0.021)$
3D-CNN + SSM	$90.23\mathrm{\%}$	$89.91\mathrm{\%}$	$90.18\mathrm{\%}$	$21.69\mathrm{\%}$	$5.91\mathrm{\%}$	$89.31\mathrm{\%}$	$88.25\mathrm{\%}$	$88.74\mathrm{\%}$	$21.92\mathrm{\%}$	$6.19\mathrm{\%}$
	$(\pm 0.104)$	$(\pm 0.201)$	$(\pm 0.065)$	$(\pm 0.081)$	$(\pm 0.109)$	$(\pm 0.061)$	$(\pm 0.048)$	$(\pm 0.045)$	$(\pm 0.051)$	$(\pm 0.063)$
DAN	$89.22\mathrm{\%}$	$90.07\mathrm{\%}$	$89.65\mathrm{\%}$	$22.87\mathrm{\%}$	$6.41\mathrm{\%}$	$88.17\mathrm{\%}$	$86.02\mathrm{\%}$	$87.61\mathrm{\%}$	$22.11\mathrm{\%}$	$6.91\mathrm{\%}$
	$(\pm 0.158)$	$(\pm 0.139)$	$(\pm 0.035)$	$(\pm 0.039)$	$(\pm 0.094)$	$(\pm 0.051)$	$(\pm 0.041)$	$(\pm 0.031)$	$(\pm 0.019)$	$(\pm 0.028)$
UA-MT	$90.14\mathrm{\%}$	$90.71\mathrm{\%}$	$89.96\mathrm{\%}$	$22.08\mathrm{\%}$	$6.12\mathrm{\%}$	$89.05\mathrm{\%}$	$87.96\mathrm{\%}$	$88.13\mathrm{\%}$	$21.87\mathrm{\%}$	$6.25\mathrm{\%}$
	$(\pm 0.29)$	$(\pm 0.21)$	$(\pm 0.041)$	$(\pm 0.029)$	$(\pm 0.103)$	$(\pm 0.063)$	$(\pm 0.037)$	$(\pm 0.024)$	$(\pm 0.022)$	$(\pm 0.031)$
SS-DTC	$91.45\mathrm{\%}$	$92.72\mathrm{\%}$	$91.07\mathrm{\%}$	$18.01\mathrm{\%}$	$5.34\mathrm{\%}$	$90.31\mathrm{\%}$	$89.15\mathrm{\%}$	$89.42\mathrm{\%}$	$22.13\mathrm{\%}$	$5.94\mathrm{\%}$
	$(\pm 0.25)$	$(\pm 0.32)$	$(\pm 0.028)$	$(\pm 0.051)$	$(\pm 0.089)$	$(\pm 0.076)$	$(\pm 0.041)$	$(\pm 0.039)$	$(\pm 0.028)$	$(\pm 0.032)$
MV-RegLP	$91.31\mathrm{\%}$	$92.16\mathrm{\%}$	$90.23\mathrm{\%}$	$18.21\mathrm{\%}$	$5.73\mathrm{\%}$	$89.12\mathrm{\%}$	$88.19\mathrm{\%}$	$88.82\mathrm{\%}$	$22.82\mathrm{\%}$	$6.08\mathrm{\%}$
	$(\pm 0.091)$	$(\pm 0.103)$	$(\pm 0.056)$	$(\pm 0.071)$	$(\pm 0.083)$	$(\pm 0.071)$	$(\pm 0.047)$	$(\pm 0.061)$	$(\pm 0.054)$	$(\pm 0.037)$
MV-HyLP	$94.04\mathrm{\%}$	$93.16\mathrm{\%}$	$92.56\mathrm{\%}$	$15.16\mathrm{\%}$	$5.12\mathrm{\%}$	$91.08\mathrm{\%}$	$89.98\mathrm{\%}$	$89.91\mathrm{\%}$	$19.67\mathrm{\%}$	$5.85\mathrm{\%}$
	$(\pm 0.052)$	$(\pm 0.051)$	$(\pm 0.026)$	$(\pm 0.028)$	$(\pm 0.034)$	$(\pm 0.047)$	$(\pm 0.023)$	$(\pm 0.018)$	$(\pm 0.025)$	$(\pm 0.011)$

The baseline PBSC method achieves the lowest overall score across all competitors $(\mathit{DSC}=(82.23\mathrm{\%},78.85\mathrm{\%}))$. However, its variant that incorporates the stacked feature descriptors offers a noticeable improvement in accuracy $(\mathit{DSC}=(84.35\mathrm{\%},\text{break}81.91\mathrm{\%}))$. Both methods build an extensive patch library ${\mathbf{P}}_{\mathcal{L}}$ for each target voxel by exhaustively extracting patches from its neighborhood. However, feature vectors comprising the patch library assume substantially different characterizations in each case: PB_SC forms the feature vector by simply stacking the voxels’ intensities, while PB_SC^(stacked), leverages three powerful descriptors HOG, LBP, HCFG (Section 3.2.3). The incorporation of information relevant to texture, shape and geometry via the above stacked features demonstrably benefits PB_SC^(stacked) allowing it to more efficiently distinguish between similarly appearing classes (Femoral & Tibial Cartilage, Cartilages & Background). Finally, the PB-JLF method seems to obtain the best overall performance among the three patch-based methods considered here. The incorporation of of an extra step of deformable registration on top of the affine one, proves to be beneficial with respect to the accuracy of the resulting segmentation maps.

The eight deep learning models exhibit relatively comparable results, with all of them achieving superior outcomes in both femoral and tibial cartilage segmentation with respect to the two baseline patch-based methods. Overall, the 3D-CNN + SSM demonstrates better performance across most measures $(\mathit{DSC}=(90.18\mathrm{\%},88.74\mathrm{\%}))$, showcasing the efficacy of combining a fully 3D convolutional network with a sophisticated post-processing regularization step, in order to finalize the segmentation maps. KCB-Net comes second, demonstrating a state-of-the-art segmentation accuracy $(\mathit{DSC}=(88.92\mathrm{\%},87.96\mathrm{\%}))$ owning to its refined training strategy incorporating both unsupervised and subsequent supervised learning steps. CAN3D reports a slightly lower performance than KCB-Net w.r.t $\mathit{DSC}$ scores $(\mathit{DSC}=(88.32\mathrm{\%},86.68\mathrm{\%}))$, showcasing the effectiveness of the Context Aggregation Network framework (CAN). SegNet $(\mathit{DSC}=(87.22\mathrm{\%},85.71\mathrm{\%}))$ ranks just below CAN3D in both femoral & tibial cartilage. DenseVoxNet $(\mathit{DSC}=(87.12\mathrm{\%},\text{break}85.05\mathrm{\%}))$ and VoxResNet $(\mathit{DSC}=(85.81\mathrm{\%},84.52\mathrm{\%}))$, via the utilization of similar strategies employing skip connections, achieve comparable results in cartilage components, respectively, albeit lower than SegNet, CAN3D and especially KCB-Net. MM-CNN-FC $(\mathit{DSC}=(85.09\mathrm{\%},84.12\mathrm{\%}))$ achieves the second lowest performance with respect to the previous deep models, albeit still superior to the patch-based ones. Finally, the Triplanar CNN method showcases the lowest overall performance with regards to the deep learning methods $(\mathit{DSC}=(84.98\mathrm{\%},83.19\mathrm{\%}))$, demonstrating the fact that decoupling the spatial cohesion of the input image patches by splitting them in three separate planes is a sub-optimal choice.

The next three methods models, utilising the semi-supervised setting in some capacity, demonstrate superior results both with regards to the patch-based methods, as well as compared to the standard deep learning ones. Specifically, DAN $(\mathit{DSC}=(88.65\mathrm{\%},87.61\mathrm{\%}))$ reports superior performance to all previous methods with the exception of KCB-Net, while UA-MT $(\mathit{DSC}=89.61\mathrm{\%},88.13\mathrm{\%})$ and SS-DTC $(\mathit{DSC}=91.07\mathrm{\%},\text{break}89.42\mathrm{\%})$ achieve better overall segmentation results against all previous reported methods. The performance of these approaches serves to strengthen the argument for semi-supervised case, demonstrating the effectiveness of incorporating additional unlabelled data into the learning process.

Regarding our proposed labelling methods, MV-RegLP $(\mathit{DSC}=(90.23\mathrm{\%},88.82\mathrm{\%}))$ outperforms both patch-based methods by a significant margin, while exhibiting slightly better, or similar results to those achieved by the deep learning-based ones (supervised or semi-supervised). The combined synergy of multi-view integration, semi-supervised learning, representation through sparse neighborhoods and aggregation of global and local information, proves to be an effective tool. MV-HyLP pushes the performance even further $(\mathit{DSC}=(92.56\mathrm{\%},89.91\mathrm{\%}))$, by introducing the consideration of an additional spatial scale in the label transfer mechanism, through the concept of LMEs, demonstrably outperforming the rest of the competing methods.

To statistically validate the observed differences across the competing methods, the Friedman statistical test is initially performed. The reported ${\chi }^2$ value with 5 degrees of freedom is $441.81$, corresponding to a $\mathit{p}$-value of less than $2.2\cdot {10}^{-16}$, rejecting the null hypothesis of statistical equivalence among the competitors. This allows for the execution of Nemenyi’s pairwise post-hoc test, to identify the statistically superior method. The reported $\mathit{p}$-values in Table 6 suggest that the performance discrepancy in PB_SC and PB_SC^(stacked) is statistically significant, highlighting the positive impact of combining multiple descriptors instead of relying on raw voxels’ intensities alone, even in the trivial case of stacking the together. Regarding the eight supervised deep learning methods employed, all are proven to be statistically superior to the patch-based ones. More notably, all five methods evaluated in this study, that in some capacity fall within the semi-supervised learning paradigm, are demonstrated to achieve superior performance with regards to their fully supervised counterparts, with the sole exception being the 3D-CNN + SSM method. Concretely, the two methods employing the teacher-student framework (DAN, UA-MT) achieve statistically superior results compared to the standard patch-based and deep learning supervised approaches, showcasing its effectiveness in this particular application. Our own MV-RegLP manages to achieve better overall results, with SS-DTC and 3D-CNN + SSM, in turn, superseding it. Finally, our proposed MV-HyLP is shown to compare favourably against all previous methods reported.

Table 6. Pairwise comparisons using the Nemenyi test for the cartilage $\mathit{DSC}$ results of the competing methodologies of Table 5. Reported are the corresponding $\mathit{p}$-values.

	PBSC	PBSC^(stacked)	PB-JLF	Triplanar CNN	MM-CNN-FC	VoxResNet	DenseVoxNet	SegNet	CAN3D	KCB-Net	DAN	UA-MT	3D-CNN+SSM	MV-RegLP	SS-DTC
PBSC^(stacked)	$4.86\times {10}^{-03}$	-	-	-	-	-	-	-	-	-	-	-	-	-	-
PB-JLF	$6.19\times {10}^{-04}$	$3.81\times {10}^{-03}$	-	-	-	-	-	-	-	-	-	-	-	-	-
Triplanar CNN	$1.92\times {10}^{-05}$	$7.04\times {10}^{-04}$	$3.19\times {10}^{-02}$	-	-	-	-	-	-	-	-	-	-	-	-
MM-CNN-FC	$7.64\times {10}^{-09}$	$3.88\times {10}^{-05}$	$6.81\times {10}^{-04}$	$2.95\times {10}^{-02}$	-	-	-	-	-	-	-	-	-	-	-
VoxResNet	$5.66\times {10}^{-12}$	$7.60\times {10}^{-06}$	$1.75\times {10}^{-05}$	$6.11\times {10}^{-04}$	$3.34\times {10}^{-03}$	-	-	-	-	-	-	-	-	-	-
DenseVoxNet	$1.26\times {10}^{-13}$	$8.32\times {10}^{-05}$	$1.92\times {10}^{-05}$	$6.62\times {10}^{-04}$	$1.24\times {10}^{-03}$	$9.32\times {10}^{-02}$	-	-	-	-	-	-	-	-	-
SegNet	$8.46\times {10}^{-14}$	$6.91\times {10}^{-08}$	$5.72\times {10}^{-07}$	$1.19\times {10}^{-05}$	$1.19\times {10}^{-05}$	$7.36\times {10}^{-02}$	$4.71\times {10}^{-02}$	-	-	-	-	-	-	-	-
CAN3D	$<2.2\times {10}^{-16}$	$2.30\times {10}^{-10}$	$9.61\times {10}^{-08}$	$5.28\times {10}^{-07}$	$4.29\times {10}^{-07}$	$3.91\times {10}^{-04}$	$6.29\times {10}^{-03}$	$1.78\times {10}^{-02}$	-	-	-	-	-	-	-
KCB-Net	$<2.2\times {10}^{-16}$	$9.60\times {10}^{-09}$	$4.94\times {10}^{-09}$	$2.61\times {10}^{-08}$	$7.59\times {10}^{-07}$	$8.22\times {10}^{-06}$	$5.73\times {10}^{-03}$	$4.29\times {10}^{-02}$	$2.78\times {10}^{-01}$	-	-	-	-	-	-
DAN	$<2.2\times {10}^{-16}$	$8.89\times {10}^{-08}$	$3.61\times {10}^{-08}$	$8.19\times {10}^{-07}$	$6.14\times {10}^{-07}$	$9.21\times {10}^{-06}$	$4.71\times {10}^{-05}$	$2.58\times {10}^{-03}$	$7.62\times {10}^{-02}$	$3.41\times {10}^{-02}$	-	-	-	-	-
UA-MT	$<2.2\times {10}^{-16}$	$4.21\times {10}^{-11}$	$1.61\times {10}^{-11}$	$7.04\times {10}^{-09}$	$8.51\times {10}^{-09}$	$1.03\times {10}^{-07}$	$7.91\times {10}^{-07}$	$2.61\times {10}^{-05}$	$5.81\times {10}^{-03}$	$1.22\times {10}^{-03}$	$4.91\times {10}^{-02}$	-	-	-	-
MV-RegLP	$<2.2\times {10}^{-16}$	$8.10\times {10}^{-12}$	$3.90\times {10}^{-12}$	$7.41\times {10}^{-10}$	$4.82\times {10}^{-09}$	$1.71\times {10}^{-07}$	$4.23\times {10}^{-05}$	$1.08\times {10}^{-05}$	$4.77\times {10}^{-04}$	$1.09\times {10}^{-03}$	$3.92\times {10}^{-02}$	$8.61\times {10}^{-01}$	-	-	-
3D-CNN + SSM	$<2.2\times {10}^{-16}$	$6.12\times {10}^{-13}$	$3.82\times {10}^{-11}$	$8.57\times {10}^{-09}$	$6.04\times {10}^{-07}$	$1.19\times {10}^{-07}$	$3.94\times {10}^{-06}$	$8.01\times {10}^{-5}$	$1.04\times {10}^{-05}$	$4.85\times {10}^{-05}$	$2.96\times {10}^{-04}$	$5.01\times {10}^{-03}$	$8.01\times {10}^{-03}$	-	-
SS-DTC	$<2.2\times {10}^{-16}$	$8.12\times {10}^{-12}$	$9.04\times {10}^{-10}$	$1.07\times {10}^{-10}$	$6.28\times {10}^{-09}$	$1.91\times {10}^{-07}$	$7.23\times {10}^{-05}$	$8.98\times {10}^{-04}$	$7.81\times {10}^{-04}$	$2.23\times {10}^{-04}$	$1.76\times {10}^{-03}$	$2.01\times {10}^{-02}$	$1.41\times {10}^{-02}$	$1.09\times {10}^{-02}$	-
MV-HyLP	$<2.2\times {10}^{-16}$	$3.61\times {10}^{-15}$	$6.43\times {10}^{-12}$	$3.19\times {10}^{-10}$	$8.81\times {10}^{-09}$	$8.21\times {10}^{-06}$	$8.87\times {10}^{-05}$	$6.29\times {10}^{-04}$	$1.64\times {10}^{-03}$	$8.71\times {10}^{-02}$	$7.15\times {10}^{-02}$	$5.87\times {10}^{-02}$	$3.41\times {10}^{-02}$	$2.03\times {10}^{-02}$	$1.71\times {10}^{-02}$

10. Conclusions & future work

The proposed MV-KCS successfully incorporates multi-view and semi-supervised learning in the multi-atlas patch-based segmentation framework. Leveraging multiple views on the data to learn a combined sparse graph, it is able to achieve state-of-the-art segmentation performance by integrating information along two distinct axes, spectral – spatial and global – local. The proposed labelling mechanisms within MV-KCS, namely, MV-RegLP and MV-HyLP, were proven to be statistically superior to other existing patch-based and deep-learning based approaches, with MV-HyLP yielding the best results ($\mathit{DSC}=92.56\mathrm{\%}(\mathit{FC})$, $\mathit{DSC}=89.91\mathrm{\%}\text{break}(\mathit{TC})$). Furthermore, a comparison with other state-of-the-art methods in the knee cartilage segmentation literature, firmly places MV-KCS among the top performing methodologies. The implementation for both MV-RegLP and MV-HyLP can be found in .²

While demonstrating satisfactory results, the proposed MV-RegLP and MV-HyLP are not without their limitations. In order for these methods to perform effectively, an atlas library has to be constructed for each target image to be segmented, significantly increasing the computational load and execution time, in case such libraries are not already available. Finally, determination of voxel similarity, which facilitates the construction of sparse graphs, in the original feature space may not be adequate in capturing the full information content of these voxels, possibly leading to sub-optimal representations. A potential remedy for this issue may lie in the use of deep models, that will be better suited to capture the varying degrees of similarity between the data.

Our future work will primarily proceed along three research directions: 1) Evaluation of alternative multi-view semi-supervised learning schemes in knee cartilage segmentation 2) Application of the suggested methods in different medical imaging domains (e.g. brain segmentation). 3) The employment of deeper models within the graph-learning framework (i.e. Graph Neural Networks).

Disclosure statement

No potential conflict of interest was reported by the author(s)

Data availability statement

Due to the nature of the research, due to [ethical/legal/commercial] supporting data is not available.

Additional information

Notes on contributors

Christos G. Chadoulos

Christos G. Chadoulos received the B.S. and M.S. degrees in electrical and computer engineering from the Aristotle University of Thessaloniki, Thessaloniki, Greece, in 2017. He is currently pursuing the Ph.D. degree at the Depart- ment of Electrical and Computer Engineering, Aristotle University of Thessaloniki. His research interests include computer vision, image analysis, and machine learning.

Dimitrios E. Tsaopoulos

Dimitrios E. Tsaopoulos graduated with a BSc from the University of Thessaly, Greece, in the field of Sports Science in 2003. He received his PhD from the Manchester Metropolitan University in 2007, focusing on in vivo human knee joint mechanics. He is currently a Research Director (Grade A) on the field of Biomechanics in CERTH/IBO. Dimitris aspires to the advancement of research and technologies in areas of human and animal-movement biology, aiming at diagnosis, assessment and treatment of locomotor dysfunction. He has participated in more than 20 National (UK and GR) and European research projects.

Serafeim P. Moustakidis

Serafeim P. Moustakidis has wide experience in computational intelligence, machine learning and data processing with more than 12 years of research experience in various fields. His research focuses in developing novel algorithms for solving important existing and/or emerging problems. His main scientific interests cover various application fields such as Deep Learning, Big Data, Biomechanics, Bio-economy, Health, Remote Sensing, Energy Optimization, Non- Destructive Testing (NDT) and machine learning-empowered imaging. He has been involved in the technical implementation, scientific or overall management of 25 R&D projects of a total budget of 35 million Euros. He has worked for several research organisations across Europe.

John B. Theocharis

John B. Theocharis (M’90) received the degree in elec- trical engineering and the Ph.D. degree from the Aristotle University of Thessaloniki, Thessaloniki, Greece, in 1980 and 1985, respectively. He is currently a Professor in the Department of Electrical and Computer Engineering, Aris- totle University of Thessaloniki. His research activities in- clude fuzzy systems, neural networks, evolutionary algo- rithms, pattern recognition and image analysis. He has published numerous papers in several application areas such as neuro-fuzzy modeling, power demand and wind speed prediction, land cover classification and segmenta- tion from remotely sensed images. Recently his research is focused on addressing challenges in medical imaging using machine learning and deep learning techniques.

Notes

1. Here, we make the concession of using ${\mathbf{x}}_i$ interchangeably to refer both to a voxel location, as well as its description.

2. https://gitlab.com/christos_chadoulos/semi-supervised-multi-view-multi-atlas-mri-segmentationMV-RegLP & MV-HyLP implementations