Deep Learning COVID-19 Features on CXR using Limited Training Data Sets

Our segmentation network aims to extract lung and heart contour from the chest radiography images. We adopted an extended fully convolutional (FC)-DenseNet103 to perform semantic segmentation [15]. The training objective is

where $\mathcal{L}(\Theta)$ is the cross entropy loss of multi-categorical semantic segmentation and $\Theta$ denotes the network parameter set, which is composed of filter kernel weights and biases. Specifically, $\mathcal{L}(\Theta)$ is defined as

where $\mathds{1}(\cdot)$ is the indicator function, $p_{\Theta}(x_{j})$ denotes the softmax probability of the $j$-th pixel in a CXR image ${\bm{x}}$, and $y_{j}$ denotes the corresponding ground truth label. $s$ denotes class category, i.e., $s\in$ {background, heart, left lung, right lung}. $\lambda_{s}$ denotes weights given to each class category.

CXR images from different dataset resources may induce heterogeneity in their bits depth, compression type, image size, acquisition condition, scanning protocol, postprocessing, etc. Therefore, we develop a universal preprocessing step for data normalization to ensure uniform intensity histogram throughout the entire dataset. The detailed preprocessing steps are as follows:

1. 1.

   Data type casting (from uint8/uint16 to float32)

2. 2.

   Histogram equalization (gray level = [0, 255.0])

3. 3.

   Gamma correction ($\gamma$ = 0.5)

4. 4.

   Image resize (height, width = [256, 256])

Using the preprocessed data, we trained FC-DenseNet103[15] as our backbone segmentation network architecture. Network parameters were initialized by random distribution. We applied Adam optimizer [16] with an initial learning rate of 0.0001. Whenever training loss did not improve by certain criterion, the learning rate was reduced by factor 10. We adopted early stopping strategy based on validation performance. Batch size was optimized to 2. We implemented the network using PyTorch library [17].

The classification network aims to classify the chest X-ray images according to the types of disease. We adopted the relatively simple ResNet-18 as the backbone of our classification algorithm for two reasons. The first is to prevent from overfitting, since it is known that overfitting can occur when using an overly complex model for small number of data. Secondly, we intended to do transfer learning with pre-trained weights from ImageNet to compensate for the small training data set. We found that these strategy make the training stable even when the dataset size is small.

The labels were divided into four classes: normal, bacterial pneumonia, tuberculosis (TB), and viral pneumonia which includes the pneumonia caused by COVID-19 infection. We assigned the same class for viral pneumonia from other viruses (e.g. SARS-cov or MERS-cov) with COVID-19, since it is reported that they have similar radiologic features even challenging for the experienced radiologists [18]. Rather, we concentrated on more feasible work such as distinguishing bacterial pneumonia or tuberculosis from viral pneumonia, which show considerable differences in the radiologic features and are still useful for patient triage.

The pre-processed images were first masked with the lung masks from the segmentation networks, which are then fed into a classification network. Classification network were implemented in two different versions: global and local approaches. In the global approach, the masked images were resized to $224\times 224$, which were fed into the network. This approach is focusing on the global appearance of the CXR data, and was used as a baseline network for comparison. In fact, many of the existing researches employs similar procedure [6, 7, 8, 9].

In the local patch-based approach, which is our proposed method, the masked images were cropped randomly with a size of $224\times 224$, and resulting patches were used as the network inputs as shown in Fig. 1(b). In contrast to the global approach, various CXR images are resized to a much bigger $1024\times 1024$ image for our classification network to reflect the original pixel distribution better. Therefore, the segmentation mask from Fig. 1(a) are upsampled to match the $1024\times 1024$ image size. To avoid cropping the patch from the empty area of the masked image, the centers of patches were randomly selected within the lung areas. During the inference, $K$-number of patches were randomly acquired for each image to represent the entire attribute of the whole image. The number $K$ was chosen to sufficiently cover all lung pixels multiple times. Then, each patch was fed into the network to generate network output, and among $K$ network output the final decision was made based on majority voting, i.e. the most frequently declared class were regarded as final output as depicted in Fig. 1(b). In this experiments, the number of random patches $K$ was set to 100, which means that 100 patches were generated randomly from one whole image for majority voting.

For network training, pre-trained parameters from ImageNet are used for network weight initialization, after which the network was trained using the CXR data. As for optimization algorithm, Adam optimizer [16] with learning rate of 0.00001 was applied. The network were trained for 100 epochs, but we adopted early stopping strategy based on validation performance metrics. The batch size of 16 was used. We applied weight decay and $L_{1}$ regularization to prevent overfitting problem. The classification network was also implemented by Pytorch library.

We investigate the interpretability of our approach by visualizing a saliency map. One of the most widely used saliency map visualization methods is so-called gradient weighted class activation map (Grad-CAM) [14]. Specifically, the Grad-CAM saliency map of the class $c$ for a given input image ${\bm{x}}\in{\mathbb{R}}^{m\times n}$ is defined by

where ${\bm{f}}^{k}({\bm{x}})\in{\mathbb{R}}^{u\times v}$ is the $k$-th feature channel at the last convolution layer (which corresponds to the layer 4 of ResNet-18 in our case), $\textsc{Up}(\cdot)$ denotes the upsampling operator from a $u\times v$ feature map to the $m\times n$ image, $\sigma(\cdot)$ is the rectified linear unit (ReLU) [14]. Here, $\alpha_{k}^{c}$ is the feature weighted parameter for the class $c$, which can be obtained by

for some scaling parameter $Z$, where $y_{c}$ is the score for the class $c$ before the softmax layer and $f_{i}^{k}$ denotes the $i$-th pixel value of ${\bm{f}}^{k}({\bm{x}})$. The Grad-CAM map ${\bm{l}}^{c}$ is then normalized to have value in $[0,1]$. The Grad-CAM for the global approach is used as a baseline saliency map.

However, care should be taken in applying Grad-CAM to our local patch-based approach, since each patch has different score for the COVID-19 class. Therefore, to obtain the global saliency map, patch-wise Grad-CAM saliency maps should be weighted with the estimated probability of the class, and their average value should be computed. More specifically, our probabilistic Grad-CAM with respect to the input image ${\bm{x}}\in{\mathbb{R}}^{m\times n}$ has the following value at the $i$-th pixel location:

where ${\bm{x}}_{k}\in{\mathbb{R}}^{p\times q}$ denotes the $k$-th input patch, ${{\mathcal{Q}}}_{k}:{\mathbb{R}}^{p\times q}\mapsto{\mathbb{R}}^{m\times n}$ refers to the operator that copies the $p\times q$-size $k$-th patch into the appropriate location of a zero-padded image in ${\mathbb{R}}^{m\times n}$, and ${\bm{l}}^{c}({\bm{x}}_{k})\in{\mathbb{R}}^{p\times q}$ denotes the Grad-CAM computed by (3) with respect to the input patch ${\bm{x}}_{k}\in{\mathbb{R}}^{p\times q}$, and $K_{i}$ denotes the number of the frequency of the $i$-th pixel in the total $K$ patches. Additionally, the class probability $r^{c}({\bm{x}}_{k})$ for the $k$-th patch can be readily calculated after the soft-max layer. Accordingly, the average probability of each pixel belonging to a given class can be taken into consideration in Eq. (5) when constructing a global saliency map.

We used public CXR datasets, whose characteristics are summarized in Table I and Table II. In particular, the data in Table I are used for training and validation of the segmentation networks, since the ground-truth segmentation masks are available. The curated data in Table II are from some of the data in Table I as well as other COVID-19 resources, which were used for training, validation, and test for the classification network. More detailed descriptions of the dataset are follows.

The following standard biomarkers from CXR image analysis are investigated.

• •

  Lung morphology: Morphological structures of the segmented lung area as illustrated in Fig. 2(b) was evaluated throughout different classes.

• •

  Mean lung intensity: From the segmented lung area, we calculated mean value of the pixel intensity within the lung area as shown in Fig. 2(c).

• •

  Standard deviation of lung intensity: From the intensity histogram of lung area pixels, we calculated one standard deviation which is indicated as the black double-headed arrow in Fig. 2(c).

• •

  Cardiothoracic Ratio (CTR): CTR can be calculated by dividing the maximal transverse cardiac diameter by the maximal internal thoracic diameter annotated repectively as red and blue double-headed arrows in Fig. 2(a). Cardiothoracic Ratio (CTR) is a widely used marker to diagnosis cardiomegaly [25, 26]. We hypothesized that if cardiothoracic boundary become blurred by rounded opacities or consolidation in COVID-19 CXR [3, 2, 4], distinct off-average CTR value can be utilized as an abnormality alarm.

Statistical analysis for the potential biomarkers was performed using MATLAB 2015a (Mathworks, Natick). Kolmogorov Smirnov test was used to evaluate the normal distribution of marker candidates. For non-normally distributed variables, Wilcoxon signed rank test was used to compare segmentation performance with identical data size, and Wilcoxon rank sum test was used to compare COVID-19 marker candidates to those of other classes with different data sizes. Statistical significance (SS) levels were indicated as asterisks; * for $p<$0.05, ** for $p<$0.01 and *** for $p<$0.001.

The performance of the classification methods was evaluated using the confusion matrix. From the confusion matrix, true positive (TP), true negative (TN), false positive (FP), and false negative (FN) values were obtained, and 5 metrics for performance evaluation were calculated as below:

1. 1.

   $\mathrm{Accuracy}=(TN+TP)/(TN+TP+FN+FP)$

2. 2.

   $\mathrm{Precision}=TP/(TP+FP)$

3. 3.

   $\mathrm{Recall}=TP/(TP+FN)$

4. 4.

   $\mathrm{F1\,score}=2(\mathrm{Precision}\times\mathrm{Recall})/(\mathrm{Precision}+\mathrm{Recall})$

5. 5.

   $\mathrm{Specificity}=TN/(TN+FP)$

Among these, the F1 score was used as the evaluation metric for early stopping. The overall metric scores of the algorithm were calculated by averaging each metric for multiple classes.

Segmentation performance of anatomical structure was evaluated using Jaccard similarity coefficient. Table V presents the Jaccard similarity coefficient of each contour on the validation dataset. The results confirmed our method provides comparable accuracy to previous works using the JSRT dataset and the NLM(MC) dataset [27, 28].

To evaluate segmentation performance on cross-database, we tested either original or preprocessed images of the NLM dataset as inputs. The result shows that our universal preprocessing step for data normalization contributes to the processing of cross-database with statistically significant improvement on segmentation accuracy (Jaccard similarity coefficients from 0.932 to 0.943, $p<$ 0.001). This result indicates that preprocessing is crucial factor to ensure segmentation performance in cross-database.

To analyze morphological characteristics in the segmented lung area, a representative CXR radiograph for each class was selected for visual evaluation. Lung contour of each class showed differentiable features and showed mild tendency. In normal and tuberculosis cases (the first and the second row of Fig. 3, respectively), overall lung and heart contour were well-segmented. In the bacterial case, however, the segmented lung area was deformed due to wide spread opacity of bacterial pneumonia as shown in the third row of Fig. 3, and both the right cardiac and thoracic borders were lost. In overall bacterial infection cases, similar findings were occasionally observed which caused degraded segmentation performance. This suggests that abnormal morphology of the segmentation masks may be a useful biomarker to differentiate the severe infections. In the fourth row of Fig. 3, viral infection caused bilateral consolidations [29], thus partial deformation of lung area was observed. In the COVID-19 case of the fifth row of Fig. 3, despite the bi-basal infiltrations [30], lung area was fully segmented. In overall cases of the viral and the COVID-19 classes, lung areas were either normally or partially-imcompletely segmented, so morphological features of the segmentation masks may not be sufficiently discriminatory markers for viral and COVID-19 classes. Based on these morphological findings in segmented lung area, we further investigated other potential COVID-19 biomarkers.

We hypothesized that CXR appearance influenced by consolidations or infiltration of COVID-19 may be reflected in intensity of the radiograph. Thus, intensity-related COVID-19 marker candidates were extracted and compared.

The classification performances of the proposed method are provided in Table X. The confusion matrices for the (a) global method and the (b) local patch-based method are shown in Fig. 6. The proposed local patch-based approach showed consistently better performance than global approach in all metrics. In particular, as depicted in Table XI, our method showed the sensitivity of $92.5\%$ for COVID-19 and viruses, which was acceptable performance as a screening method, considering the fact that the sensitivity of COVID-19 diagnosis by X-ray image is known to be $69\%$ even for clinical experts and that the current gold standard, RT-PCR, has sensitivity of $91\%$ [4]. Moreover, compared to the global approach, the sensitivity of other classes are significantly high, which confirms the efficacy of our method.

Fig. 7 and Fig. 8 illustrate the examples of visualization of saliency map. As shown in Fig. 7(a), the existing Grad-CAM method for global approach showed the limitation that it only focuses on the broad main lesion so that it cannot properly differentiate multifocal lesions within the image. On the other hand, with the probabilistic Grad-CAM, multifocal GGOs and consolidations were visualized effectively by our local patch-based approach as shown in Fig. 7(c), which was in consistent with the findings reported by clinical experts. In particular, when we compute the probabilistic Grad-CAM for the COVID-19 class using patient images from various classes (i.e., normal, bacterial, TB, and COVID-19), a noticeable activation map was observed only in the COVID-19 patient data set, whereas almost no activations were observed in patients with other diseases and conditions as shown in Fig. 8. These results strongly support our claim that the probabilistic Grad-CAM saliency map from our local patch-based approach is more intuitive and interpretable compared to the existing methods.

In the diagnosis of COVID-19, other diseases mimicking COVID-19 pneumonia should be differentiated, including community-acquired pneumonia such as streptococcus pneumonia, mycoplasma and chlamydia related pneumonia, and other coronavirus infections.

In radiological literature, most frequently observed distribution patterns of COVID-19 are bilateral involvement, peripheral distribution and ground-glass opacification (GGO) [13]. Wong et al [4] found that consolidation was the most common finding (30/64, 47%), followed by GGO (21/64, 33%). CXR abnormalities had a peripheral (26/64, 41%) and lower zone distribution (32/64, 50%) with bilateral involvement (32/64, 50%), whereas pleural effusion was uncommon (2/64, 3%).

Our statistical analysis of the intensity distribution clearly showed that the globally distributed localized intensity variation is a discriminatory factor for COVID-19 CXR images, which was also confirmed with our saliency map. This clearly confirmed that the proposed method clearly reflects the radiological findings.

In pandemic situation of infectious disease, the distribution of medical resources is a matter of the greatest importance. As COVID-19 is spreading rapidly and surpassing the capacity of medical system in many countries, it is necessary to make reasonable decision to distribute the limited resources based on the ‘triage’, which determine the needs and urgency for each patient. In general, the most common cause of community acquired pneumonia is bacterial infection [5]. Specifically, most studies reported that S. pneumoniae is the most frequent causative strain ($15-42\%$), after which H. influenza ($11-12\%$) and viral pneumonia follow as the second and the third most common causes of pneumonia, respectively [5]. In addition, depending on the geological region, substantial proportion of pneumonia may be diagnosed as tuberculosis (up to $10\%$) [5]. Summing these up, the proportion of bacterial pneumonia and tuberculosis is suspected to be still large even in this pandemic situation of COVID-19. In this respect, the disease such as bacterial pneumonia or tuberculosis as well as normal condition can be excluded primarily, to preserve limited medical resources such as RT-PCR or CT only for those who suspected to be infected with COVID-19.

The detailed triage workflow that utilizes the proposed algorithm is described in Fig. 9. Specifically, our neural network is trained to classify other viral and COVID-19 in the same class. This is not only because it is strongly correlated with the radiological findings [18], but also useful as a triage. More specifically, by excluding normal, bacterial pneumonia, and TB at the early stage, we can use RT-PCR or CT for only those patients classified as other virus and COVID-19 cases for final diagnosis. By doing this procedure, we can save limited medical resources such as RT-PCR or CT to those patients whose diagnosis by CXR is difficult even by radiologists.

In order to investigate the origin of the apparent advantages of using local patch-based training over the global approach, we investigate the training dynamics to investigate the presence of overfitting. This is especially important, given that the training data is limited due to the difficulty of systematic CXR data collection for COVID-19 cases under current public health emergency.

Fig. 10 shows the curves for accuracy and F1 score of (a) the global approach and (b) the proposed local patch-based approach for each epoch. Note that both approaches use the same number of weight parameters. Still, thanks to the increasing training data set from the random patch cropping across all image area, our local patch-based algorithm did not showed any sign of overfitting even with the small numbers of training data, while the global approach showed significant overfitting problem. This clearly indicates that with the limited data set the patch-based neural network training may be a promising direction.

Since the proposed patch-based neural network architecture is designed by considering limited training data set, we investigated any potential performance loss in comparison with other state-of-the art (SOTA) deep learning approach that has been developed without such consideration. Specifically, COVID-Net [6] is one of the most successful approaches in COVID-19 diagnosis, so we chose it as the SOTA method.

The comparison between our method and COVID-Net is shown in Table XII. With the same dataset, our method showed overall accuracy of 91.9 % which is comparable to that of 92.4 % for COVID-Net. Furthermore, our method provided significantly improved sensitivity to COVID-19 cases compared to the COVID-Net. In addition, it is also remarkable that our method uses only about 10% number of parameters (11.6 M) compared to that of COVID-Net (116.6 M), because the proposed algorithm is developed based on less complex network architecture without increasing the complexity of the model. This may bring the advantages not only in the aspect of computational time but also in the aspect of performance and stability with small-sized dataset.

We are aware that the current study has limitations due to the lack of well-curated pneumonia CXR dataset. Specifically, our CXR data set come from a single or at most two sources (see Table II). Moreover, publicly available COVID-19 dataset [23] are largely extracted from online publications, website, etc, so they are not collected in a rigorous manner.

To mitigate the issue of potential bias from the limitation of the database, we employed a universal preprocessing step for data normalization for the entire dataset as discussed before. We investigated the effects of our preprocessing step on cross-database generalization by investigating the COVID-19 dataset, which poses the most severe intra-dataset heterogeneity. As shown in Fig. 11(b), each original CXRs of the COVID-19 class showed highly-varying intensity characteristics among each segmented anatomies. Thanks to our preprocessing step, the mean pixel intensity distribution between lung and heart regions of the preprocessed COVID-19 dataset (see Fig. 11(c)) became similar to the normal class in Fig. 11(a). The problem of heterogeneity can be also mitigated as shown in the intensity histograms (see Fig. 11(d)-(f)). The results confirmed that the original COVID-19 data could be well preprocessed to have comparable intensity distribution to that of well-curated normal data.