LiSnowNet: Real-time Snow Removal
for LiDAR Point Clouds

Large-scale datasets for autonomous driving [1, 4, 3, 2] become more ubiquitous nowadays. Many of the recent ones contain LiDAR data collected during adverse weather conditions, targeting challenging scenarios like fog, rain, and snow.

[15] collected a dataset with points clouds in a weather-controlled climate chamber, simulating rain and fog in front of a vehicle. The dataset provides point-wise labels with three classes: clear, rain, and fog, where the “clear” points are the ones which are not caused by the adverse weather.

[13] published the CADC dataset which contains $32754$ point clouds under various snow precipitation levels. The point clouds were collected from a top-mounted HDL-32e and were grouped into sequences where each sequence has about $100$ consecutive frames. Though the CADC dataset has 3D bounding boxes for vehicles and pedestrians, we do not use them in our work as the proposed method is designed to be unsupervised.

Recently [14] released the WADS with only $1828$ point clouds which is significantly smaller than CADC. But unlike CADC, each point cloud in WADS has point-wise labels. There are $22$ classes, including “active falling snow” and “accumulated snow”.

Since the target application of the proposed network is snow removal, only CADC and WADS are used for training and evaluation.

Similar to the proposed method, the WeatherNet [15] projects point clouds onto a spherical coordinate as range images, which we go into detail in Section III-A. With this representation and the dataset released alongside with the WeatherNet, they formulated the de-noising problem as supervised multi-class pixel-wise classification. The network is constructed by a series of LiLaBlock [19] and is trained with point-wise class labels. The WeatherNet is able to remove rain and fog, but it is unclear how it performs on snow. Additionally, training the WeatherNet requires point-wise label. This prevents the WeatherNet from taking advantage of large-scale unlabeled point cloud datasets such as CADC.

The 2D median filter is successful in removing salt and pepper noise. Once we convert point clouds into range images, we can apply the median filter to them. As shown in Fig. 2, the pixels corresponding to snow are relatively darker than their nearby pixels in both channels. They are darker in the distance channel because the background is further away from the ego vehicle; in addition, snowflakes tend to absorb the energy of the beams which results in dark pixels in the intensity channel. An asymmetrical median filter can be applied to remove just the “peppers” (i.e. snow, isolated dark pixels) while leaving other pixels untouched.

The DROR (DROR) [11] filter removes the noise by first constructing a $k-d$ tree, and count the number of neighbors of each point using dynamic search radii. The idea is that the density of the points from the scene is inversely proportional to the distance, while the density of the points from snowflakes follows a different distribution. The search radius of each point is dynamically adjusted by its distance to the LiDAR center accordingly. If a point does not have enough number of neighbors within its search radius, it is classified as an outlier (i.e. snow).

The DSOR (DSOR)[12] filter also relies on nearest neighbor search, but it calculates the mean and variance of relative distances given a fixed number of neighbors of each point. In other words, it estimates the local density around each point, and outliers can then be filtered out. DSOR is shown to be slightly faster than DROR while having a similar performance in de-noising.

Both DROR and DSOR require nearest neighbor search which prevents them from being real-time. A moderate LiDAR like the HDL‐64E generates around $100k$ points at $10$Hz. Querying this many points with a $k-d$ tree of this size takes hundreds of milliseconds on modern hardware. On the other hand, WeatherNet [15] and our method work with image-like data, which will be discussed later in Section III-A. Common tools in image processing, such as FFT, DWT, and CNN, can be applied directly and very efficiently on GPUs. This allows us to develop a real-time de-noising algorithm that runs orders of magnitude faster than DROR and DSOR.

We compare our method with DROR, DSOR, and a $5$ median filter in Section V, since their implementations are publicly available.

The basic assumption of our work is that real-world clean signals are sparse under some transformation $\mathcal{T}$. Both DWT and FFT are known to generate sparse representations on real-world multi-dimensional grid data like audio and images, meaning that the Fourier and Wavelet coefficients of natural signals are mostly very close zero.

Typically, a signal corrupted by additive noise becomes less sparse in FFT and DWT than the underlying clean one. This is due to the fact that common noises, such as isolated spikes or Gaussian, have wide bandwidths in the frequency domain, as shown in Fig. 2c and 2d. Therefore de-noising is possible if the underlying clean signal is sufficiently sparse in Fourier or Wavelet. The de-noising can be formulated as maximizing the sparsity of the transformed signal.

However, solving Eq. 3 directly is computationally intensive for any reasonably sized signal. Thus we propose to use a CNN to approximate the solution and use Eq. 3 as the groundwork of our loss functions. More details about our formulation can be found in Section III-B and III-C.

Finally, not every pixel has a value due to the lack of points at certain directions like the sky and some transparent surfaces. The void pixels is processed with a series of operations, collectively denoted as $g_{2}:\mathbb{R}^{h\times w\times 2}_{+}\mapsto\mathbb{R}^{h\times w\times 2}_{+}$, before entering the network. In sequence, the operations are:

1. 1.

   $3\times 3$ dilation to fill isolated void pixels.

2. 2.

   Filling the remaining void pixels with the value $(\mu_{k}+\sigma_{k})$, where $\mu_{k}\in\mathbb{R}^{h\times 2}_{+}$ and $\sigma_{k}\in\mathbb{R}^{h\times 2}_{+}$ are the row-wise mean and standard deviation, respectively.

3. 3.

   Subtraction with a $3\times 3$ DoG to reduce the amplitude of isolated noisy pixels.

4. 4.

   $7\times 7$ average pooling to smooth out the image.

The proposed network, LiSnowNet $f_{\theta}:\mathbb{R}^{h\times w\times 2}_{+}\mapsto\mathbb{R}^{h\times w\times 2}$, is based on the MWCNN [21] with several key modifications. First, all convolution layers are replaced by residual blocks [22] with two circular convolution layers. In addition, a dropout layer is placed after the first ReLU activation in each residual block to further regularize the network. Lastly, the number of channels is dramatically reduced – the first level only has $8$ channels compared to the $160$ channels in MWCNN. The number of channels of the higher levels are also reduced proportionally. These modifications are critical for processing sparse representations (i.e. FFT and DWT coefficients) of panoramic range images in real-time.

The proposed network is designed to produce the residual image $\Delta_{k}:=f_{\theta}\left(\tilde{I}_{k}\right)$, satisfying

Let $\mathcal{F}:\mathbb{R}^{h\times w\times 2}\mapsto\mathbb{R}^{h\times w\times 2}$ be the real-valued FFT and $\mathcal{W}:\mathbb{R}^{h\times w\times 2}\mapsto\mathbb{R}^{h\times w\times 2}$ be the DWT with the Haar basis. With the formulation in Eq. 3 and 7, we designed three novel loss functions, two of which quantify the sparsity of the range image, while the third one quantifies the sparsity of the noise.

Note that $\mathcal{L}_{\mathcal{F}}$ uses the log-magnitude of the Fourier coefficients. Since $\log(1+\epsilon)\approx\epsilon$ for some small $\epsilon$, it prevents the network from focusing too much on the low-frequency part of the spectrum while pertaining the properties of the $L_{1}$ norm. The total loss function $\mathcal{L}$ can thus be constructed as

The network is trained on CADC and WADS with each dataset split into $80\%$ for training and $20\%$ for validation, where the minimum unit is a sequence. The batch size is $32$ and $8$ on the two datasets, respectively. Each training session runs $20$ epochs with the Adam optimizer and an initial learning rate of $0.001$. The learning rate is updated at the end of each epoch with a learning decay of $0.89$.

After the network is trained, the residual images $\Delta_{k}~{}\forall k$ are used for filtering. Let $\Delta^{d}_{k}\in\mathbb{R}^{h\times w}$ be the distance channel and $\Delta^{i}_{k}\in\mathbb{R}^{h\times w}$ be the intensity channel of the residual image $\Delta_{k}$ at time $k$. We define the decision boundary primarily in the residual space. A point $p$ with residuals $\delta^{d}_{k}\in\Delta^{d}_{k}$ and $\delta^{i}_{k}\in\Delta^{i}_{k}$ is classified as snow (i.e. positive) if it satisfies all the following conditions:

• •

  Being the foreground (i.e. $p$ is closer than nearby pixels).

  $$\delta^{d}_{k}>0$$

• •

  Absorbing most of the energy of LiDAR beams (i.e. $p$ is darker than nearby pixels).

  $$\delta^{i}_{k}>0$$

• •

  The residuals of $p$ are not sparse.

  $$\left(\delta^{d}_{k}\right)^{n_{d}}\cdot\left(\delta^{i}_{k}\right)^{n_{i}}>\bar{\delta}$$

As stated in Section II-C, the $L_{1}$ norm widely is used as a proxy to sparsity of real-world signals. To further justify the usage of $L_{1}$ norm as oppose to some other common norms, such as the $L_{2}$ norm, we train the same network but replace the $L_{1}$ norm with the $L_{2}$ norm in the loss functions. We label the trained networks LiSnowNet-$L_{1}$ and LiSnowNet-$L_{2}$ to reflect the norm used during training, as shown in TABLE I.

We see that LiSnowNet-$L_{2}$ already performs better than DROR in all metrics. Using the $L_{1}$ norm brings further improvements in precision and IoU while keeping a similar recall. The results indicate that using the $L_{1}$ norm as a proxy to sparsity is more effective than using the $L_{2}$ norm.

Both FFT and DWT are linear operators which are know to transform natural images to sparse representations in the frequency domain. The key difference is that FFT captures global features across the whole image, while DWT captures local features at various scales of the image.

We see that the variations of the performance remain reasonably small, with the maximum difference being $0.0251=0.9501-0.9250$ in recall. This number is way smaller than the variation between different methods presented in TABLE I, indicating that the proposed network is robust to the selection of sparse transformations if $\alpha$ and $\beta$ are properly tuned.