Thigh muscle segmentation using a hybrid FRFCM-based multi-atlas method and morphology-based interpolation algorithm

The volume of lower extremity muscles is affected by some diseases. Quantiﬁcation of thigh muscles in medical images can lead to an easier investigation of these diseases. Most of the previous works in thigh muscle segmentation are based on models and atlases that require manually segmented datasets in 3D. As manual segmentation of these muscles is a time-consuming task, in this work, only one initial slice is segmented by a new hybrid FRFCM-based multi-atlas method and other slices are segmented based on this slice. In the proposed method, after noise reduction, the muscle region is extracted from other tissues by the FRFCM method. Then, an initial slice of each dataset is segmented by a multi-atlas method. The segmented muscles in the initial slice are used to segment muscles in the other slices of each dataset. The proposed method was evaluated with 20 CT datasets. The average DSC, Precision, and Sensitivity of the method for individual muscle segmentation were 91 . 20 ± 2 . 37, 91 . 95 ± 3 . 54, and 90 . 71 ± 3 . 89, respectively. The quantitative and intuitive results of the proposed method show the better results of this method in comparison to other state-of-the-art thigh muscle segmentation techniques.


INTRODUCTION
Thigh muscles quantification has an important role in surgical planning, postoperative assessment in orthopedic surgery [1], and diagnosis of some diseases such as muscular dystrophies [2], sarcopenia [3], and chronic obstructive pulmonary [4]. To quantify the volume of these muscles from medical images, the first step is segmentation.
Uemura et al. have assessed thigh muscles volumes changes after total hip Arthroplasty in CT images. They segmented the muscles manually and in a few slices [5]. The manual segmentation of thigh muscles in 3D medical images is an extremely time-consuming and tedious process. Therefore, the automatic segmentation of these muscles is very important and several studies have been devoted to the automatic segmentation of thigh muscles in 3D. Some of the previous works have segmented thigh muscles as one object [6][7][8][9][10][11] and some of them segmented muscles individually [1,[12][13][14][15][16]. Most of the previous works have been segmented thigh muscles in MR images, and only Yokota et.al have segmented the pelvic and hip muscles in the CT images [1].

MRI provides higher contrast in comparison
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. IET Image Processing published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology to CT for the soft tissues but it has a limited field of view due to its long acquisition time and requirement fine-tuned protocol design for each study [1]. Also, thigh muscle segmentation is important in postoperative assessment in orthopaedic surgery and CT images are commonly used in this surgery. Therefore, in this research, thigh muscles are segmented individually in CT-scan data sets.
Andrews et al. have segmented thigh muscles in MRI images using the generalized log-ratio representation in 3D [13]. This method is fully automated but the segmentation of the muscles is started just above the knee and extended 200 mm to the upper thigh. Kemnitz et al. [15] have segmented thigh muscles using active shape and active contour models in MR images. They segmented muscles between 25-65% of the thigh in four groups and not individually.
Yokota et al. have segmented thigh muscles using a hierarchical multi-atlas method in CT images. The accuracy of this method for the large muscle segmentation is high, but its performance is lower for the smaller muscles. Also, this method requires segmented muscles in 3D to create atlases, and manually segmentation of these muscles is a very time-consuming The block diagram of the proposed method process. In this research, only an initial slice is segmented using an FRFCM (Fast and Robust Fuzzy C-Means Clustering)-based multi-atlas method and the other slices of each dataset are segmented automatically using the registration and refinement of the initial slice. Therefore, this method has a lower computation than previous methods [1,13,15]. The main steps of the proposed method are as follows: 1. After noise reduction, the muscle region is extracted by the FRFCM method. 2. An initial slice of each dataset, located at 35% of the distal to proximal femur length, is selected. 3. All of the thigh muscles in the initial slice of each dataset are segmented using a multi-atlas method. 4. Some key slices are selected and muscles are segmented in these slices using a non-rigid registration and level set method.

Segmentation of the muscles is performed based on
Morphology-Based interpolation in the slices between the key slices.
The rest of this paper laid out as follows. In the "Materials and Methods" section, the proposed method is described. Afterward, in the "Results and discussion" section, the experimental results of the method are reported. The conclusion is presented in the final section.

Dataset
In this research, 20 CT datasets taken from [17] were used. There were 8 males and 12 females with a mean age of 50.09 ± 12.44 years. In these datasets, the image size and the slice thickness were 512 × 512 and 1.25 mm, respectively. The X-ray tube voltage was 110 or 130 kVp and tube current was between 80-136 mA.

The proposed method
In the proposed method, after noise reduction, the muscle region is extracted from the other tissues. Then, an initial slice is segmented by a multi-atlas method and this slice is used to segment the other slices in each dataset. Figure 1 shows the block diagram of the proposed method and its steps are explained as follows.

2.2.1
Thigh muscle extraction from other tissues At first, in order to reduce the noise of the image and to improve the edges, the BM3D-SH3D filter proposed by Dabov et al. in [18] is used. In this method, block matching of the image is conducted and a 3D matrix is constructed. A 3D transform is applied to this matrix and the noise of the image is attenuated by hard-thresholding the 3D transform spectrum. To enhance denoised coefficients, alpha-rooting on the hard-thresholded 3D transform spectrum is applied and the 3D transform is inverted to produce filtered grouped blocks.
Alpha-rooting image enhancement is realized through the following function: where t is a transform spectrum of a signal and contains a DC coefficient termedt (0), t sh is the spectrum of the resultant signal. An exponent > 1 is used to realize coefficient amplification and enhance the image details. After image enhancement, the thigh region in each axial slice is extracted by thresholding, and skin which has a similar intensity to the muscles is removed by morphological erosion. Finally, to extract muscles region from other tissues, as there are four main classes of intensities in CT images: "dark" (background), "dark-grey" (fat), "light-grey" (muscle and skin), and "light" (bone), the FRFCM algorithm is applied to classify these four classes. This method is an improved FCM, that is significantly faster and more robust than the FCM method. Although FCM is efficient for images with simple texture, it fails to segment images with complex texture or images corrupted by noise because it only considers grey-level information without considering the spatial information. Also, the number of iterations of FCM is larger for an image corrupted by noise than the image uncorrupted by noise. In the FRFCM method, local spatial information is used with a low computational and pixels' membership modified without depending on the calculation of the distance between pixels complexity within local spatial neighbors and clustering centres. The objective function of this algorithm is defined as [19]: where u kl represents the fuzzy membership of grey value l with respect to cluster k, and m denotes the fuzzy weighted factor. l is the number of pixels with a grey value equal to l . is the image reconstructed by the bi-directional morphological filter and l is a gray level, 1 ≤ l ≤ q, q denotes the number of the gray levels contained in . v k is the centre of the kth cluster and c defines the number of clusters.
Lagrange multiplier technique is employed to minimize the objective function and u kl and v k are obtained as follows: The optimal solution of u kl is obtained iteratively. Equations (3)-(4) are repeatedly implemented until their variations will be lower than a threshold.
After clustering the image using the FRFCM method, the second class with a higher intensity is considered as the muscle region. In order to find the muscles' weak boundaries, other tissues are replaced by the mean of the fat tissue class intensity. This approach has been employed due to the intermuscular adipose tissue between muscles. The obtained image is called a "muscles mask".

Individual muscle segmentation
Muscle segmentation individually consists of three steps. At the first step, an initial slice located at 35% of the distal to proximal femur length is selected. In this slice, all of the thigh muscles are exist. To segment muscles in this slice, a multi-atlas method is performed. In the second step, segmented muscles in this initial slice are used to segment muscles in the key slices. The final step is the segmentation of the muscles in the slices between key slices by an interpolation method. These steps are explained as follows: A. Initial slice segmentation using a multi-atlas method based on FRFCM algorithm: After noise reduction and muscle region extraction by the FRFCM method, the initial slice located at 35% of the distal to proximal femur length is selected in all 20 datasets. A twostage multi-atlas method is performed to segment muscles individually in the initial slice of each dataset. In order to perform the multi-atlas method for each dataset, the other 19 datasets are used as the atlases and the muscles are segmented individually in this target image. In the multi-atlas segmentation method, if there are n atlases,{A i } n i = 1 , each atlas is registered to the target image by a transferred function {T i } n i = 1 . Then, the labels of the atlases are transferred by the corresponding function. In the target image, for each voxel v, there are n label according to the atlases labels, If there are L labels in each atlas, l = {1, 2, … , L}, to find the final label of each voxel v, the probability function of the voxel is calculated as: After calculating the probability of the labels in each voxel of the target image, the label with the maximum probability is selected for the voxel.
Muscles masks extracted by the FRFCM method in the previous section are used for multi-atlas segmentation. To this aim, the following steps are performed: 1. At first, the muscles mask of each atlas is registered to the muscles mask of the target case using affine registration, and the deformation field is computed. 2. This deformation field is applied to the labels of the muscles in each atlas. 3. Consequently, the deformed muscles masks are registered to the target muscles mask using a non-rigid registration method based on Markov random field. 4. The registered labels in step 2 are deformed using this nonrigid deformation field.
Finally, the deformed labels of all atlases are fused and the target image is segmented. Affine registration is implemented by IRTK [20] and non-rigid registration is performed by Drop [21].
B. Muscle segmentation in the key slices using non-rigid registration and level set method: Since the near-slices in each dataset are more similar to each other in comparison with the other atlases, the initial segmented slice in each dataset is considered as the atlas in order to segment the other slices. At first, the initial segmented slice is modified by the user to segment the other slices more accurately.
In our datasets, the space between the slices is low (1.25 mm). Therefore, the cross-sectional area of the muscles changes very slowly between the slices. At first, the slice located at 20 slices upper/lower than the initial slice, defined as a key slice, is segmented using a non-rigid registration method based on the Markov random field. Then to refine the segmentation results in this key slice, a level set method is performed. After that, the next slice located at 20 slices upper/lower than this slice, defined as a key slice, is segmented in the same manner. This process is repeated until the last slice located at the top/down of the thigh is segmented.
To refine the results of the registration in each key slice, the boundary of the segmented area for each muscle is defined as the initial contour of the same muscle, and the Chan-Vese active contour model (ACM) is used to find the final boundary of the muscle. If the contours of the muscles have overlapped area, this area would be allocated to the muscle that belongs to it before applying ACM.
Finally, the middle slices between the key slices are segmented by a Morphology-Based Interpolation method explained in the next section.
C. The middle slices segmentation based on Morphology-Based Interpolation: To segment slices between the key slices, a Morphology-Based interpolation method proposed by Albu et al. in [22] is applied. In this method, Shape alignment is necessary for handling cases where the planar position of anatomical structures shifts significantly between slices. In this alignment approach, the translation vector that corresponds to the minimum shape displacement achieving maximum shape overlap is computed. After Shape alignment, a shape-based interpolation is done by conditional dilation operations as follows: Let A and B be two discrete sets in the same image lattice, such that B ⊂ A. The conditional dilation operation on B with respect to A using structuring element K is defined as: In order to generate A from B, m iterative conditional dilations with respect to A is required as: In this method, two parallel deformation processes based on iterative conditional dilation are used to create a transition sequence between two segmented key slices. The first process transforms the intersection A∩B into region A using lA iterative conditional dilations, with A as reference region and K as structuring element. Similarly, the second process transforms the intersection A∩B into region B using lB iterative conditional dilations. The mathematical description of these two processes is given by: The transition sequence using these processes is created as follows: (10) The integration of two dilations in order to transform shape A into shape B is described in the first branch of Equation (10). The second and third branches of this equation represent mutually exclusive alternatives to handle the difference in the lA, lB lengths of the two dilation processes that are combined.
To reduce the consumed computation time of the proposed method, all of the muscles are interpolated concurrently in the slices between key slices.

RESULTS AND DISCUSSION
In this research, to evaluate the performance of the proposed method, 20 CT datasets were used. The proposed method was implemented on MATLAB 9.1 environment and the ground truth segmentation was done on ITK_SNAP 3.6.0 [23]. Slicer 4.10.2 (www.slicer.org) [24] was applied to visualize the segmentation results.  To evaluate the method, the Dice Similarity Coefficient (DSC), Precision, Accuracy, Sensitivity, and Specificity [25] were calculated between the manual segmentation and the result of the presented method. These metrics are computed according to the following equations: Speci ficity = TN FP + TN (15) In equation (11), A and B are the segmentation results of the proposed method and ground truth. TP, FP, TN, and FN are True Positive, False Positive, True Negative, and False Negative, respectively.
At the first stage of the method, thigh muscles are segmented using the proposed multi-atlas method based on the FRFCM algorithm in the initial slice. The Dice similarity coefficient (DSC) of this method in comparison to the hierarchical multi-atlas method proposed in [1] and U-Net for all muscles are reported in table 1.
According to table 1, the proposed multi-atlas method has better results than the previous hierarchical multi-atlas method and U-Net in the initial slice. As the deep convolutional neural networks need more data to have the desired performance, it is expected that the results of the U-Net will improve with increasing the number of data sets. However, with the limited numbers of data sets, the proposed method is more accurate. The intuitive results of the proposed method for the initial slice of some cases are shown in Figure 2.
As shown in Figure 2 the multi-atlas method has better results for the cases with more adipose tissue and for the cases that aren't very different from atlases. Consequently, the segmentation precision for muscular cases is lower than the obese cases. According to Figure 2 and Table 1, the accuracy of the proposed multi-atlas method in the segmentation of the Rectus femoris, Semimembranosus, and Semitendinosus muscles is lower than the other muscles. It is due to this fact that the Semitendinosus and Semimembranosus muscles are adjacent and due to the partial volume effect, their boundary is invisible in the muscular cases. Also, the boundary of the Rectus femoris muscle with the other quadriceps muscles is indiscernible in the muscular cases, which reduces the accuracy of the segmentation in the very different cases from atlases. As the initial slice of each dataset is used to segment other slices, the results of the multiatlas method are improved by the user. Then, the other slices are segmented with the proposed method. The initial slice is selected in the part of the thigh, where all the muscles present. It was heuristically found that all of the thigh muscles segmented in this research exist in the slices approximately located in a distance between 25% to 45% of the distal to proximal femur length. This region is between the origin of the short head of the Biceps femoris muscle and the upper end of the Quadriceps tendon. Therefore, the middle slice in this region is selected. However, other slices in this distance could be selected as the initial slice. For example, the Dice similarity coefficient of the method with the initial slice located at 30%, 35%, and 40% of the femur length are reported in Table 2. These results show the method is not very sensitive to the initial slice selection.
In this research, the distance between the key slices is chosen 20 slices. The overlap between the key slices increases if the distance between the slices reduces. Therefore, the results of the initial slice registration on the key slices and the final interpolation between the key slices would be better, but the time consumption will increase. Table 3 shows the mean of the Dice similarity coefficient and consumed computation time for the method with 10, 20, and 30 slices distance between the key slices. According to these results, by choosing a distance of 20 slices between the key slices, the performance of the algorithm does not decrease very much. However, the implementation time is reduced by 25%. The boxplots of the Dice, Precision, and Sensitivity of the proposed method for 3D segmentation of thigh muscles are shown in Figure 3(a). Figure 3(b) shows the boxplots of the specificity and accuracy.
The overall means of the DSC, Precision, and Sensitivity for all individual muscles are 91.20 ± 2.37, 91.95 ± 3.54, and 90.71 ± 3.89, respectively. The means of the Specificity and Accuracy of the method are more than 99%, which is due to a very smaller area of each muscle compared to the image size which results in a high True Negative. The intuitive results of the proposed method and manual segmentation of two cases in 3D are shown in Figure 4 from the right and left views.
According to Figures 3 and 4, the proposed method has very good performance in thigh muscle segmentation. The mean of the Dice similarity coefficient, sensitivity, specificity, accuracy, and precision of all individual muscles for some typical cases are shown in Figure 5.
Since in the proposed method only one initial slice is segmented by the multi-atlas method and the other slices are segmented based on this slice, the average consumed computation time for the proposed method is about 3 min per case which is very lower than the previous work [1] (41 min). The mean of the consumed computation time by the several steps of the method are reported in Table 4. Also, the average Dice similarity coefficient of the method is significantly higher than the previous work [1] (91.20 ± 2.37 vs 83.76 ± 5.00). Furthermore, the mean of the DSC for the proposed method is higher than the U-Net (91.20 ± 2.37 vs 85.27±5.66).
The consumed computation time of the proposed method in comparison with the previous works in thigh muscle segmentation from MRI datasets [13,15] is very low. Also, the segmented volume of the thigh in the proposed method is larger than these works. The average Dice similarity coefficient of the method is higher than the method introduced in [13]. Although the method proposed in [15] has a higher Dice similarity coefficient than the proposed method, this method has segmented muscles in four groups and not individually.  The results of the proposed method (top) and manual segmentation (down) for two cases in 3D from right and left views. The differences between the results of the method and ground truth are shown by white arrows

CONCLUSION
Thigh muscle quantification is very critical in surgical planning and postoperative assessment in the orthopedic surgery. Segmentation of these muscles from one another is a challenging task. As CT images are commonly used in orthopedic surgery, in this research, thigh muscles were segmented in these images. To this aim, a new hybrid method based on the multi-atlas method, FRFCM algorithm, and level set method has been proposed.
In this method, for all 20 datasets, an initial slice located at 35% of the distal to proximal femur length is segmented by an FRFCM-based multi-atlas method. As the variations of the muscles volume along with their length are slowly and continuously, the muscles could be segmented in the slices near to the initial slice by a non-rigid registration of the initial slice to these slices. In the proposed method, the slice located at 20 slices upper∖ lower, key slices, than the initial slice has been segmented using non-rigid registration and the initial slice is chosen as the The average results of the proposed method for all individual muscles of some typical cases atlas to segment these key-slices. Afterward, the segmentation of the muscles has been improved using the level set method in this slice and the final segmentation has been registered to the slice located 20 slices upper∖lower than this slice. The process is repeated up to the insertion of the Gluteus Maximus muscle at the top of the thigh and to the knee in the downward direction. Finally, after segmenting the key slices, the middle slices are segmented by a Morphology-based interpolation method. The quantitative and intuitive results of the proposed method show higher accuracy of the method than the previous work. Also, this method is very accurate for the small and large muscle segmentation, whereas the previous work had less accuracy in the small muscle segmentation.
The goal of this study is the visualization of muscle fibers. To this end, the next step of this research is to estimate the muscle attachment regions to the bones which are necessary for fiber arrangement estimation. Therefore, future work is devoted to estimating bone-muscle attachment in CT images.