Spatial fuzzy c‐means thresholding for semiautomated calculation of percentage lung ventilated volume from hyperpolarized gas and 1H MRI

To develop an image‐processing pipeline for semiautomated (SA) and reproducible analysis of hyperpolarized gas lung ventilation and proton anatomical magnetic resonance imaging (MRI) scan pairs. To compare results from the software for total lung volume (TLV), ventilated volume (VV), and percentage lung ventilated volume (%VV) calculation to the current manual “basic” method and a K‐means segmentation method.

R espiratory diseases affect a large portion of the population, meaning development of sensitive imaging markers for diagnostic and prognostic assessment of lung disease is of growing importance. Ventilation-weighted hyperpolarized gas (HP) and anatomical proton ( 1 H) lung magnetic resonance imaging (MRI) can be used for quantitative evaluation of lung function, including detection of early obstructive changes. 1 The most commonly used quantitative measures of lung function derived from HP and 1 H scan pairs are lung ventilated volume percentage (%VV, the ratio of ventilated lung volume in HP images to total lung volume in 1 H images) and its counterpart the ventilation defect percentage (%VD 5 100%-%VV). 1 Calculation of %VV, therefore, requires segmentation of both HP and 1 H image sets and, with manual lung segmentation taking on the order of 1-2 hours, depending on image resolution, this presents a timeconsuming barrier for routine clinical application.
Development of methods to segment these scan pairs is thus necessary for quick, accurate, and reproducible quantitative analysis for clinical uptake of the methodology, while tackling technical challenges such as partial volume effects and motion artifacts from the heart. 2,3 Previous methods for HPG ventilation image segmentation have been based on manual intensity thresholding, 1 K-means clustering, 4,5 multiple atlas labeling, 6 and globally optimal graph cuts. 7 He et al recently developed a method of characterizing the distribution of ventilation via linear binning, 8 while Zha et al added an adaptive aspect of the K-means algorithm. 9 Segmentation of 1 H anatomical scans has been proposed using a seeded region-growing algorithm, 4 active contours within a closed homogeneous region, 10 and a multiple atlas labeling approach. 11 Most of these segmentation techniques 4,6-9,12 require little to no manual input, and the methods developed in Refs. 4 and 8 also grade ventilation. However, the K-means segmentation method 4,9 can fail, with low signal-to-noise ratio (SNR) images due to its binary clustering nature and inability to differentiate noise from lung tissue. A second potential disadvantage of the method developed previously 4 is the Gaussian filtering of the 1 H images, which may lead to underestimation of the total lung volume. The performance of the method developed in Ref. 7 also has limitations, with low SNR HPG images as the segmentation of the 1 H anatomical image is dependent on the HPG image itself. Although multiple atlas labeling 6,11 may be the most automated approach, its complexity and need for suitable prior information are disadvantages.
Fuzzy C-means (FCM) clustering, 12,13 which has been applied to HPG image segmentation, 12 is based on histogram information. The spatial Fuzzy C-means (SFCM) algorithm, developed previously, 14 incorporates the use of spatial information into the calculation of the membership function and, unlike the standard FCM algorithm, it also allows use of neighboring pixel information when classifying a voxel, leading to a more robust segmentation in the presence of noise and partial volume effects.
The purpose of this work was to modify the SFCM method 14 to segment both HPG and 1 H images of the lung and incorporate it into an image processing pipeline with high resilience to noise within a graphical user interface (GUI). Furthermore, in order to quantify the algorithm's performance, the secondary aim was to compare the outputs of this novel semiautomated approach to that of the current basic segmentation and a K-means based method for %VV calculation.

MATERIALS AND METHODS
Patients scanned with 3 He were analyzed with local Research Ethics Committee approval, while patients scanned with 129 Xe gave informed consent as part of a separate research study.

Imaging
All imaging was carried out on a GE HDx 1.5T MR scanner (GE Healthcare, Milwaukee, WI). 3D 1 H anatomical (spoiled gradient echo [SPGR]) and HP 3 He ventilation-weighted (balanced steadystate free-precession [bSSFP]) images were acquired during the same breath-hold. 2,15 3D HP 129 Xe ventilation-weighted bSSFP images 16 and 3D 1 H anatomical SPGR images were acquired in separate breaths. Sequence parameters are provided in the Supplementary Information.

Observers
Three observers (O1, O2, and O3) with 6, 1, and 5 years' experience in lung image segmentation, respectively, performed the analysis independently. O1, O2, and O3 analyzed 3 He scans and O2 and O3 analyzed 129 Xe scans, with both the basic method and the semiautomated method.

Participants
Six 3 He-1 H scan pairs ( 3 He SNR range 47-72) were selected from a database of patients with respiratory conditions of various severities. Patients' ages ranged from 23-68 years (three male, three female) and forced expiratory volume in 1 second (FEV 1 ) (% predicted) ranged from 24-63%. One patient suffered from horseshoe lung, one patient had asthma, and four patients had chronic obstructive pulmonary disease (COPD).
In addition, six 129 Xe-1 H scan pairs ( 129 Xe SNR range 18-34) from patients with lung cancer were analyzed to test the applicability of the method to HP 129 Xe ventilation images. Patient ages ranged from 62-76 (three male, three female), while FEV 1 (% predicted) ranged from 36-94%.

Image Analysis
The primary quantitative measure was %VV defined as ventilated volume (VV, from HPG ventilation-weighted images) divided by total lung volume (TLV, from 1 H anatomical images) 1 : The ventilated volume was masked by the total lung volume so that only voxels included in the TLV contributed to the calculated VV for all methods. 129 Xe and separate-breath 1 H acquisitions were coregistered using ANTs registration software. 17

Basic Method
Segmentations were carried out using in-house software written in MatLab (MathWorks, Natick, MA). Thresholding was used to segment the ventilation images. 1 H images were segmented manually using the ventilation images overlaid on the 1 H images.

Semiautomated Method
SUMMARY. A summary of the semiautomated image processing workflow is shown in Fig. 1. The GUI was implemented in Mat-Lab. Briefly, images were filtered using a bilateral filter with the filtering parameters tuned for 3 He, 129 Xe, and 1 H images. They were then clustered using the SFCM algorithm, using six clusters for 1 H anatomical images and four clusters for HP gas images. A single cluster was chosen that best represented the lung, then the membership values were thresholded to create a binary mask. Spatial and intensity membership weightings were both 1, 14 however, the pipeline will modify the spatial weighting to 2 if image SNR drops below 20.
ALGORITHM DEVELOPMENT: FCM VS. SFCM. The FCM algorithm assigns N pixels to C clusters via Fuzzy memberships (l). 12 These Fuzzy memberships do not take into account the spatial information, only intensity information. The SFCM algorithm, on the other hand, makes use of a window centered on each voxel of the image, 14 which incorporates the membership information of these voxels. This spatial information will then weight the membership function towards the correct cluster (with a specific weight being placed on the intensity and spatial memberships) only if the voxel was, for example, corrupted by noise, and would be incorrectly classified by the standard FCM method.
To decide whether to use FCM or SFCM as the clustering technique in this work, both FCM and SFCM segmentations were carried out on all 3He and 129Xe gas images and then assessed qualitatively (by visual inspection of features included/excluded by the algorithm) and quantitatively by comparing the unedited ventilated volume returned by the separate algorithms (the initial ventilation mask including airways).
Using the same number of clusters, the FCM method consistently included areas of low signal intensity deemed to be defects or noise via qualitative assessment. Additionally, there was a significant difference between the FCM and SFCM ventilated volumes for both 3 He and 129 Xe images (P 5 0.0312 for both) (see Supplementary Information for detailed results; Fig. 2, for example, from one patient). Therefore, SFCM was used in the final semiautomated method.  FILTERING. To improve robustness and resilience against image noise and artifacts, a bilateral filter was used. 18 The key property of this filter in the context of image segmentation is that edges are maintained due to the use of a nonlinear combination of range and domain filtering that weights pixel values depending on spatial and intensity similarity. Filter parameters that maintained ventilation defect integrity, smoothed artifact/noise, preserved edges, and ensured TLV was within an acceptable error margin (65%) of manual segmentation were empirically determined by processing 12 sets of HP gas ventilation and 1 H anatomical image pairs. The training images (from patients with COPD and healthy volunteers) used for this optimization process were acquired in the same way, but are separate from the data used to evaluate algorithm performance.
Different filter values and mask thresholds were required for 129 Xe and 3 He images due to differing imaging resolutions and SNR. 3 He images (SNR range 25-67) filter had a window size of 3 3 3 and spatial and intensity standard deviation of 3 and 0.15, respectively, with a binary mask threshold (membership threshold) of 0.1. For 129 Xe images (SNR range 20-45) the filter intensity standard deviation was increased to 0.2 and the binary mask threshold decreased to 0.05. For 1 H images the intensity standard deviation was reduced to 0.1 and a binary mask threshold of 0.15 was used. All processing was carried out in-plane.
AUTOMATED/MANUAL EDITING. 1 H masks had data outside the lung region removed by a border-clearing algorithm (see Supplementary Information for details). Main airways and vessels, and any regions of noise or tissue misclassified as lung volume, were removed manually from the masks using the software ITK-SNAP. 19

K-means Method
All datasets were analyzed using a modified version of the method developed by Kirby et al. 4 The size of the window used in the Gaussian filter for 1 H anatomical image segmentation was reduced from 15 3 15 to 3 3 3 and the standard deviation reduced to 0.01, the radius of the closing structuring element was reduced from 15 to 7, and data outside the lung region was removed by a border-clearing algorithm. No filtering of HP gas images was applied as per Ref. 4

Performance Evaluation
Performance analysis was carried out on %VV, VV, and TLV on a slice-by-slice basis. Intraclass correlation (ICC) using the two-way mixed model for absolute agreement, Bland-Altman analysis, and Spearman's or Pearson's correlation (R) were performed using GraphPad Prism (GraphPad Software, La Jolla, CA). Spatial comparisons were carried out on VV and TLV masks using the Dice Similarity Coefficient (DSC). T-tests and Wilcoxon signed rank tests were performed to assess statistical difference between analysis metrics.  Table 1 shows the correlation, Bland-Altman limits of agreement (LOA) and ICC for %VV computed by different observers and methods. Correlation improved (P 5 0.25) between observers when using the semiautomated method (%VV mean R 5 0.984) when compared to the basic method (%VV mean R 5 0.863), while mean ICC also increased (P 5 0.25) from 0.873 using the basic method to 0.980 using the semiautomated method. LOA (P 5 0.50) and %VV bias magnitude (P 5 0.25) were reduced when using the semiautomated method (%VV mean LOA 5 7.5%, mean jbiasj 5 2.3%) when compared to the basic method (%VV mean LOA 5 14.2%, mean jbiasj 5 4.6%). These improvements were also seen in the VV and TLV measures. DSC significantly improved using the semiautomated method (VV mean DSC 5 0.973, TLV mean DSC 5 0.980) when compared to the basic method (VV mean DSC 5 0.947, TLV mean DSC 5 0.957) (P < 0.01 for both VV and TLV DSC). The K-means method underestimated TLV when compared to both other methods. %VV was overestimated when compared to the basic (mean bias 5 5.0%) and semiautomated (mean bias 5 9.7%) methods (Fig. 3, for example). The Bland-Altman plot in Fig. 4a(iii) shows poor agreement of the K-means method %VV with the basic method for O2 and is representative of the pattern seen when comparing K-means with both the basic (%VV mean LOA 5 29.8%) and semiautomated (%VV mean LOA 5 28.2%) methods for all observers.

He
On average, the semiautomated method underestimated %VV by 4.6% when compared to the basic method carried out by the same observer, with a mean LOA of 19.7%. The semiautomated method reduced average segmentation time from 1 hour (basic) to 25 minutes.

Xe
Correlation, LOA, bias magnitude, and DSC improved between observers when using the semiautomated method when compared to the basic method ( Fig. 4b(i,ii)). The Kmeans method underestimated TLV and overestimated %VV compared to the other methods to a greater extent than for the 3 He data. The semiautomatic method underestimated %VV by 2.3% compared to the basic method for O2 ( Fig. 4b(iv)) and overestimated %VV by 18.6% for O3. The mean LOA for %VV calculated by the same observer between the basic and semiautomatic methods was 26.4%.

DISCUSSION
The semiautomated image processing workflow developed reduced interobserver variability, a problem in longitudinal imaging studies when multiple observers may be required. The use of the coregistered dual 3 He-1 H image acquisition 2 in this work circumvents the need for image registration, which is commonly used in %VV analysis by other groups. 4,8,20 However, image registration was required for the 129 Xe image analysis.
The semiautomated method results were more similar to those of the basic method than those obtained by the fully automated K-means method. 4 Remaining differences in %VV values returned by the basic and semiautomated method could be explained by the basic method including pixels that correspond to tissue on the 1 H image but where 3 He signal is also present, while the proposed method considers the edge of the lung from the 1 H image as the ground truth to exclude those pixels.
In the development phase of the semiautomated method, filtering followed by SFCM clustering was found to be robust to choice of imaging sequence and parameters. The technique returned similar masks of ventilated volume from HP gas images acquired with different sequences (SSFP/SPGR) and parameters (TE, TR) from the same subjects (see supplementary information). Semiautomatic segmentations of both the HP gas and proton images analyzed here were consistently of good quality, from images acquired with a range of fields of view (FOVs) (36-40 cm) and SNR (18-72) from patients with a variety of different diseases. Beyond removal of the main airways and vessels, little manual editing of the masks was required.
The K-means method underestimated TLV when compared to both other methods due to Gaussian filtering of the 1 H images, despite the filtering being substantially reduced in this work. This effect was exacerbated by the lower resolution of the 1 H images paired with the 129 Xe images. In addition, the original ventilation masks returned by the K-means method classified regions of noise as ventilated lung tissue. Both these factors lead to a systematic overestimation of %VV by the K-means method, which has implications for its use calculating outcome measures in clinical studies.
The limitations of this technical development study are the small numbers of patients analyzed as well as the reduced number of observers who segmented the 129 Xe images and the lack of comparison to other established techniques for %VV calculation.
In conclusion, the method presented here provides a robust and repeatable means of semiautomated lung MRI segmentation and was demonstrated on both 3 He and 129 Xe ventilation images. The method proposed improves interobserver agreement, and is easier and quicker to use than the current basic segmentation used to calculate lung ventilation volumes.