Rationalising the Renormalisation Method of Kanatani.
The renormalisation technique of Kanatani is intended to iteratively minimise a cost function of a certain form while avoiding systematic bias inherent in the common method of minimisation due to Sampson. Within the computer vision community, the technique has generally proven difficult to absorb. This work presents an alternative derivation of the technique, and places it in the context of other approaches. We first show that the minimiser of the cost function must satisfy a special variational equation. A Newton-like, fundamental numerical scheme is presented with the property that its theoretical limit coincides with the minimiser. Standard statistical techniques are then employed to derive afresh several renormalisation schemes. The fundamental scheme proves pivotal in the rationalising of the renormalisation and other schemes, and enables us to show that the renormalisation schemes do not have as their theoretical limit the desired minimiser. The various minimisation schemes are finally subjected to a comparative performance analysis under controlled conditions.
Introduction
Many problems in computer vision are readily formulated as the need to minimise a
cost function with respect to some unknown parameters. Such a cost function will
often involve (known) covariance matrices characterising uncertainty of the data and
will take the form of a sum of quotients of quadratic forms in the parameters. Finding
the values of the parameters that minimise such a cost function is often difficult.
One approach to minimising a cost function represented as a sum of fractional expressions
is attributed to Sampson. Here, an initial estimate is substituted into the denominators
of the cost function, and a minimiser is sought for the now scalar-weighted
numerators. This procedure is then repeated using the newly obtained estimate until
convergence is obtained. It emerges that this approach is biased. Noting this, Kenichi
Kanatani developed a renormalisation method whereby an attempt is made at each iteration
to undo the biasing effects. Many examples may be found in the literature of
problems benefiting from this approach.
In this work, we carefully analyse the renormalisation concept, and place it in the
context of other approaches. We first specify the general problem form and an associated
cost function to which renormalisation is applicable. We then show that the
cost function minimiser must satisfy a particular variational equation. Interestingly,
we observe that the renormalisation estimate is not a theoretical minimiser of the cost
function, and neither are estimates obtained via some other commonly used methods.
This is in contrast to a fundamental numerical scheme that we present.
New derivations are given for Kanatani's first-order and second-order renormalisation
schemes, and several variations on the theme are proposed. This serves as a
rationalising of renormalisation, making recourse to various statistical concepts. Experiments
are carried out on the benchmark problem of estimating ellipses from synthetic
data points and their covariances. The renormalisation schemes are shown to
perform better than more traditional methods in the face of data exhibiting noise that
is anisotropic and inhomogeneous. None of the methods outperforms the relatively
straightforward fundamental numerical scheme.
Problem Formulation
A wide class of computer vision problems may be couched in terms of an equation of
the form
Here is a vector representing unknown parameters;
is a vector representing an element of the data (for example, the locations of a pair of
corresponding points); and is a vector with the data transformed
in such a way that: (i) each component is a quadratic form in the compound
vector one component of u(x) is equal to 1. An ancillary constraint may
also apply that does not involve the data, and this can be expressed as
for some scalar-valued function . The estimation problem can now be stated as fol-
lows: Given a collection image data, determine  6= 0 satisfying (2)
such that (1) holds with x replaced by x i for 1  i  n. When n > l and noise is
present, the corresponding system of equations is overdetermined and as such may fail
to have a non-zero solution. In this situation, we are concerned with finding  that best
fits the data in some sense. The form of this vision problem involving (known) co-variance
information was first studied in detail by Kanatani [12], and later by various
others (see, e.g., [4, 13,14,20, 21]).
Conic fitting is one problem of this kind [2, 23]. Two other conformant problems
are estimating the coefficients of the epipolar equation [6], and estimating the coefficients
of the differential epipolar equation [3, 22]. Each of these problems involves an
ancillary cubic constraint. The precise way in which these example problems accord
with our problem form is described in a companion work [4].
3 Cost Functions and Estimators
A vast class of techniques for solving our problem rest upon the use of cost functions
measuring the extent to which the data and candidate estimates fail to satisfy (1). If-
for simplicity-one sets aside the ancillary constraint, then, given a cost function
corresponding estimate b  is defined by
6=0
Since (1) does not change if  is multiplied by a non-zero scalar, it is natural to demand
that b  should satisfy (3) together with all the  b
, where  is a non-zero scalar. This is
guaranteed if J is -homogeneous:
Henceforth only -homogeneous cost functions will be considered. The assignment
of b  (uniquely defined up to a scalar factor) to x will be termed the J-based
estimator of .
Once an estimate has been generated by minimising a specific cost function, the
ancillary constraint (if it applies) can further be accommodated via an adjustment pro-
cedure. One possibility is to use a general scheme delivering an 'optimal correction'
described in [12, Subsec. 9.5.2]. In what follows we shall confine our attention to the
estimation phase that precedes adjustment.
3.1 Algebraic Least Squares Estimator
A straightforward estimator is derived from the cost function
JALS
where A
l Here each summand  T A i
is the square of the algebraic distance j T u(x i )j. Accordingly, the JALS -based estimate
of  is termed the algebraic least squares (ALS) estimate and is denoted b ALS . It
is uniquely determined, up to a scalar factor, by an eigenvector of
associated
with the smallest eigenvalue [4].
3.2 Approximated Maximum Likelihood Estimator
The ALS estimator treats all data as being equally valuable. When information about
the measurement errors is available, it is desirable that it be incorporated into the estimation
process. Here we present an estimator capable of informed weighting. It is
based on the principle of maximum likelihood and draws upon Kanatani's work on
geometric fitting [12, Chap. 7].
The measurement errors being generally unknowable, we regard the collective data
sample value taken on by an aggregate of vector-valued random
variables We assume that the distribution of exactly
specified but is an element of a collection fP  j  2 Hg of candidate distributions,
with H the set of all (n
The candidate distributions are to be such that if a distribution P  is in effect, then each
n) is a noise-driven, fluctuating quantity around x i .
We assume that the data come equipped with a collection
positive
definite k  k covariance matrices. These matrices constitute repositories of prior
information about the uncertainty of the data. We put the x i
in use by assuming that,
for each  2 H, P is the unique distribution satisfying the following conditions:
for any the random vectors x i and x j (or equivalently,
the noises behind x i and x j ) are stochastically independent;
for each normal distribution
with mean value vector x i and covariance matrix x i
, that is:
Each distribution P will readily be described in terms of a probability density
function (PDF) (~ x
Resorting to the principle of maximum
likelihood, we give the greatest confidence to that choice of  for which the
likelihood function  7! f(x attains a maximum. Using the explicit
form of the PDF's involved, one can show that the maximum likelihood estimate is the
parameter b  ML at which the cost function
attains a minimum [4, 12]. Each term in the above summation represents the squared
Mahalanobis distance between x i and x i . Note that the value of b
ML remains unchanged
if the covariance matrices are multiplied by a common scalar.
The parameter naturally splits into two parts:  1
These parts encompass the principal parameters and nuisance
parameters, respectively. We are mostly interested in the  1 -part of b  ML , which we
call the maximum likelihood estimate of  and denote b ML . It turns out that b
ML
can be identified as the minimiser of a certain cost function which is directly derivable
from JML . This cost function does not lend itself to explicit calculation. However, a
tractable approximation [4] can be derived in the form of the function
where @x vector y and any k  k
matrix , we let
and next, for each
then JAML can be simply
written as
The JAML -based estimate of  will be called the approximated maximum likelihood
(AML) estimate and will be denoted b AML .
It should be observed that JAML can be derived without recourse to principles of
maximum likelihood by, for example, using a gradient weighted approach that also
incorporates covariances. Various terms may therefore be used to describe methods
that aim to minimise a cost function such as JAML , although some of the terms may
not be fully discriminating. Candidate labels include 'heteroscedastic regression' [13],
'weighted orthogonal regression' [1, 9], and `gradient weighted least squares' [24].
3.3 Variational Equation
Since b AML is a minimiser of JAML , we have that, when the following
equation holds:
Here @  JAML denotes the row vector of the partial derivatives of JAML with respect to
. We term this the variational equation. Direct computation shows that
where X  is the symmetric matrix
Thus (7) can be written as
This is a non-linear equation and is unlikely to admit solutions in closed form.
Obviously, not every solution of the variational equation is a point at which the
global minimum of JAML is attained. However, the solution set of the equation provides
a severely restricted family of candidates for the global minimiser. Within this
set, the minimiser is much easier to identify.
4 Numerical schemes
Closed-form solutions of the variational equation may be infeasible, so in practice
b AML has to be found numerically. Throughout we shall assume that b AML lies close
to b ALS . This assumption is to increase the chances that any candidate minimiser obtained
via a numerical method seeded with b ALS coincides with b AML .
4.1 Fundamental Numerical Scheme
A vector  satisfies (9) if and only if it falls into the null space of the matrix X  .
Thus, if  k 1 is a tentative guess, then an improved guess can be obtained by picking
a vector  k from that eigenspace of X k 1
which most closely approximates the null
space of X  ; this eigenspace is, of course, the one corresponding to the eigenvalue
closest to zero. It can be proved that as soon as the sequence of updates converges,
the limit is a solution of [4]. The fundamental numerical scheme implementing the
above idea is presented in Figure (1). The algorithm can be regarded as a variant of the
Newton-Raphson method.
1. Set
2. Assuming that  k 1 is known, compute the matrix X k 1
3. Compute a normalised eigenvector of X k 1
corresponding to
the eigenvalue closest to zero and take this eigenvector for  k .
4. If  k is sufficiently close to  k 1 , then terminate the procedure;
otherwise increment k and return to Step 2.

Figure

1: Fundamental numerical scheme.
4.2 Sampson's Scheme
Let
and let
be the modification of JAML (; x which the variable  in the denominators
of all the contributing fractions is "frozen" at the value . For simplicity, we
abbreviate J 0
was the first to propose a scheme aiming to minimise a function
involving fractional expressions, such as JAML (although his cost functions did not
incorporate covariance matrices). Sampson's scheme (SMP) applied to JAML takes
b ALS for an initial guess  0 , and given  k 1 generates an update  k by minimising
the cost function  7! J 0
Assuming that the sequence f k g converges,
the Sampson estimate is defined as b . Note that each function
AML (;  k 1 ) is quadratic in . Finding a minimiser of such a function is straight-
forward. The minimiser  k is an eigenvector of M  k 1
corresponding to the smallest
eigenvalue; moreover, this eigenvalue is equal to J 0
scheme is summarised in Figure (2).
1. Set
2. Assuming that  k 1 is known, compute the matrix M  k 1
3. Compute a normalised eigenvector of M  k 1
corresponding to
the smallest (non-negative) eigenvalue and take this eigenvector
for  k .
4. If  k is sufficiently close to  k 1 , then terminate the procedure;
otherwise increment k and return to Step 2.

Figure

2: Sampson's scheme.
A quick glance shows that this scheme differs from the fundamental numerical
scheme only in that it uses matrices of the form M  instead of matrices of the form
. These two types of matrix are related by the formula X
Letting k !1 in (11) and taking into account
the equality J 0
we see that b SMP satisfies
[M  JAML ()I l
where I l is the l  l identity matrix. We call this the Sampson equation. Note that it is
different from the variational equation (9) and that, as a result, b SMP is not a genuine
minimiser of JAML .
Renormalisation
The matrices M  E  and M  JAML ()I l underlying the variational equation
and the Sampson equation (12) can be viewed as modified or normalised forms
of M  . As first realised by Kanatani [12], a different type of modification can be
proposed based on statistical considerations. The requirement is that the modified or
renormalised M  be unbiased in some sense. Using the renormalised M  , one can
formulate an equation analogous to both the variational and Sampson equations. This
equation can in turn be used to define an estimate of .
Here we rationalise the unbiasing procedure under the condition that noise in an
appropriate statistical model is small. In the next section, various schemes will be presented
for numerically computing the parameter estimate defined in this procedure. A
later section will be devoted to the derivation of an unbiasing procedure appropriate for
noise that is not necessarily small, and to the development of schemes for numerically
computing the parameter estimate defined in this more general procedure.
5.1 Unbiasing M
Regard the given data value taken on by the random vectors
introduced earlier. Suppose that
Form the following random version of M
with 'true' value
In view of (4), A(x i
0: On the other hand, since each rank-one matrix A(x i ) is non-negative definite, and
since also each B(x
) is non-negative definite 1 , M  is non-negative definite. As
the A(x i ) are independent, M is generically positive definite, with E
0:
Thus on average  T M   does not attain its 'true' value of zero, and as such is biased.
The bias can be removed by forming the matrix
The terms E
can be calculated explicitly. There is a matrix-valued function
to be specified later, such that, for each
As
Thus B(x; ) is not positive definite, but merely non-negative definite.
The unbiased M can be written as
The random matrix Y  is a raw model for obtaining a fully deterministic modification
of M  . For each
Guided by (14), we take
for a modified M  . Somewhat surprisingly, this choice turns out not to be satisfactory.
The problem is that while the A i do not change when the x i
are multiplied by a common
scalar, the D i do change. A properly designed algorithm employing a modified
should give as outcomes values that remain intact when all the x i
are multiplied
by a common scalar. This is especially important if we aim not only to estimate the
parameter, but also to evaluate the goodness of fit. Therefore further change to the
numerators of the fractions forming Y  is necessary.
The dependence of D(x;) on  is fairly complex. To gain an idea of what
needs to be changed, it is instructive to consider a simplified form of D(x;). A
first-order (in some sense) approximation to D(x;) is, as will be shown shortly, the
defined in (6). The dependence of B(x; ) on  is simple: if  is
multiplied by a scalar, then B(x;) is multiplied by the same scalar. This suggests
we introduce a compensating factor J com (; x com () in short, with the
property that if the x i
are multiplied by a scalar, then J com () is multiplied by the
inverse of this scalar. With the help of J com (), we can form, for each
renormalised numerator  T A i  J com ()  T B i  and can next set
where M  is given in (10) and N  is defined by
The numerators in (16) are clearly scale invariant. Note that J com () plays a role
similar to that played by the factors ( T A i )=( the formula for X  given
in (8). Indeed, if the x i
are multiplied by , then so are the B i , and consequently
the are multiplied by  1 . The main difference between J com ()
and the ( T A i )=( latter fractions change with the index i, while
J com () is common for all the numerators involved.
To find a proper expression for J com (), we take a look at X  for inspiration. Note
that, on account of (8),  T X  0: By analogy, we demand that  T Y  0: This
equation together with (16) implies that
T N
=n
It is obvious that J com () thus defined has the property required of a compensating fac-
tor. Moreover, this form of J com () is in accordance with the unbiasedness paradigm.
Indeed, if we form the random version of J com
insofar as E
abbreviating J com (; x
to J com (),
and further
0:
We see that Y  given by
is unbiased, which justifies the design of Y  .
Since, in view of (19), J com () is equal to 1 in the mean, the difference between
and
is blurred on average. Thus the refined renormalisation based on (16) is close in spirit
to our original normalisation based on (15).
5.2 Renormalisation Equation
The renormalisation equation
is an analogue of the variational and Sampson equations alike. It is not naturally derived
from any specific cost function, and, as a result, it is not clear whether it has
any solution at all. A general belief is that in the close vicinity of b ALS there is a
solution and only one. This solution is termed the renormalisation estimate and is denoted
REN . Since the renormalisation equation is different from the variational and
distinct both from b AML and b SMP . It should be stressed
that the difference between b REN and b AML may be unimportant as both these estimates
can be regarded as first-order approximations to b ML and hence are likely to be
statistically equivalent.
In practice, b REN is represented as the limit of a sequence of successive approximations
to what b REN should be. If, under favourable conditions, the sequence is
convergent, then the limit is a genuine solution of (20). Various sequences can be taken
to calculate b REN in this way. The simplest choice results from mimicking the fundamental
numerical scheme as follows. Take b ALS to be an initial guess  0 . Suppose that
an update  k 1 has already been generated. Form Y  k 1
, compute an eigenvector of
corresponding to the eigenvalue closest to zero, and take this eigenvector for  k .
If the sequence f k g converges, take the limit for b REN . As we shall see shortly, b REN
thus defined automatically satisfies (20). This recipe for calculating b REN constitutes
what we term the renormalisation algorithm.
6 First-Order Renormalisation Schemes
First-order renormalisation is based on the formula
as already pointed out in Subsection 5.1. To justify this formula, we retain the sequence
of independent random vectors as a model for our data
assume that distribution P for some
make a fundamental assumption to the effect that the noise driving each x i is small.
For simplicity, denote x i by x, and contract x i
to . Since the noise driving x is
small, we can replace u(x) by the first-order sum in the Taylor expansion about x,
next, taking into account that  T 0, we can write
Hence
and further, in view of (5) and (6),
which, on account of (13), establishes (21).
With the formula (21) validated, we can safely use Y  in the form given in (16)
(with J com given in (18)). The respective renormalisation estimate will be called the
first-order renormalisation estimate and will be denoted b REN1 .
6.1 The FORI Scheme
By introducing an appropriate stopping rule, the renormalisation algorithm can readily
be adapted to suit practical calculation. In the case of first-order renormalisation, the
resulting method will be termed the first-order renormalisation scheme, Version I, or
simply the FORI scheme. It is given in Figure (3).
1. Set
ALS .
2. Assuming that  k 1 is known, compute the matrix Y  k 1
using
(16).
3. Compute a normalised eigenvector of Y  k 1
corresponding to the
eigenvalue closest to zero and take this eigenvector for  k .
4. If  k is sufficiently close to  k 1 , then terminate the procedure;
otherwise increment k and return to Step 2.

Figure

3: First-order renormalisation scheme, Version I.
6.2 The FORII Scheme
The FORI scheme can be slightly modified. The resulting first-order renormalisation
scheme, Version II, or the FORII scheme, is effectively one of two schemes proposed
by Kanatani [12, Chap. 12] (the other concerns second-order renormalisation).
Introduce a function
com
Abbreviating J 0
com
com (; ), let
com (; )B i
com (; )N  :
It is immediately verified that
for each  and each  6= 0. We also have that
com
and
As previously, take b ALS to be an initial guess  0 . Suppose that an update  k 1 has
already been generated. Then
be a normalised eigenvector of Y  k 1
cor-
responding to the smallest eigenvalue  k . In view of (25), the update Y  k
is straight-forwardly
generated from the updates M  k
, and c k . It turns out that, under a
certain approximation, c k can be updated directly from c k 1 rather than by appealing
to the above formula.
We have
Substituting  k for  and  k 1 for  in (22), we obtain
The last two equations imply that
Now, assume that J 0
since both terms are close to
J com ( k ), this is a realistic assumption. Under this assumption we have c
and equation (26) becomes
This is a formula for successive updating of the c k . Defining consecutive Y k with
the help of M  k
and c k as in (25), and proceeding as in the FORI scheme, we
obtain a sequence f k g. If it converges, we take the corresponding limit for b REN1 . It
can be shown that b REN1 thus defined satisfies the renormalisation equation. The first-order
renormalisation scheme, Version II, or the FORII scheme, based on the above
algorithm is summarised in Figure (4).
6.3 The FORIII Scheme
With the help of the function J 0
com , yet another defining sequence can be constructed.
Take b ALS for an initial guess  0 . Suppose that an update  k 1 has already been
generated. Define  k to be the minimiser of the function  7! J 0
com (;  k 1
6=0
com (;  k 1
Since
com
ALS and c
2. Assuming that  k 1 and c k 1 are known, compute the matrix
3. Compute a normalised eigenvector of M  k 1
cor-
responding to the smallest eigenvalue  k and take this eigenvector
for  k . Then define c k by (27).
4. If  k is sufficiently close to  k 1 , then terminate the procedure;
otherwise increment k and return to Step 2.

Figure

4: First-order renormalisation scheme, Version II.
where
T N
com (; )N
T N
it follows that  k satisfies
Assuming that the sequence f k g converges, let b
clearly,
b REN1 satisfies
which is an equation equivalent to (20). Note that in this method b REN1 is defined as
the limit of a sequence of minimisers of cost functions. As such the algorithm is similar
to Sampson's algorithm, but the latter, of course, uses different cost functions.
The minimisers  k can be directly calculated. To see this, rewrite (28) as
com
We see that  k is an eigenvector and J 0
com (;  k 1 ) is a corresponding eigenvalue of
the linear pencil P defined by
( a real
If  is any eigenvector of P k 1 with eigenvector
necessarily, J 0
com
com
conclude that J 0
is an eigenvector of P k 1 corresponding
.
2. Assuming that  k 1 is known, compute the matrices M  k 1
and N  k 1
3. Compute a normalised eigenvector of the eigenvalue problem
corresponding to the smallest eigenvalue and take this eigenvector for  k .
4. If  k is sufficiently close to  k 1 , then terminate the procedure; otherwise
increment k and return to Step 2.

Figure

5: First-order renormalisation scheme, Version III.
to the smallest eigenvalue. This observation leads to the first-order renormalisation,
Version III, or the FORIII scheme, given in Figure 5. The matrices N  k 1
are singular,
so the eigenvalue problem for P k 1 is degenerate. A way of reducing this problem
to a non-degenerate one, based on the special form of the matrices M  and N  , is
presented in [4].
7 Second-Order Renormalisation
Second-order renormalisation rests on knowledge of the exact form of D(x). Here
we first determine this form and next use it to evolve a second-order renormalisation
estimate and various schemes for calculating it.
7.1 Calculating D(x;)
Determining the form of D(x;) is tedious but straightforward. We commence by
introducing some notation.
be the vector of variables. Append to this vector a unital
component, yielding
be the vector of carriers. Given the special form of
u(x) as described in Section 2, each u
can be expressed as
where K
In what follows we adopt
Einstein's convention according to which the summation sign for repeated indices is
omitted. With this convention, equation (29) becomes
Let x be a random Gaussian k vector with mean x and covariance matrix
Clearly, y has
for the mean, and the defined by
for the covariance matrix.
Since
we have
Now
By a standard result about moments of the multivariate normal distribution,
In view of (31),
and so
Hence
be the l  l matrices defined by
d
d 2;
With these matrices, (32) can be written as
Hence
which is the desired formula.
7.2 Redefining J com () and Y
We retain the framework of Subsection 5.1, but use the full expression for D(x;) instead
of the first-order approximation. We aim to modify, for each
T A i  into a term similar to  T A i   T D i  (recall that D
)),
remembering the need for suitable compensation for scale change. The main problem
now is that D(x;) does not change equivariantly with . Under the scale change
7! , the two components of D(x;) defined in (34) undergo two different
transformations:
a solution as follows. We introduce a compensating factor J com (; x
J com () in short, with the property that if the x i
are multiplied by , then J com () is
multiplied by  1 . We place this factor in front of D
its square in
front of D
), and form a modified numerator as follows:
This numerator is obviously invariant with respect to scale change. In analogy to (17),
we introduce
and, in analogy to (16), let
Demanding again that  T Y  we obtain the following quadratic equation for
J com ():
This equation has two solutions
com
As M  is positive definite, and N 1; and N 2; are non-negative definite, we have
com ()  0. Since the compensating factor used in the first-order
renormalisation is non-negative, we take, by analogy, J
com () to be a compensating
factor and denote it by J com (); thus
J com
Multiplying both numerator and denominator of J com by [( T N
we see that
J com
If  T N 2;    T M   is small compared to ( T N then we may readily infer
that
J com ()
T N
This expression is very similar to formula (18) for J com (), which indicates that the
solution adopted is consistent with the first-order renormalisation.
Inserting J com () given in (38) into (37), we obtain a well-defined expression for
Y  . We can now use it to define a renormalisation estimate using the renormalisation
equation (20). We call this estimate the second-order renormalisation estimate and
denote it b
REN2 .
7.3 The SORI Scheme
Mimicking the FORI scheme, we can readily advance a scheme for numerically finding
b REN2 . We call this the second-order renormalisation scheme, Version I, or the SORI
scheme. Its steps are given in Figure (6).
1. Set
2. Assuming that  k 1 is known, compute the matrix Y  k 1
using (37) and
(38).
3. Compute a normalised eigenvector of Y  k 1
corresponding to the eigenvalue
closest to zero and take this eigenvector for  k .
4. If  k is sufficiently close to  k 1 , then terminate the procedure; otherwise
increment k and return to Step 2.

Figure

renormalisation scheme, Version I.
7.4 The SORII Scheme
The SORI scheme can be modified in a similar way to that employed with the FORI
scheme. The resulting second-order renormalisation scheme, Version II, or the SORII
scheme, is effectively the second of the two schemes originally proposed by Kanatani.
Introduce
com
Abbreviating J 0
com
com (; ), let
com
It is immediately verified that equations (22) and (24) hold, as does J com
com (; );
the counterpart of (23).
Again, take b ALS to be an initial guess  0 . Suppose that an update  k 1 has already
been generated. Note that
where
Let  k be a normalised eigenvector of Y  k 1
corresponding to the smallest eigenvalue
k . We intend to find an update c k appealing directly to c k 1 . To this end, observe that
Substituting  k for  and  k 1 for  in (22), taking into account (41), assuming
explained analogously when deriving the FORII
scheme), and taking into account that, by (43), J 0
Combining this equation with (44) yields
Taking into account that c 2
we can rewrite (45) as the quadratic constraint on c k 1 given by
Let
Equation (46) has two solutions
which are real when D k 1  0.
Suppose that D k 1  0. If c k were directly defined by (43), it would be non-
negative. It is therefore reasonable to insist that c k obtained by updating c k 1 also be
non-negative. This requirement can be met by setting c next by ensuring that
k. For this reason, we select  1 to be c k 1 , obtaining
To treat the case D k 1 < 0, we first multiply the numerator and denominator of the
fractional expression in (48) by
obtaining
Next we note that if  k
k is small compared to  k
and further
This formula is very similar to (27). We use it with the equality sign instead of the
approximation sign to generate c k in the case D k 1 < 0.
In this way, we arrive at the following update formula:
The SORII scheme can now be formulated as in Figure (7).
1. Set
2. Assuming that  k 1 and c k 1 are known, compute the matrix Y  k 1 ;k 1
by using (42).
3. Compute a normalised eigenvector of Y k 1 ;k 1
corresponding to the
smallest eigenvalue  k and take this eigenvector for  k . Then define c k by
using (49).
4. If  k is sufficiently close to  k 1 , then terminate the procedure; otherwise
increment k and return to Step 2.

Figure

7: Second-order renormalisation scheme, Version II.
7.5 The SORIII Scheme
The estimate b REN2 can be represented as a limit of a sequence of minimisers of
cost functions as follows. Take b ALS for an initial guess  0 . Suppose that an up-date
has already been generated. Define  k to be the minimiser of the function
com (;  k 1
6=0
com (;  k 1
Assuming that the sequence f k g converges, take lim k!1  k for b REN2 . It can readily
be shown that b REN2 thus defined satisfies Y  It also can be shown that each
com
com
Here  k is an eigenvector and J 0
corresponding eigenvalue of the
defined by
( a real
In fact,  k is an eigenvector of P k 1 corresponding to the smallest eigenvalue. This observation
leads to the second-order renormalisation scheme, Version III, or the SORIII
scheme, given in Figure (8).
1. Set
2. Assuming that  k 1 is known, compute the matrices M  k 1
and N 2;k 1
3. Compute a normalised eigenvector of the eigenvalue problem
corresponding to the smallest eigenvalue and take this eigenvector for  k .
4. If  k is sufficiently close to  k 1 , then terminate the procedure; otherwise
increment k and return to Step 2.

Figure

8: Second-order renormalisation scheme, Version III.
The eigenvalue problem for a quadratic pencil can readily be reduced to the eigenvalue
problem for a linear pencil. Indeed,  and  satisfy
if and only if there exists  0 such that
I l
I l 0
in which case, necessarily,  . The matrices N 1 and N 2 appropriate for the
SORIII scheme are non-negative definite, but not necessarily positive definite. A way
of reducing the problem (51) to a similar problem involving positive definite matrices
2 is given in [4] This method takes advantage of the special form of the
matrices M  , N 1; , and N 2; .
8 Experimental Results
The previously derived algorithms were tested on the problem of conic fitting, which
constitutes a classical benchmark problem in the literature [2, 5, 7, 8, 10,11, 15-19, 23].
Specifically, the fitting algorithms were applied to contaminated data arising from a
portion of an ellipse.
Synthetic testing is employed here as this enables precise control of the nature of
the data and their associated uncertainties. Results obtained in real world testing, in
applications domains described earlier, are presented in subsequent work.
Our tests proceeded as follows. A randomly oriented ellipse was generated such
that the ratio of its major to minor axes was in the range [2; 3], and its major axis was
approximately 200 pixels in length. One third of the ellipse's boundary was chosen as
the base curve, and this included the point of maximum curvature of the ellipse. A set
of true points was then randomly selected from a distribution uniform along the length
of the base curve.
For each of the true points, a covariance matrix was randomly generated (using
a method described below) in accordance with some chosen average level of noise,
. The true points were then perturbed randomly in accordance with their associated
covariance matrices, yielding the data points. In general, the noise conformed to an
inhomogeneous and anisotropic distribution. Figure 9 shows a large ellipse, some selected
true points, a small ellipse for each of these points and the data points. Each
of the smaller ellipses represents a level set of the probability density function used
to generate the datum, and as such captures graphically the nature of the uncertainty
described by its covariance matrix.

Figure

9: True ellipse, data, and associated covariance ellipses
The following procedure was adopted for generating covariance matrices associated
with image points, prescribing (anisotropic and inhomogeneous) noise at a given
average level . The scale  of a particular covariance matrix was first selected from
a uniform distribution in the range [0; 2]. (Similar results were obtained using other
distributions about .) Next, a skew parameter  was generated from a uniform distribution
between 0 and 0:5. An intermediate covariance matrix was then formed by
setting
This matrix was then 'rotated' by an angle
selected from a uniform distribution between
and 2 to generate the final covariance
with
O
cos
sin
sin
cos
Let Tr A denote the trace of the matrix A. Since
is clear that the above procedure ensures that E [Tr
With the data points and their associated covariances prepared, each method under
test was then challenged to determine the coefficients of the best fitting conic. (Note,
therefore, that it was not assumed that the conic was an ellipse.) The methods were
supplied with the data points, and if a specific method was able to utilise uncertainty
information, it was also supplied with the data points' covariance matrices. Then estimates
were generated, and for each of these a measure of the error was computed
using a recipe given below. Testing was repeated many times using newly generated
data points (with the covariance matrices and true data points remaining intact). The
average errors were then displayed for each method.
The error measure employed was as follows. Assume that a particular method
has estimated an ellipse. The error in this estimate was declared to be the sum of the
shortest (Euclidean) distances of each true point from the estimated ellipse. Note that
this measure takes advantage of the fact that the underlying true points are known. Were
these unknown, an alternative measure might be the sum of the Mahalanobis distances
from the data points to the estimated ellipses.
The methods tested were as follows:
least squares scheme,
renormalisation scheme 1,
renormalisation scheme 2,
renormalisation scheme 3,
renormalisation scheme 1,
renormalisation scheme 2,
renormalisation scheme 3,

Table

1 shows the average error obtained when each method was applied to 500
sets of data points, with  varying from 1 to 10 pixels in steps of 1. Each set of data
points was obtained by perturbing shows the tabular
data in graphical form. The algebraic least squares method performs worst while some
of the renormalisation schemes and the fundamental numerical scheme perform best.
method is systematically deficient, generating average errors up to 22%
greater than the best methods. The SORI and SORII schemes are similarly deficient;
however, they are best seen as incremental developments leading to SORIII. Finally,
the FORI, FORII, FORIII and SORIII schemes are seen to trail FNS only very slightly.
NOISE LEVEL ALS SMP FOR I FOR II FOR III SOR I SOR II SOR III FNS
1:0 2:710 1:093 1:072 1:072 1:071 1:093 1:093 1:071 1:075
2:0 5:579 2:078 1:990 1:987 1:987 2:078 2:078 1:987 1:976
3:0 8:340 3:169 3:077 3:067 3:067 3:169 3:169 3:067 3:049
5:0 15:091 5:662 5:129 5:092 5:092 5:655 5:661 5:092 5:054
8:0 26:036 9:294 8:115 8:037 8:037 9:254 9:288 8:037 7:966
9:0 31:906 10:791 9:036 8:948 8:948 10:748 10:776 8:950 8:827

Table

1: Error results obtained for all methods
Acerage noise level (pixels)515
Average
error
(pixels)
ALS

Figure

results against average noise level depicted graphically. ALS refers
to the algebraic least squares method (with some errors out of range). Group 1 comprises
tabulated results for details.
9 Conclusion
The statistical approach to parameter estimation problems of Kenichi Kanatani occupies
an important place within the computer vision literature. However, a critical component
of this work, the so-called renormalisation method, concerned with minimising
particular cost functions, has proven difficult for the vision community to absorb. Our
major aim in this paper has been to clarify a number of issues relating to this renormalisation
method.
For a relatively general problem form, encompassing many vision problems, we
first derived a practical cost function for which claims of optimality may be advanced.
We then showed that a Sampson-like method of minimisation generates estimates which
are statistically biased. Renormalisation was rationalised as an approach to undoing
this bias, and we generated several novel variations on the theme.
Pivotal in the establishing of a framework for comparing selected iterative minimisation
schemes was the devising of what we called the fundamental numerical scheme.
It emerges that this scheme is not only considerably simpler to derive and implement
than its renormalisation-based counterparts, but it also exhibits marginally superior
performance.

Acknowledgements

The authors are grateful for the insightful comments of Marino Ivancic, Kenichi Kanatani,
Garry Newsam, Naoya Ohta and Robyn Owens. In addition, the authors would like to
thank two anonymous referees for providing suggestions that led to improvements in
the presentation of the paper. This work was in part funded by the Australian Research
Council and the Cooperative Research Centre for Sensor Signal and Information Processing



--R

Statistical Analysis of Measurement Error Models and Applications (Arcata
Fitting conic sections to scattered data
Determining the egomotion of an uncalibrated camera from instantaneous optical flow

Image and Vision Computing 10

Direct least square fitting of ellipses
A buyer's guide to conic fitting
Measurement

Statistical bias of conic fitting and renormalisation

Heteroscedastic regression in computer vision: problems with bilinear constraint
The role of total least squares in motion analysis
Fitting ellipses and predicting confidence envelopes using a bias corrected kalman filter
A note on the least square fitting of ellipses
Nonparametric segmentation of curves into various representations
Fitting conics sections to 'very scattered' data: An iterative refinement of the Bookstein algorithm
Estimation of planar curves
A new approach to geometric fitting


a tutorial with application to conic fitting

--TR
Measurement error models
Fitting ellipses and predicting confidence envelopes using a bias corrected Kalman filter
Estimation of Planar Curves, Surfaces, and Nonplanar Space Curves Defined by Implicit Equations with Applications to Edge and Range Image Segmentation
Ellipse detection and matching with uncertainty
Three-dimensional computer vision
A note on the least squares fitting of ellipses
A buyer''s guide to conic fitting
The Development and Comparison of Robust Methods for Estimating the Fundamental Matrix
On the Optimization Criteria Used in Two-View Motion Analysis
Direct Least Square Fitting of Ellipses
Heteroscedastic Regression in Computer Vision
Statistical Optimization for Geometric Computation
Statistical Bias of Conic Fitting and Renormalization
Nonparametric Segmentation of Curves into Various Representations
The Role of Total Least Squares in Motion Analysis
Optimal Estimation of Matching Constraints
Motion analysis with a camera with unknown, and possibly varying intrinsic parameters

--CTR
Wojciech Chojnacki , Michael J. Brooks, On the Consistency of the Normalized Eight-Point Algorithm, Journal of Mathematical Imaging and Vision, v.28 n.1, p.19-27, May       2007
N. Chernov, On the Convergence of Fitting Algorithms in Computer Vision, Journal of Mathematical Imaging and Vision, v.27 n.3, p.231-239, April     2007
N. Chernov , C. Lesort , N. Simnyi, On the complexity of curve fitting algorithms, Journal of Complexity, v.20 n.4, p.484-492, August 2004
Wojciech Chojnacki , Michael J. Brooks , Anton van den Hengel , Darren Gawley, On the Fitting of Surfaces to Data with Covariances, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.11, p.1294-1303, November 2000
N. Chernov , C. Lesort, Least Squares Fitting of Circles, Journal of Mathematical Imaging and Vision, v.23 n.3, p.239-252, November  2005
Wojciech Chojnacki , Michael J. Brooks , Anton van den Hengel , Darren Gawley, Revisiting Hartley's Normalized Eight-Point Algorithm, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.9, p.1172-1177, September
Wojciech Chojnacki , Michael J. Brooks , Anton van den Hengel , Darren Gawley, From FNS to HEIV: A Link between Two Vision Parameter Estimation Methods, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.26 n.2, p.264-268, January 2004
maximum likelihood;covariance matrix;surface fitting;conic fitting;renormalisation;statistical methods;fundamental matrix estimation