\section{Introduction}\label{sec:introduction}}
\else
\section{Introduction}
\label{sec:introduction}
\fi

Computer vision techniques, such as object detection and semantic segmentation, are used for many aerial image analysis applications. As an example, traffic can be analyzed by using object detection to find vehicles in aerial images. Semantic segmentation is the basis for creating and updating maps \cite{girard2018aligning} and can be used for tracking city growth \cite{gomez2020spatiotemporal} or tracking deforestation \cite{lee2020classification}.

However, any error in the prediction of a computer vision model, i.e., object detection or semantic segmentation, causes errors in applications based on said prediction. Precise predictions are required, which can be achieved by exploiting domain knowledge. One fact rarely exploited in the deep learning-based analysis of aerial images is that many objects and regions can be approximated well by polygons, especially in true orthophotos. By restricting a model's prediction to a limited number of polygons, the degrees of freedom in the prediction are reduced. Arbitrary region shapes, e.g., circle-shaped cars which make no sense, are no longer feasible predictions.

We created BSPSegNet\cite{gritzner2021semantic}, a model for aerial image segmentation which predicts polygons as shown in Fig. \ref{fig:teaser}.  The model partitions the input image into $8\times 8$ blocks. For each block, a binary space partitioning (BSP) tree \cite{fuchs1980visible} with depth 2 is predicted, i.e., every block is partitioned into up to four segments with a separate class prediction for each segment. The parameters of the inner nodes of the BSP tree encode parameters of a line partitioning a region into two subregions, while the parameters of the leaf nodes are class logits. The lines encoded by the inner nodes together with the block boundaries form polygons. As we demonstrated in our earlier paper, the ground truth of common aerial image datasets can be encoded extremely accurately this way ($99\%$ accuracy and $99\%$ mean Intersection-over-Union (mIoU)).

\begin{figure}[t]
	\begin{center}
		\includegraphics[width=\columnwidth]{figures/teaser.pdf}
	\end{center}
	\caption{Overview of our approach: our model predicts binary space partitioning (BSP) trees from aerial images. The inner nodes of such a tree define the shape of regions, while the leaf nodes define the content of each region. The leaf nodes effectively map shapes to classes. BSP trees can be rendered in a differentiable way into a full segmentation thus enabling end-to-end model training. Our proposed refinements improve gradient computations in the BSP renderer and parts of the model. Additionally, a novel loss function improves the predicted inner node parameters, i.e., the predicted shapes. Furthermore, we extended the approach to predict multiple trees at the same time in order to enable class-specific shape predictions.}
	\label{fig:teaser}
\end{figure}

However, despite reaching state-of-the-art performance in semantic segmentation, BSPSegNet still does not offer proven improvements over existing models. In this article, we aim to improve BSPSegNet such that it performs better than the state-of-the-art for segmentation.
We present multiple contributions:
\begin{enumerate}
    \item We replace BSPSegNet's two feature decoders with an architecture using residual connections and modify the model's differentiable BSP tree renderer to avoid using sigmoid activations. BSPSegNet tries to push the output of these activations close to $0$ or $1$. However, in these value regions the derivative of the sigmoid function is almost $0$. Vanishing gradients are avoided due to these modifications.
    \item We introduce a novel loss function punishing BSP trees which create partitions containing multiple different classes (according to the ground truth). This loss specifically affects only the inner nodes of the BSP trees, instead of all the nodes as categorical cross entropy does, thus improving the shapes of the predicted partitions.
    \item We expand BSPSegNet's prediction s.t. it is able to predict multiple BSP trees for each $8\times 8$ block, thus enabling class-specific segmentations for more precise shape predictions.
\end{enumerate}
Due to the last modification, we rename our model to SegForestNet. SegForestNet is able to reach state-of-the-art performance while using up to $60\%$ less parameters, i.e., it is less prone to overfitting to a small dataset. Small dataset size is a common problem in aerial image analysis.

With these modifications SegForestNet can be trained end-to-end without the two phase training previously recommended in \cite{gritzner2021semantic}. In the first phase, an autoencoder phase, a mapping of the ground truth to BSP trees was learned. In the second phase, this new ground truth representation was used to learn the actual semantic segmentation without using the differentiable BSP renderer and thus avoiding its effects on gradients. Our modifications allow the faster and more efficient training of SegForestNet without the autoencoder phase.

In addition to the above contributions, we also experimented with different signed distance functions (SDFs) and tree structures. SDFs are used by the inner node nodes in order to partition a segment into subsegments. By trying different SDFs, we were able to generalize the shape of the partitioning boundary from lines to other shapes such as circles. Also, our approach generalizes to different tree structures such as quadtrees \cite{finkel1974quad} and k-d trees \cite{bentley1975multidimensional}. However, we found that neither of these two changes provide any benefit for the semantic segmentation of aerial images. Since we believe that these ideas may be useful in other domains with more organic shapes, e.g., the segmentation of images captured by cameras of self-driving cars, we still present these experiments here so that they become more easily reproducible.

In this article's next section we will discuss related work, followed by a description of SegForestNet. We then evaluate our contributions by comparing SegForestNet on several datasets to multiple state-of-the-art semantic segmentation models. We finally conclude with a summary.

Our implementation can be found on GitHub\footnote{\url{https://github.com/gritzner/SegForestNet}}.

\section{Related Work}
Modern semantic segmentation models rely on deep learning and an encoder-decoder structure \cite{papandreou1502weakly,yu2015multi,jegou2017one,chen2017rethinking,wu2019fastfcn,PesYan2022}. An encoder is used to map an input image to a feature map with a reduced spatial resolution but a much higher number of channels/features. While a large number of features is computed, the spatial information is bottlenecked. This implies that downsampling, e.g., through max-pooling or strided convolutions, is performed. Using such a feature map, a decoder predicts the desired semantic segmentation. To get a segmentation of the same resolution as the input image, the decoder uses bilinear sampling and/or transpose convolutions (also called deconvolutions) in order to upsample the spatial resolution of the feature map back to the initial input resolution.

U-Net \cite{ronneberger2015u} is a popular model following this paradigm. Though designed for medical image analysis, it is also used for other types of data, including even 3D data \cite{cciccek20163d}. U-Net's encoder and decoder are almost symmetrical and skip connections at each spatial resolution are used to allow predictions with high frequency details. Other models, such as Fully Convolutional Networks for Semantic Segmentation (FCN) \cite{long2015fully}, have an asymmetrical architecture with a significantly more sophisticated encoder than decoder.
They rely on existing classification models, such as VGG \cite{simonyan2014very}, ResNet \cite{he2016deep}, MobileNetv2 \cite{sandler2018mobilenetv2} or Xception
\cite{chollet2017xception}, as encoder (also called backbone or feature extractors in these cases) by removing the final classification layers.
Instead, semantic segmentation models append a novel decoder to these encoders. The complexity of these decoders varies. While FCN uses a rather simple decoder, DeepLabv3+ \cite{chen2018encoder} uses a more sophisticated decoder and also adds atrous spatial pyramid pooling \cite{chen2017deeplab} to the encoder. By doing so, DeepLabv3+ is able to account for context information at different scales. Other modifications to the basic encoder-decoder idea are adding parallel branches to process specific features such as edges/boundaries \cite{takikawa2019gated}, processing videos instead of images \cite{zhu2019improving}, explicitly modelling long-range spatial and channel relations \cite{mou2019relation}, separating the image into fore- and background \cite{zheng2020foreground}, or point-wise affinity propagation to handle imbalanced class distributions with an unusually high number of background pixels \cite{li2021pointflow}.

More recent research also works on more complex problems, e.g., instance segmentation  \cite{he2017mask} and panoptic segmentation \cite{kirillov2019panoptic,xiong2019upsnet,yang2019deeperlab,yang2020sognet}. In instance segmentation pixel masks of detected objects, such as cars or persons, are predicted. Panoptic segmentation is a combination of semantic segmentation and instance segmentation: a class is predicted for each pixel in the input image while for all pixels of countable objects an additional instance ID is predicted. Models dealing with these more complex tasks still usually require a semantic segmentation as a component to compute the final prediction. AdaptIS \cite{sofiiuk2019adaptis} is such a model. When performing instance segmentation, it uses U-Net as a backbone while using ResNet and DeepLabv3+ when performing panoptic segmentation.

Our previous model, BSPSegNet \cite{gritzner2021semantic}, predicts BSP trees to perform a semantic segmentation. These trees are usually used in computer graphics \cite{sanglard2019game} rather than computer vision. Each inner node of a BSP tree partitions a given region into two segments while the leaf nodes define each segment's content, i.e., its class in the case of a semantic segmentation. Since BSP trees are resolution-indepedent, BSPSegNet's decoder does not have to perform any kind of upsampling. Furthermore, BSPSSegNet is the only segmentation model which inherently semantically separates features into shape and content features. Similar works using computer graphics techniques for computer vision are
PointRend \cite{kirillov2020pointrend} and
models reconstructing 3D meshes
\cite{tatarchenko2017octree,wang2018pixel2mesh,mescheder2019occupancy,gkioxari2019mesh}. PointRend iteratively refines a coarse instance segmentation into a finer representation. BSPnet \cite{chen2020bsp} also uses BSP, but instead of a segmentation, it tries to reconstruct polygonal 3D models. Additionally, rather than predicting a BSP tree hierarchy, it instead predicts a set of BSP planes.

As mentioned before, BSP trees are a data structure more common in computer graphics than in computer vision. Other commonly used data structures in computer graphics and other applications, e.g., data compression, are $k$-d trees \cite{bentley1975multidimensional}, quadtrees \cite{finkel1974quad}, and bounding volume hierarchies \cite{ericson2004real}. They differ from BSP trees in that they perform a partitioning using an axis-aligned structure which makes slopes expensive to encode compared to BSP trees.

\section{Semantic Segmentation through Spatial Partitioning}
We will first briefly describe BSPSegNet's architecture in this section in order to be able to precisely discuss our contributions. The third subsection presents the first contribution mentioned in the introduction, which reduces the risk of vanishing gradients, while the fourth subsection presents our second contribution, the novel loss function. The fifth subsection then discusses the changes necessary to enable to the prediction of multiple partitioning trees simultaneously, our third contribution. In a final subsection we show how SegForestNet's differentiable BSP rendering generalizes to different signed distance functions (SDFs) and tree structures.

\subsection{BSPSegNet}
Fig. \ref{fig:model} (top) shows a common encoder-decoder architecture used by state-of-the-art models. An encoder maps an input image to a feature map with a bottlenecked spatial resolution. This feature map is then decoded into a semantic segmentation at the same spatial resolution as the input image. Skip connections at different spatial resolutions help to maintain high frequency details.

\begin{figure}
	\begin{center}
		\includegraphics[width=\columnwidth]{figures/models.pdf}
	\end{center}
	\caption{A comparison of state-of-the-art models (top) and BSPSegNet/SegForestNet (bottom; ours). All models use an encoder-decoder architecture, however, our models semantically splits the feature map (hexagon) into shape and content features. These are decoded into the inner nodes and leaf nodes of BSP trees respectively (see Fig. \ref{fig:bsptree}). Our models also needs a differentiable BSP renderer to enable end-to-end training. The renderer is a fixed function without learnable parameters.}
	\label{fig:model}
\end{figure}

The bottom of Fig. \ref{fig:model} shows BSPSegNet in comparison. As other state-of-the-art models, it uses an encoder-decoder architecture. However, it uses two separate decoders. The feature map computed by the encoder is split along the feature dimension into two separate maps: one for shape features (blue) and one for content features (orange). Each feature map is processed by a different decoder. While the decoders share the same architecture, they both have their own unique weights. The shape features are decoded into the inner nodes of the BSP trees, which represent the shape of the regions in the final segmentation, as shown in Fig. \ref{fig:bsptree}. The content features are decoded into the leaf nodes of the BSP trees, which represent each region's class logits.

\begin{figure}
	\begin{center}
		\includegraphics[width=\columnwidth]{figures/bsptree.pdf}
	\end{center}
	\caption{A partitioning of a square region (right) defined by a BSP tree (center). Shape features are decoded
	into the parameters of the inner nodes (blue), which define lines (green) creating the partitioning. Content features are decoded into the parameters of the leaf nodes (orange), which are the class logits predicted for each partition.}
	\label{fig:bsptree}
\end{figure}

Like in other semantic segmenation models, BSPSegNet uses the feature extraction part of a classification model like MobileNetv2 or Xception as an encoder. By removing pooling operations or reducing the stride of convolution operations, the encoder is modified such that the input image is downsampled by a factor of $8$ along each spatial dimension. To keep in line with classification model design, the earlier downsampling operations are kept while the later ones are removed. This enables more efficient inference and training, as the spatial resolution is reduced early. Furthermore, the number of features in  the final feature map computed by the encoder is reduced such that there is a bottleneck.
The shape feature map has fewer features than the number of parameters required for the inner nodes, while the content feature map has fewer features than the parameters required for the leaf nodes.

BSPSegNet subdivides the input image into $8\times 8$ blocks (the same as the downsampling factor in the encoder) and predicts a separate BSP tree for each block, i.e., the spatial resolution of the feature map is such that there is exactly one spatial unit per block. The two decoders keep the spatial resolution and only modify the number of features. This is done by concatenating three blocks, each consisting of a $1\times 1$ convolution outputting $256$ features, batch normalization and a LeakyReLU activation. These blocks are followed by one final $1\times 1$ convolution to predict the final BSP tree parameters.

\subsection{BSPSegNet's Differentiable BSP Renderer}
The differentiable BSP renderer is based on the idea of a region map $\textbf{R}$. For each pixel $\textbf{p} = (x,y)$, the entry $\textbf{R}_{iyx}$ is the probability that $\textbf{p}$ is part of the $i$-th region (leaf node). This implies $\sum_{i=1}^k
\textbf{R}_{iyx} = 1$ for all pixels $\textbf{p}$ with $k$ equal to the number of leaf nodes.

Initially, BSPSegNet starts with all $\textbf{R}_{iyx}$ set to $1$ and then updates the region map for every inner node. Each inner node consists of three parameters, two parameters for the normal vector $\textbf{n}$ and one parameter $d$. These parameters together define a line in the two-dimensional image space. BSPSegNet does neither normalize $\textbf{n}$, nor does it enforce $|\textbf{n}| = 1$ through any kind of loss. Rather, it lets the network learn to predict an appropriately scaled $d$. For every inner node and pixel $\textbf{p}$, the signed distance function
$ f(\textbf{p}) = \textbf{n}\cdot\textbf{p} - d $
is computed. The equations
\begin{equation}
    g(\textbf{p}) = \sigma(\lambda_1 \cdot f(\textbf{p}))
    \label{eq:sigmoid}
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} \cdot \lambda_2 \cdot g(\textbf{p})
    \label{eq:rleft}
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} \cdot \lambda_2 \cdot (1 - g(\textbf{p}))
    \label{eq:rright}
\end{equation}
with the sigmoid function $\sigma$ are then used to update the region map $\textbf{R}$. For each inner node, equation \ref{eq:rleft} is used for all leaf nodes $i$ reachable from the left child node while equation \ref{eq:rright} is used for all $i$ reachable by the right child node.

Depending on where the pixel $\textbf{p}$ lies, the function $g$ either approaches $0$ or $1$. It approaches $1$ if $\textbf{p}$ is reachable via the left child node and it approaches $0$ if $\textbf{p}$ is reachable via the right child node. Therefore, $\textbf{R}_{iyx}$ is multiplied with a value close to $0$ for those regions $i$ which $\textbf{p}$ does not belong to. Thus, eventually, all $\textbf{R}_{iyx}$ for $i$ to which $\textbf{p}$ does not belong will be close to $0$. The one remaining entry $\textbf{R}_{iyx}$ will be close to $\lambda_2^D$, where $D$ is the depth of the BSP tree.

After iterating over all inner nodes, a Softmax operation across the region dimension $i$ is applied. This ensures that each $\textbf{R}_{iyx}$ represents the aforementioned probability. The hyperparameter $\lambda_2$ controls how close the respective final entries of $\textbf{R}$ will be to $0$ or $1$, while $\lambda_1$ controls how large $|f(\textbf{p})|$ must be in order for $\textbf{p}$ to be assigned to one specific child node ($g(\textbf{p})$ close to $0$ or $1$) instead of the model expressing uncertainty about which child region $\textbf{p}$ belongs to ($g(\textbf{p}) \approx 0.5$).

The final output $h$ of the BSP renderer is
\begin{equation}
    h(\textbf{p}) = \sum_{i=1}^k \textbf{R}_{iyx} \cdot \textbf{v}_i
    \label{eq:finalprediction}
\end{equation}
with the vectors $\textbf{v}_i$ of class logits, one vector for each leaf node. The function $h$ can be computed in parallel for all pixels $\textbf{p}$. BSP trees are resolution-independent, therefore the final output resolution is determined by how densely $h$ is sampled, i.e., for how many pixels $\textbf{p}$ the function is computed. This is also the reason why BSPSegNet, in contrast to other segmentation models, does not have to perform any kind of upsampling in its decoders. By sampling $h$ at the same resolution as the input image, a pixelwise semantic segmentation can be computed like for any other state-of-the-art model.

\subsection{Refined Gradients}
As mentioned in the introduction, a two-phase training process is recommended for the original BSPSegNet. In the second phase, the differentiable BSP renderer is skipped altogether. We hypothesize that insufficient gradients are responsible: by pushing the computation of $g$ (Eq. \ref{eq:sigmoid}) into value ranges where its derivative is almost $0$ gradients vanish. Additionally, the two feature decoders do not use residual connections, which would help gradients to propagate through the model. Our first contribution aims to improve on these two points in order to avoid the two-phase training process and to enable the use of decoders with more layers which can learn more complex mapping functions from features to BSP tree parameters. The goal is to enable true end-to-end training without the need to learn a mapping to an intermediate representation. This way, training is faster, allowing for more iterations during development, e.g., for hyperparameter optimization, or saving energy and thus reducing costs.

\begin{figure}
	\begin{center}
		\includegraphics[width=\columnwidth]{figures/decoder.pdf}
	\end{center}
	\caption{Our refined decoder architecture using residual blocks. The initial and the last convolution change the number of features $F$. The number of features are set such that there is a bottleneck with regard to the tree parameters, i.e., $F_{encoder} < F_{tree}$.}
	\label{fig:decoder}
\end{figure}

The first part of this contribution is a new feature decoder architecture. Again, we use the same architecture for both decoders, but they still do not share weights. Our new architecture, shown in Fig. \ref{fig:decoder}, starts and ends with a $1\times 1$ convolution, the first one maps the bottlenecked feature dimension coming from the encoder to an intermediate number of features, while the last convolution maps the intermediate features to the final number of tree parameters. In between these two convolutions are several residual blocks. Each block consists of a depthwise $3\times 3$ convolution with zero padding and a $1\times 1$ convolution. The number of features stays constant throughout the residual blocks. Each convolution in the decoders, except for the very last one mapping to the tree parameters, is followed by a batch normalization layer \cite{ioffe2015batch} and a LeakyReLU activation \cite{xu2015empirical}.
This architecture not only uses residual connections for improved gradient propagation, but it also includes spatial context for each block, i.e., for each BSP tree, instead of solely relying on the receptive field of the encoder.

The second part of our first contribution is an updated region map $\textbf{R}$ computation (Eqs. \ref{eq:sigmoid}-\ref{eq:rright}). We start by setting all $\textbf{R}_{iyx}$ to $0$ initially instead of $1$. For each inner node, we use the equations
\begin{equation}
    g(\textbf{p}) = \lambda \cdot f(\textbf{p})
    \label{eq:nosigmoid}
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} + \text{ReLU}(g(\textbf{p}))
    \label{eq:rleft2}
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} + \text{ReLU}(-g(\textbf{p}))
    \label{eq:rright2}
\end{equation}
for updating $\textbf{R}$ where Eq. \ref{eq:rleft2} replaces Eq. \ref{eq:rleft} and Eq. \ref{eq:rright2} replaces Eq. \ref{eq:rright}. After processing all inner nodes, there still is a Softmax operation across the region dimension $i$. Using our Eqs. \ref{eq:nosigmoid}-\ref{eq:rright2}, we avoid vanishing gradients by replacing the sigmoid function with ReLU. This necessitates moving from multiplications in Eqs. \ref{eq:rleft} and \ref{eq:rright} to additions in Eqs. \ref{eq:rleft2} and \ref{eq:rright2}. As a nice side effect, we also have only one hyperparameter $\lambda$ to optimize instead of two, thus simplifying the model. Our $\lambda$ combines the meaning of the previous $\lambda_1$ and $\lambda_2$.

\subsection{Region Map-Specific Loss}
Cross-entropy, the most common loss function used for semantic segmentation, is applied pixelwise and therefore gives the model to be trained a signal for each mispredicted pixel. However, from the point-of-view of BSPSegNet and SegForestNet, there is no distinction between wrong class logits in a leaf node and wrong line parameters in an inner node in this loss signal. We therefore present a novel loss function, our second contribution, which specifically punishes the later, i.e., wrong line parameters in inner nodes. We do so by formulating three additional losses which we combine with cross-entropy into a new SegForestNet-specific loss function. Each new loss is designed to punish a specific undesirable trait of the region map $\textbf{R}$ which may occur in the model's predictions.

First, as an intermediate result used later by two of the new losses, we compute
\begin{equation}
    \textbf{Y}_\textbf{B}^i = \sum_{(x,y) \in \textbf{B}} \textbf{Y}_{yx} \cdot \textbf{R}_{iyx}
    \label{eq:yblockregion}
\end{equation}
for every $8\times 8$ pixel block $\textbf{B}$ and region/leaf node $i$ where every $\textbf{Y}_{yx}$ is a one-hot class vector. $\textbf{Y}_\textbf{B}^i$ is therefore a weighted vector containing the number of pixels per class in region $i$ of block $\textbf{B}$. This vector is weighted by $\textbf{R}$, i.e., when a pixel belongs only partially to a region $i$ according to $\textbf{R}$ it is counted only partially for $\textbf{Y}_\textbf{B}^i$. From $\textbf{Y}_\textbf{B}^i$ the size $s_\textbf{B}^i$ of a region in pixels and the region's class probability distribution $\textbf{P}_\textbf{B}^i$ can be easily computed:
\begin{equation}
    s_\textbf{B}^i = \sum_{c=1}^{|C|} \textbf{Y}_\textbf{B}^i(c)
\end{equation}
\begin{equation}
    \textbf{P}_\textbf{B}^i = \frac{\textbf{Y}_\textbf{B}^i}{s_\textbf{B}^i}
\end{equation}
where $\textbf{Y}_\textbf{B}^i(c)$ is the number of pixels in region $i$ of block $\textbf{B}$ belonging to class $c$.

Our first new loss specifically punishes region maps $\textbf{R}$ which define segmentations with regions $i$ in blocks $\textbf{B}$ which contain more than one class according to the ground truth. We calculate
\begin{equation}
    \mathcal{L}_\textbf{Y} = \frac{1}{N} \sum_{i,\textbf{B}} H(\textbf{P}_\textbf{B}^i)
    \label{eq:lossY}
\end{equation}
with the Gini impurity $H(\textbf{P}_\textbf{B}^i) = 1 - \sum_{c=1}^{|C|} \left(\textbf{P}_\textbf{B}^i(c)\right)^2$ and $N$ equal to the number of blocks $\textbf{B}$ multiplied by the number of leaf nodes $i$ per BSP tree. $\mathcal{L}_\textbf{Y}$ is the average Gini impurity across all regions $i$ and blocks $\textbf{B}$. We use $H$ as a proxy for the entropy. Using the actual entropy instead did not produce better results, however, due to including a logarithm, it was slower and less numerically stable to calculate. $\mathcal{L}_\textbf{Y}$ ensures that every predicted region $i$ contains only a single class according to the ground truth as this loss becomes smaller the closer $\textbf{P}_\textbf{B}^i$ is to a one-hot vector, i.e., only one class is in the partition $i$ of block $\textbf{B}$.

We also compute
\begin{equation}
    \mathcal{L}_s = \frac{1}{N} \sum_{i,\textbf{B}} \max\{s_\text{min} - s_\textbf{B}^i, 0\}
    \label{eq:minregionsizeloss}
\end{equation}
with a hyperparameter $s_\text{min}$ specifying a minimum desired region size. A minimum region size ensures that the predicted lines used for partitioning intersect the blocks they belong to. This is done to make sure that the model quickly learns to utilize all available regions rather then trying to find a solution that utilizes only a subset of the regions which may result in a larger error long-term but may present itself as an undesirable local minimum during training. We observed that most $8\times 8$ blocks contain only a single class according to the ground truth. Therefore, using only a single or very few partitions/regions per block may present itself as a local minimum to the model that is hard to avoid or to get out of during training.

While we do not enforce constrains on the normal $\textbf{n}$ and the distance $d$ used to calculate the signed distance function $f$ in Eq. \ref{eq:nosigmoid}, we still want our model to favor predictions which result in sharp boundaries between regions $i$, rather than having lots of pixels which belong to multiple regions partially. The hyperparameter $\lambda$ can address this issue in theory, however, for any adjustments made to $\lambda$, the model may learn to adjust its predictions of $\textbf{n}$ and $d$ accordingly, rendering $\lambda$ useless. We therefore also compute the loss
\begin{equation}
    \mathcal{L}_\textbf{R} = \frac{1}{|\textbf{I}|} \sum_{(x,y) \in \textbf{I}} H(\textbf{R}_{yx})
\end{equation}
where $\textbf{I}$ is the set of image pixels and $\textbf{R}_{yx}$ is the probability vector defining how likely pixel $(x, y) \in \textbf{I}$ belongs to any of the regions $i$. By minimizing $H$, we make sure that our model favors predictions of $\textbf{R}_{yx}$ which are one-hot, thus defining sharp boundaries between regions. Again, the Gini impurity $H$ serves as a proxy to the entropy producing results of equal quality.

The full loss function we use for model training consists of cross-entropy $\mathcal{L}_\text{CE}$ and the three additional region map $\textbf{R}$ specific components just described:
\begin{equation}
    \mathcal{L}_\text{total} = \mu_1 \mathcal{L}_\text{CE} + \mu_2 \mathcal{L}_\textbf{Y} + \mu_3 \mathcal{L}_s + \mu_4 \mathcal{L}_\textbf{R}.
    \label{eq:totalloss}
\end{equation}
We set the loss weights $\mu_i$ s.t. $\sum_{i=1}^4 \mu_i = 1$. This constraint limits the range for hyperparameter optimization as any desired effective $\mu_i$ can be achieved by adjusting the constrained $\mu_i$ and the learning rate accordingly. $\mathcal{L}_\text{CE}$ is the most important component of the loss as it is the only component affecting the leaf nodes and therefore the predicted class logits. All four loss components affect the region map $\textbf{R}$ and therefore the inner nodes and the line parameters they contain. As initially described, $\mathcal{L}_\textbf{Y}$, $\mathcal{L}_s$ and $\mathcal{L}_\textbf{R}$ are specifically designed to improve the predicted $\textbf{R}$.

\subsection{Class-Specific Trees per Block}
In order to allow our model, SegForestNet to learn class-specific partitionings of $8 \times 8$ blocks into regions, we allow partitioning the classes $C$ into subsets $C_j$ and predicting a separate BSP tree for each subset. By creating a single subset with all classes our model can emulate the original BSPSegNet. By creating subsets which contain exactly one class we can create class-specific trees. Our approach generalizes to any partitioning of $C$ though.

For each subset $C_j$ we create a separate pair of decoders (shape and content). The encoder's output is split among all subsets s.t. there is no overlap in features going into any pair of two distinct decoders. Also, the feature dimensions and the split are specified s.t. $F_{encoder} < F_{tree}$ for all decoders (see Fig. \ref{fig:decoder}). The output dimension of the content decoder for each $C_j$ is set to $|C_j|\cdot N_\text{leaves}$ where $N_\text{leaves}$ is the number of leaf nodes per tree, i.e., each content decoder only predicts class logits for the classes in its respective $C_j$. For each $C_j$ we obtain a semantic segmentation with a reduced number of classes. To obtain the final semantic segmentation with all classes we simply concatenate the class logits predicted by all the trees.

The computation of $\mathcal{L}_\text{total}$ so far only considered a single predicted tree. However, all loss components but $\mathcal{L}_\text{CE}$ need to account for the fact that there are now multiple predicted trees. First, we calculate a separate $\textbf{Y}_\textbf{B}^i$ (Eq. \ref{eq:yblockregion}) for each $C_j$. To do so, we split the one-hot vector $\textbf{Y}_{yx}$ into $\textbf{Y}_{yx}^j$, which consists of all the classes in $C_j$, and $\textbf{Y}_{yx}^{\lnot j}$, which consists of all the other classes. From these two vectors we compute a new one-hot class vector
$\textbf{Y}_{yx}^{j'}$ by appending $\sum_{c \in C_{\lnot j}} \textbf{Y}_{yx}^{\lnot j}(c)$ to $\textbf{Y}_{yx}^j$. The one-hot vector $\textbf{Y}_{yx}^{j'}$ then consists of all classes in $C_j$ with an extra dimension representing all other classes. This vector can be used to compute tree-specific versions of $\textbf{Y}_\textbf{B}^i$ which can then be used for tree-specific loss components $\mathcal{L}_\textbf{Y}^j$, $\mathcal{L}_s^j$ and $\mathcal{L}_\textbf{R}^j$. Finally, we use the equations
\begin{equation}
    \mathcal{L}_\textbf{Y} = \sum_j \frac{|C_j|}{|C|} \cdot \mathcal{L}_\textbf{Y}^j
\end{equation}
\begin{equation}
    \mathcal{L}_s = \sum_j \frac{|C_j|}{|C|} \cdot \mathcal{L}_s^j
\end{equation}
\begin{equation}
    \mathcal{L}_\textbf{R} = \sum_j \frac{|C_j|}{|C|} \cdot \mathcal{L}_\textbf{R}^j
\end{equation}
to calculate the loss components necessary for Eq. \ref{eq:totalloss}, i.e., we weight each tree-specific loss by the number of classes in its respective subset $C_j$.

\subsection{Generalized Differentiable Rendering of Trees}
\label{sec:generalization}
The approach of using a region map $\textbf{R}$ for differentiable rendering generalizes to other signed distance functions $f$ and tree structures. In this last subsection we discuss alternatives to $f$ and how to implement other tree structures than binary space partitioning trees.

Eq. \ref{eq:nosigmoid} can be used with any signed distance function $f$, i.e., with any function which satisfies the following conditions:
\begin{enumerate}
    \item $f(\textbf{p}) = 0 \Leftrightarrow$ the point $\textbf{p}$ lies exactly on a boundary defined by $f$
    \item $\text{sgn}(f(\textbf{p}))$ specifies on which side of the boundary defined by $f$ the point $\textbf{p}$ lies
    \item $|f(\textbf{p})|$ is the distance of $\textbf{p}$ to the boundary defined by $f$
\end{enumerate}
The third condition can even be relaxed for the use in Eq. \ref{eq:nosigmoid}. It is sufficient if $|f(\textbf{p})|$ increases monotonically with the distance of $\textbf{p}$ to the boundary instead of it being the actual distance. Especially if $\textbf{p}$ is far from the boundary an approximate distance suffices.

\begin{figure*}
    \centering
    \includegraphics[width=.8\textwidth]{figures/sdfs.pdf}
    \caption{The decision boundaries created by signed distance functions based on different geometric primitives. From left to right: line, square, circle, ellipse, hyperbola, parabola. The blue area shows points for which the respective signed distance function is non-negative. As an example, in the inside of the circle $f_3$ is negative whereas it is positive on the outside.}
    \label{fig:sdfs}
\end{figure*}

So far, we only used lines as partition boundaries. Therefore, three values need to be predicted for each inner node: two values for the normal vector $\textbf{n}$ and one value for the distance $d$ to the origin. The signed distance $f_1(\textbf{p})$ of a point $\textbf{p}$ to the line can then be computed as
\begin{equation}
    f_1(\textbf{p}) = \textbf{n}\cdot\textbf{p} - d.
    \label{eq:linesdf}
\end{equation}
Fig. \ref{fig:sdfs} shows partitionings defined by signed distance functions based on other geometric primitives. The corresponding equations are as follows. To compute the approximate signed distance $f_2(\textbf{p})$ to a square, its center $\textbf{x}$ (two values) and its size $s$ (one value) need to be predicted:
\begin{equation}
    f_2(\textbf{p}) = \max\{|\textbf{x}_1 - \textbf{p}_1|,|\textbf{x}_2 - \textbf{p}_2|\} - s.
\end{equation}
In this equation $\textbf{x}_i$ and $\textbf{p}_i$ are the two components of the vectors $\textbf{x}$ and $\textbf{p}$. As with the square, the signed distance $f_3(\textbf{p})$ to a circle requires a predicted center $\textbf{x}$ and radius $r$ (one value) as well:
\begin{equation}
    f_3(\textbf{p}) = \lVert \textbf{x} - \textbf{p} \rVert_2^2 - r.
    \label{eq:circlesdf}
\end{equation}
This formulation is a pseudo signed distance which avoids the need to compute a square root and expects the model to learn to predict the square of the radius. The boundary of an ellipse can be defined by $d_1 + d_2 = c$ where $d_i$ are distances to fixed points and $c \in \mathbb{R}^+$ is a constant. To compute an approximate signed distance $f_4(\textbf{p})$, two points $\textbf{x}$ and $\textbf{y}$, each consisting of two values, need to be predicted as well as one value for the constant $c$. The signed distance then is
\begin{equation}
    f_4(\textbf{p}) = \lVert \textbf{x} - \textbf{p} \rVert_2 + \lVert \textbf{y} - \textbf{p} \rVert_2 - c.
\end{equation}
Similar to an ellipse, a hyperbola can be defined by $d_1 - d_2 = c$ with $d_i$ and $c$ as before. The corresponding approximate signed distance $f_5(\textbf{p})$ is
\begin{equation}
    f_5(\textbf{p}) = \left|\lVert \textbf{x} - \textbf{p} \rVert_2 - \lVert \textbf{y} - \textbf{p} \rVert_2\right| - c.
\end{equation}
A parabola can be defined as the set points satisfying $d_1 = d_2$ where $d_1$ is the distance to a fixed point and $d_2$ is the distance to a line. Therefore, five values need to be predicted, two for a fixed point $\textbf{x}$, two for a normal vector $\textbf{n}$ and one for a distance $d$ of the line to the origin. The approximate signed distance $f_6(\textbf{p})$ is
\begin{equation}
    f_6(\textbf{p}) = \lVert \textbf{x} - \textbf{p} \rVert_2 - f_1(\textbf{p})
\end{equation}
All (pseudo) signed distance functions $f_i$ except $f_2$ use formulations s.t. that arbitrarily rotated, scaled and translated instances of the underlying geometric primitives used for partitioning can be predicted.

\begin{figure}
    \centering
    \includegraphics[width=.6\columnwidth]{figures/kdtree.pdf}
    \caption{Visualization of a $k$-d tree (left) and a region it partitions (right). The parameters of each inner node (blue) are a fixed dimension, indicated by the orientation of the line used for partitioning, and a predicted threshold, indicated by the position of the line along the fixed dimension.}
    \label{fig:kdtree}
\end{figure}

While BSPSegNet and SegForestNet only use BSP trees, differentiable rendering of other tree structures is possible by using the same region map $\textbf{R}$ approach. A $k$-d tree \cite{bentley1975multidimensional}, as shown in Fig. \ref{fig:kdtree}, is spatial partitioning data structure similar to a BSP tree that be used to partition a space with an arbitrary number of dimensions. It is also a binary tree. Instead of using a signed distance function to decide on which side of a boundary a given point lies, each inner node of $k$-d tree uses a threshold and a dimension index to decide whether a point lies in the left or right subset. The threshold $t$ is a value that needs to predicted while the dimension index $i$ is a fixed value chosen when instantiating the model. Eqs. \ref{eq:nosigmoid}-\ref{eq:rright2} can be used for $k$-d trees as well with the signed distance function
\begin{equation}
    f_7(\textbf{p}) = t - \textbf{p}_i
\end{equation}
where $\textbf{p}_i$ is the $i$-th component of the vector $\textbf{p}$. An inner node of a $k$-d tree requires only a single parameter to be predicted as opposed to the three to five parameters required for the signed distance functions used by an inner node of a BSP tree. However, generally $k$-d trees need to be deeper, e.g., to be able to define slopes, which counteracts this advantage. This depth disadvantage can be mitigated by making $i$ a predicted parameter. We call a tree with such inner nodes a dynamic $k$-d tree. A downside of this mitigation strategy is that one or two additional parameters (depending on the implementation) need to be predicted per inner node to decide whether the first dimension $i=1$ or the second dimension $i=2$ shall be used for partitioning. Even with this extra flexibility a dynamic $k$-d tree still has a depth disadvantage in the worst case, e.g., in the aforementioned slope example.

\begin{figure}
    \centering
    \includegraphics[width=.6\columnwidth]{figures/quadtree.pdf}
    \caption{Visualization of a quadtree (left) and a region it partitions (right). The parameters of each inner node (blue) are the predicted coordinates of a point. This point partitions a given region into four quadrants which are then children of the associated inner node.}
    \label{fig:quadtree}
\end{figure}

A quadtree \cite{finkel1974quad}, shown in Fig. \ref{fig:quadtree}, is a spatial partitioning data structure specifically designed to partition two dimensional planes. In each inner node of a quadtree a point $\textbf{x}$ is used to partition the space into four quadrants. Compared to a $k$-d tree, a quadtree inner node partitions a space along both dimensions instead of only a single fixed dimension $i$ and it uses a different threshold $t$ for each dimension by using the components of $\textbf{x}$ as thresholds. The equations
\begin{equation}
    t_1 = \textbf{x}_1 - \textbf{p}_1, t_2 = \textbf{x}_2 - \textbf{p}_2
\end{equation}
\begin{equation}
    t_3 = \text{ReLU}(t_1), t_4 = \text{ReLU}(-t_1)
\end{equation}
\begin{equation}
    t_5 = \text{ReLU}(t_2), t_6 = \text{ReLU}(-t_2)
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} + \lambda \cdot t_4 \cdot t_5
    \label{eq:quadtreetl}
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} + \lambda \cdot t_3 \cdot t_5
    \label{eq:quadtreetr}
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} + \lambda \cdot t_4 \cdot t_6
    \label{eq:quadtreebl}
\end{equation}
\begin{equation}
    \textbf{R}_{iyx} := \textbf{R}_{iyx} + \lambda \cdot t_3 \cdot t_6
    \label{eq:quadtreebr}
\end{equation}
are used to update the region map $\textbf{R}$ for a quadtree at a location $\textbf{p} = (x,y)$. Note: in Eqs. \ref{eq:quadtreetl}-\ref{eq:quadtreebr} $i$ refers to indices of leaf nodes reachable from one of the child nodes. In Eq. \ref{eq:quadtreetl} $i$ needs to be set to the indices of the leaf nodes covering the top-left partition, in Eq. \ref{eq:quadtreetr} $i$ refers to the top-right partition, in Eq. \ref{eq:quadtreebl} $i$ refers to the bottom-left partition and in Eq. \ref{eq:quadtreebr} $i$ refers to the bottom-right partition.

This region map $\textbf{R}$ based framework is general enough to even allow for trees in which different inner nodes are instances of different tree types (BSP, $k$-d, quadtree) or use different signed distance functions in the case of BSP tree inner nodes. When using class-specific trees, i.e., when partitioning the set of classes $C$ into subsets $C_j$, separate region maps specific to each of the subsets have to be computed anyway. Consequently, different tree configurations can be used for each of the subsets $C_j$ if so desired.

\section{Evaluation}
In this section we evaluate our contributions. First, we describe the datasets and the test setup we used. We will then present a toy example which proves that our generalizations in section \ref{sec:generalization} may be helpful if the alternate signed distance functions or tree structures suit the data. We then explore good value ranges for the hyperparameters we introduced, followed by an ablation study. We finish the section with a comparison to state-of-the-art models.

\subsection{Datasets}
We used seven different datasets for evaluation. We used aerial images, in particular true orthophotos, from the German cities Hannover, Buxtehude, Nienburg, Vaihingen and Potsdam. The later two datasets are part of the ISPRS 2D Semantic Labeling Benchmark Challenge \cite{ISPRSBenchmark}. Additionally, we used images from Toulouse \cite{roscher2020semcity}. The Toulouse dataset also includes instance labels for buildings for performing panoptic segmentation or instance segmentation. However, we only used the semantic labels. Our last dataset was iSAID \cite{waqas2019isaid}, which is actually an instance segmentation dataset based on the object detection dataset DOTA \cite{Xia2018CVPR}. Again, we ignored the instance labels. However due to its origin, iSAID has very sparse non-background labels, i.e., it cannot be considered a true semantic segmentation datasets which have much denser non-background label coverage.

Hannover, Buxtehude and Nienburg consist of 16 patches of $2500\times 2500$ pixels each. We omitted one such patch from Hannover since it consisted almost entirely of trees. We randomly divided the image patches into subsets s.t. roughly $70\%$ of pixels were used for training with the remainder split evenly into validation and test sets. Since Toulouse consists of only four images, we used two for training, one for validation and one for testing. iSAID is already pre-split into training and validation. We used the validation images as our test images and split the training set into training ($95\%$ of pixels) and actual validation. For Vaihingen and Potsdam, we used the images originally released without ground truth as test images and used the remaining images as training (roughly $90\%$ of pixels) and validation (roughly $10\%$ of pixels).

The ground truth classes of the German city images are impervious surfaces, buildings, low vegetation, tree, car, and clutter/background, the later two being rare. The Toulouse dataset also contains sports venues and water as separate classes in addition to the classes used for segmentation of the German cities. In Toulouse the classes car, sports venues, water, and background are rare. Since the background class is so rare and also is used to express uncertainty about the true class in the case of at least Potsdam, we ignored this class for all metrics, including training losses. The exception is the iSAID. Due to its sparsity of non-background pixels, less than $3\%$ of pixels belong to a non-background class, we did not ignore the background class in this case. iSAID has 15 non-background classes including different vehicle and sports venue categories.

For data augmentation, we randomly sampled 8000 input patches from the training images. The position (random translation) of the input patches was chosen s.t. input patches from the entire image could be sampled. The size of the square area being sampled was randomly scaled between $70\%$ and $130\%$ of $224\times 224$ pixels, the actual input size, making objects randomly appear smaller or larger for augmentation. The same scaling was applied to both spatial dimensions. The random shearing was chosen independently for both axes, sampled uniformly from $[-6\degree, 6\degree]$. Random rotation depended on the dataset and was chosen from $[-50\degree, 50\degree]$ for Vaihingen and from $[-30\degree, 30\degree]$ for all other datasets. This dataset-specific augmentation for Vaihingen improved results on this dataset slightly. No flipping (horizontal or vertical) was used. All channels were normalized to zero mean and a standard deviation of $1$. After normalization, noise sampled from a normal distribution with standard deviation $0.1$ was added. The noise was sampled for each pixel and each channel individually. Furthermore, we used bilinear filtering. For the ground truth, we used the bilinear interpolation coefficients as voting weights to determine the discrete class for each pixel. Validation and test images were partitioned into axis-aligned input patches with no overlap and stride equal to the input patch size. We performed this input patch creation once for each dataset. Additionally, we computed the NDVI \cite{rouse1973monitoring} for all datasets but iSAID, which lacks the necessary infrared channel. Other exceptions made for iSAID were sampling 28000 input patches and ensuring that the entropy of the class distribution of each sample was at least $0.04$ to ensure that the patches actually contained non-background pixels.

We used all available input channels and therefore adapted the number of input channels of the very first convolution of every model to the dataset used. All datasets except Vaihingen, which has no green channel, consist of at least RGB images with a few grayscale images in iSAID which we converted to RGB. The German cities all also have infrared channels and a digital surface model ("depth"). Toulouse actually consists of multispectral images with eight channels, including RGB and infrared.

\subsection{Training and Test Setup}
We used PyTorch v1.10 with CUDA 10.2 to train models on nVidia GeForce RTX 2080 Ti GPUs. Depending on the backbone used, we set the mini-batch size to 36 (MobileNetv2 \cite{sandler2018mobilenetv2}) or 18 (all other backbones). As an optimizer we used AdamW
\cite{kingma2014adam,loshchilov2017decoupled} with a cosine annealing learning rate schedule \cite{loshchilov2016sgdr}. As loss function we either used categorical cross entropy or our novel loss function (Eq. \ref{eq:totalloss}). For the cross entropy part, we set class weights s.t. rarer classes have higher weights based on the class distribution in the training images.

We used random search for hyperparameter optimization for all models, in particular optimizing their learning rates. For this, we
computed the mean Intersection-over-Union (mIoU) on the validation subset of Buxtehude. We then used the maximum validation mIoU across all epochs of a run with a randomly chosen value for the hyperparameter being optimized to evaluate the performance of this run and thus hyperparameter value. The final hyperparameters used can be found in our published code\footnote{\texttt{core/defaults.yaml} and \texttt{cfgs/semseg.yaml}}. After hyperparameter optimization, we ran each experiment ten times and computed means and standard deviations from these ten runs. For the final comparison of our contributions to the state-of-the-art, we again used the maximum validation mIoU to identify the epoch after which each model performed best and then evaluated these particular model parameters on the test data.

We also made some small modifications to state-of-the-art models and backbones to either improve their performance slightly or simplify them while keeping their performance the same: we use LeakyReLUs \cite{xu2015empirical} instead of ReLUs and we removed all atrous/dilated convolutions. For DeepLabv3+ we modified the strides of the backbone s.t. the spatial bottleneck had a downsampling factor of $8$. Training all models from scratch, we changed the first convolution in each backbone to accept all input channels offered by each respective dataset.

For the decoder of SegForestNet we used eight residual blocks with $F_{intermediate} = 96$ (see Fig. \ref{fig:decoder}). For the shape features, $F_{encoder}$ and $F_{tree}$ were set to $8$ and $9$ respectively. For the content features, we chose $F_{encoder}=16$ and $F_{tree}=4\cdot |C|$ with $|C|$ being the number of classes in the dataset. The numbers for $F_{tree}$ are the result of using BSP trees with depth $=2$, i.e., trees with three inner nodes and four leaf nodes.

\subsection{Toy Experiment}
While the generalization to different signed distance functions and different tree structures does not provide benefits for the semantic segmentation of aerial images, we still performed two experiments using toy datasets to show that these generalization may prove useful in other contexts. In street scenes captured by cameras attached to cars or in industrial applications rounded shapes, e.g., tires, pipes or drums, may occur frequently and thus a signed distance function capturing such shapes more accurately, e.g. $f_3$ (Eq. \ref{eq:circlesdf}), may prove beneficial.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/toy0b.pdf}
    \caption{Toy dataset showing randomly created circles on a solid background. The ground truth visualization is equal to the input images, as there was a one-to-one correspondence of colors to classes. Using a signed distance function more suitable to the data (third row) produces more accurate predictions compared to using the default line-based signed distance function (second row). The restricted models in rows two and three were chosen to make the difference more visually obvious. The last row shows the prediction of a less restricted model to prove that the alternate signed distance functions can produce good results.}
    \label{fig:toysdfsample}
\end{figure}

We created the first toy dataset by generating $128\times 128$ pixel RGB images showing circles as depicted in Fig. \ref{fig:toysdfsample}. We only used the eight corners of the RGB cube as colors and assigned a different class to each corner/color. One color was used as background and circles of random translation, size and color were drawn onto the background. The background color was never used for any of the random circles. We used a model which subdivides the image into blocks of size $32\times 32$ instead of $8\times 8$. For each block a shallow BSP tree of depth one was predicted. We heavily limited the prediction capability of the model to make the visual difference more clear. We trained the model ten times using $f_1$ (Eq. \ref{eq:linesdf}) in the inner nodes and ten times using $f_3$ (Eq. \ref{eq:circlesdf}) instead. Using $f_1$, the model achieves an average validation mIoU of $63.3\%$ ($\pm 0.9\%$) across the ten trainings, while the $f_3$-based model achieves $64.9\%$ ($\pm 1.9\%$). The model with the more data-appropriate signed distance function, $f_3$, performed better.

While a tree depth of one, i.e., only a single inner node in each tree, might suggest that every $32\times 32$ block may contain only two different classes, there are blocks with three or more classes. An example is visible in the bottom-left of the $f_1$-based prediction in the last column of Fig. \ref{fig:toysdfsample}. This is due to the linear combination of the class logit predictions in the leaf nodes. To understand this phenomenon, suppose the predicted class logits are represented by the vectors $(v_1, v_2, 0)$ and $(0, v_2, v_3)$ with $v_1>v_2>0$ and $v_3>v_2>0$ with the five missing dimensions for the other classes all being $0$ in both vectors. For any pixel in the segmentation the region map consists of a vector of the form $(\alpha, 1-\alpha)$, i.e., the predicted class logits are $\alpha\cdot(v_1,v_2,0) + (1-\alpha)\cdot(0,v_2,v_3) = (\alpha v_1, v_2, (1-\alpha)v_3)$. Depending on the value of $\alpha \in (0,1)$, which itself is dependent on the signed distance of the pixel to the decision boundary, the largest class logit may be $\alpha v_1$, $v_2$, or $(1-\alpha)v_3$, thus resulting in three regions with different classes despite there being only a single decision boundary. With other class logit distributions even more distinct regions are possible.

\begin{figure}
    \centering
    \includegraphics[width=.6\columnwidth]{figures/toy12sample.pdf}
    \caption{Toy dataset for investigating the impact of different tree structures. Again, the ground truth visualization is equal to the input images. Both tested tree structures can segment the images in the dataset perfectly eventually if the model is trained sufficiently long.}
    \label{fig:toytreesample}
\end{figure}

The toy dataset for the second experiment was created by partitioning $128\times 128$ pixel RGB vertically into two segments and then partitioning each segment horizontally resulting in examples as shown in Fig. \ref{fig:toytreesample}. As before, we used the corners of the RGB cube as colors and classes. The resulting images can be perfectly segmented by a BSP tree of depth two using $f_1$ as signed distance function as well as by a $k$-d tree of depth two. However, the $k$-d tree needs fewer parameters to describe the inner nodes, thus is more efficient. For both tree structures, BSP trees and $k$-d trees, we trained ten models using this toy dataset. Both models eventually learned to segment the toy images perfectly, with the $k$-d trees having the negligible advantage of reaching this state an epoch or two sooner. The experiment still proves that region map-based differentiable rendering generalizes to other tree structures which may be desirable for types of data other than aerial images.

\subsection{New Hyperparameters}
We introduced several new hyperparameters throughout our contributions, $\lambda$ in Eq. \ref{eq:nosigmoid} and $\mu_1$ to $\mu_4$ in Eq. \ref{eq:totalloss}. To study them, we used SegForestNet with MobileNetv2 as backbone, BSP trees using $f_1$ (Eq. \ref{eq:linesdf}) as signed distance function and a tree depth of 2. As we expect $\lambda$ to have no real effect anyway, especially since we introduced $\mu_4$, we simply set $\lambda = 1$ and did not study it further. We ran random searches with 100 to 200 iterations to optimize each of the different $\mu_i$. We optimized $\mu_3$ and $\mu_4$ together in the same random search. We set $s_\text{min} = 8$ (Eq. \ref{eq:minregionsizeloss}), i.e., each of the four regions should use at least one eighth of the area of each block. We chose this as a compromise, so that each region in a $8\times 8$ block contributes a significant part without restricting the model too much by setting $s_\text{min}$ too high.

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{figures/hpo.pdf}
    \caption{Validation mIoU on Buxtehude across $\mu_i$ (Eq. \ref{eq:totalloss}). As all $\mu_i$ except $\mu_3$ and $\mu_4$ were optimized in separate random searches, the corresponding mIoU ranges differ. However, as only finding the proper settings for $\mu_i$ is important, the exact mIoU ranges do not matter anyway. The highlighted intervals show the optimal range for each $\mu_i$.}
    \label{fig:hyperparameters}
\end{figure*}

The results of our random searches are shown in Fig. \ref{fig:hyperparameters}. The impact of cross-entropy, which was controlled by $\mu_1$, was the most significant. It also had a rather large range of $0.91$ to $0.95$ in which it produced optimal results. The other $\mu_i$ produced a slight increase in performance as they were increased, especially when considering the lower bound of the performance, with rather sharp drops in performance beyond a certain threshold. For $\mu_2$ that threshold was $0.035$, while it was $0.01$ for $\mu_3$ and $\mu_4$. We therefore chose values for $\mu_2$ to $\mu_4$ which were just below that threshold: $\mu_2=0.034$ and  $\mu_3=\mu_4=0.0095$. We then computed $\mu_1 = 1 - (\mu_2+\mu_3+\mu_4) = 0.947$ which is still in the optimal range for $\mu_1$. We therefore picked these values for the hyperparameters in Eq. \ref{eq:totalloss}.

\begin{figure*}
    \centering
    \includegraphics[width=.8\textwidth]{figures/mu4.pdf}
    \caption{Visualization of the region map $\textbf{R}$ and the Gini impurity $H$ (Eq. \ref{eq:lossY}) for different loss functions and different values for $\mu_i$ (Eq. \ref{eq:totalloss}). The second to last column used optimized values for $\mu_i$, while the last column used $\mu_1 = 0.95$, $\mu_2=\mu_3=0$, and $\mu_4=0.05$. These visualizations show that our novel loss achieves the desired effects of sharp region boundaries and minimum region sizes (middle row) and pure regions containing just a single class according to the ground truth (bottom row).}
    \label{fig:regionmapvis}
\end{figure*}

In order to verify that the different loss components have the desired effects, we created a visualization of the region map $\textbf{R}$ and the Gini impurity $H$ (Eq. \ref{eq:lossY}). To create these visualizations, we used a variant of Eq. \ref{eq:finalprediction} in which we used colors or $H(\textbf{P}_\textbf{B}^i)$ in place of the class logits $\textbf{v}_i$. Using colors (red, green, blue, and cyan) instead of class logits created a visualization showing which pixel belongs to which of the four regions in a block (Fig. \ref{fig:regionmapvis}, middle row). Using $H(\textbf{P}_\textbf{B}^i)$ instead created a visualization highlighting pixels in gray or even white which belong to regions which contain more than one class according to the ground truth (Fig. \ref{fig:regionmapvis}, bottom row). Using our novel loss with optimized values for $\mu_i$ created region maps with the desired minimum region size and sharp region boundaries. Furthermore, the average $H(\textbf{P}_\textbf{B}^i)$ are smaller compared to using cross-entropy. Cross-entropy also created blurry boundaries and seemed to associate certain regions with certain classes: the green region was used for buildings (blue in the ground truth), the cyan region was used for trees (red in the ground truth) and the red region for the remaining classes. Using our novel loss but with unoptimized $\mu_i$, in particular setting $\mu_2$ (pure regions containing only class according to the ground truth) and $\mu_3$ (minimum region size) to $0$, showed that the individual losses have the desired effects. The average $H(\textbf{P}_\textbf{B}^i)$ is slightly higher (hard to see in Fig. \ref{fig:regionmapvis}) and certain regions (blue and cyan) are never used. Still, the region boundaries are sharp since $\mu_4$ is sufficiently high.

\subsection{Ablation Study}
In an ablation study, shown in Table \ref{tab:ablation}, we examined the effect of each contribution individually. As before, we used SegForestNet with MobileNetv2 as backbone. The new region map computation (Eqs. \ref{eq:nosigmoid} to \ref{eq:rright2}) had the greatest impact. Not only did it improve performance and reduce the variance, it was also required for the novel loss to actually provide a benefit. Without the new region map computation the new decoder architecture also provided a significant performance improvement but without the reduction in variance. With the new region map computation there was still a small performance improvement provided by the new decoder architecture. Furthermore, the new decoder architecture reduced the total model size by $5\%$ due to each decoder using only $85.3$k parameters instead of $137.7$k. As mentioned before, the novel loss only was beneficial in combination with the new region map computation and even then the improvements were rather small. Overall, there was an improvement of almost $1\%$ when using all contributions. We did not include class-specific trees in this ablation study as their benefit is dataset dependant which we show in the next section.

\begin{table}
\caption{Ablation study. The top-most row is SegForestNet with all our contributions, while the bottom-most row is the original BSPSegNet. The last column shows max. validation mIoU on Buxtehude.}
\label{tab:ablation}
\begin{center}
\begin{tabular}{ccc|c}
new & new region & novel & \\
decoder & map computation & loss & mIoU \\
\hline
\checkmark & \checkmark & \checkmark & $77.4\% \pm 0.2\%$ \\
\checkmark & \checkmark & & $77.3\% \pm 0.1\%$ \\
\checkmark & & \checkmark & $76.8\% \pm 1.5\%$ \\
\checkmark & & & $77.3\% \pm 0.3\%$ \\
\hline
 & \checkmark & \checkmark & $77.3\% \pm 0.2\%$ \\
 & \checkmark & & $77.1\% \pm 0.2\%$ \\
 & & \checkmark & $73.1\% \pm 1.9\%$ \\
 & & & $76.5\% \pm 1.5\%$ \\
\end{tabular}
\end{center}
\end{table}

\subsection{Comparison to State-Of-The-Art}
Finally, we compared SegForestNet to state-of-the-art models. Results are shown in Table \ref{tab:miou} and Fig. \ref{fig:miou} with samples of predicted segmentations in Fig. \ref{fig:samples}. We grouped models according to their number of parameters (Table \ref{tab:modelsize}) in order to be able to compare models of similar size and complexity with each other. We created four groups, small models (less than five million parameters), medium size models (around ten million parameters), large models (20 to 30 million parameters) and very large models (more than 30 million parameters). Since all tested models except U-Net and U-Net++ have interchangeable backbones we tested all models in a smaller variant, using MobileNetv2 as a backbone, and a larger variant, using Xception as a backbone. With these two backbones, SegForestNet appears in the small model and the large model group. Additionally, SegForestNet was tested in two variants, the base variant in which only a single BSP tree per $8\times 8$ block was predicted, and a variant, SegForestNet*, in which one BSP tree per class per block was predicted. Some older models, e.g., U-Net, still show competitive performance and thus we chose a selection of older and newer models to compare SegForestNet to: FCN \cite{long2015fully}, DeepLab v3+ \cite{chen2018encoder}, (S-)RA-FCN \cite{mou2019relation}, PFNet \cite{li2021pointflow}, FarSeg \cite{zheng2020foreground}, U-Net++ \cite{zhou2018unetpp} and U-Net \cite{ronneberger2015u}.

\begin{table*}
\caption{Average test mIoU across ten runs for various models. Our models are SegForestNet and SegForestNet*. Models marked with $\dag$ use Xception as the backbone, while all others use MobileNetv2 instead. The two U-Net variants are an exception, since they have a fixed architecture rather then exchangeable backbones. The models are grouped according to their size (Table \ref{tab:modelsize}) with the best model in each bracket highlighted in bold.}
\label{tab:miou}
\begin{center}
\begin{tabular}{r|ccccccc}
 & Hannover & Nienburg & Buxtehude & Potsdam & Vaihingen & Toulouse & iSAID \\
\hline
$\text{FCN}$ & $72.6\%\pm 0.2\%$ & $73.2\%\pm 0.4\%$ & $76.2\%\pm 0.2\%$ & $78.1\%\pm 0.2\%$ & $72.0\%\pm 0.2\%$ & $54.1\%\pm 2.9\%$ & $38.3\%\pm 0.9\%$ \\
$\text{SegForestNet}$ & $72.9\%\pm 0.5\%$ & $73.8\%\pm 0.3\%$ & $76.2\%\pm 0.3\%$ & $\textbf{78.9}\%\pm 0.3\%$ & $72.9\%\pm 0.3\%$ & $52.9\%\pm 3.2\%$ & $44.5\%\pm 0.7\%$ \\
$\text{DeepLab v3+}$ & $72.9\%\pm 0.3\%$ & $73.8\%\pm 0.3\%$ & $76.5\%\pm 0.3\%$ & $77.8\%\pm 0.8\%$ & $72.2\%\pm 0.3\%$ & $52.8\%\pm 3.8\%$ & $34.9\%\pm 1.2\%$ \\
$\text{RA-FCN}$ & $71.0\%\pm 0.4\%$ & $71.0\%\pm 0.4\%$ & $74.1\%\pm 0.8\%$ & $74.1\%\pm 1.1\%$ & $70.3\%\pm 0.3\%$ & $49.4\%\pm 2.2\%$ & $35.4\%\pm 4.1\%$ \\
$\text{SegForestNet*}$ & $\textbf{73.6}\%\pm 0.2\%$ & $74.1\%\pm 0.3\%$ & $76.2\%\pm 0.2\%$ & $78.8\%\pm 0.3\%$ & $\textbf{72.9}\%\pm 0.2\%$ & $\textbf{54.2}\%\pm 1.6\%$ & $42.8\%\pm 0.5\%$ \\
$\text{PFNet}$ & $73.0\%\pm 0.4\%$ & $\textbf{74.2}\%\pm 0.3\%$ & $\textbf{76.8}\%\pm 0.4\%$ & $\textbf{78.9}\%\pm 0.3\%$ & $72.6\%\pm 0.2\%$ & $53.9\%\pm 1.3\%$ & $\textbf{45.8}\%\pm 0.4\%$ \\
\hline
$\text{FarSeg}$ & $72.9\%\pm 0.3\%$ & $74.2\%\pm 0.2\%$ & $76.7\%\pm 0.2\%$ & $\textbf{77.6}\%\pm 0.5\%$ & $\textbf{72.2}\%\pm 0.3\%$ & $54.5\%\pm 1.9\%$ & $43.4\%\pm 0.8\%$ \\
$\text{U-Net++}$ & $\textbf{73.8}\%\pm 0.3\%$ & $\textbf{75.6}\%\pm 0.2\%$ & $\textbf{77.4}\%\pm 0.2\%$ & $77.3\%\pm 0.3\%$ & $71.7\%\pm 0.2\%$ & $\textbf{56.8}\%\pm 2.0\%$ & $\textbf{45.4}\%\pm 0.6\%$ \\
\hline
$\text{FCN}^\dag$ & $73.4\%\pm 0.1\%$ & $74.1\%\pm 0.2\%$ & $76.4\%\pm 0.2\%$ & $79.4\%\pm 0.2\%$ & $72.7\%\pm 0.2\%$ & $54.3\%\pm 1.1\%$ & $48.4\%\pm 0.5\%$ \\
$\text{SegForestNet}^\dag$ & $73.7\%\pm 0.5\%$ & $74.5\%\pm 0.3\%$ & $76.8\%\pm 0.1\%$ & $79.1\%\pm 0.3\%$ & $73.1\%\pm 0.4\%$ & $55.1\%\pm 0.7\%$ & $50.4\%\pm 0.6\%$ \\
$\text{DeepLab v3+}^\dag$ & $73.5\%\pm 0.2\%$ & $74.7\%\pm 0.2\%$ & $76.6\%\pm 0.2\%$ & $78.8\%\pm 0.2\%$ & $72.7\%\pm 0.3\%$ & $55.2\%\pm 0.6\%$ & $48.9\%\pm 0.3\%$ \\
$\text{SegForestNet*}^\dag$ & $\textbf{74.0}\%\pm 0.3\%$ & $74.3\%\pm 0.5\%$ & $76.9\%\pm 0.2\%$ & $79.5\%\pm 0.2\%$ & $72.9\%\pm 0.3\%$ & $53.8\%\pm 1.3\%$ & $49.2\%\pm 0.8\%$ \\
$\text{PFNet}^\dag$ & $73.9\%\pm 0.2\%$ & $\textbf{74.9}\%\pm 0.2\%$ & $\textbf{77.2}\%\pm 0.2\%$ & $\textbf{79.6}\%\pm 0.4\%$ & $\textbf{73.3}\%\pm 0.3\%$ & $\textbf{55.2}\%\pm 1.6\%$ & $\textbf{50.6}\%\pm 0.2\%$ \\
$\text{FarSeg}^\dag$ & $54.3\%\pm 29.6\%$ & $61.7\%\pm 25.9\%$ & $56.3\%\pm 29.6\%$ & $64.3\%\pm 27.6\%$ & $59.5\%\pm 25.9\%$ & $54.8\%\pm 1.8\%$ & $46.3\%\pm 3.2\%$ \\
\hline
$\text{U-Net}$ & $\textbf{74.4}\%\pm 0.1\%$ & $\textbf{75.8}\%\pm 0.2\%$ & $\textbf{77.2}\%\pm 0.5\%$ & $\textbf{78.2}\%\pm 0.5\%$ & $\textbf{71.5}\%\pm 0.3\%$ & $\textbf{57.3}\%\pm 0.9\%$ & $\textbf{47.8}\%\pm 0.3\%$ \\
$\text{RA-FCN}^\dag$ & $67.0\%\pm 1.5\%$ & $69.1\%\pm 0.3\%$ & $71.4\%\pm 0.8\%$ & $72.6\%\pm 0.8\%$ & $65.5\%\pm 1.5\%$ & $48.2\%\pm 1.9\%$ & $34.0\%\pm 1.4\%$
\end{tabular}
\end{center}
\end{table*}

\begin{table}
\caption{Size of the tested models in terms of number of model parameters as a proxy for model complexity. We grouped the models in four brackets: small (less than 5M), medium (about 10M), large (20-30M), and very large (more than 30M). The symbol $\dag$ is used as in Table \ref{tab:miou}.}
\label{tab:modelsize}
\begin{center}
\begin{tabular}{r|c}
 & \# of model parameters [M] \\
\hline
$\text{FCN}$ & 1.8 \\
$\text{SegForestNet}$ & 2.0 \\
$\text{DeepLab v3+}$ & 2.3 \\
$\text{RA-FCN}$ & 2.3 \\
$\text{SegForestNet*}$ & 2.9 - 4.7 \\
$\text{PFNet}$ & 4.9 \\
\hline
$\text{FarSeg}$ & 9.1 \\
$\text{U-Net++}$ & 9.2 \\
\hline
$\text{FCN}^\dag$ & 20.9 \\
$\text{SegForestNet}^\dag$ & 21.1 \\
$\text{DeepLab v3+}^\dag$ & 22.2 \\
$\text{SegForestNet*}^\dag$ & 22.2 - 24.4 \\
$\text{PFNet}^\dag$ & 26.4 \\
$\text{FarSeg}^\dag$ & 29.2 \\
\hline
$\text{U-Net}$ & 31.4 \\
$\text{RA-FCN}^\dag$ & 40.0
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption{Count of how often each model appears in the Pareto front of a given dataset (compare Fig. \ref{fig:miou}). SegForestNet' is a combination of SegForestNet or SegForestNet*: if either of the later two is in the Pareto front of a given dataset, the count for SegForestNet' is increased by 1. The symbol $\dag$ is used as in Table \ref{tab:miou}.}
\label{tab:pareto}
\begin{center}
\begin{tabular}{c|c}
\# of times in &  \\
Pareto front & models \\
\hline
7 & $\text{FCN}$ \\
6 & $\text{SegForestNet'}$ \\
5 & $\text{SegForestNet}$ \\
4 & $\text{SegForestNet*}$, $\text{U-Net++}$, $\text{SegForestNet'}^\dag$ \\
3 & $\text{PFNet}$, $\text{PFNet}^\dag$, $\text{U-Net}$ \\
2 & $\text{FCN}^\dag$, $\text{SegForestNet}^\dag$, $\text{SegForestNet*}^\dag$ \\
1 & $\text{DeepLab v3+}$, $\text{FarSeg}$ \\
0 & $\text{RA-FCN}$, $\text{DeepLab v3+}^\dag$, $\text{FarSeg}^\dag$, $\text{RA-FCN}^\dag$
\end{tabular}
\end{center}
\end{table}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/miou.pdf}
    \caption{Visualization of the columns for Potsdam and Hannover from Table \ref{tab:miou}. Dots show means and bars indicate the standard deviation. The symbol $\dag$ is used as in Table \ref{tab:miou}.}
    \label{fig:miou}
\end{figure}

\begin{figure*}
	\begin{center}
		\includegraphics[width=\textwidth]{figures/samples.pdf}
		\includegraphics[width=.6\textwidth]{figures/sampleslegend.pdf}
	\end{center}
	\caption{Samples of validation and test examples from Hannover. All models except U-Net use Xception as backbone.}
	\label{fig:samples}
\end{figure*}

We trained all models ourselves: we used optimized parameters suggested by the respective authors, except that we optimized the learning rates via random search. All models were trained ten times. To avoid overfitting, we used the validation mIoU for each training run to identify the epoch in which the model performed best and then used the parameters from this best epoch to determine the mIoU on the test set. All results shown are the average test mIoUs.

There is no clear winner between SegForestNet and its variant SegForestNet*. Sometimes, e.g, when using the MobileNetv2 backbone to train a model on Potsdam or iSAID, SegForestNet performs better while at other times, e.g., same backbone but training on Hannover or Toulouse, SegForestNet* performs better. The situation may even differ for different backbones on the same dataset, e.g., when using the Xception backbone to train a model on Potsdam, SegForestNet* outperforms SegForestNet. Therefore, no clear recommendation whether to use class-specific segmentations (SegForestNet*) or not (SegForestNet) can be given. Rather, this is has to be determined based on the backbone and dataset to be used.

On all datasets except iSAID most models performed well, even older models. The exceptions were RA-FCN and FarSeg. RA-FCN throughout performed noticably worse than all other models. It also has scaling issues, as its number of parameters increased far more when going from the MobileNetv2 backbone to the Xception backbone while at the same time its performance became worse instead of better. FarSeg's performance with MobileNetv2 was decent (when ignoring the full model's relatively large size) but its performance became very unstable when using Xception. Due to these noticable performance deficits, we excluded RA-FCN and FarSeg in Fig. \ref{fig:miou}.

Generally, the best peforming models appeared to be SegForestNet (in at least one variant), PFNet, U-Net++ and U-Net. In the small and large model brackets the best performing model always is SegForestNet, SegForestNet* or PFNet. Most of the time at least one SegForestNet variant and PFNet are very close in performance. In the small model bracket, SegForestNet is the best model once, while SegForestNet* is so three times and PFNet is the best model four times (SegForestNet and PFNet share their victory on Potsdam). In the large model segment, SegForestNet* is the best model once and PFNet is the best model on all other datasets. But again, when compared to the standard deviation, the SegForestNet variants and PFNet can be said to perform equally well most of the time. In the other model size brackets, U-Net++ and U-Net almost always performed best, however, there was little competition in these brackets.

An advantage SegForestNet has over the competition is that the model is relatively small. When using MobileNetv2 for all models and comparing to PFNet, it uses almost $60\%$ fewer parameters. SegForestNet* uses up to $40\%$ fewer parameters than PFNet, depending on the dataset since the number of predicted trees depends on the numer of classes in the dataset. These numbers decrease to $20\%$ and $16\%$ respectively when moving from MobileNetv2 to Xception as a backbone. SegForestNet's advantage is even larger when compared to the U-Net variants. This, same performance at a much smaller size, proves that the prediction of polygons instead of allowing arbitrary shapes, is a useful inductive bias to learning the semantic segmentation of aerial and satellite images. With dataset size and therefore overfitting being a common problem for aerial image analysis, smaller models are beneficial as they are less prone to overfitting.

What the best model for a given dataset is also depends on the computational budget available. We used the number of model parameters as a proxy for this budget and computed how often each model appeared in the Pareto front, see Fig. \ref{fig:miou}, in each of the seven datasets. The results are shown in Table \ref{tab:pareto}. The table again highlights that the SegForestNet variants, PFNet and the U-Net variants are the top models. The seemingly good performance of FCN with MobileNetv2 is simply due to it being the smallest model in the comparison. Which model is the best suited depends on the dataset and the available computational budget.

Fig. \ref{fig:samples} shows a qualitative comparison of the best models. All of these models provide good results, while no model is perfect. SegForestNet generally performs better for cars. The other models more often overlook cars entirely and tend to have multiple distinct car segments merge into each other. Areas in shadow proof difficult for all shown models with sealed surfaces sometimes mistaken as buildings or low vegetation/grass mistaken as sealed surface. SegForestNet tends to have straighter lines where appropriate, e.g., buildings and cars, while still being able to predict sufficiently rounded shapes, e.g., trees. However, when errors occur, e.g., when low vegetation in shadow gets mistaken as a sealed surface, SegForestNet produces visible block artefacts.

\section{Conclusion}
In this paper we present a model for the semantic segmentation of aerial images using binary space partitioning trees, as well as three modifications to this model to improve its performance. The first modificatios, a refined decoder and new region map computation strategy, are aimed at improving gradients during backpropagation, while the second is a novel loss function improving the shape of the predicted segments and the last modification is an extension which enables class-specific segmentations. Taking all modifications together, our model achieves state-of-the-art performance while using up to $60\%$ fewer model parameters when using a small backbone model or up to $20\%$ fewer model parameters when using a large backbone model. In the future, we want to investigate how to let our model learn by itself what the optimal type of tree, signed distance function in each inner node, and number of trees is for a given dataset. This would reduce the number of design decisions having to be made when applying our model. Additionally, we want to  expand our model to be able to predict instance IDs in order to perform instance and/or panoptic segmentation.

\IEEEtriggeratref{51}
\bibliographystyle{plain}