\section{Introduction}
\label{sec1}
The novel coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has emerged as one of the deadliest virus of the century, resulting in about 137 million people infected with over 2.9 million death worldwide as of April 2021. In the light of the unprecedented pandemic of COVID-19, public health systems have faced many challenges, including scarce medical resources, which are pushing healthcare providers to face the threat of infection
 \citep{ng2020covid}. Considering its ominously contagious nature, the early screening of COVID-19 infection becoming increasingly important to avert further spread of disease and thereby reduce the burden on saturated health care system.

Currently, the real-time polymerase chain reaction (RT-PCR) is considered as gold standard in diagnosis of COVID-19 for its high sensitivity and specificity \citep{tahamtan2020real}, but it takes several hours and even days depending on regions to get the  exam results due to overstressed laboratories. Since the majority of patients with confirmed COVID-19 present positive radiological findings, the radiologic examinations can be useful for rapid screening of disease \citep{shi2020radiological}. Although computed tomography (CT) scan has excellent sensitivity and specificity for COVID-19 diagnosis \citep{bernheim2020chest}, the  use of CT is a major burden because of its high cost and potential for cross-contamination in the radiology suite.
Therefore, CXR holds many practical advantages as primary screening tool in the pandemic situation.
In addition, CXR is useful for follow-up, which should be inexpensive and low in radiation exposure, to assess response to treatment.

Consequently, many studies have reported early application of CXR deep learning for  diagnosis \citep{wang2020covid, hemdan2020covidx, narin2020automatic, oh2020deep} or severity quantification of COVID-19 \citep{cohen2020predicting, sig2020covid, zhu2020deep, wong2020covidnet}, but they suffered from ineradicable drawbacks of poor generalization capability stemming from the scanty labelled COVID-19 data \citep{hu2020challenges, zech2018variable,roberts2021common}. The stable generalization performance on unseen data is indispensable for widespread adoption of the system \citep{roberts2021common}.

One of the most commonly used measures to solve this problem is to build a robust model with innumerable training data \citep{chen2020more}, but it is difficult to construct large-scale dataset with labeled COVID-19 cases under the current pandemic situation. As a result, several methods have been proposed to mitigate the problem by transfer learning \citep{apostolopoulos2020covid}, weakly supervised learning \citep{zheng2020deep,wang2020weakly}, and anomaly detection \citep{zhang2020covid}, but their performances are still suboptimal.

The previous studies mostly utilize convolutional neural network (CNN) models, which was not specially designed for manifestations of COVID-19 which can be characterized by bilateral involvement, peripheral and lower zone dominance of ground glass opacities and patchy consolidations \citep{cozzi2020chest}. Although CNN architecture has shown to be superb in many vision tasks, it may not be optimal for problems requiring high-level CXR disease classification, where the global characteristics like multiplicity, distribution, and patterns have to be considered. This is due to the
the intrinsic locality of pixel dependencies in the convolution operation.

 To overcome the similar limitation of CNN in computer vision problems that requires the integration of global relationshop between pixels, Vision Transformer (ViT) equipped with the Transformer architecture  \citep{vaswani2017attention} was proposed to model long-range dependency among pixels through the self-attention mechanism, showing the state-of-the-art (SOTA) performance in image classification task \citep{dosovitskiy2020image}.
Since the Transformer
   was originally invented for natural language processing (NLP) in order  to attend different positions of the input sequence within a corpus and compute a representation of that sequence,
the choice of an appropriate corpus is the prerequisite for  the Transformer design.

In the original paper  \citep{dosovitskiy2020image}, two ViT models were suggested utilizing either direct pixel-patch embedding or feature embedding by ResNet backbone as corpora for Transformer.
A problem occurs here, however, that neither the direct pixel-patch embedding nor feature embedding from ResNet may not be the optimal input embedding for CXR diagnosis of COVID-19. Fortunately, several large-scale CXR data sets are  constructed before COVID-19 pandemic and publicly available. For example, CheXpert \citep{irvin2019chexpert}, a large dataset that contain over 220,000 CXR images, provides labeled common low-level CXR findings (e.g. consolidation, opacity, edema, etc.), which is also useful for diagnosis of infectious disease. Moreover, an advanced CNN architecture has been suggested using the same dataset \citep{ye2020weakly}, which uses probabilistic class activation map (PCAM) pooling to leverage the class activation map to enhance the localization ability as well as classification performance. To take the maximum advantage of both the dataset and the network architecture for COVID-19, here we propose a novel ViT architecture which utilizes this advanced CNN architecture as a feature extractor for low-level CXR feature corpus, upon which Transformer is trained for downstream tasks of diagnosis  by utilizing the self-attention mechanism in Transformer.

It is worth to mention that our network basically identical to the text classification task with Transformer architecture, where the Transformer not only add up the meaning but also consider the location and relationship of words to make classification in sentence-level. Moreover, our method emulates the clinical experts who determine the final diagnosis of CXR (e.g. normal, bacterial pneumonia, COVID-19 infection, etc.) by comprehensively considering the low-level features with their pattern, multiplicity, location and distribution (e.g. \emph{Multiple} opacities and patch consolidations exist with \emph{lower lung zone dominance}: high probability for COVID-19) as illustrated in Fig. \ref{fig1}.

Another important contribution of this paper  is to show that our ViT framework can be also used for COVID-19 severity quantification and localization, enabling the serial follow-up of severity and thereby assisting the treatment decision of clinicians \citep{cohen2020covid}. The severity of COVID-19 can be determined by quantifying the extent of COVID-19 involvement.
Recently,  array-based simple severity annotations where 1 or 0 is assigned to each 6 subdivisions of lungs are proposed by \citet{toussie2020clinical}, and we are interested in utilizing this weak labeling approach for severity quantification.
As the Transformer output already incorporates the long-range relationship between regions through self-attention, we use this Transformer output  to design a light-weighted network that  can accurately
quantify and localize the COVID-19 extents from weak labels.
Specifically, we adopt  the region of interest (ROI) max-pooling of the output Transformer feature to bridge the severity map and simple array.
Consequently, in addition to the global severity score from 0 to 6, our model can create an intuitive severity level map where each pixel value explicitly means the likelihood of the presence of a COVID-19 lesion using
the weak array-based labels.
In summary,  our main contributions are as follows.
\begin{itemize}
\item A novel ViT model for COVID-19 is proposed by leveraging the low-level CXR feature corpus that contain the representations for common CXR findings with pre-built large-scale dataset.
\item We have not limited our model to diagnosis, but expanded our model to quantify severity to provide clinicians with clinical guidelines for making treatment decisions.
\item We experimentally demonstrated that our method outperforms other Transformer-based models as well as CNN-based models especially in terms of the generalization on unseen data.
\end{itemize}

The remainder of this paper is organized as follow. Section \ref{sec2} summarizes the related works. Section \ref{sec3} and Section \ref{sec4} describes the proposed framework and datasets, respectively. Experimental results are presented in Section \ref{sec5}. Finally, we conclude this work in Section \ref{sec6}.

\begin{figure*}[!t]
\centering
\includegraphics[width=0.95\textwidth]{Fig2.jpg}
\caption{Backbone network to extract low-level CXR feature corpus.}
\label{fig2}
\end{figure*}

\section{Related work}
\label{sec2}
\subsection{Vision Transformer}
Transformer \citep{vaswani2017attention}, which was originally invented for NLP, is a deep neural network based on self-attention mechanism that facilitates appreciably large receptive fields. After demonstrating its astounding performance,
not only has Transformer become a de facto standard practice in NLP, but it has also motivated the computer vision community to explore its applications in computer vision by taking advantage of the long-range dependency between pixels \citep{khan2021transformers}.

The ViT was the first major attempt to apply a pure Transformer directly to image, suggesting that it can completely replace the standard convolution operations by attaining SOTA performance. However, the experimental results showed that training vanilla ViT model requires a huge computational cost.
Therefore, the authors also suggested hybrid architecture by conjugating CNN backbone (e.g. ResNet) to Transformer. With the feature extracted by ResNet, the Transformer can mainly focus on modeling the global attention. The experimental results suggest that it was able to achieve higher performance with the hybrid approach with relatively small amount of computations.

After the introduction of ViT, the application of Transformer in computer vision has become an active area of investigation, resulting in many variant models of ViT showing SOTA performance in a variety of vision tasks including object detection \citep{zhu2020deformable}, classification \citep{dosovitskiy2020image, chen2020generative}, segmentation \citep{zheng2020rethinking}, and so on.

\subsection{Probabilistic Class Activation Map Pooling}
Class activation map (CAM) is a sort of class-specific saliency map obtained by quantifying the contribution of particular area of an image to the prediction of network. The most useful aspect of CAM is that it enables the localization of the important area only with weak labels, namely image-level supervision. Despite its excellent localization ability, most of previous works utilized CAM to generate heatmaps for lesion localization and visualization during inference. To leverage the localization ability of CAM to enhance the performance of network itself, one recent study utilized the CAM during training in CXR classification and localization tasks \citep{ye2020weakly}. They devised a novel global pooling operation that explicitly leverages the CAM in a probabilistic manner and is known as Probabilistic-CAM (PCAM) pooling. Different from standard approaches that use CAM for direct localization, they bound it with additional fully-connected layer and sigmoid function to get probabilities for each CXR findings. Then, the normalized attention weights were obtained from these output probability to make weighted feature maps containing more useful representation for each class. They showed that PCAM pooling operation can enhance both localization and diagnostic performance of  the model and achieved  first place in the 2019  CheXpert Challenge.

\subsection{COVID-19 Severity Quantification}
To build an automated algorithm for severity quantification, pixel-level annotation such as lesion segmentation label can offer a plentiful information. However, this type of labelling methods are labor-intensive
and collecting large data with this pixel-level annotated label is not feasible under the global pandemic of COVID-19. To alleviate the problem,  simplified severity annotation methods, such as score-based and array-based methods, have been proposed. For example, \citet{cohen2020predicting} suggested a geographic extent score and a lung opacity score based on a rating system of lung oedema proposed by \citet{warren2018severity}. A geographic extent score assigns scores that range from 0 to 4, while lung opacity score assigns values of 0 to 3 based on the severity of involvement in each lung area. \citet{borghesi2020covid} designed Brixia score, another array-type severity labeling method, dividing lung with anatomic landmarks and assign score of 0-3 to each subdivision. Similarly, \citet{toussie2020clinical} suggested an array-based severity score for COVID-19. After dividing both lungs into six divisions, each area is assigned a value of 0 or 1, depending on the presence of COVID-19 involvement, which adds up to an overall severity of 0 to 6. We adopted the array-based annotation method suggested by \citet{toussie2020clinical} for severity quantification of COVID-19.

\begin{figure*}[!h]
\centering
\includegraphics[width=0.9\textwidth]{Fig3.jpg}
\caption{Proposed framework and Vision Transformer model for diagnosis and severity quantification of COVID-19 on CXR for (A) disease classification, and (B) severity quantification.}
\label{fig3}
\end{figure*}

\section{The Proposed Framework}
\label{sec3}
One of the novel contributions  of our approach is to show that we can maximize the performance of the Transformer model by using the low-level CXR corpus that comes from the backbone network trained with a large well-curated public record to produce common CXR findings.
 As the backbone network is trained with large number of data, the network is less prone to overfitting and thereby the generalization capability can be improved.

\subsection{Pre-training Backbone Network for Low-level Feature Corpus}
As a backbone network to extract low-level feature, we used the modified version of network proposed by \citet{ye2020weakly}. Firstly, the backbone network was pre-trained to classify 10 common low-level findings with a large public dataset. As depicted in Fig. \ref{fig2}, feature maps in each layer can be the candidates for utilizable feature embedding for the subsequent Transformer, and we experimentally found that the common embedding before the PCAM operation comprises of  most useful information. Nevertheless, care should be exercised since the PCAM operation for specific low-level CXR findings (e.g. lung opacity, consolidation, etc.) turns out to be crucial to achieve the optimal embedding at intermediate level, as  PCAM aligns these features to obtain the better classification results. More detailed experimental results about the level of feature map will be provided using ablation studies in Section \ref{sec5}.6.2.

\subsection{Vision Transformer for Classification}
The overall framework and the architecture of our ViT model is provided in Fig. \ref{fig3}.
Specifically,  for a given $H\times W$ size input image $\boldsymbol{x} \in \mathbb{R}^{{H}\times{W}}$,
 the backbone network $\boldsymbol{\mathcal{G}}$  generates  $H'\times W'$ size feature maps $\boldsymbol{F}$
\begin{flalign}
&& \boldsymbol{F}&= \boldsymbol{\mathcal{G}}(\boldsymbol{x}) &
\end{flalign}
Here,  the feature tensor  $\boldsymbol{F} \in \mathbb{R}^{{H'}\times{W'}\times{C'}}$ is defined as
\begin{flalign}
&&\boldsymbol{F} &= \begin{bmatrix}\boldsymbol{f}_{1} &\boldsymbol{f}_{2} &\cdots & \boldsymbol{f}_{{H'}\times{W'}}\end{bmatrix} &
\end{flalign}
where
 $\boldsymbol{f}_n \in \mathbb{R}^{{C'}}$   denotes a $C$-dimensional  embedded representation of low-level features at the $n$-th encoded block.
 These feature vectors  are used
 to construct the low-level CXR feature corpora for Transformer.

Then, similar to Bidirectional Encoder Representations from Transformers  (BERT) \citep{devlin2018bert}, our Vision Transformer uses Transformer encoder layers to the input embedding. Specifically, since the Transformer encoder utilizes constant latent vector of dimension $D$, the extracted $C'$ dimension feature $\boldsymbol{f}_n \in \mathbb{R}^{{C'}}$ is first projected to a $D$ dimension feature $\tilde \boldsymbol{f}_n \in \mathbb{R}^{{D}}$ using $1$ $\times$ $1$ convolution kernel.
 We then prepended learnable \texttt{[class]} token embedding vector $\boldsymbol{f}_\texttt{cls} \in \mathbb{R}^{{D}}$ to projected feature tensor.
This leads to the following composite projected feature tensor:
\begin{align}
\tilde\boldsymbol{F} &= \begin{bmatrix}\boldsymbol{f}_\texttt{cls} & \tilde\boldsymbol{f}_{1} &\tilde\boldsymbol{f}_{2} &\cdots & \tilde\boldsymbol{f}_{{H'}\times{W'}}\end{bmatrix
\end{align}
 A positional embedding $\textbf{E}_{pos}$ that has the same shape to the projected feature tensor $\tilde\boldsymbol{F}$ is then added to encode a notion of the sequential order:
 \begin{align*}
 \boldsymbol{Z}^{(0)}
 &=\tilde\boldsymbol{F} + \textbf{E}_{pos}
 \end{align*}
This is then used as an input to a Transformer composed of $L$ successive encoder layers:
\begin{align}
\boldsymbol{Z}^{(l)} = \boldsymbol{\mathcal{T}}^{(l)}\left(\boldsymbol{Z}^{(l-1)} \right),\quad l=1,\cdots,L
\end{align}
where $\boldsymbol{Z}^{(l)}=\begin{bmatrix}\boldsymbol{z}_0^{(l)}&\boldsymbol{z}_1^{(l)}& \cdots & \boldsymbol{z}_{H'\times W'}^{(l)} \end{bmatrix}$
and $\boldsymbol{\mathcal{T}}^{(l)}$ denotes the $l$-th encoder layer.
The encoder layers used in our model are the same as standard Transformer which consists of repeated layers of multihead self-attention (MSA), multilayer perceptron (MLP), layer normalization (LN), and residual connections in each block, as shown in Fig.~\ref{fig3}.
Then, the first column $\boldsymbol{z}_0^{(L)}$ of $\boldsymbol{Z}^{(L)}$ represents
the Transformer attended feature vector with respect to
 the \texttt{[class]} token.
 Therefore, by simply adding a linear classifier as the classification head,
 we can obtain the diagnosis result of the input CXR image $\boldsymbol{x}$ (see Fig.~\ref{fig3}(A)).

\begin{figure*}[!t]
\centering
\includegraphics[width=0.85\textwidth]{Fig4.jpg}
\caption{Procedure of severity prediction and labeling. (A) Map head and ROI max-pooling of proposed framework. (B) Our severity annotation method for severity quantification on CXRs.}
\label{fig4}
\end{figure*}

For the interpretability of the classification model,
we adopted a visualization method of saliency map tailored for ViT suggested by \citep{chefer2020transformer}, which computes relevancy for Transformer network.  Specifically,
unlike the traditional approaches of gradient propagation methods \citep{selvaraju2017grad, smilkov2017smoothgrad, srinivas2019full} or attribution propagation methods \citep{bach2015pixel, gu2018understanding}, which rely on the heuristic propagation along attention graph or the obtained attention maps, the method  in \citet{chefer2020transformer} calculate the local relevance with deep Taylor decomposition, which is then propagated throughout the layers. This relevance propagation method is especially useful for models based on Transformer architecture, as it overcomes the problem of self-attention operations and skip connections.

\subsection{Vision Transformer for Severity Quantification}
In the classification task, only one Transformer output at the \texttt{[class]} token position is used.
However, the rest of the Transformer output also produces feature embedding at each block position by taking into account of long-range relations
between the blocks. Therefore, we conjecture that  this information is useful
 for the severity quantification, as severity is determined by both local and global manifestation of the disease.
Accordingly, as shown in Fig. \ref{fig3}, these outputs are combined by additional lightweighted network to produce the COVID-19 severity map.

Specifically, as shown in Figs. \ref{fig3}(b) and  \ref{fig4},
 we first extract the Transformer output $\boldsymbol{Z}^{(L)}$ except the \texttt{[class]} token position:
\begin{align}
\boldsymbol{Z}_{res}   &=  \begin{bmatrix}\boldsymbol{z}_1^{(L)} & \cdots & \boldsymbol{z}_{H'\times W'}^{(L)}\end{bmatrix}
\end{align}
which is used as an input to the map head network $\boldsymbol{\mathcal{N}}$
\begin{align}
\boldsymbol{S}=\boldsymbol{\mathcal{N}}(\boldsymbol{Z}_{res})
\end{align}
Then, the network output is multiplied pixel-wise with the segmentation mask $\boldsymbol{M}$, after which
ROI max-pooling (RMP) is applied to generate the severity mask $\boldsymbol{Y}_{sev}\in \mathbb{R}^{3\times 2}$:
\begin{align}
\boldsymbol{Y}_{sev} =  \texttt{RMP}\left(\boldsymbol{S}\otimes \boldsymbol{M}\right)
\end{align}
where $\otimes$ denotes the Hadamard product. In detail, the lung was divided into a total of six subdivisions, by dividing the right and left lungs into three subdivision (upper, middle, lower zone) with 5/12 and 2/3 line. Next, the largest values within each six subdivision were assigned as predicted values of severity array. Then, the map head network is trained by minimizing error of the estimated severity array with respect to the weakly annotated severity label as in Fig.~\ref{fig4}.

To generate the lung segmentation mask, we used method introduced by \citet{oh2021unifying}.
In contrast to the existing approaches that are prone to under-segmentation for the severely infected lung with large consolidations,
this novel approach enables the accurate segmentation of abnormal lung as well as normal lung area by learning common features
using
 a single generator with AdaIN layers.
 Since a single generator is used for all these tasks  by simply changing the AdaIN codes, the generator can synergistically learn the common features
 to improve segmentation performance for abnormal CXR data.

\begin{table*}[!t]
  \caption{Summary of dataset resources and disease classes for PA view CXRs.}
  \label{table:dataset_PA}
  \centering
  \resizebox{0.90\textwidth}{!}
  {%
  \begin{tabular}{M{2cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}}
  \Xhline{1.2pt}
  \multirow{2}{*} {\textbf{PA view}} & \multirow{2}{*} {\textbf{Total}} & \multicolumn{3}{c|}{\textbf{External test}} & \multicolumn{4}{c}{\textbf{Training and Validation}}\\
  \cline{3-9}
  & & {CNUH} & {YNU}& {KNUH}  & {AMC} & {NIH} & {Brixia} & {BIMCV} \\
  \cline{2-3} \cline{6-7}
  \hline
   \multicolumn{1}{l|}{Normal} & 13,649 & 320 & 300 & 400 & 8,861 & 3,768 & {-} & {-}\\
  \multicolumn{1}{l|}{Other infection} & 1,468 & 39 & 144 & 308 & 977 & {-} & {-} & {-} \\
   \multicolumn{1}{l|}{{COVID-19}}  & 2,431 & 6 & 8 & 80 & {-} & {-} & 1,929 & 408\\
   \cline{1-9}
   \multicolumn{1}{l|}{{Total images}}  & 17,548 & 365 & 452 & 788 & 9,838 & 3,768 & 1,929 & 408\\
   \Xhline{1.2pt}
  \end{tabular}}
\end{table*}

\begin{table*}[!t]
  \caption{Summary of dataset resources and disease classes for AP view CXRs.}
  \label{table:dataset_AP}
  \centering
  \resizebox{0.98\textwidth}{!}
  {%
  \begin{tabular}{M{2.0cm}|M{1.2cm}|M{2.0cm}|M{1.2cm}|M{1.2cm}|M{1.2cm}|M{1.2cm}|M{1.2cm}|M{1.2cm}|M{1.2cm}|M{1.2cm}}
  \Xhline{1.2pt}
  \multirow{2}{*} {\textbf{AP view}} & \multirow{2}{*} {\textbf{Total}} & \textbf{External test} & \multicolumn{8}{c}{\textbf{Training and Validation}}\\
  \cline{3-11}
  & & {CNUH} & {YNU}& {KNUH}  & {AMC} & {NIH} & {CheXpert} & {PADchest} & {Brixia} & {BIMCV} \\
  \cline{2-3} \cline{6-11}
  \hline
   \multicolumn{1}{l|}{Normal} & 14,507 & 97 & {-} & {-} & 117 & 3,390 & 9,500 & 1,310 & {-} & 93\\
  \multicolumn{1}{l|}{Other infection} & 204 & 19 & 76 & 92 & 17 & {-} & {-} & {-} & {-} & {-} \\
   \multicolumn{1}{l|}{{COVID-19}}  & 3,334 & 75 & 278 & 213 & {-} & {-} & {-} & {-} & 2,384 & 374\\
   \cline{1-11}
   \multicolumn{1}{l|}{{Total images}}  & 18,045 & 191 & 354 & 305 & 134 & 3,390 & 9,500 & 1,310 & 2,384 & 467\\
   \Xhline{1.2pt}
  \end{tabular}}
\end{table*}

\section{Dataset}
\label{sec4}
Datasets used for this study can be divided into three: dataset for pre-training backbone, dataset for classification, dataset for severity quantification.

\subsection{Dataset for Pre-training}
For the pre-training of the backbone netowrk to extract the low-level CXR features, we used CheXpert dataset containing 10 labeled CXR findings: no finding, cardiomegaly, opacity, edema, consolidation, pneumonia, atelectasis, pneumothorax, pleural effusion, and support device. With a total of 224,316 CXR images from 65,240 subjects, the 32,387 lateral view images were excluded, leaving 29,420 PA and 161,427 AP view data available. With this large number of CXRs, it was able to train the backbone network robust to the variation in subjects, which is one of the key strengths of our model.

\subsection{Dataset for Classification}
Table \ref{table:dataset_PA} and Table \ref{table:dataset_AP} summarizes dataset resources and partitioning for PA and AP views, respectively. To train and evaluate the Transformer model, we utilized both public datasets containing labeled cases of infectious disease (Valencian Region Medical Image Bank [BIMCV] \citep{de2020bimcv}, Brixia \citep{signoroni2020end}, National Institutes of Health [NIH] \citep{wang2017chestx}) and deliberately collected CXR data from four hospitals (Asan Medical Center [AMC], Seoul, Korea; Chonnam National Univerity Hospital [CNUH], Daejeon, Korea; Yeungnam University Hospital [YNU], Daegu, Korea; Kyungpook National University Hospital [KNUH], Daegu, Korea) labeled by board-certified radiologists for this study. Finally, the integrated dataset was divided into three label classes including normal, other infections (e.g. bacterial infection, tuberculosis) and COVID-19 infection, considering the application in real clinical setting. For PA view images, we put 3 institutional data (CNUH, YNU, KNUH) aside as external test datasets to evaluate the generalization capability by using data collected from independent hospitals with different devices and settings. On the other hand, CNUH data was solely used as the external test dataset for AP view images, since it was the only dataset containing all three label classes.

\subsection{Dataset for Severity Quantification}
Table \ref{table:dataset_sev} summarizes dataset resources and global severity levels. Different from diagnosis, the PA and AP view data were integrated and utilized without division for severity quantification task, since there is possibility that follow-up images may be obtained with both PA and AP view even in a single patient. Two board certified radiologists labeled the severity for three institutional datasets (CNUH, YNU, KNUH) using the array-based severity labeling method of \citet{toussie2020clinical} as in Fig. \ref{fig6}. We also utilized publicly available data, Brixia dataset, after translating its severity score the same as that of the institutional datasets. We alternately used one institutional dataset as an external testset and trained the models with two remaining datasets together with the Brixia dataset to evaluate the generalization capability in various external settings. Besides, 12 COVID-19 cases from BIMCV dataset were used to compare the severity map generated by our model to those annotated by clinical experts.

\begin{table}[!h]
  \caption{Summary of dataset resource and global severity level for severity labeled CXRs.}
  \label{table:dataset_sev}
  \centering
  \resizebox{0.4\textwidth}{!}
  {\small
  \begin{tabular}{c|c|c|c|c|c
  \Xhline{1.2pt}
  {\textbf{Severity}} & {\textbf{Total}} & {\textbf{CNUH}} & {\textbf{YNU}} & {\textbf{KNUH}} & {\textbf{Brixia}} \\
  \cline{1-6}
  \hline
   \multicolumn{1}{c|}{1}  & 361   & 26   & 63  & 25   & 247    \\
   \multicolumn{1}{c|}{2}  & 521   & 11   & 59  & 22   & 429    \\
   \multicolumn{1}{c|}{3}  & 448   & 8    & 25  & 18   & 397    \\
   \multicolumn{1}{c|}{4}  & 920   & 7    & 35  & 31   & 847    \\
   \multicolumn{1}{c|}{5}  & 774   & 12   & 18  & 29   & 715    \\
   \multicolumn{1}{c|}{6}  & 1,758 & 17   & 86  & 171  & 1,484  \\
   \cline{1-6}
   \multicolumn{1}{c|}{{Total}}  & 4,782 & 81   & 286 & 296  & 4,119 \\
   \Xhline{1.2pt}
  \end{tabular}}
\end{table}

\begin{table*}[!t]
  \caption{Diagnostic performance of the proposed model in various external test datasets from three different institutions for PA view.}
  \label{table:ext_PA}
  \centering
  \resizebox{0.9\textwidth}{!}
  {%
  \begin{tabular}{M{1.06cm}|M{1.02cm}|M{1.0cm}|M{1.04cm}|M{1.0cm}|M{1.02cm}|M{1.0cm}|M{1.04cm}|M{1.0cm}|M{1.02cm}|M{1.0cm}|M{1.04cm}|M{1.0cm}}
  \Xhline{1.2pt}
  \multicolumn{1}{l|}{\textbf{Metrics}} & \multicolumn{4}{c|}{\textbf{External dataset 1 (CNUH)}} & \multicolumn{4}{c|}{\textbf{External dataset 2 (YNU)}} & \multicolumn{4}{c}{\textbf{External dataset 3 (KNUH)}} \\
  \cline{2-13}
  & {Avg.} & {Normal} & {Others} & {COVID} & {Avg.} & {Normal} & {Others} & {COVID} & {Avg.} & {Normal} & {Others} & {COVID} \\
  \hline
   \multicolumn{1}{l|}{AUC} & 0.932 & 0.938 & 0.926 & 0.931 & 0.921 & 0.947 & 0.908 & 0.907 & 0.928 & 0.955 & 0.908 & 0.921 \\
  \multicolumn{1}{l|}{Sensitivity} & 83.4 & 84.7 & 82.1 & 83.3 & 88.4 & 92.3 & 85.4 & 87.5 & 85.4 & 87.3 & 87.7 & 81.3 \\
   \multicolumn{1}{l|}{Specificity} & 85.3 & 88.9 & 88.3 & 78.6 & 84.2 & 90.8 & 82.5 & 79.3 & 86.8 & 89.4 & 82.1 & 88.8 \\
   \multicolumn{1}{l|}{Accuracy} & 83.8 & 85.2 & 87.7 & 78.6 & 84.9 & 91.8 & 83.4 & 79.4 & 86.9 & 88.3 & 84.3 & 88.1 \\
   \Xhline{1.2pt}
  \end{tabular}}
\end{table*}

\subsection{Details of Implementation and Evaluation}
The CXR images were preprocessed via histogram equalization, Gaussian blurring with $3 \times 3$ kernel, normalization, and finally resized to $512 \times 512$. As our backbone network, the modified version of the network proposed by \citet{ye2020weakly}, which comprises DenseNet-121 baseline followed by PCAM operations. Among several layers of intermediate feature maps, we used the feature map of size $16 \times 16 \times 1024$ just before the PCAM operation. For subsequent Transformer architecture, we used standard Transformer model with 12 layers and 12 heads per each layer, while comparison with Transformer architecture with different network size is also provided in Section \ref{sec5}.6.1.

For pre-training of the backbone network, Adam optimizer with learning rate of 0.0001 was used. We trained the backbone network for 160,000 optimization steps with step decay scheduler with batch size of 8. For training of classification model, SGD optimizer with momentum 0.9 was used with learning rate of 0.001. A max gradient norm of 1 was applied to stabilize training. We trained the model for 10,000 optimization steps with cosine warm-up scheduler (warm-up steps = 500) with batch size of 16. We trained two individual classification models for PA and AP view images, respectively. For severity quantification, a map head with four upsizing convolutional blocks is used, with last block followed by sigmoid non-linearity which squashes output into [0-1] range. Training of severity quantification model was done with SGD optimizer with learning rate of 0.003 for 12,000 optimization steps with constant learning rate, and batch size of 4 was used. These optimal hyperparameters were determined experimentally.

We used the area under receiver operating characteristic curve (AUC) as the evaluation metrics for diagnostic performance of the classification model, but also calculated sensitivity, specificity, and accuracy after adjusting the thresholds to meet the sensitivity value of $\geq 80 \%$, if possible. As evaluation metrics for severity quantification, we used the Mean Squared Error (MSE) as the main metric, but the Mean Absolute Error (MAE), Correlation Coefficient (CC), and $R^2$ score were also measured and compared. For classification model, the comparison between models and ablation studies were performed with PA view data, since it usually offers better diagnostic quality and is standard position for CXR diagnosis of lung disease, while both PA and AP view CXRs were used in severity quantification model.

All experiments including preprocessing, development and evaluation of the model, was performed using Python version 3.7 and Pytorch library version 1.7 on Nvidia Tesla V100 and NVidia RTX 3090.

\section{Experimental results}
\label{sec5}

\subsection{Diagnostic Performance on External Test Datasets}
The diagnostic performances of the proposed model for PA view images are provided in Table \ref{table:ext_PA}. On average of 3 label classes (normal, other infection, COVID-19), our model showed the mean AUCs of 0.932, 0.947, 0.928, sensitivities of 83.4\%, 88.4\%, 85.4\%, specificities of 85.3\%. 84.2\%, 86.8\%, and accuracy of 83.8\%, 84.9\%, 86.9\% in three external institutions, which confirmed the stable performance (AUC $\geq 0.900$) and the generalization capability of our method in clinical situations with different devices and settings. For AP view images, our model showed mean AUCs of 0.890, 0.880, 0.828 for each 3 label classes (Table \ref{table:ext_AP}), which was slightly decreased compared to those of PA view images but still showed fair performance (AUC $\geq 0.800$) in the external test dataset, considering the fact that the diagnosis of infectious disease using only AP view image is not standard and usually deteriorates the diagnostic performance.

\begin{figure}[!t]
\centering
\includegraphics[width=0.43\textwidth]{Fig5.jpg}
\caption{Examples of visualization results for each disease classes. (A) Bacterial infection, (B) tuberculosis, and (C) COVID-19 infection.}
\label{fig5}
\end{figure}

\begin{table}[!h]
  \caption{Diagnostic performance in external test dataset for AP view.}
  \label{table:ext_AP}
  \centering
  \resizebox{0.45\textwidth}{!}
  {%
  \begin{tabular}{M{1.4cm}|M{1.4cm}|M{1.4cm}|M{1.4cm}|M{1.4cm}}
  \Xhline{1.2pt}
  \multicolumn{1}{l|}{\textbf{Metrics}} & \multicolumn{4}{c}{\textbf{External dataset (CNUH)}} \\
  \cline{2-5}
  & \multicolumn{1}{c|}{Avg.} & \multicolumn{1}{c|}{Normal} & \multicolumn{1}{c|}{Others} & \multicolumn{1}{c}{COVID} \\
  \hline
   \multicolumn{1}{l|}{AUC} & 0.866 & 0.890 & 0.880 & 0.828 \\
   \Xhline{1.2pt}
  \end{tabular}}
\end{table}

\subsection{Model Interpretability Results}
Fig. \ref{fig5} exemplifies the visualization of saliency maps for each disease classes in the external test datasets. As shown in the examples, our model well-localized a focal infected area either by bacterial infection (Fig. \ref{fig5} (a)) or tuberculosis (Fig. \ref{fig5} (b)), while it was also able to delineate the multi-focal lesions in periphery of both lower lungs in Fig. \ref{fig4} (c), which is typical findings for COVID-19 pneumonia.

\begin{figure*}[!t]
\centering
\includegraphics[width=0.90\textwidth]{Fig6.jpg}
\caption{Examples of severity quantification results of our models on CNUH external dataset.}
\label{fig6}
\end{figure*}

\subsection{Severity Quantification Results on External Test Datasets}
The results of severity quantification of our model are shown in Table \ref{table:ext_sev}. Our model showed   the MSE of 1.682, 1.677, 1.607, the MAE of 1.028, 1.102, 0.930, correlation coefficient of 0.781, 0.777, 0.682, and $R^2$ score of 0.572, 0.572, 0.432 in three external institutions. Brixia dataset contains a consensus subset of 150 CXR images labeled by five independent radiologists. Within this subset, the average MSE between the consensus severity score calculated from majority voting and each radiologist's rating is 1.683. As a result, the MSEs of 1.657, 1.696, and 1.676 in three external institutions show our model's performance comparable to those of experienced radiologists and generalization capability in the clinical environment.

Fig \ref{fig6} illustrates the examples of severity quantification, including the predicted scores, arrays, maps, and lesion contours in one of the external test datasets, which confirms that not only can our model correctly predict global severity, but it also generates an intuitive severity map that highlights the affected area, which can also be used to contour lesions.

Finally, Fig \ref{fig7} exemplifies the comparison between ground truth segmenteation label of involved area and model's prediction of involvement in BIMCV dataset. As shown in the figure, the model generally well-localized the areas of involvement.

\begin{figure}[!t]
\centering
\includegraphics[width=0.45\textwidth]{Fig7.jpg}
\caption{Comparison of localization results in BIMCV dataset. Green: radiologist's annotation. Yellow: model's prediction after thresholding.}
\label{fig7}
\end{figure}

\begin{table}[!t]
  \caption{Severity quantification performance of the proposed model in various external test datasets from different institutions.}
  \label{table:ext_sev}
  \centering
  \small
  \resizebox{0.48\textwidth}{!}
  {%
  \begin{tabular}{l|c|c|c}
  \Xhline{1.2pt}
\multicolumn{1}{l|}{\textbf{Metrics}}
& \thead{\textbf{External dataset 1}\\\textbf{(CNUH)}}
& \thead{\textbf{External dataset 2}\\\textbf{(YNU)}}
& \thead{\textbf{External dataset 3}\\\textbf{(KNUH)}} \\ \hline
\multicolumn{1}{l|}{MSE}   & 1.682   & 1.677   & 1.607 \\
\multicolumn{1}{l|}{MAE}   & 1.028   & 1.102   & 0.930 \\
\multicolumn{1}{l|}{CC}    & 0.781   & 0.777   & 0.682 \\
\multicolumn{1}{l|}{$R^2$}   & 0.572   & 0.572   & 0.432 \\
   \Xhline{1.2pt}
   \end{tabular}}
\end{table}

\begin{table*}[!t]
  \caption{Comparison with CNN and other Transformer-based Models.}
  \label{table:comparison}
  \centering
  \resizebox{\textwidth}{!}
  {%
  \begin{tabular}{M{1.06cm}|M{1.02cm}|M{1cm}|M{1.04cm}|M{1cm}|M{1.02cm}|M{1cm}|M{1.04cm}|M{1cm}|M{1.02cm}|M{1cm}|M{1.04cm}|M{1cm}}
  \Xhline{1.2pt}
   \multicolumn{1}{l|}{\textbf{Methods}} & \multicolumn{4}{c|}{\textbf{External dataset 1 (CNUH)}} & \multicolumn{4}{c|}{\textbf{External dataset 2 (YNU)}} & \multicolumn{4}{c}{\textbf{External dataset 3 (KNUH)}}
  \\
  \cline{2-13}
  & Avg. & {Normal} & {Others} & {COVID} & Avg. & {Normal} & {Others} & {COVID} &Avg. & {Normal} & {Others} & {COVID} \\
  \hline
  \multicolumn{1}{l}{\textbf{CNN-based}} \\
   \multicolumn{1}{l|}{ResNet-50} & 0.843 & 0.725 & 0.883 & 0.922 & 0.816 & 0.768 & 0.914 & 0.766 & 0.846 & 0.755 & 0.874 & 0.910 \\
  \multicolumn{1}{l|}{ResNet-152} & 0.823 & 0.755 & 0.893 & 0.821 & 0.841 & 0.811 & 0.924 & 0.789 & 0.828 & 0.795 & 0.800 & 0.890 \\
   \multicolumn{1}{l|}{DenseNet-121} & 0.852 & 0.740 & 0.907 & 0.908 & 0.880 & 0.820 & \textbf{0.937} & 0.883 & 0.835 & 0.795 & 0.792 & 0.920 \\
   \hline
   \multicolumn{1}{l}{\textbf{Transformer-based}} \\
   \multicolumn{1}{l|}{ViT} & 0.885 & 0.883 & 0.879 & 0.894 & 0.900 & 0.901 & 0.906 & 0.895 & 0.838 & 0.855 & 0.753 & 0.906 \\
   \multicolumn{1}{l|}{ViT(hybrid)} & 0.857 & 0.865 & 0.848 & 0.858 & 0.902 & 0.936 & 0.885 & 0.883 & 0.857 & 0.889 & 0.760 & \textbf{0.922} \\
   \multicolumn{1}{l|}{\textbf{Ours}} & \textbf{0.932} & \textbf{0.938} & \textbf{0.926} & \textbf{0.931} & \textbf{0.921} & \textbf{0.947} & 0.908 & \textbf{0.907} & \textbf{0.928} & \textbf{0.955} & \textbf{0.908} & 0.921 \\
   \cline{1-13}
   \Xhline{1.2pt}
  \end{tabular}}
\end{table*}

\begin{table}[!t]
  \caption{Comparison of severity quantification performance with CNN and other Transformer based Models.}
  \label{table:comparison_sev}
  \centering
  \resizebox{0.48\textwidth}{!}
  {\large
  \begin{tabular}{l|c|c|c}
  \Xhline{1.2pt}
\multicolumn{1}{l|}{\textbf{Methods}}
& \thead{\textbf{External dataset 1}\\\textbf{(CNUH)}}
& \thead{\textbf{External dataset 2}\\\textbf{(YNU)}}
& \thead{\textbf{External dataset 3}\\\textbf{(KNUH)}} \\ \hline
\multicolumn{1}{l}{\textbf{CNN-based}}  & \multicolumn{1}{l}{}& \multicolumn{1}{l}{}& \multicolumn{1}{l}{} \\
\multicolumn{1}{l|}{ResNet-50}        & \textbf{1.649}   & 2.072   & 2.024  \\
\multicolumn{1}{l|}{ResNet-152}       & 1.718            & 2.135   & 2.026   \\
\multicolumn{1}{l|}{DenseNet-121}     & 1.839            & 1.963   & 1.693   \\ \hline
\multicolumn{1}{l}{\textbf{Transformer-based}}           & \multicolumn{1}{l}{}& \multicolumn{1}{l}{}& \multicolumn{1}{l}{}\\
\multicolumn{1}{l|}{ViT}              & 3.416            & 3.918   & 3.348   \\
\multicolumn{1}{l|}{ViT(hybrid)}      & 1.894            & 2.781   & 2.099   \\
\multicolumn{1}{l|}{\textbf{Ours}}    & 1.682            & \textbf{1.677}  & \textbf{1.607} \\
   \Xhline{1.2pt}
   \multicolumn{4}{l}{\footnotesize{Mean squared error (MSE) of the severity score (0-6) is chosen as the metric for comparison.} \par}
  \end{tabular}}
\end{table}

\subsection{Comparison with CNN and Transformer-based Models}
To compare the performance with other CNN-based and Transformer-based SOTA models, ResNet-50, ResNet-152, DenseNet-121, the standard ViT and hybrid ViT models were tested with the same external test datasets. All models except ours were initialized with ImageNet pre-traiend weights as it significantly improved the performance. Other training and evaluation settings were kept the same for fair comparison. As suggested in Table \ref{table:comparison} and Table \ref{table:comparison_sev}, our model outperformed not only the CNN-base models but also the Transformer-based models in all of the external test datasets for both classification and severity quantification tasks, demonstrating superb and stable performance in real-world application. The performance improvement was not just the result of increased complexity of the model, considering that our model also considerably outperformed the model that has increased complexity (ResNet-152).

\subsection{Ablation Studies}
To get better understanding of our models, we conducted a series of ablation studies, as provided in Table \ref{table:ablation} and Table \ref{table:ablation_sev}.  More details are as follows.

\begin{table*}[!t]
  \caption{Ablation studies for classification.}
  \label{table:ablation}
  \centering
  \resizebox{1.0\textwidth}{!}
  {\large
  \begin{tabular}{l|c|c|c|c}
  \Xhline{1.2pt}
   \multicolumn{1}{l|}{\textbf{Methods}} & \multicolumn{1}{c|}{\textbf{Avg. of 3 External datasets}} & \multicolumn{1}{c|}{\textbf{External dataset 1 (CNUH)}} & \multicolumn{1}{c|}{\textbf{External dataset 2 (YNU)}} & \multicolumn{1}{c}{\textbf{External dataset 3 (KNUH)}}
  \\
  \cline{2-5}
  & \multicolumn{1}{c|}{\textbf{Avg. AUC of 3 classes}} & \multicolumn{1}{c|}{\textbf{Avg. AUC of 3 classes}} & \multicolumn{1}{c|}{\textbf{Avg. AUC of 3 classes}} & \multicolumn{1}{c}{\textbf{Avg. AUC of 3 classes}}\\
  \hline
  \multicolumn{1}{l}{\textbf{Network size}} \\
   \multicolumn{1}{l|}{2 layers and 4 heads} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.914} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.929} \\
  \multicolumn{1}{l|}{4 layers and 8 heads} & \multicolumn{1}{c}{0.915} & \multicolumn{1}{c}{0.914} & \multicolumn{1}{c}{0.916} & \multicolumn{1}{c}{0.915} \\
   \multicolumn{1}{l|}{12 layers and 12 heads \textbf{(ours)}} & \multicolumn{1}{c}{\textbf{0.927}} & \multicolumn{1}{c}{0.932} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.928} \\
   \hline
  \multicolumn{1}{l}{\textbf{Intermediate feature map}} \\
   \multicolumn{1}{l|}{feature before PCAM \textbf{(ours)}} & \multicolumn{1}{c}{\textbf{0.927}} & \multicolumn{1}{c}{0.932} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.928} \\
   \multicolumn{1}{l|}{1 feature after PCAM} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.923} & \multicolumn{1}{c}{0.908} & \multicolumn{1}{c}{0.933} \\
   \multicolumn{1}{l|}{3 features after PCAM} & \multicolumn{1}{c}{0.903} & \multicolumn{1}{c}{0.901} & \multicolumn{1}{c}{0.897} & \multicolumn{1}{c}{0.910} \\
   \multicolumn{1}{l|}{10 features after PCAM} & \multicolumn{1}{c}{0.907} & \multicolumn{1}{c}{0.903} & \multicolumn{1}{c}{0.919} & \multicolumn{1}{c}{0.900} \\
   \hline
  \multicolumn{1}{l}{\textbf{Freezing backbone weights}} \\
  \multicolumn{1}{l|}{Trainable backbone weights \textbf{(ours)}} & \multicolumn{1}{c}{\textbf{0.927}} & \multicolumn{1}{c}{0.932} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.928} \\
   \multicolumn{1}{l|}{Frozen backbone weights} & \multicolumn{1}{c}{0.919} & \multicolumn{1}{c}{0.923} & \multicolumn{1}{c}{0.919} & \multicolumn{1}{c}{0.915} \\
   \hline
  \multicolumn{1}{l}{\textbf{Image size}} \\
  \multicolumn{1}{l|}{512 $\times$ 512 \textbf{(ours)}} & \multicolumn{1}{c}{\textbf{0.927}} & \multicolumn{1}{c}{0.932} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.928} \\
   \multicolumn{1}{l|}{1024 $\times$ 1024} & \multicolumn{1}{c}{0.844} & \multicolumn{1}{c}{0.887} & \multicolumn{1}{c}{0.835} & \multicolumn{1}{c}{0.809} \\
   \hline
  \multicolumn{1}{l}{\textbf{Self-supervised pre-training}} \\
   \multicolumn{1}{l|}{ViT w/o pretrain} & \multicolumn{1}{c}{0.875} & \multicolumn{1}{c}{0.885} & \multicolumn{1}{c}{0.900} & \multicolumn{1}{c}{0.838} \\
  \multicolumn{1}{l|}{ViT w pretrain} & \multicolumn{1}{c}{0.879} & \multicolumn{1}{c}{0.872} & \multicolumn{1}{c}{0.875} & \multicolumn{1}{c}{0.889} \\
  \multicolumn{1}{l|}{ViT(hybrid) w/o pretrain} & \multicolumn{1}{c}{0.872} & \multicolumn{1}{c}{0.857} & \multicolumn{1}{c}{0.902} & \multicolumn{1}{c}{0.857} \\
  \multicolumn{1}{l|}{ViT(hybrid) w pretrain} & \multicolumn{1}{c}{0.866} & \multicolumn{1}{c}{0.865} & \multicolumn{1}{c}{0.890} & \multicolumn{1}{c}{0.841} \\
  \multicolumn{1}{l|}{Proposed w/o pretrain \textbf{(ours)}} & \multicolumn{1}{c}{0.927} & \multicolumn{1}{c}{0.932} & \multicolumn{1}{c}{0.921} & \multicolumn{1}{c}{0.928} \\
  \multicolumn{1}{l|}{Proposed w pretrain \textbf{(ours)}} & \multicolumn{1}{c}{\textbf{0.928}} & \multicolumn{1}{c}{0.951} & \multicolumn{1}{c}{0.918} & \multicolumn{1}{c}{0.915} \\
   \Xhline{1.2pt}
  \end{tabular}} \\
\end{table*}

\subsubsection{Network Size}
Since the network size can be various according to the number of Transformer encoder layers and self-attention heads, we evaluated the effect of network size by constructing models with different size consisting of 2 layers and 4 heads per layer, 4 layers and 8 heads per layer, 12 layers and 12 heads per layer. The experimental results shows that the standard Transformer architecture with 12 layers and 12 heads per layer works the best.

\subsubsection{Intermediate Feature Map}
In our model, the intermediate feature map of backbone network was used as input feature corpus for Transformer architecture. Here we tested the various options for intermediate feature maps. To determine whether the features before or after PCAM operations are better, the common feature map before PCAM operation and weighted feature maps after PCAM operation were compared. In addition, in case of weighted feature maps after PCAM operation, there might be options whether to use all 10 weighted feature maps for each CXR findings to utilize all relevant features, or to use weighted feature maps of CXR findings related to infectious disease (opacity, consolidation, and pneumonia), or just to use only the most related one feature (pneumonia) to reduce redundancy. As shown in Table \ref{table:ablation}, we can see that the common feature before PCAM operation achieved best performance.

\begin{figure*}[!t]
\centering
\includegraphics[width=0.96\textwidth]{Fig8.jpg}
\caption{Possible applications of our model in clinical scenario.}
\label{fig8}
\end{figure*}

\subsubsection{Freezing Backbone Weights}
Since our model utilizes the pre-trained backbone network with large-scaled dataset, we evaluated whether freezing or leaving the backbone weights trainable yield better result. The experimental result shows that it is better to train weights than to freeze them in all three external test data sets. This could be due to an increased capacity through a trainable backbone network, which dispels the fear of overfitting.

\subsubsection{Image Size}
Regarding the input image size, we compared the input CXR of size $1024 \times 1024$ and $512 \times 512$ with the same model and experimental settings. The size of intermeidate feature map ($16 \times 16 \times 768$) remained the same by increasing the stride of patch embedding layer by two fold. As suggested in Table \ref{table:ablation}, increasing the resolution of input CXR does rather deteriorate than improve the performance, leading to severe overfitting and poor generalization performance.

\subsubsection{Self-supervised Contrastive Pre-training}
Since previous works of Transformer-based models have suggested the possibilities that self-supervised pre-training to model structural sequence may be beneficial \citep{devlin2018bert, radford2018improving, chen2020generative}, we evaluated whether Transformer-based models (standard ViT, hybrid ViT and our model) benefit from self-supervsied pre-training. We adopted SimCLR \citep{chen2020simple}, a contrastive learning technique for visual representations with data augmentation framework, as our self-supervised pre-training method. As shown in Table \ref{table:ablation}, the benefit of self-supervised pre-training was not prominent for both hybrid ViT and our models, though it slightly improved the performance of standard ViT model. Nevertheless, our model still outperformed the other Transformer-based models pre-trained with self-supervised learning, which alludes that the low-level CXR feature corpus generated by the proposed method may be more suitable input feature embedding than direct pixel-patch embedding (standard ViT) or feature embedding by ResNet backbone (hybrid ViT) for CXR classification task.

\subsubsection{ROI Pooling Type}
 ROI max-pooling method is rooted in the annotation rule where 1 or 0 is assigned to each subdivision according to the presence or absence of the abnormality, not averaging each subdivision's opacity degree. Therefore, if there are pixels with a high probability of abnormality in a small region, although most of the pixels have a low probability, the value of 1 should be assigned into the subdivision. So, the max-pooling can be more appropriate for bridging the severity map and the array than average pooling. To demonstrate the hypothesis, we compared the performance of severity quantification between the model using average-pooling and max-pooling for converting the severity map into the array. As in Table \ref{table:ablation_sev}, the ROI max-pooling method outperforms the ROI average-pooling.

  \begin{table}[!h]
  \caption{Ablation studies for severity quantification.}
  \label{table:ablation_sev}
  \centering
  \resizebox{0.48\textwidth}{!}
  {\Large
  \begin{tabular}{l|c|c|c|c
  \Xhline{1.2pt}
\multicolumn{1}{l|}{\textbf{Methods}}
& \thead{\textbf{Avg. of}\\\textbf{3 External datasets}}
& \thead{\textbf{External dataset 1}\\\textbf{(CNUH)}}
& \thead{\textbf{External dataset 2}\\\textbf{(YNU)}}
& \thead{\textbf{External dataset 3}\\\textbf{(KNUH)}} \\ \hline

\multicolumn{1}{l}{\textbf{ROI Pooling}} & \multicolumn{1}{l}{}& \multicolumn{1}{l}{}& \multicolumn{1}{l}{}\\
\multicolumn{1}{l|}{Avg-pool}      & \multicolumn{1}{c|}{1.846}   & \multicolumn{1}{c}{1.842}  & \multicolumn{1}{c}{1.764}    & \multicolumn{1}{c}{1.933} \\
\multicolumn{1}{l|}{Max-pool \textbf{(ours)}}   & \multicolumn{1}{c|}{\textbf{1.655}}   & \multicolumn{1}{c}{1.682}  & \multicolumn{1}{c}{1.677}    & \multicolumn{1}{c}{1.607} \\
\hline
\multicolumn{1}{l}{\textbf{Training data}} & \multicolumn{1}{l}{}& \multicolumn{1}{l}{}& \multicolumn{1}{l}{} \\
\multicolumn{1}{l|}{Brixia only}          & \multicolumn{1}{c|}{2.503}   & \multicolumn{1}{c}{2.392}  & \multicolumn{1}{c}{3.310}    & \multicolumn{1}{c}{1.806} \\
\multicolumn{1}{l|}{Three datasets \textbf{(ours)}}       & \multicolumn{1}{c|}{\textbf{1.655}}   & \multicolumn{1}{c}{1.682}  & \multicolumn{1}{c}{1.677}    & \multicolumn{1}{c}{1.607} \\
\Xhline{1.2pt}
\multicolumn{4}{l}{\footnotesize{Mean squared error (MSE) of the severity score (0-6) is chosen as the metric for comparison.} \par}
  \end{tabular}}
\end{table}

\subsubsection{Training Data for Severity Quantification}
As Brixia dataset constitute the majority of data for severity quantification, we evaluated whether the model can show the stable performance only with the Brixia dataset. However, as in Table \ref{table:ablation_sev}, the performance was detrimental when the institutional datasets are all excluded from training data, addressing the necessity for including at least one institutional dataset to obtain stable performance of model. This degradation in performance may be result from a variety of reasons, including the discrepancy between details in labeling method of the Brixia dataset and the institutional datasets (e.g. anatomic landmarks determining the horizontal lines), the anatomical difference between races (e.g. Caucasian vs. Asian) and so on.

\section{Discussion and Conclusion}
\label{sec6}
In this study, we developed a novel ViT model that leverages low-level CXR feature corpus for diagnosis and severity quantification of COVID-19. The novelty of this work is to decouple the overall framework into two steps, the first is to pre-train the backbone network to classify low-level CXR findings with prebuilt large-scale dataset to embed optimal feature corpus, which was then leveraged in the second step by Transformer for high-level diagnosis of disease including COVID-19. By maximally utilizing the benefit of large-scale dataset containing more than 220,000 CXR images, overfitting problem of neural network with limited numbers of COVID-19 cases can be substantially alleviated. The extensive experimental results on various external institutions have demonstrated that our model not only outperforms other CNN and Transformer-based SOTA models but also retains outstanding generalization performance irrespective of the external settings, which is sine qua non of extensive adoption of system. In addition, we also adapted the proposed method to severity quantification problem,
demonstrating a performance similar to that of clinical experts, thereby expanding its application in the clinical setting.

In the current pandemic situation, our method holds great promise in a variety of clinical scenarios (see Fig. \ref{fig8}). Primarily, it can be used in patient triage along with the RT-PCR to isolate the suspected subjects waiting for RT-PCR results, as it was reported that positive radiological findings precede positive RT-PCR results in substantial portion (308 out of 1,014) of patients \citep{ai2020correlation}. In addition, it is possible to give guidance in treatment decision or to evaluate the response by applying our severity prediction algorithm to consecutive CXRs. Given the lack of experienced radiologists and unavailability of examination with higher sensitivity, the applications of our model could be of great value in underprivileged countries.

Finally, the concept of making higher-level diagnosis by aggregating low-level feature corpus, which is readily available with pre-built datasets, can be applied to quickly develop a robust algorithm against newly emerging pathogen, since it is expected to share the common low-level CXR features with existing diseases.

\section{Acknowledgement}
\label{sec7}
This work was supported in part by National Research Foundation of Korea under Grant NRF-2020R1A2B5B03001980, and in part by KAIST Mobile Clinic Module Project under Grant MCM-2021-N11210024.

\bibliographystyle{model2-names.bst}\biboptions{authoryear}


\section{Introduction}

Coherent X-ray source is one of the most important equipment for modern science. Its application covers various fields such as material science and biology \cite{gisriel2019membrane,xu2014single,alonso2020femtosecond,hirata2014determination,dean2016ultrafast}. Since the first hard X-ray free electron laser facility started its operation \cite{mcneil2009first}, extensive researches to provide better X-ray conditions for applications have been carried out. One of the research directions was generating two X-ray pulses having time separation or even energy separation, so called two-color X-ray \cite{PhysRevLett.84.2861,PhysRevLett.110.134801,PhysRevLett.113.254801,PhysRevAccelBeams.23.030703,hara2013two,marinelli2015high,PhysRevResearch.2.042018}. Recently these researches were extended to generate multi-color X-rays. So far, only a few different methods have been proposed \cite{lutman2016fresh,hemsing2019soft,qiang2011generation,xiang2012mode}, and further development is necessary to  explore possible performance improvements.

We introduce a new method to generate multi-color X-ray in this paper which uses beam's longitudinal phase space manipulation. This method consists of two key elements, longitudinal phase space (LPS) modulator and LPS manipulator. The modulator applies sinusoidal modulation on a single beam's LPS. Then, the manipulator controls this modulation to generate spectral and temporal bunching. Any modulation method can be used for this purpose, but here we assume to use wakefield structure. Also, we use double emittance exchange (EEX) beamline as the LPS manipulator.

The LPS modulator should provide a modulation of which wavelength is shorter than the bunch length. Here the number of periods within the bunch determines the number of bunches from temporal/spectral bunching. The modulation amplitude is one of the factors determining bunching level. In the case of wakefield structures, the modulation wavelength is determined by geometry and material \cite{PhysRevLett.61.2756,PhysRevSTAB.3.101302,BANE2012106}, and tunable frequency control is available for dielectric slab structures \cite{PhysRevE.56.4647,PhysRevLett.108.244801}. The modulation amplitude can be controlled by both the design of the structure and the charge level of the beam passing through the structure. If one uses collinear wakefield acceleration concept \cite{bane1985col,PhysRevSTAB.15.011301,PhysRevLett.120.114801}, the amplitude can be adjusted by controlling the charge level of a wakefield-driving beam.

In the case of the LPS manipulator, double EEX \cite{Ha:FEL2017-TUP054,zholents2011new} is one of possible options. The double EEX beamline enables the control of LPS via transverse manipulations. As shown in many references \cite{PhysRevAccelBeams.21.014401,seok2019suppression,seok2018sub,PhysRevSTAB.14.022801,Ha:IPAC2019-TUPGW089,PhysRevLett.105.234801,PhysRevAccelBeams.19.121301}, longitudinal properties can be manipulated using quadrupole magnets. This provides variable  ${R}_{55}$, ${R}_{56}$, ${R}_{65}$ and ${R}_{66}$,  which enables bunching of particles in both time and energy. It controls generated bunch's separations simultaneously.

The following sections provide the principle of the proposed method and experimental data showing its feasibility. Note that this experimental data is what we observed during other experiment. Thus, it does not show all aspects of this new method.

\section{Principle of new method}

\begin{figure*}[t]
\centerline{\includegraphics[width=0.5\textwidth,keepaspectratio=true]{fig2.eps}}
\caption{\textbf{Longitudinal phase spaces of numerically tracked beam.} An artificially generated beam travels through a 200 GHz dielectric structure and double EEX beamline to illustrate density and spectral bunching cases. Panel (a) shows the initial beam's LPS with uniform density distribution. Its full width is 6 mm to fit in four modulation periods. Panel (b) shows the LPS after the dielectric structure. Panels (c) and (d) correspond to spectral bunching and density bunching results, respectively. Particle tracking through a double EEX beamline is performed by using transfer matrix.}
\label{fig:stream}
\end{figure*}

In the case of conventional methods, they first introduce a sinusoidal modulation to beam's LPS. Then, proper $R_{56}$ from anisochromatic beamline such as chicane \cite{antipov-2013-a} converts this modulation to density modulation. Due to the energy difference that energy modulation introduced, particles having higher energy catch up preceding low energy particles via $R_{56}>0$. They move to zero-crossing points of the sine curve, and the density profile can be bunched when an appropriate $R_{56}$ is introduced. The new method we introduce uses the same principle but we use all $R_{55}$, $R_{56}$, $R_{65}$, and $R_{66}$ to generate either temporal or spectral bunching, and control their separation via EEX beamline. Equations given later in this section confirm this process and several simple but rough estimation of its characteristics. We assume a beam having zero longitudinal emittance with a linear chirp. Also, only its itnitial LPS was used to describe final LPS to simplify the discussion.

The first step of the method is applying sinusoidal modulation to the beam. When the beam passes through a high impedance medium such as dielectric tube, it generates wakefield. When the wakefield modulation frequency is high enough, this wakefield can introduce sinusoidal energy modulation to the beam. It can be written as,
\begin{equation} \label{eq_modulation}
\delta(z) = \frac{G_{m}L_{m}}{ E_0} \sin\left(\frac{\omega z}{c}\right) + Hz,
\end{equation}
where $G_m$ is the averaged gradient of the modulator, $L_m$ is the length of the modulator, $\omega$ is the angular frequency of the modulator's wakefield, $E_0$ is the reference beam energy, and $H$ is the initial longitudinal chirp that we call as macro-chirp. Because we know that the bunching happens near the zero-crossings, it is reasonable to observe only a half wavelength area nearby the zero-crossings (i.e. $-\frac{\lambda}{4} < z-n\lambda < +\frac{\lambda}{4}$) where the bunching happens. Then, Eq.~\eqref{eq_modulation} can be linearly approximated as,
\begin{equation}
\delta_n \simeq h\bar{z} + H(\bar{z}+n\lambda),
\end{equation}
where $|\bar{z}|<\frac{\lambda}{4}$ and $h \equiv \frac{\omega G_mL_m}{cE_0}$. We will call this $h$ as micro-chirp.

Next, this beam enters a double EEX beamline which has certain ${R}_{55}$, ${R}_{56}$, ${R}_{65}$ and ${R}_{66}$. Note that there will be quadrupole magnets in the middle of the double EEX beamline, and they can control these transfer matrix elements. Thus, the final LPS coordinates after the double EEX beamline become,
\begin{eqnarray}
z_{n,f} = R_{55}(\bar{z}+n\lambda) + R_{56}h\bar{z} + R_{56}H(\bar{z}+n\lambda), \\
\delta_{n,f} = R_{65}(\bar{z}+n\lambda) + R_{66}h\bar{z} + R_{66}H(z+n\lambda).
\end{eqnarray}
These simple equations tell us that longitudinal density can be bunched when the energy gain from modulator ($G_mL_m$) is equal to $-\frac{R_{55}+R_{56}H}{ R_{56}}\frac{cE_0}{\omega}$. Typically, the macro-chirp is much smaller than the micro-chirp, so it can be further simplified to $-\frac{R_{55}}{ R_{56}}\frac{cE_0}{\omega}$. Here, the bunch-to-bunch separation can be found from terms having no $\bar{z}$ dependency [i.e. $\lambda(R_{55}+R_{56}H)$]. Because both the bunching condition and the separation are determined by $R_{55}$ and $R_{56}$, quadrupole combinations are able to bunch the beam with an arbitrary separation. The same argument holds for the spectral bunching. Here, the bunching condition will be $-\frac{R_{65}+R_{66}H}{ R_{66}}\frac{cE_0}{\omega}$ and the energy deviation will be $\lambda (R_{65}+R_{66}H)E_0$.

This concept can be easily confirmed by tracking particles using transfer matrix (see FIG.~\ref{fig:stream}). For a numerical example, we assumed a single mode dielectric structure whose frequency is 200 GHz. Quadrupole magnets were assumed to be located in the middle of double EEX beamline to control ${R}_{55}$, ${R}_{56}$, ${R}_{65}$ and ${R}_{66}$. Their strengths were optimized to produce density and spectral bunchings. To apply the wakefield more realistically, the wake function of the structure was convoluted with beam's longitudinal density distribution. Two particle distributions having different charge were generated to represent wakefield-drive beam and target beam, and their charges were 2 nC and 200 pC, respectively. The beam energy of 50 MeV was assumed. The target beam had uniform longitudinal density distribution and its length was 4 times of the modulation wavelength (i.e., 6 mm); see FIG.~\ref{fig:stream}.

After the target bunch passes through the structure, it achieved a sinusoidal energy modulation as shown in the panel (b). This modulation is successfully converted to spectral and density bunchings as shown in panels (c) and (d).  The macro-chirp was zero, and the micro-chirp from the modulation was -10.98 m$^{-1}$. Quadrupole magnets were optimized to generate bunchings, and the double EEX beamline provided \textit{R}-terms as given in Table \ref{tab}. Because the macro-chirp is zero, the bunching condition was 11.08 m$^{-1}$ for both cases, which shows good agreement with the micro-chirp calculated from the distribution in FIG.~\ref{fig:stream}(b). Also, we expect that the density and the energy separations would be 2 mm and 100 keV, respectively. These also show good agreements with FIGS.~\ref{fig:stream}(c) and (d).

\begin{table}[]
    \caption{Matrix element used for numerical tracking. $R_{55}$ and $R_{56}$ are matrix elements for density bunching, and the other two elements for the density bunching case are not listed in this table. Similarly, $R_{65}$ and $R_{66}$ in this table are matrix elements the spectral bunching, and the other two elements are not listed.}
    \begin{ruledtabular}
    \begin{tabular}{c c c c}
    $R_{55}$ & $R_{56}$ & $R_{65}$ & $R_{66}$  \\
    \hline
    1.33 & 0.12 & 1.33 & 0.12
    \end{tabular}
    \label{tab}
\end{ruledtabular}
\end{table}

\section{Experimental data showing feasibility of the method}

Argonne wakefield accelerator facility recently performed an experiment to demonstrate the beam manipulation concept using a transverse wiggler (details can be found from \cite{Ha:IPAC2019-TUPGW089}). Although this experiment implies sinusoidal modulation to the transverse phase space and uses a single EEX beamline, its underlying principle is the same as the method introduced in this paper. A sinusoidal modulation is introduced into the phase space, and \textit{R}-terms of the EEX beamline are controlled by quadrupole magnets to bunch the beam and control its separation. While the main goal of this experiment is to demonstrate the concept of the transverse wiggler, they observed one interesting result shown in FIG.~\ref{fig:stream2} at certain quadrupole setting. Both LPS image and corresponding spectral profile show clear bunching in the energy axis. The beam was originally generated from a photocathode gun, and it had a smooth energy spectrum. However, the modulation introduced by the modulator and proper setting of the following EEX beamline provided spectral bunching as we expected from the previous section.  This is a clear evidence of the feasibility of the new method.

Note that the fractional energy separation in FIG.~\ref{fig:stream2} is about 1\%, and it corresponds to 400 keV. Each micro-bunch is not critically compressed. The energy spread of each micro-bunch can be further reduced by adjusting quadrupole setup. While the bunch is compressed in density in this data, the bunch length can also be controlled up to some level.

\begin{figure}[htb]
\centerline{\includegraphics[width=0.5\textwidth,keepaspectratio=true]{exp.eps}}
\caption{\textbf{Longitudinal phase space and corresponding spectral profile from transverse wiggler experiment.} Panel (a) shows measured longitudinal phase space after a double EEX beamline. Panel (b) shows spectral profile of the beam corresponding the panel (a).}
\label{fig:stream2}
\end{figure}

\section{Conclusions}
We introduced a new method to generate multi-energy beam that can be used to generate multi-color X-ray. This paper discussed the principle of the method, and one experimental data showing the feasibility of the method is provided. The strength of this method is its capability of bunching the single beam in both time and energy axes. Also, an arbitrary control of the separations is available.  The details of the method are not covered in this paper, but we can see that the separation of time and energy can be controlled simultaneously up to some level. We believe that this method can be one way to provide interesting beam configurations to light source users.


\section{Multi-Scale Unsupervised Regional Attentive Deep Feature}
\label{sec:method}

\begin{figure*}[!ht]
	\centering
	\includegraphics[width=0.9\linewidth]{design/MS-RMAC}
	\caption{Overall diagram of the proposed method. An input image with three different resolutions ($160\times160$, $224\times224$ and $320\times320$) is sent to a Convolutional Neural Network (CNN) to extract features from a convolutional layer. In this work, all features are extracted from ResNet50/Conv4 \cite{he2016deep}. These features are processed with the R-SMAC module and unsupervised regional attention module (URA). Finally, features are aggregated through the weighted sum. Weight factors are calculated with the URA module.}
	\label{fig:ms-rmac}
\end{figure*}

In this section, we introduce three modifications to the R-MAC (Regional-Maximum Activation of Convolutions) pipeline to improve its performance for TR. They are multi-resolution (MR), sum and max activation of convolution (SMAC), and unsupervised regional attention (URA). As shown in Fig. \ref{fig:ms-rmac}, MR is placed in the front of the R-MAC pipeline, while SMAC and URA are at the end. Before presenting details of these modifications, we briefly introduce the R-MAC pipeline.

R-MAC is an advanced version of MAC that is widely used as an aggregator for convolutional features. Several CBIR studies show that mid-level convolutional feature maps with maximum or sum outperform fully-connected layer features \cite{tolias2015particular,Jimenez_2017_BMVC,kalantidis2016cross,babenko2015aggregating}, as they contain more general features. However, MAC sacrifices local spatial information for compactness. For example, a convolution feature map of Image $I$ is $\mathcal{X}$. The shape of $\mathcal{X}$ is  $C (channel) \times W  (width) \times H (height)$. MAC will perform a spatial maximum pooling for each channel $\mathcal{X}_c$ of $\mathcal{X}$. Therefore, the MAC of $\mathcal{X}$ is,

\begin{equation}
	\mathbf{f} = [f_1, \dots, f_c, \dots, f_C], with f_c =  \max_{x\in\mathcal{X}_c} x.
\end{equation}

The MAC operation generates a compact representation $\mathbf{f}$ of size $C$, but it discards all information except maximum values. To retain important local information, Tolias \textit{et al.~} \cite{tolias2015particular} sampled multi-scale square regions from $\mathcal{X}$ in a sliding window fashion as shown in the ``region sampling" module of Fig. \ref{fig:ms-rmac}. The width of the sliding window is $2\times min(W,H)/(s+1), s = 1\dots S$ and its stride ($S$)is 60\% of it's width. The total number of regions, $N$, is decided by $S$. In our experiment, $S$ is set to 4, $N$, therefore, is $30$. Here, we use the notation $\mathbf{f_{R_i}} = [f_{R_i, 1}, \dots, f_{R_i, c}, \dots, f_{R_i, C}]$ for the MAC feature of the region $i$.

The R-MAC feature is the sum aggregation of $N$ regional MAC features. Usually, prior to sum aggregation, post-processing such as $\l_2$-normalisation and PCA-Whitening \cite{jegou2012negative} is applied. Here, we use the notation `` $\hat{}$ '' to represent this post-processing. The R-MAC feature thus defined as,
\begin{equation}
	\mathcal{F} = \sum_{i=1}^{N}\hat{\mathbf{f}}_{R_i}.
\end{equation}

In experiments, we learned the PCA-whitening on the sampled 30,000 trademarks, and the feature size after PCA is set to 256.

{\bf{Sum and Max Activation of Convolution (SMAC)}} MAC only encodes the maximum ``local'' response of each of the convolutional filters, which causes the loss of other important information embedded in the convolution features. To overcome this, we also apply sum pooling over regional features in addition to max pooling, and use the region-wise concatenation of them. This doubles the dimension size of the temporal regional features, however after post-processing the dimension size is reduced to 256.

{\bf{Multi Resolution (MR)}} Recent works \cite{gordo2017end,seddati2017towards} have improved R-MAC performance by using multiple resolutions. Inspired by these, we extract three convolutional feature maps for an input image $I$ with three resolutions ($160\times160$, $224\times224$ and $320\times320$). We therefore obtain $3\times N$ temporal regional features. The final R-MAC is the sum aggregation of these. Note all temporal regional features are post-processed. With this, $\mathcal{F}$ is equal to,

\begin{equation}
	\mathcal{F} = \sum_{j=1}^{3} \sum_{i=1}^{N}\hat{\mathbf{f}}_{R_{i,j}},
\end{equation}

where $\hat{\mathbf{f}}_{R_{i,j}}$ represents $i$th temporal regional feature of $j$th resolution.

{\bf{Unsupervised Regional Attention (URA)}} Kim \textit{et al.~} \cite{kim2018regional} note that R-MAC suffers from background clutter and varying importance of regions. They applied context-aware soft-attention to overcome this drawback of R-MAC. However, they generate the context-aware soft-attention signal via a regional attention network that is trained with ImageNet \cite{russakovsky2015imagenet}, and their method is intended for tasks where no labelled dataset is available for the target domain. However, they don't consider the context-difference between ImageNet and the target domain. For example, what is considered to be an important region in ImageNet might be unimportant for the target domain say TR. Therefore, we propose an unsupervised regional attention method to generate domain-specific context-aware soft attention. Our method learns the domain-specific context-aware soft-attention in an unsupervised manner from the gallery images that are available for every image retrieval task.

To model this domain-specific context-awareness, we applied a bag-of-words (BoW) model \cite{sivic2003video}. We build a regional deep feature dictionary by clustering sampled regional deep features extracted from the gallery images. K-means clustering (using the  FAISS library \cite{JDH17}) is applied to build the dictionary of 1,024 words. Noting the wide use of term-frequency (TF) and inverse-document-frequency (IDF) statistics \cite{salton1988term} for text-retrieval and context modelling, we calculate TF-IDF values for each word here to obtain a form of soft-attention. However, we note that TF is implicitly utilised during the sum aggregation step in R-MAC. Therefore, we only calculated the IDF for each regional deep feature,
\begin{equation}
	\mathcal{F} = \sum_{j=1}^{3}\sum_{i=1}^{N}\hat{\mathbf{f}}_{R_{i,j}}*IDF(\hat{\mathbf{f}}_{R_{i,j}}).
\end{equation}

In our experiments, we measure the similarity of two trademarks by calculating the Euclidean distance between their $l_2$ normalised modified R-MAC features.

\section{introduction}
\label{sec:intro}
A trademark (logo) is one of the most valuable intellectual properties of a company or individual. All trademarks require registration to avoid reputational and profit damages caused by trademark infringements. A trademark will be registered only if no duplication is found when it is compared with other registered trademarks. However, the exponential increase in the total number of trademark registrations and applications has made the registration process challenging. According to statistics reported by the world intellectual property office (WIPO), 11.5 million trademark applications were filed worldwide in 2019, which is a 5.8\% increase over 2018 \cite{WIPO}.

Large-scale trademark retrieval (LSTR) systems have been developed to detect and prevent trademark infringements. Early LSTR systems are text or code-based systems (\textit{i.e.} Vienna System), where each trademark is captioned by a human expert. Later, LSTR using content-based image retrieval (CBIR) algorithms have been used thanks to it's efficiency and accuracy. Hand-crafted features based-on shape, color or texture were developed for early CBIR-LSTR systems \cite{tursun2015metu,tursun2019componet}. With the rise of deep learning, off-the-shelf deep features have been applied for LSTR, demonstrating higher accuracy and efficiency compared to traditional hand-crafted features.

Due to the lack of publicly available labeled trademark datasets, recent studies focus on improving LSTR with deep features via post-processing off-the-shelf deep features, or pre-processing inputs rather than fine-tuning a pre-trained network. This type of image retrieval study is defined as ``pre-trained single-pass" by Zheng \textit{et al.~} \cite{zheng2017sift}, and our work also belongs to this category. The recent state-of-the-art (SOTA) LSTR study \cite{tursun2019componet} shows R-MAC (Regional-Maximum Activation of Convolutions) \cite{tolias2015particular} is efficient, simple and accurate compared to other post-processing techniques \cite{babenko2015aggregating,kalantidis2016cross,Jimenez_2017_BMVC}. However, its results suffer when text-components appear in trademarks as they increase the number of regions with text or background that are not essential for similarity detection. Previous work  \cite{tursun2019componet} improved R-MAC results by removing text-components that appeared in trademarks, and this approach can be considered as hard-attention. Kim \textit{et al.~} \cite{kim2018regional} also claim that R-MAC suffers in the presence of background clutter and regions of varying importance. They improve R-MAC performance by applying regional context-aware soft-attention for each regional MAC feature. However, they generate this attention via a regional attention network that requires supervised training. Although they trained the regional attention network on ImageNet \cite{russakovsky2015imagenet}, the domain difference between the target domain and ImageNet should be taken into consideration.

In this work, we also improve R-MAC's performance using soft-attention, where soft-attention is learned in an unsupervised manner on the gallery images. We apply a bag-of-words (BoW) model \cite{sivic2003video} to temporal regional features in the R-MAC pipeline, so each region is viewed as a word, and each trademark as a document. Later, we calculate the inverse document frequency (IDF) \cite{salton1988term} value for each word, which is used as soft-attention for regional features.

Additionally, we replace the MAC pooling adopted in R-MAC with a concatenation of sum and max pooling. The max pooling used by R-MAC only selects the spatially maximum activations of convolutional features, that results in a loss of other spatial information. We therefore integrate sum and max pooling. This extra step only increases the dimension of temporal regional features, however, the dimension of the final feature after aggregation remains the same as the feature dimension is reduced with a post-processing operation such as $l_2$ normalisation and PCA-whitening \cite{jegou2012negative}. Finally, to combat scale difference in TR, we introduce multi-resolution \cite{seddati2017towards,gordo2017end} inputs to the R-MAC pipeline.

We have tested our method on the challenging METU trademark dataset \cite{metudeeptursun}. All modifications show non-trivial improvements, and our final system surpasses the existing state-of-the-art methods by a large margin.

\section{Experiments}
\label{sec:exp}

\begin{figure*}[th!]
	\centering
	\includegraphics[width=0.9\linewidth]{images/qualitative/4}
	\includegraphics[width=0.9\linewidth]{images/qualitative/32}
	\includegraphics[width=0.9\linewidth]{images/qualitative/59}
	\caption{Top10 retrieved results from the METU dataset for example queries shown in the leftmost column.}
	\label{fig:qual}
\end{figure*}

\subsection{Dataset and Evaluation Protocol}
We select the METU trademark retrieval dataset \cite{metudeeptursun} as our testing dataset. It is the largest public dataset for TR. It includes nearly 1 million trademarks composed of text-only marks, figure-only marks and figure and text marks. Its evaluation set is composed of 35 similar groups, and each group includes around 10 to 14 similar trademarks. In total it includes 417 queries.

For evaluation, we follow the same evaluation protocol described in \cite{tursun2015metu}. In detail, we first return the ranking results for each query by sorting the similarity scores of the gallery images. We evaluated performance using the normalized average rank (NAR) and mean average precision (MAP) metrics. The NAR is calculated by normalizing the average ranking position of the ground-truth results of the queries. MAP values are only calculated for the top 100 results.

\subsection{Comparison with State-of-the-Art Methods}
We compare our multi-resolution R-MAC with unsupervised regional attention method with the recent SOTA TR methods that have been tested on the METU trademark dataset. Comparison results are shown in Table \ref{tab:sota}. Systems are categorized into three groups: hand-crafted features, fine-tuned off-the-shelf deep features, and pre-trained single pass \cite{zheng2017sift,kim2018regional}. Our method belongs to the last group as we don't fine-tune the feature extraction network. Our method achieves state-of-the-art NAR and MAP@100 results. In addition, our methods feature dimension is 256, which is the same as the previous SOTA.

\begin{table}[h!]
	\small
	\begin{center}
		\begin{tabular}{l c  c  c}
			\hline
			\bf Method & \bf DIM $\downarrow$ & \bf NAR $\downarrow$ & \bf MAP@100 $\uparrow$\\ \hline
			\multicolumn{4}{c}{\bf hand-crafted features} \\\hline
			Feng \textit{et al.~} \cite{feng2018aggregation} & 6,224 &  0.083 &  - \\
			Tursun \textit{et al.~} \cite{metudeeptursun}   & 10k   &  0.062 &  - \\\hline
			\multicolumn{4}{c}{\bf fine-tuned off-the-shelf deep features} \\\hline
			Perez  \textit{et al.~} (vis) \cite{perez2018trademark} & 4,096  & 0.066  & - \\
			Perez  \textit{et al.~} (con) \cite{perez2018trademark} & 4,096 & 0.063 & - \\
			Perez  \textit{et al.~} (vis, con) \cite{perez2018trademark} & 4,096  & 0.047 & - \\\hline
			\multicolumn{4}{c}{\bf pre-trained single pass \cite{zheng2017sift,kim2018regional}} \\\hline
			SPoC \cite{babenko2015aggregating,tursun2019componet} &  256 & 0.120 & 18.7\\
			CRoW \cite{kalantidis2016cross,tursun2019componet} 			& 256   &  0.140 & 19.8     \\
			R-MAC \cite{tolias2015particular}  		& 256 	& 0.072 & 24.8 \\
			MAC \cite{tolias2015particular,tursun2019componet} 	& 512  & 0.120 & 21.5 \\
			Jimenez \cite{Jimenez_2017_BMVC,tursun2019componet}   & 256 	& 0.093 & 21.0  \\
			CAM MAC \cite{tursun2019componet}  &  256  & 0.064   &22.3  \\
			ATR MAC \cite{tursun2019componet} &  512  & 0.056   &24.9 \\
			ATR R-MAC \cite{tursun2019componet} &	256  & 0.063	& 25.7 \\
			ATR CAM MAC \cite{tursun2019componet}	&  512  & 0.040   &  25.1 \\
			MR-R-MAC w/UAR (ours)  &  256  &  \textbf{0.028}  & \textbf{30.6} \\\hline
		\end{tabular}
	\end{center}
	\caption{Comparison with the previous state-of-the-art results on the METU dataset. NAR is the normalized average rank metric.}
	\label{tab:sota}
\end{table}

\begin{table}[h!]
	\small
	\begin{center}
		\begin{tabular}{l | c  c c }
			\hline
			\bf Method & \bf DIM  $\downarrow$ & \bf NAR $\downarrow$ & \bf MAP@100 $\uparrow$ \\ \hline
			R-MAC \cite{tursun2019componet}  &  256  & 0.051 & 27.3 \\
			R-SMAC \cite{tursun2019componet}  &  256 & 0.045 & 27.8 \\
			MR-R-MAC  \cite{tursun2019componet} &256 & 0.040 & 29.1 \\
			R-MAC w/URA & 256   & 0.039 & 29.6 \\
			MR-R-MAC w/URA & 256  & 0.033 &  30.6  \\
			MR-R-SMAC w/URA & 256  & \bf 0.028 &  \bf 31.0  \\
			\hline
		\end{tabular}
	\end{center}
	\label{tab:abl}
	\caption{Ablation study for the proposed approach.}
\end{table}

\subsection{Ablation Study}
We conduct an ablation study to consider the three modifications introduced to the R-MAC pipeline: multi-resolution (MR), sum and max pooling (R-SMAC), and unsupervised attention (UAR). We evaluate all possible combinations of the three modifications. Results are given in Table \ref{tab:abl}. All modifications bring non-trivial improvements. Moreover, the results show that they are complementary to one another. MR-R-SMAC with UAR (all three modifications) achieves the best result.

\subsection{Qualitative Results}
In Fig. \ref{fig:qual}, we visualized the top 10 retrieved results for three sampled queries. Here, all of the retrieved results are highly similar or related to the corresponding query. For the figure-only query displayed in the top row, figure-text marks are listed in the top 10 queries. On the other hand, for figure-text queries in the middle and bottom rows, no text-only marks are returned in the top 10 results. We therefore conclude that the proposed domain-specific context-aware unsupervised attention is beneficial.
\section{Related Studies}
\label{sec:lit}
The most recent trademark retrieval works are based on off-the-shelf deep features, as deep features are efficient and more accurate than hand-crafted features \cite{metudeeptursun,aker2017analyzing}. However, due to the lack of an annotated trademark dataset, most approaches are ``pre-trained single pass" methods. Early works \cite{metudeeptursun, aker2017analyzing} deploy pre-trained deep features from fully-connected layers for trademark retrieval. Later, aggregated deep-features from convolutional layers are studied \cite{lan2017similar,tursun2019componet}. Lan \textit{et al.~} \cite{lan2017similar} apply uniform local binary patterns (LBP) as an aggregator. They showed improved results compared to deep-features from full-connected layers, although the aggregation method is slow and returns features with high dimensionality. In comparison, Tursun \textit{et al.~} \cite{tursun2019componet} tested popular aggregation methods including SPoC, MAC, CRoW and R-MAC on the METU dataset. Deep convolutional features with these aggregation methods not only achieved improved results, but also reduced the feature size. To enable further improvements, Tursun \textit{et al.~} \cite{tursun2019componet} proposed soft and hard attention methods.

Few studies that fine-tune deep networks exist in the literature. Perez \textit{et al.~} \cite{perez2018trademark} improved deep feature performance by fine-tuning deep features with classification loss. To achieve this, they built a visual similarity dataset with 151 classes and a conceptual similarity dataset with 205 classes. Lan \textit{et al.~} \cite{lan2018similar} also improved deep feature performance by fine-tuning deep features with the triplet loss. Xia \textit{et al.~} \cite{xia2019trademark} proposed a transformation-invariant deep hashing method for efficient and transform-invariant trademark retrieval. However, the training sets of these works are either private or a subset of the METU query set. This makes their results difficult to fairly compare with the results presented here.

\section{Conclusion}
\label{sec:con}
We have introduced three modifications for the R-MAC pipeline to improve trademark retrieval performance: the use of multi-resolution input; domain-specific context-aware regional attention; and the use of both sum and max pooling for feature aggregation. With them, R-MAC features are more discriminant and robust to background-clutter and changes in scale. In the METU trademark dataset, the proposed method achieves a $5.3\%$ improvement over the state-of-the-art in MAP@100 metric and decreases the NAR score by 0.012.