21. Tong, N.; Gou, S.; Yang, S.; Ruan, D.; Sheng, K. Fully Automatic Multi-organ Segmentation for Head and Neck Cancer Radiotherapy Using Shape Representation Model Constrained Fully Convolutional Neural Networks. *Med. Phys.* **2018**, *45*, 4558–4567. [CrossRef]

22. Ronneberger, O.; Fischer, P.; Brox, T. *U-Net: Convolutional Networks for Biomedical Image Segmentation*; Springer: Berlin/Heidelberg, Germany, 2015; pp. 234–241.

23. Liu, L.; Cheng, J.; Quan, Q.; Wu, F.-X.; Wang, Y.-P.; Wang, J. A Survey on U-Shaped Networks in Medical Image Segmentations. *Neurocomputing* **2020**, *409*, 244–258. [CrossRef]

24. Wang, Y.; Zhao, L.; Wang, M.; Song, Z. Organ at Risk Segmentation in Head and Neck Ct Images Using a Two-Stage Segmentation Framework Based on 3D U-Net. *IEEE Access* **2019**, *7*, 144591–144602. [CrossRef]

25. Rundo, L.; Han, C.; Nagano, Y.; Zhang, J.; Hataya, R.; Militello, C.; Tangherloni, A.; Nobile, M.S.; Ferretti, C.; Besozzi, D.; et al. USE-Net: Incorporating Squeeze-and-Excitation Blocks into U-Net for Prostate Zonal Segmentation of Multi-Institutional MRI Datasets. *Neurocomputing* **2019**, *365*, 31–43. [CrossRef]

26. Schlemper, J.; Oktay, O.; Schaap, M.; Heinrich, M.; Kainz, B.; Glocker, B.; Rueckert, D. Attention Gated Networks: Learning to Leverage Salient Regions in Medical Images. *Med. Image Anal.* **2019**, *53*, 197–207. [CrossRef] [PubMed]

27. Yeung, M.; Sala, E.; Schönlieb, C.-B.; Rundo, L. Focus U-Net: A Novel Dual Attention-Gated CNN for Polyp Segmentation during Colonoscopy. *Comput. Biol. Med.* **2021**, *137*, 104815. [CrossRef] [PubMed]

28. Liang, X.; Nguyen, D.; Jiang, S.B. Generalizability Issues with Deep Learning Models in Medicine and Their Potential Solutions: Illustrated with Cone-Beam Computed Tomography (CBCT) to Computed Tomography (CT) Image Conversion. *Mach. Learn. Sci. Technol.* **2020**, *2*, 015007. [CrossRef]

29. Qiu, B.; van der Wel, H.; Kraeima, J.; Glas, H.H.; Guo, J.; Borra, R.J.H.; Witjes, M.J.H.; van Ooijen, P.M.A. Automatic Segmentation of Mandible from Conventional Methods to Deep Learning—A Review. *J. Pers. Med.* **2021**, *11*, 629. [CrossRef]

30. Hesse, L.S.; Kuling, G.; Veta, M.; Martel, A.L. Intensity Augmentation to Improve Generalizability of Breast Segmentation Across Different MRI Scan Protocols. *IEEE Trans. Biomed. Eng.* **2021**, *68*, 759–770. [CrossRef]

31. Bosniak, M.A. The Current Radiological Approach to Renal Cysts. *Radiology* **1986**, *158*, 1–10. [CrossRef]

32. Bluemke, D.A.; Moy, L.; Bredella, M.A.; Ertl-Wagner, B.B.; Fowler, K.J.; Goh, V.J.; Halpern, E.F.; Hess, C.P.; Schiebler, M.L.; Weiss, C.R. Assessing Radiology Research on Artificial Intelligence: A Brief Guide for Authors, Reviewers, and Readers—From the *Radiology* Editorial Board. *Radiology* **2020**, *294*, 487–489. [CrossRef]

33. Kim, H.; Shim, E.; Park, J.; Kim, Y.-J.; Lee, U.; Kim, Y. Web-Based Fully Automated Cephalometric Analysis by Deep Learning. *Comput. Methods Programs Biomed.* **2020**, *194*, 105513. [CrossRef]

34. Tao, Q.; Yan, W.; Wang, Y.; Paiman, E.H.M.; Shamonin, D.P.; Garg, P.; Plein, S.; Huang, L.; Xia, L.; Sramko, M.; et al. Deep Learning-Based Method for Fully Automatic Quantification of Left Ventricle Function from Cine MR Images: A Multivendor, Multicenter Study. *Radiology* **2019**, *290*, 81–88. [CrossRef] [PubMed]

35. B, S.; R, N. Transfer Learning Based Automatic Human Identification Using Dental Traits- An Aid to Forensic Odontology. *J. Forensic Leg. Med.* **2020**, *76*, 102066. [CrossRef] [PubMed]

36. Ghafoorian, M.; Mehrtash, A.; Kapur, T.; Karssemeijer, N.; Marchiori, E.; Pesteie, M.; Guttmann, C.R.G.; de Leeuw, F.-E.; Tempany, C.M.; van Ginneken, B.; et al. Transfer Learning for Domain Adaptation in MRI: Application in Brain Lesion Segmentation. In *International Conference on Medical Image Computing and Computer-Assisted Intervention*; Springer: Berlin/Heidelberg, Germany, 2017; pp. 516–524. [CrossRef]

37. Lee, K.-S.; Jung, S.-K.; Ryu, J.-J.; Shin, S.-W.; Choi, J. Evaluation of Transfer Learning with Deep Convolutional Neural Networks for Screening Osteoporosis in Dental Panoramic Radiographs. *J. Clin. Med.* **2020**, *9*, 392. [CrossRef] [PubMed]

38. Weiss, K.; Khoshgoftaar, T.M.; Wang, D. A Survey of Transfer Learning. *J. Big Data* **2016**, *3*, 9. [CrossRef]

39. He, K.; Zhang, X.; Ren, S.; Sun, J. Deep Residual Learning for Image Recognition. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Las Vegas, NV, USA, 27–30 June 2016; pp. 770–778.

40. Frei, S.; Cao, Y.; Gu, Q. Algorithm-Dependent Generalization Bounds for Overparameterized Deep Residual Networks. *arXiv* **2019**, arXiv:1910.02934.

41. Huang, K.; Tao, M.; Wang, Y.; Zhao, T. Why Do Deep Residual Networks Generalize Better than Deep Feedforward Networks? — A Neural Tangent Kernel Perspective. 2020, 12. *arXiv* **2020**, arXiv:2002.06262.

42. Ang, K.K.; Zhang, Q.; Rosenthal, D.I.; Nguyen-Tan, P.F.; Sherman, E.J.; Weber, R.S.; Galvin, J.M.; Bonner, J.A.; Harris, J.; El-Naggar, A.K. Randomized Phase III Trial of Concurrent Accelerated Radiation plus Cisplatin with or without Cetuximab for Stage III to IV Head and Neck Carcinoma: RTOG 0522. *J. Clin. Oncol.* **2014**, *32*, 2940. [CrossRef]

43. Nikolov, S.; Blackwell, S.; Zverovitch, A.; Mendes, R.; Livne, M.; De Fauw, J.; Patel, Y.; Meyer, C.; Askham, H.; Romera-Paredes, B. Deep Learning to Achieve Clinically Applicable Segmentation of Head and Neck Anatomy for Radiotherapy. *arXiv* **2018**, arXiv:1809.04430.

44. Clark, K.; Vendt, B.; Smith, K.; Freymann, J.; Kirby, J.; Koppel, P.; Moore, S.; Phillips, S.; Maffitt, D.; Pringle, M. The Cancer Imaging Archive (TCIA): Maintaining and Operating a Public Information Repository. *J. Digit. Imaging* **2013**, *26*, 1045–1057. [CrossRef]

45. Bosch, W.R.; Straube, W.L.; Matthews, J.W.; Purdy, J.A. Data from Head-Neck_cetuximab. *Cancer Imaging Arch.* **2015**, *10*, K9.

46. Han, X. Automatic Liver Lesion Segmentation Using a Deep Convolutional Neural Network Method. *Med. Phys.* **2017**, *44*, 1408–1419. [CrossRef] [PubMed]
47. Bulat, A.; Tzimiropoulos, G. Binarized Convolutional Landmark Localizers for Human Pose Estimation and Face Alignment with Limited Resources. In Proceedings of the IEEE International Conference on Computer Vision, Venice, Italy, 22–29 October 2017; pp. 3706–3714.

48. Szegedy, C.; Liu, W.; Jia, Y.; Sermanet, P.; Reed, S.; Anguelov, D.; Erhan, D.; Vanhoucke, V.; Rabinovich, A. Going Deeper with Convolutions. In Proceedings of the Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Boston, MA, USA, 7–12 June 2015; pp. 1–9.

49. Seo, H.; Huang, C.; Bassenne, M.; Xiao, R.; Xing, L. Modified U-Net (MU-Net) with Incorporation of Object-Dependent High Level Features for Improved Liver and Liver-Tumor Segmentation in CT Images. *IEEE Trans. Med. Imaging* **2019**, *39*, 1316–1325.
[CrossRef] [PubMed]

50. Gupta, S.; Gupta, R.; Ojha, M.; Singh, K.P. A Comparative Analysis of Various Regularization Techniques to Solve Overfitting Problem in Artificial Neural Network. In Proceedings of the Data Science and Analytics; Panda, B., Sharma, S., Roy, N.R., Eds.; Springer: Singapore, 2018; pp. 363–371.

51. He, F.; Liu, T.; Tao, D. Why ResNet Works? Residuals Generalize. *arXiv* **2019**, arXiv:1904.01367.
<figure><img src="image_35.png" /></figure>
Article

# On Unsupervised Methods for Medical Image Segmentation: Investigating Classic Approaches in Breast Cancer DCE-MRI

Carmelo Militello ¹*, Andrea Ranieri ², Leonardo Rundo ³,⁴,⁵, Ildebrando D'Angelo ⁶,⁷, Franco Marinozzi ², Tommaso Vincenzo Bartolotta ⁶,⁸, Fabiano Bini ² and Giorgio Russo ¹

¹ Institute of Molecular Bioimaging and Physiology, Italian National Research Council (IBFM-CNR), 90015 Cefalu, Italy; giorgio.russo@ibfm.cnr.it

² Department of Mechanical and Aerospace Engineering, Sapienza University of Rome, 00184 Roma, Italy; ranieri.1837003@studenti.uniroma1.it (A.R.); franco.marinozzi@uniroma1.it (F.M.); fabiano.bini@uniroma1.it (F.B.)

³ Department of Radiology, University of Cambridge, Cambridge CB2 0QQ, UK; lr495@cam.ac.uk

⁴ Cancer Research UK Cambridge Centre, Cambridge CB2 0RE, UK

⁵ Department of Information and Electrical Engineering and Applied Mathematics (DIEM), University of Salerno, 84084 Fisciano, Italy

⁶ Department of Radiology, Fondazione Istituto "G. Giglio", 90015 Cefalu, Italy; ildebrando.dangelo@hsrgiglio.it (I.D.); tommasovincenzo.bartolotta@unipa.it (T.V.B.)

⁷ Breast Unit, Fondazione Istituto "G. Giglio", 90015 Cefalu, Italy

⁸ Section of Radiology—Department of Biomedicine, Neuroscience and Advanced Diagnostics (BiND), University Hospital "Paolo Giaccone", 90127 Palermo, Italy

* Correspondence: carmelo.militello@ibfm.cnr.it

**Abstract:** Unsupervised segmentation techniques, which do not require labeled data for training and can be more easily integrated into the clinical routine, represent a valid solution especially from a clinical feasibility perspective. Indeed, large-scale annotated datasets are not always available, under-mining their immediate implementation and use in the clinic. Breast cancer is the most common cause of cancer death in women worldwide. In this study, breast lesion delineation in Dynamic Contrast Enhanced MRI (DCE-MRI) series was addressed by means of four popular unsupervised segmentation approaches: Split-and-Merge combined with Region Growing (SMRG), k-means, Fuzzy C-Means (FCM), and spatial FCM (sFCM). They represent well-established pattern recognition techniques that are still widely used in clinical research. Starting from the basic versions of these segmentation approaches, during our analysis, we identified the shortcomings of each of them, proposing improved versions, as well as developing ad hoc pre- and post-processing steps. The obtained experimental results, in terms of area-based—namely, Dice Index (DI), Jaccard Index (JI), Sensitivity, Specificity, False Positive Ratio (FPR), False Negative Ratio (FNR)—and distance-based metrics—Mean Absolute Distance (MAD), Maximum Distance (MaxD), Hausdorff Distance (HD)—encourage the use of unsupervised machine learning techniques in medical image segmentation. In particular, fuzzy clustering approaches (namely, FCM and sFCM) achieved the best performance. In fact, for area-based metrics, they obtained DI = 78.23% ± 6.50 (sFCM), JI = 65.90% ± 8.14 (sFCM), sensitivity = 77.84% ± 8.72 (FCM), specificity = 87.10% ± 8.24 (sFCM), FPR = 0.14 ± 0.12 (sFCM), and FNR = 0.22 ± 0.09 (sFCM). Concerning distance-based metrics, they obtained MAD = 1.37 ± 0.90 (sFCM), MaxD = 4.04 ± 2.87 (sFCM), and HD = 2.21 ± 0.43 (FCM). These experimental findings suggest that further research would be useful for advanced fuzzy logic techniques specifically tailored to medical image segmentation.

**Keywords:** medical image segmentation; breast cancer; pattern recognition; machine learning; clinical feasibility; magnetic resonance imaging; computer-assisted segmentation

## 1. Introduction

The use of advanced imaging technologies has significantly improved the quality of medical care delivered to patients, allowing medical imaging to be an essential part of today's healthcare system [1]. In fact, medical imaging comprises techniques for acquiring

**Citation:** Militello, C.; Ranieri, A.; Rundo, L.; D'Angelo, I.; Marinozzi, F.; Bartolotta, T.V.; Bini, F.; Russo, G. On Unsupervised Methods for Medical Image Segmentation: Investigating Classic Approaches in Breast Cancer DCE-MRI. *Appl. Sci.* **2022**, *12*, 162. https://doi.org/10.3390/app12010162

Academic Editor: Syoji Kobashi

Received: 10 November 2021

Accepted: 19 December 2021

Published: 24 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_36.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
images that convey detailed information about the anatomy and physiology of the imaged organs [2]. Moreover, many imaging-enabled tools were developed, supporting clinicians in several tasks of the care process: assisted segmentation [3,4], diagnosis support [5], treatment response assessment [6], radiomic analyses [7,8]. For these reasons, computer-assisted image analysis is considered an essential instrument in the clinical workflow [9,10].

Despite the technological progress characterizing the modern era, some medical tasks, such as image annotation, are still performed manually, often via time-consuming and operator-dependent procedures. As an example, the Gross Tumor Volume (GTV) segmentation for radiotherapy treatments is usually delineated by means of a fully manual procedure [11,12]. Considering that dozens of slices have to be contoured, this manual process is extremely time-consuming. Moreover, operator dependence is critical in terms of result reproducibility. Indeed, these manual procedures are strongly dependent on clinician's knowledge and experience: this means that a remarkable intra- and inter-operator variability can seriously affect the segmentation and quantification results.

In this scenario, computer-assisted approaches (automatic or semi-automatic) allow us to mitigate some of the typical drawbacks of manual procedures. Semi-automatic segmentation techniques involve a minimal level of user interactions and exploit automated algorithms to produce accurate and repeatable results. User interaction, for example, may involve the selection of an approximate initial ROI, which is subsequently used to segment the image [13]. As a consequence, semi-automatic approaches provide more reproducible measurements—compared to fully-manual ones—with a significant reduction of the segmentation time. In particular, fully-automatic segmentation approaches do not require any user interaction. Most existing fully-automatic approaches exploit machine learning or deep learning techniques—such as Support Vector Machines (SVMs) or deep Convolutional Neural Networks (CNNs)—to successfully handle the variability characterizing biomedical data [14,15]. In this scenario, supervised learning techniques are, in general, more complex since they require high computation times and a large amount of labeled data for training. Furthermore, it is important to point out that—from a clinical feasibility perspective—machine learning approaches, which do not require training, are advantageous: the amount of labeled data needed to adequately train and evaluate the approaches based on deep learning is not always available [16,17].

The aim of this work is to show the potential of unsupervised pattern recognition techniques, which do not require training and can be more easily integrated into care routine, especially from a clinical feasibility perspective. For this reason, supervised approaches—such as CNNs—were not treated here. As a relevant case study, we consider the segmentation of contrast-enhancing masses on DCE-MRI. It is worth noting that, in [18], classical unsupervised techniques, both automatic and semi-automatic, allowed us to obtain results comparable or superior to the deep learning approaches.

From experimental evidence, it is possible to observe that fuzzy clustering techniques significantly outperformed direct region detection approaches (i.e., split-and-merge and region growing) and crisp k-means. FCM and sFCM obtained comparable results, although the integration of spatial information into the sFCM allowed for the best performance. Therefore, the explicit management of segmentation uncertainty via multiple degrees of class memberships, along with spatial information, represented the best computational framework for the problem at hand. The main contributions of this study are:

* unsupervised segmentation methods, based on classic pattern recognition techniques, can still provide an effective solution despite the prevalence of supervised approaches, such as deep CNNs;

* the use of traditional unsupervised approaches—requiring no training—is an advantage in terms of immediate clinical feasibility;

* fuzzy clustering techniques significantly outperformed direct region detection approaches and crisp k-means;

* the integration of spatial information into the FCM algorithm achieved the best performance;
* with the goal of providing a guide for beginners, as well as possibly enabling new future extensions from other researchers, this study provides all the technical information needed to understand both the functioning of each of the studied algorithms and the implemented workflow.

The remainder of this work is organized as follows. Section 2 introduces the theoretical background about unsupervised segmentation approaches, focusing on the algorithms used in this study. A detailed description of the performed analysis and of the implemented processing steps is proposed in Section 3, where the exploited DCE-MRI dataset is also described. Section 4 formulates the metrics used to evaluate the performance of the analyzed approaches. Section 5 illustrates the experimental results, including a discussion about the comparison of the proposed techniques. Finally, discussion and conclusions are provided in Section 6.

## 2. Theoretical Background

A number of algorithms and techniques for image segmentation have been developed and implemented over the years, and a large amount of literature papers about this non-trivial task were proposed. The aim of this section is to present a brief overview about theoretical notions of literature approaches from which this work drawn inspiration. For this reason, we decided for a comprehensive description, to provide the reader with all the technical information needed to understand the functioning of each of the investigated algorithms, as well as the implemented workflow.

The segmentation involves the image partitioning into homogeneous and meaningful sub-regions. By a formal point of view, the segmentation of an image $\mathcal{I}$ involves the identification of a finite set of regions $\mathcal{R}_1, \mathcal{R}_2, \dots, \mathcal{R}_N$ as in Equation (1):

$$ \bigcup_{i=1}^{N} \mathcal{R}_i = \mathcal{I}' \qquad (1) $$

with the following constraints:

$$ \mathcal{R}_i \cap \mathcal{R}_j = \emptyset, \quad \text{for } i \neq j, \qquad (2) $$

$$ P(\mathcal{R}_i) = \text{TRUE}, \quad \text{for } i = 1, 2, \dots, N, \qquad (3) $$

$$ P(\mathcal{R}_i \cup \mathcal{R}_j) = \text{FALSE}, \quad \text{for } i \neq j. \qquad (4) $$

With more details, P is an appropriate logical predicate leading the segmentation process. Equation (1) states that the union of all the sub-regions resulting from segmentation process. Equation (2) points out that the intersection of two different sub-regions is the empty set: this means that the segmented sub-regions do not overlap each other. According to Equation (3), the result of the logical predicate P on all the pixels belonging to the same sub-region is always TRUE; in other words, all the pixels belonging to the same region share the same characteristics according to the predicate P. As a consequence, Equation (4) states that the result of P on the union of two distinct sub-regions is FALSE.

It is important to clarify that each type of medical image has a specific set of features reflecting its own properties: in fact, each image is the result of a complex interaction between the human body and the scanner (i.e., X-rays for CT, magnetic fields for MRI, radioactive decay for nuclear medicine exams) [19]. As a consequence, not all the segmentation techniques obtain the same results on all image types, but some algorithms yield better results when applied on a specific kind of image. Furthermore, it is necessary to point out that, typically, an approach that works very well with one type of image does not mean that it continues to perform well even on different images. As a matter of fact, bioimages exhibit a very high variability, as they depend on various factors, both intrinsic (e.g., patients) and extrinsic (e.g., imaging modalities, acquisition parameters). Accordingly,
ad hoc modifications might be required at the level of both the segmentation approach and
the pre-/post-processing phases.

In what follows, we outline the techniques investigated and compared for contrast-
enhancing mass segmentation on DCE-MRI.

## 2.1. Split-and-Merge Combined with Region Growing

The simplest segmentation approaches use a global threshold applied on pixel intensity to partition the original image: pixels with an intensity greater than threshold *T* are assigned to one region, while those below the threshold *T* are assigned to another one [20–23]. In this way, a binary image is created providing the segmentation of the original image with respect to the chosen threshold value. Split-and-Merge and Region Growing (SMRG) is basically a threshold-based approach that combines Split-and-Merge (SM)—composed of a first split (top-down) phase of the image followed by a merge (bottom-up) phase—with RG for the refinement of the identified regions [23].

The SM algorithm represents a valid alternative to thresholding, because it can find homogeneous regions in terms of uniformity criteria [24–27]. Unlike SM, the RG algorithm, starting from one or more seed-points, identifies an ROI through a growing procedure guided by appropriate similarity properties that describe ROI intensity features. Generally SR and RG are used individually for image segmentation, but the use of both together allows us to exploit the overall potential [28,29].

With more details, the idea behind the SM algorithm involves successive splits of the whole image into disjoint regions lead by a logical predicate P. The algorithm starts with an arbitrary partition $\mathcal{R}$ of the original image (i.e., the whole image) and yields an output composed of uniform sub-regions $\mathcal{R}_i$, for $i = 1, 2, ..., n$, according to the logical criterion expressed by P. At the generic step $t$, if $P(\mathcal{R}_i) = FALSE$, each region $\mathcal{R}_i$ is split into four sub-regions (also called 'quad-regions'): this process iteratively continues until a quad-region such that $P(\mathcal{R}_i) = TRUE$ or with an area smaller than a certain threshold is found. The logical predicate chosen for this work allows us to find the quad-regions with a mean intensity that is greater than the threshold value yielded by the Otsu's method [30]. If only splitting is used, the final partition contains adjacent regions with identical properties: this drawback can be overcome by merging only adjacent regions where the combined pixels satisfy the predicate P.

After the initial rough ROI identification obtained by means of SM, the RG algorithm expands this initial region to properly identify the lesion. According to the classic RG algorithm, each region begins its own growth from a single pixel (seed-point). Instead, when the SM is exploited, the seed-point can be obtained by using the ROI yielded by the SM algorithm (more correctly, a seed-region) [23]: the seed-region is iteratively grown by evaluating, for each pixel on the boundary, its 8-neighborhood as candidate for the growth. A stopping rule is necessary for interrupting the growing procedure if no more pixels match the membership criterion, which often refers to the proximity of the pixel intensities.

## 2.2. K-Means

Classification refers to data labeling into disjoint sets according to a common set of features. Among these, clustering algorithms can be used to determine the natural structures in the data. Clustering algorithms can use more sophisticated properties of the image: in digital imaging, this means that spatial and/or spectral features concerning pixels can be exploited [31].

K-means is an unsupervised clustering technique that aims at partitioning an input set of *N* observations into *k* clusters [32]. Let $\mathcal{X} = \{x_1, x_2, ..., x_N\}$ be a set of vector observations such that $x_i \in \mathbb{R}^n$ for $i = 1, 2, ..., N$. In image segmentation, each component of a vector *x* represents a numerical pixel attribute: if segmentation is based on gray-scale intensity alone, the *n*-dimensional observation $x_i \in \mathbb{R}^n$ degenerates into the scalar value $x_i \in \mathbb{R}$ representing the intensity of the *i*-th pixel. The final purpose of k-means is to partition the *N* observations into *k* < *N* disjoint cluster sets $C = \{C_1, C_2, ..., C_k\}$ such that the sum of the distances from each point in a set to the mean of that set is minimum. From
a mathematical point of view, k-means turns classification into an optimization problem with the following cost function in Equation (5):

$$ \mathcal{J}(\mathcal{X}, \mathcal{V}) = \underset{c}{\operatorname{argmin}} \left( \sum_{i=1}^{k} \sum_{x \in C_i} |x - v_i|^2 \right), \qquad (5) $$

where $v_i$ is the centroid of the samples in the set $C_i$ for $i = 1, 2, \dots, k$. The function $\mathcal{J}(\mathcal{X}, \mathcal{V})$ has no analytical solution: as a result, k-means proceeds by iteratively finding the minimum of its cost function. In particular, at each iteration $t$ the centroid values of the prototype set $\hat{Y}^{(t)}$ are updated according to Equation (6).

$$ \hat{\delta}_i^{(t)} = \frac{1}{|C_i^{(t)}|} \sum_{x \in C_i} x', \qquad (6) $$

where $|C_i^{(t)}|$ is the number of objects belonging to the $i$-th at the step $t$.

## 2.3. Fuzzy C-Means

Biomedical images are characterized by an intrinsic uncertainty, thus causing not well-defined regions (e.g., blurry boundaries or poor anatomical details). This aspect makes thresholding approaches and crisp approaches (such as k-means) not always suitable. Therefore, the natural fuzziness characterizing the Fuzzy C-Means (FCM) clustering approach allows to reach better segmentation results than the hard partitioning offered by k-means [33–36].

The FCM algorithm is an unsupervised clustering technique that searches for the optimal partition of an input data set. The idea leading the FCM classification process is that of minimizing the intra-cluster variance as well as maximizing the inter-cluster variance, in terms of a distance metrics between the feature vectors. Formally, the FCM technique searches for the optimal partition of an input data set $\mathcal{X} = \{x_1, x_2, \dots, x_N\}$ of $N$ objects into $C$ clusters. With respect to the k-means algorithm in which each point is assigned to one cluster only, FCM allows each object to belong to multiple clusters with different degrees of membership. This 'soft' classification allows us to define a fuzzy partition, $\mathcal{P}$ defined as a fuzzy set family $\mathcal{P} = \{Y_1, Y_2, \dots, Y_C\}$ such that each point can have a partial membership to multiple clusters. In mathematical terms, the matrix $U = [u_{ik}] \in \mathbb{R}^{C \times N}$ denotes a fuzzy C-partition of the data set $X$ by means of $C$ membership functions $u_i : X \to [0, 1]$, whose values $u_{ik} := u_i(x_k) \in [0, 1]$ represent membership grades of each element $x_k$ to the $i$-th fuzzy cluster $Y_i$, and have to hold the constraints in Equation (7):

$$ \left\{ \begin{array}{l} 0 \le u_{ik} \le 1 \\ \sum_{i=1}^{C} u_{ik} = 1, \forall k \in 1, 2, \dots, N \\ 0 < \sum_{i=1}^{N} u_{ik} < N, \forall i \in 1, 2, \dots, C \end{array} \right. . \qquad (7) $$

Although hard clustering works well on compact and well-separated groups of data, in many real-world situations clusters overlap each other: as a result, assigning them with gradual memberships by exploiting a soft computing approach may be more appropriate. Computationally, the FCM algorithm assigns to the sample $x_k$ the membership function values using the relative distance (i.e., intensity value similarity) of $x_k$ from the $C$ prototype points $\mathcal{V} = \{v_1, v_2, \dots, v_C\}$ identifying the centroids of the $C$ clusters. Such as k-means algorithm, FCM may be rewritten as an optimization problem with respect to the objective function in Equation (8).

$$ \mathcal{J}_m(U, \mathcal{V}; \mathcal{X}) = \sum_{i=1}^{C} \sum_{k=1}^{N} (u_{ik})^m |x_k - v_i|^2, \qquad (8) $$
where:

* m is the fuzzification constant (i.e., a weighting exponent such that $1 \le m < \infty$) controlling the fuzziness of the classification process. If $m = 1$, the FCM algorithm degenerates to a k-means clustering: in general, the higher the m value the greater will be the fuzziness degree (the most common value is $m = 2$);

* U is the fuzzy C-partition of the set $\mathcal{X}$;

* $|x_k - v_i|^2$ is the Euclidean distance between the elements $x_k$ and the centroid $v_i \in \mathcal{V}$.

Considering that the optimization problem described by FCM does not have a closed-form solution, the minimum of the cost function $\mathcal{J}_m(U, \mathcal{V}; \mathcal{X})$ has to be found iteratively. In particular, at each iteration t, the centroid values of the prototype set $\hat{\mathcal{V}}^{(t)}$ and the elements of the matrix $\hat{U}^{(t)}$ are updated according to Equations (9) and (10), respectively.

$$ \hat{v}_i^{(t)} = \frac{\sum_{j=1}^{N} (\hat{u}_{ij}^{(t)})^m x_j}{\sum_{j=1}^{N} (\hat{u}_{ij}^{(t)})^m}, \quad (9) $$

$$ \hat{u}_{ik}^{(t)} = \left( \sum_{j=1}^{C} \left( \frac{|x_k - \hat{v}_i^{(t)}|}{|x_k - \hat{v}_j^{(t)}|} \right)^{\frac{2}{m-1}} \right)^{-1}, \quad (10) $$

with $m > 1$ and $x_k \ne \hat{v}_j^{(t)}$, $\forall j, k$.

At each iteration, each object $x_k$ is compared with the elements of the centroid vector and is assigned to the nearest cluster. The process stops when the convergence condition (i.e., the matrix norm distance between $\hat{U}^{(t+1)}$ and $\hat{U}^{(t)}$ is less than a fixed value (i.e., minimum improvement in the objective function $\mathcal{J}$ between two consecutive iterations $\epsilon$) or the maximum number of iterations $T_{max}$ is reached. After the convergence, a defuzzification is applied to assign each pixel to the cluster with the highest membership degree, thus achieving a binary classification.

## 2.4. Spatial Fuzzy C-Means

The traditional FCM clustering does not take into account spatial relationship among neighboring pixels, making it sensitive to noise and other imaging artifacts [37,38]. Breast lesions generally tend to grow in an isotropic way, preserving a pseudo-spherical appearance [39]. Relying on those features, it is expected that neighbouring pixels in a digital image are highly correlated and that the probability that they belong to the same cluster is great. Therefore, the use of the sFCM, taking advantage of spatial relationship of neighbouring pixels, can help image segmentation [40,41]. The spatial function used by sFCM is defined in Equation (11):

$$ h_{ij} = \sum_{k \in \mathcal{N}(x_j)} u_{ik}, \quad (11) $$

where $\mathcal{N}(x_j)$ represents a square neighborhood (in the spatial domain) around the pixel $x_j$. The term $h_{ij}$ represents the probability that the pixel $x_j$ belongs to i-th cluster: as a result, the spatial function of a pixel for a cluster is large if the majority of its neighbourhood belongs to the same cluster. The contribution of the spatial function modifies the classic FCM membership function according to Equation (12):

$$ u'_{ij} = \frac{u_{ij}^p h_{ij}^q}{\sum_{k=1}^{C} u_{kj}^p h_{kj}^q}, \quad (12) $$

where p and q control the relative importance of both functions. In a homogeneous region, the spatial functions simply emphasize the original membership and the clustering remains
unchanged. On the other hand, for a noisy pixel this formula reduces the weighting of a noisy cluster by the labels of its neighboring pixels: as a result, misclassified pixels from noisy regions can easily be corrected. The sFCM clustering process involves two steps at each iteration: the former, which is the same as that in standard FCM, allows for calculating the membership function in the features domain, while the latter maps the membership information of each pixel into the spatial domain and allows to calculate the spatial functions for each pixel of the image. At this point, the FCM iteration proceeds with the new membership that is incorporated with the spatial function and stops when the maximum difference between cluster centers at two successive iterations is lower than a certain threshold. After the convergence, a defuzzification scheme (i.e., maximum membership) is applied.

### 3. Materials and Methods

As introduced in Section 1, the purpose of this work is to present an in-depth analysis of classical unsupervised segmentation algorithms, developing them appropriately to adapt them to the clinical case addressed (i.e., breast lesions detection), thus improving the overall performance.

This section, along with showing the characteristics of the DCE-MRI dataset analyzed, describes the processing pipeline of the proposed analysis. It is necessary to point out that particular attention was paid to the optimization of each step for the specific case study, basically on three different levels:

* *pre-processing*: obtaining images with similar characteristics for the downstream processing steps;

* *segmentation*: optimization of the parameters of the investigated segmentation algorithms;

* *post-processing*: definition of appropriate region properties (based on connected-components) to refine the segmentation results.

#### 3.1. MRI Dataset Description

The analysis of this study was performed on a clinical DCE-MRI dataset composed of 50 patients with breast cancer: a total of 599 slices were processed. The main details on MRI acquisition parameters are reported in Table 1, while in Figure 1 the phases related to some benign and malignant lesions are shown. The dataset includes patients with different stages of breast cancer, allowing us to cover a wide clinical scenario: as matter of fact, this dataset contains various levels of segmentation difficulty with some scans showing low contrast and large inhomogeneities. Lesions with a non-homogeneous enhancement region, masses with an irregular shape or necrotic core are also included.

<table><caption>**Table 1.** Some characteristics of the DCE-MRI dataset used for this work.</caption><thead><tr><td>Characteristic</td><td>Value</td></tr></thead><tbody><tr><td>Scanner manufacturer</td><td>General Electric (GE)</td></tr><tr><td>Scanner model</td><td>Signa HDxt</td></tr><tr><td>MRI sequence</td><td>DCE-MRI</td></tr><tr><td>Time Repetition [ms]</td><td>37.720</td></tr><tr><td>Time Echo [ms]</td><td>17.640</td></tr><tr><td>Slice thickness [mm]</td><td>2</td></tr><tr><td>Interslice spacing [mm]</td><td>1</td></tr><tr><td>Slice pixel spacing [mm]</td><td>0.6875</td></tr><tr><td>Matrix size [pixels]</td><td>512 × 512</td></tr></tbody></table>
<figure><img src="image_37.png" /><figcaption>**Figure 1.** Examples of DCE-MRI phases related to four breast cancer patients (one per row): malignant lesions are shown on rows 1 and 3, while benign lesions are shown on rows 2 and 4. Among phases, for segmentation purposes the clinician selected the strongest one. In particular, for these lesions, the 4th, 3rd, 4th, and 3rd phase was chosen, respectively.</figcaption></figure>

## 3.2. *The Proposed Analysis*

Clinical knowledge points out that breast lesions in DCE-MRI appear hyperintense compared to the adipose and muscle surrounding tissues. Furthermore, anatomical atlases and diagnosis reports also refer that breast lesions have the tendency to grow in an isotropic manner, preserving a pseudo-spherical shape [39]. Relying on these morphological features, in this study we decided to start the segmentation from the central slice, which should be the one with the largest area: adjacent slices were processed successively. The idea behind this choice is to segment initially the slice in which the lesion appears more evident, finding the centroid of the connected-component identifying the lesion. On slices different from the central one, the lesion ROI was found by identifying the connected-component with the centroid that has the minimum distance from the centroid found for the central slice. The whole segmentation process can be summarized as follows:

* *dataset loading:* all DCE-MRI images (belonging to the selected patient) and its ground-truth are loaded;

* *ROI selection:* a bounding box containing the lesion is manually drawn in the central image: in this way, the algorithm analyzes only the pixels within the box, reducing computational times and avoiding classification mismatches;

* *pre-processing:* necessary to reduce noise and provide images with similar characteristics for the next processing steps;

* *lesion segmentation:* once the images are pre-processed, the MR image stack of each patient is segmented by all the analyzed methods (i.e., SMRG, k-means, FCM, and sFCM);

* *post-processing:* performed to refine the segmentation and to properly identify the ROI from each binary mask obtained in the previous step;
* performance evaluation: the final masks are stored and performance metrics (i.e., area-based and distance-based metrics) are calculated comparing the masks against the ground-truth.

In addition to the well-known segmentation techniques, also the exploited pre- and post-processing operations are widely used by the international community in the field of medical imaging. In fact, the literature offers many works that use the same pre-processing [42,43] and post-processing [44,45] steps in the image analysis phases (before segmentation), as well as in the refinement phases (after segmentation), with excellent results.

Figure 2 shows the flow diagram of the implemented and proposed analysis. Each processing block is described in the following subsections. All the methods were developed using the MatLab environment (Natick, MA, USA). The code is available via GitHub: https://github.com/carmilitello/UnsupervisedSegmentation.git (accessed on 9 November 2021).

### 3.2.1. Dataset Loading and ROI Selection

This step selects the MRI series to be analyzed. The data loading step, even though does not represent a real processing step, is essential for each algorithm—both supervised and unsupervised—to obtain the data to process. Moreover, the concept of ‘supervision’ refers generally to the training phase of an algorithm by requiring labeled data. For these reasons, the ‘unsupervised’ nature of the algorithm used here is preserved. In order to reduce computational times and improve segmentation performance, a minimal user interaction is needed. In fact, by manually tracing a rectangular bounding-box on the image, the operator provides to approximately select the initial ROI containing the breast lesion [13]. Once traced, this ROI is used to crop all the slices of the MR image stack

The manual selection of the initial ROI containing the lesion is useful to reduce the processing time (since only a portion of the whole image is processed), as well as to exclude regions containing pixels with characteristics similar to the lesion that could complicate the algorithm performance, thus invalidating the final result. Moreover, the initial, interactive input provides the clinician with the ability of confidently controlling the entire segmentation process, which is generally preferred compared to a fully automatic process [18,46]. Moreover, considering that the method works on an ROI—selected by the clinician—it is assumed that the tumor is always present within the ROI.
<figure><img src="image_38.png" /><figcaption>**Figure 2.** Flow diagram of the proposed analysis for breast lesion segmentation.</figcaption></figure>

### 3.2.2. Pre-Processing

In order to improve segmentation results, some pre-processing operations are applied after the ROI selection (Figure 3). In particular, the aim of this preliminary phase is to perform denoising and data pre-processing, allowing us to achieve reliable results during the segmentation phase.

The first pre-processing operation deals with noise reduction [47]. The use of median or average (with Gaussian or flat kernels) filters is a common choice in MR image processing. For instance, median filtering was applied to facilitate the ROI identification by reducing the outlier introduced by anatomical peculiarities [48] or also deal with small patient shifts [49]. We used Gaussian kernels having the form in Equation (13), which are the only circular
symmetric kernels that are also separable and that allow you to reduce noise by altering the image less than the average filter.

$$G(r) = ke^{-\frac{r^2}{2\sigma^2}}, \qquad (13)$$

where r, σ and k represent the radius, the standard deviation and the normalization factor of the Gaussian function G, respectively. Kernel normalization, obtained by multiplying its coefficients by k—obtained as the inverse of the sum of all kernel coefficients—has two purposes: (i) the average value of an area of constant intensity would equal that intensity in the filtered image (as it should) and (ii) it prevents the introduction of biases during filtering (i.e., the sum of the pixels in the original and filtered images will be the same). For the purpose of this work, a Gaussian kernel with r = 5 and σ = 0.6 was chosen: this standard deviation value performed denoising without an excessive blurring of the original image.

<figure><img src="image_39.png" /><figcaption>Figure 3. Flow diagram of the pre-processing steps approach.</figcaption></figure>

The next pre-processing step subtracts the mean value from the original image. This is a very helpful step in signal processing because it emphasizes signal variations by shifting its mean to zero. For this reason, during the pre-processing phase, the algorithm removed this bias from the original image amplifying gray level variations. Once the mean value subtraction, the pixels resulting in a negative value have to be clipped to zero in order to avoid visualization problems in the following steps.

After noise reduction and mean subtraction, the last step performs a modified gamma-transformation to stretch the original image histogram. The general form of a gamma transformation is defined in Equation (14):

$$s = cr^{\gamma}, \qquad (14)$$

where c and γ are positive constants. Power-law curves with fractional values of γ map a narrow range of dark input values into a wider range of output values, and vice versa for bright input values. When c = γ = 1 the gamma transformation reduces to the identity transformation. The modified gamma transformation proposed in this work used c = 1 and γ = 0.7, by means a piece-wise function in which the final value of each pixel depends on a certain threshold value θ. In particular, the function mapping the initial value of a pixel r(x, y) into its final value s(x, y) is defined in Equation (15). The value of θ, used for the piece wise gamma selection, is the one suggested by the Otsu's thresholding method [30]. Figure 4 shows the results of each pre-processing step applied on two (one per row) DCE-MRI breast lesions.

$$s(x,y) = \begin{cases} r(x,y) & \text{if } r(x,y) \le \theta \\ r(x,y)^{\gamma} & \text{if } r(x,y) > \theta \end{cases}. \qquad (15)$$
<figure><img src="image_40.png" /><figcaption>**Figure 4.** Examples of pre-processing steps on two breast lesions (one per row): (1st column) DCE-MRI breast image (after cropping); (2nd column) Gaussian smoothing result; (3rd column) centering result; and (4th column) gamma correction result. With respect to the position in the DCE-MRI sequence, all images have been rotated 90° clockwise, in order to improve the graphic representation and allow all the images relating to a specific lesion to be displayed in a single line.</figcaption></figure>

### 3.2.3. Segmentation

After pre-processing, the segmentation involves the partitioning of a digital image into multiple sets of pixels, according to the specific clinical purpose. Nowadays, segmentation of digital images still remains one of the most challenging topics in image processing: in fact, medical image segmentation is often still performed manually, via time-consuming and operator-dependent procedures.

As indicated in Section 3.2.1, to start the segmentation process, it is necessary to select the ROI containing the lesion in the central image. The central slice—which is generally the one with the largest tumor section—is determined automatically after the operator sets the range of the initial and final slices containing the tumor. Starting from the selected ROI in the central slice, the ROIs in the adjacent slices are determined starting from the centroid of the segmented lesion in the previous slide. By doing so, the ROIs (set automatically) in the adjacent slices are centered (or almost) with the lesion section they contain. After the segmentation, a check is made to verify that the lesion does not touch the ROI boundary. In fact—considering that the ROI size is set equal to that of the ROI in the previous slice—it could happen that the selected ROI is too small and, consequently, the lesion was cropped. If so, the ROI is enlarged in order to fully contain the lesion section and repeated segmentation.

The segmentation is sequential: starting from the central slice, all the upper slices are processed first and then the lower ones. All the segmentation steps are the same for both the central slice and the other slices. The only difference lies in the different definition of the ratios—defined in Equations (18) and (19)—used in the post-processing steps, differentiated with the goal of optimizing the segmentation result.

#### SMRG Setting

There are several parameters controlling the SMRG behavior.

First of all, the splitting predicate P, a logical predicate fundamental to achieve satisfactory outputs. In order to identify breast lesions, homogeneity criteria are defined in terms of the mean value $\mu$ of each quad-regions. Here, the logical predicate in Equation (16) was used:
$$ P := \begin{cases} \text{TRUE} & \text{if } (\mu > 0.1) \land (\mu < 0.6) \\ \text{FALSE} & \text{otherwise} \end{cases}. \qquad (16) $$

The minimum block dimension $\rho_{min}$ sets the minimum quad-region size beyond which no further splitting is carried out: the best results are found with minimum block dimensions of 4 × 4 pixels, because small regions are detected too. After the splitting phase, the final partition will contain adjacent regions with identical properties: this drawback can be addressed by merging only adjacent regions whose combined pixels satisfy the predicate P.

Finally, a stopping rule interrupts the growing procedure if no more pixels match the membership criterion. For the aim of this work, a reasonable stopping rule uses a similarity criterion between the candidate pixel to be incorporated and the pixels already belonging to the region. The criterion for the stopping rule is defined in terms of absolute distance between the regional mean of each quad-region and the threshold provided by the Otsu's method [30] on the original cropped image, both of them calculated without the contribution of null pixels.

## Clustering Setting

Here, the optimal setting used for k-means, FCM and sFCM clustering algorithms is reported, in terms of the following parameters: (i) number of clusters; (ii) maximum number of iterations; (iii) minimum improvement in objective function between two consecutive iterations; (iv) exponent for the fuzzy partition matrix (only for fuzzy approaches); (v) $<p, q>$ parameters (only for fuzzy approaches).

Regarding the number of clusters, the lower the partition fuzziness, the better the segmentation result: the best clustering is achieved when $V_{pe}$ is minimal. If the partition entropy for $C = 2$ is lower than the one obtained for $C = 3$, the lesion does not include any necrotic region: if so, pixels belonging the brightest cluster are turned to 1 and the others to 0. Otherwise, the two brightest clusters are fused together and their pixels are turned to 1, allowing for the inclusion of necrotic cores into the preliminary mask. Considering that, sometimes breast lesions can include an inner necrotic region that appears darker (i.e., hypo-intense) with respect to the rest of the lesion. In those cases where a binary clustering ($C = 2$) does not allow for properly detecting a necrosis, the number of clusters might be increased ($C = 3$). The choice of the optimal number of clusters to be used in k-means, FCM and sFCM is automatically set by evaluating the partition entropy for both cases $C = 2$ and $C = 3$. The partition entropy $V_{pe}$ is a cluster validity function defined as in Equation (17):

$$ V_{pe} = \frac{\sum_{i=1}^{N} \sum_{j=1}^{C} u_{ij} \log u_{ij}^{2}}{N}. \qquad (17) $$

The $<p, q>$ parameters, dealing with FCM and sFCM clustering, control the relative importance of both functions defined in Equations (10) and (11), respectively. In particular, with $<p, q> = <1, 1>$, we sFCM is equivalent to the traditional FCM. As a matter of fact, to properly weight spatial information, the values $<p, q> = <1, 2>$ were used. The choice was guided by the analysis and results obtained in [18], tackling a similar problem on DCE-MRI images. All these parameters are reported in Table 2.

At the end of the segmentation of each slice, the algorithm also verifies whether or not the contour of the detected lesion touches the boundaries of the initial rectangular bounding-box: if so, the shape of the lesion results in a cut version of the original one. For this reason, when the lesion intersects the boundaries of the original $m \times n$ rectangular crop, the algorithm shifts and expands the bounding-box itself until the boundaries of both the lesion and the crop do not intersect each other.
<table><caption>**Table 2.** Parameters setting for k-means, FCM and sFCM clustering algorithms.</caption><thead><tr><th>Unsupervised Algorithm</th><th>Number of Clusters</th><th>Maximum Number of Iterations</th><th>Minimum Tolerance</th><th>Exponent of Fuzzy Partition Matrix</th><th>⟨p, q⟩ Parameters</th></tr></thead><tbody><tr><td>k-means</td><td>2 or 3</td><td>100</td><td>N.A.</td><td>N.A.</td><td>N.A.</td></tr><tr><td>FCM</td><td>2 or 3</td><td>100</td><td>1 × 10<sup>-5</sup></td><td>2.0</td><td>&lt;1,1&gt;</td></tr><tr><td>sFCM</td><td>2 or 3</td><td>100</td><td>1 × 10<sup>-2</sup></td><td>2.0</td><td>&lt;1,2&gt;</td></tr></tbody></table>

### 3.2.4. Post-Processing

The result of the segmentation process is a binary image in which a label is assigned to each pixel: 0 if the pixel does not belong to the mask, and 1 otherwise. The aim of the post-processing phase is to properly choose, among the connected-components resulting from the previous step, the only one representing the lesion. The post-processing itself consists of several sub-steps allowing for removing connected-components that do not meet specific morphological criteria (Figure 5).

<figure><img src="image_41.png" /><figcaption>**Figure 5.** Diagram of the post-processing steps for both central slice and other slices.</figcaption></figure>

The first post-processing phase allows us to extract all the connected-components from the raw binary mask. For each connected-component, a list of morphological features (i.e., extreme points, area, centroid, eccentricity, mean intensity and solidity) are computed that the algorithm is going to exploit in the following steps. In order to remove spurious connected-components with an area that is too small to be considered as a lesion, an area-based selection is performed: all the connected-components with an area smaller than 0.75 × meanArea of all the connected-components. In such a way, spurious regions with a very small area are removed, reducing the number of valid candidates to be analyzed in the following steps.
Clinical evidence allows us to state that breast lesions generally grows in an isotropic way preserving a pseudo-spherical appearance [39]. As a consequence, their particular shape can be described in terms of eccentricity, which is defined as the ratio of the distance between the foci of an ellipse and its major axis length. Eccentricity values lay in the range [0, 1]: the extreme values represent degenerate cases identifying a circle and a line segment, respectively. The aim of this post-processing step is to delete from the list of the lesion candidates the regions with eccentricity lower than 0.3 and higher than 0.9 (these threshold values were determined experimentally). Eccentricity values close to 1 relate to very elongated lesions, while eccentricity values close to 0 relate to almost perfectly round shapes. In the case of breast lesions, even if overall there are rounded lesions, the lesion has a 'lobed' trend for which the final value of the eccentricity never takes values below 0.2–0.4. Using 0.3 as eccentricity lower limit in the post-processing phase, allow us to eliminate those (almost perfectly) circular connected components with pixel values similar to the lesion (and which are therefore incorrectly selected), but which are instead part of the background. As a consequence of their tendency to be round-like shaped, breast lesions exhibit an high solidity, which is defined as the ratio between the region area and the (including the region) convex polygon area. Considering this morphological information, the lesion ROI identification can be turned into an optimization problem with the ultimate goal of maximizing ratioCentral, defined in Equation (18):

$$ \text{ratioCentral}_i = \frac{\mu_i s_i}{\epsilon_i}, \qquad (18) $$

where $\mu_i$ is the mean intensity, $s_i$ is the solidity and $\epsilon_i$ denote the eccentricity of the i-th connected-component, respectively.

The post-processing in the other slices differs from the central slice. In this case, the parameter to maximize ratioOther is defined according to Equation (19):

$$ \text{ratioOther}_i = \frac{\mu_i s_i}{d_i^2}, \qquad (19) $$

where $d_i^2$ represents the square of the distance between the centroid of the lesion region in the central slice and the centroids of all the connected-components in the current slice. In order to reduce false positives identification, ratioOther maximization searches for the connected-component with the centroid that is the closest to the one of the central slice mask. Recalling that breast lesions usually preserve a pseudo-spherical appearance, as the slices move away from the central one, the cross section identifying the lesion is reduced. In Equation (19) the replacement of $\epsilon_i$ with $d_i^2$ aims to avoid misclassification of spurious regions in slices different from the central one.

Figure 6 shows the results of each post-processing step applied to a segmentation mask obtained by means of sFCM clustering. Unfortunately, breast lesions might be characterized by different scenario—in terms of uniformity, contrast and well-defined boundaries—which a simple thresholding cannot properly manage. Figure 7 shows two segmentation results that allow us to appreciate the lesion non-homogeneity that only clustering approaches can properly manage, thus maximizing the result accuracy.
<figure><img src="image_42.png" /><figcaption>**Figure 6.** Examples of post-processing steps on two breast lesions (one per row): (1st column) binary mask obtained after sFCM clustering; (2nd column) mask after the small-area removal; (3rd column) mask after both the small-area removal and the eccentricity-based selection; and (4th column) mask after ratio criteria, defined in Equations (18) and (19). With respect to the position in the DCE-MRI sequence, all images have been rotated 90° clockwise, in order to improve the graphical representation and allow all the images, related to aspecific lesions, to be displayed in a single line.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Figure 7.** Two examples of segmentation obtained by sFCM: (1st row) the original DCE-MRI images; (2nd row) the corresponding segmentation. Compared to the images shown in Figures 4 and 6, a 2× zoom factor was used, that allows us to better appreciate the lesion non-homogeneity that only clustering approaches can properly manage, thus maximizing the final result accuracy.</figcaption></figure>
## 4. Segmentation Performance Evaluation

Several evaluation measures are computed to quantify segmentation performances by comparing masks—obtained from analyzed unsupervised methods—and the ground-truth—provided by a radiologist, with more than 5-year experience on breast MRI, in consensus with a consultant breast radiologist (with more than 30-year experience on breast imaging). To obtain an accurate and detailed quantification, both area-based and distance-based metrics were used. The reason behind this choice is that area-based metrics strongly depend on region size and are not always able to evaluate the precision of a segmentation approach. On the other hand, distance-based metrics take into account the distance between the boundaries of the two segmentations to be compared, ignoring the actual volume difference between the two masks.

### 4.1. Spatial Area-Based Metrics

Spatial area-based metrics compare the semi-automatic segmented regions with the manually segmented ones ($R_A$ and $R_T$, respectively) by calculating the overlapping percentage of area between the two masks obtained from the segmentation of the image $I$. Recalling some basics on statistical decision theory measures, the regions containing ‘true positives’ (TP), ‘false positive’ (FP), ‘false negatives’ (FN), and ‘true negatives’ (TN) are defined as:

$$ 
\begin{aligned}
R_{TP} &= R_A \cap R_T \\
R_{FP} &= R_A - R_T \\
R_{FN} &= R_T - R_A \\
R_{TN} &= I - R_T - R_A
\end{aligned}
 $$

When validating the segmentation results, the two most used area-based metrics are the Dice Index (DI) and the Jaccard Index (JI), defined in Equations (20) and (21), respectively. DI and JI are used to describe how much similar the manual (ground-truth) and the semi-automatic segmentations are: the greater they are, the higher is the overlapping percentage between the two masks.

$$ DI = \frac{2 \cdot R_{TP}}{R_A + R_T} \qquad (20) $$

$$ JI = \frac{R_A \cap R_T}{R_A \cup R_T} \qquad (21) $$

Sensitivity and specificity—defined in Equations (22) and (23)—represent the portion of positive pixels (foreground) and negative pixels (background) correctly detected by a segmentation method with respect to the ground-truth, respectively.

$$ Sensitivity = \frac{R_{TP}}{R_{TP} + R_{FN}} \qquad (22) $$

$$ Specificity = \frac{R_{TN}}{R_{TN} + R_{FP}} \qquad (23) $$

False Positive Ratio (FPR) and False Negative Ratio (FNR)—defined in Equations (24) and (25)—denote the presence of false positives and false negative compared to the reference region, respectively.

$$ FPR = \frac{R_{FP}}{R_{FP} + R_{TN}}, \qquad (24) $$
$$ FNR = \frac{R_{FN}}{R_{FN} + R_{TP}}. \qquad (25) $$

## 4.2. Spatial Distance-Based Metrics

Area-based metrics are susceptible to differences between the positions of segmented regions and strongly dependent on their own size. To take into account the spatial position of the pixels, it is necessary to quantify the distance between the boundaries computed by the semi-automatic methods and the ground-truth delineated by the expert. Let $A = \{a_i : i = 1, 2, \dots, K\}$ be the set of vertices belonging to the semi-automatic mask and $T = \{t_i : i = 1, 2, \dots, N\}$ the set of vertices belonging to the ground-truth, the distance between the i-th pixel in A and the set T is defined as:

$$ d(a_i, T) = \min_{j \in \{1, 2, \dots, N\}} \|a_i - t_j\|_2, \qquad (26) $$

where $||a_i - t_j||_2$ denotes the Euclidean distance between two points.

Many metrics can be defined in order to quantify the similarity/dissimilarity between two segmentations.

The Mean Absolute Distance (MAD)—defined in Equation (27)—quantifies the average error in the segmentation process. The Maximum Distance (MaxD)—defined in Equation (28)—measures the maximum difference between the two ROI boundaries. The Hausdorff Distance (HD) between the point sets A and T—defined in Equation (29)—measures the maximal distance from a point in the first set to a nearest point in the other one.

$$ MAD = \frac{1}{K} \sum_{i=1}^{K} d(a_i, T) \qquad (27) $$

$$ MaxD = \max_{i \in 1, 2, \dots, K} \{d(a_i, T)\} \qquad (28) $$

$$ HD = \max \{h(T, A), h(A, T)\}, \qquad (29) $$

where $h(T, A) = \max_{t \in T} \{\min_{a \in A} \{d(t, a)\}\}$ is the so-called 'directed Hausdorff Distance'.

It is important to point out that all the measured distances are expressed in pixels: in this way, they result will be independent from the spatial resolution among different MRI datasets (i.e., pixel spacing).

# 5. Experimental Findings

## 5.1. Area-Based Metrics Segmentation Results

Area-based metrics obtained by each segmentation algorithm are shown in Tables 3 and 4—expressed as mean ± standard deviation. As easily appreciable, the results showed that the fuzzy framework offered by FCM and sFCM reflects the intrinsic uncertainty that characterizes medical images, allowing us to achieve better segmentation results compared to the hard clustering performed by k-means. Furthermore, it is worth to note that spatial constraints taken into account by sFCM contribute to reduce the standard deviation of the final result, thus ensuring higher reliability with respect to SMRG and k-means. On the other hand, the large value in standard deviation indicates that SMRG has a high variability of the results that affects its reliability. The boxplots in Figures 8 and 9 summarize the obtained results.

Regarding sensitivity, the results showed that the FCM-based approach offered better performance in terms of both mean value and standard deviation; on the contrary, the hard partition offered by k-means did not provide satisfying results. The specificity values showed that the clustering-based approaches obtained better performances compared to SMRG. In particular, sFCM and k-means achieved slightly better results than FCM.
Altogether, the results in Tables 3 and 4 showed that segmentation approaches based on soft-clustering techniques achieved better performance compared to SMRG and k-means. In fact, by explicitly exploiting the fuzziness, both FCM and sFCM better handled the intrinsic uncertainty and the natural variability of medical images. On the other hand, the crisp k-means allowed us to reach satisfying results in specificity values, but without granting good performance in sensitivity.

<figure><img src="image_44.png" /><figcaption>**Figure 8.** Boxplots of the area-based metrics achieved by the considered unsupervised segmentation approaches. From left to right: Dice Index, Jaccard Index, Sensitivity and Specificity. The lower and upper bounds of each box represent the first and the third quartiles of the metric distribution, respectively. The median is represented by a red line, while outliers are displayed as red crosses.</figcaption></figure>

<table><caption>**Table 3.** Area-based metrics achieved by the considered unsupervised segmentation approaches: the results are expressed as average value ± standard deviation.</caption><thead><tr><th>Method</th><th>DI</th><th>JI</th><th>Sensitivity</th><th>Specificity</th></tr></thead><tbody><tr><td>SMRG</td><td>71.82 ± 10.83</td><td>58.98 ± 11.32</td><td>76.47 ± 12.41</td><td>73.65 ± 13.50</td></tr><tr><td>k-means</td><td>72.65 ± 8.26</td><td>59.87 ± 9.20</td><td>67.23 ± 10.78</td><td>86.97 ± 8.93</td></tr><tr><td>FCM</td><td>77.48 ± 6.77</td><td>65.22 ± 8.13</td><td><strong>77.84 ± 8.72</strong></td><td>81.94 ± 10.41</td></tr><tr><td>sFCM</td><td><strong>78.23 ± 6.50</strong></td><td><strong>65.90 ± 8.14</strong></td><td>74.69 ± 9.39</td><td><strong>87.10 ± 8.24</strong></td></tr></tbody></table>
<figure><img src="image_45.png" /><figcaption>**Figure 9.** False Positive Ratio and False Negative Ratio boxplots. The lower and upper bounds of each box represent the first and the third quartiles of the metric distribution, respectively. The median is represented by a red line, while outliers are displayed as red crosses.</figcaption></figure>

<table><caption>**Table 4.** False Positive and False Negative Ratios achieved by the considered unsupervised segmentation approaches: the results are expressed as average value ± standard deviation.</caption><thead><tr><td>Method</td><td>False Positive Ratio</td><td>False Negative Ratio</td></tr></thead><tbody><tr><td>SMRG</td><td>0.55 ± 1.05</td><td>0.33 ± 0.11</td></tr><tr><td>k-means</td><td>0.25 ± 0.20</td><td>0.25 ± 0.09</td></tr><tr><td>FCM</td><td>0.16 ± 0.13</td><td>0.23 ± 0.12</td></tr><tr><td>sFCM</td><td><b>0.14 ± 0.12</b></td><td><b>0.22 ± 0.09</b></td></tr></tbody></table>

Figure 10 shows some segmentation examples focusing on scenarios where the results are not particularly satisfactory in terms of FPs and FNs. In general, it is possible to note that SMRG and k-means have a greater tendency to leave out parts of the lesion (FN)—especially, when these are not uniform—and to include areas outside the lesion (FP), while the techniques based on clustering, and in particular the sFCM guarantees a better lesion detection.
<figure><img src="image_46.png" /><figcaption>**Figure 10.** Examples of interesting segmentation results showing FPs and FNs. In particular, the examples compare results yielded by the four investigated unsupervised method—namely, SMRG (cyan), k-means (yellow), FCM (red), sFCM (blue)—against the ground-truth (green).</figcaption></figure>
5.2. *Spatial Distance-Based Metrics Segmentation Results*

The need for using both area- and distance-based metrics comes out in considering
that area-based metrics do not take in account pixels' spatial distribution. This leads to the
need to quantify the distance between the boundaries computed by the semi-automatic
methods and the ground-truth. The boxplots in Figure 11 summarize the obtained results.

In terms of HD, all the presented methods share similar characteristics with a mean
value of ≈2.2 and a standard deviation of ≈0.43. This means that the boundaries of the
semi-automatic masks are quite close to the manually traced ones.

The MAD metric is quite similar across the examined techniques, except for the SMRG,
which shows the higher mean value and standard deviation. Clustering-based approaches,
on the other hand, result more stable and precise, offering a better performance on the
whole dataset.

<figure><img src="image_47.png" /><figcaption>**Figure 11.** Boxplots of the spatial distance-based metrics achieved by the considered unsupervised segmentation approaches. From left to right: Mean Absolute Distance, Maximum Distance and Hausdorff Distance. The lower and upper bounds of each box represent the first and the third quartiles of the metric distribution, respectively. The median is represented by a red line, while outliers are displayed as red crosses.</figcaption></figure>

Table 5 shows distance-based metrics obtained using the considered unsupervised
segmentation approaches: lower distance values indicate better segmentation results.
Observing the general trend, just a small deviation between the segmentations of the pro-
posed methods and those of the experienced radiologist can be denoted. Furthermore, the
achieved spatial distance-based indices are consistent with area-based metrics, confirming
that clustering-based segmentation approaches allow us to reach better results with respect
to SMRG.
<table><caption>**Table 5.** Spatial distance-based metrics achieved by the proposed segmentation approaches: the results are expressed as average value ± standard deviation.</caption><thead><tr><th>Method</th><th>MAD</th><th>MaxD</th><th>HD</th></tr></thead><tbody><tr><td>SMRG</td><td>2.00 ± 1.41</td><td>5.63 ± 3.81</td><td>2.32 ± 0.44</td></tr><tr><td>k-means</td><td>1.58 ± 0.99</td><td>4.22 ± 2.73</td><td>2.24 ± 0.42</td></tr><tr><td>FCM</td><td>1.58 ± 1.22</td><td>4.46 ± 3.12</td><td><strong>2.21 ± 0.43</strong></td></tr><tr><td>sFCM</td><td><strong>1.37 ± 0.90</strong></td><td><strong>4.04 ± 2.87</strong></td><td>2.21 ± 0.44</td></tr></tbody></table>

As shown at the top of Figure 12, FCM (with $p = 6.268 \times 10^{-6}$ and $p = 2.670 \times 10^{-7}$) and sFCM (with $p = 8.566 \times 10^{-5}$ and $p = 1.855 \times 10^{-7}$) clustering methods achieved significantly higher DI values compared to SMRG and k-means, respectively.

<figure><img src="image_48.png" /><figcaption>**Figure 12.** DI values achieved by the investigated traditional classic approaches. The bar graph and error bars denote the average value and the standard deviation DI values, respectively. The *p*-values, obtained from the statistical validation procedure, are shown at the top of the bars as brackets. The Wilcoxon rank-sum test for pairwise result comparison was used with the alternative hypothesis that the samples do not have equal medians of DI values. A significance level of $\alpha = 0.05$ with a correction using the Bonferroni-Holm method. Notation: \*\*\*\* $p < 0.0001$.</figcaption></figure>

To statistically validate the obtained results, the two-sided Wilcoxon signed rank test [50] on paired DI results was performed with the null hypothesis that the samples come from continuous distributions with equal medians (considering a significance level of 0.05). Obtained *p*-values are shown in Table 6. With more details, this test on paired results was used to statistically compare the distributions of the DI values achieved by two competing methods and identify significant differences. The *p*-values were corrected by the Bonferroni-Holm method [51] for multiple comparisons.

The values reported in Tables 3–5 represent the average value ± standard deviation over all 50 breast masses. All the values obtained on each breast mass are reported in the Supplementary Materials.

<table><caption>**Table 6.** *p*-Values obtained from the statistical validation procedure. The Wilcoxon rank-sum test for pairwise result comparison was used with the alternative hypothesis that the samples do not have equal medians of DI. A significance level of $\alpha = 0.05$ with the Bonferroni-Holm correction for multiple comparisons was used. **Boldface** indicates that the null hypothesis can be rejected.</caption><thead><tr><th></th><th>SMRG</th><th>k-Means</th><th>FCM</th></tr></thead><tbody><tr><td>k-means</td><td>0.973</td><td></td><td></td></tr><tr><td>FCM</td><td><strong>6.268 × 10<sup>-6</sup></strong></td><td><strong>2.670 × 10<sup>-7</sup></strong></td><td></td></tr><tr><td>sFCM</td><td><strong>8.566 × 10<sup>-5</sup></strong></td><td><strong>1.855 × 10<sup>-7</sup></strong></td><td>0.423</td></tr></tbody></table>
### 5.3. *Processing Times*

In order to evaluate the processing times, all segmentations were performed on the entire dataset, by calculating the average elapsed time and the corresponding standard deviation for each of the four investigated algorithms over all the analyzed images, obtaining the following values: SMRG: 2.78 ± 2.79 s; k-means: 1.35 ± 0.78 s; FCM: 1.65 ± 0.88 s; sFCM: 1.76 ± 0.83 s. These processing times were measured using the Matlab R2019b IDE (by means of the tic and toc stopwatch timer functions) running on a Windows 10 Pro general-purpose PC equipped with an Intel I7-3630QM@2.40 GHz CPU and 8 GB RAM.

As expected, SMRG had the longest processing time due to the two-stage approach and also considering the iterative processes employed during the Split-and-Merge and Region Growing executions. Interestingly, the introduction of the fuzzy logic is negligible compared to the crisp k-means, as well as the integration of the spatial constraints into the sFCM algorithm does not require a remarkable computational overhead in addition to the standard FCM clustering.

Overall, these results demonstrate the clinical feasibility of the investigated classic unsupervised methods also in terms of both computational resources and processing times, by considering that supervised CNN-based approaches for segmentation require a training phase and then an inference phase typically performed on Graphics Processing Units [18].

### 5.4. *Difficult Cases*

As previously highlighted, medical images are characterized by an intrinsic variability in which boundaries or anatomical details may be not well defined. Furthermore, noise corrupts digital images, thus affecting certain features within the original image. MRI suffers from various kinds of noise and artefacts because of the nature of the signal detection and spatial encoding [37,38]. For instance, hardware-induced errors are often caused by the complicated acquisition scheme depending on radiofrequency coils, while thermal noise can derive from transmission lines, receiver circuits and polarization magnetic field $B_0$ drift during the scan acquisitions. In addition, natural body motion (e.g., respiratory and cardiac motion) can degrade the image quality too. As a consequence, even after a proper pre-processing step, it is common to deal with low-contrast, noisy images. The proper lesion segmentation in this kind of images is not a trivial task and user interaction would be required to produce accurate results.

#### 5.4.1. Case with Low Contrast-Enhanced Mass

In some cases, MR images are difficult to segment because the lesion itself does not appear brighter with respect to the muscle and adipose surrounding tissue. This kind of images exhibit a narrow histogram located typically toward the middle of the intensity scale, implying a washed-out grey look through the whole image. As a consequence, the meaningful partition of the original image results in a very difficult task leading to imprecise results. The opposite is true for the histogram of a high-contrast image, which covers a wide range of the intensity scale and has a pixel distribution not too far from uniform. The effect is an image that shows a great deal of gray-level detail and has a high dynamic range. As a matter of fact, a segmentation process on this kind of images will produce satisfying results with a high reproducibility. As easily appreciable in Figure 13, MR images offer a very difficult scenario where the lesion boundaries are not clearly distinguishable from the rest of the image. The manual segmentation (green contour) cannot properly segment the whole lesion because of the strong uncertainty due to the low contrast characterizing the whole image. On the other hand, the FCM segmentation (red contour) correctly identifies the lesion but, because of the low percentage of overlapping area with the manual mask, it does not ensure satisfying results.
<figure><img src="image_49.png" /><figcaption>**Figure 13.** Case with a low-contrast lesion: (a) entire axial DCE-MRI slice; (b) crop of the slice shown in (a). Comparison between manual ground-truth (green) and automatic segmentation (red).</figcaption></figure>

### 5.4.2. Case with Blurred Boundary Mass

From Figure 14, it is possible to appreciate how blur severely compromises the contrast of the original image, making lesion segmentation a very challenging task. As mentioned above, because of the strong uncertainty in boundaries delineation, the manual segmentation (green contour) identifies a very simple and smooth shape into which the lesion is included. Of course, this kind of strategy allows the identification of lesion's location, but also includes into the ROI a lot of false positives. On the other hand, the semi-automatic mask (red contour) identifies a more precise region avoiding the misclassification of FP pixels. Even in this case, because of the imperfect area, overlap area-based metrics will not yield a high score.

<figure><img src="image_50.png" /><figcaption>**Figure 14.** Case with a blurred boundary lesion: (a) entire axial DCE-MRI slice; (b) crop of the slice shown in (a). Comparison between manual ground-truth (green) and automated segmentation (red).</figcaption></figure>

### 5.4.3. Case with Irregular Mass

Opposed to what is reported in the literature [39], breast lesions may sometimes exhibit irregular shapes and borders with internal divisions. In these cases, breast masses could cause difficulties during the manual segmentation process: in fact, as the lesion contour becomes more irregular, the manual tracing of the ROI becomes more challenging. As a consequence, the manual segmentation of breast lesions with unusual shapes does not always match the effective lesion contour. Figure 15 exhibits one of this cases with unusual
elongated breast lesion. The manual segmentation (green contour) completely cut the left portion in the upper part of the lesion which is, instead, properly included in the semi-automatic boundary (red contour). As a consequence, computer-assisted segmentation process results in a more precise mask that properly reproduces the narrowed-shape of lesion.

<figure><img src="image_51.png" /><figcaption>**Figure 15.** Case with an irregular lesion: (a) entire axial DCE-MRI slice; (b) crop of the slice shown in (a). Comparison between manual ground-truth (green) and automated segmentation (red).</figcaption></figure>

**6. Discussion and Conclusions**

The main objective of this work was to offer a detailed analysis of well-established classical unsupervised segmentation techniques by carefully comparing them in a real clinical application. Breast cancer is the most common cause of cancer death in women worldwide [52,53] and the second most common cancer overall [54]. Fortunately, science evolution has led to the development of medical imaging techniques, which are used to detect abnormalities in breast parenchyma. Among imaging techniques, multiparametric MRI plays a crucial role and it is widely used in clinical applications, due to its high resolution images and the ability to precisely differentiate soft tissues.

As matter of fact, as a case study, the contrast-enhancing mass delineation in breast DCE-MRI was addressed by means of four popular unsupervised segmentation methods, namely: SMRG, k-means, FCM, and sFCM. Although they represent well-known approaches in the literature, they are still widely used in clinical tools. Starting from the basic versions of these approaches, during the initial analysis, we identified the shortcomings of each of them, developing and implementing improved versions, when possible.

Nowadays, deep learning techniques represent the state-of-the-art, allowing us to achieve high performance and accurate lesion segmentation for datasets of thousands of patients [14]. Deep learning approaches for image segmentation are generally supervised techniques that require a considerable computation times and a large amount of data for training [15]. In fact, these data must be representative of all the possible scenarios in which the deep neural network could operate, not always available in small- or medium-sized hospitals, and therefore not clinically feasible. Moreover, it should be noted that the authors of [18] showed that semi-automatic approached based on classic unsupervised techniques obtained results comparable or superior to the deep CNNs (namely, SegNet and U-Net). Therefore, this study was focused on classic pattern recognition approaches with the goal of providing an in-depth analysis.

It is important to point out that, even considering all the disadvantageous aspects related to manual segmentation, the contribution of an expert radiologist still remains essential at least to validate the results obtained by means a computer-assisted approach. In fact, clinicians rely on computational approaches with interpretable results [41]. From this point of view, the classical unsupervised approaches—such as those analyzed in this
work—provide this important advantage. This aspect is even more critical in deep learning architectures, where CNNs are generally adopted as ‘black-box’, thus making it difficult to offer a physical interpretation to the features encoded in all the intermediate CNN layers [46].

The obtained experimental results, in terms of area- and distance-based metrics, encourage the use of unsupervised pattern recognition techniques in medical image segmentation. In particular, consistently with [55], clustering-based segmentation approaches achieved better performance compared to the SMRG, the only ‘thresholding-based’ approach considered. As a consequence, crisp segmentation techniques—such as k-means and SMRG—are not well-suited for medical images that are characterized by an uncertain/variability (sometimes related to the noise), yielding inaccurate boundaries and not well-defined details. Both FCM and sFCM clustering techniques—implementing fuzzy modeling that provides an intrinsic flexibility—significantly achieved the best performance. In fact, on area-based metrics, they obtained DI = 78.23% ± 6.50 (sFCM), JI = 65.90% ± 8.14 (sFCM), sensitivity = 77.84% ± 8.72 (FCM), and specificity = 87.10% ± 8.24 (sFCM), FPR = 0.14 ± 0.12 (sFCM), and FNR = 0.22 ± 0.09 (sFCM). On distance-based metrics, they obtained MAD = 1.37 ± 0.90 (sFCM), MaxD = 4.04 ± 2.87 (sFCM), and HD = 2.21 ± 0.43 (FCM).

A second segmentation of the same radiologist or the segmentation of a different radiologist would have certainly allowed us to quantify the inter-/intra-operator variability of the results. Nevertheless, as already observed in [56], the mean DI was 0.81 (range 0.19–0.96). The mean DI is higher for the ‘easy tumors’ compared to the ‘challenging tumors’ (0.83 vs. 0.75, respectively, *p* < 0.001). The mean DI for each observer combination separately, for all tumors, ranged between 0.78 and 0.83, where the segmentations of the breast radiologist and the medical student showed the highest overlap. These results confirm that the performance achieved by the best performing methods are in line with the inter-observer agreement, also in terms of metrics variability according to the lesion types.

As further developments, we plan to investigate innovative improvements to further improve the performance with fuzzy clustering, by using (i) more sophisticated membership functions, and (ii) more advanced pre- and post-processing steps. Moreover, investigating and comparing the latest machine learning techniques, such as Generative Adversarial Networks (GANs), for unsupervised detection and segmentation [57,58] would be relevant with a sufficient amount of data for training and test. Finally, the implementation of multiparametric or multimodal approaches [59], by using different types of co-registered medical images—i.e., Diffusion Weighted Imaging (DWI) and Positron Emission Tomography (PET)/MRI—probably would allow us to improve the detection performance [60,61].

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/app12010162/s1: Table S1: Dice Index (DI) values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S2: Jaccard Index (JI) values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S3: Sensitivity values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S4: Specificity values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S5: False Positive Ratio (FPR) values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S6: False Negative Ratio (FNR) values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S7: Mean Absolute Distance (MAD) values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S8: Maximum
Distance (MaxD) values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported; Table S9: Hausdorff Distance (HD) values obtained by the four investigated unsupervised methods on each of the 50 segmented breast masses on DCE-MRI. In the last row the mean value ± the standard deviation is reported.

**Author Contributions:** Conceptualization, C.M. and A.R.; methodology, C.M. A.R. and L.R.; software, A.R.; validation, C.M. and L.R.; formal analysis, C.M. and L.R.; investigation, A.R.; resources, C.M., L.R. and I.D.; data curation, C.M. and I.D.; writing—original draft preparation, C.M., A.R. and L.R.; writing—review and editing, F.M., T.V.B., F.B. and G.R.; visualization, C.M. and L.R.; supervision, F.M., T.V.B., F.B. and G.R.; project administration, T.V.B., F.B. and G.R.; funding acquisition, L.R. and G.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study has received funding by GeSeTon project, funded by Italian MISE Grant No. 489 of 21 February 2018. This study has also been partially supported by The Mark Foundation for Cancer Research and Cancer Research UK Cambridge Centre [C9685/A25177] and by the Royal Society for the International Exchanges 2020 Cost Share with the Italian CNR (project No. IEC/R2/202313). Additional support was also provided by the National Institute of Health Research (NIHR) Cambridge Biomedical Research Centre [BRC-1215-20014]. The views expressed are those of the authors and not necessarily those of the NHS, the NIHR, or the Department of Health and Social Care.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of "Azienda Ospedaliera Universitaria Policlinico P.Giaccone" of Palermo, Italy (protocol code n.1/2020-15/01/2020).

**Informed Consent Statement:** Retrospective data collection was approved by the Ethics Committee. The requirement for evidence of informed consent was waived because of the retrospective nature of our study.

**Data Availability Statement:** The data presented in this study may be available on reasonable request from the corresponding author. The data are not publicly available due to ethical and privacy restrictions.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Badr, E. Images in Space and Time: Real Big Data in Healthcare. *ACM Comput. Surv.* **2021**, *54*, 113. [CrossRef]

2. Duncan, J.S.; Ayache, N. Medical image analysis: Progress over two decades and the challenges ahead. *IEEE Trans. Pattern Anal. Mach. Intell.* **2000**, *22*, 85–106. [CrossRef]

3. Ghavami, N.; Hu, Y.; Gibson, E.; Bonmati, E.; Emberton, M.; Moore, C.M.; Barratt, D.C. Automatic segmentation of prostate MRI using convolutional neural networks: Investigating the impact of network architecture on the accuracy of volume measurement and MRI-ultrasound registration. *Med. Image Anal.* **2019**, *58*, 101558. [CrossRef]

4. Lee, B.; Yamanakkanavar, N.; Choi, J.Y. Automatic segmentation of brain MRI using a novel patch-wise U-net deep architecture. *PLoS ONE* **2020**, e0236493. [CrossRef]

5. Pasero, E.; Castagneri, C. Application of an automatic ulcer segmentation algorithm. In *Proceedings of the International Forum on Research and Technologies for Society and Industry (RTSI)*, Modena, Italy, 11–13 September 2017; pp. 1–4. [CrossRef]

6. Yankeeov, T.E.; Mankoff, D.A.; Schwartz, L.H.; Lieberman, F.S.; Buatti, J.M.; Mountz, J.M.; Erickson, B.J.; Fennessy, F.M.; Huang, W.; Kalpathy-Cramer, J.; et al. Quantitative imaging in cancer clinical trials. *Clin. Cancer Res.* **2016**, *22*, 284–290. [CrossRef] [PubMed]

7. Aerts, H.J.; Velazquez, E.R.; Leijenaar, R.T.; Parmar, C.; Grossmann, P.; Carvalho, S.; Bussink, J.; Monshouwer, R.; Haibe-Kains, B.; Rietveld, D.; et al. Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach. *Nat. Commun.* **2014**, *5*, 4006. [CrossRef] [PubMed]

8. Lambin, P.; Leijenaar, R.T.; Deist, T.M.; Peerlings, J.; de Jong, E.E.; van Timmeren, J.; Sanduleanu, S.; Larue, R.T.; Even, A.J.; Jochems, A.; et al. Radiomics: The bridge between medical imaging and personalized medicine. *Nat. Rev. Clin. Oncol.* **2017**, *14*, 749. [CrossRef]

9. Krupinski, E.A. Current perspectives in medical image perception. *Atten. Percept. Psychophys.* **2010**, *72*, 1205–1217. [CrossRef] [PubMed]

10. Rundo, L.; Pirrone, R.; Vitabile, S.; Sala, E.; Gambino, O. Recent advances of HCI in decision-making tasks for optimized clinical workflows and precision medicine. *J. Biomed. Inf.* **2020**, *108*, 103479. [CrossRef]

11. Lee, N.Y.; Lu, J.J. *Target Volume Delineation and Field Setup: A Practical Guide for Conformal and Intensity-Modulated Radiation Therapy*, 1st ed.; Springer: Berlin/Heidelberg, Germany, 2013. [CrossRef]
12. Hamamci, A.; Kucuk, N.; Karaman, K.; Engin, K.; Unal, G. Tumor-Cut: Segmentation of brain tumors on contrast enhanced MR images for radiosurgery applications. *IEEE Trans. Med. Imaging* **2012**, *31*, 790–804. [CrossRef]

13. Rundo, L.; Militello, C.; Russo, G.; Vitabile, S.; Gilardi, M.C.; Mauri, G. GTVcut for neuro-radiosurgery treatment planning: An MRI brain cancer seeded image segmentation method based on a cellular automata model. *Nat. Comput.* **2018**, *17*, 521–536. [CrossRef]

14. Litjens, G.; Kooi, T.; Bejnordi, B.E.; Setio, A.A.A.; Ciompi, F.; Ghafoorian, M.; Van Der Laak, J.A.; Van Ginneken, B.; Sánchez, C.I. A survey on deep learning in medical image analysis. *Med. Image Anal.* **2017**, *42*, 60–88. [CrossRef]

15. Tajbakhsh, N.; Jeyaseelan, L.; Li, Q.; Chiang, J.N.; Wu, Z.; Ding, X. Embracing imperfect datasets: A review of deep learning solutions for medical image segmentation. *Med. Image Anal.* **2020**, *63*, 101693. [CrossRef]

16. Tajbakhsh, N.; Shin, J.Y.; Gurudu, S.R.; Hurst, R.T.; Kendall, C.B.; Gotway, M.B.; Liang, J. Convolutional neural networks for medical image analysis: Full training or fine tuning? *IEEE Trans. Med. Imaging* **2016**, *35*, 1299–1312. [CrossRef] [PubMed]

17. Ravi, D.; Wong, C.; Deligianni, F.; Berthelot, M.; Andreu-Perez, J.; Lo, B.; Yang, G.Z. Deep learning for health informatics. *IEEE J. Biomed. Health Inform.* **2017**, *21*, 4–21. [CrossRef] [PubMed]

18. Militello, C.; Rundo, L.; Dimarco, M.; Orlando, A.; Conti, V.; Woitek, R.; D'Angelo, I.; Bartolotta, T.V.; Russo, G. Semi-automated and interactive segmentation of contrast-enhancing masses on breast DCE-MRI using spatial fuzzy clustering. *Biomed. Signal Process. Control.* **2022**, *71*, 103113. [CrossRef]

19. Acharya, R.; Wasserman, R.; Stevens, J.; Hinojosa, C. Biomedical imaging modalities: A tutorial. *Comput. Med. Imaging Graph.* **1995**, *19*, 3–25. [CrossRef]

20. Tirpude, N.N.; Welekar, R.R. Effect Of Global Thresholding On Tumor-Bearing Brain MRI Images. *Int. J. Eng. Comput. Sci.* **2013**, *2*, 728–731.

21. Militello, C.; Vitabile, S.; Rundo, L.; Russo, G.; Midiri, M.; Gilardi, M.C. A fully automatic 2D segmentation method for uterine fibroid in MRgFUS treatment evaluation. *Comput. Biol. Med.* **2015**, *62*, 277–292. [CrossRef] [PubMed]

22. Islam, M.R.; Imteaz, M.R.; Marium-E-Jannat. Detection and analysis of brain tumor from MRI by Integrated Thresholding and Morphological Process with Histogram based method. In Proceedings of the International Conference on Computer, Communication, Chemical, Material and Electronic Engineering (IC4ME2), Rajshahi, Bangladesh, 8–9 February 2018; pp. 1–5.

23. Rundo, L.; Militello, C.; Vitabile, S.; Casarino, C.; Russo, G.; Midiri, M.; Gilardi, M.C. Combining Split-and-Merge and Multi-Seed Region Growing Algorithms for Uterine Fibroid Segmentation in MRgFUS Treatments. *Med. Biol. Eng. Comput.* **2016**, *54*, 1071–1084. [CrossRef]

24. Horowitz, S.L.; Pavlidis, T. Picture Segmentation by a Tree Traversal Algorithm. *J. ACM* **1976**, *23*, 368–388. [CrossRef]

25. Manousakas, I.N.; Undrill, P.E.; Cameron, G.G.; Redpath, T.W. Split-and-Merge Segmentation of Magnetic Resonance Medical Images: Performance Evaluation and Extension to Three Dimensions. *Comput. Biomed. Res.* **1998**, *31*, 393–412. [CrossRef]

26. Saad, N.M.; Abu-Bakar, S.A.R.; Muda, S.; Mokji, M. Automated segmentation of brain lesion based on diffusion-weighted MRI using a split and merge approach. In Proceedings of the IEEE EMBS Conference on Biomedical Engineering and Sciences (IECBES), Kuala Lumpur, Malaysia, 30 November–2 December 2010; pp. 475–480.

27. Saad, N.M.; Abu-Bakar, S.A.R.; Muda, S.; Mokji, M.; Abdullah, A.R. Automated region growing for segmentation of brain lesion in diffusion-weighted MRI. In Proceedings of the International MultiConference of Engineers and Computer Scientists, IMECS 2012, Hong Kong, China, 14–16 March 2012; pp. 674–677.

28. Adams, R.; Bischof, L. Seeded region growing. *IEEE Trans. Pattern Anal.* **1994**, *16*, 641–647. [CrossRef]

29. Chang, Y.L.; Li, X. Adaptive image region-growing. *IEEE Trans. Image Process.* **1994**, *3*, 868–872. [CrossRef]

30. Otsu, N. A threshold selection method from gray-level histograms. *IEEE Trans. Syst. Man Cybern.* **1975**, *11*, 23–27. [CrossRef]

31. Joseph, R.P.; Senthil Singh, C.; Manikandan, M. Brain tumor MRI image segmentation and detection in image processing. *Int. J. Res. Eng. Technol.* **2014**, *3*, 1–5.

32. Kanungo, T.; Mount, D.M.; Netanyahu, N.S.; Piatko, C.D.; Silverman, R.; Wu, A.Y. An efficient k-means clustering algorithm: analysis and implementation. *IEEE Trans. Pattern Anal. Mach. Intell.* **2002**, *24*, 881–892. [CrossRef]

33. Bezdek, J.C. Objective function clustering. In *Pattern Recognition with Fuzzy Objective Function Algorithms*, 1st ed.; Springer: Secaucus, NJ, USA, 1981; pp. 43–93. [CrossRef]

34. Li, Y.I.; Shen, Y. An automatic fuzzy c-means algorithm for image segmentation. *Soft Comput.* **2010**, *14*, 123–128. [CrossRef]

35. Militello, C.; Rundo, L.; Vitabile, S.; Russo, G.; Pisciotta, P.; Marletta, F.; Ippolito, M.; D'Arrigo, C.; Midiri, M.; Gilardi, M. Gamma Knife treatment planning: MR brain tumor segmentation and volume measurement based on unsupervised Fuzzy C-Means clustering. *Int. J. Imaging Syst. Technol.* **2015**, *25*, 213–225. [CrossRef]

36. Rundo, L.; Militello, C.; Russo, G.; Garufi, A.; Vitabile, S.; Gilardi, M.C.; Mauri, G. Automated prostate gland segmentation based on an unsupervised fuzzy c-means clustering technique using multispectral T1w and T2w MR imaging. *Information* **2017**, *8*, 49. [CrossRef]

37. Caponetti, L.; Castellano, G.; Corsini, V. MR brain image segmentation: A framework to compare different clustering techniques. *Information* **2017**, *8*, 138. [CrossRef]

38. Rundo, L.; Militello, C.; Tangherloni, A.; Russo, G.; Vitabile, S.; Gilardi, M.C.; Mauri, G. NeXt for neuro-radiosurgery: A fully automatic approach for necrosis extraction in brain tumor MRI using an unsupervised machine learning technique. *Int. J. Imaging Syst. Technol.* **2018**, *28*, 21–37. [CrossRef]

39. Feder, J.M.; de Paredes, E.S.; Hogge, J.P.; Wilken, J.J. Unusual breast lesions: Radiologic-pathologic correlation. *Radiographics* **2019**, *19*, S11–S26. [CrossRef] [PubMed]
40. Chuang, K.S.; Tzeng, H.L.; Chen, S.; Wu, J.; Chen, T.J. Fuzzy c-means clustering with spatial information for image segmentation. *Comput. Med. Imaging Graph.* **2006**, *30*, 9–15. [CrossRef] [PubMed]

41. Rundo, L.; Beer, L.; Ursprung, S.; Martin-Gonzalez, P.; Markowetz, F.; Brenton, J.D.; Crispin-Ortuzar, M.; Sala, E.; Woitek, R. Tissue-specific and interpretable sub-segmentation of whole tumour burden on CT images by unsupervised fuzzy clustering. *Comput. Biol. Med.* **2020**, *120*, 103751. [CrossRef] [PubMed]

42. Patil, S. Preprocessing To Be Considered For MR and CT Images Containing Tumors. *IOSR J. Electr. Electron. Eng.* **2012**, *1*, 55–57. [CrossRef]

43. Vasuki, P.; Kanimozhi, J.; Devi, M.B. A survey on image preprocessing techniques for diverse fields of medical imagery. In Proceedings of the 2017 IEEE International Conference on Electrical, Instrumentation and Communication Engineering (ICEICE), Tamilnadu, India, 27–28 April 2017; pp. 1–6.

44. Behrenbruch, C.; Petroudi, S.; Bond, S.; Declerck, J.; Leong, F.; Brady, J. Image filtering techniques for medical image post-processing: An overview. *Br. J. Radiol.* **2004**, *77*, S126–S132. [CrossRef] [PubMed]

45. Seeram, E.; Seeram, D. Image Postprocessing in Digital Radiology—A Primer for Technologists. *J. Med. Imaging Radiat. Sci.* **2008**, *39*, 23–41. [CrossRef]

46. Castiglioni, I.; Rundo, L.; Codari, M.; Di Leo, G.; Salvatore, C.; Interlenghi, M.; et al. AI applications to medical images: From machine learning to deep learning. *Phys. Med.* **2021**, *83*, 9–24. [CrossRef]

47. Motwani, M.C.; Gadiya, M.C.; Motwani, R.C.; Harris, F.C. Survey of image denoising techniques. *Glob. Signal Process. Expo Conf. (GSPX)* **2004**, *27*, 27–30.

48. Marrone, S.; Piantadosi, G.; Fusco, R.; Petrillo, A.; Sansone, M.; Sansone, C. Breast segmentation using Fuzzy C-Means and anatomical priors in DCE-MRI. In Proceedings of the 23rd International Conference on Pattern Recognition (ICPR), Cancun, Mexico, 4–8 December 2016; pp. 1472–1477.

49. Piantadosi, G.; Fusco, R.; Petrillo, A.; Sansone, M.; Sansone, C. LBP-TOP for volume lesion classification in breast DCE-MRI. In *International Conference on Image Analysis and Processing; Lecture Notes in Computer Science; Springer: Cham, Switzerland, 2015; pp. 647–657*. [CrossRef]

50. Wilcoxon, F. Individual comparisons by ranking methods. *Biom. Bull.* **1980**, *1*, 196–202. [CrossRef]

51. Holm, S. A Simple Sequentially Rejective Multiple Test Procedure. *Scand. J. Stat.* **1979**, *6*, 65–70.

52. DeSantis, C.E.; Ma, J.; Gaudet, M.M.; Newman, L.A.; Miller, K.D.; Goding Sauer, A.; Jemal, A.; Siegel, R.L. Breast cancer statistics, 2019. *CA Cancer J. Clin.* **2019**, *69*, 438–451. [CrossRef]

53. Siegel, R.L.; Miller, K.D.; Fuchs, H.E.; Jemal, A. Cancer Statistics, 2021. *CA Cancer J. Clin.* **2021**, *71*, 7–33. [CrossRef] [PubMed]

54. International Agency for Research on Cancer. The Global Cancer Observatory. 2020. Available online: https://gco.iarc.fr/today/data/factsheets/cancers/20-Breast-fact-sheet.pdf (accessed on 28 July 2021).

55. Frackiewicz, M.; Koper, Z.; Palus, H.; Borys, D.; Psiuk-Maksymowicz, K. Breast lesion segmentation in DCE-MRI Imaging. In Proceedings of the 14th International Conference on Signal-Image Technology & Internet-Based Systems (SITIS), Las Palmas de Gran Canaria, Spain, 26–29 November 2018; pp. 308–313. [CrossRef]

56. Granzier, R.; Verbakel, N.; Ibrahim, A.; van Timmeren, J.; van Nijnatten, T.; Leijenaar, R.; Lobbes, M.; Smidt, M.; Woodruff, H. MRI-based radiomics in breast cancer: Feature robustness with respect to inter-observer segmentation variability. *Sci. Rep.* **2020**, *10*, 14163. [CrossRef] [PubMed]

57. Han, C.; Rundo, L.; Murao, K.; Noguchi, T.; Shimahara, Y.; Milacski, Z.Á.; Koshino, S.; Sala, E.; Nakayama, H.; Satoh, S. MADGAN: Unsupervised medical anomaly detection GAN using multiple adjacent brain MRI slice reconstruction. *BMC Bioinf.* **2021**, *22*, 31. [CrossRef]

58. Wu, X.; Bi, L.; Fulham, M.; Kim, J. Unsupervised Positron Emission Tomography Tumor Segmentation via GAN based Adversarial Auto-Encoder. In Proceedings of the 16th International Conference on Control, Automation, Robotics and Vision (ICARCV), Shenzhen, China, 13–15 December 2020; pp. 448–453. [CrossRef]

59. Rundo, L.; Stefano, A.; Militello, C.; Russo, G.; Sabini, M.G.; D'Arrigo, C.; Marletta, F.; Ippolito, M.; Mauri, G.; Vitabile, S.; et al. A fully automatic approach for multimodal PET and MR image segmentation in Gamma Knife treatment planning. *Comput. Methods Programs Biomed.* **2017**, *144*, 77–96. [CrossRef]

60. Woitek, R.; McLean, M.A.; Gill, A.B.; Grist, J.T.; Provenzano, E.; Patterson, A.J.; Ursprung, S.; Torheim, T.; Zaccagna, F.; Locke, M.; et al. Hyperpolarized 13C MRI of tumor metabolism demonstrates early metabolic response to neoadjuvant chemotherapy in breast cancer. *Radiol. Imaging Cancer* **2020**, *2*, e200017. [CrossRef]

61. Carmona-Bozo, J.C.; Manavaki, R.; Woitek, R.; Torheim, T.; Baxter, G.C.; Caracò, C.; Provenzano, E.; Graves, M.J.; Fryer, T.D.; Patterson, A.J.; et al. Hypoxia and perfusion in breast cancer: Simultaneous assessment using PET/MR imaging. *Eur. Radiol.* **2021**, *31*, 333–344. [CrossRef]
Article

# Development of Detection and Volumetric Methods for the Triceps of the Lower Leg Using Magnetic Resonance Images with Deep Learning

Yusuke Asami<fnref n="1" />, Takaaki Yoshimura<fnref n="2" />,<fnref n="3" />, Keisuke Manabe<fnref n="1" />, Tomonari Yamada<fnref n="1" /> and Hiroyuki Sugimori<fnref n="4" />,<fnref n="5" />,*

**Abstract:** Purpose: A deep learning technique was used to analyze the triceps surae muscle. The devised interpolation method was used to determine muscle's volume and verify the usefulness of the method. Materials and Methods: Thirty-eight T1-weighted cross-sectional magnetic resonance images of the triceps of the lower leg were divided into three classes, i.e., gastrocnemius lateralis (GL), gastrocnemius medialis (GM), and soleus (SOL), and the regions of interest (ROIs) were manually defined. The supervised images were classified as per each patient. A total of 1199 images were prepared. Six different datasets separated patient-wise were prepared for K-fold cross-validation. A network model of the DeepLabv3+ was used for training. The images generated by the created model were divided as per each patient and classified into each muscle types. The model performance and the interpolation method were evaluated by calculating the Dice similarity coefficient (DSC) and error rates of the volume of the predicted and interpolated images, respectively. Results: The mean DSCs for the predicted images were >0.81 for GM and SOL and 0.71 for GL. The mean error rates for volume were approximately 11% for GL, SOL, and total error and 23% for GL. DSCs in the interpolated images were >0.8 for all muscles. The mean error rates of volume were <10% for GL, SOL, and total error and 18% for GM. There was no significant difference between the volumes obtained from the supervised images and interpolated images. Conclusions: Using the semantic segmentation of the deep learning technique, the triceps were detected with high accuracy and the interpolation method used in this study to find the volume was useful.

**Keywords:** deep learning; semantic segmentation; triceps surae muscle

**Citation:** Asami, Y.; Yoshimura, T.; Manabe, K.; Yamada, T.; Sugimori, H. Development of Detection and Volumetric Methods for the Triceps of the Lower Leg Using Magnetic Resonance Images with Deep Learning. *Appl. Sci.* **2021**, *11*, 12006.
https://doi.org/10.3390/app112412006

Academic Editors: Leonardo Rundo, Carmelo Militello and Andrea Tangherloni

Received: 9 November 2021

Accepted: 14 December 2021

Published: 16 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_52.png" /></figure>

Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Deep learning technology has been widely used in recent years for automatic driving, drones, weather forecasting, and games [1–4]. The deep learning techniques include classification [5], object detection [6], and semantic segmentation. Semantic segmentation has been used to visualize the three-dimensional (3D) anatomical structures of multiple organs in 3D computed tomography images [7], to facilitate the quantitative coronary angiography-based diagnosis of major vessels in X-ray coronary angiography [8], and to quantitate whole breast image analysis in diffusion-weighted images [9]. An applied method using Conditional Generative Adversarial Networks (cGANs) has been reported [10] to be highly accurate for in segmenting peri-knee tissue using magnetic resonance imaging (MRI).

<footnote n="1">¹ Graduate School of Health Sciences, Hokkaido University, Sapporo 060-0812, Japan; yusuke12@eis.hokudai.ac.jp (Y.A.); ksk0843@eis.hokudai.ac.jp (K.M.); tomonarihandball@eis.hokudai.ac.jp (T.Y.)</footnote>

<footnote n="2">² Department of Health Sciences and Technology, Faculty of Health Sciences, Hokkaido University, Sapporo 060-0812, Japan; takaaki.ysm@med.hokudai.ac.jp</footnote>

<footnote n="3">³ Department of Medical Physics, Hokkaido University Hospital, Sapporo 060-8648, Japan</footnote>

<footnote n="4">⁴ Department of Biomedical Science and Engineering, Faculty of Health Sciences, Hokkaido University, Sapporo 060-0812, Japan</footnote>

<footnote n="5">5 Clinical AI Human Resources Development Program, Faculty of Medicine, Hokkaido University, Sapporo 060-8648, Japan</footnote>

<footnote>* Correspondence: sugimori@hs.hokudai.ac.jp; Tel.: +81-11-706-3410</footnote>
These techniques are used in various body parts and modalities. This study focused on semantic segmentation of the triceps surae muscle, which consists of the gastrocnemius and soleus (SOL) muscles. The gastrocnemius is a biceps muscle that comprises the gastrocnemius lateralis (GL) and gastrocnemius medialis (GM). The triceps is an important muscle because it makes activities, such as walking and standing, possible by its repeated contraction and relaxation. Triceps muscle weakness increases the risk of falling and prevents these essential activities [11]. The muscle volume is commonly used as a measure of the muscle size, and it is frequently measured in studies because it correlates with various functional parameters [12]. In a study by Thom et al., the volume of the triceps femoris muscle was calculated to assess the loss of muscle mass with aging [13]. The volume was calculated by the product of the cross-sectional area and slice thickness obtained by manual segmentation using an imaging software; however, this procedure takes a long time [14]. Friedberger et al. [15] showed that it is now possible to perform semi-automatic volume calculations using the Random Forest classifier for hand muscle segmentation. However, this method is still semi-automatic, and even though it reduces the effort of manual segmentation, it still requires modifications and is significantly time-consuming. In addition, since magnetic resonance imaging of the lower extremities is often performed with thicker slices and more spacing between images, volume calculations using only the acquired images are likely to differ significantly from the actual volume. Therefore, it was believed that by slice interpolating the MRI images and calculating the volume using the interpolated images, this difference from the actual image could be reduced.

Many studies have been conducted on the segmentation of quadriceps [16–21], there have been no previous studies focusing on the segmentation of triceps with deep learning. In a study conducted by Essafi et al. [22], the authors performed the segmentation of the medial gastrocnemius, one of the triceps muscles, without deep learning and reported that the average Dice similarity coefficient (DSC) in 25 subjects was 0.55. In the studies of thigh muscle segmentation using deep learning [16–19], the higher the DSC > 0.9, the higher the accuracy of muscle detection. Furthermore, in the studies of segmentation of thigh muscles without deep learning [20,21], there are examples of detection with relatively high accuracy, although not as high compared with studies that have used deep learning. Based on these results, that the triceps could be detected with high accuracy using the deep learning method.

In this study, deep learning was used to detect the triceps surae muscle. Then, the devised interpolation method was used to determine the muscle volume and verify the usefulness of the method.

The manuscript is structured as follows. Section 2 outlines the preprocessing of the dataset and the methodology and evaluations for transfer learning and the image interpolation method for semantic segmentation. Section 3 describes the created models and the evaluation of the image interpolation method. In Section 4, the results obtained in this study are compared with those of other papers, and the limitations of this study are discussed. Finally, the conclusions and future directions of this study are presented in Section 5.

## 2. Materials and Methods

### 2.1. Subjects

The 38 consecutive patients who underwent a noncontrast-enhanced lower leg MRI examination were retrospectively evaluated. The Digital Imaging and Communications in Medicine (DICOM) images were subsequently converted to anonymized DICOM files from the image server. These images have been approved by the Ethics Committee of the Hokkaido University Hospital. The details information of the obtained images is shown in Table 1.
<table><caption>**Table 1.** The details information of the obtained images (N = 38).</caption><thead><tr><th></th><th>Mean ± SD [Min–Max]</th></tr></thead><tbody><tr><td>Number of slices</td><td>17.9 ± 4.9 [3–29]</td></tr><tr><td>Field of view [mm]</td><td>342.1 × 295.0 [160 × 160–500 × 425]</td></tr><tr><td>Acquisition matrix size [pixel]</td><td>446.0 × 402.4 [320 × 224–672 × 672]</td></tr><tr><td>Pixel size [mm]</td><td>0.764 × 0.764 [0.559 × 0.559–0.928 × 0.928]</td></tr><tr><td>Slice thickness [mm]</td><td>5.8 ± 0.5 [4–6]</td></tr><tr><td>Slice gap [mm]</td><td>13.6 ± 4.7 [4.8–22.8]</td></tr><tr><td>Length of acquisition [mm]</td><td>231.0 ± 76.2 [42–342]</td></tr></tbody></table>

## 2.2. Preprocessing

For supervised image creation, the DICOM images were converted to 8-bit using the MATLAB (The MathWorks, Inc., Natick, MA, USA) "mat2gray" function because the bit depth of DICOM images is 16-bit. Those with only one leg captured were resized to 256 × 256 portable network graphics (PNG) images, whereas those with both legs captured were trimmed to the appropriate size for each leg and resized to 256 × 256. The left leg was reversed and oriented in the same way as the right leg (Figure 1). The coordinates of the trimmed area were also recorded to estimate the pixel size and the number of pixels in a trimmed image for volume calculation. Images of 67 leg muscles were obtained from 38 patients. The trimmed images were used for training and evaluation.

<figure><img src="image_53.png" /><figcaption>**Figure 1.** The conversion process from DICOM to PNG image.</figcaption></figure>

## 2.3. Dataset

These images were saved as PNG images by dividing the triceps into three classes (GM, GL, and SOL) with an in-house MATLAB software; the regions of interest (ROIs) were also manually defined (Figure 2). The ROIs were manually defined by the author alone and were then verified by two radiological technologists with 5 (T.Y.) and 20 years of experience (H.S.), respectively. A total of 1199 images were prepared from 38 patients, and the Train:Test was set to 10:2. The 38 patients were classified into six groups, and six datasets were created so that one could be used as a test (Figure 2). Six different datasets with 1199 images were prepared for K-fold cross-validation so that the number of
tests would be approximately 200. Data augmentation for improving the training [6] was
performed on the training images, and they were rotated by 5° from -45 to 45°, increasing
the number of images by 19 times.

<figure><img src="image_54.png" /><figcaption>Figure 2. Six divided datasets and 6-fold cross-validation.</figcaption></figure>

## 2.4. *Training for Creating Models*

The software was developed with MATLAB software, and a computer with NVIDIA GeForce GTX 1080Ti 12GB (NVIDIA Corporation, Santa Clara, CA, USA) was used. DeepLabv3+ was used for the architecture for implementation of image segmentation models. The model was inputted with 2D images taken by MRI. The following training parameters were used: the batch size for the number of training samples was 32, the number of epochs was 100, and the initial learning rate was 0.0001. Using this network model, six different datasets were trained to create six models.

## 2.5. Interpolation

The supervised images and the images generated using the devised model were
divided in each patient and classified into three muscle types (GM, GL, and SOL). Only the
images with muscles present in the supervised images were extracted and interpolated
by excluding even-numbered images from the predicted images in the same position
(Figures 3 and 4). The interpolation method is described in Figure 5. Interpolation was
performed in the slice direction from two images of the muscle. The contours of the two
muscles were represented by points. The centroid of the larger muscle was calculated
and connected to the centroid by a point on the larger muscle. The distance between the
connecting line and all points on the smaller muscle was calculated, and the point with
the smallest distance was adopted. Moreover, the midpoint was calculated by connecting the adopted point to a point on the larger muscle. The same process was performed for all points of the larger muscle, and by connecting all the calculated midpoints, the interpolated image of the muscle was completed.

<figure><img src="image_55.png" /><figcaption>Figure 3. Classifying muscles into three for interpolation.</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>**Figure 4.** Creating the interpolated images from the predicted images based on the presence or absence of muscles in the supervised images. (A) Muscles in supervised images; (B) Predictive muscles using segmentation model; (C) Interpolated muscles using the interpolation method by thinning out the even number of predicted images; (D) Muscles with even-numbered predicted images thinned out and slice thickness doubled without using the interpolation method; 1–9; The number of the slice where the muscle is located.</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>Figure 5. Procedure of interpolation: (1) The outline with dots and the calculated centroid of the muscle with the most dots were expressed; (2) a line was drawn connecting the red dot and the centroid; (3) measured the distance between all the blue points and the yellow line, found the nearest blue point from yellow line; (4) the midpoint of the red dot and the selected blue dot was marked; (5) a different red dot was selected and steps 2 through 4 were repeated for as many red dots as there are red dots, connecting all the yellow dots.</figcaption></figure>

## 2.6. Indicators Used for Evaluation

### 2.6.1. DSC

When the supervised images were set as A and the predicted images were set as B, the following formula was used to calculate DSC. This value is an index to check the agreement between the images, and the closer the value is to 1, the better the agreement is.

$$DSC = \frac{2 \times |A \cap B|}{|A| + |B|} \qquad (1)$$

### 2.6.2. Calculation of the Volume and Error Rate

The spacing between the slices and pixels were obtained from the DICOM data. The number of pixels used in the labeled muscle was determined, and the volume was calculated from the product of these values. The error rate was calculated to compare the volume of the supervised images with that of the predicted images. The error rate was calculated using the following formula.

$$\left| \frac{\text{Volume of predicted images} - \text{Volume of supervised images}}{\text{Volume of supervised images}} \right| \times 100 (\%) \qquad (2)$$

## 2.7. Evaluation of the Created Models

DSCs were calculated for each muscle using the abovementioned formula, and the averages were calculated for each dataset to show the overall average. All images were used for this evaluation. The volumes were calculated for each patient and muscle, and the average error rates of the datasets were calculated using the mentioned formula. For this evaluation, only the images with muscle present in both the supervised and predicted images were used to calculate the error rate.
## 2.8. Evaluation of the Interpolation Method

Only the interpolated images were used to obtain DSCs, and the average DSC for each muscle was calculated. DSCs were calculated from the interpolated image and the supervised image of the same number (Figure 4C). The volumes were calculated from the interpolated images, and the error rates were calculated for comparison with the volumes of the supervised images. To confirm that the interpolation method is useful in terms of reducing volumetric errors, the error rate was calculated by comparing the volume of the supervised image and the volume of the image with the spacing between its slices doubled without interpolation while excluding even numbered images (Figure 4D). All patients with all three muscles and all predicted images that could be interpolated were used for these evaluations. Statistical tests were performed to compare the interpolated and supervised volumes. First, the Shapiro-Wilk test was applied to evaluate normality; if normality was accepted, the *t*-test was used; otherwise, the Wilcoxon signed-ranked test, a nonparametric test, was used. The differences in the statistical analyses were considered statistically significant when *p* was <0.05.

# 3. Results

## 3.1. Evaluation of the Created Models

The DSC values per dataset are shown in Table 2. The following table shows the DSC values of all 67 leg muscles. GM and SOL were detected with relatively high accuracy. However, the DSCs of GL were lower than those of the others. The following are some representative examples of segmentation successes and failures (Figures 6 and 7)

<table><caption>**Table 2.** Average DSC values per dataset for the predicted images.</caption><thead><tr><th>Position</th><th>Dataset1</th><th>Dataset2</th><th>Dataset3</th><th>Dataset4</th><th>Dataset5</th><th>Dataset6</th><th>Mean ± SD</th></tr></thead><tbody><tr><td>GM</td><td>0.897</td><td>0.754</td><td>0.876</td><td>0.830</td><td>0.776</td><td>0.745</td><td>0.813 ± 0.064</td></tr><tr><td>GL</td><td>0.804</td><td>0.591</td><td>0.765</td><td>0.773</td><td>0.639</td><td>0.704</td><td>0.713 ± 0.084</td></tr><tr><td>SOL</td><td>0.889</td><td>0.820</td><td>0.791</td><td>0.823</td><td>0.840</td><td>0.844</td><td>0.835 ± 0.033</td></tr></tbody></table>

GM: gastrocnemius medialis, GL: gastrocnemius lateralis, SOL: soleus, SD: standard deviation.

<figure><img src="image_58.png" /><figcaption>**Figure 6.** Representative examples of accurately predicted images.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>**Figure 7.** Representative examples of partially incorrectly predicted images.</figcaption></figure>

The error rates of volume per dataset are shown in Table 3. The following table shows the error rates of the volumes of all 67 legs that were present in the supervised and predicted images. The total error represents the error rate of all three muscle volumes. The error rate of MG, SOL, and sum was 10% but that of the GL volume was large.

<table><caption>**Table 3.** The error rate of volume per dataset for the predicted images.</caption><thead><tr><th>Position</th><th>Dataset1</th><th>Dataset2</th><th>Dataset3</th><th>Dataset4</th><th>Dataset5</th><th>Dataset6</th><th>Mean ± SD</th></tr></thead><tbody><tr><td>GM</td><td>5.68</td><td>11.84</td><td>8.96</td><td>10.28</td><td>8.80</td><td>24.48</td><td>11.67 ± 6.60</td></tr><tr><td>GL</td><td>27.46</td><td>23.03</td><td>18.74</td><td>27.14</td><td>32.71</td><td>12.36</td><td>23.57 ± 7.22</td></tr><tr><td>SOL</td><td>11.71</td><td>17.22</td><td>14.60</td><td>6.87</td><td>5.92</td><td>9.24</td><td>10.93 ± 4.43</td></tr><tr><td>total error</td><td>11.07</td><td>13.38</td><td>11.87</td><td>7.27</td><td>7.84</td><td>12.43</td><td>10.64 ± 2.51</td></tr></tbody></table>

[%], GM: gastrocnemius medialis, GL: gastrocnemius lateralis, SOL: soleus, SD: standard deviation.

### 3.2. Evaluation of the Interpolation Method

The DSC values of the interpolated images are shown in Table 4. The following table shows the DSC of 60 GMs, 54 GLs, and 66 SOLs in the interpolated images. All muscles had mean DSC larger than 0.8.

<table><caption>**Table 4.** The DSC values of the interpolated images.</caption><thead><tr><th>Position</th><th>DSC</th></tr></thead><tbody><tr><td>GM</td><td>0.877 ± 0.134</td></tr><tr><td>GL</td><td>0.809 ± 0.170</td></tr><tr><td>SOL</td><td>0.867 ± 0.078</td></tr></tbody></table>

GM: gastrocnemius medialis, GL: gastrocnemius lateralis, SOL: soleus, DSC: Dice similarity coefficient.

The error rates of volume for the interpolated images are shown in Table 5. The following table shows the error rate between the volume of the supervised images (Figure 4A) and the volume calculated from the 52 legs to which all the three muscles could be interpolated (Figure 4C). For GM, SOL, and the total error, the percentage was <10%. In
fact, the percentage of the total error was as low as 7.5%. The usefulness of this interpolation method in reducing volume calculation errors was evaluated by comparing the error rate between the volume of the supervised image (Figure 4A) and the volume of the non-interpolated image (Figure 4D). In the comparison between the supervised and non-interpolated muscles, all muscles were >12%. Both volumes were calculated for each leg and were then statistically tested. *p* value results showed that there were no significant differences between the volumes of the supervised and interpolated muscles. The same sequence for the supervised and non-interpolated muscle were performed, and the *p* value results also showed that there were no significant differences.

<table><caption>**Table 5.** The error rate of volume for the interpolated images.</caption><thead><tr><th>Comparison Target</th><th>Position</th><th>Error Rate [%]</th><th>p Value</th></tr></thead><tbody><tr><td rowspan="4">Supervised muscle and<br>interpolated muscle</td><td>GM</td><td>9.41 ± 7.65</td><td>0.5052</td></tr><tr><td>GL</td><td>17.89 ± 6.37</td><td>0.3098</td></tr><tr><td>SOL</td><td>9.43 ± 5.07</td><td>0.7280</td></tr><tr><td>total error</td><td>7.69 ± 3.49</td><td>-</td></tr><tr><td rowspan="4">Supervised muscle and<br>non-interpolated muscle</td><td>GM</td><td>20.28 ± 9.68</td><td>0.1499</td></tr><tr><td>GL</td><td>33.93 ± 7.54</td><td>0.2539</td></tr><tr><td>SOL</td><td>12.70 ± 5.79</td><td>0.3442</td></tr><tr><td>total error</td><td>14.97 ± 3.98</td><td>-</td></tr></tbody></table>

GM: gastrocnemius medialis; GL: gastrocnemius lateralis; SOL: soleus. Supervised muscle: Figure 4A; interpolated muscle: Figure 4C; non-interpolated muscle: Figure 4D.

## 4. Discussion

With regards to the evaluation of the created models, DSCs were used to compare the agreement of the images and evaluate whether the detection was correct in the two dimensions. Based on the results of the evaluation of the created models, MG and SOL were detected with higher accuracy than GL, with DSCs exceeding 0.8. As for MG and SOL, the DSC of segmentation in the thigh muscle, which exceeded 0.8 in many previous studies, were found to be lower than the respective values in DSCs, but the values themselves cannot be compared because the target areas were different. However, this study could detect MGs with a higher accuracy than the reported segmentation of MGs without using deep learning [22]. For GM and SOL, there were many large muscles present, but there were some false positives, which were detected with a DSC of >0.8 on average. As for the DSCs of GL, GL was shorter in the craniocaudal direction than the other two muscles, and the number of slices in which GL exists is small. The slices at the edge of the GLs were so small that it was difficult for the human eye to recognize them, and the ratio of the number of small GLs slices was higher than that of other muscles. In small muscles where the other muscles and surrounding fat tissues are mistakenly detected as GL, false positives were more common and this was believed to be the cause of the decline in DSCs. Figure 7 is a representative example of partially incorrectly predicted images. Not only are there false detections of the SOL as the GL and GM, but there are also examples of inaccurate location. As per Ding's study [17], when false positives, which were not present in the surveillance images, appeared in the prediction images, DSC was zero and affected the average DSC. In addition, some images used in the test were blurry, and false positives were particularly common in these. The triceps surae is a long muscle that exists in a craniocaudal direction, and the shape of the muscle changes, particularly at the level of the knee and calf. In this study, the number of slices was uneven in each position; therefore, the number of slices at the knee level was small, which affected DSCs. As there were no reports on the segmentation of the triceps surae, the results were compared with those of the quadriceps. In a study by Kemnitz et al. [16], using T1-weighted 250 images of thighs, including the quadriceps, hamstrings, sutures, adductor muscles, and others, were prepared and trained with U-Net, which reported a muscle-only DSC of >0.9. However, as this was only a result for a specific anatomical location, it may not be suitable for the volume calculation and interpolation that were performed in this study. In a report by Ding et al. [17], 46 thigh
4968-slice Fat-water decomposition MRI images were trained with U-Net. The mean DSC of the quadriceps, hamstrings, sutures, and thin muscles was 0.89, and the mean percent difference in volume was 7.57%. A similar study was conducted to determine DSC and error rates for different image types and regions used; the results showed a similar trend to the present study, indicating the validity of the present results. With regards to volumetry without interpolation, the agreement was represented in three dimensions by calculating the error rates of the volume.

Based on the average error rates of volumes per dataset, the error rate of GM, SOL, and sum was 10%. The error rate in the GL volume was large, probably due to the presence of many false positives as in DSC. In addition, the predicted images were overestimated in all muscles, which might be due to the false positives because they were surrounded by other muscles and fat. Of the images used for training, the images with all three muscles present were the most common; however, there were few images at the level of the knee where only the gastrocnemius muscle was present, and the quadriceps muscles were mistakenly detected as gastrocnemius. In addition, images with contrast and roughness were also detected that made it difficult to see the boundaries between the muscles, and false positives were frequent in these images.

The interpolation method was evaluated in-plane by calculating the DSC of the supervised and interpolated images. All muscles had mean DSC larger than 0.8. The evaluation of Table 4 was DSC evaluated only with interpolated images; hence, it can be said that the interpolation was performed with high accuracy. In this interpolation method, the centroid of the larger muscle contour was discovered, and the centroid was connected to the dots that make up the contour; the dot with the shortest distance between the line and that of the smaller muscle contour were adopted. If a dot in the false positive area was adopted, the DSC was lowered. Even if they were not considered in this interpolation, all of them were >0.8, so the correct dot with a higher probability among the corresponding dots in the small contour was selected. In a study by Yap Abdullah et al. [23], the results of an interpolation to measure the intracranial volume were comparable to those of manual segmentation, and its use led to the development of a software for rapid measurement. In a study by Nordez et al. [24], to calculate the volume of the quadriceps muscle in MRI images, four different methods were used to calculate the error rate; it was reported that the interpolation method was the best. These results using the interpolation method showed the same trend as these studies.

The interpolation method was evaluated along the slice direction by calculating and comparing the error rates of muscle volumes, including the supervised and interpolated images. The volume contained interpolated images generated from the two images but with low error rates. Therefore, the results in Table 5 demonstrate that the interpolation method is effective in reducing the volumetric error. The percentages of GM, SOL, and total errors were >10%; for the total error, the percentage was as low as 7.5% because if there were no false positives in the images before and after the interpolated images, there were no false positives in the interpolated images. This was believed to lead to a decreased error rate because of the possibility of reduced false positives. Even if false positives existed, it seemed that the probability of adopting dots that constitute false positives was low. There was no significant difference between the volumes of the supervised and interpolated images, indicating that this method can be used without any problems in calculating the volume. From the results in Tables 4 and 5, interpolation was performed by excluding the even numbered images from the subsequent evaluation. As these findings underline that the interpolation was accomplished with high accuracy, while observing at the same time the effect of volume reduction, the error with the actual volume could be further reduced by interpolating the actual captured images that were found.

This study has some limitations. First, this study was unable to interpolate all images used in the test. Because muscles were present in the supervised images and the even number of those present in the predicted images were excluded, if the muscles were small in the supervised images, they could not be detected and the images could not be
complemented. Therefore, the solution is to improve the detection capability by training
more images or to prepare a test image where the number of interpolatable images exists.
Second, the error rate of the volume calculated in this study could not be calculated in
the legs without gastrocnemius or SOL muscles. Therefore, to get a more correct error
rate, it would be necessary to have all the three muscle types in the legs used for training
or testing. Third, most studies on quadriceps segmentation used U-Net [16–18]. In these
studies [25,26], network models were compared for segmentation, and deeplabv3puls
was shown to have the highest detection ability. Therefore, DeepLabv3+ was used for
segmentation instead of U-Net because Deeplabv3+, which has a higher detection capability,
was more likely to be able to detect muscles with higher accuracy than U-net. Using this
algorithm, the triceps could be detected immediately and determine its volume. However,
the recent advances in the channel and spatial attention mechanisms [27–29] have definitely
led to the evolution of semantic segmentation using U-Net. Therefore, a further study on
CNNs in the state of the art is needed. Fourth, the hyperparameters were not adjusted
because this study did not aim to tune the hyperparameters as in the study by Chieh
et al. [30]. Setting and training the optimal hyperparameters may lead to more accurate
detection with Bayesian Optimization. Fifth, the sample size of data is small. Not only
is it difficult to guarantee generalization performance with small sample size, but it is
also possible to improve accuracy by increasing sample size. Additionally, in this study,
training and test images were divided 10:2 for six-fold cross validation, but five- or ten-fold
is general [31,32]. This value was obtained by dividing the data by person, and it is thought
that cross-validation can be performed with general values by increasing the sample size.
Sixth, there are no studies of triceps surae muscle segmentation using Deep Learning.
Therefore, this study was compared with a study of quadriceps segmentation [15–20]
and a study of triceps surae muscle segmentation that did not use deep learning [21].
The results showed that the DSC of the quadriceps was higher than the DSC of this
study, which suggests that this study has room for DSC improvement. In addition, it
was more accurate than the reported segmentation of the triceps surae muscle. This
interpolation method can be used to calculate their volume without any problems found.
Although muscle perimeter length and cross-sectional area, which are correlated with
muscle volume [33,34], are sometimes used as indices, it is more accurate to calculate
volume from images of the entire muscle [35]. However, it is more accurate to calculate the
volume from the whole muscle image. Therefore, this study, which can calculate the actual
volume by immediately calculating the volume and interpolating between slices, suggests
the possibility of advancing research in the field of physical therapy. The technique of
detecting the triceps and calculating their volume can be used to other parts and organs,
and it was assumed that the interpolation method can be used to determine the volume
more accurately. The interpolation technique may also be useful in smoothing out the
images captured with a thicker slice for multiplanar reconstruction processing and 3D
visualization.

5. Conclusions

This deep learning segmentation technique could detect the triceps with relatively
high DSCs. This eliminates the time needed to manually identify the muscles, and thus
enables us to immediately calculate the respective volume. The ability to immediately
calculate the volume of the triceps muscle will provide a distinct contribution to research
studies in the field of physical therapy, including the ability to immediately perform the
quantitative evaluations of rehabilitation devices in leg muscles. The interpolation method
used in this study to determine the volume was useful. Using the interpolation method
allowed us to find that the difference between the volume calculated by the interpolation
method and the actual measurement could be smaller than the volume calculated solely
from the captured image. In the future, this interpolation method can be used to calculate
the volume of other muscles and organs as well, thereby making it possible to calculate the
respective volume in a more accurate manner.
**Author Contributions:** Y.A. contributed to the data analysis, algorithm construction, and the writing and editing of the manuscript. T.Y. (Takaaki Yoshimura), K.M. and T.Y. (Tomonari Yamada) reviewed and edited the manuscript. H.S. proposed the idea and contributed to the data acquisition, performed supervision, project administration, and reviewed and edited the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Hokkaido University Hospital.

**Informed Consent Statement:** Informed consent was waived because of the retrospective nature of the study and the analysis used anonymous clinical data.

**Data Availability Statement:** The created models in this study are available on request from the corresponding author. However, the image datasets presented in this study are not publicly available due to ethical reasons, e.g., because they contain sensitive information that could compromise the privacy of the participants. The source code of this study is available at https://github.com/MIA-laboratory/TricepsMRIsegmentation (accessed on 7 December 2021).

**Acknowledgments:** The authors would like to thank the laboratory members of the Medical Image Analysis laboratory and Yoshimura's laboratory for their help. The authors would also like to thank Masayuki Kugimoto for providing us with useful advice and information.

**Conflicts of Interest:** The authors declare that no conflict of interest exist.

## References

1. Son, S.; Jeong, Y.; Lee, B. An audification and visualization system (AVS) of an autonomous vehicle for blind and deaf people based on deep learning. *Sensors* **2019**, *19*, 5053. [CrossRef]

2. Chen, Y.; Aggarwal, P.; Choi, J.; Jay, C.C. A deep learning approach to drone monitoring. In Proceedings of the 2017 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSPIA ASC), Kuala Lumpur, Malaysia, 12–15 December 2017; pp. 686–691. [CrossRef]

3. Ghaderi, A.; Sanandaji, B.M.; Ghaderi, F. Deep Forecast: Deep Learning-based Spatio-Temporal Forecasting. In Proceedings of the International Conference on Machine Learning, Time Series Workshop, Sydney, NSW, Australia, 6–11 August 2017.

4. Silver, D.; Schrittwieser, J.; Simonyan, K.; Antonoglou, I.; Huang, A.; Guez, A.; Hubert, T.; Baker, L.; Lai, M.; Bolton, A.; et al. Mastering the game of Go without human knowledge. *Nature* **2017**, *550*, 354–359. [CrossRef]

5. Sugimori, H.; Hamaguchi, H.; Fujiwara, T.; Ishizaka, K. Classification of type of brain magnetic resonance images with deep learning technique. *Magn. Reson. Imaging* **2021**, *77*, 180–185. [CrossRef]

6. Sugimori, H.; Kawakami, M. Automatic detection of a standard line for brain magnetic resonance imaging using deep learning. *Appl. Sci.* **2019**, *9*, 3849. [CrossRef]

7. Zhou, X. Automatic Segmentation of Multiple Organs on 3D CT Images by Using Deep Learning Approaches. *Adv. Exp. Med. Biol.* **2020**, *1213*, 135–147. [CrossRef] [PubMed]

8. Yang, S.; Kweon, J.; Roh, J.H.; Lee, J.H.; Kang, H.; Park, L.J.; Kim, D.J.; Yang, H.; Hur, J.; Kang, D.Y.; et al. Deep learning segmentation of major vessels in X-ray coronary angiography. *Sci. Rep.* **2019**, *9*, 16897. [CrossRef]

9. Zhang, L.; Mohamed, A.A.; Chai, R.; Guo, Y.; Zheng, B.; Wu, S. Automated deep learning method for whole-breast segmentation in diffusion-weighted breast MRI. *J. Magn. Reson. Imaging* **2020**, *51*, 635–643. [CrossRef] [PubMed]

10. Kessler, D.A.; MacKay, J.W.; Crowe, V.A.; Henson, F.M.D.; Graves, M.J.; Gilbert, F.J.; Kaggie, J.D. The optimisation of deep neural networks for segmenting multiple knee joint tissues from MRIs. *Comput. Med. Imaging Graph.* **2020**, *86*, 101793. [CrossRef] [PubMed]

11. Cattagni, T.; Scaglioni, G.; Laroche, D.; Gremeaux, V.; Martin, A. The involvement of ankle muscles in maintaining balance in the upright posture is higher in elderly fallers. *Exp. Gerontol.* **2016**, *77*, 38–45. [CrossRef]

12. Belavý, D.L.; Miokovic, T.; Rittweger, J.; Felsenberg, D. Estimation of changes in volume of individual lower-limb muscles using magnetic resonance imaging (during bed-rest). *Physiol. Meas.* **2011**, *32*, 35–50. [CrossRef]

13. Tortorella, C.; Simone, O.; Piazzolla, G.; Stella, I.; Cappiello, V.; Antonaci, S. Role of phosphoinositide 3-kinase and extracellular signal-regulated kinase pathways in granulocyte macrophage-colony-stimulating factor failure to delay fas-induced neutrophil apoptosis in elderly humans. *J. Gerontol. Ser. A* **2006**, *61*, 1111–1118. [CrossRef]

14. Karamanidis, K.; Epro, G.; König, M.; Mersmann, F.; Arampatzis, A. Simplified Triceps Surae Muscle Volume Assessment in Older Adults. *Front. Physiol.* **2019**, *10*, 1299. [CrossRef]

15. Friedberger, A.; Figueiredo, C.; Bäuerle, T.; Schett, G.; Engelke, K. A new method for quantitative assessment of hand muscle volume and fat in magnetic resonance images. *BMC Rheumatol.* **2020**, *4*, 72. [CrossRef] [PubMed]
16. Kemnitz, J.; Baumgartner, C.F.; Eckstein, F.; Chaudhari, A.; Ruhdorfer, A.; Wirth, W.; Eder, S.K.; Konukoglu, E. Clinical evaluation of fully automated thigh muscle and adipose tissue segmentation using a U-Net deep learning architecture in context of osteoarthritic knee pain. *Magn. Reson. Mater. Phys. Biol. Med.* **2020**, *33*, 483–493. [CrossRef] [PubMed]

17. Ding, J.; Cao, P.; Chang, H.C.; Gao, Y.; Chan, S.H.S.; Vardhanabhuti, V. Deep learning-based thigh muscle segmentation for reproducible fat fraction quantification using fat-water decomposition MRI. *Insights Imaging* **2020**, *11*, 128. [CrossRef]

18. Gadermayr, M.; Li, K.; Müller, M.; Truhn, D.; Krämer, N.; Merhof, D.; Gess, B. Domain-specific data augmentation for segmenting MR images of fatty infiltrated human thighs with neural networks. *J. Magn. Reson. Imaging* **2019**, *49*, 1676–1683. [CrossRef] [PubMed]

19. Ghosh, S.; Ray, N.; Boulanger, P. A Structured Deep-Learning Based Approach for the Automated Segmentation of Human Leg Muscle from 3D MRI. In Proceedings of the 2017 14th Conference on Computer and Robot Vision (CRV), Edmonton, AB, Canada, 16–19 May 2017; pp. 117–123. [CrossRef]

20. Andrews, S.; Hamarneh, G. The Generalized Log-Ratio Transformation: Learning Shape and Adjacency Priors for Simultaneous Thigh Muscle Segmentation. *IEEE Trans. Med. Imaging* **2015**, *34*, 1773–1787. [CrossRef] [PubMed]

21. Baudin, P.-Y.; Azzabou, N.; Carlier, P.G.; Paragios, N. Automatic Skeletal Muscle Segmentation through Random Walks and Graph-Based Seed Placement. In Proceedings of the 2012 9th IEEE International Symposium on Biomedical Imaging (ISBI), Barcelona, Spain, 2–5 May 2012; pp. 1036–1039. [CrossRef]

22. Essafi, S.; Langs, G.; Deux, J.F.; Rahmouni, A.; Bassez, G.; Paragios, N. Wavelet-driven knowledge-based MRI calf muscle segmentation. In Proceedings of the 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, Boston, MA, USA, 28 June–1 July 2009; pp. 225–228. [CrossRef]

23. Abdullah, J.Y.; Rajion, Z.A.; Martin, A.G.; Jaafar, A.; Ghani, A.R.I.; Abdullah, J.M. Shape-based interpolation method in measuring intracranial volume for pre- and post-operative decompressive craniectomy using open source software. *Neurocirugia* **2019**, *30*, 115–123. [CrossRef]

24. Nordez, A.; Jolivet, E.; Südhoff, I.; Bonneau, D.; De Guise, J.A.; Skalli, W. Comparison of methods to assess quadriceps muscle volume using magnetic resonance imaging. *J. Magn. Reson. Imaging* **2009**, *30*, 1116–1123. [CrossRef]

25. Ahmed, I.; Ahmad, M.; Khan, F.A.; Asif, M. Comparison of deep-learning-based segmentation models: Using top view person images. *IEEE Access* **2020**, *8*, 136361–136373. [CrossRef]

26. Khan, Z.; Yahya, N.; Alsaih, K.; Ali, S.S.A.; Meriaveau, F. Evaluation of deep neural networks for semantic segmentation of prostate in T2W MRI. *Sensors* **2020**, *20*, 3183. [CrossRef]

27. Rundo, L.; Han, C.; Nagano, Y.; Zhang, J.; Hataya, R.; Militello, C.; Tangherloni, A.; Nobile, M.S.; Ferretti, C.; Besozzi, D.; et al. USE-Net: Incorporating Squeeze-and-Excitation blocks into U-Net for prostate zonal segmentation of multi-institutional MRI datasets. *Neurocomputing* **2019**, *365*, 31–43. [CrossRef]

28. Schlemper, J.; Oktay, O.; Schaap, M.; Heinrich, M.; Kainz, B.; Glocker, B.; Rueckert, D. Attention gated networks: Learning to leverage salient regions in medical images. *Med. Image Anal.* **2019**, *53*, 197–207. [CrossRef]

29. Yeung, M.; Sala, E.; Schönlieb, C.-B.; Rundo, L. Focus U-Net: A novel dual attention-gated CNN for polyp segmentation during colonoscopy. *Comput. Biol. Med.* **2021**, *137*, 104815. [CrossRef] [PubMed]

30. Chen, L.C.; Zhu, Y.; Papandreou, G.; Schroff, F.; Adam, H. Encoder-decoder with atrous separable convolution for semantic image segmentation. In *Lecture Notes in Computer Science*; Springer: Berlin/Heidelberg, Germany, 2018; pp. 833–851. [CrossRef]

31. Do, D.T.; Le, T.Q.T.; Le, N.Q.K. Using deep neural networks and biological subwords to detect protein S-sulfonylation sites. *Brief. Bioinform.* **2021**, *22*, bbaa128. [CrossRef]

32. Le, N.Q.K.; Huynh, T.T. Identifying SNAREs by Incorporating Deep Learning Architecture and Amino Acid Embedding Representation. *Front. Physiol.* **2019**, *10*, 1501. [CrossRef] [PubMed]

33. Miyachi, R.; Yamazaki, T.; Ohno, N.; Miyati, T. Relationship between muscle cross-sectional area by mri and muscle thickness by ultrasonography of the triceps surae in the sitting position. *Healthcare* **2020**, *8*, 166. [CrossRef] [PubMed]

34. Henninger, H.B.; Christensen, G.V.; Taylor, C.E.; Kawakami, J.; Hillyard, B.S.; Tashjian, R.Z.; Chalmers, P.N. The Muscle Cross-sectional Area on MRI of the Shoulder Can Predict Muscle Volume: An MRI Study in Cadavers. *Clin. Orthop. Relat. Res.* **2020**, *478*, 871–883. [CrossRef] [PubMed]

35. Akagi, R.; Takai, Y.; Ohta, M.; Kanehisa, H.; Kawakami, Y.; Fukunaga, T. Muscle volume compared to cross-sectional area is more appropriate for evaluating muscle strength in young and elderly individuals. *Age Ageing* **2009**, *38*, 564–569. [CrossRef]
<figure><img src="image_60.png" /></figure>
<figure><img src="image_61.png" /></figure>

Article

# Automated Breast Lesion Detection and Characterization with the Wavelia Microwave Breast Imaging System: Methodological Proof-of-Concept on First-in-Human Patient Data

Angie Fasoula <fnref n="1" />,*, Luc Duchesne <fnref n="1" />, Julio Daniel Gil Cano <fnref n="1" />, Brian M. Moloney <fnref n="2" />,<fnref n="3" />, Sami M. Abd Elwahab <fnref n="4" /> and Michael J. Kerin <fnref n="3" />,<fnref n="4" />

**Citation:** Fasoula, A.; Duchesne, L.; Gil Cano, J.D.; Moloney, B.M.; Abd Elwahab, S.M.; Kerin, M.J.
Automated Breast Lesion Detection and Characterization with the Wavelia Microwave Breast Imaging System: Methodological Proof-of-Concept on First-in-Human Patient Data. *Appl. Sci.* **2021**, *11*, 9998.
https://doi.org/10.3390/app11219998

Academic Editors: Leonardo Rundo, Carmelo Militello and Andrea Tangherloni

Received: 28 September 2021

Accepted: 19 October 2021

Published: 26 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_62.png" /></figure>

Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Microwave Breast Imaging (MBI) is an emerging non-ionizing imaging modality, with the potential to support breast diagnosis and management. Wavelia is an MBI system prototype, of 1st generation, which has recently completed a First-In-Human (FiH) clinical investigation on a 25-symptomatic patient cohort, to explore the capacity of the technology to detect and characterize malignant (invasive carcinoma) and benign (fibroadenoma, cyst) breast disease. Two recent publications presented promising results demonstrated by the device in this FiH study in detecting and localizing, as well as delineating size and malignancy risk, of malignant and benign palpable breast lesions. In this paper, the methodology that has been employed in the Wavelia semi-automated Quantitative Imaging Function (QIF), to support breast lesion detection and characterization in the FiH clinical investigation of the device, is presented and the critical design parameters are highlighted.

**Keywords:** breast cancer detection; microwave breast imaging; computer-aided diagnosis (CAD); first-in-human (FiH) study

## 1. Introduction

Microwave Breast Imaging (MBI) uses the scattering wave, or reflected wave, that arises from the contrast in dielectric properties between the various breast tissues, in the microwave frequency range [1]. The increased volume of water within the denser breast tissues is responsible for the detectable electromagnetic scattering associated with microwave imaging. The increase in sodium and water, particularly in-bound water within the tumor cells, is expected to lead to even greater conductivity and permittivity of the tumorous tissues [2,3]. Due to the dielectric contrast, back-scattered radar signals are physically generated, when the breast is illuminated with low-power electromagnetic waves in the microwave frequency range.

MBI has been investigated as a novel modality for the detection of breast disease, offering a non-ionizing, non-compressive approach [4–6] and as a potential diagnostic management strategy in the monitoring of neoadjuvant chemotherapy [7]. To date, a total of at least 10 MBI system prototypes have been employed in human subject tests, to investigate the clinical utility of MBI [8–11]. Despite encouraging clinical results being reported, several recurrent limitations, as outlined in [12], remain unresolved across most studies and justify further clinical research with alternative MBI systems, such as Wavelia.

<footnote n="1">¹ Medical Imaging Department, MVG Industries, 91140 Villejust, France; luc.duchesne@mvg-world.com (L.D.); julio_daniel.gil_cano@mvg-world.com (J.D.G.C.)</footnote>

<footnote n="2">² Department of Radiology, Galway University Hospital, Saolta University Healthcare Group, H91 YR71 Galway, Ireland; brianmoloney1@hotmail.com</footnote>

<footnote n="3">³ Discipline of Surgery, Lambe Institute for Translational Research, School of Medicine, National University of Ireland Galway, H91 TK33 Galway, Ireland; michael.kerin@nuigalway.ie</footnote>

<footnote n="4">⁴ Department of Surgery, Galway University Hospital, Saolta University Healthcare Group, H91 YR71 Galway, Ireland; sami.elwahab@hse.ie</footnote>

<footnote>* Correspondence: angie.fasoula@mvg-world.com</footnote>
Wavelia™ is an MBI system prototype, of 1st generation, which demonstrated the ability to detect dielectric contrast between tumor phantoms and synthetic fibroglandular tissue in preclinical studies [13] and has recently completed a First-In-Human (FiH) clinical investigation on a 25-symptomatic patient cohort, hosted in NUIG Clinical Research Facility Galway, Ireland. In this study (ClinicalTrials.gov NCT03475992), the Wavelia MBI system was evaluated in the clinical setting for the first time, using mammography as the reference conventional imaging modality and post-surgery histology data to assess the size of the cancers. Ultrasound, MRI and core biopsy data were also collected as a reference and were available as part of the patient’s standard of care. In this FiH study, Wavelia demonstrated the capacity to detect and approximate underlying breast abnormalities to the appropriate location, in patients with palpable biopsy-confirmed invasive carcinomas and benign breast lesions, such as cysts and fibroadenomas [12]. The device also demonstrated promising results in delineating the size and malignancy risk of the detected breast lesions [14].

The methodology that was employed in the Wavelia semi-automated Quantitative Imaging Function (QIF) during this FiH study, to support morphological breast lesion detection based on persistence, lesion sizing and lesion characterization in a low-dimensional feature space, spanning shape and texture-based features, is presented in this article.

## 2. Materials and Methods

The Wavelia MBI Quantitative Imaging Function (QIF) was initially conceived using experimental MBI datasets from anthropomorphic breast phantoms [13,15] and was further developed and configured following training on the available FiH patient datasets [12,16]. The MBI parametric radar image formation and clinical feature extraction are performed offline at this stage of development of Wavelia.

### 2.1. Wavelia MBI: Parametric Radar Image Formation

The Wavelia MBI system operates using 18 antennae arranged in a circle in a horizontal plane outside a cylinder. With the patient lying in the prone position, one breast is submerged at a time into the cylinder, which is filled with a creamy transition liquid. The liquid has dielectric properties similar to the ones of the human skin within the microwave frequency spectrum, thus favoring the penetration of the electromagnetic waves in the breast. The device illuminates the breast using low-power electromagnetic waves in the frequency range [0.5–4] GHz. The probe array moves vertically below the examination table and illuminates the breast at regular intervals of 5 mm. Coronal sections of the breast, of a given thickness (10 mm) are generated using the MBI data at each vertical scan position of the probe array. Partially overlapping consecutive coronal breast sections, formed per azimuthal sector of illumination based on multi-static radar detection technology, are integrated to form a 3D MBI image of the dielectrically contrasted interior breast tissues.

As specified in prior publications on Wavelia MBI [13,15], the multi-static radar imaging algorithm, which is employed for MBI image formation, is the Time-Reversal Multiple Signal Classification (TR-MUSIC) algorithm, which was originally conceived for the detection of obscured radar targets in heavily cluttered environments [17]. The intensity of the TR-MUSIC images gets maximized in the imaging pixels where the MBI sensor array illumination vector is more orthogonal to the noise subspace; thus, the image intensity is indicative of the probability for the presence of a dielectrically contrasted scatterer on each pixel of the image. The noise subspace is estimated at each frequency, by means of decomposition and analysis of the Multi-Static Frequency Response Matrix (MFRM) of the imaging array. The illumination vector of the imaging array, at each pixel p of the imaging scene and each frequency f, is defined as:

$$G_{\text{sec t}}(p, f) = [ \mathbf{g}_0(\mathbf{p}_{\text{TRx}_{\text{sec t},1}}, p, f) \quad \mathbf{g}_0(\mathbf{p}_{\text{TRx}_{\text{sec t},2}}, p, f) \quad \dots \quad \mathbf{g}_0(\mathbf{p}_{\text{TRx}_{\text{sec t},N_s}}, p, f) ]^T \quad (1)$$

with:

$$\mathbf{g}_0(\mathbf{p}_{\text{TRx}_{\text{sec t},1}}, \mathbf{p}, \mathbf{f}) = j \cdot \mathbf{H}_0^{(1)} \left( \frac{2\pi f}{c_0} \cdot \left( \sqrt{\mathbf{e}_{\mathbf{r},\text{trans}}(\mathbf{f})} \cdot \dot{\mathbf{d}}_{\text{OutOfBreast},i,p} + \sqrt{\mathbf{e}_{\mathbf{r},\text{InBreast}}(\mathbf{f})} \cdot \dot{\mathbf{d}}_{\text{InBreast},i,p} \right) \right) \quad (2)$$
the assumed underlying ElectroMagnetic (EM) wave propagation model for the antenna element at position $p_{TRX_{sec,i}}$, $H_0^{(1)}$ the Hankel function of 1st kind and 0th order, $c_0$ the speed of light in vacuum, $e_{r,trans}(f)$ the known permittivity of the transition liquid at the frequency $f$, $\hat{d}_{OutOfBreast,i,p}$ an estimate of the distance travelled by the EM wave in the transition liquid up to reaching the imaging pixel $p$, $\hat{d}_{InBreast,i,p}$ an estimate of the distance travelled by the EM wave within the breast up to reaching the pixel $p$, and:

$$ \hat{e}_{r,InBreast}(f) = \left( p c_{fib} \cdot \hat{e}_{r,fibroglandular}(f) + (1 - p c_{fib}) \cdot \hat{e}_{r,adipose}(f) \right) \cdot 10^{-2} \quad (3) $$

the average permittivity of the background healthy tissues of the breast, defined as a weighted average (weighting by $pc_{fib}$) of the adipose tissue and fibro-glandular tissue “mean” dielectric properties, as derived by Sugitani et al. [18].

The breast external envelope is reconstructed at first using the Wavelia MBI scan data. The geometry is exploited to split the imaging scene in “Out of breast” and “In breast” segments and further estimate $\hat{d}_{OutOfBreast,i,p}$ and $\hat{d}_{InBreast,i,p}$ for each transceiver i and each pixel p in the imaging scene.

The critical elements of the tailored implementation of TR-MUSIC in the Wavelia QIF are summarized below and in the flowchart of Figure 1.

<figure><img src="image_63.png" /><figcaption>Figure 1. Wavelia MBI radar image formation: the critical parameters.</figcaption></figure>

* Sectorization of the imaging scene: Considering the high level of heterogeneity of the breast tissues, but also the potentially irregularly shaped breast tumours, MBI image formation is performed using sectorized subsets of the circular sensor network of the Wavelia MBI system, at each vertical scan position. Both the physical size of the sensor sub-array (sector) and the number $N_s$ of sensors used for the elementary sub-image formation is critical to the achievable performance of the MBI system, in terms of unambiguous detection and valid characterization of breast lesions. $N_s = 6$ was fixed in the Wavelia QIF for the FiH clinical investigation.

* Sensor fidelity zone setting: Bounding the portion of the imaging scene, being efficiently illuminated by each sensor sub-array. In the current implementation of the
Wavelia QIF, a spatial filter activates the contribution of a given antenna to a given pixel p in the imaging scene, only if the Euclidean distance between the pixel and the phase centre of the antenna is inferior to a pre-set value $d_{max}$.

* Automated, data-driven, frequency selection for imaging: The TR-MUSIC imaging is performed per frequency. In the Wavelia QIF, a limited number of frequency points is automatically selected for imaging in each azimuthal sector of illumination, based on the information content of the MFRM at each frequency. The employed frequency selection criterion was previously defined in [15].

* pc_fib parameter setting in multiple search ranges: A large variability exists in the dielectric properties of each breast tissue type over the population, as demonstrated by multiple studies involving ex-vivo dielectric measurements of a large sample of excised breast tissues [18–21]. Considering that the full dielectric map of each breast cannot become practically available, data-driven techniques are employed in the Wavelia QIF to deduce the unknown dielectric properties of the healthy breast tissue in each breast, by assessing the pc_fib parameter. The pc_fib parameter, which is involved in the formulation of the illumination vector of the MBI sensor array, is physically associated with the percentage of fibro-glandular tissue along the propagation path within the breast, from a given transmitting antenna to the interrogated imaging pixel and back to a given receiving antenna, as defined in Equations (1)–(3). The Wavelia QIF generates a set of parametric MBI radar images under various assumptions on pc_fib. The generated set of parametric images is further evaluated in terms of focusing, using the image curvature [22,23] as a focusing quality measure. To better handle the heterogeneity of the breast and potentially better reveal the non-uniform angular response of the breast lesions to MBI, the pc_fib parameter setting is performed independently in each azimuthal imaging sector, while employing multiple search ranges. In the Wavelia QIF, X1 wide and X2 narrow pc_fib parameter search ranges are systematically employed for image formation, thus a total number of X = (X1 + X2) MBI images are formed per patient’s breast.

The wide pc_fib search ranges result in 3D MBI images including the most complete representations of the detected breast lesion shape. The narrow pc_fib search ranges are expected to lead to partial representations of the detectable breast lesions. The X = 5 pc_fib search ranges, which were systematically employed during the FiH clinical investigation of Wavelia, are listed below:

* Wide pc_fib search range #1 (W1): pc_fib ∈ [10 20 30 40 50 60]%

* Wide pc_fib search range #2 (W2): pc_fib ∈ [20 30 40 50]%

* Narrow pc_fib search range #1 (n1): pc_fib ∈ [10 20]%

* Narrow pc_fib search range #2 (n2): pc_fib ∈ [30 40]%

* Narrow pc_fib search range #3 (n3): pc_fib ∈ [50 60]%

As explained in the next subsection, the persistent presence of a Region-Of-Interest (ROI) in the set of X MBI parametric radar images of a given breast is further exploited, to support the association of automatically extracted ROIs with breast lesions and validate their reporting for clinical analysis.

## 2.2. Morphological MBI Image Post-Processing: Breast Lesion Detection Based on Persistence

Automated breast lesion detection is performed in the Wavelia QIF by means of morphological post-processing of the set of parametric radar images, which are formed with the employment of X = 5 search ranges for the pc_fib parameter. Automated segmentation of ROIs and association, or not, to a breast lesion is based on morphological properties (solidity and volume) of the ROI and its persistence on the set of parametric radar images, which is evaluated by means of spatial clustering. The persistent visibility of a ROI over multiple pc_fib search ranges is indicative of the association of the ROI with a physical object (breast lesion) in the MBI image. On the other hand, the presence of a ROI in the minority of the pc_fib search ranges under test is indicative of it being associated with an imaging artefact. This setting has been inspired by the “breast mass” definition for
mammography, as the space-occupying a 3D lesion seen in two different projections [24]. To
the authors' knowledge, no such breast lesion detection method, based on persistence in a
set of parametric images carrying redundant information content, has ever been integrated
into any of the state-of-the-art MBI systems before. A second novel element of the proposed
method is the coupling of morphological properties (solidity) with the notion of persistence
to validate a ROI detection.

The automated breast lesion detection method, as designed and integrated into the Wavelia QIF for the FiH clinical investigation of device prototype #1, is outlined below and in the block diagram in Figure 2.

<figure><img src="image_64.png" /><figcaption>**Figure 2.** Wavelia MBI: morphological image post-processing for breast lesion detection based on persistence.</figcaption></figure>

1. Iterative Image threshold setting: The following operation is first-of-all performed to set a threshold for the raw 3D MBI image, with no a priori available on how normal breast tissue is represented in this type of image:

a. Progressive increase in the image threshold, starting from the null threshold,

b. At each iteration, identification of the "connected" objects in the thresholded image,

c. Threshold setting based on the maximal accepted volumetric size of "connected" objects, potentially defining a breast mass in the image (default value 3 cm³, in this implementation).

2. Semi-automated ROI extraction, based on morphological properties: The connected objects to be retained as ROIs, are defined based on a set of user-defined characteristics, including:

a. Volumetric size: all the small objects, of volume inferior to 1 cm³, are removed from the FiH clinical data analysis, considering the status of the Wavelia system prototype #1 in terms of minimum size of detectable lesions.

b. Solidity: this structural feature measures the density (or convexity) of an object. A measure of solidity can be obtained as a ratio of the volume of the object to the volume of a convex hull of the object. A value of 100% indicates a solid object,
and a value less than 100% indicates an object having an irregular boundary or containing holes. All connected objects with solidity >30% have been ultimately retained, for the data analysis of the FiH clinical investigation of Wavelia.

c. **Intensity Contrast:** In case of ambiguity (i.e., multiple connected objects to be retained in a single 3D MBI image, based on the volume and solidity criteria), each connected object is retained as ROI, only if it is associated with the maximal intensity in the image and is minimally contrasted (by at least 5%) against the intensity of all the "competing" connected objects in the image.

3. **Refinement of the ROI segmentation:** An Active Contour segmentation module (Chan-Vese algorithm for segmentation without edges [25]) has been configured and employed to refine the contour of the extracted ROIs, in order to enable a more valid characterization of the lesions in terms of shape and texture, as defined in the following sub-section. This module may have a critical impact, especially in the case of small lesions, which are low contrasted against the background healthy breast tissue.

4. **Spatial clustering of the ROIs:** which have been extracted on the set of X MBI images, formed with varying pc_fib search ranges, is performed. In the current implementation of the Wavelia QIF, the Euclidean distance between the centroids of two ROIs is required to be shorter than $d_R = 1$ cm for the ROIs to be associated with the same lesion in the breast.

5. **Persistence over varying pc-fib search range:** Five pc-fib search ranges (2 wide and 3 narrow) are systematically used during the FiH clinical investigation of Wavelia to generate parametric 3D MBI images of each patient's breast, as earlier stated. Detection of a breast lesion in a minimum number of parametric images (3 out of the 5 pc_fib search ranges) is required for a ROI to be considered persistent and validated.

### 2.3. Combination of 3D Shape Descriptors and Texture Features for Breast Lesion Characterization with Microwave Breast Imaging (MBI)

Apart from using reflected microwave energy to reconstruct images of the breast, additional information on the size, shape, and surface texture can be extracted and potentially exploited for discrimination between benign and malignant breast lesions using microwaves [1]. Malignant tumors usually present the following characteristics: irregular and asymmetric shapes, blurred boundaries (lack of sharpness), rough and complex surfaces with spicules or micro-lobules, non-uniform permittivity, and irregular tissue density. Conversely, benign tumors tend to have the following characteristics: well-circumscribed contours, compactness, and a smooth surface. Previous research works on breast lesion characterization/classification with MBI [26–30] considered principally the MBI received signals as input to a classifier, with or without dimensionality reduction. These state-of-the-art research works [26–30] have been based on simulated datasets and/or simplified experimental setups; no evaluation of such methods on patient clinical datasets has been published to date. Among the state-of-the-art MBI prototypes which have been tested on clinical datasets, two of them published studies on breast lesion classification with MBI. Early concept work on the exploitation of the pattern of the frequency-domain Radio-Frequency (RF) responses of the ROIs representing the breast lesions in the MBI image was published in [31] for the MARIA M5 [4] MBI system. For the MammoWave MBI system, machine learning methods were employed with raw received signals in the frequency domain to classify them as healthy or non-healthy responses [32].

In the Wavelia QIF a module is integrated for the characterization of the ROIs which have been prior detected and validated based on morphological properties and their persistence, and thus associated with breast lesions. This module includes the following operations:

* **Breast Lesion sizing:** by means of fitting an ellipsoid to the ROI associated with the persistent lesion detection, in the 3D MBI images that have been generated by applying either of the two wide pc_fib search ranges. The greatest linear dimension
of the lesion is defined as the length of the longest axis of the fitted ellipsoid. This definition is compatible with the conventional method that is applied for sizing breast abnormalities based on 2D mammography and ultrasound images [24,33]. During the FiH clinical investigation, the Wavelia MBI system showed promise for measuring lesion size with a more favorable linear trend between MBI and post-surgery histological lesion size, compared to the results obtained for conventional imaging [14]. Two challenging patient cases in terms of breast lesion sizing are indicatively discussed in Results Section 3.2, to better highlight the status of the MBI lesion sizing method, as integrated into the current version of the Wavelia QIF.

* Malignant-to-benign Breast Lesion labelling is performed in the Wavelia QIF using a combination of 3 features, extracted from morphologically validated ROIs in the MBI images. The selected features include: a shape descriptor [34], a Gray-Level Co-occurrence Matrix (GLCM) texture feature [35] and a Neighborhood Gray Tone Difference Matrix (NGTDM) texture feature [36]. The 3 features which were selected in the Wavelia QIF implementation for the FiH study of the investigational device are more specifically the following:

* Shape descriptors—Solidity: This feature measures the density, or the convexity, of an object. It is computed as the ratio of the volume of the object to the volume of the convex hull of the object, as illustrated in Figure 3. Breast lesion scoring, in terms of risk for malignancy, is routinely based on visual inspection and evaluation of the shape and margins of the imaged breast lesion, as per BIRADS [24,33]. Shape descriptors have been earlier considered for breast lesion classification with mammography [34,37] and ultrasound [38].

* GLCM texture—Correlation: The GLCM texture features measure the spatial relationship between pixels per specific directions, thus highlighting the properties of uniformity, homogeneity, randomness, and linear dependency of the image [35]. More specifically, the “correlation” feature varies between 0 (uncorrelated) and 1 (perfectly correlated), showing the linear dependency of gray level values to their respective voxels, as graphically illustrated in Figure 3.

* NGTDM texture—Busyness: The NGTDM texture features measure the spatial relationship among three or more pixels neighborhood, closely approaching the human perception of the image [36], as graphically illustrated in Figure 3. More specifically, for the “busyness” feature, a high value indicates a “busy” image, with rapid changes of intensity between pixels and its neighborhood.

<figure><img src="image_65.png" /><figcaption>**Figure 3.** The 3-d feature vector employed in Wavelia for malignant-to-benign breast lesion discrimination.</figcaption></figure>

Texture-based features have been earlier considered in Radiomics Research for cancerous lesions identification on CT, PET and MRI images [39–41]. Breast lesion classification, employing texture features on multi-parametric breast MRI images has also been introduced in the state-of-the-art [42].
In contrast to the Radiomics Research studies, which suggest the employment of high-dimensional feature vectors (typical size > 30) [39–41], appropriate feature selection has been considered in the Wavelia QIF, to achieve malignant-to-benign lesion separability in a feature space of low dimensionality. To the authors’ knowledge, no shape-based or texture-based feature extraction from Microwave Breast Images (MBI) has ever been considered in the past.

The 3-dimensional (3-d) lesion feature vector data [Solidity; Correlation; Busyness] is exploited in a malignant-to-benign breast lesion classification framework in the Wavelia QIF. A 2-class discrimination problem is defined, with: (i) Class #1: Malignant breast lesions, and (ii) Class #2: Benign breast lesions. Two classifiers have been trained in this 3-d feature space. The two classifiers, i.e., a Naïve Bayesian (NB) classifier and a Quadratic Discriminant Analysis (QDA) classifier, were selected such that their decision hypersurface partitions the 3-d feature space in two disjoint and continuous manifolds (malignant lesions subspace vs. benign lesions subspace).

In the Wavelia FiH clinical investigation [12], female patients were recruited from the symptomatic unit to one of three groups: Biopsy-proven breast cancers (Group-1), unaspirated cysts (Group-2) and biopsy-proven benign breast lesions (Group-3). For the training of the 2 classifiers:

* The Group-1 patient datasets were labelled as Class #1.

* The Group-2 and Group-3 patient datasets were labelled as Class #2.

A total of 25 patients underwent MBI in this FiH study. Of these, 24 were included in the final data analysis (11 Group-1, 8 Group-2 and 5 Group-3 patients). The patient who was excluded from the final analysis was a patient who presented with a palpable lump which was determined to be normal breast tissue, and who also had small, scattered, cysts appearing in a different breast quadrant.

Given the small total number of analyzed patients, the number of training data samples which was extracted from each patient dataset equals the number of pc_fib search ranges for which the detection of each breast lesion was morphologically validated based on persistence. This implies that each detected breast lesion was represented by 3–5 points in the 3-d feature space, as depicted in Figure 4c. The confusion matrix and classification loss were estimated for the two trained classifiers by means of 10-fold cross-validation (i.e., 10 partitions of the full dataset in disjoint training and test datasets) to evaluate the potential for discrimination between malignant and benign breast lesions with Wavelia MBI. The confusion matrices and the decision surfaces are shown for the two classifiers in Figure 4a–c. This proof-of-concept FiH patient dataset suggested the good potential for discrimination between malignant and benign lesions in the defined 3-d feature space. The two classifiers demonstrated very comparable performance and associated classification loss 11.5–12.5%, as depicted in Figure 4.
<figure><img src="image_66.png" /><figcaption>**Figure 4.** Average confusion matrices for the 2 trained classifiers, estimated with 10-fold cross-validation, (a) NB classifier, (b) QDA classifier, (c) the training dataset and the decision hypersurfaces of the 2 classifiers, (d) Partitioning of the 3-d feature space illustrated with 3 cuts for the QDA classifier.</figcaption></figure>

## 3. Results

### 3.1. Semi-Automated Breast Lesion Detection Based on Persistence

The Wavelia MBI algorithm for morphological breast lesion detection based on persistence has been specified in the previous section. The lesion persistence is assessed over a set of MBI images that were generated under varying assumptions on the dielectric properties of the healthy tissue of the breast (varying pc_fib parameter search ranges). Lesions that are morphologically detected in at least 3 out of the 5 pc_fib search ranges under evaluation are considered persistent and validated. The principle of the breast lesion detection method is illustrated in Figures 5–8 on two indicative patient test cases.

* Patient 032: Group-1: 54-years old patient with an Invasive Lobular Carcinoma (ILC) of size 30 mm (MRI data) at the 12 o’clock position of the Right Breast. Breast density: BIRADS Category c, Volumetric Breast Density (VBD) = 13.3%.

* Patient 031: Group 3: 38-years old patient with a Fibroadenoma of size 19 mm (Ultrasound) in the lateral Left Breast. Breast density: BIRADS Category c, VBD = 10.8%.

The achieved persistence level of each breast lesion on MBI may vary depending on the histological type of the lesion and the density of the breast. It is interesting to note that the 30 mm ILC of Patient 032, which was not clearly visible on both mammogram and ultrasound, was persistent at 60% (i.e., 3 out of the 5 raw MBI images formed with the employment of distinct pc_fib search ranges) with MBI. More than a single dominant ROIs were visible on the raw MBI images. In Figure 6, the ROI which was extracted and
validated in terms of morphological properties, sufficient intensity contrast against the other competing ROIs in each image and persistence over varying pc_fib search ranges is presented encircled in the 3 out of the 5 raw MBI images in which it was detectable. As illustrated in Figures 7 and 8, the fibroadenoma of Patient 031 was persistent at 100% (i.e., 5 out of the 5 raw MBI images formed with the employment of distinct pc_fib search ranges) and was predominantly visible with MBI. Both patients had comparably dense breasts (P032: VBD = 13.3%, P031, VBD = 10.8%); however, the difference in terms of consistency of the two lesions may have been the principal reason for the distinct level of persistence of the response of the two lesions to MBI. The MBI scan datasets for the two patients have been processed using the same configuration of the Wavelia QIF.

In future upgraded implementations of the Wavelia QIF, both the persistence level over the varying assumption of the dielectric properties of the healthy tissue in the breast (varying pc_fib search ranges), but also the presence of a single dominant ROI or various competing ROIs in the image, may serve to define a confidence level for each MBI lesion detection, to better support the diagnosis.

<figure><img src="image_67.png" /><figcaption>**Figure 5.** Patient 032, ILC in the Right Breast: (a) Bilateral mammogram Cranio-Caudal (CC) view, (b) Bilateral mammogram Medio-Lateral Oblique (MLO) view, (c) Ultrasound scan, Right Breast, (d) MRI scan, bilateral axial image and sagittal image of the Right breast (e) MBI test results, Right Breast.</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>**Figure 6.** Patient 032, Right breast: ILC morphological detection based on persistence, with Wavelia #1 MBI.</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>**Figure 7.** Patient 031, Fibroadenoma in the Left Breast: (a) Bilateral mammogram Cranio-Caudal (CC) view, (b) Bilateral mammogram Medio-Lateral Oblique (MLO) view, (c) Ultrasound scan, Left Breast, (d) MBI test results, Left Breast.</figcaption></figure>
<figure><img src="image_70.png" /><figcaption>**Figure 8.** Patient 031, Left Breast: Fibroadenoma morphological detection based on persistence, with Wavelia MBI.</figcaption></figure>

## 3.2. Breast Lesion Sizing: Correlation with Conventional Imaging and Post-Surgery Histology

For two of the cancer patients, post-surgery histological analysis of the excised tumor demonstrated total tumor sizes which were much larger than the invasive tumor size. This was the case with Patient 002 and Patient 029, as reported in Table 1. The conventional imaging data (Mammography, Ultrasound) and the MBI imaging test results are depicted in Figure 9 for the case of Patient 002 and in Figure 10 for Patient 029.

It is interesting to observe in Table 1 and in Figure 9 that for the Patient 002 case the MBI lesion size estimate varies considerably depending on the pc_fib search range. Maximal linear dimension [34–51] mm, overestimated against the conventional imaging but better fitting to the total tumor size as confirmed with post-surgery histological analysis of the excised tumor, was retrieved with MBI for this lesion. For a subset of 3 out of the five 5 pc_fib search ranges being systematically evaluated in the Wavelia QIF, the irregularly shaped finding of the MBI system extended over a large volume (maximal linear dimension = 51 mm), including the core of the invasive tumor, as identified at triple assessment. By comparison with the patient's mammograms, it was deemed reasonable to consider that the Wavelia MBI system detected either the total tumor, or the invasive tumor and a concentration of fibro-glandular tissue adjacent to it. For the second subset of pc_fib search ranges, the size of the MBI detection was smaller, and its location seemed to correlate closely with the invasive tumor site. Due to the uncertainty on the orientation and the deformability of the patient's breast during the MBI scan, inaccuracies in the 3D reconstruction and localization of the tumor may arise when compared to conventional imaging data. This difficulty is not considered to be MBI-specific though. The registration of multi-modality imaging data of any kind, in the case of soft and deformable organs, like the breast, is a challenging task due to variations in the natural suspended position of the breast in the upright, supine and prone position.

<table><caption>**Table 1.** Breast lesion size: maximal linear dimension [mm].</caption><thead><tr><th rowspan="2">Patient ID</th><th colspan="2">Post-Surgery Histology</th><th colspan="3">MBI</th><th colspan="3">Conventional Imaging</th></tr><tr><th>Invasive Tumor</th><th>Total Tumor</th><th>W1</th><th>W2</th><th>Max.</th><th>Mammography</th><th>Ultrasound</th><th>Max.</th></tr></thead><tbody><tr><td>P002</td><td>20</td><td>40</td><td>51</td><td>34</td><td>51</td><td>25</td><td>15</td><td>25</td></tr><tr><td>P029</td><td>22</td><td>35</td><td>24.5</td><td>19.8</td><td>24.5</td><td>20</td><td>37</td><td>37</td></tr></tbody></table>
<figure><img src="image_71.png" /><figcaption>**Figure 9.** Patient 002, Invasive Ductal Carcinoma (IDC) at the 3 o'clock position of the Right breast (Breast density: BIRADS Category c, VBD = 8.5%): Bilateral Mammogram, Ultrasound and MBI test results of the Right Breast.</figcaption></figure>

Patient 029 was a patient with very dense breasts (VBD = 15.4%). It is clear in Figure 10 that the delineation of the margins of the tumor was not evident on the mammogram. In the radiology report, the presence of a 20 mm spiculated mass in the lower outer quadrant of the Left breast and calcifications extending anteriorly and medially from the mass and measuring up to (42 mm) × (47 mm) was reported. Ultrasound scan of the Left breast highlighted a 23 mm irregular hypoechoic mass in the lower outer quadrant (concurring with the invasive tumor), and a smaller node of 7 mm with indeterminate appearance, immediately superior to the mass. The total inclusive diameter of both lesions was reported to be 37 mm in the craniocaudal direction, concurring with the total tumor
size. MBI highlighted the presence of two persistent ROIs. The two ROIs are clearly visible in the raw MBI image shown in Figure 10. The ROI which was morphologically validated was located in the lower outer quadrant of the breast, had a maximal linear dimension of 24.5 mm and was associated with the invasive tumor in this analysis. The second ROI of volumetric size > 3 cm³ and rather low solidity (≈0.4) was present in the upper breast and could be, interestingly, associated with the extended zone of calcifications, as reported on the patient's mammogram. In the current version of the Wavelia QIF, the ROI definition is based on the notion of pixel connectivity, thus "discontinuous" constellations such as the one highlighted on ultrasound for this patient case and concurring with the total tumor size (as confirmed with post-surgery histology) could not be revealed. This patient case represents a limitation, which may be addressed in subsequent versions of the Wavelia semi-automated lesion sizing method.

<figure><img src="image_72.png" /><figcaption>**Figure 10.** Patient 029, IDC in the Left breast (Breast density: BIRADS Category c, VBD = 15.4%): Bilateral Mammogram, Ultrasound Scan and MBI test results of the Left Breast.</figcaption></figure>

### 3.3. Discrimination between Malignant and Benign Breast Lesions in a 3-d Feature Space (Shape-Based and Texture-Based Feature Employment)

While three narrow pc_fib search ranges are systematically used, together with the two wide pc_fib search ranges, to analyze the persistence of radar echoes for lesion detection, for the characterization of the detected lesions, i.e., sizing, shape and texture analysis, the wide pc_fib search ranges are mostly adequate to be employed, as they are expected to be associated with the most complete representations of the lesions in the available set of MBI images. In the course of the FiH clinical investigation of Wavelia, mapping of the wide pc_fib search range detections in the 3-d feature space was performed and the posterior probability for each detection to be associated to "Class #1 = Malignant lesion" (i.e., probability of malignancy) was computed, as predicted by the trained QDA classifier.

If a breast lesion was detected and validated in both wide pc_fib search range MBI images, two probabilities of malignancy were reported for the breast lesion, as depicted in
Figure 11. The maximal probability of malignancy was ultimately considered to represent a unique MBI classification score for the lesion in the data analysis. Patient 027 was the only Group-1 patient (IDC) for whom the probability of malignancy was inferior to 50%, for both wide pc_fib search ranges. Patient 029 (IDC) was an ambiguous case, with a probability of malignancy 16.5% and 55.7% for the 2 wide pc_fib search ranges, correspondingly. As depicted in Figure 10, this was a patient with very dense breasts, thus rendering the ROI delineation sensitive to the specific parameterization of the MBI morphological detector (Wavelia QIF), in its current version. The probability of malignancy was superior to 95% for all the other Group-1 patients (invasive carcinomas), inferior to 38% for all the Group-3 lesion detections (biopsied benign lesions) and inferior to 14% for all the Group-2 lesion detections (cysts).

<figure><img src="image_73.png" /><figcaption>**Figure 11.** Posterior probability of malignancy (classification score) per patient case in the Wavelia FiH study.</figcaption></figure>

In Figure 12, four patient cases are used to illustrate the impact of the three selected features (solidity, correlation, busyness) on the MBI lesion classification score. The MBI morphological images (i.e., ROIs detected and validated based on morphological properties and persistence) are superimposed with the outer surface of the breast, as reconstructed using the auxiliary Wavelia Optical Breast Contour Detection (OBCD) subsystem, in Figure 12a for the four patient test cases. The Wavelia OBCD subsystem which is employed to reconstruct the external surface of the breast with high resolution, based on optical data, was earlier introduced in [12,13,16]. This superposition serves to better highlight the location of the MBI breast lesion detection with reference to the nipple of the breast, which is visible in the OBCD reconstruction. Ultrasound images of the four patient test cases are included in Figure 12b, for a straightforward comparison with the MBI findings, both in terms of the morphology of each lesion and its localization in the breast. Mapping of the four breast lesions (1 IDC, 1 ILC, 1 fibroadenoma and 1 cyst) on the 3-d feature space of Wavelia MBI, together with the QDA decision surface, are shown in Figure 12c. The probability of malignancy, as predicted by the trained QDA classifier, is also annotated for each of the four lesions. The values of the three features and the associated probabilities of malignancy are reported in Table 2. The morphological detection (W1 or W2 pc_fib search range) which was associated with the highest probability of malignancy, has been used to represent each patient test case in Figure 12, and in Table 2.
<figure><img src="image_74.png" /><figcaption>**Figure 12.** Breast Lesion Characterization in the Wavelia 3-d feature space—Four illustrative patient test cases: (a) morphologically validated persistent ROI detections with Wavelia MBI, (b) Ultrasound images, (c) mapping in the Wavelia 3-d feature space, probability of malignancy and QDA decision surface.</figcaption></figure>

<table><caption>**Table 2.** Wavelia MBI lesion feature values and derived probability of malignancy.</caption><thead><tr><th rowspan="2">Patient ID</th><th rowspan="2">Lesion Histological Type</th><th colspan="3">Features</th><th rowspan="2">Probability of Malignancy [%]</th></tr><tr><th>Solidity</th><th>Correlation</th><th>Busyness</th></tr></thead><tbody><tr><td>P010</td><td>Invasive Ductal Carcinoma</td><td>0.529</td><td>0.496</td><td>0.041</td><td>99.95%</td></tr><tr><td>P032</td><td>Invasive Lobular Carcinoma</td><td>0.661</td><td>0.556</td><td>0.074</td><td>100%</td></tr><tr><td>P031</td><td>Fibroadenoma</td><td>0.643</td><td>0.403</td><td>0.018</td><td>34.03%</td></tr><tr><td>P040</td><td>Simple Cyst</td><td>0.816</td><td>0.422</td><td>0.02</td><td>7.55%</td></tr></tbody></table>
It is interesting to observe the following:

* the clear differentiability of the simple cyst (Patient 040) in terms of higher solidity,

* the fibroadenoma and the two cancerous lesions had similar solidity levels, however clear distinction was achieved between the benign and the two malignant lesions in terms of texture features,

* both the correlation and the busyness feature values were slightly increased in the case of the malignant lesions, with the increase being more notable on the busyness feature,

* substantially increased busyness value was associated with more heterogeneous lesion patterns, which may be interpretable as being indicative of distributed non-mass like ILC's, such as in the Patient 032 case.

This illustration highlights the physical reasoning behind the selection of the three specific features, for benign-to-malignant MBI lesion classification, based on shape and texture in the Wavelia QIF. It also serves to justify the achievable level of separability with this small proof-of-concept FiH dataset, while working with continuous subspace manifolds and very simple classifier models. The potential generalization of the above findings is intended to be confirmed with future clinical investigations, involving larger and more diverse patient datasets. Expansion of the feature space to include additional dimensions (features), supporting the generalization of the above findings on larger patient cohorts, will be also evaluated further during the development of Wavelia.

## 4. Discussion and Conclusions

In this article, the methodology that was employed in the Wavelia semi-automated Quantitative Imaging Function (QIF) during the FiH clinical investigation of the device, to support morphological breast lesion detection based on persistence, lesion sizing and lesion characterization in a low-dimensional feature space, spanning shape and texture-based features, has been outlined and the critical design parameters highlighted.

Semi-automated breast lesion detection using morphological post-processing of a set of parametric radar images, which are formed with Wavelia MBI under varying assumptions on the dielectric properties of the healthy background tissue of the breast was introduced in the Wavelia QIF. Automated segmentation of ROIs and association, or not, to a breast lesion is based on morphological properties (solidity, volume) of the ROI and its persistence on the set of parametric images, as evaluated based on spatial clustering. The novelty of the proposed method lies in the exploitation of the notion of persistence and its combination with the solidity feature to support validation of ROI detection in MBI images.

A methodology for malignant-to-benign breast lesion discrimination, based on mapping in a low-dimensional feature space, which spans both shape-based features (solidity) and texture-based features (correlation, busyness), and training of a Naïve Bayes (NB) and Quadratic Discriminant Analysis (QDA) classifier was also introduced in the Wavelia QIF. An interesting level of separability between malignant and benign breast lesions was achieved, with a classification loss of 11.5% estimated with 10-fold cross-validation for the trained QDA classifier. This is a result pending to be verified, reproduced, and validated with larger datasets in future clinical investigations of Wavelia. While extensive research work is already published on Radiomics applied to the well-established breast imaging modalities, to our knowledge, it is the first time that shape descriptors and texture-based features are computed for ROIs extracted from MBI images, to support breast lesion characterization and malignant/benign lesion labelling.

In this FiH study, which was conducted in 25 patients, the Wavelia #1 prototype system demonstrated the preliminary potential to detect and discriminate between malignant and benign palpable breast lumps, the imaging procedure had no safety issues and patients reported a favorable experience of the MBI scan. Although the number of subjects included in the FiH study was small and was not intended to permit a clinically meaningful statistical analysis, the promising findings from this study provided initial data to support the valid clinical association of the technology and warranted the preparation of further
clinical investigations, with an upgraded prototype version of the Wavelia system (Wavelia #2) and its semi-automated QIF, to progressively address the identified technological challenges. Larger and more diverse patient datasets are needed to validate these findings and delineate the cases where the Wavelia MBI modality may offer a beneficial adjunct to current diagnostic protocols.

For the first conception of the Wavelia QIF and its feasibility testing for the first time in the clinical setting, simplifications were imposed to the FiH data processing for the analysis to become more straightforward. The two most important limitations of this analysis, which are intended to be loosened in future clinical investigations for the Wavelia QIF to become sufficiently realistic and clinically meaningful, are discussed in the two following paragraphs.

Palpable breast lesions larger than 1 cm³ were only considered in the current implementation of the Wavelia QIF. All the objects with a volumetric size inferior to 1 cm³ were excluded and not morphologically validated with the Wavelia QIF. This setting was fixed in accordance with the expected minimum size of the detectable lesion with the first prototype of Wavelia and in order to avoid extreme degradation of the overall specificity of the system.

Detection of a single abnormality zone, which remained persistently contrasted against the surrounding breast tissues, was targeted with the applied MBI imaging algorithm and morphological lesion detector, at this preliminary stage of development, and in the context of the FiH feasibility study of Wavelia. Patients with bilateral breast disease were excluded and the focus was on the detection of a single (the largest) cyst, in the case of patients with multiple cysts in their breasts. This constraint is planned to be loosened in future clinical investigations, while evolving towards more realistic and generalized subsequent phases of the clinical evaluation of the Wavelia QIF.

The MBI scans were performed at least two weeks following the biopsy, in the case of Group-1 (invasive carcinoma) and Group-3 (benign biopsied lesion) patients. The two-week time lapse was considered sufficient to allow healing of the biopsy site in the breast. It is noteworthy though that, in most of these patients, a metallic biopsy clip was placed in their breast, as standard-of-care practice, to mark the tumor site. The size of the biopsy clip was small (~3 mm) compared to the targeted tumor sizes in this study (all palpable lumps), therefore its impact on the MBI images and the associated breast lesion detectability and characterization results were not considered to be significant. However, as the impact of the presence of a biopsy clip has not been characterized so far, future clinical investigations enabling the MBI examination to be performed prior to biopsy will be needed to investigate the impact of the biopsy clip. In this FiH study, there were only two Group-1 patients (P008, P043—both ILC's) and one Group-3 patient (P017—fibroadenoma) with no biopsy clip placed in their breast prior to the MBI examination and positive MBI findings nevertheless.

Quantitative evaluation, by means of a computable confidence level, is also planned to be implemented in subsequent versions of the Wavelia QIF, such that the imaging system performance can be assessed both in terms of lesion detectability rate and detection confidence level, in the case of various lesion types (solid/liquid, mass-like/non-mass-like, malignant/benign) and different breast density levels. Other factors, such as the breast size, the location of the lesion in the breast (superficial/deep, distance to the chest wall), the size of the lesion, the patient's age and breast deformability in the scanner, will also be investigated in terms of achievable lesion detectability rate and detection confidence level, as well as malignant-to-benign lesion separability. The assessment of these factors will be only feasible at a pivotal clinical investigation stage, and after sufficient stabilization of the Wavelia MBI scanning system and the associated QIF.

## 5. Patents

Two patents have been filed resulting from the work reported in this manuscript. The first patent covers the morphological breast lesion detection method, based on persistence, including sectorization of the imaging scene and MBI image reconstruction using multiple