Supervised Fitting of Geometric Primitives to 3D Point Clouds

Introduction

Recent 3D scanning techniques and large-scale 3D repositories have widened opportunities for 3D geometric data processing. However, most of the scanned data and the models in these repositories are represented as digitized point clouds or meshes. Such low-level representations of 3D data limit our ability to geometrically manipulate them due to the lack of structural information aligned with the shape semantics. For example, when editing a shape built from geometric primitives, the knowledge of the type and parameters of each primitive can greatly aid the manipulation in producing a plausible result (Figure 1). To address the absence of such structural information in digitized data, in this work we consider the conversion problem of mapping a 3D point cloud to a number of geometric primitives that best fit the underlying shape.

Representing an object with a set of simple geometric components is a long-standing problem in computer vision. Since the 1970s , the fundamental ideas for tackling the problem have been revised by many researchers, even until recently . However, most of these previous work aimed at solving perceptual learning tasks; the main focus was on parsing shapes, or generating a rough abstraction of the geometry with bounding primitives. In contrast, our goal is set at precisely fitting geometric primitives to the shape surface, even with the presence of noise in the input.

For this primitive fitting problem, RANSAC-based methods remain the standard. The main drawback of these approaches is the difficulty of finding suitable algorithm parameters. For example, if the threshold of fitting residual for accepting a candidate primitive is smaller than the noise level, over-segmentation may occur, whereas a too large threshold will cause the algorithm to miss small pieces primitives. This problem happens not only when processing noisy scanned data, but also when parsing meshes in 3D repositories because the discretization of the original shape into the mesh obscures the accurate local geometry of the shape surface. The demand for careful user control prevents RANSAC-based methods to scale up to a large number of categories of diverse shapes.

Such drawback motivates us to consider a supervised deep learning framework. The primitive fitting problem can be viewed as a model prediction problem, and the simplest approach would be directly regressing the parameters in the parameter space using a neural network. However, the regression loss based on direct measurement of the parameter difference does not reflect the actual fitting error – the distances between input points and the primitives. Such misinformed loss function can significantly limit prediction accuracy. To overcome this, Brachmann et al. integrated the RANSAC pipeline into an end-to-end neural network by replacing the hypothesis selection step with a differentiable procedure. However, their framework predicts only a single model, and it is not straightforward to extend it to predict multiple models (primitives in our case). Ranftl et al. also introduced a deep learning framework to perform model fitting via inlier weight prediction. We extend this idea to predict weights representing per-point membership for multiple primitive models in our setting.

In this work, we propose Supervised Primitive Fitting Network (SPFN) that takes point clouds as input and predicts a varying number of primitives of different types with accurate parameters. For robust estimation, SPFN does not directly output primitive parameters, but instead predicts three kinds of per-point properties: point-to-primitive membership, surface normal, and the type of the primitive the point belongs to. Our framework supports four types of primitives: plane, sphere, cylinder, and cones. These types form the most major components in CAD models. Given these per-point properties, our differentiable model estimator computes the primitive parameters in an algebraic way, making the fitting loss fully backpropable. The advantage of our approach is that the network can leverage the readily available supervisions of per-point properties in training. It has been shown that per-point classification problems (membership, type) are suitable to address using a neural network that directly consumes a point cloud as input . Normal prediction can also be handled effectively with a similar neural network .

We train and evaluate the proposed method using our novel dataset, ANSI 3D mechanical component models with $17k$ CAD models. The supervision in training is provided by parsing the CAD models and extracting the primitive information. In our comparison experiments, we demonstrate that our supervised approach outperforms the widely used RANSAC-based approach with a big margin, despite using models from separate categories in training and testing. Our method shows better fitting accuracy compared to even when we provide the latter with much higher-resolution point clouds as input.

We propose SPFN, an end-to-end supervised neural network that takes a point cloud as input and detects a varying number of primitives with different scales.

Our differentiable primitive model estimator solves a series of linear least-square problems, thus making the whole pipeline end-to-end trainable.

We demonstrate the performance of our network using a novel CAD model dataset of mechanical components.

Related Work

Among a large body of previous work on fitting primitives to 3D data, we review only methods that fit primitives to objects instead of scenes, as our target use cases are scanned point clouds of individual mechanical parts. For a more comprehensive review, see survey .

RANSAC and its variants are the most widely used methods for primitive detection in computer vision. A significant recent paper by Schnabel et al. introduced a robust RANSAC-based framework for detecting multiple primitives of different types in a dense point cloud. Li et al. extended by introducing a follow-up optimization that refines the extracted primitives based on the relations among them. As a downstream application of the RANSAC-based methods, Wu et al. and Du et al. proposed a procedure to reverse-engineer the Constructive Solid Geometry (CSG) model from an input point cloud or mesh. While these RANSAC variants showed state-of-the-art results in their respective fields, their performance typically depends on careful and laborious parameter tuning for each category of shapes. In addition, point normals are required, which are not readily available from 3D scans. In contrast, our supervised deep learning architecture requires only point cloud data as input and does not need any user control at test time.

Network-based Primitive Fitting.

Neural networks have been used in recent approaches to solve the primitive fitting problem in both supervised and unsupervised settings. However, these methods are limited in accuracy with a restricted number of supported types. In the work of Zou et al. and Tulsiani et al. , only cuboids are predicted and therefore can only serve as a rough abstraction of the input shape or image. CSGNet is capable of predicting more variety of primitives but with low accuracy, as the parameter extraction is done by performing classification on a discretized parameter space. In addition, their reinforcement learning step requires rendering a CSG model to generate visual feedback for every training iteration, making the computation demanding. Our framework can be trained end-to-end and thus does not need expensive external procedures.

Supervised Primitive Fitting Network

Notice that we do not assume a consistent ordering of ground truth primitives, so we do not assume any ordering of the columns of our predicted $\mathbf{\hat{W}}$. In Section 3.1, we describe the primitive reordering step used to handle such mismatch of orderings. In Section 3.2, we present our differential model estimator for predicting primitive parameters $\{\mathbf{\hat{A}}_{k}\}$. In Section 3.3, we define each term in our loss function. Lastly, in Section 3.4, we describe implementation details.

Inspired by Yi et al. , we compute Relaxed Intersection over Union (RIoU) for all pairs of columns from the membership matrices $\mathbf{W}$ and $\mathbf{\hat{W}}$. The RIoU for two indicator vectors $\mathbf{w}$ and $\mathbf{\hat{w}}$ is defined as follows:

The best one-to-one correspondence (determined by RIoU) between columns of the two matrices is then given by Hungarian matching . We reorder the ground truth primitives according to this correspondence, so that ground truth primitive $k$ is matched with the predicted primitive $k$. Since the set of inputs where a small perturbation will lead to a change of the matching result has measure zero, the overall pipeline remains differentiable almost everywhere. Hence we use an external Hungarian matching solver to obtain optimal matching indices, and then inject these back into our network to allow further loss computation and gradient propagation.

2 Primitive Model Estimation

We can then define $\mathbf{\hat{A}}$ as the minimizer to the weighted sum of squared distances as a function of $\mathbf{A}$:

By solving $\frac{\partial\mathcal{E}_{\text{plane}}}{\partial d}=0$, we obtain $d=\frac{\sum_{i=1}^{N}\mathbf{w}_{i}\mathbf{a}^{\mathbf{T}}\mathbf{P}_{i,:}}{\sum_{i=1}^{N}\mathbf{w}_{i}}$. Plugging this into Equation 3 gives:

where $\mathbf{X}_{i,:}=\mathbf{P}_{i,:}-\textstyle\frac{\sum_{i=1}^{N}\mathbf{w}_{i}\mathbf{P}_{i,:}}{\sum_{i=1}^{N}\mathbf{w}_{i}}$. Hence minimizing $\mathcal{E}_{\text{plane}}(\mathbf{A};\mathbf{P},\mathbf{w})$ over $\mathbf{a}$ becomes a homogeneous least square problem subject to $\|\mathbf{a}\|=1$, and its solution is given as the right singular vector $\mathbf{v}$ corresponding to the smallest singular value of matrix $\left(\text{diag}\left(\mathbf{w}\right)\right)^{\frac{1}{2}}\mathbf{X}$. As shown by Ionescu et al. , the gradient with respect to $\mathbf{v}$ can be backpropagated through the SVD computation.

Sphere.

In the sphere case (also in the cases of cylinder and cone), the squared distance is not quadratic. Hence minimizing the weighted sum of squared distances over parameters as done in the plane is only available via nonlinear iterative solvers . Instead, we consider minimizing over the weighted sum of a different notion of distance:

Solving $\frac{\partial\mathcal{E}_{\text{sphere}}}{\partial r^{2}}=0$ gives $r^{2}=\frac{1}{\sum_{i=1}^{N}\mathbf{w}_{i}}\sum_{j=1}^{N}\mathbf{w}_{j}\|\mathbf{P}_{j}-\mathbf{c}\|^{2}$. Putting this back in Equation 6, we end up with a quadratic expression in $\mathbf{c}$ as a least square:

where $\mathbf{X}_{i,:}=2\left(-\mathbf{P}_{i,:}+\textstyle\frac{\sum_{j=1}^{N}\mathbf{w}_{j}\mathbf{P}_{j,:}}{\sum_{j=1}^{N}\mathbf{w}_{j}}\right)$ and $\mathbf{y}_{i}=\mathbf{P}_{i,:}^{\mathbf{T}}\mathbf{P}_{i,:}-\textstyle\frac{\sum_{j=1}^{N}\mathbf{w}_{j}\mathbf{P}_{j,:}^{\mathbf{T}}\mathbf{P}_{j,:}}{\sum_{j=1}^{N}\mathbf{w}_{j}}$. This least square can be solved via Cholesky factorization in a differentiable way .

Cylinder.

where $\mathbf{v}=\mathbf{p}-\mathbf{c}$. As in the sphere case, directly minimizing over squared true distance is challenging. Instead, inspired by Nurunnabi et. al. , we first estimate the axis $\mathbf{a}$ and then solve a circle fitting to obtain the rest of the parameters. Observe that the normals of points on the cylinder must be perpendicular to $\mathbf{a}$, so we choose $\mathbf{a}$ to minimize:

which is a homogeneous least square problem same as Equation 4, and can be solved in the same way.

Once obtaining the axis $\mathbf{a}$, we consider a plane $\mathscr{P}$ with normal $\mathbf{a}$ that passes through the origin, and notice the projection of the cylinder onto $\mathscr{P}$ should form a circle. Thus we can choose $\mathbf{c}$ and $r$ to be the circle that best fits the projected points $\{\text{Proj}_{\mathbf{a}}(\mathbf{P}_{i,:})\}_{i=1}^{N}$, where $\text{Proj}_{\mathbf{a}}(\cdot)$ denotes the projection onto $\mathscr{P}$. This is exactly the same formulation as in the sphere case (Equation 6), and can thus be solved similarly.

Cone.

where $\mathbf{v}=\mathbf{p}-\mathbf{c},\alpha=\arccos\left(\frac{\mathbf{a}^{\mathbf{T}}\mathbf{v}}{\|\mathbf{v}\|}\right)$. Similarly with the cylinder case, we use a multi-stage algorithm: first we estimate $\mathbf{a}$ and $\mathbf{c}$ separately, and then we estimate the half-angle $\theta$.

We utilize the fact that the apex $\mathbf{c}$ must be the intersection point of all tangent planes on the cone surface. Using the predicted point normals $\mathbf{\hat{N}}$, the multi-plane intersection problem is formulated as a least square similar with Equation 7 by minimizing

where $\mathbf{y}_{i}=\mathbf{\hat{N}}_{i,:}^{\mathbf{T}}\mathbf{P}_{i,:}$. To get the axis direction $\mathbf{a}$, observe that $\mathbf{a}$ should be the normal of the plane passing through all $\mathbf{N}_{i}$ if point $i$ belongs to the cone. This is just a plane fitting problem, and we can compute $\mathbf{a}$ as the unit normal that minimizes Equation 3, where we replace $\mathbf{P}_{i,:}$ by $\mathbf{\hat{N}}_{i,:}$. We flip the sign of $\mathbf{a}$ if it is not going from $\mathbf{c}$ into the cone. Finally, using the apex $\mathbf{c}$ and the axis $\mathbf{a}$, the half-angle $\theta$ is simply computed as a weighted average:

3 Loss Function

We define our loss function $\mathcal{L}$ as the sum of the following five terms without weights:

Each loss term is described below for a single input shape.

The primitive parameters can be more accurately estimated when the segmentation of the input point cloud is close to the ground truth. Thus, we minimize $(1-\text{RIoU})$ for each pair of a ground truth primitive and its correspondence in the prediction:

Point Normal Angle Loss.

For predicting the point normals $\mathbf{\hat{N}}$ accurately, we minimize the absolute cosine angle between ground truth and predicted normals:

The absolute value is taken since our predicted normals are unoriented.

Per-point Primitive Type Loss.

We minimize cross entropy $H$ for the per-point primitive types $\mathbf{\hat{T}}$ (unassigned points are ignored):

Fitting Residual Loss.

Most importantly, we minimize the expected squared distance between $\mathbf{S}_{k}$ and the predicted primitive $k$ parameterized by $\mathbf{\hat{A}}_{k}$ across all $k=1,\ldots,K$:

where $\mathbf{p}\sim U(\mathbf{S})$ means $\mathbf{p}$ is sampled uniformly on the bounded surface $\mathbf{S}$ when taking the expectation, and $D^{2}_{l}(\mathbf{p},\mathbf{\hat{A}})$ is the squared distance from $\mathbf{p}$ to a primitive of type $l$ with parameter $\mathbf{\hat{A}}$, as defined in Section 3.2. Note that every $\mathbf{S}_{k}$ is weighted equally in Equation 17 regardless of its scale, the surface area relative to the entire shape. This allows us to detect small primitives that can be missed by other unsupervised methods.

Note that in Equation 17, we use the ground truth type $\mathbf{t}_{k}$ instead of inferring the predicted type based on $\mathbf{\hat{T}}$ and then properly weighted by $\mathbf{\hat{W}}$. We do this because coupling multiple predictions can make loss functions more complicated, resulting in unstable training. At test time, however, the type of primitive $k$ is predicted as

Axis Angle Loss.

Estimating plane normal and cylinder/cone axis using SVD can become numerically unstable when the predicted $\mathbf{\hat{W}}$ leads to degenerate cases, such as when the number of points with a nonzero weight is too small, or when the points with substantial weights form a narrow plane close to a line during plane normal estimation (Equation 4). Thus, we regularize the axis parameters with a cosine angle loss:

where $\Theta_{t}(\mathbf{A},\mathbf{\hat{A}})$ denotes $|\mathbf{a}^{\mathbf{T}}\mathbf{\hat{a}}|$ for plane (normal), cylinder (axis), and cone (axis), and $1$ for sphere (so the loss becomes zero).

4 Implementation Details

Experiments

For training and evaluating the proposed network, we use CAD models from American National Standards Institute (ANSI) mechanical components provided by TraceParts . Since there is no existing scanned 3D dataset for this type of objects, we train and test our network by generating noisy samples on these models. From 504 categories, we randomly select up to 100 models in each category for balance and diversity, and split training/test sets by categories so that training and test models are from disjoint categories, resulting in 13,831/3,366 models in training/test sets. We remark that the four types of primitives we consider (plane, sphere, cylinder, cone) cover $94.0\%$ percentage of area per-model on average in our dataset. When generating the point samples from models, we still include surfaces that are not one of the four types. The maximum number of primitives per shape does not exceed $20$ in all our models. We set $K_{\text{max}}=24$ where we add $4$ extra columns in $\mathbf{\hat{W}}$ to allow the neural net to assign a small number of points to the extra columns, effectively marking those points unassigned because of the threshold $\epsilon_{\text{discard}}$.

From the CAD models, we extract primitives information including their boundaries. We then merge adjacent pieces of primitive surfaces sharing exactly the same parameters; this happens because of the difficulty of representing boundaries in CAD models, so for instance a complete cylinder will be split into a disjoint union of two mirrored half cylinders. We discard tiny pieces of primitives (less than 2% of the entire area). Each shape is normalized so that its center of mass is at the origin, and the axis-aligned bounding box for the shape is included in $$range along every axis. In experiments, we first uniformly sample$8192$points over the entire surface of each shape as the input point cloud ($N=8192$). This is done by first sampling on the discretized mesh of the shape and then projecting all points onto its geometric surface. Then we randomly apply noise to the point cloud along the surface normal direction in$[-0.01,0.01]$range. To evaluate the fitting residual loss$\mathcal{L}_{\text{res}}$, we also uniformly sample 512 points per primitive surface for approximating$\mathbf{S}_{k}$($M=512$).

2 Evaluation Metrics

We design our evaluation metrics as below. Each quantity is described for a single shape, and the numbers are reported as the average of these quantities across all test shapes. For per-primitive metrics, we first perform primitive reordering as in Section 3.1 so the indices for predicted and ground truth primitives are matched.

Segmentation Mean IoU: $\frac{1}{K}\sum_{k=1}^{K}\text{IoU}(\mathbf{W}_{:,k},\mathcal{I}(\mathbf{\hat{W}}_{:,k}))$, where $\mathcal{I}(\cdot)$ is the one-hot conversion.

Mean point normal difference: $\frac{1}{N}\sum_{i=1}^{N}\arccos\left(|\mathbf{N}_{i,:}^{\mathbf{T}}\mathbf{\hat{N}}_{i,:}|\right)$.

When the predicted primitive numbers is less than $K$, there will be less than $K$ matched pairs in the output of the Hungarian matching. In this case, we modify the metrics of primitive type accuracy, axis difference, and $\{S_{k}\}$ residual mean/std. to average only over matched pairs.

3 Comparison to Efficient RANSAC [28]

We compare the performance of SPFN with Efficient RANSAC and also hybrid versions where we bring in predictions from neural networks as RANSAC input. We use the CGAL implementation of Efficient RANSAC with its default adaptive algorithm parameters. Following common practice, we run the algorithm multiple times ($3$ in our all experiments), and pick the result with highest input coverage. Different from our pipeline, Efficient RANSAC requires point normals as input. We use the standard jet-fitting algorithm to estimate the point normals from the input point cloud before feeding to RANSAC.

We report the results of SPFN and Efficient RANSAC in Table 1. Since Efficient RANSAC can afford point clouds of higher resolution, we test it both with the identical $8k$ input point cloud as in SPFN (row 1), and with another $64k$ input point cloud sampled and perturbed in the same way (row 2). Even compared to results from high-resolution point clouds, SPFN outperforms Efficient RANSAC in all metrics. Specifically, both $\{\mathbf{S}_{k}\}$ and $\mathbf{P}$ coverage numbers with threshold $\epsilon=0.01$ show big margins, demonstrating that our SPFN fits primitives more precisely.

We also test Efficient RANSAC by bringing in per-point properties predicted by SPFN. We first train SPFN with only $\mathcal{L}_{\text{seg}}$ loss, and then for each segment in the predicted membership matrix $\mathbf{\hat{W}}$ we use Efficient RANSAC to predict a single primitive (Table 1, row 4). We further add $\mathcal{L}_{\text{type}}$ and $\mathcal{L}_{\text{norm}}$ losses in training sequentially, and use the predicted primitive types $\mathbf{\hat{t}}$ and point normals $\mathbf{\hat{N}}$ in Efficient RANSAC (row 5-6). When the input point cloud is first segmented with a neural network, both $\{\mathbf{S}_{k}\}$ and $\mathbf{P}$ coverage numbers for Efficient RANSAC increase significantly, yet still lower than SPFN. Notice that the point normals and primitive types predicted by a neural network do not improve the $\{\mathbf{S}_{k}\}$ and $\mathbf{P}$ coverage in RANSAC.

Figure 5 illustrates $\{\mathbf{S}_{k}\}$ coverage with $\epsilon=0.01$ for varying scales of ground truth primitives. Efficient RANSAC coverage improves when leveraging the segmentation results of the network, but still remains low when the scale is small. In contrast, SPFN exhibits consistent high coverage for all scales.

4 Comparison to Direct Parameter Prediction Network (DPPN)

We also consider a simple neural network named Direct Parameter Prediction Network (DPPN) that directly predicts primitive parameters without predicting point properties as an intermediate step. DPPN uses the same PointNet++ architecture that consumes $\mathbf{P}$, but different from SPFN, it outputs $K_{\text{max}}$ primitive parameters for every primitive type (so it gives $4K_{\text{max}}$ sets of parameters). In training, the Hungarian matching to the ground truth primitives (Section 3.1) is performed with fitting residuals as in Equation 17 instead of RIoU. Since point properties are not predicted and the matching is based solely on fitting residuals (so the primitive type might mismatch), only $\mathcal{L}_{\text{res}}$ is used as the loss function. At test time, we assign each input point to the closest predicted primitive to form $\mathbf{\hat{W}}$.

The results are reported in row 7 of Table 1. Compared to SPFN, both $\{\mathbf{S}_{k}\}$ and $\mathbf{P}$ coverage numbers are far lower, particularly when the threshold is small ($\epsilon=0.01$). This implies that supervising a network not only with ground truth primitives but also with point-to-primitive associations is crucial for more accurate predictions.

5 Ablation Study

We conduct ablation study to verify the effect of each loss term. In Table 1 rows 8-11, we report the results when we exclude $\mathcal{L}_{\text{seg}}$, $\mathcal{L}_{\text{norm}}$ (use jet-fitting normals computed from $64k$ points), $\mathcal{L}_{\text{res}}$, and $\mathcal{L}_{\text{axis}}$, respectively. The coverage numbers drop the most when the segmentation loss $\mathcal{L}_{\text{seg}}$ is not used (-$\mathcal{L}_{\text{seg}}$). When using point normals computed from $64k$ input point clouds instead of predicting them (-$\mathcal{L}_{\text{norm}}$+J*), the coverage also drops despite more accurate point normals. This implies that SPFN predicts point normals in a way to better fit primitives rather than to just accurately predict the normals. Without including the fitting residual loss (-$\mathcal{L}_{\text{res}}$), we see a drop in coverage and segmentation accuracy. Excluding the primitive axis loss $\mathcal{L}_{\text{axis}}$ not only hurts the axis accuracy, but also gives lower coverage numbers (especially $\{\mathbf{S}_{k}\}$ coverage). Row 12 ($\mathbf{\hat{t}}\rightarrow$ Est.) shows results when using predicted types $\mathbf{\hat{t}}$ in the fitting residual loss (Equation 17) instead of the ground truth types $\mathbf{t}$. The results are compatible but slightly worse than SPFN where we decouple type and other predictions in training.

6 Results with Real Scans

For testing with real noise patterns, we 3D-printed some test models and scanned the outputs using a DAVID SLS-2 3D Scanner. Notice that SPFN trained on synthesized noises successfully reconstructed all primitives including the small segments (Figure 6).

Conclusion

We have presented Supervised Primitive Fitting Network (SPFN), a fully differentiable network architecture that predicts a varying number of geometric primitives from a 3D point cloud, potentially with noise. In contrast to directly predicting primitive parameters, SPFN predicts per-point properties and then derive the primitive parameters using a novel differentiable model estimator. The strong supervision we provide allows SPFN to accurately predict primitives of different scales that closely abstract the underlying geometric shape surface, without any user control. We demonstrated in experiments that this approach gives significant better results compared to both the RANSAC-based method and direct parameters prediction. We also introduced a new CAD model dataset, ANSI mechanical component dataset, along with a set of comprehensive evaluation metrics, based on which we performed our comparison and ablation studies.

The authors wish to thank Chengtao Wen and Mohsen Rezayat for valuable discussions and for making relevant data available to the project. Also, the authors thank TraceParts for providing ANSI Mechanical Component CAD models. This project is supported by a grant from the Siemens Corporation, NSF grant CHS-1528025 a Vannevar Bush Faculty Fellowship, and gifts from and Adobe and Autodesk. A. Dubrovina acknowledges the support in part by The Eric and Wendy Schmidt Postdoctoral Grant for Women in Mathematical and Computing Sciences.

References

Supplementary Material

In our differentiable model estimator, we are solving two linear algebra problems: homogeneous least square and unconstrained least square. In both problems, numerical stability issues can occur.

We solve the homogeneous least square using SVD to find $\mathbf{v}$, the right singular vector corresponding to the smallest singular value. However, when backpropagating the gradient through SVD , the gradient value goes to infinity when the singular values of the input matrix are not all distinct. In our case, such an issue happens only when the output segment (decided by membership matrix $\mathbf{\hat{W}}$) becomes degenerate. For instance, when fitting a plane to a segment via SVD, non-distinct singular values correspond to the case where the points in the segment with significant weights concentrate on a line or a single point. Hence if we get good segmentation by minimizing the segmentation loss (Section 3.3 ), then such degenerate cases should not happen. Thus, we handle the issue by simply bounding the gradient in the following way. We implemented a custom SVD layer following , and when computing $K_{ij}=\frac{1}{\sigma_{i}-\sigma_{j}}$ in Equation 13 of where $\sigma_{i}$, $\sigma_{j}$ are singular values, we instead use $K_{ij}=\frac{1}{\text{sign}(\sigma_{i}-\sigma_{j})\max(|\sigma_{i}-\sigma_{j}|,\epsilon)}$ for $\epsilon=10^{-10}$.

When solving the unconstrained least square using Cholesky factorization, numerical unstability can happen even when the segmentation is correct, but the type used in the estimator does not match with the segment. For instance, when fitting a sphere to a segment that is almost a flat plane, the optimal sphere is the one with center at infinity. To deal with such a singular case (as well as cases when the segments are degenerate), we add a $l_{2}$-regularizer to the formulation (Equation 7 ) and solve instead

with $\lambda=10^{-8}$. Even with such a modification, Cholesky factorization can still become unstable when the condition number of $\text{diag}(\mathbf{w})\mathbf{X}$ is too large, where the condition number of a matrix is defined to be the ratio of its largest singular value over its smallest singular value. To deal with this, we trivialize the least square problem when the condition number is larger than $10^{5}$ by setting $\mathbf{X}=\mathbf{0}$ to prevent gradient flow.

S.2 Training Details

We use the default hyperparameters for training PointNet++ with a batch size of $16$, initial learning rate $10^{-3}$, and staircase learning decay $0.7$. All neural network models in the experiments are trained for $100$ epochs, using Adam optimizer. The longest experiment (SPFN and its ablation studies) took $50$ hours to train on a single Titan Xp GPU, although the decay of the total loss was not substantial after $50$ epochs. We will release our source code and include a link to the code in the final version.

S.3 DPPN Architecture

The output of DPPN is simply a collection of $4K_{\text{max}}$ primitives including $K_{\text{max}}$ planes, $K_{\text{max}}$ spheres, $K_{\text{max}}$ cylinders, and $K_{\text{max}}$ cones. In order to compare with SPFN outputs, as a post-processing step, we construct auxiliary membership matrix $\mathbf{\hat{W}}$ by assigning each input point to the closest primitive among the $4K_{\text{max}}$ predicted primitives. Similarly, we construct per-point type matrix $\mathbf{\hat{T}}$ by assigning the type of each point to be the type of its closest primitive. The numbers reported in Table 1 are computed in the same evaluation pipeline as in SPFN after such post-processing step.

S.4 Primitive Correspondences

Tulsiani et al. and Sung et al. introduced a type of neural networks capable of discovering correspondences across different inputs without direct supervision. Notice in SPFN, changing the ordering of the columns in $\mathbf{\hat{W}}$ does not affect the loss. Despite such ambiguity, SPFN implicitly learns a preferred order such that the primitives represented by the same columns in $\mathbf{\hat{W}}$ in different shapes appear to be similar, resulting in rich correspondence information for primitives from different shapes (Figure S2). These results provide insight into the possible design variations for the same category of shapes.

S.5 Additional Experiments

To further study the capability of SPFN, we have conducted the following additional experiments.

PointNet++ used in our architecture has a limitation of handling high resolution point clouds during training time due to the increase of memory consumption. However, it is also known that PointNet++ is robust to the change of the resolution of point clouds at test time (See Section 3.3 in ). Hence, we can consider processing high resolution input point clouds in the test time by training the network with lower resolution point clouds. In Table S1, row 2 shows the results of testing $64k$ point clouds with the same SPFN model in Table 1 (trained with $8k$ point clouds), and it exhibits a slight improvement in nearly all metrics. We also assessed SPFN by adding not only noise in the inputs (as described in Section 4.1 ) but also outliers. Row 3 describes the results when we add $10\%$ outliers, which are uniformly sampled in space outside of the central cube $[-0.5,0.5]^{3}$, in both training and test data. The results show a little drop but still comparable performance. Lastly, we also investigated how robust SPFN can be if we change the maximum number of primitives, $K_{\text{max}}$. Row 4 illustrated the results when training the network with $K_{\text{max}}=48$, and we observed no substantial difference in performance.