Pseudo-polar based estimation of large translations rotations and
scalings in images
¤Y.Keller1;2,A.Averbuch1,M.Isreali3
1School of Computer Science
Tel Aviv University,Tel Aviv69978,Israel
2Dept.of Electrical Engineering Systems
Tel-Aviv University Tel-Aviv69978,Israel
3Faculty of Computer Science,Technion,
Haifa32000,Israel
Phone:*972-55-694-464
Fax:*972-151-55-694-464
Email:keller@post.tau.ac.il
Submitted to the IEEE Transactions on Image Processing
August2002
Abstract
One of the major challenges related to image registration is the estimation of large motions without prior knowledge.This paper presents a Fourier based approach that estimates large
translation,scale and rotation motions.The algorithm uses the pseudo-polar transform to
achieve substantial improved approximations of the polar and log-polar Fourier transforms of
an image.Thus,rotation and scale changes are reduced to translations which are estimated using
phase correlation.By utilizing the pseudo-polar grid we increase the performance(accuracy,
speed,robustness)of the registration algorithms.Scales up to4and arbitrary rotation angles
can be robustly recovered,compared to a maximum scaling of2recovered by the current state-
of-the-art algorithms.The algorithm utilizes only1D FFT calculations whose overall complexity
is signi…cantly lower than prior works.Experimental results demonstrate the applicability of
these algorithms.
Keywords:Global motion estimation,Sub-pixel registration,Gradient methods,image
alignment
EDICS Category={2-ANAL,2-MOTD}
1Introduction
Image registration plays a vital role in many image processing applications such as video compres-sion[1,2],video enhancement[3]and scene representation[4,5,6]to name a few.This problem was analyzed using various computational techniques,such as pixel domain Gradient methods[1,4,5], correlation techniques[20]and discrete Fourier(DFT)domain algorithms[10,16].Gradient meth-ods based image registration algorithms are considered to be the state-of-the-art.But,may fail unless the two images are misaligned by a moderate motion.Fourier based schemes,which are able to estimate relatively large rotation,scaling and translation,are often used as bootstrap for the more accurate gradient methods.The basic notion related to Fourier based schemes is the shift property[23]of the Fourier transform which allows robust estimation of translations using the nor-malized phase-correlation algorithm[10,12,13,14,15].Hence,in order to account for rotations and scaling,the image is transformed into a polar or log-polar Fourier grid(referred to as the Fourier-Mellin transform).Rotations and scaling are reduced to translations in these representations and are estimated using phase-correlation.
In this paper we propose to estimate the polar and log-polar DFT using the pseudo-polar FFT (PPFFT)[27,28].The PPFFT estimates the DFT on a non-Cartesian grid,which is geometri-cally similar to the polar grid.Hence,by interpolating the polar and log-polar DFT using the PPFFT a signi…cantly smaller interpolation error is achieved.Furthermore,the smaller the rela-tive rotation/scaling between the images becomes,the smaller the di¤erence between their PPFFT magnitudes and the interpolated polar/log-polar DFT magnitudes.In prior works[16,21]the relative interpolation error is not decreased by the decrease of the relative motion.
The resulting algorithm is able to robustly register images rotated by arbitrary angles and scaled up to a factor of4.It should be noted that the maximum scale factor recovered in[16]and[21] was2.0and1.8,respectively.In particular,the proposed algorithm does not result to interpolation in either spatial or Fourier domain.Only1D FFT operations are used,making it much faster and especially suited for real-time applications.
The rest of paper is organized as follows:Prior results related to FFT based image registration are given in Section2,while the pseudo-polar FFT,which is used as a basis for the proposed algo-rithm is presented in Section3.The proposed algorithm,is presented in Section4and its accuracy is analyzed in Section5.Experimental results are discussed in Section6and…nal conclusions are given in Section7.
2
Previous related work 2.1Translation estimation
The basis of the Fourier based motion estimation is the shift property [23]of the Fourier transform.Denote by
F f f (x;y )g ,b f (!x ;!y )(2.1)
the Fourier transform of f (x;y ).Then,
F f f (x +￠x;y +￠y )g =b f (!x ;!y )e j (!x ￠x +!y ￠y ):(2.2)Equation (2.2)can be used for the estimation of image translation [10,15].Assume the images I 1(x;y )and I 2(x;y )satisfy
I 1(x +￠x;y +￠y )=I 2(x;y ):
(2.3)Equation (2.3)is Fourier transformed,yielding
b I 1(!x ;!y )e j (!x ￠x +!y ￠y )=b I 2(!x ;!y )(2.4)
and
b I 2(!x ;!y )b I 1(!x ;!y )=e j (!x ￠x +!y ￠y ):(2.5)Thus,the translation parameters (￠x;￠y )can be estimated in the spatial domain [10,11]by taking the inverse FFT of Eq.(2.5)
Corr (x;y ),F ?1n e j (!x ￠x +!y ￠y )o =±(x ?￠x;y ?￠y )(2.6)
and by …nding the position of the maximum value of the correlation function Corr (x;y ):
(x;y )=arg ?max (~x ;~y )f Corr (~x ;~y )g ?:(2.7)
By normalizing the input images’magnitude,the Normalized phase correlation [10,11]is derived
]Corr (!x ;!y ),b I 1(!x ;!y )b I ¤2(!x ;!y )ˉˉˉb I 1(!x ;!y )ˉˉˉˉˉˉb I ¤2(!x ;!y )ˉˉˉ=e j (!x ￠x +!y ￠y )(2.8)
where *denotes the complex conjugate.
This estimator of ]Corr
is robust to gain and o¤set changes I 1=a ￠I 2+b
where a and b are scalars.The intensity gain a is compensated by Eq.(2.8)while the o¤set b
a¤ects only the DC coe￠cient of ]Corr
.
This scheme was proven to robustly estimate large translations where the corresponding overlap between the images to be registered is down-to 30%of the smallest image size [15].As no smoothness assumption is used,non-smooth and noisy functions (such as 2D DFT coe￠cients)can be accurately registered.Shekarforoush et-al [11]extended the phase-correlation based algorithm to subpixel
accuracy by analyzing the shape of ]Corr
(!x ;!y )in Eq.(2.8)around its maximum.A di¤erent approach to phase correlation based translation estimation [9,14,15,30]utilizes 2D linear regression in order to …t the phase values calculated in Eq.(2.8)to a two-dimensional linear function ?i ￠log 3]Corr (!x ;!y )′=!x ￠x +!y ￠y 8(!x ;!y ):(2.9)Solving Eq.(2.9)using linear regression might prove inaccurate [11]due to aliasing and phase wrapping (around 2?)of the spectra at low-frequencies.An iterative solution to phase unwrapping was suggested in [14]while Stone et-al [15]presented two approaches for modelling aliasing e¤ects in order to improve the registration accuracy.The phase domain regression approach can be naturally extended to subpixel accuracy [30],since Eq.(2.9)holds for non-integral values of (￠x;￠y ).
2.2Polar Fourier representations
The polar Fourier representation (Fourier-Mellin transform)is used to register images that are misaligned due to translation,rotation and scale [16,21].Let I 1(x;y )be a translated,rotated and scaled replica of I 2(x;y )
I 1(x;y )=I 2(s ￠x ￠cosμ0+s ￠y ￠sinμ0+￠x;?s ￠x ￠sinμ0+s ￠y ￠cosμ0+￠y )(2.10)where μ0,s and (￠x;￠y )are the rotation angle,scale factor and translation parameters,respec-tively.The DFT of Eq.(2.10)is
b I 1(r;μ)=e j (!x ￠x +!y ￠y )11s b I 2(r j s j ;μ+μ0):(2.11)
Hence,M 1and M 2,the magnitudes of b I 1and b I 2,respectively,are related by rotation and scaling
around the DC component M 2(r;μ)=12M 1μr ;μ+μ0?:(2.12)
where (r;μ)are polar coordinates !x =r cos μ;
!y =r sin μ:
Therefore,rotations and scale changes can be recovered …rst,regardless the translation https://www.doczj.com/doc/ed7387479.html,ing a polar or log-polar DFT,rotation and scaling are reduced to translations,which can
be recovered using the phase-correlation procedure described in https://www.doczj.com/doc/ed7387479.html,ing Eq.(2.12)to estimate the rotation angleμresults in an ambiguity of?[16]in the estimation of the rotation angle.This ambiguity is resolved by applying both hypothesesμandμ+?and recovering the translational motion(￠x;￠y)and the correlation peak of each hypothesis.The rotation hypoth-esis and translation values,which are related to the highest correlation peak,are chosen as the result.
Thus,the following algorithm is used to account for rotation and scaling:
1.The input images are transformed into a polar or log-polar Fourier grid using bilinear inter-
polation.
2.Rotations and scaling in these representations are reduced to translations and can be esti-
mated using phase-correlation(Eq.(2.8)).
3.The translation is estimated twice-for rotations ofμandμ+?.The result with the highest
correlation peak is chosen.
The estimation of the Fourier transform on a polar or log-polar grid can be conducted using either Image domain warping[20]and then2D FFT calculation[19],or interpolation of the2D DFT transform in the Fourier domain[16,21].
Figure1illustrates the above algorithm for rotation and translation estimation.The interpo-lation used in step#1induces signi…cant errors to the approximation of the polar and log-polar FFT,thus limiting the accuracy of the motion estimation[20].
3The pseudo-polar FFT
The proposed registration algorithm is based on the existence of a fast,accurate and invertible discrete pseudo-polar FFT(PPFFT)[27].An FFT where the evaluation frequencies lie in an oversampled set of non-angularly equispaced points which we call the pseudo-polar(PP)grid.The PPFFT includes fast forward and inverse transforms and a quasi-Parseval relation.A special feature of this approach is that it involves only1-D equispaced FFT’s.In particular,there is no need for re-gridding or interpolation.We describe below how to compute the PPFFT of a function given on a Cartesian grid with all the above features.
Section3.1presents the Fractional FFT which is used in the computation of the PPFFT in section3.2.A geometric interpretation of the PPFFT is given in section3.3.
The description of the inverse transform is omitted,since it is not used by the proposed regis-tration algorithm,a detailed description of it is given in[27].
Figure1:The?ow of the state-of-the-art FFT based image registration:1.The magnitudes of the polar DFT are approximated by interpolating the magnitudes of the2D FFT.2.The rotation ￠μis recovered using1D phase-correlation on theμaxis. 3.One of the input images is rotated by angle￠μ. 4.The translation is recovered and theμambiguity is resolved by applying a2D phase-correlation twice:using I1rotated by bothμandμ+?.
3.1Fractional FFT
The Fractional FFT(FRFT)[29](also referred to as the Chirp Fourier transform[24]),is a fast (N log N)algorithm,for the computation of any set of M equally spaced samples of the DFT on the unit circle
!k=!0+k￠!;k=0;1;:::;M?1(3.1)
where!0,￠!;M are arbitrary and N is the length of the input signal.For the speci…c case of the DFT!0=0,M=N and￠!=2?
:The frequency sampling properties of the FFT and FRFT are
N
Figure3:The construction of the pseudo-polar grid from the Z and N sub-grids.
3.2.1P P Z computation
The computation of P P Z is based on the separability of the2-D DFT[2].Hence,P P Z is computed by applying the1-D FFT in the Y direction and then computing the CFRFT in the X direction.
A varying?factor is chosen such that the Z sub-grid geometry is achieved.
Algorithm Flow
1.Let I be the input image of size(n0;m0).I is zero-padded to a size of(n;n)=?2k;2k￠k2Z,
where k is chosen such that2k?2￠Max(n0;m0):
2.A1-D FFT is applied on the Y direction I b x=F F T X(I)and the result is cyclically shifted, such that the DC component is in the center of I b x.(The shifting operation is similar to the “¤tshift”command in Matlab TM).
3.The CFRFT operator F?is applied to the rows of I b x
P P i N=F?
i?I i b x￠(3.2)
where P P i Z is the i th row of P P Z,F?
i
is the CFRFT operator with a compression ratio of?i and I i
b x is the i th row of I b x and
?i=i·n
2
?=n?2i
n
i=n
2
?=0
i>n
2
?=?n?2i
n
9>>=
>>;:(3.3)
Starting from the outer rows,each row is shrunk by a single sample more than its predecessor: the…rst row has?i=1(no compression),while the second is compressed to an interval of N?1??i=n?1n￠and the third row will be concentrated over an interval of N?1??i=n?2n￠:
4.The lower half of P P Z is ?ipped left-right.
For the construction of P P N we transpose the input image I apply the above
algorithm.
(a
)
(b )
Figure 4:Frequency samples of the FFT and Fractional FFT on the unit circle using N samples and spacing of ￠!.(a )The FFT samples the Fourier transform equally over [0;2?],where ￠!=2?N .(b )The FRFFT samples the Fourier transform equally with any spacing ￠!and initial phase !0.
3.3Geometric interpretation and properties
P P Z and P P N are matrices of size n ￡n whose elements correspond to the sampling of the DFT on the PP grid.By examining Fig.3,it follows that a column in the matrices P P Z and P P N corresponds to a ray on the PP grid with a …xed angle μ,while the rows correspond to a constant distance r .The polar and PP grids di¤er on the non-uniformity of the angular and radial spacings of the PP grid shown in Fig.5.For the polar grid we have
￠μP olar (i )=2?;￠r P olar (j )=￠r 08i while in the PP grid,￠μP P (i )and ￠r P P (j )vary smoothly as a function of i
μP P (i )=arctan μ2i ?;i =0;:::;n (3.4)
where n
(3.5)
(3.6)
Figure5:￠μP P and ￠r P P,
4The
This section FFT (PPFFT)to domains, respectively.in Section2.1.Thus,the phase-correlation algorithm is applied on three domains:spatial,polar DFT’s magnitude and log-polar DFT’s magnitude.Table1summarizes the properties of the translation estimation in each domain.Detailed implementation issues are given in section4.1.
Section4.2introduces the rotation and translation estimation algorithm illustrated in Fig.6, while the simultaneous estimation of scaling,rotation and translation is presented in section4.3 and Fig.7.
4.1Translation estimation implementation
The translation is estimated by Eq.(2.8)where several implementation issues are considered according to the domain on which the motion is conducted.The evaluation of Eq.(2.8)requires the input signals to be zero padded in order to have the same size.Furthermore,implementing Eq.
Table1:Properties of the translation estimation algorithms in various domains:spatial,polar DFT’s magnitude and log-polar DFT’s magnitude.
(2.8)using the DFT estimates the cyclic correlation[24]which di¤ers from the regular correlation needed for translation estimation.The wrap-around problem is solved by zero padding[24].Section 4.1.1presents the padding used for spatial domain translation estimation,while in sections4.1.2 and4.1.3the?!μaxis is cyclic and no padding is needed.
4.1.1Spatial domain translation estimation
Denote by length(I)the length of the1D signal I.Then,in order to avoid the cyclic wrap-around problem in the DFT computation,both signals are zero-padded such that:
length3e I1′=max(length(I1);length(I2))+￠X max
length3e I2′=length3e I1′(4.1) where
2e I1and e I2are the zero padded versions of I1and I2,respectively.
2￠X max is the maximum translation value for which the correlation function is evaluated.
In the2D case the input signals I1(i;j)and I2(i;j)are zero padded in both axes.This is done according to Eq.(4.1)using the maximal translations(￠X max;￠Y max)in both axes.
4.1.2Translation estimation in the polar domain
In the polar domain,the angular axis?!μis cyclic over the interval[0;2?].Hence,no zero padding is needed to avoid the alignment of non-corresponding samples.The images are padded to the same size.
4.1.3Translation estimation in the log-polar domain
In the log-polar domain a distinction is made between the angular axis?!μand the radial axis?!r. The motion along the angular axis?!μis estimated according to section4.1.2,while the motion along the radial axis?!r is estimated according to section4.1.1.
4.2Simultaneous estimation of rotation and translation
Let I1(i;j)and I2(i;j)be the input images to be registered.I(n)1(i;j)denotes the fact that I1 evolves throughout the registration process where initially I(0)1(i;j)=I1(i;j)(n=0).I2(i;j) remains unchanged.The rotation and translation estimation algorithm operates as follows:
1.Let(m1;l1)and(m2;l2)be the sizes of I1(i;j)and I2(i;j);respectively.Then,at iteration
n=0,I1(i;j)and I2(i;j)are zero padded such that
m1=l1=m2=l2=2k;k2Z:(4.2)
2.The magnitudes M P P
1?μi;r j￠and M P P2?μi;r j￠of the PPFFTs of I(n)1(i;j)and I2(i;j)are calculated,respectively.
3.The polar DFTs,magnitudes c M P olar1?μi;r j￠and c M P olar2?μi;r j￠of I(n)1(i;j)and I2(i;j)are
substituted by M P P
1?μi;r j￠and M P P2?μi;r j￠respectively.
4.The translation along the?!μaxis of M P P1?μi;r j￠and M P P2?μi;r j￠is estimated using the
procedure described in section4.1.2.The result is denoted by￠μn.
5.Letμn be the accumulated rotation angle estimated at iteration n
μn,
n
X i=0￠μi=μn?1+￠μn:
Then,the input image I1(i;j)is rotated byμn(around the center of the image)using fast FFT based image rotation[7].
I(n+1)
1
(μ;r)=I(0)1(μ+μn;r);n=1;:::
Rotation around the central pixel of I1(i;j)is recommend,making the bounding rectangular of the rotated image as small as possible.
6.Steps2-5are reiterated until the angular re…nement term￠μn is smaller than a prede…ned threshold"μ,i.e.j￠μn j<"μ,or a prede…ned number of iterations n max is reached.
7.The rotated images I (n )
1(i;j )and I 2(i;j )are used as inputs to the spatial domain translation estimation algorithm presented in section 4.1.1.The rotation ambiguity is resolved by the procedure presented in section 2.2.4.3Simultaneous estimation of rotation,translation and scaling
Let I 1(i;j )and I 2(i;j )be the input images to be registered.I (n )1(i;j )denotes the fact that I 1
evolves throughout the registration process,while I 2(i;j )remains unchanged.We assume that initially I (0)
1(i;j )=I 1(i;j )(n =0).The scaling,rotation and translation estimation algorithm operates as follows:
1.Let (m 1;l 1)and (m 2;l 2)be the sizes of I (n )1(i;j )and I 2(i;j );respectively.Then,at iteration n =0,I (0)1(i;j )and I 2(i;j )are zero padded such that
m 1=l 1=m 2=l 2=2k ;k 2Z :(4.3)The padding is needed to abide by the limitations of the PPFFT described in section 3.
2.The magnitudes M P P 1(μ;r )and M P P 2(μ;r )of the PPFFTs of I (n )
1(i;j )and I 2(i;j )are cal-culated,respectively.
3.Let base length (M P P 2(μi ;r j ))=length ?M P P 2?μi ;r j ￠￠(
4.4)
be the logarithmic base used to de…ne the log domain in the ?!r axis.The log-polar DFTs
of I (n )
1(i;j )and I 2(i;j )are approximated by M P P 1(μ;r )and M P P 2(μ;r );respectively,using the following procedure:
c M Log ?P olar 1(i;j )=M P P 1?i;￥base j |￠c M Log ?P olar 2(i;j )=M P P 2?i;￥base j |￠(4.5)
where b x c denotes the integral part of x (x >0).The polar axis is approximated using the same procedure as in section 4.2,while the radial axis is approximated using nearest-neighbor interpolation.
4.The relative translation between c M Log ?P olar 1(i;j )and c M Log ?P olar 2(i;j )is recovered by a 2D phase correlation on the ?!μand ?!r axes,which was described in section 4.1.3.
Figure6:The?ow of the pseudo-polar FFT based rotation and translation image registration algorithm:1.The magnitude of the pseudo-polar FFT is computed after the input images were zero padded to the same size.2.The rotation￠μis recovered by applying1D phase-correlation on theμaxis in the pseudo-polar domain.3.One of the input images is rotated by the accumulated rotation angleμn.4.Steps1-3are reiterated until a stoppage condition is met.5.The translation is recovered and theμambiguity is resolved by applying the2D phase-correlation twice:using I1 rotated byμandμ+?.
5.Let￠μn and￠r n be the rotation angle and the scaling value estimated at iteration n,re-
spectively.Then,the input image I1(x;y)is rotated(around the center of the image)[7]and then scaled using DFT domain zero padding
I(n+1)
1
(μ;r)=I(0)1(μ+μn;r￠r n)(4.6)
where
μn=
n
X i=0￠μi=μn?1+￠μn
r n=
n
Y i=0￠r i=r n?1￠￠r n:(4.7)
6.Steps1-3are reiterated until either the updated estimates are smaller than prede…ned thresh-
olds j￠μn j<"μand￠r n<"r or a prede…ned number of iterations n max is reached.
7.The rotation and scale parametersμn and r n;respectively,are used as inputs to the spatial
domain registration algorithm presented in section4.2step7.
4.4Complexity estimation
The overall complexity of the proposed algorithm is O?n2log n￠,since it is dominated by the complexity of the2D FFT:
2The pseudo-polar FFT is computed using O?n2log n￠operations(Section3).
2The fast image rotation algorithm[7]uses O?n2log n￠operations.
2O?n2log n￠operations are needed for image zooming using FFT domain zero padding. 5Interpolation accuracy analysis
This section provides performance analysis of the proposed algorithm by comparing it to prior state-of-the-art techniques.The registration accuracy and robustness of FFT based algorithms in the polar and log-polar grids are related to the interpolation accuracy of the DFT on these grids. In prior works[13,16],the DFT’s magnitudes on polar and log-polar grids were interpolated using the Cartesian2D DFT while the proposed algorithm uses the PPFFT.
The interpolation accuracy depends on the distance between the interpolated grid points(polar, log-polar)and the source grid used for the interpolation(pseudo-polar,Cartesian).We show that the accumulated distance between the pseudo-polar and polar/log-polar grids is smaller than
Figure7:The?ow of the pseudo-polar FFT based algorithm for scaling,rotation and translation estimation of image registration:1.The magnitude of the pseudo-polar FFT is computed after the input images were zero padded.2.The rotation￠μand scaling￠r are recovered by applying2D phase-correlation on theμand r axes in the pseudo-polar domain.3.One of the input images is scaled and rotated by an accumulated rotationμn and scale￠r.4.Steps1-3are reiterated until a stoppage criteria is met.5.The translation is recovered and the§?ambiguity in the value ofμis resolved by applying the2D phase-correlation twice:using I1rotated byμandμ+?.
the accumulated distance between the Cartesian and polar/log -polar grids.Thus,improving the registration accuracy and robustness of rotations and scalings.
The angular di¤erences between the PPFFT and uniform polar grids are analyzed in Section
5.1,while radial accuracy analysis (with respect to the polar domain)is given in Section 5.2.Accumulated distance calculations for several image sizes are given in Section 5.3,where it is compared to the accumulated distances of interpolations based on the Cartesian grid.
5.1Angular accuracy analysis
The polar grid shown in Fig.5is characterized by a constant angular spacing
￠μP olar =2?
(5.1)
where n is the size of the polar ?!μaxis.￠μP P ,the PPFFT’s angular spacing is given in Eq.(3.5)
and illustrated in Fig.8.It changes smoothly as a function of i the initial ray index.
Figure 8:Illustration of the non-uniformity of ￠μP P ,the angular spacing of the pseudo-polar grid.
In order to analyze the non-uniformity of ￠μP P ,we expand ￠μP P using a second order Taylor series expansion
arctan (x +h )?arctan (x )=
11+x 2h +122x (1+x 2)
2h 2+O (h )3:(5.2)
By substituting Eq.(3.5)into Eq.(5.2)?x =2i n ;h =2n ￠we get ^￠μP P (h;i ),11+?2i n ￠2
μ2?+2i n 31+?2i n ￠2′2μ2?2(5.3)=
2n n +4i +8i (n 2+4i )2n =2n 22+O μ15
?:
For common image sizes (n ?300?500);O
?15￠
=O ?10?13￠and the error term in Eq.(5.3)is negligible.Then,^￠μP P ?2n
n +4i (5.4)
and a constant angular spacing is achieved for low angles 2i ?n .
5.2Radial accuracy analysis
Following Fig.5,￠r P P the radial spacing in the pseudo-polar grid is constant along each ray,while ￠r P P changes from ray to ray.For the j th ray we have
Figure 9shows ￠r P P as a function of the ray index.On the polar grid ￠r is constant
￠r (j )=18j:
The radial spacing smoothly changes as a function of the ray angle,where the …nest spacing ￠r P P =1is achieved for low angles μ?!?
2k;k 2Z .
Figure 9:The pseudo-polar grid spacing ￠r as it changes from ray to ray.
5.3Accuracy analysis of the polar and log-polar DFT interpolation using the
pseudo-polar FFT
We compare the accuracy of approximating the polar and log–polar DFT magnitudes using the PPFFT and FFT evaluated on Cartesian grids.Denote by "P and "LP ,the approximation errors of the polar and log-polar DFTs’magnitudes,respectively
"P=LP ?μi ;r j ￠=3c M P=LP 1?μi ;r j ￠?M P=LP 1?μi ;r j
￠′2(5.5)
?@2@μ@r M P olar 1
"2Grid ?μi ;r j ￠
where M P=LP
1?μi;r j￠and c M P=LP1?μi;r j￠are the magnitude and approximated magnitude of the polar/log-polar DFT at a spatial frequency?μi;r j￠.The interpolation error"P=LP,de…ned in Eq.
(5.5),depends on the divergence of the pseudo-polar grid from the polar and log-polar pointset. Thus,the…rst order interpolation error(bilinear interpolation)can be estimated by
"P=LP?μi;r j￠?@2M P olar1"2Grid?μi;r j￠:(5.6) "2Grid?μi;r j￠is the L2distance between a point on the polar grid and the closest point on the pseudo-polar grid
"2Grid?μi;r j￠=3(X P seudo?P olar(i;j)?X P olar(i;j))2+(Y P seudo?P olar(i;j)?Y P olar(i;j))2′2;(5.7) where
f X P seudo?P olar(i;j);X P seudo?P olar(i;j)g=f cos(μi)￠r(j);cos(μi)￠r(j)g(5.8)
f X P olar(i;j);X P olar(i;j)g=?cosμ2?i?￠j;cosμ2?i?￠j?:(5.9) The pseudo-polar coordinatesμi and￠r(j)are de…ned in Eqs.(3.4)and(3.6),respectively.
The second order error(cubic interpolation)is given by
"P=LP?μi;r j￠?"3Grid?μi;r j￠:(5.10) In order to evaluate the grid resampling accuracy of the…rst and second order interpolations,the approximation errors are given by
"1Grid=P i;j"2Grid?μi;r j￠ˉˉr2M?μi;r j￠ˉˉ
"2Grid=P i;j j"Grid j3?μi;r j￠ˉˉr3M?μi;r j￠ˉˉ(5.11) where M?μi;r j￠is the magnitude of the2D FFT.We can look at Eq.(5.11)as a weighted sum of grid errors where the weights are the magnitudes’derivatives.The weighted formulation is signi…cant since most of the energy of natural image’s spectra is concentrated at low frequencies.
The testbed for evaluating Eq.(5.11)was the lena image whose log-magnitude is presented in Fig.10,where the energy concentration in low frequencies is evident.The error estimations"1Grid and"2Grid were computed in the following way:lena0s DFT magnitudes were calculated on grids of sizes64￡64;128￡128and256￡256.The derivatives of these magnitudes were computed and used to numerically evaluate Eq.(5.11).The results are presented in Table2.
Table2:The reduction of the approximation error achieved by the pseudo-polar based algorithm. The error’s decrease is computed relatively to the error of the polar/log-polar grid approxima-tion algorithms using a Cartesian DFT.The proposed pseudo-polar based algorithm improves the approximation accuracy.Better results are archived by higher order approximations,while using bigger grids in most cases does not increase the relative gain.
6Experimental results
The proposed registration algorithm was tested using the Lena and Air…eld images,which were rotated,scaled and translated(similar to[16,20])to create the test pairs shown in Figs.11and12. The images were registered using the rotation and rotation/scale algorithms presented in sections 4.2and4.3,respectively.The results are illustrated by overlaying the edges of one image over the other according to the estimated motion parameters,while the corresponding numerical results are presented in tables3and4,respectively.The translations were estimated by the phase-correlation algorithm discussed in sections2.1and4.1.1.We note that the choice of integral scaling factors for the experiments is not signi…cant,the proposed algorithm was able to recover any scaling factor within its dynamic range that is discussed later.
Table3and Fig.11present the registration results of the Air…eld image,which contains man-made objects characterized by sharp edges surrounded by smooth regions.The algorithm was able to register the images in Fig.11a having a scaling factor of2and a large rotation of145±.This test was repeated for various rotation angles,resulting in similar accurate registrations.The registration of these images was not in?uenced by the image boundaries,whose alignment corresponds to zero motion.
Figure11b shows image sets with partial alignment,which were registered by the proposed algorithm.The results are attributed to the statistical interpretation of the DFT’s magnitude as the Periodogram[23]of the images,which can be considered an approximation of the image’s