# Assessing radiomics feature stability with simulated CT acquisitions

Kyriakos Flouris  Computer Vision Lab, ETH Zurich, Zurich, Switzerland
<EMAIL_ADDRESS>Oscar Jimenez-del-Toro  University of Applied
Sciences Western Switzerland (HES-SO) Valais, Sierre, Switzerland  Christoph
Aberle  Department of Radiology, University Hospital Basel, University of
Basel, Basel, Switzerland  Michael Bach  Department of Radiology, University
Hospital Basel, University of Basel, Basel, Switzerland  Roger Schaer  Faculty
of Medicine, University of Geneva (UNIGE), Geneva, Switzerland  Markus Obmann
Department of Radiology, University Hospital Basel, University of Basel,
Basel, Switzerland  Bram Stieltjes  Department of Radiology, University
Hospital Basel, University of Basel, Basel, Switzerland  Henning Müller
University of Applied Sciences Western Switzerland (HES-SO) Valais, Sierre,
Switzerland  Faculty of Medicine, University of Geneva (UNIGE), Geneva,
Switzerland  Adrien Depeursinge  University of Applied Sciences Western
Switzerland (HES-SO) Valais, Sierre, Switzerland  Department of Nuclear
Medicine and Molecular Imaging, Lausanne University Hospital, Lausanne,
Switzerland  Ender Konukoglu  Computer Vision Lab, ETH Zurich, Zurich,
Switzerland

###### Abstract

Medical imaging quantitative features had once disputable usefulness in
clinical studies. Nowadays, advancements in analysis techniques, for instance
through machine learning, have enabled quantitative features to be
progressively useful in diagnosis and research. Tissue characterisation is
improved via the “radiomics” features, whose extraction can be automated.
Despite the advances, stability of quantitative features remains an important
open problem. As features can be highly sensitive to variations of acquisition
details, it is not trivial to quantify stability and efficiently select stable
features. In this work, we develop and validate a Computed Tomography (CT)
simulator environment based on the publicly available ASTRA toolbox
(www.astra-toolbox.com). We show that the variability, stability and
discriminative power of the radiomics features extracted from the _virtual
phantom_ images generated by the simulator are similar to those observed in a
tandem phantom study. Additionally, we show that the variability is matched
between a multi-center phantom study and simulated results. Consequently, we
demonstrate that the simulator can be utilised to assess radiomics features’
stability and discriminative power.

##  1 Introduction

Computerized quantitative analysis of medical images is emerging as a
promising approach in radiological practice and healthcare research [1, 2, 3,
4]. These methods extract measurements quantifying various aspects of the
image that include basic intensity statistics as well as more complicated
metrics quantifying spatial intensity heterogeneity. Extracted measurements
are then used as image biomarkers in predicting relevant outcomes. In recent
years, numerous researchers demonstrated the capability of this approach for
diagnosis, stratification, and prognosis [5, 6]. Moreover, since the
extraction of measurements as well as the prediction stage are all
algorithmic, quantitative analysis is an efficient approach that can
complement radiologists’ visual interpretation and analysis.

Advanced artificial intelligence techniques [7], such as deep learning, take
the quantitative analysis approach one step further [8, 9]. They remove the
need to engineer measurements to extract from images for a given task.
Instead, they optimize their parameters to extract task-optimal measurements
and predict based on them. In the respective language, quantitative
measurements are called “features”. While the optimization requires large
number of data samples, i.e., training samples, if such large datasets exist,
deep learning algorithms can provide substantial accuracy gains [10].

An important limitation of the quantitative analysis approach is its
sensitivity to variations in scanning conditions [11]. While the methods aim
to extract measurements characterising the underlying tissue composition and
microstructure, they are indeed measurements taken from the image, which is
merely a representation of the tissue. Critically, image characteristics
heavily rely on the acquisition details, e.g., resolution, radiation dose,
noise, reconstruction algorithm. Depending on the properties of the algorithm
and the measurement, the extracted quantities can be highly sensitive to
variations in the image acquisition parameters [12, 13, 14]. This sensitivity
inhibits the generalisation capabilities of such measurements. If acquisition
details are not perfectly matched, two different images, even of the same
tissue, will yield different measurements. A number of studies have reported
the impact on CT radiomics analysis caused by the variability of acquisition
parameters and post-process variables [15, 16, 17, 18]. Any algorithm or
analysis based on these measurements will therefore not be reliable for use
with unseen scanners.

The ideal way to study the sensitivity of measurements is to perform test-
retest studies [19]. This would comprise of imaging a group of subjects imaged
under different acquisition details. To study sensitivity of a measurement,
values extracted from corresponding images would be compared. When new
measurements or new algorithms to extract measurements are proposed, they
would be studied the same way. As this is not feasible for various imaging
modalities, such as Computed Tomography (CT) due to the radiation exposure of
patients in these studies, anthropomorphic printed phantoms have been proposed
for CT variability studies [20, 21, 22].

Phantom studies have been successfully used for various imaging modalities.
Especially for CT, advances in 3D printing technologies allow printing
volumetric patient images using materials with attenuation properties
comparable to human tissue. Recent work reported variability studies using
such phantoms [23, 24, 25].

While phantoms make it possible to study variability without imaging cohorts,
they still require acquiring and imaging phantoms. This can be costly as well
as resource and time consuming. In this work, we study whether sensitivity
analysis using advanced in-silico CT simulators can yield similar results to
real phantom imaging studies. To this end, a CT-scan simulator environment was
set up using the publicly available [26, 27] ASTRA toolbox (www.astra-
toolbox.com). Using a high-dose CT-image as input, the simulator outputs raw
projections, which can be manipulated accordingly. For example, stationary and
uncorrelated noise can be added. Additionally, the simulator allows for some
freedom in geometrical parameters such as the number of projections, slice
thickness, and distances. The CT-image can be reconstructed with a variety of
algorithm choices, e.g. filtered back-projection and simultaneous iterative
reconstruction technique.

The method is compared with an empirical anthropomorphic phantom variability
study published in [23]. In this unique setup, the simulated phantom study is
performed using the same original image from which the anthropomorphic phantom
was printed and the study in [23] conducted. In a sense, this can be viewed as
the theoretical replication. The simulator environment was implemented to
reconstruct images at different noise levels, reconstruction algorithms, and
number of projections. To mimic repetition and introduce variability, each
simulation parameter set was repeated via a variation of the Poisson noise
random seed. For the simulated images, radiomics features were extracted and
analysed. As the same source image is used for both the empirical phantom
study and this work, direct comparison of the results of sensitivity analyses
is possible.

The next section describes the CT simulator environment method including a
brief introduction of the anthropomorphic phantom and the phantom study. In
the results section a comprehensive validation and comparison of the simulator
with respect to the phantom study is presented. Furthermore, a stability and
discriminative power analysis and discussions can be found in the same
section. The paper is summarised in the conclusions section.

##  2 Method

First, we introduce the details of the novel anthropomorphic phantom created
for the tandem phantom study [23]. A high dose CT-scan of this phantom is used
as the simulator input. Second, the extracted radiomics features, the
principal component analysis and the simulator environment are described in
detail.

###  2.1 Anthropomorphic phantom and phantom CT acquisitions

Here, we provide brief details of the anthropomorphic phantom study presented
in [23] for completeness. For further details, we refer the reader to the
original publication.

A realistic radio-opaque three-dimensional phantom was designed from real
patient CT data. Namely, the compilation of a half-mirrored lung including a
tumor and an abdominal liver section with a metastasis from a colon carcinoma
[23]. The phantom was manufactured via stacking sheets of printed aqueous
potassium iodide solution on paper [28]. The lung tumor section is a
replication of a publicly available patient data set for radiomics phantoms,
from the Image Biomarker Standardization Initiative [29]. The lung section was
neither used in this work nor the tandem phantom CT study. Tissue equivalent
attenuation at a defined energy spectrum was calibrated at 120 kVp. The
contrast resolution of the printing technique in the phantom goes from -100 to
1000 Hounsfield units (HU). Overall, no structures can be represented whose HU
is below this paper-induced threshold. To test the contrast resolution, a
circular intensity ramp was printed in the phantom running through an HU range
of 0 to 1000. A reliable resolution of 2 HU difference was achieved.
Consequently the abdominal region was adequately depicted for a quantitative
analysis within the printed HU range.

The phantom was imaged with a Siemens SOMATOM Definition Edge CT scanner
(SSDE). To define the acquisition and image reconstruction parameters, a
survey of clinical CT protocols was performed including 9 radiological
institutes. All the CT scans in that study were acquired with the same
acquisition parameters, which resulted in an approximate CT dose index of
10mGy. Namely, a tube voltage of 120 kVp, a helical pitch factor of 1.0, a 0.5
second rotation time, and a tube current time product of 147 mAs. No automatic
tube current modulation was used.

Typical reconstruction parameter settings for clinical protocols in thoracic
and abdominal oncology were varied for the phantom study as follows:
Reconstruction algorithm, iterative reconstruction (IR) or filtered back
projection (FBP); reconstruction kernel, 2 standard soft tissue kernels per
algorithm; slice thickness in millimeter, 1, 1.5, 2, 3; and slice spacing in
millimeter, 0.75, 1, and 2. Series reconstructed with an IR algorithm used an
ADMIRE (advanced modeled iterative reconstruction) at strength level 3. In
total, 8 groups of parameter variations were selected for the phantom study to
assess their impact on classic radiomics features. Initially, 20 repetition
scans were performed without re-positioning of the phantom, followed by 10
repetitions with re-positioning between each measurement. Therefore, 30
distinct acquisitions were performed for each of 8 parameter variation groups.

In the abdominal section six 3D regions of interest (ROIs) were manually
annotated by a board-certified radiologist using a thin-section phantom series
with 2 mm slice thickness and 1 mm spacing. The ROIs were annotated
conservatively, well within the margins, thus no cross-check step of the
annotations was performed by other radiologists. A polygonal outline was used
on all slices individually to define the ROIs. The six ROI binary masks were
stored in a 3D volume NIfTI format. Two normal liver tissue regions, two
cysts, a hemangioma, and a liver metastasis from a colon carcinoma were
included during the annotation process, regions can be found in Figure 1.
Further details of the annotated regions and the 8 variation groups can be
found in Jimenez-del-Toro et al. [23].

Figure 1: Annotated regions of interest on the anthropomorphic phantom.

A multi-center phantom CT study was also carried out with 13 different
scanners at selected locations in Switzerland. The scanners used were two
Siemens SOMATOM Definition Edge, two Siemens SOMATOM Definition Flash, a
Siemens SOMATOM Edge Plus, a Siemens SOMATOM X.Cite, a Philips Brilliance iCT,
a GE BrightSpeed S, a Philips Brilliance CT 64, a GE Revolution Evo, a GE
Revolution Apex, a Canon Aquilion Prime SP and a Canon Aquilion CXL. The same
protocol was implemented (as closely as possible) in all acquisitions. A tube
voltage of 120 kVp, a helical pitch factor of 1.0 and a 0.5 second rotation
time were used. The tube current time product was adjusted accordingly to
achieve the required dose of 10mGy. The IR reconstruction algorithm was used
with slice thicknesses 2 or 2.5 millimeter and slice spacing 1 or 1.25
millimeter.

###  2.2 Radiomics feature extraction and principal component analysis

From both the simulated and phantom CT scans, a total of 86 radiomics features
were extracted in 3D from the manually segmented ROIs using the open source
Pyradiomics python toolkit [29]. Definitions for the radiomics features are
available in the Pyradiomics documentation online
(https://pyradiomics.readthedocs.io/en/latest/features.html). The 86 features
extracted include 18 first-order statistics, 22 grey level co-occurrence
matrices, 14 grey level dependence matrices, 16 grey level run length matrices
and 16 grey level size zone matrices, as described briefly in the Appendix B.
Radiomics features parameters were set to their default values. More
specifically, no filter was applied to the input image and a fixed bin width
of 25 was used for the discretisation of the image grey level. Fixed bin size
discretisation is defined such that a new bin is assigned for every intensity
interval within the bin width starting at the lowest occurring intensity.
Additionally, no normalization, no spatial resampling, no resegmentation were
performed and no HU cutoffs were used within the ROIs for the extraction. The
distance between the center voxel and the neighbor, for which angles should be
generated, was set to one pixel. Furthermore, for the first order radiomics
the voxel array shift parameter was set to zero, for the grey level co-
occurrence matrices the co-occurrences was assessed in two directions per
angle, which results in a symmetrical matrix and for the grey level dependence
matrices no cutoff value for dependence was set, i.e. a neighbouring voxel was
always considered independent.

For the phantom CT acquisitions, an analysis was carried out via the principal
component analysis (PCA). The first two principal components of the 86
radiomics features from all 240 phantom CT acquisitions are shown in Figures 6
and 7 with black markers. The ROIs can be separated into 4 distinct tissue
classes, i.e. normal liver tissue, cyst, hemangioma, and liver metastasis. The
differences between the four ROI classes (inter-class variation) are larger
than all CT parameter variations (intra-class variation). ROIs from the normal
liver tissue class are closer in the feature space than those from the other
classes. All four classes remain linearly separable despite the CT parameter
variations.

Furthermore, the Wilcoxon statistic $W$ was used to assess the stability and
discriminative power of isolated radiomics features [23]. We set a threshold
of $W$ < 1 to indicate a stable comparison. The top 10 ranked features of the
phantom CT acquisitions are shown on the right-hand of the appendix Table 3.

###  2.3 CT simulator

The simulator environment was implemented to reconstruct images at different
noise levels, with different reconstruction algorithms, and number of
projections. Each simulation parameter set was repeated ten times for
different noise random seeds to approximate repeated scans. For the simulated
images, feature values were extracted and analysed. Specific features are
explained in detail in Section 2.2.

The ASTRA toolbox CT-scan and reconstruction simulator [26] was employed for
the purpose of this study. The simulator is based on simple geometric
principles for the creation of projection data (sinograms). These sinograms
can then be manipulated to mimic more realistic scenarios, for example through
adding Poisson noise. Subsequently, the processed images are passed to the
reconstruction algorithm. To match the simulator to the phantom acquisitions,
a helical scanning sequence of pitch one was realised by explicitly specifying
a sequence of helical projection vectors. These explicit projection vectors
define the scanning frequency, i.e. the total number of projections. A conical
beam is utilized and the target and detector are placed at 500mm and 1000mm
respectively to approximate the real scanner geometry. A flat square detector
of 512 by 512 of continuous pixels (1mm) was implemented for simplicity. The
number of detector pixels is higher than for a clinical CT scanner
(approximately 1000 by 64) but is nevertheless compensated by an equivalent
decrease in the scanning frequency, making the simulations simultaneously
efficient and realistic.

Random uncorrelated noise is added at the projection level by sampling from a
Poisson distribution,

$f(k;\lambda)=\frac{\lambda^{k}e^{-\lambda}}{k!},$

where, $f(k;\lambda)$ describes the probability of k occurrences and $\lambda$
is both the expectation and the variance of the distribution. A background
intensity $I_{0}$ is used to define the noise level i.e. at each pixel of the
projection images:

$\displaystyle I_{sampled}\sim f(\lambda=I_{0}e^{-I_{image}}),$ $\displaystyle
I_{final}=-\log(I_{sampled}/I_{0}).$

$I_{image},I_{sampled},I_{final}$ represent the initial image, sampled and
final intensities respectively. Therefore the background intensity is
inversely related to the Poisson noise. Here we denote the added noise level
$I_{0}^{-1}$ as $A$. The noise is added using the “add_noise_to_sino” function
in the ASTRA toolbox.

To calibrate an appropriate noise level $A$ the average pixel-wise variance
$\sigma^{2}$ is calculated for a range of $A$s and compared to the
$\sigma^{2}$ of the phantom CT acquisitions, see Figure 3. The average
$\sigma^{2}$ of the low dose (1 mGy) and high dose (10 mGy) acquisitions are
plotted as the horizontal lines. An approximate linear relation is observed
between $\sigma^{2}$ and $A$ as seen from the linear fit. The crossing points
between the horizontal line limits and the fitted line serve as a guide for a
realistic $A$ parameter range. In the simulation study, noise levels close to
the 10 mGy were used as this was the dose level used in the tandem phantom
study.

The reconstructions are performed with the simultaneous iterative "SIRT" and
filtered back-projection "FBP" 3D algorithms as implemented in the ASTRA
toolbox. Specifically the "$SIRT3DCUDA$" with 500 iterations and "$FDKCUDA$"
were used, the reconstruction kernels are fixed by the simulator and the slice
thickness is the same as the pixel resolution, i.e. 1 mm. Furthermore, a
distinct numerical random seed is used for the Poisson noise, to imitate
repetitions as performed for the phantom CT acquisitions [23]. The method is
very efficient numerically, as total computational time on a modern GPU is in
the order of minutes per complete reconstruction.

Figure 2: Axial views of anthropomorphic radio-opaque phantom. Left, original input. Middle, filtered back-projection reconstruction, right, iterative reconstruction, both obtained by the CT simulator.  Figure 3: Average pixel-wise variance of the iterative method simulated image plotted against the arbitrary noise measure. Black and grey lines denote the average variance of the high dose and low dose acquisitions. Parameter | range | optimal
---|---|---
Noise level ($\sigma^{2}/$pixel in HU2) | 2.5$\times 10^{-3}$ \- 2.8$\times 10^{-3}$ | 2.5$\times 10^{-3}$
Number of Projections | 150-450 | 450
Table 1:  Parameter choice for the simulation environment.

##  3 Results and Discussion

Pilot simulations are carried out with the optimal set of parameters as seen
in Table 1, i.e. minimum noise and maximum number of projections for ten
repetitions. First, the procedure is verified qualitatively by visual
inspection of the reconstructions with optimal parameter choice, axial
snapshots can be seen in Figure 2. Both reconstruction methods are
sufficiently successful. The iterative reconstruction has low noise and no
artifacts are visible. The FBP method is marginally noisier and exhibits some
minor artifacts, these differences are expected as theoretically the iterative
method is superior, albeit more computationally expensive.

Additionally, the Wilcoxon statistic $W$ is employed to analyze the stability
and discriminative power of the radiomics features as extracted from the
simulated CTs. To this end, a study is carried out to mimic the phantom CT
acquisitions. Namely the simulations are separated into 8 distinct groups with
different projection number and reconstruction algorithms, see Table 2. Within
each group, repetitions are achieved via a different Poisson noise random
seed. Across the study, the same noise level was added at the projection
stage. The ROIs are separated into 4 distinct tissue classes, i.e. normal
liver tissue, cyst, hemangioma, and liver metastasis. The analysis aims to
quantify stability and discriminative power of features across parameter
groups using the class definitions.

The result is depicted in Figure 4. The stability (intra-class variation)
percentage is calculated from a pairwise comparison among the 8 parameter
variation groups. This process is repeated for all available tissue classes,
while all other CT parameters are kept constant. Expressly, for each feature
from each class, $W$ is calculated in-between the groups. To this end, a
threshold of correlation is predetermined for $W$ at $1$. I.e. the repetitions
within the two tested groups in question follow the same distribution if $W<1$
and the pairwise comparison is considered successful. The percentage is
calculated as the total fraction of the successful pairwise comparisons for
each feature. The discriminative power (inter-class variation) is calculated
via pairwise comparison in-between tissue classes for each feature and group.
In contrast to stability, here a successful comparison is achieved if $W>1$.
Again the percentage represents the fraction of successful comparisons.

Group | Reconstruction | Projections
---|---|---
1 | SIRT | 150
2 | SIRT | 200
3 | SIRT | 250
4 | SIRT | 300
5 | FBP | 150
6 | FBP | 200
7 | FBP | 250
8 | FBP | 300
Table 2:  Parameter choice for the stability and discriminative power study.
For all simulations, the noise level was set to $A=0.0001$, i.e equivalent to
approximately 10mGy dose, and 10 different random seeds were used to achieve
repetitions within the group. Figure 4: Percentage stability of features as
intra-class comparison and discriminative power inter-class comparison.

The results show that although the majority of radiomics features had low
stability for CT parameter variations, as has been previously shown in other
studies [30, 31, 32, 33], the discriminative power is high in the task of
differentiating in-between the tissue classes. This relation is again observed
in the phantom study that is mimicked [23]. The top ten features across each
axis selected by the simulation environment, i.e., virtual phantom, and the
phantom CT acquisitions are compared in Table 3. To demonstrate the ability of
the simulation environment to predict stable features, an overlap of the best
scoring features relative the phantom CT acquisitions is plotted in Figure 5.
The $x$-axis represents an ascending percentage of features that are
considered as the highest scoring group (e.g. $10\%=$ top 9 out of 86
features) and the $y$-axis the percentage within that group that overlaps with
the top features seen in the phantom study, http://links.lww.com/RLI/A632. For
both stability and discriminate power, the overlap is consistently high, i.e.,
not in a linear relationship as expected for non-correlated lists.

Figure 5: Overlap of highest scoring features between simulation and phantom
CT acquisitions plotted against ascending percentage that are considered
highest scoring. Plotted for the stability and discriminative power alike. The
grey area represents un- or negatively-correlated overlap between the two
methods.

Furthermore, the radiomics features of the simulator are compared to the
empirical phantom acquisitions in Figure 6 in a variability analysis. To this
end, the principal components are calculated to investigate the similarity and
variability of the radiomics, and we use the parameter range as shown in Table
1. As seen from Figure 6, the simulation radiomics variability is in agreement
with the empirical results. It should be noted here that the study was carried
out in a semi-blind methodology, i.e. after matching all the possible
parameters to reality, the best possible values were used to create the
optimal reconstruction. Afterwards, an appropriate noise level was chosen
using Figure 3 for the purpose of this variability study.

Figure 6: Principal component analysis. The black markers indicate the
empirical phantom study data of the region of interest with the same shape as
the colored markers. The colored markers indicate the equivalent result of the
simulator with iterative reconstruction and the parameter range shown in Table
1. The average value and variability of the two principal components are
closely matched.

The filtered back-projection method creates an inferior reconstruction as seen
in Figures 2 and 7. There is a larger discrepancy between empirical
distributions and distributions obtained through the simulation. Nevertheless,
when the PCs are plotted the results indicate that the variance is well within
the experimental result. There is a lateral shifting of the first PC. The
variability is well captured by the simulator for all six ROIs.

Figure 7: FBP principal component analysis. The black markers indicate the
empirical phantom study data of the region of interest with the same shape as
the colored markers. The colored markers indicate the equivalent result of the
simulator with the filtered back-projection reconstruction and the parameter
range shown in Table 1. The average value and variability of the two principal
components are matched up to a shift of the first principal component.

Furthermore, the radiomics features of the simulator are compared to the
multi-center empirical phantom acquisitions in Figure 8 with a PCA variability
analysis. In the simulator the projection number is fixed to 200 and 250 and
the noise range extended to 2.5$\times 10^{-3}$ \- 2.9$\times 10^{-3}$
($\sigma^{2}/$pixel in HU2). This parameter range mimics the fixed slice
reconstruction thickness and the extended noise range was used to realise the
unknown differences inherent in a multi-center study. As seen from Figure 8,
the simulation radiomics variability is in agreement with the empirical
results.

Figure 8: Multi-center principal component analysis. The black markers
indicate the multi-center empirical phantom study data of the region of
interest with the same shape as the colored markers. The colored markers
indicate the equivalent result of the simulator with iterative reconstruction.

###  3.1 Discussion

Experimental comparison showed striking similarity between sensitivity
analyses carried out with the anthropomorphic phantom and the CT simulator.
Despite the approximations, the CT simulator was able to generate images with
very similar characteristics, as quantified by the features studied here, to
real images of the phantom. This is essentially a model whose parameters can
be changed to match those observed in phantom studies. This similarity opens
different avenues for further investigation and practical opportunities. For
example, studies with multi-centre, multi-vendor data sets [34] can be
effortlessly scaled up and automated.

First, the results suggest that sensitivity analysis for new features or new
ways to extract features can be initially performed with a CT simulator. This
would substantially reduce the efforts and costs required to study
generalization properties of new radiomics features, radiomics analyses and
image-based learning techniques to new acquisition settings. This is crucial
since this generalization ends up being one of the most notorious obstacles in
front of translating new quantitative image analysis technologies to clinical
practice.

While hand-crafted features’ stability can in theory still be quantified with
phantom studies, this approach remains very limited when it comes to assessing
stability of advanced algorithms that extract features in a data-driven way,
e.g., neural networks. Phantom studies yield very limited number of images and
this inhibits using them for assessing stability of neural network-based
feature extraction methods. The simulation study we showed here is a direct
solution to this issue. The approach can use any CT image as a "phantom",
therefore yields a large number of images to perform accurate stability
analysis of such advanced algorithms.

Second, training of learning-based methods can be modified to encourage
robustness to variations of imaging characteristics during training. For
instance, through extensive data augmentation one can gain robustness to
variations in Magnetic Resonance Imaging (MRI) acquisitions of the same
contrast [35]. As the CT simulator can generate images realistic enough to
yield similar sensitivity analyses as an empirical phantom study, one can
imagine using such simulators for training of highly robust deep leaning
models.

The CT simulator we used here did not consider various details of the
acquisition due to simplifications of the system’s physical model. Our
experiments with more complete models, such as Geant4 [36], showed that using
such models is challenging due to the difficulty in replicating a given scan
and computation time. Making more accurate simulators more user friendly and
faster may improve the quality of the sensitivity analyses. In addition, a
possible extension of this work can be the application of an automatic
segmentation method. Allowing for automated and accurate determination of
ROIs, especially useful for the segmentation of liver tissue regions.

##  4 Conclusions

Based on the astratoolbox we have created an environment to reproduce
artificial variability on an initial CT-image. The environment was verified to
replicate the diversification observed from empirical acquisitions via a
principal component analysis, both for intra- and inter-scanner analyses. The
methodology and simulational tool can accelerate the creation and testing of
stable and discriminative radiomics features. More crucially, this tool can
generate realistically variable CT-image datasets for training highly robust
deep learning models.

## References

  * [1] Kim, T.-Y., Son, J. & Kim, K.-G.  The recent progress in quantitative medical image analysis for computer aided diagnosis systems.  _Healthcare informatics research_ 17, 143–149, DOI: 10.4258/hir.2011.17.3.143 (2011).  22084808[pmid].
  * [2] Gillies, R. J., Kinahan, P. E. & Hricak, H.  Radiomics: Images Are More than Pictures, They Are Data.  _Radiology_ 278, 563–577 (2016).
  * [3] Zerunian, M. _et al._ Ct based radiomic approach on first line pembrolizumab in lung cancer.  _Scientific Reports_ 11, 6633, DOI: 10.1038/s41598-021-86113-5 (2021).
  * [4] Fu, Z. _et al._ CT features of COVID-19 patients with two consecutive negative RT-PCR tests after treatment.  _Scientific Reports_ 10, DOI: 10.1038/s41598-020-68509-x (2020).
  * [5] Liu, Z. _et al._ The applications of radiomics in precision diagnosis and treatment of oncology: Opportunities and challenges.  _Theranostics_ 9, 1303–1322, DOI: 10.7150/thno.30309 (2019).  30867832[pmid].
  * [6] Court, L. E., Rao, A. & Krishnan, S.  Radiomics in cancer diagnosis, cancer staging, and prediction of response to treatment.  _Translational Cancer Research_ 5 (2016).
  * [7] Kononenko, I.  Machine learning for medical diagnosis: history, state of the art and perspective.  _Artificial Intelligence in Medicine_ 23, 89–109, DOI: https://doi.org/10.1016/S0933-3657(01)00077-X (2001).
  * [8] Bermejo-Peláez, D., Ash, S. Y., Washko, G. R., Estépar, R. S. J. & Ledesma-Carbayo, M. J.  Classification of interstitial lung abnormality patterns with an ensemble of deep convolutional neural networks.  _Scientific Reports_ 10, DOI: 10.1038/s41598-019-56989-5 (2020).
  * [9] Hu, Q., Whitney, H. M. & Giger, M. L.  A deep learning methodology for improved breast cancer diagnosis using multiparametric MRI.  _Scientific Reports_ 10, DOI: 10.1038/s41598-020-67441-4 (2020).
  * [10] Kleppe, A. _et al._ Designing deep learning studies in cancer diagnostics.  _Nature Reviews Cancer_ 21, 199–211, DOI: 10.1038/s41568-020-00327-9 (2021).
  * [11] Yip, S. S. & Aerts, H. J.  Applications and limitations of radiomics.  _Phys Med Biol_ 61, R150–166 (2016).
  * [12] Prayer, F. _et al._ Variability of computed tomography radiomics features of fibrosing interstitial lung disease: A test-retest study.  _Methods_ 188, 98–104, DOI: https://doi.org/10.1016/j.ymeth.2020.08.007 (2021).  Artificial Intelligence Approaches for Imaging Biomarkers.
  * [13] Bae, Y.-K., Lee, J.-W. & Hong, S.  Effects of image distortion and hounsfield unit variations on radiation treatment plans: An extended field-of-view reconstruction in a large bore CT scanner.  _Scientific Reports_ 10, DOI: 10.1038/s41598-020-57422-y (2020).
  * [14] Mackin, D. _et al._ Measuring computed tomography scanner variability of radiomics features.  _Investigative Radiology_ 50, 757–765, DOI: 10.1097/rli.0000000000000180 (2015).
  * [15] Shafiq-Ul-Hassan, M. _et al._ Intrinsic dependencies of CT radiomic features on voxel size and number of gray levels.  _Med Phys_ 44, 1050–1062 (2017).
  * [16] Schmidt, R. M. _et al._ Assessment of CT to CBCT contour mapping for radiomic feature analysis in prostate cancer.  _Sci Rep_ 11, 22737 (2021).
  * [17] Balagurunathan, Y. _et al._ Reproducibility and Prognosis of Quantitative Features Extracted from CT Images.  _Transl Oncol_ 7, 72–87 (2014).
  * [18] Delgadillo, R. _et al._ Repeatability of CBCT radiomic features and their correlation with CT radiomic features for prostate cancer.  _Med Phys_ 48, 2386–2399 (2021).
  * [19] Fedorov, A. _et al._ An annotated test-retest collection of prostate multiparametric mri.  _Scientific Data_ 5, 180281, DOI: 10.1038/sdata.2018.281 (2018).
  * [20] Zhang, F. _et al._ Design and fabrication of a personalized anthropomorphic phantom using 3d printing and tissue equivalent materials.  _Quantitative Imaging in Medicine and Surgery_ 9 (2018).
  * [21] Gear, J. I. _et al._ Abdo-man: a 3d-printed anthropomorphic phantom for validating quantitative SIRT.  _EJNMMI Physics_ 3, DOI: 10.1186/s40658-016-0151-6 (2016).
  * [22] Irnstorfer, N., Unger, E., Hojreh, A. & Homolka, P.  An anthropomorphic phantom representing a prematurely born neonate for digital x-ray imaging using 3d printing: Proof of concept and comparison of image quality from different systems.  _Scientific Reports_ 9, DOI: 10.1038/s41598-019-50925-3 (2019).
  * [23] Jimenez-Del-Toro, O. _et al._ The Discriminative Power and Stability of Radiomics Features With Computed Tomography Variations: Task-Based Analysis in an Anthropomorphic 3D-Printed CT Phantom.  _Invest Radiol_ (2021).
  * [24] Cheng, C. P. & Halchenko, Y. O.  A new virtue of phantom MRI data: explaining variance in human participant data.  _F1000Research_ 9, 1131, DOI: 10.12688/f1000research.24544.1 (2020).
  * [25] Jahnke, P. _et al._ A radiopaque 3D printed, anthropomorphic phantom for simulation of CT-guided procedures.  _Eur Radiol_ 28, 4818–4823 (2018).
  * [26] van Aarle, W. _et al._ The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography.  _Ultramicroscopy_ 157, 35–47 (2015).
  * [27] van Aarle, W. _et al._ Fast and flexible x-ray tomography using the astra toolbox.  _Opt. Express_ 24, 25129–25147, DOI: 10.1364/OE.24.025129 (2016).
  * [28] Jahnke, P. _et al._ Radiopaque Three-dimensional Printing: A Method to Create Realistic CT Phantoms.  _Radiology_ 282, 569–575 (2017).
  * [29] Zwanenburg, A. _et al._ The image biomarker standardization initiative: Standardized quantitative radiomics for high-throughput image-based phenotyping.  _Radiology_ 295, 328–338, DOI: 10.1148/radiol.2020191145 (2020).  PMID: 32154773, https://doi.org/10.1148/radiol.2020191145.
  * [30] Berenguer, R. _et al._ Radiomics of CT features may be nonreproducible and redundant: Influence of CT acquisition parameters.  _Radiology_ 288, 407–415, DOI: 10.1148/radiol.2018172361 (2018).
  * [31] Traverso, A., Wee, L., Dekker, A. & Gillies, R.  Repeatability and reproducibility of radiomic features: A systematic review.  _International Journal of Radiation Oncology Biology Physics_ 102, 1143–1158, DOI: 10.1016/j.ijrobp.2018.05.053 (2018).
  * [32] Meyer, M. _et al._ Reproducibility of CT radiomic features within the same patient: Influence of radiation dose and CT reconstruction settings.  _Radiology_ 293, 583–591, DOI: 10.1148/radiol.2019190928 (2019).
  * [33] van Timmeren, J. E. _et al._ Test–retest data for radiomics feature stability analysis: Generalizable or study-specific?  _Tomography_ 2, 361–365, DOI: 10.18383/j.tom.2016.00208 (2016).
  * [34] Perkonigg, M. _et al._ Dynamic memory to alleviate catastrophic forgetting in continual learning with medical imaging.  _Nature Communications_ 12, 5678, DOI: 10.1038/s41467-021-25858-z (2021).
  * [35] Zhang, L. _et al._ Generalizing deep learning for medical image segmentation to unseen domains via deep stacked transformation.  _IEEE Transactions on Medical Imaging_ 39, 2531–2540, DOI: 10.1109/TMI.2020.2973595 (2020).
  * [36] Agostinelli, S. _et al._ Geant4—a simulation toolkit.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 506, 250–303, DOI: https://doi.org/10.1016/S0168-9002(03)01368-8 (2003).
  * [37] Zwanenburg, A. _et al._ The image biomarker standardization initiative: Standardized quantitative radiomics for high-throughput image-based phenotyping.  _Radiology_ 295, 328–338, DOI: 10.1148/radiol.2020191145 (2020).  PMID: 32154773, https://doi.org/10.1148/radiol.2020191145.

##  Appendix A Best performing radiomics comparison

Stability ($\%$) - simulator | Stability ($\%$) - empirical phantom study
---|---
original glszm SmallAreaLowgreyLevelEmphasis | original gldm LargeDependenceHighgreyLevelEmphasis
original gldm LargeDependenceHighgreyLevelEmphasis | original glszm SmallAreaLowgreyLevelEmphasis
original gldm SmallDependenceLowgreyLevelEmphasis | original firstorder Median
original firstorder Kurtosis | original glcm ClusterShade
original glcm ClusterShade | original gldm SmallDependenceLowgreyLevelEmphasis
original glcm InverseDifferenceMomentNormalised | original firstorder Mean
original glszm LowgreyLevelZoneEmphasis | original glcm InverseDifferenceMomentNormalised
original glrlm ShortRunLowgreyLevelEmphasis | original glrlm ShortRunLowgreyLevelEmphasis
original glrlm LongRunHighgreyLevelEmphasis | original glszm LowgreyLevelZoneEmphasis
original glcm InverseDifferenceNormalised | original glrlm LowgreyLevelRunEmphasis
Discriminative power ($\%$) - simulator | Discriminative power ($\%$) - empirical phantom study
original firstorder 90Percentile | original firstorder Median
original firstorder Energy | original firstorder Mean
original firstorder Mean | original glszm LargeAreaHighgreyLevelEmphasis
original firstorder Median | original firstorder Energy
original firstorder RootMeanSquared | original firstorder TotalEnergy
original firstorder TotalEnergy | original glszm greyLevelNonUniformity
original gldm DependenceNonUniformity | original firstorder Minimum
original glrlm RunEntropy | original glrlm greyLevelNonUniformity
original glszm greyLevelNonUniformity | original glszm LargeAreaLowgreyLevelEmphasis
original glszm LargeAreaEmphasis | original glszm SizeZoneNonUniformity
Table 3:  Comparison table of highest stability and discriminative power
radiomics features as observed in the phantom acquisitions and simulational
studies.

##  Appendix B Radiomics extraction description

The feature definitions are described in the Imaging Biomarker Standardization
Initiative [37]. More specifically, the radiomics feature categories used are,
https://pyradiomics.readthedocs.io/en/latest/index.html, :

  * •

18 first-order statistics which describe the distribution of voxel intensities
within a region through common metrics. For example,

$\textit{energy}=\displaystyle\sum^{N_{p}}_{i=1}{(\textbf{X}(i)})^{2},$ (1)

which is a measure of the magnitude (sum of the squares) of voxel values in an
image. X is a set of $N_{p}$ voxels included in the ROI.

  * •

22 grey level co-occurrence matrices (GLCM). A GLCM of size $N_{g}\times
N_{g}$ describes the second-order joint probability function of an image
region constrained by the mask and is defined as
$\textbf{P}(i,j|\delta,\theta)$. The $(i,j)^{\text{th}}$ element of this
matrix represents the number of times the combination of levels $i$ and $j$
occur in two pixels in the image, that are separated by a distance of $\delta$
pixels along angle $\theta$. The distance $\delta$ from the center voxel is
defined as the distance according to the infinity norm. In this work we
computed symmetrical GLCM. For example,

$\textit{autocorrelation}=\displaystyle\sum^{N_{g}}_{i=1}\displaystyle\sum^{N_{g}}_{j=1}{p(i,j)ij},$
(2)

which is a measure of the magnitude of the fineness and coarseness of texture.

$\textbf{P}(i,j)$ is the cooccurence matrix for an arbitrary $\delta$ and
$\theta$. $p(i,j)$ is the normalized cooccurence matrix and equal to
$\textbf{P}(i,j)/\sum{\textbf{P}(i,j)}$. $N_{g}$ is the number of discrete
intensity levels in the image.

  * •

14 Grey level dependence matrices (GLDM). A GLDM quantifies grey level
dependencies in an image. A grey level dependency is defined as the number of
connected voxels within distance $\delta$ that are dependent on the center
voxel. A neighbouring voxel with grey level $j$ is considered dependent on
center voxel with grey level $i$ if $|i-j|\leq\alpha$. In a grey level
dependence matrix $\textbf{P}(i,j)$ the $(i,j)$ th element describes the
number of times a voxel with grey level $i$ with $j$ dependent voxels in its
neighbourhood appears in image. For example:

$\textit{grey level
variance}=\displaystyle\sum^{N_{g}}_{i=1}\displaystyle\sum^{N_{d}}_{j=1}{p(i,j)(i-\mu)^{2}}\text{,
where
}\mu=\displaystyle\sum^{N_{g}}_{i=1}\displaystyle\sum^{N_{d}}_{j=1}{ip(i,j)}$
(3)

Measures the variance in grey level in the image. $N_{g}$ is the number of
discrete intensity values in the image. $N_{d}$ is the number of discrete
dependency sizes in the image. $N_{z}$ is the number of dependency zones in
the image. $\textbf{P}(i,j)$ is the dependence matrix. $p(i,j)$ is the
normalized dependence matrix, defined as $p(i,j)=\textbf{P}(i,j)/N_{z}$.

  * •

16 grey level run length matrices (GLRLM). A GLRLM quantifies grey level runs,
which are defined as the length in number of pixels, of consecutive pixels
that have the same grey level value. In a grey level run length matrix
$\textbf{P}(i,j|\theta)$, the $(i,j)^{\text{th}}$ element describes the number
of runs with grey level $i$ and length $j$ occur in the image (ROI) along
angle $\theta$. The value of a feature is calculated on the GLRLM for each
angle separately, after which the mean of these values is returned. For
example, grey Level Variance as above.

  * •

16 grey level size zone matrices (GLSZM). A GLSZM quantifies grey level zones
in an image. A grey level zone is defined as a the number of connected voxels
that share the same grey level intensity. A voxel is considered connected if
the distance is 1 according to the infinity norm (26-connected region in a
3D). In a grey level size zone matrix $P(i,j)$ the $(i,j)^{\text{th}}$ element
equals the number of zones with grey level $i$ and size $j$ appear in image.
Contrary to GLCM and GLRLM, the GLSZM is rotation independent, with only one
matrix calculated for all directions in the ROI. For example, grey Level
Variance as above.


Annual modulations of

the angular recoil–flux/energy distributions of

WIMP–scattered target nuclei

observed at an underground laboratory

Chung-Lin Shan

Preparatory Office of
the Supporting Center for
Taiwan Independent Researchers

P.O.BOX 21 National Yang Ming Chiao Tung University,
Hsinchu City 30099, Taiwan, R.O.C.

In this article,

we compare
the target and WIMP–mass dependent “annual modulations” of
the angular distributions of
the recoil flux and energy of
WIMP–scattered target nuclei
observed at different underground laboratories.

For readers' reference,

simulation plots
with different WIMP masses
and frequently used target nuclei
for all functionable underground laboratories

can be found and downloaded
on our online (interactive) demonstration webpage

§ INTRODUCTION

While most direct Dark Matter (DM) detection experiments
measure only recoil energies
of (elastic) WIMP–nucleus scattering events
deposited in underground laboratory detectors,

“directional” direct detection experiments
aim to provide 3-dimensional information
(recoil tracks and/or head–tail senses)
of WIMP–scattered target nuclei,

as a promising experimental strategy
for discriminating WIMP signals from backgrounds
[1, 2, 3]

and some incoming–direction–known astronomical events

the very low theoretically estimated event rate
and thus
only a few (tens) of WIMP scattering events
to be observed per year,

instead of
the originally proposed observation of
the diurnal modulation
in directional DM detection experiments,

in Refs. [5, 6]
we discussed and compared
the annual modulations of
the angular distributions of
the 3-D velocity (flux) and (average) kinetic energy of
Galactic halo WIMPs,
which impinge on our detectors,
observed at different underground laboratories.

After that,

in Refs. [7, 8]
we finally accomplished
our double Monte Carlo
scattering–by–scattering simulation package
for the 3-D elastic WIMP–nucleus scattering process

and can provide
the 3-D recoil direction and then the recoil energy of
the WIMP–scattered target nuclei
event by event
in different celestial coordinate systems.

in this article,

we will demonstrate and compare
(the target and WIMP–mass dependences of)
the annual modulations of
the angular distributions of
the recoil flux
and the average recoil energy of
scattered target nuclei
observed at different underground laboratories.

§ ANNUAL MODULATIONS OF
THE ANGULAR RECOIL–FLUX/ENERGY DISTRIBUTIONS
IN THE HORIZONTAL FRAME

By using our simulation package
described in detail
in Ref. [7],

in Figs. <ref>

we show
the angular distributions of
the recoil flux of
$\rmF$ nuclei
scattered by 100-GeV WIMPs
in four advanced seasons
[9, 7]
observed in the horizontal coordinate systems of
the SUPL,
the ,
the ,
the ,
and the ANDES

Interested readers can click each row
in Figs. <ref>
to <ref>
Figs. <ref>
to <ref>

to open the webpage of
the animated demonstration
for the corresponding annual modulation

(for more considered WIMP masses
and target nuclei
as well as
other underground laboratories).

5,000 experiments with
500 total accepted events on average
in one experiment
in each 60-day observation period of four advanced seasons
have been simulated
and binned into 12 $\times$ 12 bins.

in Figs. <ref>
(and Figs. <ref>),

we consider $\rmXe$
as one heavy target nucleus
(and raise the simulated mass of incident WIMPs
to $\mchi = 200$ GeV)

in order to demonstrate
the target (and the WIMP–mass) dependence(s) of
the annual modulation of
the angular recoil–flux distribution.

in Figs. <ref>,
to <ref>

we show
the angular distributions of
the average recoil energy of
$\rmF$ and $\rmXe$ target nuclei
scattered by 100- and 200-GeV WIMPs
in four advanced seasons
observed in the horizontal coordinate systems of
five considered laboratories,

The angular distributions of
the recoil flux of
$\rmF$ nuclei
scattered by 100-GeV WIMPs
in four advanced seasons
[9, 7]
observed in the horizontal coordinate systems of
the SUPL,
the ,
the ,
the ,
and the ANDES

As in Figs. <ref>,

except that

a heavy nucleus $\rmXe$
has been considered as our target.

As in Figs. <ref>:

$\rmXe$ has been considered as the target nucleus,

except that

the mass of incident WIMPs
has been considered as heavy as $\mchi = 200$ GeV.

The angular distributions of
the average recoil energy of
$\rmF$ nuclei
scattered by 100-GeV WIMPs
in four advanced seasons
observed in the horizontal coordinate systems of
five considered laboratories,

As in Figs. <ref>,

except that

a heavy nucleus $\rmXe$
has been considered as our target.

As in Figs. <ref>:

$\rmXe$ has been considered as the target nucleus,

except that

the mass of incident WIMPs
has been considered as heavy as $\mchi = 200$ GeV.

At first,
in each row (for each laboratory),

one can observe clearly
the annual variations of
the distribution patterns of
both of
the nuclear recoil flux
and the average recoil energy.

More precisely,

the all–sky average value of
the average recoil energy
(Figs. <ref>
to <ref>)
would be the largest (smallest)
in the advanced summer (winter)
for all considered laboratories
(with all three considered target–WIMP mass combinations).

by comparing
the distribution patterns of
both of
the nuclear recoil flux
and the average recoil energy
simulated with different target–WIMP mass combinations
(for each laboratory),

one can also find that

the heavier the target nucleus
and/or the mass of incident halo WIMPs,
the wider the distribution patterns,

the larger the most frequent/energetic recoil angles
(the angle between
the nuclear recoil
and the WIMP incoming velocity),

caused by
the stronger the cross section (nuclear form factor) suppression

theoretically speaking,

the 3-D distribution patterns of
the nuclear recoil flux
in particular,
the average recoil energy
would be useful
for constraining (or even identifying) the mass (range) of halo WIMPs.

The comparisons of
the angular distribution patterns of
the 3-D recoil flux and the (average) recoil energy
observed in the
horizontal coordinate systems of
each pair of laboratory groups
(Figure from Ref. [6]).

See the text
for detailed descriptions.

the LSC laboratory
in the European (Eu) group
as the comparison standard
(the middle row of
Figs. <ref>
to <ref>),

as summarized in Fig. <ref>,
it can be found that

* the distribution patterns
observed at Kamioka
(the East–Asian (EAs) group)

is similar
but 3-month earlier;

* the distributions
observed at SNOLAB
(the North–America (NAm) group)

is also similar
but 3-month later;

* the distributions
observed at SUPL
(the South–Pacific (SPa) group)

is similar to those
observed at Kamioka
but with a half–year difference
and a 180$^\circ$–degree rotation;

* the distributions
observed at ANDES
(the South–America (SAm) group)

is similar to those
observed at SNOLAB
but also with a half–year difference
and a 180$^\circ$–degree rotation;

* the distributions
observed at Kamioka and ANDES

as well as
observed at SUPL and SNOLAB

are similar
in the same seasons
but with a 180$^\circ$–degree rotation.

Note however that,

the angular distribution patterns of
the 3-D recoil flux and the (average) recoil energy
vary from laboratory to laboratory
(due to their geographical locations),

the all–sky average values of
the average recoil energy
(simulated with each target–WIMP mass combination)
measured at different laboratories
in the same season (observation period)
should be (almost) equal!

§ CONCLUSIONS

In this article,
we have demonstrated and compared
the (target and WIMP–mass dependent) annual modulations of
the angular recoil–flux/energy distributions of
WIMP–scattered target nuclei
observed at different underground laboratories.

this work could be useful for our colleagues
working on
directional direct Dark Matter detection physics.

§ ACKNOWLEDGMENTS

The author would like to thank
the pleasant atmosphere of
the W101 Ward and the Cancer Center of
the Kaohsiung Veterans General Hospital,

where part of this work was completed.

This work
was strongly encouraged by
the “Researchers working on
e.g. exploring the Universe or landing on the Moon
should not stay here but go abroad.” speech.

§ REFERENCES

[1]
S. Ahlen et al.
The Case for a Directional Dark Matter Detector and
the Status of Current Experimental Efforts
Int. J. Mod. Phys. A25 1–51
(Preprint 0911.0323)

[2]
F. Mayet et al.
A Review of the Discovery Reach of Directional Dark Matter Detection
Phys. Rept. 627 1–49
(Preprint 1602.03781)

[3]
S. E. Vahsen, C. A. J. O'Hare and D. Loomba
Directional Recoil Detection
Ann. Rev. Nucl. Part. Sci. 71 189–224
(Preprint 2102.04596)

[4]
M. Heikinheimo, N. Mirabolfathi, K. Nordlund, S. Sassi and K. Tuominen
Solar Neutrinos and Dark Matter Detection with Diurnal Modulation
Phys. Rev. D104 063037
(Preprint 2103.08511)

[5]
Shan C-L
Annual and Diurnal Modulations of
the Angular Distribution of the 3-D WIMP Velocity
Observed at an Underground Laboratory
J. Phys.: Conf. Ser. 1468 012017
(Preprint 1910.11763)

[6]
Shan C-L
Simulations of
the Angular Kinetic–Energy Distribution of
Halo Weakly Interacting Massive Particles
for Directional Dark Matter Detection Experiments
in publication

[7]
Shan C-L
Monte Carlo Scattering–by–Scattering Simulation of
3-Dimensional Elastic WIMP–Nucleus Scattering Events
(Preprint 2103.06485)
in publication

[8]
Shan C-L
Simulations of
the Angular Recoil–Energy Distribution of
WIMP–Scattered Target Nuclei
for Directional Dark Matter Detection Experiments
(Preprint 2103.06486)
in publication

[9]
Shan C-L
Simulations of
the 3-Dimensional Velocity Distribution of
Halo Weakly Interacting Massive Particles
for Directional Dark Matter Detection Experiments
(Preprint 1905.11279)
in publication


n$. In either case, RR3 never holds, and so these are unique geodesics.

$\mathcal{C}_{2}$: The case where $w=x^{j}$ is clear, so suppose
$w=x^{\mu_{1j}}y^{\epsilon_{1j}}\dots y^{\epsilon_{n_{j}j}}x^{j}$. Then
$\mu_{ik}\in[-(2k-j-1),j]\cap\mathbb{Z}_{\neq 0}$ for all $1\leq i\leq n_{j}$
. For RR3 to hold, we would need some $\sigma>0>\delta$ where
$\sigma,\delta\in[-(2k-j-1),j]\cap\mathbb{Z}_{\neq 0}$ such that
$\sigma+|\delta|=2k$. However if that was the case, then

$\sigma+|\delta|\leq j+2k-j-1=2k-1<2k.$

Hence RR3 never holds, and so geodesics are unique.

$\mathcal{C}_{3}$: consider $x^{\alpha}y^{\epsilon_{1}}x^{\nu_{1}}\dots
x^{\nu_{n}}\in\mathcal{C}_{3}$. We have
$\alpha,\nu_{i}\in[-(k-1),k-1]\cap\mathbb{Z}$ for all $1\leq i\leq n$. Hence
if RR3 holds for some $\sigma>0>\delta$, then

$\sigma+|\delta|\leq 2(k-1)=2k-2<2k.$

Hence RR3 can never be applied, and so geodesics are unique. This completes
the set $\mathcal{C}$ of unique geodesics.

For the set $\mathcal{D}$, assume $w$ is of the form in Eq. 6. If RR3 holds
for some $\sigma>0>\delta$ such that $\sigma\neq j$ or $\delta\neq-(2k-j)$,
then since $\sigma,\delta\in[-(2k-j-1),j-1]\cap\mathbb{Z}_{\neq 0}$, we have

$\sigma+|\delta|\leq j-1+2k-j-1=2k-2<2k,$

which is a contradiction. The only other possibility to apply RR3 is using any
$x^{j}$ and $x^{-(2k-j)}$ terms. Since $j+|-(2k-j)|=2k$, we can apply RR3 to
rewrite our word in the form

$ux^{j}vx^{-(2k-j)}w=_{G}ux^{-(2k-j)}vx^{j}w,$

for all pairs of $x^{j},x^{-(2k-j)}$ terms. Hence $\mathcal{D}$ consists of
non-unique geodesics. ∎

######  Remark A.3.

In the original classification of geodesics in [10], there is a fourth type of
geodesic defined as follows when $p=2k$:

  1. Type 4.

$\left(x^{\mu_{1k}}y^{\epsilon_{1}}\dots
x^{\mu_{nk}}y^{\epsilon_{n}}\right)x^{k}\left(y^{\epsilon}x^{\nu_{1}}y^{\delta_{1}}\dots
x^{\nu_{m}}y^{\delta_{m}}\right)x^{-k}\left(y^{\gamma_{1}}x^{\lambda_{1}}\dots
y^{\gamma_{l}}x^{\lambda_{l}}\right)$ where:

     * •

$n,m,l\geq 0$,

     * •

$\epsilon_{s},\epsilon,\delta_{t},\gamma_{u}\in\mathbb{Z}_{\neq 0}\;(1\leq
s\leq n,1\leq t\leq m,1\leq u\leq l)$,

     * •

$\mu_{sk}\in[-(k-1),k]\cap\mathbb{Z}_{\neq 0}\;(1\leq s\leq n)$,

     * •

$\nu_{t}\in[-(k-1),k-1]\cap\mathbb{Z}_{\neq 0}\;(1\leq t\leq m)$,

     * •

$\lambda_{u}\in[-k,k-1]\cap\mathbb{Z}_{\neq 0}\;(1\leq u\leq l)$.

Due to an error in the original paper, specifically $\mu_{sk}\neq-k$, Type 4
is in fact equivalent to our Type 2 geodesics as defined in 4.1, when $j=k$.
This can be seen by the fact that
$\nu_{t},\lambda_{u}\in[-k,k-1]\cap\mathbb{Z}_{\neq
0}\subset[-k,k]\cap\mathbb{Z}_{\neq 0}\;(1\leq t\leq m,1\leq u\leq l)$.

## Acknowledgments

The authors would like to thank Yago Antolín, Martin Edjvet, Thomas Haettel
and Fujii Michihiko for their generous advice and helpful discussions.

## References

  * [1] Usman Ali.  Conjugacy Classes of the 3-Braid Group.  Algebra Colloquium, 17(1):829–840, 2010.
  * [2] Yago Antolín and Laura Ciobanu.  Formal Conjugacy Growth in Acylindrically Hyperbolic Groups.  International Mathematics Research Notices, 2017(1):121–157, 2017\.
  * [3] Yago Antolín and Andrew W. Sale.  Permute and conjugate: the conjugacy problem in relatively hyperbolic groups.  Bulletin of the London Mathematical Society, 48(4):657–675, 2016\.
  * [4] Mladen Bestvina and Koji Fujiwara.  A Characterization of Higher Rank Symmetric Spaces Via Bounded Cohomology.  Geometric and Functional Analysis, 19(1):11–40, 2009.
  * [5] Martin R. Bridson and André Haefliger.  Metric Spaces of Non-Positive Curvature.  Springer, 1999.
  * [6] Laura Ciobanu, Alex Evetts, and Meng-Che Ho.  The conjugacy growth of the soluble Baumslag-Solitar groups.  New York Journal of Mathematics, 26:473–495, 2020.
  * [7] Laura Ciobanu, Susan Hermiller, Derek Holt, and Sarah Rees.  Conjugacy languages in groups.  Israel Journal of Mathematics, 211(1):311–347, 2016.
  * [8] Laura Ciobanu, Susan Hermiller, and Valentin Mercier.  Formal conjugacy growth in graph products I.  Groups, Geometry, and Dynamics, 17(2):427–457, 2023.
  * [9] Michel Coornaert.  Asymptotic growth of conjugacy classes in finitely generated free groups.  International Journal of Algebra and Computation, 15(05n06):887–892, 2005.
  * [10] M. Edjvet and D. L. Johnson.  The growth of certain amalgamated free products and HNN-extensions.  Journal of the Australian Mathematical Society. Series A. Pure Mathematics and Statistics, 52(3):285–298, 1992.
  * [11] Alex Evetts.  Rational growth in virtually abelian groups.  Illinois Journal of Mathematics, 63(4):513 – 549, 2019.
  * [12] Philippe Flajolet.  Analytic models and ambiguity of context-free languages.  Theoretical Computer Science, 49(2):283–309, 1987.
  * [13] Philippe Flajolet and Robert Sedgewick.  Analytic Combinatorics.  Cambridge University Press, 2009.
  * [14] Michihiko Fujii.  A new formula for the spherical growth series of an amalgamated free product of two infinite cyclic groups.  Kodai Mathematical Journal, 41(3):475–511, 2018.
  * [15] Ilya Gekhtman and Wen‐yuan Yang.  Counting conjugacy classes in groups with contracting elements.  Journal of Topology, 15(2):620–665, 2022.
  * [16] C. P. Gill.  Growth series of stem products of cyclic groups.  International Journal of Algebra and Computation, 09(01):1–30, 1999\.
  * [17] Victor Guba and Mark Sapir.  On the conjugacy growth functions of groups.  Illinois Journal of Mathematics, 54(1):301–313, 2010.
  * [18] Thomas Haettel.  XXL type Artin groups are CAT(0) and acylindrically hyperbolic.  Annales de l’Institut Fourier, 72(6):2541–2555, 2022.
  * [19] Derek F. Holt.  Garside groups have the falsification by fellow-traveller property.  Groups, Geometry, and Dynamics, 4(4):777–784, 2010.
  * [20] Derek F. Holt and Sarah Rees.  Artin groups of large type are shortlex automatic with regular geodesics.  Proceedings of the London Mathematical Society, 104(3):486–512, 2012.
  * [21] Derek F. Holt and Sarah Rees.  Conjugacy in Artin groups of extra-large type.  Journal of Algebra, 434:12–26, 2015.
  * [22] Michael Hull and Denis Osin.  Conjugacy growth of finitely generated groups.  Advances in Mathematics, 235:361–389, 2013.
  * [23] Roger C. Lyndon and Paul E. Schupp.  Combinatorial Group Theory.  Springer, 2001.
  * [24] Jean Mairesse and Frédéric Mathéus.  Growth series for artin groups of dihedral type.  International Journal of Algebra and Computation, 16(06):1087–1107, 2006.
  * [25] Valentin Mercier.  Conjugacy growth series of some wreath products, 2017.  http://arxiv.org/abs/1610.07868.
  * [26] Walter D. Neumann and Michael Shapiro.  Automatic structures, rational growth, and geometrically finite hyperbolic groups.  Inventiones mathematicae, 120(1):259–287, 1995.
  * [27] Igor Rivin.  Growth in free groups (and other stories)—twelve years later.  Illinois Journal of Mathematics, 54(1):327–370, 2010.
  * [28] Peijun Xu.  The genus of closed 3-braids.  Journal of Knot Theory and Its Ramifications, 01(03):303–326, 1992\.

DEPARTMENT OF MATHEMATICS, HERIOT-WATT UNIVERSITY, AND THE MAXWELL INSTITUTE
FOR MATHEMATICAL SCIENCES, EDINBURGH, EH14 4AS.

Email address<EMAIL_ADDRESS><EMAIL_ADDRESS>


# Video Representation Learning and Latent Concept Mining for Large-scale
Multi-label Video Classification

Po-Yao Huang, Ye Yuan, Zhenzhong Lan, Lu Jiang, Alexander G. Hauptmann
School of Computer Science, Carnegie Mellon University
{poyaoh, yey1, lanzhzh, lujiang<EMAIL_ADDRESS>
###### Abstract

We report on CMU Informedia Lab’s system used in Google’s YouTube 8 Million
Video Understanding Challenge. In this multi-label video classification task,
our pipeline achieved 84.675% and 84.662% GAP on our evaluation split and the
official test set. We attribute the good performance to three components: 1)
Refined video representation learning with residual links and hypercolumns 2)
Latent concept mining which captures interactions among concepts. 3) Learning
with temporal segments and weighted multi-model ensemble. We conduct
experiments to validate and analyze the contribution of our models. We also
share some unsuccessful trials leveraging conventional approaches such as
recurrent neural networks for video representation learning for this large-
scale video dataset. All the codes to reproduce our results are publicly
available111https://github.com/Martini09/informedia-yt8m-release.

##  1 Introduction

Ranging from booming personal video collections, surveillance recordings, and
professional video documentary archives, we have witnessed an unprecedented
growth of a wide range of video data. Numerous methods have been invented to
understand video contents and enable searching over huge volumes of
accumulated video data. Recently released large-scale video datasets such as
Google’s YouTube 8 Million (Youtube-8M) video collection bring advancements in
video understanding tasks and create new possibilities for many emerging
applications such as personalized assistant like Google Home and Microsoft
Cortana. Youtube 8M is a multi-label video classification benchmark composed
of pre-extracted Inception-v3 features [23], labels and their hierarchy in the
knowledge graph over more than 8 million videos.

The quantity makes Youtube-8M a unique video classification test-bed. There
are 5.7 millions training videos, 1.6 millions validation videos, and 0.8
testing videos respectively. The length of videos range from 120 to 500
seconds. Frame-level features are extracted under 1 frame per second (FPS)
sampling rate. Video level features are mean-pooled from frame-level features.
The size of topic theme pool is 4,716. Each video is with 3.4 labels on
average. In comparison to other weakly-labeled datasets [11], the precision is
reasonably good ($\sim 85\%$) while recall remains poor.

Learning an effective model for video understanding at this scale is
challenging for the three reasons: First, although effort in extracting
features at a scale of 8 million is alleviated, the provided frame-level
features are prepossessed with some unknown PCA and followed by a simple mean
pooling to generate video-level representation. We propose to learn an
attentive pooling kernel followed by a refined representation learning module
to further boost model performance. Second, labels (classes/concepts) are
assumed to be independent in the Youtube 8M dataset, which fails to capture
the authentic underlying relationship (such as co-occurrence, exclusion and
hierarchy) between concepts. We address this issue by learning and
incorporating latent concepts for multi-label classification. Third, multi-
model ensemble at this scale is under-explored. We design a systematic model
ensemble scheme and quantify the importance over heterogeneous models.

Our contribution in this paper is threefold: 1) We investigate feasible neural
architectures to enhance mixture of experts (MoE) model with refined
representation learning via residual links and hypercolumns. 2) We introduce a
novel latent concept learning layer to capture relationships among concepts 3)
We incorporate temporal segment data augmentation and leave-one-out ensemble
to further boost classification accuracy.

##  2 Related work

In this section, we briefly discuss the most related previous work falling
into two categories: how to train a better individual deep video
classification model and how to fuse these individual models to achieve best
accuracy.

At the risk of oversimplification, we identify two classes of methods to
improve deep neural networks for video classification. The first class of
methods is having more data. Generally, there are two ways that deep models
can be benefit from having more data. The first approach is supervised pre-
training with external data as in [24, 25]. This approach is especially useful
when the inputs are raw video frames. However, this direction is not feasible
for the Youtube-8M dataset since we can only access pre-processed video
features. The second approach to have more data is through augmenting internal
data [24]. We design a simple way to augment data by exploiting the fact that
video information is highly redundant and the boundary of motion is often
arbitrary. In addition to apply augmentation for training, we further extend
this approach at the inference phase.

The second classes of methods is designing better network structure. Likewise,
this class of methods can be characterized into two groups: task-independent
and task-dependent. There are numerous improvements that are task-independent.
For example, dropout [20], inception structure [22], residual structure [8],
just name a few. In this work, we explore the residual structure [8] and its
variant [4] and found that these structures not only can help to learn deeper
networks but can also improve shallow networks. There are also improvements
that are more specific to video classification tasks. Most of these
improvements try to capture temporal dependencies among frames. The general
lesson is that short-term dependency are useful and easy to capture, but long-
terms ones are much more difficult to capture. For example, Simonyan et al.
[19] design a two-stream architecture to capture variances between consecutive
frames and there are significant amount of following works that try to further
improve the network structures [24] or capture longer-term temporal
information [6, 25, 18, 13, 15]. These work represent the state-of-the-art for
video classification and they find that sequence model such as LSTM often
perform worse than simple BoF models. These observations are consistent with
ours on the Youtube-8M dataset.

Since the capacity of a single model is limited, researchers often use an
ensemble of multiple models to improve video classification accuracy. Multiple
models can be learned jointly or separately followed by fusion. MoE [10] and
dropout [20] both learn a large number of single models jointly. However,
these single models are often need to be homogeneous and the size of the model
is also limited. To explore the diverse characteristics of heterogeneous
models, we train a set of different MoEs separately and fuse them using leave-
one-out [14] fusion method. Another reason for learning different models is
that different modalities of data may be harder to learn at the same time.
Luckily, we find that concatenating visual data and audio data together and
feed them into the same network (early fusion) is better than learning them
separately and fuse the results later (late fusion) for the Youtube-8M
dataset.

Figure 1: Video-level and Frame-level models

##  3 Models

Training video-level and frame-level models differ a lot in cost and
performance. Without losing choices for practical system design, in this
paper, we start with improving video-level models, which is more cost
effective, and then we treat frame-level models as models with an additional
encoding module that can be stacked upon video-level methods.

###  3.1 Video Level Models

Built upon MoE models in [2] without the null experts, we propose to improve
the model accuracy in the multi-label video classification task by refining
the representation learning with residual links and hypercolumns.

Residual Learning: While growing MoE deeper, results show that performance
gain is mediocre under the same model complexity. As shown in Fig. 1, we then
apply residual links to the middle layers of expert part of MoE to learn
better representation for multi-label classification. We name the combined
model Mixture of Residual Experts (MoRE). Residual learning is originally
designed to relieve gradient vanishing problem when training very deep neural
networks [8]. However, we find that simply adding identity mapping to the
shallow MoE network also brings significant improvement. Formally, the
residual expert network learns the refined representation
$\mathbf{x}^{\prime}$ from the original input $\mathbf{x}$ with:

$\mathbf{x}^{\prime}=F(\mathbf{x},{W_{i}})+W_{s}\mathbf{x}$ (1)

where in practice $F()$ is stack of a single layer perceptron followed by a
batch-normalization layer and dropout layer. We choose a identity matrix for
$W_{s}$. Note that the mixture network which determines the weights of
individual classes over experts still take original input instead of the
learned representation. We find that using deep structure for the mixture
network harm the performance with or without residual links.

Hypercolumn: Similar to MoRE, we also design Mixture of Hypercolumn Experts
(MoHCE). In a typical CNN, higher convolutional layers can capture high-level
global context. But they could miss low-level details, thus numerous
approaches have built predictors based on the concatenation of multi-stage
features [4] [16]. We propose to formulate a hypercolumn by concatenating
features from different layers of a MLP. Formally, the final features are
given by:

$\mathbf{x^{\prime}}=[F_{1}(\mathbf{x}),F_{2}(\mathbf{x}),...,F_{n}(\mathbf{x})]$
(2)

where $F_{i}(\mathbf{x})$ denotes the features of $i$-th layer. This
hypercolumn is then fed into mixture of experts (MoE) models to produce final
predictions. For MoHCE, we also use standard batch-normalization and dropout
techniques to regularize the network training.

###  3.2 Frame Level Models

Recurrent Neural Network Models (LSTMs). Videos as a collection of frames are
inherently rich with temporal information. Although many efforts have been
made in training recurrent neural networks (RNNs) for video classification,
finding a feasible representation remains an open question since its
challenging to 1) regularize RNNs and 2) pool representative information frame
by frame. Viewing video $v$ as a sequence of frame-level features
$x_{1:F_{v}}^{v}$, where $x_{j}^{v}$ is the features on $j$-th frame, the
video-level representation learning from frame-level features with RNN is:

$\mathbf{x}^{\prime}=P(\\{\Gamma(\mathbf{x}_{t},\mathbf{h}_{t-k})\\}_{k=1...K,t=1...T}$
(3)

where $P()$ is the pre-defined or learned pooling kernel. $\Gamma()$ is the
RNN function with $t$ as the time index and $k$ as the recurrence time index.
For example, in [2], the authors use long short term memory (LSTM) network for
$\Gamma()$ and the pooling function is a simple concatenation of the last
output and the hidden state $\mathbf{x}=o_{-1}||h_{-1}$. In this work we
compare performance of different RNNs and RNNs with the learned pooling kernel
by attending to individual frames. On the other hand, since videos usually
contain a lengthy sequence of frames, we also experiment regularization
methods for RNN such as dropout, zone-out [12], layer-wise batch-normalization
[5] to avoid the vanishing gradient problem.

Intriguingly, experimental results shows that these models with recurrent link
do not outperform simple attentive pooling for the PCA-ed features in the
Youtube 8M dataset.

Attentive Pooling Models (Deep BoF, NetVLAD): Another direction to exploit
temporal information is to use attentive pooling of the encoded input
features. We propose to remove the RNN module $\Gamma()$ from Eq. 3 and
directly learn a pooling function
$P(\\{\mathbf{x}_{t}\\})=\sum_{t}{\alpha_{t}E(\mathbf{x}_{t})}$, where $E(.)$
is some learnable encoding/transformation and $\alpha_{t}$ is the (attentive)
pooling weight. There are many possible design choices of $E(.)$. One example
is sparse encoding for deep bag-of-Frames (DBoF) in [2], where mean pooling
has been utilized to aggregate the encoded frames into video-level
representation.

A simple yet effective improvement can be achieved is to make the pooling
function learnable. As generalized in NetVLAD [3], the pooling weights are
calculated with a trained soft-attention from frames inputs to some given
cluster center. In our frame-level models, we enhance original NetVLAD with an
additional multi-layer perceptron for sparse encoding. Formally, the attentive
pooling weights are calculated as:

$\alpha_{tk}=\frac{\mathbf{w_{k}}^{T}E(\mathbf{x}_{t})}{\sum_{i}^{K}\mathbf{w_{k}}^{T}E(\mathbf{x}_{i})}$
(4)

where $K$ is number of clusters. To make the representation more condense, our
model learn another dense layer $G(.)$ to reduce the dimension after
concatenation. The final video representation
$\mathbf{x}^{\prime}\in\mathbb{R}^{d}$becomes:

$\mathbf{x}^{\prime}=G(\\{\sum_{t}\alpha_{t1}\\{E(\mathbf{x}_{t})-c_{1}\\}||...||\sum_{t}\alpha_{tK}\\{E(\mathbf{x}_{t})-c_{K}\\})$
(5)

###  3.3 Latent Concept Learning

Classes (or semantic concepts as the more general term) may share complex
relationships with each other. As can be seen from the Youtube-8M
EDA222https://www.kaggle.com/philschmidt/youtube8m-eda, one concept may
frequently co-occur with another concept. For example, chairs and tables can
usually be seen together in the indoor-scene videos. Additionally,
hierarchical structure of the concepts may also be encoded within the
relationships between concepts. For instance, ”iPhone6”, ”iPhone5” belong to
the category ”iPhone”. One approach to disentangle the hierarchical structure
among concepts is to map concept names using hand-crafted knowledge
graph333https://developers.google.com/knowledge-graph/. Unfortunately, the
direct alignments between video labels and the knowledge graph entities
remains unclear, making the built-upon graphical structure un-grounded. To
address these issues, we propose appending an additional deep neural network
at the output of MoRE to directly capture the relationship between concepts.
We also show that incorporating these latent concepts would improve the multi-
label classification performance.

Unlike the approach in [21] which uses an additional layer to absorb label
noise by twiddling the regularization term, we propose to learn the latent
concepts represented in the middle layer of deep residual networks directly
and incorporate the original MoRE outputs for final classification.
Specifically, we add an additional 2-layered residual network (we called it
the latent concept layer, $LC$-layer) with batch normalization as well as
input dropout at the output $\mathbf{y}$ of MoRE:

$\mathbf{y}^{\prime}=P(G(F(\mathbf{y},W_{i}^{LC}),W_{s}^{LC}))$ (6)

where $y^{\prime},y\in\mathbb{R}^{4716}$, $W_{i}^{LC}$ and $W_{s}^{LC}$ are
the $LC$-layer parameters, $G()$ is the aggregation function and $P()$ is the
final output layer with sigmoid function for multi-label classification. To
further alleviate error propagation of ill-classified concepts at the early
training stage, we first train without learning latent concepts for 10 epochs
and then append a randomly initialize $LC$-layer to start latent concept
learning. We find that this late-fire strategy is critical to mine reasonable
latent concepts and improve final classification performance. For the
aggregation function $G()$, we explore $add()$, $max()$, $append()$ and find
that $add()$ delivers the best performance.

###  3.4 Ensemble

We use leave-one-out method [14] to determine the fusion weights of individual
models. For a defined evaluation metric $s$ and model set $M$, the fusion
weight $w_{m},\sum_{m}w_{m}=1$ for model $m\in M$ is proportional to
$d_{m}=s_{M-m}-s_{M}$, which is the performance drop without model $m$ in
comparison to the baseline performance $s_{M}$. A typical choice of the
baseline is the result fusing all models equally. In the Youtube-8M video
understanding challenge, we use Global Average Precision (GAP) (Details in
Section 4) to calculate the weights.

##  4 Experiments

###  4.1 Training and Evaluation

Features: In the Youtube-8M dataset, raw visual features are extracted from
Google’s Inception-v3 model trained on Imagenet [23]. Raw audio features are
extracted from a CNN-inspired architecture trained for audio classification as
described in [9]. Both visual and audio features follow a PCA whitening
process (the PCA matrix is unknown) to further reduce the dimension to 1,024
and 128 respectively. The video-level features are mean-pooled from frame-
level features. For training we use the original training set and
$\frac{15}{16}$ of the validation set. We evaluated our models on the excluded
$\frac{1}{16}$ validation split. The difference in GAP between our
$\frac{1}{16}$ validation split and the real test set is less than 0.02%.

Model Training Details. All our models are trained with cross entropy loss
using the Adagrad algorithm [7] with a 0.0002 learning rate and a 0.8
exponential decrease over 8 million examples. We also experimented different
loss function such as weighted label loss as in [17], optimizer such as
momentum and RMSProp but didn’t found them beneficial for convergence speed.
Common techniques for training neural models including gradient clipping
(0.8), dropout (0.8 keep rate), and batch normalization are applied. The mini-
batch size is 1,024 (videos) for video level models and 256 (frames) for frame
level models.

For video level models, the input is a simple concatenation of l2-normalized
mean_RGB and mean_audio features in the original dataset. We vary the size of
number of mixtures while setting the unit of hidden layer in the residual
expert 4,096. Deep MoE models share the same amount of parameter as MoRE but
are without residual links. The number of latent concepts we learned in the
$LC$-layer is 4,096. The models are trained with 80,000 and 160,000 steps for
extended training split and for extended training split with data augmentation
(DA) respectively. For frame-level models (DBoF, NetVLAD), we set size of
sparse coding 4,096, sample from 80% of frames, and then compress final
representation into a 2,400 dimension vector. The dimension is also 2,400 for
LSTM models (dimension for Bi-directional LSTM is doubled). The typical step
for convergence is 360,000.

Our code is based on Tensorflow [1]. We use AWS p2-2xlarge instances (Intel
Xeon E5-2670 CPU and NVIDIA K80 GPU) to train our models. For most models, we
experiment hyper parameters of each model such that each model can maximally
fit into 12 GB GPU memory444MoRE24 and MoHCE12 share roughly the same
parameters. MoRE28 and MoRE30 can also fit but the training is sometimes
unstable. For frame level models, we use MoRE8 as the classification model.
There is a notable difference in cost of storage and computation for training
video level and frame level models. Training video level models takes around
150 GB storage and around 8 hours to converge while frame-level models takes
around 2TB and roughly 5 days to converge. The cost (computation and storage)
ratio between a video-level and a frame-level model is roughly 1:16.

Evaluation Metric. We report our results in three metrics: Mean Average
Precision (mAP) and Precision at equal recall rate (PERR) [2] and Global
Average Precision (GAP). The GAP metric takes the predicted labels that have
the highest k (k = 20) confidence scores for each video, then treats each
prediction as an individual data point in a long list of global predictions
sorted by their confidence scores. The list are then be evaluated with Average
Precision across all of the predictions and all the videos. Formally,

$AP=\sum_{i=1}^{N}p(i)\Delta r(i)$ (7)

where $N=20\times\text{ number if videos}$, $p(i)$ is the precision, and
$r(i)$ is the recall given the first $i$ predictions. For detailed definition
and interpretation of this new metric, readers may refer to Google’s metric
page555https://www.kaggle.com/c/youtube8m##evaluation.

###  4.2 Temporal Segment Data Augmentation

To alleviate missing of temporal information to boost classification
performance while constraining the data size and training time, we propose to
temporally segment videos followed by pooling to generate more temporal
informative features for training video-level models. Specifically, we
temporally split frames in a video in to $N$ segments (each with length
$\lfloor F_{v}/s\rfloor$) and then perform a fixed pooling functions $P$ with
normalization over video segments to generate video-level features:
$\\{P(x_{s_{i-1}:s_{i}}^{v})\\}_{i=1...N}$. In practice we choose $N=3$ and
apply a simple mean pooling function $\mu()$ over frame-level RGB and audio
features to generate additional video-level features for training. With 4
times larger augmented training data $30\rightarrow 120$ GB, we observe
consistent performance boost with the training time for convergence will be
roughly doubled. Additionally, as can be seen in Appendix, increase of $N$ and
other pooling functions such as standard deviation do not help final
performance. [25] proposed to randomly sample one short snippet from each
segment for training, which in our experiment cause over-fitting in the
Youtube-8M dataset.

Similar techniques can also be applied at the inference phase. By feeding the
same model with additional features mean-pooled from 3 temporal segments and
the original mean-pooled features, we merged the 4 inference results with a
fixed weight $(0.1,0.1,0.1,0.7)$ into one final prediction for videos. We name
this approach segmented inference.

###  4.3 Results of Video-level Models

Table 1: Video-level Model Performance, LC stands for latent concept learning, SI stands for segmented inference, DA stands for data augmentation Model Name |  GAP (%) |  mAP (%) |  PERR (%)
---|---|---|---
Baseline1 (MoE2) |  78.41 |  41.58 |  70.9
Baseline3 (MoE8) |  79.30 |  42.20 |  71.9
Deep MoE8 |  79.67 |  43.20 |  72.5
MoRE8 |  81.15 |  44.49 |  73.6
MoRE16 |  81.32 |  44.71 |  73.8
MoRE24 |  81.61 |  47.30 |  74.3
MoHCE12 |  82.15 |  48.76 |  74.6
MoRE24 + LC |  82.37 |  49.63 |  75.1
MoHCE12 + LC |  82.42 |  47.39 |  74.9
MoRE24 + LC + SI |  82.56 |  49.68 |  75.1
MoRE28 + LC + SI |  82.67 |  50.04 |  75.3
MoHCE12 + LC + SI |  82.62 |  47.89 |  75.1
MoRE28 + LC + DA |  82.78 |  50.47 |  75.5
MoRE30 + LC + DA |  82.79 |  50.74 |  75.6
MoRE28 + LC + SI + DA |  82.97 |  50.89 |  75.6
MoHCE12 + LC + SI + DA |  83.28 |  50.41 |  76.0

Table 1 summarizes the performance of video-level models. Models with refined
representation via residual learning (MoRE) and hypercolumn (MoHCE) deliver a
significant $1.5\sim 2\%$ gain in GAP on average. Increasing model capacity by
adding more experts results in consistent however marginal gain. Learning the
latent concept provides an additional $0.5\sim 1\%$ performance boost in GAP
and PERR. However, we notice that learning latent concept may slightly harm
mAP for some models. One possible explanation is the propagation of
classification errors from MoRE.

Temporal segment data augmentation has proved to be useful. Data augmentation
for training fuels up roughly $0.3\%$ GAP gain. Segmented inference by
predicting upon multiple temporal segments brings about $0.2\%$ GAP gain.

###  4.4 Results of Frame-level Models

As can be seen in Table 2, surprisingly, the recurrent neural network models
(with or without recently developed regularization mechanisms) do not achieve
good performance. Attentive pooling methods like deep bag-of-feature, VLAD
seems to be more feasible for frame-level representation learning. We suspect
this phenomenon is rooted from the pre-processed features in the Youtube-8M
dataset. Subtle temporal difference may be missing due to PCA. Deep NetVLAD
model with smaller codebook would achieve better performance. The decreasing
performance with larger codebook size also implies that there might be only a
few distinguishable clusters with the PCA-ed features.

Comparing performance of video and frame-level models, the mean-pooled
features along with simple data augmentation provide a strong model and
achieve competitive single model performance. The result shows that video-
level models may be more cost effective than frame-level models.

Table 2: Frame-level Model Performance Model Name |  GAP(%) |  mAP(%) |  PERR(%)
---|---|---|---
baseline1 LSTM |  79.43 |  37.43 |  72.2
baseline2 DBoW |  78.34 |  36.97 |  70.9
bi-LSTM |  79.98 |  38.51 |  72.5
LSTM + zoneout [12] |  80.15 |  40.27 |  72.4
LSTM + batch_norm [5] |  80.03 |  40.21 |  72.4
LSTM + attn pooling |  80.10 |  40.33 |  73.2
DBoW + attn pooling |  79.89 |  40.94 |  72.5
VLAD + MoRE8 |  81.71 |  46.18 |  74.1
NetVLAD + LC ($k=1$) |  82.47 |  48.37 |  75.1
NetVLAD + LC ($k=2$) |  82.11 |  48.33 |  74.9
NetVLAD + LC ($k=4$) |  81.65 |  47.61 |  74.4
NetVLAD + LC ($k=8$) |  81.29 |  45.54 |  74.1

###  4.5 Ensemble Results

We iteratively grow the ensemble set by randomly adding grouped models from 64
trained models (42 video-level and 22 frame-level models with GAP $>81.0\%$)
followed by a leave-one-out approach to determine fusion weight and remove
detrimental models. Table 3 records the growth of set. Both video and frame-
level model are required to achieve higher GAP.

Table 4 summarizes the final ensemble set with detailed model weights
determined by leave-one-out. Surprisingly, the best single model is not the
most important one for ensemble. The individual model weights are not
proportional to the original model performances. After reviewing the errors
made by different models, we hypothesize that the GAP are greatly affected by
hard examples instead of easy examples. The best model may have a wide range
of coverage but are not specialized in detecting hard examples. Ensemble of
simple models can reach those low-hanging fruits, however, it does not help
classifying hard examples. We also found out that regularization mechanisms
(e.g. dropout) share the similar phenomena. Training complementary models (for
example, ensemble with both video and frame level models) are critical for
achieving better GAP.

We provide our insight why video-level models (with mean-pooled features) and
frame-level models can be complimentary. Mean pooling is robust against noises
over shots in a video. However, in some cases, a few shots define the labels
of a video and therefore frame-level models have a better chance to localize
them. Depending on the underlying shot variants of the videos in the
Youtube-8M dataset, there might be some trade-off between the two approaches.

Finding a approach to guild training complementary models for ensemble remains
open. Readers should be careful when interpreting single model performance and
its significance for ensemble.

Our final ensemble achieve 84.68% GAP on our validation split and 84.66% on
the official test split. Leave-one-out analysis indicates that video level
models contributes 54% of ensemble weight while 46% from frame-level models.
The result shows that to some extent the simple mean-pooled features (video-
level models) are strong enough and are crucial for ensemble. Although frame-
level models are more general and critical for ensemble, they alone might not
be the best choice in sense of performance (roughly the same or slightly
worse) and cost ($\sim 16:1$) in comparison to video-level models.

Table 3: Model Ensemble Performance. V means number of video-level models and F for frame-level # Model | V | F | GAP(%) | mAP(%) | PERR(%)
---|---|---|---|---|---
10 | 10 | 0 | 83.84 | 51.65 | 76.2
12 | 10 | 2 | 84.07 | 52.00 | 76.5
14 | 12 | 2 | 84.34 | 52.86 | 76.8
18 | 12 | 6 | 84.42 | 53.28 | 76.9
18 | 11 | 7 | 84.47 | 53.39 | 76.9
21 | 12 | 9 | 84.52 | 53.64 | 77.0
26 | 17 | 9 | 84.56 | 53.79 | 77.0
34 | 18 | 16 | 84.68 | 54.01 | 77.2

##  5 Conclusion

In this paper, we share our exploration of feasible neural network
architectures for large-scale video tagging. We introduce two enhanced
versions of mixture-of-expert model for video representation learning:
mixture-of-residual-expert (MoRE) model and mixture-of-hypercolumn-expert
(MoHCE) model that boost performance for both video and frame-level models.
The proposed delayed-start layer for latent concept learning also demonstrates
its favored capability to capture the underlying relationships over concepts
and to improve multi-label video classification. The temporal segment data
augmentation provides a simple yet effective way to improve video-level
models. For frame-level models, we observed that pooling-with-attention
methods consistently outperform recurrent neural network models. Ensemble
based on leave-one-out enables us to quantify the importance of heterogeneous
models and achieve 84.662% GAP. While training complementary models seems to
be critical, a mechanism to guide complementarity is still under-explored. We
would expect a improved framework to learn and ensemble specialized models for
large-scale video understanding in the future.

## References

  * [1] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [2] S. Abu-El-Haija, N. Kothari, J. Lee, P. Natsev, G. Toderici, B. Varadarajan, and S. Vijayanarasimhan.  Youtube-8m: A large-scale video classification benchmark.  CoRR, abs/1609.08675, 2016.
  * [3] R. Arandjelovic, P. Gronat, A. Torii, T. Pajdla, and J. Sivic.  Netvlad: Cnn architecture for weakly supervised place recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5297–5307, 2016.
  * [4] A. Bansal, X. Chen, B. Russell, A. G. Ramanan, et al.  Pixelnet: Representation of the pixels, by the pixels, and for the pixels.  arXiv preprint arXiv:1702.06506, 2017.
  * [5] T. Cooijmans, N. Ballas, C. Laurent, Ç. Gülçehre, and A. Courville.  Recurrent batch normalization.  arXiv preprint arXiv:1603.09025, 2016.
  * [6] A. Diba, V. Sharma, and L. V. Gool.  Deep Temporal Linear Encoding Networks.  arXiv preprint arXiv:1611.06678, 2016.
  * [7] J. Duchi, E. Hazan, and Y. Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of Machine Learning Research, 12(Jul):2121–2159, 2011.
  * [8] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [9] S. Hershey, S. Chaudhuri, D. P. Ellis, J. F. Gemmeke, A. Jansen, R. C. Moore, M. Plakal, D. Platt, R. A. Saurous, B. Seybold, et al.  Cnn architectures for large-scale audio classification.  In Acoustics, Speech and Signal Processing (ICASSP), 2017 IEEE International Conference on, pages 131–135. IEEE, 2017.
  * [10] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton.  Adaptive mixtures of local experts.  Neural computation, 3(1):79–87, 1991.
  * [11] Y.-G. Jiang, Z. Wu, J. Wang, X. Xue, and S.-F. Chang.  Exploiting feature and class relationships in video categorization with regularized deep neural networks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2017\.
  * [12] D. Krueger, T. Maharaj, J. Kramár, M. Pezeshki, N. Ballas, N. R. Ke, A. Goyal, Y. Bengio, H. Larochelle, A. Courville, et al.  Zoneout: Regularizing rnns by randomly preserving hidden activations.  arXiv preprint arXiv:1606.01305, 2016.
  * [13] Z. Lan, Y. Zhu, and A. G. Hauptmann.  Deep Local Video Feature for Action Recognition.  arXiv preprint arXiv:1701.07368, 2017.
  * [14] Z.-Z. Lan, L. Jiang, S.-I. Yu, S. Rawat, Y. Cai, C. Gao, S. Xu, H. Shen, X. Li, Y. Wang, et al.  Cmu-informedia at trecvid 2013 multimedia event detection.
  * [15] Z. Li, E. Gavves, M. Jain, and C. G. M. Snoek.  VideoLSTM Convolves, Attends and Flows for Action Recognition.  arXiv preprint arXiv:1607.01794, 2016.
  * [16] T.-Y. Lin, P. Dollár, R. Girshick, K. He, B. Hariharan, and S. Belongie.  Feature pyramid networks for object detection.  arXiv preprint arXiv:1612.03144, 2016.
  * [17] N. Natarajan, I. S. Dhillon, P. K. Ravikumar, and A. Tewari.  Learning with noisy labels.  In Advances in neural information processing systems, pages 1196–1204, 2013.
  * [18] Z. Qiu, T. Yao, and T. Mei.  Deep Quantization: Encoding Convolutional Activations with Deep Generative Model.  arXiv preprint arXiv:1611.09502, 2016.
  * [19] K. Simonyan and A. Zisserman.  Two-Stream Convolutional Networks for Action Recognition in Videos.  NIPS, 2014.
  * [20] N. Srivastava, G. E. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research, 15(1):1929–1958, 2014.
  * [21] S. Sukhbaatar, J. Bruna, M. Paluri, L. Bourdev, and R. Fergus.  Training convolutional networks with noisy labels.  arXiv preprint arXiv:1406.2080, 2014.
  * [22] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich.  Going deeper with convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1–9, 2015.
  * [23] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna.  Rethinking the inception architecture for computer vision.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2818–2826, 2016.
  * [24] L. Wang, Y. Xiong, Z. Wang, and Y. Qiao.  Towards good practices for very deep two-stream convnets.  arXiv preprint arXiv:1507.02159, 2015.
  * [25] L. Wang, Y. Xiong, Z. Wang, Y. Qiao, D. Lin, X. Tang, and L. Van Gool.  Temporal segment networks: Towards good practices for deep action recognition.  arXiv preprint arXiv:1608.00859, 2016.

##  Appendix A Ensemble Table

Table 4: Final Ensemble Set Model Name |  GAP (%) |  mAP (%) |  PERR (%) |  Ensemble Weight
---|---|---|---|---
NetVLAD8_180f_MoRE8_2400_deep_ckpt274921 |  81.294 |  45.541 |  74.1 |  36
NetVLAD8_180f_MoRE8_2400_deep_ckpt280260 |  81.194 |  45.564 |  74.0 |  35
NetVLAD8_180f_MoRE8_1600_deep_ckpt265814 |  81.249 |  46.883 |  74.1 |  34
MoRE30_LC_deep_1st |  82.851 |  50.228 |  75.6 |  31
MoRE30_LC_deep_drop_1st |  82.004 |  49.866 |  75.1 |  28
MoRE28_LC_deep_1st |  80.427 |  48.214 |  73.7 |  28
MoRE30_LC_deep_2nd |  82.868 |  50.522 |  75.6 |  28
NetVLAD4_180f_MoRE8_2400_deep |  82.164 |  48.077 |  74.9 |  27
MoRE30_LC_deep_drop_2nd |  82.038 |  50.232 |  75.2 |  27
MoRE30_LC_deep_drop_3rd |  82.786 |  50.744 |  75.6 |  27
MoHCE16_LC_1st |  82.702 |  48.593 |  75.3 |  27
NetVLAD2_180f_MoRE8_2400_deep_dp |  82.066 |  48.380 |  74.8 |  26
MoRE28_LC_deep_2nd |  81.867 |  48.782 |  74.8 |  26
MoHCE16_LC_2nd |  81.881 |  47.614 |  74.6 |  26
NetVLAD8_180f_MoRE8 |  81.827 |  45.308 |  74.3 |  26
MoRE30_LC_deep_3rd |  81.734 |  49.663 |  74.8 |  25
MoRE30_LC_deep_l3_1st |  81.104 |  48.789 |  74.3 |  25
MoRE28_LC_deep_3rd |  82.783 |  50.472 |  75.5 |  25
MoHCE12_concat_deep_LC |  83.281 |  50.412 |  76.0 |  25
MoRE28_LC_deep_4th |  82.669 |  50.575 |  75.5 |  25
NetVLAD8_1200 |  82.005 |  45.612 |  74.4 |  25
NetVLAD4_180f_MoRE8_2400_deep_dp_ckpt304488 |  81.650 |  47.577 |  74.5 |  24
MoRE28_LC_deep_5th |  81.854 |  50.080 |  74.9 |  24
NetVLAD4_180f_MoRE8_2400_deep_dp_ckpt299007 |  81.604 |  47.610 |  74.4 |  24
NetVLAD1_180f_MoRE8_4096_deep_dp |  82.318 |  48.627 |  75.1 |  24
NetVLAD4_180f_MoRE8_2400 |  81.971 |  45.607 |  74.5 |  24
MoRE28_LC_deep_6th |  81.669 |  49.510 |  74.6 |  23
NetVLAD2_180f_MoRE8_2400_deep |  82.102 |  48.281 |  74.8 |  23
MoRE30_LC_deep_l3_2nd |  82.476 |  49.810 |  75.2 |  23
NetVLAD4_180f_MoRE8 |  81.942 |  45.330 |  74.4 |  23
MoRE28_LC_residual |  82.753 |  50.314 |  75.5 |  22
NetVLAD2_180f_MoRE8_2400_deep_dp |  82.108 |  48.327 |  74.9 |  22
NetVLAD1_180f_MoRE8_4096_deep |  82.048 |  48.309 |  74.7 |  22
NetVLAD1_180f_MoRE8_2400_deep_dp |  82.472 |  48.387 |  75.1 |  22
MoRE28_LC_deep_7th |  82.700 |  50.016 |  75.4 |  22