Fitting Cylinder to Point Cloud Data

Paláncz B,  Somogyi Á, Rehány N. and  Lovas T.

Department of Photogrammetry and Geoinformatics,
Budapest University of Technology and Economy
1521 Budapest, Hungary

Abstract

A new robust parameter estimation method taking into account the real model error distribution is presented for large size of noisy data points. Maximum likelihood technique is employed to compute the model parameters assuming that the distribution of the model errors is a Gaussian mixture corresponding to the inlier and outlier data points. The maximization is carried out by local method, where the initial guess values were computed via numerical Groebner basis. After parameter estimation based on an initially computed distribution, the real errors are determined and the corresponding Gaussian mixture is identified via expectation maximization algorithm. The iteration procedure is converging when the error distribution becomes stationary. The method is illustrated via identifying  tree stems using ground-based Lidar data.

Keywords

parameter estimation, cylinder, point cloud, Gröbner basis, outliers, maximum likelihood, Gaussian mixture, expectation maximization, Lidar.

   1 - Introduction

Fitting real circular cylinders is an important problem in the representation of geometry of man made structures such as industrial plants Vosselman et al (2004), deformation analysis of tunnel, Stal et al (2012), detecting and monitoring the deformation in deposition holes, Carrea et al (2014), positioning of femur pieces for surgical fracture reduction, Winkelbach et al (2003), estimating tree stems, Khameneh (2013) and so on. Since planes and cylinders compose up to 85% of all objects in industrial scenes, research in 3D reconstruction and modeling - see CAD - CAM applications - have largely focused on these two important geometric primitives, i.e. Petitjean (2002).

In general, to define a cylinder we need 5 parameters: 4 for the axis and 1 for the radius, so one requires at last 5 points to determine parameters. It goes without saying, that in special cases, like cylinder is parallel to an axis or to a plane, less parameter is enough, see Beder and Förstner (2006).

In case of more than 5 points, one has to find the 5 parameters of the cylinder, that the sum of the distances of data points from the cylinder surface is minimum in least square sense. Basically,  two approaches can be followed:

● find the direction vector of the cylinder center-axis and transform the data points into vertical position to the x -y plane, then the remained 3 parameters of the cylinder oriented parallel to z axis - 2 shifting parameters and the radius - can be computed, Beder and Förstner (2006).

● all of the 5 parameters are computed via least square method employing local optimization technique, i.e. Lukacs et al (1998), Lichtblau (2007).

There are  different methods to find the orientation of the cylinder center axis:

◦ one may use Principal Component Analysis (PCA),

◦ considering the general form of a cylinder as a second order surfaces, the direction vector can be partially extracted, and computed via linear least square, see Khameneh (2013),

◦ the PCA method can be modified such as that instead of single points, a local neighborhood of  randomly selected points are employed, Ruiz et al (2013).

◦ employing Hough transformation, see Su and Bethel (2010), Rabbani and van den Heuvel (2005).

All of these methods consider least square method assuming that model error  has normal Gaussian distribution with zero mean value. In this study we consider realistic model error distribution applying maximization likelihood technique for parameter estimation, where this distribution is represented by a Gaussian mixture of the in and outliers error, and identified by expectation maximization algorithm. For geometric modeling of a general cylinder, a vector algebraic approach is followed, see Lichtblau (2007).

   2 - Vector algebraic definition

Given a cylinder with axis line L in . We may parametrize it using five model parameters (a, b, c, d, r), where r stands for the radius of the cylinder. Let  P (x, y, z) is a point of the cylinder and the vector of the axis line of  L is vec = {1, a, c}. Let us consider L’ to pass through the origin of the coordinate system and parallel to L. The translation vector is offset = {0, b, d}. The result of this translation is P’ → P and L’ → L. We project this P’ point onto L’ and denote the length of the orthogonal projection perp. It is computed as follows.

The vector of P’ is projected onto L' and projection will be subtract from the vector of P’. Then for the magnitude of perp, for the norm of it, one can write

Let us apply this definition to set the equation of the cylinder. Considering a point

The vector of the locus of this point on the cylinder having model parameters  (a, b, c, d, r) can be computed as it follows, see Fig. 1.

The vector of the axis line

The vector translating the axis line to pass the origin is,

The function computing the orthogonal projection perp, carries out projection onto the translated axis line L’ and subtracts it from the vector of locus P’ is

Applying it to point P

Fig. 1 Explanation of  perp function

Clearing denominators and formulating of the equation for perp

This is practically the implicit equation of the cylinder with  model parameters  (a, b, c, d, r).

It is important to realized that creating this equation, the algebraic error definition,   has been used.

  3 - Parametrized form of the cylinder equation

In order to employ these model parameters to define the parametric equation of the cylinder, let us developed the usual parametric equation in form of x = x (u, v), y = y (u,v) and z = z (u,v) with actual scaling parameters (u, v) and using the model parameters  (a, b, c, d, r),

The locus of point P is obtained as sum of a vector on L’ plus a vector of length r perpendicular to L. Let v is a parameter along the length of axis L’. Then the projection of  P (u, v) on L  is a vector

All vector perpendicular to L are spanned by any independent pair. We can obtain an orthonormal pair {,} in the standard way by finding the null space to the matrix whose one row is the vector along the axial direction, that is vec, and then using Gram -Schmidt to orthogonalize that pair. Using parameter u, this vector having length r and perpendicular to L can be written as              

Then the locus vector of a general point of the cylinder is, see Fig. 2,

Fig.2 Explanation of general locus vector λ

Let us carry out this computation step by step in symbolic way employing Mathematica,

Then the parametric equation of a general circular cylinder using {a, b, c, d, r) parameters

where u and v the actual scaling parameters.

4 - Implicit equation from the parametric one

One can develop the implicit equation from the parametric one, too.The problem practically means to eliminate u and v or avoiding trigonometrical expression, sin(u), cos(u) and v from  algebraic system symbolically,

and

In order to carry out this computation in Mathematica, from practical reason let sin(u) = U and cos(u)= V. Then our system is

Let us clear denominators

Applying Groebner basis to eliminate U and V and v, the implicit form is

This is the same equation, which was computed in Section 1.

  5 - Computing model parameters in determined case

Let us consider five triples of points

Substituting these points into the implicit form and setting = rsqr, we get

This is a polynomial system for the parameters based on the algebraic error definition. In order to create the determined system for the geometric error model, let us consider geometric error model , namely {| perp|}  -r ,

which can be applied to the points5 data, then

  6 - Computing model parameters in overdetermined case

Let us consider algebraic error model first. This means to minimize the residual of the implicit form

where the algebraic error is, see implicit or vector,

while the geometric error is,

In order to carry out this minimization problem via local minimization, since is much faster than the global one, we may solve the determined system ( n = 5) for randomly selected data points to get initial guess values.

  7 - Application to estimation of tree stem diameter

Outdoor laser scanning measurements have been carried out in the backyard of Budapest University of Technology and Economics, see Fig.3/a. In order to get simple fitting problem rather than segmentation one, the test object, the lower part of the stem of a tree was preselected by segmentation, see Fig. 3/ b

Fig.3/a/b Test environment (to left) and test object (to right)

The experiment has been carried out with a Faro Focus 3D terrestrial laser scanner, see Fig.4.

Fig. 4 Faro Focus 3D scanner

The scanning parameters were set to ½ resolution that equals to 3mm/10m point spacing.  The test data set was cropped from the point cloud; moreover, further resampling was applied in order to reduce the data size. The final data set is composed of 16 434 points in ASCII format, and only the x, y, z coordinates were kept (no intensity values). Let us load the measured data:

Fig.5 The  points of cloud of data

We shall compute the parameters from the geometric error model employing local minimization of the sum of the square of the residual. To do that we compute the initial guess values from 5 randomly chosen points of data using algebraic error model, since the resulted polynomial system can be solved easily via numerical Groebner basis. In deterministic case the number of real solution can even number (0, 2, 4 , or 6). According to our numerical experiences, to ensure reliable initial values, one need such a 5 points which provides at least 4 real solutions.

The equations are

Employing integer coefficients in order to avoid round-off error

Solving the system via numerical Gröbner basis

The real solutions

We create the objective function with the geometric error model,

Employing all of the points

Applying the geometric error model

Then the objective function

From these solution we can select the very one, which gives the smallest residual

This solution will be used as initial guess values for the local minimization,

The result is,

Fig. 6  The fitted cylinder via geometrical error model

The distribution of the model error

Fig.7 Histogram of the model error of the standard approach assumung N (0, σ)

It is clear from Fig .7, that the assumption for the model error is not true. Therefore one should employ likelihood function for the parameter estimation. The real model error distribution as well as the model parameters can be computed  in iterative way. We can consider the distribution shown in Fig. 7. as a first guess for this iteration. In order to handle an empirical distribution in the likelihood function there are at least three ways:

◦ Gaussian kernel functions can be applied to the values of different bins,

◦ Employing a known type of distribution to approximate the empirical histogram,

◦ Considering the empirical distribution as a Gaussian mixture of errors corresponding to inliers and ouliers

Here, we considered the third approach, and employed expectation maximization algorithm to compute the parameters of the component distributions.

8 - Expectation Maximization

Let us consider a two-component Gaussian mixture represented by  the mixture model in the following form,

where

and ’s are the membership weights constrained with

We are looking for the parameter values ) and ). The log-likelihood function in case of  N  samples is

where θ = (,) the parameters of the normal densities. The problem is the direct maximization of  this function, because of the sum of terms inside the logarithm. In order to solve this problem let us introduce the following alternative log-likelihood function

Here ’s are considered as unobserved latent variables taking values 0 or 1. If belongs to the first component then = 0, so

otherwise   belongs to the second component then =1, therefore

where and are the number of the elements of the mixture, which belong to the first and to the second component,respectively.

Since the values of the ’s are actually unknown, we proceed in an iterative fashion, substituting for each its expected value,

This expression is also called as the responsibility of component two for observation i. Then the procedure called the EM algorithm for two -component Gaussian mixture is the following:

Take initial guess for the parameters: θ = (,) and for

Expectation Step: compute the responsibilities:

Maximization Step: compute the weighted means and variances for the two components:

and the mixing probability

Iterate these steps until convergence.

This algorithm is implemented for Mathematica see, Fox et al (2013) as a Mathematica Demonstration project. The code has been modified and applied here.  In the next section we illustrate how this function works.

First we should compute the parameters of the two Gaussian via EM algorithm. Considering the histogram of the model errors, Fig. 7, the guess values of the two Gaussian are

The result of the parameter values

These can be displayed in a table form

Table 1.
Parameters of the Gaussian mixture after zero iteration

Fig. 8 shows the density functions of the two component.

Fig 8. The PDF of the two components

Fig 9. The joint PDF of the mixture

Fig 10. The PDF of the two components and the normalized histogram of the data

The following figures show the convergence of the different parameters

Fig 11. The convergence of the means of the two components

Fig 12. The convergence of the standard deviations

Fig 13. The convergence of the membership weights

Now, we separate the mixture of the samples into two clusters: cluster of outliers and cluster of inliers. Membership values of the first cluster

Membership values of the second cluster

In order to get Boolean (crisp) clustering let us round the membership values

The elements in the first cluster (outliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S1)

The number of these elements (outliers)

Similarly, the elements in the second cluster (inliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S2)

and

which is the number of the inliers. Let us display the ouliers and the inliers

Fig 14. The inliers (blue) and outliers (red) as a first guess resulted standard minimization of the total residual

9 - Maximum likelihood for Gaussian mixture

The likelihood function for a two-component Gaussian mixture can be written as,

where the likelihood function for one of the components can be developed as it follows.

We have seen that the geometric error is

The probability density function for a single Gaussian

Let us substitute the expression of the geometric error, (x →)

We apply a Mathematica function developed by Rose and Smith (2000) which is entitled as SuperLog. This function utilizes pattern-matching code that enhances Mathematica’s ability to simplify expressions involving the natural logarithm of a product of algebraic terms.

Then the likelihood function is

Considering its logarithm which should be maximized

LogL = Log[L]

It is a rather complicated expression, but fortunately Mathematica has built in function to compute this expression numerically.

Now we can create the corresponding likelihood functions for the identified ouliers and inliers.

Then the residuals should be computed,

The likelihood functions can be developed by Mathematica,

Then the likelihood function for the mixture, which should be maximized

In order to carry out local maximization, the result of the standard parameter estimation (see above) can be considered

The result

Now, let us compute the model error distribution

Fig.15 Histogram of the model error employing maximum likelihood method

It means, that employing the model error distribution represented Fig.7, and employing maximum likelihood method for estimating the parameters, with these parameters we get model error distribution represented by Fig.15. If the two distribution is the same, than the parameter estimation is correct. To answer this question let us employ EM for this later distribution.

The result of the parameter values

These can be displayed in a table form

Table 2.
Parameters of the Gaussian mixture after zero iteration

Fig. 2 shows the density functions of the two component.

Fig 16. The PDF of the two components

Fig 17. The joint PDF of the mixture

Fig 18. The PDF of the two components and the normalized histogram of the data

Fig 19. The convergence of the means of the two components

Fig 20. The convergence of the standard deviations

Fig 21. The convergence of the membership weights

Now, we separate the mixture of the samples into two clusters: cluster of outliers and cluster of inliers. Membership values of the first cluster

Membership values of the second cluster

In order to get Boolean (crisp) clustering let us round the membership values

The elements in the first cluster (outliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S1)

The number of these elements (outliers)

Similarly, the elements in the second cluster (inliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S2)

and

which is the number of the inliers. Let us display the ouliers and the inliers

Fig 22. The inliers (blue) and outliers (red) as a first guess resulted standard minimization of the total residual

Since the two distribution is different, we need to carry out a new iteration step.

10 - Iteration process

Second Iteration

The result of the parameter values

These can be displayed in a table form

Now, we separate the mixture of the samples into two clusters: cluster of outliers and cluster of inliers. Membership values of the first cluster

Membership values of the second cluster

In order to get Boolean (crisp) clustering let us round the membership values

The elements in the first cluster (outliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S1)

The number of these elements (outliers)

Similarly, the elements in the second cluster (inliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S2)

and

which is the number of the inliers. Let us display the ouliers and the inliers

Third Iteration

The result of the parameter values

These can be displayed in a table form

Now, we separate the mixture of the samples into two clusters: cluster of outliers and cluster of inliers. Membership values of the first cluster

Membership values of the second cluster

In order to get Boolean (crisp) clustering let us round the membership values

The elements in the first cluster (outliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S1)

The number of these elements (outliers)

Similarly, the elements in the second cluster (inliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S2)

and

which is the number of the inliers. Let us display the ouliers and the inliers

We consider it as the last iteration step. The parametric equation of the fitted cylinder is

Let us visualize it with the inliers and outliers data points

Fig 23. The inliers (blue) and outliers (red) as a first guess resulted by the maximum likelihood technique

The results of  3 iteration steps can be seen in Table 3 and 4.

Table 3

The computed cylinder parameters

Table 4

Parameters of the components of the Gaussian mixture

Iteration
0	0.0813	-0.010	0.0536	0.0227	0.1037	0.8962	1318	15116
1	0.0826	-0.011	0.0597	0.0183	0.1188	0.8811	1577	14857
2	0.0814	-0.012	0.0631	0.0176	0.1250	0.8749	1641	14793
3	0.0837	-0.012	0.0627	0.0176	0.1229	0.8771	1637	14797

11 - Application to leafy tree

In the example above there were only  ∼10 % outliers.  Now let us see another example where the ratio of the outliers are much more higher, nearly 40% . This situation is closer to segmentation than to simple fitting problem. Here the test object is the same tree, but with foliage, see Fig.

Fig 24. The inliers (blue) and outliers (red) as a first guess resulted standard minimization of the total residual

Let us load the data

Fig. 24 The points of cloud of data

The first guess will be the result, which we have got without foliage.

Fig. 25 The first estimation of the cylinder fitting to leafy tree

Applying the geometric error model, the model error distribution of this first approach

Employing all of the points

Fig. 26 The model error distribution of the first estimation

Now let us identify the components of the Gaussian mixture

The result of the parameter values

These can be displayed in a table form

Table 5.
Parameters of the Gaussian mixture after zero iteration

The density functions of the two component.

Fig 27. The PDF of the two components

Fig 28. The joint PDF of the mixture

Fig 29. The PDF of the two components and the normalized histogram of the data

Fig 30. The convergence of the means of the two components

Fig 31. The convergence of the standard deviations

Fig 32. The convergence of the membership weights

Now, we separate the mixture of the samples into two clusters: cluster of outliers and cluster of inliers. Membership values of the first cluster

Membership values of the second cluster

In order to get Boolean (crisp) clustering let us round the membership values

The elements in the first cluster (outliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S1)

The number of these elements (outliers)

Similarly, the elements in the second cluster (inliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S2)

and

which is the number of the inliers. Let us display the ouliers and the inliers

Fig. 33 The inliers (blue) and outliers (red) as a first guess resulted by the first approach

Now let us compute the first iteration employing maximum likelihood method

Fig 34. The PDF of the two components and the normalized histogram of the data

Now the parameters of the Gaussian mixture

The result of the parameter values

These can be displayed in a table form

Table 6.
Parameters of the Gaussian mixture after zero iteration

The density functions of the two component.

Fig 35. The PDF of the two components

Fig 36. The joint PDF of the mixture

Fig 37 The PDF of the two components and the normalized histogram of the data

Fig 38. The convergence of the means of the two components

Fig 39. The convergence of the standard deviations

Fig 40. The convergence of the membership weights

Now, we separate the mixture of the samples into two clusters: cluster of outliers and cluster of inliers. Membership values of the first cluster

Membership values of the second cluster

In order to get Boolean (crisp) clustering let us round the membership values

The elements in the first cluster (outliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S1)

The number of these elements (outliers)

Similarly, the elements in the second cluster (inliers) are the corresponding elements of those having value 1 in the set containing the rounded member values (S2)

and

which is the number of the inliers. Let us display the ouliers and the inliers

Fig. 34 The inliers (blue) and outliers (red) as a first guess resulted by the first iteration

The inlier points,

Fig. 35 The inliers

We consider it as the last iteration step. The parametric equation of the fitted cylinder is

Fig. 36 The inliers (blue) and the fitted cylinder

Fig. 37 The inliers (blue) and outliers (red) and the fitted cylinder

Table 7

The computed cylinder parameters

Table 8

Parameters of the components of the Gaussian mixture

Iteration
0	0.3412	-0.011	0.2759	0.0205	0.4640	0.5360	39982	51107
1	0.3435	-0.010	0.2751	0.0208	0.4588	0.5418	39404	51685

Conclusions

The results our computation illustrate that in case of noisy measurement, the frequently assumed N (0, σ) model error distribution is not true. Therefore to estimate the model parameters, iteration is necessary with the proper representation of the non-Gaussian distribution. In this study expectation maximization algorithm was employed. In case of general cylinder geometry, the 5 model parameters can be computed via local maximization of the maximum likelihood function. The implicit equation based on algeraic error definition, and solved easily by Groebner basis, can provide a proper initil guess value for the maximization problem.

References

Beder C and Förstner W (2006) Direct Solutions for Computing Cylinders from Minimal Sets of 3D Points, Computer Vision-ECCV 2006, Lecture Notes in Computer Sci. Vol. 3951 pp. 135-146

Carrea D, Jaboyedoff M and Derron M-H (2014) Feasibility Study of Point Cloud Data from Test Deposition Holes for Deformation Analysis, Working Report 2014-01, Universite de Lausanne (UNIL)

Fox A, Smith J, Ford R and Doe J (2013) Expectation Maximization for Gaussian Mixture Distributions, Wolfram Demonstration Project.

Khameneh M (2013) Tree Detection and Species Identification using LiDAR Data, MSc. Thesis, KTH, Royal Institute of Technology, Stockholm

Lichtblau D (2007) Cylinders Through Five Points: Computational Algebra and Geometry. Automated Deduction in Geometry, Lecture Notes in Computer Sci. Vol. 4869, pp 80-97

Lukacs G, Martin R and Marshall D (1998) Faithful Least-Squares Fitting of Spheres, Cylinders, Cones and Tori for Reliable Segmentation (1998) Burkhardt H, Neumann B (eds) Computer Vision -ECCV’98 Vol.I.LNCS 1406, pp.671-686, Springer-Verlag.

Petitjean S (2002) A survey of methods for recovering quadrics in triangle meshes, ACM Computing Surveys, Vol. 34, No.2, pp. 211-262.

Ruiz O, Arroyave S and Acosta D (2013) Fitting of Analytic Surfaces to Noisy Point Clouds, American J. of Computational Mathematics, 3, pp 18-26

Rose C and Smith D (2000) Symbolic Maximum Likelihood Estimation with Mathematica, The Statistician 49, pp.229 - 240

Stal C, Timothy Nuttens T,Constales D, Schotte K, De Backer H, De Wulf A (2012) Automatic Filtering of Terrestrial Laser Scanner Data from Cylindrical Tunnels, TS09D - Laser Scanners III, 5812, FIG Working Week 2012, Knowing to manage the territory, protect the environment, evaluate the cultural heritage, Rome, Italy, 6-10 May 2012

Su Y-T and Bethel J (2010) Detection and Robust Estimation of Cylinder Features in Point Clouds, ASPRS 2010 Annual Conference, San Diego, California April 26-30.

Vosselman G, Gorte B, Sithole G and Rabbani T (2004) Recognizing structure in laser scanner point clouds. In Int. Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol. 46. pp. 33-38.

Winkelbach S, Westphal R and Goesling T (2003) Pose estimation of cylinder fragments for semi-automatic bone fracture reduction. In Pattern Recognition (DAGM 2003). Michaelis B and Krell G eds. Lecture Notes in Computer Science 2781: pp. 566 - 573, Springer- Verlag 2003