The publicly available database of the Osteoarthritis Initiative (OAI) ¹ was established to provide researchers with resources to promote the prevention and treatment of knee osteoarthritis. We use 2,399 sagittal Double Echo Steady-State (DESS) 3D MRI scans from the OAI database acquired using Siemens Trio 3.0 Tesla scanners (Peterfy et al., 2006). Additionally, 2,396 sagittal Intermediate-Weighted Turbo Spin-Echo (IW TSE) MRI scans are investigated for the same patients. The demographics of our study are shown in Table 1.

The OAI database includes multiple reading studies of respective osteoarthritis characteristics, which can be assessed in medical image data. As a gold standard, we utilize labels from MOAKS (Hunter et al., 2011) image reading studies performed by clinical experts. In the MOAKS scoring system, the menisci are divided into three anatomical sub-regions: anterior horn, body, posterior horn. We consider a sub-region as not containing a tear if the MOAKS score is “normal” or indicates a signal abnormality (which is not extending through the meniscal surface and, hence, is no tear). We considered any other type of abnormality (radial, horizontal, vertical, etc.) as a meniscal tear (c.f. Supplementary Table S1). Examples of the MRI sequences, signal abnormalities, and meniscal tears are shown in Figure 1.

In a first step of our pre-processing, the intensities of all MR images are scaled to a range of [0, 1] using min-max normalization. Following that, a standardization is applied to each MR image I i according to: I ̃ i = I i − μ σ , (1) where μ is the mean intensity and σ is the standard deviation of the training population of normalized scans. Leveraging meniscal segmentations generated by the method of Tack et al. (2018) RoIs spanning the MM and LM are created for DESS MRI data (see Figure 2). RoIs are computed by querying the minimum and maximum position of the menisci along each dimension of the binary segmentation masks: x _min, x _max, y _min, y _max, z _min, z _max. The bounding boxes are uniquely defined as the 3D center coordinate B B c e n t e r = ( x m a x − x m i n ) / 2 + x m i n ( y m a x − y m i n ) / 2 + y m i n ( z m a x − z m i n ) / 2 + z m i n , (2) and with the respective height (x _max − x _min), width (y _max − y _min), and depth (z _max − z _min). These values are represented as relative image coordinates. Hence, a bounding box is defined by 6 floating values: [ B B c e n t e r x , B B c e n t e r y , B B c e n t e r z , h e i g h t , w i d t h , d e p t h ] .

For the IW TSE data 600 segmentations are generated in a semi-automated fashion using Amira ZIB Edition ² (Reddy, 2017). These masks are defined as voxel-wise annotations of the tissue belonging to the respective meniscus. The method of Tack et al. (2018) was originally developed and evaluated on DESS MRI data. Since the DESS and IW TSE MRI sequences differ significantly in the image resolution (number of slices), that could pose an issue, we have decided to train the self-adapting nnU-net framework (Isensee et al., 2021) on these 600 training datasets. The nnU-net offers 2D and 3D architectures with 3D architectures usually yielding better results (Isensee et al., 2021). For this reason, we have used a 3D variant of the nnU-net that employs 3D convolutions in an encoder-decoder framework with skip-connections. For the IW TSE data, the nnU-net has been automatically configured to have an input size of 24 × 256 × 256 pixels and seven layers of 3D convolutions (Isensee et al., 2021). We train the nnU-net with data augmentation such as random rotations and random cropping using a dice similarity coefficient loss (Isensee et al., 2021) until convergence is reached. Hereby, the dice similarity coefficient is computed between the output of the nnU-net and the respective hand-labelled target segmentation masks. Afterwards, the nnU-net is employed to segment all 2,396 IW TSE MRI scans to yield the respective meniscal RoIs. In order to achieve this, multiple patches of the MRI with a size of 24 × 256 × 256 pixels are being processed by the nnU-net. These patches overlap by half of the patch size in each dimension. Afterwards, the nnU-net framework merges all patches to a final 3D segmentation mask employing a majority voting for every pixel.

Two distinct models, which are based on 3D counterparts of ResNet architectures (He et al., 2016; Yu et al., 2017) are introduced. ResNets have been widely applied to the medical domain and provide good properties due to the employed skip connections. In theory, the residual connections allow the design of very deep ResNets without exhibiting problems of vanishing gradients (Ide and Kurita, 2017). We have chosen 3D counterparts of 2D ResNets since 3D convolutions are able to comprehend three-dimensional context inherently. It has previously been shown in the context of musculoskeletal MRI analysis that 3D convolutions are more powerful than concatenation of 2D slices as well as a provision of multiple 2D slices as input to a CNN that employs 2D convolutions (Ambellan et al., 2019; Tack and Zachow, 2019). We adapt these 3D ResNet architectures to the three different approaches and their associated input volume sizes. Each model consists of a ResNet encoder followed by one or two Multi-Layer Perceptron (MLP) heads. The BB-crop approach has a dilation ResNet-C-26 architecture with an MLP head for the multi-label classification. The Full-scale approach has a ResNet50 encoder with a classifier MLP head, and the BB-loss approach consists of a ResNet50 encoder with two MLP heads. The performance of the classification task is improved in the BB-loss approach by solving additionally a second task, which is to learn a bounding box regression simultaneously. Again, the first MLP head is employed for multi-label classification. The second MLP head is responsible for the bounding box regression task. All ResNets comprise of a series of convolutional layers, each followed by batch normalization (Ioffe and Szegedy, 2015) and a Rectified Linear Unit (ReLU) activation function (Agarap, 2018).

Our approaches that will be presented in the following sections are designed based on (a selection of) encoders and MLP heads:

In our first and most straightforward approach, the complete 3D MRI is provided as input to the CNN. The CNN consists of a ResNet50 encoder followed by an MLP head. The outputs of the MLP after a sigmoid activation represent the probabilities for the six meniscal sub-regions to contain a tear.

The CNN is trained by minimizing the binary cross-entropy loss L B C E for a given batch of N samples. With a target matrix Y ∈ Z 2 N × C and an output matrix Y ̂ ∈ R N × C for all C meniscal sub-region labels the definition of L B C E is: L B C E = 1 N ∑ c = 1 C ∑ i = 1 N w c [ y i , c l o g ( σ ( y ̂ i , c ) ) + ( 1 − y i , c ) l o g ( 1 − σ ( y ̂ i , c ) ) ] , (3) where w _c is an inverse weighting of label frequencies and σ(⋅) is a sigmoid activation function. The Full-scale approach is visualized under A) in Figure 2.

Cropping 3D MRI data to the meniscal RoI is expected to provide two desirable properties. First, it provides smaller volumes reducing the required GPU memory as well as the run time. Second, the Full-scale 3D MR images can be considered noisy as they provide additional and unnecessary information about surrounding anatomical structures. By cropping the data to the RoI of the menisci, this unnecessary information is suppressed. Leveraging the RoI generated as described in section 2.2 the 3D MR images are cropped with a 5% margin around the menisci. Each cropped image is then resampled with trilinear interpolation to the closest multiples of 16, given the biggest bounding box in the training set. Figure 2 visualizes the cropping and resampling process. Consequently, the cropped and resampled images have a size of (64, 64, 176) for the DESS data and (16, 64, 176) for the IW TSE data. BB-crop utilizes a Dilation Resnet-C-26 encoder followed by an MLP classifier head. The CNN is trained by minimizing the L B C E as given in Eq. 3. The framework is visualized under B) in Figure 2.

The BB-crop approach requires segmentation of both menisci (or at least the determination of a meniscal region) in training and testing. Since generating segmentations is time-consuming (the method of Tack et al. (2018) requires approximately 5 min of run time), it is beneficial to avoid this step. Moreover, this approach heavily relies on high-quality bounding boxes in training and inference, which are difficult to obtain and strongly influence the performance quality. Thus, the motivation for our final BB-loss approach is to detect meniscal tears in 3D MRI data without extensive pre-processing requirements such as segmenting the menisci or computing bounding boxes for meniscal regions. Instead, the location of the menisci is added as an additional loss term for the training. The encoder is kept identical to the Full-scale approach, namely a ResNet-50 encoder. Furthermore, an identical MLP head is utilized for the meniscal tear detection. Additionally, we show that the meniscal position information helps the CNN to focus on these regions in the image yielding better results. A second MLP head is employed in the BB-loss approach to regress the coordinates of the meniscal RoI. By incorporating this knowledge as a loss in the training process, the locations of the menisci must not be explicitly provided at test time. The total loss in the BB-loss setting is computed considering the multi-label classification and the bounding box regression task. For detection of meniscal tears L B C E is employed (Eq. 3). In the bounding box regression the outputs of the MLP head are 6 coordinates d for the MM and LM, respectively. Utilizing a sigmoid activation function, these values are given as relative positions within the image in a range of [0, 1] of the respective dimension. For a detailed description of the bounding box generation procedure, we refer the reader to section 2.2. The first component of the bounding box loss is an L1-term L L 1 defined as L L 1 = ‖ B − B ̂ ‖ , (4) with a predicted bounding box B ̂ and a target bounding box B that is derived from the automated segmentation masks. These N × 2d matrices contain N rows with medial and lateral bounding box values. Where b _n,i and b ̂ n , i describe the nth element of the batch and the ith value of the concatenated bounding boxes. With this formulation the loss is given as L L 1 = 1 N ∑ n = 1 N ∑ i = 1 2 d | b n , i − b ̂ n , i | . (5)

The second component of the bounding box loss is a modified Intersection over Union (IoU) term, more specifically the Generalized-IoU (GIoU) L G I o U (Rezatofighi et al., 2019) defined as: L G I o U = 1 − I o U + | C \ \ ( B ∪ B ̂ ) | | C | , (6) where C is a convex hull enclosing the predicted and the target box. The operator |⋅| computes the box volume. The convex hull is the smallest possible region that encloses both the output and the target bounding boxes. It can be defined as a bounding box, fully characterised by the 6 coordinates elaborated above. It is computed by taking the minimum and maximum extent of both the target bounding box and the predicted bounding box coordinates along the x-, y- and z-axis. The numerator of the third term of the L G I o U is the convex hull volume subtracted by the volume of B and B ̂ , and the denominator is the volume of the convex hull. Hence, the third term of the L G I o U can be considered as the relative volume of the convex hull not covered by the union of predicted and target bounding box. The IoU is defined as | B ∩ B ̂ | | B ∪ B ̂ | , that is, the ratio of the intersecting voxels of B and B ̂ to their union. The L G I o U is computed for each meniscal RoI and averaged for the given batch. The overall loss L for the BB-loss approach is given as L = L B C E + L L 1 + L G I o U . (7)

The BB-loss approach is visualized under C) in Figure 2.

The given MRI data of the OAI are randomly split into 50% training data, 15% validation data and 35% testing data. Hence, our two experiments have 1200/359/840 and 1197/359/840 training/validation/testing scans for the DESS data and the IW TSE data, respectively. We implemented the CNNs of all approaches in PyTorch 1.9. Convolutional weights are initialized using a normal distribution as in He et al. (2015) tailored towards our deep neural networks with asymmetric ReLU activation functions. While, batch normalization weights and biases are initialized constant with 1 and 0. We train our CNNs on an Nvidia A100 GPU with 40 GB memory. Training our three ResNets, separate learning rates and dropout probabilities for the ResNet-encoders and the MLP-heads are introduced. Suitable learning rate, dropout and batch size hyper-parameters are found using the validation data of the DESS scans. The learning rate values for all parts (ResNet encoder, classifier head and bounding box head) are evaluated in an interval of [ 1 e − 5 , 0.01]. Dropout percentages are varied in an interval of [0.1, 0.9]. Further, the training batch size limited by the input size is varied from 2 to 64 for the BB-crop approach. Due to a larger input volume in approach Full-scale and BB-loss  batch size was kept constant at a value of 4. For a complete summary of our hyper-parameter values, please refer to Supplementary Table S2. Training is performed using the ADAM optimizer (Kingma and Ba, 2014) with β ₁ = 0.9, β ₂ = 0.999 and ϵ = 1 e − 08 with a learning rate decay of 0.5 every 50 epochs. Training on the IW TSE sequence is not performed from scratch, instead, both ResNet encoder and MLP weights are fine-tuned. In both DRN-C-26 and ResNet50 cases, we use the CNNs that achieve the lowest validation loss on the DESS sequence.

On-the-fly data augmentation is performed during training. Specifically, this means, random cropping around the RoI, horizontal flips, rotations, Gaussian noise, and intensity scaling are applied with 50% probability. For the Full-scale approach, we perform random cropping of up to 10% along coronal, 20% sagittal and 20% axial direction. In the BB-crop approach, random crops are performed by uniformly cropping within a 20% margin around the menisci. The BB-loss approach uniformly samples possible crops around the menisci. All cropped images are resampled with trilinear interpolation to attain consistent sizes per approach and dataset. Input images for the Full-scale and BB-loss approach are sampled for the DESS sequence data to (160, 384, 384) and for IW TSE images to (44, 448, 448). The BB-crop approach resamples to (64, 64, 176) and (16, 64, 176), respectively. The added Gaussian noise is pixel-wise sampled as ϵ ∈ N ( 0.1 , 0.5 ) . The random rotation is uniformly sampled from U ( − 5 ° , + 5 ° ) and image intensity is scaled by a uniformly sampled multiplication factor b ∈ U ( 0.9 , 1.1 ) .

For all experiments, we plot the true positive rate (TPR = sensitivity) against the false positive rate (FPR = 1–specificity) at various decision thresholds to create ROC curves (Brown and Davis, 2006). Additionally, we compute the ROC AUC to assess the quality of our classifiers. The quality of our predicted bounding boxes is assessed by computing the IoU with the target bounding boxes. We consider IoU values over 0.5 as successful localization of the menisci since this is a common value in object detection tasks (Girshick et al., 2014).

Gradient saliency maps (Simonyan et al., 2013) (otherwise called pixel attribution maps or sensitivity maps) highlight pixel regions in the input image that mostly influenced a neural network’s decision. To attain such pixel attributions, one computes the derivative of the final linear layer in a neural network with respect to the input via back-propagation. More formally, a gradient saliency map S _c for a sub-region c for which our neural network f yields a detection of meniscal tears is calculated as: S c ( I ̃ i ) = ∂ f c ( I ̃ i ) ∂ I ̃ i . (8)

For our two most promising approaches BB-crop and BB-loss, these maps are computed by applying a slight enhancement to the original mechanism - the SmoothGrad method (Smilkov et al., 2017). Similar to the SmoothGrad approach of Smilkov et al. (2017) we augmented the input image slightly, introducing noise, such that through averaging, the saliency maps of different noise levels are smoothed out. We apply Gaussian distributed noise ϵ ∈ N ( 0.1 , 0.5 ) , random horizontal flips, uniformly sampled rotations r ∈ U ( − 5 ° , + 5 ° ) and uniformly sampled pixel intensity shift with a multiplication factor b ∈ U ( 0.9 , 1.1 ) . Each image is augmented 20 times with a probability of 50% per augmentation, and the resulting maps are averaged.

Employing the Full-scale approach, the AUC values are 0.74, 0.84, 0.85 for the anterior horn, body, and posterior horn of the MM. For the LM, the AUC values are 0.94, 0.92, 0.91. The BB-crop approach usually yields higher AUC values, being 0.87, 0.89, 0.89 and 0.95, 0.93, 0.91. The BB-loss gives the highest AUC values, being 0.94, 0.93, 0.93 and 0.96, 0.94, 0.91. The ROC curves employing all three approaches are shown in Figure 3. In addition, all ROC AUC results are summarized in Table 2.

Employing the Full-scale approach, the AUC values are 0.82, 0.87, 0.82 for the anterior horn, body, and posterior horn of the MM. For the LM, the AUC values are 0.88, 0.85, 0.85. The BB-crop approach usually yields higher AUC values, being 0.84, 0.89, 0.86, and 0.92, 0.90, 0.90. The BB-loss gives similar AUC values, being 0.84, 0.88, 0.86, and 0.95, 0.91, 0.90. The ROC curves of all approaches are shown in Figure 4. Further, all AUC values are summarized in Table 3.

To investigate the bounding box regression quality of the proposed method we evaluate the distribution of the IoU values for the predicted bounding boxes (Figure 5). For the DESS dataset (our primary benchmark), we observed a very high quality of MM and LM bounding box predictions. With the values being close to normally distributed around a mean value of 0.71 (95% confidence interval (CI): 0.71–0.72) and standard deviation of 0.13. With the IoU threshold of 0.5, we conclude that 95% of the resulted bounding boxes are identified correctly. Unfortunately, we observed a clear decrease in the object detection performance in the IW TSE dataset. With a mean value of 0.58 (95% CI: 0.57–0.59) and a standard deviation of 0.14. Applying the same detection threshold as above we testify, that only around 76% of menisci were detected correctly, with the overall quality of the bounding boxes being more widely spread.

Figure 6 shows SmoothGrad saliency maps for the BB-crop and BB-loss approach overlaid to MR images. Examples are shown for randomly selected test cases, displaying different kinds of meniscal tears for DESS and IW TSE data. The RoIs for the BB-loss approach were extracted using predicted bounding boxes and the respective close-ups are shown. Red arrows point at the location of meniscal defects. Most saliency maps obtained this way display a plausible localization of the meniscal tears. The plausibility of these maps was qualitatively evaluated by their correspondence to the target labels of the regions in which the tears could also be confirmed with the help of visual inspection of the image data. SmoothGrad saliency maps are capable of highlighting more than just one affected sub-region, i.e., in the presence of defects in multiple sub-regions of one meniscus, one similarly observes these being correctly highlighted. With the Dilation ResNet-C-26 employed in the BB-crop approach, we observed that this CNN yields smoother and less noisy SmoothGrad saliency maps. However, in many cases, ResNet-50 saliency maps targeted the affected region better, but did not outline this region sharply.

Even though the occurrence of defects varies between meniscal sub-regions (see Supplementary Figure S2), we observe only minimal differences between AUC values of sub-regions in DESS MRI data (c.f. Table 2). However, we analyzed the false positive classifications and found that for all sub-regions, signal abnormalities were more often misclassified than normal menisci were (see Supplementary Figure S1). The misclassification rate of signal abnormalities is highest for the posterior horn of the lateral meniscus, the region with the least AUC for the DESS data. Conversely, the lowest signal abnormality misclassification rate is prevalent in the posterior horn of the medial meniscus, the sub-region with the highest number of signal abnormalities (Supplementary Table S1).

The least common types of tears occurring in the data are radial and vertical tears, amounting to 72 and 69, respectively. Vertical tears were most challenging for our method to detect in DESS data and led to the most false negative results (see Supplementary Figure S2). Radial meniscal tears were the ones yielding the second highest rate of misclassifications.