# Back-Propagation Optimization and Multi-Valued Artificial Neural Networks
for Highly Vivid Structural Color Filter Metasurfaces

Arthur Clini de Souza
Université Côte d’Azur, Inria, CNRS, LJAD, 06902 Sophia Antipolis Cedex,
France
Laboratory of Applied and Computational Electromagnetism (LEMAC),
School of Electrical and Computer Engineering (FEEC),
University of Campinas (UNICAMP), Campinas, São Paulo, Brazil
Solnil, 95 Rue de la République, 13002 Marseille, France
&Stéphane Lanteri
Université Côte d’Azur, Inria, CNRS, LJAD,
06902 Sophia Antipolis Cedex, France
&Hugo Enrique Hernandez-Figueroa
Laboratory of Applied and Computational Electromagnetism (LEMAC),
School of Electrical and Computer Engineering (FEEC),
University of Campinas (UNICAMP), Campinas, São Paulo, Brazil
&Marco Abbarchi
Solnil, 95 Rue de la République,
13002 Marseille, France
Université Aix Marseille, CNRS,
Université de Toulon, IM2NP, UMR 7334,
F-13397 Marseille, France
&David Grosso
Solnil, 95 Rue de la République,
13002 Marseille, France
Université Aix Marseille, CNRS,
Université de Toulon, IM2NP, UMR 7334,
F-13397 Marseille, France
&Badre Kerzabi
Solnil, 95 Rue de la République,
13002 Marseille, France
&Mahmoud Elsawy∗
Université Côte d’Azur, Inria, CNRS, LJAD,
06902 Sophia Antipolis Cedex, France
<EMAIL_ADDRESS>

###### Abstract

We introduce a novel technique for designing color filter metasurfaces using a
data-driven approach based on deep learning. Our innovative approach employs
inverse design principles to identify highly efficient designs that outperform
all the configurations in the training dataset, which consists of 810 distinct
geometries solely. By combining Multi-Valued Artificial Neural Networks and
back-propagation optimization, we overcome the limitations of previous
approaches, such as poor performance due to extrapolation and undesired local
minima. Consequently, we successfully create reliable and highly efficient
configurations for metasurface color filters capable of producing
exceptionally vivid colors that go beyond the sRGB gamut. Furthermore, our
deep learning technique can be extended to design various pixellated
metasurface configurations with different functionalities.

##  1 Introduction

Optical color filter are structures or materials designed to discriminate and
manipulate distinct light wavelengths through the selective transmission or
reflection of particular colors while simultaneously absorbing or attenuating
undesired colors [1, 2]. Conventional color filters rely on the manipulation
of chemical composition to achieve the desired optical properties, which can
lead to issues such as absorption losses, thermal effects, and alterations in
chemical characteristics [3]. An alternative approach involves the utilisation
of structural color filters, offering distinct advantages and applications in
diverse fields such as photorealistic color printing, color holography, anti-
counterfeiting devices, and much more [4, 5, 6].

Metasurfaces have emerged as a promising platform for structural color filters
[7, 8], owing to its peculiar capability of controlling all the light
properties at the nanoscale, enabling a plethora of applications [9, 10, 11,
12]. Dielectric metasurfaces play a crucial role in color filter applications,
especially within the visible spectrum range where the plasmonic conterpart
based on metals is less performing owing to intrinsic optical losses. The
limited losses of dielectrics (e.g. Si3N4, GaN, TiO2, ZrO2, HfO2) make them
highly desirable for designing efficient devices with sharp resonance
responses [13, 14, 15, 16]. Resonant dielectric metasurfaces achieve precise
control over the phase of reflected and transmitted light by leveraging
various resonant phenomena (e.g. Mie resonances) [17, 18]. Through meticulous
engineering of the resonators, selective interaction with different
wavelengths is enabled, leading to efficient and vivid color filters. Such
kind of metasurfaces offer exceptional phase control, high-quality factors,
and sharp resonances, resulting in enhanced color purity and spectral
selectivity [19, 20, 21]. Yet, the design of an ideal color filter demands
capability to selectively filter all colors across the optical spectrum. In
other words, at each desired wavelength, it is crucial to eliminate any
background resonances in order to achieve a pure color response characterized
by sharp reflection or transmission amplitudes. Given the fabrication
constraints, finding the appropriate resonator shape to achieve a desired
response, is a challenging task that has garnered significant attention in the
research community. Numerous studies explored this area, employing
sophisticated optimization algorithms including advanced Deep Learning (DL)
approaches to tackle the inherent complexity of the problem [21, 22, 23, 24,
25, 26, 27, 28, 20, 29]. However, relying on classical optimization approaches
requires several costly simulations when optimizing various color targets
simultaneously [30, 31, 32]. A viable solution for the design of vivid
metasurface color filters is one-shot optimization using Artificial Neural
Network (ANN). However, it is not straightforward owing to the presence of
several designs with similar optical response whereas regular ANN has only a
single output [33].

Here we present a novel data-driven methodology for efficiently designing
fabrication constrained color filter metasurfaces. Our approach combines the
ability to find suitable designs of Multi-Valued Artificial Neural Networks
(MVANNs) with the solution refinement of back-propagation optimization.
Thereby overcoming the fundamental limitations of relying on either the
latter, which will lead to an undesirable local minima, or solely on a MVANN
leading to poor performance associated with extrapolation [33, 23, 24]. In our
case, relying solely on 810 simulations and by varying four parameters allows
for the optimization of a continuous spectrum of objectives and the
identification of highly vivid metasurface color filters. The optimized
geometries exhibit a single sharp resonance response representing all the
colors across the visible regime. To the best of our knowledge, our research
outcomes exceed the previous findings documented in the literature,
positioning our color filter as the foremost advancement in terms of vividness
and overall performance [34, 35, 36, 8].

##  2 Geometry and surrogate model

Figure 1 represents the considered metasurface geometry composed of slanted
ridges of titanium dioxide (TiO2) on top of a silicon dioxide (SiO2) substrate
with refractive index $\mathrm{n_{s}=1.45}$. The refractive index of TiO2 is
determined through ellipsometry (available in the supplementary information).
The inclusion of slanted gratings in the metasurface design introduces
additional degrees of freedom, enabling the appearance of sharp resonances in
the reflection spectrum. This is achieved by breaking the symmetry in the
z-direction, resulting in high-quality resonance modes [37, 38]. Four
parameters, namely the metasurface period (P), resonator height (H), base
width (W), and the ratio between the top and base widths of the resonator (S),
are optimized in this study.

Figure 1: Schematic representation of the simulated structure. The light is
injected from top with normal incidence with electric field polarized along
the x-direction. The inset refers to the geometrical parameters associated to
the single unit-cell surrounded with periodic boundary conditions.

The first step to train the ANN is to generate a dataset. One important aspect
that we considered to build it is the fabrication constraint that imposes a
maximum aspect ratio of 2 between the height and width of the resonators. The
dataset consists of 810 simulations generated by uniformly sweeping the four
parameters indicated in Figure 1. The period P ranges from 250 nm and 510 nm.
For the height H, the range spans from 100 nm to 425 nm. In order to ensure
that the ratios H/W and H/(P-W) remain less than 2, the width W varies between
H/2 and (2P-H)/2, and therefore, the height H must obey the inequality H/2 <
(2P-H)/2. Additionally, the parameter S ranges from 0 to 0.8, with intervals
of 0.1. The simulations were performed on an Intel® Xeon® W-2125 Processor
operating at 4.0 GHz and allocating 8 threads, taking a total of 45 minutes
and 48 seconds. The description of the 2D Python MEEP Finite-Difference Time-
Domain (FDTD) simulations are described in the supplementary information.

Subsequently, a feedforward ANN surrogate model was trained to forecast the
reflection spectral response, considering the resonator’s geometry as input.
In this study, we leverage the surrogate model for two distinct purposes.
Firstly, it enables rapid estimation of the reflection spectrum, achieving
significant computational speed-ups compared to full wave FDTD simulations.
Secondly, the surrogate model is crucial in the inverse design process,
playing a critical role in computing optimal solutions. As a result, the
surrogate model’s performance is of paramount importance, as it must deliver
exceptional precision to provide the most accurate approximations possible. By
fulfilling these requirements, the surrogate model accuracy significantly
contributes to the effectiveness and success of this study. Our surrogate
model is a classical Multilayer Perceptron (MLP) model with fully connected
layers. In this scenario, there are several ways to configure this MLP, and we
will compare two different strategies to map the input geometries into a
spectrum.

The first approach uses the geometric parameters as inputs and generates a
multi-dimensional vector as output. Each dimension of this vector corresponds
to a distinct wavelength in the reflection spectrum. The second approach
involves incorporating the wavelength as an input parameter while outputting a
single dimension representing the reflection amplitude at that particular
wavelength (Figure 2). Although both methods initially seems to be identical,
the second is more convenient for this study. By incorporating the wavelength
as an input, each point in the spectrum becomes a unique data point for the
surrogate model. In contrast, the first method contains only 810 data points
available for training the ANN. Conversely, the second method significantly
enhances the availability of data samples, yielding a total of 405,000
individual data points resulting from sampling the spectrum across 500
different wavelengths.

Figure 2 represents the architecture of the surrogate model. In general, we
consider the Gaussian Error Linear Unit (GELU) activation function which is a
gating activation with a continuous derivative. Considering GELU implies the
continuity of the output and its derivative, providing smoother response
compared to the classical Rectified Linear Unit (ReLU) activation [39]. It is
worth mentioning that the last layer contains a single linear neuron to
compute a singular output dimension aligning with the reflection observed at
the designated wavelength as indicated by the red arrow in Figure 2.

Figure 2: The diagram represents the architecture of the surrogate model ANN.
The green arrows represent the input parameters. The surrogate model is
designed to calculate a single output dimension, which corresponds to the
reflection at the desired wavelength as shown by the red arrow. On the left
side, a typical example of the resonator’s geometry is provided, while on the
right side, an example spectrum is displayed. It is important to note that
these examples are purely for illustrative purposes and do not represent any
real simulations. All the layers in this architecture are initialized using
Glorot normal initializers as discussed in Ref. [40]. The model was trained
using the TensorFlow framework, employing a holdout methodology [41]. The
dataset is divided into two subsets: 769 simulations to update the weights of
the ANN and 41 simulations for validation. The optimizer of choice was Adam
with MSE loss and batch size of 512 [42]. The training stopped at epoch 252
due to the lack of improvement of the validation loss within the last 100
epochs. The final MSE loss was $2.1443\times 10^{-4}$ and validation loss,
$4.3698\times 10^{-4}$. The total training time was 18 min 6s on a Google
colab’s NVIDIA® T4 GPU.

Figure 3 illustrates the generalization capacity of the surrogate model.
Despite being trained on a relatively small dataset, the model is able to
reproduce the spectra associated to the geometries that lie beyond the
training dataset. This accomplishment can be attributed to the integration of
wavelength as an input variable and the meticulous training methodology
employed to mitigate overfitting. The demonstrated efficiency of this approach
underscores the potential for constructing dependable surrogate models even
when confronted with restricted training data.

Figure 3: Comparison between fullwave simulation and the prediction from the
surrogate model. Top figures indicate the reflection spectra of three randomly
selected samples from the dataset. The bottom figures corresponding to the
same top geometries while introducing a small random deviation to the
geometrical parameters. Notably, the surrogate model demonstrates robust
generalization capabilities despite the inherent sparsity of the dataset.

##  3 Inverse Design

In this particular section, an inverse design methodology is formulated for
metasurface color filter. The primary objective is to determine the
corresponding set of parameters based on the desired spectrum response. To
initiate the process, the target line shape spectrum must be established, and
in our research, a Lorentzian function shape resembling a Fano resonance-like
shape is employed:

$\centering
L\left(f\right)=\frac{\omega^{2}}{\omega^{2}+4\left(f-f_{0}\right)^{2}},\@add@centering$
(1)

where $\omega$ represents the full width at half maximum (FWHM) while $f_{0}$
corresponds to the central frequency. The units employed are in accordance
with MEEP’s configuration, where both frequency and FWHM are measured in
$\mathrm{\mu m^{-1}}$. Additionally, it is worth noting that throughout the
optimization process, the FWHM remains unchanged for all targets. This means
that the Quality Factor (QF) of the desired spectra undergoes variations as a
function of $f_{0}$, in accordance with Equation 2.

$QF:=\frac{f_{0}}{\omega}$ (2)

Figure 4: Optimization diagram for the inverse design methodology, where the
target reflection spectrum is specified in the cyan curve. To achieve the
desired outcome, two complementary paths are identified. The first path
involves exploiting the dataset scan and searching for the best design that
will be subjected to a back-propagation optimization to fine-tune the
parameters, as indicated in the top blocks. On the other hand, the second path
involves passing the target through the MVANN to obtain 20 designs. These
designs are then validated using the surrogate model, and the solution that
presented the least MSE compared to the target is chosen. Finally, a back-
propagation optimization is performed on that design. Ultimately, the best
solution among the four selected designs from the two paths is chosen. Further
details can be found in the text.

The optimization scenario illustrated in Figure 4 depicts the inverse design
process. The desired spectrum is represented by the cyan curve in the top left
corner. In this stage, we employ two complementary approaches for inverse
design. The first approach involves a straightforward search through a
dataset, represented by the arrow pointing to the right of the target
spectrum. For each simulation in the dataset, we calculate the MSE loss
between the target spectrum and the simulated spectrum. By doing so, we can
identify the design that produces the lowest MSE. Although this design is the
best within the dataset, it fails to accurately match the desired response, as
indicated in Figure 5. To address this limitation, we utilize a fully
differentiable surrogate model and perform gradient optimization to finely
adjust the parameters. As it can be seen in Figure 5 running the back-
propgation yields interesting designs. However, still this approach is not
able to produce vivid colors along all the visible regime. Further the
resonances are mainly associated with higher order modes appearing at shorter
wavelengths (see last column) that allows for non-pure color. The optimized
solution, along with the best design found within the dataset, are both stored
for later comparison.

Figure 5: Results from the dataset search using a Lorentzian lineshape in
frequency as target for the dataset search for the wavelengths of
$420\text{\,}\mathrm{n}\mathrm{m}$ $480\text{\,}\mathrm{n}\mathrm{m}$,
$500\text{\,}\mathrm{n}\mathrm{m}$, $530\text{\,}\mathrm{n}\mathrm{m}$,
$650\text{\,}\mathrm{n}\mathrm{m}$ $650\text{\,}\mathrm{n}\mathrm{m}$ with
$\omega=0.01$ from left to right, top to bottom. The metric used to evaluate
the results was the MSE (L) between the target spectrum and the closest one in
the dataset, while the geometric parameters (P, H, W and S) were described in
Figure 1. Blue and red curves represent the best design without and with back-
propagation optimization on best design found in the dataset. A comparison of
sharp resonances and broader resonances can be found in the supplementary
information section.

The second complementary approach, indicated by the downward arrow from the
target spectrum in Figure 4, involves the utilization of a MVANN to generate
robust solutions for the inverse design problem. A MVANN is essentially an ANN
capable of producing multiple solutions based on a single input. In the
context of our study, the adoption of the MVANN proves highly advantageous as
it effectively addresses the challenge of multiple parameter responses
corresponding to a given objective target [33]. To ensure reasonable accuracy,
we choose to employ a model that yields 20 outputs. The detailed architecture
and training procedure of the MVANN are provided in the supplementary
information section. After acquiring the set of 20 distinct designs, a
validation process is conducted using the surrogate model. Among these
designs, the one that exhibits the lowest MSE loss when compared to the target
is initially stored for future comparison. To further improve the results,
this best design is subjected to back-propagation optimization. By employing a
combination of the dataset search method indicated by the top blocks in Figure
4 and the MVANN approach, followed by back-propagation indicated by the bottom
blocks in Figure 4, we successfully identify four robust designs. This entire
process, encompassing both approaches, takes approximately 100 seconds to
complete. Finally, the design exhibiting the lowest MSE loss compared to the
target is selected as the final choice.

By employing the aforementioned inverse design methodology illustrated in
Figure 4, we perform a total of 100 optimization iterations, as depicted in
Figure 6 (a). These iterations involve sweeping the desired wavelength across
the range of $400\text{\,}\mathrm{n}\mathrm{m}$ to
$700\text{\,}\mathrm{n}\mathrm{m}$, thereby encompassing the entire spectrum
of colors. The target function is defined as a Lorentzian spectrum with a Full
Width at Half Maximum (FWHM) value of $\omega=0.01$. Subsequently, we
calculate the perceived color values (x and y) on the chromaticity diagram
using the CIE 1931 $\mathrm{2^{\degree}}$ standard observer color matching
functions and the "Equal energy" illuminant, which ensures that incident light
of all wavelengths possesses uniform power [43]. As depicted Figure 6 (a), we
successful optimized metasurface designs that generate colors closely
resembling our target (gray curve). Consequently, a significant number of
these optimized designs fall outside the sRGB zone, represented by the dashed
triangle in Figure 6 (a), expanding the color gamut achievable with
metasurfaces beyond the limitations of conventional display technologies.

In Figure 6 (b), we select six distinct designs to demonstrate the spectral
response of our optimized configurations, as indicated by the black arrows in
Figure 6 (a). Our inverse design approach has led to metasurface designs that
exhibit a reflection spectrum response that is consistent with the target
Lorentzian spectrum across all desired wavelengths. Moreover, the background
response in the spectrum is nearly flat, providing a vibrant color response
for our optimized configurations. A field profile for one of the optimized
geometries is presented in the supplementary information section.
Additionally, a numerical comparison between FDTD calculations and COMSOL has
been included to verify our results. Furthermore, we have conducted a full-
wave simulation with finite structure to identify the number of required
periods to retrive the results of a single unitcell simulation with periodic
boundary conditions (for more details, please refer to the supplementary
section). To the best of our knowledge, the optimized metasurface designs
yield the most vivid color filter reported so far in the literature [34, 35,
36, 8]. It is worth mentioning that the entire optimization process took 2
hours, 47 minutes, and 41 seconds. Among this time, 1 hour, 29 minutes, and 57
seconds were dedicated to simulations, 1 hour, 25 minutes, and 9 seconds were
allocated for back-propagation optimization, 8 minutes and 30 seconds were
spent on predictions using both the surrogate model and the MVANN, and 330
milliseconds were utilized for searching the best design within the dataset.

Figure 6: (a) Chromaticity diagram depicting the outcomes of the 100
optimization results. The arrows in the diagram indicate the corresponding
designs showcased in (b). The gray line represents the calculated perceived
color corresponding to the target spectrum. (b) detailed examination of 6
different optimizations, displaying their dimensions (in nm) of each
geometrical parameter and the corresponding sRGB color representation of the
simulated spectrum

The results presented in Figure 6 show that the optimized designs have varying
heights. However, considering the complexity associated with fabricating
structures of different heights, it is preferable to find designs with a fixed
height. Therefore, we introduce a constraint to the back-propagation
optimization process, requiring all designs to have equal heights. To achieve
this, we modified the inverse design workflow as follows: Initially, we
obtained the lowest 5 losses from the MVANN for each of the 6 different target
designs, each representing a different color and having a different height.
Subsequently, we identified designs that were closer in their height values.
We calculated the mean height of these selected designs for each color and
used this mean value as an initial guess for optimizing each target using the
back-propagation optimization process.

To accommodate the height constraint, we introduced a dedicated layer in the
surrogate model. This layer connected the height input solely to the initial
guess and consisted of a single weight and bias. This shared layer was
employed throughout all optimization iterations. By optimizing all 6 models
simultaneously within the same batch, we ensured that the loss for each
individual model progressively decreased with each training epoch.
Consequently, we transformed the fixed height constraint into an optimized
parameter, allowing us to determine the most suitable height value for
achieving the desired outcomes across all six colors.

The results presented in Figure 7 demonstrate the effectiveness of our
approach in achieving favorable outcomes for most of the considered colors,
with the exception of the last color. It is noteworthy that the optimization
of this particular vivid red color posed additional challenges, even in the
absence of height constraints. These challenges indicate inherent complexities
in the optimization process for achieving this specific color and necessitate
the utilization of specialized design approaches [20, 44], which are beyond
the scope of our current work. It is important to mention that our study
focuses on achieving pure colors across the entire spectrum, rather than
emphasizing a single color. Furthermore, it is worth noting that prior studies
employed half the number of simulations to optimize a single color compared to
our methodology.

Figure 7: Optimization results based on a fixed height configurations. (a)
refers to the chromaticity diagram of 6 optimized designs. The simulated
spectra is given in (b).

##  4 Conclusion

In conclusion, our research presents a novel DL methodology for optimizing
resilient designs of vivid color filter metasurfaces. By utilizing a surrogate
model constructed from a dataset of only 810 simulations, our approach
demonstrates exceptional efficiency in the optimization process. The numerical
tool developed in this study enables the cost-effective fabrication of
structural color filters by exploring a wide range of narrow line shapes that
exhibit high-quality resonances, aligning with desired spectral reflection
responses. Notably, our methodology expands the color gamut beyond the
conventional RGB colors, offering unprecedented versatility in color
generation. Furthermore, our DL approach successfully respects fabrication
constraints, ensuring practical feasibility. The achievements of our research
significantly contribute to the field of optical device design. By pushing the
boundaries of metasurface optimization, we open up new possibilities for the
development of advanced optical devices. The proposed methodology holds
promise for various applications, such as display technologies, data encoding,
and artistic expression. These notable advancements not only enhance the
understanding of metasurface design principles but also provide valuable
insights for future research endeavors.

## Author contributions statement

A.C.D developed the DL architecture and performed the numerical simulation. ME
proposed the structure and coordinated the work. M.E, A.C.D, S.L discussed the
numerical results. M.E, A.C.D, H.E.H.F, and M.A discussed the physical
results. All authors contributed to the writing and reviewed the manuscript.
The authors declare that they have no competing interests.

## Acknowledgment

Hugo EHF would like to thank the support of Brazilian government agencies
CAPES and CNPq (through project: 312714/2019-0), and the São Paulo Research
Foundation (FAPESP) agency through the projects: 2021/11380-5 (CCD/CPTEn) and
2021/06506-0.

## References

  * [1] Craig F Bohren and Donald R Huffman.  Absorption and scattering of light by small particles.  John Wiley & Sons, 2008.
  * [2] Milton Kerker.  The scattering of light and other electromagnetic radiation.  Elsevier, 2016.
  * [3] Mark A Brown and Stephen C De Vito.  Predicting azo dye toxicity.  Critical reviews in environmental science and technology, 23(3):249–324, 1993.
  * [4] Hyoki Kim, Jianping Ge, Junhoi Kim, Sung-eun Choi, Hosuk Lee, Howon Lee, Wook Park, Yadong Yin, and Sunghoon Kwon.  Structural colour printing using a magnetically tunable and lithographically fixable photonic crystal.  Nature Photonics, 3(9):534–540, 2009.
  * [5] Wei Ting Chen, Kuang-Yu Yang, Chih-Ming Wang, Yao-Wei Huang, Greg Sun, I-Da Chiang, Chun Yen Liao, Wei-Lun Hsu, Hao Tsun Lin, Shulin Sun, et al.  High-efficiency broadband meta-hologram with polarization-controlled dual images.  Nano letters, 14(1):225–230, 2014.
  * [6] Haibo Hu, Qian-Wang Chen, Jian Tang, Xian-Yi Hu, and Xu-Hui Zhou.  Photonic anti-counterfeiting using structural colors derived from magnetic-responsive photonic crystals with double photonic bandgap heterostructures.  Journal of Materials Chemistry, 22(22):11048–11053, 2012.
  * [7] Fei Ding, Anders Pors, and Sergey I Bozhevolnyi.  Gradient metasurfaces: a review of fundamentals and applications.  Reports on Progress in Physics, 81(2):026401, 2017.
  * [8] Egor Khaidarov, Damien Eschimese, Keng Heng Lai, Aihong Huang, Yuan Hsing Fu, Qunying Lin, Ramon Paniagua-Dominguez, and Arseniy I Kuznetsov.  Large-scale vivid metasurface color printing using advanced 12-in. immersion photolithography.  Scientific Reports, 12(1):14044, 2022.
  * [9] Arseniy I Kuznetsov, Andrey E Miroshnichenko, Mark L Brongersma, Yuri S Kivshar, and Boris Luk’yanchuk.  Optically resonant dielectric nanostructures.  Science, 354(6314):aag2472, 2016.
  * [10] Guoxing Zheng, Holger Mühlenbernd, Mitchell Kenney, Guixin Li, Thomas Zentgraf, and Shuang Zhang.  Metasurface holograms reaching 80% efficiency.  Nature nanotechnology, 10(4):308–312, 2015.
  * [11] Ehsan Arbabi, Seyedeh Mahsa Kamali, Amir Arbabi, and Andrei Faraon.  Full-stokes imaging polarimetry using dielectric metasurfaces.  Acs Photonics, 5(8):3132–3140, 2018.
  * [12] Andreas Tittl, Aleksandrs Leitis, Mingkai Liu, Filiz Yesilkoy, Duk-Yong Choi, Dragomir N Neshev, Yuri S Kivshar, and Hatice Altug.  Imaging-based molecular barcoding with pixelated dielectric metasurfaces.  Science, 360(6393):1105–1109, 2018.
  * [13] Yuanmu Yang, Ivan I Kravchenko, Dayrl P Briggs, and Jason Valentine.  All-dielectric metasurface analogue of electromagnetically induced transparency.  Nature communications, 5(1):5753, 2014.
  * [14] Amir Arbabi, Yu Horie, Mahmood Bagheri, and Andrei Faraon.  Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission.  Nature nanotechnology, 10(11):937–943, 2015.
  * [15] Mikhail I Shalaev, Jingbo Sun, Alexander Tsukernik, Apra Pandey, Kirill Nikolskiy, and Natalia M Litchinitser.  High-efficiency all-dielectric metasurfaces for ultracompact beam manipulation in transmission mode.  Nano letters, 15(9):6261–6266, 2015.
  * [16] Dianmin Lin, Pengyu Fan, Erez Hasman, and Mark L Brongersma.  Dielectric gradient metasurface optical elements.  science, 345(6194):298–302, 2014.
  * [17] Manuel Decker, Isabelle Staude, Matthias Falkner, Jason Dominguez, Dragomir N Neshev, Igal Brener, Thomas Pertsch, and Yuri S Kivshar.  High-efficiency dielectric huygens’ surfaces.  Advanced Optical Materials, 3(6):813–820, 2015.
  * [18] Ye Feng Yu, Alexander Y Zhu, Ramón Paniagua-Domínguez, Yuan Hsing Fu, Boris Luk’yanchuk, and Arseniy I Kuznetsov.  High-transmission dielectric metasurface with 2$\pi$ phase control at visible wavelengths.  Laser & Photonics Reviews, 9(4):412–418, 2015.
  * [19] Han-Don Um, Deokjae Choi, Amit Solanki, Emerald Huang, Kwanyong Seo, and Fawwaz Habbal.  Dynamic selection of visible wavelengths using resonant tio2 nanostructures.  Nanophotonics, (0), 2023.
  * [20] Ronghui Lin, Vytautas Valuckas, Thi Thu Ha Do, Arash Nemati, Arseniy I Kuznetsov, Jinghua Teng, and Son Tung Ha.  An efficient neural optimizer for resonant nanostructures: demonstration of highly-saturated red silicon structural color.  arXiv preprint arXiv:2304.13516, 2023.
  * [21] Xiao Shang, Jiebin Niu, He Li, Longjie Li, Huakui Hu, Cheng Lu, and Lina Shi.  Polarization-sensitive structural colors based on anisotropic silicon metasurfaces.  In Photonics, volume 10, page 448. MDPI, 2023.
  * [22] Wei Ma, Feng Cheng, Yihao Xu, Qinlong Wen, and Yongmin Liu.  Probabilistic representation and inverse design of metamaterials based on a deep generative model with semi-supervised learning strategy.  Advanced Materials, 31(35):1901111, 2019.
  * [23] John Peurifoy, Yichen Shen, Li Jing, Yi Yang, Fidel Cano-Renteria, Brendan G DeLacy, John D Joannopoulos, Max Tegmark, and Marin Soljačić.  Nanophotonic particle simulation and inverse design using artificial neural networks.  Science advances, 4(6):eaar4206, 2018.
  * [24] Takashi Asano and Susumu Noda.  Optimization of photonic crystal nanocavities based on deep learning.  Optics express, 26(25):32704–32717, 2018.
  * [25] Jiaqi Jiang and Jonathan A Fan.  Global optimization of dielectric metasurfaces using a physics-driven neural network.  Nano letters, 19(8):5366–5372, 2019.
  * [26] Fufang Wen, Jiaqi Jiang, and Jonathan A Fan.  Progressive-growing of generative adversarial networks for metasurface optimization.  arXiv preprint arXiv:1911.13029, 2019.
  * [27] Alireza Makhzani, Jonathon Shlens, Navdeep Jaitly, Ian Goodfellow, and Brendan Frey.  Adversarial autoencoders.  arXiv preprint arXiv:1511.05644, 2015.
  * [28] Fufang Wen, Jiaqi Jiang, and Jonathan A Fan.  Robust freeform metasurface design based on progressively growing generative networks.  ACS Photonics, 7(8):2098–2104, 2020.
  * [29] Peter R Wiecha and Otto L Muskens.  Deep learning meets nanophotonics: a generalized accurate predictor for near fields and far fields of arbitrary 3d nanostructures.  Nano letters, 20(1):329–338, 2019.
  * [30] Mahmoud MR Elsawy, Stéphane Lanteri, Régis Duvigneau, Gauthier Brière, Mohamed Sabry Mohamed, and Patrice Genevet.  Global optimization of metasurface designs using statistical learning methods.  Scientific reports, 9(1):1–15, 2019.
  * [31] Dong Xu, Yu Luo, Jun Luo, Mingbo Pu, Yaxin Zhang, Yinli Ha, and Xiangang Luo.  Efficient design of a dielectric metasurface with transfer learning and genetic algorithm.  Optical Materials Express, 11(7):1852–1862, 2021.
  * [32] Samad Jafar-Zanjani, Sandeep Inampudi, and Hossein Mosallaei.  Adaptive genetic algorithm for optical metasurfaces design.  Scientific reports, 8(1):11040, 2018.
  * [33] Chao Zhang, Jing Jin, Weicong Na, Qi-Jun Zhang, and Ming Yu.  Multivalued neural network inverse modeling and applications to microwave filters.  IEEE Transactions on Microwave Theory and Techniques, 66(8):3781–3797, 2018.
  * [34] Hongjun Liu, Huiren Peng, Ke Li, Liping Lu, Junhong Deng, Yanjun Liu, Chengfeng Qiu, Guixin Li, and Xing Cheng.  Transfer printing of solution-processed 3d zno nanostructures with ultra-high yield for flexible metasurface color filter.  Advanced Materials Interfaces, 9(21):2101963, 2022.
  * [35] Soumyashree S Panda and Ravi S Hegde.  Transmission-mode all-dielectric metasurface color filter arrays designed by evolutionary search.  Journal of Nanophotonics, 14(1):016014–016014, 2020.
  * [36] Yicheng Wang, Weikai Huang, Yu-Sheng Lin, and Bo-Ru Yang.  A tunable color filter using a hybrid metasurface composed of zno nanopillars and ag nanoholes.  Nanoscale Advances, 4(17):3624–3633, 2022.
  * [37] D Barton III, M Lawrence, and J Dionne.  Wavefront shaping and modulation with resonant electro-optic phase gradient metasurfaces.  Applied Physics Letters, 118(7):071104, 2021.
  * [38] Mark Lawrence, David R Barton III, Jefferson Dixon, Jung-Hwan Song, Jorik van de Groep, Mark L Brongersma, and Jennifer A Dionne.  High quality factor phase gradient metasurfaces.  Nature Nanotechnology, 15(11):956–961, 2020.
  * [39] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus).  arXiv preprint arXiv:1606.08415, 2016.
  * [40] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings, 2010.
  * [41] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  Tensorflow: a system for large-scale machine learning.  In Osdi, volume 16, pages 265–283. Savannah, GA, USA, 2016.
  * [42] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [43] János Schanda.  Colorimetry: understanding the CIE system.  John Wiley & Sons, 2007.
  * [44] Zhaogang Dong, Lei Jin, Soroosh Daqiqeh Rezaei, Hao Wang, Yang Chen, Febiana Tjiptoharsono, Jinfa Ho, Sergey Gorelik, Ray Jia Hong Ng, Qifeng Ruan, et al.  Schrödinger’s red pixel by quasi-bound-states-in-the-continuum.  Science Advances, 8(8):eabm4512, 2022.


submitted February 2022

# Coronal Mass Ejection image edge detection in Heliospheric Imager STEREO
SECCHI data

Marc D. Nichitiu  The Stony Brook School, NY, USA

###### Abstract

We present an algorithm to detect the outer edges of Coronal Mass Ejection
(CME) events as seen in differences of Heliospheric Imager STEREO SECCHI HI-1
images from either A or B spacecraft, as well as its implementation in Python.

##  1 Introduction

Studying and monitoring solar activity, and in particular, its stronger
radiation and particle emission, is important for both science and practical
purpose, for our life on and around Earth. In particular, Coronal Mass
Ejections (CME), which carry magnetized plasma, can induce geomagnetic storms
on our planet, and also endanger astronauts or space equipment. Increased
interest in studying their origins, structure, and trajectories resulted in
more and more international space missions, among which several ones are still
operational. The STEREO pair A and B of NASA spacecraft, out of which STEREO-B
was lost in 2014, is one of such latest missions, continuously imaging the Sun
and the space between Sun and Earth, with several instruments, as described in
[Eyles et al.(2009)Eyles, Harrison, Davis, Waltham, Shaughnessy, Mapson-
Menard, Bewsher, Crothers, Davies, Simnett, Howard, Moses, Newmark, Socker,
Halain, Defise, Mazy, & Rochus, Eyles et al.(2007)Eyles, Davis, Harrison,
Waltham, Halain, Mazy, Defise, Howard, Moses, Newmark, & Plunkett] by orbiting
around the Sun at roughly the same distance as the Earth. Early years of raw
and processed image data are available from both, and STEREO-A is still
operational. Two visible-light cameras HI-1 and HI-2 on each spacecraft are
thus precisely aimed in order to cover the line Sun-Earth, while being
shielded from the Sun itself. Their images show various parts of the solar
corona, receiving free-electron scattered light from the K-corona, dust
scattered light from the F-corona, as well as light from background stars.
Many studies of these images start from running differences, as described by
[Sheeley et al.(1997)Sheeley, Wang, Hawley, Brueckner, Dere, Howard, Koomen,
Korendyke, Michels, Paswaters, Socker, St. Cyr, Wang, Lamy, Llebaria, Schwenn,
Simnett, Plunkett, & Biesecker] for coronagraph images, to enhance the faint
K-corona light and remove the slowly-varying but much stronger F-corona light,
as also explained in [Davies et al.(2009)Davies, Harrison, Rouillard, Sheeley,
Perry, Bewsher, Davis, Eyles, Crothers, & Brown]. CME perturb the K-corona,
and thus appear as bright regions in these images. Successive images, obtained
from space missions such as the SOHO LASCO [Yashiro et al.(2004)Yashiro,
Gopalswamy, Michalek, St. Cyr, Plunkett, Rich, & Howard], and STEREO
coronagraph and HI cameras, show their propagation and evolution. Considerable
research studying the origin of CMEs, shape, and trajectories is in progress,
but open problems remain, as described in the very comprehensive review from
[Webb & Howard(2012)], in particular regarding propelling forces and
interaction with solar wind: such models can be found in [Xie et al.(2006)Xie,
Gopalswamy, Ofman, St. Cyr, Michalek, Lara, & Yashiro, Michalek et
al.(2006)Michalek, Gopalswamy, Lara, & Yashiro, Gopalswamy & Yashiro(2007),
Howard et al.(2008)Howard, Nandy, & Koepke]. Given the large amount of data as
well as the need for tracking, computerized pre-analysis tools can provide
important help, in particular to identify the CME boundaries and kinematics.

Computer algorithms able to identify boundaries of various objects in images
have been already studied. They can use gradient techniques such as in
[Canny(1986)], or active contour techniques: [Caselles et al.(1993)Caselles,
Catté, Coll, & Dibos, Chan & Vese(2001), Márquez-Neila et al.(2014)Márquez-
Neila, Baumela, & Alvarez]. However, as [Young & Gallagher(2008)] explain, it
is not straightforward to use these existing techniques because of the diffuse
nature of CMEs. Moreover, they also have complex inner features, and their
signal can be affected by the presence of bright stars and even the Milky Way
passing in the field of view, which adds to the noise.

Thus, several complex algorithmic methods specialized for this task have been
devised: CACTus [Robbrecht & Berghmans(2004)], SEEDS [Olmedo et
al.(2008)Olmedo, Zhang, Wechsler, Poland, & Borne], and multiscale wavelet
analysis [Young & Gallagher(2008), Byrne et al.(2009)Byrne, Gallagher,
McAteer, & Young], mainly focused on the coronagraph images. These methods
isolate the regions of interest in specific ways, and then estimate geometric
and kinematic parameters, using underlying assumptions such as constant
acceleration and elliptical shapes.

The algorithm presented here only focuses on the image segmentation stage,
identifying the outer edge of CME regions in STEREO HI images using a staged
approach composed of simple steps. The algorithm starts by filtering and
smoothing the image difference, subsequently exploring it from the center of
the bright CME region, marking its outer edges.

Through subsequent work, it can be integrated into various kinematic
estimation schemes. Given its modular nature, the algorithm itself could also
be further adapted to examine specific internal features apparent in the
images.

The rest of the paper is organized as follows. In Section 2 we describe the
algorithm outline and its subsequent processing of an image. In Section 3 we
show its results on several other STEREO HI-1 images of known CME, and in
Section 4 we discuss future lines of work.

##  2 Algorithm

Figure 1: 2008 CME event image (STEREO HI-1 Spacecraft B) used as input. The
traveling CME appears as a bright region entering the field of vision on the
left. The image is obtained by normalizing the difference of two consecutive
background-filtered images to a 0-255 range of pixel values.

In Figure 1 we see an example of what the algorithm starts with, which is an
early phase of a CME visible in the differenced image from December 12, 2008.
The Sun is behind the right edge of the image, which has $n_{c}\times
n_{r}=1024\times 1024$ pixels.

The main steps, graphically illustrated in Figure 2, are as follows, where
$I[c,r]$ represents the pixel intensity for coordinates column $c$ and row $r$
of the input image. The origin is the bottom left of the image.

  1. 1.

Find the center of mass $(x_{S},y_{S})$ of the bright area and set
$C=(c_{S},y_{S})$ where $c_{S}$ is close to the image edge $c_{0}$ towards the
Sun (right for spacecraft A, so typically $c_{0}=1023$, and left for B, i.e.
$c_{0}=0$).

  2. 2.

Going from $c_{S}$ to $c_{E}$ (1023 for B, 0 for A), on each image row $r$,
mark the farthest bright pixels seen along. Then do the same on each image
column $c$, first from $y_{S}$ up, and then from $y_{S}$ down, with
$t_{1}=240$ for 0-255 gray level images.

1:$e_{1}\leftarrow[]$

2:for each image row $r$ do

3: $q\leftarrow$False

4: $v\leftarrow c_{0}$

5: for each image column $c$ from $c_{S}$ to $c_{E}$ do

6: if $I[c,r]>t_{1}$ then

7: $v\leftarrow c$

8: $q\leftarrow$True

9: end if

10: end for

11: if $q$ then

12: $e_{1}+\leftarrow[v,r]$

13: end if

14:end for

15:for each image column $c$ do

16: for each image row $r$ from $y_{S}$ to $r_{E}$ do

17: if $I[c,r]>t_{1}$ then

18: $v\leftarrow r$

19: $q\leftarrow$True

20: end if

21: end for

22: if $q$ then

23: $e_{1}+\leftarrow[c,v]$

24: end if

25: for each image row $r$ from $y_{S}$ to $0$ do

26: if $I[c,r]>t_{1}$ then

27: $v\leftarrow r$

28: $q\leftarrow$True

29: end if

30: end for

31: if $q$ then

32: $e_{1}+\leftarrow[c,v]$

33: end if

34:end for

  3. 3.

Smooth the image $I$ into $I_{M}$, using a sum-scale-and-modulo-256 5x5 filter
to enhance contrast at the edges

  4. 4.

Mark the contrasting pixels created in step 3, with $t_{2}=120$ for 0-255 gray
level images

1:$e_{2}\leftarrow[]$

2:for each image row $r$ do

3: for each image column $c$ except the last one do

4: if $|I_{M}[c,r]-I_{M}[c+1,r]|>t_{2}$ then

5: $e_{2}+\leftarrow[c,r]$

6: end if

7: end for

8:end for

9:for each image column $c$ do

10: for each image row $r$ except the last one do

11: if $|I_{M}[c,r]-I_{M}[c,r+1]|>t_{2}$ then

12: $e_{2}+\leftarrow[c,r]$

13: end if

14: end for

15:end for

  5. 5.

Mark the centers where both kinds of previous marks are present in a
$w_{i}\times w_{i}$ window, as created in step 2 and step 4, namely from the
$e_{1}$ and $e_{2}$ lists:

1:$I_{1}=\textrm{ones}[n_{c},n_{r}]$

2:$I_{2}=\textrm{ones}[n_{c},n_{r}]$

3:for each element $[x,y]$ of $e_{1}$ do

4: $I_{1}[x,y]\leftarrow 2$

5:end for

6:for each element $[x,y]$ of $e_{2}$ do

7: $I_{2}[x,y]\leftarrow 3$

8:end for

9:$I_{3}\leftarrow I_{1}+I_{2}$

10:$I_{4}\leftarrow\textrm{filter}(I_{3},\textrm{prod},w_{i})$

11:$e_{3}\leftarrow$ list of $[x,y]$ such that $I_{4}[x,y]$ is a multiple of
$2$ and $3$, or of $5$

  6. 6.

Estimate the density of the markers from step 5 in a larger window
$w_{d}\times w_{d}$ and mark those above half into a boolean mask $I_{D}$.

  7. 7.

Circularly sweeping from the center, for each such radial half-line gather the
last mark as created in step 6, i.e. from $I_{D}$, and build the perimeter:

1:$e_{4}\leftarrow[]$

2:for each angle $\alpha$ from $0$ to $2\pi$ in steps of $\varepsilon$ do

3: $v\leftarrow$Nothing

4: for each pixel $[x,y]$ along the $\alpha$ ray from $[c_{S},y_{S}]$ to any
edge of the image do

5: if $I_{D}[x,y]$ then

6: $v\leftarrow[x,y]$

7: end if

8: end for

9: if $v$ then

10: $e_{4}+\leftarrow v$

11: end if

12:end for

  8. 8.

Remove spikes (outliers) created in step 7 in the list $e_{4}$ of pairs of
pixel coordinates, reconnect and validate the final perimeter, based on CME
geometry from [Fisher(1984), Fisher & Munro(1984), Crifo(1983)].

Figure 2: Diagram illustrating the different steps of the algorithm

The idea behind step 3 is to augment existing gradual contrasting zones by
creating an abrupt transition through the modulo operation.

Given the larger density of markers around the true edge, even if the rays in
step 7 are one or two pixels wide, they should not miss them when they cross
it outwards.

Outliers appear occasionally because of brighter background light points far
away from the CM, and can be removed in step 8. The validation of the
perimeter, to eliminate unphysical results, is based on CME general geometry,
using bubble and loop models as described in [Crifo(1983)]. Specifically, as
[Fisher(1984)] and [Fisher(1984)] point out, the width to depth ratio is
typically 3:2, so we compute and compare these diameters.

##  3 Results

We have tracked CME events using a set of 100 STEREO 11-day background-removed
L2 HI-1 a and b images from the UK Solar System Data Center,
https://www.ukssdc.ac.uk. The algorithm performed with overall high tracking
capabilities and without any fine-tuning or training set. We also conducted a
manual identification of the CME regions, and then compared it with the
automatic one. We counted the relative differences in areas between the two
perimeters as errors, and we measured by hand that the percentage error of
region identifications was below 15% for the most irregular events. We did not
find false positives. The images in Figure 3 and Figure 5 illustrate these
cases.

Figure 3: A few examples of processing results: the program recognized the CME
bright region boundary and marked it with dark pixels, forming a contour line.
We can notice that inner darker regions do not prevent the program from
finding the outermost edges.

A simplified initial calculation appears in agreement with speeds of the order
of $300$km/s from [Davis et al.(2009)Davis, Davies, Lockwood, Rouillard,
Eyles, & Harrison] for the 2008 December 12 CME event, as effectively measured
by the SWEPAM instrument aboard the Advanced Composition Explorer spacecraft
[Stone et al.(1998)Stone, Frandsen, Mewaldt, Christian, Margolies, Ormes, &
Snow]. We used the transformation (4) described in [Thompson(2006)]:

$x\simeq D_{\odot}\left(\frac{\pi}{180^{\circ}}\right)\theta_{x}$

to compute the displacement of the front, from the Helioprojective longitude
displacement of $\theta_{x}=0.27^{\circ}=1000$ seconds of arc from 18:09 to
18:49 of two algorithm-generated contours for HI-1 B, shown in Figure 4. Thus,
$\Delta x=730000$km, and $\Delta t=2400s$, so $v=300$km/s, using only two
significant digits.

Figure 4: Two-image sequence to estimate front speed, defined as the speed at
which the outermost edge advances from a snapshot to another one. The
calculation is detailed in the Results section.

Succesive images, such as in Figure 5 can allow more accurate estimations of
kinematics.

Figure 5: Longer sequence: successive images of a CME event advancing in the
field of vision of the STEREO-A camera. The front speed can be estimated from
such consecutive pairs and then averaged throughout, as explained in the
Results section.

##  4 Discussion

To improve the smoothness and completeness of the perimeter detection, we have
also tried a subsequent stage of the morphological geodesic active contour
method from [Márquez-Neila et al.(2014)Márquez-Neila, Baumela, & Alvarez],
based on level sets seeded with the interior of an initially-found contour. So
far, it did not seem to succeed, mainly because of the of the various bright-
dark alternating inner fronts, and also because of the background noise from
the light of the stars. The latter could be further reduced by carefully
realigning the images before subtraction, as suggested by [Davies et
al.(2009)Davies, Harrison, Rouillard, Sheeley, Perry, Bewsher, Davis, Eyles,
Crothers, & Brown].

A possible line of improvement would be to better adapt the active contour
methods to the already-found contour to fine-align it along less-contrasting
borders. A more precise way to distinguish the far-away background could also
help. Finally, I am planning on working to integrate it with existing models
to estimate CME velocities and 3D shape changes, and then systematically
comparing it with other methods.

##  5 Acknowledgements

We have implemented the algorithm in Python, and made it available at
https://github.com/MDNich/CME-Image-Edge-Detection-in-Python-from-STEREO-data,
being thus grateful to the community maintaining Numpy, SciPy and Matplotlib,
as well as [The SunPy Community et al.(2020)The SunPy Community, Barnes,
Bobra, Christe, Freij, Hayes, Ireland, Mumford, Perez-Suarez, Ryan, Shih,
Chanda, Glogowski, Hewett, Hughitt, Hill, Hiware, Inglis, Kirk, Konge, Mason,
Maloney, Murray, Panda, Park, Pereira, Reardon, Savage, Sipőcz, Stansby, Jain,
Taylor, Yadav, Rajul, & Dang]. We thank the reviewers for their careful
reading and helpful comments and for suggestions for applications of these
techniques to missions such as PUNCH and Solar Orbiter.

## References

  * [Byrne et al.(2009)Byrne, Gallagher, McAteer, & Young] Byrne, J. P., Gallagher, P. T., McAteer, R. T. J., & Young, C. A. 2009, A&A, 495, 325
  * [Canny(1986)] Canny, J. 1986, IEEE Transactions on Pattern Analysis and Machine Intelligence, 8, 679
  * [Caselles et al.(1993)Caselles, Catté, Coll, & Dibos] Caselles, V., Catté, F., Coll, T., & Dibos, F. 1993, Numer. Math., 66, 1
  * [Chan & Vese(2001)] Chan, T. F. & Vese, L. A. 2001, IEEE Transactions on Image Processing, 10, 266
  * [Crifo(1983)] Crifo, F. 1983, Solar Physics, 83, 143
  * [Davies et al.(2009)Davies, Harrison, Rouillard, Sheeley, Perry, Bewsher, Davis, Eyles, Crothers, & Brown] Davies, J. A., Harrison, R. A., Rouillard, A. P., et al. 2009, Geophys. Res. Lett., 36, L02102
  * [Davis et al.(2009)Davis, Davies, Lockwood, Rouillard, Eyles, & Harrison] Davis, C. J., Davies, J. A., Lockwood, M., et al. 2009, Geophys. Res. Lett., 36, L08102
  * [Eyles et al.(2007)Eyles, Davis, Harrison, Waltham, Halain, Mazy, Defise, Howard, Moses, Newmark, & Plunkett] Eyles, C., Davis, C., Harrison, R., et al. 2007, in Solar Physics and Space Weather Instrumentation II, ed. S. Fineschi & R. A. Viereck, Vol. 6689, International Society for Optics and Photonics (SPIE), 40 – 52
  * [Eyles et al.(2009)Eyles, Harrison, Davis, Waltham, Shaughnessy, Mapson-Menard, Bewsher, Crothers, Davies, Simnett, Howard, Moses, Newmark, Socker, Halain, Defise, Mazy, & Rochus] Eyles, C. J., Harrison, R. A., Davis, C. J., et al. 2009, Sol. Phys., 254, 387
  * [Fisher(1984)] Fisher, R. 1984, Adv. Space Res., 4, 163
  * [Fisher & Munro(1984)] Fisher, R. & Munro, H. R. 1984, Astrophys. J., 280, 163
  * [Gopalswamy & Yashiro(2007)] Gopalswamy, N. & Yashiro, S. 2007, Journal of Geophysical Research, 112
  * [Howard et al.(2008)Howard, Nandy, & Koepke] Howard, T., Nandy, D., & Koepke, A. 2008, Journal of Geophysical Research, 113
  * [Michalek et al.(2006)Michalek, Gopalswamy, Lara, & Yashiro] Michalek, G., Gopalswamy, N., Lara, A., & Yashiro, S. 2006, Space Weather, 4, S10003
  * [Márquez-Neila et al.(2014)Márquez-Neila, Baumela, & Alvarez] Márquez-Neila, P., Baumela, L., & Alvarez, L. 2014, IEEE Transactions on Pattern Analysis and Machine Intelligence, 36, 2
  * [Olmedo et al.(2008)Olmedo, Zhang, Wechsler, Poland, & Borne] Olmedo, O., Zhang, J., Wechsler, H., Poland, A., & Borne, K. 2008, in AGU Spring Meeting Abstracts, Vol. 2008, SP43A–02
  * [Robbrecht & Berghmans(2004)] Robbrecht, E. & Berghmans, D. 2004, A&A, 425, 1097
  * [Sheeley et al.(1997)Sheeley, Wang, Hawley, Brueckner, Dere, Howard, Koomen, Korendyke, Michels, Paswaters, Socker, St. Cyr, Wang, Lamy, Llebaria, Schwenn, Simnett, Plunkett, & Biesecker] Sheeley, N. R., Wang, Y. M., Hawley, S. H., et al. 1997, ApJ, 484, 472
  * [Stone et al.(1998)Stone, Frandsen, Mewaldt, Christian, Margolies, Ormes, & Snow] Stone, E. C., Frandsen, A. M., Mewaldt, R. A., et al. 1998, Space Sci. Rev., 86, 1
  * [The SunPy Community et al.(2020)The SunPy Community, Barnes, Bobra, Christe, Freij, Hayes, Ireland, Mumford, Perez-Suarez, Ryan, Shih, Chanda, Glogowski, Hewett, Hughitt, Hill, Hiware, Inglis, Kirk, Konge, Mason, Maloney, Murray, Panda, Park, Pereira, Reardon, Savage, Sipőcz, Stansby, Jain, Taylor, Yadav, Rajul, & Dang] The SunPy Community, Barnes, W. T., Bobra, M. G., et al. 2020, The Astrophysical Journal, 890, 68
  * [Thompson(2006)] Thompson, W. T. 2006, Astronomy and Astrophysics, 449, 791
  * [Webb & Howard(2012)] Webb, D. F. & Howard, T. A. 2012, Living Reviews in Solar Physics, 9
  * [Xie et al.(2006)Xie, Gopalswamy, Ofman, St. Cyr, Michalek, Lara, & Yashiro] Xie, H., Gopalswamy, N., Ofman, L., et al. 2006, Space Weather, 4, S10002
  * [Yashiro et al.(2004)Yashiro, Gopalswamy, Michalek, St. Cyr, Plunkett, Rich, & Howard] Yashiro, S., Gopalswamy, N., Michalek, G., et al. 2004, Journal of Geophysical Research (Space Physics), 109, A07105
  * [Young & Gallagher(2008)] Young, C. A. & Gallagher, P. T. 2008, Sol. Phys., 248, 457


# A simple, strong baseline for building damage detection on the xBD dataset

Sebastian Gerard  KTH Royal Institute of Technology, Stockholm, Sweden
<EMAIL_ADDRESS>Paul Borne-Pons  KTH Royal Institute of
Technology, Stockholm, Sweden
<EMAIL_ADDRESS>CentraleSupélec (Université Paris-Saclay),
Paris, France
<EMAIL_ADDRESS>Josephine Sullivan  KTH Royal Institute of
Technology, Stockholm, Sweden
<EMAIL_ADDRESS>
###### Abstract

We construct a strong baseline method for building damage detection by
starting with the highly-engineered winning solution of the xView2
competition, and gradually stripping away components. This way, we obtain a
much simpler method, while retaining adequate performance. We expect the
simplified solution to be more widely and easily applicable. This expectation
is based on the reduced complexity, as well as the fact that we choose
hyperparameters based on simple heuristics, that transfer to other datasets.
We then re-arrange the xView2 dataset splits such that the test locations are
not seen during training, contrary to the competition setup. In this setting,
we find that both the complex and the simplified model fail to generalize to
unseen locations. Analyzing the dataset indicates that this failure to
generalize is not only a model-based problem, but that the difficulty might
also be influenced by the unequal class distributions between events.

Code, including the baseline model:
https://github.com/PaulBorneP/Xview2_Strong_Baseline

##  1 Introduction

Figure 1: Published weights, our reproduction, our simplified model,
generalization. We compare the following models: The leftmost (brown) bar
represents the published weights of the competition winning solution. [1] The
$2^{nd}$ (dark green) bar uses the published code of the winning solution to
retrain the model on our hardware. We can see a drop in performance that is
not based on any intentional changes to the code. The $3^{rd}$ (dark blue) bar
shows our strong baseline model, derived from the winning solution by various
simplification steps. It performs slightly worse than our reproduction. The
two rightmost bar are, in order, the winning model (light green) and our
strong baseline (light blue), retrained on a data split where the test
disasters are not seen in training. Generalization proves difficult for both
models. Although the strong baseline yields worse results, this difference is
small, compared to the performance drop between the two splits. The drop is
especially steep for ’minor damage’ and ’major damage’. While the ’no damage’
and ’destroyed’ classes are easy to distinguish, it is difficult to clearly
distinguish the two damage levels in between, so seeing the performance drop
strongest in those two classes is not surprising. The results are based on
individual ResNet34-U-Net models.

With progressing climate change, extreme weather events are projected to occur
more often across many regions of the world. [2] Since these events typically
impact large areas, it is useful to incorporate large-scale observations from
remote sensing systems in the disaster response. Modern computer vision
methods can be used to process such large-scale observations automatically and
provide various analysis results in the area of humanitarian assistance and
disaster response (HADR).

As part of the xView2 competition, the xBD building damage prediction
dataset[3] was published. It contains high-resolution satellite image pairs
from before and after various types of disasters, paired with high-quality
building damage annotations for the post-disaster images. The natural starting
point for research on this kind of problem now is the solution that won the
xView2 competition. [1] However, since that solution was constructed with the
goal of winning a competition, it is highly engineered and fine-tuned in many
small ways, that make it hard to tweak and extend.

In this paper, we construct a simple, yet strong, baseline method by stripping
away components of the xView2-winning solution step by step. It retains most
of the performance of the original solution, while being much simpler, and
therefore easier to use and extend in further research. We then show that both
methods suffer greatly from generalization issues, when tested on a set of
events of xView2 that is not seen during training. This is contrary to the
original setup of the xView2 competition, in which images of each event are
seen during training and testing.

We publish the simplified baseline as a LightningModule, as well as providing
a LightningDataModule to load the xBD dataset itself. Both of these are part
of the PyTorch Lightning framework [4], that makes it easy to encapsulate all
relevant steps for training and inference into one object, instead of having
these details spread across various files and function calls, interspersed
with the training logic.

In summary, our contributions are the following:

  * •

We propose a strong baseline method based on the xView2-winning solution. We
demonstrate that the step-by-step removal of various components only has a
small influence on the performance of the method and the final method is only
about 2 percentage points worse than our reproduction of the original method.

  * •

We demonstrate that both the competition winner, as well as the strong
baseline, suffers from a strong generalization failure in two of the four
damage classes. This only becomes clear once we construct a dataset split, in
which the events do not overlap between train and test set.

  * •

We publish the strong baseline in an easily accessible form, to facilitate
further research.

Structure In section 2, we first describe the xBD dataset, including class
balance issues that might contribute to generalization problems. Then, we
describe our strong baseline model in section 3. In section 5, we show the
experimental results of simplifying the complex competition-winning model step
by step until we arrive at our baseline model. This includes the description
of the complexities of the winning solution. We move on to the question of
generalization in section 6, comparing the generalization performance of
competition-winning model and the simplified model on a non-overlapping
dataset split. section 6.

##  2 xBD: A Dataset for Assessing Building Damage from Satellite Imagery

(a) Pre-disaster

(b) Post-disaster

(c) Building damage annotation

Figure 2: xBD: Example images from the dataset The images show a location
impacted by the 2013 Moore tornado. The dataset contains one image from before
and one after the disaster. The annotation indicates the location of buildings
and the respective building damage class.

The xBD dataset [3] is a satellite image dataset for building damage
prediction. It contains images covering 22 disasters in 15 countries,
including hurricanes, tornadoes, wildfires, earthquakes, floods and volcano
outbreaks. It represents the task to detect buildings and classify them into
one of four building damage classes, given one satellite image from before and
one from after the disaster. The respective labels are high-quality
annotations created by humans. Each building is consistently labelled with one
class, even if only a part of the building is actually damaged. Figure 2 shows
example images. The xBD dataset was the basis for the xView2 challenge. For
training in the ablation experiments, we combine the train subset (2799 image
pairs) and the tier3 subset (6369 image pairs). In testing, we use the test
subset, containing 933 image pairs from events that are already part of the
training subsets. All images have a size of 1024² pixels.

###  2.1 Generalization difficulty on the dataset-level

In this section we explore dataset properties that make it harder for the
trained models to generalize well to unseen events. In Figure 3, we see how
strongly the class distribution differs between events. For example, wildfire
damages are rarely annotated as ’minor damage’ or ’major damage’, making it
very difficult to learn to predict this class for this disaster type. While
almost all distributions are dominated by the ’no damage’ class, more than
half of hurricane Matthew’s building pixels are labeled as ’minor damage’.
This is not the case for other hurricanes or tornadoes. This imbalance is
based on the dataset composition and makes it hard for the model to perform
equally well on all damage classes for all events, and also to generalize to
events that might have very different class distributions.

Figure 3: Damage class distributions per event. The distribution between
events and event categories can differ greatly. For example, wildfire damages
are mostly annotated as destroyed, while floods are mostly annotated as minor
or major damage. This imbalance adds to the difficulty of training a model to
perform well on all classes of all events.

Part of the generalization problem is also that the way that damage classes
look changes based on the geo-location. Minor damage looks different when it
affects a large, expensive house in Asia than when it affects low-income
housing in North America. This difficulty is inherent to the general task of
building damage classification itself. Additionally, the distinction between
damage classes from satellite images is no easy task in annotation either,
especially given the different appearance of these classes across geographies
and disaster types. While the ’no damage’ and ’destroyed’ classes are clearly
distinguishable, the distinction between ’minor damage’, ’major damage’ and
their respective neighboring classes is likely difficult, based only on the
satellite images. Labeling with in-situ visitations could improve upon this,
but this approach is of course completely impractical on the scale of this
dataset.

###  2.2 Assessing generalization via a new dataset split

To be able to investigate how well models generalize, we want to test them on
data that they have not seen during training. Since the original xBD dataset
splits have images from all test events in the training set, we need to create
a new split, that assigns each event exclusively to one of the subsets. At the
same time, we want to retain at least one disaster event per event type in
each of the subsets. For example, it does not make much sense to test on
wildfires if the model has never seen wildfires during training. When plotting
the images’ locations on a map, we further found that several of the disaster
events overlap spatially (see Figure 4) or are in close spatial proximity (see
Figure 5). To ensure a fair generalization test, we keep all subgroups of
spatially close events contained in the same subset. The new split is shown in
Table 1.

Figure 4: Spatial overlap between ”socal fire” and ”woolsey fire”, two disasters that can be found in our dataset and that actually describe the same fire that occurred around Los Angeles in November 2018. Map by OpenStreetMap. Figure 5: Spatial proximity between ”joplin tornado”, ”moore tornado” and ”midwestern floodings”. Those disaster occurred in geographically close areas in the state of Oklahoma and were kept in the same set as prevention. Map by OpenStreetMap. Train set |  | Test set
---|---|---
event name | count |  | event name | count
lower-puna-volcano | 291 |  | tuscaloosa-tornado | 343
palu-tsunami | 155 |  | guatemala-volcano | 23
mexico-earthquake | 159 |  | sunda-tsunami | 138
socal-fire | 1130 |  | santa-rosa-wildfire | 300
woolsey-fire | 878 |  | hurricane-matthew | 311
portugal-wildfire | 1869 |  | — | —
pinery-bushfire | 1845 |  | — | —
nepal-flooding | 619 |  | — | —
midwest-flooding | 359 |  | — | —
moore-tornado | 277 |  | — | —
joplin-tornado | 149 |  | — | —
hurricane-florence | 427 |  | — | —
hurricane-harvey | 427 |  | — | —
hurricane-michael | 441 |  | — | —
Table 1: New dataset splits To fairly evaluate the generalization performance
of the models, we rearrange the dataset to be free of overlaps between the
train and test set. This includes ensuring that disaster events that happened
in close proximity are contained in the same data subset. Those disaster are
in bold in the table

##  3 The strong baseline method

The strong baseline is a streamlined version of the xView2 challenge winning
solution [1], where we have kept the most important components of training and
architecture. The ablation experiments for the intermediate steps between
winning solution and simplified baseline are shown in section 5. The strong
baseline is built around an encoder-decoder network $f_{\text{\tiny enc-dec}}$
with skip connections and convolutional layers akin to a U-Net [5].
$f_{\text{\tiny enc-dec}}$ is parameterized by $\Theta_{1}$. This network is
applied independently to the before image $X_{\text{before}}$ and after image
$X_{\text{after}}$ of the disaster site, each of which has width $W$, height
$H$ and depth $D=3$:

$\displaystyle Z_{\text{before}}$ $\displaystyle=f_{\text{\tiny enc-
dec}}(X_{\text{before}},\Theta_{1}),$ $\displaystyle Z_{\text{after}}$
$\displaystyle=f_{\text{\tiny enc-dec}}(X_{\text{after}},\Theta_{1})$ (1)

with $Z_{\text{before}},Z_{\text{after}}\in\mathbb{R}^{W\times H\times
D_{1}}$. These outputs are concatenated
$Z_{\text{cat}}=[Z_{\text{before}},Z_{\text{after}}]$ and input into a final
small convolutional network $g$ consisting of $1\times 1\times 2D_{1}$
convolutions, parameterized by $\Theta_{2}$:

$\displaystyle Z$ $\displaystyle=g(Z_{\text{cat}},\Theta_{2})$ (2)

The final output is $Z\in\mathbb{R}^{W\times H\times\hat{C}}$, where $\hat{C}$
is the number of damage classes ($\hat{C}=4$: no damage, minor damage, major
damage, destroyed). $Z$ contains the score for each pixel of belonging to each
class. Each score in $Z$ is transformed to a scalar between 0 and 1 by
applying the sigmoid function:

$\displaystyle P$ $\displaystyle=\sigma(Z)$ (3)

where $\sigma(\cdot)$ denotes the sigmoid function. It is applied to each
entry of $Z$ independently and results in $P\in[0,1]^{W\times
H\times\hat{C}}$. One can interpret the entry $P_{ij\hat{c}}$ of $P$ as the
probability of pixel at location $(i,j)$ belonging to damage class $\hat{c}$.

Note: We represent the multi-class classification problem, where each pixel
should ultimately be assigned to just one class, as multiple binary
classification problems during training. While not strictly the correct
approach (as each pixel can have only one label), this modelling has been
shown [6, 7] to lead to better accuracy for supervised image classification.

###  3.1 Building localization

To discriminate building pixels from non-building pixels, we use a
localization mask $L\in\\{0,1\\}^{W\times H}$. We compute the final
segmentation mask $M\in\\{0,1,2,3,4\\}^{W\times H}$ as the output of our
method by masking the damage segmentation with the localization mask. Class
$0$ indicates the background class, while classes $1-4$ indicate the different
damage classes:

$\displaystyle M=(1+argmax(P_{ij\hat{c}}))\mathbb{1}_{L_{ij}=1}.$ (4)

The localization mask $L$ can be computed in a number of ways, all of which
yielded similar localization performance in our ablation studies:

  1. 1.

In the winning solution, the localization mask $L$ is created by thresholding
the prediction output $P_{loc}$ of a separate localization model, that has the
same architecture as the damage classification model.

  2. 2.

Our strong baseline predicts an additional class that distinguishes between
building and no-building, independently of the damage classification.

  3. 3.

Another approach we used on the way towards the strong baseline makes use of
the building damage predictions that are already part of the model. The
localization map can be created from the damage predictions by thresholding
them (we use OTSU’s method [8]). If at least one of the binary classifications
indicates that any level of building damage is present, it means there is also
a building at this position. This is possible, because we formulated the
problem as several independent binary classification problems, instead of
using a softmax function. We used this method for the early versions of our
solution, since a bug in the code led us to believe that the previously
mentioned method was not working well enough.

###  3.2 Loss function

The final version of our baseline predicts $C=5$ channels: four to predict
building damage and one to predict the presence or absence of buildings. We
keep the same notations as above but with $C$ instead of $\hat{C}$.

Learning the parameters, $\Theta=(\Theta_{1},\Theta_{2})$, of our complete
network from labelled training data requires optimizing a loss function w.r.t.
$\Theta$, given a training dataset. The loss in the winning solution is a
combination of several standard loss functions used in semantic segmentation.
Here, we review the mathematical definition of these losses and in turn also
introduce our notation.

Let $Y=(Y_{1},Y_{2},\ldots,Y_{C})$, where each $Y_{c}\in\\{0,1\\}^{W\times
H}$, for $1\leq c\leq C$, is the 2d binary matrix slice through the ground
truth $Y$, representing which pixel locations have label $c$. Similarly let
$P=(P_{1},P_{2},\ldots,P_{C})$ be the respective predictions.

For easier definition of the focal loss function, let $\tilde{P}_{c}$
represent the predicted probability for the presence or absence of class $c$
at each pixel location as indicated by $Y_{c}$:

$\displaystyle\tilde{P}_{c}$ $\displaystyle=(1-Y_{c})\odot(1-P_{c})+Y_{c}\odot
P_{c}$ (5)

where $\odot$ represents element-wise multiplication. For a ground truth
segmentation map $Y_{c}$ and a predicted segmentation map $P_{c}$ (and thus in
turn also $\tilde{P}_{c}$), the standard loss functions used to train the
strong baseline network for one input are:

$\displaystyle L_{\text{\tiny focal}}(P_{c},Y_{c})$
$\displaystyle=-\frac{1}{WH}\sum_{i,j}(1-\tilde{P}_{ijc})^{\gamma}\,\log(\tilde{P}_{ijc})$
Focal loss (6) $\displaystyle L_{\text{\tiny dice}}(P_{c},Y_{c})$
$\displaystyle=1-\frac{2\sum_{i,j}Y_{ijc}\,P_{ijc}}{\sum_{i,j}Y_{ijc}+\sum_{i,j}P_{ijc}}$
Soft dice loss (7)

where $\gamma\geq 0$. The combination loss, for non-negative scalars
$w_{\text{\tiny bce}},w_{\text{\tiny focal}},$ and $w_{\text{\tiny dice}}$, is
then defined as:

$\displaystyle L_{\text{\tiny combo}}(P_{c},Y_{c})=w_{\text{\tiny
focal}}L_{\text{\tiny focal}}(P_{c},Y_{c})+w_{\text{\tiny dice}}L_{\text{\tiny
dice}}(P_{c},Y_{c}).$ (8)

Finally, the overall training loss over all classes, for one labelled input
image pair is:

$\displaystyle L(P,Y)=\sum_{c=1}^{C}w_{c}\,L_{\text{\tiny
combo}}(P_{c},Y_{c}),$ (9)

where each $w_{c}$ is a non-negative scalar.

The strong baseline uses $\gamma=2$, $(w_{\text{\tiny focal}},w_{\text{\tiny
dice}})=(1,1)$ and $w_{c}=\frac{1}{f_{c}}$ for $c\in\\{1,2,3,4,5\\}$, where
$f_{c}$ is the relative frequency of class $c$ in the whole training dataset.
As channel 0 represents the presence of a building and not the presence of
background, the different classes are relatively balanced (compared to other
segmentation tasks).

###  3.3 Further implementation details

We use a single ResNet-34-based U-Net, which is the smallest among the winning
solution’s architectures. During training, the following geometric data
augmentations are applied to the whole images:

  * •

horizontal flipping

  * •

rotations of 0, 90, 180, 270 degrees;

  * •

random rotations in the range $[-10^{\circ},10^{\circ}]$ \+ random global
scaling in the range $[.9,1.1]$

For each labelled example the same augmentation is applied to before and after
images. Afterwards, random square patches of side length in the range
$[529,715]$ are cropped from the data-augmented images with a bias towards
patches that cover the building pixels from less frequent damage classes
(using inverse frequency) and resized to side length 608.

##  4 Experimental setup

The experiments are mainly based on the xView2 competition winner’s code [1].
It uses Python and PyTorch [9] to implement the models and all of the training
and evaluation code. We additionally use Weights&Biases [10] for tracking
experiments.

Hyperparameters For training, we used an AdamW optimizer, a learning rate of
2e-4 and weight decay of 1e-6. The learning rate was iteratively halved on
epoch 5, 11, 17, 23, 29 and 33. Furthermore, we used half-precision and
DataParallel training on three RTX 2080Ti GPUs with a global batch size of 14.

Evaluation metric The competition score is computed as a weighted average of
the localization F1 score $F1_{loc}$ and the harmonic mean of all building
damage F1 scores $\text{F1}_{\text{\tiny dmg}}=0.3~{}\text{F1}_{\text{\tiny
loc}}+0.7~{}\text{F1}_{\text{\tiny dmg}}$. Since localization and damage
estimation are evaluated separately, the damage predictions are only evaluated
where the ground truth labels show a building. Otherwise, localization errors
would influence the damage prediction score.

Dataset splits For the ablation studies in section 5, we use the train and
tier3 subsets as training data, and randomly split off 10% into a validation
set. The split is stratified over the disaster events, which is different from
the completely random split used in the winning solution. To compute the test
performance, we use the test set. The holdout set is not used. For the
generalization experiments, we use the non-overlapping splits described in
subsection 2.1. The validation dataset is split off in the same way as in the
ablation experiments.

##  5 Ablations: The path to a strong, but simple, baseline

Table 2: Performance of the winning solution for the xView2 challenge |  |  $F_{1}$ on main tasks |  |  $F_{1}$ score on damage classes
---|---|---|---|---
Model | Score | Local. | Damage |  | None | Minor | Major | Destroyed
Ensemble (4 architectures à 3 seeds) | 0.804 | 0.862 | 0.779 |  | 0.929 | 0.615 | 0.778 | 0.872
ResNet-34 (3 seeds) | 0.792 | 0.859 | 0.764 |  | 0.921 | 0.592 | 0.765 | 0.866
ResNet-34 (best seed) | 0.789 | 0.858 | 0.759 |  | 0.916 | 0.583 | 0.769 | 0.862
ResNet-34 (best seed, our rerun) | 0.770 | 0.854 | 0.734 |  | 0.911 | 0.544 | 0.741 | 0.860
Table 3: Ablation studies: We start from the re-run of the xView2 winner’s
solution (top row), and reduce its complexity step by step. We reach a much
simpler baseline approach, while the xView2 competition score only drops from
0.75 to 0.73. Our model does not use a separate finetuning step with a jump
back up in learning rate. We use a single model for localization and
classification, instead of two separate ones. This includes not validating the
model on the competition metric, while using the localization model, but on
the loss function instead. We use an equal-weighted combination of Focal and
soft Dice loss. We use class weights for the losses and crop-selection that
reflect the class distribution, instead of optimized weights. We only use
geometric augmentations. We do not dilate building classes in training, nor
the ’minor damage’ class in testing. The competition scores are the mean and
standard deviation over three repetitions with varying random seeds.
| Classif.
---
pretraining
| Separate
---
Finetuning
| Loc
---
method
| Tuned weights
---
in combo loss
function
| Optimized
---
loss
class weights
| Photom.
---
augs
| Dilation
---
train | test
| Optimized
---
cropping
class weights
| Competition
---
score
Localization | ✓ | Loc model | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | 0.75 ± 0.02
ImageNet | ✓ | Loc model | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | 0.75 ± 0.01
ImageNet | - | Loc model | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | 0.74 ± 0.0
ImageNet | - | OTSU | ✓ | ✓ | ✓ | ✓ | - | ✓ | 0.76 ± 0.01
ImageNet | - | OTSU | - | - | ✓ | ✓ | - | ✓ | 0.74 ± 0.01
ImageNet | - | OTSU | - | - | - | ✓ | - | ✓ | 0.74 ± 0.00
ImageNet | - | OTSU | - | - | - | - | - | ✓ | 0.73 ± 0.01
ImageNet | - | OTSU | - | - | - | - | - | - | 0.73 ± 0.01
ImageNet | - | Loc channel | - | - | - | - | - | - | 0.73 ± 0.01

We want to find out which components of the xView2 competition’s winning
solution (#1 solution) are important for its success, and which can be left
out, without losing a lot of performance.

Large ensemble The #1 solution consists of an ensemble of four architectures,
each trained on three random seeds. We investigate the effect of this ensemble
approach by using the published model weights [1] of the original author.
Table 2 shows that if we simply pick the best performing seed of the smallest
model, a ResNet-34-based U-Net, we only lose 1.5pp (percentage points) in the
competition score. This way, we greatly reduce the cost of training and
inference. In the following ablation experiments, we will retrain the model
with various changes to the initial setting.

To have a fair point of comparison, we re-run the training of the best
individual model among the published models. This model’s performance will be
used as the upper bound of performance that we compare all ablated model
versions to. Notably, the retrained model has 1.9pp lower performance than the
published weights. One factor that might cause this is the change from
randomly split train/val sets in the original code, to a stratified split, in
our experiment. Another factor is that the original hyperparameters might be
optimized for the used random seeds, which ensure reproducibility on the same
machine, but lead to different behavior when moving to our hardware.

Single model Starting from the ResNet-34 U-net, we remove different aspects of
the model, step by step. The results are shown in Table 3. The first row
represents the re-run of the best individual model as a reference point. We
will walk through the results here, to also explain the intricacies and tricks
used to squeeze out every last bit of performance, and win the xView2
competition. The resulting competition scores can be seen in Table 2, but
won’t be mentioned in the text, since the performance differences are all
rather small.

Pretraining & finetuning The original approach uses two separate models, one
for binary building localization, one for damage classification. The
localization model is trained first, and then used as a starting point for the
classification model. We use an ImageNet-pretrained model as the starting
point for classification instead. After training the classification model, the
#1 solution then performs something the original author calls tuning. This
consists of training for 3 more epochs, with a learning rate that is halved
every epoch and slightly changed augmentations. The learning rate is set to a
value that’s higher than the one at the end of the training, constituting a
slight warm restart. We remove this separate finetuning step. We increase the
number of training epochs from 20 to 40, to account for the fact that ImageNet
pretraining is less close to the building damage segmentation than the
localization pretraining. Additionally, we simplify the thresholding of the
localization prediction, which in the #1 solution is a combination of three
different threshold conditions, applied to localization and building damage
channels. We only retain the simplest threshold, which is a fixed threshold
applied only to the localization prediction.

Integrating localization and classification As a next step, we remove the
separate localization model. Instead, we use the building damage predictions,
that are modeled as four independent binary classification problems. If any of
them indicates the presence of building damage, the respective pixel must
contain a building. We apply OTSU-thresholding [8] to each building damage
channel and combine them with logical OR operations, to compute the
localization mask. At the same time, we also remove test-time dilation, which
dilates ’minor damage’ predictions to overwrite ’no damage’ predictions as a
post-processing step.

Tuned class and loss weights The #1 solution uses a weighted combination of
multiple loss functions (see subsection 3.2). For the ResNet34-based model,
these are focal and soft dice loss. For other architectures, they also mix in
a cross-entropy loss. We replace the weights for the loss functions with equal
weights for both losses. Furthermore, each per-class loss is weighted in a way
that is likely optimized via search and/or hand-tuned. We replace these
weights with the inverse relative frequency of each class in the training set.
The original class weights do not quite follow this distribution, and also
overemphasize the difficult ’minor damage’ class.

Training-time augmentations The #1 solution uses both geometric augmentations
(see subsection 3.3), as well as a range of photometric augmentations:
Channel-wise addition in RGB and HSV space, CLAHE, Gaussian noise, blurring,
and changes in saturation, brightness and contrast. We remove the photometric
augmentations and retain only the geometric ones. The #1 solution also dilates
all building damage classes during training, again preferring the ’minor
damage’ class if multiple dilated classes overlap. We remove this training-
time dilation. During random cropping, the #1 solution considers multiple
candidate crops and computes a preference score, again based on per-class
weights multiplied with the number of pixels of each class in the candidate
crop. The crop with the highest preference score is used. Oversampling
buildings makes sense, so we keep this aspect, but again replace the per-class
weights with distribution-based weights.

Localization as fifth channel Lastly, we replace the localization via OTSU
thresholding with simply predicting the presence of buildings as a separate
channel. This had been our first choice, since it is simple and clean, but due
to a problem in the code, it seemed like this option was not working well.
Instead of the adaptive OTSU threshold, we also go back to the simple fixed
threshold value set in the original solution. Both options should work in
practice. We make this last change and receive a much simpler method than the
one we started with, while only losing two percentage points in competition
score across all of the changes made.

###  5.1 Aspects we kept

During the ablation process, we attempted some modifications that resulted in
a severe drop in performance. Notably, when we removed the dice part of the
combo loss. The dice part in the loss appears to be necessary for learning
good building localization.

Another critical aspect was the importance of choosing the right batch size
during training. A batch size that is too small can negatively impact the
stability of training and ultimately affect the final score. This issue
becomes particularly problematic when using a weighted loss function, as some
classes may be underrepresented in the dataset. Consequently, computing the
loss for a few images can lead to encountering edge cases, hence making the
training unstable.

Other aspects that we did not try to change were the initial learning rate,
learning rate schedule and optimizer, as well as the fixed threshold used for
the localization prediction. The building damage predictions are aggregated
via argmax. We also did not try to change the relatively standard geometric
data augmentations.

##  6 Generalization: difficulties on the model-level

To assess the generalization performance of both models, and whether the
simpler model is maybe less overfit to the dataset, we train the model on the
new data splits shown in Table 1. We use the same hyperparameters as in
section 5, except that the class weights have been recomputed for the new
training set. The results are shown in the respective two rightmost bars in
the columns of Figure 1.

We see that both models overall perform much worse than on the original data
split. The localization, as well as the detection of the ’no damage’ and
’destroyed’ class are only a few percentage points worse than on the old
split, confirming that they are distinct enough to be easily distinguishable.
However, the performance on the ’minor damage’ class drops roughly 20% and the
’major damage’ class even drops roughly 50% for both models. These large gaps
in performance between the old, overlapping split, and the new, non-
overlapping split, demonstrates a big problem of generalization. Notably, the
’minor damage’ class is the one that received special attention in the xView2
winner’s solution. For example, it received a disproportionately high class
weight in training, and predictions on it were artificially dilated as a post-
processing step. This can explain some of the higher performance of the
winning solution. Comparing the simplified baseline with the winning solution,
the winning solution is again clearly better than the simplified model, but
the differences between the models’ performances is much smaller than the
generalization gap.

Based on the observations about class distributions differing between events
(see subsection 2.1), we also investigate if the performance on the individual
events in the test set are associated with the distance between either the
overall training set’s class distribution or the specific event type’s class
distribution in training. In both cases, we do not find a strong association.
The detailed results are not reported here.

##  7 Conclusion

We have shown that a greatly simplified version of the xView2 competition
winner’s approach still achieves an adequate performance. In the simplified
approach, various hyperparameter choices are made in clear, principled ways,
making it easy to train the model on different datasets or different splits of
the same dataset. Other parts of the complex approach are removed without
losing much performance. By retraining both models on a non-overlapping
dataset split, we find that generalization to unseen locations is a problem
for both models. Due to the simplified version, we know that this problem is
not caused by the winning model being overfit to the original dataset. After
showing that the class distributions vary greatly between events, we assume
that this imbalance in the dataset composition likely contributes to the
difficulty of generalizing to unseen events.

## References

  * [1] Victor Durnov  “xview2 1st place solution”, 2020  URL: https://github.com/vdurnov/xview2_1st_place_solution
  * [2] Intergovernmental Panel on Climate Change (IPCC)  “Weather and Climate Extreme Events in a Changing Climate”  In _Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change_ Cambridge University Press, 2023, pp. 1513–1766  DOI: 10.1017/9781009157896.013
  * [3] Ritwik Gupta, Richard Hosfelt, Sandra Sajeev, Nirav Patel, Bryce Goodman, Jigar Doshi, Eric Heim, Howie Choset and Matthew Gaston  “xBD: A Dataset for Assessing Building Damage from Satellite Imagery” arXiv:1911.09296 [cs]  arXiv, 2019  URL: http://arxiv.org/abs/1911.09296
  * [4] William Falcon and The PyTorch Lightning team  “PyTorch Lightning” doi: 10.5281/zenodo.3828935, https://www.pytorchlightning.ai, Last visited 2021-09-17, 2019  URL: https://www.pytorchlightning.ai
  * [5] Olaf Ronneberger, Philipp Fischer and Thomas Brox  “U-net: Convolutional networks for biomedical image segmentation”  In _International Conference on Medical image computing and computer-assisted intervention_ Springer, 2015, pp. 234–241
  * [6] Lucas Beyer, Olivier J. Hénaff, Alexander Kolesnikov, Xiaohua Zhai and Aäron van den Oord  “Are we done with ImageNet?” arXiv:2006.07159 [cs]  arXiv, 2020  URL: http://arxiv.org/abs/2006.07159
  * [7] Simon Kornblith, Ting Chen, Honglak Lee and Mohammad Norouzi  “Why do better loss functions lead to less transferable features?”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 28648–28662
  * [8] N. Otsu  “A Threshold Selection Method from Gray-Level Histograms” Conference Name: IEEE Transactions on Systems, Man, and Cybernetics  In _IEEE Transactions on Systems, Man, and Cybernetics_ 9.1, 1979, pp. 62–66  DOI: 10.1109/TSMC.1979.4310076
  * [9] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai and Soumith Chintala  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems 32_ Curran Associates, Inc., 2019, pp. 8024–8035  URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [10] Lukas Biewald  “Experiment Tracking with Weights and Biases”, 2020  URL: https://www.wandb.com/

##  Appendix A Contributions according to the CRediT system

Description of the CRediT system: https://credit.niso.org/

  * •

Sebastian Gerard: Conceptualization, Formal analysis, Investigation,
Methodology, Software, Writing – original draft, Writing – review & editing

  * •

Paul Borne-Pons: Formal analysis, Investigation, Methodology, Software,
Validation, Visualization, Writing – original draft

  * •

Josephine Sullivan: Conceptualization, Funding acquisition, Project
administration, Resources, Supervision, Writing – review & editing