著者原稿 - 画像情報メディア研究室

情報処理学会研究報告, 2005-CVIM-151-15, 2005/11/17,18 pp. 107–114.
107
超精度の楕円当てはめ
山田 純平
金谷 健一
岡山大学大学院自然科学研究科
平面上の点列に楕円を当てはめる問題は,従来から最尤推定が最も高精度であるとみなされ,くりこみ法,HEIV 法,
FNS 法などの計算法が提案されている.本論文ではこれを上回る「超精度」の当てはめ法が存在することを理論的お
よび実験的に検証する.これは最尤推定の摂動解析によって高次の偏差項を評価し,これを最尤推定解から差し引くも
のである.ただし,最尤推定解は既に精度に理論限界である KCR 下界をほぼ達成しているので,補正による効果は非
常にわずかである.このため,提案方法によって実際の応用で実質的な差が生じることはないが,当てはめの精度の理
論限界と最尤推定に関して,理論的な観点から意味のある結果である.
Ellipse Fitting with Hyperaccuracy
Junpei Yamada,
Kenichi Kanatani
Department of Computer Science, Okayama University, Okayama 700-8530 Japan
For fitting an ellipse to a point sequence in 2-D, maximum likelihood (ML) estimation has been regarded as
having the highest accuracy, and numerical schemes such as renormalizatin, HEIV, and FNS have been proposed
for computing the solution. In this paper, we demonstrate, theoretically and experimentally, the existence of
a “hyperaccurate” method having higher accuracy than ML. The hyperaccuracy is achieved by perturbation
analysis of the ML solution followed by evaluation of the high-order bias terms and subtraction of them from
the ML solution. However, since the ML solution almost attains the theoretical accuracy bound (the KCR lower
bound), the resulting improvement is too small to produce noticeable difference in practical applications. Yet,
our analysis has theoretical significance, illuminating the relationship between accuracy of the ML solution and
the theoretical accuracy bound.
1. まえがき
シーン中の円形や球形の物体をカメラで撮影する
と一般に楕円に投影され,その投影像からその物体
の 3 次元位置が解析できる [19].このため,画像から
抽出した点列に円や楕円を当てはめることは,視覚
ロボットを含む広範な応用の基本的な処理の一つで
ある.実際,楕円を当てはめに関してこれまで多く
の論文が書かれている.それらは次の二つに大別で
きる.
1. 画像から抽出したエッジ点列の中から,円や楕
円をなすものをどう判別するか.一つの点列に
円または楕円以外の点(アウトライア)が含ま
れているかをどう判定するか.
2. 円または楕円上にあることが既知の点列(イン
ライア)に対して,データに誤差があるとき,な
るべく正確に円または楕円の方程式を当てはめ
るにはどうすればよいか.
51, 56, 57, 58, 59, 61, 62, 63, 64]
,最近では [37] がある.
第 2 については,以前は発見的な方法(投票法,最
小二乗法など)が主であったが [4, 5, 18, 29, 41, 42, 48,
49, 52, 53, 54, 55]
,これを誤差のあるデータからの統
計的推定とみなして最適な推定量を求めようという
理論研究が増えた.しかし,多くはデータ数が増大
するときの推定量の一致性や有効性の解析を中心で
ある [1, 3, 6, 30, 31, 60] .
それに対して,筆者らは画像処理を念頭において,
データ数を固定した少数のデータに対して,画像の
誤差を小さくすると急速に精度が向上する手法を追
求した [23].それは,そのような手法は望ましい精度
を実現するのにより大きなデータの誤差を許容でき
るからである [24].
筆者らはさらに,楕円当てはめを幾何学的当ては
めという一般論に抽象して,どのような手法を用い
ても超えることのできない精度の理論限界を共分散
行列の下界という形として示した [22, 23].Chernov
第 1 については,古くからさまざまなアルゴリズ
ら [8] はこれを KCR(金谷・クラメル・ラオ)下界
ムとその効率的な処理方法が研究され [2, 7, 10, 12, 13,
と呼び,より弱い仮定から同様の結果が得られるこ
14, 15, 16, 17, 33, 35, 36, 38, 39, 40, 43, 44, 45, 46, 47, 50,
† 700-8530
岡山市津島中 3–1–1 Tel/Fax: (086)251-8173
{yamada,[email protected]
とを指摘した.
= 1, ..., N は式 (2) のよって R6 の N 点 {ξα } に埋
の高次の項を除いてこの下界を達成することが示さ め込み写像されるので,楕円当てはめは R6 での超
れる [8, 23, 26].これまでに,くりこみ法 [20, 23, 25], 平面を当てはめる問題とみなせる.
HEIV 法 [32],および FNS 法 [9] と呼ぶ反復手法が
ただし,式 (1) で表されるのは楕円とは限らず,放
提案されているが,いずれも理論的にはその精度が 物線,双曲線,およびそれらの退化も含んでいる.こ
“最尤推定” と同等であることが示される [27].この れらは総称してコニック(円錐曲線)と呼ばれるの
意味で楕円当てはめは実用的な観点からは既に解決 で,上述の定式化は「コニック当てはめ」とも呼ばれ
済みの問題であるといえる.
る [19].実際,誤差が大きいと楕円上の点列に双曲線
それに対して本論文では理論的な観点から,最尤 が当てはまったりする.それを防ぐ方法もあるが [11],
推定を上回る解法は本当に存在しないのかという問 本論文ではそのような場合も解に含めて解析する.
題を探求する.便宜上,最尤推定と同等の精度を「高
3. 最尤推定
精度」と呼び(くりこみ法,HEIV 法,FNS 法など),
上記の問題を解くために最も広く用いられている
それに及ばないものを「低精度」と呼ぶ(最小二乗
法,Taubin 法 [54] など).両者の差は歴然としてお のは次式を最小にする「最小二乗法」である.
N
り,高精度手法を用いることによって精度が著しく
JLS =
(u, ξα )2
(4)
向上することが実験的に確認されている [21, 28, 34].
α=1
これら対比して,最尤推定を上回る精度を「超精度」
これは
と呼ぶことにする.
N
M
=
ξα ξα
(5)
LS
本論文では超精度を達成する手法が実際に存在す
α=1
ることを理論的に示し,実験によって確認する.た
と定義すれば単位ベクトル u の2次形式 JLS =
だし,最尤推定は誤差の高次の項を除いて精度の理
(u, M LS u) に書き直せるので,解(「最小二乗解」)
論限界を達成しているので,精度の向上の程度は非
ˆ LS は M LS の最小固有値に対する単位固有ベクト
u
常に小さいものである.本論文では実験によってこ
ルである.しかし,最小二乗解には大きな偏差が存
れがどの程度のものかを確認する.
在して精度が低いことが知られている [23].
2. 楕円当てはめ
それに対して「最尤推定」は,各点 ξ α からその超
平面上の N 点 {(xα , yα )}, α = 1, ..., N に楕円を 平面までのマハラノビス距離の二乗和を最小にする
当てはめる問題を考える.楕円は次のように表せる. ように超平面を当てはめるものである [23].すなわち,
この枠組によると,“最尤推定” と呼ぶ方法が誤差
N
2
Ax + 2Bxy + Cy
2
+ 2f0 (Dx + Ey) + F f02
(ξ α − ξ¯α , V0 [ξ α ]− (ξ α − ξ¯α ))
J=
= 0 (1)
(6)
α=1
ただし,f0 は任意のスケール定数である.新しい係
¯ ) = 0, α = 1, ..., N のもとで最小
を拘束条件 (u, ξ
α
数ベクトル u と変数ベクトル ξ を
にする.ただし,V0 [ξ α ] はデータ点 {(xα , yα )} の誤
u=
A B
C
D
E
差から導かれる ξ α の共分散行列を正規化したもの
F
である.正規化とは,共分散行列を定数倍しても式
(6) を最小にする解は変わらないことから,誤差の大
(2) きさに依存する定数を除いて {(x , y )} の式のみで
α α
表したものである.変数 ξ は 2 自由度しかないから
と置けば,式 (1) は次のように書ける.
V0 [ξ α ] は一般にランクが 2 であり,V0 [ξα ]− はその
(u, ξ) = 0
(3) 一般逆行列を表す.ラグランジュ乗数を導入して拘
束条件を除去すれば,式 (6) は次式となる [23].
ただし,以下ベクトル a, b の内積を (a, b) と書く.
N
ベクトル u の絶対的な大きさは不定であるから,以
(u, ξα )2
(7)
J=
(u, V0 [ξα ]u)
後 u = 1 と正規化する.
α=1
ξ=
x2
2xy
y2
2f0 x 2f0 y
f02
幾何学的には式 (3) は ξ の 6 次元空間 R6 の超平
面を表している.誤差を含むデータ点 {(xα , yα )}, α
各データ点 (xα , yα ) はその真の位置 (¯
xα , y¯α ) に平均
0,標準偏差 σ の正規分布に従う誤差が独立に加わっ
108
たものであるとすると,ξ α の共分散行列は 4σ 2 V0 [ξ α ]
u
となり,V0 [ξ α ] は σ の高次の項を除いて1 ,次のよう
u
に書ける.

x
¯2α
x
¯α y¯α
0
f0 x
¯α
0
0



x
O
¯2α + y¯α2 x
¯α y¯α f0 y¯α f0 x
¯α 0 
 ¯α y¯α x

 0

2
ˆ の誤差の直交成分と平行成分.
図 1: 推定量 u
x
¯α y¯α
y¯α
0 f0 y¯α 0 

V0 [ξ α ] = 

2
 f0 x
0
f0
0 0
¯ は真の値 {ξ¯α } を用いて計算した M の
ただし,M
 ¯α f0 y¯α

 0

2
f0 x
¯α f0 y¯α 0
f0 0 

値であり,∆1 M , ∆2 M はそれぞれ {∆ξ α } の1次の
0
0
0
0
0 0
項を含む部分,および2次の項を含む部分である.
(8)
式 (11) の行列 L については次のようになる.
4. 最尤推定の誤差解析
N
N
(ξ¯α+∆ξ α ,u)2 V0 [ξα ]
(∆ξα , u)2 V0 [ξ α ]
=
(u, V0 [ξ α ]u)2
(u, V0 [ξ α ]u)2
α=1
α=1
L=
式 (7) の微分は次のようになる.
N
= ∆2 L
(16)
2(ξ α , u)ξ α
2(ξ α , u)2 V0 [ξ α ]u
−
(u, V0 [ξ α ]u) α=1 (u, V0 [ξ α ]u)2
¯ = O, ∆1 L = O である.以上より,式
α=1
すなわち,L
(9)
(10) に式 (12) を代入すると次のように書ける.
ˆ
ただし,∇u は u に関する勾配を表す.最尤推定量 u
¯ + ∆1 M + ∆∗ M + ∆2 M + ∆∗ M )(u + ∆1 u
は上式を 0 にするものであるから,次式の解である. (M
N
∇u J =
1
M u = Lu
2
+∆2 u + · · ·) = ∆2 L(u + ∆1 u + ∆2 u + · · ·)
(10)
(17)
¯ , ∆1 M の分母の u を u
N
ˆ
∆∗1 M , ∆∗2 M はそれぞれ M
ξα ξα
(ξ α , u)2 V0 [ξα ]
M=
, L=
に置き換えるために生じる摂動項であり,次のよう
(u, V0 [ξ α ]u)
(u, V0 [ξ α ]u)2
α=1
α=1
(11) になる(∆2 M の摂動項は 3 次以上の高次となる).
Chojnacki ら [9] の FNS 法も,Leedan ら [32] の HEIV
N
((∆1 u, V0 [ξ α ]u) + O(σ 2 ))ξ¯α ξ¯α
∆∗1 M = −2
法も,式 (10) を固有値問題,あるいは一般固有値問
(u, V0 [ξα ]u)2
α=1
題の反復によって解く方法である.また,くりこみ
(18)
法 [20, 23, 25] も理論的も理論的に精度が同等である
2
((∆1 u, V0 [ξ α ]u) + O(σ ))
ことが示される [27].
N
×(∆ξα ξ¯α + ξ¯α ∆ξ α )
データに誤差がないときの真の解を u とし,式 (10)
∆∗2 M = −2
(19)
(u, V0 [ξ α ]u)2
α=1
ˆ を次のように展開する.
を満たす最尤推定量 u
N
ˆ = u + ∆1 u + ∆2 u + · · ·
u
式 (17) の両辺の O(1), O(σ), O(σ 2 ) の項をそれぞ
(12)
れ取り出すと,次の結果を得る(導出は付録に示す).
ただし,∆k u は ∆ξ α の成分の k 乗を含む O(σ k ) の
¯ +∆ξ
項を表す.式 (11) の行列 M の右辺に ξ = ξ
¯ − ∆1 M u
∆1 u = −M
を代入すると次のように書ける.
¯ ∆2 M u + M
¯ ∆1 M M
¯ ∆1 M u
∆2 u = −M
α
¯ + ∆1 M + ∆2 M
M =M
α
α
∆ξα ξ¯α + ξ¯α ∆ξ α
(u, V0 [ξ α ]u)
α=1
N
∆2 M =
∆ξ α ∆ξα
(u,
V0 [ξα ]u)
α=1
(14)
(20)
−
−
¯ − ∆∗ M ∆1 u − M
¯ − ∆∗ M u
−M
1
2
(13)
N
∆1 M =
−
¯ − ∆2 Lu − M
¯ − ∆1 M u 2 u
+M
(21)
¯ u = 0 が,したがって
式 (11) の M の定義から M
−
¯ u = 0 が成り立つ.このため,式 (20), (21) の右
M
辺の項は最後のもの以外はすべて u に直交する.最
(15) 後の項 − M
¯ − ∆1 M u 2 u は u に平行であり,正規
化 u = 1 を強制するための項である(図 1).
1 後述の実験では省略した高次の項の影響がないことを確認し
ている.
109
5. 最尤推定量の最適性
ただし,上式中で各データ xα に入る誤差は互いに
独立という仮定から E[∆ξ α ∆ξ β ] = 4σ 2 δαβ V0 [ξ α ] で
2 節に述べた問題に対して,どのような推定を行っ
あることを用いた2 .
ˆ
ても,データに誤差がある限り,得られる推定量 u
¯ の定義から,式 (25) が式 (24) の右辺に一
行列 M
の共分散行列 V [ˆ
u] には,下回ることのできない下界
致していることがわかる.これに ∆2 u を加えても,
ˆ の共分散行列 V [ˆ
が存在する.推定量 u
u] は次のよ
誤差は正負に対称で ∆ξ α 奇数次の項の期待値が 0 に
うに定義する.
ˆ の共分散行列に対する寄与は O(σ 4 ) で
なるため,u
ˆ の共分散行列は O(σ 4 )
ˆ )(P u u
ˆ) ]
V [ˆ
u] = E[(P u u
(22) ある.すなわち,最尤推定量 u
を除いて KCR 下界を達成しており,この意味で最尤
P u は次のように定義する射影行列である(I は単位 推定量は最適である.
ベクトル).
P u = I − uu
(23) 6. 最尤推定量の偏差
ˆ の偏差を調べる.式 (14) より
次に最尤推定量 u
これは u に直交する平面上へ射影を表す.これを作
E[∆1 M ] = O であるから,式 (20) より E[∆1 u] =
用させるのは,パラメータ u が単位ベクトルに正規
0 である.次に ∆2 u の期待値を計算する.∆2 M の
ˆ も定
化されているためである.その結果,推定量 u
期待値は次のようになる.
義域が単位球面となるが,誤差が小さいとして真値
N
u での接平面に射影して,その接平面上で誤差を評
E[∆ξα ∆ξ α ]
E[∆2 M ] =
= 4σ 2 N
(26)
ˆ が任意の不
価するという意味である.このとき,u
(u,
V
[ξ
]u)
0
α
α=1
偏推定量に対して次の不等式が成り立つことが導か
ただし,次のように置いた.
れる [22, 23].
N
N
V [ˆ
u]
ただし,
4σ
2
ξ¯α ξ¯α
(u, V0 [ξ α ]u)
α=1
−
N=
(24)
は左辺から右辺を引いたものが半正値対
¯ − ∆1 M M
¯ − ∆1 M u]
E[M
上式の右辺を KCR 下界と呼ぶ.Chernov ら [8] はさ
¯
= E[M
ˆ が不偏推定量でなくても,σ → 0 で u
ˆ →u
らに,u
であれば下界が O(σ 4 ) を除いて上式で表されること
N
を示した.
式 (20) の誤差項 ∆u1 がこの KCR 下界に等しい
β=1
−
=
¯ − ∆1 M uu ∆1 M M
¯ −]
E[∆1 u∆1 u ] = E[M
∆ξα ξ¯α + ξ¯α ∆ξ α
uu
(u, V0 [ξ α ]u)
α=1
N
N
= 4σ 2
∆ξβ ξ¯β + ξ¯β ∆ξβ
¯ −]
M
(u, V0 [ξ β ]u)
(u, V0 [ξα ]u)2
(28)
∆∗1 M ∆1 u は次のようになる.
(u, E[∆ξ α ∆ξ β ]u)ξ¯α ξ¯β
¯−
M
(u, V0 [ξα ]u)(u, V0 [ξ β ]u)
N
∆∗1 M ∆1 u = −2
(∆1 u, V0 [ξ α ]u)ξ¯α (ξ¯α , ∆1 u)
(u, V0 [ξα ]u)2
α=1
¯ −∆1 M u,V0 [ξα ]u)(ξ¯α , M
¯ −∆1 M u)ξ¯α
(M
(u, V0 [ξα ]u)2
α=1
N
N
ξ¯α ξ¯α
¯−
M
(u,
V
[ξ
]u)
0
α
α=1
¯ −M
¯M
¯ − = 4σ 2 M
¯−
= 4σ 2 M
¯ − ξ¯α , V0 [ξα ]u)M
¯ − ξ¯α
(M
−
−
¯ ξ¯α )M
¯ V0 [ξ α ]u
+(ξ¯α , M
α=1
α,β=1
¯−
= 4σ 2 M
(u, V0 [ξ α ]u)(u, V0 [ξβ ]u)
α,β=1
N
β=1
¯−
=M
ξ¯α ∆ξ α + ∆ξα ξ¯α ¯ −
M
(u, V0 [ξ α ]u)
α=1
¯ − ξ¯α (M
¯ − ξ¯β ) E[∆ξ α ∆ξ ]u
M
β
−
¯ E[∆ξ ∆ξ ]u(ξ¯ , M
¯ − ξ¯ )
+M
α
α
β
β
N
は次のようになる.
N
ξ¯β ∆ξβ + ∆ξ β ξ¯β
u]
(u, V0 [ξβ ]u)
変動を与えることがわかる.実際,その共分散行列
N
(27)
¯ − ∆1 M M
¯ − ∆1 M u の期待値は次のようになる.
M
称行列であることを意味する.Chernov ら [8] に従い,
¯−
= E[M
V0 [ξα ]
(u, V0 [ξα ]u)
α=1
= −2
(25)
2δ
αβ はクロネッカのデルタであり,α = β のとき 1,それ以
外で 0 をとる.
110
¯ − V0 [ξ ]u, (∆1 M u)(∆1 M u) M
¯ − ξ¯ )ξ¯
(M
α
α α
2
(u,
V
[ξ
]u)
0
α
α=1
N
¯ − ξ¯ )ξ¯
(V0 [ξα ]u, M
α α
2
(u,
V
[ξ
]u)
0
α
α=1
N
= −2
= 8σ 2
(34)
(29) ∆2 L の期待値は次のようになる.
N
(∆1 M u)(∆1 M u) の期待値は次のようになる.
E[∆2 L] = E[
E[(∆1 M u)(∆1 M u) ]
N
ξ¯α (∆ξα , u)
= E[
(u, V0 [ξ α ]u)
α=1
N
=
α,β=1
N
N
β=1
ξ¯β (∆ξβ , u)
]
(u, V0 [ξ β ]u)
=
(u, E[∆ξ α ∆ξα ]u)V0 [ξ α ]
(u, V0 [ξα ]u)2
α=1
N
(u, E[∆ξα ∆ξ β ], u)ξ¯α ξ¯β
(u, V0 [ξ α ]u)(u, V0 [ξ β ]u)
= 4σ 2
¯ − ∆1 M u
M
N
ξ¯α ξ¯β
¯
= 4σ 2
= 4σ 2 M
(u,
V
[ξ
]u)
0
α
α=1
2
V0 [ξα ]
= 4σ 2 N
(u,
V
[ξ
]u)
0
α
α=1
(35)
の期待値は次のようになる.
−
¯ ∆1 M u 2 ] = E[(M
¯ − ∆1 M u, M
¯ − ∆1 M u)]
E[ M
N ¯
N ¯
ξ β (∆ξβ , u)
ξα (∆ξα , u) ¯ − 2
, (M )
)]
= E[(
(u, V0 [ξ α ]u)
(u, V0 [ξβ ]u)
α=1
(30)
したがって ∆∗1 M ∆1 u の期待値は次のようになる.
β=1
E[∆∗1 M ∆1 u]
N
¯ − )2 ξ¯β )
(u, E[∆ξ α ∆ξ β ]u)(ξ¯α , (M
(u, V0 [ξ α ]u)(u, V0 [ξβ ]u)
¯ − V0 [ξ α ]u, M
¯M
¯ − ξ¯α )ξ¯α
(M
= −8σ 2
(u, V0 [ξα ]u)2
α=1
=
¯ −M
¯M
¯ − ξ¯α )ξ¯α
(V0 [ξ α ]u, M
2
= −8σ
(u, V0 [ξα ]u)2
α=1
= 4σ 2
¯ − ξ¯α )ξ¯α
(V0 [ξ α ]u, M
= −8σ 2
(u, V0 [ξα ]u)2
α=1
= 4σ 2 tr(
N
α,β=1
¯ − )2 ξ¯α )
(ξ¯α , (M
(u, V0 [ξ α ]u)
α=1
N
N
N
ξ¯α ξ¯α
¯ − )2 )
(M
(u,
V
[ξ
]u)
0
α
α=1
N
∆∗2 M u
(∆ξα , u)2 V0 [ξα ]
]
(u, V0 [ξα ]u)2
α=1
(31)
¯ (M
¯ − )2 ) = 4σ 2 tr(M
¯ −M
¯M
¯ −)
= 4σ 2 tr(M
は次のようになる.
−
¯ )
= 4σ 2 tr(M
(36)
N
(∆1 u, V0 [ξ α ]u)ξ¯α (∆ξ α , u)
∆∗2 M u = −2
(u, V0 [ξα ]u)2
α=1
7. 超精度補正
ˆ の期待値は次のよ
前節の結果より,最尤推定量 u
¯ − ∆1 M u, V0 [ξα ]u)(∆ξ α , u)ξ¯α
(M
=2
(u, V0 [ξ α ]u)2
α=1
N
−
N
=2
うになる.
¯ V0 [ξα ]u, ∆1 M u∆ξ u)ξ¯α
(M
α
2
(u,
V
[ξ
]u)
0
α
α=1
N
N
E[∆1 M u∆ξα u] = E[
β=1
N
=
β=1
α=1
¯
= 4σ
したがって
(u, V0 [ξ α ]u)
∆∗2 M u
α
−
¯ )u + O(σ 4 )
−tr(M
(37)
ˆ から引けば超精度の当てはめが実現
て最尤推定量 u
¯ − )u は単位
できると期待される.しかし,項 tr(M
ベクトルにするためのノルムの調節であり(図 1),
推定量は単位ベクトルに正規化するのでこの項は考
慮する必要はない.したがって補正量は次のように
2¯
= 4σ ξα
(33) 表せる.
の期待値は次のようになる.
N
∆c u = 4ˆ
σ2
¯ − V0 [ξ ]u, ξ¯ )ξ¯
(M
α
α α
2
(u,
V
[ξ
]u)
0
α
α=1
N
E[∆∗2 M u] = 8σ 2
α
(u, V0 [ξ α ]u)2
この結果から,上式右辺の σ 2 に比例する項を推定し
(∆ξβ , u)ξ¯β ∆ξ α u
]
(u, V0 [ξ β ]u)
ξ¯β (u, E[∆ξ β ∆ξ α ]u)
(u, V0 [ξ β ]u)
2 ξ α (u, V0 [∆ξ α ]u)
α
E[ˆ
u] = u + 4σ 2
(32)
∆1 M u∆ξ α u の期待値は次のようになる,
¯ − ξ¯α , V0 [ξ α ]u)M
¯ − ξ¯α
(M
¯ − ξ¯ )M
¯ − V0 [ξ ]u
+(ξ¯ , M
α=1
(M − ξα , V0 [ξα ]ˆ
u)M − ξα
+(ξα , M − ξα )M − V0 [ξ α ]ˆ
u
(ˆ
u, V0 [ξα ]ˆ
u)2
(38)
111
図 2: 実験に用いた楕円とその上の 20 点.
図 4: 誤差のある点列への楕円当てはめの一例.
0.1
0.05
図 5: 誤差のある点列への楕円当てはめの一例.
対して独立に 10,000 回の実験を行い,次のように評
0
0.01
σ
価した.
0.02
図 3: データの誤差に対する楕円当てはめの誤差.
E=
ただし,式 (37) の右辺は真の値によって表されてい
るものを推定値に置き換えた.すなわち u は最尤推
¯ はデータ {(xα , yα )} に
ˆ に置き換え,行列 M
定量 u
よって定義される ξ α によって計算する行列 M に置
き換える.その計算過程に含まれる共分散行列 V0 [ξ α ]
の計算にも ξ α を用いる.そして分散 σ 2 は次のよう
に推定する.
ˆ)
(ˆ
u, M u
σ
ˆ2 =
4(N − 5)
(39)
1
10000
10000
ˆ (a)
P uu
2
(42)
a=1
ˆ (a) は a 回目の値である.ただし,符号の不
ここに u
定性があるので,(ˆ
u(a) , u) ≥ 0 となる符号を選んだ.
図中の太い実線が最尤推定量に対するものであり,
細い実線が超精度補正を加えたものである.破線は
ˆ LS に対する結果である.点線は式 (24)
最小二乗解 u
の KCR 下界に対応する
N
ˆ )/4σ 2 が自由度 N − 5 の
これは理論的には (ˆ
u, M u
tr
ξ¯α ξ¯α
(u, V0 [ξα ]u)
α=1
−
(43)
χ2 分布に従い,上式のように σ
ˆ 2 を定義すると E[ˆ
σ2 ]
= σ 2 となるからである [23].
をプロットしたものである.
ˆ
このように計算した補正量を用いて最尤推定量 u
図 3 からわかるように,最小二乗解は極めて精度
を次のように補正する.
が低い.それに対して最尤推定量は極めて高精度で
あり,誤差が小さいときは KCR 下界にほぼ一致し
˜ = N [ˆ
u
u − ∆c u]
(40)
ている.誤差がやや大きくなると精度が KCR 下界
からわずかに離れるが,超精度補正を加えると再び
ただし,N [ · ] は単位ベクトルへの正規化を表す.
KCR 下界にほぼ一致している.
8. 実験
図 4 は誤差を加えた点列に楕円を当てはめた一例
図 2 は楕円
である(σ = 0.009).点線が真の楕円であり,破線
が最小二乗解,太い実線が最尤推定量,細い実線が
y2
x2
+
=
1
(41)
超精度補正を加えたものである.この例では超精度
502
1002
補正を加えるとより真の楕円に近づいている.
上に等間隔に N = 20 個の点 {(¯
xα , y¯α )} をとったも
図 5 は別の一例である(σ = 0.009).この例では
のである.各点の x, y 座標に独立に期待値 0,標準
最尤推定量が既に真の楕円にかなり近く,超精度補
偏差 σ の正規分布に従う誤差を加えたものをデータ
正を加えるとかえって真の楕円から離れる.このよ
{(xα , yα )} として楕円の当てはめを行った.最尤推
うな場合も生じるが,全体としては図 4 のように精
ˆ の計算には Chojnacki ら [9] の FNS 法を用
定量 u
度が向上する例のほうが多く,平均すると図 3 に示
いた.
すようなわずかではあるが精度の改良となる.
図 3 は横軸に σ をとり縦軸に平方平均誤差をプロッ
いろいろな例を観察すると,図 5 のように補正に
トしたものである.ただし,平方平均誤差は各 σ に
よって精度が低下するのはほとんどの場合,最尤推
112
D = 2σ
定で当てはめた楕円が真の楕円の内側に来る場合で
あった.しかし,最尤推定では真の楕円の外側に当
[11]
てはまる場合のほうが圧倒的に多いので,全体とし
ては精度の向上となっている.
[12]
真の楕円の外側に当てはまる可能性が多いのは,楕
円の長軸または短軸方向の半径を a とすると,x2 ま
[13]
たは y 2 の係数が 1/a2 に比例し,これを未知数とし
て推定するためと思われる.実際,1/a2 が真値 1/¯
a2
より大きくなる確率と小さくなる確率が等しいとし
ても,関数 y = 1/x2 のグラフの形からわかるよう
[14]
[15]
に,a ≥ a
¯ となる確率が圧倒的に大きくなる.
[16]
9. まとめ
[17]
平面上の点列に楕円を当てはめる問題は,従来か
ら最尤推定が最も高精度であるとみなされ,くりこ
[18]
み法,HEIV 法,FNS 法などの計算法が提案されて
いる.本論文ではこれを上回る「超精度」の当ては
[19]
め法が存在することを理論的および実験的に検証し
[20]
た.これは最尤推定の摂動解析によって高次の偏差
[21]
項を評価し,これを最尤推定解から差し引くもので
ある.ただし,最尤推定解は KCR 下界をほとんど
[22]
達成しているので,超精度補正による効果はわずか
[23]
であり,実際の応用で目に見えるほどの差を生じる
ことはない.しかし,当てはめの精度の理論限界と
最尤推定の精度に関して,理論的な観点から意味の
[24]
ある結果であるといえる
[25]
謝辞: 本研究の一部は文部科学省科学研究費基盤研究 C
(No. 17500112) の助成によった.
[26]
参考文献
[27]
[1] D. A. Anderson, The circular structural model, J. R.
Statist. Soc. B-43-2 (1981), 131–141.
[2] 浅山泰祐, 塩野 充, 円弧近似を用いた任意傾斜楕円のハフ
変換による検出実験, 電子情報通信学会技術報告, PRMU98123 (1998-11).
[3] M. Berman and D. Culpin, The statistical behaviour of
some least squares estimators of the centre and radiusl
of a circle, J. R. Statist. Soc. B-48-2 (1986), 183–196.
[4] F. J. Bookstein, Fitting conic sections to scattered data,
Comput. Graphics Image Process., 9 (1979) 56–71.
[5] J. Cabrera and P. Meer, Unbiased estimation of ellipses
by bootstrapping, IEEE Trans. Patt. Mach. Intelli.,
18-7 (1996-7), 752–756.
[6] N. N. Chan, On circular functional relationships, J. R.
Statist. Soc. B-27 (1965), 45–56.
[7] Y. C. Cheng and S. C. Lee, A new method for quadratic
curve detection using K-RANSAC with acceleration
techniques, Patt. Recog, 28-5 (1995-5), 663–682.
[8] N. Chernov and C. Lesort, Statistical efficiency of curve
fitting algorithms, Comput. Stat. Data Anal., 47-4
(2004-11), 713–728.
[9] W. Chojnacki, M. J. Brooks, A. van den Hengel and D.
Gawley, On the fitting of surfaces to data with covariances, IEEE Trans. Patt. Anal. Mach. Intell., 22-11
(2000), 1294–1303.
[10] D. B. Cooper and N. Yalabik, On the computational
cost of approximating and recognizing nose-perturbed
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
113
straight lines and quadratic arcs in the plane, IEEE
Trans. Comput., 25-10 (1976-10), 1020–1032.
A. Fitzgibbon, M. Pilu and R. B. Fisher, Direct least
square fitting of ellipses, IEEE Trans. Patt. Anal.
Mach. Intell., 21-5 (1999-5), 476–480.
藤本公三, 岩田剛治, 仲田周次, θ-ρ ハフ変換平面からの2
次曲線のパラメータ抽出, 電子情報通信学会論文誌 D-II,
J74-D-II (1991-9), 1184–1191.
N. Guil and E. L. Zapata, Lower order circle and ellipse
Hough transform, Patt. Recog., 30-10 (1997-10), 1729–
1744.
C.-T. Ho and L.-H. Chen, A fast ellipse/circle detector
using geometric symmetry, Patt. Recog., 28-1 (1995),
117–124.
C. L. Huang, Elliptical feature extraction via an improved Hough transform, Patt. Recog. Lett., 10-2 (19892), 93–100.
J. Illingworth and J. Kittler, IEEE Trans. Patt. Anal.
Mach. Intell., 9-5 (1999-9), 690–698.
D. Ioannou, W. Huda, and A. F. Laine, Circle recognition through a 2D Hough Trnasform and radius histogramming, Image Vision Comput. 17-1 (1999-1), 15–
26.
S. H. Joseph, Unbiased least squares fitting of curcular
arcs, CVGIP: Models Image Process., 56-5 (1994-9),
424–432.
K. Kanatani, Geometric Computation for Machine Vision, Oxford University Press, Oxford, U.K., 1993.
金谷健一, コンピュータビジョンのためのくりこみ法, 情報
処理学会論文誌, 35-2 (1994-2), 201–209.
K. Kanatani, Statistical bias of conic fitting and renormalization, IEEE Trans. Patt. Anal. Mach. Intell., 163 (1994-3), 320–326.
金谷健一, 当てはめ問題の最適推定と精度の理論限界, 情報
処理学会論文誌, 36-8 (1995-8), 1865–1873.
K. Kanatani, Statistical Optimization for Geometric Computation: Theory and Practice, Elsevier Science, Amsterdam, The Netherlands, 1996, reprinted by
Dover, New York, 2005.
金谷 健一, 画像からの幾何学的推論はどういう統計的モデル
に基づくのか,電子情報通信学会論文誌 D-II, J86-D-II-7
(2003-7), 966-973.
金谷 健一,くりこみ法その後: 波紋と発展, 情報処理学会研
究報告, 2003-CVIM-139-5 (2003-7), 33–40.
金谷 健一, 最尤推定の最適性と KCR 下界,情報処理学会研
究報告, 2005-CVIM-147-8 (2005-1), 59-64.
金谷 健一, くりこみ法の謎を解く, 情報処理学会研究報告,
2005-CVIM-149-3 (2005-5), 15-22.
Y. Kanazawa and K. Kanatani, Optimal conic fitting
and reliability evaluation, IEICE Trans. Inf. & Sys.,
E79-D-9 (1996-9), 1323–1328.
D. Keren, D. Cooper and J. Subrahmonia, Describing complicated objects by implicit polynomials, IEEE
Trans. Patt. Anal. Mach. Intell., 16-1 (1994-1), 38–53.
A. Kukush, I. Markovsky and S. Van Huffel, Consistent
estimation in an implicit quadratic measurement error
model, Comput. Stat. Data Anal., 47-1 (2004-8), 123–
147.
A. Kukush and E. O. Maschlce, The efficiency of adjusted least squares in the linear functional relationship,
J. Multivariate Anal., 87-2 (2003-10), 261–274.
Y. Leedan and P. Meer, Heteroscedastic regression in
computer vision: Problems with bilinear constraint, Int.
J. Comput. Vision., 37-2 (2000), 127–150.
V. J. Milenkovic, Multiple resolution search techniques
for the Hough transform in high dimensinonal parameter spaces, in A. Rosenfeld (Ed.), Techniques for 3-D
Machine Perception, Elsevier Science, Amsterdam, The
Netherlands, 1986, pp. 231–254.
三島等, 太田直哉, 金谷健一, 信頼性評価を備えた最適なコニッ
ク当てはめプログラム, 情報処理学会研究報告, 98-CVIM111-4 (1998-5), 25–32.
森 克己, 渡邊栄治, 片桐重和, 弦とその共役直径に基づく
楕円弧の識別, 電子情報通信学会論文誌 D-II J84-D-II-2
(2001-2), 287–294.
[36] 永田亮一, 川口 剛, 遺伝アルゴリズムとエッジ点の線分領
域へのグルーピングを利用する楕円検出, 電子情報通信学会
論文誌 D-II, J81-D-II-9 (1998-9) 2074–2085.
[37] 岡部光生, 金谷健一, 太田直哉, 楕円成長法による円形物体
の自動検出, 電子情報通信学会論文誌 D-II, J85-D-II-12
(2002-12), 1823–1831.
[38] 恩田寿和, 藤原伸行, 阿部清秀, 森宣仁, 三次元円検出による
部品位置決めと事前のハンド干渉チェックにより実現した視
覚ベースビンピッキングシステム, 日本ロボット学会誌, 18-7
(2000-10), 995–1002.
[39] 恩田邦夫, 渡並 智, 青木由直, Hough 変換平面からの楕円
パラメータ決定法, 電子情報通信学会論文誌 D-II, J72-DII-10 (1989-10), 1760–1764.
[40] D. Pao, H. F. Li, and R. Jayakumar, A decomposable parameter space for the detection of ellipses, Patt.
Recog. 14-12 (1993), 951–958.
[41] J. Porrill, Fitting ellipses and predicting confidence envelopes using a bias corrected Kalman filter, Image Vision Comput., 8-1 (1990-2), 37–41.
[42] V. Pratt, Direct least-squares fitting of algebraic surfaces, Comput. Graphics, 21-4 (1987-7), 145–152.
[43] G. Roth and M. D. Levine, Extracting geometric primitives, CVGIP: Image Understand., 58-1(1993-7), 1–22.
[44] G. Roth and M. D. Levine, Geometric primitive extraction using a genetic alorithm IEEE Trans. Patt. Anal.
Mach. Intell., 16-9 (1994-9), 901–905.
[45] P. L. Rosin, Ellipse fitting by accumulating five-points
fits, Patt. Recog. Lett., 14-6 (1993-8), 661–669.
[46] P. L. Rosin, A note on the least squares fitting of ellipses, Patt. Recog. Lett., 14 (1993), 799–808.
[47] P. L. Rosin and G. A. W. West, Nonparametric segmentation of curves into various representations, IEEE
Trans. Patt. Anal. Mach. Intell., 17-12 (1995-12),
1140–1153.
[48] R. Safee-Rad, I. Tchoukanov, B. Benhabib and K.
C. Smith, Accurate parameter estimation of quadratic
curves from gray-level images, CVGIP: Image Understand., 54-2 (1991-9), 259–274.
[49] P. D. Sampson, Fitting conic sections to “very scattered” data: An iterative refinement of the Bookstein algorithm, Comput. Graphics Image Process., 18 (1982),
97–108.
[50] 塩野 充, 黒点ランダム抽出と Hough 曲面の交点計算によ
る円図形検出の一手法, 情報処理学会論文誌 32-2 (1991),
179–187.
[51] 塩野 充, 黒点ランダム抽出と重心を用いたハフ変換による円
弧の検出実験, 電子情報通信学会論文誌, J75-D-II (1992-7),
1195–1201.
[52] H. Sp¨
ath, Least-squares fitting of ellipses and hypoerbolas, Comput. Statistics, 12-3 (1997-9) 329–341.
[53] H. Sp¨
ath, Orthogonal distance fitting by circles and ellipses with given area, Comput. Statistics, 12-3 (1997-9)
343–354.
[54] G. Taubin, Estimation of planar curves, surfaces, and
non-planar space curves defined by implicit equations
with applications to edge and rage image segmentation,
IEEE Trans. Patt. Anal. Mach. Intell., 13-11 (199111), 1115–1138.
[55] G. Taubin, F. Cukierman, S. Sullivan, J. Ponce and
D. J. Kriegman, Parameterized families of polynomials
for bounded algebraic curve and surface fitting, IEEE
Trans. Patt. Anal. Mach. Intell., 16-3 (1994-3), 287–
303.
[56] S. Tsuji and F. Matsumoto, Detection of ellipses by a
modified Hough transform, IEEE Trans. Comput., 27-8
(1978-8), 777–781.
[57] 渡辺孝志, 畠山雅充, 木村彰男, ハフ変換を用いた接線情報
の抽出と欠損楕円の検出, 電子情報通信学会論文誌 D-II,
J82-D-II-12 (1999-12), 2221–2229.
[58] 渡辺孝志, 木村彰男, 丹波澄雄, 横山隆三, Li-Lavin-Le Master 型高速ハフ変換による欠損楕円の検出, 電子情報通信学
会論文誌 D-II, J76-D-II-12 (1993-12), 2504–2512.
[59] 渡辺孝志, 柴田俊浩, Hough 変換と階層化画像を用いた欠
損楕円の検出, 電子情報通信学会論文誌 D-II, J73-D-II-2
(1990-2), 159–166.
[60] M. Werman and D. Keren, A Bayesian method for filtering parametric and nonparametric models to noisy data,
IEEE Trans. Patt. Anal. Mach. Intell., 23-5 (2001-5),
528–534.
[61] W.-Y. Wu and M.-J. J. Wang, Elliptical object detection by using its geometric properties, Patt. Recog., 2610 (1993-10), 1449–1500.
[62] 大和淳二, 入澤和義, 石井郁夫, 牧野秀夫, 重み付け中点図
面を用いた楕円抽出アルゴリズム, 電子情報通信学会論文誌
D-II, J72-D-II-7 (1989), 1009–1016.
[63] H. K. Yuen, J. Illingworth, and J. Kittler, Detecting
partially occluded ellpises using the Hough transform,
Image Vision Comput., 7-1 (1989-2), 31–37.
[64] J. H. Yoo and I. K. Seth, An ellipse detection method
from the polar and pole definition of conics, Pattern
Recog., 26-2 (1993-2), 307–315.
付録. 最尤推定量の摂動解析
1. u は単位ベクトルであるという条件のもとに摂動する
ので,次式が成立しなければならない.
u + ∆1 u + ∆2 u + · · ·
2
= (u + ∆1 u + ∆2 u + · · · , u + ∆1 u + ∆2 u + · · ·)
= 1
(44)
O(1), O(σ), O(σ 2 ) の項をそれぞれ取り出すと次のように
なる.
(u, u) = u 2 = 1,
(u, ∆1 u) = 0
(45)
(u, ∆2 u) = −(∆1 u, ∆1 u) = − ∆1 u
2
(46)
¯
2. 式 (17) の両辺から O(1) の項を取り出すと M u = 0
となる.
3. 式 (17) の両辺から O(σ) の項を取り出すと次のように
なる.
¯ ∆1 u + ∆1 M u + ∆∗1 M u = 0
M
(47)
式 (18) より次のようになる.
N
∆∗1 M u = −2
α=1
((∆1 u,V0 [ξ α ]u)+O(σ 2 ))ξ α (ξ α ,u)
(u, V0 [ξ α ]u)2
= 0
(48)
¯
ゆえに式 (47) に左から M
−
を掛けると次式が得られる.
−
¯ ∆1 M u
P u ∆1 u = −M
(49)
¯ −M
¯
ただし,P u は式 (23) で定義した射影行列である(M
= P u であることに注意 [23]).式 (45) の第2式より ∆1 u
は u に直交するから,P u ∆1 u = ∆1 u である.したがっ
て式 (20) が得られる.
4. 式 (17) の両辺から O(σ 2 ) の項を取り出すと次のよう
になる(式 (48) より ∆∗1 M u には O(σ 2 ) の項は含まれて
いない).
¯ 2 u+∆1 M ∆1 u+∆∗1 M ∆1 u+∆2 M u +∆∗2 M u
M∆
= ∆2 Lu
(50)
¯ − を掛け,移項すると次式を得る.
両辺に左から M
¯ − ∆2 M u + M
¯ − ∆1 M M
¯ − ∆1 M u
P u ∆2 u = −M
−
−
−
∗
∗
¯ ∆1 M ∆1 u − M
¯ ∆2 M u + M
¯ ∆2 Lu (51)
−M
これは ∆2 u の u に直交する成分である.u 方向の成分は
式 (46) より − ∆1 u 2 u であり,式 (20) より ∆1 u 2 =
¯ − ∆1 M u 2 u である.以上より式 (21) を得る.
− M
114