c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

Computer-based assessment for facioscapulohumeral dystrophy diagnosis c ˇ O. Chambers a,∗ , J. Milenkovi´c b,d , A. Praznikar , J.F. Tasiˇc b a

ˇ Stefan”, Jamova cesta 39, 1000 Ljubljana, Slovenia Institute “Jozef ˇ ska cesta 25, 1000 Ljubljana, Slovenia University of Ljubljana, Faculty of Electrical Engineering, Trzaˇ c University Medical Centre of Ljubljana, Department of Neurology, Zaloˇska cesta 2, 1000 Ljubljana, Slovenia d Faculty of Medicine, Vrazov ˇ trg 2, 1000 Ljubljana,Slovenia b

a r t i c l e

i n f o

a b s t r a c t

Article history:

The paper presents a computer-based assessment for facioscapulohumeral dystrophy

Received 25 November 2014

(FSHD) diagnosis through characterisation of the fat and oedema percentages in the muscle

Received in revised form

region. A novel multi-slice method for the muscle-region segmentation in the T1-weighted

27 February 2015

magnetic resonance images is proposed using principles of the live-wire technique to find

Accepted 23 March 2015

the path representing the muscle-region border. For this purpose, an exponential cost function is used that incorporates the edge information obtained after applying the edge-

Keywords:

enhancement algorithm formerly designed for the fingerprint enhancement. The difference

MRI

between the automatic segmentation and manual segmentation performed by a medical

Hyperintensity quantification

specialists is characterised using the Zijdenbos similarity index, indicating a high accuracy

Muscle dystrophy

of the proposed method. Finally, the fat and oedema are quantified from the muscle region

Muscle-region segmentation

in the T1-weighted and T2-STIR magnetic resonance images, respectively, using the fuzzy

FSHD

c-mean clustering approach for 10 FSHD patients. © 2015 Elsevier Ireland Ltd. All rights reserved.

1.

Introduction

Facioscapulohumeral dystrophy (FSHD) is the third most common muscular dystrophy after Duchenne’s and myotonic dystrophy [1]. FSHD is characterised by progressive degenerative changes in the muscle fibres resulting in accumulation of the IntraMuscular Fat (IMF) [2,3]. IMF can be defined as a sum of the infiltrated fat within the muscle region and intermuscular fat between individual muscle groups [4]. Research studies indicated a correlation between the muscle force and level of IMF in the muscle region of the FSHD patients [5,6], which could provide an important prediction for the disease progression [7]. The T1-weighted magnetic resonance imaging

∗

Corresponding author. Tel.: +386 1 477 39 00; fax: +386 1 251 93 85. E-mail address: [email protected] (O. Chambers).

http://dx.doi.org/10.1016/j.cmpb.2015.03.006 0169-2607/© 2015 Elsevier Ireland Ltd. All rights reserved.

(MRI) serves as a non-invasive biomarker of the IMF distribution in the lower limb enabling therapy-treatment estimation [8]. Interpretation of the MRI images to determine whether the FSHD pathology is present or not is relatively easy, however, it is much more difficult to quantify the muscle changes both in terms of the affected muscle-volume percentage and disease progression. Recently, the increased water mass, or oedema, has been demonstrated as a crucial intermediate state of the change between the healthy muscle and fatty replacement [9,10]. The oedema areas characterised by an increased signal in the T2-weighted short-tau inversion recovery (T2-STIR) images in muscles not yet replaced by a fat tissue may show a normal signal in the T1-weighted images [11]. Therefore, this

38

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

complementary source of information is necessary to estimate the FSHD progression and to make therapeutic decisions. Manual assessment of the FSHD progression in the MRI images is a relatively complex and subjective task for large datasets. Therefore, automatic and semi-automatic methods of segmentation and quantification for diagnostic assessment are preferable in clinical evaluation. In the literature, the muscle-region segmentation is commonly performed on the fat-unaffected muscle structures, where a good contrast between the subcutaneous adipose tissue (SAT) and muscle tissues is observed. In those cases, the segmentation is usually performed by using a simple histogram analysis and thresholding operations [12]. More advanced techniques for the muscle-region segmentation incorporate both the spatial and intensity information of pixels. The fuzzy c-mean (FCM) clustering [13,14] and level-set approach [15–18] are the most commonly used techniques [19,17]. However, the presence of IMF in the muscle region of patients affected by FSHD makes the segmentation more difficult, because of the IMF and SAT gray-level intensity similarity. Moreover, the spatial resolution limitations blur the border between SAT and IMF making it almost invisible. Wald et al. [14] proposed a morphological opening with a circular structure element to break any connection between SAT and IMF after using the FCM clustering. However, this may lead to a loss of the muscle-region border information. By using the level-set approach, Positano et al. in [17] observe intrusion of the internal SAT border into the intermuscular fat. Moreover, the blood-vessel regions may significantly affect the effectiveness of the level-set approach. A more recent work of Commean et al. in [20] proposes the use of the second derivative of brightness in the fat-saturated MRI images for semi-automatic delineation of the lower-limb anatomical elements, where the user is required to break connections in the binary edge map between individual elements. This method is sensitive to the texture edges and IMF in muscles. The authors reported that more advanced edge-detection techniques may improve results. In this paper, an edge-enhancement algorithm designed for fingerprint enhancement is used for segmentation showing advantages over the classical edge-detection methods for the lower-limb segmentation. The proposed method uses the basic principles of the live-wire technique introduced by Mortensen and Barrett [21], where the manual interaction is substituted with an automatic seed-point placing. Originally, this technique searches the path with the smallest cost between the manually placed seed points, where the cost is calculated by incorporating the information from the edges obtained by the Laplace zero-cross filtering, gradient magnitude and gradient direction. The multi-slice muscle-region segmentation is performed in two phases. During the first phase, the muscle region of the lower-limb with the largest cross-section area is segmented and used to perform segmentation in the adjacent slices during the second phase. The obtained results using the proposed method are compared with the results of the FCM clustering and level-set approach with respect to the medical specialists’ annotations using two metrics, the root mean square and Zijdenbos similarity index. Finally, the fat and oedema percentages in the muscle region are calculated to characterise the FSHD progression.

2.

Materials

2.1.

Magnetic resonance imaging parameters

Our MRI examination of 10 patients with a genetically confirmed diagnosis of FSHD was performed on a 3.0T Philips scanner using the software version 3.2.1. Imaging was made at the calf level with a 5-mm thickness parameter. The T1-weighted images were acquired using a Fast SpinEcho (FSE) sequence with the following parameters: TE/TR = 20/500 ms, ETL = 3 (Echo Train Length), NEX = 1, 784 × 784 acquisition matrix size, FOV = 330 mm. The T2-STIR images were acquired using an FSE sequence with the following parameters: TE/TR = 30/1400 ms, ETL = 15, NEX = 2, 1152 × 1152 acquisition matrix size, FOV = 330 mm and 5 mm slice thickness.

2.2.

Dataset description

The T1-weighted and T2-STIR MRI images contain a different but complimentary information necessary for a computerbased diagnostic assessment of FSHD. The T1-weighted MRI images of the lower limbs provide information required for the muscle-region segmentation and fat quantification, while the T2-STIR images provides information required for the oedema quantification [22]. Fig. 1 shows a typical example of the T1weighted and T2-STIR images of a patient with FSHD. Anatomically, the muscle region is surrounded by SAT. In the FSHD patients, the muscle region contains also IMF [23]. Each bone, e.g. fibula and tibia, is constructed of a bone marrow (higher signal intensity) and surrounding cortical bone (lower signal intensity). In the T1-weighted MRI images the muscle tissue appears as a medium grayscale intensity level, while SAT surrounding the muscle region and IMF appear with a high grayscale intensity level similar to that of the bone marrow. The intensity inhomogeneity may be observed in the T1weighted images as signal-intensity variations. T2-STIR has a higher signal-to-noise ratio, however, it is independent of the magnetic-field inhomogeneity.

3.

Methods

The muscle-region segmentation is performed in the T1weighted MRI images incorporating information from the adjacent slices. After segmentation, quantification is performed to estimate the fat/oedema percentage in patients with FSHD.

3.1.

Muscle-region segmentation

The proposed segmentation method relies on the anatomy similarity between the adjacent slices. The segmentation is performed in two phases illustrated in Fig. 2. It can be noticed that both phases have a similar flow. The principal difference between them refers to the input parameters for the shortestpath search, where the shortest-path search technique is used to find the path along the SAT internal border.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

39

Fig. 1 – Lower-limb MRI images of a patient affected by FSHD. Top image: T1-weighted MRI image with fat infiltration; bottom image: T2-STIR image with oedema.

An edge-enhancement algorithm resulting in the Edge Map (EM) is used for the path-cost calculation in the initial slice. For other slices, EM is optimised by incorporating the border information from the adjacent slice resulting in an optimised EM or OEM. The optimisation is performed during an edge-enhancement procedure by redefining the edge flow with respect to the so-called border mask. The border mask is obtained as a tube-like area around a border copied from the nearest adjacent slice using a morphological dilation by a circular structuring element of a small diameter. The border mask is used to specify the relevant edges to reduce the overand miss-segmentation during the shortest-path search. The initial slice is identified as a slice with the largest cross-sectional area by summing the pixels belonging to the limb-region mask in each slice. Then, the fibula and tibia bones are excluded from the region limited by the calculated shortest path to obtain the muscle region. The multi-slice segmentation flow continues in both directions, interiorly and superiorly, from the initial slice until all the slices are segmented.

3.1.1. Fig. 2 – Scheme of the multi-slice muscle-region segmentation in the lower limb using the proposed method.

Edge map

The edge map is obtained by using the Hong et al. [24] edge-enhancement algorithm formerly designed for fingerprint enhancement in the fingerprint matching application [25] and being based on the tenable Gabor filter.

40

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

The edge-enhancement algorithm includes the following steps: Step 1: The orientation field is obtained by dividing the input image into square blocks, where the gradient is calculated for every pixel in each block. The orientation vector for each block is derived by averaging all the vectors orthogonal to the gradient pixels in the block. The orientation image is then smoothed using a low-pass filter. Step 2: The edge frequency image is obtained by dividing the input image into square blocks, where the oriented window [24] is calculated for each block. Then, the projection of all the gray-level values in the oriented window along a direction orthogonal to the orientation vector is performed. This projection forms a sinusoidal-shape wave and the distance between the consecutive peaks in this wave is used to estimate the frequency of the edge. Step 3: Based on the local orientation and edge frequency around each pixel, the Gabor filter is applied to each pixel location in the image. Their values are used to construct an even-symmetric Gabor filter. Gabor filter G is applied to the input image I by spatially convolving the image with the filter to obtain enhanced image E using Eq. (1), where O is the orientation field and F is the edge frequency image.

E (i, j) =



The output of the edge-enhancement algorithm is an enhanced version of the image edges that have been smoothed in the direction of these edges. By applying the global threshold of zero on image E, the binary image is obtained. The morphological thinning [26] of the obtained binary-image supplementation produces EM, where the thinning erodes the binary-object pixels until they are one pixel wide. The result of the edge-enhancement algorithm can be regulated by changing the algorithm-step parameters. It was found that the algorithm produces more accurate edges along the required structure if the orientation field is calculated directly from the mask, termed here as border mask, created as a tubelike area around the border by using morphological dilation. The OEM is obtained from Steps 1 to 3, where the orientation field in Step 1 is calculated directly from the border-mask image optimising the edges only within the border mask. Fig. 3 shows the impact of the orientation field on the edgeenhancement algorithm performance. It can be observed that OEM is more smooth and more connected, while EM contains more disconnected and less smooth edges within the border mask.

3.1.2. G (u, v; O (i, j) , F(i, j)) I (i − u, j − v)

(1)

Shortest-path search

The shortest path is calculated between the seed points along the SAT internal border, where the seed points are used to

Fig. 3 – Results of the edge-enhancement algorithm obtained by using an oriented image calculated from the original MRI (first row) and from the border-mask image (second row).

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

41

Fig. 4 – Example of the detected seed points (marked as blue dots). Left image: seed points in the initial slice (Phase I); right image: seed points in the adjacent slice (Phase II). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

regulate the path direction. The number and position of the seed points are empirically selected and tested to obtain a relative combination that gives the best overall result. Once the seed points are placed, the shortest-path search procedure starts, where the image pixels are associated with the graph nodes and connection of two 8-adjacent image pixels with the weighted graph edges. To perform this search, one of the seed points is selected as the starting point and then Dijkstra’s algorithm [27] is used to search for the paths in the clockwise direction until reaching the starting point. A local cost is assigned to each graph edge to weight their probability of being included in a path. The path with a minimal sum of costs is defined as the shortest path.

3.1.2.1. Seed points. The seed point placing was performed on EM along the SAT internal border as it has a common border with the muscle region. Anatomically, SAT circulates around the lower limb and, therefore, its internal border may be associated with a circle that has the centre in the limb-region centre of the mass. To perform the seed-point placing in the initial slice, the circle is divided into 16 equal segments, where the circle radius is selected to include the whole limb region. By using the lines between the segments endpoints and the centre of the circle, the seed points are found on the intersection of the EM subedges (edge segments between the points of bifurcation) with these lines. Firstly, the largest binary element of the muscle tissue mask is identified using the second mask of the three-class FCM and then the most outward mask point from the circle centre that has the intersection with the lines is identified as the starting seed point. The second seed point is identified right from the starting point on the intersection of the nearest line with a sub-edge that belongs to the same edge. If more than one intersection point is identified, then the point closest to the starting point is selected. If all the intersection points

belong to the disconnected sub-edges, then the seed point is placed on the intersection with the sub-edge which beginning deviates no more than the pre-defined error measure from the end of the sub-edge corresponding to the previously found seed point. If the sub-edge with the aforementioned rule is not identified, then no seed point is placed in that segment. The error measure needs to be identified by the user. Other seed points are identified similarly as the second point with respect to the previously found seed point. For segmentation on the adjacent slices, only the relevant edges within the border mask are used and, therefore, less seed points are required for the shortest-path search. Four of the 16 uniformly distributed lines are used to find four seed points on the intersection with the edges within the border mask. Fig. 4 shows the seed points placed in the initial slice and adjacent slice for one example from the research dataset.

3.1.2.2. Cost function. Since it is known that the SAT interior border is a continuous circular border, the path should be enforced to have a non-linear form. This is especially critical for the regions with no edge score. The exponential cost function defined in Eq. (2) achieves this and is found to give a robust result, where cost(u) is the cost of the graph edge in node u. cost(u) = exp(1 − EM(u))

(2)

During Phase II, the shortest path is found by incorporating only the relevant edges of OEM (see Fig. 4) using Eq. (3), where cost(u) is the cost of the graph edge in node u and B is the border mask of the nearest adjacent slice. The cost function defined in this way ensures a more accurate and robust result. cost(u) = exp(1 − OEM(u) · B)

(3)

42

3.1.3.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

Muscle-region border detection

The muscle-region border is obtained by subtracting the tibia and fibula bones from the region limited by the obtained shortest path along the SAT internal border. This subtraction is performed by using morphological operations. In the T1-weighted MRI images, the lower-limb intensity may be classified into three groups: low-signal intensity (background, vessels, cortical bone, etc.), medium-signal intensity (muscle tissue) and high-signal intensity (subcutaneous adipose tissue, muscle fat, bone marrow, etc.) Therefore, the three-class FCM clustering is applied on the intensity-inhomogeneity-corrected image to produce three binary masks corresponding to the low-, medium- and highsignal intensity, respectively. The cortical bones are identified as two binary elements with the largest area within the region limited by the shortest path along the SAT internal border in the first mask and then filled to include all the pixels corresponding to the bones. To guarantee the smoothness of the muscle-region borders, the morphological erosion and dilation (four pixels) are applied on the muscle-region mask.

3.1.3.1. Quantification. Once the muscle region is obtained, the quantification of the fat/oedema is performed, where by the fat quantification we mean the quantification of IMF. Fat quantification. The IMF quantification relies on the large difference in signals between the fat and muscle tissues in the T1-weighted images. To correct the intensity inhomogeneity in the T1-weighted images, the method proposed in [28] is used. Moreover, the four-class FCM clustering used in [28] to quantify the muscle fat in the 1.5 T MRI images is used here for the fat quantification in the 3.0 T MRI images, where an average between the maximum signal intensity in the second class and the minimum signal intensity in the fourth class is used. The clustering is applied using all the image pixels and the fat is recognised as a binary mask within the pre-defined muscle region. Oedema quantification. As the T2-STIR images have different matrix sizes, they are, firstly, rescaled on the dimensionality of the T1-weighted images to have the pixel-by-pixel correspondence between the images and, then, the iterative affine registration [29] is applied to reduce the motion artefact. Thus, the muscle-region border obtained in the T1-weighted image may be directly projected with a high accuracy on the corresponding T2-STIR image of the same patient. The muscle oedema is characterised by the presence of the muscle hyperintensity in the T2-STIR MRI images [30,31] that is almost always due to the increased intracellular or extracellular free water characterising biologically active processes. There are no methods proposed in the literature for the lowerlimb oedema quantification using the image processing tools. Therefore, we applied a technique similar to that one used for the quantification of the myocardial oedema [32]. Through experimental analysis it was found that the minimal image intensity within the fourth class (the class with the highest average intensity value) of the four-class FCM is the most optimal threshold for the oedema intensity estimation in the T2-STIR images that was verified by a medical specialists. Fig. 5 shows an example of the fat and oedema quantification results from the T1-weighted and T2-STIR images,

respectively, of the same patient. The quantification results are shown as binary masks below the original images. It can be observed that the border between SAT and IMF is easily observed in the T1-weighted images, while for the T2-STIR images this border is partially invisible. Thus, a direct segmentation of the muscle tissue from the T2-STIR images is difficult indicating the importance of the T1-weighted images as a supplemental information.

4.

Validation metrics

To validate the accuracy of the muscle-region segmentation and fat/oedema quantification, the obtained results are compared to the manual segmentation. In order to reduce the interobserver variability, two medical specialists (A.P. 14 years of experience and A.Z. 23 years of experience) were asked to perform manual segmentation independently and then their results were averaged. By using the image analysis software ImageJ version 1.34 [33], specifically utilising the ROI Manager, Threshold and Measure tools, the medical specialists were asked to segment the muscle region manually and then to find the best thresholding parameter for the fat/oedema quantification from the intensity-inhomogeneity-free images with the following manual corrections for each FSHD patient. Two types of corrections in the obtained fat/oedema quantification results were made: inclusion of the missing fat/oedema regions and exclusion of the regions that were wrongly included as part of the fat/oedema. Two metrics, the root mean squares (RMS) and Zijdenbos similarity measure (ZSI) [34], are chosen for the numerical validation. RMS measures the difference between the two curves by using Eq. (4), where d = {d1 , d2 , . . ., dN } is a set of the minimal Euclidian distances from the N equally spaced points on the shortest path to the medical annotation.

  N  RMS =  d2 i

(4)

i=1

ZSI describes the similarity between the muscle region obtained by using the proposed automatic method (A) and manual segmentation (M) using Eq. (5). The higher is the value of ZSI, better is the correspondence between the manual and automatic segmentations.

ZSI =

2 × |A ∩ M| |A| + |M|

(5)

The non-parametric Wilcoxon signed-rank test was applied to determine the significance level based on the difference in quantification results obtained by using automatic segmentation and manual interaction. The difference was considered to be statistically significant if the estimated pvalue of the statistical test was less than the risk level (␣=0.05).

43

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

Fig. 5 – Fat/oedema quantification from FSHD MRI images. First row: the axial T1-weighted (left) and T2-STIR (right) MRI images of the same patient affected by FSHD; second row: binary masks obtained by using the FCM clustering approach.

5.

Results

The proposed method was applied for the segmentation and quantification of the 10 T1-weighted FSHD MRI images. After the muscle-region segmentation, the fat and oedema percentages were calculated within the obtained muscle region in the T1-weighted and T2-STIR MRI images, respectively. The different stages of the proposed method were validated. Table 1 provides the muscle-region segmentation accuracy of the proposed method in terms of RMS and ZSI for each patient. Fig. 6 shows an example of the muscle-region segmentation results for a randomly selected FSHD MRI image from the research dataset. The three-dimensional surface model is generated by interpolating the adjacent muscle-region borders. The fat and oedema quantifications were performed for the left and right limbs separately, where the hyperintensity observed in the muscles in the T1-weighted MRI images was considered as IMF and the muscle hyperintensity observed in the T2-STIR MRI images was considered as an oedema. The obtained average accuracies of the fat and oedema quantifications for all patients in terms of ZSI are 0.89 and 0.83, respectively. There was no significant difference observed

between quantification results obtained by using automatic segmentation and manual interaction (p>0.05). For the medical diagnosis purpose, the fat percentage was calculated in the total muscle region, while the oedema percentage was calculated only within a fat-unaffected part of

Table 1 – An average root mean square (RMS) and Zijdenbos similarity index (ZSI) between the manualand automatic-segmentation results obtained by using the proposed method for the 10 FSHD MRI images from the research dataset, where the mean and standard deviations (std) are calculated for the patient MRI stack. Patient ID P1 P2 P3 P4 P5 P6 P7 P8 P9 P10

RMS (mean ± sd) 0.12 ± 0.05 0.08 ± 0.03 0.16 ± 0.03 0.12 ± 0.08 0.11 ± 0.05 0.06 ± 0.03 0.18 ± 0.04 0.02 ± 0.05 0.04 ± 0.02 0.13 ± 0.09

ZSI (mean ± sd) 0.91 ± 0.02 0.92 ± 0.01 0.88 ± 0.03 0.90 ± 0.02 0.91 ± 0.03 0.96 ± 0.02 0.90 ± 0.04 0.95 ± 0.03 0.96 ± 0.01 0.90 ± 0.05

44

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

Fig. 6 – MRI stack-segmentation results for patient P9 using the proposed method, where the resulting muscle region is marked by a red colour. The left image shows a 3-D volume reconstruction obtained by interpolating the obtained muscle regions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

the muscle region using Eq. (6), where the area of the object is the number of the pixels belonging to this object.

=

Area (hyperintensity) Area (muscle)

(6)

The quantification results are shown in Table 2 for each FSHD patient, where the result for each patient is presented as a fat/oedema percentage in the MRI stack. The patients gender, age and body mass index (BMI) are also included in the table. There are several things to be noted in Table 2. The obtained results indicate the oedema presence in the fat-unaffected muscles for all the FSHD patients that provides an important information for the muscle-dystrophy analysis in an early stage. The fat distribution asymmetry in the left and right FSHD lower legs is also observed confirming the previous reports [35,6]. In [31], the asymmetric muscle involvement is seen in 15% of the dataset including five patients, while in our dataset, the asymmetry is observed in 40%. Moreover, the asymmetry of the oedema presence is also observed, that may characterise further muscle degeneration. There is no correlation observed neither between FSHD disease progression and the patients’ age nor between FSHD disease progression and the patients’ gender. The youngest patients (P1 and P2 in Table 2) were still unaware of the disease symptoms at 38 and

35 years of age, while the oldest patient (P7 ) did not show the largest muscle degeneration. The muscular-dystrophy patients generally have a lower body weight and a smaller muscle volume compared to the control patients, because of the loss of the muscle mass [36]. Most of the patients in our study had the normal BMI, while having different stages of the muscle degeneration. Nevertheless, a correlation between BMI and muscle involvement can be observed. Thus, the patients with a higher BMI (i.e. P1 , P9 , P10 ) have a lower muscle involvement than the patients with a lower BMI (i.e. P3 , P5 ). The proposed method was implemented on a PC (Intel§Core (TM) i7 CPU 920 @ 2.67 GHz, Memory (RAM): 3 GB, Windows Vista 64-bit Operating System) with Matlab 2012b. EM was obtained by using Kovesi’s Matlab implementation [37] of the Hong et al. [24] edge-enhancement algorithm.

6.

Discussion

The study suggests the feasibility of using automated image processing for the muscle-region segmentation to quantitatively measure the muscle hyperintensity for the FSHD diagnosis. The experimental results show that the proposed method allows computation of a fairly stable and accurate set of the MRI biomarkers reported in Table 2.

Table 2 – The percentages (%) of the fat in the whole muscle region and percentages of the oedema in the fat-unaffected part of the muscle region in the left/right lower limbs of the patients affected by FSHD. Patient ID P1 P2 P3 P4 P5 P6 P7 P8 P9 P10

Gender F M M M F M M M F M

Age 38 35 52 45 51 50 60 51 49 44

BMI (kg/m2 ) 26.3 23.1 22.6 25.4 18.4 25.0 24.2 24.3 25.8 24.7

Left/right fat 4/3 4/3 75/83 22/8 59/31 37/44 39/35 52/51 11/6 7/20

Left/right oedema 4/9 11/13 16/15 16/8 20/24 17/21 12/19 16/7 5/8 14/4

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

The intermediate result of the muscle-region segmentation is the SAT internal-border detection that has a practical importance in clinical diagnosis. The performance of the SAT internal-border detection was compared by using the proposed method with the results obtained by using the most common literature techniques, such as the FCM clustering, level-set and classical implementation of the live-wire, implemented for the research dataset. Fig. 7 illustrates the results of the SAT internal-border detection in images with various fat percentages, where it can be observed that the proposed method produces a more smooth and accurate border compare to the literature methods. The boxplot below images visualise the average segmentation accuracy of these methods in terms of ZSI for research dataset indicating a better matching of result obtained using our method with the manual segmentation. Result of the FCM clustering shows the largest dispersion indicating less robust to the IMF presence. The poor results of the literature methods may be explained by their sensitivity to the fat infiltration in muscles on the border with SAT. For segmentation using the FCM clustering, three classes are used with a subsequent morphological opening to disconnect IMF and SAT as proposed in [14]. The initial contour for the level-set evaluation is obtained as a limb-region external contour as proposed in [17]. To obtain the optimal segmentation results with the level-set, the parameters such as the distance regularisation weight, contour integral weight, and pressure force are set for each image individually due to the anatomical variability. The live-wire technique is performed using the seed points found as proposed in this paper for the initial slice. This technique has already been recognised as a relevant semi-automatic segmentation method for the border detection in medical images. However, in classical implementation of live-wire, the definition of the cost function is shown to be ineffective for the research dataset. This is due to the poor results of the edge detection by the zero-cross operator presented in cost function. Most of the literature edge-based methods are sensitive to the additive noise and small textural information leading to an incorrect edge detection at the required border. The edge-enhancement algorithm used in the proposed method, formerly designed for the fingerprint application, shows to be more robust to small textural edges producing smooth edges along the SAT internal border in the FSHD MRI images. Moreover, the border mask defined in this paper using the information from the adjacent slice restricted the border search area leaving only the edges corresponding to the required structure that significantly speeds up the processing time and accuracy. Another important advantage of the proposed method is that it requires no training set or shape constraints learning, performing the segmentation for each patient independently. It has to be noted that a direct segmentation of the muscle region in the T2-STIR image is challenging due to the poor information of the musculoskeletal structure and low contrast between SAT and muscle tissue. Therefore, the T1-weighted MRI is important for the automatic T2-STIR muscle-region segmentation. This study contains some limitations that should be recognised. One limitation is that the proposed method is based on the assumption that the image properties at the SAT

45

border are similar between the adjacent slices. Therefore, the slice thickness should not be more than 5 mm to guarantee an accurate edge-map optimisation. A large slice thickness may lead to a non-correspondence of the adjacent slices anatomy and, therefore, inaccurate multi-slice segmentation. Recent research studies have extended the 2D segmentation methods into the 3D segmentation methods because 3D analysis provides much richer information. Thus, it is possible to improve the proposed segmentation method by incorporating third dimension (z-direction) in order to segment the whole MRI stack using the 3D methods as the 3D FCM clustering, 3D edge detection, etc. However, there is no trivial way to implement a 3D segmentation because the increase in the dimension increases the complexity of the seed-points placing and the exponential increase in computational cost [38]. Therefore, additional research is required to develop a 3D segmentation. The proposed multi-slice segmentation converts a 3D problem into 2D problem while using information from the adjacent slices, i.e. prior knowledge, for the segmentation task being capable of running in real-time for processing whole 3D MRI stack [39]. As there were only a limited number of patients available for analysis, the robustness against the differences in the population and acquisition parameters could not be investigated. Nevertheless, the results reported in Tables 1 and 2 indicate the similarity of the segmentation results to those obtained manually and robustness of the proposed method to various percentages of fat infiltration. For more robust validation of the proposed method, the research dataset needs to be enlarged with more FSHD MRI cases and cases of other types of muscular dystrophy. Once the muscle-region segmentation is performed, the fat/oedema is quantified from the MRI images to characterise the FSHD progression. The volumetric data are still not widely used for clinical diagnosis of disease progression. It is still unclear which measure of the muscle morphology, e.g. the muscle volume or cross-sectional area of a single-slice, is preferable. Some studies recommend the use of the muscle volume, whereas others recommend the cross-sectional area. In [40], it was determined that a cross-sectional area provides the same sensitivity to the fat change as a total muscle volume. In this study, rather than focusing on volumes, we computed the average fat percentage over all cross-sectional areas in MRI stack for each patient. The quantification accuracy may also depend on the quality of the image acquisition. Different MRI protocols affect the reproducibility and accuracy of any automatic algorithm. The T1-weighted MRI was recognised as a relevant technique for the fat quantification when using the FCM clustering [41,19]. In our previous study [28], the fat quantification by the four-class FCM was proposed for the 1.5 T MRI, where the average ZSI equal to 0.91 was obtained. This approach is also verified here for 3.0 T MRI images giving the average ZSI of 0.89. The similarity between the results indicates the accuracy of the four-class FCM clustering for the fat quantification in images acquired with different magnetic strengths, 1.5 T and 3.0 T. Our results demonstrated less accurate fat quantification compared to the results of the previous studies. In [19], the overall accuracy of an adipose tissue and muscle tissue segmentation in the T1weighted MRI images using the FCM clustering was 0.98 in

46

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

Fig. 7 – Top image: some results of the SAT internal-border detection using the FCM clustering, level-set approach, live-wire technique and the proposed method for the FSHD MRI images with various fat percentages (images in column 1, column 2 and column 3 correspond to patients P1 , P4 and P3 , respectively, which are described in Table 2). Bottom image: boxplot of the segmentation results similarity between the manual and automatic segmentation for the research dataset. The medians are denoted by solid black lines, the top and bottom box edges denote the first and third quartile, while the whiskers denote the largest and smallest values.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

terms of ZSI, while in [23], this accuracy was approximately 0.96. Nevertheless, the difference in the acquisition parameters and subjectivity of the manual delineation should be taken in consideration. In [34], it was reported that ZSI greater than 0.7 represents an “excellent” agreement in the medical image-analysis community. Therefore, our results were characterised by a medical specialist as good. Recent studies [42,43] suggest that the new MRI acquisition sequences optimised for the fat assessment, such as the least-squares estimation (IDEAL), three-point Dixon techniques, chemical-shift based water/fat separation and watersaturated sequences, may increase the contrast between the fat and non-fat tissues improving the performances of the fatquantification algorithms. Nevertheless, the left-/right-leg IMF asymmetry could be estimated fairly accurately in the presence of smaller segmentation errors (as are likely present in manual segmentations). A joint multivariate analysis of the water-suppressed and fat-suppressed T1-weighted MR images along with their non-suppressed counterparts may lead to an increased segmentation performance and may provide an improved separation between the muscle and fat voxels. The lack of knowledge of the actual oedema distribution makes it difficult to assess the actual amount of the oedema in the T2-STIR MRI images [44]. The manual delineation of the region of interest by experts in the neuroanatomical area has been the accepted standard, therefore, we believe that the obtained oedema quantification results, verified here by a medical specialists, are relevant for diagnosis. In future, the proposed segmentation method will be applied to perform segmentation of different muscle groups for more precise FSHD analysis.

7.

Conclusion

The paper introduces a computer-aided assessment for the facioscapulohumeral dystrophy diagnosis, where a fully automatic muscle-region segmentation method is proposed. The method refers to the live-wire approach for path search along the muscle-region border, where the path cost is calculated using an advanced edge-detection method. The obtained segmentation results are compared with the medical specialists manual delineations’ showing that the method is fairly insensitive to the fat percentage in the muscle region that is important for an accurate segmentation of the FSHD MRI images. We confirm that the T1-weighted MRI is a sensitive modality for the muscle-region segmentation and fat quantification for the FSHD diagnosis. Moreover, the quantification results of the oedema from the T2-STIR MRI images confirm the oedema presence in the fat-unaffected muscles that is an important parameter for the diagnosis of FSHD in an early stage. We believe that the proposed segmentation method may be of interest not only for the current problem, but also for a broad range of the related problems in the computer-based assessment during clinical diagnosis.

Acknowledgements This research was supported by the Research Program of the Slovenian Research Agency (PR-02600-1). We are also thankful

47

to the University Medical Centre of Ljubljana for providing the MRI data used in this research.

references

[1] C. Strom, P. Chan, V. Sulcova, D. Tsao, P. Walrafen, M. Braastad, C. Jocson, W. Sun, A. Anguiano, F. Boyar, DNA combing assay for detection of contraction-dependent facioscapulohumeral muscular dystrophy (FSHD1), Neurology 80 (PD3.005) (2013) 11–14. [2] G. Tasca, M. Monforte, E. Iannaccone, R. Frusciante, F. Laschena, P. Ottaviani, G.D. Lella, G. Galluzzi, F. Tiziano, A. Broccolini, M. Masciullo, G. Silvestri, E. Ricci, P2.38 lower limb muscle MRI in a large cohort of FSHD patients, Neuromuscul. Disord. 21 (2011) 671. [3] A.J. Skalsky, R.T. Abresch, J.J. Han, C.S. Shin, C.M. McDonald, The relationship between regional body composition and quantitative strength in facioscapulohumeral muscular dystrophy (FSHD), Neuromuscul. Disord. 18 (11) (2008) 873–880. [4] A.S. Gorgey, G.A. Dudley, Skeletal muscle atrophy and increased intramuscular fat after incomplete spinal cord injury, Spinal Cord 45 (4) (2007) 304–309. [5] D. Olsen, P. Gideon, T. Jeppesen, J. Vissing, Leg muscle involvement in facioscapulohumeral muscular dystrophy assessed by MRI, J. Neurol. 253 (11) (2006) 1437–1441. [6] H.E. Kan, T.W. Scheenen, M. Wohlgemuth, D.W. Klomp, I. van Loosbroek-Wagenmans, G.W. Padberg, A. Heerschap, Quantitative MR imaging of individual muscle involvement in facioscapulohumeral muscular dystrophy, Neuromuscul. Disord. 19 (5) (2009) 357–362. [7] R. Marcus, O. Addison, P. LaStayo, Intramuscular adipose tissue attenuates gains in muscle quality in older adults at high risk for falling. A brief report, J Nutr Health Aging 17 (3) (2013) 215–218. [8] B.H. Janssen, R.J.-W. Arts, N.B.-M. Voet, C.I.-H. Nabuurs, B.G.-M. Engelen, A. Heerschap, P2.37 MR as biomarker for disease evaluation in facioscapulohumeral dystrophy, Neuromuscul. Disord. 21 (9) (2011) 671. [9] S.D. Friedman, S.L. Poliachik, G.T. Carter, C.B. Budech, T.D. Bird, D.W. Shaw, The magnetic resonance imaging spectrum of facioscapulohumeral muscular dystrophy, Muscle Nerve 45 (4) (2012) 500–506. [10] S.L. Poliachik, S.D. Friedman, G.T. Carter, S.E. Parnell, D.W. Shaw, Skeletal muscle edema in muscular dystrophy: clinical and diagnostic implications, Phys. Med. Rehabil. Clin. N. Am. 23 (1) (2012) 107–122. [11] G. Tasca, M. Pescatori, M. Monforte, M. Mirabella, E. Iannaccone, R. Frusciante, T. Cubeddu, F. Laschena, P. Ottaviani, E. Ricci, Different molecular signatures in magnetic resonance imaging-staged facioscapulohumeral muscular dystrophy muscles, PLoS ONE 7 (06) (2012) 1–11. [12] A. Colin, E. Erbland, C. Datin, J.Y. Boire, A. Veyre, M. Zanca, Automatic muscle/fat quantification on MR images, in: Engineering in Medicine and Biology Society, 1995, IEEE 17th Annual Conference, September, 1995, pp. 479–480. [13] V. Barra, J.Y. Boire, Segmentation of fat and muscle from MR images of the thigh by a possibilistic clustering algorithm, Comput. Methods Programs Biomed. 68 (3) (2002) 185–193. [14] D. Wald, B. Teucher, J. Dinkel, R. Kaaks, S. Delorme, H. Boeing, K. Seidensaal, H.-P. Meinzer, T. Heimann, Automatic quantification of subcutaneous and visceral adipose tissue from whole-body magnetic resonance images suitable for large cohort studies, J. Magn. Reson. Imaging 36 (6) (2012) 1421–1434. [15] C. Wrslin, F. Springer, B. Yang, F. Schick, Compensation of RF field and receiver coil induced inhomogeneity effects in

48

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26] [27] [28]

[29]

[30]

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 2 0 ( 2 0 1 5 ) 37–48

abdominal MR images by a priori knowledge on the human adipose tissue distribution, J. Magn. Reson. Imaging 34 (3) (2011) 716–726. S.A. Sadananthan, B. Prakash, M.K.-S. Leow, C.M. Khoo, H. Chou, K. Venkataraman, E.Y. Khoo, Y.S. Lee, P.D. Gluckman, E.S. Tai, S.S. Velan, Automated segmentation of visceral and subcutaneous (deep and superficial) adipose tissues in normal and overweight men, J. Magn. Reson. Imaging 41 (4) (2014) 924–934. V. Positano, T. Christiansen, M.F. Santarelli, S. Ringgaard, L. Landini, A. Gastaldelli, Accurate segmentation of subcutaneous and intermuscular adipose tissue from MR images of the thigh, J. Magn. Reson. Imaging 29 (3) (2009) 677–684. J. Peijie, C. Yang, S. Wu, R. Xiangyan, The research of tissue segmentation method on thigh MRI, in: Multimedia Technology (ICMT), 2010 International Conference on, 2010, pp. 1–3. H. Kang, A. Pinti, A. Taleb-Ahmed, X. Zeng, An intelligent generalized system for tissue classification on MR images by integrating qualitative medical knowledge, Biomed. Signal Process. Control 6 (1) (2011) 21–26. P.K. Commean, L.J. Tuttle, M.K. Hastings, M.J. Strube, M.J. Mueller, Magnetic resonance imaging measurement reproducibility for calf muscle and adipose tissue volume, J. Magn. Reson. Imaging 34 (6) (2011) 1285–1294. E.N. Mortensen, W.A. Barrett, Intelligent scissors for image composition, in: Computer Graphics, SIGGRAPH Proceedings, 1995, pp. 191–198. F.M. Cox, M. Reijnierse, C.S.P. van Rijswijk, A.R. Wintzen, J.J. Verschuuren, U.A. Badrising, Magnetic resonance imaging of skeletal muscles in sporadic inclusion body myositis, Rheumatology 50 (6) (2011) 1153–1161. L. Urricelqui, A. Malanda, A. Villanueva, Automatic segmentation of thigh magnetic resonance images 3 (10) (2009) 964–971. L. Hong, Y. Wan, A.K. Jain, Fingerprint image enhancement: algorithm and performance evaluation, IEEE Trans. Pattern Anal. Mach. Intell. 20 (8) (1998) 777–789. S.S. Chikkerur, A.N. Cartwright, V. Govindaraju, Fingerprint image enhancement using STFT analysis, in: PROC. ICAPR, 2005, pp. 20–29. Z. Guo, R.W. Hall, Parallel thinning with two-subiteration algorithms, Commun. ACM 32 (3) (1989) 359–373. E. Dijkstra, A note on two problems in connexion with graphs, Numer. Math. 1 (1) (1959) 269–271. O. Malgina, A. Praznikar, J.F. Tasic, Inhomogeneity correction and fat-tissue extraction in MR images of facioscapulohumeral muscular dystrophy, Pattern Recogn. Lett. 34 (12) (2013) 1364–1371. G. OLaighin, B. Broderick, M. Clarke-Moloney, F. Wallis, P.A. Grace, A technique for the computation of lower leg muscle volume from MRI images in the context of venous return, in: Engineering in Medicine and Biology Society, 2007. EMBS 2007. 29th Annual International Conference of the IEEE, August, 2007, pp. 951–954. G. Frisullo, R. Frusciante, V. Nociti, G. Tasca, R. Renna, R. Iorio, A. Patanella, E. Iannaccone, A. Marti, M. Rossi, A. Bianco, M. Monforte, P. Tonali, M. Mirabella, A. Batocchi, E. Ricci, Cd8+ T cells in facioscapulohumeral muscular dystrophy patients with inflammatory features at muscle MRI, J. Clin. Immunol. 31 (2) (2011) 155–166.

[31] P. Reilich, N. Schramm, B. Schoser, P. Schneiderat, N. Strigl-Pill, J. Mller-Hcker, W. Kress, A. Ferbert, S. Rudnik-Schneborn, J. Noth, H. Lochmller, J. Weis, M. Walter, Facioscapulohumeral muscular dystrophy presenting with unusual phenotypes and atypical morphological features of vacuolar myopathy, J. Neurol. 257 (7) (2010) 1108–1118. [32] N. Baron, N. Kachenoura, F. Beygui, P. Cluze, P. Grenier, A. Herment, F. Frouin, Quantification of myocardial edema and necrosis during acute myocardial infarction, in: Computers in Cardiology, September 2008, 2008, pp. 781–784. [33] M. Abramoff, P. Magalhaes, S. Ram, Image processing with Imagej, Biophoton. Int. 11 (7) (2004) 36–42. [34] A.P. Zijdenbos, B.M. Dawant, R.A. Margolin, A.C. Palmer, Morphometric analysis of white matter lesions in MR images: method and validation, IEEE Trans. Med. Imaging 13 (4) (1994) 716–724. [35] P.J, De ziekte van Landouzy-Dejerine, University of Amsterdam, The Netherlands, 1983 (Ph.D. thesis). [36] C.A. Canapari, N. Barrowman, L. Hoey, S.W. Walker, E. Townsend, B.S. Tseng, S.L. Katz, Truncal fat distribution correlates with decreased vital capacity in Duchenne muscular dystrophy, Pediatr. Pulmonol. 50 (1) (2015) 63–70. [37] P.D. Kovesi, Matlab and Octave Functions for Computer Vision and Image Processing, Centre for Exploration Targeting, School of Earth and Environment, The University of Western Australia, 2000. [38] T. Shen, X. Huang, 3D medical image segmentation by multiple-surface active volume models, in: G.-Z. Yang, D. Hawkes, D. Rueckert, A. Noble, C. Taylor (Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI 2009, Vol. 5762 of Lecture Notes in Computer Science, 2009, pp. 1059–1066. [39] W. Hong, B. Georgescu, X.S. Zhou, S. Krishnan, D. Comaniciu, Database-guided simultaneous multi-slice 3D segmentation for volumetric data, in: ECCV, Springer-Verlag, 2006, pp. 397–409. [40] M. Hudelmaier, W. Wirth, M. Himmer, S. Ring-Dimitriou, A. Snger, F. Eckstein, Effect of exercise intervention on thigh muscle volume and anatomical cross-sectional areas-quantitative assessment using MRI, Magn. Reson. Med. 64 (6) (2010) 1713–1720. [41] S. Makrogiannis, S. Serai, K.W. Fishbein, C. Schreiber, L. Ferrucci, R.G. Spencer, Automated quantification of muscle and fat in the thigh from water-, fat-, and nonsuppressed MR images, J. Magn. Reson. Imaging 35 (5) (2012) 1152–1161. [42] T.A. Bley, O. Wieben, C.J. Franois, J.H. Brittain, S.B. Reeder, Fat and water magnetic resonance imaging, J. Magn. Reson. Imaging 31 (1) (2010) 4–18. [43] R. Kijowski, M.A. Woods, K.S. Lee, K. Takimi, H. Yu, A. Shimakawa, J.H. Brittain, S.B. Reeder, Improved fat suppression using multipeak reconstruction for ideal chemical shift fat–water separation: application with fast spin echo imaging, J. Magn. Reson. Imaging 29 (2) (2009) 436–442. [44] R.I. Johnstone, J.P. Greenwood, J.D. Biglands, S. Plein, J.P. Ridgway, A. Radjenovic, Assessment of tissue edema in patients with acute myocardial infarction by computer-assisted quantification of triple inversion recovery prepared MRI of the myocardium, Magn. Reson. Med. 66 (2) (2011) 565–574.