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 Fitting_Cylinder_1.png. 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

Fitting_Cylinder_2.png

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

Fitting_Cylinder_3.png

Fitting_Cylinder_4.png

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

Fitting_Cylinder_5.png

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

Fitting_Cylinder_6.png

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

Fitting_Cylinder_7.png

Applying it to point P

Fitting_Cylinder_8.png

Fitting_Cylinder_9.png

Fitting_Cylinder_10.gif

Fig. 1 Explanation of  perp function

Clearing denominators and formulating of the equation for perp

Fitting_Cylinder_11.png

Fitting_Cylinder_12.png

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, Fitting_Cylinder_13.png  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

Fitting_Cylinder_14.png

All vector perpendicular to L are spanned by any independent pair. We can obtain an orthonormal pair {Fitting_Cylinder_15.png,Fitting_Cylinder_16.png} 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              

Fitting_Cylinder_17.png

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

Fitting_Cylinder_18.png

Fitting_Cylinder_19.gif

Fitting_Cylinder_20.png

Fig.2 Explanation of general locus vector λ

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

Fitting_Cylinder_21.png

Fitting_Cylinder_22.png

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

Fitting_Cylinder_23.png

Fitting_Cylinder_24.png

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,

Fitting_Cylinder_25.png

Fitting_Cylinder_26.png

Fitting_Cylinder_27.png

and

Fitting_Cylinder_28.png

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

Fitting_Cylinder_29.png

Fitting_Cylinder_30.png

Let us clear denominators

Fitting_Cylinder_31.png

Fitting_Cylinder_32.png

Fitting_Cylinder_33.png

Fitting_Cylinder_34.png

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

Fitting_Cylinder_35.png

Fitting_Cylinder_36.png

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

Fitting_Cylinder_37.png

Fitting_Cylinder_38.png

  5 - Computing model parameters in determined case

Let us consider five triples of points

Fitting_Cylinder_39.png

Fitting_Cylinder_40.png

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

Fitting_Cylinder_42.png

Fitting_Cylinder_43.png

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 ,

Fitting_Cylinder_44.png

which can be applied to the points5 data, then

Fitting_Cylinder_45.png

Fitting_Cylinder_46.png

  6 - Computing model parameters in overdetermined case

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

Fitting_Cylinder_47.png

where the algebraic error is, see implicit or vector,

Fitting_Cylinder_48.png

while the geometric error is,

Fitting_Cylinder_49.png

Fitting_Cylinder_50.png

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

Fitting_Cylinder_51.gif

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.

Fitting_Cylinder_52.gif

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:

Fitting_Cylinder_53.png

Fitting_Cylinder_54.png

Fitting_Cylinder_55.png

Fitting_Cylinder_56.png

Fitting_Cylinder_57.gif

Fig.5 The  points of cloud of data

Fitting_Cylinder_58.png

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.

Fitting_Cylinder_59.png

Fitting_Cylinder_60.png

Fitting_Cylinder_61.png

The equations are

Fitting_Cylinder_62.png

Fitting_Cylinder_63.png

Employing integer coefficients in order to avoid round-off error

Fitting_Cylinder_64.png

Fitting_Cylinder_65.png

Solving the system via numerical Gröbner basis

Fitting_Cylinder_66.png

Fitting_Cylinder_67.png

The real solutions

Fitting_Cylinder_68.png

Fitting_Cylinder_69.png

We create the objective function with the geometric error model,

Fitting_Cylinder_70.png

Employing all of the points

Fitting_Cylinder_71.png

Fitting_Cylinder_72.png

Applying the geometric error model

Fitting_Cylinder_73.png

Then the objective function

Fitting_Cylinder_74.png

Fitting_Cylinder_75.png

Fitting_Cylinder_76.png

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

Fitting_Cylinder_77.png

Fitting_Cylinder_78.png

Fitting_Cylinder_79.png

Fitting_Cylinder_80.png

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

Fitting_Cylinder_81.png

Fitting_Cylinder_82.png

Fitting_Cylinder_83.png

Fitting_Cylinder_84.png

The result is,

Fitting_Cylinder_85.png

Fitting_Cylinder_86.png

Fitting_Cylinder_87.png

Fitting_Cylinder_88.png

Fitting_Cylinder_89.gif

Fig. 6  The fitted cylinder via geometrical error model

The distribution of the model error

Fitting_Cylinder_90.png

Fitting_Cylinder_91.png

Fitting_Cylinder_92.gif

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

Fitting_Cylinder_93.png

Fitting_Cylinder_94.png

Fitting_Cylinder_95.png

Fitting_Cylinder_96.png

Fitting_Cylinder_97.png

Fitting_Cylinder_98.png

Fitting_Cylinder_99.png

Fitting_Cylinder_100.png

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,

Fitting_Cylinder_101.png

where

Fitting_Cylinder_102.png

and Fitting_Cylinder_103.png’s are the membership weights constrained with

Fitting_Cylinder_104.png

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

Fitting_Cylinder_107.png

where θ = (Fitting_Cylinder_108.png,Fitting_Cylinder_109.png) 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

Fitting_Cylinder_110.png

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

Fitting_Cylinder_114.png

otherwise  Fitting_Cylinder_115.png belongs to the second component then Fitting_Cylinder_116.png=1, therefore

Fitting_Cylinder_117.png

where Fitting_Cylinder_118.png and Fitting_Cylinder_119.png 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 Fitting_Cylinder_120.png’s are actually unknown, we proceed in an iterative fashion, substituting for each Fitting_Cylinder_121.png its expected value,

Fitting_Cylinder_122.png

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: θ = (Fitting_Cylinder_123.png,Fitting_Cylinder_124.png) and for Fitting_Cylinder_125.png

Expectation Step: compute the responsibilities:

Fitting_Cylinder_126.png

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

Fitting_Cylinder_127.png

Fitting_Cylinder_128.png

Fitting_Cylinder_129.png

Fitting_Cylinder_130.png

and the mixing probability

Fitting_Cylinder_131.png

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.

Fitting_Cylinder_132.png

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

Fitting_Cylinder_133.png

The result of the parameter values

Fitting_Cylinder_134.png

Fitting_Cylinder_135.png

These can be displayed in a table form

Fitting_Cylinder_136.png

Fitting_Cylinder_137.png

Table 1.
Parameters of the Gaussian mixture after zero iteration

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

Fitting_Cylinder_138.png

Fitting_Cylinder_139.gif

Fig 8. The PDF of the two components

Fitting_Cylinder_140.png

Fitting_Cylinder_141.gif

Fig 9. The joint PDF of the mixture

Fitting_Cylinder_142.png

Fitting_Cylinder_143.gif

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

Fitting_Cylinder_144.png

Fitting_Cylinder_145.png

Fitting_Cylinder_146.gif

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

Fitting_Cylinder_147.png

Fitting_Cylinder_148.gif

Fig 12. The convergence of the standard deviations

Fitting_Cylinder_149.png

Fitting_Cylinder_150.gif

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

Fitting_Cylinder_151.png

Fitting_Cylinder_152.png

Membership values of the second cluster

Fitting_Cylinder_153.png

Fitting_Cylinder_154.png

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

Fitting_Cylinder_155.png

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)

Fitting_Cylinder_156.png

The number of these elements (outliers)

Fitting_Cylinder_157.png

Fitting_Cylinder_158.png

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)

Fitting_Cylinder_159.png

and

Fitting_Cylinder_160.png

Fitting_Cylinder_161.png

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

Fitting_Cylinder_162.png

Fitting_Cylinder_163.png

Fitting_Cylinder_164.png

Fitting_Cylinder_165.gif

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,

Fitting_Cylinder_166.png

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

We have seen that the geometric error is

Fitting_Cylinder_167.png

The probability density function for a single Gaussian

Fitting_Cylinder_168.png

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

Fitting_Cylinder_170.png

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

Fitting_Cylinder_171.png

Fitting_Cylinder_172.png

Considering its logarithm which should be maximized

LogL = Log[L]

Fitting_Cylinder_173.png

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.

Fitting_Cylinder_174.png

Then the residuals should be computed,

Fitting_Cylinder_175.png

Fitting_Cylinder_176.png

Fitting_Cylinder_177.png

Fitting_Cylinder_178.png

The likelihood functions can be developed by Mathematica,

Fitting_Cylinder_179.png

Fitting_Cylinder_180.png

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

Fitting_Cylinder_181.png

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

Fitting_Cylinder_182.png

Fitting_Cylinder_183.png

Fitting_Cylinder_184.png

Fitting_Cylinder_185.png

Fitting_Cylinder_186.png

Fitting_Cylinder_187.png

The result

Fitting_Cylinder_188.png

Fitting_Cylinder_189.png

Now, let us compute the model error distribution

Fitting_Cylinder_190.png

Fitting_Cylinder_191.png

Fitting_Cylinder_192.png

Fitting_Cylinder_193.png

Fitting_Cylinder_194.gif

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.

Fitting_Cylinder_195.png

The result of the parameter values

Fitting_Cylinder_196.png

Fitting_Cylinder_197.png

These can be displayed in a table form

Fitting_Cylinder_198.png

Fitting_Cylinder_199.png

Table 2.
Parameters of the Gaussian mixture after zero iteration

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

Fitting_Cylinder_200.png

Fitting_Cylinder_201.gif

Fig 16. The PDF of the two components

Fitting_Cylinder_202.png

Fitting_Cylinder_203.gif

Fig 17. The joint PDF of the mixture

Fitting_Cylinder_204.png

Fitting_Cylinder_205.gif

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

Fitting_Cylinder_206.png

Fitting_Cylinder_207.png

Fitting_Cylinder_208.gif

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

Fitting_Cylinder_209.png

Fitting_Cylinder_210.gif

Fig 20. The convergence of the standard deviations

Fitting_Cylinder_211.png

Fitting_Cylinder_212.gif

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

Fitting_Cylinder_213.png

Fitting_Cylinder_214.png

Membership values of the second cluster

Fitting_Cylinder_215.png

Fitting_Cylinder_216.png

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

Fitting_Cylinder_217.png

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)

Fitting_Cylinder_218.png

The number of these elements (outliers)

Fitting_Cylinder_219.png

Fitting_Cylinder_220.png

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)

Fitting_Cylinder_221.png

and

Fitting_Cylinder_222.png

Fitting_Cylinder_223.png

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

Fitting_Cylinder_224.png

Fitting_Cylinder_225.png

Fitting_Cylinder_226.png

Fitting_Cylinder_227.gif

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

Fitting_Cylinder_228.png

Fitting_Cylinder_229.png

Fitting_Cylinder_230.png

Fitting_Cylinder_231.png

Fitting_Cylinder_232.png

Fitting_Cylinder_233.png

Fitting_Cylinder_234.png

Fitting_Cylinder_235.png

Fitting_Cylinder_236.png

Fitting_Cylinder_237.png

Fitting_Cylinder_238.png

Fitting_Cylinder_239.png

Fitting_Cylinder_240.png

Fitting_Cylinder_241.png

Fitting_Cylinder_242.png

Fitting_Cylinder_243.png

Fitting_Cylinder_244.png

Fitting_Cylinder_245.png

Fitting_Cylinder_246.png

Fitting_Cylinder_247.png

Fitting_Cylinder_248.gif

Fitting_Cylinder_249.png

The result of the parameter values

Fitting_Cylinder_250.png

Fitting_Cylinder_251.png

These can be displayed in a table form

Fitting_Cylinder_252.png

Fitting_Cylinder_253.png

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

Fitting_Cylinder_254.png

Fitting_Cylinder_255.png

Membership values of the second cluster

Fitting_Cylinder_256.png

Fitting_Cylinder_257.png

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

Fitting_Cylinder_258.png

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)

Fitting_Cylinder_259.png

The number of these elements (outliers)

Fitting_Cylinder_260.png

Fitting_Cylinder_261.png

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)

Fitting_Cylinder_262.png

and

Fitting_Cylinder_263.png

Fitting_Cylinder_264.png

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

Fitting_Cylinder_265.png

Fitting_Cylinder_266.png

Fitting_Cylinder_267.png

Fitting_Cylinder_268.gif

Third Iteration

Fitting_Cylinder_269.png

Fitting_Cylinder_270.png

Fitting_Cylinder_271.png

Fitting_Cylinder_272.png

Fitting_Cylinder_273.png

Fitting_Cylinder_274.png

Fitting_Cylinder_275.png

Fitting_Cylinder_276.png

Fitting_Cylinder_277.png

Fitting_Cylinder_278.png

Fitting_Cylinder_279.png

Fitting_Cylinder_280.png

Fitting_Cylinder_281.png

Fitting_Cylinder_282.png

Fitting_Cylinder_283.png

Fitting_Cylinder_284.png

Fitting_Cylinder_285.png

Fitting_Cylinder_286.png

Fitting_Cylinder_287.png

Fitting_Cylinder_288.png

Fitting_Cylinder_289.gif

Fitting_Cylinder_290.png

The result of the parameter values

Fitting_Cylinder_291.png

Fitting_Cylinder_292.png

These can be displayed in a table form

Fitting_Cylinder_293.png

Fitting_Cylinder_294.png

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

Fitting_Cylinder_295.png

Fitting_Cylinder_296.png

Membership values of the second cluster

Fitting_Cylinder_297.png

Fitting_Cylinder_298.png

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

Fitting_Cylinder_299.png

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)

Fitting_Cylinder_300.png

The number of these elements (outliers)

Fitting_Cylinder_301.png

Fitting_Cylinder_302.png

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)

Fitting_Cylinder_303.png

and

Fitting_Cylinder_304.png

Fitting_Cylinder_305.png

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

Fitting_Cylinder_306.png

Fitting_Cylinder_307.png

Fitting_Cylinder_308.png

Fitting_Cylinder_309.gif

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

Fitting_Cylinder_310.png

Fitting_Cylinder_311.png

Let us visualize it with the inliers and outliers data points

Fitting_Cylinder_312.png

Fitting_Cylinder_313.png

Fitting_Cylinder_314.gif

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

Iteration a b c d r
0 -0.6043 -11.3773 -3.6870 146.479 0.1788
1 -0.6438 -11.3064 -3.6080 146.386 0.1752
2 -0.6556 -11.2872 -3.5783 146.343 0.1749
3 -0.6619 -11.2767 -3.5779 146.344 0.1750

Table 4

Parameters of the components of the Gaussian mixture

Iteration Fitting_Cylinder_315.png Fitting_Cylinder_316.png Fitting_Cylinder_317.png Fitting_Cylinder_318.png Fitting_Cylinder_319.png Fitting_Cylinder_320.png Fitting_Cylinder_321.png Fitting_Cylinder_322.png
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.

Fitting_Cylinder_323.gif

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

Let us load the data

Fitting_Cylinder_324.png

Fitting_Cylinder_325.png

Fitting_Cylinder_326.png

Fitting_Cylinder_327.png

Fitting_Cylinder_328.gif

Fig. 24 The points of cloud of data

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

Fitting_Cylinder_329.png

Fitting_Cylinder_330.png

Fitting_Cylinder_331.png

Fitting_Cylinder_332.png

Fitting_Cylinder_333.png

Fitting_Cylinder_334.gif

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

Fitting_Cylinder_335.png

Employing all of the points

Fitting_Cylinder_336.png

Fitting_Cylinder_337.png

Fitting_Cylinder_338.png

Fitting_Cylinder_339.png

Fitting_Cylinder_340.png

Fitting_Cylinder_341.gif

Fig. 26 The model error distribution of the first estimation

Now let us identify the components of the Gaussian mixture

Fitting_Cylinder_342.png

The result of the parameter values

Fitting_Cylinder_343.png

Fitting_Cylinder_344.png

These can be displayed in a table form

Fitting_Cylinder_345.png

Fitting_Cylinder_346.png

Table 5.
Parameters of the Gaussian mixture after zero iteration

The density functions of the two component.

Fitting_Cylinder_347.png

Fitting_Cylinder_348.gif

Fig 27. The PDF of the two components

Fitting_Cylinder_349.png

Fitting_Cylinder_350.gif

Fig 28. The joint PDF of the mixture

Fitting_Cylinder_351.png

Fitting_Cylinder_352.gif

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

Fitting_Cylinder_353.png

Fitting_Cylinder_354.png

Fitting_Cylinder_355.gif

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

Fitting_Cylinder_356.png

Fitting_Cylinder_357.gif

Fig 31. The convergence of the standard deviations

Fitting_Cylinder_358.png

Fitting_Cylinder_359.gif

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

Fitting_Cylinder_360.png

Fitting_Cylinder_361.png

Membership values of the second cluster

Fitting_Cylinder_362.png

Fitting_Cylinder_363.png

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

Fitting_Cylinder_364.png

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)

Fitting_Cylinder_365.png

The number of these elements (outliers)

Fitting_Cylinder_366.png

Fitting_Cylinder_367.png

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)

Fitting_Cylinder_368.png

and

Fitting_Cylinder_369.png

Fitting_Cylinder_370.png

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

Fitting_Cylinder_371.png

Fitting_Cylinder_372.png

Fitting_Cylinder_373.png

Fitting_Cylinder_374.gif

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

Fitting_Cylinder_375.png

Fitting_Cylinder_376.png

Fitting_Cylinder_377.png

Fitting_Cylinder_378.png

Fitting_Cylinder_379.png

Fitting_Cylinder_380.png

Fitting_Cylinder_381.png

Fitting_Cylinder_382.png

Fitting_Cylinder_383.png

Fitting_Cylinder_384.png

Fitting_Cylinder_385.png

Fitting_Cylinder_386.png

Fitting_Cylinder_387.png

Fitting_Cylinder_388.png

Fitting_Cylinder_389.png

Fitting_Cylinder_390.png

Fitting_Cylinder_391.png

Fitting_Cylinder_392.png

Fitting_Cylinder_393.gif

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

Now the parameters of the Gaussian mixture

Fitting_Cylinder_394.png

The result of the parameter values

Fitting_Cylinder_395.png

Fitting_Cylinder_396.png

These can be displayed in a table form

Fitting_Cylinder_397.png

Fitting_Cylinder_398.png

Table 6.
Parameters of the Gaussian mixture after zero iteration

The density functions of the two component.

Fitting_Cylinder_399.png

Fitting_Cylinder_400.gif

Fig 35. The PDF of the two components

Fitting_Cylinder_401.png

Fitting_Cylinder_402.gif

Fig 36. The joint PDF of the mixture

Fitting_Cylinder_403.png

Fitting_Cylinder_404.gif

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

Fitting_Cylinder_405.png

Fitting_Cylinder_406.png

Fitting_Cylinder_407.gif

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

Fitting_Cylinder_408.png

Fitting_Cylinder_409.gif

Fig 39. The convergence of the standard deviations

Fitting_Cylinder_410.png

Fitting_Cylinder_411.gif

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

Fitting_Cylinder_412.png

Fitting_Cylinder_413.png

Membership values of the second cluster

Fitting_Cylinder_414.png

Fitting_Cylinder_415.png

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

Fitting_Cylinder_416.png

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)

Fitting_Cylinder_417.png

The number of these elements (outliers)

Fitting_Cylinder_418.png

Fitting_Cylinder_419.png

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)

Fitting_Cylinder_420.png

and

Fitting_Cylinder_421.png

Fitting_Cylinder_422.png

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

Fitting_Cylinder_423.png

Fitting_Cylinder_424.png

Fitting_Cylinder_425.png

Fitting_Cylinder_426.gif

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

The inlier points,

Fitting_Cylinder_427.png

Fitting_Cylinder_428.gif

Fig. 35 The inliers

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

Fitting_Cylinder_429.png

Fitting_Cylinder_430.png

Fitting_Cylinder_431.png

Fitting_Cylinder_432.png

Fitting_Cylinder_433.gif

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

Fitting_Cylinder_434.png

Fitting_Cylinder_435.gif

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

Table 7

The computed cylinder parameters

Iteration a b c d r
0 -0.6619 -11.2767 -3.5779 146.344 0.1750
1 -0.6611 -11.2791 -3.6086 146.398 0.1753

Table 8

Parameters of the components of the Gaussian mixture

Iteration Fitting_Cylinder_436.png Fitting_Cylinder_437.png Fitting_Cylinder_438.png Fitting_Cylinder_439.png Fitting_Cylinder_440.png Fitting_Cylinder_441.png Fitting_Cylinder_442.png Fitting_Cylinder_443.png
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

Created with the Wolfram Language