*Article* **Learning the Incremental Warp for 3D Vehicle Tracking in LiDAR Point Clouds**

**Shengjing Tian 1, Xiuping Liu 1,\*, Meng Liu 2, Yuhao Bian 1, Junbin Gao <sup>3</sup> and Baocai Yin <sup>4</sup>**


**Abstract:** Object tracking from LiDAR point clouds, which are always incomplete, sparse, and unstructured, plays a crucial role in urban navigation. Some existing methods utilize a learned similarity network for locating the target, immensely limiting the advancements in tracking accuracy. In this study, we leveraged a powerful target discriminator and an accurate state estimator to robustly track target objects in challenging point cloud scenarios. Considering the complex nature of estimating the state, we extended the traditional Lucas and Kanade (LK) algorithm to 3D point cloud tracking. Specifically, we propose a state estimation subnetwork that aims to learn the incremental warp for updating the coarse target state. Moreover, to obtain a coarse state, we present a simple yet efficient discrimination subnetwork. It can project 3D shapes into a more discriminatory latent space by integrating the global feature into each point-wise feature. Experiments on KITTI and PandaSet datasets showed that compared with the most advanced of other methods, our proposed method can achieve significant improvements—in particular, up to 13.68% on KITTI.

**Keywords:** point clouds; 3D tracking; state estimation; Siamese network; deep LK

#### **1. Introduction**

Single object tracking in point clouds aims to localize the time-varying target represented by point clouds with the supervision of a 3D bounding box in the first frame. It is a challenging yet indispensable task in many real-world applications, such as autonomous driving [1,2] and mobile robot tracking [3,4]. Generally, object tracking encompasses two subtasks, target discrimination and state estimation, which are the fundamental steps for an agent to sense the surrounding environment and conduct motion planning [5]. Over the last few years, 2D single object tracking task has been explored extensively [6–9]. Inspired by that success, many RGB-D based methods refer to the pattern of 2D tracking to conduct 3D tracking [10–13]. Although working well in the conventional 2D domain, these methods rely heavily on the RGB modality. They hence may fail when color information is lowquality or even unavailable. In this work, we focus on the 3D vehicle tracking task using deep learning in point clouds, which is still in the development stage due to several factors, such as self-occlusion, disorder, density change, and the difficulty of the state estimation.

Recently, the high-end sensors for LiDAR (light detection and ranging) have attracted much attention, since they have high accuracy and are less sensitive to weather conditions than most other sensors. More importantly, it can capture the structure of a scene by producing plenty of 3D point clouds to provide reliable geometric information from far away. However, the source data constitute the unstructured representation, where the standard convolution operation is not applicable. This hampers the application of deep

**Citation:** Tian, S.; Liu, X.; Liu, M.; Bian, Y.; Gao, J.; Yin, B. Learning the Incremental Warp for 3D Vehicle Tracking in LiDAR Point Clouds. *Remote Sens.* **2021**, *13*, 2770. https:// doi.org/10.3390/rs13142770

Academic Editor: Fahimeh Farahnakian

Received: 24 June 2021 Accepted: 9 July 2021 Published: 14 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

learning models in 3D object tracking. To overcome this barrier, some studies projected point clouds onto planes from a bird's-eye view (BEV), and then discretized them into 2D images [5,14,15].

Although they could conduct tracking by detecting frame by frame, the BEV loses abundant geometric information. Consequently, starting from source point clouds, Giancola et al. [2] proposed SiamTrack3D to learn a generically template matching function which is trained by the shape completion regularization. Qi et al. [16] leveraged deep Hough voting [17] to produce potential bounding boxes. It is worth mentioning that the aforementioned approaches merely focus on determining the best proposal from a set of object proposals. In other words, they thoroughly ignore the importance of comprehensively considering both the target discrimination and the state estimation [18].

To address this problem, we elaborately designed a 3D point cloud tracking framework with the purpose of bridging the gap between target discrimination and state estimation. It is mainly comprised of two components, a powerful target discriminator and an accurate target state estimator, which realize their respective functions through the Siamese network. The state estimation subnetwork (SES) is proposed to estimate a optimal warp using the template and candidates extracted from the tracked frame. This subnetwork extends the 2D Lucas and Kanade (LK) algorithm [19] to the 3D point cloud tracking problem by incorporating it into a deep network. However, it is non-trivial, since the Jacobian matrix from the first-order Taylor expansion cannot be calculated as in the RGB image, where the Jacobian matrix can be split into two partial terms using the chain rule. The reason is that the gradients in *x*, *y*, and *z* cannot be calculated, as connections among points are lacking in 3D point clouds. To circumvent this issue, we thoughtfully present an approximation-based solution and a learning-based solution. By integrating them into a deep network in an end-to-end manner, our state estimation subnetwork can take a pair of point clouds as inputs to predict the incremental warp parameters. Additionally, we introduce an efficient target discrimination subnetwork (TDS) to remedy the deficiency of the SES. In order to project 3D shapes into a more discriminatory latent space, we designed a new loss that takes global semantic information into consideration. During online tracking, by forcing these two components to cooperate with each other properly, our proposed model could cope with the challenging point cloud scenarios robustly.

The key contributions of our work are three-fold:


#### **2. Related Work**

#### *2.1. 2D Object Tracking*

In this paper, we focus on single object tracking problem, which can be divided into two subtasks: target discrimination and state estimation [18]. Regarding 2D visual tracking, some discrimination-based methods [7,20] have recently shown outstanding performance. In particular, the family of the correlation filter trackers [8,20] have enjoyed great popularity in the tracking research community. These methods leverage the properties of circulant matrices, which can be diagonalized by discrete Fourier transformation (DFT) to learn a classifier online. With the help of the background context and the implicit exhaustive convolution in a 2D-grid, correlation filter methods achieve impressive performance. In addition, other discrimination approaches based on deep learning have achieved competitive results on 2D tracking benchmarks [21,22]. For instance, MDNet [7] first learns general feature representation in a multi-domain way, and then it captures domain-specific information via online updating. CFNet [23] creatively integrates the correlation filter into SiamFC [24]. However, the majority of these approaches put attention into developing a powerful discriminator and simply rely on brute-force multi-scale searching to adjust the target state. There also exist several special methods such as DeepLK [19] and GONTURN [25] which are merely derived from state estimation, but these obtain merely passable performances. To sum up, the situation is that most approaches for tracking the target start only with one aspect of the two subtasks.

To mitigate this situation, Danelljan et al. [18] designed a tracker called ATOM that seamlessly combines the intersection-over-union (IoU) network [26] with the fast online classifier. Taking ATOM as baseline, Zhao et al. [27] proposed an adaptive feature fusion and obtained considerable improvements. Their method takes into account both state estimation and target discrimination, thereby achieving better accuracy and robustness. Afterwards, following the same state estimation component presented in ATOM [18], Bhat et al. [28] proposed an efficient discriminator which resorts to a target predictor employing an iterative optimization technique. Our work is motivated to bridge the gap between state estimation and target discrimination via deep network, and can be thought of as 3D counterpart of them. However, utilizing a deep network to exert the potentiality of the unstructured point cloud is still challenging in 3D tracking tasks. In this work, we present a unified framework to track the target in point cloud with a dedicated state estimation subnetwork (Section 3.2) and discrimination subnetwork (Section 3.3).

#### *2.2. 3D Point Cloud Tracking*

Point cloud is a prevalent trend for representing objects in the real 3D world. A number of advanced algorithms based on point clouds have been flourishing in object classification [29,30], detection [17,31,32], registration [19], and segmentation [33,34]. Nevertheless, point cloud tracking based on deep learning has been untapped. As we all know, many algorithms dedicated to RGB-D data have been widely studied [10,12,13,35], but most of them are mainly used to boost the 2D tracking methods with the depth channel and are not good at tackling the long-range scenario. Therefore, designing an effective pattern for tracking those partial point clouds is a very promising problem.

In the past few years, there has emerged some approaches to track the target in 3D spatial data [2,5,14,36–38]. For instance, Held et al. [36] used color-augmented search alignment algorithm to obtain the separated vehicle's velocity. Subsequently, combining shape, color, and motion information, Held et al. [37] utilized the dynamic Bayesian probabilistic model to explore the state space. However, these methods is bound to segmentation and data association algorithms. Different from them, Xiao et al. [39] simultaneously detect and track pedestrian using motion prior. All these traditional methods only obtain point segments instead of 3D orientated bounding boxes. Recently, some deep-learning-based methods infused new energy into point cloud tracking. For example, AVOD [14], FaF [5], and PIXOR [15] are designed for object detection based on BEV inputs, but can be applied tracking task in a tracking-by-detection manner. Specially for 3D tracking, Giancola et al. [2] introduced completion regularization to train a Siamese network. Subsequently, in light of the limitation of candidate box generation, Qi et al. [16] designed a point-to-box network, Zou [40] reduced redundant search space using a 3D frustum, and Fang et al. extended the region proposal network into pointNet++ [41] for 3D tracking. Nevertheless, all of above methods put more emphasis on distinguishing the target from a lot of proposals. We aimed to deal with both target discrimination and state estimation, with a dedicated Siamese network and extended LK algorithm.

#### *2.3. Jacobian Matrix Estimation*

Many tasks involve the estimation of the Jacobian matrix. As we know, the visual servoing field [42,43] usually relies on approximating the inverse Jacobian to control an agent favorably. In addition, for facial image alignment, Xiong et al. [44] proposed a supervised descent method (SDM), which avoids the calculation of the Jacobian and Hessian matrices with a sequence of learned descent directions based on hand-crafted feature. Lin et al. [45] proposed the conditional Lucas-Kanade algorithm to improve the SDM. Subsequently, Han et al. [46] dealt with the image-to-image alignment problem by jointly learning the feature representation for each pixel and partial derivatives. In this work, we innovatively estimate the Jacobian of extended LK algorithm in point cloud tracking.

#### **3. Method**

#### *3.1. Overview*

The proposed 3D point cloud tracking approach not only discriminates the target from distractors but also estimates the target state in a unified framework. Its pipeline is shown in Figure 1. Firstly, the template cropped from the reference frame and the current tracked frame are fed into the target discrimination subnetwork (TDS). It can select the best candidate in terms of the confidence score. Then, such selected candidate and the template are fed into the state estimation subnetwork (SES) to produce a incremental warp parameters Δ*ρ*. These parameters are applied to the rough state of the best candidate, leading to a new state. Next, the warped point cloud extracted by the new state is sent into the SES again, producing Δ*ρ* together with the template. This procedure is implemented iteratively until the terminal condition is satisfied. We use the same feature backbone but train the TDS and SES separately. In Section 3.2, we first present our SES in detail, which extends the LK algorithm for 2D tracking to 3D point clouds. In Section 3.3, the powerful TDS is introduced. Finally, in Section 3.4, we describe an online tracking strategy that illustrates how two components cooperate with each other.

**Figure 1.** Overview of the proposed method for 3D point cloud tracking. During online tracking, the TDS first provides a rough state of the best candidate. Afterwards, provided with the template from the reference frame, the SES produces the incremental warp of the rough state. It is implemented iteratively until the terminal condition ( Δ*ρ* < ) is satisfied. The state estimation subnetwork (SES) and the target determination subnetwork (TDS) are separately trained using the KITTI tracking dataset.

#### *3.2. State Estimation Subnetwork*

Our state estimation subnetwork (SES) is designed to learn the incremental warp parameters between the template cropped from the first frame and candidate point clouds, so as to accommodate any motion variations. We took inspiration from DeepLK [19] and

extended it to the 3D point cloud tracking task. To describe the state estimation subnetwork, we briefly revisit the inverse compositional (IC) LK algorithm [47] for 2D tracking.

The IC formulation is very ingenious and efficient because it avoids the repeated computation of the Jacobian on the warped source image. Given a template image *T* and a source image *I*, the essence of IC-LK is to solve the incremental warp parameters Δ*ρ* on *T* using sum-of-squared-error criterion. Therefore its objective function for one pixel **x** can be formulated as follows:

$$\min\_{\Delta\rho} \left\| I(\mathbf{x}) - T(\mathcal{W}(\mathbf{x}; \rho + \Delta\rho)) \right\|\_{2^{\prime}}^2 \tag{1}$$

where *<sup>ρ</sup>* ∈ *<sup>R</sup>D*×<sup>1</sup> are currently known state parameters, <sup>Δ</sup>*<sup>ρ</sup>* is the number of increments the state parameters are to go through, **x** = (*x*, *y*) are the pixel coordinates, and W is the warp function. More concretely, if one considers the location shift and scale, i.e., *ρ* = (*δx*, *δy*, *δs*) , the warp function can be written as W(**x**; *<sup>ρ</sup>*)=(*δsx* + *<sup>δ</sup>x*, *<sup>δ</sup>sy* + *<sup>δ</sup>y*) ∈ *<sup>R</sup>*2×1. Using the first-order Taylor expansion at the identity warp *ρ*0, the Equation (1) can be rewritten as

$$\min\_{\mathbf{x}, \boldsymbol{\rho}} \left\| I(\mathbf{x}) - T(\mathcal{W}(\mathbf{x}; \boldsymbol{\rho}\_0)) - \nabla T \frac{\partial \mathcal{W}(\mathbf{x}; \boldsymbol{\rho}\_0)}{\partial \boldsymbol{\rho}} \Delta \boldsymbol{\rho} \right\|\_{2}^{2} \tag{2}$$

where W(**x**, *ρ*0) = **x** is the identity mapping and ∇*T* = *∂<sup>T</sup> <sup>∂</sup><sup>x</sup>* , *<sup>∂</sup><sup>T</sup> ∂y* ∈ *<sup>R</sup>*1×<sup>2</sup> represents the image gradients. Let the Jacobian *<sup>J</sup>* <sup>=</sup> <sup>∇</sup>*<sup>T</sup> <sup>∂</sup>*W(**x**;*ρ*0) *∂ρ*<sup>∈</sup> *<sup>R</sup>*1×*D*. We hence can obtain <sup>Δ</sup>*<sup>ρ</sup>* by minimizing the above Equation (2); namely,

$$
\Delta \rho = (\boldsymbol{I}^{\top} \boldsymbol{I})^{-1} \boldsymbol{I}^{\top} [\boldsymbol{I}(\mathbf{x}) - \boldsymbol{T}(\mathcal{W}(\mathbf{x}; \rho\_0))].\tag{3}
$$

Compared with 2D visual tracking, 3D point cloud tracking has an unstructured data representation and high-dimensional search space for state parameters. Let *PT* ∈ *<sup>R</sup>*3×*<sup>N</sup>* denote the template point cloud. *PI* ∈ *<sup>R</sup>*3×*<sup>N</sup>* denotes the source point cloud in the tracked frames, which is extracted by a bounding box with inaccurate center and orientation. Note that we set the quantities of both *PT* and *PI* to *N*, and when their totals of points are less than *N*, we repeat sampling from existing points. In this work, we treat the deep network *<sup>φ</sup>* : *<sup>R</sup>*3×*<sup>N</sup>* "→ *<sup>R</sup>K*×<sup>1</sup> as a learnable "image" function. In light of this, the template point cloud *PT* and the source point cloud *PI* can obtain their descriptors using the network *φ* after transforming them into the canonical coordinate system. In addition, we regard the rigid transformation *<sup>G</sup>* ∈ *<sup>R</sup>*3×<sup>4</sup> between *PT* and *PI* as the "warp" function. In this way, we can apply the philosophy of IC-LK to the 3D point cloud tracking problem. In practice, the 3D bounding box is usually utilized to represent the target state which can be parametrized by S = (*x*, *y*, *z*, *h*, *w*, *l*, *θ*) in the LiDAR system, as shown in Figure 2. Therein, (*x*, *y*, *z*) is the target center coordinate, (*h*, *w*, *l*) represents the target size, and *θ* is the rotation angle around the y-axis. Due to the target size remaining almost unchanged in 3D spatial space, it is sufficient to focus only on the state variations in the angle and x, y, and z axes. Consequently, the transformation *G* will be represented by four warping parameters *ρ* = (*x*, *y*, *z*, *θ*) . More concretely, it can be simplified as follows:

$$G(\rho) = \begin{pmatrix} \cos(\theta) & 0 & -\sin(\theta) & x \\ 0 & 1 & 0 & y \\ \sin(\theta) & 0 & \cos(\theta) & z \end{pmatrix} \tag{4}$$

Now the state estimation problem in 3D tracking can be transformed to find *G*(*ρ*) satisfying *φ*(*G*(*ρ*) ◦ *PT*) = *φ*(*PI*), where (◦) is the warp operation on the homogeneous coordinate with *G*(*ρ*). Being analogous to the aforementioned IC-LK in Equation (2), the objective of state estimation can be written as

$$\min\_{\Delta\rho} \left\| \phi(P\_I) - \phi(P\_T) - \frac{\partial\phi(G(\rho\_0)\circ P\_T)}{\partial\rho} \Delta\rho \right\|\_2^2. \tag{5}$$

**Figure 2.** Object state representation in the sensor coordinate system. An object can be encompassed with a 3D bounding box (blue). Therein, (*x*, *y*, *z*) represents the object center location in the LiDAR coordinate system. (*h*, *w*, *l*) are the height, width, and length of object, respectively. *θ* is the radian between the motion direction and x-axis. The right-bottom also exhibits different views of object point clouds produced by LiDAR.

Similar to the Equation (2), we could solve the incremental warp <sup>Δ</sup>*<sup>ρ</sup>* = (ˆ*<sup>J</sup>* <sup>ˆ</sup>*J*)−1ˆ*<sup>J</sup>* [*φ*(*PI*) <sup>−</sup> *<sup>φ</sup>*(*PT*)] with the Jacobian matrix <sup>ˆ</sup>*<sup>J</sup>* <sup>=</sup> *∂φ*(*G*(*ρ*0) ◦ *PT*)/*∂ρ*. Unfortunately, this Jacobian matrix cannot be calculated like the classical image manner. The core obstacle is that the gradients in *x*, *y*, and *z* cannot be calculated in the scattered point clouds due to the lack of connections among points or another regular convolution structure.

We introduce two solutions to circumvent this problem. One direct solution is to approximate the Jacobian matrix through a finite difference gradient [48]. Each column of the Jacobian matrix ˆ*J* can be computed as

$$f\_i = \frac{\phi(G\_i \circ P\_T) - \phi(P\_T)}{\mu\_i} \tag{6}$$

where *μ<sup>i</sup>* are infinitesimal perturbations of the warp parameters Δ*ρ*, and *Gi* is a transformation involving only one of the warp parameters. (In other words, only the *i*-th warp parameter has a non-zero value *μi*. Please refer to Appendix A for details).

On the other hand, we treat the construction of <sup>ˆ</sup>*<sup>J</sup>* as a non-linear function <sup>F</sup> with respect to *φ*(*PT*). We hence propose an alternative: to learn the Jacobian matrix using a multi-layer perceptron, which consists of three fully-connected layers and ReLU activation functions (Figure 3). In Section 4.4, we report the comparison experiments.

Based on the above extension, we can analogously solve the incremental warp Δ*ρ* of the 3D point cloud in terms of Equation (3) as follows:

$$
\Delta \varphi = J^\dagger [\phi(P\_I) - \phi(P\_T)],
\tag{7}
$$

where *J*† = (ˆ*J* ˆ*J*)−<sup>1</sup> ˆ*J* is a Moore–Penrose inverse of ˆ*J*. Afterwards, the source point cloud cropped from the coming frame can adjust its state by the following formula

$$\mathcal{S}\_I \leftarrow \mathsf{Compose}(\mathcal{S}\_I, \Delta \mathfrak{p}), \tag{8}$$

where Compose is the inverse compositional function and S*<sup>I</sup>* is the state representation of the source point cloud *PI*.

*Network Architecture.* Figure 3 summarizes the architecture of the state estimation subnetwork. Owing to the inherent complexity of the state estimation, it is non-trivial to train a powerful estimator on the fly under the sole supervision of the first point cloud scenario. We hence train the SES offline to learn general properties for predicting the incremental warp. It is natural that we opt to adopt a Siamese architecture for producing the incremental warp parameters between the template and candidate. In particular, our network contains two branches sharing the same feature backbone, each of which consists of two blocks. As shown in Figure 3, Block-1 first generates the global descriptor. Then Block-2 consumes the aggregation of the global descriptor and the point-wise features to generate the final *K*-dimensional descriptor, based on which the Jacobian matrix ˆ*J* can be calculated. Finally, the LK module jointly considers *φ*(*PI*), *φ*(*PT*), and ˆ*J* to predict Δ*ρ*. It is notable that this module theoretically provides the fusion strategy, namely, *φ*(*PI*) − *φ*(*PT*), between two features produced by the Siamese network. Moreover, we adopt the conditional LK loss [19] to train this subnetwork in an end-to-end manner. It is formulated as

$$L\_{\rm scs} = \frac{1}{M} \sum\_{m} \mathcal{L}\_1 \Big( f^{\dagger (m)} [\phi(P\_I^{(m)}) - \phi(P\_T^{(m)})], \Delta \rho\_{\rm gt}^{(m)} \Big), \tag{9}$$

where Δ*ρgt* is the ground-truth warp parameter, L<sup>1</sup> is the smooth *L*<sup>1</sup> function [49], and *M* is the number of paired point clouds in a mini-batch. This loss can propagate back to update the network when the derivative of the batch inverse matrix is implemented.

**Figure 3.** Illustration of the proposed state estimation subnetwork (SES). Firstly, the SES extracts the shape descriptors of the paired point clouds using the designed architecture. Its details are shown in the bottom dashed box. Subsequently, the Jacobian matrix is computed by one of the solutions: approximation-based or learning-based. Its details are shown in the right dashed box. Finally, the LK module generates the incremental warp parameters Δ*ρ*.

#### *3.3. Target Discrimination Subnetwork*

In the 3D search space, how to efficiently determine the presence of the target is very critical for an agent to conduct state estimation. In this section, considering that the SES lacks discrimination ability, we present the design of a target discrimination subnetwork (TDS) to realize a strong alliance with the SES. It aims to distinguish the best candidate from distractors, thereby providing a rough target state. We leverage the matching function based method [2] to track the target. Generally, its model can be written as

$$\Psi(P\_T, P\_I) = \mathcal{g}(\psi(P\_T), \psi(P\_I)),\tag{10}$$

where <sup>Ψ</sup> is the confidence score function, *<sup>ψ</sup>* : *<sup>R</sup>*3×*<sup>N</sup>* "→ *<sup>R</sup>K*×<sup>1</sup> is the feature extractor, and *g* is a similarity metric. Under this framework, the candidate with the highest score is selected as the target.

In this work, to equip the model with global semantic information, we incorporate the intermediate features generated by the first block to point-wise features, and then pass them to the second block, as shown in Figure 4. Consequently, our TDS could project 3D partial shapes into a more discriminatory latent space, which allows an agent to distinguish the target more accurately from distractors.

**Figure 4.** The scheme of our TDS. Its feature backbone is composed of two blocks as same as the SES. Particularly, the global feature generated from Block-1 is repeated and concatenated with each point-wise feature. Afterwards, the intermediate aggregation feature is further fed into Block-2. Finally, we use the combination of similarity loss, global semantic loss, and regularization completion loss to train the TDS.

*Network Architecture.* We trained the TDS offline from scratch with an annotated KITTI dataset. Based on the Siamese network, the TDS takes paired point clouds as inputs and directly produces their similarity scores. Specifically, its feature extractor also consists of two blocks the same as in the SES. As can be seen in Figure 4, Block-1 *ψ*<sup>1</sup> generates the global descriptor, and Block-2 *ψ*<sup>2</sup> utilizes the aggregation of the global point-wise features to generate the more discriminative descriptor. As for the similarity metric *g*, we conservatively utilize hand-crafted cosine function. Finally, the similarity loss, global semantic loss, and regularization completion loss are combined in order to train this subnetwork; i.e.,

$$L\_{\rm tds} = \mathcal{L}\_2(\mathcal{g}(\psi\_2(P\_T), \psi\_2(P\_I)), \mathbf{s}) + \lambda\_1 \mathcal{L}\_2(\mathcal{g}(\psi\_1(P\_T), \psi\_1(P\_T)), \mathbf{s}) + \mathcal{L}\_3(P\_T, P\_T), \tag{11}$$

where L<sup>2</sup> is mean square error loss, *s* is the ground-truth score, *λ*<sup>1</sup> is the balance factor, and <sup>L</sup><sup>3</sup> is the completion loss for regularization [2], where *<sup>P</sup>*<sup>ˆ</sup> *<sup>T</sup>* represents each template point cloud predicted via shape completion network [2].

#### *3.4. Online Tracking*

Once trained offline, the two subnetworks can be combined for online tracking. We denote our whole framework SETD. For one scenario (frame) *Ft*, *t* ∈ {1, 2, ··· , *Q*}, where *Q* is the total number of frames, we first sampled a collection of candidates C*<sup>t</sup>* in terms of Kalman filter as SiamTrack3D does, and then passed them into the TDS, which provided the rough state representation <sup>S</sup>(0) *<sup>t</sup>* of the optimal candidate. Afterwards, it was adjusted iteratively by the SES until the termination condition was satisfied, and the best state S*<sup>t</sup>* was thus determined. Finally, the template point cloud *P*(0) *<sup>T</sup>* was updated by appending the selected candidate to itself. Algorithm 1 shows the whole process in detail.

#### **Algorithm 1:** SETD online tracking.

```
Input: Frames:{Ft|t = 1, ··· , Q}, Template: P(0)
                                          T
  Output: Target state: {St}
1 for t = 1, ··· , Q do
2 Ct ← Sample(Ft)
3 S(0)
      t ← TDS(Ct, P(t−1)
                   T )
4 P(0)
      I ← Crop(Ft, S(0)
                    t )
5 for i = 1, ··· , Niter do
6 Δρ ← SES(P(i−1)
                   I , P(t−1)
                         T )
7 if 
Δρ
 ≤  then
8 Break
9 else
10 S(i)
            t ← Compose(S(i−1)
                           t , Δρ)
11 P(i)
            I ← Crop(Ft, S(i)
                          t )
12 end
13 end
14 St ← BestState({S(i)
                       t })
15 P(t)
      T ← Update(P(t−1)
                    T , Crop(Ft, St))
16 end
17 *Crop means getting the points inside the 3D bounding box parameterized by S.
```
#### **4. Experiments**

KITTI [50] is a prevalent dataset of outdoor LiDAR point clouds. Its training set of contains 21 scenes (over 27,000 frames), and each frame has about 1.2 million points. For a fair comparison, we followed [2] to divide this dataset into a training set (scene 0–16), a validation set (scene 17–18), and a testing set (scene 19–20). In addition, to validate the effectiveness of different trackers, we also evaluated them on another large-scale point cloud dataset—PandaSet [51]. It covers complex driving scenarios, including lighting conditions at day time and night, steep hills, and dense traffic. In this dataset, more than 25 scenes were collected for testing, and the tracked instances are split into three levels (easy, middle, and hard) according to the LiDAR range.

#### *4.1. Evaluation Metrics*

To evaluate the tracking results, we adopted one-pass evaluation (OPE) [21] based on the location error and the overlap. The overlap represents the intersection-over-union (IoU) between the predicted bounding box *BP* and the corresponding ground-truth bounding box *BG*, i.e., volume(*BT* ∩ *BG*)/volume(*BT* ∪ *BG*). The location error measures the Euclidean distance between the centers of *BG* and *BT*. In this paper, the success and precision metrics are utilized as evaluation metrics. Specifically, the success metric is defined as the areaunder-curve (AUC) where the x-axis denotes the overlap threshold ranging from 0 to 1 and the y-axis refers to the ratio above the threshold. The precision metric is defined as the AUC where the x-axis represents the location error threshold ranging from 0 to 2 meters and the y-axis is the ratio below this threshold.

#### *4.2. Implementation Details*

*Training.* We conducted experiments with PyTorch and Python 3.7 on a PC equipped with a GTX 2080Ti, 32 GB RAM, and 4.00 GHz Intel Core i7-4790K CPU. When training our SES, we first sampled a pair of target shapes (template and source point clouds) from the same sequence. Additionally, these two point clouds were transformed into a canonical coordinate system according to the respective bounding box. Assuming that the target motion obeyed the Gaussian distribution, we randomly produced the warp parameters

Δ*ρgt* and applied them to the source point clouds for the supervised learning. In practice, only when the IoU between the warped bounding box and its corresponding ground truth is larger than 0.1 can this paired data be fed into the SES. The mean of Gaussian distribution was set to zero, and the covariance was a diagonal matrix diag(0.5, 0.5, 5.0). The dimension *K* of shape descriptor generated by *φ* was set to 128. The network was trained from scratch using the Adam optimizer with the batch size of 32 and the initial learning rate of 1 × <sup>10</sup><sup>−</sup>3.

Regarding our TDS, the input data were the paired point clouds transformed into a canonical coordinate system. The outputs were similarity scores. The ground-truth score is the soft distance obtained by the Gaussian function. The output dimensions *K* of *ψ* were set to 128. Our proposed loss (Equation (11)) was utilized to train it from scratch using an Adam optimizer. The batch size and initial learning rate were set to 32 and 1 × <sup>10</sup>−3, respectively. As for *λ*1, we reported its performance using several metrics in Section 4.4. The learning rates of both subnetworks were reduced via multiplying by a ratio of 0.1 when the loss of the validation set reached a plateau, and the maximum number of epochs was set to 40.

*Testing.* During the online testing phase, the tracked vehicle instance was usually specified in the first frame. When dealing with a coming frame, we exhaustively drew a set of 3D candidate boxes C*<sup>t</sup>* over the search space [2]. The number of C*<sup>t</sup>* was set to 125. The number of iterations *Niter* was set to 2, and the termination parameter was set to 1 × <sup>10</sup><sup>−</sup>5. Besides, for each frame, our SETD tracker only took about 120 ms of GPU time (50 ms for the TDS and 70 ms for the SES) to determine the final state. We did not take into account the time cost of the generation and normalization of the template and candidates, which was 300 ms of CPU time, approximately.

#### *4.3. Performance Comparison*

We first compare the proposed SETD with the baseline in relation to several attributes, such as dynamics and occlusion. Then, we evaluate the tracking performances of recent related trackers on large-scale datasets.

#### 4.3.1. Comparison with Baseline

SiamTrack3D [2] was the first method made to deal with this special task using the Siamese network and gives some referential insights via adequate ablation studies. It is a strong baseline for state-of-the-art tracking performance. We first show visualization comparison results for different attributes in Figures 5–7, and then report quantitative comparison results in Table 1.

Figure 5 shows a visualization of density variation. As can be seen, "Object-1" changed from dense to sparse, and "Object-2" varied from sparse to dense. For all these scenes, our method tracked the target accurately, whereas SiamTrack3D exhibited skewing. Figure 6 presents some tracking results for dynamic scenes. Generally, the scene is treated as dynamic when the center distance between consecutive frames is larger than 0.709 [2]. As shown in Figure 6, our method performed better than SiamTrack3D when the target moved quickly, which is attributed to the seamless integration of target discrimination and state estimation. In particular, even though the target was partly occluded and moving at high speed in the scenario presented in the second row, SETD obtained satisfying results, whereas SiamTrack3D produced greatly varying results. Figure 7 plots the tracking results from when the target suffered from different degrees of occlusion. As can be seen, whether the target was visible, partly occluded, or largely occluded, our SETD performed better than SiamTrack3D.


**Table 1.** Performance comparison with the baseline [2] in terms of several attributes.

**Figure 5.** Visual results on density change. We exhibit some key frames of two different objects. Compared with SiamTrack3D (**blue**), our SETD (**black**) has a larger overlap with the ground truth (**red**). The number after # refers to the frame ID.

**Figure 6.** Visual results in terms of dynamics. We show two dynamic scenes. When the target ran at a high speed (dynamic), our SETD obtained better results, whereas SiamTrack3D resulted in significant skewing.

**Figure 7.** Visual results in terms of occlusion. The first row shows the results of a visible vehicle. The second row shows partly occluded vehicles. The last row is a largely occluded vehicle. Our SETD performed better than SiamTrack3D in all three degrees of occlusion.

> In addition, we quantitatively compared their performances according to four conditions: visible, occluded, dynamic, and static. The success and precision metrics for said attributes are shown in Table 1. Overall, the proposed SETD fully outperformed SiamTrack3D. Specifically, compared with SiamTrack3D, SETD not only significantly improved detection by 16.45%/13.61% (success/precision) in visible scenes, but also 11.31%/14.43% in occluded scenes. Meanwhile, as shown in the last two rows of Table 1, SETD also achieved great improvements of 7.68%/7.92% when detecting dynamic scenes, and 16.54%/16.73% when detecting static scenes. These significant improvements on four types of scenes thoroughly and powerfully demonstrate the effectiveness of learning incremental warp for accurate 3D point cloud tracking.

#### 4.3.2. Comparison with Recent Methods

Apart from SiamTrack3D, we compared with other methods—AVODTrack [14], P2B [16], SiamTrack3D-RPN [52], and ICP&TDS—on the testing set. AVODTrack is a tracking-bydetection method which evolved from an advanced 3D detector, AVOD [14], by equipping it with an online association algorithm. To be more precise, it consumes point cloud BEVs and RGB images to generate a 3D detection box for every frame; then the final box is the one that has the highest IoU with the previous bounding box. P2B is a new, advanced method which integrates the target feature augmentation module into deep Hough voting [17]. SiamTrack3D-RPN evolved from SiamTrack3D by jointly learning on 2D BEV images and 3D point clouds. Furthermore, to estimate the state of moving vehicles, one may intend to obtain the motion transformation using the iterative closest point (ICP) [53], and then apply this transformation to previous bounding box. In this work, we first leveraged the proposed TDS to obtain a candidate point cloud and then run the ICP algorithm between the template and this candidate in a canonical coordinate system. This method is called ICP&TDS. In SiamTrack3D, as in the Kalman filter, the ground truth of the tracked frame is also utilized to approximate dense sampling for further testing of a tracker's discrimination ability. It applies grid search centered at the tracked ground truth. We also considered this sampling strategy for comprehensive evaluation. Note that a method with the suffix "Dense" means that it adopts this dense sampling.

Table 2 summarizes the above methods' performances on KITTI. In addition to conducting OPE of the 3D bounding box, we also present the results of the 2D BEV box, which was obtained by projecting the 3D box onto a rectangle from a bird's-eye view. As shown

in this table, SETD-Dense is superior to all other methods, given its high success and precision metrics. Specifically, the success and precision metrics of our SETD constituted 13.68% and 13.48% improvements compared with SiamTrack3D, and our SETD-Dense provided 5.12% and 6.77% improvements over SiamTrack3D-Dense. This demonstrates the validity of bridging the gap between state estimation and target discrimination. ICP&TDS and ICP&TDS-Dense obtained poor performances in these specific outdoor scenes. We deem that ICP lacks strength for partial scanned point clouds. This also proves that the proposed SES plays a critical role in the point cloud tracking task. In addition, even when using multiple modalities of RGB images and LiDAR point clouds, AVODTrack was inferior to the dense sampling models (SiamTrack3D-Dense and SETD-Dense) by large margins. P2B obtained better performances than SiamTrack3D and SETD, because P2B uses a learning procedure based on deep Hough voting to generate high quality candidates, whereas SiamTrack3D and SETD only use traditional Kalman filter sampling. Hence, better candidate generation is important for the following tracking process, and we recon that integrating a learning-based candidate generation strategy into SiamTrack3D and SETD will facilitate improving their accuracy.

**Table 2.** Performance comparison with the state-of-the-art methods on KITTI. The OPE evaluations of 3D bounding boxes and 2D BEV boxes are reported.


Table 3 reports the tracking results on PandaSet. We compare the proposed method with two advanced open-source trackers: P2B and SiamTrack3D. Their performances were obtained by running their official code on our PC. As shown in the Table 3, our SETD performed considerably better than P2B in all easy, middle, and hard sets, especially in obtaining the success/precision improvements of 6.77/9.34% with the middle set. When compared with SiamTrack3D, SETD also outperformed it by a large margin on easy and middle sets. Nevertheless, on the hard set, SETD was inferior to SiamTrack3D. The reasons were that: (1) there exist some extremely sparse objects in the hard set, which makes the SES product a worse warp parameter, (2) SiamTrack3D has a better prior because it is first trained on ShapeNet and then fine-tuned on KITTI, whereas SETD is trained only from scratch.

**Table 3.** Performance comparison with the state-of-the-art methods on PandaSet. The results on three sets of different difficulty levels are reported.


#### *4.4. Ablation Studies*

We carried out five self-contrast experiments to demonstrate the necessity of each part. *SES and TDS.* In order to prove the effectiveness of the combination of the state estimation and the target discrimination, we examined the tracking performance only using TDS or SES. TDS-only tracks the target merely via the target discrimination subnetwork, which selects a candidate bounding box with the highest confidence score. Based on the state estimation subnetwork, SES-only directly rectifies the estimated bounding box of the previous frame for tracking the target. Our SETD properly combines these two components to make up for their performance gap.

The results are shown in Table 4. According to the success metric, our SETD achieved 10.41% and 11.28% improvements in comparison with SES-only and TDS-only, respectively. As for the precision metric, SETD (69.65%) also significantly surpassed both SES-only (49.70%) and TDS-only (60.19%). These improvements of SETD highlight the importance of combining these two components. In addition, we observed that SES-only performed worse than TDS-only according to the success and precision metrics. The main reason lies in that (1) TDS-only selects the best one of many candidates generated by the Kalman filter, but SES-only directly uses the previous result while lacking discrimination; (2) the previous result used by SES often drifts due to self-occlusion and density variations, leading to farfetched warp parameters. This also proves our observation mentioned in Section 3.3: that determining the presence of the target is crucial for an agent to conduct state estimation.

**Table 4.** Self-contrast experiments evaluated by the success and precision ratio. SETD achieved the best performance.


Moreover, Figure 8 presents some tracking results obtained without or with our state estimation subnetwork. As can be seen when going through the state estimation subnetwork, some inaccurate results (blue boxes), which were predicted solely via a target discrimination subnetwork, can be adjusted towards the corresponding ground truth.

*Iteration or not.* We also investigated the effect of iteratively adjusting the target state. Specifically, we designed a variant model named Iter-non, which does not apply the iterative online tracking strategy (Algorithm 1). In other words, it directly uses the first prediction of the SES as the final state increment. As shown in Table 4, Iter-non obtained a 44.67% success ratio and a 59.13% precision ratio on the KITTI tracking dataset. Compared with SETD (53.77%/69.65%), Iter-non fell short by 9% in success and precision metrics, which proves the effectiveness of our iterative online tracking strategy. In fact, the iteration process is an explicit cascaded regression that is more effective and verifiable for tasks solved in continuous solution spaces [19,48,54].

*Jacobian Approximation or Learning.* Two solutions have been provided to tackle the Jacobi matrix issue that occurred in the SES. We can approximate it via finite difference gradient or learn it using a multi-layer perceptron. To comprehensively compare these two solutions, we plotted their loss curves during the training phase in addition to reporting the success and precision metrics on the testing set. As shown in Figure 9, the learning-based solution had a far lower cost and flatter trend than the approximation-based one. Moreover, the last row of Table 5 shows the approximation-based solution achieved success/precision of 49.93%/67.15%; the learning-based solution reached 53.77%/69.65%. The reason may be that a teachable Jacobian module could be coupled with the shape descriptor *φ*(*PT*), whereas the finite difference gradient defined by a hand-crafted formula is a hard constraint. Please refer to Appendix B for more details.

**Figure 8.** Tracking results with or without the state estimation subnetwork (SES). The black bounding boxes were obtained with SES, and the blue bounding boxes without SES. As we can see, with the help of SES, a rough state (blue boxes) can be favorably meliorated. The number after # is the frame ID.

**Figure 9.** Loss curves during the training phase, where the blue and red curves correspond to the learning-based solution and the approximation-based solution, respectively. Obviously, the former had a low cost and fast convergence.

*Descriptor Using Block-1 or Block-2.* A cascaded network architecture is proposed for the 3D point cloud tracking problem. We explored the impact of using descriptors generated by different feature blocks when extending the traditional LK algorithm to the 3D point cloud tracking task. Each column of Table 5 shows that the descriptor from Block-2 is superior that from Block-1. This benefits from the novelty that we incorporate the global semantic information into point-wise features.

**Table 5.** Performance comparison between models using different solutions for the Jacobian problem. Each row shows results using a different feature block.


*Key Parameter Analysis.* In Section 3.3, in order to robustly determine the presence of the target in a point cloud scenario, we proposed a new loss that combines similarity loss, global semantic loss, and regularization completion loss. Therein, the parameter *λ*<sup>1</sup> in Equation (11) plays a key role in the global information trade-off. In Figure 10, we compared different values of *λ*1. As we can see, it obtained the best performance in success and precision metrics when *<sup>λ</sup>*<sup>1</sup> = <sup>1</sup> × <sup>10</sup><sup>−</sup>4.

**Figure 10.** Influence of the parameter *λ*1. The OPE success and precision metrics for different values of *λ*<sup>1</sup> are reported.

#### *4.5. Failure Cases*

Figure 11 shows some failure cases of our proposed model. In this figure, "Object 1" (the first row) and "Object 2" (the second row) could not be tracked accurately by our SETD (black). The reason is that the extremely spare points could not extract an explicit pattern to discriminate target or estimate state. "Object 3" (the last row) drifted to similar distractors surrounding the target. This is because the previous bounding box, when applied to the current frame, covered similar adjacent objects due to its very fast movement.

**Figure 11.** Failure cases of our proposed method. The first and second rows show failure to track the target due to extremely sparse point clouds. The last row shows failure due to similar distractors. The number after # is the frame ID.

#### **5. Conclusions and Future Work**

This paper presents a 3D point cloud tracking framework bridging the gap between state estimation and target discrimination subnetworks. Particularly, the traditional LK algorithm has been creatively extended to the case of 3D tracking for accurate state estimation. Meanwhile, a new loss method was proposed in the hopes of providing more powerful target discrimination. Experiments on the KITTI and PandaSet datasets have shown our method significantly outperforms others. Last but not least, the ablation studies fully demonstrated the effectiveness of each part and gave some key analyses of the descriptor, iteration strategy, and Jacobian matrix calculation.

SiamTrack3D [2] is the first point cloud tracker based on the Siamese network. This method starts with state estimation and target discrimination (inspired by 2D tracker[18,19]), and extends them to 3D point cloud tracking. Although achieving promising performance, it has huge room for improvement. For example, both SiamTrack3D and SETD are struggling with the proposal extraction issue. To be specific, they obtain proposals via Kalman filter or a dense sampling strategy. In light of this, in the future, it will be very important to explore an efficient proposal extraction algorithm. Despite the recent literature [52] providing better proposals via BEV, the joint learning on 2D BEV and 3D point cloud Siamese networks even drops the final discrimination ability. Besides, a new feature backbone [41,55] is also worth studying instead of using pointNet, which is used alone in SiamTrack3D. Last but not least, it will be important to study an end-to-end network, including the flow embedding layer [56], proposal generation, similarity metric, and state refinement.

**Author Contributions:** Conceptualization, S.T. and X.L.; methodology, S.T. and X.L.; validation, S.T.; formal analysis, S.T., M.L. and X.L.; investigation, S.T. and Y.B.; writing—original draft preparation, S.T. and M.L.; writing—review and editing, M.L. and J.G.; visualization, S.T. and Y.B.; supervision, X.L. and B.Y.; project administration, S.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Natural Science Foundation of China (grant number 61976040 and U1811463) and the National Key Research and Development Program of China (grant number 2020YFB1708902).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The KITTI tracking dataset is available at http://www.cvlibs.net/ download.php?file=data\_tracking\_velodyne.zip; accessed on 24 June 2021. The PandaSet dataset is available at https://scale.com/resources/download/pandaset; accessed on 24 June 2021. For the reported results, one can obtain it by request to corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Details of Approximation-Based Solution**

After extending the LK algoritm [19] designed for 2D visual tracking to 3D point cloud tracking task, we have the following objective

$$\min\_{\Delta\rho \in R^3} \left\| \phi(P\_I) - \phi(P\_T) - \frac{\left\| \phi(G(\rho\_0) \circ P\_T) \right\|\_2}{\left\| \rho \right\|\_2} \Delta\rho \right\|\_2^2. \tag{A1}$$

As the warp parameters *<sup>ρ</sup>* <sup>∈</sup> *<sup>R</sup>*3×<sup>1</sup> and *<sup>φ</sup>*(*PT*) <sup>∈</sup> *<sup>R</sup>K*×1, the Jacobian matrix <sup>ˆ</sup>*<sup>J</sup>* in the Equation (A1) belongs to *RK*×3. The formula of the finite difference gradient [48] is as follows:

$$\hat{f}\_i = \frac{\phi(G\_i \circ P\_T) - \phi(P\_T)}{\mu\_i}. \tag{A2}$$

We use infinitesimal perturbations to approximate each column ˆ*Ji* of ˆ*J*. Therein, *Gi*, *i* = 1, 2, 3, 4 corresponds to the transformation that is obtained by only perturbing the *i*-th warp parameter, which can be formulated as

$$\mathbf{G}\_1 = \begin{pmatrix} 1 & 0 & 0 & \mu\_1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}, \mathbf{G}\_2 = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & \mu\_2 \\ 0 & 0 & 1 & 0 \end{pmatrix}, \mathbf{G}\_3 = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & \mu\_3 \end{pmatrix}, \mathbf{G}\_4 = \begin{pmatrix} \cos(\mu\_4) & 0 & -\sin(\mu\_4) & 0 \\ 0 & 1 & 0 & 0 \\ \sin(\mu\_4) & 0 & \cos(\mu\_4) & 0 \end{pmatrix}. \tag{A3}$$

When training the SES by approximation-based solution, we set the infinitesimal perturbations *μ<sup>i</sup>* to 0.1.

#### **Appendix B. Analysis of the Jacobian Module**

Our loss is defined as follows:

$$L\_{\rm scs} = \frac{1}{M} \sum\_{m} \mathcal{L}\_1(f^{\dagger(m)}[\phi(P\_I^{(m)}) - \phi(P\_T^{(m)})], \Delta \rho\_{\mathcal{g}^{\dagger}}^{(m)}).\tag{A4}$$

To enable the state estimation network to be trained in an end-to-end way, the differentiation of the Moore-Penrose inverse in Equation (A4) needs to be derived as [19] did. Concretely, the partial derivative of smooth *L*<sup>1</sup> function over the feature component *φi*(*PI*) can be written as

$$\frac{\partial \mathcal{L}\_1}{\partial \phi\_i(P\_I)} = \nabla \mathcal{L}\_1 f^\dagger \delta\_{i\nu} \tag{A5}$$

where *<sup>δ</sup><sup>i</sup>* ∈ {0, 1}*K*×<sup>1</sup> is one-hot vector and ∇L<sup>1</sup> is the derivative of the smooth *<sup>L</sup>*<sup>1</sup> loss. Besides, the partial derivative of smooth *L*<sup>1</sup> function over the *φi*(*PT*) is

$$\frac{\partial \mathcal{L}\_1}{\partial \phi\_i(P\_T)} = \nabla \mathcal{L}\_1 \left( \frac{\partial f^\dagger}{\partial \phi\_i(P\_T)} [\phi(P\_I) - \phi(P\_T)] - f^\dagger \delta\_i \right). \tag{A6}$$

Therein, the key step is to obtain the differentiation of *J*† = (ˆ*J* ˆ*J*)−<sup>1</sup> ˆ*J* . By the chain rule, it can be written as

$$\frac{\partial f^{\dagger}}{\partial \phi\_{i}(P\_{T})} = \frac{\partial (f^{\top}f)^{-1}}{\partial \phi\_{i}(P\_{T})} f^{\top} + (f^{\top}f)^{-1} \frac{\partial f^{\top}}{\partial \phi\_{i}(P\_{T})},\tag{A7}$$

where

$$\frac{\partial(\hat{f}^{\top}\hat{f})^{-1}}{\partial\phi\_{i}(P\_{T})} = -(\hat{f}^{\top}\hat{f})^{-1} \left(\hat{f}^{\top}\frac{\partial\hat{f}}{\partial\phi\_{i}(P\_{T})} + \frac{\partial\hat{f}^{\top}}{\partial\phi\_{i}(P\_{T})}\hat{f}\right)(\hat{f}^{\top}\hat{f})^{-1}.\tag{A8}$$

According to the above equation, the derivative of a batch inverse matrix can be implemented in PyTorch such that the SES network can be trained in an end-to-end manner.

In this work, we present two solutions for calculating the Jacobian in our paper. Here we give their back-propagation formulae to deeply compare them with each other. When using multi-layer perceptron <sup>F</sup> (learning-based) to calculate the Jacobian, the elements of <sup>ˆ</sup>*<sup>J</sup>* are related to each component of *φ*(*PT*). We hence have

$$\frac{\partial \mathcal{f}}{\partial \theta\_{\dot{j}}^{\Phi}} = \frac{\partial \mathcal{f}}{\partial \phi(P\_T)} \frac{\partial \phi(P\_T)}{\partial \theta\_{\dot{j}}^{\Phi}},\tag{A9}$$

where *<sup>∂</sup>* <sup>ˆ</sup>*<sup>J</sup> ∂φi*(*PT*) <sup>=</sup> *<sup>∂</sup>*F(*φ*(*PT*)) *∂φi*(*PT*) is adaptively updated. When using finite difference gradient (approximition-based), we have

> *∂* ˆ*J ∂* ˆ*J*

$$\frac{\partial \hat{f}}{\partial \theta\_{\hat{j}}^{\Phi}} = \frac{\partial \hat{f}}{\partial \phi(P\_{\Gamma})} \frac{\partial \phi(P\_{\Gamma})}{\partial \theta\_{\hat{j}}^{\Phi}} + \sum\_{k} \frac{\partial \hat{f}}{\partial \phi(G\_{k} \circ P\_{\Gamma})} \frac{\partial \phi(G\_{k} \circ P\_{\Gamma})}{\partial \theta\_{\hat{j}}^{\Phi}},\tag{A10}$$

where

$$
\frac{\partial \hat{f}}{\partial \phi\_i(P\_T)} = \begin{pmatrix} 0 & \cdots & 0 \\ \vdots & & \vdots \\ -\frac{1}{\mu\_1} & \cdots & -\frac{1}{\mu\_4} \\ \vdots & & \vdots \\ 0 & \cdots & 0 \end{pmatrix},
\tag{A11}
$$

and

$$\frac{\partial \hat{\jmath}}{\partial \phi\_i (G\_k \circ P\_T)} = (a\_{mn})\_\prime a\_{mn} = \begin{cases} \frac{1}{\mu\_k} & \text{if } m = i\,, n = k;\\ 0, & \text{otherwise.} \end{cases} \tag{A12}$$

As can be seen from the above formulae, *<sup>∂</sup>* <sup>ˆ</sup>*<sup>J</sup> ∂φi*(*PT*) in the finite difference is fixed while the learning function is updated adaptively in the multi-layer perception. Thus, the learning-based solution may be more easily coupled with the feature extractor *φ*(*PT*) than the approximation-based one.

#### **References**

