Plane Fitting to Point Cloud via Gröbner Basis

Béla Paláncz, Árpád Somogyi, Tamás Lovas and Bence Molnár

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

Abstract

This notebook presents a new approach for surface reconstruction from point clouds. A fast algorithm with zero complexity to fit a plane to point clouds in Total Least Squares (TLS) sense has been developed using Gaussian-type error distribution. The application of maximum likelihood method results  a multivariate polynomial system which can be solved via symbolic computation using reduced Gröbner basis in analytic form.
In case of  plane fitting when datasets are affected heavily by outliers, this algorithm can be embedded into the RANdom SAmple Consensus (RANSAC) method. In order to increase the effectiveness of this integrated algorithm   Self-Organizing Map (SOM) can be employed as data preprocessing method to eliminate the main portion of the noisy data points by creating representative code-book vectors.  In this way the total running time can be decreased  considerably.
The study compares the performance of  the suggested algorithm to that of other methods and demonstrates its efficiency on both synthetic data adopted from the literature as well as on real, laser scanned measurements.

Keywords

Surface reconstruction, plane fitting, point cloud, LiDAR, TLS, Gröbner basis, outliers, RANSAC, SOM.

1 - Introduction

In many fields such as robotics (Poppinga  at al  2008), computer vision (Mitra and Nguyen 2003), digital photogrammetry (Yang and Förtsner 2010), surface reconstruction (Nurunnabi, Belton and West 2012), computational geometry (Lukács, Martin and Marshall 1998) as well as  in the increasing applications of laser scanning  (Stathas et al 2003) , it is a fundamental task for extracting plane features from 3D point cloud. Since the physical limitations of the sensors, the occlusions, multiple reflectance and noise can produce off-surface points, robust fitting techniques are required.  Perhaps the three most popular approaches for plane fitting in sense of Total Least Squares (TLS) are the Singular Values Decomposition  (SVD), the Principal Component Analysis (PCA) and the RANdom Sample Consensus (RANSAC) which can eliminate the so-called outliers using a noise threshold. For example Huang and Tseng (2008) employed Danish robust estimation with TLS as well as PCA. However, there are other type of methods, too e.g. that are based on minimum covariance determinant method (Russeeuw and Van Driessen, 1999),  expectation maximization method  (Lakaemper and Latecki 2006),  Bayesian techniques (Diebel et al, 2006),   region growing algorithm (Chen and Stamos, 2007) and so on. In this computational study we introduce a new efficent method based on computer algebra.

The notebook is arranged as follows:  Section 2,  we consider a toy problem and solve it using SVD and  PCA techniques. In Section 3 a new algorithm is developed based on the maximum likelihood estimation employing TLS error definition and Gaussian error distribution. The related  optimization problem is converted into a multivariate polynomial system which can be solved via reduced Gröbner basis in symbolic form. This results a non-iterative formula which has zero complexity and practically is as fast the PCA method even, when the solution is carried out -  from practical point of view-  with numerical Gröbner basis. In Section 4 these 3 methods - SVD, PCA and the Gröbner basis solution - are compared on the basis of a synthetic dataset adopted from the literature as a benchmark. In Section 5 the problem of outliers is illustrated. In Section 6 the Gröbner basis solution is integrated into RANSAC and applied to an artificial dataset corrupted with outliers. In Section 7 SOM is employed as a data preprocessing technique which reduces the amount of points to some tens of representative points, so called code-book vectors. In this way the running time of the SOM and that of the following RANSAC can be half of the running time of RANSAC with the original point of cloud. 
In Section 8 real laser scanner measurement of  grassy natural slope has been investigated. Computations with and without RANSAC as well as with SOM data preprocessing have been carried out and their results have been compared. These computations also verify the efficiency of the suggested algorithm based on Gröbner basis as well as that of the SOM preprocessing technique. In Appendix we demonstrate the known fact that PCA is not proper for RANSAC, (Mitra and Nguyen 2003) therefore our suggested algorithm is reasonable.
The computations have been carried out on HP Z 420 workstation with configuration of  Xeon E5-1650, 3.2 GHz processor and 16 GB RAM.

2 - Related works and principles

2 - 1 Total Least Squares Approach

In case of ordinary Least Squares (OLS) the residuals are the vertical distances between the points and the plane to be fitted. To overcome the bias to one direction TLS can be employed that minimizes the squared sum of the orthogonal distances between the points and the plane. This is the geometrical interpretation of the TLS. However, there is an equivalent algebraic interpretation, too. In this case one should adjust all of the measured ((Subscript[x, i]), (Subscript[y, i]) (Subscript[z, i])) coordinates in order to satisfy the plane equation, namely

(Subscript[z, i])+(Subscript[Δz, i])= α ((Subscript[x, i]) +(Subscript[Δx, i]))+β ((Subscript[y, i])+(Subscript[Δy, i]))+γ=0,   i= 1,...,ℕ

and the sum of these squared adjustment must be minimized.

Therefore the objective function is

(∑_(i =1)^ℕ)((Subscript[Δx, i])^2)+((Subscript[Δy, i])^2)+((Subscript[Δz, i])^2)

and this should be minimized under the constrains of Eq. (1). Since the constrains are linear, eliminating (Subscript[Δz, i])

(∑_(i=1)^ℕ)((Subscript[Δx, i])^2)+((Subscript[Δy, i])^2)+(((α ((Subscript[x, i]) +(Subscript[Δx, i]))+β ((Subscript[y, i])+(Subscript[Δy, i]))+γ-(Subscript[z, i])))^2)

we get an unconstrained problem with 2 ℕ+3 variables instead of  3 ℕ+3.

Let us consider a toy problem of four  3D points to illustrate the method, see  Nievergelt (2000). The coordinates of the measured points are,

XYZ={{11,45,38},{47,54,38},{17,12,14},{21,29,58}};

ℕ=Length[XYZ];

{x,y,z}=Transpose[XYZ];

The objective function is,

obj=Apply[Plus,Table[((Subscript[Δ𝓍, i])^2)+((Subscript[Δ𝓎, i])^2)+(((α (x[[i]]+(Subscript[Δ𝓍, i]))+β (y[[i]]+(Subscript[Δ𝓎, i]))+γ-z[[i]]))^2),{i,1,ℕ}]];

Now the unknown variables are not only the parameters (α, β, γ) but the adjustments ((Subscript[Δx, i]), (Subscript[Δy, i])), too.

vars=Join[Table[{(Subscript[Δ𝓍, i]),(Subscript[Δ𝓎, i])},{i,1,ℕ}],{α,β,γ}]//Flatten;

Let us employ built-in optimization method,

AbsoluteTiming[sol=NMinimize[obj,vars];]

Then the solution is

sol

Geometrically, the linear problem of total least squares amounts to fitting a hyperplane H minimizing the average squared Euclidean distance (measured perpendicularly to the fitted hyperplane) to data points (Subscript[c, 1]),...,  (Subscript[c, m])∈(ℂ^(n+1)). The problem then reduces to finding a point (Subscript[c, 0])∈H and a non-zero normal vector x orthogonal to H that minimizes the sum of the squared distances,

D (x,(Subscript[c, 0]))= (∑_(i=1)^m)(((|⟨(Subscript[c, i])-(Subscript[c, 0]),x⟩|)^2)/(⟨x,x⟩))

Introducing a matrix (Subscript[C, (Subscript[c, 0])]) with i-th rows (Subscript[c, i]) - (Subscript[c, 0]),

(Subscript[C, (Subscript[c, 0])])=(([ (Subscript[c, 1])-(Subscript[c, 0]),..., (Subscript[c, m])-(Subscript[c, 0])])^T)

we obtain that

D(x,(Subscript[c, 0]))=(((Subscript[‖(Subscript[C, (Subscript[c, 0])])x‖, 2])^2)/((Subscript[‖x‖, 2])^2))

Now let us consider two lemmas. 
Lemma 1 : The optimal hyperplane must pass through the centroid of the data. Namely

(Subscript[c, 0])=(1/m)(∑_(i=1)^m)(Subscript[c, i])

Lemma 2 : For every point (Subscript[c, 0]), let v be a singular vector corresponding to the smallest singular value σ of  matrix (Subscript[C, (Subscript[c, 0])]), then

D (x,(Subscript[c, 0]))≥  D (v,(Subscript[c, 0]))

Consequently the algorithm to find the solution of the linear total least squares problem is as follows (see Nievergelt  2000)

Consider the four data points in space,

(Subscript[c, 1])=(({{11}, {45}, {38}})),(Subscript[c, 2])=(({{47}, {54}, {38}})),(Subscript[c, 3])=(({{17}, {12}, {14}})),(Subscript[c, 4])=(({{21}, {29}, {58}}))

The centroid of the data,

(Subscript[c, 0])=(1/4)(∑_(i=1)^4)(Subscript[c, i])= (({{24}, {35}, {37}}))

c0=Mean[XYZ]

Then the matrix (Subscript[C, (Subscript[c, 0])]) is,

(Subscript[C, (Subscript[c, 0])])=(({{-13, 10, 1}, {23, 19, 1}, {-7, -23, -23}, {-3, -6, 21}}))

Cc0=Map[#-c0&,XYZ];Cc0//MatrixForm

The smallest singular value of this matrix is σ = (Subscript[σ, 3]) = 18 and the corresponding singular vector is

v = (Subscript[v, 3])=(({{(2/3)}, {-(2/3)}, {(1/3)}}))

AbsoluteTiming[{u,w,v}=SingularValueDecomposition[Cc0];]

w//MatrixForm

ν=Last[Transpose[v]]

({𝓍,𝓎,𝓏}-c0).ν==0//Simplify

Solve[%,𝓏]//Flatten

Thus the hyperplane H passes through (Subscript[c, 0]) and perpendicular to the vector v . Consequently the hyperplane is

(2/3)(x-24)-(2/3)(y-35)+(1/3)(z-37) = 0

or in explicit form

z = 15 - 2 x + 2 y

Fig. 1  The plane fitted in total least squares sense to the data points projected perpendicular to the plane

2 - 2 Statistical approach

Principal component analysis is a popular statistical technique that describes the covariance structure of data by means of a small number of components. These components are linear combinations of the original variables that  rank the variability  in the data through the variances, and produce directions using the eigenvectors of the covariance matrix. The eigenvector corresponding to the smallest eigenvalue is exactly the normal of the best-fitted plane. 

The covariance of matrix (Subscript[C, (Subscript[c, 0])]) is

AbsoluteTiming[A=Covariance[Cc0];MatrixForm[A]]

The eigenvalues of A are

Eigenvalues[A]

The smallest eigenvalue is (Subscript[λ, 3])= 108 and the corresponding eigenvector is

{nx,ny,nz}=Eigenvectors[A][[3]]

The Hesse model representing the plane via its normal vector is

(Subscript[n, x]) x +(Subscript[n, y]) y+(Subscript[n, z])z-d=0

where the normal vector of the plane is ((Subscript[n, x]),(Subscript[n, x]),(Subscript[n, x])) . In order to get the d parameter, let us consider the residual of this model

R ((Subscript[n, x]) ,(Subscript[n, y]),(Subscript[n, z]),d) =(∑_(i=1)^ℕ)((((Subscript[n, x]) (Subscript[x, i]) +(Subscript[n, y]) (Subscript[y, i])+(Subscript[n, z])(Subscript[z, i])-d))^2)

therefore considering that

((∂R)/(∂ d))=0

we get

(Subscript[n, x]) (∑_(i=1)^ℕ)(Subscript[x, i])+(Subscript[n, y]) (∑_(i=1)^ℕ)(Subscript[y, i])+(Subscript[n, z]) (∑_(i=1)^ℕ)(Subscript[z, i])- d ℕ = 0

Then

sold=Solve [{nx,ny,nz}.Total[XYZ]-𝒹 ℕ==0,𝒹]//Flatten

Consequently the surface is

z =- ((Subscript[n, x])/(Subscript[n, z])) x -((Subscript[n, y])/(Subscript[n, z])) y +(d/(Subscript[n, z])) =-2 x+ 2 y + 15

{-(nx/nz),-(ny/nz),(𝒹/nz)}/.sold

Unfortunately, both the classical variance (which is being maximized) and the classical covariance matrix (which is being decomposed) are very sensitive to anomalous observations. Consequently, the first components are often attracted toward outlying points, and may not capture the variation of the regular observations. Therefore, data reduction based on PCA becomes unreliable if outliers  are present in the data,  see Mitra and Ngugen (2003) or Hubert at al (2005). Although there are efforts to improve the robustness of PCA with adding complementary outlier identification method, e.g.. Nurunnabi et al, (2012), the TLS (SVD) has advantages to PCA since it can describe the similarity as well as the correlation amongst the features of data both together, while PCA can consider only one or the other of them selectively.

That is why, preferring to TLS approach,  a new fast algorithm with zero complexity will be developed for TLS in the next Section.

3 - TLS solution via Gröbner Basis

3 -1 Error model

Clear["Global`*"]

Let a point be (Subscript[P, i]) with coordinates  {(Subscript[𝓍, i]),(Subscript[𝓎, i]),(Subscript[𝓏, i])}  and  (Subscript[𝒹, i]) is its  distance  from a general point {x,y,z}
of the plane  z=α x+β y +γ, then

(Subscript[𝒹, i])=((((Subscript[𝓏, i])-α x-β y -γ))^2)+((((Subscript[𝓍, i])-x))^2)+((((Subscript[𝓎, i])-y))^2)

The distance of (Subscript[P, i]) from its perpendicular  projection to the plane is the shortest distance, which represents its regression error (Subscript[(Subscript[ℯ, ℳ]), i]). From the necessary conditions we get

deq1=D[(Subscript[𝒹, i]),x]//Expand

deq2=D[(Subscript[𝒹, i]),y]//Expand

xy=Solve[{deq1==0,deq2==0},{x,y}]//Flatten

Then the square of the error is

(Subscript[Δ, i])=(Subscript[𝒹, i])/.xy//Simplify

Consequently our error model is

ϵ=(((Subscript[𝓏, i])-(Subscript[𝓍, i]) α-(Subscript[𝓎, i]) β-γ)/(Sqrt[1+(α^2)+(β^2)]));

3 - 2 Application of Maximum Likelihood Estimation

In order to avoid direct minimization and to get explicit formula for the estimated parameters symbolic computation can be employed. 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.

SuperLog[Q_]:=Module[{erk,iii,nnn},
Product[iii,{iii,nnn}];(*pre-load Product*)
Which[
Q===On,
Unprotect[Log];Clear[Log];
Log[Product[x_,{k_.a_,b_}]]:=
Log[Product[Times[erk,x],{k,a,b}]]/.erk->1;
Log[Product[HoldPattern[Times[x__]],{k_,a_,b_}]]:=
Simplify[Map[Sum[#,{k,a,b}]&,Plus@@Map[Expand[PowerExpand[Log[#]]]&,List[x]]]//.
Sum[u_,w_,{kk_,aa_,bb_}]:>u Sum[w,{kk,aa,bb}]/;FreeQ[u,kk]==True];
Protect[Log];Print["SuperLog is now On"],
Q===Off,
Unprotect[Log];Clear[Log];Protect[Log];Print["SuperLog is now Off"],
True,
Print["Please use SuperLog[On] or SuperLog[Off]!"]
]]

Let us activate this function

SuperLog[On]

Then the Maximum Likelihood (ML) estimator for Gaussian-type noise is

ℒ=(∏_(i=1)^ℕ)PDF[NormalDistribution[0,σ],(Subscript[(Subscript[ℯ, ℳ]), i])]

Considering its logarithm which should be maximized

Logℒ=Log[ℒ/.(Subscript[(Subscript[ℯ, ℳ]), i])->ϵ]//Simplify

From the necessary conditions one can obtain

eq1=D[Logℒ,α]

eq2=D[Logℒ,β]

eq3=D[Logℒ,γ]

We get a multivariate algebraic system for the unknown parameters α, β and γ. In compact form

Clear[a,b,c,d,e,f,g,h,i]

eq1=ℕ α (γ^2)+2 (α^2)γ a-γ a (1+(α^2)+(β^2))+(α^3)b-α b(1+(α^2)+(β^2))+2 α β γ c+α (β^2)d+2 (α^2)β e-β e(1+(α^2)+(β^2))-2 α γ f-2 α β g+α h -2 (α^2)i+i(1+(α^2)+(β^2))//Expand

eq2=ℕ β (γ^2)+2 α β γ a+(α^2)β b+2 (β^2)γ c-γ c(1+(α^2)+(β^2))+2 α (β^2)e-α e(1+(α^2)+(β^2))+(β^3)d-β d(1+(α^2)+(β^2))-2 β γ f-2 α β i+β h-2 (β^2)g+g (1+(α^2)+(β^2))//Expand

eq3=f-ℕ γ-α a-β c

where the constants are

a=(∑_(i=1)^ℕ)(Subscript[𝓍, i]),      b=(∑_(i=1)^ℕ)((Subscript[𝓍, i])^2),     c=(∑_(i=1)^ℕ)(Subscript[𝓎, i]),     d=(∑_(i=1)^ℕ)((Subscript[𝓎, i])^2),     e=(∑_(i=1)^ℕ)(Subscript[𝓍, i])(Subscript[𝓎, i]),     f=(∑_(i=1)^ℕ)(Subscript[𝓏, i]),     g=(∑_(i=1)^ℕ)(Subscript[𝓎, i])(Subscript[𝓏, i]),     h=(∑_(i=1)^ℕ)((Subscript[𝓏, i])^2),     i=(∑_(i=1)^ℕ)(Subscript[𝓍, i])(Subscript[𝓏, i])

Similarly the compact form of objective function, which will be required later

L=-((ℕ (γ^2))/(2 (1+(α^2)+(β^2)) (σ^2)))-(1/2) ℕ Log[2]-(1/2) ℕ Log[π]-ℕ Log[σ]- a ((α γ)/((1+(α^2)+(β^2)) (σ^2)))-(((α^2) b)/(2 (1+(α^2)+(β^2)) (σ^2)))-((β γ c)/((1+(α^2)+(β^2)) (σ^2)))-((α β e)/((1+(α^2)+(β^2)) (σ^2)))-(((β^2) d)/(2 (1+(α^2)+(β^2)) (σ^2)))+((γ f)/((1+(α^2)+(β^2)) (σ^2)))+((α i)/((1+(α^2)+(β^2)) (σ^2)))+((β g)/((1+(α^2)+(β^2)) (σ^2)))-(h/(2 (1+(α^2)+(β^2)) (σ^2)))

This technique is similar to finding regression parameters of a straight line, see Paláncz (2014).

3 - 3 Solution via Symbolic Gröbner Basis

The polynomial system can be reduced to monomials of higher order. First let us eliminate α via Gröbner basis, we get

gbαβ=GroebnerBasis[{eq1,eq2,eq3},{α,β,γ},{γ},MonomialOrder->EliminationOrder]

The basis consists of 7 polynominals of variables α and β. The length of these polynomials are

Map[Length[#]&,gbαβ]

Let us consider the first and the third one and eliminate β via Sylvester resultant,

alfa=Resultant[gbαβ[[1]],gbαβ[[3]],β]

so we get a univariate  polynomial for  α.  Similarly eliminate α  from the first and second basis,

beta=Resultant[gbαβ[[1]],gbαβ[[2]],α]

which leads to a univariate polynomial of β.

The coefficients of the uni-variable polynomial for α are

ca0=Coefficient[alfa,α,0]

ca1=Coefficient[alfa,α,1]

ca2=Coefficient[alfa,α,2]

ca3=Coefficient[alfa,α,3];Short[ca3,10]

ca4=Coefficient[alfa,α,4];Short[ca4,10]

ca5=Coefficient[alfa,α,5];Short[ca5,10]

ca6=Coefficient[alfa,α,6];Short[ca6,10]

ca7=Coefficient[alfa,α,7];Short[ca7,10]

ca8=Coefficient[alfa,α,8];Short[ca8,10]

ca9=Coefficient[alfa,α,9];Short[ca9,10]

ca10=Coefficient[alfa,α,10]

Therefore our polynomial for α can be written as

A=(ca2 +ca3 α+ca4 (α^2)+ca5 (α^3)+ca6 (α^4)+ca7 (α^5)+ca8 (α^6)+ca9 (α^7)) (α^2);

Let us check these coefficients

(alfa-A)//Simplify

Now, the coefficients of the  polynomial of β

cb0=Coefficient[beta,β,0]

cb1=Coefficient[beta,β,1]

cb2=Coefficient[beta,β,2];Short[cb2,10]

cb3=Coefficient[beta,β,3];Short[cb3,10]

cb4=Coefficient[beta,β,4];Short[cb4,10]

cb5=Coefficient[beta,β,5];Short[cb5,10]

cb6=Coefficient[beta,β,6];Short[cb6,10]

cb7=Coefficient[beta,β,7];Short[cb7,10]

cb8=Coefficient[beta,β,8];Short[cb8,10]

cb9=Coefficient[beta,β,9];Short[cb9,10]

cb10=Coefficient[beta,β,10]

then

B=(cb2 +cb3 β+cb4 (β^2)+cb5 (β^3)+cb6 (β^4)+cb7 (β^5)+cb8 (β^6)+cb9 (β^7)) (β^2);

(beta-B)//Simplify

Employing the numerical values of the toy problem

XYZ={{11,45,38},{47,54,38},{17,12,14},{21,29,58}}

X=XYZ[[All,1]]; Y=XYZ[[All,2]];Z=XYZ[[All,3]];

{a,b,c,d,e,f,g,h,i}=Total/@{X,X^2,Y,Y^2,X Y,Z,Y Z,Z^2,X Z};

ℕ=Length[XYZ]

The monomial for α

AbsoluteTiming[pα=ca2 +ca3 α+ca4 (α^2)+ca5 (α^3)+ca6 (α^4)+ca7 (α^5)+ca8 (α^6)+ca9 (α^7);]

pα

The roots of this monomial are

AbsoluteTiming[solα=Reduce[pα==0,α,Reals];]

solα

ToRules[solα[[1]]]

Similarly the monomial for β and γ, and their solutions are

AbsoluteTiming[pβ=cb2 +cb3 β+cb4 (β^2)+cb5 (β^3)+cb6 (β^4)+cb7 (β^5)+cb8 (β^6)+cb9 (β^7);]

pβ

solβ=Reduce[pβ==0,β,Reals]

ToRules[solβ[[1]]]

pγ=Solve[eq3==0,γ]//Flatten

%/.ToRules[solα[[1]]]/.ToRules[solβ[[1]]]

AbsoluteTiming[gbαβ;]

3 - 4 Solution via Numerical Gröbner Basis

The key problem with the explicit Gröbner basis solution is the proper selection of the corresponding roots. Therefore direct numerical solution of the polynomial system via numerical Gröbner basis employing the built-in function NSolve has been chosen.

eq1

eq2

eq3

AbsoluteTiming[solαβγ=NSolve[{eq1,eq2,eq3},{α,β,γ}];]

The running time is acceptably short.

solαβγ

Let us select the real solutions

selsol=Select[solαβγ,Im[#[[1]][[2]]]==0&]

How can we  select the proper real solution? Let us consider the values of the formal objective function based on the different real solutions

Lαβγ=Map[L/.#&,selsol]/.σ->1//N

The real solution will be selected accordingly which maximizes the objective function. The value of σ will not influence the relative comparison of the real solutions, therefore we consider σ = 1. Then the optimal solution, which gives the highest value for the maximum likelihood function is

solopt=selsol[[Position[Lαβγ,Max[Lαβγ]]//Flatten//First]]

4 - Laser Scanning Example 1

4 - 1 Artificial Data

LiDAR (Light Detection and Ranging) point clouds contain abundant spacial (3D) spatial information. Dense distribution of scanned points on object surfaces prominently implies surface features. Particularly, plane features commonly appear in a typical LiDAR dataset of artificial structures. Plane fitting is the key process for extracting plane features from LiDAR data.
Synthetic dataset are generated from a multivariate Gauss distribution. Regular 3D points have means of  (2, 8, 6) and variances (5, 5, 0.01), see Nurunnabi, Belton and West (2012). The ideal plane is z = 6.

In this case let us consider here ℕ =20000 points.

Clear[XYZ]

ℕ=20000;

XYZ=Table[{Null,Null,Null},{i,1,ℕ}];

Do[XYZ[[i]]=RandomVariate[MultinormalDistribution[{2.,8,6},{{5.,0,0},{0,5.,0},{0,0,0.01}}]],{i,1,ℕ}];

p1=Plot3D[6,{x,-8,12},{y,-2,18},BoxRatios->{1,1,0.5}, ClippingStyle->None,Mesh->None,PlotRange->{5.5,6.5}];

p2=ListPointPlot3D[XYZ,PlotStyle->Blue];

Figure 2 shows the data points with the ideal plane.

Show[{p1,p2}]

Fig. 2 The ideal plane and the  synthetic measured points

4 - 2 Solution via SVD

c0=Mean[XYZ]

Cc0=Map[#-c0&,XYZ];

AbsoluteTiming[{u,w,v}=SingularValueDecomposition[Cc0];]

ν=Last[Transpose[v]]

({𝓍,𝓎,𝓏}-c0).ν==0//Simplify

NSolve[%,𝓏]//Expand

4 - 3 Solution via PCA

AbsoluteTiming[A=Covariance[Cc0];MatrixForm[A]]

{nx,ny,nz}=Eigenvectors[A][[3]]

sold=Solve [{nx,ny,nz}.Total[XYZ]-𝒹 ℕ==0,𝒹]//Flatten

Consequently the coefficients of surface model is

{-(nx/nz),-(ny/nz),(𝒹/nz)}/.sold

4 - 4 Solution via Numerical Gröbner Basis

X=XYZ[[All,1]]; Y=XYZ[[All,2]];Z=XYZ[[All,3]];

{a,b,c,d,e,f,g,h,i}=Total/@{X,X^2,Y,Y^2,X Y,Z,Y Z,Z^2,X Z}

ℕ=Length[XYZ];

eq1=ℕ α (γ^2)+2 (α^2)γ a-γ a (1+(α^2)+(β^2))+(α^3)b-α b(1+(α^2)+(β^2))+2 α β γ c+α (β^2)d+2 (α^2)β e-β e(1+(α^2)+(β^2))-2 α γ f-2 α β g+α h -2 (α^2)i+i(1+(α^2)+(β^2));
eq2=ℕ β (γ^2)+2 α β γ a+(α^2)β b+2 (β^2)γ c-γ c(1+(α^2)+(β^2))+2 α (β^2)e-α e(1+(α^2)+(β^2))+(β^3)d-β d(1+(α^2)+(β^2))-2 β γ f-2 α β i+β h-2 (β^2)g+g (1+(α^2)+(β^2));
eq3=f-ℕ γ-α a-β c;

AbsoluteTiming[solαβγ=NSolve[{eq1,eq2,eq3},{α,β,γ}];]

selsol=Select[solαβγ,Im[#[[1]][[2]]]==0&];
Lαβγ=Map[L/.#&,selsol]/.σ->1//N;
solopt=selsol[[Position[Lαβγ,Max[Lαβγ]]//Flatten//First]]

4 - 5 Comparison of the running times

Table 1 and Figure 3 shows the running times of the different methods in case of different dataset sizes.

Table 1.
Running times in sec

({{Number of data points, SVD, Numerical Gröbner Basis, PCA}, {(10^3), 0.218, 0.047, 0.000}, {5×(10^3), 0.250, 0.094, 0.000}, {(10^4), 1.076, 0.062, 0.016}, {1.5×(10^4), 2.309, 0.094, 0.016}, {2×(10^4), 4.072, 0.093, 0.000}, {2.5 ×(10^4), 6.599, 0.078, 0.016}})

dataSVD={{1000,0.218},{5000,0.25},{10000,1.076},{15000,2.309},{20000,4.072},{25000,6.599}};

dataGB={{1000,0.047},{5000,0.094},{10000,0.062},{15000,0.094},{20000,0.093},{25000,0.078}};

dataPCA={{1000,0.},{5000,0.0},{10000,0.016},{15000,0.016},{20000,0.0},{25000,0.016}};

ListPlot[{dataGB,dataSVD,dataPCA},Joined->True,Frame->True,FrameLabel->{"Number of data","Running Time"},Axes->False]

Fig. 3 Running time vs. Number of points

Both PCA and numerical Gröbner basis solutions have zero complexity with very low running time. PCA seems to be faster, but this fact has no practical significance.

4 - 6 Complexity of the SVD algorithms

According to Figure 3 one may approach the relation between the running time (T) and the number of data (N) in the following form (assuming polynomial complexity),

T∼ (N^η)

Therefore the complexity can be measured with η.

Fit[dataSVD//Log,{1,LogN},LogN]

Therefore  complexity of the SVD solution based on our numerical experiment  is η ∼1.1 while that of the Gröbner Basis as well as the PCA solution is  η ∼ 0!

Remark
Considering a matrix m ×n  the theoretical computation time of the SVD algorithm (for w and v only) is ∼ 4 (mn^2)+8 (n^3). Since in our case n = 3, consequently the theoretical time complexity is η = 1, see campar - wwwnavab.cs.tum.edu - 3D computer Vision Script Draft v.000 (downloaded at  28.11. 2013).

5 - Outliers

5 - 1 Problem of Outliers

A definition of the outliers provided by Zuliani (2012) is as follows: a datum is considered to be an outlier if it will not fit the “true” model instantiated by the “true” set of  parameters within some error threshold that defines the maximum deviation attributable to the effects of noise. The estimation of the parameters of a model in presence of outliers constitutes a very tricky and challenging problem although it belongs to the daily business activities of modern geodesists. There are many different techniques to solve this problem. Here we suggest a new algorithm where our numerical Gröbner basis solution is integrated into the RANdom Sample Consensus (RANSAC) algorithm. However, first let us see, how the estimators SVD, PCA and numerical Gröbner basis can handle outliers without any data preparation.

5 - 2 Artificial Data

For outliers also artificial datasets are generated from a multivariate Gauss distribution. The outliers have means of  (15, 15, 10) and variances (10, 2, 1), see Nurunnabi, Belton and West (2012). The ideal plane is z = 6, see Figure 4.

There are additional 3000 points, which represent outliers.

ℕ=20000;ℕOutLiers=3000;

Now, let us generate 3000 outliers randomly

XYZOutLiers=Table[{Null,Null,Null},{i,1,ℕOutLiers}];

Do[XYZOutLiers[[i]]=RandomVariate[MultinormalDistribution[{15.,15,10},{{10.,0,0},{0,2.,0},{0,0,1}}]],{i,1,ℕOutLiers}];

p1=Plot3D[6,{x,-8,35},{y,-2,19},BoxRatios->{1,1,0.5}, ClippingStyle->None,Mesh->None,PlotRange->{4.5,12}];

p3=ListPointPlot3D[XYZOutLiers,PlotStyle->Red];

Show[{p1,p2,p3}]

Fig. 4  The plane and the artificial measured points with outliers (red)

5 - 3 Solution via SVD

Let us compute the parameters with SVD

XYZC=Join[XYZ,XYZOutLiers];

Length[XYZC]

c0=Mean[XYZC]

Cc0=Map[#-c0&,XYZC];

AbsoluteTiming[{u,w,v}=SingularValueDecomposition[Cc0];]

ν=Last[Transpose[v]]

({𝓍,𝓎,𝓏}-c0).ν==0//Simplify

NSolve[%,𝓏]//Expand

The outliers changed the true parameters considerably, see Figure 5. Let us compute the parameters with PCA, too

p4=Plot3D[𝓏/.%[[1]],{𝓍,-8,35},{𝓎,-3,19},BoxRatios->{1,1,0.5}, ClippingStyle->None,Mesh->None,PlotRange->{3,15},ColorFunction->"RustTones"];

Show[{p4,p1}]

Fig. 5 The original plane (white) and the estimated plane (rust) in case of outliers

5 - 4 Solution via PCA

AbsoluteTiming[A=Covariance[Cc0];MatrixForm[A]]

{nx,ny,nz}=Eigenvectors[A][[3]]

sold=Solve [{nx,ny,nz}.Total[XYZC]-𝒹 ℕ==0,𝒹]//Flatten

Consequently the coefficients of surface model are

{-(nx/nz),-(ny/nz),(𝒹/nz)}/.sold

The result is quite the same, only the constant member (γ) has different value compared to the SVD solution.

5 - 5 Solution via Numerical Gröbner Basis

X=XYZC[[All,1]]; Y=XYZC[[All,2]];Z=XYZC[[All,3]];

{a,b,c,d,e,f,g,h,i}=Total/@{X,X^2,Y,Y^2,X Y,Z,Y Z,Z^2,X Z};

ℕℂ=Length[XYZC]

eq1=ℕℂ α (γ^2)+2 (α^2)γ a-γ a (1+(α^2)+(β^2))+(α^3)b-α b(1+(α^2)+(β^2))+2 α β γ c+α (β^2)d+2 (α^2)β e-β e(1+(α^2)+(β^2))-2 α γ f-2 α β g+α h -2 (α^2)i+i(1+(α^2)+(β^2));
eq2=ℕℂ β (γ^2)+2 α β γ a+(α^2)β b+2 (β^2)γ c-γ c(1+(α^2)+(β^2))+2 α (β^2)e-α e(1+(α^2)+(β^2))+(β^3)d-β d(1+(α^2)+(β^2))-2 β γ f-2 α β i+β h-2 (β^2)g+g (1+(α^2)+(β^2));
eq3=f-ℕℂ γ-α a-β c;

AbsoluteTiming[solαβγ=NSolve[{eq1,eq2,eq3},{α,β,γ}];]

selsol=Select[solαβγ,Im[#[[1]][[2]]]==0&];
Lαβγ=Map[L/.#&,selsol]/.σ->1//N;
solopt=selsol[[Position[Lαβγ,Max[Lαβγ]]//Flatten//First]]

The Gröbner basis solution gives the same result as SVD, while  PCA provides somewhat different value for the coefficient of γ (5.4.. ≠ 4.7..). This fact indicates that PCA is more sensitive to outliers than SVD or Gröbner. This feature of PCA causes that it can not be employed in RANSAC reliably, since its result is strongly effected by anomalous observations, see Appendix.

6 - Application of the Algorithm in RANSAC

6 - 1 RANSAC algorithm

Let us try to apply the RANSAC method, given in Yaniv (2010), which is proved to be successful for detecting outliers. The basic RANSAC algorithm is as follows:

1) Pick up a model type (ℳ)
2) Input data as
     -  dataQ  - data corrupted with outliers (cardinality (dataQ) = n)
     -  s     - number of data elements required per subset
     -  𝒩           -  number of subsets to draw from the data
     -  τ            -  threshold which defines if  data element, (Subscript[d, i])∈ dataQ, agrees with the model ℳ

Remarks

In general s can be the minimal number of the data which results a determined system for the unknown parameters of the model.
The number of subsets to draw from the data, 𝒩 is chosen high enough to ensure that at least one of the subsets of the random examples does not include an outliers (with the probability p, which is usually set to 0.99). Let u represent the probability that any selected data point is an inlier and v = 1 - u the probability of observing an outlier. Then the iterations 𝒩 can be computed as

𝒩 = ((log(1-p))/(log(1-(((1-v))^s))))

3) maximalConsensusSet ← ∅
4) Iterate 𝒩 times:
            a) ConsensusSet ← ∅
            b) Randomly draw a subset containing s elements and estimate the parameters of the model ℳ
            c) For each data element,  (Subscript[d, i])∈ dataQ:
                if  agree ((Subscript[d, i]),ℳ, τ), ConsensusSet ← (Subscript[d, i])
            d) if cardinality (maximalConsensusSet) < cardinality(ConsensusSet),
                maximalConsensusSet ←   ConsensusSet
                
   5) Estimate model parameters using maximalConsensusSet

6 - 2 The necessary number of iterations

In our case we have 3 parameters to determine, so we need minimum 3 equations, therefore

s=3;

In this case the parameters can be computed directly from a linear system. However, in general if  s > 3 we need numerical Gröbner solution. Let

p=0.99;

In addition let us assume that the probability of observing outlier is

ν=1/2;

Then the number of the necessary iterations is

𝒩=Round[Log[1-p]/Log[1-(((1-ν))^s)]]

6 - 3 Computation of the threshold

Let the threshold τ

τ=10;

An element  (Subscript[d, i])∈ dataQ will be accepted if the local error is less than the prespecified threshold

((|(Subscript[z, i])-(Subscript[α, s]) (Subscript[x, i])-(Subscript[β , s])(Subscript[y, i])-(Subscript[γ, s])|)/(Sqrt[1+(α^2)+(β^2)]))< τ

Remarks

(Subscript[d, i]) is represented by the triplet ((Subscript[x, i]), (Subscript[y, i]), (Subscript[z, i])) and the parameters ((Subscript[α, s]),(Subscript[β, s]),(Subscript[γ, s])) computed from a subset of s elements,

here ordinary least square ((Subscript[OLS, z])) estimation could be also applied as maximum support, see Fischler and Bolles (1981).

6 - 4 The computation model step by step

inliers=Table[{},{𝒩}] (* the consensus sets, j = 1,2, ...𝒩*)

Let us consider a subset  (j = 1,2, ...𝒩).

j=2 (*second iteration means second subset*)

Assigning noisy data to dataQ

dataQ=XYZC;

Let us select a triplet randomly

dataS=RandomSample[dataQ,s]

For further computation we need a function form of our algorithm.

6 - 5 The Mathematica function for the Numerical Gröbner basis Solution

We want to integrate numerical Gröbner basis solution into the RANSAC algorithm, therefore it is reasonable to write a function for this suggested algorithm,

PlaneTLS[XYZ_]:=Module[{X,Y,Z,a,b,c,d,e,f,g,h,i,ℕ,eq1,eq2,eq3,solαβγ,selsol,Lαβγ,L,solopt},
 X=XYZ[[All,1]]; Y=XYZ[[All,2]];Z=XYZ[[All,3]];
{a,b,c,d,e,f,g,h,i}=Total/@{X,X^2,Y,Y^2,X Y,Z,Y Z,Z^2,X Z};
    ℕ=Length[XYZ];
eq1=ℕ α (γ^2)+2 (α^2)γ a-γ a (1+(α^2)+(β^2))+(α^3)b-α b(1+(α^2)+(β^2))+2 α β γ c+α (β^2)d+2 (α^2)β e-β e(1+(α^2)+(β^2))-2 α γ f-2 α β g+α h -2 (α^2)i+i(1+(α^2)+(β^2));
eq2=ℕ β (γ^2)+2 α β γ a+(α^2)β b+2 (β^2)γ c-γ c(1+(α^2)+(β^2))+2 α (β^2)e-α e(1+(α^2)+(β^2))+(β^3)d-β d(1+(α^2)+(β^2))-2 β γ f-2 α β i+β h-2 (β^2)g+g (1+(α^2)+(β^2));
eq3=f-ℕ γ-α a-β c;
solαβγ=NSolve[{eq1,eq2,eq3},{α,β,γ}];
selsol=Select[solαβγ,Im[#[[1]][[2]]]==0&];
L=-(h/(2 (1+(α^2)+(β^2)) (σ^2)))+((i α)/((1+(α^2)+(β^2)) (σ^2)))-((b (α^2))/(2 (1+(α^2)+(β^2)) (σ^2)))+((g β)/((1+(α^2)+(β^2)) (σ^2)))-((e α β)/((1+(α^2)+(β^2)) (σ^2)))-((d (β^2))/(2 (1+(α^2)+(β^2)) (σ^2)))+((f γ)/((1+(α^2)+(β^2)) (σ^2)))-((a α γ)/((1+(α^2)+(β^2)) (σ^2)))-((c β γ)/((1+(α^2)+(β^2)) (σ^2)))-((ℕ (γ^2))/(2 (1+(α^2)+(β^2)) (σ^2)))-(1/2) ℕ Log[2]-(1/2) ℕ Log[π]-ℕ Log[σ];
Lαβγ=Map[L/.#&,selsol]/.σ->1//N;
solopt=selsol[[Position[Lαβγ,Max[Lαβγ]]//Flatten//First]]]

Now let us compute the parameters from the randomly selected subset

sols=PlaneTLS[dataS]

Then check for all elements of dataQ (i = 1,2, ..., n)  whether it can satisfy the threshold. For example let

i=1131;

dataQ[[i]]

If[((((Abs[dataQ[[i]][[3]]-α dataQ[[i]][[1]]-β dataQ[[i]][[2]]-γ])/(Sqrt[1+(α^2)+(β^2)])))/.sols)<τ,inliers[[j]]=Append[inliers[[j]],dataQ[[i]]]]

inliers

Employing the following double For cycle

For j = 1, 2, ... 𝒩
  For i = 1,2,...ℕ

we can get the maximal consensus whose size is between s and n, so mostly overdetermined. Then the computation is  straightforward.

6 - 6 Application of the RANSAC model

Now let us use for threshold  τ = 1 m

τ=1;

and let

s=3;

therefore

𝒩=Round[Log[1-p]/Log[1-(((1-ν))^s)]]

Preparation of the sets of inliers for every iteration step - at the beginning they are empty sets:

inliers=Table[{},{𝒩}] (* the consensus sets, j = 1,2, ...𝒩*);

Then let us carry out the double  For cycle to find the maximum consensus set

AbsoluteTiming[For[j=1,j≤𝒩,j++,
(*get three random points*)
dataS=RandomSample[dataQ,s];
(*determine model parameters*)
sols=PlaneTLS[dataS];
(* determine inliers*)
For[i=1,i≤ℕ,i++,
If[(((1/(Sqrt[1+(α^2)+(β^2)]))Abs[dataQ[[i]][[3]]-α dataQ[[i]][[1]]-β dataQ[[i]][[2]]-γ])/.sols)<τ,inliers[[j]]=Append[inliers[[j]],dataQ[[i]]]];
];];]

Select the consensus set which has the highest number of inliers

trimmeddata=inliers[[Ordering[inliers,-1]]][[1]];

Length[trimmeddata]

p5=ListPointPlot3D[trimmeddata,PlotStyle->Green];

Show[{p1,p5}]

Fig. 6 The original plane and the inliers computed via RANSAC  (green)

Then we can compute the parameters

PlaneTLS[trimmeddata]

7 - Application of Self-Organizing Map (SOM)

The Self-Organizing Map (SOM) is a neural network algorithm, which uses a competitive learning technique to train itself in an unsupervised manner. SOMs are different from other artificial neural networks in the sense that they use a neighborhood function to preserve the topological properties of the input space and they are used to create an ordered representation of multi-dimensional data which simplifies complexity and reveals meaningful relationships. 
In the last decade a lot of research has been carried out to develop surface reconstruction methods. A more recent approach to the problem of surface reconstruction is that of learning based methods, see e.g. DalleMole et al (2010). Learning algorithms are able to process very large and/or noisy data, such as point clouds obtained from 3D scanners and are used to construct surfaces. Following this approach some studies have been employed SOM and their variants for surface reconstruction. SOM is suitable for this problem because it can form topological maps representing the distribution of input data. In our case this mapping occurs from 3D space to 2D one , see Kohonen  (1998).                                                                                                                

Therefore we try to represent this cloud of data points by considerably less points using Kohonen-map with a lattice of  SOM of size 6 × 6 code-book vectors with symmetrical neighborhood

<<NeuralNetworks`

AbsoluteTiming[{som,fitrecord}=UnsupervisedNetFit[XYZC,36,35,SOM->{6,6},Neighbor->Symmetric];]

Fig. 7 The convergence of the SOM iteration

dataSOM=som[[1]];

p6=ListPointPlot3D[dataSOM,PlotStyle->Directive[Magenta,PointSize[0.015]]];

Show[{p1,p6}]

Fig. 8 The original plane and the code-book vectors of the SOM of 6 × 6

Now let us apply RANSAC  to these 36 code-book vectors with a low threshold

τ=0.05;

NSOM=Length[dataSOM];

inliers=Table[{},{𝒩}] (* the consensus sets, j = 1,2, ...𝒩*);

AbsoluteTiming[For[j=1,j≤𝒩,j++,
(*get three random points*)
dataS=RandomSample[dataSOM,s];
(*determine model parameters*)
sols=PlaneTLS[dataS];
(* determine inliers*)
For[i=1,i≤NSOM,i++,
If[(((1/(Sqrt[1+(α^2)+(β^2)]))Abs[dataSOM[[i]][[3]]-α dataSOM[[i]][[1]]-β dataSOM[[i]][[2]]-γ])/.sols)<τ,inliers[[j]]=Append[inliers[[j]],dataSOM[[i]]]];
];];]

Select the consensus set which has the highest number of inliers

trimmeddata=inliers[[Ordering[inliers,-1]]][[1]];

Length[trimmeddata]

p7=ListPointPlot3D[trimmeddata,PlotStyle->Directive[Magenta,PointSize[0.015]]];

Show[{p1,p7}]

Fig. 9 The original plane and the inliers of RANSAC  (magenta)

Then we can compute the parameters

PlaneTLS[trimmeddata]

The error of the method is quite acceptable. The total computation time of  SOM + RANSAC combination is less than that of the direct RANSAC with the original number of points.

8 - Laser Scanning Example 2

8 -1 Real measurement

Outdoor laser scanning measurements have been carried out in a hilly park of Budapest, see Fig.10. The test area is on a steep slope, covered with dense but low vegetation.

Fig. 10 The test area in Budapest

The experiment has been carried out with a Faro Focus 3D terrestrial laser scanner, see Fig.11. The test also aimed to investigate the tie point detection capabilities of the scanner’s processing software; different types of spheres were deployed all over the test area. In case of multiple scanning positions these spheres can be used for registering the point clouds.

Fig. 11 : Faro Focus 3D scanner

The measurement range of the scanner is 120 m, the ranging error is ± 2 mm, according to the manufacturer' s technical specification.

Fig. 12 The scanner at the top of the measured steep slope with the different sizes of white spheres as control points in the background

The scanning parameters were set to ½ resolution that equals to 3mm/10m point spacing. This measurement resulted 178.8 million points that was acquired in 5 and half minutes. 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 38 318 points in ASCII format, and only the x, y, z coordinates were kept (no intensity values). Let us load the measured data:

Clear[XYZ,ℕ]

XYZ=Import["D:\\slope.dat"];

ℕ=Length[XYZ]

p8=ListPointPlot3D[XYZ,PlotStyle->Blue]

Fig. 13 The measured LiDAR point cloud

8 - 2 Solution via SVD

c0=Mean[XYZ]

Cc0=Map[#-c0&,XYZ];

AbsoluteTiming[{u,w,v}=SingularValueDecomposition[Cc0];]

ν=Last[Transpose[v]]

({𝓍,𝓎,𝓏}-c0).ν==0//Simplify

NSolve[%,𝓏]//Expand

8 - 3 Solution via PCA

AbsoluteTiming[A=Covariance[Cc0];MatrixForm[A]]

{nx,ny,nz}=Eigenvectors[A][[3]]

sold=Solve [{nx,ny,nz}.Total[XYZ]-𝒹 ℕ==0,𝒹]//Flatten

Consequently the coefficients of surface model are

{-(nx/nz),-(ny/nz),(𝒹/nz)}/.sold

8 - 4 Solution via Numerical Gröbner Basis

solopt=PlaneTLS[XYZ]

p9=Plot3D[(γ +α 𝓍+β 𝓎)/.solopt,{𝓍,-8,6.2},{𝓎,-8,4},PlotRange->{143.2,147.5},BoxRatios->{1,1,0.5}, ClippingStyle->None,Mesh->None];

Show[{p9,p8}]

Fig. 14 The measured LiDAR point cloud

8 - 5 Error analysis

Error=Map[(((#[[3]]-α #[[1]]-β #[[2]]-γ)/(Sqrt[1+(α^2)+(β^2)])))/.solopt&,XYZ];

Min[Error]

Max[Error]

Mean[Error]

StandardDeviation[Error]

Histogram[Error]

Fig. 15 The histogram of the errors

8 - 6  Application of RANSAC

Now let us use for threshold  τ = 0.1 m

τ=0.1;

and let

ν=1/2;s=3;dataQ=XYZ;

therefore

𝒩=Round[Log[1-p]/Log[1-(((1-ν))^s)]]

Preparation of the sets of inliers for every iteration step - at the beginning they are empty sets:

inliers=Table[{},{𝒩}] (* the consensus sets, j = 1,2, ...𝒩*);

Then let us carry out the double For cycle to find the maximum consensus set

AbsoluteTiming[For[j=1,j≤𝒩,j++,
(*get three random points*)
dataS=RandomSample[dataQ,s];
(*determine model parameters*)
sols=PlaneTLS[dataS];
(* determine inliers*)
For[i=1,i≤ℕ,i++,
If[(((1/(Sqrt[1+(α^2)+(β^2)]))Abs[dataQ[[i]][[3]]-α dataQ[[i]][[1]]-β dataQ[[i]][[2]]-γ])/.sols)<τ,inliers[[j]]=Append[inliers[[j]],dataQ[[i]]]];
];];]

Select the consensus set which has the highest number of inliers

trimmeddata=inliers[[Ordering[inliers,-1]]][[1]];

Length[trimmeddata]

p10=ListPointPlot3D[trimmeddata,PlotStyle->Green]

Fig. 16 The original plane and the inliers of RANSAC  (green)

Then we can compute the parameters

soloptT=PlaneTLS[trimmeddata]

p11=Plot3D[(γ +α 𝓍+β 𝓎)/.soloptT,{𝓍,-8,6.2},{𝓎,-8,4},PlotRange->{143.2,147.5},BoxRatios->{1,1,0.5}, ClippingStyle->None,Mesh->None];

Show[{p11,p10}]

Fig. 17 The original plane and the inliers of RANSAC  (green)

8 - 7 Error analysis

Error=Map[(((#[[3]]-α #[[1]]-β #[[2]]-γ)/(Sqrt[1+(α^2)+(β^2)])))/.soloptT&,trimmeddata];

Min[Error]

Max[Error]

Mean[Error]

StandardDeviation[Error]

pHT=Histogram[Error]

Fig. 18 The histogram of the errors

8 - 8 Data preprocessing via SOM

Let us employ now a Kohonen map of size 7 × 6 code-book vectors with symmetrical neighborhood.

AbsoluteTiming[{som,fitrecord}=UnsupervisedNetFit[XYZ,42,35,SOM->{7,6},Neighbor->Symmetric];]

Fig. 19 The convergence of the SOM iteration

dataSOM=som[[1]];

p12=ListPointPlot3D[dataSOM,PlotStyle->Directive[Magenta,PointSize[0.015]]];

Show[{p11,p12}]

Fig. 20 The original plane and the code-book vectors of the SOM of 7 × 6

Applying RANSAC  to these 42 code-book vectors with threshold

τ=0.1;

NSOM=Length[dataSOM];

inliers=Table[{},{𝒩}] (* the consensus sets, j = 1,2, ...𝒩*);

AbsoluteTiming[For[j=1,j≤𝒩,j++,
(*get three random points*)
dataS=RandomSample[dataSOM,s];
(*determine model parameters*)
sols=PlaneTLS[dataS];
(* determine inliers*)
For[i=1,i≤NSOM,i++,
If[(((1/(Sqrt[1+(α^2)+(β^2)]))Abs[dataSOM[[i]][[3]]-α dataSOM[[i]][[1]]-β dataSOM[[i]][[2]]-γ])/.sols)<τ,inliers[[j]]=Append[inliers[[j]],dataSOM[[i]]]];
];];]

Select the consensus set which has the highest number of inliers

trimmeddataS=inliers[[Ordering[inliers,-1]]][[1]];

Length[trimmeddataS]

p13=ListPointPlot3D[trimmeddataS,PlotStyle->Directive[Magenta,PointSize[0.015]]];

Show[{p11,p13}]

Fig. 21 The original plane and the code-book vectors of the SOM of 7 × 6

Then we can compute the parameters

soloptS=PlaneTLS[trimmeddataS]

p14=Plot3D[(γ +α 𝓍+β 𝓎)/.soloptS,{𝓍,-8,6.2},{𝓎,-8,4},PlotRange->{143.2,147.5},BoxRatios->{1,1,0.5}, ClippingStyle->None,Mesh->None];

Show[{p11,p14}]

Fig. 22 The original plane and the plane estimated by the code-book vectors of the SOM of 7 × 6

8 - 9 Error analysis

Error=Map[(((#[[3]]-α #[[1]]-β #[[2]]-γ)/(Sqrt[1+(α^2)+(β^2)])))/.soloptS&,trimmeddata];

Min[Error]

Max[Error]

Mean[Error]

StandardDeviation[Error]

pHS=Histogram[Error,ChartStyle->Red]

Fig. 23 The histogram of the errors in case of employing data reduction via code-book vectors of the SOM of 7× 6

Show[{pHS,pHT}]

Fig. 24 The histogram of the error in case of using the total dataset (gray, see Fig.18) and the histogram of the errors in case of data preprocessing (red, see Fig. 23)

Conclusions

In this  a new algorithm for plane fitting to large cloud of data set is presented. This method is based on the Gröbner bases solution of the polynomial form of maximum likelihood equations with error definition in TLS sense. This method integrated into RANSAC algorithm is more efficient - faster and provide better fitting - than RANSAC with PCA, see Appendix.The computations carried out with artificial benchmark data as well as with real LiDAR measurements  demonstrate that employing first SOM as data preprocessing reducing the thousands of points of the cloud of data to a dozens number of  code-book vectors, and then running RANSAC with Gröbner solution, could reduce the computation time of the fitting process roughly by half.

References

Chen C C and Stamos I (2007) Range image segmentation for modeling and object detection in urban scenes, 3DIM2007.1

DalleMole V, do Rego R and Araújo A (2010) The Self-Organizing Approach for Surface Reconstruction from Unstructured Point Clouds,  in Matsopoulos G.K.(Ed.) Self-Organizing Map, In Tech, 

Diebel  J R, Thrun S and Brunig M (2006) A Bayesian method for probable surface reconstruction and decimation, ACM Transactions on Graphics (TOG). 1

Fischler M A and Bolles R C (1981) Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM, 24, pp.381-395

Huang C-M and Tseng Y-H. Plane fitting methods of Lidar point cloud, Dept. of Geomatics, National Cheng Kung Uni. Taiwan, tseng@mail.ncku.edu.tw

Hubert M , Rousseeuw P J and Van den Branden K (2005) ROBPCA: A new approach to robust principal component analysis, TECHNOMETRICS, Vol. 47. No. 1. pp.64 - 79.

Kohonen T (1998) The self-organizing map, Neurocomputing, Vol. 21,pp.1-6.

Lakaemper R and Latecki L J (2006) Extended EM for Planar Approximation  of 3D Data, IEEE Int. Conf. on Robotics and Automation (ICRa), Orlando, Florida, May 2006.

Lukács G, Martin R and Marshall D (1998) Faithful Leat-Squares Fitting of Spheres, cylinders, Cones and Tori for Reliable Segmentation, in Burkhardt H and Neumann B (Eds.) Computer Vision, ECCV ‘98, Vol I, LNCS 1406, pp.671-686, Springer-Verlag Berlin Heidelberg.

Mitra N J and Nguyen (2003) SoCG’03, June 8-10, San Diego, California, USA, ACM 1-58113-663-3/03/0006,pp. 322-328.

Nievergelt Y (2000)  A tutorial history of least square with applications to astronomy and geodesy, J. Comp. and Appl. Mathematics,121, pp. 37-72.

Nurunnabi A, Belton D and West G. (2012) Diagnostic -Robust Statistical Analysis for Local Surface Fitting in 3D Point Cloud Data, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol.I-3,2012 XXII ISPRS Congress, 25 Aug. 2012, Melbourne, Australia, pp.269-275

Paláncz B (2014) Fitting Data with Different Error Models, The Mathematica Journal (accepted for publication, to be published in April)

Poppinga J, Vaskevicius N, Birk A and Pathak K (2006)Fast Plane Detection and Polygonalization in noisy 3D Range Images, Int. Conf. on Intelligent Robots and Systems (IROS), Nice, France, IEEE Press 2008

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

Russeeuw P J and Van Driessen K (1999)  A Fast Algorithm for the Minimum Covariance Determinant Estimator, TECHNOMETRICS, Vol. 41. No. 3. pp. 212-223.

Stathas D, Arabatzi O, Dogouris S, Piniotis G, Tsini D and Tsinis D (2003) New Monitoring Techniques on the Determination of Structure Deformation, Proc. 11th FIG Symp. on Deformation Measurements, Santorini, Greece. 

 Zuliani M (2012) RANSAC for Dummies, vision.ece.ucsb.edu/~zuliani

Yang M Y and Förtsner W (2010) Plane Detection in Point Cloud Data,TR-IGG-P-2010-01, Techn.Report. Nr.1, Dept. of Photogrammetry Inst. of Geodesy and Geo-information, Uni., Bonn, Germany.

Yaniv Z (2010) Random Sample Consensus (RANSAC) Algorithm, A Generic Implementation, Georgetown University Medical Center, Washington, DC, USA, zivy@isis.georgetown.edu

Appendix

A - 1 The Mathematica function for the PCA Solution

PlanePCA[XYZ_]:=Module[{𝒩,co,Cco,nx,ny,nz,𝒹,sol𝒹},
𝒩=Length[XYZ];
co=Mean[XYZ];
Cco=Map[#-co&,XYZ];
{nx,ny,nz}=Eigenvectors[Covariance[Cco]][[3]];
sol𝒹=Solve [{nx,ny,nz}.Total[XYZ]-𝒹 𝒩==0,𝒹]//Flatten;
{α->-(nx/nz),β->-(ny/nz),γ->(𝒹/nz)}/.sol𝒹]

A - 2 Application of PCA in RANSAC

Let us consider the synthetic data, see Section 5.2,

τ=1;

s=3;

ν=1/2;

𝒩=Round[Log[1-p]/Log[1-(((1-ν))^s)]]

Preparation of the sets of inliers for every iteration step - at the beginning they are empty sets:

inliers=Table[{},{𝒩}] (* the consensus sets, j = 1,2, ...𝒩*);

AbsoluteTiming[For[j=1,j≤𝒩,j++,
(*get three random points*)
dataS=RandomSample[XYZC,s];
(*determine model parameters*)
sols=PlanePCA[dataS];
(* determine inliers*)
For[i=1,i≤ℕℂ,i++,
If[(((1/(Sqrt[1+(α^2)+(β^2)]))Abs[XYZC[[i]][[3]]-α XYZC[[i]][[1]]-β XYZC[[i]][[2]]-γ])/.sols)<τ,inliers[[j]]=Append[inliers[[j]],XYZC[[i]]]];
];];]

Select the consensus set which has the highest number of inliers

trimmeddata=inliers[[Ordering[inliers,-1]]][[1]];

Length[trimmeddata]

p15=ListPointPlot3D[trimmeddata,PlotStyle->Green];

Show[{p1,p15}]

Length

Fig. 25 The original plane and the inliers of RANSAC  (green)

Then we can compute the parameters

PlanePCA[trimmeddata]

This result illustrates that integrating PCA into RANSAC is less efficient than RANSAC with the integrated Gröbner basis solution.