According to Euler’s formula and Equation 1, the phase
shift angle (phase angle of
Q
(
u
,
v
)) represents a 2D plane
defined by shift parameters
x
0
and
y
0
in the Cartesian coordi-
nates of
u
-
v
, which can be defined as:
∠
Q
(
u,v
) =
ux
0
+
vy
0
(2)
Consequently, the subpixel shifts
x
0
and
y
0
can be obtained
by solving the plane parameters in Equation 2. However, the
phase shift angle is inevitably contaminated by various cor-
ruptions such as image noise and aliasing effect. Based on this
observation, Stone
et al.
(2001) proposed a subpixel phase
correlation method that finds the slope of plane by a least
square estimate after masking the high-frequency components
and the frequency components with small magnitude. Stone’s
method considers shift estimation as a plane fitting problem
in a theoretical sense, and the proposed frequency masking
operator can reduce the influence of corruptions to some
extent. However, the frequency masking processing is inad-
equate and relies on selection of specified thresholds control-
ling masking, which makes the least square estimates and the
shift results unreliable and unstable. The practical perfor-
mance of Stone’s method is thus questioned, especially in the
case of complicated conditions with high noise (González
et
al.
, 2010; Tzimiropoulos
et al.
, 2011; Tong
et al.
, 2015).
Therefore, an improved subpixel phase correlation method
is introduced in this study in order to ensure the robustness
and practicality. Considering the advantages and drawbacks
of the original Stone’s method, the improved method inte-
grates the plane fitting manner and frequency masking with
four additional measures, i.e., the image gradient representa-
tion, phase filtering, an effective robust estimation and robust-
ness iteration operation. These four additional measures are
described as follows:
1. Similar as Tzimiropoulos
et al.
(2011) and Ren
et al.
(2014), the image gradients are adopted to perform subpix-
el phase correlation. The original image intensity function
(e.g.,
f
(
x,y
) and
g
(
x,y
)) is replaced with a gradient-based
complex representation:
Gr
(
x,y
) =
Gr
x
(
x,y
) +
jGr
x
(
x,y
)
(3)
where
Gr
x
=
∇
x
I
and
Gr
y
=
∇
y
I
are the image gradients along
each direction. The gradients can be calculated using
first-order or second-order central difference and Savitzky-
Golay differentiators. After constructing the gradient
representations, the normalized cross-power spectrum
matrix can be calculated by Equation 1 and applied to the
following steps. The joint use of gradient magnitude and
orientation solely keeps the frequency response of the sa-
lient features and suppresses the effect of useless uniform
area, dissimilar part and noise, which makes the results
more accurate and robust.
2. Phase filtering operation is carried out to weaken the influ-
ence of the corrupted phase values on the following phase
unwrapping and plane fitting. As the wrapped phase shift
angles are discontinuous, filtering directly on the phase
shift data is infeasible. Similar as Tong
et al.
(2015), the
real and imaginary parts of the normalized cross-power
spectrum matrix are smoothed respectively. The neighbor-
hood size of the smoothing filter is selected small to not
alter the slope.
3. A robust model fitting algorithm using higher than mini-
mal subset sampling (
HMSS
) (Tennakoon
et al.
, 2016) is
applied to estimate the plane parameters. This is an im-
proved version of the random sample consensus (
RANSAC
)
algorithm, which utilizes the hypothesize-and-verify
strategy. There are two main modifications compared to
the standard
RANSAC
algorithm. On the one hand, hypoth-
eses generated by minimal subsets can be far from the
underlying model due to the influence of measurement
noise, even when the minimal subsets contain purely inli-
ers. The strategy of higher than minimal subset sampling
is thus carried out in this algorithm. On the other hand,
the least
k
-th order statistic is selected as the cost function
to make the threshold for the inlier/outlier discrimination
adaptive (Moisan
et al.
, 2012), and a new stopping criteri-
on is devised to stop the algorithm once a good estimate is
reached. Both of these points make this algorithm accurate
and efficient when compared to the state-of-the-art robust
model fitting algorithms and suitable to our phase correla-
tion method. The
HMSS
algorithm can adaptively exclude
the corrupted phase values and provide the precise and
robust estimation.
4. The robustness and reliability of the subpixel phase cor-
relation method can be improved by iterating it (Leprince
et al.
, 2007). Once the first shift estimations (
x
1
0
,
y
1
0
,) have
been achieved, the new normalized cross-power spectrum
matrix can be defined as:
Q u v Q u v
i ux vy
2
1
0
1
0
1
,
, exp
( )
=
( )
− +
(
)
{
}
(4)
The process can be seen as an image resampling by com-
pensating the achieved shifts in the frequency domain accord-
ing to the Fourier shift property. The residual shifts (
x
2
0
,
y
2
0
,)
are possible to be calculated from
Q
2
(
u,v
) by implementing
the frequency masking and filtering as well as robust estima-
tion again. The robustness iteration can be performed until
the residual shifts become lower than a preset threshold, and
the final shift estimation is obtained by summing the iterative
shift estimation. In practice, it is found that one iteration is
enough in consideration of the trade-off between the robust-
ness improved and the computational cost.
In addition, phase unwrapping is not considered in original
Stone’s method, so only subpixel shifts with range in -0.5 to
0.5 pixels can be calculated. According to Foroosh and Balci
(2004), the ill-posed 2D phase unwrapping process is theoreti-
cally separable, and can be realized using two independent
and consecutive 1D unwrapping along two directions.
The workflow of the improved method is depicted in
Figure 1 and implemented as follows: (1) Gradient representa-
tions are calculated using second-order central difference, and
a raised-cosine window is adopted to reduce the influence of
the edge effects; (2) Normalized cross-power spectrum matrix
Q
between the two images is obtained through
FFT
; (3) The
real and imaginary parts of
Q
are separately smoothed with a
small neighborhood size, and the useless spectral components
at each periphery of
Q
and with small magnitudes are masked
out. The filtered phase shift angles are then unwrapped
consecutively in two directions; (4) The
HMSS
algorithm is uti-
lized on the unwrapped data to find the best plane fitting. The
slope coefficients, which indicate the subpixel shifts between
the two images, are refined by least squares minimization
based on the inliers obtained from the robust model fitting
algorithm; and (5) The new normalized cross-power spectrum
matrix is computed as Equation 4, and the robustness itera-
tion is carried out once to further refine the shift estimation.
3D Videogrammetric Measurement System
As a developing measurement technique, there is a need to
thoroughly investigate the videogrammetric technique by
employing outstanding image processing algorithms. An
improved videogrammetric system is introduced based on the
PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING
September 2018
581
