Stable Algebraic Surfaces for 3D Object ... - Semantic Scholar
J Math Imaging Vis (2008) 32: 127–137 DOI 10.1007/s10851-008-0092-3
Stable Algebraic Surfaces for 3D Object Representation Turker Sahin · Mustafa Unel
Published online: 29 April 2008 © Springer Science+Business Media, LLC 2008
Abstract Linear fitting techniques by implicit algebraic models usually suffer from global stability problems. Ridge regression regularization can be used to improve the stability of algebraic surface fits. In this paper a Euclidean Invariant 3D ridge regression matrix is developed and applied to a particular linear algebraic surface fitting method. Utilization of such a regularization in fitting process makes it possible to globally stabilize 3D object fits with surfaces of any degree. Robustness to noise and moderate levels of occlusion has also been enhanced significantly. Experimental results are presented to verify the improvements in global stability and robustness of the resulting fits. Keywords Algebraic surfaces · Implicit polynomials · Fitting · Stability · Ridge regression
1 Introduction Representing objects in intensity and range images is a central problem in computer vision, and there has been a significant body of work devoted to 2D and 3D object representation. An important method among these is algebraic curve and surface representations, which model objects by fitting implicit polynomials (IPs) to their boundary data [1–6]. Algebraic curves/surfaces model objects by contributions from
T. Sahin Department of Computer Engineering, Gebze Institute of Technology, Gebze, Kocaeli, Turkey e-mail: [email protected] M. Unel () Faculty of Engineering and Natural Sciences, Sabanci University, Tuzla, Istanbul, Turkey e-mail: [email protected]
all points in their data sets, and result in global representations. There are also other techniques in the literature, such as parametric curves/surfaces, which model data by local polynomial patches [10, 11]; Fourier descriptors that represent shapes based on Fourier theory methods according to polar or arc-length parameterizations [13]; superquadrics which fit lower degree surfaces to sections of complicated objects of usually ordered symmetry [12]; level set methods representing objects by propagating fronts converging to their boundary data [14, 15]; and radial basis type neural network classifiers [16, 17]. Using radial basis functions, very complicated objects such as human faces can be represented with great accuracy, far more than with polynomials, whose global nature makes it almost impossible to use for fitting shapes with fine details. So why use polynomial curves/surfaces at all? We think that polynomial models have important advantages which make them preferable in many computer vision tasks. Many natural and man-made objects, which are too complicated to be represented by a simple ellipse/ellipsoid, can well be approximated by a higher degree curve/surface. These models are robust versus noise and reasonable levels of missing data. They are also very suitable for point classification problems and curve-surface intersection calculations, as this involves only a simple polynomial evaluation. New methods which enable more optimal data set representations have been emerging [5–7], making the results of these calculations more accurate than before. Algebraic curve and surface models are preferred over Fourier Descriptors as they are able to represent non-star, occluded shapes as well as unordered data; and over point based modelling approaches that suffer from noise, missing data and have significant dependence on accurate point correspondences [24]. Moreover, the geometric invariants of 2D algebraic models have been derived [8, 9], which have been
128
applied to important vision problems such as efficient pose estimation [22] and motion estimation and tracking of objects [23]. One drawback of implicit algebraic models has been considered as their containing high degree terms, making them prone to instability. Some nonlinear techniques [2, 5, 6, 19] have been attempted for IP modelling, but like any other nonlinear optimization they need a good initial guess and they were not computationally efficient. Moreover, reference [25] explains why these nonlinear methods don’t provide satisfactory fits in all cases and the need for further improvement. Thus, linear approaches such as [20, 21, 25] are preferred. The work in [25] brings substantial improvement over those previous approaches not only by providing higher quality fits, but also significantly reducing the fitting time. However, it is limited to 2D curves only. In this paper we present 3D ridge regression (RR) regularization as an effective method for overcoming global instability problem of a particular fitting technique, the 3L algorithm [21]. We first present an alternative level set generation scheme for 3L algorithm, which is not properly addressed in [21], to improve local stability and then derive a Euclidean invariant RR matrix to globally stabilize the resulting surfaces. RR for improving the global stability of algebraic curves first appeared in [25]; however, as mentioned above that solution was limited to 2D curve fitting problems. This work extends RR regularization to 3D models, which helps in fitting stable and accurate algebraic surfaces of virtually any degree to complex data. Most of the objects which are considered in our work are not simple surfaces. Some of them have complicated structures like holes etc. in their data sets. Proposed method takes these structures into account rather than avoiding them. A first draft of this paper appears in [26]. The structure of this text is as follows: in Sect. 2 we give a brief introduction to the implicit surface models followed by the general implicit polynomial fitting problem and introduction to linear fitting techniques in Sects. 3 and 4. In Sect. 5 we explain the 3L fitting method, and introduce a new level set generation technique we have employed for 3D data. Section 6 is on the 3D ridge regression method in detail. In Sect. 7 we present experimental results and discussions followed by conclusions in Sect. 8.
2 Implicit Surface Models 2.1 Algebraic Surfaces IP fitting techniques model 3D object data by nth degree algebraic surfaces with c = (n + 1)(n + 2)(n + 3)/6 coeffi-
J Math Imaging Vis (2008) 32: 127–137
cients of the form:
Fn (x, y, z) =
aij k x i y j zk
0≤i,j,k;i+j +k≤n
= a000 + a100 x + a010 y + a001 z + · · · H0
H1 (x,y,z)
+ an00 x + · · · + a00n zn n
Hn (x,y,z)
T = 1 x . . . zn a000 a100 . . . a00n mT
= mT a
a
(2.1)
where m is the vector of monomials, a is the coefficient vector and each ternary form Hr (x, y, z) is a homogeneous polynomial of degree r in x, y and z. 2.2 Data Set Normalization Data set normalization is an important technique, which first centers the mean point of data at origin and then scales the data to be as near the unit circle as possible. By this way the difference between the values of higher and lower degree monomial terms are reduced to reasonable levels, such that the resulting IP’s are more stable. We have applied the radial distance normalization, which is a linear process that improves the results obtained using other data normalizations such as [25]. Accordingly the center of mass of original data, (xi , yi , zi )N i=1 , is shifted to origin by: (xˆi , yˆi , zˆ i ) = (xi , yi , zi ) − μ
(2.2)
where μ = N1 N i=1 (xi , yi , zi ) is the mean point of original points. Then the normalized data is calculated by scaling this result according to: (xˆni , yˆni , zˆ ni ) =
(xˆi , yˆi , zˆ i ) s
2 2 2 1/2 is the average radial where s = N1 N i=1 (xˆ i + yˆi + zˆ i ) distance of the centralized data.
3 General Implicit Algebraic Model Fitting Problem The problem of representing a data set or some of its sections by an algebraic curve or surface is referred to as the polynomial fitting problem. Accordingly a cost function like the following is minimized J=
N 1 Dist2 (Pi , Z(Fn )) N i=1
(3.3)
J Math Imaging Vis (2008) 32: 127–137
where Dist(Pi , Z(F )) is the distance of the data point Pi = (xi , yi , zi )N i=1 to the zero set Z(Fn ) = {(x, y, z) : Fn (x, y, z) = 0} of the IP model Fn (x, y, z). There are two main categories of IP fitting techniques: the nonlinear methods, which use the true geometric distance or Euclidean approximation, and linear methods based on the algebraic distance. The nonlinear approaches are invariant to transformations in the Euclidean space and are not biased [3, 4]. However, except for very basic shapes, there is no available closed form expression for the Euclidean distance. Therefore, time consuming iterative optimization procedures must be carried out. For this reason linear techniques are preferred, which formulate (3.3) as a linear least squares problem resulting in faster computation times [20, 21, 25].
4 Linear Algebraic Surface Fitting Techniques The simplest linear surface fitting approach is the Least Squares method, which minimizes the algebraic distance over the data sets (Pi )N i=1 with an additional constraint for avoiding the trivial solution. The solution of this optimization problem is the unit eigenvector a associated with the smallest eigenvalue of Sa = λa, where S = M T M is the scatter matrix of monomials. Here M is the monomials matrix composed of the monomial vectors mi of data points such that:
T M = m1 m2 . . . mN (4.4) This technique is affine invariant [20], but unfortunately is very unstable, which makes it unsuitable for any practical application. Thus more satisfactory linear fitting methods are preferred over the least squares technique. These methods make use of data set normalization and local continuity measures, both of which improve accuracy and stability. One example is Gradient One algorithm, which applies gradient control by incorporating scatter matrices from local normals and tangents in addition to the original scatter matrix of monomials [25]. In this method, the normal and tangential components are imposed to be 1 and 0 for the resulting surface gradient to be as near to unity as possible. However this technique is more sensitive to changes in data normalization and in its tuning parameter than 3L method. Also the practical Euclidean invariance properties of Gradient One are weaker, which make 3L a more suitable candidate for vision tasks such as motion estimation and object tracking.
5 The 3L Fitting Method The 3L technique incorporates control for local continuity by placing two additional data layers at a distance ±ε out-
129
side and inside the original data as the optimization constraint [21]. Thus the resulting fit is constrained to be no more than ±ε away from the actual points. The solution of this problem is: a = M †b
(5.5)
where M † = (M T M)−1 M T is the pseudo-inverse of the matrix of monomials, M. Similar to the other linear approaches, this method achieves local stability about the data, but unfortunately cannot provide a solution for the global instability problem. However, it is less sensitive to changes in its tuning parameter, ε, and in applied data normalization. For this reason we have applied ridge regression regularization to this method, which gives better overall results. A crucial aspect of this method is level set generation. To achieve better fits, we have triangulated the 3D data sets, and then calculated the point normals as the average of the neighboring triangular faces for forming the inner and outer level sets. The related triangulation and level set formation methods are summarized in the next sections. 5.1 Triangulation The crust triangulation algorithm of [18] has been applied for forming the triangular tessellations of object data. This method uses Voronoi diagram and Delaunay triangulation approaches to find the topological connection of data points. The algorithm first computes the Voronoi diagram of the sample set. For this reason it forms syntectic cubic volumes called Voronoi cells for every data point. It then selects the two furthest vertices of the cells from each other, which are referred to as poles. One pole has to be inside the data surface while the other must be outside. Next the method computes the Delaunay triangulation of the union set of the original sample points and their poles. Finally, only the triangles whose all three vertices are sample points are selected as the triangle patches. The algorithm returns the list of triangular faces and the neighboring vertices for data points. This method may not be able to detect some triangles in the cases of very scarce or noisy data. However, it gives reasonable results with sufficiently dense and evenly spaced points. The object data we have employed and corresponding triangular meshes are in Figs. 1 and 2 respectively to depict the results of the employed algorithm. 5.2 Level Set Generation The triangular faces obtained from the crust algorithm are used to calculate the local normals of data points for obtaining the virtual inner and outer data layers correctly. For sufficiently accurate operation of the 3L technique, the neighboring triangular patch normals of each data point should be
130
J Math Imaging Vis (2008) 32: 127–137
Fig. 1 The 3D data used in our fits. Examples include a torus, a phone receiver, an eight shape, a sphere, a duck, a heart, an apple, a cube, a face, a twisted object, a golf stick handle, a head, a ball joint, a hip and a patella bone
Fig. 2 Triangulation of the sample data by Crust Algorithm
correctly aligned at each point to get a good overall representation of the level sets. Hence after the neighboring faces of a point have been obtained from the generated mesh, the directions of their normals are rearranged to satisfy the right hand rule for local alignment as in Fig. 3(a). Then these normals are summed and normalized to give the unit normal at the data point. There may be problems, which prevent accurate local alignment, when object data has very complex topology or is noisy. Such effects cause the neighboring triangle arrays (resulting from the algorithm) not to be of purely planar topology. Thus to avoid possible errors, these exceptional ill-structured patch arrays are checked beforehand and eliminated by the application of a single normal, which is perpendicular to the plane fit of the point and the neighboring vertices. The resulting point normals may be globally unaligned, that is they may not be all pointing to
outside the data or vice versa. Hence they should also be globally checked for having correct directions. To achieve this we have applied the local continuity of gradient principle to the neighborhood of each data point. Accordingly, the direction differences of adjacent point normals should have angles less than 90◦ . Starting with an arbitrary point, its normal can be used as a reference to align each neighboring normal by a simple dot product check. To be more robust versus the discrepancies in the local continuity of the data set, a margin can be set aside so that the normal difference angles can be kept at a significantly lower degree than acute angle for the checking process. This procedure is repeated for the adjacent point as reference, until all normals have been readjusted accurately. This global alignment is important for placing the synthetic layers correctly. In cases of correct global alignment, all points of the inner and outer
J Math Imaging Vis (2008) 32: 127–137
(a) The point normals Ni obtained from face normals NFj by right hand rule
(b) Cross-section of the level sets of a 3D free form object Fig. 3 Local alignment of surface point normals and the resulting level sets
level sets should stay ±ε inside and outside the actual data as in Fig. 3(b) of the 2D silhouette of the 3 layers of a 3D free form object. Another important factor for satisfactory 3L fits is the choice of the ε parameter. This selection is important for correct alignment of the 2 synthetic level sets, as depicted in Fig. 3(b), which prevents the least squares problem of (5.5) from being ill-posed. ε should be selected not too small for flexibility and not too high for the correct alignment of inner and outer layers at the detailed sections of the data sets. For the radial distance normalization that we have applied, choices of ε, which give better results are typically between 0.01–0.05. It is not recommended to exceed these values, as otherwise the detailed sections may not be preserved accurately.
6 Ridge Regression Regularization for Stable Algebraic Surface Fits Linear fitting methods in general have difficulty in achieving global stability. Especially in modelling complicated shapes, many extra open surface pieces are generated, which make the resulting fits nearly impossible to be interpreted. An important reason for global instability is the collinearity of modelled data, which cause the M T M matrix of the associated least squares problem to be nearly singular with some
131
eigenvalues negligible compared to the others. Such eigenvalues contribute very little to the overall shape of the fit, but are one of the main causes of instability. The ratio of the largest eigenvalue to the smallest one is referred to as the condition number and is a measure of the singularity of matrices. The work in [25] employs RR regularization as a way to counteract on such singularity effects. Fitting results are impressive, but limited to 2D. Our work will handle ridge regression regularization of algebraic surfaces in 3D. RR improves the condition numbers and the stability of implicit polynomials by adding a new κD term to M T M data scatter matrices. Accordingly, (5.5) is modified as: aκ = (M T M + κD)−1 M T b
(6.6)
where κ is the RR parameter and D is a positive definite, symmetric matrix, which will be selected as diagonal to facilitate easier tunability. For RR to achieve its purpose, κ should be increased from 0 to higher values until a stable closed bounded surface is obtained. By this way the small eigenvalues are significantly increased, while the larger eigenvalues that contribute more to the overall shape of the fit are virtually unaffected, hence reducing the condition number to reasonable levels. Thus provided that the κ term is not selected too high, unstable fits can be stabilized with no loss of shape. 6.1 The Euclidean Invariance of Ridge Regression Regularization The other important part of the RR term is the D matrix of regularization, which is selected such that Euclidean invariance is not violated while at the same time the shape representation is preserved. Accordingly, when an arbitrary Euclidean transformation E such as: ⎡ ⎤ ⎡ ⎤ ⎤ x ⎡ x R t ⎢ ⎥ ⎢y ⎥ ⎢ ⎥=⎣ ⎦ ⎢y ⎥ (6.7) ⎣z ⎦ ⎣z⎦ 0 0 0 1 1 1 E
is applied to data points where R is the 3-D rotational part and t = [tx ty tz ]T is the translational vector, the original and transformed surface monomial vectors m(x, y, z) and m (x , y , z ) will be related by a c × c square Eˆ matrix of block diagonal rotational matrices part and corresponding translational sections: ˆ m = Em
(6.8)
As the initial and transformed surfaces should be equivalent under Euclidean transformations, we can relate the two coefficient vectors by mT a = mT a
(6.9)
132
J Math Imaging Vis (2008) 32: 127–137
As Eˆ is composed of block diagonal rotational matrix sections, which all have inverses on their own, it is an invertible matrix. Thus applying (6.8) to (6.9) we get the following relation between the original and transformed coefficient vectors.
QT 1/2 DT 1/2 QT = T 1/2 DT 1/2
a = (Eˆ T )−1 a
D = T −1
(6.10)
For Euclidean invariance, we should be able to derive the relation (6.6) for the transformed coefficient vector, a , as well. Thus using (6.10) and (6.6): a = (Eˆ T )−1 (M T M + κD)−1 M T b
(6.11)
Here substituting Eˆ −1 Eˆ = I after the inverted monomial matrix, the following form is obtained: ˆ T M)Eˆ T + κ ED ˆ Eˆ T −1 EM ˆ Tb a = E(M
(6.12)
By (4.4) and (6.8), this equation can be rearranged to: ˆ Eˆ T )−1 M T b a = (M T M + κ ED
(6.13)
which should preserve Euclidean invariance when the RR part is equal to κD. This is satisfied for: ˆ Eˆ T = D ED
(6.15)
where Q is an arbitrary block diagonal rotation matrix with blocks corresponding to the ternary forms of (2.1). T is the diagonal matrix of trinomial coefficients, whose elements are: (i + j + k)! Tll = i!j !k!
(6.16)
where i, j, k are incremented as follows: (i, j, k) = (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), . . . , (n, 0, 0), . . . , (0, n, 0), . . . , (0, 0, n) and the index l of the diagonal terms are incremented from 1 to the number of IP coefficients, according to the above variation of i, j, k exponents by the formula: (j + k)(j + k + 1) 2 (i + j + k)(i + j + k + 1)(i + j + k + 2) +1 + 6
l=k+
As Q is a rotation matrix, i.e. QQT = I , a Euclidean invariant D matrix can be derived as: (6.17)
which is the inverse of the trinomial coefficients matrix. This selection of D does not maintain the local shape representation, but this can be fixed by multiplying it by the weighted sum of the diagonal elements of the monomials matrix M T M. Such terms correspond to the square of the distance of the normalized data points to the origin, thus preserve Euclidean invariance and improve the local shape properties at the same time. The terms of the resulting diagonal D matrix are: Dll = βi+j +k
1 Tll
(6.18)
Here the term βi+j +k is contributed from the data set points (Pt )N t=1 formulated as: βi+j +k =
N (xt2 + yt2 + zt2 )i+j +k t=1
(6.14)
If (6.7) is reduced to a pure rotation, which is the case for normalized data, Eˆ can be represented by the similarity transformation: Eˆ = T −1/2 QT 1/2
Substituting (6.15) into (6.14) implies that
As a simple example the D matrix that should work on quadric surface monomial matrices should be: 1 1 1 D = diag β0 , β1 , β1 , β1 , β2 , β2 , β2 , β2 , β2 , β2 2 2 2 2 2 2 where β0 = N , β1 = N t=1 (xt + yt + zt ) and β2 = N 2 2 2 2 t=1 (xt + yt + zt ) . 6.2 The Closed-Boundedness of Ridge Regression in 3D The RR regularized 3L method (3LRR) shows the important boundedness property of RR fits as the κ term is increased towards infinity in value. Accordingly provided that the modelled object data sets are closed (i.e. with no significant gaps or occlusion) and the employed IPs are of even degree, the resulting fits always converge to a bounded surface. This behavior can be proved as follows: From (6.6), the limit of the 3LRR IP vector as κ → ∞ is: aij k =
1 −1 T D M b κ
(6.19)
Here the D matrix is according to the selection in (6.18). Thus we concentrate on the M T b part, composed of the sums of monomials in ascending order, such that: T i j k i j k M b= ε xy z − xy z (6.20) +ε
−ε
c×1
J Math Imaging Vis (2008) 32: 127–137
with i, j, k in the order of (6.16) and ±ε represent the sums of monomials in the outer and inner data layers, respectively. This summation is in the form of directional derivative at each data point in the outward direction from the data set surface, i.e. (x i y j zk ) (x i y j zk )+ε − (x i y j zk )−ε = s 2ε
MT b = S
T N · ∇(x i y j zk ) c×1 dS
diag{1/Dll }
=
(i + j + k)! βi+j +k i!j !k!
(6.22)
(i + j + k)! βi+j +k i!j !k!
(6.26)
dxdydz
i>1,j >1,k>1; i+j +k=n
n(n − 1) = βn
i>1,j >1,k>1; i+j +k=n
xˆ 2
(n − 2)! (i − 2)!j !k!
ˆ j (zˆz)n−i−j × (x x) ˆ i−2 (y y) (n − 2)! (x x) ˆ i (y y) ˆ j −2 (zˆz)n−i−j i!(j − 2)!k!
(n − 2)! (x x) ˆ i (y y) + zˆ ˆ j (zˆz)n−i−j −2 i!j !(k − 2)! 2
dxdydz (6.27)
The three integrand terms are in the form of the trinomial expansion formula:
V
i(i − 1)x i−2 y j zk
(x + y + z)n =
V
aij k xˆ i yˆ j zˆ k = 0
+ zˆ 2
+ yˆ 2
n i=0
(6.24)
with V the region bounded by the surface S. From (2.1) the zero set of the IP vector a∞ is obtained to be: T a∞ m=
(i + j + k)! j (j − 1)(x x) ˆ i (y y) ˆ j −2 (zˆz)k βi+j +k i!j !k!
As i + j + k = n, (6.26) can be re-expressed as: aij k xˆ i yˆ j zˆ k
(6.23)
∇ 2 (x i y j zk )dxdydz
+ j (j − 1)x i y j −2 zk + k(k − 1)x i y j zk−2 dxdydz
+ yˆ 2
ˆ j (zˆz)k−2 ×k(k − 1)(x x) ˆ i (y y)
(i + j + k)! βi+j +k i!j !k!
ˆ j (zˆz)k ×i(i − 1)(x x) ˆ i−2 (y y)
MT b
(i + j + k)! βi+j +k i!j !k!
xˆ 2
i>1,j >1,k>1; i+j +k=n
scaled by a constant term. Applying the divergence theorem and utilizing ∇ · ∇f = ∇ 2 f where ∇ 2 represents the Laplacian operator, we get: aij k =
=
(6.21)
scaled with a constant, where N T is the local normal at each data point and ∇(x i y j zk ) is the gradient of the monomials. These terms are integrated over the complete surface of the data, S. Hence (6.19) can be brought into the form:
T (i + j + k)! aij k = N · ∇(x i y j zk ) c×1 dS βi+j +k i!j !k! S
aij k xˆ i yˆ j zˆ k
i>1,j >1,k>1; i+j +k=n
with s being the vector in the direction of the derivative. As the directional derivative of the monomial terms is related to the gradient in the outward direction (i.e. the normal) to the surface, we also express the M T b terms as:
133
(6.25)
0≤i+j +k≤n
where mT = (xˆ i yˆ j zˆ k ) is the monomial vector. If we can show that the leading form of (6.25) is strictly positive, we complete the proof of the boundedness of 3LRR fits as the RR term tends to infinity. By embedding (6.24) to (6.25) and concentrating on the leading form, we obtain:
n! x i y j z(n−i−j ) i!j !(n − i − j )!
(6.28)
therefore the leading form can be simplified to: aij k xˆ i yˆ j zˆ k i>1,j >1,k>1; i+j +k=n
=
n(n − 1)(xˆ 2 + yˆ 2 + zˆ 2 ) βn × (x xˆ + y yˆ + zˆz)n−2 dxdydz
(6.29)
Here if n is even, the integrand and all the other terms become positive. Hence the leading form of the 3LRR IP equation in 3D is strictly positive, imposing the resulting fits to converge to an ellipsoid like closed bounded surface. The proof that we have just presented for the closedboundedness of the resulting fits ensures the existence
134
J Math Imaging Vis (2008) 32: 127–137
Fig. 4 Least squares fits of 3D sample data: the torus, the phone receiver, eight, cube and the heart
Fig. 5 Stable 3L surface fits of the torus, the phone, the cube, the heart and the sphere objects
of some κ values which can stabilize the fits. In fact, as shown in the experimental section, a κ value in the range 10−4 –10−3 stabilizes the resulting fits in virtually all cases by forcing the unbounded/disconnected components to move away from the data and finally being eliminated from the fits. 6.3 Suitable Selection of κ Term for Shape Representation and Recognition The RR term κ reduces the variability of parameters and makes the IP model less sensitive to small perturbations in data. It also has the effect of smoothing the object surface as κ is increased in value. However, for an as efficient as possible shape modelling, the IP model should not be smoothed to much, and yet all extra open components should be eliminated. Thus starting from zero, the κ parameter can be iteratively increased to higher values, until the surface is closed bounded. The boundedness of the surface can be checked according to the IP leading form positive-definiteness principle introduced in Sect. 6.2. Accordingly, the smallest value of κ obeying this test should be the most suitable selection for shape modelling. For object recognition purposes, the main objective should be to minimize the error between the actual 3D data and the fitted surface polynomial. However, there is no single optimal κ value that obeys this criterion and κ for best recognition may vary between different sets of data; however, this value should be higher than the threshold value suitable for shape modelling. 7 Experimental Results and Discussion The conventional linear fitting approaches in 3D is rather limited to modelling objects of simpler topologies. Among these methods the least squares algorithm cannot provide stable and accurate fits for any data except the torus as depicted in Fig. 4. The stability of the torus is mainly the result of the perfect symmetry of this object, and is still limited to a 4th degree fit. When slightly more complicated shapes such
Fig. 6 Unstable 3L surface fits, which include an eight shape, a duck, the man face, the woman face, the golf, ball joint, hip, the twisted object and the patella in this order
as the phone and the eight shaped object are employed, accuracy of fits is lost significantly as in the same figure. Finally when more complicated structures like the cube and heart are modelled, this method fails completely. The 3L fitting is adequate for simple data shapes that can be modelled by lower degree fits. More data examples than least squares fit like the phone, the cube, the heart, and the sphere can be stably and accurately fit as depicted in Fig. 5. However, for practical applications, more complex objects have to be fit with stable IP surfaces. The remainder of our sample objects are in this category, for which the 3L method fails to give good IP models. The related fits are depicted in Fig. 6, all with extra surface pieces.
J Math Imaging Vis (2008) 32: 127–137
135
Fig. 7 The stabilization of the Eight, Face and twisted objects by RR regularization
Fig. 8 The stabilized surface fits of objects in Fig. 6 by RR
RR works to achieve global stability as the κ parameter is increased from zero to sufficiently higher values. Accordingly, it gives stable fits for the objects that the 3L method could not manage to model adequately. Some examples are presented in Figs. 7 and 8. In Fig. 7, the stabilization of surface fits of the eight shaped object, the face and the twisted object are depicted. Here the eight and the face objects are modelled by 6th degree surfaces, while the twisted object has an 8th degree fit. The leftmost sub-figures are the 3L only subplots with κ = 0. The following intermediate subplots are for higher values of κ to depict how RR moves the extra surfaces away from the zero set and eliminates
them. The stabilizations are achieved for the face and the twisted object at κ = 10−5 and for the eight shaped object at κ = 10−4 . These are typical values for the RR parameter and the fits have very good shape representation, although lower degree IP’s are employed. As in zoomed rightmost sub-figures, the eight and twisted objects are modelled with all their major details, and the face surface contains critical sections like the nose and the forehead. All the objects, 3L method could not stabilize in Fig. 6 can be modelled accurately by RR regularization as presented in Fig. 8, which verify the associated global stabilization properties. Moreover RR is robust versus perturbations and structural problems in data sets. As in Fig. 9, despite a significant section of the objects such as 15% of the golf, 12% of the duck, and 10% of the face and the eight, have been discarded, RR gives stable and accurate surface fits for all data. The only difference is that the RR parameter was increased slightly more to achieve stabilization; however, this virtually unaffects the shape of the fits. Robustness versus noise property of this method is also strong. The noise tests are presented in Fig. 10 for the eight, patella, twisted object and the face. The objects of smoother data sets with less fine detail are robust versus noise levels of 5% of object sizes, as with σ 2 = 0.0025 for the eight and σ 2 = 0.002 for the patella. It should be noted that the employed data have been normalized to unit ball and the 3L parameter has been set to ε = 0.05. Thus the inner and outer synthetic data levels are at ±5% of the object size away from the actual data surface, and at higher noise levels, they are indistinguishably independent of RR. As a result we can conclude that RR can restore the stability of least squares type fitting meth-
136
J Math Imaging Vis (2008) 32: 127–137
Fig. 9 The face, eight, golf and the duck RR fits under missing data conditions
Fig. 10 RR fits for the eight, patella, twisted and face data subjected to noise levels of variances σ 2 = 0.0025, σ 2 = 0.002, σ 2 = 0.0009 and σ 2 = 0.0004, respectively
ods while preserving accuracy up to the levels permitted by these techniques. The twisted and the face objects with more fine detail can also be fit stably and accurately up to noise levels of 0.0009 and σ 2 = 0.0004, respectively. At higher noise levels, stability is still ensured although some object fine detail cannot be preserved satisfactorily.
of the modelled objects as much as possible. Thus globally stable and accurate fits of higher degrees have been obtained for range data with fine details. The resulting fits are also robust versus noise levels of up to 5% of object size and 15% missing data. As a result this fitting technique should be suitable for modelling 3D objects with algebraic surfaces of any practical degree, and thus should find important applications in 3D modelling and computer vision.
8 Conclusions In this paper a Euclidean invariant 3D ridge regression matrix is developed and applied for the regularization of 3L fitting technique. This approach results in closed bounded algebraic surfaces, while preserving the local shape properties
References 1. Pilu, M., Fitzgibbon, A., Fisher, R.: Ellipse specific direct least squares fitting. In: Proc. IEEE International Conference on Image Processing, Lausanne, Switzerland, September 1996
J Math Imaging Vis (2008) 32: 127–137 2. Keren, D., Cooper, D., Subrahmonia, J.: Describing complicated objects by implicit polynomials. IEEE Trans. Pattern Anal. Mach. Intell. 16, 38–53 (1994) 3. Faber, P., Fischer, R.B.: Euclidean fitting revisisted. In: 4th IWVF, pp. 165–175 (2001) 4. Faber, P., Fischer, R.B.: Pros and cons of Euclidean fitting. In: Proc. Annual German Symposium for Pattern Recognition, DAGM01, Munich. LNCS, vol. 2191, pp. 414–420. Springer, Berlin (2001) 5. Keren, D., Gotsman, C.: Fitting curves and surfaces with constrained implicit polynomials. IEEE Trans. Pattern Anal. Mach. Intell. 21(1), 31–41 (1999) 6. Keren, D.: Topologically faithful fitting of simple closed curves. IEEE Trans. Pattern Anal. Mach. Intell. 26(1), 118–123 (2004) 7. Redding, N.J.: Implicit polynomials, orthogonal distance regression, and the closest point on a curve. IEEE Trans. Pattern Anal. Mach. Intell. 22(2), 191–199 (2000) 8. Unel, M., Wolovich, W.A.: On the construction of complete sets of geometric invariants for algebraic curves. Adv. Appl. Math. 24, 65–87 (2000) 9. Unel, M., Wolovich, W.A.: A new representation for quartic curves and complete sets of geometric invariants. Int. J. Pattern Recognit. Artif. Intell. 13(8), 1137–1149 (1999) 10. Huang, Z., Cohen, F.S.: Affine-invariant B-spline moments for curve matching. IEEE Trans. Image Process. 5(10), 1473–1480 (1996) 11. Cohen, F.S., Wang, J.-Y.: Modeling image curve using invariant 3-D object curve models—a path to 3-D recognition and shape estimation from image contours. IEEE Trans. Pattern Anal. Mach. Intell. 16(1), 1–12 (1994) 12. Terzopoulos, D., Metaxas, D.: Dynamic 3D models with local and global deformations: deformable superquadrics. IEEE Trans. Pattern Anal. Mach. Intell. 13(7), 703–714 (1991) 13. Wu, M.F., Shen, H.T.: Representation of 3D surfaces by twovariable Forurier descriptors. IEEE Trans. Pattern Anal. Mach. Intell. 20(8), 583–588 (0000) 14. Osher, S., Sethian, J.A.: Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. J. Comput. Phys. 79, 12–49 (1988) 15. Malladi, R., Sethian, J.A., Vemuri, B.C.: Shape modelling with front propagation: a level set approach. IEEE Trans. Pattern Anal. Mach. Intell. 17(2), 158–175 (1995) 16. Dinh, H.Q., Turk, G., Slabaugh, G.: Reconstructing surfaces by volumetric regularization using radial basis functions. IEEE Trans. Pattern Anal. Mach. Intell. 24(10), 1358–1371 (2002) 17. Zhang, M., Ma, L., Zeng, X., Wang, Y.H.: Image based 3D face modelling. In: Proceedings of the IEEE International Conference on Computer Graphics, Imaging and Visualization, CGIV04 (2004) 18. Amenta, N., Choi, S., Kolluri, R.K.: The power crust. In: Proceeding of the 6th ACM Symposium on Solid Modelling and Applications, pp. 249–266 (2001) 19. Taubin, G., et al.: Parametrized families of polynomials for bounded algebraic curve and surface fitting. IEEE Trans. Pattern Anal. Mach. Intell. 16(3), 287–303 (1994) 20. Deepak, K., Mundy, J.: Fitting affine invariant conics to curves. In: Geometric Invariance in Machine Vision. MIT Press, Cambridge (1992) 21. Blane, M., Lei, Z., et al.: The 3L algorithm for fitting implicit polynomial curves and surfaces to data. IEEE Trans. Pattern Anal. Mach. Intell. 22(3), 298–313 (2000)
137 22. Unel, M., Wolovich, W.A.: Pose estimation and object identification using complex algebraic representations. Pattern Anal. Appl. J. 1(3), 178–188 (1998) 23. Sahin, T., Unel, M.: An adaptive estimation method for rigid motion parameters of 2D curves. In: LNCS, vol. 3212, pp. 355–362. Springer, Heidelberg (2004) 24. Ramesh, J., Kasturi, R., Schunc, B.G.: Machine Vision. McGraw Hill, New York (1995) 25. Tasdizen, T., Tarel, J.-P., Cooper, D.B.: Improving the stability of algebraic curves for applications. IEEE Trans. Image Process. 9(3), 405–416 (2000) 26. Sahin, T., Unel, M.: Fitting globally stabilized algebraic surfaces to range data. In: Proeceedings of the 10th IEEE International Conference on Computer Vision, Beijing, China, October 2005
Turker Sahin has received his B.Eng. Degree in Electronic Engineering with Optoelectronics from University College London in 1999 and his M.Sc. Degree in Control Systems from Imperial College in 2001. During 2003–2007, he was a Ph.D. candiate and research assistant in the Department of Computer Engineering at Gebze Institute of Technology, Kocaeli, Turkey where he did work on algebraic and differential geometry of curves and surfaces, visual motion estimation and mobile robotics. He obtained his Ph.D. degree in January 2008. His current research interests are in the areas of vision, control and robotics. Mustafa Unel received both the BSEE degree in 1994 and the M.S. degree in Systems and Control Engineering in 1996 from Bogazici University, Istanbul, Turkey. He received a second M.S. degree in Applied Mathematics in 1998 and the Ph.D. degree in Electrical Sciences and Computer Engineering in 1999 both from Brown University, Providence, USA. During 1999–2001, he was a postdoctoral fellow at Center for Computational Vision and Control, Yale University, New Haven, USA. In Spring 2001, he joined the Department of Computer Engineering at Gebze Institute of Technology, Turkey as an assistant professor. Since Fall 2004, he has been with the Faculty of Engineering and Natural Sciences, Sabanci University, Istanbul, Turkey where he is currently an associate professor. His research interests are in the areas of computational vision, motion control systems, robotics and microsystems.
Stable Algebraic Surfaces for 3D Object Representation Turker Sahin · Mustafa Unel
Published online: 29 April 2008 © Springer Science+Business Media, LLC 2008
Abstract Linear fitting techniques by implicit algebraic models usually suffer from global stability problems. Ridge regression regularization can be used to improve the stability of algebraic surface fits. In this paper a Euclidean Invariant 3D ridge regression matrix is developed and applied to a particular linear algebraic surface fitting method. Utilization of such a regularization in fitting process makes it possible to globally stabilize 3D object fits with surfaces of any degree. Robustness to noise and moderate levels of occlusion has also been enhanced significantly. Experimental results are presented to verify the improvements in global stability and robustness of the resulting fits. Keywords Algebraic surfaces · Implicit polynomials · Fitting · Stability · Ridge regression
1 Introduction Representing objects in intensity and range images is a central problem in computer vision, and there has been a significant body of work devoted to 2D and 3D object representation. An important method among these is algebraic curve and surface representations, which model objects by fitting implicit polynomials (IPs) to their boundary data [1–6]. Algebraic curves/surfaces model objects by contributions from
T. Sahin Department of Computer Engineering, Gebze Institute of Technology, Gebze, Kocaeli, Turkey e-mail: [email protected] M. Unel () Faculty of Engineering and Natural Sciences, Sabanci University, Tuzla, Istanbul, Turkey e-mail: [email protected]
all points in their data sets, and result in global representations. There are also other techniques in the literature, such as parametric curves/surfaces, which model data by local polynomial patches [10, 11]; Fourier descriptors that represent shapes based on Fourier theory methods according to polar or arc-length parameterizations [13]; superquadrics which fit lower degree surfaces to sections of complicated objects of usually ordered symmetry [12]; level set methods representing objects by propagating fronts converging to their boundary data [14, 15]; and radial basis type neural network classifiers [16, 17]. Using radial basis functions, very complicated objects such as human faces can be represented with great accuracy, far more than with polynomials, whose global nature makes it almost impossible to use for fitting shapes with fine details. So why use polynomial curves/surfaces at all? We think that polynomial models have important advantages which make them preferable in many computer vision tasks. Many natural and man-made objects, which are too complicated to be represented by a simple ellipse/ellipsoid, can well be approximated by a higher degree curve/surface. These models are robust versus noise and reasonable levels of missing data. They are also very suitable for point classification problems and curve-surface intersection calculations, as this involves only a simple polynomial evaluation. New methods which enable more optimal data set representations have been emerging [5–7], making the results of these calculations more accurate than before. Algebraic curve and surface models are preferred over Fourier Descriptors as they are able to represent non-star, occluded shapes as well as unordered data; and over point based modelling approaches that suffer from noise, missing data and have significant dependence on accurate point correspondences [24]. Moreover, the geometric invariants of 2D algebraic models have been derived [8, 9], which have been
128
applied to important vision problems such as efficient pose estimation [22] and motion estimation and tracking of objects [23]. One drawback of implicit algebraic models has been considered as their containing high degree terms, making them prone to instability. Some nonlinear techniques [2, 5, 6, 19] have been attempted for IP modelling, but like any other nonlinear optimization they need a good initial guess and they were not computationally efficient. Moreover, reference [25] explains why these nonlinear methods don’t provide satisfactory fits in all cases and the need for further improvement. Thus, linear approaches such as [20, 21, 25] are preferred. The work in [25] brings substantial improvement over those previous approaches not only by providing higher quality fits, but also significantly reducing the fitting time. However, it is limited to 2D curves only. In this paper we present 3D ridge regression (RR) regularization as an effective method for overcoming global instability problem of a particular fitting technique, the 3L algorithm [21]. We first present an alternative level set generation scheme for 3L algorithm, which is not properly addressed in [21], to improve local stability and then derive a Euclidean invariant RR matrix to globally stabilize the resulting surfaces. RR for improving the global stability of algebraic curves first appeared in [25]; however, as mentioned above that solution was limited to 2D curve fitting problems. This work extends RR regularization to 3D models, which helps in fitting stable and accurate algebraic surfaces of virtually any degree to complex data. Most of the objects which are considered in our work are not simple surfaces. Some of them have complicated structures like holes etc. in their data sets. Proposed method takes these structures into account rather than avoiding them. A first draft of this paper appears in [26]. The structure of this text is as follows: in Sect. 2 we give a brief introduction to the implicit surface models followed by the general implicit polynomial fitting problem and introduction to linear fitting techniques in Sects. 3 and 4. In Sect. 5 we explain the 3L fitting method, and introduce a new level set generation technique we have employed for 3D data. Section 6 is on the 3D ridge regression method in detail. In Sect. 7 we present experimental results and discussions followed by conclusions in Sect. 8.
2 Implicit Surface Models 2.1 Algebraic Surfaces IP fitting techniques model 3D object data by nth degree algebraic surfaces with c = (n + 1)(n + 2)(n + 3)/6 coeffi-
J Math Imaging Vis (2008) 32: 127–137
cients of the form:
Fn (x, y, z) =
aij k x i y j zk
0≤i,j,k;i+j +k≤n
= a000 + a100 x + a010 y + a001 z + · · · H0
H1 (x,y,z)
+ an00 x + · · · + a00n zn n
Hn (x,y,z)
T = 1 x . . . zn a000 a100 . . . a00n mT
= mT a
a
(2.1)
where m is the vector of monomials, a is the coefficient vector and each ternary form Hr (x, y, z) is a homogeneous polynomial of degree r in x, y and z. 2.2 Data Set Normalization Data set normalization is an important technique, which first centers the mean point of data at origin and then scales the data to be as near the unit circle as possible. By this way the difference between the values of higher and lower degree monomial terms are reduced to reasonable levels, such that the resulting IP’s are more stable. We have applied the radial distance normalization, which is a linear process that improves the results obtained using other data normalizations such as [25]. Accordingly the center of mass of original data, (xi , yi , zi )N i=1 , is shifted to origin by: (xˆi , yˆi , zˆ i ) = (xi , yi , zi ) − μ
(2.2)
where μ = N1 N i=1 (xi , yi , zi ) is the mean point of original points. Then the normalized data is calculated by scaling this result according to: (xˆni , yˆni , zˆ ni ) =
(xˆi , yˆi , zˆ i ) s
2 2 2 1/2 is the average radial where s = N1 N i=1 (xˆ i + yˆi + zˆ i ) distance of the centralized data.
3 General Implicit Algebraic Model Fitting Problem The problem of representing a data set or some of its sections by an algebraic curve or surface is referred to as the polynomial fitting problem. Accordingly a cost function like the following is minimized J=
N 1 Dist2 (Pi , Z(Fn )) N i=1
(3.3)
J Math Imaging Vis (2008) 32: 127–137
where Dist(Pi , Z(F )) is the distance of the data point Pi = (xi , yi , zi )N i=1 to the zero set Z(Fn ) = {(x, y, z) : Fn (x, y, z) = 0} of the IP model Fn (x, y, z). There are two main categories of IP fitting techniques: the nonlinear methods, which use the true geometric distance or Euclidean approximation, and linear methods based on the algebraic distance. The nonlinear approaches are invariant to transformations in the Euclidean space and are not biased [3, 4]. However, except for very basic shapes, there is no available closed form expression for the Euclidean distance. Therefore, time consuming iterative optimization procedures must be carried out. For this reason linear techniques are preferred, which formulate (3.3) as a linear least squares problem resulting in faster computation times [20, 21, 25].
4 Linear Algebraic Surface Fitting Techniques The simplest linear surface fitting approach is the Least Squares method, which minimizes the algebraic distance over the data sets (Pi )N i=1 with an additional constraint for avoiding the trivial solution. The solution of this optimization problem is the unit eigenvector a associated with the smallest eigenvalue of Sa = λa, where S = M T M is the scatter matrix of monomials. Here M is the monomials matrix composed of the monomial vectors mi of data points such that:
T M = m1 m2 . . . mN (4.4) This technique is affine invariant [20], but unfortunately is very unstable, which makes it unsuitable for any practical application. Thus more satisfactory linear fitting methods are preferred over the least squares technique. These methods make use of data set normalization and local continuity measures, both of which improve accuracy and stability. One example is Gradient One algorithm, which applies gradient control by incorporating scatter matrices from local normals and tangents in addition to the original scatter matrix of monomials [25]. In this method, the normal and tangential components are imposed to be 1 and 0 for the resulting surface gradient to be as near to unity as possible. However this technique is more sensitive to changes in data normalization and in its tuning parameter than 3L method. Also the practical Euclidean invariance properties of Gradient One are weaker, which make 3L a more suitable candidate for vision tasks such as motion estimation and object tracking.
5 The 3L Fitting Method The 3L technique incorporates control for local continuity by placing two additional data layers at a distance ±ε out-
129
side and inside the original data as the optimization constraint [21]. Thus the resulting fit is constrained to be no more than ±ε away from the actual points. The solution of this problem is: a = M †b
(5.5)
where M † = (M T M)−1 M T is the pseudo-inverse of the matrix of monomials, M. Similar to the other linear approaches, this method achieves local stability about the data, but unfortunately cannot provide a solution for the global instability problem. However, it is less sensitive to changes in its tuning parameter, ε, and in applied data normalization. For this reason we have applied ridge regression regularization to this method, which gives better overall results. A crucial aspect of this method is level set generation. To achieve better fits, we have triangulated the 3D data sets, and then calculated the point normals as the average of the neighboring triangular faces for forming the inner and outer level sets. The related triangulation and level set formation methods are summarized in the next sections. 5.1 Triangulation The crust triangulation algorithm of [18] has been applied for forming the triangular tessellations of object data. This method uses Voronoi diagram and Delaunay triangulation approaches to find the topological connection of data points. The algorithm first computes the Voronoi diagram of the sample set. For this reason it forms syntectic cubic volumes called Voronoi cells for every data point. It then selects the two furthest vertices of the cells from each other, which are referred to as poles. One pole has to be inside the data surface while the other must be outside. Next the method computes the Delaunay triangulation of the union set of the original sample points and their poles. Finally, only the triangles whose all three vertices are sample points are selected as the triangle patches. The algorithm returns the list of triangular faces and the neighboring vertices for data points. This method may not be able to detect some triangles in the cases of very scarce or noisy data. However, it gives reasonable results with sufficiently dense and evenly spaced points. The object data we have employed and corresponding triangular meshes are in Figs. 1 and 2 respectively to depict the results of the employed algorithm. 5.2 Level Set Generation The triangular faces obtained from the crust algorithm are used to calculate the local normals of data points for obtaining the virtual inner and outer data layers correctly. For sufficiently accurate operation of the 3L technique, the neighboring triangular patch normals of each data point should be
130
J Math Imaging Vis (2008) 32: 127–137
Fig. 1 The 3D data used in our fits. Examples include a torus, a phone receiver, an eight shape, a sphere, a duck, a heart, an apple, a cube, a face, a twisted object, a golf stick handle, a head, a ball joint, a hip and a patella bone
Fig. 2 Triangulation of the sample data by Crust Algorithm
correctly aligned at each point to get a good overall representation of the level sets. Hence after the neighboring faces of a point have been obtained from the generated mesh, the directions of their normals are rearranged to satisfy the right hand rule for local alignment as in Fig. 3(a). Then these normals are summed and normalized to give the unit normal at the data point. There may be problems, which prevent accurate local alignment, when object data has very complex topology or is noisy. Such effects cause the neighboring triangle arrays (resulting from the algorithm) not to be of purely planar topology. Thus to avoid possible errors, these exceptional ill-structured patch arrays are checked beforehand and eliminated by the application of a single normal, which is perpendicular to the plane fit of the point and the neighboring vertices. The resulting point normals may be globally unaligned, that is they may not be all pointing to
outside the data or vice versa. Hence they should also be globally checked for having correct directions. To achieve this we have applied the local continuity of gradient principle to the neighborhood of each data point. Accordingly, the direction differences of adjacent point normals should have angles less than 90◦ . Starting with an arbitrary point, its normal can be used as a reference to align each neighboring normal by a simple dot product check. To be more robust versus the discrepancies in the local continuity of the data set, a margin can be set aside so that the normal difference angles can be kept at a significantly lower degree than acute angle for the checking process. This procedure is repeated for the adjacent point as reference, until all normals have been readjusted accurately. This global alignment is important for placing the synthetic layers correctly. In cases of correct global alignment, all points of the inner and outer
J Math Imaging Vis (2008) 32: 127–137
(a) The point normals Ni obtained from face normals NFj by right hand rule
(b) Cross-section of the level sets of a 3D free form object Fig. 3 Local alignment of surface point normals and the resulting level sets
level sets should stay ±ε inside and outside the actual data as in Fig. 3(b) of the 2D silhouette of the 3 layers of a 3D free form object. Another important factor for satisfactory 3L fits is the choice of the ε parameter. This selection is important for correct alignment of the 2 synthetic level sets, as depicted in Fig. 3(b), which prevents the least squares problem of (5.5) from being ill-posed. ε should be selected not too small for flexibility and not too high for the correct alignment of inner and outer layers at the detailed sections of the data sets. For the radial distance normalization that we have applied, choices of ε, which give better results are typically between 0.01–0.05. It is not recommended to exceed these values, as otherwise the detailed sections may not be preserved accurately.
6 Ridge Regression Regularization for Stable Algebraic Surface Fits Linear fitting methods in general have difficulty in achieving global stability. Especially in modelling complicated shapes, many extra open surface pieces are generated, which make the resulting fits nearly impossible to be interpreted. An important reason for global instability is the collinearity of modelled data, which cause the M T M matrix of the associated least squares problem to be nearly singular with some
131
eigenvalues negligible compared to the others. Such eigenvalues contribute very little to the overall shape of the fit, but are one of the main causes of instability. The ratio of the largest eigenvalue to the smallest one is referred to as the condition number and is a measure of the singularity of matrices. The work in [25] employs RR regularization as a way to counteract on such singularity effects. Fitting results are impressive, but limited to 2D. Our work will handle ridge regression regularization of algebraic surfaces in 3D. RR improves the condition numbers and the stability of implicit polynomials by adding a new κD term to M T M data scatter matrices. Accordingly, (5.5) is modified as: aκ = (M T M + κD)−1 M T b
(6.6)
where κ is the RR parameter and D is a positive definite, symmetric matrix, which will be selected as diagonal to facilitate easier tunability. For RR to achieve its purpose, κ should be increased from 0 to higher values until a stable closed bounded surface is obtained. By this way the small eigenvalues are significantly increased, while the larger eigenvalues that contribute more to the overall shape of the fit are virtually unaffected, hence reducing the condition number to reasonable levels. Thus provided that the κ term is not selected too high, unstable fits can be stabilized with no loss of shape. 6.1 The Euclidean Invariance of Ridge Regression Regularization The other important part of the RR term is the D matrix of regularization, which is selected such that Euclidean invariance is not violated while at the same time the shape representation is preserved. Accordingly, when an arbitrary Euclidean transformation E such as: ⎡ ⎤ ⎡ ⎤ ⎤ x ⎡ x R t ⎢ ⎥ ⎢y ⎥ ⎢ ⎥=⎣ ⎦ ⎢y ⎥ (6.7) ⎣z ⎦ ⎣z⎦ 0 0 0 1 1 1 E
is applied to data points where R is the 3-D rotational part and t = [tx ty tz ]T is the translational vector, the original and transformed surface monomial vectors m(x, y, z) and m (x , y , z ) will be related by a c × c square Eˆ matrix of block diagonal rotational matrices part and corresponding translational sections: ˆ m = Em
(6.8)
As the initial and transformed surfaces should be equivalent under Euclidean transformations, we can relate the two coefficient vectors by mT a = mT a
(6.9)
132
J Math Imaging Vis (2008) 32: 127–137
As Eˆ is composed of block diagonal rotational matrix sections, which all have inverses on their own, it is an invertible matrix. Thus applying (6.8) to (6.9) we get the following relation between the original and transformed coefficient vectors.
QT 1/2 DT 1/2 QT = T 1/2 DT 1/2
a = (Eˆ T )−1 a
D = T −1
(6.10)
For Euclidean invariance, we should be able to derive the relation (6.6) for the transformed coefficient vector, a , as well. Thus using (6.10) and (6.6): a = (Eˆ T )−1 (M T M + κD)−1 M T b
(6.11)
Here substituting Eˆ −1 Eˆ = I after the inverted monomial matrix, the following form is obtained: ˆ T M)Eˆ T + κ ED ˆ Eˆ T −1 EM ˆ Tb a = E(M
(6.12)
By (4.4) and (6.8), this equation can be rearranged to: ˆ Eˆ T )−1 M T b a = (M T M + κ ED
(6.13)
which should preserve Euclidean invariance when the RR part is equal to κD. This is satisfied for: ˆ Eˆ T = D ED
(6.15)
where Q is an arbitrary block diagonal rotation matrix with blocks corresponding to the ternary forms of (2.1). T is the diagonal matrix of trinomial coefficients, whose elements are: (i + j + k)! Tll = i!j !k!
(6.16)
where i, j, k are incremented as follows: (i, j, k) = (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), . . . , (n, 0, 0), . . . , (0, n, 0), . . . , (0, 0, n) and the index l of the diagonal terms are incremented from 1 to the number of IP coefficients, according to the above variation of i, j, k exponents by the formula: (j + k)(j + k + 1) 2 (i + j + k)(i + j + k + 1)(i + j + k + 2) +1 + 6
l=k+
As Q is a rotation matrix, i.e. QQT = I , a Euclidean invariant D matrix can be derived as: (6.17)
which is the inverse of the trinomial coefficients matrix. This selection of D does not maintain the local shape representation, but this can be fixed by multiplying it by the weighted sum of the diagonal elements of the monomials matrix M T M. Such terms correspond to the square of the distance of the normalized data points to the origin, thus preserve Euclidean invariance and improve the local shape properties at the same time. The terms of the resulting diagonal D matrix are: Dll = βi+j +k
1 Tll
(6.18)
Here the term βi+j +k is contributed from the data set points (Pt )N t=1 formulated as: βi+j +k =
N (xt2 + yt2 + zt2 )i+j +k t=1
(6.14)
If (6.7) is reduced to a pure rotation, which is the case for normalized data, Eˆ can be represented by the similarity transformation: Eˆ = T −1/2 QT 1/2
Substituting (6.15) into (6.14) implies that
As a simple example the D matrix that should work on quadric surface monomial matrices should be: 1 1 1 D = diag β0 , β1 , β1 , β1 , β2 , β2 , β2 , β2 , β2 , β2 2 2 2 2 2 2 where β0 = N , β1 = N t=1 (xt + yt + zt ) and β2 = N 2 2 2 2 t=1 (xt + yt + zt ) . 6.2 The Closed-Boundedness of Ridge Regression in 3D The RR regularized 3L method (3LRR) shows the important boundedness property of RR fits as the κ term is increased towards infinity in value. Accordingly provided that the modelled object data sets are closed (i.e. with no significant gaps or occlusion) and the employed IPs are of even degree, the resulting fits always converge to a bounded surface. This behavior can be proved as follows: From (6.6), the limit of the 3LRR IP vector as κ → ∞ is: aij k =
1 −1 T D M b κ
(6.19)
Here the D matrix is according to the selection in (6.18). Thus we concentrate on the M T b part, composed of the sums of monomials in ascending order, such that: T i j k i j k M b= ε xy z − xy z (6.20) +ε
−ε
c×1
J Math Imaging Vis (2008) 32: 127–137
with i, j, k in the order of (6.16) and ±ε represent the sums of monomials in the outer and inner data layers, respectively. This summation is in the form of directional derivative at each data point in the outward direction from the data set surface, i.e. (x i y j zk ) (x i y j zk )+ε − (x i y j zk )−ε = s 2ε
MT b = S
T N · ∇(x i y j zk ) c×1 dS
diag{1/Dll }
=
(i + j + k)! βi+j +k i!j !k!
(6.22)
(i + j + k)! βi+j +k i!j !k!
(6.26)
dxdydz
i>1,j >1,k>1; i+j +k=n
n(n − 1) = βn
i>1,j >1,k>1; i+j +k=n
xˆ 2
(n − 2)! (i − 2)!j !k!
ˆ j (zˆz)n−i−j × (x x) ˆ i−2 (y y) (n − 2)! (x x) ˆ i (y y) ˆ j −2 (zˆz)n−i−j i!(j − 2)!k!
(n − 2)! (x x) ˆ i (y y) + zˆ ˆ j (zˆz)n−i−j −2 i!j !(k − 2)! 2
dxdydz (6.27)
The three integrand terms are in the form of the trinomial expansion formula:
V
i(i − 1)x i−2 y j zk
(x + y + z)n =
V
aij k xˆ i yˆ j zˆ k = 0
+ zˆ 2
+ yˆ 2
n i=0
(6.24)
with V the region bounded by the surface S. From (2.1) the zero set of the IP vector a∞ is obtained to be: T a∞ m=
(i + j + k)! j (j − 1)(x x) ˆ i (y y) ˆ j −2 (zˆz)k βi+j +k i!j !k!
As i + j + k = n, (6.26) can be re-expressed as: aij k xˆ i yˆ j zˆ k
(6.23)
∇ 2 (x i y j zk )dxdydz
+ j (j − 1)x i y j −2 zk + k(k − 1)x i y j zk−2 dxdydz
+ yˆ 2
ˆ j (zˆz)k−2 ×k(k − 1)(x x) ˆ i (y y)
(i + j + k)! βi+j +k i!j !k!
ˆ j (zˆz)k ×i(i − 1)(x x) ˆ i−2 (y y)
MT b
(i + j + k)! βi+j +k i!j !k!
xˆ 2
i>1,j >1,k>1; i+j +k=n
scaled by a constant term. Applying the divergence theorem and utilizing ∇ · ∇f = ∇ 2 f where ∇ 2 represents the Laplacian operator, we get: aij k =
=
(6.21)
scaled with a constant, where N T is the local normal at each data point and ∇(x i y j zk ) is the gradient of the monomials. These terms are integrated over the complete surface of the data, S. Hence (6.19) can be brought into the form:
T (i + j + k)! aij k = N · ∇(x i y j zk ) c×1 dS βi+j +k i!j !k! S
aij k xˆ i yˆ j zˆ k
i>1,j >1,k>1; i+j +k=n
with s being the vector in the direction of the derivative. As the directional derivative of the monomial terms is related to the gradient in the outward direction (i.e. the normal) to the surface, we also express the M T b terms as:
133
(6.25)
0≤i+j +k≤n
where mT = (xˆ i yˆ j zˆ k ) is the monomial vector. If we can show that the leading form of (6.25) is strictly positive, we complete the proof of the boundedness of 3LRR fits as the RR term tends to infinity. By embedding (6.24) to (6.25) and concentrating on the leading form, we obtain:
n! x i y j z(n−i−j ) i!j !(n − i − j )!
(6.28)
therefore the leading form can be simplified to: aij k xˆ i yˆ j zˆ k i>1,j >1,k>1; i+j +k=n
=
n(n − 1)(xˆ 2 + yˆ 2 + zˆ 2 ) βn × (x xˆ + y yˆ + zˆz)n−2 dxdydz
(6.29)
Here if n is even, the integrand and all the other terms become positive. Hence the leading form of the 3LRR IP equation in 3D is strictly positive, imposing the resulting fits to converge to an ellipsoid like closed bounded surface. The proof that we have just presented for the closedboundedness of the resulting fits ensures the existence
134
J Math Imaging Vis (2008) 32: 127–137
Fig. 4 Least squares fits of 3D sample data: the torus, the phone receiver, eight, cube and the heart
Fig. 5 Stable 3L surface fits of the torus, the phone, the cube, the heart and the sphere objects
of some κ values which can stabilize the fits. In fact, as shown in the experimental section, a κ value in the range 10−4 –10−3 stabilizes the resulting fits in virtually all cases by forcing the unbounded/disconnected components to move away from the data and finally being eliminated from the fits. 6.3 Suitable Selection of κ Term for Shape Representation and Recognition The RR term κ reduces the variability of parameters and makes the IP model less sensitive to small perturbations in data. It also has the effect of smoothing the object surface as κ is increased in value. However, for an as efficient as possible shape modelling, the IP model should not be smoothed to much, and yet all extra open components should be eliminated. Thus starting from zero, the κ parameter can be iteratively increased to higher values, until the surface is closed bounded. The boundedness of the surface can be checked according to the IP leading form positive-definiteness principle introduced in Sect. 6.2. Accordingly, the smallest value of κ obeying this test should be the most suitable selection for shape modelling. For object recognition purposes, the main objective should be to minimize the error between the actual 3D data and the fitted surface polynomial. However, there is no single optimal κ value that obeys this criterion and κ for best recognition may vary between different sets of data; however, this value should be higher than the threshold value suitable for shape modelling. 7 Experimental Results and Discussion The conventional linear fitting approaches in 3D is rather limited to modelling objects of simpler topologies. Among these methods the least squares algorithm cannot provide stable and accurate fits for any data except the torus as depicted in Fig. 4. The stability of the torus is mainly the result of the perfect symmetry of this object, and is still limited to a 4th degree fit. When slightly more complicated shapes such
Fig. 6 Unstable 3L surface fits, which include an eight shape, a duck, the man face, the woman face, the golf, ball joint, hip, the twisted object and the patella in this order
as the phone and the eight shaped object are employed, accuracy of fits is lost significantly as in the same figure. Finally when more complicated structures like the cube and heart are modelled, this method fails completely. The 3L fitting is adequate for simple data shapes that can be modelled by lower degree fits. More data examples than least squares fit like the phone, the cube, the heart, and the sphere can be stably and accurately fit as depicted in Fig. 5. However, for practical applications, more complex objects have to be fit with stable IP surfaces. The remainder of our sample objects are in this category, for which the 3L method fails to give good IP models. The related fits are depicted in Fig. 6, all with extra surface pieces.
J Math Imaging Vis (2008) 32: 127–137
135
Fig. 7 The stabilization of the Eight, Face and twisted objects by RR regularization
Fig. 8 The stabilized surface fits of objects in Fig. 6 by RR
RR works to achieve global stability as the κ parameter is increased from zero to sufficiently higher values. Accordingly, it gives stable fits for the objects that the 3L method could not manage to model adequately. Some examples are presented in Figs. 7 and 8. In Fig. 7, the stabilization of surface fits of the eight shaped object, the face and the twisted object are depicted. Here the eight and the face objects are modelled by 6th degree surfaces, while the twisted object has an 8th degree fit. The leftmost sub-figures are the 3L only subplots with κ = 0. The following intermediate subplots are for higher values of κ to depict how RR moves the extra surfaces away from the zero set and eliminates
them. The stabilizations are achieved for the face and the twisted object at κ = 10−5 and for the eight shaped object at κ = 10−4 . These are typical values for the RR parameter and the fits have very good shape representation, although lower degree IP’s are employed. As in zoomed rightmost sub-figures, the eight and twisted objects are modelled with all their major details, and the face surface contains critical sections like the nose and the forehead. All the objects, 3L method could not stabilize in Fig. 6 can be modelled accurately by RR regularization as presented in Fig. 8, which verify the associated global stabilization properties. Moreover RR is robust versus perturbations and structural problems in data sets. As in Fig. 9, despite a significant section of the objects such as 15% of the golf, 12% of the duck, and 10% of the face and the eight, have been discarded, RR gives stable and accurate surface fits for all data. The only difference is that the RR parameter was increased slightly more to achieve stabilization; however, this virtually unaffects the shape of the fits. Robustness versus noise property of this method is also strong. The noise tests are presented in Fig. 10 for the eight, patella, twisted object and the face. The objects of smoother data sets with less fine detail are robust versus noise levels of 5% of object sizes, as with σ 2 = 0.0025 for the eight and σ 2 = 0.002 for the patella. It should be noted that the employed data have been normalized to unit ball and the 3L parameter has been set to ε = 0.05. Thus the inner and outer synthetic data levels are at ±5% of the object size away from the actual data surface, and at higher noise levels, they are indistinguishably independent of RR. As a result we can conclude that RR can restore the stability of least squares type fitting meth-
136
J Math Imaging Vis (2008) 32: 127–137
Fig. 9 The face, eight, golf and the duck RR fits under missing data conditions
Fig. 10 RR fits for the eight, patella, twisted and face data subjected to noise levels of variances σ 2 = 0.0025, σ 2 = 0.002, σ 2 = 0.0009 and σ 2 = 0.0004, respectively
ods while preserving accuracy up to the levels permitted by these techniques. The twisted and the face objects with more fine detail can also be fit stably and accurately up to noise levels of 0.0009 and σ 2 = 0.0004, respectively. At higher noise levels, stability is still ensured although some object fine detail cannot be preserved satisfactorily.
of the modelled objects as much as possible. Thus globally stable and accurate fits of higher degrees have been obtained for range data with fine details. The resulting fits are also robust versus noise levels of up to 5% of object size and 15% missing data. As a result this fitting technique should be suitable for modelling 3D objects with algebraic surfaces of any practical degree, and thus should find important applications in 3D modelling and computer vision.
8 Conclusions In this paper a Euclidean invariant 3D ridge regression matrix is developed and applied for the regularization of 3L fitting technique. This approach results in closed bounded algebraic surfaces, while preserving the local shape properties
References 1. Pilu, M., Fitzgibbon, A., Fisher, R.: Ellipse specific direct least squares fitting. In: Proc. IEEE International Conference on Image Processing, Lausanne, Switzerland, September 1996
J Math Imaging Vis (2008) 32: 127–137 2. Keren, D., Cooper, D., Subrahmonia, J.: Describing complicated objects by implicit polynomials. IEEE Trans. Pattern Anal. Mach. Intell. 16, 38–53 (1994) 3. Faber, P., Fischer, R.B.: Euclidean fitting revisisted. In: 4th IWVF, pp. 165–175 (2001) 4. Faber, P., Fischer, R.B.: Pros and cons of Euclidean fitting. In: Proc. Annual German Symposium for Pattern Recognition, DAGM01, Munich. LNCS, vol. 2191, pp. 414–420. Springer, Berlin (2001) 5. Keren, D., Gotsman, C.: Fitting curves and surfaces with constrained implicit polynomials. IEEE Trans. Pattern Anal. Mach. Intell. 21(1), 31–41 (1999) 6. Keren, D.: Topologically faithful fitting of simple closed curves. IEEE Trans. Pattern Anal. Mach. Intell. 26(1), 118–123 (2004) 7. Redding, N.J.: Implicit polynomials, orthogonal distance regression, and the closest point on a curve. IEEE Trans. Pattern Anal. Mach. Intell. 22(2), 191–199 (2000) 8. Unel, M., Wolovich, W.A.: On the construction of complete sets of geometric invariants for algebraic curves. Adv. Appl. Math. 24, 65–87 (2000) 9. Unel, M., Wolovich, W.A.: A new representation for quartic curves and complete sets of geometric invariants. Int. J. Pattern Recognit. Artif. Intell. 13(8), 1137–1149 (1999) 10. Huang, Z., Cohen, F.S.: Affine-invariant B-spline moments for curve matching. IEEE Trans. Image Process. 5(10), 1473–1480 (1996) 11. Cohen, F.S., Wang, J.-Y.: Modeling image curve using invariant 3-D object curve models—a path to 3-D recognition and shape estimation from image contours. IEEE Trans. Pattern Anal. Mach. Intell. 16(1), 1–12 (1994) 12. Terzopoulos, D., Metaxas, D.: Dynamic 3D models with local and global deformations: deformable superquadrics. IEEE Trans. Pattern Anal. Mach. Intell. 13(7), 703–714 (1991) 13. Wu, M.F., Shen, H.T.: Representation of 3D surfaces by twovariable Forurier descriptors. IEEE Trans. Pattern Anal. Mach. Intell. 20(8), 583–588 (0000) 14. Osher, S., Sethian, J.A.: Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. J. Comput. Phys. 79, 12–49 (1988) 15. Malladi, R., Sethian, J.A., Vemuri, B.C.: Shape modelling with front propagation: a level set approach. IEEE Trans. Pattern Anal. Mach. Intell. 17(2), 158–175 (1995) 16. Dinh, H.Q., Turk, G., Slabaugh, G.: Reconstructing surfaces by volumetric regularization using radial basis functions. IEEE Trans. Pattern Anal. Mach. Intell. 24(10), 1358–1371 (2002) 17. Zhang, M., Ma, L., Zeng, X., Wang, Y.H.: Image based 3D face modelling. In: Proceedings of the IEEE International Conference on Computer Graphics, Imaging and Visualization, CGIV04 (2004) 18. Amenta, N., Choi, S., Kolluri, R.K.: The power crust. In: Proceeding of the 6th ACM Symposium on Solid Modelling and Applications, pp. 249–266 (2001) 19. Taubin, G., et al.: Parametrized families of polynomials for bounded algebraic curve and surface fitting. IEEE Trans. Pattern Anal. Mach. Intell. 16(3), 287–303 (1994) 20. Deepak, K., Mundy, J.: Fitting affine invariant conics to curves. In: Geometric Invariance in Machine Vision. MIT Press, Cambridge (1992) 21. Blane, M., Lei, Z., et al.: The 3L algorithm for fitting implicit polynomial curves and surfaces to data. IEEE Trans. Pattern Anal. Mach. Intell. 22(3), 298–313 (2000)
137 22. Unel, M., Wolovich, W.A.: Pose estimation and object identification using complex algebraic representations. Pattern Anal. Appl. J. 1(3), 178–188 (1998) 23. Sahin, T., Unel, M.: An adaptive estimation method for rigid motion parameters of 2D curves. In: LNCS, vol. 3212, pp. 355–362. Springer, Heidelberg (2004) 24. Ramesh, J., Kasturi, R., Schunc, B.G.: Machine Vision. McGraw Hill, New York (1995) 25. Tasdizen, T., Tarel, J.-P., Cooper, D.B.: Improving the stability of algebraic curves for applications. IEEE Trans. Image Process. 9(3), 405–416 (2000) 26. Sahin, T., Unel, M.: Fitting globally stabilized algebraic surfaces to range data. In: Proeceedings of the 10th IEEE International Conference on Computer Vision, Beijing, China, October 2005
Turker Sahin has received his B.Eng. Degree in Electronic Engineering with Optoelectronics from University College London in 1999 and his M.Sc. Degree in Control Systems from Imperial College in 2001. During 2003–2007, he was a Ph.D. candiate and research assistant in the Department of Computer Engineering at Gebze Institute of Technology, Kocaeli, Turkey where he did work on algebraic and differential geometry of curves and surfaces, visual motion estimation and mobile robotics. He obtained his Ph.D. degree in January 2008. His current research interests are in the areas of vision, control and robotics. Mustafa Unel received both the BSEE degree in 1994 and the M.S. degree in Systems and Control Engineering in 1996 from Bogazici University, Istanbul, Turkey. He received a second M.S. degree in Applied Mathematics in 1998 and the Ph.D. degree in Electrical Sciences and Computer Engineering in 1999 both from Brown University, Providence, USA. During 1999–2001, he was a postdoctoral fellow at Center for Computational Vision and Control, Yale University, New Haven, USA. In Spring 2001, he joined the Department of Computer Engineering at Gebze Institute of Technology, Turkey as an assistant professor. Since Fall 2004, he has been with the Faculty of Engineering and Natural Sciences, Sabanci University, Istanbul, Turkey where he is currently an associate professor. His research interests are in the areas of computational vision, motion control systems, robotics and microsystems.
