Multiway Non-rigid Point Cloud Registration via Learned Functional Map Synchronization

Jiahui Huang,
Tolga Birdal, Zan Gojcic,

Leonidas J. Guibas Senior Member, IEEE
and Shi-Min Hu Senior Member, IEEE
J. Huang and S.M. Hu are with the BNRist, Department of Computer Science and Technology, Tsinghua University.

Shi-Min Hu is the corresponding author.
T. Birdal and L. J. Guibas are with Stanford University.
Z. Gojcic is with Nvidia.

Transactions on Pattern Analysis and Machine Intelligence, arxiv preprint
Huang et al.: Multiway Non-rigid Registration.

We present , a novel way to jointly register multiple non-rigid shapes by synchronizing the maps that relate learned functions defined on the point clouds.
Even though the ability to process non-rigid shapes is critical in various applications ranging from computer animation to 3D digitization, the literature still lacks a robust and flexible framework to match and align a collection of real, noisy scans observed under occlusions.
Given a set of such point clouds, our method first computes the pairwise correspondences parameterized via functional maps. We simultaneously learn potentially non-orthogonal basis functions to effectively regularize the deformations, while handling the occlusions in an elegant way.
To maximally benefit from the multi-way information provided by the inferred pairwise deformation fields, we synchronize the pairwise functional maps into a cycle-consistent whole thanks to our novel and principled optimization formulation.
We demonstrate via extensive experiments that our method achieves a state-of-the-art performance in registration accuracy, while being flexible and efficient as we handle both non-rigid and multi-body cases in a unified framework and avoid the costly optimization over point-wise permutations by the use of basis function maps.
Our code is available at <https://github.com/huangjh-pub/synorim>.

3D Point Cloud, Non-rigid Registration, Functional Map Synchronization.

§ INTRODUCTION

The prevalence of reliable 3D data capture fueled countless applications impacting from movie industry to robotics.
In a wide variety of these applications, one needs to capture (in 3D) non-rigidly moving objects from multiple angles or over time <cit.>.
This leads to a dynamic, multi-scan alignment problem, further obstructed by the presence of occlusions, ambiguities, and noise.

Solving this challenging task fostered the development of a plethora of temporal, mesh-based dynamic non-rigid registration algorithms ( <cit.>).
However, these methods suffer from two main limitations:
(1) the data provided by 3D sensors hardly come in mesh format, let alone the difficulty associated with preserving the mesh topology.
Signed distance field based reconstruction methods like <cit.> can overcome some of these problems, but (2) they still assume a streaming depth map and fail to maintain correspondences, which are critical for defining non-rigid deformations.
In addition to those nuisances, a drift-free alignment almost surely demands a global optimization step, which exploits all possible loop closure constraints in the form of a graph optimization <cit.> or bundle adjustment <cit.>.

overview. (a) Input point clouds. (b) Synchronized canonical functions $\Hfunc$ (visualized via Principal Component Analysis; see <ref> for details). (c) Scene flow estimations for all pairs (here only two sets of flow, orange$\rightarrow$caribbeangreen and orange$\rightarrow$amethyst, are shown for brevity). (d) Registered (accumulated) point clouds by gathering warped points from other inputs.

In this paper, we set off in pursuit of alleviating both of these issues. In particular, we (1) leverage unstructured point cloud representations for maximal generality and flexibility; (2) assume a fully connected graph in lieu of the sequential order, incorporating drift reduction in the early stages.

Additionally, compared to rigid motion <cit.> or pure as-rigid-as-possible <cit.> deformation, our relaxed assumptions allow for more general object dynamics.

In particular, we propose  framework, synchronization for non-rigid multiway registration.
Input to , is a collection of point cloud scans containing a potentially non-rigidly deforming object. We then seek to recover a coherent and consistent 3D scene flow, originating from any source shape to a latent target shape that is to be discovered simultaneously. We start by a motion coherence observation <cit.> that pairwise correspondences between natural shapes are smooth and band-limited, nearby points in the source point cloud should map to nearby points in the target (except for topological changes and occlusions). This allows us to substitute the point-wise matches between pairs of scans by linear maps of learned smooth functions defined on the points. Such notion of a functional map was first introduced in the geometry processing community <cit.> for estimating correspondences between isometric meshes. Unfortunately, due to the lack of proper basis functions, an extension to point cloud analysis is not trivial <cit.>. In  we propose to directly learn the bases from a large-scale training corpus, either in a fully-supervised or unsupervised fashion. Such bases learned by neural networks are found to be robust to occlusions, partiality, and errors in the initialization because an end-to-end training informs the bases about the shape context while the inherent smoothness enables meaningful extrapolations of the target coordinate functions.
A subsequent functional map synchronization algorithm coerces the maps to be globally consistent, by making use of the information that individual observations are essentially samples from a common underlying surface. Based on this, we finally refine the pairwise deformations eliminating the ambiguities. An illustrative figure is shown in <ref>.

In summary, our contributions are:

* We propose, to the best of our knowledge, the first end-to-end data-driven framework to learn consistent registration between multiple, possibly partial point cloud observations of non-rigidly moving bodies / objects.

* Our deep network can learn potentially non-orthogonal functional bases on point sets without requiring to define Laplace-Beltrami operators (LBO). Functional maps between such bases can match non-isometric shapes.

* We propose a novel functional map synchronization algorithm enforcing cycle-consistency among the pairwise deformation fields, estimated in isolation. This harmonizes the 3D flows into a coherent whole and thus enables multiway registration.

We demonstrate the efficacy of our algorithm through extensive evaluations on both rigid and non-rigid scenarios, showing superior performance on all datasets.

§ RELATED WORKS

Point-set Registration

Finding reliable alignments between point clouds plays a fundamental role for many downstream tasks.
Rigid registration estimates a single transformation matrix through either heuristic searches <cit.> or learned local/global descriptors <cit.>, while common non-rigid registration techniques aim to best align the clouds under the various data terms <cit.> and regularizations <cit.> with different deformation representations <cit.>.
In this paper we consider the general non-rigid deformation scheme and also demonstrate results on hybrid ones such as multibody <cit.>.

On the other hand, scene flow describes the transition between two point clouds using a three-dimensional vector field and is a low-level task agnostic of the deformation type.
Existing techniques <cit.> handle the task via accurate modeling of the point-wise features as well as its neighborhood context, and reach a good performance for cluttered scenes or driving scenarios, but is not robust under large deformations and ambiguities <cit.>.

Function-based Correspondence
In the field of geometry processing, the use of functional techniques is a recent trend for building reliable correspondences between 3D shapes (, discretized manifold meshes).
First introduced in <cit.>, such methods compute the LBO eigenvectors as basis functions and infer a linear transformation of a subset of bases to indicate low-rank shape mapping.
Many extensions have been vastly explored such as hierarchical matching <cit.>, partial-to-full handling <cit.> or integration into deep learning frameworks <cit.> with either supervised or unsupervised <cit.> methods.
Contrarily in the domain of point set analysis, such methods are less popular despite a few <cit.>.
Notably, <cit.> also propose to learn a linearly-invariant embedding as bases, yet their learning scheme is different from ours and they are not robust to partialities or occlusions, which are commonly observed in point cloud data.
We further highlight that the idea to project natural signals to lower dimensions has also been partially explored in the vision community <cit.>.

Synchronization on 3D Geometry
Though initially established as a theory for clock systems, synchronization has now been widely used in many vision tasks such as structure from motion <cit.>, semantic segmentation <cit.>, correspondence refinement <cit.> and multiway rigid/multibody registration <cit.>.
By enforcing cycle consistencies within a system, the relative measurements are globally harmonized thanks to the averaging of local noise.
In the field of 3D geometry analysis, analyzing and synchronizing spectral functional maps <cit.> from a mesh collection are widely applied to full shape matching, either with low-rank factorization <cit.>, limit shapes extraction <cit.>, coarse-to-fine strategy <cit.> or joint point-spectrum optimization <cit.>.
In contrast, our novel synchronization formulation is tailored for the bases learned from raw point clouds that are free of the geometric impositions like (near-)isometries and directly takes correspondences into account, robustly.

§ OVERVIEW

Problem Formulation
The input to our method is a point cloud graph $\Graph := (\Vertex, \Edge)$, whose vertices $\Vertex$ represent the input set of $K$ point clouds $\Vertex := \{\X_k \in \R^{N_k \times 3}, k \in [1,K] \}$ and the edges $\Edge := \{(k,l), \X_k \in \Vertex, \X_l \in \Vertex\}$ represent the graph connectivity. By default we assume $\Graph$ to be fully-connected.
The output of our method contains all the pairwise per-point 3D flow vectors $\Fset := \{ \flow_{kl} \in \R^{N_k \times 3}, (k,l) \in \Edge \}$.
The flow vectors naturally induce the non-rigid warp field from $\X_k$ to $\X_l$ as $\warp_{kl} (\X_k) := \X_k + \flow_{kl}$, optimally aligning the given point cloud pairs by deforming the source onto the target.

We additionally encourage the cyclic consistency of the estimates $\Fset$, defined loosely as:

\begin{equation}
\warp_{k_1 k_2} \circ \dots \circ \warp_{k_{p-1} k_p} \circ \warp_{k_p k_1} = \Id, \forall (k_1, \dots k_p)\in \mathcal{C}(\Graph),
\end{equation}

where $\mathcal{C}(\Graph)$ are the set of cycles in the $\Graph$ and $\Id$ is the identity warping.
The domain of the above composed warp map is the union region $(\X_{k_1}, \dots, \X_{k_p})$.

Method. During training, our method is supervised in a pairwise fashion. We first train 172 $\descnet$ (, $\corr_\pd$) to establish putative correspondences between each point cloud pair. We then estimate a set of basis functions $\{\basis_k\}$ for each point cloud using $\basisnet$ to obtain 173 the initial functional map $\C_{kl}^0$ before refining the 3D flow estimates with $\refnet$. During test time with multiple inputs, we estimate the map set $\{\C_{kl}^0\}$ for all pairs. 174 The estimated maps are subsequently synchronized to optimize cycle consistency among the inputs. Finally, 175 3D flows are estimated from the optimized functional maps $\{\C_{kl}^\star\}$ as our final output. The registered point cloud is a fusion of all initial point clouds warped by the estimated flows.

Method Summary
Our method begins by establishing sparse point-level putative correspondences using a correspondence generator $\corr(\cdot)$ (<ref>).
In the meantime, a basis network $\basisnet$ is applied independently to each input point cloud to generate a set of basis functions defined on the points (<ref>).
Given the sparse correspondences and bases we can compute the initial pairwise functional map matrices.
The technique to recover flows from the matrices is realized through another network $\refnet$ (<ref>).
We then feed the initial map estimates into a synchronization module that jointly optimizes for all pairwise functional mappings, considering the cycle consistency (<ref>), before utilizing the $\refnet$ module once more to get the optimized deformations.
The full pipeline is illustrated in <ref>.

§ PAIRWISE FUNCTIONAL REGISTRATION

§.§ Preliminaries on Functional Maps

Given two (abstract) shapes $\Shape_k$ and $\Shape_l$ in the form of smooth continuous manifolds, a functional map $T_{kl}: L_2(\Shape_l) \mapsto L_2(\Shape_k)$ maps from the space of square-integrable real-valued functions ($L_2$-space) defined on $\Shape_l$ to $\Shape_k$.
Such an operator is proved to be linear <cit.>, , $T_{kl} (\alpha_1 f_1 + \alpha_2 f_2) = \alpha_1 T_{kl} (f_1) + \alpha_2 T_{kl} (f_2)$, where $f_1$ and $f_2$ are functions defined on $\Shape_l$ and $\alpha_1$ and $\alpha_2$ are the coefficients.

Every function $f_k$ on $\Shape_k$ (or $f_l$ on $\Shape_l$) can be represented as a linear combination of basis functions $\{\phi_{k,m}\}$ (or $\{\phi_{l,m}\}$): $f_k=\sum_m h_{k,m} \phi_{k,m}, f_l = \sum_m h_{l,m} \phi_{l,m}$, with $h_{k,m}$ and $h_{l,m}$ being the coefficients.
Moreover, most natural functions defined on shapes can be approximated by linearly combining a finite set of $M$ bases, if correctly chosen.
A functional map matrix $\C_{kl} \in \R^{M \times M}$ can be then defined as a replacement for $T_{kl}$, and satisfies $T_{kl}(\phi_{l,m_l}) = \sum_{m_k=1}^M (\C_{kl})_{m_k,m_l} \phi_{k,m_k}$, where $(\C_{kl})_{m_k,m_l}$ denotes the element at the $m_k$-th row and the $m_l$-th column.
We can then re-write the relation $f_k = T_{kl} (f_l)$ in linear algebra as $\h_k = \C_{kl} \h_l$, where $\h_k := [h_{k,1},...,h_{k,M}]^\top$, $\h_l := [h_{l,1},...,h_{l,M}]^\top$.
This re-formulation admits various computation tools available for optimization.
We refer readers to a full tutorial provided in <cit.> for more in-depth discussions.

Discretization on Point Clouds
As points are samples from the surfaces, we can define real-valued functions on each point cloud (, $\X_k$) as a column vector in $\R^{N_k}$.
Horizontally stacking the set of $M$ basis functions $\{\phi_{k,m}\}_{m=1}^M$ gives a full-rank compact basis matrix $\basis_k \in \R^{N_k \times M}$.
The pairwise linear functional map matrix $\C_{kl}$, when left-multiplied with $\basis_k$, linearly re-combines the bases from $\X_k$ and yields the transferred set of bases $\basis_l$ from $\X_l$ to $\X_k$, , $\basis_k \C_{kl} \approx \Perm_{kl} \basis_l$, where $\Perm_{kl}$ is the point-wise permutation matrix between the two point clouds.
For other arbitrary functions $\mathbf{f} \in \R^{N_k}$, their coordinates under the bases are given by $\basis_k^+ \mathbf{f}$, where $\cdot^+$ is the Moore-Penrose pseudo-inverse operator.

§.§ Generating Putative Correspondences

We estimate deep features $\desc_k \in \R^{N_k \times F}$ for each point cloud in $\Vertex$ using a sparse-convolution-based <cit.> feature descriptor network $\descnet: \X_k \rightarrow \desc_k$.
Given a pair of descriptors $\desc_k$ and $\desc_l$, we construct the soft permutation matrix $\Perm_{kl}^\de := \mathrm{softmax} (\hat{\Perm}_{kl}^\de) \in \R^{N_k \times N_l}$ and compute the scene flow as follows:

\begin{equation}
\begin{aligned}
\flow_{kl}^{\de} & := \Perm_{kl}^\de \X_l - \X_k, \\
\rev{(\hat{\Perm}_{kl}^\de)_{ij}} & := \rev{\rev{-}\frac{1}{t^\de} \lVert (\desc_k)_{i:} - (\desc_l)_{j:} \rVert},
\end{aligned}
\end{equation}

where $t^\de$ is a trainable parameter with the initial value of 1.0 and a minimum value of 0.02.
$\mathrm{softmax}(\cdot)$ performs softmax normalizations over all the rows of $\hat{\Perm}_{kl}^\de$ so that the output $\Perm_{kl}^\de$ becomes row-stochastic (, each row sums to 1).
The flow loss $\loss_\mathrm{f}$ introduced in <ref> will be used to supervise $\descnet$ until convergence.
Remarkably, such a simple strategy, similar to <cit.>, already produces the flow $\flow_{kl}^\de$.
However, such an estimation is corrupted with noise, occlusions, and inconsistencies hence performing a lot worse than our full pipeline.
This is verified through the `Ours ($\descnet$ only)' baseline in the experiments, where we simply map each point to its nearest neighbor in the space of $\desc_k$.

After $\descnet$ is trained, it is fixed and we define the correspondence generator as $\corr_\pd: \left( (\X_k, \basis_k ), (\X_l, \basis_l ) \right) \rightarrow ( \basis_k^{(kl)}, \basis_l^{(kl)} )$, which selects corresponding rows (, matched points) from the input bases $\basis_k$ and $\basis_l$ (detailed below), and produces $\basis_{k}^{(kl)} \in \R^{I_{kl}\times M}$ and $\basis_{l}^{(kl)} \in \R^{I_{kl}\times M}$, where $I_{kl}$ is the number of matches.
The row indices are determined by a nearest neighbor search performed on the descriptors $\desc_k$ and $\desc_l$ with cross-check, and the rows are matched if the L2 distance of the descriptors is smaller than a conservative threshold of 0.3.

§.§ Computing Basis Functions

The basis functions $\{\basis_k\}$ play a critical role on the power of learned representations. For triangulated shapes where connectivity information is provided, bases can be easily formed via the standard method of decomposing the graph Laplacian <cit.>.
Similarly, for point clouds one can also build a mesh structure such as a k-NN graph or intrinsic Delaunay triangulation <cit.> before computing the bases.
However, the latter can be problematic because the constructed graph would have no notion of semantics or geodesics, leading to redundant or erroneous function approximation.
Instead, we propose to learn the basis functions directly from point clouds using a sparse-convolution-based <cit.> neural network $\basisnet: \X_k \rightarrow \basis_k$. Basis functions learned in such manner are both
(1) accurate, focus only on the flow properties we are interested in, and
(2) compact with no redundancies, as will be demonstrated in <ref>, which allows achieving higher accuracy even with a small number of bases.

In order to train $\basisnet$, we compute the optimal functional maps as

\begin{equation}
\label{eq:amin}
\C_{kl}^0 := \argmin_{\C} E_{kl}(\C),
\end{equation}

\begin{equation}
\label{eq:ekl}
\begin{aligned}
E_{kl} (\C) & := \sum_{i=1}^{I_{kl}} \rho \left( \lVert \basis_{l,i:}^{(kl)} - \basis_{k,i:}^{(kl)} \C \rVert \right), \\
( \basis_k^{(kl)}, \basis_l^{(kl)} ) & = \corr \left( (\X_k, \basis_k ), (\X_l, \basis_l ) \right).
\end{aligned}
\end{equation}

Here, $\corr(\cdot)$ is a putative correspondence generator detailed earlier.
The rows of $\basis_{k}^{(kl)}$ and $\basis_{l}^{(kl)}$ are indexed via $(\cdot)_{i:}$.
During training we use the strategy similar to <cit.> by randomly instantiating $\corr$ either with the ground-truth $\corr_\gt$ (if provided) or the predicted one $\corr_\pd$ with a probability of 50%. We empirically observe the positive impact of this on final performance.
During test time, when ground-truth information is not available, we simply let $\corr = \corr_\pd$.
$\rho(\cdot)$ is the Huber robust function <cit.> used for outlier rejection that individually treats each summand $i$.
The scale of $\rho(\cdot)$ is chosen to be 0.05 and empirically we do not observe any significant difference using more advanced adaptive scales such as the median absolute deviations <cit.>.
After obtaining $\C_{kl}^0$ we feed it into the method in <ref> to compute the flow vectors, and supervise the networks using the loss function defined later in <ref>.

Learned bases visualization. We visualize for two shapes (left) a chosen subset of its $M$ bases. The basis range is normalized.

§.§ Recovering Flows from Functional Maps

Given a computed functional map matrix $\C_{kl}$, which could either be $\C_{kl}^0$ or the optimized one $\C_{kl}^\star$ from <ref>, one can again compute the soft permutation matrix $\Perm_{kl} \rev{:= \mathrm{softmax} (\hat{\Perm}_{kl})}$ between the two clouds and extract the flow as:

\begin{equation}
\label{eq:flownn}
\begin{aligned}
\flow_{kl}^{\nn} & := \Perm_{kl} \X_l - \X_k, \\ \rev{(\hat{\Perm}_{kl})_{ij}} & := \rev{\rev{-}\frac{1}{t} \lVert (\basis_k \C_{kl})_{i:} - (\basis_l)_{j:} \rVert}.
\end{aligned}
\end{equation}

However, as pointed out before, $\Perm_{kl}$ is row-stochastic and the warped position $\Perm_{kl} \X_l$ hence never extends beyond the convex hull of $\X_l$, which is detrimental for computing flows for occluded points.
An alternative way to compute the flow is to treat the target positions $\X_l$ themselves as three functions and map their coordinates via $\C_{kl}$, evaluated as:

\begin{equation}
\label{eq:flowbs}
\flow_{kl}^{\bs} := \basis_k \C_{kl} \basis_l^+ \X_l - \X_k.
\end{equation}

Note that the ground-truth mapping already implies $\basis_k \C_{kl} \approx \Perm_{kl} \basis_l$ as noted in <ref>. However, the two methods for computing flows are fundamentally different, conceptually:
$\flow_{kl}^{\bs}$ elegantly handles the issue of occlusions because the mapping allows us to map unknown positions to valid ranges given the continuity of the bases, leading to a learned extrapolation scheme via band-limited regularization.
Nevertheless, the truncated basis may over-regularize the positions in geometrically-complicated parts and limit the representation power of the flows. In contrast, $\flow_{kl}^{\nn}$ can perform well on those regions thanks to notion of point-wise correspondences.

Hence, we combine the best of the two worlds by an additional refinement network $\refnet$ which digests the concatenation of the two, along with $\X_k$, and outputs the residual:

\begin{equation}
\label{eq:flowfin}
\flow_{kl} := \flow_{kl}^{\nn} + \refnet ([\X_k, \flow_{kl}^{\bs}, \flow_{kl}^{\nn}]).
\end{equation}

The refinement network $\refnet$ is instantiated with a sparse-convolution-based network, which works in a coarse-to-fine fashion by first aligning the general cloud structure with $\flow^\bs_{kl}$ and then fixing small-scale detailed errors by considering the information from $\flow^\nn_{kl}$.
The smoothness property of such networks also helps prune predicted outliers as already demonstrated by previous works (<cit.>).
A comparison between the different flows is shown in <ref>.

§.§ Loss and Training

We train our network in two steps, where:
(1) We train the correspondence generator $\corr$ described in <ref> using a flow loss $\loss_\mathrm{f}$ only, and
(2) we jointly train $\basisnet$ and $\refnet$ (<ref> and <ref>) end-to-end adding the consistency loss $\loss_\mathrm{c}$, :

\begin{equation}
\loss_\mathrm{f} + \lambda_\mathrm{c} \loss_\mathrm{c}.
\end{equation}

Notably, the $\mathrm{argmin}$ operator in <ref> is differentiated by computing the solution with Iteratively Reweighted Least Squares (IRLS) and unrolling the re-weighting iterations.
We provide the details of this algorithm in Appendix C.

Flow Loss $\loss_\mathrm{f}$
If ground-truth annotations $\flow_{kl}^\gt$ are provided, we can directly use:

\begin{equation}
\loss_\mathrm{f,sup} := \lVert \flow_{kl} - \flow_{kl}^{\mathrm{gt}} \rVert_\fro^2,
\end{equation}

where $\lVert \cdot\rVert_\fro$ denotes Frobenius norm.
Otherwise, we use the self-supervised loss inspired from <cit.>:

\begin{equation}
\loss_\mathrm{f,unsup} := \sum_{\mathrm{type} \in \{ \mathrm{chamfer,smooth,lap} \}} \lambda_\mathrm{type} \rev{\loss_\mathrm{type}},
\end{equation}

where for conciseness, the respective loss terms $\loss_\mathrm{type}$, the weights being $\lambda_\mathrm{type}$, are defined in Appendix A.
We evaluate both settings using either $\loss_\mathrm{f,sup}$ or $\loss_\mathrm{f,unsup}$ in our experiments (see <ref> for details).
The more general $\loss_\mathrm{f}$ refers to both for convenience.

Consistency Loss $\loss_\mathrm{c}$
We impose the pairwise consistency loss <cit.> as:

\begin{equation}
\loss_\mathrm{c} := \lVert \C_{kl}^0 \C_{lk}^0 - \Id_M \rVert_\fro^2,
\end{equation}

where $\Id_M$ is the identity matrix with size $M \times M$.
This loss regularizes the compound mapping from functions on $\X_k$ to functions on $\X_l$ and its backward forms an identity, ensuring a local consistency of the learned bases.
Additionally, we tried other regularization terms over the bases, such as Laplacian commutativity <cit.>[This means to project the Laplacian operator $\mathbf{L}$ to the current bases via $\basis^+\mathbf{L}\basis$.] or descriptor preservation <cit.> but empirically no improvement is observed.

§.§ Discussions

Further Notice of $\basis$ and $\C$
Examples of the learned basis functions are visualized in <ref>.
Remarkably, the learned bases $\basis$ and the functional map matrix $\C$ are not constrained to a particular geometric structure, , orthonormality $\basis^\top \basis = \Id$ or orthogonality $\C \C^\top = \Id$, which hold only under isometric deformations with full geometry for meshes.
The coarse-to-fine structure of $\basis$ is not enforced either: the network finds its own way of representing high-frequency bases via learning small differences.
Empirically we find that enforcing any of these regularizations harms the performance.
On the other hand, unlike most existing works (, <cit.>) that formulates <ref> by projecting probe functions $\mathbf{P}_k,\mathbf{P}_l$ via $\C_{kl}^0 := \argmin_{\C} \lVert \C \basis_l^+ \mathbf{P}_l - \basis_k^+ \mathbf{P}_k \rVert $, our method is more robust due to the per-match outlier filtering and gets rid of the trade-off between the flexibility of probe functions and the compactness of the bases for flow estimation.

Special Case: Affinity Bases
For large-scale scenes with multiple moving rigid objects as often observed in the case of autonomous driving, the resulting flow can be effectively regularized by constraining the basis structure.
Specifically, if the segmentation of $S$ rigid bodies (each indexed by $s$) is known beforehand (, acquired via <cit.> or semantic segmentations), instead of learned ones, we can manually define the basis functions as:

\begin{equation}
\forall s, (\basis_{k})_{:,4s:4s+4} = \diag(\mathbf{G}_k = s) \begin{bmatrix}
\X_k & \one
\end{bmatrix},
\end{equation}

where $\diag(\cdot)$ is the diagonal matrix of an indicator vector, $\mathbf{G}_k \in \R^{N_k}$ are the rigid body indices corresponding to $\X_k$ and the final bases $\basis_k \in \R^{N_k \times 4S}$ are split into $S$ blocks (as visualized in the inset as a 2D example).
Note that the resulting map $\C_{kl} \in \R^{4S \times 4S}$ not only shows the affine transformation matrices in its $4\times 4$ sub-blocks, but the sparsity of the sub-blocks also reveals rigid-body-level permutations.
As a demonstration of the effectiveness of such bases, we provide evaluations on driving scenes in <ref>.

§ FUNCTIONAL SYNCHRONIZATION

For a geometrically coherent outcome, we further ask the loops in $\Graph$ to be closed.
Note that this step happens only during test time after the pairwise networks are trained.
In our setting of functional mappings, this amounts to enforcing the cycle-consistency <cit.>:

\begin{equation}
\C_{k_1 k_2} \dots \C_{k_{p-1} k_p}\C_{k_p k_1} = \Id, \quad \forall (k_1, \dots k_p)\in \mathcal{C}(\Graph).
\end{equation}

It can be easily shown by construction that the above formulation is equivalent to defining a set of $V$ canonical functions $\Hfunc_k = [ \h_k^1, \dots, \h_k^V ] \in \R^{M \times V} $ on each point cloud $\X_k$, expressed as the coordinates under $\basis_k$.
By enforcing the canonical functions to evaluate the same for the corresponding regions on different shapes, we define the following cycle consistency energy:

\begin{equation}
E_{\mathrm{cycle}}(\{\Hfunc_k\}, \{\C_{kl}\}) := \sum_{(k,l)\in\Edge} \lVert \Hfunc_k - \C_{kl} \Hfunc_l \rVert_F^2.
\end{equation}

Note that a naïve solution will collapse to $\Hfunc_k = \zero\,\,\forall k$, so an additional orthonormality constraint is added as $\Hfunc^\top \Hfunc = \Id_V$, where $\Hfunc := [ \Hfunc_1^\top, \dots, \Hfunc_K^\top ]^\top \in \R^{KM \times V}$.
This constraint eliminates the trivial solution by fixing the gauge and encourages a larger functional space spanned by the canonical functions.

Additionally, we define the data term as the summation of all pairwise energies:

\begin{equation}
\label{eq:data}
E_{\mathrm{data}}(\{\C_{kl}\}) := \sum_{(k,l)\in\Edge} E_{kl} (\C_{kl}),
\end{equation}

with $E_{kl}$ defined in <ref>.
Recall that during test time we set $\corr$ (needed in $E_{kl}$) to $\corr_\pd$ without any reliance on ground-truth information.
The solution to the problem

\begin{equation}
\label{eq:sync}
\argmin_{\{\Hfunc_k\}, \{\C_{kl}\}} E_{\mathrm{cycle}}(\{\Hfunc_k\}, \{\C_{kl}\}) + E_{\mathrm{data}}(\{\C_{kl}\})
\end{equation}

gives us the synchronized functional maps $\{\C_{kl}^\star\}$ as well as one possible set of canonical functions.
After the optimization, we re-use the technique described in <ref> to recover $\Fset$ as our final output, by feeding in the optimized result $\{\C_{kl}^\star\}$ instead of $\{\C_{kl}^0\}$.

In our formulation, each summand $E_{kl}$ in $E_\mathrm{data}$ of <ref> now contains the robust function $\rho(\cdot)$ whose residual not only reflects the alignment with the current pairwise map, but also the consistency among all the maps in $\Graph$, related by $\{\Hfunc_k\}$.
Remarkably, compared to existing point-based synchronization techniques (<cit.>), the use of functional maps is lightweight reducing the number of variables from the order of $O(N^2)$ to $O(M^2)$ where $M \ll N$. Besides. its additional degree-of-freedoms endow us with more flexibility in representing non-rigid transformations compared to,  <cit.>.

Numerical Optimization
Inspired by <cit.>, we use an alternating method to solve <ref>. Specifically, we iterate over optimizing $\{\Hfunc_k\}$ with fixed $\{\C_{kl}\}$ and optimizing $\{\C_{kl}\}$ with fixed $\{\Hfunc_k\}$ until convergence. Details of the solver are deferred to Appendix B.

Basis Preconditioning
While we do not regularize the columns of $\basis_k$ to be orthogonal during the training of $\basisnet$, optimizations related to it can be highly ill-posed due to the emergence of near-parallel columns.
This is due to the use of pseudo-inverse that implicitly factors out the high-frequency bases from the network output, which by nature tends to produce over-smoothed results.

Such a phenomenon, though does not affect network training without synchronization (<ref>), jeopardizes the convergence and stability when jointly optimizing <ref>.
In order to mitigate this issue, during test time, we explicitly leverage orthogonality as a means of preconditioning by performing the singular value decomposition (SVD) as $\basis_k = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^\top$ and replace $\basis_k$ with $\mathbf{U}$ as the new set of basis functions.
Note that <ref> and <ref> are calibrated during training to engulf the original basis scale from the raw network output, hence the multiplier $\mathbf{\Sigma} \mathbf{V}^\top$ should be applied back after synchronization to recover the final 3D flow vectors.

§ EXPERIMENTS

§.§ Dataset and Settings

State-of-the-art methods for non-rigid registration or matching usually evaluate on different datasets with different settings.
In order to provide a comprehensive and fair evaluation, we experiment with a full spectrum of possible input data, ranging from partial to full point clouds, from non-rigid human/clothes/animals to multi-body articulated objects, and from synthetic to real-world scenes.
To this end, we extend the following datasets <cit.> to the task of multi-way non-rigid registration (hence the prefix Multi-Point Cloud used below).
For all the datasets we keep only the points from the foreground deforming objects.

*  is sampled from CAPE <cit.> dataset, which contains scanned sequences of clothed human using a commercial 3dMD scanner. Raw scans are cropped using a fixed-height horizontal plane to remove the static ground points. We randomly sample non-consecutive $K$-frame snippets from the sequences and obtain the ground-truth flow from the fitted template mesh skinned via a Linear Blend Skinning <cit.> scheme.

*  is from the DeformingThings4D <cit.> dataset, with hundreds of dynamic sequences of humanoids and animals whose skeletons and motions are designed by experts. We extract the partial point clouds by synthetically scanning the object under different motions. We ensure that the (instance, action) tuples never overlap when we create the splits.

*  <cit.> contains partial views captured with a real RGB-D camera. The scene flows are provided by the original dataset via non-rigidly tracking the temporally densely-sampled frames. The provided object masks are used to remove the static background points. We use the original validation set as our test set and split the original training set for training and validation, without instance overlap between the splits.

*  <cit.> consists of simulated, realistic, articulated object-level models that are commonly observed in daily life. We follow the data generation and sampling strategies implemented in <cit.>, which ensure that there is no instance overlap between the splits.

Statistics of the dataset splits are shown in <ref>.

The main metric for evaluating our flow quality is the (1) , also termed as Mean Absolute Error (MAE), or End-Point Error (EPE), , the average norm of the flow error vectors over all points.
We also include metrics from the scene flow <cit.> community but adapt their thresholds to our object-level setting.
Specifically, we add:

(2) 3D Accuracy Strict (AccS), the percentage of points whose relative error $<5\%$ or $<2$cm,
(3) 3D Accuracy Relaxed (AccR), the percentage of points whose relative error $<10\%$ or $<5$cm, and
(4) Outlier Ratio, the percentage of points whose relative error $>30\%$.

In datasets with occlusions, we additionally compute all the metrics both for non-occluded points and the full point clouds.

We subsequently collect the statistics among all the $K(K-1)$ pairs of the metrics and report the mean value and standard deviation ($\pm$), the latter of which shows the consistency of the estimates among all individual pairs, as desired in the multiway setting.

Comparison of fully-supervised . Left: Color-coded per-point errors (0cm  10cm) are shown on the registered point clouds by warping all points to a selected pose. Right: Registered point clouds using our inferred 3D flows. Dashed rectangles are the views to which we warp other point clouds on the left side.

We compare the performance of  against a number of baselines grouped by the learning strategy.
Hereinafter, we use the underlined names for brevity:

* Non-learned: Coherent Point Drift method (CPD) <cit.>,
Group-wise unnormalized information potential with Rényi's entropy (GUIP)<cit.>,
and Non-Parametric Part (NPP) <cit.> using up to 25 RJ-MCMC sample steps;

* Supervised: Consistent ZoomOut (CZO) <cit.> where we replace the original Laplacian-Beltrami basis functions with the non-manifold one in <cit.> (the number of bases is set to 80),
scene flow estimation methods including FLOT <cit.> and PointPWC-Net (PWC) <cit.>, MultiBodySync (MBS) <cit.> where for non-rigid case we remove the motion segmentation module and apply only the weighted permutation synchronization with downsampled 512 points, and EmbAlign <cit.> which also estimates point cloud bases.

* Self-supervised:

PointPWC-Net (PWC-U) <cit.> with the unsupervised loss, and RMA-Net (RMA)<cit.> with the proposed differentiable rendering loss.

For  where full human bodies are captured, we additionally compare to template-based methods 3D-CODED <cit.> and PTF <cit.>, as well as deep-functional-map-based methods GeomFMaps <cit.> (supervised) and DeepShells <cit.> (self-supervised), where bases are replaced with <cit.>.
For baselines that do not have a principled way of  handling multiple inputs, we treat each pair of the input point clouds independently.
Readers are referred to the appendix for more descriptions of the baseline settings.

Parameter Settings
The networks $\descnet$, $\basisnet$, and $\refnet$ are all trained using AdamW <cit.> optimizer with the initial learning rate of $10^{-3}$ and a 30$\%$ decay every 50k samples.
$\descnet$ and $\basisnet$ are both 4-layer U-Net <cit.> with 32, 96, 64, 192 channels and 32, 64, 128, 256 channels respectively, while $\refnet$ is a 2-layer small network with 32 and 64 channels.
All networks use InstanceNorm <cit.> as normalization.
Further details about the network architectures are illustrated in Appendix D.
The initial $t$ in <ref> is chosen to be 0.1, and the weight for $\mathcal{L}_\mathrm{c}$ is $\lambda_\mathrm{c} =5 \times 10^{-3}$.
For our self-supervised training scheme, we balance the Chamfer/Laplacian/smoothness losses with the weights 1.0/5.0/1.0, respectively.
For benchmarking we choose $K=4$ and use $M=24$ for all of our experiments except for the one in <ref>.
$V$ is set to be $M-4$ for  and $M-2$ for all the other datasets.

§.§ Evaluations on

is mainly used to measure the fitting tightness and the deformation quality of the warped point clouds.
As shown in <ref>, our method achieves a significant performance boost, reaching a $21.1\%$ and $16.2\%$ lower error compared to the nearest baseline under full supervision and self-supervised scenarios, respectively.

Comparison of different methods on the  dataset. The rightmost four columns show color-coded per-point errors (0cm  10cm) on the registered point clouds by warping all points to a selected pose.

Referring to the visualizations in <ref>,  is most effective in point clouds with large changes, where traditional methods fail either due to wrong associations caused by matching ambiguity or insufficient degrees of warps.
For learning-based methods:
3D-CODED and PTF rely on predefined canonical human shape prior, which prevents them from guaranteeing a tight fitting to the input points.
CZO builds strict point-level correspondences.
Hence, it is not capable of densifying the points and is error-prone once the number of outliers from $\descnet$ is large.
Moreover, the approximate point cloud bases, used also in GeomFMaps and DeepShells, introduce further drifts.
While PWC and FLOT are trained using dense flow annotations and are free from large systematic errors, they still yield noisy estimates and result in shrunk artifacts in the detailed region.
Comparably, we use the limited-band bases to effectively regularize the motions, while preserving local rigidity and details.
By explicitly projecting point-level matches onto smooth functions, the generated warping is more accurate and robust to noise.
This is in contrast to EmbAlign, where the descriptors are globally projected onto the bases, largely limiting the distinctiveness of the point features.
In the unsupervised scenario, we are able to outperform PWC-U and RMA thanks to the superiority of our matching representation and the proposed learning strategy.
Moreover, with the help of the synchronization module, we reach a better estimation by pruning inconsistent matches from the graph, reducing the standard deviation of the L2 error by over 5.8%.

FAUST Challenge
To reinforce our results on CAPE, we test our pairwise registration module (without synchronization) on the online benchmark FAUST <cit.>, where ground-truth information is not accessible. We perform the evaluation on the intra-subject challenge, which measures the registration error between 60 pairs of shapes, where each pair consists of two full meshes of the same subject in different poses.
We directly test our pretrained model without further finetuning on the points sampled from the input meshes. The results are shown in <ref>.
Our method ranks 3rd (as of Mar. 2022) among all published methods.
The two methods that surpass ours either require manually-labelled landmarks during test time (DHNN <cit.>), or rely on artist-designed canonicalization of a predefined human model (DVM <cit.>). We are free of the above assumptions and output the more general scene flow representation applicable to different subject categories.

Handling low-overlap samples. We warp two views ($\X_k, \X_l$) to the third view ($\X_m$) using different strategies for flow computation from the functional map matrix (See <ref> for notations and analysis). Dashed rectangles highlight the differences.

§.§ Evaluations on

Scanned by a virtual camera,  contains challenging partialities and occlusions, the handling which is an important yet often missing ingredient from many methods.
Input shapes that overlap only partially can cause catastrophic failure for some methods that are incapable of modeling the flow into empty space, especially when large deformations are present.

Somewhat surprisingly, <ref>, shows that none of the SoTA methods can effectively deal with the hard problem setting presented in this dataset.
On the one hand, methods like CZO and FLOT build the point-level correlation matrix (permutations), which results in large errors on the boundaries of the point cloud.
On the other hand, the motion regularization implicitly introduced in CPD and GUIP cannot be easily adapted to the various kinds of input shapes.
Moreover, the utilization of consistency on the level of points renders MBS inapplicable under such scenarios.

Comparably, as also shown in <ref>, we remark that in our method, the learning of the basis functions

constructs a shared pattern of similar shapes under different views, hence producing a good functional extrapolation of point coordinates in a learned manner.
By not operating on the level of points, our synchronization module bypasses the occlusion problem via aligning the globally-defined basis functions, successfully reducing the error by 5.4% for the full point cloud.
Further challenging examples that can be successfully tackled by our method are visualized in <ref>.

§.§ Evaluations on

As shown in <ref>, our method consistently outperforms scene-flow based methods (19.4% / 29.2% lower  and 2.9% / 11.3% higher AccS than FLOT / PWC, resp.), even with the presence of non-negligible real-world noise, holes, and large motions.
We show the application of dynamic reconstruction from sparse views (usually hundreds of frames apart as detailed in the appendix) in <ref> and found that our learned basis functions are smooth and robust, despite the difficulties in .

§.§ Evaluations on

Quantitative and qualitative evaluations on  are shown in <ref>.

Note,  contains no occlusions and only small local deformations. This allows <cit.> to accurately estimate the point cloud bases, making CZO the most effective out of all the competitors.

MBS is specially designed to handle multi-scan multi-body case, and already decreases the PWC backbone error by half thanks to the rigidity prior. However, the permutation and segmentation synchronization modules of MBS limit the maximum number of input points to hundreds coarsening the segmentation and causing a large overall error.
As visualized earlier in <ref>, the basis functions learned from our method already contain motion cues whose patterns are similar to moving part segmentations.
Hence, we conclude that our good performance stems mainly from the explicit usage of bases and the synchronization module, which itself contributes an 11.2% and 30.3% improvement over mean  and its variance, providing a better and more consistent flow estimation.

§.§ Affinity Bases on LidarKITTI

We now conduct an experiment on the cluttered scenes such as roads in autonomous driving to demonstrate the broad applicability of . To this end, we use the LidarKITTI <cit.> dataset with ground points removed.
Note that this dataset contains only pairwise Lidar frames, hence the synchronization module is not applied.

The correspondence estimator $\corr$ is trained on FlyingThings3D <cit.> dataset using the method described in <ref> and the bases are computed according to the affinity bases strategy in <ref> with the foreground segmentation module from <cit.>.

We also train our full pairwise registration module using <cit.> with learned bases.

Quantitative comparison on  dataset. Left and right subfigures show the empirical cumulative distribution function (ECDF) and the complementary cumulative distribution function (CCDF) of  and AccS metrics. Higher curve is better. Shaded areas denote the standard deviations of the respective metrics among views.

Dynamic reconstruction on  dataset by accumulating the input point clouds. Left: RGB image of the reference frame (), which are not taken as input. Middle: Unprojected depth point cloud from . Right: Accumulated points by warping other 3 frames to .

Comparison of fully-supervised . Left: Color-coded per-point errors (0.0  0.1) are shown on the registered point clouds by warping all points to a selected pose. Right: Registered point clouds using our inferred scene flows for all input poses.

The flow estimation results are listed in <ref>. Note that here the metric ratios use the original thresholds defined in <cit.>.
Compared to our baselines, ours with affinity bases reach a better performance, showing that the basis structure can effectively introduce meaningful priors to the scene flow estimation.
Visualized in <ref>, the main error comes from the affinely-warped background due to wrong raw correspondences estimated by $\corr$.
In fact, scene flow estimation in the wild achieves improved performance mainly due to more accurate ego-motion estimation as already pointed out in, , <cit.> and our result can be further improved if such rigidity constraint is added. Though, this is out of our current scope, as our aim is to validate the effectiveness of our affinity bases.

Nevertheless, we find that in such large-scale multi-object scenes, the learned global bases are less accurate than hand-crafted ones.
This is because in such scenes, each object moves independently and each basis should thus be concentrated on a single rigid body. However, the number of bodies vary between scenes and the boundaries are hard to fit.
These properties counteract our $\basisnet$ network and are challenging to learn.
Such problems do not exist in object-level datasets where the geometric structures are simpler.

Results on LidarKITTI dataset. (a) and (b) show the best and worst result estimated using our method, respectively (s are 1.7cm and 87.7cm). caribbeangreen Source points; orange Target points; amethyst Warped points.

§.§ Ablation Studies

Number of Bases
The number of basis functions $M$ reflects the prior knowledge on the intrinsic motion variability of the input shape.
<ref> shows that as the number of bases increases, the flow estimation error first decreases, but after passing a certain point, rises significantly.
This is because the noise in the descriptors mistakenly incorporated into the learned bases start to dominate the increased representation power brought by adding more bases.

In the meantime, computing the pseudo-inverse of a wider $\basis_k$ matrix demands stronger numerical stability hampering the convergence of the entire model.
In practice, we found that selecting $M=24$ strikes a good balance between accuracy and computation time of the synchronization, which grows exponentially as $M$ increases.

Number of Frames

Our method can be easily applied to more frames ($K>4$) without re-training.
To demonstrate this, we additionally sample $K=6,8,10$ frames (under the same time interval) for  dataset and test our trained model on them.
As shown in <ref>, the accuracy of our method only drops marginally when $K$ increases.
This is because all of our networks are trained in the pairwise setting and there is no learnable component in the multi-frame synchronization module.

Influence of frame count $K$. Changes of  and AccS the number of frames $K$ are shown separately in the two plots. Dark green thin curves show the average metric (Avg.) while light green bold curves show the standard deviation (+/-). The black arrows highlight the changes of the average metrics from $K=4$ to $K=10$.

Choice of Basis
Different basis functions defined on point clouds are available and we compare  to one of the non-learned methods <cit.>.
Despite being proven by <cit.> to be optimal in representing smooth functions on manifold (in the sense that its gradient magnitudes are bounded), such bases functions are not task-tailored, whereas ours are learned to maximize the scene flow estimation performance and can work robustly under noise.

As shown in <ref>, even with 80 bases, the flow recovered using <cit.> is worse than ours with only 24 bases.
Moreover, the basis computation from a sparse convolution backbone accelerated with GPU is highly efficient and benefits from parallel implementation.
The reduced number of bases not only makes the eigen-decomposition faster, but also boosts the efficiency of subsequent synchronization module.

Training Scheme of $\basisnet$
We validate the effectiveness of the training strategies applied for $\basisnet$, including (a) the use of the robust term $\rho(\cdot)$, (b) the addition of the consistency loss $\mathcal{L}_\mathrm{c}$ and (c) the mixing strategy of $\corr_\gt$ and $\corr$.

As demonstrated in <ref>, adding the robust term gives the framework more tolerance towards noise while the consistency loss acts as an efficient regularizer and speeds up network convergence.
Apart from the performance gain from the mixing strategy, we empirically found that using only the ground-truth $\corr$ causes singularities in the network prediction which gives ill-posed scene flow when two frames are too close.

Refinement Strategy
We show that the combination of flows in <ref> is indeed helpful, especially in the occluded regions.
As illustrated in <ref>, feeding the three pieces of information altogether helps the refinement module achieve the best result.
$\X$ acts as a regularizer by reducing gross errors. Despite its decisive effect on flow accuracy, $\flow^\nn$ is, by definition, not capable of transferring points to the empty region, where $\flow^\bs$ comes to the rescue by extending the target mapped positions to the full 3D space (<ref> for visual comparisons).

To further demonstrate how our model generalizes to unseen domains, we re-split the animals part of the  dataset into two disjoint category sets $\Cset_1$ and $\Cset_2$ (detailed in the appendix).
Results in <ref> show that our model trained only on $\Cset_1$ can reach a similar performance to the model trained on $\Cset_1 \cup \Cset_2$ when tested on $\Cset_2$ (note the actions of the training and test set never overlap, denoted by $\Aset_1, \Aset_2, \Aset_3$).
This empirically verifies that our method can generalize across different animal breeds.
However, our method fails to generalize from humanoids (denoted as $\Cset_\human$) to novel categories with drastically different geometries like $\Cset_2$. As echoed in <ref>, this probably stems from the fact that our descriptor network and basis network overfit to the global structural information.

Effect of Synchronization
Synchronization links the information from multiple input point clouds exploiting cycle consistency and reduces the overall error.
Here, we ablate different optimization strategies in <ref> using  dataset.
The use of the robust term $\rho(\cdot)$, different from the one used for pairwise training, now jointly considers multi-scan cues and prunes the inconsistent outliers which is hard to detect in the pairwise setting.
Although the basis preconditioning (Basis Precond.) strategy does not result in a large performance difference at convergence, it has a positive impact on the convergence rate.
As a direct replacement to our alternating optimization scheme (Alter. Optim.), we borrow the solver from <cit.> to jointly optimize $\Hfunc$ and $\C$ on the product manifold of Stiefel and Euclidean.
At convergence, this direct optimizer reaches a slightly lower error.
Nevertheless, a lower flow error does not necessarily indicate a better convergence of <ref> and vice versa.
We conclude, though, that our optimization scheme can reach a decent performance with much better efficiency than its counterparts.

Please see <ref> for visual comparisons.

Furthermore, to show the effect under large view differences, we create a -360 dataset, where a moving camera is circling around the performing human actor, and captures 5 frames in between. All individual views are partial, but combined, they cover the full body. After an additional post-processing step of rigidly aligning the point clouds from camera space to world space, our method can successfully align different views and assemble a complete body, as illustrated in <ref>. The synchronization module closes the loop by jointly considering multiple views and corrects the alignment errors, which are hard to tackle under the pairwise setting.

Visual comparison of the synchronization effect. Flow error is shown on the warped source point cloud. The last 2 columns are computed using additional inputs (not shown here) for synchronization.

Qualitative results on -360. Input column shows overlayed point clouds from all 5 views (different colors denote different views), while Before Sync. and After Sync. show the warped point clouds obtained from our full pipeline, before and after synchronization, respectively.

§.§ Timing and Memory Analysis

We evaluate the efficiency of  on a workstation with a 3.70GHz Intel i9 CPU with a single GeForce RTX 3090 graphics card.
On average it takes 1.22s, 1.86s, 2.80s, 3.89s to generate 12, 20, 30, 42 pairs of flow vectors for $K=4, 5, 6, 7$ (, roughly 100ms/pair), respectively.
The detailed break-down of the relative time cost by each component is shown in the pie chart above:
while the refinement module (set3_red) takes more time than $\basisnet$ (set3_green) and computing $\corr$ (set3_yellow) when $K$ gets larger (because the number of pairs grow quadratically), the synchronization module (set3_purple) is always the one that dominates the computation.

We find that the training of both $\basisnet$ and $\refnet$ usually reaches a full converge within two days, while baselines like PWC take more than 5 days.
The parameters for $\basisnet$, $\descnet$ and $\refnet$ are 8.8M, 5.0M and 0.7M, respectively.
Runtime memory footprint for the entire pipeline is around 2.8GiB for $K=4$.

Non-rigid reconstruction. We show demonstrative results on the dataset from <cit.>. Reconstructions are obtained via <cit.>.

Failure Cases.
Large topological changes and matching ambiguities fail , suffering limited cross-domain generalizability.

§ CONCLUSION

In this paper we present , a novel method for multiway non-rigid point cloud registration.
In our framework, we learn the basis functions defined on the domain of point clouds and represent the registration under the framework of functional maps, endowing us with both efficiency and flexibility.

The maps are then digested by our robust synchronization module to incorporate cycle-consistency constraints into the solution, producing more accurate scene flow estimates.
A potential application of non-rigid reconstruction is shown in <ref>.
Our method leverages ideas from computer graphics community concerned with processing clean geometry and applies these ideas to the task of computer vision where occlusions and noise are inevitable.
The good performance of  on a full spectrum of cases gives hope for future works in bridging the gap between these two schools.

Despite the state-of-the-art performance,   also has some limitations which we point out in the following and illustrate in <ref>:
(1) the estimated bases sometimes overly-regularize the motion under large displacements and topological changes, causing wrong parts to be stitched together;
(2)  is less robust to handling motion ambiguities caused by severe occlusions — in the extreme case $\C$ will be ill-conditioned if the corresponding functional coordinates amount to zero; and
(3) we find that the trained network has difficulties to generalize to new domains (from humans to animals/objects). This is partially due to the global field-of-view we employed in our basis network being easily fitted to the entire object structure.
(4) Due to the large and possibly ambiguous solution space of non-rigid registration, we limit all the datasets to a relatively small view change, in line with the weakly-supervised setting in <cit.>.
Future studies on transformation in-/equi-variant backbones may help relax such an assumption.
(5) Performance-wise, the synchronization module requires $O(K^2)$ pairwise maps as input and the computational cost grows quadratically the number of views. The root cause of this is the fully-connected graph $\mathcal{G}$ we have assumed: A principled way to build and solve on a sparsified graph $\mathcal{G}$, or to use batch-wise or incremental synchronization algorithms are good directions to take.

Future Work Besides addressing aforementioned limitations,  leaves ample room for future works, including better representation of the deformation flow and incorporating more priors to capture complicated motions.

§ ACKNOWLEDGMENTS

§ ACKNOWLEDGMENT

This paper was supported by National Key R&D Program of China (project No. 2021ZD0112902), a Vannevar Bush faculty fellowship, ARL grant W911NF2120104, NSF grant IIS-1763268, and a gift from the Autodesk Corporation.

[]Jiahui Huang
received his B.S. degree in computer science and technology from Tsinghua University in 2018.
He is currently a Ph.D. candidate at the Departemnt of Computer Science and Technology, Tsinghua University.
His research interests include computer vision, robotics and computer graphics.
His personal webpage is at <https://cg.cs.tsinghua.edu.cn/people/ huangjh/>.

[]Tolga Birdal
is a postdoctoral researcher at the Geometric Computing Group of Stanford University. He has defended his masters and Ph.D. theses at the Computer Vision Group, CAMP Chair, TU Munich. He was also a Doktorand at Siemens AG. His current foci of interest involve geometric machine learning and 3D computer vision. More theoretical work is aimed at investigating the limits in geometric computing, non-Euclidean inference and principles of deep learning.
His personal webpage is at <https://www.tbirdal.me>.

[]Zan Gojcic is currently a research scientist at the NVIDIA AI Lab, which he joined in 2021. He has received his PhD degree from ETH Zurich in 2021, under the supervision of Andreas Wieser. His research interests revolve around 3D computer vision, domain adaptation, and reducing the supervision. His personal webpage is located at <https://zgojcic.github.io/>

[]Leonidas J. Guibas
received his Ph.D. degree from Stanford University in 1976, under the supervision of Donald Knuth. His main subsequent employers were Xerox PARC, MIT, and DEC/SRC. Since 1984, he has been at Stanford University, where he is a professor of computer science. His research interests include computational geometry, geometric modeling, computer graphics, computer vision, sensor networks, robotics, and discrete algorithms. He is a senior member of the IEEE and the IEEE Computer Society.

[]Shi-Min Hu
is currently a professor in Computer Science at Tsinghua University.
He received a Ph.D. degree from Zhejiang University in 1996.
His research interests include geometry processing, image & video processing, rendering, computer animation and CAD.
He has published more than 100 papers in journals and refereed conferences.
He is Editor-in-Chief of Computational Visual Media, and on the editorial boards of several journals, including Computer Aided Design and Computer & Graphics.

§ SELF-SUPERVISED FLOW LOSS

In this section we detail the self-supervised flow loss used to train our networks, including the following three terms:

* Chamfer Loss encourages $\X_k$ to move as close as possible to $\X_l$:

\begin{align}
\loss_\mathrm{chamfer} := & \sum_{ \x_k^w \in \X_k^w } \min_{ \x_l \in \X_l } \lVert \x_k^w - \x_l \rVert_2^2 \hspace{0.5em} + \\ & \sum_{\x_l \in \X_l} \min_{\x_k^w \in \X_k^w} \lVert \x_k^w - \x_l \rVert_2^2,
\end{align}

where $\X_k^w := \X_k + \flow_{kl}$ and the elements $\x$ in $\X$ refer to the points in the point set.

* Smoothness Loss encourages the local smoothness of the predicted scene flow for nearby points:

\begin{equation}
\loss_\mathrm{smooth} := \sum_{\x_{i} \in \X_k} \frac{1}{|N(\x_{i})|} \sum_{\x_{j} \in N(\x_{i})} \lVert \flow_{k,i} - \flow_{k,j} \rVert_2^2,
\end{equation}

where $(\cdot)_i$ and $(\cdot)_j$ are the point indices and $N(\x)$ is the set of $k$ nearest neighbors of $\x$.

* Laplacian Loss enforces the geometric structure depicted by the laplacian coordinate vector is preserved after warping, giving:

\begin{equation}
\loss_\mathrm{lap} := \sum_{\x_k^w \in \X_k^w} \lVert \Delta(\x_k^w) - \Delta(\x_l^\mathrm{inter}) \rVert_2^2,
\end{equation}

where $\x_l^\mathrm{inter}$ is the interpolated position of $\x_k^w$ on $\X_l$ and the laplace operator is defined as:

\begin{equation}
\Delta(\x_i) := \frac{1}{|N(\x_i)|} \sum_{\x_j \in N(\x_i)} (\x_j - \x_i).
\end{equation}

For more details please refer to <cit.>. The only difference between the loss we used and theirs is that we remove the multi-scale supervision on each hierarchy.

§ SOLVING THE SYNCHRONIZATION PROBLEM

We use an alternating method to solve <ref> as described in the main text. The optimization starts with the pairwise initialization $\{\C_{kl}^0 \}$ obtained from <ref> and iterates over the following two steps:

* Optimizing $\{\Hfunc_k\}$ with fixed $\{\C_{kl}\}$. As demonstrated in <cit.>, after fixing $\{\C_{kl}\}$, the problem becomes a generalized Rayleigh problem with the orthonormality constraint. The spectral solution $\Hfunc := [ \Hfunc_1^\top, \dots, \Hfunc_K^\top ]^\top \in \R^{KM \times V}$ are the $V$ eigenvectors of the graph connection Laplacian matrix $\Lap \in \R^{KM \times KM}$ corresponding to the $V$ smallest eigenvalues. The matrix $\Lap$ is constructed as follows:

\begin{align}
\label{eq:wgcl}
\hspace{-3mm}\Lap = \begin{bmatrix}
\Lap_1 & -\hat{\C}_{12} & \dots & -\hat{\C}_{1K} \\
-\hat{\C}_{21} & \Lap_2 & \dots & -\hat{\C}_{2K} \\
\vdots & \vdots & \ddots & \vdots \\
-\hat{\C}_{K1} & -\hat{\C}_{K2} & \dots & \Lap_K,
\end{bmatrix},
\end{align}

\begin{equation}
\begin{aligned}
\hat{\C}_{kl} & := \C_{kl} + \C_{lk}^\top, \\
\Lap_{k} & := \sum_{(k,l)\in \Edge} \Id_M + \sum_{(l,k)\in \Edge} \C_{lk}^\top \C_{lk}.
\end{aligned}
\end{equation}

* Optimizing $\{\C_{kl}\}$ with fixed $\{\Hfunc_k\}$. With $\{\Hfunc_k\}$ being considered as constant, the full problem is then decomposed into several small independent quadratic problems optimizing each $\C_{kl}, (k,l) \in \Edge$, which is tackled using the robust IRLS technique shown in <ref>. During each inner iteration the optimal solution $\C_{kl}^{[t]}$ (with $t$ being the iteration count) is given by solving a closed-form $M^2 \times M^2$ linear system:

\begin{equation}
\label{eq:c}
\begin{aligned}
\vec(\C_{kl}^{[t]}) = & \hspace{0.3em} \mathbf{\Gamma}_{kl}^{-1} \bb_{kl}, \\
\mathbf{\Gamma}_{kl} := & \hspace{0.3em} \Id \otimes \A_{kl}^\top \A_{kl} + \Hfunc_l \Hfunc_l^\top \otimes \Id, \\
\bb_{kl} := & \hspace{0.3em} \vec(\A_{kl}^\top \B_{kl} + \Hfunc_k \Hfunc_l^\top),
\end{aligned}
\end{equation}

where $\otimes$ is kronecker product, $\vec(\cdot)$ is column-wise vectorization, $\A_{kl}:=\hwt_{kl} \basis_k^{(kl)}$, $\B_{kl}:=\hwt_{kl} \basis_l^{(kl)}$ and the computation of the diagonal weight matrix $\hwt_{kl} \in \R^{I_{kl} \times I_{kl}}$ is described in <ref>.

The convergence criteria of the above optimization is defined as the mean relative change of each element in all pairwise $\C_{kl}$ being less than $3\times10^{-4}$ or exceeding the maximum number of iterations, which we set to 20.

§ DETAILS OF THE IRLS SOLVER

The solution to the robust argmin problem is used in both the end-to-end network training (<ref>) and the synchronization module (<ref>). Hence the computation must be fast and differentiable ($\frac{\partial E_{kl}}{\partial \C_{kl}^0}$ for back propagation).
We use the IRLS algorithm shown in <ref> for the forward pass and the backward pass is achieved by auto-differentiation via unrolling the iterations.
The optional input of $\C_{kl}^{[0]}$ is provided as the result from the last iteration (or simply $\C_{kl}^0$ at the beginning of the optimization) while solving the synchronization problem.

InputInputAligned bases $\basis_k^{kl}, \basis_l^{kl}$ and optional initial $\C_{kl}^{[0]}$.
OutputOutputOptimal $\C_{kl}$.
$t\gets 1$ $T$
$t = 1$ and $\C_{kl}^{[0]}$ not exists
$\hwt_{kl} \gets \Id_{I_{kl}}$.

$r_i \gets \lVert \basis_{l,i:}^{(kl)} - \basis_{k,i:}^{(kl)} \C^{[t-1]}_{kl} \rVert, \hspace{0.5em} i \in [1, I_{kl}] $,

$\hwt_{kl} \gets \diag \big( \omega(r_{1}), ..., \omega(r_{I_{kl}}) \big)^{\frac{1}{2}} $.

$\A_{kl}\gets\hwt_{kl} \basis_k^{(kl)}, \hspace{0.5em} \B_{kl}\gets\hwt_{kl} \basis_l^{(kl)}$,

Set $\C_{kl}^{[t]} \gets \A_{kl}^+ \B_{kl} $ or apply <ref> for solving <ref> or <ref>, respectively.

$\C_{kl} \gets \C_{kl}^{[T]}$.
Solving the argmin problem with IRLS.

During the weight computation of $\hwt_{kl}$, we define $\omega(\cdot)$ as the influence function of the Huber robust kernel <cit.>, formulated as:

\begin{equation}
\omega(r) := \begin{cases}
1 & \text{if $|r| < \kappa$}\\
\kappa / |r| & \text{if $|r| \geq \kappa$}
\end{cases},
\end{equation}

where the scale factor $\kappa$ is empirically chosen as 0.05.
With such a scale, not only the robustness is maintained but the resulting energy landscape also has few local minima, hence making the convergence reasonably fast within a few iterations.
For efficiency considerations, we let $T=2$ for solving <ref> and $T=1$ for solving one iteration of <ref>.
Note that although we iterate only once for the synchronization, <ref> still needs to be run multiple times during the alternating optimization scheme of $\{\Hfunc_k\}$ and $\{\C_{kl}\}$.

Network architectures. Top-left: legends and the definition of the residual block (ResBlock). Bottom-left: illustration of the 2-layer $\refnet$. Right: illustration of the 4-layer $\basisnet$. The same structure is applied in $\descnet$, and only the number of channels change from (32, 64, 128, 256) to (32, 96, 64, 192) as described in the main text.

§ NETWORK ARCHITECTURES

The architectures of all three networks ($\basisnet$, $\descnet$ and $\refnet$) used in our experiments are adopted from the work of <cit.>, and implemented using MinkowskiEngine <cit.>.
As 3D point cloud data are sparsely distributed in space, running a 3D convolution on a dense voxel grid scattered from the sparse points usually introduces large memory consumption and waste of computation over empty regions.
To this end, a `sparse' version of convolution is introduced in <cit.> that create and track voxels only on non-empty regions, maintained using a hash table to enable fast neighborhood queries.
In practice, such sparse operators are not only significantly faster and more memory-friendly, but also capable of offering impressive performance for many downstream tasks due to the convolutional nature.

In , we use the  structure, a U-Net with skip connections and residual blocks within each hierarchy, as illustrated in <ref>.
The layers and their arguments are defined as follows:

* Conv ($a \rightarrow b$, $K$, $S$): Sparse convolution with input channels $a$, output channels $b$, kernel size $K$ and stride $S$. Note $S=2$ means spatially shrinking the feature map sizes by two.

* InstanceNorm and ReLU are instance normalization <cit.> and rectified linear unit <cit.> respectively.

* ConvT ($a \rightarrow b$, $K$, $S$): Sparse transposed convolution with input channels $a$, output channels $b$, kernel size $K$ and stride $S$. Here $S=2$ means spatially dilating the feature map sizes by two.

When creating the sparse voxels for convolution, input points have to be discretized.
We use $0.01^3$ voxel size for all datasets except for LidarKITTI, where we use $0.1^3$.
As conversions back and forth between point and voxel representations are inefficient and possibly harm performance, in practice we directly work on the level of voxels, including the putative correspondences $\corr$ and the soft permutation matrices $\Perm_{kl}$.
Final per-point flow vectors are interpolated from nearby voxels via (subscript ${kl}$ omitted):

\begin{equation}
\flow_{i:} := \frac{\sum_{j\in \knn_i} \flow_{j:}^v \lVert \x_i - \x_j^v \rVert_2^{-1} }{\sum_{j\in \knn_i} \lVert \x_i - \x_j^v \rVert_2^{-1}},
\end{equation}

where $\knn_i$ returns the indices of nearby voxels to point $\x_i \in \X_i$, and $\x^v$, $\flow_{j:}^v$ are coordinates and flow vectors on the center of the $j$-th voxel, respectively.

Still, we choose to keep the descriptions in the main text based on points for brevity and clarity, and not to break the generality of our method in the sense that the backbones can be replaced with more advanced ones free of such voxelizations.

§ MISCELLANEOUS DETAILS

EmbAlign <cit.> Baseline
The original version of EmbAlign uses a simple PointNet backbone to extract the features.
To enable a fair comparison, we attempt to replace the PointNet with the sparse-conv backbone used in our method.
However, we found that replacing both the networks (, the `embedding network $\mathcal{N}$' and the second (probe) network `$\mathcal{G}$') leads to bad convergence.
By exhaustively trying out different combinations, we find that only replacing the second (probe) network `$\mathcal{G}$' leads to the best result.
This is probably due to the design choice in EmbAlign, that requires the probe functions to be explicitly projected onto the bases, preventing it from recover high-frequency features that are beneficial for accurate correspondences.

DeepShells <cit.> Baseline
For DeepShells, the use of normal data makes the method unstable because accurate estimation of normal (especially the orientation) from point cloud is non-trivial.
In the final version we use for testing, we hence remove the normal term from the iterative optimization algorithm (Smooth Shell).
We also tried to replace the SHOT feature with FPFH but the results are worse than directly feeding the absolute coordinates into the network.

Frame Indices for <ref>
The frame indices we used for the qualitative results shown in <ref> are respectively:
172 Seq008: 400, 500, 600, 800;
173 Seq009: 200, 400, 600, 800;
174 Seq012: 300, 600, 800, 900;
175 Seq004: 600, 800, 900, 1000;
176 Seq026: 200, 300, 500, 600;
177 Seq027: 300, 400, 600, 700.
Underlined frames are shown as reference () to which other frames are warped.

Categories of
The dataset is generated including the following categories:
box, dishwasher, display, door, eyeglasses, faucet, kettle, knife, laptop, lighter, microwave, oven, phone, pliers, refrigerator, safe, scissors, stapler, storageFurniture, toilet, trashCan, washingMachine, and window.

Split of the Dataset in <ref>(<ref>)
The splits we use are:
$\Cset_1$: bear, bucks, bull, cat, crocodile, dino, doggie, dragon, elephant, elk, fox, goat, hippo, hog, huskydog, lioness, moose, pig, puma, seabird, sheep, zebra;
$\Cset_2$: bunny, canie, duck, grizz, leopard, milkcow, procy, rhino, tiger;
Validation set: cattle, cetacea, chicken, deer, rabbit, raccoon, raven.