Statistical Properties of the Hough Transform Estimator in the Presence of Measurement Errors Itai Dattner University of Haifa
Statistics Seminar The William Davidson Faculty of Industrial Engineering and Management, Technion August 4, 2008
Contents 1. Introduction and motivation • The Hough transform • Literature • Observations model 2. Asymptotic results • Consistency • Rate of convergence • Limit distribution 3. General properties and illustrations • Breakdown point • Equivariance • Choice of the parameter r 4. Conclusions
3
INTRODUCTION and MOTIVATION
4
The Hough transform - the technique
(a)
(b)
4
4 3.5
3
3
t0
2.5
1
t
Y
2
φ0
0
(φ0, t0)
2 1.5 1
−1
−2 −2
Xcosφ0+Ysinφ0=t0 −1
0
0.5
1
X
2
3
4
0 0
pi
φ
2pi
Figure 1: (a) Original domain. (b) Hough domain.
5
Example - Test image
NY.gif
Figure 2: A test image.
6
Example (cont.) - Edge data
(a)
(b)
NY2.gif 5
10
15
20
25 5
10
15
20
25
30
Figure 3: (a) Relevant data from the image of Figure 2. (b) The data after edge detection.
7
Example (cont.) - Hough space 30
25
t
20
15
10
5
0
0
pi φ
2pi
Figure 4: The Hough transform for the data points in Figure 3(b).
8
Example (cont.) - Accumulator array
25
20
H(φ, t)
15
10
5
0 40 2pi 20 pi
t
φ 0
0
Figure 5: The accumulator array corresponding to Figure 4. ˆ tˆ) = (φ24 , t79 ) = ( π , 23.0057) (φ, 4.2553
9
The Hough transform - definition • The transform (Xi , Yi ) 7→ Ci (φ, t) := {(φ, t) : Xi cos φ + Yi sin φ = t} • The rectangular parameter cell centered at θ := (φ, t) S(θ) = {θ0 = (φ0 , t0 ) : |φ0 − φ| ≤
δt δφ , |t0 − t| ≤ } 2 2
• The HT estimator ˆ tˆ) = arg θˆn = (φ,
n
max
θ∈[0,2π)×R1+
1X 1{S(θ) ∩ Ci (θ) 6= ∅} n i=1
10
Literature - Computer vision • Basic papers – Cartesian parameterization - Hough (1959) – Polar parameterization - Duda and Hart (1972)
• Algorithm enhancement – RHT/PHT - Xu et al. (1990), Kiryati et al. (1991, 2000)
• Detection of general shapes – Circular arcs - Kimme, Ballard and Sklansky (1975) – Arbitrary shapes - Ballard (1981)
• Applications of the HT – Eyes and eyelids localization - Dobes et al (2006) – Inspection of rice seed - Cheng and Ying (2004)
11
Literature - Statistics - 1 • Dual plot – m test - Daniels (1954) – Resistant line - Johnstone and Velleman (1985) – Regression depth - Rousseeuw and Hubert (1999) 3.5 3
intercept
2.5 2 1.5 1 0.5 0 −0.5
0
0.5
1
1.5
2
2.5
3
slope
Figure 6: Dual plot
12
Literature - Statistics - 2 • Mode estimation - Chernoff (1964) – X1 , . . . , Xn i.i.d. f (a unimodal density). – For a given point θ ∈ R1 and r > 0 define n
Mr,n (θ) θˆr,n
=
1X 1{Xi ∈ [θ − r/2, θ + r/2]}, n i=1
= arg max1 Mr,n (θ). θ∈R
³
´1/3 d – Then nC(r, f ) (θˆr,n − θ0 ) → W, n → ∞, where W is a maximizer of the random process Z(z) − z 2 , and Z(·) is the standard two–sided Brownian motion.
13
Literature - Statistics - 3 • HT - Goldenshluger and Zeevi (2004) – Yi = a0 + b0 Xi + ²i , i = 1, ..., n,, θ0 = (a0 , b0 ). – X is independent of ², and ² is a random variable with bounded, symmetric and strictly unimodal density f² which is continuously differentiable with bounded first derivative. – θˆr,n = arg maxθ∈R2
1 n
Pn
2 2 2 1 {|X b + a − Y | ≤ r (X i i i + 1)}. i=1
– Cube-root asymptotic.
14
Observations model - set up - 1 • Unobservable {X˜i , Y˜i }ni=1 i.i.d. random variables lying on the unknown line X˜i cos φ0 + Y˜i sin φ0 = t0 . • Observations Xi
=
Yi
=
X˜i + γi , Y˜i + ηi ,
˜i , Y˜i ). where γi and ηi are noise variables independent of (X • Objective estimate θ0 = (φ0 , t0 ) from the observed data.
15
Observations model - set up - 2 • Notation ˜ Y˜ )T , Z = (X, Y )T , θ = (φ, t)T , λ = (γ, η)T , Uφ = (cos φ, sin φ)T Z˜ = (X, • Template with line segment cell mθ (Zi ) = 1{|ZiT Uφ − t| ≤
r } 2
• The criterion function n
1X Mr,n (θ) = mθ (Zi ) n i=1 • The HT estimator θˆr,n = arg
max
θ∈Θ=[0,2π)×R1+
Mr,n (θ)
16
The Hough transform - Geometric representation Dθ = {(x, y) : |x cos φ + y sin φ − t| ≤ 2r }, θ = (φ, t) ∈ [0, 2π) × R1+ (a)
(b)
4
4
3.5
3.5
Dθ
3
2.5
2.5
2
2
t
Y
3
1.5
1.5
1
1
0.5
0.5
0 0
0.5
1
1.5
2
X
2.5
3
S(θ)
3.5
4
0 0
pi
φ
2pi
Figure 7: (a) Original domain. (b) Hough domain.
17
ASYMPTOTIC RESULTS
18
Technical tools for the study of asymptotic properties • M-estimators - Huber (1981), Van der Vaart and Wellner (1996), Chapter 3.2 – M-estimators is defined as
θˆ = arg max θ∈Θ
n X
mθ (Zi ),
i=1
– Zi i.i.d. random variables, – mθ (Z) = log f (Z, θ), – f is the joint probability density of observations Z , Then θˆ is the maximum likelihood estimator of θ0 .
19
Consistency • Theorem 1 Let
˜ n * {λi }n i=1 be i.i.d. random vectors independent of {Zi }i=1 ; * for any φ ∈ [0, 2π), the random variable λT Uφ has density f which is symmetric and strictly unimodal. Then for any fixed r > 0, p ˆ θr,n → θ0 , as n → ∞.
20
Consistency - proof scheme • Deterministic condition - θ0 is a well separated point of maximum of Mr (θ) = EMr,n (θ):
sup
Mr (θ) < Mr (θ0 ).
θ:kθ−θ0 k≥²
– The density of λT Uφ is f , thus symmetric and unimodal. Anderson (1955):
Mr (θ) =
0,
n1/3 kθˆr,n − θ0 k = Op (1), as n → ∞.
23
Cube root asymptotics - Kim and Pollard (1990) Recall mode estimation; X1 , . . . , Xn i.i.d. f a unimodal density: n
Mn (θ) M (θ)
=
1X 1{Xi ∈ [θ − 1, θ + 1]}, n i=1
= EMn (θ) = P(θ − 1 ≤ X ≤ θ + 1).
* M (θ) − M (θ0 ) ≈ −c1 (θ − θ0 )2 . R θ0 −1 R θ0 +1 1 * Mn (θ) − Mn (θ0 ) = n [Bin(n, θ−1 f (z)dz) − Bin(n, θ+1 f (z)dz)]. Thus [Mn (θ) − Mn (θ0 )] − [M (θ) − M (θ0 )] ∼ N (0,
c22 |θ−θ0 | ). n
Conclude that Mn (θ) − Mn (θ0 ) ≈ c2 n−1/2 |θ − θ0 |1/2 − c1 (θ − θ0 )2 is maximized for |θ − θ0 | of order n−1/3 .
24
Limit distribution • Theorem 3 ˜ 2 , EY˜ 2 < ∞. Let λi ∼ N2 (0, σ 2 I), i = 1, ..., n; EX Then for any fixed r > 0, d n1/3 (θˆr,n − θ0 ) → W, as n → ∞,
where W has the distribution of the maximizer of the process
1 θ 7→ G(θ) + θT V0 θ, 2 ˜ T (Z) ˜ – V0 = [f 0 ( 2r ) − f 0 (− 2r )]EA(Z)A ˜ = (Z˜ T Uφ −π/2 , 1)T , – A(Z) 0 ¡r¢ 2 ˜ – E[G(g) − G(h)] = 2f 2 E|(h − g)T A(Z)|; g, h ∈ Θ.
25
GENERAL PROPERTIES AND ILLUSTRATIONS
26
Breakdown point - Donoho and Huber (1983) • Addition breakdown point n
o k ˆ Zn ) := min ˆ n ∪ Z 0 ) − θ(Z ˆ n )k = ∞ ²add (θ; : sup kθ(Z k n + k Zk0 • Replacement breakdown point ˆ Zn ) := min ²rep (θ;
nk n
o ˆ n )k = ∞ ˆ k ) − θ(Z : sup kθ(Z n k Zn
27
Breakdown point • Theorem 4 Let * Zn = {Z1 , ..., Zn } be a sample belonging to a ball of finite radius, Then – ²add (θˆr,n ; Zn ) = – ²rep (θˆr,n ; Zn ) =
nMr,n (θˆr,n )+1 n+nMr,n (θˆr,n )+1 1 n
l
nMr,n (θˆr,n ) 2
m
Moreover a.s. – ²add (θˆr,n ; Zn ) → p(1 + p)−1 , n → ∞, a.s. – ²rep (θˆr,n ; Zn ) → p/2, n → ∞,
where p =
R r/2
−r/2
f (x)dx.
28
Breakdown point - illustration
(a)
(b)
12
7 data HT LS
data HT LS 6
10
5 8
Y
Y
4 6
3 4 2
2
0
1
0
1
2
3
4
5 X
√
6
7
8
9
10
0
0
1
2
3
2
4
5 X
6
7
8
9
10
R 0.3
Figure 8: θ0 = (3π/4, 2), f = N (0, 0.2 ), r = 0.6, −0.3 f (x)dx = 0.8664 corresponds to a 43.32% replacement BP. (a) 40% of contamination - θˆ0.6,50 = (2.3483, 1.4176). (b) 48% of contamination.
29
Rotation equivariance • Rotation matrix Rα
• Transformed sample Z ∗ = Rα Z
• Because r r T 1{|(Rα Z) Uφ+α − t| ≤ } = 1{|Z Uφ − t| ≤ } 2 2 T
Therefore
ˆ ∗ ) + α = φ(Z) ˆ tˆ(Z ∗ ) = tˆ(Z), φ(Z
30
Rotation equivariance - illustration
(a)
(b) 2
2 Original data Transformed data HT HT
1.5
1.8
1.6
1
1.4
1.2
t
Y
0.5 1
0 0.8
0.6
−0.5
0.4
−1 0.2
−1.5 0.8
0
1
1.2
1.4 X
1.6
1.8
2
0
pi φ
2pi
ˆ ∗ ) + 0.3808 = φ(Z) ˆ . Figure 9: α = π/8 = 0.3927, tˆ(Z ∗ ) = tˆ(Z), φ(Z (a) Blue cloud is the original data set. (b) Hough domain.
31
Choice of the parameter r Objective: rˆ = arg minr∈R tr{cov(W )}
(a)
(b)
1
(c)
0.9
0.07
0.8 0.8
0.06
0.7 0.05 0.6
0.4
0.5
H(φ, t)
H(φ, t)
H(φ, t)
0.6
0.4
0.04 0.03
0.3 0.02 0.2
0.2
0.01
0.1 0 1
0 1 0.8
2pi 0.6
0 1 0.8
2pi 0.6
0.4
pi
2pi 0.6
0.4
0.2 t
0.8
pi
0.4
0.2 0
0
φ
t
pi 0.2
0
0
φ
t
0
0
φ
Figure 10: (a) large r. (b) ”optimal” r. (c) small r.
32
Choice of the parameter r - experiment set up • r ∈ R = [0.276 : (0.025) : 0.976] • Sample size n = 200 • Simulations i = 1, ..., N = 200 • Experiments j = 1, ..., M = 192 i,j is the estimator computed in the i-th simulation in • θˆr,n experiment j PN ˆi,j • S(r, j) = i=1 kθr,n − θ0 k2
• rˆj = arg minr∈R S(r, j) • φ0 = π/4, t0 = 0.5, λi ∼ N2 (0, 0.22 I) for all i.
33
Choice of the parameter r - experiment results
(a)
(b)
(c)
50
0.85 45
0.8
7 6.5
40
0.75 6
35
0.7
25
20
5
0.65 S(r, j)
estimated r
30
5.5
0.6
4.5 4 3.5
0.55
3
15
0.5
10
0.45
5
0.4
2.5 2 150 100 0.976
50
0 0.3
0.4
0.5
0.6
0.7 estimated r
0.8
0.9
1
1 boxplot
0.626 0
j
0.276
Figure 11: (a) rˆj , j = 1, ..., 192. (b) Boxplot. (c) S(r, j).
r
34
CONCLUSIONS
35
Conclusions • Consistency, rate of convergence n1/3 • No analytic expression for the limit distribution • Statistically inefficiency in comparison with the least squares estimator but very good robustness properties • Numerical experience suggests r = 3σ
36
This work is based on the M.A. thesis written under the supervision of A. Goldenshluger