# PSSPNN: PatchShuffle Stochastic Pooling Neural Network for an Explainable Diagnosis of COVID-19 with Multiple-Way Data Augmentation

Authors: Shui-Hua Wang; Yin Zhang; Xiaochun Cheng; Xin Zhang; Yu-Dong Zhang
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

Aim. COVID-19 has caused large death tolls all over the world. Accurate diagnosis is of significant importance for early treatment. Methods. In this study, we proposed a novel PSSPNN model for classification between COVID-19, secondary pulmonary tuberculosis, community-captured pneumonia, and healthy subjects. PSSPNN entails five improvements: we first proposed the n-conv stochastic pooling module. Second, a novel stochastic pooling neural network was proposed. Third, PatchShuffle was introduced as a regularization term. Fourth, an improved multiple-way data augmentation was used. Fifth, Grad-CAM was utilized to interpret our AI model. Results. The 10 runs with random seed on the test set showed our algorithm achieved a microaveraged F1 score of 95.79%. Moreover, our method is better than nine state-of-the-art approaches. Conclusion. This proposed PSSPNN will help assist radiologists to make diagnosis more quickly and accurately on COVID-19 cases.


---

## Body

## 1. Introduction

COVID-19 is a type of disease caused by a new strain of coronavirus. “CO” means corona, “VI” virus, and “D” disease. Up to 22 December 2020, COVID-19 has caused more than 78.0 million confirmed cases and over 1.71 million deaths (US 326.5 k, Brazil 188.2 k, India 146.1 k, Mexico 119.4 k, Italy 69.8 k, and UK 68.3 k).To diagnose COVID-19, two methods exist: (i) real-time reverse transcriptase PCR with nasopharyngeal swab samples to test the existence of RNA fragments and (ii) chest imaging (CI) examines the evidence of COVID-19. The first type of rRT-PCR method needs to wait for a few days to get the results, while the second type of CI approach could get quick results within minutes. The CI approaches have several advantages compared to rRT-PCR. First, the swab may be contaminated [1, 2]. Second, CI can detect the lesions of lungs where “ground-glass opacity (GGO)” will be observed to distinguish COVID-19 from healthy subjects. Third, CI can provide an immediate result as soon as imaging is complete. Fourth, reports show that chest computed tomography (CCT), one type of CI approach, can detect 97% of COVID-19 infections [3].Currently, there are three types of CI approaches: chest X-ray (CXR) [4], chest CT (CCT), and chest ultrasound (CUS) [5]. Among all three types of approaches, CCT can provide the finest resolution than the other two CI approaches, allowing visualization of extremely small nodules in the lung. The additional advantage of CCT is that it can provide high-quality, three-dimensional chest data where radiologists can clearly view the COVID-19 lesions, which may be obscure in the other two CI approaches.However, manual labeling by human experts is tedious, labor-intensive, and time-consuming [6]. Besides, the labeling performances are easily affected by interexpert and intraexpert factors (e.g., emotion, tiredness, and lethargy). Moreover, the labeling throughputs of radiologists are not comparable with artificial intelligence (AI) models. For example, senior radiologists may diagnose one scanning within five minutes, but AI can analyze thousands of samples within one minute. Particularly, early lesions are small and similar to surrounding normal tissues, which make them more challenging to measure and hence can potentially be neglected by radiologists.Traditional AI methods have successfully been applied in many medical fields. For instance, Wang et al. [7] chose wavelet packet Tsallis entropy as a feature descriptor and employed a real-coded biogeography-based optimization (RCBO) classifier. Jiang and Zhang [8] proposed a 6-level convolutional neural network (6L-CNN) for therapy and rehabilitation. Their performances were improved by replacing the traditional rectified linear unit with a leaky rectified linear unit. Fulton et al. [9] used ResNet-50 (RN-50) to classify Alzheimer’s disease with and without imagery. The authors found that ResNet-50 models help identify AD patients. Guo and Du [10] utilized a ResNet-18 (RN-18) model to recognize thyroid ultrasound standard plane (TUSP), achieving a classification accuracy of 83.88%. The experiments verified the effectiveness of RN-18. The aforementioned four algorithms can be transferred to the multiclass classification task of COVID-19 diagnosis.On the COVID-19 datasets, several recent publications reported promising results. For example, Cohen et al. [11] proposed a COVID severity score network (CSSNet), which achieved a mean absolute error (MAE) of 1.14 on geographic extent score and an MAE of 0.78 on lung opacity score. Li et al. [12] developed a fully automatic model (COVNet) to detect COVID-19 using chest CT and evaluated its performance. Wang et al. [13] proposed a 3D deep convolutional neural network to detect COVID-19 (DeCovNet). Zhang et al. [14] proposed a seven-layer convolutional neural network for COVID-19 diagnosis (7L-CCD). Their performance achieved an accuracy of 94.03±0.80 for the binary classification task (COVID-19 against healthy subjects). Ko et al. [15] proposed a fast-track COVID-19 classification network (FCONet). For the sake of the page limit, the details of those methods are not described, but we shall compare our method with those state-of-the-art methods in the following sections. Wang et al. [16] presented a CCSHNet via transfer learning and discriminant correlation analysis.Our study’s inspiration is to improve recognition performances of COVID-19 infection in CCT images by developing a novel deep neural network, PSSPNN, short for PatchShuffle stochastic pooling neural network. Our contributions entail the following five angles:(i)
The “n-conv stochastic pooling module (NCSPM)” is proposed, which comprises n-times repetitions of convolution layers and batch normalization layers, followed by stochastic pooling(ii)
A novel “stochastic pooling neural network (SPNN)” is proposed, the structure of which is inspired by VGG-16(iii)
A more advanced neural network, PatchShuffle SPNN (PSSPNN), is proposed where PatchShuffle is introduced as the regularization term in the loss function of SPNN(iv)
An improved multiple-way data augmentation is utilized to help the network avoid overfitting(v)
Grad-CAM is used to show the explainable heatmap, which displays association with lung lesions


## 2. Dataset

This retrospective study was exempt by Institutional Review Board of local hospitals. Four types of CCT were used: (i) COVID-19-positive patients, (ii) community-acquired pneumonia (CAP), (iii) second pulmonary tuberculosis (SPT), and (iv) healthy control (HC). Three diseased classes (COVID-19, CAP, and SPT) were chosen since they are all infectious diseases of the chest regions. We intend to include the fifth category (chest tumors) in our future studies.For each subject,nk slices of CCT were chosen via a slice level selection (SLS) method. For the three diseased groups (COVID-19, CAP, and SPT represented as k=1,2,3), the slice displaying the largest number of lesions and size was chosen. For HC subjects (k=4), any slice within the 3D image was randomly chosen. The slice-to-subject ratio n¯ is defined as
(1)n¯k=NSkNPk,k=1,⋯,4,where NS stands for the number of slices via the SLS method and NP is the number of patients.In all, we enrolled 521 subjects and produced 1164 slice images using the SLS method. Table1 lists the demographics of the four-category subject cohort with the values of triplets n¯,NP,NS, where n¯ of the total set equals to 2.23.Table 1
Subjects and images of four categories.
Class indexClass namen¯NPNS1COVID-192.271252842CAP2.281232813SPT2.181342934HC2.20139306Total2.235211164Three experienced radiologists (two juniors:C1 and C2 and one senior: C3) were convened to curate all the images. Suppose xCCT means one CCT scan, Y stands for the labeling of each individual radiologist. The last labeling YF of the CCT scan xCCT is obtained by
(2)YFxCCT=YC1,xCCTYC1,xCCT==YC2,xCCTMAVYallxCCTotherwise,where MAV denotes majority voting and Yall the labeling of all radiologists, viz.,
(3)YallxCCT=YC1,xCCT,YC2,xCCT,YC3,xCCT.The above two equations mean that in situations of disagreement between the analyses of two junior radiologistsC1,C2, we consult a senior radiologist C3to reach a MAV-type consensus. Data is available on request due to privacy/ethical restrictions.


## 3. Methodology

Table2 gives the abbreviation and full meanings in this study. Section 3.1 shows the preprocessing procedure. Sections 3.2–3.5 offer four improvements. Finally, Section 3.6 gives the implementation, measure indicators, and explainable technology used in our method.Table 2
Abbreviation and full name.
AbbreviationFull nameAPAverage poolingCAPCommunity-acquired pneumoniaCCTChest computed tomographyCIChest imagingCXRChest X-rayCUSChest ultrasoundDPDDiscrete probability distributionFCLFully connected layerFMFeature mapFMSFeature map sizeHSHyperparameter sizeL2Pl2-norm poolingMAMicroaveragedMAVMajority votingMDAMultiple-way data augmentationMPMax poolingNCSPMn-conv stochastic pooling moduleNWLNumber of weighted layersSAPNSalt-and-pepper noiseSCStrided convolutionSLSSlice level selectionSPStochastic poolingSPNNStochastic pooling neural networkSPTSecond pulmonary tuberculosisTCMTest confusion matrix


### 3.1. Preprocessing

The original raw dataset containedV slice images vai,i=1,2,⋯,V. The size of each image was sizevai=1024×1024×3. Figure 1 presents the pipeline for preprocessing of this dataset.Figure 1
Illustration of preprocessing.First, the color CCT images of four classes were converted into grayscale by retaining the luminance channel and obtaining the grayscale data setVB:
(4)VB=OgrayVA=vb1,vb2,⋯,vbi,⋯vbV,where Ogray means the grayscale operation.In the second step, the histogram stretching (HS) was utilized to increase the contrast of all images. Take thei-th image vbi,i=1,2,⋯,V as an example; its minimum and maximum grayscale values vbli and vbhi were calculated as follows:
(5)vbli=minw=1WBminh=1HBminc=1CBvbi∣w,h,c,vbhi=maxw=1WBmaxh=1HBmaxc=1CBvbi∣w,h,c.Here, (w, h, c) means the index of width, height, and channel directions along image vbi, respectively. WB,HB,CB means the maximum values of width, height, and channel to the image set VB. The new histogram stretched image set VC was calculated as follows:
(6)VC=OHSVB=vci≝vbi−vblivbhi−vbli,where OHS stands for the HS operation.In the third step, cropping was performed to remove the checkup bed at the bottom area and eliminate the texts at the margin regions. Cropped datasetVD is yielded as
(7)VD=OcropVC,a1,a2,a3,a4=vd1,vd2,⋯,vdi,⋯,vdV,where Ocrop represents crop operation. Parameters a1,a2,a3,a4 mean pixels to be cropped from four directions of the top, bottom, left, and right, respectively (unit: pixel).In the fourth step, each image was downsampled to a size ofWE,HE, obtaining the resized image set VE as
(8)VE=ODSVD,WE,HE=ve1,ve2,⋯,vei,⋯veVwhere ODS:a↦b represents the downsampling (DS) procedure, in which b stands for the downsampled image of the raw image a.Figure2 displays example images of the four categories, in which three are diseased, and one is healthy. The original size of each image in VA is 1024×1024×3, and the final preprocessed image is 256×256×1.Figure 2
Samples of four categories (three diseased and one healthy).
(a)
COVID-19(b)
CAP(c)
SPT(d)
HC


### 3.2. Improvement I:n-conv Stochastic Pooling Module

First, stochastic pooling (SP) [17] was introduced. In the standard convolutional neural networks, pooling is an essential component after each convolution layer, which was applied to reduce the size of feature maps (FMs). SP was shown to give better performance than average pooling and max pooling in recent publications [18–21]. Recently, strided convolution (SC) is commonly used, which also can shrink the FMs [22, 23]. Nevertheless, SC could be considered a simple pooling method, which always outputs the region’s fixed-position value [24].Suppose we have a postconvolution FMF=fiji=1,⋯,M×P,j=1,⋯,N×Q. The FM can be separated into M×N blocks, in which every block has the size of P×Q. Now we focus on the block Bmn=bx,y,x=1,⋯,P,y=1,⋯,Q which stands for the m-th row and n-th column blocks.The strided convolution (SC) traverses the input activation map with the strides, which equals the size of the blockP,Q, so here its output is set to
(9)BmnSC=b1,1.The l2-norm pooling (L2P), average pooling (AP), and max pooling (MP) generate the l2-norm value, average value, and maximum value within the blockBmn, respectively.
(10)BmnL2P=∑x=1P∑y=1Qb2x,yP×Q,BmnAP=1P×Q∑x=1P∑y=1Qbx,y,BmnMP=maxx=1Pmaxy=1Qbx,y.The SP provides a solution to the shortcomings of AP and MP. The AP outputs the average, so it will downscale the largest value, where the important features may sit on. On the other hand, MP reserves the maximum value but worsens the overfitting problem. SP is a three-step procedure. First, it generates the probability map (PM) for each entry in the blockBmn.
(11)pmx,y=bx,y∑x=1P∑y=1Qbx,y,s.t.∑x=1P∑y=1Qpmx,y=1,where pmx,y stands for the PM value at pixel x,y.In the second step, create a random variabler that takes the discrete probability distribution (DPD) as
(12)Prr=1,1=pm1,1Prr=1,2=pm1,2⋯Prr=1,Q=pm1,Q,Prr=2,1=pm2,1Prr=2,2=pm2,2⋯Prr=2,Q=pm2,Q,⋯⋯⋯⋯Prr=P,1=pmP,1Prr=P,2=pmP,2⋯Prr=P,Q=pmP,Q,where Pr represents the probability.In the third step, a sampler0 is drawn from the random variable r, and the corresponding position posr0=xr0,yr0. Then, the output of SP is at location posr0, namely,
(13)BmnSP=bxr0,yr0.Figure3 presents the comparison of four different pooling techniques. The top left shows the raw FM in which the pooling will take place at a 3×3 kernel. If we take the top-right block (in an orange rectangle) as an example, the L2P outputs 4.96, while AP and MP output 4.08 and 9.3, respectively. For the SP method, it will first generate the PM and then sample a position based on the PM (see the red fonts), and thus, SP outputs the value of 6.Figure 3
Comparison of four different pooling techniques (red rectangle indicates the example discussed in Section3.2).A new “n-conv stochastic pooling module” (NCSPM) is proposed in this study based on the SP layer discussed in previous paragraphs. The NCSPM entailsn-repetitions of a conv layer and a batch normalization layer, followed by an SP layer. Figure 4 shows the schematic of the proposed NCSPM module. In this study, n=1∨2, since we experimented using n=3, but the performance using n=3 did not improve.Figure 4
Schematic of proposed NCSPM.


### 3.3. Improvement II: Stochastic Pooling Neural Network

The second improvement of this study is to propose a stochastic pooling neural network (SPNN), whose structure was inspired by VGG-16 [25]. In VGG-16, the network used small kernels instead of large kernels and always used 2×2 filters with a stride of 2 for pooling. In the end, VGG-16 has two fully connected layers (FCLs).This proposed SPNN will follow the same structure design of VGG-16 but using the NCSPM module to replace the convolution block in VGG-16. The details of SPNN are shown in Table3, where NWL means the number of weighted layers, HS is the hyperparameter size, and FMS is the feature map size.Table 3
Structure of proposed 11-layer SPNN.
IndexNameNWLHSFMS1Input256×256×12NCSPM-123×3,32×2128×128×323NCSPM-223×3,32×264×64×324NCSPM-323×3,64×232×32×645NCSPM-423×3,64×216×16×646NCSPM-513×3,128×18×8×1287Flatten1×1×81928FCL-11140×8192, 140×11×1×1409FCL-214×140, 4×11×1×410Softmax1×1×4Compared to ordinary CNN, the advantages of SPNN are two folds: (i) SPNN helps prevent overfitting; (ii) SPNN is parameter-free. (iii) SPNN can be easily combined with other advanced techniques, such as batch normalization and dropout. In total, we create this 11-layer SPNN. We have attempted to insert more NCSPMs or more FCLs, which does not show performance improvement but more computation burden. The structure of the proposed model is summarized in Table3. The c1×c1,c2×c3 related to NCSPM stands for c3 repetitions of c2 filters with size of c1×c1. For the FCL, the c1×c2,c3×c4 stands for a weight matrix is with size of c1×c2 and a bias matrix is with size of c3×c4. In the last column of Table 3, the format of c1×c2×c3 represents the feature map’s size in three dimensions: height, width, and channel. Directly using transfer learning is another alternative.In this study, we chose to create a custom neural network by designing its structure and training the whole network using our own data. The reason is some reports have shown this “built your own network from scratch” can achieve better performance than transfer learning [26, 27].


### 3.4. Improvement III: PatchShuffle SPNN

Kang et al. [28] presented a new PatchShuffle method. In each minibatch, images and feature maps undergo a transformation such that pixels with that patch are shuffled. By generating fake images/feature maps with interior order-less patches, PatchShuffle creates local variations and reduces the possibility of the AI model overfitting. Therefore, PatchShuffle is a beneficial supplement to various existing training regularization methods [28].Assume there is a matrixX of M×M entries. A random variable v controls whether the matrix X to be PatchShuffled or not. The random variable v obeys the Bernoulli distribution
(14)v~Bernoulliε.Namely,v=1 with probability ε, and v=0 with probability 1−ε. The resulted matrix X^ is written as
(15)X^=1−vX+vFPSX,where FPS is the PatchShuffle operation. Suppose the size of each patch is m×m, we can express the matrix X as
(16)X=x11x12⋯x1,M/mx21x22⋯x2,M/m⋮⋮⋱⋮xM/m,1xM/m,2⋯xM/m,M/m,where xij stands for a nonoverlapping patch at i-th row and j-th column. The PatchShuffle transformation works on every patch.
(17)FPSX=FPSxij,i=1,⋯,Mm,j=1,⋯,Mm.The shuffled patch is formulated as(18)FPSxij=eij×xij×eij′,where eij is the row permutation matrix and eij′ is the column permutation matrix [29]. In practice, a randomly shuffle operation is used to replace the row and column permutation operation. Each patch will undergo one of the m2! possible permutations.We proposed integrating PatchShuffle into our SPNN, and this new network model is named as PatchShuffle stochastic pooling neural network (PSSPNN). The PatchShuffle operation acts on both the input image layer (see grayscale image Figure5 and colorful image Figure 6 with different values of m) and the feature maps of all the convolutional layers (9 conv layers from NCSPM-1 to NCSPM-5).Figure 5
Illustration of PatchShuffle on a grayscale COVID-19 image.
(a)
Raw image(b)
m=2(c)
m=3(d)
m=4(e)
m=5(f)
m=6Figure 6
Illustration of PatchShuffle on a colorful house image.
(a)
Raw image(b)
m=2(c)
m=3(d)
m=4(e)
m=5(f)
m=6The schematics of building PSSPNN from the SPNN are drawn in Figure7, where either inputs or feature maps are randomly selected to undergo the PatchShuffle operation. To reach the best bias-variance trade-off, only a small percentage (ε) of the images/feature maps will undergo FPS operation.Figure 7
Build PSSPNN from SPNN.For simplicity, we consider the PatchShuffling images as an example, and the training loss functionℓ of the proposed PSSPNN is written as
(19)ℓPSSPNNX,y,W=1−vℓX,y,W+vℓFPSX,y,W,where ℓ stands for the ordinary loss function and ℓPSSPNN the loss function of PSSPNN. X represents the original images and FPSX the PatchShuffled images. The label is symbolized as y, and the weights are symbolized as W.Considering the extreme situations whenv=0∨1, we have
(20)ℓPSSPNNX,y,W=ℓX,y,Wv=0,ℓFPSX,y,Wv=1,which means the loss function degrades to ordinary loss function when v=0, and meanwhile, the loss function equals to training all images PatchShuffled when v=1. Taking mathematical expectation of v, equation (19) turns to
(21)11−εEvℓPSSPNNX,y,W=ℓX,y,W+ε1−εℓFPSX,y,W,where ε/1−εℓFPSX,y,W acts as a regularization term.


### 3.5. Improvement IV: Improved Multiple-Way Data Augmentation

This small four-category dataset makes our AI model prone to overfitting. In order to alleviate the overfitting and handle the low sample-size problem, the multiple-way data augmentation (MDA) [30] method was chosen and further improved. In the original 14-way MDA [30], the authors used seven different data augmentation (DA) techniques to the raw image and its horizontal image. Their seven DA techniques are as follows: noise injection, horizontal shear, vertical shear, rotation, Gamma correction, scaling, and translation.Figure8 shows a 16-way data augmentation method. The difference between the proposed 16-way DA with the traditional 14-way DA is that we add the salt-and-pepper noise (SAPN). Although the SAPN defies intuition as it never takes place in realistic CCT images, we found that it can increase performance. The same observation was reported by Li et al. [31], where the authors used salt and pepper noise for the identification of early esophageal cancer. Table 4 shows the pseudocode of this proposed 16-way improved data augmentation.
(22)vi↦Di=SviMviSk=1,..,8fkDAvi⏟WSk=1,..,8fkDAMvi⏟W,where W=30 in this study. We tested a greater value of W, but it does not bring about significant improvement. Hence, one image vi will generate Di=16W+2=482 images (including original image), as shown in Figure 8.Step 1.
Eight geometric or photometric DA transforms were utilized on raw imagevi, as shown in Figure 8. We use fkDA,k=1,⋯,8 to denote each DA operation. Note each DA operation fkDA will yield W new images. So, for a given image vi, we will produce an enhanced dataset Sk=1,..,8fkDAvi, where S stands for concatenation function.Step 2.
Horizontal mirror image is generated asMvi, where M means horizontal mirror function.Step 3.
The raw imagevi, the mirrored image Mvi, all the above 8-way DA results of raw image Sk=1,..,8fkDAvi, and 8-way DA results of horizontal mirrored image Sk=1,..,8fkDAMvi are combined. Mathematically, one training image vi will generate to a dataset Di, which contains 16W+1 new images.Figure 8
The proposed 16-way improved data augmentation method.Table 4
Pseudocode of proposed 16-way improved data augmentation.
InputRaw imageviStep 1Eight geometric or photometric DA transforms were utilized on raw imagevi.Step 2A horizontal mirror image is generated.Step 3The raw imagevi, the mirrored image, the above 8-way DA results of the raw image, and the 8-way DA results of the horizontal mirrored image are combined to form a dataset Di.OutputEnhanced dataset of raw imageDiTaking Figure2(a) as an example raw image, Figure 9 shows the 8-way DA results, i.e., fkDAvi,k=1,⋯,8. Due to the page limit, the mirror image and its corresponding 8-way DA results are not shown here.Figure 9
Proposed 16-way data augmentation results.
(a)
Gaussian noise(b)
Salt-and-pepper noise(c)
Horizontal shear(d)
Vertical shear(e)
Rotation(f)
Gamma correction(g)
Random translation(h)
Scaling


### 3.6. Implementation, Measure, and Explainability

Table5 itemizes the nontest and test sets for each category. The whole dataset V contains four nonoverlapping categories V=V1,V2,V3,V4. For each category, the dataset will be split into nontest set and test set Vk→Vkntest,Vktest,k=1,2,3,4.
(23)V=V1V2V3V4=VntestVtest=V1ntestV1testV2ntestV2testV3ntestV3testV4ntestV4test.Table 5
Dataset splitting.
Nontest (10-fold cross-validation)Test (10 runs)TotalCOVID-19V1ntest=227V1test=57V1=284CAPV2ntest=225V2test=56V2=281SPTV3ntest=234V3test=59V3=293HCV4ntest=245V4test=61V4=306Our experiment entails two phases. At phase I, 10-fold cross-validation was used for validation on the nontest set to select the best hyperparameters and best network structure. The 16-way DA was used on the training set of 10-fold cross-validation. The hyperparameter of the proposed PSSPNN was determined over the nontest setVntest. Afterward at phase II, we train our model using the nontest set Vntest 10 times with different initial seeds and attain the test results over the test set Vtest. After combining the Rtest runs, we attain a summation of the test confusion matrix (TCM) Dtest. Table 5 shows the dataset splitting, where x stands for the number of elements in the dataset x.The ideal TCM is a diagonal matrix with the form of(24)Didealtest=Rtest×V1test0000V2test0000V3test0000V4test,in which all the off-diagonal elements are zero, meaning no prediction errors. In realistic scenarios, the AI model will make errors, and the performance is calculated per category. For each class z=1,⋯,4, we set the label of that class as positive, and the labels of all the rest classes as negative. Three performance metrics (sensitivity, precision, and F1 score) per category are defined below:
(25)Senz=TPzTPz+FNz,Prcz=TPzTPz+FPz,F1z=2∗Prcz∗SenzPrcz+Senz.The performance can be measured over all four categories. The microaveraged (MA) F1 (symbolized asF1μ) is used since our dataset is slightly unbalanced:
(26)F1μ=2∗Prcμ∗SenμPrcμ+Senμ,where
(27)Senμ=∑z=14TPz∑z=14TPz+FNz,Prcμ=∑z=14TPz∑z=14TPz+FPz.Finally, gradient-weighted class activation mapping (Grad-CAM) [32] was employed to provide explanations on how our model makes the decision. It exploits the gradient of the categorization score regarding the convolutional features decided by the deep model to visualize the regions of the image that are the most vital for the image classification task [33]. The output of NCSPM-5 in Table 3 was used for Grad-CAM.


## 3.1. Preprocessing

The original raw dataset containedV slice images vai,i=1,2,⋯,V. The size of each image was sizevai=1024×1024×3. Figure 1 presents the pipeline for preprocessing of this dataset.Figure 1
Illustration of preprocessing.First, the color CCT images of four classes were converted into grayscale by retaining the luminance channel and obtaining the grayscale data setVB:
(4)VB=OgrayVA=vb1,vb2,⋯,vbi,⋯vbV,where Ogray means the grayscale operation.In the second step, the histogram stretching (HS) was utilized to increase the contrast of all images. Take thei-th image vbi,i=1,2,⋯,V as an example; its minimum and maximum grayscale values vbli and vbhi were calculated as follows:
(5)vbli=minw=1WBminh=1HBminc=1CBvbi∣w,h,c,vbhi=maxw=1WBmaxh=1HBmaxc=1CBvbi∣w,h,c.Here, (w, h, c) means the index of width, height, and channel directions along image vbi, respectively. WB,HB,CB means the maximum values of width, height, and channel to the image set VB. The new histogram stretched image set VC was calculated as follows:
(6)VC=OHSVB=vci≝vbi−vblivbhi−vbli,where OHS stands for the HS operation.In the third step, cropping was performed to remove the checkup bed at the bottom area and eliminate the texts at the margin regions. Cropped datasetVD is yielded as
(7)VD=OcropVC,a1,a2,a3,a4=vd1,vd2,⋯,vdi,⋯,vdV,where Ocrop represents crop operation. Parameters a1,a2,a3,a4 mean pixels to be cropped from four directions of the top, bottom, left, and right, respectively (unit: pixel).In the fourth step, each image was downsampled to a size ofWE,HE, obtaining the resized image set VE as
(8)VE=ODSVD,WE,HE=ve1,ve2,⋯,vei,⋯veVwhere ODS:a↦b represents the downsampling (DS) procedure, in which b stands for the downsampled image of the raw image a.Figure2 displays example images of the four categories, in which three are diseased, and one is healthy. The original size of each image in VA is 1024×1024×3, and the final preprocessed image is 256×256×1.Figure 2
Samples of four categories (three diseased and one healthy).
(a)
COVID-19(b)
CAP(c)
SPT(d)
HC


## 3.2. Improvement I:n-conv Stochastic Pooling Module

First, stochastic pooling (SP) [17] was introduced. In the standard convolutional neural networks, pooling is an essential component after each convolution layer, which was applied to reduce the size of feature maps (FMs). SP was shown to give better performance than average pooling and max pooling in recent publications [18–21]. Recently, strided convolution (SC) is commonly used, which also can shrink the FMs [22, 23]. Nevertheless, SC could be considered a simple pooling method, which always outputs the region’s fixed-position value [24].Suppose we have a postconvolution FMF=fiji=1,⋯,M×P,j=1,⋯,N×Q. The FM can be separated into M×N blocks, in which every block has the size of P×Q. Now we focus on the block Bmn=bx,y,x=1,⋯,P,y=1,⋯,Q which stands for the m-th row and n-th column blocks.The strided convolution (SC) traverses the input activation map with the strides, which equals the size of the blockP,Q, so here its output is set to
(9)BmnSC=b1,1.The l2-norm pooling (L2P), average pooling (AP), and max pooling (MP) generate the l2-norm value, average value, and maximum value within the blockBmn, respectively.
(10)BmnL2P=∑x=1P∑y=1Qb2x,yP×Q,BmnAP=1P×Q∑x=1P∑y=1Qbx,y,BmnMP=maxx=1Pmaxy=1Qbx,y.The SP provides a solution to the shortcomings of AP and MP. The AP outputs the average, so it will downscale the largest value, where the important features may sit on. On the other hand, MP reserves the maximum value but worsens the overfitting problem. SP is a three-step procedure. First, it generates the probability map (PM) for each entry in the blockBmn.
(11)pmx,y=bx,y∑x=1P∑y=1Qbx,y,s.t.∑x=1P∑y=1Qpmx,y=1,where pmx,y stands for the PM value at pixel x,y.In the second step, create a random variabler that takes the discrete probability distribution (DPD) as
(12)Prr=1,1=pm1,1Prr=1,2=pm1,2⋯Prr=1,Q=pm1,Q,Prr=2,1=pm2,1Prr=2,2=pm2,2⋯Prr=2,Q=pm2,Q,⋯⋯⋯⋯Prr=P,1=pmP,1Prr=P,2=pmP,2⋯Prr=P,Q=pmP,Q,where Pr represents the probability.In the third step, a sampler0 is drawn from the random variable r, and the corresponding position posr0=xr0,yr0. Then, the output of SP is at location posr0, namely,
(13)BmnSP=bxr0,yr0.Figure3 presents the comparison of four different pooling techniques. The top left shows the raw FM in which the pooling will take place at a 3×3 kernel. If we take the top-right block (in an orange rectangle) as an example, the L2P outputs 4.96, while AP and MP output 4.08 and 9.3, respectively. For the SP method, it will first generate the PM and then sample a position based on the PM (see the red fonts), and thus, SP outputs the value of 6.Figure 3
Comparison of four different pooling techniques (red rectangle indicates the example discussed in Section3.2).A new “n-conv stochastic pooling module” (NCSPM) is proposed in this study based on the SP layer discussed in previous paragraphs. The NCSPM entailsn-repetitions of a conv layer and a batch normalization layer, followed by an SP layer. Figure 4 shows the schematic of the proposed NCSPM module. In this study, n=1∨2, since we experimented using n=3, but the performance using n=3 did not improve.Figure 4
Schematic of proposed NCSPM.


## 3.3. Improvement II: Stochastic Pooling Neural Network

The second improvement of this study is to propose a stochastic pooling neural network (SPNN), whose structure was inspired by VGG-16 [25]. In VGG-16, the network used small kernels instead of large kernels and always used 2×2 filters with a stride of 2 for pooling. In the end, VGG-16 has two fully connected layers (FCLs).This proposed SPNN will follow the same structure design of VGG-16 but using the NCSPM module to replace the convolution block in VGG-16. The details of SPNN are shown in Table3, where NWL means the number of weighted layers, HS is the hyperparameter size, and FMS is the feature map size.Table 3
Structure of proposed 11-layer SPNN.
IndexNameNWLHSFMS1Input256×256×12NCSPM-123×3,32×2128×128×323NCSPM-223×3,32×264×64×324NCSPM-323×3,64×232×32×645NCSPM-423×3,64×216×16×646NCSPM-513×3,128×18×8×1287Flatten1×1×81928FCL-11140×8192, 140×11×1×1409FCL-214×140, 4×11×1×410Softmax1×1×4Compared to ordinary CNN, the advantages of SPNN are two folds: (i) SPNN helps prevent overfitting; (ii) SPNN is parameter-free. (iii) SPNN can be easily combined with other advanced techniques, such as batch normalization and dropout. In total, we create this 11-layer SPNN. We have attempted to insert more NCSPMs or more FCLs, which does not show performance improvement but more computation burden. The structure of the proposed model is summarized in Table3. The c1×c1,c2×c3 related to NCSPM stands for c3 repetitions of c2 filters with size of c1×c1. For the FCL, the c1×c2,c3×c4 stands for a weight matrix is with size of c1×c2 and a bias matrix is with size of c3×c4. In the last column of Table 3, the format of c1×c2×c3 represents the feature map’s size in three dimensions: height, width, and channel. Directly using transfer learning is another alternative.In this study, we chose to create a custom neural network by designing its structure and training the whole network using our own data. The reason is some reports have shown this “built your own network from scratch” can achieve better performance than transfer learning [26, 27].


## 3.4. Improvement III: PatchShuffle SPNN

Kang et al. [28] presented a new PatchShuffle method. In each minibatch, images and feature maps undergo a transformation such that pixels with that patch are shuffled. By generating fake images/feature maps with interior order-less patches, PatchShuffle creates local variations and reduces the possibility of the AI model overfitting. Therefore, PatchShuffle is a beneficial supplement to various existing training regularization methods [28].Assume there is a matrixX of M×M entries. A random variable v controls whether the matrix X to be PatchShuffled or not. The random variable v obeys the Bernoulli distribution
(14)v~Bernoulliε.Namely,v=1 with probability ε, and v=0 with probability 1−ε. The resulted matrix X^ is written as
(15)X^=1−vX+vFPSX,where FPS is the PatchShuffle operation. Suppose the size of each patch is m×m, we can express the matrix X as
(16)X=x11x12⋯x1,M/mx21x22⋯x2,M/m⋮⋮⋱⋮xM/m,1xM/m,2⋯xM/m,M/m,where xij stands for a nonoverlapping patch at i-th row and j-th column. The PatchShuffle transformation works on every patch.
(17)FPSX=FPSxij,i=1,⋯,Mm,j=1,⋯,Mm.The shuffled patch is formulated as(18)FPSxij=eij×xij×eij′,where eij is the row permutation matrix and eij′ is the column permutation matrix [29]. In practice, a randomly shuffle operation is used to replace the row and column permutation operation. Each patch will undergo one of the m2! possible permutations.We proposed integrating PatchShuffle into our SPNN, and this new network model is named as PatchShuffle stochastic pooling neural network (PSSPNN). The PatchShuffle operation acts on both the input image layer (see grayscale image Figure5 and colorful image Figure 6 with different values of m) and the feature maps of all the convolutional layers (9 conv layers from NCSPM-1 to NCSPM-5).Figure 5
Illustration of PatchShuffle on a grayscale COVID-19 image.
(a)
Raw image(b)
m=2(c)
m=3(d)
m=4(e)
m=5(f)
m=6Figure 6
Illustration of PatchShuffle on a colorful house image.
(a)
Raw image(b)
m=2(c)
m=3(d)
m=4(e)
m=5(f)
m=6The schematics of building PSSPNN from the SPNN are drawn in Figure7, where either inputs or feature maps are randomly selected to undergo the PatchShuffle operation. To reach the best bias-variance trade-off, only a small percentage (ε) of the images/feature maps will undergo FPS operation.Figure 7
Build PSSPNN from SPNN.For simplicity, we consider the PatchShuffling images as an example, and the training loss functionℓ of the proposed PSSPNN is written as
(19)ℓPSSPNNX,y,W=1−vℓX,y,W+vℓFPSX,y,W,where ℓ stands for the ordinary loss function and ℓPSSPNN the loss function of PSSPNN. X represents the original images and FPSX the PatchShuffled images. The label is symbolized as y, and the weights are symbolized as W.Considering the extreme situations whenv=0∨1, we have
(20)ℓPSSPNNX,y,W=ℓX,y,Wv=0,ℓFPSX,y,Wv=1,which means the loss function degrades to ordinary loss function when v=0, and meanwhile, the loss function equals to training all images PatchShuffled when v=1. Taking mathematical expectation of v, equation (19) turns to
(21)11−εEvℓPSSPNNX,y,W=ℓX,y,W+ε1−εℓFPSX,y,W,where ε/1−εℓFPSX,y,W acts as a regularization term.


## 3.5. Improvement IV: Improved Multiple-Way Data Augmentation

This small four-category dataset makes our AI model prone to overfitting. In order to alleviate the overfitting and handle the low sample-size problem, the multiple-way data augmentation (MDA) [30] method was chosen and further improved. In the original 14-way MDA [30], the authors used seven different data augmentation (DA) techniques to the raw image and its horizontal image. Their seven DA techniques are as follows: noise injection, horizontal shear, vertical shear, rotation, Gamma correction, scaling, and translation.Figure8 shows a 16-way data augmentation method. The difference between the proposed 16-way DA with the traditional 14-way DA is that we add the salt-and-pepper noise (SAPN). Although the SAPN defies intuition as it never takes place in realistic CCT images, we found that it can increase performance. The same observation was reported by Li et al. [31], where the authors used salt and pepper noise for the identification of early esophageal cancer. Table 4 shows the pseudocode of this proposed 16-way improved data augmentation.
(22)vi↦Di=SviMviSk=1,..,8fkDAvi⏟WSk=1,..,8fkDAMvi⏟W,where W=30 in this study. We tested a greater value of W, but it does not bring about significant improvement. Hence, one image vi will generate Di=16W+2=482 images (including original image), as shown in Figure 8.Step 1.
Eight geometric or photometric DA transforms were utilized on raw imagevi, as shown in Figure 8. We use fkDA,k=1,⋯,8 to denote each DA operation. Note each DA operation fkDA will yield W new images. So, for a given image vi, we will produce an enhanced dataset Sk=1,..,8fkDAvi, where S stands for concatenation function.Step 2.
Horizontal mirror image is generated asMvi, where M means horizontal mirror function.Step 3.
The raw imagevi, the mirrored image Mvi, all the above 8-way DA results of raw image Sk=1,..,8fkDAvi, and 8-way DA results of horizontal mirrored image Sk=1,..,8fkDAMvi are combined. Mathematically, one training image vi will generate to a dataset Di, which contains 16W+1 new images.Figure 8
The proposed 16-way improved data augmentation method.Table 4
Pseudocode of proposed 16-way improved data augmentation.
InputRaw imageviStep 1Eight geometric or photometric DA transforms were utilized on raw imagevi.Step 2A horizontal mirror image is generated.Step 3The raw imagevi, the mirrored image, the above 8-way DA results of the raw image, and the 8-way DA results of the horizontal mirrored image are combined to form a dataset Di.OutputEnhanced dataset of raw imageDiTaking Figure2(a) as an example raw image, Figure 9 shows the 8-way DA results, i.e., fkDAvi,k=1,⋯,8. Due to the page limit, the mirror image and its corresponding 8-way DA results are not shown here.Figure 9
Proposed 16-way data augmentation results.
(a)
Gaussian noise(b)
Salt-and-pepper noise(c)
Horizontal shear(d)
Vertical shear(e)
Rotation(f)
Gamma correction(g)
Random translation(h)
Scaling


## 3.6. Implementation, Measure, and Explainability

Table5 itemizes the nontest and test sets for each category. The whole dataset V contains four nonoverlapping categories V=V1,V2,V3,V4. For each category, the dataset will be split into nontest set and test set Vk→Vkntest,Vktest,k=1,2,3,4.
(23)V=V1V2V3V4=VntestVtest=V1ntestV1testV2ntestV2testV3ntestV3testV4ntestV4test.Table 5
Dataset splitting.
Nontest (10-fold cross-validation)Test (10 runs)TotalCOVID-19V1ntest=227V1test=57V1=284CAPV2ntest=225V2test=56V2=281SPTV3ntest=234V3test=59V3=293HCV4ntest=245V4test=61V4=306Our experiment entails two phases. At phase I, 10-fold cross-validation was used for validation on the nontest set to select the best hyperparameters and best network structure. The 16-way DA was used on the training set of 10-fold cross-validation. The hyperparameter of the proposed PSSPNN was determined over the nontest setVntest. Afterward at phase II, we train our model using the nontest set Vntest 10 times with different initial seeds and attain the test results over the test set Vtest. After combining the Rtest runs, we attain a summation of the test confusion matrix (TCM) Dtest. Table 5 shows the dataset splitting, where x stands for the number of elements in the dataset x.The ideal TCM is a diagonal matrix with the form of(24)Didealtest=Rtest×V1test0000V2test0000V3test0000V4test,in which all the off-diagonal elements are zero, meaning no prediction errors. In realistic scenarios, the AI model will make errors, and the performance is calculated per category. For each class z=1,⋯,4, we set the label of that class as positive, and the labels of all the rest classes as negative. Three performance metrics (sensitivity, precision, and F1 score) per category are defined below:
(25)Senz=TPzTPz+FNz,Prcz=TPzTPz+FPz,F1z=2∗Prcz∗SenzPrcz+Senz.The performance can be measured over all four categories. The microaveraged (MA) F1 (symbolized asF1μ) is used since our dataset is slightly unbalanced:
(26)F1μ=2∗Prcμ∗SenμPrcμ+Senμ,where
(27)Senμ=∑z=14TPz∑z=14TPz+FNz,Prcμ=∑z=14TPz∑z=14TPz+FPz.Finally, gradient-weighted class activation mapping (Grad-CAM) [32] was employed to provide explanations on how our model makes the decision. It exploits the gradient of the categorization score regarding the convolutional features decided by the deep model to visualize the regions of the image that are the most vital for the image classification task [33]. The output of NCSPM-5 in Table 3 was used for Grad-CAM.


## 4. Experiments, Results, and Discussions

The experiment was carried out on the programming platform of Matlab 2020b. The programs ran on Windows 10 with 16GB RAM and 2.20GHz Intel Core i7-8750H CPU. The performances are reported over the test set with 10 runs.


### 4.1. Comparison of SPNN and Other Pooling Methods

In the first experiment, we compared the proposed SPNN with four standard CNNs with different pooling methods. The first CNN uses strided convolution in five modules to replace the stochastic pooling. The second to fourth comparison CNN models use L2P, AP, and MP, respectively. Those four baseline methods are called SC-CNN, L2P-CNN, AP-CNN, and MP-CNN, respectively. The results of 10 runs over the test set are shown in Table6. The bar plot is displayed in Figure 10, where “k-S,” “k-P,” and “”k-F1” stand for the sensitivity, precision, and F1 score for category k∈1,2,3,4.Table 6
Comparison of SPNN with four standard CNNs.
ModelClassSenPrcF1SC-CNNC193.5193.3593.43C286.7991.8789.26C394.4191.0192.68C494.5993.3793.97MA92.40AP-CNNC192.8693.0292.94C295.7692.9394.32C394.2693.9594.11C493.2293.2293.22MA93.18SPNN (ours)C198.0794.5996.30C291.7995.1993.45C394.9294.9294.92C495.2595.4095.32MA95.02L2P-CNNC188.7793.3691.01C293.2193.7293.46C394.4192.6893.53C493.6190.6392.10MA92.53MP-CNNC195.4494.7795.10C294.6491.7093.15C394.0793.4393.75C491.6495.7293.63MA93.91Figure 10
3D bar plot of SPNN against other network models.The results in Table6 and Figure 10 are coherent with our expectation that SPNN obtained the best results among all FM reduction approaches. The SPNN arrives at the sensitivities of all four categories are 98.07%, 91.79%, 94.92%, and 95.25%, respectively. The precisions of all four categories are 94.59%, 95.19%, 94.92%, and 95.40%, respectively. The F1-scores of the four categories are 96.30%, 93.45%, 94.92%, and 95.32%, respectively. The overall microaveraged F1 is 95.02%.In terms of microaveraged F1, the second-best algorithm is MP-CNN, which obtains a microaveraged F1 score of 93.91%. The third best is AP-CNN, with a microaveraged F1 score of 93.18%. The two comparably worst algorithms are SC-CNN and L2P-CNN, with the microaveraged F1 scores of 92.40% and 92.53%, respectively.SPNN obtains the best results because SP can prevent overfitting [17], which is the main shortcoming of max pooling. On the other hand, AP and L2P will average out the maximum activation values, which will impair the performances of convolutional neural network models. For SC-CNN, it only uses one quarter information of the input FM, and therefore may neglect those greatest values [34]. In all, this proposed SPNN can be regarded as an improved version of vanilla CNN models, where the SP is used to replace traditional MP.


### 4.2. PSSPNN versus SPNN

In this second experiment, we compared our two proposed network models, PSSPNN against SPNN, to validate the effectiveness of PatchShuffle. The results of 10 runs over the test set with different combinations of hyperparameters are shown in Table7, and the 3D bar chart is shown in Figure 11. The optimal hyperparameter we found from the 10-fold cross-validation of the nontest set is ε=0.05 and patch size is 2×2, which are coherent with reference [28].Table 7
Hyperparameter optimization of PSSPNN in terms of microaveraged F1 score.
εPatch size1×22×22×43×30.0195.1995.5495.1594.980.0595.4995.7995.4595.110.1095.2495.2895.5495.410.1595.2495.1195.2895.060.2095.0695.1594.8594.72Figure 11
3D bar chart of microaveraged F1 against patch size and probability.In addition, the PSSPNN with optimal hyperparameter is compared with SPNN. The results are shown in Table8. From the table, we can observe that PSSPNN provides better F1 values for all four categories and the overall microaverage, which shows the effectiveness of PatchShuffle. The reason is PatchShuffle adds regularization terms [28] in the loss function, and thus can improve the generalization ability of our SPNN model.Table 8
Comparison of SPNN against PSSPNN.
ModelClassSenPrcF1SPNN (ours)C198.0794.5996.30C291.7995.1993.45C394.9294.9294.92C495.2595.4095.32MA95.02PSSPN (ours)C197.8995.0696.46C292.8696.3094.55C395.7695.4495.60C496.5696.4096.48MA95.79


### 4.3. Comparison to State-of-the-Art Approaches

We compared our proposed PSSPNN method with 9 state-of-the-art methods: RCBO [7], 6L-CNN [8], RN-50 [9], RN-18 [10], CSSNet [11], COVNet [12], DeCovNet [13], 7L-CNN-CD [14], and FCONet [15]. All the comparison was carried on the same test set of 10 runs. The comparison results are shown in Table 9.Table 9
Comparison with state-of-the-art approaches.
ModelClassSenPrcF1RCBO [7]C171.9384.1977.58C272.8672.7372.79C373.5676.4174.96C480.6668.9174.32MA74.85RN-50 [9]C187.7285.0386.36C287.6891.2689.44C393.3989.8991.60C484.9287.6586.26MA88.41CSSNet [11]C194.0492.2593.14C293.7595.1194.42C391.3693.5892.45C494.4392.7593.58MA93.39DeCovNet [13]C191.0590.5890.81C293.7590.9992.35C390.5186.9788.70C488.6995.5892.01MA90.94FCONet [15]C192.2895.6493.93C296.7994.4395.59C394.7595.8895.31C494.9292.9493.92MA94.686L-CNN [8]C172.4683.9477.78C278.9377.8278.37C381.8675.0078.28C489.8487.5488.67MA80.94RN-18 [10]C182.8182.6682.73C281.0774.4377.61C374.2476.9875.58C482.1386.3884.20MA80.04COVNet [12]C189.8286.6388.20C289.8292.6391.21C393.7390.6692.17C487.3890.9689.13MA90.177L-CCD [14]C189.4793.5891.48C293.9392.4493.18C393.7395.1894.45C495.0891.3493.17MA93.09PSSPNN (ours)C197.8995.0696.46C292.8696.3094.55C395.7695.4495.60C496.5696.4096.48MA95.79For ease of comparison, Figure12 only compares the microaveraged F1 (MA F1) score of all algorithms, from which we can observe this proposed PSSPNN achieves the best performance among all the algorithms. This experiment is a simulation-based comparison. In the future, we will apply our algorithm to rigorous clinical testing and verification.Figure 12
Algorithm comparison in terms of MA F1.


### 4.4. Explainability of Proposed PSSPNN

We take Figure2 images as examples; the heatmaps of those four images are shown in Figures 13(a)–13(d), and the manual delineation, shown in Figures 13(e)–13(h), delineates the lesions of the three disease samples. Note there are no lesions of healthy control (HC) image. The NCSPM-5 feature map in PSSPNN was used to generate heatmaps with the Grad-CAM approach.Figure 13
Delineation of three diseased samples.
(a)
Heatmap of COVID-19(b)
Heatmap of CAP(c)
Heatmap of SPT(d)
Heatmap of HC(e)
Lesion of COVID-19(f)
Lesion of CAP(g)
Lesion of SPT(h)
No lesion for HCWe can observe from Figure13 that the heatmaps via our PSSPNN model and Grad-CAM can capture the lesions effectively and meanwhile neglect those nonlesion regions. Traditionally, AI is regarded as a “black box,” which impairs its widespread usage, e.g., the black box properties of traditional AI are problematic for the FDA. Nevertheless, with the help of explainability of modern AI techniques [35], the radiologist and patients will gain confidence in our proposed AI model, as the heatmap provides a clear and understandable interpretation of how AI predicts COVID-19 and other chest infectious disease from healthy subjects, which was also stated in reference [36]. Many new AI-based anatomic pathological systems now pass through FDA approval, such as whole slide images (WSI) [37], since the doctors know the relationships between the diagnosis and the explained answer.In the future, the explainability of our proposed AI model can be used in patient monitoring [38] and health big data [39]. Some novel network improvement and signal processing techniques may help our AI model in future researches, such as filters [40, 41], fuzzy [42, 43], edge computing [44], knowledge-aid [45, 46], autofocus [47], graph integration, and cross-domain knowledge exploitation [48–50].


## 4.1. Comparison of SPNN and Other Pooling Methods

In the first experiment, we compared the proposed SPNN with four standard CNNs with different pooling methods. The first CNN uses strided convolution in five modules to replace the stochastic pooling. The second to fourth comparison CNN models use L2P, AP, and MP, respectively. Those four baseline methods are called SC-CNN, L2P-CNN, AP-CNN, and MP-CNN, respectively. The results of 10 runs over the test set are shown in Table6. The bar plot is displayed in Figure 10, where “k-S,” “k-P,” and “”k-F1” stand for the sensitivity, precision, and F1 score for category k∈1,2,3,4.Table 6
Comparison of SPNN with four standard CNNs.
ModelClassSenPrcF1SC-CNNC193.5193.3593.43C286.7991.8789.26C394.4191.0192.68C494.5993.3793.97MA92.40AP-CNNC192.8693.0292.94C295.7692.9394.32C394.2693.9594.11C493.2293.2293.22MA93.18SPNN (ours)C198.0794.5996.30C291.7995.1993.45C394.9294.9294.92C495.2595.4095.32MA95.02L2P-CNNC188.7793.3691.01C293.2193.7293.46C394.4192.6893.53C493.6190.6392.10MA92.53MP-CNNC195.4494.7795.10C294.6491.7093.15C394.0793.4393.75C491.6495.7293.63MA93.91Figure 10
3D bar plot of SPNN against other network models.The results in Table6 and Figure 10 are coherent with our expectation that SPNN obtained the best results among all FM reduction approaches. The SPNN arrives at the sensitivities of all four categories are 98.07%, 91.79%, 94.92%, and 95.25%, respectively. The precisions of all four categories are 94.59%, 95.19%, 94.92%, and 95.40%, respectively. The F1-scores of the four categories are 96.30%, 93.45%, 94.92%, and 95.32%, respectively. The overall microaveraged F1 is 95.02%.In terms of microaveraged F1, the second-best algorithm is MP-CNN, which obtains a microaveraged F1 score of 93.91%. The third best is AP-CNN, with a microaveraged F1 score of 93.18%. The two comparably worst algorithms are SC-CNN and L2P-CNN, with the microaveraged F1 scores of 92.40% and 92.53%, respectively.SPNN obtains the best results because SP can prevent overfitting [17], which is the main shortcoming of max pooling. On the other hand, AP and L2P will average out the maximum activation values, which will impair the performances of convolutional neural network models. For SC-CNN, it only uses one quarter information of the input FM, and therefore may neglect those greatest values [34]. In all, this proposed SPNN can be regarded as an improved version of vanilla CNN models, where the SP is used to replace traditional MP.


## 4.2. PSSPNN versus SPNN

In this second experiment, we compared our two proposed network models, PSSPNN against SPNN, to validate the effectiveness of PatchShuffle. The results of 10 runs over the test set with different combinations of hyperparameters are shown in Table7, and the 3D bar chart is shown in Figure 11. The optimal hyperparameter we found from the 10-fold cross-validation of the nontest set is ε=0.05 and patch size is 2×2, which are coherent with reference [28].Table 7
Hyperparameter optimization of PSSPNN in terms of microaveraged F1 score.
εPatch size1×22×22×43×30.0195.1995.5495.1594.980.0595.4995.7995.4595.110.1095.2495.2895.5495.410.1595.2495.1195.2895.060.2095.0695.1594.8594.72Figure 11
3D bar chart of microaveraged F1 against patch size and probability.In addition, the PSSPNN with optimal hyperparameter is compared with SPNN. The results are shown in Table8. From the table, we can observe that PSSPNN provides better F1 values for all four categories and the overall microaverage, which shows the effectiveness of PatchShuffle. The reason is PatchShuffle adds regularization terms [28] in the loss function, and thus can improve the generalization ability of our SPNN model.Table 8
Comparison of SPNN against PSSPNN.
ModelClassSenPrcF1SPNN (ours)C198.0794.5996.30C291.7995.1993.45C394.9294.9294.92C495.2595.4095.32MA95.02PSSPN (ours)C197.8995.0696.46C292.8696.3094.55C395.7695.4495.60C496.5696.4096.48MA95.79


## 4.3. Comparison to State-of-the-Art Approaches

We compared our proposed PSSPNN method with 9 state-of-the-art methods: RCBO [7], 6L-CNN [8], RN-50 [9], RN-18 [10], CSSNet [11], COVNet [12], DeCovNet [13], 7L-CNN-CD [14], and FCONet [15]. All the comparison was carried on the same test set of 10 runs. The comparison results are shown in Table 9.Table 9
Comparison with state-of-the-art approaches.
ModelClassSenPrcF1RCBO [7]C171.9384.1977.58C272.8672.7372.79C373.5676.4174.96C480.6668.9174.32MA74.85RN-50 [9]C187.7285.0386.36C287.6891.2689.44C393.3989.8991.60C484.9287.6586.26MA88.41CSSNet [11]C194.0492.2593.14C293.7595.1194.42C391.3693.5892.45C494.4392.7593.58MA93.39DeCovNet [13]C191.0590.5890.81C293.7590.9992.35C390.5186.9788.70C488.6995.5892.01MA90.94FCONet [15]C192.2895.6493.93C296.7994.4395.59C394.7595.8895.31C494.9292.9493.92MA94.686L-CNN [8]C172.4683.9477.78C278.9377.8278.37C381.8675.0078.28C489.8487.5488.67MA80.94RN-18 [10]C182.8182.6682.73C281.0774.4377.61C374.2476.9875.58C482.1386.3884.20MA80.04COVNet [12]C189.8286.6388.20C289.8292.6391.21C393.7390.6692.17C487.3890.9689.13MA90.177L-CCD [14]C189.4793.5891.48C293.9392.4493.18C393.7395.1894.45C495.0891.3493.17MA93.09PSSPNN (ours)C197.8995.0696.46C292.8696.3094.55C395.7695.4495.60C496.5696.4096.48MA95.79For ease of comparison, Figure12 only compares the microaveraged F1 (MA F1) score of all algorithms, from which we can observe this proposed PSSPNN achieves the best performance among all the algorithms. This experiment is a simulation-based comparison. In the future, we will apply our algorithm to rigorous clinical testing and verification.Figure 12
Algorithm comparison in terms of MA F1.


## 4.4. Explainability of Proposed PSSPNN

We take Figure2 images as examples; the heatmaps of those four images are shown in Figures 13(a)–13(d), and the manual delineation, shown in Figures 13(e)–13(h), delineates the lesions of the three disease samples. Note there are no lesions of healthy control (HC) image. The NCSPM-5 feature map in PSSPNN was used to generate heatmaps with the Grad-CAM approach.Figure 13
Delineation of three diseased samples.
(a)
Heatmap of COVID-19(b)
Heatmap of CAP(c)
Heatmap of SPT(d)
Heatmap of HC(e)
Lesion of COVID-19(f)
Lesion of CAP(g)
Lesion of SPT(h)
No lesion for HCWe can observe from Figure13 that the heatmaps via our PSSPNN model and Grad-CAM can capture the lesions effectively and meanwhile neglect those nonlesion regions. Traditionally, AI is regarded as a “black box,” which impairs its widespread usage, e.g., the black box properties of traditional AI are problematic for the FDA. Nevertheless, with the help of explainability of modern AI techniques [35], the radiologist and patients will gain confidence in our proposed AI model, as the heatmap provides a clear and understandable interpretation of how AI predicts COVID-19 and other chest infectious disease from healthy subjects, which was also stated in reference [36]. Many new AI-based anatomic pathological systems now pass through FDA approval, such as whole slide images (WSI) [37], since the doctors know the relationships between the diagnosis and the explained answer.In the future, the explainability of our proposed AI model can be used in patient monitoring [38] and health big data [39]. Some novel network improvement and signal processing techniques may help our AI model in future researches, such as filters [40, 41], fuzzy [42, 43], edge computing [44], knowledge-aid [45, 46], autofocus [47], graph integration, and cross-domain knowledge exploitation [48–50].


## 5. Conclusion

In this paper, we proposed a PSSPNN, which entails five improvements: (i) proposed NCSPM module, (ii) usage of stochastic pooling, (iii) usage of PatchShuffle, (iv) improved multiple-way data augmentation, and (v) explainability via Grad-CAM. Those five improvements enable our AI model to deliver improved performances compared to 9 state-of-the-art approaches. The 10 runs on the test set showed our algorithm achieved a microaveraged F1 score of 95.79%.There are three shortcomings of our method, which will be resolved in the future: (i) the dataset currently contains three chest infectious diseases. In the future, we shall try to include more classes of chest diseases, such as thoracic cancer. (ii) Some new network techniques and models are not tested, such as transfer learning, wide network module design, attention mechanism, and graph neural network. Those advanced AI technologies will be studied. (iii) Our model does not go through strict clinical validation, so we will attempt to release our software to hospitals and get feedback from radiologists and consultants.



---

*Source: 6633755-2021-03-10.xml*