Probabilistic 3d regression with projected huber distribution

We will investigate a setting where the model takes an input image which is produced by capturing a scene on a $S\times S$ sensor. The camera is a pin-hole camera without radial distoritions and have the known intrinsics

$$\begin{bmatrix}f&0&S/2\\ 0&f&S/2\\ 0&0&1\\ \end{bmatrix}$$

(1)

Coordinates for this camera are denoted $v=(x,y,z)$ where the $z$ axis is aligned with the principal axis of the camera and $v=\bar{0}$ correspond to the center of the camera.

The method will predict a probability distribution parameterized by $\theta$ based on the input image.

In this section we describe constraints which we argue a probability distribution estimated from camera data should conform to. We also motivate why these constaints should exist.

constraint 1: The model can express any variance for the projected coordinates on the image sensor and depth independently.

Formally:

$$\forall(a,A)\in\mathbb{R}^{+}\times\textbf{S}^{2}_{++}\exists\theta\text{ such that }Var[z|\theta]=a\text{ and }Var[x/z,y/z|\theta]=A$$

(2)

Motivation

For many tasks the error of depth estimation is inherently large due to scale/depth ambiguity. Therefore for objects which have a variation in size, but where the absolute size is hard to estimate there will be an inherent error proportional to the scale error and the distance to the object. Despite this estimating the position in projected coordinates can often be done accurately. For example assuming that the height of humans is between 1.5-1.9m and assuming that it is intrinsically hard to estimate the size of an unseen person an uncertainty of 13% of the distance to the person is reasonable to expect. Despite this it is possible to estimate the position of keypoints in the projected coordinate space of an image accurately, often within a few pixels.

constraint 2: The model should not try to estimate the coordinates directly, but instead estimate the depth divided by the focal length and the projected position of the object.

Formally:

$$p_{cam}(x,y,z|\theta(image);f)=p_{cam}(fx/z,fy/z,z/f|\theta(image))$$

(3)

Motivation:

It is difficult to separate the the quantities scale, depth and focal length, since the main effect of each variable is to change the apparent size of the object on the sensor. For this reason the focal length needs to be taken into account when estimating the depth of an object, otherwise the variation in focal length adds additional ambiguity. When scale/depth ambiguity is present it is not possible to estimate the absolute position for the x or y coordinate accurately, but it is often possible to do so for the projected position of the object. This motivates why a method should also estimate the location in image space. Methods which conform to these constrants are not uncommon. For example Zhen et al. (2020) estimates the depth after dividing by the focal length. It is also common to estimate the location of objects in image space.

constraint 3: The output probability should have support only for points in front of the camera.

Formally:

$$p_{cam}(x,y,z)=0\text{ if }z<0$$

(4)

Motivation:

Cameras only depict objects in front of them. The output distribution should reflect that.

constraint 4: The output probability should have support for all points in the field of view.

Formally: for a camera with focal length $f$ and sensor size $S$

$$\forall x,y,z\text{ such that }|fx/z|<S/2,|fy/z|<S/2\text{ and }z>0\text{ then }p_{cam}(x,y,z|\theta)>0$$

(5)

Motivation:

One standard way to optimize neural networks which output probability distributions is to minimize the negative log likelihood. Backpropagation is only possible if the predicted probability is larger than 0 for the ground truth position. Since there are not any guarantees what an untrained network predicts the probability has to be positive for all locations which are reasonable without considering the input data. Since a camera only depict objects in its field of view this is a sufficient condition.

constraint 5: The negative log likelihood of the position is convex. Formally:

$$-\log(p_{cam}(v|\theta))$$

(6)

is convex with respect to $v$

Motivation:

A common method to combine estimates which are expressed as probability distributions is by assuming that the two estimates have errors which are independent random variables and letting the output of the fusion be the maximum likelihood point under this assumption Mohlin et al. (2021). It is desirable if computing this combination is a convex optimization problem since it guarantees that it is possible to find the fusion quickly. In section 3.3 we show that a sufficient constraint for this property is that the negative log likelihood is convex with respect to position.

constraint 6: The probability for objects infinitely far away from the camera is 0 Formally

$$\lim\limits_{r\rightarrow\infty}\max\limits_{\|v\|=r}p(v)=0$$

(7)

$$\forall x,y\in\mathbb{R}^{2}\text{ }p(x,y,0)=0$$

(8)

Motivation:

Firstly objects of finite size are not visible if they are infinitely far away from the camera. Secondly we need this property to guarantee convergence when doing multi-view fusion.

constraint 7: $p_{cam}(v|\theta)$ is continuous with respect to $v$.

Motivation:

It is not reasonable that small changes in position change the density by a large amount.

Proposition 1: All distributions which are twice differentiable at least at one point and decompose into

$$p_{composed}(x,y,z|\theta)=p_{projected}(x/z,y/z|\theta)p_{depth}(z|\theta)$$

(9)

Do not fulfill all of the above constraints. Proof in supplementary B

Proposition: 2 In multi view fusion, if the field of view for the cameras have a non-empty intersection then if each camera produce probability estimates conforming to constraints 1-7 then a valid maximum likelihood fusion point exist and can be found by convex optimization. Proof in supplementary C

Proposition 3: Imposing a ground plane constraint is a convex optimization problem. If any point of the ground plane is in the field of view of the camera, then the constraints also guarantee that a solution exist. Proof in supplementary C

The distribution can be decomposed into one component $p_{proj}$ which mainly model the probability in the projected coordinates and another component $p_{depth}$which models the distribution in depth

The component which mainly models the probability over the projected coordinates is

$$p_{proj}(x,y|z;\mu,A)=\dfrac{1}{K_{depth}(A,\mu_{z})}\exp\left(-h\left(\left\|A\begin{bmatrix}x/z-\mu_{x}\\ y/z-\mu_{y}\end{bmatrix}\right\|_{2}\dfrac{z}{\mu_{z}}\right)\right)$$

(10)

By including $z$ in this distribution we can avoid the problem described in proposition 1. Note that conditioning on $z$ is necessary to get a proper distribution.

$\mu_{x}$, $\mu_{y}$ $\in\mathbb{R}^{2}$ model the mean position on the camera sensor. $\mu_{z}\in\mathbb{R}^{+}$ models the estimated depth of the object and $A\in\textbf{S}^{2}_{++}$ models the precision. The function $h(.)$ is a huber function.

The distribution which models the depth is

$$p_{depth}(z;\mu_{z},a)=\begin{cases}\dfrac{1}{K_{depth}(a,\mu_{z})}\exp(-a\max(z/\mu_{z},\mu_{z}/z))&\text{ if }z>0\\ 0&\text{ otherwise}\end{cases}$$

(11)

Where $\mu_{z}\in\mathbb{R}^{+}$ models the estimated depth. Note that this is the same parameter as for $p_{proj}$. $a\in\mathbb{R}^{+}$ roughly models the precision, that is a larger $a$ reduces the variance.

By combining these distributions we get the Projected Huber Distribution

$\displaystyle p_{combined}(x,y,z;\mu,A,a)=p_{proj}(x,y|z;A,\mu)p_{depth}(z;\mu_{z},a)$

(12)

for x,y,z in the cameras coordinate system.

The normalizing factors are

	$\displaystyle K_{depth}(\mu_{z},a)$	$\displaystyle=\mu_{z}(\exp(-a)/a+\Gamma(-1,a)a)$		(13)
	$\displaystyle K_{proj}(\mu_{z},A)$	$\displaystyle=\dfrac{\mu_{z}^{2}}{|A|}2\pi(1+\exp(-1/2))$		(14)
	$\displaystyle K_{combined}(\mu_{z},A,a)$	$\displaystyle=K_{depth}(\mu_{z},a)K_{proj}(\mu_{z},A)$		(15)

Where $K_{depth}$ is bounded by

$$\mu_{z}\dfrac{\exp(-a)}{a}\leq K_{depth}(\mu_{z},a)\leq\mu_{z}\exp(-a)(1/a+1)$$

(16)

This distribution will conform to all constraints, except constraint 2 which relates to how to predict the parameters of the model. We show how to predict the parameters in a way which conform to constraint 2 in subsection 4.4

proposition 5: the moments of this distribution are

	$\displaystyle E\left[\begin{bmatrix}x/z\\ y/z\end{bmatrix}|\mu_{p}\right]$	$\displaystyle=\mu_{p}$		(17)
	$\displaystyle Var[\begin{bmatrix}x/z\\ y/z\end{bmatrix}|A]$	$\displaystyle=A^{-2}\dfrac{4+3\exp(-1/2)}{2+2\exp(-1/2)}$		(18)
	$\displaystyle E\left[z\right]$	$\displaystyle=\mu_{z}\dfrac{\Gamma(2,a)/a^{2}+\Gamma(-2,a)a^{2}}{\Gamma(1,a)/a+\Gamma(-1,a)a}$		(19)
	$\displaystyle Var\left[z\right]$	$\displaystyle=\mu_{z}^{2}\dfrac{(\Gamma(3,a)/a^{3}+\Gamma(-3,a)a^{3})(\Gamma(1,a)/a+\Gamma(-1,a)a)-(\Gamma(2,a)/a^{2}+\Gamma(-2,a)a^{2})^{2}}{(\Gamma(1,a)/a+\Gamma(-1,a)a)^{2}}$		(20)

The $\Gamma$ terms are not very intuitive, but the fraction for the expected value quickly decrease to 1 while the fraction for the variance behaves similar to $1/a^{2}$ that is

	$\displaystyle E\left[z\right]$	$\displaystyle\approx\mu_{z}$		(21)
	$\displaystyle Var\left[z\right]$	$\displaystyle\approx\mu_{z}^{2}/a^{2}$		(22)

proof in supplementary F

constraint 1: Proof: From proposition 5 we see that $A$ models the uncertainty in projected coordinates, while $a$ models the uncertainty in depth coordinates.

constraint 3: $p(z)=0$ if $z<0$, from definition in equation 12

constraint 4: $p(z)>0$ if $z>0$ since the range of $\exp$ is the $(0,\infty)$ and the expression can be evaluated for all values of $x,y,z\in\mathbb{R}^{2}\times\mathbb{R}^{+}$. The field of view for a camera is a subset of the half plane $z>0$.

constraint 5: Proof sketch: Show that both $p_{proj}$ and $p_{depth}$ has a convex negative log likelihood. An affine basis change turns the argument of the $-\log(p_{proj})$ into a norm $h(\|q\|_{2})$ which is convex. $-\log(p_{depth})$ is also convex. Full proof in supplementary G.1

constraint 6: Proof sketch. For the region $z\leq 0$ the proof is trivial. Then we prove it for the region $0<z<=(x^{2}+y^{2})^{1/4}$. Here the $\exp(-h(.))$ term will go to 0. as $\|v\|_{2}\rightarrow\infty$ while the other factors are bounded. For the region $z>(x^{2}+y^{2})^{1/4}$ the factor $\exp(-a\max(z/\mu_{z},\mu_{z}/z))$ goes to 0 while the other factors are bounded. Full proof in supplementary G.2

constraint 7: Proof sketch for $z>0$ the function is trivially continuous. For $z<0$ the function is also trivially continuous. When z approaches 0 the factor $\exp(-a\max(z/\mu_{z},\mu_{z}/z))$ goes to 0. Full proof in G.3

When the parameters $\theta$ of Projected Huber Distribution are estimated by a neural network distribution it is necessary to turn the output of the neural network into a valid parameterization for Projected Huber Distribution . Valid in this sense refers to conforming to the constraints $A\in\textbf{S}^{++}$, $a\in\mathbb{R}^{+}$ and $\mu_{z}\in\mathbb{R}^{+}$, that is $A$ is positive definite while a and $\mu_{z}$ are positive.

Neural networks give outputs in $\mathbb{R}^{d_{out}}$ where $d_{out}$ is decided by the architecture of the network. These outputs are do not conform to our desired parameter constraints, unless a suitable activation is applied.

To construct this activation we start by doing a basis change to conform to constraint 2 from the motivation.

	$\displaystyle v_{p}$	$\displaystyle=\begin{bmatrix}x_{p}\\ y_{p}\end{bmatrix}=\begin{bmatrix}\dfrac{2fx}{zS}\\ \dfrac{2fy}{zS}\end{bmatrix}$		(23)
	$\displaystyle z_{p}$	$\displaystyle=\dfrac{z}{\mu_{z0}f}$		(24)

Where $\mu_{z0}$ is a bias term defined by

$$\mu_{z0}=\sqrt{\left(\max\limits_{z,f\in\textit{Dataset}}z/f\right)\left(\min\limits_{z,f\in\textit{Dataset}}z/f\right)}$$

(25)

For the purpose of proving that our loss has bounded gradients we also need the constant $D$ defined by

$$D=\sqrt{\left(\max\limits_{z,f\in\textit{Dataset}}z/f\right)\bigg{/}\left(\min\limits_{z,f\in\textit{Dataset}}z/f\right)}$$

(26)

These values can either be derived from known properties of the dataset or by computing these values for all samples in the training dataset.

$f/z$ proportional to the projected scale of an object. Even in the extreme case where the scale ambiguity is 40% and the projected size varies between 2 and 200 pixels $D$ would be less than 20. Furthermore $D$ is only used to prove that the loss has bounded gradients.

proposition 5 $\|v_{p}\|_{\infty}\leq 1$ and $1/D\leq z_{p}\ \leq D$. for points in the field of the camera with a depth conforming to equation 26.

proof Proof sketch $v_{p}$ correspond to the projected image coordinates, scaled to be between -1 and 1. $z_{p}$ is normalized by defininition 25-26. Full proof in supplementary H

Now $v_{p}$ correspond to the coordinates in the projected image and $z_{p}$ is the z coordinate after compensating for the scale change of the focal length and applying a logarithm to map $\mathbb{R}^{+}$ to $\mathbb{R}$

Define

	$\displaystyle\nu_{p}$	$\displaystyle=\begin{bmatrix}\nu_{x}\\ \nu_{y}\end{bmatrix}=\dfrac{f\mu_{z0}}{\mu_{z}}A\begin{bmatrix}\mu_{x}\\ \mu_{y}\end{bmatrix}$		(27)
	$\displaystyle\nu_{z}$	$\displaystyle=\dfrac{\mu_{z}}{\mu_{z0}f}$		(28)
	$\displaystyle B$	$\displaystyle=A\dfrac{Sf\mu_{z0}}{2\mu_{z}}$		(29)

Using a negative log likelihood of $p_{combined}$ in the original basis can be written as

	$\displaystyle L_{combined\_original}(\mu,A,a;x,y,z)$	$\displaystyle=L_{proj\_original}(\mu_{z},a,z)+L_{depth\_original}(\mu,A;x,y,z)$		(30)
	$\displaystyle L_{proj\_original}$	$\displaystyle=h\left(\left\|A\begin{bmatrix}x/z-\mu_{x}\\ y/z-\mu_{y}\end{bmatrix}\right\|_{2}\dfrac{z}{\mu_{z}}\right)-\log(|A|)+2\log(\mu_{z})$		(31)
	$\displaystyle L_{depth\_original}$	$\displaystyle=a\max(z/\mu_{z},\mu_{z}/z)+\log(\exp(-a)/a+\Gamma(-1,a)a)+\log(\mu_{z})$		(32)

Which in the new basis is

	$\displaystyle L_{combined}(\nu,B,a;v_{p},z_{p})$	$\displaystyle=L_{proj}(\nu_{p},A,v_{p},z_{p})+L_{depth}(\nu_{z},a;z_{p})$		(33)
	$\displaystyle L_{proj}(B,\nu_{p};v_{p},z_{p})$	$\displaystyle=h\left(\left\|Bv_{p}-\nu_{p}\right\|_{2}z_{p})\right)-\log(|B|)$		(34)
	$\displaystyle L_{depth}(a,\nu_{z};z_{p})$	$\displaystyle=a\max(z_{p}/\nu_{z},\nu_{z}/z_{p})+\log(\exp(-a)/a+\Gamma(-1,a)a)+\log(\nu_{z})$		(35)

In these losses we exclude terms not required for computing gradients such as terms only containing constants, $f,S,\mu_{z0}$.

$K_{depth}(a)$ and its gradients is computed by numerical integrals. The full method is described in supplementary J

Estimating these parameters will conform to constraint 2 since $\nu_{p}$ and $B$ models the position and precision for the projected coordinates while $\nu_{z}$ and $a$ models the depth and precision for the depth estimate.

Designing the activation which outputs $B$ and $\nu_{p}$ is done in the same way as Mohlin et al. (2021).

Estimating $a$ and $\nu_{z}$ is done by starting with two real numbers $w_{1}$, $w_{2}$

$$a(w_{1})=\begin{cases}w_{1}+1&\text{ if }w_{1}<0\\ exp(w_{1})&\text{ otherwise}\end{cases}$$

(36)

$$\nu_{z}(w_{1},w_{2})=\begin{cases}1+w_{2}/a(w_{1})&\text{if }w_{2}>0\\ 1/(1-w_{2}/a(w_{1}))&\text{otherwise}\end{cases}$$

(37)

With this parameterization the the loss is convex with respect to the network output when $a>1$. The gradients of the loss will be bounded with respect to $w$. Proof in K and I

Having a loss which is Lipschitz-continous should aid with stability during training since it avoids back-propagating very large gradients.

We construct a synthetic dataset by rendering objects in a similar way as Johnson et al. (2017). We choose our task to be to estimate the location of a rendered cylinder. The cylinder is rendered different reflection material properties, color and with random orientation.

The scene is rendered on a sensor of size $224\times 224$ pixels, the focal length of the camera is sampled from a uniform distribution $f\sim U(1200,2000)(pixels/m)$ The cylinder has a height of 0.2 meters and a radius of 0.1 meter. The depth is uniformly sampled from $z\sim U(3,5)m$ The x and y coordinates are sampled from $\dfrac{Sz}{f}U(-0.5,0.5)$ for each axis.

In this experimental setup the normalizing constants are

	$\displaystyle\mu_{z0}$	$\displaystyle=\sqrt{\left(\max\limits_{f\in(1200,2000),z\in(3,5)}z/f\right)\left(\min\limits_{f\in(1200,2000),z\in(3,5)}z/f\right)}\approx 2.5*10^{-3}$		(38)
	$\displaystyle D$	$\displaystyle=\sqrt{\left(\max\limits_{f\in(1200,2000),z\in(3,5)}z/f\right)\bigg{/}\left(\min\limits_{f\in(1200,2000),z\in(3,5)}z/f\right)}\approx 1.7$		(39)

The data is fitted by the method by using the rendered image as input to a standard Resnet-50 with an output dimension of 7. Instead of applying a softmax on the output we use the mapping described in section 4.4 to predict the parameters of Projected Huber Distribution . Applying a negative log likelihood loss gives equation 33 which is used as the loss when fitting the parameters of the network. The parameters are updated using Adam with default parameters. This setup is visualized in figure 2

To showcase different cases we also generate additional datasets where the cylinder is scaled with a factor uniformly sampled from U(0.8,1.2). This dataset is constructed to showcase how the method performs when a scale/depth ambiguity is present.

We also generate a dataset where we also add a cube and a sphere to both add visual complexity and introduce occlusions which incentive the model to predict different uncertainties for the case where the cylinder is occluded compared to when it is not. For this dataset no scale ambiguity is present. The cube and sphere have orientations, material properties and color sampled indepenently from the cylinder and each other.

Finally we generate a dataset where we render a scene with occluding objects and with scale ambiguity. The scaling factor is sampled independently for all objects.

Each dataset containts 8000 training samples and 2000 test samples.