\section{Introduction} \label{sec:intro}

A large part of the real world almost never moves.
This may be surprising, since moving entities easily attract our attention \cite{franconeri2003moving}, and because we ourselves spend most of our waking hours continuously in motion. Objects like roads and buildings, however, stay put.
Can we leverage this property of the world to learn visual features suitable for interpreting complex scenes with many moving objects?

In this paper, we hypothesize that a correspondence module learned in static scenes will also work well in dynamic scenes. This is motivated by the fact that the content of dynamic scenes is the same as the content of static scenes. We would like our hypothesis to be true, because correspondences are far cheaper to obtain in static scenes than in dynamic ones. In static scenes, one can simply deploy a Simultaneous Localization and Mapping (SLAM) module to obtain a 3D reconstruction of the scene, and then project reconstructed points back into the input imagery to obtain multiview correspondence labels. Obtaining correspondences in dynamic scenes would require taking into account the motion of the objects (i.e., tracking).

We propose to leverage multiview data of static points in arbitrary scenes (static or dynamic),
to learn a neural 3D mapping module which produces features that are correspondable across viewpoints \textit{and} timesteps. The neural 3D mapper consumes RGB-D (color and depth) data as input, and produces a 3D voxel grid of deep features as output. We train the voxel features to be correspondable across viewpoints, using a contrastive loss. At test time, given an RGB-D video with approximate camera poses, and given the 3D box of an object to track, we track the target object by generating a map of each timestep and locating the object's features within each map.

In contrast to models that represent video streams in 2D \cite{DBLP:journals/corr/WangG15a,DBLP:journals/corr/abs-1903-07593,vondrick2018tracking}, our model's 3D scene representation is disentangled from camera motion and projection artifacts.
This provides an inductive bias that scene elements maintain their size and shape across changes in camera viewpoint, and reduces the need for scale invariance in the model's parameters. Additionally, the stability of the 3D map under camera motion allows the model to constrain correspondence searches to the 3D area where the target was last seen, which is a far more reliable cue than 2D pixel coordinates. In contrast to models that use 2.5D representations (e.g., scene flow \cite{Menze2015CVPR}), our neural 3D maps additionally provide features for partially occluded areas of the scene, since the model can infer their features from context. This provides an abundance of additional 3D correspondence candidates at test time, as opposed to being limited to the points observed by a depth sensor.

Our work builds on geometry-aware recurrent neural networks (GRNNs) \cite{commonsense}. GRNNs are modular differentiable neural architectures that take as input RGB-D streams of a static scene under a moving camera and infer 3D feature maps of the scene, estimating and stabilizing against camera motion. The work of Harley et al.~\cite{adam3d} showed that training GRNNs for contrastive view prediction in static scenes helps semi-supervised 3D object detection, as well as moving object discovery. In this paper, we extend those works to learn from dynamic scenes with independently moving objects, and simplify the GRNN model by reducing the number of losses and modules.
Our work also builds on the work of Vondrick et al.~\cite{vondrick2018tracking}, which showed that 2D pixel trackers can emerge without any tracking labels, through self-supervision on a colorization task. In this work, we show that 3D voxel trackers can emerge without tracking labels, through contrastive self-supervision on static data. In the fact that we learn features from correspondences established through triangulation, our work is similar to Dense Object Nets~\cite{florence2018dense}, though that work used object-centric data, and used background subtraction to apply loss on static objects, whereas in our work we do moving-object subtraction to apply loss on \textit{anything} static. We do not assume a priori that we know which objects to focus on.

We test the proposed architectures in simulated and real datasets of urban driving scenes (CARLA~\cite{Dosovitskiy17} and KITTI~\cite{Geiger2013IJRR}). We evaluate the learned 3D visual feature representations on their ability to track objects over time in 3D. We show that the learned visual feature representations can accurately track objects in 3D, simply supervised by observing static data. Furthermore, our method outperforms 2D and 2.5D baselines, demonstrating the utility of learning a 3D representation for this task instead of a 2D one.

The main contribution of this paper is to show that learning feature correspondences from static 3D points causes 3D object tracking to emerge. We additionally introduce a neural 3D mapping module which simplifies prior works on 3D inverse graphics, and learns from a simpler objective than considered in prior works.
Our code and data are publicly available\footnote{\url{https://github.com/aharley/neural_3d_tracking}}.

\section{Related  Work} \label{sec:related}

\subsection{Learning to see by moving}

Both cognitive psychology and computational vision have realised the importance of motion for the development of visual perception \cite{Gibson1979-GIBTEA,wiskott2002slow}.
Predictive coding theories \cite{rao1999predictive,friston2003learning} suggest that the brain
predicts observations at various levels of abstraction; temporal prediction is thus of central interest. These theories currently have extensive empirical support: stimuli are processed more quickly if they are predictable \cite{mcclelland1981interactive,pinto2015expectations}, prediction error is reflected in increased neural activity \cite{rao1999predictive,brodski2015faces}, and disproven expectations lead to learning \cite{schultz1997neural}. Several computational models of frame predictions have been proposed   \cite{rao1999predictive,Eslami1204,DBLP:journals/corr/TatarchenkoDB15,commonsense,oord2018representation}.
Alongside future frame prediction, predicting some form of contextual or missing information has also been explored, such as
predicting frame ordering \cite{lee2017unsupervised}, temporal instance-level associations \cite{DBLP:journals/corr/WangG15a}
color from grayscale \cite{vondrick2018tracking}, egomotion  \cite{DBLP:journals/corr/JayaramanG15,DBLP:journals/corr/AgrawalCM15} and motion trajectory forecasting \cite{DBLP:journals/corr/WalkerDGH16}. Most of these unsupervised methods are evaluated as pre-training mechanisms for object detection or classification \cite{DBLP:journals/corr/MisraZH16,DBLP:journals/corr/WalkerDGH16,DBLP:journals/corr/WangG15a,adam3d}.

Video motion segmentation literature  explores the use of videos in unsupervised moving object discovery \cite{Pont-Tuset_arXiv_2017}.
Most motion segmentation methods operate in 2D image space, and cluster 2D optical flow vectors or 2D flow trajectories  to segment moving objects \cite{OB11,springerlink:10.1007/978-3-642-15555-0_21,Fragkiadaki:topology}, or use low-rank trajectory constraints \cite{DBLP:conf/iccv/CheriyadatR09,Tomasi:1992:SMI:144398.144403,10.1109/ICCV.1995.466815}.
Our work differs in  that we address object detection and segmentation in 3D as opposed to 2D, by estimating 3D motion of the ``imagined" (complete) scene, as opposed to 2D motion of the  pixel observation stream.

\subsection{Vision as inverse graphics}

Earlier works in Computer Vision proposed  casting visual recognition as  inverse rendering \cite{bruno,yuillek06}, as opposed to feedforward labelling. The ``blocks world" of Roberts \cite{blocksworld} had the goal of reconstructing the 3D scene depicted in the image in terms of 3D solids found in a database.
A key question to be addressed is: what representations should we use for the  intermediate latent 3D structures? Most works seek to map images to explicit 3D representations, such as 3D pointclouds \cite{Wu2016,DBLP:journals/corr/TungHSF17,tinghuisfm,sfmnet}, 3D meshes \cite{Loper:2015:SSM:2816795.2818013,DBLP:journals/corr/abs-1711-07566}, or binary 3D voxel occupancies \cite{DBLP:journals/corr/TulsianiZEM17,LSM,3Dshapenets}.
The aforementioned manually designed 3D representations, e.g., 3D meshes, 3D keypoints, 3D pointclouds, 3D voxel grids, may not be general enough to express the rich 3D world, which contains liquids, deformable objects, clutter, dirt, wind, etc., and at the same time may be over descriptive when detail is unnecessary. In this work, we opt for \textit{learning-based 3D feature representations} extracted end-to-end from the RGB-D input stream as proposed by Tung et al.~\cite{commonsense} and Harley et al.~\cite{adam3d}. We extend the architectures of those works to handle and learn from videos of dynamic scenes, as opposed to only videos of static scenes.

\section{Learning to Track with Neural 3D Mapping} \label{sec:model}

We consider a mobile agent that can move around in a 3D scene, and observe it from multiple viewpoints. The scene can optionally contain dynamic (moving and/or deforming) objects. The agent has an RGB camera with known intrinsics, and a depth sensor registered to the camera's coordinate frame.
It is reasonable to assume that a mobile agent who moves at will has access to its approximate egomotion, since it chooses where to move and what to look at \cite{patla1991visual}. In simulation, we use ground truth camera poses; in real data, we use approximate camera poses provided by an inertial navigation system (GPS and IMU). In simulation, we use random viewpoints; in real data, we use just one forward-facing camera (which is all that is available). Note that a more sophisticated mobile agent might attempt to select viewpoints intelligently at training time.

Given the RGB-D and egomotion data, our goal is to learn 3D feature representations that can correspond entities across time, despite variations in pose and appearance. We achieve this by training inverse graphics neural architectures that consume RGB-D videos and infer 3D feature maps of the full scenes, as we describe in Section~\ref{sec:neural_mapping}. To make use of data where some parts are static and other parts are moving, we learn to identify static 3D points by estimating a reliability mask over the 3D scene, as we describe in Section~\ref{sec:reliability}. Finally, we track in 3D, by re-locating the object within each timestep's 3D map, as described in Section~\ref{sec:tracking}. Figure~\ref{fig:overview} shows an overview of the training and testing setups.

\subsection{Neural 3D Mapping} \label{sec:neural_mapping}

Our model learns to map an RGB-D (RGB and depth) image to a 3D feature map of the scene in an end-to-end differentiable manner. The basic architecture is based on prior works \cite{commonsense,adam3d}, which proposed view prediction architectures with a 3D bottleneck. In our case, the 3D feature map is the \textit{output} of the model, rather than an intermediate representation.

Let $\M \in \mathbb{R}^{W \times H \times D \times C}$ denote the 3D feature map representation, where $W,H,D,C$ denote the width, height, depth and number of feature channels, respectively. The map corresponds to a large cuboid of world space, placed at some pose of interest (e.g., surrounding a target object). Every $(x,y,z)$ location in the 3D feature map $\M$ holds a $C$-length feature vector that describes the semantic and geometric content of the corresponding location of the world scene. To denote the feature map of timestep $i$, we write $\Mi$. We denote the function that maps RGB-D inputs to 3D feature maps as $f: (\I, \D) \mapsto \M$. To implement this function, we voxelize the inputs into a 3D grid, then pass this grid through a 3D convolutional network, and $L_2$-normalize the outputs.

Tung et al.~\cite{commonsense} learned the parameters of $f$ by predicting RGB images of unseen viewpoints, and applying a regression loss; Harley et al.~\cite{adam3d} demonstrated that this can be outperformed by contrastive prediction objectives, in 2D and 3D. Here, we drop the view prediction task altogether, and focus entirely on a 3D correspondence objective: if a static point $(x,y,z)$ is observed in two views $i$ and $j$, the corresponding features $m_i = \Mi_{x,y,z}, m_j = \Mj_{x,y,z}$ should be similar to each other, and distinct from other features. We achieve this with a cross entropy loss \cite{sohn2016improved,oord2018representation,he2019momentum,chen2020simple}:
\begin{equation}

    \mathcal{L}_{i,j} = -\log \frac{ \exp(  m_i^\top m_j / \tau)}{\sum_{k \neq i} \exp( m_i^\top m_k / \tau) },
\end{equation}
where $\tau$ is a temperature parameter, which we set to $0.07$, and the sum over $k \neq i$ iterates over non-corresponding features. Note that indexing correctly into the 3D scene maps to obtain the correspondence pair $m_i, m_j$ requires knowledge of the relative camera transformation across the input viewpoints; we encapsulate this registration and indexing in the notation $\M_{x,y,z}$. Following He et al.~\cite{he2019momentum}, we obtain a large pool of negative correspondences through the use of an offline dictionary, and stabilize training with a ``slow'' copy of the encoder, $f_\textrm{slow}$, which is learned via high-momentum updates from the main encoder parameters.

Since the neural mapper (a) does not know a priori which voxels will be indexed for a loss, and (b) is fully convolutional, it learns to generate view-invariant features densely in its output space, even though the supervision is sparse. Furthermore, since (a) the model is encouraged to generate correspondable features invariant to viewpoint, and (b) varying viewpoints provide varying contextual support for 3D locations, the model learns to \textit{infer} corrrespondable features from limited context, which gives it robustness to partial occlusions.

\subsection{Inferring static points for self-supervision in dynamic scenes}
\label{sec:reliability}

The training objective in the previous subsection requires the location of a static point observed in two or more views. In a scenario where the data is made up entirely of static scenes, as can be achieved in simulation or in controlled environments, obtaining these static points is straightforward: any point on the surface of the scene will suffice, provided that it projects into at least two camera views.

To make use of data where some parts are static and other parts are moving, we propose to simply discard data that appears to be moving. We achieve this by training a neural module to take the difference of two scene features $\Mi, \Mj$ as input, and output a ``reliability'' mask indicating a per-voxel confidence that the scene cube within the voxel is static: $ g : (\textbf{sg}(\Mi-\Mj)) \mapsto \mathbb{R}^{W,H,D}$, where $\textbf{sg}$ stops gradients from flowing from $g$ into the function $f$ which produces $\Mi, \Mj$. We implement $g$ as a per-voxel classifier, with 2-layer fully-connected network applied fully convolutionally. We do not assume to have true labels of moving/static voxels, so we generate synthetic labels using static data: given two maps of the \textit{same} scene $\Mi, \Mj$, we generate positive-label inputs with $\Mi - \Mj$ (as normal), and generate negative-label inputs with $\Mi - \textrm{shuffle}(\Mj)$, where the shuffle operation ruins the correspondence between the two tensors. After this training, we deploy this network on pairs of frames from dynamic scenes, and use it to select high-confidence static data to further train the encoder $f$.

Our training procedure is then: (1) learn the encoder $f$ on static data; (2) learn the reliability function $g$; (3) in dynamic data, finetune $f$ on data selected by $g$. Steps 2-3 can be repeated a number of times. In practice we find that results do not change substantially after the first pass.

\subsection{Tracking via point correspondences} \label{sec:tracking}

Using the learned neural 3D mapper, we track an object of interest over long temporal horizons by re-locating it in a map produced at each time step. Specifically, we re-locate each voxel of the object, and use these new locations to form a motion field. We convert these voxel motions into an estimate for the entire object by fitting a rigid transformation to the correspondences, via RANSAC.

We assume we are given the target object's 3D box on the zeroth frame of a test video. Using the zeroth RGB-D input, we generate a 3D scene map centered on the object. We then convert the 3D box into a set of coordinates $X_0, Y_0, Z_0$ which index into the map. Let $m_i = \Mz_{x_i \in X_0, y_i \in Y_0, z_i \in Z_0}$ denote a voxel feature that belongs to the object. On any subsequent frame $t$, our goal is to locate the new coordinate of this feature, denoted $x_i^t, y_i^t, z_i^t$. We do so via a soft spatial argmax, using the learned feature space to provide correspondence confidences:
\begin{equation}
    (x_i^t, y_i^t, z_i^t) = \sum_{x \in X, y \in Y, z \in Z}  \left( \frac{ \exp(m_i^\top \Mt_{x,y,z})}{\sum_{x_2 \in X, y_2 \in Y, z_2 \in Z} \exp( m_i^\top \Mt_{x_2,y_2,z_2} )}(x,y,z) \right),
    \label{eq:reloc}
\end{equation}
where $X, Y, Z$ denote the set of coordinates in the search region. We then compute the motion of the voxel as $(\delta_{x_i}^t, \delta_{y_i}^t, \delta_{z_i}^t) = (x_i^t, y_i^t, z_i^t) - (x_i, y_i, z_i)$. After computing the motion of every object voxel in this way, we use RANSAC to find a rigid transformation that explains the majority of the correspondences.
We apply this rigid transform to the input box, yielding the object's location in the next frame.

Vondrick et al.~\cite{vondrick2018tracking} computed a similar attention map during tracking (in 2D), but did not compute its soft argmax, nor explain the object's motion with a single transformation, but rather propagated a ``soft mask'' to the target frame, which is liable to grow or shrink. Our method takes advantage of the fact that all coordinates are 3D, and makes the assumption that the objects are rigid, and propagates a fixed-size box from frame to frame.

We empirically found that it is critical to constrain the search region of the tracker in 3D. In particular, on each time step we create a search region centered on the object's last known position. The search region is $16m \times 2m \times 16m$ large, which is half of the typical full-scene resolution. This serves three purposes. The first is: it limits the number of spurious correspondences that the model can make, since it puts a cap on the metric range of the correspondence field, and thereby reduces errors. Second: it ``re-centers'' the model's field of view onto the target, which allows the model to incorporate maximal contextual information surrounding the target. Even if the bounds were sufficiently narrow to reduce spurious correspondences, an object at the \textit{edge} of the field of view will have less-informed features than an object at the middle, due to the model's convolutional architecture. The third reason is computational: even 2D works \cite{vondrick2018tracking} struggle with the computational expense of the large matrix multiplications involved in this type of soft attention, and in 3D the expense is higher. Searching locally instead of globally makes Eq.~\ref{eq:reloc} tractable.

\section{Experiments}\label{sec:exp}

\begin{figure}[t]
 \centering
 \includegraphics[width=1.0\linewidth]{figures/eccv_fig_2d_tracking_v6.pdf}
 \caption{\textbf{Visualization of tracking inputs, inferred 2D scene features,
 and tracking outputs, for our adaptation of Dense Object Nets \cite{florence2018dense}.} This model operates in the same way as our full 3D one, but creates features in 2D, and unprojects them to a sparse 3D grid for tracking.
 }
 \label{fig:tracking2d_qual}
\end{figure}

We test our model in the following two datasets:
\begin{enumerate}
\item \textbf{Synthetic RGB-D videos of urban scenes rendered in the CARLA simulator \cite{Dosovitskiy17}}. CARLA is an open-source photorealistic simulator of urban driving scenes. It permits moving the camera to any desired viewpoint in the scene.
We obtain data from the simulator as follows. We begin by generating 10000 autopilot episodes of 16 frames each, at 10 FPS.
We define 18 viewpoints along a 40m-radius hemisphere anchored to the ego-car (i.e., it moves with the car).
In each episode, we sample 6 random viewpoints from the 18 and randomly perturb their pose, and then capture each timestep of the episode from these 6 viewpoints. We discard episodes that do not have an object in-bounds for the duration of the episode.

We treat the Town1 data as the ``training'' set, and the Town2 data as the ``test'' set, so there is no overlap between the train and test sets. This yields 4313 training videos, and 2124 test videos.

\item \textbf{Real RGB-D videos of urban scenes, from the KITTI dataset \cite{Geiger2013IJRR}.} This data was collected with a sensor platform mounted on a moving vehicle, with a human driver navigating through a variety of areas in Germany. We use the ``left'' color camera, and LiDAR sweeps synced to the RGB images.

For training, we use the ``odometry'' subset of KITTI; it includes egomotion information accurate to within 10cm.
The odometry data includes ten sequences, totalling 23201 frames.

We test our model in the validation set of the ``tracking'' subset of KITTI, which has twenty labelled sequences, totalling 8008 frames. For supervised baselines, we split this data into 12 training sequences and 8 test sequences. For evaluation, we create 8-frame subsequences of this data, in which a target object has a valid label for all eight frames. This subsequencing is necessary since objects are only labelled when they are within image bounds. The egomotion information in the ``tracking'' data is only approximate.
\end{enumerate}
We evaluate our model on its ability to track objects in 3D. On the zeroth frame, we receive the 3D box of an object to track. On each subsequent frame, we estimate the object's new 3D box, and measure the intersection over union (IOU) of the estimated box with the ground truth box. We report IOUs as a function of timesteps.

\subsection{Baselines}

We evaluate the following baselines. We provide additional implementation details for each baseline (and for our own model) in the supplementary file.

\begin{itemize}
\item \textbf{Unsupervised 3D flow \cite{adam3d}.} This model uses an unsupervised architecture similar to ours, but with a 2-stage training procedure, in which features are learned first from static scenes and frozen, then a 3D flow module is learned over these features in dynamic scenes. We extend this into a 3D tracker by ``chaining'' the flows across time, and by converting the trajectories into rigid motions via RANSAC.

\item \textbf{2.5D dense object nets \cite{florence2018dense}.} This model learns to map input images into dense 2D feature maps, and uses a contrastive objective at known correspondences across views. We train this model using static points for correspondence labels (like our own model).
We extend this model into a 3D tracker by ``unprojecting'' the learned embeddings into sparse 3D scene maps, then applying the \textit{same tracking pipeline} as our own model.

\item \textbf{2.5D tracking by colorization \cite{vondrick2018tracking,mast}.} This model learns to map input images into dense 2D feature maps, using an an RGB reconstruction objective. The model trains as follows: given two RGB frames, the model computes a feature map for each frame; for each pixel of the first frame's feature map, we compute that feature's similarity with all features of the second frame, and then use that similarity matrix to take a weighted combination of the second frame's colors; this color combination at every pixel is used as the reconstruction of the first frame, which yields an error signal for learning.
We extend this model into a 3D tracker in the same way that we extended the ``dense object nets'' baseline.

\item \textbf{3D neural mapping with random features.} This model is equivalent to our proposed model but with randomly-initialized network parameters.
This model may be expected to perform at better-than-chance levels due to the power of random features \cite{rahimi2008random} and due to the domain knowledge encoded in the architecture design.

\item \textbf{3D fully convolutional siamese tracker (supervised) \cite{bertinetto2016fully}.} This is a straightforward 3D upgrade of a fully convolutional 2D siamese tracker, which uses the object's feature map as a cross correlation template, and tracks the object by taking an argmax of the correlation heatmap at each step. It is necessary to supervise this model with ground-truth box trajectories.
We also evaluate a ``cosine windowing'' variant of this model, which suppresses correlations far from the search region's centroid \cite{bertinetto2016fully}.

\item \textbf{3D siamese tracker with random features.} This model is equivalent to the 3D siamese supervised tracker, but with randomly-initialized network parameters. Similar to the random version of 3D neural mapping, this model measures whether random features and the implicit biases are sufficient to track in this domain.

\item \textbf{Zero motion.} This baseline simply uses the input box as its final estimate for every time step. This baseline provides a measure for how quickly objects tend to move in the data, and serves as a lower bound for performance.
\end{itemize}
All of these are fairly simple trackers. A more sophisticated approach might incrementally update the object template \cite{matthews2004template}, but we leave that for future work.

\paragraph{Ablations.} We compare our own model against the following ablated versions. First, we consider a model without search regions, which attempts to find correspondences for each object point in the \textit{entire} 3D map at test time. This model is far more computationally expensive, since it requires taking the dot product of each object feature with the entire scene. Second, we consider a similar ``no search region'' model but at half resolution, which brings the computation into the range of the proposed model. This ablation is intended to reveal the effect of resolution on accuracy. Third, we omit the ``static point selection'' (via the function $g$). This is intended to evaluate how correspondence errors caused by moving objects (violating the static scene assumption) can weaken the model.

\subsection{Quantitative results}

We evaluate the mean 3D IOU of our trackers over time. Figure~\ref{fig:all_ious}-left
shows the results of this evaluation in CARLA.
As might be expected, the supervised 3D trackers perform best, and cosine windowing improves results.

Our model outperforms all other unsupervised models, and nearly matches the supervised performance.
The 2.5D dense object net performs well also,
but its accuracy is likely hurt by the fact that it is limited exclusively to points observed in the depth map. Our own model, in contrast, can match against both observed and unobserved (i.e., hallucinated or inpainted) 3D scene features. The colorization model performs under the 2.5D dense object net approach, likely because this model only indirectly encourages correspondence via the colorization task, and therefore is a weaker supervision than the multi-view correspondence objectives used in the other methods.

Random features perform worse than the zero-motion baseline, both with a neural mapping architecture and a siamese tracking architecture. Inspecting the results qualitatively, it appears that these models quickly propagate the 3D box off of the object and onto other scene elements. This suggests that random features and the domain knowledge encoded in these architectures are not enough to yield 3D trackers in this data.

\begin{figure}[t!]
 \centering
 \includegraphics[width=1.0\linewidth]{figures/tracking_iou_both.pdf}
 \caption{
 \textbf{Single-object tracking accuracy, in mean IOU across time, in CARLA (left) and KITTI (right).} Methods are sorted in the legend in decreasing order of mean IOU. }
 \label{fig:all_ious}
\end{figure}

\begin{comment}
\setlength{\tabcolsep}{4pt}
\begin{table}
\begin{center}
\caption{Single-object tracking accuracy, in mean IOU across time, in CARLA.}
\label{table:tracking_carla}
\begin{tabular}{lllll}
\hline\noalign{\smallskip}
Method & IOU@2 & IOU@4 & IOU@6 & IOU@8 \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Zero motion  & 0.63 & 0.33 & 0.21 & 0.17\\
2.5D dense object nets  & 0.66 & 0.47 & 0.39 & 0.33\\
Harley et al.~\cite{adam3d}  & 0.58 & 0.47 & 0.38 & 0.29\\

Ours  & 0.80 & 0.69 & 0.66 & 0.61\\
\hline
3D siamese (supervised)  & 0.69 & 0.65 & 0.63 & 0.61\\
\hline
\end{tabular}
\end{center}
\end{table}
\setlength{\tabcolsep}{1.4pt}
\end{comment}

We perform the same evaluation in KITTI, and show the results in Figure~\ref{fig:all_ious}-right.
On this benchmark, accuracies are lower for all models, indicating that the task here is more challenging. This is likely related to the fact that (1) the egomotion is imperfect in this data, and (2) the tracking targets are frequently farther away than they are in CARLA. Nonetheless, the ranking of methods is the same, with our 3D neural mapping model performing best. One difference is that cosine windowing actually worsens siamese tracking results in KITTI. This is likely related to the fact that egomotion stabilization is imperfect in KITTI: a zero-motion prior is only helpful in frames where the target is stationary \textit{and} camera motion is perfectly accounted for; otherwise it is detrimental.

We additionally split the evaluation on stationary vs. moving objects in CARLA. We use a threshold of $1m$ total distance (in world coordinates) across the 8-frame trajectories to split these categories. At the last timestep, the mean 3D IOU for all objects together it is $0.61$ (as shown in Figure~\ref{fig:all_ious}-left); for static objects only, the value is $0.64$; for moving objects only, it is $0.56$. This suggests that the model tracks stationary objects more accurately than moving objects, likely because their appearance changes less with respect to the camera and the background.

Finally, we evaluate the top two models in CARLA using the standard 2D tracking metrics, multi-object tracking accuracy (MOTA) and multi-object tracking precision (MOTP) \cite{bernardin2006multiple}, though we note that our task only has one target object per video. We find that the 3D siamese + cosine model (supervised) achieves a MOTA of $0.9578$ and MOTP of $0.7407$, while our model achieves MOTA $0.9125$ and MOTP $0.8384$. This suggests that the supervised tracker makes fewer misses, but our method delivers slightly better precision.

\begin{comment}
\setlength{\tabcolsep}{4pt}
\begin{table}
\begin{center}
\caption{Single-object tracking accuracy, in mean IOU across time, in KITTI.}
\label{table:tracking_kitti}
\begin{tabular}{lllll}
\hline\noalign{\smallskip}
Method & IOU@2 & IOU@4 & IOU@6 & IOU@8 \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Zero motion  & 0.59 & 0.28 & 0.18 & 0.13\\
2.5D dense object nets  & 0.68 & 0.55 & 0.4 & 0.31\\
Harley et al.~\cite{adam3d}  & 0.55 & 0.42 & 0.39 & 0.22\\

Ours  & 0.66 & 0.58 & 0.52 & 0.46\\
\hline
3D siamese (supervised)  & 0.61 & 0.60 & 0.58 & 0.58\\
\hline
\end{tabular}
\end{center}
\end{table}
\setlength{\tabcolsep}{1.4pt}
\end{comment}

\subsection{Qualitative results}

We visualize our tracking results, along with inputs and colorized visualizations of our neural 3D scene maps, in Figure~\ref{fig:tracking3d_qual}. To visualize the neural 3D maps, we take a mean along the vertical axis of the grid (yielding a ``bird's eye view'' 2D grid), compress the deep features in each grid cell to 3 channels via principal component analysis, normalize, and treat these values as RGB intensities.
Comparing the 3D features learned in CARLA vs those learned in KITTI reveals a very obvious difference: the KITTI features appear blurred and imprecise in comparison with the CARLA features. This is likely due to the imperfect egomotion information, which leads to slightly inaccurate correspondence data at training time (i.e., occasional failures by the static point selector $g$).

In Figure~\ref{fig:tracking2d_qual}, we visualize the features learned by Dense Object Nets in this data. From the PCA colorization it appears that objects are typically colored differently from their surroundings, which is encouraging, but the embeddings are not as clear as those in the original work~\cite{florence2018dense}, likely because this domain does not have the benefit of background image subtraction and object-centric losses.

In the supplementary file, we include video visualizations of the learned features and 3D tracking results.

\subsection{Ablations}

We evaluate ablated versions of our model, to reveal the effect of (1) search regions, (2) resolution, and (3) static point selection in dynamic scenes. Results are summarized in  Table~\ref{tab:ablations}.

Without search regions, the accuracy of our model drops by 20 IOU points, which is a strong impact. We believe this drop in performance comes from the fact that search regions take advantage of 3D scene constancy, by reducing spurious correspondences in far-away regions of the scene.

Resolution seems to have a strong effect as well: halving the resolution of the wide-search model  reduces its performance by 15 IOU points. This may be related to the fact that fewer points are then available for RANSAC to find a robust estimate of the object's rigid motion.

Since static point selection is only relevant in data with moving objects, we perform this experiment in KITTI (as opposed to CARLA, where the training domain is all static by design).  The results show that performance degrades substantially without the static point selection. This result is to be expected, since this ablation causes erroneous correspondences to enter the training objective, and thereby weakens the utility of the self-supervision.

\setlength{\tabcolsep}{4pt}
\begin{table}
\begin{center}
\caption{Ablations of our model.}
\label{tab:ablations}
\begin{tabular}{lc}
\hline\noalign{\smallskip}
Method & Mean IOU \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Ours in CARLA  & 0.61\\
\ldots without search regions & 0.40\\
\ldots without search regions, at half resolution & 0.25\\
Ours in KITTI  & 0.46\\
\ldots without static point selection & 0.39\\
\hline
\end{tabular}
\end{center}
\end{table}
\setlength{\tabcolsep}{1.4pt}

\subsection{Limitations}\label{sec:limitations}

The proposed model has three important limitations.
First,
our work assumes access to RGB-D data with accurate egomotion data at training time, with a wide variety of viewpoints. This is easy to obtain in simulators, but real-world data of this sort typically lies along straight trajectories (as it does in KITTI), which limits the richness of the data.
Second, our model architecture requires a lot of GPU memory, due to its third spatial dimension. This severely limits either the resolution or the metric span of the latent map $\M$.
On 12G Titan X GPUs we encode a space sized $32m \times 4m \times 32m$ at a resolution of $128 \times 32 \times 128$, with a batch size of 4; iteration time is ~0.2s/iter.
Sparsifying our feature grid, or using
points instead of voxels,
are clear areas for future work.
Third, our test-time tracking algorithm makes two strong assumptions: (1) a tight box is provided in the zeroth frame, and (2) the object is rigid. For non-rigid objects, merely propagating the box with the RANSAC solution would be insufficient, but the voxel-based correspondences might still be helpful.

\section{Conclusion} \label{sec:conclusion}
We propose a model which learns to track objects in dynamic scenes just from observing static scenes.
We show that a multi-view contrastive loss allows us to learn rich visual representations that are correspondable not only across views, but across time.
We demonstrate the robustness of the learned representation by benchmarking the learned features on a tracking task in real and simulated data. Our approach outperforms prior unsupervised 2D and 2.5D trackers, and approaches the accuracy of supervised trackers.
Our 3D representation benefits from denser correspondence fields than 2.5D methods, and is invariant to the artifacts of camera projection, such as apparent scale changes of objects.
Our approach opens new avenues for learning trackers in arbitrary environments, without requiring explicit tracking supervision: if we can obtain an accurate pointcloud reconstruction of an environment, then we can learn a tracker for that environment too.

\section*{Acknowledgements}
This material is based upon work funded and supported by the Department of Defense under Contract No. FA8702-15-D-0002 with Carnegie Mellon University for the operation of the Software Engineering Institute, a federally funded research and development center. We also acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), AiDTR, the DARPA Machine Common Sense program, and the AWS Cloud Credits for Research program.

\section{Implementation details}

\subsection{Neural 3D mapping}

Our 3D convolutional network $f$ has the following architecture. The input is a 4-channel 3D grid shaped $128 \times 32 \times 128$. The four channels are RGB and occupancy. ``Unoccupied'' voxels have zeros in all channels. This input is passed through a 3D convolutional network, which is an encoder-decoder with skip connections. There are three encoder layers: each have stride-2 3D convolutions, with kernel size $4 \times 4 \times 4$, and the output channel dimensions are $64, 128, 192$. There are three decoder layers: the first two have stride-2 transposed convolutions with stride 2 and $4 \times 4 \times 4$ kernels with $256$ output channels; after each of these we concatenate the same-resolution layer from the encoder half; the last layer has a $1 \times 1 \times 1$ kernel and $64$ channels.

At training time $f$ encodes a metric space sized $32 \times 4 \times 32$ meters. At test time we ``zoom in'' on the area of interest (i.e., the search region), and encode a space sized $16 \times 4 \times 16$ meters.

We train with the Adam optimizer for $200,000$ iterations at a learning rate of $1e-4$.

\subsection{Baselines}

\begin{itemize}
\item \textbf{Unsupervised 3D flow \cite{adam3d}.}
The original work by Harley et al.~\cite{adam3d} only estimated flow for pairs of frames. We extend this into a 3D tracker by deploying the flow module in an object-centric manner, and by ``chaining'' the flows together across time. We effectively compute the flow from the first frame \textit{directly} to every other frame in the sequence, but with an alignment step designed to keep the object within the field of view of the flow module. Specifically, given the object location in the first frame, we extract a volume of flow vectors in the object region, and run RANSAC to find a rigid transformation that explains the majority of the flow field. We apply this rigid transform to the box, yielding the object's location in the next frame. Then, we back-warp the next frame according to the box transformation, to re-center the voxel grid onto the object. Estimating the flow between the original frame and the newly backwarped frame yields the \textit{residual} flow field; adding this residual to the original flow provides the new cumulative motion of the object. We repeat these steps across the length of the input video, back-warping according to the cumulative flow and estimating the residual.

\item \textbf{2.5D dense object nets \cite{florence2018dense}.} While the original work used a margin loss for contrastive learning, we update this to a cross entropy loss with a large offline dictionary and a coupled ``slow'' encoder, consistent with state-of-the-art contrastive learning \cite{he2019momentum}. We extend this model into a 3D tracker as follows. First, we use the available depth data to ``unproject'' the learned embeddings into sparse 3D scene maps. Then, we use the same tracking pipeline as our own model (with search regions and soft argmaxes and RANSAC). The differences between this model and our own are (1) it learns a 2D CNN instead of a 3D one, and (2) its 3D output is constrained to the locations observed in the depth map, instead of producing features densely across the 3D grid.

\item \textbf{2.5D tracking by colorization \cite{vondrick2018tracking,mast}.} We use the latest iteration of this method \cite{mast}, which uses the LAB colorspace at input and output, and color dropout at training time; this outperforms the quantized-color cross entropy loss of the original work.
We extend this model into a 3D tracker in the same way that we extended the ``dense object nets'' baseline. Therefore, the only difference between the two 2.5D methods is in the training: this method can train on arbitrary data but has an objective that only indirectly encourages correspondence (through an RGB reconstruction loss), while ``2.5D dense object nets'' requires static scenes and directly encourages correspondences (through a contrastive loss).

\item \textbf{3D neural mapping with random features.} We use the same inputs, architecture, outputs, and the same hyperparameters for resolution and search regions, but do not train the parameters.

\item \textbf{3D fully convolutional siamese tracker (supervised) \cite{bertinetto2016fully}.} We use the same resolutions and feature encoders for this model as we do for our main model. The cosine window adds a zero-velocity bias into the model, and makes it less likely for the model's outputs to make large jumps \cite{bertinetto2016fully}.
\end{itemize}

\section{Detailed results}

In Tables~\ref{tab:tracking_carla} and~\ref{tab:tracking_kitti}, we provide the numerical values of the data plotted in the main paper. The notation IOU@N denotes 3D intersection over union at the Nth frame of the video. Note that IOU@0 $= 1.0$ in this task, since tracking is initialized with a ground-truth box in the zeroth frame.

\setlength{\tabcolsep}{4pt}
\begin{table}
\begin{center}
\caption{Single-object tracking accuracy, in mean IOU across time, in CARLA.}
\label{tab:tracking_carla}
\begin{tabular}{lcccc}
\hline\noalign{\smallskip}
Method & IOU@2 & IOU@4 & IOU@6 & IOU@8 \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Zero motion  & 0.63 & 0.33 & 0.21 & 0.17\\
Random 3D neural mapping  & 0.13 & 0.07 & 0.05 & 0.03\\
Random 3D siamese  & 0.05 & 0.02 & 0.01 & 0.00\\
Random 3D siamese + cosine  & 0.48 & 0.22 & 0.10 & 0.05\\
2.5D colorization~\cite{vondrick2018tracking,mast}  & 0.41 & 0.29 & 0.25 & 0.19\\
2.5D dense object nets~\cite{florence2018dense}  & 0.66 & 0.47 & 0.39 & 0.33\\
3D flow~\cite{adam3d}  & 0.58 & 0.47 & 0.38 & 0.29\\
Ours  & 0.80 & 0.69 & 0.66 & 0.61\\
\hline
3D siamese~\cite{bertinetto2016fully} (supervised)  & 0.69 & 0.65 & 0.63 & 0.61\\
3D siamese + cosine~\cite{bertinetto2016fully} (supervised)  & 0.75 & 0.71 & 0.69 & 0.65\\
\hline
\end{tabular}
\end{center}
\end{table}
\setlength{\tabcolsep}{1.4pt}

\setlength{\tabcolsep}{4pt}
\begin{table}
\begin{center}
\caption{Single-object tracking accuracy, in mean IOU across time, in KITTI.}
\label{tab:tracking_kitti}
\begin{tabular}{lcccc}
\hline\noalign{\smallskip}
Method & IOU@2 & IOU@4 & IOU@6 & IOU@8 \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Zero motion  & 0.59 & 0.28 & 0.18 & 0.13\\
Random 3D neural mapping  & 0.10 & 0.04 & 0.03 & 0.02\\
Random 3D siamese & 0.00 & 0.01 & 0.00 & 0.00\\
Random 3D siamese + cosine & 0.40 & 0.11 & 0.03 & 0.01\\
2.5D colorization~\cite{vondrick2018tracking,mast}  & 0.30 & 0.20 & 0.11 & 0.10\\
2.5D dense object nets~\cite{florence2018dense}  & 0.68 & 0.55 & 0.40 & 0.31\\
3D flow~\cite{adam3d}  & 0.55 & 0.42 & 0.39 & 0.22\\

Ours  & 0.66 & 0.58 & 0.52 & 0.46\\
\hline
3D siamese~\cite{bertinetto2016fully} (supervised)  & 0.61 & 0.60 & 0.58 & 0.58\\
3D siamese + cosine~\cite{bertinetto2016fully} (supervised) & 0.60 & 0.56 & 0.55 & 0.55\\
\hline
\end{tabular}
\end{center}
\end{table}
\setlength{\tabcolsep}{1.4pt}


\section{Introduction} \label{sec:intro}

The third generation partnership project (3GPP) has extended the support for vehicular communications as a part of the fifth generation new radio (5G-NR) standard \cite{3gpp_rel16}. The wide range of use cases supported by vehicular communications includes the exchange of short safety messages with high reliability and also video streaming services with low latency and high data rates. In order to design efficient networks that meet these stringent requirements, it is imperative to understand the additional burden on the network infrastructure due to vehicular users. The distribution of the number of users served by a cellular MBS, which is referred to as the load on the MBS, is a key metric that represents the demand of network resources.
Inspired by the randomness in the spatial distribution of vehicular users, we employ tools from stochastic geometry to characterize the load on the MBSs due to vehicular users. The prior art in this area is discussed next.

{\em Prior art.} While the distribution of load on MBSs due to users modeled by PPPs has been investigated in the literature \cite{Dhi_letter}, the results cannot be directly applied to vehicular networks because of their unique spatial geometry in which the locations of nodes are
coupled with underlying roadways. This shortcoming of the PPP model was recognized long ago in \cite{bacc_plp}, which proposed a doubly stochastic model to capture this spatial coupling by modeling the layout of roads by a PLP and the locations of vehicles on each of these roads by independent 1D PPPs. This spatial model, referred to as the PLCP, has recently been employed in several works on vehicular communication networks \cite{vishnuJ2, vishnuL1, vishnuJ4,  baccchoi}. The reliability of the typical link in a vehicular ad hoc network (VANET) has been studied in \cite{vishnuL1}. The signal-to-interference ratio (SIR)-based coverage analysis of the typical receiver in a vehicular network has been presented in \cite{vishnuJ2, baccchoi}. In \cite{vishnuJ4}, the coverage analysis of a vehicular network in the presence of shadowing has been studied by leveraging the asymptotic characteristics of the PLCP. In the same work, the authors have also provided the first moment of the load on the MBSs and computed the rate coverage of the typical receiver in a multi-tier vehicular network. However, the characterization of the PMF of the load on the MBSs due to vehicular users modeled by a PLCP remains a key open problem and is the main focus of this letter. More details of our contribution are provided next.

{\em Contributions.} In this letter, we consider a single-tier vehicular network in which the locations of vehicular users are modeled by PLCP and the locations of MBSs serving these users are modeled by a homogeneous 2D PPP. For this setup, we first compute the Laplace transform of the distribution of the total chord length in the typical Poisson Voronoi (PV) cell. We then derive the exact PMF of the load on the typical MBS in terms of the derivative of the Laplace transform of the chord length distribution. Since the evaluation of the exact expression is not computationally efficient, we also propose a simple yet accurate approximation of the PMF of the load on the typical MBS. Using a similar approach, we also derive the PMF of the load on the MBS that serves the typical vehicular user. We also discuss the application of this result in the computation of the rate coverage of the typical receiver in the network.

\section{System Model}\label{sec:sysmod}

We consider a vehicular network in which the locations of vehicular users are modeled by a motion-invariant PLCP, as illustrated in Fig. \ref{fig:sysmod}. First, we model the physical layout of the roads by a motion-invariant PLP $\Phi_l$ with line density $\mu_l$. We denote the corresponding point process of $\Phi_l$ in the representation space $\mathcal{C} : {\mathbb{R}}^+ \times [0, 2\pi)$ by $\Phi_\mathcal{C}$ which is a homogeneous 2D PPP with density $ \lambda_l = \frac{\mu_l}{\pi} $ \cite{stoyan, vishnuJ2}. Further, we model the locations of vehicular users on each line $L_i \in \Phi_l$ by a homogeneous 1D PPP $\Psi_{L_i}$ with density $\lambda_v$, thereby forming a motion-invariant PLCP $\Phi_v$. We assume that the vehicular users are served by cellular MBSs, whose locations are modeled by a homogeneous 2D PPP $\Phi_b$ with density $\lambda_b$.

We assume that all the MBSs have the same transmit power. Further, we assume a maximum average received power based association policy and as a result, the association regions of the MBSs are the cells of the Poisson Voronoi tessellation (PVT) formed by $\Phi_b$. For completeness, the PV cell for ${\rm x} \in \Phi_b$ is defined as
\begin{align}
V_{\rm x} = \{ {\rm y} \in {\mathbb{R}}^2 : \| {\rm y} - {\rm x} \| \leq \|{\rm y} - {\rm x}'\|, \forall {\rm x}' \in \Phi_b  \}.
\end{align}
The cell $V_o$ formed by the \textit{typical MBS} at the origin is termed the \textit{typical cell}.
Our goal is to characterize the number of vehicular users located within the association region of the MBS. Specifically, we will derive the PMF of the load on the typical MBS and the serving MBS of the typical vehicular user, which will be referred to as the \textit{tagged MBS}.

\begin{figure
	\centering
	\includegraphics[width=0.36\textwidth]{./Chetlur_WCL2020-0849_fig1}
	\caption{{Illustration of the system model.}}
	\label{fig:sysmod}
\end{figure}

\section{Load on the typical MBS}
{In this section, we will derive the PMF of the load on the	typical MBS due to vehicular users, which is one of the key inputs to the network dimensioning process and is also useful in the characterization of volume of traffic at the MBS.}

\subsection{Chord length distribution}

As the locations of vehicular users are confined to roadways, the number of users within the typical cell directly depends on the length of the lines inside the typical cell. So, in order to derive the PMF of the load, we will first focus on the characterization of the total length of the lines of the PLP inside the typical cell, which will be referred to as the \textit{total chord length} and denoted by $W$.
The distribution of the length of the chord segment of the line passing through the typical PV cell has already been provided in \cite[Section 4]{muche}. Therefore, we state this result in the following lemma without a proof.

\begin{lemma}\label{lem:pdfc}
	The probability density function (PDF) of the chord length of the line intersecting the typical cell $V_o$ is $ f_C (c) $
	\begin{align}
	 &= \frac{\pi}{2} \lambda_b^{\frac{3}{2}} \int_0^\infty \mspace{-12mu} \int_0^\pi \mspace{-6mu} \tau \bigg[ \lambda_b \left(\frac{\partial K_{\tau, \alpha} }{\partial c } \right)^2 \scalebox{2}[1]{-}  \frac{\partial^2 K_{\tau, \alpha} }{\partial c^2 } \bigg] e^{- \lambda_b K_{\tau, \alpha} } {\rm d} \alpha {\rm d} \tau,\\
	\notag &{where} \quad	K_{\tau, \alpha} = 2 \pi \tau^2  - 2 \pi \tau c \cos \alpha - \tau^2 \left(\alpha - 0.5 \sin 2\alpha \right)\\
	& \quad + \pi c^2 - (\tau^2 -2 \tau c \cos \alpha + c^2) \left(\phi -  0.5 \sin 2\phi \right),\\
	 &{and} \hspace{1em} \phi = \arccos \left[ \frac{2 c^2 - 2 \tau c \cos \alpha }{2 c \sqrt{\tau^2 - 2 \tau c \cos \alpha + c^2 }}\right].
	\end{align}
\end{lemma}
Note that the number of lines that intersect the typical cell is random. Using the properties of the PLP and the PDF of the chord length given in the above lemma, we will now characterize the total chord length distribution in the typical cell in the following lemma.

\begin{lemma}\label{lem:lapl_chordlength}
	The Laplace transform of the total chord length distribution is
	\begin{align}
	\notag &\mathcal{L}_W(s) = \frac{1}{1163} \lambda_b^{3.7903} \int_0^\infty \mspace{-12mu} u^{6.5806}\exp \bigg[ - \lambda_l u  \\
	&  \times  \int_0^\infty \mspace{-12mu} \left(1 - e^{-sc} \right)  f_C (c) {\rm d} c -0.1165 (\sqrt{\lambda_b} u)^{2.3361}\bigg] {\rm d} u,
	\end{align}
	where $f_C(c)$ is given in Lemma \ref{lem:pdfc}.
\end{lemma}
\begin{IEEEproof}
	First, we compute the Laplace transform of the distribution of the total chord length, denoted by $W$, conditioned on the number of lines $N_0$ that intersect the typical cell $V_o$. Thus, we have
	\begin{align}
	\notag &\mathcal{L}_W(s | N_0) = {\mathbb{E}} \left[ e^{-s \sum_{i=1}^{n_0}  C_i } \mid N_0 = n_0\right] \\
	 & \stackrel{(a)}{=} \prod_{i=1}^{n_0} {\mathbb{E}} \left[ e^{-s C_i}  \mid N_0 = n_0\right]
	 \stackrel{(b)}{=}  \left[ \int_0^\infty \mspace{-12mu} e^{-s c}f_C (c) {\rm d} c \right]^{n_0},
	\end{align}
	where $C_i$ denotes the chord length of the line $L_i$ that intersects the typical cell, and (a) and (b) follow from the independent and identical distribution of chord lengths of different lines intersecting the typical cell.

	In order to calculate the Laplace transform of the total chord length, we must now compute the expectation of the above expression w.r.t. $N_0$. From the properties of the PLP, it follows that the number of lines of a motion-invariant PLP with line density $\mu$ that intersect a convex body $K$ follows a Poisson distribution with mean $\frac{\mu}{\pi} S_K$, where $S_K$ denotes the perimeter of $K$ \cite[Section II]{vishnuJ2}, \cite[Example 8.2]{stoyan}. Note that the perimeter of the typical cell is also random and its empirical distribution is well-known in the literature \cite{tanemura}. Thus, conditioned on the perimeter $U$ of the typical cell $V_o$, the Laplace transform of the chord length distribution is obtained as
	\begin{align}
	\notag \mathcal{L}_W&(s | U )   = {\mathbb{E}}_{N_0} \left[ \left( \int_0^\infty \exp(-s c) f_C (c) {\rm d} c \right)^{n_0} \biggm| U\right] \\

	\notag & = \sum_{n_0 = 0}^\infty \frac{ \exp\left( - {\lambda_l u}\right) (\lambda_l u)^{n_0}}{n_0!}\left( \int_0^\infty \mspace{-6mu} \exp(-s c) f_C (c) {\rm d} c \right)^{n_0} \\
	\label{eq:cond_LW} & \stackrel{(a)}{=} \exp \left[ - \lambda_l u  + \lambda_l u \int_0^\infty \exp(-s c) f_C (c) {\rm d} c \right] ,
	\end{align}
	where (a) follows from the Taylor series formula for the exponential function. From \cite{tanemura}, the PDF of the perimeter of the typical cell is given by
	\begin{align}
	&f_U(u) = \frac{\sqrt{\lambda_b}}{4} g\left(2.33609, 2.97006, 7.58060, \frac{\sqrt{\lambda_b}u}{4} \right),\\
	&\text{where }g(a,b,c,x) = {a b^{c/a}}{(\Gamma\left(c/a\right))^{-1}} x^{c-1} e^{ - b x^a}.
	\end{align}

	Upon computing the expectation of the expression in \eqref{eq:cond_LW} w.r.t. $U$, we obtain the final expression.
\end{IEEEproof}

\subsection{Load distribution}

Using the Laplace transform of the distribution of the total chord length, we will now derive the PMF of the load on the typical MBS in the following theorem.

\begin{thm}\label{thm:load_typ_exact}
	The PMF of the load on the typical MBS is
	\begin{align}\label{eq:pmf_load_exact}
	\mathbb{P} (M = m) = \frac{(-\lambda_v)^m}{m!} \left[ \frac{\partial^m }{\partial s^m} \mathcal{L}_W(s) \right]_{s= \lambda_v},
	\end{align}
	where $\mathcal{L}_W(s)$ is given in Lemma \ref{lem:lapl_chordlength}.
\end{thm}
\begin{IEEEproof}
	Recall that the distribution of points on each line of the PLCP $\Phi_v$ follows a PPP with density $\lambda_v$. Therefore, due to the independent spatial distribution of points on different lines, the total number of points of $\Phi_v$ that lie inside the typical cell $V_o$ follows a Poisson distribution with mean $\lambda_v W$. Thus, conditioned on $W$, the PMF of the load is given by
	\begin{align}
	\mathbb{P} (M = m \mid W) = \frac{\exp(- \lambda_v w) (\lambda_v w )^m }{m!}.
	\end{align}
	By taking the expectation of the above expression w.r.t. $W$, we obtain the PMF of the load as
	\begin{align}
	\mathbb{P} (M =m ) = \int_0^\infty \frac{\exp(- \lambda_v w) (\lambda_v w )^m }{m!} f_W (w) {\rm d} w.
	\end{align}
	Note that it is not straightforward to obtain the PDF of $W$. However, the above expression can be expressed in terms of the Laplace transform of $W$ as $\mathbb{P} (M =m )$
		\begin{align}
		\notag &  = \int_0^\infty \frac{(-\lambda_v)^m}{m!}  \left[ \frac{\partial^m e^{-sw}}{\partial s^m}\right]_{s= \lambda_v}  \mspace{-24mu} f_W (w) {\rm d} w \stackrel{(a)}{=} \frac{(\scalebox{2}[1]{-}\lambda_v)^m}{m!} \\
		\notag  & \times \left[ \frac{\partial^m }{\partial s^m} \int_0^\infty \mspace{-12mu} e^{-s w} f_W (w) {\rm d} w \right]_{s= \lambda_v} \mspace{-24mu} \stackrel{(b)}{=} \frac{(\scalebox{2}[1]{-}\lambda_v)^m}{m!} \left[ \frac{\partial^m \mathcal{L}_W(s)}{\partial s^m}  \right]_{s= \lambda_v}\mspace{-12mu},
		\end{align}
	where (a) follows from Leibniz's rule for derivative of integral, and (b) follows from the definition of Laplace transform. {Note that the PMF for the case of $m=0$ yields the probability that there are no users in the typical cell and is simply obtained using $\mathcal{L}_W(s)$ (i.e., without computing its derivative).}

\end{IEEEproof}

\subsection{Approximation of the PMF}
While the derivation of the exact PMF of the load concretely captures the subtle dependence of load on other sources of randomness, such as the perimeter of the typical cell, the resulting expression is not very conducive for numerical evaluation. Consequently, it is much more difficult to evaluate other metrics that utilize this result. Motivated by this, we propose an accurate approximation of the PMF of the load using the following assumption.

\begin{assumption}\label{assum:disc}
	We assume that the distribution of number of users in the typical cell centered at the origin is the same as that of the number of users inside a disc $B(o, R_t)$ centered at the origin with radius $R_t$ such that the area of the disc is equal to that of the typical cell, i.e. $\pi R_t^2 = | V_o|$
\end{assumption}

{This assumption is inspired by the asymptotic result that the large PV cells are circular \cite[Theorem 4]{praful}. As will be shown in Section \ref{sec:results}, the PMF of the load computed using this approximation is remarkably accurate in the non-asymptotic regime as well.} Using this assumption, we derive the approximate PMF of the load in the following theorem.
\begin{thm}\label{thm:load_typ_approx}
	Under Assumption \ref{assum:disc}, the PMF of the load on the typical MBS is given by
	\begin{multline}
	\mathbb{P} ( M =m ) \approx  \frac{(-\lambda_v)^m}{m!}\Bigg[\frac{\partial^m}{\partial s ^m} \int_0^\infty  \exp\bigg[ - 2 \pi \lambda_l r_t\\
	\quad +2 \pi \lambda_l \int_0^{r_t}  e^{-2s \sqrt{r_t^2 - \rho^2} }  {\rm d} \rho \bigg] f_{R_t}(r_t) {\rm d} r_t \Bigg]_{s = \lambda_v} .
	\end{multline}
\end{thm}
\begin{IEEEproof}
	As we have shown earlier, the first step in the computation of the PMF of the load on the typical MBS is the characterization of the total chord length. For a line $L$ at a perpendicular distance $\rho < R_t$ from the origin, the length of the chord that intersects $ B(o,R_t) $ is $2 \sqrt{R_t^2 - \rho^2}$. From the properties of PLP, we also know that the number of lines of $\Phi_l$ that intersect the disc $B(o,R_t)$ follows a Poisson distribution with mean $ 2 \pi \lambda_l R_t $. Thus, conditioned on the radius of the disc $R_t$, the Laplace transform of the total chord length $W$ in $B(o, R_t)$ can be computed as
	\begin{align}
	\notag &\mathcal{L}_W(s | R)  \approx {\mathbb{E}} \left[ e^{ -s \sum\limits_{L \in \Phi_l}  \nu_1 (B(o,r_t) \cap L) }\right] = {\mathbb{E}} \bigg[ \mspace{-12mu} \prod_{\substack{(\rho, \theta) \in \Phi_\mathcal{C} \\ \rho \leq r_t}} \mspace{-20mu} e^{ - 2 s \sqrt{r_t^2 - \rho^2}}\bigg]\\
	\notag &  \stackrel{(a)}{=} \exp\left[ - 2 \pi \lambda_l \int_0^{r_t} 1 - \exp(-2s \sqrt{r_t^2 - \rho^2})  {\rm d} \rho \right],
	\end{align}
	where (a) follows from the PGFL of the PPP $\Phi_\mathcal{C}$. From \cite{tanemura}, the PDF of the area of the typical cell $Z$ is given by
	\begin{align}
	f_Z(z) = \lambda_b g(1.07950, 3.03226, 3.31122, \lambda_b z ).
	\end{align}

	Since $R_t = \sqrt{|V_o|/ \pi}$, the PDF of the radius of the disc is $f_{R_t}(r_t) = 2 \pi r_t f_Z(\pi r_t^2)$.
	Thus, the Laplace transform of the total chord length distribution is
	\begin{align}\label{eq:lapW_approx}
	  &\mathcal{L}_W(s) \approx \int_0^\infty \mspace{-12mu}  \exp\bigg[ \scalebox{2}[1]{-} 2 \pi \lambda_l \int_0^{r_t}  \mspace{-6mu} 1\scalebox{2}[1]{-} e^{\scalebox{2}[1]{-}2s \sqrt{r_t^2 \scalebox{2}[1]{-} \rho^2}}  {\rm d} \rho \bigg] f_{R_t}(r_t) {\rm d} r_t.
	\end{align}
	From Theorem \ref{thm:load_typ_exact}, we know that the PMF of the load can be expressed in terms of the derivative of $\mathcal{L}_W(s)$. Hence, upon substituting \eqref{eq:lapW_approx} in \eqref{eq:pmf_load_exact}, we obtain the approximate PMF of the load on the typical MBS.

\end{IEEEproof}

\section{Load on the tagged MBS}
In this section, we will derive the PMF of the load on the tagged MBS that serves the typical vehicular user. Without loss of generality, we translate the origin to the location of the typical vehicular user. Thus, under the Palm distribution of the PLCP, the resulting point process of vehicular users $\Phi_{v_0}$ is the superposition of $\Phi_v$, a 1D PPP $\Psi_{L_0}$ with density $\lambda_v$ along the line $L_0$ that passes through the origin, and an atom at the origin, i.e., $\Phi_{v_0} = \Phi_v \cup \Psi_{L_0} \cup \{o\}$. The line $L_0$ will henceforth be referred to as the typical line. The Voronoi cell of the tagged MBS that contains the origin is called the \textit{zero cell} and it is mathematically expressed as $ V(o) =$
\begin{align*}
 \{ {\rm y} \in {\mathbb{R}}^2 : \| {\rm y} - {\rm x} \| \leq \| {\rm y} - {\rm x}'\|,\  \|{\rm x}\| < \| {\rm x}'\|, \forall {\rm x}' \in \Phi_b  \}.
\end{align*}

{Owing to the Palm distribution of the PLCP and the statistically larger area of the zero cell, its total chord length distribution is different from that of the typical cell.} So, we first derive the Laplace transform of the distribution of the total chord length in the following lemma by using an assumption similar to that of Assumption \ref{assum:disc}.

\begin{lemma}\label{lem:lapl_chordlength_zerocell}
	The Laplace transform of the distribution of the total chord length in the zero cell is
	\begin{align}
	\notag&\mathcal{L}_W(s) \approx \int_0^\infty \int_0^\infty \exp \bigg[ -s c_0 - 2 \pi \lambda_l r_z + 2 \pi \lambda_l \\
	&\times  \int_0^{r_z} e^{-2s \sqrt{r_z^2 - \rho^2}}  {\rm d} \rho \bigg] f_{R_z}(r_z) f_{C_0}(c_0) {\rm d} c_0 {\rm d} r_z,
	\end{align}
	where $f_{R_z}(r_z) = 2 \pi^2 \lambda_b r_z^3 f_Z(r_z)$, and $f_{C_0} = \frac{4 \sqrt{\lambda_b}}{\pi} c_0 f_C(c_0)$.
\end{lemma}
\begin{IEEEproof}
	We know that the typical line $L_0$ intersects the zero cell and we denote the corresponding chord length by $C_0$. In addition to this, there exists a random number of lines of $\Phi_l$ that intersect the zero cell. Thus, we write the total chord length of the zero cell as $W = C_0 + C_1$, where $C_1$ denotes the sum of chord lengths of all lines that intersect $V(o)$ excluding the typical line. We will compute the Laplace transforms of the distribution of $C_0$ and $C_1$ separately.

	The PDF of the typical chord length of a line in an arbitrarily chosen PV cell is given in Lemma \ref{lem:pdfc}. However, we are interested in the length of the chord segment that contains the origin, which is expected to be longer than the length of the typical chord segment. Thus, using length-biased sampling, the PDF of $C_0$ is given by
	\begin{align}
	f_{C_0}(c_0)  =  \frac{c_0 f_{C}(c_0)}{{\mathbb{E}}[C]} = \frac{4 \sqrt{\lambda_b}}{\pi} c_0 f_C(c_0).
	\end{align}
	Using this result, the Laplace transform of the distribution of $C_0$ can be computed as
	\begin{align}\label{eq:lapC0}
	\mathcal{L}_{C_0}(s) = \int_0^\infty \exp(-s c_0) f_{C_0}(c_0) {\rm d} c_0.
	\end{align}

	We will now focus on the computation of the Laplace transform of the distribution of $C_1$. First, we assume that the number of vehicular users in the zero cell, excluding the users on the typical line, is the same as the number of vehicular users in a disc $B({\rm x}, R_z)$ centered at the location of the serving MBS ${\rm x}$ with radius $R_z$ such that $\pi R_z^2 = |V(o)| $. From area-biased sampling \cite{Dhi_letter}, the PDF of the area of the zero cell is
	\begin{align}
	f_{Z'}(z') = \frac{z' f_Z(z')}{{\mathbb{E}}[Z]} = \lambda_b z' f_Z(z').
	\end{align}

	Following the same procedure as in the proof of Theorem \ref{thm:load_typ_approx}, we obtain the approximate Laplace transform of the distribution of $C_1$ as $\mathcal{L}_{C_1}(s) $
	\begin{align}\label{eq:lapC1}
	 \approx \mspace{-6mu} \int_0^\infty \mspace{-12mu} \exp\bigg[ - 2 \pi \lambda_l \int_0^{r_z} \mspace{-12mu}  1-e^{-2s \sqrt{r_z^2 - \rho^2} }  {\rm d} \rho \bigg] f_{R_z}(r_z) {\rm d} r_z,
	\end{align}
	where $f_{R_z}(r_z) = 2 \pi r_z f_{Z'}(\pi r_z^2)$.

	 As $C_0$ and $C_1$ are independent, the Laplace transform of the distribution of total chord length can be computed as
	\begin{align}\label{eq:lapW_zerocell}
	\mathcal{L}_W(s) = \mathcal{L}_{C_0}(s)  \mathcal{L}_{C_1}(s).
	\end{align}
	Substituting \eqref{eq:lapC0} and \eqref{eq:lapC1} in the above equation yields the expression for the Laplace transform of the total chord length distribution in the zero cell.
	\end{IEEEproof}

\begin{thm}\label{thm:load_tagged}
	The PMF of the load on the tagged MBS is
	\begin{align}
		\notag &\mathbb{P} (M = m+1) = \frac{(-\lambda_v)^m}{m!} \int_0^\infty \mspace{-12mu} \int_0^\infty \Bigg[  \frac{\partial^m }{\partial s^m} \exp \bigg[-s c_0 \\
		\label{eq:pmf_zerocell} & \scalebox{2}[1]{-}2 \pi \lambda_l  \int_0^{r_z} \mspace{-12mu} 1\scalebox{2}[1]{-} e^{-2s \sqrt{r_z^2 - \rho^2}}  {\rm d} \rho \bigg]  \Bigg]_{s= \lambda_v} \mspace{-28mu}f_{R_z}(r_z) f_{C_0}(c_0) {\rm d} c_0 {\rm d} r_z,
	\end{align}
	where $ m = 0,\ 1,\ 2, \dots$.
\end{thm}
\begin{IEEEproof}
	The proof follows along the same lines as that of Theorem \ref{thm:load_typ_exact}.
\end{IEEEproof}

\section{Results and Discussion}

In this section, we will verify the accuracy of our analytical results and also discuss the application of these results in the computation of rate coverage of the typical vehicular user.
\begin{figure
	\centering
	\includegraphics[width=0.4\textwidth]{./Chetlur_WCL2020-0849_fig2}
	\caption{{PMF of the load on the typical MBS ($\mu_l = 5$ km$^{-1}$, $\lambda_v = 2$ nodes/km, and $\lambda_b = 1$ node/km$^2$).}}
	\label{fig:pmf_typ}
\end{figure}
\begin{figure
	\centering
	\includegraphics[width=0.4\textwidth]{./Chetlur_WCL2020-0849_fig3}
	\caption{{PMF of the load on the tagged MBS ($\mu_l = 5$ km$^{-1}$, $\lambda_v = 2$ nodes/km, and $\lambda_b = 1$ node/km$^2$).}}
	\label{fig:pmf_zero}
\end{figure}

\subsection{Numerical results}\label{sec:results}
In Fig. \ref{fig:pmf_typ}, we plot the PMF of the load on the typical cell evaluated from expressions given in Theorem \ref{thm:load_typ_approx} along with the empirical PMF of the load obtained through Monte-Carlo simulations. It is quite evident that the results derived in Theorem \ref{thm:load_typ_approx} match closely with the empirical results from simulations thereby demonstrating the accuracy of our analytical results. Similarly, the PMF of the load on the tagged MBS evaluated using the expressions in Theorem \ref{thm:load_tagged} matches closely with the simulations as shown in Fig. \ref{fig:pmf_zero}.

\subsection{Rate coverage}
In this subsection, we will focus on the characterization of rate coverage of the typical receiver in the network using the load distribution of vehicular users. We consider a single tier vehicular network as described in Section \ref{sec:sysmod}. Without loss of generality, we translate the origin to the location of the typical receiver. Thus, the resulting point process of vehicular users is $\Phi_{v_0} \equiv \Phi_v \cup \Psi_{L_0} \cup \{o\}$. We assume that all the MBSs have the same transmit power $P_t$ and the typical receiver connects to the closest MBS in the network. In this setup, the SIR measured at the typical receiver is
\begin{align}
\mathtt{SIR} = \frac{P_t H_0 R^{-\alpha}}{\sum_{{\rm x} \in \Phi_b} P_t H_{\rm x} \| {\rm x} \|^{-\alpha} },
\end{align}
where $R$ denotes the serving distance, and $H_0$ and $H_{\rm x}$ denote the fading gains of the links from the typical receiver to the serving MBS and interfering MBSs at ${\rm x}$, respectively, which are exponentially distributed with unit mean. Assuming that the bandwidth $B$ is equally shared by all the users served by the serving MBS, the achievable rate at the typical receiver is
\begin{align}\label{eq:rate}
{\mathcal{R}} = \frac{B}{M}\log_2(1 + \mathtt{SIR}),
\end{align}
where $M$ denotes the load on the serving MBS. We present the rate coverage of the typical receiver, which is defined as the probability with which the rate at the typical receiver exceeds a desired threshold $T$, in the following theorem.

\begin{thm}\label{thm:rate_coverage}
	The rate coverage of the typical receiver is $ \mathsf{R}_{\rm c} $
	\begin{align*
	 & \approx \sum_{m=1}^\infty P_m 2\pi\lambda_b \int_0^\infty \mspace{-12mu} r\exp \bigg[ \scalebox{2}[1]{-}\lambda_b \pi r^2   \scalebox{2}[1]{-}  \int_r^\infty \mspace{-6mu} \frac{2 \pi \lambda_b \gamma y}{\gamma  +  (y/r)^\alpha } {\rm d} y \bigg] {\rm d} r,
	\end{align*}
	where $\gamma = 2^{\frac{T m}{B}} -1$, and $P_m \equiv \mathbb{P}(M = m) $ is given by Theorem \ref{thm:load_tagged}.
\end{thm}
\begin{IEEEproof}
	Although $M$ and $\mathtt{SIR}$ in \eqref{eq:rate} are negatively correlated in general, we consider them to independent in the interest of analytical tractability. This is a popular assumption in the literature and its accuracy has been verified in many settings, e.g., see \cite{Dhi_letter}. Thus, the rate coverage of the typical receiver can be computed as
	\begin{align}
	\notag \mathsf{R}_{\rm c} &= \mathbb{P} ( {\mathcal{R}} > T) = \mathbb{P} \left( \frac{B}{M}\log_2 (1 + \mathtt{SIR}) > T\right)\\
	\label{eq:rc_step1} &=\sum_{m=1}^\infty \mathbb{P}(M = m) \mathbb{P} \left( \mathtt{SIR} > 2^{\frac{T m}{B}} -1 \right).
	\end{align}
	For this setup, {the probability that $\mathtt{SIR}$ at the typical receiver is greater than a desired threshold $\beta$, referred to as the \textit{coverage probability} of the typical receiver in an interference dominant setting,} is well-known in the literature \cite{DhiGanBacAnd} and is given by $\mathbb{P}(\mathtt{SIR}>\beta)= $
	\begin{align}\label{eq:coverage}
	 &2\pi\lambda_b \int_0^\infty \mspace{-12mu} r\exp \bigg[ -\lambda_b \pi r^2  -  \int_r^\infty \frac{2 \pi \lambda_b \beta y}{ \beta +   r^{-\alpha}y^\alpha } {\rm d} y \bigg] {\rm d} r.
	\end{align}
	We also know the expression for the PMF of the load on the serving MBS from Theorem \ref{thm:load_tagged}. So, upon substituting \eqref{eq:coverage} and \eqref{eq:pmf_zerocell} in \eqref{eq:rc_step1}, we obtain the expression for rate coverage.
\end{IEEEproof}

\begin{figure
	\centering
	\includegraphics[width=0.4\textwidth]{./Chetlur_WCL2020-0849_fig4}
	\caption{{Rate coverage of the typical receiver ($\mu_l = 5$ km$^{-1}$, $\lambda_v = 2$ nodes/km, $\alpha$ = 4, and $B$ = 10 MHz).}}
	\label{fig:rate_coverage}
\end{figure}

From Fig. \ref{fig:rate_coverage}, it can be noted that the rate coverage of the typical receiver evaluated using the analytical expression given in Theorem \ref{thm:rate_coverage} matches with the simulation results.

\section{Conclusion}
In this letter, we have derived the PMF of the load on cellular MBSs due to vehicular users modeled by PLCP in a single-tier network. Specifically, we characterized the total chord length distribution in a typical PV cell, which is an important result in its own right, and computed the PMF of the load on the typical MBS. As the evaluation of this expression is quite challenging, we proposed an accurate approximation for the PMF of the load on the typical MBS. Following a similar approach, we then derived the PMF of the load on the tagged MBS. Using this result, we also computed the rate coverage of the typical vehicular user in the network.

\bibliographystyle{IEEEtran}