Introduction

ImageCLEF (the image retrieval and analysis evaluation campaign of the Cross-Language Evaluation Forum, CLEF) has organized challenges on image classification and retrieval since 2003 [1]. Since 2004, a medical image analysis and retrieval task has been organized [2,3].

The ImageCLEF 2017 [4] challenge included a task based on tuberculosis CT (Computed Tomography) volumes, the ImageCLEF 2017 tuberculosis task [5]. In this task, a dataset of lung CT scans was provided and 2 subtasks were proposed. When tuberculosis affects the lungs, several visual patterns can be seen in a CT image. However, the final diagnosis usually required other analyses than only the images [6]. A preliminary visualization of the CT volumes available in the ImageCLEF task showed lighter regions forming patterns that could be characterized with an holistic description of the lung. In [7], a graph-model of the lungs capable to differentiate between pulmonary hypertension and pulmonary embolism patients is presented. Both diseases present somewhat similar visual defects in lung CT scans, with different shapes and distributions. The graph was based on dividing the lung into several regions and use these as nodes of a graph. The edges of the graph encoded the difference between HU distributions in the lung regions. However, Dual Energy CT (DECT) scans were used in the study and the HU distribution can be described with more detail than in a standard single energy CT. Preliminary results showed that only one energy level did not contain enough information about the HU distribution to differentiate between the pulmonary hypertension and embolism patients.

Following the same approach, more complex features are used in this work. The descriptors based on the HU distribution were replaced by 3D texture features. Moreover, a deep analysis of the edges of the graph are helpful to describe the tuberculosis patterns. Our hypothesis is that a holistic analysis of the relations between regional texture features is able to encode subtle differences between patients with different tuberculosis type and drug resistance.

The following section contains a brief overview of the subtasks and dataset of the ImageCLEF 2017 tuberculosis task. More detailed information on the task can be found in the overview article [5]. Section 3 explains the process of building the textural graph-model of the lungs and all the variations tested for this task in detail. The results obtained by this approach in both subtasks are shown in Section 4. Finally, Section 5 concludes our participation in this challenge.

Subtasks and Datasets

The ImageCLEF 2017 TB task proposed two subtasks: i) Multi-drug resistance (MDR) prediction, and ii) Tuberculosis type (TBT) detection. The MDR task is a 2-class problem and the TBT task contains 5 classes. For both subtasks volumetric chest CT images with different voxel sizes and automatic segmentations of the lungs were provided. No other lung segmentation was attempted in this work and the masks provided were used. These masks were obtained with the method described in [8].

The competition was divided into two phases. In the first phase, the organizers released for each subtask a set of patient CT volumes as training set with their lung masks and ground truth labels. In the second phase, the test set with its lung segmentations but no labels were provided. The evaluation on the test data was performed by the organizers after the scheduled deadline for all runs that were submitted in time. The number of CT volumes for each task and set are specified in Tables 1 and 2.

Methods

This section details the process for the feature vector extraction from the CT images provided by the ImageCLEF TB task. The same technique was applied for describing the patients of both subtasks. This technique consists of creating graph-models of the lungs, with nodes based on a geometrical atlas and weighted edges encoding dissimilarities between 3D texture descriptors of each atlas region.

Isometric Volumes

The approach is based on 3D texture features. These features require having isometric voxels. The first step of our approach is to make the 3D images and masks provided by the organizers isometric. After analyzing the multiple resolutions found in the dataset and inter-slice distances, we opted for a voxel size of 1 mm to capture a maximum of information.

Atlas of the Lung

To build a graph with a fixed structure over the lung physiology we first need a localization system. In this case we chose the atlas developed by Depeursinge et al. in [9]. This atlas is only based on the mask of the lungs and provide 36 geometric regions dividing the lungs as shown in Figure 1. It is not based on the lung lobes.

3D Texture Features

Two state-of-the-art 3D texture features were selected to describe the texture inside the lung. The first method is a histogram of gradients based on the Fourier transform HOG (FHOG) introduced in [10]. 28 3D directions are used for the histogram obtaining a 28-dimensional feature vector per image voxel (f H ∈ R 28 ).

The second approach is the locally-oriented 3D Riesz-wavelet transform introduced by Dicente et al. in [11]. The parameters that obtained the best results in this article were used in the approach. These are: 3rd-order Riesz transform, 4 scales and 1st-order alignment. This configuration results in 40-dimensional feature vectors for each image voxel. The feature vector for a single voxel is then 10-dimensional containing the energy of each filter along the 4 scales (f R ∈ R 10 ).

Textural Graph-Model of the Lungs

Nodes Using as a base the 36-region atlas, the texture information from the 3D textural features was embedded in a 36-node graph. With this fixed number of nodes, we defined several graphs varying the number of edges and their weights. The following notation is used: given a patient p ∈ P and its 36-region atlas A p of the lungs with regions {r 1 , . . . , r 36 } ∈ A p we let G p be the 36-node graph of the patient p with nodes {N 1 , . . . , N 36 }. Each node N a represents one region r a ∈ A p .

Edges Three graphs were defined.

-Graph Full : This is the fully connected 36-node graph. For every pair of nodes N a and N b with a = b there exists an undirected edge E a,b . The total number of edges in this case is 630 ( 36•35 2 ). -Graph 66 : Based on the region adjacency defined by the atlas A p , in this case, there exists edge E a,b between nodes N a and N b if regions r a and r b are 3D adjacent in the atlas. This graph contains 66 edges in total. -Graph 84 : The last graph has the same 66 edges as Graph 66. Moreover, it has 18 additional edges connecting each pair of nodes representing opposite regions inside the atlas.

The three graphs are shown in Figure 2.

Graph Full

Graph 66 Graph 84 Feature Vector of a Region Several feature vectors can be extracted from a region r. Section 3.3 introduced the features extracted for a single voxel, f H and f R . Given a region r a , we extracted the mean (µ a ) and standard deviation (σ a ) of the features inside the region, i.e.: µ a (f H ), σ a (f H ), µ a (f R ), and σ a (f R ).

Feature Vector of a Patient Finally, the feature vector w p of a patient p, is defined as the ordered concatenation of the weights w a,b ∈ G p . Depending on the graph used, this feature vector can be 630-, 66-, or 84-dimensional.

Feature Normalization The last step before using the feature vectors w p with the classifier is to normalize them along the set of training patients P TRN . Each component w p,i of a vector w p corresponds to the weight of a different edge in the graph. Then, the components of a feature vector can not be seen independently and the normalization is required to keep these relations. Thus, the normalization was done for all components simultaneously. Several normalizations were tested:

The first is linearly resizing the min and max values among all components between 0 and 1, i.e.:

[ min

p∈PTRN {w p }, max p∈PTRN {w p }] → [0, 1]

. The second normalization aims to remove outliers. It considers the feature vector components as elements of a Gaussian distribution, and it centers the data around 0 with a standard deviation of 1, i.e.:

[µ p∈PTRN (w p ) − σ p∈PTRN (w p ), µ p∈PTRN (w p ) + σ p∈PTRN (w p )] → [−1, 1]

. These normalizations were also applied to the set of test patients P TST using the limits computed on the training set.

Feature Concatenation Fixing a graph structure (Graph Full, Graph 66, or Graph 84 ) and a measure between features (corr, cos, euc, or sumNorm), several sets W d of feature vectors w d,p were defined when varying the features used to encode the 3D texture in the regions. The possible features are µ a (f H ), σ a (f H ), µ a (f R ), and σ a (f R ). After normalizing each set of vectors W d , concatenations of these descriptors were tested in order to better describe each patient. These were based on the underlying features. The combinations tested were:

-FHOG using the mean and the std (µ(f H ) and σ(f H )), -Riesz using the mean and the std (µ(f R ) and σ(f R )), and FHOG and Riesz using the mean, the std, and both (µ(f H ) and µ(f R ), σ(f H ) and σ(f R ), and µ(f H ) and µ(f R ) and σ(f H ) and σ(f R )),

Feature Space Reduction When using the graph Graph Full and feature concatenations, the feature space dimension was much larger than the number of patients. To avoid the known problems when using such large feature spaces, we tested two feature space reduction techniques. Both were applied in the training phase. In the first case, we considered the ground truth labels to be an step function, with values {0, 1} for the MDR subtask, and {1, . . . , 5} for the TB task. Then, we computed the correlation between each feature dimension and these step functions, and selected the feature dimensions that correlated best with them. The threshold was set up as the mean absolute correlation of the feature dimensions with the labels. The second technique is based on the standard deviation of each feature. Only the components with a standard deviation higher than the mean of the standard deviations were selected. Both techniques reduced the size of the feature space by 2 approximately.

Classification

Multi-class support vector machine (SVM) classifiers with RBF kernel were used for each run of both subtasks, particularly, 2-class for the MDR task and 5-class for the TBT task. Grid search over the RBF parameters cost C and gamma γ was applied. Since the data were normalized, both C and γ moved between [2 −10 , 2 −9 , . . . , 2 10 ]. The best C and γ combination for a run was set as the one with highest cross-validation accuracy in the training set of each subtask.

Tested Runs

The procedure explained in Section 3.4 resulted in 648 runs per subtask. Table 3 summarizes all possible options for each step using the same name coding as in the results tables. Submitted Runs A total of 10 runs could be submitted in the ImageCLEF 2017 TB task. 5 runs were submitted for each task. The 5 runs selected for each task were either one of the 5 best runs with respect to the cross-validation accuracy on the training set (Acc TRN ), or a late fusion of a few of them. The late fusion was executed using the probabilities that the SVM classifier returned and the mean probability of belonging to each class. The same procedure was applied in both tasks. Tables 4 and 5 show the identifier and run setup of the 5 best runs for each task respectively. This section shows the results obtained on the training set by the submitted runs (Acc TRN ), and the final performance in the competition (Acc TST ). The final ranking was based on the Area Under the ROC Curve (AUC) for the MDR task, and on the unweighted Cohen's Kappa coefficient (Kappa) for the TBT task. Table 6 shows the results for the MDR subtask ordered by the ranking provided by the task organizers. The run identifiers MDR TopBest3 and MDR TopBest5 were obtained by doing late fusion with the 3 and 5 best runs respectively (see Section 3.6). The results for the TBT task are shown in Table 7. Again, the run identifiers TBT TopBest3 and TBT TopBest5 correspond to the late fusion of the 3 and 5 best runs respectively. The best run of the competition is shown in the same table.

Conclusions

This work presents a new graph-model of the lung based on regional 3D texture features for describing lungs affected by tuberculosis. The participation in the ImageCLEF 2017 tuberculosis 2017 allows for an objective comparison between methods since the ground truth for the test set was never released. For the MDR task, our method participated with 5 runs and obtained the 1st, 2nd and 3rd place in the challenge. Also in the case of the TBT subtask 5 runs were submitted but the best rank obtained was 10. The results underline the difficulty of both tasks and the suitability of our approach for describing TB patients. However, the results also suggest some overfitting by our method when comparing the accuracies obtained for the training and test sets.