Novel Subtypes of Pulmonary Emphysema Based on Spatially-Informed Lung Texture Learning: The Multi-Ethnic Study of Atherosclerosis (MESA) COPD Study

Pulmonary emphysema overlaps considerably with chronic obstructive pulmonary disease (COPD), and is traditionally subcategorized into three subtypes previously identified on autopsy. Unsupervised learning of emphysema subtypes on computed tomography (CT) opens the way to new definitions of emphysema subtypes and eliminates the need of thorough manual labeling. However, CT-based emphysema subtypes have been limited to texture-based patterns without considering spatial location. In this work, we introduce a standardized spatial mapping of the lung for quantitative study of lung texture location and propose a novel framework for combining spatial and texture information to discover spatially-informed lung texture patterns (sLTPs) that represent novel emphysema subtype candidates. Exploiting two cohorts of full-lung CT scans from the MESA COPD (n = 317) and EMCAP (n = 22) studies, we first show that our spatial mapping enables population-wide study of emphysema spatial location. We then evaluate the characteristics of the sLTPs discovered on MESA COPD, and show that they are reproducible, able to encode standard emphysema subtypes, and associated with physiological symptoms.

Preliminary CT-based clinical studies suggest that regional analysis will be instrumental in advancing the understanding of multiple pulmonary diseases [19]. In the case of emphysema, it is suspected that different emphysema subtypes affect the lungs in specific anatomical regions. But the problem of how many subtypes exist, how they evolve in time and how they vary with spatial (anatomical) location is still unsolved. To date, categorization of emphysema on CT images has relied only on analysis of local textural patterns, using either grey-level co-occurrence matrix (GLCM) features [12], [16], texton features [13], [14], or local binary pattern (LBP) features [11]. All these approaches use intensity information without consideration of spatial location.
In two previous studies [17], [18], we proposed to use local textural patterns to generate unsupervised lung texture patterns (LTPs) followed by LTP-grouping based on their spatial co-occurrence in local neighborhoods. Such separate use of intensity and spatial information cannot guarantee spatial and textural homogeneity of the final LTPs.
In this study, we propose to perform discovery of LTPs via unsupervised clustering of joint spatial and textural information of local texture patterns. Spatial information can be inferred from crude partitioning of the lung with subdivisions of Cartesian coordinates or by segmenting the lung into zones (e.g. upper, lower) [4] or lobes [20]. However, such approaches have limited spatial precision and lack relative information such as peripheral versus central positioning, which is important in defining paraseptal emphysema and subpleural bullae. We introduced in [21] a new standardized lung shape spatial mapping, called Poisson distance conformal mapping (PDCM), which enables detailed, precise and standardized mapping of voxel positions with respect to the lung surfaces. This paper further refines the PDCM algorithm and exploits it for the study of emphysema spatial patterns across populations of CLE-, PLE-and PSE-predominant subjects. This paper also provides an exhaustive description of the framework for combining spatial and texture information in the unsupervised discovery of emphysema-specific texture patterns, which are called spatially-informed LTPs (sLTPs).
Exploiting a cohort of 317 full-lung CT scans from the MESA COPD study [4], and 22 longitudinal CT scans from the EMCAP study [22], the discovered sLTPs are extensively evaluated as emphysema subtype candidates in terms of reproducibility with respect to training sets, labeling task and scanner generations, ability to encode standard emphysema subtypes, and associations with respiratory symptoms. A graphical pipeline of the learning and evaluation steps is provided in the Supplementary Material.

1.
Generate spatial mapping of the lung masks: mapping voxels within the lung masks into a custom Poisson distance map (PDM) to encode the "peel to core" distance, and a conformal mapping (PDCM) to distinguish superior versus inferior, anterior versus posterior and medial versus lateral voxel positions;

2.
Encode regions of interest (ROIs) within emphysema-like lung: sampling ROIs from emphysema segmentation masks, and generating spatial features (based on spatial mapping) and texture features of each ROI;

3.
Discover an initial set of LTPs: clustering training ROIs into a large number of clusters, based on texture features, and then iteratively augment the LTPs with spatial information via regularization;

4.
Generate the final set of sLTPs: measure the similarity between LTPs in the initial set, group similar / redundant LTPs and generate the final set of sLTPs via partitioning a similarity graph.
We now detail these four steps individually.

B. Spatial Mapping of the Lung Masks
To generate spatial mapping of the lung masks, we first use the concept of Poisson distance map (PDM), introduced in [23], to encode the shape of individual lung masks V. PDM is commonly used for characterizing the silhouette of an object via continuous labeling of voxel positions with scalar field values U 3d in the range of [0, 1]. In our case, the field value U 3d encodes the "peel to core" distance between a given voxel and the external lung surface ∂V. This field is computed by solving the following Poisson equation: ΔU 3d x, y, z = − 1, for x, y, z ∈ V subject to U 3d x, y, z = 0, for x, y, z ∈ ∂V (1) where ΔU = U xx + U yy + U zz is the Laplacian operator based on 2nd-order spatial derivatives along x, y, z.
The solution for U proposed in [23] is guaranteed to be smooth according to [24]. It has the advantage of generating distance values that are sensitive to global shape characteristics, unlike other distance metrics (e.g. Euclidian or Metropolis distances) which exploit single contour points. PDM can therefore reflect rich shape properties of the lung.
The core of the PDM is the set of voxels (one or very few) where U 3d (x, y, z) = 1. The PDM generated from a lung surface generally exhibits nice star-shaped profiles when viewed in axial cuts, with maxima near the center. On the other hand, core positions can vary greatly among subjects along the superior-inferior axis, due to variable morphologies of the lungs, especially near the heart and at the base. We illustrate an example in Fig. 1 (b) where the PDM generated with Equation (II-B) has core point(s) located close to the base of the lung rather than concentrated toward the middle of the longitudinal axis. We propose the following calibration of lung PDMs to (1) prevent U 3d = 1 in most apical and basal regions, and (2) enforce U 3d = 1 in a large range of axial slices. This makes PDM values numerically more consistent between subjects over a comparable range of axial slices. We denote U 3d max i the maximal in-slice value of U 3d (.,.,i), where the i th axial slice index is counted from the apex. We denote i v % the highest slice index value such that the total lung volume sumed over all slices with lower indices is < V% of the total lung volume. A normalized version (denoted as U 2d ) of the original PDM U 3d , is then defined, per axial slice index i, as U 2d . , . , i = U 3d . , . , i /U 3d max i . We further define U mod by combining U 3d and U 2d values, as follows: U mod ( . , . , i) = U 2d ( . , . , i), ∀i u ⩽ i ⩽ i d U mod ( . , . , i) = U 3d ( . , . , i)/U 3d max i u , ∀i < i u U mod ( . , . , i) = U 3d ( . , . , i)/U 3d max i d , ∀i > i d (2) with i u (resp. i d ) the smallest (resp. highest) slice index where U 3d max reaches a local maximum. To ensure that a consistent portion of the lung is included in [i u , i d ] we further enforce: if i u > i 25% then i u = i 25% (resp. if i d < i 75% then i d = i 75% ). We illustrate in Fig. 1(c) an example where U mod = 1 over a large range of axial slice indices and exhibits decreasing values when moving toward the apex or the base of the lung.
To equip the PDM with a coordinate system, we set the final core coordinate center point as the point on the axial slice index i 50% where U mod (x, y, i 50% ) = 1 and closest to the 2D center of mass of the axial lung mask (in case of multiple candidates, we would select one abitrarily, but such situation was not encountered on our dataset.).
To uniquely encode 3D voxel positions, we define radial values r = 1 − U mod and add conformal mapping of voxels positions onto a sphere, generating a Poisson distance conformal map (PDCM). We encode superior versus inferior, anterior versus posterior and medial versus lateral voxel positioning via latitude and longitude angles (θ, ϕ) with respect to the PDM core defined above and standard image axis. The generation of the spatial PDCM mapping is illustrated in Fig. 1(d).
The PDCM spatial mapping will be exploited for sLTP learning, and also to study population-based spatial location of emphysema, as reported in Section III-B.

1) Prior Emphysema Segmentation and ROI Sampling:
Texture and spatial analysis is performed within local ROIs centered on a subset of lung voxels. Sampling ROIs from emphysema-like lung requires prior emphysema segmentation. In this study, we exploited a training cohort of full-lung CT scans and their associated emphysema masks, which are generated using both a thresholding-based voxel selection and a hidden Markov measure field (HMMF) segmentation [25]. For thresholding, voxels with attenuation below −950 HU are selected. The HMMF segmentation enforces spatial coherence of the labeled emphysematous regions, and relies on parametric modeling of intensity distributions within emphysematous and normal lung tissues to adapt to individual and scanner variability. Percent emphysema measures the proportion of emphysematous voxels within the lung region, and is denoted %emph −950 or %emph HMMF , depending on the emphysema segmentation method.
In preliminary implementations, we tested several options for ROI sampling such as keypoint sampling in [17] and regular sampling in [18]. In this study, we use the systematic uniform random sampling (SURS) strategy as suggested in [26] for use on lung CT scans. Each individual lung mask is randomly sampled via dividing the bounding box of the lung into 3D stacks, and then selecting voxels per stack with a random shift of positions. Two parameters are used for the sampling: β 1 is used for the random shift of positions and β 2 is used to set the number of sampled voxels per stack. The SURS sampling ensures even representation of all lung regions while introducing variability in the position of sampled points with the random shift parameter β 1 . Only ROIs with both percent emphysema %emph −950 > 1% and %emph HMMF > 1% are retained for training to ensure sufficient representation of emphysematous regions (i.e. each training ROI has a minimal proportion of emphysema but can be a mixture of normal and emphysematous tissues).

2) Texture Features:
We use texton-based texture features to characterize each ROI, which model textures as the repetition of a few basic primitives (called textons), and were shown to outperform other texture features in unsupervised lung texture learning in [18]. A texton codebook is constructed by retaining the cluster centers (textons) of raw pixel representations of small-sized training patches. The clustering is performed with K-means.
By projecting all small-sized patches of a ROI onto the codebook, the texton-based feature of the ROI is the normalized histogram of texton frequencies.

3) Spatial Features:
To generate spatial features of individual ROIs, we divide the lung masks into lung sub-regions by discretizing our continuous lung shape spatial mapping with a minimal granularity. We divided r ∈ [0, 1] into 3 regular intervals to distinguish pleural from mid from core regions, divided θ ∈ [0, 2π] into 4 regular intervals to distinguish anterior, medial, posterior and lateral regions, and divided ϕ ∈ [−π/2, π/2] into 3 regular intervals to distinguish inferior, mid-level and superior regions. The spatial feature of each ROI is a one-hot vector indicating the lung sub-region it belongs to. Ordering of the bins that represent the sub-regions is done via arbitrary spatial rastering as no assumption needs to be made on spatial adjacency of adjacent bins.

D. Initial Augmented LTPs
We formulate the discovery of spatially-informed lung texture patterns (sLTPs) as an unsupervised clustering problem. One key factor in unsupervised clustering is the choice of the number of clusters. The algorithm is expected to find finer-grained emphysema subtype candidates than the three standard ones. Therefore, the number of clusters should be large enough to handle the diversity of textures encountered in the lung volumes (i.e. good intra-cluster homogeneity), and small enough to avoid redundancy (i.e. good inter-cluster differences) for clinical interpretation. Fixing a priori the number of clusters may prevent the discovery of rare patterns. We therefore propose a two-stage learning strategy, where we first generate an arbitrary large number of fine-grained lung texture patterns (LTPs), and then group similar LTPs to produce the final set of sLTPs, according to a dedicated metric.
LTPs {LTP k } ({·} denotes a set of variables hereafter) are characterized by their texture and spatial feature centroids (F T LT P k , F S LT P k ), which are encoded as histograms via averaging over assigned ROIs. An initial set of LTPs is generated by clustering with texture features, and is then augmented with spatial regularizations via iterative updates of the centroids and the ROIs assignements as described in Algorithm 1 and using the following mixed χ 2 -ℓ 2 similarity metric to enforce spatial concentration of LTPs while preserving their intra-class textural homogeneity: where P 95 denotes the 95 th percentile, Λ LT P k (t) denotes the set of ROIs that are labeled as LTP k at iteration t and ∧ LT P k t λ, W , γ * denotes the optimal labeling identified with a set of parameters {λ, W, γ} and the centroids updated at iteration t − 1. Designing proper distance metrics for histograms plays a crucial role in many computer vision tasks. Two popular choices are the χ 2 and the ℓ 2 distance metrics. The latter equally weights distances of all bins and is favored to compare one-hot vectors, while the former is a weighted distance favored to compare probability distributions. For the texture feature histograms that encode distributions over textons the first distance metric χ 2 (·) measures the χ 2 distance between the textural features of a ROI x and the centroid of LTP k . For the spatial features that are sparse one-hot vectors for individual ROIs, the second distance metric ⋅ 2 2 measures the ℓ 2 distance between the spatial features of a ROI x and the centroid of LT P k . A textural penalty term is then introduced as the third term, where 1 is the indicator function. Update of LTP centroids (step 2 in Algorithm 1) is performed after relabeling each ROI with the LTP to which it has the smallest weighted feature distance without turning on the textural penalty.

Parameter W:
This parameter is used to scale contributions between textural and spatial feature distances so that λ can be tuned within a small range of values. We defined it as: where SST T and SST S are respectively the texture and spatial total sum-of-square distances, computed on the whole N training ROIs to measure the overall diversity of texture and spatial features.
Parameter λ: This parameter controls the spatial regularization which will inevitably decrease textural homogeneity of individual LTPs. The value of λ is set as follows. First we define SSW T as the initial sum-of-square intra-cluster homogeneity of texture features without spatial regularization: Then we define SSW T λ as the SSW T measured on augmented LTPs with spatial regularization enforced with λ ∈ [0, 2]. Final value of λ is set to: In the context of unsupervised discovery, we hereby spatially regularize the augmented LTPs via an empirically acceptable textural homogeneity loss with the threshold L T (set based on data observations, as reported in Section III).
Parameter γ: This parameter weights the textural penalty term which is used for ROI labeling. We set γ = ∞ to prevent a ROI from being labeled to a spatially preferred but texturally dissimilar LTP.

E. Final sLTPs
In this final step, we generate sLTPs by partitioning a weighted undirected graph G where nodes are the N LT P initial augmented LTPs. As in [21], we define the edge weight between nodes i and j as the average replacement ratio of training ROIs relabeled from label LT P i to LT P j if LT P i is removed from the set of centroids and vice versa. In the replacement task, a ROI with a textural distance to the LT P k centroid exceeding the maximal intra-cluster textural distance of LT P k is not re-labeled. To prevent weak associations of LTPs that are not easily replaceable, we remove edges with weights lower than 0.5 (i.e. 50% replacement). Indeed, graph partitioning tends to preserve nodes that are not connected, which in our case would correspond to LTPs that are not easily replaced by other ones in the labeling task, hence not redundant. We use the Infomap algorithm [27] to partition the similarity graph G. As part of its optimization process that minimizes the description length of the network, Infomap selects an optimal number of clusters of aggregated LT P k which define our final sLTPs. Final texture and spatial centroids of the sLTPs are then computed utilizing the training ROIs labeled in our final {LT P k }.

F. Labeling of CT Scans With sLTPs
In the test stage, scans in the whole dataset are labeled by extracting sample points and their ROIs {x}. Since it is computationally prohibitive to evaluate the textural and spatial features on every voxels within the lung masks, we only label centers of ROIs densely sampled using again SURS. Sampled ROIs with %emph −950 ⩽ 1% or %emph HMMF ⩽ 1% have their center labeled as no-emphysema class. Remaining sampled centers get a sLTP label, via minimization of the following cost metric: Non-sampled voxels are labeled with the sLTP index of the nearest sampled center point via nearest neighbor search within the lung mask (i.e. using a Voronoi diagram). Labeling lung scans with the discovered sLTPs generates histograms of sLTPs, which are efficient lung texture signatures exploited for several tasks, as described in the evaluation sections.

G. Visualization of the sLTPs Spatial Density
To study the spatial distribution of sLTPs, we generate spatial visualization by scatter plotting of voxels labeled with individual sLTPs in sagittal projections, as follows.
We first randomly sample an initial set of ROIs over each lung via SURS sampling. Each ROI is associated with its center point coordinates (r, θ, ϕ) in the PDCMs. To avoid artificial higher densities on the scatter plot in regions close to the core, we adapt the number of ROIs selected per radial regions. The r values are binned into N r intervals with midpoint values r 1 , … , r N r to generate isovolumetric subvolumes of the lung. We then define the sub-sampling ratio α i = r i /r N r (which approximates the ratio of areas in the scatter plot) and set the number of ROIs sampled per r bin to N IsoV i = α i ⋅ N IsoV where N IsoV is a pre-set number of ROIs sampled in the outermost part of the lung.
All ROI centers in the sub-sampled set are converted to (x, y, z) Cartesian image coordinates and accumulated in a sagittal single plane, by setting x = 0. Final density plots of sLTPs are shown in projected radial coordinates r′ = y 2 + z 2 and ϕ′ = atan z/y . We color code each point on the sagittal projection with the following density measure: where Λ (r′, ϕ′) denotes the set of ROIs at (r′, ϕ′) positions. The numerator (first term) in Equation (8) measures the probability of sLT P k at projected position (r′, ϕ′), and the denominator (second term) measures the observed overall probability of (r′, ϕ′) to host any sLT P i .

A. Data
The data used for evaluation consists of full-lung CT scans of 317 subjects. All subjects had underwent CT scanning in the MESA COPD study [4], between 2009-2011. In addition, 22 out of the 317 subjects underwent CT scanning in the EMCAP study [22], between 2008-2009.
For the MESA COPD study, all CT scans were acquired at full inspiration with either a Siemens 64-slice scanner or a GE 64-slice scanner, at 120 kVp, speed 0.5 s, and current (mA) set according to body mass index following the SPIROMICS protocol [28]. Images were reconstructed using B35/Standard kernels with axial pixel resolutions within the range [0.58, 0.88] mm, and 0.625 mm slice thickness.
For the EMCAP study, scans were acquired with a Siemens 16-slice scanner, at 120 kVp, speed 0.5 s, and a current between 169 mA and 253 mA. Images were reconstructed using the B31f kernel with axial resolutions within the range [0.49, 0.87] mm, and 0.75 mm slice thickness.
Emphysema subtypes and severity have previously been assessed visually in the MESA COPD study (details available in [4]). The raters included four experienced chest radiologists from two academic medical centers. They assessed emphysema subtypes on CT scans by assigning a percentage of the lung volume affected by CLE, PLE and PSE respectively. Based on [4], N = 205 subjects do not exhibit emphysema, and are used here as the control set of no emphysema (NE) subjects. The remaining N = 112 subjects exhibit light (N = 53) or mild-to-severe (N = 59) emphysema. For these subjects, predominant emphysema subtype is defined as the subtype affecting the greatest proportion of the lungs.
In addition, the following clinical characteristics are available for the scans in MESA COPD study (details in [4]): demographic factors (age, race, gender, height, weight); forced expiratory volume in 1 second (FEV1); MRC dyspnea scale measure (5-level scale); six-minute walking test (6MWT) total distance; pre (baseline) 6MWT pulse oximetry; post 6MWT pulse oximetry; reported post 6MWT fatigue; and reported post 6MWT breathlessness. We used these measures for evaluating the clinical significance of the discovered sLTP.

B. Population Evaluation of Emphysema Using PDCM
We first demonstrate the ability of our proposed PDCM lung shape mapping to study the spatial patterns of emphysema over a population of subjects (cf. Fig. 2). For each scan in MESA COPD study, PDCM maps of voxels inside individual lungs are generated, attributing to each voxel a coordinate (r, θ, ϕ). Voxel intensity values in PDCM maps are then averaged and visualized along two types of projections:

2.
Radial projections: intensity values averaged over all angular directions at a subset of N r = 60 regular radial positions r 1 , … , r N r .
An illustration of these two PDCM intensity projections on a sample lung are visualized in Fig. 2 (a).
Population-average PDCM angular and radial intensity projections over subjects without emphysema (NE) are displayed in Fig. 2 (b). The averaged angular projection shows a clear pattern of lower attenuations (i.e. intensity values) in the anterior versus posterior region, which agrees with the intensity gradient due to gravity-dependent regional distribution of blood flow and air [29], [30]. The averaged radial projection shows a slight gradient from core to peel regions, which is likely due to the inclusion of voxels belonging to the mediastinal and costal pleura inside the lung mask.
Population-average PDCM intensity projections over subjects with CLE-, PLE-, and PSEpredominant emphysema subtypes are visualized in Fig. 2 (c). To highlight differences with respect to the control set, we display relative values after subtraction of the values from the corresponding NE average projection in Fig. 2 (b). Color coding represents relative intensity differences with more emphysema (more negative attenuation values) corresponding to the red color.
We can see on the relative angular PDCM intensity projections that regions of normal attenuation (green to blue) are absent for PLE-predominant subjects, whereas CLE-and PSE-predominant subjects appear to have emphysema regions (red) concentrated in the superior part. The average relative radial PDCM intensity projections on emphysema subjects show systematic lower attenuation values consistent with more emphysema in the core part for CLE-predominant subjects and more emphysema in the peel part for PSE-predominant subjects.

C. Qualitative Evaluation of Discovered sLTPs
For the discovery of sLTPs, 3/4 of the total scans in MESA COPD study (N = 238) were used for training, using random stratified sampling without replacement, while the other scans (N = 79) were used for testing. We summarize the setting of pre-defined parameters for the sLTP learning in TABLE I. In addition, spatial regularization weight λ is set via empirical tuning using Eq. (II-D). Based on the relative texture homogeneity loss measure ΔSSW T , we chose L T = 1% which corresponds to λ = 1.52, above which ΔSSW T increases drastically.
A total of 12 sLTPs were discovered using the full training set, and were used to label both the training and test scans in emphysema-like lung. Each sLTP was detected (i.e. %sLT P k > 0) in at least 5% of scans both in training and test sets. In Fig. 3, we illustrate in (a) the sLTP labeling of two sample CT scans; and in (b) the characteristics of each sLTP via visual illustrations of labeled patches, average occurrence in MESA COPD scans, and spatial distribution of their occurrence within the lungs. For the patch illustrations, 9 samples were randomly selected from all available labeled ROIs (see the Supplementary Material for high-resolution illustrations). For the average occurrence, we averaged %sLT P k values over scans with %sLT P k > 0. For the spatial distributions, we generated spatial scatter plots of sLTP locations from labeled ROIs, following the method described in II-G, with N IsoV = 5, 000, and N r = 60.
We can observe that patches belonging to an individual sLTP appear to be textually homogeneous. sLTP 1 and 4 show clear spatial accumulation in superior (apical) regions, sLTP 3, 5 and 7 in anterior regions, and sLTP 10, 11 and 12 in posterior regions. The brightest LTPs (11 and 12) have very distinct visual appearance and resemble combined pulmonary fibrosis emphysema (CPFE). Since we jointly enforce spatial prevalence and textural homogeneity, some sLTP can have spatial "outliers" that are texturally favored. All sLTPs returned similar occurrences in training and test sets. Some sLTPs are rare, such as sLTP 12 which covers ~1% of the lungs when present, but is still found in 24 scans over the whole MESA COPD cohort.

D. Reproducibility of sLTPs 1) Reproducibility of sLTP Labeling Versus Training Sets:
To test the reproducibility of sLTPs learning, we first compare the N sLTP = 12 sLTPs {sLT P k } generated with the full set of training scans, to N set = 4 sLTPs sets sLT P k c (c = 1, 2, 3, 4) using subsets of training data by eliminating via stratified subsampling 25% of the training scans without overlap on the left-out scans. Reproducibility of sLTPs is evaluated on the ROI labeling task, by computing the average overlap of labeled test ROIs with the following metric: R ln = 1 N set ⋅ N sLTP c 1 N set k 1 N sLTP |Λ sLT P k Λ π(sLT P k c ) | |Λ sLT P k | (9) where Δ sLT P k denotes the set of ROIs labeled with sLT P k , and π() denotes the permutation operator on the {sLT P k c } determined by the Hungarian method [31] for optimal matching between sets {sLT P k } and {sLT P k c }.
Compared with the N sLTP = 12 sLTPs learned on the full training set, we discovered N sLTP c = 12, 12, 13, and 13 sLTPs on training subsets. We obtain an overall labeling reproducibility measure of R ln = 0.91 which corresponds to a high reproducibility level.
We then further compute the reproducibility measure, denoted as R ln ′ , among training subsets. The metric is similar to Equation 9, replacing {sLT P k } and sLT P k c with sLTPs sLT P k c1 and sLT P k c2 (c1 ≠ c2) learned on different training subsets. We obtain an overall labeling reproducibility measure of R ln ′ = 0.85 (standard deviation = 0.07).
To evaluate the contribution of spatial features in sLTP learning, we further generate sets of lung texture patterns using only texture features (i.e. using initial LTPs without spatial augmentation in Section II-D, and setting λ = 0 for the replacement test in Section II-E).
We discovered 11 patterns using the full training set, and 11, 11, 12 and 12 patterns on training subsets. The reproducibility measures R ln and R ln ′ equal to 0.84 and 0.78 (standard deviation = 0.12), are lower than the ones obtained using the proposed sLTP learning, hence confirming the benefit of adding spatial features.
2) Reproducibility of sLTP Labeling Versus ROI Sampling: As detailed in Section II-F, sLTP labeling is based on a subset of voxels setting ROI positions, using SURS-based sampling strategy, which is controlled with the parameter β 2 (number of samples per stack).
The selected ROIs have an influence on the final outline of the label map, which is hopefully minor if ROIs are sampled densely enough and if sLTPs are generic enough. In this experiment, we test this hypothesis by generating two different sets of ROIs on test scans using two different random seedings, and measure the reproducibility of the generated label masks using the {sLT P k } discovered on the full training set, while varying the β 2 parameter.
We measure labeling reproducibility using the two sets of ROIs with the following metrics: DC (sLT P k , β 2 ) = average of Dice coefficients of label masks of sLT P k over all test scans; • R 1a CC (sLT P k , β 2 ) = Spearman correlation coefficients of %sLT P k values within the lungs over all test scans.
We illustrate in Fig. 4 (a), the average, max and min values of R la * measures overall {sLT P k }, for β 2 ∈ [1,20]. Both reproducibility measures increase with β 2 in an exponential manner. We obtain an average R la DC > 0.8 when β 2 > 10, corresponding to sampling less than 0.05% points in each stack. We obtain an average R la CC > 0.9 when β 2 > 5. Minimum R la values always occur for sLTP 12, which is the rarest sLTP, as reported in Section III-C. We used the 12 sLTPs discovered on the full MESA COPD training set. Because of differences in scanner generations (axial CT in EMCAP versus spiral CT in MESA COPD) and radiation dose settings, intensity calibration was required, implemented in two steps: 1) equalizing the outside air mean intensity value (according to [25]); 2) histogram mapping of normal lung parenchyma identified with the HMMF-based emphysema masks. The sLTPs 2 to 12 were found to be present in both datasets, but sLTPs {2, 3, 4, 12} occur in less than 6 pairs of scans. We report in Fig. 4 (b) the Cohen's Kappa coefficients of sLT P k presence for sLTPs 2-12, and the Spearman correlation coefficients of %sLT P k for the frequent sLTPs only (sLTPs 5 to 11). The Cohen's Kappa coefficients and Spearman correlations are all above 0.8, which confirms robust sLTP presence and percentage labeling on the 22 subjects scanned on different scanner types in two studies.

E. sLTPs' Ability to Encode Standard Emphysema Subtypes
When generating unsupervised lung texture patterns (either sLTPs in this work or earlier generations of LTPs in previous work), we expect them to be finer-grained than the three standard emphysema subtypes used in [4], while still capable to encode them, hence linking unsupervised image-based emphysema subtyping with clinical prior knowledge.
The (s)LTPs (either LTPs or sLTPs) can correspond to a single standard subtype or a mixture of those. We hereby evaluate the ability of the generated (s)LTPs to predict the overall extent of standard emphysema subtypes. To do this, we generate, for each scan and per lung, two signature vectors: 1) a (s)LTP signature histogram composed of the percentage of non-emphysema class (obtained as in Section II-F) and the percentages of individual (s)LTPs in the emphysema-like lung. This normalized histogram is called the (s)LTP predictor signature and is of size N predictor = N s LT P + 1; 2) a ground-truth signature composed of the percentage of non-emphysema and the three standard emphysema subtypes (CLE, PLE, PSE), as visually evaluated in [4]. A constrained multivariate regression model is used on labeled training scans to learn regression coefficients between the (s)LTP and ground-truth signatures, using the following optimization: argmin A XA − Y 2 2 s.t. 0 < A (k, i) < 1 and i A k i 1 (10) where X N scan × N predictor is composed of all training (s)LTP signatures in N scan training scans, and Y N scan × 4 contains the ground-truth signatures. A N predictor × 4 is the matrix of regression coefficients {A k,i }, which measure the probability of a voxel labeled as a certain predictor belonging to one of the ground-truth classes, and are therefore constrained to be in the range of [0, 1]. Optimization of regression was solved using the CVX toolbox [32] Quality of prediction is measured with the intraclass correlation (ICC) between predicted and ground-truth exploiting the full MESA COPD dataset. We use a 4-fold cross validation (3/4 label masks used for training the regression and 1/4 used for testing and measuring prediction quality). Significance of differences in ICC values was assessed using Fisher's r-to-z transformation and a two-tailed test of the resulting z-scores.
In Fig. 5, we compare prediction quality with 7 sets of emphysema-specific (s)LTPs (re)trained on the same set of emphysematous ROIs: 1) the 12 sLTPs learned in this study; 2-3) the initial set of 100 LTPs generated in this study before (denoted as LTP init-T) and after (denoted as LTP init-TS) spatial augmentation; 4) LTPs generated by one-stage clustering (denoted as LTP TS) of the proposed texture and spatial features, by setting N LT P = 12 directly (this is to test the contribution of the proposed two-stage learning in Section II-D); 5-6) LTPs re-generated using Method A [17], discovered via graph partitioning of 100 candidates based on local spatial co-occurrence and with N LT P = 8 as in [17] or 12; 7) LTPs re-generated using Method B [18], discovered via merging 100 candidates based on texture similarity and local spatial co-occurrence, and setting N LT P = 12 for the iterative merging. Fig. 5 shows that the two sets of 100 LTP models achieve overall best prediction accuracy, and that the newly discovered 12 sLTPs have the best performance among the 5 small (s)LTP sets. Difference of ICC values between the sLTPs and the 100 LTP models was not significant for PLE emphysema subtype.

F. Clinical Associations of sLTPs
To evaluate clinical association of sLTPs, we first compute Spearman's partial correlations between %sLT P k within both lungs and the seven clinical characteristics listed in III-A, on the full MESA COPD dataset, using two models: Model 1 adjusted for demographical factors (age, race, gender, height and weight), and Model 2 further adjusted for %emph −950 .
The results are reported in Fig. 6. Correlation values for MRC dyspnea scale, post 6MWT breathlessness and post 6MWT fatigue are flipped in the figure so that more negative correlation values always correspond to more severe symptoms.
Overall, we obtained 47 and 31 significant correlations with Models 1 and 2. The sLTPs 7 and 8 are associated with less severe symptoms (positive correlations), while the other sLTPs correlate with symptoms (negative correlations). In Model 1, all clinical variables show significant correlations with 2 to 11 sLTPs. Model 2 looses significant correlations for post 6MWT breathlessness, but preserves all, or almost all, significant correlations for FEV1, 6MWT total distance, dyspnea and post-6MWT oximetry. With further adjustment for FEV1 in Model 2, sLTP 3 remains significantly correlated with baseline and post-6MWT pulse oximetry, sLTPs 2, 4 and 7 remain significantly correlated with 6MWT total distance, and sLTP 7 remains significantly correlated with MRC dyspnea scale.

IV. DISCUSSION & CONCLUSION
In this work, we propose a novel unsupervised learning framework for discovering emphysema-specific lung texture patterns and a small set of emphysema subtype candidates on the MESA COPD cohort of CT scans. The proposed method incorporates spatio-textural features via an original cost metric combining χ 2 -ℓ 2 constraints, along with data-driven parameter tuning, and Infomap graph partitioning.
Our methodological framework includes the introduction of a standardized spatial mapping of the lung shape utilizing Poisson distance map and conformal mapping to uniquely encode 3D voxel positions and enable comparison of CT scans. Our PDCM lung shape spatial mapping enables straightforward population-wide study of emphysema spatial patterns. By visualizing relative angular PDCM intensity projections on CLE-, PLE-and PSE-predominant subjects, we can see that regions of normal attenuation are absent for PLE-predominant subjects, which agrees with the definition of PLE (diffuse emphysema subtype). CLE-and PSE-predominant subjects appear to have emphysema regions concentrated in the superior part. This agrees with the observation made in [4] on the same dataset that CLE and PSE severity was greater in upper versus lower lung regions, whereas severity of PLE did not vary over the lung. By visualizing relative radial PDCM intensity projections, we can see that emphysema subjects show systematic lower attenuation values than subjects without emphysema, as expected. CLE-predominant subjects have more emphysema in the core part, whereas PSE-predominant subjects have more emphysema in the peel part. This agrees with the definitions of CLE and PSE. As a standardized tool, the proposed PDCM spatial mapping is not tied to emphysema patterns, and our future work will exploit such spatial mapping to study other pulmonary diseases.
With the proposed method, we discovered 12 spatially-informed lung texture patterns (sLTPs) in the MESA COPD Study. Qualitative visualization show that the discovered sLTPs appear to be textually homogeneous with specific average intensities and/or spatial prevalence. Using texton-based features to encode both texture and intensity is supported by [33] where "combination of both texture and densitometric measures strengthened the association with lung function" as we rely on association with physiological symptoms to evaluate our sLTPs. sLTPs (11,12) resemble CPFE studied in [34], where posterior emphysematous areas were more likely involved with interstitial lung abnormalities, which agrees with the posterior spatial prevalence seen in Fig. 3.
Extensive evaluations show that the discovered sLTPs are reproducible with respect to training sets, sampling of ROI for labeling, and certain scanner changes. The proposed incorporation of spatial and texture features obtains higher learning reproducibility compared to using texture features only, confirming the benefit of spatial regularization. The number of discovered sLTPs varies slightly between training subsets. This can be caused by a large change in the proportion of rare LTPs within these subsets, which modifies the weights in the Infomap similarity graph. A larger dataset with more diseased cases might be beneficial to solve such issue and would enable us to measure reproducibility on non-overlapping training subsets, which is a limitation of our study.
The sLTPs are able to encode the three standard emphysema subtypes, and thus link unsupervised discovery with clinical prior knowledge. Prediction quality is better than previous models, and close to the optimal level reached with 100 emphysema-specific LTPs.
While intra-cluster LTP homogeneity increases with the number of LTPs, hence leading to higher prediction performance, working with 100 LTPs leads to redundancy between subtypes which is detrimental when studying associations of individual LTPs with clinical measures. One-stage clustering leads to significantly lower prediction power for PLE and PSE subtypes compared to sLTPs, which demonstrate the benefit of the proposed two-stage learning.
Significant correlations with physiological symptoms were found for several measures.  [35] performed on a similar cohort size, but with highest COPD-prevalence, while reporting fewer significant positive correlations when proposing 7 radiological emphysema subtypes (called "factors") learned from 80 emphysema visual patterns. Correlations with standard emphysema subtypes, using similar models, were studied for the same population in [4]. Without adjusting for FEV1, CLE and PLE only showed significant associations with MRC dyspnea scale and 6MWT total distance, and only CLE showed significant associations with FEV1. With further adjustment for FEV1, only CLE and PLE showed significant associations with 6MWT total distance.
Progression patterns of the sLTPs will be investigated in the future, via sLTP labeling of longitudinal CT scans (with large time lapse). The sLTP histograms extracted in this study provide texture signatures that can be used to characterize and group CT scans. Patient grouping was found beneficial to study physiological indicators of COPD in [16], and will be considered in our future study. Further development is possible to improve the generation of image-based sLTPs with demographic and population-wide information, which would likely reveal population-specific and population-invariant patterns, but requiring a larger and more diseased cohort for training.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. posterior and medial vs. lateral positions, in addition to "peel to core" distance.       Intraclass correlation (ICC) and 95% confidence interval between predicted standard emphysema subtype scores and ground-truth. Differences with sLTP-based values are marked as ⋆ when significant (p < 0.05).  Black-boxes indicate statistically significant values (p < 0.05).  = 40, targeting 10 textons per standard emphysema subtype and normal tissue class, according to [13] Texton size 3×3×3 pixels, according to [18] # of lung sub-regions (for spatial features) = 36, according to binning of (r, θ, ϕ) in Section II-C3 N L T P : # of LTPs in initial set = 100, as suggested in [18], for sufficient diversity of the patterns and being able to discover rare emphysema types IEEE Trans Med Imaging. Author manuscript; available in PMC 2021 December 29.