Direct Fitting of Gaussian Mixture Models

TL;DR

This paper introduces a method for directly fitting Gaussian Mixture Models to 3D meshes, improving shape representation and registration accuracy over traditional point cloud approaches.

Contribution

It presents a novel formulation for fitting GMMs directly to triangular meshes, enhancing model quality and registration performance for 3D geometric data.

Findings

Improved 3D registration accuracy for meshes and RGB-D frames.

Enables fitting higher-quality GMMs with diverse initializations.

Applicable to various geometric objects and sensor data.

Abstract

When fitting Gaussian Mixture Models to 3D geometry, the model is typically fit to point clouds, even when the shapes were obtained as 3D meshes. Here we present a formulation for fitting Gaussian Mixture Models (GMMs) directly to a triangular mesh instead of using points sampled from its surface. Part of this work analyzes a general formulation for evaluating likelihood of geometric objects. This modification enables fitting higher-quality GMMs under a wider range of initialization conditions. Additionally, models obtained from this fitting method are shown to produce an improvement in 3D registration for both meshes and RGB-D frames. This result is general and applicable to arbitrary geometric objects, including representing uncertainty from sensor measurements.

Tables

Table 1. TABLE I : Results of repeating the experiments in Figures  3 and 6 on mulitple models from  [ 41 ] . All experiments used k-means initialization and K = 100 𝐾 100 K=100 . Details in Sec.  VI-C .

Model	Likelihood (larger is better)		Translation Error (% of ICP)		Rotation Error (% of ICP)
	points	mesh	points	mesh	points	mesh
Armadillo	-14.6	-12.2	127	37	161	33
Bunny	7.6	8.2	50	28	41	17
Dragon	6.9	7.6	68	25	40	19
Happy	7.3	8.2	101	27	85	27
Lucy	-21.5	-18.3	95	23	122	35

Equations

$$\mathcal{N}(x;\mu ,\Sigma )=2{\pi }^{-k/2}\det (\Sigma {)}^{-\frac{1}{2}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}$$

$$\log \mathcal{N}(x;\mu ,\Sigma )=$$

$$-\frac{1}{2}\log (\det (\Sigma ))$$

$$-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )$$

$$g(x)=\sum _{i=1}^K{\lambda }_i\mathcal{N}(x;{\mu }_i,{\Sigma }_i)$$

$$\begin{array}{rl}\ell & \cong \prod _{i=1}^{N}p(x_i)\\ \ell & \cong {\left(\prod _{i=1}^{N}p(x_i)\right)}^{\frac{1}{N}}\\ \ell & =\underset{N\to \infty }{\lim }{\left(\prod _{i=1}^{N}p(x_i)\right)}^{\frac{1}{N}}\\ \ell & =\underset{N\to \infty }{\lim }\exp \left(\log \left({\left(\prod _{i=1}^{N}p(x_i)\right)}^{\frac{1}{N}}\right)\right)\\ \ell & =\underset{N\to \infty }{\lim }\exp \left(\frac{1}{N}\sum _{i=1}^N\log (p(x_i))\right)\\ \ell & =\exp \left(\int \log (p(x))dx\right)\end{array}$$

$$\begin{array}{rl}L& =\exp \left(\sum _{j=1}^M{\int }_{△}\log \left(\sum _{i=1}^K{\lambda }_i\mathcal{N}(x;{\mu }_i,{\Sigma }_i)\right)dx\right)\\ \log (L)& =\sum _{j=1}^M{\int }_{△}\log \left(\sum _{i=1}^K{\lambda }_i\mathcal{N}(x;{\mu }_i,{\Sigma }_i)\right)dx\\ & \ge \sum _{j=1}^M\sum _{i=1}^K{\int }_{△}\log \left({\lambda }_i\mathcal{N}(x;{\mu }_i,{\Sigma }_i)\right)dx\end{array}$$

$$L\approx L_S=\prod _{j=1}^M{\left(\sum _{i=1}^K{\lambda }_i\mathcal{N}({\mu }_j;{\mu }_i,{\Sigma }_i)\right)}^{\frac{{\alpha }_j}{{\sum }_k{\alpha }_k}}$$

$$LB=\sum _{j=1}^M\sum _{i=1}^K{\eta }_{ij}\log ({\lambda }_i\mathcal{N}(x_{jk};{\mu }_i,{\Sigma }_i))$$

$$\begin{array}{rl}\frac{\partial LB}{\partial {\mu }_i}& =\frac{1}{2}\sum _{j=1}^M{\Sigma }_i^{-1}(x_j-u_i){\eta }_{ij}\\ \frac{\partial LB}{\partial {\Sigma }_i^{-1}}& =\frac{1}{2}\sum _{j=1}^M{\eta }_{ij}\left({\Sigma }_i-(x_j-{\mu }_i)(x_j-{\mu }_i{)}^T\right)\end{array}$$

$$w_{ij}={\eta }_{ij}{\alpha }_j$$

$$W_i=\sum _{j=1}^M{w}_{ij}$$

$$\log (L)$$

$$[2\log ({\lambda }_i)-k\log (2\pi )-\log (\det ({\Sigma }_i))$$

$$-({\mu }_j-{\mu }_i{)}^T{\Sigma }_i^{-1}({\mu }_j-{\mu }_i)$$

$$-\frac{1}{12}(A_j^T{\Sigma }_i^{-1}{A}_j+B_j^T{\Sigma }_i^{-1}{B}_j+C_j^T{\Sigma }_i^{-1}{C}_j$$

$$-3{\mu }_j^T{\Sigma }_i^{-1}{\mu }_j)]$$

$${\mu }_i=\frac{1}{W_i}\sum _{j=1}^M{w}_{ij}{\mu }_j$$

$${\Sigma }_i=\frac{1}{W_i}\sum _{j=1}^M{w}_{ij}\left[\underset{cov({\mu }_j,{\mu }_i)}{\underbrace{({\mu }_j-{\mu }_i)({\mu }_j-{\mu }_i{)}^T}}+\underset{cov({△}_j)}{\underbrace{\frac{1}{12}(A_j{A}_j^T+B_j{B}_j^T+C_j{C}_j^T-3{\mu }_j{\mu }_j^T)}}\right]$$

$$\log (L)$$

$$-\log (\det ({\Sigma }_i))-({\mu }_p-{\mu }_i{)}^T{\Sigma }_i^{-1}({\mu }_p-{\mu }_i)-{\Sigma }_p]$$

$${\mu }_i$$

$${\Sigma }_i$$

Full text

I Introduction

In robotics and computer vision, there exist many forms of representation for 3D geometric data. For example, some researchers use unordered point sets [1], others require points with surface normals [2] or dense volumetric representations such as signed distance fields [3]. The variation in forms of representation is related to the wide variety of sources and uses for this data, from raw depth sensor measurements to Computed-Aided Design (CAD) models.

Many researchers have found that Gaussian Mixture Models provide a powerful representation, especially for use in registration between unknown poses [4, 5, 6, 7, 8]. Producing a Gaussian Mixture Model requires only unstructured point sets from the underlying geometric data, which makes them widely applicable. Our contribution in this work is to demonstrate how Gaussian Mixture Model can be constructed, evaluated on and fit directly to geometric primitives, such as the triangles of a polygon mesh. This is done by incorporating the structural information from each primitive into the algorithm, for a visual example, see Figure 1.

Our contributions include a mathematical framework for how geometric primitives can be incorporated with probability distributions (Section II-B). We demonstrate how to obtain the structural properties for a triangular mesh (Section II) and how it can be generalized to other primitives (Section V-A). Incorporating structural information allows us to build Gaussian Mixture Models that not only converge faster and in more conditions (Section IV) but also provide a representation that produces higher quality registration results when the models are used in practice (Section VI).

The code for all methods and experiments in this paper is available at https://github.com/leonidk/direct_gmm.

II Method

II-A Gaussian Mixture Models

The Gaussian Mixture Model (GMMs) is a well studied probability distribution. It is possible to fit these models to empirical data via Maximum Likelihood Estimation (MLE) [9, 10]. The likelihood of any point $x$ in a Gaussian is given by

[TABLE]

where $\mu\in\mathbb{R}^{k}$ is the mean and $\Sigma\in\mathbb{R}^{k\times k}$ is the positive-definite covariance matrix. The log-likelihood of a given point $x$ is given by

[TABLE]

In a Gaussian Mixture Model with $K$ components, with $\lambda_{i}$ are mixing weights subject to $\sum_{i}\lambda_{i}=1$ and $\lambda_{i}\geq 0\ \forall i$, the probability of a point $x$ is given by

[TABLE]

II-B Geometric Objects in a Probability Distribution

First we must handle the general case of how to evaluate the probability of a known geometric object in a probability distribution. Consider sampling N points from the geometric object, where a notion of likelihood can be evaluated by taking their product. To account for the variable number of samples, we take a geometric mean. And to obtain the likelihood of the object, we take the limit as the number of samples grows to infinity.

[TABLE]

Equation 3 is a form of product integral [11, 12, 13], which can be used to evaluate the likelihood of a known geometric object. The extension to multiple objects is straightforward, simply adjust their sampling weights accordingly. This form of product integral has two nice properties, it is invariant to resampling and it produces a result of [math] if $p(x)=0$ anywhere along the geometric object. By invariance to resampling, we mean that one large primitive with sophisticated integration bounds gives the same answer as many small disjoint pieces of surface with their own bounds.

Following the use of Jensen’s inequality, we get the following lower bound on likelihood

[TABLE]

We also consider a simplified model, where each triangle is sampled at its center of mass ($\mu_{j}$), and has weight corresponding to its area ($\alpha_{j}$). As combining probabilities is done with multiplication, we use a weighted geometric mean over all points, obtaining the following approximation

[TABLE]

III Modifying EM maximization to account for triangles

In this section we derive the expressions for fitting a GMM to a triangular mesh. We will represent each GMM component with parameters $\mu_{i},\Sigma_{i},\lambda_{i}$ and each triangle with vertices $A_{j},B_{j},C_{j}$, centroid $\mu_{j}$ and area $\alpha_{j}$. Traditional EM minimization is possible by analyzing the lower bound of the log-likelihood. Following standard formulations (see [10, 14, 15]), we add mixture sampling probabilities $\eta_{ij}$ and move the logarithm inside the summation to obtain a valid lower-bound to minimize by using Jensen’s inequality.

To perform fitting, we need an M-step which obtains $\lambda_{i},\mu_{i},\Sigma_{i}$ by maximizing the lower bound

[TABLE]

To maximize this expression, we can take derivatives with respect to the variables of interest and obtain

[TABLE]

Now we can integrate these three expressions over the two dimensional surface of a triangle via a change of variables substitution, then set the result equal to zero and solve, thus obtaining the update equations. For clarity we will also define a weight variable $w_{ij}$ and its corresponding normalization constant $W_{i}$

[TABLE]

[TABLE]

The resulting lower-bound likelihood can be written as

[TABLE]

The new mean is obtained as simply weighted mean of centroids. This update equation is identical to the one derived for the approximation in equation 5

[TABLE]

The same technique will provide an answer to the update equation for covariance.

[TABLE]

The final update is surprisingly simple, it is the area weighted average of the covariance obtained by using centroids as point measurements plus the covariance equation for a triangle. That is, at every iteration, each mixture is updated with some fraction of the structure of the triangles associated with it. A visual example of ellipsoids showing triangle covariance structures is shown in Fig. 1. Our derived expression for the covariance of a triangle is expressed in terms of vertices, but is consistent with the standard formulation in CGAL [16].111The reference document has a typographic error in the moment matrix for triangles, which should be a 5x multiple of the one for 3D tetrahedrons. The CGAL source code correctly implements this matrix in practice. The update equation for eq. 5 is similar, simply lacking the $cov(\triangle_{j})$ term.

III-A Evaluating the derived loss function

To evaluate the validity of the expression in equation 8, we compare its fitting fidelity numerically against a large number of sampled points. The results are shown in figure 2. We can see that, the lower bound expression for triangles is equal to that obtained numerically from a large number of points. Since the lower-bound expression is all that’s needed in EM optimization [10], the equality of this expression suggests we can use it in fitting real data.

IV Results

We performed experiments fitting Gaussian Mixture Models to triangular meshes. We swept a wide range of $K$ mixture components (from 6 to 400) and evaluated two different initialization schemes. The first initialization performs k-means++ [17] clustering, and uses those clusters as initial assignments for the EM method. The second method uses simple random assignments for initialization. We run 25 iterations of EM for all methods, with a tight tolerance ($\epsilon=10^{-12}$) to prevent an early exit from the optimization.

To evaluate the converged model, we use a densely sampled point cloud of the initial mesh (Figure 5(e)) and report its likelihood according to equation 1. Since all of our experiments tend to operate on 1,000 points or triangles, the use of 50,000 points for evaluating the model should provide a good test of GMM model fidelity. In all of these cases, we focus on the Stanford Bunny. Visual examples of our input and output data is shown in Figure 5.

IV-A Mesh Input Data

The first set of experiments, shown in Figure 3, compares fitting GMM models to different input formats of the res4 Stanford Bunny. We compare our exact (eq. 8) and approximate (eq. 5) mesh loss equations against fitting a traditional point-loss to the triangle centroids and the mesh vertices. The best results came from our exact mesh expression, with the second best being its approximation. The proposed methods handled random initialization and k-means initialization. On the other hand, the point-based methods often had a preferred initialization. A qualitative look at the resulting models, shows that the mesh GMM (Fig. 5(h)) produces a fuller model of the Stanford Bunny than the center-of-mass GMM (Fig. 5(d)), even when the evaluated likelihood was numerically very similar.

IV-B Mesh Decimation

While the above experiments used the low-resolution res4 Stanford Bunny, we also evaluated which GMM fitting strategy works best when subsampling high-resolution mesh data. In our case, we try two methods point-sampling (Poisson and Random) and two methods of triangle collapse (Quadric and Clustering). Our experimental procedure involves fitting a GMM to a low-resolution mesh (1,000 faces) or point cloud (1,000 points) and evaluating the likelihood of the resulting GMM against a dense point sampling of the original shape (50,000 points). The different sampling strategies can be seen visually in Figure 5.

To generate a low-resolution mesh, we try two methods, clustering decimation [18], which is fast, and quadric error decimation [19], which is more accurate. To generate a low-resolution point cloud, we both randomly sample points on the surface on the mesh, and use Poisson Disc sampling to ensure uniform samples [20]. As before, the first is faster while the latter produces better results. Poisson Disc sampling is also used to generate the high-resolution, ”ground truth”, point cloud used for evaluation.

The results of these experiments are reported in Figure 4. As before, the mesh-based registration strategies are largely invariant to initialization method while the point-based strategies often prefer k-means initialization (exceptions discussed in Sec. IV-C). The best results often came from the use of Poisson Disc Sampling (with k-means), which ensures uniform coverage of the surface areas. In contrast, using random samples generated the lowest quality results. The mesh-based techniques proved to be reliable across all tests, regardless of mixture number and initialization.

IV-C Discussion

All of these experiments were run using the GaussianMixture code base from scikit-learn v0.20.0  [21]. We modified the code to support additional weight and covariance terms (Sec. V-A), which were general enough to allow us to implement all proposed methods.

One surprising result in Figures  3 and 4 was that random sampled points performed better with random initialization than k-means initialization (at high mixture numbers). As the EM algorithm only finds a local minima, this suggests that k-means may not always be an ideal initialization technique. We believe that this occurred due to either a bad local minima from initialization, or fairly flat cost during optimization, leading to an early exit condition being triggered (despite our tight tolerance of $\epsilon=10^{-12}$). This behavior was never observed when using our proposed exact mesh formulation.

We used 25 iterations for all of these experiments. When 100 iterations were used, the point-based methods performed better (nearly as good as our proposed method). However, our method often converged in about $\frac{1}{3}$ the number of iterations, so we picked a lower iteration number for consistency in runtime.

V Extensions

V-A Generalization to other primitives

While equations 8,9,10 were derived specifically for triangles, the update equations can be written more generally for any primitive (triangles, Gaussian mixtures, cuboids, etc.) using $\mu_{p},\Sigma_{p},\alpha_{p}$ to denote primitive’s mean, covariance and size ($\int_{S}dS$) respectively. Then the loss, mean update, and covariance update equations for a Gaussian Mixture can be written with equations 11,12,13. The previous equations can be seen as a special case of these formulas, which provide an M-step update for any set of geometric primitives $p\in\textbf{P}$ in fitting a Gaussian Mixture Model.

[TABLE]

[TABLE]

We note that these are the exact same update equations used in fitting hierarchical Gaussian Mixture Models [22]. However, while previous work applied these equations to fitting GMMs to existing GMMs, we show how this update can be used for fitting geometric data. This general form allows for easy substitution of known structural information (in the form of a second moment) about any geometric primitive. Computing covariance structures for arbitrary polyhedra is well studied area of research [23, 16].

V-B Number of Mixtures

In our later experiments, we fix the number of mixture models. Unless otherwise stated, we use $K=100$. A visual example of this mixture can be see in Figures 5(d) and 5(h). We picked $K=100$ as this matches the experimental conditions recommended for using GMMs for SLAM [8]. In practice, there are many ways to select this number, including flatness of the distribution’s KL-divergence  [24, 25], flatness of the mixture themselves [7], or by evaluating an information criterion [26]. We believe that when this technique is used in practice, this number can either be found through cross-validation [8] on a registration dataset or by using external system information such as depth sensor noise models [27].

VI Applications

The experiments in section IV showed that fitting Gaussian Mixture Models using structural information tends to produce higher quality probability distributions. Some recent work has focused solely on the efficient nature of GMMs in representing shapes [6]. Here we show that our improvements in model quality produce an appreciable performance improvement in actual 3D computer vision applications.

Gaussian Mixture Models have found wide use in the 3D registration literature. From the Normal Distance Transform [30, 31], to variants of the $L_{2}$ loss [32, 33] and even Coherent Point Drift [28], many 3D registration methods utilize Gaussian Mixture Models. Their benefits include robustness to noise, smooth variation over 3D space, speed of evaluation, and straightforward control over model complexity. These models can provide results that are state-of-the-art in both runtime and registration accuracy [7]. We show that applying our proposed mesh GMM fitting can produce an improvement in these results.

VI-A Mesh Registration

We replicate the experimental setup of a recent paper [7], demonstrating how Gaussian Mixture Models can be used for efficient 3D registration. As our experimental setup matches [7], the 20 different dozen registration methods compared in Figure 3 of that work can be directly compared against the results here. Their experiment operates on taking a large number of random deformations of the Stanford Bunny and evaluating the final quality of fitting result. Our results can be seen in Figure 6.

To perform 3D registration, we first build a GMM for the res4 Stanford Bunny, as in Section IV-A, using 100 iterations of EM with a tolerance of $10^{-5}$. We then sample vertex number of 3D points from the surface of the mesh ($N=453$). While previous work has focused on a Point-to-Distribution (P2D) technique with a polynomial approximation of the likelihood function [34], we do straightforward P2D in our experiments. Our registration process consists of finding the rigid body transformation that maximizes equation 1. We use the identity transformation for initialization and then perform gradient-based optimization to find the local minima. We compute gradients using numerical differences. For the optimizer, we tried both Conjugate Gradients [35, 36] and BFGS [37] as optimization strategies, which produced similar results and we report the BFGS results as it ran faster.

To parameterize our rigid-body transformation, we perform the optimization on $\mathbb{R}^{7}$, with a translation $t\in\mathbb{R}^{3}$ and a quaternion $q\in\mathbb{R}^{4}$. Quaternions are well studied in the context of optimization for rigid-body transformation [38, 39]. As we use numerical differences in our optimizer, we did not utilize methods the closed-form gradients for quaternions [40]. While many authors prefer the exponential map for optimizing rigid-body transformations [8], our experiments using the rotation vector $v=\theta\hat{v}\in\mathbb{R}^{3}$ (with rotation angle $\theta$ around the unit vector $\hat{v}$) produced nearly identical results in our final registration result.

We performed these experiments on our mesh and point-derived Gaussian Mixture Models, as well as two baselines. We implemented our own point-to-point Iterative Closest Point (ICP) method [1] and used an exiting implementation of Coherent Point Drift (CPD) [28] from pyCPD [29]. We adjusted pyCPD to run for 150 iterations to approximately match the run-time of our P2D GMM registration. ICP we ran for up to 50,000 iterations, or until the improvement in mean matching error was below $10^{-9}$. For consistency, all methods used in this paper were implemented in the Python programming language and only used the CPU.

VI-B Analysis of Mesh Registration

The results in Figure 6 demonstrate that our mesh-derived Gaussian Mixture Model provides improved registration results when using P2D compared to the existing baselines, ICP and CPD. Not only does our method produce better registration on average, but it also demonstrates a better distribution of errors. Specifically, the small difference between the median and mean errors shows that our method is less prone to outliers. On the other hand, the point-based GMM P2D registration results had outliers that dragged the mean towards the worst quartile of results. The randomly initialized point-based GMM had a mean that was about three times that of its median result, suggesting that some experiments produced results in the incorrect local minima.

VI-C Other 3D Models

We report results on additional 3D models in Table I. For consistency, we decimated each model to 1000 faces using [19] and then repeated our previous experiments exactly (except that the registration results are now the average of 25 runs). The likelihood column reports the per-sample average log-likelihood of ground truth, where larger numbers are better. The translation and rotation errors are reported as a percentage of the error obtained by ICP registration. In all our experiments, the mesh-based GMM always outperformed the point-based one, often significantly.

VI-D Visual Odometry

Our proposed method can be also used for improved models of partial view observations, such as those seen in simultaneous localization and mapping (SLAM). In this case, the geometric primitives used represent not the exact surface, as with our mesh experiment above, but instead an uncertainty region for each 3D measurement.

We performed experiments using a GMM distribution-to-distribution (D2D) registration method for visual odometry [8], reproducing an experiment on an RGB-D dataset sequence from the TUM dataset (freiburg3 long office household) [42]. We are able to incorporate structural information into the fitting of the Gaussian Mixture Models by adding depth uncertainty information around each 3D point and applying equation 13 during GMM fitting. The results are shown in Figure 7.

For our registration experiments, we first subsampled the depth images to $160\times 120$ resolution before performing frame-to-frame registration over the 2510 frame sequence. The ICP method used our aforementioned point-to-point ICP method over 2,500 points randomly sampled from each point cloud (selected to roughly match the run-time performance of our method). The GMM method fit a $K=100$ GMM to the point cloud using our uncertainty model, and performed registration using a D2D metric [8]. We used the determinant-free method as it was much faster in our re-implementation. Additionally, our implementation used numerical gradients and BFGS [37] as the optimizer.

In this case, our primitive model was simple one, a rectangle representing the size of each 3D measurement in the X and Y axes of the camera. This generated a trajectory with absolute translational error RMSE of $0.878m$, a small improvement ($2.4\%$) over the $0.899m$ APE produced from building GMMs without uncertainty primitives. Additionally, we found that D2D registration time was $22\%$ faster when using GMMs built from primitives.

VII Conclusion

We have shown how to build Gaussian Mixture Models by incorporating structural information into their Expectation Maximization algorithm. We demonstrate theoretical and empirical equivalence with traditional techniques, along with providing a fast approximation to our proposed method. By using the covariance structure from the triangles of a mesh, we are able to build GMM models more quickly, robustly and to higher quality. Additionally, these models provide an improved result in 3D registration. For a theoretical understanding of how geometric structures, point samples, and integrals interact, our product integral derivation provides a model that is invariant to resampling (such as triangles being merged or split while retaining the same overall 3D structure). We believe that paper demonstrates that using structural information can lead to methods that are faster, more robust, and more lead to improved performance.