# Spin-wave coupling to electromagnetic cavity fields in dysposium ferrite

M. Białek<EMAIL_ADDRESS>Institute of Physics, École Polytechnique
Fédérale de Lausanne (EPFL), 1015 Lausanne, Switzerland  A. Magrez  Institute
of Physics, École Polytechnique Fédérale de Lausanne (EPFL), 1015 Lausanne,
Switzerland  J.-Ph. Ansermet  Institute of Physics, École Polytechnique
Fédérale de Lausanne (EPFL), 1015 Lausanne, Switzerland

###### Abstract

Coupling of spin-waves with electromagnetic cavity field is demonstrated in an
antiferromagnet, dysprosium ferrite (DyFeO3). By measuring transmission at
0.2–0.35 THz and sweeping sample temperature, magnon-photon coupling
signatures were found at crossings of spin-wave resonances with Fabry-Pérot
cavity modes formed in samples. The obtained spectra are explained in terms of
classical electrodynamics and a microscopic model.

antiferromagnetism, antiferromagnetic resonance, magnons, electromagnons, spin
waves, optical cavities, cavity modes, light-matter coupling, THz, THz
spectroscopy, cavitronics

###### pacs:

71.36.+c, 76.50.+g, 78.30.-j, 75.50.Ee

Coupling of matter and electromagnetic radiation Mills and Burstein (1974) is
a topic of great interest in solid state physics research because of their
hybrid quantum nature. Khitrova  _et al._ (2006) In the THz range, phonon-
polaritons are a well-know example of a light-matter coupling. Recently,
polaritons in the THz region were shown with intersubband transitions,Anappara
_et al._ (2009); Todorov  _et al._ (2010) cyclotron resonance Bayer  _et al._
(2017) and plasmons Yu  _et al._ (2018) in two-dimensional electron gases, as
well as with intermolecular transitions in organic materials. Damari  _et al._
(2019) In these systems, a strong-coupling regime can be achieved when losses
are smaller than the exchange rate between light and matter, Paravicini-
Bagliani  _et al._ (2017) giving rise to the vacuum Rabi splitting. Polaritons
are composite particles, which are studied in basic research on quantum
optics,Khitrova  _et al._ (2006) and can be considered for use in quantum
computing and quantum memories.Wesenberg  _et al._ (2009); Schuster  _et al._
(2010); Kubo  _et al._ (2010); Huebl  _et al._ (2013); Frisk Kockum  _et al._
(2019) The coupling of electromagnetic cavity-modes to magnons was researched
intensively in ferromagnets at GHz frequencies, Huebl  _et al._ (2013); Yao
_et al._ (2015); Castel  _et al._ (2017) meeting with expectations of energy-
efficient spintronic devices. Matsukura  _et al._ (2015) Coupling of magnons
with superconducting qubits was also investigated. Tabuchi  _et al._ (2015);
Lachance-Quirion  _et al._ (2017) The Purcell enhancement and the vacuum Rabi
splitting were demonstrated in ferromagnetic materials.Zhang  _et al._ (2014);
Tabuchi  _et al._ (2014); Bai  _et al._ (2015); Maier-Flaig  _et al._ (2016)

It is interesting to investigate magnon-photon coupling in antiferromagnetic
materials Mills and Burstein (1974) in view of their high-frequency spin
dynamics,Kampfrath  _et al._ (2011); Jungwirth  _et al._ (2016); Baierl  _et
al._ (2016); Nemec  _et al._ (2018); Schlauderer  _et al._ (2019) comparing to
that of ferromagnets. As this phenomenon is readily taken into account by
classical electrodynamics, it was accounted for in the analysis of optical
investigations of antiferromagnetic materials, for instance in Ref. Balbashov
_et al._ (1985); Kozlov  _et al._ (1993); Mikhaylovskiy  _et al._ (2014,
2015a, 2015b). Some experimental reports focused on characterization of
interaction of antiferromagnetic magnons with photons in FeF2 Sanders  _et
al._ (1978); Jensen  _et al._ (1997), NiO Li  _et al._ (2011), TmFeO3
Grishunin  _et al._ (2018) and in ErFeO3. Li  _et al._ (2018); Sivarajah  _et
al._ (2019) Most of these reports use classical electrodynamics. A microscopic
picture was developed in Ref. Li  _et al._ , 2018. Here, we confront the
classical electrodynamic model with predictions of a microscopic model by
estimating, quantitatively from our data, the strength of the interaction of
antiferromagnetic spin-waves with electromagnetic cavity modes in dysprosium
ferrite DyFeO3 (DFO).

Dysprosium ferrite is an orthogonally-distorted perovskite. It shows
antiferromagnetic ordering below a Néel temperature of about 640 K White  _et
al._ (1982). Dzyaloshinskii-Moriya leads to a weak ferromagnetismIn DFO, the
spin canting allows two antiferromagnetic resonance modes to be excited: the
quasi-ferromagnetic (qFMR) and the quasi-antiferromagnetic (qAFMR),Zvezdin and
Matveev (1979); White  _et al._ (1982); Balbashov  _et al._ (1985); Koshizuka
and Hayashi (1988); Kozlov  _et al._ (1993) which are excited by the magnetic
component of radiation.Reid  _et al._ (2015)

Our measurements were performed above 400 K, where both resonances,
monotonously soften with rising temperature. We used polycrystaline samples
for their isotropic properties, as large anisotropic properties in crystals
might hinder the coupling. Białek  _et al._ (2019) We have created disk-shaped
samples of 1.0 and 0.6 mm thicknesses, thus having different cavity-mode
spectra. A sample was placed in a furnace, which allowed to control
temperature up to 700 K. By using quasioptical methods, we measured
transmission with our continuous-wave THz spectrometer based on frequency
extenders to a vector network analyzer (VNA). Caspers  _et al._ (2016) This
complex signal $S_{21}(f,T,H)$ is a function of frequency $f$, temperature $T$
and magnetic field $H$. We report its power amplitude in dB units and its
phase in degrees. In order to extract a signal related to magnetic resonances,
we measured transmission through a sample at different temperatures $T$. We
obtained temperature-differential spectra by subtracting averaged spectra
measured at subsequent temperatures:

$\frac{\partial S_{21}}{\partial T}=\Delta T^{-1}(S_{21}(f,T+\Delta
T,0)-S_{21}(f,T,0)),$ (1)

with $\Delta T=1$ K. We used this technique to measure spin-wave resonances in
bismuth ferrite at high Białek  _et al._ (2018) and low temperatures.Białek
_et al._ (2019)

In classical electrodynamics, the dispersion of electromagnetic radiation in a
material with a resonance is modified. This leads to creation of two polariton
states split by a frequency representing coupling of the resonance with the
electromagnetic radiation.Mills and Burstein (1974) In DFO, we model magnetic
susceptibility $\mu_{c}(f,T,H)$ using lorentzian distributions

$\mu_{c}(f,T,H)=1+\sum_{m=1}^{M}\frac{\Delta\mu_{m}(T,H)f_{m}^{2}(T,H)}{f_{m}^{2}(T,H)-f^{2}-if\gamma_{m}(T,H)},$
(2)

where, for the $m$-th magnetic resonance, $f_{m}$ is its frequency,
$\gamma_{m}$ its width and $\Delta\mu_{m}$ is its input to the zero-frequency
magnetic susceptibility. In the case of DFO, $M=2$, with $m=1$ corresponding
to the qFMR and $m=2$ to the qAFMR.

Ferrites have dielectric functions $\epsilon_{c}(f,T)$, which, in our
experimental frequency range, i.e. far from phonons resonant frequencies, can
be approximated as linear around a convenient point ($T_{0}=400$ K,
$f_{0}=0.3$ THz)

$\epsilon_{c}(f,T)=\epsilon_{00}+a(f-f_{0})+b(T-T_{0}).$ (3)

For our polycrystalline samples, we assume an effectiveBirchak  _et al._
(1974); Heller (1945) magnetic susceptibility as
$\sqrt{\mu}=p\sqrt{\mu_{c}}+(1-p)$ and effective dielectric function as
$\sqrt{\epsilon}=p\sqrt{\epsilon_{c}}+(1-p)$, where the factor $p=0.64$ is a
volume fraction of material to air in our pelletized samples.Białek  _et al._
(2019) This value is the maximum density of random-packed hard spheres.Song
_et al._ (2008) This assumption allows us to obtain values of $\epsilon_{c}$
and $\mu_{c}$ that are comparable with literature values for single crystal
samples.Kozlov  _et al._ (1993); Lobo  _et al._ (2007) The complex wave vector
is $k=2\pi f\sqrt{\epsilon\mu}/c$, where $c$ is the speed of light in vacuum.
Transmission of electric field $t(f,T,H)$ though a slab of a thickness $d$ and
infinite lateral dimensions is: Born and Wolf (1980)

$t(f,T,H)=\frac{(1-r^{2})e^{ikd}}{1-r^{2}e^{i2kd}},$ (4)

where
$r^{2}=(\sqrt{\epsilon}-\sqrt{\mu})^{2}/(\sqrt{\epsilon}+\sqrt{\mu})^{2}$ is
the square of the reflection coefficient at the vacuum-material interface. For
frequencies far away from resonances, Eq. 4 implies an interference pattern,
related to subsequent cavity modes of a slab, called Fabry-Pérot modes. At 300
GHz, in our polycrystalline DFO samples with the refractive index
$\sqrt{\epsilon\mu}\approx 3.6$, the wavelength is about 300 $\mu$m. With
rising temperature, $\sqrt{\epsilon\mu}$ increases, which shortens the period
of the interference pattern. The magnetic susceptibility takes into account
the effect of the magnetic resonance. When its frequency is close to a mode of
a cavity, both the magnon and the cavity lines are altered by the matter-
photon interaction. Since, the magnetic resonance frequencies have a much
stronger temperature dependence than that of the interference pattern, they
cross several interference pattern minima as temperature rises. We calculated

$\left\|\frac{\partial S_{21}}{\partial T}\right\|=\Delta
T^{-1}20log_{10}\left\|\frac{t(f,T+\Delta T,0)}{t(f,T,0)}\right\|$ (5)

to fit amplitude of temperature-differential spectra, where $\Delta T=1$ K is
the temperature step. We calculated phase of temperature-differential spectra
using

$\arg(\partial S_{21}/\partial T)=\arg(t(f,T+\Delta T,0)-t(f,T,0)).$ (6)

The experimental temperature-differential signal amplitude for the 1.0-mm-
thick DFO sample is shown in Fig. 1a.

Figure 1: Experimental temperature-differential spectra for the 1.0 mm-thick
DFO sample: (a) amplitude, (b) median-shifted phase with dashed black lines
showing coupled modes. Fit of the model to the amplitude data: (c) magnitude,
(d) phase. In the segment (d), green dashed lines show uncoupled spin-wave and
cavity modes and purple lines show coupled modes.

These data show clearly that the resonant lines are distorted when they cross
the sequence of sample-cavity modes. The amplitude and the widths of the
resonances are altered because of their interaction with the electromagnetic
standing waves. This is accounted for by Eq. 5 as shown in Fig. 1c. Thus, we
find that this model reproduces most of the important features of Fig. 1a. The
fitting parameters are parameters of the simplified dielectric function (Eq.
3) and of the magnetic susceptibility (Eq. 2). We assumed that resonance
frequencies have a temperature dependence described by a power law
$f_{m}^{*}(1-T/T_{N})^{\beta_{m}}$, applicable when approaching the Néel
temperature $T_{N}$,Eibschütz  _et al._ (1966) where $f_{m}^{*}$ has a unit of
frequency and $\beta_{m}\approx\frac{1}{3}$ in the case of DFO. To improve the
quality of our fits, we assumed that magnetic resonances have linear
dependences of their widths and amplitudes on temperature. These formulas and
values of fit parameters are given in the supplementary materials.

We measured the temperature-differential phase of transmitted electric field,
as presented in Fig. 1b. Despite high noise, features predicted in Fig. 1d are
observed in the experiment. The phase reveals clearly interactions between
electromagnetic waves and magnetization dynamics. In Fig. 2a, a spectrum
obtained at a temperature between crossings points shows that the qAFMR can be
accounted for with a harmonic model. In Fig. 2b, the spectrum obtained at the
temperature of the crossing with the $l=6$ cavity mode shows a structure that
cannot be explained using a single oscillator. It is due to two polariton
states which are not well-separated, i.e. they are in a weak coupling
regime.Khitrova  _et al._ (2006)

Figure 2: Examples of phase spectra extracted from Fig. 1b. (a) spectrum at a
temperature not showing light-matter interaction, (b) temperature close to a
crossing point. Arrows mark positions of cavity modes, qAFMR and upper and
lower polaritons. The violet lines show fits obtained using data in the entire
temperature range (Fig. 1d), thus better reflecting average properties of the
sample.

We used this phase prediction to estimate the cavity mode-magnetic resonance
coupling strength. Thus, we superimposed on the Fig. 1d predictions of the
harmonic coupling model Mills and Burstein (1974); Huebl  _et al._ (2013);
Harder  _et al._ (2016)

$f_{\pm}=\frac{1}{2}\left(f_{(l)}+f_{m}\pm\sqrt{(f_{(l)}-f_{m})^{2}+4\kappa^{2}f_{(l)}}\right),$
(7)

where $f_{\pm}$ indicates upper and lower polariton frequencies, $f_{(l)}$ is
the $l$-th cavity mode frequency, $f_{m}$ with $m=1,2$ are the resonances
frequencies. The coupling strength is given by Huebl  _et al._ (2013)

$\kappa=\frac{g_{s}\mu_{B}}{2h}\sqrt{\frac{\mu_{0}h}{2}p\rho},$ (8)

where $g_{s}=2$, $\rho$ is density of resonators and $p=0.64$ is the mass
filling factor of our polycrystalline sample. Song  _et al._ (2008); Białek
_et al._ (2019) Equation 8 assumes that a magnon is coupled to a single
electromagnetic cavity mode with a coupling strengthSoykal and Flatté (2010)
$g_{s}\mu_{B}B_{0}/2h$, where $B_{0}=\sqrt{\mu_{0}hf_{(l)}/2V_{(l)}}$ is the
magnetic component of vacuum fluctuations.Niemczyk  _et al._ (2009) Amplitude
of these is only $\approx 5\cdot 10^{-11}$ T for an electromagnetic mode at
0.3 THz taking the volume $V_{(l)}\approx 16\pi 10^{-9}$ m3 of the 1-mm-thick
sample. This is about 3 orders of magnitude smaller than the amplitude of the
THz field in our experiment. However, in an ensemble of $N$ resonators that
collectively interact with a cavity mode, coupling is increased by a factor
$\sqrt{N}$. Raizen  _et al._ (1989); Soykal and Flatté (2010); Huebl  _et al._
(2013); Harder  _et al._ (2016); Li  _et al._ (2018) Thus, the collective
coupling strength $\kappa$ depends only on oscillators density $\rho=N/V_{c}$
and the ratio $p$ of a crystalline volume $V_{c}$ to a cavity volume $V_{(l)}$
Huebl  _et al._ (2013); Li  _et al._ (2018). We take
$\rho=\frac{1}{3}\rho_{Fe}$, where $\rho_{Fe}=1.76\times 10^{28}$ m-3 is
density of iron atoms in DFO.Eibschütz (1965) This factor reflects the fact
that in a polycrystalline material, on average, only $1/3$ of magnetic dipoles
are excited by linear polarized electromagnetic wave. Under this assumption
$\kappa=1.78\times 10^{4}$ Hz1/2, which results in a splitting
$2\kappa\sqrt{f_{(l)}}\approx 19.5$ GHz with a cavity mode of $f_{(l)}=0.3$
THz.

Using Eq. 7, we calculated modes undergoing subsequent interactions with the
same coupling strength $\kappa$. We determined $f_{(l)}$ from the condition
$c(l-1/2)=2\Re{(\sqrt{\epsilon\mu})}f_{(l)}d$ for minimum of $\arg(\partial
S_{21}/\partial T)$. In the case of 1-mm-thick sample, the lowest visible mode
at $\approx 230$ GHz has $l=6$. Resonant modes frequencies were obtained from
the same fit to the experimental $\partial S_{21}/\partial T$ magnitude. For
the 1.0-mm-thick sample, the temperature dependence of both coupled and
uncoupled modes is presented in Fig. 1d.

In summary, we have observed coupling of spin-waves with electromagnetic
fields in high-temperature antiferromagnet DyFeO3. We observed the coupling of
THz-frequency magnetic resonances to modes of cavities formed by the samples
themselves. Our research demonstrates a possibility for a spin-cavitronic
system based on antiferrromagnetic resonances. To account for our results, we
applied both classical electrodynamics and a microscopic magnon-photon
coupling model,Huebl  _et al._ (2013) taking into account the polycrystalline
nature of our samples.

###### Acknowledgements.

We would like to thank prof. Can-Ming Hu of University of Manitoba for very
fruitful discussions and Claude Amendola for help in the samples preparation.
This work was partially supported by a Requip (206021150707/1) grant of the
Swiss National Science Foundation and by EPFL.

## References

  * Mills and Burstein (1974) D. L. Mills and E. Burstein, Reports on Progress in Physics 37, 817 (1974).
  * Khitrova  _et al._ (2006) G. Khitrova, H. M. Gibbs, M. Kira, S. W. Koch, and A. Scherer, Nature Physics 2, 81 (2006), review Article.
  * Anappara  _et al._ (2009) A. A. Anappara, S. De Liberato, A. Tredicucci, C. Ciuti, G. Biasiol, L. Sorba, and F. Beltram, Phys. Rev. B 79, 201303(R) (2009).
  * Todorov  _et al._ (2010) Y. Todorov, A. M. Andrews, R. Colombelli, S. De Liberato, C. Ciuti, P. Klang, G. Strasser, and C. Sirtori, Phys. Rev. Lett. 105, 196402 (2010).
  * Bayer  _et al._ (2017) A. Bayer, M. Pozimski, S. Schambeck, D. Schuh, R. Huber, D. Bougeard, and C. Lange, Nano Letters 17, 6340 (2017).
  * Yu  _et al._ (2018) Y. Yu, Z. Zheng, H. Qin, J. Sun, Y. Huang, X. Li, Z. Zhang, D. Wu, Y. Cai, B. Zhang, and V. V. Popov, Opt. Express 26, 31794 (2018).
  * Damari  _et al._ (2019) R. Damari, O. Weinberg, D. Krotkov, N. Demina, K. Akulov, A. Golombek, T. Schwartz, and S. Fleischer, Nature Communications 10, 3248 (2019).
  * Paravicini-Bagliani  _et al._ (2017) G. L. Paravicini-Bagliani, G. Scalari, F. Valmorra, J. Keller, C. Maissen, M. Beck, and J. Faist, Phys. Rev. B 95, 205304 (2017).
  * Wesenberg  _et al._ (2009) J. H. Wesenberg, A. Ardavan, G. A. D. Briggs, J. J. L. Morton, R. J. Schoelkopf, D. I. Schuster, and K. Mølmer, Phys. Rev. Lett. 103, 070502 (2009).
  * Schuster  _et al._ (2010) D. I. Schuster, A. P. Sears, E. Ginossar, L. DiCarlo, L. Frunzio, J. J. L. Morton, H. Wu, G. A. D. Briggs, B. B. Buckley, D. D. Awschalom, and R. J. Schoelkopf, Phys. Rev. Lett. 105, 140501 (2010).
  * Kubo  _et al._ (2010) Y. Kubo, F. R. Ong, P. Bertet, D. Vion, V. Jacques, D. Zheng, A. Dréau, J.-F. Roch, A. Auffeves, F. Jelezko, J. Wrachtrup, M. F. Barthe, P. Bergonzo, and D. Esteve, Phys. Rev. Lett. 105, 140502 (2010).
  * Huebl  _et al._ (2013) H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. B. Goennenwein, Phys. Rev. Lett. 111, 127003 (2013).
  * Frisk Kockum  _et al._ (2019) A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nature Reviews Physics 1, 19 (2019).
  * Yao  _et al._ (2015) B. M. Yao, Y. S. Gui, Y. Xiao, H. Guo, X. S. Chen, W. Lu, C. L. Chien, and C.-M. Hu, Phys. Rev. B 92, 184407 (2015).
  * Castel  _et al._ (2017) V. Castel, R. Jeunehomme, J. Ben Youssef, N. Vukadinovic, A. Manchec, F. K. Dejene, and G. E. W. Bauer, Phys. Rev. B 96, 064407 (2017).
  * Matsukura  _et al._ (2015) F. Matsukura, Y. Tokura, and H. Ohno, Nature Nanotechnology 10, 209 EP (2015), review Article.
  * Tabuchi  _et al._ (2015) Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Science 349, 405 (2015).
  * Lachance-Quirion  _et al._ (2017) D. Lachance-Quirion, Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, and Y. Nakamura, Science Advances 3 (2017), 10.1126/sciadv.1603150.
  * Zhang  _et al._ (2014) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Phys. Rev. Lett. 113, 156401 (2014).
  * Tabuchi  _et al._ (2014) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. Lett. 113, 083603 (2014).
  * Bai  _et al._ (2015) L. Bai, M. Harder, Y. P. Chen, X. Fan, J. Q. Xiao, and C.-M. Hu, Phys. Rev. Lett. 114, 227201 (2015).
  * Maier-Flaig  _et al._ (2016) H. Maier-Flaig, M. Harder, R. Gross, H. Huebl, and S. T. B. Goennenwein, Phys. Rev. B 94, 054433 (2016).
  * Kampfrath  _et al._ (2011) T. Kampfrath, A. Sell, G. Klatt, A. Pashkin, S. Mährlein, T. Dekorsy, M. Wolf, M. Fiebig, A. Leitenstorfer, and R. Huber, Nature Photonics 5, 31 (2011).
  * Jungwirth  _et al._ (2016) T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, Nature Nanotechnology 11, 231 EP (2016), review Article.
  * Baierl  _et al._ (2016) S. Baierl, M. Hohenleutner, T. Kampfrath, A. K. Zvezdin, A. V. Kimel, R. Huber, and R. V. Mikhaylovskiy, Nature Photonics 10, 715 EP (2016).
  * Nemec  _et al._ (2018) P. Nemec, M. Fiebig, T. Kampfrath, and A. V. Kimel, Nature Physics 14, 229 (2018).
  * Schlauderer  _et al._ (2019) S. Schlauderer, C. Lange, S. Baierl, T. Ebnet, C. P. Schmid, D. C. Valovcin, A. K. Zvezdin, A. V. Kimel, R. V. Mikhaylovskiy, and R. Huber, Nature 569, 383 (2019).
  * Balbashov  _et al._ (1985) A. M. Balbashov, A. A. Volkov, S. P. Lebedev, A. A. Mukhin, and A. S. Prokhorov, Soviet Physics - JETP 61, 543 (1985).
  * Kozlov  _et al._ (1993) G. V. Kozlov, S. P. Lebedev, A. A. Mukhin, A. S. Prokhorov, I. V. Fedorov, A. M. Balbashov, and I. Y. Parsegov, IEEE Transactions on Magnetics 29, 3443 (1993).
  * Mikhaylovskiy  _et al._ (2014) R. V. Mikhaylovskiy, E. Hendry, V. V. Kruglyak, R. V. Pisarev, T. Rasing, and A. V. Kimel, Phys. Rev. B 90, 184405 (2014).
  * Mikhaylovskiy  _et al._ (2015a) R. V. Mikhaylovskiy, E. Hendry, A. Secchi, J. H. Mentink, M. Eckstein, A. Wu, R. V. Pisarev, V. V. Kruglyak, M. I. Katsnelson, T. Rasing, and A. V. Kimel, Nature Communications 6, 8190 EP (2015a), article.
  * Mikhaylovskiy  _et al._ (2015b) R. V. Mikhaylovskiy, T. J. Huisman, A. I. Popov, A. K. Zvezdin, T. Rasing, R. V. Pisarev, and A. V. Kimel, Phys. Rev. B 92, 094437 (2015b).
  * Sanders  _et al._ (1978) R. Sanders, V. Jaccarino, and S. Rezende, Solid State Communications 28, 907 (1978).
  * Jensen  _et al._ (1997) M. R. F. Jensen, S. A. Feiven, T. J. Parker, and R. E. Camley, Journal of Physics: Condensed Matter 9, 7233 (1997).
  * Li  _et al._ (2011) J. Li, T. Higuchi, N. Kanda, K. Konishi, S. G. Tikhodeev, and M. Kuwata-Gonokami, Opt. Express 19, 22550 (2011).
  * Grishunin  _et al._ (2018) K. Grishunin, T. Huisman, G. Li, E. Mishina, T. Rasing, A. V. Kimel, K. Zhang, Z. Jin, S. Cao, W. Ren, G.-H. Ma, and R. V. Mikhaylovskiy, ACS Photonics 5, 1375 (2018).
  * Li  _et al._ (2018) X. Li, M. Bamba, N. Yuan, Q. Zhang, Y. Zhao, M. Xiang, K. Xu, Z. Jin, W. Ren, G. Ma, S. Cao, D. Turchinovich, and J. Kono, Science 361, 794 (2018).
  * Sivarajah  _et al._ (2019) P. Sivarajah, A. Steinbacher, B. Dastrup, J. Lu, M. Xiang, W. Ren, S. Kamba, S. Cao, and K. A. Nelson, Journal of Applied Physics 125, 213103 (2019).
  * White  _et al._ (1982) R. M. White, R. J. Nemanich, and C. Herring, Phys. Rev. B 25, 1822 (1982).
  * Zvezdin and Matveev (1979) A. K. Zvezdin and V. M. Matveev, Sov. Phys. JETP 50, 543 (1979).
  * Koshizuka and Hayashi (1988) N. Koshizuka and K. Hayashi, Journal of the Physical Society of Japan 57, 4418 (1988).
  * Reid  _et al._ (2015) A. H. M. Reid, T. Rasing, R. V. Pisarev, H. A. Dürr, and M. C. Hoffmann, Applied Physics Letters 106, 082403 (2015).
  * Białek  _et al._ (2019) M. Białek, T. Ito, H. Rønnow, and J.-P. Ansermet, Phys. Rev. B 99, 064429 (2019).
  * Caspers  _et al._ (2016) C. Caspers, V. P. Gandhi, A. Magrez, E. de Rijk, and J.-P. Ansermet, Applied Physics Letters 108, 241109 (2016).
  * Białek  _et al._ (2018) M. Białek, A. Magrez, A. Murk, and J.-P. Ansermet, Phys. Rev. B 97, 054410 (2018).
  * Birchak  _et al._ (1974) J. R. Birchak, C. G. Gardner, J. E. Hipp, and J. M. Victor, Proceedings of the IEEE 62, 93 (1974).
  * Heller (1945) W. Heller, Phys. Rev. 68, 5 (1945).
  * Song  _et al._ (2008) C. Song, P. Wang, and H. A. Makse, Nature 453, 629 EP (2008).
  * Lobo  _et al._ (2007) R. P. S. M. Lobo, R. L. Moreira, D. Lebeugle, and D. Colson, Phys. Rev. B 76, 172105 (2007).
  * Born and Wolf (1980) M. Born and E. Wolf, in  _Principles of Optics (Sixth Edition)_, edited by M. Born and E. Wolf (Pergamon, 1980) sixth edition ed., pp. 1 – 70.
  * Eibschütz  _et al._ (1966) M. Eibschütz, S. Shtrikman, and D. Treves, Solid State Communications 4, 141 (1966).
  * Harder  _et al._ (2016) M. Harder, L. Bai, C. Match, J. Sirker, and C. Hu, Science China Physics, Mechanics & Astronomy 59, 117511 (2016).
  * Soykal and Flatté (2010) O. O. Soykal and M. E. Flatté, Phys. Rev. Lett. 104, 077202 (2010).
  * Niemczyk  _et al._ (2009) T. Niemczyk, F. Deppe, M. Mariantoni, E. P. Menzel, E. Hoffmann, G. Wild, L. Eggenstein, A. Marx, and R. Gross, Superconductor Science and Technology 22, 034009 (2009).
  * Raizen  _et al._ (1989) M. G. Raizen, R. J. Thompson, R. J. Brecha, H. J. Kimble, and H. J. Carmichael, Phys. Rev. Lett. 63, 240 (1989).
  * Eibschütz (1965) M. Eibschütz, Acta Crystallographica 19, 337 (1965).


* This work was accepted to be presented at the IEEE/GRSS Workshop on Large Scale Computer Vision for Remote Sensing Imagery (EarthVision) to be held at the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) 2021.

# Machine-learned 3D Building Vectorization from Satellite Imagery

Yi Wang1 Stefano Zorzi2 Ksenia Bittner1
1German Aerospace Center (DLR)  2Graz University of Technology
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a machine learning based approach for automatic 3D building
reconstruction and vectorization. Taking a single-channel photogrammetric
_digital surface model (DSM)_ and a _panchromatic (PAN)_ image as input, we
first filter out non-building objects and refine the building shapes of the
input DSM with a _conditional generative adversarial network (cGAN)_. The
refined DSM and the input PAN image are then used through a semantic
segmentation network to detect edges and corners of building roofs. Later, a
set of vectorization algorithms are proposed to build roof polygons. Finally,
the height information from refined DSM is processed and added to the polygons
to obtain a fully vectorized _level of detail (LoD)_ -2 building model. We
verify the effectiveness of our method on large-scale satellite images, where
we obtain state-of-the-art performance.

##  1 Introduction

The availability of accurate 3D building models has become highly demanded in
various applications like the modeling of global urbanization process, urban
planning, disaster monitoring, . As traditional methods performed by human
operators for 3D building modeling are expensive, time-consuming and limited
to a small area, modern automatic 3D building model reconstruction methods
have drawn wide research interests.

Current automatic 3D building reconstruction methods can be generally
categorized into data-driven, model-driven and hybrid approaches. While model-
driven approaches extract the primitives of buildings and fit them to the most
appropriate models [17], data-driven methods extract geometrical components of
building roof planes from 3D point clouds or _digital surface models (DSMs)_
with point- or image-based segmentation techniques, and these components are
merged to 3D models with respect to some geometrical topology [27]. With
model-driven methods being unable to solve complex situations and data-driven
methods being commonly noisy, hybrid approaches, including this work, tend to
integrate the two types of approaches, where a data-driven approach extracts
the building components, and a model-driven approach utilizes prior knowledge
of the geometrical building models to help reconstruct 3D buildings [35].

Figure 1: Sample results of the proposed 3D building vectorization method.
$(a)$ refined DSM; $(b)$ edge and corner segmentation; $(c)$ vectorized 3D
building model.

While _light detection and ranging (LiDAR)_ point clouds and aerial images
have been the most common sources to extract 3D building information in the
past years [9, 10, 12], satellite images become more and more important as
they are convenient to acquire, cover wide areas and update frequently. Apart
from optical images, modern satellites can also provide DSMs using
photogrammetric stereo matching techniques, from which we can extract both
building objects and their height information. However, satellite DSMs show a
reasonable amount of noise and outliers because of matching errors or the
existence of non-building objects, thus refinement methods have been studied
to improve their quality. With traditional methods using filter-based
techniques like _principle component analysis (PCA)_ [20], Kalman filter [30]
and _fast Fourier transform (FFT)_ [3] to remove outliers, recent researches
have shown promising improvement by using deep learning based methods. Bittner
et al. [5] proposed firstly a _conditional generative adversarial network
(cGAN)_ based approach to filter out non-building objects and refine building
shapes of photogrammetric DSMs, which was further developed by a set of works
[6, 7, 8] to step-by-step improve the generation quality. Stucker and
Schindler [26] proposed an improvement for traditional stereo image matching
by regressing a residual correction with a convolutional neural network.

The revolutionary appearance of machine learning and deep learning techniques
has also brought significant contributions to the whole process of 3D building
reconstruction tasks. Not only building footprints can be extracted and
regularized with neural networks [28, 34, 36], but also the heights and roof
elements can be detected and predicted [1, 2], leading to constructed 3D
building models. Recent researches can be found in [21], where the authors
combined building object detection, semantic segmentation and height
prediction in a multi-task manner, and [29], where the authors proposed a deep
learning based model-driven approach to perform parametric building
reconstruction. While most of these researches focusing on _level of detail
(LoD)_ -1, LoD-2 building modeling is relatively new. One example is presented
in [24], where a hybrid 3D building reconstruction method is applied to detect
and decompose building boundaries, classify roof types, and fit predefined
building models.

Challenges for LoD-2 building reconstruction contain the requirement for
accurate building height prediction and roof element extraction, and the
complexity to form vectorized 3D roofs. Most existing methods utilize or
predict coarse height maps for detection tasks of neural networks and later
perform optimization [2, 24]. Our work, by contrast, uses network refined DSMs
to extract roof elements and proposes a corresponding vectorization pipeline
to form 3D models.

In this paper, we propose a machine learning based approach to reconstruct
LoD-2 building models from photogrammetric DSMs and _panchromatics (PANs)_
image obtained from satellites. Our contributions can be described as
following:

  * •

We improve the state-of-the-art cGAN based DSM refinement network proposed by
Bittner et al. [7] by adding a popular self-attention _convolutional block
attention module (CBAM)_ [32].

  * •

We propose an edge and corner detection network sharing the architecture of
the previous DSM refinement network.

  * •

We propose a novel vectorization pipeline to polygonize building roofs and
reconstruct 3D building models.

##  2 Methodology

As is shown in Figure 2, our multi-stage 3D building vectorization approach
starts with a cGAN architecture for photogrammetric DSM building shape
refinement. The refined DSM, together with the input PAN image, is then used
to detect building edges and corners with a semantic segmentation network that
shares the structure of the cGAN generator. The detected edges and corners are
later vectorized to building roof polygons. In the final stage, the refined
DSM and 2D polygons are combined to reconstruct 3D building models.

###  2.1 DSM building shape refinement

The proposed deep neural network for DSM refinement is an extension of the
network presented by Bittner et al. [7] based on an image-to-image translation
cGAN introduced by Isola et al. [14]. The network jointly learns a generator
and a discriminator to do the domain transfer, from a source domain, the
photogrammetric DSM, to a target domain, the refined DSM. With the
discriminator following the PatchGAN architecture proposed by Isola et al.
[14], the generator has a UNet-like structure with both long skip connections
from the encoders to the decoder and short skip connections in-between the
residual blocks inside the encoders. To enhance the feature of building
objects, we add a CBAM as presented by Woo et al. [32] before the decoder. The
CBAM is a combination of 1D channel attention and 2D spatial attention, which
are sequentially multiplied to the input feature maps. The overall generator
architecture is shown in Figure 3.

Figure 2: Overview of the proposed method. Given a photogrammetric DSM and a
PAN image as input, a cGAN based DSM refinement network and a semantic
segmentation network are sequentially applied to refine building shapes and
detect edges and corners. A set of vectorization algorithms are then applied
to reconstruct a full 3D building model. Figure 3: Generator architecture of
the proposed DSM refinement network.

Following the idea presented by Bittner et al. [7], we combine several types
of losses in a multi-task manner for optimizing the proposed DSM refinement
network:

$\begin{array}[]{c}\mathcal{L}_{\text{total}}\left(G\right)=\alpha\cdot\mathcal{L}_{\mathrm{cLSGAN}}(G,D)+\beta\cdot\mathcal{L}_{L_{1}}(G)\\\
+\gamma\cdot\mathcal{L}_{\text{normal
}}\left(\mathcal{N}^{\mathrm{t}},\mathcal{N}^{\mathrm{p}}\right)\end{array}$
(1)

where $\alpha$, $\beta$ and $\gamma$ represent the weighting parameters of
different loss terms.

GAN loss. We combine a conditional GAN [23] and a Least Squares GAN [22] for
the DSM refinement network, thus a _convolutional least squares GAN (cLSGAN)_
loss is utilized:

$\begin{array}[]{c}\underset{G}{\min}\underset{D}{\max}\mathcal{L}_{\text{cLSGAN
}}(G,D)=\mathbb{E}_{x,y\sim p_{\text{real }}(y)}\left[(D(y,x)-1)^{2}\right]\\\
+\mathbb{E}_{x,z\sim p_{z}(z)}\left[D(G(z,x),x)^{2}\right]\end{array}$ (2)

where $y\sim p_{\text{real }}(y)$ represents real samples, and $G(z)$
represents generated samples transferred from usually latent noise variables
$z\sim p_{z}(z)$. Respectively, $x$ denotes the _generative adversarial
network (GAN)_ ’s condition (the input DSM), $D(y,x)$ represents discriminator
output of real samples, and $D(G(z,x),x)$ represents discriminator output of
generated samples.

L1 loss. It is common to blend the objective functions for GANs with
traditional regression losses like $L1$ or $L2$ distances to help the
generator create images as close as possible to the given ground truth. Since
$L1$ loss encourages less blurring to the image boundaries, it is added to our
generator losses:

$\mathcal{L}_{L_{1}}(G)=\mathbb{E}_{x,y\sim p_{\text{real }}(y),z\sim
p_{z}(z)}\left[\|\mathrm{y}-G(z,x)\|_{1}\right]$ (3)

Normal vector loss. To further refine the surface of building roof planes, a
normal vector loss [13], which measures the angles between normal vectors of
generated and target DSMs, is added to the generator losses:

$\mathcal{L}_{\text{normal
}}\left(\mathcal{N}^{\mathrm{t}},\mathcal{N}^{\mathrm{p}}\right)=\frac{1}{m}\sum_{i=1}^{m}\left(1-\frac{\left\langle
n_{i}^{\mathrm{t}},n_{i}^{\mathrm{p}}\right\rangle}{\left\|n_{i}^{\mathrm{t}}\right\|\left\|n_{i}^{\mathrm{p}}\right\|}\right),$
(4)

where
$\mathcal{N}^{\mathrm{t}}=\left\\{n_{1}^{\mathrm{t}},\ldots,n_{m}^{\mathrm{t}}\right\\}$
and
$\mathcal{N}^{\mathrm{p}}=\left\\{n_{1}^{\mathrm{p}},\ldots,n_{m}^{\mathrm{p}}\right\\}$
represent normal vectors of the target and predicted DSM, and
$\langle\cdot,\cdot\rangle$ denotes the scalar product of the two vectors.
This normal vector loss emphasizes the planarity and inclination of building
roofs. The smaller the angle, the more planar the predicted surface and the
more consistent to the target surface.

The combination of different losses forms a multi-task learning problem, thus
an automatic weighting method proposed firstly by Kendal et al. [15] and
investigated in remote sensing in [18, 19] is applied to automatically tune
the loss weights considering the homoscedastic uncertainty of each separate
task:

$w_{l}=\left\\{\begin{array}[]{ll}0.5\cdot\exp\left(-\log\left(\sigma_{l}^{2}\right)\right)&\text{
for }\mathcal{L}_{L_{1}}\text{ and }\mathcal{L}_{\text{normal }}\\\
\exp\left(-\log\left(\sigma_{l}^{2}\right)\right)&\text{ for
}\mathcal{L}_{\text{cLSGAN }}\end{array}\right.$ (5)

where $\sigma_{l}^{2}$ is a learnable parameter, which represents the
variance, uncertainty of each task through the training process. In order to
avoid over-controlled parameter values, a regularization term
$0.5\cdot\log\left(\sigma_{l}^{2}\right)$ is added following each weighted
loss. As a result, the final loss of the generator of this DSM refinement
network can be formulated as:

$\mathcal{L}_{\text{total}}\left(G\right)=\sum_{l}w_{l}\cdot\mathcal{L}_{l}+\mathcal{R}_{l}$
(6)

while the discriminator loss remains the same as the cLSGAN loss:

$\begin{array}[]{c}\mathcal{L}_{\text{total }}(D)=\mathbb{E}_{x,y\sim
p_{\text{real }}(y)}\left[(D(y,x)-1)^{2}\right]\\\ +\mathbb{E}_{x,z\sim
p_{z}(z)}\left[D(G(z,x),x)^{2}\right]\end{array}$ (7)

Figure 4: Overview of the proposed vectorization pipeline. (a) detected edges
and corners. (b) edges and filtered corners. (c) vectorized edges and corners.
(d) roof polygons. (e) 3D building model.

###  2.2 Building edge and corner detection

Given the refined DSM and PAN image, a semantic segmentation network is used
to detect building edges and corners. The network architecture is identical to
the generator of the DSM refinement network (see Figure 3), except the change
of the three-channel output layer. A simple multi-class cross-entropy loss is
applied:

$\mathcal{L}_{\mathrm{CE}}(x,t)=\mathbb{E}\left[-\sum_{i=1}^{3}t_{i}\log
x_{i}\right]$ (8)

where $x_{i}$ is the predicted probability for a certain class $i$, and
$t_{i}$ is either 0 or 1 depending on the label of class $i$ for the
corresponding target. The output probability remains for further processing.

###  2.3 3D building model reconstruction

Figure 5: Examples of building edge vectorization.

In the final stage, a novel 3D building vectorization method is proposed using
the refined DSM and detected building edges and corners. Assuming building
edges are straight lines, the core idea is to step-by-step build a graph data
structure that stores points, lines, faces and their relationships for every
single building. As being a hybrid method, the proposed approach is not
limited to the complexity of different types of buildings, thus performing
well especially for large area 3D building modeling. A general workflow is
shown in Figure 4.

Corner point selection. For each corner pixel in the ground truth, multiple
surrounding pixels may be detected as corners, thus a _non-maximum suppression
(NMS)_ algorithm is implemented to filter out best fitting corner points. As
is shown in Figure 4 $(a)-(b)$, for each detected corner pixel (the
candidate), a surrounding $n\times n$ window is used as the evaluating box.
For each neighbor pixel in this window, if the pixel value (corner
probability) is no bigger than the candidate, it is set to zero; otherwise if
it exceeds the candidate, this pixel remains while the candidate is set to
zero. This process is iterated over all corner candidates and those isolated
best candidates are seen as final corner points.

Roof edge vectorization. Before we start the vectorization process, a
_connected component labeling (CCL)_ [33] algorithm is applied to label
connected pixels into building instances. Two pixels are connected when they
are neighbors and have a non-zero value. Here the neighborhood is defined in a
2-connected sense, which means every pixel has eight neighbors in eight
directions. As shown in Figure 5 $(a)$, different sets of connected pixels
would have different IDs and separate different buildings, which enables the
next steps to be performed within the scope of every single building.

Then we connect the corners to form edges based on two conditions. The first
condition is the average pixel value of a line buffer between a pair of corner
points. If the average value is above a threshold, an edge line is determined
between the corners. This condition would possibly fail when the edge is
curved in reality, thus a second condition is applied in parallel. By
utilizing the CCL algorithm again in a rectangle buffer between the pair of
corners, an edge is determined if the labels of the two corners are identical.
Two examples are shown in Figure 5 $(b)$, where both an edge with a hole and a
little curved edge can be successfully detected.

With the two conditions we can efficiently and thoroughly detect building
edges, yet still one problem needs to be considered. As it is shown in Figure
5 $(c)$, corner $A$ and $B$, corner $B$ and $C$ form two edges, but corner $A$
and corner $C$ can also form an edge which is redundant since it covers $AB$
and $BC$. To solve this issue, we again create a rectangle buffer for each
potential two-corner pair and, if other corner points exist inside this
buffer, this pair can not form an edge anymore.

Roof polygon generation. The vectorized edges are then polygonized to roof
faces ( see Figure 4 $(d)$), which can be easily done by graph search
algorithms. For each building, an undirected graph is firstly built from the
obtained edges. A simple _depth first search (DFS)_ is then applied to detect
and mark a cycle (a roof polygon) in this graph by tracing a back edge to
vertices that have been visited. This is run iteratively to extract all cycles
with corresponding different marks. To avoid face overlapping, large cycles
which cover small cycles are removed in the final step. In practice, the
polygonization process can also be directly applied with a $polygonize$
function from the open-source _shapely_ package which is popular for
manipulation and analysis of planar geometric objects [11].

Figure 6: The construction of final 3D building model. Height information from
refined DSM is processed and added to the polygons to build 3D roofs, walls
and ground face, together forming the final 3D model.

3D building modeling. In the final stage, walls and the ground face are
constructed utilizing roof polygons and the refined DSM to produce a full 3D
building model. Firstly, a _normalized digital surface model (nDSM)_ is
generated from the refined DSM with the method proposed by Qin et al. [25].
Then the adjacent roof faces are merged into a union, a polygon whose edges
are the building outlines. This gives us the footprints of the building, which
also means the 2D shape of the ground face. In the next step, the height
information from the nDSM is applied to corner points both inside roofs and on
the building boundaries. To avoid apparent height difference between endpoints
of an edge due to corner miss-matching (especially on outer boundaries where
corner height is supposed to be much bigger than neighboring ground pixels), a
small window is applied again to adjust height values. This is done by giving
the corner point the maximum height value in this small window. Though
slightly decreasing general accuracy, it can largely improve the robustness
and smoothness of resulting 3D models.

The edges of the union polygon represent both the upper and lower boundaries
of the building’s surrounding walls. With the height of upper corners already
determined with the maximum height value in the window, the height of lower
corners is determined by giving the minimum height value in the window, zero,
hence forming the building walls in 3D. Meanwhile, the lower edges form also
the ground face of the building, resulting in the final 3D building model. The
modeling process is shown in Figure 6.

##  3 Experiments and results

The proposed approach is evaluated on Worldview-1 data of Berlin, Germany. The
input consists of a space-borne photogrammetric DSM and a panchromatic image
with 0.5 $m$ spatial resolution covering a total area of 410 ${km}^{2}$. The
ground truth is generated from the public _city geography markup language
(CityGML)_ dataset following the same procedure as described in [4]. The
CityGML data for Berlin is freely available at
https://www.businesslocationcenter.de/en/economic-atlas/download-portal/. Open
datasets for some other worldwide cities can be found at
https://3d.bk.tudelft.nl/opendata/opencities/.

###  3.1 Implementation details

The DSM refinement network is based on the Coupled-UResNet cGAN architecture
proposed by Bittner et al.[7], with an additional CBAM [32] applied before the
decoder. The edge and corner detection network shares the architecture of the
generator of the DSM refinement network, while the last layer is changed to
three-channel output with a _$softmax$_ activation function.

The networks are trained on a single NVIDIA TITAN X (PASCAL) GPU with 12 GB
memory. To fit the training data into the GPU memory, the satellite images are
tiled into 21480 samples of size 256$\times$256 px. A minibatch of 4 is
applied in both networks. The samples are augmented not only by horizontal and
vertical flipping but also tiled from the original image with a random overlap
every epoch to give the model a clue about building geometries which happened
to be on the patch border in previous epochs. During the training of both
networks, the ADAM optimizer is used with an initial learning rate of
$\alpha=0.0002$ and momentum parameters ${\beta}_{1}=0.5,{\beta}_{2}=0.999$.
For the DSM refinement network, the generator is pre-trained for 100 epochs as
a warm-up and later interpolated with the cGAN’s generator. This so-called
network interpolation [31] can balance CNN’s over-smoothing and GAN’s over-
sharpening. The initial learnable weighting parameters as described in
Equation 5 are equally set to 1.

During the vectorization process, the window size for both corner point
filtering and corner height valuing is set to $5\times 5$ pix, while the width
for rectangle buffers (edge connecting and overlap elimination shown in Figure
5 $(c)$) is set to 7 pix.

###  3.2 Results and evaluation

(a) PAN image

(b) Photogrammetric DSM

(c) Ground truth DSM

(d) Refined DSM

(e) Detected edges and corners

(f) Vectorized edges and corners

Figure 7: Experimental results of a 500m $\times$ 500m testing area. Some
buildings in $(c)$ are not shown in other images because of the time
difference. Some edges are missing in $(f)$ compared to $(e)$ because they
don’t meet the requirements of vectorization process, especially for boundary
objects as they are incomplete.

Figure 7 $(c)$ shows the DSM refinement result, from which it can be seen that
the proposed network can both filter out and regularize building objects from
the photogrammetric DSM. This in parallel shows the robustness and accuracy of
our approach to detect correct buildings, as we can see from Figure 7 $(d)$
that the ground truth consists of several buildings that are not shown in
satellite images due to the time difference. _Mean absolute error (MAE)_ ,
_root mean squared error (RMSE)_ and _normalised median absolute deviation
(NMAD)_ are applied for quantitative evaluation of the DSM refinement result:

$\varepsilon_{\mathrm{MAE}}(h,\hat{h})=\frac{1}{n}\sum_{j=1}^{n}\left|\hat{h}_{j}-h_{j}\right|$
(9)

$\varepsilon_{\mathrm{RMSE}}(h,\hat{h})=\sqrt{\frac{1}{n}\sum_{j=1}^{n}\left(\hat{h}_{j}-h_{j}\right)^{2}}$
(10)

$\varepsilon_{\mathrm{NMAD}}(h,\hat{h})=1.4826\cdot{median}_{j}\left(\left|\Delta
h_{j}-m_{\Delta h}\right|\right)$ (11)

where $\hat{h}$ denotes the predicted heights, $h$ denotes the target heights,
$\Delta h$ denotes height error and $m_{\Delta h}$ denotes median height
error. As is shown in Table 1, our network improves all three metrics
evaluated over the testing area compared to Bittner et al. [7]. The RMSEs of
all DSMs are relatively large compared to the ground truth, which can be
explained by the time difference between the reference data and the given
satellite DSM. There can be cases when in one data source the buildings exist
and in the other not (due to new buildings construction or their destruction),
and vice versa.

Table 1: Quantitative metrics for refined DSM evaluated over the testing area.

| MAE (m) | RMSE (m) | NMAD (m)
---|---|---|---
Photogrammetric DSM | 3.91 | 7.14 | 1.40
Bittner et al. [7] | 1.73 | 4.02 | 0.93
Ours (with attention) | 1.42 | 3.65 | 0.60

Figure 7 $(e)$ and $(f)$ present the edge and corner detection and
vectorization results. By combining building height and shape information from
the refined DSM and intensity information from the PAN image, the results show
well-formed building skeletons with accurate corners and complete outlines. As
a result of the requirements from the vectorization process, edges which have
only one or none corner detected, or which are over-curved are unable to be
determined. However, though missing some of the expected line segments, most
of the building outer boundaries and inner edges are successfully constructed.
Meanwhile, it might be helpful to mention that during the experiments we tried
also combining the two steps (DSM refinement and edge and corner detection)
together in a multi-task manner, but the results got worse, as the edge and
corner detection network benefits more from an already refined DSM as input.

The final vectorized 3D building model is shown in Figure 8, where most of the
buildings are well reconstructed in 3D space. Even though some buildings are
not fully visible in PAN image and blurry in photogrammetric DSM, we can still
reconstruct them to a good shape. It is also seen that some buildings are
missing or incomplete, which is due to the missing of those vectorized edges
and corners whose quality doesn’t meet the vectorization process.

Figure 8: Reconstructed 3D building model of a $500m\times 500m$ testing area.

For quantitative evaluation of the height of reconstructed buildings, the
generated nDSM is compared to the ground truth. MAE, RMSE and NMAD are applied
again to evaluate the quality of the generated nDSM. The evaluation result is
shown in Table 2, from which we can see that both photogrammetric nDSM and our
genereated nDSM have better metrics than DSMs (Table 1) after removing the
height of ground surface. Meanwhile, our result presents large improvement
compared to photogrammetric nDSMs.

Table 2: Quantitative metrics for building nDSM evaluated over the testing
area.

| MAE (m) | RMSE (m) | NMAD (m)
---|---|---|---
Photogrammetric nDSM | 3.21 | 6.04 | 0.85
Ours | 0.80 | 2.28 | 0.47

Table 3: Quantitative metrics for roof orientation error evaluated over the
testing area.

| min (∘) | max (∘) | mean (∘) |  $\sigma$ (∘)
---|---|---|---|---
Photogrammetric nDSM | 0.08 | 75.84 | 22.46 | 22.28
Ours | 0.10 | 75.83 | 9.31 | 15.53

To evaluate the quality of the reconstructed 3D roofs, an orientation error is
applied to examine the inclination of the constructed roof planes. As proposed
by Koch et al. [16], the orientation error can be formulated as the angle
difference between the normal vectors of 3D planes fitted to the predicted
surface points and the given ground truth points:

$\varepsilon_{\text{orie
}}\left(G\odot\mathcal{P}\right)=\arccos\left(n_{i}^{\mathrm{t}}\cdot\check{n}_{i}^{\mathrm{p}}\right)$
(12)

where $n_{i}^{\mathrm{t}}$ and ${n}_{i}^{\mathrm{p}}$ denote the normal vector
of a certain plane on target and predicted image respectively.
$G\odot\mathcal{P}$ represents the predicted depth image $G$ masked with a
binary mask $\mathcal{P}$ containing a certain number of roof planes. Table 3
shows the average orientation error of constructed 3D roof faces compared to
corresponding ground truth, showing that the average plane angle is within
10∘, which is much better than using only the photogrammetric nDSM.

(a) PAN image

(b) Photogrammetric DSM

(c) Refined DSM

(d) Reconstructed 3D model

Figure 9: Testing results of a sub-area of Munich. Table 4: Comparison of eave
and ridge heights of the building model for selected buildings.

Building No. | Ridge (m) | Eave (m)
---|---|---
Reference | Partovi et al.[24] | Ours | Reference | Partovi et al.[24] | Ours
17 | 15 | 14.03 | 15.05 | 11 | 11.29 | 11.72
18 | 19 | 17.46 | 18.21 | 15 | 13.38 | 15.50
19 | 15 | 14.42 | 16.13 | 11 | 12.52 | 13.01
20 | 15 | 14.22 | - | 11 | 10.86 | -
21 | 15.5 | 14.08 | 15.33 | 11.9 | 12.21 | 11.54
22 | 15.6 | 15.28 | 14.87 | 11.5 | 11.87 | 11.94
23 | 20.0 | 20.76 | 21.80 | 16.5 | 17.35 | 17.58
24 | 16.2 | 15.87 | 17.03 | 12.3 | 11.03 | 13.66
25 | 17.4 | 16.21 | 18.02 | 13.6 | 13.58 | 13.77
26 | 16.8 | 16.40 | 17.19 | 12.5 | 10.54 | 11.36
27 | 15 | 14.66 | 13.88 | 10.9 | 10.49 | 11.40
28 | 16.8 | 16.27 | 17.11 | 12.5 | 10.41 | 12.96
29 | 14.7 | 13.94 | 15.54 | 10.7 | 10.46 | 11.06
30 | 16.8 | 15.66 | 16.00 | 12.5 | 10.51 | 13.79
$\mu_{|\Delta H|}$ | - | 0.79 | 0.74 | - | 0.93 | 0.99
$\sigma_{|\Delta H|}$ | - | 0.41 | 0.44 | - | 0.77 | 0.74
RMSE | - | 0.89 | 0.74 | - | 1.20 | 1.54
NMAD | - | 0.55 | 0.59 | - | 0.82 | 0.68

In addition, we compare our proposed 3D building vectorization method with the
work presented by Partovi et al. [24], who developed a multi-stage hybrid
method for 3D building reconstruction using PAN images, photogrammetric DSMs
and multi-spectral images from satellite data. Figure 9 presents the
reconstruction results of a sub-area of Munich using Worldview-2 satellite
data. The ridge and eave heights of 14 reconstructed buildings in this area
are compared with reference data from the Department of Urban Planning and
Building (DUPB) of Munich. As is shown in Table 4, $|\Delta H|$ denotes the
absolute height difference between the predicted model and reference, and
$\mu_{|\Delta H|}$ and $\sigma_{|\Delta H|}$ represent the mean and standard
deviation of the height difference, respectively. The building numbers refer
to Figure 9 $(c)$. It can be seen that both methods lead to lower accuracy in
eave heights than ridge heights because the surroundings of building
boundaries are usually more complex than inner-roof ridges both in PAN image
and photogrammetric DSM. Our method tends to get bigger values for eave
heights, which can be explained by our valuing method for the height of
building corners. In order to avoid the mismatching between DSM heights and
corner positions, we give the eave corner the maximum height value in a
surrounding window and the minimum height value for the corresponding ground
corner. This would increase the relative height of the building eaves, yet we
can see that this systematic error is within a small range. Apart from that,
the overall accuracy shows promising superiority of our method, where we get
comparative metric performance with a simpler approach. Meanwhile, as a price
of simplicity, the biggest problem remaining to be solved is the lack of
completeness of our constructed model. As can be seen from both Figure 8 and
Figure 9, some building components are lost after vectorization, which
quantitatively reduces the recall score from $0.88$ to $0.81$ (Berlin testing
area) compared to the refined DSM before vectorization.

##  4 Conclusion

In this paper, we present a multi-stage large-scale 3D building vectorization
approach. We extend the application of recent deep learning based techniques
on photogrammetric _digital surface model (DSM)_ refinement and bring it to
the application of automatic 3D building model reconstruction. With the help
of a self-attention module, we obtain promising results for both regression of
building heights and semantic segmentation of edges and corners. Based on
that, we propose a simple yet effective vectorization pipeline to reconstruct
_level of detail (LoD)_ -2 building models. We apply _non-maximum suppression
(NMS)_ to filter out best fitting corner points, define buffer connectivity
and buffer thresholds to determine edges, and polygonize them to roof faces.
By utilizing again the height information from the refined DSM, we finally
reconstruct fully vectorized 3D building models. Though limitations exist in
straight edge assumptions and the completeness of reconstructed building
models, results prove the overall robustness and accuracy of our proposed
method.

## References

  * [1] Fatemeh Alidoost, Hossein Arefi, and Michael Hahn.  Y-shaped convolutional neural network for 3d roof elements extraction to reconstruct building models from a single aerial image.  ISPRS Annals of Photogrammetry, Remote Sensing & Spatial Information Sciences, 5(2), 2020.
  * [2] Fatemeh Alidoost, Hossein Arefi, and Federico Tombari.  2d image-to-3d model: knowledge-based 3d building reconstruction (3dbr) using single aerial images and convolutional neural networks (cnns).  Remote Sensing, 11(19):2219, 2019.
  * [3] Katherine Arrell, Steve Wise, Jo Wood, and Danny Donoghue.  Spectral filtering as a method of visualising and removing striped artefacts in digital elevation data.  Earth Surface Processes and Landforms: The Journal of the British Geomorphological Research Group, 33(6):943–961, 2008.
  * [4] Ksenia Bittner, Pablo d’Angelo, Marco Körner, and Peter Reinartz.  Dsm-to-lod2: Spaceborne stereo digital surface model refinement.  Remote Sensing, 10(12):1926, 2018.
  * [5] Ksenia Bittner and Marco Korner.  Automatic large-scale 3d building shape refinement using conditional generative adversarial networks.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, pages 1887–1889, 2018.
  * [6] Ksenia Bittner, Marco Körner, Friedrich Fraundorfer, and Peter Reinartz.  Multi-task cgan for simultaneous spaceborne dsm refinement and roof-type classification.  Remote Sensing, 11(11):1262, 2019.
  * [7] Ksenia Bittner, Lukas Liebel, Marco Körner, and Peter Reinartz.  Long-short skip connections in deep neural networks for dsm refinement.  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 43(B2):383–390, 2020.
  * [8] Ksenia Bittner, Peter Reinartz, and Marco Korner.  Late or earlier information fusion from depth and spectral data? large-scale digital surface model refinement by hybrid-cgan.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pages 0–0, 2019.
  * [9] Claus Brenner.  Building reconstruction from images and laser scanning.  International Journal of Applied Earth Observation and Geoinformation, 6(3-4):187–198, 2005.
  * [10] Liang-Chien Chen, Tee-Ann Teo, J-Y Rau, Jin-King Liu, and Wei-Chen Hsu.  Building reconstruction from lidar data and aerial imagery.  In Proceedings. 2005 IEEE International Geoscience and Remote Sensing Symposium, 2005. IGARSS’05., volume 4, pages 2846–2849. IEEE, 2005.
  * [11] Sean Gillies et al.  Shapely: manipulation and analysis of geometric objects, 2007–.
  * [12] Norbert Haala and Martin Kada.  An update on automatic 3d building reconstruction.  ISPRS Journal of Photogrammetry and Remote Sensing, 65(6):570–580, 2010.
  * [13] Junjie Hu, Mete Ozay, Yan Zhang, and Takayuki Okatani.  Revisiting single image depth estimation: Toward higher resolution maps with accurate object boundaries.  In 2019 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 1043–1051. IEEE, 2019.
  * [14] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros.  Image-to-image translation with conditional adversarial networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1125–1134, 2017.
  * [15] Alex Kendall, Yarin Gal, and Roberto Cipolla.  Multi-task learning using uncertainty to weigh losses for scene geometry and semantics.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7482–7491, 2018.
  * [16] Tobias Koch, Lukas Liebel, Friedrich Fraundorfer, and Marco Korner.  Evaluation of cnn-based single-image depth estimation methods.  In Proceedings of the European Conference on Computer Vision (ECCV) Workshops, pages 0–0, 2018.
  * [17] Florent Lafarge, Xavier Descombes, Josiane Zerubia, and Marc Pierrot-Deseilligny.  Structural approach for building reconstruction from a single dsm.  IEEE Transactions on pattern analysis and machine intelligence, 32(1):135–147, 2008.
  * [18] Lukas Liebel, Ksenia Bittner, and Marco Körner.  A generalized multi-task learning approach to stereo dsm filtering in urban areas.  ISPRS Journal of Photogrammetry and Remote Sensing, 166:213–227, 2020.
  * [19] Lukas Liebel and Marco Körner.  Auxiliary tasks in multi-task learning.  arXiv preprint arXiv:1805.06334, 2018.
  * [20] Carlos López.  Improving the elevation accuracy of digital elevation models: a comparison of some error detection procedures.  Transactions in GIS, 4(1):43–64, 2000.
  * [21] Jisan Mahmud, True Price, Akash Bapat, and Jan-Michael Frahm.  Boundary-aware 3d building reconstruction from a single overhead image.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 441–451, 2020.
  * [22] Xudong Mao, Qing Li, Haoran Xie, Raymond YK Lau, Zhen Wang, and Stephen Paul Smolley.  Least squares generative adversarial networks.  In Proceedings of the IEEE international conference on computer vision, pages 2794–2802, 2017.
  * [23] Mehdi Mirza and Simon Osindero.  Conditional generative adversarial nets.  arXiv preprint arXiv:1411.1784, 2014.
  * [24] Tahmineh Partovi, Friedrich Fraundorfer, Reza Bahmanyar, Hai Huang, and Peter Reinartz.  Automatic 3-d building model reconstruction from very high resolution stereo satellite imagery.  Remote Sensing, 11(14):1660, 2019.
  * [25] Rongjun Qin, Jiaojiao Tian, and Peter Reinartz.  Spatiotemporal inferences for use in building detection using series of very-high-resolution space-borne stereo images.  International Journal of Remote Sensing, 37(15):3455–3476, 2016\.
  * [26] Corinne Stucker and Konrad Schindler.  Resdepth: Learned residual stereo reconstruction.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pages 184–185, 2020.
  * [27] Fayez Tarsha-Kurdi, Tania Landes, Pierre Grussenmeyer, and Mathieu Koehl.  Model-driven and data-driven approaches using lidar data: Analysis and comparison.  In ISPRS workshop, photogrammetric image analysis (PIA07), pages 87–92, 2007.
  * [28] Maria Vakalopoulou, Konstantinos Karantzalos, Nikos Komodakis, and Nikos Paragios.  Building detection in very high resolution multispectral data with deep learning features.  In 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), pages 1873–1876. IEEE, 2015.
  * [29] Ke Wang and Jan-Michael Frahm.  Single view parametric building reconstruction from satellite imagery.  In 2017 International Conference on 3D Vision (3DV), pages 603–611. IEEE, 2017.
  * [30] Ping Wang.  Applying two dimensional kalman filtering for digital terrain modelling.  Proceedings of International Archives of Photogrammetry, Remote Sensing, and Spatial Information Sciences, pages 649–656, 1998.
  * [31] Xintao Wang, Ke Yu, Shixiang Wu, Jinjin Gu, Yihao Liu, Chao Dong, Yu Qiao, and Chen Change Loy.  Esrgan: Enhanced super-resolution generative adversarial networks.  In Proceedings of the European Conference on Computer Vision (ECCV) Workshops, pages 0–0, 2018.
  * [32] Sanghyun Woo, Jongchan Park, Joon-Young Lee, and In So Kweon.  Cbam: Convolutional block attention module.  In Proceedings of the European conference on computer vision (ECCV), pages 3–19, 2018.
  * [33] Kesheng Wu, Ekow Otoo, and Arie Shoshani.  Optimizing connected component labeling algorithms.  In Medical Imaging 2005: Image Processing, volume 5747, pages 1965–1976. International Society for Optics and Photonics, 2005.
  * [34] Kang Zhao, Jungwon Kang, Jaewook Jung, and Gunho Sohn.  Building extraction from satellite images using mask r-cnn with building boundary regularization.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, pages 247–251, 2018.
  * [35] Yuanfan Zheng, Qihao Weng, and Yaoxing Zheng.  A hybrid approach for three-dimensional building reconstruction in indianapolis from lidar data.  Remote Sensing, 9(4):310, 2017.
  * [36] Stefano Zorzi, Ksenia Bittner, and Friedrich Fraundorfer.  Machine-learned regularization and polygonization of building segmentation masks.  arXiv preprint arXiv:2007.12587, 2020.