Orthogonal distance fitting of implicit curves and surfaces

Document Sample
Orthogonal distance fitting of implicit curves and surfaces Powered By Docstoc
					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 …a†gm :
                                                                                                     H         iˆI              …S†
                F …b; X† ˆ       bj Xkj Y lj Z mj ˆ H;                                       aH m
                                                                                                 fXi giˆI 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
                                                                                                                  iˆI
                                                                                        each given point fXi gm . For a given point xi in frame xyz
                                                                                                                iˆI

                                                                                                                  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 À V†y ˆ 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 …x†k ; xk‡I ˆ xk ‡ Áx:          …IQ†
            @x                                                                                         
                    k                                                                PI ‡ H rf        Áx
If necessary, especially in order to prevent the point update                           r„f      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
                                                                                                         iˆI
iteration cycles with most fitting tasks) and the orthogonal          feature from each given point fXi gm are found and now
                                                                                                            iˆI
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 XˆXH
                                                                                                                          i
                                                                                                   À          Á
                                                                                                             H „                            
                                                                                                     Xi À Xi            @x @RÀI        @Xo    
                                                                                                ˆÀ                RÀI       ‡    ‰xŠ ‡
                                                                                                    kXi À XHi k         @a    @a        @a xˆxH
                                                                                                                                                 i
                                                                                                   À        Á
                                                                                                           H „    
                                                                                                     xi À xi @x  
                                                                                                ˆÀ
                                                                                                    kxi À xHi k @axˆxH
                                                                                                   À          Á 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 †   == rfjxˆxH 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                         r„f   ˆ      ˆÀ
                                                                                       @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                                                     xˆxHi
                                                                                              À        Á„                                 
                                                                 …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 ;    ak‡I ˆ 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 ;             ak‡I ˆ 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 k‡I % 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
                                                                                                                                   XˆXi
   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 ; ak‡I ˆ 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 k‡I % 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 XˆXH         @a   @a         @a xˆxH                    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 Š :
                  @axˆxH              @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 ˆ di—g…wI ; F F F ; wq † ;
                f™j …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
        @f™j           À                Á                                                              PJ        PJ
  w™j         Áa ˆ Àw™j f™j …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                                                          ;                                           iˆI
                                                                                                              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                                                                                       m‡nÀ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=a†P=" ‡…jyj=b†P=" À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 t—n… =P††P ˆ 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 t—n… =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;™ir™le =m;      rP ˆ r™ir™le :
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=a†P="I ‡…jyj=b†P="I         ‡…jzj=c†P="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 iˆI 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 …m—x…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; m—x…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; m—x…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 À y†x=aP À …xi À x†y=bP
estimation accuracy is higher than 0.1 m for length unit and
0.1 r—d 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 À y†x=aP À …xi À x†y=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=a†P="    —nd    B ˆ …jyj=b†P=" :
                       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
                                                      xˆxi
                                              
                          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=a†P="I =d;     B ˆ …jyj=b†P="I =d; C ˆ …jzj=c†P="P ;
                  @f       H        H           I=aP   H
                     ˆ                                                                 D ˆ d"I = "P ;   with d ˆ …jxj=a†P="I ‡…jyj=b†P="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:
Stats:
views:50
posted:5/17/2011
language:English
pages:19