Document Sample

620 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 Orthogonal Distance Fitting of Implicit Curves and Surfaces È Sung Joon Ahn, Wolfgang Rauh, Hyung Suck Cho, Member, IEEE, and Hans-Jurgen Warnecke AbstractÐDimensional model fitting finds its applications in various fields of science and engineering and is a relevant subject in computer/machine vision and coordinate metrology. In this paper, we present two new fitting algorithms, distance-based and coordinate-based algorithm, for implicit surfaces and plane curves, which minimize the square sum of the orthogonal error distances between the model feature and the given data points. Each of the two algorithms has its own advantages and is to be purposefully applied to a specific fitting task, considering the implementation and memory space cost, and possibilities of observation weighting. By the new algorithms, the model feature parameters are grouped and simultaneously estimated in terms of form, position, and rotation parameters. The form parameters determine the shape of the model feature and the position/rotation parameters describe the rigid body motion of the model feature. The proposed algorithms are applicable to any kind of implicit surface and plane curve. In this paper, we also describe algorithm implementation and show various examples of orthogonal distance fit. Index TermsÐImplicit curve, implicit surface, curve fitting, surface fitting, orthogonal distance fitting, geometric distance, orthogonal contacting, nonlinear least squares, parameter estimation, Gauss-Newton method, parameter constraint, parametric model recovery, object segmentation, object classification, object reconstruction. æ 1 INTRODUCTION W ITH image processing, pattern recognition, and compu- ter/machine vision, dimensional model (curve and surface) fitting to a set of given data points is a very common shortest distance point on an implicit model feature from a given point. Section 3 describes two new algorithms, the distance-based and coordinate-based algorithms, for ortho- task carried out during a working project, e.g., edge detection, gonal distance fitting of implicit model features. Various information extraction from 2D-image or 3D-range image, fitting examples are given in Section 4 and Section 5, and object reconstruction. For the purpose of dimensional including examples of real data fitting. This paper ends with a summary in Section 6. model fitting, we can consider three methods, namely, moment method [15], [32], [35], [42], Hough transform [8], 1.1 Least-Squares Dimensional Model Fitting [19], [28], and least-squares method (LSM) [23]. The moment The least-squares method [23] is one of the best known, and method and Hough transform are efficient for fitting of most often applied, mathematical tools in various disci- relatively simple models, while their application to a complex plines of science and engineering. With applications of LSM object model or to an object model with a large number of to dimensional model fitting, the natural and best choice of model parameters is not encouraged. In this paper, we the error distance is the shortest distance between the given consider the LS-fitting algorithms for implicit model features. point and the model feature. This error definition is By data modeling and analysis in various disciplines of prescribed in coordinate metrology guidelines [16], [30]. A science and engineering, implicit features are very often used century ago, Pearson [32] elegantly solved the orthogonal because of their compact description in form of f a; x H distance fitting problem of line and plane in space by using and because of the possibility of a simple on-off and inside- the moment method (an art of LSM). However, except for some simple model features, computing and minimizing outside decision. the square sum of the shortest error distances are not simple We review the current LS-fitting algorithms for implicit tasks for general features. For applications, an algorithm for model features and introduce new approaches in this dimensional model fitting would be very advantageous if section. Section 2 contains two algorithms for finding the estimation parameters are grouped and simultaneously estimated in terms of form, position, and rotation para- . S.J. Ahn and W. Rauh are with the Department of Information Processing, meters [24]. Furthermore, it would be helpful if the Fraunhofer Institute for Manufacturing Engineering and Automation reliability (variance) of each parameter could be tested. (IPA), Nobelstr. 12, 70569 Stuttgart, Germany. There are two main categories of LS-fitting problems of E-mail: {sja, wor}@ipa.fhg.de. . H.S. Cho is with the Department of Mechanical Engineering, Korea implicit features, algebraic and geometric fitting, and these Advanced Institute of Science and Technology (KAIST), Kusung-dong are differentiated by their respective definition of the error 373-1, Yusung-gu, 305-701 Taejeon, South Korea. distances involved [21], [36]. E-mail: hscho@lca.kaist.ac.kr. Algebraic fitting is a procedure whereby model feature is . H.-J. Warnecke is with the Fraunhofer Society, Leonrodstr. 54, 80636 described by implicit equation F b; X H with algebraic Munich, Germany. E-mail: warnecke@zv.fhg.de. parameters b and the error distances are defined with the Manuscript received 3 Apr. 2001; revised 25 Oct. 2001; accepted 29 Jan. 2002. Recommended for acceptance by S. Sarkar. deviations of functional values from the expected value For information on obtaining reprints of this article, please send e-mail to: (i.e., zero) at each given point. If F b; Xi T H, the given tpami@computer.org, and reference IEEECS Log Number 113920. point Xi does not lie on the model feature (i.e., there is 0162-8828/02/$17.00 ß 2002 IEEE AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 621 some error-of-fit). Most publications about LS-fitting of 1.2 New Approaches implicit features have been concerned with the square sum We will now introduce two versatile and very efficient of algebraic distances or their modifications [12], [36], [40] orthogonal distance fitting algorithms for implicit surfaces Â m h m iP and plane curves. The general goal of orthogonal distance ÃP P H F b; Xi or P H F b; Xi =krF b; Xi k : fitting of a model feature to m given points is the estimation of i i the feature parameters a minimizing the performance index I P d P Pd H Q In spite of advantages in implementation and computing cost, the algebraic fitting has drawbacks in accuracy. The or estimated parameters are not invariant to coordinate P X À XH P P X À XH ; R transformation (e.g., a simple parallel shift of the given H points set causes changes not only in position but also in where form and rotation of the estimated model feature) and, generally, we cannot find a physical interpretation of an . X: Coordinate column vector of the m given points, algebraic error definition. Finally, resolving the algebraic X X ; F F F ; X , X Xi ; Yi ; Zi ; I m i parameters b into the physical parameters (e.g., in terms of . XH : Coordinate column vector of the m correspond- form, position, and rotation parameters) is not a simple task ing points on the model feature; for general features [24]. . d: Distance column vector, d dI ; F F F ; dm , In geometric fitting, also known as best fitting or orthogonal q distance fitting, the error distance is defined as the shortest di kXi À XHi k Xi À XHi Xi À XHi Y distance (geometric distance) of the given point to the model feature. However, contrary to its clear and preferable definition, the geometric distance has nonlinear nature with . P P: Weighting matrix or error covariance matrix the model parameters and, as said above, computing and (positive definitive); minimizing the square sum of geometric distances are not . P: Nonsingular symmetric matrix [17]. easy tasks for general features. Thus, including the normal- If we assume independent identically distributed (i.i.d.) ized algebraic distance F b; Xi =krF b; Xi k shown in the measurement errors, the two performance indexes (3) and (4) second formula of (1) [36], [40], various approximating have the same value up to the a priori standard deviation of measures of the geometric distance are used for dimensional the measurement errors. model fitting [14], [29], [31]. Because using approximating We are calling the new algorithms minimizing the measures of the geometric distance cannot satisfactorily fill performance indexes (3) and (4) distance-based algorithm the gap between algebraic and geometric fitting, the devel- and coordinate-based algorithm, respectively. Our algorithms opment of geometric fitting algorithms has been pursued in are the generalized extensions of an orthogonal distance the literature. While there are geometric fitting algorithms for fitting algorithm for implicit plane curves [3], [4]. Because explicit [11], and parametric features [5], [13], [21], [38], [41], no assumption is made about the mathematical handling of we are considering in this paper only fitting algorithms for implicit features, the proposed algorithms can be applied to implicit features. Sullivan et al. [39] presented a geometric any kind of implicit dimensional feature. Each of the two fitting algorithm for implicit algebraic features algorithms consists of two nested iterations: À Á q min min P fXHi agm : H iI S F b; X bj Xkj Y lj Z mj H; aH m fXi giI PF j The inner iteration of (5), to be carried out at each outer minimizing the square sum of the geometric distances iteration cycle, finds the shortest distance points on the d b; Xi kXi À XHi k, where Xi is a given point and XHi is current model feature from each given point (Section 2) and the shortest distance point on the model feature F from Xi . In the outer iteration cycle updates the model feature para- order to locate XHi , they have iteratively minimized the meters (Section 3). Once appropriate break conditions are Lagrangian function L ; X kXi À XkP F b; X. To satisfied, e.g., when no more significant improvement on obtain the partial derivatives @d b; Xi =@b (these are neces- the performance indexes (3) and (4) could be gained sary for the nonlinear iteration updating the algebraic through parameter updating, the outer iteration terminates. parameters b), they have utilized the orthogonal contacting By the new algorithms, the estimation parameters a are equations defined in machine coordinate system XYZ grouped in three categories and simultaneously estimated. First, the form parameters ag (e.g., radius r for a circle/ F cylinder/sphere, three axis lengths a, b, c for an ellipsoid) F b; Xi ; X 0: P rF Â Xi À X describe the shape of the standard model feature f defined A weakness of Sullivan et al. algorithm, the same as in the case in model coordinate system xyz (Fig. 1) of algebraic fitting, is that the physical parameters are f ag ; x H with ag aI ; F F F ; al : T combined into an algebraic parameters vector b. And, as we will show in Section 3.1, they have unnecessarily expensively The form parameters are invariant to the rigid body motion obtained the partial derivatives @d b; Xi =@b from (2). of the model feature. The second and the third parameters Helfrich and Zwick described a similar algorithm in [27]. group, respectively, the position parameters ap and the 622 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 which usually has only a few form parameters. Functional interpretations and treatments of the position/rotation parameters concerning the coordinate transformation (7) are identical for all implicit model features. 2 ORTHOGONAL CONTACTING POINT The time consuming part of the proposed orthogonal distance fitting algorithms in this paper is the inner iteration of (5), finding the shortest distance points (orthogonal contacting points) fXHi gm on a general model feature from iI each given point fXi gm . For a given point xi in frame xyz iI xi R Xi À Xo ; W we determine the orthogonal contacting point xHi on the standard model feature (6). Then, the orthogonal contacting point XHi in frame XYZ from the given point Xi will be obtained through a backward transformation of xHi into XYZ. For some implicit features (e.g., circle, sphere, cylinder, cone, Fig. 1. Implicit features and the orthogonal contacting point xHi in frame and torus), the orthogonal contacting point xHi can be xyz from the given point Xi in frame XYZ: (a) plane curve and (b) surface. determined in closed form at a low computing cost. However, in general, we must iteratively find xHi satisfying rotation parameters ar , describe the rigid body motion of the appropriate orthogonal contacting conditions. We describe model feature f in machine coordinate system XYZ two algorithms for finding the orthogonal contacting point xHi in frame xyz. The first algorithm is derived from the general X RÀI x Xo or x R X À Xo ; where U properties of the shortest distance point. The second algorithm is based on a constrained minimization method, R R R' R! rI rP rQ ; RÀI R ; the method of Lagrange multipliers [20]. Within the follow- V ing sections, the orthogonal contacting point is marked with a ap Xo Xo ; Yo ; Zo ; nd ar !; '; : prime, as XHi in frame XYZ and as xHi in frame xyz, once it has We characterize the form parameters as intrinsic parameters been determined and is ready to be used. Otherwise, it is and the position/rotation parameters as extrinsic parameters of notated with the general coordinate vector X or x. the model feature according to their context. In this paper, we intend to simultaneously estimate all these parameters 2.1 Direct Method A necessary condition for the shortest distance point x on a a ; a ; a g p r the implicit model feature (6) from the given point xi is that aI ; F F F ; al ; Xo ; Yo ; Zo ; !; '; aI ; F F F ; aq : the connecting line of xi with x should be parallel to the feature normal rf at x (Fig. 1): For plane curve fitting, we simply ignore all terms concerning z, Z, !, and '. We may try to describe the model feature F with rf Â xi À x 0 ; where r @=@x; @=@y; @=@z : algebraic parameters b in machine coordinate system XYZ IH F b ag ; ap ; ar ; X f ag ; x ap ; ar ; X H; Equation (10) is equivalent to although this has no interest for us from the viewpoint of rx f= xi À x ry f= yi À y rz f= zi À z our model fitting algorithms. and only two of three equation rows of (10) are indepen- In comparison with other dimensional model fitting dent. Equation (10) describes a space curve as the algorithms, our algorithms simultaneously estimate the intersection of two implicit surfaces defined in frame xyz physical parameters a of the model feature in terms of by two independent equation rows of (10). We can find that form ag , position ap , and rotation parameters ar , thus (10) generally satisfies the parallel condition of the vector providing very useful algorithmic features for various xi À x to the feature normal rf of the iso-feature of (6) applications, e.g., with robot vision, motion analysis, and coordinate metrology. For such applications, algebraic para- f ag ; x À onst H: II meters b must be converted into physical parameters a even though the conversion is not an easy task for general model As the constant in (11) is continuously varying, the features. trajectories of the points lying on (11) and satisfying (10) The novelty of our algorithms, as will be described in draw a curve described by (10). We name the curve (10) Section 3, is that the necessary derivative values for nonlinear centrortho-curve about the point xi to the iso-features (11). iteration are to be obtained from the standard model feature Our aim is to find the intersection of the centrortho-curve equation (6) and from the coordinate transformation equa- (10) with the implicit model feature (6). In other words, we tions (7). For application of our algorithms to a new model would like to find the point x which simultaneously feature, we need merely to provide the first and the second satisfies (6) and (10) (orthogonal contacting condition, see derivatives of the standard model feature equation (6) Fig. 2 for the case of a plane curve, an ellipse): AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 623 Fig. 3. Iterative search (direct method) for the orthogonal contacting point xHi on f ag ; x H from the given point xi . The points xi;I and xi;P are, respectively, the first and the second approximation of xHi . Fig. 2. Isocurves xP =TR yP =IT À onst H (onst > H) and the trajec- the isocurve tangent in xy-plane [3], [4]. From this fact and for tories (centrortho-curves) Qx À V Qy P IT H (thick line) and Qx À Vy H (broken line) of the orthogonal contacting points on the the sake of a common program module for curves and isocurves from the point (2, 2) and (2, 0), respectively. surfaces, we interchange the second and the fourth row of (12). Then, for the case of plane curves, we activate only the f first two rows of the modified equation (12) in the program f ag ; xi ; x 0 with xi R Xi À Xo : rf Â xi À x module (see Appendix A and Fig. 2 for an ellipse). IP 2.2 Constrained Minimization Orthogonal contacting equation (12) governs the variational Alternatively, we can find the orthogonal contacting point behavior of the orthogonal contacting point in frame xyz minimizing the Lagrangian function [20] relative to the differential change of the feature parameters a and will play an important role with the coordinate-based L ; x xi À x xi À x f: IR algorithm to be described in Section 3.2. We directly solve (12) for x by using a generalized The first order necessary condition for a minimum of the Newton method starting from the initial point of xH xi Lagrangian function (14) is (Fig. 3) (how to compute the matrix @f =@x will be shown in rL ÀP xi À x rf Section 3.2): @L 0: IS @ f @f We iteratively solve (15) for x and [27], [39] Áx Àf xk ; xkI xk Áx: IQ @x k PI H rf Áx If necessary, especially in order to prevent the point update rf H Á Áx from heading for divergence zone with an implicit IT P xi À x À rf @ feature having very high local curvature, we underweight ; where H rf: Àf @x the on-the-feature condition in (12) as below: In practice, we have observed that, once the second term of wf the Lagrangian function (14) becomes negative (the first term f ag ; xi ; x 0; with H < w I: rf Â xi À x is always nonnegative), it is occasionally prone to become In the first iteration cycle with the starting point of xH xi , more negative (i.e., diverging) in the succeeding iteration the linear system (13) and its solution (i.e., the first moving cycles of (16) to minimize (14), although the second term is step) are equivalent to commanded to be zero. Consequently, the convergence of Àf iteration (16) is very affected by the selection of the initial rf Á Áx Àf nd Áx rf : values. From this observation and with a hint at determining rf Â Áx xi 0 xi krfkP x i the first moving step size taken from the direct method It should be noted that the first moving step is the negative described in Section 2.1, we use the following initial values for of the normalized algebraic error vector shown in the x and : second formula of (1) which has been regarded in the V literature [36], [40] as a good approximation of the b f=krfkP ` X xH xi xi geometric error vector. H I X xH T xi nd f ag ; xH ! H b X For plane curves on the xy-plane (i.e., z H), we simply ÀI X xH T xi nd f ag ; xH < H: ignore the second and the third row of (12), which have terms IU of coordinate z in the cross product and have been automatically satisfied, then the fourth row of (12) without Iteration (16) with the initial values in (17) converges terms of coordinate z in the cross product is nothing other relatively fast. Nevertheless, if the convergence of iteration than the orthogonality condition of the error vector xi À x to (16) fails, we catch it using the direct method. 624 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 2.3 Verification of the Orthogonal Contacting Point If no closed form solution of the orthogonal contacting point xHi is available, thus, if we have to find xHi through iteration, we would like to start the inner iteration described in Section 2.1 and Section 2.2, preferably from the initial point of xH xi because we need not provide a special initial point (blind search). And, as another advantage, the initial point of xH xi almost always leads the iteration (13) and (16) to the global minimum distance point that is just the nearest point on the model feature from xi . However, despite this advantageous initial point selection, the global minimum distance of xHi from xi could not be generally guaranteed. Conditions (12) and (15) merely constrain the connecting line of point x on f ag ; x H with the given point xi to be parallel to the feature normal rf at x. Thus, the orthogonal contacting point xHi found by iterations (13) or (16) is, in principle, only a local extreme (maximum or minimum) distance point. For example, there is a maximum of four local extreme distance points on an ellipse for a given point (Fig. 2). While there is no general sufficient Fig. 4. Verification of the orthogonal contacting point. The point xHi;I condition for global minimum distance of xHi from xi , we satisfies (18), while xHi;P does not: (a) f ag ; xi > H. (b) f ag ; xi < H. check the correctness of the point xHi by testing multiple independent necessary conditions. of the orthogonal distance fitting if we use the current With many a well-known model feature having axis or orthogonal contacting point XHi as a reasonable initial point plane symmetry (e.g., ellipse, parabola, hyperbola, ellipsoid, for the inner iteration inside the next outer iteration cycle paraboloid, etc.), the two points xHi and xi must lie in the same half-plane/space of the model coordinate frame xyz [4]. From xH;aÁa RaÁa XHi;a À Xo;aÁa : IW this fact, we can simply check if xHi is not the global minimum distance point from xi by comparing the coordinate signs However, a pitfall of using (19) must be mentioned. If the between xHi and xi . If necessary, we selectively invert the current orthogonal contacting point XHi once becomes a local coordinate signs of xHi and repeat iteration (13) or (16) starting (not global) minimum distance point after an update of the from the modified xHi . parameters a, it can be inherited into the next outer iteration as We can also check the correctness of point xHi by long as we use (19) (Fig. 5). In order to interrupt this examining the direction of the feature normal rf at xHi . unintended development, we periodically (e.g., with every Once iteration (13) or (16) is successfully terminated, we 5-10th outer iteration cycle) start the inner iteration from the verify the direction of the feature normal by using another given point Xi during the second half-phase of the outer orthogonal contacting condition (Fig. 4): iteration rf H Á xi À xi À sign f ag ; xi H: H kx À xH k IV xH;aÁa RaÁa Xi À Xo;aÁa : PH krfk x i i i We summarize the strategy for an efficient and inexpen- Rarely, if (18) could not yet be satisfied with a value sive finding of the orthogonal contacting point xHi : (probably ÆP) other than zero, then the connecting line of the two points xi and xHi pierces the model feature an odd . Use the closed form solution of xHi , if available; number of times (i.e., at least, xHi is not the global minimum . Start the inner iteration from the initial point xH of (20) distance point). In this case, we repeat iteration (13) or (16) in the first half-phase of the outer iteration, during which the parameter update Áa is generally wild; with an adjusted starting point (e.g., xi;I on the connecting . In the second half-phase of the outer iteration, start line of the two points xi and xHi in Fig. 4). the inner iteration from the initial point xH of (19); 2.4 Acceleration of Finding the Orthogonal . Periodically, in the second half-phase of the outer Contacting Point iteration, start the inner iteration from the initial The use of the initial point of xH xi for the iterative search of point xH of (20). the orthogonal contacting point xHi is very convenient and preferable, however, it demands a relatively high computing 3 ORTHOGONAL DISTANCE FITTING cost. Fortunately, we are aware of the fact that the parameter update Áa will generally have been moderate in the second With the algorithms described in Section 2 (inner iteration), half-phase of the outer iteration (after the first 5-10 outer the shortest distance points fXHi gm on the current model iI iteration cycles with most fitting tasks) and the orthogonal feature from each given point fXi gm are found and now iI contacting point XHi in frame XYZ changes its position very available for the parameter updating (outer iteration) in the incrementally between two consecutive outer iteration cycles. nested iteration scheme (5). The inner iteration and the Then, we can dramatically save the computing cost of the parameter updating are to be alternately repeated until inner iteration and, consequently, the overall computing cost appropriate break conditions are satisfied. In this section, we AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 625 From (7) and (25), we get À Á @di Xi À XHi @X Jdi ;a À @a kXi À XHi k @a XXH i À Á H Xi À Xi @x @RÀI @Xo À RÀI x kXi À XHi k @a @a @a xxH i À Á H xi À xi @x À kxi À xHi k @axxH À Á 2 i 3 Xi À Xi H @RÀI À 0 I xHi ; Fig. 5. Orthogonal contacting points XHi;I and XHi;P on the updated feature kXi À XHi k @ar F a Áa; X H from the given point Xi . Iterative search of the orthogonal contacting point may converge to XHi;I , if the iteration started PT from the current orthogonal contacting point XHi on F a; X H (19), while starting from the given point Xi converges to XHi;P (20), being where kXi À XHi;I k > kXi À XHi;P k. a a ; a ; a ; g p r present two algorithms for updating the model feature H I x 0 parameters, the distance-based algorithm and the coordi- @R @R @R @R f FF g nate-based algorithm minimizing the performance index (3) ; nd x d F e: @ar @! @' @ and (4), respectively. 0 x 3.1 Distance-Based Algorithm Furthermore, from the facts xi À xHi == rfjxxH and i The first order necessary condition for a minimum of the @x @f @x @f performance index (3) as a function of the parameters a is rf À @a @x @a @a @f @ P @d À 0 0 y f a; x f ag ; x H; H PJ P Pd 0 ; where J : PI @ag @a @a the Jacobian matrix Jdi ;a in (26), eventually becomes We may iteratively solve (21) for a by using the Newton H I method sign xi À xHi rf @f f g Jdi ;a d 0 0 e À Á J P PJ HP Pd Áa ÀJ P Pd; @Pd where H P ; krfk @ag @a xxHi À Á PP Xi À XHi @RÀI H À 0 I xi : or, more conveniently, using the Gauss-Newton method, kXi À XHi k @ar ignoring the second derivatives term in (22) PU In comparison with Sullivan et al.'s algorithm [39], our J P PJk Áa ÀJ P Pdk ; akI ak Áa: PQ algorithm, including (24) and (27), simultaneously esti- In this paper, we iteratively estimate the parameters a by mates the parameters a in terms of form, position, and using the direct form of the Gauss normal equation (23) rotation parameters. And, Sullivan et al.'s algorithm obtains PJjk Áa ÀPdjk ; akI ak Áa; PR its Jacobian matrix unnecessarily expensively from (2). If we apply a similar procedure of (26) and (27) to F b; X with break conditions for iteration (24) H with algebraic parameters b defined in machine coordi- kJ P Pdk % H or kÁak % H or P k ÀP kI % H: H H nate system XYZ, the Jacobian matrix for the Sullivan's algorithm can be obtained as below at a lower implementa- The first break condition denotes that the error distance tion cost without using (2): vector and its gradient vectors in the linear system (24) should be orthogonal and is equivalent to condition (21). If this is sign Xi À XHi rF @F satisfied, there will be practically no more parameter update Jdi ;b ; with (kÁak % H) and improvement on the performance index of P . krF k @b H H XXi In order to complete the linear system (24), we have to r @=@X; @=@Y ; @=@Z : provide the Jacobian matrix of each distance di between the given point Xi and the orthogonal contacting point XHi on the current model feature 3.2 Coordinate-Based Algorithm In a similar way to the procedure of (21)-(24) carried out in q di kXi À XHi k Xi À XHi Xi À XHi : PS Section 3.1, we construct the Gauss-Newton iteration minimizing the performance index of (4) 626 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 @XH generality, the second and the fourth row of (12). Then, for J ; PJjk Áa P X À XH jk ; akI ak Áa; PV @a plane curve fitting, we consider only the first two rows of the with break conditions modified equation (12), in other words, we consider only the upper-left block of the FHG matrix in (31) (see Appendix A). kJ P P X À XH k % H or kÁak % H or P k ÀP kI % H H H 3.3 Comparison of the Two Algorithms and with the Jacobian matrices of each orthogonal contacting Each of the distance-based and coordinate-based algorithms point XHi on the current model feature, directly derived for orthogonal distance fitting of implicit features is versatile from (7), and very efficient from the viewpoint of application to a new model feature. We would like to stress that only the standard @X @x @RÀI @Xo model feature equation (6), without involvement of the JXHi ;a RÀI x @a XXH @a @a @a xxH position/rotation parameters, is required in (31) for imple- i i PW mentationofeitherofthetwoalgorithmstoanewmodelfeature ÀI @x @RÀI H (the second derivatives @rf=@ag is required only for coordi- R 0 I xi : @axxH @ar nate-based algorithm). The overall structure of our algorithms i remains unchanged for all dimensional fitting problems of The derivative matrix @x=@a at x xHi in (29) describes the implicit features. All that is necessary for orthogonal distance variational behavior of the orthogonal contacting point xHi in fitting of a new implicit model feature is to derive the frame xyz relative to the differential changes of the FHG matrix of (31) from (6) of the new model feature and to parameters vector a. Purposefully, we obtain @x=@a from supply a proper set of initial parameter values aH for iterations the orthogonal contacting equation (12). Because (12) has an (24) and (28). This fact makes possible a realization of versatile implicit form, its derivatives lead to and very efficient orthogonal distance fitting algorithms for implicit surfaces and plane curves. An overall schematic @f @x @f @xi @f @f @x @f @xi @f information flow is shown in Fig. 6 (coordinate-based 0 or À ; algorithm). We compare the two algorithms from the view- @x @a @xi @a @a @x @a @xi @a @a point of implementation and computing cost as follows: QH . Common features where @xi =@a is, from (9), E For implementation to a new model feature, @xi @R @Xo Xi À Xo À R only the standard model feature equation (6) is @a @a @a eventually required; @R 0 ÀR Xi À Xo : E Computing time and memory space cost are @ar proportional to the number of the data points; The other three matrices @f =@x, @f =@xi , and @f =@a in (13) and E The two algorithms demand about the same (30) are to be directly derived from (12). The elements of these computing cost because the time-consuming part of the overall procedure is the finding of three matrices are composed of simple linear combinations of the orthogonal contacting points. components of the error vector xi À x with elements of the . Different features following three vector/matrices rf, H, and G (FHG matrix): E The distance-based algorithm demands smaller @f @f @f @ @ f memory space (with plane curve fitting one half rf ; ; ; H rf; G ; QI @x @y @z @x @ag rf and with surface fitting one third of what the H I coordinate-based algorithm does); H H H E To a new model feature, the implementation of @f f yi À y f À xi À x H g g the distance-based algorithm is relatively easy f gH @x d À zi À z H xi À x e because it does not require the second deriva- H zi À z À yi À y tives @rf=@ag of the model feature; H @f @f @f I H I E With the coordinate-based algorithm, each co- H H H @x f @f @y @z g ordinate Xi ; Yi ; Zi of a measurement point can f f À @f @f H g f @y À @f @x H g g @f f @y @x g be individually weighted, while only pointwise f @f g; f @f g; weighting is possible with the distance- fÀ H @f g @xi f @z d H À @f g @x e d @z @x e based algorithm. This algorithmic feature of H @f À @f H À @f @z @f @y the coordinate-based algorithm is practical @z @y H I for measuring equipment whose measuring I H H H accuracy is not identical between measuring @f f H f yi À y À xi À x H g g axes. From this fact, we can say the coordinate- f g G 0 0 : @a d H À zi À z H xi À x e based algorithm is the more general one than H H zi À z À yi À y the distance-based algorithm; E A comparison between (26) and (29) reveals Now, (30) can be solved for @x=@a at x xHi and, conse- quently, the Jacobian matrix (29) and the linear system (28) di N Xi À XHi ; i Jdi ;a ÀN JXHi ;a ; i can be completed and solved for parameter update Áa. Xi À XHi For the sake of a common program module for surfaces where Ni : kXi À XHi k and plane curves, we have interchanged, without loss of AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 627 Fig. 6. Information flow for orthogonal distance fitting of implicit features (coordinate-based algorithm). 3.4 Parameter Constraint PJ UWV If any parameter constraint (e.g., on size, area, volume, W J Â Ã position, orientation of the model feature) with U U V V I; W dig wI ; F F F ; wq ; fj a À onstj H; j I; F F F ; n; QP and now the linear systems (34) and (35) can be solved for or, in vector form, Áa. After a successful termination of iterations (34) and (35), along with the performance index P , the Jacobian matrix H f a À const 0 provides useful information about the quality of parameter is to be additionally applied, we append any following estimations as follows: n equations (33) with the large weighting values w to the linear systems (24) and (28): . Covariance matrix: 4 5ÀI @fj À Á PJ PJ wj Áa Àwj fj a À onstj ; j I; F F F ; n: QQ gov ^ a VWÀP V Y @a W J W J Considering parameter constraints (32), we can rewrite the performance indexes (3) and (4) and the linear systems (24) . Parameter covariance: and (28) as follows: Vji Vki q Pd Pd a gov j ; ak ; j; k I; F F F ; q Y QT P H ; iI wPi W f a À const W f a À const PJ Pd QR Áa À ; . Variance of parameters: W J W f a À const @d @f P J ; J P j a H a gov j ; aj ; j I; F F F ; q Y @a @a mnÀq for the distance-based algorithm and . Correlation coefficients: P X À XH P X À XH P H ; a gov j ; ak W f a À const W f a À const j ; ak p ; a j; k I; F F F ; q: a gov j ; aj gov k ; ak a PJ P X À XH QS Áa ; W J ÀW f a À const Using the above information, we can test the reliability of @XH @f J ; J the estimated parameters ^ and the propriety of the model a @a @a selection (object classification). For example, we may try to for the coordinate-based algorithm. fit an ellipsoid to a set of measurement points of a sphere- 3.5 Parameter Test and Object Classification like object surface. Then, besides a % b % c and the existence The Jacobian matrix in (34) and (35) will be decomposed by of strong correlations between the rotation parameters and SVD [25], [33], the other parameters, we get relatively large variances of 628 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 Fig. 7. Orthogonal distance fit to the points set in Table 1: (a) Circle fit. (b) Ellipse fit. (c) Superellipse fit. (d) Convergence of the fit. Iteration number 0-6: circle, 7-21: ellipse, and 22-: superellipse fit. the rotation parameters. In this case, the rotation parameters processing and pattern recognition is very desirable because are redundant here (overparameterized) and, although the superellipse can represent various 2D-features, e.g., rectangle fitting of an ellipsoid to the points set has better (" ( I), ellipse (" I), and star (" > P). Many researchers performance than a sphere fitting according to the have tried to fit a superellipse by using the moment method performance index P , the proper (reliable) model feature H [10], [43] or using the LSM with algebraic distances [34]. for the points set is a sphere. In order to demonstrate the outstanding performance of The parameter covariance (36) with the distance-based our geometric fitting algorithms, we fit an extreme shape of a algorithm is generally a little larger than that with the superellipse, a rectangle. We obtain the initial parameter coordinate-based algorithm. If P is small enough, justifying H values with " I from a circle fitting and an ellipse fitting, ignoring of the second derivatives term in (24) and (28), successively. In detail, we use the gravitational center and rms there is practically no difference of the parameter covar- central distance of the given points set as the initial parameter iance (36) between the two algorithms. values for circle fitting, circle parameters for ellipse fitting [3], [4], and, finally, ellipse parameters for superellipse fitting (Fig. 7, Table 2). Superellipse fitting to the eight points in 4 FITTING EXAMPLES Table 1 representing a rectangle terminated after seven 4.1 Superellipse Fitting iteration cycles for kÁak U:V Á IHÀW (Fig. 7d, Table 2). A superellipse [22] is a generalization of an ellipse and it is described by TABLE 1 f a; b; "; x jxj=aP=" jyj=bP=" ÀI H; Eight Coordinate Pairs Representing a Rectangle where a, b are the axis lengths and exponent " is the shape coefficient. The use of superellipses in applications of image AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 629 TABLE 2 Results of the Orthogonal Distance Fitting to the Points Set in Table 1 4.2 Cone Fitting the constraint weighting value w is large and, thus, locks the The standard model feature of a cone in frame xyz can be parameters from changing. If we perturb this quasi-equili- described below: brium state, similarly to the random walking technique, through adding a small artifact error to the initial parameter f ; x xP yP À z tn =PP H; values, then we could expect a faster convergence. With a where , the only form parameter, is the vertex angle of a nonzero initial vertex angle value of =IH, the iteration cone. The position Xo , the origin of frame xyz, is defined at terminates after only 13 cycles for kÁak P:W Á IHÀT (Fig. 8d). the vertex. The orientation of a cone is represented by the direction cosines of the z-axis, rQ in (8). However, from the 4.3 Torus Fitting viewpoint of coordinate metrology [30], a better parameter- A torus is a quartic surface, and is described by ization of a cone is f rI ; rP ; x xR yR zR P xP yP yP zP zP xP P P P f ; r; x x y À r À z tn =P H; À P rP rP xP yP P rP À rP zP rP À rP P H; P I P I P I with a constraint on the position and rotation parameters where rI is the radius of the circular section of the tube and (see (32)) rP is the mean radius of the ring. We obtain the initial parameter values from a 3D-circle fitting [5] (Fig. 9) with f ap ; ar Xo À X rQ !; ' H; q where the second form parameter r is the sectional radius of rI P H;irle =m; rP rirle : a cone cut by the xy-plane (z H) and X is the gravitational center of the given points set. Torus fitting to the 10 points in Table 5 representing a half As one of the test data sets used in an authorized torus is terminated after nine iteration cycles for kÁak I:R Á IHÀU (Fig. 9d, Table 6). certification [18], [30] process of our algorithms, the 10 points in Table 3 representing only a quarter of a cone slice were 4.4 Superquadric Fitting prepared by the German federal authority PTB (Physikalisch- A superquadric (superellipsoid) [9] is a generalization of an Technische Bundesanstalt). We have obtained the initial ellipsoid and is described by parameter values with H from a 3D-circle fitting [5] and a cylinder fitting, successively (Fig. 8, Table 4). The cone fitting f a; b; c; "I ; "P ; x "I ="P to the points set in Table 3 with the constraint weighting jxj=aP="I jyj=bP="I jzj=cP="P ÀI H; value w I is terminated after 60 iteration cycles for kÁak P:Q Á IHÀU . The convergence is somehow slow because where a, b, c are the axis lengths and exponents "I , "P are the the initial cone parameters (i.e., the cylinder parameters) seem shape coefficients. In comparison with the algebraic fitting to already provide a local minimum and it is especially slow if algorithm [37], our algorithms can also fit extreme shapes of TABLE 3 Ten Coordinate Triples Representing a Quarter of a Cone Slice 630 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 Fig. 8. Orthogonal distance fit to the points set in Table 3: (a) 3D-circle fit. (b) Circular cylinder fit. (c) Circular cone fit. (d) Convergence of the fit. Iteration number 0-4: 3D-circle, 5-16: cylinder, and 17-:cone fit with the initial value of =IH. a superquadric (e.g., a box with "I;P ( I or a star with Superquadric fitting to the 30 points in Table 7 representing "I;P > P). After a series of experiments with numerous sets a box is terminated after 18 iteration cycles for kÁak of data points, we have concluded that our algorithms U:W Á IHÀU (Fig. 10d). The intermediate ellipsoid fitting safely converge with the parameter zone (a ! b ! c) showed a slow convergence in its second half- phase (iteration number 10 to 42 in Fig. 10d) because of a H:HHP "I nd "P IH nd "I ="P SHH: relatively large variance of the estimated major axis length a of Otherwise, there is a danger of data overflow destroying the the ellipsoid caused by the distribution of the given points. In Jacobian matrices of linear systems (13), (16), (24), and (28). other words, the performance index P is not being changed H Additionally, a sharp increase of the computing cost with very much by the variation of a (compare the standard inner iteration is unavoidable. We obtain the initial parameter values with "I "P I from a sphere fitting deviation of the intermediate ellipsoid ellip: IQ:RQWQ a a and an ellipsoid fitting, successively (Fig. 10, Table 8). a with that of the superquadric H:IHQR). TABLE 4 Results of the Coordinate-Based Orthogonal Distance Fitting to the Points Set in Table 3 AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 631 Fig. 9. Orthogonal distance fit to the points set in Table 5: (a) 3D-circle fit. (b) Initial torus. (c) Torus fit. (d) Convergence of the fit. Iteration number 0-6: 3D-circle, and 7-: torus fit. 5 REAL DATA FITTING because they distort the fitting results. Conceptually, our 5.1 Data Acquisition algorithms possess very advantageous algorithmic features In order to show how our algorithms perform on real data, a for segmentation and outlier detection of dimensional sphere and an ellipsoid will be fitted to the same set of measurement data, namely, the geometric error definition measurement points of a sphere-like object surface (Fig. 11a). and the parameters grouping in terms of form, position, and The points set is acquired by our stripe projecting rotation. In this paper, utilizing these algorithmic features, we 3D-measurement system with the Coded-Light-Approach used the following interactive procedure between segmenta- (CLA) [7], [44] using the Gray code [26]. The measuring system consists of a standard CCD camera and an LCD stripe tion, outlier detection, and fitting of the measurement points: projector of ABW [1] and is calibrated by photogrammetric 1. From the following procedure, exclude the measure- bundle adjustment using the circular coded targets for ment points which belong to the known objects, e.g., automation of the calibration procedure [2], [6]. Our measuring stage or already detected and identified 3D-measurement software determines the coordinates of objects. the object surface points along the stripe edges. The number 2. Initialize the model feature parameters through a of the measured surface points is about 29,400 (Fig. 11b), model fitting to a well-conditioned object portion or from which every eighth point (subtotal 3,675 points) is used through a fitting of a simple model feature, e.g. a for the sphere and the ellipsoid fitting in this section (Fig. 12). sphere instead of the target ellipsoid. Also, initialize 5.2 Segmentation, Outlier Detection, and the safety margin t of the point searching domain Model Fitting volume to be used in the next step, e.g., tH r=P for a sphere. With real data fitting, we have to find an effective search 3. For the current model feature, determine the (segmentation) method for the measurement points which orthogonal error distances di of each measurement potentially belong to the target model feature. And, outliers point lying inside the domain volume (see the next should be detected and excluded from the model fitting section for a detailed description) and evaluate the rms error distance TABLE 5 s Ten Coordinate Triples Representing a Half Torus I P m error d ; m iI i where m is the number of the measurement points lying inside the domain volume. 632 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 TABLE 6 Results of the Coordinate-Based Orthogonal Distance Fitting to the Points Set in Table 5 Fig. 10. Orthogonal distance fit to the points set in Table 7: (a) Sphere fit. (b) Ellipsoid fit. (c) Superquadric fit. (d) Convergence of the fit. Iteration number 0-42: ellipsoid and 43-: superquadric fit. 4. Determine the inliers set of measurement points accuracy of the measuring system, terminate the lying inside the domain volume and having error procedure. distances of di tthreshold Á error , where tthreshold is the 6. Update the safety margin t 2 tgrowing Á tthreshold Á outlier threshold factor (we used tthreshold P:S, i.e., fit error of the domain volume, where tgrowing is we assumed 1.24 percent outlier probability in the the domain growing factor (we used tgrowing I:H). domain volume according to the Gaussian error 7. Repeat again from the third step. distribution). Optionally, if there is no change in the 5.3 Domain Volume for Measurement Points inliers set, terminate the procedure. In the first step of the procedure described in the previous 5. Fit the model feature to the inliers set and save the section, 256 of 3,675 measurement points were detected and resulting rms error distance fit error . If fit error is removed from the points list as they belong to the measuring smaller than a given error level concerning the stage (Z % H). In order to not only exclude the outliers from the model fitting but also avoid unnecessary determination of TABLE 7 the orthogonal contacting points in the third step of the Thirty Coordinate Triples Representing a Box procedure, we applied a two-step domain volume criterion. First, we restrict the measurement points domain within a box with a safety margin t containing the interest portion of the standard model feature (6) defined in frame xyz. In particular, for an ellipsoid, the domain box can be defined below: jxj a t; jyj b t; nd jzj c t: With the domain box criterion, the potential measurement points of the target object surface can be screened from a AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 633 TABLE 8 Results of the Coordinate-Based Orthogonal Distance Fitting to the Points Set in Table 7 large set of measurement points at minimal computing cost grouping in terms of form, position, and rotation, makes (for each measurement point, we need only a few multi- possible a simple and inexpensive search for the measure- plications for coordinate transformation). ment point candidates which potentially belong to the model Next, of the measurement points lying inside the domain feature. The time-consuming determination of the orthogo- box, we actually determine the orthogonal contacting point nal contacting points, which is necessary for the outlier if the measurement point also lies between two iso-features detection, is to be carried out only for the measurement point (11) with onstinner < H and onstouter > H, respectively. In candidates lying inside the domain volume. particular, for an ellipsoid, we use 5.4 Fitting Results onstinner f ag ; xinner with With the ellipsoid fitting, we have applied the parameter V b mx a=P; a À t; H; H X a min a; b; c constraints of ! H and ' H with a weighting value of w ` xinner H; mx b=P; b À t; H I:H Á IHV (see Section 3.4) since the object surface seems to be X b min a; b; c b X the top face of a horizontally lying ellipsoid. The computing H; H; mx c=P; c À t X c min a; b; c; time with a Pentium III 866 MHz PC was a total of 7.8 s for and 17 rounds of sphere fitting with a final parameter update size of kÁak Q:U Á IHÀU and a total of 94.1 s for 26 rounds of onstouter f ag ; xouter with ellipsoid fitting with kÁak T:S Á IHÀU . The number of the V b a t; H; H ` X a min a; b; c outliers (including the out-of-domain points) was 727 for xouter H; b t; H X b min a; b; c error H:HSWT mm and 1,005 for error H:HQHW mm with the b X H; H; c t X c min a; b; c: sphere and the ellipsoid fitting, respectively (Fig. 12d). If we compare Fig. 12b and Fig. 12c, an ellipsoid seems to be the The above two-step criterion, defining the domain better model feature than a sphere for the measured object volume at a low computing cost by utilizing the parameters surface. However, a look at Table 9 and Table 11 reveals the Fig. 11. Sphere-like object surface: (a) Stripe projection. (b) Measured points cloud. 634 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 Fig. 12. Results of the model fitting to the points cloud (Q; RIW Q; TUS À PST points) from Fig. 11b. Distance error bars are 10 times elongated. (a) Sphere fit without and (b) with outlier elimination. (c) Ellipsoid fit with outlier elimination under the constraints of ! H and ' H. (d) Convergence of the fit with outlier elimination. Iteration number 0-51: 17 rounds of sphere fitting and 52-162: 26 rounds of ellipsoid fitting. The initial sphere parameters are taken from the gravitational center and the rms central distance of the data points. After the first round of the sphere fitting, the initial number of the outliers (including the out-of-domain points) Noutlier P; SIS is reduced to IUU and, thereafter, increased. relatively large variances and the very strong correlations of 6 SUMMARY the axis lengths of the estimated ellipsoid. From this fact, a Dimensional model fitting finds its applications in various sphere could be considered rather as the proper model fields of science and engineering and is a relevant subject in feature for the measured object surface. The maximal computer/machine vision, coordinate metrology, and com- anticorrelation between the position parameter Zo and the puter-aided geometric design. In this paper, we have sphere radius a (also the axis length c of the ellipsoid, Table 10 and Table 11) should be noted which is caused by the one- presented two new algorithms for orthogonal distance fitting sided distribution of the data points on the top face of the of implicit surfaces and plane curves which minimize the sphere/ellipsoid. On the other hand, there is no correlation square sum of the orthogonal error distances between between the first two rotation parameters (! and ') and the the model feature and the given data points. Easily other parameters of the ellipsoid, which is conditioned by the understandable optimization methods, but with a fresh idea, parameter constraints of ! H and ' H. are applied to the problem of geometric fitting, a widely TABLE 9 Results of the Coordinate-Based Orthogonal Distance Fitting to the Points Cloud of Fig. 12b (Sphere Fitting) and Fig. 12c (Ellipsoid Fitting) AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 635 TABLE 10 TABLE 11 Correlation Coefficients of the Sphere Parameters in Table 9 Correlation Coefficients of the Ellipsoid Parameters in Table 9 recognized difficult problem. The new algorithms are versatile and very efficient from the viewpoint of implementation and application to a new model feature. . Rigid body motion (7) of the model in machine Our algorithms converge very well and do not require a coordinate frame XY: necessarily good initial parameter values set, which could also be internally supplied (e.g., gravitational center and rms X C ÀS x Xo ; central distance of the given points set for sphere fitting, Y S C y Yo sphere parameters for ellipsoid fitting, and ellipsoid para- with C os ; S sin : meters for superquadric fitting, etc.). Memory space and computing time cost are proportional to the number of the . FHG matrix (31): given data points. The estimation parameters are grouped in x=aP form/position/rotation parameters and are simultaneously rf ; y=bP estimated, thus providing useful algorithmic features for H I ÀxP =aQ ÀyP =bQ various applications. With the new algorithms, the quality of I=aP H f g parameter estimation (including the propriety of model H ; G d ÀPx=aQ H e: H I=bP selection) can be simply tested by utilizing the parameter H ÀPy=bQ covariance matrix. Together with other algorithms for orthogonal distance fitting of parametric features [5], our . Orthogonal contacting equation (12): algorithms are certified by the German federal authority PTB xP =aP yP =bP À I=P [18], [30], with a certification grade that the parameter f a; b; xi ; x 0: yi À yx=aP À xi À xy=bP estimation accuracy is higher than 0.1 m for length unit and 0.1 rd for angle unit for all parameters of all tested model . Linear system equation (13) of the direct method features with all test data sets. (Section 2.1): 2 APPENDIX A H H I=aP H IMPLEMENTATION EXAMPLE WITH ELLIPSE FITTING yi À y À xi À x H I=bP 3 In this paper, the geometric fitting algorithms for implicit x=aP y=bP Áx curves and surfaces are quite generally described and, in y=bP Àx=aP Áy practice, they are implemented in a highly modular manner. P P P P x =a y =b À I=P For application to a new curve/surface, the user needs merely À : yi À yx=aP À xi À xy=bP to supply the FHG matrix (31) of the standard model feature (6). Nevertheless, in order to help the interested reader realize . Linear system equation (16) of the constrained his/her own implementation, we give, in this section, an minimization (Section 2.2): example of how some key intermediate equations really look H IH I like with a particular fitting task, the ellipse fitting. P =aP H x=aP Áx f gf g d H P =bP y=bP ed Áy e . Parameters vector of an ellipse in XY-plane: x=aP y=bP H Á a a; b; Xo ; Yo ; : H I P xi À x À x=aP . Standard model (6) in model coordinate frame xy: f g d P yi À y À y=bP e: À xP =aP yP =bP À I=P I xP yP f a; b; x P À I H: P aP b 636 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 . Linear system equation (24) for the distance-based APPENDIX B algorithm (Section 3.1): GRADIENT AND HESSIAN MATRIX OF SUPERELLIPSE H I H I AND SUPERQUADRIC JdI ;a dI f F gÁa ÀPf F g; Except for superellipse and superquadric, the derivation of Pf d F g F e f F g d F e the FHG matrix (31) is straightforward with the curves and Jdm ;a dm surfaces shown in this paper. Because the power function q and the logarithm function can cause data overflow and with di Xi À Xi P Yi À YiH P : H domain error, a careful program source coding of the FHG matrix is essential for the superellipse fitting and the superquadric fitting. . Jacobian matrix Jdi ;a in (27) for the distance-based algorithm (Section 3.1): B.1 Superellipse For the sake of a compact description and efficient source À Á Jdi ;a Àgi xHi P =aQ Àgi yHi P =bQ H H H coding of the FHG matrix, we use the following super- H ellipse function: I Xi À Xi H H I H À YiH À Yo À ; 2 P=" 3 di Yi À YiH H H H I H Xi À Xo I jxj P=" jyj f ag ; x f a; b; "; x ÀI ; P a b with À Á then the FHG matrix of this superellipse function will be: sign xi À xHi xHi =aP yi À yHi yHi =bP gi q : xHi =aP P yHi =bP P I A=x P À " A=xP H rf ; H P ; " B=y " H B=yP . Linear system equation (28) for the coordinate- H I ÀA=a ÀB=b À A ln A B ln B=P based algorithm (Section 3.2): If g G d ÀPA="ax H ÀA I ln A="x e; H I H " JXI ;a H I H ÀPB="by ÀB I ln B="y H XI À XI f JYIH ;a g f YI À YIH g f g f g with A jxj=aP=" nd B jyj=bP=" : f F g f F F g Pf F gÁa Pf F g: f F g f g From the fact limx3H x ln x H, we prevent the domain error d JXH ;a e d Xm À X H e m m H of ln H by treating H ln H H as a group. JYm ;a H Ym À Ym B.2 Superquadric . Jacobian matrix JXHi ;a in (29) for the coordinate- For superquadric, we use the following function: based algorithm (Section 3.2): f ag ; x f a; b; c; "I ; "P ; x H2 I ÀI P="I P="I 3"I ="P P="P C ÀS @f @f @xi @f Id jxj jyj jzj JXHi ;a À ÀIe: S C @x @xi @a @a H P a b c xxi H H I H À YiH À Yo In order to reduce the danger of a data overflow by the H ; H H H I Xi À Xo power function, we introduce some intermediate variables cancelling similarly sized values below: with A jxj=aP="I =d; B jyj=bP="I =d; C jzj=cP="P ; @f H H I=aP H D d"I = "P ; with d jxj=aP="I jyj=bP="I : @x yi À y À xi À x H I=bP P x=a P y=b Then, the FHG matrix of the above superquadric function ; will be: P y=b Àx=aP @f H H AD "I À "P ; Hxx PA P À "I ; @xi Ày=b x=aP P "I "P xP "P @xi H H ÀC ÀS yi BD "I À "P ; H I Hyy PB P À "I ; @a H H S ÀC Àxi AD=x "I "P yP "P Id @f I H H rf BD=y e; C "P Hzz P P P À "P ; @a C=z "P z H yi À y À xi À x H I PABD "I À "P ÀxP =aQ ÀyP =bQ H H H Hxy Hyx ; f g "I "P xy "P d ÀPx=aQ H H H H e: Hyz Hzy Hzx Hxz H; H ÀPy=bQ H H H AHN ET AL.: ORTHOGONAL DISTANCE FITTING OF IMPLICIT CURVES AND SURFACES 637 H I AD=a [12] F.L. Bookstein, ªFitting Conic Sections to Scattered Data,º Computer f g Graphics and Image Processing, vol. 9, no. 1, pp. 56-71, 1979. f BD=b g @f ÀI f f g g; [13] B.P. Butler, A.B. Forbes, and P.M. Harris, ªAlgorithms for f C=c g Geometric Tolerance Assessment,º Report no. DITC 228/94, Nat'l @ag "P f g Physical Laboratory, Teddington, U.K., 1994. d D A ln A B ln B=P e [14] X. Cao, N. Shrikhande, and G. Hu, ªApproximate Orthogonal C ln C D ln D=P Distance Regression Method for Fitting Quadric Surfaces to Range Data,º Pattern Recognition Letters, vol. 15, no. 8, pp. 781-796, 1994. [15] B.B. Chaudhuri and G.P. Samanta, ªElliptic Fit of Objects in Two H I and Three Dimensions by Moment of Inertia Optimization,º P À" A "I"P P I Pattern Recognition Letters, vol. 12, no. 1, pp. 1-7, 1991. f a g f g [16] DIN 32880-1, Coordinate Metrology; Geometrical Fundamental Prin- f P À" B "I"P P g ciples, Terms and Definitions. Berlin: Beuth Verlag, German @ @f ÀAD f f b g g Standard, 1986. f H g; @ag @x "I "P x f g [17] N.R. Draper and H. Smith, Applied Regression Analysis, third ed. f "I À"P g New York: John Wiley and Sons, 1998. f A ln A B ln B "P ln A g d e È [18] R. Drieschner, B. Bittner, R. Elligsen, and F. Waldele, ªTesting "I Coordinate Measuring Machine Algorithms: Phase II,º BCR "P I ln D Report no. EUR 13417 EN, Commission of the European Communities, Luxemburg, 1991. H I [19] R.O. Duda and P.E. Hart, ªUse of the Hough Transformation to P À" A "I"P P Detect Lines and Curves in Pictures,º Comm. ACM, vol. 15, no. 1, f a g pp. 11-15, 1972. f g f P À" B "I "P P I g [20] R. Fletcher, Practical Methods of Optimization. New York: John @ @f ÀBD f f b g g Wiley & Sons, 1987. f H g; [21] W. Gander, G.H. Golub, and R. Strebel, ªLeast-Squares Fitting of @ag @y "I "P y f g f "I À"P g Circles and Ellipses,º BIT, vol. 34, no. 4, pp. 558-578, 1994. f A ln A B ln B " ln B g [22] M. Gardiner, ªThe Superellipse: A Curve that Lies between the d P e "I Ellipse and the Rectangle,º Scientific Am., vol. 213, no. 3, pp. 222- " I ln D 234, 1965. H I P H [23] C.F. Gauss, Theory of the Motion of the Heavenly Bodies Moving about f g the Sun in Conic Sections (Theoria Motus Corporum Coelestium in f H g Sectionibus Conicis Solem Ambientum), first published in 1809, @ @f ÀC f g P f P=c g: f g translation by C.H. Davis. New York: Dover, 1963. @ag @z "P z f g [24] R.N. Goldman, ªTwo Approaches to a Computer Model for d H e Quadric Surfaces,º IEEE Computer Graphics and Applications, vol. 3, no. 9, pp. 21-24, 1983. I ln C [25] G.H. Golub and C. Reinsch, ªSingular Value Decomposition and Least Squares Solutions,º Numerische Mathematik, vol. 14, no. 5, pp. 403-420, 1970. [26] F. Gray, ªPulse Code Communication,º US Patent 2,632,058, Mar. REFERENCES 17, 1953. [1] ABW GmbH, http://www.abw-3d.de/home_e.html, Feb. 2002. [27] H.-P. Helfrich and D. Zwick, ªA Trust Region Method for Implicit [2] S.J. Ahn, ªCalibration of the Stripe Projecting 3D-Measurement Orthogonal Distance Regression,º Numerical Algorithms, vol. 5, System,º Proc. 13th Korea Automatic Control Conf. (KACC '98), pp. 535-545, 1993. pp. 1857-1862, 1998 (in Korean). [28] P.V.C. Hough, ªMethod and Means for Recognizing Complex [3] S.J. Ahn and W. Rauh, ªGeometric Least Squares Fitting of Circle Patterns,º US Patent 3,069,654, Dec. 18, 1962. and Ellipse,º Int'l J. Pattern Recognition and Artificial Intelligence, [29] G. Hu and N. Shrikhande, ªEstimation of Surface Parameters vol. 13, no. 7, pp. 987-996, 1999. Using Orthogonal Distance Criterion,º Proc. Fifth Int'l Conf. Image [4] S.J. Ahn, W. Rauh, and H.-J. Warnecke, ªLeast-Squares Orthogo- Processing and Its Applications, pp. 345-349, 1995. nal Distances Fitting of Circle, Sphere, Ellipse, Hyperbola, and [30] ISO 10360-6:2001, Geometrical Product Specifications (GPS)ÐAc- Parabola,º Pattern Recognition, vol. 34, no. 12, pp. 2283-2303, 2001. ceptance and Reverification Tests for Coordinate Measuring È [5] S.J. Ahn, E. Westkamper, and W. Rauh, ªOrthogonal Distance Machines (CMM)ÐPart 6: Estimation of Errors in Computing Fitting of Parametric Curves and Surfaces,º Proc. Int'l Symp. Gaussian Associated Features, Int'l Standard, ISO, Dec. 2001. Algorithms for Approximation IV: Huddersfield 2001, I. Anderson and [31] D. Marshall, G. Lukacs, and R. Martin, ªRobust Segmentation of J. Levesley, eds., 2002. Primitives from Range Data in the Presence of Geometric [6] S.J. Ahn, W. Rauh, and S.I. Kim, ªCircular Coded Target for Degeneracy,º IEEE Trans. Pattern Analysis and Machine Intelligence, Automation of Optical 3D-Measurement and Camera Calibra- vol. 23, no. 3, pp. 304-314, Mar. 2001. tion,º Int'l J. Pattern Recognition and Artificial Intelligence, vol. 15, [32] K. Pearson, ªOn Lines and Planes of Closest Fit to Systems of no. 6, pp. 905-919, 2001. Points in Space,º The Philosophical Magazine, Series 6, vol. 2, no. 11, [7] M.D. Altschuler, B.R. Altschuler, and J. Taboada, ªMeasuring pp. 559-572, 1901. Surfaces Space-Coded by a Laser-Projected Dot Matrix,º Proc. [33] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, SPIE Conf. Imaging Applications for Automated Industrial Inspection Numerical Recipes in C: The Art of Scientific Computing. Cambridge, and Assembly, vol. 182, pp. 187-191, 1979. U.K.: Cambridge Univ. Press, 1988. [8] D.H. Ballard, ªGeneralizing the Hough Transform to Detect [34] P.L. Rosin, ªFitting Superellipses,º IEEE Trans. Pattern Analysis and Arbitrary Shapes,º Pattern Recognition, vol. 13, no. 2, pp. 111-122, Machine Intelligence, vol. 22, no. 7, pp. 726-732, July 2000. 1981. [35] R. Safaee-Rad, K.C. Smith, B. Benhabib, and I. Tchoukanov, [9] A.H. Barr, ªSuperquadrics and Angle-Preserving Transforma- ªApplication of Moment and Fourier Descriptors to the Accurate tions,º IEEE Computer Graphics and Applications, vol. 1, no. 1, pp. 11- Estimation of Elliptical Shape Parameters,º Pattern Recognition 23, 1981. Letters, vol. 13, no. 7, pp. 497-508, 1992. [10] M. Bennamoun and B. Boashash, ªA Structural-Description-Based [36] P.D. Sampson, ªFitting Conic Sections to `Very Scattered' Data: An Vision System for Automatic Object Recognition,º IEEE Trans. Iterative Refinement of the Bookstein Algorithm,º Computer Systems, Man, and Cybernetics, Part B: Cybernetics, vol. 27, no. 6, Graphics and Image Processing, vol. 18, no. 1, pp. 97-108, 1982. pp. 893-906, 1997. [37] F. Solina and R. Bajcsy, ªRecovery of Parametric Models from [11] P.T. Boggs, R.H. Byrd, and R.B. Schnabel, ªA Stable and Efficient Range Images: The Case for Superquadrics with Global Deforma- Algorithm for Nonlinear Orthogonal Distance Regression,º SIAM J. tions,º IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 12, Scientific and Statistical Computing, vol. 8, no. 6, pp. 1052-1078, 1987. no. 2, pp. 131-147, Feb. 1990. 638 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 24, NO. 5, MAY 2002 [38] D. Sourlier, ªThree Dimensional Feature Independent Bestfit in Hyung Suck Cho received the BS degree from Coordinate Metrology,º PhD Thesis no. 11319, ETH Zurich, Seoul National University, Korea in 1971, the MS Switzerland, 1995. degree from Northwestern University Evanston, [39] S. Sullivan, L. Sandford, and J. Ponce, ªUsing Geometric Distance Illinois, in 1973, and the PhD degree from the Fits for 3-D Object Modeling and Recognition,º IEEE Trans. Pattern University of California at Berkeley in 1977. From Analysis and Machine Intelligence, vol. 16, no. 12, pp. 1183-1196, 1977 to 1978, he was a postdoctoral fellow in the Dec. 1994. Department of Mechanical Engineering, Univer- [40] G. Taubin, ªEstimation of Planar Curves, Surfaces, and Nonplanar sity of California at Berkeley. Since 1978, he has Space Curves Defined by Implicit Equations with Applications to been a professor in the Department of Production Edge and Range Image Segmentation,º IEEE Trans. Pattern Engineering, Department of Automation and De- Analysis and Machine Intelligence, vol. 13, no. 11, pp. 1115-1138, sign, Seoul Campus, and is currently with the Department of Mechanical Nov. 1991. Engineering, Korea Advanced Institute of Science and Technology [41] D.A. Turner, ªThe Approximation of Cartesian Coordinate Data by (KAIST), Korea. From 1984 to 1985, he was a visiting scholar at the Parametric Orthogonal Distance Regression,º PhD thesis, School of Fraunhofer Institute for Manufacturing Engineering and Automation (IPA), Computing and Math., Univ. of Huddersfield, U.K., 1999. Germany, where he carried out research on robot-based assembly. He È [42] K. Voss and H. Suûe, ºInvariant Fitting of Planar Objects by has been invited as a short-term visiting scholar to several universities: Primitives,º IEEE Trans. Pattern Analysis and Machine Intelligence, Ritsumeikan University, Japan, the University of Paderborn, Germany, vol. 19, no. 1, pp. 80-84, Jan. 1997. and the New Jersey Institute of Technology, in 1987, 1992, 1998, È [43] K. Voss and H. Suûe, ªA New One-Parametric Fitting Method for respectively. From 1995 to 1996, as a visiting professor, he participated in Planar Objects,º IEEE Trans. Pattern Analysis and Machine graduate school education for the Advanced Manufacturing Program Intelligence, vol. 21, no. 7, pp. 646-651, July 1999. (AMP) of the University of California, San Diego. Dr. Cho's research [44] F.M. Wahl, ªA Coded-Light Approach for 3-Dimensional (3D) interests are focused on the area of environment perception and Vision,º Research Report no. 1452, IBM Zurich Research Labora- recognition for mobile robots, machine vision and pattern classification, tory, Switzerland, Sept. 1984. and application of artificial intelligence/machine intelligence. During his research career, he has published six book chapters and more than 346 Sung Joon Ahn received the BS degree in research papers (276 papers in international journals and conferences mechanical design and production engineering and 66 papers in Korean journals). He now serves on the editorial boards from Seoul National University in South Korea in of five international journals, including Control Engineering Practice 1985. He received the MS degree in production (IFAC). In addition to his academic activities, he serves and/or has served engineering from the Korea Advanced Institute of on a few technical committees (TC) of IFAC and IMEKO: TC on Robotics, Science and Technology (KAIST) in 1987. He TC on Advanced Manufacturing (IFAC), and Measurements on Robotics worked as a junior research scientist at the (IMEKO). He has been general chair and/or cochair of several interna- research center of LG Electronics from 1987 to tional conferences, including SPIE-Opto Mechatronic Systems (2000, 1990. Since 1991, he has been working as a 2001) and IEEE/RSJ IROS (1999). For his achievements in research research scientist for the Fraunhofer Institute for works in the fields of robotics and automation, Dr. Cho has been endowed Manufacturing Engineering and Automation (IPA) in Stuttgart, Germany. with a POSCO professorship from Pohang Steel Company, Korea, since His research interests include pattern recognition, optical 3D-measure- 1995. In 1998, he was awarded the Thatcher Brothers Prize from the ment, close-range photogrammetry, camera calibration, and 3D-informa- Institution of Mechanical Engineers, United Kingdom, and he also tion processing. received the fellowship of Alexander von Humboldt of Germany in 1984. He is a member of the IEEE. Wolfgang Rauh received the Doctorate degree È Hans-Jurgen Warnecke studied mechanical in 1993 from the University of Stuttgart. He engineering at the Technical University of worked at the University of Stuttgart from 1984 Braunschweig, Germany. He was a research until 1989 after having received the Diploma engineer at the Institute for Machine Tools and degree in mechanical engineering from the Manufacturing afterward, from 1965 to 1970, a University of Karlsruhe in 1984. He then became managing director for planning and manufactur- head of the Industrial Image Processing Group ing at the Rollei-Werke, Braunschweig, a camera at the Fraunhofer Institute for Manufacturing factory. Since 1971, he has been a full professor Engineering and Automation (IPA). In 1991, he for industrial manufacturing and management, was made head of the Information Processing University of Stuttgart, and the head of the Department at the same institute. In this department, a wide variety of Fraunhofer Institute for Manufacturing Engineering and Automation projects in the area of image processing are carried out. (IPA) of the Fraunhofer Society for Applied Research. His main fields of research and development are management and organization of industrial companies, computer integrated manufacturing, computer- aided planning of material flow, simulation, production planning and control, manufacturing processes, flexible manufacturing systems, industrial robots, automation of handling and assembly, quality control, industrial measurement, cleanroom manufacturing, semiconductor manufacturing automation, and maintenance. He has published 10 books with many publications in the field of production, industrial engineering, and automation. Since 1 October 1993, he has been the president of the Fraunhofer Society, the largest institution for applied research and development in Europe with a budget of EUR 0.9 billion, of which EUR 300 million are proceeds from 3,000 private enterprises. There are 11,000 employees in 56 institutes. In the years 1995 to1997, he was the president of VDI, the Association of German Engineers, the largest institution of this kind in Europe. . For more information on this or any other computing topic, please visit our Digital Library at http://computer.org/publications/dlib.

DOCUMENT INFO

Shared By:

Categories:

Tags:
Least Squares, in Space, parametric curves, fitting algorithms, Sung Joon, Data Points, Curve fitting, Implicit Curves, Lecture Notes, J. Ahn

Stats:

views: | 50 |

posted: | 5/17/2011 |

language: | English |

pages: | 19 |

OTHER DOCS BY ert634

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.