(wacv) Package wacv Warning: Single column document - WACV requires papers to
have two-column layout. Please load document class ‘article’ with ‘twocolumn’
option

# So you think you can track? Supplementary material

The following appendices are included in this file to supplement the main
paper.

  1. 1.

Dataset Coordinate Systems

  2. 2.

Homography Re-estimation Method

  3. 3.

Homography Error Metrics

  4. 4.

Additional GPS Trajectory Plots

  5. 5.

Object Detector

  6. 6.

Experimental Details

  7. 7.

Treatment of Personally Identifiable Information

  8. 8.

Known Data Artifacts and Anomalies

##  1 Dataset Coordinate Systems

The Interstate-24 Video dataset utilizes 3 sets of coordinates utilized
throughout I-24 MOTION system [gloudemans202324] datasets:

  * •

Image Coordinates: are given in pixels. $(y^{(im)},x^{(im)})$ denotes the row
and column of the specified pixel. By convention the top left pixel is (0,0).

  * •

State Plane Coordinates: specify a rectilinear and orthogonal coordinate
system. The EPSG 2274 state plane coordinate system for Tennessee is specified
in feet relative to a known survey point. $(x^{(st)},y^{(st)})$ indicates the
coordinate (in feet) along the first (roughly east-west) and second (roughly
north-south) coordinate axis defined by the state plane coordinate system.
(Note that a common conversion from state plane coordinates to
latitude/longitude coordinates (e.g. WSG84 or NAD83) can be utilized if
desired.) A third orthogonal coordinate axis (z-axis) is defined and
corresponds to distance off the roadway, such that $z^{(st)}=0$ for all points
on the roadway plane.

  * •

Roadway Coordinates: are defined such that the primary (x) axis lies along the
median (or more precisely, midway between the two interior yellow lines for
the interstate) at all points within the instrument extents, and the secondary
(y) axis is defined locally to perpendicular to the primary axis at all points
along the roadway. All coordinates with a distance from the primary axis less
than the local radius of curvature (including all points on the roadway) have
a unique $(x^{(r)},y^{(r)})$ coordinate. By left-hand rule convention, we
define the positive y-axis to be in the direction of the eastbound roadway
lanes at all points along the roadway.

###  1.1 Notation

Throughout the rest of this appendix to disambiguate the various coordinate
systems, the following notation is used:

  * •

$x$,$y$, and $z$ refer to coordinate axes. A superscript $(im)$, $(st)$, or
$(r)$ specifies all variables corresponding to a specific coordinate system
(e.g. $x^{(st)}$).

  * •

Vectors and matrices in a specified coordinate system are denoted in bold
(e.g. $\textbf{O}^{(st)}$).

  * •

Homography matrices are also listed in bold script, without superscript (e.g
H).

  * •

A subscript indexes a specific point (e.g. $x^{(st)}_{bbl}$), and subscript
$i$ indicates an arbitrary element index from a set of elements (e.g $a_{i}$).

  * •

An $x$,$y$, or $z$ without a subscript indicates a generic variable along the
specified axis within the specified coordinate system.

Symbol | Definition
---|---
H | 3$\times$3 matrix of homography parameters $h_{ij}$
P | 3$\times$4 matrix of homography parameters $p_{ij}$
$s$ | homography scale parameter
$x^{(im)}$, $y^{(im)}$ | image coordinates (y indicates pixel row and x indicates pixel column)
$x^{(st)}$, $y^{(st)}$, $z^{(st)}$ | state plane coordinates
$x^{(r)}$, $y^{(r)}$ | roadway coordinates
$\mathbf{O}^{(st)}$ | state plane coordinates for object, equal to [$\mathbf{o}^{(st)}_{bbl}$,$\mathbf{o}^{(st)}_{bbr}$,$\mathbf{o}^{(st)}_{btl}$ ,$\mathbf{o}^{(st)}_{btr}$ ,$\mathbf{o}^{(st)}_{fbl}$ ,$\mathbf{o}^{(st)}_{fbr}$ ,$\mathbf{o}^{(st)}_{ftl}$ ,$\mathbf{o}^{(st)}_{ftr}$]
$\mathbf{o}^{(st)}_{bbl}$ | back bottom left state plane coordinate of object, equal to [$x^{(st)}_{bbl},y^{(st)}_{bbl},z^{(st)}_{bbl}$]
$\mathbf{o}^{(st)}_{c}$ | back bottom center state plane coordinate, primary reference coordinate for the object
$\mathbf{o}^{(st)}_{spl}$ | state plane coordinates of point on center-line spline ($y^{(r)}=0$) with the same $x^{(r)}$ coordinate as $\mathbf{o}^{(st)}_{c}$
$\mathbf{O}^{(r)}$ | roadway coordinates for object, $[x^{(r)_{o}},y^{(r)_{o}},l,w,h]$
$x^{(r)}$ | generic longitudinal roadway coordinate along curvilinear spline axis
$y^{(r)}$ | generic lateral roadway coordinate along axis locally perpendicular to longitudinal roadway coordinate axis
$x^{(r)}_{o}$ | object longitudinal roadway coordinate along curvilinear spline axis
$y^{(r)}_{o}$ | object lateral roadway coordinate along axis locally perpendicular to longitudinal roadway coordinate axis
$l$,$w$,$h$ | rectangular prism dimensions (length, width and height)
$F(x^{(r)})$ | spline defining state plane coordinate roadway center-line spline parameterized by $x^{(r)}$
$\tilde{G}(x^{(st)})$ | spline approximating the center-line spline in roadway coordinates $x^{(r)}$ parameterized by $x^{(st)}$
Table 1: Summary of symbols used in this section.

A list of all variables along with their descriptions is given in Table 1.
Transformations between image and state plane coordinates, and transformations
between state plane and roadway coordinates are detailed in the next two
sections.

###  1.2 Image $\leftrightarrow$ State Plane Conversion

(parts of this subsection rely on similar definitions and descriptions to the
supplement in [gloudemans2023dataset].) A homography relates two views of a
planar surface. For each camera, we provide homography information such that
the 8-corner coordinates of the stored 3D bounding-box annotation can be
projected into any camera view for which the vehicle is visible, creating a
monocular 3D bounding box within that camera field of view. For each direction
of travel in each camera view, for each scene, a homography relating the image
pixel coordinates to the state plane coordinate system is defined. (Though the
same cameras are used for different scenes, the positions of the cameras
changes slightly over time due). A local flat plane assumption is used (the
state plane coordinate system is assumed to be piece-wise flat)
[hartley2003multiple]. A series of correspondence points
$a_{i}=[x^{(im)},y^{(im)},x^{(st)},y^{(st)},z^{(st)}]$ are used to define this
relation, where $(y^{(im)},x^{(im)})$ is the coordinate of selected
correspondence point $a$ in pixel coordinates (row, column) and
$(x^{(st)},y^{(st)},z^{(st)})$ is the selected correspondence point in state
plane coordinates.

All selected points are assumed to lie on the state plane, so $z^{(st)}=0$ for
all selected correspondence points. Visible lane marking lines and other
easily recognizable landmarks on the roadway are used as correspondence points
in each camera field of view. Each correspondence point is also labeled in
global information system (GIS) software, giving the precise GPS / state-plane
coordinate system coordinates for each labeled corresponding point. The
corresponding pixel coordinates are manually selected in each camera field of
view, for each direction of travel on the roadway (see Appendix 1.5).

A perspective transform (Equation 2) is fit to these correspondence points. We
first define a 2D perspective transform which defines a linear mapping
(Equation 1) of points from one plane to another that preserves straight
lines. The correspondence points are then used to solve for the best
perspective transform H as defined in equation 2, where $s$ is a scale factor.

$s\begin{bmatrix}x^{(st)}_{i}\\\ y^{(st)}_{i}\\\
1\end{bmatrix}\sim\textbf{H}\begin{bmatrix}x^{(im)}_{i}\\\ y^{(im)}_{i}\\\
1\end{bmatrix}$ (1)

where H is a $3\times 3$ matrix of parameters:

$\textbf{H}=\begin{bmatrix}h_{11}&h_{12}&h_{13}\\\ h_{21}&h_{22}&h_{23}\\\
h_{31}&h_{32}&h_{33}\end{bmatrix}$ (2)

For each camera field of view and each direction of travel, the best
perspective transform $\textbf{H}^{*}$ is determined by minimizing the sum of
squared re-projection errors according to equation 3 as implemented in
OpenCV’s $find\mathunderscore homography$ function [bradski2000opencv]:

$\textbf{H}^{*}=\operatorname*{arg\,min}_{\textbf{H}}\sum_{i}\left(x^{(st)}_{i}-\frac{h_{11}x^{(im)}_{i}+h_{12}y^{(im)}_{i}+h_{13}}{h_{31}x^{(im)}_{i}+h_{32}y^{(im)}_{i}+h_{33}}\right)^{2}+\left(y^{(st)}_{i}-\frac{h_{21}x^{(im)}_{i}+h_{22}y^{(im)}_{i}+h_{23}}{h_{31}x^{(im)}_{i}+h_{32}y^{(im)}_{i}+h_{33}}\right)^{2}$
(3)

The resulting matrix $\textbf{H}^{*}$ allows any point lying on the plane
within the camera field of view to be converted into state plane coordinates.
The corresponding matrix $\textbf{H}_{inv}$ can easily be obtained to project
state plane coordinates on the $z=0$ plane into image coordinates. However,
since each vehicle is represented by a 3D bounding box, the top corner
coordinates of the box do not lie on the ground plane. A 3D perspective
transform P is needed to linearly map coordinates from 3D state plane
coordinate space to 2D image coordinate space, where P is a $3\times 4$ matrix
of parameters:

$\textbf{P}=\begin{bmatrix}p_{11}&p_{12}&p_{13}&p_{14}\\\
p_{21}&p_{22}&p_{23}&p_{24}\\\ p_{31}&p_{32}&p_{33}&p_{34}\end{bmatrix}$ (4)

and P projects a point in 3D space $(x^{\prime},y^{\prime},z^{\prime})$ into
the corresponding image point $(x,y)$ according to:

$\textbf{P}\begin{bmatrix}x^{(st)}\\\ y^{(st)}\\\ z^{(st)}\\\
1\end{bmatrix}\sim s^{\prime}\begin{bmatrix}x^{(im)}\\\ y^{(im)}\\\
1\end{bmatrix}$ (5)

where $s^{\prime}$ is a new scaling parameter. By observing the case where
$z^{(st)}=0$, it is evident columns 1,2, and 4 of P are equivalent to the
columns of $\textbf{H}_{inv}$ and can be fit in the same way. Thus, we need
only solve for column 3 of P. Next, we note as in [hartley2003multiple] that
$(\frac{p_{11}}{p_{31}},\frac{p_{21}}{p_{31}})$ is the vanishing point (in
image coordinates) of perspective lines drawn in the same direction as the
state plane coordinate x-axis. The same is true for the 2nd column and the
state plane coordinate y-axis, the 3rd column and the state plane coordinate
z-axis, and the 4th column and the state plane coordinate origin.

Thus, to fully determine P it is sufficient to locate the vanishing point of
the z-axis in state plane coordinates and to estimate the scaling parameter
$p_{33}$. The vanishing point is located in image coordinates by finding the
intersection point between lines drawn in the z-direction. Such lines are
obtained by manually annotating vertical lines in each camera field of view.
The scale parameter is estimated by minimizing the sum of squared reprojection
errors defined in equation 1.2 for a sufficiently large set of state plane
coordinates and corresponding, manually annotated coordinates in image space.

$\displaystyle\textbf{P}^{*}=\operatorname*{arg\,min}_{p_{33}}\sum_{i}$
$\displaystyle\left(x^{(im)}_{i}-\frac{p_{11}x^{(st)}_{i}+p_{12}y^{(st)}_{i}+p_{13}z^{(st)}_{i}+p_{14}}{p_{31}x^{(st)}_{i}+p_{32}y^{(st)}_{i}+p_{33}z^{(st)}_{i}+p_{34}}\right)^{2}+$
$\displaystyle\left(y^{(im)}_{i}-\frac{p_{21}x^{(st)}_{i}+p_{22}y^{(st)}_{i}+p_{23}z^{(st)}_{i}+h_{24}}{p_{31}x^{(st)}_{i}+p_{32}y^{(st)}_{i}+p_{33}z^{(st)}_{i}+h_{34}}\right)^{2}$
(6)

The resulting 3D perspective transform $\textbf{P}^{*}$ allows for the
lossless conversion of points in roadway coordinates to the corresponding
points in image coordinates. Observing that a lossless conversion from image
coordinates to state plane coordinates is available provided that the
converted point lies on the $z^{(st)}=0$ plane, it is possible to precisely
convert a rectangular prism from image space to state plane coordinates by i.)
converting the footprint of the prism near-losslessly into state plane
coordinates (the only source of error comes from a set of 4 image coordinates
that cannot be perfectly converted into a rectangle in state plane
coordinates), ii.) shifting the footprint in state plane coordinates along the
z-axis, iii.) re-projecting the resulting points back into the image, iv.)
comparing the reprojected “top points” to the original top of the rectangular
prism in image coordinates, and v.) adjusting the height iteratively to
minimize the re-projection error until convergence.

###  1.3 State Plane $\rightarrow$ Roadway Coordinate System Conversion

Next, we consider the conversion of points in state plane coordinates to
roadway coordinates. In most cases, we care to convert a set of state plane
coordinate points roughly in a rectangular prism (i.e. vehicle 3D bounding
box) into roadway coordinates; thus, we define this conversion for a
rectangular prism. A single point can be converted between state plane
coordinates and roadway coordinates by treating it as a rectangular prism with
zero length, width and height.

Let $\mathbf{O}^{(st)}$ be a 3D bounding box representation in state plane
coordinates, an 8$\times$3 matrix of x,y, and z coordinates for each corner of
the box. (Note that these corners need not exactly correspond to an orthogonal
rectangular prism, but the roadway coordinate equivalent will be exactly
orthogonal so some truncation will occur.) We reference, for example, the back
bottom right (from the perspective of the rear of the vehicle) of object
$\mathbf{O}^{(st)}$ as
$\mathbf{o}_{bbr}^{(st)}=[x^{(st)}_{bbr},y^{(st)}_{bbr},z^{(st)}_{bbr}]$, such
that =
$\mathbf{O}^{(st)}=[\mathbf{o}^{(st)}_{bbl},\mathbf{o}^{(st)}_{bbr},\mathbf{o}^{(st)}_{btl},\mathbf{o}^{(st)}_{btr},\mathbf{o}^{(st)}_{fbl},\mathbf{o}^{(st)}_{fbr},\mathbf{o}^{(st)}_{ftl},\mathbf{o}^{(st)}_{ftr}]$.
(For the single-point case described above, all 8 corner coordinates are
identical).

Next, Let $\mathbf{O}^{(r)}=[x^{(r)}_{o},y^{(r)}_{o},l,w,h]$ be the
corresponding object representation of $\mathbf{O}^{(st)}$ in roadway
coordinates. $x^{(r)}$ and $y^{(r)}$ are the roadway coordinate longitudinal
and lateral coordinates (in feet), and $l$, $w$,and $h$ are the length, width,
and height of the object respectively (in feet).

Let $\mathbf{o}^{(st)_{c}}$ denote the back bottom center coordinate of object
$\mathbf{O}^{(st)}$. By convention, this point is referenced as the primary
position of object $\mathbf{O}^{(st)}$. Let $\mathbf{o}^{(st)_{spl}}$ denote
the point on the center-line spline (i.e. $y^{(r)}=0$) with the same $x^{(r)}$
coordinate as $\mathbf{o}^{(st)}_{c}$.

Let $F$ be the second-order spline parameterizing the roadway center-line in
state plane coordinates. In other words, $F$ defines the longitudinal
curvilinear axis $y^{(r)}=0$ along this spline. $F$ is fit by manually
labeling a sufficiently large number of points along the interior yellow line
for both directions of travel (in state plane coordinates). A spline is fit to
each yellow line, and a third spline is fit to lie precisely halfway between
these two splines. Spline control points are selected at suitably sparse
intervals (200 foot minimum spacing) such that the spline is relatively smooth
while still capturing the roadway curvature (see Appendix 1.6).

Given $\mathbf{O}^{(st)}$, we first obtain $l$,$w$ and $h$ by computing the
average distance between points on the front and back, left and right, or top
and bottom of the vehicle respectively. Next, we obtain
$\mathbf{o}^{(st)}_{c}$ by computing the average $x^{(st)}$ and $y^{(st)}$
state plane coordinates of the 4 back rectangular prism corners.

Next, we solve for $x^{(r)}_{o}$ by solving the following optimization:

$x^{(r)}_{o}=\operatorname*{arg\,min}_{x^{(r)}}||(F(x^{(r)}),\mathbf{o}^{(st)}_{c}||_{2}$
(7)

using L2-norm (Euclidean distance) between the two points in state plane
coordinate space. In other words, determine the point on the roadway spline
closest to the back center of the rectangular prism $\mathbf{o}^{(st)}_{c}$.
This minimizing point is the corresponding roadway longitudinal coordinate
$x^{(r)}_{o}$, and the distance from the minimum distance point is roadway
lateral coordinate $y^{(r)}_{o}$.

$y^{(r)}_{o}=\min_{x^{(r)}}||(F(x^{(r)}),\mathbf{o}^{(st)}_{c}||_{2}$ (8)

Noting that the I-24 MOTION roadway segment has monotonically increasing
$x^{(st)}$ coordinate, a secondary spline $\tilde{G}(x^{(st)})$ is defined to
parameterize $x^{(r)}$ as a function of $x^{(st)}$, which yields a good
initial guess for the closest roadway longitudinal coordinate for a given
point in state plane coordinates. This optimization can then be solved to
arbitrary precision, yielding the complete roadway coordinate for the object
$\mathbf{O}^{(r)}=[x^{(r)}_{o},y^{(r)}_{o},l,w,h]$.

Constant Yellow Line Constraint: It is observed that points along the yellow
line for each roadway direction of travel have non-constant y-position due to
the varying width of the median. It is desirable that the yellow line (and by
extension each set of lane-dividing markings) has constant y position. We
finally apply the following shift:

$y^{(r)}_{o}\mathrel{+}\mathrel{=}(C-\gamma(x^{(r)}_{o}))$ (9)

where $C$ is the desired constant yellow line y-coordinate (in this case -12
for westbound-side coordinates and +12 for eastbound-side coordinates) and
$\gamma(x^{(r)}_{o})$ represents the uncorrected yellow-line coordinate at the
given x-position per roadway side. Note that this creates a discontinuity in
the coordinate system near $y=0$ on each side, but these portions of the
coordinate system are not used for objects on the roadway.

###  1.4 Roadway $\rightarrow$ State Plane Conversion

Reverse Yellow Line Constraint: First, the inverse yellow-line shift must be
applied to the y-coordinate:

$y^{(r)}_{o}\mathrel{-}\mathrel{=}(C-\gamma(x^{(r)}_{o}))$ (10)

Next, given roadway coordinates for an object $\mathbf{O}^{(r)}$, first find
the corresponding point on the roadway center-line spline in state plane
coordinates $\mathbf{o}^{(st)}_{c}$ according to:

$F(x^{(r)}_{o})=\mathbf{o}^{(st)}_{spl}$ (11)

To obtain the back center coordinate $\mathbf{o}^{(st)}_{c}$, we must offset
$\mathbf{o}^{(st)}_{spl}$ by length $y^{(r)}$ in the direction perpendicular
to the roadway centerline spline at $\mathbf{o}^{(st)}_{c}$. Let
$\overrightarrow{\mathbf{u}}_{F}$ be the unit vector in the same direction as
the derivative spline $F^{\prime}$, and let
$\overrightarrow{\mathbf{u}}_{1/F}$ be the unit vector in the perpendicular
direction (along the state plane, i.e. $z^{(st)}=0$. Note that care should be
given to ensure that the positive direction of
$\overrightarrow{\mathbf{u}}_{1/F}$ points towards the eastbound side of the
roadway with positive $y^{(r)}$.) Then, $\mathbf{o}^{(st)}_{c}$ is given by:

$\mathbf{o}^{(st)}_{c}=\mathbf{o}^{(st)}_{spl}+y^{(r)}\cdot\overrightarrow{\mathbf{u}}_{1/F}$
(12)

From here, the corner state plane coordinates for the right and left
coordinates of the rectangular prism can be obtained by offsetting
$\mathbf{o}^{(st)}_{c}$ by $\pm\frac{1}{2}$ times $w$ in the direction of
$\overrightarrow{\mathbf{u}}_{1/F}$, and the front coordinates of the
rectangular prism can similarly be obtained by offsetting by $l$ in the
direction of $\overrightarrow{\mathbf{u}}_{F}$ or in the opposite direction
for objects on the westbound or negative $y^{(r)}$ side of the roadway.
Similarly, the top coordinates can be obtained by offsetting by a factor of
$h$ in the $z^{(st)}$ direction. The direction of travel for an object can be
obtained as the sign of the $y^{(r)}$ coordinate (negative for WB, positive
for EB).For example, for an eastbound object the front top left coordinate can
be obtained as:

$\mathbf{o}^{(st)}_{ftl}=\mathbf{o}^{(st)}_{c}-\frac{1}{2}\cdot
w\cdot\overrightarrow{\mathbf{u}}_{1/F}+l\cdot\overrightarrow{\mathbf{u}}_{F}+h\cdot\mathbf{[0,0,1]}$
(13)

###  1.5 Correspondence Point Labeling

A brief description and a few pictures of overhead and camera fields of view.

Correspondence points are obtained for visible lane marking lines and other
easily recognizable landmarks on the roadway are used as correspondence points
in each camera field of view. Labeling is carried out in a custom OpenCV GUI,
with semi-automated point selection to speed this process over the 230+
cameras utilized in this work. An average frame from each camera is used to
minimize the presence of occluding vehicles during labeling. Figure 1 shows an
example of labeled points and the labeled relevant field of view for a single
camera. The corresponding pixel coordinates for each manually selected point
in each camera field of view, are stored, for each direction of travel on the
roadway, each with a unique identifier. Typically, at least 10 roadway
markings or 40 correspondence points are labeled per roadway direction of
travel. Additionally, note that this initial correspondence point selection
was performed for a day a few weeks prior to the day on which video is
recorded in this work.

Figure 1: Correspondence points (blue) and fields of view (shaded polygons)
for each roadway direction of travel, labeled in a single camera field of view
(P27C01). Around 100 correspondence points are visible for each side of the
roadway. Figure 2: Correspondence points (green for WB red for EB) labeled in
camera imagery for all cameras on one pole (P17). Fields of view (blue for WB
and orange for EB) and a mask denoting relevant portions of the image for
tracking (dotted red polygon) are also shown. Similar plots are included in a
separate file for all cameras.

Each correspondence point is also labeled in  global information system (GIS)
software, giving the precise GPS / state-plane coordinate system coordinates
for each labeled corresponding point. The same unique identifier is assigned
to each corresponding point previously labeled within a camera image.
Additionally, points along each yellow line were labeled periodically (every
$\sim$50 feet). Figure 3 shows an example of labeled points in the state plane
aerial imagery.

Figure 3: Correspondence points labeled in aerial imagery using the GIS tool
for a very small subset ($\sim$500 feet) of roadway, for one direction of
travel only.

###  1.6 Centerline Spline Fitting

The following procedure was used to fit the center-line spline for the roadway
coordinate system ($y^{(r)}=0$):

  1. 1.

For each roadway direction of travel, fit a spline parameterizing the
$(x^{(st)},y^{(st)})$ points as a function of $u$. Let $F_{(WB)}(u)$ and
$F_{(EB)}(u)$ denote these splines. All splines implemented in this work are
based on SciPy’s spline package [2020SciPy-NMeth] and are constrained to have
a minimum control point spacing of 200 feet to ensure smoothness.

  2. 2.

Moving along the EB yellow line spline $F_{(EB)}(u)$, sample points at fine (1
foot) intervals. For each sampled point $u$, find the closest point on
$F_{(WB)}(u)$ $u^{\prime}$. Store the midpoint $u*$.

  3. 3.

Fit another spline $F_{med}(u)$ to the set of midpoints from the previous
step.

  4. 4.

Sample $F_{med}(u)$ at fine (0.1 foot intervals). Compute the distance between
each consecutive pair, and compute the cumulative distance along the spline to
each point via finite difference approximation.

  5. 5.

Re-parameterize $F_{med}(u)$ such that each point on the midpoint corresponds
to the integrated distance along the spline to that point. This yields the
final spline $F(x^{(r)})$ used above.

Lastly, $\gamma(x^{(r)}))$ must be computed for each roadway direction of
travel. This is done by sampling $F(x^{(r)})$ for each direction at regular
(5-foot) intervals and recording the distance to the closest point on each of
$F_{(WB)}$ and $F_{(EB)}$. Since these offsets change very little ($\sim$1
foot per mile), the offsets are stored in a lookup table rather than fitting a
true offset function. Thus, $\gamma(x^{(r)}))$ returns the recorded offset
from $F(x^{(r^{\prime})})$ to $F_{(WB)}$ or $F_{(EB)}$ for the closest sampled
value $x^{(r^{\prime})}$.

Figure 4(a) shows the labeled points from aerial imagery (yellow and black
lines) and the points labeled in each image (blue dots), transformed into
roadway coordinates. Figure 4(b) shows the largest magnitude shift in labeled
correspondence points, when converted to the roadway coordinate system. Such
“discontinuities” are due to slight misalignments betweeen aerial images taken
during different passes of the photographing aircraft, and result in small
(less than 2 foot) errors in all cases. While these misalignments could be
corrected for by applying a smoothing to the correspondence points labeled
within the aerial imagery, this was ultimately decided against because the
smoothing would need to be performed with respect to, essentially, the roadway
coordinate system spline itself, which is in turn a product of the labeled
aerial imagery points. This would produce complex and less well-understood
artifacts in the resulting final coordinate system.

(a)

(b)

Figure 4: (a) correspondence points labeled in aerial imagery (yellow and
black) and in individual camera fields of view (blue) are projected into the
roadway coordinate system. (b) A close-up detailing a type of artifact visible
in the coordinate system as a result of misalignment in the aerial imagery.

##  2 Homography Re-estimation Methods

In ideal circumstances homographies can be estimated once (as discussed in the
previous section), and used continuously until some drastic change in the
system, (e.g. the roadway section is rebuilt or the lane markings are
repainted). Aside from these rare events there are several other factors which
significantly degrade accuracy; long term events can be the replacement of the
camera, inaccuracies in the pan-tilt mechanism during homing, settlement of
the foundation or the seasonal temperature change. Although these can be dealt
with occasional re-calibration, the sunflower effect requires constant re-
estimation of the homographies. This effect, the tilting of the metal
infrastructure poles due to the differential heating of the sun and shade-
facing sides, can cause significant homography errors (sometimes greater than
10 feet, see Figure 5 for a typical camera example), both in timescales of
hours for the daily warm up and cool down cycle, and also in the timescales of
minutes, caused by the varying cloud coverage (see Figure 6, showing how the
positional drift in homographies varies over just 12.5 minutes). A video is
included in supplementary material showing the magnitude of drift over the
course of a day. The camera movement can easily be seen by viewing the
relation between the lane markers and the ROI rectangle (which maintains
constant pixel coordinates throughout the video). For reference, a typical
dash line is 10 feet long.

(a)

(b)

(c)

Figure 5: Homography goodness of fit and uncorrected drift between reference
homography and true scene homography according to two metrics (Sub Drift and
Full Drift) for a single typical camera. Horizontal grey dashes show the two
time instances for Figure 6.

(a)

(b)

Figure 6: Full drift for the same camera at two time instances. For a given
error value $X$, this indicates that the original reference homography and the
true homography for the scene map the same image coordinate to points in the
state plane $X$ feet apart. The westbound field of view (i.e. the most
important portion of the image) is shown as a red polygon.

Note, that tilting not only degrades single-camera accuracy, but due to the
multiple pole architecture and the camera FOVs, homography errors cause
significant misalignment in state plane coordinates between poles, resulting
in vehicle tracking fragmentations . (Cameras on the same pole are less
susceptible because they are rigidly mounted and roughly ”moving” together.)
Thus, the re-estimation of the homography is necessary on relatively short
timescales, or at a minimum on a daily basis.

Since manual re-annotation is not possible for continuous operation, a feature
point re-identification method is developed, along with re-estimation methods
for a static daily and for a dynamic, time-dependent homography.

All methods and filtering steps are based on the assumption that the image is
only slightly distorted from the manually labeled ground truth. This
assumption is necessary, because the homography and the re-discovery relies on
the lane markings which are non-unique and repetitive, thus a significant
shift can cause a mis-identification and misalignment. All presented methods
are offline (they utilize both past and future information relative to the
specified time instance (i.e. use a non-casual filter)).

Automated homography generation is divided into two distinct steps. First the
feature points are re-discovered and a new homography is calculated for a
single time instance. Second, a refinement step operates on the homographies
themselves of which removes outliers, filters, and aggregates over time. This
approach provides robustness for the final estimates based on the assumption
that the re-estimation provides good results in general but are prone to
outliers and errors. To comply to this idea the implemented refinement
processes are draconian, and tuned to minimize false positives at the expense
of occasionally removing some true positives. The ”good in general” assumption
is empirically tested, by rigorously inspecting the results and fine tuning
the filter parameters of each steps. In case of faulty estimates manual
inspection always shows problems which could cause confusion even for human
observers or significant noise and drop in image quality (e.g.: at dawn).

###  2.1 Notation

The following notation convention is used throughout the remainder of this
Section. Also see Table 2 for a list of utilized symbols and their meanings.

  * •

Mathcal script (e.g. $\mathcal{I}$) is used for sets of points. Bold Text
(e.g. H) is used for vectors and matrices.

  * •

The subscript $t$ denotes a set of points or a homography for a specific time
instance (e.g. $\mathcal{I}_{t}$)

  * •

No subscript denotes the initial reference points or homography (e.g.
$\mathcal{I}$, $\mathcal{A}$, and H)

  * •

The prime symbol denotes rediscovered points (e.g. $\mathcal{I}^{\prime}_{t}$,
which corresponds to subset $\mathcal{A}_{t}$) , or homographies based on
rediscovered points (e.g.: $\textbf{H}_{t}^{\prime}$)

  * •

a superscript s denotes an image-to-image homography based on SIFT-FLANN
matching (e.g. $\textbf{H}_{t}^{s}$)

  * •

A superscript asterisk symbol $\ast$ denotes a homography estimate calculated
with methods presented in the previous appendix (e.g.: $\textbf{H}^{\ast}$).

  * •

An arrow $\mathcal{I}\xrightarrow{\textbf{H}}$ indicates a linear
transformation on the point set $\mathcal{I}$ according to homography H. (e.g.
$\mathcal{I}_{t}^{\prime}\xrightarrow[]{\textbf{H}_{t}^{\prime}}$
$\mathrel{\stackon[1.5pt]{=}{\stretchto{\scalerel*[\widthof{=}]{\wedge}{\rule{4.30554pt}{12.91663pt}}}{0.5ex}}}$
$\mathcal{A}_{t}$ denotes that the set $\mathcal{I}_{t}^{\prime}$ of
rediscovered points at time $t$, projected from image coordinates to state
plane coordinates by $\textbf{H}_{t}^{\prime}$, estimates the set of
corresponding points labeled in state plane coordinates $\mathcal{A}_{t}$).

Symbol | Definition
---|---
$\mathcal{I}$ | labeled image points (in image coordinates, unit: pixels)
$\mathcal{I}_{t}$ | subset of the labeled image points at time $t$, corresponding to the successfully rediscovered points ($\mathcal{I}_{t}\subset\mathcal{I}$)
$\mathcal{I}_{t}^{\prime}$ | rediscovered image points ($|\mathcal{I}_{t}^{\prime}|\leq|\mathcal{I}|$)
$\mathcal{I}_{t}^{s}$ | labeled image points after perspective transformation, utilizing the SIFT-FLANN matcher, such that $\mathcal{I}\xrightarrow{\textbf{H}_{t}^{s}}\mathcal{I}_{t}^{s}$
$\mathcal{I}_{t}^{\ast}$ | ”labeled” image points generated such as $A\xrightarrow{\textbf{H}_{t}^{\ast\neg}}\mathcal{I}_{t}^{\ast}$
$\mathcal{A}$ | aerial points (in state plane coordinates, unit: feet)
$\mathcal{A}_{t}$ | subset of the aerial points at time $t$, corresponding to the successfuly rediscovered points ($\mathcal{A}_{t}\subset\mathcal{A}$)
H | homography such that $\mathcal{I}\xrightarrow{\textbf{H}}\mathcal{A}$
$\textbf{H}^{\neg}$ | inverse homography such that $\mathcal{A}\xrightarrow{\textbf{H}^{\neg}}\mathcal{I}$
$\textbf{H}_{t}^{\prime}$ | homography such that $\mathcal{I}_{t}^{\prime}\xrightarrow{\textbf{H}_{t}^{\prime}}\mathcal{A}_{t}$
$\textbf{H}^{\ast}$ | new static homography estimate by method from Appendix 2.3
$\textbf{H}_{t}^{\ast}$ | new dynamic homography estimate for time $t$ by methods from Sections 2.4-2.5
$\textbf{H}_{t}^{s}$ | projection (image to image) generated with the SIFT-FLANN matcher results
Table 2: Notations for various point sets and homographies utilized in this
work.

###  2.2 Homography Generation for a Single Time Instance

This section briefly describes the steps necessary for the homography
generation for a single time instance $t$.

  1. 1.

Generate a background extracted image to remove the vehicles from the scene.
This is implemented as a 1 minute long averaging of frames, with 50% overlap
in time. Note that in case of heavy traffic (essentially for stopped vehicles)
this time period might not be sufficient for complete removal.

  2. 2.

Compute an initial alignment estimate between the frame and the reference
frame, on which the original annotations were made. This step utilize the SIFT
algorithm to find feature points and a FLANN based matcher (both implemented
in OpenCV). Note that the detected features can be anywhere on the image, not
necessarily corresponding to semantically meaningful points (e.g. grass,
trees, parking cars outside the road), and most importantly generally not lie
on the plane of the roadway.

  3. 3.

Based on the SIFT-FLANN generated corresponding points generate an image-to-
image homography ($\textbf{H}_{t}^{s}$).

  4. 4.

Examine $\textbf{H}_{t}^{s}$: transformation (e.g.: translation, rotation,
scale) should be minimal, based on the slight distortion assumption. Otherwise
discard the current time instant, because the calculated alignment is likely
from erroneously matched points. No further steps are possible.

  5. 5.

Binarize the averaged frame via OpenCV’s $threshold$ function
[bradski2000opencv] and run a contour detection algorithm. This provides good
results because lane markings are high contrast, distinct features.

  6. 6.

Transform the originally labeled feature points $\mathcal{I}$, i.e. the four
corner points of the lane markers, with the calculated $\textbf{H}_{t}^{s}$
(producing $\mathcal{I}^{s}_{t}$) and calculate the geometric center for each
marker. This provides a seed point for rediscovery.

  7. 7.

Select contour areas which include a seed point, thus creating a set of
candidate for lane marks.

  8. 8.

For each marker, select the 4 corner points on the corresponding contour that
produce the largest area, thus creating a quadrilateral.

  9. 9.

Filter the candidates by contour area, contour area to quadrilateral area
ratio, ratio of area compared to the labeled reference area. This step is
important to remove candidates which are merged, partially occluded, or
otherwise not rediscovered correctly.

  10. 10.

Select the proper labels and store the 4 corner points as feature points of
each remaining markers. With this step the rediscovery steps are complete,
yielding $\mathcal{I}^{\prime}_{t}$.

  11. 11.

Run a homography matrix estimator (which utilizes RANSAC) between the
rediscovered $\mathcal{I}^{\prime}_{t}$ and the corresponding aerial
$\mathcal{A}_{t}$ points. Remove the feature points which are considered
outliers by the algorithm. In the current implementation these are points of
which are further than 2̃ feet on the state plane. The resulting homography
for time instance $t$ is written as $\textbf{H}_{t}^{\prime}$.

Figure 7: An example of SIFT-FLANN based seed points from step 6 (red dots),
detected contours (purple outlines) and re-detected dash corner points (yellow
dots) for a 1-minute average frame. Because traffic at this time instance was
partially stopped, some vehicle artifacts are visible in the averaged frame. A
few dashes are not successfully rediscovered.

The result of the listed process is a set of rediscovered feature points
($\mathcal{I}^{\prime}_{t}$) and a new homography ($\textbf{H}_{t}^{\prime}$)
for a given time instant $t$. Note that although the steps involve heavy
filtering and outlier detection the resulting, a standalone homography should
not be used without further processing for the following reasons:

  * •

The number and positions of the rediscovered lane markings are not guaranteed,
thus it is possible to only rediscover lane markings which lie on the same
line, thus providing poor estimates for the perpendicular axis.

  * •

Lane markings ”far” from the camera are hard to properly detect due to their
small area and slight camera movements can cause misidentification.

  * •

The RANSAC based homography estimator does not guarantee that the selected
points are proper. In extreme cases it can select an erroneous subset and fit
a good homography to it. As discussed previously, the lane markings are non-
unique and repetitive: usually there are multiple subsets of points with
nearly identical relative arrangement.

###  2.3 Static Homography Generation

To counter for the long term homography errors and somewhat compensate for the
sunflower effect a single homography can be generated for a given time
interval, e.g.: the length of the recording or for the desired time range.
This method has the advantage of time-invariance (simplicity) and can utilize
all the instantaneous estimates over the specified window, resulting in ample
redundancy for outlier removal. The method includes two steps:

  1. 1.

In the preparation phase the outliers are removed, i.e. the homography
estimates of which significantly differ from the others. The current
implementations considers a homography an outlier if any of its component (in
the 3x3 matrix) deviates more than 30% from the arithmetic mean of all
estimates. Note that this process is iterative.

  2. 2.

The main step is the averaging, i.e. calculating the arithmetic mean for each
matrix component to produce $\textbf{H}^{\ast}$.

###  2.4 Dynamic Homography Generation

To counter for both long and short term homography problems a dynamic, time
dependent method is proposed. This solution is capable of modeling short-term
camera motion (e.g. from the sunflower effect), but is less trivial implement
and use compared to the static version because it requires modifying the
homography as a function of time.

The underlying idea of this method is the ”smoothing” of the original
estimates, based on the temporally nearby values. Note that this is not
straightforward because estimates are not necessarily available at all time
instances, e.g.: in case of heavy traffic the occlusion prevents homography
estimation, similarly at the beginning and at the end of the recording there
is no preceding and subsequent information available. (This is especially
problematic at dawn and dusk where the even the recorded images are noisy and
blurred thus making feature rediscovery hard or even impossible.) Thus a
method is necessary of which is capable to adapt to these situations. An
important and necessary property of a good dynamic homography is that it
provides ”smooth” variability over time, preventing discontinuities in
positional data (which could in turn cause fragmentations during object
tracking.)

The proposed solution includes three steps:

  1. 1.

In the preparation phase the outliers are removed. This step is exactly the
same as in the static case.

  2. 2.

A ”window size” parameter is computed, based on the available number of
homography estimates in a specified temporal neighbourhood around the given
time instance. This step ensures that in case of missing estimates the window
size is larger, thus accommodating more data points.

  3. 3.

The final step is the ”smoothing” of each homography matrix component, where a
Gaussian kernel function is used over time. The shape of the function, or more
precisely the variance is determined by the previous step.

Note that this method can provide a homography estimate
$\textbf{H}_{t}^{\ast}$ at any time instance, not just when instantaneous
estimates are available. In our implementation a new estimate is generated at
10 second intervals.

###  2.5 SIFT-FLANN-based homography generation (existing)

No lane marking rediscovery is performed. Instead, the set of shifted points
($\mathcal{I}_{t}^{s}$) generated from the original labels ($\mathcal{I}$) is
direcly used from Appendix 2.2 step 6 and above to fit a homography estimate
$\textbf{H}_{t}^{\ast}$ for time $t$ i.e.
$\mathcal{I}\xrightarrow[]{\textbf{H}_{t}^{s}}\mathcal{I}_{t}^{s}\xrightarrow[]{\textbf{H}_{t}^{\ast}}\mathcal{A}$.

##  3 Homography Error Metrics

Next, we turn to defining metrics for assessing the effectiveness and accuracy
of the proposed point rediscovery and homography re-estimation metrics. This
appendix provides an overview of the possible metrics and their properties.

###  3.1 Notation

In addition to the notation defined in the previous section, the comparison
operator $E$ is defined as a function of two sets of points
($\mathcal{B},\mathcal{C}\subset\mathbb{R}^{2}$). The operator calculates the
point-wise L2-norm between the common subset of points shared by $\mathcal{B}$
and $\mathcal{C}$ , (i.e. the distances between two points which correspond to
the same feature), and produces a set of values. The resulting set can be used
to calculate aggregate metrics, such as mean, maximum, and standard deviation
of error.

$E(\mathcal{B},\mathcal{C})=\\{||b_{i},c_{i}||_{2},...,||b_{n},c_{n}||_{2}\\}\>\forall
i\text{\> s.t.
\>}b_{i}\in\mathcal{B},c_{i}\in\mathcal{C},b_{i}\mathrel{\stackon[1.5pt]{=}{\stretchto{\scalerel*[\widthof{=}]{\wedge}{\rule{4.30554pt}{12.91663pt}}}{0.5ex}}}c_{i}$
(14)

###  3.2 Reference Metrics

Metrics listed in this appendix are computed by comparing the original
reference labels ($\mathcal{I}$ and $\mathcal{I}$) to rediscovered points
($\mathcal{I}_{t}^{\prime}$). In some cases the same set of points is
transformed with the original reference homography (H or $\textbf{H}^{\neg}$)
and the instantaneous homography estimate $\textbf{H}_{t}^{\prime}$. They
measure i.) the quality of the rediscovery and instantaneous homography re-
estimation process, and ii.) the movement of the cameras.

Most metrics presented here can be calculated in both image and state plane
coordinate systems; thus we present them in pairs. The resulting error values
have units of pixels or distance (feet), respectively. In practical
evaluations we utilize the state plane distance (units of feet) as it is
invariant to the image resolution and FOV, thus comparable across cameras. As
in Appendix 1, we use a superscript $(im)$ or $(st)$ to denote the coordinate
system for each metric.

  * •

Sub Drift I: This metric compares a subset of the originally labeled image
points to the rediscovered points. The resulting distances can be interpreted
as the drift caused by the camera movement. Note that in case of perfect
alignment the maximum error is zero.

$\text{SubDriftI}^{(im)}=E(\mathcal{I}_{t},\mathcal{I}_{t}^{\prime})$ (15)
$\boxed{\text{SubDriftI}^{(st)}=E(\mathcal{I}_{t}\xrightarrow{\textbf{H}},\mathcal{I}_{t}^{\prime}\xrightarrow{\textbf{H}})}$
(16)

  * •

Sub Drift II (not used): This metric compares a subset of the originally
labeled aerial points to the rediscovered points. The resulting distances are
a combined error of the drift caused by the camera movement and the homography
fitness error. Thus even in case of no movement the error is non-zero, as the
homography is not a perfectly fit mapping of $\mathcal{I}$ to $\mathcal{A}$.

$\text{SubDriftII}^{(im)}=E(\mathcal{A}_{t}^{\prime}\xrightarrow{\textbf{H}^{\neg}},\mathcal{I}_{t}^{\prime})$
(17)
$\text{SubDriftII}^{(st)}=E(\mathcal{I}_{t}^{\prime}\xrightarrow{\textbf{H}},\mathcal{A}_{t})$
(18)

  * •

Full Drift: This metric compares the full set of the originally labeled image
or aerial points transformed by the reference and the re-estimated homography.
The resulting distances can be interpreted as drift caused by the camera
movement. In the case of no camera movement, the error is zero (except for
possible stochasticity in the optimal homography $H^{\prime}_{t}$ selected by
the RANSAC algorithm utilized.)

$\text{FullDrift}^{(im)}=E(\mathcal{A}\xrightarrow{\textbf{H}^{\neg}},\mathcal{A}\xrightarrow{\textbf{H}_{t}^{\prime\neg}})$
(19)
$\boxed{\text{FullDrift}^{(st)}=E(\mathcal{I}\xrightarrow{\textbf{H}},\mathcal{I}\xrightarrow{\textbf{H}_{t}^{\prime}})}$
(20)

  * •

Sub Fitness:  This metric shows the fitness of the estimated homography,
utilizing the rediscovered points and a matching subset of aerial points. This
error is caused by a combination of e.g.: inaccurate feature point
rediscovery, lens distortion, non flat roadway.

$\text{Fitness}^{(im)}=E(\mathcal{A}_{t}\xrightarrow{\textbf{H}_{t}^{\prime\neg}},\mathcal{I}_{t}^{\prime})$
(21)
$\boxed{\text{Fitness}^{(st)}=E(\mathcal{I}_{t}^{\prime}\xrightarrow{\textbf{H}_{t}^{\prime}},\mathcal{A}_{t})}$
(22)

###  3.3 Homography Re-estimation Metrics

Metrics listed here are derived from the metrics introduced at the previous
section. Substituting the reference labels and homographies with the re-
estimated homographies and points sets derived from them. The interpretation
and purpose of these metrics are to qualify how good are the new estimations;
i.e. a perfect homography provides a perfect alignment to the detections.

  * •

Sub Drift: This metric compares a subset of points generated with the re-
estimated homography (from the originally labeled points) to the rediscovered
points. The resulting distances can be interpreted as a combined error of the
re-estimation and the homography fitness error. Therefore even in case of a
perfectly aligned homography the error is non-zero. Note that, because the new
subset of points are generated from the reference labeled points this metric
essentially merge Sub Drift I and II presented in the previous section.

$\text{SubDrift}^{(im)}=E(\mathcal{I}_{t}^{\ast},\mathcal{I}_{t}^{\prime})\Longrightarrow
E(\mathcal{A}_{t}\xrightarrow{\textbf{H}_{t}^{\ast\neg}},\mathcal{I}_{t}^{\prime})$
(23)
$\text{SubDrift}^{(st)}=E(\mathcal{I}_{t}^{\ast}\xrightarrow{\textbf{H}_{t}^{\ast}},\mathcal{I}_{t}^{\prime}\xrightarrow{\textbf{H}_{t}^{\ast}})\Longrightarrow
E(\mathcal{A}\xrightarrow{\textbf{H}_{t}^{\ast\neg}}\xrightarrow{\textbf{H}_{t}^{\ast}},\mathcal{I}_{t}^{\prime}\xrightarrow{\textbf{H}_{t}^{\ast}})\Longrightarrow
E(\mathcal{A}_{t},\mathcal{I}_{t}^{\prime}\xrightarrow{\textbf{H}_{t}^{\ast}})$
(24)

  * •

Full Drift: This metric compares the full set of the originally labeled image
or aerial points transformed by the the re-estimated and by the instantaneous
homography. In effect, these original labeled points are chosen as “proxy
points” for the real locations of the correspondence points in the image as
they are a full set known to lie near the true locations. The resulting
distances can be interpreted as the error of the homography re-estimation
process. Note that in case of a perfect alignment the error is zero.

$\text{FullDrift}^{(im)}=E(\mathcal{A}\xrightarrow{\textbf{H}_{t}^{\ast\neg}},\mathcal{A}\xrightarrow{\textbf{H}_{t}^{\prime\neg}})$
(25)
$\boxed{\text{FullDrift}^{(st)}=E(\mathcal{I}\xrightarrow{\textbf{H}_{t}^{\ast}},\mathcal{I}\xrightarrow{\textbf{H}_{t}^{\prime}})}$
(26)

  * •

Fitness:  Since the fitness metrics does not depend on the re-estimated
homography, they are equivalent to the ones discussed in the previous section.

The homography re-estimation performance can be measured by the Full and Sub
Drift metrics, but both have caveats to consider: Sub Drift is a better
measure of the  real magnitude of the error, because it incorporates
inaccuracies of both the homography and the homography fitness, e.g. errors
caused by intrinsic camera problems and the non-flat roadway. On the other
hand the Sub Drift metric only includes feature points which are re-detected,
thus hard to detect points further from the pole are often missing from the
set. This is crucial because the perspective those ”far” points account for
larger state plane errors than the more easily detectable closer ones. The
combined effect of these that the Sub Drift is better in the assessment of the
expected minimal error, this is not the case for the maximum. Also worth
noting is that Sub Drift does not explictly rely on $H^{\prime}_{t}$ being
well-fit, so produces a reliable estimate of error even when the fitting of
$H^{\prime}_{t}$ fails (though in practice this usually indicates poor point
rediscovery).

The Full Drift metric is useful to assess the performance of the homography
re-estimation itself, because in case of perfect alignment the error is
reduced to zero. In opposite to the Sub Drift this, metric includes all
labeled points, some of which are actually outside the processed ROI, thus
resulting in an apparently larger maximum error than during vehicle detection.

Based on these considerations, we select Full Drift (Equation 20) as our
primary metric and use it to report the $error$ for each method utilized as
equation 2 the main text. We also utilize Sub Fitness (Equation 22) as a
baseline measure for the error floor of each fit homography ($fitness$ in
equation (1) of the main text).

###  3.4 Additional Homography Metric Figures

Figure 8 shows the different state-plane error metrics for a camera with a
long field of view. P05C06 is chosen as it clearly illustrates both long and
short-term drift in the homography, as well as the relations between the
different metrics: On Figure 8(a) the homography Fittness is shown,
representing the minimum achievable accuracy for the vehicle detection. On
Figure 8(b) the SubDriftI error is slightly higher than FullDrift shown on
Figure 8(c), because it contains both the fitness and homography error. In
this case, the maximum SubDriftI values are higher than the maximum FullDrift
values because the rediscovery rate for the feature points is $\sim$80% for
this particular case (relatively high), thus distant points are also likely
included. If the drift were more extreme, the most distant correspondence
points might not be successfully rediscovered, and the overall reported
SubDriftI error may be driven down as a result. This dependence on the number
of rediscovered points is one reason why FullDrift is preferred to SubDrift in
this work for comparing homography re-estimation performance.

(a)

(b)

(c)

Figure 8: (Repeated from Section 2.) Goodness of fit and uncorrected drift
between reference homography and true scene homography according to two
metrics (Sub Drift and Full Drift) for a single typical camera. Horizontal
grey dashes show the two time instances for Figure 6.

Figure 9, 10 shows the mean FullDrift and SubDrift calculated and averaged
over the time of the recording, for all cameras in the system, sorted by
magnitude, for each direction of roadway travel. Note that not all cameras
have fields of view defined for both sides of the roadway (hence there are
more total westbound fields of view than eastbound fields of view). The mean
SubDrift is Lower than the mean FullDrift over all cameras, and the maximum
mean SubDrift (over all correspondence points and times) is lower than for
FullDrift. However as seen in Figure 8(a), for some time instances the maximum
SubDrift can be higher than the maximum FullDrift.

Figure 9: Mean averaged SubDrift and FullDrift calculated for all westbound
cameras in the system, sorted by error magnitude.

(a)

(b)

Figure 10: Mean averaged SubDrift and FullDrift calculated for all eastbound
cameras in the system, sorted by error magnitude.

The effectiveness of the homography re-estimation methods is illustrated on
Figure 11(b), which shows the additional error above the error floor,
utilizing the FullDrift metric, for two cameras. The reference (blue)
indicates the resulting error without any mitigation, showing both long term
(high mean) and short term (high variance) error. The SIFT-FLANN method
(orange) illustrates the performance of an optical flow based ”camera
stabilization” solution. The static, all-day average method (green) removes
the long term error, although it is mostly unable to remove the error caused
by the sunflower effect. Finally since the dynamic homography (red) utilizes
nearby (temporal) homography estimations for a given time instance, it can
cope with short-term fluctuations. Note that in most cases the all-day average
is superior to the SIFT-FLANN method, although there can be time instances
where the former can produce better results. Figure 11(b) also illustrates,
that long field of view cameras (e.g. P05C06) are more sensitive to camera
movement, compared to mainly downwards-facing cameras (e.g. P05C04), thus the
error caused by drift is more pronounced for them.

(a)

(b)

Figure 11: Error dynamics for two cameras over time with each homography re-
estimation methods. Long field of view cameras e.g. 11(a)

have more pronounced errors over down-looking cameras, e.g. 11(b).

Lastly, figure 12 compares the remaining mean average SubDrift and FullDrift
errors for each camera after (black) SIFT-FLANN feature-matching, (orange)
one-day best fit homography re-estimation, and (red) dynamic homography re-
estimation methods relative to original reference homography baseline (blue).
The overall trends and relative performance amongst methods is unchanged; in
most cases the SIFT-FLANN baseline outperforms the uncorrected reference
homography, the all-day average outperforms SIFT-FLANN, and the dynamic method
outperforms the all-day average.

(a) SubDrift

(b) FullDrift

Figure 12: Remaining SubDrift and FullDrift errors for each camera after
(black) SIFT-FLANN feature-matching, (orange) one-day best fit homography re-
estimation, and (red) dynamic homography re-estimation methods relative to
orignal reference homography baseline (blue). Cameras are grouped by position
on pole (see Figure 2 in main text.) and by side of roadway (westbound
homographies on top, eastbound on bottom).)

##  4 Additional GPS Trajectory Plots

Figure 13 shows plots for additional GPS tracks through the video scene.
Manually annotated object positions (circles) and GPS positions (lines) for
are shown for the overall scene, divided by westbound (top left) and eastbound
(bottom left) direction of travel. (right) x-position relative to the
corrected GPS track (top) and absolute y-position (bottom) are plotted for the
highlighted GPS track. It can be seen that a.) the corrected GPS trajectories
align more closely with object detections (black dots) and b.) there is
significant deviation between the corrected and uncorrected GPS trajectories.
Note that especially the Y-coordinate error in the uncorrected trajectories
varies in character across different GPS trajectories.

Figure 13: Plots for 4 individual GPS trajectories. In each sub-figure, (top
left) shows westbound X-position and (bottom left) shows eastbound X-position
for the whole scene duration. (top right) shows nearby (in relative
X-coordinates) and (bottom right) shows nearby (in Y-coordinates) detections,
manual annotations, and uncorrected GPS trajectory, relative to the corrected
GPS trajectory.

##  5 Object Detector

###  5.1 3D Bounding Box Parameterization

Object detections are produced by a Retinanet object detector with Resnet-50
FPN backbone [lin2017focal]. We extend the first convolutional layer of this
network to allow two input frames (the current and a rolling average frame).
This network has two output heads: a classification head and a regression
head, each of which outputs a set of outputs per anchor box. The shape of the
classification head output $C$ is $[n,c]$ where $n$ is the number of anchor
boxes in the network and $c$ is the number of classes. For an anchor box $i$,
the final output object class and object confidence are taken as:

$\text{class}_{i}=\operatorname*{arg\,max}_{c}(C[i,:])$ (27)
$\text{confidence}_{i}=\max(C[i,:])$ (28)

The shape of the regression head output $R$ is $[n,8]$. Let $A_{i}$
parameterize a single anchor box $i$ according to:

$A[i,:]=[x_{a},y_{a},w_{a},h_{a}]$ (29)

denoting the x-coordinate,y-coordinate,width and height of the anchor box,
respectively. The regression output corresponding to that anchor box $R[i,:]$
parameterizes a rectangular prism (ignoring the effects of perspective
foreshortening) according to:

$R[i,:]=[x_{c},y_{c},x_{l},y_{l},x_{w},y_{w},x_{h},y_{h}]$ (30)

Figure 14: Bounding box rectangular prism parameterization. Anchor box
(purple) reference coordinate system denoted as (x,y) and global frame (blue)
reference coordinate system denoted as (X,Y). Angle of $h$ is exaggerated to
show sub-components $x_{h}$ and $y_{h}$.

where ($x_{c}$,$y_{c}$) indicates the center of the rectangular prism relative
to the anchor box top left corner, $x_{l}$ and $y_{l}$ indicate the x-pixel
and y-pixel components of the rectangular prism’s length, $x_{w}$ and $y_{w}$
indicate the x-pixel and y-pixel components of the rectangular prism’s width,
and $x_{h}$ and $y_{h}$ indicate the x-pixel and y-pixel components of the
rectangular prism’s height. All components (e.g. $x_{l}$,$y_{l}$) are relative
to the anchor box dimensions (in this case $w_{a}$ and $h_{a}$). By
convention, each dimension is directional and is measured relative to the rear
bottom left corner. Figure 14 provides a visual overview of each component.

The (X,Y) pixel coordinates for back bottom left coordinate ($x_{rbl},y_{rbl}$
within the overall image can then be written as:

$\displaystyle x_{bbl}=x_{a}+(x_{c}-x_{l}/2-x_{w}/2+x_{h}/2)*w_{a}$ (31)
$\displaystyle y_{bbl}=y_{a}+(y_{c}-y_{l}/2-y_{w}/2+y_{h}/2)*h_{a}$

The 7 other corner coordinates can be similarly written. The resulting
rectangular prism can then be converted into state plane and roadway
coordinates via the methods discussed in Appendix 1. Readers familiar with 3D
detection in the autonomous vehicle context will likely question the use of 3D
bounding boxes produced explicitly in state space, as in e.g. [hu2019joint,
roddick2018orthographic]. In these contexts, a single camera is statically
mapped to the world coordinate plane (i.e. the same pixel within an image
always corresponds to the same point on the ground plane in world space). This
is not the case in our application, in which a single CNN is desired to
produce object detections for frames from all 234 cameras. Thus, a viewpoint
agnostic anchor box / bounding box formulation such as the one used here is
required. We also consider the direct regression of 3D bounding box corner
coordinates as in [li2020rtm3d] but empirically this results in poor
performance.

###  5.2 Training

We train the above object detector on the I24-3D dataset
[gloudemans2023dataset]. We use the following loss formulation:

$Loss=Loss_{2D}+\beta\times Loss_{3D}+Loss_{cls}+\gamma Loss_{vp}$ (32)

where:

  * •

$Loss_{2D}$ is intersection-over-union loss [rezatofighi2019generalized] for
the minimum enclosing 2D bounding boxes for each of the predicted and target
bounding box

  * •

$Loss_{3D}$ is MSE loss computed between the target and predicted bounding box
corners

  * •

$Loss_{cls}$ is classification focal loss [lin2017focal]

  * •

$\beta$ is a weighting coefficient, here set to 2.

  * •

$\gamma$ is a weighting coefficient, here set to 1/3

and $Loss_{vp}$ is a term designed to enforce the length, width, and height
components of the predicted rectangular prism to align closely to the
corresponding vanishing points by penalizing the angle between these vectors
and the vanishing point directions (relative to the center of the predicted
bounding box). Succinctly, $Loss_{vp}$ is 0 when each bounding box is
perfectly axis-aligned relative to the vanishing points, and is 1 when the
axes of the bounding box are perfectly orthogonal to their respective
vanishing points. In this way, the training enforces the CNN to utilize visual
cues from the image and align bounding boxes along these visual axes.

$\displaystyle Loss_{vp}=$ $\displaystyle
1/2-\alpha_{l}((y_{l}(y_{vpl}-y_{c}))+x_{l}(x_{vpl}-x_{c}))/(2(x_{vpl}-x_{c})(y_{vpl}-y_{c})y_{l}x_{l})\>+$
(33) $\displaystyle
1/2-\alpha_{w}((y_{w}(y_{vpw}-y_{c}))+x_{w}(x_{vpw}-x_{c}))/(2(x_{vpw}-x_{c})(y_{vpw}-y_{c})y_{w}x_{w})\>+$
$\displaystyle
1/2-\alpha_{h}((y_{h}(y_{vph}-y_{c}))+x_{h}(x_{vph}-x_{c}))/(2(x_{vph}-x_{c})(y_{vph}-y_{c})y_{h}x_{h})$

where ($x_{vpl}$,$y_{vpl}$) is the vanishing point, in pixel coordinates,
aligned with the roadway direction of travel, ($x_{vpw}$,$y_{vpw}$) is the
vanishing point, in pixel coordinates, perpendicular with the roadway
direction of travel, and ($x_{vph}$,$y_{vph}$) is the vanishing point, in
pixel coordinates for vertically aligned lines within the frame, $\alpha_{l}$
is -1 if the rear of the bounding box is closer to ($x_{vpl}$,$y_{vpl}$) than
the front of the bounding box and 1 otherwise, $\alpha_{w}$ is -1 if the left
side of bounding box is closer to ($x_{vpw}$,$y_{vpw}$) than the right side of
the bounding box and 1 otherwise, and $\alpha_{h}$ is -1 if the bottom of the
bounding box is closer to ($x_{vph}$,$y_{vph}$) than the front of the bounding
box and 1 otherwise.

##  6 Experimental Details

###  6.1 Evaluation Protocol

Each object tracker is run using the detection set from Appendix 5. GPS
trajectories and detections from each camera are obtained at slightly
different times. To account for this, tracking evaluation is performed at
fixed 0.1 second intervals, and each GPS trajectory and object tracklet
position is linearly interpolated at each evaluation time. GPS trajectories
extend somewhat outside of the temporal duration for which detections are
available. Evaluation is performed only over the temporal range for which both
GPS trajectories and detections exist. Moreover, each GPS trajectory is
clipped in the X-range $[0,23000]$ such that the trajectory is always visible
within the field of view of the overall camera network. Evaluation is
performed as in [bernardin2008evaluating]. For all metrics other than HOTA
metrics, a lax IOU of 0.1 is required for an object tracklet and GPS
trajectory to be matched, and matching is performed with the Hungarian
Algorithm for bipartite matching [kuhn1955hungarian]. For each object
tracklet, we use the median reported dimension ($l$,$w$, and $h$, in feet)
over all reported tracklet dimension measurements.

###  6.2 Parameter Settings

Table 3 lists relevant parameter settings for each implemented algorithm.
Kalman filter parameters were empirically fine-tuned using the I24-3D dataset
[gloudemans2023dataset].

Algorithm | Parameter | Value | Description
---|---|---|---
ALL | $t_{max}$ | 2 | Maximum time (sec) between detections before track is terminated
| $t_{min}$ | 2 | Minimum track length (sec)
| $f_{eval}$ | 10 | Evaluation time step (Hz)
SORT [bewley2016simple] | $\sigma_{high}$ | 0.5 | Required object confidence to be included in detection set
| $\phi_{nms}$ | 0.1 | Non-maximal suppression IOU threshold
| $d_{max}$ | 10 | Maximum allowable distance for a match (ft)
| $f_{track}$ | 10 | Tracking time step (Hz)
IOUT [bochinski2017high] | $\sigma_{high}$ | 0.5 | Required object confidence to be included in detection set
| $\phi_{nms}$ | 0.1 | Non-maximal suppression IOU threshold
| $\phi_{min}$ | 0.1 | Minimum IOU for a match
| $f_{track}$ | 15 | Tracking time step (Hz)
KIOU [kiout] | $\sigma_{high}$ | 0.5 | Required object confidence to be included in detection set
| $\phi_{nms}$ | 0.1 | Non-maximal suppression IOU threshold
| $\phi_{min}$ | 0.1 | Minimum IOU for a match
| $f_{track}$ | 10 | Tracking time step (Hz)
ByteTrack (L2) [zhang2022bytetrack] | $\sigma_{high}$ | 0.01 | Required object confidence to be included in detection set
| $\phi_{nms}$ | 0.1 | Non-maximal suppression IOU threshold
| $d_{max}$ | 10 | Maximum allowable distance for a match (ft)
| $\tau_{high}$ | 0.4 | Required confidence to be included in first matching step
| $f_{track}$ | 10 | Tracking time step (Hz)
ByteTrack (IOU) [zhang2022bytetrack] | $\sigma_{high}$ | 0.01 | Required object confidence to be included in detection set
| $\phi_{nms}$ | 0.1 | Non-maximal suppression IOU threshold
| $\phi_{min}$ | 0.1 | Minimum IOU for a match
| $\tau_{high}$ | 0.4 | Required confidence to be included in first matching step
| $f_{track}$ | 10 | Tracking time step (Hz)
Oracle | $\phi_{min}$ | 0.1 | Minimum IOU for a match
| $f_{track}$ | 10 | Tracking time step (Hz)
Table 3: Caption

##  7 Treatment of Personally Identifiable Information

As stated in [gloudemans202324], the primary purpose of the I-24 system is not
to produce video data, but rather is to produce vehicle trajectory data. The
data management plan for the system states that in general, video data is not
released to the public or stored for long periods of time, but occasionally
video data may be released to enable work on training, testing, and validation
of trajectory generation algorithms. Thus, the video released in this work
represents an edge case for the system,

Nevertheless, the release of video data (or any data with individuals
represented) carries with it the risk of releasing personally identifiable
information (PII) on the included individuals. In this work we make every
effort to prevent the leakage of PII to the public, and furthermore make it
infeasible to automatically extract PII such that doing so becomes extremely
onerous to potential bad actors. A three-tiered approach is used: 1.) We
release, more or less, a random hour of data, such that the potential for
capturing a discrete event of interest to a third party is near zero. 2.) We
automatically redact all license plate information from all visible vehicles.
3.) We manually redact any regions containing private property or visible
people.

Even in un-redacted video data, license plates are in almost all cases
impossible to read. Figure 15 provides a typical camera field of view, with
the license plate only about 15 pixels wide and subject to significant
blurring from vehicle motion. (The vehicle used in this image is part of the
GPS instrumented vehicle fleet and does not belong to an individual).
Nevertheless, we run an off-the-shelf license plate redaction algorithm
[silva2021flexible] on all frames and cover all detected license plates with a
black rectangle.

Figure 15: Example license plate from this dataset before redaction. License
plate information is unrecoverable; for reference, the orange numbers on the
rear window of the vehicle are about 5 times as large as license plate text
and are barely discernible. (The pictured vehicle is part of the GPS
instrumented vehicle fleet and does not belong to an individual).

We then manually inspect each video and redact the following sets of
information. These areas are blurred in the released videos. Polygons defining
the redacted regions are also released such that data users can replace these
pixel values as desired for computer vision applications. Figure 16 shows an
example of redacted regions.

  * •

All visible people

  * •

All private property

  * •

Any stopped law enforcement vehicle (we leave $\sim 1$ sec of the stopped
vehicle visible to allow for graceful handling of these vehicles by tracking
algorithms.

Figure 16: Example redacted regions (red outline, blurred) for one camera
field of view .

##  8 Known Data Artifacts and Anomalies

Here, we report a list of known anomalies and artifacts in the data as of
submission time:

  * •

Visible in some cameras on poles 3 and 4, there is a region that is blacked
out in all frames and camera fields of view. This region corresponds to a
private residence. Different from other redaction areas in this dataset, a
virtual mask is applied to this region in all cameras in which it is visible,
such that no visual informatiom from this region is ever recorded. Other
regions are redacted after recording.

  * •

All cameras from Pole 25 are missing; this pole was struck by a vehicle a few
days prior to the recording day and could not be restored in time.

  * •

Homographies on the eastbound roadway side are not defined for cameras on
poles 1-7. Construction work required that temporary solid lane markings be
painted, for which i.) corner points were not uniquely distinguishable and
ii.) no matching aerial imagery exists due to the short term nature of the
construction work.

  * •

Reference homographies for cameras on poles 1 and 2 are defined, but
homographies are not re-estimated for these cameras. The lane markings were
altered between the reference homography day and the recording day meaning
that SIFT-FLANN matching and lane marking re-detection fail.

  * •

Camera P22C04 has a black ring visible on the left portion of the frame due to
a mechanical misalignment of the camera lens and body.

  * •

Occasionally, GPS trajectories have missing recordings for time periods on the
order of $\sim 1$ sec. The onboard sensor filtering attempts to compensate for
this missing data, producing “sawtooth” artifacts in the trajectory. Figure 17
shows an example. The majority of these artifacts were removed during data
refinement, but due some artifacts may still remain.

Figure 17: “Sawtooth” artifact in uncorrected GPS trajectories (lines).
Circles depict corresponding manually annotated vehicle positions.


# DFT+U study of UO2: Correct lattice parameter and electronic band-gap

Mahmoud Payami  Email<EMAIL_ADDRESS>School of Physics & Accelerators,
Nuclear Science and Technology Research Institute, AEOI, P. O. Box 14395-836,
Tehran, Iran

###### Abstract

Hubbard-corrected density functional theory, denoted by DFT+U method, was
developed to enable correct prediction of insulating properties for strongly-
correlated electron systems. UO2 is an example having O-$2p$, U-$6d$, and
U-$5f$ incomplete electronic shells. Usually, researchers apply the Hubbard
correction only to the localized incomplete $5f$ electrons of U atoms and
succeed to predict insulating property and good geometric properties by
tweaking the Hubbard-U parameter. However, it turned out that in such a way it
was impossible to obtain reasonable values for both geometry and electronic
band-gap at the same time. In this work, we show that it is possible to
produce good values for those properties just by applying and tuning the
Hubbard corrections to all incomplete shells of O-$2p$, U-$6d$, and U-$5f$.

##  I Introduction

UO2, as a common fuel for nuclear power reactors, has attracted the interests
of researchers for a better theoretical description within DFT+U approachc.[1,
2, 3, 4] Uranium dioxide has a 3D anti-ferromagnetic (AFM) crystal structure
at temperatures less than $30~{}$K,[5, 6] but usually a simpler 1D-AFM model
is used for the description. Recent XRD experiment[7] has shown that UO2
crystallizes with a cubic space group $Pa\bar{3}$ (No. 205). However, if the
structure is modeled by a slightly different but more symmetric cubic space
group $Fm\bar{3}m$ (No. 225) with experimental lattice constant of 5.47Å,
which is shown in Fig. 1(a), then the structure can be represented by a simple
tetragonal unit cell with 6 atoms as shown in Fig. 1(b).

(a)

(b)

Figure 1: (a)- UO2 crystal structure with cubic space group $Fm\bar{3}m$ (No.
225) and lattice constant of 5.47Å; (b)- description by a simple tetragonal
crystal structure with six atoms. Gray and small red balls represent uranium
and oxygen atoms, respectively.

Experiment has shown[8] that UO2 is electrically an insulator with a gap of
2.10 eV. Ordinary approximations in density-functional theory (DFT) such as
local-density approximation (LDA)[9, 10] or semi-local approximations such as
generalized gradient approximation (GGA) [11] for the localized orbitals
usually lead to incorrect metallic behavior. One workaround is to estimate the
interactions of localized orbitals using the Hubbard model and add it to the
DFT energy functional and then subtract the double-counting contributions from
the DFT energy functional:[12, 13, 1, 14]

$E_{DFT+U}=E_{DFT}-E_{dc}+E_{Hub}.$ (1)

The interaction term in Hubbard model, when the Hamiltonian is represented in
the basis of strongly localized Wannier functions, is written as:

$E_{Hub}=U\sum_{i}n_{i\uparrow}n_{i\downarrow}$ (2)

where $U$ is a real number, $n_{i\sigma}$ with $\sigma=\uparrow,\downarrow$
denote the particle number operators, and $i$ specifies the lattice site ${\bf
R}_{i}$. For positive values of $U$, the interaction behaves as on-site
repulsion among the electrons, while on the other hand, negative values of $U$
means that there exist on-site attraction among electrons.

In previous DFT+U calculations for UO2, the on-site Hubbard correction with
positive interaction parameter $U-5f$ was applied to only $5f$ electrons of
uranium atoms which led to gap opening and thus correct insulating behavior.
However, the gap size and geometric properties such as lattice constant both
depend on the interaction parameter. By tuning this on-site parameter, it is
possible to reproduce only one of those properties: gap size or lattice
constant. As is seen from Fig. 2, for the approximation used here, the correct
band gap is reproduced by assuming $U_{\rm{U}-5f}$=3.2 eV while the correct
lattice constant is reproduced by taking $U_{\rm{U}-5f}$=4.8 eV. In this work,
we have extended the Hubbard correction to cover $6d$ of U atoms as well as
$2p$ of O atoms, and determine the relevant interaction parameter values that
reproduce both band-gap and lattice constant of the GS in very good agreement
with experiment.

The organization of this paper is as follows. In Section II, we explain the
computational details; in Section III the calculated results are presented and
discussed; and finally in Section IV we concludes this work.

##  II Computational details

The DFT+U calculations are done by solution of the KS equations using the
Quantum-ESPRESSO code package [15, 16]. Ultra-soft pseudo-potentials (USPP)
are used for U and O atoms that has been generated by the atomic code, using
the generation inputs (with small modifications for more desired results) from
the pslibrary, [17] at https://github.com/dalcorso/pslibrary. The valence
configurations of U($6s^{2},\,6p^{6},\,7s^{2},\,7p^{0},\,6d^{1},\,5f^{3}$) and
O($2s^{2},\,2p^{4}$) were used in the generation. The relativistic effects
were accounted at the level of scalar-relativistic (SR) approximation,[18]
which has been shown to give reasonable GS geometric results[4] for
$U_{\rm{U}-5f}$=4.53 eV when the Perdew-Zunger [19] (PZ) LDA approximation was
used for the exchange-correlation, and the projection on to Hubbard orbitals
were chosen to be atomic ones that were not orthonormalized. The appropriate
kinetic energy cutoffs for the plane-wave expansions were chosen as 90 and 720
Ry for the wavefunctions and densities, respectively. Also, the Methfessel-
Paxton smearing method [20] for the occupations with a width of 0.01 Ry is
used for better convergency process. For the Brillouin-zone integrations in
geometry optimizations, a $8\times 8\times 6$ grid were used. All geometries
were fully optimized for total residual pressures on unit cells to within 0.5
kbar, and residual forces on atoms to within 10-3 mRy/a.u. The occupation
matrix control (OMC)[1] is used to avoid metastable states. The starting
magnetization for oxygen atoms are set to zero values and for U atoms they are
set to $\pm 0.5$ to make anti-ferromagnetic (AFM) configuration along the $z$
direction. Since in the present work we apply Hubbard corrections to $5f$,
$6d$ localized orbitals of U atoms and $2p$ orbitals of O atoms, we have
occupation matrices of dimensions $7\times 7$, $5\times 5$, and $3\times 3$,
respectively. Applying the Hubbard correction for each of the orbitals U-$5f$,
U-$6d$ and O-$2p$ separately one at a time showed that only the Hubbard
corrections to U-$5f$ orbitals lead to metastable states and the other two are
insensitive to initial occupations. Since in the U-atom pseudo-potential the
$5f$ orbital is occupied by 3 electrons, then we have $C_{3}^{7}=35$ different
ways for occupying the diagonal elements of $7\times 7$ matrix by 3 electrons:
[1110000], [1101000], [1100100], $\cdots$, [0001011], [0000111].

##  III Results and Discussions

Examining Hubbard corrections for $5f$, $6d$ of U-atom and $2p$ of O-atom
separately one at a time shows that only correction on $5f$ is able to open an
energy gap and give reasonable lattice constant. The situation is shown in
Figs. 2-4.

(a)

(b)

Figure 2: (a)- Lattice constant $a$ in $\AA$ and band gap in eV as functions
of Hubbard correction strength $U_{\rm{U}-5f}$; (b)- Deviation form cubic
geometry, ($c-a$), in $\AA$ as function of Hubbard correction strength
$U_{\rm{U}-5f}$ of U-atom.

Inspecting Fig. 2(a) it is seen that the value for experimental lattice
constant is reproduced around $U_{\rm{U}-5f}\sim 4.8$ eV while the correct
band gap is reproduced by assuming $U_{\rm{U}-5f}$=3.2 eV. This implies that
with only one correction parameter (i.e., $U_{\rm{U}-5f}$) one fails to
reproduce reasonable values for both lattice constant and the band gap at the
same time.

In Fig. 2(b), the deviation from cubic geometry ($c-a$) is shown to be very
small, $\sim 0.01\AA$, so that modeling the system by 1D-AFM (instead of
3D-AFM) does not cause any significant error in this study. On the other hand,
Figs. 3 and 4 show that the value of lattice constant is relatively
insensitive to the values $U_{\rm{U}-6d}$ and $U_{\rm{O}-2p}$. These results
hint that one should do fine-tuning of $U_{\rm{U}-5f}$ around the value of 4.0
eV. In addition, similar to the result in Fig. 2(b), the deviations from cubic
geometries in Figs. 3-4 are negligible.

Figure 3: Lattice constant $a$ and deviation from cubic geometry ($c-a$) in
$\AA$ as functions of Hubbard correction strength $U_{\rm{U}-6d}$. Figure 4:
The same as in Fig. 3 but for Hubbard correction strength $U_{\rm{O}-2p}$.

In the next step, we apply Hubbard corrections to two orbitals at a time. In
above it was shown that the correct lattice constants were reproduced by
applying the correction to only U-$5f$ with the strength of $\sim 4.00$ eV.
So, we consider the correction to U-$5f$ with strength 4.00 eV as the main
correction and add that for U-$6d$ as a background one with different values.
The result is shown in Fig. 5. As is seen from Fig. 5, adding the background
correction for U-$6d$ almost does not change the lattice constant for
$U_{\rm{U}-6d}<3.0$ eV and so we ignore the background correction for U-$6d$.

Figure 5: The same as in Fig. 3 as functions of Hubbard correction strength
$U_{\rm{U}-6d}$ for fixed value of $U_{\rm{U}-5f}$=4.00 eV.

We now concentrate on adding the background correction of O-$2p$ orbitals. As
is seen from Fig. 6, in contrast to the case of U-$6d$, here the background
correction to O-$2p$ orbitals significantly modifies the results attained by
the correction on U-$5f$. That is, in order to maintain the reasonable value
for the lattice constant, one should use negative values for the Hubbard
correction parameter for O-$2p$ orbitals, meaning that the background
correction is as an on-site attraction. To summarize, Fig. 6 indicates that
the combination of two simultaneous corrections with $U_{\rm{U}-5f}=4.0$ eV
and $U_{\rm{O}-2p}=-3.00$ eV revives the reasonable value for the lattice
constant. But now we expect that the electronic band gap is changed from the
value 2.91 eV, obtained if only U-$5f$ correction was applied.[4] Fig. 6 also
indicates that the deviation from cubic geometry is still acceptable.

Figure 6: The same as in Fig. 5 as functions of Hubbard correction strength
$U_{\rm{O}-2p}$ for fixed value of $U_{\rm{U}-5f}$=4.00 eV.

To have a closer inspection on the effect of negative values for
$U_{\rm{O}-2p}$, we have calculated the GS lattice constants and electronic
band gaps for different values of $U_{\rm{U}-5f}$, keeping $U_{\rm{O}-2p}$
fixed at three values of -3.00, -3.50, and -4.00 eV. The results are presented
in Table 1.

Table 1: Lattice constants in $\AA$ and $E_{gap}$ in eV as functions of $U_{\rm{U}-5f}$ and fixed values of background correction $U_{\rm{O}-2p}$=-3.00, -3.50, and -4.00 eV. $U_{\rm{O}}$ (eV) |  $U_{\rm{U}}$ (eV) |  $a$ ($c$) ($\AA$) |  $E_{gap}$ (eV)
---|---|---|---
-3.00 | 3.00 | 5.4477 (5.4667) | 1.7651
| 3.10 | 5.4498 (5.4684) | 1.8397
| 3.20 | 5.4518 (5.4701) | 1.9096
| 3.30 | 5.45391(5.4718) | 1.9706
| 3.40 | 5.4559 (5.4734) | 2.0308
| 3.50 | 5.4579 (5.4751) | 2.0900
| 3.60 | 5.4599 (5.4767) | 2.1484
| 3.70 | 5.4619 (5.4784) | 2.2059
| 3.80 | 5.4644 (5.4795) | 2.2623
| 3.90 | 5.4666 (5.4817) | 2.3180
| 4.00 | 5.4682 (5.4838) | 2.3731
-3.50 | 3.00 | 5.4499 (5.4693) | 1.7387
| 3.10 | 5.4520 (5.4710) | 1.8125
| 3.20 | 5.4541 (5.4727) | 1.8789
| 3.30 | 5.4561 (5.4744) | 1.9387
| 3.40 | 5.4581 (5.4759) | 1.9977
| 3.50 | 5.4601 (5.4776) | 2.0556
| 3.60 | 5.4622 (5.4794) | 2.1129
| 3.70 | 5.4642 (5.4812) | 2.1691
| 3.80 | 5.4665 (5.4830) | 2.2245
| 3.90 | 5.4685 (5.4847) | 2.2789
| 4.00 | 5.4705 (5.4864) | 2.3320
-4.00 | 3.00 | 5.4522 (5.4720) | 1.7112
| 3.10 | 5.4543 (5.4737) | 1.7842
| 3.20 | 5.4564 (5.4754) | 1.8472
| 3.30 | 5.4584 (5.4769) | 1.9057
| 3.40 | 5.4604 (5.4786) | 1.9636
| 3.50 | 5.4626 (5.4805) | 2.0206
| 3.60 | 5.4648 (5.4823) | 2.0765
| 3.70 | 5.4669 (5.4841) | 2.1315
| 3.80 | 5.4689 (5.4858) | 2.1855
| 3.90 | 5.4708 (5.4874) | 2.2383
| 4.00 | 5.4728 (5.4891) | 2.2902

In order to estimate the proper combinations of Hubbard strengths for
$U_{\rm{U}-5f}$ and $U_{\rm{O}-2p}$ for a desired value of band gap (2.00,
2.10, 2.20 eV), we have plotted the data of Table 1 in Fig. 7.

Figure 7: Variation of band gap as function of $U_{\rm{U}-5f}$ for different
fixed values of $U_{\rm{O}-2p}$. It is seen that, in this example, to obtain
each desired values of $E_{gap}$=2.00, 2.10, and 2.20 eV one has three choices
for the Hubbard strength combinations.

As we see from Fig. 7, since here we have chosen three fixed values for
$U_{\rm{O}-2p}$, there exist three different combinations of the Hubbard
strengths $U_{\rm{U}-5f}$ and $U_{\rm{O}-2p}$ for each desired value of band
gap. To verify the validity of this guess, we have calculated the GS
properties for nine combinations of Hubbard strengths hinted by plots of Fig.
7 and presented the results in Table 2.

Table 2: Hubbard parameters $U_{\rm{U}}$, $U_{\rm{O}}$, in eV, needed to be used to obtain a desired band gap $\tilde{E}_{gap}$ along with the resulted lattice constants $a\;(c)$ and band gap $E_{gap}$.  $\tilde{E}_{gap}$ (eV) |  $U_{\rm{U}}$, $U_{\rm{O}}$ (eV) |  $a$ ($c$) ($\AA$) |  $E_{gap}$ (eV)
---|---|---|---
2.00 | 3.40, -3.00 | 5.4560 (5.4735) | 2.03
| 3.45, -3.50 | 5.4592 (5.4768) | 2.03
| 3.48, -4.00 | 5.4622 (5.4802) | 2.01
2.10 | 3.50, -3.00 | 5.4579 (5.4751) | 2.09
| 3.60, -3.50 | 5.4623 (5.4795) | 2.11
| 3.65, -4.00 | 5.4658 (5.4832) | 2.10
2.20 | 3.70, -3.00 | 5.4619 (5.4784) | 2.21
| 3.78, -3.50 | 5.4661 (5.4827) | 2.21
| 3.84, -4.00 | 5.4697 (5.4865) | 2.21

From the data in Table 2, we see that applying simultaneous Hubbard on-site
corrections on the U-$5f$ and O-$2p$ orbitals it is possible to tune both
lattice constant and band gap to their experimental values.

##  IV Conclusions

In previous theoretical studies of UO2 crystal, in order to predict correct
insulating behavior, researchers used Hubbard corrections for the U-$5f$
localized orbitals in the DFT+U approach. It was already shown that depending
on what XC functional is used and whether the Hubbard orbitals were
orthonormalized or not, for a given Hubbard-$U$ parameter (say 4.0 eV)
different results were obtained for equilibrium lattice constant and the KS
band gap. None of those results were satisfactory in predicting simultaneous
reasonable values for the lattice constant and the size of band gap. In this
work, employing LDA-PZ scheme for the XC energy functional, we have shown that
applying the on-site Hubbard corrections simultaneously to U-$5f$ and O-$2p$
orbitals one can choose certain values to obtain results for both the lattice
constant and energy band gap of the ground state in good agreement with
experiment.

## Acknowledgement

This work is part of research program in School of Physics and Accelerators,
NSTRI, AEOI.

## Data availability

The raw or processed data required to reproduce these results can be shared
with anybody interested upon sending an email to M. Payami.

## References

  * Dorado  _et al._ [2009] B. Dorado, B. Amadon, M. Freyss, and M. Bertolus, DFT+ U calculations of the ground state and metastable states of uranium dioxide, Phys. Rev. B 79, 235125 (2009).
  * Payami [2022] M. Payami, Spin-symmetry broken ground-state of UO2 in DFT+U approach: The SMC method, Iranian J. Phys. Res. 22 (2022).
  * Payami and Sheykhi [2022] M. Payami and S. Sheykhi, Relativistic effects in the study of structure and electronic properties of UO2 within DFT+U method, arXiv preprint arXiv:2209.00724 (2022).
  * Payami [2023] M. Payami, Comparison of SMC and OMC results in determining the ground-state and meta-stable states solutions for UO2 in DFT+U method, arXiv preprint arXiv:2302.04231 (2023).
  * Amoretti  _et al._ [1989] G. Amoretti, A. Blaise, R. Caciuffo, J. M. Fournier, M. T. Hutchings, R. Osborn, and A. D. Taylor, 5f-electron states in uranium dioxide investigated using high-resolution neutron spectroscopy, Phys. Rev. B 40, 1856 (1989).
  * Faber  _et al._ [1975] J. Faber, G. H. Lander, and B. R. Cooper, Neutron-diffraction study of UO2: Observation of an internal distortion, Phys. Rev. Lett. 35, 1770 (1975).
  * Desgranges  _et al._ [2017] L. Desgranges, Y. Ma, P. Garcia, G. Baldinozzi, D. Siméone, and H. Fischer, What is the actual local crystalline structure of uranium dioxide, UO2? a new perspective for the most used nuclear fuel, Inorg. Chem. 56, 321 (2017).
  * Schoenes [1978] J. Schoenes, Optical properties and electronic structure of UO2, J. Appl. Phys. 49, 1463 (1978).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * Cococcioni and De Gironcoli [2005] M. Cococcioni and S. De Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA+ U method, Phys. Rev. B 71, 035105 (2005).
  * Himmetoglu  _et al._ [2014] B. Himmetoglu, A. Floris, S. De Gironcoli, and M. Cococcioni, Hubbard-corrected DFT energy functionals: The LDA+ U description of correlated systems, Int. J. Quantum Chem. 114, 14 (2014).
  * Freyss  _et al._ [2012] M. Freyss, B. Dorado, M. Bertolus, G. Jomard, E. Vathonne, P. Garcia, and B. Amadon,  _First-principles DFT+U study of radiation damage in UO2: $f$ electron correlations and the local energy minima issue_, $\psi_{k}$ Scientific Highlight Of The Month No. 113 (2012).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2020] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum ESPRESSO toward the exascale, J. Chem. Phys. 152, 154105 (2020).
  * Dal Corso [2014] A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Comput. Mater. Sci. 95, 337 (2014).
  * Koelling and Harmon [1977] D. Koelling and B. Harmon, A technique for relativistic spin-polarised calculations, J. Phys. C: Solid State Physics 10, 3107 (1977).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981).
  * Methfessel and Paxton [1989] M. Methfessel and A. Paxton, High-precision sampling for brillouin-zone integration in metals, Phys. Rev. B 40, 3616 (1989).