BPNet: Bézier Primitive Segmentation on 3D Point Clouds

1 Introduction

Structuring and abstracting 3D point clouds via segmentation is a prerequisite for various computer vision and 3D modeling applications. Many approaches have been proposed for semantic segmentation, but the finite set of semantic classes limits their applicability. 3D instance-level segmentation and shape detection are much more demanding, while this literature lags far behind its semantic segmentation counterpart. Finding a generalized way to decompose point clouds is essential. For example, man-made objects can be decomposed into canonical primitives such as planes, spheres, and cylinders, which are helpful for visualization and editing. However, the limited types of canonical primitives are insufficient to describe objects' geometry in real-world tasks. We are looking for a generalized way of decomposing point clouds. The task of decomposing point clouds into different geometric primitives with corresponding parameters is referred to as parametric primitive segmentation. Parametric primitive segmentation is more reasonable than semantic instance segmentation for individual 3D objects, which unifies the 3D objects in the parametric space instead of forming artificially defined parts. However, the task is quite challenging as 1) there is no exhaustive repertoire of canonical geometric primitives, 2) the number of primitives and points belonging to that primitive may significantly vary, and 3) points assigned to the same primitive should belong to the same type of primitive.

∗Corresponding author

Inspired by the fact that B´ezier decomposition, where NURBS models can be divided into canonical geometric primitives (plane, sphere, cone, cylinder, etc.) and parametric surfaces into rational B´ezier patches, we propose to learn B´ezier decomposition on 3D point clouds. We focus on segmenting point clouds sampled from individual objects, such as CAD models. Departing from previous primitive segmentation, we generalize different primitive types to B´ezier primitives, making them suitable for end-to-end and batch training. To the best of our knowledge, our method is the only work to learn B´ezier decomposition on point clouds. To summarize our contributions:

1. We introduce a novel soft voting regularizer for the relaxed intersection over union (IOU) loss, improving our primitive segmentation results.
2. We design a new auto-weight embedding module to cluster point features which is free of iterations, making the network robust to real-scan data and work for axis-symmetric free-form point clouds.
3. We propose an innovative reconstruction module where we succeed in using a generalized formula to evaluate points on different primitive types, enabling our training process to be fully differential and compatible with batch operations.
4. Experiments demonstrate that our method works on the free-form point clouds and real-scan data even if we only train our model on the ABC dataset. Furthermore, we present one application of B´ezier primitive segmentation to reconstruct the full B´ezier model while preserving the sharp features. The code is available at: https://github. com/bizerfr/BPNet.

2 RelatedWork

B´ezier primitive segmentation involves parametric fitting, instance segmentation, and multi-task learning. We now provide a brief review of these related research areas.

Primitive segmentation. Primitive segmentation refers to the search and approximation of geometric primitives from point clouds. Primitives can be canonical geometric primitives, such as planes or spheres, or parametric surface patches, such as B´ezier, BSpline, or NURBS. We can classify primitive segmentation methods into two lines of approaches: geometric optimization and machine learning. Popular geometric optimization-based methods include RANSAC [Fischler and Bolles, 1981; Schnabel et al., 2007], region growing [Marshall et al., 2001] and Hough transforms [Rabbani et al., 2007]. We refer to [Kaiser et al., 2019] for a comprehensive survey. One limitation of geometric optimizationbased methods is that they require strong prior knowledge and are hence sensitive to parameters. In order to alleviate this problem, recent approaches utilize neural networks for learning specific classes of primitives such as cuboids [Zou et al., 2017; Tulsiani et al., 2017]. The SPFN supervised learning approach [Li et al., 2019] detects a wider repertoire of primitives such as planes, spheres, cylinders, and cones. Apart from the canonical primitives handled by SPFN, ParSeNet [Sharma et al., 2020] and HPNet [Yan and Yang, 2021] also detect open or closed BSpline surface patches. Nevertheless, different types of primitives are treated separately with insufficient genericity. This makes them unsuitable for batch operations, thus suffering long inference times. Deep learning-based methods are less sensitive to parameters but often support a limited repertoire of primitives. Our work extends SPFN, ParSeNet, and HPNet with more general B´ezier patches.

Instance segmentation. Instance segmentation is more challenging than semantic segmentation as the number of instances is not known a priori. Points assigned to the same instance should fall into the same semantic class. We distinguish between two types of methods: proposal-based [Yi et al., 2019; Yang et al., 2019; Engelmann et al., 2020] and proposal-free methods [Wang et al., 2018; Jiang et al., 2020; Huang et al., 2021]. On the one hand, proposal-based methods utilize an object-detection module and usually learn an instance mask for prediction. On the other hand, proposalfree methods tackle the problem as a clustering step after semantic segmentation. We refer to a recent comprehensive survey [Guo et al., 2020]. The significant difference between instance segmentation and primitive segmentation is that instance segmentation only focuses on partitioning individual objects where primitive fitting is absent.

Patch-based representations. Patch-based representations refer to finding a mapping from a 2D patch to a 3D surface. Previous works including [Groueix et al., 2018; Yang et al., 2018; Deng et al., 2020; Bednarik et al., 2020] learn a parametric 2D mapping by minimizing the Chamfer distance [Fan et al., 2017]. One issue with Chamfer distance is that it is not differentiable when using the nearest neighbor to find matched pairs. We learn the uv mapping instead. Learning uv parameters enables us to re-evaluate points from our proposed generalized B´ezier primitives, making our training process differentiable and supporting batch operations.

Multi-task learning. Multi-task learning aims to leverage relevant information contained in multiple related tasks to help improve the generalization performance of all the tasks [Zhang and Yang, 2021]. Compared to single-task learning, the architectures used for multi-task learning-see, e.g., [Zhang et al., 2014; Dai et al., 2016]-share a backbone to extract global features, followed by branches that transform the features and utilize them for specific tasks. Inspired by [Dai et al., 2016], we use a cascaded architecture for our joint optimization tasks.

3 Method

Figure 1 shows an overview of the proposed neural network. The input to our method is a 3D point cloud P = {pi|0 ≤ i ≤ N - 1}, where pi denotes the point coordinates (with or without normals). The output is the per-point patch labels {Pk| ∪k=0 Pk = P}, where each patch corresponds to a B´ezier primitive. The network will also output patch degree (du-by-dv) and weighted control points C = {ckmn = (x, y, z, w)|0 ≤ m ≤ du, 0 ≤ n ≤ dv, 0 ≤ k ≤ K - 1}, where K denotes the number of patches. We constrain the maximum degree to be Md ∗Nd. We let our network output a maximum number of K B´ezier patches for all CAD models, and we use ˆK to denote the ground-truth number of patches which is smaller than K and varies for each CAD model.

3.1 Architecture

Our architecture consists of two components: a backbone for extracting features and a cascaded structure for joint optimization. The backbone is based on three stacked EdgeConv [Wang et al., 2019] layers and extracts a 256D pointwise feature for each input point. Let P ∈ RN×Din denote the input matrix, where each row is the point coordinates (Din is three) with optional normals (Din is six). Let X ∈ RN×256 denote the 256D pointwise feature matrix extracted from the backbone. We use a cascaded structure to optimize the perpoint degree probability matrix D ∈ RN×(Md∗Nd) , the soft membership matrix W ∈ RN×K , the UV parameter matrix T ∈ RN×2 , and the weighted control points tensor C ∈ RK×(Md+1)×(Nd+1)×4 jointly. Because D, W, T, and C are coupled, it is natural to use a cascaded structure to jointly optimize them. Here, the cascaded structure is similar to [Dai et al., 2016], where the features are concatenated and transformed for different MLP branches.

3.2 Joint Optimization

We have four modules: decomposition, fitting, embedding, and reconstruction. They are coupled to optimize D, W, T and C jointly by using our proposed four modules.

Decomposition Module

Degree classification. We use B´ezier primitive with different degrees to replace classical primitives, including plane, sphere, plane, BSpline, etc. For the sake of the classification of degrees, the straightforward idea would be to use a crossentropy loss: CE = -log(pt), where pt denotes the possibility of the true degree labels. However, the degree type is highly imbalanced. For example, surfaces of degree type 1by-1 represent more than 50%, while 3-by-2 surfaces are rare. To deal with the imbalance, we utilize the multi-class focalloss [Lin et al., 2017]: FL = -(1 - pt)γ log(pt), where γ denotes the focusing parameter. Then the degree type classification loss is defined as:

Primitive segmentation. The output of primitive segmentation is a soft membership indicating per-point primitive instance probabilities. Each element wik is the probability for a point pi to be a member of primitive k. Since we can acquire pointwise patch labels from our data pre-processing, we use a relaxed IOU loss [Kr¨ahenb¨uhl and Koltun, 2013; Yi et al., 2018; Li et al., 2019] to regress theW:

where W denotes the output of the neural network and ˆW is the one-hot encoding of the ground truth primitive instance labels. The best matching pairs (k, ˆ k) between prediction and ground truth are found via the Hungarian matching [Kuhn, 1955]. Please refer to [Li et al., 2019] for more details.

Soft voting regularizer. Since we learn D and W separately, points belonging to the same primitive instance may have different degrees, which is undesirable. To favor degree consistency between points assigned to the same primitive, we propose a soft voting regularizer that penalizes pointwise degree possibilities. We first compute a score for each degree case for all primitive instances by S = WT D, where each element skd denotes the soft number of points for degree d in primitive instance k. We then perform L1-normalization to convert S into primitive degree distributions ˆS:

where the first term denotes the sum of each column and ⊙ denotes the element-wise product. Finally, we utilize a focal loss to compute the primitive degree voting loss:

where FL denotes the focal loss. The global loss for the decomposition module is defined as: Ldec = Ldeg + Lseg + Lvoting.

Fitting Module

Parameter regression. Through B´ezier decomposition we obtain the ground truth labels for the (u, v) parameters and record all parameters into matrix ˆT. We regress the uv parameters using a mean squared error (MSE) loss:

Control point regression. We select a maximum number of primitive instances K for all models. As the ground truth primitive instance ˆK varies for each model, we reuse the matching pairs directly from the Hungarian matching already computed in the primitive segmentation step. Note that as the predicted degree (du, dv) may differ from the ground truth ( ˆ du, ˆ dv), we align the degree to compute the loss via a maximum operation as (max(du, ˆ du),max(dv, ˆ dv)). The network always outputs (Md + 1) × (Nd + 1) control points for each primitive corresponding to the predefined maximum degree in U and V direction, and these control points will be truncated by the aligned degree. Furthermore, if the groundtruth degree is smaller than the prediction, we can pad "fake" control points that are zero for the ground-truth patch; otherwise, we just use the aligned degree, which is the maximum of the predicted and the ground truth. Finally, the control point loss is defined as:

where ct and ˆct denote the matched control points, and Nc is the number of matched control point pairs. Finally, we define the Lfit loss as: Lfit = Lpara +Lctrl.

Embedding Module

We use the embedding module to eliminate oversegmentation by pulling point-wise features toward their center and pushing apart different centers. Unlike ParSeNet and HPNet, 1) we do not need a mean-shift clustering step which is time-consuming; 2) we calculate the feature center in a weighted manner rather than simply averaging. The weights are chosen as W and will be automatically updated in the decomposition module; 3) W will be further optimized to improve the segmentation. Moreover, our embedding module is suitable for batch operations even though the number of primitive instances for each CAD model and the number of points for each primitive varies. Otherwise, one has to apply mean-shift for each primitive, which deteriorates timing further.

To be specific, we use W to weight X to obtain primitive features for all candidate primitive instances. Then, we reuse W to weigh all the primitive instance features to calculate a "soft" center feature for each point. We favor that each point feature embedding should be close to its "soft" center feature, and each primitive instance feature embedding should be far from each other. The primitive instance-wise feature matrix Xins is defined as:

where each row of Xins denotes the instance-wise features for each patch. We then compute the "soft" center feature matrix Xcenter as: Xcenter = WXins, where each row denotes the "soft" center for each point.

Then we define Lpull as:

and we define Lpush as:

Finally, the total embedding loss Lemb is defined as: Lemb = Lpull + Lpush.

Reconstruction Module

The reconstruction module is designed to reconstruct points from the predicted multiple B´ezier primitives, i.e., rational B´ezier patches, and further jointly optimize W. One difficulty is that each CAD model has various numbers of primitives, and the degree of each primitive is also different. Therefore, we seek a generalized formula to support tensor operations on re-evaluating points for a batch of CAD models. The straightforward approach would be to compute a synthesizing score for all degree types. Assume the maximum number of primitive instances is K, and we have Md∗Nd types of different degrees. The total number of combinations is K∗Md∗Nd. We define a synthesizing score for each case in Einstein summation form: (sw)kci = wik ∗ skc, where wik denotes the probability of point pi to belong to primitive instance k and skc denotes the degree score for degree type m-by-n indexed with c = M ∗ (m - 1) + (n - 1) for primitive instance k coming from S. Then, we need to normalize (sw)kdi such that ∑k,d,i (sw)kdi = 1. Finally, the reconstructed point coordinates pi are defined as:

where parameter (ui, vi) for point pi is shared for all combinations. Such a formulation makes extending the formula in matrix form easy and avoids resorting to loop operations. However, such an approach is too memory-intensive. We thus truncate the degree from the degree probability matrix by redefining the Bernstein basis function for degree d as:

where 0 ≤ l ≤ M, and M is the maximum degree. Then, the reconstructed point coordinates for pi for a degree m-by-n patch k is:

where cmini denotes the control point coordinates and (cw)mini denotes its weight, and wik is the element ofW.

If we also input the normal (nxi , nyi , nzi ) for point pi, we can also reconstruct the normal (n′xi , n′yi , n′zi ) by:

where × denotes the cross product.

pi denotes the input point coordinates. p∗i denotes the reconstructed point coordinates. npi denotes the input point normals. n∗pi denotes the reconstructed normals. The coordinate loss is defined as:

If we also input the normals, the normal loss is defined as:

The loss for the reconstruction module is defined as:

Total Loss

The total loss is defined as the sum of decomposition, fitting, embedding, and reconstruction losses: L = Ldec + Lfit + Lemb + Lrecon. We do not use different weights for each loss item because all point clouds are normalized into a unit sphere. Moreover, the uv parameters are outputted directly from a sigmoid layer, and the control points are outputted directly by a tanh layer. Thus, each loss item is almost at the same scale, so we do not need different weights for each loss item. Furthermore, we use different learning rates for different modules to balance the training. Specific training details are listed in section 4.2.

4 Experiments

4.1 Dataset Pre-Processing

We evaluate our approach on the ABC dataset [Koch et al., 2019]. However, the ABC dataset does not have the annotations to learn B´ezier decomposition on point clouds. Therefore, we do a pre-processing step. Specifically, we utilize the CGAL library [CGAL, 2009] and OpenCascade library [OpenCascade, 2018] to perform B´ezier decomposition on STEP files directly and perform random sampling on the surface to obtain the following labels: point coordinates, point normals, point uv parameters, surface patch indices of the corresponding points, surface patch degrees, and surface patch control points. Finally, we use 5,200 CAD models for training and 1,300 CAD models for testing. Each CAD model contains randomly sampled 8, 192 points (non-uniform) with annotations.

4.2 Training Details

We train a multi-task learning model. The learning rates differ depending on the MLP branch. The learning rate for the backbone, soft membership, and uv parameters is set to 10-3 , while the learning rate for the degree probabilities and control points is set to 10-4 . As we have several learning tasks that are not independent, we set a lower learning rate for loss items, such as degree probabilities which converges faster. We set γ as 3.0 for the focal loss, and δpull as 0 and δpush as 2.0 for the embedding losses. We employ ADAM to train our network. The model is then trained using 150 epochs.

4.3 Comparisons

We compare our algorithm with SPFN, ParSeNet, and HPNet [Li et al., 2019; Sharma et al., 2020; Yan and Yang, 2021]. We use both points and normals for training all the algorithms. Since SPFN only supports four types of canonical primitives (plane, sphere, cone, and cylinder), we consider points belonging to other primitives falling out of the supported canonical primitive types as the "unknown" type. To make fair comparisons, we modify SPFN to let the network take point coordinates and normals as input for training. For ParSeNet, we only train the segmentation module on the ABC dataset. We use their pre-trained fitting model (SplineNet) directly. For HPNet, we also use the pre-trained fitting model directly, which is the same as ParSeNet. We observed that the output of HPNet is very sensitive to the number of points. In order to use HPNet at its best, we down-sample the point clouds to 7k points for training and testing. We choose the following evaluation metrics:

1. Primitive Type Accuracy ("Acc"):

1K ∑K-1 k=0 I(tk == ˆ tk), where tk and ˆ tk are predicted primitive type and ground truth type, respectively. This is used to measure the type accuracy. Note that our primitive types differ from other baselines.

2. Rand Index ("RI"):

a+b c , where c is (N2 ) denoting the total possible pairs for all points, and a denotes the number of pairs of points that are both in the same primitive of prediction and ground truth, while b denotes the number of pairs of points that are in a different primitive of prediction and ground truth. Rand index is a similarity measurement between two instances of data clustering, and a higher value means better performance [Chen et al., 2009; Yi et al., 2018].

1. Normal Error ("Err"): 1N ∑N-1 i=0 arccos(|nTpi n∗pi |), where npi and n∗pi are ground truth and predicted unit normal, respectively.
2. Inference Time ("Time"): The inference time on the whole test dataset.

The predicted average number of primitives on the

1. Average Primitive Number ("Num"): whole test data set.

We record these evaluation metrics in table 1 and 2. Figure 2 shows visual depictions of the results. Our results show the best performance regarding primitive type accuracy, normal fitting error, and inference time. Our method is much faster for inference because it uses a general formula for different primitive types, and the embedding module is free of iterations. Other methods treat primitives with different equations, and ParSeNet and HPNet need a mean-shift step. Even though our approach may lead to more segmented primitives by the nature of B´ezier decomposition, the evaluation metrics of primitive type accuracy and normal fitting error are computed in a point-wise manner. Thus, over-segmentation and under-segmentation will not lead to smaller or bigger errors due to fewer or more segmented primitives.

We also show the performance of all the methods without normals as input. For our method and SPFN, we only input point coordinates into the neural networks but use normals as supervision. Since ParSeNet does not regress normals, we cannot use normals as supervision. We train ParSeNet without normals as input to test its performance. HPNet uses the network to regress the normals from the input and also utilizes the ground truth normals to construct an affinity matrix as a post-processing step for clustering. We modify HPNet to let the affinity matrix be constructed from the regressed normals instead of the ground-truth normals. Table 2 records the evaluation metrics of each method. From the experiments, we deduce that normals are important for the task of parametric primitive segmentation.

4.4 Ablation Studies

We first conduct experiments to verify the usefulness of the soft voting regularizer. The soft voting regularizer favors point primitive type consistency for each primitive instance, i.e., points assigned to the same primitive instance should have the same primitive type. From our experiment, we find that the soft voting regularizer not only improves the primitive type accuracy but also accelerates training relaxed IOU. Please refer to figure 3 and the last two rows of table 3.

We also verify the functionalities of each module. If we only use the decomposition module, the result is not good even though the "Acc" and "RI" are slightly higher because the decomposition module ignores the fitting, limiting the segmentation applicable to specific datasets. The reconstruction module reduces the "Err" significantly compared to the fitting module because the reconstruction module controls how "well-fitted" a predicted B´ezier primitive is to the input point clouds. In contrast, the fitting module only regresses the control points and uv parameters. The embedding module is designed to eliminate small patches that contain few points, seeing the "Num" column. Therefore, experimenting with the embedding module results in fewer patch numbers than its counterpart. To conclude, training with all the modules yields the best results.

4.5 Stress Tests

To test whether our algorithm can work in real-world scenarios, we show more results from the real-scan data from the Aim@Shape dataset [Falcidieno, 2004]. The sampling is non-uniform, with missing data and measurement noise compared to the ABC dataset. Besides, We cannot train the network on those data directly because they lack groundtruth labels. Instead, we use the models trained on the ABC dataset and test the performance on real-scan data. Our algorithm still works, while other methods are sensitive. Another positive aspect is that our algorithm could decompose the axis-symmetric free-form point clouds with much smoother boundaries of different patches. Please refer to figure 4.

We also test the performance of our network by adding Gaussian white noise. Specifically, we apply different scales of Gaussian white noise to the point coordinates after normalizing them into a unit sphere. The noise scale denotes the standard deviation of the Gaussian white noise. It ranges from 0.01 to 0.05. We train our network on noise-free data but test the network with Gaussian white noise. Please refer to table 4.

4.6 Applications

We can reconstruct the full B´ezier model from the B´ezier primitive segmentation. We do not follow ParSeNet to pretrain a model that outputs a fixed control point size. Instead, we reuse the rational B´ezier patch to refit the canonical B´ezier patch. We treat the degrees of the canonical B´ezier patch the same as the rational B´ezier patch. As a result, we fetch the segmentation and degrees of each patch predicted from the network. Then, we use the parameterization [L´evy et al., 2002] to recompute uv parameters and least squares to refit control points for each patch. Each patch is expanded by enlarging the uv domain to guarantee intersections with its adjacent patches. After that, we use the CGAL co-refinement package [CGAL, 2009] to detect intersecting polylines for adjacent tessellated patches and trim the tessellated patch with the intersected polylines. Our reconstructed full B´ezier model can preserve the sharp features, while the boundaries of ParSeNet for different primitives are jaggy and thus fail to preserve the sharp features. Please refer to figure 5.

5 Conclusion

This paper presents an end-to-end method to group points by learning B´ezier decomposition. In contrast to approaches treating different geometric primitives separately, our method uses a general formulation for different primitive types. Regarding limitations, B´ezier decomposition may naturally generate overly complex segmentations. In addition, we choose the rational B´ezier patch as the primitive type. As the formulation is not linear, fitting the parametric patch is not direct. In future work, we wish to use the neural network to directly regress the canonical B´ezier patch.

Acknowledgements

This research is part of a project that has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 860843. The work of Pierre Alliez is also supported by the French government, through the 3IA Cˆote d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.