楕円フーリエ解析
06 Mar 2023
輪郭形状の解析のための楕円フーリエ解析について、理論から実装までのメモです。
楕円フーリエ解析とは
輪郭 contour の解析のために用いられる解析方法です。 形態計測学 Morphometricsで用いられます。
- 輪郭がそれぞれのサンプルでどのくらい似ているのか
- サンプル同士の輪郭形状の違いがどのようなものか
といったことを定量化および可視化することができます。 研究例としては、以下のような研究が挙げられます。 本ページでは、理論からRによる実装、解析を含めた流れをできる限り簡単に説明することを試みます。
チェインコード
文献によって方向や開始点などが若干異なりますが、方向に数字を結びつけたものをチェインと呼び、その数字を並べたベクトルをチェインコードと呼びます。
\[V = a_1 a_2 a_3 \cdots a_K\]ここで、\(V\)がチェインコードで\(a_1 a_2 a_3 \cdots a_K\)がそれぞれのチェインになります。 おそらく最初に提唱されたのはFreeman(1974)だそうです。 0度から360度まで、45度ごと反時計回りに矢印を伸ばします。 矢印は8本になるため、それぞれの矢印に0から7の番号を振ります。 図示すると、以下のようになります。
簡単な例を挙げて輪郭形状をチェインコードで表してみます。
ここでは、すべての方向を含むような単純な八角形を考えます。
チェインコードは以下のように表されます。
\[V = 07654321\]図示すると以下のようになります。 開始点は基本的には左上になります。
この八角形は正八角形ではないことに注意が必要です。
すなわち、すべての辺の長さが同じではないということです。
縦横に移動する場合は長さが\(1\)になりますが、斜めに移動する場合は長さは\(\sqrt{2}\)になります。
それぞれのチェインの\(x\)座標と\(y\)座標の移動とその長さの対応は以下の表のようになります。
| チェイン | \(x\)方向変位 | \(y\)方向変位 | 長さ |
|---|---|---|---|
| 0 | 1 | 0 | 1 |
| 1 | 1 | 1 | \(\sqrt{2}\) |
| 2 | 0 | 1 | 1 |
| 3 | -1 | 1 | \(\sqrt{2}\) |
| 4 | 0 | -1 | 1 |
| 5 | -1 | -1 | \(\sqrt{2}\) |
| 6 | 0 | -1 | 1 |
| 7 | 1 | -1 | \(\sqrt{2}\) |
一本の矢印の長さ(大きさ)が輪郭をなぞる時間と対応すると考えます。 例えば、0, 2, 4, 6の矢印は1、1, 3, 5, 7の矢印は\(\sqrt{2}\)の時間がかかります。 つまり、\(p\)番目までのチェインまでをなぞる時間\(t_p\)は以下の式で表されます。
\[t_p = \sum^p_{i=1} \Delta t_i\]ここで、\(\Delta t_i\)は一つの矢印をなぞる時間です。 また、\(p\)番目までのチェインまでの\(x\)座標と\(y\)座標の変位はそれぞれ以下の式で表されます。
\[x_p = \sum^p_{i=1} \Delta x_i\] \[y_p = \sum^p_{i=1} \Delta y_i\]輪郭におけるフーリエ級数展開
ここから、フーリエ級数展開などの話になります。 フーリエ級数展開、フーリエ変換などについての基本的な理論についてはここでは省略します。 本ページの最終目標は輪郭の解析をすることなので、難しい場合は読み飛ばしていただけたらと思います。
\(x\)座標におけるフーリエ級数展開 Fourier series expansionは以下の式で表されます。
\[x(t) = A_0 + \sum^\infty_{n=1} a_n \cos \frac{2n \pi t}{T} + b_n \sin \frac{2n \pi t}{T}\]ここで、\(T\)は周期で、輪郭形状の場合はなぞる時間と言えます。 \(A_0, a_n, b_n\)はそれぞれ以下の式で表されます。
\[A_0 =\]