# Classical versus quantum calculation of radiative electric quadrupole
transition rates for hydrogenic states

Michael Horbatsch  Department of Physics and Astronomy, York University,
Toronto, Ontario M3J 1P3, Canada  Marko Horbatsch<EMAIL_ADDRESS>Department
of Physics and Astronomy, York University, Toronto, Ontario M3J 1P3, Canada

###### Abstract

The semiclassical Kepler-Coulomb problem and the quantum-mechanical
Schrödinger-Coulomb problem are compared for their predictions of quadrupole
E2 transitions. The semiclassical treatment involves an extension of previous
work for the electric dipole transitions (Physical Review A 71, 020501), and
rates are derived for $\Delta\ell=0,\pm 2$ transitions on the basis of the
multipolar properties of the emitted radiation. For the quantum case a
derivation is presented within the Schrödinger framework without reference to
spin. Comparison of the E2 rates shows reasonable agreement, but not as good
as was found for the electric dipole case.

##  I Introduction

The electromagnetic radiation problem for an accelerated charged particle is
known in the physics curriculum as the Larmor problem Landau and Lifshitz
(1975); Jackson (1999a). The dominant dipolar contribution is treated in many
textbooks, as is the associated problem of the radiation reaction force which
can be introduced to account for the transfer of mechanical energy due to the
radiating charge. In the context of quantum mechanics the radiation that is
associated with transitions from excited to lower states can be motivated by
the classical radiation problem Corney (2006), but once one introduces matrix
elements that connect initial and final states to calculate transition rates
for particular pairs of states the connection to the classical problem appears
to be lost or at least deeply hidden.

For the dominant electric dipole transitions from an initial state
characterized by principal and orbital angular momentum numbers
$(n_{i},\ell_{i})$ to a final state labeled by $(n_{f},\ell_{f})$ a reasonably
accurate quantum-classical correspondence has been established in Ref.
Horbatsch  _et al._ (2005a). In this work the semiclassical Kepler-Coulomb
trajectories for the initial state were represented through their Fourier-
Bessel series, in principle, as was done originally by Heisenberg and Kramers,
just before the advent of modern quantum mechanics. One modification, however,
was the use of semiclassical quantization in the sense of the WKB
approximation. The equidistant frequency spectrum arising from the Fourier-
Bessel series was not interpreted literally, but a re-scaling of the
contributions from (infinitely) many integer orders was performed to a non-
integer index which took into account the discrete $n$-dependent Balmer
spectrum of final available states for the mechanical system. This re-scaled
approach, in which one uses the fact that $n_{f}=1$ represents the ground
state which is not allowed to decay further resulted in quite accurate
transition rates as compared to the Schrödinger results, if one made the
restriction that $\ell_{i}>0$. In fact, it was possible to derive branching
ratios associated with angular-momentum-decreasing
($\Delta\ell\equiv\ell_{f}-\ell_{i}=-1$) and also angular-momentum-increasing
($\Delta\ell=+1$) transitions.

In this way it was possible to establish with more confidence than was
provided by the general statements in the literature that using radiated power
as predicted by the Larmor expression and dividing by the photon energy,
accurate transition rates could be obtained for electric dipole (E1) allowed
transitions. This work was also extended to the case of hydrogenic diamagnetic
states Horbatsch  _et al._ (2005b), for which Hamilton-Jacobi perturbation
theory had been developed previously Delos  _et al._ (1983). More recently
this work was extended to the Stark states Horbatsch and Horbatsch (2021). In
both cases the Kepler-Coulomb ellipses are perturbed by the external field,
yet the classical-quantum correspondence remains largely intact.

Electric quadrupole (E2) transitions of hydrogenic states (and magnetic dipole
M1, in particular) are usually addressed within a relativistic approach based
on solutions to the Dirac equation Jitrik and Bunge (2004). For our purposes
of connecting with the semiclassical Coulomb-Kepler problem it is, however,
sufficient to deal with them at the level of Schrödinger-Pauli theory, and if
we ignore the magnetic dipole case, we should really compare to Schrödinger
solutions, i.e., ignore fine structure. Electric quadrupole transition rates
are of interest for cases where electric dipole transitions are not allowed,
e.g., by a sequence of two such transitions via some intermediate state. The
E2 selection rules, which are discussed, e.g., in Refs. Bethe and Salpeter
(1957); Weissbluth (1978), then require us to look at the dominant case of
$\Delta\ell=-2$, the case of $\Delta\ell=0$, and $\Delta\ell=+2$.
Interestingly, the calculation of an electric quadrupole transition does come
up in an undergraduate textbook as a problem (Problem 11.31 in Ref. Griffiths
and Schroeter (2018)). Some discussion of how to calculate these transition
rates is found in the classic book of Bethe and Salpeter Bethe and Salpeter
(1957), but it comes up short in providing a complete treatment, since it only
considers the case of radiation emitted in particular directions. Thus, our
section on how to calculate these rates in quantum mechanics can fill a gap on
its own.

The classical quadrupole radiation problem plays a dominant role in
gravitational physics Peters and Mathews (1963), but the situation is similar
to the classical electromagnetic problem, yet not identical. In the
electromagnetic Kepler problem the appearance of opposite charges leads to a
time-dependent electric dipole moment whose acceleration results in E1 Larmor
radiation. In the gravitational case where the mass plays the role of charge
the center of mass does not accelerate in the two-body problem, and thus
dipole radiation is not allowed Landau and Lifshitz (1981). Quadrupole
radiation can be associated with the third time derivative of the quadrupole
moment, e.g., as represented by the moment of inertia, and it arises in both
Kepler problems. In the Kepler-Coulomb problem it is just a correction to
dipole radiation, in the gravitational Kepler problem it becomes the dominant
radiative process. Motivation for studying the subject of the present paper is
of interest in this context, since eventually a quantum theory of gravity will
need to connect quantum and classical gravitational radiation for objects or
situations where the classical description may be an approximation only.

The electric dipole case with its two possible decay paths ($\Delta\ell=\pm
1$) was handled by considering a decomposition of the Kepler motion in such a
way that one could identify the parts that correspond to radiation carrying
away angular momentum of definite sign, but more detailed work is required to
separate the quadrupole transitions into the three contributions
($\Delta\ell=-2,0,2$). Thus, one has to extract information from a multipole
decomposition of the emitted electromagnetic fields in the radiation zone.
Fortunately, Chapter 9 in Jackson’s textbook Jackson (1999b) prepares the
ground for all required calculations. For recent pedagogical work beyond the
dipole approximation we also refer to Ref. Rovenchak and Krynytskyi (2018).

Concerning experimental determination of E2 transitions in hydrogenic systems
we can make the following remarks. The measurement of electric quadrupole
lifetime in cesium for the $\rm 6S-6D$ and $\rm 6S-5D$ transitions using
various methods has been reported in Refs. Zimmermann  _et al._ (1974); Glab
and Nayfeh (1981); Tojo  _et al._ (2005) reaching eventually good precision.
Atomic hydrogen spectroscopy Horbatsch and Hessels (2016) has reached such
high precision that quantum electrodynamics and proton size can be tested, and
this involves spectroscopy which relies on dipole-forbidden transitions (and
often two-photon transitions) Grinin  _et al._ (2020); Fleurbaey  _et al._
(2018). Sum rules for electric quadrupole transitions in two-electron systems
were derived and tested in Ref. Rosenthal  _et al._ (1974). E2 transition
rates for atoms with effectively two active electrons in even-parity states
were treated non-relativistically in Ref. Wilson (1991) using the formalism
provided in Ref. Shortley (1940) and the book Condon and Odabasi (1980), and
reasonable agreement was found with experimental multi-photon ionization rates
for neutral silicon atoms Ioannidou-Phili and Christodoulides (1990). For the
alkalis transition rates were calculated within model potential approaches in
Refs. Çelik  _et al._ (2014); Majumder, Sonjoy  _et al._ (2004), and
expressions for the E2 transition rate for LSJ coupled states can be found
there, as well as in Ref. Condon  _et al._ (1967). There does not, however,
appear to be a direct comparison of experiment versus theory for these
transitions in atomic Rydberg levels, and thus, the present work focuses
purely on the comparison of Schrödinger and Coulomb-Kepler results. An
interesting recent development, however, is the observation of big changes to
normally forbidden transitions in hydrogen when the atoms are placed in a
nano-scale environment with external fields Kim and Kim (2018).

The paper is organized as follows. In Sect. II.1 the semiclassical E2 rates
for $(n_{i},\ell_{i})\to(n_{f},\ell_{f})$ transitions are derived from general
expression given by Jackson Jackson (1999b). In Sect. II.2 we provide the
corresponding calculation for the quantum-mechanical E2 transition rates. In
Sect. III we provide some comparisons for a few cases, indicating strengths
and weaknesses of the semiclassical approach. We present graphical results for
transitions from $(n_{i}=30,\ell_{i}=15)$ to illustrate the decay
possibilities for $\Delta\ell=-2,0,2$, in terms of their branching ratios as a
function of emitted photon frequency, and show how the semiclassical results
compare with the quantum ones. This is followed up by a table which compares
absolute transition rates for a number of cases. The paper ends with
concluding remarks in Sect. IV.

##  II Theory

###  II.1 Semiclassical approach

Here we analyze the electric quadrupole radiation emitted by a classical model
of a hydrogenic atom, i.e. a point particle of charge $-e$ moving along a
Kepler ellipse. The theory of the multipole expansion of radiation fields is
developed in Ref. Jackson (1999b). From Jackson’s equation (9.155), the
radiated power from the $\ell=2$ electric multipole is given by

$P=\int_{0}^{\infty}d\omega\frac{c}{2\epsilon_{0}\omega^{2}}\sum_{m=-2}^{2}\left|a_{E}^{(\omega)}(2,m)\right|^{2}\,.$
(1)

Note that the discussion in Ref. Jackson (1999b) considers only a single
Fourier mode, whereas here we need multiple Fourier modes to describe ellipses
with non-zero eccentricity, hence the integral over $\omega$.

The electric quadrupole coefficients $a_{E}^{(\omega)}(2,m)$ appearing in (1)
may be computed using equations (9.169)-(9.170) from Ref. Jackson (1999b):

$a_{E}^{(\omega)}(2,m)=\frac{\omega^{4}}{15ic^{3}}\sqrt{\frac{3}{2}}\int|\vec{x}|^{2}\,Y_{2,m}^{*}\,\rho_{\omega}\,d^{3}\vec{x}\,,$
(2)

where the relevant spherical harmonics $Y_{\ell,m}$ are given explicitly on
pages 108-109 of Ref. Jackson (1999a), and

$\rho_{\omega}(\vec{x})=\frac{1}{\pi}\int_{-\infty}^{\infty}\rho(\vec{x},t)e^{i\omega
t}dt$ (3)

are the Fourier components of the electric charge density $\rho(\vec{x},t)$.
For a point particle with charge $-e$ and trajectory $\vec{r}(t)$, we have

$\rho(\vec{x},t)=-e\delta^{(3)}(\vec{x}-\vec{r}(t))\,.$ (4)

In principle, in the classical model of an atom, the above equation should
contain an additional term describing the presence of the nucleus, but we do
not write it down, because it does not directly contribute to the radiation as
long as reduced-mass effects are neglected.

In order to evaluate the electric quadrupole coefficients (2), we need an
explicit description of the Kepler trajectory, which is known as a Fourier
series:

$\displaystyle\vec{r}(t)$ $\displaystyle=$
$\displaystyle\Bigl{(}x(t),y(t),0\Bigr{)}\,,$ (5)
$\displaystyle\left[x(t)\right]^{2}$ $\displaystyle=$ $\displaystyle
a^{2}\sum_{k=0}^{\infty}{A}_{k}\cos(k\Omega t)\,,$ (6)
$\displaystyle\left[y(t)\right]^{2}$ $\displaystyle=$ $\displaystyle
a^{2}\sum_{k=0}^{\infty}{B}_{k}\cos(k\Omega t)\,,$ (7) $\displaystyle
x(t)\cdot y(t)$ $\displaystyle=$ $\displaystyle
a^{2}\sum_{k=1}^{\infty}{C}_{k}\sin(k\Omega t)\,,$ (8)

where $a$ and $\Omega$ are the semimajor axis and angular frequency of the
Kepler orbit, respectively, and

$\displaystyle{A}_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{k}\left[J_{k-2}(k\epsilon)-J_{k+2}(k\epsilon)-2\epsilon(J_{k-1}(k\epsilon)-J_{k+1}(k\epsilon))\right]\,,$
(9) $\displaystyle{B}_{k}$ $\displaystyle=$
$\displaystyle\frac{1-\epsilon^{2}}{k}\left[J_{k+2}(k\epsilon)-J_{k-2}(k\epsilon)\right]\,,$
(10) $\displaystyle{C}_{k}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{1-\epsilon^{2}}}{k}\left[J_{k+2}(k\epsilon)+J_{k-2}(k\epsilon)-\epsilon(J_{k+1}(k\epsilon)+J_{k-1}(k\epsilon))\right]\,,$
(11) $\displaystyle{A}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+4\epsilon^{2})\,,$ (12) $\displaystyle{B}_{0}$
$\displaystyle=$ $\displaystyle\frac{1}{2}(1-\epsilon^{2})\,,$ (13)

are the Fourier coefficients Maggiore (2007). The standard notation $J_{k}(z)$
is used for the Bessel function of order $k$ Milton and Stegun (1965). The
constant terms ${A}_{0}$ and ${B}_{0}$ do not enter into any final radiation
formulas, and are only provided for completeness. The quantity $\epsilon$ is
the eccentricity of the Kepler ellipse. By means of semiclassical quantization
Langer (1937), properties of the Kepler orbit may be related to the principal
quantum number $n$, and the orbital quantum number $\ell$:

$a=\frac{n^{2}\hbar}{Z\alpha
m_{e}c}\,,\qquad\Omega=\frac{Z^{2}\alpha^{2}m_{e}c^{2}}{n^{3}\hbar}\,,\qquad\epsilon=\sqrt{1-\frac{(\ell+\tfrac{1}{2})^{2}}{n^{2}}}\,.$
(14)

Here $Z$ is the atomic number (so the nucleus has charge $+Ze$), $m_{e}$ is
the electron mass, and $\alpha=e^{2}/(4\pi\epsilon_{0}\hbar c)\sim 1/137$ is
the fine-structure constant.

With the Fourier series in hand, equation (2) may now be explicitly evaluated,
yielding

$\displaystyle a_{E}^{(\omega)}(2,0)$ $\displaystyle=$
$\displaystyle-\frac{iea^{2}\omega^{4}}{8c^{3}}\sqrt{\frac{2}{15\pi}}\sum_{k=1}^{\infty}({A}_{k}+{B}_{k})\delta(\omega-k\Omega)\,,$
(15) $\displaystyle a_{E}^{(\omega)}(2,\pm 2)$ $\displaystyle=$
$\displaystyle\frac{iea^{2}\omega^{4}}{8\sqrt{5\pi}c^{3}}\sum_{k=1}^{\infty}({A}_{k}-{B}_{k}\pm
2{C}_{k})\delta(\omega-k\Omega)\,,$ (16) $\displaystyle a_{E}^{(\omega)}(2,\pm
1)$ $\displaystyle=$ $\displaystyle 0\,.$ (17)

In the above equations, zero-frequency terms involving $\delta(\omega)$ have
been dropped, since they do not contribute to the radiated power (1). Also,
the vanishing of the $m=\pm 1$ components is a consequence of choosing the
coordinate system such that the Kepler ellipse lies in the $xy$ plane.

Substituting equations (15)-(17) into the power formula (1), we obtain

$\displaystyle P_{m=0}$ $\displaystyle=$
$\displaystyle\frac{e^{2}a^{4}\Omega^{6}}{960\pi\epsilon_{0}c^{5}}\sum_{k=1}^{\infty}k^{6}({A}_{k}+{B}_{k})^{2}\,,$
(18) $\displaystyle P_{m=+2}$ $\displaystyle=$
$\displaystyle\frac{e^{2}a^{4}\Omega^{6}}{640\pi\epsilon_{0}c^{5}}\sum_{k=1}^{\infty}k^{6}({A}_{k}-{B}_{k}+2{C}_{k})^{2}\,,$
(19) $\displaystyle P_{m=-2}$ $\displaystyle=$
$\displaystyle\frac{e^{2}a^{4}\Omega^{6}}{640\pi\epsilon_{0}c^{5}}\sum_{k=1}^{\infty}k^{6}({A}_{k}-{B}_{k}-2{C}_{k})^{2}\,.$
(20)

From the calculation and discussion in sections 9.8 and 9.11 of Ref. Jackson
(1999b) it follows that:

  * •

The $k^{\rm th}$ term in equation (18), divided by $\hbar k\Omega$, is the
rate for $\Delta\ell=0$ transitions.

  * •

The $k^{\rm th}$ term in equation (19), divided by $\hbar k\Omega$, is the
rate for $\Delta\ell=-2$ transitions.

  * •

The $k^{\rm th}$ term in equation (20), divided by $\hbar k\Omega$, is the
rate for $\Delta\ell=+2$ transitions.

Explicitly, and in terms of quantities familiar to atomic physicists, the
classical rates (before harmonic index rescaling) are given by

$R_{k,\Delta\ell}=\frac{Z^{6}\alpha^{7}m_{e}c^{2}}{60n^{7}\hbar}\tilde{R}_{k,\Delta\ell}\,,$
(21)

with

$\displaystyle\tilde{R}_{k,\Delta\ell=0}=\frac{k^{5}}{4}({A}_{k}+{B}_{k})^{2}\,,$
(22) $\displaystyle\tilde{R}_{k,\Delta\ell=\pm
2}=\frac{3k^{5}}{8}({A}_{k}-{B}_{k}\mp 2{C}_{k})^{2}\,.$ (23)

These expressions complement the previously obtained semiclassical rates for
E1 transitions Horbatsch  _et al._ (2005a):

$R_{k,\Delta\ell=\pm
1}=\frac{2Z^{4}\alpha^{5}m_{e}c^{2}}{3n^{5}\hbar}k\left(J_{k}^{\prime}(k\epsilon)\mp\sqrt{\epsilon^{-2}-1}J_{k}(k\epsilon)\right)^{2}\,.$
(24)

The given rates for E1 or E2 transitions are in terms of Fourier order $k$,
i.e., correspond to an equidistant photon energy spectrum $\hbar k\Omega$. In
order to compare with transitions to allowed hydrogenic states with energies
$E_{n^{\prime}}$ in accord with the Balmer formula they can be rescaled to
non-integer order $k_{\Delta n}$ and weighted appropriately, as explained in
Ref. Horbatsch  _et al._ (2005a), and in Section III (cf. eq. (48) and Fig.
1). While the Fourier rates depend only on the properties of the initial state
$(n,\ell)$ whose spontaneous decay is being considered, the rescaling
introduces dependence on the final state, and makes comparison with quantum
transition rates to states $(n^{\prime},\ell^{\prime})$ possible.

###  II.2 Quantum mechanical approach

As mentioned in the introduction most publications reporting electric
quadrupole transition rates quote formulae given in handbooks for the case of
total angular momentum eigenstates which include spin-orbit coupling. These
can be derived using the formalism of the Wigner-Eckart theorem Weissbluth
(1978); Shore (1990). Given that we are interested in quantum-classical
correspondence we would like to work only with orbital angular momentum
eigenstates. The calculation can be carried out directly, and this approach
will be followed here.

Following Ref. Griffiths and Schroeter (2018) (p. 556, eq. 11.130, and
including an apparently omitted factor of $1/4$), the quadrupole transition
rate is given by

$R_{n\ell m\to
n^{\prime}\ell^{\prime}m^{\prime}}=\frac{e^{2}\omega^{5}}{4\pi\epsilon_{0}\hbar
c^{5}}\,\left|\sum_{i,j}{\epsilon_{i}k_{j}\langle
n^{\prime}\ell^{\prime}m^{\prime}|r_{i}r_{j}|n\ell m\rangle}\right|^{2}\,.$
(25)

Here $r_{i}$ are the components of the electron position vector $\vec{r}$,
$\omega$ is the circular photon frequency, $\epsilon_{i}$ are the components
of the (unit) photon polarization vector $\hat{\epsilon}$, and $k_{i}$ are the
components of the (unit) photon propagation direction vector $\hat{k}$. The
latter two vectors must be orthogonal, i.e.

$\hat{k}\cdot\hat{\epsilon}=0\,.$ (26)

The first thing that we need to do, is average the quadrupole rate (25) over
all photon directions $\hat{k}$ and polarizations $\hat{\epsilon}$, while
enforcing the constraint (26). To do this, begin by introducing polar and
azimuthal angles for $\hat{\epsilon}$:

$\epsilon_{x}=\sin\theta\cos\phi\,,\qquad\epsilon_{y}=\sin\theta\sin\phi\,,\qquad\epsilon_{z}=\cos\theta\,.$
(27)

For fixed values of $\theta$ and $\phi$, a geometrical argument involving
rotation matrices can be used to show that the vector $\hat{k}$ is constrained
to lie on a unit circle that may be parametrized according to

$\displaystyle k_{x}$ $\displaystyle=$
$\displaystyle\cos\phi\cos\theta\cos\psi-\sin\phi\sin\psi\,,$ (28)
$\displaystyle k_{y}$ $\displaystyle=$
$\displaystyle\sin\phi\cos\theta\cos\psi+\cos\phi\sin\psi\,,$ (29)
$\displaystyle k_{z}$ $\displaystyle=$ $\displaystyle-\sin\theta\cos\psi\,,$
(30)

where $0\leq\psi<2\pi$. The direction- and polarization-averaged rate is now
given by

$R_{n\ell m\to
n^{\prime}\ell^{\prime}m^{\prime}}=\frac{e^{2}\omega^{5}}{4\pi\epsilon_{0}\hbar
c^{5}}\,\frac{1}{4\pi}\int_{0}^{\pi}d\theta\int_{0}^{2\pi}\sin\theta
d\phi\,\frac{1}{2\pi}\int_{0}^{2\pi}d\psi\,\Bigl{|}\sum_{i,j}{\epsilon_{i}k_{j}Q_{ij}}\Bigr{|}^{2}\,,$
(31)

where we have introduced the notation

$Q_{ij}=Q_{ji}=\langle n^{\prime}\ell^{\prime}m^{\prime}|r_{i}r_{j}|n\ell
m\rangle$ (32)

for the quadrupole matrix elements. Evaluating the angular integrals using a
computer algebra program we obtain

$\displaystyle R_{n\ell m\to n^{\prime}\ell^{\prime}m^{\prime}}$
$\displaystyle=$ $\displaystyle\frac{e^{2}\omega^{5}}{4\pi\epsilon_{0}\hbar
c^{5}}\Biggl{[}\frac{1}{15}\left(|Q_{xx}|^{2}+|Q_{yy}|^{2}+|Q_{zz}|^{2}\right)$
(33) $\displaystyle\phantom{\frac{e^{2}\omega^{5}}{\pi\epsilon_{0}\hbar
c^{5}}\Biggl{[}}+\frac{1}{5}\left(|Q_{xy}|^{2}+|Q_{xz}|^{2}+|Q_{yz}|^{2}\right)$
$\displaystyle\phantom{\frac{e^{2}\omega^{5}}{\pi\epsilon_{0}\hbar
c^{5}}\Biggl{[}}-\frac{1}{30}\left(Q_{xx}Q^{*}_{yy}+Q_{xx}^{*}Q_{yy}\right)$
$\displaystyle\phantom{\frac{e^{2}\omega^{5}}{\pi\epsilon_{0}\hbar
c^{5}}\Biggl{[}}-\frac{1}{30}\left(Q_{xx}Q^{*}_{zz}+Q_{xx}^{*}Q_{zz}\right)$
$\displaystyle\phantom{\frac{e^{2}\omega^{5}}{\pi\epsilon_{0}\hbar
c^{5}}\Biggl{[}}-\frac{1}{30}\left(Q_{yy}Q^{*}_{zz}+Q_{yy}^{*}Q_{zz}\right)\Biggr{]}\,.$

Now we compute the six quadrupole matrix elements. Using explicit formulas for
the spherical harmonics $Y_{\ell,m}$ with $\ell=2$, we obtain

$\displaystyle Q_{xx}$ $\displaystyle=$
$\displaystyle\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}}\left\\{\sqrt{\frac{2\pi}{15}}\left[\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{2}{}^{\ell}_{m}}+\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{-2}{}^{\ell}_{m}}\right]-\frac{1}{3}\sqrt{\frac{4\pi}{5}}\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{0}{}^{\ell}_{m}}+\frac{1}{3}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\right\\}\,,$
(34) $\displaystyle Q_{yy}$ $\displaystyle=$
$\displaystyle\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}}\left\\{-\sqrt{\frac{2\pi}{15}}\left[\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{2}{}^{\ell}_{m}}+\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{-2}{}^{\ell}_{m}}\right]-\frac{1}{3}\sqrt{\frac{4\pi}{5}}\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{0}{}^{\ell}_{m}}+\frac{1}{3}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\right\\}\,,$
(35) $\displaystyle Q_{zz}$ $\displaystyle=$
$\displaystyle\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}}\left\\{\frac{2}{3}\sqrt{\frac{4\pi}{5}}\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{0}{}^{\ell}_{m}}+\frac{1}{3}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\right\\}\,,$
(36) $\displaystyle Q_{xy}$ $\displaystyle=$
$\displaystyle\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}}\left\\{\frac{1}{i}\sqrt{\frac{2\pi}{15}}\left[\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{2}{}^{\ell}_{m}}-\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{-2}{}^{\ell}_{m}}\right]\right\\}\,,$
(37) $\displaystyle Q_{xz}$ $\displaystyle=$
$\displaystyle\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}}\left\\{\frac{1}{2}\sqrt{\frac{8\pi}{15}}\left[\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{-1}{}^{\ell}_{m}}-\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{1}{}^{\ell}_{m}}\right]\right\\}\,,$
(38) $\displaystyle Q_{yz}$ $\displaystyle=$
$\displaystyle\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}}\left\\{\frac{i}{2}\sqrt{\frac{8\pi}{15}}\left[\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{-1}{}^{\ell}_{m}}+\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{2}_{1}{}^{\ell}_{m}}\right]\right\\}\,.$
(39)

Here, we have introduced the notation

$\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}}=\int_{0}^{\infty}drr^{4}R^{*}_{n^{\prime}\ell^{\prime}}(r)R_{n\ell}(r)$
(40)

for the radial integrals ($R_{n\ell}(r)$ are the real-valued hydrogenic radial
wavefunctions), and

$\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{\ell_{1}}_{m_{1}}{}^{\ell_{2}}_{m_{2}}}=\int_{0}^{\pi}d\theta\int_{0}^{2\pi}\sin\theta
d\phi\,Y^{*}_{\ell^{\prime},m^{\prime}}(\theta,\phi)Y_{\ell_{1},m_{1}}(\theta,\phi)Y_{\ell_{2},m_{2}}(\theta,\phi)$
(41)

for the angular integrals. These angular integrals may be expressed in terms
of Clebsch-Gordan coefficients, see, e.g., Ref. Sakurai and Napolitano (2021)
(p. 218, eq. 3.393):

$\tensor{*}{\Lambda}{*^{\ell^{\prime}}_{m^{\prime}}{}^{\ell_{1}}_{m_{1}}{}^{\ell_{2}}_{m_{2}}}=\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)}{4\pi(2\ell^{\prime}+1)}}\tensor{*}{C}{*^{\ell_{1}}_{0}{}^{\ell_{2}}_{0}{}^{\ell^{\prime}}_{0}}\tensor{*}{C}{*^{\ell_{1}}_{m_{1}}{}^{\ell_{2}}_{m_{2}}{}^{\ell^{\prime}}_{m^{\prime}}}\,.$
(42)

To proceed further we observe some properties of the transition rate (33):
summing over the magnetic quantum numbers $m^{\prime}$ for the final state we
find numerically an answer that is independent of the initial-state quantum
number $m$. Usually the initial $m$-value is not specified, and one averages
over $m$, but needs to sum over all final $m^{\prime}$. This leads to the idea
of defining the strength for $(n,\ell)\to(n^{\prime},\ell^{\prime})$ E2
transitions, in analogy to the detailed discussion for electric dipole
transitions in Ref. Weissbluth (1978) (chapter 23):

$S_{n\ell\to n^{\prime}\ell^{\prime}}=\sum_{m,m^{\prime}}{\tilde{R}_{n\ell
m\to n^{\prime}\ell^{\prime}m^{\prime}}}\ ,$ (43)

where $\tilde{R}_{n\ell m\to n^{\prime}\ell^{\prime}m^{\prime}}$ is given by
(33) with the prefactor omitted.

Using angular momentum algebra it is possible to simplify the result such that

$S_{n\ell\to
n^{\prime}\ell^{\prime}}=\frac{(2\ell+1)(2\ell^{\prime}+1)}{15}\begin{pmatrix}\ell^{\prime}&2&\ell\\\
0&0&0\end{pmatrix}^{2}(\tensor{*}{I}{*^{n^{\prime}}_{\ell^{\prime}}{}^{n}_{\ell}})^{2}\
,$ (44)

where the term in round brackets is a Wigner 3j symbol and the squares of the
radial integrals (40) need to be computed. In order to compare with the
semiclassical transition rates one needs to include the pre-factor given in
(33), and divide by the multiplicity of the initial state in order to account
for statistical averaging.

In order to compare with literature values obtained within relativistic
quantum mechanics Jitrik and Bunge (2004) one has to convert the transition
strength to states with defined total angular momentum in accord with the LSJ
coupling scheme. This can be accomplished as outlined by Weissbluth Weissbluth
(1978) for the dipole case (eq. 23.1-25), but using the Wigner-Eckart theorem
(eq. 6.3-39) for the tensor rank $q=2$. This results in the transition
strength expression

$S_{n\ell j\to
n^{\prime}\ell^{\prime}j^{\prime}}=(2j+1)(2j^{\prime}+1)\begin{Bmatrix}j^{\prime}&2&j\\\
l&\tfrac{1}{2}&l^{\prime}\end{Bmatrix}^{2}S_{n\ell\to
n^{\prime}\ell^{\prime}}\ ,$ (45)

where generally $j=\ell\pm\tfrac{1}{2}$, and the term in curly brackets is a
Wigner 6j symbol.

The final expressions then become the following: for comparison with the
semiclassical rates we have

$R_{n\ell\to
n^{\prime}\ell^{\prime}}=\frac{e^{2}\omega^{5}}{4\pi\epsilon_{0}\hbar
c^{5}}\frac{1}{2\ell+1}S_{n\ell\to n^{\prime}\ell^{\prime}}\ ,$ (46)

while for comparison with the relativistic calculations we use

$R_{n\ell j\to
n^{\prime}\ell^{\prime}j^{\prime}}=\frac{e^{2}\omega^{5}}{4\pi\epsilon_{0}\hbar
c^{5}}\frac{1}{2j+1}S_{n\ell j\to n^{\prime}\ell^{\prime}j^{\prime}}\ .$ (47)

##  III Results

We consider results for hydrogen only, i.e., $Z=1$. We have checked that the
quantum LSJ coupling results based on the Schrödinger equation treatment,
i.e., eq. (45) agrees to several digits with the Dirac equation based results
from Ref. Jitrik and Bunge (2004), which are available in tabular form for
many states Jitrik and Bunge (2003). For higher values of $Z$ relativistic
effects due to the energy differences between initial and final states, as
well as the character of the radial functions will become important. Thus, our
comparison between semiclassical and quantum results will be based solely on
the $(n,\ell)\to(n^{\prime},\ell^{\prime})$ transitions. Electric quadrupole
E2 transitions are stated in the literature to be five orders of magnitude
weaker than E1 transitions, which can be estimated from the prefactors given,
e.g., in eqs. (21) vs. (24). Therefore, the most important comparisons will be
for those cases where a dipole transition cascade will not connect two levels.
It would be an interesting exercise in its own right to check how a sequence
of two E1 transitions (e.g., when $n-n^{\prime}\geq 2$) would compete against
a direct E2 transition.

We consider first a high-$n$, intermediate-$\ell$ initial state to allow for
pathways involving $\Delta\ell=-2,0,2$ transitions, and begin with a
comparison of the overall transition rates to lower principal quantum number
$n^{\prime}=n-\Delta n$ ($\Delta n$ is positive). In the classical limit
electromagnetic radiation at Fourier order $k$ is approximately treated as
corresponding to a transition to a lower state according to $\Delta n=k$. The
Balmer spectrum is not equidistant, however, and therefore the re-scaling
approach of Ref. Horbatsch  _et al._ (2005a) was introduced to re-scale the
Fourier coefficients to non-integer order

$k_{\Delta n}=\frac{n}{2}\left[(1-\Delta n/n)^{-2}-1\right]$ (48)

which also implies a weighting of this effective harmonic transition rate by a
factor of $(1-\Delta n/n)^{-3}$, (cf. Eqs. (22,23) in Ref. Horbatsch  _et al._
(2005a)).

$\begin{array}[]{cc}\leavevmode\resizebox{216.81pt}{}{\includegraphics{Fig2a.pdf}}&\leavevmode\resizebox{216.81pt}{}{\includegraphics{Fig2b.pdf}}\\\
\leavevmode\resizebox{216.81pt}{}{\includegraphics{Fig2c.pdf}}&\leavevmode\resizebox{216.81pt}{}{\includegraphics{Fig2d.pdf}}\end{array}$

Figure 1: Branching ratios for E2 quadrupole transitions from the state
$(n=30,\ell=15)$ to states $n^{\prime}=n-\Delta n$ as a function of emitted
radiation frequency $f$. The red dots represent quantum-mechanical ratios, the
green dots show the classical ratios, and the blue dots show the re-scaled
classical branching ratios according to eq. (48). The top left panel shows the
quantum and re-scaled classical branching ratios as a function of photon
frequency only, i.e., irrespective of $\Delta\ell$, while the other three
panels show the $(\Delta n,\Delta\ell)$-dependent ratios separately for
$\Delta\ell=-2$, $\Delta\ell=0$, and $\Delta\ell=+2$ transitions.

The top left panel of Fig. 1 shows the $(\Delta n)$-dependent branching
fractions in percent irrespective of whether angular momentum is decreasing,
remains the same, or is increasing for the transitions. The points correspond
to transitions with $\Delta n=1,2,...$ from left to right, and are only shown
for the quantum and re-scaled classical results. A peak occurs at $\Delta
n=13$, or about $6-8\ \rm THz$ in frequency. At this point the rates exceed
the low-$\Delta n$ results by almost two orders of magnitude. The minimum at
$\Delta n=3$ which is obtained both in the quantum and in the re-scaled
classical calculation can be explained when looking at the $\Delta\ell=-2$
transitions alone.

Considering now the branching fractions associated with angular momentum
decreasing transitions (top right panel) we observe a number of things. First
of all, the green points show the classical Fourier spectrum of frequencies,
which is equidistant (harmonic). At small $k\approx\Delta n$ it follows the
quantum results quite well, also exhibiting a minimum. There is a near-
cancellation at this Fourier order, i.e., at $k=3$. This is happening in the
combination of Fourier coefficients in eq. (23), caused by the fact that there
is a sign change in the combination $({A}_{k}-{B}_{k}+2{C}_{k})$ for some
(non-integer) value of $k=k_{c}$. In the quantum calculation it is the
expectation value $\langle n,\ell|r^{2}|n^{\prime},\ell-2\rangle$ that changes
sign between two subsequent values of $n^{\prime}=n-\Delta n$. The result is
that for given $(n,\ell)$ at particular $\Delta n$ a minimum occurs. For
higher values of $n$ this minimum is shifted to higher $\Delta n$, or
classical Fourier order $k$.

With increasing $\Delta n$ the quantum transition energies spread apart
according to the Balmer formula, thus causing a distinctive pattern in the
frequency spectrum. The transition rates increase with $\Delta n$ due to the
phase space factor, i.e., higher photon energies are favored. Eventually,
however, there is a turnover, because the radial matrix elements become
unfavorable. The classical Fourier series (green points) also displays a
maximum (at around 5 THz). Many equally-spaced Fourier orders contribute, and
thus the branching fractions are smaller. They need to be added up by a
binning procedure to take into account the allowed final states. This is
accomplished in a continuous way by the re-scaling given in eq. (48). The re-
scaling makes use of non-integer index $k_{\Delta n}$, and a weight is applied
to the result. The blue points show that this is effective in the sense that a
similar distribution is obtained as for the quantum case.

While the quantum calculation depends on the initial- and final-state
wavefunctions, the re-scaled semiclassical calculation makes use of the
initial-state Fourier series, but then applies the Balmer formula to sum up
the contributions from multiple Fourier orders that correspond to a transition
to final $n^{\prime}=n-\Delta n$. At moderate values of $\Delta n$ there is an
overestimation of the transition rates (blue points), but at high $\Delta n$ a
very similar maximum as for the quantum results can be observed, and then a
cross-over, where the highest-allowed frequency contribution is overestimated.

For the $\Delta\ell=0$ transitions shown in the bottom left panel no
particularly interesting feature arises, the transition rates rise to form a
maximum at $\Delta n=3-4$, and then fall off. The classical Fourier rates
(green points) when binned effectively to correspond to quantized final-state
energies using eq. (48) underestimate the quantum rates at the maximum, and
then overestimate them strongly at the highest allowed frequencies. A similar
behavior is found for the $\Delta\ell=+2$ transitions shown in the bottom
right panel.

The decays from a high-$n$, intermediate-$\ell$ initial state discussed above
allowed us to contrast the classical Fourier radiation spectrum from the
elliptic orbit that represents the initial state with the available space of
quantum levels. Through an effective continuous binning process (eq. (48)) we
demonstrated how quantum transitions can be made compatible with a
semiclassical picture. For the case shown it works reasonably well for the
total $n\to n^{\prime}$ transition rates at not too large values of $\Delta
n=n-n^{\prime}$. It is based mostly on a reasonable treatment of
$\Delta\ell=-2$ transitions, and a moderately successful representation of the
case of $\Delta\ell=0$, while the overall less important case of
$\Delta\ell=+2$ was seen to be described more poorly in the semiclassical
analysis.

In order to show results for cases where the absolute transition rates are
larger we present in table 1 some numerical comparisons for the cases of
$\Delta\ell=-2$, $\Delta\ell=0$, and $\Delta\ell=+2$. A classic textbook case
is of course the $3d\to 1s$ transition (this is a case where a sequence of E1
transitions, namely $3d\to 2p\to 1s$ will dominate the overall transition
strength). The semiclassical calculation based on re-scaling does not fare too
well, being too small by a factor of about 3. It is interesting to observe,
however, that the agreement becomes quite good as the initial-state principal
quantum number $n$ is increased. Moving to higher initial values of $\ell$
shows that the semiclassical treatment yields respectable results.

For the $\Delta\ell=0$ transitions we observe that for $\ell=\ell^{\prime}=1$,
i.e., p-state to p-state transitions, factor-of-two level agreement is
attained, and that the situation improves considerably as one moves to
d-states, and higher angular momentum, making the semiclassical result eqs.
(21-23) a practical contender for evaluating transition strengths, which could
be useful for estimates of absorption of electromagnetic radiation.

$(n,\ell)\to(n^{\prime},\ell^{\prime})$ | QM | SCL | $(n,\ell)\to(n^{\prime},\ell^{\prime})$ | QM | SCL
---|---|---|---|---|---
$3d\to 1s$ | 594. | 201. | $4f\to 2p$ | 61.8 | 70.1
$4d\to 1s$ | 327. | 256. | $4f\to 3p$ | 5.79 | 6.42
$5d\to 1s$ | 185. | 167. | $5g\to 3d$ | 11.6 | 13.4
$6d\to 1s$ | 112. | 107. | $5g\to 4d$ | 1.00 | 0.85
$3p\to 2p$ | 23.9 | 14.7 | $4d\to 3d$ | 1.19 | 1.06
$4p\to 2p$ | 10.3 | 6.29 | $5d\to 3d$ | 0.574 | 0.521
$5p\to 2p$ | 5.27 | 3.22 | $5f\to 4f$ | 0.153 | 0.150
$6p\to 2p$ | 3.05 | 1.86 | $6f\to 5f$ | 0.053 | 0.050
$5p\to 4f$ | 0.047 | 0.158 | $6d\to 5g$ | 5.2(-3) | 2.1(-2)
$6p\to 4f$ | 0.032 | 0.124 | $7f\to 6h$ | 9.0(-4) | 4.0(-3)
$7p\to 4f$ | 0.021 | 0.089 | $8g\to 7i$ | 2.0(-4) | 9.5(-4)
$8p\to 4f$ | 0.014 | 0.064 | $9h\to 8k$ | 5.5(-5) | 2.7(-4)
Table 1: Some quantum versus semiclassical E2 transition rates in decays per
second. The column labeled QM shows quantum rates obtained from eq. (46),
while the column labeled SCL shows the semiclassical rates according to eqs.
(21-23). The top third of the table shows $\Delta\ell=-2$ transition rates,
while the middle is for $\Delta\ell=0$, and the bottom third contains results
for $\Delta\ell=+2$. The numbers in the bottom right part of the table are in
a notation where $x.y(-i)$ represents $x.y\times 10^{-i}$.

The case of $\Delta\ell=+2$, however, is somewhat disappointing, in that one
finds at best order-of-magnitude agreement. When comparing the semiclassical
and quantum calculations we note that there are two aspects to consider. On
the one hand there is the radial matrix element of $r^{2}$, which enters
equally for all three cases, $\Delta\ell=0,\pm 2$. The Fourier series of the
elliptic orbit somehow encodes the right physics into the coefficients, such
that a correspondence with the squared quantum matrix element becomes
possible. How do the calculations separate the transition strength into the
parts that correspond to $\Delta\ell=0,\pm 2$ transtions respectively? This is
accomplished in the semiclassical case by forming in eqs. (21-23) particular
linear combinations of the coefficients $A_{k},B_{k},C_{k}$. In the quantum
calculation, on the other hand, it is angular momentum coupling that dictates
how the transition strength splits up. The common aspect of these calculations
is, of course, that the total angular momentum is conserved when taking into
account the emitted electromagnetic radiation. While we have no detailed
understanding as to why the orbital angular momentum decreasing transitions
are described more favorably in the semiclassical approach, we note the
similarity with the E1 dipole case Horbatsch  _et al._ (2005a), where the
semiclassical $\Delta\ell=-1$ branching fractions agreed much better with the
quantum ones than those for $\Delta\ell=+1$.

##  IV Conclusions

This work expands on previous work for electric dipole transitions Horbatsch
_et al._ (2005a) in hydrogenic atoms demonstrating that quantum-classical
correspondence can be observed in a more complicated process connected with
emission and absorption of electromagnetic waves with higher angular momentum.
The semiclassical calculation involves a combination of mechanics (description
of the time evolution of the Kepler orbit by a Fourier series), and its
coupling to electromagnetic radiation. The main basis for evaluating the
coupling of the orbital motion and the properties of the emitted radiation is
provided in the graduate text Ref. Jackson (1999a). The quantum calculation
was done mostly without sophistication by brute force and computer algebra,
but final expressions required the use of the Wigner-Eckart theorem in general
form, as presented in the graduate text Ref. Weissbluth (1978). While we
provide the relevant expression for LSJ coupling we leave it as an exercise
for students to perform the comparison with the relativistic rates Jitrik and
Bunge (2003), which are available for higher $Z$, and based on eq. (47) in the
present work. Observation of what happens as $Z$ is increased would form a
useful exercise for graduate students. Finally, we may consider the present
semiclassical calculation a useful warm-up exercise before attempting a study
of gravitational radiation in a post-Newtonian framework.

###### Acknowledgements.

Financial support from the Natural Sciences and Engineering Research Council
of Canada (NSERC) is acknowledged. Marko Horbatsch acknowledges contributions
made by Nick Hamer during the early stages of this project, and his support
through an NSERC undergraduate research award.

## References

  * Landau and Lifshitz (1975) L. D. Landau and E. M. Lifshitz,  _The Classical Theory of Fields, 4th ed._ (Butterworth-Heinemann, Oxford, U.K., 1975).
  * Jackson (1999a) J. D. Jackson,  _Classical Electrodynamics, 3rd ed._ (Wiley, New York, 1999).
  * Corney (2006) A. Corney,  _Atomic and Laser Spectroscopy_ (Oxford University Press, Oxford, U.K., 2006).
  * Horbatsch  _et al._ (2005a) M. W. Horbatsch, E. A. Hessels, and M. Horbatsch, Phys. Rev. A 71, 020501 (2005a).
  * Horbatsch  _et al._ (2005b) M. W. Horbatsch, E. A. Hessels, and M. Horbatsch, Phys. Rev. A 72, 033405 (2005b).
  * Delos  _et al._ (1983) J. B. Delos, S. K. Knudson, and D. W. Noid, Phys. Rev. A 28, 7 (1983).
  * Horbatsch and Horbatsch (2021) M. Horbatsch and M. Horbatsch, Journal of Physics B: Atomic, Molecular and Optical Physics 54, 065004 (2021).
  * Jitrik and Bunge (2004) O. Jitrik and C. F. Bunge, Physica Scripta 69, 196 (2004).
  * Bethe and Salpeter (1957) H. A. Bethe and E. E. Salpeter,  _Quantum Mechanics of One- and Two-Electron Atoms_ (Springer-Verlag, Berlin-Heidelberg, 1957).
  * Weissbluth (1978) M. Weissbluth,  _Atoms and Molecules_ (Academic Press, New York, 1978).
  * Griffiths and Schroeter (2018) D. J. Griffiths and D. F. Schroeter,  _Introduction to Quantum Mechanics, 3rd ed._ (Cambridge University Press, Cambridge, U.K., 2018).
  * Peters and Mathews (1963) P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963).
  * Landau and Lifshitz (1981) L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics: Non-Relativistic Theory, 3rd ed._ (Butterworth-Heinemann, Oxford, U.K., 1981).
  * Jackson (1999b) J. D. Jackson,  _Classical Electrodynamics, 3rd ed._ (Wiley, New York, 1999) Chap. 9.
  * Rovenchak and Krynytskyi (2018) A. Rovenchak and Y. Krynytskyi, American Journal of Physics 86, 727 (2018), https://doi.org/10.1119/1.5052427 .
  * Zimmermann  _et al._ (1974) P. Zimmermann, T. W. Ducas, M. G. Littman, and D. Kleppner, Optics Communications 12, 198 (1974).
  * Glab and Nayfeh (1981) W. Glab and M. H. Nayfeh, Optics Communications 38, 262 (1981).
  * Tojo  _et al._ (2005) S. Tojo, T. Fujimoto, and M. Hasuo, Phys. Rev. A 71, 012507 (2005).
  * Horbatsch and Hessels (2016) M. Horbatsch and E. A. Hessels, Phys. Rev. A 93, 022513 (2016).
  * Grinin  _et al._ (2020) A. Grinin, A. Matveev, D. C. Yost, L. Maisenbacher, V. Wirthl, R. Pohl, T. W. Hänsch, and T. Udem, Science 370, 1061 (2020), https://science.sciencemag.org/content/370/6520/1061.full.pdf .
  * Fleurbaey  _et al._ (2018) H. Fleurbaey, S. Galtier, S. Thomas, M. Bonnaud, L. Julien, F. m. c. Biraben, F. m. c. Nez, M. Abgrall, and J. Guéna, Phys. Rev. Lett. 120, 183001 (2018).
  * Rosenthal  _et al._ (1974) D. Rosenthal, R. P. McEachran, and M. Cohen, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 337, 365 (1974), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1974.0055 .
  * Wilson (1991) M. Wilson, Zeitschrift für Physik D Atoms, Molecules and Clusters 21, 7 (1991).
  * Shortley (1940) G. H. Shortley, Phys. Rev. 57, 225 (1940).
  * Condon and Odabasi (1980) E. U. Condon and H. Odabasi,  _Atomic Structure_ (Cambridge University Press, Cambridge, U.K., 1980).
  * Ioannidou-Phili and Christodoulides (1990) A. Ioannidou-Phili and A. A. Christodoulides, Zeitschrift für Physik D Atoms, Molecules and Clusters 16, 271 (1990).
  * Çelik  _et al._ (2014) G. Çelik, Y. Gökçe, and M. Yıldız, Atomic Data and Nuclear Data Tables 100, 792 (2014).
  * Majumder, Sonjoy  _et al._ (2004) Majumder, Sonjoy, Gopakumar, G., Chaudhuri, R. K., Das, B. P., Merlitz, H., Mahapatra, U. S., and Mukherjee, D., Eur. Phys. J. D 28, 3 (2004).
  * Condon  _et al._ (1967) E. U. Condon, E. W. Condon, and H. Odishaw,  _Handbook of Physics_ (McGraw-Hill, New York, 1967) Chap. 9, p. 177.
  * Kim and Kim (2018) H. Y. Kim and D. S. Kim, Nanophotonics 7, 229 (2018).
  * Maggiore (2007) M. Maggiore,  _Gravitational Waves, Volume 1: Theory and Experiments_ (Oxford University Press, Oxford, U.K., 2007).
  * Milton and Stegun (1965) A. Milton and I. Stegun,  _Handbook of Mathematical Functions_ (Dover Publications, New York, 1965).
  * Langer (1937) R. E. Langer, Phys. Rev. 51, 669 (1937).
  * Shore (1990) B. W. Shore,  _The Theory of Coherent Atomic Excitation_ (Wiley, New York, 1990) Chap. 20.
  * Sakurai and Napolitano (2021) J. J. Sakurai and J. Napolitano,  _Modern Quantum Mechanics, 3rd ed._ (Cambridge University Press, Cambridge, U.K., 2021).
  * Jitrik and Bunge (2003) O. Jitrik and C. F. Bunge, “Atomic transition probabilities for hydrogen-like atoms,” (2003), https://www.fisica.unam.mx/research/tables/spectra/1el/, Last accessed on 2021-07-13.


# User-Guided Free-Form Asset Modelling

Daniel Beale University of Bath. Email<EMAIL_ADDRESS>Supported by the
EPSRC grant EP/K02339X/1 - ”Acquiring Complete and Editable Outdoor Models
from Video and Images”

In this paper a new system for piecewise primitive surface recovery on point
clouds is presented, which allows a novice user to sketch areas of interest in
order to guide the fitting process. The algorithm is demonstrated against a
benchmark technique for autonomous surface fitting, and, contrasted against
existing literature in user guided surface recovery, with empirical evidence.
It is concluded that the system is an improvement to the current documented
literature for its visual quality when modelling objects which are composed of
piecewise primitive shapes, and, in its ability to fill large holes on
occluded surfaces using free-form input.

##  1 Introduction

Creating models of existing buildings and man made structures is a problem for
a range of different industries such as visual effects, gaming and
architecture, and is one of the most important challenges for computer
graphics. It is possible for highly skilled artists to reconstruct an existing
building if reference photographs and plans are available. This process can be
expensive, motivating the search for alternatives that build models directly
from image or video data with little or no demand placed on the user.

The need for 3D models from the real world has lead to a wide variety of
research and systems. Previous work in this field is not far removed from the
paradigm of traditional geometric modelling packages such as
Maya111http://www.autodesk.co.uk/products/maya/overview and
Modo222https://www.thefoundry.co.uk/products/modo/; the user, at a
conventional workstation, is able to perform standard operations such as
selecting specific points and drawing precise curves. While this workflow
leads to models which are visually pleasing from the point of view of an
artist or director, they do not easily model the geometry of a real existing
scene, leaving room for techniques for speeding up the process through
automation.

Standard industrial reconstruction pipelines recreate surfaces from
photographs by first building a point cloud, which can be obtained from a
number of different technologies such as photographs or range scanners. A
common problem among all of these technologies is the undesired production of
holes in regions which are occluded, for example, or exhibit a high degree of
specular reflection or, indeed, no reflection at all. Autonomous methods for
filling holes rely on prediction and prior knowledge of the expected surface
shape. User guided techniques are able to draw on human prior experience for
identifying surfaces, but also require a way to interact with the system in a
way which is not time consuming or difficult.

Rather than motivating the problem through the need for precise models,
attention is drawn to the necessity of accurate and fast modelling using
current point cloud technology. The demands of modern film production, and
increasingly television and gaming, requires pre-
visualisation333https://www.smpte.org/sections/hollywood/previsualization-
current-filmmaking444http://digitalfilmfarmworkshops.com/pre-visualize-
directing-video-film/ ; this is the task of rapidly prototyping and
visualising the individual shots of a film before actual production. While
pre-visualization has clear demands from a technical point of view; including
the need for logistical planning, safety, rigging and setup, and managing
costs; there are also many artistic advantages, for example, ensuring that you
get the shots you want during production. With the rise of virtual production,
the turn-around times between pre-visualisation and production are getting
smaller. In the area of 3D modelling, these time pressures give rise to a need
for rapid prototyping of geometry that can be used to visualise and plan shots
before subsequent detailed modelling for final post-production and composting.

The need for methods which allow a user to quickly build surface models from
point clouds are identified, while maintaining robustness to large holes and
noise. The new system introduced in this paper allows a user to guide fitting
through a collection of free form curves. The workflow first requires the user
to highlight required points on the cloud, and then fit a collection of
maximum-likelihood primitive surfaces over a number of user interactions. The
final surfaces are closed using a well-known and contemporary iso-surface
algorithm. It is concluded that the method presented in this paper produces
high quality surfaces, and is an improvement beyond contemporary methods when
compared on a dataset of man-made objects.

##  2 Background

Reconstructing point clouds and surfaces has motivated the fields of computer
vision and graphics for many years, and remains an area of active interest.
Methods are developed which remain accurate in the presence of noise and
occlusion, maintain a low algorithmic complexity and also allow creative but
precise user input.

Since surface reconstruction has been successfully used in industry for some
time now, the common types of software pipeline are detailed with allusion to
creative contributions in the area. A large number of algorithms are
documented in order to solve specific problems that a general reconstruction
algorithm cannot. To begin with, major, fully autonomous systems for surface
recovery from images are detailed, and also, alternative hardware such as
range scanners. Methods for improving the process through user guidance are
introduced towards the end of the section.

There is a large body of research in to surface reconstruction from point
clouds, with algorithms that are used in modelling software such as The
Foundry’s Nuke555https://www.thefoundry.co.uk/products/nuke/. The area has a
well formed and standard pipeline, which is useful for the range of different
surfaces that it can deal with, and also for its ability to process point
clouds which are generated from alternative hardware. The pipeline consists of
the following processes,

  1. 1.

Point cloud reconstruction. A dense collection of points is reconstructed from
images using available software e.g.[Changchang(2011a), Changchang(2011b),
Furukawa(2010), Kolmogorov(2004)]. Alternatively one can use a different
sensor such as a range finder, radar or sonar with varying reconstruction
qualities.

  2. 2.

Normal estimation. A consistent set of normals are estimated for each point in
the cloud. This indicates both where the surface is and also which side of the
surface a point is on [Mitra and Nguyen(2003), Hoppe et al.(1992)Hoppe,
DeRose, Duchamp, McDonald, and Stuetzle].

  3. 3.

Surface fitting. Finally, an implicit surface; such as a Poisson surface
[Kazhdan et al.(2006)Kazhdan, Bolitho, and Hoppe, Kazhdan and Hoppe(2013),
Bolitho et al.(2009)Bolitho, Kazhdan, Burns, and Hoppe], a Gaussian Process
[Williams and Fitzgibbon(2007), Dragiev et al.(2011)Dragiev, Toussaint, and
Gienger] or Delaunay triangulation [Gopi et al.(2000)Gopi, Krishnan, and
Silva]; or a parametric surface such as a B-Spline [Krishnamurthy and
Levoy(1996), Wang et al.(2006)Wang, Pottmann, and Liu] is fit to the data.

A main alternative to the method described above is to reconstruct a volume or
surface directly from a collection of images, without the intermediate point
cloud reconstruction. These methods produce excellent results when there are a
large number of cameras viewing the object (e.g. [Furukawa et
al.(2009)Furukawa, Curless, Seitz, and Szeliski, Campbell et
al.(2010)Campbell, Vogiatzis, Hernández, and Cipolla, Starck et
al.(2006)Starck, Miller, and Hilton]), but will not work when there are only a
few images. They are a good candidate for reconstructing models of objects
which are small enough to be viewed from a large number of angles, but are
generally unsuitable for medium to large scale objects such as buildings and
landscapes.

The primary concern of this paper is with man made objects, on which the
assumption of ‘piecewise primitive’ surface structure is reasonable. It is
beneficial to assume that the surface is piecewise primitive if there are
large holes in the point cloud or it is extremely noisy, since the assumptions
places a large enough constraint that it is possible to predict with a high
degree of accuracy in unknown areas. There are methods which are able to fit
primitives such as toroids, quadrics, or super-quadrics without user
interaction [Debevec et al.(1996)Debevec, Taylor, and Malik, Musialski et
al.(2013)Musialski, Wonka, Aliaga, Wimmer, Gool, and Purgathofer, Bo et
al.(2012)Bo, Ling, and Wang]. For example, [Monszpart et al.(2015)Monszpart,
Mellado, Brostow, and Mitra] fit a collection of planes to the surface using
an energy minimisation framework; [Schnabel et al.(2007)Schnabel, Wahl, and
Klein] use probabilistic RanSaC to fit primitives including toroids and
cylinders; [Li et al.(2011)Li, Wu, Chrysathou, Sharf, Cohen-Or, and Mitra]
optimise using a formulation which minimises global and pairwise energy; while
[Labatut et al.(2009)Labatut, Pons, and Keriven] minimise energy using a
binary space partitioning tree. For man made objects this approach is
reasonable, and well studied, but complicated natural phenomena, such as tree
bark or moss, break the assumptions and so produce poor quality results.

Incorporating the use of a small amount of human interaction can be a large
benefit to the process (e.g. [van den Hengel et al.(2007)van den Hengel, Dick,
Thormählen, Ward, and Torr, Portelli et al.(2010)Portelli, Ganovelli, Tarini,
Cignoni, Dellepiane, and Scopigno]). One can provide extra points and normals
to surface by selecting a plane which is orthogonal to a point on the medial
axis and sketching on it [Yin et al.(2014)Yin, Huang, Zhang, Gong, Cohen-Or,
and Chen], before sweeping along the skeleton to fit the surface. Methods
which fit a collection of planes to the surface are useful since they easily
snap together along the edges [Sinha et al.(2008)Sinha, Steedly, Szeliski,
Agrawala, and Pollefeys, Nan et al.(2010)Nan, Sharf, Zhang, Cohen-Or, and
Chen]. Using quadric fitting with interaction has been studied in [Andrews and
Séquin(2014)]. The authors allow a user to identify a collection of points on
a mesh from a single stroke and fit a quadric, displaying triangles which are
near the original mesh. [Cao et al.(2014)Cao, Ju, Fu, and Hu] consider objects
which are made by extruding 2D profiles, such as furniture and buildings. The
work described in [Habbecke and Kobbelt(2012)] use software to apply
constraints to designs that are made against pictures, but do not use computer
vision in the fitting process.

Alternative to the multiple image / 3D-reconstruction pipeline, shape and
volume can be quickly recovered from a single image, without any knowledge of
the scene geometry [Lau et al.(2010)Lau, Saul, Mitani, and Igarashi, Chen et
al.(2013)Chen, Zhu, Shamir, Hu, and Cohen-Or, Jiang et al.(2009)Jiang, Tan,
and Cheong]. In 3-Sweep [Chen et al.(2013)Chen, Zhu, Shamir, Hu, and Cohen-Or]
a user can quickly outline contours on objects and recover depth using high
level knowledge about object ordering, the resulting models can then be
recomposed into different photographs. Methods in this area can quickly
recover plausible looking geometry. Without the addition of strong, and
potentially unavailable, priors, they do not have the same metric accuracy as
a point cloud reconstruction.

It is worth mentioning that editable assets can also be generated from scratch
using interaction alone, without any computer vision (e.g. [Chen et
al.(2008)Chen, Kang, Xu, Dorsey, and Shum, Lipp et al.(2014)Lipp, Wonka, and
Müller, Xu et al.(2014)Xu, Chang, Sheffer, Bousseau, Mccrae, and Singh,
Iarussi et al.(2015)Iarussi, Bommes, and Bousseau]). These family of
approaches are suited to modelling novel buildings from an architectural
standpoint, and belong to a large class of CAD modelling software. The prior
use of point-based reconstruction provides a cue to speed up the extraction of
surfaces on existing buildings. It should be viewed as a potential aid for
current CAD modelling tools, removing the need for intensive user interaction
where a 3D model already exists.

| Fit type | Input
---|---|---
Sinha et al. | Planes | Polygon tool
VideoTrace | B-Splines | Spline tool
Morphit | B-Splines | Free-form
This system | Quadric forms | Free-form
Table 1: A comparison of some relevant, and current, interactive techniques
for surface reconstruction. Figure 1: A schematic diagram of the system.

The system presented in this paper follows the standard modelling pipeline
described at the beginning of this section. A point cloud is reconstructed
first, before any user interaction, to extract points for fitting geometric
primitives. This places the work in to a mainstream research area, allowing it
to be compared to fully autonomous point cloud techniques, and also
interactive multiple view surface reconstruction algorithms.

The literature arranged in Table 1 is a representative sample of the current
state of the art in interactive surface recovery. It compares the system
presented in this paper to the most recent and most similar alternative
systems in the area. It can be seen that while [Sinha et al.(2008)Sinha,
Steedly, Szeliski, Agrawala, and Pollefeys] and VideoTrace [van den Hengel et
al.(2007)van den Hengel, Dick, Thormählen, Ward, and Torr] are able to fit
primitive parametric models to the data, they cannot deal with free-form
input. The method in Morphit [Yin et al.(2014)Yin, Huang, Zhang, Gong, Cohen-
Or, and Chen] allows a user to improve implicit surface fitting using free-
form input, which works well on complicated surfaces but cannot predict when
there are large holes in the surface. While there does not appear to be any
previous work that uses rough, or loose, annotations for primitive geometry
recovery, there are other applications that have similar interfaces for
related editing tasks such as colorization [Levin et al.(2004)Levin,
Lischinski, and Weiss] and segmentation [Subr et al.(2013)Subr, Paris, Soler,
and Kautz].

The limitations of contemporary work is given a more detailed review in the
results section, with some specific examples and a more in depth discussion.
It is concluded that the method presented in this paper is a contribution to
the area of free-form modelling with geometric primitives.

##  3 System Overview

|  Point selection |  |  Surface fitting
---|---|---|---
Quadric |  |  $\rightarrow$ |
LVM |  |  $\rightarrow$ |
Figure 2: An example of some sketch input to the system and the corresponding output from the primitive fitting algorithms. The first row contains an example of fittinAn example of someg an ellipsoid to a column of a building, and the second row is an example of a perpendicular extrusion using a pair of latent variable models. Primitive geometry | Poisson reconstruction
---|---
|
Figure 3: After iteratively selecting points in the original cloud, system
produces a collection of primitives which form the surface of the required
object. The left image shows an example of the extracted geometry, and the
right image shows the mesh after it has been closed using Poisson
reconstruction.

The system assumes the existence of a point cloud with a collection of
calibrated images. The user is presented with two panes, one which displays an
interface with the image and allows a frame number to be selected, and the
other containing the point cloud and tools to view the cloud from an arbitrary
pose and scale. Sketches can be drawn on the image pane with colour selection
representing curves which correspond across each of the views.

A schematic diagram for the system is presented in Figure 1. The point clouds
are generated through image based reconstruction, using the well known and
freely available software VisualSFM [Changchang(2011a), Furukawa(2010)]. The
rest of this paper describes the interactive modelling, which makes up the
principle contribution. The meshes are finalised with a Poisson reconstruction
[Kazhdan et al.(2006)Kazhdan, Bolitho, and Hoppe] and then textured.

The system is a free form, colour coded, tool for highlighting corresponding
trajectories. The trajectories enable a sub-selection of the point cloud to be
extracted for more advanced fitting. The types of surface that can be fit are
first classified into ‘whole surface fitting’ and ‘combining curves’. The
former reduces the point cloud to a single point sub-cloud, from which a
surface derived from a quadric is fitted, either an Ellipsoid or Cylinder. The
latter takes two point sub-clouds representing two curves in three dimensions
and combines them through perpendicular extrusion or interpolation. The
categories are labelled below,

  1. 1.

Whole surface fitting. A single point cloud is derived from the user sketches,
from which the algorithm fits, either,

    1. (a)

an Ellipsoid,

    2. (b)

or a Cylinder.

The method for fitting a quadric is described in Section 5, with a discussion
on how to obtain the primitive surfaces from the parameters.

  2. 2.

Combining curves. A pair of point clouds are derived from the user sketches,
representing a pair of curves in 3D. The curves are then combined to form a
surface through;

    1. (a)

Perpendicular extrusion, one of the curves is copied along the other to form
the vertices of a surface;

    2. (b)

or, Interpolation, where the curves form the ends of the surface and
intermediate vertices are interpolated.

The latent variable model for the curves is described in Section 6, from which
the surfaces defined above are created.

Figure 4: The graphical user interface.

The graphical user interface is displayed in Figure 4. The user is presented
with a panel for sketching on camera views on the left, and with the point
cloud on the right. Some examples of fitting primitives to the point cloud
using the interface are given in Figure 2. The top row shows an ellipsoid
being fit to the data from two curves on two views, and the bottom row shows a
perpendicular extrusion. The process is refined over a number of iterations to
produce an output such as on the left of Figure 3. The surface is saved as a
collection of meshes which are then processed by a Poisson surface
reconstruction algorithm, in order to fill the gaps between the surfaces;
displayed on the right hand side of Figure 3.

The method of free-form point selection is first described in Section 4,
before the various fitting algorithms. The curves are modelled using a latent
variable model, this is described in Section 6 along with the methods for
reconstructing surfaces from curves. The ellipsoid and cylinder fitting
algorithms are described in Section 5.

##  4 Free-Form Point Selection

The problem of selecting points is difficult if both the input and the point
cloud is contains noise. It is natural to appeal to Bayesian statistics, since
it is a well established framework for these kinds of problem. Although there
are many different types of noise, it is assumed that the user-input and point
cloud exhibit additive Gaussian noise, and the corresponding points that were
selected in the cloud are extracted by thresholding the probability
distribution.

The interface takes a pen width and colour, and allows a user to draw a
profile across the image. By marking the image the user is supplying the
likelihood that a pixel belongs to a group of 3D points that are of interest
in the point cloud, either because the pixels represent the boundary between
surfaces or fall within the interior of the surface. It is shown how to
combine these distributions in order to derive a probability that each point
in the cloud belongs to a sketch. A threshold can then be applied to the
density in order to extract points of interest.

Let the collection of all images be denoted $F=\\{f_{1},\ldots,f_{N}\\}$. The
subset of images that have been marked are denoted $U$, so that $i\in U$
refers to image $f_{i}$. Typically only a few images are required for marking
and so the size of $U$ is small compared to $F$, particularly in the case of a
video.

Points along the user sketch are denoted $\bm{\mu}_{ij}$; this is the $j^{th}$
point in image $f_{i}$. Each point has an associated covariance
$\Sigma_{ij}=\sigma^{2}_{ij}I$; the $\sigma_{ij}$ represents the uncertainty
in the estimated edge location and are chosen to be the brush width in the
user interface. If we now consider an arbitrary pixel
$\bm{y}\in\mathbb{N}^{2}$, in an image. the likelihood that this pixel belongs
to the mark drawn by the user is modelled as a Gaussian mixture,

$\displaystyle
p(\bm{y}|f_{i})=\frac{1}{K_{i}}\sum_{j=1}^{K_{i}}\mathcal{N}(\bm{y}|\bm{\mu}_{ij},\Sigma_{ij})$
(1)

Now suppose that we have an arbitrary 3D point $\bm{z}\in\mathbb{R}^{3}$, and
a collection of cameras $\Pi=\\{\pi_{1},\ldots,\pi_{|U|}\\}$,

$\displaystyle\pi_{i}:\mathbb{R}^{3}\rightarrow$ $\displaystyle\mathbb{R}^{2}$
(2) $\displaystyle\bm{z}\mapsto$ $\displaystyle\frac{\left[\begin{array}[]{c c
c}1&0&0\\\ 0&1&0\end{array}\right]P_{i}\left[\begin{array}[]{c}\bm{z}\\\
1\end{array}\right]}{(P_{i}\left[\begin{array}[]{c}\bm{z}\\\
1\end{array}\right])_{4}}.$ (9)

where $P_{i}\in\mathbb{R}^{3\times 4}$ is the $i^{th}$ camera matrix.

We wish to find the overall probability of the point in the cloud
$\bm{z}_{k}$, where $0<k<N$. In order to do it we assume that the sketches are
conditionally independent given the projection functions and the 3D point, in
which case the distribution can be written as follows,

$\displaystyle
p(\bm{z}_{k}|\Pi)=\frac{1}{\mathcal{C}}\prod_{i=1}^{|U|}p(\pi_{i}(\bm{z}_{k})|f_{i})$
(10)

The normalising factor $\mathcal{C}$ is then the sum over each point in the
cloud,

$\displaystyle\mathcal{C}=\sum_{k=1}^{N}\left(\prod_{i=1}^{|U|}p(\pi_{i}(\bm{z}_{k})|f_{i})\right)$
(11)

We use this probability to identify 3D points among the million or so points
in the cloud. Every point in the cloud is assigned a probability, many of them
very low. We simply keep all of the points in the cloud which have probability
larger than the mean probability.

It would be possible to include further interaction to identify the surface
normals. This is not necessary for this method since the positions of the
cameras are known. The interior of the surface is set so that all of the
normals are pointing towards a candidate camera.

##  5 Fitting a Quadric

The first category of point extraction is termed ‘whole surface fitting’, in
which the user selects a region which contains points that belong to a whole
surface. Either an ellipsoid or a cylinder is then fitted. Both of these
surface types can be neatly represented as a quadric. A full and detailed
explanation of the method can be found in [Beale et al.(2016)Beale, Yang,
Cosker, and Campbell]. The salient points are described in this section, along
with the method that is used for trimming the final surfaces.

Quadrics are implicit functions, specifically, solutions of the equation,

$\displaystyle\bm{z}^{T}A\bm{z}+\bm{b}^{T}\bm{z}+c=0$ (12)

with $A\in\mathbb{R}^{3\times 3}$ matrix which is symmetric positive definite,
$\bm{b}\in\mathbb{R}^{3}$ a vector and $c\in\mathbb{R}$ a scalar.

In three dimensions $A$ and $\bm{b}$ can be written as follows,

$\displaystyle A=\left[\begin{array}[]{c c c}a_{11}&a_{12}&a_{13}\\\
a_{21}&a_{22}&a_{23}\\\
a_{31}&a_{32}&a_{33}\end{array}\right]\bm{b}=\left[\begin{array}[]{c}b_{1}\\\
b_{2}\\\ b_{3}\end{array}\right]$ (19)

We let,

$\displaystyle\bm{\psi}=[a_{11},a_{22},a_{33},a_{12},a_{13},a_{23},b_{1},b_{2},b_{3},c]^{T}$
(20)

Suppose that the matrix of data points are denoted
$Z=[\bm{z}_{1},\ldots,\bm{z}_{N}]$ with the $k^{th}$ column denoted
$\bm{z}_{k}\in\mathbb{R}^{3}$, and let the $j^{th}$ element of $\bm{z}_{k}$ be
$z_{jk}$. We create a vector of monomials for each data point as follows,

$\displaystyle\bm{x}_{i}$ $\displaystyle=[$ $\displaystyle
z_{1i}^{2},z_{2i}^{2},z_{3i}^{2},2z_{1i}z_{2i},2z_{1i}z_{3i},2z_{2i}z_{3i},$
(21) $\displaystyle z_{1i},z_{2i},z_{3i},1]^{T}$

If we define a monomial data matrix $X=[\bm{x}_{1},\ldots,\bm{x}_{N}]$, the
maximum a-priori solution for the parameters is,

$\displaystyle\bm{\psi}=\left(\sigma^{2}X^{T}X+\left[\begin{array}[]{cc}J&\bm{0}\\\
\bm{0}&\bm{0}\end{array}\right]\right)^{-1}\left[\begin{array}[]{c}\bm{1}\\\
\bm{0}\end{array}\right]$ (26)

Where $\sigma$ is a prior parameter which controls how close to a sphere the
quadric is. The smaller it is, the more likely that the quadric will take the
shape of a sphere. This method is useful since it prevents the quadric from
forming a hyperboloid which is made of two disconnected surfaces. Usually when
fitting primitive geometry we are looking for a single closed surface, and the
best fit to be a contiguous collection of points.

In practise, the surfaces must be trimmed so that they do not go too far
beyond the point cloud. This is done by considering the eigen-decomposition of
the matrix $A$. Considering the expansion,

$\displaystyle(\bm{x}-\bm{\mu})^{T}A(\bm{x}-\bm{\mu})=\tau$ (27)
$\displaystyle\iff$
$\displaystyle\bm{x}^{T}A\bm{x}-2\bm{\mu}^{T}A\bm{x}+\bm{\mu}^{T}A\bm{\mu}$
(28)

and taking $\bm{b}=-2A\bm{\mu}$ and $c=\bm{\mu}^{T}A\bm{\mu}$, it can be seen
that Equations 27 and 12 are equivalent. Removing the mean and projecting the
points onto the principal vectors provides a way to determine how much of the
quadric is covered by the points. To trim the final surface we simply remove
vertices and faces which are outside the maximum and minimum projected values.

|  |
---|---|---
(a) Point cloud |  (b) Points selected by user sketches |  (c) Recovered surface
Figure 5: An example of fitting a quadric to a point cloud. (a) The point
cloud. (b) Points selected by user sketches are shown in green. (c) The final
quadric surface fit to the point cloud. The surface extends to the extremal
points selected by the users sketches or curves.

Figure 5 shows an example of fitting a quadric surface to a set of 3D points
selected from the point cloud.

##  6 Fitting a Latent Variable Model

The second categories of point extraction is termed ‘combining curves’. In
this mode the user selects a pair of curves in the interface which are then
combined either through interpolation of perpendicular extrusion in order to
create the surface. This section describes the method used for fitting a
parametric function to the curves and then combining them in to a single
surface.

One of the principal difficulties in fitting a parametric model to a point
cloud is that the points are noisy and unordered. We suppose that the points
relate to the output of a continuous function, of a single variable. In other
words given a value $t$ we would like to be able to compute an arbitrary point
on the curve. If we know how the value of $t$ relates to the output, the
problem is one of standard regression, if we do not, the $t$ values are
unknown or hidden, and are termed ‘Latent Variables’. A number of different
methods for extracting latent variables and fitting parameters are documented
in [Lawrence(2004)]. Following the work of [Bishop et al.(1998)Bishop,
Svensén, and Williams], we make use of the Generative Topographic Map (GTM),
since it is flexible to a range of different types of function, requiring only
a finite set of basis vectors.

There are several different types of curve, such as B-Splines, Bezier curves,
or linear combinations of sinusoids or exponential kernels. A multidimensional
polynomial is chosen since it is simple and portable. Letting
$\bm{\psi}_{d}\in\mathbb{R}^{L}$ denote the basis vectors for $0<d<D$, and the
matrix $\Phi=[\bm{\phi}_{0},\ldots,\bm{\phi}_{D}]$. For a polynomial basis
this is written as follows,

$\displaystyle\Phi=\left[\begin{array}[]{c c c c c}1&1&1&\ldots&1\\\
0&1&2&\ldots&D\\\ \vdots&\vdots&\vdots&&\vdots\\\
0^{L}&1^{L}&2^{L}&\ldots&D^{L}\end{array}\right]$ (33)

Given a weight matrix $W\in\mathbb{R}^{3\times L}$, the product $W\Psi$ is
then a collection of polynomials evaluated at each of the points $0,\ldots,D$;
with coefficients in $W$. The value of $D$ can be chosen by the user, but is
generally fixed at about $50$. The values of $W\Psi$ are mapped in to the 3D
point cloud space, defining the equidistant centres of a Gaussian mixture
along the polynomial. The model is best represented by Figure 6 (right) where
the black points are unordered 2D data distributed along the curve, and the
curve itself is shown in blue.

Figure 6: An example of multivariate polynomial regression using both
Principle Component Analysis (left) and a Generative Topographic Map (right).
The black points are the noisy point cloud in 2D, the blue ellipses represent
Gaussians.

The constrained Gaussian mixture model, for the collection of points
$Z=[Z_{1},\ldots,Z_{N}]$ is given by,

$\displaystyle
p(Z|\Phi,W)=\prod_{k=1}^{N}\left(\sum_{j=1}^{D}\omega_{k}\mathcal{N}(Z|W\bm{\phi}_{j},\sigma^{2}I)\right)$
(35)

The mixture is solved for the parameters $W$ and $\sigma$ by finding the
derivatives of the log-likelihood of Equation 35. The solution is a fixed
point iteration, which can be formulated as an Expectation-Maximisation
algorithm as follows. We first compute the expected probability of a point
given the parameters and the data,

$\displaystyle
R_{ij}=\frac{\omega_{i}\mathcal{N}(Z_{j}|W\bm{\phi}_{i}.\sigma^{2}I)}{\sum_{m=1}^{D}\omega_{m}\mathcal{N}(Z_{j}|W\bm{\phi}_{m},\sigma^{2}I)}$
(36)

The maximisation step finds the maximum likelihood values for the parameters
given the expected probabilities,

$\displaystyle W^{T}_{new}=$ $\displaystyle(\Phi G\Phi^{T})^{-1}\Phi RZ^{T}$
(37) $\displaystyle\sigma^{2}=$
$\displaystyle\sum_{i=1}^{D}\sum_{j=1}^{K}R_{ij}||W_{new}\bm{\phi}-Z_{j}||^{2}$
(38)

where G is a diagonal matrix with elements $G_{ii}=\sum_{n}R_{in}$.

Since this procedure is iterative it requires an initialisation of the
variables $W$ and $\sigma$. This is done using principle component analysis.
Consider a line defined by the function $f(t)=\bm{n}t+\bm{a}$, and a Gaussian
which has been fit to the data defined by by parameters $\bm{\mu}$ and
$\Sigma$. We would like the line to pass through the mean of the data and
along the largest principle component of $\Sigma$. This is shown graphically
on the left of Figure 6, the unordered data is represented by the black
circles and we would like to fit the line represented by the blue dashes. The
parametrisation of $f$ makes makes this easy: set $\bm{a}=\bm{\mu}$ and let
$\bm{n}$ be the vector corresponding to the largest eigenvalue of $\Sigma$. In
the notation of the above the weight matrix is then taken to be
$W=[\bm{a},\bm{n},\bm{0},\ldots,\bm{0}]$, and the variance is initialised as
the smallest eigenvalue of the covariance matrix.

Generally speaking, the ends of the curve do not need to be computed since the
algorithm converges to fit the points precisely. It is possible that the shape
of the curve causes the Gaussians to be centres beyond the points. In this
case the GTM is trimmed by removing components at either end which have a
value of $\frac{1}{K}\sum_{k=1}^{K}R_{ij}$ that falls below a threshold, or in
other words, Gaussian components are removed from the ends which do not have
any points associated with them. Doing this keeps the resulting curve tight to
the original data, even after convergence.

The final surface is built from the curves. Either the surface is an
interpolation of two drawn edge curves, or one curve is extrapolated in the
direction of a perpendicular second curve. In practise, drawing a profile and
extrapolating it along an edge gives very good results; a similar approach was
taken in 3-sweep [Chen et al.(2008)Chen, Kang, Xu, Dorsey, and Shum].

Letting $\bm{q}_{i},\bm{p}_{i}\in\mathbb{R}^{3}$ denote two extracted curves,
and $\bm{s}_{kj}$ be a set of vertices on the desired surface, for
$0<i,j,k\leq D$. We first reorder the curves so that they are oriented in the
same direction; if $||\bm{q}_{0}-\bm{p}_{D}||<||\bm{q}_{0}-\bm{p}_{0}||$ then
the curve $(\bm{p}_{j})_{j=1}^{D}$ is reversed. The interpolation is then
performed using the following equation,

$\displaystyle\bm{s}_{kj}=\left(\frac{j}{D}\right)\bm{q}_{k}+\left(1-\frac{j}{D}\right)\bm{p}_{k}$
(39)

If the curves are perpendicular then the curves are reordered in the same way
as above. One side of the surface is initialised with a trajectory and then
copied and translated it according to the velocities of the second. This
method avoids the requirement of computing curve intersections, which may be
difficult if the data is noisy. The surface vertices are as follows,

$\displaystyle\bm{s}_{kj}=\left\\{\begin{array}[]{r l}\bm{q}_{k}&:j=0\\\
\bm{q}_{k}+\bm{p}_{j}-\bm{p}_{j-1}&:j>0\end{array}\right.$ (42)

|  |
---|---|---
(a) First edge |  (b) Second edge |  (c) Final surface
Figure 7: An example of fitting a perpendicular GTM surface from a pair of
curves. A single collection of sketches define each of the edges, and the
first edge is extruded along the second edge to create the final surface.

An example of fitting a perpendicular GTM surface is shown in Figure 7.

##  7 Finalising the surface

The system described in previous sections provides a way for a user to fit
geometric primitives to a point cloud. It remains to fill in the holes in the
surfaces.

To close the holes between the primitive surfaces the Poisson reconstruction
algorithm is used, as described in [Kazhdan et al.(2006)Kazhdan, Bolitho, and
Hoppe]. An open source implementation is available in the software package
MeshLab [Cignoni et al.(2008)Cignoni, Corsini, and Ranzuglia]. Since the
output from the software is a collection of meshes with normals no pre-
processing is required, and a closed surface is fit to the geometry. The mesh
is completed by removing all of the vertices and faces which have long edges.
This is possible since the Poisson reconstruction tends to produce large
triangles in regions where it is predicting.

##  8 Results

The results in this section present comparisons against state of the art
alternatives for reconstructing surfaces from point clouds. It is difficult to
assess the quality of the output without first establishing a high quality
ground truth. The problem of overfitting in statistical regression motivates
an argument that, even with a good quality ground truth, the quality of a
result is difficult to measure, and ultimately lies with the perceiver.
Another difficulty in assessing the results is in the human interface. The
competing methods require user interaction in order to fit the surfaces.
Rather than make arguments based on the effort or complexity of the user
input, a single working example is chosen with rational arguments about the
limitations of the methods. It is argued that each of the competing softwares
are fit for the purpose that they were designed, but fail on the dataset that
we have, making the method presented in this paper a reasonable replacement.

The interactive methods that compared against are, VideoTrace [van den Hengel
et al.(2007)van den Hengel, Dick, Thormählen, Ward, and Torr] and Morfit [Yin
et al.(2014)Yin, Huang, Zhang, Gong, Cohen-Or, and Chen], both of which are
contemporary, and present an alternative approach to surface reconstruction.
The comparison is made on a point cloud created from photographs of a building
with columns and a non-planar facade. A large portion of the entrance is
occluded and there are holes in the point cloud as a result of it. These
problems make it very difficult to automatically reconstruct the surface, and
a user must provide some idea of where the surfaces should be. This is the
example that is used to make the initial comparisons in Sections 8.1 and 8.2.

The principal argument for using the method presented in this paper is that a
large percentage of man-made objects are composed of a collection primitive
surfaces. The final set of results show the performance on a dataset
containing 6 examples of buildings and outdoor furniture, the comparisons
provide a visual assessment of the improvement beyond the standard method of
normal estimation and Poisson reconstruction.

###  8.1 Comparison against Morfit

The Morfit [Yin et al.(2014)Yin, Huang, Zhang, Gong, Cohen-Or, and Chen]
software package generates a surface based on the morphological skeleton of
the object. The software allows a user to edit the skeleton, and then sketch
extra vertices and normals on a plane which is orthogonal to an arbitrary
point on the skeleton, before fitting surface. So long as the skeleton is a
good enough prediction the surface will allow a user to specify the surface
shape in regions where there are large holes.

Point cloud | Skeleton
---|---
|
Figure 8: An example of a skeleton fitted to point cloud of a building with a
non-planar and partially occluded facade.

The principle limitation of Morfit is that it requires a morphological
skeleton in order to compute the surface. Consider the point cloud shown in
Figure 8. The building has 6 columns arranged in a cylinder, the occluded
entrance is cylindrical and each of the walls are planar. The right hand side
of Figure 8 shows the automatically extracted skeleton which correctly
identified the two front columns, but was unable to find a skeleton for other
features such as the back wall. In fact, a skeleton does exist for this
building. It is part of a larger closed structure which is only partially
visible. Using the Morfit interface to fix the skeleton is extremely time
consuming, and running the sweeper on the resulting branches yields poor
results.

The dataset that Morfit was tested on is primarily made up of models which are
small and closed, such as the model Gekko featured throughout the paper. The
algorithm works well on this data, since the skeleton is much easier to
compute and the reconstructed points are evenly distributed along the object.

###  8.2 Comparison against VideoTrace

VideoTrace is an interactive software package for reconstructing 3D surfaces
from a collection of photographs of an object from multiple views. It first
reconstructs a sparse point cloud and then allows a user to draw a polygon or
non-uniform rational b-spline in three dimensions using the interface. The
curves can be joined to create a surfaces, which can then either be extruded
or mirrored two create a more complicated closed surface.

VideoTrace
---
View 1 | View 2 | Reconstruction
|  |
|  |

This system
---
View 1 | View 2 | Reconstruction
|  |

Figure 9: Example output from video trace. The first row shows a successful
surface reconstruction of an occluded face, and the second row shows an
unsuccessful reconstruction of the same face. The final row is the output from
the system presented in this paper.

One of the principal limitations of this system is the way in which the user
identifies the 3D surface. Assuming that the camera geometry and point cloud
have been accurately reconstructed, the system takes a polygon from a single
view and uses the point cloud to infer where the polygon lies in 3D space. The
user can then refine the fit by moving control points from different views.
This presents a problem if there is a surface which is occluded, such as the
post office canopy in Figure 9. The column has been highlighted for
comparison. On the top row the occluded surface was identified by creating a
polygon on the left hand image, the computed geometry was correct in the
second view, and also in the 3D reconstruction on the right hand side. The
second row shows this method to fail when the user wants to extract the larger
segment on the occluded face. Again, the polygon was identified in the left
hand column but is significantly erroneous in both the second view and the 3D
reconstruction on the right hand side. The method presented in this paper
correctly identifies the surface, even though we have sketched over an
occlusion, and is displayed on the bottom row.

###  8.3 Qualitative results

We provide some final comparisons on a dataset of objects, all of which are
man-made objects in an urban environment. The dataset contains building which
have a surface structure which is not restricted to planes alone, making a
large number of traditional methods inapplicable.

The results can be visually compared in Figure 10. It can be seen that the
Poisson reconstruction fails where there are large holes in the surface such
as the floor of the post-office, or on the roof of the phone booths. The final
textured meshes are shown in Figure 11.

| Point cloud | Poisson reconstruction | Primitives | Closure
---|---|---|---|---
The cross |  |  |  |
The post office |  |  |  |
Mantova Basilica |  |  |  |
Phone booth |  |  |  |
Mantova Piazza |  |  |  |
Post box |  |  |  |
Figure 10: Some results from the method presented in this paper after fitting the primitives and applying surface closure. The far left shows the original point cloud, the second column shows the results of normal estimation and Poisson reconstruction, and the right hand side shows the extracted primitives, and the result of applying a Poisson reconstruction to the output.  |  |
---|---|---
|  |
Figure 11: The final textured results.

The software is demonstrated in the supplementary material. A candidate user
is recorded fitting primitives using only free-form input. The videos
demonstrate the ease of use for both latent variable models and quadric
surfaces, on a touch screen display. The use of Bayesian statistics makes the
system robust to noise in the input trajectories, resulting in geometry which
is a good fit to the data.

##  9 Conclusion

The system that has been presented in this paper enables a user to quickly fit
parametric surfaces to point clouds. It is demonstrated that we outperform
alternative methods on a dataset of medium to large scale objects such as
buildings and street furniture, which are man-made and comprise surfaces which
are piecewise primitive. Alternative interactive methods either do not make
the correct assumptions about the data or do not correctly deal with occluded
surfaces, making them laborious for desktop use by an artist.

The system performs well on data which exhibits a large amount of noise, and
has holes in the point cloud; a limiting artefact for modern autonomous
reconstruction techniques such as Poisson reconstruction.

While the system does perform well under the assumptions, it is not expected
to work well if they are broken. Highly complicated natural objects such as
flowing water, vegetation and so on, require an algorithm which does not
constrain the smoothness of the surface to the same degree. To some extent,
simplicity is a matter of scale, and so at a close range one might expect
better results; but it is an extension that we do not consider in this work.

We conclude that the system provides a useful interface for reconstructing
surfaces from urban scenes, and produces high quality and visually pleasing
results.

## References

  * [Andrews and Séquin(2014)] James Andrews and Carlo H Séquin.  Type-constrained direct fitting of quadric surfaces.  _Computer-Aided Design and Applications_ , 11(1):107–119, 2014.
  * [Beale et al.(2016)Beale, Yang, Cosker, and Campbell] Daniel Beale, Yongliang Yang, Darren Cosker, and Neill Campbell.  Fitting quadrics with a bayesian prior.  _Computational Visual Media_ , 2016.
  * [Bishop et al.(1998)Bishop, Svensén, and Williams] Christopher M Bishop, Markus Svensén, and Christopher KI Williams.  Gtm: The generative topographic mapping.  _Neural computation_ , 10(1):215–234, 1998.
  * [Bo et al.(2012)Bo, Ling, and Wang] Pengbo Bo, Ruotian Ling, and Wenping Wang.  A revisit to fitting parametric surfaces to point clouds.  _Computers & Graphics_, 36(5):534–540, 2012\.
  * [Bolitho et al.(2009)Bolitho, Kazhdan, Burns, and Hoppe] Matthew Bolitho, Michael Kazhdan, Randal Burns, and Hugues Hoppe.  Parallel poisson surface reconstruction.  In _Advances in Visual Computing_ , pages 678–689. Springer, 2009\.
  * [Campbell et al.(2010)Campbell, Vogiatzis, Hernández, and Cipolla] Neill DF Campbell, George Vogiatzis, Carlos Hernández, and Roberto Cipolla.  Automatic 3d object segmentation in multiple views using volumetric graph-cuts.  _Image and Vision Computing_ , 28(1):14–25, 2010\.
  * [Cao et al.(2014)Cao, Ju, Fu, and Hu] Yan-Pei Cao, Tao Ju, Zhao Fu, and Shi-Min Hu.  Interactive image-guided modeling of extruded shapes.  In _Computer Graphics Forum_ , volume 33, pages 101–110. Wiley Online Library, 2014.
  * [Changchang(2011a)] W. Changchang.  Visualsfm: A visual structure from motion system.  _http://ccwu.me/vsfm/_ , 2011a.
  * [Changchang(2011b)] W. Changchang.  Multicore bundle adjustment.  _Computer Vision and Pattern Recognition (CVPR)_ , 2011b.
  * [Chen et al.(2013)Chen, Zhu, Shamir, Hu, and Cohen-Or] Tao Chen, Zhe Zhu, Ariel Shamir, Shi-Min Hu, and Daniel Cohen-Or.  3-sweep: extracting editable objects from a single photo.  _ACM Transactions on Graphics (TOG)_ , 32(6):195, 2013.
  * [Chen et al.(2008)Chen, Kang, Xu, Dorsey, and Shum] Xuejin Chen, Sing Bing Kang, Ying-Qing Xu, Julie Dorsey, and Heung-Yeung Shum.  Sketching reality: Realistic interpretation of architectural designs.  _ACM Transactions on Graphics (TOG)_ , 27(2):11, 2008.
  * [Cignoni et al.(2008)Cignoni, Corsini, and Ranzuglia] Paolo Cignoni, Massimiliano Corsini, and Guido Ranzuglia.  Meshlab: an open-source 3d mesh processing system.  _Ercim news_ , 73(45-46):6, 2008.
  * [Debevec et al.(1996)Debevec, Taylor, and Malik] Paul E Debevec, Camillo J Taylor, and Jitendra Malik.  Modeling and rendering architecture from photographs: A hybrid geometry-and image-based approach.  In _Proceedings of the 23rd annual conference on Computer graphics and interactive techniques_ , pages 11–20. ACM, 1996.
  * [Dragiev et al.(2011)Dragiev, Toussaint, and Gienger] Stanimir Dragiev, Marc Toussaint, and Michael Gienger.  Gaussian process implicit surfaces for shape estimation and grasping.  In _Robotics and Automation (ICRA), 2011 IEEE International Conference on_ , pages 2845–2850. IEEE, 2011.
  * [Furukawa(2010)] J. Furukawa, Y. Ponce.  Accurate, dense and robust multi-view stereopsis.  _IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI)_ , 32(8):1362–1378, 2010.
  * [Furukawa et al.(2009)Furukawa, Curless, Seitz, and Szeliski] Yasutaka Furukawa, Brian Curless, Steven M Seitz, and Richard Szeliski.  Reconstructing building interiors from images.  In _Computer Vision, 2009 IEEE 12th International Conference on_ , pages 80–87. IEEE, 2009.
  * [Gopi et al.(2000)Gopi, Krishnan, and Silva] M Gopi, Shankar Krishnan, and Cláudio T Silva.  Surface reconstruction based on lower dimensional localized delaunay triangulation.  In _Computer Graphics Forum_ , volume 19, pages 467–478. Wiley Online Library, 2000.
  * [Habbecke and Kobbelt(2012)] Martin Habbecke and Leif Kobbelt.  Linear analysis of nonlinear constraints for interactive geometric modeling.  In _Computer Graphics Forum_ , volume 31, pages 641–650. Wiley Online Library, 2012.
  * [Hoppe et al.(1992)Hoppe, DeRose, Duchamp, McDonald, and Stuetzle] Hugues Hoppe, Tony DeRose, Tom Duchamp, John McDonald, and Werner Stuetzle.  _Surface reconstruction from unorganized points_ , volume 26.  ACM, 1992.
  * [Iarussi et al.(2015)Iarussi, Bommes, and Bousseau] Emmanuel Iarussi, David Bommes, and Adrien Bousseau.  Bendfields: Regularized curvature fields from rough concept sketches.  _ACM Transactions on Graphics (TOG)_ , 34(3):24, 2015.
  * [Jiang et al.(2009)Jiang, Tan, and Cheong] Nianjuan Jiang, Ping Tan, and Loong-Fah Cheong.  Symmetric architecture modeling with a single image.  _ACM Transactions on Graphics (TOG)_ , 28(5):113, 2009.
  * [Kazhdan and Hoppe(2013)] Michael Kazhdan and Hugues Hoppe.  Screened poisson surface reconstruction.  _ACM Transactions on Graphics (TOG)_ , 32(3):29, 2013.
  * [Kazhdan et al.(2006)Kazhdan, Bolitho, and Hoppe] Michael Kazhdan, Matthew Bolitho, and Hugues Hoppe.  Poisson surface reconstruction.  _Proceedings of the fourth Eurographics symposium on Geometry processing_ , 2006.
  * [Kolmogorov(2004)] Vladimir Kolmogorov.  _Graph based algorithms for scene reconstruction from two or more views_.  PhD thesis, Citeseer, 2004.
  * [Krishnamurthy and Levoy(1996)] Venkat Krishnamurthy and Marc Levoy.  Fitting smooth surfaces to dense polygon meshes.  In _Proceedings of the 23rd annual conference on Computer graphics and interactive techniques_ , pages 313–324. ACM, 1996.
  * [Labatut et al.(2009)Labatut, Pons, and Keriven] Patrick Labatut, J-P Pons, and Renaud Keriven.  Hierarchical shape-based surface reconstruction for dense multi-view stereo.  In _Computer Vision Workshops (ICCV Workshops), 2009 IEEE 12th International Conference on_ , pages 1598–1605. IEEE, 2009.
  * [Lau et al.(2010)Lau, Saul, Mitani, and Igarashi] Manfred Lau, Greg Saul, Jun Mitani, and Takeo Igarashi.  Modeling-in-context: user design of complementary objects with a single photo.  In _Proceedings of the Seventh Sketch-Based Interfaces and Modeling Symposium_ , pages 17–24. Eurographics Association, 2010.
  * [Lawrence(2004)] Neil D Lawrence.  Gaussian process latent variable models for visualisation of high dimensional data.  _Advances in neural information processing systems_ , 16:329–336, 2004.
  * [Levin et al.(2004)Levin, Lischinski, and Weiss] Anat Levin, Dani Lischinski, and Yair Weiss.  Colorization using optimization.  _ACM Trans. Graph._ , 23(3):689–694, August 2004\.
  * [Li et al.(2011)Li, Wu, Chrysathou, Sharf, Cohen-Or, and Mitra] Yangyan Li, Xiaokun Wu, Yiorgos Chrysathou, Andrei Sharf, Daniel Cohen-Or, and Niloy J Mitra.  Globfit: Consistently fitting primitives by discovering global relations.  In _ACM Transactions on Graphics (TOG)_ , volume 30, page 52. ACM, 2011.
  * [Lipp et al.(2014)Lipp, Wonka, and Müller] Markus Lipp, Peter Wonka, and Pascal Müller.  Pushpull++.  _ACM Transactions on Graphics (TOG)_ , 33(4):130, 2014.
  * [Mitra and Nguyen(2003)] Niloy J Mitra and An Nguyen.  Estimating surface normals in noisy point cloud data.  In _Proceedings of the nineteenth annual symposium on Computational geometry_ , pages 322–328. ACM, 2003.
  * [Monszpart et al.(2015)Monszpart, Mellado, Brostow, and Mitra] Aron Monszpart, Nicolas Mellado, Gabriel J Brostow, and Niloy J Mitra.  Rapter: rebuilding man-made scenes with regular arrangements of planes.  In _ACM SIGGRAPH_ , volume 2015, 2015.
  * [Musialski et al.(2013)Musialski, Wonka, Aliaga, Wimmer, Gool, and Purgathofer] Przemyslaw Musialski, Peter Wonka, Daniel G Aliaga, Michael Wimmer, L Gool, and Werner Purgathofer.  A survey of urban reconstruction.  In _Computer Graphics Forum_ , volume 32, pages 146–177. Wiley Online Library, 2013.
  * [Nan et al.(2010)Nan, Sharf, Zhang, Cohen-Or, and Chen] Liangliang Nan, Andrei Sharf, Hao Zhang, Daniel Cohen-Or, and Baoquan Chen.  Smartboxes for interactive urban reconstruction.  _ACM Transactions on Graphics (TOG)_ , 29(4):93, 2010.
  * [Portelli et al.(2010)Portelli, Ganovelli, Tarini, Cignoni, Dellepiane, and Scopigno] Davide Portelli, Fabio Ganovelli, Marco Tarini, Paolo Cignoni, Matteo Dellepiane, and Roberto Scopigno.  A framework for user-assisted sketch-based fitting of geometric primitives.  2010\.
  * [Schnabel et al.(2007)Schnabel, Wahl, and Klein] Ruwen Schnabel, Roland Wahl, and Reinhard Klein.  Efficient ransac for point-cloud shape detection.  In _Computer graphics forum_ , volume 26, pages 214–226. Wiley Online Library, 2007.
  * [Sinha et al.(2008)Sinha, Steedly, Szeliski, Agrawala, and Pollefeys] Sudipta N Sinha, Drew Steedly, Richard Szeliski, Maneesh Agrawala, and Marc Pollefeys.  Interactive 3d architectural modeling from unordered photo collections.  In _ACM Transactions on Graphics (TOG)_ , volume 27, page 159. ACM, 2008.
  * [Starck et al.(2006)Starck, Miller, and Hilton] Jonathan Starck, Gregor Miller, and Adrian Hilton.  Volumetric stereo with silhouette and feature constraints.  In _British machine vision conference_ , pages 1189–1198, 2006.
  * [Subr et al.(2013)Subr, Paris, Soler, and Kautz] Kartic Subr, Sylvain Paris, Cyril Soler, and Jan Kautz.  Accurate binary image selection from inaccurate user input.  _Computer Graphics Forum_ , 32(2pt1):41–50, 2013\.
  * [van den Hengel et al.(2007)van den Hengel, Dick, Thormählen, Ward, and Torr] Anton van den Hengel, Anthony Dick, Thorsten Thormählen, Ben Ward, and Philip HS Torr.  Videotrace: rapid interactive scene modelling from video.  volume 26, page 86. ACM, 2007.  ISBN 0730-0301.
  * [Wang et al.(2006)Wang, Pottmann, and Liu] Wenping Wang, Helmut Pottmann, and Yang Liu.  Fitting b-spline curves to point clouds by curvature-based squared distance minimization.  _ACM Transactions on Graphics (ToG)_ , 25(2):214–238, 2006.
  * [Williams and Fitzgibbon(2007)] Oliver Williams and Andrew Fitzgibbon.  Gaussian process implicit surfaces.  _Gaussian Proc. in Practice_ , 2007.
  * [Xu et al.(2014)Xu, Chang, Sheffer, Bousseau, Mccrae, and Singh] Baoxuan Xu, William Chang, Alla Sheffer, Adrien Bousseau, James Mccrae, and Karan Singh.  True2form: 3d curve networks from 2d sketches via selective regularization.  _ACM Transactions on Graphics_ , 33(4), 2014.
  * [Yin et al.(2014)Yin, Huang, Zhang, Gong, Cohen-Or, and Chen] Kangxue Yin, Hui Huang, Hao Zhang, Minglun Gong, Daniel Cohen-Or, and Baoquan Chen.  Morfit: Interactive surface reconstruction from incomplete point clouds with curve-driven topology and geometry control.  _ACM Transactions on Graphics_ , 33(6), 2014.