Alignment Using Distributions of Local Geometric Properties.
AbstractWe describe a framework for aligning images without needing to establish explicit feature correspondences. We assume that the geometry between the two images can be adequately described by an affine transformation and develop a framework that uses the statistical distribution of geometric properties of image contours to estimate the relevant transformation parameters. The estimates obtained using the proposed method are robust to illumination conditions, sensor characteristics, etc., since image contours are relatively invariant to these changes. Moreover, the distributional nature of our method alleviates some of the common problems due to contour fragmentation, occlusion, clutter, etc. We provide empirical evidence of the accuracy and robustness of our algorithm. Finally, we demonstrate our method on both real and synthetic images, including multisensor image pairs.
Introduction
Image alignment (variously known as registration, positioning etc.) amounts to establishing a common
frame of reference for a set of images. Alignment is an essential step for many tasks such
as data fusion, change detection, pose recovery etc. and has been widely investigated in various
contexts. A good survey of the existing literature is [4]. Traditionally, the transformation required
to achieve image alignment is computed using either feature matching [13], or a search strategy
that optimises some meaningful image similarity measure (eg. mutual information [19], normalised
cross-correlation [10]). While feature matching methods can give very accurate solutions, obtaining
correct matches of features is a hard problem especially in the case of images acquired using different
sensors or widely different viewing positions. Therefore most methods that use feature matching
for such scenarios either use some specific domain knowledge or assume that the features are well
preserved under the two different views. On the other hand, while the methods using similarity
maximisation are relatively insensitive to variations in sensor characteristics, they can be computationally
expensive and typically need good initial guesses to ensure correct convergence. Therefore
these methods may converge to local minima in the case of images with large transformations between
them.
Recently "direct methods" [2, 11, 17] have had a lot of success. These methods assume that the
images that are to be aligned have a small transformation between them. A multi-scale approach
is used and a measure of image intensity similarity is minimised using optimisation routines. The
multi-scale approach ensures convergence of the optimisation method and the transformation model
is gradually made more complex (from simple translation, to Euclidean to affine) so as to best
describe the transformation between images. However, such methods have the limitation that they
cannot be guaranteed to converge in the cases where there is a large transformation between the
two images (eg. if the scale change is large). Also if we have images that have been taken from very
different view-points and with different radiometric properties (eg. a pair of visible and infra-red
images, or images taken under different lighting conditions), the direct method is not applicable
(In [10], the case of multisensor images is handled by a nonlinear filtering of the images which
emphasises the high frequency components. This enhances the discontinuities in images. However
such an aligment scheme can only be applied in a limited multisensor context). Other methods
that are correspondenceless are restrictive, they either solve for only translation [1], assume that
the image consists of a single, unfragmented curve [12, 15], assume that the images have textured
patterns [16, 3, 21] or cannot handle noise and/or occlusion in a robust fashion [6]. In this paper
we propose an alignment method that overcomes some of these limitations. Our interest is in
developing a method that is robust to the above mentioned variations and can handle large changes
in viewing positions as well as small ones.
(a) Contour C (b) Contour ~

Figure

1: The contour ~
C is a rotated version of contour C. The tangent lines at corresponding
points on the two curves are indicated.
While it is difficult to directly relate the intensity information across images, the geometric properties
of image discontinuities remain relatively stable, ie. the boundaries between regions do not
change under illumination changes etc. The geometric information contained in these stable image
contours is often sufficient to determine the transformation between images. This utilisation of the
geometric properties of image contours to estimate the relevant transformation between images is
the key idea developed in this paper. We assume that the scene is approximately planar and that
the relative geometry between images can be adequately captured by an affine transformation. Fur-
thermore, we assume that the images contain atleast one contour that can be extracted by standard
techniques. For our treatment of the alignment problem for discrete features like points and lines,
see [18].
To develop an intuitive understanding of our method, we will use a simple illustrative example.
Consider the contours shown in Fig. 1. The contour in (b) is a rotated version of the one in (a). Let
the rotation angle between the two contours be '. Aligning the contours is equivalent to estimating
the rotation angle.
Now consider a point p on contour C and its corresponding point on the second contour, ~
cated by the circles). The slope angles of the tangents to points p and ~
are denoted by \Psi and ~
respectively. We can observe that ~
'. Now this relationship holds for every corresponding
point pair on the contours, and hence the rotation angle can be easily recovered from such measurements
on given point pairs. However we do not need to establish point correspondences to extract
(a) Slope distribution of C
(b) Slope distribution of ~
Rotation angles
(c) Rotation estimate

Figure

2: (a) and (b) show the distributions of slope angles for the two contours. (c) shows the
cross-correlation between the distributions. The peak occurs at the correct rotation angle.
the rotation angle. Instead, we note that since the above relation holds for every point pair, it also
holds for the distribution of the measures on each image. In other words if we denote by P (:) the
distribution of slope angles of a contour, we have the simple relationship,
and hence given the two distributions of the slope angles in the two images, we can recover the rotation
angle ' in a simple fashion (This can be done by minimising the norm jjP ( ~
which is equivalent to maximising the cross-correlation between the two distributions. See Fig. 2.
Note that the probability distributions can be computed independently on each image without any
need to establish feature correspondences.
The above idea can be extended to other transformation parameters. Any transformation model
T that consists of n parameters can be reparametrised as
are
independent parameters. This parametrisation is chosen in such a way that each parameter g i is
computable from the observed images using the method we develop in this section. For example,
under a Euclidean transformation, the parameters used would be the rotation angle ' and translations
in x and y directions. Each of these new parameters can be estimated in a manner similar to
that described above for the case of the rotation angle. These notions are formally developed below
and the different transformation models are described in Section 3. In the discussions that follow,
we consider a point p on curve C in image I 1 and its corresponding point ~
on curve ~
C in image
I 2 . The local geometric properties measured at points p and ~
are called geometric descriptors and
are denoted by D and ~
D respectively. D and ~
are chosen in a manner such that for a pair of
corresponding points (p; ~
p),
where g is a parameter of the transformation T . In other words, the pair of operators (D; ~
D) are
chosen in such a manner that for a given corresponding pair of points, the value of ~
D can be related
to that of D given the value of the parameter g.
Intuitively, the above notation implies that for the correct value of the parameter g, the functions
D(:) and ~
D(:) have identical values when computed on corresponding points p and ~
p. This notion
is key to the estimation process since it enables us to choose the parameter that satisfies this
relationship. For an n parameter transformation T , we would need to establish n such relationships
using the different transformation parameters g i 's to estimate the transformation T . We denote
a parameter g i to be "observable" if its value can be recovered from the descriptor pair. In the
case of rotation,we can choose
\Psi. Hence ~
which implies that
is observable, since its value can be obtained from those of D and ~
D. Thus for an n parameter
transformation model, we can recover the transformation given n independent descriptors that make
the parameters observable.
However, since we do not have explicit correspondences available, we can convert observability
of parameters through descriptor pairs to observability through distributions of these descriptors.
Given a descriptor D, we can easily determine its probability density function P(D) in an image by
computing its value on points along the contours and computing the histograms of the computed
values. We denote the probability functions thus obtained from images I and ~
I as P(D) and P( ~
D)
respectively. Now for the correct value of g i , we have the relationship
Hence for
observation of the distribution under noise, we can estimate g i by maximising the similarity between
P(D) and P( ~
In practice, the observed distributions would not be identical to P( ~
to errors introduced by image noise, discretisation of the contours etc. Since the true distributions
of these errors are not analytically tractable, we will assume that the errors are well described by a
monotonically decreasing noise process with mode at zero, like Gaussian, Laplacian etc.
We can now describe the estimator for g i given the pdf's (ie. P(D) and P( ~
D)) as follows:
Theorem 1 Given an image pair I and ~
I, and descriptors D(p) and ~
D(~ p; g), the Maximum Likelihood
Estimator (MLE) for g is
arg min g
P( ~
where P and -
P are the observed probability measures of D and ~
D respectively and
is the L n
norm for n ? 0. The probability distribution function (pdf) of the observation noise is assumed to
be monotonically decreasing with a mode at zero.
For convenience, we shall henceforth denote g(T ) by g.
We denote the true probability distributions of D and ~
D as P and ~
The observed distribution of D is P and that of ~
D is modelled as :
P( ~
Now we know that for every possible distribution pair f P , ~
P g the relationship
P( ~
holds. For notational convenience, we denote the conditional probability measures ~
P( ~
Djg) and
P( ~
Djg) as ~
respectively.
Therefore, using the observation model for estimating
the MLE is defined as
where IP denotes the probability of the observations -
P. 2 In other words, we maximise the probability
of jointly observing the two probability distributions P and -
P on the images I and ~
I respectively.
For the given observation model, we get :
Due to the monotonicity of the noise model we observe that:
By the triangle inequality and Eqn. (2), we know that:
Therefore,
Pg ;Pg
Since the right hand side of equation (4) is equal to k ~
by the assumption of
independence of noise, we have
1 Note that here the true distributions P and ~
P are not random variables since for a given image pair, the
distribution of the geometric properties is deterministic, hence the distributions should be treated as general functions.
P are the observations we extract from the images.
arg min g
Hence the Maximum Likelihood Estimator for the parameter g is given by Eqn. (1).
In effect, the above theorem states that we can determine the optimal estimate of a transformation
parameter by maximising the similarity of the two distributions of the descriptor sets. So in the
case of a Gaussian noise model, our estimation process amounts to picking that parameter which
minimises the least squares error between the two distributions on the two images.
3 Transformation models
In the following subsections, we examine specific transformation models and show how we can
estimate the relevant parameters. We shall denote the points in the first image by
those in the second image by ~
. The transformation model adopted is
~
denotes the 2-D translation vector and the matrix T denotes a 2 \Theta 2 invertible
matrix. Henceforth we shall refer to matrix T as the transformation matrix.
3.1 Euclidean
The Euclidean transformation is parametrised by three parameters, the rotation angle ' and the
two translation parameters t x and t y . Here we have
~
To compute rotation we choose the descriptors D(p) and ~
D(~ p) to be the slope angles of points
on the curves in the two images. From the example shown in Fig. 1, the rotation value can be
computed using the MLE defined above. Having compensated for the rotation between the images,
we can compute the x-direction translation t x between the two images using
~
is the x-coordinate of point p. The y-component of translation can also
be computed in a similar fashion. Note that this estimation process is different from the traditional
method of taking the difference of the centers of mass of the two sets of contours. This is so because
our method of comparing the distributions is a robust estimator unlike the mean (See Section 4 for
an experimental comparision).
3.2 Similarity
To compute the similarity transformation, we need to compute an additional scale parameter s.
Since the radius of curvature (R) of a point is directly proportional to the scaling parameter, we
have p). From this we can deduce that
Hence
we have the simple additive relationship
Also since the radius of curvature is independent of the slope angle , we can compute the scaling
and rotation parameters independently.
3.3 Quasi-affine
We now consider the case of a "quasi-affine" transformation which is defined to be of the form
~
In this case, curvature and slope angles are no longer independent. However we can reparametrise
Eqn. (8) as
~
aeA@ cos ' sin '
sy
To compute the determinant of the transformation (jT we consider the following
The curves C and ~
C are parametrised by the indices s and ~ s respectively. The 2 \Theta 2 matrices P
and ~
defined as
p(s)] and ~
~
p), denote the derivatives
with respect to the appropriate parametrisation s or ~ s of the curves.
We have
x
~
x
From the above we can note that
~
ds
d~s
(See [5] for details). If the parametrisation is chosen such that ds
d~s
we have the simple relationship

~
Therefore we use
p]j and ~
~
to compute the determinant of the
transformation T. By scaling the curves we get the new relationship
~
aeA@ cos ' sin '
sin ' cos 'A -
The rotation angle ' can be computed in a manner invariant to the ratio of the scales ae 2 . To do
this we use the relationship
~
y
The two sides of Eqn. (11) can be equated to D(p) and ~
can be solved for.
The parameter ae, now satisfies the relationship
y
y
y
which implies that
x
y
y
y
Thus we can also solve for the parameter ae and hence we can recover the transformation matrix T
3 . Translation can be recovered as before.
3.4 Affine
Here we consider the more general case of an affine transformation between two images, described
by the equation
~
where
T =@ a b
c dA (15)
is a non-singular matrix. In this case, we need to compute four independent parameters of the
transformation matrix T , which can be accomplished in the following manner.
The determinant jT j can be computed in the manner described in the previous subsection. Now we
use the simple relationship
x
~
x
a -
y
3 Note that in Eqn. 12, ' is known since it is estimated using Eqn. 11.
By a simple reparametrisation of the parameters (a;
in the above
equation, we get
x
~
x
y
from which we can solve for the ratio ( a
using the left-hand side of Eqn. (16) for D(p) and the
right-hand side for ~
By similar analysis, we can recover the ratio ( c
d
using the relationship for - ~
y
y
. Finally, since we now
know the ratios ( a
d
), we can observe that
y
d
a
y
c
d
y
from which the ratio ( b
d
can be computed by taking the logarithm on both sides of the above
equation.
Now we have 4 unknowns and 4 different relationships between these parameters. We denote a
c
d
d
which when substituted into the relationship jT
Using this relationship we can solve for the value of a and then by using the other relationships, we
can derive the value of the matrix T . The sign ambiguity can be easily resolved by considering the
correctness of transformations applied to the image contours.
It may be noted that the reparametrisation of (a; b) and (c; d) results in a finite range of possible
estimate values of the new parameters (ie, range of OE is [\Gamma-]). It may also be noted that using
the above equations we could compute either of the ratios ie. a
or b
a
, or c
d
or d
c
(See Section 4 for
details on how this is used for robust estimation). The translational component can be recovered
in a manner identical to those in the previous cases. It is to be emphasised that in all of the above
cases we have chosen D and ~
D such that the computations on the two sides of the equality of a
relationship can be carried out independently on two different images. This choice enables us to
eliminate the need for correspondences.
4 Implementation and Evaluation
In this section we describe how our method is implemented in practice, following which we consider
various issues that arise out of the implementation and practice of our method. We also present an
empirical evaluation of its performance.
4.1 Implementation
In this subsection we detail the practical issues in the implementations of our method. All the results
detailed in Section 5 are obtained by applying the same fully automatic techniques. A Canny-like
edge detector is used to extract curves. Segments of curves that are smaller than 20 pixels are
discarded. Subsequently, the relevant differential properties are estimated along the curves in both
the images and the distributions computed to estimate the transformation parameters.
ffl Derivative Computation : An important issue in any such implementation is the accurate
computation of the differential properties which are in general sensitive to noise. As shown in [20],
using a Gaussian smoothing kernel to fit curves for computing the derivative of a discrete contour
is incorrect since such an estimator is very biased. To overcome this problem, we use a method
for robust computation of derivatives detailed in [14]. This method involves least squares fitting
of a continuous function (using polynomial orthogonal bases) to a neighbourhood centered on the
discrete point of interest. For such computations, closed form solutions exist in the form of convolution
kernels and the derivative can be easily computed by convolving the curve with these kernels.
We have found the estimates of derivatives using this method to be stable even for the cases where
noise was added to the curves.
ffl Discrete Distribution Representations : Since we use a finite number of bins to compute the
distributions of the geometric properties, we need to choose the bin size to be that of the accuracy
in estimation that we are interested in. For example, if we want to compute the angle OE which
represents the estimate in the affine case for values a
or c
d
to within an accuracy of 0:5 ffi , we need to
use at least 720 bins for the appropriate distribution. Also since the distributions are computed by
assigning each descriptor measure to an appropriate bin, the computational complexity is directly
related to the number of points on the curves for which the geometric descriptors are computed. In
fact if we have N points on the curves and we seek a resolution accuracy of \Delta, the computational
complexity is O( N
). However this computational load is substantially reduced by adopting a multi-resolution
technique. For example, in the case of the estimation of angle, we can estimate the angle
at a coarser resolution (say 2 ffi ) and then refine the estimate by computing the correct angle at a
finer resolution around the current estimate. This is important in the case where the parameter is
not separable (eg., a
for the affine case), in which case the estimation of the MLE amounts to a
search through the space of possible solutions. In our non-optimised implementation of the affine
parameter estimation (in MATLAB), the computation time was typically on the order of tens of
seconds on a Sparc Ultra 1. Faster implementations are easily conceivable. It may also be noted
that in the case where the parameter is separable, ie. when it can be expressed as a function of the
geometric descriptors alone (eg. rotation under the Euclidean model or scale under the Similarity
model), the estimation of the parameter is a straight forward maximisation of correlation which can
Two different parametrisations

Figure

3: The points shown on the two curves are the ones that would be the closest correspondences
for different parametrisation of the curves.
be achieved using convolution. Hence the estimation is considerably faster in such cases and the
computational time is of the order of seconds.
With regard to accuracy it may be noted that in the estimation of the fractions a
or c
d
we use the
tangent function which could lead to inaccurate solutions if the estimated angle OE is close to 90 ffi .
However in such cases, it may be noted that if
is close to 90 ffi , the ratio b
a
would be
close to 0 and can be estimated accurately without any instabilities. This is a strategy we adopt in
our implementation, where the angle tan \Gamma1 ( a
is estimated and if required we switch the representation
to the parametrisation that uses the ratio b
a
instead of a
. The same applies for the ratio of
c and d.
ffl Parametrisation of curves : The estimation accuracy is also dependent on the manner in
which the curves are parametrised. As noted in Section 3.3, if the curves are parametrised such
that ds
d~s
sampling the geometric properties at uniform intervals on the curve would be correct
since for every point we sample on a curve, we would ensure that we sample its corresponding point
on its transformed version. However for our method it is not necessary to parametrise the curves in
a manner that they satisfy the above criteria. Since we have a finite number of bins for representing
any distribution, points that have values within a fraction of the bin size would all fall within the
same bin.
To show that the above is true, consider the following. In Fig. 3, we show a curve with two different
parametrisations. We also indicate a series of points on the curves according to two independent
parametrisations, s and ~
s. Let the step size in both the cases be \Delta. Obviously, if s and ~
s satisfy the
relationship ds
d~s
then the sets of points would be in correspondence. And hence any measurement
of the geometric properties would be equivalent upto the parameter that we want to estimate. Now
consider the case when ds
d~s
1. Let the geometric property we measure be denoted F . Then the set
of geometric descriptor values measured in the first instance, (fFg) would be values computed at
points according to the parametrisation s and those in the second case (f ~
be values computed at points according to parametrisation ~ s which would be
measured on a different set of points than those according to parametrisation s. Now every point
in the first set would have a neighbour within distance \Delta in the second set. Therefore if we assume
that dF
ds
is bounded we have,
ds
where dH
F) is the Hausdorff distance [9] between sets F and ~
F. Thus,
lim
\Delta!0
\Delta!0
Therefore if we sample the curves densely enough we would ensure that for every point sampled in
the first image, we would sample a point close enough to its corresponding point in the second image.
This implies that we can effectively capture the statistical nature of the two distributions without
needing the correct parametrisation. We have found this to be true in course of our experimen-
tation. A simple arc-length parametrisation followed by uniform sampling was found to be sufficient.

Figure

4: The images on the left have 25% overlap. The correctly registered images are shown in
the right half of the image.
ffl Iterative Refinement : One of the underlying assumptions for our method is that the scenes
being imaged are the same, ie. the contours in the two images arise from the same area and hence
minimising the metric would result in the correct estimate for the transformation parameter. How-
ever, in the case where the overlapping areas of the image are small (ie. say or less), this
assumption cannot be true and hence the parameter estimates may not be correct. However by
applying the estimated transformation to the images, we can increase the percentage of the overlap
and hence better approximate the assumptions of a common scene. We have found that by extracting
the overlapping areas after applying this initial estimate and recomputing the transformation
using these areas we can get the correct transformation estimate. We illustrate this in the example
Fig. 4. The two images on the left were extracted from a larger image. The images are of size
300 \Theta 300 and the relative translation between the images is 150 pixels in both the x and y directions
(overlap is 25%). The initial estimate for the translation was (151,124) pixels. The estimated
translation was applied to the images and the overlapping areas of the two images were extracted
(The common area is now about 80%). Hence the new images are a better approximation of the
underlying assumption of a common scene. The required transformation is recomputed and the
new estimate is found to be correct, ie (150,150) pixels. The registered images are shown in the
right half of Fig. 4. We have found that in similar cases with little overlap, correct registration can
be achieved by applying one or two iterations of the above method.
In general, each parameter estimate assumes that the previous estimates are correct. Hence the
stagewise errors can accumulate. However at each stage, since the error is bounded, the maximum
possible error is also bounded. In practice, this problem has not been found to be significant.
Moreover, it can be easily corrected for by either increasing the resolution of our estimates at each
stage (by increasing the number of bins) or by recomputing the transformation using the idea of
iterative refinement as described above.
ffl Multiple peaks : In our method the estimation of any parameter amounts to the maximisation
of a given metric. Therefore we need to detect the largest peak of a given function. In the case
of the computation of rotation, we can see that there are multiple peaks (See Fig. 2) of which
the correct one is ideally the most prominent. However this uniqueness of the correct solution
may be violated in many practical scenarios, eg. when there is a lot of clutter, high amounts of
fragmentation or when there is a small amount of image overlap. In such cases, since the underlying
distributions no longer arise from curves that have a one-to-one correspondence across images, we
can get spurious peaks. A particular case is when we have rectangular buildings in aerial images.
Since buildings typically have strong edges that are oriented 90 ffi apart, we would get peaks that are
90 ffi away from the true solution. We tackle such situations by using a simple verification process to
eliminate spurious peaks. In the case of rotation, it is easy to see that any estimate can be verified
for its accuracy in registering the curves (edges) observed. Choosing the wrong peak would result in
totally incorrect alignment results. Thus in cases where there are multiple competing solutions for
a particular parameter, we maintain them as possible solutions and as we progress in our estimation
process, we eliminate the spurious solutions to arrive at a unique solution. In practice, we have
observed that there are at most 3 different estimates that we need to consider before choosing the
correct one.
4.2 Evaluation
In characterising the accuracy and robustness of any algorithm, typically additive white Gaussian
noise assumptions are made on the observed data. However in our method, it is not possible to
propagate such assumptions to arrive at appropriate noise models for the final representations used
in the estimation process. For example, even if we assume that the image noise is Gaussian, it is
extremely difficult and cumbersome to model the error in the resultant distribution of slope angle
or any other geometric property that we measure, since a number of processing steps are involved.
Moreover, a particularly strong assumption that many image alignment methods (eg. moment
based methods) make is that the curves are complete, ie. not fragmented. In practice this is seldom
the case and this fact has been one of the prime motivations for the development of our method
that uses a distributional representation to alleviate the problem of missing data due to curve frag-
mentation. Under such circumstances it is important to determine the robustness of the alignment
methods. However since it is not possible to accurately model the process of curve fragmentation,
we need to take recourse to an empirical evaluation.
To evaluate the accuracy and robustness of the proposed algorithm, we carried out the following
evaluations.
ffl Translation accuracy and comparision with mean estimation
ffl Rotation accuracy and comparision with moment based methods
ffl Robustness of affine estimation to noise
ffl Robustness of affine estimation to fragmentation
To standardise the evaluation, we used the same test pattern. This test pattern is shown in Fig. 5.
To test the accuracy of our translation method, we subjected the test pattern to a translation
of (10; 20) pixels to form the second curve. Subsequently both the test and translated patterns
were randomly fragmented to the extent that 10% (or 25%) of their total lengths was lost and
we estimated the translation using our method. To compare the accuracy of our method we also
computed the translation by the standard method of taking the difference of the centroids of the
two curves. This experiment was repeated a 1000 times and the results are tabulated in Table. 1 in
the form of the means and standard deviations of the two estimators. As can be easily noted, our
method is more accurate and more robust than the simple process of taking the mean.
10% fragmentation 25% fragmentation
Method Mean Std Mean Std
True Value (10.00,20.00) - (10.00,20.00) -
Our Method
Mean based (9.58,21.25) 2.58 (8.37,22.56) 7.86

Table

1: Comparision of translation estimation accuracy. 1000 experiments each. The true translation
is (10,20) and the table shows the mean estimates in position and the standard deviation of
our method and the mean-based method for two different fragmentation levels.
Method Mean error Std Dev of error
Our method 0:17
Moment based 1:24

Table

2: Comparision of rotation estimation accuracy. 1000 experiments were conducted with 25%
fragmentation
To compare the accuracy of rotation estimation, we used an ellipse and tested the recovery of its
rotation. An ellipse (( rotated by an angle chosen from the range [\Gamma-]. Both
the original ellipse and test ellipse were subjected to fragmentation of 25% of their total lengths.
We then recovered the rotation angle using our method. For the sake of comparision, we also used
a moment based method (by solving from rotation angle using the moment equations) to recover
the angle of rotation. This experiment was repeated 1000 times for different rotation angles and
the absolute error was computed for all the cases. Table. 2 shows the results for the two cases. The
mean and the standard deviation of the absolute error is shown. As can be noted, our method of
rotation estimation is more accurate than that of a moment based one, which like in the case of
translation is a non-robust estimator due to the integrative nature of the moment functions. In such
a case, once the moments are computed, we cannot "separate" out the contributions from spurious
curve segments etc. In this context we would like to point out that most traditional image alignment
techniques that use moments do so in a different manner. In such cases, moment invariants are used
for matching features and then algebraic equations are solved to calculate the required alignment
transformation ( [7]). This is quite different from using a moment based approach to compute the
transformation without using explicit matching. The few cases that do use the moments to calculate
the transformation (eg. [15]) make very strong assumptions that the scene has a single contour and
that there is no fragmentation. Such assumptions imply that there is explicit matching and do not
constitute a general correspondenceless technique.

Figure

5: The test pattern used for evaluation of the affine alignment accuracy.
(a)
level
Average
registration
error
(b)

Figure

The figure on the left shows an example of noise contamination. The noise is white
Gaussian with a standard deviation of 1 pixel. The graph on the right shows the average error for
different noise levels. Notice that the alignment is quite reasonable even for extremely high levels
of noise!
To evaluate the estimation accuracy of our affine parameter estimation method, we test for the
following cases.
ffl Estimation accuracy under noise
ffl Estimation accuracy under fragmentation
For the case of estimation under noise, we apply a known affine transformation to the test pattern
shown in Fig. 5. Subsequently we add white Gaussian noise of different standard deviations of 0:25
to every point on the test pattern. One such instance is shown in Fig. 6 a). To evaluate the accuracy
of the estimated transformation, we use the estimate to warp back the transformed contour onto
the original test pattern and we measure the root mean square error (RMSE) obtained between the
(a) (b)

Figure

7: The circle in the figure on the left is an outlier and does not appear in the second image.
The dotted curve in the figure on the right shows the registration achieved in the presence of the
contaminating outlier. As can be noticed the alignment is quite reasonable given the high levels of
contamination due to the outlier. This can be easily refined to get the correct estimate.

Figure

8: The fragmented, transformed version of the test pattern.
test pattern and the estimated aligned pattern. The average RMSE errors for different noise levels
are shown in Fig 6 b). As can be observed, the performance of our estimator degrades gracefully
and gives reliable estimates even under the severe amounts of noise shown in Fig. 6 a).
To study the effects of outliers on our estimation process, we consider the scenario in Fig. 7. The
large circle is not present in the second image, hence it will corrupt the distributions of different
geometric properties. In this case, the effect of the presence of the circle is severe since it contributes
large numbers of samples to the distributions. However as can be observed in Fig. 7, the estimate
of the transformation is reasonable considering the severe amounts of contamination. Given this
estimate, it is easy to reject the outliers and we can reestimate the transformation to get the correct
solution for the affine transformation between the two images.
To evaluate our method under fragmentation of contours, we used the same test pattern shown in
10% fragmentation 25% fragmentation
Measure (in pixels) (in pixels)
Median 0.45 0.97
Mean 0.58 1.11
Std. Dev 0.27 0.72

Table

3: RMSE for 1000 experiments with fragmentation
Fig. 5 and subjected it to the same affine transformation used in the evaluation for robustness to
noise. Then we fragmented both the contours and estimated the transformation. We considered
fragmentations of 10% and 25% of the total length and repeated the experiments a 1000 times in
each case (One instance of the fragmented pattern is shown in Fig. 8). The results are tabulated in

Table

. 3. As can be noted, our estimation process is fairly robust to even severe fragmentations of
the order of 25% of pixels in both contours.
In this section we present the results obtained by applying our method to a variety of images. The
same technique is used for all the examples. The implementation is detailed in Section 4.

Figure

Figure
An image and its affine transformed version

Figure

9: An image and its affine transformed version.
Fig. 9 shows a synthetic image and its affine transformed version. The resulting registration achieved
is shown in Fig. 10 and has sub-pixel accuracy (The root mean squared error (RMSE) is about 0.5
pixels). In Fig. 11 we show the registration achieved when about 25% of the image is occluded. As
can be observed from the result, the registration accuracy is not affected since the nature of the
Registered images

Figure

10: Alignment result of the images in Fig. 9. The transformed contours of the second image
are shown as dotted curves.
distributions are not altered to the extent that the MLE's are significantly perturbed.
Fig. 12 shows a mosaic constructed by aligning images from a video sequence in a common frame of
reference using a Euclidean model. The sequence was obtained by a hand-held camera and involved
translation and rotation. In this case a Euclidean model was sufficient for achieving accurate
alignment.
Fig. 13 shows the alignment of aerial images of the Mojave desert obtained from a balloon flying over
the area. The alignment achieved is accurate as is evident from the alignment of the image features
like the roads, rock outcrops etc. In Fig. 14 we show the results for aligning a pair of images from
the Landsat Thematic Mapper (TM). The images are from different bands of the electro-magnetic
spectrum and therefore have different radiometric properties. The images are correctly aligned as
can be observed from the continuity of the coastline and the alignment of features that run across
the two images. Fig. 15 shows the results obtained for another pair of Landsat TM images. We
used the quasi-affine transformation model to achieve the alignment for Figs. 13, 14, 15.
As an illustration of our algorithm for affine transformation estimation, we show two examples
from different application domains. In 16 we show the registration achieved for two MRI images of
different modalities. The image on the left in Fig. 16 is a proton density MRI image and the one in
the middle is the corresponding T2 weighted image that has been subjected to an arbitrary affine
transformation. The image on the right shows the alignment of a single contour from the different
modalities. It may be noted that the correct alignment is achieved inspite of the difference in the
photometric properties of the images. The alignment can be further refined if necessary using any
of the standard energy minimisation techniques.
Original Image
Occluded Transformed Image
Registered image

Figure

11: Alignment of occluded image. Only the transformed version of the occluded image is
shown. The registration achieved is to within subpixel accuracy.
Finally, we demonstrate the applicability of our method to alignment of multi-sensor images that
have large transformations between them. In Fig. 17, we show two satellite images of a river taken
under the SPOT and TM modalities. As can be observed, the image intensity patterns are quite
different and therefore any intensity based "direct method" is not directly applicable. Moreover,
the relative scaling between the images is large (about 1:4). Such large scalings are seldom dealt
with using optimisation schemes. Our alignment of the two images is shown in Fig. 18. We use a
checkerboard pattern to illustrate the alignment of the two images from different sensors. Minor
discrepancies do exist in the alignment (probably due to the deviation of the images from the affine
model), but the overall alignment is very good especially given the large scaling between the two
images. If desired, the final alignment can be easily refined using different techniques.
All the above examples illustrate the ability of our algorithm in achieving the correct alignment for a
variety of possible scenarios. We have also demonstrated that our method gives good performance
for a broad range of cases especially in the case where the transformations between images are
large and when there is a small amount of overlap between the images, something that energy
minimisation methods would fail to achieve.
6 Discussions
An advantage of our method is the fact that we use the same method for computing the alignment
for all cases. This is possible since we use the same distributional framework of local geometric
(a) Some frames (b) Mosaic

Figure

12: Mosaic created out of a video sequence. Some of the frames used are shown in (a).
properties to estimate the transformation parameters. As a result, the method works just as well
for images with a large transformation between them as with images with small transformations.
This is not always the case with energy minimisation methods which typically require the solution
to be started fairly close to the true solution to avoid being trapped in local minima. Also any
explicit need for domain dependent knowledge is avoided since the image primitives we use are
contours that are usually extractable using standard edge-detection techniques. We believe that
this is a significant advantage in the case of multi-sensor images since most methods that deal with
multi-sensor image alignment utilise specific knowledge about the imaging geometry and radiometry
to tune their algorithms to work on the specific domains of application.
The use of a distributional method has other significant advantages too. Fragmentation of contours
is easily handled since the local geometric properties can still be computed on the fragments without
any significant loss of information. This is possible since if we have a contour, say, of length 200
pixels, fragmented into many segments that add up to say 180 pixels in total length, then we still
have 180 points on which we can compute the geometric properties. The probability distribution
that we estimate is not significantly changed. Therefore, we can observe a graceful degradation of
the method with increasing loss of information due to occlusion or fragmentation. As demonstrated
in the evaluations, such robustness would not be possible with the use of traditional global techniques
like moments unless special care is taken to handle the change in shape due to occlusion or
fragmentation. It is easy to note that the results of a moment-based method would be perturbed
significantly due to small amounts of fragmentation of contours, while our method would still give
the correct estimates.
In most voting based schemes, there is a combinatorial explosion of possible solutions that need
to be checked as we increase the dimensionality of the search space or more significantly, as the

Figure

13: Alignment of a set of aerial images
number of observations are increased. This is so since most voting schemes, use all possible combinations
of "hypothesised" feature matches to populate the space of possible solutions and then
search for a maxima in this space. In our case, we have two advantages in this regard. Firstly, all
computations are carried out independently on each image and its only in the final stage that the
two distributions are compared. Secondly, each stage of parameter estimation is one-dimensional
since we parametrise each transformation into a set of independently estimated parameters. As a
result, the computational load is reduced. Also detection of peaks in one-dimensional functions is
easier than in higher dimensions especially under noisy scenarios.
As demonstrated by our examples, the theoretical framework we have developed works well for
many real-life examples. While the existence of clutter, occlusion or the fragmentation of image
contours does violate the underlying assumptions, our method is robust enough to be able to easily
handle these problems. However it would be advantageous to modify the metric that we minimise
so as to take into account knowledge about the scene. If we can determine which parts of the image
contain areas that are visible in both images, which parts are occluded etc. then we should modify
the metric to exploit this information to make our method more robust. One limitation of our
framework is that the scene is assumed to be approximately planar which may not be the case for
scenarios where the perspective effects in the images are dominant. However we can still use our
method to get a reasonable estimate of the true transformation.

Figure

14: Alignment of LANDSAT TM images from different bands.

Figure

15: Alignment of two LANDSAT TM images

Figure

The image on the left is a proton density image, the one in
the middle is T2 weighted and the image on the right shows the registration of a single contour of
the two modalities
7 Conclusion
We have described a framework for image alignment that does not use explicit feature correspon-
dences. We have demonstrated the effectiveness and explained the advantages of using a distributional
framework for computation of the parameters required for image alignment. This framework
is robust and is more general than many existing methods.

Acknowledgments

We would like to acknowledge E. Rignot , B. S. Manjunath and maintainers of the UCSB Image
Registration web site for providing some of the data sets used in this paper. We would also like to
thank R. Chellappa, Z. Duric and J. Oliensis for comments on the paper.



--R

"A Robust, Correspondenceless, Translation-Determining Algo- rithm"
"Hierarchical Model-Based Motion Estimation"
"Shape from Texture : Estimation, Isotropy and Moments"
"A Survey of Image Registration Techniques"
"Invariant Signatures for Planar Shape Recognition Under Partial Occlusion"
"Recovering 3D Rigid Motion Without Correspondence"
"A Moment-Based Approach to Registration of Images with Affine Geometric Distortion"
"Image Registration Without Explicit Point Correspondences"
Comparing Images Using the Hausdorff Distance"
"Robust Multi-Sensor Image Alignment"
"Recovery of Ego-motion Using Image Stabilisation"
"Recovery of Global Nonrigid Motion - A Model Based Approach without Point Correspondences"
"A Contour-Based Approach to Multisensor Image Regis- tration"
"Smooth Differentiation Filters For Images"
"Using moments to acquire the motion parameters of a deformable object without correspondences"
"Extracting the affine transformation from textured moments"
"Compact Representations of Videos through Dominant and Multiple Motion Estimation"
"Multisensor Image Registration Using Feature Consensus"
"Alignment by Maximization of Mutual Information"
"Noise Resistant Invariants of Curves"
"Recovering surface shape and orientation from texture"
--TR

--CTR
Rujirutana Srikanchana , Jianhua Xuan , Matthew T. Freedman , Charles C. Nguyen , Yue Wang, Non-Rigid Image Registration by Neural Computation, Journal of VLSI Signal Processing Systems, v.37 n.2-3, p.237-246, June-July 2004
Kenneth Nilsson , Josef Bigun, Localization of corresponding points in fingerprints by complex filtering, Pattern Recognition Letters, v.24 n.13, p.2135-2144, September
Christophe Doignon , Dominique Knittel, Detection of noncircularity and eccentricity of a rolling winder by artificial vision, EURASIP Journal on Applied Signal Processing, v.2002 n.1, p.714-727, January 2002
Bogdan Georgescu , Peter Meer, Point Matching under Large Image Deformations and Illumination Changes, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.26 n.6, p.674-688, June 2004
geometric properties;statistical distributions;pose estimation;correspondenceless image alignment