Zhengyang Shen
email: zyshen@cs.unc.edu

Jean Feydy
email: jfeydy@ic.ac.uk

Peirong Liu
email: peirong@cs.unc.edu

Ariel Hernán Curiale
email: acuriale@bwh.harvard.edu

Ruben San
email: rubensanjose@bwh.harvard.edu

José Estépar

Raúl San

Marc Niethammer

Accurate Point Cloud Registration with Robust Optimal Transport

This work investigates the use of robust optimal transport (OT) for shape matching. Specifically, we show that recent OT solvers improve both optimization-based and deep learning methods for point cloud registration, boosting accuracy at an affordable computational cost. This manuscript starts with a practical overview of modern OT theory. We then provide solutions to the main difficulties in using this framework for shape matching. Finally, we showcase the performance of transport-enhanced registration models on a wide range of challenging tasks: rigid registration for partial shapes; scene flow estimation on the Kitti dataset; and nonparametric registration of lung vascular trees between inspiration and expiration. Our OT-based methods achieve state-of-the-art results on Kitti and for the challenging lung registration task, both in terms of accuracy and scalability. We also release PVT1010, a new public dataset of 1,010 pairs of lung vascular trees with densely sampled points. This dataset provides a challenging use case for point cloud registration algorithms with highly complex shapes and deformations. Our work demonstrates that robust OT enables fast pre-alignment and fine-tuning for a wide range of registration models, thereby providing a new key method for the computer vision toolbox. Our code and dataset are available online at: https://github.com/uncbiag/robot.

Introduction

Shape registration is a fundamental but difficult problem in computer vision. The task is to determine plausible spatial correspondences between pairs of shapes, with use cases that range from pose estimation for noisy point clouds [START_REF] Bouaziz | Sparse iterative closest point[END_REF] to the nonparametric registration of high-resolution medical images [START_REF] Brunn | Fast GPU 3D diffeomorphic image registration[END_REF]. As illustrated in Fig. 1, most existing approaches to this problem consist of a combination of three steps, possibly fused together by some deep learning (DL) methods: [START_REF] Ahrens | Paraview: An end-user tool for large data visualization[END_REF] feature extraction;

(2) feature matching; and (3) regularization using a class of acceptable transformations that is specified through a parametric or nonparametric model. This work discusses how tools derived from 1. Feature extraction. To establish spatial correspondences, one first computes descriptive local features. When dealing with (possibly annotated) point clouds, a simple choice is to rely on Cartesian coordinates (x, y, z) [START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF][START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF]. Going further, stronger descriptors capture local geometric and topological properties: examples include shape orientation and curvatures [START_REF] Charon | The varifold representation of nonoriented shapes for diffeomorphic registration[END_REF][START_REF] Roussillon | Kernel metrics on normal cycles and application to curve matching[END_REF], shape contexts [START_REF] Belongie | Shape context: A new descriptor for shape matching and object recognition[END_REF], spectral eigenvalues [START_REF] Lombaert | Spectral log-demons: diffeomorphic image registration with very large deformations[END_REF][START_REF] Ovsjanikov | Functional maps: a flexible representation of maps between shapes[END_REF] and annotations such as color [START_REF] Men | Color point cloud registration with 4D ICP algorithm[END_REF] or chemical fingerprints [START_REF] Gainza | Deciphering interaction fingerprints from protein molecular surfaces using geometric deep learning[END_REF][START_REF] Sverrisson | Fast end-to-end learning on protein surfaces[END_REF]. Recently, expressive feature representations have also been learned using deep neural networks (DNN): see [START_REF] Bronstein | Geometric deep learning: going beyond Euclidean data[END_REF] and subsequent works on geometric deep learning. Generally, feature extractors are designed to make shape registration as unambiguous as possible. In order to get closer to the ideal case of landmark matching [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF], we associate discriminative features to the salient points of our shapes: this increases the robustness of the subsequent matching and regularization steps.

2. Feature matching. Once computed on both of the source and target shapes, feature vectors are put in correspondence with each other. This assignment is often encoded as an explicit mapping between the two shapes; alternatively, the vector field relating the shapes can be defined implicitly as the gradient of a geometric loss function that quantifies discrepancies between two distributions of features [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF]: a) A first major approach is to rely on nearest neighbor projections and the related chamfer [START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF] and Hausdorff distances [START_REF] Bouaziz | Modern techniques and applications for real-time non-rigid registration[END_REF], as in the Iterative Closest Point (ICP) algorithm [START_REF] Besl | Method for registration of 3-D shapes[END_REF]. This method can be softened through the use of a softmax (log-sum-exp) operator as in the many variants of the Coherent Point Drift (CPD) method [START_REF] Gao | Filterreg: Robust and efficient probabilistic point-set registration using Gaussian filter and twist parameterization[END_REF][START_REF] Ma | Non-rigid visible and infrared face registration via regularized Gaussian fields criterion[END_REF][START_REF] Ma | Non-rigid point set registration by preserving global and local structures[END_REF][START_REF] Myronenko | Point set registration: Coherent point drift[END_REF], or made robust to outliers in the specific context of rigid and affine registrations [START_REF] Bouaziz | Sparse iterative closest point[END_REF][START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF][START_REF] Yang | A polynomial-time solution for robust registration with extreme outlier rates[END_REF][START_REF] Yang | Teaser: Fast and certifiable point cloud registration[END_REF]. b) Alternatively, a second approach is to rely on convolutional kernel norms such as the Energy Distance [START_REF] Rizzo | Energy distance[END_REF], which are also known as Maximum Mean Discrepancies (MMD) in statistics [START_REF] Gretton | A kernel two-sample test[END_REF]. These loss functions are common in imaging science [START_REF] Pluim | Mutual information matching in multiresolution contexts[END_REF] and computational anatomy [START_REF] Charon | The varifold representation of nonoriented shapes for diffeomorphic registration[END_REF][START_REF] Vaillant | Surface matching via currents[END_REF] but are prone to vanishing gradients [START_REF] Feydy | Interpolating between optimal transport and MMD using Sinkhorn divergences[END_REF][START_REF] Feydy | Global divergences between measures: from Hausdorff distance to optimal transport[END_REF]. c) Finally, a third type of approach is to rely on optimal transport (OT) theory [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF] and solutions of the earth mover's problem [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF]. This method is equivalent to a nearest neighbor projection under a global constraint of bijectivity that enforces consistency in the matching. On the one hand, OT has been known to provide reliable correspondences in computer vision for more than two decades [START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF][START_REF] Gold | New algorithms for 2D and 3D point matching: pose estimation and correspondence[END_REF][START_REF] Kosowsky | The invisible hand algorithm: Solving the assignment problem with statistical physics[END_REF]. On the other hand, it has often faced major issues of scalability and robustness to outliers on noisy data. As detailed below, the main purpose of this work is to overcome these limitations and enable the widespread use of OT tools for challenging registration problems.

3. Regularization with a deformation model. The output of the two steps above is a non-smooth vector field that may not be suitable for downstream tasks due to e.g. tears and compression artifacts. As a third step, most registration methods thus rely on regularization to obtain plausible deformations. This process is task-specific, with applications that range from rigid registration [2,[START_REF] Deng | Ppfnet: Global context aware local features for robust 3D point matching[END_REF][START_REF] Gojcic | The perfect match: 3D point cloud matching with smoothed densities[END_REF][START_REF] Wang | Deep closest point: Learning representations for point cloud registration[END_REF][START_REF] Wang | PRNet: Self-supervised learning for partial-to-partial registration[END_REF][START_REF] Zeng | 3DMatch: Learning local geometric descriptors from RGB-d reconstructions[END_REF] to free-form motion estimation [START_REF] Gu | HPLFlowNet: Hierarchical permutohedral lattice FlowNet for scene flow estimation on large-scale point clouds[END_REF][START_REF] Liu | FlowNet3D: Learning scene flow in 3D point clouds[END_REF][START_REF] Puy | FLOT: Scene flow on point clouds guided by optimal transport[END_REF][START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]. In Sec. 3, we address the interaction of OT matching layers with a varied collection of regularization strategies -from optimization-based spline and diffeomorphic models to DNNs.

Recent progresses. Research works on shape registration combine ideas from the three paragraphs above to best fit the characteristics of computer vision problems [START_REF] Donati | Deep geometric functional maps: Robust feature learning for shape correspondence[END_REF][START_REF] Lüthi | Gaussian process morphable models[END_REF][START_REF] Shen | Networks for joint affine and non-parametric image registration[END_REF]. Over the past few years, significant progress has been made on all fronts. On the one hand, (geometric) deep learning networks have been used to define data-driven feature maps [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF][START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF][START_REF] Wang | Dynamic graph CNN for learning on point clouds[END_REF] and multiscale regularization modules [START_REF] Litany | Deep functional maps: Structured prediction for dense shape correspondence[END_REF][START_REF] Shen | Region-specific diffeomorphic metric mapping[END_REF][START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF], sometimes fused within end-to-end architectures [START_REF] Deng | Ppfnet: Global context aware local features for robust 3D point matching[END_REF][START_REF] Puy | FLOT: Scene flow on point clouds guided by optimal transport[END_REF][START_REF] Yew | RPM-Net: Robust point matching using learned features[END_REF][START_REF] Zeng | 3DMatch: Learning local geometric descriptors from RGB-d reconstructions[END_REF]. On the other hand, nearest neighbor projections, kernel convolutions and transport-based matching strategies have all been generalized to take advantage of these modern descriptors: they can now be used in high-dimensional feature spaces [START_REF] Feydy | Fast geometric learning with symbolic matrices[END_REF][START_REF] Johnson | Billion-scale similarity search with GPUs[END_REF].

Challenges. Nevertheless, state-of-the-art (SOTA) methods in the field still have important limitations. First, modern deep learning pipelines are often hard to train to "pixel-perfect" accuracy on nonsmooth shapes, with diminishing returns in terms of model size and training data [2]. Second, scaling up point neural networks to finely sampled shapes (N > 10k points) remains a challenging research topic [START_REF] Feydy | Fast geometric learning with symbolic matrices[END_REF][START_REF] Gu | HPLFlowNet: Hierarchical permutohedral lattice FlowNet for scene flow estimation on large-scale point clouds[END_REF][START_REF] Zhou | Fully convolutional mesh autoencoder using efficient spatially varying kernels[END_REF]. Third, the impact of the choice of a specific feature matching method on the performance of deep learning models remains only partially understood [START_REF] Huang | A comprehensive survey on point cloud registration[END_REF]. (a) RobOT is equivalent to a nearest neighbor projection subject to mass distribution constraints that make it robust to translations and small deformation. We demonstrate that RobOT is now ready to be part of the standard toolbox in computer vision with extensive numerical experiments for 3D scene flow estimation (b) and lung registration (c). Rendering done with Paraview [START_REF] Ahrens | Paraview: An end-user tool for large data visualization[END_REF] and PyVista [START_REF] Sullivan | PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (vtk)[END_REF].

Related works. Following major progress on computational OT in the mathematical literature [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for L2 semi-discrete optimal transport in 3D[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] Schmitzer | Stabilized sparse scaling algorithms for entropy regularized transport problems[END_REF], improved modules for feature matching have attracted interest as a possible solution to these challenges. Works on sliced partial OT [START_REF] Bonneel | SPOT: sliced partial optimal transport[END_REF] and dustbin-OT [START_REF] Dang | Learning 3D-3D correspondences for one-shot partial-to-partial registration[END_REF] have shown that outliers can be handled effectively by OT methods for rigid registration, beyond the classic Robust Point Matching method (RPM) [START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF][START_REF] Gold | New algorithms for 2D and 3D point matching: pose estimation and correspondence[END_REF]. Going further, the Sinkhorn algorithm for entropy-regularized OT [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Kosowsky | The invisible hand algorithm: Solving the assignment problem with statistical physics[END_REF][START_REF] Léonard | From the Schrödinger problem to the Monge-Kantorovich problem[END_REF][START_REF] Schrödinger | Sur la théorie relativiste de l'électron et l'interprétation de la mécanique quantique[END_REF] has been studied extensively for shape registration in computational anatomy [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF][START_REF] Gerber | Exploratory population analysis with unbalanced optimal transport[END_REF] and computer graphics [START_REF] Eisenberger | Deep shells: Unsupervised shape correspondence with optimal transport[END_REF][START_REF] Mandad | Variance-minimizing transport plans for inter-surface mapping[END_REF][START_REF] Pai | Fast Sinkhorn filters: Using matrix scaling for non-rigid shape correspondence with functional maps[END_REF]. The Gromov-Wasserstein distance has also been used for shape analysis [START_REF] Solomon | Entropic metric alignment for correspondence problems[END_REF][START_REF] Vayer | Sliced Gromov-Wasserstein[END_REF], albeit at a higher computational cost. These applications have driven interest in the development of a complete theory for Robust Optimal Transport (RobOT), outlined in Sec. 2, which handles both sampling artifacts and outliers [START_REF] Chizat | Scaling algorithms for unbalanced optimal transport problems[END_REF][START_REF] Chizat | Unbalanced optimal transport: Dynamic and Kantorovich formulations[END_REF][START_REF] Liero | Optimal entropy-transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF][START_REF] Mukherjee | Outlier-robust optimal transport[END_REF][START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF][START_REF] Séjourné | The unbalanced Gromov Wasserstein distance: Conic formulation and relaxation[END_REF]. Most recently, this framework has started to be used in shape analysis with applications to shape matching [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF], the segmentation of brain tractograms [START_REF] Feydy | Fast and scalable optimal transport for brain tractograms[END_REF] and deep deformation estimation with the FLOT architecture [START_REF] Puy | FLOT: Scene flow on point clouds guided by optimal transport[END_REF].

Contributions. We build upon the work above to tackle challenging point cloud registration problems for scene flow estimation and computational anatomy. Our key contributions are:

1. Accurate feature matching with scalable OT solvers. For the first time, we scale up RobOT for deep feature matching to high-resolution shapes with more than 10k points. To this end, we leverage the latest generation of OT solvers [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF][START_REF] Feydy | Interpolating between optimal transport and MMD using Sinkhorn divergences[END_REF] and overcome significant issues of memory usage and numerical stability. This allows us to handle fine-grained details effectively, which is key for e.g. most medical applications. 2. Interaction with task-specific regularization strategies. We show how to interface RobOT matchings with advanced deformation models. This is in contrast with e.g. the FLOT architecture, which focuses on the direct prediction of a vector field and cannot be used for applications that require guarantees on the smoothness of the registration. 3. Challenging new dataset. We release a large dataset of lung vascular trees that should be registered between inhalation and exhalation. This relevant medical problem involves large and complex deformations of high-resolution 3D point clouds. As a new benchmark for the community, we provide two strong baselines that rely respectively on global feature matching and on deep deformation estimation. 4. Consistent SOTA performance. Our proposed models achieve SOTA results for scene flow on Kitti [START_REF] Menze | Object scene flow for autonomous vehicles[END_REF][START_REF] Menze | Joint 3D estimation of vehicles and scene flow[END_REF] and for point-cloud-based lung registration on DirLab-COPDGene [START_REF] Castillo | A reference dataset for deformable image registration spatial accuracy evaluation using the COPDgene study archive[END_REF]. Notably, we show that RobOT is highly suited to fine-tuning tasks: it consistently turns "good" matchings into nearly perfect registrations at an affordable numerical cost.

Main experimental observations. Is OT relevant in the deep learning era? To answer this question decisively, we perform extensive numerical experiments and ablation studies. We fully document "underwhelming" results in the Supplementary Material and distill the key lessons that we learned in the Deep-RobOT architecture (Section 3.2). This model relies on fast RobOT layers to cover for the main weaknesses of point neural networks for shape registration. It is remarkably easy to deploy and generalizes well from synthetic training data to real test samples. We thus believe that it will have a stimulating impact on both of the computer vision and medical imaging literature.

Robust optimal transport

This section introduces the mathematical foundations of our work. After a brief overview of Robust Optimal Transport (RobOT) theory, we discuss the main challenges that one encounters when using this framework for computer vision. To avoid memory overflows and numerical stability issues, we introduce the weighted "RobOT matching": a vector field that summarizes the information of a full transport plan with a linear memory footprint. As detailed in the next sections, this representation lets us scale up to high-resolution shapes without compromising on accuracy.

Mathematical background

The assignment problem. If A = (x 1 , . . . , x N ) and B = (y 1 , . . . , y M ) are two point clouds in R 3 with N = M, the assignment problem between A and B reads:

Assignment(A, B) = min s:[[1,N]]→[[1,N]] 1 2N N i=1 x i -y s(i) 2 R 3
, where s is a permutation. [START_REF] Ahrens | Paraview: An end-user tool for large data visualization[END_REF] This problem generalizes sorting to R 3 : if the points x i and y j all belong to a line, the optimal permutation s * corresponds to a non-decreasing re-ordering of the point sets A and B [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF].

Robust optimal transport. Further, OT theory allows us to consider problems where N = M. Non-negative weights α 1 , . . . , α N , β 1 , . . . , β M 0 are attached to the points x i , y j and account for variations of the sampling densities, while feature vectors p 1 , . . . , p N and q 1 , . . . , q M in R D may advantageously replace raw point coordinates x i and y j in R 3 . Following [START_REF] Chizat | Unbalanced optimal transport: Dynamic and Kantorovich formulations[END_REF][START_REF] Liero | Optimal entropy-transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF], the robust OT problem between the shapes A = (α i , x i , p i ) and B = (β j , y j , q j ) reads:

OT σ,τ (A, B) = min (πi,j )∈R N×M 0 N i=1 M j=1 π i,j • 1 2 p i -q j 2 R D (2)
+ σ 2 KL π i,j || α i ⊗ β j
Entropic blur at scale σ.

+ τ 2 KL j π i,j || α i π should match A. . . + τ 2 KL i π i,j || β j . . . onto B.
, for any choice of the regularization parameters σ > 0 and τ > 0. In the equation above, the Kullback-Leibler divergence KL (a i ||b i ) = a i log (a i /b i ) -a i + b i is a relative entropy that penalizes deviations of a non-negative vector of weights (a i ) to a reference measure (b i ).

Parameters. The first regularization term is scaled by the square of a blur radius σ. This characteristic length quantifies the fuzziness of the probabilistic assignment (π i,j ) between points x i and y j [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF]. The last two regularization terms promote the matching of the full distribution of points A onto the target shape B: they generalize the constraints of injectivity and surjectivity of Eq. (1) to the probabilistic setting. They are scaled by the square of a maximum reach distance τ : this parameter acts as a soft upper bound on the distance between feature vectors p i and q j that should be matched with each other [START_REF] Feydy | Fast and scalable optimal transport for brain tractograms[END_REF][START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF].

For shape registration, we use simple heuristics for the values of these two characteristic scales: the blur σ should be equal to the average sampling distance in feature space R D while the reach τ should be equal to the largest plausible displacement for any given feature vector p i . These rules are easy to follow if point features correspond to Cartesian coordinates x i and y j in R 3 but may lead to unexpected behaviors if features are output by a DNN. In the latter case, we thus normalize our feature vectors so that p i R D = q j R D = 1 and pick values for σ and τ between 0 and 2.

Working with a soft, probabilistic transport plan. As detailed in [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF], scalable OT solvers for Eq. ( 2) return a pair of dual vectors (f i ) ∈ R N and (g j ) ∈ R M that encode implicitly an optimal transport plan (π i,j ) ∈ R N×M with coefficients:

π i,j = α i β j • exp 1 σ 2 f i + g j -1 2 p i -q j 2 R D 0 . (3)
In the limit case where σ tends to 0 and τ tends to +∞, for generic point clouds (x i ) and (y j ) with N = M and equal weights α i = β j = 1/N, π i,j is a permutation matrix [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF]. We retrieve the simple assignment problem of Eq. ( 1): π i,j = 1/N if j = s * (i) and 0 otherwise. However, in general the transport plan must be understood as a probabilistic map between the point distributions A and B that assigns a weight π i,j to the coupling "x i ↔ y j ". For shape registration, this implies that the main difficulties for using robust OT are two-fold: first, the coupling π is not one-to-one, but one-to-many; second, the lines and columns of the transport plan π do not sum up to one. Notably, this implies that when τ < +∞, the gradient of the OT cost with respect to the point positions x i is not homogeneous: we observe vanishing and inflated values across the domain [START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF].

RobOT: a convenient representation of the optimal transport plan

The weighted RobOT matching. To work around these issues, we introduce the vector field:

v i = M j=1 π i,j • (y j -x i ) M j=1 π i,j ∈ R 3 with confidence weights w i = M j=1 π i,j 0 . ( 4
)
This object has the same memory footprint as the input shape A and summarizes the information that is contained in the N-by-M transport plan (π i,j ) -a matrix that is often too large to be stored and manipulated efficiently. This "weighted RobOT matching" is at the heart of our approach and generalizes the standard Monge map from classical OT theory [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF] to the setting of (deep) shape registration. In practice, the weighted vector field (w 1 , v 1 ), . . . , (w N , v N ) is both convenient to use and easy to compute on GPUs. Let us briefly explain why.

Fast implementation. Our differentiable RobOT layer takes as input the two shapes A = (α i , x i , p i ) and B = (β j , y j , q j ), with feature vectors p i and q j in R D that have been computed upstream using e.g. a point neural network. It returns the N vectors v i with weights w i that map the source points x i onto the targets y j in R 3 . Starting from the input point features p i , q j and weights α i , β j , we first compute the optimal dual vectors f i and g j using the fast solvers of the GeomLoss library [START_REF] Feydy | Interpolating between optimal transport and MMD using Sinkhorn divergences[END_REF]. We then combine Eq. (3) with Eq. ( 4) to compute the RobOT vectors v i and weights w i with O(N + M) memory footprint using the KeOps library [START_REF] Charlier | Kernel operations on the GPU, with autodiff, without memory overflows[END_REF][START_REF] Feydy | Fast geometric learning with symbolic matrices[END_REF] for PyTorch [START_REF] Paszke | Automatic differentiation in PyTorch[END_REF] and NumPy [START_REF] Van Der Walt | The NumPy array: a structure for efficient numerical computation[END_REF].

We use a log-sum-exp formulation to ensure numerical stability. Remarkably, our implementation scales up to N, M = 100k in fractions of a second. Unlike common strategies that are based on dense or sparse representations of the transport plan π, our approach is perfectly suited to a symbolic implementation [START_REF] Feydy | Fast geometric learning with symbolic matrices[END_REF] and streams well on GPUs with optimal, contiguous memory accesses.

Comparison with nearest neighbor projections. We use our RobOT layer as a plug-in replacement for closest point matching [START_REF] Besl | Method for registration of 3-D shapes[END_REF]. The blur (σ) and reach (τ ) scales play similar roles to the standard deviation (σ) and weight of the uniform distribution (w) in the Coherent Point Drift (CPD) method [START_REF] Myronenko | Point set registration: Coherent point drift[END_REF]: they allow us to smooth the matching in order to increase robustness to sampling artifacts.

The main difference between projection-based matching and RobOT is that the latter enforces a mass distribution constraint between the source and the target. This prevents our matching vectors from accumulating on the boundaries of the distributions of point features (p 1 , . . . , p N ) and (q 1 , . . . , q M ) in R D [START_REF] Feydy | Global divergences between measures: from Hausdorff distance to optimal transport[END_REF]. This property is most desirable when the shapes to register are fully observed, with a dense sampling: as detailed in Suppl. A.5, enforcing the global consistency of a matching is then a worthwhile registration prior.

Partial registration. On the other hand, we must also stress that OT theory has known limitations [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF]. First of all, the RobOT matching cannot guarantee the preservation of remarkable points or of the shapes' topologies "on its own": it should be combined with relevant feature extractors and regularizers. Going further, partial registration is not a natural fit for standard OT formulations which assume that all points from both shapes must be put in correspondence with each other.

To mitigate this issue, RobOT leverages the theory of unbalanced optimal transport [START_REF] Chizat | Scaling algorithms for unbalanced optimal transport problems[END_REF][START_REF] Chizat | Unbalanced optimal transport: Dynamic and Kantorovich formulations[END_REF][START_REF] Liero | Optimal entropy-transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF][START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF][START_REF] Séjourné | The unbalanced Gromov Wasserstein distance: Conic formulation and relaxation[END_REF]: we rely on soft Kullback-Leibler penalties to enforce a matching between the shapes A and B in Eq. (2). In practice, the RobOT confidence weights w i of Eq. ( 4) act as an attention mechanism: they vanish when no target feature vector q j can be found in a τ -neighborhood of the source vector p i in R D , where τ is the reach scale that is associated to Eq. ( 2). This lets our registration method focus on reliable matches between similar features, without being fooled by strong constraints of bijectivity.

As detailed in Suppl. A.3, combining standard FPFH features [START_REF] Rusu | Fast point feature histograms (FPFH) for 3D registration[END_REF] with the rigid projection of Eq. ( 5) allows us to register partially observed shapes that have little overlap with each other.

3 Regularization and integration with a deep learning model 3.1 Smooth-RobOT: algebraic and optimization-based regularization Notations. We now detail how to interface the weighted RobOT matching with regularization models and feature extractors that may be handcrafted [START_REF] Rusu | Fast point feature histograms (FPFH) for 3D registration[END_REF][START_REF] Rusu | Aligning point cloud views using persistent feature histograms[END_REF][START_REF] Tombari | Unique signatures of histograms for local surface description[END_REF] or learnt using a deep neural network [START_REF] Deng | Ppf-foldnet: Unsupervised learning of rotation invariant 3D local descriptors[END_REF][START_REF] Deng | Ppfnet: Global context aware local features for robust 3D point matching[END_REF][START_REF] Gojcic | The perfect match: 3D point cloud matching with smoothed densities[END_REF][START_REF] Zeng | 3DMatch: Learning local geometric descriptors from RGB-d reconstructions[END_REF]. Recall that we intend to register a source point cloud x 1 , . . . , x N onto a target y 1 , . . . , y M in R 3 . Non-negative weights α 1 , . . . , α N and β 1 , . . . β M 0 let us take into account variations in the sampling densities and we assume that point features p 1 , . . . , p N , q 1 , . . . , q M in R D have been computed upstream by a relevant feature extractor. For every source point x i , the RobOT matching layer then provides a desired displacement v i in R 3 with influence weight w i 0.

Smoothing in closed form. Standard computations let us derive closed-form expressions for rigid and affine registration [START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF][START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF][START_REF] Myronenko | Point set registration: Coherent point drift[END_REF]. These respectively correspond to transformations:

(Rigid RobOT) x ∈ R 3 → (x -x c ) U V + x c + v c ∈ R 3 , (5)
(Affine RobOT) x ∈ R 3 → (x -x c ) ( X W X) -1 ( X W Ŷ ) + x c + v c ∈ R 3 , (6)
where W = Diag(w i ) ∈ R N×N is the diagonal matrix of influence weights,

x c = i w i x i / i w i ∈ R 3 is the barycenter of the source shape, v c = i w i v i / i w i ∈ R 3 is the average desired displacement, X = (x i -x c ) ∈ R N×3 is the centered matrix of source positions, Ŷ = (x i + v i -x c -v c ) ∈ R N×3
is the centered matrix of desired targets and U SV is the singular value decomposition of X W Ŷ ∈ R 3×3 . This corresponds to a weighted Kabsch algorithm [START_REF] Kabsch | A solution for the best rotation to relate two sets of vectors[END_REF]. Likewise, we implement free-form spline registration using the KeOps library [START_REF] Charlier | Kernel operations on the GPU, with autodiff, without memory overflows[END_REF][START_REF] Feydy | Fast geometric learning with symbolic matrices[END_REF]. A Nadaraya-Watson interpolator with kernel k : (x, y) ∈ R 3 × R 3 → k(x, y) > 0 [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF] induces a transformation:

(Spline RobOT) x ∈ R 3 → x + N i=1 w i k(x i , x)v i / N i=1 w i k(x i , x) ∈ R 3 . (7)
Black-box deformation models. Going further, we interface RobOT matchings with arbitrary deformation modules Morph : (θ, x i ) → ŷi ∈ R 3 parameterized by a vector θ in R P . If Reg(θ) denotes a regularization penalty on the parameter θ (e.g. a squared Euclidean norm), we use standard optimizers such as L-BFGS-B [START_REF] Zhu | Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[END_REF] and Adam [START_REF] Kingma | A method for stochastic optimization[END_REF] to find the optimal deformation parameter:

θ * = arg min θ∈R P Reg(θ) + N i=1 w i x i + v i -Morph(θ, x i ) 2 R 3 . (8)
This optimization-based method is especially relevant in the context of computational anatomy, where smooth and invertible deformations are commonly defined through the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Bône | Deformetrica 4: an open-source software for statistical shape analysis[END_REF][START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF]. We stress that different transformation models may result in different registration results and refer to Suppl. A.3 for further details.

Optimization-based approaches provide strong geometric guarantees on the final matching. But unfortunately, these often come at a high computational price: to register complex shapes, quasi-Newton optimizers require dozens of evaluations of the deformation model Morph(θ, x) and of its gradients. In practice, fitting a complex model to a pair of high-resolution shapes may thus take several minutes or seconds [START_REF] Brunn | Fast GPU 3D diffeomorphic image registration[END_REF]. This precludes real-time processing and hinders research on advanced deformation models.

Deep-RobOT: registration via deep deformation prediction

Deformation prediction. In this context, there is growing interest in fast learning methods that avoid the use of iterative optimizers. The idea is to train a deep neural network Pred : (x i , y j ) → θ that takes as input two point clouds A = (x 1 , . . . , x N ), B = (y 1 , . . . , y M ) and directly predicts the optimal vector of parameters θ for a transformation model Morph(θ, • ) that should map A onto B. Assuming that the prediction network Pred has been trained properly, this strategy enables real-time processing while leveraging the task-specific geometric priors that are encoded within the deformation model.

Over the last few years, numerous authors have worked in this direction for rigid registration [2,[START_REF] Deng | Ppfnet: Global context aware local features for robust 3D point matching[END_REF][START_REF] Gojcic | The perfect match: 3D point cloud matching with smoothed densities[END_REF][START_REF] Wang | Deep closest point: Learning representations for point cloud registration[END_REF][START_REF] Wang | PRNet: Self-supervised learning for partial-to-partial registration[END_REF][START_REF] Zeng | 3DMatch: Learning local geometric descriptors from RGB-d reconstructions[END_REF] and scene flow estimation [START_REF] Gu | HPLFlowNet: Hierarchical permutohedral lattice FlowNet for scene flow estimation on large-scale point clouds[END_REF][START_REF] Liu | FlowNet3D: Learning scene flow in 3D point clouds[END_REF][START_REF] Puy | FLOT: Scene flow on point clouds guided by optimal transport[END_REF][START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]. Comparable research on diffeomorphic models has focused on images that are supported on a dense grid, with successful applications to e.g. the registration of 3D brain volumes [START_REF] Balakrishnan | Voxelmorph: a learning framework for deformable medical image registration[END_REF][START_REF] Shen | Networks for joint affine and non-parametric image registration[END_REF][START_REF] Yang | QuickSilver: Fast predictive image registration-a deep learning approach[END_REF]. As of 2021, prediction-based approaches have thus become standard methods for 3D shape registration.

1. Affine-RobOT pre-alignment. The D-RobOT model. In practice though, prediction methods still face three major challenges:

1. Common architectures may not be equivariant to rigid or affine transformations.

2. Dense 3D annotation is expensive, especially in a medical context. As a consequence, most predictors are trained on synthetic data and have to overcome a sizeable domain gap. We propose to address these issues using RobOT layers for pre-alignment and post-processing. Our Deep RobOT (D-RobOT) model is an end-to-end deep learning architecture that is made up of three consecutive steps that we illustrate in Figure 2 and detail in Suppl. A.4:

1. OT-based pre-alignment. We use the rigid or affine S-RobOT models of Eqs. [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Belongie | Shape context: A new descriptor for shape matching and object recognition[END_REF] to normalize the pose of the source shape. This is a fast and differentiable pre-processing. 2. Deep registration module. We combine a deep predictor with a task-specific deformation model to register the pre-aligned source onto the target. For the prediction network Pred : (x i , y j ) → θ, we use a multiscale point neural network that is adapted from the PointPWC-Net architecture [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]. We refer to Suppl. A.4 for a full description of our architecture and training loss. 3. OT-based post-processing. In order to reach "pixel-perfect" accuracy, we use the spline S-RobOT deformation model of Eq. ( 7) with a task-specific kernel k(x, y).

Complementary strengths and weaknesses. We apply these three steps successively, which brings the moving source A = (x 1 , . . . , x N ) increasingly close to the fixed target B = (y 1 , . . . , y M ). Remarkably, each step of our method covers for the weaknesses of the other modules: the RobOT-based pre-alignment makes our pipeline robust to changes of the 3D acquisition parameters; our multiscale neural predictor is able to match corresponding key points quickly, even in complex situations; the domain-specific deformation model acts as a regularizer and improves the generalization properties of the deep registration module; the RobOT-based fine-tuning improves accuracy and helps our model to overcome the domain gap between synthetic and real shape data.

As detailed below, the D-RobOT model generalizes well outside of its training dataset and outperforms state-of-the-art methods on several challenging problems. We see it as a pragmatic architecture for shape registration, which is easy to deploy and tailor to domain-specific requirements. As discussed in Suppl. A.4, we found that introducing sensible geometric priors through our RobOT layers and the deformation model Morph : (θ, x i ) → ŷi results in a "forgiving" pipeline: our model produces accurate results, even when trained on synthetic data that is not very realistic. In a context where generating plausible 3D deformations is easier than developing custom registration models for every single task (e.g. in computational anatomy), we believe that this is an important observation.

Scene flow estimation

Benchmark. We now evaluate our method on a standard registration task in computer vision: the estimation of scene flow between two successive views of the same 3D scene. We follow the same experimental setting as in [START_REF] Gu | HPLFlowNet: Hierarchical permutohedral lattice FlowNet for scene flow estimation on large-scale point clouds[END_REF][START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF], with full details provided in Suppl. A.4.3:

1. We train on the synthetic Flying3D dataset [START_REF] Mayer | A large dataset to train convolutional networks for disparity, optical flow, and scene flow estimation[END_REF], which is made up of multiple moving objects that are sampled at random from ShapeNet. We take 19,640 pairs of point clouds for training, with dense ground truth correspondences. 2. We evaluate on 142 scene pairs from Kitti, a real-world dataset [START_REF] Menze | Object scene flow for autonomous vehicles[END_REF][START_REF] Menze | Joint 3D estimation of vehicles and scene flow[END_REF]. We conduct experiments using 8,192 and 30k points per scan, sampled at random from the original data.

Performance metrics. We evaluate all methods as in [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]: EPE3D is the average 3D error, in centimeters; Acc3DS is the percentage of points with 3D error < 5 cm or relative error < 5%; Acc3DR is the percentage of points with 3D error < 10 cm or relative error < 10%; Outliers3D is the percentage of points with 3D error > 30 cm or relative error > 10%; EPE2D is the average 2D error obtained by projecting the point clouds onto the image plane, measured in pixels; Acc2D is the percentage of points with 2D error < 3 px or relative error < 5%. As detailed in Suppl. A.6, all run times were measured on a single GPU (24GB NVIDIA Quadro RTX 6000).

Methods. We study a wide range of methods and report the relevant metrics in Fig. 4:

In the upper third of the table, we report results for unsupervised methods that do not require ground truth correspondences for training. This includes the "raw" RobOT matching of Eq. ( 4), computed on (x, y, z) coordinates in R 3 with a blur scale σ = 1 cm and a reach scale τ = +∞. Please also note that PWC refers to an improved version of PointPWC-Net, released on GitHub (https: //github.com/DylanWusee/PointPWC) after the publication of [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF] with a self-supervised loss.

In the central third of the table, we benchmark a collection of state-of-the-art point neural networks.

In the lower third of the table, we study the influence of our RobOT-based layers. The methods "Pre + FLOT/PWC + Post" correspond to the FLOT and PointPWC-Net architectures, with the additional pre-alignment and post-processing modules of Sec. 3.2. The last two lines correspond to the full D-RobOT architecture (with a spline deformation model) whose training is detailed in Suppl. A.4.3.

Results. We make three major observations:

1. Without any regularization or training, a simple RobOT matching on high-resolution data outperforms many deep learning methods in terms of speed, memory footprint and accuracy (line 6 of the table). This surprising result is strong evidence that geometric methods and baselines deserve more attention from the computer vision community. 2. In the lower third of the table, RobOT-enhanced methods consistently outperform stateof-the-art methods by a wide margin. 3. As shown in Fig. 3, these improvements are most significant on high-resolution data.

Overall, as detailed in Suppl. A.4.3 and A.5, we observe that optimal transport theory is especially well suited to scene flow estimation. Assuming that ground points have been removed from the 3D frames (as a standard pre-processing), most object displacements can be explained as translations and small rotations: this is an ideal setting for our robust geometric method.

5 Registration of high-resolution lung vessel trees PVT1010: a new dataset for lung registration. Going further, we introduce a new dataset of 1,010 pairs of pulmonary vessel trees that must be registered between expiration (source) and inspiration (target). Due to the intricate geometry of the lung vasculature and the complexity of the breathing motion, the registration of these shapes is a real challenge.

As detailed in Suppl. A.1, we encode our 1,010 × 2 vessel trees as high-resolution 3D point clouds (N = M = 60k points per tree). For each point, we also provide a local estimate of the vessel radius that we use as an additional point feature or as a weight α i or β j in Eq. (2). Our first 1,000 pairs of 3D point clouds are provided without ground truth correspondences; for all of our experiments, we randomly sample 600 training and 100 validation cases from this large collection of unannotated patients. The last 10 cases correspond to the 10 DirLab COPDGene pairs [START_REF] Castillo | A reference dataset for deformable image registration spatial accuracy evaluation using the COPDgene study archive[END_REF]: they come with 300 expert-annotated 3D landmarks per lung pair, that we use to test our methods. Overall, as discussed in Suppl. A.4.2, we found that supervised training on synthetic deformations is both easier and more efficient than unsupervised training on real lung pairs. From the chamfer and Wasserstein distances [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF] to local Laplacian penalties [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF], none of the unsupervised loss functions that we experimented with was able to deal with the complex geometry of our lung vascular trees.

Methods. We benchmark a wide range of methods in Fig. 5 and Tab. 1. In the upper half of the table, we evaluate geometric approaches that require no training; in the lower half of the table, we benchmark scalable point neural networks that we trained on our synthetic dataset. We evaluate three types of RobOT-based approaches: a simple RobOT matching computed using (x, y, z) coordinates as point features, that may either be "raw" as in Eq. ( 4) or regularized using Eq. ( 6); an S-RobOT matching that we compute using the deep features of Suppl. A.3.3 and regularize with the spline smoothing of Eq. ( 7) or the LDDMM optimization of Eq. ( 10); a D-RobOT architecture that we pair with three deformation models Morph(θ, x i ) → ŷi and describe in depth in Suppl. A.4.2.

Results. We provide additional experiments in Suppl. A.4.2 and make three major observations:

1. Explicit regularization with a spline or LDDMM deformation model is key. Model-free architectures that predict raw 3D correspondences produce non-smooth results that are not anatomically plausible, even when they are trained entirely on smooth deformations. 2. The D-RobOT architecture combines a high acccuracy with fast run times. 3. Most remaining errors occur at the boundary of the lungs, where acquisition artifacts prevent the thinnest vessels from being sampled reliably in our point cloud representation.

Conclusion, limitations and future work

Our work builds upon a decade of active research in the field of computational optimal transport. We leverage major advances on RobOT solvers to define a new matching layer which is a plugand-play replacement for nearest neighbor projection. This operation has two major uses in 3D shape registration: first, it provides a very strong geometric baseline for e.g. scene flow estimation; second, it increases the accuracy and generalization abilities of point neural networks on finely sampled 3D shapes. We see D-RobOT as a mature and versatile architecture for shape registration which is easy to train and adapt to task-specific requirements in e.g. medical imaging.

Going forward, we see three main ways of improving this work. First, we still have to investigate in depth the important problem of occlusions and partial acquisitions. Second, integrating taskspecific features beyond the (x, y, z) point coordinates is often key to perfect results. In the specific setting of lung registration, working with image-based features or focusing on branching points could be a way of improving performance at the cost of portability: recent works such as [START_REF] Hansen | Deep learning based geometric registration for medical images: How accurate can we get without visual features[END_REF][START_REF] Hansen | GraphRegNet: Deep graph regularisation networks on sparse keypoints for dense registration of 3d lung CTs[END_REF] are an excellent source of inspiration. Finally, we believe that high-quality software packaging is an important part of research in our field. We intend to keep working on the topic and distribute our methods through a user-friendly Python library for widespread use by the scientific community.

A Supplementary Material

We provide additional information on our RobOT layer and the two derived methods:

• S-RobOT: feature matching with RobOT, followed by a smoothing that is implemented in closed form or relies on an optimization loop -as detailed in Sec. 3.1.

• D-RobOT: S-RobOT as a pre-and post-processing tool for a deep deformation estimatoras detailed in Sec. 3.2.

More specifically:

1. Sec. A.1 provides details on PVT1010, our dataset of pulmonary vascular trees. 2. Sec. A.2 describes our synthetic deformations for augmenting PVT1010. 3. Sec. A.3 contains more information on global feature matching with S-RobOT. 4. Sec. A.4 provides details on our deep deformation prediction approach D-RobOT, with additional results on the PVT1010 and Kitti datasets. 5. Sec. A.5 discusses the differences between RobOT and nearest neighbor projection. 6. Sec. A.6 details our computational resources. 7. Finally, Sec. A.7 discusses the societal impact of our work.

A.1 Pulmonary vascular tree dataset

Introducing the PVT1010 dataset. We introduce a new pulmonary vascular tree dataset for point cloud registration. The PVT1010 dataset includes 1,010 pairs of inhale/exhale lung vascular trees extracted from 3D computed tomography (CT) images; 10 of these correspond to the 10 cases of the public DirLab-COPDGene [START_REF] Castillo | A reference dataset for deformable image registration spatial accuracy evaluation using the COPDgene study archive[END_REF] dataset which includes, for each pair, 300 expert annotated landmarks that are in correspondence with each other and that we use to validate our results. We extracted the lung vascular geometry in both inspiratory and expiratory CT scans using a scale-space particle system [START_REF] Estépar | Computational vascular morphometry for the assessment of pulmonary vascular disease based on scale-space particles[END_REF][START_REF] Kindlmann | Sampling and visualizing creases with scale-space particles[END_REF] that is implemented in the Teem library [61]. We used the pipeline that is defined in the chest imaging platform [START_REF] Nardelli | Generative-based airway and vessel morphology quantification on chest CT images[END_REF][START_REF] San | Chest imaging platform[END_REF].

Legal and regulatory information. The vascular tree reconstructions that are used in this study were part of the COPDGene study (NCT00608764). This study has been IRB approved and participants have provided their consent. The investigators from the Brigham and Women's Hospital (Harvard Medical School) only had access to de-identified CT images to perform vascular reconstructions. Since we are performing secondary analysis using de-identified data, the work under consideration is not considered human subjects research and did not imply additional risks to the participants. Risks related to ionizing radiation exposure were described in the primary IRB-approved study. Study identifiers were re-coded for our release of PVT1010 to preserve anonymity. PVT1010 is released under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License.

Point cloud representation.

We now detail how we extracted the lung vascular trees as highresolution 3D point clouds from the raw CT images. To perform this geometric segmentation task, we rely on a system of 4D particles that are defined by three spatial coordinates plus one scale parameter that corresponds to the local radius of the lung vessel -these radii are used by RobOT as point weights α i and β j in Eq. (2). Our segmentation method starts from a point cloud that is initialized using the Frangi filter [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. Then, we fit this point cloud representation to our 3D CT volumes by minimizing iteratively a system energy that is expressed as the sum of:

• An inter-particle regularization energy that ensures convenient sampling properties. We rely on a sum of quartic polynomials of the pairwise point distances, with a tunable potential well that is chosen to induce regular sampling at a fixed distance between the points.

• A particle-image data fidelity term, which is computed using the image Hessian at the current particles' locations.

As a final result, we obtain points that are approximately equally distributed along the vessel centerlines. The vessel radii are first approximated as the image scales at which the middle eigenvalues of the Hessian are locally minimized, and then refined using the generative approach of [START_REF] Nardelli | Generative-based airway and vessel morphology quantification on chest CT images[END_REF].

The resulting dataset has the following properties, which result in a challenging registration task:

1. The vascular trees have a complex structure and exhibit complex, large motions between inhalation and exhalation. 2. To capture the complex anatomical structure of the lungs at millimeter scale, registration methods need to focus on branching points or rely on high-resolution point clouds. 3. Due to the fixed image resolution of the raw CT volume and to acquisition differences between the inhale and exhale scans, the extracted inhale and exhale vascular trees are not fully consistent with each other. This is especially true for tiny structures at the lung boundary, that may not be visible on the smaller lung images at exhalation time.

Sampling density, resolution. The original CT images from which we extract our point clouds are acquired with a uniform resolution on the x, y and z axes: depending on the patients, the side length of our voxels varies between 0.60 mm and 0.65 mm. Using the processing above, we turn these volumetric images into 3D point clouds with 60k samples per lung vascular tree: depending on the subject, the average sampling distance (from each point to its nearest neighbor in the 3D point cloud) ranges between 0.6 mm and 1.0 mm.

We note that for our 10 test cases, the Dirlab annotations were performed on a down-sampled volume with a resolution of 2.5 mm on the z axis. To guarantee a fair comparison between image-based and point-based methods, we follow standard practice for this dataset (https://www.dir-lab.com/ Results.html) and report our results in Tab. 1 with a "snap-to-voxel" post-processing: we quantize our 3D lung registrations on the original grids (with spacing ∼ 0.625 mm × 0.625 mm × 2.5 mm) before computing the average errors and percentiles. In practice, we note that this quantization slightly lowers the 25% percentile of the registration errors (as we "snap" many displacements to the correct voxel or slice) but increases the 50% and 75% percentiles (some landmarks get "snapped" to the wrong slice). Please note that in the Supplementary Material, we do not include this quantization step for e.g. ablation studies: this allows us to study more precisely the impact of each layer in our architecture.

Volume vs point cloud representation. We stress that our point clouds contain much less information than the original 3D volumes from which they have been sampled. We discard all the intensity (grayscale) values and only retain the sparse geometric support of the lung vascular tree -60k points out of 100M+ voxels. As a consequence, our registration task on 3D point clouds is significantly harder than the original DirLab benchmark (https://www.dir-lab.com/Results.html). Whereas optimization-based methods on the full CT volumes reach a nearly perfect accuracy of 0.60 mm to 1.00 mm with run times on the order of the minute [START_REF] Castillo | Computing global minimizers to a constrained B-spline image registration problem from optimal L1 perturbations to block match data[END_REF][START_REF] Castillo | A reference dataset for deformable image registration spatial accuracy evaluation using the COPDgene study archive[END_REF][START_REF] Heinrich | Estimating large lung motion in COPD patients by symmetric regularised correspondence fields[END_REF][START_REF] Hermann | Evaluation of scan-line optimization for 3D medical image registration[END_REF][START_REF] Polzin | Memory efficient LDDMM for lung CT[END_REF][START_REF] Polzin | Combining automatic landmark detection and variational methods for lung CT registration[END_REF][START_REF] Rühaak | Highly accurate fast lung CT registration[END_REF][START_REF] Rühaak | Estimation of large motion in lung CT by integrating regularized keypoint correspondences into dense deformable registration[END_REF][START_REF] Vishnevskiy | Isotropic total variation regularization of displacements in parametric image registration[END_REF][START_REF] Vishnevskiy | Isotropic total variation regularization of displacements in parametric image registration[END_REF], our point neural networks reach an average accuracy of 2 mm to 4 mm in one or two seconds.

The main purpose of our work on the PVT1010 dataset is to show that fast and accurate registration is now at hand, even on very degraded anatomical data. This is of significant interest for clinical practice: our method is suitable for real-time processing and has built-in robustness to changes of the CT acquisition parameters that may affect the intensities of the raw image volumes. Going forward, as detailed in the conclusion of our manuscript, we intend to work on improving the accuracy of our method with a better sampling strategy and image-based features. Packaging our method as an accessible Python toolbox will also open the door to a genuine multi-center evaluation of our trained models.

A.2 Augmentation of the training dataset for vascular tree registration

We now describe how to augment the PVT1010 dataset with synthetic deformations in order to create a large training set with dense ground truth annotations for lung registration. We proceed in four steps: voxel-grid sampling; local deformation; global deformation; local property distortion and degradation using an inconsistent sub-sampling. Our efficient implementation lets us generate synthetic training pairs online -just like a standard data augmentation layer. We showcase our local and global deformations in Fig. 6.

1. Voxel-grid sampling. To start, we use a voxel-grid strategy to re-sample the raw point clouds with a standard sampling density. First, we subdivide the volume space into coarse 3D blocks with spacing s voxelgrid = 0.03 mm. Second, we sort all points into these cubic bins according to their (x, y, z) coordinates. Third, we compute one barycenter per cubic cell to down-sample the original point cloud. Since we work with weighted point clouds, these local centroids are associated to the sums of the weights of all points in the corresponding cells.

Local deformation.

We then sample control points from the vessel trees and generate random displacements that we smooth using an anisotropic spline model:

1. We first sample C = 1, 000 spline control points x c uniformly at random from the point cloud (x 1 , . . . , x N ). 2. Second, we compute local covariance matrices for the distribution of points x i in a neighborhood of each control point x c , using an isotropic Gaussian kernel window. 3. We compute the three eigenvalues (e 1 c , e 2 c , e 3 c ) and unit eigenvectors (v 1 c , v 2 c , v 3 c ) of each local covariance matrix to determine the main direction of the lung vessel. To avoid rank deficiency, we use a lower threshold of 0.2 on the eigenvalues e k c . For each control point x c , we then normalize the vector of three eigenvalues (e

Σ c = 3 k=1 (s local e k c ) 2 v k c v k c
, where s local = 4 mm is a positive scaling factor. 5. To obtain a robust estimation of the local covariance structure of our point cloud, we re-run steps 2-4 with neighborhoods that are defined using an anisotropic Gaussian kernel window of covariance Σ c . 6. For every control point x c , we generate a random displacement vector ∆x c ∈ R 3 such that ∆x c d local , drawn uniformly in the ball of center 0 and radius d local = 3 mm.

7. We smooth and interpolate the displacement vector field ∆x c from the control points x c to the full point cloud {x i } using an anisotropic Nadaraya-Watson kernel interpolator:

x i ← x i + C c=1 k Σc (x c , x i )∆x c C c=1 k Σc (x c , x i ) , (9)
where k Σc (x c , •) is a Gaussian kernel with local covariance Σ c . This anisotropic formula ensures that the local connectivity structure of the lung vessel tree is preserved: the displacement of points that belong to the same vessel are strongly correlated with each other.

3. Global deformation. The step above simulates local relative displacements between lung vessels.

To take large-scale breathing movements into account, we apply a second spline deformation which is smoother but has a larger magnitude. In practice, we use a voxel-grid sampling with spacing resolution s global = 90 mm to obtain control points. We then generate random displacements of magnitude at most d global for every control point, and interpolate them to the full point cloud using a Nadaraya-Watson estimator, parameterized by an isotropic Gaussian kernel with standard deviation σ global .

4.a. Radius distortion.

Having altered the 3D coordinates of our points using the deformations above, we add random noise to the local estimates α i of the vessel radii -which are encoded as additional point features as detailed in Sec. A.1 and used in our RobOT layer as point weights. This additive noise is scaled by a positive parameter s radius = 0.1 and drawn at random in [-s radius α i , s radius α i ].

4.b. Inconsistent sampling. By construction, the steps above let us create an arbitrary number of pairs of (real, simulated) lungs with known pairwise correspondence between all points. In order to simulate acquisition artifacts and introduce challenging inconsistencies, we sample N = M = 60k points at random from the source and the synthetic (target) point clouds as a last generation step. We note that for training, the ground-truth flow is computed based on the source sampling and is not affected by this last degradation: it may or may not point to a sample in the synthetic target. We now discuss a toy example, illustrated in Fig. 7, that showcases the benefits of unbalanced optimal transport theory for partial rigid registration. We normalize and resample the partial Stanford scans using a voxel-grid sampling with 0.005 spacing. We then rotate every source shape by (120°,10°,10°) degrees to create a target shape. We illustrate two methods:

• As a baseline competitor, we use the FilterReg [START_REF] Gao | Filterreg: Robust and efficient probabilistic point-set registration using Gaussian filter and twist parameterization[END_REF] implementation of https://github. com/neka-nat/probreg (a feature-based CPD method [START_REF] Myronenko | Point set registration: Coherent point drift[END_REF]) with noise ratio set to 0.7 and an automatic annealing strategy based on the Expectation-Maximization (EM) algorithm.

• For our rigid S-RobOT registration, we first compute a 33-dimensional FPFH feature vector [START_REF] Rusu | Fast point feature histograms (FPFH) for 3D registration[END_REF] for each point using a radius of 0.02 for the normal search and a radius of 0.05 for the feature search. We then rely on Eqs. (4-5) to compute the weighted RobOT matching and project it onto the space of rigid transformations. For our robust OT problem, we set the blur parameter to 0.1 and the reach parameter to 2 units in R 33 .

Results. Given standard FPFH features [START_REF] Rusu | Fast point feature histograms (FPFH) for 3D registration[END_REF], our RobOT approach successfully registers the lowoverlap pair while FilterReg [START_REF] Gao | Filterreg: Robust and efficient probabilistic point-set registration using Gaussian filter and twist parameterization[END_REF] fails. The RobOT confidence weights w i of Eq. ( 4) naturally act as an attention mechanism.

Experiment 2: ModelNet40. Going further, we evaluate our partial registration strategy S-RobOT on the standard dataset ModelNet-40 [START_REF] Wu | 3d shapenets: A deep representation for volumetric shapes[END_REF] with comparisons to state-of-the-art methods. Instead of relying on FPFH features, we use a self-supervised deep feature learning strategy to increase the robustness and accuracy of the registration. The ModelNet40 dataset contains 12,311 CAD (Computer-Aided Design) models, from 40 object categories. We follow the same experimental setup as in [START_REF] Dang | Learning 3D-3D correspondences for one-shot partial-to-partial registration[END_REF][START_REF] Wang | PRNet: Self-supervised learning for partial-to-partial registration[END_REF]:

Partial-to-Partial registration with unseen categories. We then follow [START_REF] Wang | PRNet: Self-supervised learning for partial-to-partial registration[END_REF] and test the generalization ability of our model between object categories: we train on the first 20 categories of ModelNet40 and test on the remaining 20 categories. We report these results in Tab. 3.

Results. Overall, we observe that the unsupervised Rigid S-RobOT method performs much better than traditional approaches and is close to end-to-end methods. Since our main focus is on freeform registration, we do not push these experiments further: we use Affine and Rigid S-RobOT as pre-alignment steps and are satisfied with this level of performance.

We note that the best-performing method Partial-OT also relies on an optimal transport layer: as discussed in [START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF], partial OT is a very close cousin of the theory of unbalanced OT that we leverage in our RobOT layer. We understand the Partial-OT method as a supervised and end-to-end version of our S-RobOT baseline. Both [START_REF] Dang | Learning 3D-3D correspondences for one-shot partial-to-partial registration[END_REF] and this manuscript thus show that optimal transport theory is ready to be part of the standard toolbox in our field, with state-of-the-art performance in varied and complementary application settings.

A.3.2 Diffeomorphic registration

Background on LDDMM. Going beyond rigid and affine transformation models, the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework captures large, smooth and invertible deformations in a principled way [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]. This fluid-based model is standard in computational anatomy, especially for applications to neuroimaging where the preservation of shape topology is a key registration prior [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF].

In the LDDMM model, deformations are encoded via the integration of a time-varying velocity field v t = d dt x t over the unit time interval t ∈ [0, 1]. Given any two shapes A and B to register with each other, we seek a geodesic path (v t ) t∈[0,1] for a chosen Riemannian metric v 2 K that is induced by a convolution kernel K x over the space of vector fields v : R 3 → R 3 .

Following standard derivations from optimal control theory [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF][START_REF] Hart | An optimal control approach for deformable registration[END_REF], we know that plausible deformations are fully parameterized by the momentum m 0 : R 3 → R 3 at time t = 0, a vector field that acts as the parameter "θ" of the LDDMM deformation model. In the LDDMM framework, the optimization problem of Eq. ( 8) reads:

θ * = m 0 * = arg min m 0 :R 3 →R 3 λ reg m 0 , K x * m 0 L 2 (R 3 ,R 3 ) Regularization. + N i=1 w i x i + v i -Morph(m 0 , x i ) 2 R 3
Fidelity to the data. , [START_REF] Bookstein | Principal warps: Thin-plate splines and the decomposition of deformations[END_REF] where the deformation model Morph : (m 0 , x 0 i ) → x 1 i is computed through the integration of the geodesic shooting equation:

d dt x t = + ∂H ∂m x t , m t , d dt m t = -∂H ∂x x t , m t (11)
from time t = 0 to time t = 1 for the Hamiltonian

H(x, m) = 1 2 m, K x * m L 2 (R 3 ,R 3
) . We refer to Chapter 5 of [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF] for a presentation and implementation of these equations on point clouds with the PyTorch and KeOps libraries.

Black-box deformation models. Different transformation models may result in different projection results. In Fig. 8, we show several smooth RobOT (S-RobOT) results for spline and LDDMM models on a toy registration task.

Experiment. The source sphere (left) and the target cube (right) are sampled with N = 1, 922 and M = 1, 538 points respectively. We normalize their coordinates to be contained in [-1, 1] 3 and normalize the point cloud weights to sum up to one (α i = 1/1, 922 and β j = 1/1, 538 so that 1,922

i=1 α i = 1,538 j=1 β j = 1)
. For the initial RobOT matching, we set the blur parameter to 0.005 units in R 3 and the reach parameter to +∞ (balanced OT). For the following spline and LDDMM regularizations, we use a multi-Gaussian-kernel with standard deviations {0.05, 0.2, 0.3} and weights {0.2, 0.3, 0.5}, i.e. a kernel function:

k(x, y) = 0.2 • exp -x -y 2 R 3 /(2 • 0.05 2 ) + 0.3 • exp -x -y 2 R 3 /(2 • 0.2 2 ) (12)
+ 0.5 • exp -x -y 2 R 3 /(2 • 0.3 2 ) .
Results. We make two important observations:

• As evidenced by the disappearance of the polar sampling pattern in the fourth column, the raw RobOT matching x i → x i + v i provides a "perfect fit" to the target but does not preserve the topology of the source point cloud. The spline and the LDDMM approximations alleviate this problem: they smooth the weighted RobOT matching to combine accuracy with topology preservation and, in the case of the LDDMM model, guarantees of invertibility.

• The spline and the LDDMM model both result in smooth deformations. But crucially, the LDDMM model is more suited to large deformations and provides a better fit to the corners of the target cube. We note that workarounds exist for the over-smoothing of the second (spline) column: at the cost of an increased sensitivity to noise, singular kernels [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF] or ridge regression methods [START_REF] Bookstein | Principal warps: Thin-plate splines and the decomposition of deformations[END_REF] allow spline models to get a closer fit to the target. For medical applications, the key benefit of LDDMM is that it provides strong guarantees on the invertibility of the deformation: we refer to Chapter 5 of [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF] for a detailed discussion.

In the remainder of this work, we benchmark both spline and LDDMM deformation models whenever relevant. We observe that diffeomorphic LDDMM registrations perform better on e.g. lung data, but stress that such comparisons are task-dependent.

A.3.3 Deep feature learning

Feature learning on 3D point clouds. As discussed above, (x, y, z) coordinates or standard FPFH features can be good enough to handle simple shapes and deformations. On challenging settings however, learning task-specific features is often key to high performance.

In Sec. 3.2 and Suppl. A.4, we present our best-performing solution to this problem: the D-RobOT architecture. In this section, we discuss an alternative approach that decouples feature learning and matching. In practice, we did not obtain competitive results with this method. Nevertheless, using a separate feature extractor may increase interpretability and ease deployment issues in medical scenarios: we believe that this line of work is worth pursuing and provide full details on our experiments.

Let us process a synthetic pair of point clouds with N points in correspondence with each other -N = 60, 000 in our experiments for lung registration. Instead of sampling positive and negative samples, which is common for contrastive loss functions, we choose to rely on N × N correspondence matrices that take all possible point pairs into account. We use the KeOps library [START_REF] Feydy | Fast geometric learning with symbolic matrices[END_REF] to manipulate these objects efficiently, with extremely fast run times and without memory overflows.

Specifically, we train our feature extraction network as follows:

1. Input data. We assume that we are given two point clouds (x 1 , . . . , x N ) and (y 1 , . . . , y N ) that are in pairwise correspondence with each other. In our experiments, these are typically the output of a synthetic "data augmentation" procedure: the Flying3D objects for the Kitti benchmark and our synthetic lung pairs for the Dirlab benchmark. 2. Feature extraction. We apply the feature extractor (a trainable point neural network) on both point clouds, independently from each other. We retrieve point features p i and q i in R D that are respectively associated to the points x i and y i . 3. Feature normalization. As discussed in Sec. 2.1, we normalize the feature vectors so that

p i R D = q i R D = 1.
This prevents the feature extractor from converging to degenerate solutions and ensures that our hyper-parameters for the RobOT problem of Eq. ( 2) can be interpreted as sensible scales in R D . 4. Source self-similarities. We compute the N × N correspondence matrix for the point positions in the source point cloud:

c (xi,xj ) = softmax N j=1 (-x i -x j 2 R 3 /2κ 2 ) = exp(-x i -x j 2 R 3 /2κ 2 ) N j=1 exp(-x i -x j 2 R 3 /2κ 2 ) , (13)
where the Softmax denotes a mirrored exponential followed by a normalization over the rows of the correspondence matrix while κ > 0 is a scaling factor. Each row of the matrix c (xi,xj ) then refers to a probability distribution, a position heatmap with a peak at point x i whose radius is proportional to κ. We adapt this feature pyramid network from PointPWC-Net [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]. We implement a four-scale U-net [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] architecture using: point-wise multi-layer perceptrons on the feature embeddings (with LeakyReLU non-linearities); PointConv [START_REF] Wu | PointConv: Deep convolutional networks on 3D point clouds[END_REF] for downsampling and gathering information from the neighbors; K-Nearest Neighbor interpolation for upsampling. Please note that this architecture relies on the point coordinates (x, y, z) to define the PointConv convolution and the upsampling layer. In our experiments for lung registration, we use the local vessel radius as our only input feature: this corresponds to using a single input channel on the left-most layer of the figure above.

5. Feature correspondences. Similarly, we compute a correspondence matrix for the source and target features:

c (pi,qj ) = softmax N j=1 (-p i -q j 2 R D ) . (14)
Each row of c (pi,qj ) is a feature heatmap that indicates how well p i corresponds to q j . 6. Training loss. Our total loss is the sum over the cross entropies for each row:

CE(c (xi,xj ) , c (pi,qj ) ) = - N i=1 N j=1 c (xi,xj ) log c (pi,qj ) . (15)
We optimize it by stochastic gradient descent over the parameters of the feature extraction network (step 2).

Hyper-parameters. For feature learning, we use a U-net structured PointConv [START_REF] Wu | PointConv: Deep convolutional networks on 3D point clouds[END_REF] architecture with 7M parameters that is illustrated in Fig. 9. At the feature learning stage, we set κ (described above) to √ 2 mm. For each point we learn a 15 dimensional feature vector of unit length. For the RobOT feature matching, we set the blur parameter to σ = 0.01 units in R 15 and the reach parameter to ρ = +∞ (balanced OT). As discussed in Suppl. A.4.2, we benchmark two types of S-RobOT regularization:

1. The spline smoothing of Eq. ( 7), using a Gaussian (RBF) kernel with a standard deviation of 5 mm. 2. An LDDMM deformation model that we optimize as in Eq. ( 8) using an SGD solver. We use a multi-Gaussian-kernel with standard deviations {5, 8, 10} mm and weights {0.2, 0.3, 0.5}.

A.4 Additional material on deep deformation prediction

We now provide full details on our D-RobOT architecture and additional experiments on scene flow estimation and lung registration.

A.4.1 Deep registration module

Architecture of the prediction network. In Sec. 3.2, we rely on a modified PointPWC-Net to act as a predictor Pred : (x i , y j ) → θ. This multiscale point neural network takes as input the source and target point clouds. It returns a high-dimensional parameter θ for the deformation model Morph : (θ, x i ) → ŷi . In this work, the predicted θ is always a 3D vector field that is supported by the points x i (for the raw displacements model) or by a collection of control points c i that have been generated by farthest point sampling (for the spline and LDDMM models).

We describe our modifications to the original PointPWC-Net architecture in Fig. 10. In order to define a network that can predict spline and LDDMM parameters:

1. We replace the flow prediction layer of PointPWC-Net by a suitable prediction layer for registration parameters. 2. We use an asymmetric hierarchical architecture that outputs registration parameters for the control points c i . Since the number of control points is often much smaller than the number of points in the input point cloud, this reduces the memory footprint of our network architecture.

Losses. We use synthetic data pairs for training: for scene flow estimation on Kitti, we rely on the synthetic Flying3D dataset; for lung registration, we rely on the two-scales simulator of Suppl. A.2.

In both cases, we thus have access to ground truth deformations and can rely on a supervised learning strategy. Let us consider a PointPWC-Net with L scales, and denote by (x l 1 , . . . , x l N l ) the subsampled input for the l-th scale. For training, we compute a multi-scale similarity loss as:

Loss(ŷ i ) = L-1 l=0 W l i w i • ŷl i -y l i 2 2 , (16)
where l ∈ {0, . . . , L} is a scale (l = 0 corresponds to the raw point clouds), W l is a scalar hyperparameter that we use as a total weight for the l-th scale, ŷl i is the output of the network that corresponds to the flowed x l i at the l-th scale and y l i denotes the ground truth target that corresponds to the same point x l i with weight w i . Note that at the finest scale, we compute ŷ0 using a task-specific deformation model, Morph : (θ, x i ) → ŷi .

A.4.2 Experiments on lung vascular trees

We now provide more details on our experiments for lung registration, which are illustrated in Fig. 11 with a visualization of landmark errors in Fig. 13. Notably, we discuss both supervised and unsupervised training strategies as well as the influence of the deformation model Morph : (θ, x i ) → ŷi on the D-robot architecture.

Synthetic data vs unsupervised learning. To overcome the lack of dense 3D annotations on real shape data, we advocate the use of simulated deformations and synthetic training datasets. As detailed in Sec. 3.2, we observe that D-RobOT networks are easy to train to a high level of accuracy: RobOT-based post-processing and fine-tuning help our models to bridge the domain gap between the synthetic and real distributions of shapes.

We would like to stress that our focus on synthetic training datasets to the detriment of e.g. unsupervised approaches results from careful and extensive experiments. In the lead-up to the publication of this work, we tried several competing approaches to train our registration networks: supervised learning with dense correspondences on synthetic data; "unsupervised" learning on unannotated point clouds using geometric loss functions between point sets; a mix of both approaches. In practice, supervised learning on synthetic data clearly emerged as the most practical option for challenging registration tasks. Let us briefly explain why.

Unsupervised loss functions. When dense pointwise correspondences are not available, a popular strategy to train registration methods is to rely on permutation-invariant loss functions between point clouds. We refer to Chapters 3 and 4 of [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF] for an introduction to the topic.

Since our PVT1010 dataset contains 1,000 pairs of lung vessel trees that are in correspondence with each other (inspiration/expiration) but for which no expert-annotated landmarks are available, we are in a perfect situation to try out these tools. We thus attempted to train our networks using local Laplacian matching [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF], Maximum Mean Discrepancies, Gaussian Mixture Models (GMM) as well as Wasserstein distances [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF] on (x, y, z) coordinates. Unfortunately, we never succeeded in converging to a competitive level of accuracy. We believe that this is due to the complex geometric structure of the lung vascular trees, which are significantly more intricate than the clean point clouds and surface meshes on which these methods are usually tested [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF][START_REF] Feydy | Global divergences between measures: from Hausdorff distance to optimal transport[END_REF].

Combining synthetic (supervised) and real (unsupervised) data. Going further, we also tried to use a mixed strategy: in our training dataset, we combined synthetic deformations of real source shapes (that can be handled using a mean square error) with genuine target point clouds (that can be Figure 10: Architecture of our deep prediction network Pred : (x i , y j ) → θ, which is adapted from PointPWC-Net [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]. The original PointPWC-Net architecture is a multi-scale network that predicts 3D scene flow (a raw displacement field) in a coarse-to-fine fashion. We modify this stateof-the-art architecture to make it suitable for registration with a task-specific deformation model Morph : (θ, x i ) → ŷi . a) First of all, we compute feature pyramids at scales L0 to L4 using farthest point sampling and PointConv layers [START_REF] Wu | PointConv: Deep convolutional networks on 3D point clouds[END_REF] on the source and target shapes (left and right columns) -see Fig. 5 in [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF].

In our experiments, we always use C = 4, 096 or C = 8, 192 control points. The key design decision behind the PointPWC-Net architecture is to run through this pyramid from the coarsest scale (L4) to the finest one (L0, the original point clouds) in order to predict the final 3D scene flow: we represent this coarse-to-fine prediction as a stack of PointPWC blocks (central column). b) Those blocks share the same architecture (but not the same neural weights): an upsampling layer that interpolates our 3D vector fields from the coarser to the finer scale using an inverse distance spline kernel; another upsampling layer that interpolates the pyramid features at the coarsest scale; a warp layer that deforms the source shape according to the up-sampled 3D flow, allowing the network to focus on residual deformations at each scale; a cost "volume" layer, detailed in Fig. 2 of [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF], which relies on K-NN neighborhoods and PointConv layers to merge the source and target features into a vector of "patchwise" discrepancies for every point in the (subsampled and warped) source shape; a prediction PointConv and Multi-Layer Perceptrons (MLP), detailed in Fig. 6 of [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF], that turns this vector of discrepancies into a predicted correction for the up-sampled 3D flow. Please note that this architecture relies on the point coordinates (x, y, z) in the PointConv and upsampling layers. For our scene flow experiments, we use the (x, y, z) coordinates as input to the network: this corresponds to "Source" and "Target" arrays that have 3 input channels. For our lung registration experiments, we also add the local vessel radius as a fourth feature and thus use 4 input channels. We stress that in our RobOT-based pre-alignment and post-processing (steps 1 and 3 of the D-RobOT architecture), we rely on the (x, y, z) coordinates as point features p i and q j in Eq. ( 2). When available, the vessel radii are only used as point weights α i and β j in the RobOT problem. In all our experiments, the parameter θ is a 3D vector field that has the same memory footprint as a "raw" 3D scene flow supported by our control points: this corresponds to using 3 output channels on the C control points. Our main modification to the original PointPWC-Net architecture is to handle this registration parameter (purple) in parallel with the usual scene flow (sky blue): the final prediction module of every PointPWC block runs two estimations in parallel, which share the same PointConv layer and are equivalent to Fig. 6 in [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF].

handled using e.g. the Wasserstein distance). We expected that this combination would narrow the domain gap between our synthetic deformations and the real breathing movement. In practice, this strategy produced indecisive results:

• On the one hand, assuming that we only have access to a simple deformation model (e.g. a single-scale random field), this mixed training strategy improves the accuracy of our registrations. This is especially true when the synthetic deformations are not diverse enough. • On the other hand, assuming that we have access to the more realistic and expressive deformation model of Suppl. A.2, the introduction of real but not annotated pairs in the training loop proves slightly detrimental to performance. We observe a small increase of about 0.2mm in landmark Root Mean Squared Error (RMSE).

For the sake of simplicity, we thus trained our final D-RobOT network entirely on synthetic deformations. We found that using RobOT as a final post-processing layer in the D-RobOT pipeline is enough to address prediction errors effectively and compensate for a small domain gap between synthetic and real deformations.

Raw displacements vs spline and LDDMM regularizations. In Fig. 12, we further compare the D-RobOT results for three different deformation models. Our main observations are that:

• A deep prediction module that returns raw displacements ∆x i with a naive deformation model Morph : (θ = ∆x i , x i ) → ŷi = x i + ∆x i may produce non-smooth registration results. This holds even when the network is trained entirely on smooth, synthetic deformations. On the other hand, the regularizing spline and LDDMM models always result in smooth deformations. This vindicates the use of explicit regularizing models in safetycritical applications: we cannot trust our networks to "learn smoothness properties" from the data. • According to Fig. 12 and Tab. 4, the LDDMM model captures large deformations better than the spline model. However, we obtain similar performance with both spline and LDDMM models after the (affordable) RobOT post-processing step: fine-tuning with RobOT alleviates the need for complex and expensive deformation models.

Experiments. We now provide detailed hyperparameters for the D-RobOT architecture in our lung registration experiments. As detailed in Sec. 3.2, the full model applies successively: an affine S-RobOT pre-alignment; a deep non-parametric registration; and a spline S-RobOT post-processing.

1. Affine registration. For affine pre-alignment, we use the S-RobOT model of Eq. ( 6) with raw (x, y, z) coordinates as input features in R 3 . The vessel radii are taken into account as point weights α i and β j . In the RobOT problem of Eq. ( 2), we set the blur parameter to σ = 1 mm and the reach parameter to τ = +∞ (balanced OT with strong constraints on the marginals).

2. Deep non-parametric registration. We evaluate three deformation models Morph : (θ, x i ) → ŷi :

1. Raw displacements correspond to the case where the parameter θ = ∆x i is a 3D vector field that is supported by the point x i and the deformation model is the addition:

Morph : (∆x i , x i ) ∈ R N×3 × R N×3 → x i + ∆x i ∈ R N×3 . ( 17
) As a regularization penalty, we use the squared Euclidean norm:

Reg(∆x i ) = 1 N N i=1 ∆x i 2 R 3 .
(18) 2. The spline model corresponds to the case where the parameter θ = ∆c l is a 3D vector field that is supported by a set of control points (c 1 , . . . , c C ) in R 3 . As a deformation model, we use a simple kernel smoothing:

Morph : (∆c l , x i ) ∈ R C×3 × R N×3 → x i + C l=1 k(x i , c l )∆c l ∈ R N×3 . ( 19
) As a regularization penalty, we use the squared Euclidean norm:

Reg(∆c l ) = 1 C C l=1 ∆c l 2 R 3 .
(20) For our lung experiments, we always use a multi-Gaussian kernel with standard deviations {3, 6, 9} mm and weights {0.2, 0.3, 0.5}. With the coordinates of points x and y in millimeters, this corresponds to the positive definite kernel function:

k(x, y) = 0.2 • exp -x -y 2 R 3 /(2 • 3 2 ) + 0.3 • exp -x -y 2 R 3 /(2 • 6 2 ) (21) + 0.5 • exp -x -y 2 R 3 /(2 • 9 2 ) . 29
3. The diffeomorphic LDDMM model is equivalent to a spline deformation with continuous time. Notably, we use the same kernel k(x, y) as in the spline model above. As detailed in Suppl. A.3.2 and Chapter 5 of [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF], the parameter θ = m 0 l is a 3D vector field that is supported by a set of control points (c 1 , . . . , c C ) in R 3 . As a deformation model, we use an iterative kernel smoothing and deformation layer that corresponds to a numerical integration scheme for the Ordinary Differential Equation [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF]. More specifically, we use the differentiable Dopri-5 layer that is provided by the TorchDiffEq library [START_REF] Chen | Learning neural event functions for ordinary differential equations[END_REF][START_REF] Chen | Neural ordinary differential equations[END_REF], which implements the Dormand-Prince or Runge-Kutta 4(5) integrator [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF] with 20 time steps from time t = 0 to time t = 1. For a simpler implementation, we refer to the Euler integration scheme that is detailed in [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF], Algorithm 5.6. As a regularization term, we use the squared kernel norm:

Reg(∆c l ) = 1 C C l=1 C s=1 k(c l , c s ) m 0 l , m 0 s R 3 . (22)
For the deep prediction network Pred : (x i , y j ) → θ, we use the 4-scale hierarchical architecture that is described in Suppl. A.4.1 and set the scale weights W l to {1.0, 0.8, 0.4, 0.2}. We compute the 4-scale decompositions of the point clouds using Farthest Point Sampling. For all models, we promote a close fit to the target shape and weight the regularization term Reg(θ) by 1/100 in the training loss:

Loss(x i , y i , θ) = 1 100 Reg(θ) + L-1 l=0 W l i w i • ŷl i (θ, x i ) -y l i 2 R 3
Multiscale loss of Eq. ( 16)

.

We recall that in the expression above:

1. The weights w i are attached to the points x i and are proportional to the vessel radius.

2. The weights W 0 , W 1 , W 2 and W 3 put a strong emphasis on the finest scales.

3. Each point y l i corresponds to the ground truth target for x i at scale l. 4. Each point ŷl i (θ, x i ) corresponds to the output of our prediction and deformation networks at scale l for x i . As detailed in Fig. 10, the finest-scale output ŷ0 i (θ, x i ) = Morph(θ, x i ) corresponds to our model-based deformation; the larger-scale predictions ŷ1 i (θ, x i ), ŷ2 i (θ, x i ) and ŷ3 i (θ, x i ) correspond to raw displacements and are only used in Eq. ( 23) to make the training easier [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF].

3. Postprocessing. To fine-tune our registration, we use the spline S-RobOT deformation of Eq. ( 7) with (x, y, z) coordinates as input features in R 3 . We set the blur parameter to σ = 0.1 mm and the reach parameter to τ = 10 mm. For smoothing, we use the vessel-preserving anisotropic kernel k(x, y) of Sec. A.2 with a kernel scale s local = 8 mm. In order to get a closer fit to the target, we apply this post-processing step twice for all lung registration results. Table 4: Ablation study on the post-processing module. We assess the influence of the fine-tuning step in the D-RobOT model. As an addendum to Tab. 1, we display the Root Mean Squared Error in millimeters (median of the 10 subjects in parentheses) for expert-annotated landmarks on the 10 Dirlab lung pairs, that we use as a test set for PVT1010. The first row corresponds to the output of the deep registration module, without any fine-tuning. The second row corresponds to a nearest neighbor projection on the target point cloud. The third row corresponds to the "raw" RobOT matching of Eq. 4, extrapolated to the landmarks using a Nadaraya-Watson spline with an isotropic Gaussian kernel of deviation σ = 0.5 mm. In the fourth and fifth rows, we use a vessel-preserving anisotropic kernel to smooth a nearest neighbor projection (fourth row) and our RobOT matching (fifth row). Table 6: Ablation study on the modules of D-RobOT (with spline deformations) and the influence of the point sampling rate, performed on the 57 largest scenes from the Kitti dataset. As detailed in Suppl. A.4.3, these results correspond to average values on the 57 (out of 142) pairs of 3D scenes that are sampled with more than 30k points per frame in the original Kitti dataset. This allows us to focus our evaluation on the most challenging 3D scenes, with under-sampling artifacts: as expected, performance is lower than in the "full" Kitti benchmark of Fig. 4.

D-RobOT

  Source and target frames, sampled with 30k points each. (c) Exhale and inhale trees, sampled with 60k points each. (a) Source and target points in dimension D = 2.

Figure 1 :

 1 Figure 1: Robust Optimal Transport (RobOT) generalizes sorting to spaces of dimension D 1.(a) RobOT is equivalent to a nearest neighbor projection subject to mass distribution constraints that make it robust to translations and small deformation. We demonstrate that RobOT is now ready to be part of the standard toolbox in computer vision with extensive numerical experiments for 3D scene flow estimation (b) and lung registration (c). Rendering done with Paraview[START_REF] Ahrens | Paraview: An end-user tool for large data visualization[END_REF] and PyVista[START_REF] Sullivan | PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (vtk)[END_REF].

  2.a. Deep prediction network. 2.b. Smooth deformation model. 3. Spline-RobOT post-processing.

Figure 2 :

 2 Figure 2: The D-RobOT architecture. Left: We apply three successive registration modules that bring the moving source shape increasingly close to the fixed target point cloud. On the one hand, the RobOT-based pre-alignment (1) and fine-tuning (3) steps take the (x, y, z) coordinates as input features and do not require any training. On the other hand, our deep registration module (2) relies on a multi-scale point neural network "Pred : (x i , y j ) → θ" (2.a) and a task-specific deformation model "Morph(θ, x i ) → ŷi " (2.b). We train it end-to-end on a dataset of synthetic pairs of shapes with known ground truth correspondences. Right: To generate these pairs, we apply random deformations to real source shapes. For lung registration, we apply successively a vessel-preserving local perturbation and a smooth global deformation -as detailed in Suppl. A.2.

3 .

 3 Predicted registrations may be less accurate than optimization-based solutions. Training registration networks to pixel-perfect accuracy is notoriously hard, with diminishing returns in terms of model size and number of training samples.

Figure 3 :

 3 Figure 3: Influence of the sampling density on a detail (top car) of the Kitti frames of Fig. 1.b. First row: On sub-sampled 3D scenes, the regularizing priors of the D-RobOT architecture prevent over-fitting to the random sampling patterns of the target point cloud (blue). The D-RobOT output (last column, red) is very close to the ground truth scene flow (green). Second row: Increasing the number of points per frame reduces the influence of sampling artifacts. The simple RobOT baseline (fourth column, red) still over-fits to the target (blue) but becomes remarkably competitive.

Figure 4 :

 4 Figure 4: Evaluation on the Kitti dataset for 3D scene flow. Black numbers and crosses (×) correspond to results on scene pairs that are sampled with 8,192 points per frame; red numbers and plus signs (+) correspond to scene pairs that are sampled with 30,000 points per frame.

  (a) Source. (b) Target. (c) RobOT. (d) S-RobOT and target. (e) Deep features (source, target). (f) After pre-alignment, target. (g) After deep registration, target. (h) After fine-tuning, target.

Figure 5 :

 5 Figure 5: Registration of lung vascular trees. First row: S-RobOT registration with the deep features of Suppl. A.3.3. We display: (a) the source and (b) target shapes; (c) the "raw" RobOT registration; (d) a smoother S-RobOT registration with spline regularization, an overlap of this result with the target shape; (e) a visualization of the deep features on a pair of lung vascular trees, with colors that correspond to a t-SNE embedding of the point features in color space [116]. Second row: D-RobOT registration. We display the successive steps of our model: (f) pre-alignment with affine S-RobOT; (g) deep registration with an LDDMM model; (h) fine-tuning with spline S-RobOT.

  Augmentation of the the target and source point clouds. In our experiments, we generate our target point clouds using d global = 25 mm, σ global = 25 mm. Additionally, we also perform data augmentation for the source point cloud itself using the smaller values of d global = 8 mm and σ global = 15 mm. (a) Real source. (b) Local kernels. (c) Coarse control points. (d) Synthetic target. (e) Local perturbation, source. (f) Global deformation, source.

Figure 6 :

 6 Figure 6: Local and global deformations that we use to generate our synthetic training dataset. (a) Original (real) vascular tree. (b) Ellipsoids that represent anisotropic kernels of size 2 mm that we use to generate vessel-preserving local deformations. Note that in our experiments, we use larger kernels of size 4 mm that induce a stronger regularization but are harder to display cleanly. (c) Spline control points that we sample using a voxel-grid scheme with 90 mm spacing and use to generate a global deformation. (d) Synthetic vascular tree, the output of the process that we use as a target for training. (e) Source point cloud after local deformation. (f) After global deformation.

Figure 7 :

 7 Figure 7: Partial rigid registration of the Stanford scans based on the matching of FPFH featuresas discussed in Sec. 3.1 and Suppl. A.3.1. The feature-based CPD model FilterReg falls in a local minimum while the rigid S-RobOT of Eq. (5) results in a successful registration. In the last column, we display the attention weights w i that are derived from unbalanced OT.

Figure 8 :

 8 Figure 8: Smooth RobOT (S-RobOT) on a toy registration problem, the deformation of a sphere (left) onto a cube (right). From left to right, we display: (a) the source shape A = (x 1 , . . . , x N ); (b) the output of Spline RobOT from Eq. (7); (c) the output of LDDMM RobOT, solution of the optimization problems of Eq. (8,10); (d) the output of the "raw" RobOT matching x i → x i + v i from Eq. (4); (e) the target shape B = (y 1 , . . . , y M ).

Figure 9 :

 9 Figure9: Deep feature extractor for our S-RobOT experiments. We adapt this feature pyramid network from PointPWC-Net[START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]. We implement a four-scale U-net[START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] architecture using: point-wise multi-layer perceptrons on the feature embeddings (with LeakyReLU non-linearities); PointConv[START_REF] Wu | PointConv: Deep convolutional networks on 3D point clouds[END_REF] for downsampling and gathering information from the neighbors; K-Nearest Neighbor interpolation for upsampling. Please note that this architecture relies on the point coordinates (x, y, z) to define the PointConv convolution and the upsampling layer. In our experiments for lung registration, we use the local vessel radius as our only input feature: this corresponds to using a single input channel on the left-most layer of the figure above.

  Multiscale prediction architecture. (b) Architecture of a PointPWC block.

Figure 11 :Figure 12 :

 1112 Figure 11: Additional results for lung vascular trees. Columns refer to: (a) the source shape (expiration); (b) the registration result from D-RobOT using the LDDMM deformation model; (c) the target shape (inspiration); (d) an overlap between the D-RobOT registration and the target. Each row corresponds to a patient, with names that refer to case IDs in the original DirLab dataset.

Figure 13 :

 13 Figure13: Distributions of the DirLab-COPDGene errors for the 300 expert-annotated corresponding landmarks. Each image refers to a patient, with a title that refers to the case ID in the original DirLab dataset. We use color to display the registration error with our best performing model (D-RobOT with LDDMM deformations) on the source landmarks. The COPD 2 patient is a clear outlier in our experiment as in all previous studies on the original volumetric data (https://www.dir-lab.com/Results.html) -this may be due to annotation errors. On the other patients, most of the large errors occur in boundary regions where acquisition artifacts induce inconsistencies between the source and target point clouds: small vessels are harder to catch in the original 3D volume.

Figure 14 :

 14 Figure 14: Additional results for scene flow estimation on Kitti. From left to right, the columns refer to: the source shape; the registration result for D-RobOT with a spline deformation model (the source point cloud is displayed in gray in the background, flows are displayed as brown arrows); the target shape (the source point cloud is displayed in gray in the background); and the D-RobOT result overlapped with the target. Each row corresponds to a single pair of 3D frames.

Figure 15 :

 15 Figure 15: RobOT matching vs Nearest Neighbor projection on a synthetic local vessel tree.From left to right, the columns refer to: the source shape; the registration result for NN projection and RobOT followed by a local smoother (the source point cloud is displayed in gray in the background, flows are displayed as gray arrows); the target shape (the source point cloud is displayed in gray in the background); and the D-RobOT result overlapped with the target.

Table 1 :

 1 3D registration errors on the 3,000 expert-annotated DirLab landmarks. The † symbol denotes methods that are based on image keypoints[START_REF] Hansen | Learning deformable point set registration with regularized dynamic graph CNNs for large lung motion in COPD patients[END_REF][START_REF] Hansen | Tackling the problem of large deformations in deep learning based medical image registration using displacement embeddings[END_REF] and are evaluated on the original DirLab image dataset; all other approaches are tested on our point clouds.Data augmentation. The imbalance between our large training set and the small collection of 10 test cases reveals a fundamental challenge in computational anatomy: annotating pairs of 3D medical shapes with dense correspondences is prohibitively expensive. To work around this problem, we train our networks on synthetic deformations of the 600 × 2 lung shapes that make up our training set. As detailed in Suppl. A.2, we use a two-scales random field to generate a wide variety of deformations. This allows us to create a suitable training set with dense "ground truth" correspondences.

				Average		Percentiles
				error	25%	50%	75%	Time
			Method	mm ↓ mm ↓ mm ↓ mm ↓	s ↓
		No training	Input data ICP (affine) [7] CPD (non-rigid) [81] RobOT (affine) RobOT (raw)	23.30 13.18 22.22 31.65 15.05 9.60 14.06 20.01 9.30 5.95 8.60 11.83 332.60 -0.52 10.45 6.01 9.83 13.97 0.18 9.41 4.89 8.35 13.04 0.15
			DGCNN-CPD † [50]	4.30	-	-	-	-
	Due to its large memory requirements, CPD is tested on clouds of 20k points (in-	Supervised	DispEmd † [51] S-RobOT (spline) S-RobOT (LDDMM) D-RobOT (raw) D-RobOT (spline)	3.42 5.72 5.48 3.40 2.95	-3.19 2.86 1.40 1.30	-5.04 4.44 2.58 2.50	-7.35 7.14 3.69 3.19	-2.77 42.30 1.26 1.87
	stead of 60k).		D-RobOT (LDDMM)	2.86	1.25	2.23	3.11	1.92

  For each control point x c , we create a numerically stable anisotropic covariance matrix

	1 c , e 2 c , e 3 c ) as
	(e 1 c , e 2 c , e 3 c )/ (e 1 c ) 2 + (e 2 c ) 2 + (e 3 c ) 2 .
	4.

Table 2 :

 2 Partial-to-Partial registration with unseen objects on ModelNet40. ↓ degrees ↓ degrees ↓ 3D units 2 ↓ 3D units ↓ 3D units ↓

				Rotation			Translation
			MSE	RMSE	MAE	MSE	RMSE	MAE
	Method ICP FGR [134] Go-ICP [130] degrees Matching 1134.552 126.288 195.985 S-RobOT (rigid) 2.411	33.683 11.238 13.999 5.814	25.045 2.832 3.165 1.245	0.0856 0.0009 0.0011 0.0002	0.293 0.030 0.033 0.014	0.250 0.008 0.012 0.009
	End-to-end	PointNetLK [2] DCP(v2) [121] PRNet [122] Partial-OT [28]	280.044 45.005 10.235 0.107	16.735 6.709 3.199 0.328	7.550 4.448 1.454 0.052	0.0020 0.0007 0.0003 3.384e-06	0.045 0.027 0.016 0.00183	0.025 0.020 0.010 0.0003
				Rotation			Translation
			MSE	RMSE	MAE	MSE	RMSE	MAE
	Method ICP FGR [134] Go-ICP [130] degrees 2 Matching 1217.618 98.635 157.072 S-RobOT(rigid) 50.997	34.894 9.932 12.533 7.142	25.455 1.952 2.940 4.012	0.086 0.0014 0.0009 0.0005	0.293 0.038 0.031 0.022	0.251 0.007 0.010 0.013
	End-to-end	PointNetLK [2] DCP(v2) [121] PRNet [122] Partial-OT [28]	526.401 95.431 24.857 0.127	22.943 9.769 4.986 0.357	9.655 6.954 2.329 0.069	0.0037 0.0010 0.0004 3.953e-06	0.061 0.034 0.021 0.002	0.033 0.025 0.015 0.0004

2

↓ degrees ↓ degrees ↓ 3D units 2 ↓ 3D units ↓ 3D units ↓

Table 3 :

 3 Partial-to-Partial registration with unseen categories on ModelNet40.

	A.3 Additional material on S-RobOT and global feature matching
	A.3.1 Partial registration
	Experiment 1: Stanford scans.

  (raw) D-RobOT (spline) D-RobOT (LDDMM)

	Post-processing	mm ↓	mm ↓	mm ↓
	Without post-processing	3.46 (3.13)	4.45 (4.10)	3.33 (3.09)
	NN projection	3.33 (2.96)	3.32 (3.05)	3.31 (2.96)
	RobOT	3.35 (3.04)	3.18 (3.01)	3.19 (3.01)
	NN projection + Smoothing	3.32 (2.95)	3.08 (2.82)	2.95 (2.64)
	RobOT + Smoothing	3.33 (2.99)	2.94 (2.71)	2.83 (2.40)

Table 5 :

 5 Ablation study for the D-RobOT model (with LDDMM deformations) on the lung registration task. We benchmark the Root Mean Squared Error on the 10 × 300 pairs of DirLab landmarks for a wide range of configurations. We investigate the influence of affine S-RobOT prealignment (column 2); spline S-RobOT fine-tuning (column 4); the sampling rate for the input data (column 5); and the number of control points for the LDDMM deformation model (column 6). As a backbone architecture for the deep prediction network, we use a FlowNet3D (rows 3-5) and the modified PointPWC-Net of Fig.10(rows 6-10).

		Method	Prealign Post EPE3D Acc3DS Acc3DR Outliers3D EPE2D Acc2D
			cm ↓	% ↑	% ↑	% ↓	px ↓	% ↑
	8192 points	PWC [126] PWC PWC D-RobOT D-RobOT	6.78 4.62 4.76 4.63 4.55	78.58 83.24 82.16 79.13 80.27	90.30 94.91 94.53 95.30 94.85	23.61 19.61 20.10 21.43 20.65	2.6484 82.36 2.1411 86.82 2.0716 87.59 2.3930 82.93 2.1404 85.24
	30k points	PWC [126] PWC PWC D-RobOT D-RobOT	7.90 5.94 3.99 2.94 2.15	61.33 72.62 86.44 92.86 95.74	89.75 92.26 95.02 98.59 98.96	29.36 25.79 18.59 16.05 12.93	4.0558 65.35 3.1817 72.82 1.9592 86.09 1.5733 91.86 1.1118 95.66

Acknowledgments and Disclosure of Funding

Research reported in this publication was supported by the National Heart, Lung, and Blood Institute of the National Institutes of Health under award numbers R01HL149877 and R01HL116473. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The authors would also like to thank the anonymous reviewers for their most valuable advice.

A.4.3 Experiments on Kitti

We now provide details for our experiments on scene flow estimation, which are illustrated in Fig. 14. Notably, we discuss the influence of the sampling rate for the source and target point clouds and provide an extensive ablation study for the modules of the D-RobOT model.

Influence of the sampling rate. We follow the pre-processing strategy of PointPWC-Net [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF]:

1. We train on points whose depth with respect to the acquisition device is smaller than 35 m. 2. We sample points from both of the source and target frames at random, in a noncorresponding manner.

In the main manuscript, we report results when the source and target point clouds are sampled with either 8,192 or 30k points at a time. We note that in a similar situation, when sampling 32,768 points from each frame, HPLFlowNet [START_REF] Gu | HPLFlowNet: Hierarchical permutohedral lattice FlowNet for scene flow estimation on large-scale point clouds[END_REF] reported an average EPE3D value of 10.87 cm over all 142 pairs of the Kitti dataset. In comparison, our D-RobOT model reaches a much higher accuracy (EPE3D = 2.23 cm) with only 30,000 points per frame.

Working with the 57 largest scenes from Kitti. We note that the number of points per 3D frame varies widely in the Kitti dataset. To investigate the potential negative effects from under-sampling, we evaluate our models separately on pairs where both frames contain more than 30k points: 57 out of the 142 Kitti scene pairs meet this requirement. When dealing with those 57 "most populated" frames, sub-sampling the scene to 30k points induces a net loss of geometric information.

We display our results for this subset of Kitti in Tab. 6: unsurprisingly, performance is lower than in Fig. 4, which reports results on the full Kitti dataset and thus includes the 85 pairs of frames for which 30k points are enough to work at full resolution. Nevertheless, our conclusions remain consistent: the D-RobOT method outperforms PointPWC-Net [START_REF] Wu | PointPWC-Net: A coarse-to-fine network for supervised and self-supervised scene flow estimation on 3D point clouds[END_REF] for both sampling rates.

Ablation study. Going further, we analyze the contributions of the pre-alignment and post-processing modules in the D-RobOT model. Our results are summarized in Tab. 6 and can be described as follows:

• The pre-alignment module improves results for both sampling rates. Since rigid transformations have few degrees of freedom, the estimation of Eq. ( 5) does not rely on fine details.

• In sharp contrast, the post-processing module only improves results for high-resolution points. In practice, our RobOT-based post-processing behaves as a fast local fitting to the target point cloud. In the upper half of the table, the target point cloud and the set of "ground truth" final positions for the source points have small overlap due to the small number of sampling points for these complex scenes: a local "pixel-perfect" fine-tuning is detrimental to performance. In the lower half of the table, we increase the sampling rate to 30k points per frame: the target point cloud and the set of ground-truth final positions for the source points have a much higher overlap. This allows the RobOT post-processing to increase the accuracy of the full model.

Experiments. We now provide detailed hyperparameters for the D-RobOT architecture in our Kitti experiments. As detailed in Sec. 3.2, the full model applies successively: a rigid S-RobOT pre-alignment; a deep non-parametric registration; and a spline S-RobOT post-processing.

1. Rigid registration. For rigid pre-alignment, we use the S-RobOT model of Eq. ( 5) with raw (x, y, z) coordinates as input features in R 3 . The vessel radii are taken into account as point weights α i and β j . In the RobOT problem of Eq. ( 2), we set the blur parameter to σ = 1 m and the reach parameter to τ = +∞ (balanced OT with strong constraints on the marginals).

2. Deep non-parametric registration. We use the prediction architecture and deformation models of Suppl. A.4.2. Our spline deformation model is based on a multi-Gaussian-kernel with standard deviations {20, 40, 60} cm as and weights {0.2, 0.3, 0.5}. With the coordinates of points the x and y in centimeters, this corresponds to the positive definite kernel function:

3. Postprocessing. To fine-tune our registration, we use the spline S-RobOT deformation of Eq. ( 7) with (x, y, z) coordinates as input features in R 3 . We set the blur parameter to σ = 5 cm and the reach parameter to τ = 80 cm. For smoothing, we use the vessel-preserving anisotropic kernel k(x, y) of Sec. A.2 with a kernel scale s local = 80 cm.

A.5 RobOT matching vs nearest neighbor projection.

As detailed in Sec. 2.2, we use the RobOT matching of Eq. ( 4) as a plug-in replacement for nearest neighbor (NN) projection. At an affordable computational cost, our RobOT layer enforces a local conservation of the point density: this geometric prior is relevant for many registration tasks and mitigates accumulation artifacts that are commonly found in nearest neighbor projections [START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF][START_REF] Feydy | Global divergences between measures: from Hausdorff distance to optimal transport[END_REF].

In accordance with the theory, we thus observe that RobOT matching outperforms nearest neighbor projection in a wide range of settings.

Local translations. When the point features p i and q j correspond to the (x, y, z) coordinates in R 3 , a key property of the Monge-Brenier map that is defined by Eq. ( 4) is that it perfectly retrieves translations and dilations [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Feydy | Geometric data analysis, beyond convolutions[END_REF][START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF]. Scene flow estimation is often close to this best case scenario for RobOT theory: after we remove points that correspond to the ground (a common pre-processing for this task), most of the displacements can be explained as local translations and small rotations of solid objects (e.g. cars or trees) in the frame of reference of the acquisition device.

On the Kitti benchmark of Fig. 4, we observe that a simple RobOT matching already performs very well for 3D scene flow estimation. Remarkably, and without any training, vanilla RobOT matching on high-resolution point clouds (30k points per frame) outperforms most pre-existing deep learning methods on medium-resolution point clouds [START_REF] Bône | Deformetrica 4: an open-source software for statistical shape analysis[END_REF]192 points per frame) in terms of speed, memory usage and accuracy. This is strong evidence that geometric methods deserve more attention from the computer vision community.

Complex structures. Going further, the mass distribution constraint of Eq. ( 2) is useful to register complex shapes with appendices and branches. This property has been discussed in depth in previous works: we refer to e.g. [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF] for the matching of five-fingered hand surfaces. The (soft) constraints on the marginals of the transport plan promote a bijective matching of the fingers, preventing the thumb and the index of the source shape from being both projected on the thumb of the target.

In Fig. 15, we provide a similar example for the local registration of branching vessels. We create a synthetic pair of source and target vessels with N = 1, 000 and M = 1, 200 points respectively. We normalize their (x, y, z) coordinates to be contained in the unit cube [0, 1] 3 and assign uniform weights α i = 1 and β j = 1 to each point -we work in a realistic unbalanced scenario. We compute the raw NN matching using a simple nearest neighbor projection from the source onto the target in R 3 .

For the RobOT matching, we also use raw (x, y, z) coordinates and the squared Euclidean metric; we set the blur parameter to 0.0005 units in R 3 and the reach parameter to 1 unit (unbalanced OT).

As detailed in Suppl. A.2, we regularize both of the NN and RobOT matchings using an anisotropic muti-Gaussian-kernel with standard deviations {0.03, 0.05, 0.07} and weights {0.2, 0.3, 0.5}.

As a post-processing step. As discussed above, RobOT matching is good at handling local translations, dilations and small free-form deformations. We propose to use it as a fast post-processing step in our D-RobOT architecture, presented in Sec. 3.2. In order to validate its utility, we perform an ablation study on the PVT1010/DirLab dataset for lung registration.

We benchmark increasingly suitable fine-tuning layers: no fine-tuning; nearest neighbor projection; RobOT matching; nearest neighbor projection with vessel-preserving smoothing; RobOT matching with vessel-preserving smoothing. Our results are detailed in Tab. 4: when used in combination with a smooth deformation model, these layers result in an increasingly higher accuracy. This confirms the main findings of our work: including sensible geometric priors in a modular deep learning architecture is the key to reliable state-of-the-art performance.

A.6 Computational Resources

For the evaluation of time and memory cost in Fig. 4, we compare all models on a Ubuntu server with a single 24GB NVIDIA Quadro RTX 6000 graphics processing unit (GPU) and a 10-core Intel(R) Xeon(R) Silver 4114 CPU @ 2.20GHz. We trained all of our networks on a single 24GB NVIDIA RTX 3090 GPU, except for the training of PointPWC-Net with 30k sampled points that was performed on a 48GB NVIDIA RTX A6000 GPU.

A.7 Potential for negative societal impact

Point cloud registration is a fundamental low-level task in computer vision and computer graphics. As a consequence, negative societal impact of our work may arise in a wide range of use cases. A first example is that of point cloud registration for autonomous driving: algorithm failure could lead to incorrect driving decisions. Similar concerns may arise when using these approaches to register facial scans for the purpose of person identification. In general, while our manuscript demonstrates improved performance and robustness over competing approaches, the possibility for registration failures should always be considered when applying our models in safety critical environments.